INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS FOR NON-LINEAR FINITE ELEMENT ANALYSIS OF SHELLS P. NORDLUND,* A. E. GIANNAKOPOULOS AND B. HA®GGBLAD Department of Solid Mechanics, Royal Institute of ¹echnology (K¹H), Stockholm, Sweden Department of Materials Science and Engineering, MI¹, Cambridge, MA 02139, º.S.A. ABB Corporate Research, »ästeras s, Sweden ABSTRACT We consider methods for adaptive updating of computational meshes during incremental loading of non-linear shell and solid structures. In particular, we focus on updating methods where the initial topology of the mesh is maintained. The movement of the mesh (the convective step) is decoupled from the Lagrangian solution step and is done by using separately computed displacement fields, obtained by solving auxiliary thermoelastic problems. The properties of this mesh-updating procedure are investigated and demonstrated on several examples. 1998 John Wiley & Sons, Ltd. KEY WORDS: adaptivity; node relocation; ALE; non-linear finite element analysis 1. INTRODUCTION A major challenge in engineering computations is the analysis of non-linear shells. Many structures in automotive, naval, aerospace and civil engineering contain components in the form of plates and shells, that may undergo large deformations during their service lives. Even with elastic material response, complications in shape, loading and contact conditions might put great demands on the computing power that is needed for analysing different design alternatives. For large deformations the computational mesh will eventually become severely distorted and a complete remeshing might be needed to complete the analysis. However, by moving the grid during the incremental loading process in such a way that excessive mesh distortions are avoided and refined parts of the mesh are forced to follow the regions with high deformation gradients or large changes in the local geometry, the accuracy and reliability of the model can be maintained for a larger range of loading and the need for complete remeshing is postponed for as long as possible. Updating the mesh according to this arbitrary Lagrangian Eulerian (ALE) procedure is normally computationally much more efficient than performing complete remeshing steps repeatedly. *Correspondence to: P. Nordlund, Department of Solid Mechanics, Kungl Tekniska Högskolan, S-100 44 Stockholm, Sweden. E-mail: pern@hallf.kth.se Contract/grant sponsor: ABB Contract/grant sponsor: Swedish Research Council for Engineering Sciences CCC 0029—5981/98/081523—22$17.50 1998 John Wiley & Sons, Ltd. Received 6 March 1997 Revised 7 April 1998 1524 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD The ALE approach was originally used in fluid mechanics simulations as a complement to the pure Eulerian approach in order to handle free surfaces and interacting surfaces between fluids and structures effectively, see for example the early works of Hughes et al. and Donea et al. In solid mechanics there is the additional need for tracing the material particles (the Lagrangian view) in order to capture history or path-dependent material properties, and the development of free surfaces. This means that the introduction of ALE into the Lagrangian analyses might lead to a substantial increase of computations per analysis step, since remapping of many mechanical quantities (mass, momentum, stresses, strains, etc.) must be done between different meshes. Hence, much research has been devoted to developing methods for effective mapping (projection, advection, global and local smoothing of involved fields), strategies for less frequent updating, and effective procedures for mesh-updating. Clearly, the use of ALE methods in solid mechanics has made it possible to solve many complex non-linear problems, such as impact and metalforming processes (wire drawing, extrusion, rolling and sheet metal-forming). For some pioneering works, see References 3—6. A non-trivial task in the ALE approach is the construction of an updated mesh that has desirable properties, such as no inverted elements, appropriate directional grading, good overall approximation properties, etc. Currently, a satisfactory and fully automated adaptive procedure is still lacking. Updating by relocation of nodes (the r-method) has been studied throughout the development of the finite element method, often as a stand alone remeshing procedure outside the scope of traditional ALE, where the size of the movement of the nodes is usually kept small and the movement itself might be performed directly coupled to the solution step. Many techniques have been proposed in the context of node relocation algorithms. Early examples are various types of direct optimization of the potential energy with respect to both nodal displacements and nodal co-ordinates. In strongly non-linear problems (e.g. in metalforming) it is important to take away accumulated mesh distortions, which can be done, for example, by averaging or using smoothing methods based on Laplace or Poisson equations.\ However, if smoothing is enforced by depressing the mesh gradients towards zero, these methods tend to eliminate any existing mesh grading. The perhaps most successful r-adaptivity formulation uses an error-weighted nodal averaging formulation. This method is simple and fast, but it allows no direct control over element distortion. Furthermore, special techniques are necessary to handle curved boundaries and it works well only when a structured mesh is used. Torigaki used a similar approach, which suffers from essentially the same limitations as Diaz’ method. In the fluid mechanics community, a common formulation is based on a variational principle similar to the one of Brackbill and Saltzman together with optimization techniques for the relocation of nodes. Important quantities in adaptive mesh-updating techniques are error indicators. These can be constructed using geometric mesh distortions or a posteriori estimations of local errors based on the differential equation residuals or on the difference between the finite element solution and a post-processed smoothed solution. Procedures for adaptive updating of shell meshes are not so common in the literature. Works by Belytschko et al., Stein et al., Bonet and Okstad and Mathisen focus mainly on enrichment procedures. In the mesh updating there are additional difficulties for shells, because of the strong influence of the local geometry on the behaviour. Special precautions must be taken to preserve the shell geometry during the movement of the mesh. In this work we introduce a decoupled mesh-updating technique that permits ‘aggressive’ node relocation for obtaining meshes with strong directional grading. We strain the current mesh by 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS 1525 prescribing a fictitious thermal strain field in such a way that it will reshape the mesh towards a desired form. By assigning appropriate thermoelastic material properties to the body described by the mesh and solving the thermoelastic boundary value problem with appropriate boundary conditions we obtain a smooth approximation of the desired global mesh. The smoothing and integrity (no overlaps) of the mesh are secured by the elasticity of the fictitious material that the mesh is thought to be made up of. Since we are free to choose the material parameters and their spatial variations this mesh-updating procedure can easily be adjusted to special demands. The geometry of the structure must be retained during this node relocation analysis. This is handled by constructing frictionless contact surfaces surrounding the structure. Since the contact problem is of a simple nature the additional computational costs are in most cases negligible. The application of this mesh-updating procedure to solids is rather straightforward. In this paper, we concentrate on the harder case with geometrically non-linear shells. Similar approaches to mesh-updating can be found in Reference 24 for solids and in References 25 and 26 for fluids. However, in our approach we use a scalar field (temperature) for monitoring the mesh-updating which makes the method very flexible in complex applications, involving for example shells. We outline the basic idea behind our procedure in Section 2, and Section 3 is devoted to the theory of the updating with thermoelastic displacement fields. In Section 4 we describe some of the strategies for controlling mesh density and mesh quality. Section 5 discusses the general setting of the decoupled ALE approach and in Section 6 we perform several numerical experiments. 2. THE PROPOSED MESH-UPDATING ALGORITHM In this section we outline the basic idea behind our mesh-updating algorithm. We want to update the finite element mesh of a given non-linear problem by displacing the nodes to new positions so that an improved mesh (in some sense) will result. For that purpose we think of the mesh as describing a separate elastic body which we can deform by applying for instance forces and prescribed deformation. A prescribed strain field is applied to the current mesh such that the mesh will reshape towards a desired (target) form both with respect to mesh density and element shapes. The prescribed strain is imposed elementwise as a thermal strain tensor e "a ¹, where a is an (anisotropic) GH GH GH tensor describing the thermal expansion and ¹ the increase in temperature. The tensor a governs GH the relation between volume and deviatoric strain with the scalar field ¹ as a scaling factor. The thermal strain field is constructed in such a way that nodes are moved into regions where a higher mesh density is needed, and also so that the shape of individual elements is improved. The intuitive physical picture is of course that ‘warm’ regions will swell and others will contract, producing a node relocation which will redistribute the mesh according to, say, an error estimation. A simple example is shown in Figure 1 where the boundaries are straight and the treatment of the boundary conditions is trivial. In a general case with curved boundaries, the nodes should slide along the curved boundaries. This is accomplished by defining rigid contact surfaces (or lines in two-dimensional problems) surrounding the structure and permitting frictionless sliding without separation. By imposing the volume strain e "2a¹"2xy!1/2 on the mesh in Figure 1(a) we obtain the mesh in Figure 1(b). In this example, we could equally well have used the volume strain e "2a¹"2xy, but by subtracting the mean value of the thermal strain over the mesh (1/2) we can easily identify the curve xy"1/4, around which the elements will retain their volumes (see the 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1526 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD Figure 1. Nodal relocation algorithm applied to a square 2D mesh dashed lines in Figure 1). This kind of partition of the mesh into expanding and contractive areas will always be present in our mesh-updating method. Additional control of the element shapes is obtained by using orthotropic thermoelastic properties. The proposed mesh-updating method displays several good features, for example: 1. 2. 3. 4. it is generally applicable to solid and shell problems, mesh integrity is built in (no overlapping of elements is possible), the basic algorithms needed are available in most non-linear finite element codes, a good intuitive physical interpretation of the mechanism behind the movement of the nodes is available. 3. THEORY In this section we outline the basic thermoelastic relations for solids and shells and reinterpret the energy functionals as smoothing and mesh redistributing functionals. For shells we also discuss the boundary and curvature constraints that must be supplied. 3.1. Solids For a linear, thermoelastic solid the stresses are obtained as E l p" e !a ¹# (e !a ¹) V 1#l V V 1!2l E l p" e !a ¹# (e !a ¹) W 1#l W W 1!2l (1) E l p" e !a ¹# (e !a ¹) X 1#l X X 1!2l q "G(c !a ¹), q "G(c !a ¹), q "G(c !a ¹) VW VW VW VX VX VX WX WX WX 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS 1527 where E and l are Young’s modulus and Poisson’s ratio, respectively, a , a and a the thermal V W X expansion coefficients in the extensional directions and a , a and a the apparent coefficients of VW VX WX thermal shear present in reference systems that are different from the principal axis system of the tensor a. ¹ is the temperature and e "e #e #e the volume strain. For brevity, we have also V W X introduced the volumetric thermal expansion coefficient a "a #a #a V W X The increment in strain energy is given by (2) d¼" (p de #p de #p de #q dc #q dc #q dc ) d» V V W W X X VW VW VX VX WX WX (3) Inserting relation (1) into equation (3) and integrating gives E ¼" 2(1#l) l (e !a ¹)#(e !a ¹)#(e !a ¹)# (e !a ¹) V V W W X X 1!2l #(c !a ¹)#(c !a ¹)#(c !a ¹) d» VW VW VX VX WX WX (4) or, choosing the local co-ordinate axes along the principal directions of the a-tensor, E ¼" 2(1#l) (e !a ¹)#(e !a ¹)#(e !a ¹) V V W W X X l # (e !a ¹)#c #c #c d» VW VX WX 1!2l (5) If there are no other external forces, the resulting mesh deformation (movement) is determined by the requirement that the strain energy ¼ should be minimum over the set of functions that satisfies appropriate geometric boundary conditions (zero normal displacements and free tangential displacements). a ¹ are the imposed thermal strains, and the minimization of the functional GH will give the solution which has strains closest, in the energy sense, to the thermal ones. Note that for c "0 and l"0 we obtain a pure least-squares approximation. GH Since l (6) a ¹"a ¹#1!1" V !1 V V lI V where l and lI are deformed and undeformed side lengths, respectively, we can use a thermoelasV V tic solution for adjusting the shape of the mesh elements towards a prescribed mesh (given in terms of thermal strains) in a least-squares sense. Directional grading might be obtained by setting, for example, a O0, a "a "0 and l"0. Clearly, by using different values of the V W X components of the a-tensor it is possible to control the mesh distortion. The global behaviour of a ¹ will govern the mesh density distribution (the element sizes) and the local orthotropy of the a-tensor the shape changes. 3.2. Membrane shells Consider a membrane shell with the differential arc-length ds"A(dm )#A (dm ) 1998 John Wiley & Sons, Ltd. (7) Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1528 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD where m and m are the curvilinear co-ordinates of the surface. A dm and A dm are the differential arc-lengths along the local co-ordinate axes. w is the displacement normal to the shell surface. The strains are given by the expressions u *A w 1 *u # # e " A *m A A *m R 1 *u u *A w # # e " (8) A *m A A *m R A * u A * u # c " A *m A A *m A where u and u are the displacements in the shell surface and R and R are the radii of curvature. Supposing thermoelastic orthotropy, Hooke’s law can be written E p " [e !a ¹#l(e !a ¹)] 1!l E p " [e !a ¹#l(e !a ¹)] 1!l (9) q "G(c !a ¹) where E, l and G are the conventional moduli and a , a and a the thermal expansion coefficients. The thermoelastic functional for shells that corresponds to relation (5) for solids can be written E ¼" 2(1#l) l (e !a ¹)#(e !a ¹)# (e !a ¹)#(c !a ¹) A A dm dm 1!l (10) with e "e #e and a "a #a . The desired thermoelastic displacement field will minimize the strain energy of the membrane shell within the set of fields that satisfy the geometric boundary conditions including the out-of-plane constraint w"0. Remark 1. By introducing the new co-ordinates xN "x#u yN "y#v (11) zN "z#w we obtain *xN *yN e " !1, e " !1 etc. V *x W *y (12) Introducing this in (5), the energy ¼ can be expressed in terms of the ‘smoothing’ integrals that were used by Brackbill and Saltzman in their approach. In an approximate way our thermoelastic approach can then be reinterpreted in the terminology of Brackbill and Saltzman. 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS 1529 Remark 2. An obvious extension of the present approach is to use a hyperelastic material model instead of Hooke’s law. However, if the incremental loading steps are kept small the corresponding movement of the mesh will also be small and the linear elastic approach will be accurate enough. 3.3 Boundary conditions and curvature constraints Of great importance for the success of the method is the proper treatment of the boundaries in the auxiliary thermal stress analysis. The current positions of the boundary nodes define the boundary to be used in the non-linear thermoelastic analysis which is performed with updating of geometry (material or Eulerian formulations). This boundary must be preserved during the node relocation. The simplest choice in two dimensions is to let the boundary be a rigid slide line, defined by straight lines between the boundary nodes. This works well for boundaries that remain overall straight, but has the disadvantage that the nodes are hindered in their motion along a curved boundary due to the step-wise change in the line normals. In the method applied here we create interpolating Non-Uniform Rational B-Spline (NURBS) lines using the nodes between corner points. The NURBS lines are then directly used in the analysis as rigid slide lines. This will have a smoothing effect on the boundary which in some cases may be unwanted. In most cases, however, the smoothing gives a more accurate description of the boundary. In analyses that involve contact there might be problems. A node relocation step may create gaps or overlaps between the contacting bodies of such magnitude that the subsequent analysis will be influenced. A possible remedy to this problem would be to project the new mesh onto the old boundary, but this does not guarantee the elimination of overlaps. Another option is to enforce proper contact conditions afterwards by an equilibrium correction. This will induce spurious stresses, but if the gaps and overlaps are reasonably small these stresses will be small and lose significance in subsequent steps. This is the approach used in this work. An even better solution would be to smooth the boundary locally, based on neighbouring segments, and use this smooth boundary for both the contact surface in the main analysis and in the rigid sliding surface in the node relocation analysis. The problem with spurious gaps and overlaps is then completely avoided. In a shell structure, all nodes must be constrained to stay on the surface defined by the positions of the shell elements before relocation of the nodes. Analogous to the approach in two dimensions, one or more NURBS surfaces are interpolated from the nodal positions. The NURBS surfaces are interpolated from local patches of neighbouring element segments. This is done in order to allow general unstructured meshes in the model. In order to simplify the implementation, only larger regular patches are used in the examples below. The free edges of the structure (if any) need special treatment. Sliding surfaces perpendicular to the elements at the free edges are constructed, constraining the edge nodes to stay on two perpendicular sliding surfaces. Internal edges (fold lines), appearing, for instance, if a thin-walled box is modelled, are treated automatically with the current approach. Nodes lying on an edge are constrained to stay on the intersection between two rigid surfaces. The definition of edges and the distinction between free and internal edges is done at model definition time. If there is an angle between elements at a boundary between two surface patches due to discretization, it should be allowed to be smoothed. On the other hand, if it is a sharp edge, it should be retained by letting the nodes only slide along the edge. The mesh is supposed to slide along the contact surfaces and thus no bending stiffness of the elements is needed or even wanted. Therefore, membrane elements instead of full 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1530 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD shell elements are used for the node relocation analysis, which reduces the cost of this analysis substantially. The extension to three-dimensional structures is straightforward: corners are fixed and contact surfaces are defined in terms of NURBS surfaces interpolated from boundary nodes. In fact, the treatment of three-dimensional structures is usually simpler than the handling of shell elements with their special kind of free edges. 4. CONTROL OF NODE RELOCATION 4.1. Thermal strain fields The previous sections showed how to reshape the mesh by solving a thermoelastic problem. The control of the mesh density distribution is obtained with an appropriate choice of the thermal strain field. In each element we set up a thermal strain tensor in such a way that the resulting nodal displacement will change the element towards a state with desired volume and shape. In our previous analysis we also concluded that the global distribution of the thermal volume strain a ¹ (see equation (6)) tends to impose elemental volume changes and the relation among the components of the a-tensor tends to govern the changes in elemental shapes. In this work we will consider volume strain dependencies of the type a ¹ "f (w , » ) (13) G G G where f (w , » ) is a monotonic function of the product w » . Here w is an elemental weight G G G G G function (local error estimator, stress, strain, curvature, etc.) and » is the current volume of G element i. A typical form of the function f can be derived in the following way. In the most common cases, we want to adjust the current element sizes so that the sum of the products of the weight functions w (the local error) and the corresponding adjusted element sizes (» #*» ) will be as small as G G G possible. This will obviously minimize the global error. The optimal solution is characterized by (w #*w ) (» #*» )"const. (14) G G G G (no sum over i). As a first approximation we use w (» #*» )"const."c G G G (15) where c can be defined as c" » (16) 1/w I with » denoting the total mesh volume. Thus, *» 1 » 1# G" » w » 1/w G G G I which means that a thermal strain field that varies as c a ¹" !1 G w» G G 1998 John Wiley & Sons, Ltd. (17) (18) Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS 1531 Figure 2. The volume strain functions will produce displacements that tend to adjust the current mesh towards the desired (optimal) one. Hence, in this case we obtain c f (w , » )"f (w , » )" !1 (19) G G G G w» G G The behaviour of this function is shown in Figure 2. The parameter c plays the role of a ‘target’ quantity, the elements with f '0 will expand and for f (0 they will shrink. It can be computed according to equation (16) (a kind of mean value) or be supplied by the user (see next subsection). The alternative w» f (w , » )"! G G #1 G G c (20) is a local approximation to f in the vicinity of the target point (c, 0) . In practice, the usable interval for a ¹ is (!1, 1). The functional form f has the disadvantage that it gives infinite strain values for w "0, a situation that is not uncommon. Therefore, the form f is the one that is G utilized in this work. However, this form of f is only useful near the target point. Used globally it might lead to annihilation of volumes ( f "!1) . We could use a normalized version of f : w» 1! G G c f " (21) 1#"w » /c" G G that is well-behaved for all w » /c. If several weights are to be considered we set up a resultant G G weight wN according to G w w wN "k #k G #k G #) ) ) (22) G w w with k #k #k #) ) )"1. The values w, w, . . . are the chosen comparison weights (e.g. mean values). The first term on the right-hand side will promote equal element sizes (weight function w "1) and the other terms consider additional element quantities (error estimations, G 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1532 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD Figure 3. The mapping of an isoparametric element stresses, strains, curvatures etc.) that should govern the equidistribution. The factors k can be G used to emphasize different parts in equation (22). In addition to monitoring the volume strain for controlling the element sizes we need a procedure for controlling the element shapes. Consider the isoparametric mapping from the parent rs-domain to the actual xy-domain for a four-noded quadrilateral plane element as depicted in Figure 3. It is obvious that the local properties (local extensions and shears) can be obtained from the Jacobian matrix that relates corresponding ‘fibres’ in the two domains. The determinant of the Jacobian gives the ratio between the mapped and the original differential volumes at the chosen point in the element. However, for the special case of a four-noded quadrilateral element a simple geometric construction can be used to obtain relevant information about the element distortion. Consider Figure 4, where the parallelogram ABCD has been constructed by drawing the side line AB in parallel with R OR , BC in parallel with S OS , etc., through the midpoints of the sides in the given quadrilateral ABCD. It is easily seen that the two quadrilateral elements ABCD and ABCD have equal total areas. The lines R OR and S OS between the midpoints of the sides divide the parallelogram into four equal subareas. The corresponding subareas in the quadrilateral element ABCD might be larger (like OS DR ) or smaller (like OS BR ) than the fourth of the total area. The reshaping of the quadrilateral ABCD can now be thought of as being done in two steps. In the first step, the element ABCD is reshaped towards the parallelogram ABCD by imposing appropriate thermal volume strains at the integration points that are located in the different subareas. For example, at the integration point Q (see Figure 4) the thermal strain is set equal to a superposition of a local ‘reshaping’ part (*e ) and a global ‘resizing’ part (e ): e "*e #e (23) where area (OR AS ) !1 *e " (24) area (OR AS ) and e is given by e.g. equation (20). In the second step, the shearing of the parallelogram ABCD (and the element ABCD) is counteracted by imposing the volume strains at the integration points with the aid of orthotropic thermal strains in appropriate principal directions. Using the notations in Figure 4 we obtain the following components for the vectors P and P pointing in the principal directions: P "(1, 1!k /2) (25) 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS 1533 Figure 4. Deformations of a general parallelogram Figure 5. An application of element shape control: (a) original mesh; (b) element shape control and P orthogonal to P (see Figure 4), where k is the amount of shear, k "tan . in Figure 4 is supposed to be much less than n/4. If the imposed thermal strains in these directions are set to e e " !k /2 2 (26) e e " #k /2 2 (27) and respectively, the shearing is counteracted under the constraint of keeping the volume strain at the level of e . An example of this shape control strategy is given in Figure 5. In Table I we have summarized the algorithm for controlling the element shapes that is used in the node relocation procedure presented in this work. Remark 3. Obviously, reduced integration of the four-noded quadrilateral element does not allow this type of element shape control and should therefore be avoided in the auxiliary thermal analysis. The extension to other element types such as higher-order elements as well as threedimensional elements is in principle straightforward, but some minor modifications of the element shape control procedure might be needed. 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1534 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD Table I. Algorithm for controlling element shapes For each element j compute wN "k #k w/w (resultant weight) H H eH "(!wN » /c)#1 (mean volume strain) H H k "tan (amount of shear, ;n/4) H H H P "(1, 1!k /2), P NP (principal directions) H In each of the four integration points (i), impose in the directions P and P eHG"eHG/2!k /2 and H eHG"eHG/2#k /2 where H eHG"e #*eHG with *eG"(4» /» )!1 (» total element area and » subarea i) HG H H HG As an alternative approach to element shape control, one could consider to apply nodal forces directly. This was used by Schreurs et al. and can be applied in the current node relocation formulation as well. 4.2. Some examples of volume strain fields An often used error indicator is the one due to Zienkiewicz and Zhu. It measures the error in stress by evaluating the difference between a smoothed stress, p*, and an unsmoothed stress, p. An error measure e88 is formulated as w » "e88"#p!p*# (28) G G G G with some suitable norm. This error is, for instance, evaluated for each stress component, for an effective stress or expressed as an energy norm. In order to obtain an equidistribution of this error indicator, a volume strain field based on the element based error e88 can be formulated as e88 a ¹ "! G #1 (29) G e88 valid around e88, where e88 is a proper comparison quantity, for instance a mean value of the error indicator. In many non-linear situations, however, ‘real’ error estimates like the one above are not so useful. For example, in elastoplastic problems, equal distribution of some other solution quantity may be more appropriate. Such quantities are the plastic strain energy or the accumulated plastic strain evaluated at each element integration point. The plastic strain energy was successfully used for adaptive mesh refinement (essentially in an h-method) by Cuitin o and Ortiz. An example of a volume strain field based on effective stress is given in Section 6.1 below. Using few elements on a curved surface gives a less accurate description of it, especially if low-order elements are used. We may require equal distribution of elements with respect to the 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS 1535 curvature or compute geometry description error estimates based on, for instance, the angle between neighbouring elements. The algorithm for evaluating the angle, h, between neighbouring shell elements is as follows. First, an element-centred normal is calculated for each element by taking the cross product of the diagonals. These normals are averaged at the nodes, and for each element an angle error is calculated as the maximum angle between its centred normal and the averaged normals at the element nodes. The averaged normals at symmetry boundaries are projected onto the symmetry plane. A volume strain field can be set up as h» a ¹ "! G G #1 (0)h)n/2) G (h») (30) 5. GENERAL STRATEGIES The principal usage of the proposed node relocation algorithm is thought to be in a decoupled ALE procedure. The steps of such a procedure are typically (a) (b) (c) (d) (e) Lagrangian analysis, estimation of solution errors and mesh distortions, grid displacement calculation, mesh movement and mapping, equilibrium correction computation. During the incremental loading our algorithm will produce a sequence of meshes with corresponding sequence of finite element solutions. The adaptive procedure for the design of a new updated mesh in this sequence is based on the previous mesh with its corresponding solution. A crucial question is how the accuracy is affected by the change of meshes. For good accuracy it is in general necessary that the differences between the response curves belonging to different meshes are small, in particular when the response is path dependent. As a rule, response of a mesh will become too stiff if the mesh is loaded and distorted too much. A change to a new mesh will soften the mesh and give an improved solution accuracy so the loading can continue. An equilibrium correction step is used for re-establishing the equilibrium directly after a change of mesh. Alternatively, a more accurate procedure is used when the response is strongly path dependent. Here, the updated nodal positions are introduced in the geometry at the start of the increment and a complete re-analysis of the current step is performed. In exceptional cases this can be done iteratively until the mesh is accepted. In this work, however, the incremental displacements are assumed to be small enough so that a re-analysis of a load increment is not necessary. 6. NUMERICAL INVESTIGATIONS The node relocation algorithm has been implemented in the general purpose finite element code MarcR as a set of user supplied subroutines. The node relocation calculation is performed as an RMARC Analysis Research Corporation, 260 Sheridan Avenue, Palo Alto, CA 94306, U.S.A. 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1536 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD Figure 6. The elastic spherical indentation problem analysis separate from the main analysis. The standard routines in Marc were used for most features of the method, including the remapping of variables and the interpolation of NURBS surfaces and lines. The remapping (transfer of solution quantities from one mesh to another) used in the Marc program is based on interpolation. When a node relocation analysis is to be performed, the main analysis is halted, the thermal elastic node relocation analysis is performed and the updated mesh is put back into the main analysis which automatically continues, starting with the remapping. The node relocation analysis itself is in general performed in more than one increment (step), and the thermal strain field is updated at the start of each increment. In addition to the general approach described in previous sections we have experimented with several alternative volume strain fields essentially based on ad hoc quantities. 6.1. Solution based control of the volume strain As an illustration of how a solution quantity like the effective stress can be utilized, an elastic spherical indentation problem is studied. A large elastic block is indented with a rigid frictionless sphere, see Figure 6. The problem is modelled as an axisymmetric problem, with fixed displacements on the edges away from the indentor. The analytical Hertz solution reveals that the Mises effective stress has a maximum value at the symmetry line, below the surface at a distance of about half the contact radius. The block is first indented 0)025 m and following this the node relocation analysis is performed. Thus, the main problem is actually not solved with an improved mesh, it is only used as an example of how a stress based mesh relocation can be utilized. A volume strain field of the form proposed in Section 4 with the Mises effective stress, p, as an G ‘error’ indicator is defined as p» a ¹ "1! G G G (p») (31) for each element i. The scaling factor (p») is here chosen to be the maximum value of that product in the mesh. The top part of Figure 7 shows the original mesh used and the bottom part shows results using the above volume strain field. Element shape control was also used which 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS 1537 Figure 7. Original mesh and mesh after node relocation for the indentation problem clearly can be seen in the figures. The elements get clustered in the region of highest stress, and it is interesting to note that the mesh tends to get aligned with the von Mises stress field. The dashed line in the right part is a line of constant von Mises stress (p"100 MPa). 6.2. Geometric based control of the volume strain As a demonstration of how the node relocation algorithm distributes a geometric error distribution, the structure in Figure 8(a) is studied. The angle between normals of adjacent elements is used as an error measure for the discretization, as described in a previous section. Applying the node relocation algorithm to this model with originally uniform element distribution and the volume strain formulation (30) results in a node distribution similar to that shown in Figure 8(b). The calculated error (h » ) is now equally distributed over the surface. Note that the G G angle error itself is not equidistributed but the product h» is. In order to obtain equal distribution of h one can use a volume strain field defined simply as h a ¹ "1! G (32) G h This formulation was used in Figure 8(b) where the angle error is constant in the whole structure. The first case, which distributed h», gives a less focused mesh. The scaling factors h and » influence the convergence of the nodal motion, i.e. how fast (in number of increments) h» or h become equally distributed. In general it does not influence the 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1538 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD Figure 8. A parabolic cylindrical shell: (a) original mesh; (b) angle error control Figure 9. A structure with flat regions: (a) original mesh; (b) mesh after node relocation with the combined weight final result. They can be chosen as an averaged value as in equation (16) above or simply as maximum values. The previous examples demonstrate how the method works for cases where the error distribution is smooth and the equidistribution scheme gives a good mesh. One situation where the results may not be satisfactory is when the chosen error estimate becomes zero in a region. This often occurs, for instance, in elastoplastic analyses if an error estimate based on error in plastic strain is used. In elastic regions the error is zero. Another example is shell structures with the angle measure above when the shell is flat in some regions, a common situation in sheet metal-forming. Figure 9 shows such a case with a uniformly distributed mesh. Two flat regions are connected with a cylindrical part. With the standard volume strain field, all elements (except one row) would gather at the curved region since the error is zero in the flat parts. One remedy to 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1539 ADAPTIVE MESH-UPDATING METHODS this problem is to accept a maximum angle error, h . The volume strain field for this case is expressed as h» G G a ¹ " 1!h » , hG'h (33) G 0, h )h G When the angle error is below h everywhere, the nodal motion stops. Another approach is to use a sum of two terms in the weight factor (cf. equation (22)), where the second term is equal to one. This volume strain field can be formulated as h» » G #(1!k ) G (34) a ¹ "1! k h » » G Now the elements will tend to have equal volume in the region where h is zero. Near areas of non-zero error the mesh density will gradually change. Results are shown in Figure 9(b). In this case, the value of » will influence the size of elements in the curved region. Elements in this region will have a volume close to » . A more detailed example with this volume strain field is given below. 6.3. Pinching of a cylinder This example consists of a regular cylinder that is clamped at both ends and with diametrical opposite point forces applied. The upper-half of the structured is modelled, taking advantage of one symmetry plane. The cylinder is modelled with shell elements and a linear elastic material model is used. The problem is solved taking large displacements into account in an updated Lagrangian formulation. The structure is modelled with a regular mesh consisting of 800 four-noded shell elements. The loading is applied as forced displacement at the point where the point force would be applied. A solution without adaptivity (regularly distributed elements) is shown in Figure 10, which also shows the geometry and material parameters. Taking the angle between elements as an error measure, the region of highest error is seen in this figure for a prescribed displacement of 0)34. We first apply the node relocation algorithm with the ‘standard’ formulation for the volume strain field: h» a ¹ "1! G G (35) G h » Here the maximum angle is chosen for h and » is equal to the original element volume. The mesh is optimized after each load step. The result after 7 load steps is shown in Figure 11. The node relocation is repeated until the difference between minimum and maximum h » is less that G G 5 per cent. In this case, the maximum angle between elements is about 8° as compared with the 11° without adaptivity. If a more pronounced node clustering is wanted, one could use the angle measure alone in the expression for the volume strain field: h a ¹ "1! G G h 1998 John Wiley & Sons, Ltd. (36) Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1540 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD Figure 10. A solution without adaptivity Figure 11. Standard formulation for the volume strain field Figure 12. Results with the angle measure used alone With the same strategy as above, the mesh in Figure 12 is obtained. Now the nodes are more clustered, and there is also a higher degree of clustering towards the point where the load is applied. The node lines describe what seems to be lines of constant curvature. This case could benefit from a formulation where elements with a bad aspect ratio are split into two elements. It often happens that certain elements get stretched due to the node relocation. If 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1541 ADAPTIVE MESH-UPDATING METHODS Figure 13. A sheet drawn over a cylinder this reduces the accuracy of the solution, one could consider augmenting the current approach with element splitting. This approach is not directly applicable with the currently used formulation for obtaining the rigid contact surfaces, but with the local smoothing approach it is. A node relocation step is less expensive than a step in the main analysis. Both analyses are non-linear, the main analysis due to large deformations and the relocation analysis due to contact. The contact problem is very simple (no new contact, separation or friction) and membrane elements are used instead of more expensive shell elements. The difference is much more pronounced in, for instance, a sheet metal-forming analysis which includes contact with friction, large deformations and large strain plasticity. 6.4. Sheet drawn over a cylinder As a precursor to later applications in sheet metal-forming, we here consider an elastic shell which is drawn over a rigid cylinder with zero friction as shown in Figure 13. The upper and lower edges of the shell are pinned and other two are free. The cylinder is first moved 0)04m in the x-direction and then rotated 10° around the x-axis. For the volume strain field we utilize the combined form h» » G #(1!k ) G a ¹ "1! k G h » » (37) since the regions not in contact with the cylinder will have h close or equal to zero. Element shape control is also used in this example. The mesh is relocated after each of the loading steps, using a sufficient number of increments in the auxiliary problem so that no further relocation takes place (which took about eight increments the first time it was applied and two or three increments 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1542 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD Figure 14. Results from sheet drawn over cylinder: (a) original mesh after first load step; (b) mesh after node relocation; (c) relocated mesh after second load step in subsequent steps). » was chosen as 60 per cent of the original element volume and for h the largest angle in the mesh was used. The results in Figure 14 clearly show how the mesh is optimized with respect to the angle measure. Of particular interest is the mesh after the second load step. The mesh becomes aligned along the new cylinder axis direction, although slightly hindered by the element shape control, especially near the edges. The behaviour of the proposed node relocation algorithm in this example makes it suitable for later applications to sheet metal-forming. 7. CONCLUSIONS In this work we have presented a method for continuous adaptive updating of computational meshes, designed primarily for use in incremental loading of non-linear shell and solid structures. The method is based on relocating the nodes of a pre-existing mesh with the topology preserved. This is achieved by formulating and solving an auxiliary (separate) thermoelastic problem where the temperature field is used to impose a proper distribution of thermal strains to obtain the desired mesh properties. The thermal strain field might be based on any local error indicator including other solution features like effective stress, contact tractions or purely geometric element characteristics. In particular, the angle equidistribution between elements was found to be effective and could be directly connected to the development of characteristic lines on the shell surface that relate to constant curvature bending lines. A limitation in the current formulation is that the mesh must be locally structured because of the algorithm used for constructing the rigid contact surfaces from the current geometry. However, this is relatively easily replaced by a local smoothing procedure which allows unstructured meshes. The method can be implemented easily in almost all non-linear finite element codes and can be extended to elastoplastic analysis as well (will be addressed in a forthcoming paper). Other methods of mesh adaptivity like the p- and the h-method can readily supplement the presented one. Our investigation has shown several good features of the proposed mesh-updating method. For example, 1. it is generally applicable to solid and shell problems 2. mesh integrity and robustness are built in (no overlapping elements possible) 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) ADAPTIVE MESH-UPDATING METHODS 1543 3. the smoothing properties are easily monitored 4. the basic algorithm is available in most non-linear finite element codes 5. a good intuitive physical interpretation of the mechanism behind the movement of the nodes. ACKNOWLEDGEMENT The financial support from ABB and the Swedish Research Council for Engineering Sciences is gratefully acknowledged. REFERENCES 1. T. J. R. Hughes, W. K. Liu and T. K. Zimmermann, ‘Lagrangian—Eulerian finite element formulation for incompressible viscous flows’, Comput. Meth. Appl. Mech. Engng., 29, 329—349 (1981). 2. J. Donea, P. Fasoli-Stella and S. Giuliani, ‘Lagrangian and Eulerian finite element techniques for transient fluidstructure interaction problems’, ¹rans. 4th Int. Conf. on SMIR¹, Vol. B, 1977, pp. 1—12. 3. W. K. Liu, T. Belytschko and H. Chang, ‘An arbitrary Lagrangian—Eulerian finite element method for pathdependent materials’, Comput. Meth. Appl. Mech. Engng., 58, 227—245 (1986). 4. W. K. Liu, H. Chang, J.-S. Chen and T. Belytschko, ‘Arbitrary Lagrangian—Eulerian Petrov—Galerkin finite elements for nonlinear continua’, Comput. Meth. Appl. Mech. Engng., 68, 259—310 (1988). 5. S. Ghosh and N. Kikuchi, ‘An arbitrary Lagrangian—Eulerian finite element method for large deformation analysis of elasto-viscoplastic solids’, Comput. Meth. Appl. Mech. Engng., 86, 127—188 (1991). 6. J. Huétnik and P. T. Vreede, ‘Progress in mixed Eulerian—Lagrangian finite element simulation of forming processes’, Int. J. Numer. Meth. Engng., 30, 1441—1457 (1990). 7. C. A. Felippa, ‘Optimization of finite element grids by direct energy search’, Appl. Math. Modelling, 1, 93—96 (1976). 8. D. J. Turcke and G. M. McNeice, ‘Guidelines for selecting finite element grids based on an optimization study’, Comput. Struct., 4, 499—519 (1974). 9. A. M. Winslow, ‘Numerical solution of the quasilinear Poisson equation in a nonuniform triangular mesh’, J. Comput. Phys., 2, 149—172 (1967). 10. J. U. Brackbill and J. S. Saltzman, ‘Adaptive zoning for singular problems in two dimensions’, J. Comput. Phys., 46, 342—368 (1982). 11. A. E. Giannakopoulos and A. J. Engel, ‘Directional control in grid generation’, J. Comput. Phys., 74, 422—439 (1988). 12. A. R. Diaz, N. Kikuchi and J. E. Taylor, ‘A method of grid optimization for finite element methods’, Comput. Meth. Appl. Mech. Engng., 41, 29—45 (1983). 13. N. Kikuchi, ‘Adaptive grid-design methods for finite element analysis’, Comput. Meth. Appl. Mech. Engng., 55, 129—160 (1985). 14. T. Torigaki, ‘Contact and impact problems using adaptive finite element methods’, Ph.D. Dissertation, The University of Michigan, 1989. 15. A. Saouab and D. Vandromme, ‘Application of a variational method in compressible flow computations’, in A. S. Arcilla et al. (eds.), Numerical Grid Generation in Computational Fluid Dynamics and Related Fields, Elsevier, North-Holland, 1991. 16. J. Cabello, R. Löhner and O.-P. Jacquotte, ‘A variational method for the optimization of two- and three-dimensional unstructured meshes’, AIAA-92-0450, 30th Aerospace Sciences Meeting & Exhibit, January 6—9 Reno, NV, U.S.A., 1992. 17. I. Babus\ ka and W. Rheinboldt, ‘Error estimates for adaptive finite element computations’, SIAM J. Numer. Anal., 15, 736—754 (1978). 18. I. Babus\ ka and W. Gui, ‘Basic principles of feedback and adaptive approaches in the finite element method’, Comput. Meth. Appl. Mech. Engng., 55, 27—42 (1986). 19. O. C. Zienkiewicz and J. Z. Zhu, ‘A simple error estimator and adaptive procedure for practical engineering analysis’, Int. J. Numer. Meth. Engng., 24, 337—357 (1987). 20. T. Belytschko, B. L. Wong, E. J. Plaskacz, ‘Fission—fusion adaptivity in finite elements for nonlinear dynamics of shells’, Comput. Struct., 33, 1307—1323 (1989). 21. E. Stein, B. Seifert, S. Ohnimus and C. Carstensen, ‘Adaptive finite element analysis of geometrically non-linear plates and shells, especially buckling’, Int. J. Numer. Meth. Engng., 37, 2631—2655 (1994). 22. J. Bonet, ‘Error estimators and enrichment procedures for the finite element analysis of thin sheet large deformation processes’, Int. J. Numer. Meth. Engng., 37, 1573—1591 (1994). 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998) 1544 P. NORDLUND, A. E. GIANNAKOPOULOS AND B. HA®GGBLAD 23. K. M. Okstad and K. M. Mathisen, ‘Towards automatic adaptive geometrically non-linear shell analysis, Part I: implementation of an h-adaptive mesh refinement procedure’, Int. J. Numer. Meth. Engng., 37, 2657—2678 (1994). 24. P. J. G. Schreurs, F. E. Veldpaus and W. A. M. Brekelmans, ‘Simulation of forming processes, using the arbitrary Eulerian—Lagrangian formulation’, Comput. Meth. Appl. Mech. Engng., 58, 19—36 (1986). 25. R. Löhner and C. Yang, ‘Improved ALE mesh velocities for moving bodies’, Comm. Numer. Meth. Engng., 12, 599—608 (1996). 26. Spectrum Solver, Solution Methods Theory, Centric Engineering Systems, Inc., 1993. 27. I. Babus\ ka and W. Rheinboldt, ‘A posteriori error analysis of finite element solutions for one dimensional problems’, SIAM J. Numer. Anal., 18, 565—589 (1981). 28. L.-Y. Li and P. Bettess, ‘Notes on mesh optimal criteria in adaptive finite element computations’, Comm. Numer. Meth. Engng., 11, 911—915 (1995). 29. A. M. Cuitin o and M. Ortiz, ‘Computational modelling of single crystals’, Modelling Simul. Mater. Sci. Engng., 1, 225—263 (1992). 30. C. J. M. Gelten and A. W. A. Konter, ‘Application of mesh-rezoning in the updated Lagrangian method to metal forming analyses’, in J. F. T. Pittman et al. (eds.), Numerical Methods in Industrial Forming Processes, Pineridge Press, Swansea, U.K., 1982, pp. 511—521. 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 43, 1523—1544 (1998)

1/--страниц