Comput Mech (2008) 41:257–268 DOI 10.1007/s00466-007-0183-9

ORIGINAL PAPER

Study on sub-cycling algorithm for flexible multi-body system—integral theory and implementation flow chart J. C. Miao · P. Zhu · G. L. Shi · G. L. Chen

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

Abstract A sub-cycling integration algorithm (or named multi-time-steps integration algorithm), which has been successfully applied to FEM dynamical analysis, was firstly presented by Belytschko et al. (Comput Methods Appl Mech Eng 17/18:259–275, 1979). However, the problem of how to apply this type of algorithm to flexible multi-body dynamics (FMD) problems still lacks investigation up to now. Similar to the region-partitioning method used in FEM, this paper presents a central-difference-based sub-cycling integral method by decomposing the variables of an FMD equation into several groups and adopting different integral step sizes to each group of the variables. Based on the condensed form of an FMD equation, a group of common update formulae and a sub-step update formula, which constitute the sub-cycling together, are established in the paper. Furthermore, an implementation flowchart of the sub-cycling is presented. Stability of the sub-cycling will be analyzed and numerical examples will be performed to verify availability and precision of the sub-cycling in part II of the paper. Keywords Flexible multi-body system (FMS) · Flexible multi-body dynamics (FMD) · Region-partitioning · Sub-cycling algorithm · Implementation procedure

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]

1 Introduction Flexible multi-body systems (FMS) can be applied in various domains such as design and control of aviation and space flight, vibration reduced of automobile, design and control of robots. In these domains, accurate and efficient computation of the flexible bodies’ deformation during a large overall motion is an important issue for design and control of the system. In the past 30 years, significant developments have been achieved in studies on the FMS. Researchers are interested in three domains of studies on the FMS: modeling methods, computational strategies, and experiment investigation. There are various methods for modeling flexible multibody dynamics (FMD), which can be classified by means of modeling principles, selections of reference frames and flexible body discrete methods. No matter what means is used, the final formulae of FMD are similar, which can be written in a combined form. Mq¨ + Kq + CqT λ = QF + QV

(1)

C(q, t) = 0

(2)

Thereinto, q is a general variable vector, which need be estimated by solving Eqs. (1) and (2). M is a general mass matrix. K is a general stiffness matrix. λ is the Lagrange multiplier vector. QF and QV are external forces and internal forces respectively. C(q,t) is constrained function of the system and constituted by a group of algebraic equations. Cq is the Jaccobi matrix of C(q,t). Combination of Eqs. (1) and (2) is a kind of differential and algebraic mixed equations, generally called differentialalgebraic equations (DAEs). There are two ways to solve the DAEs. One is to transform the DAEs to an augmented form and solve it as an algebraic equation group. The largest number of variables of an augmented FMD equation needs

123

258

to be assessed at each step. However, the equivalent stiffness matrix is a sparse matrix and this characteristic is an advantage for computation of the equation [2]. The integral process is performed at the acceleration level and easily meets a constraints violated problem. The Baumgarte stabilizing technique is a widely used method to deal with this case [3,4]. Due to the system constraints, the general variables in FMD equations are not entirely independent. The general variables can be separated into one group of dependent variables and another group of independent variables. The dependent variables are expressed as a function of the independent variables. Therefore, the augmented FMD equation can be transformed to a condensed form and a pure differential equation can be obtained [5]. Both types of the FMD equations mentioned above can be computed by means of various temporal domain integral techniques such as explicit integral, implicit integral and explicit-implicit mixed integral methods [6–17]. In addition, researchers attempted some other methods, such as recursive solution procedures, parallel computational strategies, object-oriented strategies, computerized symbolic manipulation and adaptive approximation strategies, to improve computational efficiency of the FMD problems. The principle of the sub-cycling, which was firstly proposed by Belytschko et al. [1], is to separate the elements or the nodes of an FEM into different region in terms of their sizes or stiffness and adopt different integral step sizes for each domain. Neal and Belytschko [18] developed it into a more general form, of which the ratio of the sub-step-size and the major-step-size need not be an integer, and proposed a numerical example to prove the efficiency enhancement of the sub-cycling. However, Daniel [19] pointed out that it is unstable in a classical sense. Belytschko and Lu [20] proposed a central-difference-based sub-cycling based on a constant acceleration assumption, which is proved to be unstable in reference [21]. And a modified version developed by Daniel [22] based on the Belytschko–Lu method obtained a so-called statistical stability. The sub-cycling algorithms proposed by Smolinski [23,24] and Daniel [25,26] are absolutely stable in the sense of energy conservation but are unstable in the sense of momentum balance. Gravouil [27,28] presented a more complicated sub-cycling based on the Newmark integral formula and was proved to be stable but energy dissipation. Daniel [29] developed a sub-cycling, which is similar to that of Smolinski [30]. It was realized based on the central difference method and the generalized Alpha method, respectively. Thereinto, the sub-cycling based on the central difference is very simple for code realization and holds a good computational precision. But it is still only statistical stable. The sub-cycling based on the generalized Alpha method is absolutely stable but is a little sacrifice of the precision and

123

Comput Mech (2008) 41:257–268

is less elegant for code realization. Recently, Prakash and Hjelmstad [31] adopted an FETI method to decompose the FEM mesh region. The algorithm was modified based on the Gravouil-Combescure method. It reduces computational cost than the GC method significantly and holds several technique characteristics such as unconditional stability, energy conservation and no dissipation. Cavin and Gravouil [32] proposed a sub-cycling related mesh refinement method, which is called space-time automatic refinement (STAR) method. For the STAR method, the domain, where the time steps varies, changes with time to adhere to the physical phenomena. Since both the sub-cycling and the parallel computational strategies involve into decomposition of the variables, people are interesting in combination of these two techniques. Krysl and Bittnar [33], Kale and Lew [34] proposed parallel computational algorithms related with sub-cycling and expected to obtain higher computational efficiencies by means of utilizing the advantages of the two methods adequately. Besides being applied in FEM structure dynamic analysis, the sub-cycling algorithms have been applied to other disciplines. Zhen Yangming et al. [35,36] have attempted to analyze a microwave circuit by means of the sub-cycling and enhanced the computational efficiency nearly 100% with properly accuracy. Wasfy and Noor [37] pointed out that although the area of FMD maybe benefited from general sub-cycling implicit and explicit solution procedures, which include an algorithm for modeling general constraints, such procedures have not yet been presented in literatures. Similar to the region-partitioning method in FEM [38], this paper presents a central-difference-based sub-cycling method by partitioning the variables into several groups in terms of their frequencies of oscillation and proposes its implementation flowchart in detail. This paper is organized in the order as following. The principle of a central-differencebased sub-cycling algorithm for FEM is detailedly explained in Sect. 2. Then a sub-cycling method, which can be used to analyze the FMD problem, is developed in Sect. 3. Its flowchart is proposed in Sect. 4. Stability analysis and numerical examples will be presented to validate availability and conditional stability property of the sub-cycling algorithm for the FMD in part II of the paper.

2 The basic principle of sub-cycling in FEM A sub-cycling integral method involves different time step sizes in different zones of the model, which is depicted in Fig. 1. Different from the time mixed methods [39], which consist of different zones with different time integral techniques applied, but with a single time step size defined for the whole structure, an FEM model can be partitioned into several regions in terms of different sizes of elements in a sub-cycling process and the integral efficiency can be

Comput Mech (2008) 41:257–268

259

Fig. 1 Region partition principle of an FEM model

significantly enhanced with acceptable accuracy. The principle of a conventional central-difference-based sub-cycling algorithm for FEM is explained in detail. A typical structure dynamics formula can be written as following. M¨x + C˙x + Kx = F(t)

(3)

In Eq. (3), M is a general mass matrix. C is a general mass matrix. K is a general mass matrix. F(t) is the external forces. x are the unknown variables and the super-dot represent differential operations to the time variable. If element sizes or their stiffness in an FEM model are extremely uneven, which means that there are some local tiny elements or very rigid elements in the FEM model, the critical step sizes of different parts are greatly uneven. In this case, the sub-cycling algorithm is a desirable selection for an efficient integral process. There are two types of cycles within one integral step in a sub-cycling procedure, one of which can be called the major step and the other the sub-step. The major-step is usually of a larger step size and the sub-step is of a smaller step size. Corresponding to these two types of cycles, two kinds of update operations need to be performed within one integral step, one is named the common update and the other the sub-step update. The sub-cycling algorithm is usually developed by means of modification of one kind of the conventional algorithms such as the central difference method [40]. In order to simplify the problem, we assume that there are only two types of regions in an FEM model. Then the unknown variables can be set to two types of difference expressions based on the original central difference formulae. x¨ ts = x˙ ts = x¨ tL = x˙ tL =

 1
s s x − 2xts + xt+t t 2 t−t  1
s s −xt−t + xt+t 2t  1  L L L x − 2x + x t t+T T 2 t−T  1  L L −xt−T + xt+T 2T

(4) (5) (6) (7)

Thereinto, x S and x L are called rapidly changing variables and slowly changing variables, respectively. t is step size of x S and T is step size of x L . And T is an integer multiplier of t as following. T = mt

(8)

Corresponding to Eqs. (4)∼(7), Eq. (3) can also be written in a block matrix pattern as following. To simplify the problem, the damping term is omitted.    S    S  S  K SS K S L x¨ x F (t) M SS 0 + = L L T F L (t) x¨ x KS L KL L 0 ML L (9) Substituting Eqs. (4) ∼(7) into Eq. (9), a common update formula can be developed as following.  s  1 M SS s xt+t = F S (t) − M SS xt−t − 2xts 2 2 t t s −K SS xt − K S L xtL

1 ML L L L L L x x = F (t) − M − 2x L L t t+T t−T T 2 T 2 T s −K S L xt − K L L xtL

(10)

(11)

After the common update, the sub-cycling procedure enters the sub-step update stages. As variable x S changed from s s xt+( j−1)t to xt+ jt , the sub-step update formula is similar to Eq. (10).

M SS s M SS s S s xt+( j−2)t −2xt+( xt+jt = Ft+( j−1)t − j−1)t 2 2 t t s L (12) −K SS xt+( j−1)t − K S L xt+( j−1)t Hereinto, 2 ≤ j ≤ m. It is easy to note that the slowly changing components, L xt+( j−1)t , is unknown at this sub-step. The known results L , respectively. Thus, from the second are only xtL and xt+T step of the sub-steps to the end of this major-step, interaction of the nodes on interface are always undiscovered. This problem can be settled in terms of the well-known trapezoid rule. Assuming x˙ L within a major-step is a constant. Then, The following interpolation relationship can be obtained based

123

260

on the trapezoid rule.  j j L L = 1 − xt+( x L + xt+T j−1)t m t m

Comput Mech (2008) 41:257–268

(13)

Substituting Eq. (13) into Eq. (12), the sub-step update formula can be now modified as following. M SS s M SS s S xt+( j−2)t x = F − t+ jt t+( j−1)t 2 t 2

t s s −2xt+( j−1)t − K SS xt+( j−1)t  j j L L 1− • xt + • xt+T (14) −K S L m m Formulae (10), (11) and (14) constitute the common update formulae and the sub-step update formula. This type of sub-cycling can be applied to a structure dynamic analysis problem.

3 Sub-cycling algorithm for flexible multi-body dynamics Compared with a structural dynamics model, the flexible multi-body dynamic formula holds several features as following: (1) The rigid motion variables of a flexible body change slowly while the elastic deformation variables of the flexible body change rapidly. (2) There are strong coupling between the rigid motion variables and the elastic deformation variables. (3) Different external forces maybe excite elastic deformations with different frequencies. Just like an FEM model with local extremely tiny elements or very hard elements, the FMD is naturally suitable for sub-cycling because of co-existence of the rapidly changing variables and the slowly changing variables. We can expect to integrate the rigid motion variables in larger step sizes and update the elastic deformation variables in smaller step sizes. Assuming the deformation of the flexible body is under the elastic limitation, a model of the FMD can be established by adopting floating frames and general mode coordinate basis. To perform a sub-cycling integral, the DAEs formula of the FMD equation need be condensed to a pure differential equation in advance. Detailed description of modeling the FMD based on the floating frame and the mode coordinate basis and transformation of the DAEs equation to a condensed FMD equation will be described in Sect. 3.1. 3.1 Flexible multi-body dynamic model In flexible multi-body dynamics, four types of methods were frequently used, which named hypothetic mode method, finite segment method, finite element method (FEM) and

123

model synthesis method, to describe the deformation of a flexible body undergoing large overall motion. Along with progress of FEM, combination of finite element method and model synthesis method is the most widely used method for modeling an FMD. In addition, an inertial frame serves as a global reference frame for describing the motion of the multi-body system. Intermediate reference frames that are attached to each flexible component and follow the average local rigid body motion are often used. Two main types of intermediate frames are often used: floating frames and corotational frames. The floating frame follows an average rigid body motion of the entire flexible component or substructure. The floating frame formulation along with modal reduction offer the most efficient method for the simulation of flexible multi-body systems undergoing small elastic deformations and slow rotational speeds. The co-rotational frame follows an average rigid body motion of an individual finite element within the flexible component. The co-rotational formulation can handle flexible multi-body systems undergoing large deflections and large high-speed rigid body motion. In addition, if used in conjunction with an explicit solution procedure, high-speed wave propagation effects can be accurately modeled [37]. Characteristics of the floating frame and the co-rotational frame can be easily explained in Fig. 2. In this paper, the floating frame along with model synthesis method is used for modeling an FMD. Thus the sub-cycling algorithm proposed in this paper can work for linear elastic deformation of flexible bodies. For a multi-body system composed by N bodies, one of the flexible bodies is depicted in Fig. 3. we can mesh the ith (i = 1, 2, . . . , N ) flexible body, Bi , as l elements by means of FEM with lumped mass type. ρik is the vector of the kth node in the floating frame and its un-deformed locak . The translation deformation vector is tion is denoted as ρio k denoted as ui . The general coordinates of the flexible body can be defined as following.
q = Ro θ

qf

T

(15)

Thereinto,
Ro = Rx

Ry

Rz

T

(16)

is the vector of origin of the floating frame.
θ = θx

θy

θz

T

(17)

is the Eulerian angles vector of the floating frame gesture and is a pseudo-coordinate. T
q f = q1 q2 · · · qs

(18)

Comput Mech (2008) 41:257–268

261

Fig. 2 Characteristics of the floating frame and the co-rotational frame

⎤ M R R M Rθ M R f Thereinto M = ⎣ Mθ R Mθθ Mθ f ⎦ is the general mass Mf R Mfθ Mff matrix. ⎡

Fig. 3 Deformation description of an arbitrary flexible body

⎡

⎤ 0 0 0 K = ⎣0 0 0 ⎦ 0 0 Kff

is the general mode coordinates for describing the elastic deformation of a flexible body. Elastic deformation of the kth element of the flexible body can be described based on the general mode basis as following. uik = φik q f

(19)


T Here, φ = φ1 φ2 · · · φs is constituted by S ranks of natural modes, which are solution of free vibration problem of the flexible body, and regarded as the mode basis of the general deformation coordinates. Based on definition of the coordinate system and the general variables mentioned above, the free motion equation of the arbitrary flexible body can be deduced from the Lagrange equation and can be written in a matrix expression below [2]. ⎤ ⎡ ⎤⎡ R ⎤⎡ R ⎤ ¨o o 0 0 0 M R R M Rθ M R f ⎥ ⎣ ⎢ ⎥ ⎦ ⎣ Mθ R Mθθ Mθ f ⎦ ⎢ ¨ + 0 0 0 θ θ ⎣ ⎦ ⎣ ⎦ Mf R Mfθ Mff 0 0 Kff qf q¨ f ⎤ ⎡ ⎤ ⎡ Qv R QR ⎥ ⎢ ⎥ ⎢ (20) = ⎣ Qθ ⎦ + ⎣ Qvθ ⎦ ⎡

Qq

Qv f

is the general stiffness matrix. T
is the general external force Qext = Q R Qθ Qq vector.
T is the quadratic velocity term QV = Qv R Qvθ Qv f related general force. In an FMS, the bodies are not completely free; there are various kinds of constraints between bodies. Thus associated with the free motion equation of the bodies, a general function formula, which used to describe the constraints between bodies, can be written as following. C(q, t) = 0

(21)

The constraints can be inserted into the free motion equation by means of a Lagrange multiplier vector. Then we can obtain the integrated flexible multi-body system governing equation based on the general coordinates. Mq¨ + Kq + CqT λ = QF + Qv

(22)

C(q, t) = 0

(23)


T Thereinto, C(q, t) = C1 C2 Cnc . Ci (i = 1, 2, . . . , n c ) are linear independent on each other and they have two orders continuous differential coefficients of the general variables, q. λ is the Lagrange multiplier vector. Cq is the Jacobean matrix of the constraints. Equations (22) and (11) are called differential-algebraic equations (DAEs) and generally cannot be solved directly.

123

262

Comput Mech (2008) 41:257–268

Performing one differential action to the time variable t on Eq. (21), we can obtain the constraint relationship at velocity level.

U is a full rank S × S matrix. Hence qi can be sought in advance and qd can be obtained from qi . The relationship of qd and qi are guaranteed by the constraint equation.

Cq q˙ + Ct = 0

qd = f (qi , t)

(24)

Performing once more differential action to the time variable, t, the constraint relationship can be obtained at acceleration level from Eq. (24).   Cq q¨ + Ctt + 2Cqt + Cq q˙ q q˙ = 0 (25) Define symbols below.

Then we can obtain equations as following. ν = Cq q˙

(26)

γ = Cq q¨

(27)

Right now, we obtain the augmented form of FMD governing equation by combining Eqs. (22) and (27). Mrr

Mr f

CqTr

⎤⎡

q¨ r

⎤

⎡

QFr + Qvr

⎤

⎢ ⎥ ⎥⎢ ⎥ ⎢ ⎣ M f r M f f CqT f ⎦ ⎣ q¨ f ⎦ = ⎣ QF f + Qv f − Kq ⎦ γ λ Cqr Cq f 0 (28) Regarding q¨ and λ as the unknown variables, the augmented FMD equation can be numerically solved as an algebraic equation group but the entire variables should be assessed at each step simultaneously [2]. Due to the system constraints, some of the general variables in the DAEs equation are not independent. By separating the variables into dependent variables, qd , and independent variables, qi , one can obtain the condensed form of the FMD equation. Detailed inducing process will be shown below. Assuming the virtual displacements of the general coordinates satisfy the constraint equation, we can obtain the following equation. Cq δq = 0

(29)

Equation (29) is an algebraic equation group. The Jacobean matrix, Cq , can be decomposed into two sub-matrices by means of a Gauss elimination method or by other matrixdecomposing method. These two sub-matrices are denoted U and R, respectively. The variables, q, can now be separated into independent variables, qi , and dependent variables, qd . Thus Eq. (29) turns into the form below. Uδqd + Rδqi = 0

123

Then Eq. (28) can be decomposed as the form below: Mdd q¨ d + Mdi q¨ i + CqTd λ = Q∗Fd

(32)

Mid q¨ d + Mii q¨ i + CqTi λ = Q∗Fi

(33)

Cqd q¨ d + Cqi q¨ i = γ

(34)

It can be proved that Cqd is a nonsingular matrix. In terms of Eqs. (27) and (34), we can obtain the following relationship.   γ − Cqi q¨ i (35) q¨ d = Cq−1 d

  ν = −Ct , γ = −Ctt − 2Cqt − Cq q˙ q q˙

⎡

(31)

(30)

Also, from Eq. (26) and using matrix decomposition technique, an equation can be established. Cqd q˙ d + Cqi q˙ i = ν Then we can obtain a similar relationship of q˙ d and q˙ i as following.   ˙i (36) q˙ d = C−1 qd ν − Cqi q The Lagrange multiplier can be calculated from Eq. (32).

T   λ = C−1 (37) · Q∗Fd − Mdd q¨ d − Mdi q¨ i qd Substituting Eqs. (35) and (36) into Eq. (33), we obtain the final condensed FMD model, which is a pure differential equation, based on the independent variables as following [41]. ˆ i (q˙ i , qi , q˙ d , qd , t) ˆ i (qi , qd , t)q¨ i = Q M

(38)

Hereinto,

T T −1 ˆ i = Mii − Mid C−1 M qd Cqi − Cqd Cqd   × Mdi − Mdd C−1 C q qd i

T   T −1 ∗ −1 ˆ i = Q∗ − Mid C−1 Q C Q γ − C − M C γ dd qd qd qi qd Fi Fd Eq. (38) can be solved by seeking qi at each discrete time step by numerical integral. The corresponding variables q˙ i and q¨ i can be obtained after qi is known. Then the dependent variables, qd , q˙ d and q¨ d can be calculated in terms of Eqs. (31), (35) and (36). The equivalent system matrices in a condensed FMD equation are denseness. The matrix decomposition operations need to be performed at each step. The unknown variables in Eq. (38) include rigid body motion components and elastic deformation components. Usually, the rigid motion components are slowly changing

Comput Mech (2008) 41:257–268

variables and the elastic deformation components are rapidly changing variables. Coexistence and strong coupling of these two types of variables make Eq. (38) a typical stiff differential equation. Hence integral methods with a uniform step size will meet two types of problems such as computational inefficiency and accumulative integral errors of slowly changing variables.

263

S, represents a smaller step size. Thus, q¨ is expressed as a T
decomposed formation, q¨ = q¨ L q¨ S . And q¨ L is corresponding to the retained rigid body motion variables, q¨ S is corresponding to the elastic deformation variables. Then the corresponding general mass matrix and the general external force vector can be decomposed to block formations as following. 

3.2 Sub-cycling algorithm for FMD In order to solve an FMD problem by sub-cycling, the original FMD equation, which is a DAEs formula, needs to be transformed to a condensed formation in advance. This can be done by means of the process mentioned in Sect. 3.1 and the final FMD equation is expressed in Eq. (38). Since the elastic deformation variables have been transformed to the general mode coordinates and are independent to each other, Eq. (38) maybe still keeps partial rigid body motion variables and all the elastic deformation variables, which are all independent. This type of FMD formula is a pure differential equation and suitable to the sub-cycling algorithm. Assuming elastic deformation of the flexible body is under the elastic limitation, the deformation variables can be expressed based on the general natural modes basis coordinates and uncoupled. Figure 4 shows that there is no interface among the elastic deformation variables. However, two kinds of coupling phenomena still exist, one of which is coupling of the rigid body motion variables and the elastic deformation variables, the other one is coupling of the rigid body translation variables and the rigid body rotation variables. For simplification, we separate the entire variables in Eq. (38) into two groups. The retained rigid body motion variables are assigned in one group and the corresponding step size is set to a larger size, t L . The elastic deformation variables are parked in another group and the corresponding step size is set to a smaller size, t S . The subscript symbol, L, represents a larger step size and the subscript symbol,

ˆ = M

ˆt M LL

ˆt M LS

ˆt M SL

ˆt M SS



 ˆ = Q

ˆt Q L



ˆt Q S

Here, the superscript symbol, t, illustrates that the terms in the matrices are time-dependent. Thus, the condensed FMD formula can be decomposed as a block matrix formation as following. 

ˆt M LL

ˆt M LS

ˆt M SL

ˆt M SS



q¨ L q¨ S



 =

ˆt Q L



ˆt Q S

(39)

In order to perform sub-cycling integral of Eq. (39), the larger step size components, q L , and the smaller step size components, q S , are described based on the central difference rule in Eqs. (40)∼(43).  1  t−t L t+t L t q − 2q + q L L t L2 L   1 t+t L L q˙ tL = −qt−t + q L L 2t L   1 t−t S t+t S t q¨ tS = q − 2q + q S S S t S2  1  t−t S S −q S q˙ tS = + qt+t S 2t S

q¨ tL =

(40) (41) (42) (43)

Hereinto, t L = mt S , m is a positive integer. It means that there are m sub-steps within one major-step.

Fig. 4 Principle of region partitioning and coupling between variables

123

264

Comput Mech (2008) 41:257–268

Substituting Eqs. (40)∼(43) into Eq. (39), we can obtain a common update formula in the sub-cycling procedure below. ˆt  M LL t L2

L L qt−t − 2qtL + qt+t L L



 ˆt  M S S ˆt =Q − 2qtS + qt+t + L2S qt−t L S S t S  ˆ t  t−t M SL L q L L − 2qtL + qt+t L 2 t L  ˆt  M S S ˆt =Q − 2qtS + qt+t + SS2 qt−t S S S t S

t L2

L L qt−t − 2qtL + qt+t L L

(44)

(45)

(46)

(47)

Initial approximate values of the variables at time t + t L and t + t S should be predefined to start up the Newton– Raphson iteration. One way is to take the previous variables value as the initial approximate values. (1) t+t L qL (1) t+t S qS

 ˆ LL  M (k+1) t+t L = qL 2 t L

 ˆ  ˆ t − M L S (k) qt+t S Q L S 2 t S   ˆ LS M t S − 2 qt−t − 2q S S t S  ˆ LL  M t L − 2 qt−t − 2q L L t L   ˆ SS  ˆ SL  M M (k+1) t+t S (k) ˆ t (k) t+t L = Q q − q S S L t S2 t L2  ˆ SL  M t L − 2 qt−t − 2q L L t L  ˆ SS  M t S − 2q − 2 qt−t S S t S L (k+1) qt+t = L S (k+1) qt+t S

= qtL

(48)

=

(49)

123

(51)

=

(k)

(k+1) t+t L qL (k+1) t+t S qS

(52)

(53)

L − (k) qt+t L

(54)

S − (k) qt+t S

(55)

S L The terms (k+1) qt+t and (k+1) qt+t in Eqs. (54) and L S (55) are the increments of q L and q S at the kth iteration step. The criterion for stopping the iteration is expressed as following.   (k+1) t+t L (k) t+t L  qL − qL (56)  ≤ εL    (k+1) t+t S (k) t+t S  qS − qS (57)   ≤ εS

Here, ε L and ε S are the allowable error critical values of the iteration. After the common update of all variables at time t, the sub-cycling procedure enters the sub-step cycles for updating the information of the rapidly changing variables, q S . The difference formations are the same as those described in Eqs. (40)∼(43). As a result, the second sub-step update formula is obtained below. ˆ t+t S M SS t S2

qtS

(50)

Hence the implicit iteration formulae of Eqs. (46) and (47) are shown below.



 ˆ t  t−t M t+t S LS t S ˆt =0 q −Q − 2q + q S L S t S2 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 =0 −Q − 2q + q + SS2 qt−t S S S t S S +

1 = qtL + t L q˙ tL + t L2 q¨ tL 2 1 (1) t+t S qS = qtS + t S q˙ tS + t S2 q¨ tS 2

(1) t+t L qL

In Eqs. (44) and (45), the superscript t, t − t and t + t represent the current step, the previous step and the next step, respectively. In the mass matrix of a condensed FMD equation, there is a complicated coupling relationship among the larger step size variables and the smaller step size variables, q L and q S , which coexist in the general mass matrix. This is a typical inertia-coupling phenomenon. Due to this coupling, S L and qt+t coexist in the left the unknown variables qt+t L S section of Eqs. (44) and (45). Thus Eqs. (44) and (45) need to be solved in a synchronous manner. In other words, evaluation of a variable at time t +t depends on the values of other variables at time t + t. An implicit iteration method can be used to deal with this type of common update problem. Thus, the explicit central difference operation is turned into an implicit iterative operation at the common update stage. This can be done by means of the Newton–Raphson iteration method. In this case, Eqs. (44) and (45) are re-written in the form below. ˆt  M LL

In Eqs. (48) and (49), the superscript in left side of the variables denote the number of the iteration. Or we can take the two orders Taylor series expansion as the initial approximate values of the iteration.

 ˆ t+t S  t−t M t+t L SL t L q −2q +q L L L t L2  ˆ t+t S  M S − SS 2 qtS − 2qt+t (58) S t S

S S qt+2t = Qt+t − S S

Comput Mech (2008) 41:257–268

265

In Eq. (58), the general mass matrices and the general forces are time-dependant and state-dependant. It means that the current values of the general matrices and the general forces, ˆ t+t S and Qt+t S , depend on the current values ˆ t+t S , M M SS SL S S S and qt+t , at the current time step of the variables, qt+t S L and their differential coefficients. The current values of the S , and their differential coefficients were obvariables, qt+t S tained at the previous time step. However, the current values S , and their differential coefficients are of the variables, qt+t L still unknown at this stage. In order to evaluate the current values of the general mass matrices and the general forces, the constant acceleration method is used to deal with this problem. In terms of the known values of the variables, qtL L and qt+t , and their differential coefficients, we can assume L S that the current value of the valuables, qt+t , are changed L as follows. S q¨ t+t = q¨ tL L  1 1 t+t S L q˙ tL + q˙ t+t q˙ L = 1− m m L

(59) (60)

Then we can obtain the interpolation formula as following.  S qt+t = qtL + t S L

1−

 1 1 L q˙ tL + q˙ t+t m m L

Then Eq. (62) can be modified to the form below. ˆ t+( j−1)t S M SS t S2

t+ j·t S

qS

t+( j−1)·t S

= QS

−

ˆ t+( j−1)·t S M SL t L2

 M  ˆ t+( j−1)·t S t+t L SS t L − × qt−t − 2q + q L L L t S2   t+( j−2)t S t+( j−1)t S × qS − 2q S

(64)

Hereinto, the general mass matrix and the general external force can be described as following. ˆ t+( j−1)t S = M ˆ SS qt+( j−1)·t S, qt+( j−1)·t S , M SS i d

(65) t + ( j − 1)t S ˆ t+( j−1)t S = M ˆ S L qt+( j−1)·t S, qt+( j−1)·t S , M SL i d

(66) t + ( j − 1)t S ˆ t+( j−1)·t S = Q ˆ S q˙ t+( j−1)·t S, qt+( j−1)·t S , q˙ t+( j−1)·t S , Q S i i d

t+( j−1)·t S qd , t +( j − 1)t S

(61)

(67)

Finally, Eq. (58) can be modified as the form below. In ˆ t+t S , M ˆ t+t S and Qt+t S , Eq. (62), the general matrices, M SS SL S are functions of the variables.

Equation (64) is the sub-step update formula for modifying q S at the arbitrary step within one major-cycle. Equations (52)∼(55) and (64) constitute the common-update formulae and the sub-step update formula together.

ˆ t+t S (qt+t S , qt+t S , q˙ t+t S , q˙ t+t S ) M SS L S L S

3.3 Start-up the sub-cycling algorithm

S · qt+2t S t S2 S S S S S = Qt+t (qt+t , qt+t , q˙ t+t , q˙ t+t ) S L S L S t+t t+t t+t t+t t+t S ˆ (q L S , q S S , q˙ L S , q˙ S S ) M − SL t L2   t+t L t L · qt−t − 2q + q L L L

−

ˆ t+t S (qt+t S , qt+t S , q˙ t+t S , q˙ t+t S ) M SS L S L S

t 2  S  S S · qt+t − 2qt+t S S

(62)

Furthermore, at the arbitrary jth(2 ≤ j ≤ m) sub-step, we can write the interpolation formula below.

t+( j−1)·t S

qL

= qtL +

j−1  i=1

  i i L q˙ tL + q˙ t+t t S 1 − m m L (63)

Startup of the central difference method requires the known values at two previous steps of the variables. Therefore, besides the initial condition values, the values of the variables at −t L and −t S should be sought out in advance. In order to obtain values of these variables, we can calculate the initial acceleration firstly by solving Eq. (39) at time t. Then the values of the variables at time −t can be evaluated in terms of the assumption of the central difference formulae described in Eqs. (40)∼(43). According to Eq. (39), we can obtain the initial acceleration, q¨ i , below.  −1 ˆt ˆt Q q¨ it = M i i

(68)

Then the value of the variables at time −t can be evaluated by means of the hypothesis in Eqs. (40)∼(43). qit−t = qit − q˙ it t −

t 2 t q¨ 2 i

(69)

123

266

Comput Mech (2008) 41:257–268

Thus, we obtain the final formulae for startup of the central-difference-based sub-cycling algorithm. L qt−t L S qt−t S

t 2 t q¨ = qtL − q˙ tL t + 2 L t 2 t q¨ = qtS − q˙ tS t + 2 S

Fig. 5 Flow-chart of the new sub-cycling method for FMD with two sub-regions

(70) (71)

4 Implementing procedure of the sub-cycling algorithm In terms of the sub-cycling algorithm for FMD described in Sect. 3.2, we can present the program flowchart of the subcycling procedure in Fig. 5. The method can be easily generalized to S sub-regions with the assumption.

(1) Set the initial condition of the unknown variables. (2) Establish the condensed form equation of a FMS model at current time step. See equation (3-24) Calculate the general mass block matrices, M ii , M id and M dd . See equations (3-18)~(3-20) Calculate the Jacobean matrices, C q i and C q d . Calculate the values of

and

.

See equations (3-12)~(3-13) *

*

Calculate the general external forces Q Fi and Q Fd . (3) Decompose the variables into the larger step size terms and the smaller step size terms and express the system matrices in blocks form. See equation (3-25) (4) Calculate the initial value of the variables at time t. (5) Set the major-step size and the sub-step size. (6) Start up the sub-cycling procedure.

See equations (3-54)~(3-57)

(7) Common update. Set the initial values of the Newton-Raphson iteration. See equations (3-34)~(3-37) Calculate the iterative values of each variable.

See equations (3-38)~(3-41)

Convergent judgment.

See equations (3-42)~(3-43)

(8) Enter sub-step cycles for 2

j

m.

Calculate the general mass block matrices, M ii , M id and M dd . See equations (3-18)~(3-20) Calculate the Jacobean matrix C qi and C qd . Calculate the values of

and

.

See equations (3-12)~(3-13) *

*

Calculate the general external forces Q Fi and Q Fd . Calculate the interpolation values of q L at the jth sub-step. See equations (3-45)~(3-47) Calculate the values of q S at the jth sub-step. See equation (3-50) (9) Repeat step (7) to step (8) until end of the current major cycle and enter the next major cycle. (10) Repeat step (2) ~ (9) until end of the simulation time.

123

Comput Mech (2008) 41:257–268

267

t1 = t S t2 = m 2 t S ......

(72)

tm = m S t S In Eq. (72), the step size, t1 , is corresponding to the group of variables with the highest frequencies. The step size, t2 , is corresponding to the group of variables with the hypo-highest frequencies. Finally, The step size, tm , is corresponding to the group of variables with the lowest frequencies. In an FMD equation, the maximum step size, tm , is usually the step size of the rigid body motion variables and regarded as the major step size of the sub-cycling. The flowchart of the generalized sub-cycling algorithm is ignored here.

5 Conclusions This paper presents a central-difference-based sub-cycling computational strategy by partitioning the unknown variables into a group of larger step size variables and a group of smaller step size variables in terms of their frequencies of oscillation. It can be used for analyzing FMD problems with significantly enhanced efficiency. The characteristics of the parameters of an FMD model indicate that it is naturally suitable for a sub-cycling algorithm. In numerical integral of FMD equations, one can set different integral step sizes corresponding to the variables with different frequencies of oscillation. In this new algorithm for FMD, The rigid body motion variables are updated using a larger step size while the elastic deformation variables are updated using a smaller step size. The conventional central difference method is an explicit algorithm. However, due to the inertia coupling in FMD equation, the common update in the sub-cycling procedure should be calculated by implicit iteration. In the sub-steps within a major-step cycle, the current values of the slowly changing variables can be coarsely evaluated by means of the trapezoid rule. Implementation procedure of the sub-cycling algorithm was introduced in details in Sect. 4. Stability analysis and numerical examples will be proposed to verify the availability and conditional stability of the algorithm in part II of this paper. Acknowledgment 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.

References 1. Belytschko T, Yen HJ, Mullen R (1979) Mixed methods for time integration. Comput Methods Appl Mech Eng 17/18:259–275

2. Lu YF (1996) The flexible multi-body dynamics. The higher education press, Beijing 3. Baumgrate E (1972) Stabilization of constraints and integrals of motion in dynamic systems. Comput Methods Appl Mech Eng 1:1–16 4. Baumgrate E (1982) A new method of stabilization for holonomic constraints. J Appl Mech 50:869–870 5. Haug EJ (1989) Computer-aided kinematics and dynamics of mechanical systems. Allyn and Bacon, Boston, London, Sydney, Toronto 6. Wang JM (1999) Theoretic study on model of rigid-flexible coupling dynamics of flexible body and a FEM algorithm for dynamical stiffness phenomenon, Postdoctoral Thesis. Shanghai Jiao Tong University, Shanghai 7. Hong JZ, Liu YZ (1989) The state of the arts of computational dynamics of discrete systems. Adv Mech 19(2):205–210 8. Belytschko T, Hsieh BJ (1973) Non-linear transient finite element analysis with convected coordinates. Int J Numer Methods Eng 7:255–271 9. Hughes TJR, Winget J (1980) Finite rotation effects in numerical integration of rate constitutive equations arising in large deformation analysis. Int J Numer Methods Eng 15(12): 1862–1867 10. Flanagan DP, Taylor LM (1987) An accurate numerical algorithm for stress integration with finite rotations. Comput Methods Appl Mech Eng 62:305–320 11. Rice DL, Ting EC (1993) Large displacement transient analysis of flexible structures. Int J Numer Methods Eng 36:1541–1562 12. Wasfy TM (1994) Modeling continuum multibody systems using the finite element method and element convected frames. In: Proceedings of 23rd ASME Mechanisms Conference, Minneapolis, pp 327–336 13. Iura M, Atluri SN (1995) Dynamic analysis of planar flexible beams with finite rotations by using inertial and rotating frames. Comput Struct 55(3):453–462 14. Hsiao KM, Jang J (1991) Dynamic analysis of planar flexible mechanisms by co-rotational formulation. Comput Methods Appl Mech Eng 87:1–14 15. Elkaranshawy HA, Dokainish MA (1995) Co-rotational finite element analysis of planar flexible multibody systems. Comput Struct 54(5):881–890 16. Ibrahimbegovic A, AlMikdad M (1996) On dynamics of finite rotations of 3D beams. Comput Methods Appl Sci. In: 3rd ECCOMAS Comput Fluid Dyn Conf and the 2nd ECCOMAS Conf on Number Methods in Eng, pp 447–453 17. Housner JM, Wu SC, Chang CW (1988) A finite element method for time varying geometry in multibody structures. In: Proceedings of 29th Struct, Struct Dyn and Materials Conference AIAA Paper No 88-2234 18. Neal MO, Belytschko T (1989) Explicit–explicit subcycling with non-integer time step ratios for structural dynamic systems. Comput Struct 31(6):871–880 19. Daniel WJT (1998) A study of the stability of subcycling algorithms in structural dynamics. Comput Methods Appl Mech Eng 156:1–13 20. Belytschko T, Neal MO, Lu YY (1992) An explicit multi-time step integration for parabolic and hyperbolic systems. New Methods Trans Anal ASME N Y 143:25–39 21. Klisinski M, Mostrom A (1998) On stability of multitime-step integration procedures. J Eng Mech 124(7):783–793 22. Daniel WJT (1997) Analysis and implementation of a new constant acceleration algorithm. Int J Numer Methods Eng 40: 2841–2855 23. Smolinski P, Sleith S, Belytschko T (1996) Stability of an explicit multi-timestep integration algorithm for linear structural dynamics equations. Comput Mech 18:236–244

123

268 24. Wu YS, Smolinski P (2001) Multi-time step integration algorithm for structural dynamics based on the modified trapezoidal rule. Comput Methods Appl Mech Eng 187(3):641–660 25. Daniel WJT (1998) Subcycling first- and second-order generalizations of the trapezoidal rule. Int J Numer Methods Eng (42): 1091–1119 26. Daniel WJT (2001) Multi-timestep integration in computational dynamics. In: Ambrosio JAC, Kleiber MComputational aspects of nonlinear structural systems with large rigid body motion, NATO Science Series. IOS Press, Amsterdam, 27. Gravouil A, Combescure A (2001) Multi-time-step explicitimplicit method for non-linear structural dynamics. Int J Numer Methods Eng 50:199–225 28. Combescure A, Degayffier A, Gravouil A, Greffet NA (1998) A Lagrange multiplier based domain decomposition method for timedependent problems involving several time-scales. IV World Congress on Computational Mechanics 29. Daniel WJT (2003) A partial velocity approach to subcycling structural dynamics. Comput Methods Appl Mech Eng 192:375–394 30. Smolinski P, Wu YS (1998) Stability of explicit subcycling time integration with linear interpolation for first-order finite element discretizations. Comput Methods Appl Mech Eng 151:311–324 31. Prakash A, Hjelmstad KD (2004) A FETI based multi-step coupling method for newmark schemes in structural dynamics. Int J Numer Methods Eng 61:2183–2204

123

Comput Mech (2008) 41:257–268 32. Cavin P, Gravouil A, Lubrecht AA, Combescure A (2005) Automatic energy conserving space-time refinement for linear dynamic structural problems. Int J Numer Methods Eng 64:304–321 33. Petr K, Zdenek B (2001) Parallel explicit finite element solid dynamics with domain decomposition and message passing: dual partitioning scalability. Comput Struct 79:345–360 34. Kedar GK, Adrian JL (2006) Parallel asynchronous variational integrators. Int J Numer Methods Eng (published online in http://www.interscience.wiley.com). DOI: 10.1002/nme.1880. 2006 35. Zheng Y, Chu Q (2004) Multi-time-step finite difference time domain method. ACTA Electric SINICA 32(9):1504–1506 36. Zheng Y, Chu Q (2005) Analysis of microwave circuits using multitime-step finite difference time domain method. J Microw 21:11–14 37. Tamer MW, Ahmed KN (2003) Computational strategies for flexible multibody systems. Appl Mech Rev 56(6):553–613 38. Liu WK, Belytschko T (1982) Mixed-time implicit–explicit finite elements for transient analysis. Comput Struct 15(4):445–450 39. Hughes TJR, Liu WK (1978) Implicit–explicit finite elements in transient analysis: implementation and numerical examples. J Appl Mech 45:375–378 40. Bathe KJ, Wilson EL (1985) Numerical methods in finite element analysis. Scientific publishing company, Beijing 41. Haug EJ (1989) Computer-aided kinematics and dynamics of mechanical systems. Prentice-Hall, Englewood Cliffs