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.

23

7.6 NEWTON“RAPHSON ITERATION AND SOLUTION PROCEDURE

7.6.2 LINE SEARCH METHOD

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

that,

|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)

24 DISCRETIZATION AND SOLUTION

R

R

R(1)

·1 1

·

·

·1

1

R(0)

R(0) R(1)

(a) (b)

FIGURE 7.3 Quadratic line search

which yields a value for · as,

2

± ± R(0)

·= ± ’ ±; ±= (7.72)

2 2 R(1)

If ± < 0, the square root is real and the current value for · emerges as,

2

± ±

·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.

7.6.3 ARC LENGTH METHOD

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

25

7.6 NEWTON“RAPHSON ITERATION AND SOLUTION PROCEDURE

sphere

F ’

‚1F ’

F

s

0 ’

A A ∆ ‚1F

B s

’

‚0 F

u

∆ x1

0

B

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-

26 DISCRETIZATION AND SOLUTION

straint equation is,

T

∆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,

∆x

s= (7.78)

∆»ψF

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 ;

(7.80a,b,c)

In (7.79) u and γ must be such that the constraint Equation (7.77) remains

satis¬ed, hence,

T

(∆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)

where,

T

a1 = uT uF + ψ 2 F F

F

T

a2 = 2uT (∆xk + uR ) + 2∆»k ψ 2 F F

F

a3 = uT (2∆xk + uR )

R

27

7.6 NEWTON“RAPHSON ITERATION AND SOLUTION PROCEDURE

There are two solutions γ (1) and γ (2) to (7.82), which when substituted

(1)

into (7.80a), (7.79c,e), and (7.77) give two revised generalized vectors sk+1

(2)

and sk+1 . The correct parameter, γ (1) or γ (2) , is that which gives the mini-

(j)

mum “angle” θ between sk and sk+1 where θ is obtained from,

(j)

sT sk+1 ∆xk + u(j)

∆xk (j)

k

(j)

cos θ = ; sk = ; =

sk+1

(∆»k + γ (j) )ψF

s2 ∆»k ψF

(7.83)

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 .

Exercises

1. Prove that the equivalent internal nodal forces can be expressed with

respect to the initial con¬guration as,

Ta = 0 Na dV

FS

(e)

V

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,

(e)

K σ,ab = ( 0 Na 0 Nb )I dV

·S

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,

3

‚Na sym ‚Nb

(e)

= FiI C FjL dV

K c,ab

‚XJ IJKL ‚XK

ij (e)