Balance laws

Elasticity

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}\).

Linear poroelasticity

For an infinitesimal volume element \(dv\), the particles of the fluid phase and solid phase are assumed to exist together. In other words, \(dv\) is a superimposed continua, which is defined as the sum of the partial volume elements \(dv_s\) and \(dv_f\),

(27)\[ dv = dv_s + dv_f \]

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

(28)\[ \phi^f = \frac{dv_f}{dv}, \quad \phi^s = \frac{dv_s}{dv} \]

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

(29)\[ \phi^s + \phi^f = 1 \]

Note that \(\phi^f\) is the commonly used concept porosity \(\phi\) in soild mechanics. If it’s not mentioned, 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.

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}_f = \bm v_f - \bm v_s = \bm v_f - \bm v \]

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

(31)\[ \dot{\bm w} = \phi^f \tilde{\bm v}_f = \phi^f \left( \bm v_f - \bm v \right) \]

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

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 mixture theory is given by

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

where

(35)\[ \rho^b = (1 - \phi^f) \rho^s + \phi^f \rho^f, \]

is the bulk density of the porous medium. 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_f + \rho^f \ddot {\bm u} + \rho^w \ddot{\bm w} - \bm f^f \right), \]

where \(\bm f^f\) is the fluid body force,

\[ \varkappa = \frac{\varkappa_{\circ}}{\eta_f}, \]

is the permeability coefficient, with \(\varkappa_{\circ}\) the intrinsic permeability, and \(\eta_f\) the dynamic viscosity of the fluid phase, and

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

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

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

where \(\alpha_{\infty}\) is the tortuosity.

Poroelasticity at finite strain

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

(37)\[ \begin{aligned} m_{\alpha} &= \int_{\Omega^{\alpha}} dm_{\alpha} = \int_{\Omega^{\alpha}} \rho^{\alpha R}dv_{\alpha} = \int_{\Omega} \rho^{\alpha R}\phi^{\alpha} dv = \int_{\Omega} \rho^{\alpha} dv = \int_{\Omega^{\alpha}_0} \rho^{\alpha} J_{\alpha} dV_{\alpha}, \\ \rho^{\alpha R}(\bm x, t) &= \frac{dm_{\alpha}}{dv_{\alpha}}, \\ \rho^{\alpha}(\bm x, t) &= \frac{dm_{\alpha}}{dv} = \rho^{\alpha R} \phi^{\alpha}(\bm x, t), \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^{s} \rho^{\alpha}}{D t} + \tilde{\bm v}_{\alpha} \cdot \nabla_x \bm{\rho}^{\alpha}, \]

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

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

Using barotropic flow relation (for constant phase temperature \(\theta^{\alpha}\))

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

replacing \( D^f/D t\) by \(D^s/D t\) and \(\bm{v}_f = \tilde{\bm v}_f + \bm{v}_s\), then add both phases we get

(42)\[ \frac{\phi^s}{\bulk_s}\frac{D^{s} p_s}{D t} + \frac{\phi^f}{\bulk_f}\frac{D^{s} p_f}{D t} + \nabla_x\cdot \bm{v}_{s} + \nabla_x\cdot (\phi^f \tilde{\bm v}_f) + \frac{\phi^f \tilde{\bm v}_f}{\bulk_f}\cdot \nabla_x p_f = \frac{\hat{\rho}^{s}}{\rho^{sR}} + \frac{\hat{\rho}^{f}}{\rho^{fR}} \]

Assume the solid phase is nearly incompressible and there is no supply of solid mass and use the definition (31), we get simplified balance of mass for solid-fluid mixture as

(43)\[ \frac{\phi^f}{\bulk_f}\frac{D^{s} p_f}{D t} + \nabla_x\cdot \bm{v}_{s} + \nabla_x\cdot (\dot{\bm w}) + \frac{\dot{\bm w}}{\bulk_f}\cdot \nabla_x p_f = \frac{\hat{\rho}^{f}}{\rho^{fR}}. \]

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

(44)\[ \frac{J_s \phi^f}{\bulk_f}\frac{D^{s} p_f}{D t} + \frac{D^{s} J_s}{D t} + J_s \nabla_{X_s}(\dot{\bm w})\tcolon \bm{F}_s^{-T} + \frac{J_s \dot{\bm w}}{\bulk_f} \cdot \bm{F}_s^{-T}\nabla_{X_s} p_f = \frac{J_s\hat{\rho}^{f}}{\rho^{fR}}. \]

It should be noted that from (41) as \(\bulk \to \infty \), we get \(D^s \rho^{sR}/D t =0\), and from (38) and (39) we can derive

(45)\[ \phi^s = \phi^s_0/ J_s, \]

where \(\phi^s_0\) is the initial value of \(\phi^s\) at \(t = 0\). Furthermore, in the Darcy’s law (36), the permeability coefficient can be a function of porosity as

\[ \varkappa = \frac{\varkappa_{\circ}}{\eta_f} \frac{\mathcal{F}(\phi^f)}{\mathcal{F}(\phi^f_0)}, \]

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

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