close

Вход

Забыли?

вход по аккаунту

?

738

код для вставкиСкачать
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).
Документ
Категория
Без категории
Просмотров
2
Размер файла
695 Кб
Теги
738
1/--страниц
Пожаловаться на содержимое документа