EXPLICIT SIMULATIONS WITH REDUCED INTEGRATION SOLID-SHELL ELEMENTS : STABILIZATION AND SELECTIVE MASS SCALING

Explicit approaches are usually preferred for the simulation of thin walled structural problems, which are often highly nonlinear due to large deformations and possible material inelasticity. Solid-shell elements can describe the correct thickness geometry and are therefore more suitable than standard shell elements for the implementation of complex 3D material models. However, they include a too simple kinematic formulation leading to artificial stiffening phenomena called locking. To overcome this problem computationally expensive corrections, e.g. introducing enhanced strains, are required which suggest the use of reduced integration with hourglass stabilization. Furthermore, a high element maximum eigenfrequency is implied by the small thickness, leading to overly small stable time-steps. These two issues are addressed in this paper where a stabilized reduced integration solidshell element and a selective mass scaling technique for the reduction of the maximum eigenfrequency are proposed.


INTRODUCTION
The finite element simulation of shell or thin walled structures is of interest in many industrial applications, such as naval and aerospace engineering, metal forming processes and automotive industries.Due their morphology, this type of structures may often undergo large inelastic deformations for which explicit simulation approaches become preferable, especially in the presence of other sources of nonlinearities (like e.g.contact, crack propagation, delamination).In recent years, considerable effort has been devoted to the development of 3D continuum-based elements for shell applications, the so called solid-shell elements (see, for instance, [1], [2] and [3]).These elements have only translational degrees of freedom and are characterized by a fully three-dimensional stress state which makes them more suitable for the implementation of complex 3D material models.However, while the use of classical shell elements in explicit dynamics simulations is well established, solid-shell elements have been predominantly used in the framework of implicit formulations.In fact, the kinematic formulations of compatible solid-shell elements is too simple and requires corrections e.g.introduced via the enhanced or assumed natural strain method.However, these extended element formulations may become computationally too expensive for explicit approaches where very small time steps are required for stability reasons.A possible remedy consists in adopting a reduced integration which, in turn, requires an hourglass stabilization.As a further problem, the small thickness of solid-shell elements leads to a very small travel time for propagating elastic stress waves, requiring, as a consequence, overly small stable time-steps in explicit time-integrations, with prohibitive computational costs.This problem can be circumvented making use of a selective mass scaling technique, whereby only masses related to some of the degrees of freedom are increased, leaving the inertia connected with the element rigid body modes unaltered.The goal of the present paper is to propose an explicit finite element methodology for the simulation of the large inelastic deformation of shell structures.The solid-shell finite element recently proposed in [2] in the framework of an implicit approach is here reconsidered in an explicit context and a new selective mass scaling technique, specifically conceived for this type of elements is presented.

ELEMENT FORMULATION
Low-order continuum elements, from which solid-shell elements are usually derived, exhibit several types of locking behavior.In the reduced integration solid-like shell element Q1STs proposed in [2], transverse shear and curvature thickness locking are cured by means of the assumed natural strain concept while the enhanced strain method is used to avoid volumetric and Poisson thickness locking.The formulation of the element is briefly summarized below using a total Lagrangian approach which does not require the definition of a co-rotational coordinate system to maintain material objectivity.

Variational formulation
The explicit framework of solid-shell element Q1STs is derived from the dynamic variational functional: where δ denotes a variation, ρ the material density, u the displacement field and S(E) is the second Piola-Kirchhoff stress tensor.It is a function of the Green-Lagrange strain tensor E, which is additively decomposed into the displacement-based part E c and the enhanced strain part E e according to the EAS concept.Using the orthogonality between stress and enhanced strain fields, the stress field can be eliminated from the variational functional (1), (2).The Q1STs element has eight nodes, one in-plane integration point and the ability to accommodate several integration points along the thickness direction of the element.Fig. 1 shows the element topology with the node numbering order (related to the isoparametric coordinate system of the element, defined by the convective natural coordinates ξ i , i = 1, 2, 3).The covariant base vectors in the undeformed configuration are written as where X(ξ) = [X 1 , X 2 , X 3 ] T denotes the position vector in the reference configuration and J i represent the columns of the Jacobian matrix J = [J 1 , J 2 , J 3 ].The covariant base vectors in the deformed configuration are given by where x denotes the position vector in the current configuration, Ji are the columns of the Jacobian matrix with respect to the current configuration and D = ∂u/∂ξ represents the displacement gradient with respect to the vector of natural coordinates ξ.Accordingly, the contravariant base vectors in the reference and current configurations are defined as where j ij = (J −1 ) ij and jij = ( J−1 ) ij denote the coefficients of the inverse Jacobian matrix.
Starting from the definition, the deformation gradient F has the following form:

Strain tensor
The Green-Lagrange strain tensor E is used as a suitable strain measure in large deformation problems.We consider an additively decomposition of the total Green-Lagrange strain tensor E into the compatible part E c and the enhanced strain part E e .

Compatible strain
The compatible Green-Lagrange strain tensor has the following form: where I is the second order identity tensor.Applying (4), the covariant components E ij can be derived from where J i • D j and D i • D j represent the linear and nonlinear part of the strain.From the finite element point of view, the strain tensor is often represented in Voigt notation, storing the coefficients E ij into the 6 × 1 vector with Γ cξ i ξ j := 2E cξ i ξ j .

ANS concept
The transverse shear locking occurs in bending situations when the element thickness tends to zero.One of the most efficient approaches to cure this pathology is the assumed natural strain (ANS) method.This concept can be used also to cure the curvature thickness locking which occurs only in initially curved shell structures.The idea of the ANS method is to define a strain field interpolated independently of the displacement field, in terms of strains computed at locking free sampling points.In order to avoid the curvature thickness locking, the covariant compatible transverse normal strain E cζζ is evaluated at the corners of the element midplane (sampling points K = A, B, C, D) and interpolated by means of bilinear ansatz functions.On the other hand, to cure the shear locking the shear terms Γ cηζ and Γ cξζ are evaluated respectively at sampling points K = E, F, G, H and K = J, K, L, M (see Fig. 2) as suggested in [3].The compatible covariant strain terms assume the following form:

Hourglass stabilization
One key point of Q1STs element is the treatment of the hourglass modes.In fact, the use of a reduced integration scheme within the shell plane leads to the necessity of a stabilization procedure to avoid spurious singular modes.In the considered reduced integration scheme, all integration points are located along the normal through the center of the element defined by ξ z = (0, 0, ζ) T , so that the compatible Green-Lagrange strain can be split into where Ēz c = Ēc (ζ z ) is the part related to the integration points while Ēhg c = Ēc − Ēz c denotes the hourglass strain.The part Ēz c enters the constitutive law for the stress computation while the hourglass strain Ēhg c intervenes in the hourglass stabilization.Since this term does not have a direct physical meaning, it is desirable that it is computed as efficiently as possible.One possibility consists of analytically integrating the hourglass stabilization over the element domain.This could be easily achieved if the integrands were polynomials.For this purpose, a Taylor expansion of the compatible Green-Lagrange strain tensor is carried out with respect to the center of the element We can single out two contributions: one related to the integration points ( Ēz c ), constant within the shell plane, and another one bi-linear within the shell plane ( Ēhg c ), which is used as a basis for the construction of the hourglass stabilization.To correctly study problems involving thick shells, an expansion of Ēz c to a quadratic function is required.

Enhanced strain
Bending of a shell-like structure leads to a non-constant distribution of the transverse normal strain along the thickness.Consequently, in order to avoid the Poisson thickness locking in the element formulation, the transverse normal strain component has to be modeled at least linearly.Since E AN S cζζ is modeled as being independent of ζ, its polynomial order is not high enough.Furthermore, volumetric locking occurs when the material approaches incompressibility.To overcome this problem, the strain terms E ξξ , E ηη and E ζζ must be represented by polynomials of the same order.Again the term E ζζ has to be linear.It has been shown that the introduction of only enhanced-degree-of-freedom W e is sufficient to avoid these types of locking.The resulting covariant enhanced strain field E e , defined by the relation is constant within the shell plane and contributes only to the covariant strain field evaluated at the normal through the element center.

Transformation tensor
The strain components can be expressed in the global reference frame, starting from the covariant components, as where H a are the contravariant base vectors and e i represent the base vectors of a Cartesian coordinate system.Adopting Voigt notation for the representation of strain tensor, Eq. ( 18) can be written in matrix form as where T represents the transformation matrix.

2nd Piola Kirchhoff stress tensor
As done for the Green-Lagrange strain tensor, it is possible to carry out a Taylor expansion of the stress field with respect to the normal through the center of the element (ξ z = (0, 0, ζ) T ) in order to define a convenient hourglass stress Sghg which does not influence the physical element response and can be efficiently computed.
The hourglass stress Sghg is a function of the tangent C and of the cartesian compatible hourglass strains.To improve the efficiency of the hourglass stress, the non-linear tangent C can be replaced by the linear-elastic material tangent C hg , which is modified to yield a deviatoric hourglass stress.The latter assumption is necessary to overcome volumetric locking.The computation of stabilization forces at each time-step is computationally expensive and may be not necessary for the overall analysis stability [4].Criteria for the definition of the updating frequency of the stabilization force are currently under study.

SELECTIVE MASS SCALING
Explicit time integration is stable only if the used time-step is smaller than a critical threshold, which can be shown to depend on the smallest geometrical dimension of the finite elements in the mesh.This aspect is particularly critical when solid-shell elements are used for the analysis of thin walled structures, since the small thickness can lead to unacceptably small time-steps.To overcome this problem, we adopt a selective mass scaling technique, based on a linear transformation of the element degrees of freedom to increase the size of the critical time-step without affecting the dynamical response.It is well-known that the structural response in inertia dominated dynamic problems depends predominantly on inertia properties related to rigid body modes of individual finite elements.Consequently, in this kind of problems, better results can be obtained by selectively scaling element masses, in such a way that masses associated to element rigid body modes are not modified.The solid-shell dofs are linearly transformed to segregate middle plane displacement dofs controlling translational rigid body motions.Masses pertinent to relative dofs, controlling higher order modes, are selectively scaled.It is important to underline that the scaled mass matrix remains diagonal for effective use in explicit methods while the scale factor is estimated so as to obtain a critical time-step size of the same order of magnitude of the in-plane "traversal time".Introducing the definition of upper and lower surfaces as depicted in Fig. 3, we introduce the following variable transformation: where u upper and u lower indicate respectively the dofs of nodes belonging to the upper and lower surfaces.On the other hand ũMiddle and ũDiff represent the dofs of nodes on the element middle plane.In matrix notation: The lumped mass matrix associated to the middle surfaces dofs is given by It should be noted that T T M lump T is not diagonal even when the original mass matrix is so.
A further lumping procedure is therefore necessary to make M diagonal.For a parallepiped it takes the form M = ρV 4 where V is the element volume and α is the nondimensional mass scaling factor.A conceptually almost identical scaling procedure was presented in [5] where, rather than to the masses, the scaling was applied directly to nodal accelerations.For solid-shell elements, where the thickness dimension is significantly smaller than the in-plane ones, the highest element eigenfrequency always turns out to be given by the square root of the eigenfrequency ω 2 corresponding to the thickness vibration mode.In the case of a regular parallelepiped, in [6] it has been shown that the critical time step resulting from this eigenfrequency can be analytically computed as a and b being the element in-plane semi-dimensions, ν the Poisson coefficient, ρ the material density and E the Young modulus.The parameter η is computed starting from the cubic equation: resulting from the element eigenvalue problem.The coefficient γ and λ are defined as γ = c/a, λ = b/a while c represents the half-thickness.After suitable variable transformations and simplifications, the following solution of Eq. ( 28) is found: The scaling parameter α, appearing in the coefficients C 1 , C 2 and C 3 , has to be optimized so as to maximize the time-step size without significant accuracy losses.The variation of the critical time step ∆t with respect to the scaling parameter exhibits an asymptotic behavior (see Fig. 4).Since selective mass scaling has the effect to decrease the eigenfrequency corresponding to the thickness vibration mode, there is a critical value of α above which the critical time-step size becomes determined by the in-plane traversal time.As a consequence, values of α above this limit would be useless since they would not bring an increase of the critical time step, while leading to inaccurate solutions.A possible strategy for the determination of the optimal value of α consists first in the determination of the asymptotic value ∆t asy of the critical time step, which can be derived letting α → ∞ in (28).In this way, the following quadratic equation is obtained: The roots of (37) are given by The corresponding asymptotic value of the critical time step ∆t asy is then given by ∆t asy = 24 abρ(1 + ν) Eη asy η asy = max(η 1 , η2 ) As a trade-off between accuracy and computational efficiency, the optimal value α opt of the scaling parameter is chosen as the one which produces an estimate of the critical time step equal to 0.9∆t asy , since values of α higher than this limit do not produce significant increase of the time step.The value of α to be used for mass scaling is then obtained by iteratively solving Eq. ( 27) for α, having set ∆t = 0.9∆t asy .The above discussed procedure provides a viable and effective strategy for an optimized and fully automatic scaling procedure for solid-shell elements.

NUMERICAL RESULTS
A cantilever beam subject to a uniform pressure load is used to test the accuracy properties of the selective mass scaling procedure and the geometrically nonlinear response of the Q1STs solid-shell element.The cantilever is clamped on one side and subject to a uniform pressure load applied instantly and kept constant throughout the rest of the simulation.Geometrical and material data are referred to the example described in [7].The plate material is assumed to be linear elastic with Young's modulus E = 12000 psi (82.7 MPa), Poisson coefficient ν = 0.2 and material density ρ = 0.1224 × 10 −5 lb s 2 /in 4 .The beam, depicted in Fig. 5, is 10 in long, 1 in wide.A convergence study has been conducted varying the thickness of the plate (0.1 and 0.05 in) and the number of solid-shell elements (5, 10, 20 and 40 elements) along the plate length.For each simulation, the optimal value of α is computed according to the above described strategy.Table 1 summarizes the computed scaling parameter values and  the corresponding critical time step size.The time step value for the standard case (without mass scaling), computed by means of Gershgorin estimation, is also reported.Fig. 6 shows the time histories of the tip displacement for different meshes in the case of thickness t = 0.1 in and t = 0.05 in.It can be observed that also in the case of the coarsest mesh, requiring a high scaling parameter α, the dynamic response is not corrupted and the results are quite good.The present mass scaling procedure is based on the hypothesis that inertia properties related to deformation modes defined by the "difference" dofs have negligible effects on the structural response.In order to validate this assumption, Fig. 7 shows the trend of the kinetic energy associated to middle plane and difference dofs.As expected, the latter contribution results to be significantly smaller than the one due to middle plane dofs confirming that modifications of this contribution will have only minor effects on the overall dynamic behavior.

CONCLUSIONS
Solid-shell elements have attractive features for the simulation of thin-walled structures with complex material models.In view of applications in highly nonlinear problems, there is strong interest in their usage in the context of explicit approaches.The present paper considers the implementation in an explicit code of a recently proposed reduced integration solid-shell element, originally formulated in an implicit framework.An inexpensive selective mass scaling procedure is proposed to overcome the problem of the excessively small critical time-step size resulting as a consequence of the element small thickness.An automatic procedure for the optimal definition of the scaling parameter is also discussed.

Figure 1 .
Figure 1.Node numbers and Gauss points for the solid-shell element.

Figure 2 .
Figure 2. Sampling points for the application of the ANS concept.

Figure 3 .
Figure 3. Definition of lower and upper surfaces.

Figure 4 .
Figure 4. Trend of ∆t for varying value of α.