A hyperelastic extended Kirchhoff–Love shell model with out-of-plane normal stress: I. Out-of-plane deformation

This is the first part of a two-part article on a hyperelastic extended Kirchhoff–Love shell model with out-of-plane normal stress. We present the derivation of the new model, with focus on the mechanics of the out-of-plane deformation. Accounting for the out-of-plane normal stress distribution in the out-of-plane direction affects the accuracy in calculating the deformed-configuration out-of-plane position, and consequently the nonlinear response of the shell. The improvement is beyond what we get from accounting for the out-of-plane deformation mapping. By accounting for the out-of-plane normal stress, the traction acting on the shell can be specified on the upper and lower surfaces separately. With that, the new model is free from the “midsurface” location in terms of specifying the traction. We also present derivations related to the variation of the kinetic energy and the form of specifying the traction and moment acting on the upper and lower surfaces and along the edges. We present test computations for unidirectional plate bending, plate saddle deformation, and pressurized cylindrical and spherical shells. We use the neo-Hookean and Fung’s material models, for the compressible- and incompressible-material cases, and with the out-of-plane normal stress and without, which is the plane-stress case.

We are introducing a hyperelastic extended Kirchhoff-Love shell model with out-of-plane normal stress. In the first part of a two-part article, we present the derivation of the model, with focus on the mechanics of the out-of-plane deformation. To determine the out-of-plane stress, we solve the linear-momentum-balance equation in the out-of-plane direction. Accounting for the out-of-plane normal stress distribution in the out-of-plane direction affects the accuracy in calculating the deformed-configuration out-of-plane position, and consequently the nonlinear response of the shell. The improvement is beyond what was achieved with the new model's precursor [5] by accounting for the out-of-plane deformation mapping.
A good number of shell models were presented earlier in the finite element context (see, for example, [25][26][27][28][29][30][31]), with significant effort in bending representation. The model in [28] is based on a mixed formulation. The model in [31] is based on a discontinuous-Galerkin type approximation to weakly enforce C 1 continuity. The model in [30] is a TUBA family element, which has displacement derivatives as unknowns to attain C 1 continuity in the displacement. The model we are introducing here is similar to the model in [25], which uses only one parameter to represent the outof-plane deformation. Most of the other shell formulations, including some based on the Reissner-Mindlin theory, use the plane-stress assumption. The models in [26,27], based on the Reissner-Mindlin theory, are, however, without the plane-stress assumption, in the finite element context.
Continuing what was started with its precursor, the model introduced here is extending the range of applicability of the 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 material model has different versions. In the version used in [13], the first invariant of the Cauchy-Green deformation tensor appears in a squared form. In the version used in this article, as in [5], it appears without being squared, and this version has been used in a number of arterial FSI computations [32][33][34][35][36][37][38][39] with the continuum model.
By accounting for the out-of-plane normal stress, the traction acting on the shell can be specified on the upper and lower surfaces separately. This enables not only more accuracy in the linear-momentum balance in the out-of-plane direction, but also representation of the moment the shear tractions on the upper and lower surfaces generate around the midsurface. With separate out-of-plane tractions on the upper and lower surfaces, for example, we can accurately model cases that might have nonzero net force even when those out-of-plane tractions have equal magnitudes and opposite directions. The net force would be nonzero because the upper and lower surfaces would have different areas due to the curvature. To accurately account for the moment generated by the separate shear tractions on the upper and lower surfaces, we improve the rotational kinematics in the model.
We note that accounting for the out-of-plane stress improves the out-of-plane deformation mapping also in cases with no traction on the upper or lower surfaces and no body force. Those would be the cases when the shell deformation is driven by the displacements and slopes specified along the edges of the shell.
We also would like to note that the level of accuracy we are striving for in representing the tractions on the upper and lower surfaces would be meaningful in an FSI computation only if the flow solution method can deliver those tractions with a comparable level of accuracy. That level of flow solution accuracy, especially in representing the shear stress, requires moving-mesh methods [9], where the high mesh resolution near solid surfaces follows the fluid-solid interface as it moves. That is now possible even in flow computations with actual contact between solid surfaces or some other topology change. The Space-Time Topology Change method [40] enabled that. We can both represent the actual contact and have high-fidelity, moving-mesh flow solution near the solid surfaces.
Our test computations are based on solving the linearmomentum-balance equation in the out-of-plane direction with finite element discretization. The computed problems are unidirectional plate bending, plate saddle deformation, and pressurized cylindrical and spherical shells. We use the neo-Hookean and Fung's material models, for the compressible-and incompressible-material cases, and with the out-of-plane normal stress and without, which is the plane-stress case.
In Sect. 2, we provide the definitions and concepts used in the shell model, including the notations and main assumptions. In Sect. 3, we derive the weak form from the virtual-work principle and express the corresponding strong form. In Sect. 4, we describe the solution technique for the out-of-plane deformation, which is meant to be used only for the test computations in Sect. 5. The concluding remarks are given in Sect. 6. In the Appendix, we provide some supplemental derivations and the constitutive models.

Kinematics
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. Here, we assume n sd = 3. The subscript t indicates the timedependence of the domain.
We split the domain as Ω t = Γ t × (h th ) t , where Γ t represents the midsurface, which is parametrized by n pd = n sd −1, with n pd being the number of parametric dimensions. The remaining parametric direction is the out-of-plane direction, which we will explain a little later in this section. The symbol (h th ) t will represent both the shell thickness and the domain in the out-of-plane direction, depending on the context. With the position x ∈ Γ t , we define a natural coordinate system: where α = 1, . . . , n pd and ξ α represents the parametric space. Figure 1 shows a schematic shell domain. We note that this parametric space is only for representing the neighborhood around a point. With a range −1 ≤ ξ α ≤ 1 and without loss of generality, it can be seen as a parent finite element domain. The out-of-plane direction is Fig. 1 A schematic shell domain Ω t at current configuration and the midsurface Γ t . The black-framed region is for showing the parametric space used in representing the neighborhood around a point. With a range −1 ≤ ξ α ≤ 1, the parametric space can be seen as a parent finite element domain. The basis vectors of the natural coordinate system are also shown in the figure The components of the metric tensor are and this is known as the first fundamental form. We also work with the contravariant components of the metric tensor g αβ and the contravariant basis vectors g α (see Appendix A for their relationship to g αβ and g α ).
Remark 1 The symbols κ αβ and ω αβ are for notational convenience.
We note that the curvature tensor is defined by κ κ κ = κ αβ g α g β , (15) and the square of the curvature tensor is the tensor version of ω αβ : We also use the dual basis system with g αβ and g α (see Appendix A for their relationship to g αβ and g α ).
We now provide similar definitions and derivations for the undeformed configuration X ∈ Γ 0 . We start with the basis vectors where ξ α 0 = ξ α , and We also again work with the contravariant components of the metric tensor G αβ and the contravariant basis vectors G α (see Appendix A for their relationship to G αβ and G α ).
A position X ∈ Ω 0 is expressed as where . Similar to what we had for the current configuration, along ξ 3 0 , the basis vectors are represented as where The lines between Eqs. (23) and (24) are the deformedconfiguration counterpart of the lines between Eqs. (8) and (9). The metric tensor components in 3D space are where We also again use the dual basis system with G αβ and G α (see Appendix A for their relationship to G αβ and G α ).

Deformation gradient tensor
The deformation gradient tensor F is defined from which implies Using the relationship ξ α = ξ α 0 and introducing we can write

Cauchy-Green deformation tensor
The Cauchy-Green deformation tensor is defined as which is and the determinant of C gives the square of J = det F:

Green-Lagrange strain tensor
The Green-Lagrange strain tensor is defined as where I is the identity tensor, and one of the ways to express it is We can then write E as The covariant components of the in-plane strain tensor are The symbol ε αβ here represents what is typically identified as the membrane strain. That is only part of the membrane strain in our model. The normal component associated with the out-of-plane direction is

Remark 2
The term of Eq. (46) that is quadratic in ξ 3 and ξ 3 0 was omitted in [5] as a higher-order term. We retain that here.

The strain-energy density function
We express the strain-energy density function as The second Piola-Kirchhoff stress tensor, is obtained from and from that and Eq. (41), The Cauchy stress tensor, defined in the current configuration, can be obtained from

Variational formulation
The variations of the internal and kinetic energies are expressed as where ρ 0 is the density at the undeformed configuration. The principle of virtual work can then be written as where δW ext is the external-force virtual work, explained more later.

Remark 3
We note that in [9] (Section 1.2.2), −δW int is used in place of δU .

Remark 4
We also note that in [9] (Section 1.2.2), the variation of the kinetic energy, expressed in terms of the acceleration, is part of the external-force virtual work. Here we keep that separate. With this separation, Eq. (55) can be seen as d'Alembert's principle.

Remark 5
With the displacement expressed as y = x − X, we can do the substitutionsẋ =ẏ andẍ =ÿ.

Admissible variations
The variation of x is obtained by taking the variation of Eq. : The variation δx represents the virtual displacement of the midsurface, δn represents the rotation of the midsurface (see Appendix B.2), and δξ 3 represents the virtual displacement in the out-of-plane direction, relative to the midsurface. From the variation of Eqs. (8) and (32), we obtain The variation of E can be obtained by taking the variation of Eqs. (46) and (47):

Remark 6
In [5], only the first two terms of Eq. (59) were present.

Variation of the strain energy
In one of the ways to express the variation of the strain energy, we start with Because S is a symmetric tensor, we simply write Substituting Eq. (33), its variation, and Eq. (49) into this, we get Substituting Eq. (57) into that, we obtain δE : S = δg α · S αβ g β + δn ,α · S αβ g β ξ 3 We perform the integration in Eq. (53) with the integrand given by Eq. (73): We define p α as its integration along (h th ) 0 aŝ and its first moment aŝ With these definitions, we write Eq. (74) in a simpler form: We also express δU in an alternative form by using Eqs.
where the zeroth, first, and second moments of the contravariant components of S are given aŝ

Remark 9
The fourth, fifth, and sixth integrals in Eq. (79) are the contributions to the virtual work in the out-of-plane direction, which will be explained more in Sect. 3.5.

Variation of the kinetic energy
Performing the integration in Eq. (54), we get Substituting Eq. (56) into this, we obtain Remark 10 We note that where and n ×m 1 represents the time derivative of the angular momentum. The derivation of Eq. (87) is given in Appendix B.3. This rearrangement may make it easier to see the fundamental mechanics in the term.

Remark 11
The second integral in Eq. (84) is related to the time derivative of the angular momentum. It was included in [41], but we did not see other papers that explicitly mention the inclusion of the term.
The acceleration is written as Remark 12 We omit the time derivatives of ξ 3 : Otherwise they need to be stored. We see this approximation as a quasi-steady-state assumption.
With that, approximated forms of Eqs. (85) and (86) are written aŝ wherê In the third term of Eq. (84), the integration along (h th ) 0 can be written, with Eq. (91), as where which can be obtained from n ·ṅ = 0.

Remark 13
The second term in Eq. (98) represents the effect of the centripetal acceleration.

The external virtual work
We separate δW ext into three parts: where Here S t (see Fig. 2) is the edge line that traverses the edge surface as ξ 3 varies from (ξ 3 ) − to (ξ 3 ) + . The symbol h represents the traction, for all surfaces, including the edge surfaces.
Substituting Eq. (56) into Eq. (105), we obtain and the normal tractions aŝ We note that we can seep as the pressure brought to the midsurface by using the area ratios.

Remark 15
In Eq. (108), there are three integrals. In [5], only the first integral was considered, and even in that, the evaluations were not based on the actual upper and lower surfaces as how it was done with Eqs. (109) and (110). The tractions acting on the upper and lower surfaces were treated as if they were acting on the midsurface. Not being able to use the correct surfaces was a consequence of the plane-stress assumption.

Remark 16
The terms in the second line of Eq. (108) can be written as δn ·ĥ − 1 = δr · (n ×ĥ − 1 ) and δn ·ĥ + 1 = δr · (n ×ĥ + 1 ) (see Appendix B.3). We note that n×ĥ − 1 and n×ĥ + 1 represent the moment. Again these rearrangements may make it easier to see the fundamental mechanics in the terms. In a typical shell formulation, the distance between the upper and lower surfaces is not taken into account. Therefore, the tractions on those surfaces, even when they have shear components, do not produce moment.
For the representation at the edges, we introduce a unit vector T along S 0 (see Fig. 3). The contravariant components are denoted as We map T from the midsurface to the ξ 3 0 surface with edge S 0 , and make a unit vector from that: Using Eqs. (24), (26) and (27), T can be expressed as Fig. 3 The cut piece in Fig. 2 (middle) with the unit vectors N, T, and T and its contravariant components are denoted as We note that the numerator of Eq. (119) is coming from the numerator Eq. (116). For notational convenience, we define From Eqs. (121) to (122), we used Eq. (119). We use the symbol L 0 to denote L 0 at the midsurface, obtained by setting ξ 3 0 = 0: This can be seen as the length counterpart of A 0 given in Eq. (40). Using Eqs. (119), (122) and (123), we get Similarly, we introduce unit vectors t and t along S t and S t (see Fig. 4). The stretch λ T along S 0 is obtained from

Fig. 4
The cut piece in Fig. 2 (bottom) with the unit vectors n, t, and t Taking the inner product of g α with both sides of Eq. (126), we obtain Because t is a unit vector, g αβ t α t β = 1 can be used to calculate λ T as Similar to what we did for the undeformed configuration, we define which is the length counterpart of A given in Eq. (39). Using Eqs. (129), (127) and (121), we get and at the midsurface, We now define a unit outward normal vector for S 0 : and for S 0 : (133) Figure 5 shows the unit vectors B and B and the rest of the unit vectors in the undeformed configuration.
Then, corresponding to B, we introduce a unit vector b in the deformed configuration, and the stretch in the direction of B is obtained from With Eq. (134) and b = 1, the stretch can be written as We use the symbols b and λ B to denote b and λ B at midsurface: and similarly write the stretch λ B as Remark 17 We note that, unlike B, the unit vector b is not necessarily normal to S t , and therefore t·b does not have to be zero. To overcome the difficulty with using non-orthogonal basis vectors, we introduce the dual basis vectors: We refer to Appendix C.1 for more on these. Figure 6 shows the vectors t and b , the unit vector b, and some of earlier-defined unit vectors in the deformed configuration.
With the ratio of the lengths along S t and S t being L/L, we change the order of the integration in Eq. (106): With the ratio of the lengths along S t and S 0 being L/L 0 and with Eq. (130), we transform the first integration from the Fig. 6 The cut piece in Fig. 4 with the unit vector b and the dual basis vectors t and b current to undeformed configuration: Using Eq. (32), we transform the second integration from the current to undeformed configuration: Substituting Eq. (56) into this, we get We defineĥ e 0 aŝ and its first moment aŝ With these definitions, we write Eq. (143) in a simpler form: We omit the variation of ξ 3 :

Remark 18
The reasons for omitting the variation of ξ 3 in Eq. (146) are, (i) the out-of-plane shear stress is not meant to be represented in this model, and (ii) we already omitted the time derivatives of ξ 3 (see Remark 12).
(55), exclude the variations associated with the midsurface, and obtain We substitute Eq. (58) into this and consider only the out-ofplane variations at a given midsurface point: Applying integration by parts to the last integral on the lefthand side, we obtain We substitute Eqs. (99) and (100) into this, and the equation holds for all admissible δξ 3 . From that, we obtain the strong form: and Integrating Eq. (151) from (ξ 3 0 ) − to (ξ 3 0 ) + and using Eqs. (152) and (153), we obtain We rearrange this equation by substituting forp − andp + from Eqs. (113) and (114)

The virtual work for the midsurface deformation and the corresponding strong form
wherê Equation (156) is our weak form for the midsurface representation. To understand it better, we substitute into it, apply integration by parts to the first and third integrals, and obtain We note that, as shown in Fig. 5, the unit outward normal vector for S 0 is B. With δn involving derivatives, we next apply integration by parts to the second and fourth lines of Eq. (161). We do that by using two expressions derived in Appendix C, given by Eqs. (306) and (316). After that, separating δx into its normal and tangential components, we obtain We note thatp α 0 andp 1 have only in-plane components; n · p α 0 = 0 and n ·p α 1 = 0 (see Eqs. (75), (76) and (77)). With that, we get Using Eq. (164), we rewrite Eq. (162) as From that, the strong form in Γ 0 can be written as We note that we split the strong form into its parts associated with the normal and in-plane directions. Along the boundary S 0 , we have Here Eqs. (168) and (169) are the parts associated with the normal and in-plane directions.
Remark 20 We can show that the first two terms of Eq. (166) constitute an alternative form of Eq. (155) and therefore give zero. For that, we rewritep α 0 by using Eqs. (76), (75) and (8): take the inner product of both sides of that with −n ,α , and obtain We rearrange that by using Eqs. (264) and (13): Recognizing the terms on the right-hand side from Eqs.
With that and Eq. (157), Eq. (155) can be rewritten as Combining that with Eq. (166), we obtain The equation above represents the rotational balance around the midsurface.

Remark 21
and recognizing the right-hand side as the moment acting along S 0 . Using Eq. (139), Eq. (170) can be written as which can be rearranged as This means that Eq. (170) represents the tangential component of the rotation boundary condition along S 0 . We note that we cannot have a component in the n direction because n ×ĥ e 1 · n = 0, and we cannot have a component in the b direction because of the Kirchhoff-Love assumption.

Stress and strain distributions in the out-of-plane direction
We solve for the stress and strain distributions along the outof-plane direction at a given midsurface point. We use the symbol y 3 to denote the displacement in the out-of-plane direction: and write The point on the midsurface is like an integration point in a typical finite element formulation. For simplicity, we assume steady state and exclude the body force. In the test computations, which will be reported in Sect. 5, we found out that determining in advance the midsurface variables needed in the equation associated with the out-of-plane direction is not a good path. Therefore, in this section, we treat also ε αβ and κ αβ as unknowns, with total six components. Then, in addition to Eq. (148), we use the first three integrals of δU given by Eq. (79), which are associated with the midsurface deformation, with δε αβ and δκ αβ generating six additional equations. We note that other virtual-work terms in Eq. (55) play no role in the midsurface deformation. That is because δT = 0 from assuming steady state, and in δW ext , δW exbody = 0 from excluding the body force, δW exsurf = 0 from being for the upper and lower surfaces, and δW exedge = 0 from being for the edge surfaces. We also note that δω αβ appearing in the third integral of Eq. (79) can be represented in terms of δε αβ and δκ αβ , as described in Remark 7. Some of the six additional equations will be removed for consistency with the constraints coming from the boundary conditions or geometry, depending on the problem setup.

Remark 22
In the actual computations, which will be presented in the second part of this two-part article, we plan to use a different way of numerically solving the strong form of the equations in the out-of-plane direction. That will involve change of unknown variable. We expect that coupling that solution process with solving for the midsurface deformation will not be computationally burdensome.
At this point, we will change the symbol we use for representing the variation associated with the virtual work, and call it "δ a ." That is because in solving the nonlinear equations with the Newton-Raphson method, we will need to take another variation for the linearization in the iterations, and that variation will be denoted by the symbol δ b . In other words, the linearization for the Newton-Raphson iterations will be in terms of δ b . Using the test functions of δ a ξ 3 , δ a ε αβ , and δ a κ αβ , we form the equation system as where N (δ a ξ 3 , δ a ε αβ , δ a κ αβ , y 3 , ε γ δ , κ γ δ ) We note that the undeformed shape is represented by G αβ and K αβ , and the current configuration by ε αβ and κ αβ .

Remark 23
The idea of representing the undeformed shape in terms of the metric tensor is from [42], and in that context the configuration is called the integration-point based zerostress state (IPBZSS). In the shell model, full representation of the IPBZSS requires also the curvature tensor.
We also note that which follows from Eq. (45), and which follows from Eq. (254) in Appendix A.
In describing the linearization for the Newton-Raphson iterations, we will first derive the variations needed for that. We start with the variations of S αβ and S 33 : where The elastic moduli C αβγ δ , C αβ33 , C 33γ δ , and C 3333 can be derived from ϕ. In Eqs. (59) and (60), using δ b in place of δ, recognizing that δξ 3 and δ y 3 are in interchangeable, and reordering the terms, we obtain and Taking the variation of Eq. (181), we obtain From the variations of Eqs. (80), (81) and (82), we obtain Now that we have all the variations we need, we write δ b N by taking the variation of Eq. (183). After taking the variation, we substitute into it Eqs.

Remark 25
The choices made between using different combinations of variable in Eq. (203) were driven by having visible symmetry in δ a and δ b and reducing the number of terms.

Remark 26
The term δ a δ b ω αβ is not zero, because, as explained in Remark 23, δ a ω αβ is not an independent variation. We replace δ with δ a in Eq. (68), take its δ b variation, multiply that with 1 2Ŝ αβ ω , and obtain

Remark 27
When the traction conditions on the upper and lower surfaces are coming from the pressure on those surfaces, with the conditions expressed as We note that A depends on both the midsurface deformation and the displacement in the out-of-plane direction. We derive the expression for δ b A based on that dependence and use it in the computations.

Test problems
With the method described in Sect. 4, for a set of test problems, we determine the displacement in the out-of-plane direction. In all cases, we use a uniform finite element mesh made of 200 elements with linear functions 1 . The lower or inner surface serves 2 as the "midsurface." Together with the out-of-plane displacement, the midsurface variables ε αβ and κ αβ are solved for, with different boundary or geometry conditions depending on the problem setup. For comparison, we include in the tests a plane-stress shell model from [5]. The constitutive models, provided in Appendix D, are neo-Hookean and Fung's models with compressible and incompressible materials. The stress-related variables are scaled by the shear modulus μ 0 at the undeformed configuration. The Poisson's ratio ν is defined based on the bulk modulus at the undeformed configuration, with ν = 0.50 indicating incompressible material. We use separate models for compressible and incompressible materials, but the compressible-material model converges to the incompressible-material model in the limit ν → 0.50 (see [5]). The scaled stress components we report will be from the Cauchy stress, and the indices will correspond to the coordinate frame defined separately in each problem, but the index 3 will always imply the out-of-plane direction: For the length scale, we use the undeformed shell thickness (h th ) 0 = (ξ 3 0 ) + − (ξ 3 0 ) − . We note that, with the shell model and method sections behind us, now we are using the symbol (h th ) 0 for the shell thickness, as we stated in Sect. 2.1 that we would. We also note that the current thickness is expressed We use the principal curvatures at the midsurface to describe the individual problems. The principal curvatures κ 1 and κ 2 are obtained from the generalized eigenvalue problem where t β represents the contravariant components of the principal curvature direction. From Eq. (212), we can get the principal curvatures as 2 We note that, as can be inferred from the overall formulation, in the shell model proposed here, the location of the midsurface does not matter. However, it does in the model proposed in [5]. We will be comparing our results to those obtained with that model. This is the reason why the lower or inner surface serves as the midsurface.
Their nondimensional versions are

Unidirectional plate bending
The bending is only in ξ 1 direction (see Fig. 7). No traction is applied on the upper or lower surface. The bending moment M is applied along the midsurface edges. We select G α to be unit orthogonal basis vectors, and consequently G 11 = G 22 = 1, and G 12 = 0. We note that, with this selection, we no longer have −1 ≤ ξ α ≤ 1. In this problem, because it is a plate, K 11 = K 22 = K 12 = 0. We constrain the midsurface deformation with ε 22 = ε 12 = 0 and κ 22 = κ 12 = 0. With these conditions, from Eq. (213) we get and its nondimensional version is For any value of κ * 1 , we set κ 11 based on that relationship, noting that g 11 is a function of ε 11 . We solve for ε 11 together with the displacement in the out-of-plane direction. After that, we calculate the required moment around ξ 3 = 0: (220) In this problem, L L = 1 and g 21 = 0, which follow from the constraints listed earlier. Figures 8 and 9 show, for the neo-Hookean and Fung's materials, M * as a function of κ * 1 . For both materials, the model with the plane-stress assumption (σ 33 = 0) is behaving stiffer than the model with the out-of-plane normal stress (σ 33 = 0). For the Fung's material, the difference becomes smaller at higher ν values.
We report, for κ * 1 = 0.5 and 1, ξ 3 0 profiles of different variables. We report λ 3 , J , σ * 11 , and σ * 33 in Figs. 10, 11, 12 and 13 for the neo-Hookean material, and in Figs. 14, 15, 16 and 17 for the Fung's material. For the neo-Hookean material, λ 3 obtained with the σ 33 = 0 model is more than it is with the σ 33 = 0 model. Assuming that the σ 33 = 0 model is giving the true solution, σ 33 ≤ 0. This means that the σ 33 = 0 model  is behaving in a fashion that inhibits the compressive stress that should be there. The larger λ 3 makes h th higher, and that could be the source of the stiffer response mentioned when Figs. 8 and 9 were discussed. Overall, the Fung's material behaves similar to how the neo-Hookean material behaves. However, because of the strong nonlinearity, the λ 3 profile is more complex. For example, there is a local maximum close to the inner surface. Also, at high values of ν, σ 11 becomes very high at both the lower and upper surfaces. That is making σ 33 higher.
To have a better understanding of the significance of σ * 33 , we scale it by using, in general, σ * 11 max and σ * 22 max , where "max" denotes the maximum over both ξ 3 0 and ν. To be on the conservative side in not overrating that significance, we scale σ * 33 with σ * MAX12 = max σ * 11 max , σ * 22 max . We note that in this case σ * 22 max = 0. Tables 1 and 2 show, for the neo-Hookean and Fung's materials and for κ * 1 = 0.5 and 1, σ * 11 max , σ * 33 max , and σ * 33 max σ * MAX12 . We see that the  significance, as measured by σ * 33 max σ * MAX12 , is more than 10 % for the neo-Hookean material and more than 5.9 % for the Fung's material.

Plate saddle deformation
Bending moments are applied along the midsurface edges both in ξ 1 and ξ 2 directions to create a saddle deformation. No traction is applied on the upper or lower surface. Again, we select G 11 = G 22 = 1, G 12 = 0, and K 11 = K 22 = K 12 = 0 to represent a plate. We constrain the midsurface deformation with ε 12 = 0 and κ 12 = 0. With these conditions, from Eqs. (213) and (214) we get and their nondimensional versions are For any values of κ * 1 and κ * 2 , we set κ 11 and κ 22 based on these relationships, noting that g 11 is a function of ε 11 and g 22 is a function of ε 22 . We solve for ε 11 and ε 22 together with the displacement in the out-of-plane direction.
The saddle deformation is obtained by specifying the two principal curvatures to be κ * 1 > 0 and κ * 2 < 0. The curvature at the upper surface can be written as and this needs to be finite. That requirement brings the restriction Considering that, we set κ * 2 = − 1 2 . We report, for the two test cases given in Table 3, ξ 3 0 profiles of λ 3 , J , σ * 11 , σ * 22 , and σ * 33 . Figures 18, 19, 20, 21 and 22   show those profiles for the neo-Hookean material, and Figs. 23, 24, 25, 26 and 27 for the Fung's material. Case 2 is closer to having κ * 1 = −κ * 2 at ξ 3 0 / (h th ) 0 = 0.5. This is the reason behind the near symmetry between the σ * , there is more difference between what we get from the σ 33 = 0 and σ 33 = 0 models than between what we get from different ν values. In other words, the shell models being different is more significant than the ν values being different. Tables 4 and 5 show, for the neo-Hookean and Fung's materials and for Case 1 and Case 2, σ * 11 max , σ * 22 max , σ * 33 max , and σ * 33 max σ * MAX12 . This time σ * 33 max σ * MAX12 is more than 16 % for the neo-Hookean material and more than 6.2 % for the Fung's material. These numbers are higher compared to what we had in the unidirectional plate bending.

Pressurized cylindrical and spherical shells
The inner and outer pressures are p I and p O , and the inner and outer radii are R I and R O in the undeformed configuration   Table 6 Pressurized cylindrical and spherical shells. Three methods, characterized by selecting between the σ 33 = 0 and σ 33 = 0 models and whether p O is applied at the outer or inner surface and r I and r O in the deformed configuration (see Fig. 28). We note that R O = R I + (h th ) 0 and r O = r I + h th . We set (h th ) 0 2R I = 0.1, the same value used for the undeformed configurations in [5].
We compute with three different methods, characterized by selecting between the σ 33 = 0 and σ 33 = 0 models and whether p O is applied at the outer or inner surface (see Table  6). In Method 1, by applying both p I and p O at the inner surface, which serves as the midsurface, the difference between the inner and outer surface areas is not taken into account in calculating the net force acting on the shell. That is taken into account in Method 2 by applying p O at the outer surface. Method 3 also takes that into account, but with the σ 33 = 0 model. We measure the deformation by λ 1 . For the cylindrical shell, we set λ 2 = 1. We set p I = 2 p and p O = p, and p is scaled as p * = p μ 0 .

Pressurized cylindrical shell
We select the undeformed configuration as and constrain the midsurface deformation with In representing the mutual dependence between κ 11 and ε 11 based on the constraints above, we select the form κ 11 being a function of ε 11 . We solve for ε 11 together with the displacement in the out-of-plane direction, with the conditions p * O = p * and p * I = 2 p * (see Remark 27 for the way the pressure conditions are specified). Figures 29 and 30 show, for both the neo-Hookean and Fung's materials, p * as a function of λ 1 . We first compare Method 1 and Method 2. With Method 1, because p O acts on the inner surface, the net force in the inflation direction is more. Consequently, in all cases, for a given p * , Method 1 has higher λ 1 . Next we compare Method 2 and Method 3. For the neo-Hookean material, in all cases, for a given λ 1 , Method 2 has higher p * , meaning that it is stiffer than Method 3. For the Fung's material, we do not see that in all cases. When ν = 0.5, Method 2 is stiffer for all values of λ 1 , but at lower ν values, beyond a high-enough λ 1 value, Method 2 switches to being less stiff. We will discuss that more when we report the ξ 3 0 profiles. Remark 28 It is clear from Figs. 29 and 30 that when ν = 0.5, for both the neo-Hookean and Fung's materials, there is no λ 1 value that would make Method 2 less stiff. In other words, if the behavior is incompressible, Method 2, i.e. the σ 33 = 0 model, is stiffer.
We report, for λ 1 = 1.3, ξ 3 0 profiles of λ 3 , J , σ * 11 , and σ * 33 . Figures 31, 32, 33 and 34 show those profiles for the neo- Hookean material, and Figs. 35, 36, 37 and 38 for the Fung's material. Because the only difference between Method 1 and Method 2 is how the net force acting on the midsurface is calculated, when we look at the ξ 0 3 profiles coming from the two methods at a given λ 1 , meaning at the same midsurface deformation, we are already beyond how the net force acting on the midsurface was calculated. The midsurface deformations are the same, and therefore the ξ 0 3 profiles will be the same. For that reason, the ξ 3 0 profiles are reported under the labels "σ 33 = 0" and "σ 33 = 0." When ν = 0.5, the deformation patterns obtained with the σ 33 = 0 and σ 33 = 0 models are the same (see Figs. 31,32,35,36). This is expected because when ν = 0.5, the constraint J = 1 determines the profile for a given midsurface deformation.
For the Neo-Hookean material, as can be seen in Fig. 33, the σ 33 = 0 model yields higher σ * 11 values at all ν values and for the full range of ξ 3 0 . Remembering that these profiles are for a given value of λ 1 , higher σ * 11 means that σ 33 = 0 model stiffer. This is consistent with the observation we made when we discussed Fig. 29.  values for the full range of ξ 3 0 , and therefore is less stiff than the σ 33 = 0 model. This can be explained by the shifting balance between the bulk and shear moduli as λ 1 varies. At ν values not far from 0.5, when λ 1 is not so high, the bulk modulus is dominant, the material behavior is closer to being incompressible, and therefore, from Remark 28, the σ 33 = 0 model is stiffer. However, even at ν values not far from 0.5, due to the exponential form of the constitutive model, when λ 1 is high enough, the shear modulus is dominant and not the bulk modulus, the material behavior is not close enough to being incompressible, and therefore the σ 33 = 0 model is less stiff.

Pressurized spherical shell
We select the undeformed configuration as and constrain the midsurface deformation with In representing the mutual dependence between κ 11 , κ 22 , ε 22 , and ε 11 based on the constraints above, we select the form κ 11 , κ 22 , and ε 22 being functions of ε 11 . We solve for ε 11 together with the displacement in the out-of-plane direction, with the conditions p * O = p * and p * I = 2 p * . Figures 39 and  40 show, for both the neo-Hookean and Fung's materials, p * as a function of λ 1 . Our observations are essentially the same as those we made for the cylindrical shell, except the differences are more pronounced, because, with the same radius, the sphere has overall higher curvature effects.

Concluding remarks
This was the first part of a two-part article on a hyperelastic extended Kirchhoff-Love shell model with out-of-plane normal stress. We have presented the derivation of the new model, with focus on the mechanics of the out-of-plane deformation. To determine the out-of-plane stress, we solve the linear-momentum-balance equation in the out-of-plane direction. Accounting for the out-of-plane normal stress distribution in the out-of-plane direction affects the accuracy in calculating the deformed-configuration out-of-plane position, and consequently the nonlinear response of the shell. The improvement is beyond what was achieved with the new model's precursor [5] by accounting for the out-of-plane deformation mapping. Continuing what was started with the precursor, the new model is extending the range of applicability of the 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.
By accounting for the out-of-plane normal stress, the traction acting on the shell can be specified on the upper and lower surfaces separately. This enables more accuracy in the linear-momentum balance in the out-of-plane direction. For example, we can accurately model cases that might have nonzero net force even when those out-of-plane tractions have equal magnitudes and opposite directions. Accounting for the out-of-plane normal stress also enables more accuracy in representation of the moment the shear tractions on the upper and lower surfaces generate around the midsurface. To accurately account for the moment generated by the separate shear tractions on the upper and lower surfaces, we have also improved the rotational kinematics in the model. Accounting for the out-of-plane stress improves the out-ofplane deformation mapping also in cases with no traction on   the upper or lower surfaces and no body force. Those would be the cases when the shell deformation is driven by the displacements and slopes specified along the shell edges. We presented test computations for unidirectional plate bending, plate saddle deformation, and pressurized cylindrical and spherical shells. We tested the neo-Hookean and Fung's material models, for the compressible and incompressible materials, and with the out-of-plane normal stress and without, which is the plane-stress case. The test computations show that the differences between the shell models with and without the out-of-plane normal stress are more pronounced i) for plate saddle deformation than unidirectional plate bending, ii) when we have higher curvature effects, and iii) for cylindrical of spherical shells at high ν values. We leave it to the reader to judge how significant these differences are.
The level of accuracy we are striving for in representing the tractions on the upper and lower surfaces would be meaningful in an FSI computation only if the flow solution method can deliver those tractions with a comparable level of accuracy. That level of flow solution accuracy, especially in representing the shear stress, requires moving-mesh methods, where the high mesh resolution near solid surfaces follows the fluidsolid interface as it moves. That is now possible even in flow computations with actual contact between solid surfaces or some other topology change. The space-time computational methods introduced in the last decade enable that, as can be seen, for example, in [43,44]. tract FA520921C0010. This work was also supported (first author) in part by Pioneering Research Program for a Waseda Open Innovation Ecosystem (W-SPRING). The mathematical model and computational method parts of the work were also supported (fourth author) in part by ARO Grant W911NF-17-1-0046, Contract W911NF-21-C-0030 and Top Global University Project of Waseda University.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Dual basis system
The description given here for the relationship between a basis system and its dual is applicable to all the basis systems we work with, where we use the symbols g, g, G, and G.
For the covariant basis vectors g α , the components of the metric tensor are We define the contravariant basis vectors as g α · g β = δ α β .
This implies where the brackets denote the matrix formed from its components, and the contravariant components of the metric tensor are given as We note that, with n sd = 3 and n pd = 2, one of the ways to express the unit tensor is where n = g 1 × g 2 g 1 × g 2 .
(256) B Derivative and variation of the normal vector in the shell model

B.2 Variation of the normal vector
We recall the Lagrange's identity (a × b) · (c × d) = a · cb · d − a · db · c.

C.1 Dual basis system for the midsurface edges
We have a set of orthonormal basis vectors: B, T, and N. After the deformation: With the set of basis vectors and stretches, we can express the deformation gradient tensor and its inverse. The deformation gradient tensor can be expressed as The inverse of F can be written as and we have two alternative ways of verifying that. We can show that by remembering that B, T, and N are orthonormal and using Eq. (285), or show that + Γ 0 δx · n q · g β G β ,γ · G γ dΓ . (306)
Taking the inner product of g γ with Eq. (269): δn · g γ = −δg γ · n, recognizing the right-hand side of that in two places in Eq.

D Constitutive models
We test two constitutive models: neo-Hookean and Fung's materials. The strain-energy density functions are where μ is the shear modulus, and D 1 and D 2 = 8.365 are the coefficients of the Fung's material model. The shear modulus at the undeformed configuration is μ 0 = μ for the neo-Hookean material and μ 0 = 2D 1 D 2 for the Fung's material. We determine the bulk modulus from ν and μ 0 as

D.1 Incompressible material
For incompressible material, we use and p can be calculated by using the constraint J = 1 in the equation associated with the out-of-plane direction.
We use the label "ν = 0.5" in reporting the results for the incompressible-material cases, but we use the forms above in the tests.

D.2 Compressible material
For compressible material, we use where and we use β B = −2. This form with β B was introduced in [45].