Raphson iterative process whereby given a solution estimate xk at iteration

k, a new value xk+1 = xk + u is obtained in terms of an increment u by

establishing the linear approximation,

R(xk+1 ) ≈ R(xk ) + DR(xk )[u] = 0 (1.34)

This directional derivative is evaluated with the help of the chain rule as,

d

DR(xk )[u] = R(xk + u)

d =0

d R1 (x1 + u1 , x2 + u2 )

=

R2 (x1 + u1 , x2 + u2 )

d =0

= Ku (1.35)

where the tangent matrix K is,

‚Ri

K(xk ) = [Kij (xk )]; Kij (xk ) = (1.36)

‚xj xk

If we substitute Equation (1.35) for the directional derivative into (1.34), we

obtain a linear set of equations for u to be solved at each Newton“Raphson

iteration as,

K(xk )u = ’R(xk ); xk+1 = xk + u (1.37a,b)

For equations with a single unknown x, such as Equation (1.11) for the

truss example seen in Section 1.3.2 where R(x) = T (x) ’ F , the above

Newton“Raphson process becomes,

R(xk )

u=’ ; xk+1 = xk + u (1.38a,b)

K(xk )

This is illustrated in Figure 1.12.

18 INTRODUCTION

T R

K

R(xk )

0

F x

u

x

xk+1 xk

FIGURE 1.12 Newton“Raphson iteration.

In practice the external load F is applied in a series of increments as,

l

F= ∆Fi (1.39)

i=1

and the resulting Newton“Raphson algorithm is given in Box 1.1 where bold-

face items generalize the above procedure in terms of column and square

matrices. Note that this algorithm re¬‚ects the fact that in a general FE

program, internal forces and the tangent matrix are more conveniently eval-

uated at the same time. A simple FORTRAN program for solving the one-

degree-of-freedom truss example is given in Box 1.2. This program stops

once the sti¬ness becomes singular, that is, at the limit point p. A tech-

nique to overcome this de¬ciency is dealt with in Section 7.5.3. The way

in which the Newton“Raphson process converges toward the ¬nal solution

is depicted in Figure 1.13 for the particular choice of input variables shown

in Box 1.2. Note that only six iterations are needed to converge to values

within machine precision. We can contrast this type of quadratic rate of

convergence with a linear convergence rate, which, for instance, would re-

sult from a modi¬ed Newton“Raphson scheme based, per load increment,

on using the same initial sti¬ness throughout the iteration process.

19

1.4 DIRECTIONAL DERIVATIVE AND LINEARIZATION

BOX 1.1: NEWTON“RAPHSON ALGORITHM

INPUT geometry, material properties, and solution parameters

r

INITIALIZE F = 0, x = X (initial geometry), R = 0

r

FIND initial K (typically (1.13))

r

LOOP over load increments

r

r FIND ∆F (establish the load increment)

r SET F = F + ∆F

r SET R = R ’ ∆F

r DO WHILE ( R / F > tolerance )

SOLVE Ku = ’R (typically (1.38a))

r

UPDATE x = x + u (typically (1.38b))

r

FIND T (typically (1.12)) and K (typically (1.13))

r

FIND R = T ’ F (typically (1.11))

r

r ENDDO

r ENDLOOP

1

Residual norm

Modified Newton’Raphson

Newton’Raphson

’10

1e

2 4 6 8 10 12 14 16 18 20 Iterations

FIGURE 1.13 Convergence rate.

20 INTRODUCTION

BOX 1.2: SIMPLE TRUSS PROGRAM

c------------------------------------------------------

program truss

c------------------------------------------------------

c

c Newton-Raphson solver for 1 D.O.F. truss example

c

c Input:

c

c d ---> horizontal span

c x ---> initial height

c area ---> initial area

c e ---> Young modulus

c nincr ---> number of load increments

c fincr ---> force increment

c cnorm ---> residual force convergence norm

c miter ---> maximum number of Newton-Raphson

c iterations

c

c------------------------------------------------------

c

implicit double precision (a-h,o-z)

double precision l,lzero

c

data d,x,area,e,f,

resid /2500.,2500.,100.,5.e5,0.,0./

data nincr,fincr,cnorm,miter /1,1.5e7,1.e-20,20/

c

c initialize geometry data and stiffness

c

lzero=sqrt(x**2+d**2)

vol=area*lzero

stiff=(area/lzero)*e*(x/lzero)*(x/lzero)

c

c starts load increments

c

do incrm=1,nincr

f=f+fincr

resid=resid-fincr

21

1.4 DIRECTIONAL DERIVATIVE AND LINEARIZATION

BOX 1.3:

c

c Newton-Raphson iteration

c

rnorm=cnorm*2

niter=0

do while((rnorm.gt.cnorm).and.(niter.lt.miter))

niter=niter+1

c

c find geometry increment

c

u=-resid/stiff

x=x+u

l=sqrt(x*x+d*d)

area=vol/l

c

c find stresses and residual forces

c

stress=e*log(l/lzero)

t=stress*area*x/l

resid=t-f

rnorm=abs(resid/f)

print 100, incrm,niter,rnorm,x,f

c

c find stiffness and check stability

c

stiff=(area/l)*(e-2.*stress)*(x/l)*(x/l)

+(stress*area/l)

if(abs(stiff).lt.1.e-20) then

print *, ™ near zero stiffness - stop™

stop

endif

enddo

enddo

stop

100 format(™ increment=™,i3,™ iteration=™,i3,™

residual=™,g10.3,

& ™ x=™,g10.4,™ f=™,g10.4)

end

CHAPTER TWO

MATHEMATICAL PRELIMINARIES

2.1 INTRODUCTION

In order to make this text su¬ciently self-contained, it is necessary to include

this chapter dealing with the mathematical tools that are needed to achieve

a complete understanding of the topics discussed in the remaining chapters.

Vector and tensor algebra is discussed, as is the important concept of the

general directional derivative associated with the linearization of various

nonlinear quantities that will appear throughout the text.

Readers, especially with engineering backgrounds, are often tempted to

skip these mathematical preliminaries and move on directly to the main

text. This temptation need not be resisted, as most readers will be able to

follow most of the concepts presented even when they are unable to under-

stand the details of the accompanying mathematical derivations. It is only

when one needs to understand such derivations that this chapter may need

to be consulted in detail. In this way, this chapter should, perhaps, be ap-

proached like an instruction manual, only to be referred to when absolutely

necessary. The subjects have been presented without the excessive rigors of

mathematical language and with a number of examples that should make

the text more bearable.

2.2 VECTOR AND TENSOR ALGEBRA

Most quantities used in nonlinear continuum mechanics can only be de-

scribed in terms of vectors or tensors. The purpose of this section, however,

is not so much to give a rigorous mathematical description of tensor alge-

bra, which can be found elsewhere, but to introduce some basic concepts

1

2 MATHEMATICAL PRELIMINARIES

x3

v3

e3

v

e2