. 42
( 54 .)


solution is independent of the manner in which the load increments are
applied. If, however, the material is not hyperelastic this conclusion may
not hold.
An outline of the complete solution algorithm is shown in Box 7.1.

Remark 7.6.1. By comparing Equation (7.63) with the example given in
Section 2.3.4 relating to the linearization of a system of algebraic equations,
it is evident that the tangent sti¬ness matrix can be found directly as,

‚R ‚Ri
K= ; Kij = ; i, j = 1, ndof (7.65)
‚x ‚xj

where ndof is the number of degrees of freedom in the problem. For some sit-
uations, such as ¬nite deformation thin shell analysis, where the discretiza-
tion of the kinematic quantities is very much algorithm-dependent, such a
direct approach, though tedious, may be the only method of obtaining the
tangent matrix coe¬cients.

The Newton“Raphson process is generally capable of reaching the conver-
gence of the equilibrium equations in a small number of iterations. Nev-
ertheless, during the course of complex deformation processes, situations
may occur where the straight application of the Newton“Raphson method
becomes insu¬cient. A powerful technique often used to improve the conver-
gence rate is the line search method. This technique consists of interpreting
the displacement vector u obtained from Equation (7.63) as an optimal di-
rection of advance toward the solution but allowing the magnitude of the
step to be controlled by a parameter · as,
xk+1 = xk + ·u (7.66)
The value of · is normally chosen so that the total potential energy,
Π(·) = Π(xk + ·u), at the end of the iteration is minimized in the direction
of u. This is equivalent to the requirement that the residual force at the
end of each iteration, that is, R(xk + ·u), is orthogonal to the direction of
advance u. This yields a scalar equation for · as (see Figure 7.3),
R(·) = uT R(xk + ·u) = 0 (7.67)
Because of the extreme nonlinearity of the function R(·), Condition (7.67)
is too stringent and in practice it is su¬cient to obtain a value of · such
|R(·)| < ρ|R(0)| (7.68)
where, typically, a value of ρ = 0.5 is used. Under normal conditions the
value · = 1, which corresponds to the standard Newton“Raphson method,
automatically satis¬es Equation (7.68), and therefore few extra operations
are involved. Occasionally, this is not the case, and a more suitable value of
· must be obtained. For this purpose it is convenient to approximate R(·)
as a quadratic in ·. To achieve this requires the knowledge of the value
R(0), together with the derivative dR/d· at · = 0, which, recalling Remark
1, is obtained from Equation (7.67) as,
dR ‚R
= uT u = uT K(xk ) u = ’uT R(xk ) = ’R(0) (7.69)
d· ‚x x=xk
and, ¬nally, a third value, which typically is the standard value of the resid-
ual force for which · = 1, as,
R(1) = uT R(xk + u) (7.70)
The quadratic approximation thus obtained with these coe¬cients gives,
R(·) ≈ (1 ’ ·)R(0) + R(1)· 2 = 0 (7.71)



·1 1
R(0) R(1)

(a) (b)

FIGURE 7.3 Quadratic line search

which yields a value for · as,
± ± R(0)
·= ± ’ ±; ±= (7.72)
2 2 R(1)
If ± < 0, the square root is real and the current value for · emerges as,
± ±
·1 = + ’± (7.73)
2 2
Alternatively, if ± > 0 (see Figure 7.3b), then · can be simply obtained
by using the value that minimizes the quadratic function, that is, ·1 =
±/2. Within each (DO WHILE) iteration, the quadratic approximation
procedure is repeated with the three values, R(0), R (0), and R(·1 ), until
Equation (7.68) is satis¬ed.

Although the line search method will improve the convergence rate of the
Newton“Raphson method, it will not enable the solution to traverse the so-
called limit points on the equilibrium path. Figure 7.4(a) shows two such
limit points A and B on an example of snap-back behaviour. If the load is
incremented the solution will experience convergence problems in the neigh-
bourhood of A and may jump to the alternative equilibrium position A . In
many cases the equilibrium path can be followed beyond A by prescribing
a characteristic displacement and calculating the load as the correspond-
ing reaction. In Figure 7.4(a) this technique would enable the solution to

F ’
‚1F ’
0 ’
A A ∆ ‚1F

B s

‚0 F

∆ x1
x x
x0 x1 x2
(a) (b)

FIGURE 7.4 Snap-back and (b) spherical arc length method

progress to the neighborhood of the limit point B, but, again, the solution
is likely to jump to point B .
Various ad hoc schemes combining load and displacement incrementation
have been devised to deal with limit point problems, but these have been
superseded by arc length methods that constrain the iterative solution to
follow a certain route toward the equilibrium path. This is achieved by
controlling the magnitude of the external loads by a parameter » so that
the equilibrium equations become,

R(x, ») = T(x) ’ »F = 0 (7.74)

where F is a representative equivalent nodal load. The value of » is now
allowed to vary during the Newton“Raphson iteration process by imposing
an additional constraint equation. As a consequence of the proportional
loading, load increment i is de¬ned by an increment in the value of » over
the value at the end of the previous load increment i ’ 1 as,

∆Fi = ∆»F; ∆» = » ’ »i’1 (7.75a,b)

Similarly, the total change in position over the load increment is denoted
∆x, that is,

∆x = x ’ xi’1 (7.76)

A number of arc length methods have been proposed by imposing dif-
ferent constraint equations to evaluate the additional unknown ». A robust
candidate is the spherical arc length method in which the additional con-

straint equation is,
∆xT ∆x + ∆»2 ψ 2 F F = s2 (7.77)
Figure 7.4(b) illustrates this equation where the constant s can be thought
of as the magnitude of the generalized vector,
s= (7.78)
which de¬nes a generalized sphere which is the surface upon which the iter-
ative solution is constrained to move as convergence toward the equilibrium
path progresses. The term ψ 2 is a scaling factor that, in principle at least,
renders (7.77) dimensionally consistent.
The Newton“Raphson process is now established by linearizing Equa-
tion (7.74), taking into account the possible change in », to give a set of
linear equations at iteration k as,
R(xk , »k ) + K(xk )u ’ γF = 0 (7.79a)
where u represents the iterative change in position and γ denotes the iterative
change in » as,
xk+1 = xk + u; ∆xk+1 = ∆xk + u
»k+1 = »k + γ; ∆»k+1 = ∆»k + γ
Solving the above equations gives the iterative displacements u in terms of
the, as yet, unknown parameter change γ and the auxiliary displacements
uR and uF as,

uR = K(xk )’1 R(xk , »k ); uF = K(xk )’1 F
u = uR + γuF ;
In (7.79) u and γ must be such that the constraint Equation (7.77) remains
satis¬ed, hence,
(∆xk + u)T (∆xk + u) + (∆»k + γ)2 ψ 2 F F = s2 (7.81)
Expanding this equation using (7.77) and (7.80) gives a quadratic equation
for the unknown iterative load parameter change γ as,
a1 γ 2 + a2 γ + a3 = 0 (7.82a)
a1 = uT uF + ψ 2 F F
a2 = 2uT (∆xk + uR ) + 2∆»k ψ 2 F F
a3 = uT (2∆xk + uR )

There are two solutions γ (1) and γ (2) to (7.82), which when substituted
into (7.80a), (7.79c,e), and (7.77) give two revised generalized vectors sk+1
and sk+1 . The correct parameter, γ (1) or γ (2) , is that which gives the mini-
mum “angle” θ between sk and sk+1 where θ is obtained from,
sT sk+1 ∆xk + u(j)
∆xk (j)
cos θ = ; sk = ; =
(∆»k + γ (j) )ψF
s2 ∆»k ψF
In practice, where many degrees of freedom exist, the scaling factor ψ
is taken as zero, or, alternatively, if the representative load F is normalized
and ψ = 1 the second term in the constraint Equation (7.77) reduces to
∆»2 .
1. Prove that the equivalent internal nodal forces can be expressed with
respect to the initial con¬guration as,

Ta = 0 Na dV
and then validate this equation by recalculating the equivalent nodal
forces found in Example 7.3.
2. Prove that the initial stress matrix can be expressed with respect to the
initial con¬guration as,
K σ,ab = ( 0 Na 0 Nb )I dV
V (e)
and then validate this equation by recalculating the initial stress matrix
K σ,12 found in Example 7.5.
3. Show that the constitutive component of the tangent matrix can be ex-
pressed at the initial con¬guration as,

‚Na sym ‚Nb
= FiI C FjL dV
K c,ab
ij (e)


. 42
( 54 .)