NONLINEAR ANALYSES OF

STRUCTURAL SYSTEMS

Computer Methods for the

Solution of Nonlinear Problems
D
Spring 2015

Prof.Dr. Bülent AKBAŞ

Departmant of Earthquake and Structural Engineering
 Outline

• Solution of Nonlinear Problems

• Solution Strategies

• Event-to-event method

• Iterative methods
 Solution of Nonlinear Problems

Nonlienar analysis seeks for finding the state of equilibrium of a

structure/element/component corresponding to the applied nodes, e.g. the
same as in linear analysis. Equation of equilibrium can be stated as:

Pt – Ft = 0 (1)

Pt =external forces applied at a point at time t

Ft =internal forces developed at that point

Equation 1 is satisfied in the deformed geometry of the structure and should

include all nonlinearities. Note that it will include inertia and damping
forces in case of a dynamic analysis.
 Solution of Nonlinear Problems

Nodal internal forces and displacements at certain load levels or at certain

times are needed in solving a nonlinear problem. Equation 1 should be
solved for the whole time range including nonlinear geometric or material
nonlinearity. The response of a structure can be calculated using a solution
strategy for the problem. If there is no nonlinearity, the solution is only a
one-step analysis.
State-of-the-practice of nonlinear solution strategies applies an incremental
solution with a number of load steps to finally reach the total applied load.
 Solution Strategies

Nonlinear problems can be solved by iteration of a series of linearized

analysis as discussed before. Many iteration strategies have been used in
literature. Some of them are:
a. Event-to-event strategy
b. Newton-Raphson Iteration
c. Modified Newton-Raphson Iteration
d. Initial Stiffness Iteration
e. Secant Stiffness Iteration
 Solution Strategies

A. Event-to-event strategy
i. Limitations of the method:
-all elements have multilinear force-deformation (F-D) relationship.
-small displacements are assumed (linear kinematics and
equlibrium in the undeformed configuration)

F P

k4 K4
k3 K3
k5 K5
k2 K2

k1 K1

D u

Element/component Structure
 Solution Strategies
A. Event-to-event strategy
ii. Definition of ‘’event’’
Each element stiffness change is called an Event. Structure’s stiffness
changes when an element stiffness changes. The stiffness of an element and
structure are constant between events. The obvious method of analysis is to
go from event to event solving a series of linear problems.
iii. Advantages of event-to-event strategy
-conceptually simple and reliable
-gets traces of the complete force-deformation and load-
displacement paths and simulates what happens in the actual
structures P

-linear solution 1 3

-converged DP
2

-can be scaled to event

Du
u
 Solution Strategies
A. Event-to-event strategy
iv. Disadvantages of event-to-event strategy
-only limited to multilinear element models

F F

k4
k3
k5
k2

k1

D D

Multilinear Element/component model Parabolic Element/component model

 Solution Strategies
A. Event-to-event strategy
iv. Disadvantages of event-to-event strategy
-if there are many elements with many possible element events, the
number of structure events can become large. Since a complete
linear analysis is required for each event, the computational cost
increase rapidly as the size of the structure increases.

Brace behavior

Structural behavior

We have to keep tracking the “events”

very often and have to scale back the
forces and update the stiffness matrix.
This can become computationally very
cumbersome.
 Solution Strategies
A. Event-to-event strategy
v. Basic steps in event-to-event strategy
1.initialize
2.define loading (e.g. Gravity, impact, dynamic)
3.linearize and solve
4.calculate event factor
5.scale to event
6.repeat steps 3-5 until all load is applied
7.define additional loading (e.g. seismic loading)
8.repeat steps 3 through 6.
 Solution Strategies
A. Event-to-event strategy
vi. Event factor calculation

Consider a typical element in a structural model as shown.

F
Solution to linear problem
Fev
DFL
k DFL  kDDL
Fcur State at the beginning of step (current state)

Dcur Dev D

DDL
Obtained from calculated Du by kinematics

Dev  Dcur Fev  Fcur

Event factor is defined as: f  or f 
DDL DFL

The governing event factor is the smallest event factor calculated

considering all elements.
 Solution Strategies
A. Event-to-event strategy
vi. Detailed algorithm using load increments
Initialize
F=0, D=0 for all elements
P=0, u=0 for structure

Set up DP, load increment

Form tangent stiffness (current

state) for each element, k

Assemble tangent stiffness for the

structure, K

Solve KDu=DP for Du, displacement increment

State determination and event factor calculations, f

Get smallest value of f

Replace Du by fDu

Update state of each element (see Note 1)

Replace DP by (1-f)DP
Yes No
f<1
(event occurs) (f=1, no event occurs)
 Solution Strategies
A. Event-to-event strategy
vii. Detailed algorithm using total load (used more often than load increments)
Initialize
F=0, D=0 for all elements
P=0, PI=0, u=0 for structure

Set up DP, load increment and add to P (total)

Form tangent stiffness (current

state) for each element, k

Assemble tangent stiffness for the

structure, K

Solve KDu=P-PI for Du, displacement increment

State determination and event factor calculations, f

Get smallest value of f

Replace Du by fDu

Update state of each element (see Note 1)

Replace PI
Yes No
f<1
(event occurs) (f=1, no event occurs)
 Solution Strategies
A. Event-to-event strategy
Note 1: It is important that the updated state of the critical element (the one
with smallest event factor, f) be just beyond the event. Because of round-of-
error, the state based on fDP could be just before the event (consider also the
case with f=0). This can be taken into account either:
-by adding a small number to the calculated f, or
-by calculating the event factor for a point just beyond the event. For
example:
F
Solution to linear problem
Fev DDtol=a small tolerance
DDtol
DFL
k

Fcur State at the beginning of step (current state)

Dev  DDtol   Dcur

D f 
DDL
Dcur Dev
DDL
Obtained from calculated Du by kinematics
Load increments can be used when expected displacement is large and structural
behavior is complex. Convergency within each load increment is quick. However,
number of steps decides the computing time. Sometimes total load can be applied and
iteration can be used if unbalanced load appears.
 Solution Strategies
A. Event-to-event strategy
viii. Some key computational steps in algorithm using total load
-Linearization KDu=P-PI
-given current states for all elements
-form tangent stiffness matrices for elements, k
-assmeble tangent stiffness matrix for the structure, K
-Event factor calculation
-given current states for all elements
-given Du
-calculate event factor, f, for each element, get the smallest.
-State determination
-given current states for all elements
-given Du
-calculate new states for all elements
-Resisting force calculation
-given current states for all elements
-calculate nodal forces for all elements
-assemble PI
 Example :
The rigid block with four springs is subjected to lateral force of R=25 kN. F-D
relationships of the springs are given as well. Plot P-u graph and find the
maximum displacement.
a. Assume linear force-deformation relation for the element
b. Assume nonlinear force-deformation relation for the element

u F (kN)
Rigid
block A B C D
A

P=25 kN B

D 10 10 10 10

0.5 1.0 1.5 D (m)

0.5m 0.5m 0.5m

 Solution :
Let D1=D2=D3=D4=u
Force-deformation relationhips for the springs:

For spring A : F=10D, DF=10DD

For spring B : D < 0.5m  F=0D, DF=0DD
D ≥ 0.5m  F=10(D-0.5), DF=10DD
For spring C : D < 1.0m  F=0D, DF=0DD
D ≥ 1.0m  F=10(D-1.0), DF=10DD
For spring D : D < 1.5m  F=0D, DF=0DD
D ≥ 1.5m  F=10(D-1.5), DF=10DD
 Solution :

Iteration using total load algorithm:

PI K
State u FA FB FC FD P-Pi kA kB kC kD Du fA fB fC fD f fDu
(=∑Fi) (=∑ki)

0 0 0 0 0 0 0 25 10 0 0 0 10 2.5 ∞ 0.2 0.4 0.6 0.2 0.5

1 0.5 5 0 0 0 5 20 10 10 0 0 20 1.0 ∞ ∞ 0.5 1.0 0.5 0.5

2 1.0 10 5 0 0 15 10 10 10 10 0 30 0.33 ∞ ∞ ∞ 1.5 1.0 0.33

3 1.33 13.33 8.33 3.33 0 25 0

PI   Fi Dev  Dcur
K .Du  P  PI f 
DDL

u j 1  u j  fDu j ( j : number of state) f  min( f i )

 Solution :

P (kN)

25
Linearization about State 1
Linearization about State 2 Linearization about State 0
(initial State)
20

15

10

0.5 1.0 1.33 1.5 2.0 2.5 u (m)

 Solution :

Calculation of event factors:

For State 0:

0 0.5  0 1.0  0 1.5  0

fA   fB   0.20 fC   0.40 fD   0.60
2.5 2.5 2.5 2.5

For State 1:

  0.5   0.5 1.0  0.5 1.5  0.5

fA   fB   fC   0.50 f D   1.0
1.0 1.0 1.0 1.0

For State 2:
  1.0   1.0   1.0 1.5  1.0
fA   fB   fC    fD   1.5  1.0
0.33 0.33 0.33 0.33
 Solution Strategies
B. Step-by-Step Solution
Incremental step-by-step solution assumes that the solution for the discrete
time t is known and that the solution for the discrete time t+Dt is required.
Recall the equation of equlibrium at time t:

Pt – Ft = 0 (1)

Pt =external forces applied at a point at time t

Ft =internal forces developed at that point

And the equation of equlibrium at time t+Dt becomes:

Pt+Dt – Ft+Dt = 0 (2)

Pt+Dt =external forces applied at a point at time t+Dt
(independent of the deformation-assumed)

Since we know the solution at time t:

Ft+Dt = Ft + F (3)
F =increment in internal forces developed at that point
 Solution Strategies
B. Step-by-Step Solution

F ≈ Kt u (4)

u =vector of incremental nodal point displacements

F t
K  t
t
(5)
u tangent stiffness matrix
(corresponding to the geometric and material conditions
at time t – approximately)

Substituting Equations 4 and 3 into 2, we get:

Kt u= Pt+Dt - Ft (6)

Solving for u, we can get the approximate displacements at time t+Dt:

ut+Dt = ut + u (7)

The exact displacements at t+Dt are those that correspond to the applied loads
Pt+Dt . Equation 6 is only an approximation because Equation 5 is used.
 Solution Strategies
B. Step-by-Step Solution

Knowing the approximate displacements corresponding to time t+Dt,

approximate stresses and corresponding nodal point forces at time t+Dt,
Ft+Dt can be solved, and then proceeded to the next time increment
calculations. Because of the assumption of Kt in Equation 5, such a solution
may be subject to very significant errors, and depending on the time of load
step size used, may indeed be unsuitable. In practice, it is necessary to
iterate until the solution of Equation 2 is obtained with sufficient accuracy.
 Solution Strategies
C. Newton-Raphson Iteration

Newton-Raphson method is widely used as an iteration method. This method

is an extension of the simple incremental technique givein in Equations 6
and 7. Since we already obtained an increment in the nodal point
displacements, we can now define a new total displacement vector and the
incremental solution presented above can be repeated using the currently
total displacements instead of the displacements at time t.

Newton-Raphson iteration equations for i=1,2,3,...

Kit1Dt Dui  Pt  Dt  Fi t 1Dt
DPi 1 (9)
uit  Dt  uit1Dt  Dui
And the initial conditions:
(10)
u0t  Dt  u t K 0t  Dt  K t F0t  Dt  F t
Equations 9 are obtained by linearizing the response of the structural model
about the conditions at time t+Dt, iteration (i-1).
 Solution Strategies
C. Newton-Raphson Iteration
P (external load)
This error should be within the
specified tolerance
K 0t  Dt
K1t  Dt
Applied load
t  Dt
P
Pt  Dt  F2t  Dt
Pt  Dt  F1t  Dt

P t  Dt  F0t  Dt
Ki= tangent stiffness

Pt
Exact solution Nodal
Displacement
u t Du Du2 u t  Dt
1

Figure . Newton-Raphson Iteration

 Solution Strategies
D. Modified Newton-Raphson Iteration
Modified Newton-Raphson method is also widely used as an iteration method
and lies in between the full Newton-Raphson iteration method and the
initial stiffness method. The Modified Newton-Raphson iteration equations
for i=1,2,3,...

K t Dui  Pt  Dt  Fi t 1Dt (11)

DPi 1
And the initial conditions:
u0t  Dt  u t F0t  Dt  F t
t corresponds to one of the accepted equlibrium configuration at times 0, Dt,
2Dt, or t. This method involves fewer stiffness updates than the full
Newton-Raphson iteration and bases the stiffness matrix update on an
accepted equlibrium configuration. The choice of time time steps when the
stiffness matrix should be updated depends on the degree of nonlinearity in
the structure. The more nonlinear the response, the more often the
updating should be performed. Convergence in this method is slower,
because the same stiffness matris is used throughour the same time step.
 Solution Strategies
D. Modified Newton-Raphson Iteration
P (external load)

K 0t  Dt t  Dt
K 0 t  Dt
K0
Applied load
t  Dt
P
Pt  Dt  F1t  Dt

Pt  Dt  F2t  Dt
t  Dt t  Dt K t  K 0t  Dt
P F0
Ki= tangent stiffness

Pt
Exact solution Nodal
Displacement
u t Du Du2 u t  Dt
1

Figure . Modified Newton-Raphson Iteration

 Solution Strategies
E. Initial Stiffness Iteration Method
In this method, only the initial stiffness K0 needs to be calculated, thus avoiding
the expense of recalculating and factorizing many times the equation 5. This
method corresponds to a linearization of the response about the initial
configuration of the ssytem may converge very slowly or even diverge. The
Initial Stiffness Iteration equations for i=1,2,3,...

K 0 Dui  Pt  Dt  Fi t 1Dt
DPi 1
And the initial conditions:
u0t  Dt  u t F0t  Dt  F t
 Solution Strategies
E. Initial Stiffness Method
P (external load)

K0 K 0K 0 K 0
Applied load
t  Dt
P
Pt  Dt  F1t  Dt

Pt  Dt  F2t  Dt

P t  Dt  F0t  Dt
K0 = initial stiffness

Pt
Exact solution Nodal
Displacement
uDt u Du u t  Dt
1 2

Figure . Initial Stiffness Method

 Solution Strategies
F. Secant Stiffness Iteration Method
Secant stiffness iteration is an approximation of the Newton-Raphson Iteration.

F F

Fcur
Fcur

ksec ksec

D D
Dcur Dcur
Multilinear Element/component model Parabolic Element/component model
Figure . F-D relationships that can be used in secant stiffness iteration

Linearization  F = ksecD
 Solution Strategies
F. Secant Stiffness Iteration Method
Advantages of secant stiffness iteration:
-computationally simple
-can be applied to any F-D relationship
-can be applied to F-D relationship with strain softening behavior
F
Tangant k is negative (may be a problem)

Fcur

Ksec (positive)

D
Dcur
Parabolic Element/component model

Disadvantages of secant stiffness iteration:

-may be slow to converge (or may not converge at all)
-it is difficult to construct secant stiffness for complex elements
 Solution Strategies
Discussion on Convergence:
The basic iteration strategies which are discussed in the previous parts are not
guaranteed to converge. The following examples show the failures of
solution strategies for P-u space.

Case 1: Failure of Secant Stiffness Iteration:

P Applied Load

1 1
K sec

Event-to-event and Newton-Raphson work for this case.

 Solution Strategies
Discussion on Convergence:
Case 2: Failure of Newton-Raphson Iteration:
P 1 K1

P
Applied Load

Goes to zero K0

K0 and K1 are used repeatedly here and iteration is appearantly not going
anywhere.
 Solution Strategies
Discussion on Convergence:
Case 3: Failure of Newton-Raphson and Event-to-Event Iterations:
P

1 6 3 Applied Load
P
Real solution

4
2 7
8
5

Event-to-event iteration:
-advance to 1
-back to 2
-advance to 3
-back to 4
-advance to 6 (should have gone to 5)

The iteration repeats itself in the cycle and never gets to the real solution.
 Solution Strategies
Discussion on Convergence:

Newton-Raphson iteration:
-advance to 1
-back to 8
-advance from 8 to 3
-back to 5
-advance to 6 (should have gone to the right)

The real P-u relationship is unknown in particular for comlex structures.

Therefore, which iteration is used to get convergence is more experimental
than analytical. In practice, different strategies should be tried. More
choices in strategies in a pregram are desirable.
 Solution Strategies
Discussion on Event-to-Event and Iteration Methods:

Event-to-Event Method:
- tries to track down the “events” in each member
-error between internal and external forces at every step is zero
-can give accurate results, but may be computationally very cumbresome
for large structures

Iterative Methods:
-are based on error minimization between applied and resisting force
-a small error always exists (a small tolerance should always be set)
-may not be as accurate as the event-to-event method, but computationally
takes definitely less time for large structures.
 References

• Shen,J., Seismic Retrofit and Earthquake Hazard Reduction,

Class Notes.
• Graham Powell, A Presentation of CSI Educational Services,
Performance Based Design Using Nonlinear Analysis: Perform
3D Özer.
• Bathe, KJ, Introduction to Finite Element Analyses.