An Accurate, Controllably Dissipative, Unconditionally Stable Three-Sub-Step Method For Nonlinear Dynamic Analysis of Structures

See discussions, stats, and author profiles for this publication at: https://www.researchgate.
net/publication/366516894
An Accurate, Controllably Dissipative, Unconditionally Stable Three-Sub-

Step Method for Nonlinear Dynamic Analysis of Structures
Article in International Journal of Structural Stability and Dynamics · December 2022

DOI: 10.1142/S0219455423501389
CITATION READS
1 75
3 authors, including:
Yi Ji Yang Wu
Harbin Institute of Technology Institute of Applied Physics and Computational Mathematics
18 PUBLICATIONS 159 CITATIONS 12 PUBLICATIONS 105 CITATIONS
SEE PROFILE SEE PROFILE
All content following this page was uploaded by Yi Ji on 09 March 2023.
The user has requested enhancement of the downloaded file.

International Journal of Structural Stability and Dynamics
(2023) 2350138 (34 pages)
#.c World Scienti¯c Publishing Company
DOI: 10.1142/S0219455423501389
An Accurate, Controllably Dissipative, Unconditionally

by 114.246.204.185 on 03/01/23. Re-use and distribution is strictly not permitted, except for Open Access articles.
Stable Three-Sub-Step Method for Nonlinear Dynamic

Analysis of Structures
Yi Ji*, Yang Wu†,‡,¶ and Yufeng Xing§

Int. J. Str. Stab. Dyn. Downloaded from www.worldscientific.com
*MOE Key Laboratory of Dynamics and Control of Flight Vehicle

School of Aerospace Engineering
Beijing Institute of Technology, Beijing 100081, P. R. China
†
CAEP Software Center for High Performance Numerical Simulation
Beijing 100088, P. R. China
‡
Institute of Applied Physics and Computational Mathematics
Beijing 100088, P. R. China
§Institute of Solid Mechanics
Beihang University, Beijing 100083, P. R. China

¶
wuyangrjzx@163.com
Received 22 October 2022

Accepted 17 December 2022
Published 9 February 2023
An implicit truly self-starting time integration method for nonlinear structural dynamical
systems is developed in this paper. The proposed method possesses unconditional stability,
second-order accuracy, and controllable dissipation, and it has no overshoots. The well-known
BN-stability theory is employed in the design of algorithmic parameters, ensuring that the
proposed method can stably solve nonlinear structural dynamical systems without restricting
the time step size. The spectral analysis shows that compared to existing second-order accurate
time integration methods, the proposed method enjoys a considerable advantage in low-fre-
quency accuracy. For nonlinear problems where the currently popular Generalized- method
and 1 -Bathe method fail, the proposed method shows strong stability and accuracy. Further,
for nonlinear problems in which all methods' results are convergent, the proposed method has
greater accuracy, e±ciency, and energy-conservation capability.
Keywords: Nonlinear structural dynamics; unconditional stability; second-order accuracy;

controllable dissipation; zero-order overshoot.
1. Introduction
Dynamical problems can be solved most e®ectively with time integration methods
(TIMs), which are divided into explicit and implicit methods.1 The computational
e±ciency of explicit methods is higher, but the time step size is limited by their
¶ Corresponding author.
2350138-1
Y. Ji, Y. Wu & Y. Xing
intrinsic conditional stability. In contrast, most implicit methods for linear systems
are unconditionally stable, hence they are more stable for solving nonlinear
dynamical systems.
The Newmark method is a widely-used TIM, which includes some well-known
schemes, such as the Trapezoidal Rule (TR) and the central di®erence method
(CDM). For linear dynamical systems, the implicit TR is unconditionally stable and
strictly energy-conserving, but for simple nonlinear problems such as the softening
spring equation,2 the TR fails. Due to this phenomenon, in recent decades, TIMs
have emerged constantly with better stability.
In addressing nonlinear dynamical systems, the introduction of numerical dissi-
pation is a straightforward method for improving the stability of TIMs. The New-
mark method containing numerical dissipation is only ¯rst-order accurate,
motivating the development of the -type methods.3–7 With the addition of extra
parameters to the equilibrium equation, the -type methods3–6 can obtain second-
order accuracy as well as numerical dissipation. However, due to the discontinuous
external force, the acceleration calculated by the -type methods based on weighted
conception is only ¯rst-order accurate.7 To obtain dissipative TIMs with higher low-
frequency accuracy, linear multistep methods8–13 and composite methods14–21 have
been developed. In the linear multistep methods,8–13 the states of the current time
step are expressed by those of several previous time steps, resulting in the inability to
automatically start the computational process. The composite methods14–21 com-
bined with the TR, and the backward di®erence formula (BDF) or backward in-
terpolation formula (BIF) are self-starting, and since 2005,21 composite methods
have attracted wide attention. The dissipative TIMs,3–23 including -type methods,
linear multistep methods, and composite methods, are more stable when applied to
nonlinear systems, but their accuracy is lower compared to the non-dissipative
methods. Further, numerical dissipation does not prevent divergence but only delays
it, and it cannot guarantee that in simulations, the dissipative TIMs are always
stable.
To provide stable numerical solutions for nonlinear dynamical systems, energy-
conserving methods24–38 based on the Belytschko principle24 have been developed.
Owing to desirable stability for nonlinear problems, this type of method has been
widely employed in structural dynamics,24–33 and has been extended to multibody
dynamics.34,35 As applied to general nonlinear dynamical systems, existing energy-
conserving methods may need additional information, increasing the computations.
Based on the Lagrange multiplier method, the constrain energy method (CEM)36
and the constraint energy momentum algorithm (CEMA)37 were developed, but the
development of this kind of method is limited by the extra computation of discrete
energy and Lagrange multiplier. Simo et al. constructed the energy-momentum
method (EMM)38 based on the modi¯cation of the mid-point rule, and the form
of internal force in the EMM needs to be reevaluated by the modi¯ed algorithmic
strain at the element level, increasing the computational cost. In addition, most
2350138-2
An Accurate, Controllably Dissipative, Unconditionally Stable Three-Sub-Step Method
energy-conserving methods, such as the Krenk method30 and the Orden method,35
can stably solve dynamical systems with geometric nonlinearity, and for damping
nonlinear systems, few energy-conserving methods26 are applicable. It can be con-
cluded that energy-conserving methods do not lose stability when applied to solve
nonlinear problems, but the e±ciency and feasibility of this kind of method need to
be further improved and developed.
To improve the stability, accuracy, and e±ciency of TIMs in solving nonlinear

dynamical systems, an advanced unconditionally stable two-sub-step method39 was
developed recently by the authors of this work. The essential idea of the two-sub-step
method is based on the theory of BN-stability proposed by Butcher in 1975.40 The
BN-stability40 characterizes the stability of a TIM in coping with nonlinear, non-
autonomous problems. It claims that in the process of the step advance of a TIM, if
the norm of the state vector yn is less than or equal to the norm of yn1 for any time
step number n, the TIM is proven to be unconditionally stable (also called BN-
stable40). In other words, the BN-stability theory in fact reveals a condition that is
su±cient to construct an unconditionally stable TIM for nonlinear, non-autonomous
problems. Based on the theory, the developed BN-stable two-sub-step method39
exhibits superior stability and accuracy in solving nonlinear problems in which
the widely used TR and the 1 -Bathe method19 fail. Owing to their remarkable
performance for nonlinear systems, more accurate BN-stable TIMs deserve further
investigation.
From the above review, one can conclude that the energy-conserving methods and
the BN-stable methods solve the issue that TIMs may lose stability for solving
nonlinear systems, and dissipative methods delay the moment in which TIMs lose
stability. Additionally, it can be found that for general nonlinear dynamics, only BN-
stable methods can obtain unconditional stability without accuracy and e±ciency
loss. However, the two-sub-step method with BN-stability39 does not exhibit
advantages in accuracy and e±ciency compared to the widely-used second-order
TIMs, such as the TR and 1 -Bathe method.19
In this context, for nonlinear structural dynamical systems, we develop a new
TIM consisting of three sub-steps, which have BN-stability, second-order accuracy,
controllable dissipation, zero-order overshoots, and truly self-starting capability. The
newly developed three-sub-step method with BN-stability can stably address non-
linear structural dynamical problems, and it has higher accuracy, higher e±ciency,
and stronger energy-conserving capability than the existing two-sub-step scheme
with BN-stability.39
The main contents of this work are organized as follows. The construction process
of the proposed BN-stable three-sub-step method is presented in Sec. 2. Then, the
spectral radius, low-frequency accuracy, overshoot characteristic, and convergence
rate of the proposed method are analyzed in Sec. 3. Some numerical examples that
can con¯rm the superiorities of the proposed method in solving nonlinear dynamical
systems are carried out in Sec. 4. Lastly, Sec. 5 presents the conclusions of this work.
2350138-3
2. Formulation
The need to develop a more accurate BN-stable TIM motivates this work. For
this reason, we construct a BN-stable TIM that consists of three sub-steps.
For convenience, in the rest of this paper, the BN-stable method with multiple
sub-steps and controllable dissipation is referred to as the 1 -multi-sub-step-BN-
stable-n method (1 -MSSBNnÞ, in which 1 and n stand for the amount of high-
frequency dissipation and number of sub-steps, respectively. Therefore, the

existing two-sub-step39 and the newly developed three-sub-step BN-stable methods
are, respectively, named as 1 -MSSBN2 and 1 -MSSBN3. In the following,
the updating formulations and algorithmic parameters of the 1 -MSSBN3 method
are given.
2.1. Updating equations

The general structural dynamical systems are governed by
:: :
Mx þ F ðx; x; tÞ ¼ 0 ð1aÞ
with the given initial conditions

: :
xð0Þ ¼ x0 ; xð0Þ ¼ x0 ; ð1bÞ
where M is the mass matrix, F is the vector containing nonlinear internal force and
external load, x is the displacement vector, t is the time, and the superposed dot
indicates di®erentiation with respect to time. For linear systems, F ðx; x; _ tÞ ¼
Kx þ C x_ QðtÞ where K is the sti®ness matrix, C is the damping matrix, and QðtÞ
is the external force vector.
The three-sub-step 1 -MSSBN3 method is used to solve Eq. (1), in which each
time element [t, t þ t] is classi¯ed into three intervals, [t, t þ c1 t], [t þ c1 t,
t þ c2 t] and [t þ c2 t, t þ c3 t] (0 < c1 < c2 < c3 < 1). The updating equations of
the 1 -MSSBN3 method at t þ c1 t have the forms
8 :: :
< Mxtþc1 t þ F ðxtþc1 t ; xtþc1 t ; t þ c1 tÞ ¼ 0;
>
:
xtþc1 t ¼ xt þ c1 txtþc1 t ; ð2Þ
>
: x: : ::
tþc1 t ¼ x t þ c1 tx tþc1 t :
::
The vectors xtþc1 t , x_ tþc1 t , and xtþc1 t at t þ c1 t can be solved from Eq. (2),
and they are known information in the calculations of the second-sub-step [t þ c1 t,
t þ c2 t]. The state vectors at t þ c2 t are updated by
8 :: :
>
: :
xtþc2 t ¼ xt þ c2 t½ð1 Þxtþc1 t þ xtþc2 t ; ð3Þ
>
: x: : :: ::
tþc2 t ¼ x t þ c2 t½ð1 Þx tþc1 t þ x tþc2 t :
2350138-4
Fig. 1. The diagram of the 1 -MSSBN3 method.
::
from which we can solve for xtþc2 t , x_ tþc2 t , and xtþc2 t . The updates for the last
sub-step [t þ c2 t, t þ c3 t] are given by

8 :: :
>
: : :
xtþc3 t ¼ xt þ c3 t½ð1 Þxtþc1 t þ xtþc2 t þ xtþc3 t ; ð4Þ
>
:x : : :: :: ::
tþc3 t ¼ xt þ c3 t½ð1 Þxtþc1 t þ xtþc2 t þ xtþc3 t :
Once the displacements, velocities, and accelerations at t þ c1 t, t þ c2 t, and

t þ c3 t are solved from Eqs. (2)–(4), the state vectors at the end of the time step are
obtained from
8 : : :
< xtþt ¼ xt þ t½b1 xtþc1 t þ b2 xtþc2 t þ ð1 b1 b2 Þxtþc3 t ;
>
: : :: :: ::
xtþt ¼ xt þ t½b1 xtþc1 t þ b2 xtþc2 t þ ð1 b1 b2 Þxtþc3 t ; ð5Þ
>
: :: :
xtþt ¼ M 1 F ðxtþt ; xtþt ; t þ tÞ:
The solution process of the 1 -MSSBN3 method is described in Fig. 1.

Here we provide some remarks for the 1 -MSSBN3 method's formulas. The initial
::
acceleration vector x0 for most TIMs is necessary, such as the TR and the Gener-
alized- method, which can be solved from the equilibrium equation at the initial
time. However, one can observe from Eq. (2) that in the computation of xtþc1 t ,
: ::
xtþc1 t , and xtþc1 t , the initial acceleration information is not needed, which illus-
trates the 1 -MSSBN3 method being truly self-starting. Additionally, the compu-
tation of Eq. (5) requires only pure vector addition (or subtraction) operations,
meaning that the computations of the 1 -MSSBN3 method are the same as those of
other three-sub-step methods, such as the Wen method16 and the 1 -OTTBIF.18
2.2. Determination of parameters

It can be seen in Eqs. (2)–(5) that the 1 -MSSBN3 method includes eight free
algorithmic parameters, c1 , c2 , c3 , , , , b1 , and b2 . In this section, these free
parameters are designed for second-order accuracy, controllable dissipation, BN-
stability, and the same e®ective sti®ness matrix.
2350138-5
Since the di®erence formulas of displacement and velocity are the same in the
1 -MSSBN3 method, Eqs. (2)–(5) can be reformulated as
:
ztþc1 t ¼ zt þ c1 tztþc1 t ; ð6Þ
: :
ztþc2 t ¼ zt þ c2 t½ð1 Þztþc1 t þ ztþc2 t ; ð7Þ
: : :
ztþc3 t ¼ zt þ c3 t½ð1 Þztþc1 t þ ztþc2 t þ ztþc3 t ; ð8Þ
: : :
ztþt ¼ zt þ t½b1 ztþc1 t þ b2 ztþc2 t þ ð1 b1 b2 Þztþc3 t ; ð9Þ
where z T ¼ ½x T , x_ T . One can ¯nd from Eqs. (6)–(9) that the 1 -MSSBN3 method
can be regarded as a special case of the multi-stage Runge–Kutta methods,40
and hence, the 1 -MSSBN3 method's parameters can be written into Butcher's
tableau,40 as follows
(10)
To reduce the number of free parameters involved in the 1 -MSSBN3 method,

some assumptions are arti¯cially introduced here. First, we require that the diagonal
elements of matrix A given in Eq. (10) are the same, leading to
c1 ¼ c2 ¼ c3 ; ð11Þ
which can ensure that all sub-steps for linear structural dynamical systems share the
identical e®ective sti®ness matrix, referring to the Appendix. For linear structural
dynamical systems, the computational cost of the 1 -MSSBN3 method satisfying the
relations given in Eq. (11) is the same as that of the single-step time integration method,
such as the TR and the Generalized- method. In addition, the time discretization point
t þ c1 t is assumed to be symmetrical with the time discretization point t þ c3 t
about the midpoint of the time step [t, t þ t], yielding the following relation
c1 þ c3 ¼ 1: ð12Þ
It follows that in the last formulation (5), the weights corresponding to the time
discretization points t þ c1 t and t þ c3 t are required to be the same, yielding
b1 ¼ 1 b1 b2 : ð13Þ
Then, together with the four relationships provided in Eqs. (11)–(13), the number
of free parameters involved in the 1 -MSSBN3 method changes from eight to four.
The remaining relations between free parameters are established for second-order
accuracy, tunable dissipation, and BN-stability.
A standard single-degree-of-freedom (SDOF) test equation1 is introduced, and it
is written as
:
z z ¼ 0; ð14Þ
2350138-6
where is the eigenvalue. Applying the formulations (6)–(9) of the 1 -MSSBN3 method
together with Eqs. (11)–(13) to Eq. (14) can yield the following recursive equation
ztþt ¼ AðÞzt ; ¼ t; ð15Þ
where A represents the ampli¯cation factor, and it can be explicitly written as
3 ð2b1 c2 c2 c3 b1 þ b1 c3 þ c2 c3 c 23 þ c 33 b1 c3
2b1 c2 c3 þ b1 c 23 þ b1 c2 c3 þ 1Þ
þ 2 ð3c 23 3c3 þ b1 þ c2 2b1 c2 Þ þ ð3c3 2Þ þ 1
AðÞ ¼ : ð16Þ
3 ðc3 1Þ 3 þ 3 2 ðc3 1Þ 2 þ 3ðc3 1Þ þ 1
To design the convergence order of the 1 -MSSBN3 method, the local truncation
error 1 is considered ¯rst, and its de¯nition is
¼ AðÞ Aexact ðÞ; ð17Þ

where Aexact ðÞ ¼ expðÞ. If a TIM is nth-order accurate, its truncation error satis¯es
¼ Oð nþ1 Þ, which requires that up to nth-derivatives of A at ¼ 0 are all equal to
1, i.e.
Að0Þ ¼ A ð1Þ ð0Þ ¼ ¼ A ðnÞ ð0Þ ¼ 1: ð18Þ
From Eq. (16), the 1 -MSSBN3 method's A(0), A ð1Þ (0) and A ð2Þ (0) can be
derived, and their expressions are
Að0Þ ¼ 1; A ð1Þ ð0Þ ¼ 1; A ð2Þ ð0Þ ¼ 2b1 þ 2c2 4b1 c2 : ð19Þ
It can be observed that the 1 -MSSBN3 method can achieve second-order accuracy
when we require that
A ð2Þ ð0Þ ¼ 2b1 þ 2c2 4b1 c2 ¼ 1: ð20Þ
Equation (20) holds when
1
c2 ¼ : ð21Þ
2
Hence, it is used in the 1 -MSSBN3 method. Substituting Eq. (21) into Eq. (16)
can generate the updated ampli¯cation factor, which has the form
ð2c 33 2c 23 c3 b1 c3 þ 2b1 c 23 þ 1Þ 3
þ ð6c 23 6c3 þ 1Þ 2 þ 2ð3c3 2Þ þ 2
AðÞ ¼ : ð22Þ
2ðc3 þ 1Þ 3
Then, the dissipation behavior of the 1 -MSSBN3 method is designed. The
controllable dissipation ability is expected for damping out the unwanted or spurious
high-frequency information; hence, the ampli¯cation factor at frequency limit !
1 is considered here, which can be expressed as
ð2c 33 2c 23 c3 b1 c3 þ 2b1 c 23 þ 1Þ
AðÞj!1 ¼ : ð23Þ
2ð1 c3 Þ 3
2350138-7
To describe the degree of numerical dissipation, an extra parameter 1 is intro-

duced in the 1 -MSSBN3 method. The value of 1 can range from 1 to 0, which is
given by the users. By solving for AðÞj!1 ¼ 1 , for the case of 0 1 < 1,
ðc3 1Þð21 þ 41 c3 21 c 23 2c 23 þ 1Þ
b1 ¼ ð24aÞ
c3 ð2c3 1Þ
and for the case of 1 ¼ 1,

2c 23 þ 3c3 1
b1 ¼ : ð24bÞ
c3
Last, the 1 -MSSBN3 method' stability characteristics for nonlinear dynamical
systems is designed. If the elements bi > 0 (i ¼ 1; 2; . . . ; n) and the matrix N ¼
diagðbÞA þ A T diagðbÞ bb T is positive semi-de¯nite, the TIM is algebraically sta-

ble.40 The algebraic stability is the su±cient condition of the BN-stability.40 One can
¯nd that the algebraic stability formulates the BN-stability as a series of algebraic
equations; thus, the algebraic stability is used to determine the last two free para-
meters, and c3 . From Butcher's tableau of the 1 -MSSBN3 method given in
Eq. (10), we can obtain that
2 3 2 3
b1 c1
b¼4 b2 5; A ¼ 6 4 c2 ð1 Þ c2
7
5: ð25Þ
1 b1 b2 c3 ð1 Þ c3 c3
Then, in terms of the de¯nition of algebraic stability, the matrix N of the

1 -MSSBN3 method reads
N ¼ diagðbÞA þ A T diagðbÞ bb T
2 32 3
b1 c1
6 76 7
¼4 b2 54 c2 ð1 Þ c2 5
ð1 b1 b2 Þ c3 ð1 Þ c3 c3
2 32 3
c1 c2 ð1 Þ c3 ð1 Þ b1
6 76 7
þ4 c2 c3 54 b2 5
c3 ð1 b1 b2 Þ
2 3
b1
6 7
4 b2 5½b1 b2 ð1 b1 b2 Þ
ð1 b1 b2 Þ
2 3
b1 ð2c1 b1 Þ b2 ðc2 ð1 Þ b1 Þ ð1 b1 b2 Þðc3 ð1 Þ b1 Þ
6 7
6 b2 ðc2 ð1 Þ b1 Þ b2 ð2c2 b2 Þ ð1 b1 b2 Þðc3 b2 Þ 7
¼6
6
7 ð26Þ
7
4 ð1 b1 b2 Þ ð1 b1 b2 Þ ð1 b1 b2 Þ 5
ðc3 ð1 Þ b1 Þ ðc3 b2 Þ ð2c3 ð1 b1 b2 ÞÞ
2350138-8
To ensure that the N is positive semi-de¯nite, we require that

8
> S1 ¼ b1 ð2c1 b1 Þ 0;
>
>
>
> S2 ¼ b2 ð2c2 b2 Þ 0;
>
>
>
>
>
> S3 ¼ ð1 b1 b2 Þ½2c3 ð1 b1 b2 Þ 0;
>
>
>
> S4 ¼ S1 S2 b 22 ½c2 ð1 Þ b1 2 0;
>
>
>
>
>
< S5 ¼ S1 S3 ð1 b1 b2 Þ ½c3 ð1 Þ b1 0;
2 2
S6 ¼ S2 S3 ð1 b1 b2 Þ 2 ðc3 b2 Þ 2 0; ð27Þ
>
> 8 9
>
>
>
>
> < b2 ½c2 ð1 Þ b1 S3
> >
=
>
>
>
> S ¼ S S b ½c ð1 Þ b ðc b Þð1 b b Þ 2
>
>
7 1 6 2 2 1
>
:
3 2 1 2
>
;
>
> ½c ð1 Þ b
>
> 3 1
>
> b ðc b2 Þ½c2 ð1 Þ b1
>
: þ ð1 b1 b2 Þ 2 ½c3 ð1 Þ b1 2 3 0
S2 ½c3 ð1 Þ b1
where S1 , S2 and S3 are the diagonal elements of N; S4 , S5 and S6 represent the
second-order principal minor of N; and S7 is the determination of N. In the
1 -MSSBN3 method, we assume that S7 ¼ 0, leading to
2ð2c 23 3c3 þ 1Þð21 þ 41 c3 21 c 23 2c 23 þ 1Þ
¼ ; 1 2 ½0; 1Þ; ð28aÞ
c3 ð41 6c3 þ 121 c3 121 c 23 þ 41 c 33 þ 4c 33 þ 3Þ
2c3 1
¼ ; 1 ¼ 1: ð28bÞ
2c3
When the above relation (28) is adopted, for the dissipative scheme (0 1 < 1), S4 ,
S5 and S6 given in Eq. (27) of the 1 -MSSBN3 method become
S4 ¼ S5 ¼ S6
2 3
21 ð48c 63 272c 53 þ 640c 43 800c 33 þ 560c 23
6 7
6 208c3 þ 32Þ 1 ð96c 63 þ 320c 53 7
6 7
2ð1 1Þðc3 1Þ 6
4
6 312c 3 64c 3 þ 302c 3 188c3 þ 38Þ 7
4 3 2
7
6 7
4 þ ð48c 63 48c 53 72c 43 þ 96c 33 5
þ 2c 23 36c3 þ 11Þ
¼ ð29Þ
ð2c3 1Þ 6
and for the non-dissipative scheme (1 ¼ 1), S4 , S5 and S6 of the 1 -MSSBN3
method become
ð24c 23 38c3 þ 15Þ 2
S4 ¼ ; S5 ¼ S6 ¼ ð6c 23 11c3 þ 5Þ 2 : ð30Þ
4
To determine the last free parameter c3 , we require that S4 = S5 ¼ S6 ¼ 0,
yielding
c3 ¼ fð1 Þ; 1 2 ½0; 1Þ; ð31aÞ
2350138-9
Fig. 2. The implementation process of the 1 -MSSBN3 method for nonlinear structural dynamical systems.
5
c3 ¼ ; ¼ 1; ð31bÞ
6 1
where fð1 Þ is the function about 1 . All algorithmic parameters of the
1 -MSSBN3 method are clear. The values of Si and bj ði ¼ 1; 2; 3; 4; 5; 6; 7; j ¼ 1; 2)
of the 1 -MSSBN3 method are provided in Table 1. It follows that the values of
Si ði ¼ 1; 2; 3) are more than or equal to zero; the values of Si ði ¼ 4; 5; 6; 7) are equal
Table 1. The values of Si and bj (i ¼ 1; 2; 3; 4; 5; 6; 7; j ¼ 1; 2) of the 1 -MSSBN3 method.
1 S1 S2 S3 S4 S5 S6 S7 b1 b2
0 0.004361 0.017445 0.004361 0 0 0 0 0.348330 0.303340
0.1 0.003818 0.015274 0.003818 0 0 0 0 0.346210 0.307581
0.2 0.003306 0.013225 0.003306 0 0 0 0 0.344288 0.311424
0.3 0.002822 0.011287 0.002822 0 0 0 0 0.342533 0.314933
0.4 0.002361 0.009447 0.002361 0 0 0 0 0.340921 0.318159
0.5 0.001923 0.007693 0.001923 0 0 0 0 0.339430 0.321140
0.6 0.001505 0.006020 0.001505 0 0 0 0 0.338045 0.323910
0.7 0.001105 0.004421 0.001105 0 0 0 0 0.336753 0.326495
0.8 0.000722 0.002888 0.000722 0 0 0 0 0.335543 0.328915
0.9 0.000354 0.001416 0.000354 0 0 0 0 0.334405 0.331190
1 0 0 0 0 0 0 0 1/3 1/3
2350138-10
Table 2. The values of the algorithmic parameters of the 1 -MSSBN3 method.
1 c3 b1 c2 c1 b2
0 0.819575 2ð2c 23 3c3 þ 1Þð21 ðc3 1Þð21 1/2 1 c3 c1 =c2 c1 =c3 1–2 b1
þ 41 c3 21 c 23 þ 41 c3 21 c 23
2c 23 þ 1Þ 2c 23 þ 1Þ
c3 ð41 6c3 c3 ð2c3 1Þ
þ 121 c3 121 c 23
þ 41 c 33 þ 4c 33 þ 3Þ
0.1 0.821381
0.2 0.823054
0.3 0.824614
0.4 0.826076
0.5 0.827452
0.6 0.828751
0.7 0.829983
0.8 0.831153
0.9 0.832268
1 5/6 (2 c3 1Þ=ð2c3 Þ ð2c 23 3c3 þ 1Þ=ðc3 Þ
to zero; and the values of bj ðj ¼ 1,2) are more than zero. Therefore, when relations
(28) and (31) are employed, the 1 -MSSBN3 method is algebraically stable as well as
BN-stable. By using Eqs. (11), (12), (13), (21), (24), (28) and (31), the values of c1 ,
c2 , c3 , b1 , b2 , , , and are completely determined, referring to Table 2.
The designed parameters can ensure that the 1 -MSSBN3 method is second-
order accurate, controllably dissipative, and BN-stable, in which all sub-steps have
the same e®ective sti®ness matrices, refer to Appendix. In Fig. 2, we provide the
implementation process of the 1 -MSSBN3 method for general nonlinear structural
dynamical systems (1).
3. Numerical Properties
The spectral characteristics, overshoot characteristics, and convergence rate of the
1 -MSSBN3 method are analyzed in this section. In the discussion of spectral
characteristics, the 1 -MSSBN3 method's stability behavior, low-frequency accu-
racy, and high-frequency dissipation are deliberately investigated.
3.1. Spectral characteristics

The design of the algorithmic parameters provided in Sec. 2 can ensure that the
1 -MSSBN3 method strictly satis¯es BN-stability. The time integration methods
satisfying BN-stability are unconditionally stable for nonlinear, non-autonomous
dynamic problems, and the related proofs have been provided in detail in Ref. 40.
The current sub-section focuses only on the stability of the 1 -MSSBN3 method in
solving linear systems. It is well known that the numerical dissipation ability in the
high-frequency range and the low-frequency accuracy of a TIM are generally de-
scribed via spectral characteristics.1 Concerning a SDOF system x€ þ ! 2 x ¼ 0 in
2350138-11
which ! represents the natural frequency, the compact recursive scheme of the
1 -MSSBN3 method has the form

xtþt A11 A12 xt xt
: ¼ : ¼A : ð ¼ !tÞ; ð32Þ
txtþt 2 A12 A11 txt txt
where
6 ð2c 63 8c 53 þ 11c 43 4c 33 4c 23 þ 4c3 þ 2b1 c 53 7b1 c 43

þ 9b1 c 33 5b1 c 23 þ b1 c3 1Þ
þ 4 ð6c 43 20c 33 þ 24c 23 12c3 6b1 c 33 þ 9b1 c 23
3b1 c3 þ 2Þ þ 2 ð6c 23 12c3 þ 5Þ þ 2
A11 ¼ 2 2 ð33aÞ
2½ ðc3 1Þ þ 1 3
and
4 ð6c 43 22c 33 þ 30c 23 18c3 þ 6b1 c 43 15b1 c 33 þ 12b1 c 23 3b1 c3 þ 4Þ
þ 2 ð8c 23 14c3 2b1 c 23 þ b1 c3 þ 6Þ þ 2
A12 ¼ 2 2 :
2½ ðc3 1Þ þ 1 3
ð33bÞ
The characteristic equation of the matrix A has the form as

2 A1 þ A2 ¼ 0; ð34Þ
where A1 ¼ trðAÞ and A2 ¼ detðAÞ. The eigenvalues of A determine the spectral
characteristics of a TIM, including the spectral radius, numerical damping ratio, and
period elongation. The de¯nition of the spectral radius is
¼ maxfj1 j; j2 jg; ð35Þ
which performs a TIM's stability behavior for linear systems and high-frequency
dissipation capability. The numerical damping ratio and period elongation (PE)1
can evaluate the amplitude and phase accuracy of TIMs in the low-frequency range,
respectively, and they are de¯ned as
¼ lnðÞ ; ð36Þ

2
T T
PE ¼ ¼ 1; ð37Þ
T
where ¼ arctanðb=aÞ; a and b represent the real and imaginary parts of the
eigenvalues (34), respectively. The spectral radius versus of the 1 -MSSBN3
method is plotted in Fig. 3, and one can see that it is unconditionally stable, satis-
fying 0 ðÞ 1, and the degree of the numerical dissipation can be smoothly
adjusted by 1 . Hence, it can be concluded from Fig. 3 that the 1 -MSSBN3 method
has AN-stability, meaning it is unconditionally stable for linear, non-autonomous
dynamic systems.40 If a time integration method is BN-stable for nonlinear dynamical
2350138-12
Fig. 3. Variation of spectral radii versus of the 1 -MSSBN3 method.
problems, then it is AN-stable for linear ones,40 meaning that time integration
methods adopting the algorithmic parameters designed according to BN-stability
theory are unconditionally stable both for linear and nonlinear dynamical systems.
Since the 1 -MSSBN2 method's spectral characteristics are the same as those of
the currently popular 1 -Bathe method,39 this section compares only the di®erences
between the 1 -MSSBN2 method and the 1 -MSSBN3 method in the low-frequency
accuracy. The two BN-stable methods' amplitude and phase errors versus are
drawn in Fig. 4, respectively, wherein ¼ /n (n is the number of sub-steps) can
ensure that the two methods have the same computations. It follows that the
1 -MSSBN2 method and 1 -MSSBN3 method share the same accuracy perfor-
mances when 1 ¼ 1, while for the dissipative scheme (0 1 < 1), the
1 -MSSBN3 method enjoys a considerable accuracy advantage. Additionally, one
can observe that as 1 increases, the two BN-stable methods' low-frequency accu-
racy, including amplitude and phase, are both simultaneously enhanced.
3.2. Overshoot
In the ¯rst several time steps, the overshooting phenomenon39 may occur for un-
conditionally stable methods. Generally, only the case of ! 1 needs to be con-
sidered for convergent methods, considering that they have no overshoots as ! 0.
For the case of ! 1, the recursive scheme (32) of the 1 -MSSBN3 method at the
¯rst time step becomes

x1 ¼ d 1 x0 ;
: : ð38Þ
tx1 ¼ d2 x0 þ d1 tx0 ;
2350138-13
(a) numerical damping ratio
(b) period elongation
Fig. 4. Amplitude decays and period errors versus of the 1 -MSSBN n method (n ¼ 2: open; n ¼ 3:
solid).
where
2c 63 8c 53 þ 11c 43 4c 33 4c 23 þ 4c3 þ 2b1 c 53 7b1 c 43
þ 9b1 c 33 5b1 c 23 þ b1 c3 1
d1 ¼ ð39aÞ
2ðc3 1Þ 6
and
6c 43 22c 33 þ 30c 23 18c3 þ 6b1 c 43 15b1 c 33 þ 12b1 c 23 3b1 c3 þ 4
d2 ¼ :
2ðc3 1Þ 6
ð39bÞ
2350138-14
Fig. 5. Displacement at the ¯rst time step of the 1 -MSSBN3 method versus t=T .
Fig. 6. Velocity at the ¯rst time step of the 1 -MSSBN3 method versus t=T .
From Eqs. (39a) and (39b), one can ¯nd that when ! 1, the displacement and
velocity both have zero-order overshoots with a large time step size.
To test the overshoot behavior of the 1 -MSSBN3 method, the case x€ þ ! 2 x ¼ 0
(x0 ¼ 1 and x_ 0 ¼ 0) is employed. Figures 5 and 6 plot the displacement and
velocity at the ¯rst time step of the 1 -MSSBN3 method versus t=T , respectively.
It follows that the 1 -MSSBN3 method does not have overshoots both in displace-
ment and velocity, and the displacement and velocity tend to be constant when
t=T ! 1.
2350138-15
Fig. 7. Convergence rates of the 1 -MSSBN3 method (linear case).
3.3. Convergence rate

Convergence rate test is considered here to check the 1 -MSSBN3 method' accuracy
order in displacement, velocity, and acceleration, and its de¯nition is that at a
certain moment, the decrease rate of the absolute error as the step size decreases.
The following linear SDOF equation39 is ¯rst carried out
:: : 57 : 2
x þ 4x þ 5x ¼ sin 2t; xð0Þ ¼ ; xð0Þ ¼ ð40Þ
65 65
and its analytical solution is
1
xðtÞ ¼ expð2tÞðcos t þ 2 sin tÞ ð8 cos 2t sin 2tÞ: ð41Þ
65
Absolute errors at time t ¼ 1 are shown in Fig. 7, and one can observe that for
all state variables, the 1 -MSSBN3 method is strictly second-order accurate.
Additionally, with the increase in 1 , the 1 -MSSBN3 method's accuracy can be
improved.
Then, the accuracy order of the 1 -MSSBN3 method for nonlinear systems is
checked via the Van der Pol equation, i.e.
:: : :
ð1 þ x 2 Þx þ ð2 þ x 2 Þx þ x ¼ 0; xð0Þ ¼ 0:02; xð0Þ ¼ 0 ð42Þ
and the reference solution is obtained by the TR with t ¼ 10 6 . In Fig. 8, the
second-order convergence rate can be observed from the absolute errors of the
1 -MSSBN3 method.
4. Numerical Experiments
The theoretical analysis given in Sec. 3 showed the 1 -MSSBN3 method's perfor-
mances for linear structural dynamical problems, therefore, this section only
considers the 1 -MSSBN3 method's performance in solving nonlinear systems.
2350138-16
Fig. 8. Convergence rates of the 1 -MSSBN3 method (nonlinear case).
To validate the advantages of the 1 -MSSBN3 method in stability, accuracy, and

energy conservation as applied to nonlinear problems, some representative nonlinear
numerical experiments are implemented in this section.
The performances of the 1 -MSSBN3 method are compared with the widely-used
single-step Generalized- method,6 the two-sub-step 1 -Bathe method,19 and the
two-sub-step 1 -MSSBN2 method.39For ensuring that the above mentioned methods
are compared under the roughly same computations in simulations, we require that
t(Generalized-Þ ¼ tð1 -Bathe)/2 = tð1 -MSSBN2)/2 = tð1 -MSSBN3)/3.
Thus, in all examples, we only provide the Generalized- method's size, and the
results obtained from the Generalized- method with a smaller time step size are
used as the reference solution.
4.1. The SDOF equation with the internal force of the fractional function
To perform the advantages of the BN-stable methods in terms of accuracy and
stability, a simple nonlinear SDOF equation26 is ¯rst considered, which has the form
:: x :
xþ ¼ 0; xð0Þ ¼ 0; xð0Þ ¼ 1 ð43Þ
1 þ x2
and its total energy is
1 :2 1
E¼ x þ logð1 þ x 2 Þ: ð44Þ
2 2
With the initial conditions, the system's total energy should hold at 0.5 in theory,
that is, the energy E0 at the initial moment.
For this conservative system, the non-dissipative scheme (1 ¼ 1) is suggested for
the Generalized- method, 1 -Bathe method, 1 -MSSBN2 method, and 1 -MSSBN3
method to avoid energy loss. The results for the case of t(Generalized-Þ ¼ 0.5 are
shown in Fig. 9, in which one can see that all methods give convergent results, but the
2350138-17
Fig. 9. Numerical results in Sec. 4.1 for the case of t (Generalized-Þ ¼ 0:5.
Fig. 10. Numerical results in Sec. 4.1 for the case of t (Generalized-Þ ¼ 2.
two BN-stable methods have higher phase accuracy for all state variables, including
displacement, velocity and acceleration, and smaller energy errors. With the increase
in the time step size, one can ¯nd from Fig. 10 that the Generalized- method and
the 1 -Bathe method both give incorrect predictions, while the results calculated by
the two BN-stable methods are stable. Then, it can be concluded that for the BN-
stable methods, the time step size is selected based on accuracy requirement rather
than stability requirement. For the case of t (Generalized-Þ ¼ 2, the energy of the
Generalized- method, the 1 -Bathe method and the two BN-stable methods
2350138-18
Fig. 11. Energy of the Generalized- method (G-), 1 -Bathe method and the 1 -MSSBN n method for
the case of 1 < 1.
involving numerical dissipation is presented in Fig. 11, from which one can ¯nd that as
the amount of numerical dissipation increases, the 1 -Bathe method and Generalized-
method can provide convergent predictions, and the corresponding energy quickly
declines to zero. Besides, one can observe that the 1 -MSSBN2 method with 1 ¼ 0.99
can keep the energy of the system. Therefore, it can be concluded from Fig. 11 that the
introduction of numerical dissipation can improve the stability of a TIM in dealing with
nonlinear problems, but it cannot ensure globally stable solutions.
4.2. Cantilever beam

The second example considers a cantilever beam, as shown in Fig. 12, in which a
mass is attached at the free end26 to show stability and e±ciency advantages of the
BN-stable methods in solving nonlinear problems. The model is discretized by 10
two-node elements, and its physical parameters are assumed as EA ¼ 70 10 4 N,
EI ¼ 10 4 N m2, A ¼ 0:027 kg/m, L ¼ 4 m, and m ¼ AL=10 ¼ 0:0108 kg, respec-
tively. The model is motivated by an initial vertical velocity v0 .
In the computations, t (Generalized-Þ ¼ 0.002 s, and 1 ¼ 1 is used in all
compared methods. Figures 13 and 14 perform displacements and velocities in the y
direction of the free end with v0 ¼ 20 m/s and 40 m/s, respectively. From Fig. 13, one
can ¯nd that all methods' results are convergent when v0 ¼ 20 m/s. With the in-
crease in v0 , more high-frequency modes are stimulated. Therefore, for the case of
v0 ¼ 40 m/s, the 1 -Bathe method yields divergent results at the early simulation
stage, while the two BN-stable methods still achieve reliable results in the entire
simulation, further supporting their superiority in solving nonlinear dynamical sys-
tems. The e±ciency of these methods is compared in Table 3, and it can be found
2350138-19
Fig. 12. Cantilever beam model.

Fig. 13. Numerical results in Sec. 4.2 for the case of v0 ¼ 20 m/s.
Fig. 14. Numerical results in Sec. 4.2 for the case of v0 ¼ 40 m/s.
2350138-20
Table 3. CPU time and number of iterations in Sec. 4.2 (Tolerance error " ¼ 10e-12).
v0 ¼ 20 m/s, 1 ¼ 1, v0 ¼ 40 m/s, 1 ¼ 1,
Time step total time ¼ 0:5 s total time ¼ 0:5 s
Method size (s) Number of iterations CPU (s) Number of iterations CPU(s)
Generalized- 0.002 24 276 139.1286 25 179 170.6578
1 -Bathe 0.004 24 278 122.3186 — —
1 -MSSBN2 0.004 24 116 114.8487 25083 124.2971

1 -MSSBN3 0.006 23 810 110.1680 25 036 123.6134
that among these methods, the computational cost of the 1 -MSSBN3 method is
the lowest.
From the above two nonlinear numerical experiments, one can ¯nd that compared
with the currently popular Generalized- method and 1 -Bathe method, the
1 -MSSBNn methods (n ¼ 2; 3) have a considerable stability advantage in
solving nonlinear structural dynamical problems. It can be seen that for the
geometric nonlinearity problems, the accuracy of the 1 -MSSBN2 method and the
1 -MSSBN3 method are the same for the non-dissipative scheme (1 ¼ 1). In ad-
dition, the newly developed 1 -MSSBN3 method enjoys a slight e±ciency advantage
compared with the 1 -MSSBN2 method in solving nonlinear problems.
4.3. Spring pendulum

A spring pendulum model39 is considered here to validate the superiority of the
1 -MSSBN3 method for addressing nonlinear damping systems. The diagram of this
mode is described in Fig. 15, and its governing equation is written as
8 ::
> _2
< mr mðL0 þ rÞ mg cos þ kr ¼ 0;
: ð45Þ
> mð2r _ þ g sin Þ
: m€þ ¼ 0:
L0 þ r
Fig. 15. Spring pendulum model.
2350138-21
Table 4. Average absolute errors of these methods in Sec. 4.3.
Generalized- 1 -Bathe 1 -MSSBN2 1 -MSSBN3

Displacement r (1 ¼ 1) 2.5253e-02 2.6534e-02 2.5255e-02 2.4838e-02
(1 ¼ 1) 3.3906e-02 4.2147e-02 3.3934e-02 3.3852e-02
r (1 ¼ 0) 1.1241e-01 5.1329e-02 4.8233e-02 3.3073e-02
(1 ¼ 0) 1.3467e-01 7.8491e-02 6.1289e-02 4.5199e-02
:
¼ 1)
Velocity r (1 2.6727e-01 2.8011e-01 2.6685e-01 2.6686e-01

_ð1 ¼ 1) 2.0770e-01 2.3914e-01 2.0710e-01 2.0244e-01
:
r (1 ¼ 0) 1.1744eþ00 5.4279e-01 5.1293e-01 3.5110e-01
_ð1 ¼ 0) 8.4242e-01 4.3750e-01 3.7911e-01 2.6547e-01
::
Acceleration r (1 ¼ 1) 2.8338eþ00 2.9644eþ00 2.8278eþ00 2.7983eþ00
€ð1 ¼ 1) 2.2643eþ00 2.5731eþ00 2.2571eþ00 2.2758eþ00
::
r (1 ¼ 0) 1.3999eþ01 5.7100eþ00 5.4107eþ00 3.7014eþ00
€ð1 ¼ 0) 7.9549eþ00 4.6810eþ00 4.0919eþ00 2.9935eþ00
In this model, the initial radial and circumferential displacements are assumed to
be rð0Þ ¼ 0:25 m and ð0Þ ¼ 2 /9, and the initial velocities are equal to zero. The
parameters of this model are chosen as m ¼ 1 kg, g ¼ 9:81 m/s2, L0 ¼ 0:5 m, and
k ¼ 98:1 N/m.
The displacements, velocities, and accelerations within [0, 10]s are plotted in Figs.
16–18, in which t(Generalized-Þ ¼ 0.01 s. It follows that all methods have no
considerable numerical errors at the beginning, while with the increase in time their
accuracy decreases. It can also be found that the accuracy of the non-dissipative
scheme (1 ¼ 1) is higher than that of the dissipative scheme (1 ¼ 0). To distin-
guish the di®erences between these methods, the absolute errors of these methods are
provided in Table 4, wherein one can ¯nd that the accuracy of the 1 -MSSBN3
(a) 1 ¼ 1
Fig. 16. Numerical displacements in Sec. 4.3.
2350138-22
(b) 1 ¼ 0
Fig. 16. (Continued )
method is the highest, followed by the 1 -MSSBN2 method. This example illustrates
that the 1 -MSSBN3 method enjoys accuracy advantages in solving nonlinear
systems, including physical damping. Considering that physical damping is helpful
to improve the stabilities of time integration methods, the non-dissipative scheme
(1 ¼ 1) is used to obtain higher accuracy for the dynamical problems containing
nonlinear damping.
(a) 1 ¼ 1
Fig. 17. Numerical velocities in Sec. 4.3.
2350138-23
(b) 1 ¼ 0
4.4. Rigid linkage

This example considers a rigid linkage,15 as shown in Fig. 19, in which the Jacobian
matrix is not positive, hence the dissipative schemes are suggested that can improve the
condition number of Jacobian matrix. The model is motivated by the initial conditions
:
’ð0Þ ¼ =6, xð0Þ ¼ 1, ’ð0Þ ¼ 0, and pffiffiffi xð0Þ
_ ¼ 0. The physical parameters of this model
are assumed as mass m ¼ 1, a ¼ 3, b ¼ c ¼ 1, and gravity acceleration g ¼ 9:81.
(a) 1 ¼ 1
Fig. 18. Numerical accelerations in Sec. 4.3.
2350138-24
(b) 1 ¼ 0
Fig. 19. Rigid linkage.
The results of displacement, velocity, and acceleration of the Generalized-

method, 1 -Bathe method, 1 -MSSBN2 method, and 1 -MSSBN3 method are
shown in Figs. 20–22, in which t(Generalized-Þ ¼ 0.01, and 1 ¼ 0 is employed in
all compared methods. Figures 20–22 show that in addition to the 1 -MSSBN3
method, other methods exhibit obvious phase shifts after a while. In theory, the
mechanical energy of this model should hold at 16.9914 for the initial conditions.
Figure 23 shows the energy of these methods, and one can see that (1) the
1 -MSSBN3 method's energy loss is the lowest; (2) with the increase in 1 , the
energy loss of the Generalized- method and the 1 -MSSBN3 method declines,
2350138-25
Fig. 20. Numerical displacements in Sec. 4.4.
Fig. 21. Numerical velocities in Sec. 4.4.
while other methods' energy loss becomes larger. For this nonlinear system in
which the Jacobian matrix is not positive, the conclusion obtained from spectral
characteristics39 that with the increase in 1 , the amplitude errors get lower does not
hold. Therefore, the 1 -MSSBN2 method and the 1 -Bathe method both perform
considerable energy loss for the case of 1 ¼ 0.8. It can be concluded from this
example that among these dissipative schemes (1 < 1), the 1 -MSSBN3 method
enjoys superiority in accuracy and energy conservation.
2350138-26
Fig. 22. Numerical accelerations in Sec. 4.4.
Fig. 23. Numerical energy in Sec. 4.4.
4.5. Elastoplastic plate

As shown in Fig. 24, a clamped plate with end impact load is considered to further
compare the performances of the BN-stable methods with the Generalized- method
and the 1 -Bathe method in solving complex nonlinear problems. Both geometry
and material nonlinearity are involved in the problem since ¯nite deformation ele-
ments with the isotropic hardening J2 bilinear Hypoelastoplastic material are used in
the model. The material parameters are: Young's modulus E ¼ 200 GPa, Poisson's
2350138-27
(a) model information (b) load information
Fig. 24. Clamped plate model: (a) model information; (b) load information.
ratio ¼ 0:3, the density ¼ 7000 kg/m3, the elastoplastic tangent modulus Et = 50
GPa, and the initial yield stress y0 ¼ 200 MPa. The geometric parameters are as-
sumed to be a ¼ 3 m, b ¼ 1 m, and h ¼ 0:055 m. The impact forces are applied at each
node at the end of the plate, and the impact load curve is described in Fig. 24(b), in
which one can see that the external force rapidly increases from 0 N to the maximum
value of 50 000 N at a very short time of 0.1 s and is removed soon after. This structure
has 1 296 elements, as shown in Fig. 24(a), and 5 772 degrees of freedom.
The deformed con¯gurations of the elastoplastic plate at di®erent moments are
drawn in Fig. 25, in which one can see that the plate produces the irreversible
deformation. In the calculation, the Generalized- method's step size is chosen to be
0.01 s. Figures 26 and 27 show the displacements in the x and z direction of point A
(the central point of the free end), respectively, in which one can ¯nd that at the
Fig. 25. Contours on deformed con¯gurations at di®erent moments in Sec. 4.5.
2350138-28
Fig. 26. Displacements of point A in the x direction in Sec. 4.5 (1 ¼ 1).
Fig. 27. Displacements of point A in the z direction in Sec. 4.5 (1 ¼ 1).
beginning, all methods provide accurate results, while they perform the phase errors
with the increase in time. Additionally, one can ¯nd that after the external force is
removed, the responses tend to be constant (the plastic deformations) due to the
presence of physical damping. It can be concluded from Figs. 26 and 27 that among
these methods, the 1 -MSSBN3 method achieves a considerable accuracy advantage.
2350138-29
Fig. 28. Variation of displacement versus ! 2 =! 20 in Sec. 4.6 (1 ¼ 1).
4.6. Du±ng equation

To analyze the performances of the BN-stable methods in the calculations of fre-
quency and acceleration, the last example considers the following Du±ng equation
::
x þ ð! 20 þ "x 2 Þx ¼ 0; !0 ¼ 1: ð46Þ
In classi¯cation, this model belongs to the sti®ness hardening system when " > 0 and
the sti®ness softening system when " < 0. In solving the problem, the time step size of
the Generalized- method is 0.2. Figure 28 plots the variations of displacement
versus the frequency ratio ! 2 =! 20 , in which the variation frequency ! 2 approximately
Fig. 29. Variation of acceleration versus time in Sec. 4.6 (1 ¼ 1).
2350138-30
equals to (! 20 þ "x 2 Þ.41 It can be seen from Fig. 28 that (1) for the case of " ¼ 0:4, the
system sti®ness increases with the increase of jxj; (2) for the case of " ¼ 0:4, the
system sti®ness decreases with the increase of jxj. For the weak nonlinear sti®ness
problems, we also provide the acceleration results as shown in Fig. 29. One can ¯nd
that their accuracy is very close, and among them, the accuracy of the 1 -Bathe
method is slightly lower.
5. Conclusions
For structural dynamical systems, we construct a second-order accurate and °exibly
dissipative TIM with BN-stability, called the 1 -MSSBN3 method. In addition, the
1 -MSSBN3 method has no overshoots in displacement and velocity. Unlike most
existing self-starting time integration methods, the information of initial acceleration is
not necessary for the 1 -MSSBN3 method, meaning that it is a truly self-starting TIM.
The theoretical analysis given in Sec. 3 con¯rmed that the newly developed
1 -MSSBN3 method achieves higher low-frequency accuracy than the existing
1 -MSSBN2 method.39 In other words, as the number of sub-steps increases, the
accurate low-frequency range can be widened. The superiority of the 1 -MSSBN3
method as applied to nonlinear systems has been validated by nonlinear numerical
examples, from which it can be concluded that (1) the 1 -MSSBNn methods
(n ¼ 2; 3) are unconditionally stable for nonlinear structural dynamical systems, and
they can provide stable predictions for nonlinear problems wherein the 1 -Bathe
method and Generalized- method fail, seeing Secs. 4.1 and 4.2; (2) for nonlinear
systems containing physical damping, the 1 -MSSBN3 method achieves higher ac-
curacy than the currently popular second-order methods, referring to Secs. 4.3
and 4.5; (3) for the case of 0 1 < 1, compared to some advanced methods, the
1 -MSSBN3 method enjoys considerable advantages in low-frequency accuracy and
energy conservation, referring to Secs. 4.3 and 4.4; (4) the 1 -MSSBN3 method has
slight advantage in computational e±ciency, and it has less CPU time and number of
iterations, refer to Sec. 4.2.
By increasing the number of sub-steps, a more desirable BN-stable method con-
sisting of three sub-steps was developed in this work. The theoretical analysis and the
limited numerical experiments illustrate that the 1 -MSSBN3 method is a good
candidate for solving structural dynamical systems, especially for nonlinear pro-
blems. On this basis, it is expected to further increase the computational e±ciency of
the 1 -MSSBN3 method. Therefore, the solution strategy based on this work and
domain decomposition technology for large-scale complex dynamics problems will be
studied in the future.
Acknowledgments
This work was supported by the National Natural Science Foundation of China
(12202058, 12172023, 11872090), the China Postdoctoral Science Foundation
2350138-31
(2022M710386), and the Outstanding Research Project of Shen Yuan Honors

College BUAA (230121104).
Competing Interest
The authors declare that they have no known competing ¯nancial interests or
personal relationships that could have appeared to in°uence the work reported in
this paper.
Appendix A. E®ective Sti®ness Matrices and Load Vectors of the

1 -MSSBN3 Method
The time-stepping equations for linear structural dynamics can be written as
8
> b b
< K 1 xtþc1 t ¼ R1 ;
>
Kb 2 xtþc t ¼ R
b 2; ðA:1Þ
>
>
2
:b b
K 3 xtþc t ¼ R3 ; 3
where the e®ective sti®ness matrices are
b1 ¼ 1 1
K Mþ C þ K; ðA:2Þ
c 21 t 2 c1 t
b2 ¼ 1 1
K Mþ C þ K; ðA:3Þ
c 22 2 t 2 c2 t
b3 ¼ 1 1
K Mþ CþK ðA:4Þ
c 23 2 t 2 c3 t
and the load vectors have the forms as

b 1 ¼ Qðt þ c1 tÞ þ M 1 1 : 1
R x þ x þ Cxt ; ðA:5Þ
c 21 t 2 t c1 t t c1 t

b 1 1 : 1 :
R2 ¼ Qðt þ c2 tÞ þ M 2 2 2 xt þ xt þ ð1 Þxtþc1 t
c 2 t c2 t c2 2 t

1 :: 1 1 :
þ ð1 Þxtþc1 t þ C x þ ð1 Þxtþc1 t ; ðA:6Þ
c2 t t

b 3 ¼ Qðt þ c3 tÞ þ M 1 1 : ð1 Þ :
R xt þ x þ xtþc1 t
c 23 2 t 2 c3 t t c3 2 t

: ð1 Þ :: ::
þ x þ x þ x
c3 2 t tþc2 t tþc1 t
tþc2 t

1 ð1 Þ : :
þC xt þ xtþc1 t þ xtþc2 t : ðA:7Þ
c3 t
2350138-32
References
1. T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element
Analysis (Prentice-Hall, New Jersey, 1987).
2. K. C. Park, An improved sti®ly stable method for direct integration of nonlinear
structural dynamic equations, J. Appl. Mech.-Trans. ASME 42 (1975) 464–470.
3. H. M. Hilber, T. J. R. Hughes and R. L. Taylor, Improved numerical dissipation for
time integration algorithms in structural dynamics, Earthq. Eng. Struct. Dyn. 5 (1977)
283–292.
4. W. L. Wood, M. Bossak and O. C. Zienkiewicz, An alpha modi¯cation of Newmark's
method, Int. J. Numer. Methods Eng. 15 (1980) 1562–1566.
5. H. P. Shao and C. W. Cai, A three parameters algorithm for numerical integration of
structural dynamic equations, Chin. J. Appl. Mech. 5 (1988) 76–81, (In Chinese).
6. J. Chung and G. M. Hulbert, A time integration algorithm for structural dynamics with
improved numerical dissipation: The generalized- method, J. Appl. Mech.-Trans.
ASME 60 (1993) 371–375.

7. H. M. Zhang and Y. F. Xing, A three-parameter single-step time integration method for
structural dynamic analysis, Acta Mech. Sin. 35 (2019) 112–128.
8. X. Zhou and K. K. Tamma, Design, analysis, and synthesis of generalized single step
single solve and optimal algorithms for structural dynamics, Int. J. Numer. Methods Eng.
59 (2004) 597–668.
9. Y. Ji and Y. F. Xing, A two-step time integration method with desirable stability for
nonlinear structural dynamics, Eur. J. Mech. A-Solids 94 (2022) 104582.
10. J. Zhang, A-stable linear two-step time integration methods with consistent starting and
their equivalent single-step methods in structural dynamics analysis, Int. J. Numer.
Methods Eng. 122 (2021) 2312–2359.
11. Y. Ji and Y. F. Xing, A two-sub-step generalized central di®erence method for general
dynamics, Int. J. Struct. Stab. Dyn. 20 (2020) 2050071.
12. H. M. Zhang, R. S. Zhang and P. Masarati, Improved second-order unconditionally stable
schemes of linear multi-step and equivalent single-step integration methods, Comput.
Mech. 67 (2021) 289–313.
13. S. Dong, BDF-like methods for nonlinear dynamic analysis, J. Comput. Phys. 229 (2010)
3019–3045.
14. J. Z. Li, K. P. Yu and T. Hong, Further assessment of three Bathe algorithms and
implementations for wave propagation problems, Int. J. Struct. Stab. Dyn. 21 (2021)
2150073.
15. Y. Ji and Y. F. Xing, Optimization of a class of n-sub-step time integration methods for
structural dynamics, Int. J. Appl. Mech. 13 (2021) 2150064.
16. T. H. Liu, F. L. Huang, W. B. Wen, X. H. He, S. Y. Duan and D. N. Fang, Further
insights of a composite implicit time integration scheme and its performance on linear
seismic response analysis, Eng. Struct. 241 (2021) 112490.
17. Y. F. Xing, Y. Ji and H. M. Zhang, On the construction of a type of composite time
integration methods, Comput. Struct. 221 (2019) 157–178.
18. Y. Ji and Y. F. Xing, An optimized three-sub-step composite time integration method
with controllable numerical dissipation, Comput. Struct. 231 (2020) 106210.
19. G. Noh and K. J. Bathe, The Bathe time integration method with controllable spectral
radius: The 1-Bathe method, Comput. Struct. 212 (2019) 299–310.
20. W. Kim, An improved implicit method with dissipation control capability: The simple
generalized composite time integration algorithm, Appl. Math. Model. 81 (2020) 910–930.
21. K. J. Bathe and M. M. I. Baig, On a composite implicit time integration procedure for
nonlinear dynamics, Comput. Struct. 83 (2005) 2513–2524.
2350138-33
22. M. Rezaiee-Pajand, S. A. H. Esfehani and H. Ehsanmanesh, An e±cient weighted residual

time integration family, Int. J. Struct. Stab. Dyn. 21 (2021) 2150106.
23. S. Y. Chang, A family of matrix coe±cient formulas for solving ordinary di®erential
equations, Appl. Math. Comput. 418 (2022) 126811.
24. T. Belytschko and D. F. Schoeberle, On the unconditional stability of an implicit
algorithm for nonlinear structural dynamics, J. Appl. Mech.-Trans. ASME 42 (1975)
865–869.
25. S. Lopez, A dissipative momentum-conserving time integration algorithm for nonlinear

structural dynamics, Int. J. Struct. Stab. Dyn. 21 (2021) 2150116.
26. H. M. Zhang, Y. F. Xing and Y. Ji, An energy-conserving and decaying time integration
method for general nonlinear dynamics, Int. J. Numer. Methods Eng. 121 (2020) 925–944.
27. D. Kuhl and E. Ramm, Generalized energy-momentum method for nonlinear adaptive
shell dynamics, Comput. Methods Appl. Mech. Eng. 178 (1999) 343–366.
28. B. Wu, T. L. Pan, H. W. Yang, J. Z. Xie and B. F. Spencer, Energy-consistent integration
method and its application to hybrid testing, Earthq. Eng. Struct. Dyn. 49 (2020) 415–433.
29. R. Zhang and H. Z. Zhong, A quadrature element formulation of an energy-momentum
conserving algorithm for dynamic analysis of geometrically exact beams, Comput. Struct.
165 (2016) 96–106.
30. S. Krenk, Global format for energy–momentum based time integration in nonlinear
dynamics, Int. J. Numer. Methods Eng. 100 (2014) 458–476.
31. O. Gonzalez, Exact energy and momentum conserving algorithms for general models in
nonlinear elasticity, Comput. Methods Appl. Mech. Eng. 190 (2000) 1763–1783.
32. S. Mamouri, F. Hammadi and A. Ibrahimbegovic, Decaying/conserving implicit scheme
and non-linear instability analysis of 2D frame structures, Int. J. Non-Linear Mech. 67
(2014) 144–152.
33. H. Jahromi and B. A. Izzuddin, Energy conserving algorithms for dynamic contact
analysis using Newmark methods, Comput. Struct. 118 (2013) 74–89.
34. J. H. Luo, X. G. Feng, X. M. Xu, H. J. Peng and Z. G. Wu, A parameter-preadjusted
energy-conserving integration for rigid body dynamics in terms of convected base vectors,
Int. J. Numer. Methods Eng. 121 (2020) 4921–4943.
35. J. C. Orden, Energy and symmetry-preserving formulation of nonlinear constraints and
potential forces in multibody dynamics, Nonlinear Dyn. 95 (2019) 823–837.
36. T. J. R. Hughes, T. K. Caughey and W. K. Liu, Finite-element methods for nonlinear
elastodynamics with conserve energy, J. Appl. Mech. - Trans. ASME 45 (1978) 366–370.
37. D. Kuhl and E. Ramm, Constraint energy momentum algorithm and its application to
non-linear dynamics of shells, Comput. Methods Appl. Mech. Eng. 136 (1996) 293–315.
38. J. C. Simo and K. K. Wong, Unconditionally stable algorithms for rigid body dynamics
that exactly preserve energy and momentum, Int. J. Numer. Methods. Eng. 31 (1991) 19–52.
39. Y. Ji, Y. F. Xing and M. Wiercigroch, An unconditionally stable time integration method
with controllable dissipation for second-order nonlinear dynamics, Nonlinear Dyn. 105
(2021) 3341–3358.
40. J. C. Butcher, Numerical Methods for Ordinary Di®erential Equations (John Wiley &
Sons, Ltd, Chichester, 2016).
41. S. Y. Chang, An explicit method with improved stability property, Int. J. Numer. Meth.
Engng. 77 (2009) 1100–1120.
2350138-34
View publication stats

An Accurate, Controllably Dissipative, Unconditionally Stable Three-Sub-Step Method For Nonlinear Dynamic Analysis of Structures

Uploaded by

Copyright:

Available Formats

An Accurate, Controllably Dissipative, Unconditionally Stable Three-Sub-Step Method For Nonlinear Dynamic Analysis of Structures

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

An Accurate, Controllably Dissipative, Unconditionally Stable Three-Sub-Step Method For Nonlinear Dynamic Analysis of Structures

Uploaded by

Copyright:

Available Formats

See discussions, stats, and author profiles for this publication at: https://www.researchgate.

An Accurate, Controllably Dissipative, Unconditionally Stable Three-Sub-

Article in International Journal of Structural Stability and Dynamics · December 2022

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Yi Ji on 09 March 2023.

The user has requested enhancement of the downloaded file.

An Accurate, Controllably Dissipative, Unconditionally

Stable Three-Sub-Step Method for Nonlinear Dynamic

Yi Ji*, Yang Wu†,‡,¶ and Yufeng Xing§

*MOE Key Laboratory of Dynamics and Control of Flight Vehicle

Beihang University, Beijing 100083, P. R. China

Received 22 October 2022

Keywords: Nonlinear structural dynamics; unconditional stability; second-order accuracy;

To improve the stability, accuracy, and e±ciency of TIMs in solving nonlinear

frequency dissipation and number of sub-steps, respectively. Therefore, the

2.1. Updating equations

with the given initial conditions

Fig. 1. The diagram of the 1 -MSSBN3 method.

sub-step [t þ c2 t, t þ c3 t] are given by

Once the displacements, velocities, and accelerations at t þ c1 t, t þ c2 t, and

The solution process of the 1 -MSSBN3 method is described in Fig. 1.

2.2. Determination of parameters

To reduce the number of free parameters involved in the 1 -MSSBN3 method,

 ¼ AðÞ Aexact ðÞ; ð17Þ

To describe the degree of numerical dissipation, an extra parameter 1 is intro-

and for the case of 1 ¼ 1,

diagðbÞA þ A T diagðbÞ bb T is positive semi-de¯nite, the TIM is algebraically sta-

Then, in terms of the de¯nition of algebraic stability, the matrix N of the

To ensure that the N is positive semi-de¯nite, we require that

Table 1. The values of Si and bj (i ¼ 1; 2; 3; 4; 5; 6; 7; j ¼ 1; 2) of the 1 -MSSBN3 method.

Table 2. The values of the algorithmic parameters of the 1 -MSSBN3 method.

3.1. Spectral characteristics

6 ð2c 63 8c 53 þ 11c 43 4c 33 4c 23 þ 4c3 þ 2b1 c 53 7b1 c 43

The characteristic equation of the matrix A has the form as

Fig. 3. Variation of spectral radii versus of the 1 -MSSBN3 method.

(a) numerical damping ratio

(b) period elongation

Fig. 7. Convergence rates of the 1 -MSSBN3 method (linear case).

3.3. Convergence rate

Fig. 8. Convergence rates of the 1 -MSSBN3 method (nonlinear case).

To validate the advantages of the 1 -MSSBN3 method in stability, accuracy, and

4.2. Cantilever beam

Fig. 12. Cantilever beam model.

1 -MSSBN2 0.004 24 116 114.8487 25083 124.2971

4.3. Spring pendulum

Fig. 15. Spring pendulum model.

Table 4. Average absolute errors of these methods in Sec. 4.3.

Generalized- 1 -Bathe 1 -MSSBN2 1 -MSSBN3

Velocity r (1 2.6727e-01 2.8011e-01 2.6685e-01 2.6686e-01

Fig. 16. Numerical displacements in Sec. 4.3.

Fig. 16. (Continued )

Fig. 17. Numerical velocities in Sec. 4.3.

Fig. 17. (Continued )

4.4. Rigid linkage

Fig. 18. Numerical accelerations in Sec. 4.3.

Fig. 18. (Continued )

Fig. 19. Rigid linkage.

The results of displacement, velocity, and acceleration of the Generalized-

¼ AðÞ Aexact ðÞ; ð17Þ