Efficient isogeometric shell element with through-thickness stretch: application to incremental sheet forming

An isogeometric shell element with through-thickness stretch is applied to a two-point incremental forming problem. The shell element supports full three-dimensional constitutive laws and therefore does not make the plane stress assumption. An anisotropic material model is implemented to account for the sheet rolling direction. Automatically adjusting penalty stiffness is proposed for modeling the contact between the stylus tool and the sheet, whereas the die contact algorithm uses traditional constant penalty stiffness. A comparison is made between experimental results as well as results from a conventional shell formulation.


Introduction
Highly nonlinear problems such as automobile crash dynamics or sheet metal forming simulations often utilize the four-node, bilinearly interpolated Reissner-Mindlin shell elements combined with explicit time integration. These Reissner-Mindlin elements imposing the plane stress condition have dominated the finite element shell analysis for decades. The use of bilinear interpolation functions has resulted in efficient formulations and larger stable time step sizes compared to higher-order Lagrange shell elements. However, the interest in more advanced formulations such as thickness stretchable shell elements and solid-shells has increased as the plane stress assumption is not always acceptable, particularly if thick plates or problems involving doublesided contact are considered. Many recent contributions show a great interest in using solidshell elements in explicit dynamic simulations. However, the computational speed has yet to be improved for demanding problems such as incremental sheet forming. Isogeometric analysis (IGA) first proposed by Hughes et al. [1] in 2005 has recently attracted growing interest among the computational mechanics community. This concept adopts splinebased interpolation functions resulting in several potential advantages over conventional finite elements, such as increased accuracy per degree of freedom, increased stable time step size, and continuous stress and strain fields, among others [2,3,4,5,6]. Isogeometric solid-shell elements have been proposed by Hosseini et al. [7], Bouclier et al. [8,9], and Caseiro et al. [10]. However, most of the isogeometric solid-shell element formulations proposed in the literature do not consider dynamic problems nor the implications of using explicit time integration.
The isogeometric shell element introduced by Hokkanen et al. [11] supports full three-dimensional constitutive laws and is capable of simulating geometrically nonlinear dynamic problems. The purpose of this paper is to evaluate the suitability of the aforementioned shell element to incremental sheet forming (ISF) using explicit time integration.

Shell formulation
The implemented formulation is based on the degenerated solid approach with only displacement degrees of freedom. Reduced integration is used to alleviate several locking problems and to reduce computational cost. The reduced Gauss 2 × 2 quadrature preserves full rank of the stiffness matrix and therefore no ad hoc hourglass control procedures are required. The fiber mass scaling prevents the eigenfrequencies related to the transverse deformations from lowering the larger stable time step size associated with the spline-based in-plane interpolation. The implementation relies on CUDA and most computations are performed parallel on a modern consumer level GPU. A detailed description of the shell element is given by Hokkanen et al. [11]. The BWC shell element [12] based on the plane stress assumption is used as reference.

Material model
The simulation assumes hypoelastic-plastic material behavior. The material model is based on Yld2004-18p yield criterion proposed by Barlat et al. [13] to consider the anisotropy introduced by the process of rolling the sheet. The required parameters describing the constitutive relations for aluminum 7075-O are the density ρ = 2.81×10 −9 Ns 2 mm 4 , Young's modulus E = 71.7×10 3 MPa, Poisson's ratio ν = 0.33, and the initial yield stress σ y = 261 MPa. The hardening coefficient and the hardening exponent required by the Swift type hardening law are given as C sw = 333 MPa and n sw = 0.16, respectively. Moreover, the 18 parameters required by Yld2004-18p yield criterion are given by Esmaeilpour et al. [14]. The return mapping uses the explicit closest point projection method [15].

Contact
The isogeometric formulation uses a non-mortar Gauss-point-to-surface (GPTS) penalty contact algorithm [16,17] for the stylus contact. The linear pressure-overclosure relationship is given as where ε is the penalty stiffness, g is the overclosure and s is the shift which causes the contact pressure to start increasing slightly before the actual contact occurs. The penalty stiffness ε (i.e. the slope of the pressure-overclosure relationship) adjusts automatically such that zero maximum penetration is maintained (i.e. g 0). However, a minimum value for ε is specified. The described contact formulation is found robust for the stylus contact in many ISF simulations as the contact between the sheet and the spherical-headed tool comprises mostly of a single contact region. The die contact algorithm uses non-mortar Gauss-point-to-surface (GPTS) penalty contact, a linear pressure-overclosure relationship, and traditional constant penalty stiffness.
Moreover, the traditional BWC shell formulation uses a simple node-to-surface (NTS) contact formulation and traditional constant penalty stiffness for the stylus and the die.

Results
The die geometry and the toolpath for the ISF simulation are given in Figures 1 and 2, respectively. The square sheet has the thickness of 0.063 inches and the edge length of 290 mm. The tool has a hemispherical head, 20 mm in diameter, which makes sliding contact with the sheet. The stylus tool proceeds in the counterclockwise direction and steps down 2 mm after each full revolution. The assumed tool speed is 4000 mm/min. Mass scaling of 4000 is applied to increase the stable time step size and to speed up the simulation. The control point (or node) spacing is 2 mm.  The formed part is shown in Figure 3. The smooth results from the isogeometric solver are given in Figures 5 and 7. The simulation accounts for the unclamping step, i.e. the release of the clamps after the forming process has finished. After the unclamping step, the springback effect bends the edges of the sheet slightly upwards. The edges perpendicular to the longer sides of the die become more curved than the other edges. The springback predicted by the simulation matches very well with the experimental results. Furthermore, Figure 4 shows the actual thickness determined by 3D scanning both surfaces and aligning them to produce the thickness contours. The simulated thinning shown in Figure 7 is in a good agreement with the experimental results (same scales). The results from the bilinearly interpolated BWC shell element that is based on the plane stress assumption are given as a reference in Figures 6  and 8. The solution is rough in comparison with IGA as the displacement field is now only C 0 continuous. The thickness plot is even rougher as the thickness is only available at the integration points. Nevertheless, the results are still good in comparison with the experimental results.  Unfortunately, both element types overestimate the averaged tool forces by ∼50%. However, this is a common problem which has been associated with many FEM formulations in the past [18]. Furthermore, the isogeometric shell formulation overestimates these forces slightly more than the traditional BWC element. The overestimation is suspected to originate from the material model (mostly due to lack of a damage model), but also the contact formulation and element locking may play a role. Further investigation is required.
In the isogeometric analysis, the boundary Bézier elements pose stricter requirements for the stable time step size compared to the interior elements as shown by Adam et al. [6]. One way to overcome this limitation is to increase the boundary element size, which, indeed, may be acceptable for many ISF simulations.
However, it is noticed that preserving the uniform control point spacing, but simply clamping all the edges also effectively prevents instabilities originating from the boundaries.
The implemented explicit isogeometric shell formulation is found visually stable with up to ∼50% larger maximum time step size compared to the traditional BWC shell. This increase in the time step size makes close to no difference to the final results. However, it is found that the large time step size cannot be preserved if all edges are not clamped (with uniform control point spacing). Therefore, the time step size must be reduced right before the unclamping step at the end of the simulation.

Conclusions
The isogeometric shell formulation proposed by Hokkanen et al. [11] is applied to a challenging ISF problem. A comprehensive anisotropic material model is implemented to take into account the sheet rolling direction. As the chosen problem requires a die, double-sided contact must be considered. The automatically adjusting penalty stiffness is proposed to model the contact between the stylus tool and the sheet, whereas the die contact algorithm uses traditional constant penalty stiffness.
Qualitatively good results are obtained by the isogeometric shell as well as the traditional BWC shell formulation. The displacements, springback effects, and thinning agree very well with the experimental results. However, the forming forces are significantly overestimated by both shell types which requires further investigation. The isogeometric shell formulation is found very competitive even against low-order reduced integration shell elements in speed due to the increased accuracy per degree of freedom and larger stable time step size.