{'fraction_non_alphanumeric': 0.03708192518354504, 'fraction_numerical': 0.02153595152080177, 'mean_word_length': 4.239467578458909, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new anisotropic bending model for nonlinear shells: Comparison with existing models and isogeometric finite element implementation

A new nonlinear hyperelastic bending model for shells formulated directly in surface form is presented, and compared to four prominently used bending models. Through an essential set of elementary nonlinear bending test cases, the stresses and moments of each model are examined analytically. Only the proposed bending model passes all the test cases while the other bending models either fail or only pass the test cases for small deformations. The proposed new bending model can handle large deformations and initially curved surfaces. It is based on the principal curvatures and their directions in the initial configuration, and it thus can have different bending moduli along those directions. These characteristics make it flexible in modeling a given material, while it does not suffer from the pathologies of existing bending models. Further, the bending models are compared computationally through four classical benchmark examples and one contact example. As the underlying shell theory is based on Kirchhoff-Love kinematics, isogeometric NURBS shape functions are used to discretize the shell surface. The linearization and efficient finite element implementation of the proposed new model are also provided.


Introduction
Shells are curved thin-walled structures appearing in nature and engineering designs. The curvature enables shells to be designed with high load-bearing capacity at minimal use of materials. This high strength-to-weight ratio makes them ubiquitous in many applications. Extensive efforts have been made to accurately describe the load-carrying behavior of shells. The finite element method (FEM) is the predominantly used numerical technique to solve shell problems. It is common practice in FEM formulations of slender structures such as plates, membranes, and shells, to reduce the dimension from volume to surface. This significantly simplifies the numerical discretization, reduces the degrees-of-freedom (dofs), and condenses the 3D kinematics to 2D. Savings in computational time naturally follow.
The three mainly used ways of dimensionality reduction are derived, degenerate, and direct surface approaches (Bischoff et al., 2004). In the derived approach, an approximate shell theory is derived by asymptotic (analytical or numerical) integration of the 3D equations. Some examples of the early works, which mainly focused on linear theories, are those by Gol'denveizer (1963), Reissner (1963) and Cicala (1965). In the degenerate solid approach, the shell behavior is obtained by reducing or degenerating 3D continuum mechanics through kinematic assumptions. This was first developed by Ahmad et al. (1970) and a comprehensive presentation of the methodology can be found in Hughes (2012). In the direct surface approach, only the shell surface is considered ab initio and well-defined constitutive laws are proposed to obtain stresses and moments. The underlying theory of this method goes back to Cosserat (1909). This was followed by prominent works of Ericksen and Truesdell (1957), Green et al. (1965), Naghdi (1973) and Simo et al. (1990).
In addition to the dimensionality reduction, the underlying shell kinematics also plays a vital role in constructing shell theories. The two widely used shell theories are Reissner-Mindlin (RM) (Reissner, 1945;Mindlin, 1951) and Kirchhoff-Love (KL) (Kirchhoff, 1850;Love, 1888). In KL theory the cross-section remains normal to the mid surface during deformation, whereas in RM shell theory, a shear angle can appear. Therefore, RM theory has both displacement and rotational degrees-of-freedom (dofs), and C 0 -continuous shape functions suffice for discretization. On the other hand, KL theory has only displacement dofs. It therefore accommodates bending in the governing equations through displacements itself. This results in a fourth-order strong form equation for KL shells. The corresponding principle of virtual work contains second order derivatives which necessitates C 1 -continuous shape functions for discretization. There are several methods developed to enforce C 1 -continuity for Lagrange shape functions like rotation-free elements (Oñate and Zárate, 2000;Brunet and Sabourin, 2006), discontinuous Galerkin formulation (Noels and Radovitzky, 2008) and mesh-free methods (Krysl and Belytschko, 1996). However, these methods are usually complex and/or expensive. Therefore RM based shells are commonly used in commercial finite element (FE) codes as they can be used with simple classical Lagrange shape functions.
A new approach to obtain C 1 -continuous shape function called Isogeometric analysis (IGA) was introduced by Hughes et al. (2005). There the spline-based basis functions used to create the geometry, such as B-splines or Non-Uniform Rational B-Splines (NURBS), itself are used for the FE analysis. The Bézier extraction operator developed by Borden et al. (2011) enabled a NURBS surface to be decomposed into Bézier elements to seamlessly incorporate isogeometric analysis into existing FE code. The advantages of IGA over conventional approaches in shell problems were shown by Kiendl et al. (2009) for KL and Benson et al. (2010) for RM shells. A KL shell formulation doesn't just benefit from having fewer dofs compared to a RM shell, but it is also shear locking free. However, both KL and RM shells suffer from membrane-bending locking. This type of locking arises due to the undesirable coupling of membrane and bending modes. A variety of solutions exists for alleviating membrane-bending locking in IGA shells, such as reduced quadrature (Adam et al., 2015;Zou et al., 2021), mixed formulations (Bieber et al., 2018) and assumed natural strain methods (Caseiro et al., 2014). Kiendl et al. (2015) presented an isogeometric KL shell formulation such that an arbitrary 3D constitutive models can be used. Along similar lines, Duong et al. (2017) proposed a formulation that admits constitutive laws obtained by thickness-integration of 3D material models as well as those constructed directly in surface energy form. IGA shell formulations have been applied to the study of rubbers (Taylor, 2011;Elguedj and Hughes, 2014), tissues (Tepole et al., 2015;Roohbakhshan and Sauer, 2017), red blood cells (Casquero et al., 2017;Bartezzaghi et al., 2019), lipid bilayers Auddya et al., 2021) and viscoelastic materials (Dortdivanlioglu and Javili, 2021;Paul and Sauer, 2022). IGA shells have also been used in inverse problems such as shape optimization  and material reconstruction (Borzeszkowski et al., 2022).
One of the material models for the direct surface approach was proposed by Koiter (1966). It is a linear relation between certain strains and stresses that can be derived from thickness integration of the 3D St.Venant-Kirchhoff material model. The Koiter model was later extended by Steigmann (2012) to materials exhibiting arbitrary symmetry. A comprehensive derivation can also be found in Steigmann (2013). Two other direct surface material model were proposed by Canham (1970a), to study red blood cells, and Helfrich (1973), to study the elastic properties of lipid bilayers. The Helfrich model, which includes the Canham model as a special case, is one of most widely used material models in morphological studies of the vesicles, which are closed bilayer films.
KL shell formulations are primarily used with bending models derived from 3D or directly proposed for the surface such as the Koiter and Helfrich model. There aren't many examples of direct bending models besides those of Koiter and Helfrich. As is shown here, these suffer from an inability to handle nonlinear deformations or initially stress-free curved surfaces. We address these shortcomings by proposing a new nonlinear anisotropic bending model for the direct surface formulation of shells. It is motivated by eliminating the physical stretch-dependency on bending that is affecting existing bending models at large deformations. This is different from treating membrane-bending locking, which is caused by a numerical stretch-dependency of bending of very thin shells that is already a problem at small deformations. Membrane-bending locking is due to shortcomings in the underlying finite element discretization and therefore needs to be treated at that level, as was noted above. Our concern is purely physical and therefore needs to be treated at the constitutive level. Our description uses a direct surface formulation for KL shells in curvilinear coordinates together with an isogeometric surface discretization, which is very general and accurate, yet straightforward to implement .
The salient features of the proposed new bending model are: • It is objective, admits large deformations and captures anisotropic bending.
• It allows to describe initially curved stress-free shells.
• It offers great flexibility in modeling a given material.
• It is compared to existing bending models and it agrees with those at small deformations.
• In contrast to existing models, it passes a proposed set of nonlinear bending test cases.
• It allows for an efficient implementation within isogeometric shell FE.
The remainder of this paper is organized as follows: Sec. 2 provides an overview of general thin shell theory in curvilinear coordinates. In Sec. 3, existing bending models are presented along with the proposed new bending model. In Sec. 4, these models are then compared and their shortcomings are illustrated using five analytical test cases. In Sec. 5, the new bending model is examined in several numerical examples. The paper concludes with Sec. 6.

Thin-shell formulation
This section summarizes the Kirchhoff-Love thin-shell theory formulation of  and Duong et al. (2017) that is based on curvilinear coordinates and isogeometric finite elements. First, the essential kinematical relations, governing strong and weak form equations, and hyperelastic constitutive equations are introduced in Sec. 2.1-2.3. Then, the linearization and finite element approximation of the resulting nonlinear equations are presented in Sec. 2.4 and 2.5.

Surface kinematics
The current and initial configuration of a shell surface embedded within 3D space can be parametrized respectively as with parameters ξ 1 and ξ 2 . A basis in the curvilinear setting can then be defined by a pair of covariant tangent vectors (a α , α = 1, 2) and normal (n) as A similar set of basis vectors is defined for the initial configuration as Since the covariant tangent vectors are not orthogonal, i.e. the covariant surface metric a αβ = a α · a β = 0, its dual vectors are introduced by such that a α · a β = δ β α and [a αβ ] := [a αβ ] −1 . The surface stretch is given by J = J a /J A where J a := a 1 × a 2 and J A := A 1 × A 2 is the area of the parallelogram enclosed by the covariant tangents in the current and initial configuration, respectively. The components of the surface curvature can be defined as Analogous curvature definitions, denoted B αβ , B α β and B αβ , follow for the initial configuration. The two invariants of the curvature tensor, b = b αβ a α ⊗ a β , called mean and Gaussian curvature are calculated respectively as The principal curvature of the surface can then be calculated as Additionally, we define two symmetric tensors: The Green-Lagrange strain tensor, and the relative curvature tensor, A more detailed description of thin shell kinematics and its variation can be found for example in Sauer (2018).

Strong form and weak form
The quasi-static shell boundary value problem governs the displacement field u of the surface S through the fourth order partial differential equation where f is a body force on S, and the prescribed displacement (ū), traction (t), and bending moment (m τ ) boundary conditions In Eq. (11), is the stress vector. Here, N αβ and S α are the in-plane membrane and the out-of-plane shear stress components defined via Cauchy's theorem. From angular momentum balance follows where M αβ are the stress couples caused by out-of-plane bending. Tab. 1 gives an overview of the different stress components. As shown, they can be either expressed per current surface area or per reference surface area.
per current area per reference area Multiplying strong form Eq. (11) with a suitable variation δx ∈ V and integrating it over surface S, leads to the weak form where Eq. (18) contains the virtual work of the in-plane membrane deformations and out-of-plane bending deformations. It can be extended to out-of-plane shear and strain deformations (Simo et al., 1990) and in-plane bending deformations (Steigmann, 2018;Duong et al., 2020) in the context of more general shell theories. The last part of Eq. (19) represents the virtual work of moment m τ ν, with ν being the normal to the boundary where bending moment m τ is applied.

Surface Hyperelasticity
For hyperleastic materials, the membrane stress and the moment components in weak form (17) are calculated from the surface energy density W (with units [J/m 2 ]) by It is convenient to decompose the surface energy density into two parts, associated with membrane and bending deformations. Ideally, the former should only generate membrane stresses, while the latter only generates bending stresses. However this is not possible in the case of coupled membrane-bending material behavior. We proceed by introducing some simple hyperelastic membrane models below. Hyperelastic bending models are then discussed in detail in Sec. 3.
An example of a linear elastic membrane strain energy is the Koiter model. Its membrane strain energy is of the form where the fourth order tensor C is defined as Here, Λ and µ are the surface Lamé parameters (with units [N/m]). Inserting Eq. (23) into (22) leads to where tr E = E αβ A αβ , and E αβ : The Kirchhoff stress and in-plane moment contributions of the Koiter membrane model are obtained from Eq. (20) as τ αβ mKo = Λ tr EA αβ + 2µE αβ , M αβ 0 mKo = 0 .
Another commonly used membrane model, which can be derived from 3D elasticity, is the Neo-Hookean surface model. Its strain energy function can be written as Here, I 1 = A αβ a αβ is the first invariant of the surface Cauchy-Green tensors. The Neo-Hookean model only produces the membrane stress i.e. M αβ 0 mNH = 0.

Linearization
In order to solve the nonlinear Eq. (17) using the Newton-Raphson method, its linearization is necessary. This leads to the increment for the internal virtual work where the material tangents are defined as Examples for these, along with ∆δε αβ and ∆δb αβ can be found in Sauer (2018).

Isogeometric FE approximation
The geometry within an undeformed finite element Ω e 0 and its deformed counterpart Ω e is interpolated from the positions of control points X e and x e , respectively, as with Here {N A (ξ α )} ne A = 1 are the n e C 1 −continuous NURBS basis functions of Ω e . Such an arrangement results in a direct correspondence between the discretized equations and their implementation in a computer code. Discretizing Eq. (17) gives where n el is the number of elements. Based on Eq. (18) and (31), the internal forces due to membrane stress τ αβ and the bending moment M αβ 0 will be where and Here, N e,α and N e,αβ denote the first and second derivatives of N e w.r.t. parameter ξ α . Discretizing Eq. (19), the external virtual work follows as where the external FE force vectors are for the body force f = f 0 + p n. Here f 0 denotes dead loading, while p denotes an external pressure always acting normal to the surface, and ν = ν a a α denotes the boundary normal on ∂S.
Discretizing Eq. (29) gives where k e τ τ := are the material stiffness matrices and with are the geometric stiffness matrices. The discretization of Eq. (38) leads to where k e extp and k e extm can be found in Sauer et al. (2014) and Duong et al. (2017), respectively.

Bending constitution
This section presents hyperelastic bending models for thin shells. We first summarize existing bending models and then propose a new nonlinear, anisotropic bending model. In each subsection we state the surface energy per unit area and then define the stress and moment calculated using Eq. (20).

Koiter bending model
The surface energy density for the Koiter bending model (Ciarlet, 2005;Steigmann, 2013) is where the fourth order tensor F is defined w.r.t. C from Eq. (23) as Here, T is the shell thickness. Inserting F into Eq. (45) leads to where tr K = K αβ A αβ , and K αβ : while τ αβ bKo = 0. The surface Lamé parameters can be derived from Young's modulus (E) and Poisson's ratio (ν) as (Ciarlet, 2005) that admit the special case ν = 0.5. From these, one can identify the bending stiffness (Landau and Lifshitz, 1986)

Canham bending model
The surface energy density for the bending model of Canham (1970b) is of the form where c is the bending modulus with units [Nm]. According to Eq. (20), the Canham bending model contributes to both stress and moment as As seen, M αβ 0 Can is proportional to the total curvature b αβ instead of the relative curvature K αβ as for the Koiter model in Eq. (49). This makes the Canham model unsuitable for initially curved stress-free shells. Comparing their stiffness f αβγδ , further shows that the Canham and Koiter model are equivalent for small deformations of initially flat shells when implying according to Eq. (51).

Helfrich bending model
The surface energy density for the bending model of Helfrich (1973) is defined by whereH 0 is the spontaneous curvature of the material, which can be an externally prescribed quantity, or taken as the initial mean curvature of the shell. Further, k andk are material constants. The stress and moment follow from Eq. (20) as with ∆H := H −H 0 . This model reduces to the Canham model forH 0 = 0, k = 2c, and k = −c.
Comparing the stiffness f αβγδ , shows that the Helfrich and Koiter bending model are equivalent in the linear regime if which then implies c Hel = c Koi .

Analytically projected Neo-Hooke bending model (apH)
The classical 3D Neo-Hooke material model can be analytically integrated through the shell thickness as presented in Duong et al. (2017). Assuming incompressibility (ν = 0.5) this leads to the stress and moment Here, H 0 refers to the mean curvature in the initial configuration. Model (59) becomes a pure membrane model if M αβ 0 apH = 0, i.e. for T 2 → 0. Otherwise it is a complete shell model, meaning, it has both membrane and bending parts unlike the other models considered in this section. The bending moment parameter, c apH , is the same as c Koi in Eq. (51) as long as ν = 0.5.

Proposed new bending model
The preceding bending models all have drawbacks, as will be shown in Sec. 4 and 5. This motivates the following new bending model defined by the surface energy density per reference area where the kinematical quantities characterize the relative curvature, and c i , c 12 and c 3 are bending moduli. The stretch (λ i ) and curvature (κ i , κ 0i , κ 12 and κ 012 ) measures are defined as (without summation on i) and It is emphasized that these kinematic quantities do not depend on the choice of surface parameterization, hence making the new model frame invariant. Here κ 0i , for i = 1, 2, refers to the two principal surface curvatures in the initial configuration, while refers to the corresponding principal curvature directions in the initial configuration. As a consequence of Eqs. (63) and (64), λ i and κ i are generally not equal to the current principal stretches and curvatures. The proposed bending model is a generalization of the 1D fiber bending model of Duong et al. (2020) that appears for c 2 = c 12 = c 3 = 0. In contrast to the Canham model, the proposed bending model of Eq. (60) can be employed for initially stress-free curved shells.
The components of the membrane stress and bending moment follow from Eqs. (20) and (60) as The fourth order material tangents and an efficient FE implementation for the new bending model are provided in Appendices A and B, respectively. Appendix C shows that for small deformations the proposed model is equivalent to the Koiter model, if c 1 and c 2 then play the same role as c Koi .

Elementary bending test cases
In this section, three analytical test cases are investigated, and the results for the proposed new bending material model are juxtaposed with the results for the other bending models given in the previous section. In all the test cases considered, the initial configuration of a tube with radius R and length L = π R is considered. The initial surface of the tube can be parametrized using to represent any point on it as with e r = cos θ e 1 + sin θ e 2 , The tangents and normal for this surface then follow from (3) as and Using these, we can obtain the further quantities and In the elementary test cases that follow, we use Λ = 0 and µ T 2 /6 = c for all the models except the proposed new model. For the new bending model, the bending parameters are considered to be arbitrary and c 1 = c 2 = c. This generalization is used so that all the bending models can be compared with the proposed model through a common material parameter. Further, for the Helfrich model, the corresponding material constants are taken from Eq. (58) and the spontaneous curvature is set to the initial mean curvature valueH 0 = −1/(2R).

Test case 1 -bending vs. rigid rotation
As shown in Fig. 1, the first test case compares the stresses and moments for two similarly looking but differently obtained final configurations. One is obtained by applying a rigid body rotation and thus no stress or moment should be induced. The other is obtained by counter bending, wherein stress and moment are induced. This test examines whether the principal curvatures are assigned properly during deformation. This is essential for bending models directly defined on the principal curvatures, such as our proposed model.

Rigid body rotation
Given the parametrization of Eq. (70), the current configuration of the half tube rotated by 90 • around the e 2 -axis (Fig. 1a) is described by The tangents, normal, surface metric and curvature tensor components thus are a 1 = −R sin θ e 1 − R cos θ e 3 , a 2 = R e 2 , n = cos θ e 1 − sin θ e 3 , and into the equations of Sec. 3. The results are enumerated in Table 2. As seen, for the Koiter, apH and proposed new model no stress or moment are introduced. However, for the Canham model we obtain both non-zero stress and moment and for the Helfrich model, a non-zero moment. This is because for the considered initial configuration, these two models are not completely stress-free. The Helfrich bending model only satisfies this test fork = 0 as otherwise the bending stress M αβ 0 Hel becomes non-zero even when ∆H = 0 as seen in Eq. (57.2).

Counter bending
A similarly looking final configuration as before can be obtained by the counter bending shown in Fig. 1b. In this case the current configuration is given by The required kinematic quantities then are and Unlike the rigid body rotation considered before, the bending described by Eq. (88) generates non-zero bending moment components. Given the nature of the final configuration in this test case, the bending moment induced here should satisfy M 1 1 = −M 2 2 , while the Cauchy stress components, N α β , should still be zero. The stress and moment components of the five bending models are shown in Tab. 3. The Canham model neither provides accurate bending moment nor Cauchy stress components. While the Helfrich model gives accurate Cauchy stress and erroneous bending moment, the apH and Koiter models give accurate bending moments but not Cauchy stresses. Only if κ 2 is small (corresponding to small deformations) κ 2 2 becomes negligible, and hence satisfactory Cauchy stresses are obtained. Thus the Koiter and apH model can be argued to partially pass the test. But only the new model fully passes the test. The test also illustrates its anisotropic nature: If the bending parameter c 1 = c 2 then M 1 1 /(c 1 − c 12 ) = − M 2 2 /(c 2 − c 12 ).

Test case 2 -bending vs. stretching
Test case 2 also examines the behavior of two similarly looking yet different deformations. They are now obtained either by inflation or bending, as Fig. 2 shows.

Inflation
If the half tube is inflated to have radius r (Fig. 2a), the current configuration is given by based on Eqs. (70) and (72). While the components of the curvature tensor are the same as in Eq. (78), but with R replaced by r, the surface metric now is The stretch and curvature terms are then of the form and The resulting stress and moment components for test case 2(a) are listed in Tab. 4. The deformation considered here is a pure membrane action, and hence any stress contribution should only stem from the membrane part but not from the bending part of the material model. This is achieved by our new bending model as it does not introduce any stress or bending moment even though the surface curvature changes. Also the apH model achieves this for τ αβ (which comes from the membrane stiffness µ), but not for M α β and N α β (which are affected by bending stiffness c). All other models show unphysical stresses in all components. In case of the Koiter and apH model, those vanish in the limit r → R , i.e., for small ∆κ 1 := κ 1 − κ 01 . In case of the Canham and Helfrich bending models non-zero stresses remain even for small ∆κ 1 .  Table 4: Stress components for test case 2(a) -inflation -according to various bending models. In this case, all bending induced stresses (τ α β , M α β and N α β ) should be zero, which is only achieved by the new bending model. ∆κ 1 := κ 1 − κ 01 is introduced for simplification.

Pure bending
In the second case, the half tube undergoes pure bending to have the radius r and length l without any change in the arc length (Fig. 2b). The current configuration, parametrized by Eq. (70), is now described by The corresponding tangents and normal vectors then follow as Further, the surface metric and curvature tensor components become With this the stretch and curvature can be calculated as and The resulting stress and moment components for test case 2(b) are listed in Tab. 5. The components of Cauchy stress N α β should be zero, and only our new bending model accurately captures this.
(a) Initial configuration (b) Current configuration Here γ is a constant that defines the extent of twisting. The tangent vectors and normal follow from Eq. (110) as while the kinematic quantities become and We emphasize that the principal curvatures in the current configuration, κ * i , which follow from Eq. (8) as will be different from κ i .
Helfrich apH new = 0 and M 22 0 = 0, which is only achieved by the Koiter and the proposed new model (when c 12 = 0). The undesired membrane stress contributions from all bending models are much smaller than the stress contributions from a membrane model, even at large deformation.
From the Fig. 3, we know that M 22 0 should be non-zero, as the fibers along ξ 2 = φ are being bent. In contrast, M 11 0 should be equal to zero, as there is no change in curvature of the fibers along ξ 1 = θ. As shown in Tab. 6, only the Koiter model has M 11 0 = 0 and M 22 0 = 0. For the proposed new model, M 22 0 = 0 and if we consider Λ = 0, which is also used for all the other bending models, M 11 0 will be zero according to Eq. (69). As torsion, in contrast to the previous test cases, causes in-plane deformation, we need to compare any stresses coming from the bending model to those of the membrane part of the material model in order to further asses the accuracy of the model. For this, we consider the Neo-Hookean membrane model defined in Eq. (27). Using Eq. (28), the Neo-Hookean's stress contribution will be Let Λ = 0, γ = 1 and slenderness ratio R/T = 100. Then, using Eq. (69), (122) and τ αβ new in Tab. 6, the maximum absolute value of the effective membrane stress component of the new bending model (τ 11 new ) and the Neo-Hookean membrane model (τ 11 mNH ) can be compared. We find The membrane stress contribution coming from the new material model is thus negligible (10 5 times smaller) in comparison to the stress from the Neo-Hookean membrane model. The Eq. (123) also holds true for the new bending model's Cauchy stress N αβ new . The Canham bending model's membrane and Cauchy stress contributions are also of the same order as the new model, as Tab. 6 shows. The stress contribution from the apH model is equal to Eq. (122), which is acceptable as it is from its membrane part Eq. (59.1). Additionally, the Cauchy stress of the Koiter bending model (N αβ bKo ) will also be 10 5 times smaller than the Cauchy stress of the Koiter membrane model (N αβ mKo = τ αβ mKo ) according to Eq. (26).

Summary
The preceding results are categorized as passing, failing, or partial passing depending on the results of M α β and N α β . This categorization is shown in Tab. 7. As seen, only the proposed new bending model passes all test cases. In the table, an additional row is added to report if the models allow for initially stress-free curved surfaces. This is not the case for the Canham and Helfrich model. As a result of not being initially stress-free, the Canham and Helfrich models also fail to provide satisfying results for M α β or N α β in any of the elementary test cases where the body undergoes large deformation. The Koiter and apH models perform similarly well for all the test cases apart from torsion. Both the models pass test case 1(a), but only pass test case 1(b), when κ 2 is small, such that errors in N α β become negligible. They also only pass test case 2(a), when ∆ κ 1 is negligible, due to errors in M α β . In test case 2(b), the Koiter and apH models fail due to errors in both M α β and N α β . Test case 3 is passed by the Koiter, but not the apH model due to incorrect bending moments. Koiter

Simply supported plate under pressure loading
A simply supported square plate of length L = 12 L 0 and thickness T = 0.375 L 0 under the sinusoidal pressure p(x, y) = sin π x L sin π y L E 0 , is analysed. The plate is considered to have Young's Modulus E = 480 E 0 and Poisson's ratio ν = 0.38. Tab. 8 shows the list of corresponding parameters for the different bending models. As the problem is symmetric, only 1/4 of the plate is modeled and the symmetry boundary conditions are enforced using the penalty method of Duong et al. (2017) with penalty parameter = 4.8 n q−1 E 0 L 0 . Here, n refers to the number of elements per side and q is the order of shape functions used. The maximum vertical displacement is compared to the analytical solution from Ugural (2009) and the relative error plot is shown in Fig. 4d for the new bending model. All the other bending models show matching error rates and are therefore omitted from the error plot. In this plot, the results obtained with a regular mesh (Fig. 4a) and skew mesh (Fig. 4b) are compared using solid and dashed lines, respectively. The two results are very close, confirming that the new bending model works also when L i is not aligned with A α . Fig. 4c shows the deformed configuration obtained with the skew mesh.

Pinched cylindrical shell -linear case
This example analyzes a cylinder with rigid diaphragms at its ends and subjected to point forces as shown in Fig. 5. The cylinder is of dimension R = 300L 0 , L = 600L 0 and has thickness T = 3L 0 . The rigid diaphragm is realized by fixing the x-and z-degrees-of-freedom of the nodes lying at the cylinder ends. The material parameters are Young's modulus E = 3 E 0 and Poisson's ratio ν = 0.3. The corresponding material parameters for the different bending models are shown in Tab. 9. The magnitude of the two opposing pinching forces is F = 10 −6 E 0 L 2 0 . Owing to symmetry, the problem is solved for only 1/8 of the cylinder as shown in Fig. 5a. All symmetry conditions are enforced using a penalty method with parameter = 6×10 2 n q−1 l E 0 L 2 0 for axial symmetry and = 6 × 10 2 n q−1 t E 0 L 2 0 for circumferential symmetry, where n l and n t are the number of elements in axial and circumferential directions, respectively. The finite element solution is verified by the analytical solution for the displacement beneath the force from Flügge (1962) and Duong et al. (2017). The proposed bending model converges to the accurate solution as shown in Fig. 5b. Identical convergence behavior is observed for the Koiter and apH model. As the Helfrich bending model is not initially stress-free due to the rear term in Eq. (57.2), highly inaccurate results are obtained for coarse meshes (Fig. 5c).

Bending model
Parameter set 3D linear elasticity

Pinched cylindrical shell -nonlinear case
We consider the same problem as in the previous example (see Fig. 5a), but now undergoing large deformations. The radius, length, and thickness of the cylinder now is R = 100 L 0 , L = 200 L 0 , and T = 1.0 L 0 , respectively. The cylinder is considered to have a Young's modulus E = 30 E 0 , Poisson's ratio ν = 0.3. The corresponding bending parameters following from this are tabulated in Tab. 10. Fig. 6a shows the deformed configuration obtained with the point force of F = 12 E 0 L 2 0 applied in 40 loading steps. By employing Lagrange multiplierbased symmetric boundary conditions as in Duong et al. (2017), only 1/8 of the cylinder is used in the FE computations. Fig. 6c presents the force vs. displacement curve (measured at points A and B shown in Fig. 5a) (58) and (69). The following figures show the membrane stress N α β = N αγ a γβ (Fig. 7) and bending stress M α β = M αγ 0 a γβ /J (Fig. 8), obtained from post-processing based on L2-projection with lumped mass matrix. For N α β , only the result of the new bending model is shown as the two other models give similar results. However, for M α β the Helfrich model shows distinct differences to both the Koiter and new model. This again shows that the Helfrich model is not able to capture initially stress-free shells properly.

Angioplasty simulation
The last example considers contact between an expanding balloon and a tube, which mimics the angioplasty procedure used to clear blocked coronary arteries. The tube has the dimensions R t ×L t ×T t = 5 mm× 30 mm × 0.5 mm and represents a portion of an artery. The balloon within the artery is initially spherical and has the initial radius R b = R t and volume V 0 = 4π/3R 3 b . The tube is modeled using the incompressible Neo-Hookean membrane model from Eq. (59.1) together with either the Koiter or the proposed new bending model using the parameters µ t = 8 N/m (corresponding to E = 48 kPa and ν = 0.5 according to Eq. (50)). The balloon is also modeled by Eq. (59.1) but with µ b = 5 µ t . The tube and balloon are discretized using n el = 60 × 55 and n el = 64 × 32 quadratic NURBS elements, respectively. The tube and balloon are supported by fixing all normal dofs on the x-, y-and z-symmetry planes. This leaves the tube free to contract longitudinally. The balloon is inflated up to the volume V = 3V 0 . Frictionless contact is enforced using the two-half-pass penalty contact formulation of Sauer and De Lorenzis (2013) with penalty parameter = 5 E/mm. A similar angioplasty example was studied in Sauer (2017, 2018) using the material model of Gasser et al. (2006).  Fig. 10a, the deformed configurations are colored by the circumferential stretch λ 1 . Though this distribution looks similar, the balloon has slightly higher longitudinal deformation for the Koiter bending model than for the new model, which indicates that the Koiter model is stiffer. This is confirmed by the larger contact pressure p c and circumferential bending stress component M 1 1 of the Koiter model shown in Fig. 10b and c, respectively. The longitudinal bending stress component M 2 2 , on the other hand, is again similar for both models as Fig. 10d shows. The large differences between the contact pressure p c and moment component M 1 1 can also be seen in Fig. 11. At V = 3V 0 , p c is almost twice as large and M 1 1 almost three times as large for the Koiter than for then new bending model. Up to V ≈ 1.2V 0 , however, the differences are small, illustrating once more the equivalency of the two models for small deformations.
Contrary to the previous examples, the angioplasty example thus shows that there are major differences between the Koiter bending model and the proposed new bending model. Those are due to the circumferential stretch in the tube caused by the expanding balloon. In the Koiter model this stretch contributes to the bending moments and stiffness, while it does not in the proposed new model. The effect of bending and stretching can therefore be properly separated in the new model. This was also seen in test case 2 (Sec. 4.2). Since the Koiter bending model reacts to stretches it can be expected to overestimate the bending moments.

Conclusion
This work proposes a new bending model for Kirchhoff-Love shells based on the direct surface approach. The proposed bending model eliminates the spurious influence of membrane strains on bending at large deformation. It is objective, can handle initially curved surfaces and large deformations.
The new bending model is introduced in Sec. 3 along with several existing bending models. The relation between the material parameters of the different bending models is provided, and it is shown that all the bending models become equivalent for initially planar shells at small deformations. In contrast to existing bending models, the proposed new model passes an essential set of five elementary bending test cases as was seen in Sec. 4. Test case 1(a) (Sec. 4.1) shows that the Canham and Helfrich models are not initially stress-free. The Helfrich model will be initially stress-free only when one of its bending parameter (k) is zero. Test case 1(b) considers a final configuration similar to that of a rigid rotation but obtained by counter bending. This test case highlights the anisotropic nature of the new bending model in addition to providing accurate stresses. The Koiter and apH models only partially pass this test due to errors in the Cauchy stress. The Helfrich model fails due to incorrect bending stress and the Canham model fails due to incorrect Cauchy as well as bending stresses. In test case 2(a) (Sec 4.2), a cylinder is inflated, which is a pure membrane action, and hence stresses shouldn't arise from the bending models. The Koiter and apH models partially pass this test with small errors in the bending stresses. The Canham and Helfrich models fail again due to incorrect bending and Cauchy stresses. In test case 2(b), pure bending is considered and all the bending models except the proposed model fail to give zero Cauchy stresses. In test case 3 (Sec 4.3), torsion is considered for which only the Koiter and new bending model are able to provide accurate bending stresses to pass the test.
These observations are confirmed by the five numerical test cases in Sec. 5. In the simply supported linear plate problem (Sec. 5.1), all the bending models considered are identical. This problem is also solved with skew meshes to exhibit that the proposed model gives accurate results also when the principal curvature directions are not aligned with the curvilinear coordinates. The Helfrich model fails in the second linear problem -the pinched cylindrical shell (Sec. 5.2) -due to non-zero initial bending stress. The Helfrich model also deviates from the reference solution for cylinder spreading (Sec. 5.4) for the same reason. Though the Helfrich model gives accurate force vs. displacement for the nonlinear pinching problem (Sec. 5.3), its bending stress components (M α β ) are different from those obtained for the Koiter and the new bending model. Finally in the angioplasty example, the difference between the Koiter and the new model is illustrated. The example shows that bending and stretching can be properly separated in the new model, which is not the case for the Koiter model where it results in much larger bending moments. The proposed new model is able to circumvent problems in existing bending models by using stretch-invariant bending quantities. It serves as a bending model with an union of desired features for shell models in the direct surface approach. As the model is based on the principal curvature directions and has four different material parameters, it provides flexibility in modeling new materials. In future, it would be interesting to combine the new bending model with anisotropic membrane models based on the Mooney-Rivlin (Mooney, 1940;Rivlin and Saunders, 1951) and Gasser-Ogden-Holzapfel model for biological tissues (Gasser et al., 2006).
The tensor components corresponding to the bending models of Koiter, projected Neo-Hooke, and Canham and Helfrich can be found in Duong et al. (2017), Roohbakhshan and Sauer (2017) and , respectively.

B Efficient FE implementation
If the initial principal curvature directions align with surface tangents (A α ), then the quantity defined in Eq. (66) will be The stretches along those direction will then simply be, Plugging this into the quantities of Eq. (68) where A 1 := A 11 a 11 , A 2 := A 22 a 22 .
Using this simplification, the stress and moment can be directly calculated as τ 22 = − c 12 k 1 + c 2 k 2 + c 3 √ λ 1 κ 12 k 12 2 √ λ 2 τ 12 = τ 21 = 0 , M 11 0 = (c 1 k 1 + c 12 k 2 ) M 22 0 = (c 12 k 1 + c 2 k 2 ) As described in Duong et al. (2017), for the efficient computation of the FE stiffness matrices we can exploit the symmetries to rearrange the fourth order tensor, f αβγδ as Based on the simplification in Eq (131), this rearrangement simplifies to Similar rearrangement can be applied to c αβγδ , d αβγδ and e αβγδ . Further, we define the auxiliary terms L a αβ := N T ,α a β , L n α := N T ,α n , G n αβ :=Ñ T ;αβ n , where each term is an array of size (3n × 1), with n being the number of control points per element. This is further reorganized aŝ L a = [L a 11 , L a 22 , L a 12 + L a 21 ] , G n = [G n 11 , G n 22 , G n 12 + G n 21 ] .
We can then rewrite the equations for force, material stiffness, and geometric stiffness, k e M1 = − Ω e 0 b M a 11 L n 1 L nT 1 + a 22 L n 2 L nT 2 + a 12 L n 1 L nT 2 + L n 2 L nT 1 dA , (152)

C Bending moduli extraction
In order to calculate the bending moduli of the new bending model we directly compare the values of f αβγδ in (40) Choosing a parametrization with A 12 = A 21 = 0 and noting that f αβγδ has major and minor symmetries, the tensor can be represented only by the following elements  1 2 A 2 2 ) .
The corresponding fourth order tensor components for the new model are of the form The way L α i is calculated in the numerical examples is by using for a particular choice of parametrization. Thus only L 1 1 and L 2 2 will be non-zero. So by considering λ i = 1,