Comput Mech (2008) 41:269–277 DOI 10.1007/s00466-007-0214-6

ORIGINAL PAPER

Study on sub-cycling algorithm for flexible multi-body system: stability analysis and numerical examples J. C. Miao · P. Zhu · G. L. Shi · G. L. Chen

Received: 11 September 2006 / Accepted: 12 April 2007 / Published online: 16 August 2007 © Springer-Verlag 2007

Abstract Numerical stability is an important issue for any integral procedure. Since sub-cycling algorithm has been presented by Belytschko et al. (Comput Methods Appl Mech Eng 17/18: 259–275, 1979), various kinds of these integral procedures were developed in later 20 years and their stability were widely studied. However, on how to apply the sub-cycling to flexible multi-body dynamics (FMD) is still a lack of investigation up to now. A particular sub-cycling algorithm for the FMD based on the central difference method was introduced in detail in part I (Miao et al. in Comp Mech doi: 10.1007/s00466-007-0183-9) of this paper. Adopting an integral approximation operator method, stability of the presented algorithm is transformed to a generalized eigenvalue problem in the paper and is discussed by solving the problem later. Numerical examples are performed to verify the availability and efficiency of the algorithm further. Keywords Flexible multi-body system (FMS) · Flexible multi-body dynamics (FMD) · Sub-cycling algorithm · Amplification matrix · General eigenvalue problem

1 Introduction In part I [2] of this paper, a central-difference-based subcycling for flexible multi-body dynamics (FMD) problems

was developed. We introduced its basic ideas and established its formulae in the previous paper. Stability is one of the most important properties for any applied algorithm. It means that as a perturbation was introduced in during an integral process or at the start of the process, the error can be depressed or limited in a boundary at the following steps. This is a pledge for convergence of a practical integral procedure. In addition, Stability analysis is desirable for seeking the appropriate step size of an integral procedure. There are two ways widely used in analyzing the stability of a sub-cycling algorithm. One is called the energy balance method, which adopted in references [3–8]; the other one is called the integral approximation operator method [9–14]. In the integral approximation operator method, a transform matrix, which reflects the relationship of the value of the variables at n + 1 step, Xn+1 , and the value of the variables at n step, Xn , should be established in advance. It can be written as the following form [15]. ˆ t + Lrt+v ˆ t+
t = AX X

(1)

Here, the matrix, A, is called the amplification matrix of the integral. It is a function of
t, the general mass matrix and the general stiffness matrix of an FMD formula. Then the stability problem is equivalent to an eigenvalue problem. AU = λU

(2)

Or a generalized eigenvalue problem. J. C. Miao · P. Zhu · G. L. Shi · G. L. Chen School of Mechanical Engineering, Shanghai JiaoTong University, Shanghai, China J. C. Miao (B) Department of Mechanical Engineering, Shazhou Institute of Technology, Zhangjiagang, China e-mail: [email protected]

AU = λBU

(3)

Here, λ is the eigenvalues of Eq. (3) and is constituted by λi . The relationship of (
t/T ) and λi can be gained after Eq. (3) was solved. Stability of the integral process can be assessed in terms of that whether the maximum modulus of the eigenvalues, or called the spectrum radius, is less than 1.

123

270

Comput Mech (2008) 41:269–277

If the maximum modulus of the eigenvalues is less than 1, namely max|λi | < 1, the algorithm is proved to be a stable process in an exponential sense. Otherwise, the perturbation error will be linear increase as max|λi | = 1, or exponential increase as max|λi | > 1. We will apply this method to analyze the stability of the sub-cycling algorithm, which was presented in part I [2] of the paper, in the following sections.

q L represent the rapidly changing variables and the slowly changing variables respectively. In a structural dynamic analysis, if the damping is small, popularly does not change the stability of an integral procedure. In the same sense, we will ignore the damping in Eq. (6). Then it can be simplified as following.  
 
ˆt ˆt ˆt 0 M M K q¨ L qL LL LS L + =0 (7) t t t ˆ q¨ S qS ˆ ˆ 0 KS M M SL

2 Amplification matrices and stability analysis To obtain the formula of the amplification matrix, we rewrite the condensed equation of the FMD as follows. ˆ i (q˙ i , qi , q˙ d , qd , t) ˆ i (qi , qd , t) q¨ i = Q M

(4)

ˆ i is the general ˆ i is the general mass matrix and Q Herein, M force. There are three components included in the general force, which are the external force, the general force related with the quadratic term of the velocity, and the general internal force respectively. The general internal force is introduced by means of the elastic deformations of the flexible bodies. The general force related with the quadratic term of the velocity is corresponding to the centrifugal force and the Coriolis force. We can regard these two components as similar to the damping in a structural dynamic analysis. And if the deformation of a flexible body is transformed to the general mode coordinates, the general internal forces are uncoupling each other. In other words, the general stiffness matrix of an FMD formula will be a diagonal matrix. Omitting the general external forces contained in the general force, we obtain an equivalent expanded form of Eq. (4): ˆ i (q˙ i , qi , q˙ d , qd , t) + K ˆ i qi = 0 ˆ i (qi , qd , t)q¨ i + C M

(5)

ˆ i (q˙ i , qi , q˙ d , qd , t) − ˆ i (q˙ i , qi , q˙ d , qd , t) = Q There into, C ˆ Ki qi . For stability analysis of an FMD sub-cycling procedure, the mass matrix and the stiffness matrix of a linear elastic model need be partitioned into sub-matrices corresponding to the larger step size region and the smaller step size region. Using these sub-matrices, Eq. (5) can be written in a decomposed matrices form below. 

t 

t ˆ ˆt ˆt ˆ q¨ L C q˙ L M C M LL LS LL LS + ˆt ˆt ˆt ˆt q¨ S q˙ S M M C C SL SS SL SS
t 
 ˆ K 0 qL L + =0 (6) ˆt qS 0 K S where, the subscripts L and S are corresponding to the variables changed in a major cycle and the variables changed in a sub-cycle in the sub-cycling process, respectively. q S and

123

SS

In terms of the central difference formulae in reference [14], the simplified common update Eq. (7) can now be written below.  ˆ t  t−
t M t+
t L t LL L q − 2q + q L L
t L2 L  ˆt  M t+
t S t S ˆ t qt = 0 +K − 2q + q + L2S qt−
t S L L S
t S S  ˆ t  t−
t M t+
t L SL t L q − 2q + q L L
t L2 L  ˆt  M t+
t S t S ˆ t qt = 0 +K − 2q + q + SS2 qt−
t S S S S
t S S

(8)

(9)

The matrix formation of formulae (8)–(9) is as following. ⎤ ⎡ t ˆt ˆ ⎡ t+
t L ⎤ M M LS LL qL ⎢
t L2
t S2 ⎥ ⎥⎣ ⎢ ⎦ ⎣ Mˆ t ˆt ⎦ M t+
t S SL SS qS
t L2
t 2 ⎡ St ⎤ ˆt ˆ M M 1 ˆt
 LS LL − K 2 L ⎢
t L2 ⎥ qtL
t S2 ⎥ = −2 ⎢ ⎣ Mˆ t ⎦ qt ˆt M 1 ˆt SL SS S − 2 KS
t L2
t S2 ⎤ ⎡ t ˆt ˆ ⎡ t−
t ⎤ M M LS LL L ⎢
t L2
t S2 ⎥ q L ⎥ ⎢ ⎣ ⎦ +⎣ t (10) ˆ ˆt ⎦ M M t−
t S SL SS q S 2 2
t L


t S

Also, the simplified sub-step updating formula can be written as following, with 2
j
m, m is the number of sub-steps within a major cycle.  ˆ t+( j−1)
t S  t+( j−2)
t M t+( j−1)
t S t+ j·
t S SS S q − 2q + q S S S
t S2  ˆ t+( j−1)·
t S  t+( j−2)
t M t+( j−1)
t S t+ j·
t S S qL + SL 2 −2q L + qL
t L t+
t ˆ qt+
t = 0 (11) +K S

S

The stability will be evaluated by means of the eigenvalues of a generalized eigenvalue problem linking all the states during the current major cycle to that during the previous major cycle. This can be written in several ways to include

Comput Mech (2008) 41:269–277

271

displacements, velocities and accelerations. Similar to references [13,14], only the eigenvectors of displacements will be used here and the case of two-step sub-cycles is used as an illustration. Thus we can rewrite Eqs. (10) and (11) as follows: ⎤⎡ ⎡ t ⎤ ˆt ˆ M M LS LL qt+2
t L 2 ⎥ ⎢ 4
t 2
t ⎣ ⎦ ⎣ ˆt ˆt ⎦ MS L M t+
t SS qS 2 4
t 2 ⎤⎡ ⎤ ⎡
tt ˆt ˆ M M 1 ˆt LS LL qt − K 2 L ⎥⎣ L ⎦ ⎢ 4
t 2
t 2 = −2⎣ ˆ t ⎦ ˆt MS L M SS − 21 Kˆ St qtS 2 2 4
t
t ⎡ t ⎤⎡ ⎤ ˆt ˆ L M M LS LL qt−
t L 2 ⎢ 4
t 2 ⎥
t ⎣ ⎦ +⎣ ˆ t (12) ˆt ⎦ MS L M t−
t S SS q S 4
t 2
t 2 t+
t   ˆ M SS t+
t t+2
t t q − 2q + q S S S
t S2  ˆ t+
t  M + S L 2 qt−2
t − 2qtL + qt+2·
t L L 4
t ˆ + Kt+
t qt+
t =0 (13) S S

We can obtain ⎡ 4I 4 A − 8I ⎢0 4D ⎢ ⎣0 4I 0 0 ⎡ 4I ⎢ −8D =⎢ ⎣ 4 A − 8I 0

In Eqs. (12) and (13), The general mass matrix and the stiffness matrix are time-dependent and state-dependent. Hence ˆ t and K ˆ t+
t
= K ˆ t . For convenience, in terms ˆ t+
t
= M M of a widely used instantaneous structure method, the genˆ and K, ˆ can be regarded as constant matrices eral matrices, M within a major cycle. According to Eqs. (12) and (13), and assuming there are only two sub-cycles in one major step, we can obtain the relationship of the displacements in the current major step and the displacements in the previous major step as following.

Thus, stability analysis of the integral process is equivalent to a generalized eigenvalue problem as following.

⎡

ˆ SS M
t 2

⎢ ⎢ ⎢0 ⎢ ⎢ ⎢ ⎢0 ⎣ 0

ˆ LS M
t 2 ˆ SS M
t 2

ˆ SL M 4
t 2 ˆ LL M 4
t 2 ˆ SL M 4
t 2

ˆL − K

0

0

I

ˆS− K

⎡

ˆ SS 2M
t 2

ˆ SS M ⎢
t 2 ⎢ ⎢ − 2Mˆ L S ⎢
t 2

=⎢ ⎢ ˆ − ⎢K ⎣ S 0

ˆ SS 2M
t 2

0 ˆ LS M
t 2 ˆ SS M
t 2

0

LL

ˆ LL M 2
t 2

ˆ

MS L − 2
t 2

0 0 0 I

ˆ SL M 4
t 2 ˆ LL M 4
t 2 ˆ SL M 4
t 2

0

Define four symbols below. 2 ˆ −1 K ˆ A=M SS S
t ˆ L
t 2 ˆ −1 K C=M

⎤

ˆ

MS L − 2
t 2

ˆ SL ˆ −1 M B=M SS ˆ LS ˆ −1 M D=M LL

⎤ ⎥⎡q ⎥ S+2 ⎥⎢ ⎥ ⎢ q S+1 ⎥ ⎥ ⎥⎣ ⎥ q L+2 ⎦ ⎥ ⎦ qL

⎤ ⎤ ⎥⎡q ⎥ S ⎥⎢ ⎥ ⎢ q S−1 ⎥ ⎥ ⎥⎣ ⎥ qL ⎦ ⎥ ⎦ q L−2

(14)

⎤ ⎤⎡ −2B q S+2 ⎥ ⎢ 4C − 2I ⎥ ⎥ ⎢ q S+1 ⎥ ⎦ ⎣ −2B q L+2 ⎦ qL I ⎤ ⎤⎡ 0 0 B qS ⎥ ⎢ 4D 0 I ⎥ ⎥ ⎢ q S−1 ⎥ ⎦ ⎣ qL ⎦ 4I 0 B q L−2 0 I 0

B I B 0

(15)

Define symbols below. ⎡ ⎤ 4I 4 A − 8I B −2B ⎢0 4D I 4C − 2I ⎥ ⎥ AM = ⎢ ⎣0 ⎦ 4I B −2B 0 0 0 I ⎡ ⎤ 4I 0 0 B ⎢ −8D 4D 0 I ⎥ ⎥ BM = ⎢ ⎣ 4 A − 8I 4I 0 B⎦ 0 0 I 0 T

q2 = q S+2 q S+1 q L+2 q L

T q0 = q S q S−1 q L q L−2

AM q = λBM q

(16)

The amplification matrix is useful for studying stability of the algorithm applied to small models with a few degrees of freedom. We consider firstly the case of only one degree of freedom in the S region and one degree of freedom in the L region in this section. Then the symbols, A, B, C and D, become scalars. Expressions of the eigenvalues have been obtained by symbolic algebra for the two sub-cycles case. These expressions are extremely lengthy and will not be reproduced here. Computational experience suggests that the case of two degrees of freedoms system can be easily generalized to be applied to the case of multi-degrees of freedoms system [11]. A simple flexible multi-body system, which is shown in Fig. 1, is used as an example for stability analysis of the sub-cycling. In Fig. 1, a rigid rod, O A, is linked to the origin of the inertial coordinate system, O, at one end. It rotates around the origin in the X-Y plane and its general coordinates is the rotating angle denoted as θ . This variable will be integrated with a larger step size during the sub-cycling process. A lubricous sliding block, B, is mounted on the rod and linked to the origin by means of a spring and it can only oscillate along the axis of the rod. Its general coordinate is the oscillation displacement denoted as xi . This variable will be integrated with a smaller step size within a major cycle during the sub-cycling process. This model will be used as

123

272

Comput Mech (2008) 41:269–277

Fig. 1 Rotational rod-spring-sliding block system with two degrees of freedoms Fig. 3 Stability of the sub-cycling process with three sub-cycles in one major cycle (0 = stable. 1 = unstable.)

Figure 3 describes the stability of the sub-cycling with three sub-cycles in a major cycle of the same model mentioned above. We can see that the boundary of the instability is larger than that in Fig. 2. The plot also shows that along with the number of sub-steps in a major cycle is increased, the rapidly changing variables will gradually participate to affect the stability of the sub-cycling. Such a phenomenon becomes more obvious as the number of the sub-steps in a major cycle is continuously increased. 3 Numerical examples

Fig. 2 Stability of the sub-cycling process with two sub-cycles in one major cycle (0 = stable. 1 = unstable.)

a numerical example in the following section and its parameters will be described there. Stability of this model in the case of two sub-cycles is plotted in Fig. 2. The axes defined in Fig. 2 varies with courant numbers  S = ωn S
t and  L = ωn L 2
t, where ωn S and ωn L are the natural frequencies of each sub-region with scale of Hz, 2 ˆ ˆ −1 K ˆ −1 ˆ ωn2 S = M SS S , ωn L = M L L K L . A value of zero on the vertical axis represents that the process is stable, and one represents that the process is unstable. The plot shows a narrow ridge of unstable time-steps above the central difference stability limitation. It is interesting to note that the instability is mainly determined by the major step size. Although we still cannot conclude that the major cycle’s step size must be the step size of the rigid motion variables in a complicated FMS, we can certainly confirm that the instability is mainly affected by the slowly rotational variable in this simple system.

123

In this section, several numerical examples will be performed to validate the availability and efficiency of the presented sub-cycling algorithm. Readers will find that the computation process can hold obvious enhancement of the computational efficiency despite a little decrease of the computational accuracy. 3.1 Simulation of a rotational rod-spring-sliding block system The model was already displayed in Fig. 1. The rod is assumed as a lubricious rigid body. Its length is l with mass of m A ; the sliding block, B, can be regarded as a particle with mass of m B ; stiffness of the spring is k without mass and damping. τ is the torque to drive the rod rotating around the origin. The rotating angle of the rod and the dynamic response of the sliding block, θ and xi , need be estimated. As shown in Fig. 1, set XOY as the inertial coordinate system. The floating frames Oi xi yi and O j x j y j are linked to the static equilibrium position of the sliding block. Oi xi and O j x j always keep consistent with the axis of the rod. Thus

Comput Mech (2008) 41:269–277

273

ˆ i (q˙ i , qi , q˙ d , qd , t) ˆ i (qi , qd , t)q¨ i = Q M

(17)

There into   j j T qi = [θ xi ]T qd = Rxi R iy Rx R y  T   T −1 ˆ i = Mii − Mid C−1 Mdi −Mdd C−1 M qd Cqi − Cqd Cqd qd Cqi   T  T −1 ˆ i = Q∗ − Mid C−1 Q Q∗Fd − Mdd C−1 qd γ − Cqi Cqd qd γ Fi 
1  2 2 0 3 m A l + m B xi Mii = 0 mB
 −m B xi sin θ m B xi cos θ 0 0 Mid = m B cos θ m B sin θ 0 0

Q∗Fi

 =

T 0




 −


0

500 400 300 200 100 0

0

1

2

3

4

5 6 Time(s)

7

8

9

10

Fig. 4 Results of θ with and without sub-cycling (five sub-steps in one major cycle)

0.2 no subcycling subcycling

0.15 0.1

⎤T 0 ⎥ 0⎥ ⎥ 1⎥ ⎦ 0

0.05 0 -0.05 -0.1 -0.15 -0.2 0

1

2

3

4

5

6

7

8

9

10

Time(s)



kxi

no subcycling subcycling

600

Displacement(m)

T Mdi = Mid ⎤ ⎡ 0 0 mB 0 0 ⎥ mB 0 ⎢0 Mdd = ⎣ 0 0 mA 0 ⎦ 0 0 0 mA ⎡l 2 sin θ ⎢ l ⎢ − 2 cos θ Cqi = ⎢ ⎢ −(R j − R i ) sin θ + (R j − R i ) cos θ x y ⎣ x y j j (Rx − Rxi ) cos θ + (R y − R iy ) sin θ ⎡ ⎤T 0 0 1 0 0 0 1 ⎢0 ⎥ Cqd = ⎣ − cos θ − sin θ cos θ sin θ ⎦ − sin θ cos θ sin θ − cos θ

700

Rotational angle(degree)

the rod, the spring and the sliding block constitute a simple flexible multi-body system. In terms of the method mentioned in reference [16], the condensed equations can be established as following with that qi represents the independent variable vector and qd represents the dependent variable vector.



R˙ xi θ˙ xi cos θ + R˙ iy θ˙ xi sin θ + R˙ xi x˙i sin θ − R˙ iy x˙i cos θ R˙ xi θ˙ sin θ + R˙ iy θ˙ cos θ − θ˙ 2 xi
     −xi cos θ θ˙ − x˙i sin θ R˙ xi + −xi sin θ θ˙ + x˙i cos θ R˙ iy + 2xi x˙i θ˙ − mB − sin θ θ˙ R˙ xi + cos θ θ˙ R˙ iy  ⎤ ⎡ −xi cos θ θ˙ − x˙i sin θ θ˙ − sin θ θ˙ x˙i ⎢ −x sin θ θ˙ + x˙ cos θ θ˙ + cos θ θ˙ x˙ ⎥ i i i ⎥ ⎢ = −m B ⎢ ⎥ ⎦ ⎣0

Fig. 5 Results of xi with and without sub-cycling (five sub-steps in one major cycle)

−2m B

Q∗Fd

⎡

0 l ˙˙ 2 θ θ cos θ l ˙˙ 2 θ θ sin θ

⎤

⎢ ⎥ ⎢ ⎥ ⎥ γ = −⎢ ⎢ −θ˙ sin θ R˙ i + θ˙ cos θ R˙ i + θ˙ sin θ R˙ j − θ˙ cos θ R˙ j + S θ˙ ⎥ 1 x y ⎣ ⎦ x y j j −θ˙ cos θ R˙ xi − θ˙ sin θ R˙ iy + θ˙ cos θ R˙ x + θ˙ sin θ R˙ y + S2 θ˙ j j S1 = R˙ xi sin θ − R˙ iy cos θ − R˙ x sin θ + R˙ y cos θ   j j − θ˙ · −(Rx − Rxi ) cos θ − (R y − R iy ) sin θ j j S2 = − R˙ xi cos θ − R˙ iy sin θ + R˙ x cos θ + R˙ y sin θ   j j + θ˙ · −(Rx − Rxi ) sin θ + (R y − R iy ) cos θ

From Eq. (17), one can easily note that both the general mass matrix and the general force are time-dependent and state-dependent. Thus the mass matrix and the general force are changing along with the states of the variables during the sub-steps within a major cycle. This can be dealt with by means of the instantaneous structure method mentioned before. The parameters used for simulation are described below. The initial positions of center of mass of the rod and the sliding block are (0,0) relative to the float frame. Mass of the rod is 2.0 kg and mass of the sliding block is 5.0 kg. The driving torque is 100 Nm/s. The stiffness of the spring is 1,000 N/m. length of the rigid rod is 2 m. Results of θ and xi computed by means of the presented sub-cycling algorithm and the central difference algorithm

123

274

Comput Mech (2008) 41:269–277

700 no subcycling subcycling

Rotational angle(degree)

600 500 400 300 200

Fig. 8 A beam with one end jointed at the origin point of the inertial coordinate system

100 0

0

1

2

3

4

5

6

7

8

9

10

Time(s) Fig. 6 Results of θ with and without sub-cycling (ten sub-steps in one major cycle) 0.2 no subcycling subcycling

0.15

Displac ement(m)

0.1

a slowly changing variable in the system, This phenomenon is reasonable because the slowly changing variables are only updated at the start of each major cycle. In Figs. 5 and 7, we can see that although the variation trend of xi with five sub-steps in one major cycle are same as that with ten substeps in one major cycle, they are different at each discrete step. Compare the results with sub-cycling to that without sub-cycling, it is obvious that the precision of the results of the rapidly changing variable is a little decreasing along with that the sub-step number in one major cycle is increasing.

0.05

3.2 Simulation of a rotational beam

0 -0.05 -0.1 -0.15 -0.2

0

1

2

3

4

5

6

7

8

9

10

Time(s)

Fig. 7 Results of xi with and without sub-cycling (ten sub-steps in one major cycle)

are shown in Figs. 4 to 7 below. The discontinuous lines in Figs. 4 to 5 are results carried out by sub-cycling with five sub-steps in one major cycle (5:1). The discontinuous lines in Figs. 6 to 7 are results carried out by sub-cycling with 10 sub-steps in one major cycle (10:1). And the continuous lines are results carried out by the central difference method without sub-cycling. The time cost of the 5:1 sub-cycling is 3.57 s and the time cost of the central difference method is 7.98 s. The ratio of these two time costs is 44.7%. The time cost of the 10:1 sub-cycling is only 2.89 s. The computational efficiency is enhanced over double relative to the central difference method. Comparing Figs. 4 to 6, we can see that results of θ with five sub-steps in one major cycle are almost same as that with ten sub-steps in one major cycle. Since θ is regarded as

123

Figure 8 shows a flexible beam with elastic modulus EX = 7.22e2, Poisson’s ratio PRXY = 0.34, and density DENS = 2e-9. Its section is 10 × 25 mm2 and length is 300 mm. One end of the beam is linked at the origin of the inertial coordinate system by means of a revolute joint. A torque drives it rotating around the origin at an invariable angle velocity in the X-Y plane. The FMD equation of this model is established based on the mode basis coordinate system. Natural frequencies and modes of the beam can be obtained by means of an FEM analysis procedure. Generally speaking, the main contribution to deformation of the flexible body of an FMD problem is usually produced by the first several orders of modes. Various methods can be adopted to bite off the higher frequency modes [17,18]. For reducing the dimension of an FMD model, only the first five orders natural modes in X-Y plane are applied to be the mode basis of the FMD formulae in this paper. Such a truncation of the higher frequency modes will not take excessive infection of the comparison results of the sub-cycling and that without sub-cycling. Table 1 shows the first five orders of the natural frequencies of the beam. Trajectory of center of mass of the elastic beam, which are carried out by means of the methods with sub-cycling and without sub-cycling respectively, are depicted together in Fig. 9. The velocities and the angular velocities of mass of the elastic beam are depicted in Figs. 10 and 11 respectively. And these are measured referring to the

Comput Mech (2008) 41:269–277

275

Table 1 The natural frequencies corresponding to the first five orders modes in X-Y plane Order of the natural modes

Natural frequencies (Hz)

1

131.2

2

349.3

3

653.9

4

1, 025.3

5

1, 444.5

Fig. 11 Angular velocity of center of mass of the elastic beam in rotating

Fig. 9 Displacement of center of mass of the elastic beam in rotating

To improve the computational precision, local refinement is done at the bottom of the blades. Parameters of the model are as following (Units: tonne/N/mm/s). The material of the propeller is the aluminum alloy. E X = 7e4, P R = 0.3, D E N S = 2.78e − 6. Diameter of the airscrew is D = 900 mm, rotate speed is 8,000 rpm. The elements number of the propeller is 20,700. Thereinto, the elements number of the joints is 5,760, the elements number of the sleeve is 1,980, and the elements number of the blades is 12,960. The dynamic equation of the propeller can be established based on the free mode information obtained from mode analysis of FEM. Then we can simulate the large range overall motion of the airscrew by means of the sub-cycling and the central difference method respectively. The dynamic stress and the motion trajectory of an element on the joint are described in Figs. 13 and 14. The time cost of the sub-cycling is 167 s and the time cost of the central difference algorithm is 283 s. The computational efficiency is enhanced about 70%. The compared results show the suitable precision of the subcycling.

Fig. 10 Velocity of center of mass of the elastic beam in rotating

origin point of the inertial coordinate system. The discontinuous line in Fig. 9 is the results carried out by sub-cycling and the continuous line is the results carried out by the central difference algorithm without sub-cycling. The time cost of the integral process with sub-cycling and without sub-cycling are 5.5 and 12.7 s, respectively. In other words, the time cost of the sub-cycling integral is only about 43.3% of that of a single-step-cycle integral. 3.3 Simulation of propeller of a jet engine The FEM model and the local mesh of the propeller of a jet engine are displayed in Fig. 12. Experience indicates that the dangerous position is the joint of the blades and the sleeve.

4 Conclusions This paper attempts to analysis the stability of the centraldifference-based sub-cycling method for FMD. Omitting the external forces, an FMD formula can be expressed in the form as Eq. (5). Then regard the general force related with the quadratic term of the velocity as the damping, and partition variables to a group of slowly changing variables and a group of rapidly changing variables, the simplified form of Eq. (5) can be written as the block-matrix form in Eq. (7). On the basis of the integral approximation operator method introduced in reference [15], an amplification matrix, which represents the relationship of the variables, xn+1 and xn , can be established. Subsequently the stability analysis can be solved as a generalized eigenvalue problem.

123

276

Comput Mech (2008) 41:269–277

Fig. 12 FEM model of the airscrew and the local mesh

80 no subcycling subcycling

Dy namic stress(Mpa)

75 70 65 60 55 50 45

0

0.005

0.01

0.015

0.02

Time(s) Fig. 13 Dynamic stress result of one element at the joint during rotation

Acknowledgments The authors gratefully acknowledge the financial support of pre-research project from the 14th academy of the electronic science and technology corporation group, People’s Republic of China.

9 no subcycling subcycling

Rotational angle(radian)

8

limitation. And we can conclude that the instability is mainly determined by the step size of the major cycle in this case. Figure 3 shows a larger boundary of the instability than that in Fig. 2. This phenomenon indicates that along with increase of the number of sub-steps in one major-cycle, the rapidly changing variables in flexible multi-body system will gradually participate in the influential factors of the stability of the sub-cycling process. Such a phenomenon is more obvious along with increase of the number of the sub-steps in one major-cycle. Several numerical examples are performed to validate the availability and efficiency of the presented sub-cycling algorithm. The computational efficiency can be greatly enhanced although the computational accuracy is decreased a little. And the accuracy of rapidly changing variables is decreasing while the sub-step number in one major-cycle is increasing.

7

References

6 5 4 3 2 1 0 0

0.005

0.01

0.015

0.02

Time(s) Fig. 14 Rotational angle result of one element at the joint during rotation

Based on a simple flexible multi-body system, the stability of the sub-cycling algorithm was numerically calculated and displayed in Figs. 2 and 3. Figure 2 shows a narrow ridge of unstable time-steps above the central difference stability

123

1. Belytschko T, Yen HJ, Mullen R (1979) Mixed methods for time integration. Comput Methods Appl Mech Eng 17(18):259–275 2. Miao JC, Zhu P, Shi GL, Chen GL; Study on sub-cycling algorithm for flexible multi-body system—part I: integral theory and implementation flow chart. Comp. Mech. doi: 10.1007/s00466007-0183-9 3. Gravouil A, Combescure A (2001) Multi-time-step explicit– implicit method for non-linear structural dynamics. Int J Numer Methods Eng 50:199–225 4. Hughes TJR, Liu WK (1978) Implicit–explicit finite elements in transient analysis: stability theory. J Appl Mech 45:371–374 5. Belytschko T, Mullen R (1978) Stability of explicit–implicit mesh partitions in time integration. Int J Numer Methods Eng 12: 1575–1586 6. Belytschko T, Lu YY (1992) Stability analysis of elemental explicit–implicit partitions by Fourier methods. Comput Methods Appl Mech Eng 95:87–96 7. Liu WK, Lin J (1982) Stability of mixed time integration schemes for transient thermal analysis. Numer Heat Transf 5:211–222 8. Belytschko T, Smolinski P, Liu WK (1985) Stability of multi-time step partitioned integrators for first-order finite element systems. Comput Methods Appl Mech Eng 49:281–297

Comput Mech (2008) 41:269–277 9. Smolinski P, Sleith S, Belytschko T (1996) Stability of an explicit multi-time step integration algorithm for linear structural dynamics equations. Comput Mech 18:236–244 10. Farhat C, Ivelli L, Geradin M (1993) On the spectral stability of time integration algorithms for a class of constrained dynamics problems. In: AIAA 34th Structural Dynamics Meeting 11. Daniel WJT (1997) The subcycled Newmark algorithm. Comput Mech 20:272–281 12. Daniel WJT (1998) A study of the stability of subcycling algorithms in structural dynamics. Comput Methods Appl Mech Eng 156:1–13 13. Daniel WJT (1997) Analysis and implementation of a new constant acceleration subcycling algorithm. Int J Numer Methods Eng 40:2841–2855

277 14. Daniel WJT (2003) A partial velocity approach to sub-cycling structural dynamics. Comput Methods Appl Mech Eng 192: 375–394 15. Bathe KJ, Wilson EL (1976) Numerical methods in finite element analysis. Prentice Hall, Englewood Cliffs 16. Haug EJ (1989) Computer-aided kinematics and dynamics of mechanical systems. Allyn and Bacon, Boston 17. Hughs PC (1980) Model identities for elastic bodies with application to vehicle dynamics and control. J Appl Mech 47(1):177–184 18. Skelton RE, Hughs PC (1980) Modal cost analysis for linear matrix-second-order system. J Dynam Syst Meas Control 102: 151–158

123