Nonlinear elasticity of incompatible surface growth

TL;DR

This paper extends the linear elasticity theory of incompatible surface growth to include nonlinear effects, revealing limitations of the linear model in describing confined growth and elastic instabilities in elastic solids.

Contribution

It introduces a nonlinear elasticity framework for incompatible surface growth, addressing the shortcomings of previous linearized models.

Findings

Linear theory cannot describe kinematically confined growth.

Nonlinear effects are crucial for understanding growth-induced instabilities.

The new model captures complex post-growth stress states.

Abstract

Surface growth is a crucial component of many natural and artificial processes from cell proliferation to additive manufacturing. In elastic systems surface growth is usually accompanied by the development of geometrical incompatibility leading to residual stresses and triggering various instabilities. In a recent paper (PRL, 119, 048001, 2017) we developed a linearized elasticity theory of incompatible surface growth which quantitatively linked deposition protocols with post-growth states of stress. Here we extend this analysis to account for both physical and geometrical nonlinearities of an elastic solid. The new development reveals the shortcomings of the linearized theory, in particular, its inability to describe kinematically confined surface growth and to account for growth-induced elastic instabilities.

Equations

$$\left\{\begin{array}[]{lll}\text{Div}\textbf{S}(\textbf{\emph{x}})+\textbf{\emph{f}}(\textbf{\emph{x}})=\textbf{0}&&\mathring{\mathcal{B}}\\ \textbf{S}(\textbf{\emph{x}})\mathring{\textbf{\emph{n}}}(\textbf{\emph{x}})=\textbf{\emph{s}}(\textbf{\emph{x}})&&\partial\mathring{\mathcal{B}}\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\text{Div}\textbf{S}(\textbf{\emph{x}},t)+\textbf{\emph{f}}(\textbf{\emph{x}},t)=\textbf{0}&&\mathring{\mathcal{B}}_{t}\\ \textbf{S}(\textbf{\emph{x}},t)\mathring{\textbf{\emph{n}}}(\textbf{\emph{x}},t)=\textbf{\emph{s}}(\textbf{\emph{x}},t)&&\partial\mathring{\mathcal{B}}_{t}\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\text{Div}\textbf{S}(\textbf{\emph{x}},t)=0&&\mathring{\mathcal{B}}_{t}\\ \textbf{S}(\textbf{\emph{x}},t)\mathring{\textbf{\emph{n}}}(\textbf{\emph{x}},t)=0&&\partial\mathring{\mathcal{B}}_{t}.\end{array}\right.$$

$$\overline{\text{S}}(\text{x})={\text{S}}_a(\text{x})+{\text{S}}_p(\text{x})\text{x}\in \partial {\mathring{\mathcal{B}}}_t$$

$$\nabla \bar{\text{y}}-\overline{\nabla \text{y}}=||\nabla \mathring{\vartheta }||\overline{\dot{\text{y}}}\otimes \mathring{\text{n}}$$

$$e(\text{G},\mathring{\text{G}})=\frac{E}{2(1+\nu )}\left(\text{tr}({\mathbf{ϵ}}_e^2)+\frac{\nu }{1-2\nu }(\text{tr}{\mathbf{ϵ}}_e{)}^2\right).$$

$$\begin{array}[]{ccc}\displaystyle\nabla\textbf{\emph{y}}(\textbf{\emph{x}},t)=\overline{\nabla\textbf{\emph{y}}}(\textbf{\emph{x}})+\int_{\mathring{\vartheta}(\textbf{\emph{x}})}^{t}\nabla\dot{\textbf{\emph{y}}}(\textbf{\emph{x}},s)\,ds\\ \displaystyle\textbf{S}(\textbf{\emph{x}},t)=\overline{\textbf{S}}(\textbf{\emph{x}})+\int_{\mathring{\vartheta}(\textbf{\emph{x}})}^{t}\dot{\textbf{S}}(\textbf{\emph{x}},s)\,ds.\end{array}$$

$$\text{f}(\text{x},t)=\bar{\text{f}}(\text{x})+{\int }_{\mathring{\vartheta }(\text{x})}^t\dot{\text{f}}(\text{x},s)ds$$

$$\left\{\begin{array}[]{lcc}\text{Div}\dot{\textbf{S}}+\dot{\textbf{\emph{f}}}=\textbf{\emph{0}}&&\mathring{\mathcal{B}}_{t}\\ \dot{\textbf{S}}\mathring{\textbf{\emph{n}}}=\mathring{D}\left(\text{Div}\overline{\textbf{S}}+\overline{\textbf{\emph{f}}}\right)&&\partial\mathring{\mathcal{B}}_{t}.\end{array}\right.$$

$$\mathring{\mathbf{\eta }}=\overline{\mathbf{\eta }}-\text{Curl}\text{Curl}(\text{sym}(\overline{\dot{\text{u}}}\otimes \nabla \mathring{\vartheta })),$$

$$\left\{\begin{array}[]{lll}\text{Div}\textbf{S}+\textbf{\emph{f}}=\textbf{\emph{0}}&&\mathring{\mathcal{B}}_{t}\\ \text{Curl}\text{Curl}(\mathbb{C}^{-1}\textbf{S})+\mathring{\boldsymbol{\eta}}=\textbf{\emph{0}}&&\mathring{\mathcal{B}}_{t}\\ \textbf{S}\textbf{\emph{n}}=\textbf{\emph{s}}&&\partial\mathring{\mathcal{B}}_{t}.\end{array}\right.$$

$$\nabla \text{y}(\text{x},t)={\lambda }_r(r,t)\text{e}(\text{x})\otimes \text{e}(\text{x})+{\lambda }_{\theta }(r,t)\mathbb{P}(\text{x})$$

$$\mathring{\text{G}}(\text{x})={\gamma }_r^2(r)\text{e}(\text{x})\otimes \text{e}(\text{x})+{\gamma }_{\theta }^2(r)\mathbb{P}(\text{x}).$$

$$S^d=2K=2\frac{{\gamma }_{\theta }{\gamma }_r^{\prime }+r{\gamma }_r^{\prime }{\gamma }_{\theta }^{\prime }-{\gamma }_r(2{\gamma }_{\theta }^{\prime }+r{\gamma }_{\theta }^{\prime \prime }))}{r{\gamma }_r^3{\gamma }_{\theta }}.$$

$${\gamma }_r=k(r{\gamma }_{\theta }{)}^{\prime }$$

$$\left\{\begin{array}[]{lll}R_{1}^{s}=2\frac{\gamma_{\theta}\gamma_{r}^{\prime}+r\gamma_{r}^{\prime}\gamma_{\theta}^{\prime}-\gamma_{r}(2\gamma_{\theta}^{\prime}+r\gamma_{\theta}^{\prime\prime})}{r\gamma_{r}\gamma_{\theta}}\\ R_{2}^{s}=\frac{\gamma_{r}^{3}+r\gamma_{\theta}\gamma_{r}^{\prime}(r\gamma_{\theta})^{\prime}-\gamma_{r}(\gamma_{\theta}^{2}+r^{2}\gamma_{\theta}^{\prime 2}+r\gamma_{\theta}(4\gamma_{\theta}^{\prime}+r\gamma_{\theta}^{\prime\prime}))}{\gamma_{r}^{3}}\end{array}\right.$$

$$\begin{array}[]{cc}\displaystyle\lambda_{r/\theta}(r,t)=\bar{\lambda}_{r/\theta}(r)+\int_{\mathring{\vartheta}(r)}^{t}\dot{\lambda}_{r/\theta}(r,s)\,ds\\ \displaystyle s_{r/\theta}(r,t)=\overline{s}_{r/\theta}(r)+\int_{\mathring{\vartheta}(r)}^{t}\dot{s}_{r/\theta}(r,s)\,ds\end{array}$$

$$\left\{\begin{array}[]{lll}\bar{\lambda}_{r}(r)=\bar{\chi}^{\prime}(r)-\mathring{\vartheta}^{\prime}(r)\,\overline{\dot{\chi}}(r)\\ \bar{\lambda}_{\theta}(r)=\bar{\chi}(r)/r.\end{array}\right.$$

$$\left\{\begin{array}[]{llll}\displaystyle\dot{s}_{r}={\mathcal{A}}_{rr}\dot{\lambda}_{r}+{\mathcal{A}}_{r\theta}\dot{\lambda}_{\theta}\\ \displaystyle\dot{s}_{\theta}={\mathcal{A}}_{\theta r}\dot{\lambda}_{r}+{\mathcal{A}}_{\theta\theta}\dot{\lambda}_{\theta}.\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\displaystyle\gamma_{r}=\sqrt{\lambda_{r}^{2}-2s_{r}/(E\lambda_{r})+4s_{\theta}\nu/(E\lambda_{\theta}))}\\ \displaystyle\gamma_{\theta}=\sqrt{\lambda_{\theta}^{2}+2s_{r}\nu/(E\lambda_{r})+2s_{\theta}(\nu-1)/(E\lambda_{\theta})}.\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\displaystyle\gamma_{r}=\lambda_{r}\exp\left(\left(2\nu s_{\theta}\lambda_{\theta}-s_{r}\lambda_{r}\right)/E\right)\\ \displaystyle\gamma_{\theta}=\lambda_{\theta}\exp\left(\left(s_{\theta}\lambda_{\theta}(\nu-1)+\nu s_{r}\lambda_{r}\right)/E\right).\end{array}\right.$$

$$\left\{\begin{array}[]{llll}\displaystyle\tilde{\mathcal{A}}_{rr}=E\lambda_{r}^{2}(1-\nu)/((1+\nu)(1-2\nu))+s_{r}/\lambda_{r}\\ \displaystyle\tilde{\mathcal{A}}_{r\theta}=\tilde{\mathcal{A}}_{\theta r}=E\nu\lambda_{r}\lambda_{\theta}/((1+\nu)(1-2\nu))\\ \displaystyle\tilde{\mathcal{A}}_{\theta\theta}=E\lambda_{\theta}^{2}(1-\nu)/((1+\nu)(1-2\nu))+s_{\theta}/\lambda_{\theta}\end{array}\right.$$

$$\left\{\begin{array}[]{llll}\displaystyle\hat{\mathcal{A}}_{rr}=E(1-\nu)/(\lambda_{r}^{2}(1+\nu)(1-2\nu))-s_{r}/\lambda_{r}\\ \displaystyle\hat{\mathcal{A}}_{r\theta}=\hat{\mathcal{A}}_{\theta r}=E\nu/(\lambda_{r}\lambda_{\theta}(1+\nu)(1-2\nu))\\ \displaystyle\hat{\mathcal{A}}_{\theta\theta}=E(1-\nu)/(\lambda_{\theta}^{2}(1+\nu)(1-2\nu))-s_{\theta}/\lambda_{\theta}.\end{array}\right.$$

$$\frac{\partial {\dot{s}}_r}{\partial r}+\frac{n-1}{r}({\dot{s}}_r-{\dot{s}}_{\theta })=0.$$

$${\dot{s}}_r=\mathring{D}\left(\frac{d{\overline{s}}_r}{dr}+\frac{n-1}{r}({\overline{s}}_r-{\overline{s}}_{\theta })\right).$$

$${\dot{s}}_r-\frac{1}{r}{\overline{\sigma }}_{\theta }\overline{\dot{\chi }}=D\left({\overline{\sigma }}_r^{\prime }{\bar{\lambda }}_{\theta }+\frac{1}{r}({\overline{\sigma }}_r-{\overline{\sigma }}_{\theta })\right)$$

$${\dot{s}}_r-\frac{2}{r}{\overline{\sigma }}_{\theta }{\bar{\lambda }}_{\theta }\overline{\dot{\chi }}=D\left({\overline{\sigma }}_r^{\prime }{\bar{\lambda }}_{\theta }+\frac{2}{r}({\overline{\sigma }}_r-{\overline{\sigma }}_{\theta })\right){\bar{\lambda }}_{\theta }$$

$${{\mathring{\epsilon }}_{\theta }}^{\prime \prime }+\frac{2{\mathring{\epsilon }}_{\theta }^{\prime }-{\mathring{\epsilon }}_r^{\prime }}{r}=0$$

$${\sigma }_{\theta }^{\infty }=\frac{E({\mathring{\epsilon }}_r-{\mathring{\epsilon }}_{\theta })}{2(1-{\nu }^2)}$$

$${\sigma }_{\theta }^a=-\frac{E{\mathring{\epsilon }}_{\theta }}{1-{\nu }^2}.$$

Full text

I Introduction

A variety of natural and artificial processes rely on active mass deposition on the surface of a solid body. The associated class of phenomena is quite broad including such diverse processes as growth of plants Archer ; Dumais2001 , cell motility Dafalias ; johnpre14 , construction of retaining walls GoodmanSlopes ; Labuz , formation of planets Kadish2005 , crystallisation from solution Wildeman ; Fink ; Schwerdtfeger and 3D printing Ge ; LindNature16 . Surface growth is understood in this context as a continuous addition of new layers of mass on the external boundary of a solid. From the perspective of elasticity theory, the fundamental interest here is in the fact that the accreeted mass points arrive with their own reference state, which must emerge as an outcome of the manufacturing process.

Despite considerable interest of compatible surface growth Skalak ; Epstein ; DiCarloSurf ; MauginCiarletta ; Ganghoffer2010 ; Tomassetti ; MoultonGoriely ; Jenkins , here we focus on the case when the manufactured reference state is * incompatible* (non-Euclidean) in the sense that it cannot be realized in 3D without generating residual stresses. The underlying “geometric frustration” kondo49 ; Efrati ; RazPNAS ; Sharon , which is ultimately shaped by the deposition process, may be beneficial (as in growing plants GorielyPlants ) or detrimental (as in civil engineering structures Hossain ), as it was already exemplified in the early attempts to understand incompatible surface growth motivated by the necessity to explain the built-up of “growth stresses” in trees Martley , to optimize the concrete pouring Arutyunyan and to improve the quality of industrial winding Southwell .

The first systematic theoretical study of the effects of incompatibility in surface growth was conducted by the Russian school Rashba ; Trincher ; Naumov ; Manz14 ; Drozdov98 ; Lychev17 with the largely parallel development and subsequent extension of the theory in the West Goodman ; Fletcher ; Zabaras ; Gambarotta ; SozioYavari ; GuptaTrees ; GanghofferGoda ; Papadopoulos ; Goriely19 ; Bulging . In particular, these studies have raised an awareness of the presence of a “historical element” in the incompatible surface growth problems, which implies that accumulated inelastic strains keep a detailed memory of the deposition process. In defiance of these efforts, however, the relation between the accretion protocol and the ensuing state of geometric frustration remains poorly understood. The problem is that this relation is inherently nonlocal in both, space and time due to the long range character of elastic interactions, and the unavoidable coupling between the incremental adjustments of the elastic state to the advances of the accretion front.

This problem was addressed in the two recent papers ZTprl ; ZTMaugin where we developed a linearized theory of incompatible surface growth focused on the deposition-protocol-dependence of the resulting state of residual stress. The analysis in the linear case turned out to be relatively simple, primarily, because the geometry could be decoupled from elasticity. In the present paper we generalize this theory by taking into account both geometrical (e.g. finite rotations) and physical (e.g. finite stretches) nonlinearities of an elastic solid.

Geometric nonlinearities are particularly relevant when the kinematics of the growth process cannot be linearized, as, for instance, in the case of actin polymerization against a solid wall johnpre14 ; Misbah2016 , where the identification of the reference and the actual states is incompatible with the very presence of the incoming mass flux. Physical nonlinearities are crucial, for instance, when due to extreme dependence of the elastic moduli on pre-stress MacKintosh , growth induced inhomogeneity can lead to dramatic spatial heterogeneity of the elastic response. The use of a nonlinear theory becomes, of course, imperative when one deals with soft solids like biological tissues or synthetic gels Nardinocchi ; Teresi .

In this paper we formulate a general theory of large-strain incompatible surface growth. To illustrate the specific effects of nonlinearity we present a systematic analysis of the case of spherical symmetry, where partial differential equations reduce to ordinary differential equations and at least some of the computations can be performed analytically. Several case studies were chosen to show the details of the implementation of the general approach.

Our first example deals with winding protocols producing disclination-type incompatibility. The second example shows how such singular incompatibility reveals itself if a 2D disk is allowed to relax into the 3D space. The third example illustrates mass accretion under prescribed external pressure, emphasizing extreme sensitivity of the embedded incompatibility to small changes in the deposition protocol. In the last two examples we consider outward and inward surface growth against rigid constraints. Such problems are out of reach for geometrically linearized theory, moreover, we show that the account of physical nonlinearities in this context allows one to predict the emergence of growth-induced material instabilities.

The paper is organized as follows. In Sec. II we develop the concept of inelastic surface growth and discuss the geometrical meaning of incompatibility. In Sec. III we compare various deposition protocols and formulate the corresponding conditions on the growth surface. In Sec.IV we present the incremental formulation of the accretion problem in the general nonlinear setting. In Sec. V we use the nonlinear theory to obtain its linearized version developed in ZTprl ; ZTMaugin . The specialization of the general theory for the case of spherical symmetry is discussed in Section Sec. VI. The case studies, illustrating various effects of physical and geometrical nonlinearity, are presented in Sec. VII. Finally, Sec. VIII contains our conclusions.

II Preliminaries

Consider a 3D body $\mathring{\mathcal{B}}$ equipped with an arbitrary set of (Lagrangian) coordinates x and a metric tensor $\textbf{G}(\textbf{\emph{x}})$ allowing one to measure distances in $\mathring{\mathcal{B}}$. Suppose that G is compatible in the sense that there exists the deformation (embedding) of $\mathring{\mathcal{B}}$ in the Eucledian space $\mathcal{B}=\textbf{\emph{y}}(\mathring{\mathcal{B}})$ such that $\textbf{G}=\nabla\textbf{\emph{y}}^{\textsf{T}}\nabla\textbf{\emph{y}}$ where $\textbf{\emph{y}}=\textbf{\emph{y}}(\textbf{\emph{x}})$. Suppose also that $\mathring{\mathcal{B}}$ is equipped with a reference metric $\mathring{\textbf{G}}(\textbf{\emph{x}})$ which does not have to be compatible and may carry information about the local configuration of defects.

The elastic response is calibrated by the distance between G and $\mathring{\textbf{G}}$, so we can introduce an elastic energy density $e(\textbf{G},\mathring{\textbf{G}})\geq 0$, such that $e(\mathring{\textbf{G}},\mathring{\textbf{G}})=0$. In equilibrium the (Piola-Kirchhoff) stress tensor $\textbf{S}=\partial e/\partial\nabla\textbf{\emph{y}}$ must satisfy the conditions

[TABLE]

where $\mathring{\textbf{\emph{n}}}$ is the normal to $\partial\mathring{\mathcal{B}}$. We assumed here that the body $\mathring{\mathcal{B}}$ is loaded by surface tractions $\textbf{\emph{s}}(\textbf{\emph{x}})$ and body forces $\textbf{\emph{f}}(\textbf{\emph{x}})$.

Suppose now that the body $\mathring{\mathcal{B}}_{t}$ is growing due to addition of points on the evolving boundary $\partial\mathring{\mathcal{B}}_{t}$, see Fig.1. This process can be described by prescribing a function $\mathring{\vartheta}(\textbf{\emph{x}})=t$ that defines the placement of the referential surface at a time instant $t$. The normal to this boundary is $\mathring{\textbf{\emph{n}}}=||\nabla\mathring{\vartheta}||^{-1}\nabla\mathring{\vartheta}$ and its normal velocity is $\mathring{D}=||\nabla\mathring{\vartheta}||^{-1}$. Given that inertial terms can be neglected, the deformation of the growing body $\textbf{\emph{y}}(\textbf{\emph{x}},t)$ must satisfy a one-parametric family of equilibrium equations

[TABLE]

where $\textbf{\emph{s}}(\textbf{\emph{x}},t)$ and $\textbf{\emph{f}}(\textbf{\emph{x}},t)$ can be now time dependent. If the reference metric $\mathring{\textbf{G}}(\textbf{\emph{x}},t)$ is known and the energy density $e$ satisfies suitable conditions Ball , the system (2) can be solved at each instant.

To simplify the problem we can assume that the reference metric $\mathring{\textbf{G}}(\textbf{\emph{x}})$ is time independent, which means that it is fixed for each material particle at the stage of deposition and is not evolving afterwords. The challenge is to link physically realizable “controls” on the evolving surface with particular distributions of $\mathring{\textbf{G}}(\textbf{\emph{x}})$ in the grown body Vanel .

The reference metric $\mathring{\textbf{G}}$ carries information about how the body is assembled, which is contained in the Ricci tensor R$(\textbf{\emph{x}},\textbf{\emph{y}})=\tsum_{i}\langle\mathcal{R}(\textbf{\emph{x}},\textbf{\emph{e}}_{i})\textbf{\emph{y}},\textbf{\emph{e}}_{i}\rangle$, a contraction of the Riemann-Christoffel tensor $\mathcal{R}(\textbf{\emph{x}},\textbf{\emph{y}})\textbf{\emph{z}}=\nabla_{\textbf{\emph{x}}}\nabla_{\textbf{\emph{y}}}\textbf{\emph{z}}-\nabla_{\textbf{\emph{y}}}\nabla_{\textbf{\emph{x}}}\textbf{\emph{z}}-\nabla_{[\textbf{\emph{x}},\textbf{\emph{y}}]}\textbf{\emph{z}}$, where $\textbf{\emph{x}},\textbf{\emph{y}},\textbf{\emph{z}}$ are arbitrary vectors, $\left\{\textbf{\emph{e}}_{i}\right\}$ is an ortho-normal basis in $\mathbb{R}^{3}$ and $\nabla$ the Levi-Civita connection induced by $\mathring{\textbf{G}}$ (see Kuhnel for details).

If $\textbf{R}(\mathring{\textbf{G}})\equiv 0$ there exists a map $\mathring{\textbf{\emph{g}}}$ such that $\mathring{\textbf{G}}=\nabla\mathring{\textbf{\emph{g}}}^{\textsf{T}}\nabla\mathring{\textbf{\emph{g}}}$ Ciarlet and the body is unstressed in the absence of loading ($\textbf{\emph{s}}=\textbf{\emph{0}},\textbf{\emph{f}}=\textbf{\emph{0}}$). Indeed, the equilibrium equations in this case

[TABLE]

have a homogeneous solution $\textbf{S}(\nabla\textbf{\emph{y}},\nabla\mathring{\textbf{\emph{g}}}^{\textsf{T}}\nabla\mathring{\textbf{\emph{g}}})=\textbf{0}$ with $\textbf{\emph{y}}\equiv\mathring{\textbf{\emph{g}}}$. Such growth is compatible and the map $\mathring{\textbf{\emph{g}}}$ only affects the stress-free shape of the body $\mathring{\textbf{\emph{g}}}(\mathring{\mathcal{B}}_{t})$, see for instance Nardinocchi .

If, instead, $\textbf{R}(\mathring{\textbf{G}})\neq 0$, then it is not possible to find a smooth deformation y such that $\textbf{S}(\nabla\textbf{\emph{y}},\mathring{\textbf{G}})=\textbf{0}$ and the unloaded body is prestressed. Since the tensor R$(\mathring{\textbf{G}})$ identically satisfies a set of 3 scalar differential constraints (Bianchi identities), $\text{Div}(\textbf{R}-S\,\mathring{\textbf{G}}/2)=0$, where $S=\tsum_{i,j}\langle\mathcal{R}(\textbf{\emph{e}}_{i},\textbf{\emph{e}}_{j})\textbf{\emph{e}}_{j},\textbf{\emph{e}}_{i}\rangle$, only three independent components of $\textbf{R}(\mathring{\textbf{G}})$ characterize the distribution of growth related defects Davini ; RazPNAS ; Hsu68 ; Teresi . The presence of such defects is a sign of geometric frustration, which can be viewed as an embedded “information”, revealed through residual stresses.

III Deposition protocols

Suppose that the presence of an active agent on the growing surface can be modeled by nonstandard boundary conditions defining the “growth protocol”. One possibility is to directly prescribe the functions $\mathring{\textbf{G}}(\textbf{\emph{x}})$ on the growth surface. Three of its “incompatible” components will then describe the accumulated defects and acquired residual stresses, while the other three “compatible” components will characterize the shape of the body at the end of the deposition process; the latter is relevant in many biological problems Nardinocchi and in engineering problems related to residual actuation Danescu2018 .

However, the case when the functions $\mathring{\textbf{G}}(\textbf{\emph{x}},t)$ can be directly controlled on the growth surface is not very realistic. Instead, the reference state is “manufactured” during deposition, as a result of the physical actions that one can control through the deposition machinery. Typical examples of natural processes where $\mathring{\textbf{G}}(\textbf{\emph{x}})$ is not fixed a priori, include layered manufacturing Gambarotta , 3D printing Ge and biological accretion Skalak ; Skalak82 .

III.1 Stress control

In view of the mechanical nature of an “active agent” involved in the execution of a deposition protocol, it is natural to formulate the corresponding boundary conditions in terms of the deposition stress $\overline{\textbf{S}}(\textbf{\emph{x}})=\textbf{S}(\textbf{\emph{x}},\mathring{\vartheta}(\textbf{\emph{x}}))$. Here and in what follows, the overline bar will indicate the restriction of a function to $\partial\mathring{\mathcal{B}}_{t}$ NoteBar .

The tensor function $\overline{\textbf{S}}(\textbf{\emph{x}})$ is a priori partially constrained by the applied tractions $\overline{\textbf{S}}\mathring{\textbf{\emph{n}}}=\textbf{\emph{s}}$ and by the fact that the angular momentum is balanced, $\overline{\textbf{S}}(\overline{\nabla\textbf{\emph{y}}})^{\textsf{T}}=(\overline{\nabla\textbf{\emph{y}}})\overline{\textbf{S}}\,^{\textsf{T}}$. This implies that out of nine scalar components of $\overline{\textbf{S}}$, six are essentially known. The remaining three can be still used to at least partially constrain the distribution of embedded inhomogeneity, and will constitute our first set of active controls.

To be more specific, consider the decomposition of the deposition stress

[TABLE]

where the “passive” contribution $\textbf{S}_{p}(\textbf{\emph{x}})$ is defined by the applied surface tractions $\textbf{S}_{p}=\textbf{\emph{s}}\otimes\mathring{\textbf{\emph{n}}}+\mathring{\textbf{\emph{n}}}\otimes\textbf{\emph{s}}-(\textbf{\emph{s}}\cdot\mathring{\textbf{\emph{n}}})\mathring{\textbf{\emph{n}}}\otimes\mathring{\textbf{\emph{n}}}$ and the “active” component is restricted to the surface $\textbf{S}_{a}\mathring{\textbf{\emph{n}}}=0$, see Fig.2. Fixing the three independent components of $\textbf{S}_{a}(\textbf{\emph{x}})$ can be then considered as a stress-control part of the deposition protocol.

As we show in our examples, the stress control can be performed not only through the Piola stress $\textbf{S}(\textbf{\emph{x}})$ but also through the Cauchy stress $\boldsymbol{\sigma}(\textbf{\emph{y}})$. The two are related through $\textbf{S}=\boldsymbol{\sigma}\text{cof}(\nabla\textbf{\emph{y}})$, where $[\text{cof}(\textbf{A})]^{\textsf{T}}=\textbf{A}^{-1}\det\textbf{A}$, and the difference between them is an effect of geometrical nonlinearity.

III.2 Kinematic control

The other three controls can be, for instance, of kinematic origin and involve constraints on the current position of the growing surface, which amounts to prescribing the functions $\bar{\textbf{\emph{y}}}(\textbf{\emph{x}})=\textbf{\emph{y}}(\textbf{\emph{x}},\mathring{\vartheta}(\textbf{\emph{x}}))$. Such constraints can be imposed passively as, in the case of growth against a fixed wall, or actively, as in the case of growth under an oscillating piston.

Another important case when the controls may involve the field $\bar{\textbf{\emph{y}}}(\textbf{\emph{x}})$ is that of solidification, where the built-up of pressure at the growth interface depends on the current position of the solidification front Fletcher . In those cases one can use the fact that the map $\bar{\textbf{\emph{y}}}$ completely characterizes the areal stretch $\alpha=||\text{cof}\,(\nabla\bar{\textbf{\emph{y}}})\mathring{\textbf{\emph{n}}}||$ and the normal $\textbf{\emph{n}}=\alpha^{-1}\text{cof}\,(\nabla\bar{\textbf{\emph{y}}})\mathring{\textbf{\emph{n}}}$ of the current surface $\partial\mathcal{B}_{t}$.

Note that the control of the actual position of the growth surface $\textbf{\emph{y}}(\textbf{\emph{x}},\mathring{\vartheta}(\textbf{\emph{x}}))$ implies the knowledge of the reference position of this surface defined by the function $\mathring{\vartheta}(\textbf{\emph{x}})$. The knowledge of the latter is equivalent to fixing the rate of mass delivery at the referential growth surface $\dot{m}=\mathring{\varrho}\mathring{D}$, where $\mathring{\varrho}$ is the referential density.

We can also express the arriving mass flux in terms of the actual parameters $\dot{m}=\bar{\varrho}\alpha(D-\bar{\dot{\textbf{\emph{y}}}}\cdot\textbf{\emph{n}})$, where $\bar{\varrho}=\mathring{\varrho}(\det\overline{\nabla\textbf{\emph{y}}})^{-1}$ is the actual density and $D=\mathring{D}(\nabla\bar{\textbf{\emph{y}}})\mathring{\textbf{\emph{n}}}\cdot\textbf{\emph{n}}$ is the Eulerian velocity of the growing surface. Using the (commutation) relation

[TABLE]

we may rewrite the expression for $D$ in the form $\bar{J}\alpha^{-1}\mathring{D}+\bar{\dot{\textbf{\emph{y}}}}\cdot\mathring{\textbf{\emph{n}}}$ where $\bar{J}=\det\overline{\nabla\textbf{\emph{y}}}$. In the special case when the position of the growing surface is fixed in the actual configuration $(D=0)$ the Eulerian growth rate takes a simpler form $\dot{m}/\alpha=-\bar{\varrho}(\bar{\dot{\textbf{\emph{y}}}}\cdot\mathring{\textbf{\emph{n}}})$.

Finally we remark that if the functions $\bar{\textbf{\emph{y}}}(\textbf{\emph{x}})$ are not constrained directly on the growth surface, we can always use the freedom of choosing the Lagrangian coordinates of the arriving material particles to set $\bar{\textbf{\emph{y}}}(\textbf{\emph{x}})=\textbf{\emph{x}}$. This will correspond to the choice of the reference configuration coinciding at the moment of deposition with the instantaneous actual configuration and will partially specify the reference metric. Such “pseudo-linearization” of the deformation geometry would be of course impossible in many realistic cases, for instance, in problems involving growth against rigid obstacles.

IV Solving the growth problem

We now show that if the function $\mathring{\vartheta}(\textbf{\emph{x}})$ is given, the six controls $(\bar{\textbf{\emph{y}}}(\textbf{\emph{x}}),{\textbf{S}}_{a}(\textbf{\emph{x}}))$, combined with the classical mechanical conditions on the growth surface, are sufficient to ensure a particular reference metric distribution $\mathring{\textbf{G}}(\textbf{\emph{x}})$ in the grown body.

Suppose, for simplicity, that the constitutive relation $\textbf{S}=\textbf{S}(\nabla\textbf{\emph{y}},\mathring{\textbf{G}})$ can be inverted, allowing one to express functionally the (time independent) reference metric in terms of the (time dependent) current values of stress and strain, to obtain $\mathring{\textbf{G}}=\mathring{\textbf{G}}(\nabla\textbf{\emph{y}},\textbf{S})$. For instance, this will be the case when the elastic energy is quadratic in strain, which in the case of isotropy implies that

[TABLE]

Here the elastic strain $\boldsymbol{\epsilon}_{e}$ can be chosen in different ways, for instance, in the “linear” form $\boldsymbol{\epsilon}_{e}^{\text{lin}}=(\textbf{G}-\mathring{\textbf{G}})/2$ defining Kirchhoff-S.Venant material GreenNaghdi ; Efrati , or in the “logarithmic” form $\boldsymbol{\epsilon}_{e}^{\text{log}}=\left(\log(\textbf{G}\mathring{\textbf{G}}^{-1})\right)/2$ defining Hencky material Bruhns ; Auricchio ; Anand . Note that in the case of linear strain $\boldsymbol{\epsilon}_{e}^{\text{lin}}$, a finite energy is needed to squeeze a volume element into a point, while this energy is infinite when we use the logarithmic strain $\boldsymbol{\epsilon}_{e}^{\text{log}}$.

The obtained inversion $\textbf{S}=\textbf{S}(\nabla\textbf{\emph{y}},\mathring{\textbf{G}})$, however, does not solve our problem directly, because the projection of such relation on the growth surface involves $\overline{\nabla\textbf{\emph{y}}}$ rather than $\nabla\bar{\textbf{\emph{y}}}$. Indeed, the (commutation) relation (5) indicates an essential coupling of the equilibrium problems (2) at different time instants. The functions $\overline{\dot{\textbf{\emph{y}}}}(\textbf{\emph{x}})$ will then depend on the whole deposition protocol, making the relation between $\mathring{\textbf{G}}(\textbf{\emph{x}})$ and the controls $(\bar{\textbf{\emph{y}}}(\textbf{\emph{x}}),\textbf{S}_{a}(\textbf{\emph{x}}))$ fundamentally nonlocal. The simplest way to deal with such path dependence is to consider an incremental reformulation of the problem.

We start with the straightforward representations

[TABLE]

In (7) we can use the relation $\dot{\textbf{S}}=\mathbfcal{A}\,\dot{\nabla}\textbf{\emph{y}}$ where $\mathbfcal{A}=\partial_{\nabla\textbf{\emph{y}}}^{2}e(\textbf{G},\mathring{\textbf{G}})$ is the tangential elasticity tensor, which depends on the reference metric. We can also write $\text{Div}\textbf{S}=\text{Div}\bar{\textbf{S}}+\int_{\mathring{\vartheta}(\textbf{\emph{x}})}^{t}\text{Div}\dot{\textbf{S}}\,ds-\dot{\textbf{S}}\nabla\mathring{\vartheta}$, and since body forces can be written as

[TABLE]

the equilibrium condition $\text{Div}\textbf{S}+\textbf{\emph{f}}=\textbf{\emph{0}}$ transforms into the incremental field equation for the stress increment, together with a boundary condition on incremental tractions Trincher ; ZTprl

[TABLE]

The boundary condition (9)2 shows that if the new mass points are deposited on $\partial\mathring{\mathcal{B}}_{t}$ in a state of mechanical unbalance ($\text{Div}\overline{\textbf{S}}+\overline{\textbf{\emph{f}}}\neq 0$), they are instantly re-equilibrated by activating the incremental elastic displacements $\dot{\textbf{\emph{y}}}(\textbf{\emph{x}},t)$.

The knowledge of the functions $\nabla\textbf{\emph{y}}(\textbf{\emph{x}},t)$ and $\textbf{S}(\textbf{\emph{x}},t)$ at each time step allows one to use the relation $\mathring{\textbf{G}}=\mathring{\textbf{G}}(\overline{\nabla\textbf{\emph{y}}},\overline{\textbf{S}})$ to compute the reference metric of the newly adhered layers. This information is needed to update the tangential elasticity tensor $\mathbfcal{A}$, which can be then used to solve the next incremental equilibrium problem.

V Linearization

It is instructive to compare the nonlinear theory presented above with its linear counterpart ZTprl ; ZTMaugin . Suppose that the displacement $\textbf{\emph{u}}(\textbf{\emph{x}},t)=\textbf{\emph{y}}(\textbf{\emph{x}},t)-\textbf{\emph{x}}$ is small in the usual sense of linear elasticity. Then the analysis of the current metric $\textbf{G}(\textbf{\emph{x}},t)$ can be replaced by the study of the linearized strain $\boldsymbol{\epsilon}(\textbf{\emph{x}},t)=\text{sym}\nabla\textbf{\emph{u}}(\textbf{\emph{x}},t)$, and instead of the reference metric $\mathring{\textbf{G}}(\textbf{\emph{x}})$ one can consider the linear reference strain $\mathring{\boldsymbol{\epsilon}}(\textbf{\emph{x}})$. In this approximation, the elastic constitutive relation takes the simple form $\textbf{S}=\mathbb{C}\boldsymbol{\epsilon}_{e}$, where $\boldsymbol{\epsilon}_{e}=\boldsymbol{\epsilon}-\mathring{\boldsymbol{\epsilon}}$ is the elastic strain. The Hookean elasticity tensor $\mathbb{C}$ will be taken positive definite and space/time independent.

While in the nonlinear theory the function $\bar{\textbf{\emph{y}}}(\textbf{\emph{x}})$ can be appreciably different from x, in the linear setting we must have $\bar{\textbf{\emph{y}}}(\textbf{\emph{x}})\simeq\textbf{\emph{x}}$. This implies that out of our six controls, three are automatically fixed and we have only the surface stress ${\textbf{S}}_{a}(\textbf{\emph{x}})$ to work with. This is, of course, consistent with the fact that in the linear theory only the three independent components of the incompatibility of the reference strain $\mathring{\boldsymbol{\eta}}=\text{Curl}\text{Curl}\mathring{\boldsymbol{\epsilon}}$ affect the solution of the elasticity problem VanGoethem . A fundamental advantage of the linear formulation is that the relation between the incompatibility in the newly accreted points and the solution of the incremental problem can be made explicit ZTprl

[TABLE]

where the tensor $\overline{\boldsymbol{\eta}}=-\text{Curl}\text{Curl}(\mathbb{C}^{-1}\overline{\textbf{S}})$ can be viewed as a measure of the incompatibility in the arriving material. If instead the reference incompatibility $\mathring{\boldsymbol{\eta}}$ is known, the stress distribution at each moment of time $t$ can be found directly by solving the boundary value problem

[TABLE]

Such a formulation is, of course, not possible in the nonlinear setting, where the knowledge of the full tensor $\mathring{\textbf{G}}(\textbf{\emph{x}})$ is required to find the stress distribution $\textbf{S}(\textbf{\emph{x}},t)$ at each stage of mass accretion.

In the incremental problem the equilibrium equations (9) remain the same with $\dot{\textbf{S}}(\textbf{\emph{x}},t)=\mathbb{C}\nabla\dot{\textbf{\emph{u}}}(\textbf{\emph{x}},t)$. An important difference, though, is that the incremental displacement field is defined on a fixed configuration of the body, which is controlled at each moment of time by the function $\mathring{\vartheta}(\textbf{\emph{x}})$. Otherwise the procedure of solving a one parametric family of the incremental problems remain the same as in the nonlinear case with prescription of the tensor $\overline{\textbf{S}}_{p}(\textbf{\emph{x}})$ which defines the tractions $\textbf{\emph{s}}(\textbf{\emph{x}})$ and fixes the three extra components of surface stress $\overline{\textbf{S}}_{a}(\textbf{\emph{x}})$ controlling the residual incompatibility $\mathring{\boldsymbol{\eta}}(\textbf{\emph{x}})$.

While the linearized theory preserves some of the complexity of the full nonlinear formulation, it under-represents several important effects. For instance, in the linear theory the reference and actual configurations are identified, which makes it impossible to deal with problems involving confined growth in fixed domains. Another problem is that in the linearized theory we assume that the elasticity moduli are fixed, while in the nonlinear theory, the whole deposition history is encoded both in $\mathring{\textbf{G}}$, and in the tensor field $\mathbfcal{A}(\textbf{\emph{x}},t)$. The implied elastic inhomogeneity, accumulated during accretion, is particularly relevant for biological materials where the hardening nonlinearity may be extremely strong MacKintosh . Strong inhomogeneity of elastic properties can also become the origin of material instabilities FosdickEdelstein , that can, of course, serve by themselves as potential targets in the design of surface growth protocols.

VI Radial symmetry

To illustrate the above general theory, consider now in more detail the accretion process under the condition of radial symmetry (both in two and in three dimensions).

Assume that the evolving reference configuration is a part of a ball $\mathring{\mathcal{B}}_{t}=\left\{\textbf{\emph{x}}\,|\,\psi_{0}\leq r\leq\psi_{t}\right\}$ in $\textbf{R}^{n}$ where $r=||\textbf{\emph{x}}||$. The function $\psi_{t}=\psi_{0}+\mathring{D}\,t$ will define the position of the growth surface at time $t$, where $\psi_{0}$ is its initial position. For simplicity, the referential growth velocity $\mathring{D}$ will be taken as constant. We can then write $\mathring{\vartheta}(r)=(r-\psi_{0})/\mathring{D}$, where $\mathring{D}=1/\mathring{\vartheta}^{\prime}$.

Denote the unit vector pointing in radial direction by $\textbf{\emph{e}}=\textbf{\emph{x}}/r$. Then the radially symmetric deformation can be written as $\textbf{\emph{y}}(\textbf{\emph{x}},t)=\chi(r,t)\textbf{\emph{e}}(\textbf{\emph{x}})$. The ensuing deformation gradient is

[TABLE]

where $\lambda_{r}(r,t)=\partial_{r}\chi(r,t)$ and $\lambda_{\theta}(r,t)=\chi(r,t)/r$ are the stretches in radial and azimuthal directions, and $\mathbb{P}=\textbf{I}-{\textbf{\emph{e}}}\otimes{\textbf{\emph{e}}}$ is the projection tensor. In view of radial symmetry, the most general form of the time-independent reference metric is

[TABLE]

For the case of a disk ($n=2$), the Ricci tensor of $\mathring{\textbf{G}}$ reduces to the Ricci scalar (twice the Gaussian curvature $K$) Kuhnel

[TABLE]

The requirement $S^{d}=0$ produces then a single differential constraint

[TABLE]

on the functions $(\gamma_{r},\gamma_{\theta})$ where $k$ is an arbitrary constant. For a 3D sphere ($n=3$) the Ricci tensor reduces to two independent components,

[TABLE]

Such Ricci tensor vanishes if and only if $\gamma_{r}=\gamma_{\theta}=\text{const}$, which corresponds to a homogeneous conformal dilatation of the reference sphere.

To determine the unknown functions $\gamma_{r}(r)$ and $\gamma_{\theta}(r)$ we need to prescribe two auxiliary conditions on the growth surface. In the case of radial symmetry (for both $n$= 2, 3) the vector function $\bar{\textbf{\emph{y}}}$ reduces to a scalar function $\bar{\chi}(r)$ and, in view of the representation $\bar{\textbf{S}}=\bar{s}_{r}(r)\textbf{\emph{e}}\otimes\textbf{\emph{e}}+\bar{s}_{\theta}(r)\mathbb{P}$, the active stress $\textbf{S}_{a}$ reduces to its hoop component $\bar{s}_{\theta}(r)$. The additional conditions prescribing the growth protocol, may then take the form of restrictions imposed on the functions $\bar{\chi}(r)$ and $\bar{s}_{\theta}(r)$ with the other component of the deposition stress $\bar{s}_{r}(r)$ fixed passively. Note that in the problem of inward accretion against a fixed wall we have $\mathring{\textbf{\emph{n}}}=-\textbf{\emph{e}}$, so that the mass flux $\dot{m}=\bar{\varrho}\alpha\bar{\dot{\chi}}$; for outward accretion, $\mathring{\textbf{\emph{n}}}=\textbf{\emph{e}}$ and $\dot{m}=-\bar{\varrho}\alpha\bar{\dot{\chi}}$.

To formulate the incremental problem when the functions $\bar{\chi}(r)$ and $\bar{s}_{\theta}(r)$ are known, we first need to specialize (7) for the case of spherical symmetry

[TABLE]

where the deposition stretches can be written as

[TABLE]

To compute the rates $\dot{\lambda}_{r}=\partial\dot{\chi}$, $\dot{\lambda}_{\theta}=\dot{\chi}/r$ we need to introduce the incremental moduli

[TABLE]

These moduli can be expressed in terms of the current values of stretches $\displaystyle\lambda_{r/\theta}(r,t)$ and stresses $\displaystyle s_{r/\theta}(r,t)$ which, as we have mentioned before, depend on the whole accretion history.

Since $\mathring{\textbf{G}}(\textbf{\emph{x}})$ is time-independent, we can time differentiate the inverse constitutive equation $\mathring{\textbf{G}}(\nabla\textbf{\emph{y}},\textbf{S})$ to obtain the relation $\partial_{\nabla\textbf{\emph{y}}}\mathring{\textbf{G}}[\nabla\dot{\textbf{\emph{y}}}]+\partial_{\textbf{S}}\mathring{\textbf{G}}[\dot{\textbf{S}}]=0$. Then $\dot{\textbf{S}}=\tilde{\mathbfcal{A}}\nabla\dot{\textbf{\emph{y}}}$, where $\tilde{\mathbfcal{A}}=-(\partial_{\nabla\textbf{\emph{y}}}\mathring{\textbf{G}})^{-1}(\partial_{\textbf{S}}\mathring{\textbf{G}})$. In the case of a sphere, we can first use the elastic model (6) with linear strains $\boldsymbol{\epsilon}_{e}^{\text{lin}}$ to invert the elastic constitutive relation

[TABLE]

If, instead, we use the logarithmic strain $\boldsymbol{\epsilon}_{e}^{\text{log}}$, we obtain

[TABLE]

Since the inelastic strains are time-independent, by time differentiation of (20) and (21) we obtain the incremental constitutive equations. In the case when we use the linear strain measure $\boldsymbol{\epsilon}_{e}^{\text{lin}}$, the tangential elasticity tensor $\tilde{\mathbfcal{A}}$ has components

[TABLE]

while if we use the logarithmic strain $\boldsymbol{\epsilon}_{e}^{\text{log}}$,

[TABLE]

We now have all we need to formulate the sequence of incremental equilibrium problems. Specialization of the equilibrium equations (9)1 to the case of radial symmetry gives

[TABLE]

The boundary condition (9)2 on the growth surface takes the form

[TABLE]

Note that often it may be more convenient to specify the deposition protocol in Eulerian rather than Lagrangian coordinates. In this case we can still formulate the incremental problem in terms of Piola-Kirchhoff components of stress, however, the boundary conditions (25) would have to be modified. Thus, if we can control at deposition the components of Cauchy rather than Piola-Kirchhoff stress, the incremental boundary condition describing the radially symmetric growth will take the form

[TABLE]

in the case of a disk ($n=2$), and

[TABLE]

in the case of a sphere ($n=3$). Here we used the Eulerian velocity of the growing surface $D=\bar{\chi}^{\prime}\mathring{D}$, and we denoted $\bar{\sigma}_{r}^{\prime}(y)=d\sigma_{r}(y)/dy.$ In the case of a disk ($n=2$) the relation between Cauchy and Piola-Kirchhoff stresses reads $s_{r}=\sigma_{r}\lambda_{\theta}$ and $s_{\theta}=\sigma_{\theta}\lambda_{r}$, while in the case of a sphere ($n=3$), $s_{r}=\sigma_{r}\lambda_{\theta}^{2}$ and $s_{\theta}=\sigma_{\theta}\lambda_{r}\lambda_{\theta}$. Note that the boundary conditions (27), (26) couple the controlled stress $\overline{\sigma}_{r},\overline{\sigma}_{\theta}$ and the controlled shape $\bar{\chi}(r)$, which is a purely nonlinear effect.

VII Case studies

To illustrate different aspects of the developed general theory, we now elaborate few explicit case studies, which are all relevant for applications. Our main goal here is to highlight the effects of physical and geometrical nonlinearities, and to illustrate path dependence of the surface growth process. Our examples show that the deposition protocol is effectively remembered by the grown body through the accumulated incompatibility, the inhomogeneity of the incremental moduli and the shape of the body liberated from the constraints. In all examples we assume that the function $\mathring{\vartheta}(\textbf{\emph{x}})$ is prescribed, leaving aside the problem of the feedback received by the growth mechanism from the current state of stress.

VII.1 Winding of 2D disks

Consider a process of winding of an infinitesimally thin tape on a rigid mandrel. Suppose that the growth process takes place in 2D so that the emerging hollow disk is constrained to remain in plane. Assume that the winding is accomplished by pulling the adhering tape with a controlled tangential force. If the adhering layers have thickness $h$ and the tangential force is $F$, we assume that in the double limit $h\to 0$, $F\to 0$ the hoop stress $\bar{\sigma}_{\theta}=F/h$ remains finite, see Fig.3a. The growing disk $\mathcal{D}=\{\textbf{\emph{x}}\,|\,a\leq||\textbf{\emph{r}}||\leq\psi\}$ (with $\psi\leq r_{e}$) can be parameterized by the referential radius $\psi$ of the external circle where the tape deposition takes place.

We begin with the simplest inverse problem: given that the acquired reference metric $\mathring{\textbf{G}}=\text{diag}(\gamma_{r},\gamma_{\theta})$ is homogeneous, find the corresponding deposition protocol. Somewhat counter-intuitively, this metric is nontrivial, being in general locally compatible and globally incompatible. Indeed, in view of (15) the Gaussian curvature $K$ is equal to zero everywhere except for the origin, where the polar coordinate system is singular and the curvature has a delta function type singularity as long as $k=\gamma_{r}/\gamma_{\theta}\neq 1$. Physically such target configuration corresponds to a disk with removed ($k>1$) or inserted ($k<1$) wedge (with opening angle $\Delta\theta=2\pi(k^{-1}-1)$) RazPRE ; RazARMA .

Observe next that the prescription of a reference metric implicitly constrains the choice of the Lagrangian coordinates for the reference configuration of the growing body. If we assume, in addition, that the normal tractions on the growth surface are absent ($s_{r}(\psi,\psi)=0$), we can solve the sequence of incremental problems and compute both, the current hoop stress $\bar{s}_{\theta}(\psi)=s_{\theta}(\psi,\psi)$ and the current radius of the disk $\bar{\chi}(\psi)$. With this information at hand, we can express the Cauchy hoop stress $\sigma_{\theta}=s_{\theta}/\lambda_{r}$ in terms of the current position of the grow surface, to obtain the function $\sigma_{\theta}(\chi)$ controlling the winding process. Our computational results are illustrated in Fig.3b for the case of 2D Hencky material with logarithmic strain $\boldsymbol{\epsilon}_{e}^{\text{log}}$.

Note that to obtain discs with positive embedded Gaussian curvature ($\gamma_{r}>\gamma_{\theta}$) the controlled traction force must be positive, whereas to embed negative curvature ($\gamma_{r}<\gamma_{\theta}$), it must be negative. This is consistent with the fact that disks with positive curvature are obtained by removing a wedge: indeed, by wrapping through pulling we deposit less mass per length than in the case of zero force wrapping. Instead, pushing produces an opposite result and leads to disks with negative curvature.

To understand these numerical results more deeply, we can compare them with the analytically transparent linearized theory. To this end we set $\gamma_{r/\theta}=1+\mathring{\varepsilon}_{r/\theta}$ and assume that the inelastic strains $\mathring{\varepsilon}_{r/\theta}$ are small. The linearization of the scalar Ricci curvature (14) gives the compatibility condition in the bulk,

[TABLE]

which can be also rewritten as $(r\mathring{\varphi})^{\prime}=0$, where $\mathring{\varphi}={\mathring{\varepsilon}_{\theta}}^{\prime}+({\mathring{\varepsilon}}_{\theta}-{\mathring{\varepsilon}}_{r})/r$. Therefore $\mathring{\varphi}=c/r$ where $c$ an arbitrary constant and if $c\neq 0$ there will be a curvature singularity at the origin (signifying the presence of an inserted or removed wedge).

We reiterate that even though the inelastic strains $\mathring{\varepsilon}_{r/\theta}$ are homogeneous, the residual stresses in this case will be different from zero. Indeed, if $u$ is a radial displacement, the total strains $\varepsilon_{r}=u^{\prime}$ and $\varepsilon_{\theta}=u/r$ can be additively decomposed as $\varepsilon_{r/\theta}=e_{r/\theta}+\mathring{\varepsilon}_{r/\theta}$, where the elastic strains $e_{r/\theta}$ are constitutively related to stresses $\varepsilon_{r/\theta}=(\sigma_{r/\theta}-\nu\sigma_{\theta/r})/E$, where $\nu$ is the Poisson ratio and $E$ is the Young modulus. Since in equilibrium $\sigma_{r}^{\prime}+(\sigma_{r}-\sigma_{\theta})/r=0$, the residual stresses in a traction-free disk can be determined by solving the equation $(r^{3}\sigma_{r}^{\prime})^{\prime}=-Er^{2}\mathring{\varphi}$ with the right hand side $\mathring{\varphi}=({\mathring{\varepsilon}}_{\theta}-{\mathring{\varepsilon}}_{r})/r$. If ${\mathring{\varepsilon}}_{\theta}\neq{\mathring{\varepsilon}}_{r}$ the solution will be obviously nonzero even for a hollow disk with zero boundary conditions $\sigma_{r}(a)=\sigma_{r}(r_{e})=0$.

We can now use the linearized theory to address the peculiar behavior of the winding tension $\sigma_{\theta}(\psi)$ in the nonlinear theory, see Fig.3a. Solving the linear equilibrium problem with boundary conditions $u(a,\psi)=0$ and $\sigma_{r}(\psi,\psi)=0$, we obtain $\sigma_{\theta}(\psi)=E((\mathring{\varepsilon}_{r}-\mathring{\varepsilon}_{\theta})(\psi^{2}-2a^{2}(1-2\nu)\log(a/\psi))-a^{2}(\mathring{\varepsilon}_{r}+\mathring{\varepsilon}_{\theta}(3-4\nu)))/(2(a^{2}(1-2\nu)+\psi^{2})(1-\nu^{2}))$. In particular, at $\psi\rightarrow\infty$ we have

[TABLE]

whereas in the limit $\psi\rightarrow a$ we obtain

[TABLE]

These analytical results confirm the general trends observed in the numerical solution of the nonlinear problem and show to what extent the linearized theory overestimates the necessary winding tension, see Fig.3.

VII.2 Isometric embedding of 2D disks into 3D

To illustrate the growth-induced incompatibility produced by winding, we now allow the grown disk to detach from its (imaginary) 2D substrate and take a relaxed shape in 3D. Such an isometric embedding $\boldsymbol{\psi}:\mathbb{R}^{2}\mapsto\mathbb{R}^{3}$ leads to the full relaxation of elastic energy and must therefore satisfy the system of equations $\mathring{\textbf{G}}=\nabla\boldsymbol{\psi}^{\textsf{T}}\nabla\boldsymbol{\psi}$, where $(\gamma_{r},\gamma_{\theta})$ is the homogeneous diagonal reference metric acquired in the process of controlled growth discussed above. We remark that such relaxed shapes may serve by themselves as the target of controlled surface manufacturing Sharon ; Danescu .

Note first that the relaxed shape of the disk in 3D depends crucially on the sign of its Gaussian curvature which is proportional to the effective wedge opening $\Delta\theta=2\pi(\gamma_{\theta}/\gamma_{r}-1)$. As we show below, for $\gamma_{r}>\gamma_{\theta}$ the isometric embedding takes the form a circular cone, whereas for $\gamma_{\theta}>\gamma_{r}$ it is known as an excess cone (or *anti *cone) BenAmarConicalDefects ; Warner .

We begin with the trivial case $\gamma_{r}=\gamma_{\theta}$, when the isometric embedding is confined to 2D and $\boldsymbol{\psi}(\textbf{\emph{x}})=g(r)\textbf{\emph{e}}(\textbf{\emph{x}}):\mathbb{R}^{2}\mapsto\mathbb{R}^{2}$. Here e is the unit vector in the radial direction, and $g(r)=r\gamma_{\theta}=r\gamma_{r}$. Clearly, in this case, the relaxed configuration of the disk is just another flat disk.

If $\gamma_{r}>\gamma_{\theta}$, the Gaussian curvature is positive and, even if the functions $\gamma_{r}$ and $\gamma_{\theta}$ are $r$ dependent, we can write a universal embedding $\boldsymbol{\psi}(\textbf{\emph{x}})=\rho(r)\textbf{\emph{e}}(\textbf{\emph{x}})+\zeta(r)\textbf{\emph{k}}$ where the unit vector k is perpendicular to the reference plane. The unknown functions $\rho(r)$ and $\zeta(r)$ must satisfy $\gamma_{r}^{2}={\rho^{\prime 2}+\zeta^{\prime 2}}$ and $\gamma_{\theta}=\rho/r$, so that for homogeneous metric $\rho=r\gamma_{\theta}$ and $\zeta=(\gamma_{r}^{2}-\gamma_{\theta})^{2})^{1/2}$; the latter expression shows that $\zeta$ is real only if $\gamma_{r}>\gamma_{\theta}$. The corresponding relaxed shape is shown in Fig. 4a.

When $\gamma_{r}<\gamma_{\theta}$ the embedded curvature is negative and the reconstruction of the 3D relaxed surface is not straightforward BenAmarConicalDefects ; Sharon ; HanHong . An approximate construction can be built based on the auxiliary surface $\boldsymbol{\psi}(\textbf{\emph{x}})=r\rho(\theta)\textbf{\emph{e}}(\textbf{\emph{x}})+r\zeta(\theta)\textbf{\emph{k}}$ (see Warner ), whose induced metric is also diagonal $(\tilde{\gamma}_{r}^{2}(\theta),\tilde{\gamma}_{\theta}^{2}(\theta))$ with

[TABLE]

Generically, this metric is not homogeneous and is therefore, strictly speaking, incompatible with the homogeneous reference metric generated by the winding process. However, we can demand that one of its components is homogeneous $\tilde{\gamma}_{r}(\theta)=\gamma_{r}$ and the other one meets the original anzats in average so that $\frac{1}{2\pi}\int_{0}^{2\pi}\tilde{\gamma}_{\theta}(\theta)\,d\theta=\gamma_{\theta}.$ To meet this constraint, whose physical meaning is that the image of a constrained reference circle has a prescribed perimeter in the relaxed state Warner , we can, for instance, assume that $\zeta(\theta)=A\,\text{sin}(m\theta)$, where $m>1$ is an integer.

The integral condition provides a link between $A$ and $m$ and for the case of homogeneous $\gamma_{r},\gamma_{\theta}$, the amplitude of the oscillation decreases with increasing wave number $m$, see Fig. 5. The resulting smooth embedding is, of course, only an approximation of the relaxed surface which can be expected to be rather rough Nash . However, lower order modes, like the one shown in Fig. 4b, will be energetically favored if, for instance, small bending stiffness is taken into consideration.

VII.3 Outward growth under controlled pressure

Our next example concerns outward accretion of a hollow sphere under controlled pressure, which is physically relevant in the context of pressurized freezing Fletcher . The goal of this example is to illustrate the possibility of extreme path dependence of the accretion process.

We assume that the reference configuration is a hollow sphere $\mathring{\mathcal{B}}_{t}=\left\{\textbf{\emph{x}}\,|\,A\leq||\textbf{\emph{x}}||\leq\psi_{t}\leq B\right\}$, where $A$ and $B$ are the initial and final radii. The current configuration is $\mathcal{B}_{t}=\left\{\textbf{\emph{y}}\,|\,a\leq||\textbf{\emph{y}}||\leq\bar{\chi}(\psi_{t})\right\}$, which shows that the internal, non-growing surface is forced to remain on a rigid bead of radius $a$.

We further assume that the growth protocol is characterized by the following conditions on the accretion surface

[TABLE]

The second condition (32) is just the simplest assumption that the reference configuration for the arriving material particles coincides with the actual configuration. Instead, on the non-growing part of the boundary we impose the hard constraint $\chi(A,t)=A$.

For determinacy, we assume that

[TABLE]

where $\beta$ and $m$ are prescribed constants. All members of this “family” begin ($\psi_{0}\equiv A$) and end ($\psi_{T}\equiv B$) at the same level of pressure while exhibiting super or sub-linear growth depending on whether $m>1$ or $m<1$.

Our numerical simulations for the Hencky material (with strain $\boldsymbol{\epsilon}_{e}^{\text{log}}$) show appreciable protocol (history) sensitivity, which can be seen in the final stress distribution, see Fig. 6(a) and in the final inhomogeneity of tangential elastic moduli, see Fig. 6(c). The most striking effect concerns the final distribution of incompatibility. As we see in Fig. 6(b), a seemingly insignificant change in the control of pressure during deposition, can lead to either divergence of the two components of curvature at the inner boundary of the body (for $m<1$) or their convergence to zero ( for $m>1$). In other words our computations show that the embedded curvature lacks continuity with respect to the parameter $m$. The origin of this dramatic effect can be explained already in the framework of a more transparent linearized theory.

Note first that in 3D linearized theory the metric is compatible if two conditions are satisfied simultaneously: $(r\mathring{\varphi})^{\prime}=0$ and $(r^{3}\mathring{\varphi})^{\prime}=0$ where again $\mathring{\varphi}=\mathring{\varepsilon}_{\theta}^{\prime}+(\mathring{\varepsilon}_{\theta}-\mathring{\varepsilon}_{r})/r$. The solution of these equations is $\mathring{\varphi}=0$ and since the equilibrium problem for the stress distribution in the unloaded 3D body is

[TABLE]

hollow spheres behave quite differently than hollow disks, where there can be residual stresses even if the regular part of curvature vanishes.

The 3D linearized problem of surface accretion under a general pressure distribution $p(z)$ can be solved explicitly. We now assume that $a=A$. If we again replace time by the radius $\psi$ of the growing surface, the solution to the problem (24) with incremental boundary conditions $\dot{u}(a,\psi)=0$ and $\dot{\sigma}_{r}(\psi,\psi)=-p^{\prime}(\psi)$ is $\dot{u}(r,\psi)=c_{1}(\psi)r+c_{2}(\psi)/r^{2}$ with $c_{1}(\psi)=(1+\nu)(2\nu-1)\psi^{3}p^{\prime}(\psi)/(E(2(1-2\nu)a^{3}+(1+\nu)\psi^{3}))$ and $c_{2}(\psi)=-a^{3}c_{1}(\psi)$. The resulting stress distribution takes the form

[TABLE]

In particular, these expressions show that the case $a=0$ is trivial in the sense that the stress in the growing sphere is necessarily hydrostatic.

The distribution of residual stresses in a hollow sphere of internal radius $r_{i}$ and external radius $r_{e}$ is given by (34). The expression for $\mathring{\varphi}(r)$ can be found by solving the linearized problem of surface accretion with $r_{i}=a$:

[TABLE]

This expression confirms that no incompatibility can result from a growth process unless the sphere is hollow ($a\neq 0$) and the applied pressure varies along the deposition process. If we substitute (33) into (36) and expand the result in small $\delta=(r-a)/a$, we obtain $\mathring{\varphi}(\delta)\sim\delta^{m-1}$ which shows that the curvature either diverges or tends to zero in proximity of the inner radius depending on whether $m<1$ or $m>1$. This is exactly what we found numerically in the case of nonlinear Hencky material (with strain $\boldsymbol{\epsilon}_{e}^{\text{log}}$), see Fig.6b.

VII.4 Inward growth on a rigid bead

Consider now accretion on a rigid spherical surface, mimicking the process of inwards actin polymerization stimulated biochemically on a specially treated spherical bead Tomassetti ; Misbah2016 ; SykesProst , see Fig.7b. A similar process with cylindrical symmetry would be the growth of a tree where new mass is deposited between the existing trunk and the bark Archer ; GuptaTrees . We use this example to highlight the essential role of geometrical nonlinearity because the linearized setting of such a problem is not even meaningful.

Indeed, in our problem the reference configuration $\mathring{\mathcal{B}}_{t}=\left\{\textbf{\emph{x}}\,|\,\psi_{t}\leq||\textbf{\emph{x}}||\leq\psi_{0}\right\}$ and the actual configuration $\mathcal{B}_{t}=\left\{\textbf{\emph{y}}\,|\,\chi(\psi_{t},t)\leq||\textbf{\emph{y}}||\leq\chi(\psi_{0},t)\right\}$ are fundamentally different (see Fig.7a), because the current radius of the growing surface is fixed $\chi(\psi_{t},t)=a$ while its reference radius $\psi_{t}$ continuously evolves. As a result, the reference and actual domains cannot coincide, which is the natural starting assumption of any geometrically linear elasticity theory.

In physical terms, we assume that the growth process starts at a surface of a rigid spherical bead and that the arriving material is being continuously “squeezed” between the emerging grown body and the original rigid surface (the actual mechanism of mass delivery is obviously disregarded). While the growing body is expanding away from the bead and its external radius $\chi(\psi_{0},t)$ is an increasing function of time, the reference radius of the growth surface $\psi_{t}$ is a decreasing function of time.

To determine both, the current state of stress and the two unknown components $(\gamma_{r}(r),\gamma_{\theta}(r))$ of the reference metric, we impose three boundary conditions on the growth surface

[TABLE]

The first of this conditions states that the attachment stress is hydrostatic with time dependent pressure which is controlled externally. In deviation from the previous example, we assume that the pressure control is not direct but is rather an outcome of the control of the Eulerian velocity of the arriving material $\bar{\dot{\chi}}(\psi_{t})=\dot{m}/(\bar{\varrho}\alpha)$. Essentially this means the control of a volumetric inflow rate which in the context actin polymerization appear to be more realistic than the full control of the attachment stress Tomassetti ; Misbah2016 . On the exterior surface of the growing body we assume the no-tractions condition, $s_{r}(\psi_{0},t)=0$.

The succession of incremental problems can be solved numerically and we illustrate the final stress distribution in Fig.8b. Note that the deviatoric stress is maximal at the external (non-growing) surface of the body which, in principle, should lead to surface instabilities DeSimone . The parameters were tuned to match the numerical results obtained for a growing network of actin rods, biased to polymerize on contact with a spherical bead Misbah2016 . The pressure build up at the rigid surface, the final distribution of the components of the Ricci tensor and the residual inhomogeneity of the elastic moduli, are shown in Fig.9.

VII.5 Growth induced material instabilities

The aim of our last example is to show that incompatible mass accretion can lead to material instabilities. Since the letter are not related directly to geometrical effects, this example stresses the importance of physical nonlinearity in surface growth problems.

Suppose that the reference configuration is given by $\mathring{\mathcal{B}}_{t}=\left\{\textbf{\emph{x}}\,|\,\psi_{0}\leq||\textbf{\emph{x}}||\leq\psi_{t}\right\}$ where now it is the outward surface of a hollow sphere that is growing. The current configuration is $\mathcal{B}_{t}=\left\{\textbf{\emph{y}}\,|\,\chi(\psi_{0},t)\leq||\textbf{\emph{y}}||\leq\bar{\chi}(\psi_{t})\right\}$, and we assume again that the growing surface is constrained by a rigid wall $\chi(\psi_{t},t)=a$. In other wards, we assume that the mass is continuously “squeezed” between the existing surface of the body and the rigid spherical cavity, see Fig.11c. Once again, the actual and the reference configurations can not coincide at any time which makes the geometrical linearization of the problem meaningless.

On the growth surface we maintain the same conditions as in the previous example

[TABLE]

and again, instead of prescribing pressure directly, we assume that the growth is controlled through the Eulerian velocity of the arriving material $\bar{\dot{\chi}}(\psi_{t})=-\dot{m}/(\bar{\varrho}\alpha)$. Also, similarly to the case of inward growth, we assume zero tractions condition on the internal surface, $s_{r}(\psi_{0},t)=0$.

The ensuing sequence of incremental problems can be solved numerically and our focus now will be on the material stability of the grown body. More specifically, we consider the case of linear strain tensor $\boldsymbol{\epsilon}_{e}^{\text{lin}}$ (Kirchhoff-St.Venant material) and check whether the strong ellipticity condition KnopsPayne may be violated during accretion.

We recall that the strong ellipticity condition requires that $\mathbfcal{A}(\textbf{\emph{m}}\otimes\mathring{\textbf{\emph{k}}})\cdot(\textbf{\emph{m}}\otimes\mathring{\textbf{\emph{k}}})>0$ where $\mathbfcal{A}=\partial\textbf{S}/\partial\textbf{F}$ is the instantaneous elasticity tensor, and $\textbf{\emph{m}},\mathring{\textbf{\emph{k}}}$ are unit vectors. To ensure this condition it is sufficient to check that the acoustic tensor $\textbf{Q}(\textbf{\emph{m}})$, whose action on a vector $\mathring{\textbf{\emph{a}}}$ is defined by $\textbf{Q}(\textbf{\emph{m}})\mathring{\textbf{\emph{a}}}=\mathbfcal{A}(\mathring{\textbf{\emph{a}}}\otimes\textbf{\emph{m}})\textbf{\emph{m}}$, is positive definite GurtinLinearTh . To locate the domain of potential material instability we must then search for the possibility that the determinant $\det\textbf{Q}(\textbf{\emph{m}})$ is negative for some m Ogden .

The results of a typical numerical experiment are shown in Fig.11, where $\psi_{t}$ is the increasing radius of the growth surface. Before a critical radius of the growth surface is reached, strong ellipticity is ensured on the whole domain. At the critical threshold the strong ellipticity is first lost at the non-growing (traction free) surface. As the growth surface advances, the domain where the strong ellipticity is lost, spreads behind getting progressively closer to the growth surface. It should be mentioned, however, that while strong ellipticity may be formally lost, our incremental equilibrium problems are constrained to maintain radial symmetry. For this reason, our numerical tests do not show instability even after the unconstrained problem becomes ill-posed.

VIII Conclusions

A large-strain theory of incompatible surface growth is formulated here in terms of the new type of boundary conditions. These conditions are prescribed on a free surface which can move along the reference coordinates. The novelty of this formulation is in the account of geometric and constitutive nonlinearities of an elastic body, which are both necessary to adequately capture finite rotations and finite stretches. The theory is presented in a general 3D setting, with the main focus on path dependence of the incompatible accretion process, and its goal is to develop the means to control the acquired incompatibility by tailoring the deposition protocol.

We showed that if both, stresses and displacements, can be independently controlled at the growing surface, one can ensure point-wise manipulation of the resulting reference metric which keeps a detailed memory of the deposition process. The relation between the embedded Ricci tensor of the reference metric and the time dependence of the control parameters was shown to be nonlinear and nonlocal. It was shown to emerge from a solution of a one parametric family of incremental equilibrium problems for linear elastic bodies with evolving geometry and varying elastic inhomogeneity.

The problem simplifies in the case of small deformations when geometry decouples from elasticity, which becomes linear with elastic moduli both homogeneous and time independent. We showed that even in this case the incremental elastic problem is characterized by nonstandard boundary conditions. However, the acquired incompatibility can be then expressed analytically in terms of the incremental strain rates which simplifies considerably the dependence of the deposited reference metric on the surface controls. The main shortcoming of the linearized theory is its inability to deal with kinematically confined growth, and to account for the possibility of elastic instabilities induced by the growth process.

We illustrated the general theory by a series of examples emphasizing the role of finite strains in the surface growth of soft solids. Our examples highlight the inherent path dependence of the incompatible growth and emphasize the effects of geometrical and physical nonlinearities in ensuring particular outcomes of the physically realizable growth protocols. Through these examples we showed that geometrical frustration developing during deposition can be indeed fine-tuned and that a particular behavior of a system in physiological or industrial conditions can be engineered by embedding into the material a particular incompatibility. The proposed general theoretical framework allowing one to handle such information-rich solids can be used to guide bio-mimetic design and additive manufacturing.

The fact that one can engineer the incompatibility in a manufactured solid body and regulate in this way the associated distribution of residual stresses is of crucial importance for the understanding of biological processes where solids are typically functionally pre-stressed. Our analysis also highlights the possibility to manufacture non-Euclidean elastic solids artificially, not only with a tailored pre-stress, but also with a particular stiffness distribution. Such solids can be designed to undergo specific elastic instabilities and to exhibit specific patterning in the technologically relevant conditions Gao12 ; Sharon11 ; Bangham ; Roman ; Burgueno ; 2Dprinting .

IX Acknowledgments

The authors thank V.Balbi, D.Barilari, R.Kupferman, E.Sharon, B.Shoikhet and G.Tomassetti for helpful discussions. G.Z. acknowledges the PMMH-ESPCI for hospitality during 2018, the sabbatical support from NUI-Galway and the support of the Italian National Group of Mathematical Physics (INDAM-GNFM). L.T. was supported by the French government under the Grant No. ANR-10-IDEX-0001-02 PSL.