PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
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.007
PMCID: PMC3627422
NIHMSID: NIHMS443042

On the mechanics of continua with boundary energies and growing surfaces

Abstract

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.

Keywords: boundary energy, surface growth, volume growth, instability, finite elements, airway wall remodeling, thin films

1. Introduction

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].

Figure 1
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 ...

1.1. The mechanical origin of morphological instabilities

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].

1.2. Continuum modeling of mechanically-induced biological growth

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].

1.3. Continuum modeling of growth-induced mechanical instabilities

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.

1.4. Continuum modeling of material surfaces

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.

1.5. Organization of this manuscript

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.

Kinematics of growth

2.1. Kinematics of volume growth

Let B0 [set membership] R3 denote the material placement of a continuum body with material particles X [subset or is implied by] B0. We denote the time domain by T = [0, T ] [subset or is implied by] R+ 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 [var phi] : B0 × TR3 as the smooth motion of the material placement X onto its spatial placement x at any time t [set membership] T,

equation M1
(1)

and label the current placement of the body at time t as Bt = [var phi](B0), 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, t) as

equation M2
(2)

where I is the material identity tensor, and {o}|t, and {o}|X indicate fixed temporal and spatial positions. Using the above definitions, we introduce the deformation gradient F : TB0TBt, which maps material line elements dX [set membership] TB0 onto spatial line elements dx [set membership] TBt as dx = F · dX. To account for volumetric growth, we multiplicatively decompose the deformation gradient into an elastic part Fe and a growth part Fg [49],

equation M3
(3)

Figure 2
Kinematics of material and spatial configurations B0 and Bt of material body with material and spatial surfaces S0 and St and multiplicative decomposition of the volume and surface deformation gradients F = Fe · ...

In contrast to the deformation gradient F, however, its individual contributions Fe and Fg can, in general, be incompatible [32]. For notational convenience, we also introduce the inverses of these second order tensors as f = [F]−1, fe = [F]e − 1, and fg = [F]g − 1.

The Jacobian J of the deformation gradient F relates material volume elements dV0 [set membership] B0 and spatial volume elements dVt [set membership] Bt as dVt = JdV0. Similar to the deformation gradient, we can decompose it multiplicatively into an elastic part Je and a growth part Jg,

equation M4
(4)

We can explicitly express the elastic Jacobian Je = det (Fe) and the growth Jacobian Jg = det (Fg), where we introduce the following explicit representation in terms of the I = 1,2,3 covariant base vectors GI in TB0,

equation M5
(5)

Last, we introduce the right Cauchy Green tensor C in the material configuration and its elastic counterpart Ce in the intermediate configuration as characteristic deformation measures,

equation M6
(6)

which are related through the identity Ce = [fg]t · C · fg. With the above considerations, the pull back of the spatial velocity gradient l to the intermediate configuration

equation M7
(7)

obeys an additive decomposition into the elastic velocity gradient Le = fe · DtFe and the growth velocity gradient Lg = DtFg · fg.

2.2. Kinematics of surface growth

Let S0 = [partial differential]B0 denote the surface of B0, assumed smooth, such that B0 = B0 [union or logical sum] S0 is the closure of B0 and S0. We label points on the surface of the reference configuration S0 as X = X|s0, and the unit outward normal to S0 as N, see Figure 2. We can then introduce the surface deformation map [var phiv with circumflex] : S0 × TR3 as the smooth motion of points on the surface of the material configuration X onto points on the surface of the spatial configuration x̂ = x|st at any time t [set membership] T [54],

equation M8
(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} (X, t) defined on the surface,

equation M9
(9)

where Î = IN [multiply sign in circle] N is the second order surface unit tensor, which serves as surface projection operator. With the above definitions, we introduce the surface deformation gradient F : S0St, which maps tangential line elements to the material surface dX [set membership] S0 onto tangential line elements to the spatial surface dx̂ [set membership] St. Similar to the volume deformation gradient F, the surface deformation gradient F obeys a multiplicative decomposition into an elastic part Fe and a growth part Fg,

equation M10
(10)

In contrast to the surface deformation gradient F, however, neither Fe nor Fg are necessarily restricted to result from the projection of the corresponding volume quantities Fe and Fg. Although the surface deformation gradient and its elastic and growth part are non-invertible, they possess inverses in the generalized sense, f · F = Î, fe · Fe = Î, and fg · Fg = Î in terms of the surface unit tensor Î. In analogy to the volume Jacobian J, we can introduce the Jacobian Ĵ of the surface deformation gradient F, which relates material surface area elements dA0 [set membership] S0 and spatial surface area elements dAt [set membership] St as dAt = ĴdA0. Similarly, we decompose the surface Jacobian multiplicatively into an elastic part equation M11 and a growth part equation M12,

equation M13
(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 equation M14 and the growth surface Jacobian as equation M15 in terms of the covariant base vectors Aα in TS0 [17],

equation M16
(12)

Finally, we introduce the right Cauchy Green surface tensor Ĉ in the material configuration and its elastic counterpart Ĉe in the intermediate configuration as characteristic surface deformation measures,

equation M17
(13)

which are related through the identity, Ĉe = fgt · Ĉ · fg. In analogy to the volume considerations, we introduce the pull back of the spatial surface velocity gradient l to the intermediate configuration,

equation M18
(14)

which obeys an additive decomposition into the elastic surface velocity gradient Le = fe · DtFe and the growth surface velocity gradient Lg = DtFg · fg.

Remark 1. The material and spatial second order surface unit tensors Î = IN [multiply sign in circle] N and î = In [multiply sign in circle] n serve as projection tensors to map the volume deformation gradient and its inverse onto their surface counterparts,

equation M19
(15)

equation M20

where N and n = N · f/|N · f|are the unit outward normals to the material and spatial surfaces. Since the surface unit tensors Î and î are rank deficient, the surface deformation gradient is non-invertible. However, it possesses a generalized inverse according to the following singular value decomposition,

equation M21
(16)

equation M22

where the diagonal entries of Σ and σ correspond to the singular values of F and Ft, the columns of U and V are the left- and right-singular vectors associated with these singular values, and Σ+ and σ+ are the pseudoinverses of Σ and σ, which are formed by replacing every non-zero diagonal entry by its reciprocal value.

3. Balance equations of growth

3.1. Balance equations of volume growth

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

equation M23
(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,

equation M24
(18)

In the transformation above, we have utilized the time derivative of Jg as DtJg = [partial differential]FgJg : DtFg = Jg[fg]t : DtFg = Jg tr Lg. The combination of equations (17.1) and (18.2) renders the following explicit representation of the mass source R0

equation M25
(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 b,

equation M26
(20)

The Piola stresses P and volume forces b have the dimensions force per unit area and force per unit volume. We obtain the weak form of the balance of linear momentum (20) through the standard multiplication with the test function δ[var phi], integration over the volume B0, and integration by parts,

equation M27
(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 σ,

equation M28
(22)

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

equation M29
(23)

3.2. Balance equations of surface growth

The balance of mass for growing surfaces balances the rate of change of the material area density [rho with circumflex]0 with the area mass source R0 and the area mass flux R, which we again assume to vanish in the sequel, R = 0, such that

equation M30
(24)

We can then relate the surface densities in the material, spatial, and intermediate configurations [rho with circumflex]0, [rho with circumflex]t, and [rho with circumflex]g with the help of the surface Jacobians (11)

equation M31
(25)

where equation M32. The mass source then follows directly from the combination of equations (24.1) and (25.2) as

equation M33
(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 P and the surface tractions bP · N,

equation M34
(27)

Here, the surface stresses P and the surface tractions bP · N have the dimensions force per unit length and force per unit area. The surface tractions consist of two contributions, the prescribed surface tractions b and the surface tractions imposed by the underlying volume through the projected volume Piola stress P · N [54]. To transform the balance of linear surface momentum into its weak form, we multiply it with the test function δ[var phi], integrate it over the surface S0, and apply the integration by parts,

equation M35
(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 C0. Again, we can introduce the surface Piola Kirchhoff stress Ŝ and the surface Cauchy stress [sigma with hat] using the surface Piola transforms,

equation M36
(29)

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

equation M37
(30)

4. Constitutive equations

4.1. Constitutive equations of volume growth

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 Fe, or, alternatively, in terms of the volume deformation gradient F and the volume growth tensor Fg,

equation M38
(31)

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

equation M39
(32)

where Pe denotes the classic elastic volume Piola stress. The total derivative of the Piola stress P with respect to the deformation gradient F introduces the fourth order tensor of tangent volume moduli A,

equation M40
(33)

where Ae denotes the classic elastic volume tangent moduli. To characterize the growth process, we need to specify a functional form for the volume growth tensor Fg and its evolution in time. For the simplest case of isotropic volume growth, the volume growth tensor Fg takes the following functional form, parameterized in terms of a single scalar-valued growth multiplier [theta],

equation M41
(34)

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

equation M42
(35)

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

Figure 3
Temporal evolution of growth multiplier [theta] for varying growth velocities τ at fixed [theta]max, left, and varying maximum growth [theta]max at fixed τ, right.

Remark 2 (Isotropic volume growth). The special case of isotropic volume growth, Fg = [theta]I, has three immediate consequences. First, the volume growth tensor Fg can be inverted explicitly, and the volume deformation gradient F, the volume Piola stress P, and the volume tangent moduli A become equal to their growth-weighted elastic counterparts,

equation M43
(36)

Second, the growth tensor Fg always remains isotropic, and so does the volume growth velocity gradient Lg,

equation M44
(37)

Third, volume changes upon growth Jg are simply equivalent to the third power of the growth multiplier [theta],

equation M45
(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], Dtρg = 0. This implies that the volume density source R0 required to maintain a constant volume density in the intermediate configuration takes the following form,

equation M46
(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.

4.2. Constitutive equations of surface growth

To characterize surface growth, we equip the surface with its own Helmholtz free surface energy [psi], which we parameterize in terms of the elastic contribution to the surface deformation gradient Fe, or, alternatively, in terms of the surface deformation gradient F and the surface growth tensor Fg,

equation M47
(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 P as thermodynamically conjugate to the surface deformation gradient F,

equation M48
(41)

where Pe denotes the elastic surface Piola stress. Additionally, we introduce the fourth order tensor of tangent surface moduli A as the total derivative of the Piola stress P with respect to the deformation gradient F,

equation M49
(42)

such that Ae denotes the elastic surface tangent moduli. Again, we consider the simplest possible case of growth, isotropic surface growth, for which the surface growth tensor Fg can be parameterized in terms of a single scalar-valued growth multiplier [theta],

equation M50
(43)

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

equation M51
(44)

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

Remark 4 (Isotropic surface growth). The special case of isotropic surface growth, Fg = [theta]Î, has three immediate consequences. First, we can explicitly invert the surface growth tensor Fg 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 A through to their growth-weighted elastic counterparts,

equation M52
(45)

Second, the surface growth tensor Fg and the surface growth velocity gradient Lg always remain isotropic,

equation M53
(46)

Third, surface area changes upon growth equation M54 are simply equivalent to the growth multiplier [theta] squared,

equation M55
(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, Dt[rho with circumflex]g = 0, we can express the surface density source R0 in the following form,

equation M56
(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.

5. Discretization

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 [theta]g and [theta]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),

equation M57
(49)

To discretize the weak form (49) in space, we partition the domain of interest equation M58 into nbe finite volume elements equation M59, and the surface of interest equation M60 into nse finite surface elements equation M61. 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

equation M62
(50)

Figure 4
Discretization of material body with finite volume elements equation M79 and equation M80 and of its surface with finite surface elements equation M81 and equation M82. Surface elements share their nodes with the corresponding volume elements, however, they are equipped with their own energies [psi], ...

Herein, the operator A symbolizes the assembly of all element residuals at the i = 1, ‥, nvn volume element nodes to the global residual at the global node points I = 1, ‥, ngn. To solve the above equation, we suggest an incremental iterative Newton algorithm based on the consistent linearization of the residual RI with respect to the nodal vector of unknowns [var phi]J. This linearization introduces the global stiffness matrix at all global nodes I, J = 1, ‥, nng,

equation M63
(51)

For each global Newton iteration step, we iteratively update the current deformation state equation M64 until we achieve algorithmic convergence. Upon convergence, we store the volume and surface growth multipliers [theta]g and [theta]g at the integration points of the corresponding volume and surface elements.

5.1. Discretization of volume terms

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

equation M65
(52)

where Ni and Nj are the element shape functions in the volume and i, j = 1, ‥, nvn are nodes of the volume element. The gradients of the test and trial functions δF = [nabla]δ[var phi] and F = [nabla][var phi] then follow naturally in terms of the gradients of the shape functions [nabla]Ni and [nabla]Nj,

equation M66
(53)

To determine the gradients of the shape functions [nabla]Ni, we use the isoparameteric concept for the approximation of the nodal coordinates X from (52.3) to calculate the Jacobi matrix J and its transposed inverse J−t,

equation M67
(54)

5.2. Discretization of surface terms

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

equation M68
(55)

where Ni and Nj are the element shape functions on the surface and i, j = 1, ‥, nsn are nodes of the surface element. The gradients of the test and trial functions δF = ∇̂δ[var phi] and F = ∇̂[var phi] on the surface follow in terms of the covariant spatial base vectors δaα and aα multiplied with the contravariant material base vectors Aα, for α = 1,2, such that

equation M69
(56)

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

equation M70
(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

equation M71
(58)

The gradients of the shape functions ∇̂Ni follow in analogy to the volume, however, here, they are explicit functions of the contravariant base vectors Aα,

equation M72
(59)

6. Results

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.

Figure 5
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. ...

6.1. Volume growth inside a cylindrical tube

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,

equation M73

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),

equation M74

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

equation M75

along with the assumption of isotropic volume growth, Fg = [theta]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/mm2 and μ = 385 N/mm2, a maximum growth of [theta]max = 4.0, and a time constant of τ = 10.0, see Figure 3. For the thick elastic outer layer, the thickness is 0.9 cm, and we choose Lamé constants of λ = 0.577 N/mm2 and μ = 0.385 N/mm2. We initiate an inhomogeneous growth pattern through a pentagonal perturbation of five sinus waves applied to the inner surface. We discretize the tube 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. Figure 5, left, shows the finite element discretization for volume growth, with the growing inner layer displayed in gray and the elastic outer layer displayed in white.

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 Δtmin = 0.003 at a volume growth multiplier of [theta] = 1.38 corresponds exactly to the onset of buckling, when the five perturbations transition from an initially convex to an inward folded concave state.

Figure 6
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 ...
Table 1
Algorithmic performance. Quadratic convergence of the volume growth model.

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 [theta] = 1.87 for square, at [theta] = 2.00 for pentagonal, at [theta] = 2.19 for hexagonal, at [theta] = 2.28 for heptagonal, and at [theta] = 2.36 for octagonal perturbations.

Figure 7
Volume growth of inner layer of a cylindrical tube with perturbed inner surface. Spatial distribution of radial stresses, top, and deviatoric stresses, bottom. Square, pentagonal, hexagonal, heptagonal, and octagonal perturbations transition into four, ...

6.2. Surface growth inside a cylindrical tube

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 [psi]0 of Neo Hookean type in analogy to its volumetric counterpart [22],

equation M76

where [lambda with circumflex] and [mu] are the surface Lamé constants. To evaluate the discrete residual (50) and its consistent linearization (51), we calculate the surface Piola stress P using the general definition (41),

equation M77

and the surface tangent A using the definition (42),

equation M78

along with the assumption of isotropic surface growth, Fg =[theta]Î, according to equation (43). In the above equation, we have used the abbreviation i[perpendicular] for the second order normal projection tensor i[perpendicular] = iî = n [multiply sign in circle] n. For the growing surface, we choose surface Lamé constants of [lambda with circumflex] = 0.577 N/mm2 and [mu] = 0.385 N/mm2, a maximum growth of [theta]max = 4.0, and a time constant of [tau] = 10.0, see Figure 3. For the elastic bulk, we choose volume Lamé constants of λ = 0.577 N/mm2 and μ = 0.385 N/mm2. We initiate an inhomogeneous growth pattern through a pentagonal perturbation of five sinus waves applied to the inner growing surface. We discretize the tube 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. In addition, we discretize the inner surface with 125 surface elements. Since these surface elements share their nodes with the corresponding volume elements, see Figure 4, the addition of surface elements does not change the overall number of degrees of freedom. Figure 5, right, shows the finite element discretization for surface growth, with the growing inner surface displayed in gray and the elastic bulk displayed in white.

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 Δtmin = 0.020 at a surface growth multiplier of [theta] = 1.13 corresponds exactly to the onset of buckling, when the five perturbations transition from an initially convex to an inward folded concave state.

Figure 8
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 ...
Table 2
Algorithmic performance. Quadratic convergence of the surface growth model.

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 [theta] = 2.75 for all five perturbations.

Figure 9
Surface growth of inner surface of a cylindrical tube with perturbed inner surface. Spatial distribution of radial stresses, top, and deviatoric stresses, bottom. Square, pentagonal, hexagonal, heptagonal, and octagonal perturbations transition into four, ...

7. Discussion

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 [theta]. 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.

Figure 10
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 ...
Figure 11
Algorithmic performance. Evolution of growth multiplier [theta] 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].

8. Conclusion

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.

Acknowledgements

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.

Footnotes

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.

References

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]