Composite Newton-PCG and quasi-Newton iterations for nonlinear consolidation

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 86 (1991) 27-60 NORTH-HOLLAND

COMPOSITE NEWTON-PCG AND QUASI-NEWTON ITERATIONS FOR NONLINEAR CONSOLIDATION Ronaldo I. BORJA Department of Civil Engineering, Stanford University, Stanford, CA 94305, USA Received 8 September 1987 Revised manuscript received 15 December 1989

Newton's method is theoretically attractive for solving a system of nonlinear equations, but it is difficult to use in practice because it requires the solution of a linear system at each iteration. For large-scale problems where the dimension of the system may be several thousand, the solution of such linear system may be extremely expensive. The burden of equation-solving is more evident in two- and three-dimensional consolidation problems due to the coupling of soil skeletal displacements and pore fluid pressure which tends to increase the size of the matrix equation significantly. Two iterative techniques for solving such system of nonlinear equations are thus investigated in this paper and compared with the performance of the standard Newton's method in which the linear system is solved by direct triangular factorization. The first of these two techniques is a combination of Newton iteration for the nonlinear system and preconditioned conjugate gradients (PCG) with a global elastic tangent operator preconditioner for the linearized system. The PCG iterative linear equation solver is shown to dramatically improve the efficiency of Newton's method, particularly when the system of linear equations is large. The second technique involves a secant-Newton iteration in which the rank-two BFGS inv~,rse update is employed. For large symmetric systems it is shown that this technique can also be competitive with the standard Newton's method, although it appears to be incapable of handling nonsymmetric systems well. A simple line search algorithm based on the majorization principle is employed to force all of these iterative techniques to produce a norm-reducing solution. Numerical results employing realistic constitutive models for soils are then presented to compare the performance of each of these iterative techniques in the analysis of typic,'d nonlinear consolidation problems.

1. Introduction Nonlinear equations in finite dimensions are generally solved by iteration. The prototype of a majority of iterative techniques is the famous Newton's method which allows the tangent operator to vary from iteration to iteration. For systems where the number of unknown degrees of freedom n~q is large, a major impediment to finite element analysis of complex, nonlinear structures by Newton's method is the deficiency of direct equation solving techniques associated with the solution of the linearized problem. For example, an evaluation of the residual force vector and stiffness matrix necessary to carry out this iteration entails O(neq ) and O(n~q) scalar function evaluations, respectively, while the cost associated with direct equation solving techniques is O(n~q) arithmetic operations, which can well dominate the total computational expense for large systems [1, 2]. Much research effort in the past has been devoted to developing less CPU-intensive 0045-7825/91/$03.50 © 1991- Elsevier Science Publishers B.¥. (North-Holland)

_~

R. !. Borja, Composite Newton-PCG and quasi-Newton iterations

iterative algorithms for application to nonlinear problems in structural engineering and solid mechanics (see, e.g., [2-8]). However, it must be emphasized that the need for such alternative iterative strategies is far more pressing in the analysis of coupled or 'mixed' systems such as elastoplastic soil consolidation because of the coupling of various types of degrees of freedom which tends to increase the size of the problem dramatically. The enormity of the problem is demonstrated in Table 1 which compares the computing costs associated with triangular factorization of nonsymmetric stiffness matrices for meshes consisting of 100, 200 and 400 elements of the following types: (11 4-node fluid conduction element with one D O F per node, (2) 4-node quadrilateral continuum element with two D O F ' s per node, and (3) 9-node displacement/4-node pore pressure 'mixed' element with three D O F ' s at the corner nodes and two D O F ' s elsewhere. In three-dimensional coupled analyses, one can expect that the cost of triangular factorization will grow at an inconceivably faster rate. In contrast to Newton's method, an iteration method called modified Newton requires a reduced number of assembly and factorization of the global tangent operator. However, this method often leads to very slow convergence, if not total divergence, particularly when the degree of nonlinearity is high [5]. An alternative procedure is represented by a class of algorithms called quasi-Newton [1, 8] which involve rank-one or rank-two updates of the coefficient matrix to force it to satisfy a global secant equation at each iteration. This procedure often yields better search directions than does the modified Newton iteration, but at the expense of additional storage demands and computing effort associated with matrix updating. Table ! Total CPU required (in seconds) for pertorming C~rout triangular factorization of stiffness matrix (NEO = nunlber of equations): computing was done on a system called Macbeth (Section 5 ) I,

Ft~ur-node quadrilateral fluid conduction elements Mesh

NEQ

CPU, sec.

111x 1()elenlents !11x 2(Ielements 211x 2!1elements

121 231 441

11,2 (I.3 1,7

2. Four-node quadrilateral elasto-static continuum elements Mesh

NEO

CPU, sec,

!11x 10dements I() x 20 elements 211x 20 elements

242 462 882

i,0 2. ! 11,7

3. Nine-node displacement/four-node pore pressure mixed elements Mesh

NEO

CPU, see,

III x 10elements III x 211elements 211x 211elements

I,f1113 !,953 3,803

12331 458,2 4,892,0

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

29

In this paper we re-focus our investigation on Newton's method and look at alternative ways of improving its efficiency when used for nonlinear elasto-plastic soil consolidation analysis. In lieu of the conventional LDU-factorization of the global tangent operator [9], we investigate an iterative solution strategy based on preconditioned conjugate gradients (PCG) with an elastic tangent operator preconditioner for solving large systems of equations associated with the linearized problem. In the context of nonlinear consolidation, it is emphasized that the resulting global tangent operator is generally indefinite, i.e., its determinant may change in sign from step to step, and even from iteration to iteration [10]. In the context of computational plasticity, it is further emphasized that the tangent operator may be asymmetric even with an elasto-plastic soil model obeying the associative flow rule, since not all stress integration algorithms for elasto-plasticity are symmetry-preserving [11-13]. These two factors apparently make the implementation of the PCG technique less straightforward, but we show that one can simply choose appropriate search parameters and effectively employ the PCG technique essentially in its unaltered form. The resulting algorithm is now a composite Newton-PCG technique which, when tested in conjunction with a simple line search, shows remarkable efficiency and robustness. We also investigate the performance of a popular quasi-Newton technique, the rank-two BFGS inverse update, for the solution of nonlinear consolidation. Although it is well-known that accurate line searches enhance the convergence of this technique, we caution that line searches are generally more expensive in computational geomechanics since the constitutive models employed are usually more elaborate. We thus perform a line search only in situations where a unit steplength, in conjunction with the search direction obtained from the BFGS iteration, does not reduce the length of the residual vector during the iteration (thus a line search simply forces the iteration to be norm-reducing). Sufficient conditions exist to guarantee such norm-reducing solutions, and we propose a simple line search algorithm based on the majorization principle [14] to keep the required number of residual function evaluations to minimum. The BFGS iteration is then shown to be competitive with Newton's method when the system is symmetric, although it does not perform as well otherwise. As for notations and symbols, bold-face letters denote matrices and vectors, (.,.) is the inner product, and I1"11= (',')"" is the vector norm.

2. Matrix equation: Newton's method with line search

We begin this section with a review of the general recurrence formula for nonlinear consolidation based on Biot's original linear theory [15] which accommodates the one-step generalized trapezoidal family of methods and the linear multis'tep (LMS) kth order backward differentiation formula (BDF-k) methods as formulated in [16, 17]: NS(~n+l) + NW(O,+,) = (FExT),,+,,

G'

(8,

k

in which

)

, + , - ~, a,,,8,,+._,,, + ~,, Attb(O,,+O)= ~,, At(HEx-r),,+O , m=!

(2.1a) (2.1b)

30

R.i. Borja, Composite Newton-PCG and quasi-Newton iterations

~,,+~ -/3¢k(O,,+ ~) + (1 -/3)¢k(0,,), (HExT),,+~ -- ~[~(HExT)n+ 1 + (1

-

(2.2a)

fl)(HExT)

(2.2b)

n ,

where N"(/$) and N*(O) are the vectors of internal nodal forces due to intergranular stresses and pore pressures, respectively, G' is the global connectivity array which couples the displacement and pore pressure degrees of freedom, ~k(0) is the nodal vector of internal flow, 8 is the nodal displacement vector describing the motion of the soil matrix, 0 is the nodal vector interpolating the pore fluid pressure, FEXr and HEXx are the applied nodal force and fluid supply vectors, respectively, At = t,,+. - t,, is the step size, 0 ~
(y~×~).,+,-#"(8.+,)-#*(e.÷,) t

f = f(d. +, ) =

}

(2.3a)

¢,, At(Hsx T),,+/~ - ]30 Atck(O. + o ) - ¢o G 0(8. +, )

where

8.+,

d,, + i

0. + i

0(8,,+s)= I

- ~ a,,8, ., ~ |

(2.3b)

In (2.3), f E R"~, is the global residual vector, 0 is the displacement difference operator [18, 19] and neq is the number of equations. The problem is then reduced to finding d,,+~ E R"eq such that

f(a.+,)=0.

(2.4)

Table 2 Coefficients of the generalized trapezoidal and constant step size backward differentiation formula methods (see [17] for the appropriate variable step size BDF coefficients) Method

/3

/3.

aI

a.

a~

Forward Euler Crank-Nicolson BDF-1 BDF-2 BDF-3

0 1/2 1 1 1

1 1 1 2/3 6/11

1 1 1 4/3 18/11

0 0 0 -1/3 -9/11

0 0 0 0 2/11

R.I. Borja, Composite Newton-PCG and quasi.Newton iterations

31

Throughout this paper, we will be concerned with an iteration of the form d,n+l +1

_

_

dl ' +1 + Ai Ad i '

(2.5)

Ad' = AT, tf(d~ +, ) ,

i where AtE L(Rneq). For example, with Newton's method, A i = - f t (d,+~) is the negative of the Jacobian o f f at the configuration d~+ ~. The scalar quantity Ai > 0 is a steplength parameter which forces the iteration (2.5) to converge uniformly. With Newton's method, sufficient conditions exist to ensure that this uniform convergence is achieved (see Section 2.2).

2.1. Newton's method and cons&tent tangent operator

Newton's method will play a central role throughout this paper and is, therefore, summarized in Box 1. In order to carry out the steps shown in Box 1, it is essential that f b e linearized. Linearizing f about the configuration d~+ t yields the consistent tangent operator [16, 17] __f

( d i, . l ) = A i

p

_.

I N s~(6,,+ i ~)

G

1

(2.6)

In general, NS'(8 ,,+ ~ !) is positive definite while ~'(0~, +~) is negative definite. Consequently, A i is indefinite. If N" (8,,+1) and ~'(0~,+~) are both symmetric, then A,. will also be symmetric. and ~ ' (0,+l) are symmetric or not, however, A~ will always Regardless of whether N~'(6 ,,+1) ~ have a symmetric profile associated with the connectivity of the finite elements. The symmetry of a tangent operator is central to numerous theorems and propositions in mathematical analysis. Furthermore, it is known that the factorization cost is significantly higher for a nonsymmetric matrix compared to its symmetric counterpart. Thus, it is important to determine conditions under which the tangent operator will remain symmetric so that appropriate algorithms may be employed. In the following we shall assume a symmetric tangent operator @'(0) and focus our attention exclusively on the operator N"(8). In the context of computational plasticity using finite elements, the use of a non-associated flow rule results in a non-symmetric tangent operator NS'(8). On the other hand, the use of an associative flow rule automatically guarantees a symmetric N~'(8) provided that an elastoplastic definition of tangent operator is employed [20]. However, it is now well recognized that such elasto-plastic tangent operator is known to destroy the important property of asymptotic quadratic convergence in Newton's method particularly when the step size is large [21]. Box 1 Newton's method with line search for nonlinear consolidation 1.

Initialize

i=0 f°ffif(d,,) Compute Ai = -f'(de. ÷ t ) If i > 1, line search for As; otherwise, ~1+ i + A~A[tf(di.÷~) -.+1I ~d,,+i

2. d ° ÷ t f d . , 3. 4. 5.

A~ = 1.0

6. fi÷i -'Jk ~d,÷! n+l) 7.

Update

i~i+l

8. If IIf'll/llflll <~,,, exit;

else,

go to Step 3.

32

R.I. Borja, Composite Newton-PCG attd quasi-Newton iterations

In [13.21] it was shown that a class of implicit stress update algorithms of the form o',,+ ~ = ~r(e,,+~; e,,, o'., X.; At),

(2.7)

where the o"s are stresses, e's are strains and the X's are some set of plastic variables, can be linearized consistently in the sense that C = Oo'.+~/Oe.+~,

(2.8)

so that, when C is used in the finite element expression f

N' '(8,,+,)= |J,, BtCB dO

(2.9)

the Newton-Raphson iteration converges quadratically. In (2.9), C is called the matrix of consistent tangential moduli which generally depends on the choice of update procedure. Note that such update procedure is not necessarily symmetry-preserving even if the plastic flow is associative [ 11-13]. In [13] it was shown that for perfect plasticity, only the fully implicit stress update scheme is symmetry-preserving (e.g., radial return [22] and closest point projection [23]). For hardening plasticity, additional restrictions must be imposed on the constitutive equation to ensure a symmetric matrix of consistent tangential moduli, among which are [13] (a) that the direction of the hardening moduli be evaluated from the same potential function used for determining the direction of plastic flow, and (b) that the set of plastic variables be updated in the additive form X,,. ~= ,Y,, + ~X, Although this scenario appears rather restrictive, it is affirmed by the fact that the matrix of consistent tangential moduli is not symmetric even for a relatively simple modified Cam-Clay plasticity model which obeys the associative flow rule but does not satisfy restrictions (a) and (b) above [11, 12]. We thus conclude that when applied to more realistic elasto-plastic constitutive models, the consistent tangent operator N~'(6) is in general nonsymmetric. However, it is important that we distinguish between two types of asymmetry: (a) strong asymmetry, which is due to the use of a non-associated flow rule, and (b) weak asymmetry, which is due to exact linearization of the update algorithm in conjunction with an associative flow rule. Such distinction follows from the fact that as e,+m---, e, (i.e., as the step size becomes small), C-.-,C ~p (i.e., the consistent tangential moduli approach the elastoplastic moduli [13, 21]), and thus a weakly asymmetric operator eventually becomes symmetric. The significance of these definitions will become apparent when iterative methods for symmetric systems are employed to solve weakly asymmetric systems of linear equations. 2.2. A simple line search: motivation and algorithm The steplength parameter A~was introduced in Step 5 of Box 1 to force Newton's method to produce a norm-reducing solution, and thus ensure convergence. The above statement is motivated by the following well-known lemma, which is central to all descent methods [24] (for brevity, we will omit the subscripts and superscripts in various terms and take on the following definitions: d = d l , + ~, Ad = Ad z and A = At).

R.i. Borja, Composite Newton-PCG and quasi-Newton iterations

33

L E M M A 2.1. Assume that a scalar function g: D C R"~,---> R ~is G-differentiable at d ~ int(D) and that (g'(d), Ad) < 0 for some Ad ~ R"~,; then there is a ,~ > 0, such that g(d + A Ad) < g(d)

VA U (O, i ) .

(,2.10)

PROOF. By definition of G-derivative, we have lim 1 [ g(d + ;t Ad) - g(d) - ,~(g'(d), Ad) ] = O. ^-.o -~ Since ( g ' ( d ) , A d ) < 0 , it follows that we may choose a ,~>0 so that (2.10) holds. Now, define the scalar function g =

Ilfll

(2.11) []

and employ (2.10) to obtain the following result.

C O R O L L A R Y 2.1. Assume that f :R""q---->R"~. is G-differentiable at d and that f ' ( d ) is invertible. If Ad = - [ f ' ( d ) ] - ' f ( d ) (Newton's method), then there exist A > 0 so that

(2.12)

IIf(a + ,~ Aa')ll ~< IIf(d)ll •

PROOF. Define g(d + A A d ) = Ilf(d + ,~ Ad)ll and g(d) = IIf(d)ll. It follows that we may use Lemma 2.1 if we can show that (g'(d), Ad) < 0 unconditionally. Using Newton's method to define Ad, we obtain

(g'(d), Ad) = ( f ( d ) , f ' ( d ) Ad) /IIf(d)ll = - I I f ( d ) l l < 0 , provided that f ~ 0 , and thus (2.12) is proved,

(2.13)

r-1

Equation (2.12) is an important result in that it proves that Newton's method can produce a norm-reducing solution regardless of the lack of symmetry of the tangent operator f ' ( d ) . In fact this result further shows that the tangent operator does not even have to be positivedefinite for (2.12) to hold, and thus any steplength algorithm based on Corollary 2.1 can be applied directly to nonlinear consolidation problems. Corollary 2.1 thus ensures that Ascan be ,ll÷l = d,,+l i + A~Ad ~, satisfies (2.12). On the other hand, if A~is chosen so that the next iterate, -,,+1 equal to unity, then (2.12) may fail and the Newton iterate is said to 'overshoot'. In this case, it is possible that the solution will never recover and the whole iteration process may diverge. The above 'damped' Newton-Raphson iteration then prevents such divergence to occur at all times, provided that the inverse of f ' ( d ) exists. The bound on the steplength A, A, is of importance in actual numerical analysis and can be determined, at least in principle, using the following lemma.

L E M M A 2.2. Assume that g: D C R"~q'->R ~ is G.differentiable at d E int(D) and that

IIg'(d + o A d ) - g'(d)ll <~ K011Adll VO~[0, 11 , where A d = - [ f ' ( d ) ] - l f ( d )

i - 211fll/gllAdll ~.

(2.14)

(Newton's method), then (2.12) holds for all A •[0, ,~], where

R. !, Borja, Composite Newton-PCG and quasi-Newton iterations

34

PROOF. We invoke the majorization principle [14] which states that if there exists a majorization function 0(A)" R t --> R t such that

g(d + A Ad) ~< 0(A) < g(d) VA • (0, ,~),

(2.15)

then d + A Ad decreases g. Using the mean-value theorem [24] and the Cauchy-Schwarz inequality [25], we have

g(d + A Ad) = g(d) +

fo'

A ( g'(d + 0A Ad), Ad) d0

- g(d) + A(g'(d), Ad) + ~

[ / g ' ( d + O~ Ad) - g'(d), a d ) l dO

fo'II a2llAdll'K fo'

<~g(d) + A(g'(d), Ad) + ;~llAdll <.g(d)+ A(g'(d),Ad) +

g'(d + OAAd) 0d0.

g'(d)ll dO (2.16)

Thus, choosing

¢(A)= g(d) + A(g'(d),Ad) +

½A2KIIAdII2

(2.17)

guarantees that

g(d + A Ad) ~< ¢(A).

(2.18)

To determine the majorization interval [0, ,~], we impose the constraint O(A)= g(d) + A(g'(d), &d) +

½A"/flIAdlI" < g(d),

(2.19)

which yields ,~ = - 2 ( g ' ( d ) , Ad)/KIIAdll ~ and since (g'(d), Ad)

--

-Ilfll

(2.20)

,

from (2,13), we obtain

X - 211fll/t,:lladll ~ •

(2.21)

Of course, one can choose, as the 'best estimate' for the steplength,

A,-- ,~/2 - I l f l l / g l l A d l l

~,

which gives the minimum of the majorization curve 0(A).

(2.22) [:]

The obvious difficulty with the formulation above lies in the fact that the coefficient K in (2.14) cannot be readily determined, Hence, in general, one cannot obtain 0(A) directly from

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

35

(2.17). However, if we assume a trial majorization curve 0(A) of the form (2.17), then a quadratic interpolation polynomial can be uniquely defined by forcing this trial curve to satisfy the following boundary conditions: (a) 0 ( 0 ) = g(d)= II$(d)ll- g,; (b) ~ ' ( 0 ) = (g'(d), Ad)= -IIf( )ll = - g , ; and (c) 0(A 2 > 0) = II$(d + A2 Ad)ll --- g2. The resulting polynomial is (cf, (2.17)) 0(A) = g . - g,A + ((g2 + g~A2- g,)/A~) A2-

(2.23)

To implement (2.23), we assume an initial trial value of ~2-" 1.0 and compare g2 with g~. If the steplength A 2 = 1 . 0 decreases g, then we accept the iterate and proceed with the next Newton-Raphson iteration; otherwise, we compute a new A2 which gives the minimum of 0(A) from (2.23), i.e., A2

*-g,A~/2(g2 + g,A 2

-

-

(2.24)

g,).

m

If qJ(A) is indeed a majorization curve, the next value of A2 will lead to a descent in g, by virtue of Lemma 2.2, and we proceed with the next Newton-Raphson iteration; otherwise, we continue searching for A that decreases g based on the idea derived from the majorization principle. A flowchart of the steps necessary to carry out these calculations is shown in Box 2. Note that the line search algorithm of Box 2 is very inexpensive. For example, Fig. l(a) shows that if A = 1.0 does not lead to a descent in g and the first trial quadratic interpolation polynomial is a majorization curve, then only one additional residual function evaluation will be required to determine the norm value at Ai, g(A,), which results in a decrease in g. Even if the first trial polynomial is not a majorization curve, Fig. l(b) shows that the line search algorithm can converge rapidly, and that no computed information is wasted as the previous value of A is now used to define the other end-point of the trial curve. However, this algorithm a MINIMUMOF TRIAL MAJORIZATIONCURVE o POINT ON RESIDUAL NORM CURVE

g

g().)

gi~) g(1)

g g(o)

g(O)

=4

~.~

(a)

I

(b)

Fig. 1. Line search algorithm for Newton's method: trial majorization curve finds A, after (a) one iteration~ (b) two iterations.

36

R. !. Borja, C o m p o s i t e N e w t o n - P C G

a n d q u a s i - N e w t o n iterattons

Box 2 Line search algorithm employing a quadratic interpolation polynomial for use with Newton-Raphson iteration 1. Initialize g, = IIf(d)ll, g2 = Ilf(d + Ad)ll, 2. I f g , < g ~ , return: else, continue 3.

A*--A'-gl/2(g,_ + gl A -

A - 1.0

gl)

4. g, = [[f(d + AAd)[I and go to Step 2.

can also cause the majorization interval to shrink rapidly causing the solution to converge to a point that does not necessarily represent the zero of g [14], since A in (2.12) can be made arbitrarily small. In this case, more elaborate steplength algorithms must be employed such as that based on Goldstein's principle [26] to ensure that the solution eventually converges to g=0.

3. Composite N e w t o n - P C G iteration

Let x = Ad z, b = f ,td , ,~, ~) and A = - f ' (d,,+ ' i ), and hence, the system of linear equations A x = b. Instead of solving these equations directly, say, by Crout elimination [9], which entails O(n:~. ) arithmetic operations dominated by triangular factorization, we now consider an iteratlve algorithm based on preconditioned conjugate gradients (PCG). Recall that the solution of A x = b at each iteration is required to carry out Newton's method (Box 1). In this approach, we thus end up with a composite Newton-PCG iteration, with Newton's method as the primary iteration and PCG as the secondary iteration. The preconditioned conjugate gradient technique is widely used for the minimization of quadratic functionals, although generalizations to nonquadratic functionals are also possible (e.g., see [6, 8, 27-29]). However, for nonquadratic functionals the cost associated with the one-dimensional minimization for the evaluation of the steplength parameter may require successive residual function evaluations which can become prohibitively expensive particularly in geomechanics where the constitutive models employed are generally more elaborate. We thus take full advantage of the benefits derived from Newton's method (such as quadratic convergence) and investigate in this section the utility of the composite Newton-PCG technique instead. The flowchart af Box 3 summarizes the steps necessary to carry out the solution of a linear system arising from Newton's iteration by the PCG technique. Box 3 Flowchart of the preconditioned conjugate gradient algorithm I. Initialize k--O, xt'--0, r~'=b, p O = z O _ B - i r , 2. c~=(rk, z~}l(p~,Ap ~) 4,

r k "~ = r x - or ~Ap ~

5, If [[r'~'[[/[[r°[[ <~TOL, 6, Zk''=B-Irk*'

returm

7. i3~= - ( z k " , A p k) /~pk Ap* ) 8.

p~~' = z k+' + ~"p"

9. k , - k + l

and go to Step2.

else continue

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

37

If the matrix A is symmetric and positive-definite, then the solution of A x = b is equivalent to minimizing the quadratic functional (3.1)

+ ½(x, A x ) .

However, for a non-symmetric matrix A the functional (3.1) does not exist. In order to preserve the essential features of the PCG technique for application to nonsymmetric linear systems, we attempt to reformulate the PCG algorithm based upon modified definitions of the PCG parameters a k and/3 k.

3.1. Derivation o f the P C G parameters a k and {3'

Let r k ' - b - A x * E R "~q be the residual vector associated with the secondary (PCG) iteration and pk E R'% the descent direction. If A is symmetric and positive definite, then the steplength c~k can be obtained by the minimization principle [14, 28]. If A is nonsymmetric, on the other hand, then a k can be determined by insisting that A x k÷ ~ = b, i.e., (3.2)

Ax* +1. = A(x* + a ' p * ) = b ,

which implies that a , = (r*, p* ) / ( p * , A p * ) .

(3.3)

It follows that r k+l = b -

A x k+l = b -

(3.4)

A(x* + akp *) = r* - a*Ap * .

Note that cek is unaffected by the symmetry, or lack of symmetry, of A. A direct consequence of (3.3) is given by the following identity. IDENTITY

3.1

( r k, p , - i } = 0 . PROOF ( r k, p k - l } = ( r , - , _ a k - l A p k - I ' p * - , ) = (r,-~

p k - ~ } _ [(rk-~, p,-~}/(pk-~, A p k - ~ ) ] ( A p k - ~

=0. I--I

Thus, we obtain the following alternative expression for a k'

pk-~ I

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

38

k = (An-

tr* + 13'- tp,- 1>/(p*,ap*> _ ( r k B - l r k ) / ( p k A p k > + ~ k - l ( r k ' pk-

1> /(pk,

Apk )

= /(pk, Apk>.

(3.5)

This is the form for k that appears in Step 2 of Box 3. In the classical conjugate gradient method the parameter flk is chosen so that pk and pk+ are mutually conjugate with respect to A in the sense that (pk, pk+~)A=(pkApk+~)= (Ap* pk +~) = 0, thus allowing the vector pJ to be generated simultaneously with the vector x j. For a nonsymmetric A, however, such an inner product does not exist since (p*,Ap *+~) # (Ap, p,+l). On the other hand, still insisting this type of conjugacy requirement in spite of the lack of symmetry of A will result in two possible expressions for/3~: (1) If we enforce (Ap k, pk+l) = 0 , then we obtain

(Apk, B-lr *+l + ~kp*> = 0 ,

(3.6)

from which

~*= -(B-'rk+',Ap *) /(Ap*, pk> . (2) If we require that

(3.7)

(p*,Ap *+~) =0, on the other hand, then we obtain

( p * , A ( B = ' r *+' +

~*p*)> =0,

(3.8)

from which

B*-- -(AB='r *+', P~) / (Ap*, p*) .

(3.9)

Between (3.7) and (3.9), we favor the former since it is simpler and does not require an extra matrix multiplication for the evaluation of the vector A l l - t r k+l (which the latter equation does). Thus, the form (3.7) is employed in Step 7 of Box 3. It is important to note that for symmetric systems, the expression for/3* usually takes on a different form. The more widely used expression for $ , generally results in simpler implementation than does either (3.7) or (3.9). To elaborate this point, assume that A is now symmetric and consider the following identity:

IDENTITY 3.2 (rk+',B-lr k) = 0 . PROOF (rk+',B-'rk) = (r k - (~kApk B-'rk ~ = (r k, B-,rk) _ [(r k, p , ) / ( p k , Apk)](Apk, B - t r k ) .

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

39

But

=
][~k-lpk-l>

= - t3~-~ = , since A = A'. Therefore,

._ _

= ( r k, B - Irk --

pk > = _ flk -,


l

> = O.

El It is natural to choose a symmetric preconditioner B if A is symmetric. For this case we obtain another identity.

I D E N T I T Y 3.3 ( B - l r k+l, r k} = 0 . PROOF (B-trk+i,r*> = (B-irk+l, Bp k - ~k-iBpk-' >

= (B-irk+i, Bpk> - flk-i(B-irk+l

Bpk-'>

= (rk+,, pk> _ flk-t = flk-i = 0 . El As a consequence of identities (3.2) and (3.3), we can derive an alternative form for O k as follows. Consider

(rk+l, pk+,) = (rk+t B-trk+l + ~kpk) = ( r ~ + l , B - t r k + ' ) + ~k(rk+t, p k ) = (rk+' B-irk+l}.

(3.10)

Now, consider

( pk+ ', rk +' > = ( pk+ ', rk -- akAp k ) = (pk°l, rk) -- ak(pk+t, Apk) = (B-'rk+' + ~kpk, rk} = (B-'rk+', r k) + flk (pk, rk} = flk ( B - ' r k + flk-tpk-', rk) = flk (B-Irk, rk} . (3.11) Since inner products commute, we can compare (3.10) with (3.11) and obtain the following alternative expression for ok: (3.12)

40

R. !. Borja, Composite Newton-PCG and quasi-Newton iterations

For symmetric systems, (3.12) has the advantage over either (3.7) or (3.9) that only r, and not A, is needed to form the next pk. This explains the wide usage of the form (3.12). For a nonsymmetric system, on the other hand, (3.12) tends to significantly slow down the convergence rate of the PCG solution because of its inherent assumption on the symmetry of A and B. In light of the fact that in practice we may be dealing with nonsymmetric systems as well, the form (3.7) for/3 k is used in lieu of (3.12) in Step 7 of Box 3. 3.2. Convergence and choice of preconditioner The fundamental result for the conjugate gradient algorithm when A is symmetric and positive definite is that the iterates converge to x ' " = A-~b in at most n~q steps. It must be emphasized that this important result is lost when one is dealing with a nonsymmetric system since the directions p". . . . . pm no longer form a conjugate basis with respect to A in the sense that =(Ap',pJ>=O,

i#j,

O<~i,j<~m,

(3.13)

for some m < n~q (for the conjugate gradient algorithm, a proof of (3.13) based on an inductive hypothesis when A is symmetric is given in [29]). To elaborate this point, note that (3.13) is central to the results of the following lemma which describes the convergence characteristics of the conjugate gradient algorithm. L E M M A 3.1. Let A E L(R '°~'a) be symmetric and positive definite and assume that (3.13) holds. Then, x'" satisfies x'" = A ~ tb for some m <<.n~q. PROOF. Write

(r *+,, p,) _= (r k, pJ) _ ak (Apk pi> = (r k, pi) _[
/ (AI

pk)l(APk ' It>,

assuming B = 1, the identity matrix. 1he A-orthogonality of the pJ from (3.13) implies that

(r k, p'>, {rk+l pi) = O,

j#k , j=k,

(3.14)

for j = 0 , . . . , n - 1. Since the pJ are linearly independent, we have r" ffi b - Ax"' = 0 for some m <~ n~,~. Note that this important result is not recovered when the matrix A is not symmetric since the inner product (3.13) does not exist for such system. El With finite precision machines the convergence properties of the conjugate gradient algorithm described in Lemma 3.1 may not be fully realized even when A is symmetric because of the accumulation of error. To reduce the accumulation of error and achieve faster convergence, it is usually necessary to introduce a preconditioning matrix B. For nonsymmetric systems the introduction of this preconditioning matrix plays an even more crucial role since these systems do not satisfy the conjugacy properties in the sense of (3.13). Hence, we seek an effective preconditioner which can be used not only for symmetric systems but for nonsymmetric systems as well.

R.i. Borja, Composite Newton-PCG and quasi-Newton iterations

41

In general, an effective preconditioning matrix B is one which satisfies the following criteria: (a) the equivalent system B - l a x = B - ~ b is easy to solve, and (b) the condition number of the iterative matrix B-~A is close to unity. These criteria can be satisfied by a matrix that is easy to factor (an alternative is to use a matrix that has already been factored), and which is a reasonably good approximation of A. The more closely B approximates A, the faster the convergence. In the context of nonlinear (elasto-plastic) consolidation analysis~ the tangent operator 2 A = - f F(d,,+~) generally varies not only from step to step (index = n) but also from iteration to iteration (index = i) as a result of the changing step size, progressive plastification and geometric mesh distortion. For example, the small-strain tangent operator for elasto-plastic consolidation varies according to (2.6). In a typical consolidation analysis the step size At can change by several orders of magnitude, and it is thus logical to choose a B that contains information on the present value of the step size so that the condition number of the matrix product B-~4 can be made as close to unity as possible. Although it is difficult to guess the influence of local yielding for purposes of deriving an even better approximation for A, we find the following preconditioner robust: ° B • = A o = - ~f ' ( d,-,,k+' ) = -f'(d,,k )'

(3.15)

where d~,l~+~-d,, k is the approximation to d at time t = t,,k, and t,,,, t , , . , . . . , t,,~ are time stations at which the step sizes are changed. We thus choose B as the re-formed (and re-factored) tangent eperator A during the first iteration of the time step at which the step size is changed. This preconditioner will generally contain the elastic component of the tangent operator N"'(Sq,l,+ ~) provided that homogeneous essential boundary conditions are prescribed (as non-zero prescribed displacements could cause initial plastification of elements near the movin$ boundaries during the first iteration). If the step size is changed every time step, then B is the tangent operator during the first iteration at each time step. In Section 5 we show that this preconditioner is extremely effective in dealing with symmetric and weakly asymmetric systems under a variety of loading conditions, in addition to the fact that this algorithm is very easy to code in practice. 3.1. Looking at the algorithm of Box 3, one finds the need to store five vectors of length n~q, the elements below the 'skyline' of the preconditioner B (in factored form), and the elements below the same 'skyline' of the tangent operator A. If storage space is unavailable, then one can either access B and A in turn (note that they do not have to reside in the main core at the same time), or re-form A every subiteration. The tradeoff is between I / O cost and CPU. In the examples of Section 5, both A and B were stored in the main core.

REMARK

4. Quasi.Newton iteration In this section we briefly describe an update procedure for the solution of nonlinear elasto-plastic consolidation equations. It is a quasi-Newton method in which the search direction Ad i is obtained from the expression

42

R.i. Borja, Composite N e w t o n - P C G and quasi-Newton iterations

Ad i = A/-t f(d,,+ i t) ,

(4.1)

where/[i satisfies the secant equation A iS = y ,

$ = d i,,+~ _

d ~l ,,+~,

y = f ( d , , +i -~1) - f ( n + , d) .i

(4.2)

The displacement update formula is then given by (4.3)

di,,+l ÷i = d , i, + l + Ai A d i

where At is a line search parameter. Conditions for the existence of Ai > 0 necessary to produce norm-reducing solutions in the context of nonlinear consolidation analyses are presented in Section 4.2. 4.1. B F G S iteration

We shall particularly focus our attention on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) inverse update [30-33], considered one of the best update formulas available [1], in which the inverse matrix A[ I is derived from its previous value, A~_t~, through updates or corrections t~f rank two:

t

A [ ' --" I

"') ( . t ) At_ z 1



+ (y,s)

'

(4.4)

where ! is the identity matrix. This formula can be phrased more elegantly in the factored form 13, 5] "

t

A-I

Af' --(! + w,v,)Af~,(l +

t

v,w,).

(4.5)

The vectors v~ and w~ are calculated from the formulas

v, - -< ( s. y ) / ( s. A,_,s) )""A,_,s - y

(4.6)

w,=s/(s,y)

(4.7)

and .

"-I Note that (4.5) is equivalent to (4.4) if and only if A~_~ is symmetric. Furthermore, since " -I i-I s = A~_~ A d ~-~ = Ai_~Ai_~f(d,,+~), we can also write (4.6) and (4.7) in the following form:

v, = - ( A i - , ( A d ' - ~ , Y ) / ( A d ' - I ,

ttdi-~,,x ,, ÷,)) )~'~f(d'-+ ~,, , ) - Y

(4.8)

and w, = Ad'-~/ ( A d '-~, y ) ,

so that vi and w/can be generated conveniently from vectors already computed.

(4.9)'

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

43

The form for/[~-~ in (4.4) and (4.5) allows the search direction Ad' to be computed without the costly factorization required in Newton's method. However, additional storage is required for storing historical information in the vectors v, and w~. Between (4.4) and (4.5), we choose the latter form because it is easier to implement, although it must be noted that this equation does not generally satisfy the secant equation (4.2) if A, is not symmetric, since, for this condition, (4.5) is not the same as (4.4). However, note that both (4.4) and (4.5) preserve the symmetry of the tangent operator, and if one starts the iteration with a symmetric A o = A 0 (such as the tangent operator which contains the symmetric elasticity matrix), neither form for / i will approach Newton's tangent operator for the non-symmetric case, so that there is really no strong argument about choosing one specific form for A over the other except to make the implementation simpler. As with the PCG technique, we shall re-form and re-factor the tangent operator every step size change. It may also be necessary to re-form and re-factor the tangent operator in instances where

(Ad'-', y) ( Adi-l, f(d,,+,)) i-1

> > 1.0

or

(Ad'-', y)

(Ad' ', f(d,,+,)) '-'

<0.

(4.10)

The first of (4.10) corresponds to a nearly singular matrix update due to a large condition number (in this paper the criteria used for avoiding numerically dangerous updates are the same as those discussed in [3]), while the second inequality gives rise to a square-root of a negative number in (4.8). Not performing the updating in the latter situation does not generally solve the problem since this condition also cautions one that a line search parameter A, necessary to decrease g might not exist (see Remark 4.2).

4.2. Extension of the simple line search algorithm to quasi-Newton iteration We now discuss the existence of a steplength parameter Ai>0 in (4.3) to force a quasi-Newton method to produce a norm-reducing solution. The basic idea is the same as that employed for Newton's method and can be summarized in the following lemma (again, omitting the time step and iteration counters for simplicity in notation).

LEMMA 4.1. Assume that f:l~"~"~---~l~'''', is G.differentiable at d and that A and ,4 are invertible, if Ad = ,4- If(d), where A satisfies the secant equation (4.2), then there exist A > 0 so that

IIf(d +

(4.11)

Ad)ll IIftd)ll,

provided that det(A)det(/[)>0. PROOF. Define g(d + A Ad)= IIf(d + A Ad)ll and g(d)= IIf(d)ll. It follows that we may use Lemma 2.1 if we can show that (g'(d), Ad) < 0. Using a quasi-Newton method to define Ad, we obtain

( g'(d), Ad) = - ( f(d), A Ad) / llf(d)ll - - ( f(d), Aft-tf(d) ) / llf(d)ll <0. since det (A) det(,4) > 0 implies that det(AA-t) > 0, and thus (4.11) is proved.

D

(4.12)

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

44

Now, g(d + A A d ) > 0 so that there is a value of A corresponding to a critical point of g(d + A Ad), if (4.12) holds. The minimum value of this A is the Curry steplength [34], which is the minimum A that satisfies the equation g'(A) = ( f ( d + A Ad),

f'(d + A A d ) A d ) / l l f ( d

+ a

ad)ll = o .

(4.13)

Equation (4.13) thus requires that the vectors f ( d + A Ad) and Ad be orthogonal with respect to the Jacobian f ' ( d + A Ad). An accurate determination of A from the orthogonality condition (4.13) requires expensive calculations and is best avoided in practice. However, in view of the results derived from the majorization principle of Lemma 2.2, we can also generate (as in Newton iteration) a quadratic interpolation polynomial 0(A) which approximates the majorization curve by forcing 0 to satisfy the following boundary conditions: (a) 0 ( 0 ) = g ( d ) = II/td)ll- g,; (b) 0'(o) = (g'(d), Ad) = - ( f ( d ) , A AM)/llf(d)ll = -go; and (c) ~(A 2 > 0 ) = IIf(d + a: ad)ll = g,,. The resulting parabola is O(A) = g, - g,,A + (g2 + goA2 - gt )/A ~ ) A : ,

(4.14)

which has its critical point at A~, = g,,A~/2(g 2 + g,,A..- g~).

(4.15)

We see that this algorithm, which is summarized in Box 4, is similar to the algorithm of Box 2. The graphical representation of this extended line search algorithm is similar to Fig. 1 except for the fact that the tangent line at A = 0 does not necessarily intersect the line g = 0 at A - 1.0. R E M A R K 4.1. The choice of g(A)= [[f(d + A Ad)[ I as the function to be minimized differs slightly from the conventional line search algorithm in which the function to be minimized is cq, whose gradient is f, assuming that such function exists. The relationship between q3 and g is

~d--{~d.R"~'~R'[V~--f}

so that g =

IIv ll.

(4.16)

It can bc seen that a critical point of qd occurs at V ~ = f = O, which gives the minimum (and the zero) of g. Now, for given d and Ad the conventional line search algorithm entails choosing a steplength A such that (see, e.g., [1]) ~d'(A) = ( f ( d + it Ad), Ad) = 0, Box 4

Line search algorithm employing a quadratic interpolation polynomial for use with quasi-Newton iteration I,

2. 3. 4. 5.

Initialize g , = II~d)ll, g=--II~d + ad)ll, A= 1.0 ifg,
(4.17)

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

45

i.e., A cancels the projection o f f ( d + A Ad) in the direction Ad. If the tangent operator A for Newton's method (or ,4 for quasi-Newton iteration) is positive definite, then the steplength parameter determined from (4.17) will accelerate the convergence of the iteration toward a local maximum of qd since (~'(d), Ad) = (f(d), A-~f(d)) > 0 and det(A)> 0. If A is negativedefinite, on the other hand, the convergence of the iteration will be accelerated toward a local minimum of qd. Since it is known that for consolidation problems A is generally indefinite and that qd may not exist after all, we find that the norm g = [Ill[ provides a more meaningful function to minimize for the line search, not to mention that this function is in fact used in practice as a measure of how well an iterative method is converging, hence, the algorithm of Box 4.

REMARK 4.2. Line searches require repeated evaluation of the residual vector corresponding to trial steplengths and can become very expensive when one uses very elaborate constitutive models such as those commonly employed in geomechanics. To reduce the cost associated with elaborate line searches, we employ the idea of Section 2 for Newton iteration and perform a line search only to guarantee that the BFGS iterates reduce g. As Lemma 4.1 indicates, the condition that ensures the existence of a steplength A, that decreases g is that det(A) dot(,4)> 0, i.e., the two determinants have the same sign. For small Ad this condition also avoids the possibility of the second inequality in (4.10) since

-

f(d/n+lt )-f(d~.÷l)) (Ad'-'. f(d.+,))'-t

(Ad '-~,

(4.18)

and if we expandf(d~,,+l) about d,.+ ,-i 1 using the Taylor series expansion and ignore the second order terms, we get

(Ad i-t, y) (Ad~O t, f(d~+t)) ~I

(Ad'-t, A,_ i Ad '-i ) >0, (Ad '~ 1, A,_t Ad t~ t )

(4.19)

provided that det(Al_ !) dot(,41_,) > 0. Now, when dot(A) dot(A) < 0 the BFGS approximation of A is considered ill.conditioned and we thus form and factor the tangent operator anew.

REMARK 4.3. As in the composite Newton-PCG technique, the BFGS iteration requires storage for A 0 for the evaluation of the search direction, A for the evaluation of the steplength, and 2naFGS vectors of size neq representing the v,'s and wi's, in addition to auxiliary vectors necessary to compute s and y. Since the nBFGs, the number of BFGS vectors, is linked to the number of iterations required for convergence, and is not known a priori, the BFGS technique will generally require more storage space than the composite Newton-PCG technique.

$. Numerical results

In addition to the full Newton method, the composite Newton-PCG and BFGS algorithms have been incorporated into the FE program SPIN2D [35], and test problems were solved to investigate the performance of each of these iterative solutions. Two constitutive models were

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

46

used in the study to simulate the deformation behavior of the soil skeleton: (1) the ideal Von Mises elasto-plastic model [36], and (2) the more realistic modified Cam-Clay model for soils [37-39]. Incremental stresses predicted by these models were computed using the radial return algorithm [22] and the fully implicit secant elastic stress integration algorithm described in [12]. The conditioning of the tangent operator depends crucially on the behavior of the material comprising the mesh, and so we solve problems involving not only homogeneous soils but layered soil materials as well. To further obtain a better impression of the performance of each of these iterative techniques, we also consider test problems which involve both symmetric and non-symmetric loading conditions, and then investigate the influence of mesh size on the CPU demand. Since the rate of convergence generally varies from technique to technique, we also compare the computing effort (in seconds CPU) required to satisfy (a) a nominal error tolerance, and (b) a stringent error tolerance. The convergence criterion is based on the residual norm and is given by

IIf*ll/Jlf'll

(5.1)

<

where a nominal error tolerance is set at eR = 10 -2 and a stringent error tolerance is set at e R = 10 -'~. The error tolerance for the PCG subiteration is TOL = 10 -2 (cf. Box 3). All computations were performed in (72-bits) double precision on an SC-30M DEC-20 workalike system called Macbeth (which is roughly 2.5 times as fast as a PDP-10 Model 2065) at Stanford University.

CENTERLINE WATER TABLE w

I A

BEDROCK(IMPERVIOUS) LOAD &

w~

~ .....

V ,To

SCALE 10 METERS • TIME

Fig. 2. Elasto-plastic consolidation: finite element mesh for the rectangular footing problem.

R.I, Borja, Composite Newton-PCG and quasi-Newton iterations

47

5.1. Rectangular footing on a Von Mises soil In this problem we assume that the soil is an elasto-plastic Von Mises material with the following material constants: Young's modulus E = 10~ kPa, Poisson's ratio v =0.0, uniaxial yield stress tr r = 5 0 0 k P a , and hardening parameter H ' =0.0; soil permeabilities are k~ = k 2 2 - - 0 . 0 1 m / d a y , kl2 = k E ! - - 0 . 0 (isotropic permeability, Darcy's law assumed valid). The consistent tangent operator A i is then symmetric in this example. The finite element mesh used in this example consists of one hundred D8P4 mixed elements (eight-node serendipity for displacements and 4-node element for pore pressure) all of which were assumed to be initially stress-free, and is shown in Fig. 2. This mesh generates a total of r/eq = 729 unknown degrees of freedom, a symmetric tangent operator containing 74,148 (single precision) words below the 'skyline', and a mean half-bandwidth of 101 words. A strip load of w - 103 kPa was then applied linearly with time over a period of T 0 = 1.0 day in one time step, and then the pore fluid was allowed to diffuse using 6 progressively increasing time increments chosen in such a way that TAF = temporal amplification factor = At,,+ ~/At,, = 2.0, i.e., the step size was doubled each time step. An effective method for solving the ordinary differential equation associated with this plane strain consolidation problem is provided by the two-step backward differentiation formula method (BDF-2) with a Crank-Nicolson start [17]. This method provides a combination of high frequency numerical dissipation and second-order accuracy not possible with any of the one-step generalized trapezoidal methods. In Fig. 3 we plot the ground settlement-time curve for point A at the center of the footing, and in Fig. 4 the variation of excess pore pressure at a centerline point B located 15 meters below the footing, along with a second set of solutions obtained from the same BDF-2 method employing a T A F = 1.5 to confirm the accuracy and stability of the time-stepping scheme. The MandeI-Cryer effect is evident in Fig. 4 (the initial

~

4

0

s

W

O TAF [] 1.5

10

1

I

I

10

100

1000

TIME,days Fig. 3. Consolidation of a Von Mises soil: settlement-time curve for node A by the BDF-2 method (TAF = temporal amplification factor).

R. !. Borja, Composite Newton-PCG and quasi-Newton iterations

48 m

0.8

O 0.6 Q

~

0.4

nO

0.2

X

LU

0.0

,

1

T A F = 2.0

I

!

10

100

1000

TIME, days

Fig. 4. Consolidation of a Von Mises soil: excess pore pressure ratio at B (#BIw) versus time by the BDF-2 method.

increase of excess pore pressure followed by the expected pore pressure dissipation), a characteristic feature of coupled phenomena which will also manifest in the remaining examples. Figure 5 shows the deformed mesh and the resulting zone of plastification 254 days after the full value of the footing load was applied (full plastification means all Gauss points yielded, partial plastification means at least one, but not all, Gauss point yielded). In Table 3 we summarize the performano: of (a) the full Newton method employing Crout triangular factorization for direct equation solving, (b) the composite Newton-PCO iteration, and (c) the BF(3S iteration, based on a global error tolerance e e = 10-". The required CPU CENTERLINE

~~11111111

/ i!iiii!ii !i~ill!iiiiiiill!/ii!i i!!i!ilil FULLPL~S~RO~ON ........ ~ PARTIAL PLASTIFIC&TION

SCALE: 10 METERS ! I

Fig, 5. Consolidation of a Von Mises soil: deformed mesh and zone of plastification.

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

49

Table 3 Consolidation of an elasto-plastic (Von Mises) soil subjected to a rectangular strip load by the BDF-2 method (e R = 1.0e - 2) Criteria

Newton-Raphson

No. load steps No. global iterations No. global iter/step No. factorizations No. subiterations No. subiter/global iter. No. residual eval. (LS) Total CPU, sec.

Composite NR-PCG

7 30 4.3 30 0 0 0 1,163

BFGS Iteration

7 30 4.3 7 256 11.1 0 1,020

7 78 11.1 13 0 0 24 1,305

are about the same for the first two methods showing that a mesh of this size, in conjunction with a symmetric tangent operator A~, is probably the 'break-even' point between the full Newton and the composite Newton-PCG iterations. Note that the PCG error tolerance of TOL = 10 -2 did not affect the global convergence rate of the Newton-PCG iteration. Of the three techniques the BFGS iteration took the longest time to complete the solution.

5.2. Triangular footing on a Von Mises soil We now consider a mesh twice as large as that of Example 5.1 and apply a linearly distributed (triangular) strip load as shown in Fig. 6. The mesh of Fig. 6 generates a total of neq = 1,449 unknown degrees of freedom, a symmetric tangent operator containing 232,332 words below the 'skyline', and a mean half-bandwidth of 160 words. The maximum intensity of the load is w = 1.5 x 103 kPa applied linearly with time over a period of T~ - 1.0 day on a MESH I CENTERLINE I

14) WATER TABLE

I A SOiL 1

is

SOIL 2

BEDROCK (IMPERVIOUS)

SCALE 10 METERS I t

Fig. 6. Elasto-plastic consolidation: finite element mesh for the triangular footing problem.

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

50

homogeneous Von Mises soil (i.e., soil #1 is the same as soil #2) possessing the same characteristics as the soil of the preceding example. In Figs. 7 and 8 we show the displacement-time history of a surface node A on the mesh centerline (note the leftward-rightward lateral movement of the point during consolidation) and the excess pore pressure-time history for node B situated 7 meters below nodal point A, again employing the same two temporal amplification factors, TAF = 2.0 and TAF = 1.5, to demonstrate the accuracy and stability of 0.0

ii

JT=0.0

4) LOADING

E ).: Z

T=1.0

-1.0

uJ

:S uJ (J ,,J n_ (n o _1

CONSOLIDATION

-2.0 r-I TAF = 2.0

2?

O T A F = 1.5

r~ LU =) • 3.0

,

,

• 0.1

,

0,0

,

,

0.1

),2

HORIZONTAL DISPLACEMENT, meters

Fig. 7. Consolidation of a Von Mises soil: displacement of node A versus time (in days) by the BDF-2 method. 0,8

- -

m

0.8 D TAF : 2,0 ,5

0,4

0,2

W

0,0

I

I

10

100

1000

TIME, days

Fig, 8, Consolidation of a Von Mises soil: excess pore pressure ratio at B (Oa/w) versus time by the BDF-2 method.

51

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

the BDF-2 solutions. Figure 9 shows the resulting deformed mesh and zone of plastification 254 days after the soil was allowed to consolidate. Table 4 summarizes the performance of each of the three iterative techniques considered and shows that the composite N e w t o n - P C G and the B F G S solutions required about 40% less CPU compared with the full Newton method. Note that the number of subiterations per global iteration required by the composite N e w t o n - P C G technique is about the same for this mesh as for the (smaller) mesh of Example 5.1 (approximately 11 subiterations per global iteration, average). A striking result and interesting feature of this composite N e w t o n - P C G solution is that it required fewer global iterations to converge than the full Newton method, suggesting that an exact linear equation solving (by direct method) is in fact not essential when one deals with a nominal convergence tolerance (e.g., e R = 10-2). We also see that the number of global iterations required by the BFGS technique remained essentially the same as W

I:¸,:,¸%

~- .... D

_

'

i: !:!~ | ! !!!

: i :if:

I

BEDROCK (IMPERVIOUS) ~ m

FULL PI.ASTIFICATION

....

PARTIAl. PIºASTIFICATION

SCAI,E: 10 METERS I I

Fig. 9. Consolidation of a Von Mises soil: deformed mesh and zone of plastification.

Table 4 Consolidation of an elasto-plastic (Von Mises) soil subjected to a triangular strip load by the BDF-2 method (e R = 1.0e - 2)

Criteria No. load steps No. global iterations No. global iter / step No. factorizations No. subiterations No. subiter/global iter. No. residual eval. (LS) Total CPU, sec.

Newton-Raphson 7 33 4.7 33 0 0 0 5,528

CompositeNR-PCG 7 30 4.3 7 260 11.3 0 F 3,464 ..

BFGS Iteration

,,,

7 79 11.3 12 0 0 4 3,742

52

R. !. Borja, Composite Newton-PCG and quasi-Newton iterations

that of Example 5.1 (approximately 11 iterations per step, average), which suggests considerable potential of both the composite Newton-PCG and BFGS techniques for symmetric large-scale finite element equation solving particularly with application to nonlinear consolidation problems. 5.3. Consolidation of a homogeneous Undisturbed San Francisco Bay Mud

As a third example, consider the same FE mesh of Fig. 6 and assume that the soil is Bay Mud whose Cam-Clay parameters are [40]: virgin compression index A = 0,37, swell-recompression index K = 0.054, slope of the critical state line M = 1.40 and reference void ratio e, = 2.52; permeability: k~ = k22 = 0.02 m / d a y , kl2 = k21 =0.0 (isotropic permeability, Darcy's law holds). In this example, the version of the modified Cam-Clay model used to represent the behavior of the Bay Mud is described in [37-39]. The numerical integration algorithm employed in the framework of associative plasticity is discussed in [12]. It must be noted that the stress integration algorithm of [12] gives rise to a weakly asymmetric consistent tangent operator which must be considered accordingly if one is to achieve the (optimum) quadratic convergence rate in Newton's method. The first step in the analysis is to establish the initial stress condition using the gravity load imposition procedure described in [17]. This part of the analysis is necessary since the Cam-Clay model treats the elastic moduli as stress-dependent quantities. A triangular strip load was then applied in five equal time increments using a step size of At = 2.0 days during which time the soil was also allowed to consolidate. The transient analysis was carried out using the constant step size three-step backward differentiation formala method (BDF-3) with a Crank-Nicolson/BDF-2 start [17]. With sufficiently accurate initial starting values this method produces third-order accurate, unconditionally stable solutions and is capable of effectively damping high frequency modes [17, 41]. The accuracy of the BDF-3 method is

1001..0

.o.:,.

o

~ -lO0

.,oo

--%

;: .,oo -500

t 0

,

* 10

,

I

,

20

30

HORIZONTAL DISPLACEMENT OF A, mm. Fig, 10, Consolidation of Undisturbed Bay Mud: displacement of node A versus time by the BDF-3 method.

R.I. Borja, Composite Newton-PC(J and quasi-Newton iterations

53

demonstrated in Figs. 10 and 11 which show nearly identical displacement-time history for node A, and nearly the same excess pore pressure profiles along the mesh centerline obtained from five and ten time-step analyses. Figure 12 shows the resulting deformed mesh and zone of plastification, the latter represented by the shaded area within which overconsolidation ratio (OCR) is equal to unity (for three-dimensional applications OCR is defined in [38]). Comparing the final yield zone with the initial yield zone also shown in Fig. 12, note that the soil was unloaded on both sides of the strip load during this period. 10 ORAINAGE

-10

m !,._

E

AT = 2.0 AT = 1.0

-20

:If ~ -30 W 0 -40

-50

~

0.0

r: .o 0.1

r: oo 0.2

0.3

0.4

0.5

0.6

EXCESS PORE PRESSURE RATIO Fig. 11. Profile of excess pore pressure ratio #~,/w along the mesh centcrline by the BDF-3 method (T in days). le9

~////////,~/I////~./.F~

I

FULL PLASTIFICATION PARTIAL PLASTIFICATION

;~ / , ~ ' , F / / / ' / . / / ~ / / / / / : MESH SCALE: 10 METERS I I

DISPL. SCALE: 2 METERS I I

Fig. 12. Consolidation of homogeneous Bay Mud: deformed mesh and zone of plastification.

54

R.I. Borja, Conlposite Newton-PCG and quasi-Newton iterations

Table 5 Consolidation of a homogeneous Undisturbed Bay Mud subjected to a triangular strip load by the BDF-3 method (e~ = l.Oe- 2) Criteria No. load steps" No. global iterations No. global iter/step No. factorizations No. subiterations No. subiter/global iter. No. residual eval. (LS) Total CPU, see.

Newton-Raphson 6 24 4.0 24 0 0 0 7,820

Composite N R - P C G 6 24 4.0 2 170 7.7 0 2,530

BFGS Iteration 6 79 13.2 11 0 0 73 7,255

" Includes the gravity load imposition step.

In Table 5 we summarize the performance of each of the three techniques and see that the composite Newton-PCG procedure again resulted in significant savings (about 70% less CPU) compared with either the full Newton or the BFGS techniques. The BFGS technique required large number of costly re-formation and re-factorization of the global tangent operator as a result of ill-conditioning during the iteration process (see Section 4.2). Consequently, the CPU demand of this technique was high, and the performance is hence poor.

5.4. Consolidation of layered porous media We now consider a more realistic problem involving two soils having different material characteristics: a 15-meter thick granular landfill (soil #1 in Fig. 6) underlain by a 35-meter thick Bay Mud (soil #2 in the same figure), which is a typical soil profile in the San Francisco Bay Area. The Bay Mud considered is the same as that of Example 5.3 while the granular landfill, also modeled as a Cam-Clay material, is assumed to have the following parameters [42]: A=0.025, •--0.005, M--1.50, e,,=2.52; permeabilities are: ktt=k,,=2.0m/day, kt~ = k~ =0.0 (isotropic permeability, Darcy's law valid). Note that the landfill is more permeable, but less compressible, than the Bay Mud. We again apply the triangular strip load in (a) 5 equal increments and (b) 10 equal increments over a period of T,~ = 10.0days. Figures 13 and 14 show the respective timevarying displacement of node A and excess pore pressure profile along the centerline obtained using the BDF-3 method. For comparison, Fig. 13 also shows the theoretical movement of the same point A if drainage of the soil were fully suppressed (undrained loading), thus demonstrating the significant consolidation-induced ground movement occurring during the loading stage. The deformed mesh and zone of plastification are shown in Fig. 15 which indicates that movement of the foundation is smaller, although plastification is more extensive than those of the preceding example. In Table 6 we evaluate the performance of the same three iterative techniques against the backdrop of a poorly-conditioned tangent operator A~ and see that the composite NewtonPCG iteration again yields the best result (it is 64% cheaper than the full Newton method and 57% cheaper than the BFGS technique). Considering that the global preconditioner B defined in (3.15) keeps the average number of subiterations required for iterative linear equation

R.i. Borja, Composite Newton-PCG and quasi-Newton iterations

55

10 DRAINAGE

0 days / CONSOLIDATION: O 5-STEPSOLUTION 10-STEP SOLUTION

E

~-20 ILl I¢I

0 [] -

u) -10

-

AT=2.0days AT=l.0day

~ LANDFILL

G) G) E -20 I,,.

:g

0 :S .4o 0

~-60

I--

6.0 UNDRAINED IIOA,m ,~ 0 ~ [] 5-STEP SOLUTION ---- 10-STEPSOLUTION

Q. UJ Q

-30 -40

8.0~ ~)

:>

-50 T=10.0 days

-8O

I

0

I

1

I

2

I

3

I

4

-60 0.0

I

,

,

0.1 0.2 EXCESS PORE PRESSURE R A T I O

H O R I Z O N T A L M O V E M E N T , ram.

Fig. 13. Consolidation of Bay Mud-landfill: displacement of node A versus time by the BDF-3 method.

I

,

0.3

Fig. 14. Profile of excess pore pressure ratio 01~/w along the mesh centerline by the BDF-3 method.

w

FULL PLASTIFICATION PARTIAL PLASTIFICATION

MESHSCALE: 10METERS I

,I

DISPL.SCALE: 200MM I

I

Fig. 15. Consolidation of layered media (Bay Mud-landfill)' deformed mesh and zone of plastification.

solving by the PCG technique constant with respect to the mesh size, we can conclude from this example that for weakly asymmetric systems the Newton-PCG technique also possesses considerable potential for large-scale FE equation-solving. In Table 7 we re-analyze the same layered consolidation problem and prescribe an error tolerance of e R = 10 -~, with the same PCG subiteration tolerance of TOL = 10-". Compared to

R.i. Borja, Composite Newton-PCG and quasi-Newton iterations

56

Table 6 Consolidation of a layered elasto-plastic (modified Cam-Clay) soil subjected to a triangular strip load by the BDF-3 method (e R = 1 . 0 e - 2)

Criteria No. No. No. No. No. No. No.

load steps '~ global iterations global iter/step factorizations subiterations

subiter/global iter. residual eval. (LS) Total CPU. sec.

Newton-Raphson 6 22 3.7 22 0 0 0 7,122

Composite N R - P C G 6 22 3.7 2 165 8.3 0 2,533

BFGS Iteration 6 75 12.5 10 0 0 29 5,881

" Includes the gravity load imposition step. Table 7 Consolidation of a layered elasto-plastic (modified Cam-Clay) soil subjected to a triangular strip load by the BDF-3 method (e, = 1 . 0 e - 5) Criteria No. No. No. No. No. No. No.

load steps" global iterations global iter/step factorizations subiterations

subiter/global iter, residual evai, (LS) Total CPU, sec,

Newton-Raphson 6 27 4.5 27 0 0 0 8,763

Composite N R - P C G 6 31 5.1 2 239 8,2 0 3,365

BFGS Iteration 6 134 22.3 10 0 0 49 8,622

"Includes the gravity load imposition step,

the full Newton method, the composite Newton-PCG technique generally required one extra global iteration per time step to achieve this stringent convergence requirement, although it still ended up 62% cheaper than the full Newton method due to small number of costly factorizations. The slow global convergence of the BFGS iteration makes it less competitive than any of the two Newton-based iterations when the accuracy demand is high, as Table 8 shows. On another note, one can further improve the performance of the Newton-PCG technique by starting with a loose convergence requirement associated with the secondary (sub)iteration, and then gradually tightening this up as the global error associated with the primary iteration diminishes, as suggested in [8]. This will result in more elaborate error control which may not be warranted after all, considering the already superior performance of the composite Newton-PCG technique in its present form. REMARK 5.1. As can be gleaned from Tables 3-8, Newton's method does not generally require a line search as often as, say, the BFGS technique. However, if it does and a line search routine is unavailable, the solution will usually diverge and one may not finish an otherwise saveable solution. This is exemplified in Table 9 which demonstrates the usefulness of the simple line search algorithm described in this paper. The problem involves the mesh of

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

57

Table 8 R e l a t i v e residual n o r m values for step n u m b e r 2 of the layered elasto-plastic soil consolidation p r o b l e m (e R = 1 . 0 e - 5) I t e r a t i o n no.

Newton-Raphson

Composite NR-PCG

1

1.O00e + 0

1.O00e + 0

l.O00e + 0

2 3 4 5 6 7

1.615e1.735e 1.946e 1.445e --

1.599e 1.767e 1.892e 1.476e 9.653e --

6.880e 5.404e 1.785e 1.110e 5.836e 3.565e

1 2 3 7

-

BFGS Iteration

1 2 3 5 8

22 23

-

1 1 1 1 2 2

1.514e - 5 8.517e - 6

Table 9 R e l a t i v e residual n o r m values for step n u m b e r 1: elasto-plastic c o n s o l i d a t i o n o f Bay M u d s u b j e c t e d to a r e c t a n g u l a r strip load by the N e w t o n - P C G t e c h n i q u e with a n d w i t h o u t line search (E R = 1.0e - 5)

Iteration no.

N o line search

Line search

LS p a r a m e t e r ,~,

1

1.000e + 0

1.000e + 0

1.00

2 3 4 5

1.164e + 0 1.428e + 0 diverged

6.606e 7.970e 1.738e 1.544e

1 2 3 5

0.42 1.00 i.00 1.00

1.051 e - 7

1.00

6

-

Fig. 2 and assumes that the soil is a Bay Mud subjected to a rectangular strip load of w = 30 kPa. The solution was carried out by the composite Newton-PCG technique, with and without a line search, using a single time increment of At = 2.0 days. Without a line search, Table 9 indicates that the solution diverged. This type of 'numerical breakdown' is typical in analyses employing critical state models [11, 12]. Also observe that the quadratic rate of convergence in Newton's method was lost when an approximate (PCG-based) linear equation solver with an error tolerance of TOL = 10-: was employed, and that, like the convergence rate of the Newton-PCG technique shown in Table 8, the maximum rate of dissipation of the residual norm appears to be governed by this same error tolerance TOL. R E M A R K 5.2. Finally, we show in Table 10 a comparison between two Newton-PCG solutions, one employing (3.7) for the determination of/3 k and another employing the more widely used expression (3.12). We see that form (3.7) required fewer subiterations to satisfy a convergence tolerance of TOL = l 0 - 2 than (3.12), and thus, (3.7) was used in Step 7 of Box 3.

58

R. !. Borja, Composite Newton-PCG and quasi-Newton iterations

Table 10 Convergence of PCG iterations employing two expressions for/3 k for the solution of weakly asymmetric system of linear equations

Criteria No. load steps ~' No. global iterations No. global iter/step No. factorizations No. subiterations No. subiter/global iter. No. residual eval. (LS) Total CPU, sec.

Eq. (3.7) 6 31 5.1 2 239 8.2 0 3,365

Eq. (3.12) 6 31 5.1 2 285 9.8 0 3,660

" includes the gravity load imposition step.

6. Closure We have pointed out the advantages and limitations of Newton's method for the solution of nonlinear consolidation problems. For the linearized problem the task of triangular factorization when the system is large can be immense, and thus we sought alternative equation solving techniques which can perform the same task at a more reasonable cost. While our account and investigation are by no means complete, we find that for symmetric and weakly asymmetric systems the composite Newton-PCG technique possesses considerable potential for usefulness in large-scale computations. The performance of the BFGS technique is not quite as encouraging for nonsymmetric systems as it is for symmetric systems. Finally, we note that all of the numerical techniques we have described are merely a few of the most promising methods selected from a long list of nonlinear equation solving techniques that could also be applied to nonlinear consolidation problems.

Acknowledgment The author is grateful to anonymous reviewers for their constructive reviews and comments. Financial support for this research was provided in part by the National Science Foundation under Contract No. MSS-8910219, Research Initiation Award.

References

[11

J.E. Dennis, Jr. and J,J. Mord, Quasi-Newton methods, motivation and theory, SIAM Rev. 19 (1) (1977) 46-89. [2] J.M. Winget and T.J.R. Hughes, Solution algorithms for nonlinear transient heat conduction analysis employing element-by-element iterative strategies, Comput, Methods Appl. Mech. Engrg. 52 (1985) 711-815. [3] K.J. Bathe and A,E Cimento, Some practical procedures for the solution of nonlinear finite element equations, Comput, Methods Appl, Mech. Engrg, 22 (1980) 59-85,

R.I. Borja, Composite Newton-PCG and quasi-Newton iterations

59

[4] T.J.R. Hughes, I. Levit and J. Winget, An element-by-element solution algorithm for problems of structural and solid mechanics, Comput. Methods Appl. Mech. Engrg. 36 (1983) 241-254. [5l H. Matthies and G. Strang, The solution of nonlinear finite element equations, Internat. J. Numer. Methods Engrg. 14 (1979) 1613-1626. [61 A.K. Noor and J.M. Peters, Preconditioned conjugate gradient technique for the analysis of symmetric anisotropic structures, Internat. J. Numer. Methods Engrg. 24 (1987) 2057-2070. [7] M. Ortiz, P.M. Pinsky and R.L. Taylor, Unconditionally stable element-by-element algorithms for dynamic problems, Comput. Methods Appl. Mech. Engrg. 36 (1983)223-239. [81 M, Papadrakakis and C.J. Gantes, Preconditioned conjugate- and secant-Newton methods for non-linear problems, Internat. J. Numer. Methods Engrg. 28 (1989) 1299-1316. [9] R.L. Taylor, Computer procedures for finite element analysis, in: O.C. Zienkiewicz, ed., The Finite Element Method, 3rd Ed. (McGraw-Hill, New York, 1977) Chapter 24. 1101 R.I. Borja, The analysis of consolidation by a quasi-Newton technique, Internat. J. Numer. Anal. Methods Geomech. 12 (1988) 221-229. 1111 R.I. Borja and S.R. Lee, Cam-Clay plasticity, Part l: Implicit integration of elasto-plastic constitutive relations, Comput. Methods Appl. Mech. Engrg. 78 (1990) 49-72. [12] R.I. Borja, Cam-Clay plasticity, Part II: Implicit integration of constitutive equation based on a nonlinear elastic stress predictor, Comput. Methods Appl. Mech. Engrg. (1991) to appear. [13] M. Ortiz and J.B. Martin, Symmetry-preserving return mapping algorithms and incrementally extremal paths: A unification of concepts, lnternat. J. Numer. Methods Engrg. 28 (1989) 1839-1853. [141 J.M. Ortega and W,C. Rheinboldt, lterative Solution of Nonlinear Equations in Several Variables (Academic Press, New York, 1970). [lSl M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (1941) 155-164. ll61 R.I. Borja, Linearization of elasto-plastic consolidation equations, Engrg. Comput. 6 (2) (1989) 163-168. 1171 R.I. Borja, One-step and linear multistep methods for non-linear consolidation, Comput. Methods Appi. Mech. Engrg. 85 (1991)239-272. !181 C.W. Gear, B. Leimkuhler and G.K. Gupta, Automatic integration of Euler-Lagrange equations with constraints, J. Comput. Appi. Math. 12, 13 (1985) 77-90. II91 P. L(itstedt and L.R. Petzold, Numerical solution of nonlinear differential equations with algebraic constraints l: Convergence results for backward differentiation fonnuhis, Math. Comp, 46 (174) (198t~) 491-51h. izoi D.R. Owen and E, Hinton, Finite Elements in Plasticity: Theory and Practice (McGraw-Hill, New York, 1980). I211 J.C, Simo and R.L. Taylor, Consistent tangent operators for rate-independent ehlstoplasticity. Comput. Methods Appl. Mech. Engrg. 48 (1985) 101-118. 1221 M.L. Wilkins, Calculation of elasto-plastic flow, Methods of Computational Physics 3 (Academic Press. New York, 1964). [231 J,C, Simo, J,G. Kennedy and S. Govindjee, Non-smooth multisurface plasticity and viscoplasticity. Loading/ unloading ~onditions and numerical algorithms, lnternat. J. Numer. Methods Engrg. 26 (1988) 2161-2185. [241 J. Dieudonn6, Foundations of Modern Analysis (Academic Press, New York, 1969). [251 J.L. Nowinski, Applications of Functional Analysis in Engineering (Plenum, New York, 1981). [261 A. Goldstein, Minimizing functionals on normed linear spaces, SIAM J. Control 4 (1966) 81-89. 1271 J. Daniel, The conjugate gradient method for linear and nonlinear operator equations, SIAM J. Numer. Anal. 4 (1967) 10-26. 1281 R. Fletcher and C. Rezves, Function minimization by conjugate gradients, Comput. J. 7 (1964) 149-154. 129] M. Hestenes, The conjugate gradient method for solving linear systems, Proc. Sixth Syrup. Appl. Math., Amer. Mathematical See., Providence, RI (1956)83-102. 13o1 C.G. Broyden, A new double-r'~nk minimization algorithm, Notices Amer. Math. Soc. 16 (1969) 670. 1311 R. Fletcher, A new approach to variable metric algorithms, Comput. J. 13 (1970) 317-322. [321 D. Goldfarb, A family of variable-metric methods derived by variational means, Math. Comp. 24 (1970J 23-26. 1331 D.F. Shanno, Conditioning of quasi-Newton methods for function minimization, Math. Comp. 24 (1970) 647-656.

60

R. !. Borja, Composite Newton-PCG and quasi-Newton iterations

[341 H. Curry, The method of steepest descent for nonlinear minimization problems, Quart. Appl, Math. 2 (1944) 258-261. 1351 R.I. Borja and E. Kavazanjian, Jr., Finite element analysis of time-dependent behavior of soft clays, Geotech. Engrg. Res. Report No. GT1, Stanford University, Stanford, CA 94305, USA, 1984. 1361 R. Hill, The Mathematical Theory of Plasticity (Oxford Univ. Press, London, 1967). [371 J.H. Atkinson, Foundations and Slopes: An Introduction to Applications of Critical State Soil Mechanics (Halsted Press, New York, 1981). 1381 R.I. Borja and E. Kavazanjian, Jr., A constitutive model for the stress-strain-time behaviour of 'wet' clays, Gdotechnique 35 (3) (1985) 283-298. 1391 K.H. Roscoe and J.B. Burland, On the generalized stress-strain behaviour of 'wet' clays, in: J. Heyman and F.A. Leckie, eds., Engineering Plasticity (Cambridge Univ. Press, Cambridge, 1968) 535-609. [401 R. Bonaparte and J.K. Mitchell, The properties of San Francisco Bay Mud at Hamilton Air Force Base, California, Geotech. Engrg. Res. Repolt, University of California, Berkeley, 1979. [41] C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice-Hail, Englewood Cliffs, NJ, 1971). [42] T.W. Lambe and R.V. Whitman, Soil Mechanics (Wiley, New York, 1969).