Thermodynamics in Mono and Biphasic Continuum Mechanics

1. Introduction

This chapter applies the laws of thermodynamics to problems in continuum mechanics. Initially these are applied to a monophasic medium. The case of a biphasic porous medium is then treated with the aim of illustrating how a framework may be established for capturing possible couplings in the pertinent constitutive relationships. This approach is founded on the two laws of Thermodynamics. The first law expresses the conservation of energy when considering all possible forms while the second law postulates that the quality of energy must inevitably deteriorate in relation to its transformability into efficient mechanical work.

2. The principles of thermodynamics in the case of monophasic media

In order to simplify matters so that the reader can have a good intuitive understanding on the fundamental principles, in particular their physical contents, we begin with the simplest case of a monophasic continuous media.

Consider a solid body in movement, with mass density $\rho$ and a velocity field $v$ (figure 1). Our attention will be focused on an arbitrarily chosen part of this body, which occupies a volume ${\text{Ω}}_{\text{t}}$ at time $t$ . For ordinary problems of solid mechanics, we are concerned with mechanical and thermal energies. We therefore suppose that the body inside ${\text{Ω}}_{\text{t}}$ is subject to a distributed body force $f$ (for example gravity) and surface tractions $t$ on its boundary surface, noted $\partial {\text{Ω}}_{\text{t}}$ . At the same time, the body is subject to a heat flux $q$ on $\partial {\text{Ω}}_{\text{t}}$ and an internal heat source $r_q$ .

To begin with, we consider the energy and entropy balance of all the matter inside the volume ${\text{Ω}}_{\text{t}}$ , using the two principles of thermodynamics.

3. The first principle of thermodynamics

The first principle stipulates that energy must be conserved under its different forms. Limiting our study here to thermal and mechanical energies, we can write:

Figure 1.

Formula: Eqn013.wmf>is a generic part of a body in movement, with distributed body forces ${\text{Ω}}_{\text{t}}$ , surface tractions $f$ inward heat flux $t$ , and distributed heat source $q$ , $r_q$ is the outward unit normal.

$$n$$ E1

In the above equation, $\frac{d}{dt}\left(E+K\right)=P_x+Q$ and $E$ are the global internal and kinetic energies, while $K$ and $P_x$ are the total external supply of mechanical and thermal power for all matters inside $Q$ . The time derivative refers to the rate of increase of the energy content by following the same ensemble of material particles in their movement. This equation simply states that heat and mechanical energies received by a body which are not converted into kinetic energy become the internal energy. In continuum mechanics, physical quantities vary spatially from one point to another. The global quantities can be expressed in terms of the sum of local quantities:

$${\text{Ω}}_{\text{t}}$$ E2

where $E=\underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}\rho ed{\text{Ω}}_{\text{t}}K=\underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}\frac{1}{2}\rho v\Delta vd{\text{Ω}}_{\text{t}}{P}_x=\underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}f\Delta vd{\text{Ω}}_{\text{t}}+\underset{\partial }{\overset{\text{Ω}}{{\displaystyle \int }}}t\Delta vd\text{S}Q=\underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}{r}_{q}d{\text{Ω}}_{\text{t}}-\underset{\partial }{\overset{\text{Ω}}{{\displaystyle \int }}}q\Delta nd\text{S}$ , the specific internal energy is defined as the internal energy per unit of mass. Substitution of equation (2) into (1) and on account of the classic equation $e$ relating the surface traction $t=\sigma \Delta n$ to the second order symmetric stress tensor $t$ , we get after some simplifications:

$$\sigma$$ E3

where $i$ denotes the strain tensor and a dot above a variable denotes the material derivative (i.e. total derivative with respect to time) by following the movement of an elementary solid particle. Internal energy is the energy content within a given mass of material. This includes the (A) kinetic energy due to the disordered thermal agitation and the (B) interaction, or potential, energy between molecules due to their relative positions (for example the elastic strain energy). It is the macroscopic description of (A) that leads to the introduction of the absolute temperature. The internal energy can also be the energy stored due to concentration of solutes (osmotic potential), but is outside the scope of this presentation. However, it should be noted that the following energies are not counted as internal energy:

1. Kinetic energy due to the macroscopic (ordered) movement of a material body

2. Potential energy due to the position of a body relative to an external field such as gravity

The last form of energy, namely the macroscopic potential energy, is accounted for by considering conservative body forces derivable from a potential, such as the gravity force per unit volume $\rho \dot{e}=\sigma :\dot{\epsilon }+r_q-divq$ , in the term $\rho g$ in the definition of $f$ . Note that relative to the first principle, all forms of energy have an equal status.

4. The second principle of thermodynamics

The second principle of Thermodynamics confers a special status to heat, and distinguishes it from all other forms of energy, in that:

1. Once a particular form of energy is transformed into heat, it is impossible to back transform the entire amount to its original form without compensation.

2. To convert an amount of heat energy $P_x$ into useful work, a necessary condition is to have at least two reservoirs with two different (absolute) temperatures $\text{Δ}Q$ and $T_1$ (suppose $T_2$ to fix ideas).

3. Moreover, the above conversion can at best be partial in that the amount of work $\mathrm{e}W$ extractable from a given quantity of heat $T_1{T}_2$ admits a theoretical upper bound depending on the two temperatures:

$$\text{Δ}Q$$ E4

Figure 2.

The heat engine represented by the circle takes a quantity of heat $\frac{\text{Δ}W}{\text{Δ}Q}\le \frac{T_1-T_2}{T_2}$ from the hotter reservoir $\text{ΔQ}$ and rejects $T_1$ to the colder reservoir $\text{ΔQ'}$ , while it performs an amount of useful work $T_2$ . The first principle requires $\text{Δ}W$ and the second principle sets a theoretical upper bound on the efficiency $\text{Δ}W=\text{ΔQ}-\text{ΔQ'}$ attainable by heat engines.

Note that real efficienciesobtainable in practical cases are far less than that suggested by equation (4) due to unavoidable frictional losses. In the limit when the temperature becomes uniform, no mechanical work can be extracted anymore and this corresponds to some kind of thermal-death.In technical terms, when a particular form of energy is transformed into heat, the energy is degraded and becomes less available to perform useful work. The second principle gives a systematic and consistent account of why heat engines have theoretical upper limits of efficiency, and why certain phenomena can never occur spontaneously. Forexample, we cannot extract sea water at 20 C, cool it down to 0 C by extracting heat from it, and use that heat to drive the turbine and advance a ship!The theoretical formulation of the second principle viathe concept of entropy derives its basis from a very large quantity of observations. The counter-part of the generality of its validity is the high level of abstraction, making it difficult to understand.Classical irreversible thermodynamics formulated directly at the macroscopic scale has an axiomatic appearance. The entropy change is defined axiomatically with respect to heat exchange and production. To understand its molecular original requires investigations at the microscopic scale. This is not necessary if the objective is to apply thermodynamic principles to build phenomenological models, although such investigations do contribute to a better understanding of the physical origin of the phenomena.Clausius (1850) invented the thermodynamic potential - the entropy - to describe this uni-directional and irreversible degradation of energy. Formulated in terms of entropy, the second principle of thermodynamics says that whenever some form of energy is transformed into heat, the global entropy increases. It can at best stay constant for reversible processes but can never decrease. If we denote $s$ the specific entropy (per unit mass), the second principle writes:

$$\text{Δ}W/\text{ΔQ}$$ E5

In other words, for a fixed quantity of matter, the entropy increase must be greater than (resp. equal to) external heat supply divided by the absolute temperature in irreversible (resp. reversible) processes. The difference is due to other forms of energy being transformed into heat via dissipative processes. In our study here, this corresponds to internal frictional processes transforming mechanical energy into heat. Once this occurs, the process becomes irreversible. The previous inequality can be simplified to the following local form using Gauss' theorem:

$$\frac{d}{dt}\underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}\rho sd{\text{Ω}}_{\text{t}}\ge \underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}\frac{r_q}{T}d{\text{Ω}}_{\text{t}}-\underset{\partial }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}\frac{q\Delta n}{T}d\text{S}$$ E6

As a macroscopic theory, irreversible thermodynamics does not give any explanation on the origin of entropy. Similarly to the case of plastic strains, the manipulation of entropy and other thermodynamic potentials will rely on postulated functions, valid over finite domains and containing coefficients to be determined by experiments.

5. Clausius-Duheim (CD) inequality

Combining the first and the second principle, we obtain the classic Clausius-Duhem (CD) inequality in the context of solid mechanics (electric, magnetic, chemical or osmotic terms etc. can appear in more general problems):

$$\rho \dot{s}+div\frac{q}{T}-\frac{r_q}{T}\ge 0$$ E7

In the limiting case when the temperature field is uniform and the process is reversible, the above inequality becomes equality:

$$\text{Φ}=\sigma :\dot{\epsilon }+\rho \left(T\dot{s}-\dot{e}\right)-\frac{1}{T}q\Delta gradT\ge 0$$ E8

Since the specific internal energy is a state function and is supposed to be entirely determined by the state variables, we conclude from the differential form in (8) that $e$ depends naturally on $\frac{1}{\rho }\sigma :\dot{\epsilon }+T\dot{s}-\dot{e}=0\text{or d}e=\frac{1}{\rho }\sigma :\text{d}\epsilon +Td\text{s}$ and $\epsilon$ (i.e. $s$ ) and that the following state equations hold:

$$e=e\left(\epsilon ,s\right)$$ E9

However, the specific entropy $s$ is not a convenient independent variable as it is intuitively difficult to comprehend and practically difficult to control. The classical approach consists of introducing another state function, the specific Helmholtz's free energy, via the Legendre transform:

$$\sigma =\rho \frac{\partial e}{\partial \epsilon }T=\frac{\partial e}{\partial s}$$ E10

to recast inequality (7) to the following form:

$$\psi =e-Ts$$ E11

Again, in the absence of dissipative phenomena and a uniform temperature field, we have:

$$\text{Φ}=\sigma :\dot{\epsilon }-\rho \left(\dot{\psi }+s\dot{T}\right)-\frac{1}{T}q\Delta gradT\ge 0$$ E12

via the same reasoning as previously, we deduce that the specific free energy $\frac{1}{\rho }\sigma :\dot{\epsilon }-s\dot{T}-\dot{\psi }=0\text{or d}\psi =\frac{1}{\rho }\sigma :\text{d}\epsilon -sd\text{T}$ depends naturally on $\psi$ and $\epsilon$ and satisfies the following state equations:

$$T$$ E13

The Legendre transform (10) thus allows one to define a thermodynamic potential with natural independent variables which are more accessible ( $\sigma =\rho \frac{\partial \psi }{\partial \epsilon }s=-\frac{\partial \psi }{\partial T}$ instead of $T$ in the present case). The quantity $s$ , having the unit of energy per unit volume per unit time, is called total dissipation. It represents the transformation of non-thermal energy into heat via frictional processes, which then becomes less available.

6. How to use the second principle

There are two ways to make use of the second thermodynamic principle. We can first of all verify the consistency or the inconsistency of a given model with respect to the 2^nd principle, in an a posteriori manner, in the sense that the construction of the model does not rely in any way on the 2^nd principle. On the other hand, we can actually construct a model, starting from the Clausius-Duhem inequality, by specifying appropriate functional forms for the Helmholtz's free energy and the dissipation. Naturally, there is no unique way to achieve this goal since thermodynamics does not supply any information on the specific behavior of a particular material under study. This process must therefore integrate experimental data so that the model predictions are consistent with the reality. Among different representations (or models) consistent with thermodynamic principles, the best is the one with a clear logical structure and comprising a minimum number of parameters (simplicity). This last criterion allows to minimise the amount of experimental work necessary to identify these parameters, which is always a very time-consuming task.

7. Implicit but essential assumptions

All classic developments based on irreversible thermodynamics assume implicitly that the process does not deviate significantly from thermodynamic equilibrium. In consequence, despite the fact the system is in evolution therefore in non-equilibrium, the state equation expressing the condition of thermodynamic equilibrium can still be used to reduce the number of independent state parameters by one in complex problems (for example, the density, pressure and temperature of the pore fluid transiting a porous solid is related by a state equation). This is strictly speaking an approximation. Its efficiency can only be assessed a posteriori by the results.

In a heterogeneous system, the thermodynamic state hence the state parameters are position-dependent. This heterogeneity (hence non-equilibrium) is the driving force which tends to restore the system back to thermodynamic equilibrium. However, it is assumed that the (spatial) variation is sufficiently mild so that every elementary particle can be considered as under thermodynamic equilibrium. Its state parameters aretherefore linked by the stateequation expressing this equilibrium requirement. This assumption is called the “hypothesis of local equilibrium”. This assumption excludes the treatment of fast processes (for example explosions) under the framework of classic irreversible thermodynamics.

8. Applications to plasticity and viscoplasticity: General equations

To illustrate how thermodynamic principles can be used to formulate physical laws, let us consider the particular case of the inelastic behaviour of solids. The classic partition:

$$\text{Φ}$$ E14

Is assumed, where $\epsilon ={\epsilon }^e+{\epsilon }^p$ is the elastic strain and ${\epsilon }^e$ denotes for the time being all forms of irreversible (i.e. inelastic) strains. In order to satisfy the CD inequality (11), a common practice is to assume that ${\epsilon }^p$ , so that $\psi =\psi ({\epsilon }^e,T,V_k)$ . The scalar variables grouped into a tensor $\dot{\psi }=\frac{\partial \psi }{\partial \epsilon }\dot{{\epsilon }^e}+\frac{\partial \psi }{\partial T}\dot{T}+\frac{\partial \psi }{\partial V_k}{\dot{V}}_{\text{k}}$ are internal variables introduced to account for the state-dependent non-linear inelastic behaviour. In practice, this is often the irreversible strains or their scalar invariants. The CD inequality then becomes:

$$V_k$$ E15

Consider the particular case of elastic (reversible) evolution corresponding to stationary values of the internal variables $\text{Φ}=\left(\sigma -\rho \frac{\partial \psi }{\partial \epsilon }\right):{\dot{\epsilon }}^e+\sigma :{\dot{\epsilon }}^p-\rho \left(s+\frac{\partial \psi }{\partial T}\right)\dot{T}-\rho \frac{\partial \psi }{\partial V_k}\Delta {\dot{V}}_{\text{k}}-\frac{q\Delta gradT}{T}\ge 0$ and plastic strains, with uniform temperatures. We then have zero dissipation, retrieving the classic state equations (13). In the sequel it will be assumed that these state equations remain valid even under irreversible inelastic evolutions, so that the CD inequality becomes:

$$V$$ E16

Under a simplified framework, we require the mechanical and thermal dissipations to be separately non-negative (this reduces the amount of coupling to account for in the model):

$$\text{Φ}={\text{Φ}}_{\text{M}}+{\text{Φ}}_{\text{T}}=\sigma :{\dot{\epsilon }}^p-A_k\Delta {\dot{V}}_{\text{k}}-\frac{q\Delta gradT}{T}\ge 0$$ E17

The thermodynamic force ${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^p-A_k\Delta {\dot{V}}_{\text{k}}\ge 0;{\text{Φ}}_{\text{T}}=-\frac{q\Delta gradT}{T}\ge 0$ , the conjugate variable to the thermodynamic flux $A_k$ , is “defined” as:

$$V_k$$ E18

In practice, $A_k=\rho \frac{\partial \psi }{\partial V_k}$ is often the variable which determines the size (isotropic hardening) or the amount of translation (kinematic hardening) of the yield surface and represents in a simplified manner all the effects of the loading history. One particular example is the pre-consolidation pressure which determines the current yield envelope of clays (as in Camclay model).

The non-negativity of the thermal dissipation can be satisfied by the classic Fourier Law:

$$A_k$$ E19

where the thermal conductivity tensor $q=-K\Delta gradT$ must be symmetric and strictly positive, so that:

$$K$$ E20

It remains to satisfy the non-negativity of the mechanical (or intrinsic) dissipation:

$${\text{Φ}}_{\text{T}}=\frac{gradT\Delta K\Delta gradT}{T}\ge 0$$ E21

The non-negativity of the mechanical dissipation forms the basis for the construction of the material behavioral laws. Note that the equation ${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^p-A_k\Delta {\dot{V}}_{\text{k}}\ge 0$ only “defines” the variable $A_k=\rho \frac{\partial \psi }{\partial V_k}$ but does not contain any rule to calculate its evolution. Similarly, we need a rule to calculate the plastic strain rate $A_k$ .

9. Onsager’s principle

In many physical problems, the total dissipation can be written as the sum of the products between a set of thermodynamic forces ${\dot{\epsilon }}^p$ and theirs conjugates, the thermodynamic flux $X$ :

$$x$$ E22

Onsager, based on theoretical studies at molecular scales where all phenomena are reversible, suggested when the physical process only deviates slightly from the thermodynamic equilibrium, the thermodynamic forces and flux can be related by a set of phenomenological coefficients:

$$\text{Φ}=X\Delta x=X_i{x}_i\ge 0$$ E23

Onsager showed theoretically that the coefficients $X_i=L_{ij}{x}_j$ must be symmetrical. To ensure the non-negativity of the dissipation, it suffices to require $L_{ij}$ to be definite positive, other than being symmetrical. The off-diagonal coefficients allow to account for cross-couplings. This formulation seems to be better suited to moderately non-linear problems. For example, it cannot lead to the classical plastic flow rule in solids.

10. Dissipation potentials

Another, more general, way to satisfy automatically the non-negativity of $L_{ij}$ is to introduce dissipation potentials. This can also handle more general non linear behaviours.

In the case of inelastic behaviour, we define a scalar function called the dissipation potential ${\text{Φ}}_{\text{M}}$ , convex and continuously differentiable with respect to both its arguments, positive everywhere and null at the origin, such that:

$$\phi ({\dot{\epsilon }}^p,{\dot{V}}_{\text{k}})$$ E24

We get immediately:

$$\sigma =\frac{\partial \phi }{\partial {\dot{\epsilon }}^p}{A}_k=-\frac{\partial \phi }{\partial {\dot{V}}_{\text{k}}}$$ E25

In general, it is more convenient to work with ${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^p-A_k\Delta {\dot{V}}_{\text{k}}=\frac{\partial \phi }{\partial {\dot{\epsilon }}^p}:{\dot{\epsilon }}^p+\frac{\partial \phi }{\partial {\dot{V}}_{\text{k}}}\Delta {\dot{V}}_{\text{k}}\ge \phi \ge 0$ , the Legendre transform of ${\phi }^{*}(\sigma ,A_k)$ , also convex and positive definite with respect to its arguments, zero at origin, with:

$$\phi$$ E26

So that:

$${\dot{\epsilon }}^p=\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }{\dot{V}}_{\text{k}}=-\frac{\partial {\phi }^{\text{*}}}{\partial A_k}$$ E27

Theoretically, once the free energy and the dissipation function are specified, the stress-strain relation is fully defined. This is therefore one possible way to construct a constitutive model. However the above reasoning does not work for plasticity.

11. Hardening plasticity for “standard” materials

In plasticity, the dissipation potential is not differentiable. Classically, the usual way to satisfy the dissipation inequality is to define a yield function:

$${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^p-A_k\Delta {\dot{V}}_{\text{k}}=\sigma :\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }+A_k\Delta \frac{\partial {\phi }^{\text{*}}}{\partial A_k}\ge {\phi }^{\text{*}}\ge 0$$ E28

(1) convex with respect to its arguments

(2) the “elastic domain” $F=F(\sigma ,A_k)$ contains the origin, and that:

$$F\left(\sigma ,A_k\right)\le 0$$ E29

where ${\dot{\epsilon }}^p=\dot{\lambda }\frac{\partial F}{\partial \sigma }{\dot{V}}_{\text{k}}=-\dot{\lambda }\frac{\partial F}{\partial A_k}\dot{\lambda }\ge 0$ is the classic plastic multiplier, which obeys the conditions that:

$$\lambda$$ E30

The first condition says if either the stress point is strictly inside the yield surface or if it is currently on the yield surface but moves inwards, the plastic multiplier, hence the plastic strain rate is null. The second condition expresses the condition of plastic loading when the current stress point is on the yield surface and it moves outwards. In this latter case, we have:

$$\dot{\lambda }=0\text{if}F0or\dot{F}0\dot{\lambda }\ge 0\text{if}F=0\text{and}\dot{F}=0$$ E31

The non-negativity of the term between the parenthesis, namely:

$${\text{Φ}}_{\text{M}}=\dot{\lambda }\left(\sigma :\frac{\partial F}{\partial \sigma }+A_k\Delta \frac{\partial F}{\partial A_k}\right)\ge 0$$ E32

stems from geometric arguments (figure 3). This, together with $\left\{\begin{array}{c}\sigma \\ A_k\end{array}\right\}\Delta \left\{\begin{array}{c}\frac{\partial F}{\partial \sigma }\\ \frac{\partial F}{\partial A_k}\end{array}\right\}\ge 0$ , allows to ensure the non-negativity of $\dot{\lambda }\ge 0$ .

Figure 3.

The convex elastic domain contains the origin. Hence the position vector of a point on the boundary ${\text{Φ}}_{\text{M}}$ and the normal vector at the same point $\{\sigma ,A_k\}$ give a positive scalar product.

To construct an elastoplastic model, we need to define a hardening rule:

$$\{{\partial }_{\sigma }F,{\partial }_{A_k}F\}$$ E33

The plastic multiplier $A_k=A_k(V_k)$ can then be determined by the classic consistency condition:

$$\dot{\lambda }$$ E34

For stress-controlled evolutions, this yields, after a little substitution:

$$\dot{F}=\frac{\partial F}{\partial \sigma }\Delta \dot{\sigma }+\frac{\partial F}{\partial A_k}\Delta {\dot{A}}_{\text{k}}=0$$ E35

$\dot{\lambda }=\frac{\frac{\partial F}{\partial \sigma }\Delta \dot{\sigma }}{H}H=\frac{\partial F}{\partial A_k}\Delta \frac{\partial A_k}{\partial V_k}\Delta \frac{\partial F}{\partial A_k}$ is known as the hardening or plastic modulus. To relate the stress increment directly to the strain increment via the tangent stiffness tensor, we substitute:

$$H$$ E36

in the above to get:

$$\dot{\sigma }=D_e\Delta (\dot{\epsilon }-{\dot{\epsilon }}^p){\dot{\epsilon }}^p=\dot{\lambda }\frac{\partial F}{\partial \sigma }$$ E37

Restarting with $\dot{\lambda }=\frac{\frac{\partial F}{\partial \sigma }\Delta D_e\Delta \dot{\epsilon }}{H+\frac{\partial F}{\partial \sigma }\Delta D_e\Delta \frac{\partial F}{\partial \sigma }}$ and after some manipulation leads to:

$$\dot{\sigma }=D_e\Delta (\dot{\epsilon }-{\dot{\epsilon }}^p)=D_e\Delta \left(\dot{\epsilon }-\dot{\lambda }\frac{\partial F}{\partial \sigma }\right)$$ E38

Note that the associative flow rule $\dot{\sigma }=D_{ep}\Delta \dot{\epsilon };D_{ep}=\left(D_e-\frac{D_e\Delta \frac{\partial F}{\partial \sigma }\otimes \frac{\partial F}{\partial \sigma }\Delta D_e}{H+\frac{\partial F}{\partial \sigma }\Delta D_e\Delta \frac{\partial F}{\partial \sigma }}\right)$ renders the tangent matrix ${\dot{\epsilon }}^p=\dot{\lambda }\frac{\partial F}{\partial \sigma }$ symmetric. This relation is also essential in the model construction to ensure the non-negativity of the dissipation. If we replace $D_{ep}$ by ${\dot{\epsilon }}^p=\dot{\lambda }\frac{\partial F}{\partial \sigma }$ with ${\dot{\epsilon }}^p=\dot{\lambda }\frac{\partial g}{\partial \sigma }$ (non-associative flow rule), the CD inequality will no longer be automatically verified. This means that thermodynamic principles may then be violated in some evolutions.Note that in order to describe isotropic and kinematic hardening, the thermodynamic flux $g\ne F$ is often decomposed into a tensor $V_k$ and a scalar $\alpha$ , associated with thermodynamic forces $r$ and $X$ . We would then have to write:

$$R$$ E39

A common example is to identify $F=F\left(\sigma ,X,R\right){\dot{\epsilon }}^p=\dot{\lambda }\frac{\partial F}{\partial \sigma }\dot{\alpha }=-\dot{\lambda }\frac{\partial F}{\partial X}\dot{r}=-\dot{\lambda }\frac{\partial F}{\partial R}$ with the cumulated plastic deviatoric strain $r$ , defined as:

$${\gamma }^p$$ E40

where

$$r={\gamma }^p=\underset{0}{\overset{t}{{\displaystyle \int }}}{\left(\frac{2}{3}{\dot{e}}^p\left(\tau \right):{\dot{e}}^p\left(\tau \right)\right)}^{1/2}d\tau$$ E41

.

12. Viscoplasticity

We start with:

$${\dot{e}}^p=\text{dev}({\dot{\epsilon }}^p)$$ E42

then go through the same procedure as for plasticity:

$$\epsilon ={\epsilon }^e+{\epsilon }^{vp}$$ E43

and:

$$\text{Φ}=\sigma :\dot{\epsilon }-\rho \left(\dot{\psi }+s\dot{T}\right)-\frac{q\Delta gradT}{T}\ge 0$$ E44

We end up with the same dissipation inequality:

$$\psi =\psi ({\epsilon }^e,T,V_k)$$ E45

the same state equations:

$$\text{Φ}=\left(\sigma -\rho \frac{\partial \psi }{\partial \epsilon }\right):{\dot{\epsilon }}^e+\sigma :{\dot{\epsilon }}^{vp}-\rho \left(s+\frac{\partial \psi }{\partial T}\right)\dot{T}-\rho \frac{\partial \psi }{\partial V_k}\Delta {\dot{V}}_{\text{k}}-\frac{q\Delta gradT}{T}\ge 0$$ E46

the same intrinsic dissipation (we discard the thermal part here):

$$\sigma =\rho \frac{\partial \psi }{\partial \epsilon }s=-\frac{\partial \psi }{\partial T}$$ E47

the same definition for the thermodynamic force ${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^{vp}-A_k\Delta {\dot{V}}_{\text{k}}\ge 0$ conjugate to the thermodynamic flux $A_k$ :

$$A_k=?\frac{??}{?V_k}$$ E48

However, a fundamental difference with plasticity intervenes here. In viscoplasticity, a continuously differentiable dissipation potential, definite positive, convex and contains the origin, can be defined:

$$V_k$$ E49

so that the non-negativity condition can be a priori satisfied:

$${\phi }^{\text{*}}={\phi }^{\text{*}}\left(\sigma ,A_k\right){\dot{\epsilon }}^{vp}=\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }{\dot{V}}_{\text{k}}=-\frac{\partial {\phi }^{\text{*}}}{\partial A_k}$$ E50

As for plasticity, in order to describe isotropic and kinematic hardening, the internal variable ${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^{vp}-A_k\Delta {\dot{V}}_{\text{k}}=\sigma :\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }+A_k\Delta \frac{\partial {\phi }^{\text{*}}}{\partial A_k}\ge {\phi }^{\text{*}}\ge 0$ is often decomposed into a tensor $V_k$ and a scalar $\alpha$ , associated with thermodynamic forces $r$ and $X$ :

$$R$$ E51

The mechanical dissipation inequality then becomes:

$$\psi =\psi \left({\epsilon }^e,T,\alpha ,r\right)X=\frac{\partial \psi }{\partial \alpha }R=\frac{\partial \psi }{\partial r}$$ E52

with the corresponding dissipation potential :

$${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^{vp}-X\Delta \dot{\alpha }-R\dot{r}\ge 0$$ E53

We and up with:

$${\phi }^{\text{*}}={\phi }^{\text{*}}\left(\sigma ,X,R\right){\dot{\epsilon }}^{vp}=\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }\dot{\alpha }=-\frac{\partial {\phi }^{\text{*}}}{\partial X}\dot{r}=-\frac{\partial {\phi }^{\text{*}}}{\partial R}$$ E54

For example, Lemaitre's model with isotropic hardening is based on the following dissipation potential:

$${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^{vp}-X\Delta \dot{\alpha }-R\dot{r}=\sigma :\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }+X\Delta \frac{\partial {\phi }^{\text{*}}}{\partial X}+R\frac{\partial {\phi }^{\text{*}}}{\partial R}\ge 0$$ E55

Where ${\phi }^{\text{*}}\left(\sigma ,R\right)=\frac{K}{N+1}{\left(\frac{{\sigma }_{eq}-R}{K}\right)}^{N+1}\frac{1}{\zeta }$ is considered as a parameter independent of the stress tensor, with:

$$\zeta$$ E56

A differentiation gives:

$${\sigma }_{eq}=\sqrt{\frac{3}{2}s:s};s=\text{dev}\left(\sigma \right)=\sigma -\frac{1}{3}tr\left(\sigma \right)I$$ E57

and:

$${\dot{\epsilon }}^{vp}=\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }={\left(\frac{{\sigma }_{eq}-R}{K}\right)}^N\frac{1}{\zeta }\left(\frac{3}{2}\frac{s}{{\sigma }_{eq}}\right)$$ E58

where we have used the identity $\dot{r}=-\frac{\partial {\phi }^{\text{*}}}{\partial R}={\left(\frac{{\sigma }_{eq}-R}{K}\right)}^N\frac{1}{\zeta }$ . Note that the viscoplastic strain rate is purely deviatoric, in other words $\frac{\partial {\sigma }_{eq}}{\partial \sigma }=\frac{3}{2}\frac{s}{{\sigma }_{eq}}$ . Using the classic definition of the equivalent deviatoric viscoplastic strain rate:

$$tr({\dot{\epsilon }}^{vp})=0$$ E59

It can easily be verified that:

$${\dot{\gamma }}_{vp}=\sqrt{\frac{2}{3}{\dot{\epsilon }}^{vp}:{\dot{\epsilon }}^{vp}}$$ E60

${\dot{\gamma }}_{vp}=\dot{r}={\left(\frac{{\sigma }_{eq}-R}{K}\right)}^N\frac{1}{\zeta }$ is an intermediate variable to ensure the consistency of the relations. A particular choice of $\zeta$ can be $\zeta$ which is consistent with the text of Lemaitre & Chabouche (1990). In view of the above identity on $\zeta =r^{N/M}$ and ${\dot{\gamma }}_{vp}$ , we can also write:

$$\dot{r}$$ E61

To define completely the model, we still need a (hardening) relation between ${\dot{\epsilon }}^{vp}={\dot{\gamma }}_{vp}\left(\frac{3}{2}\frac{s}{{\sigma }_{eq}}\right)=\dot{r}\left(\frac{3}{2}\frac{s}{{\sigma }_{eq}}\right)$ et $R$ . This can either be defined explicitly $r$ or by specifying a specific Helmholtz free energy $R=R(r)$ and then uses $\psi$ .

13. Case of biphasic porous media

13.1. Fundamental hypotheses and definitions

In a macroscopic description, a biphasic medium is considered as the superposition of 2 continua. At a given time $R=\frac{\partial \psi }{\partial r}$ and at a given position $t$ , 2 particles, one representing the solid and the other, the fluid, occupy simultaneously the same spatial region $x$ around the geometric point $d{\text{Ω}}_{\text{t}}$ . In order to access separately the mass of each phase, we define the Eulerian porosity $x$ (resp. the Lagrangian porosity $n$ ) so that $\varphi$ (resp. $nd{\text{Ω}}_{\text{t}}$ ) represents the current volume of fluid inside $\varphi d{\text{Ω}}_0$ We have to deal with the macroscopic strain and porosity variations of the solid skeleton. Following Coussy (2004), we split the strain and porosity variation into a elastic and a plastic part:

$$d{\text{Ω}}_{\text{t}}\text{.}$$ E62

We denote by $?$ and ${?}_s$ the volumetric component of the skeleton strain and that of the solid matrix (i.e. $\epsilon ={\epsilon }^e+{\epsilon }^p\dot{\varphi }={\dot{\varphi }}_e+{\dot{\varphi }}_p\text{Δ}\varphi =\varphi -{\varphi }_0={\varphi }_e+{\varphi }_p$ , etc.), which admit the same decomposition:

$$ϵ=tr\left(\epsilon \right)$$ E63

The global volume change comes from those of the solid matrix and of the porous space. It can be proved that:

$$ϵ={ϵ}^e+{ϵ}^p;{ϵ}_s={ϵ}_s^e+{ϵ}_s^p$$ E64

Extending equation (5) to include the contributions of the fluid, we write:

$$ϵ=\left(1-{\varphi }_0\right){ϵ}_s+\varphi -{\varphi }_0{ϵ}^e=\left(1-{\varphi }_0\right){ϵ}_s^e+{\varphi }_e{ϵ}^p=\left(1-{\varphi }_0\right){ϵ}_s^p+{\varphi }_p$$ E65

where $\begin{array}{l}\frac{d^s}{dt}\underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}\left(1-n\right){\rho }_s{s}_{s}d{\text{Ω}}_{\text{t}}+\frac{d^f}{dt}\underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}n{\rho }_f{s}_{f}d{\text{Ω}}_{\text{t}}\ge \underset{\Delta }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}\frac{r_q}{T}d{\text{Ω}}_{\text{t}}-\underset{\partial }{\overset{{\text{Ω}}_{\text{t}}}{{\displaystyle \int }}}\frac{q\Delta n}{T}d\text{S}\\ \end{array}$ express the kinematics of the solid skeleton and fluid phases respectively while $\frac{d^s}{dt}\left(\Delta \right),\frac{d^f}{dt}\left(\Delta \right)$ denote the respective density and entropy of the solid and fluid phases. The Clausius-Duhem inequality corresponding to deformable porous thus admits the following:

$${\rho }_s,s_s,{\rho }_f,s_f$$ E66

where $\text{Φ}={\text{Φ}}_M+{\text{Φ}}_F+{\text{Φ}}_T\ge 0$ are as before the ïntrinsic mechanical and thermal dissipations while ${\text{Φ}}_M,{\text{Φ}}_T$ is the fluid dissipation. Going through the same procedure as in the case of monophasic media, but considering the contributions of both the solid and fluid phases, each with an independent kinematic field, the Clausius-Duhem inequality can be derived:

$${\text{Φ}}_F$$ E67

$${\text{Φ}}_{\text{M}}=\sigma :\dot{\epsilon }+p\dot{\varphi }-{\text{<mover accent="true"> <mi>Ψ</mi> <mo>˙</mo> </mover>}}_{\text{s}}\ge 0$$ E68

where ${\text{Φ}}_F=\left(-gradp+{\rho }_f\left(f^{}-{\gamma }^f\right)\right)\Delta \mathcal{V}\ge 0$ represents the body and inertia forces of the fluid; ${\rho }_f\left(f^{}-{\gamma }^f\right)$ is the filtration vector and $\mathcal{V}=n\left(V^f-V^s\right)$ is the velocity of the fluid phase relative to the solid phase. Introduce the Gibb's free energy $\left(V^f-V^s\right)$ leads to:

$$G_s={\text{Ψ}}_s-p(\varphi -{\varphi }_0)={\text{Ψ}}_s-p({\varphi }_e+{\varphi }_p)$$ E69

Restricting to the case of reversible behaviour where the plastic components and the intrinsic dissipation ${\text{Φ}}_{\text{M}}=\sigma :\dot{\epsilon }-({\varphi }_e+{\varphi }_p)\dot{p}-{\text{<mover accent="true"> <mi>G</mi> <mo>˙</mo> </mover>}}_{\text{s}}\ge 0$ vanish, so that the above inequality becomes an equality, we deduce that ${\text{Φ}}_{\text{M}}$ , and get the state equations:

$$G_s=G_s\left({\epsilon }^e,p\right)$$ E70

Differentiating the above leads to the following constitutive equations:

$$\sigma =\frac{\partial G_s}{\partial {\epsilon }^e}{\varphi }_e=-\frac{\partial G_s}{\partial p}$$ E71

with:

$$d{\sigma }_{ij}=C_{ijmn}d{\epsilon }_{mn}^e-b_{ij}dp;d{\varphi }_e=b_{ij}d{\epsilon }_{ij}^e+\frac{1}{N}dp$$ E72

For isotropic behaviour, we have:

$$C_{ijmn}=\frac{{\partial }^2{G}_s}{\partial {\epsilon }_{ij}^e\partial {\epsilon }_{mn}^e}{b}_{ij}=-\frac{{\partial }^2{G}_s}{\partial {\epsilon }_{ij}^e\partial p}\frac{1}{N}=-\frac{{\partial }^2{G}_s}{\partial p^2}$$ E73

The first of the above equations can be rewritten to introduce an elastic effective stress $a_{\mathit{ij}}^{\text{'}}$ which determines entirely the strain increments under elastic behaviour:

$$d{\sigma }_{ij}=\left(K-\frac{2}{3}G\right)d{\epsilon }_{kk}^e{\delta }_{ij}+2Gd{\epsilon }_{ij}^e-bdp{\delta }_{ij}d{\varphi }_e=bd{ϵ}^e+\frac{1}{N}dp$$ E74

Recalling the following relation resulting from fluid mass conservation and the definition of fluid bulk modulus $d{\sigma }_{ij}^{\text{'}}=\left(K-\frac{2}{3}G\right)d{\epsilon }_{kk}^e{\delta }_{ij}+2Gd{\epsilon }_{ij}^e{\sigma }_{ij}^{\text{'}}={\sigma }_{ij}+bp{\delta }_{ij}$ :

$$K_f$$ E75

Recalling the definition of fluid volume content (neglecting 2nd order terms) $\text{d}\varphi =\frac{d{m}_f}{{\rho }_f}-\varphi \frac{dp}{K_f}$ and combining with the 2nd state equation, we obtain:

$$d{v}_f=\frac{d{m}_f}{{\rho }_f}$$ E76

To introduce a simple non linear skeleton behaviour, we restart with $d{v}_f=d{\varphi }_p+bd{ϵ}^e+\frac{1}{M}dp\frac{1}{M}=\frac{1}{N}+\frac{\varphi }{K_f}$ , and postulates that:

$${\text{Φ}}_{\text{M}}=\sigma :\dot{\epsilon }+p\dot{\varphi }-{\text{<mover accent="true"> <mi>Ψ</mi> <mo>˙</mo> </mover>}}_{\text{s}}\ge 0$$ E77

Where ${\Psi }_s\left({\epsilon }^e,{\varphi }_e,V_k\right)=W_s\left({\epsilon }^e,{\varphi }_e\right)+U\left(V_k\right)$ represents the trapped energy due to hardening, depending only on the internal state parameters $U\left(V_k\right)$ . Substituting this into the Clausius-Duhem inequality and simplifying leads to:

$$V_k$$ E78

with:

$${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^p+p{\dot{\varphi }}_p+A_k\Delta {\dot{V}}_k\ge 0$$ E79

The above inequality can also be rewritten as:

$$\sigma =\frac{\partial {\text{Ψ}}_s}{\partial {\epsilon }^e}=\frac{\partial {\text{W}}_s}{\partial {\epsilon }^e}p=\frac{\partial {\text{Ψ}}_s}{\partial {\varphi }_e}=\frac{\partial {\text{W}}_s}{\partial {\varphi }_e}{A}_k=-\frac{\partial {\text{Ψ}}_s}{\partial V_k}=-\frac{\partial U}{\partial V_k}$$ E80

Hence $\delta {\text{Φ}}_{\text{M}}=\delta W^p-dU\ge 0\delta W^p=\sigma :d{\epsilon }^p+pd{\varphi }_{p}dU=\frac{\partial U}{\partial V_k}d{V}_k$ represents that part of the plastic work which is not dissipated into heat. Returning to (65), it is observed that the non-negativity of the dissipation $dU$ leads to Darcy’s law as the constitutive equation of flow, which is defined for the isotropic case as:

$${\text{Φ}}_F$$ E81

where $\text{n}\left(V^f-V^s\right)={\text{λ}}_{\text{h}}\left(-gradp+{\rho }_f\left(f-{\gamma }^f\right)\right)$ is the hydraulic conductivity or coefficient of permeability of the medium. It is interesting to note that the thermodynamic approach confirms Darcy’s law governs fluid flow relative to the solid matrix, and not with respect to a stationary observer.

13.2. Poroplastic behaviour

As for monophasic media, the dissipation potential is not differentiable in plasticity. To satisfy the non-negativity of the intrinsic dissipation, we postulate an elastic domain defined by a convex function ${\lambda }_h$ :

$$f$$ E82

The domain contains the origin, in other words:

$$F\left(\sigma ,p,A_k\right)\le 0$$ E83

Introducing the classic standard material behavioural law:

$$F\left(0,0,0\right)0$$ E84

we have:

$$d{\epsilon }^p=d\lambda \frac{\partial F}{\partial \sigma }d{\varphi }_p=d\lambda \frac{\partial F}{\partial p}d{V}_k=d\lambda \frac{\partial F}{\partial A_k}d\lambda \ge 0;F\le 0$$ E85

The quantity between square brackets represents the scalar product between the position vector ${\text{Φ}}_{\text{M}}=d\lambda \left[\sigma :\frac{\partial F}{\partial \sigma }+p\frac{\partial F}{\partial p}+A_k\frac{\partial F}{\partial A_k}\right]\ge 0$ and the outward normal vector $\left(\sigma ,p,A_k\right)$ which is perpendicular to the boundary of the elastic domain $\left(\frac{\partial F}{\partial \sigma },\frac{\partial F}{\partial p},\frac{\partial F}{\partial A_k}\right)$ . Its positivity comes from the geometric convexity of the domain $f=0$ and the fact that the domain contains the origin. In the above formulation, the yield criterion is supposed to depend both on the total stress and the fluid pressure. This can be simplified if the plastic porosity change is related to the plastic volumetric strain:

$$F\le 0$$ E86

so that:

$${\dot{\varphi }}_p=b^{\text{'}}{\dot{ϵ}}^p=b^{\text{'}}{\dot{\epsilon }}^p:\text{I}$$ E87

Mechanical stress and fluid pressure then intervene in the yield function only via a plastic effective stress ${\text{Φ}}_{\text{M}}=\sigma \text{''}:{\dot{\epsilon }}^p-{\text{<mover accent="true"> <mi>Ψ</mi> <mo>˙</mo> </mover>}}_{\text{s}}\ge 0\sigma \text{''}=\sigma +b^{\text{'}}p\text{I}$ :

$$\sigma \text{'}\text{'}$$ E88

with:

$$F\left(\sigma \text{''},A_k\right)\le 0$$ E89

However, there are two effective stresses $d{\epsilon }^p=d\lambda \frac{\partial F}{\partial {\sigma }^{\text{''}}}d{V}_k=d\lambda \frac{\partial F}{\partial A_k}d\lambda \ge 0F\le 0$ and $\sigma \text{'}$ , which is confusing. The situation will be optimum if we can assume either $\sigma \text{'}\text{'}$ , hence $b^{\text{'}}=b$ , or matrix incompressibility which implies ${\sigma }^{\text{'}}={\sigma }^{\text{'}\text{'}}$ and that $b^{\text{'}}=b=1$ . The last case is of particular importance and corresponds to the majority of cases in soils.The above flow rule is known as associative since the strain rate is normal to the yield surface, with the advantage that the non-negativity of the dissipation is always satisfied. Geomaterials exhibit complex volumetric behaviours and sometimes call for non associative flow rules:

$${\sigma }^{\text{'}}={\sigma }^{\text{'}\text{'}}=\sigma +pI$$ E90

However, the non-negativity of the dissipation is not always satisfied in this last case.

13.3. Poroviscoplastic behaviour

Recall that we have to satisfy:

$$d{\epsilon }^p=d\lambda \frac{\partial g}{\partial \sigma \text{''}}d{V}_k=d\lambda \frac{\partial g}{\partial A_k}d\lambda \ge 0F\le 0$$ E91

The dissipation potential is in this case differentiable so that we can write:

$${\text{Φ}}_{\text{M}}=\sigma :{\dot{\epsilon }}^p+p{\dot{\varphi }}_p+A_k{\stackrel{\Delta V̇}{\Delta }}_k\ge 0$$ E92

Hence:

$${\phi }^{\text{*}}={\phi }^{\text{*}}\left(\sigma ,p,A_k\right){\dot{\epsilon }}^p=\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }{\dot{\varphi }}_p=\frac{\partial {\phi }^{\text{*}}}{\partial p}{\dot{V}}_k=\frac{\partial {\phi }^{\text{*}}}{\partial A_k}$$ E93

Similar to the case of plasticity, we can simplify by supposing ${\text{Φ}}_{\text{M}}=\sigma :\frac{\partial {\phi }^{\text{*}}}{\partial \sigma }+p\frac{\partial {\phi }^{\text{*}}}{\partial p}+A_k\Delta \frac{\partial {\phi }^{\text{*}}}{\partial A_k}\ge 0$ and ${\dot{\varphi }}_p=b{\dot{ϵ}}^p=b{\dot{\epsilon }}^p:I$ . We then require the dissipative potential to satisfy:

$$\sigma \text{'}=\sigma +bpI$$ E94

For example, if we take:

$${\phi }^{\text{*}}={\phi }^{\text{*}}\left(\sigma \text{'},A_k\right){\dot{\epsilon }}^p=\frac{\partial {\phi }^{\text{*}}}{\partial \sigma \text{'}}{\dot{V}}_k=\frac{\partial {\phi }^{\text{*}}}{\partial A_k}$$ E95

We get:

$${\phi }^{\text{*}}=\frac{1}{2\eta }F{\left(\sigma \text{'},A_k\right)}^2$$ E96

14. Applications

14.1. Example 1 – Hardening plasticity – EPS geofoam

In the following example we illustrate the first type of use of the second thermodynamic principle discussed in Section 6, namely, by verifying a constitutive model of EPS geofoam a posteriori for thermodynamic consistency. This model was developed by the authors (Wong and Leo, 2006) based on experimental results from a series of standard “drained” triaxial tests. It initially adopted the Mohr-Coulomb yield function used widely in soil mechanics but upon further testing with a true triaxial apparatus (Leo et al., 2008), a Drucker-Prager type yield function was subsequently preferred. This is written as:

$${\dot{\epsilon }}^p=\frac{\partial {\phi }^{\text{*}}}{\partial \sigma \text{'}}=\frac{1}{\eta }F\frac{\partial F}{\partial \sigma \text{'}}{\dot{V}}_k=\frac{\partial {\phi }^{\text{*}}}{\partial A_k}=\frac{1}{\eta }F\frac{\partial F}{\partial A_k}$$ E97

i.e.

$$F\left(\sigma ,a\right)=\sqrt{3J_2}-b{I}_1-a=0$$ E98

where $dF=\frac{\partial F}{\partial \sigma }\Delta \sigma +\frac{\partial F}{\partial a}\Delta da=0$ is the first stress invariant, $I_1=tr\left(\sigma \right)$ is the second stress invariant and b is a material constant. Here $J_2=\frac{1}{2}s:s$ is the hardening law accounting for the isotropic hardening effects; $a\left(r\right)=a_0+\beta r$ are material constants and $r$ is an internal variable chosen as the equivalent deviatoric plastic strain defined by:

$$a_0,\beta$$ E99

Referring to the discussion in Section 11, we observe that equation (94) is a particular form of (27), $r=\underset{0}{\overset{t}{{\displaystyle \int }}}\sqrt{\frac{1}{2}{\dot{e}}^p:{\dot{e}}^p}d\tau$ of $a\left(r\right)$ , and (96) is the equivalent of (39). Geometrically, the surface of equation (94) corresponds to a conical surface, with the symmetry axis coinciding with the hydrostatic axis. The apex angle is governed entirely by the constant b, whereas a, together with b, determines the distance separating the cone tip from the origin. According to the laws of thermodynamics, an associative flow rule should have been adopted for the plastic strain (i.e. $A_k(V_k)$ in equation (28)) for this constitutive model, but we chose a non-associative flow rule instead where,

$${\dot{\epsilon }}^p=\dot{\lambda }\frac{\partial F}{\partial \sigma }$$ E100

c is a rheological parameter which depends on the initial stress. This is because experimental measurements suggest that the plastic volumetric strain is better represented by the plastic potential given in (97) rather than the yield function of (94). As discussed earlier, this means that the thermodynamics principle in terms of the non-negativity of the dissipation may possibly be violated in some evolutions since the normality rule (plastic strain increment being normal to the yield function) is not being followed. The associative flow rule, however, has been a problem with some geomaterials such as soils and rocks in that it tends to erroneously predict plastic volumetric strain. This is one instance where the insight provided by thermodynamics into post yielding volumetric behavior is seemingly at odds with experimental evidence. In these cases it is widely accepted that the plastic volumetric behavior would be better captured using a non-associative flow rule. These cases also demonstrate that while thermodynamics insights provide useful guidance to help engineers focus on important aspects of the constitutive relationships in continuum mechanics, it is necessary that these insights should ultimately be supported by experimental evidence.

14.2. Example 2 – Poroelasticity: closure of a spherical cavity

This example dealing with the closure of a deeply embedded cavity in poroelastic medium was previously studied by the authors (Wong et al. 2008). Here we illustrate the second type of use of the second thermodynamic principle discussed in Section 6, where the thermodynamics concepts from Section 13.1 are applied to formulate the constitutive relationships that lead, importantly, to the analytical solutions for the closure of a spherical cavity. The closure constitutes part of a life cycle of an underground mining cavity idealised by four stages. Initially, the ground is in a state of hydro-mechanical equilibrium. The cavity is then excavated and an internal support is provided to maintain its stability. Various techniques of support exist. For example, it can be evenly spaced steel bolts or a layer of shotcrete or a combination of them. For modelling purposes, this support can be assimilated to a layer of elastic material lining the cavity walls. At the end of its service life, the cavity is backfilled with a poro-elastic material before being abandoned. We were interested in the long term evolution of the hydro-mechanical fields in the surrounding medium and in the backfill after the its abandonment, when the support starts to deteriorate. This problem deals with a special case of thereversible behaviour where the intrinsic dissipation vanishes, namely ${\dot{\epsilon }}^p=\dot{\lambda }\frac{\partial G}{\partial \sigma }G\left(\sigma \right)=\sqrt{3J_2}-c{I}_1$ (as opposed to the more general case of irreversible behaviour for materials with plasticity and/or viscosity), leading to the state equation (67) and the constitutive equations (70) for isotropic poroelastic material. Limiting ourselves to small strains, we define:

$${\text{Φ}}_M=0$$ E101

where $d{\sigma }_{ij}={\sigma }_{ij}-{\sigma }_{ij}^{0}d{\text{ε}}_{ij}^{\text{e}}={\epsilon }_{ij}^e-{\epsilon }_{ij}^{0}dp=p-p^0$ denotes the initial stress, strain, pore pressure respectively. We make further assumptions that the solid grains of the medium are incompressible, and it thus holds that the skeletal volumetric change ${\sigma }_{ij}^0,{\epsilon }_{ij}^0,p^0$ must be the same as the change in the porosity $d{ϵ}^e=d{\epsilon }_{kk}^e$ , that is:

$$d{\varphi }^e$$ E102

By comparing (99) to the second equation of (70), it is evident that the values of Biot coefficients must be: $d{ϵ}^e=d{\varphi }^e$ and $$b=\text{ 1}$$ . Taking initial strain $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$N$}\right.=0$ , equation (70) thus yields the following constitutive relationships for a linear isotropic poroelastic material:

$${\epsilon }_{ij}^0=0$$ E103

Since we are clearly dealing with a poroelastic medium, the superscript ‘e ‘ denoting elastic strain shall be omitted in Example 2 without ambiguity, for the sake of brevity. For the fluid phase of the porous material, the constitutive equation follows from the thermodynamically consistent Darcy’s law, equation (78). Here, after neglecting inertia effects but not body forces due to gravity g, the fluid mass flux, ${\sigma }_{ij}-{\text{σ}}_{\text{ij}}^0=\left(K-\frac{2}{3}G\right){ϵ}^{\text{e}}{\delta }_{ij}+2G{\epsilon }_{ij}^e-\left(p-p^0\right){\delta }_{ij}\varphi -{\varphi }_0={ϵ}^e$ is related to the thermodynamic forces as: $w_f={\rho }_{f}n\left(V^f-V^s\right)$ . At $\raisebox{1ex}{$w_f$}\!\left/ \!\raisebox{-1ex}{${\rho }_f$}\right.={\lambda }_h\left(-gradp+{\rho }_{f}g\right)$ , the fluid is assumed to be in hydraulic equilibrium, implying that: $$t=0$$ . The difference between these two equations yields:

$$0={\lambda }_h\left(-grad{p}_0+{\rho }_{f}g\right)$$ E104

As shown above, insights from thermodynamics principles have lead to constitutive equations (100) and (101). These equations thus allow us to develop a set of governing equationswhich is applicable to the cavity closure problem. These equations are then solved with respect to the initial and boundary conditions for a spherical cavity to obtain a set of analytic solutions, of which a detailed discussion is given in Wong et al. (2008).

14.3. Example 3 – Poroviscoelasticity: closure of long cylindrical tunnel

Example 3 illustrates the use of thermodynamics principles in formulating constitutive equations for a poro-viscoelastic medium. The ultimate purpose here is also to develop solutions for a long horizontally aligned tunnel with a circular cross-section embedded in a poro-viscoelastic massif. The setting of the problem is similar to Example 2 discussed above except that the spherical cavity is replaced by a long lined tunnel (Dufour et al. 2009). We start by restricting to small strain problems where the strain tensor of a viscoelastic material can be decomposed into an elastic part (denoted by superscript ‘e ’) and a viscoelastic part (superscript ‘ν’):

$$\raisebox{1ex}{$w_f$}\!\left/ \!\raisebox{-1ex}{${\rho }_f$}\right.=-{\lambda }_{h}grad\left(p-p_0\right)$$ E105

The strain and stress tensors are separated into isotropic and deviatoric parts as follow:

$${\epsilon }_{ij}={\epsilon }_{ij}^e+{\epsilon }_{ij}^{\nu }$$ E106

where ${\epsilon }_{ij}=\frac{1}{3}ϵ{\delta }_{ij}+e_{ij}{\sigma }_{ij}=\sigma {\delta }_{ij}+s_{ij}$ are the mean and deviatoric strains defined previously; $ϵ,e_{ij}$ is the mean stress and $\sigma =\raisebox{1ex}{${\sigma }_{ij}$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$ is the deviatoric stress tensor. It is noted that the decomposition into elastic and viscoelastic parts in (102) apply separately to $s_{ij}={\sigma }_{ij}-\sigma {\delta }_{ij}$ , $ϵ$ and the porosity as well such that:

$$e_{ij}$$ E107

Correspondence between volumetric strain and porosity change holds for each of the elastic and viscoelastic components:

$$ϵ={ϵ}^e+{ϵ}^{\nu }{e}_{ij}=e_{ij}^e+e_{ij}^{}\varphi -{\varphi }_0={\varphi }_e+{\varphi }_{\nu }$$ E108

14.3.1. Poroviscoelastic constitutive equations

Following (74), we postulate the existence of trapped energy due to viscosity that depends on viscous strains only and write the free energy of the skeleton as:

$$ϵ=\varphi -{\varphi }_0{ϵ}^e={\varphi }_e{ϵ}^{\nu }={\varphi }_{\nu }$$ E109

where the following relationships are considered for the functions ${\Psi }_s\left({ϵ}^e,e_{ij}^e,{\varphi }_{\text{e}},{ϵ}^{\text{ν}},e_{ij}^{\nu }\right)=W_s\left({ϵ}^e,e_{ij}^e,{\varphi }_{\text{e}}\right)+U\left({ϵ}^{\text{ν}},e_{ij}^{\nu }\right)$ :

$$W_s,U$$ E110

$$W_s\left({ϵ}^e,e_{ij}^e,{\varphi }_{\text{e}}\right)=\frac{1}{2}{K}_0{ϵ}^e{}^2+{\mu }_0{e}_{ij}^e{}^2+\frac{1}{2}{N}_0{\left({ϵ}^e-{\varphi }_e\right)}^2$$ E111

Specialising to a linear isotropic porous material, after substituting (107) into (76) and taking into consideration the decomposition into the mean and deviatoric parts, and the initial stresses we get:

$$U\left({ϵ}^{\text{ν}},e_{ij}^{\nu }\right)=\frac{1}{2}\xi {ϵ}^{\nu }{}^2+\chi e_{ij}^{\nu }{}^2$$ E112

$$\left(\sigma -{\sigma }_0\right)+\left(p-p_0\right)=K_0\left(ϵ-{ϵ}^{\nu }\right)$$ E113

$$s_{ij}-s_{ij}^0=2{\mu }_0\left(e_{ij}-e_{ij}^{\nu }\right)$$ E114

K ₀ , μ ₀, N ₀are the initial or “short term” analogues of K, μ, Nrespectively. Further substitution of (106) – (111) into (64) yields:

$$p-p_0=-N_0\left({ϵ}^e-{\varphi }_e\right)$$ E115

In the next step, a convex dissipative potential $\left\{\left(\sigma -{\sigma }_0\right)+\left(p-p_0\right)-\xi {ϵ}^{\nu }\right\}{\dot{ϵ}}^{\nu }+\left\{\left(s_{ij}-s_{ij}^0\right)-2\chi e_{ij}^{\nu }\right\}{\dot{e}}_{ij}^{\nu }\ge 0$ is introduced so that based on (112):

$$\mathcal{D}\left({\dot{ϵ}}^{\nu },{\dot{e}}_{ij}^{\nu }\right)$$ E116

which leads to:

$$\frac{\partial \mathcal{D}}{\partial {\dot{ϵ}}^{\nu }}=\left(\sigma -{\sigma }_0\right)+\left(p-p_0\right)-\xi {ϵ}^{\nu }\frac{\partial \mathcal{D}}{\partial {\dot{e}}_{ij}^{\nu }}=\left(s_{ij}-s_{ij}^0\right)-2\chi e_{ij}^{\nu }$$ E117

where positivity of $\mathcal{D}\left({\dot{ϵ}}^{\nu },{\dot{e}}_{ij}^{\nu }\right)=\frac{1}{2}\zeta {\dot{ϵ}}^{\nu }{}^2+\eta {\dot{e}}_{ij}^{\nu }{}^2$ ensures the convexity of $\zeta \ge 0;\eta \ge 0$ . From the above developments, the constitutive equations relating stresses to strains for an isotropic poro-viscoelastic material may thus be defined by equations (109)-(111) as well as by the following equations.

$$\mathcal{D}\left({\dot{ϵ}}^{\nu },{\dot{e}}_{ij}^{\nu }\right)$$ E118

$$\left(\sigma -{\sigma }_0\right)+\left(p-p_0\right)=\xi {ϵ}^{\nu }+{\dot{\zeta ϵ}}^{\nu }$$ E119

where $s_{ij}-s_{ij}^0=2\chi e_{ij}^{\nu }+2\eta {\dot{e}}_{ij}^{\nu }$ are rheological constants. Note that these equations have been formulated based on the thermodynamics approach while adopting the convex dissipative potential, $\xi ,\zeta ,\chi ,\eta$ , in equation (114). Before proceeding further, we will now introduce the Laplace transform, defined for a typical function $\mathcal{D}\left({\dot{ϵ}}^{\nu },{\dot{e}}_{ij}^{\nu }\right)$ as follows:

$$f\left(r,t\right)$$ E120

where s is the Laplace transform parameter and $\begin{array}{l}\overline{f}\left(r,s\right)=L\left\{f\left(r,t\right)\right\}=\underset{0}{\overset{\infty }{{\displaystyle \int }}}f\left(r,t\right)e^{-st}dt\\ \overline{f}\left(r,s\right)=L^{-1}\left\{\overline{f}\left(r,t\right)\right\}=\frac{1}{2\pi i}\underset{\Delta }{\overset{\Delta \text{Γ}-i\infty \text{Γ}+i\infty }{{\displaystyle \int }}}\overline{f}\left(r,t\right)e^{st}ds\end{array}$ . In the notations adopted here, the bar over the symbol denotes the transformed function represented by the symbol. The value Γ is chosen such that all poles in the s-plane lie to the left of the vertical line $$i^{\text{2}}=-\text{1}$$ . Taking the Laplace transform of (109), (110), (113) and (116) and solving for the viscous volumetric and deviatoric strains give,

$$\text{Re}\left(\text{s}\right)=\Gamma$$ E121

The constitutive equations (115), (116), (118) are then used to developed governing equations for the closure of a long cylindrical tunnel in poroviscoelastic massif. Laplace transform solutions have been developed and discussed in detail in Dufour et al. (2009) to which interested readers may refer.