Thermodynamically consistent variational theory of porous media with a breaking component

Abstract

If a porous media is being damaged by excessive stress, the elastic matrix at every infinitesimal volume separates into a ‘solid’ and a ‘broken’ component. The ‘solid’ part is the one that is capable of transferring stress, whereas the ‘broken’ part is advecting passively and is not able to transfer the stress. In previous works, damage mechanics was addressed by introducing the damage parameter affecting the elastic properties of the material. In this work, we take a more microscopic point of view, by considering the transition from the ‘solid’ part, which can transfer mechanical stress, to the ‘broken’ part, which consists of microscopic solid particles and does not transfer mechanical stress. Based on this approach, we develop a thermodynamically consistent dynamical theory for porous media including the transfer between the ‘broken’ and ‘solid’ components, by using a variational principle recently proposed in thermodynamics. This setting allows us to derive an explicit formula for the breaking rate, i.e., the transition from the ‘solid’ to the ‘broken’ phase, dependent on the Gibbs’ free energy of each phase. Using that expression, we derive a reduced variational model for material breaking under one-dimensional deformations. We show that the material is destroyed in finite time, and that the number of ‘solid’ strands vanishing at the singularity follows a power law. We also discuss connections with existing experiments on material breaking and extensions to multi-phase porous media.

Appendices

Variational formulation for a viscous medium with heat conduction

In this Appendix, we review the variational theory of a single-component fluid or elastic medium with heat conduction and viscosity, based on the variational thermodynamics approach developed in [33,34,35]. This work was generalized for porous media in [32], in the case where the elastic medium consists of a single component. We use this variational thermodynamics framework in Sect. 3 to treat the case of a two-component (‘solid’ and ‘broken’) elastic media which includes the additional irreversible processes of transitions between the components, where both components are embedded in the fluid.

1.1 Lagrangian (material) description

In the material description, the variational formulation of thermodynamics is an extension of the Hamilton principle of continuum mechanics. We recall the material description here for pedagogical reasons since the variational formulation takes its simplest form. The variational formulation in the spatial (Eulerian) description is then derived from it, see the next section, and turns out to be a useful tool for porous media modeling.

For a single-component fluid or elastic medium, we take the Lagrangian L to be a function of the configuration diffeomorphism $\varvec{\varphi }(t, {\textbf{X}} )$, see Sect. 2.1.2, the material velocity ${\varvec{V}}(t, {\textbf{X}} )= \dot{\varvec{\varphi }}(t, {\textbf{X}} )$, and the entropy density $S(t, {\varvec{X}})$. It also depends parametrically on the density $\varrho _0({\varvec{X}})$ and on the Riemannian metric $G_0({\varvec{X}})$ on the reference configuration which can be chosen as the Euclidean one when ${\mathcal {B}} \subset {\mathbb {R}} ^3$. The metric $G_0({\varvec{X}})$ is needed to introduce the Eulerian deformation tensors.

For each fixed $\varrho _0( {\textbf{X}} )$ and $G_0( {\textbf{X}} )$, we can write the Lagrangian as a function $L: T {\text {Diff}}_0( {\mathcal {B}}) \times {\text {Den}}( {\mathcal {B}} ) \rightarrow {\mathbb {R}}$, where ${\text {Diff}}_0( {\mathcal {B}} )\ni \varvec{\varphi }$ is the group of diffeomorphisms of ${\mathcal {B}}$ keeping $\partial {\mathcal {B}}$ pointwise fixed, corresponding to no-slip boundary conditions, $T{\text {Diff}}_0( {\mathcal {B}} )\ni (\varvec{\varphi }, \dot{\varvec{\varphi }})$ is its tangent bundle, and ${\text {Den}}( {\mathcal {B}} )\ni S$ is the space of densities on ${\mathcal {B}}$. In the material description, a standard expression is given by

$$\begin{aligned} L( \varvec{\varphi }, \dot{ \varvec{\varphi } }, S)= \int _ {\mathcal {B}} \left[ \frac{1}{2} \varrho _0 | \dot{\varvec{\varphi } }|- \varrho _0 \,{\mathcal {E}} ( \nabla \varvec{\varphi }, \varrho _0, S, G_0) \right] \textrm{d}^3 {\varvec{X}}, \end{aligned}$$

(97)

with ${\mathcal {E}}$ the internal energy in the material description, subject to the usual covariance assumptions, [36, 49]. In particular, due to the material covariance assumption, the internal energy can be written in terms of the Eulerian variables $\rho , s, b$ see Sect. A.2, as

$$\begin{aligned} {\mathcal {E}} ( \nabla \varvec{\varphi }, \varrho _0, S, G_0)= e( \rho , s, b) \circ \varvec{\varphi }, \end{aligned}$$

(98)

with e the specific internal energy in the Eulerian description.

The critical action principle for thermodynamics needs the introduction of two additional variables: $\Sigma (t,{\varvec{X}})$ which is identified with the entropy generated by the irreversible processes (unlike $S(t, {\textbf{X}} )$, which is the total entropy), and $\Gamma (t,{\varvec{X}})$ identified with the thermal displacement [34]. With that notation, the critical action principle for a heat conducting viscous continuum reads [33, 34]:

$$\begin{aligned} \delta \int _0^ T \left[ L (\varvec{\varphi }, \dot{ \varvec{\varphi } }, S) + \int _ {\mathcal {B}} (S- \Sigma ) {\dot{\Gamma }} \,\textrm{d}^3 {\varvec{X}} \right] \textrm{d}t=0\,, \end{aligned}$$

(99)

subject to the phenomenological constraint

$$\begin{aligned} \frac{\delta L}{ \delta S}{\dot{\Sigma }} = - {\varvec{P}}: \nabla \dot{\varvec{\varphi }} + {\varvec{J}}_S \cdot \nabla \dot{ \Gamma } \end{aligned}$$

(100)

and with respect to variations $\delta \varvec{\varphi }$, $\delta S$, $\delta \Sigma$, $\delta \Gamma$ subject to the variational constraint

$$\begin{aligned} \frac{\delta L}{ \delta S}\delta \Sigma = - {\varvec{P}}: \nabla \delta \varvec{\varphi } + {\varvec{J}}_S \cdot \nabla \delta \Gamma . \end{aligned}$$

(101)

The tensor ${\varvec{P}}$ is the Piola–Kirchhoff viscous stress tensor and ${\varvec{J}}_S$ is the entropy flux density in Lagrangian representation, see [34]. They are the Lagrangian objects corresponding to the more widely used Eulerian viscous stress tensor $\varvec{\sigma }$ and Eulerian entropy flux density ${\varvec{j}}_s$, see below.

An application of (99)–(101) yields, in addition to the equations of motion in the Lagrangian frame that we do not present here, the conditions

$$\begin{aligned} {\dot{\Sigma }} = \dot{S} + {\text {DIV}} {\varvec{J}}_S \;\; \text{ and } \;\;{\dot{\Gamma }} = - \frac{\delta L}{\delta S} = T \end{aligned}$$

(102)

with T being the temperature, see [34], with the notation ${\text {DIV}}$ indicating the divergence operator in the reference frame, as opposed to the divergence operator ${\text {div}}$ in the spatial frame. These conditions attribute to $\Sigma$ and $\Gamma$ their physical meaning mentioned above.

Since it is $\Sigma$ that describes the entropy of irreversible processes, from the second law of thermodynamics, we must have

$$\begin{aligned} {\dot{\Sigma }} \ge 0 \,, \end{aligned}$$

(103)

whereas $\dot{S}$ does not necessarily has to have a particular sign.

Remark A.1

(Structure of the variational formulation) Observe that (99)–(101) is an extension of the Hamilton principle $\delta \int _0^T L( \varvec{\varphi }, \dot{\varvec{\varphi }} ) \textrm{d} t=0$ of continuum mechanics which involves two types of constraints: the constraint (100) on the critical curve and the constraint (101) on the variations to be used when computing this critical curve. One passes from (100) to (101) by formally replacing the time rate of changes by $\delta$-variations for each irreversible processes, such as $F_i \dot{x}^i \leadsto F_i \delta x^i$ in the case of an irreversible process due to a friction force for a finite dimensional system, see [33,34,35]. This variational formulation is reminiscent from the Lagrange-d’Alembert principle used in nonholonomic mechanics. A finite dimensional version of the variational formulation (99)–(101) will be used as a modeling tool in Sect. 5.

1.2 Eulerian (spatial) description

The Eulerian variables are the velocity ${\varvec{u}}(t, {\varvec{x}})$, mass density $\rho (t, {\varvec{x}})$, entropy density $s(t, {\varvec{x}})$, and Finger deformation tensor $b(t, {\varvec{x}})$ defined from the material variables $\varvec{\varphi } (t, {\varvec{X}})$, $\dot{ \varvec{\varphi }} (t, {\varvec{X}})$, $\varrho _0 (t, {\varvec{X}})$, $S (t, {\varvec{X}})$, and $G_0 (t, {\varvec{X}})$ in the usual way. In particular, the Finger deformation tensor is the push-forward of the inverted metric by the configuration diffeomorphism $\varvec{\varphi }$:

$$\begin{aligned} b= \varvec{\varphi } _* G_0 ^{-1}, \end{aligned}$$

(104)

see [36, 49]. In the case when ${\mathcal {B}}$ is a subset of ${\mathbb {R}}^3$ and $G_0$ is the identity tensor, the Finger deformation tensor b is written in terms of the deformation gradient tensor ${\mathbb {F}}= \nabla \varvec{\varphi }$ as

$$\begin{aligned} b = {\mathbb {F}} \cdot {\mathbb {F}}^T \,, \quad b^{ij}(t, {\varvec{x}}) = \frac{\partial x^i}{\partial X^k}\frac{\partial x^j}{\partial X^k}( \varvec{\varphi }^{-1}(t, {\varvec{x}})). \end{aligned}$$

(105)

Note that from $\varvec{\varphi } \in {\text {Diff}}_0( {\mathcal {B}})$, we have the no-slip boundary conditions ${\varvec{u}}=0$ on $\partial {\mathcal {B}}$. We also consider the thermal displacement $\gamma (t, {\varvec{x}})$, internal entropy density $\sigma (t, {\varvec{x}})$, Cauchy stress $\varvec{\sigma } (t, {\varvec{x}})$ and entropy flux ${\varvec{j}}_s (t, {\varvec{x}})$ defined as the Eulerian quantities associated to $\Gamma (t, {\varvec{X}})$, $\Sigma (t, {\varvec{X}})$, ${\varvec{P}} (t, {\varvec{X}})$, and ${\varvec{J}}_S (t, {\varvec{X}})$, see [34].

With these definitions, the spatial version of (99)–(101) gives the variational formulation

$$\begin{aligned} \delta \int _0^ T \left[ \ell ({\varvec{u}}, \rho , s, b) + \int _ {\mathcal {B}} (s - \sigma ) D_t \gamma \,\textrm{d}^3 {\varvec{x}} \right] \textrm{d}t=0\,, \end{aligned}$$

(106)

subject to the phenomenological constraint

$$\begin{aligned} \frac{\delta \ell }{ \delta s} {\bar{D}}_t \sigma = - \varvec{\sigma }: \nabla {\varvec{u}} + {\varvec{j}}_s \cdot \nabla D_t \gamma \end{aligned}$$

(107)

and with respect to variations

$$\begin{aligned} \delta {\varvec{u}}= \partial _t \varvec{\eta } + {\varvec{u}} \cdot \nabla \varvec{\eta } - \varvec{\eta } \cdot \nabla {\varvec{u}}, \quad \delta \rho = - {\text {div}}( \rho \varvec{\eta } ), \quad \delta b= - \pounds _ { \varvec{\eta } }b, \end{aligned}$$

(108)

$\delta s$, $\delta \sigma$, and $\delta \gamma$ subject to the variational constraint

$$\begin{aligned} \frac{\delta \ell }{ \delta s}{\bar{D}}_ \delta \sigma = - \varvec{ \sigma }: \nabla \varvec{\eta } + {\varvec{j}}_s \cdot \nabla D_ \delta \gamma . \end{aligned}$$

(109)

We recall that $\varvec{\eta }(t, {\varvec{x}})= \delta \varvec{\varphi } (t, \varvec{\varphi } ^{-1} (t, {\varvec{x}}))$ denotes the variation of the fluid trajectories in the Eulerian frame. In the expression of $\delta b$, $\pounds _{ \varvec{\eta }} b$ denotes the Lie derivative of the symmetric contravariant tensor b in the direction $\varvec{\eta }$, which follows from (104). We have introduced the notations

$$\begin{aligned} \begin{aligned}&D_t f= \partial _t f + {\varvec{u}} \cdot \nabla f&\qquad&D_ \delta f= \delta f + \varvec{\eta } \cdot \nabla f \\&{\bar{D}}_t f = \partial _t f + {\text {div}}(f {\varvec{u}})&\qquad&{\bar{D}}_ \delta f = \delta f + {\text {div}}(f \varvec{\eta }) \end{aligned} \end{aligned}$$

(110)

for the Lagrangian time derivative and Lagrangian variations of a scalar function and a density.

A direct application of the variational principle (106)–(109) yields the general equations of motion for a heat conducting viscous continuum with Lagrangian $\ell ( {\varvec{u}}, \rho , s, b)$ in Eulerian coordinates as

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \partial _t \frac{\delta \ell }{\delta {\varvec{u}} }+\pounds _{{\varvec{u}}}\frac{\delta \ell }{\delta {\varvec{u}}} =\rho \nabla \frac{\delta \ell }{\delta \rho } +s\nabla \frac{\delta \ell }{\delta s}- \frac{\delta \ell }{\delta b}:\nabla b- 2{\text {div}} \left( \frac{\delta \ell }{\delta b}\cdot b \right) + {\text {div}} \varvec{\sigma } \\ \displaystyle \frac{\delta \ell }{\delta s}( {\bar{D}}_t s + {\text {div}} {\varvec{j}}_s) = - \varvec{\sigma }: \nabla {\varvec{u}}- {\varvec{j}}_s \cdot \nabla \frac{\delta \ell }{\delta s} \end{array} \right. \end{aligned}$$

(111)

together with the conditions

$$\begin{aligned} {\bar{D}}_t \sigma = {\bar{D}}_t s+ {\text {div}} {\varvec{j}}_s \qquad \text {and}\qquad {\bar{D}}_t \gamma = - \frac{\delta \ell }{\delta s}. \end{aligned}$$

(112)

From these two conditions, ${\bar{D}}_t \sigma$ is interpreted as the total entropy generation rate density and $D_t \gamma$ is the temperature, hence $\gamma$ is the thermal displacement. The equations for $\rho$ and b are

$$\begin{aligned} \partial _t \rho + {\text {div}}( \rho {\varvec{u}} )=0 \quad \text {and}\quad \partial _t b + \pounds _{ {\varvec{u}}} b=0, \end{aligned}$$

which follow from the definition of $\rho$ and b in terms of $\varrho _0$ and $G_0$. Also, the variations $\delta \gamma$ at the boundary yields the insulated boundary conditions ${\varvec{j}} _s \cdot {\varvec{n}}=0$ on $\partial {\mathcal {B}}$. In the fluid momentum equations above, $\pounds _ {{\varvec{u}}} \frac{\delta \ell }{\delta {\varvec{u}}}$ denotes the Lie derivative of the fluid momentum density $\frac{\delta \ell }{\delta {\varvec{u}}}$ in the direction ${\varvec{u}}$, explicitly given as $\pounds _ {{\varvec{u}}} {\varvec{m}} = {\varvec{u}} \cdot \nabla {\varvec{m}} + \nabla {\varvec{u}}^{\textsf{T}} {\varvec{m}} + {\varvec{m}} {\text {div}} {\varvec{u}}$.

By using the standard expression of the Lagrangian

$$\begin{aligned} \ell ({\varvec{u}}, \rho , s, b) = \int _{{{\mathcal {B}}}} \left[ \frac{1}{2} \rho |{\varvec{u}}|^2 - \rho e(\rho ,s/ \rho ,b)\right] \textrm{d}^3 {\varvec{x}}\,, \end{aligned}$$

(113)

which is the Eulerian form of (97), see [36], the equations of motion (114) give the following system of equations for a viscous and heat conducting continuum

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \rho (\partial _t {\varvec{u}} + {\varvec{u}} \cdot \nabla {\varvec{u}}) = - \nabla p + {\text {div}} \varvec{\sigma }_{\textrm{el}}+ {\text {div}} \varvec{\sigma }\\ \displaystyle T({\bar{D}}_t s + {\text {div}} {\varvec{j}}_s) = \varvec{\sigma }: \nabla {\varvec{u}}- {\varvec{j}}_s \cdot \nabla T, \end{array} \right. \end{aligned}$$

(114)

with $p= \rho ^2 \frac{\partial e}{\partial \rho }$ the pressure, $T= \frac{\partial e}{\partial \eta }$ the temperature, and $\varvec{\sigma } _{\textrm{el}}=2 \rho \frac{\partial e}{\partial b}\cdot b$ the elastic stress. We refer to [31, 33,34,35] for the statement of the variational formulation in both the material and spatial description as well as the detailed computations and several applications and extensions.

Remark A.2

(Structure of the variational formulation) The variational formulation (106)–(109) inherits from its Lagrangian counterpart (99)–(101) the same structure, in which the variational constraint (109) follows from the phenomenological constraint (107) by formally replacing the time rate of changes by $\delta$-variation, now in the Eulerian setting, see [33,34,35]. The same structure underlies the variational approach to porous media developed in Sect. 3.1.

Alternative expressions for heat and mass transfer coefficients involving cross-phenomena

We describe here possible cross-effects in the reversible processes, based on the form of the entropy production equation (30).

Regarding mass and heat transfer, cross-effects can occur for $k=s,\ell =b$, while still preserving (38), which are the discrete analogues to the Soret and Dufour effects. We can thus consider

$$\begin{aligned} \begin{bmatrix} J_{sb} \frac{T_s-T_b}{T_s T_b} \\ q_{sb} \end{bmatrix} = {\mathcal {L}} _{sb} \begin{bmatrix} T_b - T_s \\ \frac{1}{T_s}\big ( \frac{1}{2} | {\varvec{u}}_s| ^2 - g_s \big )- \frac{1}{T_b}\big ( \frac{1}{2} | {\varvec{u}}_b| ^2 - g_b\big ) \end{bmatrix} \end{aligned}$$

(115)

with the symmetric part of the $2 \times 2$ matrix ${\mathcal {L}} _{sb}$ positive, and with $J_{sf}$ and $J_{bf}$ satisfying (39) or, equivalently (41). We shall note that the conditions (42) and (115) for the mass and heat transfer coefficients seem to be novel and not encountered in previous literature on the subject. On the other hand, since $q_{sb}$ can be interpreted physically as the rate of damage to the elastic matrix, these equations do have some aspects of the phenomenological equations for the evolution of damage parameters used previously in the literature, see, e.g., [48]. The equations for transfer rate q described there, in our notation, would be an affine function of the internal pressure p. Our expression has these pressure terms, but also involves additional terms related to the elastic energy, as we show below.

Besides the cross-phenomena between the scalar processes associated to mass and heat transfer, another cross-effect can be postulated, mathematically similar to the cross-phenomena between bulk viscosity and chemistry, see [22]. Ignoring the friction forces and assuming that the stresses $\varvec{\sigma } _{k\ell }$ are symmetric, we can rewrite the entropy production expression (30) as

$$\begin{aligned}&\sum _{k<\ell }\Bigg [ \varvec{\sigma } _{k\ell }^\circ : \left( \frac{{\mathbb {F}} _k^\circ }{T_k}- \frac{{\mathbb {F}} _\ell ^\circ }{T_\ell } \right) + \frac{1}{3} {\text {Tr}}( \varvec{\sigma } _{k\ell }) \left( \frac{{\text {div}} {\textbf{u}} _k}{T_k} - \frac{{\text {div}} {\textbf{u}} _\ell }{T_\ell } \right) \\&\qquad + q_{k\ell }\left( \frac{\frac{1}{2} | {\varvec{u}}_k| ^2 - g_k}{T_k}-\frac{\frac{1}{2} | {\varvec{u}}_\ell | ^2 - g_\ell }{T_\ell } \right) + J_{k\ell }\left( \frac{1}{T_\ell } - \frac{1}{T_k}\right) (T_\ell -T_k)\Bigg ]\,. \end{aligned}$$

Therefore, one is naturally led to postulate the following cross-effects of scalar processes:

$$\begin{aligned} \left[ \begin{array}{c} \frac{1}{3} {\text {Tr}}( \varvec{\sigma } _{k\ell })\\ \displaystyle q_{k\ell }\\ J_{k\ell } \frac{T_k-T_\ell }{T_kT_\ell } \end{array} \right] = {\mathcal {L}}_{k\ell } \left[ \begin{array}{c} \frac{1}{T_k}{\text {div}} {\textbf{u}} _k - \frac{1}{T_\ell }{\text {div}} {\textbf{u}} _\ell \\ \frac{1}{T_k} (\frac{1}{2} | {\varvec{u}}_k| ^2 - g_k)-\frac{1}{T_\ell }(\frac{1}{2} | {\varvec{u}}_\ell | ^2 - g_\ell ) \\ \displaystyle T_\ell -T_k \end{array} \right] , \end{aligned}$$

(116)

where the symmetric part of the $3 \times 3$ matrices ${\mathcal {L}}_{k\ell }$ must be positive from the second law. We shall not endeavor to consider the more general conditions (115) and (116) involving the cross-effects, and only concentrate on the simpler conditions (40) as the one leading to the most simple mathematical expressions treatable analytically. In spite of the relative mathematical simplicity compared to (115) and (116), expressions (40) lead to physically relevant systems providing physically valid quantitative predictions.