Accretion and ablation in deformable solids with an Eulerian description: examples using the method of characteristics

Abstract

Accretion and ablation, i.e., the addition and removal of mass at the surface, are important in a wide range of physical processes, including solidification, growth of biological tissues, environmental processes, and additive manufacturing. The description of accretion requires the addition of new continuum particles to the body, and is therefore challenging for standard continuum formulations for solids that require a reference configuration. Recent work has proposed an Eulerian approach to this problem, enabling side-stepping of the issue of constructing the reference configuration. However, this raises the complementary challenge of determining the stress response of the solid, which typically requires the deformation gradient that is not immediately available in the Eulerian formulation. To resolve this, the approach introduced the elastic deformation as an additional kinematic descriptor of the added material, and its evolution has been shown to be governed by a transport equation. In this work, the method of characteristics is applied to solve concrete simplified problems motivated by biomechanics and manufacturing. Specifically, (1) for a problem with both ablation and accretion in a fixed domain and (2) for a problem with a time-varying domain, the closed-form solution is obtained in the Eulerian framework using the method of characteristics without explicit construction of the reference configuration.

Keywords

Accretion kinematics method of characteristics Eulerian formulation growth

1. Introduction

Accretion or surface growth involves the addition of mass at the boundary of a body. In a continuum setting, this requires the introduction of new material particles to the body [1], and poses unusual challenges. Surface growth occurs in a range of settings, including in biological tissue [2], planetary formation [3, 4], solidification and casting processes [5], etching processes in silicon wafers [6], additive manufacturing processes [7], the construction of masonry structures [8], and cell motility through the assembly and disassembly of actin networks [9].

The addition and removal of material points in surface growth causes difficulties in the usual procedure of defining a fixed set of points as the reference configuration. Dealing with this in linear elasticity is somewhat easier [3, 8, 10 –12]. However, in non-linear descriptions of surface growth, the definition of the reference configuration is more subtle. Pioneering work by Skalak, Hoger, and coworkers provided a systematic approach to study the kinematics of surface growth [1, 13]. However, these studies did not consider deformation, and assumed that the body was rigid. Recent work has extended these ideas to account for deformation and stress [14 –23]. These studies used a Lagrangian formulation, i.e., the kinematic description is roughly based on reference particles in an evolving reference configuration.

To avoid the difficulties of dealing with an evolving reference configuration, our previous work [24] developed an Eulerian description of surface growth. This has the advantage that we do not need to calculate the reference configuration explicitly, e.g. [25 –27]. Conversely, in a solid, the stress response requires knowledge of the deformation gradient $F$ , and the latter is not straightforward to calculate in an Eulerian setting. Prior work, notably [28 –30], has developed approaches to find the deformation gradient in an Eulerian setting; specifically, Liu and Walkington [28] use an evolution equation for $F$ that does not require explicit calculation of the reference configuration.

While the Eulerian approach works well for the standard setting of a fixed set of material particles without growth, there is the challenge of defining the reference state and kinematic information of the added material particles during surface growth. Specifically, the specification of $F$ for the added material must ensure that kinematic compatibility is not violated. To avoid dealing with the issue of kinematic compatibility, we introduced two additional kinematic descriptors: the relaxed zero-stress deformation $F_{relax}$ and the elastic deformation $F_{e}$ . These are related to the deformation gradient through the relation $F_{e} = {FF}_{relax}^{- 1}$ . They enable us to eliminate $F$ —as well as $F_{relax}$ , using the fact that it is transported with the material particles—and thereby formulate a model that is posed in terms of $F_{e}$ as the sole kinematic descriptor for the stress response.

The evolution of $F_{e}$ is governed by a transport equation, and transparently handles both accretion and ablation: accretion corresponds to inflow, and requires the specification of boundary conditions, while ablation corresponds to outflow and requires no boundary conditions. Further, $F_{e}$ has no requirements of kinematic compatibility, enabling a simple and direct approach to the specification of boundary conditions. The boundary conditions are specified in terms of $F_{e}$ for the added particles, which is directly and solely related to their given stress state. We can think of $F_{e}$ and $F_{relax}$ as providing a mechanism for the added material to define its own reference configuration that satisfies kinematic compatibility.

We apply this formulation to examine two simplified concrete problems of recent interest. In the first problem, motivated by the assembly and disassembly of actin microfilaments during polymerization and depolymerization on a spherical bead [31, 32], and building on the model of this phenomenon from [17], we examine simultaneous accretion at the fixed inner surface and ablation at the time-dependent outer surface of a hollow sphere. In the second problem, motivated by the wound roll process of additive manufacturing, which has been applied to processes ranging from the packaging of flexible films in paper and magnetic tape industries [33] to tissue engineered blood vessels [34], we study the deformation and stress distribution when stress-free particles are deposited at the outer layer of a hollow cylinder with a time-dependent pressure applied at the inner surface. We solve these problems in closed form with the method of characteristics, using the fact that the evolution of $F_{e}$ is governed by a transport equation.

1.1. Organization

Section 2 summarizes the approach proposed in [24]. We start with a discussion on balance laws and jump or boundary conditions, using the view of surface growth as a localized source on the growing surface. We then turn to the kinematics of the deformation in the Eulerian setting with accretion. Sections 3 and 4 solve specific examples of (1) simultaneous accretion at the fixed inner surface and ablation at the time-dependent outer surface of a hollow sphere and (2) accretion on the outer surface of a long hollow cylinder with a time-dependent pressure applied at the inner surface. We use the method of characteristics to solve the transport equation that governs the evolution of the elastic deformation.

2. Physical model and governing equations

2.1. Notation

The motion is denoted $x = χ (X, t)$ , where $X$ is the location in the reference configuration of a point that is currently at the spatial location $x$ . Also, the inverse of the motion $χ - 1 (x, t)$ maps the current spatial location of a particle at $x$ to the location $X$ in the reference configuration. The deformation gradient is $F = \partial χ (X, t) / \partial X$ . For short, we write $F (x, t) = F (χ - 1 (X, t), t)$ .

Because most of our work is in the Eulerian setting, all differential operators (e.g. $\nabla, div$ ) imply derivatives with respect to $x$ , except where explicitly stated. Similarly, the argument of various field quantities will implicitly be $x$ , except where explicitly stated.

2.2. Balance laws

Following [24], we model surface growth through the introduction of sources that are localized on the growing surface. The source terms can be understood as a coarse-grained approach to treat the complex process of growth without considering the fine details of the processes in the ambient environment outside the growing body; the growth is defined only in terms of net mass and linear momentum transfer. To define the source terms, we specify the rate of mass addition per unit area, $M$ , and the vectorial rate of linear momentum addition per unit area, $P$ . We can then infer the velocity (or specific momentum) of the added material at the instant of attachment, $υ_{a} : = P / M$ ; and the mass density, $ρ : = M / (υ_{a} \cdot \hat{n})$ . Accretion is modeled by $M > 0$ and ablation by $M < 0$ . We assume that the added particles do not carry angular momentum, such as that due to individual particle spins, i.e., there is no source of angular momentum owing to growth.

The source terms appear in the jump conditions that represent the balances of mass and momentum at the growing surface:

Mass : [[ρ (V_{b} - υ) \cdot \hat{n}]] = M,

(1)

Momentum : [[ρ υ ((V_{b} - υ) \cdot \hat{n})]] + [[σ \hat{n}]] = M υ_{a} .

(2)

We have used the jump operator [[·]] to denote the difference between the limits when approaching the surface of discontinuity from either side; $V_{b}$ is the velocity of the surface and $\hat{n}$ is the surface normal.

We make the assumption that the ambient environment outside of the growing body has a negligible effect on the mechanics of the body. This leads to the significant simplification that the domain of the problem now involves only the solid body, and the jump conditions become boundary conditions:

Mass : ρ (V_{b} - υ) \cdot \hat{n} = M \Leftrightarrow V_{b} \cdot \hat{n} = υ \cdot \hat{n} + \frac{M}{ρ},

(3)

Momentum : ρ υ ((V_{b} - υ) \cdot \hat{n}) + σ \hat{n} - t_{b} = M υ_{a} \Leftrightarrow σ \hat{n} = M (υ_{a} - υ) + t_{b},

(4)

where $t_{b}$ is the traction owing to the external forces at the boundary of the growing body. We further assume that the inertial terms $ρ υ ((V_{b} - υ) \cdot \hat{n})$ and $M υ_{a}$ are negligible compared with $σ$ . Then, equation (4) reduces to $σ \hat{n} = t_{b}$ .

The balance laws in the bulk are standard:¹

Mass : \frac{\partial ρ}{\partial t} + div (ρ υ) = 0,

(5)

Momentum : \frac{\partial}{\partial t} (ρ υ) + div (ρ υ \otimes υ) = ρ b + div (σ),

(6)

where $σ$ is the Cauchy stress; $ρ$ is the mass density; $υ$ is the particle velocity; $b$ is the body force; $t_{b}$ is the traction vector owing to the external forces at the boundary of the growing body; and $V_{b}$ is the boundary velocity. From equation (3), the boundary velocity $V_{b}$ has contributions from $M / ρ$ and $υ$ (Figure 1).

Figure 1.

Evolution of a growing body in a small interval of time $Δ t$ . The spatial location of the growing boundary depends on the growth velocity and the continuum particle velocity at the boundary.

2.3. Kinematics

The stress response of the material generally requires knowledge of the deformation gradient

F = \frac{\partial χ (X, t)}{\partial X} = {(\frac{\partial χ - 1 (x, t)}{\partial x})}^{- 1};

we notice that it requires knowledge of the motion or its inverse. However, we adapt the approach from [28] that formulates the evolution of $F$ without requiring us to solve for the deformation map.

We can write the material time derivative of $F$ as

\frac{\partial F (X, t)}{\partial t} = \frac{\partial}{\partial t} \frac{\partial χ (X, t)}{\partial X} = \frac{\partial}{\partial X} \frac{\partial x}{\partial t} = \frac{\partial υ}{\partial X} = \frac{\partial υ}{\partial x} \frac{\partial x}{\partial X} = (\nabla υ) F \Rightarrow \overset{\cdot}{F} = (\nabla υ) F,

(7)

where $\overset{\cdot}{g}$ represents the material time derivative of a quantity $g$ . It follows that the evolution of $F$ satisfies the transport equation:

\frac{\partial F (x, t)}{\partial t} + (υ \cdot \nabla) F = (\nabla υ) F .

(8)

While this equation preserves the kinematic compatibility of $F$ for a fixed set of material particles, the setting of surface growth requires the imposition of boundary conditions—corresponding to the deformation state of the added material—that must also be consistent with kinematic compatibility. Assume that the stress state of the added particles is given as $σ^{*}$ . Denote the stress response function of the added particles by $\hat{σ} (\cdot)$ , and assume for simplicity and without loss of generality that $\hat{σ} (I) = 0$ . We notice immediately that $σ^{*} = \hat{σ} (F)$ ; assuming for simplicity that the stress response is invertible² implies that $F = {\hat{σ}}^{- 1} (σ^{*})$ is completely determined at the time of attachment. Since $σ^{*}$ is specified independently of the state of the body, this could potentially violate kinematic compatibility at the time of attachment.

We therefore introduce two additional kinematic descriptors: $F_{relax}$ , which quantifies the relaxed (i.e., zero stress) shape of the added particles, and $F_{e} : = {FF}_{relax}^{- 1}$ , which quantifies the elastic deformation.³ We highlight that neither $F_{e}$ nor $F_{relax}$ need satisfy any requirements of kinematic compatibility. We next redefine the stress response function to be $\hat{σ} (g) \mapsto \hat{σ} ({gF}_{relax}^{- 1})$ . We have the obvious interpretation of $F_{e}$ as the part of $F$ that causes stress, i.e. $σ = \hat{σ} (F_{e})$ . This corresponds to redefining the reference of the added material. Prior to attachment, the deformation state is irrelevant, and post-attachment deformation is governed by the standard equations of continuum mechanics. For an unattached particle, the reference can be changed in arbitrary ways; in particular, it can be changed to ensure that we satisfy kinematic compatibility. In this regard, we highlight that the stress state $σ^{*}$ of the added material is a well-defined and physically meaningful quantity, whereas the deformation can be chosen for convenience by appropriately defining the reference at the time of attachment.

Substituting $F = F_{e} F_{relax}$ in equation (8), we find:

\frac{\partial F_{e} (x, t)}{\partial t} + (υ \cdot \nabla) F_{e} = (\nabla υ) F_{e} .

(9)

We have assumed that $F_{relax}$ is constant in time for each material particle, i.e.,

\frac{\partial F_{relax}}{\partial t} + (υ \cdot \nabla) F_{relax} = 0 .

Reference [24] discusses the more general case where $F_{relax}$ can evolve in time.

To compute the stress response of the body, we only need to know $F_{e}$ , so we solve only equation (9). This is a transport equation, and requires boundary conditions only at inflow boundaries [35], i.e., where the velocity $υ$ is inwards with respect to the moving boundary. Using the given stress state $σ^{*}$ of the added particles, we specify $F_{e} = {\hat{σ}}^{- 1} (σ^{*})$ at the accretion boundaries. While $F_{e}$ has nine components that must be prescribed at the inflow boundary, we highlight that $F_{e}$ is prescribed by the stress only up to a free rotation. Loosely, the stress provides a prescription of only six components, broadly similar to the approach by Zurlo and Truskinovsky [22, 23] that projects stresses.

No boundary condition on $F_{e}$ are needed for the ablation boundaries [35].

Remark 1. Properties of the evolution equation. The form of equation (9) suggests the method of characteristics. Introducing the parametric variable $l$ , we can write the solution of this partial differential equation (PDE) system through the solution of the ordinary differential equations (ODEs)

\frac{d t}{d l} = 1, t (l = 0) = 0 \Rightarrow t = l

(10)

\frac{d x_{i}}{d l} = v_{i}, x_{i} (l = 0) = c_{i}, i = 1, 2, 3

(11)

\frac{d {(F_{e})}_{ij}}{d l} = \frac{\partial v_{i}}{\partial x_{k}} (F_{e})_{kj}, (F_{e})_{ij} (x_{k} (= c_{k}), l = t (= 0)) = initial or boundary condition i, j, k = 1, 2, 3

(12)

in a fixed Cartesian basis, and we have used index notation with implied summation.

The equations of the characteristic curves are defined by equations (10) and (11):

\hat{x} = ({\hat{x}}_{1} (l, c_{1}, c_{2}, c_{3}), {\hat{x}}_{2} (l, c_{1}, c_{2}, c_{3}), {\hat{x}}_{3} (l, c_{1}, c_{2}, c_{3}), \hat{t} (l))

(13)

and equation (12) governs the behavior of $F_{e}$ along the characteristic curves.

Remark 2. Path lines of the motion. Comparing the transport equation (equation (9)) with the transport equation of $χ - 1$ provided in [29], the characteristic lines of both equations are the same and the characteristic curves of the latter are the path lines of the motion. Therefore, the characteristic curves of equation (9) are the path lines of the motion. We will use this property together with equation (3) in the following examples to determine the shape of the body.

2.4. Summary of the proposed approach

Our approach requires the solution of the balances of mass and momentum, and the transport of the elastic deformation. In summary, we solve the three coupled equations given by

\frac{d}{d t} {\begin{matrix} ρ \\ υ \\ F_{e} \end{matrix}} = {\begin{matrix} - ρ div υ \\ \frac{1}{ρ} div \hat{σ} (F_{e}) \\ (\nabla υ) F_{e} \end{matrix}},

(14)

where $d / d t$ is the material time derivative.

3. Hollow sphere with accretion at the inner and ablation at the outer surface

In this section, we consider a hollow sphere with accretion and ablation. The accretion occurs on the inner fixed boundary and the ablation occurs on the outer boundary simultaneously. This example is motivated by the assembly and disassembly of actin microfilaments during polymerization and depolymerization on a spherical bead [31, 32], and builds on the model of this phenomenon proposed in [17]. While we find the general solution, in the special case that there is no initial body, our solution reduces to that reported in [17], where it was solved using a Lagrangian approach with a four-dimensional manifold as the reference configuration.

For the stress response, we use an incompressible neo-Hookean form:

\hat{σ} (F) = - p I + \frac{\partial W (F)}{\partial F} F^{T},

(15)

where $p$ is the pressure and $W$ is the stored energy density function. We note that the incompressible assumption violates the invertibility assumption of the stress response function. For this particular problem, we do not use this assumption; in the general setting, it is sufficient to have invertibility for the part of the stress derived from $W$ .

A schematic of the problem is given in Figure 2. We denote quantities related to the inner and outer surfaces using the subscripts $0$ and $1$ , respectively. We work in the spherical coordinate system:

{\hat{e}}_{r} = \frac{\partial x}{\partial r}, {\hat{e}}_{θ} = \frac{1}{r} \frac{\partial x}{\partial θ}, {\hat{e}}_{ϕ} = \frac{1}{r \sin θ} \frac{\partial x}{\partial ϕ} .

(16)

Figure 2.

Hollow spherical body with accretion at the inner fixed radius and ablation at the outer time-dependent radius.

3.1. Boundary and initial conditions

The inner radius is fixed and is equal to $r_{0}$ , which means that the boundary velocity ${V_{b}}_{0}$ in equation (3) is equal to zero. Particles with given time-dependent velocity $υ_{a_{0}} = - {\overset{\cdot}{Z}}_{0} (t) {\hat{e}}_{r}$ , $({\overset{\cdot}{Z}}_{0} < 0, Z_{0} (0) = 0)$ approach the body at the inner surface and attach to it. Therefore, $M_{0} = - ρ {\overset{\cdot}{Z}}_{0} (t) > 0$ . Then, using equation (3), it is clear that the normal velocity of the material particles right after attachment at $r = r_{0}$ is equal to ${\overset{\cdot}{Z}}_{0} (t)$ , i.e., $υ = - {\overset{\cdot}{Z}}_{0} (t) {\hat{e}}_{r}$ . This will provide a boundary condition for computing the velocity field later on.

Also, mass is removed from the outer surface at the ablation rate $M_{1} < 0$ . In the treadmilling regime, the boundary velocity is zero and the geometry of the body is constant, consequently $M_{1} = - ρ υ (r = r_{1} (t)) \cdot {\hat{e}}_{r}$ .

The outer boundary is traction-free.

These velocity relations correspond to those given by Tomassetti et al. [17]. However, Tomassetti et al. [17] obtain $M_{0} / ρ = - {\overset{\cdot}{Z}}_{0}$ and $M_{1} / ρ = {\overset{\cdot}{Z}}_{1} (t) < 0$ , using a sophisticated coupled analysis based on diffusion and kinetics; here, we will simply assume them as given.

We next turn to the deformation and stress states of the added particles. Since the added particles are composed of the same material as the body, they have the same stress response. We assume that no tractions are imposed on the inner surface and the attaching particles are stress-free, up to an interior hydrostatic pressure $p_{0}$ that must be determined, i.e., $σ (r = r_{0}) = - p_{0} I$ . From the assumption of incompressibility, we can choose $F_{e} = F = F_{relax} = I$ . The outer boundary is an outflow boundary and requires no boundary condition on $F_{e}$ .

Finally, for the initial conditions, we use $F_{e}$ , which is consistent with the boundary conditions and equilibrium equation.

3.2. Reduced governing equations

As the accretion velocity is radial and there is no shear traction, we assume that the motion is independent of $θ$ and $ϕ$ . The velocity $υ$ then simplifies to

υ = v_{r} (r, t) {\hat{e}}_{r} .

(17)

Based on the definition of the problem, the governing equations (equation (14)), together with the stress response (equation (15)), for a neo-Hookean incompressible body reduce to the following:

div (υ) = \frac{\partial v_{r}}{\partial r} + \frac{2 v_{r}}{r} = 0 on r_{0} < r < r_{1} (t)

(18)

υ \cdot {\hat{e}}_{r} = \frac{M_{0}}{ρ} = - {\overset{\cdot}{Z}}_{0} (t), V_{b_{0}} = 0 at r = r_{0}

(19)

\frac{M_{1}}{ρ} = {\overset{\cdot}{Z}}_{1} (t) at r = r_{1} (t) \Rightarrow V_{b_{1}} \cdot {\hat{e}}_{r} = v_{r} (r = r_{1}, t) + {\overset{\cdot}{Z}}_{1} (t)

(20)

div (σ) = (\frac{\partial σ_{rr}}{\partial r} + \frac{2}{r} (σ_{rr} - σ_{θ θ})) {\hat{e}}_{r} = 0 on r_{0} < r < r_{1} (t)

(21)

t_{b} = 0 at r = r_{1} (t)

(22)

\frac{\partial}{\partial t} [\begin{matrix} F_{e_{rr}} & F_{e_{r θ}} & F_{e_{r ϕ}} \\ F_{e_{θ r}} & F_{e_{θ θ}} & F_{e_{θ ϕ}} \\ F_{e_{ϕ r}} & F_{e_{ϕ θ}} & F_{e_{ϕ ϕ}} \end{matrix}] + v_{r} \frac{\partial}{\partial r} [\begin{matrix} F_{e_{rr}} & F_{e_{r θ}} & F_{e_{r ϕ}} \\ F_{e_{θ r}} & F_{e_{θ θ}} & F_{e_{θ ϕ}} \\ F_{e_{ϕ r}} & F_{e_{ϕ θ}} & F_{e_{ϕ ϕ}} \end{matrix}] = [\begin{matrix} \frac{\partial v_{r}}{\partial r} & 0 & 0 \\ 0 & \frac{v_{r}}{r} & 0 \\ 0 & 0 & \frac{v_{r}}{r} \end{matrix}] [\begin{matrix} F_{e_{rr}} & F_{e_{r θ}} & F_{e_{r ϕ}} \\ F_{e_{θ r}} & F_{e_{θ θ}} & F_{e_{θ ϕ}} \\ F_{e_{ϕ r}} & F_{e_{ϕ θ}} & F_{e_{ϕ ϕ}} \end{matrix}]

(23)

F_{e} (r, t = 0) given, F_{e} (r = r_{0}, t) = I

(24)

\hat{σ} (F_{e}) = - p I + G F_{e} F_{e}^{T}, G is the shear modulus .

(25)

3.3. Solution

Using the boundary condition (equation (19)), the continuity equation (equation (18)) gives us the velocity field:

υ = - {\overset{\cdot}{Z}}_{0} \frac{r_{0}^{2}}{r^{2}} {\hat{e}}_{r} .

(26)

We next notice that we can write equation (23) as

(\frac{\partial}{\partial t} + v_{r} \frac{\partial}{\partial r} - A) F_{e_{ij}} = 0, A = {\begin{matrix} \frac{\partial v_{r}}{\partial r} & if ij = rr, r θ, r ϕ, \\ \frac{v_{r}}{r} & otherwise . \end{matrix}

These differential equations are homogeneous, i.e., there is no source term on the right, and therefore those components with zero initial and boundary conditions remain zero during the growth process. Hence, the off-diagonal components of $F_{e}$ are zero throughout the growth process, assuming that $F_{e}$ is diagonal at the initial time. We can then represent $F_{e}$ as:

F_{e} = F_{e_{rr}} {\hat{e}}_{r} \otimes {\hat{e}}_{r} + F_{e_{θ θ}} {\hat{e}}_{θ} \otimes {\hat{e}}_{θ} + F_{e_{ϕ ϕ}} {\hat{e}}_{ϕ} \otimes {\hat{e}}_{ϕ} .

(27)

The evolution equation (equation (23)) and the initial (assuming $F_{e_{θ θ}} (t = 0) = F_{e_{ϕ ϕ}} (t = 0)$ ) and boundary conditions (equation (24)) are the same for $F_{e_{θ θ}}$ and $F_{e_{ϕ ϕ}}$ , as we would expect from spherical symmetry, so we solve only for $F_{e_{θ θ}}$ . Thus, only two of the nine scalar equations for the components of $F_{e}$ are independent and non-trivial. Finally, after substituting the velocity (equation (26)), we obtain two decoupled equations:

\frac{\partial F_{e_{rr}}}{\partial t} - {\overset{\cdot}{Z}}_{0} \frac{r_{0}^{2}}{r^{2}} \frac{\partial F_{e_{rr}}}{\partial r} = 2 {\overset{\cdot}{Z}}_{0} \frac{r_{0}^{2}}{r^{3}} F_{e_{rr}},

(28)

\frac{\partial F_{e_{θ θ}}}{\partial t} - {\overset{\cdot}{Z}}_{0} \frac{r_{0}^{2}}{r^{2}} \frac{\partial F_{e_{θ θ}}}{\partial r} = - {\overset{\cdot}{Z}}_{0} \frac{r_{0}^{2}}{r^{3}} F_{e_{θ θ}} .

(29)

3.3.1. The characteristic curves

To solve equations (28) and (29), we use the method of characteristics, similar to the discussion in Remark 1, because of the convective nature of the equations. The parametric representation of the characteristic lines for both of these equations is

\frac{d t}{d l} = 1 \Rightarrow \hat{t} (l) = l + C_{0},

(30)

\frac{d r}{d l} = - {\dot{Z}}_{0} \frac{r_{0}^{2}}{r^{2}} \Rightarrow {\hat{r}}^{3} (l) = - 3 Z_{0} (\hat{t}) r_{0}^{2} + C_{1},

(31)

where $l$ is the parametric variable along the characteristic curves.

Assuming that all the characteristic curves start from $t = 0$ , we have, for all curves, that $C_{0} = 0$ and $\hat{t} = l$ , although the information about $F_{e}$ might be provided at the middle of the curve (at the attachment time $τ$ of each particle). However, different characteristic curves have different initial locations $\hat{r} (l = 0)$ , so $C_{1}$ is the quantity that differentiates characteristic curves. For the special case of constant ${\overset{\cdot}{Z}}_{0}$ , different characteristic curves are plotted in Figure 3. The information about $F_{e}$ is carried along these characteristic curves, which are also path lines of the motion.

Figure 3.

Characteristic curves that convect initial and boundary deformation through out the space–time domain. (1) The solid non-bold curves are the characteristic curves; (2) the bold black curve distinguishes the initial body (green region) from the added particles (non-green region), and is also the outer radius when there is no ablation and no initial body; (3) the black dashed curve shows the outer radius in the case of no ablation and the existence of an initial body; (4) the gray shaded area shows the body when there is ablation and no initial body; (5) the gray dot-dashed curve shows the outer radius in the case of both ablation and existence of the initial body.

Instead of distinguishing different characteristic curves by their intersection with the line $t = 0$ , one could alternatively determine them by their attachment time $τ$ to the body. In other words, assuming that all characteristic curves start from $r = r_{0}$ at $l = 0$ , then $C_{1} = r_{0}^{3}$ for all characteristic curves but they are distinguished by parameter $C_{0}$ , which is the attachment time $τ$ that was defined by [13] and [17]. We note that the extrapolation of characteristic curves to the $t = 0$ axis and representing them by the attachment time $τ$ provides a four-dimensional manifold, which is denoted the “placement map” in [17]. In the case that there is no initial body at $t = 0$ and all the particles were added to the body by accretion, these two choices of distinguishing characteristic curves are equivalent from the continuum mechanics perspective. However, if there is a stressed or an unstressed initial body at $t = 0$ with accretion and ablation occurring at its boundaries, the approach of using the attachment time cannot as easily predict the deformation of the body during the growth process. In other words, there is no distinction between the particles that were initially in the body (i.e., with $τ = 0$ ) at different radial locations, so their motion can only be considered as a rigid body motion at any later time, as in [13], or extra effort is needed.

The information about $F_{e}$ is carried along these characteristic curves. It can be seen from Figure 3 that, based on the origin of each curve, they may carry either initial condition or boundary condition information. So the space–time domain splits into two regions: (1) the region influenced by the initial condition (the green region), and (2) the region influenced by the boundary condition (the non-green region). The deformation gradient is affected by the initial condition in the lower region of the black bold curve (information comes from the line $t = 0$ ), and by the boundary condition in the upper region of the black bold curve (information comes from the line $r / r_{0} = 1$ ).

Each characteristic curve in the lower region shows the location of a particle that was at specific non-dimensional radius $r / r_{0}$ at $t = 0$ . Similarly, each characteristic curve in the upper region shows the location of a particle that was added to the body at time $τ$ ; the non-dimensional time $- {\overset{\cdot}{Z}}_{0} τ / r_{0}$ is given by the intersection of the characteristic curve with the line $r / r_{0} = 1$ . It should be noted that these path lines show the motion of the particles inside the growing body, so they are only valid within the growing body, although we have extended them outside the body. The time-dependent geometry of the growing body depends on the initial shape of the body and the boundary velocity; the latter depends on the rate of accretion or ablation and the particle velocity at the boundary, or it might be given, as in the case of the inner fixed boundary.

When there is no ablation ( $M_{1} = 0$ ) and the growth starts from a substrate with radius $r_{0}$ , i.e., there is no body at $t = 0$ , the black solid bold curve in Figure 3 shows the outer radius of the body at each time because only the motion of the particles at the outer boundary determines the velocity and location of the boundary. Then, the time-dependent growing body is a hollow sphere with a constant non-dimensional inner radius of one and a non-dimensional outer radius given by the black solid bold curve. The distinguishing constant $C_{1}$ in equation (31) for the outer radius in the special case of no ablation and no initial body is obtained by substituting $\hat{r} = r_{0}$ , $\hat{t} = l = 0$ and solving for $C_{1}$ , which gives

r^{3} = r_{0}^{3} - 3 r_{0}^{2} Z_{0} (t) .

(32)

This is the equation of the black solid bold curve.

However, if there exists an initial body with inner and outer radius $r_{0}$ and $r_{1} (t = 0)$ , and there is no ablation, the outer radius of the hollow sphere is defined by the characteristic line that was at $r = r_{1} (t = 0)$ at $t = 0$ (dashed bold black curve). For the general case with ablation and an initial body, the outer radius of the growing body (denoted by $r_{1} (t)$ in this problem) increases, but at a smaller rate than the case with no ablation, and the rate depends directly on the given function ${\overset{\cdot}{Z}}_{1}$ . In this case, the outer radius is depicted (schematically) by the gray dot-dashed curve in Figure 3. If this radius reaches a steady-state constant value, the body will be “treadmilling” [17].

Finally, we consider the case when there is ablation but no initial body, following [17]. The non-dimensional outer radius of $r_{1} (t) / r_{0}$ , which depends on ${\overset{\cdot}{Z}}_{1}$ is shown (schematically) by the solid gray curve and the growing body is shown by the gray area in Figure 3. As the inner radius is fixed, the lower limit of the growing body is the vertical line at non-dimensional radius $r / r_{0} = 1$ , and the upper limit of the area is the non-dimensional outer radius $r_{1} (t) / r_{0}$ at each time (the gray curve). As can be seen from the figure, $F_{e}$ in this case is only affected by the boundary conditions, since there can be no initial conditions in the absence of an initial body.

3.3.2. Evolution of the deformation gradient

To compute the non-zero components of deformation gradient, we consider the scalar equations (equations (28) and (29)) separately, as their right-hand sides are different.

Solving equation (28) for $F_{e_{rr}}$ . The differential equation of $F_{e_{rr}}$ along characteristic curves, which is parametrized with $l$ is

\frac{d F_{e_{rr}}}{d l} = 2 {\overset{\cdot}{Z}}_{0} \frac{r_{0}^{2}}{r^{3}} F_{e_{rr}} .

Combining this with equation (31) gives

\frac{d F_{e_{rr}}}{F_{e_{rr}}} = - 2 \frac{d r}{r} \Rightarrow F_{e_{rr}} = \frac{f (C_{1})}{r^{2}} = \frac{f (r^{3} + 3 Z_{0} r_{0}^{2})}{r^{2}} .

The unknown function $f$ is the constant of integration, which is different for each characteristic curve, and is a function of $C_{1}$ , which distinguishes different characteristic curves. It is to be determined using initial and boundary conditions. Consider the initial condition of $F_{e_{rr}} (r, t = 0) = F_{e_{r r_{0}}} (r)$ , so $f$ is as follows for the region influenced by the initial condition:

F_{e_{r r_{0}}} (r) = \frac{f (r^{3})}{r^{2}} \Rightarrow f (x) = x^{2 / 3} F_{e_{r r_{0}}} (x^{1 / 3}) .

The boundary condition is also $F_{e_{rr}} (r = r_{0}, t) = 1$ , so $f$ is as follows for the region influenced by the boundary condition:

1 = \frac{f (r_{0}^{3} + 3 Z_{0} (t) r_{0}^{2})}{r_{0}^{2}} \Rightarrow f (x) = const . = r_{0}^{2} .

The function $f$ is constant because the argument depends on $t$ but it is equal to constant value $r_{0}^{2}$ .

Putting these together, we have

F_{e_{rr}} (r, t) = {\begin{matrix} \frac{{(r^{3} + 3 Z_{0} r_{0}^{2})}^{\frac{2}{3}} F_{e_{r r_{0}}} ((r^{3} + 3 Z_{0} r_{0}^{2})^{\frac{1}{3}})}{r^{2}} & if r^{3} < r_{0}^{3} - 3 r_{0}^{2} Z_{0} (t), \\ {(\frac{r_{0}}{r})}^{2} & if r^{3} > r_{0}^{3} - 3 r_{0}^{2} Z_{0} (t) . \end{matrix}

Therefore, in the problem defined in [17], where there is no body at initial state but there is ablation at the outer surface, $F_{e_{rr}}$ is ${(r_{0} / r)}^{2}$ at all times. □

Solving equation (29) for $F_{e_{θ θ}}$ . The differential equation for $F_{e_{θ θ}}$ along the characteristic curves that are parametrized by $l$ is

\frac{d F_{e_{θ θ}}}{d l} = - {\overset{\cdot}{Z}}_{0} \frac{r_{0}^{2}}{r^{3}} F_{e_{θ θ}} .

Combining this with equation (31) gives

\frac{d F_{e_{θ θ}}}{F_{e_{θ θ}}} = \frac{d r}{r} \Rightarrow F_{e_{θ θ}} = rg (C_{1}) = rg (r^{3} + 3 Z_{0} r_{0}^{2}) .

The unknown function $g$ is the constant of integration, which is different for each characteristic curve, so it is a function of $C_{1}$ , which distinguishes different characteristic curves. Using the initial condition of $F_{e_{θ θ}} (r, t = 0) = F_{e_{θ θ_{0}}} (r)$ and the boundary condition, similar to the approach used in computing $F_{e_{rr}}$ , we have

F_{e_{θ θ}} (r, t) = {\begin{matrix} \frac{r F_{e_{θ θ_{0}}} ((r^{3} + 3 Z_{0} r_{0}^{2})^{\frac{1}{3}})}{{(r^{3} + 3 Z_{0} r_{0}^{2})}^{\frac{1}{3}}} & if r^{3} < r_{0}^{3} - 3 r_{0}^{2} Z_{0} (t), \\ \frac{r}{r_{0}} & if r^{3} > r_{0}^{3} - 3 r_{0}^{2} Z_{0} (t) . \end{matrix}

For the problem defined in [17], $F_{e_{θ θ}}$ is $r / r_{0}$ at all times. □

3.4. Solution for the special case with no initial body and with ablation at the outer surface

For the example defined in [17], where there is no initial body and $r^{3} > r_{0}^{3} - 3 r_{0}^{2} Z_{0} (t)$ , our calculation of $F_{e}$ is

F_{e} = {(\frac{r_{0}}{r})}^{2} {\hat{e}}_{r} \otimes {\hat{e}}_{r} + \frac{r}{r_{0}} ({\hat{e}}_{θ} \otimes {\hat{e}}_{θ} + {\hat{e}}_{ϕ} \otimes {\hat{e}}_{ϕ}) .

(33)

As the process of addition of material is smooth, there is no jump in traction boundary conditions to change $F_{e}$ instantaneously, and the added particles and the traction boundary condition are in equilibrium. Therefore, $F_{relax} = I$ .

Our calculation of $F$ matches with that computed in [17]; a superficial difference is that they represent the deformation gradient in the mixed basis with the first leg in the current and the second leg in the reference. Further, our calculation of the equations of motion (characteristic curves) derived in equation (31) matches with the current location of particles that was added to the body at time $t_{0} (Z)$ derived in [17] with $C_{1} = r_{0}^{3} + 3 r_{0}^{2} Z$ , where $t_{0} (Z)$ is the attachment time of a particle that was at location $Z$ in the reference configuration at the time of attachment. Therefore, essentially, both $Z$ and $\hat{r} (l = 0)$ are the same and the unit vectors in $r$ and $Z$ direction are identical. Further, the problem is independent of $θ$ and $ϕ$ and the coordinates are orthonormal; hence, the coordinate vectors defined here and [17] are the same, implying that the representations of $F$ are also identical.

Finally, substituting equation (33) into equation (25), the Cauchy stress tensor is

σ = (G {(\frac{r_{0}}{r})}^{4} - p) {\hat{e}}_{r} \otimes {\hat{e}}_{r} + (G {(\frac{r}{r_{0}})}^{2} - p) ({\hat{e}}_{θ} \otimes {\hat{e}}_{θ} + {\hat{e}}_{ϕ} \otimes {\hat{e}}_{ϕ}) .

To determine the unknown Lagrange multiplier function $p$ , the balance of linear momentum (equation (21)) and the boundary condition (equation (22)) are used. Thus, the stress field in the growing body is

\begin{matrix} \frac{σ}{G} = [\frac{1}{2} ({(\frac{r_{0}}{r})}^{4} - {(\frac{r_{0}}{r_{1} (t)})}^{4}) + {(\frac{r}{r_{0}})}^{2} - {(\frac{r_{1} (t)}{r_{0}})}^{2}] {\hat{e}}_{r} \otimes {\hat{e}}_{r} + \\ [- \frac{1}{2} ({(\frac{r_{0}}{r})}^{4} + {(\frac{r_{0}}{r_{1} (t)})}^{4}) + 2 {(\frac{r}{r_{0}})}^{2} - {(\frac{r_{1} (t)}{r_{0}})}^{2}] ({\hat{e}}_{θ} \otimes {\hat{e}}_{θ} + {\hat{e}}_{ϕ} \otimes {\hat{e}}_{ϕ}) \end{matrix}

(34)

which agrees with the Cauchy stress tensor derived in [17].

4. Hollow cylinder with accretion at the outer surface

In the previous example, accretion occurred at a surface with a known location. Here we will consider the case that accretion occurs at an unknown time-dependent surface. Specifically, we consider a long hollow cylinder made of a general incompressible material with time-dependent inner and outer radii, with accretion occurring on the outer surface and traction applied on the inner surface. The formulation is motivated by the wound roll process of additive manufacturing, which has been applied to processes ranging from the packaging of flexible films in paper and magnetic tape industries [33] to tissue engineered blood vessels [34]. The material is an incompressible and isotropic hyperelastic solid with an energy function $W = W (I_{1}, I_{2})$ .

We assume that we can neglect the end effects and assume plain strain. Therefore, we use a two-dimensional polar coordinate system:

{\hat{e}}_{r} = \frac{\partial x}{\partial r}, {\hat{e}}_{θ} = \frac{1}{r} \frac{\partial x}{\partial θ} .

(35)

This problem was solved in [14] using a geometrical Lagrangian approach with two different choices of metric for the initial and added parts of the body. We discuss next the fact that the choice of the growth velocity does not change the Eulerian formulation, but it can make the Lagrangian formulation more complicated as the time of attachment has to be related to the location in the reference configuration. We also show that one can solve an Eulerian growth process using the evolution of inverse of the motion via the method discussed in [29], rather than computing the deformation gradient directly using equation (9).

4.1. Boundary and initial conditions

We denote the time-dependent inner and outer radii by $r_{in} (t)$ and $r_{out} (t)$ , their initial values by $r_{in} (t = 0) = R_{0}$ and $r_{out} (t = 0) = R_{1}$ , and the radial velocities by ${\overset{\cdot}{r}}_{out}$ and ${\overset{\cdot}{r}}_{in}$ , respectively.

The body is initially stress-free, and an internal pressure of $p_{i} (t)$ ( $p_{i} (t = 0) = 0$ ) applies a traction at the inner time-dependent surface of the cylinder, shown in Figure 4.

Figure 4.

Very long hollow cylinder body with outer accretion: the gray body is the initial body and the dashed lines show the time-dependent boundary of the growing body

There is no ablation, and the accretion occurs only at the outer surface. The mass rate of growth $M$ at the outer surface per unit area is equal to $ρ u_{g} (t)$ . The added particles are stress-free and the growing surface is traction-free. We aim to find the time-dependent geometry of the body, i.e., the expression for $r_{in} (t)$ and $r_{out} (t)$ , and the stress distribution in the body.

4.2. Reduced governing equations

From the symmetry of the problem, we assume that no quantities of interest depend on $θ$ , and that the velocity $υ$ is radial and can be written $υ = v_{r} (r, t) {\hat{e}}_{r}$ . Therefore, the governing equations, written in the coordinate system defined in equation (35), reduce to the following:

div (υ) = \frac{\partial v_{r}}{\partial r} + \frac{v_{r}}{r} = 0 on r_{in} (t) < r < r_{out} (t)

(36)

{\overset{\cdot}{r}}_{out} = V_{b} \cdot {\hat{e}}_{r} = v_{r} (r = r_{out} (t), t) + u_{g} (t) at r = r_{out} (t)

(37)

{\overset{\cdot}{r}}_{in} = V_{b} \cdot {\hat{e}}_{r} = v_{r} (r = r_{in} (t), t) at r = r_{in} (t)

(38)

div (σ) = (\frac{\partial σ_{rr}}{\partial r} + \frac{1}{r} (σ_{rr} - σ_{θ θ})) {\hat{e}}_{r} = 0 on r_{in} (t) < r < r_{out} (t)

(39)

t_{b} = 0 at r = r_{out} (t)

(40)

t_{b} = p_{i} (t) {\hat{e}}_{r} at r = r_{in} (t)

(41)

\frac{\partial}{\partial t} [\begin{matrix} F_{e_{rr}} & F_{e_{r θ}} \\ F_{e_{θ r}} & F_{e_{θ θ}} \end{matrix}] + v_{r} \frac{\partial}{\partial r} [\begin{matrix} F_{e_{rr}} & F_{e_{r θ}} \\ F_{e_{θ r}} & F_{e_{θ θ}} \end{matrix}] = [\begin{matrix} \frac{\partial v_{r}}{\partial r} & 0 \\ 0 & \frac{v_{r}}{r} \end{matrix}] [\begin{matrix} F_{e_{rr}} & F_{e_{r θ}} \\ F_{e_{θ r}} & F_{e_{θ θ}} \end{matrix}]

(42)

F_{e} = I at t = 0, R_{0} < r < R_{1}, and F_{e} = I at r = r_{out} (t)

(43)

\begin{matrix} \hat{σ} (F_{e}) = - p I + \frac{\partial W (F_{e})}{\partial F_{e}} F_{e}^{T} = (- p + 2 I_{2} \frac{\partial W}{\partial I_{2}}) I + 2 \frac{\partial W}{\partial I_{1}} F_{e} F_{e}^{T} - 2 \frac{\partial W}{\partial I_{2}} {(F_{e} F_{e}^{T})}^{- 1}, \\ where I_{1} = tr (B_{e}) and I_{2} = \frac{tr {(B_{e})}^{2} - tr (B_{e}^{2})}{2} are invariants of B_{e} = F_{e} F_{e}^{T} . \end{matrix}

(44)

As the initial and inlet boundary conditions of $F_{e}$ are consistent, there is no jump in traction boundary conditions to change $F_{e}$ instantly. Also, the added particles and the traction boundary conditions have consistent stresses. So, we also keep the boundary and initial conditions for $F$ equal to $I$ , and the relaxed shape of the added particles is equal to their shape at the time of attachment; then $F_{relax} = I$ .

4.3. Solution

The solution of the continuity equation (equation (36)) using the boundary condition of equation (38) is

υ = \frac{r_{in}}{r} {\overset{\cdot}{r}}_{in} {\hat{e}}_{r}

(45)

and $r_{in} (t)$ , and consequently ${\overset{\cdot}{r}}_{in} (t)$ , will be computed later. Substituting this in equation (37), we obtain a differential equation that can be later solved for the outer radius:

r_{in} {\overset{\cdot}{r}}_{in} = r_{out} {\overset{\cdot}{r}}_{out} - r_{out} u_{g} (t), with r_{out} (t = 0) = R_{1} and r_{in} (t = 0) = R_{0} .

(46)

Similar to the previous example, $\nabla υ$ is diagonal, so the four PDEs for the evolution of the components of $F_{e}$ are decoupled and homogeneous. As both the boundary and initial conditions for $F_{e}$ are identity, the off-diagonal components of $F_{e}$ remain zero during the growth process. Using (42) and (45), the evolution equations for the diagonal components of $F_{e}$ are:

\frac{\partial F_{e_{rr}}}{\partial t} + {\overset{\cdot}{r}}_{in} \frac{r_{in}}{r} \frac{\partial F_{e_{rr}}}{\partial r} = - {\overset{\cdot}{r}}_{in} \frac{r_{in}}{r^{2}} F_{e_{rr}},

(47)

\frac{\partial F_{e_{θ θ}}}{\partial t} + {\overset{\cdot}{r}}_{in} \frac{r_{in}}{r} \frac{\partial F_{e_{θ θ}}}{\partial r} = {\overset{\cdot}{r}}_{in} \frac{r_{in}}{r^{2}} F_{e_{θ θ}} .

(48)

We solve these using the method of characteristics, similar to the previous example. The parametric equations of the characteristic curves are

\frac{d t}{d l} = 1 \Rightarrow \hat{t} = l + C_{0},

(49)

\frac{d r}{d l} = {\overset{\cdot}{r}}_{in} \frac{r_{in}}{r} .

(50)

As in the previous example, we assume that all the characteristic curves start at $t = 0$ , so $\hat{t} (l = 0) = 0 \Rightarrow C_{0} = 0$ for all characteristic curves. Further, substituting $d l = d t$ in equation (50) and using the identity

{\overset{\cdot}{r}}_{in} r_{in} = \frac{1}{2} \frac{d}{d t} r_{in}^{2},

the function $\hat{r} (l)$ along the characteristic curves is

r d r = \frac{1}{2} ({\overset{\cdot}{r}}_{in}^{2}) d t \Rightarrow {\hat{r}}^{2} (l) + C_{1} = r_{in}^{2} (\hat{t} (l)) .

(51)

The integration constant $C_{1}$ distinguishes different characteristic curves. Obviously, the path line with $C_{1} = 0$ shows the location of the inner boundary. However, the equation of $r_{in} (t)$ has not been computed yet, and will be computed later, using the traction boundary conditions. Further, given $r_{in} (t)$ , we can compute $r_{out} (t)$ using equation (46). Figure 5 shows the path lines and body shape.

Figure 5.

Path lines of the motion in space–time. The outer radius is marked by the bold dot-dashed line, and the inner radius by the solid bold line; the time-dependent growing body lies between these radii. The characteristic curve starting at $r = R_{0}$ and $t = 0$ is shown by the solid lines. The colors mark the space–time regions corresponding to the initial body and the accreted mass.

We can also substitute for ${\overset{\cdot}{r}}_{in} r_{in}$ in equation (51) from equation (46) to find the equation for the characteristic curves:

{\hat{r}}^{2} (l) + C_{2} = r_{out}^{2} (\hat{t} (l)) - 2 \int_{0}^{\hat{t} (l)} u_{g} (t) r_{out} (t) d t,

(52)

where $C_{2}$ is a constant of integration, and is different for different characteristic curves.

Using equations (47) and (48), $F_{e_{rr}}$ and $F_{e_{θ θ}}$ along the characteristic curves are governed by

\frac{d F_{e_{rr}}}{d l} = - {\overset{\cdot}{r}}_{in} \frac{r_{in}}{r^{2}} F_{e_{rr}}, \frac{d F_{e_{θ θ}}}{d l} = {\overset{\cdot}{r}}_{in} \frac{r_{in}}{r^{2}} F_{e_{θ θ}} .

Combining this with equation (51), the general forms of $F_{rr}$ and $F_{θ θ}$ are

\hat{r} (l, C) F_{e_{rr}} = g_{1} (C),

(53)

\frac{F_{e_{θ θ}}}{\hat{r} (l, C)} = g_{2} (C),

(54)

where $g_{1}$ and $g_{2}$ are constants of integration, which are different for each characteristic line. Therefore, $g_{1}$ and $g_{2}$ are functions of $C$ which is either $C_{1}$ or $C_{2}$ , depending on whether we use equation (51) or equation (52) to distinguish the different characteristic curves.

Considering the characteristic curve starting at $r = R_{1}$ at $t = 0$ (marked by a dashed bold line in Figure 5), we substitute $(r = R_{1}, l = 0)$ to find $C_{1}$ in equation (51); hence, the equation of this characteristic curve is

R_{0}^{2} - R_{1}^{2} = r_{in}^{2} (\hat{t} (l)) - {\hat{r}}_{R_{1}}^{2} (l),

(55)

where ${\hat{r}}_{R_{1}}^{2} (l)$ is the function $\hat{r}$ for this specific characteristic curve. This curve divides the body into two regions. As can be seen from Figure 5, the characteristic curves to the left of this specific curve (green region) originate from the line $t = 0$ ; therefore, this part of the body is governed by the initial conditions of the body. Conversely, the characteristic curves to the right of this specific curve (red region) originate from the outer boundary; therefore, this part of the body is governed by the boundary condition on the outer surface. Consequently, $F_{e}$ and $σ$ behave differently in each of these regions, and they must be treated separately.

4.3.1. The region that was initially in the body

Considering the region that formed part of the initial body, the characteristic curves from equation (51) can be used to substitute $C$ in equations (53) and (54) to find the general forms of $F_{e_{rr}}$ and $F_{e_{θ θ}}$ :

\hat{r} F_{e_{rr}} (\hat{r}, \hat{t}) = g_{1} (C_{1}) = g_{1} (r_{in}^{2} (\hat{t}) - {\hat{r}}^{2}),

(56)

\frac{F_{e_{θ θ}} (\hat{r}, \hat{t})}{\hat{r}} = g_{2} (C_{1}) = g_{2} (r_{in}^{2} (\hat{t}) - {\hat{r}}^{2}) .

(57)

Using the initial condition of $F_{e_{rr}} = F_{e_{θ θ}} = 1$ at $t = 0$ to compute $g_{1}$ and $g_{2}$ , and substituting $\hat{r}$ and $\hat{t}$ with $r$ and $t$ , we can find $F_{e}$ in this region:

F_{e} = \frac{\sqrt{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}}{r} {\hat{e}}_{r} \otimes {\hat{e}}_{r} + \frac{r}{\sqrt{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}} {\hat{e}}_{θ} \otimes {\hat{e}}_{θ} .

(58)

Using equation (51), for a particle that was initially ( $l = 0$ ) at $R$ , the spatial location at any later time $t$ is

r^{2} = R^{2} + r_{in}^{2} (t) - R_{0}^{2} \Rightarrow \sqrt{R_{0}^{2} - r_{in}^{2} (t) + r^{2}} = R .

(59)

So, we can write $F_{e}$ as

F_{e} = \frac{R}{r} {\hat{e}}_{r} \otimes {\hat{e}}_{r} + \frac{r}{R} {\hat{e}}_{θ} \otimes {\hat{e}}_{θ} .

We notice that this agrees with [14], because the initial body is stress-free, and hence the initial configuration is also the reference configuration for this part of the body.

Substituting $F_{e}$ from equation (58) into equation (44), the Cauchy stress tensor in this region is

\begin{matrix} σ = [- p + 2 \frac{\partial W}{\partial I_{2}} (1 - \frac{r^{2}}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}) + 2 \frac{\partial W}{\partial I_{1}} (\frac{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}{r^{2}})] {\hat{e}}_{r} \otimes {\hat{e}}_{r} \\ + [- p + 2 \frac{\partial W}{\partial I_{2}} (1 - \frac{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}{r^{2}}) + 2 \frac{\partial W}{\partial I_{1}} (\frac{r^{2}}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}})] {\hat{e}}_{θ} \otimes {\hat{e}}_{θ} . \end{matrix}

(60)

To compute the unknown pressure field $p$ , we use equation (39) and the appropriate boundary condition (equation (41)):

\begin{matrix} \frac{\partial σ_{rr}}{\partial r} + \frac{1}{r} (σ_{rr} - σ_{θ θ}) = 0 \Rightarrow \\ \frac{\partial σ_{rr}}{\partial r} = \frac{2}{r} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (\frac{r^{4} - {(R_{0}^{2} - r_{in}^{2} (t) + r^{2})}^{2}}{(R_{0}^{2} - r_{in}^{2} (t) + r^{2}) r^{2}}) \\ = \frac{2 r}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{{(R_{0}^{2} - r_{in}^{2} (t) + r^{2})}^{2}}{r^{4}}) . \end{matrix}

(61)

As $σ_{rr}$ is a function of $r$ and $t$ , we have

σ_{rr} = \int_{r_{in} (t)}^{r} [\frac{2 r}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{{(R_{0}^{2} - r_{in}^{2} (t) + r^{2})}^{2}}{r^{4}})] d r + h (t) .

Using the boundary condition (equation (41)), $h (t)$ is equal to $- p_{i} (t)$ , thus

σ_{rr} = - p_{i} (t) + \int_{r_{in} (t)}^{r} [\frac{2 r}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{{(R_{0}^{2} - r_{in}^{2} (t) + r^{2})}^{2}}{r^{4}})] d r

(62)

and

p (r, t) = p_{i} (t) - \int_{r_{in} (t)}^{r} [\frac{r}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{{(R_{0}^{2} - r_{in}^{2} (t) + r^{2})}^{2}}{r^{4}})] d r

(63)

+ 2 \frac{\partial W}{\partial I_{2}} (1 - \frac{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}{r^{2}}) + 2 \frac{\partial W}{\partial I_{1}} (\frac{r^{2}}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}) .

(64)

Using $p$ from equation (63) in equation (60) gives the complete Cauchy stress tensor in this region of the body.

To compare this result with that given in [14], we need to write equation (62) based on the reference configuration and in Lagrangian form. Using equation (59), we have $r d r = R d R$ ; and substituting this relation in equation (62), we can write $σ_{rr}$ in the following form that agrees with [14]:

σ_{rr} = - p_{i} (t) + \int_{R_{1}}^{R} [\frac{2}{R} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{R^{4}}{r^{4}})] d R .

4.3.2. The region that is added to the body through the growth process

In the region that is governed by the outer boundary conditions, it is more convenient to use equation (52) to compute the characteristic curves as it involves the outer part of the body, in contrast to equation (51). So, for $C$ we use $C_{2}$ in equations (53) and (54). Hence, the general forms of $F_{e_{rr}}$ and $F_{e_{θ θ}}$ are:

\hat{r} F_{e_{rr}} (\hat{r}, \hat{t}) = g_{1} (C_{2}) = g_{1} (r_{out}^{2} (\hat{t}) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t)

(65)

\frac{F_{e_{θ θ}} (\hat{r}, \hat{t})}{\hat{r}} = g_{2} (C_{2}) = g_{2} (r_{out}^{2} (\hat{t}) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t)

(66)

Using the boundary condition from equation (43) at $r = r_{out} (t)$ , we have

r_{out} (\hat{t}) = g_{1} (- 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t),

(67)

\frac{1}{r_{out} (\hat{t})} = g_{2} (- 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t) .

(68)

To compute $F_{e_{rr}}$ and $F_{e_{θ θ}}$ , we need to know the functions $g_{1}$ and $g_{2}$ when their arguments are

r_{out}^{2} (t) - {\hat{r}}^{2} - 2 \int_{0}^{t} u_{g} (t) r_{out} (t) d t .

Using equation (52) for a particle that is initially at the outer radius of $R_{1}$ , substituting $(r = r_{out} (t = 0) = R_{1}, t = 0)$ , the constant $C_{2}$ for this characteristic curve is computed to be zero. Hence, another form of ${\hat{r}}_{R_{1}} (l)$ from equation (55) is as follows:

{\hat{r}}_{R_{1}}^{2} (l) = r_{out}^{2} (\hat{t}) - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t .

(69)

In Figure 5, ${\hat{r}}_{R_{1}} (l)$ is indicated by the bold dashed line, and $r_{out} (t)$ by the bold dotted dashed line.

Further, for a particle that is at $r = r_{out} (τ)$ at $t = τ$ , substituting $(r = r_{out} (t = τ), t = τ)$ , the constant $C_{2}$ is $- 2 \int_{0}^{τ} u_{g} (t) r_{out} (t) d t$ . Thus, the path line of a particle that attaches to the body at the outer radius $r_{out} (τ)$ (at time $t = τ$ ) at any time after the attachment time $t > τ$ is

r_{out}^{2} (\hat{t}) - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t = {\hat{r}}^{2} - 2 \int_{0}^{τ} u_{g} (t) r_{out} (t) d t .

(70)

Using this equation and the boundary conditions (equations (67) and (68)), we get

r_{out}^{2} (\hat{t}) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t = - 2 \int_{0}^{τ} u_{g} (t) r_{out} (t) d t \Rightarrow

(71)

\hat{r} F_{e_{rr}} (\hat{r}, \hat{t}) = g_{1} (r_{out}^{2} (\hat{t}) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t) = g_{1} (- 2 \int_{0}^{τ} u_{g} (t) r_{out} (t) d t) = r_{out} (τ),

(72)

\frac{F_{e_{θ θ}} (\hat{r}, \hat{t})}{\hat{r}} = g_{2} (r_{out}^{2} (\hat{t}) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t) = g_{2} (- 2 \int_{0}^{τ} u_{g} (t) r_{out} (t) d t) = \frac{1}{r_{out} (τ)} .

(73)

So, the deformation gradient in this region is

F_{e} = \frac{r_{out} (τ (\hat{r}, \hat{t}))}{r} {\hat{e}}_{r} \otimes {\hat{e}}_{r} + \frac{r}{r_{out} (τ (\hat{r}, \hat{t}))} {\hat{e}}_{θ} \otimes {\hat{e}}_{θ} .

(74)

Now, we need to compute $τ$ for a particle at a given spatial location $r$ at given time $t$ . Substituting equation (70) into equation (69), the following relation implicitly provides an equation to compute $τ$ for a particle that is at the spatial radius $r$ at time $t$ with ${\hat{r}}_{R_{1}} (t) < r < r_{out} (t)$ :

{\hat{r}}^{2} (l) = {\hat{r}}_{R_{1}}^{2} + 2 \int_{0}^{τ} u_{g} (t) r_{out} (t) d t .

(75)

The (as yet unknown) functions $r_{out} (\hat{t})$ and ${\hat{r}}_{R_{1}}$ are required to apply this equation to compute $τ$ . The functions $r_{out} (\hat{t})$ and ${\hat{r}}_{R_{1}}$ are related to each other by equation (69), so indeed we have only one unknown function. For now, we will continue with these unknown functions, and they will be computed using the traction boundary conditions next.

Although we want to compute $τ$ for a known spatial location $r$ at known time $t$ because of the Eulerian approach, Sozio and Yavari [14] used a similar equation to compute the unknown location $r$ at known time $t$ of a particle that has reference position $R$ in the Lagrangian approach. In that approach, one needs to also know the relation between $τ$ and $R$ as the attachment time affects the current location of the body, and the integral in the last term of equation (75) is written with respect to $d R$ in [14]. For simplicity, Sozio and Yavari [14] assume that $u_{g} (t)$ is constant and equal to $u_{0}$ , and then $τ = (R - R_{1}) / u_{0}$ . In general, the relation between $τ$ to $R$ is more complicated. The Eulerian approach does not refer to the reference configuration explicitly; hence, the choice of $u_{g} (t)$ does not affect the formulation of the problem.

We notice that we could alternatively find $F_{e}$ from equations (65) and (66) using the following argument. The boundary condition for this region is that the particles are added to the body at the outer surface with radius $r_{out} (\hat{t})$ without experiencing any deformation, i.e. $F_{e} = I$ . In other words, $F_{e_{rr}} = F_{e_{θ θ}} = 1$ at the attachment radius of the particle $r = r_{attach} (\hat{r}, \hat{t})$ ; this is different for each characteristic curve. Also, we know that, along path lines, $\hat{r} F_{e_{rr}}$ and $F_{e_{θ θ}} / \hat{r}$ are constants and equal to $g_{1} (C_{2})$ and $g_{2} (C_{2})$ , respectively. Using the boundary condition, $r_{attach} (\hat{r}, \hat{t}) = g_{1} (C_{2}) = 1 / g_{2} (C_{2})$ , and we can then write $F_{e}$ in this region as

F_{e} = \frac{r_{attach} (\hat{r}, \hat{t})}{r} {\hat{e}}_{r} \otimes {\hat{e}}_{r} + \frac{r}{r_{attach} (\hat{r}, \hat{t})} {\hat{e}}_{θ} \otimes {\hat{e}}_{θ}

And obviously the attachment radius of each particle is the outer radius of the body at the time of attachment $τ$ , so $r_{attach} (\hat{r}, \hat{t}) = r_{out} (τ)$ ; equation (75) gives us an equation to compute $τ$ .

Now that $F_{e}$ is known (equation (74)), we can compute the Cauchy stress using equation (44):

\begin{matrix} σ = [- p + 2 \frac{\partial W}{\partial I_{2}} (1 - \frac{r^{2}}{r_{out}^{2} (τ)}) + 2 \frac{\partial W}{\partial I_{1}} (\frac{r_{out}^{2} (τ)}{r^{2}})] {\hat{e}}_{r} \otimes {\hat{e}}_{r} \\ + [- p + 2 \frac{\partial W}{\partial I_{2}} (1 - \frac{r_{out}^{2} (τ)}{r^{2}}) + 2 \frac{\partial W}{\partial I_{1}} (\frac{r^{2}}{r_{out}^{2} (τ)})] {\hat{e}}_{θ} \otimes {\hat{e}}_{θ} . \end{matrix}

(76)

To find the unknown pressure field $p$ , we use equation (39) to compute $σ_{rr}$ :

σ_{rr} = \int_{{\hat{r}}_{R_{1}}}^{r} [\frac{1}{r} (σ_{θ θ} - σ_{rr})] d r = \int_{{\hat{r}}_{R_{1}}}^{r} [\frac{2 r}{r_{out}^{2} (τ)} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{r_{out}^{4} (τ)}{r^{4}})] d r + h (t) .

Using traction continuity at the interface between the initial and added material (located at ${\hat{r}}_{R_{1}}$ ), we have for $σ_{rr}$ and the function $h (t)$ :

\begin{matrix} h (t) = - p_{i} (t) + \int_{r_{in} (t)}^{r (R_{1}, t)} [\frac{2 r}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{{(R_{0}^{2} - r_{in}^{2} (t) + r^{2})}^{2}}{r^{4}})] d r, \\ σ_{rr} = - p_{i} (t) + \int_{r_{in} (t)}^{r (R_{1}, t)} [\frac{2 r}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{{(R_{0}^{2} - r_{in}^{2} (t) + r^{2})}^{2}}{r^{4}})] d r \\ + \int_{r (R_{1}, t)}^{r} [\frac{2 r}{r_{out}^{2} (τ)} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{r_{out}^{4} (τ)}{r^{4}})] d r . \end{matrix}

We note that $r_{in} (t)$ , $r_{out} (t)$ , and ${\hat{r}}_{R_{1}}$ are as yet unknown functions. Based on equations (46) and (55), these functions are not independent, and we can write both ${\hat{r}}_{R_{1}}$ and $r_{out} (t)$ in terms of $r_{in} (t)$ . Another boundary condition is that the outer surface $r = r_{out} (t)$ is stress-free. We therefore get the following equation for computing these three (non-independent) unknown functions:

\begin{matrix} p_{i} (t) = \int_{r_{in} (t)}^{r (R_{1}, t)} [\frac{2 r}{R_{0}^{2} - r_{in}^{2} (t) + r^{2}} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{{(R_{0}^{2} - r_{in}^{2} (t) + r^{2})}^{2}}{r^{4}})] d r \\ + \int_{r (R_{1}, t)}^{r_{out} (t)} [\frac{2 r}{r_{out}^{2} (τ)} (\frac{\partial W}{\partial I_{1}} + \frac{\partial W}{\partial I_{2}}) (1 - \frac{r_{out}^{4} (τ)}{r^{4}})] d r . \end{matrix}

(77)

This is in agreement with [14].

4.4. Computing the deformation gradient using the inverse of the motion

An alternative Eulerian approach to find $F$ is to work with the inverse of the motion $ξ = ξ_{r} {\hat{e}}_{r} + ξ_{θ} {\hat{e}}_{θ}$ and equation (78) provided in [29]; once we have $ξ$ , we can compute $F$ using $F = (\nabla ξ)^{- 1}$ . And we note that as $F_{relax} = I$ , $F_{e} = F$ ,

\frac{\partial ξ^{- 1}}{\partial t} + (υ \cdot \nabla) ξ^{- 1} = 0 .

(78)

From symmetry, we have $ξ_{θ} = const .$ For $ξ_{r}$ , we have

\frac{\partial ξ_{r}}{\partial t} + {\overset{\cdot}{r}}_{in} \frac{r_{in}}{r} \frac{\partial ξ_{r}}{\partial r} = 0,

(79)

which is the $r$ -component of equation (78), after substituting the velocity field from equation (45).

Using the method of characteristics, the characteristic curves are the same as those from equations (51) and (52). However, as the right-hand side of this equation is zero, it means that $ξ_{r}$ is constant along the characteristics, giving

ξ_{r} = g (C) = g_{1} (r_{in}^{2} (\hat{t}) - {\hat{r}}^{2}) = g_{2} (r_{out}^{2} (\hat{t}) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t) .

(80)

We can then compute $F$ in the initial region:

F = (\nabla ξ)^{- 1} = {[\begin{matrix} \frac{\partial ξ_{r}}{\partial r} & 0 \\ 0 & \frac{ξ_{r}}{r} \end{matrix}]}^{- 1} = {[\begin{matrix} \frac{1}{- 2 r g_{1'} (m)} & 0 \\ 0 & \frac{r}{g_{1} (m)} \end{matrix}]}_{m = r_{in}^{2} (t) - r^{2}} .

(81)

This region is influenced by the initial condition, which is $F = I$ . Also, the intersection between the line $t = 0$ and a characteristic curve that passes through $\hat{r}$ at time $\hat{t}$ is $R$ (the initial location of the particle), so

g_{1} (R_{0}^{2} - R^{2}) = R \Rightarrow g_{1} (m) = \sqrt{R_{0}^{2} - m} \Rightarrow g_{1'} (m) = \frac{- 1}{2 \sqrt{R_{0}^{2} - m}} .

(82)

Therefore, the explicit expression for $F$ is

F = [\begin{matrix} \frac{\sqrt{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}}{r} & 0 \\ 0 & \frac{r}{\sqrt{R_{0}^{2} - r_{in}^{2} (t) + r^{2}}} \end{matrix}],

(83)

which is, as expected, identical to equation (58).

In the region governed by the boundary condition, we have from equation (80) that $ξ_{r} = g_{2} (C_{2})$ . Combined with the boundary condition $F = I$ at $r = r_{out} (t)$ , we have

F = {[\begin{matrix} \frac{1}{- 2 r g_{2'} (m)} & 0 \\ 0 & \frac{r}{g_{2} (m)} \end{matrix}]}_{m = (r_{out}^{2} (\hat{t} (l)) - {\hat{r}}^{2} (l, f) - 2 \int_{0}^{\hat{t} (l)} u_{g} (t) r_{out} (t) d t)} .

(84)

At $r = r_{out} (\hat{t}),$ we have $F_{θ θ} = 1$ , so $r_{out} (\hat{t}) = g_{2} (- 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t)$ . Using equation (70):

g_{2} (m) |_{m = (r_{out}^{2} (\hat{t}) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t)} = g_{2} (- 2 \int_{0}^{τ} u_{g} (t) r_{out} (t) d t) = r_{out} (τ (\hat{r}, \hat{t})) .

Now, we need $g_{2'} (m)$ to define $F$ completely; however, the general form of $g_{2} (m)$ is not readily available to differentiate. We only know that $g_{2} (m) = r_{out} (τ)$ at $m = r_{out}^{2} (\hat{t} (l)) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t} (l)} u_{g} (t) r_{out} (t) d t$ . Therefore, we consider $g' (m)$ as an unknown function, and apply the boundary condition $F_{rr} = 1$ at $r = r_{out}$ to compute $g_{2'} (m)$ directly:

\begin{matrix} - 2 r {g'}_{2} (- 2 \int_{0}^{\hat{t} (l)} u_{g} (t) r_{out} (t) d t) = 1 \Rightarrow {g'}_{2} (- 2 \int_{0}^{\hat{t}} u_{g} (t) r_{out} (t) d t) = \frac{1}{- 2 r_{out} (\hat{t})} \Rightarrow \\ {g'}_{2} (r_{out}^{2} (\hat{t} (l)) - {\hat{r}}^{2} - 2 \int_{0}^{\hat{t} (l)} u_{g} (t) r_{out} (t) d t) = {g'}_{2} (- 2 \int_{0}^{τ} u_{g} (t) r_{out} (t) d t) = \frac{1}{- 2 r_{out} (τ)} \Rightarrow \\ F_{rr} = \frac{r_{out} (τ)}{r} . \end{matrix}

Therefore, $F$ in the region governed by the boundary conditioned region is

F = [\begin{matrix} \frac{r_{out} (τ)}{r} & 0 \\ 0 & \frac{r}{r_{out} (τ)} \end{matrix}],

(85)

which is identical to equation (74).

Instead of using the initial condition on $F$ , we could use that the region governed by the initial condition is initially stress-free to set $ξ_{r} = R$ , and use this to compute the function $g_{1} (C)$ in the initially existing body:

g_{1} (R_{0}^{2} - R^{2}) = R \Rightarrow ξ_{r} = g_{1} (C_{2}) = \sqrt{R_{0}^{2} - r_{in}^{2} (\hat{t}) - {\hat{r}}^{2}} .

However, in the region governed by the boundary condition, converting the boundary condition on $F$ to a boundary condition on $ξ_{r}$ requires the initial location of the added particles, and hence has the same difficulty as computing the reference configuration. This is the reason that we compute $F$ using equation (8) rather than through computing the inverse of the motion using equation (78).

5. Discussion

In this paper, we model the surface growth of solid bodies using an Eulerian approach. An important advantage of this approach is that we do not need to compute the time-evolving reference configuration explicitly. Conversely, to obtain the stress, we need to solve for the deformation gradient, which is not straightforward in the Eulerian setting. This is addressed by introducing the elastic deformation field and the associated transport equation that governs its evolution. To solve this time-dependent transport equation, we use the method of characteristics that enables us to find closed-form solutions in two examples with spherical and cylindrical symmetry with free boundaries.

The primary motivation for formulating the Eulerian approach is to enable future work on numerical solutions. A potential advantage of Eulerian methods that have been developed in the fluid–structure interaction (FSI) literature (e.g. [36]) is that we can solve such problems on a fixed mesh for an evolving body, rather than an evolving reference configuration and deforming mesh that needs to be updated or re-meshed frequently. A numerical scheme would enable the careful study of various interesting systems identified recently, e.g. [22, 23, 37], which focus on the problem of near-net-shape manufacturing in the presence of incompatibility and residual stress, as well as [38], which discusses incompatibility in both surface and bulk growth.

An important challenge in the Eulerian formulation is to know the location of the boundaries to be able to apply the boundary conditions. We refer to the brief discussion in Section 3.3.1 of [24], as well as [29, 30], and, noticing that the velocity field can be used to locate the boundary of the body for the parts of the boundary without growth, we follow the characteristic lines of particles located on the boundary; for the parts with growth, we use equation (3) together with the path lines to determine the location of the growing boundaries.

Briefly, we also highlight that bulk growth has received much attention, e.g., using growth tensors [39 –42], geometrical approaches [43], and mixture theory [44]. A key difference between bulk and surface growth is that the former can typically be studied by assuming that the set of the material points is unchanged, and that growth occurs through a change in the referential density, owing to a distributed mass source.

Footnotes

Acknowledgements

We thank Rohan Abeyaratne, Timothy Breitzman, Tal Cohen, and Tony Rollett for useful discussions.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Science Foundation (grant numbers CMMI MOMS 1635407, DMREF 2118945, DMS 1729478, DMS 2012259, DMS 2108784), the Army Research Office (grant number MURI W911NF-19-1-0245), the Office of Naval Research (grant number N00014-18-1-2528), and the Binational Science Foundation (grant number 2018183).

ORCID iDs

S Kiana Naghibzadeh

Kaushik Dayal

Notes

References

Skalak

Dasgupta

Moss

, et al. Analytical description of growth. J Theor Biol 1982; 94(3): 555–577.

Taber

. Biomechanics of growth, remodeling, and morphogenesis. Appl Mech Rev 1995; 48(8): 487–545.

Brown

Goodman

. Gravitational stresses in accreted bodies. Proc R Soc London, Ser A 1963; 276(1367): 571–576.

Kadish

Barber

Washabaugh

, et al. Stresses in accreted planetary bodies. Int J Solids Struct 2008; 45(2): 540–550.

Schwerdtfeger

Sato

Tacke

. Stress formation in solidifying bodies: Solidification in a round continuous casting mold. Metall Mater Trans B 1998; 29(5): 1057–1068.

Rao

Hughes

. On modelling thermal oxidation of silicon I: Theory. Int J Numer Methods Eng 2000; 47(1–3): 341–358.

Drozdov

. Continuous accretion of a composite cylinder. Acta Mech 1998; 128(1–2): 117–135.

Bacigalupo

Gambarotta

. Effects of layered accretion on the mechanics of masonry structures. Mech Based Des Struct Mach 2012; 40(2): 163–184.

Papadopoulos

Hodge

. On surface growth of actin networks. Int J Eng Sci 2010; 48(11): 1498–1506.

10.

Ong

O’Reilly

. On the equations of motion for rigid bodies with surface growth. Int J Eng Sci 2004; 42(19–20): 2159–2174.

11.

Naumov

. Mechanics of growing deformable solids: A review. J Eng Mech 1994; 120(2): 207–220.

12.

Kadish

Barber

Washabaugh

. Stresses in rotating spheres grown by accretion. Int J Solids Struct 2005; 42(20): 5322–5334.

13.

Skalak

Farrow

Hoger

. Kinematics of surface growth. J Math Biol 1997; 35(8): 869–907.

14.

Sozio

Yavari

. Nonlinear mechanics of surface growth for cylindrical and spherical elastic bodies. J Mech Phys Solids 2017; 98: 12–48.

15.

Sozio

Shojaei

Sadik

, et al. Nonlinear mechanics of thermoelastic accretion. Z Angew Math Phys 2020; 71: 87.

16.

Sozio

Yavari

. Nonlinear mechanics of accretion. J Nonlinear Sci 2019; 29(4): 1813–1863.

17.

Tomassetti

Cohen

Abeyaratne

. Steady accretion of an elastic body on a hard spherical surface and the notion of a four-dimensional reference space. J Mech Phys Solids 2016; 96: 333–352.

18.

Abi-Akl

Abeyaratne

Cohen

. Kinetics of surface growth with coupled diffusion and the emergence of a universal growth path. Proc R Soc London, Ser A 2019; 475(2221): 20180465.

19.

Abi-Akl

Cohen

. Surface growth on a deformable spherical substrate. Mech Res Commun 2020; 103: 103457.

20.

Abeyaratne

Puntel

Tomassetti

. Treadmilling stability of a one-dimensional actin growth model. International Journal of Solids and Structures. 2020; 198: 87–98.

21.

von Streng

Abi-Akl

Giovanardi

Cohen

. Morphogenesis and proportionate growth: A finite element investigation of surface growth with coupled diffusion. Journal of the Mechanics and Physics of Solids. 2021; 146: 104211.

22.

Zurlo

Truskinovsky

. Printing non-Euclidean solids. Phys Rev Lett 2017; 119(4): 048001.

23.

Truskinovsky

Zurlo

. Nonlinear elasticity of incompatible surface growth. Phys Rev E 2019; 99(5): 053001.

24.

Naghibzadeh

Walkington

Dayal

. Surface growth in deformable solids using an Eulerian formulation. J Mech Phys Solids 2021; 154: 104499.

25.

Clayton

. Nonlinear Eulerian thermoelasticity for anisotropic crystals. J Mech Phys Solids 2013; 61(10): 1983–2014.

26.

Clayton

. Nonlinear elastic and inelastic models for shock compression of crystalline solids (Shock Wave and High Pressure Phenomena). Cham: Springer, 2019.

27.

Freed

Zamani

Szabo

, et al. Laplace stretch: Eulerian and Lagrangian formulations. Z Angew Math Phys 2020; 71(5): 157.

28.

Liu

Walkington

. An Eulerian description of fluids containing visco-elastic particles. Arch Ration Mech Anal 2001; 159(3): 229–252.

29.

Kamrin

Nave

. An Eulerian approach to the simulation of deformable solids: Application to finite-strain elasticity. arXiv 2009; arXiv:09013799.

30.

Kamrin

Rycroft

Nave

. Reference map technique for finite-strain elasticity and fluid–solid interaction. J Mech Phys Solids 2012; 60(11): 1952–1969.

31.

Pantaloni

Le Clainche

Carlier

. Mechanism of actin-based motility. Science 2001; 292(5521): 1502–1506.

32.

Noireaux

Golsteyn

Friederich

, et al. Growing an actin gel on spherical surfaces. Biophys J 2000; 78(3): 1643–1654.

33.

Yagoda

. Resolution of a core problem in wound rolls. J Appl Mech 1980; 4: 847–854.

34.

Konig

McAllister

Dusserre

, et al. Mechanical properties of completely autologous human tissue engineered blood vessels compared to human saphenous vein and mammary artery. Biomaterials 2009; 30(8): 1542–1550.

35.

Trangenstein

. Numerical solution of hyperbolic partial differential equations. Cambridge: Cambridge University Press, 2009.

36.

Sugiyama

Takeuchi

, et al. A full Eulerian finite difference approach for solving fluid–structure coupling problems. J Comput Phys 2011; 230(3): 596–627.

37.

Zurlo

Truskinovsky

. Inelastic surface growth. Mech Res Commun 2018; 93: 174–179.

38.

Swain

Gupta

. Biological growth in bodies with incoherent interfaces. Proc R Soc London, Ser A 2018; 474(2209): 20170716.

39.

Rodriguez

Hoger

McCulloch

. Stress-dependent finite growth in soft elastic tissues. J Biomech 1994; 27(4): 455–467.

40.

Chenchiah

Shipman

. An energy-deformation decomposition for morphoelasticity. J Mech Phys Solids 2014; 67: 15–39.

41.

Reina

Conti

. Incompressible inelasticity as an essential ingredient for the validity of the kinematic decomposition F = F^eFⁱ. J Mech Phys Solids 2017; 107: 322–342.

42.

Garikipati

. The kinematics of biological growth. Appl Mech Revs 2009; 62(3): 030801.

43.

Yavari

. A geometric theory of growth mechanics. J Nonlinear Sci 2010; 20(6): 781–830.

44.

Humphrey

Rajagopal

. A constrained mixture model for growth and remodeling of soft tissues. Math Models Methods Appl Sci 2002; 12(03): 407–430.

45.

Sadik

Yavari

. On the origins of the idea of the multiplicative decomposition of the deformation gradient. Math Mech Solids 2017; 22(4): 771–772.