**3. Numerical procedures for calculating periodic solutions of linear dynamic systems with time-periodic coefficients**

#### **3.1. Numerical procedure based on Runge-Kutta method**

Now we consider only the periodic vibration of a dynamic system which is governed by a set of linear differential equations with periodic coefficients. As already mentioned in the previous section, these differential equations can be expressed in the compact matrix form

$$
\dot{\mathbf{x}} = \mathbf{P}(t)\mathbf{x} + f(t) \tag{12}
$$

where *x* is the vector of state variables, matrix *P*(*t*) and vector *f* (*t*) are periodic in time with period *T*. The system of homogeneous differential equations corresponding to Eq. (12) is

$$
\dot{\mathbf{x}} = \mathbf{P}(t)\mathbf{x} \tag{13}
$$

As well known from the theory of differential equations, if Eq. (13) has only non-periodic solutions except the trivial solution, then Eq. (12) has an unique *T*-periodic solution. This pe‐ riodic solution can be obtained by choosing the appropriate initial condition for the vector of variables x and then implementing numerical integration of Eq. (12) within interval 0, *T* . An algorithm is developed to find the initial value for the periodic solution [18, 19]. Firstly, the *T*-periodic solution must satisfy the following condition

$$\mathbf{x}(\mathbf{0}) = \mathbf{x}(T) \tag{14}$$

The interval 0, *T* is now divided into *m* equal subintervals with the step-size *h* =*ti* −*ti*−<sup>1</sup> =*T* / *m*. At the discrete times *ti* and *ti*+1, *x<sup>i</sup>* = *x*(*ti* ) and *xi*+1 = *x*(*ti*+1) represent the states of the system, respectively. Using the fourth-order Runge-Kutta method, we get a nu‐ merical solution [5]

$$\mathbf{x}\_{i} = \mathbf{x}\_{i-1} + \frac{1}{6} \mathbf{f} \mathbf{k}\_{1}^{(i-1)} + 2\mathbf{k}\_{2}^{(i-1)} + 2\mathbf{k}\_{3}^{(i-1)} + \mathbf{k}\_{4}^{(i-1)}\tag{15}$$

where

Expansion of Eq. (9) yields a *n*-order algebraic equation

to Floquet multipliers can be expressed as follows.

equation, the solution *x* =0 of Eq. (6) is stable.

ution *x* =0 of Eq. (6) is also stable.

*<sup>λ</sup><sup>k</sup>* <sup>=</sup> <sup>1</sup>

(10). Floquet exponents are given by

304 Advances in Vibration Engineering and Structural Dynamics

ing larger than 1.

*<sup>ρ</sup> <sup>n</sup>* <sup>+</sup> *<sup>a</sup>*1*<sup>ρ</sup> <sup>n</sup>*−<sup>1</sup> <sup>+</sup> *<sup>a</sup>*2*<sup>ρ</sup> <sup>n</sup>*−<sup>2</sup> <sup>+</sup> .... <sup>+</sup> *an*−1*<sup>ρ</sup>* <sup>+</sup> *an* =0 (10)

*<sup>T</sup>* ln*ρ<sup>k</sup>* , (*<sup>k</sup>* =1, ..., *<sup>n</sup>*) (11)

where unknowns *ρ<sup>k</sup>* (*k* =1, ..., *n*), called Floquet multipliers, can be determined from Eq.

When the Floquet multipliers or Floquet exponents are known, the stability conditions of solutions of the system of linear differential equations with periodic coefficients can be easi‐ ly determined according to the Floquet theorem [17–20]. The concept of stability according

If |*ρ<sup>k</sup>* |1, the trivial solution *x* =0 of Eq. (6) will be asymptotically stable. Conversely, the solution *x* =0 of Eq. (6) becomes unstable if at least one Floquet multiplier has modulus be‐

If |*ρ<sup>k</sup>* | ≤1 and Floquet multipliers with modulus 1 are single roots of the characteristic

If |*ρ<sup>k</sup>* | ≤1 and Floquet multipliers with modulus 1 are multiple roots of the characteristic equation, and the algebraic multiplicity is equal to their geometric multiplicity, then the sol‐

**3. Numerical procedures for calculating periodic solutions of linear**

Now we consider only the periodic vibration of a dynamic system which is governed by a set of linear differential equations with periodic coefficients. As already mentioned in the previous section, these differential equations can be expressed in the compact matrix form

**dynamic systems with time-periodic coefficients**

**3.1. Numerical procedure based on Runge-Kutta method**

$$\begin{aligned} \mathbf{k}\_1^{(i-1)} &= \mathbb{I}\left[\mathbf{P}(t\_{i-1})\mathbf{x}\_{i-1} + f(t\_{i-1})\right] \\ \mathbf{k}\_2^{(i-1)} &= \mathbb{I}\left[\mathbf{P}(t\_{i-1} + \frac{h}{2})(\mathbf{x}\_{i-1} + \frac{1}{2}\mathbf{k}\_1^{(i-1)}) + f(t\_{i-1} + \frac{h}{2})\right] \\ \mathbf{k}\_3^{(i-1)} &= \mathbb{I}\left[\mathbf{P}(t\_{i-1} + \frac{h}{2})(\mathbf{x}\_{i-1} + \frac{1}{2}\mathbf{k}\_2^{(i-1)}) + f(t\_{i-1} + \frac{h}{2})\right] \\ \mathbf{k}\_4^{(i-1)} &= \mathbb{I}\left[\mathbf{P}(t\_i)(\mathbf{x}\_{i-1} + \mathbf{k}\_3^{(i-1)}) + f(t\_i)\right]. \end{aligned} \tag{16}$$

Substituting Eq. (16) into Eq. (15), we obtain

$$\mathbf{x}\_{i} = \mathbf{A}\_{i-1}\mathbf{x}\_{i-1} + \mathbf{b}\_{i-1} \tag{17}$$

where matrix *Ai*−1 is given by

$$\begin{aligned} A\_{i-1} &= \mathbf{I} + \frac{1}{6} \Big[ h\Big[ \mathbf{P}(t\_{i-1}) + 4\mathbf{P}(t\_{i-1} + \frac{h}{2}) + \mathbf{P}(t\_i) \Big] \\ &+ h^2 \Big[ \mathbf{P}(t\_{i-1} + \frac{h}{2}) \mathbf{P}(t\_{i-1}) + \mathbf{P}^2(t\_{i-1} + \frac{h}{2}) + \frac{1}{2} \mathbf{P}(t\_i) \mathbf{P}(t\_{i-1} + \frac{h}{2}) \Big] \\ &+ \frac{h^3}{2} \Big[ \mathbf{P}^2(t\_{i-1} + \frac{h}{2}) \mathbf{P}(t\_{i-1}) + \frac{1}{2} \mathbf{P}(t\_i) \mathbf{P}^2(t\_{i-1} + \frac{h}{2}) \Big] \\ &+ \frac{h^4}{4} \mathbf{P}(t\_i) \mathbf{P}^2(t\_{i-1} + \frac{h}{2}) \mathbf{P}(t\_{i-1}) \Big] (i = 1, \dots, m), \end{aligned} \tag{18}$$

**3.2. Numerical procedure based on Newmark integration method**

represent the solution of Eq. (2) at discrete times *ti*

tion must satisfy the following conditions

ing approximation formulas [6-7]

over the time interval [6, 7].

(*Mi*+1 <sup>+</sup> *<sup>γ</sup><sup>h</sup> <sup>C</sup>i*+1 <sup>+</sup> *<sup>β</sup><sup>h</sup>* <sup>2</sup>

ments at time *ti*+1

subintervals with the step-size *h* =*ti* −*ti*−<sup>1</sup> =*T* / *m*. We use notations *q<sup>i</sup>* =*q*(*ti*

*<sup>q</sup>i*+1 <sup>=</sup>*q<sup>i</sup>* <sup>+</sup> *<sup>h</sup> <sup>q</sup>***˙***<sup>i</sup>* <sup>+</sup> *<sup>h</sup>* 2( <sup>1</sup>

*<sup>q</sup>***˙***i*+1 <sup>=</sup>*q***˙***<sup>i</sup>* <sup>+</sup> (1−*γ*)*<sup>h</sup> <sup>q</sup>***¨**

*β* =1 / 6 leads to linear interpolation of accelerations in the time interval [ *ti*

From Eq. (2) we have the following iterative computational scheme at time *ti*+1

*<sup>i</sup>*+1 <sup>=</sup>*di*+1 <sup>−</sup>*Ci*+1 *<sup>q</sup>***˙***<sup>i</sup>* <sup>+</sup> (1−*γ*)*<sup>h</sup> <sup>q</sup>***¨**

( 1 <sup>2</sup> <sup>−</sup>*β*)*q***¨** *i* , *q***˙***i*+1

The use of Eqs. (23) and (24) leads to the prediction formulas for velocities and displace‐

*<sup>M</sup>i*+1*q***¨**

*<sup>K</sup>i*+1)*q***¨**

*qi*+1

Eq. (27) can be expressed in the matrix form as

\* <sup>=</sup>*q<sup>i</sup>* <sup>+</sup> *<sup>h</sup> <sup>q</sup>***˙***<sup>i</sup>* <sup>+</sup> *<sup>h</sup>* <sup>2</sup>

where *Mi*+1 =*M* (*ti*+1), *Ci*+1 =*C*(*ti*+1), *Ki*+1 = *K*(*ti*+1) and *di*+1 =*d*(*ti*+1).

In the next step, substitution of Eqs. (23) and (24) into Eq. (25) yields

The procedure presented below for finding the *T*-periodic solution of Eq. (2) is based on the Newmark direct integration method. Firstly, the interval 0, *T* is also divided into *m* equal

Based on the single-step integration method proposed by Newmark, we obtain the follow‐

<sup>2</sup> <sup>−</sup>*β*)*q***¨**

Constants *β*, *γ* are parameters associated with the quadrature scheme. Choosing *γ* =1 / 4 and

way, choosing *γ* =1 / 2, *β* =1 / 4 corresponds to considering the acceleration average value

*<sup>i</sup>* <sup>+</sup> *<sup>β</sup><sup>h</sup>* <sup>2</sup> *q***¨**

*<sup>i</sup>* <sup>+</sup> *<sup>γ</sup><sup>h</sup> <sup>q</sup>***¨**

) and *qi*+1 =*q*(*ti*+1) to

http://dx.doi.org/10.5772/51157

307

, *ti*+1 ]. In the same

*<sup>i</sup>* . (26)

and *ti*+1 respectively. The *T*-periodic solu‐

*<sup>i</sup>*+1, (23)

*<sup>i</sup>*+1, (24)

*q*(0)=*q*(*T* ), *q***˙**(0)=*q***˙**(*T* ), *q***¨**(0)=*q***¨**(*T* ). (22)

Parametric Vibration Analysis of Transmission Mechanisms Using Numerical Methods

*<sup>i</sup>*+1 + *Ci*+1*q***˙***i*+1 + *Ki*+1*qi*+1 =*di*+1, (25)

*<sup>i</sup>* <sup>−</sup> *<sup>K</sup>i*+1 *<sup>q</sup><sup>i</sup>* <sup>+</sup> *<sup>h</sup> <sup>q</sup>***˙***<sup>i</sup>* <sup>+</sup> *<sup>h</sup>* <sup>2</sup>

*i*

\* <sup>=</sup>*q***˙***<sup>i</sup>* <sup>+</sup> (1−*γ*)*<sup>h</sup> <sup>q</sup>***¨**

( 1 <sup>2</sup> <sup>−</sup>*β*)*q***¨**

. (27)

and vector *bi*−1 takes the form

$$\begin{split} & \mathbf{b}\_{i-1} = \frac{1}{6} \Big| h\left[ f(t\_{i-1}) + 4f\left( t\_{i-1} + \frac{h}{2} \right) + f\left( t\_i \right) \right] \\ & + h^2 \Big[ \mathbf{P}\left( t\_{i-1} + \frac{h}{2} \right) f\left( t\_{i-1} \right) + \mathbf{P}\left( t\_{i-1} + \frac{h}{2} \right) f\left( t\_{i-1} + \frac{h}{2} \right) + \frac{1}{2} \mathbf{P}\left( t\_i \right) f\left( t\_{i-1} + \frac{h}{2} \right) \Big] \\ & + \frac{h^3}{2} \Big[ \mathbf{P}^2(t\_{i-1} + \frac{h}{2}) f\left( t\_{i-1} \right) + \mathbf{P}(t\_i) \mathbf{P}\left( t\_{i-1} + \frac{h}{2} \right) f\left( t\_{i-1} + \frac{h}{2} \right) \Big] \\ & + \frac{h^4}{4} \mathbf{P}(t\_i) \mathbf{P}^2(t\_{i-1} + \frac{h}{2}) f\left( t\_{i-1} \right) \Big]. \end{split} \tag{19}$$

Expansion of Eq. (17) for *i* =1 *to m* yields

$$\begin{aligned} \mathbf{x}\_1 &= \mathbf{A}\_0 \mathbf{x}\_0 + \mathbf{c}\_1\\ \mathbf{x}\_2 &= \mathbf{A}\_1 \mathbf{A}\_0 \mathbf{x}\_0 + \mathbf{c}\_2\\ \mathbf{x}\_3 &= \dots + \mathbf{x}\_{\text{new}} + \dots + \dots \\ \mathbf{x}\_m &= \left(\prod\_{i=m-1}^0 \mathbf{A}\_i\right) \mathbf{x}\_0 + \mathbf{c}\_m \end{aligned} \tag{20}$$

where *c*<sup>0</sup> =**0**, *c*<sup>1</sup> = *A*0*c*<sup>0</sup> + *b*0, *c*<sup>2</sup> = *A*1*c*<sup>1</sup> + *b*1 ,..., *c<sup>m</sup>* = *Am*−1*cm*−<sup>1</sup> + *bm*−1. Using the boundary condi‐ tion according to Eq. (14), the last equation of Eq. (20) yields a set of the linear algebra‐ ic equations

$$\left(I - \prod\_{i=w-1}^{0} A\_i\right) \mathbf{x}\_0 = \mathbf{c}\_m. \tag{21}$$

The solution of Eq. (21) gives us the initial value for the periodic solution of Eq. (12). Finally, the periodic solution of Eq. (12) with the corresponding initial value can be calculated using the computational scheme according to Eq. (15).

### **3.2. Numerical procedure based on Newmark integration method**

*Ai*−<sup>1</sup> = *I* +

+ *h* 3 <sup>2</sup> *<sup>P</sup>* <sup>2</sup>

+ *h* 4 <sup>4</sup> *<sup>P</sup>*(*ti*

and vector *bi*−1 takes the form

<sup>+</sup> *<sup>h</sup>* <sup>2</sup> *<sup>P</sup>*(*ti*−<sup>1</sup> <sup>+</sup>

*<sup>b</sup>i*−<sup>1</sup> <sup>=</sup> <sup>1</sup>

+ *h* 3 <sup>2</sup> *<sup>P</sup>* <sup>2</sup>

+ *h* 4 <sup>4</sup> *<sup>P</sup>*(*ti*

ic equations

<sup>+</sup> *<sup>h</sup>* <sup>2</sup> *<sup>P</sup>*(*ti*−<sup>1</sup> <sup>+</sup>

306 Advances in Vibration Engineering and Structural Dynamics

1

*h*

(*ti*−<sup>1</sup> + *h*

)*P* <sup>2</sup> (*ti*−<sup>1</sup> + *h*

<sup>6</sup> {*<sup>h</sup> <sup>f</sup>* (*ti*−1) <sup>+</sup> <sup>4</sup> *<sup>f</sup>* (*ti*−<sup>1</sup> <sup>+</sup>

*h*

(*ti*−<sup>1</sup> + *h*

)*P* <sup>2</sup> (*ti*−<sup>1</sup> + *h*

Expansion of Eq. (17) for *i* =1 *to m* yields

the computational scheme according to Eq. (15).

<sup>6</sup> {*<sup>h</sup> <sup>P</sup>*(*ti*−1) <sup>+</sup> <sup>4</sup>*P*(*ti*−<sup>1</sup> <sup>+</sup>

<sup>2</sup> )*P*(*ti*−1) <sup>+</sup> *<sup>P</sup>* <sup>2</sup>

<sup>2</sup> )*P*(*ti*−1) <sup>+</sup>

*h* <sup>2</sup> ) <sup>+</sup> *<sup>f</sup>* (*ti*

<sup>2</sup> ) *<sup>f</sup>* (*ti*−1) <sup>+</sup> *<sup>P</sup>*(*ti*−<sup>1</sup> <sup>+</sup>

<sup>2</sup> ) *<sup>f</sup>* (*ti*−1) <sup>+</sup> *<sup>P</sup>*(*ti*

<sup>2</sup> ) *<sup>f</sup>* (*ti*−1)}.

*h* <sup>2</sup> ) <sup>+</sup> *<sup>P</sup>*(*ti*

(*ti*−<sup>1</sup> + *h* <sup>2</sup> ) <sup>+</sup> 1 <sup>2</sup> *<sup>P</sup>*(*ti*

<sup>2</sup> )*P*(*ti*−1)}(*<sup>i</sup>* =1, ..., *<sup>m</sup>*),

)

<sup>2</sup> ) *<sup>f</sup>* (*ti*−<sup>1</sup> <sup>+</sup>

*Ai*)*x*<sup>0</sup> + *c<sup>m</sup>*

where *c*<sup>0</sup> =**0**, *c*<sup>1</sup> = *A*0*c*<sup>0</sup> + *b*0, *c*<sup>2</sup> = *A*1*c*<sup>1</sup> + *b*1 ,..., *c<sup>m</sup>* = *Am*−1*cm*−<sup>1</sup> + *bm*−1. Using the boundary condi‐ tion according to Eq. (14), the last equation of Eq. (20) yields a set of the linear algebra‐

The solution of Eq. (21) gives us the initial value for the periodic solution of Eq. (12). Finally, the periodic solution of Eq. (12) with the corresponding initial value can be calculated using

*h*

*h* <sup>2</sup> ) <sup>+</sup> 1 <sup>2</sup> *<sup>P</sup>*(*ti*

<sup>2</sup> ) *<sup>f</sup>* (*ti*−<sup>1</sup> <sup>+</sup>

*h* 2 )

*h*

)*P*(*ti*−<sup>1</sup> +

*x*<sup>1</sup> = *A*0*x*<sup>0</sup> + *c*<sup>1</sup> *x*<sup>2</sup> = *A*1*A*0*x*<sup>0</sup> + *c*<sup>2</sup> ................................

*<sup>x</sup><sup>m</sup>* =( ∏ *i*=*m*−1 0

(*I* − ∏ *i*=*m*−1 0

1 <sup>2</sup> *<sup>P</sup>*(*ti* )*P* <sup>2</sup> (*ti*−<sup>1</sup> + *h* 2 )

)

)*P*(*ti*−<sup>1</sup> +

*h* 2 )

) *f* (*ti*−<sup>1</sup> +

*Ai*)*x*<sup>0</sup> =*cm*. (21)

*h* 2 ) (18)

(19)

(20)

The procedure presented below for finding the *T*-periodic solution of Eq. (2) is based on the Newmark direct integration method. Firstly, the interval 0, *T* is also divided into *m* equal subintervals with the step-size *h* =*ti* −*ti*−<sup>1</sup> =*T* / *m*. We use notations *q<sup>i</sup>* =*q*(*ti* ) and *qi*+1 =*q*(*ti*+1) to represent the solution of Eq. (2) at discrete times *ti* and *ti*+1 respectively. The *T*-periodic solu‐ tion must satisfy the following conditions

$$
\dot{q}(0) = q(T), \ \dot{q}(0) = \dot{q}(T), \ \ddot{q}(0) = \ddot{q}(T). \tag{22}
$$

Based on the single-step integration method proposed by Newmark, we obtain the follow‐ ing approximation formulas [6-7]

$$q\_{i+1} = q\_i + h \, \dot{q}\_i + h \, ^2 \left(\frac{1}{2} - \beta\right) \ddot{q}\_i + \beta h \, ^2 \dddot{q}\_{i+1} \tag{23}$$

$$
\dot{q}\_{i+1} = \dot{q}\_i + (1 - \gamma^\prime) h \stackrel{\cdots}{q}\_i + \gamma h \stackrel{\cdots}{q}\_{i+1^\prime} \tag{24}
$$

Constants *β*, *γ* are parameters associated with the quadrature scheme. Choosing *γ* =1 / 4 and *β* =1 / 6 leads to linear interpolation of accelerations in the time interval [ *ti* , *ti*+1 ]. In the same way, choosing *γ* =1 / 2, *β* =1 / 4 corresponds to considering the acceleration average value over the time interval [6, 7].

From Eq. (2) we have the following iterative computational scheme at time *ti*+1

$$\mathbf{M}\_{i+1}\ddot{\mathbf{q}}\_{i+1} + \mathbf{C}\_{i+1}\dot{\mathbf{q}}\_{i+1} + \mathbf{K}\_{i+1}\mathbf{q}\_{i+1} = \mathbf{d}\_{i+1} \tag{25}$$

where *Mi*+1 =*M* (*ti*+1), *Ci*+1 =*C*(*ti*+1), *Ki*+1 = *K*(*ti*+1) and *di*+1 =*d*(*ti*+1).

In the next step, substitution of Eqs. (23) and (24) into Eq. (25) yields

$$\left[\mathbf{M}\_{i+1} + \gamma h \, \mathbf{C}\_{i+1} + \beta h \, \, ^2 \mathbf{K}\_{i+1}\right] \ddot{\boldsymbol{q}}\_{i+1} = \mathbf{d}\_{i+1} - \mathbf{C}\_{i+1} \left[\dot{\boldsymbol{q}}\_{i} + (1 - \gamma) h \, \ddot{\boldsymbol{q}}\_{i}\right] - \mathbf{K}\_{i+1} \left[\boldsymbol{q}\_{i} + h \, \dot{\boldsymbol{q}}\_{i} + h \, ^2 \left(\frac{1}{2} - \beta\right) \ddot{\boldsymbol{q}}\_{i}\right].\tag{26}$$

The use of Eqs. (23) and (24) leads to the prediction formulas for velocities and displace‐ ments at time *ti*+1

$$
\dot{q}\_{i+1}^\* = q\_i + h \, \dot{q}\_i + h \, ^2 \Big( \frac{1}{2} - \beta \Big) \stackrel{\cdots}{\dot{q}}\_{i\prime} \cdot \dot{q}\_{i+1}^\* = \dot{q}\_i + (1 - \gamma) h \, \stackrel{\cdots}{\dot{q}}\_{i\prime} \tag{27}
$$

Eq. (27) can be expressed in the matrix form as

$$
\begin{bmatrix}
\mathbf{q}\_{\iota+1}^{\*} \\
\dot{\mathbf{q}}\_{\iota+1}^{\*}
\end{bmatrix} = \mathbf{D} \begin{bmatrix}
\mathbf{q}\_{\iota} \\
\dot{\mathbf{q}}\_{\iota} \\
\ddot{\mathbf{q}}\_{\iota}
\end{bmatrix} \tag{28}
$$

The combination of Eqs. (28), (33) and (34) yields a new computational scheme for determin‐

( ) ( )

1 1 1 <sup>1</sup> 1 1

+ + - <sup>+</sup> + + é ù é ù é ù é ù ê ú ê ú ê ú = + ê ú ê ú ê ú ê ú ê ú ë û - ê ú ê ú ê ú ë û ë û ê ú ë û **<sup>q</sup> q 0 <sup>D</sup> q T qT 0 S HD q q S d**

In this equation, the iterative computation is eliminated by introducing the direct solution

( ) ( )

1 1 1

+ + -

1 1 1

é ù é ù é ù ê ú ê ú = = ê ú <sup>=</sup> ê ú ê ú ê ú ë û - ê ú ê ú ë û ê ú ë û **q 0 <sup>D</sup> xq A T bT 0 S HD q S d**

*i i <sup>i</sup> i i* + + -

> *x*<sup>1</sup> = *A*1*x*<sup>0</sup> + *c*<sup>1</sup> *x*<sup>2</sup> = *A*2*A*1*x*<sup>0</sup> + *c*<sup>2</sup> ................................

*<sup>x</sup><sup>m</sup>* =(∏ *i*=*m* 1

> (*I* −∏ *i*=*m* 1

*Ai*)*x*<sup>0</sup> + *c<sup>m</sup>*

Using the condition of periodicity according to Eq. (22), the last equation of Eq. (39) yields a

The solution of Eq. (40) gives us the initial value for the periodic solution of Eq. (2). Finally, the periodic solution of Eq. (2) with the obtained initial value can be calculated without diffi‐

, ,,

.

Parametric Vibration Analysis of Transmission Mechanisms Using Numerical Methods

1 1

+ +

*x<sup>i</sup>* = *Aixi*−<sup>1</sup> + *b<sup>i</sup>* (*i* =1, 2, ..., *m*). (38)

*Ai*)*x*<sup>0</sup> =*cm*. (40)

& (37)

(36)

309

http://dx.doi.org/10.5772/51157

(39)

ing the solution of Eq. (2) at the time *ti*+1 in the form

1

+

*i*

Eq. (36) can then be rewritten in the following form

By setting

1 1

+ -

for each time step. Note that *T* and *D* are matrices of constants.

*i i i i i i i i i i*

& & && &&

*ii i i*

Expansion of Eq. (38) for *i* =1 *to m* yields the same form as Eq. (20)

where *c*<sup>0</sup> =**0**, *c*<sup>1</sup> = *A*1*c*<sup>0</sup> + *b*1, *c*<sup>2</sup> = *A*2*c*<sup>1</sup> + *b*2 ,..., *c<sup>m</sup>* = *Amcm*−<sup>1</sup> + *bm*.

culties using the computational scheme in Eq. (36).

set of the linear algebraic equations

with

$$D = \begin{bmatrix} I & h \ I & h^2 (0.5 - \beta) I \\ \mathbf{0} & I & (1 - \gamma) h \ I \end{bmatrix} \tag{29}$$

where **0** represents the *n* ×*n* matrix of zeros. Eq. (26) can then be rewritten in the matrix form as

$$\ddot{\mathbf{q}}\_{\ell+1} = \left(\mathbf{S}\_{\ell+1}\right)^{-1} \mathbf{d}\_{\ell+1} - \left(\mathbf{S}\_{\ell+1}\right)^{-1} \mathbf{H}\_{\ell+1} \begin{bmatrix} \mathbf{q}\_{\ell+1}^{\*} \\ \mathbf{q}\_{\ell+1}^{\*} \end{bmatrix},\tag{30}$$

where matrices *Si*+1 and *Hi*+1 are defined by

$$\mathbf{S}\_{i+1} = \mathbf{M}\_{i+1} + \gamma h \, \mathbf{C}\_{i+1} + h \, ^2 \beta \mathbf{K}\_{i+1'} \tag{31}$$

$$\mathbf{H}\_{i+1} = \begin{bmatrix} \mathbf{K}\_{i+1} & \mathbf{C}\_{i+1} \end{bmatrix} . \tag{32}$$

By substituting relationships (28) into (30) we find

$$\boldsymbol{\ddot{\mathbf{q}}}\_{i+1} = \left(\mathbf{S}\_{i+1}\right)^{-1} \mathbf{d}\_{i+1} - \left(\mathbf{S}\_{i+1}\right)^{-1} \mathbf{H}\_{i+1} \mathbf{D} \begin{bmatrix} \mathbf{q}\_{i} \\ \mathbf{q}\_{i} \\ \ddot{\mathbf{q}}\_{i} \end{bmatrix} \tag{33}$$

From Eqs. (23), (24) and (27) we get the following matrix relationship

$$
\begin{bmatrix}
\mathbf{q}\_{\ell+1} \\
\dot{\mathbf{q}}\_{\ell+1} \\
\dot{\mathbf{q}}\_{\ell+1}
\end{bmatrix} = \mathbf{T} \begin{bmatrix}
\dot{\mathbf{q}}\_{\ell+1}^\* \\
\dot{\mathbf{q}}\_{\ell+1}^\*
\end{bmatrix},
\tag{34}
$$

where matrix *T* is expressed in the block matrix form as

$$\mathbf{T} = \begin{bmatrix} \mathbf{I} & \mathbf{0} & \mathbf{I}\beta h^2 \\ \mathbf{0} & \mathbf{I} & \mathbf{I}\gamma h \\ \mathbf{0} & \mathbf{0} & \mathbf{I} \end{bmatrix} \tag{35}$$

The combination of Eqs. (28), (33) and (34) yields a new computational scheme for determin‐ ing the solution of Eq. (2) at the time *ti*+1 in the form

$$
\begin{bmatrix}
\mathbf{q}\_{i+1} \\
\dot{\mathbf{q}}\_{i+1} \\
\ddot{\mathbf{q}}\_{i+1}
\end{bmatrix} = \mathbf{T} \begin{bmatrix}
\mathbf{D} \\
\mathbf{q}\_{i} \\
\dot{\mathbf{q}}\_{i} \\
\ddot{\mathbf{q}}\_{i}
\end{bmatrix} + \mathbf{T} \begin{bmatrix}
\mathbf{0} \\
\mathbf{0} \\
\left(\mathbf{S}\_{i+1}\right)^{-1} \mathbf{d}\_{i+1}
\end{bmatrix}.
\tag{36}
$$

In this equation, the iterative computation is eliminated by introducing the direct solution for each time step. Note that *T* and *D* are matrices of constants.

By setting

(28)

(33)

(34)

\* 1 \* 1

*i*

+ + é ù é ù ê ú ê ú <sup>=</sup> ê ú ë û ê ú

*i*

*<sup>D</sup>* <sup>=</sup> *<sup>I</sup> <sup>h</sup> <sup>I</sup> <sup>h</sup>* 2(0.5−*β*)*<sup>I</sup>*

with

as

*i*

**q <sup>q</sup> D q q**

*i*

*i*

**<sup>0</sup>** *<sup>I</sup>* (1−*γ*)*<sup>h</sup> <sup>I</sup>* (29)

ë û

**q** & & &&

where **0** represents the *n* ×*n* matrix of zeros. Eq. (26) can then be rewritten in the matrix form

( ) ( ) \* 1 1 <sup>1</sup> 1 11 1 1 \*

é ù = - ê ú


*i ii i i*

*<sup>S</sup>i*+1 <sup>=</sup>*Mi*+1 <sup>+</sup> *<sup>γ</sup><sup>h</sup> <sup>C</sup>i*+1 <sup>+</sup> *<sup>h</sup>* <sup>2</sup>

( ) ( ) 1 1 1 11 1 1

From Eqs. (23), (24) and (27) we get the following matrix relationship

where matrix *T* is expressed in the block matrix form as

*i ii i ii*

**q S d S H Dq**

&& &&&

*i i i i i i*

**q q q Tq q q** & & && &&

é ù é ù ê ú ê ú <sup>=</sup> ê ú ê ú ê ú ê ú ë û ë û

+ + + + + +

**I 0I T 0I I 00I**


ê ú = - ê ú

\* 1 1 \* 1 1 1 1

,

2 *h h* é ù ê ú <sup>=</sup> ê ú ê ú ë û

b

g

where matrices *Si*+1 and *Hi*+1 are defined by

308 Advances in Vibration Engineering and Structural Dynamics

By substituting relationships (28) into (30) we find

+ ++ + +

**<sup>q</sup> q Sd SH**

1 , *<sup>i</sup>*

*i*

é ù

**q**

*i*

ê ú ë û

**q**

**<sup>q</sup>** && &(30)

*<sup>H</sup>i*+1 <sup>=</sup> *<sup>K</sup>i*+1 *<sup>C</sup>i*+1 . (32)

*βKi*+1, (31)

(35)

+

*i*

ë û

$$\mathbf{x}\_{i} = \begin{bmatrix} \mathbf{q}\_{i} \\ \dot{\mathbf{q}}\_{i} \\ \dot{\mathbf{q}}\_{i} \end{bmatrix}, \qquad \mathbf{A}\_{i+1} = \mathbf{T} \begin{bmatrix} \mathbf{D} \\ -\left(\mathbf{S}\_{i+1}\right)^{-1} \mathbf{H}\_{i+1} \mathbf{D} \end{bmatrix}, \qquad \mathbf{b}\_{i+1} = \mathbf{T} \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \left(\mathbf{S}\_{i+1}\right)^{-1} \mathbf{d}\_{i+1} \end{bmatrix}, \tag{37}$$

Eq. (36) can then be rewritten in the following form

$$\mathbf{x}\_{i} = \mathbf{A}\_{i}\mathbf{x}\_{i-1} + \mathbf{b}\_{i} \ \text{ (} i = 1, \ \text{ 2, ..., } \ m \text{)}. \tag{38}$$

Expansion of Eq. (38) for *i* =1 *to m* yields the same form as Eq. (20)

$$\begin{aligned} \mathbf{x}\_1 &= A\_1 \mathbf{x}\_0 + \mathbf{c}\_1\\ \mathbf{x}\_2 &= A\_2 A\_1 \mathbf{x}\_0 + \mathbf{c}\_2\\ \vdots &= \vdots + \vdots + \vdots + \vdots + \vdots\\ \mathbf{x}\_m &= \left(\prod\_{i=m}^1 A\_i\right) \mathbf{x}\_0 + \mathbf{c}\_m \end{aligned} \tag{39}$$

where *c*<sup>0</sup> =**0**, *c*<sup>1</sup> = *A*1*c*<sup>0</sup> + *b*1, *c*<sup>2</sup> = *A*2*c*<sup>1</sup> + *b*2 ,..., *c<sup>m</sup>* = *Amcm*−<sup>1</sup> + *bm*.

Using the condition of periodicity according to Eq. (22), the last equation of Eq. (39) yields a set of the linear algebraic equations

$$\left(I - \prod\_{i=m}^{1} A\_i\right) \mathbf{x}\_0 = \mathbf{c}\_m. \tag{40}$$

The solution of Eq. (40) gives us the initial value for the periodic solution of Eq. (2). Finally, the periodic solution of Eq. (2) with the obtained initial value can be calculated without diffi‐ culties using the computational scheme in Eq. (36).

Based on the proposed numerical procedures in this section, a computer program with MATLAB to calculate periodic vibrations of transmission mechanisms has been developed and tested by the following application examples.

the combined stiffness of the return spring and the support of the output link. The rotating components are modeled by two rotating disks with moments of inertia *I*0 and *I*1. Let us in‐ troduce into our dynamic model the nonlinear transmission function *U* (*φ*1) of the cam mechanism as a function of the rotating angle *φ*1 of the cam shaft, the driving torque from

Parametric Vibration Analysis of Transmission Mechanisms Using Numerical Methods

http://dx.doi.org/10.5772/51157

311

The kinetic energy, the potential energy and the dissipative function of the considered sys‐

<sup>2</sup> *k*2(*y*<sup>2</sup> − *y*1)<sup>2</sup> +

<sup>2</sup> *c*2(*y*˙ <sup>2</sup> − *y*˙ 1)<sup>2</sup> +

1

1

∑ *<sup>δ</sup>A*=*<sup>M</sup>* (*t*)*δφ*<sup>0</sup> <sup>−</sup> *<sup>F</sup>* (*t*)*<sup>δ</sup> <sup>y</sup>*<sup>3</sup> (44)

2 (41)

<sup>2</sup> *k*3(*y*<sup>3</sup> − *y*2)<sup>2</sup> (42)

<sup>2</sup> *c*3(*y*˙ <sup>3</sup> − *y*˙ 2)<sup>2</sup> (43)

the motor *M*(*t*) and the external load *F*(*t*) applied on the system.

**Figure 2.** Dynamic model of the cam mechanism.

tem can be expressed in the following form

*<sup>T</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> *I*0*φ*˙ <sup>0</sup> 2 + 1 <sup>2</sup> *I*1*φ*˙ <sup>1</sup> 2 + 1 <sup>2</sup> *m*2*y*˙ <sup>2</sup> 2 + 1 <sup>2</sup> *m*3*y*˙ <sup>3</sup>

<sup>2</sup> *k*1(*φ*<sup>1</sup> −*φ*0)<sup>2</sup> +

<sup>2</sup> *c*1(*φ*˙ <sup>1</sup> −*φ*˙ 0)<sup>2</sup> +

The virtual work done by all non-conservative forces is

1

1

*<sup>Π</sup>* <sup>=</sup> <sup>1</sup>

*<sup>Φ</sup>* <sup>=</sup> <sup>1</sup>
