Products: ABAQUS/Standard ABAQUS/Explicit
This section describes the formulation of the quadrilateral finite-membrane-strain element S4R, the triangular element S3R and S3 obtained through degeneration of S4R, and the fully integrated finite-membrane-strain element S4.
At a given stage in the deformation history of the shell, the position of a material point in the shell is defined by
and other lowercase Greek subscripts which describe the quantities in the reference surface of the shell range from 1 to 2. In the above equation 
is the normal to the reference surface of the shell. The gradient of the position is 
with respect to 
. Note that in the above 
are local surface coordinates that are assumed to be orthogonal and distance measuring in the reference state. 
is the coordinate in the thickness direction, distance measuring and orthogonal to in the reference state. The thickness increase factor is assumed to be independent of .In the deformed state we define local, orthonormal shell directions 
such that
is the Kronecker delta and 
is the identity tensor of rank 2. Summation convention is used for repeated subscripts. The in-plane components of the gradient of the position are obtained as 
In the original (reference) configuration we denote the position by 
(
for the reference surface) and the direction vectors by 
, which yields
The position of the points in the shell reference surface is described in terms of discrete nodal positions with parametric interpolation functions 
. The functions are 
continuous, and 
are nonorthogonal, nondistance measuring parametric coordinates. For the reference surface positions one, thus, obtains
are 
Now consider the original configuration. The unit normal to the shell reference surface is readily obtained as
and distance measuring coordinates along these vectors. The derivatives of these coordinates with respect to follow from 
with respect to is readily obtained by inversion: 
with 
the original curvature of the reference surface.It is convenient to define the inverse of the reference surface deformation gradient
along the base vectors 
, since 
: 
The equations given in the earlier sections are valid for any local coordinate system defined in the current state. The vectors at the beginning of the analysis are determined following the standard ABAQUS conventions. In this section, we outline the way in which the in-plane coordinates are made corotational.
To obtain the updated version of , we follow a two-step approach. First, we construct orthogonal vectors 
tangential to the surface (following ABAQUS conventions). Subsequently, we calculate
to the vectors: : 
is to be determined such that the resulting deformation tensor is symmetric, as 
We assume that the nodal spin will be interpolated with the interpolation functions . During an increment the nodal spin is assumed to be constant; consequently, the value of the spin at each material point will be constant. Hence, we can use the same interpolation functions for the incremental finite rotation vector 
:
and a rotation axis 
: 
, which is defined by 
, which will again be chosen perpendicular to the reference surface. The updated normal used here will be approximately orthogonal to the reference surface, depending upon the amount of transverse shear deformation. The gradient of the updated shell normal can be obtained by differentiation: 
is equal to its conjugate (
). Hence, we can write 
follows 
, again with use of the incremental quaternion definition 
and follows 
and some algebra one obtains 
when 
.For the gradient 
of the updated shell normal we obtain
: 
is independent of 
.Calculation of involves taking the gradient with respect to the reference configuration. It is more convenient to use the reference surface curvature tensor
We already have obtained an expression for the deformation gradient in the reference surface, and we have assumed that the thickness change is constant:
, which yields the simplified relation 
The membrane strain increment follows from the incremental stretch tensor 
, whose components follow from the incremental deformation gradient 
by the polar decomposition 
.
Let 
and 
be the deformation gradient at the beginning and the end of the increment, respectively. By definition 
. The incremental deformation gradient follows as
are the components of an orthogonal matrix, the square of the incremental stretch tensor can be obtained by 
Following Koiter-Sanders shell theory, and compensating for the rotation of the base vectors relative to the material, we define the physical curvature increment 
as
relative to 
, this expression can be rewritten as 
, does not appear in this equation. Hence, there is no need to calculate the initial curvature 
, and we can assume 
. The deformation gradient can, hence, also be simplified to 
The virtual work contribution of the stresses is
and the bending moments 
: 
.The variation in the curvature is obtained by taking variations in the incremental curvature, which yields
and also terms of order 
, which yields 
with respect to the current state (at the end of the increment). Hence for the evaluation we can assume 
. Moreover, we neglect terms of the order 
since they are proportional to . Hence, we obtain 
yields 
To obtain an expression for the rate of virtual work, we first write the virtual work equation in terms of the reference volume
is the Kirchhoff stress tensor, related to the Cauchy or true stress tensor via 
indicates that the rates are taken in a material, corotational coordinate system. The terms involving stress rates are related to the material behavior. We assume constitutive equations of the form 
and transformed back to the current configuration, this yields 
. Hence, the rate of virtual work can be written as 
It remains to determine 
and 
. From the first variation 
follows
is the inverse of 
, it follows that 
yields 
We need to calculate the second variation of the curvature to calculate the initial stress contribution from the curvature:
are the components expressed in the orthonormal coordinate system (
) transformed by 
.Denoting derivatives with respect to the isoparametric coordinates as 
, the second variation of the curvature is
and 
, we find that 
indicates the skew-symmetric tensor with axial vector .Several interpolation schemes have been proposed to avoid shear-locking, which typically arises as the thickness of a plate or shell goes to zero. Here we employ an assumed strain method based on the Hu-Washizu principle. This scheme derives from that by MacNeal (1978), subsequently extended and reformulated in Hughes and Tezduyar (1981) and MacNeal (1982) and revisited in Bathe and Dvorkin (1984). Computational aspects of the nonlinear theory are investigated in Simo, Fox, and Rifai (1989) for fully integrated quadrilateral shell elements. For reduced integration quadrilateral and triangular shell elements that can be used for both implicit and explicit integration, this assumed strain method needs to be modified. We summarize below the assumed strain method used with fully integrated elements, followed by the modifications required for the one-point integration plus stabilization used in ABAQUS.
Consider a typical isoparametric finite element, as depicted in Figure 3.6.5–1, and denote by 
the set of midpoints of the element boundaries.
The following assumed transverse shear strain field is used:
, for 
, are the reference surface position vectors of the element nodes.By making use of the assumed strain field along with the update formulae for the director field, the assumed covariant transverse shear field can be written concisely in matrix notation. Recall the director field update equation and the corresponding linearized director field:
A St. Venant-Kirchhoff constitutive model for the Kirchhoff curvilinear components of the resultant transverse shear force is written in terms of the transverse shear strains as
is the transverse shear stiffness in curvilinear coordinates. For a single isotropic layer, 
is the inverse of the metric 
, where metric components in the reference configuration 
are defined by the inner product 
are calculated with the coordinate transformation 
as 
The calculation of the initial stress stiffness matrix requires the second variation of the assumed transverse strain field. This calculation can be summarized in matrix notation as follows. Define vectors of variations of the nodal displacement quantities:
is the area measure in the current configuration and 
is the (symmetric) transverse shear contribution to the initial stress, defined as follows. Let be the 
identity matrix; then define the symmetric matrices 
be the zero matrix. Then is written 
For reduced-integration elements the transverse shear force components need to be evaluated at the center of the elements. Consider 
the transverse shear contribution to the internal energy:
is written in terms of the isoparametric coordinates as 
, where 
and are the components of the reference surface metric in the undeformed configuration.This transverse shear energy can be approximated in many ways to produce a one point integration at the center of the element plus hourglass stabilization. It is important that this treatment yield accurate representation of transverse shear deformation in thick shell problems and provide robust performance for skewed elements. The treatment should collapse smoothly to a triangle, which should be insensitive to the node numbering during collapse; that is, the triangle's response should not depend on the nodal connectivity. For an entire mesh of triangular elements, the treatment should give convergent results (that is, the element should not lock). Furthermore, the high frequency response of the transverse shear treatment should be controlled so that transverse shear response does not dominate the stable time increment for explicit dynamic analysis (including for skewed triangular or quadrilateral geometries). All of these requirements are embodied in the following transverse shear treatment.
Define the transverse shear strain at the center of the element (the homogeneous part) and the “hourglass” transverse shear strain vectors as
and 
are constants determined by the element reference geometry. For geometries with constant Jacobian transformation, 
. The components of the hourglass strain vector 
are defined in terms of the edge strains as 
, 
, 
, and 
are constants determined from the reference geometry of the element. For rectangular elements 
, 
, 
, 
and 
can be identified as the strain associated with the rotational “butterfly” deformation pattern. We call 
the “crop circle” mode strain since it corresponds to a deformation pattern that resembles the sweeping over the element normals in a circular pattern.The inclusion of the crop circle strain in the homogeneous part of the transverse shear strain 
has two important consequences. First, it makes the transverse shear response insensitive to the nodal connectivity for a triangular element. That is, when a side of a quadrilateral element is collapsed to form a triangle, the element's response is independent of the choice of node numbering on the element. Second, for explicit dynamic analyses the coefficients and are chosen to minimize the highest frequencies associated with the homogeneous part of the transverse shear response.
To illustrate the crop circle and butterfly transverse shear patterns, consider a square, initially flat element. Furthermore, consider plate theory kinematics; that is, two rotations and a vertical deflection at the nodes. The crop circle pattern has zero vertical deflection at the nodes and a nodal rotation vector pattern as illustrated in Figure 3.6.5–2.
The butterfly pattern has vertical deflections that correspond to cross-diagonal bending; that is, two equal deflections at two nodes across a diagonal, with equal and opposite deflections at the remaining two nodes. The nodal rotations develop in a way that opposes the bending motion of the reference surface; that is, the rotations are opposite the rotations that would develop for this displacement pattern to produce pure bending. The butterfly mode's nodal vertical deflection and rotation vector pattern are illustrated in Figure 3.6.5–3.Let the reference element area be
. The transverse shear energy can be approximated as a center point value plus a stabilization term: 
is the transverse shear stiffness evaluated at the center of the element and the hourglass stiffness 
is the diagonal matrix 
is the average direct component of the transverse shear stiffness, 
.The formulation of the homogeneous part of the transverse shear has two contributions: the average edge strain across the element, plus the element distortion term. The average strain treatment is essentially the same as that for the assumed strain formulation of MacNeal and others presented earlier, with expressions evaluated at the center of the element (
and 
). The details of this part are omitted; only the element distortion term is presented in detail. The variation of the homogeneous transverse shear strain can be written
and 
are 
and 
evaluated at the center of the element, 
The stabilization term has a similar formulation. The variation of the hourglass strain is
The hourglass force components 
and 
are given by the constitutive relations
(1) The butterfly mode is applied with a “large” or physical hourglass stiffness. For a reference geometry with constant Jacobian, the butterfly stabilization term can be derived from an exact integration of the assumed strain formulation of the transverse shear energy. It is important to apply this constraint with a high stiffness to prevent overly flexible response for quadrilateral elements. The crop circle mode is applied with a “small” or weak stiffness. Although this mode can propagate, it is rarely problematic and is often prevented with boundary conditions.
(2) As the quadrilateral element is degenerated to a triangle, the two hourglass constraints converge into a single constraint: the crop circle constraint. However, as is well-known, for a constant strain triangle the element will lock for certain meshes with three transverse shear constraints per element. Therefore, in the case of a triangular element, the (strong) butterfly mode stabilization is not applied. Only the (weak) crop circle mode stabilization is applied. Thus, in addition to the two homogeneous transverse shear strains, the triangle has a weak constraint to prevent spurious zero energy modes, yet avoids locking in most situations.
The initial stress contribution from the stabilization terms takes the following form:
is the (symmetric) transverse shear stabilization contribution to the initial stress. Define the symmetric matrices 
is written 
are defined, is filled just as .The initial stress contribution from the homogeneous part consists of two terms, one from the assumed strain formulation (evaluated at 
) as detailed earlier, and the other from the crop circle mode addition. These two terms can be written
and 
are the shear force and matrix evaluated at the element center. The matrix expression for 
is analogous to from the stabilization terms.The in-plane displacement hourglass control is applied in the same way as in the ABAQUS membrane elements. The hourglass strains are defined by
is the hourglass mode. This mode is obtained by making the “regular” hourglass mode 
orthogonal to the homogeneous deformation mode in the undeformed shape of the element. This last condition can be written as 
is a small number chosen to be 0.005 in ABAQUS/Standard and 0.05 in ABAQUS/Explicit. When the hourglass control is based on assumed enhanced strain, the artificial stiffness factor is replaced by coefficients derived from a three-field variational principle. The hourglass force Z conjugate to z is then equal to 
. The second variation is needed for the Jacobian. From the first variation follows right away 
.The expressions for the curvature change, the transverse shear constraints, and the drilling mode constraints still leave three nonhomogeneous rotational modes unconstrained. These modes correspond to zero rotation at the midedges and zero gradient at the centroid. Hence, they correspond to the familiar hourglass pattern. To pass curvature patch tests exactly, it is necessary to use orthogonalized hourglass patterns as derived for in-plane hourglass control.
This last aspect implies that the rotational hourglass mode corresponds to the mixed derivative of the rotation at the centroid:
We cannot use the above formulation directly in a formulation suitable for multiple finite rotation increments. Hence, we use the same approach as for the calculation of the curvature change. For the purpose of the calculation we define the updated shell direction vectors
and 
can be neglected. Consequently, 
can be expressed in terms of the hourglass strain tensor with 
. Subsequent values of 
are obtained incrementally. From the expression for 
we obtain 
and . Hence, the above expression simplifies to 
is equal to its conjugate (); hence, we can write 
follows, while neglecting the products of and derivatives: 
again with use of the incremental quaternion definition 
and follows, again neglecting the products of first derivatives 
and some algebra one obtains 
when . For the updated hourglass tensor one readily obtains 
and obtain 
and 
with as final result 
In general meshes it will be desirable to collapse at least some of the quadrilateral elements to triangles or to use the triangular element S3 or S3R, which is in fact an internally collapsed S4R element. For this case the calculation of the membrane strains and the curvature changes proceeds along the same lines as before. The transverse shears will now be zero at the degenerate edges. Finally, calculation of all hourglass constraints will be omitted.
For numerical efficiency in explicit dynamic analysis, it is desirable to have the stable time increment determined by the membrane response of the structure. For this reason scaling of the rotary inertia based on the element's reference geometry is included in ABAQUS/Explicit.
In explicit dynamic analyses the stable time increment is proportional to the inverse of the highest frequency of the element. Therefore, we must ensure that the highest frequency associated with the transverse shear response does not exceed the highest frequency associated with the membrane response. For thick elements (that is, for elements whose thickness is order unity relative to a characteristic length in the element), the membrane frequencies are dominant. The primary consideration in choosing appropriate scalings is that in the limit as the element's thickness goes to zero, the transverse shear frequencies remain below the membrane ones. Recall that for a one-dimensional spring-mass oscillator, the natural frequency 
can be written in terms of the stiffness K and the mass M as
For thin elements the rotary inertia is small (negligible) relative to the rotational inertia of the mass at the nodes rotating about an axis through the element center. Therefore, we choose a scaling on the rotary inertia such that it never becomes smaller than a fixed (small) percentage of the rotational inertia of the mass at the nodes rotating about an axis through the center of the element.
Let R be the nondimensional rotary inertia scaling, where 
. When 
, the true rotary inertia is used. Consider a lumped mass matrix for a 4-node element, and let the element be flat and square. For rotations about an axis in the plane of the element, parallel to an element edge, passing through the center, the contribution to the rotational inertia of the mass at the nodes is
is the mass density. The sum of the rotary inertia at the four nodes is 
—of the mass contribution as the thickness goes to zero, asymptotically R must be proportional to 
; that is, 
For planar geometries with element directors along the normal direction, closed-form expressions are possible for the highest membrane and transverse shear frequencies. In such cases the length parameter L is interpreted as a characteristic element length that depends on the element distortion. To handle arbitrarily shaped curved elements, exact calculation of the element frequencies becomes difficult. However, we can safely bound the frequencies by an appropriate choice of L in the following scaling:
. For triangular elements the characteristic length, 
is the minimum element edge length; and for quadrilateral elements, 
For the displacement degrees of freedom, bulk viscosity introduces damping associated with volumetric straining. Linear bulk viscosity or truncation frequency damping is used to damp the high frequency ringing that leads to unwanted noise in the solution or spurious overshoot in the response amplitude. For the same reason, in shells we need to damp the high frequency ringing in the rotational degrees of freedom with linear bulk viscosity acting on the mean curvature strain rate. This damping generates a bulk viscosity “pressure moment,” m, which is linear in the mean curvature strain rate:
is the original thickness, is the mass density, 
is the current dilatational wave speed, L is the characteristic length used for rotary inertia scaling, and 
is twice the increment in mean curvature. The dilatational wave speed is given in terms of the effective Lamé's constants as 
, where h is the current thickness, is added to the direct components of the moment resultant.Element S4 is a fully integrated finite-membrane-strain shell element. Since the element's stiffness is fully integrated, no spurious membrane or bending zero energy modes exist and no membrane or bending mode hourglass stabilization is used. Drill rotation control, however, is required. Element S4 uses the same drill stiffness formulation as used for element S4R. Similarly, element S4 assumes that the transverse shear strain (and force, since the transverse shear treatment is elastic based on the initial elastic modulus of the material) is constant over the element. Therefore, all four stiffness integration locations will have the same transverse shear strain, transverse shear section force, and transverse shear stress distribution. The transverse shear treatment for S4 is identical to that for S4R.
It is well known that a standard displacement formulation will exhibit shear locking for applications dominated by in-plane bending deformation. However, a standard displacement formulation for the out-of-plane bending stiffness is not subject to similar locking response. Hence, S4 uses a standard displacement formulation for the element's bending stiffness, and the theory presented above for the rotation kinematics and bending strain measures applies to S4. The primary difference between the element formulations for S4 and S4R is the treatment of the membrane strain field. This formulation is the topic of the following discussion.
The membrane formulation used for S4 does not rely on the fact that S4 is a shell element. Hence, the discussion below details the formulation from the point of view that the membrane response is governed by the equilibrium for a three-dimensional body in a state of plane stress.
Consider an enhancement 
to the rate of deformation tensor 
. We introduce the enhanced rate of deformation tensor, 
, as
is defined subsequently.Admissible variations in the rate of deformation are also introduced as
We now introduce constraints on the enhancements and 
:
is the specified traction on 
and 
on 
. 
is an arbitrary stress field, and the constitutive equation 
() is enforced pointwise.In the modified virtual work statement all kinematic quantities and corresponding variations (, , 
, and 
) are known functions of 
, 
, and the reference configuration. A fundamental requirement for the validity of the formulation is that the modified virtual work statement leads to the proper equilibrium equations. If is arbitrary, the constraint equations can be rewritten as
. We recognize this variational statement as the usual virtual work equation, and a straightforward application of the divergence theorem leads to the standard equilibrium equations.In the actual implementation we choose to satisfy the constraints only for piecewise constant stress fields . Hence, over the element domain 
we require
and 
are chosen such that they eliminate the shear locking for in-plane bending. In addition, the direct strain field is enhanced to approximate the strains caused by Poisson's effect in bending.To pass the patch test, the choice of enhancements cannot be arbitrary. A sufficient condition for the satisfaction of the patch test is that for homogeneous deformations we have 
or 
pointwise. In that case 
and the stress is homogeneous. Since the stress is homogeneous, it can be moved outside the volume integral in the modified virtual work statement. The volume integral condition on 
implies that the expression is independent of the enhancement and leads to the standard displacement formulation, which is known to satisfy the patch test.
