close

Вход

Забыли?

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

?

373

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