This article was downloaded by: [National Institute of Technology - Rourkela]

On: 13 February 2014, At: 09:17

Publisher: Taylor & Francis
Informa Ltd Registered in England and Wales Registered Number: 1072954 Registere
d office: Mortimer
House, 37-41 Mortimer Street, London W1T 3JH, UK
Journal of the Chinese Institute of Engineers
Publication details, including instructions for authors and subscription informa
tion:
http://www.tandfonline.com/loi/tcie20
Studies of Newmark method for solving nonlinear
systems: (I) basic analysis
Shuenn-Yih Chang a
a Department of Civil Engineering , National Taipei University of Technology , T
aipei,
Taiwan 106, R.O. C. Phone: 886227712171 ext. 2653 Fax: 886227712171 ext. 2653 Email:
Published online: 04 Mar 2011.
To cite this article: Shuenn-Yih Chang (2004) Studies of Newmark method for solv
ing nonlinear systems: (I) basic analysis,
Journal of the Chinese Institute of Engineers, 27:5, 651-662, DOI: 10.1080/02533
839.2004.9670913
To link to this article: http://dx.doi.org/10.1080/02533839.2004.9670913
PLEASE SCROLL DOWN FOR ARTICLE
Taylor & Francis makes every effort to ensure the accuracy of all the informatio
n (the Content) contained
in the publications on our platform. However, Taylor & Francis, our agents, and
our licensors make no
representations or warranties whatsoever as to the accuracy, completeness, or su
itability for any purpose of
the Content. Any opinions and views expressed in this publication are the opinio
ns and views of the authors,
and are not the views of or endorsed by Taylor & Francis. The accuracy of the Co
ntent should not be relied
upon and should be independently verified with primary sources of information. T
aylor and Francis shall
not be liable for any losses, actions, claims, proceedings, demands, costs, expe
nses, damages, and other
liabilities whatsoever or howsoever caused arising directly or indirectly in con
nection with, in relation to or
arising out of the use of the Content.
This article may be used for research, teaching, and private study purposes. Any
substantial or systematic
reproduction, redistribution, reselling, loan, sub-licensing, systematic supply,
or distribution in any
form to anyone is expressly forbidden. Terms & Conditions of access and use can
be found at http://
www.tandfonline.com/page/terms-and-conditions
Journal of the Chinese Institute of Engineers, Vol. 27, No. 5, pp. 651-662 (2004
) 651
STUDIES OF NEWMARK METHOD FOR SOLVING NONLINEAR
SYSTEMS: (I) BASIC ANALYSIS
Shuenn-Yih Chang
ABSTRACT
The performance of the Newmark method in the solution of linear elastic systems
can be reliably predicted by the currently available evaluation techniques.
However, its performance for solving nonlinear systems is still not clear. There
fore,
its numerical characteristics in the solution of nonlinear systems are explored
in this
study by using a newly developed technique. Only a specific implementation of th
e
Newmark method is studied in this paper although there are several possible impl
ementations
for this integration method. It seems numerical properties of the Newmark
method in the solution of linear elastic and nonlinear systems can be entirely c
aptured
by the newly developed technique. Although this technique is only aimed at a spe
cific
time step, it is still indicative for the whole step-by-step integration procedu
re
since this procedure consists of each time step.
Key Words: nonlinear system, degree of nonlinearity, evaluation, stability, accu
racy.
S. Y. Chang is with the Department of Civil Engineering, National
Taipei University of Technology, Taipei, Taiwan 106, R.O.
C. (Tel: 886-2-27712171 ext. 2653; Fax: 886-2-27814518; Email:
[email protected])
I. INTRODUCTION
Step-by-step integration methods (Bathe and
Wilson, 1973; Chang, 1996; 1997; 2001a; 2002;
Goudreau and Taylor, 1972; Hilber et al., 1977;
Newmark, 1959) are widely used for solving equations
of motion so that dynamic behaviors of the system
under a specific loading can be obtained. In fact,
it is well recognized that the step-by-step integration
method is the most powerful technique for nonlinear
dynamic analysis. The performance of integration
methods is very important in practical application and
thus they have been thoroughly studied. However,
numerical properties of the integration method in the
solution of equations of motion are usually evaluated
for linear elastic systems only (Belytschko and
Hughes, 1983; Chang, 2001b; 2003a; Goudreau and
Taylor, 1972; Krieg, 1973) and those for nonlinear
systems are not assessed (Belytschko and Schoeberle,
1975). The Newmark method is a widely used stepby-
step integration method in dynamic response
analysis and thus is chosen for this study. Its performance
in the solution of nonlinear systems is assessed
by a newly developed technique (Chang, 2003b). In
fact, its numerical characteristics, such as local truncation
error, order of accuracy, overshooting, stability
and accuracy, in the solution of nonlinear systems
are evaluated by this technique. The details of using
this technique to evaluate the Newmark method are
presented. Meanwhile, analytical results are discussed
and correlated to numerical phenomena
experienced.
II. EQUATIONS OF MOTION
In dynamic analysis, a structural system is often
idealized as a discrete model and thus the equations
of motion for the discrete structural model are
expressed as:
Mu+Cu+Ku=f (1)
where M, C and K are the mass, viscous damping and
stiffness matrices; u, u and u are the displacement,
velocity and acceleration vectors, respectively, and f
is an external force vector. In general, M and C are
constant matrices while the stiffness matrix K is varied
with time for a nonlinear system.
The Newmark method can be used to obtain the
approximate solutions of Eq. (1) and its general formulation
is
Downloaded by [National Institute of Technology - Rourkela] at 09:17 13 February
2014
652 Journal of the Chinese Institute of Engineers, Vol. 27, No. 5 (2004)
Mai+1+Cvi+1+Kdi+1=fi+1
di+1=di+(?t)vi+(?t)2[(12
-)ai+ai+1]
vi+1=vi+(?t)[(1-?)ai+?ai+1] (2)
where di, vi and ai are approximations to u(ti), u(ti)
and u(ti), respectively, and and ? are specified coefficients
which govern the stability, accuracy and numerical
dissipation of the integration method. To obtain
the basic properties, it is convenient to consider the
single degree of freedom analogue of Eq. (2)
mai+1+cvi+1+kdi+1=fi+1
di+1=di+(?t)vi+(?t)2[(12
-)ai+ai+1]
vi+1=vi+(?t)[(1-?)ai+?ai+1] (3)
where di, v i and ai are approximations of
displacement, velocity and acceleration at the i-th
time step; m, c and k are the mass, viscous damping
and stiffness.
The coefficients for the four well-known members
of Newmark method (Newmark, 1959), which
are the Newmark explicit method, Fox-Goodwin
method, linear acceleration method and constant average
acceleration method, are summarized in Table
1 for the subsequent discussions and comparisons. In
addition, their stability limits for linear elastic systems
are also listed for comparison. In this table, the
relationship O=?(?t) is defined and ?= k/ m is the
natural frequency of the system.
III. AN IMPLEMENTATION OF NEWMARK
METHOD
In order to evaluate the numerical properties of
Newmark method in the solution of a nonlinear system,
the step-by-step integration details within a time step
are described and then are realistically reflected in
the proposed evaluation procedure. In the following
derivations, ki is used to represent the stiffness at the
end of the i-th time step or at the beginning of the
(i+1)-th time step. Similarly, ki+1 is used to denote
the stiffness at the end of (i+1)-th time step or at the
beginning of the (i+2)-th time step. This implies that
the stiffness during the (i+1)-th time step is varied
from the beginning of ki to the end of ki+1.
Immediately following the end of the i-th time
step, computation details for the (i+1)-th time step
are described next. The first step is to find the displacement
for the next time step. In fact, using the
second line of Eq. (3), the acceleration ai+1 can be
rewritten as a function of di+1. Meanwhile, the velocity
vi+1 can be expressed in terms of di+1 after substituting
the obtained ai+1 into the third line of Eq.
(3). As a result, the displacement di+1 can be obtained
after the substitutions of the derived vi+1 and ai+1 into
the first line of Eq. (3) and are found to be
[ 1
(?t)2m + ?
(?t)c + k i]di + 1
= f i + 1 + [ m
(?t)2 ][di + (?t)vi + (12
)(?t)2ai]
c [ vi + (1 ? )(?t)ai] + c[ ?
(?t) ][di + (?t)vi
+ (12
)(?t)2ai] (4)
In this equation, it is important to note that the stiffness
at the end of the i-th time step ki is used to determine
the displacement di+1 for the (i+1)-th time step.
Next, the restoring force ri+1 corresponding to this
displacement can be computed by an assumed material
model, which mathematically describes a stressstrain
relation or a load-displacement relation. In
general, the load-displacement relation for the (i+1)-
th time step can be expressed as ri+1=ki+1di+1 where
ki+1 is the stiffness at the end of the (i+1)-th time step.
After updating the stiffness ki+1, the velocity can be
computed based on the new determined stiffness ki+1
and are found to be
vi + 1 = {[ 1
? (?t) ]m + c} 1{ fi + 1 k i + 1di + 1
+ [ m
? (?t) ][vi + (1 ? )(?t)ai]} (5)
Table 1 Four well-known members of Newmark method
Method ? Stability condition
Newmark explicit method 0 12
0<O=2
Fox-Goodwin method 1
12
12
0<O= 6
Linear acceleration method 16
12
0<O=2 3
Average acceleration method 14
12
0<O=8
Downloaded by [National Institute of Technology - Rourkela] at 09:17 13 February
2014
S. Y. Chang: Studies of Newmark Method for Solving Nonlinear Systems: (I) Basic
Analysis 653
Finally, the acceleration can be computed by
ai + 1 = (1m
)( fi + 1 cvi + 1 k i + 1di + 1) (6)
In this computing procedure, the stiffness ki is used
to determine the displacement di+1 while the velocity
and acceleration are computed by using the stiffness
ki+1. This reveals that the use of different stiffnesses
to compute the displacement, velocity and acceleration
are the major differences between the linear elastic
and nonlinear systems.
IV. ANALYSIS
It is often found that for a nonlinear system the
stiffness is ki at the beginning of the (i+1)-th time
step while it becomes ki+1 at the end of this time step.
Since ki+1 is generally different from ki for a nonlinear
system, a new parameter i+1 is used to describe
the change of the stiffness during the (i+1)-th time
step. For this purpose, the degree of nonlinearity i+1
is defined as the ratio of ki+1 over ki
ki+1=i+1ki (7)
Obviously, there is no stiffness change for i+1=1 during
the (i+1)-th time step. Meanwhile, the case of
i+1>1 simulates the stiffness hardening, and the stiffness
softening can be mimicked by 0<i+1<1. Thus,
i+1 can be considered as a measure of the degree of
nonlinearity for the (i+1)-th time step. In general,
for i+1>1 the larger i+1 is, the higher the degree of
hardening nonlinearity is while for 0<i+1<1 the
smaller i+1 is, the higher the degree of softening
nonlinearity is.
After defining the degree of the nonlinearity i+1
for the (i+1)-th time step, the integration procedure
as shown in Eq. (4) to Eq. (6) can be rewritten in a
matrix form (Bathe and Wilson, 1973; Belytschko and
Hughes, 1983) as
Xi+1=Ai+1Xi+Li+1fi+1 (8)
where Xi=[di, (?t)vi, (?t)2ai]T; Ai+1 is the amplification
matrix and Li+1 is the load vector for the (i+1)-th
time step. In this study, the viscous damping is assumed
to be zero in order to simplify the algebraic
manipulations. The explicit expression of the matrix
Ai+1 is
Ai + 1 = 1
1 +Oi
2
1 1 12

i + 1?Oi
2 1 + ( i + 1? )O
i
2 (1 ? )(1 +Oi
2) (12 )(i + 1?Oi
2)
i + 1Oi
2 i + 1Oi
2 (12
)i + 1Oi
2
(9)
where Oi=?i(?t); and ?i= k i /m is the natural frequency
of the system determined from the structural
properties at the end of the i-th time step. The three
eigenvalues of the amplification matrix Ai+1 can be
found from the characteristic equation (Belytschko
and Hughes, 1983)
det (Ai+1-?I)=?3-A1?2+A2?-A3=0 (10)
where
A1 =
2 + ( + i + 1 i + 1? 12
i + 1)O
i
2
1 +Oi
2
A2 =
1 +i + 1( ? + 12
)O
i
2
1 +Oi
2
A3=0 (11)
After solving the characteristic equation, numerical
properties of the (i+1)-th time step for the Newmark
method, such as the local truncation error,
consistency, stability, accuracy and overshoot, can be
further evaluated for both linear and nonlinear
systems.
1. Local Truncation Error and Order of Accuracy
Since the object of numerical integration is to
determine sufficiently accurate approximations with
minimum effort, a means to evaluate the efficiency
of various approximation methods is required. Local
truncation error can be used as a criterion for this
purpose. The local truncation error is defined as the
error committed in each time step by replacing the
differential equation by the difference equation
(Belytschko and Hughes, 1983). In general, the local
truncation error for the (i+1)-th time step can be
expressed as
Ei + 1 = S [Hpu(p)(?t) p 2] p = 0
n
+ O[(?t) p 1] (12)
Downloaded by [National Institute of Technology - Rourkela] at 09:17 13 February
2014
654 Journal of the Chinese Institute of Engineers, Vol. 27, No. 5 (2004)
where u(p) denotes the p-th derivative of the displacement
u with respect to time and the coefficient Hp is
found to be
H0=1-A1+A2
H p =
1 + ( 1)pA2
p! p=1, 2, 3, ... (13)
The difference equation with local truncation error
Ei+1(?t) is said to be consistent with the differential
equation if
lim
?t? 0
( max
1 = i + 1 = N
Ei + 1(?t) ) = 0 (14)
( max
1 = i + 1 = N
u(ti + 1) di + 1 ) = 0 (15)
where di+1 is the approximation obtained from the difference
equation and u(ti+1) is the exact solution for
the differential equation.
After the substitutions of A1 and A2 into Eq. (12),
the local truncation error for the (i+1)-th time step of
the Newmarkwhere N is the total number of integration time steps.
If Ei+1=O[(?t)q] where q>0, then q is named the order
of accuracy (Belytscwhere N is the total number of integration time steps.
If Ei+1=O[(?t)q] where q>0, then q is named the order
of accuracy (Belytschko and Hughes, 1983). It is apparent
that consistency implies that Ei+1?0 as ?t?0,
and the order of accuracy is the rate at which the local
truncation error approaches zero as the time step
tends to zero. Furthermore, the difference equation
is convergent with the differential equation if
limwhere N is the total number of integration time steps.
If Ei+1=O[(?t)q] where q>0, then q is named the order
of accuracy (Belytschko and Hughes, 1983). It is apparent
that consistency implies that Ei+1?0 as ?t?0,
and the order of accuracy is the rate at which the local
truncation error approaches zero as the time step
tends to zero. Furthermore, the difference equation
is convergent with the differential equation if
lim
?t? 0where N is the total number of integration time steps.
If Ei+1=O[(?t)q] where q>0, then q is named the order
of accuracy (Belytschko and Hughes, 1983). It is apparent
that consistency implies that Ei+1?0 as ?t?0,
and the order of accuracy is the rate at which the local
truncation error approaches zero as the time step
tends to zero. Furthermore, the difference equation
is convergent with the differential equation if
lim
?t? 0
?t? 0hko and Hughes, 1983). It is apparent
that consistency implies that Ei+1?0 as ?t?0,
and the order of accuracy is the rate at which the local
truncation error approaches zero as the time step
tends to zero. Furthermore, the difference equation
is convergent with the differential equation if
lim
?t? 0 method is found to be
Ei + 1 = i + 1 1
1 +Oi
2?1
2u +
i + 1( ? + 12
)
1 +Oi
2 ?
i
Oi
u
+
1
12 12
[ + i + 1( ? + 12
)]
1 +Oi
2 ?
i
2Oi2u
+ O [(?t)3] (16)
In addition, for i+1=1, it reduces to be
Ei + 1 =
? 12
1 +Oi
2?
i
Oi
u +
1
12 + 12
(? 12
)
1 +Oi
2 ?
i
2Oi2u
+ O [(?t)3] (17)
This equation clearly shows that for i+1=1, the order
of accuracy is one if ??12
, two if ?=12
and three if =
1
12 in addition to the case of ?=12
. Based on Eq. (14),
consistency cannot be achieved for any degree of
hardening or softening nonlinearity since i+1?1. In
fact, the order of accuracy degenerates to 0 in the
solution of a nonlinear system for any member of the
Newmark method and thus the accuracy significantly
deteriorates when compared to that for the solution
of a linear system. However, consistency and nonzero
order of accuracy might be achieved if i+1 is
close enough to 1. This is because the first and/or
second terms on the right hand side of Eq. (16) is
smaller than the subsequent term.
2. Overshooting Effect
Goudreau and Taylor (1972) discovered a very
peculiar property of the Wilson ? method. Although
this method is unconditionally stable for linear systems,
numerical experiments showed a tendency to overshoot
exact solutions significantly in the early response.
Therefore, the tendency to overshoot is an important
factor, which should be considered in an evaluation
of an implicit method. To assess the tendency of an
integration method to overshoot the exact solutions,
one can consider a single degree of freedom system
subject to any arbitrary initial conditions di and vi and
then compute di+1 and vi+1 as a function of Oi. Since
only convergent integration methods are considered
there will be no overshooting effect as Oi?0.
However, the behavior as Oi?8 will provide an indication
of the behavior of the high frequency modes
in a structural system. This is a typical condition when
an unconditionally stable implicit method is used to
solve a large multiple degree of freedom system.
The following two equations are directly obtained
from Eq. (8) after using the equation of motion
to eliminate the acceleration term ai
di + 1 = 1
1 +Oi
2{[1 + ( 12
)O
i2]di + (?t)vi}
vi + 1 = 1
1 +Oi
2{ {1 ? + ? i + 1 + [(1 ? )
+ ( 12
)? i + 1]O
i2}O
i2(
di?t
)
+ [1 + ( i + 1? )O
i2]vi} (18)
Results for the limiting condition of Oi?8 for these
two equations are found to be
di + 1 (1 1
2 )di
v i + 1 [? 1 + ( 1
2 1)? i + 1]O
i2(
di?t
) + (1 i + 1
?
)vi
(19)
The first line of this equation reveals that there is no
overshoot in displacement for either linear or nonlinear
systems using the Newmark method. However,
it has a tendency to overshoot quadratically in Oi in
Downloaded by [National Institute of Technology - Rourkela] at 09:17 13 February
2014
S. Y. Chang: Studies of Newmark Method for Solving Nonlinear Systems: (I) Basic
Analysis 655
the velocity equation due to the displacement term
(di/?t) for any degree of nonlinearity. It is worth noting
that the constant average acceleration method is
the only member of the Newmark method that shows
no overshoot in velocity for linear elastic systems.
This is because =14
, ?= 12
and i+1=1, (di/?t) in the
second line of Eq. (19) disappears. However, the
overshooting in velocity is still present if it is used to
solve a nonlinear system.
3. Stability
There is a zero eigenvalue for the Newmark
method since A3=0 in Eq. (12). Thus, this integration
method is spectrally stable (Belytschko and Hughes,
1983) if the other two eigenvalues of the amplification
matrix Ai+1 satisfy the following two inequalitie
-(A2+1)=A1=(A2+1) if -1=A2<1
-2<A1<2 if A2=1 (20)
In this derivation, the two eigenvalues can be real or
complex. Thus, the response is in a bounded oscillatory
form or an exponential decay form. Although a
stable computation can be guaranteed by Eq. (20), it
is advisable to satisfy more stringent stability
conditions, which restricts the two eigenvalues to be
complex conjugate, and thus a bounded oscillatory
response is obtained. As a result, two complex conjugates
and |?1,2|=1 leads to
A1 2 -4A2<0
0<A2=1 (21)
If any member of the Newmark method can satisfy
the stability conditions either defined in Eq. (20) or
(21) for any value of Oi then it possesses unconditional
stability for the (i+1)-th time step. On the other
hand, if the stability conditions are met only for a
certain interval of Oi then the integration method can
only have the conditional stability. To explore the
stability conditions of the Newmark method, Eqs. (20)
and (21) will be thoroughly studied next.
(i) Two Real or Two Complex Eigenvalues - Eq. (20)
(1) Unconditional stability
For the Newmark method to have the response
in a bounded oscillatory form or an exponential decay
form, unconditional stability conditions can be
easily determined after the substitution of A1 and A2
into Eq. (20) and are found to be:
i+1(-?)+=0 and -=i+1(-?+12
)= (22)
It is found that the second inequality of this equation
will be satisfied only if the degree of nonlinearity i+1
is smaller than or equal to 1 for the case of ?=12
and
?0 while for ?=12
and =0 the first inequality is not
met for any i+1>0. Thus, in the sub-family of the
Newmark method with ?=12
an unconditionally stable
algorithm may exist for solving softening or linear
systems but not for solving hardening systems.
(2) Conditional stability
On the other hand, the conditional stability limits
for the Newmark method can be also found from
Eq. (20) and are:
i+1(-?)+<0, i+1(-?+12
)< and
0<Oi= 2
i + 1(? )
(23)
where the second inequality is met only if i+1<1 for
?=12
and ?0. This implies there is no conditionally
stable method in the sub-family of ?=12
and ?0 for
the system with i+1=1. In fact, the Newmark explicit
method is the only method in the sub-family of ?=12
that can have conditional stability for any i+1 and its
stability limits are found to be 0<O
i=2/ i or 0<O
i+1=2.
(ii) Two Complex Conjugate Eigenvalues - Eq. (21)
(1) Unconditional stability
Using the inequality of A1 2 -4A2<0, the lower and
upper bounds of are found to be:
a = if i+1=1
a==b if i+1?1 (24)
where
a = 14
(? + 12
)2
a = 1
(1 i + 1)2{i + 1[i + 1(? + 12
) ? + 32
]
2 i + 1 i + 1(? + 12
) ? + 12
}
b = 1
(1 i + 1)2{i + 1[i + 1(? + 12
) ? + 32
]
+ 2 i + 1 i + 1(? + 12
) ? + 12
} (25)
Downloaded by [National Institute of Technology - Rourkela] at 09:17 13 February
2014
656 Journal of the Chinese Institute of Engineers, Vol. 27, No. 5 (2004)
where a is proved to be greater than or equal to zero
if the result of the square root is a real number. It is
clear that the result of the square root in Eq. (25) is a
real number for i+1>1 while for the case of 0<i+1<1
it must satisfy the following inequality:
? = 12
(
1 +i + 1
1 i + 1
) (26)
As a result, for a=0 and a =0 Eq.(24) implies =0.
Thus, the following results can be further obtained
from the constraint of 0<A2=1 based on =0
c==d and ?=12
if i+1>1
?=12
if i+1=1
=?-12
if 0<i+1<1 (27)
where
c=?-12
and d = (
i + 1
i + 1 1
)(? 12
) (28)
As a further step, the conditions to have unconditional
stability for the Newmark method can be summarized
from Eqs. (24) and (27) and are
a==d, ?=12
and ?=i + 1 + 1
4 if i+1>1
a = and ?=12
if i+1=1
a==b and ?=12
( 1 +i + 1
1 i + 1
) if 0<i+1<1 (29)
In this equation, the inequality ?=(i+1+1)/4 in the first
case of i+1>1 implies that the degree of nonlinearity
i+1 must be smaller than or equal to 1 if ?=12
.
Obviously, this is contradictory to the assumption of
i+1>1 for this case. Thus, it is concluded that any
member of the Newmark method with ?=12
cannot
possess the desired unconditional stability for any
degree of hardening nonlinearity.
To see the variations of a and d with respect to
i+1 for a given value of ?, the first line of Eq. (29) is
plotted in Fig. 1. The line with an arrow denotes the
possible range of , which is 14
=<8, to have unconditional
stability for i+1=1 as ?=12
. In addition, the
solid circles represent the lower bounds of as i+1=1
for the cases of ?=12
, 34
, 1, 54
and 32
; and their correspondingiod
shrinkage for any degree of nonlinearity. For each
curve, the relative period error is gradually increased
on
of Numerical Integration Methods in Elastodynamics,
Computer Methods in Applied Mechanics40, No. 3, pp. 417-421.
Newmark, N. M., 1959, A Method of Computation
for Structural Dynamics, Journal of Engineering
Mechanics Division, ASCE, Vol. 85, pp. 67-94.
Manuscript Received: May 02, 2003
Revision Received: Dec. 10, 2003
and Accepted: Jan. 20, 2004
Downloaded by [National Institute of Technology - Rourkela] at 09:17 13 February
2014