A finite element method for modeling surface growth and resorption of deformable solids

This work explores a continuum-mechanical model for a body simultaneously undergoing finite deformation and surface growth/resorption. This is accomplished by defining the kinematics as well as the set of material points that constitute the domain of a physical body at a given time in terms of an evolving reference configuration. The implications of spatial and temporal discretization are discussed, and an extension of the Arbitrary Lagrangian–Eulerian finite element method is proposed to enforce the resulting balance laws on the grown/resorbed body in two spatial dimensions. Representative numerical examples are presented to highlight the predictive capabilities of the model and the numerical properties of the proposed solution method.

the motility of certain classes of cells [14], the erosion of tissue [9], as well as engineering processes, such as electrodeposition [1], epitaxial growth of thin films [30], construction of built-up systems [7], metal solidification [24], and additive manufacturing [31]. In contrast to volumetric growth, where addition of mass is manifested through a global balance law that modifies the density of material points, surface accretion or resorption uniquely defines the evolution of the boundary in terms of the addition or removal of layers of mass [12].
Unlike elastically loaded bodies, the stresses and deformations that a growing body experiences are inherently tied to the manner in which material is deposited. In one of the earliest works to examine this dependence [5], where the authors analyzed the stress and displacements of a linearly elastic sphere subjected to an incrementally increasing self-weight as material is added onto its surface. The history-dependence of the fields on the accreted sphere is a feature which inherently arises when attaching layers of material which must equilibrate with the previously loaded body without introducing any material gaps or overlaps. The fact that material continuity is maintained throughout the accretion process implies that the unloaded sphere contains residual stresses, as schematically illustrated in Fig. 1.
To extend the notion of surface growth/resorption to bodies undergoing finite deformation, the concept of a mapping between a material and spatial configuration must be re- interpreted since here material constantly enters or exits the body through the exterior boundary. A local reference state is defined in [18] for each deposited material point based on its initial time of existence and coordinates of attachment. Therefore, the growth surface is a collection of points that each have a unique spatio-temporal marker which identifies its initial state. The concept of a reference configuration for a growing body is extended in [19] by defining an associated manifold as an arbitrary abstract representation of both grown and ungrown material at any given time in the three-dimensional Euclidean point space. The motion of the growing and deforming body is thus characterized by a homeomorphic transformation which maps points in the associated manifold to the physical space. In [3], the notion of a oneto-one mapping between a deforming body and multiple generations of unique reference configurations is formulated and used to simulate interstitial growth of biological tissues. Similarly, in [15] the deformation of a growing body is described relative to evolving "natural" configurations which are defined for a given point by a locally stress-free state. An intermediate configuration is introduced in [13] by way of a surface growth/resorption transformation which generates or removes material from a deformed configuration. A separate mapping relates material points from a given reference state with an associated reference time onto the intermediate configuration at the current time. In this setting, the deformation and surface/growth are effectively decoupled between these two times. Recent works [16,17] formalize the notion of an evolving material manifold based on a non-Euclidean material connection. In the same spirit, a material manifold is introduced in [28,29] endowed with a growth-dependent metric tensor that characterizes the coupling between surface accretion and deformation.
Although a a computationally tractable kinematic description of surface growth/resorption for deformable bodies remains largely elusive, numerous recent works have attempted to address the computational challenges in tracking material and its history as it enters or leaves the body. In [1,25,30], numerical methods are developed that utilize an Eulerian grid to solve a global set of boundary-value problems while tracking the growth and resorption front using level-set methods. Other notable works that address different aspects of the computational modeling of surface growth/resorption include [2,[7][8][9][10][11]22,28]. Of particular importance to the method developed in this paper is the work in [13], which introduces a kinematic description of discrete accretion and resorption increments on the surface of deformable bodies based on the Arbitrary Lagrangian-Eulerian (ALE) finite element method, as described, e.g., in [6]. This technique is implemented in [14,23] to simulate a one-dimensional treadmilling cell. A critical advantage of ALE is that it involves a finite element mesh with fixed connectivity whose boundary can naturally define the interface of the growing/resorpting domain. Consequently, ALE does not require additional algorithmic infrastructure to track a moving interface on a fixed background grid such as node activation/deactivation, and the interpolation of boundary conditions between elements.
In the formulation and implementation presented here, "spurts" of growth/resorption are applied instantaneously in time onto the deforming body at prescribed intervals. The deformation of the body at a given time is tracked relative to an updated configuration (akin to the updated Lagrangian finite element method) which contains the set of material points existing after the prior growth/resorption spurt. Imposing balance of linear momentum thus requires keeping track of the displacement relative to this updated configuration, rather than tracking the material displacement of each point relative to its undeformed state. Within this framework, an ALE-like method is introduced to accommodate the boundary motion induced by surface growth/resorption on a domain without the need to constantly remesh. Interior elements deform based on a separate balance law irrespective of locations of prior growth interfaces, thus eliminating restrictions on the geometry of the domain that exist in the prior implementations [13,14]. This allows for numerical modeling of arbitrarily complex surface growth/resorption in conjunction with a Lagrangian description of a body undergoing finite deformation.
In summary, the primary focus of this paper is threefold: (a) to formulate a general kinematic description of a body undergoing concurrent surface growth/resorption and finite deformation, (b) to introduce a discretization strategy that obtains approximate solutions to the physical balance laws for a two-dimensional body, and, (c) to quantify the effect of surface growth/resorption on the numerical accuracy and robustness for the proposed algorithm. The organization of the article is as follows: In Sect. 2, the governing balance laws are derived considering an arbitrary growth/resorption increment. The space and time discretization of the balance laws are presented in Sect. 3, while the numerical implementation of discrete surface growth and resorption is introduced in Sect. 4. Algorithmic and physical properties of discrete growth/resorption are explored by means of two examples in Sect. 5. Lastly, concluding remarks are presented in Sect. 6.

Kinematic description
In conventional continuum mechanics, the Lagrangian description of a deformable body is rooted in the notion of an isomorphism between an initial configuration and the current configuration. The definitions of these configurations give meaning to kinematic quantities, such as the deformation gradient, used to determine the stress for, say, hyperelastic materials or, more generally, simple materials in the sense of Noll [20]. For a deforming body undergoing surface growth/resorption, the conventional definition of an initial configuration must be reinterpreted since new material may enter or existing material may exit through its boundary. In this work, the kinematics is presented for a body undergoing "spurts" of surface growth/resorption (applied instantaneously at discrete times) while continuously deforming.
Consider the domain R t 0 of a body B embedded in the three-dimensional Euclidean point space E 3 at some fixed initial time t 0 . In the present setting, it is assumed, for simplicity, that this body is entirely stress-free and deformation-free at t 0 . Suppose the body deforms continuously in time while undergoing a series of surface growth and resorption spurts. These spurts occur simultaneously across the entire boundary at discrete times t i , i = 1, 2, . . . Any subsequent growth occurs on top of the outermost layer of current growth regions, hence there is no interior (volumetric) growth. Let the domain of the body at the current time t be denoted R t . It is assumed that the body is in existence throughout the time interval [t 0 , t] to exclude the possibility that it "grows out of nothing" or that it is fully resorbed to where no material exists. Explicit reference to the initial time t 0 is typically omitted hereafter for brevity, with the assumption that all elapsed times are expressed with respect to t 0 .

Description of growth/resorption
Consider the i-th growth/resorption spurt occurring at some elapsed time t i . The domain of the body just before the spurt at t i is denoted R t i , with boundary ∂R t i . The transformation of the boundary relative to a fixed frame of reference at the instant that material is grown or resorbed at t i is characterized by the bijective mapping χ t i g : ∂R t i → ∂R t i , as in Fig. 2. The new surface ∂R t i encloses a volumeR t i termed the updated configuration at time t i . The primary purpose of this updated configuration is to track the set of material points existing after the application of surface growth/resorption at t i .
The parts of the boundary ∂R t i that undergo surface growth and resorption are denoted t i g and t i a , respectively. In addition, the part of the boundary undergoing no surface growth/resorption is denoted t i . Since this last part of the boundary remains unchanged by the spurt at t i , χ t i g maps points on this surface to their identical positions. It is assumed here that t i a and t i g are both smooth and have a unique orientation at each point. In view of the preceding definitions, the boundary at t i is formally defined as ∂R t i = t i g t i a t i , where denotes disjoint union.
As a consequence of the surface transformation, the interior of the updated configuration at t i contains a subset of material points that did not exist in R t i , which constitute the domain of the growth region G t i , see Fig. 2. The subset of points that were contained within the interior of R t i and are now outside the boundary ofR t i defines the resorption domain A t i . Also, the part of R t i that remains unaltered by χ t i g is denoted M t i . The surface t i g separates the unaltered domain M t i from the growth domain G t i . The new growth and resorption boundaries generated by χ t i g are denoted˜ t i g and˜ t i a , respectively. The portion of the boundary which remains unaltered is simply denoted˜ t i . Hence, the updated configuration at t i is formally defined as Fig. 2 Kinematics of surface growth/resorption and its boundary is defined as The transformation χ t i g can be determined by the imposed growth and resorption displacements u t i g and u t i a (respectively) of the boundaries˜ t i g and˜ t i a relative to their positions just prior to the i-th spurt as where X t i and x t i represent the initial and current position vector of a given material point (respectively) relative to a fixed basis. Growth and resorption can be formally defined by placing the conditions on the displacements of t i g and t i a where n ∂R t i is the unit outward-facing normal on ∂R t i .

Description of deformation and material motion
The configuration R t i contains regions which may have undergone deformation prior to the addition/removal of material at t i . However, the subsequent accretion on the boundary of R t i will always result in a continuous body which is free of gaps and/or material overlaps. Suppose each material point in the newly grown body (which is represented as the updated configurationR t i ) can be uniquely mapped to a configuration R t i which is free of deformation and stress .
This configuration will contain regions which are geometrically incompatible at the growth interface t i g , as well as other prior growth interfaces. It is reasonable to assume that each subregion within R t i contains material as it appears when it initially attaches to the body. Therefore, R t i is defined as where G t j ;t i denotes the part of G t j which continues to be in existence at time t i (≥ t j ). Also, by definition, G t 0 ;t i represents the portion of the original body that remains in existence at t i . The subscript j will be used hereafter to index each growth spurt at or prior to the i-th spurt. The existence of these subregions which have geometrically incompatible deformation-free initial configurations arises as a direct consequence of the fact that accretion is discontinuous in time. Therefore, material which is attached in spurts along the boundary of a continuously deforming body will depend only on position of the boundary at the instant of accretion. This essentially means that material points can have initial positions which overlap or have gaps in space, yet are distinguished by the different times at which they are grown onto the body. The geometric incompatibility between each of the accreting layers leads to residual stresses, as will be illustrated in the examples of Sect. 5.
Next, define the total deformation of the j-th growth incre- is the image of the initially stress-free domain of material which came into existence at t j mapped onto the continuous and deformed body R t in its current state, as shown in Fig. 3. The two-term superscript in χ t j ;t d and R t j ;t is intended to emphasize the dependence of this transformation and subregion at its current state t on its initial configuration at the time of attachment t j . It can be reasonably assumed that this mapping is bijective since the configurations G t j ;t i and R t j ;t are both occupied at different times and share the same material points.
The deformation mapping χ t j ;t d ensures that the body remain free of gaps and overlaps. This constraint of material continuity takes the form, where denotes the difference between the value of χ d on either side of the interface t j g . This constraint applies for all points on t j g and for all times at or after t j . The total material displacement at t ∈ [t i , t i+1 ) is defined for the j-th growth increment as In Eq. (7), the initial positions X t j are assumed to be fixed in time once the material comes into existence at that location. The total material displacement as defined in Eq. (7) is discontinuous across each growth interface when the body deforms between growth spurts due to the fact that its initial configuration has gaps and/or overlaps. The fact that each growth spurt has an initial configuration which is not continuously connected to any other spurt has important consequences when solving the fully discretized set of balance laws on the growing/resorbing body, which is discussed in R t j ;t and takes the form Note that there are no continuity requirements on the gradient of χ Fig. 3 Schematic of the geometrically incompatible stress-and deformation-free configurations R t2 (each subdomain's time slab shown at bottom inset) and the simply connected current body in its deformed state R t for a two-dimensional hollow ellipse after undergoing two accretion spurts g tinuous along each growth interface. Such a scenario will be explored in Sect. 5.2.
The material velocity at time t ∈ [t i , t i+1 ) is defined for a given material point accreting at t j as where d j /dt denotes the material time derivative (that is, keeping X t j fixed). It is assumed that the velocity is continuous along each growth interface, namely, An alternative kinematic description of the deformed body at the current time t can be constructed based on the updated configuration at t i , assuming that both R t andR t i share the same material points. To this end, the incremental deformation is defined as the transformation wherex t i are the coordinates in the updated configuration at t i . The kinematic constraints on the incremental deformation are expressed as Additionally, the incremental velocity defined as is continuous along the i-th growth interface, i.e., The incremental deformation gradient is defined for the entire updated configuration at t i as the mapping F t i ;t : R t . This two-point tensor can be computed based on the preceding incremental displacements as where I is the two-point identity tensor. Assuming the jth growth increment is deformed at t i , the total deformation gradient for some time t ∈ [t i , t i+1 ) is expressed in terms of the incremental deformation gradient above as This incremental approach defined above closely resembles the updated Lagrangian finite element method [4] and will be used here to simulate surface growth/resorption. To summarize, the kinematics of body undergoing finite deformation and spurts of surface growth/resorption can be described with respect to either its individual deformationand stress-free subregions which together comprise the body as a whole as in Eqs. (7) and (8) (akin to a total Lagrangian approach), or to an updated configuration containing all of the material points existing at a given time as in Eqs. (11) and (15) (akin to an updated Lagrangian approach). The former approach is useful in characterizing the growing body relative to the multiple geometrically incompatible reference configurations of each accreted layer, which ultimately leads to residual stress. The latter approach is useful as a tool to simplify the tracking of each growth/resorption interface in computation.
The relevant transformations and domains of interest described in this section are succinctly summarized in Tables 1 and 2 below.
The configuration of the body at time t i just prior to the application of the i-th growth/resorption increment ∂R t i The boundary of the body at time t i just prior to the application of the i-th growth/resorption increment The updated configuration at the i-th growth/resorption increment which connects the deformed and resorbed body as it appears at t i (M t i ) to the grown region G t i ∂R t i The boundary of the updated configuration at the i-th growth/resorption increment which includes the newly generated surface upon instantaneous growth/resorption (˜ t i g and˜ t i a respectively), as well as the portion of the boundary which does not grow or resorb ( t i ) The set of geometrically incompatible stress-and deformation-free configurations just after the i-th growth/resorption increment (equivalent to the disjoint union of G t j ;t i for j = 0 to i) which together define all the material existing at t i R t The configuration of the body at time t, which for i prior growth increments, consists of the union of R t j ;t for j = 0 to i The i-th resorption surface on the updated configuratioñ The resorbed domain at the i-th resorption increment The grown domain at the i-th growth increment The portion of the j-th growth increment remaining at t i in its initial configuration The subdomain ofR t i which contains material existing prior to t i Generates the i-th growth/resorption domain by mapping the boundary ∂R t i to the updated con- Defines the deformation of the j-th growth increment in its occupied subregion of the current configuration R t j ;t relative to its deformation-and stress-free state G t j (assuming j ≤ i for i growth increments) Defines the deformation of the i-th growth increment in the current configuration R t relative to the updated configurationR t i

Balance laws
The body described in Sect. 2.1 is composed of subregions, each of which comes into existence at different times. The subregions themselves are constrained to form a continuous body pieced together along each of the discrete growth interfaces. In addition to these continuity constraints, each subregion is individually subject to physical balance laws as described below.
The mass density ρ t j 0 is defined as a local measure at a given initial position X t j for a given stress-free region G t j ;t i through the limit where and vol(P t j δ ) is the volume of material enclosed by the sphere P t j δ ⊂ E 3 centered at X t j with radius δ > 0. The mass density in Eq. (17) characterizes the mass per initial volume of a given material point in G t j ;t i . Note that ρ t j 0 is not well-defined on the growth boundary t j g since its value can be discontinuous along this interface.
Mass is conserved in each individual subregion after it attaches onto the body, hence where ρ is the mass density in the current configuration and J t j ;t = det(F t j ;t ). The total mass gained/lost due to the i-th growth/resorption spurt is expressed as Assuming quasi-static loading conditions, the point-wise equilibrium equation is expressed in the current configuration subject to boundary conditions In Eqs. (22) and (21), T t j ;t is the Cauchy stress tensor and b is the body force per unit mass. Also, the divergence operator in (21) is taken with respect to the coordinates of the body in its current/spatial configuration. Additionally, the a-th component of the displacement or traction is prescribed on t a and u b , respectively. These subdomains satisfy the classical requirements In addition, the constraint enforces equilibrium on each of the growth interfaces (with j = 1, 2, . . . , i). The kinematic conditions (12) on the deformation and (14) on material velocity, together with the continuity-of-traction condition (24), are the three constraints on each growth interface required to achieve a unique and well-posed incremental displacement field that satisfies the equilibrium equations (21) with the proper boundary conditions in Eq. (22). It is assumed that angular momentum holds at a given material point which implies symmetry of the Cauchy stress tensor, i.e., Moreover, energy balance under purely mechanical conditions leads to a relation between the Cauchy stress T t j ;t and the internal energy density at time t, which is stated as For a hyperelastic material, Eq. (26) reduces to the familiar relation between the Cauchy stress and deformation gradient, In addition to the balance laws and boundary/initial conditions in Eqs. (21) to (23) and (26) to (22), the surface growth/resorption displacement is imposed on the portions of the boundary where material enters or leaves the domain, according to Eq. (3) for the most recent growth increment at time t i . Material which accretes at t i must have a prescribed initial state on the entire growth domain G t i when it first comes into existence. Here, it is assumed for simplicity that the accreting material does not inherit any history of stress/deformation, hence the conditions hold for all material points in G t i . The subscript "0" emphasizes the initial values of the displacement, deformation gradient, and Cauchy stress in Eq. (28), and not their values once the accreting material equilibrates with the existing body. It is assumed for simplicity that the growth/resorption surfaces remain disjoint to the Neumann and Dirichlet boundaries. This condition can be expressed as g a ∩ u t = ∅. The formal initial/boundary-value problem for a deforming body undergoing surface/growth resorption thus consists of finding the incremental material displacement u t i ;t which satisfies equilibrium as stated in Eq. (21), the boundary conditions in Eq. (22), and the constraints given by Eqs. (12), (14) and (24) given a constitutive law for the Cauchy stress which satisfies Eqs. (27) and (25), an imposed growth/resorption displacement as in Eq. (3), and growth extensions satisfying Eq. (28).

Finite element approximation
In this section, the discretized weak forms of the governing balance laws are presented in two spatial dimensions. To this end, let the closure of the current configuration R t be approximated by finite elements with domains t,e , such that where the superscript e denotes a single two-dimensional element subdomain. At any time t ∈ [t i , t i+1 ), the classical Bubnov-Galerkin approximation for the incremental material displacement u t i ;t,e where nen denotes the number of nodes per element. In addition, û t i ;t,e andξ e are the nodal values of the respective fields ordered in vector form.
Also, the element interpolation matrix is classically defined as in terms of the standard element interpolation functions N e n (x), n = 1, 2, . . . , nen.
The spaces of admissible functions are In Eq. (32), H 1 ( t,e ) denotes the Sobolev space of order 1 over all two-dimensional vectors functions in t,e . Using the interpolated variables, the discrete weak form of the equilibrium equations for element e is The global counterpart of Eq. (34) is computed by the standard assembly procedure.

Discretization of surface growth and resorption
Consider a representative time t i at which growth/resorption takes place across the boundary. It is assumed here, without substantial loss of generality that the incremental displacement fields are computed at this time relative to the configuration at time t i−1 using Eq. (11), such that they satisfy the global counterpart of the equilibrium equations (33). These are used to define the deformation gradient and Cauchy stress according to Eqs. (15) and (27), respectively. The incremental displacements, deformation gradient, and Cauchy stress at t i pertain to the deformed body t i just before the i-th growth/resorption spurt takes place.
To effect growth/resorption at time t i , the regions G t i and A t i are first created based respectively on the surface growth/resorption displacements u g and u a of Eq. (3), as sampled at t i . In particular, the nodal points lying on the original boundary are placed at the positions of the newly grown/resorbed boundary, while the placement of interior nodal points is optimized to maintain a desirable quality. This optimization consists of solving the discretized equilibrium Eq. (33) for the incremental mesh displacements u m rather than the incremental material displacements. In this setting, the boundary mesh displaces by its prescribed growth/resorption (the portion of the boundary not undergoing surface growth/resorption remains fixed along its normal direction), while the interior mesh deforms as an imaginary solid with an assumed constitutive law and user-specified mesh material parameters (in this setting, linear elasticity), as shown in Fig. 4. This mesh motion implementation closely resembles procedures commonly used in production codes, e.g., [21]. The mesh optimization strategy does not constrain element edges to align with prior growth interfaces, thus allowing the element size and shape to be independent of the growth/resorption rate and time increment between spurts. Nodes that are "flagged" for lying in the growth region G t i are initialized with a zero-valued incremental material displacement. Additionally, integration points that lie inside G t i are assigned a deformation gradient equal to the identity tensor [according to Eq. (28)]. On the other hand, the deformation gradient in M t i is initialized based on its known values at t i prior to growth/resorption. This is accomplished through a global L 2 -projection which maps each component of the deformation gradient from the mesh defining R t i onto the subdomain M t i of the newly meshed and updated configurationR t i . The Cauchy stress is computed and stored for the updated configurationR t i , based on the constitutive law in Eqs. (25) and (27) as well as the initial Fig. 4 Schematic of the boundary-value problem used to determine the mesh motion u t i m based on the imposed surface growth and resorption displacement at time t i condition (28) in the growth region. The procedure highlighted above is related to the arbitrary Lagrangian-Eulerian (ALE) finite element method [6], whereby the mesh motion is prescribed rather than deforming with the body at each point (Lagrangian) or being fixed (Eulerian). However, unlike traditional ALE methods the mesh boundary moves along the growth/resorption front, thus defining a constantly changing physical domain of material. The methodology described above is summarized in Algorithm 1. In the approach described here, the choice is made to solve for the incremental material displacement in Eq. (11), which tracks material motion relative to the updated configuration at t i rather than the total material displacement [Eq. (7)] which accounts for the entire motion of each growth spurt relative to its undeformed state. In this regard, the proposed method is also related to the classical updated Lagrangian formulation [4]. The dependence on the updated configuration rather than a globally deformation-free (geometrically incompatible) configuration only requires knowledge of the most recent growth/resorption interface. This eliminates the need to store and update the locations of all prior growth interfaces as they deform with the body. Additionally, the placement of the newest accretion increment onto the deformed body forming the updated configuration at a given time t i naturally results in a fully continuous body, hence Eqs. (12) and (14) are automatically satisfied as long the subsequent incremental deformation and incremental velocity after t i are continuous. By extension, Eqs. (6) and (10) are also satisfied.

Algorithm 1 Surface growth/resorption algorithm for quasistatic (non-inertial) bodies with no body forces
The realignment of element edges after each growth/ resorption spurt occurs independently of prior edge locations. Therefore, discontinuities in the material displacement and deformation gradient at prior growth interfaces must be approximated as piecewise polynomials constrained by the newly aligned element edges, as shown in Fig. 5. The examples presented in the following section illustrate that the errors associated with this smearing process decrease as the element size and the time-step are reduced.

Representative simulations
In this section, two idealized examples are presented under the assumption of plane strain. These examples demonstrate the numerical consistency and spatial/temporal convergence of the proposed surface growth/resorption algorithm. The first concerns an ellipsoidal cylinder undergoing simultaneous surface growth/resorption and rigid body motion, while the second involves a hollow elliptical cylinder simultaneously undergoing finite elastic deformation due to an imposed time-dependent traction on its inner boundary and surface growth/resorption on its outer boundary. The material response is assumed compressible neo-Hookean, such that for a given material point which comes to existence at time t j , the Cauchy stress at time t ∈ [t i , t i+1 ) is expressed as where B t j ;t = F t j ;t (F t j ;t ) T is the left Cauchy-Green tensor, i is the spatial rank-two identity tensor, and λ, μ are material parameters akin to the Lamé constants of linear elasticity with values selected as 7.1E+02 MPa and 1.8E+02 MPa, respectively, as is typical for various plastics. With these parameters, the material is compressible. All of the meshes used in these examples consist of bilinear quadrilaterals, and numerical integration of the weak forms is performed using standard 2 × 2 Gauss-Legendre quadrature.

Surface growth and resorption under rigid-body motion
The initial domain in this example is an ellipse with diameters of 80 mm on the major axis and 40 mm on the minor axis. The discretized domain consists of 1044 elements, as shown in Fig. 6. Two cases are examined with imposed growth/resorption and incremental material displacements, as depicted in Fig. 8. For both cases, growth spurts are assumed to occur every t = 1 s for t ∈ (0, 20] s. In Case 1, a time-dependent Fig. 7 Mesh motion of the ellipsoidal geometry for Case 1 (left) and Case 2 (right) at t = t 3 . Initial and newly grown/resorpted mesh are shown in gray and black, respectively (the motion due to growth/resorption is scaled for clarity) uniform normal growth/resorption displacement is imposed, such that for the i-th growth spurt, Case 2 involves constant growth u t i gy = 1.0e−3 m on the top boundary and, likewise, constant resorption u t i ay = 1.0e−3 m on the bottom boundary for each spurt, as well as a constant incremental material displacement u y = 5.0e−4 m. In this setting, the assumed boundary displacements due to growth/resorption preserve the initial geometry of the ellipse, resulting in a rigid "treadmilling" motion.
In both cases, the elliptical cylinder experiences addition or removal of material without experiencing deformation. In the first case, the domain expands and contracts as the surface first grows and subsequently resorbs. In the second case, material is deposited onto the top-half of the boundary and simultaneously removed from the bottom-half, while the body is treadmilling with a constant velocity in the ydirection. The Lamé constants used to solve for the mesh motion are λ = μ = 1.0E+8 MPa. This motion is illustrated at a sample time in Fig. 7. In both cases, the entire body remains stress-free for the duration of the simulation, as highlighted in Fig. 8. These simple test cases confirm that for an initially unstressed body experiencing surface growth/resorption, elastic deformation is solely dependent on externally applied forces, and is not generated by stressfree surface growth or resorption alone. In this regard, these two cases may be thought of as providing a test of consistency for the algorithmic implementation.

Growth and resorption of elliptical cylinder undergoing finite deformation
Here, results are presented for a hollow elliptical cylinder experiencing a series of constant-magnitude normal growth/resorption spurts along its outer boundary and a temporally increasing pressure along its inner boundary. In the reference configuration, the cross-section of the cylinder shows outer major and minor axes of 2 m and 3 m, respectively, and an inner radius of 1.5 m. The pressure is prescribed along the inner radius as p(t) = 10 7 t Pa, within the range t = (0, 1.0] s, while the outer surface remains traction-free. Two cases are considered with regard to growth. In Case 1, the body undergoes surface growth , while in Case 2, it experiences resorption. The geometry and loading for the two cases are shown in Fig. 9. Appealing to double symmetry, only one quadrant of an ellipse is modeled . The mesh consists of 28,441 elements and is gradually refined near the outer surface exposed to growth or resorption. Additionally, a constant t = 0.025 s is used in the analysis and the prescribed normal growth and resorption displacements for the i-th spurt are 0.5 t m (Case 1) and −0.1 t m (Case 2). The effects of surface growth and resorption on the final state of the deformation and stress are highlighted in Fig. 10 at the end-time t = 1 s . The accretion of material leads to an lower overall von Mises stress and absolute pressure defined as | p| = 1/3|tr(T )|. The results in Fig. 10 also indicate that the majority of the internal force is sustained by the original material with relatively little redistribution in the growth region. The abrupt change in the von Mises stress and pressure is most pronounced along the minor radii, where steep gradients between the original boundary and the accreted regions exist. No material is added to the cylinder in the second case where it only undergoes resorption. The points of maximum and minimum stress therefore remain on the outer major and minor axes (respectively) throughout the resorption process. The spike in the von Mises stress and absolute pressure for both the grown and resorbed cylinder occurs at the top/bottom since these are the thinnest regions of the body.
The material is assumed to be deposited at each discrete time of the analysis, hence the time-step size dictates the magnitude of the growth displacement. Unsurprisingly, selecting a coarser step-size, here t = 0.1 s, leads to visibly discrete stress discontinuity. This is illustrated in Fig. 11 for the circumferential, radial, and longitudinal components of the Cauchy stress along a vertical and horizontal section of the ellipse at t = 1 s. The observed "stair-stepping" pattern in the stress occurs due to the large differences in certain components of stress and deformation between the pre-existing and newly added material stress or deformation. Note that these discontinuities are effectively "smeared" between element edges when prior growth interfaces cut through elements, thus resulting in a steep (continuous) These gradients are more pronounced in the circumferential and longitudinal directions than in the radial direction. In contrast, the radial stress is zero along the outer surface due to the traction free boundary conditions, hence it maintains a smooth profile. Therefore, once each finite layer attaches to the cylinder, it primarily resists subsequent loading by deforming along the circumferential and longitudinal direction. Note that the mismatch in stresses between the grown and ungrown regions is more prominent along the minor rather than the major axis, since the circumferential stresses are greatest in the former region. In addition to the discontinuity at each growth increment, the slope of the circumferential stress within the oldest growth increment nearly matches the region near the original boundary at t = 0 s whereas the newer increments consist of sequentially lower slopes. Here, the slope is the change in a component of stress per distance for a given growth region, discounting the "smoothed" discontinuities. The outer-most layer of material exhibits near-zero slope in the circumferential component of the stress due to the assumed stress-free initial condition applied uniformly within the newly grown material. These trends highlight that the slopes of the circumferential and longitudinal stress across the growth interface nearly match, despite the fact that the stresses themselves are discontinuous.
As a consequence of the history-dependent surface growth/ resorption process, the cylinder maintains residual stresses upon unloading, as highlighted in Fig. 12. These fields are obtained by removing the traction in a single load step and solving Eq.  body is unloaded, thus requiring the radial stress to be zero on these surfaces. However, the radial component of the Cauchy stress increases quadratically through the thickness of the cylinder cross-section, reaching a maximum near the center of the section thickness. The circumferential and longitudinal stresses are tensile along the inner surface and compressive at the outer surface thus implying that the unloaded section along the radius has the tendency to bend.
The numerical results presented in [13] highlighted a similar stair-stepping effect to the one exhibited in Figs. 11 and 12 for one-dimensional spatial domains with the distinction that each growth interface align precisely with element boundaries. In this setting, the sharp discontinuities that appear at each growth interface are exactly captured. Although this distinction enhances the overall accuracy of the solution,  16 Hollow elliptical cylinder: element length scale metrics for meshes used for surface growth (top) and resorption (bottom) to determine spatial refinement convergence rates. Plotted values are averaged over all time steps ( t = 0.1, t = 0.05, and t = 0.025) for clarity the method in [13] does not readily extend to a general two-dimensional body. With the assumption that the discontinuities are approximated by steep gradients, the mesh motion is performed irrespective of prior locations of the growth interface. In this setting, it is assumed that the fields near each growth interface can be accurately represented as the mesh element size approaches zero.
Although the choice of time step affects the size of the growth increments and character of the solution, the discontinuities lose prominence and the solution converges to a single smooth curve as the size of the growth increments become infinitesimal relative to the scale of the physical domain. This convergence is highlighted in Figs. 13 and 14 for three time Uniform h-refinement of the mesh leads to a spatially convergent solution in the norm of the stored elastic energy of the error, as shown in Fig. 15. Solutions from six meshes were obtained with 741, 1711, 3081, 4851, 7021, and 28441 elements. The errors are taken relative to the solution of the finest mesh for each corresponding time step, and are computed based on the difference in displacement gradients. As expected, the spatial convergence rates are higher at smaller time steps (smaller growth increments) for the grown cylinder since the deformation gradient and stress become smoother as the time step decreases. In contrast, the resorpted cylinder exhibits nearly equal errors for a given mesh regardless of time step, thus indicating that temporal discretization does not play a significant role in the spatial convergence rate for the case of surface resorption. Note that in the case of resorption, the error for a given mesh size is nearly identical for each of the selected time-step sizes, therefore only the errors for a single time-step size are visible in Fig. 15. The convergence rate for surface growth increases with a decrease in time-step size, which is likely higher than expected since the errors are taken with respect to the fields computed on a slightly finer mesh (as opposed to an exact analytical solution which may be more accurate, though this is non-trivial to obtain).
The convergence rates for the mesh used in these simulations (Fig. 9) depends on the metrics used to characterize the level of refinement of the spatial discretization. Table 3 contains a comparison of the convergence rates of the grown cylinder using four element size metrics: minimum, average, volume-averaged, and maximum element lengths. The relevant values of the element size metrics are shown in Fig. 16. All length scales are taken from the deformed mesh at the final simulation time. The average element length is approx-imated as h avg= 1 nel where A t f is the area of the meshed ellipse at the final simulation time t f . The volume-averaged element length is taken as where h e is the average element length. The difference in the norm of the stored energy of a neo-Hookean material asymptotically approaches that of a linear elastic material (also referred to as the energy norm) when the deformations are considered infinitesimal. Therefore, the theoretical convergence rate for a uniform grid and bilinear elements is linear. In an average sense, the convergence rates shown in Table 3 tend to the theoretical value as the time discretization is refined.

Conclusion
This article demonstrates that suitably adapted Lagrangian methods can be used to describe the surface growth and resorption of solids within the framework of classical continuum mechanics. This is enabled by the introduction of discrete time-evolving spurts which create the updated reference configuration relative to which one may describe the kinematics and balance laws. The resulting computational problem is tackled readily by an ALE-like method which circumvents the need for constant remeshing. The mesh motion is determined independently of locations of prior growth increments, which prevents significant mesh distortion. As a consequence, elements can contain material that has come into existence at different times. Numerical examples of a rigidly growing and resorbing ellipse and a growing/resorbing elliptical cylinder undergoing finite deformation demonstrate that the proposed algorithm is both robust and convergent. An extension to three-dimensional geometries is feasible and will be pursued in the future.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the