Isogeometric hyperelastic shell analysis with out-of-plane deformation mapping

We derive a hyperelastic shell formulation based on the Kirchhoff–Love shell theory and isogeometric discretization, where we take into account the out-of-plane deformation mapping. Accounting for that mapping affects the curvature term. It also affects the accuracy in calculating the deformed-configuration out-of-plane position, and consequently the nonlinear response of the material. In fluid–structure interaction analysis, when the fluid is inside a shell structure, the shell midsurface is what it would know. We also propose, as an alternative, shifting the “midsurface” location in the shell analysis to the inner surface, which is the surface that the fluid should really see. Furthermore, in performing the integrations over the undeformed configuration, we take into account the curvature effects, and consequently integration volume does not change as we shift the “midsurface” location. We present test computations with pressurized cylindrical and spherical shells, with Neo-Hookean and Fung’s models, for the compressible- and incompressible-material cases, and for two different locations of the “midsurface.” We also present test computation with a pressurized Y-shaped tube, intended to be a simplified artery model and serving as an example of cases with somewhat more complex geometry.

In this article, we start with the formulation from [4] and derive, based on the Kirchhoff-Love shell theory and isogeometric discretization, a hyperelastic shell formulation that takes into account the out-of-plane deformation mapping. Accounting for that mapping affects the curvature term. It also affects the accuracy in calculating the deformedconfiguration out-of-plane position, and consequently the nonlinear response of the material. We are extending the range of applicability of Kirchhoff-Love shell theory to the situations where the Kirchhoff-Love shell kinematics is still valid, yet the thickness or the curvature change is significant enough to make a difference in the response. Fung's model has different versions. In the version used in [12], the first invariant of the Cauchy-Green deformation tensor appears in a squared form. In the version we use in this article, it appears without being squared, and this version has been used in a number of arterial FSI computations [24][25][26][27][28][29][30][31] with the continuum model.
In FSI analysis, when the fluid is inside a shell structure, the shell midsurface is what it would know. That would be physically wrong, especially when the thickness is significant, because the inner surface is the one that the fluid should really see. In this article we also propose, as an alternative, shifting the "midsurface" location in the shell analysis to the inner surface. Furthermore, in performing the integrations over the undeformed configuration, we take into account the curvature effects, and consequently integration volume does not change as we shift the "midsurface" location.
To evaluate the performance of the shell formulation presented, we do test computations with pressurized cylindrical and spherical shells, with Neo-Hookean and Fung's models, for the compressible-and incompressible-material cases, and for two different locations of the "midsurface." We compare the results to near-analytical reference solutions. We also do test computation with a pressurized Y-shaped tube, intended to be a simplified artery model. This serves as an example of cases with somewhat more complex geometry.
In Sect. 2 we provide the governing equations. The hyperelastic shell model is presented in Sect. 3. Test computations with the cylindrical and spherical geometries and Y-shaped tube are presented in Sect. 4, and the concluding remarks in Sect. 5. In the Appendix, we provide some derivations used in Sect. 3, and the constitutive laws.

Governing equations
Let Ω t ⊂ R n sd be the spatial domain with boundary Γ t at time t ∈ (0, T ), where n sd is the number of space dimensions. The subscript t indicates the time-dependence of the domain. The equations governing the structural mechanics are then written, on Ω t and ∀t ∈ (0, T ), as where ρ, y, f and σ σ σ are the density, displacement, body force and Cauchy stress tensor. The essential and natural boundary conditions are represented as y = g on (Γ t ) g and n · σ σ σ = h on (Γ t ) h , where n is the unit normal vector, and g and h are given functions. The Cauchy stress tensor can be obtained from where F and J are the deformation gradient tensor and its determinant, and S is the second Piola-Kirchhoff stress tensor. It is obtained from the strain-energy density function ϕ as follows: where E is the the Green-Lagrange strain tensor: C is the Cauchy-Green deformation tensor: and I is the identity tensor. From Eqs. (3) and (4),

Hyperelastic shell model
We split the domain as where Γ t represents the midsurface of the domain, which is parametrized by n pd = n sd − 1, where n pd is the number of parametric dimensions. With the position x ∈ Γ t , we define a natural coordinate system: where α = 1, . . . , n pd , and the third direction is based on The components of the metric tensor are and this is known as the first fundamental form. Similarly, we define the components of the metric tensor for the contravariant basis vectors as and obtain the tensor components and contravariant basis vectors from and We define and with that, components of the covariant curvature tensor are and this is known as the second fundamental form.
A position x ∈ Ω t is represented as where −1 ≤ ξ α ≤ 1 and ξ 3 ∈ (h th ) t . The basis vectors are represented as = g α + n ,α ξ 3 (20) See Appendix A.1 for the lines between Eqs. (20) and (21). Because g α and g α are on parallel planes (from the Kirchhoff-Love shell theory), With that, the metric tensor components in 3D space are

Remark 1
The quadratic term may be omitted. However, if the metric tensor is obtained from the basis vectors, the term will automatically be included.
We now provide similar definitions and derivations for the undeformed configuration. We start with the basis vectors: and A position X ∈ Ω 0 is expressed as where −1 ≤ ξ α 0 ≤ 1 and ξ 3 0 ∈ (h th ) 0 . The basis vectors are represented as The metric tensor components in 3D space are and B αβ is the second fundamental form for the midsurface of the undeformed configuration. On the midsurface the parametric coordinates indicate the same material points, and therefore, ξ α = ξ α 0 . In the third direction, however, because of the normalization, the coordinates may not be the same. The relationship is where λ 3 is the stretch in the third direction.

Kinematics
We obtain F from the following relationship: This means that Then we can write F as and J as From Eq. (5), we can write C as and the determinant of C gives the square of J : From Eq. (4), we can write E as The covariant components of the in-plane strain tensor are We write ξ 3 ξ 3 0 as ξ 3 ξ 3 0 =λ 3 ξ 3 0 ξ 3 0 . From the Taylor expansion ofλ 3 around ξ 3 0 = 0, we obtain We note that λ 3 is the stretch at ξ 3 0 = 0, which isλ 3 (0). With that,

Constitutive equations
The total differential of the second Piola-Kirchhoff stress tensor is where I , J , K , L, M, N = 1, . . . , n sd . From Eq. (4), the following expression can be used: For shells, because dE 3α = dE α3 = 0. From Eq. (65), we can write and from Eqs. (4) and (66), we can write From the plane stress condition S 33 = 0, dS 33 = 0, and consequently which makes and therefore we introducê In computing C 33 , we have different methods for incompressible and compressible materials. In the case of incompressible material, from Eq. (49) we can write C 33 = λ 2 3 . Because J = 1, and therefore In the case of compressible material, as can be found in [32], C 33 can be calculated by Newton-Raphson iterations that would make S 33 = 0. Because C γ δ does not change during the iterations, the iteration increment is From Eq. (67) and remembering that dC γ δ = 0 during the iterations, The update takes place as where superscript i is the iteration counter, and as the initial guess we have the following three options: The option given by Eq. (78) comes from the constitutive law for zero bulk modulus. To preclude C 33 being negative, we introduce an alternative update method based on the logarithm of C 33 :

Variational formulation
The variation of the in-plane components of the Green-Lagrange tensor is The variation of ξ 3 can be dropped (see Appendix B), and we obtain With that, which means Remark 2 Evaluation of S requires a material point correspondence in the third direction. We take that into account by integrating Eq. (40) with the 4th order Runge-Kutta method, and λ 3 can be obtained from the constitutive law given in Sect. 4.1. Figure 1 illustrates the deformation mapping. In general, stretch at a convex side is less than the stretch at the concave side, which results in a nonuniform λ 3 .
Now we derive what we need: where and the variation of the normal vector (see Appendix A.2) is

Linearization for the Newton-Raphson iterations
The linearization for the Newton-Raphson iterations is doneas The variation with subscript a is associated with the variational formulation, and the variation with subscript b is associated with the iteration linearization. Again, the variation of ξ 3 is dropped. In this part too, we derive what we need: Here, we used and the proof for this can be found in Appendix C.

Test problems
We test the formulation given in Sect. 3 by using pressurized cylindrical and spherical shells, with Neo-Hookean and Fung's models, for the compressible-and incompressiblematerial cases, and for two different locations of the "midsur- face." For the compressible-material cases, we include tests with the continuum model. We compare the results to nearanalytical reference solutions. We also do test computation with a pressurized Y-shaped tube.

Constitutive models
We test two constitutive models: neo-Hookean and Fung's. The elastic-energy density functions are where μ is the shear modulus, and D 1 and D 2 are the coefficients of the Fung's model. For incompressible material, we use where p is the pressure, which can be eliminated by the plane stress condition. For compressible material, we use where and κ is the bulk modulus.

Test computations
The pressure, applied at r = r p , is normalized by the shear modulus (at the zero-stress state): for the Fung's model, and we use D 2 = 8.365. We determine the bulk modulus from the Poisson's ratio ν as follows: for the neo-Hookean model, and for the Fung's model. How to deal with pressure acting on the inner surface is not easy because the midsurface is the geometry we are using in the computation. Here we propose two ways. In the first one, "midsurface model," the pressure is applied on the midsurface of the current configuration. In the second one,"inner-surface model," the structure "midsurface" is moved to the inner surface and the pressure is applied there.

Remark 3
In applying the pressure, the midsurface model is physically wrong, especially when the thickness is significant. The inner-surface model will have larger absolute value for ξ 3 , which would lead to larger discretization errors.
In the test cases, we use the inner and outer radii R I and R O , and the thickness H = R O − R I . The condition used here is H 2R I = 0.1, which is slightly thinner than most arteries. To have a reference solution to compare the results to, we provide in Appendix D the second Piola-Kirchhoff stress tensor expressed in terms of the principal stretches. The results are compared by inspecting pressure as a function of stretch. The stretch is λ 1 ≡ r p R for the midsurface model, where r p = r , and λ 1 ≡ r p R I for the inner-surface model, where r p = r I . In the computations, we increase the pressure gradually in obtaining the solution and calculate the stretch. In obtaining the reference solutions, we use numerical integrations, which are explained in the following subsections.

Pressurized cylinder
We use orthogonal basis vectors: the first basis vector is in the radial direction, the second one is in the cylinder axis direction, and the third one is normal to the surface. The force equilibrium gives the following relationship: Because the cylinder height does not change, λ 2 = 1, and we obtain See Appendix D.1 for S 11 . Figures 2 and 3 show the reference solutions for the neo-Hookean and Fung's models.

Remark 4
For the Fung's model, to show convergence to incompressible-material response with increasing Poisson's ratio, for the midsurface model we add Fig. 4, with two more Poisson's ratios beyond those in Fig. 3.
We compute, in 2D, with uniform, periodic cubic Bsplines with 8 elements. For comparison purposes, we also compute with the continuum model, using 128 uniform, periodic cubic B-spline elements in the circumferential direction, and 1 element in the radial direction.

Remark 5
To compare the solutions from the shell models with and without the out-of-plane deformation mapping, we use near-analytical solutions to represent the shellmodel solutions (see Appendix E). We do this for the incompressible-material case, with the midsurface model. To have some sense of scale for the difference between the solutions, we include in the comparison the solution we get when we use 1-point ξ 3 integration in the model without the out-ofplane deformation mapping. Figure 11 shows the comparison for the Fung's model. For the neo-Hookean model, there is no

Pressurized sphere
We use orthogonal basis vectors: the first two vectors are on the surface, and the third vector is normal to the surface. The force equilibrium gives the following relationship: Because of the symmetry between the two basis vector directions on the surface, λ 1 = λ 2 , and we obtain We compute, in 3D, with a cubic T-spline mesh, which consists of 296 control points and 534 Bézier elements (see Fig. 15). For comparison purposes, we also compute with the continuum model, which is extruded in the thickness direction with 1 element.

Remark 7
The number of elements used in the integration is the number of Bézier elements, which is 534 in this case. The mesh was generated by a commercial software, Rhinoceros with the T-splines plug-in. It actually has, in the finite element sense, 294 elements.

Remark 8
In the same way described in Remark 5 for the pressurized cylinder, we compare the solutions from the shell models with and without the out-of-plane deformation mapping. Figure 22 shows the comparison for the Fung's model. For the neo-Hookean model there is no visible difference between the solutions. We again note that the curvature changes are very small in the test problem.

Pressurized Y-shaped tube
The undeformed configuration of the tube is shown in Fig. 23. The end diameters of the tube are 20, 14 and 10 mm.   Figure 25 shows the maximum principal curvature for the undeformed configuration. We use the incompressiblematerial Fung's model with D 1 = 2.6447×10 3 Pa and D 2 = 8.365. The "midsurface" location in the shell analysis is the inner surface, and the pressure applied is 12.3×10 3 Pa. The thickness distribution for the undeformed configuration is shown in Fig. 26. This smooth thickness distribution is outcome of solving the Laplace's equation over the inner surface of the tube, with Dirichlet boundary conditions at the tube ends, where the value specified is 0.146 times the

Concluding remarks
We have presented a hyperelastic shell formulation based on the Kirchhoff-Love shell theory and isogeometric discretization. The formulation takes into account the out-of-plane deformation mapping. Accounting for that mapping affects the curvature term. It also affects the accuracy in calculating the deformed-configuration out-of-plane position, and consequently the nonlinear response of the material. In FSI analysis, when the fluid is inside a shell structure, the shell midsurface is what it would know. That would be physically wrong, especially when the thickness is significant, because the inner surface is the one that the fluid should really see. For that reason, in this article we also proposed, as an alternative, shifting the "midsurface" location in the shell analysis to the inner surface. The way we perform the integrations over the undeformed configuration takes into account the curvature effects, and consequently integration volume does not change if we shift the "midsurface" location. We presented test computations with pressurized cylindrical and spherical shells, with Neo-Hookean and Fung's models, for the compressible-and incompressible-material cases, and for two different locations of the "midsurface." We compared the results to near-analytical reference solutions, and in all cases we see a good match. We also presented test computation with a pressurized Y-shaped tube, intended to be a simplified artery model and serving as an example of cases with somewhat more complex geometry.

Appendix B: Variation of 3 is a second-order term
Taking the variation of both sides of Eq. (58), we obtain For a representative value of the variation, we take the average over the thickness: Thus, the variation of ξ 3 is a second-order term.