InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Chemical Engineering » "Thermodynamics - Systems in Equilibrium and Non-Equilibrium", book edited by Juan Carlos Moreno-Piraján, ISBN 978-953-307-283-8, Published: October 10, 2011 under CC BY 3.0 license. © The Author(s).

Chapter 4

Thermodynamics in Mono and Biphasic Continuum Mechanics

By Henry Wong, Chin J Leo and Natalie Dufour
DOI: 10.5772/20235

Article top


Formula: Eqn013.wmf>is a generic part of a body in movement, with distributed body forces
						, surface tractions 
						inward heat flux
						, and distributed heat source
						is the outward unit normal.
Figure 1. Formula: Eqn013.wmf>is a generic part of a body in movement, with distributed body forces Ω t , surface tractions f inward heat flux t , and distributed heat source q , r q is the outward unit normal.
The heat engine represented by the circle takes a quantity of heat 
						from the hotter reservoir 
						and rejects 
						to the colder reservoir
						, while it performs an amount of useful work
						. The first principle requires 
						and the second principle sets a theoretical upper bound on the efficiency
						 attainable by heat engines.
Figure 2. The heat engine represented by the circle takes a quantity of heat Δ W Δ Q ≤ T 1 − T 2 T 2 from the hotter reservoir ΔQ and rejects T 1 to the colder reservoir ΔQ' , while it performs an amount of useful work T 2 . The first principle requires Δ W and the second principle sets a theoretical upper bound on the efficiency Δ W = ΔQ − ΔQ' attainable by heat engines.
The convex elastic domain contains the origin. Hence the position vector of a point on the boundary 
						 and the normal vector at the same point 
						 give a positive scalar product.
Figure 3. The convex elastic domain contains the origin. Hence the position vector of a point on the boundary Φ M and the normal vector at the same point { σ , A k } give a positive scalar product.

Thermodynamics in Mono and Biphasic Continuum Mechanics

Henry Wong1, Chin J. Leo2 and Natalie Dufour1

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 ρ and a velocity field v (figure 1). Our attention will be focused on an arbitrarily chosen part of this body, which occupies a volume Ω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 Ωt is subject to a distributed body force f (for example gravity) and surface tractions t on its boundary surface, noted Ωt . At the same time, the body is subject to a heat flux q on Ωt and an internal heat source rq .

To begin with, we consider the energy and entropy balance of all the matter inside the volume Ω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 Ωt , surface tractions f inward heat flux t , and distributed heat source q , rq is the outward unit normal.

In the above equation, ddt(E+K)=Px+Q and E are the global internal and kinetic energies, while K and Px 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:

where E=ΔΩtρedΩtK=ΔΩt12ρvΔvdΩtPx=ΔΩtfΔvdΩt+ΩtΔvdSQ=ΔΩtrqdΩtΩqΔndS , 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=σΔn to the second order symmetric stress tensor t , we get after some simplifications:

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:

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

  • 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 ρe˙=σ:ε˙+rqdivq , in the term ρ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 Px into useful work, a necessary condition is to have at least two reservoirs with two different (absolute) temperatures ΔQ and T1 (suppose T2 to fix ideas).

  3. Moreover, the above conversion can at best be partial in that the amount of work eW extractable from a given quantity of heat T1T2 admits a theoretical upper bound depending on the two temperatures:


Figure 2.

The heat engine represented by the circle takes a quantity of heat ΔWΔQT1T2T2 from the hotter reservoir ΔQ and rejects T1 to the colder reservoir ΔQ' , while it performs an amount of useful work T2 . The first principle requires ΔW and the second principle sets a theoretical upper bound on the efficiency ΔW=ΔQΔ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:

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:


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):


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


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 1ρσ:ε˙+Ts˙e˙=0or de=1ρσ:dε+Tds and ε (i.e. s ) and that the following state equations hold:

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:


to recast inequality (7) to the following form:

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


via the same reasoning as previously, we deduce that the specific free energy 1ρσ:ε˙sT˙ψ˙=0or dψ=1ρσ:dεsdT depends naturally on ψ and ε and satisfies the following state equations:

The Legendre transform (10) thus allows one to define a thermodynamic potential with natural independent variables which are more accessible ( σ=ρψεs=ψ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 2nd principle, in an a posteriori manner, in the sense that the construction of the model does not rely in any way on the 2nd 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:

Is assumed, where ε=εe+εp is the elastic strain and ε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 εp , so that ψ=ψ(εe,T,Vk) . The scalar variables grouped into a tensor ψ˙=ψεεė+ψTT˙+ψVkV˙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:

Consider the particular case of elastic (reversible) evolution corresponding to stationary values of the internal variables Φ=(σρψε):ε˙e+σ:ε˙pρ(s+ψT)T˙ρψVkΔV˙kqΔgradTT0 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:

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):


The thermodynamic force ΦM=σ:ε˙pAkΔV˙k0;ΦT=qΔgradTT0 , the conjugate variable to the thermodynamic flux Ak , is “defined” as:

In practice, Ak=ρψVk 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:

where the thermal conductivity tensor q=KΔgradT must be symmetric and strictly positive, so that:

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


The non-negativity of the mechanical dissipation forms the basis for the construction of the material behavioral laws. Note that the equation ΦM=σ:ε˙pAkΔV˙k0 only “defines” the variable Ak=ρψVk but does not contain any rule to calculate its evolution. Similarly, we need a rule to calculate the plastic strain rate Ak .

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 ε˙p and theirs conjugates, the thermodynamic flux X :

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:


Onsager showed theoretically that the coefficients Xi=Lijxj must be symmetrical. To ensure the non-negativity of the dissipation, it suffices to require Lij 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 Lij 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 ΦM , convex and continuously differentiable with respect to both its arguments, positive everywhere and null at the origin, such that:


We get immediately:


In general, it is more convenient to work with ΦM=σ:ε˙pAkΔV˙k=φε˙p:ε˙p+φV˙kΔV˙kφ0 , the Legendre transform of φ*(σ,Ak) , also convex and positive definite with respect to its arguments, zero at origin, with:

So that:


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:


(1) convex with respect to its arguments

(2) the “elastic domain F=F(σ,Ak) contains the origin, and that:

where ε˙p=λ˙FσV˙k=λ˙FAkλ˙0 is the classic plastic multiplier, which obeys the conditions that:

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:


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


stems from geometric arguments (figure 3). This, together with {σAk}Δ{FσFAk}0 , allows to ensure the non-negativity of λ˙0 .


Figure 3.

The convex elastic domain contains the origin. Hence the position vector of a point on the boundary ΦM and the normal vector at the same point {σ,Ak} give a positive scalar product.

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


The plastic multiplier Ak=Ak(Vk) can then be determined by the classic consistency condition:

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

λ˙=FσΔσ˙HH=FAkΔAkVkΔFAk 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:

in the above to get:


Restarting with λ˙=FσΔDeΔε˙H+FσΔDeΔFσ and after some manipulation leads to:


Note that the associative flow rule σ˙=DepΔε˙;Dep=(DeDeΔFσFσΔDeH+FσΔDeΔFσ) renders the tangent matrix ε˙p=λ˙Fσ symmetric. This relation is also essential in the model construction to ensure the non-negativity of the dissipation. If we replace Dep by ε˙p=λ˙Fσ with ε˙p=λ˙gσ (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 gF is often decomposed into a tensor Vk and a scalar α , associated with thermodynamic forces r and X . We would then have to write:

A common example is to identify F=F(σ,X,R)ε˙p=λ˙Fσα˙=λ˙FXr˙=λ˙FR with the cumulated plastic deviatoric strain r , defined as:



12. Viscoplasticity

We start with:


then go through the same procedure as for plasticity:


We end up with the same dissipation inequality:


the same state equations:


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


the same definition for the thermodynamic force ΦM=σ:ε˙vpAkΔV˙k0 conjugate to the thermodynamic flux Ak :

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:

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


As for plasticity, in order to describe isotropic and kinematic hardening, the internal variable ΦM=σ:ε˙vpAkΔV˙k=σ:φ*σ+AkΔφ*Akφ*0 is often decomposed into a tensor Vk and a scalar α , associated with thermodynamic forces r and X :

The mechanical dissipation inequality then becomes:


with the corresponding dissipation potential :


We and up with:


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


Where φ*(σ,R)=KN+1(σeqRK)N+11ζ is considered as a parameter independent of the stress tensor, with:

A differentiation gives:


where we have used the identity r˙=φ*R=(σeqRK)N1ζ . Note that the viscoplastic strain rate is purely deviatoric, in other words σeqσ=32sσeq . Using the classic definition of the equivalent deviatoric viscoplastic strain rate:

It can easily be verified that:

γ˙vp=r˙=(σeqRK)N1ζ is an intermediate variable to ensure the consistency of the relations. A particular choice of ζ can be ζ which is consistent with the text of Lemaitre & Chabouche (1990). In view of the above identity on ζ=rN/M and γ˙vp , we can also write:

To define completely the model, we still need a (hardening) relation between ε˙vp=γ˙vp(32sσeq)=r˙(32sσeq) et R . This can either be defined explicitly r or by specifying a specific Helmholtz free energy R=R(r) and then uses ψ .

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=ψ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Ωt . In order to access separately the mass of each phase, we define the Eulerian porosity x (resp. the Lagrangian porosity n ) so that ϕ (resp. ndΩt ) represents the current volume of fluid inside ϕdΩ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:

We denote by ? and ?s the volumetric component of the skeleton strain and that of the solid matrix (i.e. ε=εe+εpϕ˙=ϕ˙e+ϕ˙pΔϕ=ϕϕ0=ϕe+ϕp , etc.), which admit the same decomposition:

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


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


where dsdtΔΩt(1n)ρsssdΩt+dfdtΔΩtnρfsfdΩtΔΩtrqTdΩtΩtqΔnTdS express the kinematics of the solid skeleton and fluid phases respectively while dsdt(Δ),dfdt(Δ) denote the respective density and entropy of the solid and fluid phases. The Clausius-Duhem inequality corresponding to deformable porous thus admits the following:


where Φ=ΦM+ΦF+ΦT0 are as before the ïntrinsic mechanical and thermal dissipations while ΦM,Φ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:

ΦM=σ:ε˙+pϕ˙ Ψ˙ s0

where ΦF=(gradp+ρf(fγf))ΔV0 represents the body and inertia forces of the fluid; ρf(fγf) is the filtration vector and V=n(VfVs) is the velocity of the fluid phase relative to the solid phase. Introduce the Gibb's free energy (VfVs) leads to:


Restricting to the case of reversible behaviour where the plastic components and the intrinsic dissipation ΦM=σ:ε˙(ϕe+ϕp)p˙ G˙ s0 vanish, so that the above inequality becomes an equality, we deduce that ΦM , and get the state equations:

Differentiating the above leads to the following constitutive equations:




For isotropic behaviour, we have:


The first of the above equations can be rewritten to introduce an elastic effective stress aij' which determines entirely the strain increments under elastic behaviour:


Recalling the following relation resulting from fluid mass conservation and the definition of fluid bulk modulus dσij'=(K23G)dεkkeδij+2Gdεijeσij'=σij+bpδij :

Recalling the definition of fluid volume content (neglecting 2nd order terms) dϕ=dmfρfϕdpKf and combining with the 2nd state equation, we obtain:

To introduce a simple non linear skeleton behaviour, we restart with dvf=dϕp+bdϵe+1Mdp1M=1N+ϕKf , and postulates that:

ΦM=σ:ε˙+pϕ˙ Ψ˙ s0

Where Ψs(εe,ϕe,Vk)=Ws(εe,ϕe)+U(Vk) represents the trapped energy due to hardening, depending only on the internal state parameters U(Vk) . Substituting this into the Clausius-Duhem inequality and simplifying leads to:



The above inequality can also be rewritten as:


Hence δΦM=δWpdU0δWp=σ:dεp+pdϕpdU=UVkdVk 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:

where n(VfVs)=λh(gradp+ρf(fγf)) 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 λh :

The domain contains the origin, in other words:


Introducing the classic standard material behavioural law:

we have:


The quantity between square brackets represents the scalar product between the position vector ΦM=dλ[σ:Fσ+pFp+AkFAk]0 and the outward normal vector (σ,p,Ak) which is perpendicular to the boundary of the elastic domain (Fσ,Fp,FAk) . 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:

so that:


Mechanical stress and fluid pressure then intervene in the yield function only via a plastic effective stress ΦM=σ'':ε˙p Ψ˙ s0σ''=σ+b'pI :



However, there are two effective stresses dεp=dλFσ''dVk=dλFAkdλ0F0 and σ' , which is confusing. The situation will be optimum if we can assume either σ'' , hence b'=b , or matrix incompressibility which implies σ'=σ'' and that b'=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:


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:


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




Similar to the case of plasticity, we can simplify by supposing ΦM=σ:φ*σ+pφ*p+AkΔφ*Ak0 and ϕ˙p=bϵ˙p=bε˙p:I . We then require the dissipative potential to satisfy:

For example, if we take:


We get:


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:




where dF=FσΔσ+FaΔda=0 is the first stress invariant, I1=tr(σ) is the second stress invariant and b is a material constant. Here J2=12s:s is the hardening law accounting for the isotropic hardening effects; a(r)=a0+βr are material constants and r is an internal variable chosen as the equivalent deviatoric plastic strain defined by:

Referring to the discussion in Section 11, we observe that equation (94) is a particular form of (27), r=0t12e˙p:e˙pdτ of a(r) , 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. Ak(Vk) in equation (28)) for this constitutive model, but we chose a non-associative flow rule instead where,


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 ε˙p=λ˙GσG(σ)=3J2cI1 (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:

where dσij=σijσij0dεije=εijeεij0dp=pp0 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 σij0,εij0,p0 must be the same as the change in the porosity dϵe=dεkke , that is:

By comparing (99) to the second equation of (70), it is evident that the values of Biot coefficients must be: dϵe=dϕe and b= 1 . Taking initial strain 1N=0 , equation (70) thus yields the following constitutive relationships for a linear isotropic poroelastic material:

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, σijσij0=(K23G)ϵeδij+2Gεije(pp0)δijϕϕ0=ϵe is related to the thermodynamic forces as: wf=ρfn(VfVs) . At wfρf=λh(gradp+ρfg) , the fluid is assumed to be in hydraulic equilibrium, implying that: t= 0 . The difference between these two equations yields:


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 ‘ν’):


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


where εij=13ϵδij+eijσij=σδij+sij are the mean and deviatoric strains defined previously; ϵ,eij is the mean stress and σ=σij3 is the deviatoric stress tensor. It is noted that the decomposition into elastic and viscoelastic parts in (102) apply separately to sij=σijσδij , ϵ and the porosity as well such that:

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


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:


where the following relationships are considered for the functions Ψs(ϵe,eije,ϕe,ϵν,eijν)=Ws(ϵe,eije,ϕe)+U(ϵν,eijν) :


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:


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


In the next step, a convex dissipative potential {(σσ0)+(pp0)ξϵν}ϵ˙ν+{(sijsij0)2χeijν}e˙ijν0 is introduced so that based on (112):


which leads to:


where positivity of D(ϵ˙ν,e˙ijν)=12ζϵ˙ν2+ηe˙ijν2 ensures the convexity of ζ0;η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.


where sijsij0=2χeijν+2ηe˙ijν are rheological constants. Note that these equations have been formulated based on the thermodynamics approach while adopting the convex dissipative potential, ξ,ζ,χ,η , in equation (114). Before proceeding further, we will now introduce the Laplace transform, defined for a typical function D(ϵ˙ν,e˙ijν) as follows:

where s is the Laplace transform parameter and f¯(r,s)=L{f(r,t)}=0f(r,t)estdtf¯(r,s)=L1{f¯(r,t)}=12πiΔΔΓiΓ+if¯(r,t)estds . 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 i2= 1 . Taking the Laplace transform of (109), (110), (113) and (116) and solving for the viscous volumetric and deviatoric strains give,

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.


1 - Biot, M.A.: General theory of three-dimensional consolidation. Journal of Applied Physics, 12, 155 164 , 1941
2 - Clausius, Rudolf 1850 On the Motive Power of Heat, and on the Laws which can be deduced from it for the Theory of Heat. Poggendorff’s Annalen der Physik, (Dover Reprint). 0-48659-065-8
3 - O. Coussy, Poromechanics. John Wiley & Sons Ltd.; 2004
4 - N. Dufour, C. J. Leo, F. Deleruyelle, H. Wong, Hydromechanical responses of a decommissioned backfilled tunnel drilled into a poro-viscoelastic medium. Soils and Foundations 4 49 495 507 2009
5 - J. Lemaitre, J. Chaboche, Mech of solid materials, Cambridge University Press, 1990EO
6 - C. J. Leo, M. Kumruzzaman, K. Wong, J. H. Yin, of. E. P. S. Behaviour, in. geofoam, triaxial. true, tests’. compression, Geotextiles, Geomembranes, 2008 26(2), 175 180 .
7 - K. Wong, C. J. A. Leo, elastoplastic. simple, constitutive. hardening, for. E. P. S. model, Geotextiles. geofoam, Geomembranes, 2006 24, 299 310 .
8 - H. Wong, M. Morvan, F. Deleruyelle, C. J. Leo, Analytical study of mine closure behaviour in a poro-elastic medium, Computers and Geotechnics, 2008 35(5), 645 654 .