A high‐order shell finite element for the large deformation analysis of soft material structures

This work proposes a higher‐order unified shell finite element for the analysis of cylinders made of compressible and nearly incompressible hyperelastic materials. The nonlinear governing equations are derived employing the Carrera unified formulation (CUF), thanks to which it is possible to build shell elements with the capability to capture three‐dimensional (3D) transverse and out‐of‐plane effects. The material and geometric nonlinearities are expressed in an orthogonal curvilinear reference system and the coupled formulation of hyperelastic constitutive law is considered. The principle of virtual work and a total Lagrangian approach is used to derive the nonlinear governing equations, which are solved by a Newton–Raphson scheme. The numerical investigations deal with a curved arch and both thick and thin cylinders subjected to line and point loadings. The obtained results are validated by comparing them with those from the literature. They demonstrate the reliability of the proposed method to analyze compressible and incompressible hyperelastic shell structures.

related numerical implementation of shells was a challenging task for scientists.Large deformations most likely appear in thin shells, and a geometrical linear approach is not able to describe such phenomena.Therefore, thin shells are often described using a geometrical nonlinear approach.Many studies have applied the St. Venant-Kirchhoff constitutive model, which leads to reliable results when shells are subjected to large displacements, but not in the case of large strains.Thus, for thin shells that undergo both large displacements and strains, the hyperelastic constitutive model must be included to capture their static response.
In recent decades, a remarkable effort on soft and hyperelastic shell structures was dedicated by scientists and engineers, due to their potential application to biomaterials.A few works are mentioned hereafter.Basar and Ding, 6 starting from a quadratic approximation over the thickness of thin shells and neglecting transverse shear effects, developed a nonlinear model for incompressible hyperelastic material.Campello et al. 15 introduce the Rodrigues rotation vector to carry out the analysis of nonlinear hyperelastic shells using a neo-Hookean material model.Basar and Itskov 16 propose a numerical implementation of the Ogden material model via a variational procedure.The large strain and finite rotation fields are analyzed by introducing a six parametric shells kinematic.The absolute nodal coordinate formulation was adopted by Luo et al. 17 to carry out the static and dynamic nonlinear analysis of thin shells.The proposed approach is based on the Kirchhoff-Love hypotheses.Spherical and circular cylindrical shells were investigated by Song and Dai, 18 using a finite-strain shell theory.They employed a high-order expansion from the 3D equations for compressible hyperelastic materials.An extensible director kinematic is used by Betsch et al. 19 and Li et al. 20 for the analysis of compressible and incompressible hyperelastic shells, respectively.A nonlinear Lagrangian formulation was employed.
To the best of the authors' knowledge, much more attention was focused on thin shells made of hyperelastic materials, rather than their thick counterparts.The latter would require the adoption of more refined theories, which could eventually describe the through-the-thickness (both normal and shear) strains employing a refined approximation of the displacement field.More in general, when the thickness deformation plays a crucial role in the mechanical behavior of a structure, the adoption of higher-order theories is demanding.This is the case of thick shells and, in particular, of biological and rubber-soft materials, as outlined in many review articles ( [21][22][23][24] ).Moreover, for biological applications, multilayered structures are often adopted, and they represent typical cases when an accurate description of the transverse components is mandatory.An important contribution for the analysis of shells with higher-order models is given by Arbind and Reddy 25 which employed a general polynomial interpolation of the displacement field to develop of a general refined shell hyperelastic model.Moreover, Amabili et al. 26 developed a geometrically nonlinear model for incompressible cylindrical hyperelastic structures.By employing a 9-parameter theory, a higher-order description of the in-plane and out-of-plane strain components was carried out.
It can be said that the main difficulties when dealing with hyperelastic shells are due to the complexity of the formulation for a hyperelastic material and the derivation of the nonlinear governing equations for shells.In this context, the Carrera unified formulation (CUF) 27 represents a good framework.Due to its hierarchical nature, it is possible to build higher-order shell models by means of the so-called fundamental nuclei which, opportunely expanded, compose the stiffness matrix.The form of those does not depend on the theory order, which can be chosen by the user as input for the analysis.CUF has already been applied for the analysis of shell structures in the geometrical nonlinear field (see References 28 and 29) and for hyperelastic structures (see Reference 30).In this article, CUF is further extended to deal with hyperelastic thin and thick shell structures.
This work is organized as follows.An introduction to the geometric relation of shell structures is given in Section 2 in the curvilinear coordinate system.The hyperelastic constitutive law in the coupled form (following the Holzapfel 31 terminology) is described in Section 3. Section 4 describes the kinematics of the proposed nonlinear shell formulation and the derivation of the governing equations introducing the finite element method (FEM).The numerical results, proposed in Section 5, involve the analysis of cylindrical thick and thin shells made of compressible and incompressible hyperelastic materials.The article concludes by drawing the main conclusions.Also, appendix sections are provided, and they give the components of the linear and nonlinear differential operators and the material Jacobian tensor.

Reference configuration of a shell
Let us consider an arbitrarily curved shell as shown in Figure 1.The structure is considered to have a constant thickness h.The shell lays within a curvilinear coordinate system  1 ,  2 ,  3 , where  3 lays along n, which is the unit vector normal to Geometry of an arbitrary shell in an orthogonal curvilinear coordinate system.
the surface that is defined by the two in-plane parametric curves  1 and  2 .It is clear that  3 ∈ [−h∕2; h∕2].The relation between the position vectors X (point P (0) ) and  (point P (0) 0 ) can be expressed as follows: Starting from the position vector X, the covariant base vectors g i and it dual counterpart base vectors g i can be derived as follows: where (i = 1, 2, 3).Moreover, the partial derivative within a curvlinear coordinate system of a generic vector A reads as follows: where the spatial covariant derivative is conventionally expressed by the semicolon ";", and it is defined as follows: Equation ( 4) introduces the Christoffel symbols, whose components can be expressed in the following: If the directions  1 and  2 are considered to lay along the curvatures of the shell structure, then the coordinate system results to be orthogonal, that is the considered case in the present study.Thus, g i and g i result to be orthogonal, and their modules are defined as follows: with H i being the Lamé parameters of an arbitrary parallel surface of the shell.They are defined as follows: where R 1 and R 2 are the principal curvature radii, whereas the mid-surface Lamé parameters are defined by A and B. In this research, cylindrical shells are examined, and, as a result, constant curvatures are presumed.Thus, both R 1 and R 2 are constant.

Motion of a shell in the orthogonal curvilinear system
Consider the motion of a shell between two configurations, as depicted in Figure 2, where the position vector of the point in the final configuration is given by x.The displacement vector u can be expressed in the total Lagrangian description as follows: Therefore, the deformation gradient tensor F is defined as follows: where the ⊗ and ∇ symbols stand for the dyad operators and the gradient operator with respect to the reference configuration, respectively.A set of orthogonal normalization vectors e i is here introduced to simplify the analysis of physical problems, since g i and g i are not necessarily unit vectors in an orthogonal curvilinear coordinate system (see Equation ( 6)).
Then, the displacement vector can be re-written, in terms of e i , as follows: F I G U R E 2 Motion of an arbitrary shell in an orthogonal curvilinear coordinate system.
where ũi express the physical components.Introducing Equations ( 8) and ( 10) into (11), one has Thus, considering Equation (3), Equation ( 9) can be expressed as follows: In this work, we consider  1 and  2 as the arc-length coordinates, and in this work they are referred to as  and , respectively (z for  3 ), and consequently, A = B = 1.The 3 × 3 deformation gradient tensor F of Equation ( 9) is given by The Green-Lagrange strain tensor E can be expressed as follows: where the transpose is indicated by the superscript "T."The linear E l and nonlinear E nl components of the strain E in Equation ( 9) can be written as follows: and Subsequently, Equations ( 16) and ( 17) can be rearranged in a matrix form as follows: where u and E can be defined as The 6 × 3 linear and nonlinear differential operators b l and b nl in Equation ( 18) are given in Appendix A.

HYPERELASTIC CONSTITUTIVE LAW
In the present work, isotropic hyperelastic behavior is described adopting a neo-Hookean strain-energy function Ψ.In general, the strain energy function Ψ can be written in terms of the invariants (I 1 , I 2 , I 3 ) of the C = F T F right Cauchy-Green tensor as follows: where the invariants are defined as follows: where tr(⋅) and det(⋅) represent the trace and the determinant operators.
The constitutive law is expressed as a differential relation between Ψ and the second Piola-Kirchhoff (PK2) stress tensor as follows: where Introducing Equation (23) into Equation ( 22), one has the final expression of the PK2 tensor, The material Jacobian tensor is given in Appendix B.

Unified shell finite element
In this work, the 2D shell finite element is built recalling the CUF.According to CUF, the 3D displacement field of Equation ( 19) can be derived as a set of in-plane functions over the  and  directions and thickness functions over the thickness domain z.,  and z are the physical coordinates.Specifically, u can be written as follows: where  = 0, 1, … , N, with N being the order of expansion functions F  .The repeated subscript  indicates the summation.In this article, cubic and quadratic expansion functions, based on three-point and four-point Lagrange polynomials, are employed.Following the notation introduced by Carrera, 32 they are recalled as Lagrange Displacement LDo, which stands for Lagrange expansion with Displacement unknowns of order o (2 and 3 for quadratic and cubic expansion, respectively).For the sake of completeness, the case of a LD2, namely a quadratic expansion, is reported hereafter where where s vary from 1 to -1, whereas s 1 and s 3 correspond to the position of the Lagrange points in the natural coordinate.
The relation between the natural and physical coordinates can be found in many books, see Reference 27.On the other hand, for the - domain discretization, FEM is recalled.After discretizing the shell mid-plane in planar finite elements, the generalized displacement vector u  (, ) can be expressed in the following: where i = 1, 2, … , p + 1, with p representing the order of the shape functions N i .The vector of the FE nodal parameters q sj is written as Note that the finite element interpolation is expressed by using global coordinates in Equation ( 28).However, stiffness terms as detailed in Section 4.2 are derived by integrating both the shape functions N i and the theory approximation expansions F  in the natural domain, as in classical FEM.In this work, nine-node bi-quadratic FEs (Q9) are adopted as shape functions.Figure 3 depicts the discretization of an arbitrary shell structure, where the expansion function F  and the shape function N i are highlighted in red and blue, respectively.
Note that, introducing CUF Equation ( 27) and FEM (28) into Equation ( 18), the resultant strain vector can be expressed as follows: where B i l and B i nl are the linear and nonlinear algebraic matrices.Interested readers can find more details about B i l and B i nl in Reference 28.

Nonlinear governing equations
In the present work, the derivation of the nonlinear governing equation is conducted recalling the PVDs.It can be written as follows: where L int and L ext are the virtual variation of the internal strain energy and the external loads, respectively, with  denoting the variation.The virtual variation of the internal work can be expressed as follows: where V is the volume of the body, and Introducing Equation (33) into Equation ( 32), one has where F i int is the 3 × 1 fundamental nucleus of the internal forces vector: The expression of the external load vector is computed from the definition of the virtual variation of the work made by external forces.Omitting some mathematical derivations, we have Finally, introducing Equations ( 34) and ( 36) into (31), the nonlinear equation holds Equation ( 37) represents a geometrically nonlinear systems, and it is typically computed adopting a linearization technique.In this article, the employed scheme is the Newton-Raphson method, according to which, the nonlinear governing equations are expressed as follows: where the residual nodal forces vector is expressed in  res .One can use a known (q, p) solution to linearize Equation ( 38) by expanding  res in Taylor's series.Therefore, where  res q = K T represents the tangent stiffness matrix.The external load is assumed to change directly with the vector of the reference loadings p ref with a variation rate expressed by , defined as the load-scaling parameter, i.e. p =  p ref .
Since  is a variable, an additional constraint equation is required and this is given by a relation constraining both q and ).Finally, one has In this work, a path-following technique is adopted as the constraint equation.In particular, an arc-length method is utilized in this work, see Criesfield 33,34 and Carrera. 35It is important to underline that K T is derived from the linearization of the equilibrium equations. 36This corresponds to the linearization of the virtual variation of the work made by internal forces in the case of conservative cases, as follows where K ijs T represents the 3 × 3 fundamental nucleus, that is, the basic building block for the formulation of the total tangent stiffness matrix.

NUMERICAL RESULTS
The numerical results concern three examples taken from the literature.Each of them is analyzed considering a neo-Hookean strain energy function.In particular, the study cases involve a curved arch and both thick and thin cylinders subjected to line and point load.For the thin cylinder subjected to a point load, both compressible and incompressible material configurations are considered, and the obtained results are compared to those from literature, when available.

Curved arch
As a first example, the behavior of a compressible soft curved arch is investigated.Figure 4 shows the geometrical characteristics and boundary condition of the structure, with r i and r e equal to 9 and 10 mm, respectively, and a = 1 mm.Then, the structure can be considered as thick (r e /a = 10), and can be modeled using the proposed shell finite element thanks to its higher-order capability.The structure is clamped and subjected to an external pressure oriented by 45 • with respect to the normal of the surface.The study case is taken from Hassani et al., 37 and the neo-Hookean material model is described in the following strain energy function: where  and k are the shear and bulk modulus, respectively.Their values are  = 80.194 N/mm 2 and k = 120.291N/mm 2 , so that the material is compressible.
A preliminary convergence analysis is carried out, in order to establish a reliable mathematical model.Figure 5 shows the results.Clearly, the 120 Q9 and 2 LD3 approximation can be considered a converged mathematical model.Its total number of Degrees of Freedom (DOFs) is 11907.The nonlinear static equilibrium curve is reported in Figure 6, along with some of the deformed configurations.A good match between the proposed model solution and the reference results is guaranteed.

Cylindrical shell subjected to line loading
The hyperelastic cylindrical structure shown in Figure 7   to 2 and 0.2 cm, respectively.A neo-Hookean material model is adopted with the same strain energy function Φ adopted in Reference17 with shear modulus  = 6000 kN∕cm 2 and bulk modulus k = 280,000 kN∕cm 2 .For symmetry reasons, only one quarter of the structure is considered (depicted in blue in 7).A 4 × 10 Q9 mesh is adopted, along with one LD2 on the thickness direction, resulting in a total of 1701 DOFs.The thick cylinder was analyzed as the first case, and the results are compared to reference ones from Luo et al.17 The nonlinear equilibrium curve is shown in Figure 8, proving a perfect match with the reference solution.Moreover, the thin cylinder is considered.The available reference solution is from Büchter et al., 38 where a single equilibrium state in the force vs displacement graph is provided (see Figure 9).The solution obtained with the present model is reported in the same graph, and a consistent match with literature results is demonstrated.In the same figure, the deformed configuration at maximum vertical displacement equals 16 cm is provided.

Cylindrical shell subjected to pinching loading
As the final case, a cylindrical shell with rigid end diaphragms is subjected to pinching point load P, as shown in Figure 10.This example was studied in References 39 and 40 by using different approaches, whereas the results chosen

F I G U R E 11
Convergence analysis of the cylinder under pinching point load.For each model one LD2 is adopted over the thickness direction.
here as reference are from Luo et al. 17 The radius, length, and thickness of the cylindrical shell are r = 100 mm, L = 200 mm, and t = 1 mm, respectively.The boundary condition of rigid diaphragms is enforced on the two edges of the cylindrical shell.This means that the displacements of the edges of the cylindrical shell are constrained except along the longitudinal direction.Thanks to the symmetrical property of the problem, only one octant of the shell is studied.A convergence analysis is carried out as a preliminary investigation to assess the mathematical model, and the results are shown in Figure 11.Clearly, the 40 × 40 Q9 mathematical model can be considered as converged, and is used for the subsequent investigations.The number of DOFs is 59049.Both compressible and nearly incompressible neo-Hookean material models are considered here, with Poisson's ratios  = 0.3 and  = 0.499, respectively.
The obtained results in terms of nonlinear equilibrium curve are shown in Figure 12 and, in the compressible case, they match well with the reference solution from Luo et al. 17 Moreover, Figure 13 depicts some of the deformed shapes of the cylindrical shell at the nonlinear static load steps of maximum vertical displacement equals to 57, 76, 86, and 96 mm, respectively.18) are reported here. and shell model.Shape function N j in blue and expansion function F s in red.

F G U R E 4 37 F I G U R E 5 37 F I G U R E 7 F I G U R E 8
is considered as the second analysis case.The geometric properties are r = 9 cm and L = 15 cm.Both thick and thin configurations are considered, with thicknesses t and t 2 equal Geometric properties of the curved arch.Study case taken from Hassani et al.Convergence analysis of the arch subjected to an external pressure.F I G U R E 6 Nonlinear static equilibrium curve of the arch subjected to an external pressure.Some deformed configurations are reported.Reference results from Hassani et al.Study case of the cylinder under line load.Nonlinear analysis of the cylinder under line load.Thickness equals to 2 cm.Reference result from Luo et al.17

F I G U R E 9
Nonlinear analysis of the cylinder under line load.Thickness equals to 0.2 cm.Reference result from Büchter et al.38

F I G U R E 10
Study case of the cylinder under pinching point load.

F I G U R E 12
Nonlinear analysis of the cylinder under pinching point load.Reference solution from Luo et al.17 for compressible case.

F I G U R E 13
Deformed shapes of the cylinder under pinching point load.

6 CONCLUSIONS
How to cite this article: Pagani A, Augello R, Carrera E. A high-order shell finite element for the large deformation analysis of soft material structures.Int J Numer Methods Eng.2024;125(7):e7417.doi: 10.1002/nme.7417APPENDIXA. LINEAR AND NONLINEAR DIFFERENTIAL OPERATORSThe 6 × 3 linear and nonlinear differential operators b l and b nl in Equation (