Balance laws

Elasticity (single-phase)

Mathematical models for continuum mechanics rely upon governing laws that ensure the balance of mass and momentum in the system being modeled. Below we provide the equations for balance of mass and linear momentum used in the material methods section.

Lagrangian view

Law	Formulation
Balance of Mass	\(\rho J - \rho_0 = 0\)
Balance of Linear Momentum	\(\rho_0 \frac{\partial}{\partial t} \bm{V} - \nabla_X \cdot \bm{P} - \rho_0 \bm{g} = 0\)

\(\bm{P}\) is the first Piola-Kirchhoff stress tensor, \(\bm{V}\) is the Lagrangian velocity defined in (3), and \(\bm{g}\) is the gravitational acceleration. \(\rho\) and \(\rho_0\) denote the mass density and initial mass density, respectively.

Eulerian view

Law	Formulation
Balance of Mass	\(\frac{\partial}{\partial t} \rho + \nabla \cdot (\rho \bm{v}) = 0\)
Balance of Linear Momentum	\(\rho \frac{D}{D t} \bm{v} - \nabla \cdot \bm{\sigma} - \rho \bm{g} = 0\)

where \(\bm{v}\) is the Eulerian velocity (4).

Note

The balance of linear momentum can be written in the equivalent form

(26)\[ \frac{\partial}{\partial t} \left(\rho \bm{v} \right) + \nabla \cdot \left( \rho \bm{v} \otimes \bm{v} \right) - \nabla \cdot \bm{\sigma} - \rho \bm{g} = 0. \]

by expanding \(\rho \frac{D}{D t} \bm{v}\) using (6) as

\[ \begin{aligned} \rho \frac{D}{D t} \bm{v} &= \rho \frac{\partial \bm{v}}{\partial t} + \rho \left( \nabla \bm{v} \right) \bm{v} = \frac{\partial \left( \rho \bm{v} \right)}{\partial t} - \frac{\partial \rho}{\partial t} \bm{v} + \left( \nabla \bm{v} \right) \left(\rho \bm{v} \right) \\ &= \frac{\partial \left( \rho \bm{v} \right)}{\partial t} + \bm{v} \nabla \cdot \left( \rho \bm{v} \right) + \left( \nabla \bm{v} \right) \left(\rho \bm{v} \right) \\ &= \frac{\partial \left( \rho \bm{v} \right)}{\partial t} + \nabla \cdot \left( \rho \bm{v} \otimes \bm{v} \right) \end{aligned} \]

where we have used the identity \(\nabla \cdot \left( \bm{u} \otimes \bm{v} \right) = \bm{u} \nabla \cdot \bm{v} + \left( \nabla \bm{u} \right) \bm{v}\).

Poromechanics (multiphase)

General concepts

For an infinitesimal volume element \(dv\), we assume a biphasic mixture comprising of a fluid (f) and solid (s). Herein, the particles of the fluid phase and solid phase are assumed to exist together. In other words, \(dv\) is a superimposed continua of constituents \(\varphi^\alpha\), which is defined as the sum of the partial volume elements \(dv_\rs\) and \(dv_\rf\),

(27)\[ dv = \sum_\alpha dv_\alpha = dv_\rs + dv_\rf\, . \]

The volume fractions of the fluid phase \(\phi^\rf\) and solid phase \(\phi^\rs\) are defined via

(28)\[ \phi^\rf = \frac{dv_\rf}{dv}, \quad \phi^\rs = \frac{dv_\rs}{dv}\, . \]

The saturation condition is obtained by a combination of above equations as

(29)\[ \sum_\alpha \phi^\alpha = \phi^\rs + \phi^\rf = 1 \]

Note that \(\phi^\rf\) is the commonly used concept porosity \(\phi\) in soild mechanics. If it’s not mentioned herein, porosity \(\phi\) is used to represent the volume fraction of the fluid phase, and \(1-\phi\) to represent the volume fraction of the solid phase.

Also, we define the real mass density and partial mass density of constituent \(\varphi^\alpha\) by

\[ \rho^{\alpha{\rm{R}}} = \frac{dm_{\alpha}}{dv_{\alpha}}, \quad \rho^{\alpha} = \frac{dm_{\alpha}}{dv} = \phi^{\alpha}\rho^{\alpha{\rm{R}}}, \]

A variety of terminologies have been introduced in engineering practice to describe the velocities of the solid and the fluid phase of the mixture. Fluid velocities are usually expressed relative to the motion of the solid phase. For example, the relative velocity vector of the fluid with respect to the solid skeleton motion (true seepage velocity) is

(30)\[ \tilde{\bm v}_\rf = \bm{v}_\rf - \bm{v}_\rs = \bm{v}_\rf - \bm v \]

The other more commonly used relative velocity is the superficial (Darcy seepage) fluid velocity given by

(31)\[ \dot{\bm w} = \phi^\rf \tilde{\bm v}_\rf = \phi^\rf \left( \bm{v}_\rf - \bm v \right) \]

Darcy velocity represents the relative volumetric rate of discharge per unit area of the fluid-solid mixture.

Linear theory

For linear theory (Biot theory), the porosity \(\phi^\rf\) is assumed to include disconnected pores, i.e., porous regions within the solid itself [SD03]. Only if the pores are all interconnected would we recover the non-linear theory (Non-linear theory (finite strain)). Furthermore, we assume \(\phi^\rf\) is constant from which it follows that for a closed system, constituent densities \(\rho^{\alpha\rm{R}}\) are also constant, and therefore current mixture density \(\rho^b\) is also constant.

The continuity equation of the fluid phase is given by [Che16, DCC13]

(32)\[ \dot{\zeta} + \nabla\cdot \bm{q} = \dot{\gamma}, \]

where \(\zeta\) is the variation in fluid content, \(\gamma\) is a fluid volumetric source and \(\bm q = \dot{\bm w}\) is the flux vector. Integrating (32) with respect to time (assumed zero initial value) we obtain

(33)\[ \zeta = -\nabla \cdot \bm w + \gamma. \]

The balance of linear momentum for the linear mixture theory is given by

(34)\[ \nabla \cdot \bm \sigma + \bm f^\rs = \rho^b \ddot{\bm u} + \rho^{\fR} \ddot{\bm w}, \]

where

(35)\[ \rho^b = (1 - \phi^\rf) \rho^{\sR} + \phi^\rf \rho^{\fR}, \]

is the bulk density (mixture density) of the porous medium, which we take to be constant for constant volume fractions and real mass densitities of \(\varphi^\alpha\). For the fluid phase, the balance of linear momentum, or generalized Darcy’s law, is

(36)\[ \bm q = \dot{\bm w} = - \varkappa \left(\nabla p_\rf + \rho^{\fR} \ddot {\bm u} + \rho^w \ddot{\bm w} - \bm f^\rf \right), \]

where \(\bm f^\rf\) is the body force acting on the pore fluid (e.g., gravity),

\[ \varkappa = \frac{\varkappa_{\circ}}{\mu_\rf}, \]

is the permeability coefficient (assuming isotropic permeability), with \(\varkappa_{\circ}\) the intrinsic permeability (units m\(^2\)), and \(\mu_f\) the dynamic shear viscosity (units Pa-s) of the fluid phase, and

\[ \rho^w = \frac{\rho^a}{(\phi^\rf)^2} + \frac{\rho^{\fR}}{\phi^\rf}, \]

\(\rho^a\) is apparent mass density. In some references, \(\rho^w\) is defined as

\[ \rho^w = \alpha_{\infty} \frac{\rho^{\fR}}{\phi^\rf}, \]

where \(\alpha_{\infty}\) is the tortuosity of the solid skeleton (solid frame).

Non-linear theory (finite strain)

Balance of Mass

To derive the balance of mass for porous material at large deformation, we start by defining differential mass, real mass and partial mass densities of constituent \(\varphi^\alpha\) as

(37)\[ \begin{aligned} m_{\alpha} &= \int_{\Omega^{\alpha}} dm_{\alpha} = \int_{\Omega^{\alpha}} \rho^{\alpha{\rm{R}}}dv_{\alpha} = \int_{\Omega} \rho^{\alpha{\rm{R}}}\phi^{\alpha} dv = \int_{\Omega} \rho^{\alpha} dv = \int_{\Omega^{\alpha}_0} \rho^{\alpha} J_{\alpha} dV_{\alpha}, \\ \rho^{\alpha{\rm{R}}} &= \frac{dm_{\alpha}}{dv_{\alpha}}, \\ \rho^{\alpha} &= \frac{dm_{\alpha}}{dv} = \phi^{\alpha}\rho^{\alpha{\rm{R}}}, \end{aligned} \]

where we have used (28). Then, the material time derivative and balance of mass is

(38)\[ \begin{aligned} \frac{D^{\alpha} m_{\alpha}}{D t} &= \int_{\Omega^{\alpha}_0} \frac{D^{\alpha} \rho^{\alpha} J_{\alpha}}{D t} dV_{\alpha} = \int_{\Omega} \hat{\rho}^{\alpha} dv, \\ &= \int_{\Omega^{\alpha}_0} \left(\frac{D^{\alpha} \rho^{\alpha}}{D t} J_{\alpha} + \rho^{\alpha} \frac{D^{\alpha} J_{\alpha}}{D t} \right)dV_{\alpha} = \int_{\Omega} \hat{\rho}^{\alpha} dv, \\ &= \int_{\Omega} \left(\frac{D^{\alpha} \rho^{\alpha}}{D t} + \rho^{\alpha} \nabla_x\cdot \bm{v}_{\alpha} \right)dv = \int_{\Omega} \hat{\rho}^{\alpha} dv, \\ &\Longrightarrow \frac{D^{\alpha} \rho^{\alpha}}{D t} + \rho^{\alpha} \nabla_x\cdot \bm{v}_{\alpha} = \hat{\rho}^{\alpha}, \end{aligned} \]

where we localized the integral at the last step, \(\hat{\rho}^{\alpha}\) is the mass supply of \(\rho^{\alpha}\) per unit total current volume, and we have used

(39)\[ \frac{D^{\alpha} J_{\alpha}}{D t} = J_{\alpha} \trace \left( \bm{F}_{\alpha}^{-1} \frac{D^{\alpha}{\bm{F}}_{\alpha}}{D t}\right) = J_{\alpha} \trace \left(\frac{\partial \bm{X}_{\alpha}}{\partial \bm{x}} \frac{\partial \bm{v}_{\alpha}}{\partial \bm{X}_{\alpha}}\right) = J_{\alpha} \nabla_x \cdot \bm{v}_{\alpha}. \]

Using material time derivative of phase \(\alpha\) (13), we can write the derivative of \(\rho^{\alpha}\) solely in terms of the solid phase motion as

(40)\[ \frac{D^{\alpha} \rho^{\alpha}}{D t} = \frac{D^\rs \rho^{\alpha}}{D t} + \tilde{\bm v}_{\alpha} \cdot \nabla_x \rho^{\alpha}, \]

where \(\tilde{\bm v}_{\alpha}\) is the relative velocity vector of constituent \(\varphi^\alpha\) with respect to the solid phase motion (see (30)). Writing (38) for solid and fluid phase s, f, and using \(\rho^\rs = \phi^\rs\rho^\sR, \rho^\rf = \phi^\rf\rho^\fR\) we arrive at

\[ \begin{aligned} \frac{\phi^\rs}{\rho^{\sR}}\frac{D^{\rs} \rho^{\sR}}{D t} + \frac{D^{\rs} \phi^{\rs}}{D t} + \phi^{\rs} \nabla_x\cdot \bm{v}_{\rs} &= \frac{\hat{\rho}^{\rs}}{\rho^{\sR}}, \\ \frac{\phi^\rf}{\rho^{\fR}}\frac{D^{\rf} \rho^{\fR}}{D t} + \frac{D^{\rf} \phi^{\rf}}{D t} + \phi^{\rf} \nabla_x\cdot \bm{v}_{\rf} &= \frac{\hat{\rho}^{\rf}}{\rho^{\fR}}, \end{aligned} \]

Using barotropic flow relation (i.e., assume no thermomechanical coupling)

(41)\[ \frac{D^{\alpha} p_{\alpha}}{D t} = \frac{\bulk_{\alpha}}{\rho^{\alpha{\rm{R}}}}\frac{D^{\alpha} \rho^{\alpha{\rm{R}}}}{D t}, \]

replacing \( D^\rf(\cdot)/Dt\) by \(D^\rs(\cdot)/D t + \nabla_x(\cdot)\cdot\tilde{\bm{v}}_\alpha\) and \(\bm{v}_\rf = \tilde{\bm v}_\rf + \bm{v}_\rs\), and summing over both phases, we get

(42)\[ \frac{\phi^\rs}{\bulk_\rs}\frac{D^{\rs} p_\rs}{D t} + \frac{\phi^\rf}{\bulk_\rf}\frac{D^{\rs} p_\rf}{D t} + \nabla_x\cdot \bm{v}_{\rs} + \nabla_x\cdot (\phi^\rf \tilde{\bm v}_\rf) + \frac{\phi^\rf \tilde{\bm v}_\rf}{\bulk_\rf}\cdot \nabla_x p_\rf = \frac{\hat{\rho}^{\rs}}{\rho^{\sR}} + \frac{\hat{\rho}^{\rf}}{\rho^{\fR}} \]

Assuming the solid phase is nearly incompressible (i.e., \(D^\rs \rho^\sR/Dt \rightarrow 0\) such that \(D^\rs p_\rs/Dt = 0\)) and that there is no supply of solid mass, and using the definition (31), we get simplified balance of mass for the solid-fluid mixture as

(43)\[ \frac{\phi^\rf}{\bulk_\rf}\frac{D^{\rs} p_\rf}{D t} + \nabla_x\cdot \bm{v}_{\rs} + \nabla_x\cdot (\dot{\bm w}) + \frac{\dot{\bm w}}{\bulk_\rf}\cdot \nabla_x p_\rf = \frac{\hat{\rho}^{\rf}}{\rho^{\fR}}. \]

Note that (43) is written in current configuration, we can derive the balance of mass in initial configuration as

(44)\[ \frac{J_\rs \phi^\rf}{\bulk_\rf}\frac{D^{\rs} p_\rf}{D t} + \frac{D^{\rs} J_\rs}{D t} + J_\rs \nabla_{X_\rs}(\dot{\bm w})\tcolon \bm{F}_\rs^{-T} + \frac{J_\rs \dot{\bm w}}{\bulk_\rf} \cdot \bm{F}_\rs^{-T}\nabla_{X_\rs} p_\rf = \frac{J_\rs\hat{\rho}^{\rf}}{\rho^{\fR}}. \]

From the assumption of mechanical incompressiblity of the solid phase (equivalent to \(\bulk_\rs \rightarrow \infty\)), we get

(45)\[ D^\rs \rho^{\sR}/D t =0 \Longrightarrow \rho^{\sR} = \rho^{\sR}_0, \]

and from (38) for \(\alpha = \rs\) and (39) we can derive, via integration in time,

(46)\[ \phi^\rs = \phi^\rs_0/ J_\rs, \]

where \(\phi^\rs_0\) is the initial value of \(\phi^\rs\) at \(t = 0\). Furthermore, in the generalized Darcy’s law

(47)\[ \dot{\bm{w}} = \phi^\rf\tilde{\bm{v}}_\rf = -\varkappa(\nabla_x p_\rf + \rho^\fR(\ddot{\bm{u}}_\rf - \bm{f}^\rf))\, , \]

the permeability coefficient can be a function of porosity:

\[ \varkappa = \frac{\varkappa_{\circ}}{\mu_\rf} \frac{\mathcal{F}(\phi^\rf)}{\mathcal{F}(\phi^\rf_0)}, \]

where \(\mathcal{F}\) is a nonlinear function of porosity \(\phi^\rf\) accounting for change in permeability due to change in porosity. E.g., the Kozeny-Carman relation is

(48)\[ \mathcal{F}(\phi^\rf) = \frac{(\phi^\rf)^3}{1 - (\phi^\rf)^2}. \]

Note that in the \((\bm{u}\text{-}p_\rf)\) formulation derived later, we assume \(\ddot{\bm{u}}_\rf \approx \ddot{\bm{u}}_\rs = \ddot{\bm{u}}\) (accelerations of \(\varphi^\alpha\) are equivalent) such that the generalized Darcy’s law may be written as

(49)\[ \dot{\bm{w}} = \phi^\rf\tilde{\bm{v}}_\rf = -\varkappa(\nabla_x p_\rf + \rho^\fR(\ddot{\bm{u}} - \bm{f}^\rf))\, . \]

Balance of Linear Momentum

The balance of linear momentum for \(\varphi^\alpha\) may be written in localized form as

(50)\[ \nabla_x\cdot\bm{\sigma}^\alpha + \rho^{\alpha}\bm{f}^\alpha + \bm{h}^\alpha = \rho^\alpha\ddot{\bm{u}}_\alpha + \hat{\rho}^\alpha\dot{\bm{u}}_\alpha\, , \]

where \(\bm{\sigma}^\alpha\) is the partial Cauchy stress, such that the total Cauchy stress \(\bm{\sigma} = \bm{\sigma}^\rs + \bm{\sigma}^\rf\), \(\hat{\rho}^\alpha\dot{\bm{u}}_\alpha\) is the mass supply momentum (which we assume is negligible), and \(\bm{h}^\alpha\) is the interaction body force from all other constituents on \(\varphi^\alpha\). Usually, the interaction body forces between constituents are due to drag and will sum to zero because they are equal and opposite. Thus, these forces do not affect the mixture as a whole such that \(\bm{h}^\rs + \bm{h}^\rf = \bm{0}\). Via the balance of angular momentum it can be shown that the partial Cauchy stresses for each \(\varphi^\alpha\) are symmetric, i.e., \(\bm{\sigma}^\alpha = (\bm{\sigma}^\alpha)^T\).

(50) can be mapped back to the initial configuration of the solid skeleton for \(\varphi^\alpha\) as

(51)\[ \nabla_{X_\rs}\cdot\bm{P}_\rs^\alpha + \rho_{0(\rs)}^\alpha \bm{f}^\alpha + J_\rs\bm{h}^\alpha = \rho_{0(\rs)}^\alpha\ddot{\bm{u}}_\alpha + \hat{\rho}^\alpha_{0(\rs)}\dot{\bm{u}}_\alpha\, . \]

From (51) we can derive the balance of linear momentum equation for the biphasic mixture in the initial configuration \(\Omega_0^\rs\) using the following equations and assumptions:

• Total Cauchy stress, and first Piola stress with respect to \(\Omega_0^\rs\):

  \[ \bm{\sigma} = \bm{\sigma}^\rs + \bm{\sigma}^\rf \ , \ \bm{P}_\rs = \bm{P}_\rs^\rs + \bm{P}_\rs^\rf\, , \]

• Decomposition of the solid partial Cauchy stress into an effective (extra) stress, hereafter referred to as the effective stress (solid skeleton extra stress, or solid extra stress):

  (52)\[ \bm{\sigma}^\rs = \bm{\sigma}^\rs_E - \phi^\rs p_\rf \bm{I}\, , \]

• Assume all mass supplies are negligible: \(\hat{\rho}^\alpha = 0\),

• Assume body forces per unit mass are only due to gravity: \(\bm{f}^\alpha = \bm{g}\), where \(\bm{g}\) is the acceleration vector of gravity.

We may then write the balance of linear momentum equation for the biphasic mixture in the initial configuration \(\Omega_0^\rs\) as

(53)\[ \begin{aligned} &\nabla_{X_\rs}\cdot\bm{P}_\rs + \rho^b_{0(\rs)}\bm{g} = \rho^\rs_{0(\rs)}\ddot{\bm{u}}_\rs + \rho^\rf_{0(\rs)}\ddot{\bm{u}}_\rf\, ,\\ &\bm{P}_\rs = \bm{P}^\rs_{E(s)} + \bm{P}^\rf_{E(s)} - J_\rs p_\rf\bm{F}_\rs^{-T}\, , \end{aligned} \]

with \(\rho^b\) defined in (35) for evolving volume fractions and pore fluid real mass density. For the \((\bm{u}\text{-}p_\rf)\) formulation shown later, we will further assume that \(\ddot{\bm{u}}_\rf \approx \ddot{\bm{u}}_\rs = \ddot{\bm{u}}\), which, in theory, should be valid only for slower dynamic loadings. With that assumption in mind, the pore fluid viscous stress tensor \(\bm{P}_{E(s)}^\rf\) drops from (54)\(_2\) because the constitutive law that defines it in includes spatial dependence on pore fluid velocity, which cannot be determined without a third governing equation for the pore fluid displacement. In this sense, the pore fluid is nearly-inviscid, since drag due to pore fluid viscosity is accounted for via Darcy’s law.

Thus the balance of linear momentum of the mixture is written as

(54)\[ \begin{aligned} &\nabla_{X_\rs}\cdot\bm{P}_\rs + \rho^b_{0(\rs)}\bm{g} = \rho^b_{0(\rs)}\ddot{\bm{u}}\, ,\\ &\bm{P}_\rs = \bm{P}^\rs_{E(s)} - J_\rs p_\rf\bm{F}_\rs^{-T}\, . \end{aligned} \]

Biot theory at finite strain

Recall the Biot modulus is defined as

\[ B \:= 1 - \frac{k_d}{k_\rs}\, , \]

such that we may define the total Cauchy stress of the mixture in terms of an effective stress (of the solid skeleton/drained mixture) and a pore fluid pressure component,

\[ \bm{\sigma} = \bm{\sigma}_E^\rs - B p_\rf\bm{I}\, . \]

For a nearly-inviscid pore fluid with \(\bm{\sigma}^\rf = -p^\rf\bm{I} = \phi^\rf p_\rf\bm{I}\), this implies that the solid stress in (52) is actually

(55)\[ \begin{aligned} &\bm{\sigma}^\rs = \bm{\sigma} - \bm{\sigma}^\rf = \bm{\sigma}_E^\rs - p_\rf(B - \phi^\rf)\bm{I}\, ,\\ &\Rightarrow \bm{\sigma}_E^\rs = \bm{\sigma}^\rs + p_\rf(B - \phi^\rf)\bm{I}\, . \end{aligned} \]

Such a formulation introduces unnecessary complexity when deriving constitutive equations for the solid and fluid phases, which will be shown as follows.

For ease of notation, we will denote Eulerian descriptions of the constituent velocity and acceleration as follows, respectively:

\[ \frac{\partial\bm{\chi}_\alpha(\bm{X}_\alpha, t)}{\partial t} = \bm{v}_\alpha\, ,\quad \frac{\partial^2\bm{\chi}_\alpha(\bm{X}_\alpha, t)}{\partial t^2} = \bm{a}_\alpha\, . \]

Starting from the dissipation inequality (refer to [Ehl02, Irw24]) for \(\varphi^\alpha\),

\[ \sum_\alpha \frac{1}{\theta^\alpha}\Big(\rho^\alpha\Big[\frac{D^\alpha\psi^\alpha}{Dt} + \eta^\alpha\frac{D^\alpha\theta^\alpha}{Dt}\Big] - \bm{\sigma}^\alpha:\bm{l}_\alpha - \hat{e}^\alpha + \hat{\rho}^\alpha\Big[\psi^\alpha - \frac{1}{2}\bm{v}_\alpha\cdot\bm{v}_\alpha\Big] + \bm{h}^\alpha\cdot\bm{v}_\alpha\\ + \frac{1}{\theta^\alpha}\nabla_x(\theta^\alpha)\cdot\bm{q}^\alpha\Big) \geq 0\, , \]

where for \(\varphi^\alpha\), \(\theta^\alpha\) is the temperature (units K), \(\eta^\alpha\) is the entropy per unit mass (units J/kg-K), \(\hat{e}^\alpha\) is the power supply to \(\varphi^\alpha\) from other phases, and \(\bm{q}^\alpha\) is the heat flux vector, we will assume locally homogeneous temperatures such that \(\theta^\alpha = \theta\) and no mass supplies (\(\hat{\rho}^\alpha = 0\)).

Then, the dissipation inequality is written as (note since the constituent temperature \(\theta^\alpha = \theta\), it can be brought outside the sum and both sides can be multiplied by \(\theta\) to eliminate the scaling by temperature)

\[ \sum_\alpha \Big(\rho^\alpha\Big[\frac{D^\alpha\psi^\alpha}{Dt} + \eta^\alpha\frac{D^\alpha\theta}{Dt}\Big] - \bm{\sigma}^\alpha :\bm{l}_\alpha - \hat{e}^\alpha + \bm{h}^\alpha\cdot\bm{v}_\alpha + \frac{1}{\theta}\nabla_x(\theta)\cdot\bm{q}^\alpha\Big) \geq 0\, . \]

In deriving constitutive equations for individual constituents, we begin by writing the dissipation inequality summed over our two constituents/phases (solid and fluid):

\[ \rho^\rs\Big[\frac{D^\rs\psi^\rs}{Dt} + \eta^\rs\frac{D^\rs\theta}{Dt}\Big] - \bm{\sigma}^\rs : \bm{l}_\rs - \hat{e}^\rs + \bm{h}^\rs\cdot\bm{v}_\rs\\ + \rho^\rf\Big[\frac{D^\rf\psi^\rf}{Dt} + \eta^\rf\frac{D^\rf\theta}{Dt}\Big] - \bm{\sigma}^\rf : \bm{l}_\rf - \hat{e}^\rf + \bm{h}^\rf\cdot\bm{v}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, . \]

For a closed system, we must have

\[ \sum_\alpha\bm{h}^\alpha = \bm{0},\, \quad \sum_\alpha\hat{e}^\alpha = 0\, , \]

and we know from balance of angular momentum for \(\varphi^\alpha\) that the partial Cauchy stresses are symmetric (\(\bm{\sigma}^\alpha = (\bm{\sigma}^\alpha)^T\)). Thus we may futher simplify the dissipation inequality of the mixture:

(56)\[ \rho^\rs\Big[\frac{D^\rs\psi^\rs}{Dt} + \eta^\rs\frac{D^\rs\theta}{Dt}\Big] - \bm{\sigma}^\rs : \bm{d}_\rs + \bm{h}^\rf\cdot\tilde{\bm{v}}_\rf\\ + \rho^\rf\Big[\frac{D^\rf\psi^\rf}{Dt} + \eta^\rf\frac{D^\rf\theta}{Dt}\Big] - \bm{\sigma}^\rf : \bm{d}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, . \]

Substituting in the assumed solid stress from (55)\(_1\) and the pore fluid stress, we have

\[ \begin{aligned} \rho^\rs\Big[\frac{D^\rs\psi^\rs}{Dt} + \eta^\rs\frac{D^\rs\theta}{Dt}\Big] - \bm{\sigma}_E^\rs : \bm{d}_\rs + p_\rf(B - \phi^\rf)\nabla_x\cdot\bm{v}_\rs + \bm{h}^\rf\cdot\tilde{\bm{v}}_\rf\\ + \rho^\rf\Big[\frac{D^\rf\psi^\rf}{Dt} + \eta^\rf\frac{D^\rf\theta}{Dt}\Big] + \phi^\rf p_\rf\nabla_x\cdot\bm{v}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, . \end{aligned} \]

This can be simplified further:

\[ \begin{aligned} \rho^\rs\Big[\frac{D^\rs\psi^\rs}{Dt} + \eta^\rs\frac{D^\rs\theta}{Dt}\Big] - \bm{\sigma}_E^\rs : \bm{d}_\rs + Bp_\rf\nabla_x\cdot\bm{v}_\rs + \phi^\rf p_\rf\nabla_x\cdot{\tilde{\bm{v}}}_\rf + \bm{h}^\rf\cdot\tilde{\bm{v}}_\rf\\ + \rho^\rf\Big[\frac{D^\rf\psi^\rf}{Dt} + \eta^\rf\frac{D^\rf\theta}{Dt}\Big] + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, . \end{aligned} \]

Following the approach of [dB05], the following response functions \(\mathcal{R}\) must be determined:

\[ \mathcal{R} := \lbrace \psi^\rs, \psi^\rf, \bm{\sigma}_E^\rs, \bm{q}^\rs, \bm{q}^\rf\rbrace \]

which depend on a set of \(\mathcal{S}\) variables, i.e.,

\[ \mathcal{R}(\mathcal{S})\, , \]

where \(\mathcal{S}\) is a to-be-determined subset of the fundamental constitutive variables \(\mathcal{V}\) for a biphasic continuum with an elastic solid constituent, i.e.,

\[ \mathcal{V} := \lbrace \theta^\alpha, \nabla_x\theta^\alpha, \phi^\rf, \nabla_x\phi^\rf, \rho^{\alpha\rm{R}}, \bm{F}_\alpha, \nabla_{X_\alpha}\bm{F}_\alpha, \nabla_x\cdot\bm{v}_\rs, \bm{v}_\rf, \nabla_{X_\rf}\bm{v}_\rf, \bm{X}_\alpha\rbrace\, . \]

For an isotropic pore fluid, it can be shown [Cro73] that the deformation gradient \(\bm{F}_\rf\) is a function of \({\rm{det}}(\bm{F}_\rf)\), which can be determined through porosity \(\phi^\rf\) and real fluid mass density \(\rho^\fR\) in the absence of fluid mass supply \(\hat{\rho}^\rf\). Using the principle of frame indifference, we may also substitute the pore fluid velocity \(\bm{v}_\rf\) by the seepage velocity \(\tilde{\bm{v}}_\rf\). The pore fluid velocity gradient \(\nabla_{X_\rf}\bm{v}_\rf\) may also be substituted by the deformation rate tensor \(\bm{d}_\rf\).

Assuming an incompressible solid, it can be shown that \(\phi^\rs = \phi^\rs_{0(\rs)}({\rm{det}}\bm{F}_\rs)^{-1}\) and thus the real mass density of the solid and its gradient, and porosity and its gradient, may be eliminated from the set of respionse functions. We have also assumed \(\theta^\alpha = \theta\). Thus the constitutive variables pertaining to the problem at hand is

\[ \mathcal{S} = \lbrace \theta, \nabla_x\theta, \rho^\fR, \nabla_x\rho^\fR, \bm{F}_\rs, \nabla_{X_\rs}\bm{F}_\rs, \nabla_x\cdot\bm{v}_\rs, \tilde{\bm{v}}_\rf, \bm{d}_\rf\rbrace\, . \]

Thus from \(\mathcal{R} = \mathcal{R}(\mathcal{S})\) we have

\[ \mathcal{R} = \lbrace \psi^\rs, \psi^\rf, \bm{\sigma}_E^\rs, \bm{q}^\rs, \bm{q}^\rf\rbrace = \mathcal{R}(\theta, \nabla_x\theta, \rho^\fR, \nabla_x\rho^\fR, \bm{F}_\rs, \nabla_{X_\rs}\bm{F}_\rs, \nabla_x\cdot\bm{v}_\rs, \tilde{\bm{v}}_\rf, \bm{d}_\rf)\, . \]

Based on the principle of phase separation, that the Helmholtz free energy of \(\varphi^\alpha\) should depend only on the \(\varphi^\alpha\) variables (in the case of locally homogeneous temperatures, we can assume cross-dependence of temperature), we may write

\[ \begin{aligned} &\psi^\rs := \psi^\rs(\theta, \nabla_x\theta, \bm{C}_\rs, \nabla_{X_\rs}\bm{C}_\rs, \nabla_x\cdot\bm{v}_\rs)\, ,\\ &\psi^\rf := \psi^\rf(\theta, \nabla_x\theta, \rho^\fR, \nabla_x\rho^\fR, \tilde{\bm{v}}_\rf, \bm{d}_\rf)\, , \end{aligned} \]

where we have used the right Cauchy-Green tensor in place of the deformation gradient (of the solid). Then, the material time derivatives of the Helmholtz free energy functions are as follows:

\[ \begin{array}{ll} &\begin{aligned} \frac{D^\rs\psi^\rs}{Dt} = \frac{\partial\psi^\rs}{\partial\theta}\frac{D^\rs\theta}{Dt} + \frac{\partial\psi^\rs}{\partial(\nabla_x\theta)}\frac{D^\rs(\nabla_x\theta)}{Dt} + \frac{\partial\psi^\rs}{\partial\bm{C}_\rs}\frac{D^\rs\bm{C}_\rs}{Dt} + \frac{\partial\psi^\rs}{\partial(\nabla_{X_\rs}\bm{C}_\rs)}\frac{D^\rs(\nabla_{X_\rs}\bm{C}_\rs)}{Dt}\\ + \frac{\partial\psi^\rs}{\partial(\nabla_x\cdot\bm{v}_\rs)}\frac{D^\rs(\nabla_x\cdot\bm{v}_\rs)}{Dt}\, , \end{aligned}\\ &\begin{aligned} \frac{D^\rf\psi^\rf}{Dt} = \frac{\partial\psi^\rf}{\partial\theta}\frac{D^\rf\theta}{Dt} + \frac{\partial\psi^\rf}{\partial(\nabla_x\theta)}\frac{D^\rf(\nabla_x\theta)}{Dt} +\frac{\partial\psi^\rf}{\partial\rho^\fR}\frac{D^\rf\rho^\fR}{Dt} + \frac{\partial\psi^\rf}{\partial(\nabla_x\rho^\fR)}\frac{D^\rf(\nabla_x\rho^\fR)}{Dt}\\ + \frac{\partial\psi^\rf}{\partial\tilde{\bm{v}}_\rf}\frac{D^\rf\tilde{\bm{v}}_\rf}{Dt} + \frac{\partial\psi^\rf}{\partial\bm{d}_\rf}\frac{D^\rf\bm{d}_\rf}{Dt}\, . \end{aligned} \end{array} \]

Substitution into the dissipation inequality of the mixture yields

\[ \begin{aligned} \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_x\theta)}\frac{D^\rs(\nabla_x\theta)}{Dt} + \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_{X_\rs}\bm{C}_\rs)}\frac{D^\rs(\nabla_{X_\rs}\bm{C}_\rs)}{Dt} + \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_x\cdot\bm{v}_\rs)}\frac{D^\rs(\nabla_x\cdot\bm{v}_\rs)}{Dt} + \rho^\rs\Big[\frac{\partial\psi^\rs}{\partial\theta}+ \eta^\rs\Big]\frac{D^\rs\theta}{Dt}\\ + \Big[\rho^\rs\frac{\partial\psi^\rs}{\partial\bm{C}_\rs} - \frac{1}{2J_\rs}\bm{S}_{E(\rs)}^\rs\Big]\frac{D^\rs\bm{C}_\rs}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial(\nabla_x\theta)}\frac{D^\rf(\nabla_x\theta)}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial(\nabla_x\rho^\fR)}\frac{D^\rf(\nabla_x\rho^\fR)}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial\tilde{\bm{v}}_\rf}\frac{D^\rf\tilde{\bm{v}}_\rf}{Dt}\\ + \rho^\rf\frac{\partial\psi^\rf}{\partial\bm{d}_\rf}\frac{D^\rf\bm{d}_\rf}{Dt} + \rho^\rf\Big[\frac{\partial\psi^\rf}{\partial\theta} + \eta^\rf\Big]\frac{D^\rf\theta}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial\rho^\fR}\frac{D^\rf\rho^\fR}{Dt}\\ + Bp_\rf\nabla_x\cdot\bm{v}_\rs + \phi^\rf p_\rf\nabla_x\cdot{\tilde{\bm{v}}}_\rf + \bm{h}^\rf\cdot\tilde{\bm{v}}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, . \end{aligned} \]

For the velocity divergence terms, we will attempt to simplify using the balance of mass of the mixture,

\[ \begin{aligned} &\phi^\rs\nabla_x\cdot\bm{v}_\rs + \frac{\phi^\rf}{\rho^\fR}\frac{D^\rf\rho^\fR}{Dt} + \phi^\rf\nabla_x\cdot\bm{v}_\rf + \nabla_x(\phi^\rf)\cdot\tilde{\bm{v}}_\rf = 0\\ &\Rightarrow \nabla_x\cdot\bm{v}_\rs + \frac{\phi^\rf}{\rho^\fR}\frac{D^\rf\rho^\fR}{Dt} + \nabla_x\cdot(\phi^\rf\tilde{\bm{v}}_\rf) = 0\\ &\Rightarrow Bp_\rf\nabla_x\cdot\bm{v}_\rs = -\frac{B\phi^\rf p_\rf}{\rho^\fR}\frac{D^\rf\rho^\fR}{Dt} - Bp_\rf\nabla_x\cdot(\phi^\rf\tilde{\bm{v}}_\rf)\, , \end{aligned} \]

such that the dissipation inequality of the mixture becomes

\[ \begin{aligned} \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_x\theta)}\frac{D^\rs(\nabla_x\theta)}{Dt} + \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_{X_\rs}\bm{C}_\rs)}\frac{D^\rs(\nabla_{X_\rs}\bm{C}_\rs)}{Dt} + \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_x\cdot\bm{v}_\rs)}\frac{D^\rs(\nabla_x\cdot\bm{v}_\rs)}{Dt} + \rho^\rs\Big[\frac{\partial\psi^\rs}{\partial\theta}+ \eta^\rs\Big]\frac{D^\rs\theta}{Dt}\\ + \Big[\rho^\rs\frac{\partial\psi^\rs}{\partial\bm{C}_\rs} - \frac{1}{2J_\rs}\bm{S}_{E(\rs)}^\rs\Big]\frac{D^\rs\bm{C}_\rs}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial(\nabla_x\theta)}\frac{D^\rf(\nabla_x\theta)}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial(\nabla_x\rho^\fR)}\frac{D^\rf(\nabla_x\rho^\fR)}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial\tilde{\bm{v}}_\rf}\frac{D^\rf\tilde{\bm{v}}_\rf}{Dt}\\ + \rho^\rf\frac{\partial\psi^\rf}{\partial\bm{d}_\rf}\frac{D^\rf\bm{d}_\rf}{Dt} + \rho^\rf\Big[\frac{\partial\psi^\rf}{\partial\theta} + \eta^\rf\Big]\frac{D^\rf\theta}{Dt} + \rho^\rf\Big[\frac{\partial\psi^\rf}{\partial\rho^\fR} - \frac{B p_\rf}{(\rho^\fR)^2}\Big]\frac{D^\rf\rho^\fR}{Dt}\\ + (1 - B)\phi^\rf p_\rf\nabla_x\cdot{\tilde{\bm{v}}}_\rf + (\bm{h}^\rf - B p_\rf\nabla_x(\phi^\rf))\cdot\tilde{\bm{v}}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, . \end{aligned} \]

Using the [CN63] argument, i.e., that the free parameters \(D^\rs\theta/Dt\), \(D^\rs(\nabla_x\theta)/Dt\), \(D^\rs(\nabla_{X_\rs}\bm{C}_\rs)/Dt\), \(D^\rs(\nabla_x\cdot\bm{v}_\rs)/Dt\), \(D^\rf\theta/Dt\), \(D^\rf(\nabla_x\theta^\rf)/Dt\), \(D^\rf\rho^\fR/Dt\), \(D^\rf(\nabla_x\rho^\fR)/Dt\), \(D^\rf\tilde{\bm{v}}_\rf/Dt\) and \(D^\rf\bm{d}_\rf/Dt\) maintain arbitrary values, the following constitutive relations must hold:

(57)\[ \begin{alignedat}{3} \Big(\rho^\rs\frac{\partial\psi^\rs}{\partial\theta} + \rho^\rs\eta^\rs\Big) &= 0 && \Rightarrow \rho_{0(\rs)}^\rs\eta^\rs = -\frac{\partial(\rho_{0(\rs)}^\rs\psi^\rs)}{\partial\theta}\, ,\\ \frac{\partial\psi^\rs}{\partial(\nabla_x\theta)} &= 0\, , &&\phantom{\Rightarrow}\\ \Big(\rho^\rs\frac{\partial\psi^\rs}{\partial\bm{C}_\rs}\bm{C}_\rs - \frac{1}{2J_\rs}\bm{S}^\rs_{E(\rs)}\Big) &= \bm{0} &&\Rightarrow \bm{S}^\rs_{E(\rs)} = 2\frac{\partial(\rho^\rs_{0(\rs)}\psi^\rs)}{\partial\bm{C}_\rs}\, ,\\ \frac{\partial\psi^\rs}{\partial(\nabla_{X_\rs}\bm{C}_\rs)} &= \bm{0}\, , &&\phantom{\Rightarrow}\\ \frac{\partial\psi^\rs}{\partial(\nabla_x\cdot{\bm{v}}_\rs)} &= \bm{0}\, , &&\phantom{\Rightarrow}\\ \Big(\rho^\rf\frac{\partial\psi^\rf}{\partial\theta} + \rho^\rf\eta^\rf\Big) &= 0 &&\Rightarrow \eta^\rf = -\frac{\partial\psi^\rf}{\partial\theta}\, ,\\ \frac{\partial\psi^\rf}{\partial(\nabla_x\theta)} &= 0\, , &&\phantom{\Rightarrow}\\ \Big(\rho^\rf\frac{\partial\psi^\rf}{\partial\rho^\fR} - \frac{B\phi^\rf}{\rho^\fR}\Big) &= 0 &&\Rightarrow p_\rf = \frac{(\rho^\fR)^2}{B}\frac{\partial\psi^\rf}{\partial\rho^\fR}\, ,\\ \frac{\partial\psi^\rf}{\partial(\nabla_x\rho^\fR)} &= 0\, , &&\phantom{\Rightarrow}\\ \frac{\partial\psi^\rf}{\partial\tilde{\bm{v}}_\rf} &= \bm{0}\, , &&\phantom{\Rightarrow}\\ \frac{\partial\psi^\rf}{\partial\bm{d}_\rf} &= \bm{0}\, , &&\phantom{\Rightarrow} \end{alignedat} \]

where we have made use of the identity

\[ \rho^\rs\frac{D^\rs\psi^\rs}{Dt} = \frac{1}{J_\rs}\frac{D^\rs\left(\rho^\rs_{0(\rs)}\psi^\rs\right)}{Dt}\, . \]

Examining (57)\(_8\), we know from thermodynamic principles of a compressible pore fluid [Cla11, Cla19, Dav08] that

(58)\[ \frac{\partial\psi^\rf}{\partial v^\rf} = -p_\rf \Rightarrow \frac{\partial\psi^\rf}{\partial\rho^\fR}\frac{\partial\rho^\fR}{\partial v^\rf} = -p_\rf \Rightarrow (\rho^\fR)^2\frac{\partial\psi^\rf}{\partial\rho^\fR} = p_\rf\, , \]

and therefore we are already in violation of thermodynamic admissibility from the result in (57)\(_8\).

Proceeding anyway, we are left with the reduced dissipation inequality

\[ \mathcal{D} := (1 - B)\phi^\rf p_\rf\nabla_x\cdot{\tilde{\bm{v}}}_\rf + (\bm{h}^\rf - B p_\rf\nabla_x(\phi^\rf))\cdot\tilde{\bm{v}}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, , \]

to determine the remaining constitutive laws. By inspection, it is evident that the heat flux terms can be conjugate to the mixture temperature gradient by assuming a Fourier’s law, i.e.,

\[ \bm{q}^\alpha = -\bm{k}^{\theta^\alpha}\nabla_x\theta\, , \]

where \(\bm{k}^{\theta^\alpha}\) is a materially-dependent isotropic (or anisotropic) heat conduction tensor for \(\varphi^\alpha\).

It is trickier to define thermodynamic conjugacy for the seepage velocity terms. For the second term, let

\[ \bm{h}^\rf_B = \bm{h}^\rf - B p_\rf\nabla_x(\phi^\rf) = \bm{h}^\rf - p_\rf\nabla_x(\phi^\rf) + p_\rf\nabla_x(\phi^\rf)(1 - B)\, . \]

In order to define thermodynamic conjugacy for \(\bm{h}^\rf_B\) (a drag-like term), it should be proportional to \(\tilde{\bm{v}}_\rf\), i.e.,

\[ \bm{h}^\rf_B = -\bm{S}_w\tilde{\bm{v}}_\rf\, , \]

where the permeability tensor \(\bm{S}_w := (\phi^\rf)^2(\bm{K}^\rs)^{-1}\) and \(\bm{K}^\rs\) is the intrinsic permeability tensor of the drained mixture (solid skeleton), such that for isotropic permeability \(\bm{K}^\rs := \varkappa\bm{I}\).

We will use the balance of linear momentum for the pore fluid, assuming \(\bm{\sigma}^\rf = -\phi^\rf p_\rf\bm{I}\), in place of \(\bm{h}^\rf\) above:

\[ \begin{aligned} &\nabla_x\cdot\bm{\sigma}^\rf + \rho^\rf\bm{f}^\rf + \bm{h}^\rf = \rho^\rf\bm{a}_\rf &\Rightarrow \bm{h}^\rf = \rho^\rf(\bm{a}_\rf - \bm{f}^\rf) + \phi^\rf\nabla_x(p_\rf)) + p_\rf\nabla_x(\phi^\rf)\, . \end{aligned} \]

Then,

\[ \bm{h}^\rf_B = \rho^\rf(\bm{a}_\rf - \bm{f}^\rf) + \phi^\rf\nabla_x(p_\rf) + p_\rf\nabla_x(\phi^\rf)(1-B)\, . \]

Using the definition \(\bm{h}^\rf_B = -\bm{S}_w\tilde{\bm{v}}_\rf\) for isotropic permeability, we have

(59)\[ -\varkappa(\rho^\fR[\bm{a}_\rf - \bm{f}^\rf] + \nabla_x(p_\rf) + p_\rf\nabla_x(\phi^\rf)[1-B]) = \phi^\rf\tilde{\bm{v}}_\rf = \dot{\bm{w}}\, . \]

This is akin to (47), except that the \((1-B)\) term is left over. The remaining term in \(\mathcal{D}\),

\[ (1 - B)\phi^\rf p_\rf\nabla_x\cdot{\tilde{\bm{v}}}_\rf\, , \]

cannot maintain thermodynamic conjugacy since \(B, \phi^\rf, p_\rf\) are already defined. Only in the limit \(B \rightarrow 1\) is the model thermodynamically admissible (and do we recover Darcy’s law from (59)).

If instead we begin to derive constituive models when assuming nothing of the form of the Cauchy stresses, then, following the approach of [Ehl02], we introduce the saturation constraint:

\[ \phi^\rs + \phi^\rf = 1, \quad \frac{D^\rs\phi^\rs}{Dt} + \frac{D^\rs\phi^\rf}{Dt} = 0\, , \]

which may also be expressed as

\[ \frac{D^\rs\phi^\rs}{Dt} + \frac{D^\rf\phi^\rf}{Dt} - \nabla_x(\phi^\rf)\cdot\tilde{\bm{v}}_\rf = 0\, , \]

and with utilization of the balance of mass of the mixture for incompressible solid, it may also be written as

\[ \phi^\rs\nabla_x\cdot\bm{v}_\rs + \frac{\phi^\rf}{\rho^\fR}\frac{D^\rf\rho^\fR}{Dt} + \phi^\rf\nabla_x\cdot\bm{v}_\rf + \nabla_x(\phi^\rf)\cdot\tilde{\bm{v}}_\rf = 0\, . \]

The saturation constraint may be added to the dissipation inequality of the mixture with the introduction of a Lagrange multiplier \(\Lambda\):

\[ \Lambda\Big(\phi^\rs\nabla_x\cdot\bm{v}_\rs + \frac{\phi^\rf}{\rho^\fR}\frac{D^\rf\rho^\fR}{Dt} + \phi^\rf\nabla_x\cdot\bm{v}_\rf + \nabla_x(\phi^\rf)\cdot\tilde{\bm{v}}_\rf\Big) = 0\, , \]

such that the dissipation inequality of the mixture,

\[ \rho^\rs\Big[\frac{D^\rs\psi^\rs}{Dt} + \eta^\rs\frac{D^\rs\theta}{Dt}\Big] - \bm{\sigma}^\rs : \bm{d}_\rs + \bm{h}^\rf\cdot\tilde{\bm{v}}_\rf\\ + \rho^\rf\Big[\frac{D^\rf\psi^\rf}{Dt} + \eta^\rf\frac{D^\rf\theta}{Dt}\Big] - \bm{\sigma}^\rf : \bm{d}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, , \]

can be written as

\[ \begin{aligned} \rho^\rs\Big[\frac{D^\rs\psi^\rs}{Dt} + \eta^\rs\frac{D^\rs\theta}{Dt}\Big] - [\bm{\sigma}^\rs + \Lambda\phi^\rs\bm{I}] : \bm{d}_\rs + [\bm{h}^\rf - \Lambda\nabla_x(\phi^\rf)]\cdot\tilde{\bm{v}}_\rf + \rho^\rf\Big[\frac{D^\rf\psi^\rf}{Dt} + \eta^\rf\frac{D^\rf\theta}{Dt}\Big]\\ - [\bm{\sigma}^\rf + \Lambda\phi^\rf\bm{I}] : \bm{d}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) - \Lambda\frac{\phi^\rf}{\rho^\fR}\frac{D^\rf \rho^\fR}{Dt}\geq 0\, . \end{aligned} \]

Introduce the following terms for ease of notation:

(60)\[ \begin{aligned} &\bm{\sigma}^\rs_E := \bm{\sigma}^\rs + \Lambda\phi^\rs\bm{I}\, ,\\ &\bm{\sigma}^\rf_E := \bm{\sigma}^\rf + \Lambda\phi^\rf\bm{I}\, ,\\ &\bm{h}^\rf_E := \bm{h}^\rf - \Lambda\nabla_x(\phi^\rf)\, ,\\ \end{aligned} \]

and substitute back into the dissipation inequality:

\[ \begin{aligned} \rho^\rs\Big[\frac{D^\rs\psi^\rs}{Dt} + \eta^\rs\frac{D^\rs\theta}{Dt}\Big] + \rho^\rf\Big[\frac{D^\rf\psi^\rf}{Dt} + \eta^\rf\frac{D^\rf\theta}{Dt}\Big] - \bm{\sigma}^\rs_E : \bm{d}_\rs - \bm{\sigma}^\rf_E : \bm{d}_\rf + \bm{h}^\rf_E\cdot\tilde{\bm{v}}_\rf\\ + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) - \Lambda\frac{\phi^\rf}{\rho^\fR}\frac{D^\rf \rho^\fR}{Dt}\geq 0\, . \end{aligned} \]

Proceeding as before, denote the response function

\[ \mathcal{R} := \lbrace \psi^\rs, \psi^\rf, \bm{\sigma}_E^\rs, \bm{\sigma}_E^\rf, \bm{h}^\rf_E, \bm{q}^\rs, \bm{q}^\rf\rbrace\, . \]

From the same arguments for the independent set of constitutive variables \(\mathcal{S}\), we have now

\[ \mathcal{R} = \lbrace \psi^\rs, \psi^\rf, \bm{\sigma}_E^\rs, \bm{\sigma}_E^\rf, \bm{h}^\rf_E, \bm{q}^\rs, \bm{q}^\rf\rbrace = \mathcal{R}(\theta, \nabla_x\theta, \rho^\fR, \nabla_x\rho^\fR, \bm{F}_\rs, \nabla_{X_\rs}\bm{F}_\rs, \nabla_x\cdot\bm{v}_\rs, \tilde{\bm{v}}_\rf, \bm{d}_\rf)\, . \]

Then the dissipation inequality may be written as follows:

\[ \begin{aligned} \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_x\theta)}\frac{D^\rs(\nabla_x\theta)}{Dt} + \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_{X_\rs}\bm{C}_\rs)}\frac{D^\rs(\nabla_{X_\rs}\bm{C}_\rs)}{Dt} + \rho^\rs\frac{\partial\psi^\rs}{\partial(\nabla_x\cdot\bm{v}_\rs)}\frac{D^\rs(\nabla_x\cdot\bm{v}_\rs)}{Dt} + \rho^\rs\Big[\frac{\partial\psi^\rs}{\partial\theta} + \eta^\rs\Big]\frac{D^\rs\theta}{Dt}\\ + \Big[\rho^\rs\frac{\partial\psi^\rs}{\partial\bm{C}_\rs} - \frac{1}{2J_\rs}\bm{S}_{E(\rs)}^\rs\Big]\frac{D^\rs\bm{C}_\rs}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial(\nabla_x\theta)}\frac{D^\rf(\nabla_x\theta)}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial(\nabla_x\rho^\fR)}\frac{D^\rf(\nabla_x\rho^\fR)}{Dt} + \rho^\rf\frac{\partial\psi^\rf}{\partial\tilde{\bm{v}}_\rf}\frac{D^\rf\tilde{\bm{v}}_\rf}{Dt}\\ + \rho^\rf\frac{\partial\psi^\rf}{\partial\bm{d}_\rf}\frac{D^\rf\bm{d}_\rf}{Dt} + \rho^\rf\Big[\frac{\partial\psi^\rf}{\partial\theta} + \eta^\rf\Big]\frac{D^\rf\theta}{Dt} + \rho^\rf\Big[\frac{\partial\psi^\rf}{\partial\rho^\fR} - \frac{\Lambda}{(\rho^\fR)^2}\Big]\frac{D^\rf\rho^\fR}{Dt}\\ - \bm{\sigma}^\rf_E : \bm{d}_\rf + \bm{h}^\rf_E\cdot\tilde{\bm{v}}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, . \end{aligned} \]

The [CN63] argument states that

(61)\[ \begin{alignedat}{3} \Big(\rho^\rs\frac{\partial\psi^\rs}{\partial\theta} + \rho^\rs\eta^\rs\Big) &= 0 && \Rightarrow \rho_{0(\rs)}^\rs\eta^\rs = -\frac{\partial(\rho_{0(\rs)}^\rs\psi^\rs)}{\partial\theta}\, ,\\ \frac{\partial\psi^\rs}{\partial(\nabla_x\theta)} &= 0\, , &&\phantom{\Rightarrow}\\ \Big(\rho^\rs\frac{\partial\psi^\rs}{\partial\bm{C}_\rs}\bm{C}_\rs - \frac{1}{2J_\rs}\bm{S}^\rs_{E(\rs)}\Big) &= \bm{0} &&\Rightarrow \bm{S}^\rs_{E(\rs)} = 2\frac{\partial(\rho^\rs_{0(\rs)}\psi^\rs)}{\partial\bm{C}_\rs}\, ,\\ \frac{\partial\psi^\rs}{\partial(\nabla_{X_\rs}\bm{C}_\rs)} &= \bm{0}\, , &&\phantom{\Rightarrow}\\ \frac{\partial\psi^\rs}{\partial(\nabla_x\cdot{\bm{v}}_\rs)} &= \bm{0}\, , &&\phantom{\Rightarrow}\\ \Big(\rho^\rf\frac{\partial\psi^\rf}{\partial\theta} + \rho^\rf\eta^\rf\Big) &= 0 &&\Rightarrow \eta^\rf = -\frac{\partial\psi^\rf}{\partial\theta}\, ,\\ \frac{\partial\psi^\rf}{\partial(\nabla_x\theta)} &= 0\, , &&\phantom{\Rightarrow}\\ \Big(\rho^\rf\frac{\partial\psi^\rf}{\partial\rho^\fR} - \frac{\Lambda}{(\rho^\fR)^2}\Big) &= 0 &&\Rightarrow \Lambda = (\rho^\fR)^2\frac{\partial\psi^\rf}{\partial\rho^\fR}\, ,\\ \frac{\partial\psi^\rf}{\partial(\nabla_x\rho^\fR)} &= 0\, , &&\phantom{\Rightarrow}\\ \frac{\partial\psi^\rf}{\partial\tilde{\bm{v}}_\rf} &= \bm{0}\, , &&\phantom{\Rightarrow}\\ \frac{\partial\psi^\rf}{\partial\bm{d}_\rf} &= \bm{0}\, , &&\phantom{\Rightarrow} \end{alignedat} \]

where

(62)\[ \mathcal{D} := - \bm{\sigma}^\rf_E : \bm{d}_\rf + \bm{h}^\rf_E\cdot\tilde{\bm{v}}_\rf + \frac{1}{\theta}\nabla_x(\theta)\cdot(\bm{q}^\rs + \bm{q}^\rf) \geq 0\, . \]

From (58) and (61)\(_8\) it follows that \(\Lambda = p_\rf\). Thus, from (60),

(63)\[ \begin{aligned} &\bm{\sigma}^\rs_E = \bm{\sigma}^\rs + \phi^\rs p_\rf\bm{I}\, ,\\ &\bm{\sigma}^\rf_E = \bm{\sigma}^\rf + \phi^\rf p_\rf\bm{I}\, ,\\ &\bm{h}^\rf_E = \bm{h}^\rf - p_\rf\nabla_x(\phi^\rf)\, . \end{aligned} \]

It is then clear that (note the lack of \(B\) now)

\[ \bm{\sigma} = \bm{\sigma}^\rs_E + \bm{\sigma}^\rf_E - p_\rf\bm{I}\, . \]

From (61)\(_3\), it is apparent that \(\bm{\sigma}^\rs_E\) is defined from the Helmholtz free energy function of the solid skeleton/drained mixture, in other words, \(\bm{\sigma}^\rs_E = \bm{\sigma}^\rs_E(\theta, \bm{C}_\rs)\) and not a function of the Biot parameter.

By inspecting the \(\bm{\sigma}^\rf_E\) term in (62), it is reasonable to assume thermodynamic conjugacy through a fourther order tensor \(\overset{4}{\bm{Z}}_\rf\) such that

\[ \bm{\sigma}_E^\rf = \overset{4}{\bm{Z}}_\rf\bm{d}_\rf\, . \]

A common choice for \(\overset{4}{\bm{Z}}_\rf\) is a simple Newtonian fluid law (see [Hol00] p. 203), i.e.,

\[ \overset{4}{\bm{Z}}_\rf := \phi^\rf\kappa_\rf(\bm{I} \otimes \bm{I})^{\overset{23}{T}} + 2\phi^\rf\mu_\rf(\bm{I} \otimes \bm{I})\, , \]

where \(\kappa_\rf\) and \(\mu_\rf\) are the bulk and shear dynamic viscosities of the pore fluid, respectively. For a nearly-inviscid pore fluid, we take \(\bm{\sigma}^\rf_E = \bm{0}\). Determing conjugacy for \(\bm{q}^\alpha\) proceeds as before.

For conjugacy of \(\bm{h}^\rf_E\), the derivation is similar to that for \(\bm{h}^\rf_B\), and it should be proportional to \(\tilde{\bm{v}}_\rf\), i.e.,

\[ \bm{h}^\rf_E = \bm{h}^\rf - p_\rf\nabla_x(\phi^\rf) = -\bm{S}_w\tilde{\bm{v}}_\rf\, , \]

where the permeability tensor \(\bm{S}_w := (\phi^\rf)^2(\bm{K}^\rs)^{-1}\) and \(\bm{K}^\rs\) is the intrinsic permeability tensor of the drained mixture (solid skeleton), such that for isotropic permeability \(\bm{K}^\rs := \varkappa\bm{I}\).

We will use the balance of linear momentum for the pore fluid, assuming \(\bm{\sigma}^\rf = -\phi^\rf p_\rf\bm{I}\), in place of \(\bm{h}^\rf\) above:

(64)\[ \begin{aligned} &\nabla_x\cdot\bm{\sigma}^\rf + \rho^\rf\bm{f}^\rf + \bm{h}^\rf = \rho^\rf\bm{a}_\rf &\Rightarrow \bm{h}^\rf = \rho^\rf(\bm{a}_\rf - \bm{f}^\rf) + \phi^\rf\nabla_x(p_\rf) + p_\rf\nabla_x(\phi^\rf)\, . \end{aligned} \]

Then the porosity gradient terms cancel, such that,

(65)\[ \bm{h}^\rf_E = \rho^\rf(\bm{a}_\rf - \bm{f}^\rf) + \phi^\rf\nabla_x(p_\rf)\, . \]

Using the definition \(\bm{h}^\rf_E = -\bm{S}_w\tilde{\bm{v}}_\rf\) for isotropic permeability, we have

(66)\[ -\varkappa(\rho^\fR[\bm{a}_\rf - \bm{f}^\rf] + \nabla_x(p_\rf)) = \phi^\rf\tilde{\bm{v}}_\rf = \dot{\bm{w}}\, , \]

where we have recovered (47).

In conclusion, the Biot theory is inappropriate for rigorous application of finite strain theory as it does not “proceed from the basic equations of the single constituents but directly from an equivalent equation for the whole body – solid skeleton and water.” [dBE88]