The Föppl–von Kármán equations of elastic plates with initial stress

Initially, stressed plates are widely used in modern fabrication techniques, such as additive manufacturing and UV lithography, for their tunable morphology by application of external stimuli. In this work, we propose a formal asymptotic derivation of the Föppl–von Kármán equations for an elastic plate with initial stresses, using the constitutive theory of nonlinear elastic solids with initial stresses under the assumptions of incompressibility and material isotropy. Compared to existing works, our approach allows us to determine the morphological transitions of the elastic plate without prescribing the underlying target metric of the unstressed state of the elastic body. We explicitly solve the derived FvK equations in some physical problems of engineering interest, discussing how the initial stress distribution drives the emergence of spontaneous curvatures within the deformed plate. The proposed mathematical framework can be used to tailor shape on demand, with applications in several engineering fields ranging from soft robotics to four-dimensional printing.


Introduction
The Föppl-von Kármán (FvK) equations are a set of nonlinear partial differential equations describing the large deflection of linear elastic plates [1,2]. They can be derived as a formal asymptotic expansion of the three-dimensional filed theory of linear elasticity in the limit of large displacements and small strains, and associated with specific boundary conditions [3][4][5][6][7].
The FvK equations are notably difficult to solve, but they proved to be very useful to give theoretical insights into many physical problems where in-plane and out-of-plane unknowns can be decoupled [8][9][10][11][12]. The research interest in FvK equations has been recently reinvigorated by the technological possibility to fabricate shape-morphing devices using soft active materials [13][14][15][16][17][18][19]. These morphable plates have applications in several engineering fields, ranging from soft robotics [20,21] to the design of biomimetic structures [16]. In particular, morphological transitions can be realized by controlling the geometric frustration of a soft plate by swelling [22][23][24], surface accretion [25], optothermal stimuli in nematic elastomers [26] and surface tension in nano-plates [27]. For these purposes, FvK equations have been derived in cases where geometrical incompatibilities arise and the undeformed configuration of the elastic plates is no longer free of initial stresses [28].
The formal asymptotic expansion leading to the FvK model has been rigorously derived as the G-limit of the three-dimensional elastic problem [29,30]. This analysis has been recently extended to the case of a pre-strained plate [31][32][33][34]. The existing approaches account for the geometrical frustration using additive or multiplicative decomposition of the deformation gradient for describing the spatial distribution of residual strains given by the underlying non-Euclidean metric [35][36][37][38]. However, the FvK equations with pre-strains require the prescription of the incompatible metric of the virtual relaxed state, while in many practical cases, the distribution of residual strains remains unknown. Since a stress-free configuration cannot be physically attained by non-invasive techniques, it is more suitable to consider a more general theoretical framework where the elastic energy explicitly depends on the spatial distribution of the initial stress in the reference configurations, without the need to prescribe a stressfree state [39].
In the following, we present a formal asymptotic derivation of the FvK equations for an elastic plate with initial stresses. In §2, we introduce the constitutive theory of nonlinear elastic solids with initial stresses under the assumptions of incompressibility and material isotropy. In §3, we introduce the scaling assumption for the geometrical parameters and the initial stress components and we derive the FvK equations for an initially stressed elastic plate. In §4, we explicitly solve the derived FvK equations in some physical problems, discussing how the initial stress concentration may drive the emergence of spontaneous curvatures within the deformed plate. Concluding remarks are finally summarized in §5.

Constitutive theory for initially stressed elastic materials
Let us consider a body that occupies a simply connected domain B t in its reference configuration. We use Cartesian unit vectors (E X , E Y , E Z ) and (e x , e y , e z ) in the reference and spatial configurations, respectively. Let X = XE X + YE Y + ZE Z be the material position vector. In this undeformed configuration, the body has an initial stress, meaning that its Cauchy stress tensor is not vanishing. We denote the initial stress tensor by t, where t : B t ! SðR 3 Þ, and SðR n Þ is the set of the self-adjoint L : R n ! R n , where L [ LðR n Þ is a linear application. In order to enforce the material balance of linear and angular momentum, this initial stress is such that where Div is the material divergence operator. We remark that if the body is residually stressed, i.e. the zero-traction boundary condition tN = 0 applies to the whole boundary @B t with material unit normal N, t must be inhomogeneous and have zero average over the volume in B t in force of the mean value theorem [40]. Let x = w(X) = xe x + ye y + ze z be the spatial material vector, so that w : B t ! B be the one-to-one mapping to the deformed configuration B and F ¼ @w=@X be the deformation gradient. In the following, we deal with incompressible materials, imposing the constraint J ¼ det F ¼ 1. We further assume that the body possesses a perfectly elastic response, defining a strain energy density C per unit of reference volume as [41,42] C ¼ CðF, tÞ: ð2:2Þ By standard arguments, the second Piola-Kirchhoff S and Cauchy σ stress tensors reads where p is the Lagrange multiplier enforcing the incompressibility constraint. In particular, when F is equal to the identity tensor I, the Cauchy stress must be equal to the initial stress tensor t, so that where p t is the value of p in B t . If we further assume that the elastic response is invariant after the application of a rigid-body motion, the strain energy can be expressed as a function of the three invariants of the right Cauchy-Green tensor C ¼ F T F, the three invariants of the initial stress tensor t, plus their four mixed invariants [42,43]. Since we are interested in developing a FvK theory, we are interested in considering the minimal constitutive response that takes into account geometric nonlinearities. Accordingly, we consider in the following the constitutive response of a pre-stressed Neo-Hookean material, that is given by [44] CðF, tÞ ¼ where E is the Young modulus of the unstressed material and r is the real root of with I t1 ¼ trt, I t2 ¼ 1 2 ½ðI 2 t1 À trðt 2 Þ, I t3 ¼ det t. Equation (2.6) is obtained by inverting the constitutive relation (2.3) and imposing the incompressibility constraint after a repeated application of the Cayley-Hamilton theorem. Enforcing the compatibility condition (2.4) for the initial stress, we also find that p t = r. Such an energy guarantees that the material properties are not affected by the initial stress distribution, for details see [45].
From (2.3), the constitutive equation for the initial stressed Neo-Hookean material with strain energy given by (2.5) is where B ¼ FF T is the left Cauchy-Green tensor.
Neglecting the presence of volume bulk forces, the equilibrium conditions in the reference configuration read while the zero-traction conditions at the boundary give SN ¼ 0 at @B t . In the following, we will use these constitutive assumptions to derive a FvK theory of elastic plates with initial stress.

The Föppl-von Kármán equations for initially stressed plates
In this section, we first make the geometric assumptions of the FvK theory, using dimensional analysis to justify the asymptotic development. We later derive the constitutive equations using a variational argument, assuming different scaling for the initial stress components.

Geometric assumptions and asymptotic analysis
We consider an elastic plate with reference configuration B t ¼ S m Â ½ÀH, H, where S m is a closed subset in R 2 , as sketched in figure 1. In particular, we require that there exists a point P ∈ S m such that Figure 1. Sketch of the reference (a) and current (b) configuration of the elastic plate with initial stress. The blue surfaces represent the mid-section of the plate.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 and L is much larger than the thickness 2H of the plate, so that e = H/L ≪ 1, where B L ðPÞ ¼ fY [ R 2 j kY À Pk , Lg. The dimensionless parameter e will be used to perform an asymptotic expansion. In this thin geometry, we introduce the following dimensionless variables Similarly, let r be the gradient operator with respect to the dimensionless variables. Let u = x − X be the displacement vector describing the deformation of the plate. We denote by u k the projection of u on the plane containing S m , while w indicates the orthogonal component of u. More explicitly where P ¼ I À E Z E Z is the projection operator, while u k and w are the dimensionless counterpart of u k and w, respectively. From now on, for a generic quantity f, we denote by f its dimensionless counterpart. We assume the classical Kirchhoff hypothesis for the in-plane displacement u k where U m ðX k Þ ¼ HU m ðX k Þ ¼ u k ðPXÞ represents the in-plane displacement of the point X k of the middle surface, obtained through the projection of X on S m , while r k is the gradient operator with respect to X k , namely in Cartesian components r k ¼ ð@=@X, @=@Y, 0Þ, while r ? ¼ ð0, 0, @=@ZÞ. From the nondimensionalization proposed in (3.1), we get The out-of-plane component of the displacement vector reads where j m ðX k Þ ¼ Hj m ðX k Þ ¼ wðPXÞ is the out-of-plane displacement of the point X k of the middle surface, obtained through the projection of X on S m and W : B t ! R 3 . The field ξ governs the out-of-plane deformation up to an order O(e 2 ).
The assumption made in (3.3) and (3.5) imposes that the displacement is a pure bending deformation in the Z-direction, so that sections that are initially perpendicular to the middle surface undergo a rotation that is driven by the local curvature of the deformed middle surface. For an initially stressed plate, the validity of this scaling requires a further constraint on the scaling of the initial stress tensor, that will be discussed later.
We further assume that the middle surface components can be expanded in powers of e as follows: Using (3.4), the deformation gradient can be written as We can compute r k u and r ? u, obtaining where comma denotes the partial derivative. Thus, the right Cauchy-Green strain tensor becomes where the tensor A is defined as for the sake of simplicity. At the leading order, the incompressibility constraint reads tr A ¼ 0, that imposes where D ¼ r Á r is the Laplace operator. The tensor A affinely depends on Z. More explicitly, we write In appendix A, we report the explicit expression of the tensors in (3.14) using the canonical basis (E X , E Y , E Z ). Finally, for future convenience, we also recall that the asymptotic expansion of the inverse of C reads

Scaling assumptions on the initial stress
Under the constitutive assumption given by (3.5), we define the dimensionless initial stress tensor as t ¼ t=E. We assume the following scaling for the initial stress tensor: and If the top and bottom surfaces are free of traction, i.e. t iZ ðX, Y, +1Þ ¼ 0, the previous equation can be further simplified by integration in the Z-direction 19Þ royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 We define the averaged planar initial stress tensor t m as The leading order of (3.19) can be automatically fulfilled by introducing the dimensionless initial Airy stress function x 0 : S m ! R, so that where r m is the dimensionless gradient operator in R 2 , cofA is the cofactor of the tensor A. For the sake of clarity, using Cartesian coordinates, we get Under these assumptions, we get . From (2.6), the leading order expression of the dimensionless term r ¼ r=E is Using (2.7), (3.15) and (3.16), the components of the dimensionless second Piola-Kirchhoff stress tensor S ¼ S=E scales as follows: where the latter scaling is enforced by imposing that the dimensionless Lagrange multiplier p ¼ p=E reads In particular, we remark that (3.24) corresponds to the classical FvK scaling for the stress tensor, and justifies the mathematical soundness of the Kirchhoff assumption for the displacement field in (3.3) and (3.5) in the presence of initial stress components which scale as in (3.16). The leading-order components of the stress tensor read We remark that both S 0 and S 1 are independent of Z, since they depend on X k through j and U.

Variational formulation
Let E ¼ E=ðEHL 2 Þ be the dimensionless counterpart of the elastic energy. From now on, we will refer only to dimensionless variables (unless differently specified) and we drop the use of the superposed line to indicate the non-dimensional quantities. We now perform a variational derivation of the generalized FvK theory for an initially stressed material by imposing the stationary conditions for the total elastic energy functional E. Neglecting traction loads at the boundary and body forces for the sake of simplicity, the first variation of the dimensionless energy reads where dS = dXdY and A : B ¼ trðA T BÞ. From (3.10), (3.13) and (3.14), the increment of the right Cauchy-Green tensor C reads royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 Since S 0 and S 1 are symmetric and S 0 E Z ¼ S 1 E Z ¼ 0, we can introduce the tensors S m0 and S m1 which are their projection onto SðR 2 Þ. Substituting (3.26) and (3.29) into (3.28), we get where dE s and dE b are the increments due to average planar stretching and bending of the plate, respectively. They are defined as where M is the dimensionless tensor representing the average bending torques imposed by the initial stress, defined as which vanishes if t αβ are even functions of Z, with α, β ∈ {X, Y}.
In the following, we detail such contributions, deriving the corresponding FvK equations as the necessary conditions for their extremal values.

Föppl-von Kármán equation for the average planar stretch
Integrating by parts (3.31) and neglecting the remainder, we get where D m ¼ r m Á r m is the Laplace operator in R 2 and the bracket operator is the Monge-Ampére bilinear form

Föppl-von Kármán equation for the average bending stretch
Using (3.27) in (3.32) and neglecting the remainder, we integrate by parts, obtaining By substituting (3.35) into (3.38) and imposing dE b ¼ 0 for each admissible variation dξ, we finally get where D m C M ¼ 9 4 r m Á ðr m Á MÞ, so that C M represents the spontaneous mean curvature imposed by the 1 The existence of the Airy stress function can be shown by using some standard theorems of the theory of distributions (see [46], p. 59, Théoréme VI). For multiply connected domains, one could proceed following the procedure proposed in ( [7], pp. 61-66).
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 initial stress distribution. In appendix A, we show the explicit expressions of the FvK equations in Cartesian and polar coordinates.

Physical examples
In

Planar initial stress
Let us start by considering a planar distribution of the initial stress. Since t does not depend on Z, the corresponding average bending torques vanish, i.e. M ¼ 0. Accordingly, such an initial stress distribution may impose only a non-zero spontaneous Gaussian curvature C G in (3.36), while C M = 0 in (3.39). For the sake of clarity, we discuss in the following some physical examples exhibiting different spontaneous Gaussian curvatures.

Positive spontaneous Gaussian curvature
Let R, Q be the polar coordinates of the generic point X k [ S m . Let us first consider an initial stress distribution having an Airy stress function given by with c being a characteristic dimensionless stress parameter, and the planar initial stress tensor t m reads From equation (3.36), the spontaneous Gaussian curvature is positive and constant, being We look for a solution such that the second Piola-Kirchhoff stress tensor vanishes. In such a case, equations (3.36)-(3.39) simplify as where the former is a Monge-Ampére equation. Using polar coordinates, we get A solution to equation (4.4) such that ξ depends only on R is given by which corresponds to a buckling of the plate, as illustrated in figure 2.

Negative spontaneous Gaussian curvature
We now consider the following Airy stress function with c being a characteristic dimensionless stress parameter. Accordingly, the planar initial stress reads royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 From (3.36), the spontaneous Gaussian curvature is negative, being It is easy to check that these equations admit a solution for each n [ N of the form jðR, QÞ ¼ cR nþ2 sinððn þ 2ÞQÞ n 2 þ 3n þ 2 : ð4:11Þ In particular, the solution for n = 0 corresponds to a twisting of the plate, as illustrated in figure 3.

Concentrated Gaussian curvature in thin plates
In this subsection, we consider the case of large deflections ξ, so that the bending term D 2 m j in equation (3.39) can be neglected [36]. In this case, if χ = 0 the second FvK equation is automatically satisfied while the first FvK equation (3.36) becomes the Monge-Ampére equation We now consider the initial Airy stress function which corresponds to the planar initial stress tensor royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 9 If the second Piola-Kirchhoff stress tensor vanishes, the first FvK equation (3.36) becomes

Inhomogeneous spontaneous Gaussian curvature
Let us now consider an initial Airy stress function given by with c being a characteristic dimensionless stress parameter and k > 0 a characteristic dimensionless decay length. It corresponds to a uniaxial compression along the Y-axis that is inhomogeneous along the X-axis, and the only non vanishing initial stress component is Assuming variable separation, we set ξ = f (Y ) e −kX so that (4.22) can be transformed in the following nonlinear ordinary differential system: 16c 2 k 2 À f 02 þ f f 00 ¼ 0 and k 4 f þ 2k 2 f 00 þ f 0000 ¼ 0: ð4:23Þ A solution of (4.23) is given by f = 4c sin (kY), so that the vertical deflection is j ¼ 4c sinðkYÞ e ÀkX , ð4:24Þ that corresponds to a sinusoidal oscillation along the Y-axis that decays exponentially along the X-axis with a decay length k and an amplitude proportional to c. The resulting morphology of the middle plane of the initially stressed plate is depicted in figure 5.

Initial stress varying along the plate thickness
We finally consider the case in which the initial stress is allowed to vary along the plate thickness. For the sake of simplicity, we take t ¼ tðX, ZÞ, with t YY = t XY = t YZ = 0. In this case, from (3.18), we can define an Airy stress function F ¼ FðX, ZÞ such that t XX ¼ F ,ZZ , t XZ ¼ ÀF ,XZ and t ZZ ¼ F ,XX : ð4:25Þ In particular, assuming that t XZ | Z±1 = t ZZ | Z±1 = 0, a leading order expansion of F along the vertical direction gives where a 0 , a 1 and b are dimensionless stress parameters. The only non-zero components of the average royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 11 initial stress tensor t m and torque tensor M are the following: Accordingly, the parameter b determines the characteristic amplitude of the average initial stress, while a 0 , a 1 determine the magnitude of the initial bending torque. Since t m XX is a constant, the initial stress distribution has zero spontaneous Gaussian curvature, i.e. C G = 0. Assuming ξ = ξ(X ), since [ξ, ξ] = 0 equation (3.36) is automatically satisfied if χ = χ 0 = bY 2 . The torque distribution in (4.27) imposes a spontaneous average curvature C M ¼ 3 2 ða 1 X þ a 0 Þ, and equation (3.39) reads 2 9 j 00 À bj À 1 3 ða 1 X þ a 0 Þ ¼ 0: ð4:28Þ If the average initial stress is zero, i.e. if b = 0, then the solution corresponds to a bending deflection having exactly the curvature C M , with where C 1 , C 2 are two constants of integration that must be fixed by boundary conditions. If the average initial stress does not vanish, i.e. b ≠ 0, the deflection reads ð4:30Þ If the average uniaxial stress is tensile, i.e. b > 0, than the deflection has an exponential trend; if it is compressive, i.e. b < 0, than the deflection has the characteristic sinusoidal behaviour of an Euler buckling where the symmetry is broken by a spontaneous average curvature. Two possible bending morphologies of the middle line of the plate are depicted in figure 6. Thus, in the limit of narrow plates, the governing equation (4.28) recovers the Euler beam theory for elastic rods, where the average initial stress t m XX appears as a traction load, and M XX as a distributed torque provoking a spontaneous curvature.

Concluding remarks
In this work, we have derived the FvK equations for an elastic plate with initial stress. The reference configuration of the plate is a parallelepiped, whose thickness 2H along the z-direction is much smaller than the characteristic length L of its edges. We have identified with e = H/L the dimensionless parameter describing the thinness of the plate. Adopting the theoretical framework of initially stressed materials developed in [41], the elastic energy of the plate is obtained by an asymptotic expansion with respect to the small parameter e. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 More explicitly, we have assumed that the plate is composed of an incompressible neo-Hookean material, whose strain energy depends on both the deformation gradient and the initial stress, see (2.5). Assuming the classical Kirchhoff hypothesis for the displacement of the plate (see (3.3)-(3.5)), we used the scaling for the initial stress components reported in (3.16).
Under these assumptions, we have obtained the balance equation through a variational approach, by enforcing that the energy functional be stationary. The equilibrium equations (3.36) and (3.39) generalize the FvK equations for the average planar and bending stretch, respectively, to the case of initially stressed plates.
We have solved the FvK equations in some physical examples of engineering interest. By tuning the initial stress distribution within the plate, we have obtained buckled configurations exhibiting a constant positive or negative Gaussian curvature. Furthermore, it is possible to obtain surfaces with non-constant Gaussian curvature by properly tuning the initial stress of the plate. In particular, we propose some new solutions with a radially inhomogeneous Gaussian curvature. We have also recovered some known explicit solutions of the FvK equations, such as the conical solutions proposed by Ben Amar & Pomeau [8] for thin plates, where the initial stress exhibits logarithmic singularities, and the solution for a geometrically frustrated FvK plate reported in [28]. Finally, we analysed the effect of an initial stress that varies along the plate thickness. The plate spontaneously bends when the initial stress along the X-direction is non-negative, while it exhibits a wrinkling pattern if it is compressive, similarly to Euler buckling of elastic beams. In the case of more complex initial stress distributions, royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220421 e.g. exhibiting checkerboard or labyrinth patterns, it is unlikely that analytical solutions can be found and approximated solutions of the proposed equations may be found exploiting numerical techniques, such as the finite-element method.
In conclusion, we have derived the FvK equations for an initially stressed plate using a formal dimensional reduction of a nonlinear strain energy function depending explicitly on both the deformation gradient and the initial stress tensor. The main advance of the proposed approach is to unravel the effects of the initial stress distribution on the spontaneous average and Gaussian curvatures of the plate without the need to prescribe incompatible pre-strains, as required in earlier works [28,35]. In fact, unlike pre-straining, the initial stress distribution within the body can be measured by means of non-destructive techniques, such as ultrasound elastography [48] or photoelasticity [49]. Moreover, a target initial stress distribution can be physically realized using modern digital fabrication techniques, such as four-dimensional printing [18,25] and UV lithography [50]. Thus, the results of the proposed model may be used to design residually stressed objects in proximity to an elastic bifurcation [51], with the aim of fabricating shape-shifting plates able to adapt their morphology in the presence of external stimuli, such as a chemical potential [52] or an electric field [53].
Data accessibility. This article has no additional data.