On a consistent finite-strain shell theory for incompressible hyperelastic materials

Abstract

In this paper, a consistent finite-strain shell theory for incompressible hyperelastic materials is formulated. First, for a shell structure made of an incompressible material, the three-dimensional (3D) governing system is derived through the variational approach, which is composed of the mechanical field equation and the constraint equation. Then, series expansions of the independent variables are conducted about the bottom surface and along the thickness direction of the shell. The recursive relations of the coefficient functions in the series expansions can be derived from the original 3D governing system. Further from the top surface boundary condition, a 2D vector shell equation is obtained, which represents the local force-balance of a shell element. The associated edge boundary conditions are also proposed. It is verified that shell equation system is consistent with the 3D variational formulation. The weak formulation of the shell equation is established for future numerical calculations. To show the validity of the shell theory, the axisymmetric deformations of a spherical and a circular cylindrical shell made of incompressible neo-Hookean materials are studied. By comparing with the exact solutions, it is shown that the asymptotic solutions obtained from the shell theory attain the accuracy of O(h²).

Keywords

Shell theory nonlinear elasticity incompressible hyperelastic materials asymptotic analysis dimension reduction

1. Introduction

The development of shell theories can be dated back to the pioneering work of Love [1] in 1888. During the past 120 years, shell theories have attracted extensive research interests and the literature in this field is extremely plentiful (cf. [2 –6] and the references therein). The present study is on a derived shell theory for constitutively nonlinear materials with incompressible constraint, so reviews on derived shell theories for constitutively linear materials are only selective. The classic shell theories were usually proposed based on an a priori hypothesis. For example, the Love’s shell theory [1] was established based on the Kirchhoff’s assumptions [7] and the thin shell approximation; Reissner [8] formulated a geometrically nonlinear shell theory based on the Kirchhoff–Love hypothesis but without making a small finite deflection assumption; Naghdi [9] derived the stress–strain relations and the appropriate boundary conditions under the condition of small deformation of thin elastic isotropic shells with uniform thickness; Koiter [10] formulated the first approximation theory of thin elastic shells that is consistent with the Euler–Bernoulli hypothesis. These classical shell theories have been applied successfully in engineering problems. However, their consistencies with the three-dimensional (3D) theory, as well as their ranges of applicability, need to be verified more strictly through mathematical approaches.

Gol’denveizer [11] first developed an approach of asymptotic integration of the elastic differential equations to derive a general theory of shells. This approach was generalized by Kaplunov et al. [12] to deal with cases other than long-wave and low-frequency approximation, such as long-wave and high-frequency approximation and short-wave approximations. Some other extension/generalization of Gol’denveizer’s approach for anisotropic shells can be found in the English version of the book by Aghalovyan [13], in which some special attention is paid on the boundary layer solutions. The proposed method presented in this book can handle the face boundary conditions with prescribed values of displacements or mixed boundary conditions of elasticity, which leads to progress in the analysis of shells resting on elastic foundations and enables novel formulations for the modulus of a foundation for layered and inhomogeneous foundations. In addition, the approach is rather efficient for tackling dynamic boundary value problems for thin walled structures, including these interacting with various physical-mechanical fields.

Vorovich and his group developed another major research direction associated with the asymptotic method. The asymptotic behavior of the solution to the problem for a hollow, finite-length cylinder subjected to an axisymmetric load and that for a thin-walled spherical shell subjected to a symmetric load were studied in [14, 15], respectively. They drew the conclusion concerning the second-order accuracy of technical shell theories such as the Vlasov theory and the Darevskii theory, and then explored a practical approach to the construction of more accurate theories. For a compressible Saint-Venant strain energy function (which contains geometric nonlinearity), Hamdouni and Millet [16, 17] deduced the models from the level of applied forces and the geometric properties of the middle surface, and obtained a constructive method of classification of shell models whose domain of validity is specified. One feature of their work is the scale of the displacements is not a data of the problem but is deduced from the scale of applied forces. For disc problems for linearized elastic materials, Kienzler and Schneider [18, 19] adopted another approach based on the series expansion of the deformed position but with no explicit kinematic assumptions. They derived the governing equations of all independent coefficients in the expansion from the 2D variational principle by first integrating out the Z variable and conducting a truncation.

An alternative asymptotic method based on the 3D weak formulation with the thickness as the small parameter has been proposed by Ciarlet and Destuynder [20] for plates and by Destuynder [21] for shells. This approach stands out by its amenability to a rigorous asymptotic analysis which provides convergence results for linear shell theories as the thickness approaches zero. The $Γ$ -convergence method was applied by Le Dret and Raoult [22, 23] and Friesecke et al. [24] to derive rigorously the shell theory for specimens with zero thickness and under appropriate hypotheses on the applied loads, where the nonlinear energy is obtained by computing the $Γ$ -limit of the sequence of 3D energies. Some high-order shell theories with thickness stretching have been proposed by Simo et al. [25] and Büchter et al. [26], who used complete 3D material laws within shell analysis and were able to approximately represent the 3D effects. Reddy and Arciniega [27] introduced a third-order shell theory based on a cubic expansion of the thickness coordinate around the middle surface, which takes into account the shear deformation and can be applied to laminated structures. For nearly traction-free cases, Steigmann [28, 29] extended the Koiter’s shell theory to the case of arbitrary material symmetry and provided an optimal third-order approximation for the 3D potential energy.

Nowadays, the unusual mechanical properties of biological soft tissues have motivated more and more research interests. In general, soft biological materials are assumed to be incompressible materials that can undergo large deformations. Furthermore, samples with thin plate or shell forms are commonly observed in biological organs or tissues, such as skins, blood vessels, leaves, and petals. To study the mechanical behaviors of these thin biological samples, one needs to adopt suitable plate or shell theories for incompressible materials. In the literature, the existing works on derived shell theory for incompressible hyperelastic materials are not so many. Makowski and Stumpf [30] formulated a finite-strain shell theory assuming that the material fibers initially normal to the shell reference surface remain straight in the deformation process. Itskov [31] proposed an orthotropic hyperelastic constitutive model, which is coupled with incompressible shell kinematics to deal with large strains and finite rotations of shell structures. Chapelle et al. [32] analyzed whether the stress assumption or the asymptotic limits of thinness can be commuted with the incompressibility condition, which justified the usages of classical shell models and a modified 3D shell models under incompressibility condition. By using the $Γ$ -convergence method, Li and Chermisi [33] rigorously derived the von Kármán shell theory for incompressible materials. Despite the existence of these representative works, there still lacks a derived shell theory for incompressible hyperelastic materials that incorporates both bending and stretching under general loadings with no ad hoc assumptions and inherits the constitutive relations from the original 3D ones naturally.

In addition to derived shell theories, another type is direct shell theories, which model the shell as a material surface with several directors modeling the deformation of a material fiber through the thickness of the shell. Based mainly on the concept of a Cosserat surface, Naghdi [34] constructed a general theory of shells by a direct approach and made some comparison with that from the 3D theory. Zhilin [35] also applied the direct approach in the formulation of the shell theory based on Cosserat’s ideas. Within the framework of the direct approach the theory of Cosserat shells can also handle incompressibility constraint, and we refer to Rubin [36] and the review by Altenbach et al. [37] for details. More recently, in the lecture notes by Altenbach and Eremeyev [38] the classical theories are summarized and various advanced theories such as direct approach to viscoelastic plates are introduced. Another is a mixed-type (between derived and direct) shell theory developed by Simmonds and co-authors (see Libai and Simmonds [39]), and such an approach can also handle the incompressibility constraint without much difficulty. One difference between a derived shell theory and a direct (or mixed-type) shell theory is that the former can directly use 3D material laws and the latter needs to use a 2D geometrically dependent constitutive law (which may not be easy to determine).

In the recent works of Dai and Song [40, 41], some consistent plate and shell theories have been proposed within the framework of 3D nonlinear elasticity. Here, the “consistency” means that the derived plate or shell equations, together with the associated boundary conditions, ensure each term in the 3D variational formulation attains a required asymptotic order, say, $O (h^{4})$ . These plate and shell theories can be applied in general loading conditions and can incorporate bending and stretching effects simultaneously. In Wang et al. [42], the consistent plate theory proposed in [40] was extended for incompressible hyperelastic materials, where the constraint equation of incompressibility and some additional variables were considered in the derivation procedure. In the current work, we shall further extend the shell theory proposed in [41] for compressible hyperelastic materials to one for incompressible hyperelastic materials. The incompressibility constraint induces an extra unknown and causes more complexity for the derivation procedure. Nevertheless, since most soft materials are incompressible, such a theory is needed. First, the 3D governing system for a shell structure made of an incompressible material will be derived through the variational approach, which is composed of the mechanical field equation and the constraint equation. Then, based on bottom surface of the shell, series expansions of the independent variables will be conducted along the thickness direction. The recursive relations of the coefficient functions in the series expansions can be derived from the 3D governing equations. Further from the top surface boundary condition, a 2D vector shell equation can be derived. The associated position and traction boundary conditions will also be proposed. It can be verified that the derived shell equation, together with the boundary conditions, ensures the required asymptotic order of the terms in the 3D variational formulation. The weak formulation of the shell equation is also established for future numerical calculations. The validity of the current shell theory is shown by studying the axisymmetric deformations of a spherical and a circular cylindrical shell made of incompressible neo-Hookean materials.

This manuscript is arranged as follows. In Section 2, the kinematics of shell samples are recalled and the 3D governing system for incompressible materials is derived. In Section 3, the 2D vector shell equation and the associated edge boundary conditions are obtained, whose consistency with the 3D variational formulation is also verified. The weak formulation of the shell equation is derived in Section 4, which can be simplified for different types of boundary conditions. In Section 5, the validity of the shell theory is shown by studying the two common shell structures made of incompressible neo-Hookean materials. Finally, some conclusions are drawn.

2. Kinematics and the 3D governing system

In the current work, we consider a homogeneous thin shell of constant thickness and composed of an incompressible hyperelastic material. In the reference configuration $K_{0}$ , the shell occupies the region $Ω \times [0, 2 h]$ , where the thickness of the shell $2 h$ is small compared with the dimensions of the top (or bottom) surface $Ω$ as well as the local radius of curvature. After the deformation, the current configuration of the shell becomes $K_{t}$ . The position of a material point in the shell is denoted by $X$ in $K_{0}$ and $x$ in $K_{t}$ . In the following notation, the Greek letters $(α, β, γ \dots)$ run from 1 to 2, whereas the Latin letters $(i, j, k \dots)$ run from 1 to 3. The repeated summation convention will be employed and a comma preceding indices means differentiation.

Following the conventional approach [5, 28], we use curvilinear coordinates $θ^{α}$ to parameterize the base surface $Ω$ of the shell. A generic point on $Ω$ is denoted by $r = r (θ^{α})$ . The tangent vectors along the coordinate lines can be calculated through $g_{α} = r_{, α} = \partial r ∕ \partial θ^{α}$ , which form the covariant basis of the tangent plane to $Ω$ . Another two vectors $g^{α}$ in the tangent plane, which satisfy the relations $g^{α} \cdot g_{β} = δ_{β}^{α}$ , form the contravariant basis of the tangent plane. Based on ${g_{α}}$ , the unit normal vector $n$ to $Ω$ can be defined according to $n = g_{1} \land g_{2} ∕ | g_{1} \land g_{2} |$ , where “ $\land$ ” means the cross product. Formally, we denote $g_{3} = g^{3} = n$ , then ${g_{i}}$ and ${g^{i}}$ ( $i = 1, 2, 3$ ) form two sets of right-handed bases.

With the above preparations, the position of a material point in $K_{0}$ can be decomposed into

X = r (θ^{α}) + Z n (θ^{α}), 0 \leq Z \leq 2 h,

(1)

where Z is the coordinate along the normal vector $n$ . At the different material points, the normal vector $n$ has different orientations. Simple calculation yields that

d n = n_{, α} d θ^{α} = (n_{, α} \otimes g^{α}) (g_{β} d θ^{β}) : = - κ d r,

(2)

where $κ$ is the curvature tensor with the following component form

κ = - n_{, α} \otimes g^{α} = κ_{α}^{β} g_{β} \otimes g^{α} = κ_{β α} g^{β} \otimes g^{α} .

(3)

Based on $κ$ , the mean and Gaussian curvatures can be defined by

H = \frac{1}{2} tr (κ) = \frac{1}{2} κ_{α}^{α}, K = \det (κ) = [{(2 H)}^{2} - κ_{α}^{β} κ_{β}^{α}] ∕ 2 .

(4)

Accordingly, the differential of the reference position vector yields that

\begin{matrix} d X = d r + Zd n + n dZ = μ d r + n dZ ≜ {\hat{g}}_{α} d θ^{α} + n dZ, \\ μ = I - Z κ, \end{matrix}

(5)

where $I$ is the rank-two identity tensor and the following covariant and contravariant base vectors have been adopted

{\hat{g}}_{α} = μ g_{α}, {\hat{g}}^{α} = μ^{- T} g^{α} .

(6)

Here ${\hat{g}}_{α}$ and ${\hat{g}}^{α}$ are also orthogonal to $n$ . Note that the assumption on the thickness of the shell implies that $| 2 h κ_{α}^{β} | < 1$ , hence the inverse $μ^{- 1}$ is well defined. The area element on $Ω$ and the volume element can then be calculated through

\begin{matrix} dA = | g_{1} \land g_{2} | d θ^{1} d θ^{2} = \sqrt{g} d θ^{1} d θ^{2} \\ dV = μ (Z) dA dZ, μ (Z) = \det (μ) = 1 - 2 HZ + K Z^{2}, \end{matrix}

(7)

where $g$ is the determinant of the metric tensor $(g_{α β})$ associated with the base vectors $g_{α}$ .

The position vector of a material point in $K_{t}$ can be written as $x = x (θ^{α}, Z)$ . The deformation gradient tensor $F$ can then be obtained from (5) and the relation $d x = Fd X$ , which is given by

F = x_{, α} \otimes {\hat{g}}^{α} + \frac{\partial x}{\partial Z} \otimes n = (\nabla x) μ^{- 1} + \frac{\partial x}{\partial Z} \otimes n,

(8)

where $\nabla = g^{α} \partial ∕ \partial θ^{α}$ is the in-plane 2D gradient on the base surface $(\nabla x = x_{, α} \otimes g^{α})$ .

For an incompressible material, we have the following constraint equation

R (F) = Det (F) - 1 = 0, in Ω \times [0, 2 h] .

(9)

Further suppose that the material has the strain-energy density function $Φ (F)$ , then the associated elastic moduli can be defined by

A^{n} (F^{(0)}) = \frac{\partial^{n + 1} Φ}{\partial F^{n + 1}} |_{F = F^{(0)}}, n = 1, 2, \dots .

(10)

It is assumed that the strain-energy function satisfies the strong-ellipticity condition

a \otimes b : A^{1} (F) [a \otimes b] > 0, for all a \otimes b \neq 0,

(11)

where the colon between second-order tensors represents the trace of their tensor product $A : B = tr (A B)$ , and the square bracket after a modulus tensor represents the operations

{A^{1} [A]}_{ij} = A_{ijkl}^{1} A_{lk}, {A^{2} [A, B]}_{ij} = A_{ijklmn}^{2} A_{lk} B_{nm} .

For the case of dead-loading, suppose that $q^{\pm}$ are the external loads applied on the top and the bottom surfaces of the shell. The edge boundary $\partial Ω$ of the base surface is composed of two parts: $\partial Ω_{0}$ with the prescribed position and $\partial Ω_{q}$ subject to the applied traction $q$ . The 3D potential energy E of the shell is then given by

\begin{matrix} E = \bar{Φ} + \bar{V}, \bar{Φ} = \int_{Ω} \int_{0}^{2 h} Φ (F) μ (Z) dZ dA, \\ \bar{V} = - \int_{Ω} {q^{-} (θ^{α}) \cdot x (θ^{α}, 0) + μ (2 h) q^{+} (θ^{α}) \cdot x (θ^{α}, 2 h)} dA \\ - \int_{\partial Ω_{q}} \int_{0}^{2 h} q (s, Z) \cdot x (s, Z) da . \end{matrix}

(12)

On $\partial Ω$ , let s be the arclength variable, and $τ$ and $ν$ be the unit tangent and outward normal vectors (which satisfy $ν = τ \land n$ ). Then, we have $d r = τ ds$ along $\partial Ω$ . Further from (5), it can be obtained that $d X$ consists of $d X_{τ} = μ τ ds$ along the tangent direction and $d X_{n} = n dZ$ along the outward normal direction. Through some conventional derivations, we have

\begin{matrix} da = \sqrt{1 - 2 κ_{τ} Z + c_{τ} Z^{2}} dZds = : \sqrt{g_{τ}} dZds, \\ κ_{τ} = τ \cdot (κ τ), c_{τ} = (κ τ) \cdot (κ τ) . \end{matrix}

(13)

To derive the equilibrium state of the shell, we need to calculate the minimum of the potential energy functional E under the constraint condition (9). Thus, the following generalized potential energy functional is considered

Ψ (x, p; X) = \bar{Φ} + \bar{V} - \int_{Ω} \int_{0}^{2 h} p (X) R (F) μ (Z) dZdA,

(14)

where $p (X)$ plays the role of a Lagrangian multiplier. The governing system of the shell model can then be derived by calculating the variations of $Ψ$ with respect to the independent variables $x$ and p.

First, the variation of $Ψ$ with respect to $x$ yields that

\begin{matrix} \frac{δ Ψ}{δ x} = - \int_{Ω} \int_{0}^{2 h} Div S \cdot δ x μ (Z) dZdA - \int_{Ω} (S^{T} n |_{Z = 0} + q^{-}) \cdot δ x (θ^{α}, 0) dA \\ + \int_{Ω} (S^{T} n |_{Z = 2 h} - q^{+}) \cdot δ x (θ^{α}, 2 h) μ (2 h) dA + \int_{\partial Ω_{0}} \int_{0}^{2 h} S^{T} N \cdot δ x (s, Z) da \\ + \int_{\partial Ω_{q}} \int_{0}^{2 h} (S^{T} N - q) \cdot δ x (s, Z) da, \end{matrix}

(15)

where

S = \frac{\partial Φ}{\partial F} - p \frac{\partial R}{\partial F}

(16)

is the nominal stress tensor of the incompressible material [43]. Owing to the arbitrariness of $δ x$ in (15), the following 3D equilibrium equation together with the boundary conditions can be obtained

\begin{matrix} Div S = 0, & in Ω \times [0, 2 h], \\ S^{T} n |_{Z = 0} = - q^{-}, & on Ω, \\ S^{T} n |_{Z = 2 h} = q^{+}, & on Ω, \\ x = b (s, Z), & on \partial Ω_{0} \times [0, 2 h], \\ S^{T} N = q (s, Z), & on \partial Ω_{q} \times [0, 2 h], \end{matrix}

(17)

where $b (s, Z)$ is the prescribed position on $\partial Ω_{0}$ .

Moreover, the variation of $Ψ$ with respect to p yields that

\frac{δ Ψ}{δ p} = - \int_{Ω} \int_{0}^{2 h} R (F) δ p μ (Z) dZdA,

(18)

from which one can obtain the constraint equation (9). Equations (9) and (17) formulate the 3D governing system of the current shell model, which contains two independent variables $x (X)$ and $p (X)$ .

3. The 2D vector shell equations

In this section, a consistent shell theory will be derived from the 3D governing system (9) and (17), which includes the 2D shell equation and some appropriate boundary conditions. Here, the criterion of consistency requires that [40], each term in the first variations of the energy functional (i.e., (15) and (18)) should be of a required asymptotic order (say, $O (h^{4})$ ) separately for the approximation. The derivation procedure is similar to that for compressible materials [41]. However, owing to incompressibility, one now has to handle an extra constraint and an additional unknown. The major steps of the derivation procedure are given in the following subsections.

3.1. Derivation of the 2D vector shell equation

To eliminate the variable $‘ Z ’$ from the 3D governing system, we consider the following series expansions of the independent variables $x (X)$ and $p (X)$ with respect to Z (sufficient smoothness is assumed):

\begin{matrix} x (X) = x^{(0)} (θ^{α}) + Z x^{(1)} (θ^{α}) + \frac{1}{2} Z^{2} x^{(2)} (θ^{α}) + \frac{1}{6} Z^{3} x^{(3)} (θ^{α}) + \frac{1}{24} Z^{4} x^{(4)} (θ^{α}) + O (Z^{5}), \\ p (X) = p^{(0)} (θ^{α}) + Z p^{(1)} (θ^{α}) + \frac{1}{2} Z^{2} p^{(2)} (θ^{α}) + \frac{1}{6} Z^{3} p^{(3)} (θ^{α}) + \frac{1}{24} Z^{4} p^{(4)} (θ^{α}) + O (Z^{5}), \end{matrix}

(19)

where ${(\cdot)}^{(n)} = \partial^{n} (\cdot) ∕ \partial Z^{n} |_{Z = 0} (n = 0, 1, \dots, 4)$ . Accordingly, the deformation gradient tensor can be expanded as

F (X) = F^{(0)} (θ^{α}) + Z F^{(1)} (θ^{α}) + \frac{1}{2} Z^{2} F^{(2)} (θ^{α}) + \frac{1}{6} Z^{3} F^{(3)} (θ^{α}) + O (Z^{4}) .

(20)

By substituting (19)₁ and (20) into (8) and comparing the coefficients of Z, we obtain the relations

F^{(m)} = \sum_{i = 0}^{m} \frac{m!}{i!} (\nabla x^{(i)}) κ^{m - i} + x^{(m + 1)} \otimes n, m = 0, 1, 2, 3 .

(21)

It can be seen from (21) that the dependence of $F^{(m)}$ on $x^{(m + 1)}$ is linearly algebraic. We also consider the series expansion of the nominal stress tensor $S$ in the following form

S (X) = S^{(0)} (θ^{α}) + Z S^{(1)} (θ^{α}) + \frac{1}{2} Z^{2} S^{(2)} (θ^{α}) + \frac{1}{6} Z^{3} S^{(3)} (θ^{α}) + O (Z^{4}) .

(22)

By further substituting (19) and (22) into (16) and comparing the coefficients of Z, we obtain the following expressions for $S^{(n)}$

\begin{matrix} S^{(0)} (θ^{α}) = A^{(0)} - p^{(0)} ℝ^{(0)}, \\ S^{(1)} (θ^{α}) = \bar{A} [F^{(1)}] - p^{(1)} ℝ^{(0)}, \\ S^{(2)} (θ^{α}) = \bar{A} [F^{(2)}] - p^{(2)} ℝ^{(0)} + A^{(2)} [F^{(1)}, F^{(1)}] - p^{(0)} R^{(2)} [F^{(1)}, F^{(1)}] - 2 p^{(1)} R^{(1)} [F^{(1)}], \\ S^{(3)} (θ^{α}) = \bar{A} [F^{(3)}] - p^{(3)} ℝ^{(0)} + 3 A^{(2)} [F^{(1)}, F^{(2)}] + A^{(3)} [F^{(1)}, F^{(1)}, F^{(1)}] \\ - 3 p^{(0)} R^{(2)} [F^{(1)}, F^{(2)}] - p^{(0)} R^{(3)} [F^{(1)}, F^{(1)}, F^{(1)}] - 3 p^{(1)} R^{(1)} [F^{(2)}] \\ - 3 p^{(1)} R^{(2)} [F^{(1)}, F^{(1)}] - 3 p^{(2)} R^{(1)} [F^{(1)}], \end{matrix})

(23)

where $R^{(n)} (n = 1, 2, 3)$ are defined by replacing $Φ (F)$ with the constraint $R (F)$ in (10), $\bar{A} = A^{(1)} - p^{(0)} R^{(1)}$ is a new modulus which resembles that of $A^{(1)}$ in the compressible case, and

A^{(0)} (F^{(0)}) = \frac{\partial Φ}{\partial F} |_{F = F^{(0)}}, ℝ^{(0)} (F^{(0)}) = \frac{\partial R}{\partial F} |_{F = F^{(0)}} = Det (F^{(0)}) {(F^{(0)})}^{- 1} .

From (23) and (21), we can observe that $S^{(n)}$ depends linearly on $x^{(n + 1)} (n = 1, 2, 3)$ as well. Owing to the series expansions (19), we obtain a total of 19 unknowns (including five vectors $x^{(m)} (m = 0, \dots, 4)$ and four scalars $p^{(m)} (m = 0, \dots, 3)$ ), thus we also need to find 19 equations to form a closed partial differential equation (PDE) system.

First, by substituting (22) into the bottom traction condition (17)₂, we obtain

\begin{matrix} {S^{(0)}}^{T} n & = {A^{(0)} (F^{(0)}) - p^{(0)} ℝ^{(0)} (F^{(0)})}^{T} n = - q^{-}, \end{matrix}

(24)

which provides three algebraic equations for the unknowns $x^{(1)}$ ( $F^{(0)}$ depends on $x^{(1)}$ linearly) and $p^{(0)}$ . To ease derivations in the following, we use the notation

\begin{matrix} y (x^{(0)}) : = ℝ^{(0) T} n = Det (F^{(0)}) {(F^{(0)})}^{- T} n = F^{{(0)}^{*}} n = F^{{(0)}^{*}} (g_{1} \land g_{2}) ∕ \sqrt{g} \\ = (F^{(0)} g_{1}) \land (F^{(0)} g_{2}) ∕ \sqrt{g} = x_{, 1}^{(0)} \land x_{, 2}^{(0)} ∕ \sqrt{g}, \end{matrix}

where “ $*$ ” represents the adjugate tensor [44]. Immediately, (24) reduces to

A^{(0) T} (F^{(0)}) n - p^{(0)} y (x^{(0)}) = - q^{-} .

(25)

Next, by substituting (22) into the mechanical field equation (17)₁ and equating the coefficients of $Z^{m} (m = 0, 1, 2, \dots)$ to be zero, one obtains

\sum_{i = 0}^{m} \frac{m!}{i!} [{(κ^{m - i})}^{T} g^{α}] \cdot S_{, α}^{(i)} + S^{(m + 1) T} n = 0 m = 0, 1, 2 .

(26)

In particular, the first three equations can be expressed explicitly as

\begin{matrix} \nabla \cdot S^{(0)} + S^{(1) T} n = 0, \\ (κ^{T} g^{α}) \cdot S_{, α}^{(0)} + \nabla \cdot S^{(1)} + S^{(2) T} n = 0, \\ 2 ({(κ^{2})}^{T} g^{α}) \cdot S_{, α}^{(0)} + 2 (κ^{T} g^{α}) \cdot S_{, α}^{(1)} + \nabla \cdot S^{(2)} + S^{(3) T} n = 0, \end{matrix}

(27)

where $\nabla \cdot S^{(0)} = g^{α} \cdot S_{, α}^{(0)}$ denotes the 2D divergence. In general, the three relations in (27) provide nine linear algebraic equations for the unknowns $x^{(i + 1)}$ and $p^{(i)} (i = 1, 2, 3)$ .

Third, by substituting (20) into the constraint equation (9) and equating the coefficients of $Z^{m} (m = 0, 1, 2, 3)$ to be zero, one obtains

\begin{matrix} R (F^{(0)}) = (x_{, 1}^{(0)} \land x_{, 2}^{(0)}) \cdot x^{(1)} - 1 = \sqrt{g} y \cdot x^{(1)} - 1 = 0, \\ ℝ^{(0)} : F^{(1)} = y \cdot x^{(2)} + ℝ^{(0)} : ((\nabla x^{(0)}) κ + \nabla x^{(1)}) = 0, \\ ℝ^{(0)} : F^{(2)} + F^{(1)} : R^{(1)} [F^{(1)}] = 0, \\ ℝ^{(0)} : F^{(3)} + 3 F^{(2)} : R^{(1)} [F^{(1)}] + F^{(1)} : R^{(2)} [F^{(1)}, F^{(1)}] = 0, \end{matrix}

(28)

which provide four additional linear algebraic equations for the unknowns $x^{(m + 1)} (m = 0, 1, 2, 3)$ .

From (27)_1,2 and (28)_2,3, we can obtain the expressions of $x^{(m + 1)}$ and $p^{(m)} (m = 1, 2)$ in terms of $x^{(0)}$ , $x^{(1)}$ and $p^{(0)}$ , which are given by

\begin{matrix} x^{(2)} = - B^{- 1} f^{(2)} + p^{(1)} B^{- 1} y, \\ p^{(1)} = \frac{1}{y} (y \cdot B^{- 1} f^{(2)} - ℝ^{(0)} : [(\nabla x^{(0)}) κ + \nabla x^{(1)}]), \\ x^{(3)} = - B^{- 1} f^{(3)} + p^{(2)} B^{- 1} y, \\ p^{(2)} = \frac{1}{y} (y \cdot B^{- 1} f^{(3)} - ℝ^{(0)} : [\sum_{i = 0}^{2} \frac{2}{i!} (\nabla x^{(i)}) κ^{2 - i}] - F^{(1)} : R^{(1)} [F^{(1)}]), \end{matrix}

(29)

where $y = y \cdot B^{- 1} y$ and

\begin{matrix} B x = {\bar{A} [x \otimes n]}^{T} n, \Rightarrow (B) |_{ij} = {\bar{A}}_{3 i 3 j}, \\ f^{(2)} = {\bar{A} [(\nabla x^{(0)}) κ + \nabla x^{(1)}]}^{T} n + \nabla \cdot S^{(0)}, \\ f^{(3)} = {\bar{A} [\sum_{i = 0}^{2} \frac{2}{i!} (\nabla x^{(i)}) κ^{2 - i}]}^{T} n + (κ^{T} g^{α}) \cdot S_{, α}^{(0)} + \nabla \cdot S^{(1)} \\ + {(A^{(2)} - p^{(0)} R^{(2)}) [F^{(1)}, F^{(1)}] - 2 p^{(1)} R^{(1)} [F^{(1)}]}^{T} n . \end{matrix}

(30)

Note that the strong-ellipticity condition (11) implies that $B$ is invertible and positive-definite. The explicit expressions of $x^{(4)}$ and $p^{(3)}$ are omitted here since they are not needed in the following derivations. To derive the explicit expressions of $x^{(1)}$ and $p^{(0)}$ , one needs to consider Equations (25) and (28)₁. However, these two equations are nonlinear algebraic equations, which can only be solved when the concrete form of the strain-energy function is given. Nevertheless, the strong-ellipticity condition and the implicit function theorem guarantee that they can be uniquely solved (see [42]).

Finally, substitution of (22) into the top traction condition (17)₃ yields that

S^{(0) T} n + 2 h S^{(1) T} n + 2 h^{2} S^{(2) T} n + \frac{4}{3} h^{3} S^{(3) T} n + O (h^{4}) = q^{+} .

(31)

By multiplying the factor $μ (2 h) = 1 - 4 h H + 4 h^{2} K$ on the two sides of (31) and through some delicate manipulations [41], we can derive the following 2D vector shell equation

\nabla \cdot (1 (\tilde{S} + h κ \hat{S})) = - \tilde{q},

(32)

where

\begin{matrix} \tilde{q} & = \frac{1}{2 h} [(1 - 4 h H + 4 h^{2} K) q^{+} + q^{-}], \\ \tilde{S} & = (1 - 2 h H) S^{(0)} + h (1 - \frac{8}{3} h H) S^{(1)} + \frac{2}{3} h^{2} S^{(2)}, \\ \hat{S} & = S^{(0)} + \frac{4}{3} h S^{(1)} . \end{matrix}

(33)

The quantity $\tilde{S} + h κ \hat{S}$ is considered as the averaged stress and $\tilde{q}$ is the weighted average of surface traction (the weight depends on the surface area). Formally, this 2D shell equation can also be deduced by an integration of the mechanical equilibrium equation (17)₁ multiplied by $μ (2 h)$ over Z, thus it preserves the 3D force equilibrium in the sense of through-thickness average. By substituting the recursion relations (29) and the expressions of $x^{(1)}$ and $p^{(0)}$ into (32), one obtains a fourth-order differential equation for $x^{(0)}$ . Once the shell equation is solved, an up to $O (h^{2})$ -correct result for $x^{(0)}$ can be obtained. Correspondingly, the other components $x^{(m)} (m = 1, 2, 3)$ can also be derived. Therefore, the distribution of 3D current position vector can be recovered. In addition, as mentioned in [41], the shell equations are actually local force equilibrium equations for a shell element and can be easily adapted to some non-dead-loading cases, when the base surface is subjected to a pressure or sitting on a foundation.

3.2. Boundary conditions

In addition to the 2D vector shell equation, appropriate boundary conditions should also be deduced from the original 3D lateral surface conditions. On both the position boundary $\partial Ω_{0}$ and the traction boundary $\partial Ω_{q}$ , two conditions regarding $x^{(0)}$ or its derivatives are required since the shell equation is of fourth order. In fact, the scheme for proposing the boundary conditions has been introduced in the authors’ previous works [40, 41]. Here, to be self-contained, we only give a simple introduction.

Case 1. Prescribed position in the 3D formulation. Suppose that on $\partial Ω_{0} \times [0, 2 h]$ the position $b$ is prescribed, we define

\bar{x} (s) = \frac{1}{2 h} \int_{0}^{2 h} x (s, Z) dZ = x^{(0)} + h x^{(1)} + \frac{2}{3} h^{2} x^{(2)} + O (h^{3}) .

(34)

Then, the following two boundary conditions can be proposed

\begin{matrix} x^{(0)} = b^{(0)} (s), \\ \bar{x} = \bar{b} on \partial Ω_{0} \Leftrightarrow x^{(1)} + \frac{2}{3} h x^{(2)} + \frac{1}{3} h^{2} x^{(3)} + O (h^{3}) = \frac{1}{h} (\bar{b} - b^{(0)}), \end{matrix}

(35)

where $b^{(0)} = b |_{Z = 0}$ and $\bar{b}$ is defined in the same way as $\bar{x}$ . The second condition contains up to the third-order derivatives of $x^{(0)}$ upon using the recursion relations.

Case 2. Prescribed traction in the 3D formulation. Suppose that on $\partial Ω_{q} \times [0, 2 h]$ , the traction $q$ is specified and is sufficiently smooth. Denote the first four coefficients of its series expansion in terms of Z by $q^{(i)} (i = 0, 1, 2, 3)$ . For future convenience, we denote

\begin{matrix} p = \sqrt{g_{τ}} S^{T} N = {[S^{(0)} + Z S^{(1)} + \frac{1}{2} Z^{2} S^{(2)} + \dots]}^{T} [ν + Z ν_{1}] \\ = S^{(0) T} ν + Z [S^{(1) T} ν + S^{(0) T} ν_{1}] + \frac{1}{2} Z^{2} [S^{(2) T} ν + 2 S^{(1) T} ν_{1}] + O (Z^{3}) \\ ≜ p^{(0)} + Z p^{(1)} + \frac{1}{2} Z^{2} p^{(2)} + O (Z^{3}), \end{matrix}

(36)

and the averaged traction

\begin{matrix} q_{0} = \frac{1}{2 h} \int_{0}^{2 h} q (s, Z) \sqrt{g_{τ}} dZ \\ = \frac{1}{2 h} \int_{0}^{2 h} [q^{(0)} + Z q^{(1)} + \frac{1}{2} Z^{2} q^{(2)} + \dots] \sqrt{1 - 2 κ_{τ} Z + c_{τ} Z^{2}} dZ \\ = q^{(0)} + h [q^{(1)} - κ_{τ} q^{(0)}] + \frac{2}{3} h^{2} [q^{(2)} - 2 κ_{τ} q^{(1)} + κ_{τ ν}^{2} q^{(0)}] + O (h^{3}) \cdot \end{matrix}

(37)

Then, the following two boundary conditions can be proposed

\begin{matrix} S^{(0) T} ν = q^{(0)}, on \partial Ω_{q}, \\ \bar{p} ≜ p^{(0)} + h p^{(1)} + \frac{2}{3} h^{2} p^{(2)} + O (h^{3}) = q_{0} \cdot \end{matrix}

(38)

Alternatively, we can define the resultant bending moment about the middle line of the lateral surface as

\begin{matrix} m (s) = \frac{1}{2 h} \int_{0}^{2 h} q (s, Z) \sqrt{g_{τ}} (Z - h) dZ \\ = \frac{1}{3} h^{2} [q^{(1)} - κ_{τ} q^{(0)}] + \frac{1}{3} h^{3} [q^{(2)} - 2 κ_{τ} q^{(1)} + κ_{τ ν}^{2} q^{(0)}] + O (h^{4}) \cdot \end{matrix}

(39)

Then, the first condition in (38) can be replaced by

\frac{1}{2 h} \int_{0}^{2 h} (Z - h) p dZ = m (s) \Leftrightarrow \frac{1}{3} p^{(1)} + \frac{1}{3} h p^{(2)} + \frac{1}{5} h^{2} p^{(3)} + O (h^{3}) = \frac{m (s)}{h^{2}} .

(40)

3.3. Examination of the consistency

With the 2D vector shell equation and the associated boundary conditions, according to the consistency criterion introduced before, we need to analyze the asymptotic orders of the terms in the variations (15) and (18).

For the first term on the right-hand side (r.h.s.) of (15), we consider the series expansions of $Div S$ in terms of Z. As the first three terms in the series expansion (26) have been used together with (28)_(2-4) to obtain the recursion relations of $x^{()} (i = 2, 3, 4)$ and $p^{(i)} (i = 1, 2, 3)$ , we have $Div S = O (Z^{3})$ . Thus, the first term in (15) is of $O (h^{4})$ . The second term on the r.h.s. of (15) is exactly equal to zero because (24) together with (28)₁ has been used to derive the expressions of $x^{(1)}$ and $p^{(0)}$ . In (31), $S^{T} n$ on the top surface has been expanded to $O (h^{4})$ , which implies that the third term in (15) is $O (h^{4})$ .

Next, we examine the asymptotic order of the fourth term in (15). Based on the definition of $p$ , we obtain

\begin{matrix} \int_{0}^{2 h} S^{T} N \cdot δ x \sqrt{g_{τ}} dZ = \int_{0}^{2 h} p \cdot δ x dZ \\ = \int_{0}^{2 h} p^{(0)} \cdot δ x dZ + \int_{0}^{2 h} Z p^{(1)} \cdot [δ x^{(0)} + Z δ x^{(1)}] dZ + \int_{0}^{2 h} \frac{1}{2} Z^{2} p^{(2)} \cdot δ x^{(0)} dZ + O (h^{4}) . \end{matrix}

(41)

From the position boundary conditions (35), it can be obtained that $δ x^{(0)} = 0$ , $δ x^{(1)} = O (h)$ and $δ x^{(1)} + 2 ∕ 3 h δ x^{(2)} = O (h^{2})$ . By substituting these results into (41), it is easy to check that all three terms on the r.h.s. of (41) are of $O (h^{4})$ . Thus, the fourth term in (15) also satisfies the consistency condition.

To examine the asymptotic order of the fifth term in (15), we denote $\hat{q} = \sqrt{g_{τ}} (S^{T} N - q) = p - \sqrt{g_{τ}} q$ . The coefficients of the series expansion of $\hat{q}$ are represented by ${\hat{q}}^{(i)} (i = 0, \dots, 3)$ . Then, the integration of the fifth term in (15) can be rewritten as

\begin{matrix} \int_{0}^{2 h} \hat{q} \cdot δ x dZ = \int_{0}^{2 h} \hat{q} \cdot δ x^{(0)} dZ + \int_{0}^{2 h} Z [{\hat{q}}^{(0)} + Z {\hat{q}}^{(1)}] \cdot δ x^{(1)} dZ \\ + \int_{0}^{2 h} \frac{1}{2} Z^{2} {\hat{q}}^{(0)} \cdot δ x^{(2)} dZ + O (h^{4}) . \end{matrix}

(42)

From the traction boundary conditions (38), it can be obtained that ${\hat{q}}^{(0)} = 0$ , ${\hat{q}}^{(1)} = O (h)$ and ${\hat{q}}^{(1)} + 2 ∕ 3 h {\hat{q}}^{(2)} = O (h^{2})$ . By substituting these results into (42), it is easy to check that all three terms on the r.h.s. of (42) are of $O (h^{4})$ . Thus, the fifth term in (15) also satisfies the consistency condition.

For the variation (18), we have considered the series expansion of $R (F)$ in terms of Z, where the coefficients of $Z^{m} (m = 0, 1, 2, 3)$ are set to be zero to derive the recursion relations of $x^{(m + 1)} (m = 0, 1, 2, 3)$ . Thus, the variation (18) is also of $O (h^{4})$ , which satisfies the consistency condition.

In summary, the 2D vector shell equation (32) together with the edge boundary conditions (35) or (38) ensures the consistency of the 3D weak formulation of the governing system to an asymptotic order of $O (h^{4})$ . Note that no higher-order stress resultants are involved in above boundary conditions, in contrast to some high-order shell theories.

4. Associated weak formulation

To be prepared for future numerical calculations, we shall derive the associated weak formulation for the previous 2D vector shell system in this section. Furthermore, based on the weak formulation, suitable boundary conditions can also be proposed for some practical loading cases.

First, by multiplying both sides of the shell equation (32) with $ξ = δ x^{(0)}$ and calculating the integration over the region $Ω$ , we obtain

\int_{Ω} [\nabla \cdot (1 (\tilde{S} + h κ \hat{S}))} \cdot ξ dA = - \int_{Ω} \tilde{q} \cdot ξ dA .

(43)

By adopting the 2D divergence theorem and through some calculations, (43) can be rewritten as

\int_{Ω} (\tilde{S} + h κ \hat{S}) : \nabla ξ dA - \int_{Ω} \tilde{q} \cdot ξ dA = \int_{\partial Ω} [{\tilde{S}}^{T} ν + h {(κ \hat{S})}^{T} ν] \cdot ξ ds = \int_{\partial Ω} \bar{p} \cdot ξ ds

(44)

where the term in square brackets corresponds to $\bar{p}$ .

In general, the weak formulation associated with the fourth-order shell equation (32) should only contain up to the second-order derivatives of $x^{(0)}$ . However, the weak formulation (44) involves the third-order derivatives, which originate from the terms $F^{(2)}$ and $p^{(2)}$ in $S^{(2)}$ . Next, we intend to eliminate these third-order derivative terms.

To extract the higher-order derivative terms from $F^{(2)}$ and $p^{(2)}$ , we first decompose $f^{(3)}$ in ${(30)}_{3}$ as follows

\begin{matrix} f^{(3)} = (κ^{T} g^{α}) \cdot S_{, α}^{(0)} + \nabla \cdot S^{(1)} + {\bar{A} [\sum_{i = 0}^{2} \frac{2}{i!} (\nabla x^{(i)}) κ^{2 - i}]}^{T} n \\ + {(A^{(2)} - p^{(0)} R^{(2)}) [F^{(1)}, F^{(1)}] - 2 p^{(1)} R^{(1)} [F^{(1)}]}^{T} n \\ = \nabla \cdot (κ S^{(0)}) - (\nabla \cdot κ) \cdot S^{(0)} + \nabla \cdot (1 S^{(1)}) - 2 H S^{(1) T} n \\ + {(A^{(2)} - p^{(0)} R^{(2)}) [F^{(1)}, F^{(1)}] - 2 p^{(1)} R^{(1)} [F^{(1)}]}^{T} n \\ + {\bar{A} [\sum_{i = 0}^{2} \frac{2}{i!} (\nabla x^{(i)}) κ^{2 - i}]}^{T} n \\ = {\nabla \cdot [κ S^{(0)} + 1 S^{(1)} - 2 H 1 S^{(0)}] + {(\bar{A} [\nabla x^{(2)}])}^{T} n} \\ + {{(A^{(2)} - p^{(0)} R^{(2)}) [F^{(1)}, F^{(1)}] - 2 p^{(1)} R^{(1)} [F^{(1)}]}^{T} n \\ - (\nabla \cdot κ) \cdot S^{(0)} - 2 H S^{(1) T} n + \nabla \cdot (2 H 1 S^{(0)}) \\ + {\bar{A} [2 (\nabla x^{(1)}) κ + 2 (\nabla x^{(0)}) κ^{2}]}^{T} n} \\ ≜ f_{2}^{(3)} + f_{1}^{(3)}, \end{matrix}

where we use the identity $\nabla \cdot S^{(1)} = \nabla \cdot (1 S^{(1)}) - 2 H S^{(1) T} n$ in the second equation and the two braces define $f_{2}^{(3)}$ and $f_{1}^{(3)}$ , respectively. In the above derivation, the term $\nabla \cdot (2 H 1 S^{(0)})$ is added and subtracted to obtain the desired form. Based on ${(29)}_{3, 4}$ , we introduce

p^{(2)} = p_{1}^{(2)} + p_{2}^{(2)}, F^{(2)} = F_{1}^{(2)} + F_{2}^{(2)},

(45)

where

\begin{matrix} p_{1}^{(2)} & = \frac{1}{y} {y \cdot B^{- 1} f_{1}^{(3)} - ℝ^{(0)} : [2 (\nabla x^{(1)}) κ + 2 (\nabla x^{(0)}) κ^{2}] - F^{(1)} : R^{(1)} [F^{(1)}]}, \\ p_{2}^{(2)} & = \frac{1}{y} {y \cdot B^{- 1} f_{2}^{(3)} - ℝ^{(0)} : \nabla x^{(2)}}, \\ F_{1}^{(2)} & = 2 (\nabla x^{(1)}) κ + 2 (\nabla x^{(0)}) κ^{2} + p_{1}^{(2)} B^{- 1} y \otimes n - B^{- 1} f_{1}^{(3)} \otimes n, \\ F_{2}^{(2)} & = \nabla x^{(2)} - B^{- 1} f_{2}^{(3)} \otimes n + p_{2}^{(2)} B^{- 1} y \otimes n . \end{matrix}

(46)

It can be found that the third-order derivative terms of $x^{(0)}$ are only involved in $p_{2}^{(2)}$ and $F_{2}^{(2)}$ . To go further, we rewrite

(\tilde{S} + h κ \hat{S}) : \nabla ξ = W_{1} (\nabla \nabla x^{(0)}, \nabla ξ) + W_{2} (\nabla \nabla \nabla x^{(0)}, \nabla ξ),

(47)

where

\begin{matrix} W_{1} = ((1 - 2 h H) S^{(0)} + h (1 - \frac{8}{3} h H) S^{(1)} + h κ \hat{S}) : \nabla ξ \\ + \frac{2}{3} h^{2} {(A^{(2)} - p^{(0)} R^{(2)}) [F^{(1)}, F^{(1)}] - 2 p^{(1)} R^{(1)} [F^{(1)}]} : \nabla ξ \\ + \frac{2}{3} h^{2} {\bar{A} [F_{1}^{(2)}] - p_{1}^{(2)} ℝ^{(0)}} : \nabla ξ, \\ W_{2} = \frac{2}{3} h^{2} {\bar{A} [F_{2}^{(2)}] - p_{2}^{(2)} ℝ^{(0)}} : \nabla ξ . \end{matrix}

To eliminate the third-order derivative terms of $x^{(0)}$ , we then substitute the expressions of $F_{2}^{(2)}$ and $p_{2}^{(2)}$ into (47). Further manipulations yield that

W_{2} = \frac{2}{3} h^{2} {S_{0} : \nabla x^{(2)} + η \cdot [\nabla \cdot (1 S^{(1)} - 2 H 1 S^{(0)} + κ S^{(0)})]},

(48)

and

\begin{matrix} S_{0} & = \bar{A} [\nabla ξ + η \otimes n] - ζ ℝ^{(0)}, η = - B^{- 1} {\bar{A} [\nabla ξ]}^{T} n + ζ B^{- 1} y, \\ ζ & = \frac{1}{y} (\bar{A} [B^{- 1} y \otimes n] - ℝ^{(0)}) : \nabla ξ ≜ ℙ : \nabla ξ . \end{matrix}

(49)

In fact, it can be proved that

δ p^{(0)} = ζ, δ x^{(1)} = η, δ S^{(0)} = S_{0} .

(50)

By denoting

W_{3} = - \frac{2}{3} h^{2} {[\nabla \cdot (1 S_{0})] \cdot x^{(2)} + (S^{(1)} - 2 H S^{(0)} + κ S^{(0)}) : \nabla η}

(51)

and through the integration by parts, we obtain

\begin{matrix} \int_{Ω} (W_{2} - W_{3}) dA = \frac{2}{3} h^{2} \int_{\partial Ω} (S_{0}^{T} ν) \cdot x^{(2)} + [{(S^{(1)} - 2 H S^{(0)} + κ S^{(0)})}^{T} ν] \cdot η ds \\ = \frac{2}{3} h^{2} \int_{\partial Ω} (S_{0}^{T} ν) \cdot x^{(2)} + p^{(1)} \cdot η ds . \end{matrix}

(52)

It can be seen from (52) that the third-order derivatives of $x^{(0)}$ have been eliminated. By substituting (47) and (52) into (44), the following 2D weak formulation can be derived

\int_{Ω} (W_{1} + W_{3} - \tilde{q} \cdot ξ) dA = \int_{\partial Ω} \bar{p} \cdot ξ - \frac{2}{3} h^{2} (p^{(1)} \cdot η + (S_{0}^{T} ν) \cdot x^{(2)}) ds .

(53)

By considering the boundary conditions, the weak formulation (53) can be further simplified. Next, we shall consider the boundary conditions in the following two distinct cases and derive the reduced weak formulations.

Case 1. Edge position and traction in the 3D formulation are known

In this case, from (35), it is easy to deduce that $ξ = δ x^{(0)} = 0$ and $η = δ x^{(1)} = O (h)$ on $\partial Ω_{0}$ , which together with (49) yields that $\nabla ξ = O (h)$ and $S_{0} = O (h)$ . Consequently, the integration on $\partial Ω_{0}$ in (53) is of $O (h^{3})$ and can be neglected and we can readily replace $\partial Ω$ by $\partial Ω_{q}$ in (53).

On $\partial Ω_{q}$ , it follows from (38) and (40) that $S_{0}^{T} ν = δ [S^{(0) T} ν] = O (h^{2})$ . Thus, the third term in the boundary integral can be neglected. In addition, replacing $p^{(1)}$ by the condition in (40) only causes a higher-order correction. Then, the 2D weak formulation (53) reduces to

\int_{Ω} (W_{1} + W_{3} - \tilde{q} \cdot ξ) dA = \int_{\partial Ω_{q}} (\bar{p} \cdot ξ - \frac{2}{3} h^{2} p^{(1)} \cdot η) ds = \int_{\partial Ω_{q}} (q_{0} \cdot ξ - 2 m \cdot η) ds,

(54)

where (38)₂ and (40) has been used to derive the last equality. If the two conditions in (38) are adopted, we should replace $m$ by a combination of $q_{0}$ and $q^{(0)}$ .

Case 2. Edge position and traction in the 3D formulation are unknown

In many practical situations, the exact traction distribution (e.g., a pinned edge) or displacement distribution (e.g., a clamped edge) on the lateral surface of the shell is unknown. In this case, we should resort to the weak formulation to propose the so-called natural boundary conditions. For that purpose, we need to recast the boundary integral in (53) in terms of $ξ$ and its normal derivative $ξ_{, ν}$ .

First, we consider the last two terms in the boundary integral of (53). For convenience, we introduce a third-order tensor $M$ through

\begin{matrix} - \frac{2}{3} h^{2} (p^{(1)} \cdot η + (S_{0}^{T} ν) \cdot x^{(2)}) = \frac{2}{3} h^{2} {\bar{A} [t \otimes n - x^{(2)} \otimes ν] + (x^{(2)} \cdot ℝ^{(0) T} ν - t \cdot y} ℙ} : \nabla ξ \\ ≜ {(M [ν])}^{T} : \nabla ξ, \end{matrix}

(55)

where

t = B^{- 1} (p^{(1)} + B_{1} x^{(2)}), {(B_{1})}_{ij} = {\bar{A}}_{3 i α j} ν_{α} .

Furthermore, we introduce the decomposition

\nabla ξ = ξ_{, s} \otimes τ + ξ_{, ν} \otimes ν,

(56)

where $τ$ is the unit tangent vector and $ξ_{, s}$ is tangent derivative on $\partial Ω$ . Substituting (55) and (56) into the boundary integral in (54) and a simple integration by parts leads to

\int_{Ω} (W_{1} + W_{3} - \tilde{q} \cdot ξ) dA = \int_{\partial Ω} {\bar{p} - {(M [ν] τ)}_{, s}} \cdot ξ + {M [ν] ν} \cdot ξ_{, ν} ds,

(57)

where the terms before $ξ$ and $ξ_{, ν}$ are the generalized traction and moment on the edge, respectively. Corresponding to the different boundary conditions (e.g., clamped, pinned, simply supported), the conjugate variables ${x^{(0)}, \bar{p} - {(M [ν] τ)}_{, s}}$ and ${x_{, ν}^{(0)}, M [ν] ν}$ should be specified properly [40].

5. Examples

To show the validity of the shell theory established in the previous sections, we shall use it to study the deformations of two commonly used shell structures, i.e., the spherical shell and the circular cylindrical shell, in this section. The shells are supposed to be made of an incompressible neo-Hookean material. The approximate solutions obtained from the shell theory will be compared with the exact solutions to show their accuracies.

5.1. The spherical shell

For a spherical shell, it is convenient to use the spherical coordinates $(R, Θ, Φ)$ within the domain

A \leq R \leq B, 0 \leq Θ \leq π, 0 \leq Φ \leq 2 π,

(58)

which are related to rectangular Cartesian coordinates by

X_{1} = R \sin Θ \cos Φ, X_{2} = R \sin Θ \sin Φ, X_{3} = R \cos Θ .

(59)

In the notation of the current shell theory, we have the corresponding relations

θ^{1} = Θ, θ^{2} = Φ, Z = R - A .

The inner surface of the spherical shell is selected as the base surface $Ω$ . For any point on the base surface, the base vectors can be calculated by

\begin{matrix} g_{1} & = A (\cos Θ \cos Φ e_{1} + \cos Θ \sin Φ e_{2} - \sin Θ e_{3}) = A e_{Θ} = A^{2} g^{1}, \\ g_{2} & = A \sin Θ (- \sin Φ e_{1} + \cos Φ e_{2}) = A \sin Θ e_{Φ} = A^{2} \sin^{2} Θ g^{2}, \\ g_{3} & = \sin Θ \cos Φ e_{1} + \sin Θ \sin Φ e_{2} + \cos Θ e_{3} = e_{R} = g^{3} = n, \end{matrix}

(60)

where $e_{i} (i = 1, 2, 3)$ are the Cartesian base vectors and $e_{Θ}$ , $e_{Φ}$ , $e_{R}$ are the commonly adopted physical base vectors along the coordinates. From the equations of Gauss and Weingarten [5], we can define the following constants

Γ_{α β}^{σ} = \frac{\partial g_{α}}{\partial θ^{β}} \cdot g^{σ}, κ_{α}^{β} = - Γ_{3 α}^{β}, κ_{α β} = Γ_{α β}^{3} .

(61)

For the current example, the non-zero versions of the above constants are given by

\begin{matrix} Γ_{12}^{2} = Γ_{21}^{2} = \cot Θ, Γ_{22}^{1} = - \cos Θ \sin Θ, \\ κ_{11} = - A, κ_{22} = - A \sin^{2} Θ, κ_{1}^{1} = κ_{2}^{2} = - \frac{1}{A} . \end{matrix}

(62)

From (4) and (62), it can be obtained that

H = - \frac{1}{A}, K = \frac{1}{A^{2}} .

(63)

To satisfy the assumptions of the shell theory on geometry and curvature, we need to require that $2 h ∕ A = (B - A) ∕ A$ is small.

Through the manipulations introduced in Section 3, the shell equation (32) takes the following form in the current example

\begin{matrix} {\bar{S}}_{, α}^{α 1} + {\bar{S}}^{11} \cot Θ + \frac{1}{A} ({\bar{S}}^{13} + {\bar{S}}^{23}) - \cos Θ \sin Θ {\bar{S}}^{22} = - {\tilde{q}}^{1}, \\ {\bar{S}}_{, α}^{α 2} + (2 {\bar{S}}^{12} + {\bar{S}}^{21}) \cot Θ + \frac{1}{A} ({\bar{S}}^{13} + {\bar{S}}^{23}) = - {\tilde{q}}^{2}, \\ {\bar{S}}_{, α}^{α 3} + {\bar{S}}^{13} \cot Θ - A ({\bar{S}}^{11} + {\bar{S}}^{22} \sin^{2} Θ) = - {\tilde{q}}^{3}, \end{matrix}

(64)

where ${\bar{S}}^{α j} ≜ {\tilde{S}}^{α j} - h Ŝ^{α j} ∕ A$ .

Next, we consider the case that the shell structure undergoes the axisymmetric deformation

x (X) = r (R) e_{R}, A \leq R \leq B .

(65)

For an incompressible neo-Hookean material, the strain-energy density function is given by

Φ (F) = \frac{1}{2} μ (tr (F F^{T}) - 3),

(66)

then the nominal stress tensor can be calculated through

S (F) = μ F^{T} - p (Z) F^{- 1} .

(67)

On the inner and the outer surfaces of the spherical shell, we consider the traction boundary conditions

q^{-} = q^{-} e_{R} q^{+} = q^{+} e_{R} .

(68)

To apply the present shell theory, we consider the following series expansions of r(Z) and p(Z) in terms of Z

\begin{matrix} r (Z) = r_{0} + r_{1} Z + \frac{1}{2} r_{2} Z^{2} + \frac{1}{6} r_{3} Z^{3} + \frac{1}{24} r_{4} Z^{4} + O (Z^{5}), \\ p (Z) = p_{0} + p_{1} Z + \frac{1}{2} p_{2} Z^{2} + \frac{1}{6} p_{3} Z^{3} + \frac{1}{24} p_{4} Z^{4} + O (Z^{5}), \end{matrix}

(69)

where the subscripts are used instead of the superscripts for clarity. The non-zero components of $F_{ij}$ are given by

F_{11} = \frac{A^{2} r}{A + Z}, F_{22} = \frac{A^{2} r \sin^{2} Θ}{A + Z}, F_{33} = r^{'} .

(70)

Since we have

F = F_{ij} g^{i} \otimes g^{j} = F^{ij} g_{i} \otimes g_{j} = F_{e}^{11} e_{Θ} \otimes e_{Θ} + F_{e}^{22} e_{Φ} \otimes e_{Φ} + F_{e}^{33} e_{R} \otimes e_{R},

(71)

the non-zero components of $F^{ij}$ and $F_{e}^{ij}$ should be

\begin{matrix} F^{11} = \frac{r}{A^{2} (A + Z)}, F^{22} = \frac{r}{A^{2} \sin^{2} Θ (A + Z)}, F^{33} = r^{'}, \\ F_{e}^{11} = \frac{r}{A + Z}, F_{e}^{22} = \frac{r}{A + Z}, F_{e}^{33} = r^{'} . \end{matrix}

(72)

The non-zero components of the nominal stress tensor can then be calculated as

\begin{matrix} S^{11} = \frac{μ r}{A^{2} (A + Z)} - \frac{p (A + Z)}{A^{2} r}, S^{22} = S^{11} ∕ \sin^{2} Θ, S^{33} = μ r^{'} - \frac{p}{r^{'}} \cdot \end{matrix}

(73)

Through the series expansions of $F$ and $S$ , the recursion relations can then be obtained in the following form

\begin{matrix} r_{1} = \frac{A^{2}}{r_{0}^{2}}, p_{0} = \frac{μ A^{4}}{r_{0}^{4}} + \frac{A^{2} q^{-}}{r_{0}^{2}}, \\ r_{2} = \frac{2 A r_{0}^{3} - 2 A^{4}}{r_{0}^{5}}, p_{1} = - 2 μ (\frac{A^{6}}{r_{0}^{7}} - \frac{2 A^{3}}{r_{0}^{4}} + \frac{1}{r_{0}}), \\ r_{3} = \frac{2 r_{0}^{6} - 12 A^{3} r_{0}^{3} + 10 A^{6}}{r_{0}^{8}}, p_{2} = \frac{3 q^{-} r_{0}^{8} + 14 A^{8} μ - 28 A^{5} r_{0}^{3} μ + 14 A^{2} r_{0}^{6} μ}{r_{0}^{10}} . \end{matrix}

(74)

The non-trivial shell equation (64)₃ has the following form

- 2 A {\tilde{S}}^{11} + 2 h Ŝ^{11} = - {\tilde{q}}^{3},

(75)

which, by substituting all the recursion relations, provides an algebraic equation for r₀. As an illustrative example for comparison, we consider the case that $q^{-} = P {(r_{0} ∕ A)}^{2}$ and $q^{+} = 0$ . Although this is not a dead-loading case, we can follow the same procedure with only slight modifications to obtain the final equation

\begin{matrix} 3 ϱ λ_{a}^{13} + 12 λ_{a}^{6} - 12 λ_{a}^{12} + \frac{h}{A} (12 ϱ λ_{a}^{10} - 36 λ_{a}^{3} + 72 λ_{a}^{6} - 36 λ_{a}^{9}) \\ + \frac{h^{2}}{A^{2}} (- 12 ϱ λ_{a}^{7} + 24 ϱ λ_{a}^{10} + 176 - 400 λ_{a}^{3} + 304 λ_{a}^{6} - 80 λ_{a}^{9}) \\ + \frac{h^{3}}{A^{3}} (32 ϱ λ_{a}^{4} - 48 ϱ λ_{a}^{7} + 40 ϱ λ_{a}^{10}) = 0, \end{matrix}

(76)

where the scales $P = ϱ h μ ∕ A$ and $λ_{a} = r_{0} ∕ A$ have been used. Equation (76) can be solved through a regular perturbation method and the asymptotic solution is given by

r_{0} = r_{0}^{(0)} + r_{0}^{(1)} h + \frac{1}{2} r_{0}^{(2)} h^{2} + O (h^{3}),

(77)

where

\begin{matrix} ϱ λ_{a}^{{(0)}^{7}} = 4 λ_{a}^{{(0)}^{6}} - 4, r_{0}^{(1)} = \frac{7 - 6 λ_{a}^{{(0)}^{3}} - λ_{a}^{{(0)}^{6}}}{- 7 λ_{a}^{{(0)}^{2}} + λ_{a}^{{(0)}^{8}}}, \\ r_{0}^{(2)} = \frac{2}{3 A λ_{a}^{{(0)}^{5}} {(7 - λ_{a}^{{(0)}^{6}})}^{3}} (686 - 1666 λ_{a}^{{(0)}^{3}} + 525 λ_{a}^{{(0)}^{6}} + 672 λ_{a}^{{(0)}^{9}} \\ - 192 λ_{a}^{{(0)}^{12}} - 90 λ_{a}^{{(0)}^{15}} + 61 λ_{a}^{{(0)}^{18}} + 4 λ_{a}^{{(0)}^{21}}), \end{matrix}

(78)

and $λ_{a}^{(0)} = r_{0}^{(0)} ∕ A$ .

On the other hand, for the current example, the exact solution to the original 3D field equation takes the form [43]

r (R) = R {(1 + \frac{r_{0}^{3} - A^{3}}{R^{3}})}^{1 ∕ 3} .

(79)

Recall $R = A + Z$ , and then by Taylor’s expansion we have

r (Z) = r_{0} + \frac{A^{2}}{r_{0}^{2}} Z + \frac{A r_{0}^{3} - A^{4}}{r_{0}^{5}} Z^{2} + \frac{r_{0}^{6} - 6 A^{3} r_{0}^{3} + 5 A^{6}}{3 r_{0}^{8}} Z^{3} + O (Z^{4}),

(80)

which gives the same recursion relations for $r_{i} (i = 1, 2, 3)$ as the shell theory. For the above specific boundary conditions, one has (see [43])

P = \int_{λ_{b}}^{λ_{a}} \frac{{\hat{Φ}}^{'} (λ) d λ}{λ^{3} - 1},

(81)

where

\begin{matrix} λ_{b} = \sqrt[3]{{(\frac{A}{A + 2 h})}^{3} (λ_{a}^{3} - 1) + 1} = \frac{\sqrt[3]{{(A + 2 h)}^{3} - A^{3} + r_{0}^{3}}}{A + 2 h}, \end{matrix}

and for a neo-Hookean material

\begin{matrix} \hat{Φ} (λ) = \frac{1}{2} μ (2 λ^{2} + λ^{- 4} - 3) . \end{matrix}

The same scale $P = ϱ h μ ∕ A$ yields that

\frac{ϱ h}{2 A} = - \frac{1}{λ_{a}} - \frac{1}{4 λ_{a}^{4}} + \frac{1}{λ_{b}} + \frac{1}{4 λ_{b}^{4}} .

(82)

By substituting $λ_{a}$ and $λ_{b}$ into (82), we obtain an equation for r₀. After expanding r₀ like (77) and comparing the coefficients of $h^{n}$ ( $n = 0, 1, 2$ ) between two sides, we find that the first three terms are the same as (78), which implies that r₀ given in (77) and (78) is correct to $O (h^{2})$ . From the recursion relations (74), the other coefficients $r_{i}$ ( $i = 1, 2, 3$ ) can also be recovered. Since the exact solution (79) gives the same recursion relations as the shell theory, we conclude that all $r_{i}$ ( $i = 1, 2, 3$ ) have the correct terms up to $O (h^{2})$ . Therefore, for this problem, the shell theory produces $O (h^{2})$ -correct results for the deformation r(Z).

5.2. The circular cylindrical shell

For a circular cylindrical shell, the cylindrical coordinates $(R, Θ, X_{3})$ within the domain

A \leq R \leq B, 0 \leq Θ \leq 2 π, 0 \leq X_{3} \leq L

(83)

are adopted, which are related to the Cartesian coordinates by

X_{1} = R \cos Θ, X_{2} = R \sin Θ, X_{3} = X_{3} .

To apply the shell theory, we identify the cylindrical coordinates as

θ^{1} = Θ, θ^{2} = X_{3}, Z = R - A .

(84)

The inner surface of the cylinder is selected as the base surface $Ω$ . For any given point on $Ω$ , the base vectors are given by

\begin{matrix} g_{1} = A (- \sin Θ e_{1} + \cos Θ e_{2}) = A e_{Θ} = A^{2} g^{1}, \\ g_{2} = e_{3} = g^{2}, \\ g_{3} = \cos Θ e_{1} + \sin Θ e_{2} = e_{R} = g^{3} = n . \end{matrix}

(85)

The non-zero versions of the coefficients $Γ_{α β}^{γ}$ , $κ_{α β}$ , and $κ_{α}^{β}$ in the equations of Gauss and Weingarten are

κ_{11} = - A, κ_{1}^{1} = - \frac{1}{A},

(86)

which imply that $H = - 1 ∕ 2 A$ and K=0. For the current example, the shell equation (32) takes the following form

\begin{matrix} {\tilde{S}}_{, α}^{α 1} - \frac{h}{A} Ŝ_{, 1}^{11} + \frac{1}{A} ({\tilde{S}}^{13} - \frac{h}{A} Ŝ^{13}) = - {\tilde{q}}^{1}, \\ {\tilde{S}}_{, α}^{α 2} - \frac{h}{A} Ŝ_{, 1}^{12} = - {\tilde{q}}^{2}, \\ {\tilde{S}}_{, α}^{α 3} - \frac{h}{A} Ŝ_{, 1}^{13} - A {\tilde{S}}^{11} + h Ŝ^{11} = - {\tilde{q}}^{3} . \end{matrix}

(87)

We consider the axisymmetric deformation of the shell defined by

r = r (R), θ = Θ, x_{3} = λ X_{3},

(88)

where $(r, θ, x_{3})$ denote the coordinates in the current configuration. We suppose that the shell is made of an incompressible neo-Hookean material with the strain-energy density function (66). The non-zero components in the deformation gradient tensor and the nominal stress tensor are

\begin{matrix} F_{11} = \frac{A^{2} r}{A + Z}, F_{22} = λ, F_{33} = r^{'}, \\ F^{11} = \frac{r}{A^{2} (A + Z)}, F^{22} = λ, F^{33} = r^{'}, \\ F_{e}^{11} = \frac{r}{A + Z}, F_{e}^{22} = λ, F_{e}^{33} = r^{'}, \end{matrix}

(89)

and

S^{11} = \frac{μ r}{A^{2} (A + Z)} - \frac{p (A + Z)}{A^{2} r}, S^{22} = μ λ - \frac{p}{λ}, S^{33} = μ r^{'} - \frac{p}{r^{'}} .

(90)

Corresponding to the series expansions of $r = r (Z)$ and $p = p (Z)$ in (69), the non-zero components given in (89) and (90) can also be expanded. For the case $λ = 1$ , we can obtain the following recursive formulas for $r_{i + 1}$ and $p_{i}$ ( $i = 0, 1, 2$ )

\begin{matrix} r_{1} = \frac{A}{r_{0}}, r_{2} = \frac{r_{0}^{2} - A^{2}}{r_{0}^{3}}, r_{3} = \frac{3 A^{3} - 3 A r_{0}^{2}}{r_{0}^{5}}, \\ p_{0} = \frac{μ A^{2}}{r_{0}^{2}} + \frac{q^{-} A}{r_{0}}, p_{1} = \frac{- A^{4} μ + 2 A^{2} r_{0}^{2} μ - r_{0}^{4} μ}{A r_{0}^{4}}, \\ p_{2} = \frac{4 A^{6} μ - 7 A^{4} r_{0}^{2} μ + 2 A^{2} r_{0}^{4} μ + r_{0}^{6} μ}{A^{2} r_{0}^{6}} . \end{matrix}

(91)

Then, the non-trivial shell equation (87)₃ becomes

- A {\tilde{S}}^{11} + h Ŝ^{11} = - {\tilde{q}}^{3},

(92)

which provides an algebraic equation for $r_{0}$ . Here we suppose that the cylindrical shell undergoes inflation subject to the internal pressure, then the boundary conditions can be set as $q^{-} = q^{-} e_{R} = (P r_{0} ∕ A) e_{R}$ and $q^{+} = 0$ . The non-trivial shell equation (92) can be rewritten as

\begin{matrix} 3 ϱ λ_{a}^{8} + 6 λ_{a}^{4} - 6 λ_{a}^{8} + \frac{h^{3}}{A^{3}} (12 ϱ λ_{a}^{2} - 12 ϱ λ_{a}^{4}) \\ + \frac{h^{2}}{A^{2}} (- 6 ϱ λ_{a}^{4} + 6 ϱ λ_{a}^{6} + 36 - 64 λ_{a}^{2} + 36 λ_{a}^{4} - 8 λ_{a}^{8}) \\ + \frac{h}{A} (6 ϱ λ_{a}^{6} - 12 λ_{a}^{2} + 18 λ_{a}^{4} - 12 λ_{a}^{6} + 6 λ_{a}^{8}) = 0, \end{matrix}

(93)

where we still use the scale $P = ϱ h μ ∕ A$ and $λ_{a} = r_{0} ∕ A$ . We concentrate on the O(1) solution to $λ_{a}$ . For $ϱ$ close to 2, the exact solution (see (98)) goes to infinity. Here we consider $ϱ = 2 - ϖ$ where $ϖ < 2$ is a positive O(1) quantity. Following a similar procedure as that in the previous subsection, the first three terms in the perturbation expansion solution of r₀ are given by

\begin{matrix} λ_{a}^{(0)} = {(\frac{2}{2 - ϱ})}^{1 ∕ 4}, r_{0}^{(1)} = \frac{8 - 3 ϱ + (2 ϱ - 4) \sqrt{4 - 2 ϱ}}{2 \times 2^{3 ∕ 4} {(2 - ϱ)}^{5 ∕ 4}}, \\ r_{0}^{(2)} = \frac{1}{24 \times 2^{3 ∕ 4} A {(2 - ϱ)}^{9 ∕ 4}} {- 24 ϱ^{3} + 159 ϱ^{2} - 320 ϱ + 256 - 44 \sqrt{2} ϱ^{2} \sqrt{2 - ϱ} + (152 ϱ - 128) \sqrt{4 - 2 ϱ}} . \end{matrix}

(94)

On the other hand, the current example has the following exact solution [43]

r (R) = \sqrt{R^{2} - A^{2} + r_{0}^{2}},

(95)

with $R = A + Z$ . It can be verified that, if this exact solution is expanded as a Taylor series, the recursion relations for $r_{i}$ ( $i = 1, 2, 3$ ) are the same as the shell theory. For the above-used specific boundary conditions, one has [43]

P = \int_{λ_{b}}^{λ_{a}} \frac{{\hat{Φ}}^{'} (λ) d λ}{λ^{2} - 1},

(96)

where

\begin{matrix} λ_{b} = \sqrt{\frac{A^{2}}{{(A + 2 h)}^{2}} (λ_{a}^{2} - 1) + 1} = \frac{\sqrt{{(A + 2 h)}^{2} - A^{2} + r_{0}^{2}}}{A + 2 h}, \\ \hat{Φ} (λ) = \frac{1}{2} μ (λ^{2} + λ^{- 2} - 2) \cdot \end{matrix}

(97)

Using the same scales as before, we have

\frac{ϱ h}{A} = \ln λ_{a} - \frac{1}{2 λ_{a}^{2}} - \ln λ_{b} + \frac{1}{2 λ_{b}^{2}} .

(98)

From (96)–(98), one can seek a perturbation expansion solution of the form (77), and it turns out that the first three terms are just the same as (94). It can also be verified that the other coefficients $r_{i}$ ( $i = 1, 2, 3$ ) are all correct to $O (h^{2})$ . Hence, for this internally pressurized case, the shell theory also produces the $O (h^{2})$ -correct results for the deformation.

From the illustrative examples given in this section, we can draw a conclusion that the current shell theory can provide asymptotically correct solutions. Furthermore, the theory can deal with not only the dead-loading cases, but also the pressurized cases. We expect that the shell theory introduced in the current paper will have more applications in practical cases.

6. Conclusions

In this article, a finite-strain shell theory for incompressible hyperelastic materials has been proposed, which is consistent with the 3D variational formulation and is applicable in general loading conditions. Starting from the 3D governing system and through the series expansion method, the constraint equation of incompressibility can be naturally taken into account without any special hypotheses. Once the solution to the shell equation is obtained, the 3D displacement and stress fields in the thin shell structures can also be recovered. The illustrative examples have revealed that the current shell theory can provide asymptotically $O (h^{2})$ correct solutions. Furthermore, the weak formulation of the 2D shell equation has also been derived, which facilitates the future numerical calculations. In our future work, the current shell theory will be further extended for studying the growth-induced pattern formations from the morphogenesis of soft biological tissues.

To end this paper, we discuss some advantages and disadvantages of this shell theory and the approach in the framework of derived shell theories. First, compared with shell theories for incompressible hyperelastic materials based on ad hoc hypotheses, the advantage of the present shell theory is that one does not need to worry about the validity of the hypotheses for a particular application. In addition, it can deal with the case of pure shear traction on the top/bottom surface, for which most (if not all) hypotheses-based shell theory will not work. The authors are not aware of works on using asymptotic methods to derive shell models for general incompressible hyperelastic materials. For compressible Saint-Venant materials, either membrane-type or bending-type shell models were deduced by using asymptotic methods for certain particular scales (in terms of the thickness) of applied tractions (e.g., both shear traction and normal traction of the order of the thickness; see Hamdouni and Millet [16, 17]). The $Γ$ -convergence approach also depends on specific a priori scalings, which also yields either membrane-type or bending-type shell models. The shell models for scalings other than those specified may not be deduced by such two approaches (e.g., shear traction is of the order of thickness and normal traction is of the order of thickness square). From a physical point of view, applied data is external, which should not be directly related to the shell thickness. In addition, for a quasi-static process, the amplitude of applied data changes and cannot be a fixed scale of the thickness. Thus, despite the mathematical elegance of these two approaches, they still leave a number of open issues in constructing shell models. The present approach does not need a priori scalings on applied data and displacements, so it does not lead to the aforementioned defects associated with them. One important aspect is that it finishes a single shell model that incorporates both stretching and bending effects. Of course, the present theory and approach also have certain disadvantages. Owing to the highly nonlinear nature of incompressible hyperelastic materials, the recursive relations, obtainable, but are rather complicated. In particular, for the second expansion coefficients one needs to solve three nonlinear algebraic equations and one linear algebraic equation (for a general strain energy function), which may have to be done numerically in practice. A finite element modeling (FEM) scheme for this shell theory is implementable (the aforementioned four algebraic equations can be treated as 2D constraints), but putting the lengthy recursive relations into a program requires very tedious work (unfortunately Matlab still does not have a strong symbolic capacity at this stage; otherwise the programming is not difficult). A disadvantage of this approach is that it cannot handle the boundary layer region and the resulting equations are only good for the interior region (existing shell theory for incompressible hyperelastic materials have the same defects). However, the edge boundary conditions proposed agree with the 3D weak formulation (for which 3D FEM is based) to proper order, which gives confidence on using them. In short, whereas the present shell theory for incompressible hyperelastic materials has a few of advantages over existing ones, efforts are still needed to ease the programming of the FEM scheme. Numerical computations based on this shell theory for particular applications on biological materials will be for future studies.

Footnotes

Acknowledgements

The authors thank Professor P.G. Ciarlet for some valuable advice.

Funding

The author disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This paper was supported by a GRF grant (project no. CityU 11303015) from the Research Grants Council of Hong Kong SAR, China. Jiong Wang is supported by a grant from the National Nature Science Foundation of China (project no. 11602086) and the Guangdong Nature Science Funds for Distinguished Young Scholar (project no. 2015A030306009).

ORCID iD

Yuanyou Li .

References

Love

AEH.

On the small free vibrations and deformation of thin elastic shells. Philos Trans R Soc Lond Ser A 1888; 179: 491–546.

Timoshenko

Woinowsky-Krieger

Theory of Plates and Shells, Vol. 2. New York: McGraw-Hill, 1959.

Koiter

Simmonds

JG.

Foundations of Shell Theory. New York: Springer, 1973.

Ciarlet

PG.

Mathematical Elasticity, Vol. III: Theory of Shells. Amsterdam: North-Holland, 2000.

Ciarlet

PG.

An introduction to differential geometry with applications to elasticity. J Elast 2005; 78: 1–215.

Pietraszkiewicz

Gorski

Shell Structures: Theory and Applications, Vol. 3. Boca Raton, FL: CRC Press, 2013.

Kirchhoff

GR.

ber das Gleichgewicht und Bewegungen einer elastischen Scheibe. J reine angew Math 1850; 40: 51–88.

Reissner

On the Theory of Thin Elastic Shells. Ann Arbor, MI: JW Edwards, 1949, pp. 231–247.

Naghdi

PM.

On the theory of thin elastic shells. Q Appl Math 1957; 14: 369–380.

10.

Koiter

. A consistent first approximation in the general theory of thin elastic shells. In: Theory of Thin Elastic Shells. Amsterdam: North-Holland, 1960, pp. 12–33.

11.

Gol’denveizer

AL.

Derivation of an approximate theory of shells by means of asymptotic integration of the equations of the theory of elasticity. J Appl Math Mech 1963; 27: 903–924.

12.

Kaplunov

Kossovitch

Nolde

EV.

Dynamics of thin walled elastic bodies. New York: Academic Press, 1998.

13.

Aghalovyan

Asymptotic theory of anisotropic plates and shells. Singapore: World Scientific, 2015.

14.

Bazarenko

Israilevich

VI.

Asymptotic solution of the elasticity problem for a hollow, finite length, thin cylinder. J Appl Math Mech 1965; 29: 1219–1238.

15.

Vilenskaia

Israilevich

VI.

Asymptotic behavior of the solution of the problem of elasticity for a thin spherical shell. J Appl Math Mech 1967; 30: 342–361.

16.

Hamdouni

Millet

Classification of thin shell models deduced from the nonlinear three-dimensional elasticity. Part I: the shallow shells. Arch Mech 2003; 55: 135–176.

17.

Hamdouni

Millet

Classification of thin shell models deduced from the nonlinear three-dimensional elasticity. Part II: the strongly curved shells. Arch Mech 2003; 55: 177–220.

18.

Kienzler

On consistent plate theories. Arch Appl Mech 2002; 72: 229–247.

19.

Kienzler

Schneider

Consistent theories of isotropic and anisotropic plates. J Theoret Appl Mech 2012; 50: 755–768.

20.

Ciarlet

Destuynder

A justification of a nonlinear model in plate theory. Comput Meth Appl Mech Eng 1979; 17: 227–258.

21.

Destuynder

Sur une Justification des Modeles de Plaques et de Coques par les Methodes Asymptotiques. PhD thesis, Universite Pierre et Marie Curie, Paris, 1980.

22.

Le Dret

Raoult

. The nonlinear membrane model as a variational limit of nonlinear three-dimensional elasticity. J Math Pures Appl 1995; 74: 549–578.

23.

Le Dret

Raoult

. The membrane shell model in nonlinear elasticity: a variational asymptotic derivation. J Nonlin Sci 1996; 6: 59–84.

24.

Friesecke

James

Mora

. Derivation of nonlinear bending theory for shells from three-dimensional nonlinear elasticity by Γ -convergence. C R Math. 2003; 336: 697–702.

25.

Simo

Rifai

Fox

DD.

On a stress resultant geometrically exact shell model. Part IV: Variable thickness shells with through-the thickness stretching. Comput Meth Appl Mech Eng 1990; 81: 91–126.

26.

Büchter

Ramm

Roehl

Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept. Int J Numer Meth Eng 1994; 37: 2551–2568.

27.

Reddy

Arciniega

Shear deformation plate and shell theories: from stavsky to present. Mech Adv Mater Struct 2004; 11: 535–582.

28.

Steigmann

DJ.

Extension of Koiter’s linear shell theory to materials exhibiting arbitrary symmetry. Int J Eng Sci 2012; 51: 216–232.

29.

Steigmann

DJ.

Koiter’s shell theory from the perspective of three-dimensional nonlinear elasticity. J Elast 2013; 111: 91–107.

30.

Makowski

Stumpf

Finite strains and rotations in shells. In: Finite Rotations in Structural Mechanics. Berlin: Springer, 1986, pp. 175–194.

31.

Itskov

A generalized orthotropic hyperelastic material model with application to incompressible shells. Int J Numer Meth Eng 2001; 50: 1777–1799.

32.

Chapelle

Mardare

Münch

Asymptotic considerations shedding light on incompressible shell models. J Elast 2004; 76: 199–246.

33.

Chermisi

The von Kármán theory for incompressible elastic shells. Calculus Variations Partial Differ Eqns 2013; 48: 185–209.

34.

Naghdi

. The theory of plates and shells. In: Handbuch der Physik, Vol. VIa/2. Heidelberg: Springer, 1972, pp. 425–640.

35.

Zhilin

PA.

Mechanics of deformable directed surfaces. Int J Solids Struct 1976; 12: 635–648.

36.

Rubin

MB.

Cosserat theories: shells, rods and points. Dordrecht: Springer Science and Business Media, 2000.

37.

Altenbach

Eremeyev

VA.

On generalized Cosserat-type theories of plates and shells: a short review and bibliography. Arch Appl Mech 2010; 80: 73–92.

38.

Altenbach

Eremeyev

Thin-walled structural elements: classification, classical and advanced theories, new applications. In: Shell-like Structures: Advanced Theories and Applications. Cham: Springer, 2017, pp. 1–62.

39.

Libai

Simmonds

JG.

The Nonlinear Theory of Elastic Shells. Cambridge: Cambridge University Press, 1998.

40.

Dai

Song

ZL.

On a consistent finite-strain plate theory based on a 3-D energy principle. Proc R Soc Lond Ser A 2014; 470: 2171.

41.

Song

Dai

HH.

On a consistent finite-strain shell theory based on 3-D nonlinear elasticity. Int J Solids Struct 2016; 97: 137–149.

42.

Wang

Song

Dai

HH.

On a consistent finite-strain plate theory for incompressible hyperelastic materials. Int J Solids Struct 2016; 78: 101–109.

43.

Ogden

Non-linear Elastic Deformations. New York: Ellis Horwood, 1984.

44.

Chadwick

Continuum Mechanics: Concise Theory and Problems (2nd edn). New York: Dover Publications, 1999.