Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3627422

Formats

Article sections

- Abstract
- 1. Introduction
- Kinematics of growth
- 3. Balance equations of growth
- 4. Constitutive equations
- 5. Discretization
- 6. Results
- 7. Discussion
- 8. Conclusion
- References

Authors

Related links

J Mech Phys Solids. Author manuscript; available in PMC 2014 June 1.

Published in final edited form as:

J Mech Phys Solids. 2013 June 1; 61(6): 1446–1463.

Published online 2013 February 8. doi: 10.1016/j.jmps.2013.01.007PMCID: PMC3627422

NIHMSID: NIHMS443042

phone: +1.650.450.0855, fax: +1.650.725.1587, url: http://biomechanics.stanford.edu, corresponding author

See other articles in PMC that cite the published article.

Many biological systems are coated by thin films for protection, selective absorption, or transmembrane transport. A typical example is the mucous membrane covering the airways, the esophagus, and the intestine. Biological surfaces typically display a distinct mechanical behavior from the bulk; in particular, they may grow at different rates. Growth, morphological instabilities, and buckling of biological surfaces have been studied intensely by approximating the surface as a layer of finite thickness; however, growth has never been attributed to the surface itself. Here, we establish a theory of continua with boundary energies and growing surfaces of zero thickness in which the surface is equipped with its own potential energy and is allowed to grow independently of the bulk. In complete analogy to the kinematic equations, the balance equations, and the constitutive equations of a growing solid body, we derive the governing equations for a growing surface. We illustrate their spatial discretization using the finite element method, and discuss their consistent algorithmic linearization. To demonstrate the conceptual differences between volume and surface growth, we simulate the constrained growth of the inner layer of a cylindrical tube. Our novel approach towards continua with growing surfaces is capable of predicting extreme growth of the inner cylindrical surface, which more than doubles its initial area. The underlying algorithmic framework is robust and stable; it allows to predict morphological changes due to surface growth during the onset of buckling and beyond. The modeling of surface growth has immediate biomedical applications in the diagnosis and treatment of asthma, gastritis, obstructive sleep apnoea, and tumor invasion. Beyond biomedical applications, the scientific understanding of growth-induced morphological instabilities and surface wrinkling has important implications in material sciences, manufacturing, and microfabrication, with applications in soft lithography, metrology, and flexible electronics.

Just like skin covers and protects the outer surface of our body, the mucous membrane covers and protects our inner body surfaces. The mucous membrane is a thin, moist layers of epithelial cells, which lines the air-organ interface of our respiratory, digestive, and urogenital tracts, and of our mouth and nose. It is attached to a soft submucosal layer supported by an outermost layer of smooth muscle cells. Surface wrinkling of the mucous membrane is critical for healthy biological function, and abnormalities in wrinkling patterns are a classical hallmark of disease [60]. A typical example is severe asthma, a chronic disease of the respiratory system, for which there is currently no cure. As illustrated in Figure 1, during chronic asthma, the airway wall undergoes an irreversible structural remodeling, associated with an increased deposition of collagen, a permanent thickening of the smooth muscle layer, growth of the mucous membrane, and mucosal folding [19]. Collectively, these changes induce an excessive airway wall narrowing, resulting in progressively reduced lung function [8].

Chronic airway wall remodeling and characteristic structural changes. Normal small airway (left). Acutely inflamed airway with thickened airway wall, in which the lumen is partly filled with an inflammatory exudate of mucus and cells (middle). Chronically **...**

The mechanical origin of mucosal folding in healthy and asthmatic airways was first studied more than a decade ago using an analytical linear elastic model [60]. Linear stability analysis of bilayered tubes can explain to which extent the ratio between thicknesses and elastic moduli of the mucosal and submucosal layers dictates the morphology of surface buckling [37, 43]. A recent study confirmed the hypothesis that the wrinkling mode, the number of folds, is highly sensitive to the thickness ratio of these two layers, with more folds in the thin mucosal layer of the airway and fewer folds in the thick mucosal layer of the esophagus [36]. Beyond mucosal folding during airway narrowing in asthma [60] and surface wrinkling in the esophagus [36], morphological instabilities have been thoroughly investigated in the context of buckling during crypt formation in the intestinal wall [45], mucosal folding in the gastrointestinal tract [39], lumen narrowoing during chronic obstructive pulmonary disease [19], rule formation in cyclically loaded arteries [33], and buckling in the pharnyx during obstructive sleep apnoea [25]. In addition to these common diseases, understanding the morphogenesis and origin of shape has been identified as one of the key challenges in developmental biology [61, 62].

While the onset of mucosal buckling was initially attributed to an elevated external pressure, current hypotheses suggest that physiological and pathological thickening of the individual layers, as illustrated in Figure 1, might play a critical role in the generation of morphological instabilities [46, 44]. A recent model experiment of a cultured growing epithelial monolayer on a thin elastic substrate confirmed that growth of biological surfaces can, indeed, induce significant substrate buckling [45]. Another model experiment of a swelling polyacrylamide gel displayed similar buckling phenomena at the gel-water interface [10]. Only now, we are beginning to recognize the striking similarities between growing biological tissues and swelling gels, which display almost identical surface phenomena. These findings initiated a thorough mathematical analyzis of crease formation in growing soft matter [24] and of surface wrinkling in core-shell soft cylinders [4] using continuum models of isotropic volumetric growth. An excellent recent overview article illustrates various phenomena of morphological instabilities including applications beyond biological systems, such as soft lithography, metrology, and flexible electronics [38].

Conceptually, volumetric growth models can be classified into two categories, models of mechanically-induced biological growth and models of growth-induced mechanical instabilities. The first pioneering continuum model for mechancially-induced biological growth was proposed almost two decades ago [49]. Since then, it has been widely adopted to model muscle growth [56, 66], vascular growth [18, 57], skin growth [7, 53], and cardiac growth [14, 26] in response to elevated stresses or strains. Initially, these models were merely analytical [11, 12], for example focussing on volume growth of idealized rotationally symmetric cylindrical tubes [57, 59]. Now, computational growth models allow us to simulate growth in more realistic, fully three-dimensional, subject-specific geometries [29, 65]. While the first generation of growth models was rather phenomenological and isotropic in nature [1, 26], the current trend goes towards incorporating the underlying mechanobiology [64] and the microstructural origin of growth [13]. Depending on the particular tissue architecture, microstructurally-motivated models are therefore unsually transversely isotropic [47, 66], orthotropic [14], or generally anisotropic [41, 58]. Various different biomedical applications of mechanically-induced growth are summarized comprehensively in two recent overview articles [2, 42].

The first continuum model for growth-induced mechanical instabilities was proposed a decade ago to study tumor growth in confined geometries [1]. In the context of surface layer folding, the problem we are interested in here, recent studies have addressed circumferential buckling instabilites inside a growing cylindrical tube [36, 43]. Again, anisotropic models have been proposed to allow for a more sophisticated analysis using independent growth multipliers in the circumferential and radial directions [40, 44]. These models attribute buckling to the progression of growth, whereas other studies report buckling instabilities as a result of cyclic pressure changes [33] and mechanical unloading [6]. While mechanical instabilities have been investigated intensely in growing solids, the interplay between growing surfaces and morphological instabilities remains severely understudied [15]. Since growing surfaces may undergo finite deformations, analytical solutions are usually either unavailable or incredibly complex, and only a handful of mathematical models exist to model a growing surface by means of nonlinear plate [9] or shell [3] theories. Using shell theories of surface growth, a recent study has found that growth of biological membranes does not necessarily have to be detrimental: In some diseases such as mitral valve leakage, an adaptive surface growth of more than 30% has been reported as a compensatory mechanism to chronically reduce mitral valve insufficiency and deleterious backflow [48]. Here, we adopt the kinematical framework of growing biological surfaces; however, rather than allowing the surface to move freely, we attach it kinematically to the underlying solid, but equip it with its own potential energy.

The concept of surface energies is by no means new; in fact is has first been recognized more than two centuries ago, formalized through the famous Young Laplace equation [31, 63]. Relating the pressure difference across a fluid surface to surface tension and mean curvature, the Young Laplace equation is known to explain a variety of different phenomena ranging from minimal surfaces, soap bubbles, and capillary pressure to cardiovascular and respiratory physiology [55]. More than three decades ago, the familiar concept of scalar-valued surface tension was generalized to the tensorial notion of surface stress in the first continuum theory of elastic material surfaces [16]. Since then, the concept of material surfaces has gained increasing popularity in crystallography [34], phase transformation [35, 51], mineralogy [30], micro- and nanofabrication, and soft lithography; broadly speaking, whenever the surface displays distinct characteristic properties [17, 54]. Although propagations of internal surfaces [52] and morphological surface instabilities [35] are widely studied in metallurgy, material sciences, and mechanics today, the computational modeling of surfaces equipped with their own energies is still in its infancy. For fluids, a few finite element approaches exist to simulate droplets and free surfaces with scalar-valued surface tension [50]. For solids, however, a generic finite element approach towards elastic surfaces with tensorial surface stresses has only been proposed very recently [20]. Conceptually speaking, this approach treats the surface as a hyperelastic membrane of zero thickness glued to the underlying bulk. This concept is mathematically elegant and easily generalizable to anisotropic surfaces [21], thermomechanical surfaces [23] and surfaces with diffusion [5]. Here, we adopt this framework and model the growing surface as a growing shell, kinematically constrained to move with the solid body, but equipped with its own potential energy.

The manuscript is organized as follows. In Sections 2, 3, and 4, we introduce the kinematic equations, the balance equations, and the constitutive equations for finite growth. In all three sections, we first illustrate growth of the solid body, and then adopt the underlying concept to growth of the surface. In Section 5, we summarize the weak form of the governing equations and discretize it in space using a combination of volume elements for the volume terms and surface elements for the surface terms. In Section 6, we illustrate the features of our algorithm using the model problem of constrained growth in a cylindrical tube. We systematically compare a growing volume layer of finite thickness to a growing surface of zero thickness. In Section 7, we discuss the differences between volume and surface growth, before we close with a conclusion in Section 8.

Let _{0} ^{3} denote the material placement of a continuum body with material particles *X*_{0}. We denote the time domain by = [0, *T* ] _{+} and consider the quasi-static case, for which time provides a history parameter to order the sequence of events. We introduce the orientation preserving map : _{0} × *T* → ^{3} as the smooth motion of the material placement ** X** onto its spatial placement

$$\mathit{x}=\mathit{\phi}\phantom{\rule{thinmathspace}{0ex}}(\mathit{X},t),$$

(1)

and label the current placement of the body at time *t* as _{t} = (_{0}), see Figure 2. Here and in the sequel, the subscripts 0 and *t* designate material and spatial quantities. We further introduce the gradient, divergence, and material rate of any field {o} (** X**,

$$\nabla \{\u25cb\}={\partial}_{\mathit{X}}\{\u25cb\}{|}_{t}\text{Div}\phantom{\rule{thinmathspace}{0ex}}\{\u25cb\}=\nabla \{\u25cb\}:\mathit{I}\text{}{\mathrm{D}}_{t}\phantom{\rule{thinmathspace}{0ex}}\{\u25cb\}={\partial}_{t}\{\u25cb\}{|}_{\mathit{X}},$$

(2)

where ** I** is the material identity tensor, and {o}|

$$\mathit{F}=\nabla \mathit{\phi}\text{with}\mathit{F}={\mathit{F}}^{\mathrm{e}}\xb7{\mathit{F}}^{\mathrm{g}}.$$

(3)

Kinematics of material and spatial configurations _{0} and _{t} of material body with material and spatial surfaces _{0} and _{t} and multiplicative decomposition of the volume and surface deformation gradients *F* = *F*^{e} · **...**

In contrast to the deformation gradient ** F**, however, its individual contributions

The Jacobian *J* of the deformation gradient ** F** relates material volume elements d

$$J=\text{det}(\mathit{F})0\text{with}J={J}^{\mathrm{e}}{J}^{\mathrm{g}}.$$

(4)

We can explicitly express the elastic Jacobian *J*^{e} = det (*F*^{e}) and the growth Jacobian *J*^{g} = det (*F*^{g}), where we introduce the following explicit representation in terms of the *I* = 1,2,3 covariant base vectors * G_{I}* in

$$\text{det}({\mathit{F}}^{\mathrm{g}})=\frac{[{\mathit{F}}^{\mathrm{g}}\xb7{\mathit{G}}_{1}]\xb7[[{\mathit{F}}^{\mathrm{g}}\xb7{\mathit{G}}_{2}]\times [{\mathit{F}}^{\mathrm{g}}\xb7{\mathit{G}}_{3}]]}{{\mathit{G}}_{1}\xb7[{\mathit{G}}_{2}\times {\mathit{G}}_{3}]}.$$

(5)

Last, we introduce the right Cauchy Green tensor ** C** in the material configuration and its elastic counterpart

$$\mathit{C}={\mathit{F}}^{\mathrm{t}}\xb7\mathit{F}\text{and}{\mathit{C}}^{\mathrm{e}}={\mathit{F}}^{\text{et}}\xb7{\mathit{F}}^{\mathrm{e}},$$

(6)

which are related through the identity *C*^{e} = [*f*^{g}]^{t} · ** C** ·

$$\mathit{l}={\mathrm{D}}_{t}\mathit{F}\xb7\mathit{f}\text{with}{\mathit{f}}^{\mathrm{e}}\xb7\mathit{l}\xb7{\mathit{F}}^{\mathrm{e}}={\mathit{L}}^{\mathrm{e}}+{\mathit{L}}^{\mathrm{g}}$$

(7)

obeys an additive decomposition into the elastic velocity gradient *L*^{e} = *f*^{e} · D_{t}*F*^{e} and the growth velocity gradient *L*^{g} = D_{t}*F*^{g} · *f*^{g}.

Let _{0} = _{0} denote the surface of _{0}, assumed smooth, such that _{0} = _{0} _{0} is the closure of _{0} and _{0}. We label points on the surface of the reference configuration _{0} as ** = **** X**|

$$\mathit{x\u0302}=\mathit{\phi \u0302}(\mathit{X\u0302},t).$$

(8)

The above considerations implicitly assume kinematic slavery, i.e., points on the surface of the body remain on the surface at all times. In analogy to the volume, we formally introduce the surface gradient and surface divergence of any field {o} (**, ***t*) defined on the surface,

$$\nabla \u0302\{\u25cb\}=\nabla \{\u25cb\}\xb7\mathit{\xce}\text{}\widehat{\text{Div}}\phantom{\rule{thinmathspace}{0ex}}\{\u25cb\}=\nabla \{\u25cb\}:\mathit{\xce},$$

(9)

where ** Î** =

$$\mathit{F\u0302}=\nabla \u0302\mathit{\phi \u0302}=\mathit{F}\xb7\mathit{\xce}\text{with}\mathit{F\u0302}={\mathit{F\u0302}}^{\mathrm{e}}\xb7{\mathit{F\u0302}}^{\mathrm{g}}.$$

(10)

In contrast to the surface deformation gradient **, however, neither **^{e} nor ^{g} are necessarily restricted to result from the projection of the corresponding volume quantities *F*^{e} and *F*^{g}. Although the surface deformation gradient and its elastic and growth part are non-invertible, they possess inverses in the generalized sense, ** · **** = **** Î**,

$$\u0134=\widehat{\text{det}}\phantom{\rule{thinmathspace}{0ex}}(\mathit{F\u0302})=\Vert \text{cof}\phantom{\rule{thinmathspace}{0ex}}(\mathit{F})\xb7\mathit{N}\Vert >0\text{with}\u0134=\widehat{{J}^{\mathrm{e}}}\widehat{{J}^{\mathrm{g}}}.$$

(11)

Here cof (o) = det(o) (o)^{−t} denotes the cofactor of the second order tensor (o). We can then express the elastic surface Jacobian as $\widehat{{J}^{\mathrm{e}}}=\u0134/\widehat{{J}^{\mathrm{g}}}$ and the growth surface Jacobian as $\widehat{{J}^{\mathrm{g}}}=\widehat{\text{det}}({\mathit{F\u0302}}^{\mathrm{g}})$ in terms of the covariant base vectors *A*_{α} in *T*_{0} [17],

$$\widehat{\text{det}}\phantom{\rule{thinmathspace}{0ex}}({\mathit{F\u0302}}^{\mathrm{g}})=\frac{|[{\mathit{F\u0302}}^{\mathrm{g}}\xb7{\mathit{A}}_{1}]\times [{\mathit{F\u0302}}^{\mathrm{g}}\xb7{\mathit{A}}_{2}]|}{|{\mathit{A}}_{1}\times {\mathit{A}}_{2}|}.$$

(12)

Finally, we introduce the right Cauchy Green surface tensor ** Ĉ** in the material configuration and its elastic counterpart

$$\mathit{\u0108}={\mathit{F\u0302}}^{\mathrm{t}}\xb7\mathit{F\u0302}\text{and}{\mathit{\u0108}}^{\mathrm{e}}={\mathit{F\u0302}}^{\text{et}}\xb7{\mathit{F\u0302}}^{\mathrm{e}}$$

(13)

which are related through the identity, *Ĉ*^{e} = ^{gt} · ** Ĉ** ·

$$\mathit{l\u0302}={\mathrm{D}}_{t}\mathit{F\u0302}\xb7\mathit{f\u0302}\text{with}{\mathit{f\u0302}}^{\mathrm{e}}\xb7\mathit{l\u0302}\xb7{\mathit{F\u0302}}^{\mathrm{e}}={\mathit{L\u0302}}^{\mathrm{e}}+{\mathit{L\u0302}}^{\mathrm{g}},$$

(14)

which obeys an additive decomposition into the elastic surface velocity gradient ^{e} = ^{e} · D_{t}^{e} and the growth surface velocity gradient ^{g} = D_{t}^{g} · ^{g}.

**Remark 1.**
*The material and spatial second order surface unit tensors Î* =

$$\mathit{F\u0302}=\mathit{F}\xb7\mathit{\xce}\text{}\mathit{\xce}=\mathit{I}-\mathit{N}\otimes \mathit{N}$$

(15)

$$\mathit{f\u0302}=\mathit{f}\xb7\mathit{\xee}\text{}\mathit{\xee}=\mathit{i}-\mathit{n}\otimes \mathit{n},$$

*where N and n* =

$$\mathit{f\u0302}\xb7\mathit{F\u0302}=\mathit{\xce}\text{}\mathit{F\u0302}=\mathit{U}\xb7\mathit{\Sigma}\xb7{\mathit{V}}^{\mathrm{t}}\text{}\mathit{f\u0302}=\mathit{V}\xb7{[{\mathit{\Sigma}}^{+}]}^{-1}\xb7{\mathit{U}}^{\mathrm{t}}$$

(16)

$${\mathit{f\u0302}}^{\mathrm{t}}\xb7{\mathit{F\u0302}}^{\mathrm{t}}=\mathit{\xee}\text{}{\mathit{F\u0302}}^{\mathrm{t}}=\mathit{U}\xb7\mathit{\sigma}\xb7{\mathit{V}}^{\mathrm{t}}\text{}{\mathit{f\u0302}}^{\mathrm{t}}=\mathit{V}\xb7{[{\mathit{\sigma}}^{+}]}^{-1}\xb7{\mathit{U}}^{\mathrm{t}},$$

*where the diagonal entries of*
**Σ**
*and*
**σ**
*correspond to the singular values of and*
^{t}, *the columns of U and V are the left- and right-singular vectors associated with these singular values, and*

The balance of mass for growing continua balances the rate of change of the material volume density ρ_{0} with the volume mass source _{0} and the volume mass flux ** R** [27, 28]. For simplicity, we assume the mass flux to be negligibly small,

$${\mathrm{D}}_{t}{\rho}_{0}={\mathcal{R}}_{0}\text{or}{\rho}_{0}(t)={\rho}_{0}({t}_{0})+{\displaystyle {\int}_{{t}_{0}}^{t}}{\mathcal{R}}_{0}(\tau )\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\tau .$$

(17)

By using the volume Jacobians introduced in equation (4), we immediately obtain the following relations between the volume densities in the material, spatial, and intermediate configurations ρ_{0}, ρ_{t}, and ρ_{g},

$${\rho}_{0}=J{\rho}_{t}={J}^{\mathrm{g}}{\rho}_{\mathrm{g}}\text{such that}{\mathrm{D}}_{t}{\rho}_{0}={J}^{\mathrm{g}}{\mathrm{D}}_{t}{\mathit{F}}^{\mathrm{g}}:{\mathit{f}}^{\mathrm{g}}{\rho}_{\mathrm{g}}+{J}^{\mathrm{g}}{\mathrm{D}}_{t}{\rho}_{\mathrm{g}}={J}^{\mathrm{g}}\phantom{\rule{thinmathspace}{0ex}}[{\rho}_{\mathrm{g}}\text{tr}{\mathit{L}}^{\mathrm{g}}+{\mathrm{D}}_{t}{\rho}_{\mathrm{g}}].$$

(18)

In the transformation above, we have utilized the time derivative of *J*^{g} as D_{t}*J*^{g} = _{F}^{g}*J*^{g} : D_{t}*F*^{g} = *J*^{g}[*f*^{g}]^{t} : D_{t}*F*^{g} = *J*^{g} tr *L*^{g}. The combination of equations (17.1) and (18.2) renders the following explicit representation of the mass source _{0}

$${\mathcal{R}}_{0}={J}^{\mathrm{g}}\phantom{\rule{thinmathspace}{0ex}}[{\rho}_{\mathrm{g}}\text{tr}{\mathit{L}}^{\mathrm{g}}+{\mathrm{D}}_{t}{\rho}_{\mathrm{g}}].$$

(19)

The balance of linear momentum balances the rate of change of the linear momentum, which we assume to vanish here for the quasi-static case, with the divergence of the volume Piola stresses ** P** and the volume forces

$$\mathbf{0}=\text{Div}\mathit{P}+\mathit{b}.$$

(20)

The Piola stresses ** P** and volume forces

$${\int}_{{\mathcal{B}}_{0}}}\nabla \delta \mathit{\phi}:\mathit{P}\mathrm{d}{V}_{0}={\displaystyle {\int}_{{\mathcal{B}}_{0}}}\delta \mathit{\phi}\xb7\mathit{b}\mathrm{d}{V}_{0}+{\displaystyle {\int}_{{\mathcal{S}}_{0}}}\delta \mathit{\phi}\xb7\mathit{P}\xb7\mathit{N}\mathrm{d}{A}_{0}.$$

(21)

Here, the lefthand side characterizes the internal forces, and the righthand side terms represent the external volume forces and the external surface forces. Alternatively, we could rephrase the weak form in terms of the Piola Kirchhoff stress ** S** or the Cauchy stress

$$\mathit{S}=\mathit{f}\xb7\mathit{P}\text{}\mathit{\sigma}=j\mathit{P}\xb7{\mathit{F}}^{\mathrm{t}},$$

(22)

such that the internal force term of the weak form (21) takes the following alternative representations,

$${\int}_{{\mathcal{B}}_{0}}}\nabla \delta \mathit{\phi}:\mathit{P}\mathrm{d}{V}_{0}={\displaystyle {\int}_{{\mathcal{B}}_{0}}}[{\mathit{F}}^{\mathrm{t}}\xb7\nabla \delta \mathit{\phi}]:\mathit{S}\mathrm{d}{V}_{0}={\displaystyle {\int}_{{\mathcal{B}}_{\mathrm{t}}}}[\nabla \delta \mathit{\phi}\xb7\mathit{f}]:\mathit{\sigma}\mathrm{d}{V}_{\mathrm{t}}.$$

(23)

The balance of mass for growing surfaces balances the rate of change of the material area density _{0} with the area mass source _{0} and the area mass flux **, which we again assume to vanish in the sequel, **** = ****0**, such that

$${\mathrm{D}}_{t}{\mathrm{\rho \u0302}}_{0}={\mathrm{\mathcal{R}\u0302}}_{0}\text{or}{\mathrm{\rho \u0302}}_{0}(t)={\mathrm{\rho \u0302}}_{0}({t}_{0})+{\displaystyle {\int}_{{t}_{0}}^{t}}{\mathrm{\mathcal{R}\u0302}}_{0}(\tau )\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\tau .$$

(24)

We can then relate the surface densities in the material, spatial, and intermediate configurations _{0}, _{t}, and _{g} with the help of the surface Jacobians (11)

$${\mathrm{\rho \u0302}}_{0}=\u0134{\mathrm{\rho \u0302}}_{t}=\widehat{{J}^{\mathrm{g}}}{\mathrm{\rho \u0302}}_{\mathrm{g}}\text{such that}{\mathrm{D}}_{t}{\mathrm{\rho \u0302}}_{0}=\widehat{{J}^{\mathrm{g}}}{\mathrm{D}}_{t}{\mathit{F\u0302}}^{\mathrm{g}}\xb7{\mathit{f\u0302}}^{\mathrm{g}}{\mathrm{\rho \u0302}}_{\mathrm{g}}+{J}^{\mathrm{g}}{\mathrm{D}}_{t}{\rho}_{\mathrm{g}}=\widehat{{J}^{\mathrm{g}}}\phantom{\rule{thinmathspace}{0ex}}[{\mathrm{\rho \u0302}}_{\mathrm{g}}{\mathit{L\u0302}}^{\mathrm{g}}+{\mathrm{D}}_{t}{\mathrm{\rho \u0302}}_{\mathrm{g}}],$$

(25)

where ${\mathrm{D}}_{t}\widehat{{J}^{\mathrm{g}}}={\partial}_{{\mathit{F\u0302}}^{\mathrm{g}}}\widehat{{J}^{\mathrm{g}}}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}{\mathrm{D}}_{t}{\mathit{F\u0302}}^{\mathrm{g}}=\widehat{{J}^{\mathrm{g}}}{[{\mathit{f\u0302}}^{\mathrm{g}}]}^{\mathrm{t}}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}{\mathrm{D}}_{t}{\mathit{F\u0302}}^{\mathrm{g}}=\widehat{{J}^{\mathrm{g}}}\text{tr}{\mathit{L\u0302}}^{\mathrm{g}}$. The mass source then follows directly from the combination of equations (24.1) and (25.2) as

$${\mathrm{\mathcal{R}\u0302}}_{0}=\widehat{{J}^{\mathrm{g}}}\phantom{\rule{thinmathspace}{0ex}}[{\mathrm{\rho \u0302}}_{\mathrm{g}}\text{tr}{\mathit{L\u0302}}^{\mathrm{g}}+{\mathrm{D}}_{t}{\mathrm{\rho \u0302}}_{\mathrm{g}}].$$

(26)

The balance of linear momentum balances the rate of change of the linear surface momentum, which we assume to be negligible in the quasi-static case, with the divergence of the surface stresses ** and the surface tractions **** − **** P** ·

$$\mathbf{0}=\widehat{\text{Div}}\phantom{\rule{thinmathspace}{0ex}}\mathit{P\u0302}+[\mathit{b\u0302}-\mathit{P}\xb7\mathit{N}].$$

(27)

Here, the surface stresses ** and the surface tractions **** − **** P** ·

$${\int}_{{\mathcal{S}}_{0}}}\nabla \u0302\delta \mathit{\phi}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\mathit{P\u0302}\mathrm{d}{A}_{0}={\displaystyle {\int}_{{\mathcal{S}}_{0}}}\delta \mathit{\phi}\xb7\mathit{b\u0302}\mathrm{d}{A}_{0}-{\displaystyle {\int}_{{\mathcal{S}}_{0}}}\delta \mathit{\phi}\xb7\mathit{P}\xb7\mathit{N}\mathrm{d}{A}_{0}+{\displaystyle {\int}_{{\mathcal{C}}_{0}}}\delta \mathit{\phi}\xb7\mathit{P\u0302}\xb7\mathit{N}\mathrm{d}{\mathit{L}}_{0}.$$

(28)

The lefthand side characterizes the internal surface forces, and the righthand side terms represent the prescribed external surface forces, the projected external surface forces imposed by the underlying volume, and the external line forces along the boundary curve to the surface _{0}. Again, we can introduce the surface Piola Kirchhoff stress ** Ŝ** and the surface Cauchy stress using the surface Piola transforms,

$$\mathit{\u015c}=\mathit{f\u0302}\xb7\mathit{P\u0302}\text{}\mathit{\sigma \u0302}=\u0135\mathit{P\u0302}\xb7{\mathit{F\u0302}}^{\mathrm{t}},$$

(29)

to rephrase the internal forces of the weak form (28) through the following expressions,

$${\int}_{{\mathcal{S}}_{0}}}\nabla \u0302\delta \mathit{\phi}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\mathit{P\u0302}\mathrm{d}{A}_{0}={\displaystyle {\int}_{{\mathcal{S}}_{0}}}[{\mathit{F\u0302}}^{\mathrm{t}}\xb7\nabla \u0302\delta \mathit{\phi}]\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\mathit{\u015c}\mathrm{d}{A}_{0}={\displaystyle {\int}_{{\mathcal{S}}_{\mathrm{t}}}}[\nabla \u0302\delta \mathit{\phi}\xb7\mathit{f\u0302}]\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\mathit{\sigma \u0302}\mathrm{d}{A}_{\mathrm{t}}.$$

(30)

To characterize volume growth in a hyperelastic medium, we introduce a Helmholtz free volume energy ψ, parameterized in terms of the elastic contribution to the volume deformation gradient *F*^{e}, or, alternatively, in terms of the volume deformation gradient ** F** and the volume growth tensor

$$\psi =\psi (\mathit{F},{\mathit{F}}^{\mathrm{g}})=\psi ({\mathit{F}}^{\mathrm{e}}).$$

(31)

The corresponding volume Piola stress ** P** follows from thermodynamic considerations as stress measure conjugate to the deformation gradient

$$\mathit{P}=\frac{\partial \psi}{\partial \mathit{F}}=\frac{\partial \psi}{\partial {\mathit{F}}^{\mathrm{e}}}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\frac{\partial {\mathit{F}}^{\mathrm{e}}}{\partial \mathit{F}}={\mathit{P}}^{\mathrm{e}}\xb7{[{\mathit{f}}^{\mathrm{g}}]}^{\mathrm{t}}\text{with}{\mathit{P}}^{\mathrm{e}}=\frac{\partial \psi}{\partial {\mathit{F}}^{\mathrm{e}}},$$

(32)

where *P*^{e} denotes the classic elastic volume Piola stress. The total derivative of the Piola stress ** P** with respect to the deformation gradient

$$\mathbb{A}=\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{F}}=\frac{\mathrm{d}\mathit{P}}{\mathrm{d}{\mathit{F}}^{\mathrm{e}}}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\frac{\mathrm{d}{\mathit{F}}^{\mathrm{e}}}{\mathrm{d}\mathit{F}}={\mathbb{A}}^{\mathrm{e}}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}[{[{\mathit{f}}^{\mathrm{g}}]}^{\mathrm{t}}\otimes \u0305{[{\mathit{f}}^{\mathrm{g}}]}^{\mathrm{t}}]\text{with}{\mathbb{A}}^{\mathrm{e}}=\frac{\mathrm{d}{\mathit{P}}^{\mathrm{e}}}{\mathrm{d}{\mathit{F}}^{\mathrm{e}}},$$

(33)

where ^{e} denotes the classic elastic volume tangent moduli. To characterize the growth process, we need to specify a functional form for the volume growth tensor *F*^{g} and its evolution in time. For the simplest case of isotropic volume growth, the volume growth tensor *F*^{g} takes the following functional form, parameterized in terms of a single scalar-valued growth multiplier ,

$${\mathit{F}}^{\mathrm{g}}=\vartheta \phantom{\rule{thinmathspace}{0ex}}\mathit{I},$$

(34)

such that = 1 in the initial ungrown state, > 1 characterizes volume growth, and < 1 characterizes volume shrinkage. We propose the following simple exponential evolution equation for the growth multiplier, which we can integrate explicitly in time,

$${\mathrm{D}}_{t}\vartheta =[{\vartheta}^{\text{max}}-{\vartheta}_{0}]\phantom{\rule{thinmathspace}{0ex}}[\text{exp}(-t/\tau )]/\tau \text{thus}\vartheta ={\vartheta}_{0}+[{\vartheta}^{\text{max}}-{\vartheta}_{0}]\phantom{\rule{thinmathspace}{0ex}}[1-\text{exp}(-t/\tau )].$$

(35)

Here, τ characterizes the velocity of volume growth, _{0} = 1.0 is the initial growth value, and ^{max} limits the maximum amount of growth towards which the growth multiplier converges with progressing time *t* → ∞. Figure 3 illustrates the temporal evolution of the growth multiplier , with an initial value _{0} = 1.0, for varying growth velocities τ at fixed ^{max} = 4.0, left, and varying maximum growth ^{max} at fixed τ = 10, right.

Temporal evolution of growth multiplier for varying growth velocities τ at fixed ^{max}, left, and varying maximum growth ^{max} at fixed τ, right.

**Remark 2 (Isotropic volume growth).**
*The special case of isotropic volume growth, F*

$${\mathit{f}}^{\mathrm{g}}=\frac{1}{\vartheta}\mathit{I}\mathit{\text{such that}}\mathit{F}=\vartheta \phantom{\rule{thinmathspace}{0ex}}{\mathit{F}}^{\mathrm{e}}\mathit{\text{and}}\mathit{P}=\frac{1}{\vartheta}{\mathit{P}}^{\mathrm{e}}\mathit{\text{and}}\mathbb{A}=\frac{1}{{\vartheta}^{2}}{\mathbb{A}}^{\mathrm{e}}.$$

(36)

*Second, the growth tensor F*

$${\mathrm{D}}_{t}{\mathit{F}}^{\mathrm{g}}={\mathrm{D}}_{t}\vartheta \phantom{\rule{thinmathspace}{0ex}}\mathit{I}\mathit{\text{such that}}{\mathit{L}}^{\mathrm{g}}=\frac{{\mathrm{D}}_{t}\vartheta}{\vartheta}\mathit{I}.$$

(37)

*Third, volume changes upon growth J*^{g}
*are simply equivalent to the third power of the growth multiplier* ,

$${J}^{\mathrm{g}}={\vartheta}^{3}\mathit{\text{such that}}{\rho}_{0}={\vartheta}^{3}{\rho}_{\mathrm{g}}.$$

(38)

*Overall, the assumption of isotropic growth significantly simplifies the general equations of volumetric growth*.

**Remark 3 (Growth at constant volume density).**
*A common additional assumption is that the newly grown material has the same volume density as the original substrate* [42], D_{t}ρ_{g} = 0. *This implies that the volume density source* _{0}
*required to maintain a constant volume density in the intermediate configuration takes the following form*,

$${\rho}_{\mathrm{g}}=\mathit{\text{const}}.\mathit{\text{such that}}{\mathcal{R}}_{0}={\rho}_{0}\text{tr}{\mathit{L}}^{\mathrm{g}}={J}^{\mathrm{g}}{\rho}_{\mathrm{g}}{\mathrm{D}}_{t}{\mathit{F}}^{\mathrm{g}}\xb7{\mathit{f}}^{\mathrm{g}}=3\phantom{\rule{thinmathspace}{0ex}}{\rho}_{\mathrm{g}}\phantom{\rule{thinmathspace}{0ex}}\vartheta {\mathrm{D}}_{t}\vartheta .$$

(39)

*For this special case, we do not need to discretize the balance of mass* (17) *explicitly, but can simply evaluate overall mass changes in a straightforward post-processing step*.

To characterize surface growth, we equip the surface with its own Helmholtz free surface energy , which we parameterize in terms of the elastic contribution to the surface deformation gradient ^{e}, or, alternatively, in terms of the surface deformation gradient ** and the surface growth tensor **^{g},

$$\mathrm{\psi \u0302}=\mathrm{\psi \u0302}(\mathit{F\u0302},{\mathit{F\u0302}}^{\mathrm{g}})=\mathrm{\psi \u0302}({\mathit{F\u0302}}^{\mathrm{e}}).$$

(40)

For simplicity, we assume the surface to behave isotropically [54], and do not account for an explicit dependence on the deformed surface plane normal ** n**. Similar to the volume, we introduce the surface Piola stress

$$\mathit{P\u0302}=\frac{\partial \mathrm{\psi \u0302}}{\partial \mathit{F\u0302}}=\frac{\partial \mathrm{\psi \u0302}}{\partial {\mathit{F\u0302}}^{\mathrm{e}}}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\frac{\partial {\mathit{F\u0302}}^{\mathrm{e}}}{\partial \mathit{F\u0302}}={\mathit{P\u0302}}^{\mathrm{e}}\xb7{[{\mathit{f\u0302}}^{\mathrm{g}}]}^{\mathrm{t}}\text{with}{\mathit{P\u0302}}^{\mathrm{e}}=\frac{\partial \psi}{\partial {\mathit{F\u0302}}^{\mathrm{e}}},$$

(41)

where ^{e} denotes the elastic surface Piola stress. Additionally, we introduce the fourth order tensor of tangent surface moduli as the total derivative of the Piola stress ** with respect to the deformation gradient ****,
**

$$\mathrm{\mathbb{A}\u0302}=\frac{\mathrm{d}\mathit{P\u0302}}{\mathrm{d}\mathit{F\u0302}}=\frac{\mathrm{d}\mathit{P\u0302}}{\mathrm{d}{\mathit{F\u0302}}^{\mathrm{e}}}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\frac{\mathrm{d}{\mathit{F\u0302}}^{\mathrm{e}}}{\mathrm{d}\mathit{F\u0302}}=\widehat{{\mathbb{A}}^{\mathrm{e}}}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}[{[{\mathit{f\u0302}}^{\mathrm{g}}]}^{\mathrm{t}}\otimes \u0305{[{\mathit{f\u0302}}^{\mathrm{g}}]}^{\mathrm{t}}]\text{with}{\mathrm{\mathbb{A}\u0302}}^{\mathrm{e}}=\frac{\mathrm{d}{\mathit{P\u0302}}^{\mathrm{e}}}{\mathrm{d}{\mathit{F\u0302}}^{\mathrm{e}}},$$

(42)

such that ^{e} denotes the elastic surface tangent moduli. Again, we consider the simplest possible case of growth, isotropic surface growth, for which the surface growth tensor ^{g} can be parameterized in terms of a single scalar-valued growth multiplier ,

$${\mathit{F\u0302}}^{\mathrm{g}}=\mathrm{\vartheta \u0302}\mathit{\xce}.$$

(43)

Similar to the case of volume growth, = 1 characterizes the initial ungrown state, > 1 corresponds to surface growth, and < 1 to surface shrinkage. Again, we suggest a simple exponential evolution equation for the growth multiplier, and integrate it explicitly in time,

$${\mathrm{D}}_{t}\mathrm{\vartheta \u0302}=[{\mathrm{\vartheta \u0302}}^{\text{max}}-{\mathrm{\vartheta \u0302}}_{0}]\phantom{\rule{thinmathspace}{0ex}}[\text{exp}(-t/\mathrm{\tau \u0302})]/\mathrm{\tau \u0302}\text{thus}\mathrm{\vartheta \u0302}={\mathrm{\vartheta \u0302}}_{0}+[{\mathrm{\vartheta \u0302}}^{\text{max}}-{\mathrm{\vartheta \u0302}}_{0}]\phantom{\rule{thinmathspace}{0ex}}[1-\text{exp}(-t/\mathrm{\tau \u0302})].$$

(44)

Here, characterizes the velocity of surface growth, _{0} = 1.0 is the initial growth value, and ^{max} limits the maximum amount of surface growth towards which the growth multiplier converges with progressing time *t*. Figure 3 illustrates the temporal evolution of the growth multiplier , with an initial value _{0} = 1.0, for varying growth velocities at fixed ^{max} = 4.0, left, and varying maximum growth ^{max} at fixed = 10, right.

**Remark 4 (Isotropic surface growth).**
*The special case of isotropic surface growth, *^{g} = **Î**, has three immediate consequences. First, we can explicitly invert the surface growth tensor ^{g}
*using the surface unit tensor Î. This allows us to express the surface deformation gradient F, the surface Piola stress P, and the surface tangent moduli*

$${\mathit{f\u0302}}^{\mathrm{g}}=\frac{1}{\mathrm{\vartheta \u0302}}\mathit{\xce}\mathit{\text{such that}}\mathit{F\u0302}=\mathrm{\vartheta \u0302}{\mathit{F\u0302}}^{\mathrm{e}}\mathit{\text{and}}\mathit{P\u0302}=\frac{1}{\mathrm{\vartheta \u0302}}{\mathit{P\u0302}}^{\mathrm{e}}\mathit{\text{and}}\mathrm{\mathbb{A}\u0302}=\frac{1}{{\mathrm{\vartheta \u0302}}^{2}}{\mathrm{\mathbb{A}\u0302}}^{\mathrm{e}}.$$

(45)

*Second, the surface growth tensor F*

$${\mathrm{D}}_{t}{\mathit{F\u0302}}^{\mathrm{g}}={\mathrm{D}}_{t}\mathrm{\vartheta \u0302}\mathit{\xce}\mathit{\text{such that}}{\mathit{L\u0302}}^{\mathrm{g}}=\frac{{\mathrm{D}}_{t}\mathrm{\vartheta \u0302}}{\mathrm{\vartheta \u0302}}\mathit{\xce}.$$

(46)

*Third, surface area changes upon growth*
$\widehat{{J}^{\mathrm{g}}}$
*are simply equivalent to the growth multiplier* *squared*,

$$\widehat{{J}^{\mathrm{g}}}={\mathrm{\vartheta \u0302}}^{2}\mathit{\text{such that}}{\mathrm{\rho \u0302}}_{0}={\mathrm{\vartheta \u0302}}^{2}{\mathrm{\rho \u0302}}_{\mathrm{g}}.$$

(47)

*As the above considerations indicate, the assumption of isotropic growth simplifies the governing equations of surface growth*.

**Remark 5 (Growth at constant surface density).**
*If we additionally assume that the newly grown surface has the same surface density as the original surface*, D_{t}_{g} = 0, *we can express the surface density source* _{0}
*in the following form*,

$${\mathrm{\rho \u0302}}_{\mathrm{g}}=\mathit{\text{const}}.\mathit{\text{such that}}{\mathrm{\mathcal{R}\u0302}}_{0}={\mathrm{\rho \u0302}}_{0}\text{tr}{\mathit{L\u0302}}^{\mathrm{g}}=\widehat{{J}^{\mathrm{g}}}{\mathrm{\rho \u0302}}_{\mathrm{g}}{\mathrm{D}}_{t}{\mathit{F\u0302}}^{\mathrm{g}}\xb7{\mathit{f\u0302}}^{\mathrm{g}}=2{\mathrm{\rho \u0302}}_{\mathrm{g}}\mathrm{\vartheta \u0302}{\mathrm{D}}_{t}\mathrm{\vartheta \u0302}.$$

(48)

*Similar to the case of constant volume density, we do not need to discretize the balance of mass* (24) *explicitly, but can simply evaluate overall mass changes in a straightforward post-processing step*.

The governing equations for finite volume and surface growth are complex and highly nonlinear. In this section, we illustrate their computational solution within an incremental iterative nonlinear finite element framework. To characterize the state of the growth process at each instant in time, we introduce the volume and surface growth multipliers ^{g} and ^{g} as internal variables, and solve their evolution equations (35) and (44) locally at the integration point level. We begin by summarizing the residual of the weak forms of the balance of linear momentum in the volume (21) and on the surface (28),

$${\int}_{{\mathcal{B}}_{0}}}\nabla \delta \mathit{\phi}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\mathit{P}\mathrm{d}{V}_{0}-{\displaystyle {\int}_{{\mathcal{B}}_{0}}}\delta \mathit{\phi}\xb7\mathit{b}\mathrm{d}{V}_{0}+{\displaystyle {\int}_{{\mathcal{S}}_{0}}}\nabla \u0302\delta \mathit{\phi}\phantom{\rule{thinmathspace}{0ex}}:\phantom{\rule{thinmathspace}{0ex}}\mathit{P\u0302}\mathrm{d}{A}_{0}-{\displaystyle {\int}_{{\mathcal{S}}_{0}}}\delta \mathit{\phi}\xb7\mathit{b\u0302}\mathrm{d}{A}_{0}-{\displaystyle {\int}_{{\mathcal{C}}_{0}}}\delta \mathit{\phi}\xb7\mathit{P\u0302}\xb7\mathit{N}\mathrm{d}{\mathit{L}}_{0}\doteq 0.$$

(49)

To discretize the weak form (49) in space, we partition the domain of interest ${\mathcal{B}}_{0}={\displaystyle {\cup}_{\mathrm{e}=1}^{{\mathrm{n}}_{\text{be}}}}{\mathcal{B}}_{0}^{\mathrm{e}}$ into n_{be} finite volume elements ${\mathcal{B}}_{0}^{\mathrm{e}}$, and the surface of interest ${\mathcal{S}}_{0}={\displaystyle {\cup}_{\mathrm{e}=1}^{{\mathrm{n}}_{\text{se}}}}{\mathcal{S}}_{0}^{\mathrm{e}}$ into n_{se} finite surface elements ${\mathcal{S}}_{0}^{\mathrm{e}}$. Although it is in principle possible to select arbitrary surface elements with independent nodes, it proofs efficient to discretize the surface in consistency with the bulk. This implies that each surface element shares its nodes with its corresponding volume elements as illustrated in Figure 4. The surface element itself than acts like a shell element, which moves in consistency with the bulk; however, it is equipped with its own independent free energy function. In the following subsections, we demonstrate the finite element discretizations of the volume and surface terms of equation (49), to obtain the discrete residual

$${\mathbf{R}}_{\mathit{I}}=\underset{\mathrm{e}=1}{\overset{{\mathrm{n}}_{\text{el}}}{\mathbf{A}}}\phantom{\rule{thinmathspace}{0ex}}{\displaystyle {\int}_{{\mathcal{B}}_{0}^{\mathrm{e}}}}\nabla {N}^{i}\xb7\mathit{P}\mathrm{d}{V}_{\mathrm{e}}-{\displaystyle {\int}_{{\mathcal{B}}_{0}^{\mathrm{e}}}}{N}^{i}\mathit{b}\mathrm{d}{V}_{\mathrm{e}}+{\displaystyle {\int}_{{\mathcal{S}}_{0}^{\mathrm{e}}}}\nabla \u0302{N}^{i}\xb7\mathit{P\u0302}\mathrm{d}{A}_{\mathrm{e}}-{\displaystyle {\int}_{{\mathcal{S}}_{0}^{\mathrm{e}}}}{\mathrm{N\u0302}}^{i}\mathit{b\u0302}\mathrm{d}{A}_{\mathrm{e}}-{\displaystyle {\int}_{{\mathcal{C}}_{0}^{\mathrm{e}}}}{\mathrm{N\u0302}}^{i}\mathit{P\u0302}\xb7\mathit{N}\mathrm{d}{\mathit{L}}_{\mathrm{e}}\doteq \mathbf{0}.$$

(50)

Discretization of material body with finite volume elements ${\mathcal{B}}_{0}^{\mathrm{e}}$ and ${\mathcal{B}}_{t}^{\mathrm{e}}$ and of its surface with finite surface elements ${\mathcal{B}}_{0}^{\mathrm{e}}$ and ${\mathcal{S}}_{t}^{\mathrm{e}}$. Surface elements share their nodes with the corresponding volume elements, however, they **...**

Herein, the operator **A** symbolizes the assembly of all element residuals at the *i* = 1, ‥, n_{vn} volume element nodes to the global residual at the global node points *I* = 1, ‥, n_{gn}. To solve the above equation, we suggest an incremental iterative Newton algorithm based on the consistent linearization of the residual **R**_{I} with respect to the nodal vector of unknowns _{J}. This linearization introduces the global stiffness matrix at all global nodes *I*, *J* = 1, ‥, n_{ng},

$${\mathbf{K}}_{\mathit{\text{IJ}}}=\frac{\mathrm{d}{\mathbf{R}}_{I}}{\mathrm{d}{\mathit{\phi}}_{J}}=\underset{\mathrm{e}=1}{\overset{{\mathrm{n}}_{\text{el}}}{\mathbf{A}}}{\displaystyle {\int}_{{\mathcal{B}}_{0}^{\mathrm{e}}}}[\mathit{I}\xb7\nabla {N}^{i}]:\mathbb{A}\xb7\nabla {N}^{j}\mathrm{d}{V}_{\mathrm{e}}+{\displaystyle {\int}_{{\mathcal{S}}_{0}^{\mathrm{e}}}}[\mathit{\xce}\xb7\nabla \u0302{N}^{i}]:\mathrm{\mathbb{A}\u0302}\xb7\nabla \u0302{N}^{j}\mathrm{d}{A}_{\mathrm{e}}.$$

(51)

For each global Newton iteration step, we iteratively update the current deformation state ${\mathit{\phi}}_{J}\leftarrow {\mathit{\phi}}_{J}-{\mathbf{K}}_{\mathit{\text{IJ}}}^{-1}\xb7{\mathbf{R}}_{I}$ until we achieve algorithmic convergence. Upon convergence, we store the volume and surface growth multipliers ^{g} and ^{g} at the integration points of the corresponding volume and surface elements.

To approximate the test functions δ, trial functions , and nodal coordinates ** X** in the volume, we apply an isoparametric Bubnov-Galerkin type finite element interpolation,

$$\delta \mathit{\phi}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{vn}}}}{N}^{i}\phantom{\rule{thinmathspace}{0ex}}\delta {\mathit{\phi}}_{i}\text{and}\mathit{\phi}={\displaystyle {\sum}_{j=1}^{{\mathrm{n}}_{\text{vn}}}}{N}^{j}\phantom{\rule{thinmathspace}{0ex}}{\mathit{\phi}}_{j}\text{and}\mathit{X}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{vn}}}}{N}^{i}\phantom{\rule{thinmathspace}{0ex}}{\mathit{X}}_{i}$$

(52)

where *N ^{i}* and

$$\delta \mathit{F}=\nabla \delta \mathit{\phi}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{vn}}}}\delta {\mathit{\phi}}_{i}\otimes \nabla {N}^{i}\text{and}\mathit{F}=\nabla \mathit{\phi}={\displaystyle {\sum}_{j=1}^{{\mathrm{n}}_{\text{vn}}}}{\mathit{\phi}}_{j}\otimes \nabla {N}^{j}.$$

(53)

To determine the gradients of the shape functions *N ^{i}*, we use the isoparameteric concept for the approximation of the nodal coordinates

$$\nabla {N}^{i}=\frac{\mathrm{d}{N}^{i}(\mathit{\xi})}{\mathrm{d}\mathit{X}}=\frac{\mathrm{d}{N}^{i}(\mathit{\xi})}{\mathrm{d}\mathit{\xi}}\xb7\frac{\mathrm{d}\mathit{\xi}}{\mathrm{d}\mathit{X}}={\mathbf{J}}^{-\mathrm{t}}\xb7\frac{\mathrm{d}{N}^{i}(\mathit{\xi})}{\mathrm{d}\mathit{\xi}}\text{with}\mathbf{J}=\frac{\mathrm{d}\mathit{X}}{\mathrm{d}\mathit{\xi}}\text{from}\mathit{X}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{vn}}}}{N}^{i}(\mathit{\xi})\phantom{\rule{thinmathspace}{0ex}}{\mathit{X}}_{i}.$$

(54)

To approximate the test functions δ, trial functions , and nodal coordinates ** X** on the surface, we apply a Bubnov-Galerkin type finite element interpolation,

$$\delta \mathit{\phi}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{sn}}}}{\mathrm{N\u0302}}^{i}\delta {\mathit{\phi}}_{i}\text{and}\mathit{\phi}={\displaystyle {\sum}_{j=1}^{{\mathrm{n}}_{\text{sn}}}}{\mathrm{N\u0302}}^{j}{\mathit{\phi}}_{j}\text{and}\mathit{X}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{sn}}}}{\mathrm{N\u0302}}^{i}{\mathit{X}}_{i}$$

(55)

where * ^{i}* and

$$\delta \mathit{F}=\nabla \u0302\delta \mathit{\phi}=\delta {\mathit{a}}_{\alpha}\otimes {\mathit{A}}^{\alpha}\text{and}\mathit{F}=\nabla \u0302\mathit{\phi}={\mathit{a}}_{\alpha}\otimes {\mathit{A}}^{\alpha}.$$

(56)

Accordingly, we first calculate the covariant spatial and material base vectors δ*a*_{α}, *a*_{α}, and *A*_{α}, see Figure 4, with

$$\delta {\mathit{a}}_{\alpha}=\delta {\mathit{\phi}}_{,\xi \alpha}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{sn}}}}\delta {\mathit{\phi}}_{i}{\mathrm{N\u0302}}_{,\xi \alpha}^{i}\text{and}{\mathit{a}}_{\alpha}={\mathit{\phi}}_{,\xi \alpha}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{sn}}}}{\mathit{\phi}}_{i}{\mathrm{N\u0302}}_{,\xi \alpha}^{i}\text{and}{\mathit{A}}_{\alpha}={\mathit{X}}_{,\xi \alpha}={\displaystyle {\sum}_{i=1}^{{\mathrm{n}}_{\text{sn}}}}{\mathit{X}}_{i}{\mathrm{N\u0302}}_{,\xi \alpha}^{i}.$$

(57)

Then, we determine the contravariant material base vectors *A*^{α} using the contravariant and covariant material metric coefficients *A*^{αβ} such that *A*_{αβ} for α,β = 1, 2, such that

$${\mathit{A}}^{\alpha}={A}^{\alpha \beta}{\mathit{A}}_{\beta}\text{with}{A}^{\alpha \beta}={\mathit{A}}^{\alpha}\xb7{\mathit{A}}^{\beta}={[{A}_{\alpha \beta}]}^{-1}\text{and}{A}_{\alpha \beta}={\mathit{A}}_{\alpha}\xb7{\mathit{A}}_{\beta}.$$

(58)

The gradients of the shape functions ∇̂*N ^{i}* follow in analogy to the volume, however, here, they are explicit functions of the contravariant base vectors

$$\nabla \u0302{N}^{i}=\frac{\mathrm{d}{\mathrm{N\u0302}}^{i}(\mathit{\xi})}{\mathrm{d}\mathit{X}}=\frac{\mathrm{d}{\mathrm{N\u0302}}^{i}(\mathit{\xi})}{\mathrm{d}\mathit{\xi}}\xb7\frac{\mathrm{d}\mathit{\xi}}{\mathrm{d}\mathit{X}}={\mathbf{\u0134}}^{-\mathrm{t}}\xb7\frac{\mathrm{d}{\mathrm{N\u0302}}^{i}(\mathit{\xi})}{\mathrm{d}\mathit{\xi}}\text{with}\mathbf{\u0134}=\frac{\mathrm{d}\mathit{X}}{\mathrm{d}\mathit{\xi}}=[{\mathit{A}}_{1},{\mathit{A}}_{2}]\text{and}{\mathbf{\u0134}}^{-1}=\frac{\mathrm{d}\mathit{\xi}}{\mathrm{d}\mathit{X}}=\left[\begin{array}{c}\hfill {\mathit{A}}^{1}\hfill \\ \hfill {\mathit{A}}^{2}\hfill \end{array}\right].$$

(59)

To systematically compare the features of volume and surface growth, we simulate constrained growth of the inner lining of a cylindrical hyperelastic tube. The tube is 5 cm long, and has an initial inner and outer radius of 1 cm and 2 cm, see Figure 5. Throughout the simulation, we apply homogeneous Dirichlet boundary conditions to the entire outer surface and to the axial deformation at the longitudinal boundaries. To initiate an inhomogeneous growth pattern, we perturb the inner surface with slight sinus-type perturbations of an amplitude of ± 0.005 cm.

Finite element discretization of cylindrical tube discretized with ten elements in the radial, 125 elements in the circumferential, and one element in the longitudinal direction, resulting in a total of 1,250 volume elements and 8,250 degrees of freedom. **...**

To illustrate the features of volume growth, we allow a thin inner layer of a cylindrical tube to grow, while the thick outer layer remains entirely hyperelastic. We characterize the hyperelastic response of both layers through an isotropic Helmholtz volume energy ψ_{0} of Neo-Hookean type,

$${\psi}_{0}=\frac{1}{2}\mu \phantom{\rule{thinmathspace}{0ex}}[{\mathit{F}}^{\mathrm{e}}:{\mathit{F}}^{\mathrm{e}}-3-2\text{ln}{J}^{\mathrm{e}}]+\frac{1}{2}\lambda {\text{ln}}^{2}{J}^{\mathrm{e}},$$

where λ and μ are the standard volume Lamé constants. To evaluate the discrete residual (50) and its consistent linearization (51), we calculate the volume Piola stress ** P** using the general definition (32),

$$\mathit{P}=\frac{\partial {\psi}_{0}}{\partial \mathit{F}}=\frac{1}{\vartheta}{\mathit{P}}^{\mathrm{e}}\text{with}{\mathit{P}}^{\mathrm{e}}=\frac{\partial {\psi}_{0}}{\partial {\mathit{F}}^{\mathrm{e}}}=\mu \phantom{\rule{thinmathspace}{0ex}}{\mathit{F}}^{\mathrm{e}}+[\lambda \text{ln}{J}^{\mathrm{e}}-\mu ]{[{\mathit{f}}^{\mathrm{e}}]}^{\mathrm{t}},$$

and the volume tangent using the general definition (33),

$$\mathbb{A}=\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{F}}=\frac{1}{{\vartheta}^{2}}{\mathbb{A}}^{\mathrm{e}}\text{with}{\mathbb{A}}^{\mathrm{e}}=\frac{\mathrm{d}{\mathit{P}}^{\mathrm{e}}}{\mathrm{d}{\mathit{F}}^{\mathrm{e}}}=\mu \phantom{\rule{thinmathspace}{0ex}}\mathit{I}\otimes \u0305\mathit{I}+[\mu -\lambda \text{ln}{J}^{\mathrm{e}}]{[{\mathit{f}}^{\mathrm{e}}]}^{\mathrm{t}}\otimes \u0332{\mathit{f}}^{\mathrm{e}}+\lambda \phantom{\rule{thinmathspace}{0ex}}{[{\mathit{f}}^{\mathrm{e}}]}^{\mathrm{t}}\otimes {[{\mathit{f}}^{\mathrm{e}}]}^{\mathrm{t}},$$

along with the assumption of isotropic volume growth, *F*^{g} = ** I**, according to equation (34). Typically, the inner layer is considered significantly thinner and stiffer than the outer layer [60]. For the thin growing inner layer, we choose a thickness of 0.1 cm, Lamé constants of λ = 577 N/mm

Table 1 demonstrates the algorithmic performance of the volume growth model. The five individual columns correspond to the snapshots displayed in Figure 6. The algorithm converges quadratically in five to six iteration steps for a convergence criterion of 10E−10. To accurately capture the growth process, we apply an ad hoc time adaptive scheme, which automatically decreases the time step size by 50% once the number of Newton iterations exceeds seven, and increases the time step size by 10% otherwise. As the growth process progresses, the time step size decreases initially and then increases. The minimum time step size of Δ*t*^{min} = 0.003 at a volume growth multiplier of = 1.38 corresponds exactly to the onset of buckling, when the five perturbations transition from an initially convex to an inward folded concave state.

Volume growth of inner layer of a cylindrical tube with pentagonally perturbed inner surface. Spatio-temporal evolution of radial stresses, top, and deviatoric stresses, bottom. The growing inner layer successively increases in volume, pulling the hyperelastic **...**

Figure 6 illustrates volume growth of the inner tubular layer at the discrete time points summarized in Table 1. The upper and lower rows display the spatio-temporal evolution of the radial and deviatoric stresses. As time progresses, the growing inner layer successively increases in volume, pulling the hyperelastic outer layer inward. During this growth process, the pentagonal perturbation grows into five distinct inner folds. Both radial and deviatoric stresses are concentrated in the growing inner layer, while the hyperelastic outer layer remains virtually stress free.

Figure 7 provides a fully three-dimensional view of the volume growth process. To explore the impact of the degree of perturbation, we systematically vary the number of perturbations applied to the inner surface. Square, pentagonal, hexagonal, heptagonal, and octagonal perturbations, displayed from left to right, gradually transition into four, five, six, seven, and eight distinct inner folds. The upper and lower rows display the spatial distribution of the radial and deviatoric stresses. Both stress patterns are conceptually similar for different degrees of perturbation, displaying maximum stresses in the growing inner, while the hyperelastic outer layer remains virtually stress free. The illustrated snapshots correspond to the time point of first self contact, which occurs at = 1.87 for square, at = 2.00 for pentagonal, at = 2.19 for hexagonal, at = 2.28 for heptagonal, and at = 2.36 for octagonal perturbations.

To illustrate the features of surface growth, we allow the inner surface of a cylindrical tube to grow, while the entire bulk remains hyperelastic. We characterize the hyperelastic response of the inner surface through an isotropic Helmholtz surface energy _{0} of Neo Hookean type in analogy to its volumetric counterpart [22],

$${\mathrm{\psi \u0302}}_{0}({\mathit{F\u0302}}^{\mathrm{e}})=\frac{1}{2}\mathrm{\mu \u0302}[{\mathit{F\u0302}}^{\mathrm{e}}:{\mathit{F\u0302}}^{\mathrm{e}}-2-2\text{ln}\widehat{{J}^{\mathrm{e}}}]+\frac{1}{2}\mathrm{\lambda \u0302}{\text{ln}}^{2}\widehat{{J}^{\mathrm{e}}},$$

where and are the surface Lamé constants. To evaluate the discrete residual (50) and its consistent linearization (51), we calculate the surface Piola stress ** using the general definition (41),
**

$$\mathit{P\u0302}=\frac{\partial {\psi}_{0}}{\partial \mathit{F\u0302}}=\frac{1}{\mathrm{\vartheta \u0302}}{\mathit{P\u0302}}^{\mathrm{e}}\text{with}{\mathit{P\u0302}}^{\mathrm{e}}=\frac{\partial {\psi}_{0}}{\partial {\mathit{F\u0302}}^{\mathrm{e}}}=\mathrm{\mu \u0302}{\mathit{F\u0302}}^{\mathrm{e}}+[\mathrm{\lambda \u0302}\text{ln}\widehat{{J}^{\mathrm{e}}}-\mathrm{\mu \u0302}]{[{\mathit{f\u0302}}^{\mathrm{e}}]}^{\mathrm{t}},$$

and the surface tangent using the definition (42),

$$\mathrm{\mathbb{A}\u0302}=\frac{\mathrm{d}\mathit{P\u0302}}{\mathrm{d}\mathit{F\u0302}}=\frac{1}{{\vartheta}^{2}}{\mathbb{A}}^{\mathrm{e}}\text{with}{\mathrm{\mathbb{A}\u0302}}^{\mathrm{e}}=\frac{\mathrm{d}{\mathit{P\u0302}}^{\mathrm{e}}}{\mathrm{d}{\mathit{F\u0302}}^{\mathrm{e}}}=\mathrm{\mu \u0302}\mathit{I}\otimes \u0305\mathit{\xce}+[\mathrm{\mu \u0302}-\mathrm{\lambda \u0302}\text{ln}\widehat{{J}^{\mathrm{e}}}]\phantom{\rule{thinmathspace}{0ex}}{[\phantom{\rule{thinmathspace}{0ex}}[{\mathit{f\u0302}}^{\mathrm{e}}]}^{\mathrm{t}}\otimes \u0332{\mathit{f\u0302}}^{\mathrm{e}}-{\mathit{\xee}}^{\perp}\otimes \u0305[{\mathit{f\u0302}}^{\mathrm{e}}\xb7{[{\mathit{f\u0302}}^{\mathrm{e}}]}^{\mathrm{t}}]]+\mathrm{\lambda \u0302}\phantom{\rule{thinmathspace}{0ex}}{[{\mathit{f\u0302}}^{\mathrm{e}}]}^{\mathrm{t}}\otimes {[{\mathit{f\u0302}}^{\mathrm{e}}]}^{\mathrm{t}},$$

along with the assumption of isotropic surface growth, ^{g} =** Î**, according to equation (43). In the above equation, we have used the abbreviation

Table 2 documents the algorithmic performance of the surface growth model. The five individual columns correspond to the snapshots displayed in Figure 8. The algorithm converges quadratically in five iteration steps for a convergence criterion of 10E−10. To accurately capture the growth process, we apply an ad hoc time adaptive scheme, which automatically decreases the time step size by 50% once the number of Newton iterations exceeds seven, and increases the time step size by 10% otherwise. As the growth process progresses, the time step size decreases initially and then increases. The minimum time step size of Δ*t*^{min} = 0.020 at a surface growth multiplier of = 1.13 corresponds exactly to the onset of buckling, when the five perturbations transition from an initially convex to an inward folded concave state.

Surface growth of inner surface of a cylindrical tube with pentagonally perturbed inner surface. Spatio-temporal evolution of radial stresses, top, and deviatoric stresses, bottom. The growing inner surface successively increases in area, pulling the **...**

Figure 8 illustrates surface growth of inner tubular surface at the discrete time points summerized in Table 2. The upper and lower rows display the spatio-temporal evolution of the radial and deviatoric stresses. As time progresses, the growing inner surface successively increases in area, pulling the elastic bulk inward. During this growth process, the pentagonal perturbation grows into five distinct inner folds. Both radial and deviatoric stresses display significant regional variations. While radial stresses are primarily concentrated in convex regions where the volume is stretched most, deviatoric stresses are concentrated mainly in concave regions where the volume is compressed most.

Figure 9 provides a fully three-dimensional view of the surface growth process. To explore the impact of the degree of perturbation, we systematically vary the number of perturbations applied to the inner surface. Square, pentagonal, hexagonal, heptagonal, and octagonal perturbations, displayed from left to right, gradually transition into four, five, six, seven, and eight distinct inner folds. The upper and lower rows display the spatial distribution of the radial and deviatoric stresses. Both stress patterns are conceptually similar for different degrees of perturbation, displaying maximum radial stresses in convex regions where the volume is stretched most, and maximum deviatoric stresses in concave regions where the volume is compressed most. The illustrated snapshots correspond to the time point of first self contact, which occurs consistent at = 2.75 for all five perturbations.

Growth in constrained geometries is an important phenomenon of equal relevance to material sciences, biology, and medicine. Here, we have compared two conceptually different types of growth: classical volume growth, which has been explored intensely within the past decade; and surface growth, which has been realized using the concept of continua with boundary energies. While volume growth can reasonably well model biological systems of finite thickness, surface growth is appropriate to characterize phenomena at the zero thickness limit. Our new continuum framework for growing surfaces is therefore particularly appealing for growing thin films and growing biological linings of several cell layers.

In this manuscript, we have established the kinematic equations, the balance equations, and the constitutive equations for continua with boundary energies and growing surfaces. To motivate this novel concept, we have derived all sets of equations in complete analogy to continua with growing volumes. We have discretized the underlying weak forms and illustrated their consistent algorithmic linearization. To systematically compare growing volumes and growing surfaces, we have simulated constrained growth of finite-thickness and zero-thickness linings of a hyperelastic cylindrical tube.

Figures 6 through through99 have shown that both volume and surface growth are capable of predicting extreme growth of the inner cylindrical lining, which more then doubled its initial area throughout the growth process. Both algorithms are robust and stable to simulate growth during the onset of buckling and beyond. Simulating the initiation of buckling instabilities is numerically highly challenging. To address this issue, we have adopted an ad hoc time adaptive scheme, which automatically reduces the time step size during phases when the Newton Raphson iteration does not converge within a limited number of iteration steps. This time adaptation has proven crucial to accurately track the initial stage of the instability and efficiently simulate growth beyond this critical point. However, this ad hoc adaptation scheme critically relies on a correct algorithmic linearization of the underlying highly nonlinear system of equations, which we have demonstrated through the characteristic quadratic convergence of the algorithmic residual in Table 1 and Table 2.

To quantitatively compare volume and surface growth, Figure 10 demonstrates the algorithmic performance of both formulations in terms of the adaptive time step size Δ*t*. The blue, turquoise, green, orange, and red curves correspond to the square, pentagonal, hexagonal, heptagonal, and octagonal perturbations displayed in Figures 7 and and9.9. With our particular time adaptive scheme, the required time step sizes for volume growth, left, were 0.002, 0.003, 0.004, 0.006, and 0.009 for the five simulations. The required time step sizes for surface growth, right, were 0.022, 0.020, 0.030, 0.030, and 0.025, i.e., they were in average an order of magnitude larger. Figure 11 compares both formulations in terms of the growth multiplier . The algorithm for volume growth, left, required 187, 154, 131, 114, and 98 simulation steps towards the point of first self contact, shown in Figure 7, while the algorithm for surface growth, right, required only 47, 47, 40, 41, and 42 steps towards this point, shown in Figure 9. Although these numbers strongly depend on the particular choice and parameterization of the time adaptation scheme, these studies indicates that surface growth is computationally more robust and more efficient than volume growth.

Algorithmic performance. Evolution of adaptive time step size Δ*t* for volume growth, left, and surface growth, right. The blue, turquoise, green, orange, and red curves correspond to the square, pentagonal, hexagonal, heptagonal, and octagonal **...**

Algorithmic performance. Evolution of growth multiplier for volume growth, left, and surface growth, right. The blue, turquoise, green, orange, and red curves correspond to the square, pentagonal, hexagonal, heptagonal, and octagonal perturbations **...**

The present formulation for continua with boundary energies and growing surfaces is particularly suited to model phenomena, in which the thickness of the growing membrane is significantly smaller than the overall dimensions of the system [16]. The thickness of the mucous membrane is typically less than a millimeter [44], the thickness of the cerebral cortex of the brain is of the order of millimeters [61], the thickness of the epidermis of skin is of the order of several cell layers [6]. During growth, neither the cortical [62] nor the epidermal [65] thicknesses undergo relevant changes, while their surfaces may grow significantly. While the proposed theory is well suited for thin biological membranes close to the zero thickness limit, which do not grow considerably in the thickness direction, volume growth with different growth multipliers in plane and out of plane [15, 43] might be more appropriate for thickening membranes. The role of thickness changes has been studied intensely in the context of mucosal folding [60]. Parametric studies focus on identifying critical thickness ratios [24], stiffness ratios [37], and growth multipliers [36] at the onset of wrinkling. Another group of studies attributes the onset of buckling to an increased muscular pressure at the outer surface and focuses on identifying critical pressure-to-luminal-areal [60] or critical pressure-to-growth [43] curves. While we have limited the current examples to growth-induced instabilities, we plan to investigate the impact of pressure-induced instabilities in the near future. Another limitation of the proposed approach is that the evolution of growth is prescribed constitutively through equations (35) and (44) as illustrated in Figure 3. Replacing the morphogenetically driven growth law by a mechanically driven growth law is conceptually possible. However, the corresponding linearization will become slightly more cumbersome, as we have recently shown in the context of volume growth [13, 64].

Many biological systems are coated by thin films, which display a fundamentally different behavior than the bulk. Here we have explored the physics and the numerics for growing thin biological surfaces. To model surface growth and the associated kinematic instabilities, we have compared two alternative concepts, classical volume growth for systems with a finite thickness and surface growth for systems with a zero thickness. Although the growing surfaces of both formations may look kinematically similar, aside from the underlying physics, there are three essential differences between the two schemes: First, and most obvious, surface growth is kinematically more localized than volume growth. It affects only the area of the inner lining, whereas volume growth also affects the thickness of the growing layer. Second, surface growth induces smoother stress profiles than volume growth. Its stresses are not just concentrated in the growing layer, but are distributed smoothly throughout the entire hyperelastic bulk. This implies that the peak stresses during surface growth are significantly smaller. This difference might become relevant in surface debonding, since elevated peak stresses are a hallmark for interface failure and an important indicator for layer separation. Third, surface growth is algorithmically more robust than volume growth. Surface growth accurately captures the onset of instabilities and traces the equilibrium path efficiently. Under the same conditions, the simulation of growing surfaces is about an order of magnitude faster than the simulation of growing volumes.

We believe that the concept of continua with boundary energies and growing surfaces will be widely applicable to predict growth-induced morphological instabilities in thin biological membranes, and, ultimately, support the design of new diagnostic tools and therapies for various types of respiratory, digestive, and urogenital disorders such as asthma, bronchitis, gastritis, obstructive sleep apnoea, and tumor invasion. In addition, the fundamental scientific understanding of growth-induced morphological instabilities has important implications in material sciences, manufacturing, and microfabrication, and may support the rational design of smart materials with tunable patterns and functional surfaces in soft lithography, metrology, and flexible electronics.

This study was supported by the National Science Foundation CAREER award CMMI 0952021 and INSPIRE grant 1233054 and by the National Institutes of Health Grant U54 GM072970.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1. Ambrosi D, Mollica F. On the mechanics of a growing tumor. Int J Eng Sci. 2002;40:1297–1316.

2. Ambrosi D, Ateshian GA, Arruda EM, Cowin SC, Dumais J, Goriely A, Holzapfel GA, Humphrey JD, Kemkemer R, Kuhl E, Olberding JE, Taber LA, Garikipati K. Perspectives on biological growth and remodeling. J Mech Phys Solids. 2011;59:863–883. [PMC free article] [PubMed]

3. Ben Amar M, Goriely A. Growth and instability in elastic tissues. J Mech Phys Solids. 2005;53:2284–2319.

4. Cao YP, Li B, Feng XQ. Surface wrinkling and folding of core-shell soft cylinders. Soft Matter. 2012;8:556–562.

5. McBride AT, Javili A, Steinmann P, Bargmann S. Geometrically nonlinear continuum theories with surface energies coupled to diffusion. J Mech Phys Solids. 2011;59:2116–2133.

6. Buganza Tepole A, Ploch CJ, Wong J, Gosain AK, Kuhl E. Growing skin - A computational model for skin expansion in reconstructive surgery. J Mech Phys Solids. 2011;59:2177–2190. [PMC free article] [PubMed]

7. Buganza Tepole A, Gosain AK, Kuhl E. Stretching skin: The physiological limit and beyond. Int J Nonlin Mech. 2012;47:938–949. [PMC free article] [PubMed]

8. Davies DE, Wicks J, Powell RM, Puddicombe SM, Holgate ST. Airway wall remodeling in asthma: New insights. J Allergy Clin Immunol. 2003;111:215–225. [PubMed]

9. Dervaux J, Ciarletta P, Ben Amar M. Morphogenesis of thin hyperelastic plates: A constitutive theory of biological growth in the Föppl-von Karman limit. J Mech Phys Solids. 2009;57:458–471.

10. Dervaux J, Ben Amar M. Buckling considerations in constrained growth. J Mech Phys Solids. 2011;59:538–560.

11. Epstein M, Maugin GA. Thermomechanics of volumetric growth in uniform bodies. Int J Plast. 2000;16:951–978.

12. Garikipati K. The kinematics of biological growth. Appl Mech Rev. 2009;62:030801.1–030801.7.

13. Göktepe S, Abilez OJ, Parker KK, Kuhl E. A multiscale model for eccentric and concentric cardiac growth through sarcomerogenesis. J Theor Bio. 2010;265:433–442. [PubMed]

14. Göktepe S, Abilez OJ, Kuhl E. A generic approach towards finite growth with examples of athlete’s heart, cardiac dilation, and cardiac wall thickening. J Mech Phys Solids. 2010;58:1661–1680.

15. Goriely A, BenAmar M. Differential growth and instability in elastic shells. Phys Rev Letters. 2005;94:198103. [PubMed]

16. Gurtin ME, Murdoch A. A continuum theory of elastic material surfaces. Arch Rational Mech Anal. 1975;57:291323.

17. Gurtin ME, Struthers A. Multiphase thermomechanics with interfacial structure. Arch Rational Mech Anal. 1990;112:97160.

18. Himpel G, Kuhl E, Menzel A, Steinmann P. Computational modeling of isotropic multiplicative growth. Comp Mod Eng Sci. 2005;8:119–134.

19. Hogg J. Pathophysiology of airflow limitation in chronic obstructive pulmonary disease. Lancet. 2004;364:709–721. [PubMed]

20. Javili A, Steinmann P. A finite element framework for continua with boundary energies. Part I: The two-dimensional case. Comp Meth Appl Mech Eng. 2009;198:2198–2208.

21. Javili A, Steinmann P. A finite element framework for continua with boundary energies. Part II: The three-dimensional case. Comp Meth Appl Mech Eng. 2010;199:755–765.

22. Javili A, Steinmann P. On thermomechanical solids with boundary structures. Int J Solids Struct. 2010;47:3245–3253.

23. Javili A, Steinmann P. A finite element framework for continua with boundary energies. Part III: The thermomechanical case. Comp Meth Appl Mech Eng. 2011;200:1963–1977.

24. Jin L, Cai S, Suo Z. Creases in soft tissues generated by growth. EPL Frontiers of Physics. 2011;95:64002.1–64002.6.

25. Kairaitis K. Pharyngeal wall fold influences on the collapsibility of the pharynx. Med Hypotheses. 2012;79:327–376. [PubMed]

26. Kroon W, Delhaas T, Arts T, Bovendeerd P. Computational modeling of volumetric soft tissue growth: Application to the cardiac left ventricle. Biomech Model Mechanobio. 2009;8:310–309. [PubMed]

27. Kuhl E, Steinmann P. Mass- and volume specific views on thermodynamics for open systems. Proc Royal Soc. 2003;459:2547–2568.

28. Kuhl E, Steinmann P. On spatial and material settings of thermohyperelstodynamics for open systems. Acta Mech. 2003;160:179–217.

29. Kuhl E, Maas R, Himpel G, Menzel A. Computational modeling of arterial wall growth: Attempts towards patient-specific simulations based on computer tomography. Biomech Mod Mechanobio. 2007;6:321–331. [PubMed]

30. Kuhl E, Schmid DW. Computational modeling of mineral unmixing and growth: An application of the Cahn-Hilliard equation. Comp Mech. 2007;39:439–451.

31. Laplace PS. Oeuvres complète. Suppl. 1. Vol. 4. Paris: Gauthier-Villars; Traité de Mécanique Céleste; p. 1805. Livre X:771–777.

32. Lee EH. Elastic-plastic deformation at finite strains. J Appl Mech. 1969;36:1–6.

33. Lee MML, Chien S. Morphologic effects of pressure changes on canine carotid artery endothelium as observed by scanning electron microscopy. Anat Rec. 1978;194:1–14. [PubMed]

34. Leo PH, Sekerka RF. The effect of surface stress on crystal-melt and crystal-crystal equilibrium. Acta Metall. 1989;37:3119–3138.

35. Leo PH, Sekerka RF. The effect of elastic fields on the morphological stability of a precipitate grown from solid solution. Acta Metall. 1989;37:3139–3149.

36. Li B, Cao YP, Feng XQ, Gao H. Surface wrinkling of mucosa induced by volumetric growth: Theory, simulation and experiment. J Mech Phys Solids. 2011;59:758–774.

37. Li B, Cao YP, Feng XQ, SW Yu. Mucosal wrinkling in animal antra induced by volumetric growth. Appl Phys Letters. 2011;98:153701.

38. Li B, Cao YP, Feng XQ, Gao H. Mechanics of morphological instabilities and surface wrinkling in soft materials: A review. Soft Matter. 2012;8:5728–5745.

39. Liao D, Zhao J, Gregersen H. The oesophageal zero-stress state and mucosal folding from a giome perspective. World J Gastroenterol. 2007;13:1347–1351. [PubMed]

40. McMahon J, Goriely A. Spontaneous cavitation in growing elastic membranes. Math Mech Solids. 2010;15:57–77.

41. Menzel A. Modelling of anisotropic growth in biological tissues - A new approach and computational aspects. Biomech Model Mechanobiol. 2005;3:147–171. [PubMed]

42. Menzel A, Kuhl E. Frontiers in growth and remodeling. Mech Res Comm. 2012;42:1–14. [PMC free article] [PubMed]

43. Moulton DE, Goriely A. Circumferential buckling instability of a growing cylindrical tube. J Mech Phys Solids. 2011;59:525–537.

44. Moulton DE, Goriely A. Possible role of differential growth in airway wall remodeling in asthma. J Appl Physiol. 2011;110:1003–1012. [PubMed]

45. Nelson MR, Howard D, Jensen OE, King JR, Rose FRAJ, Waters SL. Growth-induced buckling of an epithelial layer. Biomech Model Mechanobio. 2011;10:883–900. [PubMed]

46. Pare PD, McParland BE, Seow CY. Structural basis for exaggerated airway narrowing. Can J Physiol Pharmacol. 2007;85:653–658. [PubMed]

47. Rausch MK, Dam A, Göktepe S, Abilez OJ, Kuhl E. Computational modeling of growth: Systemic and pulmonary hypertension in the heart. Biomech Model Mechanobiol. 2011;10:799–811. [PMC free article] [PubMed]

48. Rausch MK, Tibayan FA, Miller DC, Kuhl E. Evidence of adaptive mitral leaflet growth. J Mech Behavior Biomed Mat. 2012;15:208–217. [PMC free article] [PubMed]

49. Rodriguez EK, Hoger A, McCulloch AD. Stress-dependent finite growth in soft elastic tissues. J Biomech. 1994;27:455–467. [PubMed]

50. Saksono PH, Peric D. On finite element modelling of surface tension. Variational formulation and applications. Part I: Quasistatic problems. Comp Mech. 2006;38:251–263.

51. Simha NK, Bhattacharya K. Kinetics of phase boundaries with edges and junctions. J Mech Phys Solids. 1998;46:2323–2359.

52. Simha NK, Bhattacharya K. Kinetics of phase boundaries with edges and junctions in a three-dimensional multi-phase body. J Mech Phys Solids. 2000;48:2619–2641.

53. Socci L, Rennati G, Gervaso F, Vena P. An axisymmetric computational model of skin expansion and growth. Biomech Model Mechanobiol. 2007;6:177–188. [PubMed]

54. Steinmann P. On boundary potential energies in deformational and configurational mechanics. J Mech Phys Solids. 2008;56:772–800.

55. Suo Z. Motions of microscopic surfaces in materials. Advances Appl Mech. 33:193–294. 19978.

56. Taber LA. Biomechanics of growth, remodeling and morphogenesis. Appl Mech Rev. 1995;48:487–545.

57. Taber LA, Humphrey JD. Stress-modulated growth, residual stress, and vascular heterogeneity. J Biomech Eng. 2001;123:528–535. [PubMed]

58. Tsamis A, Cheng A, Nguyen TC, Langer F, Miller DC, Kuhl E. Kinematics of cardiac growth - In vivo characterization of growth tensors and strains. J Mech Behavior Biomed Mat. 2012;8:165–177. [PMC free article] [PubMed]

59. Vandiver R, Goriely A. Differential growth and residual stress in cylindrical elastic structures. Phil Trans R Soc A. 2009;367:3607–3630. [PubMed]

60. Wiggs BR, Hrousis CA, Drazen JM, Kamm RD. On the mechanism of mucosal folding in normal and asthmatic airways. J Appl Physiol. 1997;83:1814–1821. [PubMed]

61. Wyczalkowski MA, Chen Z, Filas BA, Varner VD, Taber LA. Computational models for mechanics of morphogenesis. Birth Defects Res C. 2012;96:132–152. [PMC free article] [PubMed]

62. Xu G, Knutsen AK, Dikranian K, Kroenke CD, Bayly PV, Taber LA. Axons pull on the brain, but tension does not drive cortical folding. J Biomech Eng. 2010;132(7) 071013. [PMC free article] [PubMed]

63. Young T. An essay on the cohesion of fluids. Phil Trans R Soc London. 1805;95:65–87.

64. Zöllner AM, Buganza Tepole A, Kuhl E. On the biomechanics and mechanobiology of growing skin. J Theor Bio. 2012;297:166–175. [PMC free article] [PubMed]

65. Zöllner AM, Buganza Tepole A, Gosain AK, Kuhl E. Growing skin - Tissue expansion in pediatric forehead reconstruction. Biomech Mod Mechanobio. 2012;11:855–867. [PMC free article] [PubMed]

66. Zöllner AM, Abilez OJ, Böl M, Kuhl E. Stretching skeletal muscle: Chronic muscle lengthening through sarcomerogenesis. PLoS ONE. 2012;7(10):e45661. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |