A Mechanical Theory of Growth

Abstract

A theory of growth is developed, utilizing the notion of a directional density function that captures the number and distribution of the material particles and their changes in time. A spatial (or Eulerian) description of kinematics is adopted, and the constitutive theory for a growing body is developed that relates the stress to the directional density function. The equation that governs the evolution of the directional density function is derived. An example of internal surface growth is presented.

1 Introduction

Growth is the process of addition (or resorption) of material particles¹ and change of mass in a body. Its importance in biological systems and many man-made material systems has attracted the attention of the mechanics community, and the mechanics of growth has been an area of intense research over the past half century. Many theories and extensive analyses have been developed and published in a vast literature. Much of the progress up to 1995 can be found in a review article by Taber [16]. One of the most important developments of using continuum mechanics to model growth (Rodriguez, et al. [10]) is based on the multiplicative decomposition of the deformation gradient tensor into an elastic part and a growth part. A review article by Kuhl [8] provides a discussion of the microenvironmental cues that drive changes of the growth component of the deformation gradient. An excellent monograph in mathematics and mechanics of growth by Goriely [5] contains the contributions up to 2017. Further developments include an article by Goodbrake et al. [4] that establishes the existence of an intermediate configuration in the multiplicative decomposition of the deformation gradient as well as an article by Zhuan and Luo [18] that proposes a new growth law in which the growth part of the deformation gradient decomposition depends on fields evaluated in the current configuration.

Classical theories of mechanics have provided powerful tools for studying the mechanical aspects of growth. The physical and physiological process of growth, in turn, presents unique challenges to the theoretical foundations of mechanics, and promotes the growth of mechanics itself. While continuum mechanics has greatly advanced our understanding of growth, some fundamental issues have not been fully resolved.

Almost all existing growth theories use the basic concepts and formulations developed in the classical continuum mechanics for solids. A starting point is the introduction of a fixed reference configuration and a deformation function, which establishes a one-to-one correspondence between the material particles in a deformed configuration and the material particles in the reference configuration. Strictly speaking, such a fixed reference configuration does not exist for a growing material body, of which the number of material particles is ever changing as new material particles are added to the body and/or some existing material particles are removed from the body. While the notions of a fixed reference configuration and a deformation function are effective in describing the deformation of a non-growing material, they cannot truthfully describe a growth process.

While this difficulty may be alleviated, to some extent, for a volumetric growth process in which the material particles in a fixed reference configuration remain dense in a grown state, it is insurmountable for a surface growth in which a material surface may grow into a volume. Here we must emphasize that the term “surface growth” used in the present work is not limited to a material surface growing into another material surface (also termed area growth by some authors), but encompasses a material surface growing into a volume. In this latter case, the changes in geometry and topology make the concepts of a fixed reference configuration and a deformation function completely irrelevant. Indeed, for a volumetric growth with dense material particles, a small material region in a grown state always contains material particles that existed before the growth. In contrast, in a surface growth, the entirety of material particles in a region need not exist in an earlier state. It is therefore not possible to find a surjective mapping from a material region onto the subsequent material regions. There have been attempts to overcome this difficulty. Among others, DiCarlo [2] treats surface growth as an infinitely intense bulk growth confined in a layer of vanishingly small thickness. Sozio and Yavari [15] treat an accreting body, that is formed by surface growth on the boundary of a bulk material, as a material manifold with a time of attachment map, that serves to specify the time at which a material layer is accreted.

Another issue in the existing growth theories pertains to the concept of a growth tensor. Rodriguez, et al. [10], proposed a multiplicative decomposition of the deformation gradient into an elastic tensor and a growth tensor. This revolutionary idea has led to tremendous progress and success in developing constitutive theories for growing elastic materials. Nevertheless, as the deformation gradient tensor is the gradient of the deformation function, the said decomposition lacks a foundation if a deformation function cannot be defined in the first place. Even if a deformation gradient can be properly defined, the concept of a growth tensor appears to be peculiar in describing a growth process, which embodies the change of the number of material particles. A scalar quantity such as number density seems more appropriate.

In the present work, we formulate a growth theory that does not require a fixed reference configuration or a deformation function. For kinematics, we use the common formulation in fluid mechanics where it is difficult, if not impossible, to establish a one-to-one correspondence between the material particles in two different states. Instead, the motion of a fluid is completely described by the velocity field that specifies the velocity of individual material particles that happen to occupy a place at a time. It is also noteworthy that a conventional fluid has no naturally identifiable reference configuration.

To describe anisotropic growth, we introduce the concept of a scalar-valued directional density function which gives the number density of material particles in any assigned direction. A growth is subsequently described by the rate of change of the directional density function.

The constitutive theory for a growing/deforming material is developed by using the value of the directional density function in the current configuration as the primary constitutive variable. This completely eliminates the need for a fixed reference configuration and is thus suitable for describing the mechanical behavior of a material undergoing growth, especially surface growth. We also derive the evolution equation for the directional density function in a simultaneous growth and deformation process.

We note that there have been considerable efforts in developing the Eulerian (or spatial) description of various theories in continuum mechanics, as exemplified by the works of Eckart [3], Leonov [9], Rubin [11], Rubin and Attia [12], Rubin et al. [13], and Safadi and Rubin [14]. The independent variables in the Eulerian description are space and time, as opposed to (in the terminology adopted here) material points and time in the Lagrangian (or referential) description. However, the formulation of the constitutive theories of the above works still requires a reference configuration. For example, for Eckart [3], the elastic stress is a function of a “strain metric tensor” that specifies the state of strain from a “relaxed” state. In other works, the stress is a function of some kinematic quantities that are determined by integrating certain evolution equations from a reference state to the current state. The present work is aimed toward developing a constitutive theory that does not require a fixed reference configuration and, thus, is completely free of the need to connect the current configuration to such a reference configuration.

In Sect. 2, we discuss the kinematics and some primary variables in the present theory, including the velocity, the mass density, and the stress. In Sect. 3, we introduce the concept of directional number density function and growth rate function. In Sect. 4, we derive the equation of balance of mass, which, in comparison with the equation of conservation of mass for non-growing materials, has an additional term resulted from growth. In Sect. 5, we derive the equation of balance of linear momentum for a growing material. Again, an additional term emerges due to growth. Sect. 6 is devoted to developing the constitutive equations for materials that can undergo growth and deformation simultaneously. We first show that for elastic materials, the stress in the current state depends on the directional density function in a particular manner. We then generalize this result to a larger family of materials by assuming that the stress in the current state depends on the directional density function in a general manner. This leads to a constitutive theory that is completely different in nature from those encountered in classical continuum mechanics, including various models used in existing growth theories. An evolution equation for the directional density function in a simultaneous growth and deformation process is also derived. In the concluding Sect. 7, an internal surface growth problem, which serves as an idealized model for the growth of a tree trunk, is discussed. This problem exposes some special features associated with surface growth, such as discontinuous velocity fields and singular growth rate functions. The explicit expressions for the components of the stress in the wood and in the tree bark in terms of the growth rate function are presented.

2 Kinematics

In this work, we consider a body of material that can undergo deformation, as well as growth during which new material particles are added. The process of resorption, in which some existing material particles are removed, can be described as well.

The basic independent variables of various functions in the present work are place $\boldsymbol{x}$ and time $t$. Specifically, $\boldsymbol{x}\in \mathbb{R}^{3}$ is the position of a point in space,² which is not associated with a material particle, although it may be occupied by a material particle at a given instant of time. Three state variables commonly used in fluid mechanics, as well as in the present work, are velocity vector $\boldsymbol{v}(\boldsymbol{x}, t)$, mass density $\rho (\boldsymbol{x}, t)$, and stress tensor $\boldsymbol{T}(\boldsymbol{x}, t)$, which give, respectively, the velocity, the mass density, and the Cauchy stress of a material particle that happens to occupy the place $\boldsymbol{x}$ at the time $t$. These functions need not be smooth. For example, the velocity may be discontinuous for a singular surface growth. However, for ease of analysis, we shall proceed by assuming that all the functions possess smoothness to the extent sufficient to carry out the analysis. For non-smooth functions, the argument presented in this work can be made rigorous by using the theory of generalized functions (or distributions), and the resulting equations are interpreted in the sense of distributions.

In continuum mechanics, the kinematic theory used in this work is often referred to as spatial (or Eulerian) description, in which $\boldsymbol{x}$ and $t$ are taken as the independent variables, and the velocity ${\boldsymbol{v}(\boldsymbol{x}}, t)$ is taken as the primary kinematic quantity. On the other hand, in the referential (or Lagrangian) description, one of the independent variables is $\boldsymbol{X}$ which gives the place of a material point in a reference configuration. The motion of a body in the referential description is given by a deformation function ${\bar{\boldsymbol{x}}(\boldsymbol{X}}, t)$, which gives the place of the material point $\boldsymbol{X}$ at the time $t$. The velocity function $\boldsymbol{v}(\boldsymbol{x},t)$ in the spatial description and the deformation function $\bar{\boldsymbol{x}}(\boldsymbol{X},t)$ in the referential description are related through the equation

$$ \frac{\partial {\bar{\boldsymbol{x}}(\boldsymbol{X}}, t)}{\partial t} = \boldsymbol{v}(\bar{\boldsymbol{x}}(\boldsymbol{X}, t), t). $$

(1)

For a non-growing material, the evolution equation (1), which constitutes a system of first-order ordinary equations, can be used to find the deformation function ${\bar{\boldsymbol{x}}(\boldsymbol{X}}, t)$ from the velocity function $\boldsymbol{v}(\boldsymbol{x}, t)$, with appropriate initial conditions. For a growing material, this process may not be feasible, even when proper jump conditions for $\boldsymbol{v}(\boldsymbol{x}, t)$ are specified. Indeed, as discussed in the introduction, the deformation function $\bar{\boldsymbol{x}}(\boldsymbol{X}, t)$ need not be well-defined for a growing material.

The spatial description of kinematics is commonly used in fluid mechanics, where it is not always possible to know where a fluid particle was in the past, and where the particle will be in the future. Instead, the focus is on the current state of the fluid and the rate of change of the state variables at present without referring to a fixed reference configuration. Of course, for a growing solid, the lack of a fixed reference configuration is due to the understanding that a material particle may not exist in the past. However, the spatial description appears to be suitable for both types of materials.

3 Particle Distribution, Density, Growth

The existing theories of growth in the mechanics literature, for which the referential description is exclusively used, are built on a basic premise of multiplicatively decomposing the deformation gradient tensor into a growth tensor and an elastic deformation tensor. In addition to the unsettled issue of the existence of a fixed reference configuration, the concept of growth tensor obscures the physical nature of growth. Indeed, growth is a process during which material particles are added to the body. An intuitively appealing and natural idea is to use a scalar state variable to describe the distribution of the material particles, and to relate the growth to the change of that state variable.

A biological tissue at the cellular level is a discrete system, and a growth process can be described by specifying the distribution (the number and the locations) of cells and the change of the distribution in time. One of the simplest descriptions of distribution is to specify the volumetric number density (the number per unit volume) of cells. This description, however, cannot account for orientation dependent distributions, and therefore cannot capture an anisotropic growth process for which the change of the number of cells is orientation dependent. To describe anisotropic growth, we introduce the concept of directional number density function, or directional density function for short, $n: \overline{\mathbb{R}^{3}} \times \mathbb{R}^{3} \times \mathbb{R} \rightarrow \mathbb{R}_{+}$, where $\overline{\mathbb{R}^{3}}$ denotes the set of unit vectors, and $\mathbb{R}_{+}$ denotes the set of non-negative real numbers. The value of the directional density function, $n(\boldsymbol{e}, \boldsymbol{x}, t)$, gives the number density of cells, or material particles for a general material system,³ in the direction of $\boldsymbol{e}$ at the place $\boldsymbol{x}$ and the time $t$. Precisely, it gives the number of the material particles, per unit length, that intersect a straight-line segment in the direction of $\boldsymbol{e}$. For illustrative purposes, consider a schematic representation of a block of cells shown in Fig. 1. When taking a cell with the surrounding extracellular matrix as a material particle, the value of the directional density function in the direction of the specific $\boldsymbol{e}$ as shown in the figure is six. For a system with a random distribution of material particles, the concept of the directional density function must be taken in the statistical sense. For a system with a regular distribution of material particles, it is possible to find the explicit expression for the directional density function. For example, the directional density function for a primitive orthorhombic crystal system with mutually orthogonal lattice vectors $\boldsymbol{a}_{1}$, $\boldsymbol{a}_{2}$, and $\boldsymbol{a}_{3}$ is given by

$$ n(\boldsymbol{e}, \boldsymbol{x}, t) = \sqrt{\sum _{i=1}^{3} \left ( \frac{\boldsymbol{e}\cdot \boldsymbol{a}_{i}}{|\boldsymbol{a}_{i}|^{2}} \right )^{2}}. $$

(2)

Fig. 1

Illustration of the directional density function

For a growing material body, the directional density will change in time. We introduce a directional density growth rate function, or growth rate function for short, $g: \overline{\mathbb{R}^{3}} \times \mathbb{R}^{3} \times \mathbb{R} \rightarrow \mathbb{R}$, which, in the absence of deformation, gives the rate of change of the directional density function:

$$ \frac{\partial n(\boldsymbol{e},\boldsymbol{x}, t)}{\partial t} = g( \boldsymbol{e},\boldsymbol{x},t), $$

(3)

or, equivalently,

$$ n(\boldsymbol{e}, \boldsymbol{x}, \tau) = n(\boldsymbol{e}, \boldsymbol{x}, t)+ \int_{t}^{\tau} g(\boldsymbol{e}, \boldsymbol{x}, s)\, d s. $$

(4)

The dimension of $g$ is number per unit length per unit time. A positive value of $g$ indicates growth, while a negative value of $g$ indicates resorption. We note that the value of $g$ depends on the direction vector $\boldsymbol{e}$, in addition to place $\boldsymbol{x}$ and time $t$. It is possible, for given place $\boldsymbol{x}$ and time $t$, to have growth in one direction, and resorption in another direction. We also note that in a controlled process of growth, such as the growth of a crystalline material, the growth rate function may be controlled by a person (or device) that conducts the operation. In a natural growth process, such as the growth of a biological tissue, the growth rate function may be regulated by factors outside the material, such as gene expressions, extrinsic signals, etc. For example, Ambrosi and Mollica [1] have argued, in the context of the growth tensor associated with the multiplicative decomposition of the deformation gradient, that tumor growth strongly depends upon the availability of nutrients and on the presence of chemical signals.

As shall be discussed in Sect. 6, changes of directional density function can also be induced by a deformation, during which the locations of the material particles change. In this case, equations (3) and (4) for the change of the directional density function $n$ need to be properly modified.

The concept of directional number density function leads to other useful quantities. We define the material particle directional diameter $d(\boldsymbol{e}, \boldsymbol{x}, t)$ by

$$ d(\boldsymbol{e}, \boldsymbol{x}, t) \equiv \frac{1}{n(\boldsymbol{e}, \boldsymbol{x}, t)}, $$

(5)

which gives the effective diameter, in the direction of $\boldsymbol{e}$, of a material particle at place $\boldsymbol{x}$ and time $t$. The integral of the cubic power of the material particle directional diameter over the orientation space gives the volume $\nu(\boldsymbol{x}, t)$ of the material particle at place $\boldsymbol{x}$ and time $t$:

$$ \nu(\boldsymbol{x}, t) = \frac{1}{4 \pi} \int _{\omega} d^{3}( \boldsymbol{e}, \boldsymbol{x}, t) \; d \omega (\boldsymbol{e}). $$

(6)

Here, $d \omega$ is the differential angle, and $\omega$ is the angular space (orientation space). In the standard spherical coordinate system, we have

$$ {\boldsymbol{e}} = \sin \theta \cos \phi \,\boldsymbol{e}_{1} + \sin \theta \sin \phi \, \boldsymbol{e}_{2} + \cos \theta \, \boldsymbol{e}_{3}, $$

(7)

$$ d\omega (\boldsymbol{e}) = \sin \theta \, d\theta \, d\phi , $$

(8)

and

$$ \omega = \{(\theta , \phi ) \in \mathbb{R}^{2}: 0 \leq \theta \leq \pi , 0 \leq \phi \leq 2 \pi \}, $$

(9)

$(\boldsymbol{e}_{1}, \boldsymbol{e}_{2}, \boldsymbol{e}_{3})$ being an orthonormal basis. Furthermore, we define the volumetric number density function $n_{\nu}(\boldsymbol{x}, t)$ by

$$ n_{\nu}(\boldsymbol{x}, t) \equiv \frac{1}{\nu(\boldsymbol{x}, t)}, $$

(10)

which gives the number of material particles per unit volume. Due to deformation or growth, the volumetric number density $n_{\nu}$ may change in time, as shall be discussed in the next section.

4 Balance of Mass

For a non-growing material body, the principle of conservation of mass states that the mass of any part of the body remains unchanged during deformation. For a growing body, the above statement is replaced by the balance of mass, which states that the change of mass is balanced by the mass of additional material particles.

For simplicity, we assume that each material particle has the same mass $\mu$. As such, the balance of mass amounts to the balance of the number of material particles. Consider a part of material (a set of material particles) that occupies region $P(t)$ at time $t$. We denote by $N_{P}(t)$ the number of the material particles in $P(t)$. Then,

$$ N_{P}(t) = \int_{P(t)} n_{\nu}(\boldsymbol{x}, t) \, dV. $$

(11)

By the Reynolds’ transport theorem, the rate of change of $N_{P}(t)$ is

$$ N_{P}'(t) = \int_{P(t)} \Big(\frac{\partial n_{\nu}}{\partial t} +\operatorname{div} (n_{\nu} \boldsymbol{v})\Big)\, dV. $$

(12)

Since $P(t)$ is a material part, the number of the material particles in $P(t)$ changes if and only if growth occurs. The rate of change of $N_{P}(t)$ can be calculated from

$$ N_{P}'(t) = \int_{P(t)} g_{\nu}(\boldsymbol{x}, t) \, dV, $$

(13)

where $g_{\nu}$, called volumetric density growth rate function, is the rate of change of the volumetric number density function $n_{\nu}$ in the absence of deformation. Equating (12) and (13), and invoking the localization theorem, we obtain the local form of the balance for the number of particles:

$$ \frac{\partial n_{\nu}}{\partial t} +\operatorname{div} (n_{\nu} \boldsymbol{v})= g_{\nu}. $$

(14)

In the current theory, the directional density growth rate function $g$ is assigned as given a priori information, which allows one to calculate the directional density function $n$ using an extended version of (4), to be derived in Sect. 6. Then, one can determine the volumetric number density function $n_{\nu}$ using (5), (6), and (10). Furthermore, equation (14) serves to determine the volumetric density growth function $g_{\nu}$ to ensure the balance of the number of material particles.

Since the mass density $\rho$ is the mass per unit volume and the volumetric number density $n_{\nu}$ is the number of the material particles per unit volume, we have

$$ \rho(\boldsymbol{x}, t) = \mu\, n_{\nu}(\boldsymbol{x}, t). $$

(15)

Combining (14) and (15) gives the balance of mass:

$$ \rho_{t} +\operatorname{div} (\rho \boldsymbol{v})= \mu \, g_{\nu}. $$

(16)

In the absence of growth, we have $g_{\nu}(\boldsymbol{x}, t) = 0$, and the balance of mass (16) reduces to the usual conservation of mass in continuum mechanics:

$$ \rho_{t} + \operatorname{div}(\rho \boldsymbol{v}) = 0. $$

(17)

5 Balance of Linear Momentum

In continuum mechanics, the balance of linear momentum states that the rate of change of the linear momentum of a part of a material equals the resultant force on the part. The linear momentum of a material particle is the product of its mass and velocity, and the linear momentum of a part of a material is the sum of the linear momentum of all material particles in the part. For a growing material body, the change of linear momentum results from changes of both velocity and number of particles.

We again consider a part of material that occupies region $P(t)$ at time $t$. The linear momentum $\boldsymbol{l}(t)$ of the material part is given by

$$ {\boldsymbol{l}}(t) = \int _{P(t)} \rho (\boldsymbol{x}, t) \boldsymbol{v}(\boldsymbol{x}, t)\, dV. $$

(18)

By the Reynolds’ transport theorem, the rate of change of the linear momentum of the material part is

$$ {\boldsymbol{l}}'(t) = \int _{P(t)} \big( (\rho \boldsymbol{v})_{t} + \operatorname{div} (\rho \boldsymbol{v}\otimes \boldsymbol{v})\big)\, dV. $$

(19)

The resultant force on the material part $P(t)$ is given by

$$ {\boldsymbol{f}}(t) = \int _{P(t)} \rho (\boldsymbol{x}, t) \boldsymbol{b}(x, t)\, dV + \int _{\partial P(t)} \boldsymbol{t}( \boldsymbol{x}, t)\, dA, $$

(20)

where $\boldsymbol{b}$ is the body force per unit mass, and $\boldsymbol{t}$ the surface traction (force per unit surface area). By Cauchy’s representation theorem $\boldsymbol{t}(\boldsymbol{x}, t) = \boldsymbol{T}(\boldsymbol{x}, t) \boldsymbol{n}(\boldsymbol{x}, t)$ for the traction and the divergence theorem, we have

$$ \boldsymbol{f}(t) = \int _{P(t)} \big(\rho (\boldsymbol{x}, t) \boldsymbol{b}(\boldsymbol{x}, t) + \operatorname{div} \boldsymbol{T}( \boldsymbol{x}, t)\big)\, dV. $$

(21)

Substituting (19) and (21) into the balance of linear momentum $\boldsymbol{l}'(t) = \boldsymbol{f}(t)$, and invoking the localization theorem, we find that

$$ (\rho \boldsymbol{v})_{t} +\operatorname{div} (\rho \boldsymbol{v}\otimes \boldsymbol{v}) = \rho \boldsymbol{b}+ \operatorname{div} \boldsymbol{T}. $$

(22)

The mass density $\rho$ for a growing/deforming material must satisfy the balance of mass (16), which implies that

$$ (\rho \boldsymbol{v})_{t} +\operatorname{div} (\rho \boldsymbol{v}\otimes \boldsymbol{v})= \rho \boldsymbol{v}_{t} + \rho (\text{grad}\, \boldsymbol{v})\boldsymbol{v}+ \mu g_{\nu} \boldsymbol{v}. $$

(23)

Substituting (23) into (22) yields the balance of linear momentum for a growing/deforming material body:

$$ \rho (\boldsymbol{x}, t) \big(\boldsymbol{v}_{t}(\boldsymbol{x}, t) + (\text{grad} \, \boldsymbol{v}(\boldsymbol{x}, t))\boldsymbol{v}( \boldsymbol{x},t) \big) + \mu g_{\nu}(\boldsymbol{x}, t) \boldsymbol{v}( \boldsymbol{x}, t) = \rho (\boldsymbol{x}, t) \boldsymbol{b}( \boldsymbol{x}, t) + \operatorname{div} \boldsymbol{T}(\boldsymbol{x}, t). $$

(24)

In the absence of growth, i.e., $g_{\nu} = 0$, equation (24) reduces to the usual equation of motion in continuum mechanics. It is noted that the first two terms in (24) correspond to, respectively, the convective and local accelerations in the spatial description.

While the balance of linear momentum for a growing material contains a new term $\mu g_{\nu} \boldsymbol{v}$ in (24), a detailed analysis for the balance of angular momentum shows that no additional term emerges for a growing material. The equation for the balance of angular momentum of a growing material again states that the stress tensor $\boldsymbol{T}$ is symmetric. The reason for this contrast lies in the fact that the added material particles due to growth bring a change of the linear momentum, with $\mu g_{\nu} \boldsymbol{v}$ being precisely the rate of change, but bring no change of the angular momentum since a material particle is assumed to move with the surrounding material particles without an individual rotational motion.

On the other hand, the presence of growth necessitates fundamental changes to the constitutive equations, which relate the stress to kinematic variables and other state variables. We now turn our attention to this matter.

6 Constitutive Theory for Growing/Deforming Bodies

The current constitutive theory for growing/deforming bodies is based on the fundamental premise that the stress at a point in the current state of a material depends only on the distributions of the material particles in a neighbourhood of the point, and, in particular, on the directional density function.

This is in contrast to the classical elasticity theory, which assumes that the stress at a point in the current configuration depends on the deformation gradient from a reference configuration. We shall show that each deformation gradient corresponds to a directional density function. Therefore, for an elastic material, the stress depends on the directional density function.

The proposed constitutive theory is also in contrast to those in the existing growth theories, which are almost exclusively based on multiplicatively decomposing the deformation gradient tensor into a growth tensor and an elastic tensor. Such an approach to account for the effect of the growth on the stress through deformation is indirect to say the least, as growth and deformation are fundamentally different physical processes in nature. The proposed constitutive theory is based on the understanding that growth and deformation both cause changes of the distribution of the material particles, which, in turn, influence the stress.

6.1 Reformulation of the Elasticity Theory Using the Directional Density Function

In this subsection, we use elasticity theory as a stepping-stone to the development of the proposed constitutive theory.

In classical continuum mechanics, the deformation (motion) of a non-growing material body is described by a function $\bar{\boldsymbol{x}}: \Omega _{0} \times \mathbb{R} \rightarrow \mathbb{R}^{3}$ so that the place $\boldsymbol{x}$ of a material particle $\boldsymbol{X}$ at the time $t$ is given by

$$ \boldsymbol{x}= \bar{\boldsymbol{x}}(\boldsymbol{X}, t). $$

(25)

Here the material particle is labelled by its position vector $\boldsymbol{X}$ in a fixed reference configuration $\Omega _{0} \subset \mathbb{R}^{3}$.

Denoting the set of linear transformations (or tensors) of $\mathbb{R}^{3}$ by Lin, and the set of symmetric tensors by Sym, a homogeneous elastic material is characterized by a response function $\hat{\boldsymbol{T}}$: Lin → Sym, so that the Cauchy stress $\boldsymbol{T}$ at the place $\boldsymbol{x}$ and the time $t$ is given by the constitutive equation

$$ \boldsymbol{T}(\bar{\boldsymbol{x}}(\boldsymbol{X}, t), t) = \hat{\boldsymbol{T}}(\boldsymbol{F}(\boldsymbol{X}, t)), $$

(26)

where $\boldsymbol{F}(\boldsymbol{X}, t) \equiv \partial \bar{\boldsymbol{x}}(\boldsymbol{X}, t)/\partial \boldsymbol{X}$ is the deformation gradient. If the elastic material is isotropic in the reference configuration, the constitutive equation can be written as⁴

$$ \boldsymbol{T}(\bar{\boldsymbol{x}}(\boldsymbol{X}, t), t) = \bar{\boldsymbol{T}}(\boldsymbol{B}(\boldsymbol{X}, t)), $$

(27)

where

$$ \boldsymbol{B}(\boldsymbol{X}, t) \equiv \boldsymbol{F}( \boldsymbol{X},t) \boldsymbol{F}^{T}(\boldsymbol{X}, t) $$

(28)

is the left Cauchy–Green tensor, and the reduced response function $\bar{\boldsymbol{T}}$: Sym → Sym is an isotropic function. A homogeneous isotropic elastic material is thus characterized by a response function $\bar{\boldsymbol{T}}$.

An observation from (28) is that while the deformation gradient $\boldsymbol{F}$ is a linear transformation from the vector space associated with the reference configuration $\Omega _{0}$ to the vector space associated with the current (deformed) configuration $\Omega (t) \equiv \bar{\boldsymbol{x}}(\Omega _{0}, t)$, the left Cauchy–Green tensor $\boldsymbol{B}$ is a linear transformation from the vector space associated with $\Omega (t)$ to itself. It is thus convenient to use $\boldsymbol{B}$ and write the constitutive equation (26) in a spatial form

$$ \boldsymbol{T}(\boldsymbol{x}, t) = \bar{\boldsymbol{T}}( \bar{\boldsymbol{B}}(\boldsymbol{x}, t)), $$

(29)

where $\bar{\boldsymbol{B}}$ is defined by

$$ \bar{\boldsymbol{B}}(\bar{\boldsymbol{x}}(\boldsymbol{X}, t), t) \equiv \boldsymbol{B}(\boldsymbol{X}, t). $$

(30)

We note that although the constitutive equation (29) is presented in the current configuration $\Omega (t)$, a reference configuration $\Omega _{0}$ is still needed in the derivation of (29).

In classical continuum mechanics, an elastic material is treated as a continuum consisting of infinitely many material points. In reality, a material is a discrete system consisting of finite number of material particles, examples of which include atoms, molecules, or biological cells. The continuum hypothesis requires that the number of the material particles in a material body be very large, and that the distribution of the material particles in a neighbourhood of any point be dense. To make a connection between elasticity theory and the present theory, we extend the deformation function of a discreet material body to a smooth function defined on a continuum material body. Here, we consider regular deformations that exclude material penetration, material fragmentation, and other irregular behavior. This ensures that the extended deformation function possesses as much smoothness as desired.

We consider a homogeneous isotropic elastic material that consists of a finite number of material particles uniformly distributed in $\Omega _{0}$. Denote the centroid of the $i$-th material particle by $\boldsymbol{X}_{i}$, which is also regarded as the label for the $i$-th material particle. Define a finite subset $\hat{\Omega}_{0}$ of $\Omega _{0}$ by

$$ \hat{\Omega}_{0} \equiv \{\boldsymbol{X}_{i} \in \Omega _{0}: i = 1, 2, \ldots m\}, $$

(31)

$m$ being the number of material particles. The directional density function of a homogeneous isotropic elastic material in the reference configuration $\Omega _{0}$ is then

$$ n(\boldsymbol{e}, \boldsymbol{X}_{i}) = n_{0}, $$

(32)

where $n_{0}$ is a non-negative constant. The deformation (motion) of an elastic material is described by a discrete deformation function $\hat{\boldsymbol{x}}: \hat{\Omega}_{0} \times \mathbb{R} \rightarrow \mathbb{R}^{3}$. The place $\boldsymbol{x}$ of the material particle $\boldsymbol{X}_{i}$ at the time $t$ is given by

$$ \boldsymbol{x}= \hat{\boldsymbol{x}}(\boldsymbol{X}_{i}, t). $$

(33)

We now extend the function $\hat{\boldsymbol{x}}$ smoothly from $\hat{\Omega}_{0}$ to $\Omega _{0}$ and denote the extended function by $\bar{\boldsymbol{x}}$. As such, $\bar{\boldsymbol{x}}(\boldsymbol{X}_{i}, t) = \hat{\boldsymbol{x}}( \boldsymbol{X}_{i}, t)$, for each $\boldsymbol{X}_{i} \in \hat{\Omega}_{0}$, and $\bar{\boldsymbol{x}}$ is of C¹.

Let $\boldsymbol{\iota }$ be a differential line element in $\Omega _{0}$. At the time $t$, the line element $\boldsymbol{\iota }$ is mapped by $\bar{\boldsymbol{x}}$ to a line element $\boldsymbol{F}\boldsymbol{\iota }$ in $\Omega (t)$, where $\boldsymbol{F}$ is the gradient of the extended deformation function $\bar{\boldsymbol{x}}$. Since the elastic material considered here is non-growing, the number of the material particles associated with the line element $\boldsymbol{\iota }$ in $\Omega _{0}$ is the same as the number of the material particles associated with the line element $\boldsymbol{F}\boldsymbol{\iota }$ in $\Omega (t)$. That is,

$$ |\boldsymbol{F}\boldsymbol{\iota }|\, n( \frac{\boldsymbol{F}\boldsymbol{\iota }}{|\boldsymbol{F}\boldsymbol{\iota }|}, \boldsymbol{x}, t) = |\boldsymbol{\iota }|\, n_{0}. $$

(34)

Letting

$$ \boldsymbol{e}= \frac{\boldsymbol{F}\boldsymbol{\iota }}{|\boldsymbol{F}\boldsymbol{\iota }|} $$

(35)

in (34) leads to an expression for the directional density function in $\Omega (t)$:

$$ n(\boldsymbol{e}, \boldsymbol{x}, t)= \frac{n_{0} |\boldsymbol{\iota }|}{|\boldsymbol{F}\boldsymbol{\iota }|} = n_{0} \, |\boldsymbol{F}^{-1} \boldsymbol{e}| = n_{0} \sqrt{ \boldsymbol{e}\cdot \bar{\boldsymbol{B}}^{-1}(\boldsymbol{x}, t) \boldsymbol{e}}. $$

(36)

We have used (28) and (30) in the last step of (36). It is worth pointing to an analogy between (2) and (36). It indicates that a homogeneous deformation of a material can be regarded as the deformation of a primitive cubic crystal system to the orthorhombic crystal system aforementioned in (2). The cubic crystal system possesses mutually orthogonal lattice vectors $\boldsymbol{b}_{1}$, $\boldsymbol{b}_{2}$, and $\boldsymbol{b}_{3}$ of equal magnitude $|\boldsymbol{b}_{1}|=|\boldsymbol{b}_{2}|=|\boldsymbol{b}_{3}|= n_{0}$. The corresponding deformation gradient is given by

$$ \boldsymbol{F}=\sum _{i=1}^{3} \boldsymbol{a}_{i} \otimes \boldsymbol{b}_{i}. $$

(37)

Equation (36) establishes a relation between the directional density function $n$ in the deformed configuration $\Omega (t)$ and the left Cauchy–Green tensor $\bar{\boldsymbol{B}}$, granted that the material particles before deformation are uniformly distributed and that the discrete deformation function $\hat{\boldsymbol{x}}$ has been smoothly extended to $\Omega _{0}$. An important observation from (36) is that the directional density function $n$ contains all the information contained in the left Cauchy–Green tensor. As a result, we can use $n$, instead of $\bar{\boldsymbol{B}}$, as the primary variable of the response function of a homogeneous isotropic elastic material. We can formally write

$$ \boldsymbol{T}(\boldsymbol{x}, t) = \tilde{\boldsymbol{T}}(n( \boldsymbol{e}, \boldsymbol{x}, t)). $$

(38)

Equation (38) states that the stress depends on the directional density function, and the elastic material is characterized by a new response function $\tilde{\boldsymbol{T}}$, which is related to the response function $\bar{\boldsymbol{T}}$ by

$$ \tilde{\boldsymbol{T}}(n_{0} \sqrt{\boldsymbol{e}\cdot \bar{\boldsymbol{B}}^{-1} \boldsymbol{e}}) = \bar{\boldsymbol{T}}( \bar{\boldsymbol{B}}). $$

(39)

6.2 Toward a General Constitutive Theory Using the Directional Density Function

The preceding subsection serves as a motivation for developing a constitutive theory that uses the directional density function as the primary constitutive variable, a theory which we now describe.

To distinguish the special role of the direction vector $\boldsymbol{e}$ from that of $\boldsymbol{x}$ and $t$ in the directional density function, we define the set ℱ of continuous functions from unit vectors to non-negative real numbers,

$$ {\mathcal{F}}\equiv C^{0}(\overline{\mathbb{R}}^{3}; \mathbb{R}_{+}), $$

(40)

and treat the directional density function as a functional from the space-time, $\mathbb{R}^{3} \times \mathbb{R}$, to ℱ. Specifically, a directional density function induces a function $F \in {\mathcal{F}}$, such that

$$ F(\boldsymbol{x}, t) = n(\boldsymbol{e}, \boldsymbol{x}, t). $$

(41)

With this specification, we propose a constitutive equation

$$ \boldsymbol{T}(\boldsymbol{x}, t) = \tilde{\boldsymbol{T}}(F( \boldsymbol{x}, t)), $$

(42)

where the response function $\tilde{\boldsymbol{T}}$ is regarded as a functional from ℱ to Sym.

The constitutive equation (42) asserts that the stress at the place $\boldsymbol{x}$ and the time $t$ is completely determined by the directional density function at $\boldsymbol{x}$ and $t$. This assertion is supported by various theories in mechanics, such as molecular dynamics, where the stress defined in continuum mechanics is related to the intermolecular forces which depend on the positions of the material particles in the current state.

In spite of the apparent resemblance between (38) and (41)–(42), the constitutive theory proposed here differs significantly from the elasticity theory that led to (38). In fact, it is fundamentally different from most constitutive theories in classical continuum mechanics. Indeed, one of the basic principles in continuum mechanics is the principle of determinism,⁵ which asserts that the stress in a body is determined by the history of the motion of that body. In particular, the stress in an elastic material is determined by the deformation gradient, as indicated in (26). In contrast, the proposed constitutive equation (41)–(42), with the directional density function as the primary constitutive variable, makes no reference of a motion or deformation of the material at all. One of the consequences of this important character of the proposed constitutive theory is that a reference configuration is not required.

While the proposed constitutive theory can describe the behavior of an elastic material, it can describe the behavior of other types of materials studied in continuum mechanics. An example is provided by higher gradient materials, for which a strain-energy function depends on the deformation gradient and higher order deformation gradients. The derivatives of the strain-energy function give stress and higher order stresses (the second order stress is often called couple stress). Such a material can again be described by the constitutive equation (41)–(42), with (36) being replaced by an equation that relates the directional density function to the higher order deformation gradients. In this work, we shall not pursue the classification of the materials that can be described by the constitutive equation (41)–(42).

In (38), the directional density function is related to the left Cauchy–Green tensor through (36) for a homogeneous isotropic elastic material. No such relations are required for the general constitutive equation (41)–(42). However, if the response function of a material under consideration is known, we can use the known response function to find the response function $\tilde{\boldsymbol{T}}$ in (42). We have seen an example of this derivation for a homogeneous isotropic elastic material. If the response function $\bar{\boldsymbol{T}}$ in (29) is known, we can use (39) to find the response function $\tilde{\boldsymbol{T}}$ upon which the constitutive theory developed here is based. Another example will be provided in Sect. 7 for an orthotropic linear elastic material.

6.3 Change of the Directional Density Function in Simultaneous Growth and Deformation Processes

Suppose that at the current time $t$, a growing/deforming material occupies the region $\Omega (t)$. We take $\Omega (t)$ as a reference configuration in which the velocity vector $\boldsymbol{v}(\boldsymbol{x}, t)$ and the directional density function $\boldsymbol{n}(\boldsymbol{e}, \boldsymbol{x}, t)$ are known. We wish to find the directional density function when the material is undergoing a simultaneous growth and deformation process from the current state.

To this end, we consider a deformation function $\boldsymbol{y}: \Omega (t) \times \mathbb{R} \rightarrow \mathbb{R}^{3}$ so that at the time $\tau \geq t$, the place $\boldsymbol{y}$ of a material particle, which occupies $\boldsymbol{x}$ at the time $t$, is given by

$$ \boldsymbol{y}=\boldsymbol{y}(\boldsymbol{x}, \tau ). $$

(43)

Again, for a discrete material system, $\boldsymbol{y}$ is the smooth extension of a discrete deformation function. Obviously, the deformation function must satisfy

$$ \boldsymbol{y}(\boldsymbol{x}, t) = \boldsymbol{x}. $$

(44)

Let $\boldsymbol{\iota }$ be a differential material line element at the material particle located at $\boldsymbol{x}$ at the time $t$. At the time $\tau$, the material particle moves to $\boldsymbol{y}(\boldsymbol{x}, \tau )$, and the material line element $\boldsymbol{\iota }$ is mapped to $\boldsymbol{\iota }(\tau )$, given by

$$ \boldsymbol{\iota }(\tau ) = \boldsymbol{y}_{\boldsymbol{x}}( \boldsymbol{x}, \tau ) \boldsymbol{\iota }, $$

(45)

where $\boldsymbol{y}_{\boldsymbol{x}}$ is the gradient of the smoothly extended deformation function $\boldsymbol{y}$. Let $d \tau$ be an infinitesimal time increment. At the time $\tau + d \tau$, the material particle moves to the place $\boldsymbol{y}(\boldsymbol{x}, \tau + d \tau )$, and the material line element $\boldsymbol{\iota }$ is mapped to

$$ \boldsymbol{\iota }(\tau + d \tau ) = \boldsymbol{y}_{\boldsymbol{x}}( \boldsymbol{x}, \tau + d \tau ) \boldsymbol{\iota }. $$

(46)

Denote by $N(\tau )$ the number of the material particles on the line element $\boldsymbol{\iota }(\tau )$ at the time $\tau$. Then,

$$ N(\tau ) = |\boldsymbol{\iota }(\tau )| n( \frac{\boldsymbol{\iota }(\tau )}{|\boldsymbol{\iota }(\tau )|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ). $$

(47)

Due to the growth, the number of the material particles on the line element $\boldsymbol{\iota }(\tau )$ at the time $\tau + d \tau$ is

$$ N(\tau + d \tau ) = |\boldsymbol{\iota }(\tau )| \Big(n( \frac{\boldsymbol{\iota }(\tau )}{|\boldsymbol{\iota }(\tau )|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) + g( \frac{\boldsymbol{\iota }(\tau )}{|\boldsymbol{\iota }(\tau )|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) d \tau \Big). $$

(48)

Here, the definition of the growth rate function $g$, as stated in (3), has been invoked with the understanding that $d \tau$ is infinitesimally small. On the other hand, $N(\tau + d \tau )$ is also the number of the material particles on the line element $\boldsymbol{\iota }(\tau + d \tau )$, which is given by

$$\begin{aligned} N(\tau + d \tau ) &= |\boldsymbol{\iota }(\tau + d\tau )|\; n( \frac{\boldsymbol{\iota }(\tau +d \tau )}{|\boldsymbol{\iota }(\tau + d \tau )|}, \boldsymbol{y}(\boldsymbol{x}, \tau +d \tau ), \tau +d\tau ) \\ &= |\boldsymbol{\iota }(\tau )|\; n( \frac{\boldsymbol{\iota }(\tau )}{|\boldsymbol{\iota }(\tau )|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau )+ \frac{\partial}{\partial \tau}\Big[ |\boldsymbol{\iota }(\tau )|\; n( \frac{\boldsymbol{\iota }(\tau )}{|\boldsymbol{\iota }(\tau )|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau )\Big] d \tau . \end{aligned}$$

(49)

It then follows from (48) and (49) that

$$ |\boldsymbol{\iota }(\tau )|\; g( \frac{\boldsymbol{\iota }(\tau )}{|\boldsymbol{\iota }(\tau )|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) = \frac{\partial}{\partial \tau}\Big[ |\boldsymbol{\iota }(\tau )|\; n( \frac{\boldsymbol{\iota }(\tau )}{|\boldsymbol{\iota }(\tau )|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau )\Big], $$

(50)

which can be written, by (45), as

$$ |\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{\iota }| \; g( \frac{\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{\iota }}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{\iota }|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) = \frac{\partial}{\partial \tau}\Big( |\boldsymbol{y}_{\boldsymbol{x}}( \boldsymbol{x}, \tau ) \boldsymbol{\iota }|\; n( \frac{\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{\iota }}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{\iota }|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau )\Big). $$

(51)

Integrating (51) over the time interval $[t, \tau ]$ and using (44), we find that

$$\begin{aligned} \int _{t}^{\tau}& |\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{\iota }| \, g( \frac{\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{\iota }}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{\iota }|}, \boldsymbol{y}(\boldsymbol{x}, s), s) ds \\ & = |\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{\iota }|\, n( \frac{\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{\iota }}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{\iota }|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) - |\boldsymbol{\iota }| \, n(\frac{\boldsymbol{\iota }}{|\boldsymbol{\iota }|}, \boldsymbol{x}, t). \end{aligned}$$

(52)

Letting $\boldsymbol{e}\equiv \boldsymbol{\iota }/|\boldsymbol{\iota }|$, we write (52) as

$$\begin{aligned} n(\boldsymbol{e}, \boldsymbol{x}, t)+\int _{t}^{\tau}& | \boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{e}| \, g( \frac{\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{e}}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{e}|}, \boldsymbol{y}(\boldsymbol{x}, s), s) ds \\ & = |\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}|\, n( \frac{\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}|}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) . \end{aligned}$$

(53)

Equation (53) dictates how the directional density function $n$ evolves in a growing/deforming process, and is termed the evolution equation for the directional density function. We note that in the mechanics literature, an evolution equation for a state variable is often written in a rate form, which in the current theory is essentially (51). This rate-form evolution equation is a first-order partial differential equation for the direction density function $n$, whose solution is simply (53). We also note that if a material grows without deformation from the current configuration, we will have $\boldsymbol{y}(\boldsymbol{x}, \tau )=\boldsymbol{x}$ and recover (4) from (53).

7 Growth of a Tree Trunk, an Example of Internal Surface Growth

One of the primary growth mechanisms of a tree trunk is the nucleation of cells near the surface. Between the wood and bark is a thin layer of dividing meristematic cells called the vascular cambium, which produces new cells that move inward toward the wood and outward toward the tree bark. This represents a physical process that can be idealized as an internal surface growth. In this section, we use the theory developed in the preceding sections to analyze tree truck growth.

7.1 Growth Rate Function and the Corresponding Motion of the Material Particles

The primary input quantity (prescribed quantity) is the growth rate function $g(\boldsymbol{e}, \boldsymbol{x}, t)$. For a given $g$, we wish to compute other quantities of interest, such as the motion and the stress of the material particles in the tree.

To this end, we prescribe a growth rate function by a distribution function that is zero except on a moving circular cylindrical surface. Let $\boldsymbol{e}_{r}$, $\boldsymbol{e}_{\theta}$, and $\boldsymbol{e}_{z}$ be the unit base vectors in a cylindrical coordinate system, with the $z$-coordinate axis being along the centerline of the tree trunk, which is modeled as an evolving circular cylinder. For the present problem, it suffices to specify the growth rate function in the directions of the base vectors. We assume that

$$ g(\boldsymbol{e}_{r}, \boldsymbol{x}, t) = g_{r}(t) \delta (r( \boldsymbol{x})-R(t)), \; g(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t) = g_{\theta}(t) \delta (r(\boldsymbol{x})-R(t)), \; g(\boldsymbol{e}_{z}, \boldsymbol{x}, t)=0, $$

(54)

where $g_{r}$, $g_{\theta}$, and $R$ are prescribed functions from ℝ to ℝ, $\delta$ is the Dirac $\delta$-function, and

$$ r(\boldsymbol{x}) \equiv \boldsymbol{x}\cdot \boldsymbol{e}_{r}. $$

(55)

We assume that $R$ is a continuous and strictly increasing function of $t$, which implies that the circular cylindrical surface defined by $r(\boldsymbol{x})=R(t)$ is moving outward. The growth rate function given in (54) is a mathematical idealization of an anisotropic surface growth which occurs only on the moving surface $r(\boldsymbol{x})=R(t)$. We note that the surface $r(\boldsymbol{x})=R(t)$ is not a material surface, meaning that it is not associated with any set of material particles. It is just a moving locus where new material particles are deposited. Also, for simplicity, we choose not to consider, by (54)₃, the axial growth of the tree trunk, which can nevertheless be incorporated if needed.

We look for a solution for which the velocity is of the form

$$ \boldsymbol{v}(\boldsymbol{x}, t) = R'(t) H(r(\boldsymbol{x})- R(t)) \boldsymbol{e}_{r}, $$

(56)

where $H$ is the Heaviside function:

$$ H(x) \equiv \left \{ \textstyle\begin{array}{l@{\quad}l} 0, \vspace{5mm} \hspace{5mm} & {\mathrm{if}} \, x < 0, \vspace{-4mm} \\ 1, & {\mathrm{if}}\, x \geq 0. \end{array}\displaystyle \right . $$

(57)

The velocity field (56) asserts that the material particles behind the moving surface $r(\boldsymbol{x})=R(t)$ remain stationary, and that the material particles ahead of the moving surface move with the same velocity as the surface.

The deformation function $\boldsymbol{y}(\boldsymbol{x}, \tau )$ associated with the velocity field (56) can be found by solving the velocity-deformation equation (1). For the current spatial description, this amounts to solving the initial value problem

$$\begin{aligned} \frac{\partial \boldsymbol{y}(\boldsymbol{x}, \tau )}{\partial \tau} &= \boldsymbol{v}(\boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) = R'( \tau ) H(r(\boldsymbol{y}(\boldsymbol{x}, \tau )) - R(\tau )) \boldsymbol{e}_{r}, \end{aligned}$$

(58)

$$\begin{aligned} \boldsymbol{y}(\boldsymbol{x}, t)&=\boldsymbol{x}. \end{aligned}$$

(59)

The solution to this initial value problem is

$$ \boldsymbol{y}(\boldsymbol{x}, \tau )=\boldsymbol{x}+ \big(R(\tau )-R(t) \big) \,\big(1- H(R(t)-r(\boldsymbol{x}))\big) \boldsymbol{e}_{r}. $$

(60)

It follows from (60) that

$$\begin{aligned} \boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) = & \boldsymbol{I}+ \big(R(\tau )-R(t)\big) \delta (R(t)-r(\boldsymbol{x})) \boldsymbol{e}_{r} \otimes \boldsymbol{e}_{r} \\ &+ \frac{R(\tau )-R(t)}{r(\boldsymbol{x})}\big(1-H(R(t)-r( \boldsymbol{x}))\big) \boldsymbol{e}_{\theta} \otimes \boldsymbol{e}_{ \theta}. \end{aligned}$$

(61)

7.2 Evolution of the Directional Density Function

We can now find the directional density function in the growth process. By (61), we find that

$$ \frac{\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}|} = \boldsymbol{e}\quad \mathrm{{for} \quad \boldsymbol{e}= \boldsymbol{e}_{r}},\mathrm{ \boldsymbol{e}_{\theta}},\mathrm{ \boldsymbol{e}_{z}. } \mathrm{ }$$

(62)

Also, it follows from (54)₁, (55), and (60) that

$$\begin{aligned} g(\boldsymbol{e}_{r}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) &= g_{r}(\tau ) \delta (r(\boldsymbol{y}(\boldsymbol{x}, \tau )) - R( \tau )) \\ &=g_{r}(\tau ) \delta (r(\boldsymbol{x})-(R(\tau )-R(t)) H(R(t)-r( \boldsymbol{x})) - R(t)) \\ &=g_{r}(\tau ) H(R(t)-r(\boldsymbol{x})) \delta (r(\boldsymbol{x})-R( \tau )). \end{aligned}$$

(63)

By substituting (61), (62), and (63) into the evolution equation (53), we find that

$$\begin{aligned} n(\boldsymbol{e}_{r}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) &= \frac{n(\boldsymbol{e}_{r}, \boldsymbol{x}, t)}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}_{r}|}+ \int _{t}^{\tau} \frac{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{e}_{r}|}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}_{r}|} \, g(\boldsymbol{e}_{r}, \boldsymbol{y}(\boldsymbol{x}, s), s)\, ds \\ &= \frac{n(\boldsymbol{e}_{r}, \boldsymbol{x}, t)}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}_{r}|}+ \int _{t}^{\tau} \frac{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{e}_{r}|}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}_{r}|} \, g_{r}(s) H(R(t)-r(\boldsymbol{x})) \delta (r(\boldsymbol{x})-R(s))\, ds \\ & = \left \{ \textstyle\begin{array}{l@{\quad}l} n(\boldsymbol{e}_{r}, \boldsymbol{x}, t), \hspace{11mm} & \text{if}\; r(\boldsymbol{x}) \neq R(t), \\ \displaystyle \frac{g_{r}(t)}{R'(t)}, & \text{if}\; r(\boldsymbol{x}) = R(t). \end{array}\displaystyle \right . \end{aligned}$$

(64)

Some observations from (64) are in order. At the time $t$, the moving surface, where growth occurs, occupies the points $\boldsymbol{x}$ for which $r(\boldsymbol{x}) = R(t)$. The region for which $r(\boldsymbol{x})< R(t)$ is occupied by the material particles behind the moving surface in the wood of the tree trunk. These material particles, in accord with (56), are motionless. Therefore, the arrangement of the material particles remains unchanged. As a result, the directional number density remains unchanged in the subsequent growth process. That is, $n(\boldsymbol{e}_{r}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau )=n( \boldsymbol{e}_{r}, \boldsymbol{x}, \boldsymbol{t})$, as stated in (64)₁. Note that (64)₁ also indicates that $n(\boldsymbol{e}_{r}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau )=n( \boldsymbol{e}_{r}, \boldsymbol{x}, \boldsymbol{t})$ if $r(\boldsymbol{x}) > R(t)$. Indeed, the region for which $r(\boldsymbol{x})>R(t)$ is occupied by the material particles ahead of the moving surface in the bark of the tree trunk. These material particles, in accord with (56), move outward at the same velocity $R'(t)$. Therefore, the radial arrangement of the material particles remains unchanged: hence the claimed result.

On the other hand, (64)₂ indicates that $n(\boldsymbol{e}_{r}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau )=g_{r}(t)/R'(t)$ if $r(\boldsymbol{x}) = R(t)$, which is at the moving surface. By (54)₁, $g_{r}(t)/\epsilon$ is the rate of change of the radial directional number density as the moving surface passes through an infinitesimal region of width $\epsilon$. Also, $\epsilon/R'(t)$ is the time required for the moving surface to pass the infinitesimal region. Therefore, $g_{r}(t)/R'(t)$ is the increase of the radial directional number density $n(\boldsymbol{e}_{r}, \boldsymbol{x}, t)$ due to surface growth. A seemingly unintuitive result of (64)₂ is that the value of the radial directional density increases from zero to $g_{r}(t)/R'(t)$ as the moving growth surface passes. This is because the infinitesimal region instantaneously occupied by the moving surface grows into a finite region and the material particles within the infinitesimal region spread out radially over the finite region, leading to zero density. To put it differently, the radial directional number density would be zero as the surface $r(\boldsymbol{x})=R(t)$ moves forward, had $g_{r}$ been zero. Of course, this scenario cannot actually happen, as the surface growth is the factor that drives the moving surface of the velocity discontinuity. If $g_{r}$ were zero, the velocity field (56) would not be realized.

In a similar manner, we can find the change of the circumferential directional density function $n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)$. It follows from (54)₂, (55), and (60) that

$$ g(\boldsymbol{e}_{\theta}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) =g_{\theta}(\tau ) H(R(t)-r(\boldsymbol{x})) \delta (r( \boldsymbol{x})-R(\tau )). $$

(65)

By substituting (61), (62), and (65) into (53), we find that

$$\begin{aligned} n(\boldsymbol{e}_{\theta}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau ) &= \frac{n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}_{\theta}|}+ \int _{t}^{\tau} \frac{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{e}_{\theta}|}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}_{\theta}|} \, g(\boldsymbol{e}_{\theta}, \boldsymbol{y}(\boldsymbol{x}, s), s)\, ds \\ &= \frac{n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}_{\theta}|} \\ &\quad {}+ \int _{t}^{\tau} \frac{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, s) \boldsymbol{e}_{\theta}|}{|\boldsymbol{y}_{\boldsymbol{x}}(\boldsymbol{x}, \tau ) \boldsymbol{e}_{\theta}|} \, g_{\theta}(s) H(R(t)-r(\boldsymbol{x})) \delta (r(\boldsymbol{x})-R(s))\, ds \\ & = \left \{ \textstyle\begin{array}{l@{\quad}l} n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t), \hspace{19mm} & \text{if}\; r(\boldsymbol{x}) < R(t), \\ \displaystyle n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t) + \frac{g_{\theta}(t)}{R'(t)}, & \text{if}\; r(\boldsymbol{x}) = R(t), \\ \displaystyle \frac{R(t)}{R(\tau )} n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t), & \text{if}\; r(\boldsymbol{x}) > R(t). \end{array}\displaystyle \right . \end{aligned}$$

(66)

We next make some observations from (66). The statement of (66)₁, that is, $n(\boldsymbol{e}_{\theta}, \boldsymbol{y}(\boldsymbol{x}, \tau ), \tau )=n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)$ if $r(\boldsymbol{x}) < R(t)$, can be justified again by the argument that the material particles behind the moving surface are motionless, which implies that the arrangement of the material particles and the directional number density remain unchanged. However, unlike the radial directional number density, the circumferential directional number density changes ahead of the moving surface, as indicated by (66)₃. Indeed, the material particles ahead of the moving surface move outward at the same velocity. While this motion does not change the distribution of the material particles in the radial direction, it does change the distribution of the material particles in the circumferential direction. Precisely, as the radius of the moving surface changes from $R(t)$ to $R(\tau )$, the circumferential directional density changes from $n(\boldsymbol{e}_{\theta},\boldsymbol{x}, t)$ to $n(\boldsymbol{e}_{\theta},\boldsymbol{x}, t)\, R(t)/R(\tau )$, as stated in (66)₃. Moreover, (66)₂ can be justified by a similar argument previously employed for justifying (64)₂. The term $g_{\theta}(t)/R'(t)$ is again the increase of the circumferential directional density $n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)$ due to the surface growth. However, unlike the situation with the radial directional density function, the value of the circumferential directional density function increases from $n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)$ to $n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t) + g_{\theta}(t)/R'(t)$. This is because while the materials particles within the infinitesimal region instantaneously occupied by the moving surface spread out in the radial direction over a finite region as the surface moves, they do not spread in the circumferential direction.

The independent variables in (64) and (66) are $\boldsymbol{x}$ and $\tau$. By using (60), we can rewrite (64) and (66) in a form with $\boldsymbol{y}$ and $\tau$ being the independent variables:

$$\begin{aligned} n(\boldsymbol{e}_{r}, \boldsymbol{y}, \tau ) &= \left \{ \textstyle\begin{array}{l@{\quad}l} n(\boldsymbol{e}_{r}, \boldsymbol{y}, t), \hspace{38mm} & \text{if}\; r(\boldsymbol{y}) < R(t), \\ \displaystyle \frac{g_{r}(\hat{t}(r(\boldsymbol{y})))}{R'(\hat{t}(r(\boldsymbol{y})))}, & \text{if}\; R(t) < r(\boldsymbol{y}) < R(\tau ), \\ n(\boldsymbol{e}_{r}, \boldsymbol{y}-(R(\tau ) - R(t)) \boldsymbol{e}_{r}, t), & \text{if}\; r(\boldsymbol{y}) > R(\tau ). \end{array}\displaystyle \right . \end{aligned}$$

(67)

$$\begin{aligned} n(\boldsymbol{e}_{\theta}, \boldsymbol{y}, \tau ) &= \left \{ \textstyle\begin{array}{l@{\quad}l} n(\boldsymbol{e}_{\theta}, \boldsymbol{y}, t), \hspace{38mm} & \text{if}\; r(\boldsymbol{y}) < R(t), \\ \displaystyle \frac{R(t)}{r(\boldsymbol{y})} n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t) + \frac{g_{\theta}(\hat{t}(r(\boldsymbol{y})))}{R'(\hat{t}(r(\boldsymbol{y})))}, & \text{if}\; R(t) < r(\boldsymbol{y}) < R(\tau ), \\ \displaystyle \frac{R(t)}{R(\tau )} n(\boldsymbol{e}_{\theta}, \boldsymbol{y}-(R(\tau ) - R(t)) \boldsymbol{e}_{r}, t), & \text{if}\; r( \boldsymbol{y}) > R(\tau ). \end{array}\displaystyle \right . \end{aligned}$$

(68)

Here $\hat{t}$ is the inverse function of $R$, that is, $\hat{t}(R(t)) = t$. Since $R$ is a strictly increasing function of $t$, its inverse is well-defined. We note that the region $r(\boldsymbol{y}) < R(t)$ at time $\tau$ is occupied by the same material particles in the region $r(\boldsymbol{x})< R(t)$ at time $t$. These material particles are stationary in the time interval $[t, \tau ]$ and are not affected by growth, and thus the directional density function remains unchanged. On the other hand, the region $r(\boldsymbol{y})>R(\tau )$ at the time $\tau$ is occupied by the same material particles in the region $r(\boldsymbol{x}) > R(t)$ at time $t$, that is, the material particles in the tree bark. These material particles move outward in the time interval $[t, \tau ]$. The change of the directional density function is due to the change of places of the material particles. Most interesting is the region $R(t) < r(\boldsymbol{y}) < R(\tau )$, which is occupied by the new material particles from the growth in the time interval $[t, \tau ]$. This three-dimensional region at time $\tau$ corresponds to the two-dimensional surface $r(\boldsymbol{x}) = R(t)$ at time $t$. Such a change of topology is a special feature of a surface growth in which a surface grows into a volume. As expected, the directional density function in this region depends on the growth rate function $g$.

7.3 Determination of the Stress Field

The development in the preceding two subsections, concerning the motion of the material particles and the directional density function, makes no reference to a constitutive theory and the change of the stress field due to growth. In the remainder of this work, we undertake the task of finding the components of the stress in the tree trunk during growth. This is motivated by the observation that wood may wrap after being cut from the tree trunk, which is due to the existence of residual stress in the tree trunk,⁶ in addition to other conditions such as uneven moistures.

To obtain definite results, we assume that the tree trunk is modeled as an orthotropic linear elastic material, for which the material symmetry group contains reflections on three mutually perpendicular planes defined by $\boldsymbol{e}_{r}$, $\boldsymbol{e}_{\theta}$, and $\boldsymbol{e}_{z}$ for the present problem. The constitutive equations for this material are given⁷ by

$$ \left [ \textstyle\begin{array}{c} T_{rr} \\ T_{\theta \theta} \\ T_{zz} \end{array}\displaystyle \right ] = \left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c} C_{rr} & C_{r\theta} & C_{zr} \\ C_{r\theta} & C_{\theta \theta} & C_{\theta z} \\ C_{z r} & C_{\theta z} & C_{z z} \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c} E_{rr} \\ E_{\theta \theta} \\ E_{zz} \end{array}\displaystyle \right ] \quad \mathrm{{and} \quad }\left [\mathrm{ \textstyle\begin{array}{c} T_{r\theta} \\ T_{\theta z} \\ T_{z r} \end{array}\displaystyle }\right ]\mathrm{ = }\left [\mathrm{ \textstyle\begin{array}{c@{\quad}c@{\quad}c} D_{r \theta} & 0 & 0 \\ 0 & D_{\theta z} & 0 \\ 0 & 0 & D_{z r} \end{array}\displaystyle }\right ] \left [\mathrm{ \textstyle\begin{array}{c} E_{r \theta} \\ E_{\theta z} \\ E_{z r} \end{array}\displaystyle }\right ], $$

(69)

where $T_{rr}$, etc. are the components of the stress in cylindrical coordinates, $E_{rr}$, etc. are the components of the infinitesimal strain tensor $\boldsymbol{E} \equiv \frac{1}{2}(\boldsymbol{F}+\boldsymbol{F}^{T})-\boldsymbol{I}$, and $C_{rr}$, etc. and $D_{r\theta}$, etc are nine material constants. For the present problem, we have $E_{zz} = E_{r \theta} = E_{\theta z} = E_{z r} =0$. Hence, the relevant part of the constitutive equations involves two components of the stress, two strain components, and three material constants:

$$ \left [ \textstyle\begin{array}{c} T_{rr} \\ T_{\theta \theta} \end{array}\displaystyle \right ] = \left [ \textstyle\begin{array}{c@{\quad}c} C_{rr} & C_{r\theta} \\ C_{r\theta} & C_{\theta \theta} \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c} E_{rr} \\ E_{\theta \theta} \end{array}\displaystyle \right ]. $$

(70)

We note that $(C_{rr} C_{\theta \theta}-C^{2}_{r \theta})/C_{\theta \theta}$ and $C_{r \theta}/C_{\theta \theta}$ correspond to the Young’s modulus and the Poisson’s ratio for uniaxial tension in the $\boldsymbol{e}_{r}$ direction, and $(C_{rr} C_{\theta \theta}-C^{2}_{r \theta})/C_{r r}$ and $C_{r \theta}/C_{r r}$ to the Young’s modulus and the Poisson’s ratio in the $\boldsymbol{e}_{\theta}$ direction.

Using the same argument that has led to (36), we find that the components of the strain are related to the directional density function as

$$ E_{r r}(\boldsymbol{x}, t) = \frac{n_{0}}{n(\boldsymbol{e}_{r}, \boldsymbol{x}, t)}-1 \quad \text{and} \quad E_{\theta \theta}(\boldsymbol{x}, t) = \frac{n_{0}}{n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)}-1. $$

(71)

Upon substituting (71) into (70), we find that

$$ T_{rr}(\boldsymbol{x}, t) = C_{rr}\big( \frac{n_{0}}{n(\boldsymbol{e}_{r}, \boldsymbol{x}, t)}-1\big) + C_{r \theta}\big( \frac{n_{0}}{n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)}-1\big) $$

(72)

and that

$$ T_{\theta \theta}(\boldsymbol{x}, t) = C_{r \theta}\big( \frac{n_{0}}{n(\boldsymbol{e}_{r}, \boldsymbol{x}, t)}-1\big) + C_{ \theta \theta}\big( \frac{n_{0}}{n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)}-1\big). $$

(73)

We note that equations (72) and (73) are the relevant part of the response function $\tilde{\boldsymbol{T}}$ appearing in (38) for the orthotropic linear elastic material.

We can now draw some definite conclusions regarding the relation between the state of stress and the process of growth.

As evident from (64)₂ and (66)₂, or from (67)₂ and (68)₂, the directional density function $n$ increases in the growth rate function $g$. Under the usual assumption that the material constants $C_{rr}$, $C_{\theta \theta}$, and $C_{r \theta}$ are positive, or that the Young’s moduli and the Poisson’s ratios are positive, the radial stress $T_{rr}$ and the circumferential stress $T_{\theta \theta}$, given by (72) and (73) respectively, decrease in the values of both the radial directional density function $n(\boldsymbol{e}_{r}, \boldsymbol{x}, t)$ and the circumferential directional density function $n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)$, and therefore, decrease in the values of both the radial and circumferential growth rate functions $g(\boldsymbol{e}_{r}, \boldsymbol{x}, t)$ and $g(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t)$. While it is known that a radial growth tends to create a compressive radial stress, it may be less recognized that a circumferential growth also has the effect of creating a compressive radial stress. A similar statement can be made of the compressive circumferential stress.

For a prescribed growth rate function $g$, the directional density function in the region $R(t) < r(\boldsymbol{y}) < R(\tau )$, which is generated by the surface growth over the time interval $[t, \tau ]$, can be obtained from (67)₂ and (68)₂. Then, the stresses in this region can be obtained from (72) and (73). The control of such residual stresses in wood is of great importance to the lumber industry. As an ideal situation, suppose that the wood is stress free at time $t$, that is, $n(\boldsymbol{e}_{r}, \boldsymbol{x}, t) = n(\boldsymbol{e}_{\theta}, \boldsymbol{x}, t) = n_{0}$. It then follows from (67)₂ and (68)₂ that if we can control the growth rate function in such a way that

$$ g_{r}(s) = n_{0} R'(s) \quad \text{and} \quad g_{\theta}(s) = n_{0} R'(s) \big(1-\frac{R(t)}{R(s)}\big) $$

(74)

for all $s \in [t, \tau ]$, the wood in the grown trunk will also be stress-free.

The above discussion has been focused on the region behind the moving growth surface, where the directional density function, and therefore the stresses, depend heavily on the growth rate function $g$. In the region ahead of the moving growth surface $r(\boldsymbol{y}) > R(\tau )$, where the tree bark lies, the value of the radial directional density function remains constant, as stated in (67)₃, and the value of the circumferential directional density function decreases in a factor of $R(t)/R(\tau )$, as stated in (68)₃. This decrease of the circumferential density results from the increase of diameter of the tree trunk during the growth process, and in turn, results in an increase of the circumferential stress $T_{\theta \theta}$, as can be seen from (73). This inevitably leads to a tensile circumferential stress in the tree bark, which may eventually break the tree bark. The patterns of the tree bark in various tree species are likely dependent on the circumferential stress, the axial stress (which is due to the axial growth), and the fracture properties of the bark material.

As a concluding remark, we note that in principle the intensity of the radial growth rate $g_{r}$ and the intensity of the circumferential growth rate $g_{\theta}$ can be arbitrarily prescribed. However, they must satisfy certain conditions for the velocity field (56) to hold. Such conditions can be derived from the balance of linear momentum (24). In the region behind the moving growth surface, the velocity is zero and the balance of linear momentum (24) reduces to the equilibrium equation

$$ \mathrm{div} \,\boldsymbol{T}= \mathbf{0. } \mathbf{ }$$

(75)

For the current problem, the equilibrium equation takes the following component form:

$$ \frac{d T_{rr}}{d r} + \frac{T_{rr} - T_{\theta \theta}}{r}=0. $$

(76)

In the region $R(t) < r(\boldsymbol{y}) < R(\tau )$ which is created by the surface growth over the time interval $(t, \tau )$, the radial and circumferential directional density functions are given, respectively, by (64)₂ and (66)₂. Upon substituting (64)₂, (66)₂, (72), and (73) into (76), we find that

$$\begin{aligned} \Big[& \frac{1}{R'(s)} \frac{d}{d s}\big( \frac{C_{rr} n_{0} R(s) R'(s)}{g_{r}(s)}+ \frac{C_{r\theta} n_{0} R^{2}(s) R'(s)}{R(t) R'(s)+R(s) g_{\theta}(s)} \Big)-\frac{C_{r \theta} n_{0} R'(s)}{g_{r}(s)} \\ &- \frac{C_{\theta \theta} n_{0} R(s) R'(s)}{R(t) R'(s)+R(s) g_{\theta}(s)} \Big]_{s=\hat{t}(r(\boldsymbol{y}))} =C_{rr}-C_{\theta \theta}. \end{aligned}$$

(77)

Equation (77) is the condition that the radial growth rate $g_{r}$ and the circumferential growth rate $g_{\theta}$ must satisfy in order to realize the velocity field (56). In other words, if (77) is not satisfied, the material particles behind the moving growth surface will be set in motion. Actually, this may well be so in reality. However, such a motion is literally negligible since the velocity $R'(t)$ of the moving growth surface is so slow (in the order of millimeter per year for most trees), and the velocity field (56) is practically valid.