A NEWMARK FMD SUB-CYCING ALGORITHM
J. C. Miao
(1,2)
, P. Zhu
(1)
, G. L. Shi
(1)
and G. L. Chen
(1)
School of Mechanical Engineering, Shanghai JiaoTong University, Shanghai, China
(1)
Department of Mechanical Engineering, Shazhou institute of technology, Zhangjiagang, China
(2)
Keywords: Flexible multi-body dynamics, Condensed FMD equation, Sub-cycling algorithm, Energy balance.
Abstract: This paper presents a Newmark method based sub-cycling algorithm, which is suitable for solving the
condensed flexible multi-body dynamic (FMD) problems. Common-step update formulations and sub-step
update formulations for quickly changing variables and slowly changing variables of the FMD are
established. Stability of the procedure is checked by checking energy balance status. Examples indicate that
the sub-cycling is able to enhance the computational efficiency without dropping results accuracy greatly.
1 INTRODUCTION
Flexible multi-body system (FMS) can be applied in
various domains such as space flight, automobiles
and robots. In these domains, accurate and efficient
computation of the flexible bodies undergoing large
overall motion is important for design and control of
the system (Huang and Shao, 1996).
Conventional integration methods, such as the
Newmark algorithm, the Runge-Kutta algorithm and
the Gear algorithm and etc (Dan Negrut, et al, 2005),
were widely applied to solve FMD equations.
Sub-cycling was proposed by Belytschko T. et al
(Belytschko T.,1978). Mark et al (Mark, 1989)
applied the method to simulate an impact problem
and computational cost of the sub-cycling was only
15% of that without sub-cycling. Gao H et al. (Gao,
2005) used sub-cycling to simulate auto-body
crashworthiness and declared that the cost is only
39.3% of no sub-cycling. Tamer et al (Tamer, 2003)
pointed out that the FMD sub-cycling methods have
not yet been presented in literatures.
2 SUB-CYCLING FOR FMD
A sub-cycling is constituted by two types of cycles,
main-cycle and sub-cycle. The key for sub-cycling is
to appropriately treat with interactions of nodes on
the interface (Daniel, 1997).
2.1 Condensed Model of FMD
A flexible body is displayed in figure 1. Definition
of variables can be referenced in the literature (Lu,
1996).
Figure 1: A spatial arbitrary flexible body.
The FMD equation can be established by means
of the Lagrange formation.
()
)22(0,
)12(
−=
−+=++
""""""""
"

t
T
qC
QQλCKqqM
vFq
Thereinto, M is a general mass matrix, K is a
general stiffness matrix, C(q,t) is the constrains, Q
F
and Q
V
are general external load and general
centrifugal load. Λ is the Lagrange multiplier. By
means of two times of differentials of the constrain
equation, an augmented FMD equation can be
obtained as below.
)32( −
⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
∗
"""

C
F
q
qr
Q
Q
λ
q
0C
CM
T
Due to the constraints, variables in (2-3) are
independent. By decomposition, a condensed FMD
equation can be established as following (Haug,
1989).
207
C. Miao J., Zhu P., L. Shi G. and L. Chen G. (2007).
A NEWMARK FMD SUB-CYCING ALGORITHM.
In Proceedings of the Fourth International Conference on Informatics in Control, Automation and Robotics, pages 207-211
DOI: 10.5220/0001641402070211
Copyright
c
 SciTePress

()
(
)
)42(,,,,
ˆ
,,
ˆ
−= "

tt
ddiiiidii
qqqqQqqqM
Thereinto:
() []
() []
γCMQCCγCMQQ
CCMMCCCCMMM
qFqqqF
qqqqqq
1*11*
111
ˆ
ˆ
−−−
−−−
−−−=
−−−=
ddid
idddid
ddd
T
T
idii
dddi
T
T
idiii
Equation（2-4） is a pure differential formation.
It is suitable for sub-cycling (Lu, 1996).
2.2 Newmark Integration for FMD
The Newmark integration is as following. There, q
t
denotes value of the general variable at time t. Δ t is
the step size.
β and γ are adjustable parameters.
[]
[]
2
(1 2 ) 2 (2 5)
2
(1 ) (2 6)
tt t t t tt
tt t t tt
t
t
t
ββ
γγ
+Δ +Δ
+Δ +Δ
Δ
=+Δ+ − + −
=+Δ − + −
qqq q q
qq qq

"
 
""""""
Define symbols below.
()
tata
t
aa
a
t
a
t
a
t
a
Δ=−Δ=
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
Δ
=−=
−=
Δ
=
Δ
=
Δ
=
γγ
β
γ
β
γ
βββ
γ
β
7654
321
2
0
12
2
1
1
2
111
The dynamic equation at time t+ Δ t can be
established below.
()
(
)
[]
()
)72(,,,
ˆ
,
ˆ
)1()1()1()1(
320
−=
−−−
Δ+Δ+Δ+Δ+
Δ+Δ+Δ+
"


tt
d
tt
d
tt
i
tt
ii
t
i
t
i
t
i
tt
i
tt
d
tt
ii
aaa
qqqqQ
qqqqqqM
Finally, we can get the results.
[]
)82(
ˆ
ˆ
0
3
0
2
1
0
−+++=
Δ+
−
Δ+Δ+
"

t
i
t
i
t
i
tt
i
tt
i
tt
i
a
a
a
a
a qqqQMq
Equation (2-8) need be solved iteratively. The
iteration process is as following.
1
() ( 1) ( 1)
0
3
2
00
ˆ
ˆ
(2 9)
ktt k tt k tt
iiii
ttt
iii
a
a
a
aa
−
+Δ − +Δ − +Δ
⎡⎤
=
⎣⎦
++ + −
qMQ
qqq

"
Thereinto, the top left corner marks represent the
number of the iteration step. Start-up initialisation of
iteration can be set-up below.
[]
)102(2)21(
2
)1(
2
)1(
−+−
Δ
+Δ+=
Δ+Δ+
"

ttttttt
t
t qqqqq
ββ
2.3 Newmark Sub-cycling for FMD
For simplification, all variables are separated into
two categories. The smaller step size is set to be Δ t.
The larger step size is set to be nΔ t. n is a positive
integer. Thus,
q

is expressed as a decomposed
format. And the condensed FMD formula can be
decomposed as a block matrix format as following.
)122(
ˆ
ˆ
ˆˆ
ˆˆ
−
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
"


t
S
t
L
S
L
t
SS
t
SL
t
LS
t
LL
Q
Q
q
q
MM
MM
The subscript symbol, L and S, represent the
larger step size and the smaller step size. M and Q
are the general mass matrix and the general external
force. According to balance status of the FMS at t+
Δ t and t+nΔ t, two groups of equations can be
obtained as following.
11 11 1
ˆ
ˆˆ
(2 13)
ˆ
ˆˆ
(2 14)
LLn Ln LSn Sn Ln
SL L SS S S
+= −
+= −
Mq Mq Q
Mq Mq Q
 
"
 
"
Thereinto,
Ln
q

and
Sn
q

are accelerations at
tnt
Δ
+
.
1L
q

and
1S
q

are accelerations at tt
Δ
+
.
The general mass matrix and the general external
force are defined similarly. Define symbols below.
1
333
1
3
1
2
22
1
2
2
1
0
00
1
0
,,,,, aaaa
n
a
aaa
n
a
aaa
nnn
======
In order to compute interaction between the
coupling variables at the common update, q
Sn
and
q
L1
need be estimated simultaneously. Similar to the
method (Daniel, 1997), q
L1
can be linearly
interpolated and q
Sn
can be linearly extrapolated by
means of the trapezoid rule. The formats of the
interpolations are expressed below.
()
10
00
11
(2 15)
(2 16)
LLLn
Sn Sn S S
n
nn
n
−
=+ −
=−+ −
qqq
qqqq
""""
"""
ICINCO 2007 - International Conference on Informatics in Control, Automation and Robotics
208

Hence action of the slowly changing variables to
the rapidly changing variables is calculated below.
() )172(
ˆ
0
1
30
1
20
1
0
11
−
⎥
⎦
⎤
⎢
⎣
⎡
−−−=− "

LLLLnSLSL
aa
n
a
qqqqMF
The action of the rapidly changing variables to
the slowly changing variables is calculated below.
() []
)182(
ˆ
0001
320
−−−−=− "

S
n
S
n
SS
n
LSnLSn
aaqna qqqMF
Imposing equations (2-17) and (2-18) into
equations (2-13) and (2-14), and co-multiplying a
number, n, to the left hand and the right hand of
equation (2-13), the common update formula of the
slowly changing variables and the rapidly changing
variables can be obtained below.
)192(
ˆˆ
ˆˆ
ˆˆ
ˆˆ
ˆˆ
ˆˆ
ˆ
ˆ
ˆˆ
ˆˆ
0
011
3
0
011
2
0
0
11
0
1
1
11
0
−
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
=
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
"




L
S
LLnLSn
SlSS
L
S
LLnLSn
SlSS
L
S
LLnLSn
SlSS
Ln
S
Ln
S
LLnLSn
SlSS
nn
aa
n
a
n
n
a
q
q
MM
MM
q
q
MM
MM
q
q
MM
MM
Q
Q
q
q
MM
MM
In order to get the sub-step update formula, the
following equation need be calculated at t+（i+1）
Δ t，i=1,2…(n-2).
)202(
ˆ
ˆˆ
)1()1()1()1()1(
−=+
+++++
"

iSiSiSSiLiSL
QqMqM
Also, according to equation (2-7), the following
equation can be obtained.
(
)
()
111
(1) 0 (1) 2 3
111
(1) 0 (1) 2 3 (1)
ˆ
ˆ
ˆ
(2 21)
SL i L i Li Li Li
SSi SiSi SiSiSi
aaa
aaa
++
++ +
⎡⎤
−− − +
⎣⎦
⎡⎤
−− − = −
⎣⎦
Mqqqq
MqqqqQ


"
In equation (2-21), the slowly changing variables,
which are used to compute the interaction of the
coupling variables, can be linearly interpolated in
terms of the trapezoid rule.
)232(
)222(
0
0
−+
−
=
−+
−
=
"""""

"""""
LnLL
LnLLi
n
i
n
in
n
i
n
in
i
qqq
qqq
In terms of equations (2-22) and (2-23), the
action of the slowly changing variables to the
rapidly changing variables can be approximately
assessed.
()
[
()()
()
()
()
(1) 0 0 2 0
1
1
30
1
ˆ
(
(2 24)
SL i Ln L L Ln
SL i
LLn
aanii
n
ani i
+
+
−= −−−+
⎤
−−+ −
⎦
FM qq qq
qq

 
"""" """
Imposing equation (2-24) into equation (2-21),
the sub-step update format can be got as following.
(
)
() ()
()
()
()
0(1)(1) (1) (1)0 2 3
0
2
(1) 0 (1) 0
1
3
(1) 0
ˆ
ˆˆ
ˆˆ
ˆ
(2 25)
SS i S i S i SS i Si Si Si
SL i Ln L SL i L Ln
SL i L Ln
aaaa
a
a
ni i
nn
a
ni i
n
++ + +
++
+
=+ ++
−−+−+
+−+ −
Mq Q M qqq
Mqq M qq
Mqq


 
""""" "
The energy balance status can be calculated by
means of the equation below (Mark and Belytschko,
1989).
)262(
int
−≤−+ ""WWWW
ext
nn
kin
n
δ
Thereinto, W
n
ext
is the work of the external forces
at nΔ t. W
n
int
is the internal energy at nΔ t. W
n
kin
is
the kinetic energy of the system at nΔ t.
δ is the
available error coefficient.
3 NUMERICAL EXAMPLES
In this section, two numerical examples will be
performed to validate availability and efficiency of
the present sub-cycling algorithm.
3.1 A Bar-slider System
A bar-slider system is shown in figure 2. Mass of the
rod is 2.0 kilograms and mass of the slider is 5.0
kilograms. The driving torque is 100 Nm/s. Stiffness
of the spring is 1000 N/m. length of the rigid rod is 2
meters.
Results of rotational angle of the rod, the
vibration amplification of the sliding block and the
energy balance status computed by means of the
sub-cycling and without sub-cycling are shown in
figure 3 to figure 6. The scale values in brackets of
the figure captions, such as (5:1), represent a sub-
cycling with 5 sub-steps in one major step.
A NEWMARK FMD SUB-CYCING ALGORITHM
209

Figure 2: Rotational rod-spring-slider system.
Comparing the results in figure 3 to figure 6, we can
see that no matter which scale of the sub-cycling is
adopted, the results are very similar. The error of
sub-cycling with scale 10:1 is a little larger than that
of sub-cycling with scale 5:1. Yet all these two
errors are very small if compare with the original
results. The time cost of sub-cycling with scale 5:1
is 118 seconds and the time cost of the original
algorithm without sub-cycling is 218 seconds. The
proportion of time cost of these two algorithms is
54%. The time cost of sub cycling with scale 10:1 is
44 seconds and the proportion of time cost is only
20.2%. The results of balance checking illustrate that
the sub-cycling is stable during the integral process.
0 1 2 3 4 5 6 7 8 9 10
0
20
40
60
80
100
120
140
Time(s)
Rotational angle(radian)
no subcycling
subcycling
Figure 3: Rotational angle of the bar (5:1).
0 1 2 3 4 5 6 7 8 9 10
0
20
40
60
80
100
120
140
Time(s)
Rotational angle(radian)
no subcycling
subcycling
Figure 4: Rotational angle of the bar (10:1).
0 1 2 3 4 5 6 7 8 9 10
0
1
2
3
4
Time(s)
Amplification(m)
no subcycling
subcycling
Figure 5: Vibration amplification of the slider (5:1).
0 1 2 3 4 5 6 7 8 9 10
0
1
2
3
4
Time(s)
Amplification(m)
no subcycling
subcycling
Figure 6: Vibration amplification of the slider (10:1).
3.2 Airscrew of a Jet Engine
FEM model of the airscrew of a jet engine are
displayed in figure 7. Parameters of the model are as
following (Units: kg/N/mm/s). The material of
airscrew is aluminium alloy. EX=7e6, PR=3,
DEN=2.78e-6. Diameter of the airscrew is D=900
mm, rotate speed is 8000 rpm.
We simulate the large range overall motion of
the airscrew by means of the Newmark sub-cycling
and the original Newmark respectively. The
dynamic stress at the blade root is described in
figure 8. The time cost of the Newmark sub-cycling
is 1660 seconds and that of the Newmark is 2810
seconds. The computational efficiency is enhanced
about 70%. The compared results show the good
precision and stability of the Newmark sub-cycling.
Figure 7: The FEM model and the local mesh.
ICINCO 2007 - International Conference on Informatics in Control, Automation and Robotics
210

0 0.005 0.01 0.015 0.02
40
50
60
70
80
Time(s)
Dynamic stress(MPa)
no subcycling
subcycling
Figure 8: Dynamic stress of the blade root during rotation.
0 0.005 0.01 0.015 0.02
0
2
4
6
8
10
Time(s)
Rotational angle(radian)
no subcycling
subcycling
Figure 9: Rotational angle of the blade root.
4 CONCLUSIONS
This paper firstly presents a Newmark-based FMD
sub-cycling algorithm. By modifying the Newmark
integral formula to be fitted for sub-cycling of the
FMD problems, not only the integral efficiency can
be greatly improved, but also be more easy for
convergence of the integral process.
Because of that different integral step sizes are
adopted during the sub-cycling, the integral process
can be more efficient and easier for convergence. At
the same time, Unconditional stability of the original
Newmark are still kept in the Newmark sub-cycling.
The number of the sub-step in one major cycle
can be a little infection of the integral precision of a
sub-cycling process. However, it is unobvious as the
number is within a range. Generally speaking, the
enhancement of the integral efficiency is more
significant when the number is under a limitation.
By checking the energy balance status of the
integral process real time and adjusting the step size
when necessary, the sub-cycling procedure can keep
a well convergence property and obtain the
reasonable numerical computation results.
REFERENCES
Huang Wenhu, Shao Chengxun, 1996. Flexible multi-body
dynamics, The science Press.
Dan Negrut, Jose L, Ortiz., 2005. On an approach for the
linearization of the differential algebraic equations of
multi-body dynamics. Proceedings of IDETC/MESA.
DETC2005-85109. 24-28.Long Beach, USA.
Dan Negrut, Gisli Ottarsson, 2005. On an Implementation
of the Hilber-Hughes-Taylor Method in the Context of
Index 3 Differential-Algebraic Equations of Multibody
Dynamics. DETC2005-85096.
Belytschko T., Mullen, Robert, 1978. Stability of explicit-
implicit mesh partitions in time integration.
International Journal for Numerical Methods in
Engineering, 12(10): 1575-1586.
Neal, Mark O., Belytschko T., 1989. Explicit-explicit sub-
cycling with non-integer time step ratios for structural
dynamic systems. Computers and Structures, 31(6):
871-880.
Gao H.,Li G. Y.,Zhong Z. H., 2005. Analysis of sub-
cycling algorithms for computer simulation of
crashworthiness. Chinese journal of mechanical
engineering, 41(11): 98-101.
Tamer M Wasfy, Ahmed K Noor, 2003. Computational
strategies for flexible multi-body systems. Applied
Mechanics Reviews. 56(6): 553-613.
W. J. T. Daniel, 1997. Analysis and implementation of a
new constant acceleration sub-cycling algorithm.
international journal for numerical methods in
engineering. 40: 2841-2855.
Haug, E. J., 1989. Computer-Aided Kinematics and
Dynamics of Mechanical Systems. Prentice-Hall.
Englewood Cliffs, NJ.
Lu Youfang, 1996. Flexible multi-body dynamics. The
higher education Press.
A NEWMARK FMD SUB-CYCING ALGORITHM
211