INTERNATIONALJOURNAL FOR NUMERICAL METHODS I N ENGINEERING,VOL.39,2635-2646 (1996) A GENERALIZED NEWTON METHOD FOR HIGHER-ORDER FINITE ELEMENT APPROXIMATIONS I N NON-LINEAR ELASTICITY PANAYIOTIS PAPADOPOULOS Department of Mechanical Engineering, University of California at Berkeley, Berkeley, C A 94720-1 740, U.S.A. ROBERT L. TAYLOR Department of Civil Engineering, University of Calijornia at Berkeley, Berkeley, C A 94720-1 710, U.S.A. SUMMARY A generalized Newton method is proposed in conjunction with a higher-order Lagrangian finite element discretization of bodies undergoing finite elastic deformations. The method is based on a gradient-like modification of the Newton method, designed to suppress the sensitivity of higher-order elements during the early iterations, thus allowing for solutions to be obtained using moderately large step-sizes. KEY WORDS non-linear elasticity; finite element method; Newton method 1. INTRODUCTION Parametric Lagrangian finite element formulations are used extensively in the numerical simulation of the mechanical response of elastic bodies that may undergo finite deformations. A wellestablished limitation of such formulations stems from the inability of Lagrangian elements to stretch and rotate beyond limits set by the assumed invertibility of their motion. Higher-order elements, such as the 8- and 9-noded quadrilaterals in two dimensions or the 20- and 27-noded hexahedra in three dimensions, although very desirable for their accuracy, are substantially more sensitive to distortions than their respective low-order counterparts. This shortcoming is clearly evident when the discrete equations of motion are solved iteratively using the classical Newton method, even at relatively small strains. Indeed, the step-size of an incremental analysis may have to be severely restricted in order for a solution path to be successfully traced. The present article advocates a generalized Newton method which can yield convergent Lagrangian solutions with appreciably larger step-sizes than the classical Newton method. The Levenberg-Marquardt algorithm,' originally suggested for the computation of non-linear least-squares parameters, provides the motivation for the proposed method. The main idea is to appropriately combine the classical Newton method with a variable parameter gradient-based (or steepest-descent) method, so as to obtain an iterative solution scheme which selectively suppresses the sensitivity associated with higher-order elements, especially during the early iterations within each time step. Elements of finite elasticity are presented in Section 2 by way of background to the ensuing developments. The proposed generalized Newton method is introduced and analysed in Section 3, while the results of two numerical paradigms using this method are presented in Section4. Concluding remarks are given in Section 5. 9 CCC OO29-598 1/96/152635-12 0 1996 by John Wiley & Sons, Ltd. Received 18 M a y 1995 Revised 23 February 1996 2636 P. PAPADOPOULOS AND R. L. TAYLOR 2. A BRIEF O N FINITE ELASTICITY Consider a deformable body 98 identified with an open set R in linear space R3, equipped with a fixed orthonormal basis (El, E2,E3). A material point of 98 in the reference configuration is algebraically specified by a vector X. The boundary of 93 is represented by aR and possesses a unique outer unit normal N at each of its points. The motion of the body is described by the smooth and invertible mapping x on so that the position vector x of a material point in the current configuration R, is given by x = x(X, t ) at each time t, and the respective displacement vector u is defined according to u(X, t):= x(X,t ) - X. Invertibility of the motion at all times implies that J := det F > 0, where the deformation gradient F is defined as a, Neglecting the effect of inertia, the equations of motion are written in the reference configuration as Div P + pob = 0 in R u=U PN = p on Tu on (1) r4 where Dirichlet and Neumann boundary conditions are enforced on ru# 8 and Ts, respectively, and T,wT, = an. In the above equations, P is the first Piola-Kirchhoff stress tensor, po is the mass density in the reference configuration, b is the body force per unit mass, ii is the prescribed boundary displacement on I-,,, and p is the prescribed traction vector on r4 measured in the reference configuration. Application of a standard weighted-residual method results in a weak form associated with such that equations (l), which determines a solution u E @ := {u E H'(R) 1 u = U on r,,}, A[w,u]:= I Gradw:PdV - I* pow.bdV - 6, w-pdr=0 (2) for all w ~ 7 ( l r : = {wEH'(R)Iw = 0 on r,,}. Within the context of classical hyper-elasticity, admit the existence of a strain energy function Y = Y (X,F) per unit volume, so that, following standard arguments, the first Piola-Kirchhoff stress tensor is given by ay p=aF (3) The motion of body 99 is assumed to be infinitesimally super-stable for the given boundary data. in the sense that jnGradw:-GradwdV aP dF >0 (4) for all non-trivial infinitesimal w E W ,see the review article of Beatty (Reference 3, Section 19).As concluded in Reference 4, inequality (4) is a sufficient condition for uniqueness of an incremental boundary-value problem, where a solution is sought starting from a known equilibrium state, and imposing incremental changes to the boundary conditions (lb, c). HIGHER-ORDER FINITE ELEMENT APPROXIMATIONS 2637 It should be noted that condition (4) is generally very restrictive and should not necessarily hold for all elastic motions. In the context of this work, condition (4) is introduced only for the purpose of proving certain convergence properties of the proposed solution method. 3. A GENERALIZED NEWTON METHOD An approximate solution of (2) under constitutive assumption (3) and (when appropriate) stability condition (4) can be obtained by admitting a spatial (finite-element-based) discretization in order to construct finite-dimensional subspaces h%! and wh of 92 and *w, respectively. Subsequently, U h E %h is determined so that A[Wh. uh] (5) =0 for all wh E wh. Let a piecewise polynomial interpolation of the displacement field over the domain of analysis be written as N L U h ( x , f) = ih(X, f) + ih(X, t ) = 1 @I(X)dI(t)+ 1 @I(x)aI(t) I= 1 I=L+1 := N(X)d(t) + A(X)a(t):= N(X)d(t) (6) where degrees of freedom &, I = L + 1, . . . ,N, are added hierarchically to an initially assumed low-order polynomial approximation 6, (Reference 5, p. 139).In the above equation, OIdenote the standard and hierarchic interpolation functions for 1 6 I < L and L + 1 < I 5 PJ, respectively. Consequently, the global displacement vector dE R3N is partitioned as d = [d a]', where 2 E R3L and d E R 3 ( N - L )are the standard and hierarchic displacement sub-vectors, respectively. Similarly, write N L + := N ( X ) i ( t ) A(X)i(f):= N(X)v(t) (7) where i, and i, are arbitrary vector functions that vanish identically where Dirichlet boundary conditions are enforced, and V E R3N is the global weighting vector. It follows from (7) that L Grad wh = + N Grad(Qr(X)i1(t)) I= 1 Grad(@I(X)tr(t)) I=L+1 := B(X)i ( t ) + B(X) t ( t ):= B(X) v(t) so that equation ( 5 ) may be written as Taking into account the arbitrariness of w, equation (9) yields a system of 3N non-linear equations expressed as f(d):= Sn [f TIT jQBTP(d)dV where f may be partitioned as f = - p,NTbdV - s, NTpdr= 0 in analogous fashion with d. (10) 2638 P. PAPADOPOULOS AND R. L. TAYLOR Consider an incremental formulation based on ( 5 ) and restrict attention to the Jinite time interval (r,,, t,+,]. Determination of the displacement field and all other kinematic and kinetic quantities at time t = t n + l requires the solution of a system in the form of (10). A generalized Newton method may determine the solution iteratively according to the recursive formula d(i+ 1) - (i) i = 0, 1, . . . 1) n + l - dn+1 - A n- +1 l (i)f(d:\l):= ~ ( d ; l + (11) where A,+, is a differentiable mapping from to R3N with differentiable inverse, and for notational convenience (.),,+ := ( . ) ( r n + l). The classical Newton method is recovered as a special case from (1 1) by setting A!L = Df(d!! 1), where Df denotes the (Frtchet) derivative off, defined with the aid of equations (6),(8) and (10) as , Recalling stability condition (4) and equation (8), it follows that ap Gradwh:-GradwhdT/ > 0 aF for all non-trivial v. Thus, Df is a positive-definite operator. The remainder of this article is devoted to a generalized Newton method which is specifically designed to enhance the robustness of higher-order element approximations, while maintaining certain desirable characteristics of the classical Newton method. 3.I . A modijied Levenberg-Marquardt algorithm The failure of higher-order Lagrangian element formulations to converge using the classical Newton method for moderately large step-sizes is largely due to the high sensitivity of these elements to relative nodal motions during the iterative solution process. The above observation justifies the central objective for the proposed generalized Newton method, which is to control selected nodal motions during the early iterations, before gradually giving way to the classical Newton method. Thus, the resulting method is intended to be robust, while maintaining the quadratic rate of convergence near the equilibrium state. By way of background, the well-known Levenberg-Marquardt algorithm used extensively in non-linear least-squares parameter estimation is first reviewed. According to this algorithm, a gradient-based iteration of the form d(i+ n + l1) - d,+, (0 - {diag[wDf(dP!l)]}-'f(d','!l) (12) is combined with a classical Newton iteration, thus yielding #+I) -d:\, - - (Df(d:\,) + diag[wDf(d$\l)]}-'f(d!,?+l) (13) In the above two equations, diag[ .] is an operator defined on a 3N x 3N matrix as 3N diag[.]:= 1 (eJeT)(.)(eJeJ) J= 1 where eJ, J = 1,. . . , 3 N , are the standard orthonormal basis vectors of R3N. The choice of parameter w > 0 at each iteration is governed by a separate strategy of which details can be found 2639 HIGHER-ORDER FINITE ELEMENT APPROXIMATIONS in Reference 2. Clearly, the classical Newton method is asymptotically recovered from (13) as w+o. In the present context, a modification of the Levenberg-Marquardt algorithm is introduced, so that parameter w affects only the hierarchically interpolated displacements d. To this end, equation (13) gives way to d::," - {Df(d:! = 1!d: 1) + n[wDf(d!,$l)]}-'f(dt! 1) (14) where IC is a projection mapping on R3N x R3N, defined by 3N nC.l:= 1 (eJeT)(-)(eJeT) J=3L+1 The modified Levenberg-Marquardt method of equation (14) can be regarded as a special case of the generalized Newton method defined by equation (1l), where A:!+ 1 = Df(dP! 1) + wn[Df(df! (15) I)] It can be immediately verified that A:; 1,as defined in (15), is positive-definite, given that Df(d:! 1) is assumed positive-definite. Moreover, for large values of w, the changes to the hierarchic displacements are selectively suppressed during iterations performed according to (14). To justify this point, introduce the notation A(.)(i):= ( - ) ( i + l ) - (-)(i),and use (15) to rewrite (14) as A!$IAd$1 Omitting, for brevity, subscript n - f(dt! 1) + 1, express the above system in expanded form as fi>ii(i) fif(i) f#i) fiT(i) + [ = ] y(i)][ diag[:fiT(i)] Ai(i) = - "fd")) a, respectively. Since DfCi)is where D and 6 denote derivatives in the direction of d and necessarily non-singular, elimination of Ad(') in the above system yields A(i)A& = - f ( i ) (16) where - := - Df(i)fjf(i)- 1Df(i)+ fiT(i) + -1 diag[fjT(i)] and f(i) :=1- DT(i)@(i) - 1 f The reduced matrix %(i) is positive-definite, and its diagonal dependence on parameter o necessitates that all of its eigenvalues Imin < . . < I,,, be of order o(w). It follows that the condition number K ( & ( ~ ) )of A(i) with respect to the spectral matrix norm 11.112, which is given by 2640 P. PAPADOPOULOS AND R. L. TAYLOR so that from (16) it is seen that hence An alternative least-squares interpretation of the proposed generalized Newton method is obtained by defining the direction in which a solution is sought with reference to a modified form of equation (5). In particular, define A , according to Aw[Wh, uh] := A[Wh, uh] +0 J* (wh - 6'h).diag[6f](Uh - fih)dV (18) where the added integral on the right-hand side of (18) is clearly derivable from a strictly convex functional in uh - i,.With the aid of (6) and (7) it is seen that the derivative of the above integral along the direction Ad yields the term Using nodal quadrature to evaluate (19) and noting that fi vanishes by construction at nodes I = 1 , . . . , L, it follows that the contribution of the above term to the generalized Newton method is identical to the proposed modified Levenberg-Marquardt method. The preceding interpretation is particularly attractive for implementation to finite element approximations that do not make explicit use of hierarchic interpolation. 3.2. Convergence analysis Convergence of the proposed generalized Newton method is established in the sense of Ostrowski's theorem (Reference 6, p. 300). The main idea is to write a typical iteration of the method in the standard fixed-point from (1l), and establish that, given a fixed point d* E R3N, p(DG(d*)) -= 1 (20) where p(DG(d*)) is the spectral radius of D G at d*. Ostrowski's theorem states that d* is a point of attraction for the iterative method, hence there exists an open neighbourhood Y c R3N of d*, such that iterative process (11) converges to d* from any starting point d'O'E9. In order to show that (20) holds, use (11) to conclude that, since d* is a fixed point, DG(d*) = I - A - '(a*) Df(d*) (21) where I is the 3N x 3N identity matrix. Omitting, for brevity, the explicit reference to dependence on d*, define a symmetric matrix Q as Q := Df - DGTDf DG and use equation (21) to write Q = (A- 'Df)T(2A - D f ) (A- 'Df) HIGHER-ORDER FINITE ELEMENT APPROXIMATIONS 2641 Since, by construction, 2A - Df is symmetric positive-definite and A-'Df is non-singular, it follows that Q is a symmetric positive-definite matrix. Let (A, z) be a (generally complex) eigenpair of DG, and note that, since Df is symmetric positive-definite, ZTDfz > 0 where Z is the complex conjugate of z. Invoking the positive-definiteness of Q, it follows that ZT(Df- DGTDf DG)z = ZT Df z - itTDGTDf DG z = iiTDfZ- IAI2ZTDfZ> 0 so that IAI < 1, and condition (20) is met. Thus, the proposed generalized Newton method is convergent for constant w, when interpreted as a fixed-point iteration method. The above convergence analysis is motivated by the work of Weissinger' on iterative methods. 3.3. Implementation of the algorithm The proposed solution method is applicable to all higher-order Lagrangian finite element interpolations. For the purpose of demonstrating its implementation, consider a two-dimensional elastic body ~8,which is discretized by means of 9-noded isoparametric quadrilateral elements. The proposed method will be employed to selectively suppress the relative motions of the four mid-edge nodes and the centre node, during the early iterations of every solution step. Therefore, within a typical element Re write the piecewise polynomial displacement field as 4 9 a::= where 5 and 9 are co-ordinates in the parent element ((5, 9 ) I < 1, (91< 1}, shown in Figure 1. Also, ieand fie are standard and hierarchic nodal displacement vectors, respectively. In addition, the element shape functions N,' are written as + ;(I tr<)(I + qrq) +(I - [')(I + qIq) j ( 1 + trt)(l - q 2 ) (1 - t2)(1- q 2 ) for for for for I = 1,2,3,4 I = 5,7 I = 6,8 I =9 where ((1, q r ) are the co-ordinates of node 1. Figure 1. A 9-noded quadrilateral element 2642 P.PAPADOPOULOS AND R. L. TAYLOR S Figure 2. Representative transition functions In order to provide a smooth transition from a lower-order dominated solution to the full-order solution during the iterative process, it is important to vary parameter w in (14) following a suitable strategy. In particular, w should be 'large' when d'') is away from the solution, so as to restrict the motion of mid-edge and centre nodes relative to the corner nodes. Likewise, w should decrease as d") approaches the solution, so as to allow for the classical Newton method to provide asymptotic quadratic rate of convergence. In order to facilitate the above transition, define the residual norm at iteration i of the current time step n + 1 as given by and introduce the scalar parameter s!,? RIP! 1 s A 1 := log,,- R'" n+ 1 Subsequently, let o be prescribed according to where wois the initial value of w and 0 is a smooth transition function which decreases with s, e.g. e(s) = - as') if s < O if s > 0 where a > 0 is a parameter that controls the rate at which w decreases as the global solution is approached, see Figure 2. An alternative implementation of the proposed solution algorithm can be obtained without partitioning the degrees of freedom into standard and hierarchic. To achieve this, one may merely amend the element stiffness matrix of a standard 9-noded isoparametric quadrilateral with a contribution due to the term in (19). HIGHER-ORDER FINITE ELEMENT APPROXIMATIONS 2643 4. NUMERICAL PARADIGMS The proposed formulation has been implemented within the environment of FEAP, a generalpurpose non-linear finite element program partially documented in References 5 and 8. Selected numerical simulations have been conducted using a compressible isotropic neo-Hookean constitutive model. The strain energy function Y of the model is defined as Y = f p ( z C - 3) - p In ZZZ:” + f A(ZZZ;’~ - where C is the right Cauchy-Green deformation tensor, Zc = tr C,ZZZc = det C and A, p are the Lame constants. Here, 8- and 9-noded quadrilateral isoparametric elements are used for the discretization of all deformable bodies, and the particular finite element formulation is based on the work of Simo and Tay10r.~ 4.1. A single element test This problem is specifically designed to test the convergence of the proposed generalized Newton method and its dependence on the choice of parameter o.A single 8-noded isoparametric element is used in the analysis of an elastic body in the state of plane strain. The reference Xz,X,)I ( X II < 1, IX2I < I}, and its deformaconfiguration of the body is given by R = {(XI, tion is induced by prescribed boundary conditions U(- 1, & 1, x3;t ) = 0 and 13(1 - X$)t U(1, X,,X3;t) = 0-5(1 - X : ) t [ O I The properties of the isotropic elastic material are 1 = 16644.43and p = 33.36 (which correspond to E = 1000 and v = 0.499). The current configuration at time t = 1.0, as obtained by the finite element approximation, is shown in Figure 3. Current configuration Referenceconfiguration Figure 3. A single element test problem; reference and current configuration 2644 P. PAPADOPOULOS AND R. L. TAYLOR Figure 4. A single element test problem; binary representation of convergence. in (fi, wc)-space A one-parameter family of initial configurations for the generalized Newton method is defined by a mapping x ( O ) on the reference configuration, so that x:oyx,, t ) = XI P X , ( l x:, = X:"'(X,, t ) = x, + jX,(l x;), x:"' = x:"' x(o) = 3 - - - A = 1, 2, 3 (24) (0) x3 (XA?t ) = x3 where jis a scalar parameter. It can be readily verified that invertibility of 2"' requires IflI < 1. A solution to the above boundary-value problem is subsequently attempted in a single time step of the proposed generalized Newton method, where the solution parameters which appear in equations (22) and (23) are set to ct = 0.2 and 0 6 wo d 1, and the initial configuration is defined according to equations (24) with - 0.5 < p < 0.5. The results of the parametric analysis are illustrated in Figure 4. In this figure, convergence or divergence of the iterative process is indicated by white or gray shading, respectively. The analysis clearly demonstrates that the radius of attraction of the proposed solution method expands significantly with increasing values of coo, as expected. For comparison, it is seen that the classical Newton method (wo = 0) is unable to provide a convergent solution for the vast majority of the prescribed initial configurations. 4.2. Pressing of a square block A square block in plane strain is pressed uniformly on its top surface using displacement control, while its bottom surface is kept fixed. The block, which has linear dimensions of 24 by 24 2645 HIGHER-ORDER FINITE ELEMENT APPROXIMATIONS Figure 5. Pressing of a square block; spatial discretization iia = -5.0 uz = -7.5 Figure 6. Pressing of a square block; deformed configurations units, is spatially discretized by 9-noded isoparametric elements exploiting symmetry, as shown in Figure 5. The properties of the elastic material are as in the previous paradigm. Figure 6 shows the deformed configurations obtained in the first two convergent steps, where U2 = - 5.0 and ti2 = - 75, respectively. The solution parameters that appear in equations (22) and (23) are set to w o = 10.0 and ct = 0-2. In contrast, the classical Newton method requires several solution steps to reach the configurations depicted in Figure 6. 2646 P. PAPADOPOULOS AND R. L. TAYLOR 5. CONCLUSIONS A generalized Newton method is presented in the context of finite elasticity. The method is shown to be convergent and appears to suppress the tendency of higher-order Lagrangian elements to fail due to excessive distortion during the early iterations of a classical Newton solution. Thus, the proposed method allows the use of higher-order Lagrangian elements with time steps well beyond the limits imposed on meshes by a traditional Newton-based solution procedure. REFERENCES 1. K. Levenberg, ‘A method for the solution of certain non-linear problems in least squares’, Q. Appl. Math., 2, 164-168 ( 1944). 2. D. W. Marquardt, ‘An algorithm for least-squares estimation of non-linear parameters’, J. SOC.Indust. Appl. Math., 11, 431441 (1963). 3. M . F. Beatty, ‘Topics in finite elasticity: hyperelasticity of rubber, elastomers, and biological tissues-with examples’, ASME Appl. Mech. Rev., 40, 1699-1734 (1987). 4. R. W. Ogden, Non-Linear Elastic Deformations, Ellis Horwood, Chichester, 1984. 5. 0.C. Zienkiewicz and R. L. Taylor, The Finite Element Method; Basic Formulation and Linear Problems, Vol. 1,4th edn, McGraw-Hill, London, 1989. 6. J. M. Ortega and W. C. Rheinboldt, Iterative Solution o/ Nonlinear Equations in Several Variables, Academic Press, New York, 1970. 7. J. W eissinger, ‘Verallgemeinerungendes Seidelschen Iterationsverfahrens’, Z . Angew. Math. Mech., 33, 155-163 (1953). 8. 0. C. Zienkiewicz and R. L. Taylor, The Finite Element Method; Solid and Fluid Mechanics, Dynamics and NonLinearity, Vol. 2, 4th edn., Mc-Graw-Hill, London, 1991. 9. J . C. Simo and R. L. Taylor, ‘Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms’, Comput. Methods Appl. Mech. Eng., 85, 273-310 (1991).

1/--страниц