Computational Method for Dynamic Analysis of Multibody Systems with Deformable Elements

1. Introduction

The elastodynamic analysis of mobile mechanical systems is a highly debated research topic, especially in the context of the dynamics of mobile mechanical systems with elastic elements. Most of the research conducted among time, has aimed at accurately capturing the vibrations and flexibility of the components within the structure of the studied mechanisms, as well as experimentally validating these findings, as seen in [1,2].

The studies carried out on mobile mechanical systems have enabled the analysis of elastodynamic phenomena and vibrations under high-speed operating conditions. The use of the finite element method in the context of these analyses provides a detailed representation of deformations, fully integrated with the kinematic and dynamic models of the studied mechanisms. Furthermore, by complementing theoretical research with experimental validations for direct comparison with simulations and numerical processing, the obtained results have only increased the credibility and relevance of the research, as observed in [3,4].

Addressing to the elastodynamic analysis, and targeting mobile mechanical systems through an integrated concept of theoretical (mathematical) modeling – numerical processing – virtual prototyping – experimental validation, as exemplified in [5], these provides a classical methodological extensibility. This is compatible with formalisms such as Lagrange, Gibbs-Appell, and Hamilton, which, when combined with the finite element method or multibody dynamics, it generates modern research trends aimed at ensuring continuity in the field.

On the other hand, the studies presented in [6,7,8], which follow the aforementioned approach, are highly complex, as the use of the finite element method and the comparison of modal-dynamic analysis results, and these involves a large volume of numerical processing with high computation times, especially for extended models. Another particularly important aspect is the validation of the methods and tools, which are specific to the elastodynamic analysis. These were performed on only one mechanism, implying a limited validation. In this context, validating other mechanism configurations requires adapting the methods and newly validating the undertaken research.

If we consider the research in [9,10], which highlights measurement errors and inconsistencies in the obtained results, it becomes evident that experimental determinations are influenced by vibrations and the conditions under which experimental testing was conducted. Consequently, the results reveal incomplete alignments or discrepancies compared to the numerical processing of the mathematical models and, implicitly, with those obtained through virtual prototyping.

Moreover, the research in [10] also addresses to another highly important aspect in the dynamic behaviour of mobile mechanical systems, namely when the kinematic joints of a spatial mechanism exhibit a mechanical clearance. Thus, the analysis targets an advanced dynamic modelling of a spatial mechanism that includes impact forces and friction inside the joints, specifically adapted for microgravity conditions. As for the experimental analysis, it places greater emphasis on the joint clearance, which in microgravity can significantly degrade the mechanism motion accuracy. It is also noted that the methods used in the mechanism dynamic analysis have a limited applicability, and for generalizing them to be applied to more complex mechanisms, it requires substantial adaptations.

Similar to the studies presented in [10], the research in [11] also highlights several study methods applied to planar mechanisms, but with ambiguous experimental validation, which can only be certified through virtual prototyping.

It should be noted that a crucial influence on an elastodynamic analysis involving a mobile mechanical system is the dependence on the mechanical properties of the materials from which the system components are made. For this purpose, studies such as [12,13] have considered components made from materials with linear elasticity. However, the continuation of research in this direction also targets elastodynamic analyses of mobile mechanical systems with large component deformations or the use of materials with nonlinear mechanical properties, such as those described in [14,15,16].

In the research from [17], the authors develop a hybrid elastodynamic modelling, in which the mobile mechanical system is decomposed into substructures, some of them simulated numerically, others tested experimentally. The overall model of the mechanical system is then assembled using frequency-based substructuring methods. The aim of this research is to model the mechanical system in order to produce more realistic predictions regarding vibrations, system stability, and dynamic loads. Thus, the scientific approach involves a high methodological complexity through numerical modelling and experimental testing, which may lead to potential errors in the results and demands rigorous experimental procedures.

Therefore, the innovative aspect that ensures the continuity of research involving mobile mechanical systems elastodynamic analyses, lies in the development of mathematical models that also take into account the impact forces within joints where are present the potential deviations. These deviations significantly influence the quality of the results through which the system is rigorously validated, especially from an experimental viewpoint , as it can be seen in [18,19].

Considerring the research highlighted above and the weaknesses that an elastodynamic analysis of a mobile mechanical system may have, this paper aims to achieve the following objectives: Coupling the motion of a rigid body with one of a deformable body and identifying the static and dynamic components of the matrices in the equations of motion – interfacing the Kane method with the finite element method; Modal dynamic analysis through the modal superposition method with the nodal displacements identification; Elastodynamic analysis by determining displacements, velocities, and accelerations in deformable body motion; Developing mathematical models on a crank-connecting rod mechanism used to drive a three-cylinder inline internal combustion engine; Elastodynamic analysis of the mechanism using MSC Adams software and comparing the results with those obtained from mathematical model processing; Experimental analysis of the mechanism by monitoring the variation of the kinematic and dynamic parameters in both time and frequency domains.

These objectives are achieved according to the research workflow presented in Figure 1.

2. Mathematical Models for Dynamic Analysis of Mechanisms with Deformable Elements

There will be considered a kinematic linkage made of ,,n–1” solid–rigid, connected together through kinematic joints as it can be remarked in Figure 2. Thus, there will be accepted the solid “n” is a deformable solid.

The following notations will be introduced:

$R_i\left({\overrightarrow{x}}_i,{\overrightarrow{y}}_i,\overrightarrow{z_i}\right)$ – is the reference coordinate system attached to the link “i” with versor base of $W_i\left({\overrightarrow{i}}_i,{\overrightarrow{j}}_i,\overrightarrow{k_i}\right)$;

$R_0\left({\overrightarrow{x}}_0,{\overrightarrow{y}}_0,\overrightarrow{z_0}\right)$ – is the global reference coordinate system with a versor base of$W_0\left({\overrightarrow{i}}_0,{\overrightarrow{j}}_0,\overrightarrow{k_0}\right)$;

${\overrightarrow{\delta }}_i$– is the relative translational vector between ,,i” and ,,i–1”, depending on ,,Ri-1” tried, which exists, if between ,,i” and ,,i-1” links, there is a joint;

${\stackrel{-}{r}}_i$ – position vector of the ,,O’i” point, expressed in the ,,Ri–1” reference coordinate system, when the relative translation begins.

Thus, the following equations can be written:

$${\overrightarrow{r}}_i={\left\{r_i^x,r_i^y,r_i^z\right\}}_{i-1}={\left\{r_i\right\}}^T\left\{{\overrightarrow{W}}_{i-1}\right\}$$

(1)

$${\overrightarrow{\delta }}_i={\left\{{\delta }_i^x,{\delta }_i^y,{\delta }_i^z\right\}}_{i-1}={\left\{{\delta }_i\right\}}^T\left\{{\overrightarrow{W}}_{i-1}\right\}$$

(2)

There will be introduced the transformation matrices for crossing from a reference coordinate system to another, by having in sight Equation (3).

$$\left\{{\overrightarrow{W}}_{i-1}\right\}=\left[A_{oi-1}\right]\left\{{\overrightarrow{W}}_o\right\}$$

(3)

Thus, the Equations (1) and (2) becomes:

$$\overline{r_i}={\left\{r_i\right\}}^T\left\{{\overrightarrow{W}}_{i-1}\right\}={\left\{r_i\right\}}^T\left[A_{oi-1}\right]\left\{{\overrightarrow{W}}_o\right\}$$

(4)

$$\overline{{\delta }_i}={\left\{{\delta }_i\right\}}^T\left\{{\overrightarrow{W}}_{i-1}\right\}={\left\{{\delta }_i\right\}}^T\left[A_{oi-1}\right]\left\{{\overrightarrow{W}}_o\right\}$$

(5)

By having in sight the Figure 3, onto E_n link, there will be considered the M_n point and u(M,t) displacement. The position vector of the P point is:

$$\overrightarrow{r}=\overrightarrow{r_{n-1}}+\overrightarrow{r_u}$$

(6)

$$\overrightarrow{r_u}=\overrightarrow{r_n}+\overrightarrow{u}$$

(7)

$$\overrightarrow{r_{n-1}}={\displaystyle \sum _{i=1}^n\left({\left\{r_i\right\}}^T+{\left\{{\delta }_i\right\}}^T\right)\left[A_{0\hspace{0.17em}i-1}\right]}\left\{\overrightarrow{W_0}\right\},$$

(8)

where: $\left[A_{0i-1}\right]$ represents the coordinates transformation matrix.

The elastic displacement vector is:

$$\overrightarrow{u}={\left\{u\right\}}^T\left\{\overrightarrow{w_n}\right\}$$

(9)

$$\left\{u\right\}=\left[\text{}N\right]\left\{q\right\},$$

(10)

where: $\left[N\right]$– represents the form functions matrix or interpolating polynomial matrix; $q$– is the nodal displacements vector.

$$\left\{\overrightarrow{w_n}\right\}=\left[A_{0,n}\right]\left\{\overrightarrow{w_0}\right\}$$

(11)

Based on previous equations, there will be obtained:

$$\overrightarrow{u}={\left\{q\right\}}^T{\left[N\right]}^T\left[A_{0,n}\right]\left\{\overrightarrow{w_0}\right\}$$

(12)

$$\overrightarrow{r_n}={\left\{r_n\right\}}^T\left\{\overrightarrow{w_n}\right\}={\left\{r_n\right\}}^T\left[A_{0,n}\right]\left\{\overrightarrow{w_0}\right\}$$

(13)

$$\begin{array}{l}\overrightarrow{r_u}={\left\{r_n\right\}}^T\left[A_{0,n}\right]\left\{\overrightarrow{w_0}\right\}+{\left\{q\right\}}^T{[N]}^T\left[A_{0,n}\right]\left\{\overrightarrow{w_0}\right\}=\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\left({\left\{r_n\right\}}^T+{\left\{q\right\}}^T{[N]}^T\right)\left[A_{0,n}\right]\left\{\overrightarrow{w_0}\right\}\end{array}$$

(14)

There will be introduced the Equations (8) and (14) in Equation (6):

$$\overrightarrow{r}=\left({\displaystyle \sum _{i=1}^n\left({\left\{r_i\right\}}^T+{\left\{{\delta }_i\right\}}^T\right)}\left[A_{0\hspace{0.17em}i-1}\right]+\left({\left\{r_n\right\}}^T+{\left\{q\right\}}^T{[N]}^T\right)\left[A_{0,n}\right]\right)\left\{\overrightarrow{w_0}\right\}$$

(15)

Thus, there will be considered the following generalized coordinates:

$${\left\{p_i\right\}}^T={\left\{\left\{r_i\right\}\hspace{0.17em}\hspace{0.17em}\left\{{\delta }_i\right\}\right\}}^T$$

(16)

It will be obtained:

$$\overrightarrow{r}=\left({\displaystyle \sum _{i=1}^n{\left\{p_i\right\}}^T}\left[A_{0\hspace{0.17em}i-1}\right]+\left({\left\{r_n\right\}}^T+{\left\{q\right\}}^T{[N]}^T\right)\left[A_{0,n}\right]\right)\left\{\overrightarrow{w_0}\right\}$$

(17)

For the speeds calculus, there will be differentiated the (44) in raport with time :

$$\stackrel{•}{\overrightarrow{r}}=\stackrel{•}{\overrightarrow{r_{n-1}}}+\stackrel{•}{\overrightarrow{r_u}}$$

(18)

Also, the Equation (14) it will be differentiated in raport with the time, and there will be obtained :

$$\stackrel{•}{\overrightarrow{r_u}}=\stackrel{•}{\overrightarrow{r_n}}+\stackrel{•}{\overrightarrow{u}}=\stackrel{•}{\overrightarrow{r_u}}+{\displaystyle \frac{\partial u}{\partial t}}+\overrightarrow{\omega }\times \overrightarrow{u}$$

(19)

By reversing at the workflow with the use of matrices, we will differentiate the Equation (15), depending on time and we will obtain :

$$\begin{array}{l}\stackrel{•}{\overrightarrow{r}}=\hspace{0.17em}\left({\displaystyle \sum _{i=1}^n{\left\{\stackrel{•}{{\delta }_i}\right\}}^T\left[A_{0,i-1}\right]}+{\displaystyle \sum _{i=1}^n\left({\left\{r_i\right\}}^T+{\left\{{\delta }_i\right\}}^T\right)\left[{\tilde{\omega }}_{0,i-1}\right]\left[A_{0,i-1}\right]}+{\displaystyle \sum _{i=1}^n{\left\{\stackrel{•}{q}\right\}}^T{\left[N\right]}^T\left[A_{0,n}\right]}+\right.\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left.+{\displaystyle \sum _{i=1}^n\left({\left\{r_n\right\}}^T+{\left\{q\right\}}^T{[N]}^T\right)\hspace{0.17em}\left[{\tilde{\omega }}_{0n}\right]\left[A_{0n}\right]}\right)\left\{\overrightarrow{w_0}\right\}\end{array}$$

(20)

By taking into account the Equation (16), the Equation (20) becomes :

$$\begin{array}{l}\overrightarrow{v}={\displaystyle \frac{d\overrightarrow{r}}{dt}}=\left({\displaystyle \sum _{i=1}^n{\left\{\stackrel{•}{p_i}\right\}}^T\left[A_{0,i-1}\right]}+{\left\{p_i\right\}}^T\left[{\tilde{\omega }}_{0,i-1}\right]\left[A_{0,i-1}\right]+\right.{\left\{\stackrel{•}{q}\right\}}^T{\left[N\right]}^T\left[A_{0,n}\right]+\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left.\left({\left\{r_n\right\}}^T+{\left\{q\right\}}^T{[N]}^T\right)\hspace{0.17em}\left[{\tilde{\omega }}_{0n}\right]\left[A_{0n}\right]\right)\left\{\overrightarrow{w_0}\right\}\end{array}$$

(21)

Also we will obtain the accelerations mathematical expression:

$$\begin{array}{l}\overrightarrow{a}={\displaystyle \frac{d^2\overrightarrow{r}}{d{t}^2}}=\left({\displaystyle \sum _{i=1}^n{\left\{\stackrel{••}{{\delta }_i}\right\}}^T\left[A_{0,i-1}\right]}+2{\displaystyle \sum _{i=1}^n{\left\{{\stackrel{•}{\delta }}_i\right\}}^T\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}}+\right.\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{i=1}^n\left({\left\{r_i\right\}}^T+{\left\{{\delta }_i\right\}}^T\right)\hspace{0.17em}\left[{\tilde{\epsilon }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}}+\\ +{\displaystyle \sum _{i=1}^n\left({\left\{r_i\right\}}^T+{\left\{{\delta }_i\right\}}^T\right)\hspace{0.17em}\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}+}\\ +{\displaystyle \sum _{i=1}^n{\left\{\stackrel{••}{d}\right\}}^T{\left[N\right]}^T\left[A_{o,n}\right]}+2{\displaystyle \sum _{i=1}^n{\left\{\stackrel{•}{d}\right\}}^T{\left[N\right]}^T\left[{\tilde{\omega }}_{o,n}\right]\left[A_{o,n}\right]}+\\ +{\displaystyle \sum _{i=1}^n\left({\left\{r_n\right\}}^T+{\left\{d\right\}}^T{\left[N\right]}^T\right)}\hspace{0.17em}\hspace{0.17em}\left[{\tilde{\epsilon }}_{0,n}\right]\hspace{0.17em}\hspace{0.17em}\left[A_{0,n}\right]+\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left.{\displaystyle \sum _{i=1}^n\left({\left\{r_n\right\}}^T+{\left\{d\right\}}^T{\left[N\right]}^T\right)}\hspace{0.17em}\left[{\tilde{\omega }}_{0,n}\right]\hspace{0.17em}\left[{\tilde{\omega }}_{0,n}\right]\hspace{0.17em}\hspace{0.17em}\left[A_{0,n}\right]\right)\hspace{0.17em}\left\{\overrightarrow{w_0}\right\}\end{array}$$

(22)

Based on the Equation (16), the Equation (22) becomes :

$$\begin{array}{l}\overrightarrow{a}={\displaystyle \frac{d\overrightarrow{v}}{dt}}=\left({\displaystyle \sum _{i=1}^n{\left\{\stackrel{••}{p_i}\right\}}^T\left[A_{0,i-1}\right]}+2{\displaystyle \sum _{i=1}^n{\left\{{\stackrel{•}{p}}_i\right\}}^T\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}}+\right.\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{i=1}^n{\left\{p_i\right\}}^T\hspace{0.17em}\left[{\tilde{\epsilon }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}}+{\displaystyle \sum _{i=1}^n{\left\{p_i\right\}}^T\hspace{0.17em}\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}+}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{i=1}^n{\left\{\stackrel{••}{q}\right\}}^T{\left[N\right]}^T\left[A_{o,n}\right]}+2{\displaystyle \sum _{i=1}^n{\left\{\stackrel{•}{q}\right\}}^T{\left[N\right]}^T\left[{\tilde{\omega }}_{o,n}\right]\left[A_{o,n}\right]}+\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{i=1}^n\left({\left\{r_n\right\}}^T+{\left\{q\right\}}^T{\left[N\right]}^T\right)}\hspace{0.17em}\hspace{0.17em}\left[{\tilde{\epsilon }}_{0,n}\right]\hspace{0.17em}\hspace{0.17em}\left[A_{0,n}\right]+\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left.{\displaystyle \sum _{i=1}^n\left({\left\{r_n\right\}}^T+{\left\{q\right\}}^T{\left[N\right]}^T\right)}\hspace{0.17em}\left[{\tilde{\omega }}_{0,n}\right]\hspace{0.17em}\left[{\tilde{\omega }}_{0,n}\right]\hspace{0.17em}\hspace{0.17em}\left[A_{0,n}\right]\right)\hspace{0.17em}\left\{\overrightarrow{w_0}\right\}\end{array}$$

(23)

Thus, there will be elaborated the proper mathematical models for the motion equations in Kane formalism, for the holonomic case are [22]:

$$Q+Q^{*}=0,$$

(24)

where: Q represents the generalized forces associated to the exterior loadings and connection forces; Q^* represents the generalized forces associated to inertial torsors.

$$Q={\displaystyle \sum _{k=1}^p\overrightarrow{F_k}}\cdot {\overrightarrow{v}}_{rk}+{\displaystyle \sum _{i=1}^n{\overrightarrow{M}}_i}\cdot {\overrightarrow{\omega }}_{ri}$$

(25)

where: ${\overrightarrow{F}}_k$ is the resultant force associated to M_k point; ${\overrightarrow{v}}_{rk}$ is the partial speed associated to $\overrightarrow{v_k}$ speed of Mk point; $\overrightarrow{M_i}$is the resulted torque associated to the “i” kinematic element. $\overrightarrow{{\omega }_{ri}}$ is the partial angular speed, associated to the “i” element and to v_r generalized speed.

$$Q^{*}={\displaystyle \sum _{k=1}^p\overrightarrow{F_{i{n}_k}}\overrightarrow{v_k}}+{\displaystyle \sum _{i=1}^s\overrightarrow{M_{in}}\overrightarrow{{\omega }_{ri}}}+{\displaystyle \sum _{i=1}^s\overrightarrow{F_{i{n}_i}}}\overrightarrow{v_{r0}}$$

(26)

For the “i” element, by taking into account the “q” nodal displacements, the generalized force associated to the inertial torsors, has the following form:

$$Q_i^{*}=\overrightarrow{F_{i{n}_i}}{\displaystyle \frac{\partial \overrightarrow{v_i}}{\partial \stackrel{•}{q}}}+\overrightarrow{M_{i{n}_i}}{\displaystyle \frac{\partial \overrightarrow{{\omega }_i}}{\partial \stackrel{•}{q}}}$$

(27)

For the “i” element, it can be written:

$$\overrightarrow{F_{i{n}_i}}=-m_i\overrightarrow{a_i}$$

(28)

$$\overrightarrow{M_{i{n}_i}}{\displaystyle \frac{\partial \overrightarrow{{\omega }_i}}{\partial \stackrel{•}{q}}}\approx 0$$

(29)

The one from (29), by considering that the ω speed has small values, this it can be neglected.

From the (21) equation, it can be obtained:

$${\displaystyle \frac{\partial \overrightarrow{v}}{\partial \stackrel{•}{q}}}={\left[N\right]}^T\left[A_{o,n}\right]\hspace{0.17em}\left\{{\overrightarrow{w}}_0\right\}$$

(30)

Based on the Eqns (27) and (23), we will obtain the generalized forces attached to the inertial ones, as it follows:

$$\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{Q}^{*}=\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\hspace{0.17em}m\left(\begin{array}{l}{\displaystyle \sum _{i=1}^n\left({\left\{\stackrel{••}{p_i}\right\}}^T\right.\left[A_{0,i-1}\right]}\hspace{0.17em}{\left[A_{0,n}\right]}^T\left[N\right]+\\ +2{\left\{{\stackrel{•}{p}}_i\right\}}^T\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}{\left[A_{0,n}\right]}^T\left[N\right]+\\ \hspace{0.17em}+{\left\{p_i\right\}}^T\hspace{0.17em}\left[{\tilde{\epsilon }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}{\left[A_{0,n}\right]}^T\left[N\right]+\\ +{\left\{p_i\right\}}^T\hspace{0.17em}\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[{\tilde{\omega }}_{0,i-1}\right]\hspace{0.17em}\left[A_{0,i-1}\right]\hspace{0.17em}{\left[A_{0,n}\right]}^T\left.\left[N\right]\right)+\\ +{\left\{\stackrel{••}{q}\right\}}^T{\left[N\right]}^T\left[N\right]+2{\left\{\stackrel{•}{q}\right\}}^T{\left[N\right]}^T\left[{\tilde{\omega }}_{o,n}\right]\left[N\right]+\\ +{\left\{r_n\right\}}^T\hspace{0.17em}\hspace{0.17em}\left[{\tilde{\epsilon }}_{0,n}\right]\hspace{0.17em}\hspace{0.17em}\left[N\right]+{\left\{r_n\right\}}^T\hspace{0.17em}\left[{\tilde{\omega }}_{0,n}\right]\hspace{0.17em}\left[{\tilde{\omega }}_{0,n}\right]\hspace{0.17em}\hspace{0.17em}\left[N\right]+\\ +{\left\{q\right\}}^T{\left[N\right]}^T\hspace{0.17em}\left[{\tilde{\omega }}_{0,n}\right]\hspace{0.17em}\left[{\tilde{\omega }}_{0,n}\right]\hspace{0.17em}\hspace{0.17em}\left[N\right]\end{array}\right)$$

(31)

The generalized force attached to the exterior forces, by supplementary considering the stiffness and dumping is:

$$Q=F_e-F_k-F_a$$

(32)

The exterior forces and the exterior torque – T, which actuates on the considered element:

$$F_e=\overrightarrow{F}{\displaystyle \frac{\partial \overrightarrow{v}}{\partial \stackrel{•}{q}}}+\overrightarrow{T}{\displaystyle \frac{\partial \overrightarrow{\omega }}{\partial \stackrel{•}{q}}}=\overrightarrow{F}{\displaystyle \frac{\partial \overrightarrow{v}}{\partial \stackrel{•}{q}}}$$

(33)

Thus, we consider the angular speed ω at small values and it will result: ${\displaystyle \frac{\partial \overrightarrow{\omega }}{\partial \stackrel{•}{q}}}=0$By considering the F – force actuates on the “n” element, it can be written:

$$\overrightarrow{F}={\left\{F\right\}}^T\left\{\overrightarrow{w_n}\right\}={\left\{F\right\}}^T\left[A_{o,n}\right]\left\{\overrightarrow{w_0}\right\}$$

(34)

$$\overrightarrow{F_e}=\overrightarrow{F}{\displaystyle \frac{\partial \overrightarrow{v}}{\partial \stackrel{•}{q}}}={\left\{F\right\}}^T\left[N\right]$$

(35)

The force resulted from the “n” element stiffness is:

$$F_k={\left\{q\right\}}^T\left[K\right]\left[20\right]$$

(36)

where : [K] is the stiffness matrix, given by:

$$\left[K\right]={\left[B\right]}^T\left[D\right]\left[B\right]\left[20\right]$$

(37)

In (37), we identify: [B] is the connection matrix between deformations and displacements; [D] is the material properties matrix.

Thus, the force resulted from the dumping has the following form:

$$F_c={\left\{\stackrel{•}{q}\right\}}^T\left[C\right]$$

(38)

where: [C] is the dumping matrix given by:

$$\left[C\right]={\left[N\right]}^T\left[\mu \right]\hspace{0.17em}\left[N\right]\left[20\right]$$

(39)

where [μ] is the dumping coefficients matrix.

$$Q=F_e-F_k-F_c={\left\{F\right\}}^T\left[N\right]-{\left\{q\right\}}^T\left[K\right]-{\left\{\stackrel{•}{q}\right\}}^T\left[C\right]$$

(40)

By turning at the Equation (24), which means: Q+Q^*=0, and introducing the Equation (40) in (24) we obtain:

$$\begin{array}{l}{\left\{F\right\}}^T\left[N\right]-{\left\{q\right\}}^T\left[K\right]-{\left\{\stackrel{•}{q}\right\}}^T\left[C\right]-m\left\{{\displaystyle \sum _{i=1}^n\left({\left\{{\stackrel{••}{p}}_i\right\}}^T\right.}\left[A_{0,i-1}\right]{\left[A_{0,n}\right]}^T\left[N\right]\right.+\\ +2{\left\{{\stackrel{•}{p}}_i\right\}}^T\left[{\tilde{\omega }}_{0,i-1}\right]\left[A_{0,i-1}\right]{\left[A_{0,n}\right]}^T\left[N\right]+{\left\{p_i\right\}}^T\left[{\tilde{\epsilon }}_{0,i-1}\right]\left[A_{0,i-1}\right]{\left[A_{0,n}\right]}^T\left[N\right]+\\ \left.+{\left\{p_i\right\}}^T\left[{\tilde{\omega }}_{0,i-1}\right]\left[{\tilde{\omega }}_{0,i-1}\right]\left[A_{0,i-1}\right]{\left[A_{0,n}\right]}^T\left[N\right]\right)+{\left\{\stackrel{••}{q}\right\}}^T{\left[N\right]}^T\left[N\right]+\\ +2{\left\{\stackrel{•}{q}\right\}}^T{\left[N\right]}^T\left[{\tilde{\omega }}_{0,n}\right]\left[N\right]+{\left\{r_n\right\}}^T\left[{\tilde{\epsilon }}_{0,n}\right]\left[N\right]+{\left\{q\right\}}^T{\left[N\right]}^T\left[{\tilde{\epsilon }}_{0,n}\right]\left[N\right]+\\ +{\left\{r_n\right\}}^T\left[{\tilde{\omega }}_{0,n}\right]\left[{\tilde{\omega }}_{0,n}\right]\left[N\right]+{\left\{q\right\}}^T{\left[N\right]}^T\left[{\tilde{\omega }}_{0,n}\right]\left[{\tilde{\omega }}_{0,n}\right]\left[N\right]=0\end{array}$$

(41)

For a finite element “e” the motion is described by the following equation system:

$${\left\{\stackrel{••}{q}\right\}}^T\left[M_e\right]+{\left\{\stackrel{•}{q}\right\}}^T\left[C_e\right]+{\left\{q\right\}}^T\left[K_e\right]={\left[F_e\right]}^T\left[21\right]$$

(42)

From the Equation (41), we identify the terms coefficients ${\left\{\stackrel{••}{q}\right\}}^T$, ${\left\{\stackrel{•}{q}\right\}}^T$, ${\left\{q\right\}}^T$ namely $\left[M_e\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left[C_e\right]$, $\hspace{0.17em}\hspace{0.17em}\left[K_e\right]$ as it follows:

The mass matrix is:

$$\begin{array}{c}\left[M^{(e)}\right]={\displaystyle \underset{V^{(e)}}{\iiint }{\rho }_e{m}_e}{\left[N\right]}^T\left[N\right]dV\end{array}$$

(43)

The dumping matrix with both components, respectively static and dynamic is:

$$\begin{array}{c}\left[C^{(e)}\right]={\displaystyle \underset{V^{(e)}}{\iiint }{\left[N\right]}^T\left[q\right]\left[N\right]}dV+2{\displaystyle \underset{V^{(e)}}{\iiint }{\rho }_e{\left[N\right]}^T}\left(\left[{\tilde{\omega }}_{0n}\right]+\left[{\tilde{\omega }}_{0n}\right]\left[{\tilde{\omega }}_{0n}\right]\left[N\right]dV=\left[{C_s}^{(e)}\right]+\left[{C_d}^{(e)}\right]\right)\end{array}$$

(44)

The elemental nodal force can be determined with the following relation:

$$\begin{array}{c}\left\{F^{(e)}\right\}={\displaystyle \underset{V^{(e)}}{\iiint }{\left\{N\right\}}^T\left[F\right]}dV-{\displaystyle \underset{V^{(e)}}{\iiint }{\rho }_e{\left[N\right]}^T}\left\{\left(\left[{\tilde{\epsilon }}_{0,n}\right]+\left[{\tilde{\omega }}_{0,n}\right]\left[{\tilde{\omega }}_{0,n}\right]\right)\right.\left\{r_n\right\}+\\ +{\displaystyle \sum _{i=1}^n\left[A_{n-1}\right]\left(\left\{\stackrel{••}{p_i}\right\}\right.}+2\left[{\tilde{\omega }}_{0,i-1}\right]\left\{\stackrel{•}{p_i}\right\}+\left[{\tilde{\epsilon }}_{0,i-1}\right]\left\{p_i\right\}\left.+\left.\left[{\tilde{\omega }}_{0,i-1}\right]\left[{\tilde{\omega }}_{0,i-1}\right]\left\{p_i\right\}\right)\right\}dV\end{array}$$

(45)

The stiffness matrix with both components, respectively static and dynamic is:

$$\begin{array}{c}\left[K^{(e)}\right]={\displaystyle \underset{V^{(e)}}{\iiint }{\left[B\right]}^T\left[D\right]\left[B\right]}dV+{\displaystyle \underset{V^{(e)}}{\iiint }{\rho }_e{\left[N\right]}^T}\left(\left[{\tilde{\omega }}_{0n}\right]+\left[{\tilde{\omega }}_{0n}\right]\left[{\tilde{\omega }}_{0n}\right]\left[N\right]dV=\left[{K_s}^{(e)}\right]+\left[{K_d}^{(e)}\right]\right)\end{array}$$

(46)

3. Elastodynamic Analysis of the Mechanism from an Internal Combustion Engine

There will be analyzed, under kinematic and dynamic aspects, a mechanism from an internal combustion engine, shown in Figure 4. This mechanism consists of the following components: crankshaft 1, connecting rod 2, piston 3. The analyzed motion is transmitted from the crankshaft through connecting rods at the three in–line pistons. By knowing the functional role and the mechanism destination, it is well known the fact that the functional cycle involves rapid motions, with an accentuated dynamic character. Thus, it is necessary to develop a kinematic and kinetostatic calculus, which assumes the variation motion laws during the time for the kinematic parameters and the connection forces of the kinematic joints. The crankshaft rotational speed is 6000rpm.

The mechanism motion is monitored on a time interval equal to 0.04 seconds when the crankshaft develops two complete rotations.

Thus, by having in sight the rapid motions of the links that generate high inertial forces, there will be imposed a mechanism elastodynamic analysis. For this, the mechanism connecting rods is analyzed as deformable bodies.

For the inverse dynamic analysis processing, it is needed to know the crankshaft angular speed, which is ω=628.318 radians/second, crank slider length r= 40millimeters, connecting rod length l= 137 millimeters, mass characteristics of the mechanism links (m1= 0.952 kilograms, m2= 0.153 kilograms, and m3= 0.214 kilograms), and the gas pressure developed on the piston head. There will be determined the positions, speeds, and accelerations of A, and B characteristic points, respectively angular speed and angular acceleration of the AB connecting rod.

The force developed by the gas pressure onto the piston head is:

$$F_g=p_g{A}_p$$

(47)

where: p_g– gas pressure, established experimentally; A_p – piston head surface.

There were experimentally determined the gas pressure developed onto the piston head, depending on the crankshaft angle as can be observed in Figure 5.

In Figure 5, there are presented the main parameters variation of an internal combustion engine, depending on the crankshaft rotational angle (Crank Angle, deg). Thus, it can be remarked on the four distinctive phases of the engine stroke: intake, compression, power, and exhaust. The diagrams main elements reported in Figure 5 are as follows:

- Black curve (cylinder gas pressure): Two major peaks are visible, corresponding to successive combustion events. The maximum values exceed 90 bar, indicating the fuel expansion phase. During compression, pressure increases before ignition, followed by a rapid drop after the exhaust phase.

- Orange curve (fuel ignition): This shows the exact ignition point of the fuel, correlated with the peak pressure. A delay between injection and ignition is noticeable, which can be analyzed for optimizing engine performance.

- Green curve (reaction force in the piston, Fx): A sharp increase in force occurs at ignition, followed by a controlled decrease. This confirms that gas expansion generates a maximum force on the piston, influencing the mechanical power output.

- Red and blue curves (intake and exhaust valve lift): Intake valves (blue curve) open at the beginning of the intake stroke and close at its end. The exhaust valves (red curve) open at the end of the expansion stroke and remain open until the beginning of the intake stroke.

Also, in the superior part of Figure 5, the pressure variation diagram was developed onto the piston head, in bar as measuring unit, depending on the crankshaft rotational angle. According to the Figure 5, it was obtained the force developed by the gas pressure, by knowing the piston head surface, and this variation was adapted as a function of time represented in Figure 6. Thus, the mechanism was modeled through beam elements type, as it can be seen in Figure 7 and Figure 8.

Each mechanism kinematic link will be considered as two finite elements, beam type, with two nodes at each end, respectively 3DoF for each node (Figure 9 and Figure 10). As loads, there were considered the finite elements weights and the force developed by the gas pressure onto the piston head, namely F_b. For finite element no. 2, the piston component can be replaced with a concentrated mass – m_c.

The mechanical system motion can be reported at the following reference systems:

• R1(x1,y1,z1) – mobile reference system attached to the finite element no. 1;
• R2(x2,y2,z2) – mobile reference system attached to the finite element no. 2;
• R0(X0,Y0,Z0) – global reference system.

The proper mathematical model for ,,e’’ finite element displacements is:

$$\left\{\overrightarrow{U^{(e)}}\right\}=\left\{\begin{array}{c}u(x)\\ v(x)\end{array}\right\}=\left[N^{(e)}\right]\left\{{\overrightarrow{q}}^{(e)}\right\},$$

(48)

where: $u(x)$– axial displacements; $v(x)$,$\theta (x)$– transversal displacements and angular displacements. From Equation (48), it can be identified:

$$\left[N^{(e)}\right]=\left[\begin{array}{cccccc}{N}_{1x}& 0& 0& N_{2x}& 0& 0\\ 0& N_1& N_2& 0& N_3& N_4\end{array}\right]$$

(49)

$\left[N^{(e)}\right]$– is the interpolation functions matrix;

$$\left\{{\overrightarrow{q}}^{(e)}\right\}={\left[\begin{array}{cccccc}{q}_1& q_2& q_3& q_4& q_5& q_6\end{array}\right]}^T$$

(50)

$\left\{{\overrightarrow{q}}^{(e)}\right\}$– is the nodal displacements vector for a finite element,,e “ with e=1,2– finite elements.

In this vector structure it will have axial and transversal displacements, as follows:

$$\left\{\overrightarrow{q_x}\right\}={\left[\begin{array}{cc}{q}_1& q_4\end{array}\right]}^T–\mathrm{axial}\mathrm{nodal}\mathrm{displacements};$$

(51)

$$\left\{\overrightarrow{q_{xy}}\right\}={\left[\begin{array}{cc}{q}_2& q_5\end{array}\right]}^T–\mathrm{transversal}\mathrm{nodal}\mathrm{displacements};$$

(52)

$$\left\{\overrightarrow{q_{\theta }}\right\}={\left[\begin{array}{cc}{q}_3& q_6\end{array}\right]}^T–\mathrm{angular}\mathrm{nodal}\mathrm{displacements}.$$

(53)

For the axial displacements, there will be adopted the linear interpolating model, and the form function matrix components become:

$$N_{1x}=1-{\displaystyle \frac{x}{l}};N_{2x}={\displaystyle \frac{x}{l}};\left[21\right]$$

(54)

For the transversal displacements there will be adopted the cubic interpolating model, and the form function matrix components will be:

$$\begin{gathered} N_1={\displaystyle \frac{2x^3-3x^2+l^3}{l^3}};N_2={\displaystyle \frac{x^3-2l{x}^2+l^{2}x}{l^2}}; \\ {N}_3=-{\displaystyle \frac{2x^3-3l{x}^2}{l^3}};N_4={\displaystyle \frac{x^3-l{x}^2}{l^2}}\left[21\right] \end{gathered}$$

(55)

For the mechanism assembled structure will have:

$$\left\{\overrightarrow{q}\right\}={\left[\begin{array}{ccccccccc}{q}_1& q_2& q_3& q_4& q_5& q_6& q_7& q_8& q_9\end{array}\right]}^T$$

(56)

$\left\{\overrightarrow{q}\right\}$– nodal displacements vector in the local coordinates system.

The motion equations, in a matrix form for the mechanism, are:

$$\left[M\right]\left\{\stackrel{··}{{\displaystyle u}}\right\}+\left[C\right]\left\{\stackrel{·}{{\displaystyle u}}\right\}+\left[K\right]\left\{\stackrel{·}{{\displaystyle u}}\right\}=\left\{F\right\},$$

(57)

where:

$\left[M\right]$– mass matrix, assembled for the entire analyzed mechanical system; $\left[K\right]$– assembled stiffness matrix; $\left[C\right]$ – assembled dumping matrix; $\left[F\right]$ – nodal force vector, assembled for the entire mechanical system.

For the presented mechanical system, made of two finite elements (beam type) with 3DoF on each node, there will be determined the matrices from the motion equations structure.

3.1. Nodal Force Calculus

The nodal forces corresponding to the local coordinates system will be obtained by having in sight the Equation (45). Thus, the nodal force for the ,,e” finite element can be determined with the following relations:

– static components:

$$\left\{F^{(e)}\right\}={\displaystyle \underset{V^{(e)}}{\iiint }{\left\{N\right\}}^T\left[F\right]}dV$$

(58)

– dynamic components:

$$\begin{array}{c}\left\{{F_d}^{(e)}\right\}={\displaystyle \underset{V^{(e)}}{\iiint }{\rho }^{(e)}}\left\{{\left\{r_n\right\}}^T\left(\left[{\tilde{\epsilon }}_{0,n}\right]+\left[{\tilde{\omega }}_{0,n}\right]\left[{\tilde{\omega }}_{0,n}\right]\right)\right.\left[N\right]+\\ +{\displaystyle \sum _{i=1}^n\left({\left\{\stackrel{••}{p_i}\right\}}^T\right.}+2{\left\{\stackrel{•}{p_i}\right\}}^T\left[{\tilde{\omega }}_{0,i-1}\right]+{\left\{p_i\right\}}^T\left[{\tilde{\epsilon }}_{0,i-1}\right]\left.+{\left\{p_i\right\}}^T\left[{\tilde{\omega }}_{0,i-1}\right]\left[{\tilde{\omega }}_{0,i-1}\right]\left[A_{n-1}\right]\left[N\right]\right\}dV\end{array}$$

(59)

Finite element no.1:

The nodal force static component, which corresponds to finite element no. 1, is:

$$F_s^{(1)}={\rho }^{(1)}g{A}^{(1)}{\displaystyle \underset{0}{\overset{l_1}{\int }}{\left[N_1\right]}^{T}dx},$$

(60)

where:

${\rho }^{(e)}$ – material density; g – gravitational acceleration; $A^{(e)}$ – transversal section of the finite element „e” with e=1..2

$$N_1=\left[0\begin{array}{cccccc}1-{\displaystyle \frac{x}{l_1}}& 0& 0& {\displaystyle \frac{x}{l_1}}& 0& 0\\ 0& {\displaystyle \frac{2x^3-3l_1{x}^2+l_1^3}{l_1^3}}& {\displaystyle \frac{x^3-2l_1{x}^2+l_1^{2}x}{l_1^2}}& 0& -{\displaystyle \frac{2x^3-3l_1{x}^2}{l_1^3}}& {\displaystyle \frac{x^3-l_1{x}^2}{l_1^2}}\\ 0& 0& 0& 0& 0& 0\end{array}\right]$$

(61)

The nodal force static component expressed by taking into account the local reference system is:

$$f_s^{(1)}=F_s^{(1)}\left[\begin{array}{c}sin{\phi }_1\\ -\cos {\phi }_1\\ 0\end{array}\right]$$

(62)

The nodal force dynamic component is:

$$\left[F_d^{\left(1\right)}\right]={\rho }^{\left(1\right)}{A}^{(1)}{\displaystyle {\int }_0^{l1}\left\{{\left\{r_1\right\}}^T\left[\stackrel{\sim }{{\omega }_{01}}\right]\cdot \left[\stackrel{\sim }{{\omega }_{01}}\right]\left[N_1\right]\right\}dx}$$

(63)

The angular speed can be written as a anti-symmetrical matrix as follows:

$$\left[\stackrel{\sim }{{\omega }_{01}}\right]=\left[\begin{array}{ccc}0& -{\omega }_{01}^z& {\omega }_{01}^y\\ {\omega }_{01}^z& 0& -{\omega }_{01}^x\\ -{\omega }_{01}^y& {\omega }_{01}^x& 0\end{array}\right]$$

(64)

The force which actuates on the finite element no. 1, in the local reference system is:

$$\left[f_1\right]=\left[f_s^{\left(1\right)}\right]+\left[F_d^{\left(1\right)}\right],$$

(65)

where: $\left[T_e\right]$ – is the transformation matrix of the coordinates for crossing from a reference system attached to finite element „e”, to a global reference system (e=1,2):

$$\left[T_e\right]=\left[\begin{array}{cccccc}\cos {\phi }_e& \sin {\phi }_e& 0& 0& 0& 0\\ -\sin {\phi }_e& \cos {\phi }_e& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& \cos {\phi }_e& \sin {\phi }_e& 0\\ 0& 0& 0& -\sin {\phi }_e& \cos {\phi }_e& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$$

(66)

The force that actuates onto finite element no. 1 in the global reference system is:

$$\left[F_1\right]={\left[T_1\right]}^T\cdot \left[f_1\right],$$

(67)

where: $\left[T_1\right]$ – is the coordinates transformation matrix, for crossing from a reference system attached to finite element no. 1, to the global reference system.

Finite element no. 2:

The static components can be obtained in accordance with the global reference system by considering the nodal mass force as:

$$\begin{gathered} \left[{F_s}^{(2)}\right]={\displaystyle \underset{V^{(2)}}{\iiint }{\left[N_2\right]}^T}{G}_{2}dv \\ \left[{F_s}^{(2)}\right]={\rho }^{(2)}g{A}^{(2)}{\displaystyle \underset{0}{\overset{l_2}{\int }}{\left[N_2\right]}^{T}dx}, \end{gathered}$$

(68)

where: N₂ – is the form functions matrix which corresponds to the finite element no. 2.

The nodal force static component, depending on the local reference system is:

$$\left[f_s^{(2)}\right]=\left[{F_s}^{(2)}\right]\left[\begin{array}{c}-\sin {\phi }_2\\ \cos {\phi }_2\\ 0\end{array}\right]+f_g$$

(69)

The produced force by the gas pressure in the local reference system, which actuates in node no. 3, is given by:

$$f_g={\left[\begin{array}{cccccc}0& 0& 0& -F_g\sin {\phi }_2& F_g\cos {\phi }_2& 0\end{array}\right]}^T$$

(70)

The dynamic components are:

$$\begin{array}{c}\left\{F_d^{(2)}\right\}=-{\displaystyle \underset{V^{(2)}}{\iiint }{\rho }^{(2)}{\left[N_2\right]}^T}\left\{\left(\left[{\tilde{\epsilon }}_{02}\right]+\left[{\tilde{\omega }}_{02}\right]\left[{\tilde{\omega }}_{02}\right]\right)\right.\left\{r_2\right\}+\\ +{\displaystyle \sum _{i=1}^n\left[A_{i-1,i}\right]\left(\left\{\stackrel{••}{p_i}\right\}\right.}+2\left[{\tilde{\omega }}_{0i-1}\right]\left\{\stackrel{•}{p_i}\right\}+\left[{\tilde{\epsilon }}_{0i-1}\right]\left\{p_i\right\}\left.+\left.\left[{\tilde{\omega }}_{0i-1}\right]\left[{\tilde{\omega }}_{0i-1}\right]\left\{p_i\right\}\right)\right\}dV\end{array}$$

(71)

$$\left[F_d^{(2)}\right]=-{\displaystyle \underset{V^{(l)}}{\iiint }{\rho }^{(2)}{\left[N_2\right]}^T\left\{\left(\begin{array}{l}\left[\stackrel{\sim }{{\epsilon }_{02}}\right]+\left[\stackrel{\sim }{{\omega }_{02}}\right]\cdot \left[\stackrel{\sim }{{\omega }_{02}}\right]\cdot \left\{r_2\right\}+\\ +\left[A_{12}\right]\cdot \left(\begin{array}{l}\left\{{\stackrel{··}{p}}_2\right\}+2\left[{\tilde{\omega }}_{01}\right]\left\{\stackrel{•}{p_2}\right\}+\\ +\left[\stackrel{\sim }{{\epsilon }_{01}}\right]+\left[\stackrel{\sim }{{\omega }_{01}}\right]\cdot \left[\stackrel{\sim }{{\omega }_{01}}\right]\end{array}\right)\left\{p_2\right\}\end{array}\right)\right\}}dv,$$

(72)

where:

$\left\{r_2\right\}-$ Is the B point vector in the local coordinates system, namely R₂(x₂,y₂,z₂):

$$\left\{r_2\right\}={\left[\begin{array}{ccc}{l}_2& 0& 0\end{array}\right]}^T$$

(73)

$\left\{p_2\right\}-$Is the generalized coordinates vector for the finite element no. 2 (mass center position and finite element no. 2 orientation):

$$\left\{p_2\right\}={\left[\begin{array}{ccc}{x}_2& y_2& {\phi }_2\end{array}\right]}^T$$

(74)

$\left\{{\stackrel{·}{p}}_2\right\}-$ Is the vector of generalized speeds:

$$\left\{p_2\right\}={\left[\begin{array}{ccc}{\stackrel{·}{x}}_2& {\stackrel{·}{y}}_2& {\stackrel{·}{\phi }}_2\end{array}\right]}^T$$

(75)

$\left\{{\stackrel{··}{p}}_2\right\}-$Is the vector of generalized accelerations:

$$\left\{{\stackrel{··}{p}}_2\right\}={\left[\begin{array}{ccc}{\stackrel{··}{x}}_2& {\stackrel{··}{y}}_2& {\stackrel{··}{\phi }}_2\end{array}\right]}^T$$

(76)

Point C₂ position is:

$$\left\{\begin{array}{l}{x}_2=l_1\cos {\phi }_1+{\displaystyle \frac{l_2}{2}}\cos {\phi }_2\\ y_2=l_1\sin {\phi }_1+{\displaystyle \frac{l_2}{2}}\sin {\phi }_2\end{array}\right.$$

(77)

Point C₂ speed is:

$$\left\{\begin{array}{l}{\stackrel{·}{x}}_2=-l_1\stackrel{·}{{\phi }_1}\sin {\phi }_1+{\displaystyle \frac{l_2}{2}}\stackrel{·}{{\phi }_2}sin{\phi }_2\\ {\stackrel{·}{y}}_2=l_1\stackrel{·}{{\phi }_1}\cos {\phi }_1+{\displaystyle \frac{l_2}{2}}\stackrel{·}{{\phi }_2}\cos {\phi }_2\end{array}\right.$$

(78)

Point C₂ acceleration is:

$$\left\{\begin{array}{l}{\stackrel{··}{x}}_2=-l_1{\epsilon }_1\sin {\phi }_1-l_1{\omega }_1^2\cos {\phi }_1-{\displaystyle \frac{l_2}{2}}{\epsilon }_2\sin {\phi }_2-{\displaystyle \frac{l_2}{2}}{\omega }_2^2\cos {\phi }_2\\ {\stackrel{··}{y}}_2=l_1{\epsilon }_1\cos {\phi }_1-l_1{\omega }_1^2\sin {\phi }_1+{\displaystyle \frac{l_2}{2}}{\epsilon }_2\cos {\phi }_2-{\displaystyle \frac{l_2}{2}}{\omega }_2^2\sin {\phi }_2\end{array}\right.$$

(79)

The angular speed of the finite element “i” is:

$$\left\{{\omega }_i\right\}={\displaystyle \frac{d{\phi }_i}{dt}}={\stackrel{·}{\phi }}_i;\mathrm{i}=1,2$$

(80)

The angular acceleration of the finite element “i” is:

$$\left\{{\epsilon }_i\right\}={\stackrel{··}{\phi }}_i;\mathrm{i}=1,2$$

(81)

The coordinates transformation matrix for crossing from a local axis reference system R_i(x_i,y_i,z_i) to a global reference system R₀(X₀,Y₀,Z₀) is:

$$\left[A_{0i}\right]=\left[\begin{array}{ccc}\cos ({\phi }_i)& -\sin ({\phi }_i)& 0\\ \sin ({\phi }_i)& \cos ({\phi }_i)& 0\\ 0& 0& 1\end{array}\right];\mathrm{i}=1,2$$

(82)

Thus, in this case, the coordinates transformation matrix for crossing from a local axis reference system R₁(x₁,y₁,z₁) to a global reference system R₂(x₂,y₂,z₂) is:

$$\left[A_{12}\right]=\left[\begin{array}{ccc}\cos ({\phi }_1-{\phi }_2)& -\sin ({\phi }_1-{\phi }_2)& 0\\ \sin ({\phi }_1-{\phi }_2)& \cos ({\phi }_1-{\phi }_2)& 0\\ 0& 0& 1\end{array}\right]$$

(83)

The force which actuates on the finite element no. 2, in the local axis system is:

$$\left[f_2\right]=\left[f_s^{(2)}\right]+\left[F_d^{(2)}\right]$$

(84)

The force which actuates on the finite element no. 2, in the global axis system is:

$$\left[F_2\right]={\left[T_2\right]}^T\cdot \left[f_2\right],$$

(85)

where: $\left[T_2\right]$ – is the coordinates transformation matrix, for crossing from a reference system attached to finite element no. 1, to a global reference system.

The assembled nodal force for the engine mechanism, according to the global reference system, is:

$$\left[F\right]={\left[\begin{array}{ccccccccc}{F}_{11}^{(1)}& F_{21}^{(1)}& F_{31}^{(1)}& F_{41}^{(1)}+F_{11}^{(2)}& F_{51}^{(1)}+F_{21}^{(2)}& F_{61}^{(1)}+F_{31}^{(2)}& F_{41}^{(2)}& F_{51}^{(2)}& F_{61}^{(2)}\end{array}\right]}^T$$

(86)

3.2. Stiffness Matrix

$$\begin{array}{l}\left[K^{(e)}\right]={\displaystyle \underset{V^{(e)}}{\iiint }{\left[B\right]}^T\left[D\right]\hspace{0.17em}\left[B\right]}dV+{\displaystyle \underset{V^{(e)}}{\iiint }{\rho }_e}{\left[N\right]}^T\left(\left[{\tilde{\omega }}_{0n}\right]\left[{\tilde{\omega }}_{0n}\right]+\left[{\epsilon }_{0,n}\right]\right)\left[N\right]dV=\\ =\left[K_s^{(e)}\right]+\left[K_d^{(e)}\right]\end{array},$$

(87)

where: $\left[K_s^{(e)}\right]={\displaystyle \underset{V^{(e)}}{\iiint }{\left[B\right]}^T\left[D\right]\left[B\right]}dV$– represents the static component; $\left[{K_d}^{(e)}\right]={\displaystyle \underset{V^{(e)}}{\iiint }{\rho }_e}{\left[N\right]}^T\left(\left[{\tilde{\omega }}_{0n}\right]\left[{\tilde{\omega }}_{0n}\right]+\left[\stackrel{\sim }{{\epsilon }_{0,n}}\right]\right)\left[N\right]dV$– represents the dynamic component; $\left[D\right]=E$– material elasticity modulus.

Finite element no. 1:

Static component. The stiffness matrix in a local coordinate system is:

$$\left[K_s^{\left(l_1\right)}\right]={\displaystyle \underset{v^{\left(1\right)}}{\iiint }{\left[B\right]}^T\left[D\right]}\left[B\right]dv=E^{(1)}{A}^{(1)}{\displaystyle {\int }_0^{l_1}{\left[B\right]}^T\left[B\right]dx},$$

(88)

where: $\left[B\right]=\left[\begin{array}{cccc}{\displaystyle \frac{12x-6l_1}{{l_1}^3}}& {\displaystyle \frac{6x-4l_1}{{l_1}^2}}& -{\displaystyle \frac{12x-6l_1}{{l_1}^3}}& {\displaystyle \frac{6x-2l_1}{{l_1}^2}}\end{array}\right]$– for transversal displacements; $\left[B\right]=\left[\begin{array}{cc}-{\displaystyle \frac{1}{l_1}}& {\displaystyle \frac{1}{l_1}}\end{array}\right]$ for axial displacements.

The stiffness matrix for the transversal displacements:

$$\left[K_{ts}^{\left(l_1\right)}\right]={\displaystyle \frac{E{I}_{z1}}{{l_1}^3}}\left[\begin{array}{cccc}12& 6l_1& -12& 6l_1\\ 6l_1& 4l^2& -6l& 2{l_1}^2\\ -12& -6l_1& 12& -6l_1\\ 6l_1& 2{l_1}^2& -6l_1& 4{l_1}^2\end{array}\right]$$

(89)

The stiffness matrix for the axial displacements:

$$\left[K_{as}^{\left(l_1\right)}\right]={{\displaystyle \frac{E^{(1)}A}{l_1}}}^{(1)}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$$

(90)

Dynamic component: For the finite element no. 1, in the local axis reference system:

$$\begin{array}{l}\left[K_d^{\left(l_1\right)}\right]={\displaystyle \underset{v^{\left(e\right)}}{\iiint }{\rho }^{\left(1\right)}}{\left[N_1\right]}^T\left(\left[\stackrel{•}{\stackrel{~}{{\omega }_{01}}}\right]+\left[\widetilde{{\omega }_{01}}\right]\left[\widetilde{{\omega }_{01}}\right]\right)\left[N_1\right]dv\\ \hspace{0.17em}\end{array}$$

(91)

The stiffness matrix for the finite element no. 1, according to the local reference system is:

$$\left[K^{\left(l_1\right)}\right]=\left[K_s^{\left(l_1\right)}\right]+\left[K_d^{\left(l_1\right)}\right]$$

(92)

The stiffness matrix for the finite element no. 1, according to the global reference system is:

$$\left[K^{(1)}\right]={\left[T\right]}^T\left[K^{(l_1)}\right]\left[T\right]$$

(93)

Finite element no. 2:

The static component of the stiffness matrix in the local axis reference system is:

$$\left[K_s^{\left(l_2\right)}\right]={\displaystyle \underset{v^{\left(2\right)}}{\iiint }{\left[B\right]}^T\left[D\right]}\left[B\right]dv=E^{(2)}{A}^{(2)}{\displaystyle {\int }_0^{l_2}{\left[B\right]}^T\left[B\right]dx}$$

(94)

The stiffness matrix which corresponds to transversal displacements is:

$$\left[K_{ts}^{\left(l_2\right)}\right]={\displaystyle \frac{E{I}_{z2}}{{l_1}^3}}\left[\begin{array}{cccc}12& 6l_2& -12& 6l_2\\ 6l_2& 4{l_2}^2& -6l_2& 2{l_2}^2\\ -12& -6l_2& 12& -6l_2\\ 6l_2& 2{l_2}^2& -6l_2& 4{l_2}^2\end{array}\right]$$

(95)

The stiffness matrix which corresponds to axial displacements is:

$$\left[K_{as}^{\left(l_2\right)}\right]={{\displaystyle \frac{E^{(2)}A}{l_2}}}^{(2)}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$$

(96)

The assembled stiffness matrix for the finite element no. 2, according to the local reference axis system:

$$\left[{K_s}^{(l2)}\right]=\left[\begin{array}{cccccc}{\displaystyle \frac{A^{(2)}{E}^{(2)}}{l_1}}& 0& 0& -{\displaystyle \frac{A^{(2)}{E}^{(2)}}{l_1}}& 0& 0\\ 0& {\displaystyle \frac{12E^{(2)}{I}_{z2}}{{l_2}^3}}& {\displaystyle \frac{6E^{(2)}{I}_{z2}}{{l_2}^2}}& 0& -{\displaystyle \frac{12E^{(2)}{I}_{z2}}{{l_2}^3}}& {\displaystyle \frac{6E^{(2)}{I}_{z2}}{{l_2}^2}}\\ 0& {\displaystyle \frac{6E^{(2)}{I}_{z2}}{{l_2}^2}}& {\displaystyle \frac{4E^{(2)}{I}_{z2}}{l_2}}& 0& -{\displaystyle \frac{6E^{(2)}{I}_{z2}}{{l_2}^2}}& 2{\displaystyle \frac{E^{(2)}{I}_{z2}}{l_2}}\\ -{\displaystyle \frac{A^{(2)}{E}^{(2)}}{l_2}}0& 0& 0& {\displaystyle \frac{A^{(2)}{E}^{(2)}}{l_2}}& 0& 0\\ 0& -{\displaystyle \frac{12E^{(2)}{I}_{z2}}{{l_2}^3}}& -{\displaystyle \frac{6E^{(2)}{I}_{z2}}{{l_2}^2}}& 0& {\displaystyle \frac{12E^{(2)}{I}_{z2}}{{l_2}^3}}& -{\displaystyle \frac{6E^{(2)}{I}_{z2}}{{l_2}^2}}\\ 0& {\displaystyle \frac{6E^{(2)}{I}_{z2}}{{l_2}^2}}& {\displaystyle \frac{2E^{(2)}{I}_{z2}}{l_2}}& 0& -{\displaystyle \frac{6E^{(2)}{I}_{z2}}{{l_2}^2}}& {\displaystyle \frac{4E^{(2)}{I}_{z2}}{l_2}}\end{array}\right]$$

(97)

The dynamic component of the stiffness matrix in a local axis reference system is:

$$\left[K_d^{\left(l2\right)}\right]={\rho }^{\left(2\right)}{A}^{(2)}{\displaystyle \underset{0}{\overset{l_2}{\int }}{\left[N_2\right]}^T\left[\begin{array}{ccc}-{{\omega }_2}^2& -{\epsilon }_2& 0\\ {\epsilon }_2& -{{\omega }_2}^2& 0\\ 0& 0& 0\end{array}\right]\left[N_2\right]dx}\hspace{0.17em}\hspace{0.17em}$$

(98)

The stiffness matrix for the finite element no. 2, according to the local reference system is:

$$\left[K^{\left(l_2\right)}\right]=\left[K_s^{\left(l_2\right)}\right]+\left[K_d^{\left(l_2\right)}\right]$$

(99)

The stiffness matrix for the finite element no. 2, according to the global reference system is:

$$\left[K^{(2)}\right]={\left[T\right]}^T\left[K^{(l_2)}\right]\left[T\right]$$

(100)

3.3. Mass Matrix

The mass matrix in the local coordinates system, can be written as follows:

$$M^{\left(e\right)}={\displaystyle \underset{V^{\left(e\right)}}{\iiint }{\rho }^{\left(e\right)}}{\left[N\right]}^T\left[N\right]dv={\rho }^{\left(e\right)}{A}^{(e)}{{\displaystyle {\int }_0^l\left[N\right]}}^T\left[N\right]dx,$$

(101)

where: ${\rho }^{\left(e\right)}$ – ,,e” finite element mass; $\left[N\right]$ – form function matrix (interpolating ones).

Finite element no. 1:

$$M^{\left(l_1\right)}={\displaystyle \underset{V^{\left(e\right)}}{\iiint }{\rho }^{\left(1\right)}}{\left[N_1\right]}^T\left[N_1\right]dv={\rho }^{\left(e\right)}{A}^{(e)}{{\displaystyle {\int }_0^l\left[N_1\right]}}^T\left[N_1\right]dx$$

(102)

Finite element no. 2:

$$\begin{array}{l}{M}^{\left(l_2\right)}={\displaystyle \underset{V^{\left(2\right)}}{\iiint }{\rho }^{\left(2\right)}}{\left[N_2\right]}^T\left[N_2\right]dv={\rho }^{\left(2\right)}{A}^{(2)}{{\displaystyle {\int }_0^{l_2}\left[N_2\right]}}^T\left[N_2\right]dx+\left[m_p\right]\end{array},$$

(103)

where: m_p is the piston mass.

The piston was considered a concentrated mass, which was attached to the finite element no. 2, respectively in node no. 3.

The mass matrix in a global coordinates system is:

Finite element no. 1:

$$\left[M^{(1)}\right]={\left[T_1\right]}^T\left[M^{(l_1)}\right]\left[T_1\right]$$

(104)

Finite element no. 2:

$$\left[M^{(2)}\right]={\left[T_2\right]}^T\left[M^{(l_2)}\right]\left[T_2\right]$$

(105)

The equilibrium equation expressed in the local coordinates system has the following form:

$$k\overrightarrow{q}=\overrightarrow{f}$$

(106)

Also: $\overrightarrow{Q}=T\overrightarrow{q}$ – are the nodal displacements in the global coordinates system; $\overrightarrow{F}=T\overrightarrow{f}$ – represent the nodal forces in a global coordinates system.

The assembled equilibrium equations, in a global coordinates system have the following form:

$$\begin{array}{l}\left[K\right]\left[Q\right]=\left[F\right]\\ \left[Q\right]={\left[K\right]}^T\left[F\right]\end{array},$$

(107)

where: $\left[K\right]$– is the assembled stiffness matrix; $\left[Q\right]$ – is the assembled nodal displacements vector; $\left[F\right]$ – represents the assembled nodal forces vector.

3.4. The Assembly of the Equations and Introducing the Contour Conditions

There will be introduced contour conditions (node no. 1 has a single DoF, which allows only the Q₃ rotation, and node no. 3, has two DoFs which for this was canceled the Q₈ displacement):

Finite element no. 1:

$$x=0\to \left\{\begin{array}{l}{Q}_1=0\\ Q_2=0\end{array}\right.$$

(108)

Finite element no. 2:

$$x=l_2\to Q_8=0$$

(109)

By considering the contour conditions, from the assembled matrix which interact in the motion equation, respectively$\left[M\right],\left[K\right],\left[C\right]$, there will be canceled line 1, column 1, line 2, column 2, and line 8, column 8. Thus, in these conditions, the remaining matrices will have 6x6 dimensions.

The reduced mass matrix has the following form:

$$\left[M_r\right]\left[\begin{array}{cccccc}{M}_{33}^{(1)}& M_{34}^{(1)}& M_{35}^{(1)}& M_{36}^{(1)}& 0& 0\\ M_{43}^{(1)}& M_{44}^{(1)}+M_{11}^{(2)}& M_{45}^{(1)}+M_{12}^{(2)}& M_{46}^{(1)}+M_{13}^{(2)}& M_{14}^{(2)}& M_{16}^{(2)}\\ M_{53}^{(1)}& M_{54}^{(1)}+M_{21}^{(2)}& M_{55}^{(1)}+M_{22}^{(2)}& M_{56}^{(1)}+M_{23}^{(2)}& M_{24}^{(2)}& M_{26}^{(2)}\\ M_{63}^{(1)}& M_{64}^{(1)}+M_{31}^{(2)}& M_{65}^{(1)}+M_{32}^{(2)}& M_{66}^{(1)}+M_{33}^{(2)}& M_{34}^{(2)}& M_{36}^{(2)}\\ 0& M_{41}^{(2)}& M_{42}^{(2)}& M_{43}^{(2)}& M_{44}^{(2)}& M_{46}^{(2)}\\ 0& M_{61}^{(2)}& M_{62}^{(2)}& M_{63}^{(2)}& M_{64}^{(2)}& M_{66}^{(2)}\end{array}\right]$$

(110)

The reduced stiffness matrix has the following form:

$$\left[K_r\right]\left[\begin{array}{cccccc}{K}_{33}^{(1)}& K_{34}^{(1)}& K_{35}^{(1)}& K_{36}^{(1)}& 0& 0\\ K_{43}^{(1)}& K_{44}^{(1)}+K_{11}^{(2)}& K_{45}^{(1)}+K_{12}^{(2)}& K_{46}^{(1)}+K_{13}^{(2)}& K_{14}^{(2)}& K_{16}^{(2)}\\ K_{53}^{(1)}& K_{54}^{(1)}+K_{21}^{(2)}& K_{55}^{(1)}+K_{22}^{(2)}& K_{56}^{(1)}+K_{23}^{(2)}& K_{24}^{(2)}& K_{26}^{(2)}\\ K_{63}^{(1)}& K_{64}^{(1)}+K_{31}^{(2)}& K_{65}^{(1)}+K_{32}^{(2)}& K_{66}^{(1)}+K_{33}^{(2)}& K_{34}^{(2)}& K_{36}^{(2)}\\ 0& K_{41}^{(2)}& K_{42}^{(2)}& K_{43}^{(2)}& K_{44}^{(2)}& K_{46}^{(2)}\\ 0& K_{61}^{(2)}& K_{62}^{(2)}& K_{63}^{(2)}& K_{64}^{(2)}& K_{66}^{(2)}\end{array}\right]$$

(111)

The nodal force assembled for the considered engine mechanism, by having in sight the global reference system, is:

$$\left[F_r\right]\left[\begin{array}{c}{F}_{31}^1\\ F_{41}^1+F_{11}^2\\ F_{51}^1+F_{21}^2\\ F_{61}^1+F_{31}^2\\ F_{41}^2\\ F_{61}^2\end{array}\right]$$

(112)

The elastic displacements vector for the element no. 1 is:

$$\left[\overrightarrow{U^{(1)}}\right]=\left[N_1\right]\left[\overrightarrow{Q^{(1)}}\right]$$

(113)

The elastic displacements vector for the element no. 2 has the following form:

$$\left[\overrightarrow{U^{(2)}}\right]=\left[N_2\right]\left[\overrightarrow{Q^{(2)}}\right]$$

(114)

or:

$$\left[\overrightarrow{U}\right]=\left[N\right]\left[\overrightarrow{Q}\right],$$

(115)

where:

$$\left[N\right]=\left[\begin{array}{cccccc}1-{\displaystyle \frac{x}{l}}& 0& 0& {\displaystyle \frac{x}{l}}& 0& 0\\ 0& {\displaystyle \frac{2x^3-3x^2+l^3}{l^3}}& {\displaystyle \frac{x^3-2l{x}^2+l^{2}x}{l^2}}& 0& -{\displaystyle \frac{2x^3-3l{x}^2}{l^3}}& {\displaystyle \frac{x^3-l{x}^2}{l^2}}\\ 0& 0& 0& 0& 0& 0\end{array}\right]$$

(116)

$$\left[\overrightarrow{Q}\right]=\left[K_r^T\right]\left[\overrightarrow{F_r}\right]$$

(117)

The nodal displacements vector for the finite element no. 1, in the global reference system, is:

$$\left[\overrightarrow{Q^{(1)}}\right]={\left[\begin{array}{cccccc}{Q}_1& Q_2& Q_3& Q_4& Q_5& Q_6\end{array}\right]}^T,$$

(118)

where: Q1=0; Q2=0 (by introducing the contour conditions).

The nodal displacement vector for the finite element no. 2, in the global reference system has the following form:

$$\left[\overrightarrow{Q^{(2)}}\right]={\left[\begin{array}{cccccc}{Q}_4& Q_5& Q_6& Q_7& Q_8& Q_9\end{array}\right]}^T$$

(119)

where: Q7=0, (by introducing the contour conditions).

3.5. Modal Dynamic Analysis

The motion equations are:

$$\left[M\right]\left\{\stackrel{··}{{\displaystyle q}}\right\}+\left[C\right]\left\{\stackrel{·}{{\displaystyle q}}\right\}+\left[K\right]\left\{q\right\}=\left\{F\right\}$$

(120)

The motion equations, without dumping, are:

$$\left[M\right]\left\{\stackrel{··}{q}\right\}+\left[K\right]\left\{q\right\}=\left\{F\right\}$$

(121)

The solution of the differential Equation (121), where the coefficients are considered constant ones, has the following form:

$$\left\{q\right\}=\left\{x\right\}\sin \left({\omega }_t+\phi \right)$$

(122)

By differentiating the above equation will be obtained:

$$\left\{\stackrel{·}{q}\right\}=\left\{x\right\}\omega \cdot \cos \left({\omega }_t+\phi \right)$$

(123)

$$\left\{\stackrel{··}{q}\right\}=-\left\{x\right\}{\omega }^2\cdot \sin \left({\omega }_t+\phi \right)$$

(124)

The differential equation without dumping and external loads will be:

$$\left[M\right]\left\{\stackrel{··}{q}\right\}+\left[K\right]\left\{q\right\}=0$$

(125)

$$\left(\left[K\right]-{\omega }^2\left[M\right]\right)\left\{x\right\}\cos \left({\omega }_t+\phi \right)=0$$

(126)

There will be obtained:

$$\left(\left[K\right]-{\omega }^2\left[M\right]\right)\left\{x\right\}=0$$

(127)

$${\left[K\right]}^{-1}\left[M\right]\left\{x\right\}=\left({\displaystyle \frac{1}{{\omega }^2}}\right)\left\{x\right\}$$

(128)

This is a proper value problem for the matrix ${\left[K\right]}^{-1}\left[M\right]$.

For a system with ,,n” degrees of freedom we have ,,n” proper values ,,λ” and ,,n” proper vector sets $\left\{{\varphi }_j\right\},j=1,\overline{n}$.

For each ,,j” we have:

$${\left[K\right]}^{-1}\left[M\right]\left\{{\varphi }_j\right\}={\lambda }_j\left\{{\varphi }_j\right\}=\left({\displaystyle \frac{1}{{\omega }^2}}\right)\left\{{\varphi }_j\right\},$$

(129)

where: ${\omega }_j$ = is the natural circular frequency:

$${\omega }_j={\displaystyle \frac{\pi n_j}{30}}\left[rot/\min \right]\left[rad/s\right]$$

(130)

$f_j=$ natural frequency:

$$f_j={\displaystyle \frac{{\omega }_j}{2\pi }}\left(cicluri/\sec \right)$$

(131)

$\left\{{\varphi }_j\right\}-$There are proper vectors which represent the vibration modes, respectively the assumed structure form when this vibrates with a frequency of $f_j$${\omega }_j$ and $\left\{{\varphi }_j\right\}$ which satisfies the equation:

$${\omega }_j^2\left[M\right]\left\{{\varphi }_j\right\}=\left[K\right]\left\{{\varphi }_j\right\}$$

(132)

For solving the motion equations system (120), there will be used the nodal superposition method. This method allows the user to detach the motion equations through a linear transformation as follows:

$$q_i\left(t\right)={\displaystyle \sum _{j=1}^n{\varphi }_{ij}{a}_j\left(t\right),i=\overline{1,n}},$$

(133)

where: ${\varphi }_{ij}-$ is the modal form coordinate which represents the ,,i” mass position in a ,,j” module; $a_j\left(t\right)-$is the generalized coordinate which represents the ,,j” module response during time; ${\varphi }_{ij}-$depends on the position, not on the time; $a_j\left(t\right)-$ depends on time, and not on the position.

Under a matrix form:

$$\left\{q\right\}=\left[\varphi \right]\cdot \left\{a\left(t\right)\right\},$$

(134)

where: $\left[\varphi \right]-$modal matrix, where the modal forms will be written in columns.

We accept that the modal forms vector is orthogonal in rapport with the mass matrix, which means:

$${\left[\varphi \right]}^T\left[M\right]\left[\varphi \right]=I–\mathrm{unit}\mathrm{matrix}$$

(135)

We reconsider the equation:

$$\left[M\right]\left\{\stackrel{··}{q}\right\}+\left[K\right]\left[q\right]=\left\{F\right\}$$

(136)

We differentiate twice the mathematical expression:

$$\left\{\stackrel{··}{q}\right\}=\left[\varphi \right]\left\{\ddot{a}\left(t\right)\right\}$$

(137)

We introduce (137) in (136):

$$\left[M\right]\cdot \left[\varphi \right]\cdot \left\{\ddot{a}\left(t\right)\right\}+\left[K\right]\cdot \left[\varphi \right]\cdot \left\{a\left(t\right)\right\}=\left\{F\left(t\right)\right\}$$

(138)

We multiply on the left side of the (138) relation with the ${\left[\varphi \right]}^T$ term and taking into account the (135), there will be obtained:

$$\left\{\ddot{a}\left(t\right)\right\}+{\left[\varphi \right]}^T\left[K\right]\left[\varphi \right]\left\{a\left(t\right)\right\}={\left[\varphi \right]}^T\left\{F\left(t\right)\right\}$$

(139)

The (132) relation, can be written as:

$${\omega }^2\left[M\right]\left[\varphi \right]=\left[K\right]\left[\varphi \right]$$

(140)

In the mathematical expression (140) we multiply on the left side with${\left\{\varphi \right\}}^T$, and taking into account the (135), we will obtain:

$$\left[{\omega }^2\right]={\left[\varphi \right]}^T\left[K\right]\left[\varphi \right]$$

(141)

where:

$$\left[{\omega }^2\right]=diag\left[{\omega }_1^2,{\omega }_2^2,\mathrm{...}{\omega }_n^2\right]$$

(142)

We introduce (141) in (139) and we will obtain:

$$\left\{\ddot{a}\left(t\right)\right\}+\left[{\omega }^2\right]\left\{a\left(t\right)\right\}={\left[\varphi \right]}^T\left\{F\left(t\right)\right\}$$

(143)

The (143) mathematical expression represents the uncoupled differential equation set which are independent ones. Thus, the modal frequencies will be:

$$\left[{\omega }^2\right]=diag\left[{\omega }_1^2,{\omega }_2^2,\mathrm{...}{\omega }_6^2\right]$$

(144)

The natural frequencies are:

$$f_i={\displaystyle \frac{{\omega }_i}{2\pi }},i=\overline{1,6}$$

(145)

where: f₁= 4,957487532Hz; f₂= 542,9548523Hz; f₃= 156,8132854Hz; f₄= 3230,572390Hz; f₅= 978,5360786Hz; f₆= 93,29474415Hz.

We reconsider the following equation:

$$\left(\left[K\right]-{\omega }^2\left[M\right]\right)\left\{\varphi \right\}=0$$

(146)

where: $\varphi -$ is the modal matrix form:

$$\varphi =\left[\begin{array}{cccccc}{\varphi }_{11}& {\varphi }_{12}& {\varphi }_{13}& {\varphi }_{14}& {\varphi }_{15}& {\varphi }_{16}\\ {\varphi }_{21}& {\varphi }_{22}& {\varphi }_{23}& {\varphi }_{24}& {\varphi }_{25}& {\varphi }_{26}\\ {\varphi }_{31}& {\varphi }_{32}& {\varphi }_{33}& {\varphi }_{34}& {\varphi }_{35}& {\varphi }_{36}\\ {\varphi }_{41}& {\varphi }_{42}& {\varphi }_{43}& {\varphi }_{44}& {\varphi }_{45}& {\varphi }_{46}\\ {\varphi }_{51}& {\varphi }_{52}& {\varphi }_{53}& {\varphi }_{54}& {\varphi }_{55}& {\varphi }_{56}\\ {\varphi }_{61}& {\varphi }_{62}& {\varphi }_{63}& {\varphi }_{64}& {\varphi }_{65}& {\varphi }_{66}\end{array}\right]$$

(147)

Based on numerical processing, there will be obtained the diagrams reported in Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16, where it can be identify six proper vibration modes of the analyzed mechanism.

The generalized coordinates vector is:

$$\left\{a\left(t\right)\right\}={\left[\begin{array}{cccccc}{a}_1\left(t\right)& a_2\left(t\right)& a_3\left(t\right)& a_4\left(t\right)& a_5\left(t\right)& a_6\left(t\right)\end{array}\right]}^T$$

(148)

This vector corresponds to the accelerations vector, which has the following form:

$$\left\{\stackrel{··}{a}\left(t\right)\right\}={\left[\begin{array}{cccccc}{\stackrel{··}{a}}_1\left(t\right)& {\stackrel{··}{a}}_2\left(t\right)& {\stackrel{··}{a}}_3\left(t\right)& {\stackrel{··}{a}}_4\left(t\right)& {\stackrel{··}{a}}_5\left(t\right)& {\stackrel{··}{a}}_6\left(t\right)\end{array}\right]}^T$$

(149)

The force vector will be:

$$\left\{F\left(t\right)\right\}={\left[\begin{array}{c}{F}_{31}^1\\ F_{41}^1+F_{21}^2\\ F_{51}^1+F_{21}^2\\ F_{61}^1+F_{31}^2\\ F_{41}^2\\ F_{61}^2\end{array}\right]}^T$$

(150)

$$\left\{\ddot{a}\left(t\right)\right\}+\left[{\omega }^2\right]\left\{a\left(t\right)\right\}-{\left[\varphi \right]}^T\left\{F\left(t\right)\right\}=0$$

(189)

By solving the differential equations system (189), the equation system solution (121), is:

$$\left\{q\right\}=\left[\varphi \right]\left\{a\left(t\right)\right\}$$

(152)

where $\left\{q\right\}$ - represents the nodal displacements. The elastic displacements can be obtained by solving the equation system:

$$\left\{u\right\}=\left[\text{}N\right]\left\{q\right\}$$

(153)

3.6. Results of Numerical Processing

For the numerical processing of the mathematical models, a program was developed in the Maple programming environment, which allows simultaneous inverse dynamic analysis—when the elements of the mechanism are rigid bodies—as well as elastodynamic and modal dynamic analysis when the elements are deformable bodies.

Virtual prototyping involved building a 3D model of the engine mechanism. An interface was created for importing this model into the virtual environment of the MSC Adams software.

In the first stage, a numerical simulation of the mechanism kinematic model was performed, with the crankshaft rotation speed (6000 [rpm]) as input data, along with the geometric characteristics, mass properties of the mechanism’s kinematic links, and the local reference systems attached to the considered links and kinematic joints. The mechanism motion was studied over a time interval of 0.04 seconds when the engine crankshaft performed four complete rotations.

In the second stage, the dynamic model was developed, in which the gas pressure applied on the piston head is known, along with the geometry, material characteristics, and inertial properties of the kinematic links.

The torque variation law acting on the crankshaft was determined, and consequently, the engine crankshaft angular speed variation.

The pressure developed by the gas on the piston head was experimentally determined, and its variation, as a function of the crankshaft rotation angle, was shown in Figure 5.

In order to determine the crankshaft torque (engine torque), a function was created that enables the interface with the ADAMS programming environment, as follows:

M_m = M₀ *(1-WX(MARKER_I, MARKER_J, MARKER_J)/ω₀)

(154)

Where: M_m is the engine torque; M₀ is the initial torque; ω₀ is the initial angular speed; WX(MARKER_I, MARKER_J is the X-axis component corresponding to the crankshaft angular speed, between marker I and marker J, measured in accordance with the coordinate system attached to marker J.

For the analyzed case study, the torque function is given by:

M_m = -170000*(1-WX(MARKER_1,MARKER_2,MARKER_2)/100)

(155)

Figure 17 represents snapshots from the virtual simulation carried out using the Adams software.

The three connecting rods were analyzed as deformable bodies, with the transformation from rigid to deformable bodies being based on the finite element method.

The connecting rod linear and angular deformations were determined in the Adams programming environment, based on the theory developed by Craig & Bampton [23]. The results were compared with those obtained from the processing of the mathematical models in the modal dynamic analysis, as it can be remarked in Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22.

By processing the relation (153), there will be obtained the time variation diagrams for the transversal and longitudinal displacements of the connecting rod as it can be seen in Figure 23, Figure 24, Figure 25, Figure 26, Figure 27 and Figure 28.

4. Experimental Analysis

To determine the piston stroke, an HBK WA-L type inductive displacement transducer was used, mounted in an opening made in the engine cylinder head. This transducer is attached to the piston and allows precise measurement of its stroke during the mechanism functionality.

For piston stroke measurement, it was used an inductive displacement transducer type HBK WA-L, mounted an fixed through a hole especially created in the engine cylinder head. This transducer has the mobile component fixed to the piston and allow the user to measure with a high fidelity the piston displacement during engine mechanism functionality.

The transducers’ characteristics are the following:

• Transducer type: Inductive, linear;
• High precision: Compatible with dynamic measurements;
• Variable measuring domain: Selected for covering the entire piston stroke;
• Rigid coupling: This assures a measurement without mechanical interferences.

This transducer is used to analyze the piston kinematics, providing experimental data for validating the theoretical models and the virtual simulations performed in MSC ADAMS.

By using a precision transducer, it is possible to compare the experimental results with the predictions obtained through numerical modeling, thus ensuring the validity of the data and identifying deviations from the impact theory applied in the simulations.

For the mechanism dynamic analysis, there were used piezoelectric acceleration transducers, mounted on the connecting rod and the cylinder head for measuring longitudinal and transversal accelerations.

The configuration of the experimental setup:

• The connecting rod accelerometers were attached by using magnets to the flat surface of the connecting rod, near the strain gauge areas, to ensure accurate measurement of longitudinal and transversal acceleration.
• The cylinder head accelerometers were mounted above the analyzed piston, capturing the vibrations and accelerations transmitted through the engine.

The accelerometer characteristics are as follows:

• Type: Piezoelectric;
• Measurement range: ±50 g;
• Sensitivity: 10 mV/g
• Mounting method: Magnets for the connecting rod, screws for the cylinder head;
• Frequency domain: 0.5 Hz – 10 kHz

The purposes of using the mentioned accelerometers are:

• To determine the vibrations induced during the connecting rod and piston motions.
• To compare the accelerations simulated in MSC.ADAMS with those experimentally measured.
• To validate the dynamic model through frequency spectrum analysis of the measured accelerations.

By using these accelerometers, valuable data were obtained on the system’s dynamic behavior, allowing for an accurate correlation between simulations and experimental measurements. The following parameters were recorded:

-AccBiela_L(m/s2) – connecting rod longitudinal acceleration;

-AccBiela_T(m/s2) – connecting rod transversal acceleration;

-AccBloc_V(m/s2) – cylinder head vertical acceleration;

-AccBloc_T(m/s2) – cylinder head transversal acceleration;

-CursaPiston(mm) – piston stroke;

-FortaBiele(N) – connecting rod longitudinal force.

For a single connecting rod setup, experimental tests were carried out successively at three different excitation frequencies. The electric motor’s operating frequency was adjusted using a frequency converter, with values around: 0.7 Hertz, 1 Hertz, and 1.3 Hertz.

For each excitation frequency, the operating time was approximately 30 seconds in order to allow for a frequency analysis over a reasonably wide bandwidth. It is important to note that the frequency resolution is inversely proportional to the analysis signal length:

Df=1/Dt

(194)

Considering the relation (194), in the frequency analyses performed, the frequency displayed on the screen represents the center of the frequency band, while the actual frequency lies within that band. For this reason, it was preferred that the excitation frequency of the piston be calculated from the time-domain representations as the inverse of the piston stroke period.

The experimental recordings were carried out using the LAN-XI data acquisition system, model 3053-B120, manufactured by Brüel & Kjær, along with the Pulse Labshop software from the same manufacturer.

Post-processing analyses were performed using the Pulse Reflex and TestPoint software packages, which complement each other, and results of these are shown in Figure 30 and Figure 31.

Figure 29. Experimental test rig for the dynamic analysis of the engine mechanism: Detailed view for the dynamic experimental analysis of the piston – connecting rod mechanism.

Figure 29. Experimental test rig for the dynamic analysis of the engine mechanism: Detailed view for the dynamic experimental analysis of the piston – connecting rod mechanism.

For the connecting rod, as it can be remarked in Figure 30, the accelerations predominantly exhibit a sinusoidal pattern, with two impact zones occurring on either side of the piston stroke’s peak moment, corresponding to the points where the direction of the force acting through the connecting rod changes on the two lateral faces of the piston.

From the time-domain analysis of the connecting rod vibrations, it is observed that the longitudinal and transversal accelerations predominantly exhibit a sinusoidal pattern, with two impact zones on either side of the maximum point of the piston stroke. These correspond to the moments when the direction of the force transmitted by the connecting rod changes to the two lateral faces of the piston, as it can be remarked Figure 31.

The longitudinal force in the connecting rod shows inflection points at the maximum and minimum points of the piston stroke. These inflection points should be clearly delimited in the case of a connecting rod without clearances in the joints.

5. Conclusions

The dynamic response of the multibody system consists of two segments: one due to the initial conditions, which damps out quickly, and one due to the disturbing forces.

By applying the finite element method, the continuous system was replaced with a discrete system with a finite number of degrees of freedom. The problem unknowns are no longer displacement functions u(x,y,z,t), but rather the nodal displacements q(t).

The system of partial differential equations (the Lame equations from elasticity theory) is transformed into a system of differential equations. The matrices involved in the general equation of motion were identified based on the idea of superimposing rigid body motion over deformable body motion. The rigid body motion was defined by a set of reference coordinates that define the location and orientation of the local reference system of each kinematic link.

By eliminating the rigid body motion, a system of linear differential equations is obtained, this being clearly evidentiate due to the assumptions made, which include: the small deformations hypothesis; the small displacements hypothesis; the linear material hypothesis (Hooke law).

Also, assuming that the damping force is proportional to the velocity, we consider a system of linear differential equations with constant coefficients ([M], [K], and [C]).

By combining the two motions, a complete dynamic analysis is proposed with multiple applications in multibody system dynamics. The elastodynamic analysis involves two stages: the analytical determination of natural frequencies and vibration modes, and the determination of kinematic parameters in the deformable body motion of the kinematic links.

The mathematical models obtained for the multibody systems elastodynamic analysis were applied to a crank-connecting rod mechanism used to drive a three-piston in-line engine. Each kinematic link was considered as a finite element “beam” type, with two nodes and three degrees of freedom per node, corresponding to the nodal displacements vector. The shape function matrices for axial and transversal displacements, as well as the generalized coordinate vector for the entire assembled system, were determined. For each finite element, the nodal force static and dynamic components, the stiffness matrix, and the damping matrix were determined relative to the local axis system and the global axis system, using coordinate transformation matrices when moving from one reference system to another. This was possible by superimposing rigid body motion with elastic solid motion, using Kane’s algorithm [22].

Using the modal superposition method, there were determined the linear, transversal, and angular nodal displacements for the connecting rod of the mechanism.

Thus, by using the finite element method, the variation of longitudinal and transversal elastic displacements in 3D format was obtained along with the entire length of the connecting rod, and in a 2D format at different points in the connecting rod plane.

In order to verify the results obtained from the processing of the mathematical models, a dynamic analysis was performed in the Adams programming environment based on multibody system dynamic theory.

The system’s motion is driven by the force developed by the gas on the piston head as a result of the gas pressure which was experimentally determined. The same force was considered in the mathematical models.

A comparative analysis was performed for the linear, transversal, and angular nodal displacements, where it was observed that the maximum and minimum values are close, including the shape of the variation curves.

From the 3D diagrams, it is observed that the maximum value of the longitudinal elastic displacements is 0.003 millimeters, while the transversal elastic displacements reach a value of 0.015 millimeters.

At the midpoint of the connecting rod length, the variation curves for the two types of longitudinal and transversal displacements register the same extreme values with similar variation shapes. The recorded differences mainly occur due to the fact that in the dynamic analysis with the Adams program, the engine torque exhibits a jump during startup over a period of a few milliseconds, then stabilizes around 150 [Nm].

The matrix approach to the mathematical models allowed the development of a flexible calculation program with rapid applications for the dynamic analysis of multibody systems with deformable elements.

The experimental analysis enabled the time and frequency monitoring of kinematic and dynamic parameters, mainly for each of the three connecting rods of the mechanism (forces, longitudinal, and transversal accelerations). The dynamic response monitoring of the three connecting rods was performed at lower speeds and without load on the piston, on a longitudinally sectioned engine.