Stretch-Based Hyperelastic Material Formulations for Isogeometric Kirchhoff-Love Shells

A stretch-based hyperelastic material formulation for isogeometric Kirchhoff-Love shells is presented, particularly useful when experimental material data fitting is involved to capture the model parameters, for instance using the Ogden material model. Complementing the existing invariant-based formulation, the stretch-based stress and material tensors are expressed in spectral form and transformed to the covariant curvilinear bases for consistency with the variational formulation. For both compressible and incompressible material formulations, analytical and numerical benchmarks show excellent results. In particular, the response of a collapsing conical shell a complex collapsing mechanism was revealed when using the arc-length method.


Introduction
To model thin shells or membranes that undergo large strains, nonlinear material relations are required, rather than the commonly used Saint Venant-Kirchhoff model. Nonlinear material models (i.e. constitutive relations) are often called (finite) hyperelastic models. In general, one distinguishes between phenomenological continuum mechanics approaches and molecular statistical approaches [1], which can both be described in terms of invariants or principal stretches of the deformation. Phenomenological material models include the well-known neo-Hookean, Mooney-Rivlin and Ogden material models, and statistical ones are for example the Wang-Guth and Arruda-Boyce models [1]. Numerical implementation of invariant-based material models require a less extensive computational procedure than stretch-based formulations, since principal stretches and their directions do not have to be computed. On the other hand, stretch-based descriptions can be very useful when experimental data fitting is required to capture the model The Einstein summation convention is adopted to represent tensor operations. In this notation, the trace and determinant of a tensor are defined as follows for tensor A = A i j a i ⊗ a j as [12,24,25] tr A = A i j a i j and det{ A} = A i j Where A i j denotes the determinant of the matrix A. The inverse of a tensor A is denoted by A −1 orĀ. The derivative of the inverse and the determinant of a tensor, with respect to one of its components become:

Shell Coordinate System
The Kirchhoff-Love shell element formulation is based on the Kirchhoff Hypothesis, that is, the cross-section does not shear and orthogonal vectors in the undeformed configuration remain orthogonal after deformation. As a consequence, any point in the shell can be represented by a point on the mid-surface and a contribution in normal direction: with the shell mid-surface by r(θ 1 , θ 2 ) and the unit normal direction n(θ 3 ) for the deformed configuration x(θ 1 , θ 2 , θ 3 ). For the undeformed configurationx, the same relation holds with all quantities decorated with a·. The parametrization utilizes surface coordinates θ α and the through-thickness coordinate θ 3 . Derivatives with respect to these coordinates are denoted by (·) ,i = ∂(·)/∂θ i .
The covariant basis of the mid-surface is represented by a i a α = ∂r ∂θ α , a 3 = n = and the first fundamental form is a αβ = a α · a β . The curvature tensor b αβ is represented by the second fundamental form of surfaces, which can be obtained using the Hessian of the surface a α,β or the derivative of the normal vector n ,α b αβ = n · a a,b = −n ,β · a α .
The derivative of the normal vector can be obtained by Weingarten's formula n ,α = −b β α a β with b β α = a αγ b γβ as the mixed curvature tensor [27]. Taking the derivative of Eq. (3), the covariant basis of the shell coordinate system x can be formulated as follows: The metric coefficients are constructed by taking the inner-product of these basis vectors, i.e.

Remark 1.
In the isogeometric Kirchhoff-Love shell formulations [10,12], the last term in Eq. (7) is neglected because of the thin shell assumption, meaning (θ 3 ) 2 takes small values. However, the co-and contravariant basis vectors (g α and g α , respectively) are used in the mapping of the stretch-based material matrix onto the contravariant undeformed basis (Section 5.3). To enable an accurate comparison of the invariant-based and stretch-based formulation, we do not neglect the O((θ 3 ) 2 ) term, contrary to previous works [12,15].

Shell Kinematics
The deformed and undeformed configurations (x andx, respectively) are related to each other by the mid-plane deformation vector u by r =r + u and n = n(r + u). However, in both the invariant-based and stretch-based forms that are described in this paper, the deformations are defined using the undeformed and deformed geometries. In continuum mechanics, we define the deformation gradient F and the deformation tensor C are defined as [12,26]: Note that the deformation tensor is defined in the contravariant undeformed basis g i ⊗ g j . For shells, it is known that g α3 = g 3α = 0, hence this implies C α3 = C 3α = 0. Since g 33 = 1, which implies C 33 to be unity and meaning that the thickness remains constant under deformation. In hyperelastic Kirchhoff-Love shell formulations, the contribution of C 33 is therefore incorporated by static condensation (incompressible materials) or iteratively (compressible materials). Therefore, we denote the deformation tensor C and its inverseC as denoted as: From Eqs. (10) and (11), it can be observed that the thickness-contribution (index 3) is decoupled from the in-plane contributions (Greek indices α, β). This decoupling is a consequence of the Kirchhoff Hypothesis and therefore is only valid for Kirchhoff-Love shells. As a consequence, 4 using the definitionC = g αβg α ⊗g β , the trace and determinant of C can be simplified accordingly [24,25]: where J denotes the Jacobian determinant and J 0 is its in-plane counterpart. Furthermore, the tensor invariants of C simplify to: where λ i denote the principal stretches of the shell and their squares λ 2 i are the eigenvalues of the deformation tensor C. The squares of the eigenvalues are the roots of the characteristic polynomial: ( Corresponding eigenvectors are denoted by v i , which are normalized to have unit-length. The eigenvalue decomposition (or spectral decomposition) of the deformation tensor C can be written as [24,25]: Where the Einstein summation convention is used. Since C 33 is decoupled by construction, one can immediately see from Eqs. (10) and (18) that λ 3 = √ C 33 and v 3 =n.
For the sake of completeness, we recall the definition of the Green-Lagrange strain tensor E = E αβg α ⊗g β from [10,12] and its decomposition to membrane and bending contributions (ε and κ, respectively): Remark 2. Following up on Remark 1; the contribution of the O θ 3 2 term in Eq. (7) is neglected in the strain tensor and its derivatives. The O θ 3 2 term is only included in Eq. (7) to ensure equivalence in comparison of the stretch-and invariant-based formulations.

Variational Formulation
The shell internal and external equilibrium equations in variational form are derived by the principle of virtual work [10,12]. The variations of internal and external work are defined as: with δu being the virtual displacement, δε and δκ the virtual strain components, Ω the midsurface and dΩ = å αβ dθ 1 dθ 2 the differential area in the undeformed configuration, mapped 5 to the integration domain Ω * = [0, 1] 2 using the undeformed mid-plane measure. Furthermore, with slight abuse of notation, the tensors n = n αβg α ⊗g β and m = m αβg α ⊗g β denote the shell normal force and bending moment tensors, respectively, represented by: Here, S αβ denotes the coefficients of the stress tensor following from the constitutive relations that will be derived in Sections 3 and 4 and t stands for the shell thickness. The total differentials of the stress resultants are: Discretizing the equations using known formulations from previous publications [10,12,27], the solution u is represented by a finite sum of weighted basis functions and the tensors n, m, ε and κ are linearized around the weights using Gateaux derivatives. The linearized tensors are denoted by (·) = ∂(·) ∂u r in the following, where u r are individual weights of the solution vector. Note that u denotes the basis functions [12]. Using the discretized system, the residual vector is defined by: and must be equal to the zero vector for the weights u corresponding to the exact solution. To solve the residual equation R = 0, another linearization is performed, yielding the Jacobian matrix or tangential stiffness matrix K: (24) Note that the matrix contains a contribution for the external load depending on the solution vector (f(u)). For instance, follower-pressures are defined by f(u) = pn(u), where n is the surface normal. In order to solve for nonlinear equation, Newton iterations are performed for solution u and increment ∆u by solving

Invariant-Based Constitutive Relations
Adopting an invariant-based formulation [12], the constitutive relations for hyperelastic shells are derived and illustrated for the Mooney-Rivlin material model. That is, the constitutive laws are derived using derivatives of the strain energy density function Ψ(C) with respect to the components C i j of the deformation tensor C. The aim of this section is to keep the paper selfcontained and to show the relation between the state-of-the-art [12] and the novel formulations that will be presented in Section 4.

6
The section is structured as follows: Section 3.1 provides general formulations and background that are valid for both compressible and incompressible material formulations. Thereafter, Section 3.2 and Section 3.3 provide the derivations for incompressible and compressible material models, respectively, in the invariant-based formulations. As mentioned before, these formulations are based on [12].

General Relations
Firstly, the second Piola-Kirchhoff stress tensor S = S i jg i ⊗g j and the (tangent) material tensor C = C i jklg i ⊗g j ⊗g k ⊗g l are defined by the following well-established relations [12,[24][25][26]: These equations are valid for 3D continua and hence need to be modified for the shell equations to incorporate the through-thickness stress components. Reading Eq. (10), C αβ = g αβ but C 33 g 33 = 1 to avoid violation of the plane stress condition. To correctly incorporate the plane-stress condition (S 33 = 0), the material tensor C is modified using static condensation, which implies that the material tensorĈ corrected for plane-stress is defined by [12]: For incompressible materials, this term is derived analytically using the incompressibility condition (J = 1) whereas for compressible materials, it is corrected iteratively.

Incompressible Material Models
For incompressible materials, the incompressibility condition (J = 1) is enforced using a Lagrange multiplier p in the strain energy density function [12,24]: where Ψ el is the original strain energy density function. Using Eq. (28), Eq. (26) reduces to [12]: The plane stress condition (S 33 = 0) together with the incompressibility condition (J = 1) lead to the expression for the Lagrange multiplier p, from which its derivative with respect to C i j can easily be obtained [12]: Substituting the derivatives of the Jacobian determinant J (see Eqs. (A.2) and (A.4) in Appendix A), Eqs. (29) and (30) reduce to the following equations, which are valid for incompressible materials in general with S 33 = 0 enforced [12]: For a shell, note that only the in-plane components S αβ and C αβγδ are only considered and that C −1 αβ =C αβ = g αβ , C 33 = J −2 0 . The stress-tensor and the material tensor for incompressible Kirchhoff-Love shells become, by substituting Eqs. (10) and (11): Example 1 (Incompressible Mooney-Rivlin Model). The strain energy density function of the incompressible Mooney-Rivlin material model is defined as follows [28,29]: where c 1 and c 2 are the material parameters such that the Lamé parameter µ satisfies µ = c 1 + c 2 .
Using the derivatives of the invariants I 1 and I 2 and the Jacobian determinant J with respect to C i j , the derivatives of Ψ can be found using the identities in Appendix A (see Eqs. (A.2) and (A.4) to (A.9)): and all other derivatives are equal to zero. As a consequence, Eqs. (35) and (37) reduce to: C αβγδ = c 1 + c 2 trC 2g αβ g γδ + g αγ g βδ + g αδ g βγ − 2c 2 J −2 0 g αβ g γδ + g αβgγδ (42) Note that for c 1 = µ and c 2 = 0, this model reduces to the incompressible Neo-Hookean material model (see Eq. (96)) [12].

Compressible Material Models
For compressible models, the Jacobian determinant J is not necessarily equal to 1. As a consequence, the deformation gradient F and deformation tensor C are modified such that F and C are a multiplicative decomposition of a volume-changing (dilational) part depending on J and a volume preserving (distortional) part depending on the modified deformation gradient and deformation tensors,Ċ andḞ, respectively [30]: The modified deformation gradient and deformation tensor have determinants which are equal to 1 (corresponding to volume preservation), meaning: where the modified principal stretchesλ i are defined as: As a consequence, the invariants of the modified deformation tensorĊ become: with I i the invariants of the deformation tensor C. WithḞ,Ċ andİ k defined above, the strain energy density function Ψ(C) for a compressible material can be described in a decoupled form, separating the response in a isochoric (or distortional) elastic part Ψ iso (Ċ) and an volumetric (or dilational) elastic part Ψ vol (J) [24,25,30]: The volumetric elastic part Ψ vol is required to be strictly convex and equal to zero if and only if J = 1 andĊ = I [24].
For compressible materials, the plane stress condition is incorporated by iteratively solving S 33 = 0 for C 33 using Newton linearizations [12,31]: where C 33 is incrementally updated by C (n+1) 33 = C (n) 33 + ∆C (n) 33 with the increment on iteration n: In each iteration, the updated stress tensor S and material tensor C can be computed and iterations are continued until the plane stress condition is satisfied within a certain tolerance, i.e. S 33 < tol. When converged, static condensation can be performed for the material tensor using Eq. (27). Rather than using C (0) 33 = 0 [12], C (0) 33 = J −2 0 is used for incompressible materials, although the difference for the two approaches is negligible.
where c 1 and c 2 are the material parameters such that the Lamé parameter µ satisfies µ = c 1 + c 2 .
The bulk modulus of the material is K = 2µ(1 + ν)/(3 − 6ν) and the function contribution KG(J) is the volumetric contribution function of which G is defined as [32]: The parameter β is an empirical coefficient and a value of −2 can be adopted [33]. Note that for the volumetric strain energy density function Ψ vol other functions can be selected as well [34].
Using the derivatives of the invariants I 1 and I 2 and the Jacobian determinant J with respect to C i j , the derivatives of Ψ can be obtained: where for the sake of brevity, reference is made to Eqs. (A.5) to (A.9) in Appendix A for the derivatives of the invariants and to Eq. (2) for the derivative of the inverse of the deformation tensor. Note that from Eq. (A.6) it follows that the second derivative of I 1 is equal to zero. The other derivatives of the invariants and the derivative of the inverse of the deformation tensor are not substituted for conciseness. Substituting Eqs. (53) and (54) into Eq. (26) yields expressions for the stress and material tensors: Again, for c 1 = µ and c 2 = 0, this material model reduces to the Neo-Hookean model (see Eq. (95)) [12].

Stretch-Based Constitutive Relations
Invariant-based (in)compressible material model formulations have been obtained for the strain energy density functions Ψ(C) in Section 3 based on [12]. However, when experimental material data fitting is involved a formulation in terms of stretches (i.e. in terms of the eigenvalues of C, Ψ(λ with λ = (λ 1 , λ 2 , λ 3 )) might be preferred, meaning the invariant-based formulations (Section 3) are not directly applicable. Therefore, this section provides the main contribution of this paper: the generalized formulations for the implementation of stretch-based material models in the isogeometric Kirchhoff-Love shell model.
The section is structured as follows: Section 4.1 provides the basics for the derivation of the stretch-based constitutive relations. Thereafter, Section 4.2 and Section 4.3 provide the derivations for incompressible and compressible material models, respectively, in the stretch-based formulations. These formulations are the novelty of the present paper.

General Relations
Assuming Ψ(λ), relations for the stress and material tensor will be derived in terms of the (normalized) eigenvector bases (Eq. (18)): When S and C are known, a transformation to the basesg i ⊗g j andg i ⊗g j ⊗g k ⊗g l for the stress and material tensors, respectively. This will be discussed in Section 5.3.
The derivative of any scalar function with respect to the deformation tensor C can be written as a derivative with respect to the stretch by applying the chain rule [24]: Consequently, the second derivatives are: From this, it follows that: Furthermore, it can also be shown that for the material tensor, the following holds [5,6,8,24,35]: The first part of C i jkl represents the normal components (diagonal elements) and the second part denotes the shear components (off-diagonal elements). In the formulation of the invariant-based counterpart of this equation (Eq. (26)) these parts are not explicitly visible, since the spectral form by definition uses the principal directions of the deformation tensor, whereas shear and normal contributions are mixed in the curvilinear form of the material tensor. Note that for the second part of this equation, the case λ i = λ j results in an undefined result. Hence, using L'Hopital's rule, this limit case can be identified: 11 Similar to the derivations in terms of the components (see Section 3) static condensation is also applied to enforce plane-stress. For incompressible material models, the stretch-based static condensation term will be derived. For compressible models, static condensation will be performed iteratively in spectral form and afterwards it will be transformed to the covariant undeformed basis. Since J = λ 1 λ 2 λ 3 , the derivatives of J are:

Incompressible Material Models
Similar to the invariant-based formulation (Section 3.2), the incompressibility condition J = 1 is imposed using the Lagrange multiplier p in the strain energy density function for incompressible materials, i.e.
Using Eq. (60), the stress tensor becomes: Comparing S i with the invariant-based stress tensor in Eq. (29) shows that all components can easily be obtained using Eq. (58). To comply with the plane-stress condition S 33 = 0, the equation to be solved for p using J = 1 denotes: which implies, using the derivative of J from Eq. (A.13): It can easily be shown that Eq. (67) is similar to the expression of p in the invariant-based form (Eq. (31)) using Eq. (59) and using λ 2 3 = C 33 . The derivative of the stress tensor with respect to the stretch is required to find the material tensor, as observed in Eq. (61). From Eq. (65) it follows that: where the incompressibility condition (J = 1) is used again. Note that the Kronecker delta δ j i covers the case when i = j. The derivative of p follows from Eq. (67) and yields: 12 Again, this result can be compared to its invariant-based counterpart in Eq. (32) and using Eqs. (58) and (59) Comparison with the invariant-based formulation shows that λ −1 i in front of the second term in Eq. (70) translates toC i j in Eq. (34). Using these identities, the material tensor can be derived from Eq. (61). For the static condensation term, reference is made to Eq. (27), hence the components C αβ33 , C 33αβ and C 3333 need to be evaluated. From Eq. (61) it follows that: such that the static condensation term becomes: Using this result, the in-plane incompressible material tensor can be evaluated as: where the second term needs to be replaced by Eq. (62) if λ α = λ β .
Example 3 (Incompressible Ogden Material Model). The (incompressible) Ogden material model is an example of a (class of) material model(s) that is expressed in terms of stretches. The strain energy density function reads [24,25,32]: The last equality signs show that the function can be written as a sum of three separate functions ω(λ i ) depending on different stretches. This separated form is known as the Valanis-Landel hypothesis [24,36]. Using this form, the derivatives of Ψ with respect to the stretches are: From which again the stress and material tensors can be derived using Eqs. (70), (71) and (77).

Compressible Material Models
Similar to the compressible invariant-based formulation (Section 3.3), S and C for compressible material models will be derived for stretch-based material models. This means that the thickness-correction is applied iteratively. On the stretch-based terms of the stress and material tensors, that will be derived in this section, a transformation from the eigenvector space to the undeformed covariant tensor basis is required (Section 5.3) is applied to transform the material model to the curvilinear basis. Using the transformed stretch-based quantities, the same iterative procedure as in Section 3.3 can be used to converge towards the compressible stress and material tensors. Static condensation (Eq. (27)) is performed before transforming the material tensor.
For computation of the stress and material tensors, the relations in Eqs. (60) and (61) are used. As described in Section 3.3, the strain energy density function Ψ for compressible materials can be decomposed in an isochoric and volumetric part (see Eq. (48)) [30]. For stretch-based forms, this implies: The stretchesλ i are the modified principal stretches, defined as: In this equation, Ψ vol can be chosen according to Eq. (52) but is not specified in the further 14 derivation. Using this form, the derivatives of Ψ with respect to the stretches are: where Λ = λ

Implementation Aspects
In order to provide some implementation aspects for the presented formulations (Section 4), the assembly of the nonlinear system for isogeometric Kirchhoff-Love shells will be recalled (Section 5.1) as well as the computation of the eigenvalues and eigenvectors of the deformation tensor C (Section 5.2). Lastly, we will provide details about the transformation of the stress and material tensors S and C from spectral to curvilinear bases (Section 5.3).

System Assembly
For the implementation of Kirchhoff-Love shells recall that the vector of internal forces and the tangential stiffness matrix read [10,12]: Here, we note that the matricesD k , k = 0, 1, 2, are k th thickness moments of the material tensor represented as a 3×3 matrix andn andm are the zero-th and first thickness moments of the stress tensor, see [12]. The thickness integrals are, in the present paper and in [12], computed using numerical through-thickness integration with four Gaussian points. As discussed before by [15], the matricesD 1 can differ in the variations of the normal force tensorn and the moment tensor m depending the analytic projected or directly decoupled alternatives for thickness integration.

Eigenvalue Computation
The eigenvalues a tensor quantity can be computed by solving Eq. (17) or, alternatively, by computing the eigenvalues of the matrix that results from symbolic computation of C = C i jgi ⊗g j . Since λ 2 3 = √ C 33 is decoupled by construction, the other two eigenvalues λ 2 1 and λ 2 2 are to be computed. Here, an basic routine for the eigenvalues and eigenvectors for 3 × 3 matrices is employed on the matrix that follows from explicit computation of C = C αβgα ⊗g β .

Tensor Transformation
Since the stretch-based stress and material tensor are derived in spectral form (i.e. in the eigenvector space) a transformation towards the curvilinear basis needsis required in order to use these entities in further computations. Recall that the spectral forms of S and C are: The the invariant-based stress and material tensors are defined in the curvilinear basis, as follows: Since the strain tensors (c.f. Eq. (19)) are defined in the curvilinear basis, it is convenient to define the quantities in the variational form (c.f. Eq. (20)) defined in the curvilinear basis. Hence, the stretch-based stress and material tensors have to be transformed to the undeformed covariant curvilinear basis by: whereS i j andC i jkl are the coefficients of the stress and material tensors in the curvilinear basis.
Obviously, the tensor transformation only needs to be computed for non-zero components of C pqrs . For incompressible material models, static condensation is applied analytically, which implies that the transformations only need to be applied for indices ranging from α, β, γ, δ = 1, 2, thus the transformation consists of mapping 2 4 = 16 entries. However, it is known that for hyperelastic materials the contravariant components of the material tensor, C i jkl , posses minor and major symmetry [24,25], i.e.
C abcd = C bacd = C abdc minor symmetry, = C cdab major symmetry, So that only six unique components exists for the 2 × 2 × 2 × 2 tensor. Furthermore, Eq. (61) implies that the non-zero components of C i jkl are of the form C iiii , C ii j j , C i ji j and C i j ji of which the last two are equal by virtue of the minor symmetry property. This implies that the 2 × 2 × 2 × 2 tensor has only four uniquely defined components, namely C 1111 , C 1122 , C 2222 and C 1212 .
For compressible material models, static condensation is applied in the spectral basis, i.e. on the tensor C before it is transformed to the covariant undeformed tensor basis. From Eq. (27) we see that the static condensation procedure requires the computation of C 3333 , C αβ33 and C 33αβ , where the last two are equal by virtue of the major symmetry property. Using the minor and 16 major symmetries again, this implies computation of four extra (unique) components, namely C 1133 , C 2233 , C 1233 and C 3333 .
Accordingly, it can be concluded that for incompressible materials four and for compressible materials eight unique components of the spectral material tensor need to be computed, when exploiting minor and major symmetry, as well as the nature of Eq. (61). Hence, the additional costs of the transformation are existent, but limited due to symmetries of the spectral form.

Numerical experiments
For benchmarking purposes, the results of five numerical experiments have been used for verification and validation of the presented formulations for incompressible and compressible material models. For the uniaxial tension and pressurized balloon benchmarks (Sections 6.1 and 6.3, respectively), analytical solutions are available, meaning the results will serve as verification of the stretch-based material model formulations. The restrained sheet elongation (Section 6.2) will serve as validation of the stretch-based material model formulations because of the involved numerical (IGA and FEM) and experimental data. The benchmark of the pinched cylinder (Section 6.4) provides additional comparison of selected material models with reference works. In order to verify the presented isogemetric Kirchhoff-Love formulation for a stretch-based Ogden material with its FEM couterpart, the conical shell collapse (Section 6.5) is incorporated. The described models have been implemented in the open-source library G+Smo (Geometry + Simulation Modules) [37].
In the numerical experiments, compressible and incompressible formulations of the Neo-Hookean (NH), Mooney-Rivlin (MR) and Ogden (OG) material models have been used. The Neo-Hookean models are given by: For the MR and OG models, the established formulations (respectivelyEqs. (38) and (51) and Eqs. (78) and (83)) are used. For all models, the following volumetric part of the strain energy density function is adopted (see Eq. (52)): For the NR and MR models, the invariants can be replaced by Eqs. (14) to (16) to obtain stretchbased forms (see also Appendix B). Unless stated otherwise, for the compressible models β = −2, and for the Mooney-Rivlin model c 1 /c 2 = 7 [29] is used. For the Ogden model the coefficients from [38] are re-scaled to the value of µ: L W x y σt where µ 0 = 4.225. Lastly, for obtaining the analytical solutions, we note that the principal Cauchy stresses are defined as :

Uniaxial Tension
The first benchmark case is uniaxial tension of a material block. In this case, a block with dimensions L×W ×t = 1×1×0.001[m 2 ] is considered. The shear modulus is µ = E/(2(1+ν)) where E and ν are the Young's modulus and Poisson ratio, respectively, such that µ = 1.5 · 10 6 [N/m 2 ]. The block is modeled by shell elements, which means that the L × W plane is considered and all edges are restrained in vertical direction (z = 0). The left edge (x = 0) is restrained in x direction and on the right edge (x = L) a distributed load σt is applied. The bottom edge (y = 0) is restrained in y direction and the top edge (y = B) is free to move (see Fig. 1.
For incompressible materials, the analytical solution is obtained by the incompressibility constraint J = λ 1 λ 2 λ 3 = 1 and λ 1 = λ as the stretch of the block. As a consequence, λ 2 = λ 3 = 1/ √ λ are the changes in width and thickness of the block and the Cauchy-stresses are σ 2 = σ 3 = 0 and Neo-Hookean: For compressible materials, the incompressibility constraint is not valid, hence J = λ 1 λ 2 λ 3 1, and λ 1 = λ and λ 2 = λ 3 = √ J/λ. The plane-stress condition σ 3 = 0 using Eq. (100) is adopted to obtain J for any value of λ. Using this result, the Cauchy-stress for any λ can then be obtained Table 1: Residual norms per iteration for the 10 th load-step for uniaxial tension for all material models in compressible and incompressible forms. For the Neo-Hookean and Mooney-Rivlin models, the iteration residuals are provided for the stretch-based and invariant-based implementations. For the Ogden model, only the results for the stretch-based implementation are given, since no invariant-based formulation exists. For the Neo-Hookean and Mooney-Rivlin models, results are only observed in the last iteration, due to machine precision of the arithmetic.

Iteration
Neo using: Neo-Hookean: Ogden: In Fig. 2 the results for uniaxial tension are depicted. The numerical and analytical solutions for incompressible and compressible materials show a perfect match for all quantities studied (thickness decrease λ 3 , axial Cauchy stress σ and Jacobian determinant J). Note that the Jacobian determinant for incompressible materials is equal to 1 and hence not shown. The residual norms of the non-linear iteration convergence for the invariant-based and stretch-based Neo-Hookean and Mooney-Rivlin models as well as the stretch-based Ogden model are equal in all iterationsTable 1, showing that the present formulation provides exactly the same rates of convergence as the invariant-based method. Last but not least, the iterations converge with the second order, as can be expected when using Newton iterations.

Restrained Sheet Elongation
Related to the first benchmark in the work of [15] and on the experiments of [39], a tensile load is applied on a strip of which the short edges are fixed and the long edges are free (see Fig. 3). Focus is on the non-domensional load versus end-point displacement in longitudinal (load and displacement) direction. For incompressible materials, comparison is made with the numerical results of [15] and for compressible materials, comparison is made with the experimental and numerical results of [39]. In both cases, the geometric parameters of the sheet are L = 157.895 [mm], W = 78.947 [mm] and t = 0. 15 [mm] leading to L/W = 2 and t/W = 1.90·10 −3 as in [39]. Furthermore, for comparison with [15] Young's modulus E = 1[Pa]], and a Poisson's ratio ν = 0.5 (due to incompressibility) are used, which leads to µ = 1/3 [Pa]. For the Mooney-Rivlin model, c 1 /c 2 = 1/2 such that c 1 = 1/9 and c 2 = 2/9. Scaling according to A good match with the results of [15] for the incompressible Neo-Hookean and Mooney-Rivlin models can be seen in Fig. 4. Note that the forces in the reference paper are normalized by E = 3c 1 for both the Neo-Hookean and Mooney-Rivlin models, whereas in the present simulations, the are normalized by E = 3µ (since ν = 0.5 in the comparison with [15]). Furthermore, the deflection contour at the sheet sides (Fig. 5) shows a good agreement with the numerical/experimental results of [39] as well. The curves are presented for different strains which shows that dispite strain-normalization on the vertical axis strain dependence is still visible, meaning that there is non-linear dependence of the strains in the scaling laws.

Pressurized Balloon
The response affected pressure of a spherical balloon is used for benchmarking purposes as well. The analytical pressure formulation denotes [24]: where h and r are the undeformed thickness and radius of the balloon, respectively, σ 1 = σ 2 = σ since λ 1 = λ 2 = λ and λ 3 = 1/λ by the incompressibility condition. Substituting in Eq. , using the compressible Neo-Hookean model. The horizontal axis denotes the non-dimensional location on the length of the sheet and the vertical axis denotes the contraction over the width-axis, made nondimensional using the applied strain ε and the width of the sheet. The markers denote the experimental data from Ref. [39] and the markers are grouped on strain level (ε ∼ 0.09 − 0.30) and colored using the same colors as in their work. The black dashed line denotes the result of the FEM model used by [39]. The other lines represent present results for different meshes on ε ∼ 0.1. material model, this yields: The numerical model results are based on follower pressures, i.e. f = pn where n is the unit normal in the current configuration. The balloon is modeled as a quarter sphere, of which the bottom point is fixed in all directions, and on the sides a symmetry condition is applied by clamping the sides in normal direction and restriction deflections orthogonal to the symmetry boundary (see Fig. 6). The geometry is modelled by 4 elements over the height and 2 element over the quartercircumference, both of quadratic order.

Pinched Cylinder
The response of a pinched cylinder (Fig. 8) [12,[40][41][42] is also benchmarked. Although the response is bending dominated -meaning the influence of Ogden materials is expected to be minimal -the present model is benchmarked for a variety of material models for the sake of completeness. The geometry of the problem is depicted in Fig. 8. The dimensions include   4) show marginal differences with the references and the influence of higher-order nonlinear material models (MR and OG) show minor differences as well; this is expected behaviour due to the bending dominated response. The minor difference with the isogeometric Kirchhoff-Love shell model of [12] is explained by the fact that the mesh size of the present study is slightly finer (16 vs. 12 elements over the half-circumference).    [12] uses 16 × 12 quartic isogeometric Kirchhoff-Love shell elements and models half of the cylinder. Ref [41] uses (a) a 7-parameter model with exact strains; and (b) 6-parameter model with exact strains and an multiplicative strain-decomposition used in the incompatible mode method. Ref [42] employ solid-shell finite elements with (a) 1 element over the tickness; and with (b) 2 elements over the thickness. Ref [40] represents solutions obtained using finite elements with quadratic strains and a 7-parameter theory with condensation.

Conical Shell Collapse
A collapsing conical shell (or frustrum) is presented as a benchmark for modelling of strong non-linearities [8].  Fig. 9 is considered. The geometry is modelled with 32 quadratic elements over the height and one quadratic element over the circumference to represent axial symmetry. The corresponding material model is of the Ogden type and has the following parameters:  1 ) is either kept rigid (no x and y displacements) or free, referred to as constant or variable radius, respectively [8]. On the top edge, a uniform load p is applied, providing a uniform displacement ∆. Due to symmetry, only one quarter of the geometry is modelled, which means that symmetry boundary conditions are applied on the x = 0 and y = 0 planes (Γ 3 , Γ 4 , see Fig. 9). The quarter-conical shell is modelled with 32 quartic shell elements over the width.
Loads are applied using displacement-control (DC) or arc-length control. In the former case, displacements are applied on the top-side of the cone and the deformation of the cone as well as the corresponding load on the top-boundary are computed. In the latter case, Crisfield's spherical arc-length procedure is used with extensions for resolving complex roots. If this method does not converge to an equilibrium point, the step size is bisected until a converged step is found. After this step, the step size is reset to its original value [43].
Figs. 11 and 12 present the result of the collapsing conical shell (constant and variable radius, respectively) of the present study and the reference results from [8]. The results for the displacement-controlled (DC) solution procedure shows that the different between the used material models are negligible, since the actual strains are relatively small. The results also agree with the displacement-controlled reference results of [8], and minor differences between the results might be a result of FE shear locking as involved for the reference results. Since more steps have been used for the displacement-controlled calculations, sharp corners in the curve can be observed for ∆ ∼ 1.9 for constant radius and ∆ ∼ 1.8 for variable radius.
An arc-length based calculation was used as well. From the results, one can observe revelation of the collapsing mechanism of the conical shell. For both cases (constant and variable radius) an almost anti-symmetric pattern in the load-deflection space can be observed, which initates and finishes at the kinks in the curve that was found with the displacement-control procedure. For the constant-radius shell, Fig. 11a shows two loops of large magnitude. In Figs. 11b and 12b it can be seen that collapsing behavior of the conical shell consists of states in which multiple waves in radial direction occur. For both cases, it can be seen that after the loops with the highest force-amplitude, the shell and its collapse-path invert and continue on the path that can be obtained with displacement-control. Since variation between the material models is rather small for the DC solutions, only the results for the OG material model are given. The reference results are obtained from [8]. A movie of the collapse is available as supplementary material to this paper (video 2).
To the best of the author's knowledge, the collapsing of a collapsing conical shell was not investigated before. Complex load-displacement paths from Figs. 11a and 12a show that displacement-controlled simulations in this case ignore the collapsing behaviour of the shell with multiple limit points. The authors highly encourage further investigations on this benchmark for verification and validation.

Conclusions and recommendations
A stretch-based hyperelastic material formulation for isogeometric Kirchhoff-Love shells is presented, particularly useful when experimental material data fitting is involved to capture the model parameters. Complementing the existing invariant-based formulation [12], the implementation is based on a spectral decomposition of the stress and material tensors, which are transformed to the covariant curvilinear basis for consistency in the variational formulation. The eigenvalue computation and the transformation of the spectral entities imply additional computational costs, although limited due spectral material tensor symmetries. At the same time, the possibility to translate the stretch-based formulations into invariant-based ones is considered to be an advantage since it is not possible the other way around. There is still room for optimization of the eigenvalue computation, for instance by only evaluating the eigenvalues at the mid-plane of the shell.
For benchmarking purposes, incompressible and compressible Neo-Hookean, MOney-Rivlin and Ogden models have been considered; the latter being only defined in terms of stretches. Identical iteration residuals and correct convergence rates have been obtained. The analytical benchmarks (Sections 6.1 and 6.3) have shown a perfect agreement for the Ogden material model using the stretch-based formulation. Very good agreement was obtained for the numerical benchmarks. Benchmarking the response of a collapsing cylindrical shell, the displacement-controlled results are in good agreement. At the same time, a complex collapsing mechanism -not observed in literature before -was revealed when using the arc-length method.
Inspired by the results of the uniaxially stretched restrained sheet, we aim to further apply the present model on the modeling of wrinkling in our future work. The presented stretchbased formulation performance could be improved applying an analytical projection and direct decoupling [15] of the constitutive equations in order to prevent for numerical through-thickness integration (i.e. eigenvalue computations for all through-thickness Gaussian points). Firstly, the derivatives of J are derived based on the derivative of J 2 [24], that is: which implies the derivative of J is: and the derivative of an arbitrary power k ∈ R of J is found by the chain rule: The second derivative of J is simply evaluated by taking the second derivative of Eq. (A.2). This gives: where we used the derivative of an inverse of a tensor, from Eq. (2). From the definitions of the invariants in Eqs. (14) to (16)