Stimuli-responsive shell theory

Abstract

Soft matter mechanics generally involve finite deformations and instabilities of structures in response to a wide range of mechanical and non-mechanical stimuli. Modeling plates and shells is generally a challenge due to their geometrically nonlinear response to loads; however, non-mechanical loads further complicate matters as it is often not clear how they modify the shell’s energy functional. In this work, we demonstrate how to form a mechanical interpretation of these non-mechanical stimuli, in which the standard shell strain energy can be augmented with potentials corresponding to how a non-mechanical stimulus acts to change the shell’s area and curvature via the natural stretch and curvature. As a result, the effect of non-mechanical stimuli to deform shells is transformed into effective external loadings, and this framework allows for the application of analytical and computational tools that are standard within the mechanics community. Furthermore, we generalize the effect of mass change during biological growth to account for its effect on the stress constitution. The theory is formulated based on a standard, stress-free reference configuration which is known a priori, meaning it can be physically observed, and only requires the solution of a single-field equation, the standard mechanical momentum or equilibrium equation, despite capturing the effects of non-mechanical stimuli. We validate the performance of this model by several benchmark problems, and finally, we apply it to complex examples, including the snapping of the Venus flytrap, leaf growth, and the buckling of electrically active polymer plates. Overall, we expect that mechanicians and non-mechanicians alike can use the approach presented here to quickly modify existing computational tools with effective external loadings calculated in this novel theory to study how various types of non-mechanical stimuli impact the mechanics and physics of thin shell structures.

1. Introduction

Soft matter deforms in response to a variety of stimuli, which may include biochemical signaling during growth [1–11], voltages applied to electroactive polymers [12–21], residual stresses that build up during three-dimensional (3D) printing [22–26], externally applied magnetic fields [27–34], interactions with surrounding granular matter [35–40], and swelling with fluids [41–53]. Modeling such interactions and deformations presents a variety of challenges. One may write down the field equations corresponding to a given stimulus; however, simultaneously solving multiple field equations to capture tightly coupled multiphysics phenomena generally requires significant computational expense, as well as complex and specialized numerical formulations. Adding to that challenge are the finite deformations common in soft materials and the geometric nonlinearity and instabilities of slender structures. Taken together, this presents a formidable challenge for modeling the behavior of stimuli-responsive materials.

We know that the equilibrium configurations of a structure, and the stability of these equilibria, can be determined through the analysis of its total potential energy. The total potential energy is composed of the strain energy stored in the deforming body, along with any potentials due to external mechanical forces and moments that are acting on it. Here, we are concerned with thin shells—curved structures whose thickness is small compared to their other two spatial dimensions. Koiter’s foundational contributions to shell theory throughout the middle of the 20th century showed that the strain energy stored within a deforming shell may be decomposed into separate contributions from stretching and bending the shell’s middle surface [54–58]. The principle of minimum potential energy ensures that equilibrium configurations of a shell are found by extremizing the total potential energy, a process aided by the calculus of variations, where equilbria are found when the first variation of the total potential energy vanishes. The stability of these equilibrium configurations is then determined by the character of high-order variations of the total potential energy—it is commonly sufficient to check the second variation, with stable equilibria having a second variation that is positive-definite. A common challenge is encountered when a shell is subjected to non-mechanical stimuli, for instance, thermal, swelling, electrical, or growth—namely, how do we write the total potential energy for shells subjected to these stimuli? These problems are particularly important in the study of soft matter, where non-mechanical stimuli can cause extremely large deformations and may result in elastic instabilities.

For multi-physics systems, many theories have been actively developed to describe the behavior of 3D solid bodies subject to non-mechanical stimuli. They can be classified into three types: diffusion models [59–68], continuum elastic models [69–75], and mixture models [76–79]. Despite this ongoing research, a relatively small number of attempts have been presented to apply those modeling theories for thin bodies [80–88] and to develop modeling theories of plates and shells subject to non-mechanical stimuli [89–95]. Meanwhile, great strides were made in modeling stimuli-responsive plates through the introduction of the so-called non-Euclidean Plate (NEP) theory [90, 96–100]. In this approach, the strain energy of an elastic plate is calculated by considering the difference between the deformed configuration of the plate and the “stress-free” configuration of the plate. This intermediate stress-free configuration is distinct from the initial configuration and is in general incompatible—meaning it cannot be physically observed. This presents a practical challenge for modeling stimuli-responsive plates and shells: the intermediate configuration is very complex and difficult to conceptualize even for plates exposed to simple stimuli. As we will show, this work recasts NEP theory in a more traditional mechanical framework, that is, strains are measured with respect to the reference configuration that is known a priori, rather than the intermediate configuration. However, we are not altogether done with the intermediate configuration. In what we will put forth as a stimuli-responsive shell (SRS) theory, the intermediate configuration enters in the potential of the effective external loads, which are conjugate to the non-mechanical stimuli that act to stretch and bend the shell. In this framework, it becomes far easier to conceptualize the intermediate configuration as it is only necessary to consider how the stimuli deform a beam, whose intermediate configuration is compatible, rather than a shell, whose intermediate configuration is not.

At every point on the shell, one may imagine cutting two orthogonal beams from the surface as shown in Figure 1. Removed from the geometric constraints of the shell’s surface, these beams are free to deform in response to stimuli that act to stretch or bend them. A non-mechanical stimulus may cause the neutral axis of these beams to stretch or bend, deformations that we refer to as the natural stretch and natural curvature, respectively. This way of conceptualizing the effect of a stimulus is not merely a gedanken experiment, but can in fact be used to experimentally measure the natural stretch and natural curvature that a shell experiences [48, 51, 101] and is analogous to common biological measurements of residual stress within soft tissues [102–109]. As such, these quantities will form the basis for the SRS theory that follows. Beyond stretching and bending, non-mechanical stimuli may add mass to the shell—an effect that occurs in problems involving biological growth [110–113] and when a soft material is swollen by a macromolecule, such as a polymer chain [47, 114, 115]. As we will show, accounting for the addition of material to a shell requires a subtle accounting of how that mass changes the shell’s free energy, and we put forth a generalized framework for adding mass—either similar or dissimilar—to the existing shell.

Figure 1.

Schematic of two orthogonal beams ( $b1$ and $b2$ in the figure) cut from a surface. (a) Beam cut from the initial configuration with initial length ( $L_{b1}$ and $L_{b2}$ ) and radius of curvature ( $R_{b1}$ and $R_{b2}$ ). (b) Beam cut from the deformed configuration upon application of a stimulus, which shows the emergence of natural stretch ( ${\mathrm{\Lambda }}_{b1}$ and ${\mathrm{\Lambda }}_{b2}$ ) and natural curvature ( ${\kappa }_{b1}$ and ${\kappa }_{b2}$ ) under the assumption that the natural stretch and curvature are homogeneous along each beam cut, which itself is sufficiently narrow to be free of any constraints. The dotted shapes of the isolated beams denote the original beam cut before the geometric constraints are released.

In this work, we put forth a general shell theory to account for non-mechanical stimuli. We consider three generic forms of stimuli: (1) one that adds mass to the shell, (2) one that acts to increase the area of the shell (locally or globally), and (3) one that acts to change the curvature of the shell (locally or globally). We will show in Section 3 on thermodynamics and conservation of mass that when mass is added to an elastic body, it alters the internal energy density and the state of stress. As a result, the second Piola–Kirchhoff stress ${\mathrm{\Sigma }}^{\mathrm{ij}}$ —which is typically found by differentiating the strain energy density per unit volume with respect to the metric tensor—is modified to account for a change in mass $\mathrm{\Delta }{M}_g$ relative to the shell’s original volume $V_0$ according to

$${\mathrm{\Sigma }}^{\mathrm{ij}}=\underset{\mathrm{Classical}}{\underbrace{2\frac{\partial {\mathrm{\Psi }}_0}{\partial g_{\mathrm{ij}}}}}+\underset{\mathrm{Added}\mathrm{Mass}}{\underbrace{2\frac{\mathrm{\Delta }{M}_g}{V_0}\frac{\partial {\psi }_g}{\partial g_{\mathrm{ij}}}}},$$

(1)

where ${\mathrm{\Psi }}_0$ is the energy density form with no mass addition using the material properties of the initial body, and $g_{\mathrm{ij}}$ is the 3D metric tensor of the current (deformed) configuration. On $\psi$ denoting the free energy, the subscript $g$ indicates that the added mass may or may not be the same material as the mass of the original shell—a different material will also have a different free energy, and therefore different mechanical properties. In general, when mass is added locally to a shell, for instance, through biological growth, this added mass will increase the shell’s state of stress. Furthermore, if the material added during growth is different from what is already present in the body, this will change the material properties [111, 116, 117], which has significant ramifications for deformation and material stability. We note that the formulation we present is more general than previous work [72], which always gives the same stress amplification for all types of stimuli due to its geometric dependence on changes in the area or volume of the material.

The primary contribution of this paper is the mechanical interpretation of non-mechanical stimuli, which will be shown in Section 5. Therefore, one may take the strain energy $\mathcal{U}$ of a plate or shell (e.g. Föppl–von Kàrmàn; Donnell–Mushtari–Vlasov; Koiter) and add to it a potential corresponding to how a non-mechanical stimulus acts to change the shell’s area, $P_{\mathrm{ext}}^s$ , and a potential corresponding to a stimulus that acts to change the shell’s curvature, $P_{\mathrm{ext}}^b$ , that is,

$$\mathcal{V}=\mathcal{U}-\int P_{\mathrm{ext}}^s+P_{\mathrm{ext}}^{b}d\stackrel{°}{a},$$

(2)

as is standard for work from mechanical loads, with the initial area element $d\stackrel{°}{a}$ . Non-mechanical stimuli that act to stretch an area of the shell are characterized by a natural stretch $\mathrm{\Lambda }$ , which is work conjugate to the stretching strains of the shell’s middle surface which correspond to an intermediate stress-free configuration, ${\overline{E}}_{\alpha \beta }$ . The natural stretch will be a scalar $\mathrm{\Lambda }$ if it is uniformly acting over the shell, but in general, it will appear as a second-order tensor ${\mathrm{\Lambda }}_{\gamma \delta }$ , and as we will show, the potential of natural stretch is

$$P_{\mathrm{ext}}^s=h{\AA }^{\alpha \beta \gamma \delta }{\overline{E}}_{\alpha \beta }{\mathrm{\Lambda }}_{\gamma \delta }.$$

(3a)

As you may expect, the thickness $h$ enters linearly, and the strain tensor and natural stretch are associated with the elasticity tensor ${\AA }^{\alpha \beta \gamma \delta }$ in the reference configuration, as neglecting changes in the energy due to differences in metric tensors corresponds to the largest principal extension in the stress-free state from the initial body. An example of a non-mechanical stimulus that imparts a natural stretch would be a shell that is growing uniformly through its thickness. If a shell is experiencing differential growth, where the outer surface of the shell is growing more or less than the inner surface of the shell, this will generate a natural curvature $\kappa$ that is a generalization of Timoshenko’s heating of a bimetallic strip. This stimulus acts to locally change the curvature of the shell and is in general a second-order tensor ${\kappa }_{\gamma \delta }$ . The natural curvature is work conjugate to the curvature strains of the shell’s middle surface which correspond to an intermediate stress-free configuration, ${\overline{K}}_{\alpha \beta }$ , such that the potential of natural curvature is

$$P_{\mathrm{ext}}^b=\frac{h^3}{12}{\AA }^{\alpha \beta \gamma \delta }{\overline{K}}_{\alpha \beta }{\kappa }_{\gamma \delta },$$

(3b)

where, as you may expect, the potential of natural curvature scales as the cube of the shell’s thickness.

Our intent is to make analytical and computational modeling of shells exposed to non-mechanical stimuli generic and approachable with the tools and techniques common to theoretical and applied mechanics. As such, we also describe at length how to calculate the natural stretch ${\mathrm{\Lambda }}_{\delta \gamma }$ and natural curvature ${\kappa }_{\delta \gamma }$ from experimental measurements. The scalar forms of $\mathrm{\Lambda }$ and $\kappa$ will be somewhat familiar to readers well-versed in Timoshenko’s bimetallic beam theory. As we will show for a bilayer shell, both $\mathrm{\Lambda }$ and $\kappa$ depend on the relative expansion factor of each layer, the modulus ratio, the thickness ratio, and the stress amplification due to mass addition. The expansion factor is perhaps the least familiar, but it simply involves measuring how much the distance between material points has changed in response to a non-mechanical stimulus.

As is common with shell theories, our efforts focus on projecting the effects of these stimuli onto the middle surface of the shell. While we discuss $\mathrm{\Lambda }$ and $\kappa$ in the context of bilayer shells, this SRS theory can be applied to a shell with an arbitrary number of layers through the thickness. In Section 4, we describe how to model a generic multilayer shell as an equivalent monolayer shell, where the stresses that result from stretching and bending the shell are projected onto the shell’s middle surface. We utilize the nonlinear kinematics of Koiter’s shell theory throughout this paper, as it allows us to validate our model on a variety of well-studied problems involving stimuli-responsive shells. In Section 5, the SRS theory is developed by mapping non-mechanical stimuli to effective mechanical stretching and curvature potentials. Using a variational principle, we derive the first variation of the total potential energy, which may serve as the weak form to be solved using the finite element method, as well as the second variation of the total potential energy, which may be used to check the stability of a given equilibrium configuration. The performance of this proposed shell theory is illustrated by several benchmark problems and verified by the available analytic solutions in Section 6. A benefit of mapping non-mechanical stimuli to effective mechanical potentials is that the resulting total potential energy (i.e. equation (2)) is amenable to linear stability analysis, perturbation methods, and other tools from analysis and applied mathematics. We demonstrate this in Section 7, where we look at several example problems analytically and computationally. These include the growth of a leaf (which involves natural stretching), the rapid closure of the Venus flytrap’s leaves (which involves natural curvature), and the buckling of a plate made of a dielectric elastomer (which involves linear stability analysis). This SRS theory utilizes standard Kirchhoff–Love kinematics while accounting for nonlinearity. Because such details have been widely presented in previous works [54–58], we provide relevant background information in Appendices 1–3 while focusing the discussion in the main text on how to account for non-mechanical stimuli as equivalent mechanical stretching and curvature potentials, and how those potentials can be evaluated using common experimental measurements.

A brief note about notation. In what follows, the reader will encounter what Cartan [118] referred to as a “debauch of indices.” To make matters worse, in order to describe a deforming shell we must keep track of several curved configurations simultaneously. We begin with a reference configuration (which may also be referred to as the original, initial, or undeformed configuration). Quantities associated with the reference configuration will be denoted with an overcircle, for example, ${\stackrel{°}{a}}_{\alpha \beta }$ . Quantities associated with the current or deformed configuration of the shell will have no decoration above them, for example, $a_{\alpha \beta }$ . The differential geometry of surfaces tells us that shells often cannot directly adopt the shape prescribed by a given stimulus, meaning that the deformed configuration—the one that will be observed in practice—will possess a residual state of stress. This motivates the consideration of a third configuration that describes the stress-free (or natural) state of the shell, one that may be impossible to physically realize but nonetheless provides an essential means for performing calculations. Terms associated with this surface will be denoted with an overbar, for example, ${\overline{a}}_{\alpha \beta }$ . To make matters worse, a general surface containing multiple layers through its thickness requires us to keep track of the stress-free configuration of each layer, which we will denote with an over-tilde, for example, ${\tilde{a}}_{\alpha \beta }$ . Shells by definition are not flat and as such require curvilinear coordinates, as opposed to Cartesian coordinates, to describe material points in space. Switching between different shell surfaces, for example, ${\stackrel{°}{a}}_{\alpha \beta }\to {\overline{a}}_{\alpha \beta }\to a_{\alpha \beta }$ , requires us to change our basis, which requires to keep track of how components transform under a change of basis, with covariant components (lower indices) changing in the same way as the basis, for example, a gradient operator, and contravariant components (upper indices) changing inversely with a change in basis, for example, a velocity vector. We employ Einstein’s summation convention, where we sum over repeated indices—a contraction between a covariant component and a contravariant component. For a primer on the underlying differential geometry, the reader is referred to Needham [119]; for a brief on tensor analysis, the reader is referred to Simmonds [120]; and for the governing shell theory, the reader is encouraged to examine Koiter’s [54–58] work and Niordson’s [121] book. In this paper, Latin indices take 1,2,3, and the Greek indices are restricted to 1,2.

2. Kinematics with the multiplicative decomposition of deformation gradient

There has recently been significant interest in using the multiplicative decomposition of the deformation gradient as a means for accounting for the effect of non-mechanical stimuli, such as heating, swelling, and biological growth, on the mechanical response of structures and materials [59, 69, 71, 110, 112, 122, 123]. In this approach, the 3D deformation gradient, which is a one-to-one mapping of a material point $\stackrel{°}{\mathit{x}}$ in the initial configuration $\stackrel{°}{\mathrm{\Omega }}$ to a material point $\mathit{x}$ in the current (deformed) configuration $\mathrm{\Omega }$ , is composed as $\mathit{F}=\partial \mathit{x}/\partial \stackrel{°}{\mathit{x}}={\mathit{F}}_m{\mathit{F}}_n$ , where ${\mathit{F}}_m$ and ${\mathit{F}}_n$ are the mechanical (elastic) and non-mechanical parts, respectively. The non-mechanical part of the deformation gradient constructs a conceptual intermediate (stress-free) state as the mapping of the initial material points $\stackrel{°}{\mathit{x}}$ to the intermediate material points $\tilde{\mathit{x}}$ . Here, we use the term “state” rather than using “configuration” as the intermediate configuration $\tilde{\mathrm{\Omega }}$ cannot be realized if the deformation gradient ${\mathit{F}}_n$ does not uniquely define a one-to-one mapping, that is, incompatible. Meanwhile, the elastic part of the deformation gradient describes elastic deformations that satisfy the linear momentum balance on the deformed configuration, for which stresses are calculated via a relation between the current configuration and its corresponding intermediate state at local points. Then, the internal term of mechanical potential energy is calculated by means of the elastic part of the deformation gradient, that is, deformation from the intermediate stress-free state.

To apply this approach to the thin-shell kinematics, we begin with the 3D covariant tangent vector that is obtained via the parametric derivative with respect to a material point on the surface of a 3D body, for example, ${\mathit{g}}_i=\partial \mathit{x}/\partial {\xi }^i$ for the current configuration, based on the differential geometry with the general coordinate system using the curvilinear coordinates ${\xi }^i$ embedded in the body. The covariant components of the 3D metric tensor, including information about the distance between the material points, are calculated as $g_{\mathrm{ij}}={\mathit{g}}_i·{\mathit{g}}_j$ for the current configuration, from which the dual vector is ${\mathit{g}}^i=g^{\mathrm{ij}}{\mathit{g}}_j$ with $g^{\mathrm{ij}}=[g_{\mathrm{ij}}{]}^{-1}$ such that ${\mathit{g}}^i·{\mathit{g}}_j={\delta }_j^i$ is the Kronecker delta. All differential geometry calculations can be also applied for the initial configurations, denoted with the over-circle notation, for example, $\stackrel{°}{\mathit{x}}$ , ${\stackrel{°}{\mathit{g}}}_i$ , ${\stackrel{°}{\mathit{g}}}^i$ , ${\stackrel{°}{g}}_{\mathrm{ij}}$ , and ${\stackrel{°}{g}}^{\mathrm{ij}}$ . Meanwhile, the 3D metric tensor of the intermediate stress-free state can be calculated as ${\tilde{g}}_{\mathrm{ij}}={\tilde{\mathit{g}}}_i·{\tilde{\mathit{g}}}_j$ , where ${\tilde{\mathit{g}}}_i={\mathit{F}}_n{\stackrel{°}{\mathit{g}}}_i$ based on the definition of the non-mechanical part of deformation gradient as ${\mathit{F}}_n=\partial \tilde{\mathit{x}}/\partial \stackrel{°}{\mathit{x}}={\tilde{\mathit{g}}}_i\otimes {\stackrel{°}{\mathit{g}}}^i$ . This indicates that the stress-free 3D metric tensor is determined by the initial geometry and the specific form of ${\mathit{F}}_n$ which depends on the stimulus type and profile. The dual vector is calculated as ${\tilde{\mathit{g}}}^i={\tilde{g}}^{\mathrm{ij}}{\tilde{\mathit{g}}}_j$ with ${\tilde{g}}^{\mathrm{ij}}=[{\tilde{g}}_{\mathrm{ij}}{]}^{-1}$ such that ${\tilde{\mathit{g}}}^i·{\tilde{\mathit{g}}}_j={\delta }_j^i$ . Here, we use the over-tilde notation for the differential geometry calculation in terms of the intermediate state, for example, $\tilde{\mathit{x}}$ , ${\tilde{\mathit{g}}}_i$ , ${\tilde{\mathit{g}}}^i$ , ${\tilde{g}}_{\mathrm{ij}}$ , and ${\tilde{g}}^{\mathrm{ij}}$ .

As the elastic deformation is represented through the mechanical part of the deformation gradient, the elastic Green–Lagrange strain tensor is obtained as

$$\mathit{ϵ}=\frac{1}{2}({\mathit{C}}_m-\mathit{I})={\mathit{ϵ}}_{\mathrm{ij}}{\tilde{\mathit{g}}}^i\otimes {\tilde{\mathit{g}}}^j,{ϵ}_{\mathrm{ij}}=\frac{1}{2}(g_{\mathrm{ij}}-{\tilde{g}}_{\mathrm{ij}})$$

(4)

with the mechanical (elastic) right Cauchy–Green tensor ${\mathit{C}}_m={\mathit{F}}_m^T{\mathit{F}}_m=g_{\mathrm{ij}}{\tilde{\mathit{g}}}^i\otimes {\tilde{\mathit{g}}}^j$ and the second-rank identity tensor $\mathit{I}$ .

For thin surfaces, as the transverse shear deformations are negligible, the Kirchhoff–Love shell assumption that the normal vector of the mid-surface remains normal to the mid-surface after deformations is applicable to describe the geometry of thin shells [124]. Hence, Kirchhoff–Love shell kinematics are employed in this paper to describe our target geometry, that is, soft, thin shells. When a general thin surface is considered in 3D space, its material points are described by a general mapping from the parametric domain [125]. With the curvilinear coordinates ${\xi }^i$ , the general mapping of the current surface’s mid-surface is expressed as ${\mathit{r}=\mathit{x}|}_{{\xi }^3=0}=\mathit{r}({\xi }^{\alpha })$ so that a material point on the surface of a 3D body can be described as $\mathit{x}({\xi }^i)=\mathit{r}+{\xi }^3{\lambda }_3\mathit{n}$ with the thickness-direction coordinate ${\xi }^3:[-h/2,h/2]$ where ${\lambda }_3$ is the thickness-change ratio between the current and initial surface, $\mathit{n}$ is the surface unit normal vector, and $h$ is the initial surface thickness. In this paper, a change in thickness during elastic deformation is allowed for the sake of generality and calculated using the plane-stress condition owing to the thin thickness of shells. We note that most thin-shell theories assume a constant thickness during the elastic deformation using the plane-strain condition and the plane-stress condition. Based on the parameterization of the surface, the covariant tangent vector on the current mid-surface is obtained using the parametric derivative as

$${\mathit{a}}_{\alpha }={\mathit{r}}_{,\alpha }=\frac{\partial \mathit{r}}{\partial {\xi }^{\alpha }}.$$

(5)

Using the covariant tangent vector, the metric tensor (first fundamental form) of the current mid-surface, which has information about the lateral distance between material points on the mid-surface, is obtained as

$$a_{\alpha \beta }={\mathit{a}}_{\alpha }·{\mathit{a}}_{\beta }.$$

(6)

The dual vector is calculated as ${\mathit{a}}^{\alpha }=a^{\alpha \beta }{\mathit{a}}_{\beta }$ with $a^{\alpha \beta }=[a_{\alpha \beta }{]}^{-1}$ such that ${\mathit{a}}^{\alpha }·{\mathit{a}}_{\beta }={\delta }_{\beta }^{\alpha }$ . The unit normal vector to the current mid-surface is obtained as

$$\mathit{n}=\frac{{\mathit{a}}_1\times {\mathit{a}}_2}{||{\mathit{a}}_1\times {\mathit{a}}_2||}$$

(7)

which is called director in the Kirchhoff–Love shell theory. A projection of this surface unit normal vector on the mid-surface gives information about the local curvature of the current mid-surface through the curvature tensor (second fundamental form) as

$$b_{\alpha \beta }=\mathit{n}·{\mathit{a}}_{\alpha ,\beta }=-{\mathit{a}}_{\alpha }·{\mathit{n}}_{,\beta }.$$

(8)

The third fundamental form is also defined as $c_{\alpha \beta }={\mathit{n}}_{,\alpha }·{\mathit{n}}_{,\beta }$ .

Similar to above, all differential geometry calculations on the mid-surface can be also applied for the initial configuration and will be, in this paper, denoted with the over-circle notation, for example, $\stackrel{°}{\mathit{x}}$ , $\stackrel{°}{\mathit{r}}$ , ${\stackrel{°}{\mathit{a}}}_{\alpha }$ , ${\stackrel{°}{\mathit{a}}}^{\alpha }$ , $\stackrel{°}{\mathit{n}}$ , ${\stackrel{°}{a}}_{\alpha \beta }$ , ${\stackrel{°}{a}}^{\alpha \beta }$ , and ${\stackrel{°}{b}}_{\alpha \beta }$ .

From the shell kinematics, the 3D metric tensor for the current and initial surface body can be expressed using the fundamental forms on the mid-surface as

$$g_{\alpha \beta }=a_{\alpha \beta }-2{\xi }^3{\lambda }_3{b}_{\alpha \beta }+O({\xi }^{3^2}),g_{\alpha 3}=0,g_{33}={\lambda }_3^2$$

(9a)

$${\stackrel{°}{g}}_{\alpha \beta }={\stackrel{°}{a}}_{\alpha \beta }-2{\xi }^3{\stackrel{°}{b}}_{\alpha \beta }+O({\xi }^{3^2}),{\stackrel{°}{g}}_{\alpha 3}=0,{\stackrel{°}{g}}_{33}=1,$$

(9b)

where $O({\xi }^{3^2})$ is negligible due to the thin thickness of the shell. Assuming that the intermediate stress-free 3D metric tensor has a similar form for its in-plane components as ${\tilde{g}}_{\alpha \beta }={\tilde{a}}_{\alpha \beta }-2{\xi }^3{\tilde{b}}_{\alpha \beta }+O({\xi }^{3^2})$ yields

$$\{\begin{array}{c}{ϵ}_{\alpha \beta }=E_{\alpha \beta }-{\xi }^3{K}_{\alpha \beta }\hfill \\ \multicolumn{1}{c}{{ϵ}_{\alpha 3}=0}\\ \multicolumn{1}{c}{{ϵ}_{33}=\frac{1}{2}({\lambda }_3^2-{\mathit{F}}_n\stackrel{°}{\mathit{n}}·{\mathit{F}}_n\stackrel{°}{\mathit{n}}),}\end{array}$$

(10)

where $E_{\alpha \beta }=(1/2)(a_{\alpha \beta }-{\tilde{a}}_{\alpha \beta })$ and $K_{\alpha \beta }={\lambda }_3{b}_{\alpha \beta }-{\tilde{b}}_{\alpha \beta }$ are, respectively, the membrane and bending strain on the mid-surface, with the fundamental forms on the mid-surface of the intermediate stress-free state from the differential geometry given by

$${\tilde{a}}_{\alpha \beta }={\left[{\tilde{g}}_{\alpha \beta }\right]}_{{\xi }^3=0}$$

(11a)

$${\tilde{b}}_{\alpha \beta }=-\frac{1}{2}{\left[\frac{\partial {\tilde{g}}_{\alpha \beta }}{\partial {\xi }^3}\right]}_{{\xi }_3=0}.$$

(11b)

3. Material constitutive response accounting for mass change during biological growth

In biological growth problems, an important phenomenon that should be considered in the growth modeling is mass change, as the growth increases the structural size by means of mass addition. That is, the deformation during growth is determined by the mass uptake at growing parts in the body, and thus, mass is a variable that impacts the strain energy density per volume. As discussed in Lee et al. [113], since the internal stress of the body is calculated as the derivative of the strain energy density with respect to strain, it is clear that the mass uptake should impact the internal stress state. In this section, we will elaborate on the concept of how mass addition during growth alters the stress state in a growing body and generalize it to represent the resulting internal stress for all types of stimuli, regardless of whether the mass is added or not.

From thermodynamics, the stress constitutive equation for the second Piola–Kirchhoff stress is derived based on the energy balance and the dissipation inequality (see Appendix 2 for detailed derivation) as

$$\mathbf{\Sigma }={\mathrm{\Sigma }}^{\mathrm{ij}}{\stackrel{°}{\mathit{g}}}_i\otimes {\stackrel{°}{\mathit{g}}}_j\dots {\mathrm{\Sigma }}^{\mathrm{ij}}=2(J\rho )\frac{\partial \psi }{\partial g_{\mathrm{ij}}}$$

(12)

with the basis vector ${\stackrel{°}{\mathit{g}}}_i$ on the initial configuration, in which $\psi$ is the Helmholtz free energy per mass, $g_{\mathrm{ij}}$ is the 3D metric tensor of the current configuration, $J=det\mathit{F}$ is the Jacobian for the deformation gradient $\mathit{F}$ , and $\rho$ is the current mass density per current volume. If there is no mass uptake, $J\rho$ is equal to the initial mass density ${\rho }_0$ from mass conservation, so the second Piola–Kirchhoff stress is calculated as ${\mathrm{\Sigma }}^{\mathrm{ij}}=2{\rho }_0(\partial \psi /\partial g_{\mathrm{ij}})=2(\partial {\mathrm{\Psi }}_0/\partial g_{\mathrm{ij}})$ , where ${\mathrm{\Psi }}_0={\rho }_0\psi$ is the energy density form with no mass addition using the material properties of the initial body. In contrast, if mass uptake occurs, for example, during biological growth, $J\rho \ne {\rho }_0$ . Thus, the interpretation of $J\rho \psi$ has to be established to quantify the effect of mass uptake on the internal energy density and stress.

We consider an infinitesimal element at a material point in a body at which the current mass density $\rho$ can be expressed as $M/V$ , where $M$ and $V$ are the mass and volume of the infinitesimal element on the current configuration. Since the current mass $M$ is $M=M_0+\mathrm{\Delta }{M}_g$ , where $M_0$ is the initial mass and $\mathrm{\Delta }{M}_g$ is the mass increase due to biological growth, $J\rho$ can be calculated as $J\rho ={\rho }_0+\mathrm{\Delta }{M}_g/V_0$ , where $V_0$ is the volume of the corresponding infinitesimal element in the initial configuration. Then, the stress constitutive equation can be rewritten as

$${\mathrm{\Sigma }}^{\mathrm{ij}}=2(J\rho )\frac{\partial \psi }{\partial g_{\mathrm{ij}}}=2\frac{\partial {\mathrm{\Psi }}_0}{\partial g_{\mathrm{ij}}}+2\frac{\mathrm{\Delta }{M}_g}{V_0}\frac{\partial {\psi }_g}{\partial g_{\mathrm{ij}}},$$

(13)

where the second term is the stress modification that arises by accounting for a mass increase during biological growth. We denote the free energy corresponding to the newly added material as ${\psi }_g$ instead of just using $\psi$ because the newly grown material may have different characteristics than the existing material. The Helmholtz free energy that is a quantity “per mass” involves the mass density in the stress constitution. This stress constitution is closely related to material characteristics such as behavior types (e.g. hyperelasticity and viscoelasticity) and material properties (e.g. Young’s modulus and Poisson’s ratio). This implies that the form of free energy is determined by constituents composing the mass and indicates that if the newly grown materials consist of different constituents from the existing material, it requires the use of different free energy forms to describe their different material characteristics. Note that a negative value of $\mathrm{\Delta }{M}_g$ implies shrinkage due to mass decrease. If the added mass is assumed to have the same mechanical behavior and properties as the existing material, the free energy of the new material has the same form as the free energy of the existing material, that is, ${\psi }_g=\psi$ . Then, the effect of mass uptake on the stress can be expressed through the density ratio ${\rho }_r=J\rho /{\rho }_0$ , called stress amplification factor, as

$${\mathrm{\Sigma }}^{\mathrm{ij}}=2{\rho }_r\frac{\partial {\mathrm{\Psi }}_0}{\partial g_{\mathrm{ij}}}.$$

(14)

Thus, the key result of equation (14) is that mass change due to biological growth results in a stress amplification factor of ${\rho }_r$ that modifies the state of stress in the body compared to standard continuum mechanics of $2\partial {\mathrm{\Psi }}_0/\partial g_{\mathrm{ij}}$ that does not account for mass change, assuming that the added material is the same as that of the original body. Note that the density ratio ${\rho }_r$ also physically means the mass ratio between the initial and current mass after growth while we here keep using the density ratio instead of the mass ratio as density is a more standard variable in the mechanics field. We also note that the stress amplification factor $({\rho }_r=J\rho /{\rho }_0)$ derived in this section also intrinsically accounts for the stimulus types without mass addition because $J\rho ={\rho }_0$ for those cases. Hereafter, in this paper, it is assumed that biological growth occurs with mass addition that has the same material behavior and properties as the existing body.

The derived stress amplification in this paper is similar to the work in Ben Amar and Goriely [72] in which the internal stress level is amplified with an extra factor of $J=det\mathit{F}$ which only depends on the area or volume change of the material under consideration. However, as this extra factor is computed using elastic strains from the grown (intermediate) state, this approach gives the same stress amplification regardless of the stimulus type (i.e. heating, swelling, and biological growth), which can be problematic. For a simple example, when a material is swelling by water, the stiffness of the material will be decreased and the material will be softer. On the other hand, if the material undergoes growth by mass addition using the same material as the existing one through density-preserving growth, the stiffness of this material will be the same as before. Nevertheless, the extra factor of Ben Amar and Goriely [72] leads to the same amplification for these very different physical situations. In contrast to this, as the derived stress amplification in this paper is represented by the density ratio $(J\rho /{\rho }_0)$ , it can describe the difference of the stiffness between the stimulus types concerning the mass addition—the case of no mass addition with comparable mechanical properties, such as heating and swelling, gives $J\rho ={\rho }_0$ , that is, ${\rho }_r=1$ . We note that despite the different physical interpretations and mathematical origins of the stress amplification between this paper and Ben Amar and Goriely [72], our stress amplification factor generalizes the work of Ben Amar and Goriely [72] as the extra factor in Ben Amar and Goriely [72] is the same as our stress amplification for a specific situation where (1) the added material during growth is the same as the existing one; (2) the growth is density-preserving; and (3) elastic deformation is incompressible.

This stress modification due to mass addition during biological growth is clearly reconciled with the method of minimum potential energy. The internal potential energy $\mathcal{W}$ is calculated as

$$\mathcal{W}=\int \frac{M_0}{V}\psi +\frac{\mathrm{\Delta }{M}_g}{V}{\psi }_g\mathrm{dV},$$

(15)

where $\mathrm{dV}$ is the volume element of the current body. For the assumption in this paper that the added material during biological growth is the same as the existing material, the potential energy calculation can be expressed using the stress amplification factor as

$$\mathcal{W}=\int \rho \psi \mathrm{dV}=\int J\rho \psi d\stackrel{°}{V}=\int {\rho }_r{\mathrm{\Psi }}_{0}d\stackrel{°}{V},$$

(16)

where $d\stackrel{°}{V}$ is the volume element of the initial body. By virtue of the variation of the internal potential energy with respect to the mechanical part, the method of minimum potential energy yields

$$\begin{array}{c}\multicolumn{1}{c}{\delta \mathcal{W}=[\frac{\partial \mathcal{W}}{\partial g_{\mathrm{ij}}}+\frac{\partial \mathcal{W}}{\partial g_{\mathrm{ji}}}]\frac{1}{2}\delta g_{\mathrm{ij}}}\\ \multicolumn{1}{c}{=\int [\frac{\partial ({\rho }_r{\mathrm{\Psi }}_0)}{\partial g_{\mathrm{ij}}}+\frac{\partial ({\rho }_r{\mathrm{\Psi }}_0)}{\partial g_{\mathrm{ji}}}]\frac{1}{2}\delta g_{\mathrm{ij}}d\stackrel{°}{V}}\\ \multicolumn{1}{c}{\equiv \int {\mathrm{\Sigma }}^{\mathrm{ij}}\frac{1}{2}\delta g_{\mathrm{ij}}d\stackrel{°}{V}}\end{array}$$

(17)

which gives the same stress constitution, relying on the symmetric characteristic of the stress, as equation (14).

Two cases frequently considered in biological growth problems are density-preserving growth and densification (volume-preserving growth). Building on the concept of the multiplicative decomposition of deformation gradient, $\mathit{F}={\mathit{F}}_m{\mathit{F}}_n$ , into the mechanical $({\mathit{F}}_m)$ and non-mechanical $({\mathit{F}}_n)$ parts, the density-preserving growth means $J_m\rho ={\rho }_0$ where $J_m=det{\mathit{F}}_m$ is the Jacobian for the mechanical part of deformation gradient, so the stress amplification factor is obtained as ${\rho }_r=J\rho /{\rho }_0=J_m{J}_n\rho /{\rho }_0=J_n$ where $J_n=det{\mathit{F}}_n$ is the Jacobian for the non-mechanical part of the deformation gradient. Here, $J=J_m{J}_n$ is used from the multiplicative decomposition of deformation gradient. The physical meaning of this stress amplification factor is that ignoring the amplification factor induces a decreased Young’s modulus in the grown body as $1/J_n$ of the initial Young’s modulus. Meanwhile, the densification means $J=1$ , so the stress amplification factor has ${\rho }_r=J\rho /{\rho }_0=\rho /{\rho }_0$ .

In the case of homogeneous biological growth in the body, the stress amplification factor is applied uniformly on the whole body, so there is no effect on the equilibrium configuration. However, if a body grows locally, which is the typical reason for shape-morphing due to biological growth, the mass change may have a significant effect, particularly for processes involving large amounts of growth, such as the embryonic development process as shown in our previous work [113].

3.1. Constitutive response of shell mid-surface

In this paper, the deformation of stimuli-responsive shells is characterized by the resultant stresses on the shell mid-surface, which are used to satisfy the balance of linear momentum. However, those resultant mid-surface stresses are also impacted by the effect of added mass from growth. For general multi-layer surfaces in which each layer may have different characteristics (either mechanical properties or stimulus profile), the stress projection approach can be employed to calculate the resultant mid-surface stresses as follows.

The constitutive equations for the resultant mid-surface membrane $({\tau }_{\sigma }^{\alpha \beta })$ and bending $({\tau }_M^{\alpha \beta })$ stresses are given as [125, 126]

$${\tau }_{\sigma }^{\alpha \beta }=2\frac{\partial {\mathrm{\Psi }}_{2d}}{\partial a_{\alpha \beta }}$$

(18a)

$${\tau }_M^{\alpha \beta }=\frac{\partial {\mathrm{\Psi }}_{2d}}{\partial b_{\alpha \beta }}$$

(18b)

with ${\mathrm{\Psi }}_{2d}=\int J\rho \psi d{\xi }^3$ for the thickness-direction coordinate, ${\xi }^3$ in which $a_{\alpha \beta }$ and $b_{\alpha \beta }$ are, respectively, the metric and curvature tensors on the current, deformed mid-surface, called mid-surface fundamental forms [126]. For thin shells, the Kirchhoff–Love shell kinematics enable us to calculate the in-plane components of the current 3D metric tensor only using the mid-surface metric and curvature tensors as $g_{\alpha \beta }=a_{\alpha \beta }-2{\lambda }_3{\xi }^3{b}_{\alpha \beta }$ , with the thickness-change ratio ${\lambda }_3$ that can be analytically calculated via the plane-stress condition owing to the thin thickness of the shell. Then, using the chain rule, the resultant mid-surface stresses can be calculated for the general multi-layer surface by projecting the second Piola–Kirchhoff stress [127] as

$${\tau }_{\sigma }^{\alpha \beta }=\sum _l\left(\underset{h_b}{\overset{h_t}{\int }}2{\rho }_r\frac{\partial {\mathrm{\Psi }}_0}{\partial g_{\alpha \beta }}d{\xi }^3\right),$$

(19a)

$${\tau }_M^{\alpha \beta }=\sum _l(\underset{h_b}{\overset{h_t}{\int }}-{\lambda }_3{\xi }^{3}2{\rho }_r\frac{\partial {\mathrm{\Psi }}_0}{\partial g_{\alpha \beta }}d{\xi }^3),$$

(19b)

where the subscript $l$ denotes the layer number for considering the general multi-layer surface, and $h_t$ and $h_b$ are the top and bottom thickness coordinates of each layer, respectively. Here, it is assumed that the layers are perfectly bonded. The energy model form for ${\mathrm{\Psi }}_0$ depends on the layer as a surface may consist of layers with different materials. As a specific example, the Saint-Venant energy model, which is a quadratic function of the Green–Lagrange strain ${ϵ}_{\mathrm{ij}}=\frac{1}{2}(g_{\mathrm{ij}}-{\tilde{g}}_{\mathrm{ij}})$ measuring the difference between the current $(g_{\mathrm{ij}})$ and stress-free $({\tilde{g}}_{\mathrm{ij}})$ 3D metric tensors, is given with the plane-stress condition (see Appendix 1 for detailed derivation) by

$${\mathrm{\Psi }}_0=\frac{1}{2}{A}^{\alpha \beta \gamma \delta }{ϵ}_{\alpha \beta }{ϵ}_{\gamma \delta },$$

(20a)

$${\mathrm{\Sigma }}^{\alpha \beta }=2{\rho }_r\frac{\partial {\mathrm{\Psi }}_0}{\partial g_{\alpha \beta }}={\rho }_r{A}^{\alpha \beta \gamma \delta }{ϵ}_{\gamma \delta }$$

(20b)

with the elasticity tensor $A^{\alpha \beta \gamma \delta }={\lambda }_{\mathrm{ps}}{\tilde{g}}^{\alpha \beta }{\tilde{g}}^{\gamma \delta }+\mu ({\tilde{g}}^{\alpha \gamma }{\tilde{g}}^{\beta \delta }+{\tilde{g}}^{\alpha \delta }{\tilde{g}}^{\beta \gamma })$ and ${\tilde{g}}^{\mathrm{ij}}=[{\tilde{g}}_{\mathrm{ij}}{]}^{-1}$ , in which ${\lambda }_{\mathrm{ps}}=\frac{E\nu }{1-{\nu }^2}$ and $\mu =\frac{E}{2(1+\nu )}$ are Lamé constants using Young’s modulus $E$ and Poisson’s ratio $\nu$ of the initial body. Here, the subscript $ps$ of ${\lambda }_{\mathrm{ps}}$ denotes that this Lamé constant is a modified one based on the plane-stress condition from which the thickness-direction strain is given for the Saint-Venant model by ${\tilde{g}}^{33}{ϵ}_{33}=-\left(\frac{\lambda }{\lambda +2\mu }\right){\tilde{g}}^{\alpha \beta }{ϵ}_{\alpha \beta }$ , which results in the modification of the Lamé constant, where the original one is $\lambda =\frac{E\nu }{(1+\nu )(1-2\nu )}$ . Note that the modified Lamé constant differs from this original one for common values of $\nu$ , that is, between 0.3 and 0.5. Then, the constitutive equations for the resultant mid-surface stresses (19) are obtained as

$${\tau }_{\sigma }^{\alpha \beta }=\sum _l\left(\underset{h_b}{\overset{h_t}{\int }}{\rho }_r{A}^{\alpha \beta \gamma \delta }{ϵ}_{\gamma \delta }d{\xi }^3\right),$$

(21a)

$${\tau }_M^{\alpha \beta }=\sum _l(\underset{h_b}{\overset{h_t}{\int }}-{\lambda }_3{\xi }^3{\rho }_r{A}^{\alpha \beta \gamma \delta }{ϵ}_{\gamma \delta }d{\xi }^3).$$

(21b)

4. Equivalent monolayer system

Many thin-shell structures that are actuated by non-mechanical stimuli comprise multiple thin layers, such as the layered aorta growth [128], the pH- or thermo-responsive hydrogel actuators [129, 130], and the bio-inspired nanofiber vascular graft [131]. As differential non-mechanical stimuli induce a change in the rest length, that is, natural stretch, and curvature, that is, natural curvature, of the surface’s mid-surface, we model such multi-layer structures through an equivalent monolayer shell model. Furthermore, we demonstrate how the effects of non-mechanical stimuli, which can be measured experimentally or calculated theoretically, can be accounted for through fundamental forms ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ of the mid-surface in the rest state. To accomplish this, we adopt a conceptual stress-free (rest) mid-surface, which we call the natural mid-surface, to represent the shell energy of general multi-layer surfaces via an equivalent monolayer surface. We define for clarity that overbar symbols, for example, ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ , represent stress-free values for the equivalent monolayer system, whereas tilde symbols, for example, ${\tilde{g}}_{\alpha \beta }$ , represent stress-free values for each layer of the general multi-layer shell system under consideration. The metric and curvature tensors of the natural mid-surface, that is, ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ , give information about rest length and curvature of the mid-surface so that the equivalent shell energy density, adopting the Koiter energy form [56], can be found based on the deviation between the current and natural mid-surface as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Psi }}_{2d}=\int J\rho \psi d{\xi }^3=\int {\rho }_r{\mathrm{\Psi }}_{0}d{\xi }^3}\\ \multicolumn{1}{c}{\stackrel{\mathrm{eq}}{=}\frac{h}{2}{\overline{A}}^{\alpha \beta \gamma \delta }{\overline{E}}_{\alpha \beta }{\overline{E}}_{\gamma \delta }+\frac{h^3}{24}{\overline{A}}^{\alpha \beta \gamma \delta }{\overline{K}}_{\alpha \beta }{\overline{K}}_{\gamma \delta }}\end{array}$$

(22)

such that, in the absence of mechanical external loads, the potential energy $\mathcal{V}$ is calculated as $\mathcal{V}=\int {\mathrm{\Psi }}_{2d}d\stackrel{°}{a}$ with the initial area element $d\stackrel{°}{a}$ , where ${\overline{E}}_{\alpha \beta }=(1/2)(a_{\alpha \beta }-{\overline{a}}_{\alpha \beta })$ and ${\overline{K}}_{\alpha \beta }=b_{\alpha \beta }-{\overline{b}}_{\alpha \beta }$ are the equivalent membrane and bending strains in this equivalent monolayer shell system, respectively, and the elasticity tensor of the natural mid-surface is

$${\overline{A}}^{\alpha \beta \gamma \delta }={\lambda }_{\mathrm{eq}}{\overline{a}}^{\alpha \beta }{\overline{a}}^{\gamma \delta }+{\mu }_{\mathrm{eq}}({\overline{a}}^{\alpha \gamma }{\overline{a}}^{\beta \delta }+{\overline{a}}^{\alpha \delta }{\overline{a}}^{\beta \gamma }),$$

(23)

where ${\lambda }_{\mathrm{eq}}$ and ${\mu }_{\mathrm{eq}}$ are the equivalent Lamé constants. Consequently, the equivalent mid-surface stress resultants are calculated as

$${\tau }_{\sigma }^{\alpha \beta }=2\frac{\partial {\mathrm{\Psi }}_{2d}}{\partial a_{\alpha \beta }}\stackrel{\mathrm{eq}}{=}h{\overline{A}}^{\alpha \beta \gamma \delta }{\overline{E}}_{\gamma \delta },$$

(24a)

$${\tau }_M^{\alpha \beta }=\frac{\partial {\mathrm{\Psi }}_{2d}}{\partial b_{\alpha \beta }}\stackrel{\mathrm{eq}}{=}\frac{h^3}{12}{\overline{A}}^{\alpha \beta \gamma \delta }{\overline{K}}_{\gamma \delta }.$$

(24b)

In general, the stimulus profile may be non-homogeneous over the entire body, which requires that the change in rest length and curvature appears as a second-order tensor, so each component of the natural fundamental forms, that is, ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ , should be, respectively, evaluated according to the stimulus type and profile. To solve for the unknown quantities on the natural mid-surface, that is, fundamental forms ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ , we build on its physical meaning that the stress resultants vanish at $a_{\alpha \beta }={\overline{a}}_{\alpha \beta }$ and $b_{\alpha \beta }={\overline{b}}_{\alpha \beta }$ . While imposing this condition gives a trivial result for equation (24), if used in equation (19), with the Saint-Venant energy model of equation (21), this leads to the following system of equations as

$$\begin{array}{c}\left(\sum _l\underset{h_b}{\overset{h_t}{\int }}\frac{1}{2}{\rho }_r{A}^{\alpha \beta \gamma \delta }d{\xi }^3\right){\overline{a}}_{\gamma \delta }-\left(\sum _l\underset{h_b}{\overset{h_t}{\int }}{\lambda }_3{\xi }^3{\rho }_r{A}^{\alpha \beta \gamma \delta }d{\xi }^3\right){\overline{b}}_{\gamma \delta }\hfill \\ =\left[\sum _l\underset{h_b}{\overset{h_t}{\int }}\frac{1}{2}{\rho }_r{A}^{\alpha \beta \gamma \delta }{\tilde{g}}_{\gamma \delta }d{\xi }^3\right],\hfill \end{array}$$

(25a)

$$\begin{array}{c}\left(\sum _l\underset{h_b}{\overset{h_t}{\int }}\frac{1}{2}{\lambda }_3{\xi }^3{\rho }_r{A}^{\alpha \beta \gamma \delta }d{\xi }^3\right){\overline{a}}_{\gamma \delta }-\left(\underset{h_b}{\overset{h_t}{\int }}{({\lambda }_3{\xi }^3)}^2{\rho }_r{A}^{\alpha \beta \gamma \delta }d{\xi }^3\right){\overline{b}}_{\gamma \delta }\hfill \\ =\left[\sum _l\underset{h_b}{\overset{h_t}{\int }}\frac{1}{2}{\lambda }_3{\xi }^3{\rho }_r{A}^{\alpha \beta \gamma \delta }{\tilde{g}}_{\gamma \delta }d{\xi }^3\right],\hfill \end{array}$$

(25b)

resulting in a system of six equations and six unknowns, from which the fundamental forms of the natural mid-surface (i.e. the three components for each of ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ ) can be obtained relying on that fundamental forms of the mid-surface are symmetric.

Equation (25) is an important result because it enables a direct linkage between the natural fundamental forms and the layer’s stress-free 3D metric tensor $({\tilde{g}}_{\alpha \beta })$ , which can be experimentally measured or theoretically calculated for the non-mechanical stimuli under consideration. Therefore, from the prescribed stress-free 3D metric tensor of each layer for the general multi-layer surface, the fundamental forms ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ of the natural mid-surface, which govern the deformation of the equivalent monolayer shell system, can be obtained using equation (25).

For the special case where the non-mechanical stimulus homogeneously acts on some segment of the shell, the natural mid-surface can be represented by scalar values of $\mathrm{\Lambda }$ (natural stretch) and $\kappa$ (natural curvature) as ${\overline{a}}_{\alpha \beta }={\mathrm{\Lambda }}^2{\stackrel{°}{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }={\mathrm{\Lambda }}^2{\stackrel{°}{b}}_{\alpha \beta }+{\mathrm{\Lambda }}^2\kappa {\stackrel{°}{a}}_{\alpha \beta }$ [48, 99, 113], which enables a simpler system of equations from equation (25) and which may lead to explicit forms of the natural stretch and curvature in some cases. For example, if initially flat bilayer surfaces with thickness $h$ undergo an isotropic area expansion as ${\mathit{F}}_n=e\mathit{I}+(1-e)\stackrel{°}{\mathit{n}}\otimes \stackrel{°}{\mathit{n}}$ with expansion factor $e$ that varies as a function of the thickness but remains homogeneous within each layer, the stress-free 3D metric tensor can be expressed as ${\tilde{g}}_{\alpha \beta }=e_1^2{\stackrel{°}{g}}_{\alpha \beta }$ for the upper layer (i.e. denoted by subscript 1) and ${\tilde{g}}_{\alpha \beta }=e_2^2{\stackrel{°}{g}}_{\alpha \beta }$ for the lower layer (i.e. denoted by subscript 2), with the initial 3D metric tensor ${\stackrel{°}{g}}_{\alpha \beta }={\stackrel{°}{a}}_{\alpha \beta }-2{\xi }^3{\stackrel{°}{b}}_{\alpha \beta }$ using the initial mid-surface fundamental forms ${\stackrel{°}{a}}_{\alpha \beta }$ and ${\stackrel{°}{b}}_{\alpha \beta }$ . Substituting these stress-free metric tensors into the simpler system of equations from equation (25), analytical equations for $\mathrm{\Lambda }$ and $\kappa$ can be derived under the assumption that each layer has similar Poisson’s ratio with an assumption of small thickness change as

$$\begin{array}{c}{\mathrm{\Lambda }}^2=\frac{e_2^6{m}^4{n}^2{\rho }_{r1}^2+e_1^2{e}_2^2\mathrm{mn}(e_1^2(3m(m+1)+1)+e_2^2(m(m+3)+3)){\rho }_{r1}{\rho }_{r2}+e_1^6{\rho }_{r2}^2}{{(e_1{e}_2)}^{-2}[e_2^8{m}^4{n}^2{\rho }_{r1}^2+2e_1^4{e}_2^{4}m(m(2m+3)+2)n{\rho }_{r1}{\rho }_{r2}+e_1^8{\rho }_{r2}^2]},\hfill \end{array}$$

(26)

$$\kappa =-\frac{1}{h}\frac{3e_1^4{e}_2^4(e_1^2-e_2^2)m{(m+1)}^{2}n{\rho }_{r1}{\rho }_{r2}}{(e_2^8{m}^4{n}^2{\rho }_{r1}^2+2e_1^4{e}_2^{4}m(m(2m+3)+2)n{\rho }_{r1}{\rho }_{r2}+e_1^8{\rho }_{r2}^2){\mathrm{\Lambda }}^2},$$

(27)

where subscripts 1 and 2 denote the upper and lower layer, respectively, $m=h_1/h_2$ is the thickness ratio, and $n=E_1/E_2$ is the modulus ratio. A few equations for the scalar natural quantities have been proposed, for example, Pezzulla et al. [99] via the linear geometrical projection and Lucantonio et al. [132] via one-dimensional (1D) beam theory. However, both Pezzulla et al. [99] and Lucantonio et al. [132] did not consider the multidimensional features of shells as well as the mass change effect from biological growth.

As ${\lambda }_{\mathrm{eq}}$ and ${\mu }_{\mathrm{eq}}$ effectively represent the material properties of a general multi-layer surface, one can easily derive those values via its definition for the Saint-Venant energy model with the plane-stress condition as

$$\begin{array}{c}\multicolumn{1}{c}{\int J\rho \psi d{\xi }^3=\int {\rho }_r{\mathrm{\Psi }}_{0}d{\xi }^3}\\ \multicolumn{1}{c}{=\int {\rho }_r\frac{1}{2}{A}^{\alpha \beta \gamma \delta }{ϵ}_{\alpha \beta }{ϵ}_{\gamma \delta }d{\xi }^3}\\ \multicolumn{1}{c}{=\int {\rho }_r\frac{1}{2}[{\lambda }_{\mathrm{ps}}{\tilde{g}}^{\alpha \beta }{\tilde{g}}^{\gamma \delta }+\mu ({\tilde{g}}^{\alpha \gamma }{\tilde{g}}^{\beta \delta }+{\tilde{g}}^{\alpha \delta }{\tilde{g}}^{\beta \gamma })]{ϵ}_{\alpha \beta }{ϵ}_{\gamma \delta }d{\xi }^3}\\ \multicolumn{1}{c}{=\int \frac{1}{2}[{\lambda }_{\mathrm{eq}}{\tilde{g}}^{\alpha \beta }{\tilde{g}}^{\gamma \delta }+{\mu }_{\mathrm{eq}}({\tilde{g}}^{\alpha \gamma }{\tilde{g}}^{\beta \delta }+{\tilde{g}}^{\alpha \delta }{\tilde{g}}^{\beta \gamma })]{ϵ}_{\alpha \beta }{ϵ}_{\gamma \delta }d{\xi }^3,}\end{array}$$

(28)

which shows the equivalent Lamé constants encompass the different material properties and the effect of mass addition in terms of biological growth via the density ratio ${\rho }_r$ . Then, the equivalent Lamé constants can be calculated as

$${\lambda }_{\mathrm{eq}}=\frac{\sum _l[{\int }_{h_b}^{h_t}{\rho }_r{\lambda }_{\mathrm{ps}}{\tilde{g}}^{\alpha \beta }{\tilde{g}}^{\gamma \delta }({\overline{g}}_{\alpha \beta }-{\tilde{g}}_{\alpha \beta })({\overline{g}}_{\gamma \delta }-{\tilde{g}}_{\gamma \delta })d{\xi }^3]}{\sum _l[{\int }_{h_b}^{h_t}{\tilde{g}}^{\alpha \beta }{\tilde{g}}^{\gamma \delta }({\overline{g}}_{\alpha \beta }-{\tilde{g}}_{\alpha \beta })({\overline{g}}_{\gamma \delta }-{\tilde{g}}_{\gamma \delta })d{\xi }^3]},$$

(29a)

$${\mu }_{\mathrm{eq}}=\frac{\sum _l[{\int }_{h_b}^{h_t}{\rho }_r\mu {\tilde{g}}^{\alpha \gamma }{\tilde{g}}^{\beta \delta }({\overline{g}}_{\alpha \beta }-{\tilde{g}}_{\alpha \beta })({\overline{g}}_{\gamma \delta }-{\tilde{g}}_{\gamma \delta })d{\xi }^3]}{\sum _l[{\int }_{h_b}^{h_t}{\tilde{g}}^{\alpha \gamma }{\tilde{g}}^{\beta \delta }({\overline{g}}_{\alpha \beta }-{\tilde{g}}_{\alpha \beta })({\overline{g}}_{\gamma \delta }-{\tilde{g}}_{\gamma \delta })d{\xi }^3]},$$

(29b)

with higher order terms $O(g_{\alpha \beta }-{\overline{g}}_{\alpha \beta })$ neglected, in which ${\overline{g}}_{\alpha \beta }={\overline{a}}_{\alpha \beta }-2{\xi }^3{\overline{b}}_{\alpha \beta }$ as following the Kirchhoff–Love shell kinematics. Note that for the special case where the actual surface is composed of a single layer, that is, ${\overline{g}}_{\alpha \beta }={\tilde{g}}_{\alpha \beta }$ , the equivalent Lamé constants are ${\lambda }_{\mathrm{eq}}={\rho }_r{\lambda }_{\mathrm{ps}}$ and ${\mu }_{\mathrm{eq}}={\rho }_r\mu$ . The mass increase during growth implies that more energy is needed to make a body deform as it grows, which is clearly captured by the appearance of the stress amplification factor in the equivalent Lamé constants.

5. Stimuli-responsive shell theory

The previous derivations demonstrate how a multi-layer surface system can be expressed as an equivalent monolayer system with natural fundamental forms ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ that prescribe how the monolayer shell will deform in response to non-mechanical stimuli. In this section, we will derive the effective strain energy density while accounting for the fundamental forms and will demonstrate how the effects of non-mechanical stimuli can be captured using standard body forces, tractions, and torques that act along the area or perimeter of the shell by treating the initial configuration as a reference stress-free configuration. The resultant weak form then can be solved numerically as a standard mechanical boundary value problem, in which the standard equilibrium equation is solved accounting for the standard internal stresses due to deformations from the initial body and external forces while still accounting for the effects of the non-mechanical stimuli.

We begin by re-expressing the equivalent membrane strain in equation (22) by splitting it as

$${\overline{a}}^{\alpha \gamma }{\overline{E}}_{\gamma \beta }={\overline{a}}^{\alpha \gamma }\frac{1}{2}(a_{\gamma \beta }-{\overline{a}}_{\gamma \beta })={\stackrel{°}{E}}_{\beta }^{\alpha }+L_{\beta }^{\alpha },$$

(30a)

$${\stackrel{°}{E}}_{\beta }^{\alpha }={\stackrel{°}{a}}^{\alpha \gamma }{\stackrel{°}{E}}_{\gamma \beta }={\stackrel{°}{a}}^{\alpha \gamma }\frac{1}{2}(a_{\gamma \beta }-{\stackrel{°}{a}}_{\gamma \beta }),$$

(30b)

$$L_{\beta }^{\alpha }=\frac{1}{2}({\overline{a}}^{\alpha \gamma }-{\stackrel{°}{a}}^{\alpha \gamma })a_{\gamma \beta },$$

(30c)

where equation (30b) is the standard mechanical term, and equation (30c) is the stimulus-induced term. The equivalent bending strain in equation (22) is also re-expressed as

$${\overline{a}}^{\alpha \gamma }{\overline{K}}_{\gamma \beta }={\overline{a}}^{\alpha \gamma }(b_{\gamma \beta }-{\overline{b}}_{\gamma \beta })={\stackrel{°}{K}}_{\beta }^{\alpha }+Z_{\beta }^{\alpha },$$

(31a)

$${\stackrel{°}{K}}_{\beta }^{\alpha }={\stackrel{°}{a}}^{\alpha \gamma }{\stackrel{°}{K}}_{\gamma \beta }={\stackrel{°}{a}}^{\alpha \gamma }(b_{\gamma \beta }-{\stackrel{°}{b}}_{\gamma \beta }),$$

(31b)

$$Z_{\beta }^{\alpha }=({\overline{a}}^{\alpha \gamma }-{\stackrel{°}{a}}^{\alpha \gamma })b_{\gamma \beta }-({\overline{a}}^{\alpha \gamma }{\overline{b}}_{\gamma \beta }-{\stackrel{°}{a}}^{\alpha \gamma }{\stackrel{°}{b}}_{\gamma \beta })$$

(31c)

with its standard mechanical term of equation (31b) and stimulus-induced term of equation (31c). Then, inserting these re-expressed strains into equation (22) yields a new form of the shell energy density which is fully decomposed into an internal strain energy term ${\mathrm{\Psi }}_{2d}^{\mathrm{int}}$ , which depends on the deformation with respect to the initial configuration, and effective external terms $P_{\mathrm{ext}}$ in terms of the effect of non-mechanical stimuli, given as

$${\mathrm{\Psi }}_{2d}={\mathrm{\Psi }}_{2d}^{\mathrm{int}}-(P_{\mathrm{ext}}^s+P_{\mathrm{ext}}^b),$$

(32a)

$${\mathrm{\Psi }}_{2d}^{\mathrm{int}}=\frac{h}{2}{\AA }^{\alpha \beta \gamma \delta }{\stackrel{°}{E}}_{\alpha \beta }{\stackrel{°}{E}}_{\gamma \delta }+\frac{h^3}{24}{\stackrel{°}{A}}^{\alpha \beta \gamma \delta }{\stackrel{°}{K}}_{\alpha \beta }{\stackrel{°}{K}}_{\gamma \delta },$$

(32b)

$$P_{\mathrm{ext}}^s=\frac{h}{2}[-{\lambda }_{\mathrm{eq}}(2{\stackrel{°}{E}}_{\alpha }^{\alpha }{L}_{\gamma }^{\gamma }+L_{\alpha }^{\alpha }{L}_{\gamma }^{\gamma })-2{\mu }_{\mathrm{eq}}(2{\stackrel{°}{E}}_{\beta }^{\alpha }{L}_{\alpha }^{\beta }+L_{\beta }^{\alpha }{L}_{\alpha }^{\beta })],$$

(32c)

$$P_{\mathrm{ext}}^b=\frac{h^3}{24}[-{\lambda }_{\mathrm{eq}}(2{\stackrel{°}{K}}_{\alpha }^{\alpha }{Z}_{\gamma }^{\gamma }+Z_{\alpha }^{\alpha }{Z}_{\gamma }^{\gamma })-2{\mu }_{\mathrm{eq}}(2{\stackrel{°}{K}}_{\beta }^{\alpha }{Z}_{\alpha }^{\beta }+Z_{\beta }^{\alpha }{Z}_{\alpha }^{\beta })]$$

(32d)

in which we note that ${\mathrm{\Psi }}_{2d}^{\mathrm{int}}$ has the same form as the standard Koiter shell energy model, which depends on deformations between the current and initial mid-surface [58], with its elasticity tensor ${\AA }^{\alpha \beta \gamma \delta }={\lambda }_{\mathrm{eq}}{\stackrel{°}{a}}^{\alpha \beta }{\stackrel{°}{a}}^{\gamma \delta }+{\mu }_{\mathrm{eq}}({\stackrel{°}{a}}^{\alpha \gamma }{\stackrel{°}{a}}^{\beta \delta }+{\stackrel{°}{a}}^{\alpha \delta }{\stackrel{°}{a}}^{\beta \gamma })$ , while superscripts $s$ and $b$ on $P_{\mathrm{ext}}$ denote the stretching- and bending-related terms, respectively. Then, we can fully expand the effective external terms via Taylor series about $a_{\alpha \beta }={\overline{a}}_{\alpha \beta }$ and $b_{\alpha \beta }={\overline{b}}_{\alpha \beta }$ as

$$P_{\mathrm{ext}}^s=\frac{h}{2}[{\AA }^{\alpha \beta \gamma \delta }{\mathrm{\Lambda }}_{\alpha \beta }{\mathrm{\Lambda }}_{\gamma \delta }+2{\AA }^{\alpha \beta \gamma \delta }{\overline{E}}_{\alpha \beta }{\mathrm{\Lambda }}_{\gamma \delta }+({\AA }^{\alpha \beta \gamma \delta }-{\overline{A}}^{\alpha \beta \gamma \delta }){\overline{E}}_{\alpha \beta }{\overline{E}}_{\gamma \delta }],$$

(33a)

$$P_{\mathrm{ext}}^b=\frac{h^3}{24}[{\AA }^{\alpha \beta \gamma \delta }{\kappa }_{\alpha \beta }{\kappa }_{\gamma \delta }+2{\AA }^{\alpha \beta \gamma \delta }{\overline{K}}_{\alpha \beta }{\kappa }_{\gamma \delta }+({\AA }^{\alpha \beta \gamma \delta }-{\overline{A}}^{\alpha \beta \gamma \delta }){\overline{K}}_{\alpha \beta }{\overline{K}}_{\gamma \delta }]$$

(33b)

with natural stretch ${\mathrm{\Lambda }}_{\alpha \beta }=\frac{1}{2}({\overline{a}}_{\alpha \beta }-{\stackrel{°}{a}}_{\alpha \beta })$ and natural curvature ${\kappa }_{\alpha \beta }={\overline{b}}_{\alpha \beta }-{\stackrel{°}{b}}_{\alpha \beta }$ as general tensorial forms for local, non-homogeneous stimuli, which measure the change in the rest length and curvature of the mid-surface under the non-mechanical stimuli, respectively. We note that as these natural stretch and curvature have general tensorial forms, the anisotropic characteristics of non-mechanical stimuli can be represented in this theory, such as orthotropic area expansion of soft tissues containing microfibrils induced by swelling [133, 134]. Here, the first terms on the right-hand side (RHS) of equations (33a) and (33b) will not affect the variation of shell energy because they are independent of displacement. Furthermore, in the third terms on the RHS of equations (33a) and (33b), the deviation of elasticity tensors introduces a change in the energy of $O(ϵ)$ , which corresponds to the difference in the metric tensors in the natural state ${\overline{a}}^{\alpha \beta }$ and the undeformed, initial state ${\stackrel{°}{a}}^{\alpha \beta }$ , that is, ${\stackrel{°}{a}}^{\alpha \beta }-{\overline{a}}^{\alpha \beta }=O(ϵ)$ where $ϵ$ is the largest principal extension in the stress-free state from the initial body. In this paper, we will thus neglect the first terms and the third terms associated with the deviation that comes from ${\AA }^{\alpha \beta \gamma \delta }-{\overline{A}}^{\alpha \beta \gamma \delta }$ , as a similar deviation between the fundamental (equilibrium) and initial configuration is ignored in standard continuum mechanics [135], allowing us to finally write the effective external terms of the shell energy density as

$$P_{\mathrm{ext}}^s=h{\AA }^{\alpha \beta \gamma \delta }{\overline{E}}_{\alpha \beta }{\mathrm{\Lambda }}_{\gamma \delta },$$

(34a)

$$P_{\mathrm{ext}}^b=\frac{h^3}{12}{\AA }^{\alpha \beta \gamma \delta }{\overline{K}}_{\alpha \beta }{\kappa }_{\gamma \delta }.$$

(34b)

In this final form of the effective external terms for the effect of non-mechanical stimuli, the variable is the displacement in the current metric and curvature tensor ( $a_{\alpha \beta }$ and $b_{\alpha \beta }$ ), while the natural fundamental forms ( ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ ) are the free input parameters that can be calculated through equation (25) using any available experimental and/or theoretical data that describe how each layer of general multi-layer surface responds to the non-mechanical stimuli. Therefore, we have established a relationship between the effect of non-mechanical stimuli, such as its intensity and direction that are characterized from experimental measurements or theoretical calculation, and displacement, which is similar to the standard work-conjugate relationship between displacement and mechanical loads.

5.1. Variational approach

Having established the novel SRS energy density of equation (32a) with the internal term (32b) and effective external terms (34), we now derive the first and second variations of the potential energy. These are essential for understanding how the non-mechanical stimuli work on the body and enabling a numerical (i.e. finite element-based) solution to boundary value problems (first variation) and also for enabling examinations into the stability of the shells subject to the non-mechanical stimuli (second variation). To do so, we will closely follow the mathematical treatment by Koiter and Van Der Heijden [135], and Thompson and Hunt [136, 137] by starting with equations (32b) and (34) to obtain integral formulations in terms of the equilibrium equation and the stability criteria. For this novel SRS theory, the functional of mechanical potential energy is calculated as

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{V}[\mathit{u}({\xi }^{\alpha })]=\int {\mathrm{\Psi }}_{2d}d\stackrel{°}{a}}\\ \multicolumn{1}{c}{=\int {\mathrm{\Psi }}_{2d}^{\mathrm{int}}[\mathit{u}({\xi }^{\alpha })]d\stackrel{°}{a}-\mathcal{P}[\mathit{u}({\xi }^{\alpha })]}\end{array}$$

(35)

in which the total potential energy is the internal strain energy using equation (32b) minus the external potential using equation (34) in $\mathcal{P}[\mathit{u}({\xi }^{\alpha })]=\int P_{\mathrm{ext}}^s+P_{\mathrm{ext}}^{b}d\stackrel{°}{a}$ , where $\mathit{u}({\xi }^{\alpha })$ is the displacement vector between the current and initial configurations, and $d\stackrel{°}{a}=\sqrt{det{\stackrel{°}{a}}_{\alpha \beta }}d{\xi }^{1}d{\xi }^2$ is the area element of the initial shell with the in-plane coordinates ${\xi }^1$ and ${\xi }^2$ . We will investigate a shell that is in a state of equilibrium by expanding each term in this energy functional via Taylor series about the metric and curvature tensors at the fundamental (equilibrium) state. Then, we find the energy density in the internal strain energy as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Psi }}_{2d}^{\mathrm{int}}={\left[{\mathrm{\Psi }}_{2d}^{\mathrm{int}}\right]}_{*}}\\ \multicolumn{1}{c}{+{\left[2\frac{\partial {\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial a_{\alpha \beta }}\right]}_{*}\frac{1}{2}\delta a_{\alpha \beta }+{\left[\frac{\partial {\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial b_{\alpha \beta }}\right]}_{*}\delta b_{\alpha \beta }}\\ \multicolumn{1}{c}{+\frac{1}{2}{\left[4\frac{{\partial }^2{\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial a_{\gamma \delta }\partial a_{\alpha \beta }}\right]}_{*}\frac{1}{4}\delta a_{\alpha \beta }\delta a_{\gamma \delta }+\frac{1}{2}{\left[\frac{{\partial }^2{\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial b_{\gamma \delta }\partial b_{\alpha \beta }}\right]}_{*}\delta b_{\alpha \beta }\delta b_{\gamma \delta },}\end{array}$$

(36)

where $\delta a_{\alpha \beta }=a_{\alpha \beta }-a_{\alpha \beta }^{*}$ and $\delta b_{\alpha \beta }=b_{\alpha \beta }-b_{\alpha \beta }^{*}$ are the variation in the metric and curvature tensor, respectively, with the superscript and subscript (*) denoting a quantity related to the fundamental state here and hereafter in this section. The second and third terms on the RHS can be further examined by virtue of the variation with the stress constitution as

$$\begin{array}{c}\multicolumn{1}{c}{\delta {\mathrm{\Psi }}_{2d}^{\mathrm{int}}=\frac{1}{2}{[\frac{\partial {\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial a_{\alpha \beta }}+\frac{\partial {\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial a_{\beta \alpha }}]}_{*}\delta a_{\alpha \beta }+\frac{1}{2}{[\frac{\partial {\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial b_{\alpha \beta }}+\frac{\partial {\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial b_{\beta \alpha }}]}_{*}\delta b_{\alpha \beta }}\\ \multicolumn{1}{c}{\equiv {\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }\frac{1}{2}\delta a_{\alpha \beta }+{\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }\delta b_{\alpha \beta },}\end{array}$$

(37)

where the symmetric mid-surface membrane ${\stackrel{°}{\tau }}_{\sigma }^{\alpha \beta }$ and bending ${\stackrel{°}{\tau }}_M^{\alpha \beta }$ stresses, in which the initial configuration is the stress-free configuration, are introduced, that is,

$${\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }={\left[2\frac{\partial {\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial a_{\alpha \beta }}\right]}_{*}={\left[h{\AA }^{\alpha \beta \gamma \delta }{\stackrel{°}{E}}_{\gamma \delta }\right]}_{*},$$

(38a)

$${\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }={\left[\frac{\partial {\mathrm{\Psi }}_{2d}^{\mathrm{int}}}{\partial b_{\alpha \beta }}\right]}_{*}={\left[\frac{h^3}{12}{\AA }^{\alpha \beta \gamma \delta }{\stackrel{°}{K}}_{\gamma \delta }\right]}_{*}.$$

(38b)

Next, we can also expand the integrand of the external potential as

$$P_{\mathrm{ext}}^s+P_{\mathrm{ext}}^b={[P_{\mathrm{ext}}^s+P_{\mathrm{ext}}^b]}_{*}+{\left[\frac{\partial P_{\mathrm{ext}}^s}{\partial a_{\alpha \beta }}\right]}_{*}\delta a_{\alpha \beta }+{\left[\frac{\partial P_{\mathrm{ext}}^b}{\partial b_{\alpha \beta }}\right]}_{*}\delta b_{\alpha \beta }.$$

(39)

Now, we write the functional of mechanical potential energy as

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{V}[\delta \mathit{u}({\xi }^{\alpha })]=\int {\left[{\mathrm{\Psi }}_{2d}^{\mathrm{int}}\right]}_{*}d\stackrel{°}{a}-\int {[P_{\mathrm{ext}}^s+P_{\mathrm{ext}}^b]}_{*}d\stackrel{°}{a}}\\ \multicolumn{1}{c}{+\int {\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }\frac{1}{2}\delta a_{\alpha \beta }+{\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }\delta b_{\alpha \beta }d\stackrel{°}{a}}\\ \multicolumn{1}{c}{+\int \frac{1}{2}\left(h{\AA }^{\alpha \beta \gamma \delta }\right)\frac{1}{4}\delta a_{\alpha \beta }\delta a_{\gamma \delta }+\frac{1}{2}\left(\frac{h^3}{12}{\AA }^{\alpha \beta \gamma \delta }\right)\delta b_{\alpha \beta }\delta b_{\gamma \delta }d\stackrel{°}{a}}\\ \multicolumn{1}{c}{-\int \left(h{\AA }^{\alpha \beta \gamma \delta }{\mathrm{\Lambda }}_{\gamma \delta }\right)\frac{1}{2}\delta a_{\alpha \beta }+\left(\frac{h^3}{12}{\AA }^{\alpha \beta \gamma \delta }{\kappa }_{\gamma \delta }\right)\delta b_{\alpha \beta }d\stackrel{°}{a},}\end{array}$$

(40)

where $\delta \mathit{u}({\xi }^{\alpha })$ is the variation in displacement vector (i.e. a change between the current and fundamental configurations) such that $\delta a_{\alpha \beta }$ and $\delta b_{\alpha \beta }$ are functions of $\delta \mathit{u}$ . For convenience, we introduce $\delta a_{\alpha \beta }=\delta a_{\alpha \beta }^{\mathrm{lin}}+\delta a_{\alpha \beta }^{\mathrm{nl}(2)}+\delta a_{\alpha \beta }^{\mathrm{nl}(h)}$ and $\delta b_{\alpha \beta }=\delta b_{\alpha \beta }^{\mathrm{lin}}+\delta b_{\alpha \beta }^{\mathrm{nl}(2)}+\delta b_{\alpha \beta }^{\mathrm{nl}(h)}$ whose superscripts $\mathrm{lin}$ , $\mathrm{nl}(2)$ , and $\mathrm{nl}(h)$ denote the first-order (linear), second-order, and higher-order terms with respect to the displacement $\delta \mathit{u}$ , respectively.

Since the fundamental state for our investigation is an equilibrium state, the first-order part of the energy functional with respect to the displacement variation $\delta \mathit{u}$ , called first variation of energy functional, must be stationary for all admissible displacement fields. As the first variation only contains the linear terms in $\delta \mathit{u}$ of equation (40), this leads to

$${\delta }^1\mathcal{V}=\int {\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }\frac{1}{2}\delta a_{\alpha \beta }^{\mathrm{lin}}+{\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }\delta b_{\alpha \beta }^{\mathrm{lin}}d\stackrel{°}{a}-\int f_{\mathrm{\Lambda }}^{\alpha \beta }\frac{1}{2}\delta a_{\alpha \beta }^{\mathrm{lin}}+f_{\kappa }^{\alpha \beta }\delta b_{\alpha \beta }^{\mathrm{lin}}d\stackrel{°}{a}=0,$$

(41)

where ${\delta }^1\mathcal{V}$ is the first variation of energy functional with superscript 1 denoting the linear-order part with respect to $\delta \mathit{u}$ , with

$$f_{\mathrm{\Lambda }}^{\alpha \beta }=h{\AA }^{\alpha \beta \gamma \delta }{\mathrm{\Lambda }}_{\gamma \delta },$$

(42a)

$$f_{\kappa }^{\alpha \beta }=\frac{h^3}{12}{\AA }^{\alpha \beta \gamma \delta }{\kappa }_{\gamma \delta }.$$

(42b)

Equation (41) represents the balance of linear momentum in the integral form.

Now, we will further manipulate equation (41) to explicitly determine how the non-mechanical stimuli act to deform the shell. If we denote the displacement field $\delta \mathit{u}=\delta u_i{\stackrel{°}{\mathit{a}}}^i=\delta u_{\alpha }{\stackrel{°}{\mathit{a}}}^{\alpha }+\delta w\stackrel{°}{\mathit{n}}$ , that is, the in-plane displacement $\delta u_{\alpha }$ and the normal displacement $\delta u_3=\delta w$ , we can approximate [126]

$$\frac{1}{2}\delta a_{\alpha \beta }^{lin}=\frac{1}{2}({\stackrel{°}{\nabla }}_{\alpha }\delta u_{\beta }+{\stackrel{°}{\nabla }}_{\beta }\delta u_{\alpha })-{\stackrel{°}{b}}_{\alpha \beta }\delta w,$$

(43a)

$$\delta b_{\alpha \beta }^{lin}={\stackrel{°}{\nabla }}_{\alpha }{\stackrel{°}{\nabla }}_{\beta }\delta w+{\stackrel{°}{b}}_{\alpha \gamma }{\stackrel{°}{\nabla }}_{\beta }\delta u^{\gamma }+{\stackrel{°}{b}}_{\beta \gamma }{\stackrel{°}{\nabla }}_{\alpha }\delta u^{\gamma }+\delta u^{\gamma }{\stackrel{°}{\nabla }}_{\beta }{\stackrel{°}{b}}_{\gamma \alpha }-{\stackrel{°}{b}}_{\beta \gamma }{\stackrel{°}{b}}_{\alpha }^{\gamma }\delta w,$$

(43b)

where ${\stackrel{°}{\nabla }}_{\alpha }$ is the covariant derivative with respect to the initial configuration. Here and hereafter in this section, we use the initial metric tensor to raise and lower indices. Then, inserting equation (43) only into the external potential term, that is, the second integral term of equation (41), and using the Stokes theorem give

$$\begin{array}{c}\multicolumn{1}{c}{{\delta }^1\mathcal{V}=\int {\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }\frac{1}{2}\delta a_{\alpha \beta }^{\mathrm{lin}}+{\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }\delta b_{\alpha \beta }^{\mathrm{lin}}d\stackrel{°}{a}}\\ \multicolumn{1}{c}{-\int \underset{\genfrac{}{}{0ex}{}{\mathrm{In}-\mathrm{plane}\mathrm{body}\mathrm{force}}{\equiv {\mathcal{B}}_{\mathrm{in}}^{\beta }}}{\underbrace{[-{\stackrel{°}{\nabla }}_{\alpha }{f}_{\mathrm{\Lambda }}^{\alpha \beta }-2{\stackrel{°}{\nabla }}_{\alpha }{f}_{\kappa }^{\alpha \theta }{\stackrel{°}{b}}_{\theta }^{\beta }-f_{\kappa }^{\alpha \theta }{\stackrel{°}{\nabla }}_{\theta }{\stackrel{°}{b}}_{\alpha }^{\beta }]}}\delta u_{\beta }d\stackrel{°}{a}}\\ \multicolumn{1}{c}{-\int \underset{\genfrac{}{}{0ex}{}{\mathrm{Out}-\mathrm{of}-\mathrm{plane}\mathrm{body}\mathrm{force}(\mathrm{pressure})}{\equiv {\mathcal{B}}_{\mathrm{out}}}}{\underbrace{[-f_{\mathrm{\Lambda }}^{\alpha \beta }{\stackrel{°}{b}}_{\alpha \beta }+{\stackrel{°}{\nabla }}_{\beta }{\stackrel{°}{\nabla }}_{\alpha }{f}_{\kappa }^{\alpha \beta }-f_{\kappa }^{\alpha \beta }{\stackrel{°}{b}}_{\beta \gamma }{\stackrel{°}{b}}_{\alpha }^{\gamma }]}}\delta wd\stackrel{°}{a}}\\ \multicolumn{1}{c}{-\oint \underset{\genfrac{}{}{0ex}{}{\mathrm{In}-\mathrm{plane}\mathrm{traction}}{\equiv {\mathcal{T}}_{\mathrm{in}}^{\beta }}}{\underbrace{[(f_{\mathrm{\Lambda }}^{\alpha \beta }+f_{\kappa }^{\alpha \theta }{\stackrel{°}{b}}_{\theta }^{\beta }){\stackrel{°}{m}}_{\alpha }]}}\delta u_{\beta }+\underset{\genfrac{}{}{0ex}{}{\mathrm{Out}-\mathrm{of}-\mathrm{plane}\mathrm{traction}}{\equiv {\mathcal{T}}_{\mathrm{out}}}}{\underbrace{[-{\stackrel{°}{\nabla }}_{\alpha }{f}_{\kappa }^{\alpha \beta }{\stackrel{°}{m}}_{\beta }]}}\delta w+\underset{\genfrac{}{}{0ex}{}{\mathrm{Torque}}{\equiv t^{\beta }}}{\underbrace{\left[f_{\kappa }^{\alpha \beta }{\stackrel{°}{m}}_{\alpha }\right]}}\delta q_{\beta }d\stackrel{°}{s}}\\ \multicolumn{1}{c}{=0,}\end{array}$$

(44)

where ${\stackrel{°}{m}}_{\alpha }=\stackrel{°}{\mathit{m}}·{\stackrel{°}{\mathit{a}}}_{\alpha }$ is the unit normal vector to the initial boundary curve, and $d\stackrel{°}{s}$ is the initial boundary curve element. Here, the rotation vector $\delta q_{\alpha }={\stackrel{°}{\nabla }}_{\alpha }\delta w+{\stackrel{°}{b}}_{\alpha }^{\gamma }\delta u_{\gamma }$ is used such that $\delta \mathit{n}=-\delta q_{\alpha }{\stackrel{°}{\mathit{a}}}^{\alpha }$ [126].

Equation (44) shows how the non-mechanical stimuli, which change the rest length and curvature of the mid-surface via ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ , act to deform the shell. Specifically, non-mechanical stimuli do not impact the internal potential energy (the first line of the RHS in equation (44)). Instead, the non-mechanical stimuli act as effective external loads to deform the initial, undeformed body as the external terms are functions of the shell displacements and their derivatives (the second, third, and fourth lines of the RHS in equation (44)). Then, based on the well-known mechanical work conjugates, we can see that the effect of non-mechanical stimuli plays a similar role as mechanical, external loadings do in classical mechanics as shown in Figure 2—for body forces, in-plane body force ${\mathcal{B}}_{\mathrm{in}}={\mathcal{B}}_{\mathrm{in}}^{\alpha }{\stackrel{°}{\mathit{a}}}_{\alpha }$ (conjugated with $\delta u_{\alpha }$ ) and out-of-plane body force ${\mathcal{B}}_{\mathrm{out}}={\mathcal{B}}_{\mathrm{out}}\stackrel{°}{\mathit{n}}$ (conjugated with $\delta w$ ), that is, pressure; for boundary loads, in-plane traction ${\mathit{T}}_{\mathrm{in}}=T_{\mathrm{in}}^{\alpha }{\stackrel{°}{\mathit{a}}}_{\alpha }$ (conjugated with $\delta u_{\alpha }$ ), out-of-plane traction ${\mathit{T}}_{\mathrm{out}}=T_{\mathrm{out}}\stackrel{°}{\mathit{n}}$ (conjugated with $\delta w$ ), and torque $\mathit{t}=t^{\alpha }{\stackrel{°}{a}}_{\alpha }$ (conjugated with $\delta q_{\alpha }$ ).

Figure 2.

Transformation of the effect of non-mechanical stimulus into effective external loads according to the mechanical work conjugates. Blue arrows denote that those are body force-type effective loads, while red arrows mean boundary forces and torque.

Furthermore, the weak form of equation (44) can be utilized as the starting point to use well-established numerical techniques, such as the finite element method, to solve mechanical boundary value problems where all external loads, due to either mechanical loads or non-mechanical stimuli, are applied as standard mechanical, external loadings. We do note that the finite element approaches for the shell analysis must satisfy at least $C^1$ continuity of the shape functions due to the presence of the second derivative in the curvature tensor. However, this offers the distinct advantage that mechanicians and non-mechanicians alike can easily modify existing computational tools with additional external forces to study how various types of non-mechanical stimuli impact the mechanics and physics of thin-shell structures.

In addition, inserting equation (43) into both integral terms in equation (41) and using the Stokes theorem give the Euler–Lagrange equations as

$${\stackrel{°}{\nabla }}_{\alpha }({\stackrel{°}{\tau }}_{{\sigma }_{\ast }}^{\alpha \beta }-f_{\mathrm{\Lambda }}^{\alpha \beta })+2{\stackrel{°}{\nabla }}_{\alpha }({\stackrel{°}{\tau }}_{M_{\ast }}^{\alpha \theta }-f_{\kappa }^{\alpha \theta }){\underset{\theta }{\overset{°}{b}}}^{\beta }+({\stackrel{°}{\tau }}_{M_{\ast }}^{\alpha \theta }-f_{\kappa }^{\alpha \theta }){\stackrel{°}{\nabla }}_{\theta }{\underset{\alpha }{\overset{°}{b}}}^{\beta }=0,$$

(45a)

$$({\stackrel{°}{\tau }}_{{\sigma }_{\ast }}^{\alpha \beta }-f_{\mathrm{\Lambda }}^{\alpha \beta })\underset{\alpha \beta }{\overset{°}{b}}-{\stackrel{°}{\nabla }}_{\beta }{\stackrel{°}{\nabla }}_{\alpha }({\stackrel{°}{\tau }}_{M_{\ast }}^{\alpha \beta }-f_{\kappa }^{\alpha \beta })+({\stackrel{°}{\tau }}_{M_{\ast }}^{\alpha \beta }-f_{\kappa }^{\alpha \beta })\underset{\beta \gamma }{\overset{°}{b}}{\underset{\alpha }{\overset{°}{b}}}^{\gamma }=0$$

(45b)

with boundary conditions as

$$[({\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }-f_{\mathrm{\Lambda }}^{\alpha \beta })+2({\stackrel{°}{\tau }}_{M_{*}}^{\alpha \theta }-f_{\kappa }^{\alpha \theta }){\stackrel{°}{b}}_{\theta }^{\beta }]{\stackrel{°}{m}}_{\alpha }^{°}=0,$$

(46a)

$$({\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }-f_{\kappa }^{\alpha \beta }){\stackrel{°}{m}}_{\alpha }=0,$$

(46b)

$${\stackrel{°}{\nabla }}_{\alpha }({\stackrel{°}{\tau }}_{M_{\ast }}^{\alpha \beta }-f_{\kappa }^{\alpha \beta }){\stackrel{°}{m}}_{\beta }=0.$$

(46c)

These equations and boundary conditions of equations (45) and (46) form a boundary value problem. Note that equations (45a), (46a), and (46b) have a free index of $\beta =1,2$ as vector forms, while equations (45b) and (46c) are scalar forms.

Then, accounting for equation (41), the functional of mechanical potential energy becomes

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{V}[\delta \mathit{u}({\xi }^{\alpha })]=\int {\left[{\mathrm{\Psi }}_{2d}^{\mathrm{int}}\right]}_{*}d\stackrel{°}{a}-\int {[P_{\mathrm{ext}}^s+P_{\mathrm{ext}}^b]}_{*}d\stackrel{°}{a}}\\ \multicolumn{1}{c}{+\int {\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }\frac{1}{2}(\delta a_{\alpha \beta }^{\mathrm{nl}(2)}+\delta a_{\alpha \beta }^{\mathrm{nl}(h)})+{\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }(\delta b_{\alpha \beta }^{\mathrm{nl}(2)}+\delta b_{\alpha \beta }^{\mathrm{nl}(h)})d\stackrel{°}{a}}\\ \multicolumn{1}{c}{+\int \frac{1}{2}\left(h{\AA }^{\alpha \beta \gamma \delta }\right)\frac{1}{4}\delta a_{\alpha \beta }\delta a_{\gamma \delta }+\frac{1}{2}\left(\frac{h^3}{12}{\AA }^{\alpha \beta \gamma \delta }\right)\delta b_{\alpha \beta }\delta b_{\gamma \delta }d\stackrel{°}{a}}\\ \multicolumn{1}{c}{-\int f_{\mathrm{\Lambda }}^{\alpha \beta }\frac{1}{2}(\delta a_{\alpha \beta }^{\mathrm{nl}(2)}+\delta a_{\alpha \beta }^{\mathrm{nl}(h)})+f_{\kappa }^{\alpha \beta }(\delta b_{\alpha \beta }^{\mathrm{nl}(2)}+\delta b_{\alpha \beta }^{\mathrm{nl}(h)})d\stackrel{°}{a}.}\end{array}$$

(47)

The equilibrium state is stable if the functional of mechanical potential energy is positive definite, which can be tested via the character of the second-order part in the energy functional with respect to the displacement variation $\delta \mathit{u}$ , called second variation of energy functional, that is, ${\delta }^2\mathcal{V}\ge 0$ with ${\delta }^2\mathcal{V}$ denoting the second variation [135–137]. By virtue of variational calculus, the second variation only contains the second-order terms in $\delta \mathit{u}$ , which correspond to terms containing $\delta a_{\alpha \beta }^{\mathrm{nl}(2)}$ and $\delta b_{\alpha \beta }^{\mathrm{nl}(2)}$ as well as the linear terms in the product $\delta a_{\alpha \beta }\delta a_{\gamma \delta }$ and $\delta b_{\alpha \beta }\delta b_{\gamma \delta }$ , that is, $\delta a_{\alpha \beta }^{\mathrm{lin}}\delta a_{\gamma \delta }^{\mathrm{lin}}$ and $\delta b_{\alpha \beta }^{\mathrm{lin}}\delta b_{\gamma \delta }^{\mathrm{lin}}$ , as

$$\begin{array}{c}\multicolumn{1}{c}{{\delta }^2\mathcal{V}=\int {\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }\frac{1}{2}\delta a_{\alpha \beta }^{\mathrm{nl}(2)}+{\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }\delta b_{\alpha \beta }^{\mathrm{nl}(2)}d\stackrel{°}{a}-\int f_{\mathrm{\Lambda }}^{\alpha \beta }\frac{1}{2}\delta a_{\alpha \beta }^{\mathrm{nl}(2)}+f_{\kappa }^{\alpha \beta }\delta b_{\alpha \beta }^{\mathrm{nl}(2)}d\stackrel{°}{a}}\\ \multicolumn{1}{c}{+\int \frac{1}{2}\left(h{\AA }^{\alpha \beta \gamma \delta }\right)\frac{1}{4}\delta a_{\alpha \beta }^{\mathrm{lin}}\delta a_{\gamma \delta }^{\mathrm{lin}}+\frac{1}{2}\left(\frac{h^3}{12}{\AA }^{\alpha \beta \gamma \delta }\right)\delta b_{\alpha \beta }^{\mathrm{lin}}\delta b_{\gamma \delta }^{\mathrm{lin}}d\stackrel{°}{a}.}\end{array}$$

(48)

Therefore, the stability of the stimuli-responsive shells under investigation can be determined by the character of equation (48).

6. Validation

In this section, the performance of the proposed SRS theory is illustrated by several benchmark problems under plane-strain conditions, considering different types of stimuli and their corresponding impact on the natural stretch and curvature of shells. As the natural fundamental forms ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ are prescribed, the shell deformations are obtained analytically as well as numerically using the derived weak form (44) with the isogeometric analysis (IGA), which is a modern finite element method well-suited for shell problems due to its shape function’s higher continuity [138], and verified by available reference solutions. Details of the IGA formulation, including discretization of the weak form, can be found in Appendix 3.

6.1. Plates

During differential swelling [101] and growth [113] of thin shells, these non-mechanical stimuli can induce isotropic area expansion that is homogeneous in the in-plane direction but varies in the thickness direction. This results in a change in not only the rest length (natural stretch) but also the rest curvature (natural stretch) in the mid-surface. These changes can be represented by scalar values of $\mathrm{\Lambda }$ and $\kappa$ via ${\overline{a}}_{\alpha \beta }={\mathrm{\Lambda }}^2{\stackrel{°}{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }={\mathrm{\Lambda }}^2{\stackrel{°}{b}}_{\alpha \beta }+{\mathrm{\Lambda }}^2\kappa {\stackrel{°}{a}}_{\alpha \beta }$ as stated in Section 4. If considering the action of the scalar natural stretch and curvature on plates with thickness $h$ , it has been shown that the resulting deformation is bending-dominant and that the plates bend into cylindrical shapes with Germain (mean) curvature $H$ , which can be found analytically for incompressible plates as [48, 99]

$$\mathrm{Hh}=\frac{3}{4}\kappa h,$$

(49)

which indicates that the mean curvature change is solely due to the natural curvature. Thus, we shall test our proposed SRS theory under pure natural curvature based on this plate problem.

By first recalling that the shell deformation is dominated by bending because the shells are thin, we can analytically obtain the mean curvature by minimizing the bending-related parts in the mechanical potential energy of equation (35) by assuming that the stretching-related terms vanish as $a_{\alpha \beta }={\overline{a}}_{\alpha \beta }$ , meaning the shell deformation is entirely due to bending. Then, since the Gauss–Codazzi–Mainardi equations for the surface constraint admit homogeneous second fundamental forms for a flat homogeneous metric, the solution can be calculated by minimizing the bending-related parts with a constraint on the null Gaussian curvature via a Lagrange multiplier. As aforementioned, the homogeneous natural quantities allow us to use ${\overline{a}}_{\alpha \beta }={\mathrm{\Lambda }}^2{\stackrel{°}{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }={\mathrm{\Lambda }}^2{\stackrel{°}{b}}_{\alpha \beta }+{\mathrm{\Lambda }}^2\kappa {\stackrel{°}{a}}_{\alpha \beta }$ , and inserting these relations with ${\stackrel{°}{b}}_{\alpha \beta }=0$ for plates into equation (35) gives its mechanical potential energy as

$$\begin{array}{c}\mathcal{V}=\int \frac{1}{2}{B}_{\mathrm{eq}}({b_1^1}^2+2b_1^1{b}_2^2\nu +2{b_2^1}^2(1-\nu )+{b_2^2}^2)\hfill \\ -B_{\mathrm{eq}}(1+\nu )\kappa {\mathrm{\Lambda }}^2(b_1^1+b_2^2-2\kappa {\mathrm{\Lambda }}^2)+q_m({b_2^1}^2-b_1^1{b}_2^2)d\stackrel{°}{a}\hfill \end{array}$$

(50)

with $b_{\beta }^{\alpha }={\stackrel{°}{a}}^{\alpha \gamma }{b}_{\gamma \beta }\ne 0$ and the bending rigidity $B_{\mathrm{eq}}$ , in which $q_m$ is the Lagrange multiplier to impose the constraint of the null Gaussian curvature. Here, we assume that Poisson’s ratio $\nu$ is homogeneous in the body. If we denote the integrand of this equation as $\mathcal{L}$ , since all variables in the integrand are homogeneous, the solution of minimum potential energy can be simply obtained via $\partial \mathcal{L}/\partial b_1^1=\partial \mathcal{L}/\partial b_2^2=\partial \mathcal{L}/\partial b_2^1=\partial \mathcal{L}/\partial q_m=0$ , which gives

$$\{\begin{array}{c}{b}_1^1=(1+\nu )\kappa {\mathrm{\Lambda }}^2,b_2^2=0,b_2^1=0\mathrm{with}{q}_m=-B_{\mathrm{eq}}(1-\nu )\hfill \\ \multicolumn{1}{c}{b_1^1=0,b_2^2=(1+\nu )\kappa {\mathrm{\Lambda }}^2,b_2^1=0\mathrm{with}{q}_m=-B_{\mathrm{eq}}(1-\nu )}\end{array}$$

(51)

in which both solutions have dimensionless mean curvature $\mathrm{Hh}=\frac{1}{2}{a}^{\alpha \beta }{b}_{\alpha \beta }h=\frac{3}{4}\kappa h$ for incompressible materials, which is consistent with the previous work of equation (49), while also matching our simulation results with $E=1$ MPa and $\nu =0.5$ as shown in Figure 3.

Figure 3.

Change in mean curvature of an incompressible plate as a function of natural curvature. Symbols denote IGA simulation results obtained via the proposed stimuli-responsive shell theory based on equation (44), while the solid line is the analytical solution of equation (49) for which our proposed theory gives the same solution as the previous works [48, 99].

6.2. Cylinders

The previous example for flat plates considered a scalar natural stretch $\mathrm{\Lambda }$ and curvature $\kappa$ that are induced by a stimulus like differential swelling in the thickness direction [99, 139], which results in area expansion that is isotropic in-plane but varies in the thickness direction. In contrast, cylinders have an initial curvature, and as a result, a different deformation feature is observed for the same stimulus profile. Cylinders unroll and snap like a tape-spring, or “snap bracelet” [140], when exposed to curvature stimulus whose sign is the opposite of its initial curvature. This bending-dominant deformation of the cylinder with radius $R$ and thickness $h$ has been well studied in Pezzulla et al. [99] where the deformed mean curvature is analytically solved as

$$\mathrm{Hh}=\frac{\mathrm{hf}(\kappa ;\mathrm{\Lambda })}{2R}+\frac{(1+\nu )\kappa h}{2},f(\kappa ;\mathrm{\Lambda })=\{\begin{array}{c}\frac{1}{\mathrm{\Lambda }}\mathrm{for}\kappa >{\kappa }_s\hfill \\ \multicolumn{1}{c}{\frac{\nu }{\mathrm{\Lambda }}\mathrm{for}\kappa <{\kappa }_s}\end{array}$$

(52)

in which ${\kappa }_s=-\frac{1}{2\mathrm{\Lambda }R}$ is the critical natural curvature when the cylinder curvature changes sign. In this case, both the natural stretch and curvature appear to affect the mean curvature change during the deformation. However, it has been shown that the thickness-dependent stimulus results in minimal stretching of the shell’s mid-surface [51], so the effect of natural stretching is weak during the cylinder deformation, that is, $\mathrm{\Lambda }\simeq 1$ . Thus, we shall test our proposed novel SRS theory under a combination of natural curvature and weak stretching based on this cylinder problem. This problem also demonstrates the applicability of our proposed theory for the initially curved geometry and its ability to capture the deformation preference for a bi-stable problem [99].

Our numerical solutions via the proposed SRS theory accurately capture the evolution of the dimensionless mean curvature as well as the critical point of the sign change of curvature as shown in Figure 4. In addition, the critical point is analytically evaluated through the proposed shell theory. As in the above plate problem, for thin cylindrical shells, the shell’s bending-dominated deformation in which the stretching-related parts of the potential energy vanish, that is, $a_{\alpha \beta }={\overline{a}}_{\alpha \beta }$ , leads to homogeneous second fundamental forms for the zero Gaussian curvature of initial cylinder shapes. Then, we have for initial cylindrical shapes with a cylindrical coordinate system

$${\stackrel{°}{b}}_{\beta }^{\alpha }={\stackrel{°}{a}}^{\alpha \gamma }{\stackrel{°}{b}}_{\gamma \beta }=\left(\begin{array}{cc}1/R\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right)$$

(53)

which gives the mechanical potential energy of equation (35) using ${\overline{a}}_{\alpha \beta }={\mathrm{\Lambda }}^2{\stackrel{°}{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }={\mathrm{\Lambda }}^2{b}_{\alpha \beta }+{\mathrm{\Lambda }}^2\kappa {\stackrel{°}{a}}_{\alpha \beta }$ by assuming a homogeneous Poisson’s ratio as

$$\begin{array}{c}\mathcal{V}=\int \frac{B_{\mathrm{eq}}({\mathcal{X}}_{c1}+{\mathcal{X}}_{c2}+{\mathcal{X}}_{c3})}{2R^2}+q_m({b_2^1}^2-b_1^1{b}_2^2)d\stackrel{°}{a}\hfill \\ \dots \{\begin{array}{c}{\mathcal{X}}_{c1}=R^2[b_1^{1^2}+2b_1^1{b}_2^2\nu +2b_2^{1^2}(1-\nu )+b_2^{2^2}]\hfill \\ \multicolumn{1}{c}{{\mathcal{X}}_{c2}=-2{\mathrm{\Lambda }}^2+2{\mathrm{\Lambda }}^4[2\kappa R(1+\nu )(\kappa R+1)+1]+1}\\ \multicolumn{1}{c}{{\mathcal{X}}_{c3}=-2{\mathrm{\Lambda }}^{2}R[\kappa (1+\nu )\{R(b_1^1+b_2^2)+1\}+b_1^1+b_2^2\nu ]}\end{array}\hfill \end{array}$$

(54)

with $b_{\beta }^{\alpha }={\stackrel{°}{a}}^{\alpha \gamma }{b}_{\gamma \beta }\ne 0$ . Since the integrand $\mathcal{L}$ of this equation only includes homogeneous variables like the above plate problem, the solution of minimum potential energy is obtained via $\partial \mathcal{L}/\partial b_1^1=\partial \mathcal{L}/\partial b_2^2=\partial \mathcal{L}/\partial b_2^1=\partial \mathcal{L}/\partial q_m=0$ as

$$\{\begin{array}{c}{b}_1^1={\mathrm{\Lambda }}^2[\frac{1}{R}+(1+\nu )\kappa ],b_2^2=0\mathrm{with}{q}_m=-\frac{B_{\mathrm{eq}}R(1-{\nu }^2)\kappa }{1+R(1+\nu )\kappa }\hfill \\ \multicolumn{1}{c}{b_1^1=0,b_2^2={\mathrm{\Lambda }}^2[\frac{\nu }{R}+(1+\nu )\kappa ]\mathrm{with}{q}_m=-\frac{B_{\mathrm{eq}}(1-{\nu }^2)(1+R\kappa )}{\nu +R(1+\nu )\kappa }.}\end{array}$$

(55)

Figure 4.

Change in mean curvature of an incompressible cylinder due to natural stretch and curvature, with the critical transition point of equation (56), which is analytically obtained via the proposed SRS theory, denoted using the black vertical dashed line. Symbols denote IGA simulation results obtained via the proposed stimuli-responsive shell theory. The black solid line is the reference analytical solution of equation (52) [99], while the colored dashed line is the analytical solution of equation (57) using the proposed stimuli-responsive shell theory.

As both these solutions provide the stationary states of the cylinder under scalar natural quantities, the cylinder deforms in a way to have smaller potential energy among the two stationary solutions, that is, bi-stable state—one keeps the original bending orientation, while the other has a rotated configuration. Hence, inserting these two solutions into equation (54) gives two energies corresponding to a deformed shape with the original bending direction and the rotated one after snapping, and equating those two energies gives the transition point at which the cylinder switches its bending orientation with the critical natural curvature as

$${\kappa }_s=-\frac{1}{2R}.$$

(56)

Then, the dimensionless mean curvature is analytically obtained as

$$\mathrm{Hh}=\frac{\mathrm{hf}(\kappa )}{2R}+\frac{(1+\nu )\kappa h}{2},f(\kappa )=\{\begin{array}{c}1\mathrm{for}\kappa >{\kappa }_s\hfill \\ \multicolumn{1}{c}{\nu \mathrm{for}\kappa <{\kappa }_s}\end{array}.$$

(57)

Here, an error of $O(1-\mathrm{\Lambda })$ exists between the previous work of equation (52) and the analytical solutions of equation (57) obtained via the proposed SRS theory. This error comes from the fact that the solution obtained by the proposed theory stems from the natural curvature tensor written as ${\overline{b}}_{\alpha \beta }={\mathrm{\Lambda }}^2{\stackrel{°}{b}}_{\alpha \beta }+{\mathrm{\Lambda }}^2\kappa {\stackrel{°}{a}}_{\alpha \beta }$ , while the previous work of equation (52) utilized ${\overline{b}}_{\alpha \beta }=\mathrm{\Lambda }{\stackrel{°}{b}}_{\alpha \beta }+{\mathrm{\Lambda }}^2\kappa {\stackrel{°}{a}}_{\alpha \beta }$ for their solution. For the isotropic area expansion induced by non-mechanical stimuli, the correct expression is the one used in our solution. However, the natural curvature for thin bodies undergoing differential stimulus through the thickness direction will barely or not stretch the shell’s mid-surface [51], so this error is negligible. Therefore, the analytical solutions via our proposed novel SRS theory have good agreement with the previous work as shown in Figure 4.

6.3. Annulus

The purpose of this verification example is to both examine the performance of the proposed SRS theory for a stretching-dominated problem and demonstrate the ability to examine the onset of buckling instability through a linear stability analysis. When a monolayer plate annulus with radius $R\in [R_{\mathrm{in}},R_{\mathrm{out}}]$ and thickness $h$ is subject to an inhomogeneous stimulus profile such that the natural metric tensor has a spherical form with the radial distance $r$ from the center as

$$\left[{\overline{a}}_{\alpha \beta }\right]=\left[\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {(\sin r)}^2\hfill \end{array}\right]$$

(58)

along with a flat natural curvature tensor, the circular plate annulus remains flat for thicknesses above a critical value for a buckling transition ( $h_b$ ), while the plate buckled into a spherical shape for thicknesses below the critical buckling transition. Numerical solutions to this problem were obtained by employing the NEP theory that involves minimizing the following potential energy [90]:

$${\mathcal{V}}_{\mathrm{NEP}}=\int \frac{h}{2}{A}_{\mathrm{NEP}}^{\alpha \beta \gamma \delta }{\overline{E}}_{\alpha \beta }{\overline{E}}_{\gamma \delta }+\frac{h^3}{24}{A}_{\mathrm{NEP}}^{\alpha \beta \gamma \delta }{\overline{K}}_{\alpha \beta }{\overline{K}}_{\gamma \delta }d\overline{a}$$

(59)

in which $A_{\mathrm{NEP}}^{\alpha \beta \gamma \delta }=\frac{E}{2(1+\nu )}[\frac{2\nu }{1-\nu }{\overline{a}}^{\alpha \beta }{\overline{a}}^{\gamma \delta }+({\overline{a}}^{\alpha \gamma }{\overline{a}}^{\beta \delta }+{\overline{a}}^{\alpha \delta }{\overline{a}}^{\beta \gamma })]$ is the elasticity tensor for NEP theory, and $d\overline{a}=\sqrt{det{\overline{a}}_{\alpha \beta }}d{\xi }^{1}d{\xi }^2$ is the area element of a conceptually stress-free plate.

While the integrand in the potential energy (59) of the NEP theory is similar to our equivalent monolayer shell model with equation (22) using the natural fundamental forms ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ , there are several key differences. First, the potential energy of the NEP theory is obtained through integration over the conceptual stress-free configuration as can be seen in equation (59), which implies that for a total Lagrangian formulation, the resultant mid-surface stresses of this theory have a pre-factor ${\overline{J}}_a=\sqrt{\frac{det{\overline{a}}_{\alpha \beta }}{det{a}_{\alpha \beta }}}$ as ${\tau }_{\sigma }^{\alpha \beta }={\overline{J}}_a{\mathrm{hA}}_{\mathrm{NEP}}^{\alpha \beta \gamma \delta }{\overline{E}}_{\gamma \delta }$ and ${\tau }_M^{\alpha \beta }={\overline{J}}_a\frac{h^3}{12}{A}_{\mathrm{NEP}}^{\alpha \beta \gamma \delta }{\overline{K}}_{\gamma \delta }$ . This pre-factor ${\overline{J}}_a$ always exists regardless of the stimulus type, which may be incorrect as we discussed in Section 3 in terms of the mass change effect. In contrast, our equivalent monolayer shell model utilizes the equivalent material properties that account for the effect of mass addition in terms of biological growth on the stress constitution, such that the stiffness difference between the stimulus types with respect to the mass addition can be described. Second, while our proposed theory utilizes the initial configuration as a stress-free state, the NEP potential energy is formulated using the fictitious configuration as its stress-free state, and this makes it hard to derive physical interpretations of the role of non-mechanical stimuli in mechanical deformations and instabilities. Moreover, since the integration of potential energy in the NEP theory and its elasticity tensor are based on the natural metric tensor ${\overline{a}}_{\alpha \beta }$ , interpretation of stress states and a linear stability analysis with the NEP theory are challenging even if the natural metric tensor has a simple form.

Based on the NEP theory and assuming axisymmetry, the value of the buckling transition was numerically found to be $h_b\approx 0.3$ for incompressible materials when $\pi E/[4(1-{\nu }^2)]=1$ and the annulus geometry was $R_{\mathrm{in}}=0.1$ and $R_{\mathrm{out}}=1.1$ [90]. Here, we will test our SRS theory under pure natural stretch based on this annulus problem, for both the performance of the proposed novel SRS theory and its ability to examine the onset of buckling instability through a linear stability analysis.

We solve this problem via our proposed novel SRS theory by inserting the prescribed natural fundamental forms of equation (58) and ${\overline{b}}_{\alpha \beta }=0$ in equation (44). Our numerical solutions give the value of buckling transition as $h_b=0.28$ , which is in good agreement with the previous work. In addition, comparing the energies of equation (59) using the deformation solutions obtained via our proposed SRS theory, it is clear that the proposed theory accurately captures the evolution of all energies depending on the deformation as shown in Figure 5.

Furthermore, the solution of the membrane and bending stresses can be analytically obtained by solving equation (45) with equation (46) under axisymmetric conditions using the polar coordinate system relying on the initial geometry as

$$\begin{array}{c}\frac{{\stackrel{°}{\tau }}_{{\sigma }_{*}}^{11}}{Y_{\mathrm{eq}}}=\frac{\nu [0.25(1-\cos (2{\xi }^1))-{\xi }^1\sin ({\xi }^1)\cos ({\xi }^1)-\frac{0.255}{\nu }{\xi }^{1^2}+{\xi }^{1^2}\mathrm{sinc}(2{\xi }^1)]}{{\xi }^{1^2}}\hfill \\ \frac{{\stackrel{°}{\tau }}_{{\sigma }_{*}}^{22}}{Y_{\mathrm{eq}}}=\frac{{\nu }^2[0.25(1-\cos (2{\xi }^1))-{\xi }^1\sin ({\xi }^1)\cos ({\xi }^1)-\frac{0.255}{\nu }{\xi }^{1^2}+{\xi }^{1^2}\mathrm{sinc}(2{\xi }^1)]}{{\xi }^{1^4}}\hfill \\ {\stackrel{°}{\tau }}_{{\sigma }_{*}}^{12}={\stackrel{°}{\tau }}_{{\sigma }_{*}}^{21}={\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }=0,\hfill \end{array}$$

(60)

with the stretching rigidity $Y_{\mathrm{eq}}$ , where ${\xi }^1$ is the radial coordinate, and $\mathrm{sinc}(·)$ is the sine cardinal function. Here, as the annulus is a monolayer surface, we denote Poisson’s ratio in equation (60) without the subscript $\mathrm{eq}$ while keeping it on the membrane rigidity for the sake of generality in terms of the effect of mass addition (i.e. for stimulus types without mass addition, $Y_{\mathrm{eq}}\to Y$ also). We approximate the buckled shape by assuming that the components of $\delta \mathit{u}=\delta u_{\alpha }{\stackrel{°}{\mathit{a}}}^{\alpha }+\delta w\stackrel{°}{\mathit{n}}$ have forms as

$$\delta u_{\alpha }=0,\delta w=C\cos n_1\frac{\pi }{2}\frac{{\xi }^1}{R_{\mathrm{out}}},$$

(61)

where $C$ is a constant real number and $n_1$ is an integer. Then, the Donnell–Mushtari–Vlasov theory leads to the linear and second-order terms of the variations $\delta a_{\alpha \beta }$ and $\delta b_{\alpha \beta }$ [126] as

$$\delta a_{\alpha \beta }^{\mathrm{lin}}=\frac{1}{2}({\stackrel{°}{\nabla }}_{\alpha }\delta u_{\beta }+{\stackrel{°}{\nabla }}_{\beta }\delta u_{\alpha }),$$

(62a)

$$\delta b_{\alpha \beta }^{\mathrm{lin}}={\stackrel{°}{\nabla }}_{\alpha }{\stackrel{°}{\nabla }}_{\beta }\delta w,$$

(62b)

$$\delta a_{\alpha \beta }^{\mathrm{nl}(2)}=\frac{1}{2}{\nabla }_{\alpha }\delta w{\stackrel{°}{\nabla }}_{\beta }\delta w,$$

(63a)

$$\delta b_{\alpha \beta }^{\mathrm{nl}(2)}={\stackrel{°}{\nabla }}^{\gamma }\delta u_{\gamma }{\stackrel{°}{\nabla }}_{\beta }{\stackrel{°}{\nabla }}_{\alpha }\delta w-{\stackrel{°}{\nabla }}^{\gamma }\delta w{\stackrel{°}{\nabla }}_{\beta }{\stackrel{°}{\nabla }}_{\alpha }\delta u_{\gamma }$$

(63b)

such that the resulting balance equation obtained via the first variation of equation (48) with the integration domain ${\xi }^1:[0.1,1.1]$ and ${\xi }^2:[0,\pi /2]$ leads to the buckling transition thickness with $n_1=1$ as $h_b\approx 0.27$ which agrees well with both our numerical result and the previous NEP-based work.

Figure 5.

Annulus disk under in-plane differential natural stretch. $\mathcal{V}{}_{NEP}^s$ is the stretching energy part of the first term in equation (59), whereas $\mathcal{V}{}_{NEP}^b$ is the bending energy part of the second term. Symbols denote IGA simulation results obtained via the proposed stimuli-responsive shell theory, while solid lines mean the reference numerical solutions in Efrati et al. [90]. Vertical dashed lines at $h_b$ denote the buckling transition thickness obtained by the proposed stimuli-responsive shell theory.

The assumption of axisymmetry in this problem is made for the purpose of validation as the reference solution of Efrati et al. [90] artificially requires the axisymmetric constraint; however, we note that our numerics indicate that if the axisymmetric constraint is relaxed, the dome-like shape is no longer observed below a critical thickness.

7. Examples

In the previous section, we validated the proposed SRS theory and demonstrated that the effective external terms representing the effects of non-mechanical stimuli on deformation can capture fundamental cases involving either stretching or curvature (bending) dominated deformations. In this section, we apply the SRS theory to more complex examples considering shape changes induced by various types of stimuli, including biological growth, swelling, and applied electric voltage. One key goal of this section is to demonstrate how, using experimentally measured data for shells deformed by non-mechanical stimuli, we can obtain the corresponding stress-free 3D metric tensor of each layer, and thus the specific forms of the natural fundamental forms ${\overline{a}}_{\alpha \beta }$ and ${\overline{b}}_{\alpha \beta }$ to determine the resulting effective external loadings that deform the body. This solution process is shown in the flowchart in Figure 6.

Figure 6.

Flowchart for the process to obtain effective external loadings induced by non-mechanical stimuli, based on experimentally measured or theoretically calculated data. These external loadings can then be used to solve deformation problems via the proposed stimuli-responsive shell theory.

7.1. Long leaf growth

In the first example problem, we consider the biological growth of long leaves. We obtain numerical and analytical predictions of the shape-morphing process through the SRS theory and demonstrate the ability to capture large shape changes resulting from geometric nonlinearity. Morphogenesis in plants is primarily driven by non-homogeneous growth of tissues. Especially for long leaves, which are thin, the role of excessive cell growth within their marginal region is dominant in the shape formation in leaves. Here, we focus on the helical growth that can be commonly seen in nature, such as Dendrobium helix and Prosthechea cochleata.

Because helical growth is induced by twisting during growth, we impose a non-homogeneous growth stimulus on a monolayer thin leaf of dimensions $L=5$ cm and $W=0.5$ cm with thickness $h=200m$ (see Figure 7(a)) by implementing a width-direction differential expansion on the leaf, with its material properties of $E=350$ kPa and $\nu =0.25$ [141]. To obtain the amount of expansion induced by growth, we follow a previous work [141] that experimentally measured the width-direction differential expansion by dissecting the twisting croton mammy of Codiaeum variegatum leaf into several length-direction strips, in which the length change after the dissection was measured as the functional form $(l_{\mathrm{diss}}-L)/L=\beta (y/W{)}^2$ where $l_{\mathrm{diss}}$ is the strip length after the dissection and $\beta$ provides the maximum growth value. Therefore, the twisting leaf growth is modeled with the experimentally derived expansion factor $e(y)=1+\beta (y/W{)}^2$ , based on the concept of expansion factor $e=l_{\mathrm{diss}}/L$ , via an isotropic area expansion induced by growth as

$${\tilde{g}}_{\alpha \beta }=e^2{\stackrel{°}{g}}_{\alpha \beta }$$

(64)

which is plugged into equations (25) and (29) to calculate the natural fundamental forms and material properties. Here, we account for mass change from growth and its effect on the internal stress by considering the leaf growth as density-preserving in which the newly added mass during growth has the same material properties as the existing material, that is, ${\rho }_r=J_n=e^2$ as discussed in Section 3. Then, we numerically solve equation (44) based on the calculated natural fundamental forms and material properties. The current approach thus differs from the previous work, which modeled the leaf growth process as an equivalent thermal expansion problem and did not account for the effect of mass addition during growth [141]. Specifically, the stress amplification due to mass addition induces a higher contrast in the effect of the width-direction differential growth, which eventually leads to more twisted leaf shapes as a result.

Figure 7.

Growing long leaf. (a) Initial geometry with length $L$ and width $2W$ . (b) Shape change process as the growth amount increases, from the initial flat shape to twisted and helical shapes. (c) Twist angle per unit length. Symbols denote IGA simulation results obtained by the proposed SRS theory. The blue solid line is the analytical solution of equation (94) via our proposed SRS theory, whereas the black solid line is the analytical solution as ${\alpha }_{t}W=\sqrt{2\beta +{\beta }^2}$ in Huang et al. [141] solely from geometrical considerations.

For progressive leaf growth with increasing $\beta$ , the shell exhibits two distinct shape-shifting events as shown in Figure 7(b), showing the ability to capture significant geometric nonlinearity with the proposed SRS theory. First, the leaf morphs from an initially flat shape to one that becomes gradually more twisted. As the amount of growth increases, the twisted leaf shape morphs to the helical shape around $\beta =0.25$ , which is consistent with the previous study of Huang et al. [141].

During twisting, the twisting angle per unit length is analytically calculated depending on the amount of growth. Assuming that the deformed shape has the position vector of the twisted shape as $r=({\xi }^1,{\xi }^2\cos {\alpha }_t{\xi }^1,{\xi }^2\sin {\alpha }_t{\xi }^1)$ , we obtain the current fundamental forms as

$$\left[a_{\alpha \beta }\right]=\left[\begin{array}{cc}1+{\xi }^{2^2}{\alpha }_t^2\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right],\left[b_{\alpha \beta }\right]=\left[\begin{array}{cc}0\hfill & \frac{{\alpha }_t}{1+{\xi }^{2^2}{\alpha }_t^2}\hfill \\ \multicolumn{1}{c}{\frac{{\alpha }_t}{1+{\xi }^{2^2}{\alpha }_t^2}}& 0\hfill \end{array}\right]$$

(65)

in which ${\xi }^1\in [0,L]$ , ${\xi }^2\in [-W,W]$ , and ${\alpha }_t$ is the twisting angle per unit length. Then, inserting the current and natural fundamental forms into equation (35) gives the mechanical potential energy as

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{V}=\int \frac{{{\alpha }_t}^4{W}^8+2\beta (1+\nu )({\xi }^{2^2}\beta +2W^2)[2\beta ({\xi }^{2^2}\beta +2W^2)-{{\alpha }_t}^2{W}^4]}{{\left[{\xi }^{2^4}\mathrm{Eh}{({\xi }^{2^2}\beta +W^2)}^2\right]}^{-1}8(1-{\nu }^2)W^{12}}d\stackrel{°}{a}}\\ \multicolumn{1}{c}{=\frac{\mathrm{EhLW}({\mathcal{X}}_{l1}+{\mathcal{X}}_{l2}+{\mathcal{X}}_{l3})}{180180(1-{\nu }^2)}}\\ \multicolumn{1}{c}{\dots \{\begin{array}{c}{\mathcal{X}}_{l1}=143[5\beta (7\beta +18)+63]{{\alpha }_t}^4{W}^4\hfill \\ \multicolumn{1}{c}{{\mathcal{X}}_{l2}=-26\beta [5\beta \{7\beta (9\beta +44)+495\}+1386]{{\alpha }_t}^2(1+\nu )W^2}\\ \multicolumn{1}{c}{{\mathcal{X}}_{l3}=4{\beta }^2[5\beta \{7\beta (99{\beta }^2+702\beta +1859)+15444\}+36036](1+\nu )}\end{array}}\end{array}$$

(66)

with $O(h^3)$ , which leads to the minimum potential energy via $\partial \mathcal{V}/\partial {\alpha }_t=0$ with respect to the twisting angle per unit length as

$${\alpha }_{t}W=\sqrt{\frac{(1+\nu )\beta (1386+5\beta (495+7\beta (44+9\beta )))}{693+55\beta (18+7\beta )}}.$$

(67)

This analytical solution agrees well with the numerical simulation results as shown in Figure 7(c). Note that our analytical equation for the twisting angle ${\alpha }_t$ via the proposed SRS theory gives similar results to the one proposed of ${\alpha }_{t}W=\sqrt{2\beta +{\beta }^2}$ in Huang et al. [141], in which the best fit of Huang et al. [141] compared to equation (67) occurs for $\nu =0$ since the previous work derived the proposed equation solely from geometric considerations, while our solution is obtained via full elasticity with material properties.

7.2. Swelling-induced closure of Venus flytrap

For our next example, we examine the ability to capture the effects of a different stimulus, that is, swelling, using the proposed SRS theory. Specifically, we examine the swelling-induced closure of a flytrap, which also demonstrates that our proposed theory can accurately capture elastic instabilities that are exploited in nature to realize rapid motion and shape changes. A plant’s leaf blades or flower petals are mainly composed of mesophyll cells sandwiched between two epidermis layers. The differential expansion between the cells on the two different sides of the leaf blades in response to non-mechanical stimuli such as light and humidity dictates the opening and closure of plant tissues. Whereas the flower petal can open only once, many insectivorous plants including the Venus flytrap can repeatedly open and close their leaves during their lifetime. Furthermore, the leaf closure occurs rapidly via snap-through instability induced by differential swelling. Importantly, as we detail below, we model the snap-through process using experimentally-measured data [134] for the differential expansion (induced by swelling) occurring on opposite sides of the flytrap leaves to determine the resulting natural stretch and curvature that the flytrap leaves are subject to. This is in contrast to previous studies [133] where only the natural curvature is prescribed without linkage between the natural curvature and actual response in flytrap leaves to swelling.

In the cross-sectional cut perpendicular to the midrib of the flytrap leaves, the presence of microfibrils running perpendicular to the midrib has been observed in the upper epidermis. Due to the reinforcement of the cell wall by the microfibrils, the expansion of the upper side of flytrap induced by a stimulus is prevented during the flytrap leaf movement. Following this, a biological experiment was carried out to obtain the expansion factor for the upper and lower sides of the flytrap leaf during the closure movement, for the directions parallel and perpendicular to the midrib [134]. In this experiment, ink dots were drawn on the uppermost and lowermost superficial surfaces of flytrap leaves, and the distance change between ink dots during the closure movement was measured as $(d_{\mathrm{ink},a}-d_{\mathrm{ink},b})/d_{\mathrm{ink},b}$ , where $d_{\mathrm{ink},b}$ and $d_{\mathrm{ink},a}$ are the distance between ink dots before and after the closure movement, respectively, as listed in Table 1 [134].

Table 1.

Distance change between adjacent ink dots on the Dionaea leaf surface during the closure movement [134].

		Distance change
Upper epidermis	Perpendicular	–0.01
	Parallel	0
Lower Epidermis	Perpendicular	+0.07
	Parallel	+0.04

Building on the experimentally observed data, we imposed the thickness-direction differential swelling on a tri-layer surface of thickness $h$ , consisting of the mesophyll and the upper and lower epidermis of thickness $h_{\mathrm{ep}}$ as shown in Figure 8(a), to mimic the flytrap closure using an orthotropic area expansion induced by swelling as

$$\left[{\tilde{g}}_{\alpha \beta }\right]=\left[\begin{array}{cc}{e}_{\perp }^2{\stackrel{°}{g}}_{11}\hfill & e_{\perp }{e}_{\parallel }{\stackrel{°}{g}}_{12}\hfill \\ \multicolumn{1}{c}{e_{\perp }{e}_{\parallel }{\stackrel{°}{g}}_{21}}& e_{\parallel }^2{\stackrel{°}{g}}_{22}\hfill \end{array}\right].$$

(68a)

Figure 8.

Snapping-induced closure of flytrap. (a) Trilayer surface of the flytrap leaf. (b) Barrel-like initial shape of the flytrap model with $h=0.5$ mm, $L_f=3.46$ cm, $R_f=2$ cm, and ${\theta }_f={60}^{°}$ . (c) IGA simulation results of the change in distance $d_{\mathrm{tip}}$ between two tips of the flytrap leaves normalized by its initial value $d_{\mathrm{tip}0}$ . (d) Average mean curvature $H_{\mathrm{avg}}$ normalized by its initial value $H_{\mathrm{avg}0}$ , compared with experimental data in Forterre et al. [133].

Each layer is allowed different expansion factors ( $e_{\perp }$ and $e_{\parallel }$ ) in the perpendicular and parallel directions to the midrib, where the expansion factor for each layer is prescribed based on the experimentally measured distance changes listed in Table 1 as

$$\{\begin{array}{c}{e}_{\perp }=1-0.01t_n,e_{\parallel }=1\mathrm{for}\mathrm{upper}\mathrm{epidermis}\hfill \\ \multicolumn{1}{c}{e_{\perp }=1,e_{\parallel }=1\mathrm{for}\mathrm{mesophyll}}\\ \multicolumn{1}{c}{e_{\perp }=1+0.07t_n,e_{\parallel }=1+0.04t_n\mathrm{for}\mathrm{lower}\mathrm{epidermis}}\end{array}$$

(68b)

in which $t_n=t/{\tau }_s:[0,1]$ is the normalization of time $t$ with respect to the experimentally observed flytrap closure time ${\tau }_s=0.8$ s [133]. Here, we assume that the uppermost and lowermost superficial surfaces are almost stress-free so that the experimental data of Table 1 can represent the stress-free 3D metric tensor ${\tilde{g}}_{\alpha \beta }$ of each upper and lower layer. Thus, these expansion factors in the perpendicular and parallel direction are plugged into equations (25) and (29) to calculate the natural fundamental forms and material properties and then to numerically solve equation (44). As the type of non-mechanical stimulus is swelling, there is no effect of added mass in this case, i.e. ${\rho }_r=1$ , which can also be interpreted by considering that the fluid flowing through the body during swelling does not have mechanical stiffness. As we simplify the initial geometry of flytrap using a barrel-like shape like the actual flytrap leaves as shown in Figure 8(b), with $E=10$ MPa and $\nu =0.5$ for all layers [133], our simulation results show that the swelling-induced flytrap closure indeed occurs as a fast leaf movement via snapping as shown in Figure 8(c). Because we are able to capture the experimentally measured timescale and change in mean curvature as shown in Figure 8(d), this indicates our SRS theory can accurately capture swelling-induced snapping of plant leaves.

7.3. Voltage-induced buckling of dielectric film

Our final example demonstrates the ability of the proposed novel SRS theory to capture a different stimulus, that is, the electric field–induced actuation of thin structures, including predicting the onset of electromechanical instabilities via a stability analysis. As the dielectric elastomer is composed of a film, which is electrically insulated and can be polarized via applying an electric field, between two electrodes, a Maxwell stress is developed in the film to the surface normal direction as the electrodes approach each other under an applied electric field [142, 143]. During this process, the film isotropically expands in the in-plane direction due to the repulsion of like charges along each surface of the film, which generates significant natural stretching, which motivates the application of the proposed novel SRS theory for this problem. Clamping the film generates compressive stresses by preventing the film from expanding in-plane, which leads to buckling of the film when the compressive stress exceeds a critical value. This buckling of dielectric elastomer has been used to pump fluids within microfluidic systems in response to an applied electric field as the dielectric elastomer can be actively and reversely buckled out of plane [19–20]. To enhance the functionality of microfluidic systems, controlling and directing fluid flow are essential [144–146]. Thus, estimation of the critical voltage that induces the critical compressive stress is crucial to predicting the onset of instabilities for the dielectric elastomer within microfluidic and other systems.

Our analysis of this problem is thus motivated by the significant natural stretching that occurs, and also because while the estimation of the critical voltage has been widely studied using standard solid continuum models [19, 20, 147, 148], the estimation of the critical buckling voltage has not, to the best of our knowledge, been derived based on a shell theory. Furthermore, we demonstrate the ability to analyze the stability of complex multi-physics phenomena using the natural stretch- and curvature-based SRS theory, in which the natural stretch for this problem is obtained from a theoretically calculated expansion factor based on the Maxwell stress. Since the buckling of dielectric elastomers occurs due to the developing compressive stress in the film, we focus on the positive-definite characteristic of the second variation (48) of the film’s mechanical potential energy.

The expansion factor depending on the applied voltage can be theoretically calculated using the Maxwell stress in the dielectric film with thickness $h$ , while the Maxwell stress induced by an applied voltage is calculated as ${\mathrm{\Sigma }}_{\mathrm{Max}}=-\frac{1}{2(1-{\nu }^2)}{ϵ}_{V0}{ϵ}_V{\left(\frac{V_e}{h}\right)}^2$ , where $\nu$ is the Poisson’s ratio of the film, $V_e$ is the applied voltage, and ${ϵ}_{V0}$ and ${ϵ}_V$ are, respectively, the free space permittivity and the relative dielectric constant [142, 149]. When a free-standing dielectric film is isotropically expanding to the in-plane direction while simultaneously reducing its thickness due to an applied voltage, the second Piola–Kirchhoff stress in the film with respect to deformations from the initial body can be calculated using the Saint-Venant energy model [150] as

$${\mathrm{\Sigma }}^{\mathrm{ij}}=[\lambda {\stackrel{°}{g}}^{\mathrm{ij}}{\stackrel{°}{g}}^{\mathrm{kl}}+\mu ({\stackrel{°}{g}}^{\mathrm{ik}}{\stackrel{°}{g}}^{\mathrm{jl}}+{\stackrel{°}{g}}^{\mathrm{il}}{\stackrel{°}{g}}^{\mathrm{jk}})]\frac{1}{2}(g_{\mathrm{kl}}-{\stackrel{°}{g}}_{\mathrm{kl}}),$$

(69a)

where $\lambda =\frac{E\nu }{(1+\nu )(1-2\nu )}$ and $\mu =\frac{E}{2(1+\nu )}$ are the Lamé constants of the dielectric film with the film’s Young’s modulus $E$ . Since the dielectric film is homogeneously expanding without any constraint, the current 3D metric tensor can be assumed with the expansion factor $e$ as

$$\left[g_{\mathrm{kl}}\right]=\left[\begin{array}{cc}{e}^2{\left[{\stackrel{°}{g}}_{\alpha \beta }\right]}_{2\times 2}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\lambda }_3^2\hfill \end{array}\right]$$

(69b)

in which ${\lambda }_3$ has the information about a change in thickness as stated in Section 2. Relying on the definition of Maxwell stress, equating ${\mathrm{\Sigma }}^{33}={\mathrm{\Sigma }}_{\mathrm{Max}}$ with ${\mathrm{\Sigma }}^{\alpha \beta }=0$ for the free-standing homogeneous expansion of the film gives the expansion factor that depends on the applied voltage as

$$e=\sqrt{1+\frac{\nu {ϵ}_{V0}{ϵ}_V{V}_e^2}{E{h}^2(1-{\nu }^2)}}.$$

(70)

The dielectric film geometry is that of a circular plate with radius $R$ and thickness $h$ . Building on the polar coordinate system of radial ${\xi }^1:[0,R]$ and polar ${\xi }^2:[0,\pi /2]$ coordinates, the solution of stress states on the fundamental (equilibrium) configuration of the clamped film can be calculated via equation (45) with equation (46) using the natural fundamental forms and material properties obtained by plugging

$${\tilde{g}}_{\alpha \beta }=e^2{\stackrel{°}{g}}_{\alpha \beta }$$

(71)

into equations (25) and (29), based on the isotropic area expansion induced by the applied voltage. As the type of non-mechanical stimulus is the voltage, there is no mass addition in this case, that is, ${\rho }_r=1$ . Then, solving equation (45) with equation (46) under clamped boundary conditions, while exploiting the axisymmetry, gives the solution of the membrane and bending stresses in the film caused by the applied voltage as

$$[{\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }-f_{\mathrm{\Lambda }}^{\alpha \beta }]=\left[\begin{array}{cc}-\frac{1}{2}\frac{\nu {ϵ}_{V0}{ϵ}_V{V}_e^2}{h{(1-\nu )}^2(1+\nu )}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -\frac{1}{2}\frac{\nu {ϵ}_{V0}{ϵ}_V{V}_e^2}{h{(1-\nu )}^2(1+\nu ){\xi }^{1^2}}\hfill \end{array}\right],{\stackrel{°}{\tau }}_{M_{*}}^{\alpha \beta }-f_{\kappa }^{\alpha \beta }=0.$$

(72)

In this example problem, following Timoshenko’s and Bulson’s [151, 152] works, we solve this stability problem under the fixed boundary condition, and the effect of clamped boundary condition on the buckling is dealt with by employing the compensating factor $k_{\mathrm{bc}}=14.68$ which is a constant for a clamped plate exhibiting the first-mode buckling, as $({\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }-f_{\mathrm{\Lambda }}^{\alpha \beta })\to ({\stackrel{°}{\tau }}_{{\sigma }_{*}}^{\alpha \beta }-f_{\mathrm{\Lambda }}^{\alpha \beta })/k_{\mathrm{bc}}$ . Similar to the annulus benchmark problem above, we approximate the buckled shape by assuming that the components of $\delta \mathit{u}=\delta u_{\alpha }{\stackrel{°}{\mathit{a}}}^{\alpha }+\delta w\stackrel{°}{\mathit{n}}$ have forms as

$$\delta u_{\alpha }=0,\delta w=C\cos n_1\frac{\pi }{2}\frac{{\xi }^1}{R},$$

(73)

where $C$ is a constant real number and $n_1$ is an integer. With the Donnell-Mushtari–Vlasov theory of equations (62) and (63), the resulting balance equation obtained via the first variation of equation (48) with the integration domain ${\xi }^1:[0,R]$ and ${\xi }^2:[0,\pi /2]$ leads to the critical voltage for buckling with $n_1=1$ as

$$V_{e_c}=\frac{h^2\pi \sqrt{\left(\frac{{\pi }^2-4}{6(4+{\pi }^2)}\right)\left(\frac{E{k}_{\mathrm{bc}}(1-\nu )}{\nu {ϵ}_V{ϵ}_{V0}}\right)}}{2R}.$$

(74)

This critical voltage gives $V_{e_c}=2.96$ kV for incompressible films with $E=1.2$ MPa, $R=400$ μm, $h=60$ μm, ${ϵ}_V=3.21$ , and ${ϵ}_{V0}=8.85\times {10}^{-12}$ N/V2, which is in agreement with the previous work of $V_{e_c}=3.0$ kV in Tavakol and Holmes [20] that solved this buckling problem using the existing closed-form of critical buckling stresses.

8. Conclusion

We have presented a general shell theory to account for non-mechanical stimuli. We show that regardless of the specific type of stimuli, they can be generalized into three particular forms: those that add mass to the shell, those that increase the area of the shell through the natural stretch $\mathrm{\Lambda }$ , and those that change the curvature of the shell through the natural curvature. $\kappa$ Considering these effects, we demonstrated how to form a mechanical interpretation of non-mechanical stimuli, in which the standard shell strain energy can be augmented by a potential corresponding to how a non-mechanical stimulus acts to change the shell’s area and a potential corresponding to a stimulus that acts to change the shell’s curvature. We also generalized the effect of mass added during biological growth on the stress constitution based on a mass balance during the growth process. By developing these potentials, we enable the usage of standard tools for the mechanics community, that is, taking variations to derive a weak form that can be used for finite element analysis and also to use second variations to perform linear stability analyses. The theory is formulated based on a standard, stress-free reference configuration that is known a priori, meaning that it can be physically observed, and only requires the solution of a single-field equation, the standard mechanical momentum or equilibrium equation, despite capturing the effects of non-mechanical stimuli. Furthermore, we derived analytical equations for an equivalent monolayer system to calculate its natural fundamental forms and material properties, which prescribe the shell deformation. For these equations, we can use experimental measurements or theoretical calculations of deformations under various stimuli in general multi-layer surfaces. Overall, we expect that mechanicians and non-mechanicians alike can use the approach presented here to modify existing computational tools with the effective external loadings calculated in this novel theory in order to study how various types of non-mechanical stimuli impact the mechanics and physics of thin-shell structures.