2.10 Constrained Reactive Mixture of Solids

A solid material may consist of a heterogeneous mixture of various solid constituents that are constrained to move together. In this section we describe fundamental considerations for such constrained mixtures.

Mixture Kinematics

Consider a constrained mixture of multiple solid constituents \(\sigma\). The motion of each constituent is given by \(\boldsymbol{\chi}^{\sigma}\left(\mathbf{X}^{\sigma},t\right)\), where \(\mathbf{X}^{\sigma}\) denotes a material point in the reference configuration of that constituent. At the current time \(t\), various constituents \(\sigma\) which occupy an elemental region with a spatial position \(\mathbf{x}\) may have originated from distinct referential positions \(\mathbf{X}^{\sigma}\). We often label a convenient constituent as the master constituent \(s\) (e.g., the oldest constituent in a reactive mixture with evolving composition) and call the reference configuration \(\mathbf{X}^{s}\) the master reference configuration. All the referential mass densities and mass density supplies (see below) are evaluated relative to the master reference configuration \(\mathbf{X}^{s}\). The kinematics of each constituent \(\sigma\) may be related to the kinematics of the master constituent \(s\) through

\[ \begin{equation} \mathbf{x}=\boldsymbol{\chi}^{s}\left(\mathbf{X}^{s},t\right)=\boldsymbol{\chi}^{\sigma}\left(\mathbf{X}^{\sigma},t\right)\,.\label{eq:general-mixture-kinematics} \end{equation} \]

Taking the material time derivative of this relation in the material frame and recognizing that this relation must hold for all \(t\) in the case of a constrained mixture establishes that all constituents share the same velocity \(\mathbf{v}^{\sigma}=\mathbf{v}^{s}\). However, as detailed previously ¹², constituents may have distinct deformation gradients \(\mathbf{F}^{\sigma}=\partial\boldsymbol{\chi}^{\sigma}/\partial\mathbf{X}^{\sigma}\). When a reaction converts a reactant \(\sigma=a\) into a product \(\sigma=b\), these constituents may have distinct reference configurations. The deformation gradient of the master constituent \(\mathbf{F}^{s}\), which also serves as the total deformation gradient, may be related to the relative deformation gradient \(\mathbf{F}^{\sigma}\) of constituent \(\sigma\) by applying the chain rule to eq.\eqref{eq:general-mixture-kinematics}, producing

\[ \begin{equation} \mathbf{F}^{s}\left(\mathbf{X}^{s},t\right)=\frac{\partial\boldsymbol{\chi}^{s}\left(\mathbf{X}^{s},t\right)}{\partial\mathbf{X}^{\sigma}}\cdot\frac{\partial\mathbf{X}^{\sigma}\left(\mathbf{X}^{s}\right)}{\partial\mathbf{X}^{s}}=\mathbf{F}^{\sigma}\left(\mathbf{X}^{s},t\right)\cdot\mathbf{F}^{\sigma s}\left(\mathbf{X}^{s}\right)\,.\label{eq:general-total-def-grad} \end{equation} \]

In this expression, \(\mathbf{F}^{\sigma s}\left(\mathbf{X}^{s}\right)\) is the deformation gradient of \(\sigma\) relative to \(s\), which must be postulated by constitutive assumption. The relationship between \(\mathbf{X}^{\sigma}\) and \(\mathbf{X}^{s}\) is time-invariant; consequently, \(\mathbf{F}^{\sigma s}\) is a time-invariant spatial mapping. It follows that only one deformation gradient represents an independent state variable in a constrained mixture framework, whereas all others are related to it via eq.\eqref{eq:general-total-def-grad}; any of the \(\mathbf{F}^{\sigma}\)'s may be selected, based on convenience. In eq.\eqref{eq:general-total-def-grad} the spatio-temporal arguments have been written explicitly for clarity. These dependencies are implied in the forthcoming sections and henceforth those arguments may be selectively suppressed. Taking the determinant of eq.\eqref{eq:general-total-def-grad} produces a relation between the volume ratios \(J^{\sigma}=\det\mathbf{F}^{\sigma}\) and \(J^{s}=\det\mathbf{F}^{s}\),

\[ \begin{equation} J^{s}=J^{\sigma}J^{\sigma s}\,,\label{eq:general-jacobian-relation} \end{equation} \]

where \(J^{\sigma s}=\det\mathbf{F}^{\sigma s}\).

Mixture Composition

Each constituent \(\sigma\) has an apparent mass density \(\rho^{\sigma}\) which may evolve due to deformation, or due to reactive processes which alter the mixture composition. Following ³⁴, we define the referential apparent mass density \(\rho_{r}^{\sigma}\) of each constituent as

\[ \begin{equation} \rho_{r}^{\sigma}=J^{s}\rho^{\sigma}\,.\label{eq:general-referential-mass-density} \end{equation} \]

Equation \eqref{eq:general-referential-mass-density} expresses the mass of constituent \(\sigma\) per referential volume of the master constituent \(s\); thus \(\rho_{r}^{\sigma}\) may only evolve if the mass content changes via reactions, making it a suitable state variable for tracking composition in a reactive framework. The axiom of mass balance for each constituent \(\sigma\) may be written as

\[ \begin{equation} \dot{\rho}_{r}^{\sigma}=\hat{\rho}_{r}^{\sigma}\,,\label{eq:mass-balance-alpha-1} \end{equation} \]

where the dot operator represents the material time derivative and \(\hat{\rho}_{r}^{\sigma}\) is the referential mass supply density for constituent \(\sigma\), representing the rate at which mass (per referential volume) is added to \(\sigma\) due to reactions with all other mixture constituents ¹². A constitutive relation must be provided for \(\hat{\rho}_{r}^{\sigma}\) for various types of reactions. The mixture referential mass density \(\rho_{r}\) is given by

\[ \begin{equation} \rho_{r}=\sum_{\sigma}\rho_{r}^{\sigma}\,.\label{eq:general-mass-balance} \end{equation} \]

This summation is carried out over all constituents. When a constrained mixture of solid constituents represents a closed system, \(\rho_{r}\) remains constant over time. Taking the material time derivative of eq.\eqref{eq:general-mass-balance} and using eq.\eqref{eq:mass-balance-alpha-1} shows that the referential mass density supplies must satisfy

\[ \begin{equation} \sum_{\sigma}\hat{\rho}_{r}^{\sigma}=0\,.\label{eq:mass-balance-constraint} \end{equation} \]

However, if a constrained mixture of solid constituents represents an open system (i.e., if fluid constituents are present in the mixture, either explicitly or implicitly), then \(\rho_{r}\) is no longer necessarily constant, since there may be mass exchange between fluid and solid constituents.

Mixture Free Energy and Stress

The referential free energy density of the mixture is obtained as

\[ \begin{equation} \Psi_{r}=\rho_{r}\psi=\sum_{\sigma}\rho_{r}^{\sigma}\psi^{\sigma}\,,\label{eq:general-free-energy} \end{equation} \]

where \(\psi^{\sigma}\) is the specific free energy of constituent \(\sigma\) and \(\psi\) is the specific free energy of the mixture. An important function of state which arises later in our treatment is the chemical potential of constituent \(\sigma\), given by

\[ \begin{equation} \mu^{\sigma}=\frac{\partial\Psi_{r}}{\partial\rho_{r}^{\sigma}}\,.\label{eq:general-chemical-potential} \end{equation} \]

The mixture Cauchy stress is given by

\[ \begin{equation} \boldsymbol{\sigma}=\frac{1}{J^{s}}\frac{\partial\Psi_{r}}{\partial\mathbf{F}^{s}}\cdot\left(\mathbf{F}^{s}\right)^{T}=\frac{2}{J^{s}}\mathbf{F}^{s}\cdot\frac{\partial\Psi_{r}}{\partial\mathbf{C}^{s}}\cdot\left(\mathbf{F}^{s}\right)^{T}\,,\label{eq:general-mixture-stress} \end{equation} \]

where \(\mathbf{C}^{s}=\mathbf{F}^{s}\cdot\left(\mathbf{F}^{s}\right)^{T}\) is the right Cauchy-Green tensor. The spatial elasticity tensor may be evaluated from

\[ \begin{equation} \boldsymbol{\mathcal{C}}=\frac{4}{J^{s}}\left(\mathbf{F}^{s}\oslash\mathbf{F}^{s}\right):\frac{\partial^{2}\Psi_{r}}{\partial\mathbf{C}^{s}\partial\mathbf{C}^{s}}:\left(\left(\mathbf{F}^{s}\right)^{T}\oslash\left(\mathbf{F}^{s}\right)^{T}\right)\,.\label{eq:general-mixture-elasticity} \end{equation} \]

Closed System of Solid Constituents

In reactive frameworks where \(\rho_{r}^{\sigma}\) evolves according to eq.\eqref{eq:mass-balance-alpha-1}, and where \(\rho_{r}\) remains constant due to the fact that the solid mixture represents a closed system, it may be convenient to define the mass fraction

\[ \begin{equation} w^{\sigma}=\frac{\rho_{r}^{\sigma}}{\rho_{r}}\,,\label{eq:mass-fraction} \end{equation} \]

in which case we may rewrite eq.\eqref{eq:general-free-energy} as \(\Psi_{r}=\sum_{\sigma}w^{\sigma}\Psi_{r}^{\sigma}\) where \(\Psi_{r}^{\sigma}\equiv\rho_{r}\psi^{\sigma}\) is the referential strain energy density of solid \(\sigma\) under the assumption that it is the sole mixture constituent, normalized by the referential volume of the master constituent \(s\). Based on eq.\eqref{eq:general-mass-balance} the mass fractions satisfy \(\sum_{\sigma}w^{\sigma}=1\). In this case, when \(\mathbf{X}^{\sigma}\ne\mathbf{X}^{s}\) and \(\psi^{\sigma}\) is most conveniently expressed as a function of \(\mathbf{F}^{\sigma}\), we may use eq.\eqref{eq:general-total-def-grad} to evaluate \(\partial\mathbf{F}^{s}/\partial\mathbf{F}^{\sigma}=\mathbf{I}\oslash\left(\mathbf{F}^{\sigma s}\right)^{-T}\) and calculate the mixture stress using the alternative form

\[ \begin{equation} \boldsymbol{\sigma}=\frac{1}{J^{s}}\sum_{\sigma}w^{\sigma}\frac{\partial\Psi_{r}^{\sigma}}{\partial\mathbf{F}^{\sigma}}\cdot\left(\mathbf{F}^{\sigma}\right)^{T}\,.\label{eq:general-mixture-stress-redux} \end{equation} \]

This expression shows that the mixture stress may evolve not only due to temporal changes in the state of strain but also due to reactive changes in the mass fractions \(w^{\sigma}\).

In FEBio the referential strain energy density for any solid mixture constituent \(\sigma\) is evaluated from the same library of constitutive models used in single-constituent solids. In this library the calculation of the referential strain energy density is based on the assumption that the reference configuration corresponds to the configuration when the deformation gradient passed to those functions is equal to the identity tensor. Thus, passing \(\mathbf{F}^{s}\) to those functions returns the correct \(\Psi_{r}^{\sigma}\). However, when passing \(\mathbf{F}^{\sigma}\) as an argument to those functions, the referential volume is based on the configuration \(\mathbf{X}^{\sigma}\). Let the referential free energy density returned by FEBio for an argument \(\mathbf{F}^{\sigma}\) be denoted by \(\Psi_{0}^{\sigma}\). We may similarly denote the corresponding Cauchy stress and spatial elasticity tensors by \(\boldsymbol{\sigma}_{0}^{\sigma}\) and \(\boldsymbol{\mathcal{C}}_{0}^{\sigma}\). Here, the subscript \(0\) has two meanings: First it emphasizes that the corresponding function is evaluated using \(\mathbf{F}^{\sigma}\) for the mixture constituent \(\sigma\); second it emphasizes that the calculation returns the corresponding measure under the assumption that the mixture consists entirely of that constituent, so that its multiplication by the scale factor \(w^{\sigma}\) returns the actual contribution of that measure to the entire mixture. Based on eq.\eqref{eq:general-jacobian-relation} we find that \(\Psi_{r}^{\sigma}=J^{\sigma s}\Psi_{0}^{\sigma}\) so that the mixture free energy density may be evaluated from

\[ \begin{equation} \Psi_{r}=\sum_{\sigma}w^{\sigma}J^{\sigma s}\Psi_{0}^{\sigma}\left(\mathbf{F}^{\sigma}\right)\,,\label{eq:febio-mixture-sed} \end{equation} \]

whereas the mixture Cauchy stress and spatial elasticity are given by

\[ \begin{equation} \boldsymbol{\sigma}=\sum_{\sigma}w^{\sigma}\boldsymbol{\sigma}_{0}^{\sigma}\,,\quad\boldsymbol{\mathcal{C}}=\sum_{\sigma}w^{\sigma}\boldsymbol{\mathcal{C}}_{0}^{\sigma}\,,\label{eq:febio-mixture-stress} \end{equation} \]

where

\[ \begin{equation} \boldsymbol{\sigma}_{0}^{\sigma}=\frac{1}{J^{\sigma}}\frac{\partial\Psi_{0}^{\sigma}}{\partial\mathbf{F}^{\sigma}}\cdot\left(\mathbf{F}^{\sigma}\right)^{T}=\frac{1}{J^{s}}\frac{\partial\Psi_{r}^{\sigma}}{\partial\mathbf{F}^{\sigma}}\cdot\left(\mathbf{F}^{\sigma}\right)^{T}\label{eq:febio-stress-function} \end{equation} \]

and

\[ \begin{equation} \boldsymbol{\mathcal{C}}_{0}^{\sigma}=\frac{4}{J^{\sigma}}\left(\mathbf{F}^{\sigma}\oslash\mathbf{F}^{\sigma}\right):\frac{\partial^{2}\Psi_{0}^{\sigma}}{\partial\mathbf{C}^{\sigma}\partial\mathbf{C}^{\sigma}}:\left(\left(\mathbf{F}^{\sigma}\right)^{T}\oslash\left(\mathbf{F}^{\sigma}\right)^{T}\right)\,.\label{eq:febio-elasticity-function} \end{equation} \]

These relations show that the Cauchy stress and spatial elasticity tensor of constituent \(\sigma\) in eq.\eqref{eq:febio-mixture-stress} may be evaluated using existing FEBio functions without needing to adjust for the choice of reference configuration \(\mathbf{X}^{s}\) or \(\mathbf{X}^{\sigma}\), as evidenced by eq.\eqref{eq:febio-stress-function} in the case of the stress. However the referential strain energy density needs to be properly scaled by \(J^{\sigma s}\) as shown in eq.\eqref{eq:febio-mixture-sed}.

The calculation of the 2nd Piola-Kirchhoff stress \(\mathbf{S}\) for each generation \(\sigma\) is more elaborate. When using the master reference configuration \(\mathbf{X}^{s}\), this stress is given by

\[ \mathbf{S}^{\sigma}\left(\mathbf{F}^{s}\right)=2\frac{\partial\Psi_{r}^{\sigma}\left(\mathbf{F}^{s}\right)}{\partial\mathbf{C}^{s}}\,. \]

When using the reference configuration \(\mathbf{X}^{\sigma}\), the stress is evaluated from a similar relation

\[ \mathbf{S}_{0}^{\sigma}\left(\mathbf{F}^{\sigma}\right)=2\frac{\partial\Psi_{0}^{\sigma}\left(\mathbf{F}^{\sigma}\right)}{\partial\mathbf{C}^{\sigma}}\,. \]

It can be shown that these stresses are related according to

\[ \mathbf{S}^{\sigma}\left(\mathbf{F}^{s}\right)=J^{\sigma s}\left(\mathbf{F}^{\sigma s}\right)^{-1}\cdot\mathbf{S}_{0}^{\sigma}\left(\mathbf{F}^{\sigma}\right)\cdot\left(\mathbf{F}^{\sigma s}\right)^{-T}\,. \]

FEBio does not use this calculation for solid mixtures, as all internal calculations employ the Cauchy stress.

Open System of Solid Constituents

When \(\rho_{r}\) as defined in eq.\eqref{eq:general-mass-balance} does not remain constant due to the implicit or explicit presence of fluid constituents, we may choose to define the mass fraction \(\omega^{\sigma}\) of constituent \(\sigma\) based on the true density \(\rho_{0}^{\sigma}\) of the solid constituent,

\[ \begin{equation} \omega^{\sigma}=\frac{\rho_{r}^{\sigma}}{\rho_{0}^{\sigma}}\,,\label{eq:mass-fraction-1} \end{equation} \]

In this case, \(\omega^{\sigma}=1\) when the solid mixture consists entirely of constituent \(\sigma\). Unlike \(w^{\sigma}\) in eq.\eqref{eq:general-mass-balance}, the summation of \(\omega^{\sigma}\) over all \(\sigma\) is meaningless here, since the denominator in eq.\eqref{eq:mass-fraction-1} is not common to all \(\sigma\). Nevertheless, the mixture free energy density in eq.\eqref{eq:general-free-energy} may be rewritten conveniently as \(\Psi_{r}=\sum_{\sigma}\omega^{\sigma}\Psi_{0}^{\sigma}\) where \(\Psi_{0}^{\sigma}\equiv\rho_{0}^{\sigma}\psi^{\sigma}\) represents the free energy density of pure solid \(\sigma\). Here again, \(\psi^{\sigma}\) is most conveniently expressed as a function of \(\mathbf{F}^{\sigma}\), so that eq.\eqref{eq:general-total-def-grad} may be use to evaluate \(\partial\mathbf{F}^{s}/\partial\mathbf{F}^{\sigma}=\mathbf{I}\oslash\left(\mathbf{F}^{\sigma s}\right)^{-T}\) and calculate the mixture stress using the alternative form

\[ \begin{equation} \boldsymbol{\sigma}=\frac{1}{J^{s}}\sum_{\sigma}\omega^{\sigma}\frac{\partial\Psi_{0}^{\sigma}}{\partial\mathbf{F}^{\sigma}}\cdot\left(\mathbf{F}^{\sigma}\right)^{T}\,.\label{eq:general-mixture-stress-redux-1} \end{equation} \]

This expression shows that the mixture stress may evolve not only due to temporal changes in the state of strain but also due to reactive changes in the mass fractions \(\omega^{\sigma}\). In this type of open-system formulation it becomes the user's responsibility to ensure that all other \(\omega^{\sigma}\) values reduce to zero when one of the \(\omega^{\sigma}\) reaches unity.

Frame Indifference

To maintain frame indifference for constrained mixtures we must satisfy \(\boldsymbol{\sigma}\left(\mathbf{Q}\cdot\mathbf{F}^{s}\right)=\mathbf{Q}\cdot\boldsymbol{\sigma}\left(\mathbf{F}^{s}\right)\cdot\mathbf{Q}^{T}\), where \(\mathbf{Q}\) is an orthogonal transformation that maintains the symmetry group of \(\boldsymbol{\sigma}\) ⁵. Here, we abbreviated the list of state variables for \(\boldsymbol{\sigma}\) to just \(\mathbf{F}^{s}\), for notational simplicity. This frame indifference is automatically satisfied for \(\boldsymbol{\sigma}\) based on the hyperelasticity relation of \eqref{eq:general-mixture-stress} and the invariance of \(\psi\) to \(\mathbf{Q}\). By the same argument, \(\boldsymbol{\sigma}_{0}^{\sigma}\) as given in \eqref{eq:febio-stress-function} satisfies \(\boldsymbol{\sigma}_{0}^{\sigma}\left(\mathbf{Q}\cdot\mathbf{F}^{\sigma}\right)=\mathbf{Q}\cdot\boldsymbol{\sigma}_{0}^{\sigma}\left(\mathbf{F}^{\sigma}\right)\cdot\mathbf{Q}^{T}\) for transformations that maintain the symmetry group of \(\psi^{\sigma}\). Since the symmetry transformations \(\mathbf{Q}\) of \(\psi\) must also belong to the symmetry groups of all \(\psi^{\sigma}\), we need to ensure that \(\mathbf{F}^{s*}=\mathbf{Q}\cdot\mathbf{F}^{s}\) and \(\mathbf{F}^{\sigma*}=\mathbf{Q}\cdot\mathbf{F}^{\sigma}\) also satisfy the kinematic constraint of \eqref{eq:general-total-def-grad} for those transformations \(\mathbf{Q}\), namely \(\mathbf{F}^{s*}=\mathbf{F}^{\sigma*}\cdot\mathbf{F}^{\sigma s}\). This can be achieved if and only if \(\mathbf{F}^{\sigma s}\) is invariant under a change of frame. It follows from this argument that the right Cauchy-Green tensors are invariant, \(\mathbf{C}^{s*}=\left(\mathbf{F}^{s*}\right)^{T}\cdot\mathbf{F}^{s*}=\left(\mathbf{F}^{s}\right)^{T}\cdot\mathbf{F}^{s}=\mathbf{C}^{s}\) and similarly \(\mathbf{C}^{\sigma*}=\mathbf{C}^{\sigma}\), while also satisfying \(\mathbf{C}^{s*}=\mathbf{C}^{s}=\left(\mathbf{F}^{\sigma s}\right)^{T}\cdot\mathbf{C}^{\sigma}\cdot\mathbf{F}^{\sigma s}\). Right stretch tensors \(\mathbf{U}^{s}\) and \(\mathbf{U}^{\sigma}\) are also invariant.

Simple Solid Mixtures

In the simplest type of non-reactive constrained mixtures of solids all constituents \(\sigma\) share the same reference configuration \(\mathbf{X}^{\sigma}=\mathbf{X}^{s}\), in which case \(\mathbf{F}^{\sigma}=\mathbf{F}^{s}\) and \(J^{\sigma s}=1\). For this type of mixture where \(\rho_{r}^{\sigma}\) does not evolve, it is convenient to set \(w^{\sigma}=1\) for all \(\sigma\) and scale the material properties of \(\Psi_{r}^{\sigma}\) to properly reflect the contribution of each constituent \(\sigma\) to the mixture. Thus, the relation of eq.\eqref{eq:febio-mixture-stress} reduces to

\[ \begin{equation} \boldsymbol{\sigma}=\sum_{\sigma}\boldsymbol{\sigma}_{0}^{\sigma}\,.\label{eq:simple-mixture-stress} \end{equation} \]

This type of non-reactive mixture of solids is represented in FEBio using the material “solid mixture”.

When each constituent of a solid mixture is modeled using an uncoupled formulation to enforce nearly isochoric responses, the strain energy density of this type of mixture is given by

\[ \begin{equation} \Psi_{r}=U\left(J\right)+\sum\limits_{\sigma}\tilde{\Psi}_{r}^{\sigma}\left(\mathbf{\tilde{F}}^{s},\rho_{r}^{\sigma}\right)\,,\label{eq140} \end{equation} \]

where \(U\left(J\right)\) is the volumetric energy component, \(\tilde{\Psi}_{r}=\sum\nolimits_{\sigma}\tilde{\Psi}_{r}^{\sigma}\) is the distortional energy component, and \(\mathbf{\tilde{F}}^{s}\) is the distortional part of the deformation gradient, as described in Section Nearly-Incompressible Hyperelasticity. This type of non-reactive mixture of uncoupled solids is represented in FEBio using the material “uncoupled solid mixture”.

Multigenerational Interstitial Growth

Multigenerational interstitial growth mechanics may be modeled using a reactive mixture of constrained solids as described in ⁶. In this framework it is assumed that a porous solid matrix may gain mass via interstitial growth, such that the porosity of the solid decreases with increasing solid mass content. The history of growth is discretized temporally into generations \(\sigma\) such that the mass added in the time interval \(t^{\sigma}\le t<t^{\sigma+1}\) has a reference configuration \(\mathbf{X}^{\sigma}\). This model assumes that each generation \(\sigma\) gets deposited into the existing mixture in a stress-free state. This can be achieved by adopting the constitutive assumption that \(\mathbf{X}^{\sigma}=\boldsymbol{\chi}^{s}\left(\mathbf{X}^{s},t^{\sigma}\right)\). In other words, the reference configuration of generation \(\sigma\) is the current configuration of the solid mixture at time \(t^{\sigma}\). An alternative form of this constitutive model, which satisfies frame-indifference, is that \(\mathbf{F}^{\sigma}\left(\mathbf{X}^{s},t^{\sigma}\right)=\mathbf{R}^{s}\left(\mathbf{X}^{s},t^{\sigma}\right)\) at the time the new generation \(\sigma\) is being deposited, where \(\mathbf{R}^{s}\) is the rotation tensor in the polar decomposition of \(\mathbf{F}^{s}=\mathbf{R}^{s}\cdot\mathbf{U}^{s}\). Equivalently, \(\mathbf{F}^{\sigma s}\left(\mathbf{X}^{s}\right)=\mathbf{U}^{s}\left(\mathbf{X}^{s},t^{\sigma}\right)\) according to eq.\eqref{eq:general-total-def-grad}. Since \(\mathbf{F}^{\sigma s}\) is a right-stretch tensor, it satisfies frame indifference.

The underlying assumption of multigenerational growth is that \(\mathbf{U}^{s}\left(\mathbf{X}^{s},t^{\sigma}\right)\) evolves for each generation \(\sigma\) either due to load-induced deformations or deformations produced by swelling processes, such as cell growth or Donnan swelling.

This type of multigenerational growth material is implemented in FEBio as “multigeneration” for mixtures of elastic solids, and as “multiphasic-multigeneration” for reactive multiphasic mixtures whose solid constituent is a multigenerational growth material. In this type of material the user needs to prescribe the generation birth times \(t^{\sigma}\) and the properties of the material of each generation \(\sigma\). The code automatically prescribes \(\mathbf{F}^{\sigma s}\) based on the constitutive model given above, for each generation \(\sigma\). For example, one may use a material model whose response depends on the evolving composition \(\rho_{r}^{\sigma}\) of that generation. The composition \(\rho_{r}^{\sigma}\) may evolve either due to chemical reactions modeled in a multiphasic framework, or by associating a user-defined load curve with \(\rho_{r}^{\sigma}\). Thus, the material properties (and the stress response) of generation \(\sigma\) need not remain constant over the generation time interval \(t^{\sigma}\le t<t^{\sigma+1}\), even though \(\mathbf{F}^{\sigma s}\) remains constant during that generation.

When the growth process is negative (\(\hat{\rho}_{r}^{\sigma}<0\)), it implies that the solid constituent \(\sigma\) is losing mass (solid resorption); this loss of mass terminates when \(\rho_{r}^{\sigma}=0\). A material model that depends on \(\rho_{r}^{\sigma}\) may exhibit evolving material properties during this resorption process until that generation produces zero stress when \(\rho_{r}^{\sigma}=0\).

As a result of these evolving growth processes, this multigeneration mixture may exhibit residual stresses.

Prescribed Pre-Stretch

An alternative approach to multigenerational growth is to prescribe \(\mathbf{F}^{\sigma s}\) as a user-defined function. For isotropic stretching we may define

\[ \begin{equation} \mathbf{F}^{\sigma s}=\lambda^{\sigma s}\mathbf{I}\,,\label{eq:prescribed-defgrad-isotropic} \end{equation} \]

where \(\lambda^{\sigma s}>0\) is the stretch ratio. A value greater than unity causes swelling whereas a value less than unity causes contraction; \(\lambda^{\sigma s}=1\) produces \(\mathbf{F}^{\sigma s}=\mathbf{I}\), which recovers the simple solid mixture described in Section Simple Solid Mixtures. In FEBio a load curve may be associated with \(\lambda^{\sigma s}\) to ramp up the prescribed deposition stretch.

For orthotropic stretching we define mutually orthogonal symmetry planes with unit normals \(\mathbf{a}_{i}\) (\(i=1,2,3\) and \(\mathbf{a}_{i}\cdot\mathbf{a}_{j}=\delta_{ij}\)). Then,

\[ \begin{equation} \mathbf{F}^{\sigma s}=\sum_{i=1}^{3}\lambda_{i}^{\sigma s}\mathbf{a}_{i}\otimes\mathbf{a}_{i}\,,\label{eq:prescribed-defgrad-orthotropic} \end{equation} \]

where \(\lambda_{i}^{\sigma s}\) is the stretch ratio for the prescribed stretch along \(\mathbf{a}_{i}\); in general, \(\lambda_{1}^{\sigma s}\ne\lambda_{2}^{\sigma s}\ne\lambda_{3}^{\sigma s}\), though this model can be specialized to transversely isotropic stretch by setting two of these stretch ratios equal to each other.

At lower symmetries the constitutive model for \(\mathbf{F}^{\sigma s}\) must necessarily combine stretch and rotation. For monoclinic materials, deformations may be prescribed along three unit vectors \(\mathbf{a}_{i}\) that satisfy \(\mathbf{a}_{1}\cdot\mathbf{a}_{2}=0\), \(\mathbf{a}_{1}\cdot\mathbf{a}_{3}=0\), and \(\mathbf{a}_{2}\cdot\mathbf{a}_{3}\ne0\), such that \(\mathbf{a}_{1}\) defines the single plane of symmetry. In this case,

\[ \begin{equation} \begin{aligned}\mathbf{F}^{\sigma s} & =\lambda_{1}^{\sigma s}\mathbf{a}_{1}\otimes\mathbf{a}_{1}+\frac{\lambda_{2}^{\sigma s}}{1-\alpha_{23}^{2}}\mathbf{a}_{2}\otimes\left(\mathbf{a}_{2}-\alpha_{23}\mathbf{a}_{3}\right)\\ & +\frac{\lambda_{3}^{\sigma s}}{1-\alpha_{23}^{2}}\mathbf{a}_{3}\otimes\left(\mathbf{a}_{3}-\alpha_{23}\mathbf{a}_{2}\right) \end{aligned} \,,\label{eq:prescribed-defgrad-monoclinic} \end{equation} \]

where \(\lambda_{i}^{\sigma s}\) is the stretch along \(\mathbf{a}_{i}\) and \(\alpha_{23}=\mathbf{a}_{2}\cdot\mathbf{a}_{3}\). Under general conditions (i.e., when \(\lambda_{1}^{\sigma s}\ne\lambda_{2}^{\sigma s}\ne\lambda_{3}^{\sigma s}\) and \(\alpha_{23}\ne0\)), the polar decomposition theorem shows that this deformation gradient is not a pure stretch, as it also involves a rotation.

The deformation gradient for expansion of a triclinic material may be similarly constructed by finding \(\mathbf{F}^{\sigma s}\) such that \(\mathbf{F}^{\sigma s}\cdot\mathbf{a}_{i}=\lambda_{i}^{\sigma s}\mathbf{a}_{i}\) (no sum) for non-orthogonal and non-collinear unit vectors \(\mathbf{a}_{i}\),

\[ \begin{equation} \begin{aligned}\mathbf{F}^{\sigma s} & =\frac{\lambda_{1}^{\sigma s}}{d}\mathbf{a}_{1}\otimes\left(\left(1-\alpha_{23}^{2}\right)\mathbf{a}_{1}-\left(\alpha_{12}-\alpha_{13}\alpha_{23}\right)\mathbf{a}_{2}-\left(\alpha_{13}-\alpha_{12}\alpha_{23}\right)\mathbf{a}_{3}\right)\\ & +\frac{\lambda_{2}^{\sigma s}}{d}\mathbf{a}_{2}\otimes\left(-\left(\alpha_{12}-\alpha_{13}\alpha_{23}\right)\mathbf{a}_{1}+\left(1-\alpha_{13}^{2}\right)\mathbf{a}_{2}-\left(\alpha_{23}-\alpha_{12}\alpha_{13}\right)\mathbf{a}_{3}\right)\\ & +\frac{\lambda_{3}^{\sigma s}}{d}\mathbf{a}_{3}\otimes\left(-\left(\alpha_{13}-\alpha_{12}\alpha_{23}\right)\mathbf{a}_{1}-\left(\alpha_{23}-\alpha_{12}\alpha_{13}\right)\mathbf{a}_{2}+\left(1-\alpha_{12}^{2}\right)\mathbf{a}_{3}\right) \end{aligned} \label{eq:prescribed-defgrad-triclinic} \end{equation} \]

where

\[ \begin{equation} d=1-\alpha_{12}^{2}-\alpha_{13}^{2}-\alpha_{23}^{2}+2\alpha_{23}\alpha_{13}\alpha_{12}\label{eq:triclinic-defgrad-denominator} \end{equation} \]

and \(\alpha_{ij}=\mathbf{a}_{i}\cdot\mathbf{a}_{j}\) is the cosine of the angle between \(\mathbf{a}_{i}\) and \(\mathbf{a}_{j}\). All of these constitutive models for \(\mathbf{F}^{\sigma s}\) satisy frame indifference since they depend on material vectors \(\mathbf{a}_{i}\).

This type of elastic solid with prescribed pre-stretch is modeled in FEBio using the materials “prestretch elastic” and “prestretch uncoupled elastic”. In the uncoupled version each solid constituent in the mixture has a deformation gradient \(\mathbf{F}^{\sigma}\) which produces nearly isochoric responses (\(J^{\sigma}=1\)), whereas no constraint is placed on the volumetric strain of the pre-stretch \(\mathbf{F}^{\sigma s}\).

---

1. Ateshian, G. A.; Weiss, J. A.. "Anisotropic hydraulic permeability under finite deformation." Journal of biomechanical engineering, vol. 132, pp. 111004 (2010). ↩↩

2. Nims, Robert J; Ateshian, Gerard A. "Reactive constrained mixtures for modeling the solid matrix of biological tissues." Journal of Elasticity, vol. 129, pp. 69--105 (2017). ↩↩

3. Bowen, Ray M. "The thermochemistry of a reacting mixture of elastic materials with diffusion." Archive for Rational Mechanics and Analysis, vol. 34, pp. 97--127 (1969). ↩

4. Ateshian, G. A.. "On the theory of reactive mixtures for modeling biological growth." Biomech Model Mechanobiol, vol. 6, pp. 423-45 (2007). ↩

5. Ateshian, Gerard A; Zimmerman, Brandon K. "Continuum Thermodynamics of Constrained Reactive Mixtures." J Biomech Eng, vol. 144 (2022). ↩

6. Ateshian, G. A.; Ricken, T.. "Multigenerational interstitial growth of biological tissues." Biomech Model Mechanobiol, vol. 9, pp. 689-702 (2010). ↩