ADV COMPUTNL FLOWS & TRANSPORT
ADV COMPUTNL FLOWS & TRANSPORT CAM 397
Popular in Course
Popular in Computational
SOCI 101S - 02
verified elite notetaker
verified elite notetaker
verified elite notetaker
verified elite notetaker
verified elite notetaker
This 165 page Class Notes was uploaded by Frida Senger on Sunday September 6, 2015. The Class Notes belongs to CAM 397 at University of Texas at Austin taught by Staff in Fall. Since its upload, it has received 17 views. For similar materials see /class/181767/cam-397-university-of-texas-at-austin in Computational at University of Texas at Austin.
Reviews for ADV COMPUTNL FLOWS & TRANSPORT
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 09/06/15
v7 5 Clue n WM e y 4 v 1 31 iawacvs W39C I firmwaaue lt w AM LAN k 5 N i w p I 5 d t e v C lt Q v A a K L if f J g 7 l quot 6 r39 17 v Chapter 1 Review Formulation of Linear Static Problems Suggested ni 1 1 Preliminaries Let nscg 2 or denote the number of space dimensions of the problem under consideration Let Q C Rn be an open set with piecewise smooth boundary F A general point in Rn is denoted by 13 We will identify the point 1 with its position vector emanating from the origin of Rnsd The unit outward normal vector to F is denoted by n We shall employ the following alternative representations for a and n 336122 2 where mi and m 1 g i g nsd are the Cartesian components of a and 12 respectively see Fig Ll Unless otherwise speci ed we shall work in terms of Cartesian components of vectors and tensors 2 CHAPTER 1 REVIEW FORl leLrTIOzI OF LINEAR STATIC PROBLEMS Figure 11 Figure 221 We assume that T admits the decomposition rerJm am Where r jmzn gm and F9 and Th are open sets in F The notations are de ned as follows U is the set union symbol Thus Pg UPh means the set of all points 1 contained in either Pg or Th Also is the set intersection symbol Thus Pg Th means the set of all points contained in both IQ and Th The empty set is denoted by 0 Thus means that there is no point x contained in both F9 and Th ie F9 and Th do not intersect or overlap A bar above a set means set Closure7 ie7 the union of the set with its boundary Thus anr up Reference linear For further discussion see book We shall assume throughout that Pg D but allow for the case Th 2 B Let the indices i j k I run over the values 1 Jam Differentiation is denoted by a comma eg 213i 2 1th flu822 and repeated indices imply summation eg in R3 uriz wii 1ng 35333 2 Wei8172 Wu9y2 ya822 The summation convention only applies to the indices i j k and 1 and only to two repeated indices If there are three 0 more repeated indices in an t2 I LASSICAL LINEAR HEAT COND UCTlON S I RONG AND WEAK FORALS39 3 expression then the summation convention is not in effect This is in keeping with the usual convention Divergence theorem Let fz R be Cl Then Lfidnanidr 15 Integration by parts Let f be as above and also let 925 R be C1 Then fgdo fgidrfgnidr 16 Q Q P 12 Classical Linear Heat Conduction Strong and Weak Forms Equivalence Let 1 denote Cartesian components of the heat ux vector let u be the temperature and let f be the heat supply per unit volume Assume the heat ux vector is de ned in terms of the temperature Use f a gradient by the generalized Fourier law f for 3 t load a qz 2 Mind Iiij is symmetry 1 f and f load vec Where the conductivities mj s are given functions of 3 If the mj s are constant throughout Q the body is said to be homogeneous The conductivity matriac n Hijli is assumed positive de nite namely pTxp 2 0 for all nSd vectors p5 and anp z 0 implies p 2 0 The most common situation Def post is in practice is the isotropic case in which Isl ax Hiw zj where dais the Kronecker delta 394 CRAP TER i R E VI W FO R M U LATiO N O F L I NEA R S TATHT P R0 B L EMS A formal statement of the strong form of the boiinlaryvalue problem is as follows Given f Q R gle 4 R and h fh v R nd uz gt R such that Gm f in 9 heat equation 18 S u 2 g on Pg 19 qmz h on Th 110 Where q is de ned by 17 Use The functions 9 and h are the prescribed boundary temperature and heat ux respectively This 2 for be and h problem possesses a unique solution for appropriate restrictions on the given data f esh size 6 In a more mathematical terminology 18 is a generalized Poisson equation 19 is a Dirichlet boundary condition and 110 is a Neumann boundary condition We shall now recall the weak formulation of the boundary value problem 19 and 110 will be treated as essential and natural boundary conditions respectively Let 8 denote the trial Del asbefore solution space and V the variation space S and V consist of realvalued functions de ned on 5 Del this time satisfying certain smoothness requirements such that all members of S satisfy 19 Whereas if in E V then 100 on Pg 111 No boxes thru The weak formulation of the problem goes as folloyvs 12 CLASSICAL LINEAR HEAT CONDUCTION S i RONC AND W EA Ii FORMS 1 Given fzfl 115121 a R and h H R nd a E S such that for all an E V W7 wygqid zf wfdQ lj whdl 112 V 39 Q Q Fh where q is de ned by 17 Theorem If all functions involved are suf ciently smooth then a solution of S is a solution of W and vice versa It is convenient to introduce an abstract version of 112 Let awu magmaude 113 wf Awfdn 114 mm frwhdl 115 Then 112 may be written as am u w f hlr 116 It is useful to also adopt an indexfree notation for 112 Aside from stemming the proliferation of indices this formulation is conductive to the computer implementation of the element arrays especially in more complicated situations such as elasticity In introducing our index free notation we shall assume for de niteness that nsd 2 Let V 6 CHAPTER 1 REVIEW FORMULATION OF LINEAR STATIC PROBLEMS denote the gradient operator thus uyl Vu 1152 2 117 2 7131 Va 2 w 11s w2 In the case of two space dimensions the conductivity matrix may be written as 511 K12 n I Lij symmetric 119 N21 1522 In the isotropic case 119 simpli es to 1 0 K 13527 2 R 0 1 In terms of the above expressions the integrand of 113 may be written in indexfree fashion wwuw V39wTI V39u 121 Thus in piece of 113 we may write awu f9VwTnVu d9 122 13 HEAT COND UCTION GALERKIN FORMULA HON quot Figure 12 Figure 241 13 Heat Conduction Galerkin Formulation Symmetry and Positivede niteness of K Let Sh and V11 be nite dimensional approximations to S and V respectively We assume all members of Vb vanish or vanish approximately on Pg and that each member of 5h admits the representation uh vh gh 123 Where vb E Vb and gh results in satisfaction or at least approximate satisfaction of the boundary condition it g on Pg The Galerkin formulation is given as follows Del cf Given f g and h as in WU nd uh vb gh E Sh such that for all wh E Vb G atwhwh W what mm 124 We now View our domain as discretized into element domains 26 1 g e g neg in two dimensions the element domains might be simply triangles and quadrilaterals see Fig 12 Nodal points may exist anywhere on the domain but most frequently appear at the element vertices and interelement boundaries and less often in the interiors Del rel 8 CHAPTER Z REVIEW FORMULATION Of LINEAR STATIC PliOBLEil IS in simple 0nedimensional models the global nodal ordering and ordering of equations in the ma trix system may coincide ln multidimensional applications this would prove to be an inconvenient restriction with regard to data preparation in what follows a more exible scheme is described Let 77 2 12 nnp the set of global node numbers where mm is the number of nodal points By the terminology g node we shall mean a node A at which it is prescribed that uh 9 Let 77g C 77 be the set of g nodes 7 The complement of 779 in 17 denoted by 77 97g is the set of nodes at which uh is to be determined The number of nodes in 77 174 equals neg the number of equations A typical member of Vb is assumed to have the form whm Z Man 125 AEW Wg where NA is the shape function associated with node number A and CA is a constant We assume throughout that wh 0 if and only if CA 0 for each A E 77 77g vhm 2 Z imam 126 AEW WQ where dA is the unknown at node A ie temperature and ghee 2 warns gt914 Aeng 94 gtmi 127 From 127 we see that 911 has been chosen to be the nodal interpolate of g by way of the shape functions Consequently gh will be generally only an approximation ofg See Fig 13 Additional sources of error are the use of approximations fh and hh in place of f and h respectively and 2 domain approximations in which the element boundaries do not exactly coincide with F Analyses 13 HEAT CONE TCTION GALERKIN FORMULATION 9 Figure 13 Figure 242 of these approximations are presented in Strang and Fix Substituting 125 127 into 124 and making use of the fact that the 058 are arbitrary see 1 for more details results in Change mi 2 1NANBdB NArf NAh Z aNANBgB A E 7 77g 128 BEnng 136779 To de ne the global stiffness matrix and force vector we need to rst specify the global ordering of equations For this purpose we introduce the ID array sometimes called the destination army which assigns to node A the corresponding global equation number viz Global equation number A p P if A E 3923 77g lDA 129 0 if A E 779 where 1 g P g neg The dimension of ID is nnp As may be seen from 129 nodes at which 9 is prescribed are assigned equation number zero An example of the setup of ID and other important dataprocessing arrays is presented in Sec 15 The matrix equivalent of 128 is given as follows Kd F r 130 10 CHAPTER Z REVIEW FORMULATION OF LINEAR STATIC PROBLEMS K K1962 ddq FFp 13Paneq 131 Km 2 mivs PIDltAgt ozmw 132 Fp 39NAfgtNAhr ZazvwNBgB 133 86779 The main properties of K are symmetry and positivede niteness 14 Heat Conduction Element Stiffness Matrix and Force Vec tor We can break up the global arrays into sums of elemental contributions nel K 2K6 K5K3Q 134 61 6 F 2 FC FeFfa 135 621 Where KgQ 1NAN B QVNATnVNBdQ 136 F NAf NAh Z aNANBegB 86779 2 fVAfdflf iVAhdFquot Z aJVA7YBSQB 137 9 3 FE Belg P Fh fe PID4 QIDB 138 See Fig 14 for an illustration of Iquot The element stiffness ice and element force vector F may be deduced from these equations 5b 22m Lgmgnm 139 14 HEAT BOND UCTION ELEAHfh S I IF FNESS MATRIX AND FORCE VEQ IquotOR ii Figure 14 Figure 251 3 auvwvbyzf VlVaT sVNbdQ 140 96 f 2 Nafdnf NahdF Zkgbgf 141 95 r bl where nen is the number of element nodes and g gm of g is prescribed at node number 1 Del recall and equals zero otherwise An implicit assumption in localizing the gterm is that if 934 is not a node attached to element 6 then NAW 2 0 for all a E T29 Otherwise the last term in 137 may involve gdata of nodes not attached to element 6 which is not accounted for in 141 The global arrays K and F may be formed from the element arrays Ice and F respectively by way of an assembly algorithm as described in Change mi The element stiffness matrix can be written in a standard form convenient for programming Ice BTDBdQ 142 Qe where in the present case D 2 K 143 nsd gtlt nsd B 311an 144 Tlsd X Hen Ba 2 Viva 1 40 v I 12 CHAPTER 1 REVIEW FORMULATION OF LINEA R STATIC PROBLEMS The component version of 142 is 35 APBEDBdeZ 146 Let 4 d6 d6 g 397 V a 14 nan X l diam Where e h e da u 148 d8 is called the element temperature vector It can be shown that the heat ux at point x E 2 3 can be calculated from the formula qm DwBmde Dmf Bad 149 11 15 Heat Conduction Data Processing Arrays ID IEN and LM The element nodal data is stored in the array IEN the element nodes army which relates local node numbers to global node numbers viz lEN 7 L50 Local Element Global node number node number number 16 CLASSICAL LINEAR ELASTOSTATICS STRONG AND WEA K FORMS 13 The relationship between global node numbers and global equation numbers as well as nodal boundary condition information is stored in the ID array see 129 In practice the EN and D arrays are set up from input data The LM array is de ned from the relation LMae IDlENae 151 Because of the previous relationship we often think of LM as the element localization of ID Strictly speaking the LM array is redundant However it is generally convenient in computing to set up LM once and for all rather than make use of 151 repeatedly We postpone further description of these arrays until we consider elasticity which is more general than the present case In terms of the IEN and LM arrays a precise de nition of the 93 s may be given see 141 e 0 if LMae 0 9A if LMae 0 Where A lENae This de nition may be easily programmed 16 Classical Linear Elastostatics Strong and Weak Forms Classical elastostatics is a rich subject in its own right See for background and references to the literature These references range from the very physical to the very mathematical The most physical book is the one by Timoshenko and Goodier In ascending order of mathematical content are Sokolnikoff Gurtin DuvantLions F ichera The reader is reminded that indices 2 j k and 5 take on values 13nsd where nsd is the number of spatial dimensions and the summation convention applies to repeated indices 239 j k and 1 only Del 6 Del 6 141 CHAPTER 1 REVIEW FORMULATION OF LINEAR STATIC PROBLEMS Let or denote Cartesian components of the Cauchy strese tensor let ui denote the dis placement vector and let fz be the prescribed body force per unit volume The in niteez mal strain tensor is de ned to be the symmetric part of the displacement gradients Viz def um up nag strain displacement equations 153 The stress tensor is de ned in terms of the strain tensor by the generalized Hooke s law Uij CijlekJ 154 Where the czjkz s the elastic coe icients are given functions of m If the cijkz s are constants throughout the body is called homogeneous The elastic coef cients are assumed to satisfy the following properties Symmetry Cijkl CW major symmetry 155 6239ij sz kl mmor symmetmes 156 Cijkl Cijlk Posit ive de niteness CemWWW m CD H m 1 Cijkl1bij1rbkl 0 implies i10 399 t OI 00 V for all m E 3239 and all sz 2 w Note The positivede niteness condition is in terms of symmetric arrays 16 L LASSICAL LINEAR ELASTOSTATICS STRONG AND lVEAK FORMS 15 A consequence of the major symmetry 155 is that K is symmetric The first minor symmetry Change per implies the symmetry of the stress tensor ie m a The positive de niteness condition when combined with appropriate boundary conditions on the displacement leads to the positive de niteness of K In the present theory the unknown is a vector pie the displacement vector A formal statement of the strong form of the boundarywalue problem goes as follows Given fill gt R gizl g R and hizl h gt R find uz gt R such that 015139 fr 2 0 in 2 equilibium equations 159 53 ui 2 g on F9 160 o zjnj 2 on Th 161 where aij is de ned in terms of ui by 153 and 154 Remark The functions g and hi are called the prescribed boundary displacements and tractions respec tively is sometimes referred to as the mixed boundaryvalue problem of linear elastostatz39cs Under appropriate hypotheses on the data possesses a unique solution see Fichera The additional complexity of the present theory when compared with heat conduction is that the un known ie u 2 is vectorsvalued rather than scalar valued In practice it is important to be able to deal with somewhat more complicated boundarycondition speci cation See 1 for further CllSCllSSlOll Reference ti book Let Si denote the trial solution space and Vi the variation space Each member ng E 8 satis es ence linear book l6 CHAPTER 1 REVIEWquot FORMULATION OF LINEAR STATIC PROBLEMS ui gi on Pg whereas each wi E Vi satis es mi 0 on Pg Equation 161 will be a natural boundary condition The weak formulation goes as follows Given fi 9 gt R g Pg R and hi l h w R nd ui 6 Si such that for all 39wi 6 Vi w wtjaijdnszifidnqu wihidl 162 9 Q Fr where 025 is de ned in terms of uz by 153 and 154 In the solid mechanics literature is sometimes referred to as the principle of virtual work or principle of virtual displacements 20 being the virtual displacements It can be shown that a solution of S is a solution of and viceversa See The abstract notation for the present case is awu Aw dpijklwhnd w fgwifida 164 amp 2 wihidl 165 Th Let S 2 u I m 6 Si and let V w f 20339 6 Vi Then the weak form can be concisely written in terms of163165 as follows 16 CLASSICAL LINEAR ELASTOSTATICS STRONG AND WEAK FORMS Given f g and h in which the components are de ned in WW nd u E S such that W39 for all w E V awuwfwhr 156 As discussed in Sec 12 it is desirable to construct indexfree counterparts of the expression on the righthand sides of 163165 For concreteness we shall assume that nsd 2 thus 1 g ijkl S 2 According to conventional notation e uh03 is the strain tensor However we will have no need for this matrix7 and consequently we reserve 6 for the vector of strain components U11 Iu U2 2 1012 U2 711 60 61w 158 1112 1012 1021 This will considerably simplify subsequent writing Note that factors of one half have been elimi nated from the shearing components Le last components of 167 and 168 Compare 167 and 168 with 153 A similar notational rede nition occurs with respect to the vector of stress components 169 Modify Modify 18 CHAPTER 1 REVIEW FORMULATION OF LINEAR STATIC PROBLEMS Table Ll Table 271 I z j J k l l l l 3 l 2 3 392 l 2 2 2 Let D11 D12 D13 D DU 2 D22 023 170 symmetric D33 Where Dij Ciij 171 Change at add in which the indices are related by Table 11 diag s As should be clear we have collapsed pairs of indices 239 j and 19 into single indices I and J respectively taking account of the symmetries of Cz39jkz WM and wad Observe that the indices I and J take on values 1 2 and 3 In nsd dimensions the I and J indices will take on values 12 nsdmsd 12 It can be shown that w j6ijkzucz ElwlTD Wl 172 and so awuQewTDeudQ quot 173 17 ELASTOSTATICS GALERKIN FORMULATION l 17 Elastostatics Galerkin Formulation Symmetry and Positive de niteness of K Let Sh and V11 be nite dimensional approximations to 5 and V respectively We assume members 11211 6 Vb result in satisfaction or approximate satisfaction of the boundary condition quot111g 2 0 on F9 and members of SI2 admit the decomposition uh vh gb 174 Where vb E Vb and 912 results in satisfaction or approximate satisfaction of the boundary condition ui 2 g on F9 The Galerkin formulation of our problem is given as follows Given f g and h as in nd uh vb g11 E S11 such that for all wh E Vb h arm vh th what arwhah 175 To de ne the global stiffness matrix and force vector for elasticity it is necessary to introduce the ID array This entails a generalization of the de nition given in Sec 15 since in the present case there will be more than 1 degree of freedom per node For elasticity there are nsd degrees of freedom per node but in order to include in our de nition cases such as heat conduction we shall take the fully general situation in which it is assumed that there are ndof degrees of freedom In 20 CHAPTER 1 REVIEW FORMULATION OF LINEAR STATIC PROBLEMS this case Global equation number A P HAen ID i A 2 176 V V 39 Degree of Global 0 if A E 77g freedom DOde number number where 1 lt 239 lt ndof Thus ID has dimensions ndof gtlt Tldof if ndof 2 1 we reduce to the case considered previously in Sec 15 ie ID239 A 2 IDA Recall that 77 2 1 2 nnp denotes the set of global node numbers Let 779 C 77 be the set of nodes at which gi and let 77 ng be the complement of 779 For each node in 77 779 the nodal value of is to be determined The explicit representations of v and 9 in terms of the shape functions and nodal values are H s vb 2 NA dz 177 AET 779Degree Of freedom node number number 9 5 Piagt4 178 A6779 Let ez denote the ith Euclidean basis vector for Rn ei has a 1 in slot i and zeros elsewhere For example 715d 2 392 e1 2 62 2 179 The vector versions of 177 and 178 may be de ned with the aid of 23 viz v v ei L80 17 a ciA ELASTOSTATICS GALERKIN FORMULATION 21 9h 2 931 181 Likewise a typical member 11211 E V has the representation wh wglei w Z NAciA 182 AEnmg Substituting 1777 178 and 180 182 into 175 and making use of the fact that the s are arbitrary results in verify Deli refi ndoj Z Z 0NAeiaNBejij NAeivf NAeuhh i1 BEH 779 ndof Z Z 01VAej1VBEjij A E 77 779 1 S i g ndof 183 j1 86779 This is equivalent to the matrix equation 239 to ndof Kd F 184 where K IS pQ 185 d z 4Q 186 F Fp 187 IX PQ CLUVAehNBej 188 ndof FP NAenfNAeuhP Z Z 0NAeiaNBejng 189 j1 BE VIQ 22 CHAPTER 1 REVIEW FORMULATION OF LINEAR STATIC PROBLEMS in which P 21131111 Q mugs 190 Equation 188 may be written in more explicit form by using 173 and noting that see 16 and 168 emAei BAei 191 Where for example NA1 0 M 2 BA o NA2 192 NA 2 NA1 With these 188 becomes Degree of freedom numbers KPQ ZeganTADBB Cler 193 V Global Global node numbers equation numbers and the indices are related by 190 Equation 189 is also amenable to explication Note that by 164 NM 1 LNAfZ d9 194 and likewise by 165 JVTAehhh 2 N Ahidf 195 Fr 18 ELASTOSTATICS ELEMENT STIFFNESS MATRIX AND FORCE VECTOR 215 Thus 189 may be written as ndof Fp 21VAfidQl VAllidll Z Z aJVA y1vgej39gj5 9 W 11 Bang Now that we have de ned K we can establish its fundamental properties K is symmetric and if appropriate displacement boundary conditions are prescribed then X is also positive de nite Positivede niteness of K is based upon two requirements a positive de niteness condition on the constitutive coef cients and suitable boundary conditions being incorporated into Vb See 1 for details Reference book 18 Elastostatics Element Stiffness Matrix and Force Vector As usual K and F may be decomposed into sums of elemental contributions These results will be omitted here as the reader should now be familiar with the ideas involved cf Sec 14 We will proceed directly to the de nitions of Ice and F as k8 lkiqla f3 2 f5 1 S paq nee name 197 kgq eKQeBZDBbdoej p new 1 i 198 I ned5 1j i a 0 Ba 2 0 Nag nsd 2 2 199 1 an Nva l 24 CHAPTER 1 REVIEW FORMULATION OF LINEAR STATIC PROBLEMS and f z Naf do Nahidl giggly r r rtquot 1100 Where g gjb gja if 93 is prescribed at node 6 and equals zero otherwise It is useful for programming purposes to de ne the nodal submatrix 25 BZDBtdQ 1101 V 9e nedxned From 198 we see that kg e3quotkbej 1102 This means the pqcomponent of Ice is the ijcomponent of the submatrix kgb By 197 through 199 we see that ksquot3 may be written as k5 BTDB d9 1103 Qe Where BB1BgBnm 1104 For example in the case of a two dimensional 1e nsd z ned 2 fournoded element Ice takes the form kil 162 166133 194 e 1631 332 933 54 y d k 2 110 8 X 8 k 1 1932 1633 63 ki 1632 1333 1934 Add dash In practice the submatrices above the dashed lines are computed and those below if required are determined by symmetry 19 ELASTOSTATICS DATA PROCESSING ARRAYS ID IEN AND LM 25 The global arrays K and F may be formed from the element arrays k6 and 6 respectively by way of an assembly algorithm as outlined in Change ref Let di d 3 d 3 d 1 l v a quot l 39 neexl dim dis ned 2 a 2 1107 dig Where d umg 1108 d8 is called the element displacement vector The vector of stress components 169 at point a 6 QB can be calculated from the formula Modify and nan am DmBmde DmZBamd 1109 41 19 Elastostatics Data Processing Arrays ID IEN and LM we have already noted that the de nition of the ID array must be generalized for the present case as indicated in Sec 17 We must also generalize our de nition of the LM array However the lEN array remains the same as before In the present and fully general cases the LM array is threezdimensional with dimensions 26 CHAPTER 1 REVIEW FORMULA I I39OEV OF LINEAR S l t39l l PROBLEMS neg gtlt rim X 7561 and is de ned by LMiaeJlDilENae 11110 Degree of Local Element freedom node number number number Alternatively it is sometimes convenient to think of LM as two dimensional with dimensions n68 gtlt n81 viz LM p 7 6 l LMiae p nadi a 1 gt2 1111 Local Element equation number number To see how everything works in practice it is helpful to run through a simple example Example 1 Consider the mesh of four node rectangular elements illustrated in Fig 15 We assume that the local node numbering begins in the lower left hand corner for each element and proceeds in counterclockwise fashion This is shown for element 4 which is typical Four displacement ie gtype boundary conditions are speci ed namely the horizontal displacement is speci ed at nodes 1 and 10 and the vertical displacement is speci ed as nodes 1 and 3 Since nnp 1 1 2 ndof 2 ned 2 2t and 4 displacement degrees of freedom are speci ed we have neg 20 We adopt the convention that the global equation numbers run in ascending order with respect to the ascending order of global node numbers However in practice equation numbers are often renumbered internally to minimize the bandwidth of K and thus decrease storage and solution effort The ID EN and LM arrays are given in Fig 16 The reader is urged to verify the results 19 ELASTOSTATICS DATA PROCESSING ARRAYS ID EN AND LM 10 I Figure 15 Mesh 0f four node rectangular elasticity elements global and local node numbers element numbers and displacement boundary conditions 28 ID array IEN array LM array CHAPTER 1 REVIEW FORMULATION OF LINEAR STATIC PROBLEMS Figure 16 ID IEN and LM arrays for the mesh of Fig 15 19 ELASTOSTATICS DATA PROCESSING ARRAYS ID IEN AND LM 29 In terms of the IEN and LM arrays a precise de nition of the gg s may be given see 1100 0 if LM ae 75 0 g 9 1112 gik Where A IENae ifLMiae 0 This de nition may be easily programmed 30 CHAPTER 1 REVIEW FORMULATION OF LINEAR STATIC PROBLEMS Bibliography 1 TJR Hughes The Finite Element Method Linear Static and Dynamic Finite Element Anal ysis PrenticeHall Englewood Cliffs New Jersey 1987 2 G Strang and GJ Fix An Analysis of the Finite Element Method PrenticeHall Englewood Cliffs New Jersey 1973 3 G Duvaut and JL Lions Les In quations en M canique et en Physique Dunod Paris 1972 4 G Fichera Existence theorems in elasticity in C Truesdell ed Handbuch der Physik Vol VlaQ lVIechanics of Solids II Springer Verlag New York 1972 5 M Gurtin The linear theory of elasticity in C Truesdell ed Handbuch der Physik Vol Via2 Mechanics of Solids II Springer Verlag New York 1972 6 LS Sokolnikoff Nlathematical Theory of Elasticity 2nd ed McGrawHiLl New York 1956 7 S Timoshenko and JN Goodier Theory ofElasticity 3rd ed McGraw Hill New York 1969 31 Finite Element Methods for Nonlinear Problems by Thomas JR Hughes Robert M Ferencz Division of Applied Mechanics Durand Building Stanford University Stanford California Copyright 1990 TJR Hughes and RM Ferencz Hughes and Ferencz Nonlinear FENI Chapter 2 Draft Chapter 2 Solution Algorithms for Nonlinear Problems All linear problems are trivial and all nonlinear problems are impossible R Sachs Berkeley 1970 Outline 21 Basic Features of Linear and Nonlinear Finite Element Formulations 22 Newton type Methods 221 Consistent Tangents 222 The Newton Raphson Method 223 Developing an Understanding of Some Aspects of Nonlinear Equation Solving 23 Practical Enhancements for Nevvtontype Methods 231 Line Search 232 A QuasiNewton Method BFGS 24 Arclength Methods 241 Initialization 242 Convergence Criterion Hughes and Ferencz Nonlinear FEE Chapter 392 Draft 21 Basic Features of Linear and Nonlinear Finite Element Formulations Linear Statics Recall from Chapter 1 the conceptual algorithm used to develop a finite element method 5 ltgt W a G gt M V V V V Strong form Weak Galerkin Matrix of BVP The first step begins with a statement of differential equations and boundary conditions the strong form of the boundary value problem S To develop a nite element method we need to convert S into a weak form ie an integral formulation lV incorporating socalled natural boundary conditions Discretization of the weak form via the Galerkin method leads to a nite dimensional problem Introduction of a speci c basis for the nite element functions results in a linear algebraic7 or matrix problem llI K d F 211 For linear elastostatics the terms in the matrix equation are denoted the stiffness K the nodal displacements d and applied external loads F The nature of typical FE functions leads to a bandpro le coef cient matrix which for linear elasticity and heat conduction is symmetric and positive de nite Fig 211 Bandpro le structure of the matrix arising from typical nite element functions Symmetric go to Hughes and Ferencz Nonlinear FENI Chapter 2 Draft Linear Dynamics In linear dynamics we proceed along similar lines time LSZ ltgt iiV 2 G ltgt llI a stepping Strong form tcontinuous semi ODE S algoathm 0f IBVP discretization linear The strong form is now an initial boundary value problem the differential equations and boundary conditions being augmented with initial conditions The classical framework for FE formulations uses a time continuous weak form derived from 5 through integration of the spatial domain only The discretization is performed similar to the static case the FE bases functions are spatially de ned but the nodal variables remain time continuous This is denoted a semidiscretization and adoption of speci c FE function spaces leads to a matrix system of ordinary differential equations Time discretization is then accomplished through a timestepping algorithm Which results in a sequence of linear algebraic problems e g NPAa R 212 Where M is the effective mass matrix Aa is a vector of unknown acceleration increments and R is the residual or outofbalance force The algebraic problem is accompanied by update formulae which emanate from the particular timestepping algorithm Both linear static and dynamic problems discretized by Galerkin s method are reduced to linear algebraic equations which may be generically denoted as the system An b 213 Again for many common applications A is symmetric and positive de nite The Lin ear Text presents in detail a reasonably ef cient data structure and algorithm for direct solution of such problems Nonlinear Statics With the above as background we examine the case of nonlinear statics The con ceptual background is as before but we defer consideration of the details until Chapter 3 v N Nonlinear S 4 W G gt problem Hughes and Ferencz Nonlinear FER Chapter 2 Draft We start with a strong form of the nonlinear boudary value problem and develop an equiv alent weak formulation Introduction of typical nite element functions to approximate the weak form with a Galerkin method leads to a nonlinear problem Roughly speaking for typical cases we seek unkown internal forces to balance speci ed external forces ie Fint z Fext As an example consider elastostatics in which the internal forces may be represented as nonlinear algebraic functions of the unknown displacement vector F m Nd 215 Where N is a vectorvalued function N an IR and neg is the number of unknown degrees of freedom In components N Np 1 lt P g neg 216 where P is the equation number In contrast with the linear cases 211 and 212 the present case requires solution of a nonlinear algebraic problem Notice however that there is no matrix apparent in 214 the lefthand side only contains a vectorvalued function Any matrix brought into the solution of 214 is the result of an algorithmic strategy The typical stategy for solving nonlinear algebraic problems is to reduce the solution to a sequence of linear algebraic problems by way of incremental iterative strategies These methodologies will be discussed in the next section Nonlinear Dynamics It is also the case in nonlinear dynamics that within each time step we will have to solve a nonlinear algebraic problem The details of nite element formulations in nonlinear dynamics have not been treated yet and will not be discussed until Chapter Here we will only present a sketch which follows along the lines of the conceptual algorithms presented previously just to illustrate this point N 1 time 3 39 lV z gt Grim z stepping ODES algorithm Hughes and Ferencz Nonlinear FERI Chapter 2 Draft Again we begin with a strong form of the problem which is a nonlinear system of dif ferential equations in space and time appended with boundary and initial conditions We develop an equivalent weak formulation involving spatial integrals This is approximated with a semi discrete Galerkin method which results in typically nonlinear ordinary di 39er ential equations For example in elastodynamics the Galerkin method would yield the nonlinear system of secondorder ODEs Md39 Fint Fm 217 where Fint z N d is a nonlinear function of the displacements The introduction of a time stepping algorithm will result in a nonlinear algebraic problem having similar struc ture to the static case eg NAa R 218 but where here the unknown might be for example the acceleration increment Aa 22 Newtontype Methods Let us assume the external load is parameterized by time or equivalently a load param eter which we denote by t For each value of t of interest we wish to achieve a balance between the external and internal forces For elastostatics we would seek the solution dt such that Ndt E Fina FEX U 991 The difference between the external and internal forces if dt is only approximated is called the residual force RU Rt 13 Pena Fi quott 222 As a rst example of a simple algorithm we introduce the so called incremental load method We mention that every nonlinear nite element solver is incremental with respect to loading Application of the entire external load in one step is not a feasible solution strategy as it would likely lead to excessive computational costs or simply a lack of converegence for signi cantly nonlinear problems to C Hughes and Ferencz Nonlinear FENI Chapter 392 Draft The incremental load method is summmarized in Box 221 Note that in the right hand side we have applied the external force at time tn1 and subtracted the internal force at time 73 based upon the approximate solution dn to form the residual Rn 1 A matrix K is introduced but left momentarily unde ned and is a vehicle for solving for an increment Ad This is a linear algebraic problem of the type discussed previously It enables us to compute a solution increment which is used to update the solution Once updated we go on to a new step with a new approximate solution but in going to the next step we increment the external load and so forth Box 221 Algebraic de nition of the incremental load method n 0 Step loop KAd Rn 223 def Fez1 Fint dn1 d Ad 994 n e n 1 continue step loop The weakness of such a solution strategy is shown in Fig 221 The residual is never driven to zero and this error entails not obtaining the exact solution curve Rather some approximate solution is obtained which for typical nonlinear problems may poorly represent the exact solution In order to drive the residual to zero in each step we need to introduce some type of iteration Thus the incremental load method although simplicity personi ed motivates the need to iterate within load steps to drive the residual to zero and thereby obtain accurate solutions to the nonlinear algebraic problem However before pursuing that point we must examine the definition of the matrix K This will be the focus of the next subsection Hughes and Ferencz Nonlinear FEZVI i Chapter 2 Draft Feast approx I Nd exact Kt Fig 221 Use of the incremental load method leads to divergence from the true solution path 221 Consistent Tangents We now discuss the consistent tangent matrix arising from two nonlinear model problems Presentation of the associated consistent element tangent matrices and other element quan tities will occur in Chapter 3 Nonlinear Elastostatics Let us consider as a rst example nonlinear elastostatics In this case the internal force is a nonlinear algebraic function Fint N The consistent tangent stiffness is de ned here to be dig 8Nd 9 9 quot39 ltquotquot 39a or in components N Np 226 3N GNP r 9 7 K ad ital quot 1 M where P and Q are equation number indices 1 S P Q S neg As just shown the consistent tangent is de ned via a differentiation operator This has implications with respect to the symmetry of this matrix Socalled hyperelasticity assumes the existence of a scalarvalued 9 4 Hughes and Ferencz Nonlinear FERI Chapter 2 Draft strain energy function U 1 such that 8U N 22 ad 8 ie the internal forces are the gradient of a potential In components 9U V 1 P adp 2 2 9 This implies the symmetry of the consistent tangent as follows 8N1 82U K K 22 9 adQ 8dQ8dp Qquot l 0 due to symmetry of mixed partial derivatives Consequently K K T for hyperelasticity which is a valuable computational property The consistent tangent is so far unmotivated It can be motivated by linearizing the nonlinear equations We proceed as follows for elasticity Fe 2 FintidnH Ndn1 2211 Expanding in a Taylor series 3Ndn Ndn ad Ad 2212 In order to develop an equation which can be used to de ne Ad we ignore higherorder terms in the Taylor series and rearrange in the following fashion 8Ndn 6d So the linearized operator ie matrix engendered by the Taylor expansion leads to the Ad Fext N01 2213 consistent tangent 6N d K M 2214 Now the fact that we have used the Taylor expansion suggests the nature of the approximation This is illustrated in Fig 222 for the singledegreeof freedom case Note that the residual which is de ned to be Rn1 KdnAd is a reasonable approximation to the difference of the internal forces Ndn Ad NM This illustrates in a simple setting the important property of a consistent tangent it gives a good local approximation Also from the gure it is clear that the exact solution dn1 is approximated by dn plus the Ad 28 Hughes and Ferencz Nonlinear FEM Chapter 2 Draft Fact 1 Kd z slope Nd d dt Fig 222 The consistent tangent provides the best local approximation to the loaddisplacement curve computed via this linear matrix equation T his fact holds true in general for according to Taylor s theorem the consistent tangent is the best local approximation in the sense IlKdAd N01 f Ad NdH 0 as gt O In fact K d is the unique matrix achieving this property So far we have seen that consistent tangent matrices can be advantageous computa tionally in that they may be symmetric and that they also have superior local approxima tion properties However it is important to realize that consistent tangents are not always symmetric In fact for many typical eld theories the tangent will not be symmetric An illustrative case is that of nonlinear heat conduction Nonlinear Heat Conduction As we will discuss in Chapter 3 in nonlinear heat conduction the internal force vector can be written in the particular form F 2 Kat 2215 29 Hughes and Ferencz Nonlinear FEBI Chapter 2 Draft where K is a symmetric matrix and d is a vector of nodal temperatures This is a man ifestation of the Galerkin method and the Fourier law of heat conduction see Chapter 3 for details The consistent tangent in the spirit of a Taylor expansion is the gradient of the internal force with respect to its argument namely aFim 817 22 84 l at l lt 16 Let us compute this quantity as follows 8KPR CIR A aKPR d 22 quot39 ado KPQ adQ R 11 where the repeated subscript R implies summation Note that two terms are present the first clearly being symmetric Symmetry of the second term would rely on the condition allpg l afi QR 397 r 8d asp 3918 However it will be veri ed in Chapter 3 for the simplest possible case that 2218 does not hold whenever K exhibits any temperature dependence ie whenever the problem is nonlinear rather than linear In this eXample the band pro le structure for a typical finite element system would look as shown in Fig 223 Although the pro le is symmetric the corresponding entries within the pro le are not equal ie it has a symmetric architecture but it is a non symmetric matrix This results in a doubling of storage and solution costs These costs can be extremely high for large systems and one naturally asks if this extra storage and computation is worth the potentially superior performance of the consistent tangent No unequivocal answer can be given In simple cases when the nonsymmetric contributions to the tangent are small it is probably adequate to iterate with only the symmetric part This algorithm is adopted in many codes However if the nonsymmetric contributions are sub stantial one would anticipate that their inclusion would lead to signi cant computational efficiencies So in many cases the extra storage and solution eitort is justi ed Remark As previously shown in elasticity the basic material component of the stiffness is symmetric Here we comment on the inclusion of so called follower forces in which the external force is dependent upon the deformed geometry ie Fe Fex d In this 210 Hughes and Ferencz Nonlinear FEi lI Chapter 2 Draft Fig 223 Band pro le structure of the consistent tangent matrix for nonlinear heat conduction Components nonsymmetnc 39 case the gradient of Fe needs to be included in the consistent tangent Recasting the linearization the statement of equilibrium FeXtdn1 Ndnl is expanded in the Taylor series aFext dn 6d leading to the linear algebraic system amen arena 6Ndn ext F dn ad AdNdn ext Ad 2 221 Typically the contribution from the follower force is nonsymmetric Consequently even in simple cases such as elasticity a nonsymmetric equation solver can be useful for nonlinear analysis I 0 w FM Hughes and Ferencz Nonlinear FEM Chapter 392 Draft 222 The NewtonRaphson Method Now we come to an algorithm whose superior convergence properties are based upon the use of the consistent tangent the Newton Raphson method A basic Newton Raphson scheme is presented in Box 2221 where we suppress the step number subscript n and n 1 to simplify the writing The steps are as follows Here we iterate within each step to drive the residual to a suf ciently small value thereby ensuring an accurate approximate solution to the nonlinear algebraic problem in each step We introduce an iteration counter 239 initialized to O The initial iterate d is set to the solution from the previous step dn The linearized problem involves the consistent tangent K 115 and we solve for an increment in the solution Ad The righthand side is the residual calculated from the current iterated value In the rst iteration we use Rm 2 Rdn which is subsequently re evaluated as iteration proceeds and the displacement updated dm ll d Ad Now we check if dii l l is an accurate solution This convergence test can be quite elaborate in practice and involve several different quantities However for present purposes we only present perhaps the most basic criterion which compares the current residual with the initial value for the step We require that some norm eg the 2onorm of the residual be reduced to a preassigned fraction of the initial residual This fraction is de ned by the userinput parameter a which might have a value of 10 3 10 etc If the convergence test is not satis ed we typically continue the iteration loop solve the linearized problem with respect to the updated residual update the solution and so forth If convergence is attained the solution at step n 1 is de ned to be the latest iterate which satis ed the convergence test At this point we update the the step number and begin the next step Remark We Wish to emphasize that in the ertonRaphson method the matrix K is indeed the consistent tangent In practice the terminology NewtonRaphson method is often used simply to denote algorithms in which a new lefthand side matrix is formed for each iteration irrespective of that matrix s analytic de nition However strictly speaking the term NewtonRaphson implies use of the consistent tangent I The drawback of the NewtonRaphson method is the expense of formulation and factorization of the consistent tangent matrix at each iteration As mentioned this is a largescale computation and can entail great expense This motivates the practical expediency of replacing the consistent tangent with something simpler For example in the modi ed Newton method K may be taken to be the rst tangent formed for the 212 Hughes and Ferencz Nonlinear FEM39 Chapter 2 Draft Box 2221 Algebraic de nition of the NewtonRaphson algorithm n 0 Step loop i O di dn Iteration loop Kdlt gt Ad 2 R 2222 d f new d 1 d Adm 2223 Convergence test RH1 5 120 2224 No i lt i 1 continue iteration loop Yes and dlt51gt 2225 n n 1 continue step loop entire analysis or the rst tangent formed in each step Subsequently K would be frozen for the remaining iterations in the analysis or in the step In general a modi ed Newton method implies that some convenient approximation to the consistent tangent is employed in place of forming K in each iteration In order to appreciate the nature of convergence of the NewtonRaphson method let us examine the following concept the order of convergence of an iterative scheme To illustrate this concept consider a single degree offreedom example Let d be the exact solution and d be the approximation at iterate i The error is denoted a 2 d d 2226 We attempt to obtain an expression governing the evolution of successive iterates of 213 Hughes and Ferencz Nonlinear FERI Chapter 2 Draft the form I lead lt c 18ml 2227 Where c and k are constants Now if the error is small at iteration i obviously the greater the power k the smaller the error will be after the succeeding iterate For the Newton Raphson method k 2 2 This is referred to as quadratic convergence and represents an exceptional property of this algorithm The interpretation of quadratic convergence is as follows Assume the hypothetical data lam 9 and 621 Then according to the inequality 2227 subsequent interations behave as given in Table 2221 Note that as we get close to the solution the property of quadratic convergence implies that the number of leading Zeroes in the error doubles with each succeeding iteration Thus for the NewtonRaphson method the convergence is quite rapid once we obtain an approximation close to the solution Table 2221 An example of secondorder convergence It erate Error bound i am lt c 1501ki 936 Q13 000016 4030 quothector H CO In the modi ed Newton method the order of convergence is only k 1 The implica tions for a hypothetical case are illustrated in Figs 224 and 225 Note the effect of the consistent tangent in the NewtonRaphson method It more faithfully represents the slope of the curve and quickly leads to convergence In contrast the modi ed Newton method employing the initial tangent in this example requires many more iterations to converge It is recommended to the reader that 5 he carefully examine the figures and verify that the illustrations do graphically represent the solution strategies Now clearly in this example the NewtomRaphson method attains more rapid convergence However there is a tradeoff due to the savings in computational effort engendered by the modi ed Newton method as 214 Hughes and Ferencz Nonlinear FEM Chapter 2 Draft the number of formulations and factorizations of K is signi cantly reduced Again a rigid discrimination between the two methods is dif cult without considerable computational experience It Feat Nd1 0 do do d2d Fig 224 An example interation path with the Newton Raphson method Facet III N dm 7 NEdmg N 39n in O d0 d1 C12 d Fig 225 An example interation path with the modi ed Newton method to L4 On Hughes and Ferencz Nonlinear FENI Chapter 2 Draft 223 Developing an Understanding of Some Aspects of Nonlinear Equation Solving In practice with large systems it is often impractical to update the tangent at every iteration ie use the NewtonRaphson method We can assume then that often we will only use a matrix which is an approximation to the consistent tangent Consequently in some sense there is a possibility this matrix will be soft or stiff relative to the consistent tangent What are the implications of using too soft or too stiff an approximate tangent matrix To gain an appreciation for this behavior it is possible to consider a very simple model problem that is in fact linear This illustrates the important point that often it is possible to understand nonlinear behavior by a properly constructed linear model problem Consider the scalar equation Kd F 2228 where K is a positive stiffness or spring coef cient A Newton Raphsontype algorithm for this problem is Characterized by the following equation fr Ad F Kai 2229 The righthand side is the residual at iterate i Rm and the lefthand side coef cient represents some approximation of the true stiffness K Consequently we can write K aK where a is a nondimensional parameter we assume a is positive If a lt 1 then the approximate stiffness Ii is too soft If a gt 1 then Ii is too stiff Let us analyze the behavior of this method We can write the algorithmic equation as follows aKd 1gt dm F ltd 2230 We wish to subtract from the algorithmic equation the identity 0 0 written as aKd d 2 F Kd 2231 Subtracting these equations leads to a recursion relationship for the error aKeltquot1gt em Ks 2232 where we recall em 2 d d Note that the exact tangent cancels out of this recursion which may be rewritten as elti1gt 44 2233 216 Hughes and Ferencz Nonlinear FEM Chapter 2 Draft Where the ampli cation factor A is de ned is ma quot 1 z 1 i 2234 0 O A it is clear that the error will be reduced in each iteration if lt l implying the algorithm will converge Conversely gt 1 means the error increases each iteration ie the algorithm diverges We plot A as a function of a in Fig 226 The ideal value of A is 0 meaning that eff O and the exact solution has been obtained This value of A corresponds to a 1 which denotes use of the exact stiffness Furthermore any value of A close to 0 ie values of or close to 1 represent good approximations to the exact stiffness and lead to rapid convergence Now if a is increased this represents a stiff approximation to the true stiffness and A asymptotically approaches 1 Thus for extremely stiff approximations convergence While attained Will be quite slow because A R 1 implies little reduction in the error for each iteration On the other hand if or is reduced below 1 the algorithmic behavior is somewhat different At the value a the ampli cation factor A 1 A too soft an approximate stiffness If lt K results in divergence of Newton Raphsontype Further reduction in a ie o E 0 leads to A lt 1 and gt 1 Consequently using algorithms Again we emphasize that using too stiff an approximate stiffness still results in convergence but may require many more iterations than the optimal case using the exact stiffness Remark The interpretation given above motivates a comment on plasticity Fig 227 schematically illustrates typical behavior of an elasto plastic material such as a metal The material begins by deforming elastically then flows plastically subsequent unloading will also be initially elastic Now for example in a plasticity calculation a consistent tangent may be formed while the material is plastically owing This tangent is quite soft relative to the elastic responSe If during the iterative process the material elastically unloads the consistent tangent would become stiff again However if the soft plastic tangent is retained for subsequent iterations we can anticipate the behavior illustrated by the previous singledegreaof freedom example in other words divergence This type of behavior is observed in practice for precisely this case Similar results have also been recorded for some implementations of contact algorithms where using excessively soft tangent stiffnesses leads to catastrophic divergence I 217 Hughes and Ferencz Nonlinear FERI Chapter 2 Draft A n 1 A gtlasa oo good f Stiff 1 Soft Fig 226 Newton Raphson method s error ampli ca tion factor as a function of a E Fig 227 An example stressstrain response with elasticplastic loading and elastic unloading Hughes and Ferencz Nonlinear FERI Chapter 2 Draft Is this singledegree offreedom model problem representative of the behavior of large scale systems The following example illustrates that indeed it is Consider a linear multi degreeof freedom problem in which the analogous iterative solution strategy is applied KAdU F th39 2235 dlti1 do Adm 2236 This time we introduce the symmetric positive definite matrix K as an algorithmic ap proximation to the true stiffness matrix K which is also assumed to be syrmnetic positive de nite To analyze this system we construct the associated eigenproblem K ARM 2 0 2237 The eigenvalues and eigenvectors Awi q satisfying 2237 may be used to uncouple the multidegreeof freedom system into a set of scalar equations ME E A dga for each I lneq 2238 It is important to emphasize that the use of the eigenproblem is not to be construed to have any algorithmic signi cance relevance It is simply a device which permits analysis of the matrix problem with regards to the previous results for the scalar equation Dropping the modal subscript Z we can derive a recursion for the modal error di1 5di F Ad 2239 d d F Ad 2240 H 1 1 0 We 13 Ag 2241 where A is interpreted as a modal ampli cation factor Note that d is de ned to be the exact solution of the uncoupled mode ie F 5 Ad We can denote an equivalence between 2241 and the previous single degreeof freedom case 2233 by de ning a modal a 5 1 This implies that we have reduced each mode of the coupled system to the previous analysis Consequently a good approximate left hand side matrix K is one such that each is nearly unity implying a z 1 for each mode 1 being the ideal value Furthermore the clear implication from the scalar case is 219 Hughes and Ferencz Nonlinear FEM Chapter 2 Draft that soft modal components in the sense of the eigenproblem 2237 will cause the system iteration process to diverge Finally the linear system is interpreted to characterize the behavior of the nonlinear system locally These analyses provide some insight into the important behavior of stiff and soft com ponents of the system One conclusion we draw is that if we are working with a soft tangent it is likely that it is less than half as soft as the consistent tangent in some modes Consequently an unsophisticated use of a modi ed Newton method could lead to catastrophic divergence as the Fourier coef cients of certain softly represented components grow This is not satisfactory ie not suf ciently robust for practical nonlinear equation solving strategies for largescale engineering analysis This motivates the inclusion of addi tional algorithmic strategies intended to preclude catastrophic divergence and enhance the probability of rapid convergence Without engendering unacceptable computational costs Ix lx39 O Hughes and Ferencz Nonlinear FERI Chapter 2 Draft 23 Practical Enhancements for Newtontype Methods 231 Line Search As we have seen previously the use of a tangent matrix which is not a reasonable ap proximation to the correct one can lead to slow convergence or actual divergence of Newton type algorithms In large scale problem solving the cost of the NewtonRaphson method is often prohibitive Consequently when using a modi ed Newton method there is a high probability that we are sometimes iterating with a tangent which is in some sense locally too soft and will lead to divergence We desire a means of preventing such divergence and the line search algorithm has proved effective in this role The line search algorithm may be viewed as an insurance policy for Newtonlike methods If the nonlinear iteration is behaving reasonably then line search adds relatively little computational cost and has little effect upon convergence of the Newton process For instance if we are near the solution and the Newton algorithm is converging quadratically then line search would not interfer However should the nonlinear iteration attempt to diverge then line search would intervene to control the solution and prevent divergence Linear Case To illustrate the line search algorithm we rst examine the linear case We consider again a Newtonlike iterative strategy for solving the linear algebraic problem K d F The iterative driver is written as before KAdl F th39 R 231 where K is an approximation to K but the update formula d 1 d 5 Ad 232 now employs Ad as a search direction That is rather than simply updating di1 d Ad the search direction is scaled by the line search parameter 5 prior to updating the previous solution The objective is to nd a value of s such that the solution is in some sense improved relative to the nonscaled update We will present two derivations which under the current circumstances lead to identical definitions of s The rst derivation follows from the observation that the potential energy 1 1F ddi ef 3Frat dTF 233 221 Hughes and Ferencz Nonlinear FERI Chapter 2 Draft is minimized by the exact solution d K IF Introducing the function p m as F W t in w 234 we Wish to choose 3 such that p is minimized The stationary values of p are derived from the condition dp ds leading to the explicit formula for the search parameter 0 i Adm nd W o S Adli KAdi 35 or equivalently Ad Rh 0 9 S Adm KAdi 36 Note that both the numerator and denominator are quadratic forms involving matrices assumed to be symmetric and positive de nite Consequently if Ad 0 then s gt 0 Possible addition Furthermore under these assumptions d2 2 i 19 2 Ad KAd gt 0 237 so the stationary value of 3 de ned by 236 is indeed a local minima The second derivation follows from a different point of view The new residual Ra may be written in terms of s 12 F Kd 1 F K 01 stud 238 Select 3 such that RU has zero component in the direction Ad viz Ad 120 0 239 This design condition leads to the same formula for 3 as given in 236 Although the results for the present case are identical this alternative derivation is attractive due to its generality It may be utilized to de ne s in cases where no potential energy function exists lx O 3 Hughes and Ferencz Nonlinear FEBI Chapter 2 Draft Nonlinear Case We now write the Newtonlike iterative strategy of 231 and 232 in the form appropriate for the nonlinear algebraic problem Fext Fintdi di1 dh 32 Adi 2311 Again we use the residual on the righthand side of 2310 and assume that K is some approximate tangent matrix Employing the philosophy of the second derivation above de ne the scalar valued nonlinear function C39 2 Ad Rslt gt 2312 Ad Fext Fi d sltigtAdltigt 2313 The desired design condition for 3 is Gsi 0 Finding a scalar root of a scalarvalued nonlinear function may not appear to be a computationally intensive operation However note that each evaluation of 63 for a trial value of 3 involves computing both the external force and the internal force contributions from every element in the mesh The latter may require the evaluation of complex nonlinear constitutive relationships and these operations can certainly be nontrivial Thus if an extended iterative process is required to compute s such that C30 z 0 the line search algorithm would be unacceptably expensive Consequently less restrictive criteria are typically employed for the line search parameter Matthies and Strang suggest accepting 3 if GHQ salsa Where G0 is based on the residual derived from the external forces at the current step minus the internal forces at the previous step For details of the procedure advocated by Matthies and Strang we refer the reader to their paper This procedure is encoded in the software included in Chapter in subroutine SEARCH and similar algorithms are included in many nonlinear nite element codes available to the general public Remarks 1 The algorithm just described is a general procedure for nonlinear problems No assumptions have been made concerning the structure of the internal force vector nor the o 93 Hughes and Ferencz Nonlinear FEAI Chapter 2 Draft lefthand side coef cient matrix If there exists a potential function Ud such that rmt z Nd 2 2314 then the nonlinear version of the rst derivation of line search leads to the same nonlinear scalar equation 2313 given for s This is the case for nonlinear elastostatics Where Pd Ud dTFe 2315 Note that this assumes the external loads are not a function of d a situation commonly denoted as dead loading If Fe F e 1 it may still be possible to de ne a more complicated potential which would admit the application of the rst derivation 2 For situations in which the consistent tangent is not symmetric and positive de nite a minimization of the residual is frequently used to define 3 through the condition min Rd sAd Rd sAd 2316 Procedure of this kind are frequently found in fluid dynamics The Generalized Minimum Residual GMRES algorithm embodies this concept and has shown great success in such applications ref Shakib Saad amp Schultz etc 3 Possible addition We reiterate the perspective that line search is intended to insure against catastrophic divergence of Newtondike methods If we envision iterating with an extremely soft approximate tangent this will lead to the accumulation of displace ment increments which are much too large leading to divergence In this case the line search intervenes by scaling back these large displacement increments and hopefully pre eluding divergence In contrast line search may in some instances enhance convergence but is not employed with this expectation For example if K is a stiff approximation to the consistent tangent the line search may return values of s gt 1 However at least one developer ref Hallquist has observed that such overrelaxation can lead to poor con vergence in the presence of highly nonlinear phenomena such as contact Thus the NIKE codes limit the acceptable value of the line search parameter to s S 10 h Ix A Hughes and Ferencz Nonlinear FEDI Chapter 2 Draft 232 Quasi Newton Methods BFGS As discussed previously the line search algorithm cannot in general enhance the rapidity of convergence of Newtonlike procedures The scaledback displacement increments returned by the line search may prevent divergence yet this does not imply that they are optimal increments The only way to capture these increments accurately is to have an effective representation of the tangent This can of course be accomplished through the Newton Raphson method but this entails the unattractive computational costs discussed earlier Certain amounts of tangent information can be inferred from the solution and residual iterates computed within any Newtonvlike method Such socalled update procedures are now used in a number of nonlinear nite element codes and represent an economical means of instilling the coef cient matrix with some approximate tangent information it Fext Ndlti1gt a ARU Ndltigt dbl dil di2 d Fig 231 A singledegree offreedom example of the quasiNewton update Quasi Newton updates are at times also labeled secant updates which is perhaps more descriptive of their motivating principle We rst examine the singledegreeof freedom case Which is illustrated by Fig 231 It is assumed that iterates i and i l have been completed ie d and 19 and the corresponding internal forces are known The ultimate goal of the Newtonlike iteration is to compute the displacement at corre spending to the applied external load Fe Having computed 03 the immediate goal 225 Hughes and Ferencz Nonlinear FEBI Chapter 2 Draft is to establish a linearized problem which gives a reasonable approximation to the tangent and whose solution permits us to update 550 to dm39g Observe that knowing the inter nal forces at iterates i and i 1 is equivalent to knowing the change in the residual over the ith iterate This may be seen as follows AB 120 R 2317 Fart 1Vdi1 Fart w Ndigt Nd 1 Ndlt gt 2319 Question Should we mention that 2318 and 2319 are only strictly true in the presence of dead loading Assuming the cost of computing the true tangent K di1 to be excessive we seek a good approximation From the diagram we see that the secant K connecting states z39 and i 1 is a reasonable choice Let us de ne the secant as follows def AR5 and then use this secant in the next iteration to de ne di2 KM dlti1gt R0 5 Fe Ndltquot1gt 2321 This is a viable strategy and it develops a useful update method For a single equation the secant condition 2320 uniquely de nes Iii Our real interest is in systems of equations and the generalization of this concept of a secant to the multi degreeoffreedom case We write an equation analogous to 2320 which is applicable to systems as follows K dquot1gt N ARlt 3 2322 Again we assume that d and dm39l are known The righthand side can be computed as shown for the SDOF case thus 2322 is a constraint upon the unknown K It is called the quasiNewton equation Obviously for multidegree offreedom systems it does not uniquely de ne K it is simply a restriction upon the definition of K The discipline of quasi Newton updates provides de nitions for R such that the quasi Newton equation is satis ed Many such procedures have been proposed eg Broyden 226 Hughes and Ferencz Nonlinear FEMquot Chapter 2 Draft Davidon BFGS etc The reader may consult lE Dennis and JJ More JE Dennis and RB Schnabel Luenberger Murray Gill and Wright for surveys of these updates The process of quasiNewton updating embodies some approximate tangent informa tion which can have a bene cial effect upon convergence Nevertheless the possibility of catastrophic divergence still exists consequently quasiNewton updates are often combined with line searches Introducing the search parameter into the update formula at d SmAdm 2323 the quasi Newton equation becomes 1230mm AR 2324 In the presence of line search this equation is the appropriate design constraint for the de nition of K The BFGS Inverse Update Perhaps the most widely used quasiNewton update in nite element analysis is the BFGS update This update is named for its creators Broyden Fletcher Goldfarb and Shanno and has the form K l d f I va I I va 2325 where v and w are update vectors Whose de nitions will be made precise subsequently We caution the reader that the use of the inverse matrices in the update equation is merely formal As is typical in numerical linear algebra the actual computation of the inverse matrix is not implied Rather we View 12 1 as a representation of the ability to solve a linear system of equations with K as the coef cient matrix If a direct equation solver is to be employed as we will assume for the remainder of this section then I might signify a factored form of K stored in a pro lecolumn data structure Due to the product structure of the BFGS update if K 1 is symmetric then K 1 will inherit that symmetry Another property of this update is that if I 1 is positive de nite and v w 1 then K 1 is also positive de nite If we multiply through the righthand side of 2325 it may be found that K l is the sum of I and two rank one matrices For this reason the BFGS update is denoted a ranktwo update These design conditions ie the product form 2325 the requirement 2 w 1 and satisfaction of 2227 Hughes and Ferencz Nonlinear FEM Chapter 2 Draft the quasiNewton equation lead to unique de nitions of v and w Other quasi Newton updates will of course satisfy the same quasi Newton equation but are derived from other auxiliary design conditions It is important to understand the data structure and sequence of computations in volved when performing the BFGS or any other quasi Newton update within a nite element code Obstensibly the update formula implies that the pro le structured I 1 is transformed in general into a full matrix However we emphasize that 3 1 is never explicitly calculated We only need a means of replicating its action upon any vector Consider the need to solve the system KAdW I RU Viewing K 1 as an operator we write A620 124126 2326 and substitute the update formula 3 I va K 1I 10275120 2327 The essential computational feature of this update is that the operations proceed from right to left Ada I vaflt1Rlt U v Rlti1gtw 2328 or V I w39vTIquot1R 2329 After solving the system RAJ R with the factored form of if AdUH I vaTmJ 2330 An w Ail1 2331 The factored 1239 is never altered during these calculations It is simply recalled as required to perform the forward reduction and backsubstitution with R Note that these operations with the factored K remain the dominant cost of solving the system KAdltiU z 12 Each of the left and right updates eg R Ri1vRi1w only entails one dot product and one SAXPY see remark below requiring 4n operations In contrast the forward reduction and backsubstitution use 43726 operations Where B is the mean bandwidth of K Furthermore this righttoleft strategy preserves the sparsity of Again a full neg gtlt neg matrix is not stored rather the two neqolength vectors 1 and 228 Hughes and Ferencz Nonlinear FERI Chapter 2 Draft 10 are simply appended to the existing storage Neither the storage scheme for K nor the value of its coef cients are in any way altered Remark The term SAXPY is derived from the LINPACK 3 subroutine of that same name which performs the operations 2 ax y ie a scalar multiple of a vector summed to another vector AX Plus Y The S denotes the single precision version of the subroutine but also makes it easier to pronounce This sequence of operations is highly ef cient on vector processors so numerical linear algebra algorithms are often expressed and implemented with the SAXPY as one of the basic building blocks I As the nonlinear solution strategy iterates a sequence of quasi Newton updates can be derived from iterate to iterate In this case I represents the most recently formed and factored coef cient matrix If a total of i updates have been derived since the last matrix formation then the BFGS update formula talces the form 21 I name I ve 1wi1T x I v1w1T 3 1 I w1voT X I w2v2T I w v T 2332 The solution of a linear system incorporating multiple BF GS updates again occurs through a series of operations ordered from right to left Consider such a system written in operator form as Ad i K 1R 1 39 2333 which we expand schematically to clarify the structure 239 239 1 l 1 2 i lt gtlt gtlt gt gt gt gtRlti1gt 2334 leftside updates rightside updates W direction of operations The rightside updates proceed from the outside in Ru Ri1 10 Ri1wi 2335 Rm RU vi 1R1lwi 2336 9 09 di Hughes and Ferencz Nonlinear FEftI Chapter 2 Draft R 2 RV 121 Rltquot1gtwlt1gt 2337 The factored form of K is then used to solve the intermediate system a3 2 13112 2338 Finally the leftside updates proceed from the inside out a3 mi 1111 adT gtvlt1gt 2339 mm Mm 102 Mame 2340 A30 Age 1 wh Agi 1gtvigt 2341 whereby the desired solution is attained ie7 Adm 2 Ad 2342 This multipleupdate solution strategy is summarized in Box 2321 With K unal tered the only storage required is that engendered by the set of BFGS update vectors vce wuwi k1 which occupies i X 2 X neg units of floatpoint storage Typically a strategy is adapted which limits the maximum number of updates i to some practical value eg imax 15 Attaining this limit is then viewed as a signal to the nonlinear solution strategy that K should be reformed and factorized in turn the update procedure is reinitialized and the pool of BFGS vectors overwritten as new pairs are generated Our discussion of the BFGS update has focused upon its form and operational struc ture We have not assigned explicit de nitions to the update vectors 2 and w The derivation of these definitions may be found in the previously referenced works which lead to the formulae vk d Ad 3346 mac dgf ARM altkgtRltH 2347 kgtARk Ade k def a VW 3348 230 where Hughes and Ferencz Nonlinear FEM Chapter 2 Draft Box 2321 Solution procedure for linear systems with multiple BFGS updates Given a matrix ff 3 set of BFGS update vector pairs 170 19mm1 and a residual vector Ra Initialize Rm 5 Rlti1gt Rightside update loop do k 1 z39 RUC Rk 1vi k1 Rk1wik1 2343 Solve intermediate system with factored matrix A3 1211 2344 Leftside update loop do k 1 i 133k AC kl 1004 AJk 1pk 2345 Terminate MU Adm Note the appearance of the search parameter so in the de nition of do assumes that a line search precedure is incorporated within the nonlinear solution strategy Alternative de nitions of the BFGS vectors are often used computationally These may be derived from the line search function Cs de ned in the previous section We reiterate that de nition Gltkgtslt gt Adm R dW 5kAdlk 2349 and subtract 0000 Ad R 1 2350 to find the identity Gquotgtsltkgt away Ad an r 2351 231 Hughes and Ferencz Nonlinear FERI Chapter 2 Draft This relationship may be substituted into the previous de nitions to nd I vac z 43 9 3 5 ak3kgt 3000 39 212k AR 040012quot 2353 where 31Gk3k Gquot0 am z 6mm 23a4 Further information on these formulae may be found in the paper of Strang and Matthies These alternative de nitions are the basis for the implementation found in NLEARN go OJ N Hughes and Ferencz Nonlinear FEM Chapter 2 Draft 24 ArcLength Methods So far in this chapter external loading has been assumed to be a given function of time Now we wish to trace loaddisplacement processes for which there is no a priori knowledge of the load time function The dif culties in such situations can be described with the aid of Fig 241 Suppose we insist upon solving this problem as an applied force problem in which the external load is monotonically increased This strategy should be adequate between the initial state A and the local maximum B However as the load increases beyond qB the only solutions exist well to the right ie the solution path would have to jump from B to D Solutions might then be obtained along the segment D F but the local maximum at F could lead to another jump If F is in fact a global maximum then no additional solution states could be identi ed under load control Clearly just applying the load monotonically without any knowledge of the nature of the nonlinear system may lead only to a partial solution or none at all One apparent means of obtaining a unique solution beyond point B is the application of a monotonically increasing displacement This would permit the solution path to be traced from B to C C to D and so forth Note that the solution is not unique between displacement levels dH and 18 Furthermore under displacement control G is a local maximum and subsequent increase in the displacement is dependent upon a discontinuous jump from G to I Liam painquot undo E disohcemm control quotx G I Dynamic I unto undev l displacement cannot a Dynamic 39xnm under load mmval Limit points under load control C D39 ectson a Fig 241 An example of snap buckling after Cris eld Clearly displacement control is no more reliable than load control for the solution of nonlinear systems exhibiting the local extrema typi ed by the previous example More effective solution methods are needed In structures these enhanced algorithms are called 233 Hughes and Ferencz Nonlinear FENI Chapter 2 Draft arclength methods The basic ideas are due to Riks 5 Wempner l Ramm l l Bergan ll and Cris eld 23 Other useful papers are also cited in the references In the mathematical literature similar procedures are referred to as continuation methods The essential idea is to increment on some quantity which unlike pure load or displace ment is increasing monotonically throughout the evolution of the problem In arc length strategies the quantity employed to increment the solution is the arc length of the solution curve in load displacement space Clearly this measure is monotonic as the solution path evolves Furthermore it will automatically accommodate reversals in the load or displace ment as is needed to trace the solution curve The arc length methods append a scalar arolength constraint equation to the system of equilibrium equations This introduces the variable A the loading parameter as an additional unknown into the system of nonlinear equations The load parameter automatically controls the amplitude of the applied load Fe The continuous version of the problem may be stated as follows Given a desired rate of change it of the arclength function f nd dt and t such that W at W 241 fldy A d where La 242 The arclength function f should satisfy the de nition of a norm on the neg 1 dimensional space consisting of pairs d 39 The discrete counterpart of 241 is seen more frequently in the literature It may be stated as follows Given a desired increment 5a of the arc length function f nd dn1 and An such that Fint dn 2 n ext 1 A 1F f6d ti 6a where E fidn39il quot39 dn nl quot A112 a 234 Hughes and Ferencz Nonlinear PEAI Chapter 2 Draft Let us introduce a candidate arclength function A form for f which is frequently given in the literature is ma 6A 2 0 MT diagK g6d 6 6X2 245 where c 2quot 246 qT dlasKo q q d IZO IFEX 247 b 6 01 248 and K0 denotes the first formedfactorized leftvhand side matrix Let us attempt to understand this function The parameter b scales between the relative weigth of displacement and load control as represented by the rst and second terms under the radical in 245 The coef cient c incorporates this scaling and non dimensionalizes the displacement term using the auxiliary solution vector 1 Note that the displacement vector 6d need not be dimensionally homogeneous A common example of this occurs in structural mechanics where the unknown degrees of freedom typically include both translations and rotations Due to this possible inhomogeneity the diagonal entries of K0 are employed as a metric that renders the inner product dimensionally consistent Returning to the example of structural mechanics the inner product 601T diag R0 6d has dimensions of energy To gain an appreciation for this function let us consider some special cases The arc length constraint is more easily visualized if we de ne the scalar incremental displacement variable 6dT diagKg 6d 6d d qT diagKo q 249 This permits us to place a local 6d 6A coordinate system whose origin corresponds to the last converged state attained on the solution path in A space ie dn An we rst assume I O in 245 ie socalled displacement control The arc length constraint plotted in 6d 6 space is shown in Fig 242 The function f de nes the boundaries of a vertical cylinder with radius 6a centered about the 6 axis Enforcing 235 Hughes and Ferencz Nonlinear F EllI Chapter 2 Draft equality of the arc length constraint requires the solution lie on the boundary of the shaded region thus the two points Where the solution curve intersects that boundary are the only admissible converged solutions for state dn1n1 The constraint in the discrete problem 243 can be loosened to f6d 6A 3 6a In this case any point on the solution curve within the shaded region is admissible as the next converged state The assumption 1 1 in 245 leads to socalled load control In this case the constraint region is a horizontal cylinder of radius 6a centered about the 6d axis as illustrated in Fig 243 6A Fig 242 Arc length constraint region for pure displacement con trol 6a 6A 6a Fig 243 Arcslength constraint region for pure load control Hughes and Ferencz Nonlinear FEM Chapter 2 Draft Now consider intermediate cases which blend displacement and load control Selecting b a where 0 lt 5 ltlt 1 leads to an elliptic constraint region as shown in Fig 244 The ellipse has a major axis which coincides with the 6A axis and a minor radius of 2 6a This value of b emphasizes displacement control over load control but controls both quantities simultaneously In contrast setting 6 1 e emphasizes load control Then the elliptic constraint region has a major axis coinciding with the 6d axis While retaining a minor radius of 2 6a see Fig 245 Finally equally weighting displacement and load control 6A Fig 244 Arc length constraint region when displacement control dominates Fig 245 Arc length constraint region when load control domi nates 237 Hughes and Ferencz Nonlinear FEM Chapter 2 Draft Fig 246 Arcolength constraint region for equally weighted dis placement and load control by selecting b leads to the circular constraint region pictured in Fig 246 with radius 2 6a In practice computational experience favors pure displacement control ie b 0 Fig 247 illustrates why the arc length constraint can actually enhance convergence when iterating with an approximate tangent This is important when we consider going through local maxima or minima points of the solution curve where the tangent stiffness is singular In such cases one can iterate with a nonsingular stiffness and still get decent convergence The modi edNewton process shown in the gure uses the initial tangent Kn The arclength constraint ensures the intermediate iterates lie on the surface of an ellipsoid centered about the point dn An Therefore the converged solution dn An must be located at an itersection 0 1 that ellipsoid with the solution curve Without this constraint the modi edNewton iteration would converge more slowly Indeed it might not converge at all if the target load level exceeded some local maxima of the solution CUI VE We now consider a solution procedure for the constrained nonlinear system 243 As with the unconstrained case discussed in Sec 221 this is based upon a consistent linearization However now we linearize with respect to both unknowns ie dn1 and An1 We begin by substituting an e d l Ad 2410 An1 i A331 AA 2411 into 243 and linearize with respect to Adm and AM It is important to bear in mind 238 Hughes and Ferencz Nonlinear FERI Chapter 2 Draft n1 6A 7quot 6d d Fig 247 Modified Newton iteration path controlled by arclength constraint the di 39erences implied by the A and 6 quantities The A variables denote the difference between individual iterates eg 0131 d213 Ad The 6 variables represent the total change over the step to the current iterate eg dm l dn 6di Using Taylor s n1 formula on both parts of 243 Fintd 31 Kd 11Adlt gt 3 Fi d H Adm 2412 A311 AA F 2413 and z faz j3r1 Adm dn A5511 AA Aquot 6a 2415 where 331 fwd 5W 2416 if Y 3 5dltigt 5W 2417 aod n1 Bad x 239 Hughes and Ferencz Nonlinear FEZVI Chapter 2 Draft 5f M W quoti m t 2 85971 MAW tax g 418 M 1511 dn 2419 5A0quot ASL in 2420 These linearized relations may be arranged in the nu 1 equation matrix problem def i c Rivii W K 1 Fext Adh AgilFext Fintdr11 T 2421 3f 6f 2 i am 86A m 6 39 fn1 Remarks 1 K can be replaced by any convenient lefthand side matrix say K accounting as an example for quasiNewton updates etc 2 Line searches may also be accounted for by using the line search update formula c1533 dfflrl Mud 24222 Where 3 is an approximate root of the line search function ie 0 so 43 Ad ASLIFW Fiwd l 043510 2423 3 It is not expensive to calculate the lower partition of the matrix in 2421 so normally the latest updated values are used even when R is not being updated Using the form of the arclength function given in 245 explicit formulae for these entries may be stated a i c 86 fa diagK05di 2424 n1 n1 a i t i 5 2425 86 n1 240 Hughes and Ferencz Nonlinear FEM Chapter 2 Draft Now consider the solution of 2421 Note that 58 5quot may be zero This would obviously be the case in displacement control when the coef cient 13 is set equal to zero Given that the augmented matrix is only positive semide nite7 some care must be exercized in solving the system ie the entire ncq 1 x neg 1 matrix cannot simply be factorized Conceptually we solve the upper partition for Adm as a function of the unknown AA This result is substituted into the lower partition and solved for AA which in turns permits evaluation of Ad This leads to the explicit formulae i a i O a ffhll 52 1Ad AW 3 W i W 2 2426 353 quotH as Adgl Aqu M 2427 where the intermediate quantities A3 i S K lR l 2428 q i Klpex 2429 represent the computationally intensive portions of the solution Thus the solution of 245 entails two forward reductionback substitutions with K This assumes that K can be factorized ie is nonsingular This apparently contradicts a prime motivation for adopting arclength algorithms However it appears that in practice the tangent stiffness is rarely re formed precisely at the singular point and that this strategy has been sufficiently robust to solve many problems 241 Initialization Solution of the linear system 2421 is seen to be straightforward An essential in gredient in the success of an arclength method is the initialization of the algorithm in each load step The strategy employed can determine the success or failure of the partic ular method We rst consider a simple initialization strategy which suggests some basic elements common to most algorithms Subsequently we will describe a more sophisticated strategy that has been found to perform well in practice 241 Hughes and Ferencz Nonlinear FEM Chapter 2 Draft Starting the arclength method in each step requires the speci cation of 6A which in turns determines 6510 Remember unlike the customary assumption of loadcontrol algorithms employed in FE analysis we have no a priori knowledge of an acceptable mag nitude of external load to be applied in the current step Rather some reasonable estimate is required To this end we define q d f K 1dnFe 2430 and set 6010 6A q 2431 We think of K 1dn as a tangent matrix formed at the beginning of this current step but in fact it may be any convenient approximation based upon quasiNewton updates of an earlier tangent etc With 6110 determined up to a scalar factor all that is required is to de ne 6M0 The arclength constraint determines the magnitude of 6M viz 6a fwd 5Alt gt 2432 f5lt q 6AM 2433 16Alt gtfti1deq 1 2434 implying aim f 71 2435 Remark It is important to emphasize that an adaptive strategy is required to select the arclength constraint in each step In regions of large curvature in the load de ection curve a small increment in arclength 6a is needed if convergence is to be obtained In contrast in regions of small curvature ie nearly linear response much larger increments are practical Consequently an important ingredient in an arclength method is a strategy to select 6a One such strategy is illustrated by the driver in NLEARN While somewhat ad hoc in nature the algorithm presented expands the arc length increment when convergence is rapidly attained in the previous step and contracts 6a when slow convergence is detected I In addition to determining the magnitude of 6M0 E Skill the sign of SAM is also an algorithmic variable The simplest strategy proceeds as follows If the following two conditions are true we reverse the sign of 6M relative to the initial value for the previous step ie sign 01310 sign6n 242 Hughes and Ferencz Nonlinear FEW Chapter 2 Draft The relevant conditions are i The determinant of the left hand side matrix changes ie det Kdn det dud ii In changing Sign the determinant has passed through zero as shown for example in the following sketch detK If either of the above is not true we retain the sign from the previous step ie Sign 6M3 sign 6A3 We mention that the calculation of detK39 is trivial when using a direct equation solver based upon the Crout factorization K UTDU In this case the determinant of K is equal to the product of the diagonal entries of D Remark To understand the rationale behind the strategy consider the scalar case The condition that the determinant passes through zero corresponds to a local maximum or a minimum in the loaddisplacement diagram Consider the left hand sketch On the left ascending branch the slope is positive and passes through Zero at the maximum Subsequently we would want to begin with a negative load increment to initiate descent along the right branch of the curve Analgous logic dictates the change in sign as we pass through the minimum in the righthand sketch The 243 Hughes and Ferencz Nonlinear FEZW Chapter 2 Draft important thing to note is that the determinant must have passed through zero for these to be valid pictures Yet the determinant can also change sign by passing through in nity as illustrated in this sketch deth In the scalar case this corresponds to turning points in the load de ection curve Hughes and Ferencz Nonlinear FEDI Chapter 2 Draft The desired behavior at turning points is quite different from that at maximaminima In this case we want to maintain the current direction of loading unloading as we pass through the turning point rather than reversing the load increment direction as discussed for the local extrema The behavior of detK can be monitored to differentiate between these cases by storing successive values and testing whether they pass monotonically through zero I Arclength algorithms have not yet progressed to the point where they are totally reliable Nevertheless their use expands the domain of solvable problems well beyond that attainable with Newtonlike strategies and some impressive practical calculations have been performed In practice it has been found that the following initialization strategy leads to a more robust arc length algorithm for engineering simulations In this strategy the initial solution AS1011 dial is extrapolated from recently converged states using the arclength parameter as an indepedent variable The extrapolation is based upon the quadratic Lagrange polynomials 2 a a artIX quot an ll an 2 anlan 2 quot39 an 13m a an 2a an 2437 an1 39 anan1 aft 2 1301 a and GM 2438 an quot an2an an1l These lead to the initializations 0 n1vZ anl n 2 Igian1plnl Zi9n1lln 2439 245 Hughes and Ferencz Nonlinear FENI Chapter 2 Draft digi1 liianildn2 liiarHildm i liianildn 2440 where liianH m 2441 13m 5 1quot 1 160quot1 2442 130111 6an 1 6a 6an16an 6421 2443 611 6an6an Note that 6an1 and 6an are converged arclength increments stored from previous steps While 6an1 represents the desired arclength increment selected by the algorithm at the beginning of this new step Only upon convergence of this step will the actual value of 6an1 be knowri Remark This extrapolation eliminates the need for the previously discussed logic to assign signSt01 The assumed quadratic form is capable of automatically setting signx 01 signn as the solution passes through a limit point but will also maintain sign t ll signn near turning points 242 Convergence Criteria Various convergence criteria may be employed in conjuction with the arclength method For example we can testagainst the residual norm For this purpose we de ne the residual norm to be HR is RTdiagI lt0 1R 2444 Again the inverses of diagonal entries of the initial left hand side matrix are employed as a metric to render the inner product dimensionally homogeneous cf the discusion on the arc length constraint function The residual convergence criterion may be stated in terms of a user specified tolerance 63 Where we insist that the current step satisfy the condition 112311 3 ea W 2445 ie the nal residual in the current step must be a preassigned fraction of the initial residual 246 Hughes and Ferencz Nonlinear FERI Chapter 2 Draft In addition we must ensure the arclength constraint is maintained We put this in the form of an inequality constraint corresponding to the interior of the shaded regions illustrated in Figs 241 247 viz frail 39 6an1 lt 5f 5an1 2446 The parameter 5f is another userde ned tolerance The geometric interpretation of 2446 is that the solution diquot1 A3111 must lie within an ellipse as pictured in the following sketch 6A not 6d OK Trial solutions outside the ellipse are not accepted even though they may in fact be valid solutions lying on the loaddisplacement curve Some programs also employ a displacement constraint Where the displacement norm employed is as de ned previously for the arc length function 245 Hde d deiagKo d 2447 with the same motivation for dimensional homegeneity We insist in this case that consec utive displacement increments Ad 2 111 satisfy the constraint Ad lt c act d d 9 4 48 H H d H H Ed H n1 nll where 54 is another user de ned tolerance It is worth noting that as long as b l in 245 ie not pure load control the displacement criterion is redundant with the arcclength constraint In situations Where 247 Hughes and Ferencz Nonlinear FEW Chapter 2 Draft other algorithms are employed eg where using Newtonlike methods the displacement criterion should generally be added to the residual convergence criterion As an application of the arclength method we present Fig 248 from Stanley This gure illustrates the buckling of an axially loaded cylindrical shell with various degrees of geometric imperfection It is clear from the solutions that the adaptive control of the arc length constraint has decreased 50 around the sharper portions of the loaddisplacement curve and increased 6a elsewhere 39 Buckling Ez39gemolution percc J C 1h IMP 0J 10 oh lMP C 07 n U l L 3 n g 95 o a5 50 4 1MP g 04 3 r E b E 03 quot z i a I 39 I 39 Lugond a 39 Z Z V quot3 quotquot 4513 9 x 9 39 900 901 533 05 992 MAL DEREC TION Fig 248 Postbuckling of aidally compressed cylinder with arclength method imper fection sensitivity with xed mesh from Stanley Hughes and Ferencz Nonlinear FERi Chapter 392 Draft References Section 22 1 TJR Hughes The Finite Element Method Linear Static and Dynamic Finite Element Analysis PrenticeHall Englewood Cliffs 1987 Section 22 ll J M Ortega and WC Reinholdt Iterative Solution of Nonlinear Equations in Several Variables Academic Press New York 1970 Section 23 I JE Dennis and J J More QuasiNewton Methods Motivation and Theory SIAM Review Vol 19 No 1 1977 2 JE Dennis and RB Schnabel Numerical Methods for Unconstrained Optimization and Nonlinear Equations PrenticeHall Englewoocl Cliffs 1983 3 JJ Dongarra JR Bunch CB Moler and GW Stewart LINPACK Users Guide SIAM Philadelphia 1979 4 H Matthies and G Strang The Solution of Nonlinear Equations International Journal for Numerical Methods in Engineering Vol 14 16131626 1979 Section 24 1 PG Bergan Automated IncrementalIterative Solution Schemes pp 291305 in Numerical AIethods for Nonlinear Problems Vol 1 eds C Taylor and E Hinton Pineridge Press UK 1980 2 MA Cris eld A Fast IncrementalIterative Solution Procedure that Handles Snap through Computers and Structures Vol 13 5562 1981 3 MA Cris eld An Arclength Method Including Line Searches and Accelerations In ternational Journal for Numerical Methods in Engineering Vol 19 1269 1289 1983 4 CA Felippa Solution of Nonlinear Static Equations to appear in the NorthHolland Handbook Series on Computational Methods in Mechanics Volume on Large De ection and Stability of Structures 5 E Riks An Incremental Approach to the Solution of Snapping and Buckling Prob lems International Journal for Numerical Methods in Engineering Vol 15 524 551 1979 61 KH Schweizerhof and P Wriggers Consistent Linearization for Path Following Meth ods in Nonlinear RE Analysis Computer M39ethods in Applied M39echanics and Eng39 neering Vol 59 261 279 1986 7 G M Stanley Continuum based Shell Elements PhD Thesis Stanford University Stanford California 1985 g l quot3 Hughes and Ferencz Nonlinear FEDI Chapter 3 Draft Chapter 3 Formulation of Some Elementary Nonlinear Problems Outline 31 Nonlinear Heat Conduction 311 A 1D Example of Nonlinear Heat Conduction 32 Nonlinear Elastostatics 321 A 1D Example of Nonlinear Elastostatics 31 Nonlinear Heat Conduction We begin with a domain Q E anm having the boundary I1 Which is decomposed into two parts m see Fig 311 Pg denotes the region on which temperature essential boundary conditions are speci ed While 17 denotes the region on which heat flux natural boundary conditions are speci ed More general boundary conditions can be considered but these will suf ce for present purposes Ph L Fig 311 Nomenclature for the problem domain 31 Hughes and Ferencz Nonlinear FEB139 Chapter 3 Draft The strong or classical form of the boundary value problem is stated as follows Given the internal heat source f Q IR the temperature boundary conditions g I IR and the heat flux boundary conditions ii Pg lR nd the temperature eld a i2 1R such that gm 2 f on Q heat equation 311 SH u 9 on Pg 312 9239 m FL on ID 313 where u satis es the generalized Fourier law q Qi mug 314 2 ng14 uJ 315 l and n is a unit outward normal vector on T The generalized Fourier law is a constitutive equation which completes S by relating the ux vector q to the temperature u In this example we assume the conductivity K is a nonlinear function of u but that q is still linearly proportional to the temperature gradient Substitution of 315 into 311 defines a nonlinear PDE If the material is inhomogeneous ie the conductivities vary with the position a E Q the constitutive equations are replaced by qi qiuujc 316 m uauj 317 Statement of the weak form corresponding to S requires the de nition of two function spaces The trial solution space S is the set of all functions satisfying any inhomogeneous essential boundary conditions 5 Sulu g on F9 The weighting function space V in turn contains all functions satisfying homogeneous essential boundary conditions on F9 V Vww O on F9 An introductory discussion of the technical conditions 5 and V must satisfy may be found in the Linear Text The weak form can now be written in abstract notation as follows Hughes and Ferencz Nonlinear FEBI Chapter 3 Draft Given f g and r ias in 3 nd u E S such that for all to E V mum w f w a 318 where awuu dz jwjizsijuu j d9 319 9 W 770011 3110 wf 13 fwfdo 3111 2 def 1011 2 fw dl 3112 l Pi Note that the de nition of awu u incorporates the constitutive relation stated in Indeed the third argument in awuu is to indicate the nonlinear dependence of the conductivities upon temperature While the second argument denotes the linear temper ature gradient term This segregation of it into nonlinear and linear arguments of a is not always possible and motivates the introduction of the more general operator notation 77w u Which Will not be a symmetric bilinear form as previously encountered for linear problems We emphasize however that for the formulations considered herein the weak form operators will always be linear with respect to the weighting function 211 The basic idea of the Galerkin formulation is to restate the weak form in the context of nitedimensional subspaces of 5 and V We indicate these subspaces With a superscript h ie 5quot C 5 and Vquot C V The superscript h represents a measure of some characteristic dimension of the finite element mesh employed The nite dimensional subspaces are in some sense parameterized by h A ne mesh would have subspaces with a small h while a coarser mesh would have a comparatively larger h Using the operators de ned in EV the Galerkin problem may be succinctly written 33 Hughes and Ferencz Nonlinear FEM Chapter 3 Draft Given f g and Ass in 5 find it E Sh such that for all wh E 12quot G awh7uhuh whvf wha lr 3113 Where the operators a and 3 f are as defined in W Substitution of expansions for the trial solutions and weighting functions in terms of nite element basis functions and coef cients as was previously done for the linear case7 now leads to the nonlinear algebraic system of equations Md F 3114 where N Md 51 new arts F Fnodal a 1 3 3116 N d represents thermal forces arising from temperature gradients and is a function of the nodal temperatures d The assembly operator A maps the element equation numbers onto the global equation numbers and adds the contributions of n6 to the global vector N Detailed discussion of the assembly operator may be found in the Linear Text In the linear case the element contributions 71 3 simply equals kede F represents heat sources and consists of nodal constributions and sums of element contributions The nodal forces may be considered to be 6 functions or concentrated sources located coincidentally at nodes in the mesh All other sources would enter N through the element contributions Recapping problem N requires the nonlinear vectorvalued function N d to balance the externally supplied heat coming from f and FL as manifest in F Viewing the element contributions to N in component form quot6 71 f3 ff 1 lt a S nan 3117 Where a is the local node number and made qu do 3113 96 34 Hughes and Ferencz Nonlinear FEMquot Chapter 3 Draft BZ Duh2ZBbd an 3119 03 f rVafdQ1Wa dr 3120 96 I The vector Ba is the gradient of the shape function Na so in turn 22 B 5d represents the gradient of the temperature The conductivity matrix D is a function of uh Nedi In the linear case 3120 simpli es to the familiar result nen rigor 2 foDBb an d Z kgbdg 3121 5198 b1 Where It8 is the element contribution to the linear global conductivity matrix Remark In the expansion for the trial solution certain nodal temperatures will be located on F9 ie these values Will be speci ed essential boundary conditions These known values are components of the nodal temperature vector d used in N This contrasts with the linear case Where these essential boundary conditions are immediately segregated from the unknown degrees of freedom and applied to the righthand side of the balance equation through a correction to the element contribution f5 cf Chapter 1 Such a separation is not possible when N d is a nonlinear function of I thus the absence of any correction term in 3121 for g type boundary conditions The segregation of d into active unknown and inactive known degrees of freedom cannot occur until N has been linearized In the present nonlinear context to derive an anolog of the linear conductivity matrix we take the variational derivative of N d quot61 DNd A Dnede 3122 61 which is the assembly of the element tangent matrices ans 8 a lt lt 8 3123 Dn 8 1abn ie the partial derivatives of the element internal force With respect to the nodal tem 35 iiughes and Ferencz Nonlinear FEM Chapter 3 Draft peratures Substituting the de nition of made ane a T en en a z V Bede 2 6 amp aDanJ m win 9e 11 31 and using the chain and product rules for di 39erentiation leads to two terms 5D quotaquot n T Ba may lVdddb Z BedC d9 Qe 111 c1 foDBbdn 3125 28 The second term is reminiscent of the linear conductivity matrix 138 except now the matrix of material conductivities D is dependent on the temperature Given the symmetry of 13 the second term is clearly symmetric The first term is entirely new and it will be shown explicitly in the 1D example to follow that it is nonsymmetric even under the simplest circumstances As discussed in the earlier remark the components of is used to evaluate ned6 and Dnede will consist of both unknown and known nodal temperatures The mapping between global and element degrees of freedom is accomplished through the LilI and I EN arrays de d a 9247 LAIQJ e 0 P LilIQ a e Lila e O A IENa e The de nitions of the L114 and I EN arrays are reviewed in Chapter 1 and fully developed in the Linear Text BL We now have at hand de nitions for all element contributions to the linear matrix problem M MWAdFNMm BL where Ad consists solely of components corresponding to unknown degrees of freedom A Newton type iteration stategy may be constructed around 3 as discussed in Chapter 392 I to obtain a solution to IV 36 Hughes and Ferencz Nonlinear FEM Chapter 3 Draft 311 A 1D Example of Nonlinear Heat Conduction We wish to make the previous discussion more concrete by considering a simple one dimensional heat conduction problem and deriving the element quantities emanating from the choice of piecewise linear shape functions Assume that the problem domain is the unit interval 9 10 l as shown in Fig 312 A temperature boundary condition is imposed at x 1 and a heat ux boundary condition is speci ed at x O The strong form of the nonlinear heat conduction problem is now written g f on Q 101 3128 151 z 9 3129 1q0 5 3130 qiuvu cix Kuauz 3131 The second negative sign in 3130 re ects that the outward normal at x 0 is n 1 see Fig 312 The constitutive equation acknowledges that the domain is inhomogeneous as the conductivity is in part a function of position at Fig 312 Problem domain for the 1D heat conduction example 3 E 0 1 Hughes and Ferencz Nonlinear FEM Chapter 3 Draft The corresponding weak form is Find u E S 9 such that for all w E V wlw1 0 Where 1 W 77wu awu u d f wznuzu72 d2 3133 0 1 def wf fwf do 3134 0 def w p w06 3135 and the Galerkin form follows immediately by restricting W to nitedimensional trial and weighting spaces Find uh E 8quot uhluh1 2 9 such that for all 10quot E Vquot whlwhl O G nwhuh whfwha1 3136 where the operators a 3 and r are as de ned in W We state the element constituents of the nonlinear balance equation N nde Naqua 3137 96 f szafdx 35M 11L 3133 QC Linearization of N leads to a matrix problem A involving the variational derivative of 38 amp Hughes and Ferencz Nonlinear FEM Chapter 3 Draft N d The components of the associated element matrix Dned are 8718 85uhm quotm a Ivar 1 AN 1 c e 39 39 ad au 6 EN dc dm 3 i 39 QC c1 NawuhxNtzdx 3140 90 The respective nonsymmetry and symmetry of these two terms will be explicitly illustrated by the following example 0 Twonode Linear Element We now Will evaluate the previously presented element matrices in detail for a par ticular choice of element shape functions the familiar piecewiselinear hat functions pictured in Fig 313 along with the element s parent domain 5 E 1ll N2 Fig 313 Parent domain and linear shape functions for a two node element To evaluate the integrals we will employ one point Gaussian quadrature An integral is rst transformed to the parent domain 1 dxjd 93 1 3141 where j 2 235 is the 1D version of the jacobian determinant The integral is then approx imated by sampling the integrand at the midpoint of the element and scaling it with the 39 D Hughes and Ferencz Nonlinear F 133 appropriate weighting factor VV2 1 jd 3 W 3150 1 For this simple element the jacobian is a constant j lac2 therefore dm z 1711 50 Q6 The shape functions and their derivatives may be concisely written 1 1 I a Na 31 tial 31Tquot1l l5 1 2 Then from the isoparametric representation a 2 221 Nam it follows Va g Z 392 15 Z Nay xz 21 6 6 e 2 1 2 2 so by the chain rule 1 a Ivan lva5x5 1 he The components of the internal force are de ned nd eNa qua 2 he Naizq 0 25 423 117 40 1 10 Expanding the arguments of the flux term qltogt quhlt08gtufoegt am 7137 nen nan q lVadZ Z Nbadi Z N i 61 Cl 111 3110 Chapter 3 Draft 3142 3143 3144 3145 3146 3147 3148 3149 39Hughes and Ferencz Nonlinear FEM Chapter 3 Draft d8 l d de de 8 e q1 2h61712gt 3153 Upon evaluation 10 can be used with 3150 to yield 6 e 1 n d q0 1 3104 ie the ux is equal and opposite at the two nodes A convenient approximation typically employed is to interpolate the source function by its nodal values7 ie 2 f e fh 2M f 3155 61 Where ff f Assuming the FL boundary condition is not active on this element f Nafda z 1Vafh dcc 3156 m m a he Nam Loft 3157 therefore 6 he 1 1 ff f z l Mfg 318 3159 In essence the nodal values of f are averaged and placed at each of the nodes The components of the tangent matrix may be calculated as e h quotan and an u ad 2 he Nagma If m Nczd 60 he ayzs uhxe Nm 50 3160 heel 6 Mama 3 1 a he an 2 he 81a e e 1 W he e W0 MU D be 3151 311 3 Hughes and Ferencz Nonlinear FEM Chapter 3 Draft or in matrix form 6 c quot391 1 dg dflz BMW 5327b 1 1 1 1 0 1 1 The two terms arise from the two dependencies on temperature as previously discussed 3162 Clearly the rst term is nonsymmetric and the second is symmetric The rst term involves the derivative of the conductivity with respect to temperature as scaled by the temperature gradient across the element The rst term will collapse to zero only if at least one of these two derivatives is zero otherwise D776 and the assembled DN Will be nonsymmetric 312 Hughes and Perencz Nonlinear FEDI Chapter 3 Draft 32 Nonlinear Elastostatics A discussion of the general theory of nonlinear elastostatics will be deferred until the next chapter For the present we limit ourselves to a restricted case which is nonethe less illustrative of the finite element techniques used in the general case without being encumbered by the mechanics required to treat the full theory Let us make the following assumptions 1 The deformations in particular strains are small 2 The body is initially unstressed 3 Material behavior is the only source of nonlinearities The first assumption is reminiscent of the linear theory but we now admit with the third assumption that the stress is a nonlinear function of the linear strains through some generalized Hooke s law We assume the same problem domain considered for the nonlinear heat conduction problem and illustrated in Fig 311 The strong form of the boundary value problem for an unspeci ed number of space dimensions nsd may be stated I Given the body force distribution f 2 1R7 the displacement bound ary conditions g I an d and the boundary tractions 17 1R7 nd the displacement eld u 5 gt n and the stress or 52 151quot X an d such that S lt 0ng 1 0 on Q 321 u r g on l g 322 at nj a on 13 323 These equations are form identical to the linear case cf Chapter 1 The difference be tween these theories will arise from our choice of a constitutive equation relating the deformations u and the stress a In particular we assume the body is hyperelastic ie that there exists a given strain energy density 5 3 such that 8 2 72 8515 3 4 The applicable strain measure to be employed is the symmetric part of the deformation 313 leughes and Ferencz Nonlinear FEM Chapter 3 Draft gradient d f 1 5 2 3013 uji 320 The material moduli are in turn de ned to be def 9033 52 2 2 326 6 96 653585 where for a nonlinear inhomogeneous body ngH ciJkz 52 In the linear case I has the simple quadratic form 1 P 195 239Cijk15ij5kl 321 The symmetry of mixed partial derivatives and the de nition 326 immediately implies the major symmetry of the material moduli tensor Cijkz ckuj The geometrical inter pretation of c is that it provides tangent information about the stressstrain diagram as shown in Fig 321 0a Fig 321 Example of a nonlinear elastic stressstrain constitutive relation The weak form corresponding to S is known as the principle of virtual work in the mechanics literature Again we must de ne the trial solution space 5 whose members satisfy the essential boundary conditions on F9 and the weighting function space W whose members satisfy homogeneous essential boundary conditions on F9 The Weak form is then written 31 Hughes and Ferencz Nonlinear FERI Chapter 3 Draft Given 1 g and Ii as in 5 nd u E S g on F9 such that for all wEV wlwizOonl g mm w J w or me Where def nwu 2 10 3031 6 d9 329 W j 9 I 11 f 2 wifidQ 3210 2 def ma wiadr 3211 L 1391 Again the operator n is linear in the weighting function 10 but nonlinear in to via the consitutive relation de ning 05 The notation wad denotes the symmetric part of the gradient of the weighting function This restriction can be applied because the stress is symmetric therefore wijaij wad0 The body force and traction boundary condition terms are familiar from the linear case The Galerkin approximation to W derives from restricting the trial solution and weighting function spaces to nitedimensional subspaces 5quot C S and Vquot C V Given f g and 5 as in 539 find 21quot E 5quot 92 on Pg such that for all toquot E Vquot z 0 on F9 G nwhuhgt w f w a 3212 Where the operators 7 0 and 1 are as de ned in HIV 4 Hughes and Ferencz Nonlinear FEIV39I Chapter 3 Draft The nonlinear algebraic equation N resulting from G has the standard form Fint a Nd F 5 Fe 3213 where N Md 51 new 3214 F Fmdai El 3quot 3 3215 The element contributions to N may be expressed in component form 71646 712015 fa f5 1 S p lt nee 3216 The component index p is the local equation number p neda 1 i de ned by the local node number a and the direction number i The parameter ned represents the number of element degrees of freedom per node which for the present case would equal nsd Then the element internal force is written nude ifEf e d9 3217 98 which in the linear case only reduces to e83D Z Bbdz d9 J BIDS d9 3218 9e bzl 52 In this linear case the constitutive matrix D is a constant The notation of 3217 and 3218 are made more precise by de ning its constituents assuming the case of plane strain The vector 8 is simply a cartesian basis vector eg wt mm The relevant components of the stress and strain tensors are arranged in a convenient 0 11 511 0 05 0 22 8 Z 522 0 12 2512 3219 vector form 3220 and the straindisplacement matrix corresponding to 5 is 316 I A Hughes and Ferencz Nonlinear FEBI Chapter 3 Draft man The conventions for a e and B in the instance of plane stress torsionless axisyrnmetry or three dimensions are fully detailed in Chapter 2 of the Linear T ext The components of the element external force are identical to the linear case zM mMan 8 Qc TI and consist of contributions from the body forces and surface tractions Solution of N requires its linearization with respect to the unknown nodal displace ments A consistent linearization gives rise to the tangent matrix net DNd A Dnede 3223 621 which is assembled from the element contributions 39 8776 ang e Ii m 4 M in a lt3 2gt Where the second component index q 2 med 1 1 j Taking the derivative of the internal force ans 8 T T 2 8d adgb e jBaa39e d9 32 5 J 96 we note that only the strains are a function of the nodal displacements e ET 6 as do 3226 1 a 99 Substituting the expansion for the strain vector if d EB e 89m 8 2 B e 2 C 7 c1 C C c1 dgc 31 A J Hughes and Ferencz Nonlinear FEBI Chapter 3 Draft and then using the chain rule ne and a quotequot 3 T T 6 22 ad 61 Ba as adj ZBCdC an 3 a Q cl 83quotfo Qf gjiwc aend 3229 Qe alt33135 BbdQ ej 3230 96 This element tangent matrix is formidentical to the linear case however now the coeffi cients of the constitutive matrix 001 551 must be evaluated for the current strain state The contraction between the J and ijkl DIJ Cijkl5 subscripts is discussed in the Linear Text Remark It is worthwhile to underscore the contrast between the tangent matrices aris ing in nonlinear heat conduction and nonlinear elasticity As just seen in elasticity the tangent matrix is formidentical to the linear stiffness and hence typically possesses the same properties of symmetry and positive de niteness As shown earlier for the simplest example of nonlinear heat conduction there the resulting tangent consists of a symmetric part which is formddentical to the linear conductivity matrix plus a second nonsymmet ric contribution Thus even for familiar continuum theories the admission of nonlinearity can lead to nonsymmetric tangent matrices and this possibility must be examined on a casebycase basis I The evaluation of nede and Dnede requires the localization of both active degrees of freedoms for Which the displacements are being approximated as well as essential bound ary conditions which are speci ed This mapping from global to local degrees of freedom LAIia e 774 0 P LiMz39a e is written d6 dP 3 9 31 P T gig LMia e O A IENa e quotquot The reader may consult Chapter 1 or the Linear T ext for background on these mapping arrays We have now defined all element contributions to the linear algebraic problem 31 needed for a Newton like solution of the nonlinear algebraic problem 318 AHughes and Ferencz Nonlinear FEW Chapter 3 Draft 1D Example of Nonlinear Elastostatics We now examine a simple onedimensional elasticity problem and derive the element quantities arising from the use of piecewiselinear shape functions The problem domain is the unit interval 9 10 1 previously illustrated in Fig 312 A displacement boundary condition is enforced at a 1 while the traction is speci ed at x 0 The strong form of this nonlinear elasticity problem is aquot f 0 on Q 10 ll 3232 539 u1 9 3233 0390 5 3234 Physically these equations govern the longitudinal deformations of a rod or the trans verse deformations of a string This strong form is completed by the speci cation of a scalar constitutive relationship a 05 2 0uz Where the nonlinearity arises because now do C C 3235 d5 5 lt gt ie the modulus is a function of the current strain state The associated weal form is f Find it E S 9 such that for all 212 E V wlw1 O 77w7u w f w r 3236 where 1 W lt 77wu d ws due dx 3237 0 1 unf d wf d2 3238 1 L w rap 2 want 3239 339 Hughes and Ferencz Nonlinear FERI Chapter 3 Draft and the Galerkin approximation follows immediately by restricting EV to nite dimensional trial and weighting spaces as denoted by a superscript h Find uh E Squot uhluhU g such that for all inquot E 1 quot whlth 0 may 2 Muff wk 4 3240 Where the operators 77 O and 3 p are as defined in WV The element contributions to the nonlinear algebraic problem N may be Written much more simply because the tensor indices contract out of the one dimensional problem eg the straindisplacement matrix simpli es to N Explicating ni de 2 lez 0uf zdx 3241 QC 2 Nma Zdeg d1 3242 95 61 and f3 Nf dx al 61 3943 P a z 0 e gt 1 39w 9 The components of an element s contribution to the variational derivative DNd are ne 2 V NI ad 1 ce 5 dx 3244 QC Again this is formddentical to the linear case but the modulus is now a function of the current strain 0 Two node Linear Elements To evaluate the previously de ned element quantities for the case of piecewiselinear shape functions see Fig 313 we will again employ onespoint Gauss quadrature to com pute the integrals Hence 7 fNa zadx 3245 QC 320 7s qughes and Ferencz Nonlinear FEEI Chapter 3 Draft a he i 0 a go 3246 2 he Naz03 U50 324 he 0uyz06 3248 Computing the strain at the center of the element quoten 3 e 612 d z a V edZ 2 u 0 6421 0 he 3 49 then 6 e de de 1 n d a 1 3250 The element external force is identical that of the heat conduction case e he ff ff 1 A f 39339 1 681 kD 3201 Finally we have the consistent tangent 871a adb 2 NM CNM d9 3252 QC z he 139 a chJ 3253 e 1039 e 1b h he ce0 he 3204 So in matrix form d5 d 1 1 1 C 7239 21 1 which is familiar from the linear case cf Chapter 2 of the Linear Text Chapter 4 Finite Deformation HyperElastostatics 41 Continuum Mechanics Background As a point of departure recall the strong form of the boundaryvalue problem of linear abstract notation elastostatics Find the displacement vector u 51lch that diva f 0 in Q equilibium equations 5 u g on l g 41 t 2 039 n 2 h on I The stress tensor 039 is de ned in terms of the strain tensor 6 by a constitutive equation the generalized Hooke s law Gij CWW 42 1 i 2 CHAPTER 4 FINITE DEFORMATION HYPERELASTOSTATICS lrepeated indices imply summation unless otherwise indicated where Cijk GigMac are the elastic coefficients and 6k 7101 43 This is a special restricted version of smalldeforrnation nonlinear elasticity Where there exists a given potential 03 the strain energy density function such that 91 I r 44 a 5g k l and the tangent moduli Cijk GigHOB 6 form its Hessian matrix 80 Cijkl max 4 5 862566 The major symmetry of the tangent moduli thus follows from the existence of a strain energy function The weak formulation of the boundaryvalue problem sometimes referred to as the prin ciple of Virtual work is W Find u E S satisfying essential boundary conditions such that for all m 6 V satisfying the homogeneous counterpart of essential boundary conditions QuitJog A wgfi A wihi Where in the nonlinear case mj is no longer defined by the generalized Hooke s law 42 CONTINCCM MECHANICS BACKGROUND 3 Figure 41 1 Fig 12 p 3 In the case of nite deformations e is not a bona de nonlinear strain measure To Linden stand the situation we need to introduce some elementary notions from nonlinear continuum mechanics The following employs notation from A reference con guration 13 is chosen for the body Points in B called material points are described in Lagrangian coordinates X X1 X2 X3 A con guration of 8 is a mapping qsz gt R3 representing a deformed state of the body Points in R3 called spatial points i a are described in Eulerian coordinates 93 z 331 12 33 Parameterizing the con guration by time t leads to a motion a X i MX 47 Where the second line refers to t xed We assume that the initial con guration coincides with the reference con guration ie 450 1 48 so that X 450Xl 49 4 CHAPTER 4 FINITE DEFRMATION HYPER EELASTOSTATICS The displacement of material point X is defined by UdX MIC X 410 W a The initial displacement is thus Uolxl QMX X X 2 0 411 In dynamic applications often U0 Og ie the initial con guration differs from the reference con guration A general rigid motion of B C R3 is a3 Xt 2 max C25 412 Where Q is a proper orthogonal matrix representing a rigid rotation 6 0 for such a motion thus disqualifying it as a measure of strain and c is a rigid translation Kinematics The components of the motion are set aiXit 413 The matrix of partial derivatives of 5 is called the deformation gradient 414 EX w 8 15 m 36239 F iaXIl 4 Z CONTIN l M ECH ANI CS BA CKGRO UNI 5 i Lower case indices i j k and 1 denote spatial quantities eg a Ii and capital I J K and L are used for material quantities for example X XI The material time derivative denoted by a superposed dot is the derivative with respect to 73 holding X xed Thus the velocity of a material point X is WX z X zg mm regarded as a vector emanating from point 13 The material acceleration of a motion A 2 AlXt is AV2 aw The spatial velocity of the motion is Magi VlX t where m X t We may thus de ne v by v tXt VXt 417 ie vt t 2 VJ which de nes the composition TM 0 I5 Vt 418 v is the Eulerian representation of the velocity and V is its Lagrangian representation To obtain the spatial velocity gradient V we rst consider the material derivative of the 6 CHAPTER 4 FINITE DEFORMATION HYPER ELASTOSTATECS deformation gradient By 418 employing the chain rule inwg ax EEaX av 31 This may be written in components as V211 MOM IMF 3391 with implied summation on the index j Since z and 3 2 V W FF 1 assuming F 1 exists 419 420 421 41 COLV39TINCUM MECHANICS BACKGROUND The Sacobian determinant of the motion is J det F 423 9 72 Note that X If gt 0 corresponding to impenetrability a property of real motions The material or Lagrangian Green Lagmnge strain tensor is de ned by 2E 2 PW 1 424 so that E 2 ET and in components 21911 FHEJ 6J 425 In order to express the Green strain in terms of displaCements recall 8gb F 2 8X 2 X U X 8Xlt lt gtgt lt4 26gt where the second line follows from the de nition of the displacement 410 In components 9 Xi 1 8X T 1239 6 221 427 8 CHAPTER 4 FINITE BEFORAXE ATION HYPERELASTQSTATICS Figure 42 Current con guration with outward normal since Xi ogJXJ Components of the Green strain are thus 21511 H 551 UiJXCSiJ Um 517 H 617 UL U1 UMUN 5U so that 1 1 EN 5Uw in EZUMUN With implied summation on i Stress The rst HoleKirchhoff stress tensor is P 1ch4 or in components 8X1 Pi Jm v I J Where J 2 det F and the third term is FquotT I X1 339 428 A 429 i 430 431 The symmetric true Cauchy stress tensor 0 is uniquely determined from Cauchy trace tion vector t 2 an which represents the force per unit current area exerted on a surface element with outward normal n Consider an in nitesimal volume element in the current 4 E CONTIN U 5M MECHANICS BA CKGRO END 9 Figure 43 Mapping of reference con guration with N to current with n con guration dt 2 tdV not to be confused with velocities Similarly an in nitesimal surface element in the current con guration is do 2 gthl not to be confused with ac celeration We claim that there exists a rst PlOlaijlI Chl39lO traction vector T 2 PN representing force per unit reference area such that t da TdA 432 so that the two traction vectors are parallel To see this recall from calculus and differential geometry that for a volume element dv 1 IV 433 and for a surface element nda JF TN dA 434 Therefore 0n da a JF TNdA JUF TNdA H H PN M 435 10 CHAPTER FINITE DEFORMATION HiPER ELASTOSTATICS Where the last line follows from the de nition of the rst KolaKirchhoff stress tensor 430 or in components min da 2 PiINI We have thus Veri ed 432 Consider integration over a closed surface For the rst Piola Kirchhoff traction vector jndA z jHIJVIdA A A 2 PHJdV V Where the second line follows from the divergence theorem Likewise for the Cauchy traction vector ftgda awrydd 039ij3j d U 102 ij H H Where the last line follows from the change of variables These integral relations hold for arbitrary closed surfaces which are the boundaries of subsets of v and V such as U C V and 411 C 2 Thus lt2quotan Pm 439 Remark P is not generally symmetric Returning to the concept of strain energy density a given function of Green strain 42 CONTINE39UM MECHANICS BACKGROUND 11 Figure 44 Mapping of closed surface A in 23 to a in 343 we have from thermodynamic considerations a symmetric stress tensor such that 8Q s 5E 440 called the second Pioletrchzo stress tensor In components this material tensor is 8G 41 SH aEU 4 l The relation to the previously de ned stress tensors is P 2 F3 2 JaF T 442 01 039 J1PFT z J 1FSFT 443 The material tangent moduli are de ned by 8311 u 821 CIJKLE 8EKL 444 representing the slopes of the SE diagram The major symmetry of the tangent moduli thus follows from the existence of the strain energy function and the minor symmetries from 12 CHAPTER 4 FINITE DEFORMATION HYPERjELASTOSTATICS 39 the symmetry of Green strain Similarly the analog of the Cijkl in the smalldeformation nonlinear theory considered previously is CijkllE 3 J1EIFJFkKFlLCIJKL 445 but there will be another term in the consistent tangent to be seen later 42 Variational Weightedresidual Formulations The weak form survives intact widowd1 fwifi wihida 39U U 1h where v QSJV lls the role of Q the current con guration but the tangent is calculated from the material description which hence must be derived Consider the change of variables we max WXX 447 which is thus de ned by the composition VVZ 102 0 QB 448 42 VHRIATIONAL WEIGHTEDRESID U53 L FORMULATIONS l3 Figure 45 Mapping of boundary partition into Ag and AH toog and 15 We assume that the boundary conditions can be stated in terms of material coordinates and rede ned as follows The displacement boundary condition7 which is built into the de nition of the set of trial solutions u 2 g on a9 449 Where ag iAg lls the role of F9 is equivalent to U 2 g on Ag 450 Similarly the traction boundary condition t h on ah 451 Where ah tAH lls the role of Th is equivalent to T 2 H on AH 452 Note that t T since tda TdA Weighting functions satisfy the homogeneous analog of 14 CHAPTER 4 FINITE DEFORMATION HYPERELASTOSTATICS Figure 46 Mapping of L C B to gm essential boundary conditions 100 onag and W20 onAg Due to the change of variables integration becomes jgdvVJdv and for surface integration mida F39T 39 NIdA a A lt gt21 453 454 I 455 456 We also need to employ conservation of mass Consider a set of material points 2 C B d E pg 2 3 0 754m 4 2 iCL XRIATIONAL W El GHTED RESID ML FORMULATIONS 15 iwhere p w is the mass density of the deformed body at time t integrating p52 dvtz pm 61722151 458 45W calm Taking 752 t and t1 2 0 corresponding to the reference con guration and noting that duo JG dV dV since J0 2 det 2 1 leads to ptdvt fpodV 459 43w 1 Employing the formula for change of variables dvt Jt dV we obtain Apnth Up0dl 460 This holds for any set of material points 2 leading to the local form of conservation of mass tht Z 00 461 in terms of the mass density of the reference con guration pgX The material description of 446 requires that the integrand on the left hand side be written in a more useful form Recall that the material description of weighting functions is de ned by the composition 448 Thus 8W 0 a lwiod 8w 8a 462 16 CHAPTER 5 FiNITE DEFORMATION HYPERJELASTQSTATICS where the second line follows from the chain rule in abstract notation we write GRADW grad39w F or gradw GRADWF1 which in components is 811 2 We 463 Returning to the weak form 446 by change of variables the left hand side becomes fvwmaijdv sz JO39z delf 8X1 aquot d An axjaj1 V fi l39i IPudV V V smug dV 464 Where the third and fourth lines follow from the de nitions of the rst and second Piola Kirchhoff stress tensors 430 and 442 respectively The lefthand side of the weak form is left in this form at this stage We now turn to the righthand side of the weak form 446 The body force per unit current volume f can be expressed as f 2 pb Where b is the body force per unit mass eg gravity We express the body force B 2 b o I as a function of X and assume that z b o ie B does not depend on the motion this is referred to as dead loading The contribution of the distributed loading is obtained by the change of variables H 13quotquotde LwWX pb XJdV H V W 300 dV 465 42 tiaRIATIONAL WEIGHTEDRESIDUtah FORMULATIONS l since oj 2 p0 by conservation of mass 461 Note that nothing in the second line depends on the motion 1 The surface traction term is treated in similar fashion Recall the rst PiolavKirchho quot traction vector T was de ned so that tda 2 TdA and so for the traction boundary condition it da 2 HdA The contribution of the traction boundary condition is Lhwihida Li ill iHidA 466 Where we again assume no dependence of the loading on the motion Hz HAX ie no follower forces etc Employing equations 464 466 to the residual ie the balance of internal and exter nal forces of the Weak form 446 o H dv wihida wg j0 gj do 11 oh 1 A manaw Wimda WiJPZIdV 467 f I Where v 2 2 Note that the two descriptions are identical Discretizing in the usual manner W39ZWX EGANAX and ZciANAm or locally lV X Z chiVAX and Z ciaNaz Within each element We need to calculate ti and then E1313 511 SIAE 039 and P and CJKLE By change of variables to the parent domain element integrals become fuemdvzfamjd 468 18 CHAPTER 4 FINITE DEFORMATION HYPERELASTOSTATICS Figure 47 Mapping of to ue by q Jh where 226 2 Q6 and j det and in the material description mdvzjmjodtj 469 V5 I where jg det Note that V8 and v8 are sets of the same material points This is the Lagrangian description of the motion in terms of sets of material points In other applications such as fluid dynamics the Eulerian description is used in terms of sets of spatial points T we alternative descriptions of motion are the updated Lagrangian Where integrals are calculated in the current con guration and the total Lagrangian Where integrals are calculated in the reference con guration Both are identical 43 Linearized Operators The variational equations are linearized to obtain the left hand side operator The consistent tangent 8N K 2 83 470 is obtained by discretization prior to linearization Alternatively the variational equations can be directly differentiated assuming the constitutive algorithm is exact which is true in 43 LINEARIZED OPERATORS lg elasticity The left hand side is considered a function of the motion H fl l 10143094 do L39 WUsz at V H H V WMFz JSJi dV BWAXWMX 7 I v ax WSMEWWV an Where v 2 Q and V 00 and the Green strain E depends on the motion by its de nition 424 and the de nition of the deformation gradient F 414 Linearizing this function eav l D iAU an lezo Where 5 is a real parameter and the right hand side represents the action of a linear operator on the displacement increment AU de ned by the change of variables AUX Auq5m or the composition AU Au 0 d Discretizing the displacement increment as ZNaAdZa Within each element leads to K dAd Linearizing the left hand side of the variational equation 473 520 dV 50 n sav atAUz39l quot V W47 aXJ i f quot d W1I W8XJ 511E eAUDdVU 39 ast EM N titan an lt d5 Q 3920 CHAPTER 4 FINITE DEFORMATION HYPER ELASTOSTATICS Where lL CHIS L Employing the de nition of Green strain 35KL 8 5AU T a 5AU 2E 5AUlt Wax gt Max 1 leads to d2E eAU aw T93 93 Tam d5 50 8X 8X ax 8X and in components 2 KL AUiI39 Llt5iKAUiL d5 50 H AUM39FM FiKAUiL Thus d A 1 EM 211 50 H V VVz JAUiJSJI I ViJFiJCJIKLAUjKFjL 03V 474 475 1 V mIAUiJSJI WgFiJC JIKL AUjKFjL Fj t lUggd dV 477 Where the second line follows from the minor symmetry on indices KL of the tangent moduli Pushing forward to the current con guration we employ do J dV Where J 2 det F Wit1 me Where F g xle and AU 2 ALIMP The rst term is LEVEVIAULJSdeV 1 A U i ijAui kaJSJdV H 1 J39wi ijjSJIijAugkj d H wijaijuzvk do I 478 p LINEARIZED OPERATORS 21 i Where the last line follows from 443 This is thus an initial stress contribution The second term is fvVVLIFiJC39JIKLAiBKFjLdV Wier JFjLC JrKLAUjjrCW H l fwiJcFkIFiJFjLCJIKLAujJFZKj dv l fwjiFiIFjJELFkKCJIKLAUlk j d0 ulwancz jkgz ua m where dummy indices are exchanged on the third line and the last line is obtained by the de nition of the spatial tangent moduli 445 and their symmetries This term is similar to the one leading to stiffness in the linear theory giving rise to the standard f BTDB The straindisplacement matrix B is identical to the one in linear theory The matrix D is derived from CMMVF Combining these results and exploiting minor symmetries of the second term yields 1 Efltb 6 AU 2 wi jajkAui k wig0251611321 d1 50 The second expression may be written in terms of strain due to the minor symmetries and the major symmetry leads to symmetry of the stianess matrix We Wish to investigate the symmetries of the rst term starting with the major one For this purpose we write wiJO AuiJ 23 039jz6akl uzci 481 J N CHAPTER 4 FINITE DEFORMATION HYPERJELASTOSTATICS The expression in the parentheses represents efiective moduli due to current stress leading to initial stress stiffnessa as mentioned The rst pair of indices ij may indeed be replaced by the second pair 393 in this expression since O jgdik Z O39gjdk so that major symmetry does exist and symmetry of the stiffness matrix follows This is expected by the existence of a potential the strain energy density function Looking at minor symmetries we need to ascertain whether i and j can be switched and Whether k and l can be switched However 0j15ki Uij kj 483 O39jl ki O39J MSH Consider for example i z k 2 1 and j 2 l 2 This leads to 611 1 on the left hand side in both cases whereas on the righthand side we have 612 621 0 so that equality cannot hold This expression does not possess minor symmetries so that the additional term depends on all the information contained in the deformation gradient Fa and cannot be written in terms of strain the symmetric part alone To demonstrate this point consider the polar decomposition 1 pp 5155l of the deformation gradient F RU 484 43 LINEARIZED OPERATGRS Where R is a proper orthogonal tensor ie RTR RRT 1 and det R 2 1 representing local rotation and U is a symmetric positive definite tensor representing local stretching not to be confused With displacement of material points FTF z RUFRU UTRTRU z UTU U2 485 For the Lagrangian strain measure E 213 FTF 1 Uquot2 1 486 so that it does not contain information related to rotation The linearized lefthand side in its spatial and material representations is w5j6ik039jl CHICOAWJCZV L mJ6ikSJL FiIFkKCIJKLAUkL dV 487 Both versions are identical The expressions in the parentheses represent the effective moduli which due to the rst term in each which leads to initialstress stiffness do not possess minor symmetries 24 CHAPTER 4 FINITE DEFORMATION HYPERJELASTOSTATECS 44 Discretization Coding As previously mentioned the term containing material moduli fu w jCijklAukyf in leads to the standard stillness term fut BTDB dc The strain displacement matrix B is identical to the one in the case of small deformations The matrix D is derived from cijde The geometric stiiiness term which results from initial stresses and may not always be included in the tangent stiffness W 435 lt1 E M to 514 fw j ika zhuzzdv Z Ci 4ixrAj6ikUjZAdkBJNTB Zdv leads to the following contribution to the element stiffness matrix kieakb Naj6tk0ijbzdv 489 in which the Kronecker delta can be taken outside of the integration Consider for example a fournoded linear tetrahedron in three dimensional analysis Fig 48 with three degrees of freedom per node for a total of 12 element equations The 12 X 12 element stiffness matrix is made up of 4 X 4 nodal blocks of size 3 x 3 each Due to the Kronecker delta initial stresses contribute only diagonal terms in each nodal block in z fNaj0ijbzdv no sum on i 16 H LVNaTaV1Vbjdij 490 where j 2 det 33 In each block these three diagonal entries are equal This is similar 44 DISCRETIZATION 3 31 Figure 48 Linear tetrahedron Figure 49 Block diagonal contribution of initial stress to stiffness to the case of heat conduction The stress matrix 039 need not be positive with the effect of destabilizing the stiffness as in the case of buckling under compressive loads The small deformation nonlinear case may be recaptured from this general nonlinear formulation by 1 Not adding displacements to the initial coordinates ie not updating the geometry 2 Not adding the quadratic terms in the strain calculation 3 Not computing and accounting for the initial stress stiffness Furthermore the linear elasticity case may be similarly obtained if in addition the elastic coef cients CW are xed as constants To recapitulate the nite deformation case mffil X czin 491 Where X 30 is used to de ne the current element geometry Displacements are updated 26 CHAPTER 4 FINITE DEFORMATION HYPER ELASTOSTATICS KM z Fm Nd jl1 492 Where the righthand side fag Nafi dvfai Na h da ezr fue 330 do is identical to the small deformation case and the left hand side stiffness K 2 Kc Ka is obtained by contributions from material moduli stiffness and initial stress stiffness Element arrays Given d311 and X 7 the above quantities are localized Shape functions and derivatives are calculated at the quadrature points Derivatives Nay are needed for spatial integrals containing Ba To calculate 0 derivatives 1 a are needed for F since 92 8X 8n 1 l 493 Where displacements are approximated in the usual way inside an element new uh Z Nada 494 a1 The Jacobian determinants jg gg and j 2 are also obtained from these shape function derivatives The calculation of F is thus completed Next the Green strain is calculated E FTF 1 495 It 44 DISCRETIZATION 2 The second PioiaKirchhoff stress tensor 39 81 4 96 8E r l is calculated by a formula that is encoded and from it Cauchy stress 1 T 7 Where J z det F The material moduli C are also obtained by an encoded formula and then c is calculated from 1 02w J iIFjJFka FlLCIJKL 498 Principle of Minimum Potential Energy Consider the functional of a motion 45 under dead loads ma 2 A ltIE dV V as adv A H qbdA 499 H Where I is the strain energy density per unit undeformed volume B is the body force per unit mass in the initial con guration and 71 is the traction per unit undeformed reference surface area Consider variations of 475 ie weighting functions 1 that satisfy IV 2 0 on Ag The stationary value is obtained by forming H 5W Where 5 is a real parameter and taking 3928 CHAPTER 4 FINITE DEFORMATION HYPERELASTOSTATICS the variational derivative see 2 p 188 0 2 611 up A DIHng39W d H W d5 50 H ltL El EWdV Vp03 5WdV LH H 5WdAgtgt 81gtE 8E 5W 8E1 i 85 W 31 50 dv VpOBWdt AHHWdA 4100 50 Recall from the calculation of the linear operator BE A 1 z EmmaJ AuwFu 4101 0 and by symmetry of this term the stationary value is found from SUMile av pr de f H WdA 0 4102 V V AH This is the equation of Virtual work or the weak form in the material con guration For veri cation push forward to the current con guration employing W7 2 31quotng and H dA hda mm 531511307 dv fpb39U dv quot h wda 0 4103 U U 3h 0395 and since p pgJ the spatial version is obtained Remark The Virtual work equation f 0 is a nonlinear variational equation that 44 DISCRETIZATION 29 is written in this form for linearization by taking q q5 5W 1W5 DEW W 4104 is linear in W Df Au 2 D2Hq W Au 4105 is linear in both W and Au and is symmetric by the symmetry of the second derivative leading to K 2 KT 30 CHAPTER 4 FINITE DEFORMATION HYPER EELASTOSTATICS Bibliography 1 JE Marsden and TJR Hughes Mathematical Foundations of Elasticity PrenticeHall Englewood Cliffs New Jersey7 1983 2 TJR Hughes The Finite Element Method Linear Static and Dynamic Finite Element Analysis Prentice Hall7 Englewood Cliffs New Jersey 1987 31 Chapter 5 Finite Deformation Elastodynamics 51 Continuum Mechanics The strong form of the boundaryvalue problem of linear elastostatics O divaf inQ u g on TE 5391 tzzc39rn h onlh survived in its entirety the generalization to the spatial description of nitedeformation elasticity in the domain 2 replacing Q with f pb l 39 2 CHAPTER 5 FINITE DEFORMATION ELASTODYNAMICS The equations of linear elastodynamics are pa divo f in Qx0T 39u 2 g on ngl0T cr n h on PhXOT 52 um00 u0m0 in Q 11a300 2 210830 in Q where a superposed dot denotes the material time derivative the derivative with respect to 15 holding X xed and mo 2 X These equations too survive the generalization to the spatial description of nite deformation but for the initial conditions V replaces Q Recall the material displacement velocity and accelerations U u o qb V 39v o 5 53 A a 0 gb Where V 2 E 11 since at qbt 1 Thus the material and spatial descriptions of acceleration are A V U 54 152 SEMLDISCRETE lHRIATIONAL WElGHTEDRESID HAL FORMLEATIONS 3 Remark Recall that 1123 1 2 v X t t By the chain rule the acceleration is given by 817 3v 1 on since vi 2 This is the form often used in fluid dynamics The rst term on the right hand side represents convective acceleration In the spatial description pi 2 divo39 pb 56 Recall that DlVP Jdiva39 or7 in components PM Jam Substituting in the equation of motion pa 2 DIVP pb 57 since it 2 a and multiplying by J yields the material description mAzDWPmB ma Where p0 Jp 52 Semidiscrete Variational Weightedresidual For mulations Recall that the variational equations of dynamics are constructed by taking the static version and invoking D Alambert s principle namely pb lt pb pa in the spatial description and in the material description p03 4 003 pOA Thus the spatial description of the weak form 4 CHAPTER 5 FINITE DEFORMATION ELASTODYNARHCS fwgpagdv39zvgjaij dv fwipbidv f f wihida 59 U 1 v ah which in the static case no acceleration after discretizatiou leads to Nd F and the material description is V WipoAi W WMPH dv jv WipOBz czv AH new dA 510 Where p0 pJ The inertial term uwipiiidv V Wzpo ridV 511 is semidiscretiZed in the usual way by employing WT 2 W39ZWX and U UHX 315 Ah w phizg M W wiped M 512 in the integral over the xed domain Vh Substituting the standard nitedimensional ex pansions W Z CiANAX A 513 72h 2 Z di1itNAX 1 52 SEMI DISCRETE VARIATIONAL WEIGHTED RESID UAL FORM TLATIGNS 5 into the approximation of the inertia 2 lV hpg ijl div 2 2 C2 6 V pervpvt dV 333 as v AItAJB leading to a constant in time mass matrix M 2 illmm since p0 p0Xl which is identical to that of the linear case The inertial forces are added to the balance of forces from nonlinear statics yielding a system of nonlinear ordinary differential equations Ma Nd Fm 515 with initial conditions as before 10 2 do 516 10 2 do which de nes the matrix problem for fully nonlinear elastodynamics 6 CHAPTER 5 FINITE DEFORMATION ELASTODYNAMICS 53 Classical Timestepping algorithms The next step is to discretize with respect to time Recall Newmarkgs family of algorithms 1 p 491 with predictors At dn1 Z dn 39i Atvn quotiquot 17n1 390 m Where At is the time step 3 and 397 are free parameters and dn um and an are approximations of dam and El n respectively and correctors dn H Z ani 53L5t20n1 518 quotn1 ni7Atani Where the predictors are known from the previous step and an1 is calculated from the equation of motion Man1 2 F221 Where F 2 FeXtUn in which the only di erence from the linear case is the term Ndn1 This nonlinear algebraic problem is solved in each time step To do this an iterative strat egy modeled on the predictor multicorrector algorithms described in 1 p 562 is introduced CLASSICAL TIMESTEPPING ALGORITHMS I First the iteration counter is initialized Then there is a predictor phase 1534 ginH vn1 nd1 C D pa Next the residual APSE in is formed 14 R3311 N 1511 Matt 529 by assembling elementlevel quantities At this stage a convergence Check is performed If Bill is sufficiently small quit the iteration ie go to the next time step otherwise continue The linear system MAa 2 R531 523 is solved for Au the matrix M is de ned later A corrector phase is performed at a lwzsa 21 2 vim1 5 vnl7Atai 39 8 CHAPTER 5 FINITE DEFORMATION ELASTODYNAMICS uglmsma 5 24 13 139 gamma 2 625 ammo W Ad The counter is incremented 239 i l 525 and calculations continue by forming a new residual To solve the linear system 523 the effective mass M must be de ned The nonlinear equation is Man1 Nldn1l I F5511 526 and so in each iteration we want to satify MaSINd f11 lt52 Where 3 new 528 11551 dlfll m2na I CLASSICAL TIMESTEPPING ALGORITHMS 9 in which Act is unknown Linearizmg about the state t 2 N 1311 N 1511 M 3At2Aa OH ama12 529 27 Where K is the consistent tangent Neglecting higher order terms and substituting for in ternal forces yields the linear system M AtZK Aa F3531 n N 4191 Mag 530 V M 1211 which de nes the effective mass matrix the consistent tangent for the algorithm The problem Mum 12331 531 may be viewed in the same way as KAd 2 12331 532 of statics in Which the residual contains no intertial term The consistent tangent M could be replaced by F including BFGS updates etc Generalizations 1 Other transient rithms could be considered such as the amethod 1 p 532 The predictor and corrector formulas of Newmark s method 521 and 524 are retained but treatment of outo balance forces in the nonlinear timediscrete equation of motion 3910 dropped g n1a CHAPTER 5 FINITE DEFORMATION ELASTODYNA HCS is modi ed Man1 1 allFi i quot Nldn1ll alFZXt 39 The Nevvmark method is the special case of a 0 To solve this system of nonlinear algebraic equations we linearize as before obtaining M i a At2K Aa Eff 534 h F d M3 Where R231 1 was Nltdn1gtgt M2 Ndn Matti 535 Explicit predictorcorrector algorithms are obtained by replacing correctors With pre dictors in the timediscrete equation of motion 1 p 553 The variant of Newmarkis method Mam N 2 F3511 536 is solved With a diagonal mass matrix as the effective mass in one pass since R2031 F3521 N avail quot Man 0 537 The generalization to explicit predictormulticorrector algorithms is obtained by em ploying displacements from the previous iteration in the time discrete equation of mo CLASSICAL TIMESTEPPING ALGORFI HMS 11 tion For example in the Newmark algorithm Matti N 811 538 A second pass can increase time accuracy in the presense of damping 1 p 563 The update is Obtained by solving the standard linear system Where again MM 4 In the presense of damping eg nonlinear Viscoelasticity the timec0ntinuous equa tion of motion is mg Nd d I ext 539 Discretizing in time by the predictormulticorrector algerithm based on the Newmark method leads to n1 Mags N dffjff v81 Iv t 540 Linearizing as usual about the state I 8N dg Squot N am as N as v521 WA T Ad NV 1011 vl ll quot 7 n t 2 2 8 pm 0 AdH lAvH 541 quot39 V quotquot39quot39 A D C Where C is the consistent tangent damping matrix The standard linear system 523 4 12 CHAPTER 5 FINITE BEFORMATION ELASTODYNAMICS is solved for Aa where in this case the effective mass matrix is Mx M VAZC39 3At2K 542 and the residual is R311 F3531 N dgllavSil Magi1 5 Employing implicitexplicit nite element mesh partitions lg pp 552 564 requires as in the linear case nothing more than substituting into the effective mass matrix the assembled damping and stiffness matrices of the implicit group for the full damping and stiffness matrices respectively namely M M 7301 6At2KI 544 where 8N1 I C 81 545 8N1 I K 39 8d and N I is taken only over the group of implicit elements Nothing else in the formu lation including the right hand side changes 54 STABILITY 13 Figure 51 top of 1 Fig 941 p 555 Generalizations of this kind to the amethod including implicit explicit mesh parti tions are derived in 54 Stability Does the concept of numerical stability work There are various notions of stability in the nonlinear case in linear analysis numerical stability implies that a conservation law or decay inequality imposes limits on the growth of solutions This notion is combined with consistency based on the local truncation error and related to order of accuracy to prove convergence 1 pp462 468 A weaker notion of stability is actually suf cient to prove convergence In the nonlinear case a notion of stability should be at least as strong as the one used in linear analysis Thus analogous results to those in the linear case may be obtained This leads to the notion of linearized stability which is Widely used in nonlinear analysis However7 this concept does not necessarily impose limits on growth of solutions to nonlinear problems It provides only necessary conditions for stability and therefore it is not suf cient on its own7 to guarantee success in practice Consider for example the Newmark method with predictors an and 5M1 given in 518 correctors dn1 and on given in 519 and the nonlinear time discrete equation of V 14 CHAPTER 5 FINITE DEFORMATION ELASTODYNAJVHCS iiiotion Man Ngdn1 vn1l Fifi To analyze the lineariZed stability of the time stepping scheme which is a separate issue from the iterative strategy used to solve the nonlinear algebraic equations the reaction of the algorithm to small perturbations is examined if the algorithm is well behaved in this sense it is called stable this is a very weak notion The algorithm is written in terms of a state vector z as fz 0 547 where T zT lta 17vg13d3 17vf1 dn17ag lv 7d gt 53948 and an an Am egg 25m 611 Un w fz an an mam l 549 vn1 n1 7At0n1 Man1 Ntdn1e quotmil Fiji Consider perturbations of the state vector 6z such as small changes of initial data The algorithm becomes fz 6z 0 550 l 54 SIABILITY 15 Where T T T T MT T 6 quot lt6GR16vn16dn17ovn1 6dn17 Sag vgadg 551 To determine how the perturbations 6z evolve from one time step to the next fz 62 is expanded about z 0 2 fz5z Dfz 6z 0 15412 552 0 This leads to linear equations for 6z Dfz 6z 2 0 553 where Dfz6z 2 z i e 554 80 621W 6d mm egg1 m 2mmquot 617ml 6390n At 39y6an mm m 63n1 5At26an1 6Un1 6gn1 7At6an1 8Ndnla vn1 a1Vdn1 1 n1 M60n1 6d 6dn1 av 517n1 0 K C All the linear terms are replaced by the corresponding perturbed quantities This is the Newmark algorithm for a linear system with mriable stiffness and mass matrices and no 16 CHAPTER 5 FINITE DEFORMATION ELASTODYNAMECS loading In the linear case the algorithm is modally reduced to single degreeoffreedom form with necessary assumptions on the damping and then cast in rstorder form as a one step multivalue method yn1 2 Ayn Where A is the ampli cation matrix 1 p 492 which leads to yn1 2 An yo 556 Spectral analysis is then employed to determine stability in terms of the eigenvalues of A 1 pp 497498 In the present case A varies from one time step to the other so that yn1 An l lAn A1310 557 which is much more complicated For mild nonlinearities it can be assumed that K and C are slowly varying leading to the usual spectral conditions on each An The natural frequency a and the damping ratio 6 in the Newmark stability results are replaced by current values and These are necessary conditions Remarks 1 These results obviously work for the linear case in fact they are identical 2 In the nonlinear case the time step restriction for conditionally stable explicit schemes seems to be suf cient 54 STABILITY N For implicit schemes that are unconditionally stable in linear analysis maintaining time accuracy is suiiicient for stability in the nonlinear case Poor time accuracy may cause biased accumulation of truncation errors leading to pathological results This is examined in more detail in the following in summary linearized stability results are Viewed as necessary conditions but caution is required with implicit schemes in which automatic time stepping strategies are important 5 For general follow on reading see 4 Chap 2 and in particular Sec C which contains convergence proofs Energy Methods Energy concepts may be invoked to assess growth and decay of solutions of nonlinear elastodynamics Consider the timecontinuous case with no damping Premultiplying the 39 T equation of motion by d yields d and d dTFe 558 Where the total energy E Tm Ud 559 is composed of kinetic energy T dTMcl 560 and strain energy U its gradient is Integrating from time 0 to 25 leads to the 39 l8 CHAPTER 5 FINETE DEFORMATION ELASTODYNAMICS fundamental energy identity f 39r Erasmus z Ed0d0O d mech 561 which solutions of the initialvalue problem must satisfy In particular if F 2 0 total energy is conserved Emmett Eidmmon 562 at all time as in the linear case 1 p 457 The total energy plays the role of a norm squared in the displacement velocity phase space Keeping track of this quantity during a calculation gives an indication of the stability of evolution Example Consider small def0rmations in a materially nonlinear elastic 0nedimensi0nal rod model 100d Id 3 2 Nd 5 63 200 sgnd 1 gt 2 3 The material law is representing a bilinear elastic since there is no hysteresis spring relation whose tangent is zero for d gt 2 Fig 52 The strain energy is the area under the Nd curve Consider the trapezoidal rule often refered to as the average acceleration method in structural dynamics which is Newmark s method with 8 14 and 7 12 It has been shown that in the absence of dissipative effects and external forces total energy is asymptotically conserved for small and large time steps and that attainable values of En z Edmvn for the trapezoidal rule are bounded as qualitatively shown in Fig 53 39 54 STABILITY 19 Figure 52 Bilinear spring with zero tangent material law Figure 53 Fig 2 Consider the initial value problem ciNd 0 d0 0 564 am 2 25 This problem With the material law 563 exhibits interesting numerical phenomena Taking At 02 is neither a small time step nor a large one for both of which total energy is conserved but in the intermediate range In the linear branch 02 k m 100 The period is T 27rw 2 0628 and thus AtT z 03 so that the selected time step is not suf ciently small for time accuracy The solution of the trapeziodal rule in exact arithmetic with a time step of At 02 exactly reproduces itself every 24 time steps with Iain S 5 12 lt 30 and 3125 E0 S En 3 9125 Numerical calculations initially manifest this behavior accurately but for longer times the numerical solution deteriorates and a persistent pattern of total energy growth appears This drift from the expected behavior is due to truncation i 28 CHAPTER 5 FINITE DEFORMATION ELASTODYNAMICS Figure 54 bilinear hardening law errors It is observed that a certain amount of damping stabalizes the calculation but a thorough analytical understanding of these results is still lacking The important conclusion from this example is that the unconditional stability of the trapezoidal rule in linear analysis is not automatic in the nonlinear regime and thus caution should be exercised in its use Certain types of nonlinearity may have to be approached with greater care than others and success for one class of nonlinear problems is no guarantee of success for another These statements hold for many methods of temporal integration in nonlinear analysis Example Also presented in 3 is the initialvalue problem cZNd 0 em 4 5 65 d0 0 in one case with a bilinear hardening spring law Fig 54 as in foamtype materials and in the other bilinear softening Fig 55 For the hardening spring law total energy is conserved for small time steps as expected but at larger time steps a periodic pattern of energy growth and decay appears see Fig 56 5 54 STABILITY 21 Figure 55 bilinear softening law Figure 56 Fig 4 the period in this case is approximately As the time step is increased further the response begins to ip op and En approaches E0 Recall the behavior of the trapezoidal rule in linear analysis In the absence of damping and external loading the total energy is conserved Furthermore as At gt 00 on gt 1 v0 566 We see that this holds in the nonlinear case as well Energy conservation for large time steps and a pattern of periodic growth and decay for an intermediate time step are observed in the softening spring law as well Fig 57 In this case however the numerically computed energy is bounded from above by the initial value Whereas the opposite is true for the hardening spring law Remedies HOW may these pathologies be avoided 0 CHAPTER 5 FINITE DEFORMATION ELASTODYNAMICS Figure 57 Fig 6 One may employ better algorithms than the trapezoidal rule such as the amethod which incoporates numerical damping Alternatively the algorithm may be modi ed so that a discrete counterpart of the energy identity is satis ed as presented in the following For implicit schemes that are unconditionally stable in linear analysis time steps that are suf ciently small to provide time accuracy must employed thus assuring stability in the nonlinear regime Adaptive strategies based on a priori estimates of the solution error as a function of the time step would be preferred but unfortunately are usually valid only for linear problems In nonlinear analysis a posteriori estimation of the error is employed Bil8 In this approach a time step is selected the resulting solution computed and then employed to estimate the error made and nally if need be the time step is modi ed accordingly Energy conserving Algorithms Stronger notions than linearized stability require different classes of algorithms than those Which are derived directly from linear analysis One such algorithm may be obtained by a modi cation to the trapezoidal rule for nonlinear applications resulting in physically appropriate energy growth characteristics In particular when external forces are absent energy is identically conserved in nonlinear elastodynamics a property shared by the exact 54 STABILITY 23 continuum equations As a consequence unconditional stability is automatically attained The following exposition applies to linear elastic systems for which N 2 The trapzoidal rule Newmark s method With 3 14 and 7 1 may be written as Mam Nidn1l F311 At dn H dn 17n1 557 At vn1 vn l Tquotan quotl NIH 1 2 Where do and 3900 are given and a0 1 1F0 Ndg In practice vn and an1 are eliminated from the equations resulting in a nonlinear algebraic problem for Lin 4 4 EMdml Mean 2 F331 M an mm mm 568 We Wish to impose energy conservation as a constraint condition on the solution For this purpose the nonlinear algebraic problem is viewed as the EulenLagrange equation of a functional The Constraint is then enforced by employing the Lagrange multiplier method see eg 1 p 195 The functional is 4 2 an d 1Mdn1Udn1dg1 F3311 M an EM Am 569 Then 8fdn1 T 4 ex 4 saf m 6dn1 Ej Mazn1 Ndn1 m1 M an Em Aim 570 24 CHAPTER FINITE DEFORMATION ELASTODYNAMICS Thus the nonlinear aigebraic problem is equivalent to 50 being identically Zero for arbi trary gvariatiousquot 6dn1 The trapezoidal rule is modi ed so that the following discrete energy identity is satis ed 1 a En1 2 En 3mm dnTF E1 F2 571 This is the algorithmic equivalent of the fundamental energy identity 561 When F 0 this implies that the total energy of the approximate solution is identically conserved To achieve the desired modi cation let 1 gldnu En1 En dn1 dnTF f1 F2 572 Where En1 E dn1 dn1 dn 22 and the second argument is the expression substituted for vn1 from 567 Thus the equality GRILL H 0 573 implies satisfaction of the discrete energy identity The modi ed algorithm is based on 567 with the nonlinear timediscrete equation of motion replaced by the EulerLagrange equation of the augmented functional fdn1 Agdn1 Where A is a Lagrange multiplier Differentiating the augmented functional yields adnu adnu l WWW i 55 CONSISTENCY ACCURACY For this expression to vanish for arbitrary 6dn1 and 6A requires satisfaction of 573 and 4 i 1 1 M15de Mains t1 3A F3251 3mg 575 4 1 Mn 1dn Atlt1 gt H aAt2 2Avn 0 This is the modi ed equation of motion which accounts for the forces induced by the Lagrange multiplier in order to enforce the discrete energy identity The solution of equations 573 and quot75 for unknowns dn1 and A by an iterative algorithm is not substantially more expensive than the unmodi ed trapezoidal rule An effective procedure is discussed in 55 Consistency Accuracy CHAPTER FINITE BEFORMATION ELASTOD YNAMICS Bibliography 1 TJR Hughes The Finite Element Method Linear Static and Dynamic Finite Element Analysis PrenticeHall Englewood Cliffs New Jersey 1987 2 l Miranda TJR Hughes and RM Ferencz Earthquake Engrg Struct Dynamics 1989 3 TJR Hughes Stability convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics Comput amp Structures 6 1976 313324 4 T Belytschko and TJR Hughes Computational Methods for Transient Analysis North Holland Amsterdam 1983 5 HD Hibbit and BI Karlsson Analysis of pipe whip Pressure Vessels and Piping Conference ASME 79PVP 122 San Francisco California 1979 6 C Johnson Numerical Solution of Partial Differential Equations by the Finite Element l vlethods Cambridge University Press Cambridge 1987 M7 2
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'