\chapter{An Eulerian perspective}
\label{eulerian-perspective}

As detailed at the outset of Chapter~\ref{lagrangian-perspective}, the
continuum treatment presented thus far has stemmed from classical
theories for solid continua, which are traditionally formulated in a
Lagrangian setting. Since our continuum idealisation of tissues
includes fluid components, and material coordinates are, in general,
not known in fluid mechanics, this chapter revisits the derivation of
the governing field equations of growing tissues following an Eulerian
approach.

In this {\em spatial description}, attention is turned to a region of
space coinciding with the {\em current configuration} of the tissue,
where the evolution of field variables of interest is studied.
Remarkably, this dissimilar approach results in a set of balance
equations which is completely equivalent to that deduced in Section%
~\ref{balance-laws}, just {\em pushed-forward} to the current
configuration. But more significantly, the spatial approach presented
below naturally leads to the identification of a different set of
primitive variables more suitable to physically relevant boundary
value problems.

This chapter is divided into two main parts. The first part
(Section~\ref{eu-balance-laws}) briefly recapitulates the fundamental
quantities characterising the tissue, pointing out noteworthy
differences from Section~\ref{balance-laws}, and derives the balance
laws governing finite deformation growth in terms of spatial
quantities. The second (Section~\ref{eu-entropy-inequality}), closes
this set of balance laws by deducing a set of physiologically relevant
constitutive relationships which ensure that the entropy production
inequality is satisfied a priori.

\section{Balance laws for an open mixture}
\label{eu-balance-laws}

We initiate this discussion with the introduction of $\Omega_{t}$, a
temporally-varying closure of an open set in $\mathbb{R}^{3}$ with a
piecewise smooth boundary, which we define to be our region of
interest. This region of interest is constructed to coincide with the
current position of the solid component of the tissue, and our primary
interest lies in the evolution of various field variables inside
$\Omega_{t}$, as observed from an inertial reference frame
(\citet{newton1726}, Corollary~V, p.~423).

\begin{figure}
  \centering
  \psfrag{A}{$\bX$}
  \psfrag{B}{$\Omega_0$}
  \psfrag{C}{$\Bvarphi_{t}$}
  \psfrag{D}{$\bx$}
  \psfrag{E}{$\Omega_t$}
  \psfrag{F}{$\bX^{\iota}$}
  \psfrag{G}{$\Omega_{0}^{\iota}$}
  \psfrag{H}{$\Bvarphi^{\iota}_{t}$}
  \includegraphics[width=0.8\textwidth]{images/elucidation/%
    continuum-potato-current-configuration}
  \caption{An Eulerian point of view.}
  \label{continuum-potato-current-configuration}
\end{figure}

It is assumed that the there exists a {$\mathit{C}^{2}$} (in space and
time), bijective and orientation preserving map $\Bvarphi_{t}(\bX):
\mathbb{R}^3\times\mathbb{R}^{+}\cup\{0\} \rightarrow\mathbb{R}^{3}$
such that $\Omega_{t} = \Bvarphi_{t} (\Omega_{0})$ for some
convenient, {\em fixed} subset of $\mathbb{R}^{3}$, $\Omega_{0}$. This
ensures that $\Omega_{t}$ evolves in a well-behaved manner,
disallowing non-physical deformations from being imposed on the
tissue, and furthermore, affords the application of Reynolds'
Transport Theorem (Appendix~\ref{reynolds-transport}). The map
$\Bvarphi_{t}$ is visualised in Figure~\ref{continuum-potato%
  -current-configuration}.

Denoting a point in $\Omega_{t}$ by $\bx$, $\bv(\bx,t) =
\frac{\partial \Bvarphi_{t}}{\partial t}$ defines the {\em spatial
  velocity} of the system domain. In contrast to
Section~\ref{balance-laws}, it is recognised at the outset that each
species of the tissue is capable of undergoing its own motion
independent of the tissue's solid component; and so we introduce the
spatial velocity of an arbitrary species~$\iota$,
$\bv^{\iota}(\bx,t)$, which is assumed to be $\mathit{C}^{1}$ in space
and time.\footnote{It is important to note that these quantities are
  primitive variables in themselves. While these species velocities
  can be formally understood as \mbox{$\bv^{\iota}(\bx,t) =
    \frac{\partial \Bvarphi^{\iota}_{t}}{\partial t}$}, where
  $\Bvarphi^{\iota}_{t}$ is the deformation map of each
  species~$\iota$ from an arbitrary reference configuration to
  $\Omega_{t}$ (shown in the portion of Figure~\ref{continuum-potato%
    -current-configuration} constructed using dashed-lines), the
  quantities $\Bvarphi^{\iota}_{t}$ are neither explicitly defined nor
  tracked.}  Unlike the corresponding quantities defined in
Section~\ref{balance-of-linear-momentum}, these velocities are defined
to be the {\em total} velocities of each species, not velocities
relative to the solid component of the tissue.\footnote{And, with the
  redefinition of these quantities as total velocities, the velocity
  of the solid phase ceases to be special. This is manifest in the
  forms of the balance laws deduced in Sections~\ref{eu-balance-of%
    -mass}--\ref{eu-balance-of-energy}.}

The following discussion of balance laws is carried out entirely in
terms of an arbitrary species~$\iota$, and specialisation to the solid
collagenous, fluid and solute phases (i.e. \mbox{$\iota=\mathrm{c},
  \mathrm{f}, \mathrm{s}$}) is reintroduced only in Section~\ref{eu%
  -entropy-inequality} when discussing constitutive relationships
specific to these different components of the tissue.

\subsection{Balance of mass}
\label{eu-balance-of-mass}

We now turn our attention to the evolution of the first of our field
variables of interest, the current concentration of an arbitrary
species~$\iota$ constituting the system, $\rho^{\iota}(\bx,t):
\mathbb{R}^3\times\mathbb{R}^{+}\cup\{0\} \rightarrow
\mathbb{R}$. These are defined as the mass of species~$\iota$ per unit
       {\em system volume}, $\Omega_{t}$, and are assumed to be
       $\mathit{C}^{1}$ in time and $\mathit{C}^{2}$ in space. The
       total {\em spatial density} of the tissue can be obtained by
       their summation, $\sum\limits_{\iota}\rho^\iota = \rho$.

With these quantities defined, we have from the conservation of matter
for species $\iota$ over $\Omega_{t}$,

\begin{equation}
\underbrace{\frac{\mathrm{d}} {\mathrm{d}t}\left(\int_{\Omega_{t}}
  \rho^{\iota} \ \mathrm{d}v \right)}_{\text{Rate of change of mass}}
= \underbrace{\int_{\Omega_{t}} \pi^{\iota} \ \mathrm{d}v}_{\text{Mass
    being created}} - \underbrace{\int_{\partial \Omega_{t}}
  \rho^{\iota} \left(\bv^{\iota} - \bv\right) \cdot
  \bn\ \mathrm{d}a,}_{\text{Mass leaving the domain}}
\label{eu-globalbalanceofmass}
\end{equation}

\noindent where $\pi^{\iota}(\bx,t)$ is the volumetric source (or
sink) of species~$\iota$, which specifies the species' mass production
rate per unit system volume, and $\bn$ is the outward normal vector
over $\partial \Omega_{t}$, the boundary of $\Omega_{t}$. At this
point, the only restriction on the source terms, $\pi^{\iota}(\bx,t)$,
is that they be integrable.

On applying Reynolds' Transport Theorem
(Appendix~\ref{reynolds-transport}) to the left hand-side,

\begin{equation*}
\int_{\Omega_{t}} \frac{\partial \rho^{\iota}}{\partial t}
\ \mathrm{d}v + \cancel{\int_{\partial \Omega_{t}} \rho^{\iota} \bv
  \cdot \bn\ \mathrm{d}a} = \int_{\Omega_{t}} \pi^{\iota}
\ \mathrm{d}v - \int_{\partial \Omega_{t}} \rho^{\iota}
\left(\bv^{\iota} - \cancel{\bv}\right) \cdot \bn\ \mathrm{d}a,
\end{equation*}

\noindent Gauss' Divergence Theorem (Appendix~\ref{gauss-divergence})
to the area integral,

\begin{equation*}
\int_{\Omega_{t}} \frac{\partial \rho^{\iota}}{\partial t}
\ \mathrm{d}v = \int_{\Omega_{t}} \pi^{\iota} \ \mathrm{d}v -
\int_{\Omega_{t}} \mathrm{div} \left( \rho^{\iota} \bv^{\iota}\right)
\ \mathrm{d}v,
\end{equation*}

\noindent and localising, we arrive at the final form of the balance
of mass of species $\iota$,

\begin{equation}
\frac{\partial \rho^{\iota}}{\partial t} = \pi^{\iota} -\ \mathrm{div}
\left(\rho^{\iota} \bv^{\iota}\right) \quad \mathrm{in}\ \Omega_{t},
\label{eu-localbalanceofmass}
\end{equation}

\noindent where $\mathrm{div} (\bullet)$ denotes the spatial
divergence operator. This result is consistent with classical mixture
theory \citep{TruesdellToupin:60} and is the current configuration
analogue of Equation~(\ref{localbalanceofmass}).

As in Section~\ref{balance-of-mass}, recall that for an external
observer, the rate of change of mass of the entire system, affected
only by external agents, is independent of interconversion between
species. From the viewpoint of such an observer, the balance of mass
for the tissue as a whole reads,

\begin{equation}
\frac{\mathrm{\mathrm{d}}}{\mathrm{\mathrm{d}}t} \sum\limits_{\iota}
\left(\int_{\Omega_{t}} \rho^{\iota} \ \mathrm{d}v \right) =
-\sum\limits_{\iota} \int_{\partial \Omega_{t}} \rho^{\iota}
\left(\bv^{\iota} - \bv\right) \cdot \bn\ \mathrm{d}a.
\label{eu-systembalanceofmass}
\end{equation}

\noindent Comparing Equation~(\ref{eu-systembalanceofmass}) to a
summation of Equation~(\ref{eu-globalbalanceofmass}) over all species,
both being valid statements of the balance of mass of the tissue as a
whole, it is clear that the sources and sinks satisfy,

\begin{equation}
\sum\limits_{\iota}\pi^{\iota} = 0,
\label{eu-summationrelationmass}
\end{equation}

\noindent which specifies that, during biochemical interactions
between species, the rate of consumption of the reactants equals the
rate of creation of the products.

An identical relation is obtained by \citet{TruesdellToupin:60} and
others who follow their ideas (see, for e.g., \citet{bowen76} and
\citet{passmanetal}), but their deduction of this result, and other
similar results involving quantities internal to the system
(Equations~(\ref{eu-summationrelationmomentum}) and (\ref{eu-sum%
  mationrelationenergy})), stems from a formally stated assumption
that ``the mean response of a heterogeneous mixture obeys the ordinary
equations of a continuum.'' While we do not explicitly employ this
assumption, it is implicit in our assertion that external observers
can quantify field variables characterising species within a system
without being aware of phenomena internal to the system.

\subsection{Balance of linear momentum}
\label{eu-balance-of-linear-momentum}

We now look at how the momentum of a species evolves under the action
of external agents, accounting for mass sources and allowing for the
possibility that matter can leave the domain.

As observed from an inertial reference frame, the balance of momentum
of a species~$\iota$ over $\Omega_{t}$ requires,

\begin{equation}
\begin{split}
\underbrace{\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Omega_{t}}
  \rho^{\iota} \bv^{\iota}\ \mathrm{d}v \right)}_{\text{Rate of change
    of momentum}} = & \underbrace{\int_{\Omega_{t}} \rho^{\iota}
  \left(\bg^{\iota} +
  \bq^{\iota}\right)\ \mathrm{d}v}_{\text{Resultant body force}} +
\underbrace{\int_{\partial \Omega_{t}}
  \Bsigma^{\iota}\bn\ \mathrm{d}a}_{\text{Boundary traction}}\\ + &
\underbrace{\int_{\Omega_{t}}
  \pi^{\iota}\bv^{\iota}\ \mathrm{d}v}_{\text{Momentum being created}}
- \underbrace{\int_{\partial \Omega_{t}} \left(\rho^{\iota}
  \bv^{\iota} \right) \left(\bv^{\iota} - \bv\right) \cdot
  \bn\ \mathrm{d}a,}_{\text{Momentum leaving the domain}}
\end{split}
\label{eu-globalbalanceofmomentum}
\end{equation}

\noindent where, in addition to the quantities introduced previously,
$\bg^{\iota}(\bx,t)$ is the resultant body force of {\em external}
origin acting on species $\iota$, $\bq^{\iota}(\bx,t)$ is the
resultant body force on species $\iota$ from all other species {\em in
  the mixture}, and $\Bsigma^{\iota}(\bx,t)$ is the partial Cauchy
stress on species $\iota$. The partial Cauchy stress tensor
corresponding to species~$\iota$ is the portion of the total stress
borne by the species. All of these quantities are assumed to be
sufficiently smooth. In particular, $\bg^{\iota}$ and $\bq^{\iota}$
are assumed to be integrable, and $\Bsigma^{\iota}$ is assumed to be
$\mathit{C}^{1}$ in space and time.

Application of Reynolds' Transport Theorem (Appendix~\ref{reynolds%
  -transport}) to the left hand-side of Equation~(\ref{eu-%
  globalbalanceofmomentum}) yields,

\begin{equation*}
\begin{split}
\int_{\Omega_{t}} \frac{\partial
  \left(\rho^{\iota}\bv^{\iota}\right)}{\partial t} \ \mathrm{d}v
+\cancel{\int_{\partial \Omega_{t}}
  \left(\rho^{\iota}\bv^{\iota}\right) \bv \cdot \bn\ \mathrm{d}a} = &
\int_{\Omega_{t}} \rho^{\iota} \left(\bg^{\iota}+\bq^{\iota}\right)
\ \mathrm{d}v + \int_{\partial \Omega_{t}}
\Bsigma^{\iota}\bn\ \mathrm{d}a\\ + & \int_{\Omega_{t}}
\pi^{\iota}\bv^{\iota} \ \mathrm{d}v - \int_{\partial \Omega_{t}}
\left(\rho^{\iota} \bv^{\iota} \right)\left(\bv^{\iota} -
\cancel{\bv}\right) \cdot \bn\ \mathrm{d}a.
\end{split}
\end{equation*}

\noindent Using the product rule, Gauss' Divergence Theorem
(Appendix~\ref{gauss-divergence}), and the balance of
mass{\footnote{Recognising that the balance of mass need not be
    satisfied exactly, pointwise, in a numerical implementation.}}
(Equation~\ref{eu-localbalanceofmass}), we have,

\begin{equation*}
\begin{split}
\cancel{\int_{\Omega_{t}} \frac{\partial \rho^{\iota}} {\partial t}
  \bv^{\iota} \ \mathrm{d}v} + \int_{\Omega_{t}} \rho^{\iota}
\frac{\partial \bv^{\iota}} {\partial t} \ \mathrm{d}v = &
\int_{\Omega_{t}} \rho^{\iota} \left(\bg^{\iota}+\bq^{\iota}\right)
\ \mathrm{d}v + \int_{\Omega_{t}}
\mathrm{div}\left(\Bsigma^{\iota}\right)\ \mathrm{d}v\\ + &
\cancel{\int_{\Omega_{t}} \pi^{\iota}\bv^{\iota} \ \mathrm{d}v} -
\int_{\Omega_{t}}\left(\cancel{\mathrm{div}
  \left(\rho^{\iota}\bv^{\iota}\right) \bv^{\iota}} +
\mathrm{grad}\left(\bv^{\iota}\right)
\rho^{\iota}\bv^{\iota}\right)\ \mathrm{d}v,
\end{split}
\end{equation*}

\noindent where $\mathrm{grad} (\bullet)$ denotes the spatial gradient
operator. Upon localisation, we obtain the final form of the balance
of momentum of species $\iota$,

\begin{equation}
\rho^{\iota} \frac{\partial \bv^{\iota}} {\partial t} = \rho^{\iota}
\left(\bg^{\iota}+\bq^{\iota}\right) +
\mathrm{div}\left(\Bsigma^{\iota}\right) -
\mathrm{grad}\left(\bv^{\iota}\right) \rho^{\iota}\bv^{\iota} \quad
\mathrm{in}\ \Omega_{t},
\label{eu-localbalanceofmomentum}
\end{equation}

\noindent which is a result consistent with classical mixture theory
\citep{TruesdellToupin:60} and is the current configuration analogue
of Equation~(\ref{localbalanceofmomentum}).

Neglecting the interaction terms as an external observer, the balance
of momentum for the entire system can be written as follows,

\begin{equation}
\begin{split}
\sum\limits_{\iota}
\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Omega_{t}} \rho^{\iota}
\bv^{\iota} \ \mathrm{d}v \right) = & \sum\limits_{\iota}
\left(\int_{\Omega_{t}} \rho^{\iota} \bg^{\iota} \ \mathrm{d}v +
\int_{\partial \Omega_{t}} \Bsigma^{\iota}\bn\ \mathrm{d}a\right)\\ -
&\sum\limits_{\iota}\int_{\partial \Omega_{t}} \left(\rho^{\iota}
\bv^{\iota} \right) \left(\bv^{\iota} - \bv\right) \cdot
\bn\ \mathrm{d}a.
\end{split}
\label{eu-systembalanceofmomentum}
\end{equation}

\noindent Comparing Equation~(\ref{eu-systembalanceofmomentum}) to a
summation of Equation~(\ref{eu-globalbalanceofmomentum}) over all
species, it is clear that the sources and interaction forces satisfy
the relation,

\begin{equation}
\sum\limits_{\iota}\left( \rho^{\iota}\bq^{\iota} + \pi^{\iota}
\bv^{\iota} \right) = 0,
\label{eu-summationrelationmomentum}
\end{equation}

\noindent which states that the momentum being introduced at a point
due to the creation of matter has to be negated by momentum
interactions with other species; ensuring that there is no mechanism
for momentum production internal to the system.

\subsection{Balance of angular momentum}
\label{eu-balance-of-angular-momentum}

Consider the position vector $\bp(\bx)$ of a point on the tissue
relative to a fixed point\footnote{Which may or may not be the origin
  of the system's Euclidean space.} in space. The balance of angular
momentum about $\bp$, as observed from an inertial reference frame, of
a species~$\iota$ over $\Omega_{t}$ requires,

\begin{equation}
\begin{split}
\underbrace{\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Omega_{t}} \bp
  \times \rho^{\iota} \bv^{\iota}\ \mathrm{d}v\right)}_{\text{Rate of
    change of angular momentum}} = & \underbrace{\int_{\Omega_{t}} \bp
  \times \rho^{\iota} \left(\bg^{\iota} +
  \bq^{\iota}\right)\ \mathrm{d}v}_{\text{Moment from body forces}} +
\underbrace{\int_{\partial \Omega_{t}} \bp \times
  \left(\Bsigma^{\iota}\bn\right)\ \mathrm{d}a}_{\text{Moment from
    traction}}\\ + & \underbrace{\int_{\Omega_{t}} \bp \times
  \pi^{\iota}\bv^{\iota}\ \mathrm{d}v}_{\text{Angular momentum being
    created}} \\ - & \underbrace{\int_{\partial \Omega_{t}} \left(\bp
  \times \rho^{\iota} \bv^{\iota}\right) \left(\bv^{\iota} -
  \bv\right) \cdot \bn\ \mathrm{d}a,}_{\text{Angular momentum leaving
    the domain}}
\end{split}
\label{eu-globalbalanceofangularmomentum}
\end{equation}

\noindent since it is reasonable to assume that the material
comprising the tissue is not a {\em polar material}.\footnote{Some
  materials, such as liquid crystals, become polarised under the
  presence of electric fields and consequently have additional global
  torque contributions to their balance of momentum equations
  \citep{TruesdellNoll:65}.}

On applying Reynolds' Transport Theorem (Appendix~\ref{reynolds%
  -transport}), Equation~(\ref{eu-globalbalanceofangularmomentum})
reduces to,

\begin{equation*}
\begin{split}
\int_{\Omega_{t}} \frac{\partial}{\partial t}\left( \bp \times
\rho^{\iota} \bv^{\iota}\right)\ \mathrm{d}v \quad = &
\int_{\Omega_{t}} \bp \times \rho^{\iota} \left(\bg^{\iota} +
\bq^{\iota}\right)\ \mathrm{d}v + \int_{\partial \Omega_{t}} \bp
\times \left(\Bsigma^{\iota}\bn\right)\ \mathrm{d}a\\ + &
\int_{\Omega_{t}} \bp \times \pi^{\iota}\bv^{\iota}\ \mathrm{d}v -
\int_{\partial \Omega_{t}} \left(\bp \times \rho^{\iota}
\bv^{\iota}\right) \bv^{\iota} \cdot \bn\ \mathrm{d}a.
\end{split}
\end{equation*}

Using Gauss' Divergence Theorem (Appendix~\ref{gauss-divergence}) and
the product rule, we have the following relations:

\begin{displaymath}
\int_{\partial \Omega_{t}} \bp \times
\left(\Bsigma^{\iota}\bn\right)\ \mathrm{d}a = \int_{\Omega_{t}} \bp
\times \mathrm{div}\left(\Bsigma^{\iota}\right)\ \mathrm{d}v
+\int_{\Omega_{t}} \Bepsilon\colon\Bsigma^{\iota^{\mathrm{T}}}
\ \mathrm{d}v,
\end{displaymath}

\noindent where $\Bepsilon$ is the permutation symbol (introduced in
Section~\ref{balance-of-angular-momentum}),

\begin{displaymath}
\int_{\Omega_{t}} \frac{\partial}{\partial t}\left( \bp \times
\rho^{\iota} \bv^{\iota}\right)\ \mathrm{d}v = \int_{\Omega_{t}}
\frac{\partial \rho^{\iota}} {\partial t} \bp \times
\bv^{\iota}\ \mathrm{d}v + \int_{\Omega_{t}} \rho^{\iota} \bp \times
\frac{\partial \bv^{\iota}} {\partial
  t}\ \mathrm{d}v,\quad\mathrm{and}
\end{displaymath}

\begin{displaymath}
\int_{\partial \Omega_{t}} \left(\bp \times \rho^{\iota}
\bv^{\iota}\right) \bv^{\iota} \cdot \bn\ \mathrm{d}a =
\int_{\Omega_{t}}\bp \times \bv^{\iota} \mathrm{div}
\left(\rho^{\iota}\bv^{\iota}\right)\ \mathrm{d}v +
\int_{\Omega_{t}}\bp \times
\left(\mathrm{grad}\left(\bv^{\iota}\right)
\rho^{\iota}\bv^{\iota}\right)\ \mathrm{d}v,
\end{displaymath}

\noindent since $\rho^{\iota} \bv^{\iota}\times\bv^{\iota} =
0$. Substituting these relations above and invoking the balance of
mass (\ref{eu-localbalanceofmass}) and the balance of linear momentum
(\ref{eu-localbalanceofmomentum}), we obtain on localisation that,

\begin{equation}
\Bepsilon\colon\Bsigma^{\iota^{\mathrm{T}}} = 0,
\label{eu-localbalanceofangularmomentum}
\end{equation}

\noindent or the partial Cauchy stress tensor, $\Bsigma^{\iota}$, is
symmetric. This classical result is the pushed-forward form of the
synonymous result derived earlier (Section~\ref{balance-of-angular%
  -momentum}) in terms of the partial first Piola-Kirchhoff stress
tensor.

When the balance of angular momentum of the entire system, deduced by
neglecting the interaction terms, is compared to the form of the
equation obtained by summation of the individual balance of angular
momenta (\ref{eu-globalbalanceofangularmomentum}) over all species
present, one obtains a relationship between the interaction forces and
interconversion terms that is identical to Equation~(\ref{eu%
  -summationrelationmomentum}).

\subsection{Balance of energy}
\label{eu-balance-of-energy}

Recall from Section~\ref{balance-of-energy} that the internal energy
per unit mass of species~$\iota$ is denoted $e^\iota$, the external
heat supply to species~$\iota$ per unit mass of that species is
$r^\iota$ and the interaction energy, $\tilde{e}^\iota$, accounts for
the energy transferred to $\iota$ by all other species, also per unit
mass of $\iota$. We now denote the partial heat flux vector of $\iota$
defined on $\Omega_{t}$ as $\bh^\iota$, and for notational simplicity,
introduce the total energy of each species~$\iota$ per unit mass:
$E^{\iota} = e^\iota + \frac{1}{2} \Vert\bv^\iota\Vert^2$.

With these quantities defined, the rate of change of internal and
kinetic energies of species~$\iota$ under the action of mechanical
loads, processes of mass production and transport, and heating and
energy transfer, in $\Omega_{t}$, is,

\begin{equation}
\begin{split}
\underbrace{\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Omega_{t}}
  \rho^{\iota} E^{\iota}\ \mathrm{d}v \right)}_{\text{Rate of change
    of energy}} = & \underbrace{\int_{\Omega_{t}} \rho^{\iota}
  \left(\bg^{\iota} +
  \bq^{\iota}\right)\cdot\bv^{\iota}\ \mathrm{d}v}_{\text{Work done by
    body forces}}\quad + \underbrace{\int_{\partial \Omega_{t}}
  \left(\Bsigma^{\iota}\bn\right)
  \cdot\bv^{\iota}\ \mathrm{d}a}_{\text{Work done by boundary
    traction}}\\ + & \underbrace{\int_{\Omega_{t}}
  \pi^{\iota}E^{\iota}\ \mathrm{d}v}_{\text{Energy being created}} -
\underbrace{\int_{\partial \Omega_{t}} \left(\rho^{\iota} E^{\iota}
  \right) \left(\bv^{\iota} - \bv\right) \cdot
  \bn\ \mathrm{d}a}_{\text{Energy lost due to mass flux}} \\ +
&\underbrace{\int_{\Omega_{t}} \rho^{\iota} \left(r^{\iota} +
  \tilde{e}^{\iota} \right)\ \mathrm{d}v}_{\text{Energy supplied}} -
\underbrace{\int_{\partial \Omega_{t}}
  \bh^{\iota}\cdot\bn\ \mathrm{d}a.}_{\text{Heat outflux}}
\end{split}
\label{eu-globalbalanceofenergy}
\end{equation}

On applying Reynolds' Transport Theorem (Appendix~\ref{reynolds%
  -transport}) to the left hand-side and Gauss' Divergence Theorem
(Appendix~\ref{gauss-divergence}) to the area integrals, we have,

\begin{equation*}
\begin{split}
\int_{\Omega_{t}} \frac{\partial \left(\rho^{\iota}
  E^{\iota}\right)}{\partial t} \ \mathrm{d}v = & \int_{\Omega_{t}}
\rho^{\iota} \left(\bg^{\iota} +
\bq^{\iota}\right)\cdot\bv^{\iota}\ \mathrm{d}v + \int_{\Omega_{t}}
\left(\mathrm{div}\left(\Bsigma^{\iota}\right)\cdot\bv^{\iota} +
\Bsigma^{\iota}\colon\mathrm{grad}\left(\bv^{\iota}\right)
\right)\ \mathrm{d}v\\ + & \int_{\Omega_{t}}
\pi^{\iota}E^{\iota}\ \mathrm{d}v - \int_{\Omega_{t}}
\mathrm{div}\left(E^{\iota} \rho^{\iota}
\bv^{\iota}\right)\ \mathrm{d}v \\ + &\int_{\Omega_{t}} \rho^{\iota}
\left(r^{\iota} + \tilde{e}^{\iota} \right)\ \mathrm{d}v -
\int_{\Omega_{t}} \mathrm{div}\left(\bh^{\iota}\right)\ \mathrm{d}v.
\end{split}
\end{equation*}

\noindent Upon further simplification with the product rule and the
balance of mass (\ref{eu-localbalanceofmass}), we are left with:

\begin{equation*}
\begin{split}
\int_{\Omega_{t}} \rho^{\iota} \frac{\partial E^{\iota}}{\partial
  t}\ \mathrm{d}v = & \int_{\Omega_{t}} \rho^{\iota} \left(\bg^{\iota}
+ \bq^{\iota}\right)\cdot\bv^{\iota}\ \mathrm{d}v + \int_{\Omega_{t}}
\left(\mathrm{div}\left(\Bsigma^{\iota}\right)\cdot\bv^{\iota} +
\Bsigma^{\iota}\colon\mathrm{grad}\left(\bv^{\iota}\right)
\right)\ \mathrm{d}v\\ - & \int_{\Omega_{t}} \rho^{\iota}\bv^{\iota}
\cdot \mathrm{grad}\left( E^{\iota} \right)\ \mathrm{d}v \\ +
&\int_{\Omega_{t}} \rho^{\iota} \left(r^{\iota} + \tilde{e}^{\iota}
\right)\ \mathrm{d}v - \int_{\Omega_{t}}
\mathrm{div}\left(\bh^{\iota}\right)\ \mathrm{d}v.
\end{split}
\end{equation*}

We now expand the total energy of each species~$\iota$, $E^{\iota} =
e^\iota + \frac{1}{2} \Vert\bv^\iota\Vert^2$, and take its spatial and
temporal derivatives to give,

\begin{equation*}
\begin{split}
\int_{\Omega_{t}} \rho^{\iota} \left(\frac{\partial
  e^{\iota}}{\partial t} + \bv^{\iota}\cdot\frac{\partial
  \bv^{\iota}}{\partial t}\right)\ \mathrm{d}v = & \int_{\Omega_{t}}
\rho^{\iota} \left(\bg^{\iota} +
\bq^{\iota}\right)\cdot\bv^{\iota}\ \mathrm{d}v\\ + &
\int_{\Omega_{t}}
\left(\mathrm{div}\left(\Bsigma^{\iota}\right)\cdot\bv^{\iota} +
\Bsigma^{\iota}\colon\mathrm{grad}\left(\bv^{\iota}\right)
\right)\ \mathrm{d}v\\ - & \int_{\Omega_{t}} \rho^{\iota}\bv^{\iota}
\cdot \left(\mathrm{grad}\left(e^{\iota}\right)
+\mathrm{grad}(\bv^{\iota})^{\mathrm{T}}
\bv^{\iota}\right)\ \mathrm{d}v \\ + &\int_{\Omega_{t}} \rho^{\iota}
\left(r^{\iota} + \tilde{e}^{\iota} \right)\ \mathrm{d}v -
\int_{\Omega_{t}} \mathrm{div}\left(\bh^{\iota}\right)\ \mathrm{d}v.
\end{split}
\end{equation*}

\noindent Finally, we apply the balance of linear momentum~%
(\ref{eu-localbalanceofmomentum}) to the equation above and, upon
localisation, arrive at:

\begin{equation}
\rho^{\iota} \frac{\partial e^{\iota}}{\partial t} =
\Bsigma^{\iota}\colon\mathrm{grad}\left(\bv^{\iota}\right) +
\rho^{\iota} \left(r^{\iota} + \tilde{e}^{\iota} \right) -
\mathrm{div}\left(\bh^{\iota}\right) - \rho^{\iota}\bv^{\iota} \cdot
\mathrm{grad} \left(e^{\iota}\right) \quad
\mathrm{in}\ \Omega_{t},
\label{eu-localbalanceofenergy}
\end{equation}

\noindent the form of the balance of energy of a species~$\iota$ in
$\Omega_{t}$ which is most convenient for combining with the entropy
inequality, leading to the Clausius-Duhem form of the dissipation
inequality~(\ref{eu-clausiusduhemform}) in the following section.
This result is consistent with classical mixture theory \citep{%
  TruesdellToupin:60} and is the current configuration analogue of
Equation~(\ref{localbalanceofenergy}).

In order to obtain a relationship between the momentum transfer, mass
interconversion and inter-species energy transfer terms within the
system, we first express the rate of change of energy of the system
interacting with its environment from the point of view of an external
observer unaware of these internal interactions, as:

\begin{equation}
\begin{split}
\sum\limits_{\iota}
\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Omega_{t}} \rho^{\iota}
E^{\iota}\ \mathrm{d}v \right) = & \sum\limits_{\iota}
\int_{\Omega_{t}} \rho^{\iota} \bg^i\cdot\bv^{\iota}\ \mathrm{d}v +
\sum\limits_{\iota} \int_{\partial \Omega_{t}}
\left(\Bsigma^{\iota}\bn\right)\cdot\bv^{\iota}\ \mathrm{d}a\\ -&
\sum\limits_{\iota} \int_{\partial \Omega_{t}} \left(\rho^{\iota}
E^{\iota} \right) \left(\bv^{\iota} - \bv\right) \cdot
\bn\ \mathrm{d}a \\ + & \sum\limits_{\iota} \int_{\Omega_{t}}
\rho^{\iota} r^{\iota}\ \mathrm{d}v - \sum\limits_{\iota}
\int_{\partial \Omega_{t}} \bh^{\iota}\cdot\bn\ \mathrm{d}a.
\end{split}
\label{eu-energyexternal}
\end{equation}
 
\noindent Since the equation above, and a summation of Equation~%
(\ref{eu-globalbalanceofenergy}) over all $\iota$ are equivalent
statements of the balance of energy for the entire system, it follows
upon inspection and localisation that,

\begin{equation}
 \sum\limits_{\iota} \left( \rho^{\iota}\bq^{\iota}\cdot\bv^{\iota} +
 \pi^{\iota} E^{\iota} + \rho^{\iota}\tilde{e}^{\iota} \right) = 0,
\label{eu-summationrelationenergy}
\end{equation}

\noindent which ensures that there is no net energy production
mechanism internal to the system.

\section{Thermodynamically-consistent constitutive framework}
\label{eu-entropy-inequality}

In order to close the system of balance equations deduced in
Sections~\ref{eu-balance-of-mass}--\ref{eu-balance-of-energy}, and
make thermodynamically-valid constitutive choices pertinent to
biological tissues, we turn to the principles of entropy growth and
material frame-indifference. The following treatment builds upon the
same fundamental assumptions underlying our system of interest as
those stated in Section~\ref{entropy-inequality}:

\begin{enumerate}
\item[(\romannumeral 1)] The entropy production inequality is assumed
  to hold at a continuum point for all species as a whole, but, in
  general, not for each individual species.\footnote{Enforcing the
    entropy inequality separately for each of the individual
    constituents of the mixture imposes unrealistic constraints on the
    mixture \citep{BedfordDrumheller:1983}.}
\item[(\romannumeral 2)] All species occupying a continuum point in
  the tissue have the same absolute temperature, $\theta$.
\end{enumerate} 

\noindent Additionally, in the derivation of the Clausius-Duhem form
of the Second Law of Thermodynamics in
Section~\ref{eu-clausius-duhem-form} immediately below, we assume the
viewpoint of an observer external to the system unaware of any
internal interactions.

\subsection{The Clausius-Duhem form}
\label{eu-clausius-duhem-form}

With the assumptions introduced above, and denoting by $\eta^\iota$
the entropy per unit mass of species~$\iota$, the entropy inequality,
when written out for the entire system in $\Omega_{t}$ reads,

\begin{equation}
\begin{split}
\sum\limits_{\iota}
\underbrace{\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{\Omega_{t}}
  \rho^{\iota} \eta^{\iota}\ \mathrm{d}v \right)}_{\text{Rate of
    change of entropy}} \geq & \sum\limits_{\iota}
\underbrace{\int_{\Omega_{t}} \frac{ \rho^{\iota}
    r^{\iota}}{\theta}\ \mathrm{d}v}_{\text{Entropy supplied}} -
\sum\limits_{\iota} \underbrace{\int_{\partial \Omega_{t}}
  \frac{\bh^{\iota}\cdot\bn}{\theta}\ \mathrm{d}a.}_{\text{Entropy
    outflux}} \\ - & \sum\limits_{\iota} \underbrace{\int_{\partial
    \Omega_{t}} \left(\rho^{\iota} \eta^{\iota} \right)
  \left(\bv^{\iota} - \bv\right) \cdot
  \bn\ \mathrm{d}a.}_{\text{Entropy lost due to mass flux}}
\end{split}
\label{eu-globalbalanceofentropy}
\end{equation}

On applying Reynolds' Transport Theorem (Appendix~\ref{reynolds%
  -transport}) to the left hand-side and Gauss' Divergence Theorem
(Appendix~\ref{gauss-divergence}) to the area integrals, we have,

\begin{equation*}
\begin{split}
\sum\limits_{\iota} \int_{\Omega_{t}} \frac{\partial
  \left(\rho^{\iota} \eta^{\iota}\right) }{\partial t} \ \mathrm{d}v
\geq & \sum\limits_{\iota} \int_{\Omega_{t}} \frac{ \rho^{\iota}
  r^{\iota}}{\theta}\ \mathrm{d}v - \sum\limits_{\iota}
\int_{\Omega_{t}} \left(
\frac{\mathrm{div}\left(\bh^{\iota}\right)}{\theta} -
\frac{\bh^{\iota}\cdot\mathrm{grad}\left(\theta\right)}{\theta^{2}}
\right)\ \mathrm{d}v \\ - & \sum\limits_{\iota} \int_{\Omega_{t}}
\left( \eta^{\iota}\ \mathrm{div} \left(\rho^{\iota}\bv^{\iota}
\right) +
\rho^{\iota}\bv^{\iota}\cdot\mathrm{grad}\left(\eta^{\iota}\right)
\right)\ \mathrm{d}v.
\end{split}
\end{equation*}

\noindent Applying the product rule and the balance of mass%
~(\ref{eu-localbalanceofmass}),

\begin{equation*}
\begin{split}
\sum\limits_{\iota} \int_{\Omega_{t}} ( \rho^{\iota} \frac{\partial
  \eta^{\iota}}{\partial t} +
\cancelto{\pi^{\iota}\eta^{\iota}}{\eta^{\iota} \frac{\partial
    \rho^{\iota}}{\partial t}}) \ \mathrm{d}v \geq &
\sum\limits_{\iota} \int_{\Omega_{t}} \frac{ \rho^{\iota}
  r^{\iota}}{\theta}\ \mathrm{d}v\\ -& \sum\limits_{\iota}
\int_{\Omega_{t}} \left(
\frac{\mathrm{div}\left(\bh^{\iota}\right)}{\theta} -
\frac{\bh^{\iota}\cdot\mathrm{grad}\left(\theta\right)}{\theta^{2}}
\right)\ \mathrm{d}v \\ - & \sum\limits_{\iota} \int_{\Omega_{t}}
\left( \cancel{\eta^{\iota}\ \mathrm{div}
  \left(\rho^{\iota}\bv^{\iota} \right)} +
\rho^{\iota}\bv^{\iota}\cdot\mathrm{grad}\left(\eta^{\iota}\right)
\right)\ \mathrm{d}v,
\end{split}
\end{equation*}

\noindent and rearranging terms and localising, we have the following
form of the entropy inequality for the entire system:

\begin{equation}
\begin{split}
\sum\limits_{\iota}\ \left( \rho^{\iota} \frac{\partial
  \eta^{\iota}}{\partial t} + \pi^{\iota}\eta^{\iota}\right) \geq &
\sum\limits_{\iota} \left( \frac{ \rho^{\iota} r^{\iota}}{\theta} -
\rho^{\iota}\bv^{\iota}\cdot\mathrm{grad}\left(\eta^{\iota}
\right)\right)\ \\ -& \sum\limits_{\iota} \left(
\frac{\mathrm{div}\left(\bh^{\iota}\right)}{\theta} -
\frac{\bh^{\iota}\cdot\mathrm{grad}\left(\theta\right)}{\theta^{2}}
\right).
\end{split}
\label{eu-localentropyinequality}
\end{equation}

Now, multiplying Equation~(\ref{eu-localentropyinequality}) by the
temperature field, $\theta$, and subtracting it from a summation of
the balance of energy~(\ref{eu-localbalanceofenergy}) over all
species~$\iota$, we obtain,

\begin{equation*}
\begin{split}
\sum\limits_{\iota}\ \left(\rho^{\iota} \dot{e^{\iota}} - \rho^{\iota}
\dot{\eta^{\iota}} \theta - \pi^{\iota}\eta^{\iota}\theta \right) \leq
& \sum\limits_{\iota} \left(
\Bsigma^{\iota}\colon\mathrm{grad}\left(\bv^{\iota}\right) -
\frac{\bh^{\iota}\cdot\mathrm{grad}\left(\theta\right)}{\theta}
\right)\ \\ +& \sum\limits_{\iota}\bigg( \rho^{\iota}
\tilde{e}^{\iota} - \rho^{\iota}\bv^{\iota} \cdot\left( \mathrm{grad}
\left(e^{\iota}\right) - \mathrm{grad}\left(\eta^{\iota}
\right)\theta\right)\bigg),
\end{split}
\end{equation*}

\noindent where $\dot{\bullet}$ denotes the partial derivative with
respect to time. Using Equation~(\ref{eu-summationrelationenergy}) to
eliminate the species interaction energy, $\tilde{e}^{\iota}$, we
deduce that

\begin{equation*}
\begin{split}
\sum\limits_{\iota}\ \left(\rho^{\iota} \dot{e^{\iota}} - \rho^{\iota}
\dot{\eta^{\iota}} \theta - \pi^{\iota}\eta^{\iota}\theta \right) \leq
& \sum\limits_{\iota} \left(
\Bsigma^{\iota}\colon\mathrm{grad}\left(\bv^{\iota}\right) -
\frac{\bh^{\iota}\cdot\mathrm{grad}\left(\theta\right)}{\theta}
\right)\ \\ -& \sum\limits_{\iota} \rho^{\iota}\bv^{\iota} \cdot\left(
\mathrm{grad} \left(e^{\iota}\right) - \mathrm{grad}\left(\eta^{\iota}
\right)\theta\right) \\ -& \sum\limits_{\iota} \left(
\rho^{\iota}\bq^{\iota}\cdot\bv^{\iota} + \pi^{\iota} \left(e^\iota +
\frac{1}{2} \Vert\bv^\iota\Vert^2\right) \right).
\end{split}
\end{equation*}

Introducing the Helmholtz free energy per unit mass of
species~$\iota$, \mbox{$\psi^\iota = e^\iota - \theta\eta^\iota$}, and
regrouping terms, the relation above reduces to,

\begin{equation}
\begin{split}
\sum\limits_{\iota}\ \bigg(\rho^{\iota} \dot{e^{\iota}} - \rho^{\iota}
\dot{\eta^{\iota}} \theta
-\Bsigma^{\iota}\colon\mathrm{grad}\left(\bv^{\iota}\right) +
\frac{\bh^{\iota}\cdot\mathrm{grad}\left(\theta\right)}{\theta} \bigg)
& \\ + \sum\limits_{\iota} \left( \rho^{\iota}
(\bq^{\iota}+\mathrm{grad} \left(e^{\iota}\right) -
\mathrm{grad}\left(\eta^{\iota} \right)\theta) \cdot\bv^{\iota} +
\pi^{\iota} \left(\psi^{\iota}+\frac{1}{2}
\Vert\bv^\iota\Vert^2\right) \right) & \leq 0.
\end{split}
\label{eu-clausiusduhemform}
\end{equation}

\noindent the Clausius-Duhem (or reduced dissipation) inequality for
the growth process; a rule which any prescribed constitutive
relationship must satisfy \citep{TruesdellToupin:60}.

The form of the Clausius-Duhem inequality arrived at in (\ref{eu%
  -clausiusduhemform}) is equivalent to the forms in recent work on
mixture theory-based models for biological growth \citep{loret05,
  ateshian07}. However, subsequently varying choices made in the
different works, including this one, for the constitutive independent
variables result in altered constitutive specification. I believe it
is significant, and must be emphasised at this point that not only do
the constitutive choices detailed in the following sections ensure
that the Clausius-Duhem inequality is satisfied {\em a priori}, but
also that they are general enough to handle a fairly large class of
physics, and most significantly, have been implemented in a coupled
formulation retaining much of their rich detail, as evidenced by the
computational examples in Chapter~\ref{numerical-simulations-2}.

\subsection{Duhamel's law of heat conduction}
\label{eu-duhamel-law}

A suitable constitutive relation motivated by experiment which relates
the Cauchy heat flux vector $\bh^{\iota}(\bx,t)$ to the spatial
temperature gradient $\mathrm{grad}\left(\theta\right)(\bx,t)$ is the
classical heat conduction law of Duhamel:

\begin{equation}
\bh^{\iota} = -
\Bkappa^{\iota}\ \mathrm{grad}\left(\theta\right),\;\forall\,\iota,
\label{eu-heat-conduction}
\end{equation}

\noindent where $\Bkappa^{\iota}(\bx,t)$, the spatial thermal
conductivity tensor, is positive semi-definite. The relation above
states that heat flows down a temperature gradient, and using it
ensures that the following portion of the Clausius-Duhem inequality
(\ref{eu-clausiusduhemform}),

\begin{equation*}
\sum\limits_{\iota} \left(\frac{\bh^{\iota} \cdot
  \mathrm{grad}\left(\theta\right)}{\theta} \right) \leq 0,
\end{equation*}

\noindent is satisfied since the heat conduction law (\ref{eu-heat%
  -conduction}) is valid for all species in the mixture.

In fact, the constitutive form~(\ref{eu-heat-conduction}) specified
above was motivated directly by virtue of it being sufficient to
ensure that the corresponding portion of the Clausius-Duhem inequality
had the appropriate sign. Another approach that can be followed to
arrive at Duhamel's law (and, similarly, the constitutive relations
that follow) is to first define a general set of independent variables
upon which the various field variables can depend (such as in
\citet{TruesdellToupin:60}). One then proceeds to consider a scenario
where the functional forms of some of the terms in the dissipation
inequality are fixed, and examine the implications of entropy
production on the remaining terms, assuming they can be varied
independently.

In practice, the two approaches outlined above are mathematically
equivalent, and the former is employed in the following treatment.

While Duhamel's law (\ref{eu-heat-conduction}) is indeed useful in a
general setting, our primary interest lies in specialising concepts
from classical thermodynamics to biological systems, where the
temperature field under most physiological scenarios is relatively
uniform and constant in time. In what follows, we will restrict our
attention to this case. Additionally, the general nature of the
analysis thus far has resulted in it being presented in terms of an
arbitrary species~$\iota$. Recall, from Section~\ref{balance-laws},
that the various species constituting the tissue are grouped into a
solid species, consisting of solid \emph{collagen fibrils} and
\emph{cells}, denoted by $\mathrm{c}$, an extra-cellular \emph{fluid}
species, denoted by $\mathrm{f}$, consisting primarily of water, and
\emph{solute} species, consisting of precursors to reactions,
byproducts, nutrients, and other regulatory chemicals, denoted by
$\mathrm{s}$. We now return to this earlier classification, and
motivated by specific portions of (\ref{eu-clausiusduhemform}),
proceed to prescribe constitutive relations pertinent to the mechanics
and biochemistry of growing biological tissue for each of these
species groups.

\subsection{Energy-dependent mass source terms}
\label{eu-energy-dependent-source}

We first focus on the following term of the Clausius-Duhem form
(\ref{eu-clausiusduhemform}),

\begin{equation*}
\sum\limits_{\iota=\mathrm{c},\mathrm{f},\mathrm{s}} \pi^{\iota}
\left(\psi^{\iota}+\frac{1}{2} \Vert\bv^\iota\Vert^2\right) \leq 0.
\end{equation*}

Recall that the extra-cellular fluid species, $\mathrm{f}$, consists
primarily of water and does not take part in reactions, i.e.,
$\pi^\mathrm{f} = 0$. Therefore, from the summation relation for mass
sources and sinks~(\ref{eu-summationrelationmass}), we have that
$\pi^\mathrm{s} = - \pi^\mathrm{c}$. Using this fact, the inequality
above reads,

\begin{equation*}
 \pi^{\mathrm{c}} \left(\psi^{\mathrm{c}}+\frac{1}{2}
 \Vert\bv^{\mathrm{c}}\Vert^2\right) + \left(-\pi^{\mathrm{c}}\right)
 \left(\psi^{\mathrm{s}}+\frac{1}{2}
 \Vert\bv^{\mathrm{s}}\Vert^2\right) \leq 0.
\end{equation*}

\noindent The above relation provides a broad guideline that allows
for the specification of a very general class of source terms for
collagen. An obvious form that it suggests is linear in terms of the
energy (Helmholtz free energy plus the kinetic energy) difference
between the collagen phase and the solute phase:

\begin{equation}
 \pi^{\mathrm{c}} = -\kappa^{\mathrm{c}}
 \left(\psi^{\mathrm{c}}+\frac{1}{2} \Vert\bv^{\mathrm{c}}\Vert^2 -
 \psi^{\mathrm{s}}-\frac{1}{2} \Vert\bv^{\mathrm{s}}\Vert^2\right),
\label{eu-energy-source-linear}
\end{equation}

\noindent where $\kappa^{\mathrm{c}}$ is a non-negative rate constant.

Thus, the thermodynamics suggests that it is an energetic difference
between the reactants and products of a chemical reaction that drive
the reaction forward; a well-established concept in chemistry. More
interestingly, since all that the thermodynamics requires is that the
source term for collagen be positive (negative) when the solute is
more (less) energetic than collagen, this opens up the selection to a
variety of forms, such as the nonlinear (in terms of the energy
difference) example shown below,

\begin{equation}
\pi^{\mathrm{c}}=\epsilon\ \kappa^{\mathrm{c}}\left( \mathrm{exp}
\left[-\epsilon\ U \left(\psi^{\mathrm{c}}+\frac{1}{2}
  \Vert\bv^{\mathrm{c}}\Vert^2 - \psi^{\mathrm{s}} - \frac{1}{2}
  \Vert\bv^{\mathrm{s}}\Vert^2\right) \right] - 1\right)
\label{eu-energy-source-nonlinear}
\end{equation}

\noindent (where $\epsilon =
\mathrm{sign}\left(\psi^{\mathrm{c}}+\frac{1}{2}
\Vert\bv^{\mathrm{c}}\Vert^2 - \psi^{\mathrm{s}}-\frac{1}{2}
\Vert\bv^{\mathrm{s}}\Vert^2\right)$, $U>0$ and
$\kappa^{\mathrm{c}}\geq 0$), and can be tailored to represent
different sorts of biochemistry. The evolution of this nonlinear
collagen source term with the energy difference between the collagen
and solute phases is shown in Figure~\ref{energy-based-source}; which
clearly shows the source for collagen having a sign opposite to that
of the energy difference between collagen and the solute.

\begin{figure}
  \centering
  \includegraphics[width=0.6\textwidth,angle=270]{images/elucidation/%
    energy-based-source}
  \caption{A thermodynamically-motivated collagen source.}
  \label{energy-based-source}
\end{figure}

With the introduction of an energy-dependent mass source (such as
(\ref{eu-energy-source-linear}) or (\ref{eu-energy-source%
  -nonlinear})), and drawing upon the assumption that the temperature
field is uniform and constant\footnote{Which allows one to write
  $\dot{e^{\iota}} - \dot{\eta^{\iota}} \theta = \dot{\psi^{\iota}}$
  and $\mathrm{grad} \left(e^{\iota}\right) -
  \mathrm{grad}\left(\eta^{\iota} \right)\theta =
  \mathrm{grad}\left(\psi^{\iota} \right)$.} during the course of a
biological experiment, a portion of the Clausius-Duhem form (\ref{eu%
  -clausiusduhemform}) is satisfied a priori. The remainder now reads:

\begin{equation}
\sum\limits_{\iota=\mathrm{c},\mathrm{f},\mathrm{s}}\ \bigg(\rho^{\iota}
\dot{\psi^{\iota}}
-\Bsigma^{\iota}\colon\mathrm{grad}\left(\bv^{\iota}\right)+
\rho^{\iota} (\bq^{\iota}+\mathrm{grad}\left(\psi^{\iota}\right))
\cdot\bv^{\iota} \bigg) \leq 0.
\label{eu-reduced-dissipation-1}
\end{equation}

In the treatment that follows, recognising that the fluid, solid and
solute are fundamentally different substances, we turn to the
specification of different constitutive independent variables most
suited to each of them in deriving relevant constitutive
relationships. Since our primary interest lies in the mechanics
response of soft collagenous tissues, we will be paying particular
attention to the constitutive forms of the solid collagen and fluid
stress tensors, as well as the interaction forces between the collagen
and fluid phases.

It is well established that collagenous tissues demonstrate a
rate-dependent response to applied loads (see, for e.g.,
\citet{Provenzanoetal:2001}). This behaviour may be attributed to the
viscoelastic deformation of the collagen network (which is the focus
of Section~\ref{eu-viscoelastic-solid}), the viscosity of the
extra-cellular fluid (as seen in Section~\ref{eu-newtonian-fluid}),
frictional effects arising from stress-driven fluid flow through the
network (as discussed in Section~\ref{eu-interaction-forces}), or
combinations thereof. In order to be able to model each of these
cases, the following sections derive the relevant constitutive
relations entirely from thermodynamic considerations.

\subsection{A viscoelastic solid}
\label{eu-viscoelastic-solid}

For the terms arising from the solid collagen phase, the reduced
dissipation inequality (\ref{eu-reduced-dissipation-1}) requires that,

\begin{equation*}
\rho^{\mathrm{c}} \dot{\psi^{\mathrm{c}}}
-\Bsigma^{\mathrm{c}}\colon\mathrm{grad}\left(\bv^{\mathrm{c}}\right)+
\rho^{\mathrm{c}} \mathrm{grad}\left(\psi^{\mathrm{c}}\right)
\cdot\bv^{\mathrm{c}} +\rho^{\mathrm{c}}
\bq^{\mathrm{c}}\cdot\bv^{\mathrm{c}} \leq 0.
\end{equation*}

\noindent We will return to the last term, $\rho^{\mathrm{c}}
\bq^{\mathrm{c}}\cdot\bv^{\mathrm{c}}$ when discussing interaction
forces between the solid and fluid phases in Section~(\ref{eu%
  -interaction-forces}). Currently, we turn our attention to the
remaining terms in the inequality above:

\begin{equation}
\rho^{\mathrm{c}} \dot{\psi^{\mathrm{c}}}
-\Bsigma^{\mathrm{c}}\colon\mathrm{grad} \left(\bv^{\mathrm{c}}\right)
+\rho^{\mathrm{c}} \mathrm{grad}\left(\psi^{\mathrm{c}}\right)
\cdot\bv^{\mathrm{c}} \leq 0.
\label{eu-reduced-dissipation-solid}
\end{equation}

Looking ahead to the numerical examples presented in Chapter~\ref{%
  numerical-simulations-2}, and noting that they solve the momentum
balance equation (\ref{localbalanceofmomentum}) for the solid phase in
the reference configuration,\footnote{It is instructive to note here
  that since $\Omega_{t}$ is constructed to coincide with the current
  position of the solid component of the tissue (see
  Section~\ref{eu-balance-laws}), we are aware of the deformation
  gradient $\bF^{\mathrm{c}}=\bF=\frac{\partial \Bvarphi_{t}}{\partial
    \bX}$ of this phase. This allows us to consistently pose the
  balance laws in either configuration.} we deduce the constitutive
response of the solid phase in terms of its partial first and second
Piola-Kirchhoff stress tensors, $\bP^{\mathrm{c}}$ and
$\bS^{\mathrm{c}}$, respectively.

From the balance of angular momentum for the solid (\ref{eu%
  -localbalanceofangularmomentum}), we know that its partial Cauchy
stress is symmetric and thus, $\Bsigma^{\mathrm{c}} \colon
\mathrm{grad} \left(\bv^{\mathrm{c}}\right) =\Bsigma^{\mathrm{c}}
\colon (\dot{\bF}\bF^{-1}) =\frac{1}{J} \bP^{\mathrm{c}} \colon
\dot{\bF}$. Substituting this in (\ref{eu-reduced-dissipation-solid})
above, recalling that \mbox{$\rho^{\mathrm{c}}\ J =
  \rho_0^{\mathrm{c}}$} and $\rho^{\mathrm{c}}
\dot{\psi^{\mathrm{c}}}+ \rho^{\mathrm{c}}
\mathrm{grad}\left(\psi^{\mathrm{c}}\right) \cdot\bv^{\mathrm{c}}$
provides the material time derivative of $\psi^{\mathrm{c}}$, we have,

\begin{equation*}
\rho_{0}^{\mathrm{c}} \dot{\psi^{\mathrm{c}}}
-\bP^{\mathrm{c}}\colon\dot{\bF}\leq 0,\ \mathrm{in}\ \Omega_{0}.
\end{equation*}

\noindent Introducing the elasto-growth kinematic split
$\bF=\bF^{\mathrm{e}}\bF^{\mathrm{g}}$ discussed in Section%
~\ref{kinematics-of-growth} for the solid phase, c, and applying the
properties of the contraction operation, the above inequality reads,

\begin{equation*}
\rho_{0}^{\mathrm{c}} \dot{\psi^{\mathrm{c}}} - \bP^{\mathrm{c}}
\bF^{\mathrm{g}^{\mathrm{T}}} \colon \dot{\bF^{\mathrm{e}}} -
\bF^{\mathrm{e}^{\mathrm{T}}} \bP^{\mathrm{c}} \colon
\dot{\bF^{\mathrm{g}}} \leq 0.
\end{equation*}

\noindent The term $-\bF^{\mathrm{e}^{\mathrm{T}}} \bP^{\mathrm{c}}
\colon \dot{\bF^{\mathrm{g}}}$ is the focus of the following section
on stress-dependent growth deformation gradients
(\ref{eu-stress-dependent-growth}). Currently, we consider only the
following two terms,

\begin{equation}
\rho_{0}^{\mathrm{c}} \dot{\psi^{\mathrm{c}}} - \bP^{\mathrm{c}}
\bF^{\mathrm{g}^{\mathrm{T}}} \colon \dot{\bF^{\mathrm{e}}} \leq 0.
\label{eu-reduced-dissipation-solid-2}
\end{equation}

As in Section~\ref{constitutive-framework}, if we were to assume at
this point that the Helmholtz free energy of the solid had the form,
$\psi^{\mathrm{c}} = \frac{1}{\tilde{\rho_{0}^{\mathrm{c}}}}
\hat{\psi^{\mathrm{c}}} (\bF^{\mathrm{e}})$, where
$\tilde{\rho_{0}^{\mathrm{c}}}$ is the intrinsic density of the solid
collagen in the reference configuration, i.e., the free energy depends
only upon the elastic portion of the deformation gradient, then a
sufficient condition to satisfy the inequality above is to assume a
hyperelastic model of the form:

\begin{equation}
\bP^{\mathrm{c}} =
\frac{\rho_{0}^{\mathrm{c}}}{\tilde{\rho_{0}^{\mathrm{c}}}}
\frac{\partial \hat{\psi^{\mathrm{c}}}}{\partial \bF^{\mathrm{e}}}
\bF^{\mathrm{g}^{\mathrm{-T}}}.
\label{eu-solid-hyperelastic}
\end{equation}

\noindent We have already seen one specific form of this model used in
the computations in Section~\ref{constitutive-framework}.

However, we are currently interested in allowing for an inelastic
response in the solid, and for this, we turn to a body of established
work on continuum formulations for viscoelastic materials undergoing
finite strains (see, for e.g., \citet{simo86}, \citet{holzapfel96},
and \citet{SimoHughes:98}). The treatment below follows in the same
vein.

We begin by assuming a Helmholtz free energy for the solid collagen of
the form: $\psi^{\mathrm{c}} = \frac{1}{\tilde{\rho_{0}^{\mathrm{c}}}}
\hat{\psi^{\mathrm{c}}} (\bC^{\mathrm{e}}, \BGamma_{1}, \ldots,
\BGamma_{m})$, where $\bC^{\mathrm{e}}= \bF^{\mathrm{e}^{\mathrm{T}}}
\bF^{\mathrm{e}}$ is the elastic right Cauchy-Green tensor and
$\BGamma_{1}, \ldots, \BGamma_{m}$ are a set of second order tensorial
internal history variables. Substituting this form into
(\ref{eu-reduced-dissipation-solid-2}), and rewriting the partial
first Piola-Kirchhoff stress tensor in terms of the partial second
Piola-Kirchhoff stress tensor using the relation $\bP^{\mathrm{c}} =
\bF \bS^{\mathrm{c}}$, we have,\footnote{The computations presented in
  Chapter~\ref{numerical-simulations-2} use the second partial
  Piola-Kirchhoff stress tensor,~$\bS^{\mathrm{c}}$, since it is a
  symmetric quantity, and thus requires less memory for storage.}

\begin{equation}
\frac{\rho_{0}^{\mathrm{c}}}{\tilde{\rho_{0}^{\mathrm{c}}}}
\frac{\partial \hat{\psi^{\mathrm{c}}}}{\partial \bC^{\mathrm{e}}}
\colon (2 \bF^{\mathrm{e}^{\mathrm{T}}} \dot{\bF^{\mathrm{e}}}) +
\sum_{\alpha=1}^{m}
\frac{\rho_{0}^{\mathrm{c}}}{\tilde{\rho_{0}^{\mathrm{c}}}}
\frac{\partial \hat{\psi^{\mathrm{c}}}}{\partial \BGamma_{\alpha}}
\colon \dot{\BGamma}_{\alpha} - (\bF^{\mathrm{e}}\bF^{\mathrm{g}})
\bS^{\mathrm{c}} \bF^{\mathrm{g}^{\mathrm{T}}} \colon
\dot{\bF^{\mathrm{e}}} \leq 0.
\label{eu-reduced-dissipation-solid-3}
\end{equation}

\noindent A sufficient condition to satisfy (\ref{eu-reduced%
  -dissipation-solid-3}) is to specify that the partial second
Piola-Kirchhoff stress tensor has the form,

\begin{equation}
\bS^{\mathrm{c}} = \bF^{\mathrm{g}^{\mathrm{-1}}}
\frac{\rho_{0}^{\mathrm{c}}}{\tilde{\rho_{0}^{\mathrm{c}}}} 2
\frac{\partial \hat{\psi^{\mathrm{c}}}}{\partial \bC^{\mathrm{e}}}
\bF^{\mathrm{g}^{\mathrm{-T}}},
\label{eu-solid-viscoelastic}
\end{equation}

\noindent and provide a suitable evolution equation for the internal
variables, $\BGamma_{\alpha}$. Motivated by the fact that some
compressible materials exhibit dissimilar bulk and shear response, we
proceed to decompose the free energy function into volumetric and
isochoric parts:

\begin{equation}
\hat{\psi^{\mathrm{c}}} (\bC^{\mathrm{e}}, \BGamma_{1}, \ldots,
\BGamma_{m}) = \sW_{\mathrm{vol}}(J^{\mathrm{e}}) +
\sW_{\mathrm{iso}}(\bar{\bC^{\mathrm{e}}}) +
\sum_{\alpha=1}^{m}\Bgamma_{\alpha}
(\bar{\bC^{\mathrm{e}}},\BGamma_{\alpha}),
\label{eu-solid-energy-split}
\end{equation}

\noindent where $J^{\mathrm{e}}$ is the determinant of the elastic
portion of the deformation gradient tensor and $\bar{\bC ^
  {\mathrm{e}}} = J ^ {\mathrm{e}^{-2/3}}\bC^{\mathrm{e}}$. The first
two terms in the decomposition above characterise the volumetric and
isochoric equilibrium response of the solid phase, and the last term
is the dissipative potential which contributes to the viscoelastic
response. The equilibrium response (that of a purely-elastic material)
is recovered during infinitely slow processes.

Using the decomposition (\ref{eu-solid-energy-split}) in the stress
constitutive relation (\ref{eu-solid-viscoelastic}), we see that the
solid stress takes the form,

\begin{equation}
\bS^{\mathrm{c}} = \bF^{\mathrm{g}^{\mathrm{-1}}}
\frac{\rho_{0}^{\mathrm{c}}}{\tilde{\rho_{0}^{\mathrm{c}}}}\left(
\bS^{\mathrm{c}}_{\mathrm{vol}}+\bS^{\mathrm{c}}_{\mathrm{iso}}+
\sum_{\alpha=1}^{m}\bQ_{\alpha} \right)
\bF^{\mathrm{g}^{\mathrm{-T}}},
\label{eu-solid-visco-stress}
\end{equation}

\noindent where, denoting by $\mathbb{I}$ the fourth order unit
tensor,

\begin{equation*}
\bS^{\mathrm{c}}_{\mathrm{vol}} = J^{\mathrm{e}} \frac{\partial
  \sW_{\mathrm{vol}}(J^{\mathrm{e}})}{\partial J^{\mathrm{e}}} \quad
\mathrm{and} \quad \bS^{\mathrm{c}}_{\mathrm{iso}} =
J^{\mathrm{e}^{-2/3}} \left(\mathbb{I} - \frac{1}{3}\ \bC^{-1} \otimes
\bC \right) \colon 2
\frac{\partial\sW_{\mathrm{iso}}(\bar{\bC^{\mathrm{e}}})}{\partial
  \bar{\bC^{\mathrm{e}}}}.
\end{equation*}

\noindent Equation (\ref{eu-solid-visco-stress}), along with the
following evolution equations for the {\em non-equilibrium stresses},
$\bQ_{\alpha}$,\footnote{Which are work-conjugate variables to
  $\BGamma_{\alpha}$.} in agreement with the dissipation inequality
motivated by a generalised Maxwell model,

\begin{equation}
\dot{\bQ}_{\alpha}+\frac{\bQ_{\alpha}}{\tau_{\alpha}} =
\dot{\bS}^{\mathrm{c}}_{\mathrm{iso}_{\alpha}},\ \alpha=1,\ldots,m,
\label{eu-linearviscoelasticity}
\end{equation}

\noindent where each $\tau_{\alpha}$ is a characteristic relaxation
time and $\bS^{\mathrm{c}}_{\mathrm{iso}_{\alpha}} =
\beta_{\alpha}\ \bS^{\mathrm{c}}_{\mathrm{iso}}
(\bar{\bC^{\mathrm{e}}})$, where $\beta_{\alpha}$ is a non-negative
strain-energy factor associated with $\tau_{\alpha}$, completes the
specification of a linear viscoelastic solid. This can be extended in
a straightforward manner to a nonlinear model by introducing the
notion of a {\em modified relaxation time} \citep{Eyring:36}.

\subsection{Effects of the stress state on tissue growth}
\label{eu-stress-dependent-growth}

One important influence that the local stress state has on tissue
growth directly relates to its regulation of species production (and
consumption) rates. An example of this fact lies in oft-cited work of
\citet{wolff1892} who found that bone is deposited and resorbed in
accordance with the stresses placed upon it. In this formulation, this
effect is modelled through the use of the strain-energy dependent
source terms~(\ref{strain-energy-based-source}) introduced in
Section~\ref{nature-of-sources}.

There is experimental evidence to suggest that there is another
important influence of the stress state on the development of tissues,
relating to the spatial alignment of deposited matter.  An example of
this is found in \citet{Provenzanoetal:2003}, where it is observed
that during wound healing, newly-de\-po\-sit\-ed collagen fibres are
found to align with the applied stress, whereas under unstressed
conditions, they are deposited isotropically. Momentarily, will see
that the thermodynamics naturally motivates a form for the rate of
change of the growth deformation gradient tensor that reflects this
observation.

During the derivation of the constitutive relationship for the solid
stress in the preceding section, the implication of the following
inequality,

\begin{equation*}
- \bF^{\mathrm{e}^{\mathrm{T}}} \bP^{\mathrm{c}} \colon
\dot{\bF^{\mathrm{g}}} \leq 0,
\end{equation*}

\noindent had not been explored. Turning our attention to it now, it
is clear that a sufficient condition to satisfy this inequality is,

\begin{equation}
\dot{\bF^{\mathrm{g}}} = \lambda\ \bF^{\mathrm{e}^{\mathrm{T}}}
\bP^{\mathrm{c}},
\label{eu-stressgrowthtensorrate}
\end{equation}

\noindent where $\lambda$ is a non-negative scalar.

Equation~(\ref{eu-stressgrowthtensorrate}) specifies that incremental
changes in the growth deformation gradient in time have to be aligned
along the partial first Piola-Kirchhoff stress tensor, suitably
transformed by the elastic portion of the deformation gradient.

The treatment presented in Sections~\ref{eu-viscoelastic-solid} and
\ref{eu-stress-dependent-growth} involve the notion of deformation
gradients and are thus only applicable to the solid collagen phase for
reasons discussed. In contrast, the analyses for the fluid and solute
phases below are carried out in $\Omega_{t}$, and we work with their
respective velocities as primitive variables characterising their
motion.
 
\subsection{A Newtonian fluid}
\label{eu-newtonian-fluid}

As alluded to toward the end of Section~\ref{eu-energy-dependent%
  -source}, one of the mechanisms underlying the rate-dependent
behaviour of soft collagenous tissues is the inherent viscosity of the
extra-cellular fluid. In this section, we derive the constitutive
relationship for a Newtonian (viscous) fluid from thermodynamic
considerations.

For the terms arising from the fluid phase, the reduced dissipation
inequality (\ref{eu-reduced-dissipation-1}) requires that,

\begin{equation*}
\rho^{\mathrm{f}} \dot{\psi^{\mathrm{f}}}
-\Bsigma^{\mathrm{f}}\colon\mathrm{grad}\left(\bv^{\mathrm{f}}\right)+
\rho^{\mathrm{f}} \mathrm{grad}\left(\psi^{\mathrm{f}}\right)
\cdot\bv^{\mathrm{f}} +\rho^{\mathrm{f}}
\bq^{\mathrm{f}}\cdot\bv^{\mathrm{f}} \leq 0.
\end{equation*}

\noindent We will return to the last term, $\rho^{\mathrm{f}}
\bq^{\mathrm{f}} \cdot \bv^{\mathrm{f}}$, when considering interaction
forces between the solid and fluid phases in
Section~(\ref{eu-interaction-forces}). Now, we turn our attention to
the remaining portion of the dissipation inequality for the fluid:

\begin{equation}
\rho^{\mathrm{f}} \dot{\psi^{\mathrm{f}}} -\Bsigma^{\mathrm{f}} \colon
\mathrm{grad}\left(\bv^{\mathrm{f}}\right)+ \rho^{\mathrm{f}}
\mathrm{grad} \left(\psi^{\mathrm{f}}\right) \cdot\bv^{\mathrm{f}}
\leq 0.
\label{reduced-dissipation-fluid}
\end{equation}

Assuming that the Helmholtz free energy of the fluid depends only on
its current concentration, i.e, $\psi^{\mathrm{f}} = \frac{1}
{\tilde{\rho^{\mathrm{f}}}} \hat{\psi^{\mathrm{f}}}
(\rho^{\mathrm{f}})$, where $\tilde{\rho^{\mathrm{f}}}$ is the
intrinsic density of the fluid, and substituting this form into
(\ref{reduced-dissipation-fluid}), we have,

\begin{equation*}
\frac{\rho^{\mathrm{f}}}{\tilde{\rho^{\mathrm{f}}}} \frac{\partial
  \hat{\psi^{\mathrm{f}}}}{\partial \rho^{\mathrm{f}}}
\dot{\rho^{\mathrm{f}}} -\Bsigma^{\mathrm{f}}\colon\mathrm{grad}
\left(\bv^{\mathrm{f}}\right) +
\frac{\rho^{\mathrm{f}}}{\tilde{\rho^{\mathrm{f}}}} \frac{\partial
  \hat{\psi^{\mathrm{f}}}}{\partial \rho^{\mathrm{f}}}
\mathrm{grad}\left(\rho^{\mathrm{f}}\right) \cdot\bv^{\mathrm{f}} \leq
0.
\end{equation*}

\noindent Invoking the balance of mass (\ref{eu-localbalanceofmass})
for the fluid, with $\pi^\mathrm{f} = 0$, we obtain,

\begin{equation*}
\frac{\rho^{\mathrm{f}}}{\tilde{\rho^{\mathrm{f}}}} \frac{\partial
  \hat{\psi^{\mathrm{f}}}}{\partial \rho^{\mathrm{f}}} \left(-
\rho^{\mathrm{f}} \mathrm{div}\left(\bv^{\mathrm{f}} \right)\right)
-\Bsigma^{\mathrm{f}}\colon\mathrm{grad} \left(\bv^{\mathrm{f}}\right)
\leq 0.
\end{equation*}

\noindent Finally, recalling the definition of pressure in terms of
the Helmholtz free energy at fixed temperature from classical
thermodynamics, we have the following form of the reduced dissipation
inequality for the fluid,

\begin{equation*}
\bigg( - \cancelto{p^{\mathrm{f}}(\rho^{\mathrm{f}})}
     {\frac{(\rho^{\mathrm{f}})^2}{\tilde{\rho^{\mathrm{f}}}}
       \frac{\partial \hat{\psi^{\mathrm{f}}}}{\partial
         \rho^{\mathrm{f}}}} \bone -\Bsigma^{\mathrm{f}} \bigg) \colon
     \mathrm{grad}\left(\bv^{\mathrm{f}} \right) \leq 0,
\end{equation*}

\noindent since \mbox{$\mathrm{div}(\bullet) = \bone \colon
  \mathrm{grad}(\bullet)$}, where $\bone$ is the second order identity
tensor.

A suitable form for the partial Cauchy stress tensor of the fluid
phase motivated by the above inequality is: $\Bsigma^{\mathrm{f}} = -
p^{\mathrm{f}} \bone + 2 \mu^{\mathrm{f}} \mathrm{grad}\left(
\bv^{\mathrm{f}} \right)$, where $\mu^{\mathrm{f}}$, the viscosity of
the fluid, is a non-negative scalar. Setting it to zero results in the
case of an ideal fluid. Furthermore, recognising from the balance of
angular momentum that $\Bsigma^{\mathrm{f}}$ is symmetric, we rewrite
the above constitutive relationship in terms of the fluid rate of
deformation tensor,~$\bd^{\mathrm{f}}$:

\begin{equation}
\Bsigma^{\mathrm{f}} = - p^{\mathrm{f}} \bone + 2
\mu^{\mathrm{f}}\bd^{\mathrm{f}},
\label{eu-stress-rel-fluid}
\end{equation}

\noindent which defines the behaviour of a classical Newtonian fluid.

\subsection{Frictional interaction forces}
\label{eu-interaction-forces}

We now focus our attention on the following two terms that were left
unaccounted for when exploring the implications of the terms arising
from the solid and fluid phases in the reduced dissipation inequality
(\ref{eu-reduced-dissipation-1}) in
Sections~\ref{eu-viscoelastic-solid} and \ref{eu-newtonian-fluid}:

\begin{equation}
\rho^{\mathrm{c}} \bq^{\mathrm{c}}\cdot\bv^{\mathrm{c}} +
\rho^{\mathrm{f}} \bq^{\mathrm{f}}\cdot\bv^{\mathrm{f}} \leq 0.
\label{eu-interactionremaining}
\end{equation}

In the time-scales of experiments studying the mechanical response of
tissues, growth is not usually significant, i.e. $\pi^\mathrm{c} =
\pi^\mathrm{f} = 0$. Applying this simplification to the summation
relationship between the interaction forces
(\ref{eu-summationrelationmomentum}), we see that,

\begin{equation}
\rho^\mathrm{c} \bq^\mathrm{c} = - \rho^\mathrm{f} \bq^\mathrm{f}.
\label{eu-simplifiedsummationrelation}
\end{equation}

\noindent Substituting Equation~(\ref{eu-simplifiedsummationrelation})
in (\ref{eu-interactionremaining}), we obtain,

\begin{equation}
\rho^{\mathrm{c}} \bq^{\mathrm{c}}\cdot \left(\bv^{\mathrm{c}} -
\bv^{\mathrm{f}}\right) \leq 0.
\label{eu-frictionresult}
\end{equation}

Relationships (\ref{eu-simplifiedsummationrelation}) and
(\ref{eu-frictionresult}) indicate that a suitable form for the
interaction forces between the solid collagen and fluid phases are:

\begin{equation}
\rho^{\mathrm{c}} \bq^{\mathrm{c}} = - \rho^{\mathrm{f}}
\bq^{\mathrm{f}} = -\bD^{\mathrm{fc}} \left(\bv^{\mathrm{c}} -
\bv^{\mathrm{f}}\right),
\label{eu-interactionforces}
\end{equation}

\noindent where $\bD^{\mathrm{fc}}$ is a positive semi-definite
frictional coefficient tensor.\footnote{This is only the simplest form
  possible. See \citet{massoudi03} for a thorough review of various
  thermodynamically and phenomenologically motivated forms for this
  interaction force.} This frictional interaction between the phases
is the basis for the time-dependent mechanical response observed in
the biphasic computations presented in
Section~\ref{biphasic-examples-2}.

\subsection{Diffusive solute fluxes}
\label{eu-solute-velocity}

Since the solute species are present only in low concentrations in
solution, they are not capable of exerting significant forces upon
other phases, and do not bear appreciable stress. Consequently, for
the solutes, the reduced dissipation inequality~(\ref{eu-reduced%
  -dissipation-1}) requires only that,

\begin{equation*}
\sum \left( \rho^{\mathrm{s}} \dot{\psi^{\mathrm{s}}} +
\rho^{\mathrm{s}} \mathrm{grad}\left(\psi^{\mathrm{s}}\right)
\cdot\bv^{\mathrm{s}} \right) \leq 0,
\end{equation*}

\noindent where the summation above is carried out over all solute
species. A sufficient condition to satisfy the above inequality is to
specify the following form for the mass flux of each species s,

\begin{equation}
\rho^{\mathrm{s}} \bv^{\mathrm{s}} = -\bD^{\mathrm{s}}
\mathrm{grad}\left(\psi^{\mathrm{s}}\right)
\label{eu-diffusive-flux}
\end{equation}

\noindent and only allow for solute interactions which decrease their
overall free energy density. Here, $\bD^{\mathrm{s}}$ is a positive
semi-definite mobility, and Equation~(\ref{eu-diffusive-flux}) results
in classical diffusive flux of the solute species when their Helmholtz
free energies are purely functions of their concentrations.

%

% Local Variables:
% TeX-master: "thesis"
% mode: latex
% mode: flyspell
% End:
