# A universal phase-field mixture representation of thermodynamics
and shock wave mechanics in porous soft biologic continua

J.D. Clayton<EMAIL_ADDRESS>Terminal Effects Division, DEVCOM
ARL, Aberdeen Proving Ground, MD 21005-5066 USA

###### Abstract

A continuum mixture theory is formulated for large deformations, thermal
effects, phase interactions, and degradation of soft biologic tissues. Such
tissues consist of one or more solid and fluid phases and can demonstrate
nonlinear anisotropic elastic, viscoelastic, thermoelastic, and poroelastic
physics. Under extremely large or rapid deformations, for example impact or
shock loading, tissues may fracture, tear, or rupture. Mechanisms are
encompassed in a universal, thermodynamically consistent formulation that
combines the continuum theory of mixtures with phase-field mechanics of
fracture. A metric tensor of generalized Finsler space supplies geometric
insight on effects rearrangements of microstructure, for example degrading
collagen fibers. Governing equations are derived, and energy potentials and
kinetic laws posited, for generic soft porous tissues with solid and liquid or
gas phases. Shock waves are modeled as singular surfaces; Hugoniot states and
shock decay are studied analytically. Suitability of the framework for
representing blood, skeletal muscle, and liver is demonstrated. Insight into
physics presently unresolved by experiments is obtained.

Suggested keywords

###### pacs:

87.10.Pq, 87.19.R-, 62.20.D-, 46.50.+a, 46.40.Cd

##  I  Introduction

Constitutive models describe the complex physics of soft biological tissue—for
example, skin, muscle, connective tissue, blood vessels, and the internal
organs—when subject to mechanical and thermal stimuli. These materials are
often simultaneously nonlinear anisotropic elastic, viscoelastic,
thermoelastic, and poroelastic [1, 2]. Large deformations manifesting
nonlinear mechanical response occur even under normal physiological activity.
Medical events involving disease and surgery incur cutting or tearing (i.e.,
fracture) of the tissue. For traumatic scenarios involving dynamic impact or
blast, even more extreme deformations and rapid loading rates arise [1, 3, 4,
5, 6].

Microstructures are often hierarchical and highly complex [1, 2]. Materials
consist of multiple solid- and fluid-like phases. The liver includes the
spongy parenchyma, connective tissue, blood vessels and ducts, plus blood and
bile [7, 8]. Skeletal and cardiac muscles contain fibers (i.e., cells),
connective tissue, blood and interstitial fluid [2, 9]. Skin includes dermal
tissue layers, cells, and interstitial fluid [10]. Lung includes parenchyma
(e.g., alveoli and ducts), air, and surfactant fluid, as well as stiffer
elements of the bronchiole structure [1]. Membraneous tissues (e.g., in skin)
contain ground substance, elastin, and collagen fibers [11]. Blood has cells
immersed in the extracellular plasma [2, 12, 13]. Cells themselves include
solid-like walls and internal fluids [14, 15].

Most constitutive treatments of mechanics of soft tissues focus on their
nonlinear elastic response [16, 9, 17]. Effects of fluids are included
implicitly in the energy densities and parameters, perhaps augmented with
viscoelasticity or other dissipation elements [18, 19]. Even in the non-
dissipative case, sophisticated models are needed to account for anisotropy
and nonlinearity, for example due to collagen fiber distributions [20, 21,
22]. If tearing or degradation occur, continuum damage mechanics theories are
fairly standard [23, 24], whereby a phenomenological kinetic equation is
prescribed for an internal state variable measuring local loss of stress-
bearing capacity.

In contrast to the former, the phase-field approach has been more recently
advocated for modeling tearing of soft biologic tissues [25, 26, 27]. Diffuse-
interface modeling has witnessed use over a broad application space; models
with elasticity include phase changes, twinning, and dislocations in crystals
[28, 29, 30, 31]. The method appeared for brittle fracture around 25 years ago
[32, 33, 34, 35, 36], typically for small-strain, isotropic elasticity theory.
Finite-strain approaches for fracture appeared thereafter [37, 38]. Fluid
cavitation [39] and fluid injection [40] have been modeled.

Continuum mixture theories and porous media theories emerged in the mid-20th
century [41, 42, 43, 44]. Microstructure details are smeared, but these models
distinguish stresses and deformations of each constituent and capture
transfers of mass, momentum, and energy between phases. Mixture- and porous-
media theories introduce length and time scales that are absent in single-
phase viscoelasticity [45, 46]. Although porous media models have been used
for soft tissues [45, 47, 7, 48, 49, 8], they have not been combined with the
phase-field fracture approach for the class of soft biologic materials. Finite
strains [50] and thermal exchange effects [51] have, however, been considered
previously with phase-field fracture of brittle porous media.

The theory of porous media classically models a two-phase system of one solid
and one fluid, though models with multiple fluids exist. The fully dense solid
and fluid phases are usually treated as incompressible, but the mixture
overall is compressible as fluid is locally squeezed out [45, 47, 7]. Solid
and fluid are usually, but not always [8], treated as inviscid individually,
but interactions of viscous origin between phases are captured by permeability
and dissipation from fluid transport (e.g., Darcy’s law). An effective stress
principle [52] decomposes total stress into solid skeleton stress and pore
pressure. Constitutive equations are supplied for the solid skeleton (i.e.,
drained material) and fully dense fluid rather than for partial stresses of
mixture components. This method has proven successful for soil and rock
mechanics in which the distinction between solid and fluid(s) is clear (e.g.,
water flowing through sand or porous rock), and it has been used elsewhere for
porous tissues [45, 47, 53].

The effective stress principle with solid skeleton concept is not implemented
here for several reasons. First, tissue may consist of multiple solid and
fluid phases. For example, the liver has the parenchyma, blood vessels, and
connective tissue (which might be tied or displace independently) and distinct
fluids of blood and bile. Second, designation of a constituent as a “solid” or
“fluid” may be ambiguous, for example, the extracellular matrix of tissue or
blood demonstrating both viscous (fluid-like) and viscoelastic (solid-like)
physics [12, 13, 54]. With no single true solid or fluid phase, the soil
mechanics analogy is tenuous. Cracking in ambiguously soft materials also
shows richer physics than in brittle elastic solids due to stochastic bond
opening and closure (i.e., healing) [55].

A mixture theory similar to that of Refs. [56, 43, 57, 58] is advocated here
instead, with constituent energies formulated on a per-unit-mass basis.
Complexity in thermodynamic derivations is reduced since partial stresses
rather than effective stress and pore pressure are used. However, potential
difficulty can arise in ascertaining properties of isolated phases as
experiments might measure responses of only some (e.g., fluid) phase(s) and
the mixture as a whole. This issue is rectified by example later. The typical
assumption of incompressibility of individual phases is abandoned to resolve
bulk sound waves and longitudinal shocks. Incompressibility is also
unrealistic if materials dilate or undergo cavitation (i.e., tensile damage).

A unified framework is formulated in Secs. II and III, newly synthesizing
continuum mixture theory, nonlinear anisotropic thermoelasticity, viscosity,
viscoelasticity, and phase-field fracture mechanics for biologic systems. A
universal, thermodynamically consistent theory containing all such features
relevant to soft-tissue mechanics for intense rapid loading does not seem to
exist. Another unique aspect is generalized Finsler geometry describing the
material manifold with evolving microstructure [59].

Traditional continuum models are couched in ambient Euclidean 3-space.
Residual stresses from growth and remodeling have been considered using
Riemannian geometry with a metric of non-vanishing curvature tensor [60, 61].
A recent approach [62, 63] accounts for effects of microstructure using a
generalized Finsler metric [64] as opposed to a Riemannian metric. In addition
to residual stress and growth, the Finsler-geometric approach on the material
manifold has been used to describe changes in local configurations of collagen
fibers as tissues degrade and tear [59]. A Finsler metric was used elsewhere
to quantify effects of fiber orientations on mechanical responses of soft
tissues in a discrete bond-based model [65] rather than a continuum physics
approach as herein. The Finsler-geometric theory has similarities to
micropolar theory [66, *ikeda1981, 68], but the former is more general [62].
Much past work focused on hard crystalline solids [69, 70].

Shock waves in soft tissues are analyzed here with the mixture theory and
constitutive frameworks. Hugoniot solutions and evolution equations are
derived for compressive shocks, extending prior works [71, 72, 73, 74] to
simultaneously account for internal state variables (e.g., order parameters)
and the Finsler metric that, in the present application, is transformed to an
osculating Riemannian metric [75, 76, 59]. Dissipation from viscoelasticity,
momentum and energy transfer between phases, and degradation due to shear-
induced tearing all potentially affect shock amplitudes over time. Treatment
of shocks as singular surfaces of velocity differs from those of continuous
waves in nonlinear materials [77]. Shock structures (e.g., shapes of
continuous wave forms) have been analyzed elsewhere for nonlinear solids [78]
and fluids [79, 58].

Modeled in Sec. IV are three soft-tissue systems: skeletal muscle with
interstitial fluid, liver with blood, and lung with air. Solutions to physical
problems involving tension, compression, or shock-wave loading grant new
understanding of the physics demonstrated by these materials. Concluding
remarks give closure in Sec. V.

##  II  Governing Equations

The present theory builds significantly on the finite-strain mixture theories
of Refs. [43, 44, 74]. Two enhancements are incorporated here. First, time-
dependent general internal state vectors are introduced. Elements of these
vectors are later associated with history dependent mechanisms in the material
microstructure, namely viscoelasticity, active tension, and damage, all
generally anisotropic. A dependence of energy potentials on internal state and
the gradient of internal state is permitted, with terms associated with
internal state gradients properly incorporated in the balance of energy and
boundary conditions. This enables a phase-field type representation when the
variables are interpreted as order parameters [80, 28, 30], suitable for
modeling regularized fracture [37, 38, 50, 26]. Second, metric tensors on
spatial and material manifolds are permitted to depend on internal state and
can be time-dependent. Distances measured in the material can include remnant
strains from dissipative processes, or biologic growth and remodeling if the
latter physics are resolved by state variables. If internal-state dependence
of metrics is explicit and distinct from coordinate dependence, a generalized
Finsler representation emerges [62, 63]. If internal-state dependence is
implicit and state vectors are (time-dependent) functions of position [81],
then an osculating Riemannian geometry [76, 59] is obtained. In the Riemannian
setting, similarities with Refs. [60, 61] are noted.

A mixture of $N\geq 1$ constituents is considered, where at time $t$, these
occupy a shared infinitesimal control volume $\rm{d}\Omega$ centered at
spatial position $\bf{x}$. Greek superscripts $\alpha=1,\ldots,N$ denote
constituents having reference coordinates $\bf{X}^{\alpha}$; constituents
coincident at $\bf{x}$ may occupy different $\bf{X}^{\alpha}$ due to diffusion
processes. Motions are

${\bf x}={\bm{\chi}}^{\alpha}({\bf X}^{\alpha},t).$ (1)

A spatial manifold comprising the material body, parameterized by coordinate
chart(s) $\\{x^{k}\\}$, is $\mathfrak{m}$. Referential manifold(s)
$\mathfrak{M}^{\alpha}$ are parameterized by coordinate chart(s)
$\\{(X^{\alpha})^{K}\\}$ corresponding to reference positions of constituent
$\alpha$. Denote by $\\{{\bm{\xi}}^{\alpha}({\bf x},t)\\}$ and
$\\{\bm{\Xi}^{\alpha}({\bf X}^{\alpha},t)\\}$ sets of internal state variables
viewed as auxiliary coordinates over $\mathfrak{m}$ and
$\mathfrak{M}^{\alpha}$, respectively. Metric tensors $\bf{g}$ and $\bf
G^{\alpha}$ with components $g_{ij}$ and $G^{\alpha}_{IJ}$ are permitted to be
coordinate- and time-dependent, of the forms

$\displaystyle{\bf g}={\bf g}({\bf x},t)=\tilde{\bf g}({\bf
x},\\{{\bm{\xi}}^{\alpha}({\bf x},t)\\}),$ (2) $\displaystyle{\bf
G}^{\alpha}={\bf G}^{\alpha}({\bf X}^{\alpha},t)=\tilde{\bf G}^{\alpha}({\bf
X}^{\alpha},\\{{\bm{\Xi}}^{\alpha}({\bf X}^{\alpha},t)\\}).$ (3)

The $\sim$ notation denotes the generalized Finsler description [64, 62, 63]
of metrics on $\mathfrak{m}$ and $\mathfrak{M}^{\alpha}$, whereas unadorned
versions are interpreted as osculating Riemannian metric tensors [75, 76, 59]
at any fixed time $t$. Determinants are written $g=\det{\bf g}$ and
$G^{\alpha}=\det{\bf G}^{\alpha}$. The partial time derivative at fixed $\bf
x$ is $\partial_{t}(\cdot)$; the material time derivative at fixed ${\bf
X}^{\alpha}$ is $D^{\alpha}_{t}(\cdot)$, related by

$D^{\alpha}_{t}(\cdot)=\partial_{t}(\cdot)+\nabla(\cdot)\cdot{\bm{\upsilon}}^{\alpha},\quad({\upsilon}^{\alpha})^{k}=D^{\alpha}_{t}({\chi}^{\alpha})^{k}.$
(4)

Particle velocity is ${\bm{\upsilon}}^{\alpha}$, and $\nabla(\cdot)$ is the
covariant derivative with respect to $\bf x$ where Christoffel symbols are
those of the Levi-Civita connection derived from $g_{ij}$. The covariant
derivative with respect to $\bf X^{\alpha}$ on $\mathfrak{M}^{\alpha}$ is
$\nabla_{0}^{\alpha}(\cdot)$. Spatial and referential gradient operators obey

$\nabla(\cdot)=\partial(\cdot)/\partial x^{k}\otimes{\bf
g}^{k},\quad\nabla_{0}^{\alpha}(\cdot)=\partial(\cdot)/\partial(X^{\alpha})^{K}\otimes{\bf
G}^{K}.$ (5)

When used on typical scalars, vectors, and tensors, $\partial_{t}(\cdot)$ and
$D^{\alpha}_{t}(\cdot)$ are performed with natural basis vectors ${\bf
g}_{k}=\partial{\bf x}/\partial x^{k}$ and ${\bf G}^{\alpha}_{K}=\partial{\bf
X}/\partial(X^{\alpha})^{K}$ held fixed with respect to $t$ at ${\bf x}$ and
$\bf X^{\alpha}$, respectively, so $\partial_{t}{\bf g}_{k}$ and
$D_{t}^{\alpha}{\bf G}^{\alpha}_{K}$ vanish in such circumstances. This
produces commutation rules:

$\nabla[\partial_{t}(\cdot)]=\partial_{t}[\nabla(\cdot)],\quad\nabla_{0}^{\alpha}[D^{\alpha}_{t}(\cdot)]=D_{t}^{\alpha}[\nabla_{0}^{\alpha}(\cdot)].$
(6)

The deformation gradient ${\bf F}^{\alpha}$ and Jacobian determinant
$J^{\alpha}$ are defined as follows, giving a relation between reference and
spatial gradient operators:

$\displaystyle(F^{\alpha})^{i}_{J}=\frac{\partial(\chi^{\alpha})^{i}}{\partial(X^{\alpha})^{J}},\quad
J^{\alpha}=\det[(F^{\alpha})^{i}_{J}]\sqrt{g/G^{\alpha}},$ (7)
$\displaystyle{\bf F}^{\alpha}=(F^{\alpha})^{i}_{J}{\bf g}_{i}\otimes{\bf
G}^{J},\quad\nabla^{\alpha}_{0}(\cdot)=\nabla(\cdot){\bf F}^{\alpha}.$ (8)

The velocity gradient ${\bf l}^{\alpha}$ and its trace are

${\bf l}^{\alpha}=\nabla{\bm{\upsilon}}^{\alpha}=D^{\alpha}_{t}{\bf
F}^{\alpha}({\bf
F^{\alpha}})^{-1},\quad\nabla\cdot{\bm{\upsilon}}^{\alpha}={\rm tr}{\bf
l}^{\alpha}.$ (9)

Spatial and material volume elements, ${\rm d}\Omega$ and ${\rm
d}\Omega^{\alpha}_{0}$, obey

${\rm d}\Omega({\bf x}({\bf X}^{\alpha},t),t)=J^{\alpha}({\bf
X}^{\alpha},t){\rm d}\Omega_{0}^{\alpha}({\bf X}^{\alpha},t).$ (10)

Time derivatives of local volume elements allow for time dependence of metric
tensors, extending Ref. [61] as

$\displaystyle\partial_{t}({\rm d}\Omega)={\textstyle{\frac{1}{2}}}{\rm
tr}(\partial_{t}{\bf g}){\rm d}\Omega=\partial_{t}(\ln\sqrt{g}){\rm d}\Omega,$
(11) $\displaystyle D^{\alpha}_{t}({\rm
d}\Omega_{0}^{\alpha})={\textstyle{\frac{1}{2}}}{\rm tr}(D^{\alpha}_{t}{\bf
G}^{\alpha}){\rm
d}\Omega_{0}^{\alpha}=D^{\alpha}_{t}(\ln\sqrt{G^{\alpha}}){\rm
d}\Omega_{0}^{\alpha},$ (12) $\displaystyle D^{\alpha}_{t}({\rm
d}\Omega)=(D^{\alpha}_{t}J^{\alpha}+J^{\alpha}D^{\alpha}_{t}\ln\sqrt{G^{\alpha}}){\rm
d}\Omega_{0}^{\alpha}$
$\displaystyle\qquad\quad\,\,\,=(\nabla\cdot{\bm{\upsilon}}^{\alpha}+\partial_{t}\ln\sqrt{g}){\rm
d}\Omega.$ (13)

In (13), $D^{\alpha}_{t}\sqrt{g/G^{\alpha}}$ is included in
$D_{t}^{\alpha}J^{\alpha}$, and $\partial_{t}g=D^{\alpha}_{t}g$ since $\nabla
g$ vanishes for a Levi-Civita connection [82].

Define the following, all generally fields of $({\bf x},t)$: partial Cauchy
stress tensor ${\bm{\sigma}}^{\alpha}$, traction vector ${\bf
t}^{\alpha}={\bm{\sigma}}^{\alpha}\cdot{\bf n}$, body force per unit mass
${\bf b}^{\alpha}$, partial internal energy per unit mass $u^{\alpha}$, heat
source per unit mass $r^{\alpha}$, heat flux vector ${\bf q}^{\alpha}$, mass
supply rate $c^{\alpha}$, linear momentum exchange ${\bf h}^{\alpha}$, and
energy exchange rate $\epsilon^{\alpha}$. The spatial mass density of
constituent $\alpha$ is $\rho^{\alpha}({\bf x},t)$. Referential mass density
at a fixed time $t=t_{0}$ is $\rho_{0}^{\alpha}({\bf X^{\alpha}})$ with metric
${\bf G}_{0}^{\alpha}({\bf X^{\alpha}})$ and mass element ${\rm
d}m_{0}^{\alpha}({\bf
X}^{\alpha})=\rho_{0}^{\alpha}\sqrt{G_{0}^{\alpha}/G^{\alpha}}{\rm
d}\Omega^{\alpha}_{0}$. At arbitrary $t$, the mass element is ${\rm
d}m^{\alpha}({\bf X}^{\alpha},t)=\rho^{\alpha}{\rm
d}\Omega=\rho^{\alpha}J^{\alpha}{\rm d}\Omega_{0}^{\alpha}$.

###  II.1 Continuous processes

A global energy balance for each constituent is posited on a region of
$\mathfrak{m}$ occupying control volume $\Omega$, enclosed by oriented
boundary $\partial\Omega$ with unit outward normal ${\bf n}$. Velocity and
heat flux normal to $\partial\Omega$ are
$\upsilon^{\alpha}_{n}={\bm{\upsilon}^{\alpha}}\cdot{\bf n}$ and
$q^{\alpha}_{n}={\bf q^{\alpha}}\cdot{\bf n}$. The global balance per
constituent extends that of Ref. [43] to account for working of generalized
tractions $\\{\bf z^{\alpha}\\}=\\{{\bm{\zeta}^{\alpha}}\cdot{\bf n}\\}$ work-
conjugate to rates of auxiliary variables
$\\{D_{t}^{\alpha}{\bm{\xi}}^{\alpha}\\}$ on $\partial\Omega$, similarly to
Ref. [80]:

$\displaystyle\partial_{t}\int_{\Omega}\rho^{\alpha}(u^{\alpha}+{\textstyle{\frac{1}{2}}}|{\bm{\upsilon}}^{\alpha}|^{2}){\rm{d}\Omega}$
$\displaystyle\qquad\qquad+\oint_{\partial\Omega}\rho^{\alpha}(u^{\alpha}+{\textstyle{\frac{1}{2}}}|{\bm{\upsilon}}^{\alpha}|^{2})\upsilon^{\alpha}_{n}{\rm{d}\partial\Omega}$
$\displaystyle=\oint_{\partial\Omega}({\bf
t}^{\alpha}\cdot{\bm{\upsilon}}^{\alpha}+\\{{\bf
z}^{\alpha}\\}\cdot\\{D^{\alpha}_{t}{\bm{\xi}}^{\alpha}\\}-q_{n}^{\alpha}){\rm{d}\partial\Omega}$
$\displaystyle\qquad\qquad+\int_{\Omega}[\rho^{\alpha}({\bf
b}^{\alpha}\cdot{\bm{\upsilon}}^{\alpha}+r^{\alpha})+{\bf
h}^{\alpha}\cdot{\bm{\upsilon}}^{\alpha}+\epsilon^{\alpha}$
$\displaystyle\qquad\qquad\qquad\,\,+c^{\alpha}(u^{\alpha}+{\textstyle{\frac{1}{2}}}|{\bm{\upsilon}}^{\alpha}|^{2})]{\rm{d}\Omega}.$
(14)

Angular momentum exchange between constituents [43] is omitted herein. Field
variables in (14) are assumed sufficiently smooth on $\Omega$ such that the
divergence theorem applies on ${\mathfrak{m}}$. Energy conservation under
rigid motions (i.e., various time-dependent translations and rigid-body
rotations [83, 84]) of the entire mixture then furnishes local conservation
laws for mass, momenta, and energy:

$\displaystyle\partial_{t}\rho^{\alpha}+\nabla\cdot(\rho^{\alpha}{\bm{\upsilon}}^{\alpha})+\rho^{\alpha}\partial_{t}\ln\sqrt{g}=c^{\alpha},$
(15) $\displaystyle\nabla\cdot{\bm{\sigma}}^{\alpha}+\rho^{\alpha}{\bf
b}^{\alpha}+{\bf
h}^{\alpha}=\rho^{\alpha}D^{\alpha}_{t}{\bm{\upsilon}}^{\alpha},\quad{\bm{\sigma}}^{\alpha}=({\bm{\sigma}}^{\alpha})^{\mathsf{T}},$
(16)
$\displaystyle\rho^{\alpha}D_{t}^{\alpha}u^{\alpha}={\bm{\sigma}}^{\alpha}:\nabla{\bm{\upsilon}}^{\alpha}+\nabla\cdot(\\{{\bm{\zeta}}^{\alpha}\\}\cdot\\{D_{t}^{\alpha}{\bm{\xi}}^{\alpha}\\})$
$\displaystyle\qquad\qquad\qquad\qquad\quad\,-\nabla\cdot{\bf
q}^{\alpha}+\rho^{\alpha}r^{\alpha}+\epsilon^{\alpha}.$ (17)

Denote by $\eta^{\alpha}$ the local entropy per unit mass and
$\theta^{\alpha}>0$ the absolute temperature of constituent $\alpha$. As in
Refs. [43, 74], an entropy inequality for each constituent is avoided in lieu
of an inequality for the whole mixture:

$\displaystyle\partial_{t}\int_{\Omega}\sum_{\alpha}\rho^{\alpha}\eta^{\alpha}{\rm{d}\Omega}+\oint_{\partial\Omega}\sum_{\alpha}\rho^{\alpha}\eta^{\alpha}\upsilon_{n}^{\alpha}{\rm{d}\partial\Omega}$
$\displaystyle\qquad\qquad\geq\int_{\Omega}\sum_{\alpha}\frac{\rho^{\alpha}r^{\alpha}}{\theta^{\alpha}}{\rm{d}\Omega}-\oint_{\partial\Omega}\sum_{\alpha}\frac{q_{n}^{\alpha}}{\theta^{\alpha}}{\rm{d}\partial\Omega}.$
(18)

From the divergence theorem, (15), and localization 111Misprints in eqs.
(2.15), (4.1), (4.19)1, (4.30)3, (4.45), (4.63), (5.2)2, (5.13)2 of Ref. [74]
are corrected herein.,

$\sum_{\alpha}[\rho^{\alpha}D^{\alpha}_{t}\eta^{\alpha}+\frac{\nabla\cdot{\bf
q}^{\alpha}}{\theta^{\alpha}}-\frac{{\bf
q^{\alpha}}\cdot\nabla\theta^{\alpha}}{(\theta^{\alpha})^{2}}-\frac{\rho^{\alpha}r^{\alpha}}{\theta^{\alpha}}+c^{\alpha}\eta^{\alpha}]\geq
0.$ (19)

Let $\psi^{\alpha}$ be Helmholtz free energy per unit mass, whereby
substitution of (II.1) into (19) yields

$\displaystyle\psi^{\alpha}=u^{\alpha}-\theta^{\alpha}\eta^{\alpha},$ (20)
$\displaystyle\sum_{\alpha}({1}/{\theta^{\alpha}})[{\bm{\sigma}}^{\alpha}:\nabla{\bm{\upsilon}}^{\alpha}+\nabla\cdot(\\{{\bm{\zeta}}^{\alpha}\\}\cdot\\{D_{t}^{\alpha}{\bm{\xi}}^{\alpha}\\})$
$\displaystyle\qquad\qquad-({{\bf
q^{\alpha}}\cdot\nabla\theta^{\alpha}})/{\theta^{\alpha}}+\epsilon^{\alpha}+c^{\alpha}\theta^{\alpha}\eta^{\alpha}$
$\displaystyle\qquad\qquad\qquad-\rho^{\alpha}(D^{\alpha}_{t}\psi^{\alpha}+\eta^{\alpha}D^{\alpha}_{t}\theta^{\alpha})]\geq
0.$ (21)

###  II.2 Singular surfaces

Now consider a propagating singular surface $\Sigma(t)$ in $\mathfrak{m}$,
with image $\Sigma^{\alpha}(t)$ in $\mathfrak{M}^{\alpha}$. The Eulerian
function $\phi$ defines this surface, from which the unit normal $\bf n$ and
Eulerian speed $\mathcal{U}>0$ in the direction of $\bf n$ follow:

$\phi({\bf x},t)=0,\quad{\bf
n}=\nabla\phi/|\nabla\phi|,\quad\mathcal{U}=-\partial_{t}\phi/|\nabla\phi|.$
(22)

Let $(\cdot)^{+}$ and $(\cdot)^{-}$ label limiting values of $(\cdot)$ as
$\Sigma$ is approached from either side; $\bf n$ is directed from the
$(\cdot)^{-}$ side (behind $\Sigma$) to the $(\cdot)^{+}$ side. The jump
across $\Sigma$ is

$\llbracket(\cdot)\rrbracket=(\cdot)^{-}-(\cdot)^{+}.$ (23)

Across a shock front $\Sigma$, each ${\bm{\chi}}^{\alpha}$ is continuous, but
space-time derivatives of ${\bm{\chi}}^{\alpha}$ are not necessarily so. Nor
always are other field variables such as $\rho^{\alpha}$,
${\bm{\sigma}}^{\alpha}$, $\theta^{\alpha}$, etc. With $\bf n$ defined per
(22), the normal velocity, heat flux, and tractions are
$\upsilon_{n}^{\alpha}={\bm{\upsilon}}^{\alpha}\cdot{\bf n}$,
$q_{n}^{\alpha}={\bf q}^{\alpha}\cdot{\bf n}$, ${\bf
t}^{\alpha}={\bm{\sigma}}^{\alpha}\cdot{\bf n}$, and $\\{\bf
z^{\alpha}\\}=\\{{\bm{\zeta}^{\alpha}}\cdot{\bf n}\\}$. Conservation laws for
mass, linear momentum, and energy, and the entropy imbalance, across $\Sigma$
are derived using principles set forth in Refs. [42, 86]. Specifically, a
closed region $\Omega$ of $\mathfrak{m}$ is partitioned by $\Sigma$ at an
instant in time into two sub-bodies, within which the local continuum balance
laws (15)–(II.1) and (19) hold. These continuum laws are integrated over
$\Omega$ and then over each sub-body, the latter accounting for possible
boundary contributions on $\Sigma$ from tractions and heat flux. Integrals are
assumed to be physically meaningful even if $\mathfrak{m}$ is non-Euclidean.
The total mass rate, linear momentum rate, energy rate, and entropy production
rate are required to be equal for $\Omega$ and the summed contributions of
each sub-body with boundary $\Sigma$. Application of a form of Reynolds’
transport theorem [86, 87] derived from the divergence theorem and (13) then
gives analogs of (15), (16), and (II.1) across $\Sigma$ for each $\alpha$:

$\displaystyle\llbracket\rho^{\alpha}(\upsilon_{n}^{\alpha}-\mathcal{U})\rrbracket=0,$
(24)
$\displaystyle\llbracket\rho^{\alpha}{\bm{\upsilon}^{\alpha}}(\upsilon_{n}^{\alpha}-\mathcal{U})\rrbracket=\llbracket{\bf
t}^{\alpha}\rrbracket,$ (25)
$\displaystyle\llbracket\rho^{\alpha}(u^{\alpha}+{\textstyle{\frac{1}{2}}}|{\bm{\upsilon}}^{\alpha}|^{2})(\upsilon_{n}^{\alpha}-\mathcal{U})\rrbracket$
$\displaystyle\qquad\quad=\llbracket{\bf
t}^{\alpha}\cdot{\bm{\upsilon}}^{\alpha}+\\{{\bf
z}^{\alpha}\\}\cdot\\{D^{\alpha}_{t}{\bm{\xi}}^{\alpha}\\}-q_{n}^{\alpha}\rrbracket.$
(26)

A jump equation for angular momentum can be derived, but it encompasses
nothing beyond (25). Similarly, arguments in Ref. [86] applied to (19) furnish
a local entropy inequality across $\Sigma$ for the mixture:

$\sum_{\alpha}\llbracket\rho^{\alpha}\eta^{\alpha}(\mathcal{U}-\upsilon_{n}^{\alpha})-q_{n}^{\alpha}/\theta^{\alpha}\rrbracket\geq
0.$ (27)

Now consider one-dimensional (1-D) loading conditions: $n_{k}\rightarrow
n_{1}=1$, $x^{k}\rightarrow x^{1}=x$,
$(\chi^{\alpha})^{k}\rightarrow(\chi^{\alpha})^{1}=\chi^{\alpha}$,
$(F^{\alpha})^{i}_{J}\rightarrow(F^{\alpha})^{1}_{1}=\partial\chi^{\alpha}/\partial
X^{\alpha}=F^{\alpha}$,
$(\upsilon^{\alpha})^{k}\rightarrow(\upsilon^{\alpha})^{1}=\upsilon^{\alpha}_{n}=\upsilon^{\alpha}$,
$(t^{\alpha})_{k}\rightarrow(t^{\alpha})_{1}=(\sigma^{\alpha})^{1}_{1}=t^{\alpha}$,
$\\{(z^{\alpha})_{k}\\}\rightarrow\\{(z^{\alpha})_{1}\\}=\\{(\zeta^{\alpha})^{1}_{1}\\}=\\{z^{\alpha}\\}$,
$\\{(\xi^{\alpha})^{k}\\}\rightarrow\\{(\xi^{\alpha})^{1}\\}=\\{\xi^{\alpha}\\}$,
and $(q^{\alpha})^{k}\rightarrow(q^{\alpha})^{1}=q^{\alpha}_{n}=q^{\alpha}$.
The shock is planar, with Eulerian speed $\mathcal{U}$. Eulerian forms of
Rankine-Hugoniot equations (24)–(27) reduce to

$\displaystyle\llbracket\rho^{\alpha}(\upsilon^{\alpha}-\mathcal{U})\rrbracket=0,$
(28)
$\displaystyle\llbracket\rho^{\alpha}{\upsilon^{\alpha}}(\upsilon^{\alpha}-\mathcal{U})\rrbracket=\llbracket{t}^{\alpha}\rrbracket,$
(29)
$\displaystyle\llbracket\rho^{\alpha}(u^{\alpha}+{\textstyle{\frac{1}{2}}}|{\upsilon}^{\alpha}|^{2})(\upsilon^{\alpha}-\mathcal{U})\rrbracket$
$\displaystyle\qquad\quad=\llbracket{t}^{\alpha}{\upsilon}^{\alpha}+\\{{z}^{\alpha}\\}\\{D^{\alpha}_{t}{\xi}^{\alpha}\\}-q^{\alpha}\rrbracket,$
(30)
$\displaystyle\sum_{\alpha}\llbracket\rho^{\alpha}\eta^{\alpha}(\mathcal{U}-\upsilon^{\alpha})-q^{\alpha}/\theta^{\alpha}\rrbracket\geq
0.$ (31)

A Lagrangian description of a singular surface, denoted $\Sigma^{\alpha}(t)$
for each constituent $\alpha$, is, analogously to (22),

$\displaystyle\Phi^{\alpha}({\bf X^{\alpha}},t)=0,\qquad{\bf
N}^{\alpha}=\nabla^{\alpha}_{0}\Phi^{\alpha}/|\nabla^{\alpha}_{0}\Phi^{\alpha}|,$
$\displaystyle\mathcal{U}^{\alpha}=-D^{\alpha}_{t}\Phi^{\alpha}/|\nabla_{0}^{\alpha}\Phi^{\alpha}|.$
(32)

In 1-D, Eulerian and Lagrangian shock speeds are [73]

${\mathcal{U}}(t)={\rm d}\Sigma(t)/{\rm
d}t,\quad{\mathcal{U}}^{\alpha}(t)={\rm d}\Sigma^{\alpha}(t)/{\rm d}t.$ (33)

From 1-D inverse motion
$X^{\alpha}(\chi^{\alpha},t)=(\chi^{\alpha})^{-1}(x,t)$ [73],

$\Sigma^{\alpha}(t)=(\chi^{\alpha})^{-1}(\Sigma(t),t)\Rightarrow{\mathcal{U}}^{\alpha}=(F^{\alpha})^{-1}({\mathcal{U}}-\upsilon^{\alpha}),$
(34)

with $(F^{\alpha})^{-1}=\partial(\chi^{\alpha})^{-1}/\partial x$. Now assume a
continuous referential density field $\rho_{0}^{\alpha}(X^{\alpha})$ exists
and can be related in 1-D to any other spatial density field
$\rho^{\alpha}(x,t)$ via $\rho_{0}^{\alpha}=F^{\alpha}\rho^{\alpha}$.
Sufficient conditions for the latter consistent with mass conservation are
$c^{\alpha}=\rho^{\alpha}\partial_{t}\ln\sqrt{g}$ with
$g(\chi^{\alpha}(X^{\alpha},t),t)=G^{\alpha}(X^{\alpha},t)$. In this case,
(28) with (34) and $\llbracket\rho^{\alpha}_{0}\rrbracket=0$ produces
$\llbracket\mathcal{U}^{\alpha}\rrbracket=0$. Thus, noting that
$\mathcal{U}^{\alpha}$ is intrinsic and (34) should hold at all
$(\cdot)^{\pm}$ states, the following Lagrangian forms of the 1-D Rankine-
Hugoniot equations are derived:

$\displaystyle\rho^{\alpha}_{0}{\mathcal{U}}^{\alpha}\llbracket
1/\rho^{\alpha}\rrbracket=-\llbracket\upsilon^{\alpha}\rrbracket,$ (35)
$\displaystyle\rho^{\alpha}_{0}{\mathcal{U}}^{\alpha}\llbracket\upsilon^{\alpha}\rrbracket=-\llbracket{t}^{\alpha}\rrbracket,$
(36) $\displaystyle\rho^{\alpha}_{0}{\mathcal{U}}^{\alpha}\llbracket
u^{\alpha}+{\textstyle{\frac{1}{2}}}|{\upsilon}^{\alpha}|^{2}\rrbracket$
$\displaystyle\qquad\quad=-\llbracket{t}^{\alpha}{\upsilon}^{\alpha}+\\{{z}^{\alpha}\\}\\{D^{\alpha}_{t}{\xi}^{\alpha}\\}-q^{\alpha}\rrbracket,$
(37)
$\displaystyle\sum_{\alpha}(\rho^{\alpha}_{0}{\mathcal{U}}^{\alpha}\llbracket\eta^{\alpha}\rrbracket-\llbracket
q^{\alpha}/\theta^{\alpha}\rrbracket)\geq 0.$ (38)

Routine algebra produces velocity jumps and Lagrangian shock speeds in terms
of jumps in stress and mass density:

$\llbracket\upsilon^{\alpha}\rrbracket=(\llbracket
t^{\alpha}\rrbracket\llbracket
1/\rho^{\alpha}\rrbracket)^{1/2},\quad\rho^{\alpha}_{0}\mathcal{U}^{\alpha}=(\llbracket
t^{\alpha}\rrbracket/\llbracket 1/\rho^{\alpha}\rrbracket)^{1/2}.$ (39)

Also, since $\rho_{0}^{\alpha}/\rho^{\alpha}=F^{\alpha}$ by construction,

$\llbracket\upsilon^{\alpha}\rrbracket=-{\mathcal{U}}^{\alpha}\llbracket
F^{\alpha}\rrbracket,\qquad\llbracket
t^{\alpha}\rrbracket=\rho_{0}^{\alpha}({\mathcal{U}}^{\alpha})^{2}\llbracket
F^{\alpha}\rrbracket.$ (40)

The Rankine-Hugoniot energy balance follows by eliminating particle and shock
velocities from (II.2):

$\llbracket u^{\alpha}\rrbracket=\langle t^{\alpha}\rangle\llbracket
1/\rho^{\alpha}\rrbracket-\frac{\llbracket\\{{z}^{\alpha}\\}\\{D^{\alpha}_{t}{\xi}^{\alpha}\\}-q^{\alpha}\rrbracket}{(\llbracket
t^{\alpha}\rrbracket/\llbracket 1/\rho^{\alpha}\rrbracket)^{1/2}}.$ (41)

where $\langle(\cdot)\rangle=\frac{1}{2}[(\cdot)^{+}+(\cdot)^{-}]$ is the
average across $\Sigma^{\alpha}$.

The displacement derivative $\delta_{t}(\cdot)$ is defined as follows in 1-D
[42], consistently with (4) and (33):

$\delta_{t}(\cdot)=\partial_{t}(\cdot)+{\mathcal{U}}\partial(\cdot)/\partial
x=D^{\alpha}_{t}(\cdot)+{\mathcal{U}}^{\alpha}\partial(\cdot)/\partial
X^{\alpha}.$ (42)

This is the time derivative of a quantity measured by an observer moving with
the shock front.

Jump equations can be derived for conservation of mass, momentum, and energy
between any two points in a structured steady wave form [88, 89]. In the 1-D
Lagrangian description, let $\mathcal{D}^{\alpha}$ be a constant steady wave
speed, such that for a differentiable function $f(X^{\alpha},t)$, within the
steady waveform,

$\displaystyle
f(X^{\alpha},t)=f(X^{\alpha}-{\mathcal{D}}^{\alpha}t)=f(Y^{\alpha}),$ (43)
$\displaystyle\partial f/\partial X^{\alpha}={\rm d}f/{\rm d}Y^{\alpha},\quad
D^{\alpha}_{t}f={-\mathcal{D}}^{\alpha}{\rm d}f/{\rm d}Y^{\alpha}.$ (44)

Applying (44) to 1-D equations
$D^{\alpha}_{t}F^{\alpha}=\partial\upsilon^{\alpha}/\partial X^{\alpha}$ and

$\rho_{0}^{\alpha}D_{t}^{\alpha}\upsilon^{\alpha}=\partial t^{\alpha}/\partial
X^{\alpha}+\rho^{\alpha}_{0}b^{\alpha}+(\partial\chi^{\alpha}/\partial
X^{\alpha})h^{\alpha},$ (45)

which is the first of (16), gives

$\displaystyle{\rm d}\upsilon^{\alpha}/{\rm d}Y={-\mathcal{D}}^{\alpha}{\rm
d}F^{\alpha}/{\rm d}Y,$ (46) $\displaystyle{\rm d}t^{\alpha}/{\rm
d}Y=-\rho_{0}^{\alpha}{\mathcal{D}}^{\alpha}{\rm d}\upsilon^{\alpha}/{\rm
d}Y-\rho^{\alpha}_{0}b^{\alpha}-F^{\alpha}h^{\alpha}.$ (47)

Select two points at steady-wave coordinates $Y^{\pm}$, and define the jump in
a quantity between these points as in (23): $\llbracket
f(Y)\rrbracket=f(Y^{-})-f(Y^{+})$. Direct integration of (46) over
$Y^{-}\rightarrow Y^{+}$ and substitution into (47) gives

$\displaystyle\llbracket\upsilon^{\alpha}\rrbracket=-{\mathcal{D}}^{\alpha}\llbracket
F^{\alpha}\rrbracket,$ (48) $\displaystyle\llbracket
t^{\alpha}\rrbracket=\rho_{0}^{\alpha}({\mathcal{D}}^{\alpha})^{2}\llbracket
F^{\alpha}\rrbracket+\int_{-}^{+}(\rho^{\alpha}_{0}b^{\alpha}+F^{\alpha}h^{\alpha}){\rm
d}Y.$ (49)

Using the same procedure for 1-D continuum laws of energy conservation and
entropy production, (II.1) and (19),

$\displaystyle\rho_{0}^{\alpha}D_{t}^{\alpha}u^{\alpha}=t^{\alpha}\partial\upsilon^{\alpha}/\partial
X^{\alpha}+\partial(\\{z\\}^{\alpha}\\}\\{D_{t}^{\alpha}{\xi}^{\alpha}\\})/\partial
X^{\alpha}$ $\displaystyle\qquad\qquad\quad-\partial q^{\alpha}/\partial
X^{\alpha}+\rho_{0}^{\alpha}r^{\alpha}+(\partial\chi^{\alpha}/\partial
X^{\alpha})\epsilon^{\alpha},$ (50)
$\displaystyle\sum_{\alpha}[\rho_{0}^{\alpha}D^{\alpha}_{t}\eta^{\alpha}+\partial(q^{\alpha}/\theta^{\alpha})/\partial
X^{\alpha}-\rho_{0}^{\alpha}r^{\alpha}/\theta^{\alpha}$
$\displaystyle\qquad\qquad\qquad\qquad\quad+(\partial\chi^{\alpha}/\partial
X^{\alpha})c^{\alpha}\eta^{\alpha}]\geq 0,$ (51)

produces jumps between two points $Y^{\pm}$ in a steady wave:

$\displaystyle\rho^{\alpha}_{0}{\mathcal{D}}^{\alpha}\llbracket
u^{\alpha}+{\textstyle{\frac{1}{2}}}|{\upsilon}^{\alpha}|^{2}\rrbracket=-\llbracket{t}^{\alpha}{\upsilon}^{\alpha}+\\{{z}^{\alpha}\\}\\{D^{\alpha}_{t}{\xi}^{\alpha}\\}-q^{\alpha}\rrbracket$
$\displaystyle\qquad+\int_{-}^{+}\\{\rho_{0}^{\alpha}(r^{\alpha}+b^{\alpha}\upsilon^{\alpha})+F^{\alpha}(\epsilon^{\alpha}+h^{\alpha}\upsilon^{\alpha})\\}{\rm
d}Y,$ (52)
$\displaystyle\sum_{\alpha}(\rho^{\alpha}_{0}{\mathcal{D}}^{\alpha}\llbracket\eta^{\alpha}\rrbracket-\llbracket
q^{\alpha}/\theta^{\alpha}\rrbracket)$
$\displaystyle\qquad\qquad-\int_{-}^{+}(\rho^{\alpha}_{0}r^{\alpha}/\theta^{\alpha}-F^{\alpha}c^{\alpha}\eta^{\alpha}){\rm
d}Y\\}\geq 0.$ (53)

Relation (48) is identical to the mass conservation law for a singular surface
in the first of (40) when
${\mathcal{D}}^{\alpha}\rightarrow{\mathcal{U}}^{\alpha}$. Relation (49) is
identical to the momentum conservation law in the second of (40) when the
integral on the right of (49) vanishes (e.g., constant body force and null
drag). Relation (52) is identical to the energy conservation law in (II.2)
when its integral terms vanish, and (53) is identical to entropy inequality
(38) when its integral terms involving heat and mass supplies vanish. In such
cases, the Eulerian jump conditions in (28) and (29) can be recovered for a
steady wave form of Eulerian speed ${\mathcal{D}}$ where
$f(x,t)=f(x-{\mathcal{D}}t)$ when
${\mathcal{D}}=F^{\alpha}{\mathcal{D}}^{\alpha}+\upsilon^{\alpha}={\rm
constant}$.

###  II.3  Total expressions for mixture

The spatial mass density of the mixture $\rho$, mean velocity of the mixture
${\bm{\upsilon}}$ , and diffusion velocities ${\bm{\mu}}^{\alpha}$ are

$\displaystyle\rho=\sum_{\alpha}\rho^{\alpha},\quad{\bm{\upsilon}}=\frac{1}{\rho}\sum_{\alpha}\rho^{\alpha}{\bm{\upsilon}}^{\alpha},\quad{\bm{\mu}}^{\alpha}={\bm{\upsilon}}^{\alpha}-{\bm{\upsilon}}.$
(54)

Note $\sum_{\alpha}\rho^{\alpha}{\bm{\mu}}^{\alpha}={\bf 0}$. Define the
material time derivative of $\square$ with respect to the mixture as

$\dot{\square}=\partial_{t}(\square)+\nabla(\square)\cdot{\bm{\upsilon}}\,\Rightarrow\,D^{\alpha}_{t}(\square)=\dot{\square}+(\nabla\square)\cdot{\bm{\mu}}^{\alpha}.$
(55)

Summing (15) over $\alpha$ then gives the total mass balance for the mixture
at space-time location $({\bf x},t)$:

$\dot{\rho}+\rho\nabla\cdot{\bm{\upsilon}}=\hat{C},\qquad\hat{C}=\sum_{\alpha}c^{\alpha}-\rho\,\partial_{t}\ln\sqrt{g}.$
(56)

In classical mixture theory [43], rates $c^{\alpha}$ account for exchange of
mass between phases but not mass production from extrinsic sources, thus
locally summing to zero. To accommodate mass production from biologic (e.g.,
cellular growth) processes even in a single-phase material,
$\sum_{\alpha}c^{\alpha}({\bf x},t)$ could be nonzero [61, 90].

Define, respectively, the total Cauchy stress tensor, total body force vector,
total internal energy density, total entropy density, total heat supply, total
heat flux, total internal state vector, and total conjugate force vector:

$\displaystyle{\bm{\sigma}}=\sum_{\alpha}({\bm{\sigma}}^{\alpha}-\rho^{\alpha}{\bm{\mu}}^{\alpha}\otimes{\bm{\mu}}^{\alpha}),\quad{\bf
b}=\frac{1}{\rho}\sum_{\alpha}\rho^{\alpha}{\bf b}^{\alpha},$ (57)
$\displaystyle
u=\frac{1}{\rho}\sum_{\alpha}\rho^{\alpha}\left[{u}^{\alpha}+{\textstyle{\frac{1}{2}}}|{\bm{\mu}}^{\alpha}|^{2}\right],$
(58)
$\displaystyle\eta=\frac{1}{\rho}\sum_{\alpha}\rho^{\alpha}{\eta}^{\alpha},\qquad
r=\frac{1}{\rho}\sum_{\alpha}\rho^{\alpha}{r}^{\alpha},$ (59)
$\displaystyle{\bf q}=\sum_{\alpha}({\bf
q}^{\alpha}-{\bm{\sigma}}^{\alpha}\cdot{\bm{\mu}}^{\alpha}+\rho^{\alpha}u^{\alpha}{\bm{\mu}}^{\alpha}+\frac{\rho^{\alpha}}{2}|{\bm{\mu}}^{\alpha}|^{2}{\bm{\mu}}^{\alpha}),$
(60)
$\displaystyle\\{{\bm{\xi}}\\}=(\\{{\bm{\xi}}\\}^{1},\ldots,\\{{\bm{\xi}}\\}^{N}),\quad\\{{\bm{\zeta}}\\}=(\\{{\bm{\zeta}}\\}^{1},\ldots,\\{{\bm{\zeta}}\\}^{N}).$
(61)

Define $\hat{c}^{\alpha}=c^{\alpha}-\rho^{\alpha}\partial_{t}\ln\sqrt{g}$ and
impose the constraints

$\displaystyle\sum_{\alpha}\hat{c}^{\alpha}=0,\qquad\sum_{\alpha}({\bf
h}^{\alpha}+\hat{c}^{\alpha}{\bm{\mu}^{\alpha}})={\bf 0},$ (62)
$\displaystyle\sum_{\alpha}[\epsilon^{\alpha}+{\bf
h}^{\alpha}\cdot{\bm{\mu}}^{\alpha}+\hat{c}^{\alpha}(u^{\alpha}+{\textstyle{\frac{1}{2}}}|{\bm{\mu}}^{\alpha}|^{2})]=0.$
(63)

Ensuring that the second of (62) holds, net energy rates from biologic growth
processes can be included in $r^{\alpha}$. Given (54)–(61) with constraints
(62) and (63), the total balances of mass, linear momentum, angular momentum,
and energy and the dissipation inequality are obtained by accumulating (15),
(16), (II.1), and (19) over $\alpha=1,\ldots,N$, following steps detailed by
Bowen [43], with the addition of work contributions from (61) whose individual
entries are defined as orthogonal:

$\displaystyle\dot{\rho}+\rho\nabla\cdot{\bm{\upsilon}}=0,\quad\nabla\cdot{\bm{\sigma}}+\rho{\bf
b}=\rho\dot{{\bm{\upsilon}}},\quad{\bm{\sigma}}={\bm{\sigma}}^{\mathsf{T}},$
(64)
$\displaystyle\rho\dot{u}={\bm{\sigma}}:\nabla{\bm{\upsilon}}+\nabla\cdot(\\{{\bm{\zeta}}\\}\cdot\\{{\bm{\xi}}^{\prime}\\})-\nabla\cdot{\bf
q}+\rho r$ $\displaystyle\qquad+\sum_{\alpha}\rho^{\alpha}{\bf
b}^{\alpha}\cdot{\bm{\mu}}^{\alpha},\quad\\{{\bm{\xi}}^{\prime}\\}=\\{\dot{\bm{\xi}}\\}+\sum_{\alpha}\\{\nabla{\bm{\xi}}^{\alpha}\cdot{\bm{\mu}}^{\alpha}\\},$
(65) $\displaystyle\rho\dot{\eta}+\nabla\cdot\sum_{\alpha}(\frac{{\bf
q}^{\alpha}}{\theta^{\alpha}}+\rho^{\alpha}\eta^{\alpha}{\bm{\mu}}^{\alpha})-\sum_{\alpha}\frac{\rho^{\alpha}r^{\alpha}}{\theta^{\alpha}}\geq
0.$ (66)

For the particular case when $\theta^{\alpha}=\theta$ is uniform among
constituents at $({\bf x},t)$, then (66) simplifies to

$\rho\dot{\eta}+\nabla\cdot(\hat{\bf q}/\theta)-\rho r/\theta\geq
0,\quad\hat{\bf q}=\sum_{\alpha}({\bf
q}^{\alpha}+\theta\rho^{\alpha}\eta^{\alpha}{\bm{\mu}}^{\alpha}).$ (67)

Jump conditions across singular surfaces can be derived from (64), (II.3), and
(67) using the methods already discussed for constituent $\alpha$. Let
$\mathcal{U}$ now be the Eulerian shock speed for the mixture as a whole.
Analogs of (24)–(27) are, with an obvious change of notation,

$\displaystyle\llbracket\rho(\upsilon_{n}-\mathcal{U})\rrbracket=0,\qquad\llbracket\rho{\bm{\upsilon}}(\upsilon_{n}-\mathcal{U})\rrbracket=\llbracket{\bf
t}\rrbracket,$ (68)
$\displaystyle\llbracket\rho(u+{\textstyle{\frac{1}{2}}}|{\bm{\upsilon}}|^{2})(\upsilon_{n}-\mathcal{U})\rrbracket=\llbracket{\bf
t}\cdot{\bm{\upsilon}}+\\{{\bf
z}\\}\cdot\\{{\bm{\xi}}^{\prime}\\}-\hat{q}_{n}\rrbracket,$ (69)
$\displaystyle\llbracket\rho\eta(\mathcal{U}-\upsilon_{n})-\hat{q}_{n}/\theta\rrbracket\geq
0.$ (70)

Eulerian equations for the 1-D case follow trivially. Lagrangian equations for
the 1-D case are obtained by defining Lagrangian speed $\mathcal{U}_{0}$,
deformation mapping $F$, and space-time continuous reference mass density
$\rho_{0}$ to obey

$\mathcal{U}_{0}=F^{-1}(\mathcal{U}-\upsilon),\qquad F=\rho_{0}/\rho.$ (71)

Given $\rho_{0}$ and $\mathcal{U}$, $F$ and $\mathcal{U}_{0}$ are well-
defined. Analogs of (35)–(38) are then derived as, with $t_{k}\rightarrow
t_{1}=\sigma$,

$\displaystyle\rho_{0}{\mathcal{U}}_{0}\llbracket
1/\rho\rrbracket=-\llbracket\upsilon\rrbracket,\qquad\rho_{0}{\mathcal{U}}_{0}\llbracket\upsilon\rrbracket=-\llbracket{\sigma}\rrbracket,$
(72) $\displaystyle\rho_{0}{\mathcal{U}}_{0}\llbracket
u+{\textstyle{\frac{1}{2}}}{\upsilon}^{2}\rrbracket=-\llbracket{\sigma}{\upsilon}+\\{{z}\\}\\{{\xi}^{\prime}\\}-\hat{q}\rrbracket,$
(73)
$\displaystyle\rho_{0}{\mathcal{U}}_{0}\llbracket\eta\rrbracket-\llbracket\hat{q}/\theta\rrbracket\geq
0.$ (74)

Equations fully analogous to (39)–(41) can be derived summarily with the
obvious substitutions.

Presuming invertibility and integrability of the second of (71), deformation
gradient $F=\partial x/\partial X$ and Lagrangian coordinate $X$ for the
mixture exist whereby

$\displaystyle\dot{F}=(\partial\upsilon/\partial x)F=\partial\upsilon/\partial
X,$ (75) $\displaystyle
X(x,t)=\chi^{-1}(x,t)=\int_{x_{0}}^{x}F^{-1}(\tilde{x},t){\rm d}\tilde{x},$
(76)

then the analog of the displacement derivative (42) for the mixture is

$\delta_{t}(\square)=\partial_{t}(\square)+{\mathcal{U}}\partial(\square)/\partial
x=(\dot{\square})+{\mathcal{U}_{0}}\partial(\square)/\partial X.$ (77)

Given (64), (II.3), (67), (75), and (76), the treatment of steady Lagrangian
waves can be applied to the mixture as a whole. Let
$f(X,t)=f(X-{\mathcal{D}}t)=f(Y)$ in a steady wave, with $\mathcal{D}$ the
constant speed. Using relations akin to (43) and (44), the local compatibility
and linear momentum equations (75) and (64) (in 1-D) are

$\displaystyle{\rm d}\upsilon/{\rm d}Y={-\mathcal{D}}{\rm d}F^{\alpha}/{\rm
d}Y,$ (78) $\displaystyle{\rm d}\sigma/{\rm d}Y=-\rho_{0}{\mathcal{D}}{\rm
d}\upsilon/{\rm d}Y-\rho_{0}b.$ (79)

Integrating from $Y^{+}\rightarrow Y^{-}$ gives conditions like (72):

$\displaystyle\llbracket\upsilon\rrbracket=-{\mathcal{D}}\llbracket
F\rrbracket=-{\mathcal{D}}\rho_{0}\llbracket 1/\rho\rrbracket,$ (80)
$\displaystyle\llbracket\sigma\rrbracket=-\rho_{0}{\mathcal{D}}\llbracket\upsilon\rrbracket+\int_{-}^{+}\rho_{0}b\,{\rm
d}Y.$ (81)

Energy conservation and entropy production in (II.3) and (66) can be addressed
similarly, omitted here for brevity.

Returning to the general 3-D case, diffusion problems are often analyzed using
dimensionless measures of local amounts of each constituent. Recall
$\rho^{\alpha}$ is the local mass of $\alpha$ per unit total spatial volume of
mixture. The spatial volume fraction $n^{\alpha}$ is the ratio of volume
occupied by $\alpha$ to that of the mixture, while the real mass density
$\rho^{\alpha}_{\rm R}$ is the local mass of $\alpha$ per unit spatial volume
occupied by $\alpha$ (i.e., $\rho^{\alpha}_{\rm R}$ is mass density of the
isolated, fully dense constituent). The spatial mass concentration
$m^{\alpha}$ is the mass of $\alpha$ per total mass of the mixture. Equations
are

$\displaystyle n^{\alpha}({\bf x},t)=\frac{\rho^{\alpha}({\bf
x},t)}{\rho^{\alpha}_{\rm R}({\bf x},t)},\qquad\sum_{\alpha}n^{\alpha}=1,$
(82) $\displaystyle m^{\alpha}({\bf x},t)=\frac{\rho^{\alpha}({\bf
x},t)}{\rho({\bf x},t)},\qquad\sum_{\alpha}m^{\alpha}=1.$ (83)

At a reference $t=t_{0}$ when all particles occupy positions ${\bf
X^{\alpha}}$, reference volume and mass fractions are $n^{\alpha}_{0}({\bf
X}^{\alpha})=\rho^{\alpha}_{0}({\bf X}^{\alpha})/\rho^{\alpha}_{{\rm R}0}({\bf
X}^{\alpha})$ and $m^{\alpha}_{0}({\bf X}^{\alpha})=\rho^{\alpha}_{0}({\bf
X}^{\alpha})/\rho_{0}({\bf X}^{\alpha})$.

##  III  Constitutive Theory

Thermodynamic identities are derived in Sec. III.1 by appealing to the local
balance of energy and entropy inequality. Pragmatic and thermodynamically
admissible energy functions and kinetic relations are posited in respective
Secs. III.2 and III.3.

###  III.1 Thermodynamics

Helmholtz free energy per unit mass and entropy density are of the following
functional forms, each depending only on state variables for its particular
constituent $\alpha$ and not others $\beta$ with $\beta\neq\alpha$:

$\displaystyle\psi^{\alpha}=\psi^{\alpha}({\bf
F}^{\alpha},\theta^{\alpha},\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla^{\alpha}_{0}{\bm{\xi}}^{\alpha}\\},{\bf
X}^{\alpha}),$ (84) $\displaystyle\eta^{\alpha}=\eta^{\alpha}({\bf
F}^{\alpha},\theta^{\alpha},\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla^{\alpha}_{0}{\bm{\xi}}^{\alpha}\\},{\bf
X}^{\alpha}).$ (85)

Partial stress consists of elastic $\bar{\bm{\sigma}}^{\alpha}$ and viscous
$\hat{\bm{\sigma}}^{\alpha}$ parts:

$\begin{split}{\bm{\sigma}}^{\alpha}&=\bar{\bm{\sigma}}^{\alpha}({\bf
F}^{\alpha},\theta^{\alpha},\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla^{\alpha}_{0}{\bm{\xi}}^{\alpha}\\},{\bf
X}^{\alpha})\\\ &\qquad+\hat{\bm{\sigma}}^{\alpha}({\bf
F}^{\alpha},\theta^{\alpha},\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla^{\alpha}_{0}{\bm{\xi}}^{\alpha}\\},D^{\alpha}_{t}{\bf
F}^{\alpha},{\bf X}^{\alpha}).\end{split}$ (86)

Arguments in (2) and (3) are linked at time $t$ via [63, 59]

$\\{{\bm{\Xi}}^{\alpha}\\}(t)=\\{{\bm{\Xi}}^{\alpha}(\\{\bm{\xi}^{\alpha}\\},{\bf
X}^{\alpha},{\bf x})\\}(t).$ (87)

Thus, all dependence of response functions on metrics $({\bf
g},\bf{G}^{\alpha})$ and states $\\{{\bm{\Xi}}^{\alpha}\\}$ is implicitly
included via arguments $\\{{\bm{\xi}}^{\alpha}\\}$. Kinetic equations for heat
flux, internal state, and interphase mass, momentum, and energy exchange are
more general, allowing for dependence on all constituents $\beta=1,\ldots,N$
including $\alpha=\beta$ and $\alpha\neq\beta$:

$\displaystyle{\bf q}^{\alpha}={\bf q}^{\alpha}({\bf F}^{\beta},\nabla{\bf
F}^{\beta},\theta^{\beta},\nabla\theta^{\beta},\bm{\xi}^{\beta},\nabla^{\beta}_{0}{\bm{\xi}}^{\beta},{\bm{\upsilon}}^{\beta},{\rho}^{\beta}),$
(88) $\displaystyle{\bf h}^{\alpha}={\bf h}^{\alpha}({\bf
F}^{\beta},\nabla{\bf
F}^{\beta},\theta^{\beta},\nabla\theta^{\beta},\bm{\xi}^{\beta},\nabla^{\beta}_{0}{\bm{\xi}}^{\beta},{\bm{\upsilon}}^{\beta},{\rho}^{\beta}),$
(89) $\displaystyle{{c}}^{\alpha}={c}^{\alpha}({\bf F}^{\beta},\nabla{\bf
F}^{\beta},\theta^{\beta},\nabla\theta^{\beta},\bm{\xi}^{\beta},\nabla^{\beta}_{0}{\bm{\xi}}^{\beta},{\bm{\upsilon}}^{\beta},{\rho}^{\beta}),$
(90) $\displaystyle{{\epsilon}}^{\alpha}={\epsilon}^{\alpha}({\bf
F}^{\beta},\nabla{\bf
F}^{\beta},\theta^{\beta},\nabla\theta^{\beta},\bm{\xi}^{\beta},\nabla^{\beta}_{0}{\bm{\xi}}^{\beta},{\bm{\upsilon}}^{\beta},{\rho}^{\beta}),$
(91)
$\displaystyle{D^{\alpha}_{t}\\{{\bm{\xi}}^{\alpha}\\}}={D^{\alpha}_{t}\\{{\bm{\xi}}^{\alpha}\\}}({\bf
F}^{\beta},\nabla{\bf F}^{\beta},\theta^{\beta},\nabla\theta^{\beta},\cdots$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\cdots\bm{\xi}^{\beta},\nabla^{\beta}_{0}{\bm{\xi}}^{\beta},{\bm{\upsilon}}^{\beta},{\rho}^{\beta}).$
(92)

Notation $\\{\cdot\\}$ on ${\bm{\xi}}^{\alpha}$ and admissible explicit
dependence on ${\bf X}^{\alpha}$ for heterogeneous phases $\alpha$ are omitted
in arguments of (88)–(92) for brevity. Particular forms of (88)–(91) must
satisfy principles of spatial invariance for objective spatial vectors ${\bf
q}^{\alpha}$ and ${\bf h}^{\alpha}$ and scalars $c^{\alpha}$ and
$\epsilon^{\alpha}$. Evolution equations (92) must also be objective.
Invariance under rigid translation of the mixture as a whole necessitates that
dependence on velocities ${\bm{\upsilon}}^{\alpha}$ is at most only on the
$N-1$ velocity differences
${\bm{\upsilon}}^{1}-{\bm{\upsilon}}^{N},\ldots,{\bm{\upsilon}}^{N-1}-{\bm{\upsilon}}^{N}$.
See Ref. [43]. Diffusion velocities
${\bm{\mu}}^{\alpha}={\bm{\mu}}^{\alpha}({\bm{\upsilon}}^{\beta},\rho^{\beta})$
of (54) fulfill this requirement. Spatial invariance of (84) and (85) is
obtained via dependence on ${\bf F}^{\alpha}$ through symmetric deformation
tensor ${\bf C}^{\alpha}$:

$\displaystyle{\bf C}^{\alpha}=({\bf F}^{\alpha})^{\mathsf{T}}{\bf
F}^{\alpha},\,(C^{\alpha})^{K}_{J}=(G^{\alpha})^{KI}(F^{\alpha})^{i}_{I}g_{ij}(F^{\alpha})^{j}_{J};$
(93) $\displaystyle\partial\psi^{\alpha}/\partial{\bf F}^{\alpha}=2{\bf
F}^{\alpha}\partial\psi^{\alpha}/\partial{\bf C}^{\alpha},\quad
J^{\alpha}=\sqrt{\det{\bf C}^{\alpha}};\quad$ (94) $\displaystyle
D^{\alpha}_{t}{\bf C}^{\alpha}=2({\bf F}^{\alpha})^{\mathsf{T}}{\bf
d}^{\alpha}{\bf F}^{\alpha},\quad 2{\bf d}^{\alpha}={\bf l}^{\alpha}+({\bf
l}^{\alpha})^{\mathsf{T}}.$ (95)

The spatial deformation rate is ${\bf d}^{\alpha}$, and $D^{\alpha}_{t}{\bf
C}^{\alpha}$ is taken with $(G^{\alpha})^{IJ}$ and $g_{ij}$ fixed with respect
to $t$ in (95).

Expanding $D^{\alpha}_{t}\psi$ using the chain rule on (84) and inserting the
result and (86) into (II.1) gives

$\displaystyle\sum_{\alpha}\frac{1}{\theta^{\alpha}}\\{[\bar{\bm{\sigma}}^{\alpha}({\bf
F}^{\alpha})^{-\mathsf{T}}-2\rho^{\alpha}{\bf
F}^{\alpha}\frac{\partial\psi^{\alpha}}{\partial{\bf
C}^{\alpha}}]:D_{t}^{\alpha}{\bf F}^{\alpha}$
$\displaystyle-\rho^{\alpha}[\eta^{\alpha}+\partial\psi^{\alpha}/\partial\theta^{\alpha}]D_{t}^{\alpha}\theta^{\alpha}$
$\displaystyle+[\\{({\bf
F}^{\alpha})^{-1}{\bm{\zeta}}^{\alpha}\\}-\rho^{\alpha}\frac{\partial\psi^{\alpha}}{\partial\\{\nabla_{0}{\bm{\xi}}^{\alpha}\\}}]:\\{D^{\alpha}_{t}(\nabla_{0}{\bm{\xi}}^{\alpha})\\}$
$\displaystyle+[\\{({\bf
F}^{\alpha})^{-1}:\nabla_{0}^{\alpha}{\bm{\zeta}}^{\alpha}\\}-\rho^{\alpha}\frac{\partial\psi^{\alpha}}{\partial\\{{\bm{\xi}}^{\alpha}\\}}]\cdot\\{D^{\alpha}_{t}{\bm{\xi}}^{\alpha}\\}$
$\displaystyle+\hat{\bm{\sigma}}^{\alpha}:{\bf d}^{\alpha}-{({\bf
q^{\alpha}}\cdot\nabla\theta^{\alpha})}/{\theta^{\alpha}}+\epsilon^{\alpha}+c^{\alpha}\theta^{\alpha}\eta^{\alpha}\\}\geq
0.$ (96)

Identities from (6) and (8) have been used to obtain

$\begin{split}\nabla\cdot(\\{{\bm{\zeta}}^{\alpha}\\}&\cdot\\{D_{t}^{\alpha}{\bm{\xi}}^{\alpha}\\})=(\nabla\cdot\\{{\bm{\zeta}}^{\alpha}\\})\cdot\\{D_{t}^{\alpha}{\bm{\xi}}^{\alpha}\\}\\\
&\qquad\qquad+\\{({\bf
F}^{\alpha})^{-1}{\bm{\zeta}}^{\alpha}\\}:\\{D^{\alpha}_{t}(\nabla^{\alpha}_{0}{\bm{\xi}}^{\alpha})\\}.\end{split}$
(97)

From standard arguments [91, 92, 80, 38] and (84)–(92), the first three sets
of terms in (III.1) should vanish for admissibility under general
thermodynamic processes, leading to the following constitutive equalities:

$\displaystyle\bar{\bm{\sigma}}^{\alpha}=2\rho^{\alpha}{\bf
F}^{\alpha}\frac{\partial\psi^{\alpha}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}},\qquad\eta^{\alpha}=-\frac{\partial\psi^{\alpha}}{\partial\theta^{\alpha}},$
(98) $\displaystyle\\{{\bm{\zeta}}^{\alpha}\\}=\rho^{\alpha}\\{{\bf
F}^{\alpha}{\partial\psi^{\alpha}}/{\partial\nabla_{0}^{\alpha}{\bm{\xi}}^{\alpha}\\}}=\rho^{\alpha}\frac{\partial\psi^{\alpha}}{\partial\\{\nabla{\bm{\xi}}^{\alpha}\\}},$
(99)
$\displaystyle\\{{\bm{\pi}}^{\alpha}\\}=\rho^{\alpha}\partial\psi^{\alpha}/\partial\\{{\bm{\xi}}^{\alpha}\\},$
(100)

where (100) defines a conjugate force to internal state variables or order
parameters. Then (III.1) reduces to

$\displaystyle\sum_{\alpha}\frac{1}{\theta^{\alpha}}[(\\{\nabla\cdot{\bm{\zeta}}^{\alpha}\\}-\\{{\bm{\pi}}^{\alpha}\\})\cdot\\{D^{\alpha}_{t}{\bm{\xi}}^{\alpha}\\}$
$\displaystyle+\hat{\bm{\sigma}}^{\alpha}:{\bf d}^{\alpha}-{({\bf
q}^{\alpha}\cdot\nabla\theta^{\alpha})}/{\theta^{\alpha}}+\epsilon^{\alpha}+c^{\alpha}\theta^{\alpha}\eta^{\alpha}]\geq
0.$ (101)

Applying the Legendre transformation from (20) with

$\displaystyle u^{\alpha}=u^{\alpha}({\bf
F}^{\alpha},\eta^{\alpha},\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla_{0}^{\alpha}{\bm{\xi}}^{\alpha}\\},\bf{X}^{\alpha}),$
(102) $\displaystyle\theta^{\alpha}=\theta^{\alpha}({\bf
F}^{\alpha},\eta^{\alpha},\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla_{0}^{\alpha}{\bm{\xi}}^{\alpha}\\},\bf{X}^{\alpha}),$
(103)

in conjunction with (98) and (99), gives

$\displaystyle\bar{\bm{\sigma}}^{\alpha}=2\rho^{\alpha}{\bf
F}^{\alpha}\frac{\partial u^{\alpha}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}},\qquad\theta^{\alpha}=\frac{\partial
u^{\alpha}}{\partial\eta^{\alpha}},$ (104)
$\displaystyle\\{{\bm{\pi}}^{\alpha}\\}=\rho^{\alpha}\frac{\partial
u^{\alpha}}{\partial\\{{\bm{\xi}}^{\alpha}\\}},\quad\\{{\bm{\zeta}}^{\alpha}\\}=\rho^{\alpha}\frac{\partial
u^{\alpha}}{\partial\\{\nabla{\bm{\xi}}^{\alpha}\\}}.$ (105)

Define specific heat per unit mass at constant strain $c_{\epsilon}^{\alpha}$,
thermal stress coefficients ${\bm{\beta}}^{\alpha}$, and Grüneisen tensor
${\bm{\gamma}}^{\alpha}$:

$\displaystyle
c_{\epsilon}^{\alpha}=\theta^{\alpha}\partial\eta^{\alpha}/\partial\theta^{\alpha}=-\theta^{\alpha}\,{\partial^{2}\psi^{\alpha}}/{\partial(\theta^{\alpha})^{2}},$
(106)
$\displaystyle{\bm{\beta}}^{\alpha}=\rho^{\alpha}c^{\alpha}_{\epsilon}{\bm{\gamma}}^{\alpha}=-2\rho^{\alpha}\,{\partial^{2}\psi^{\alpha}}/{\partial\theta^{\alpha}\partial{\bf
C}^{\alpha}}.$ (107)

Define the intrinsic dissipation for constituent $\alpha$:

${\mathfrak{D}}^{\alpha}=(\\{\nabla\cdot{\bm{\zeta}}^{\alpha}\\}-\\{{\bm{\pi}}^{\alpha}\\})\cdot\\{D^{\alpha}_{t}{\bm{\xi}}^{\alpha}\\}+\hat{\bm{\sigma}}^{\alpha}:{\bf
d}^{\alpha}.$ (108)

Expand the rate of $\eta^{\alpha}$ using (85), (98), (106), and (107):

$\begin{split}\rho^{\alpha}\theta^{\alpha}D_{t}^{\alpha}\eta^{\alpha}&=\rho^{\alpha}c_{\epsilon}^{\alpha}D^{\alpha}_{t}\theta^{\alpha}+{\textstyle{\frac{1}{2}}}\theta^{\alpha}{\bm{\beta}}^{\alpha}:D_{t}^{\alpha}{\bf
C}^{\alpha}\\\
&-\rho^{\alpha}\theta^{\alpha}[(\partial^{2}\psi/\partial\theta^{\alpha}\partial\\{{\bm{\xi}}^{\alpha}\\})\cdot\\{D_{t}^{\alpha}{\bm{\xi}}^{\alpha}\\}\\\
&+(\partial^{2}\psi/\partial\theta^{\alpha}\partial\\{\nabla{\bm{\xi}}^{\alpha}\\}):\\{\nabla(D_{t}^{\alpha}{\bm{\xi}}^{\alpha})\\}].\end{split}$
(109)

From (II.1), time differentiation of (20), and (98) and (99):

$\begin{split}\rho^{\alpha}\theta^{\alpha}D_{t}^{\alpha}\eta^{\alpha}={\mathfrak{D}}^{\alpha}-\nabla\cdot{\bf
q}^{\alpha}+\rho^{\alpha}r^{\alpha}+\epsilon^{\alpha}.\end{split}$ (110)

Temperature rates then are, combining (109) and (110),

$\begin{split}\rho^{\alpha}c_{\epsilon}^{\alpha}D^{\alpha}_{t}\theta^{\alpha}=&\,{\mathfrak{D}}^{\alpha}-{\textstyle{\frac{1}{2}}}\theta^{\alpha}{\bm{\beta}}^{\alpha}:D_{t}^{\alpha}{\bf
C}^{\alpha}\\\
&+\rho^{\alpha}\theta^{\alpha}[(\partial^{2}\psi/\partial\theta^{\alpha}\partial\\{{\bm{\xi}}^{\alpha}\\})\cdot\\{D_{t}^{\alpha}{\bm{\xi}}^{\alpha}\\}\\\
&+(\partial^{2}\psi/\partial\theta^{\alpha}\partial\\{\nabla{\bm{\xi}}^{\alpha}\\}):\\{\nabla(D_{t}^{\alpha}{\bm{\xi}}^{\alpha})\\}]\\\
&-\nabla\cdot{\bf
q}^{\alpha}+\rho^{\alpha}r^{\alpha}+\epsilon^{\alpha}.\end{split}$ (111)

Decompositions of metrics of (2) and (3) into symmetric position-dependent
(i.e., classical) and dimensionless, invertible, space-time dependent parts
are [59, 61]

$\displaystyle{\bf g}({\bf x},t)=\bar{\bf g}({\bf x})\,\hat{\bf
g}(\\{{\bm{\xi}}^{\alpha}({\bf x},t)\\}),$ (112) $\displaystyle{\bf
G}^{\alpha}({\bf X}^{\alpha},t)=\bar{\bf G}^{\alpha}({\bf
X}^{\alpha})\,\hat{\bf G}^{\alpha}(\\{{\bm{\Xi}}^{\alpha}({\bf
X}^{\alpha},t)\\}).$ (113)

A deformation $\bar{\bf C}^{\alpha}$ and Jacobian $\bar{J}^{\alpha}$ based on
$(\bar{\bf g},\bar{\bf G}^{\alpha})$ are

$\displaystyle(\bar{C}^{\alpha})^{K}_{J}=\bar{G}^{KI}(F^{\alpha})^{i}_{I}\bar{g}_{ij}(F^{\alpha})^{j}_{J}({\bf
G}^{\alpha})_{K}\otimes({\bf G}^{\alpha})^{J},$ (114)
$\displaystyle\bar{J}^{\alpha}=\sqrt{\det\bar{\bf
C}^{\alpha}}=J^{\alpha}\sqrt{\hat{G}^{\alpha}/\hat{g}},$ (115)

with dimensionless $\hat{g}=\det\hat{\bf g}$ and
$\hat{G}^{\alpha}=\det\hat{\bf G}^{\alpha}$. Alternative constitutive
equations [59], also energetically objective, are obtained by positing
dependence of $\psi^{\alpha}$, $\eta^{\alpha}$, $u^{\alpha}$, and
$\theta^{\alpha}$ through $\bar{\bf C}^{\alpha}(\bf{F}^{\alpha})$ rather than
${\bf C}^{\alpha}$, whereby

$\begin{split}&\partial\psi^{\alpha}/\partial{\bf F}^{\alpha}=2{\bf
F}^{\alpha}\partial\psi^{\alpha}/\partial\bar{\bf C}^{\alpha},\\\
&\quad\leftrightarrow\,\partial\psi^{\alpha}/\partial(F^{\alpha})^{i}_{J}=2\bar{g}_{ik}(F^{\alpha})^{k}_{L}\partial\psi^{\alpha}/\partial(\bar{C}^{\alpha})_{JL}.\end{split}$
(116)

Derivations in (III.1)–(111) continue to apply for $\psi^{\alpha}(\bar{\bf
C}^{\alpha},\cdot)$, $\eta^{\alpha}(\bar{\bf C}^{\alpha},\cdot)$, etc. with
several changes manifested by (116):

$\displaystyle(\bar{\sigma}^{\alpha})^{j}_{i}=2{\rho}^{\alpha}\bar{g}_{ik}(F^{\alpha})^{k}_{L}(F^{\alpha})^{j}_{J}\partial\psi^{\alpha}/\partial(\bar{C}^{\alpha})_{JL}$
$\displaystyle\qquad\,\,\,=2{\rho}^{\alpha}\bar{g}_{ik}(F^{\alpha})^{k}_{L}(F^{\alpha})^{j}_{J}\partial
u^{\alpha}/\partial(\bar{C}^{\alpha})_{JL},$ (117)
$\displaystyle\bar{\bm{\beta}}^{\alpha}=\rho^{\alpha}c^{\alpha}_{\epsilon}\bar{\bm{\gamma}}^{\alpha}=-2\rho^{\alpha}\,{\partial^{2}\psi^{\alpha}}/{\partial\theta^{\alpha}\partial\bar{\bf
C}^{\alpha}},$ (118)
$\displaystyle\rho^{\alpha}c_{\epsilon}^{\alpha}D^{\alpha}_{t}\theta^{\alpha}=\,{\mathfrak{D}}^{\alpha}-{\textstyle{\frac{1}{2}}}\theta^{\alpha}\bar{\bm{\beta}}^{\alpha}:D_{t}^{\alpha}\bar{\bf
C}^{\alpha}$
$\displaystyle\qquad\qquad+\rho^{\alpha}\theta^{\alpha}[(\partial^{2}\psi/\partial\theta^{\alpha}\partial\\{{\bm{\xi}}^{\alpha}\\})\cdot\\{D_{t}^{\alpha}{\bm{\xi}}^{\alpha}\\}$
$\displaystyle\qquad\qquad+(\partial^{2}\psi/\partial\theta^{\alpha}\partial\\{\nabla{\bm{\xi}}^{\alpha}\\}):\\{\nabla(D_{t}^{\alpha}{\bm{\xi}}^{\alpha})\\}]$
$\displaystyle\qquad\qquad-\nabla\cdot{\bf
q}^{\alpha}+\rho^{\alpha}r^{\alpha}+\epsilon^{\alpha}.$ (119)

Mixed-variant $\bar{\bm{\sigma}}^{\alpha}$ in (III.1) excludes $\hat{\bf
g}({\bf x},t)$. Contravariant stress defined as
$(\bar{\sigma}^{\alpha})^{ij}=\frac{1}{2}[g^{ik}(\bar{\sigma}^{\alpha})^{j}_{k}+g^{jk}(\bar{\sigma}^{\alpha})^{i}_{k}$]
or $(\bar{\sigma}^{\alpha})^{ij}=\bar{g}^{ik}(\bar{\sigma}^{\alpha})^{j}_{k}$
must be symmetric. The former depends on $\hat{\bf g}({\bf x},t)$ implicitly
from $g^{ik}$. This choice presumes, a priori, that skew contributions from
(III.1) perform no work in the energy balance so can thus be redefined as
zero. The latter prescription either redefines raising/lowering indices on
Cauchy stress or presumes that ${\bf d}^{\alpha}$ is defined in covariant form
by lowering of ${\bf l}^{\alpha}$ with $\bar{g}_{ij}$, rather than the typical
$g_{ij}$ prior to symmetrization.

###  III.2 Energy functions

Internal state variables $\\{{\bm{\xi}}^{\alpha}\\}$ consist of three sets:
configurational variables associated with viscoelastic processes
$\\{{\bm{\Gamma}}^{\alpha}\\}$, damage variables associated with degradation
processes $\\{{\bf D}^{\alpha}\\}$, and electrochemical activation (e.g.,
muscle contraction) variables $\\{{\bm{\Delta}}^{\alpha}\\}$ [23, 19, 26, 59]:

$\\{{\bm{\xi}}^{\alpha}\\}({\bf x},t)=(\\{{\bm{\Gamma}}^{\alpha}\\},\\{{\bf
D}^{\alpha}\\},\\{{\bm{\Delta}}^{\alpha}\\})({\bf x},t).$ (120)

In the present work, as in phase-field and related theories [26, 59], free and
internal energy functions can depend on spatial gradients of damage variables,
which are viewed as order parameters, but not on spatial gradients of
viscoelastic and tissue activation variables. Denote by
$\varsigma^{\alpha}_{\rm V}$ and $\varsigma^{\alpha}_{\rm S}$ degradation
functions associated with loss of strength due to changes in bulk and
deviatoric strain energies, respectively. These scalar functions obey

$\displaystyle\varsigma^{\alpha}_{\rm V}=\varsigma^{\alpha}_{\rm V}(\\{{\bf
D}^{\alpha}\\},{\bf C}^{\alpha})\in[0,1],\,\,\varsigma^{\alpha}_{\rm
S}=\varsigma^{\alpha}_{\rm S}(\\{{\bf D}^{\alpha}\\})\in[0,1],$ (121)
$\displaystyle\partial\varsigma^{\alpha}_{\rm V}/\partial{\bf
C}^{\alpha}(\\{{\bf D}^{\alpha}\\},{\bf C}^{\alpha})={\bf 0}\,\forall\,{\bf
C}^{\alpha}\neq{\bf 1}.$ (122)

A degradation operator for fibrous energy contributions with similar
properties is $\varsigma^{\alpha}_{\rm F}(\\{{\bf D}^{\alpha}\\})$. Let
$\Psi^{\alpha}=\rho^{\alpha}_{{\rm R}0}\psi^{\alpha}$ and
$U^{\alpha}=\rho^{\alpha}_{{\rm R}0}u^{\alpha}$ be free and internal energies
per unit reference volume of individual phases. Pragmatic functional forms
consist of the following sums:

$\displaystyle\Psi^{\alpha}({\bf C}^{\alpha},$
$\displaystyle\,\theta^{\alpha},\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla_{0}^{\alpha}{\bf
D}^{\alpha}\\})=$ $\displaystyle\quad\varsigma^{\alpha}_{\rm V}(\\{{\bf
D}^{\alpha}\\},{\bf C}^{\alpha})\Psi^{\alpha}_{{\rm
V}}(J^{\alpha},\theta^{\alpha})$ $\displaystyle+\varsigma^{\alpha}_{\rm
S}(\\{{\bf D}^{\alpha}\\})[\Psi^{\alpha}_{{\rm S}}({\bf
C}^{\alpha})+\Psi^{\alpha}_{\Gamma}({\bf
C}^{\alpha},\\{{\bm{\Gamma}}^{\alpha}\\})]$
$\displaystyle+\varsigma^{\alpha}_{\rm F}(\\{{\bf
D}^{\alpha}\\})\circ[\Psi^{\alpha}_{\rm F}({\bf
C}^{\alpha})+\Psi^{\alpha}_{\Phi}({\bf
C}^{\alpha},\\{{\bm{\Gamma}}^{\alpha}\\})]$ $\displaystyle+\Psi^{\alpha}_{\rm
A}({\bf
C}^{\alpha},\\{{\bm{\Delta}}^{\alpha}\\})+\Psi^{\alpha}_{\rm\theta}(\theta^{\alpha})$
$\displaystyle+\Psi^{\alpha}_{\sigma}(J^{\alpha})+\Psi^{\alpha}_{\rm
D}(\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla_{0}^{\alpha}{\bf D}^{\alpha}\\}),$
(123) $\displaystyle U^{\alpha}({\bf C}^{\alpha},$
$\displaystyle\,\eta^{\alpha},\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla_{0}^{\alpha}{\bf
D}^{\alpha}\\})=$ $\displaystyle\quad\varsigma^{\alpha}_{\rm V}(\\{{\bf
D}^{\alpha}\\},{\bf C}^{\alpha})U^{\alpha}_{{\rm
V}}(J^{\alpha},\eta^{\alpha})$ $\displaystyle+\varsigma^{\alpha}_{\rm
S}(\\{{\bf D}^{\alpha}\\})[U^{\alpha}_{{\rm S}}({\bf
C}^{\alpha})+U^{\alpha}_{\Gamma}({\bf
C}^{\alpha},\\{{\bm{\Gamma}}^{\alpha}\\})]$
$\displaystyle+\varsigma^{\alpha}_{\rm F}(\\{{\bf
D}^{\alpha}\\})\circ[U^{\alpha}_{\rm F}({\bf
C}^{\alpha})+U^{\alpha}_{\Phi}({\bf
C}^{\alpha},\\{{\bm{\Gamma}}^{\alpha}\\})]$ $\displaystyle+U^{\alpha}_{\rm
A}({\bf
C}^{\alpha},\\{{\bm{\Delta}}^{\alpha}\\})+U^{\alpha}_{\rm\theta}(\eta^{\alpha})$
$\displaystyle+U^{\alpha}_{\sigma}(J^{\alpha})+U^{\alpha}_{\rm
D}(\\{{\bm{\xi}}^{\alpha}\\},\\{\nabla_{0}^{\alpha}{\bf D}^{\alpha}\\}).$
(124)

In $\Psi^{\alpha}$, volumetric equilibrium free energy for the entire
constituent $\alpha$ is $\Psi^{\alpha}_{\rm V}$, including isotropic
thermoelastic coupling. Deviatoric equilibrium energy of the isotropic matrix
is $\Psi^{\alpha}_{\rm S}$. Viscoelastic configurational energy of the
isotropic matrix is $\Psi^{\alpha}_{\Gamma}$. Anisotropic deviatoric
equilibrium free energy from fibrous microstructures is $\Psi^{\alpha}_{\rm
F}$. Configurational energy, often but not always anisotropic, from fibers is
$\Psi^{\alpha}_{\Phi}$. Energy from fiber activation is $\Psi^{\alpha}_{\rm
A}$. Thermal energy of specific heat is $\Psi^{\alpha}_{\theta}$. Energy from
a nonzero reference pressure is $\Psi^{\alpha}_{\sigma}$. Surface energy from
fractures, tears, and other damage is contained in $\Psi^{\alpha}_{\rm D}$.
Fully analogous descriptors apply to internal energy contributions in
$U^{\alpha}$. Forms (III.2) and (III.2) are not the most mathematically and
physically general, but they are sufficient for soft tissue materials of
present interest given the scope of available data on their properties and
response. A few, or even most, terms vanish for certain classes of materials
(e.g., isotropic solids, viscoelastic fluids, gas phases, and so forth). All
functions in (III.2) and (III.2) in materials with heterogeneous properties
can further depend explicitly on ${\bf X}^{\alpha}$, omitted in the arguments
for brevity. Dependence of state-dependent metric tensors is implicit, for
example in $J^{\alpha}$ and scalar functions of certain vectors and tensors.

Ideal gas EOS. For gaseous fluids such as air in the lung, an ideal gas model
[74] is sufficient. For the ideal gas, $\psi^{\alpha}=(\Psi^{\alpha}_{\rm
V}+\Psi^{\alpha}_{\theta})/\rho^{\alpha}_{{\rm R}0}$ and
$u^{\alpha}=(U^{\alpha}_{\rm V}+U^{\alpha}_{\theta})/\rho^{\alpha}_{{\rm
R}0}$. At a reference state, the following conditions hold: $J^{\alpha}=1$,
$\theta^{\alpha}=\theta^{\alpha}_{0}$, $\eta^{\alpha}=\eta^{\alpha}_{0}$,
$\rho^{\alpha}_{0}=n^{\alpha}_{0}\rho^{\alpha}_{{\rm R}0}$, $p^{\alpha}_{\rm
V}=p^{\alpha}_{0}=n^{\alpha}_{0}p^{\alpha}_{{\rm R}0}$, and
$c^{\alpha}_{\epsilon}=c^{\alpha}_{\epsilon 0}$. Quantities with zero
subscripts are constants; $p^{\alpha}_{\rm V}=-\frac{1}{3}{\rm
tr}\bar{\bm{\sigma}}^{\alpha}$ is the partial inviscid pressure. From the
identity $\partial J^{\alpha}/\partial{\bf
C}^{\alpha}=\frac{1}{2}J^{\alpha}({\bf C}^{\alpha})^{-1}$, the stress
contribution is spherical: $\bar{\bm{\sigma}}^{\alpha}=-p^{\alpha}_{\rm V}{\bf
1}$. The ideal gas constant is ${\mathfrak{R}}^{\alpha}$. Equation of state
(EOS) and internal energy function are

$p^{\alpha}_{\rm V}=\rho^{\alpha}{\mathfrak{R}}^{\alpha}\theta^{\alpha},\qquad
u^{\alpha}=c^{\alpha}_{\epsilon 0}\theta^{\alpha}.$ (125)

From (125) and
$\psi^{\alpha}=u^{\alpha}+\theta^{\alpha}(\partial\psi^{\alpha}/\partial\theta^{\alpha})$,
it follows that

$\displaystyle\psi^{\alpha}(J^{\alpha},\theta^{\alpha})=-{\mathfrak{R}}^{\alpha}\theta^{\alpha}\ln
J^{\alpha}-c^{\alpha}_{\epsilon
0}\theta^{\alpha}[\ln(\theta^{\alpha}/\theta^{\alpha}_{0})-1],$ (126)
$\displaystyle u^{\alpha}(J^{\alpha},\eta^{\alpha})=c^{\alpha}_{\epsilon
0}\theta^{\alpha}_{0}(J^{\alpha})^{-\gamma^{\alpha}_{0}}\exp(\eta^{\alpha}/c^{\alpha}_{\epsilon
0}).$ (127)

noting $\eta^{\alpha}_{0}=0$ and
$\gamma^{\alpha}_{0}={\mathfrak{R}}^{\alpha}/c^{\alpha}_{\epsilon 0}$. Thermal
stress tensor is
${\bm{\beta}}^{\alpha}=\rho^{\alpha}{\mathfrak{R}}^{\alpha}({\bf
C}^{\alpha})^{-1}$ in (107), and $c^{\alpha}_{\epsilon}=c^{\alpha}_{\epsilon
0}$ in (106).

Condensed matter EOS. For solid and liquid tissue phases, an EOS combining the
third-order logarithmic form used for high-pressure physics [93, 94] with an
exponential form for tissue mechanics [95] is sufficiently general for the
present applications. Thermoelastic coupling is linear and isotropic with
constant volumetric expansion coefficient $A^{\alpha}$, and specific heat
$c^{\alpha}_{\epsilon}$ is constant. Reference temperature is
$\theta^{\alpha}_{0}$, and reference pressure is $p^{\alpha}_{{\rm R}0}$. The
reference isothermal bulk modulus is $B^{\alpha}_{\theta}$, and the pressure
derivative of the isothermal bulk modulus in the reference state is
$B^{\alpha}_{\theta{\rm p}}$. Analogously, the isentropic bulk modulus and
pressure derivative are $B^{\alpha}_{\eta}$ and $B^{\alpha}_{\eta{\rm p}}$.
Denote a constant controlling exponential stiffening by $k^{\alpha}_{\rm V}$.
Free energies per unit initial volume of constituent $\alpha$ are

$\displaystyle\Psi^{\alpha}_{\rm
V}=\frac{B^{\alpha}_{\theta}}{2}\left[\frac{\exp\\{k^{\alpha}_{\rm V}(\ln
J^{\alpha})^{2}\\}-1}{k^{\alpha}_{\rm V}}-\frac{(B^{\alpha}_{\theta{\rm
p}}-2)(\ln J^{\alpha})^{3}}{3}\right]$
$\displaystyle\quad-A^{\alpha}B^{\alpha}_{\theta}(\theta^{\alpha}-\theta^{\alpha}_{0})\ln
J^{\alpha},\quad\Psi^{\alpha}_{\sigma}=-p^{\alpha}_{{\rm R}0}\ln J^{\alpha},$
(128) $\displaystyle\Psi^{\alpha}_{\theta}=-\rho^{\alpha}_{{\rm
R}0}c^{\alpha}_{\epsilon}[\theta^{\alpha}\ln(\theta^{\alpha}/\theta^{\alpha}_{0})-(\theta^{\alpha}-\theta^{\alpha}_{0})].$
(129)

The contribution to stress $\bar{\bm{\sigma}}^{\alpha}$ from
$\Psi^{\alpha}_{\rm V}$ and $\Psi^{\alpha}_{\sigma}$ is spherical, with Cauchy
pressure

$\displaystyle p^{\alpha}_{\rm V}=-\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\frac{\partial(\Psi^{\alpha}_{\rm V}+\Psi^{\alpha}_{\sigma})}{\partial\ln
J^{\alpha}}$ $\displaystyle=-\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}B^{\alpha}_{\theta}\ln J^{\alpha}[\exp\\{k^{\alpha}_{\rm V}(\ln
J^{\alpha})^{2}\\}-{\textstyle{\frac{1}{2}}}(B^{\alpha}_{\theta{\rm p}}-2)\ln
J^{\alpha}]$ $\displaystyle\qquad+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}A^{\alpha}B^{\alpha}_{\theta}(\theta^{\alpha}-\theta^{\alpha}_{0})+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}p^{\alpha}_{{\rm R}0}.$ (130)

From (15), if $c^{\alpha}=\rho^{\alpha}\partial_{t}\ln\sqrt{g}$,
$\partial_{t}\ln\sqrt{g}=D^{\alpha}_{t}\ln\sqrt{G^{\alpha}}$, and
$\rho_{0}^{\alpha}=\rho_{0}^{\alpha}({\bf X}^{\alpha})$, then
$\rho^{\alpha}_{0}=\rho^{\alpha}J^{\alpha}\Rightarrow{\rho^{\alpha}}/{\rho^{\alpha}_{{\rm
R}0}}=n^{\alpha}_{0}/J^{\alpha}$. The thermal stress tensor and Grüneisen
tensor are

${\bm{\beta}}^{\alpha}=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}A^{\alpha}B^{\alpha}_{\theta}({\bf
C}^{\alpha})^{-1},\quad{\bm{\gamma}}^{\alpha}=\frac{A^{\alpha}B^{\alpha}_{\theta}}{\rho^{\alpha}_{{\rm
R}0}c^{\alpha}_{\epsilon}}({\bf C}^{\alpha})^{-1}.$ (131)

The scalar Grüneisen constant is
$\gamma^{\alpha}_{0}=A^{\alpha}B^{\alpha}_{\theta}/(\rho^{\alpha}_{{\rm
R}0}c^{\alpha}_{\epsilon})$. Internal energy complementary to (III.2) and
(129) is

$\displaystyle U^{\alpha}_{\rm
V}=\frac{B^{\alpha}_{\eta}}{2}\left[\frac{\exp\\{k^{\alpha}_{\rm V}(\ln
J^{\alpha})^{2}\\}-1}{k^{\alpha}_{\rm V}}-\frac{(B^{\alpha}_{\eta{\rm
p}}-2)(\ln J^{\alpha})^{3}}{3}\right]$
$\displaystyle\qquad-\rho^{\alpha}_{{\rm
R}0}\theta_{0}^{\alpha}\gamma_{0}^{\alpha}\eta^{\alpha}\ln J^{\alpha},\quad
U^{\alpha}_{\sigma}=-p^{\alpha}_{{\rm R}0}\ln J^{\alpha},$ (132)
$\displaystyle U^{\alpha}_{\theta}=\rho^{\alpha}_{{\rm
R}0}\theta^{\alpha}_{0}\eta^{\alpha}[1+\eta^{\alpha}/(2c^{\alpha}_{\epsilon})].$
(133)

Pressure and temperature from (III.2) and (133) are

$\displaystyle p^{\alpha}_{\rm V}=-\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\frac{\partial(U^{\alpha}_{\rm V}+U^{\alpha}_{\sigma})}{\partial\ln
J^{\alpha}}$ $\displaystyle=-\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}B^{\alpha}_{\eta}\ln J^{\alpha}[\exp\\{k^{\alpha}_{\rm V}(\ln
J^{\alpha})^{2}\\}-{\textstyle{\frac{1}{2}}}(B^{\alpha}_{\eta{\rm p}}-2)\ln
J^{\alpha}]$
$\displaystyle\qquad\qquad+{\rho^{\alpha}}\theta^{\alpha}_{0}\gamma^{\alpha}_{0}\eta^{\alpha}+\rho^{\alpha}p^{\alpha}_{{\rm
R}0}/\rho^{\alpha}_{{\rm R}0},$ (134)
$\displaystyle\theta^{\alpha}=\frac{1}{\rho^{\alpha}_{{\rm
R}0}}\frac{\partial(U^{\alpha}_{\rm
V}+U^{\alpha}_{\theta})}{\partial\eta^{\alpha}}=\theta_{0}^{\alpha}[1+\frac{\eta^{\alpha}}{c^{\alpha}_{\epsilon}}-\gamma_{0}^{\alpha}\ln
J^{\alpha}].$ (135)

Bulk moduli $B^{\alpha}_{\theta}$ and $B^{\alpha}_{\eta}$ are nonnegative, and
$k^{\alpha}_{\rm V}$ should be nonnegative for stiffening under large strain
typical of soft tissues, compressive or tensile. If $B^{\alpha}_{\theta{\rm
p}}>2$, the material stiffens in compression and softens in tension, and vice-
versa for $B^{\alpha}_{\theta{\rm p}}<2$. Similar statements hold for
$B^{\alpha}_{\eta{\rm p}}$. Energy functions (III.2) and (III.2) are not
(poly)convex in $J^{\alpha}$. Polyconvexity is appealing for existence of
unique solutions to boundary value problems [96] but is not essential. If
$\varsigma^{\alpha}_{\rm V}<1$, $p^{\alpha}_{\rm V}$ contributions from
$\Psi^{\alpha}_{\rm V}$ and $U^{\alpha}_{\rm V}$ (i.e., all terms except
rightmost in (III.2) and (III.2) with $p^{\alpha}_{{\rm R}0}$) require
multiplication by $\varsigma^{\alpha}_{\rm V}$, as do ${\bm{\beta}}^{\alpha}$,
${\bm{\gamma}}^{\alpha}$, and $\gamma_{0}^{\alpha}$.

Deviatoric matrix equilibrium. Deviatoric deformation gradient and deformation
tensor are, with $f=f(\tilde{\bf C}^{\alpha})$ a generic differentiable
function of its argument,

$\displaystyle\tilde{\bf F}^{\alpha}=(J^{\alpha})^{-1/3}{\bf
F}^{\alpha},\qquad\tilde{\bf C}^{\alpha}=(J^{\alpha})^{-2/3}{\bf C}^{\alpha},$
(136) $\displaystyle\frac{\partial f}{\partial{\bf
C}^{\alpha}}=(J^{\alpha})^{-2/3}\left[\frac{\partial f}{\partial\tilde{\bf
C}^{\alpha}}-\frac{1}{3}\left(\frac{\partial f}{\partial\tilde{\bf
C}^{\alpha}}:{\bf C}^{\alpha}\right)({\bf C}^{\alpha})^{-1}\right].$ (137)

Let $\mu^{\alpha}_{\rm S}\geq 0$ be a shear modulus. Energy is [96]

$\Psi^{\alpha}_{\rm S}=U^{\alpha}_{\rm
S}={\textstyle{\frac{1}{2}}}\mu^{\alpha}_{\rm S}({\rm tr}\,\tilde{\bf
C}^{\alpha}-3).$ (138)

From (137), the contribution of (138) to Cauchy stress is

$\displaystyle{\bm{\sigma}}^{\alpha}_{\rm
S}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
F}^{\alpha}\frac{\partial\Psi^{\alpha}_{\rm S}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
F}^{\alpha}\frac{\partial U^{\alpha}_{\rm S}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}$
$\displaystyle\quad=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\mu^{\alpha}_{\rm S}[\tilde{\bf
B}^{\alpha}-{\textstyle{\frac{1}{3}}}({\rm tr}\tilde{\bf B}^{\alpha}){\bf
1}],\quad\tilde{\bf B}^{\alpha}=\tilde{\bf F}^{\alpha}(\tilde{\bf
F}^{\alpha})^{\mathsf{T}}.$ (139)

This contribution is linear in spatial deformation tensor $\tilde{\bf
B}^{\alpha}$, traceless, and ultimately scaled by $\varsigma^{\alpha}_{\rm
S}$. Further nonlinearity can be furnished by $\Psi^{\alpha}_{\rm F}$ and
$U^{\alpha}_{\rm F}$. Function (138) is polyconvex [96] and isotropic.

Fiber equilibrium. Let index $k$ denote a fiber family of reference alignment
by unit vector ${\bm{\iota}}^{\alpha}_{k}$. Let
$\kappa^{\alpha}_{k}\in[0,\frac{1}{3}]$ be dispersion constants. Structure
tensors [20] are

${\bf H}_{k}^{\alpha}=\kappa^{\alpha}_{k}{\bf
1}+(1-3\kappa^{\alpha}_{k}){\bm{\iota}}^{\alpha}_{k}\otimes{\bm{\iota}}^{\alpha}_{k}.$
(140)

Strain energy contributions are of functional forms

$\displaystyle\Psi^{\alpha}_{\rm F}=\Psi^{\alpha}_{\rm F}(\tilde{\bf
C}^{\alpha},{\bf H}_{k}^{\alpha}({\bf X}^{\alpha}))=U^{\alpha}_{\rm
F}(\tilde{\bf C}^{\alpha},{\bf H}_{k}^{\alpha}({\bf X}^{\alpha})),$ (141)

with ${\bf H}_{k}^{\alpha}$ time-independent at ${\bf X}^{\alpha}$ (i.e., not
transient state variables). Depending on the number of fiber families $k$ and
their orientations, different scalar invariants entering (141) are possible.
For the current presentation, one invariant per fiber family is sufficient:
$I^{\alpha}_{k}=\tilde{\bf C}^{\alpha}:{\bf H}^{\alpha}_{k}$. The particular
form of (141) is polyconvex [16, 96, 9, 59]:

$\displaystyle\Psi^{\alpha}_{\rm F}=\sum_{k}\Psi^{\alpha}_{{\rm F}k}$
$\displaystyle\quad=\sum_{k}\frac{\mu^{\alpha}_{k}}{4k^{\alpha}_{k}}\\{\exp[k_{k}^{\alpha}(I^{\alpha}_{k}-1)^{2}]-1\\}{\rm
H}(I^{\alpha}_{k}-1).$ (142)

A fiber modulus and stiffening coefficient are $\mu^{\alpha}_{k}\geq 0$ and
$k^{\alpha}_{k}>0$. Optional right-continuous Heaviside function is ${\rm
H}(\cdot)$; this disables fiber stiffness for buckling in compression along
${\bm{\iota}}^{\alpha}_{k}$. Contributions to stress are traceless:

$\displaystyle{\bm{\sigma}}^{\alpha}_{\rm
F}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
F}^{\alpha}\frac{\partial\Psi^{\alpha}_{\rm F}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
F}^{\alpha}\frac{\partial U^{\alpha}_{\rm F}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}$
$\displaystyle\quad=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}\mu^{\alpha}_{k}(I^{\alpha}_{k}-1)\exp[k_{k}^{\alpha}(I^{\alpha}_{k}-1)^{2}]{\rm
H}(I^{\alpha}_{k}-1)\tilde{\bf h}^{\alpha}_{k},$ $\displaystyle\tilde{\bf
h}^{\alpha}_{k}=\tilde{\bf F}^{\alpha}{\bf H}^{\alpha}_{k}(\tilde{\bf
F}^{\alpha})^{\mathsf{T}}-{\textstyle{\frac{1}{3}}}{\rm tr}[\tilde{\bf
F}^{\alpha}{\bf H}^{\alpha}_{k}(\tilde{\bf F}^{\alpha})^{\mathsf{T}}]{\bf 1}.$
(143)

Family $k$ is isotropic as $\kappa_{k}^{\alpha}\rightarrow\frac{1}{3}$. Fiber
compressibility is encompassed by the condensed matter EOS for phase $\alpha$
rather than distinct energetic terms. Stress contributions from (III.2) are
affected by $\varsigma^{\alpha}_{\rm F}$ if fibers are damaged.

Matrix viscoelasticity. The viscoelastic formulation combines features from
prior works [97, 98, 99, 95, 100] in a thermodynamically consistent manner.
Let $\\{{\bm{\Gamma}}^{\alpha}\\}\rightarrow\\{{\bm{\Gamma}}^{\alpha}_{{\rm
V}l},{\bm{\Gamma}}^{\alpha}_{{\rm S}m},{\bm{\Gamma}}^{\alpha}_{{\Phi}k,n}\\}$
be internal strain-like configurational variables for constituent $\alpha$.
Index $l$ spans a set of discrete relaxation time constants
$\tau^{\alpha}_{{\rm V}l}=\tau^{\alpha}_{{\rm V}1},\ldots$ for viscoelastic
relaxation processes associated with volumetric deformation of the matrix.
Index $m$ spans times $\tau^{\alpha}_{{\rm S}m}$ associated with deviatoric
(shear) deformation of the matrix. Index $n$ spans times
$\tau^{\alpha}_{{\Phi}k,n}$ associated with fiber family $k$ discussed in the
next subsection. Internal stresses $\\{{\bf Q}^{\alpha}_{{\rm V}l},{\bf
Q}^{\alpha}_{{\rm S}m}\\}$ conjugate to the matrix internal strains, in
coordinates referred to ${\mathfrak{M}}^{\alpha}$, obey [98, 99]

$\displaystyle{\bf Q}^{\alpha}_{{\rm
V}l}=-\frac{\partial\Psi^{\alpha}_{\Gamma}}{\partial{\bm{\Gamma}}^{\alpha}_{{\rm
V}l}}=2\frac{\partial\Psi^{\alpha}_{{\rm V}l}}{\partial{\bf
C}^{\alpha}},\quad{\bf Q}^{\alpha}_{{\rm
S}m}=-\frac{\partial\Psi^{\alpha}_{\Gamma}}{\partial{\bm{\Gamma}}^{\alpha}_{{\rm
S}m}}=2\frac{\partial\Psi^{\alpha}_{{\rm S}m}}{\partial{\bf C}^{\alpha}},$
(144)
$\displaystyle\Psi^{\alpha}_{\Gamma}=U^{\alpha}_{\Gamma}=\sum_{l}\Psi^{\alpha}_{{\rm
V}l}({\bm{\Gamma}}^{\alpha}_{{\rm V}l},{\bf
C}^{\alpha})+\sum_{m}\Psi^{\alpha}_{{\rm S}m}({\bm{\Gamma}}^{\alpha}_{{\rm
S}m},{\bf C}^{\alpha})$
$\displaystyle\quad=\sum_{l}\int{\textstyle{\frac{1}{2}}}{\bf
Q}^{\alpha}_{{\rm V}l}:{\rm d}{\bf
C}^{\alpha}+\sum_{m}\int{\textstyle{\frac{1}{2}}}{\bf Q}^{\alpha}_{{\rm
S}m}:{\rm d}{\bf C}^{\alpha}.$ (145)

Indefinite integrals (III.2) are not needed explicitly. Evolution equations
for internal stresses are

$\displaystyle D^{\alpha}_{t}{\bf Q}^{\alpha}_{{\rm V}l}+{\bf
Q}^{\alpha}_{{\rm V}l}/\tau^{\alpha}_{{\rm
V}l}=2D^{\alpha}_{t}(\partial\hat{\Psi}^{\alpha}_{{\rm V}l}/\partial{\bf
C}^{\alpha}),$ (146) $\displaystyle D^{\alpha}_{t}{\bf Q}^{\alpha}_{{\rm
S}m}+{\bf Q}^{\alpha}_{{\rm S}m}/\tau^{\alpha}_{{\rm
S}m}=2D^{\alpha}_{t}(\partial\hat{\Psi}^{\alpha}_{{\rm S}m}/\partial{\bf
C}^{\alpha}),$ (147) $\displaystyle\hat{\Psi}^{\alpha}_{{\rm
V}l}={\textstyle{\frac{1}{2}}}\beta^{\alpha}_{{\rm V}l}B^{\alpha}_{\theta}(\ln
J^{\alpha})^{2},$ (148) $\displaystyle\hat{\Psi}^{\alpha}_{{\rm
S}m}={\textstyle{\frac{1}{2}}}\beta^{\alpha}_{{\rm S}m}\mu^{\alpha}_{\rm
S}({\rm tr}\,\tilde{\bf C}^{\alpha}-3).$ (149)

Dimensionless factors are $\beta^{\alpha}_{{\rm V}l}\geq
0,\beta^{\alpha}_{{\rm S}m}\geq 0$. Initial conditions and convolution
solutions to (146) and (147) are

$\displaystyle{\bf Q}^{\alpha}_{{\rm V}l0}=2\partial\hat{\Psi}^{\alpha}_{{\rm
V}l}/\partial{\bf C}^{\alpha},\quad{\bf Q}^{\alpha}_{{\rm
S}m0}=2\partial\hat{\Psi}^{\alpha}_{{\rm S}m}/\partial{\bf C}^{\alpha},$ (150)
$\displaystyle{\bf Q}^{\alpha}_{{\rm V}l}(t)={\bf Q}^{\alpha}_{{\rm
V}l0}\exp[{-t}/{\tau^{\alpha}_{{\rm V}l}}]$
$\displaystyle\qquad+\int_{0+}^{t}\exp[{-(t-s)}/{\tau^{\alpha}_{{\rm
V}l}}]D^{\alpha}_{s}(2\partial\hat{\Psi}^{\alpha}_{{\rm V}l}/\partial{\bf
C}^{\alpha}){\rm d}s,$ (151) $\displaystyle{\bf Q}^{\alpha}_{{\rm S}m}(t)={\bf
Q}^{\alpha}_{{\rm S}m0}\exp[{-t}/{\tau^{\alpha}_{{\rm S}m}}]$
$\displaystyle\qquad+\int_{0+}^{t}\exp[{-(t-s)}/{\tau^{\alpha}_{{\rm
S}m}}]D^{\alpha}_{s}(2\partial\hat{\Psi}^{\alpha}_{{\rm S}m}/\partial{\bf
C}^{\alpha}){\rm d}s.$ (152)

Cauchy stress contributions are sums over $l,m$:

$\displaystyle{\bm{\sigma}}^{\alpha}_{\Gamma}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}{\bf F}^{\alpha}\frac{\partial\Psi^{\alpha}_{\Gamma}}{\partial{\bf
C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
F}^{\alpha}\frac{\partial U^{\alpha}_{\Gamma}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}$
$\displaystyle=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}\sum_{l}{\bf
F}^{\alpha}{\bf Q}^{\alpha}_{{\rm V}l}({\bf
F}^{\alpha})^{\mathsf{T}}+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{m}{\bf F}^{\alpha}{\bf Q}^{\alpha}_{{\rm S}m}({\bf
F}^{\alpha})^{\mathsf{T}}.$ (153)

As $t/\tau^{\alpha}_{{\rm V}l}\rightarrow 0$ and $t/\tau^{\alpha}_{{\rm
S}m}\rightarrow 0$, ${\bm{\sigma}}^{\alpha}_{\Gamma}$ in (III.2) sums to the
instantaneous (glassy) viscoelastic stresses

$\displaystyle 2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}\sum_{l}{\bf
F}^{\alpha}\frac{\partial\hat{\Psi}^{\alpha}_{{\rm V}l}}{\partial{\bf
C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}+2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{m}{\bf F}^{\alpha}\frac{\partial\hat{\Psi}^{\alpha}_{{\rm
S}m}}{\partial{\bf C}^{\alpha}}({\bf F}^{\alpha})^{\mathsf{T}}$
$\displaystyle=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{l}\beta^{\alpha}_{{\rm V}l}B^{\alpha}_{\theta}(\ln J^{\alpha}){\bf
1}$ $\displaystyle\qquad+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{m}\beta^{\alpha}_{{\rm S}m}\mu^{\alpha}_{{\rm S}}[\tilde{\bf
B}^{\alpha}-\frac{1}{3}({\rm tr}\tilde{\bf B}^{\alpha}){\bf 1}].$ (154)

As $t/\tau^{\alpha}_{\rm{V}l}\rightarrow\infty$ and
$t/\tau^{\alpha}_{\rm{S}m}\rightarrow\infty$, ${\bf Q}^{\alpha}_{{\rm
V}l}\rightarrow{\bf 0}$ and ${\bf Q}^{\alpha}_{{\rm S}m}\rightarrow{\bf 0}$ so
that ${\bm{\sigma}}^{\alpha}_{\Gamma}\rightarrow{\bf 0}$ in (III.2) for
relaxed equilibrium response.

Fiber viscoelasticity. Dissipative response of fiber families $k=1,\ldots$ is
dictated by internal variables ${\bm{\Gamma}}^{\alpha}_{{\Phi}k,n}$ each with
$n=1,\ldots$ relaxation times $\tau^{\alpha}_{{\Phi}k,n}$ and conjugate
internal stresses ${\bf Q}^{\alpha}_{{\Phi}k,n}$. Internal stresses and
energies are

$\displaystyle{\bf
Q}^{\alpha}_{{\Phi}k,n}=-{\partial\Psi^{\alpha}_{\Phi}}/{\partial{\bm{\Gamma}}^{\alpha}_{{\Phi}k,n}}=2{\partial\Psi^{\alpha}_{{\Phi}k,n}}/{\partial{\bf
C}^{\alpha}},$ (155)
$\displaystyle\Psi^{\alpha}_{\Phi}=U^{\alpha}_{\Phi}=\sum_{k}\Psi^{\alpha}_{{\Phi}k}=\sum_{k}\sum_{n}\Psi^{\alpha}_{{\Phi}k,n}({\bm{\Gamma}}^{\alpha}_{{\Phi}k,n},{\bf
C}^{\alpha})$
$\displaystyle\qquad\qquad=\sum_{k}\sum_{n}\int{\textstyle{\frac{1}{2}}}{\bf
Q}^{\alpha}_{{\Phi}k,n}:{\rm d}{\bf C}^{\alpha}.$ (156)

Evolution equations and stored viscoelastic energies are

$\displaystyle D^{\alpha}_{t}{\bf Q}^{\alpha}_{{\Phi}k,n}+{\bf
Q}^{\alpha}_{{\Phi}k,n}/\tau^{\alpha}_{{\Phi}k,n}=2D^{\alpha}_{t}(\partial\hat{\Psi}^{\alpha}_{{\Phi}k,n}/\partial{\bf
C}^{\alpha}),$ (157)
$\displaystyle\hat{\Psi}^{\alpha}_{{\Phi}k,n}=\frac{\beta^{\alpha}_{{\Phi}k,n}\mu^{\alpha}_{k}}{4k^{\alpha}_{k}}\\{\exp[k_{k}^{\alpha}(I^{\alpha}_{k}-1)^{2}]-1\\}{\rm
H}(I^{\alpha}_{k}-1),$ (158)

with $\beta^{\alpha}_{{\Phi}k,n}\geq 0$. Initial conditions and solutions
(convolution integrals), followed by Cauchy stress terms, are

$\displaystyle{\bf
Q}^{\alpha}_{{\Phi}k,n0}=2\partial\hat{\Psi}^{\alpha}_{{\Phi}k,n}/\partial{\bf
C}^{\alpha},$ (159) $\displaystyle{\bf Q}^{\alpha}_{{\Phi}k,n}(t)={\bf
Q}^{\alpha}_{{\Phi}k,n0}\exp[{-t}/{\tau^{\alpha}_{{\Phi}k,n}}]$
$\displaystyle\qquad+\int_{0+}^{t}\exp[{-(t-s)}/{\tau^{\alpha}_{{\Phi}k,n}}]D^{\alpha}_{s}(2\partial\hat{\Psi}^{\alpha}_{{\Phi}k,n}/\partial{\bf
C}^{\alpha}){\rm d}s,$ (160)
$\displaystyle{\bm{\sigma}}^{\alpha}_{\Phi}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}{\bf F}^{\alpha}\frac{\partial\Psi^{\alpha}_{\Phi}}{\partial{\bf
C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
F}^{\alpha}\frac{\partial U^{\alpha}_{\Phi}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}$
$\displaystyle=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}\sum_{n}{\bf F}^{\alpha}{\bf Q}^{\alpha}_{{\Phi}k,n}({\bf
F}^{\alpha})^{\mathsf{T}}.$ (161)

As $t/\tau^{\alpha}_{{\Phi}k,n}\rightarrow 0$, ${\bm{\sigma}}^{\alpha}_{\Phi}$
in (III.2) becomes the glassy stress

$\displaystyle 2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}\sum_{n}{\bf
F}^{\alpha}\frac{\partial\hat{\Psi}^{\alpha}_{{\Phi}k,n}}{\partial{\bf
C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}\sum_{n}\\{\beta^{\alpha}_{{\Phi}k,n}$
$\displaystyle\quad\times\mu^{\alpha}_{k}(I^{\alpha}_{k}-1)\exp[k_{k}^{\alpha}(I^{\alpha}_{k}-1)^{2}]{\rm
H}(I^{\alpha}_{k}-1)\tilde{\bf h}^{\alpha}_{k}\\}.$ (162)

As $t/\tau^{\alpha}_{{\Phi}k,n}\rightarrow\infty$, ${\bf
Q}^{\alpha}_{{\Phi}k,n}\rightarrow{\bf 0}$ leading to
${\bm{\sigma}}^{\alpha}_{\Phi}\rightarrow{\bf 0}$ in (III.2) .

Active tension. Electrochemistry of soft tissue cellular activation [101] is
beyond the present scope. A phenomenological approach is instead, generalizing
other continuum models [102, 23, 103]. Let
$\\{{\bm{\Delta}}^{\alpha}\\}\rightarrow\\{\Delta^{\alpha}_{k}\\}({\bf
X}^{\alpha},t)$ be a set of scalar internal variables associated with
potentially active fiber families $k$ in phase $\alpha$. These variables can
include internal strains in contractile elements and time-dependent switching
functions. Define the fiber orientation tensors ${\bf H}^{\alpha}_{k}$ as in
(140). In many models, cells are fully aligned such that $k=1$ and
$\kappa^{\alpha}_{1}=0$ [102, 23], but this is inessential [103]. Stretch in
the fiber direction is $\lambda^{\alpha}_{k}=\sqrt{I^{\alpha}_{k}}$. A generic
energy function is

$\Psi^{\alpha}_{\rm A}=U^{\alpha}_{\rm A}=\sum_{k}\Psi^{\alpha}_{{\rm
A}k}=\sum_{k}[\Lambda^{\alpha}_{k}(\lambda^{\alpha}_{k},\\{\Delta^{\alpha}_{k}\\})+\chi^{\alpha}_{k}(\\{\Delta^{\alpha}_{k}\\})].$
(163)

Strain energy functions $\Lambda^{\alpha}_{k}$ furnish active stress terms

$\displaystyle{\bm{\sigma}}^{\alpha}_{\rm
A}=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
F}^{\alpha}\frac{\partial\Psi^{\alpha}_{\rm A}}{\partial{\bf C}^{\alpha}}({\bf
F}^{\alpha})^{\mathsf{T}}=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}\frac{1}{\lambda^{\alpha}_{k}}\frac{\partial\Lambda^{\alpha}_{k}}{\partial\lambda^{\alpha}_{k}}\tilde{\bf
h}^{\alpha}_{k}.$ (164)

Energy functions $\chi^{\alpha}_{k}$ ensure nonnegative net dissipation.
Stresses ${\bm{\sigma}}^{\alpha}_{\rm A}$ should vanish for passive
conditions.

Damage. Internal variables $\\{{\bf
D}^{\alpha}\\}\rightarrow\\{\bar{D}^{\alpha},D^{\alpha}_{k}\\}$. Scalar damage
measures in the isotropic matrix are $\bar{D}^{\alpha}\in[0,1]$, and
$D^{\alpha}_{k}\in[0,1]$ are scalar functions for each fiber family $k$. All
are akin to order parameters in phase-field fracture theory [37, 26].
Degradation functions in (121)–(III.2) are, with
$\bar{\vartheta}^{\alpha}\in[0,\infty),\vartheta^{\alpha}_{k}\in[0,\infty)$
constants,

$\displaystyle\varsigma^{\alpha}_{\rm{V}}=[1-\bar{D}^{\alpha}{\rm H}(\ln
J^{\alpha})]^{\bar{\vartheta}^{\alpha}},\quad\varsigma^{\alpha}_{\rm{S}}=(1-\bar{D}^{\alpha})^{\bar{\vartheta}^{\alpha}},$
(165) $\displaystyle\varsigma^{\alpha}_{\rm
F}\circ(\cdot)=\varsigma^{\alpha}_{\rm
F}\circ\sum_{k}(\cdot)_{k}=\sum_{k}\varsigma^{\alpha}_{{\rm
F}k}(\cdot)_{k}=\sum_{k}(1-D^{\alpha}_{k})^{\vartheta^{\alpha}_{k}}(\cdot)_{k},$
$\displaystyle\varsigma^{\alpha}_{{\rm
F}k}=(1-D^{\alpha}_{k})^{\vartheta^{\alpha}_{k}}.$ (166)

The Heaviside function in $\varsigma^{\alpha}_{\rm V}$ prevents degradation in
compression so the bulk modulus is maintained [37, 89]. Operator
$\varsigma^{\alpha}_{\rm F}$ in (III.2) and (III.2) is applied to a sum of
energetic contributions (hyperelastic, viscoelastic, and active) over all
families $k$ via (III.2). Partial stress less viscous stress of constituent
$\alpha$ is
$\bar{\bm{\sigma}}^{\alpha}={\bm{\sigma}}^{\alpha}-\hat{\bm{\sigma}}^{\alpha}$.
For an ideal gas,
$\bar{D}^{\alpha}=\bar{\vartheta}^{\alpha}=0\rightarrow\varsigma^{\alpha}_{\rm
V}=1$, so
$\bar{\bm{\sigma}}^{\alpha}=-\rho^{\alpha}{\mathfrak{R}}^{\alpha}\theta^{\alpha}{\bf
1}$ by (125). For solid and liquid $\alpha$, applying (III.2), this stress is
the sum of (III.2), (III.2), (III.2), (III.2), (III.2), and (164) scaled by
one or more functions in (165) and (III.2):

$\displaystyle\bar{\bm{\sigma}}^{\alpha}=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}[1-\bar{D}^{\alpha}{\rm H}(\ln
J^{\alpha})]^{\bar{\vartheta}^{\alpha}}B^{\alpha}_{\theta}$
$\displaystyle\qquad\quad\times\\{\ln J^{\alpha}[\exp\\{k^{\alpha}_{\rm V}(\ln
J^{\alpha})^{2}\\}-{\textstyle{\frac{1}{2}}}(B^{\alpha}_{\theta{\rm p}}-2)\ln
J^{\alpha}]$
$\displaystyle\qquad\qquad\quad-A^{\alpha}(\theta^{\alpha}-\theta^{\alpha}_{0})\\}{\bf
1}-({\rho^{\alpha}}/{\rho^{\alpha}_{{\rm R}0})p^{\alpha}_{{\rm R}0}}{\bf 1}$
$\displaystyle\qquad+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}(1-\bar{D}^{\alpha})^{\bar{\vartheta}^{\alpha}}\mu^{\alpha}_{\rm
S}[\tilde{\bf B}^{\alpha}-{\textstyle{\frac{1}{3}}}({\rm tr}\tilde{\bf
B}^{\alpha}){\bf 1}]$
$\displaystyle\qquad+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}(1-{D}^{\alpha}_{k})^{\vartheta^{\alpha}_{k}}\mu^{\alpha}_{k}(I^{\alpha}_{k}-1)$
$\displaystyle\qquad\qquad\qquad\qquad\times\exp[k_{k}^{\alpha}(I^{\alpha}_{k}-1)^{2}]{\rm
H}(I^{\alpha}_{k}-1)\tilde{\bf h}^{\alpha}_{k}$
$\displaystyle\qquad+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}(1-\bar{D}^{\alpha})^{\bar{\vartheta}^{\alpha}}\sum_{l}{\bf
F}^{\alpha}{\bf Q}^{\alpha}_{{\rm V}l}({\bf F}^{\alpha})^{\mathsf{T}}$
$\displaystyle\qquad+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}(1-\bar{D}^{\alpha})^{\bar{\vartheta}^{\alpha}}\sum_{m}{\bf
F}^{\alpha}{\bf Q}^{\alpha}_{{\rm S}m}({\bf F}^{\alpha})^{\mathsf{T}}$
$\displaystyle\qquad+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}[(1-{D}^{\alpha}_{k})^{\vartheta^{\alpha}_{k}}\sum_{n}{\bf
F}^{\alpha}{\bf Q}^{\alpha}_{{\Phi}k,n}({\bf F}^{\alpha})^{\mathsf{T}}]$
$\displaystyle\qquad+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}\frac{1}{\lambda^{\alpha}_{k}}\frac{\partial\Lambda^{\alpha}_{k}}{\partial\lambda^{\alpha}_{k}}\tilde{\bf
h}^{\alpha}_{k}.$ (167)

To preclude damage induced by normal muscle contraction, the final term is
decoupled from $\\{{\bf D}^{\alpha}\\}$ consistently with (III.2) and (III.2).
If (III.2) is used instead of (III.2), spherical terms on the right in (III.2)
should appeal to (III.2) rather than (III.2).

Let $\Psi^{\alpha}_{\rm D}=U^{\alpha}_{\rm D}$ comprise cohesive and surface
energies of fracture per unit referential volume scaled by contributions of
dimensionless Finsler-type metric $\hat{\bf G}^{\alpha}$ in (113) [63, 59]
where $\hat{G}^{\alpha}=\hat{G}^{\alpha}(\\{\bm{\xi}^{\alpha}\\},{\bf
X}^{\alpha},{\bf x})$ via (87). As in phase-field theory, quadratic forms for
matrix and fiber $k$ contributions are

$\displaystyle\hat{\Psi}^{\alpha}_{\rm D}=\Psi^{\alpha}_{\rm
D}/\sqrt{\hat{G}^{\alpha}}=\bar{E}^{\alpha}_{\rm
C}|\bar{D}^{\alpha}|^{2}+\bar{\Upsilon}^{\alpha}\bar{l}_{\rm
R}^{\alpha}|\nabla^{\alpha}_{0}\bar{D}^{\alpha}|^{2}$
$\displaystyle\qquad\qquad\quad+\sum_{k}[E^{\alpha}_{{\rm
C}k}|D^{\alpha}_{k}|^{2}+{\Upsilon}^{\alpha}_{k}{l}^{\alpha}_{{\rm
R}k}|\nabla^{\alpha}_{0}{D}^{\alpha}_{k}|^{2}].$ (168)

In (III.2), cohesive energies per unit volume are $\bar{E}^{\alpha}_{\rm C}$
and $E^{\alpha}_{{\rm C}k}$, surface energies are $\bar{\Upsilon}^{\alpha}$
and $\Upsilon^{\alpha}_{k}$, and gradient regularization lengths are
$\bar{l}^{\alpha}_{\rm R}=\bar{\alpha}^{\alpha}\bar{l}^{\alpha}$ and
${l}^{\alpha}_{{\rm R}k}=\alpha^{\alpha}_{k}l^{\alpha}_{k}$, all non-negative
constants. For solids, typically $\bar{E}^{\alpha}_{\rm
C}=\bar{\Upsilon}^{\alpha}/\bar{l}^{\alpha}$ and $E^{\alpha}_{{\rm
C}k}={\Upsilon}^{\alpha}_{k}/{l}^{\alpha}_{k}$. Dimensionless factors
$\bar{\alpha}^{\alpha}\in[0,\infty)$ and
$\alpha^{\alpha}_{k}l^{\alpha}_{k}\in[0,\infty)$ allow independent cohesive
energies and gradient regularization lengths [89]. For cavitation of a fluid,
gradient terms can be dropped [39] (i.e., $\bar{\alpha}^{\alpha}\rightarrow
0$); cohesion energy $\bar{E}^{\alpha}_{{\rm C}}$ will capture fracture of the
fluid for tensile pressure. Isotropic surface energies are assumed for
gradient terms of (III.2). These could be extended to anisotropic
contributions [26] if data exist. However, $\\{D^{\alpha}_{k}\\}$ furnish
stress-damage anisotropy regardless.

###  III.3 Kinetics

Viscous stress. Isotropic Newtonian behavior is usually adequate for each
constituent $\alpha$, with $\hat{B}^{\alpha}(\theta^{\alpha})\geq 0$ and
$\hat{\mu}^{\alpha}(\theta^{\alpha})\geq 0$ possibly temperature-dependent
bulk and shear viscosities. Viscous stresses and dissipation are

$\displaystyle\hat{\bm{\sigma}}^{\alpha}=[\hat{B}^{\alpha}-{\textstyle{\frac{2}{3}}}\hat{\mu}^{\alpha}]{\rm
tr}({\bf d}^{\alpha}){\bf 1}+2\hat{\mu}^{\alpha}{\bf d}^{\alpha}$ (169)
$\displaystyle\hat{\mathfrak{D}}^{\alpha}=\hat{\bm{\sigma}}^{\alpha}:{\bf
d}^{\alpha}$ $\displaystyle\quad\,\,=\hat{B}^{\alpha}|{\rm tr}({\bf
d}^{\alpha})|^{2}+2\hat{\mu}^{\alpha}|{\bf
d}^{\alpha}-{\textstyle{\frac{1}{3}}}{\rm tr}({\bf d}^{\alpha}){\bf
1}|^{2}\geq 0.$ (170)

Viscous pressure and shear are
$\hat{p}^{\alpha}=-\hat{B}^{\alpha}\nabla\cdot{\bm{\upsilon}}^{\alpha}$ and
$\hat{\bm{\sigma}}^{\alpha}_{\rm S}$, whereby
$\hat{\bm{\sigma}}^{\alpha}=-\hat{p}^{\alpha}{\bf
1}+\hat{\bm{\sigma}}^{\alpha}_{\rm S}$.

Viscoelasticity. Viscoelastic internal state variables are the subset of
$\\{{\bm{\xi}}^{\alpha}\\}$ consisting of $\\{{\bm{\Gamma}}^{\alpha}_{{\rm
V}l},{\bm{\Gamma}}^{\alpha}_{{\rm S}m},{\bm{\Gamma}}^{\alpha}_{{\Phi}k,n}\\}$.
Conjugate forces entering (108) are a subset of $\\{{\bm{\pi}}^{\alpha}\\}$:

$\begin{split}{\bm{\pi}}^{\alpha}_{{\rm V}l}&=-\varsigma_{\rm
S}^{\alpha}\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
Q}^{\alpha}_{{\rm V}l},\qquad{\bm{\pi}}^{\alpha}_{{\rm S}m}=-\varsigma_{\rm
S}^{\alpha}\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
Q}^{\alpha}_{{\rm S}m},\\\ {\bm{\pi}}^{\alpha}_{{\Phi}k,n}&=-\varsigma_{{\rm
F}k}^{\alpha}\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}{\bf
Q}^{\alpha}_{{\Phi}k,n}.\end{split}$ (171)

Kinetic laws for internal variables [97, 98, 95] are

$\begin{split}D^{\alpha}_{t}{\bm{\Gamma}}^{\alpha}_{{\rm V}l}&=\frac{{\bf
Q}^{\alpha}_{{\rm V}l}}{\beta^{\alpha}_{{\rm
V}l}B^{\alpha}_{\theta}\tau^{\alpha}_{{\rm V}l}},\quad
D^{\alpha}_{t}{\bm{\Gamma}}^{\alpha}_{{\rm S}m}=\frac{{\bf Q}^{\alpha}_{{\rm
S}m}}{\beta^{\alpha}_{{\rm S}m}\mu^{\alpha}_{\rm S}\tau^{\alpha}_{{\rm
S}m}},\\\ D^{\alpha}_{t}{\bm{\Gamma}}^{\alpha}_{{\Phi}k,n}&=\frac{{\bf
Q}^{\alpha}_{{\Phi}k,n}}{\beta^{\alpha}_{{\Phi}k,n}\mu^{\alpha}_{k}\tau^{\alpha}_{{\Phi}k,n}}.\end{split}$
(172)

Dissipation from viscoelasticity is nonnegative in (108):

$\begin{split}{\mathfrak{D}}^{\alpha}_{\Gamma}=&\quad\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{l}\frac{\varsigma_{\rm S}^{\alpha}{\bf Q}^{\alpha}_{{\rm V}l}:{\bf
Q}^{\alpha}_{{\rm V}l}}{\beta^{\alpha}_{{\rm
V}l}B^{\alpha}_{\theta}\tau^{\alpha}_{{\rm V}l}}\\\
&+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm R}0}}\sum_{m}\frac{\varsigma_{\rm
S}^{\alpha}{\bf Q}^{\alpha}_{{\rm S}m}:{\bf Q}^{\alpha}_{{\rm
S}m}}{\beta^{\alpha}_{{\rm S}m}\mu^{\alpha}_{\rm S}\tau^{\alpha}_{{\rm
S}m}}\\\ &+\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sum_{k}\sum_{n}\frac{\varsigma_{{\rm F}k}^{\alpha}{\bf
Q}^{\alpha}_{{\Phi}k,n}:{\bf
Q}^{\alpha}_{{\Phi}k,n}}{\beta^{\alpha}_{{\Phi}k,n}\mu^{\alpha}_{k}\tau^{\alpha}_{{\Phi}k,n}}\geq
0.\end{split}$ (173)

Initial conditions for state variables are ${\bm{\Gamma}}^{\alpha}_{{\rm
V}l0}={\bf 0}$, ${\bm{\Gamma}}^{\alpha}_{{\rm S}m}={\bf 0}$, and
${\bm{\Gamma}}^{\alpha}_{{\Phi}k,n}={\bf 0}$. For energies,
$\Psi^{\alpha}_{{\rm V}l}({\bm{0}},{\bf C}^{\alpha})=\hat{\Psi}^{\alpha}_{{\rm
V}l}({\bf C}^{\alpha})$, $\Psi^{\alpha}_{{\rm S}m}({\bm{0}},{\bf
C}^{\alpha})=\hat{\Psi}^{\alpha}_{{\rm S}m}({\bf C}^{\alpha})$,
$\Psi^{\alpha}_{{\Phi}k,n}({\bm{0}},{\bf
C}^{\alpha})=\hat{\Psi}^{\alpha}_{{\Phi}k,n}({\bf C}^{\alpha})$.
Configurational energies are integrated over time as

$\displaystyle\Psi^{\alpha}_{\Gamma}=\sum_{l}\left[\hat{\Psi}^{\alpha}_{{\rm
V}l}-\int_{0}^{t}{\bf Q}^{\alpha}_{{\rm
V}l}:D^{\alpha}_{s}{\bm{\Gamma}}^{\alpha}_{{\rm V}l}\,{\rm d}s\right]$
$\displaystyle\qquad\quad+\sum_{m}\left[\hat{\Psi}^{\alpha}_{{\rm
S}m}-\int_{0}^{t}{\bf Q}^{\alpha}_{{\rm
S}m}:D^{\alpha}_{s}{\bm{\Gamma}}^{\alpha}_{{\rm S}m}\,{\rm d}s\right],$ (174)
$\displaystyle\Psi^{\alpha}_{\Phi}=\sum_{k}\sum_{n}\left[\hat{\Psi}^{\alpha}_{{\Phi}k,n}-\int_{0}^{t}{\bf
Q}^{\alpha}_{{\Phi}k,n}:D^{\alpha}_{s}{\bm{\Gamma}}^{\alpha}_{{\Phi}k,n}{\rm
d}s\right].$ (175)

Active tension. Internal state variables in $\\{{\bm{\xi}}^{\alpha}\\}$ are
the scalar functions $\\{\Delta^{\alpha}_{k}\\}$ with $k$ the fiber family
number. Kinetic equations with initial conditions are imposed directly [102,
23, 103] rather than by more sophisticated electrochemical physics [101]
outside the present scope:

$\begin{split}D^{\alpha}_{t}\\{D^{\alpha}_{t}\Delta^{\alpha}_{k}\\}&=\\{D^{\alpha}_{t}\Delta^{\alpha}_{k}\\}({\bf
X}^{\alpha},t),\\\ \\{\Delta^{\alpha}_{k}\\}({\bf
X}^{\alpha},0)&=\\{\Delta^{\alpha}_{k0}\\}({\bf X}^{\alpha}).\end{split}$
(176)

Evolution equations (176) should implicitly be affected by local states; for
example, the history of fiber damage $\\{D^{\alpha}_{k}\\}$, if severe, should
limit maximum contractile stress. Conjugate forces in (108) are the following
parts of $\\{{\bm{\pi}}^{\alpha}\\}$:

$\begin{split}&\\{\pi^{\alpha}_{{\rm
A}k}\\}(\lambda^{\alpha}_{k},\\{\Delta^{\alpha}_{k}\\})=\\\
&\qquad\varsigma_{{\rm F}k}^{\alpha}\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\frac{\partial}{\partial\\{\Delta^{\alpha}_{k}\\}}[\Lambda^{\alpha}_{k}(\lambda^{\alpha}_{k},\\{\Delta^{\alpha}_{k}\\})+\chi^{\alpha}_{k}(\\{\Delta^{\alpha}_{k}\\})].\end{split}$
(177)

Dissipation from activation or passivation should be nonnegative, to be
ensured by storage-release functions $\chi^{\alpha}_{k}$:

${\mathfrak{D}}^{\alpha}_{\rm A}=-\sum_{k}\\{\pi^{\alpha}_{{\rm
A}k}\\}\cdot\\{D^{\alpha}_{t}\Delta^{\alpha}_{k}\\}\geq 0.$ (178)

Damage. From (III.2) and (III.2), conjugate forces to damage measures for the
matrix in $\\{{\bm{\pi}}^{\alpha}\\}$ and $\\{{\bm{\zeta}}^{\alpha}\\}$ are

$\displaystyle\bar{\pi}^{\alpha}_{\rm
D}=\rho^{\alpha}\frac{\partial\psi^{\alpha}}{\partial\bar{D}^{\alpha}}=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\frac{\partial}{\partial\bar{D}^{\alpha}}[\sqrt{\hat{G}^{\alpha}}\hat{\Psi}^{\alpha}_{\rm
D}+$ $\displaystyle\qquad\qquad\qquad\qquad+\varsigma^{\alpha}_{\rm
V}\Psi^{\alpha}_{\rm V}+\varsigma^{\alpha}_{\rm S}(\Psi^{\alpha}_{\rm
S}+\Psi^{\alpha}_{\Gamma})]$
$\displaystyle\quad=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}[2\sqrt{\hat{G}^{\alpha}}\bar{E}^{\alpha}_{\rm
C}\bar{D}^{\alpha}+\Psi^{\alpha}_{\rm
D}\frac{\partial}{\partial\bar{D}^{\alpha}}\ln\sqrt{\hat{G}^{\alpha}}]$
$\displaystyle\quad\,\,-\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\bar{\vartheta}^{\alpha}[1-\bar{D}^{\alpha}{\rm H}(\ln
J^{\alpha})]^{\bar{\vartheta}^{\alpha}-1}{\rm H}(\ln
J^{\alpha})\Psi^{\alpha}_{\rm V}$
$\displaystyle\quad\,\,-\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\bar{\vartheta}^{\alpha}[1-\bar{D}^{\alpha}]^{\bar{\vartheta}^{\alpha}-1}(\Psi^{\alpha}_{\rm
S}+\Psi^{\alpha}_{\Gamma}),$ (179)
$\displaystyle\bar{\bm{\zeta}}^{\alpha}_{\rm D}=\rho^{\alpha}{\bf
F}^{\alpha}\frac{\partial\psi^{\alpha}}{\partial\nabla^{\alpha}_{0}\bar{D}^{\alpha}}=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sqrt{\hat{G}^{\alpha}}{\bf
F}^{\alpha}\frac{\partial\hat{\Psi}^{\alpha}_{\rm
D}}{\partial\nabla_{0}\bar{D}^{\alpha}}$
$\displaystyle\quad=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sqrt{\hat{G}^{\alpha}}\bar{\Upsilon}^{\alpha}\bar{l}^{\alpha}_{\rm
R}{\bf F}^{\alpha}[(\nabla\bar{D}^{\alpha}){\bf F}^{\alpha}].$ (180)

Define total conjugate force to damage in the matrix:

$\bar{\digamma}^{\alpha}_{\rm D}=-\bar{\pi}^{\alpha}_{\rm
D}+\nabla\cdot\bar{\bm{\zeta}}^{\alpha}_{\rm D}.$ (181)

Define viscosity $\bar{\nu}^{\alpha}_{\rm D}\geq 0$. A Ginzburg-Landau kinetic
law and nonnegative dissipation for the matrix in (108) are

$\displaystyle\bar{\nu}^{\alpha}_{\rm
D}D^{\alpha}_{t}\bar{D}^{\alpha}=\bar{\digamma}^{\alpha}_{\rm D},$ (182)
$\displaystyle\bar{\mathfrak{D}}^{\alpha}_{\rm D}=\bar{\digamma}^{\alpha}_{\rm
D}D^{\alpha}_{t}\bar{D}^{\alpha}=\bar{\nu}^{\alpha}_{\rm
D}|D^{\alpha}_{t}\bar{D}^{\alpha}|^{2}\geq 0.$ (183)

To render the damage rate always nonnegative, (182) can be modified to
$\bar{\nu}^{\alpha}_{\rm
D}D^{\alpha}_{t}\bar{D}^{\alpha}=\bar{\digamma}^{\alpha}_{\rm D}{\rm
H}(\bar{\digamma}^{\alpha}_{\rm D})$. Damage kinetics are suppressed for
$\bar{\nu}^{\alpha}_{\rm D}\rightarrow\infty$ and rate-independent for
$\bar{\nu}^{\alpha}_{\rm D}\rightarrow 0$ with equilibrium condition
$\bar{\digamma}^{\alpha}_{\rm D}=0$. For rate insensitivity, irreversibility
is enforced by setting $\bar{D}^{\alpha}(t)$ to the maximum of the argument of
$\bar{\digamma}^{\alpha}_{\rm D}(\bar{D}^{\alpha};\cdot)(t)=0$ and
$\bar{D}^{\alpha}(s)\,\forall\,s<t$, where the latter renders
$D^{\alpha}_{t}\bar{D}^{\alpha}(t^{-})\rightarrow 0$. Usual, but inessential,
initial conditions are $\bar{D}^{\alpha}_{0}=0$.

Damage in fiber families $k$ is treated analogously. Viscosities are
${\nu}^{\alpha}_{{\rm D}k}\geq 0$. Conjugate thermodynamic forces and
dissipative Ginzburg-Landau kinetics are

$\displaystyle{\pi}^{\alpha}_{{\rm D}k}$
$\displaystyle=\rho^{\alpha}\frac{\partial\psi^{\alpha}}{\partial{D}^{\alpha}_{k}}=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\frac{\partial}{\partial{D}^{\alpha}_{k}}[\sqrt{\hat{G}^{\alpha}}\hat{\Psi}^{\alpha}_{\rm
D}+\varsigma^{\alpha}_{{\rm F}k}(\Psi^{\alpha}_{{\rm F}k}+\Psi^{\alpha}_{\Phi
k})]$ $\displaystyle=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}[2\sqrt{\hat{G}^{\alpha}}{E}^{\alpha}_{{\rm
C}k}{D}^{\alpha}_{k}+\Psi^{\alpha}_{\rm
D}\frac{\partial}{\partial{D}^{\alpha}_{k}}\ln\sqrt{\hat{G}^{\alpha}}]$
$\displaystyle\qquad\quad\,\,-\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}{\vartheta}^{\alpha}_{k}[1-{D}^{\alpha}_{k}]^{{\vartheta}^{\alpha}_{k}-1}(\Psi^{\alpha}_{{\rm
F}k}+\Psi^{\alpha}_{\Phi k}),$ (184) $\displaystyle{\bm{\zeta}}^{\alpha}_{{\rm
D}k}$ $\displaystyle=\rho^{\alpha}{\bf
F}^{\alpha}\frac{\partial\psi^{\alpha}}{\partial\nabla^{\alpha}_{0}{D}^{\alpha}_{k}}=\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sqrt{\hat{G}^{\alpha}}{\bf
F}^{\alpha}\frac{\partial\hat{\Psi}^{\alpha}_{\rm
D}}{\partial\nabla_{0}{D}^{\alpha}_{k}}$
$\displaystyle=2\frac{\rho^{\alpha}}{\rho^{\alpha}_{{\rm
R}0}}\sqrt{\hat{G}^{\alpha}}{\Upsilon}^{\alpha}_{k}{l}^{\alpha}_{{\rm R}k}{\bf
F}^{\alpha}[(\nabla{D}^{\alpha}_{k}){\bf F}^{\alpha}],$ (185)
$\displaystyle{\digamma}^{\alpha}_{{\rm D}k}=-{\pi}^{\alpha}_{{\rm
D}k}+\nabla\cdot{\bm{\zeta}}^{\alpha}_{{\rm D}k},$ (186)
$\displaystyle{\nu}^{\alpha}_{{\rm
D}k}D^{\alpha}_{t}{D}^{\alpha}_{k}={\digamma}^{\alpha}_{{\rm D}k},$ (187)
$\displaystyle{\mathfrak{D}}^{\alpha}_{\rm
DF}=\sum_{k}{\mathfrak{D}}^{\alpha}_{{\rm
D}k}=\sum_{k}{\digamma}^{\alpha}_{{\rm D}k}D^{\alpha}_{t}{D}^{\alpha}_{k}$
$\displaystyle\qquad\qquad\qquad\quad=\sum_{k}{\nu}^{\alpha}_{{\rm
D}k}|D^{\alpha}_{t}{D}^{\alpha}_{k}|^{2}\geq 0.$ (188)

To forbid healing, ${\nu}^{\alpha}_{{\rm
D}k}D^{\alpha}_{t}{D}^{\alpha}_{k}={\digamma}^{\alpha}_{{\rm D}k}{\rm
H}({\digamma}^{\alpha}_{{\rm D}k})$ in lieu of (187). For rate insensitivity,
${\nu}^{\alpha}_{{\rm D}k}\rightarrow 0\Rightarrow{\digamma}^{\alpha}_{{\rm
D}k}=0$ with possible irreversibility constraints analogous to those for
$\bar{D}^{\alpha}(t)$. Usual initial conditions are $D^{\alpha}_{k0}=0$.

Heat conduction. Fourier conduction is usually sufficient for each bulk
constituent $\alpha$, with isotropic conductivity
$\kappa_{\theta}^{\alpha}(\theta,\\{{\bm{\xi}}^{\alpha}\\})\geq 0$, often
temperature and internal-state dependent. It could degrade with damage via,
e.g., $\kappa^{\alpha}_{\theta}\approx\varsigma^{\alpha}_{\rm
V}\kappa^{\alpha}_{\theta 0}$. Heat flux and entropy production are

$\displaystyle{\bf
q}^{\alpha}=-\kappa^{\alpha}_{\theta}\nabla\theta^{\alpha},\quad{\mathfrak{D}}_{\rm
q}^{\alpha}=\frac{-{\bf
q}^{\alpha}\cdot\nabla\theta^{\alpha}}{\theta^{\alpha}}=\frac{\kappa^{\alpha}_{\theta}|\nabla\theta^{\alpha}|^{2}}{\theta^{\alpha}}\geq
0.$ (189)

Momentum transfer. Momentum exchange includes Darcy-like contributions from
velocity differences ${\bm{\upsilon}}^{\alpha}-{\bm{\upsilon}}^{\beta}$ =
${\bm{\mu}}^{\alpha}-{\bm{\mu}}^{\beta}$ [49, 58, 74] and mass exchange to
satisfy (62):

$\displaystyle{\bf
h}^{\alpha}=-\sum_{\beta}[\lambda^{\alpha\beta}({\bm{\mu}}^{\alpha}-{\bm{\mu}}^{\beta})]-\hat{c}^{\alpha}{\bm{\mu}}^{\alpha}.$
(190)

The inverse hydraulic-type conductivity matrix is
$\lambda^{\alpha\beta}=\lambda^{\beta\alpha}$; entries
$\lambda^{\alpha\beta}\geq 0$ can depend on temperature and volume of each
phase (e.g., to account for changes in interphase viscosity with temperature
and in permeability with porosity [49]) and degrade with damage [50, 51]:

$\displaystyle\lambda^{\alpha\beta}(J^{\alpha},J^{\beta},\theta^{\alpha},\theta^{\beta},\bar{D}^{\alpha},\bar{D}^{\beta})$
$\displaystyle\qquad\qquad=\bar{\lambda}^{\alpha\beta}(J^{\alpha},J^{\beta},\theta^{\alpha},\theta^{\beta})\sqrt{\varsigma^{\alpha}_{\rm
V}\varsigma^{\beta}_{\rm V}}.$ (191)

Energy transfer. Energy exchange includes heat transfer from temperature
differences $\theta^{\alpha}-\theta^{\beta}$ as well as momentum and mass
exchange terms to satisfy (63):

$\displaystyle\epsilon^{\alpha}=-\sum_{\beta}[\omega^{\alpha\beta}(\theta^{\alpha}-\theta^{\beta})]-m^{\alpha}\sum_{\beta}{\bf
h}^{\beta}\cdot{\bm{\mu}}^{\beta}$
$\displaystyle\qquad\qquad\qquad\quad-m^{\alpha}\sum_{\beta}\hat{c}^{\beta}(u^{\beta}+{\textstyle{\frac{1}{2}}}|{\bm{\mu}}^{\beta}|^{2}),$
(192)

recalling mass concentration $m^{\alpha}=\rho^{\alpha}/\rho$ from (83). The
matrix of heat transfer coefficients
$\omega^{\alpha\beta}=\omega^{\beta\alpha}\geq 0$ can depend on state
variables and damage [51] like (III.3):

$\displaystyle\omega^{\alpha\beta}(J^{\alpha},J^{\beta},\theta^{\alpha},\theta^{\beta},\bar{D}^{\alpha},\bar{D}^{\beta})$
$\displaystyle\qquad\qquad=\bar{\omega}^{\alpha\beta}(J^{\alpha},J^{\beta},\theta^{\alpha},\theta^{\beta})\sqrt{\varsigma^{\alpha}_{\rm
V}\varsigma^{\beta}_{\rm V}}.$ (193)

In compression, contact among fully broken constituents permits momentum and
heat transfer in respective (III.3) and (III.3) even if
$\bar{D}^{\alpha}\rightarrow 1$ or $\bar{D}^{\beta}\rightarrow 1$. Damage
dependence differing from basic illustrations (III.3) and (III.3) (e.g., [50,
51]) can be substituted if more appropriate.

Mass transfer. Terms $\hat{c}^{\alpha}$ sum to zero in (62); they account for
mass transfer rates between constituents. In biology, these could originate
from growth and remodeling, for example, exchange of nutrients dissolved in a
fluid phase to support growth of new solid tissue. Detailed constitutive
equations for $\hat{c}^{\alpha}$ are beyond the present scope. Thermodynamic
constraints emerge from (III.1) with (190) and logical stipulation that the
rightmost two sums in (III.1) should be nonnegative in concert:

$\displaystyle\sum_{\alpha}\left[\frac{\epsilon^{\alpha}}{\theta^{\alpha}}+c^{\alpha}\eta^{\alpha}\right]$
$\displaystyle=\sum_{\alpha}\sum_{\beta}\frac{\omega^{\alpha\beta}(\theta^{\alpha}-\theta^{\beta})^{2}}{2\theta^{\alpha}\theta^{\beta}}$
$\displaystyle+\sum_{\alpha}\frac{m^{\alpha}}{2\theta^{\alpha}}\sum_{\beta}\lambda^{\alpha\beta}|{\bm{\mu}}^{\alpha}-{\bm{\mu}}^{\beta}|^{2}$
$\displaystyle+\sum_{\alpha}\frac{m^{\alpha}}{\theta^{\alpha}}\sum_{\beta}\hat{c}^{\beta}\left[\frac{|{\bm{\mu}}^{\beta}|^{2}}{2}-u^{\beta}\right]$
$\displaystyle+\sum_{\alpha}\hat{c}^{\alpha}\eta^{\alpha}+\rho\eta\,\partial_{t}\ln\sqrt{g}\geq
0.$ (194)

The first two double sums in (III.3) are always nonnegative. When all
$\epsilon^{\alpha}=0$ and $\hat{c}^{\alpha}=0$ (i.e., no phase interactions),
(III.3) becomes, with (112), $\eta\,\partial_{t}\hat{g}\geq 0$. Then when
$\rho\eta=\sum_{\alpha}\rho^{\alpha}\eta^{\alpha}>0$, metric $\hat{\bf g}$
should only be dilating.

###  III.4 Stress, dissipation, and boundary conditions

Total stress ${\bm{\sigma}}^{\alpha}$ for constituent $\alpha$ is the sum of
(III.2) and (169). Total stress for the mixture ${\bm{\sigma}}$ is (57). Total
constituent dissipation $\mathfrak{D}^{\alpha}$ entering (108) is the sum of
(III.3), (173), (178), (183), (III.3), and (189), each individually
nonnegative:

$\displaystyle{\mathfrak{D}}^{\alpha}-({\bf
q}^{\alpha}\cdot\nabla\theta^{\alpha})/\theta^{\alpha}$
$\displaystyle=\hat{\mathfrak{D}}^{\alpha}+{\mathfrak{D}}^{\alpha}_{\Gamma}+{\mathfrak{D}}^{\alpha}_{\rm
A}$ $\displaystyle\quad+\bar{\mathfrak{D}}^{\alpha}_{\rm
D}+{\mathfrak{D}}^{\alpha}_{\rm DF}+{\mathfrak{D}}^{\alpha}_{\rm q}\geq 0.$
(195)

The total dissipation inequality of (III.1) is then

$\displaystyle\sum_{\alpha}\left[\frac{{\mathfrak{D}}^{\alpha}+{\mathfrak{D}}^{\alpha}_{\rm
q}}{\theta^{\alpha}}\right]+\sum_{\alpha}\left[\frac{\epsilon^{\alpha}}{\theta^{\alpha}}+c^{\alpha}\eta^{\alpha}\right]\geq
0.$ (196)

From (14), boundary conditions are required for each constituent
$\alpha=1,\ldots,N$. Mechanical conditions are prescribed histories of
traction ${\bf t}^{\alpha}$ or velocity ${\bm{\upsilon}}^{\alpha}$ on
$\partial\Omega$. Thermal conditions are histories of flux $q_{n}^{\alpha}$ or
temperature $\theta^{\alpha}$ on $\partial\Omega$. For internal variables
$\\{{\bm{\xi}}^{\alpha}\\}$ with gradient energetic dependence (e.g., order
parameters for damage), histories of fluxes $\\{{\bf
z}^{\alpha}\\}=(\rho^{\alpha}\partial\psi^{\alpha}/\partial\\{\nabla{\bm{\xi}}\\}^{\alpha})\cdot{\bf
n}$ or conjugate rates $D^{\alpha}_{t}{\bm{\xi}}^{\alpha}$ are needed on
$\partial\Omega$. Histories of body force ${\bf b}^{\alpha}$ and heat source
$r^{\alpha}$ are prescribed over $\Omega$.

###  III.5 Finsler metrics

Coordinate forms of (112) and (113) are

$\displaystyle{\bf g}({\bf x},t)=\bar{g}_{ik}({\bf
x})\,\hat{g}^{k}_{j}(\\{{\bm{\xi}}^{\alpha}({\bf x},t)\\}){\bf
g}^{i}\otimes{\bf g}^{j}$ (197) $\displaystyle{\bf G}^{\alpha}({\bf
X}^{\alpha},t)=(\bar{G}^{\alpha})_{IK}({\bf X}^{\alpha})$
$\displaystyle\qquad\qquad\times(\hat{G}^{\alpha})^{K}_{J}(\\{{\bm{\Xi}}^{\alpha}({\bf
X}^{\alpha},t)\\})({\bf G}^{\alpha})^{I}\otimes({\bf G}^{\alpha})^{J}.$ (198)

Canonical transformations [63, 70, 59] between representations of state
variables on $\mathfrak{m}$ and ${\mathfrak{M}}^{\alpha}$ are used for (87):

$\\{{\bm{\Xi}}^{\alpha}({\bf X}^{\alpha},t)\\}=\\{{\bm{\xi}}^{\alpha}({\bf
x},t)\\}\circ{\bm{\chi}}^{\alpha}({\bf X}^{\alpha},t).$ (199)

Though other relationships are admissible, the analogous transformation law
[59] between components of $\hat{\bf G}^{\alpha}$ and $\hat{\bf g}$ is
prescribed here, with $\delta^{i}_{J}$ Kronecker’s delta symbols:

$\displaystyle(\hat{G}^{\alpha})^{I}_{J}({\bf
X}^{\alpha},t)=\delta^{I}_{i}\delta^{j}_{J}\hat{g}^{i}_{j}({\bm{\chi}}^{\alpha}({\bf
X}^{\alpha},t),t),$ (200)

Only $\\{{\bm{\xi}}^{\alpha}\\}$ and
$\hat{g}^{i}_{j}(\\{{\bm{\xi}}^{\alpha}\\})$ are defined constitutively, with
(199) and (200) yielding $\\{{\bm{\Xi}}^{\alpha}\\}$ and
$(\hat{G}^{\alpha})^{I}_{J}(\\{{\bm{\Xi}}^{\alpha}\\})$, or vice-versa if
referential versions are defined instead.

Dependence of $\hat{\bf g}$ on $\\{{\bm{\xi}}^{\alpha}\\}$ is henceforth
restricted to dependence on damage parameters
$(\\{\bar{D}^{\alpha}\\},\\{D^{\alpha}_{k}\\})$; (199) is
$\bar{D}^{\alpha}({\bf
X}^{\alpha},t)=\bar{D}^{\alpha}({\bm{\chi}}^{\alpha}({\bf
X}^{\alpha},t),t)\circ{\bm{\chi}}^{\alpha}_{t}$ and so on for
$\\{D^{\alpha}_{k}\\}$. If the geometric framework is extended to describe
biologic growth and remodeling, then $\\{{\bm{\xi}}^{\alpha}\\}$ can be
expanded with internal state variable(s) associated with such processes, for
which kinetic equations in (92) are needed. For a theory as in Ref. [61],
$\hat{\bf G}^{\alpha}$ should depend on these additional (e.g., growth)
function(s), and $\hat{g}^{i}_{j}\rightarrow\delta^{i}_{j}$, so that
$\mathfrak{m}$ is Euclidean but $\mathfrak{M}^{\alpha}$ need not be.

In the current application, as tears and commensurate fiber rearrangements
arise in constituents of the mixture, the body manifold can expand and shear
[59]. Mixed-variant tensor $\hat{\bf g}$ is a product of matrix and fiber
terms:

$\displaystyle\hat{g}^{i}_{j}(\\{\bar{D}^{\alpha}\\},\\{D^{\alpha}_{k}\\})=\bar{\gamma}^{i}_{k}(\\{\bar{D}^{\alpha}\\})\tilde{\gamma}^{k}_{j}(\\{D^{\alpha}_{k}\\}).$
(201)

Contributions from isotropic matrix damage $\\{\bar{D}^{\alpha}\\}$ in
$\bar{\gamma}^{i}_{j}$ are assumed spherical (e.g., Weyl-type scaling [62]),
measured by determinants
$\bar{\gamma}^{\alpha}=\bar{\gamma}^{\alpha}(\bar{D}^{\alpha})$. Spherical
contributions of phases $\alpha=1,\ldots,N$ are merged multiplicatively in
$\hat{\bf g}$ since their sequence is irrelevant. Forms are

$\displaystyle\bar{\gamma}^{i}_{j}=\delta^{i}_{j}\prod_{\alpha}(\bar{\gamma}^{\alpha})^{1/3},\quad\bar{\gamma}^{\alpha}=\exp\left[\frac{2n_{0}^{\alpha}\bar{\kappa}^{\alpha}}{\bar{r}^{\alpha}}(\bar{D}^{\alpha})^{{\bar{r}}^{\alpha}}\right].$
(202)

Recall $n_{0}^{\alpha}({\bf X}^{\alpha})\in[0,1]$ is a reference volume
fraction of phase $\alpha$, and $\bar{r}^{\alpha}>0$ and
$\bar{\kappa}^{\alpha}$ are constants, the latter positive for dilatant
damage. Remnant volumetric strain [59] at $\bar{D}^{\alpha}=1$ is the ratio of
constants
$\bar{\epsilon}^{\alpha}=n_{0}^{\alpha}\bar{\kappa}^{\alpha}/\bar{r}^{\alpha}$.
Fiber contributions from constituents $\alpha$ and fiber families $k$ are
merged additively into $\hat{\bf g}$ since these terms are generally
anisotropic [63, 59]. Defining
$(H^{\alpha}_{k})^{i}_{j}=\delta^{i}_{I}\delta^{J}_{j}(H^{\alpha}_{k})^{I}_{J}$
with
$(H^{\alpha}_{k})^{I}_{J}=\kappa^{\alpha}_{k}\delta^{I}_{J}+(1-3\kappa^{\alpha}_{k})(\iota^{\alpha}_{k})^{I}(\iota^{\alpha}_{k})_{J}$
from (140),

$\displaystyle\tilde{\gamma}^{i}_{j}=\delta^{i}_{j}+\sum_{\alpha}\sum_{k}(H^{\alpha}_{k})^{i}_{j}\\{\exp\left[\frac{2n_{0}^{\alpha}\tilde{\kappa}^{\alpha}_{k}}{\tilde{r}_{k}^{\alpha}}({D}_{k}^{\alpha})^{\tilde{r}_{k}^{\alpha}}\right]-1\\}.$
(203)

Constants $\tilde{r}_{k}^{\alpha}>0$ and $\tilde{\kappa}_{k}^{\alpha}$ measure
the logarithmic remnant strain contributions
$\tilde{\epsilon}^{\alpha}_{k}=n_{0}^{\alpha}\tilde{\kappa}^{\alpha}_{k}/\tilde{r}^{\alpha}_{k}$
at $D^{\alpha}_{k}=1$. Noting $\hat{G}^{\alpha}=\det\hat{\bf
G}^{\alpha}=\det\hat{\bf g}=\hat{g}$ are related through (200), derivatives in
conjugate forces (III.3) and (III.3) are found from (202) and (203) as

$\displaystyle{\partial}(\ln\sqrt{\hat{G}^{\alpha}})/{\partial\bar{D}^{\alpha}}=n^{\alpha}_{0}\bar{\kappa}^{\alpha}(\bar{D}^{\alpha})^{{\bar{r}}^{\alpha}-1},$
(204)
$\displaystyle\frac{\partial(\ln\sqrt{\hat{G}^{\alpha}})}{\partial{D}^{\alpha}_{k}}=\exp\left[\frac{2n_{0}^{\alpha}\tilde{\kappa}^{\alpha}_{k}}{\tilde{r}_{k}^{\alpha}}({D}_{k}^{\alpha})^{\tilde{r}_{k}^{\alpha}}\right]$
$\displaystyle\qquad\qquad\quad\quad\times(\tilde{\gamma}^{-1})^{i}_{j}(H^{\alpha}_{k})^{j}_{i}n^{\alpha}_{0}\tilde{\kappa}^{\alpha}_{k}({D}_{k}^{\alpha})^{{\tilde{r}}^{\alpha}_{k}-1}.$
(205)

Values $\bar{r}^{\alpha}\in(0,1)$ and $\tilde{r}^{\alpha}_{k}\in(0,1)$ produce
non-singular $\hat{\bf g}$ and admit solutions to some equilibrium problems
[59]. However, these ranges can result in singularities at
$\bar{D}^{\alpha}=0$ and $D^{\alpha}_{k}=0$ in (204) and (III.5). Such
singularities can be avoided by choosing $\bar{r}^{\alpha}\geq 1$ and
$\tilde{r}^{\alpha}_{k}\geq 1$. Stronger conditions $\bar{r}^{\alpha}>1$ and
$\tilde{r}^{\alpha}_{k}>1$ usefully ensure (204), (III.5) vanish at (e.g.,
initial) states having $\bar{D}^{\alpha}=0$, ${D}_{k}^{\alpha}=0$.

##  IV  Soft-tissue physics

###  IV.1  Shock Hugoniot response

Fluids. The theory is first exercised for three fluids that comprise the
majority of soft tissues: water, extracellular fluid (ECF, representative of
blood plasma [2] and interstitial fluid in skeletal muscle and skin), and
whole blood, the latter with a realistic hematocrit of 0.4. Shock responses
are modeled via theory of Sec. II.2 with each fluid a single-phase material
($\alpha=N=1$, superscript henceforth suppressed). Planar (1-D) impact is
along the $x=x^{1}$ direction. Fluid is quiescent upstream, at density
$\rho_{0}$, temperature $\theta_{0}$, and ambient pressure $p_{0}=p_{{\rm
R}0}=1$ atm. Eulerian and Lagrangian shock speeds are identical, labeled
${\mathcal{U}}$. Particle velocity in Hugoniot states is
$\upsilon^{-}=\upsilon$. In Hugoniot (i.e., downstream shocked) $(\cdot)^{-}$
states,

$\displaystyle\rho=\rho_{0}/J,\quad J=F^{1}_{1}=F=\partial{\chi}/\partial X,$
(206) $\displaystyle u=U/\rho_{0},\quad t_{1}=\sigma^{1}_{1}=-P=-p,$ (207)

where $U$ is energy per unit reference volume and $P$ is longitudinal Cauchy
stress, positive in compression. Since volume fraction $n^{1}_{0}=n^{1}=1$,
$\rho_{{\rm R}0}=\rho_{0}$, and $\rho_{\rm R}=\rho$.

To calculate the material response, $J$ is decremented from unity. At each
decrement, the constitutive model, here the condensed matter EOS for $U_{\rm
V}$ and $p_{\rm V}$ in (III.2)–(III.2), is solved concurrently with Hugoniot
energy equation (41) for $P$ and $\eta$. With macroscopically adiabatic
conditions assumed and $\\{z\\}=0$, the latter reduces to

$\displaystyle U={\textstyle{\frac{1}{2}}}(P+p_{0})(1-J)$ (208)

since $U_{0}=0$ and $J=1$ upstream. Then $\upsilon$, $\mathcal{U}$, and
$\theta$ are found from (39) and (135). For single-phase materials,
${\bm{\mu}}^{\alpha}={\bf 0}$, ${\bf h}^{\alpha}={\bf 0}$,
$\epsilon^{\alpha}=0$, and here $\hat{c}^{\alpha}=0$. For compression, fluid
cavitation is omitted, so no damage is modeled. Hence, all metrics are
Euclidean: $g=G=\hat{g}=\hat{G}=1$.

Table 1: Physical properties or model parameters for water, extracellular fluid (ECF), human blood, and porcine skeletal muscle (cells and matrix, $\alpha=k=1$, vol. fract. $n_{0}^{1}=0.9$)  Property | Water | ECF | Human Blood | Porcine Muscle
---|---|---|---|---
$\rho^{\alpha}_{{\rm R}0}$ [g/cm3] | 1.00 | 1.03 | 1.06 | 1.10
$B^{\alpha}_{\eta}$ [GPa] | 2.10 | 2.20 | 2.64 | 3.28
$c^{\alpha}_{\epsilon}$ [J/g$\cdot$K] | 4.15 | 3.96 | 3.58 | 3.25
$\gamma_{0}^{\alpha}$ [-] | 0.120 | 0.132 | 0.160 | 0.313
$B^{\alpha}_{\eta{\rm p}}$ [-] | 6.96 | 6.96 | 12.0 | 8.0
$k^{\alpha}_{\rm V}$ [-] | 7.0 | 7.0 | 0.0 | 6.0
$\hat{B}^{\alpha}$ [mPa$\cdot$s] | 2.1 | $\ldots$ | $\ldots$ | $\ldots$
$\hat{\mu}^{\alpha}$ [mPa$\cdot$s] | 0.8 | 1.2 | 5.0 | $\ldots$
$\mu^{\alpha}_{\rm S}$ [kPa] | $\ldots$ | $\ldots$ | $\ldots$ | 1.0
$\mu^{\alpha}_{k}$ [kPa] | $\ldots$ | $\ldots$ | $\ldots$ | 1.0
$k^{\alpha}_{k}$ [-] | $\ldots$ | $\ldots$ | $\ldots$ | 10.0
$\beta^{\alpha}_{\rm S}$ [-] | $\ldots$ | $\ldots$ | $\ldots$ | $1.0\times 10^{5}$
$\beta^{\alpha}_{\Phi k}$ [-] | $\ldots$ | $\ldots$ | $\ldots$ | $1.0\times 10^{5}$
${\mathsf{J}}_{\rm C}$ [kJ/m2] | $\ldots$ | $\ldots$ | $\ldots$ | 0.84
$l^{\alpha}$ [mm] | $\ldots$ | $\ldots$ | $\ldots$ | 0.88
$\vartheta^{\alpha}$ [-] | $\ldots$ | $\ldots$ | $\ldots$ | 2.0
${\epsilon}^{\alpha}_{\rm r}$ [-] | $\ldots$ | $\ldots$ | $\ldots$ | 0.2
$\bar{r}^{\alpha}=\tilde{r}^{\alpha}_{k}$ [-] | $\ldots$ | $\ldots$ | $\ldots$ | 2.0

Properties entering the EOS are listed in Table 1. All are obtained or
estimated from experimental literature [104, 2, 105, 106, 107, 108] with the
exception of nonlinear bulk stiffening parameters $B_{\eta{\rm p}}$ and
$k_{\rm V}$ that are fit to the experimental Hugoniot data. The ambient bulk
modulus is related to the bulk sound velocity $c_{\rm
B}=\sqrt{(B_{\eta}+p_{0})/\rho_{0}}$. For water, the usual relationship
$B_{\eta{\rm p}}=4S-1$ is used, where $S$ is the slope of a linear fit to the
$\mathcal{U}$-$\upsilon$ Hugoniot [88]. Since shock compression data are
apparently unavailable for ECF, $B_{\eta{\rm p}}$ and $k_{\rm V}$ for water
are assigned. Newtonian viscosities $\hat{B},\hat{\mu}$ are listed for
completeness [2, 105], where $\hat{\mu}$ for blood is for high rates [2] and
bulk viscosity is supplied only for water [105]. These do not enter the
present analysis: Newtonian viscosity (169) and Fourier conduction (189) are
necessarily omitted as both are incompatible with treatment of shocks as
singular surfaces [109, *morro1980b].

Compared in Fig. 1(a) are $P$-$\rho$ Hugoniot predictions and experimental
data on water ($\theta_{0}=297$ K)[106] and human blood ($\theta_{0}\approx
310$ K) [111]. Shock data on ECF do not seem to exist; predictions are for
$\theta_{0}=310$ K. Compared in Fig. 1(b) are $\mathcal{U}$-$\upsilon$ data
for water [106] and predictions for all three fluids; shock velocity data were
not reported for blood in Ref. [111]. The $c_{\rm B}$ value for blood,
${\mathcal{U}}\rightarrow c_{\rm B}$ as $\upsilon\rightarrow 0$, is from Ref.
[108]. Results in Fig. 1 confirm validity of the EOS for water and blood. The
latter is stiffer than ECF, which is stiffer than water.

(

Figure 1:  Model results and shock data [106, 111] for (a) Hugoniot stress vs.
mass density and (b) shock velocity vs. particle velocity in water,
extracellular fluid (ECF), and human blood.

a)
(b)

Skeletal muscle. Planar shock response of skeletal muscle is predicted next.
Properties and loading conditions replicate impact experiments of Wilgeroth et
al. [112, 5] on porcine muscle tissue. The material is modeled as a mixture of
two coexisting phases ($\alpha=1,2$): a “solid” tissue of initial volume
fraction $n_{0}^{1}=0.9$ [113] and an interstitial fluid depicted by the ECF,
comprising remaining fraction $n_{0}^{2}=0.1$. The first phase consists of the
muscle cells (i.e., fibers), collagenous connective tissues, and ground
substance between and encasing the cells (i.e., the extracellular matrix).
Muscle cells contain significant internal fluid whose physical properties are
included implicitly in the constitutive model for the first phase.

Experiments [112] show a single-wave structure with steep shock front (rise
time on the order of $\mu{\rm s}$ or smaller) rather than multiple waveforms
that would be expected if shock and particle velocities among the phases
differed significantly [73, 74]. This suggests inverse hydraulic conductivity
is very large (e.g., $\lambda^{\alpha\beta}\tilde{\rightarrow}\infty$) at
these high-pressure dynamic conditions, and diffusion velocities
${\bm{\mu}}^{\alpha}$ are negligible. It is thus assumed particle velocity
histories match in each phase:
$\upsilon^{\alpha}(x,t)=\upsilon^{1}(x,t)=\upsilon^{2}(x,t)=\upsilon(x,t)$.
Therefore, $J(x,t)=F(x,t)=\partial\chi/\partial X$ is identical in each
constituent. Microsecond scales are too brief for biologic mass exchange:
$\hat{c}^{\alpha}=0$. Shock compression is adiabatic: $q^{\alpha}\rightarrow
0$). Heat transfer in $\epsilon^{\alpha}$ of (III.3) is likewise assumed null
in Hugoniot states: $\theta^{\alpha}=\theta^{1}=\theta^{2}=\theta$.

The solution procedure is similar to that for single fluids, but the
constitutive model is now much more complex. The shock response is that of the
mixture, where governing and jump equations are in (68)–(76). Both phases are
quiescent and at reference $\theta_{0}=310$ K and $p_{{\rm R}0}=1$ atm
upstream; $\theta_{0}$ was unreported in Ref. [112], but model results are
insensitive to $\theta_{0}$. From (57) and (58) with ${\bm{\mu}}^{\alpha}={\bf
0}$, mixture stress and internal energy are

$\displaystyle\bm{\sigma}=\sum_{\alpha}{\bm{\sigma}}^{\alpha},\quad
U=\sum_{\alpha}n_{0}^{\alpha}\rho^{\alpha}_{{\rm
R}0}u^{\alpha}=\sum_{\alpha}n^{\alpha}_{0}U^{\alpha}.$ (209)

Because all phases deform equally without mass exchange, mixture density is
$\rho=\rho_{0}/J$. As explained later in the context of (213), $\\{{\bf
z}^{\alpha}\\}\rightarrow\\{z^{\alpha}\\}=0$. The analog of (41) for the
mixture reduces to (208) with $U_{0}=0$ and
$p_{0}=\sum_{\alpha}n^{\alpha}_{0}p^{\alpha}_{{\rm R}0}$. In calculations, $J$
is reduced incrementally from unity. In each decrement, energy equation (208)
and constitutive equations for each phase are solved simultaneously and
summed, if appropriate, to give mixture values $P$, $U$, and $\theta$. Given
$\theta$ and $J$, entropy $\eta^{\alpha}$ of each phase is found by inversion
of (135). Particle and shock velocities are found from mixture analogs of
(39). The response of the fluid phase (ECF, $\alpha=2$) is calculated as
before; its energy and pressure contributions are given fully by
$U^{\alpha}_{\rm V}$ and $p^{\alpha}_{\rm V}$.

The tissue phase (including intracellular fluids), $\alpha=1$, has a total
internal energy per unit reference volume $U^{\alpha}$:

$\displaystyle U^{1}=U_{\rm V}^{1}+\varsigma^{1}_{\rm S}\cdot(U^{1}_{\rm
S}+U^{1}_{\Gamma})+\varsigma^{1}_{\rm F}\circ(U^{1}_{\rm
F}+U^{1}_{\Phi})+U^{1}_{\rm D}.$ (210)

The first term on the right is the EOS (noting $\varsigma_{\rm V}=1$ for
compression), second and third are deviatoric matrix elasticity and
viscoelasticity, third and fourth are deviatoric fiber elasticity and
viscoelasticity, and the last is surface energy of soft-tissue degradation
(i.e., damage). Only passive states are modeled: $U^{1}_{\rm A}=0$. Thermal
variables $\theta^{1}$ and $\eta^{1}$ only enter $U_{\rm V}^{1}$, which fully
specifies the partial pressure $p^{1}$. Notation $U$ and $\Psi$ is
interchangeable for remaining terms that only affect deviatoric response.

Matrix viscoelasticity is limited to the shear response, following typical
assumptions for nearly incompressible soft materials [97, 98, 19]. For very
rapid loading modeled here, viscous relaxation for all ($m$) configurational
variables ${\bm{\Gamma}}^{1}_{{\rm S}m}$ is assumed negligible with
$t/\tau^{1}_{{\rm S}m}\rightarrow 0$, so

$\displaystyle U^{1}_{\rm S}+U^{1}_{\Gamma}\rightarrow U^{1}_{\rm
S}+\sum_{m}\hat{\Psi}^{1}_{{\rm S}m}=U^{1}_{\rm S}(1+\sum_{m}{\beta}^{1}_{{\rm
S}m}).$ (211)

Thus, $\check{\mu}^{1}_{\rm S}=\mu^{1}_{\rm S}(1+\sum_{m}\beta^{1}_{{\rm
S}m})$ is the glassy shear modulus of the matrix, with static energy and
modulus in (138).

Muscle fibers comprise one family, $k=1$, of direction
${\bm{\iota}}^{\alpha}_{k}={\bm{\iota}}$ with $\kappa^{1}_{1}=0$ in (140).
Viscous relaxation for all ($n$) configurational variables
${\bm{\Gamma}}^{1}_{\Phi k,n}$ is assumed negligible:

$\displaystyle U^{1}_{\rm F}+U^{1}_{\Phi}\rightarrow U^{1}_{\rm
F}+\sum_{n}\hat{\Psi}^{1}_{{\Phi}1,n}=U^{1}_{\rm
F}(1+\sum_{n}{\beta}^{1}_{{\Phi}1,n}).$ (212)

Static strain energy of fibers $U^{1}_{\rm F}$ is (III.2) with
$\mu^{\alpha}_{k}=\mu^{1}_{1}$ and ${\rm H}(\cdot)$ omitted, supporting
compressive stress. A dynamic fiber modulus is
$\check{\mu}^{1}_{1}={\mu}^{1}_{1}(1+\sum_{n}{\beta}^{1}_{{\Phi}1,n})$. Fiber
directions relative to $x=x^{1}$ are ambiguous in Ref. [112]. Calculations
apply loading parallel or transverse to ${\bm{\iota}}$, both pure mode
directions. In the former, the longitudinal sound speed ${\mathcal{C}}$ obeys
$\rho_{0}{\mathcal{C}}^{2}=B_{\eta}+p_{0}+\frac{4}{3}n_{0}^{1}(\check{\mu}_{\rm
S}^{1}+\frac{2}{3}\check{\mu}^{1}_{1})$. In the latter,
$\rho_{0}{\mathcal{C}}^{2}=B_{\eta}+p_{0}+\frac{4}{3}n_{0}^{1}(\check{\mu}_{\rm
S}^{1}+\frac{1}{6}\check{\mu}^{1}_{1})$.

Damage order parameters for the matrix,
$\bar{D}^{\alpha}=\bar{D}^{1}=\bar{D}\in[0,1]$, and fibers,
$D^{\alpha}_{k}=D^{1}_{1}=D_{1}\in[0,1]$, degrade respective deviatoric stress
contributions from matrix and fibers in (III.2). They also supply surface
energy $U^{1}_{\rm D}=\Psi^{1}_{\rm D}$ in (III.2). Jumps in $\bar{D}$ and
$D_{k}$ are allowed across the shock front. This necessitates
$\bar{\alpha}^{\alpha}=\alpha^{\alpha}_{k}=0\Rightarrow\bar{l}^{\alpha}_{\rm
R}={l}^{\alpha}_{{\rm R}k}=0$ in (III.2) to avoid infinite energy in the
front. Gradient energies and conjugate forces $\bar{\bm{\zeta}}^{\alpha}_{\rm
D},{\bm{\zeta}}^{\alpha}_{{\rm D}k}$ vanish identically, as do $\\{{\bf
z}^{\alpha}\\}$. Treatment of shocks as singular surfaces mandates viscosities
$\bar{\nu}^{\alpha}_{\rm D}=\nu^{\alpha}_{{\rm D}k}=0$ for fracture to avoid
infinite dissipation in the shock front if $\bar{D}$ and $D_{k}$ are