# Temperature and Pressure Dependence of the Effective Thermal Conductivity of Geomaterials: Numerical Investigation by the Immersed Interface Method

Authors: Duc Phi Do; Dashnor Hoxha
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The present work aims to study the nonlinear effective thermal conductivity of heterogeneous composite-like geomaterials by using a numerical approach based on the immersed interface method (IIM). This method is particularly efficient at solving the diffusion problem in domains containing inner boundaries in the form of perfect or imperfect interfaces between constituents. In this paper, this numerical procedure is extended in the framework of non linear behavior of constituents and interfaces. The performance of the developed tool is then demonstrated through the studies of temperature- and pressure-dependent effective thermal conductivity of geomaterials with imperfect interfaces.


---

## Body

## 1. Introduction

The estimation of the effective transfer properties of heterogeneous media such as geomaterials is still nowadays a challenging research field despite perpetual advances of research and an increasing number of published works. For this class of materials, the overall properties depend not only on the properties of the matrix, the shape and spatial distribution of inclusions, an the volume fractions and properties of constituents, but also on the transfer field’s distribution in the medium and, in many situations, on the interfacial properties between phases.An example of this type of problems is that of the effective thermal conductivity of geomaterials with thermal conductivity of constituents being a function of temperature which lead to an overall nonlinear thermal conductivity. The study of the temperature dependence of thermal conductivity of porous medium like rocks and soils thas been intensively carried out in the last decade due to the large number of applications such as geothermal reservoirs, underground storage of nuclear waste, and petroleum and natural gas geology. Various researches conducted on different types of rocks, soils, and minerals have shown that the thermal conductivity of geomaterials decreases when the temperature increases [1–11] even if in some exceptional cases the conductivity slightly increases with temperature.In addition to the nonlinearity of constituent’s law, the behavior of the interfaces and their state could also play an important role on the nonlinearity of effective properties. Very often these interfaces are favorite places of cracking by debonding that generally leads to a modification of transport properties of the interface.From a numerical point of view, the effect of cracking at the boundary of grains to the overall thermal conductivity can be accounted for by considering the fully or partially debonded inclusions embedded in the homogeneous matrix of materials. This last feature is an important problem not only in gesocience materials, but also in general engineering science materials, and various modeling techniques have been developed to deal with it such as perturbation expansion method [12, 13] or adaptive finite element method [14, 15]. In particular, the immersed interface method (IIM), under such conditions performs quite well. This method, considered as an extension of the finite difference method for the case of media with discontinuities in uniform Cartesian grid points stencil, allows in the case of imperfect interfaces to capture the jump of state and/or flux variables across the interface. Uniform Cartesian grids avoid mesh regeneration and allow for fast flow solvers which contribute to the simplicity and efficiency of the IIM method. Otherwise, the IIM can be used in conjunction with the level set method to treat various problems involving moving interface, non linear interface and free boundary problems (such as Stefan problems and crystal growth, the incompressible flows with moving interface modeled by Navier-Stokes equations to mention a few). Interested readers could refer to [16, 17] and references cited therein for an overview of the advantages and applications of this method.In a previous work of the present authors [18], the IIM was used to take into account the interfacial resistance between linear constituent phases. In this context, the IIM allows calculating the effective transfer properties and enables demonstrating the role of various factors such as shape, size, spatial distribution, and volume fraction of constituents as well as their properties of constituents and those of interface on overall effective properties.In the present paper, this numerical approach is extended to the evaluation of effective thermal conductivity of geomaterials considering the nonlinearity of constituents and the presence of imperfect interfaces. The paper is organized as follows. Firstly, we describe the principle of IIM used to solve the non linear heat transfer problem with the presence of contact resistance. Then, the application of the numerical procedure to evaluate the effective thermal conductivity of non linear composite-like geomaterials is studied by counting for the interfacial resistance between phases. In the whole paper, a lower underlined symbol indicates a vector while a bold one represents a second-order tensor.


## 2. Solving Nonlinear Heat Transfer Problem by the Immersed Interface Method

### 2.1. Mathematical Model

The problem considered here is that of a heterogeneous material constituted by a matrix containing inclusions of smooth shapes whose properties are different from that of the matrix. Moreover, the interfaces between the matrix and any inclusion (Figure1) are considered to be thin layers (no thickness) having their own properties and sufficiently smooth to assure the existence of all derivatives involved in equations developed later. The interfaces are not necessarily perfect so that a jump in some variables (state or/and flux variables) could be observed across them. For the sake of simplicity, the presentation of the method is limited here on 2D steady state heat transfer problem without source. The extension of the method for more general cases presents no conceptual difficulties. Local behavior of each constituent phase of heterogeneous media is characterized by the conservation and Fourier law written as
(1)div(q_)=0,q_=-K(u)·grad(u)_,
where q_ and K stand for the flux and the second order tensor of thermal conductivity. As the behavior of constituents is assumed non linear, K is supposed to be a known function of the temperature u.Figure 1
Schematic presentation of the problem of interface in a uniform Cartesian grid for the interface problem.Without losing the generality, in what follows, the isotropic behavior will be considered; that is,K(u)=K(u)·1 with 1 being unity second-order tensor. The combination of these two equations leads to the following non linear partial differential equation with respect to the unknown u:
(2)div(K(u)·grad(u)_)=0.
The technique of resolution of this type of elliptical equation by IIM method was first used by [19] in the context of magnetorheological fluids problem of perfect interface (i.e., assuming the continuity of flux q and of the solution u across the interface). In the extended mathematical framework adopted here, a jump of solution [u] on the contact between matrix and an inclusion is considered:
(3)[u]=u+-u-,[K(u)∂u∂n]=K+(u)un+-K-(u)un-=0
with n_ being the outward oriented unitary normal vector to the interface.From a physical point of view, these conditions describe an imperfect contact between phases (e.g., due to the presence of the roughness or air film at the interface). This imperfection is often stated in thermal or hydraulic transport problems even if it is very difficult to be measured experimentally.A common hypothesis used by many authors (see, e.g., [12, 13, 18] and references cited therein) is to consider the jump of the solution u across the interface being proportional to the normal component of flux across the interface, with the coefficient of proportionality being a characteristic interface parameter, α, called the resistance of interface:
(4)[u]=u+-u-=-α(q_·n_)=αK+(u)un+=αK-(u)un-.
It easy to observe that the perfect interface problem is a special case of the problem with resistance of interface equal to zero (α=0).Hereafter, the superscripts (+) and (−) present values evaluated at two opposed sides of a given interface (Figure1).


### 2.2. Iterative Resolution of Nonlinear Elliptical Equation by Immersed Interface Method

Like finite difference method from which it has been derived, the IIM uses a uniform Cartesian grid which does not need to coincide with the interfaces or boundaries of internal domains to obtain discretized form of equation:(5)∂∂x(K∂ui,j∂x)+∂∂y(K∂ui,j∂y)=∑m=1ns‍γmui+im,j+jm+Ci,j=0.
This equation is written at each node of the grid with coordinates (xi, yj), identified in the 2D case by the couple of indexes (i, j). Since the behavior of constituents is supposed non linear, then K=K(ui,j). To resolve the system of (2) obtained by the discretization of differential equation, we use the so-called substitution method, similar to that proposed by [19]. Following this scheme, from an initial guess of solution u0, the (k+1)th iterative step consists to obtain a new estimation of solution uk+1 of partial derivative equation employing the known information of the previous step uk.Bearing in mind spatial and temporal discretization, the solutionu and the thermal conductivity K must be written as ui,j=uk+1(xi,yj) and K=K(ui,jk). However, in order to keep simple formulas, this dependency is omitted from notation but it is implicitly assumed for the rest of the paper. Other parameters ns and γm represent, respectively, the number of grid points and the weighting coefficient of mth node involved in the evaluation of the solution at point (xi, yj), and indexes im and jm take values in the set {0,±1,±2,…}. Since γm coefficients, indexes im and jm, and the correction term Ci,j depend on point (i, j), it needs to be written γijm. Again for the sake of simplicity, this dependency is omitted from notation and is assumed implicitly. Furthermore, a constant number of stencils is kept for all nodes of the same type but different for regular and irregular points.Aregular point is a grid node away from the interface (see Figure 1) for which the centered FDM differentiation scheme of differential equation to be solved is used. For such nodes, the stencil contains five points so that ns=5. All nodes where the approximation of solution involves only points on the same side of the interface are regular points and the approximation of solution at these points coincides with classical formula of a standard five points stencil. For these points the approximation of the partially derivative equation to the second-order accuracy demonstrates a vanishing correction term (Ci,j=0) and is given as follows:
(6)∂∂x(K∂ui,j∂x)+∂∂y(K∂ui,j∂y)=Ki,j-1/2h2ui,j-1+Ki-1/2,jh2ui-1,j+Ki+1/2,jh2ui+1,j+Ki,j+1/2h2ui,j+1-Ki,j-1/2+Ki-1/2,j+Ki+1/2,j+Ki,j+1/2h2ui,j+O(h3),
where
(7)Ki,j=K(uk(xi,yj)),ui,j=uk+1(xi,yj),Ki±1/2,j=12(Ki,j+Ki±1,j),Ki,j±1/2=12(Ki,j+Ki,j±1).
If the grid point is near the interface, the points involved in the approximation of the solution would be from materials in both sides of interface. These points are called irregulars since the use of the standard five points stencil does not insure any more the second-order accuracy of the solution. As discussed in a [18, 20] and [16, 17, 19] a good choice in order to maintain the second-order accuracy at irregular points is to take a stencil of nine points, that is, for these irregular nodes ns=9.The procedure of constructing a second order accurate solution at irregular points is detailed in the previously mentioned references and is briefly described here. Let (x*, y*) be the projection on the interface of an irregular node (xi, yj). A local coordinates system attached to the projected point is assumed with two axis (ξ, η) coinciding with normal and tangential directions of the interface at this point (Figure 1):
(8)ξ=(x-x*)cosθ+(y-y*)sinθ,η=-(x-x*)sinθ+(y-y*)cosθ,
with θ being the angle between x and ξ axes.The relationship between the global and local coordinates systems gives(9)uξ=uxcosθ+uysinθ;uη=-uxsinθ+uycosθ;Kξ=Kxcosθ+Kysinθ;Kη=-Kxsinθ+Kycosθ,
where ux and Kx indicate, respectively, the derivative of u and K with respect to x.It is then possible to replace the original discretized equation (5) by one written in the local coordinate of (ξ, η):
(10)∂∂x(K∂ui,j∂x)+∂∂y(K∂ui,j∂y)=∂∂ξ(K∂ui,j∂ξ)+∂∂η(K∂ui,j∂η)=∑m=1ns‍γmui+im,j+jm+Ci,j.
Further, in respect of local coordinates system, the solution ui+im,j+jm at a grid point could be written as a Taylor expansion at (x*, y*) so that
(11)u(xi+im,yj+jm)=u(ξm,ηm)=u±+ξmuξ±+ηmuη±+12ξm2uξξ±+12ηm2uηη±+ξmηmuξη±+12ξmηm2uξηη±+O(h3).
The reasons of keeping third-order derivative of solution uξηη in Taylor’s development (11) are related to their presence on the interface conditions relations as it will be clear later (see (15)). Additionally, the values u±, uξ±, uη±, uξξ±, uηη±, uξη±, and uξηη± in (11) represent the solution and its derivatives determined at projection point (x*, y*) of the interface. The (+) or (−) sign is chosen depending on whether (ξm,ηm) belongs to the (+) or (−) side of the interface.By replacing in (10) each component ui+im,j+jm as given by (11) and rearranging terms, the following expression of the truncation error at an irregular point is obtained (i, j):
(12)κi,j=a1u-+a2u++a3uξ-+a4uξ++a5uη-+a6uη++a7uξξ-+a8uξξ++a9uηη-+a10uηη++a11uξη-+a12uξη++a13uξηη-+a14uξηη+-Ci,j+O(maxm|γm|h3).The coefficients{a1,a2,…,a14} of (12) depend on the position of the stencil point relative to the interface:
(13)a1=∑(m∈M-)‍γm;a2=∑(m∈M+)‍γm;a3=∑(m∈M-)‍γmξm;a4=∑(m∈M+)‍γmξm;a5=∑(m∈M-)‍γmηm;a6=∑(m∈M+)‍γmηm;a7=12∑(m∈M-)‍γmξm2;a8=12∑(m∈M+)‍γmξm2;a9=12∑(m∈M-)‍γmηm2;a10=12∑(m∈M+)‍γmηm2;a11=∑(m∈M-)‍γmξmηm;a12=∑(m∈M+)‍γmξmηm;a13=12∑(m∈M-)‍γmξmηm2;a14=12∑(m∈M+)‍γmξmηm2,
where the sets M± are defined as
(14)M±={m:(ξm,ηm)∈Ω±}.
At this end it is important to specify the relations of the solution and its derivatives at both sides of the interface. After some algebraic manipulation, from the contact conditions (3) the following relations are deduced:
(15)u+=u-+αK-uξ-,uξ+=ρuξ-,uη+=(1-αχ′′K-)uη-+αKη-uξ-+αK-uξη-,uξξ+=ρuξξ--[Kη+uη+-Kη-uη-+Kξ+uξ+-Kξ-uξ-+K+uηη+-K-uηη-]K+;uηη+=uηη--χ′′(uξ+-uξ-)+α(χ′′Kξ-+Kηη-)uξ-+αKη-(uξη--χ′′uη-)+αKη-(uξη--χ′′uη-)+αK-(χ′′uξξ-+uξηη-2χ′′uηη-)-α(χ′′)2K-uξ-,uξη+=(Kη-uξ--Kη+uξ+)K++ρ(uξη--χ′′uη-)+χ′′uη+,uξηη+=ρuξηη-+χ′′(ρuξξ--uξξ+)+[uξ-(χ′′Kξ-+Kηη-)-uξ+(χ′′Kξ++Kηη+)]K++2[Kη-uξη--Kη+uξη+-χ′′(Kη-uη--Kη+uη+)]K+-2χ′′(ρuηη--uηη+),
with ρ=K-/K+. The curvature χ′′ of the interface at (x*, y*) is the second order derivative of function ξ=χ(η) with respect to η. We note also that at (x*, y*) we have χ(0)=0 and for a smooth interface (as assumed here) χ′(0)=0.The interface relations (15) are used to recast (12) in a compact form using only quantities from one side of the interface (the Ω- side, e.g.):
(16)κi,j=B1u-+B2uξ-+B3uη-+B4uξξ-+B5uηη-+B6uξη-+B7uξηη-+O(maxm|γm|h3),
where Bi (i=1,7) are expressions of coefficients {a1,a2,…,a14} and other quantities used in (15). Since we are looking for a second order accuracy, the coefficients Bi of κi,j expression must vanish in the same time with the correction term Ci,j. These constraints lead to following system of equations with respect to variables ai:
(17)a1+a2=0,a2αK-+a3+a4ρ+a6αKη-+a8c8(2)+a10c10(2)+a12c12(2)+a14c14(2)=Kξ-,a5+a6(1-αχ′′K-)+a8c8(3)-2αχ′′Kη-a10+a14c14(3)+χ′′(1-ρ-αχ′′K-)a12=Kη-,a7+(ρ-αχ′′K-)a8+αχ′′K-a10+3α(χ′′)2K-a14=K-,a9+(1-2αχ′′K-)a10+(ρ-1+2αχ′′K-)a8+3χ′′(1-ρ-2αχ′′K-)a14=K-,αK-a6-(2αKη--αρKη+)a8+2αKη-a10+a11+(ρ+αχ′′K-)a12+a14c14(6)=0,-αK-a8+αK-a10+a13+(ρ+3αχ′′K-)a14=0,
where
(18)c8(2)=-(αKη-Kη++αK+Kηη--Kξ-+ρKξ+)/K+-(1+αKξ--ρ)χ′′+αK-(χ′′)2,c10(2)=χ′′(1-ρ)+αKηη-+αχ′′(Kξ--χ′′K-),c12(2)=(Kη--ρKη+)K++αχ′′Kη-,c14(2)=(-2Kη-Kη++2ρ(Kη+)2)(K+)2+(αχ′′Kη-Kη++Kηη--ρKηη+)K++3χ′′(αKηη-+χ′′+αχ′′Kξ--ρχ′′(1+αχ′′K+)),c8(3)=(Kη--Kη+)K++αχ′′(2Kη-+ρKη+),c14(3)=-χ′′(3Kη--Kη+-2ρKη+)K+-(χ′′)2(6αKη-+αρKη+),c14(6)=2(Kη--ρKη+)K++αχ′′(6Kη-+ρKη+).
Since all coefficients ai are defined as functions of γm (see (13)), one could write the system (17) in a compact form as
(19)[A]·γ_=b_
with the vector γ_=[γ1,γ2,…,γns]t and b_=[0,Kξ-,Kη-,K-,K-,0,0]t.Thus, the central point to obtain the solution of the problem (2) by using IIM is to construct and resolve the system (19). Note that, at irregular points, the set of (19) could be obtained using at least six points stencil as proposed in the original IIM method [21]. However, as discussed in some recent works [16, 17, 19], nine-point stencil (ns=9) is preferred because the resulting linear system of equations is in this case a block tridiagonal one.However, using of a nine points stencil (ns=9) leads to an underdetermined system (19) whose solution could be obtained by using the maximum principle [16, 17, 19] and an optimization approach. For that, the following constrained quadratic optimization problem is considered:
(20)minγ{12∥γ-g∥22},
with
(21)γm≥0if(im,jm)≠(0,0),γm<0if(im,jm)=(0,0).
In this equation, the vector g_=[g1,g2,…,gns]t has the following components:
(22)gm=Ki+im/2,j+jm/2h2,if(im,jm)∈{(-1,0),(1,0),(0,-1),(0,1)},gm=-Ki,j-1/2+Ki-1/2,j+Ki+1/2,j+Ki,j+1/2h2,if(im,jm)=(0,0),gm=0otherwise.
This vector g_ represents in effect the coefficients of the standard five points stencil to whom the solution of the optimization problem has to coincide in case of continuous properties of constituent phases (i.e., when [K]=0).The existence of the solution of the optimization problem (20) as well as the convergence of the substitution method in that case has been demonstrated in [16, 17] and will not be discussed here. We just mention that, as in a previous work of these authors [18], this optimization problem is solved by using the optimization toolbox available in MATLAB.


### 2.3. Computing Solution, Its Derivatives at Interface, and Effective Properties

As previously outlined, the non linear transfer problem is treated here by using an iterative procedure consisting for each iteration to solve a linearized partial differential equation. The resolution of this later by using the IIM method needs to determine the set of weighting coefficientsγm. At each irregular point, these coefficients have to be evaluated by solving an undetermined system via an optimization approach. The coefficients K±, Kξ±, and Kη± of this system of (17) are in reality defined as functions of the solution u(x*,y*) at interface. For the substitution scheme followed here, these coefficients are calculated at each step from the solution of the previous iteration uk at projection point (x*, y*).The solution value at this pointu(x*,y*) is interpolated from solution at grid points using a least squares scheme. For this interpolation, grid points situated at the same side of the considered point (xi,yj) are used. As discussed in chapter 6 of [17] as well as in the contribution of [19], using this approach by involving in least squares approximation at least six points of the same side avoids to use the interface relations and maintains the second order accuracy. For example, if the irregular point (xi,yj) belongs to Ω-, one can interpolate the solution at its projection point (x*,y*) from the values of points in the same side Ω- as follows:
(23)u-(x*,y*)=∑(xm,yn)∈Ω-,|(xm,yn)-(x*,y*)|≤Rε‍γm,num,n,
where Rε can be chosen between 4.1h and 6.1h.The same method of least squares scheme is also used to calculate the derivatives of the solutionux and uy using uξ by applying relation (9). Note that once the value of the solution u- is known (and so are the values of K-=K(u-)) as well as its derivative uξ- at Ω- side of the interface, it is possible to evaluate the respective values at the other side Ω+ of interface by using the following interface relation:
(24)u+=u-+αK-uξ-;K+=K(u+).The same approach can be used to determine the coefficientsKx- and Ky- at (x*, y*) by involving values of Km,n=K(um,n) determined at grid points situated at the same side Ω- (or Ω+). We obtain finally Kξ± and Kη± from the relation written in (9).Once all coefficients of the finite difference scheme are known, it is possible to solve the system of (5) and obtain the solution (uk+1) of the (k+1)th step at all grid nodes. The iterative calculation is finished when the convergence criterion is reached. For all numerical simulations presented in the next sections, the norm of residue of the iterative scheme is fixed to 10-6. These results numerically confirm the convergence of the iterative procedure which is obtained after only a small number of iteration (~5 iterations).From the solution of the transfer problem, one could use it to calculate the overall conductivity of the heterogeneous media from the well-known relation of volume average quantities:(25)〈q_〉=-Keff·{〈grad(u)_〉+1|Ω|∫Γ‍α[q_(s)·n_(s)]n_(s)ds},
where 〈·〉 denotes the average over all elements of the representative volume.


## 2.1. Mathematical Model

The problem considered here is that of a heterogeneous material constituted by a matrix containing inclusions of smooth shapes whose properties are different from that of the matrix. Moreover, the interfaces between the matrix and any inclusion (Figure1) are considered to be thin layers (no thickness) having their own properties and sufficiently smooth to assure the existence of all derivatives involved in equations developed later. The interfaces are not necessarily perfect so that a jump in some variables (state or/and flux variables) could be observed across them. For the sake of simplicity, the presentation of the method is limited here on 2D steady state heat transfer problem without source. The extension of the method for more general cases presents no conceptual difficulties. Local behavior of each constituent phase of heterogeneous media is characterized by the conservation and Fourier law written as
(1)div(q_)=0,q_=-K(u)·grad(u)_,
where q_ and K stand for the flux and the second order tensor of thermal conductivity. As the behavior of constituents is assumed non linear, K is supposed to be a known function of the temperature u.Figure 1
Schematic presentation of the problem of interface in a uniform Cartesian grid for the interface problem.Without losing the generality, in what follows, the isotropic behavior will be considered; that is,K(u)=K(u)·1 with 1 being unity second-order tensor. The combination of these two equations leads to the following non linear partial differential equation with respect to the unknown u:
(2)div(K(u)·grad(u)_)=0.
The technique of resolution of this type of elliptical equation by IIM method was first used by [19] in the context of magnetorheological fluids problem of perfect interface (i.e., assuming the continuity of flux q and of the solution u across the interface). In the extended mathematical framework adopted here, a jump of solution [u] on the contact between matrix and an inclusion is considered:
(3)[u]=u+-u-,[K(u)∂u∂n]=K+(u)un+-K-(u)un-=0
with n_ being the outward oriented unitary normal vector to the interface.From a physical point of view, these conditions describe an imperfect contact between phases (e.g., due to the presence of the roughness or air film at the interface). This imperfection is often stated in thermal or hydraulic transport problems even if it is very difficult to be measured experimentally.A common hypothesis used by many authors (see, e.g., [12, 13, 18] and references cited therein) is to consider the jump of the solution u across the interface being proportional to the normal component of flux across the interface, with the coefficient of proportionality being a characteristic interface parameter, α, called the resistance of interface:
(4)[u]=u+-u-=-α(q_·n_)=αK+(u)un+=αK-(u)un-.
It easy to observe that the perfect interface problem is a special case of the problem with resistance of interface equal to zero (α=0).Hereafter, the superscripts (+) and (−) present values evaluated at two opposed sides of a given interface (Figure1).


## 2.2. Iterative Resolution of Nonlinear Elliptical Equation by Immersed Interface Method

Like finite difference method from which it has been derived, the IIM uses a uniform Cartesian grid which does not need to coincide with the interfaces or boundaries of internal domains to obtain discretized form of equation:(5)∂∂x(K∂ui,j∂x)+∂∂y(K∂ui,j∂y)=∑m=1ns‍γmui+im,j+jm+Ci,j=0.
This equation is written at each node of the grid with coordinates (xi, yj), identified in the 2D case by the couple of indexes (i, j). Since the behavior of constituents is supposed non linear, then K=K(ui,j). To resolve the system of (2) obtained by the discretization of differential equation, we use the so-called substitution method, similar to that proposed by [19]. Following this scheme, from an initial guess of solution u0, the (k+1)th iterative step consists to obtain a new estimation of solution uk+1 of partial derivative equation employing the known information of the previous step uk.Bearing in mind spatial and temporal discretization, the solutionu and the thermal conductivity K must be written as ui,j=uk+1(xi,yj) and K=K(ui,jk). However, in order to keep simple formulas, this dependency is omitted from notation but it is implicitly assumed for the rest of the paper. Other parameters ns and γm represent, respectively, the number of grid points and the weighting coefficient of mth node involved in the evaluation of the solution at point (xi, yj), and indexes im and jm take values in the set {0,±1,±2,…}. Since γm coefficients, indexes im and jm, and the correction term Ci,j depend on point (i, j), it needs to be written γijm. Again for the sake of simplicity, this dependency is omitted from notation and is assumed implicitly. Furthermore, a constant number of stencils is kept for all nodes of the same type but different for regular and irregular points.Aregular point is a grid node away from the interface (see Figure 1) for which the centered FDM differentiation scheme of differential equation to be solved is used. For such nodes, the stencil contains five points so that ns=5. All nodes where the approximation of solution involves only points on the same side of the interface are regular points and the approximation of solution at these points coincides with classical formula of a standard five points stencil. For these points the approximation of the partially derivative equation to the second-order accuracy demonstrates a vanishing correction term (Ci,j=0) and is given as follows:
(6)∂∂x(K∂ui,j∂x)+∂∂y(K∂ui,j∂y)=Ki,j-1/2h2ui,j-1+Ki-1/2,jh2ui-1,j+Ki+1/2,jh2ui+1,j+Ki,j+1/2h2ui,j+1-Ki,j-1/2+Ki-1/2,j+Ki+1/2,j+Ki,j+1/2h2ui,j+O(h3),
where
(7)Ki,j=K(uk(xi,yj)),ui,j=uk+1(xi,yj),Ki±1/2,j=12(Ki,j+Ki±1,j),Ki,j±1/2=12(Ki,j+Ki,j±1).
If the grid point is near the interface, the points involved in the approximation of the solution would be from materials in both sides of interface. These points are called irregulars since the use of the standard five points stencil does not insure any more the second-order accuracy of the solution. As discussed in a [18, 20] and [16, 17, 19] a good choice in order to maintain the second-order accuracy at irregular points is to take a stencil of nine points, that is, for these irregular nodes ns=9.The procedure of constructing a second order accurate solution at irregular points is detailed in the previously mentioned references and is briefly described here. Let (x*, y*) be the projection on the interface of an irregular node (xi, yj). A local coordinates system attached to the projected point is assumed with two axis (ξ, η) coinciding with normal and tangential directions of the interface at this point (Figure 1):
(8)ξ=(x-x*)cosθ+(y-y*)sinθ,η=-(x-x*)sinθ+(y-y*)cosθ,
with θ being the angle between x and ξ axes.The relationship between the global and local coordinates systems gives(9)uξ=uxcosθ+uysinθ;uη=-uxsinθ+uycosθ;Kξ=Kxcosθ+Kysinθ;Kη=-Kxsinθ+Kycosθ,
where ux and Kx indicate, respectively, the derivative of u and K with respect to x.It is then possible to replace the original discretized equation (5) by one written in the local coordinate of (ξ, η):
(10)∂∂x(K∂ui,j∂x)+∂∂y(K∂ui,j∂y)=∂∂ξ(K∂ui,j∂ξ)+∂∂η(K∂ui,j∂η)=∑m=1ns‍γmui+im,j+jm+Ci,j.
Further, in respect of local coordinates system, the solution ui+im,j+jm at a grid point could be written as a Taylor expansion at (x*, y*) so that
(11)u(xi+im,yj+jm)=u(ξm,ηm)=u±+ξmuξ±+ηmuη±+12ξm2uξξ±+12ηm2uηη±+ξmηmuξη±+12ξmηm2uξηη±+O(h3).
The reasons of keeping third-order derivative of solution uξηη in Taylor’s development (11) are related to their presence on the interface conditions relations as it will be clear later (see (15)). Additionally, the values u±, uξ±, uη±, uξξ±, uηη±, uξη±, and uξηη± in (11) represent the solution and its derivatives determined at projection point (x*, y*) of the interface. The (+) or (−) sign is chosen depending on whether (ξm,ηm) belongs to the (+) or (−) side of the interface.By replacing in (10) each component ui+im,j+jm as given by (11) and rearranging terms, the following expression of the truncation error at an irregular point is obtained (i, j):
(12)κi,j=a1u-+a2u++a3uξ-+a4uξ++a5uη-+a6uη++a7uξξ-+a8uξξ++a9uηη-+a10uηη++a11uξη-+a12uξη++a13uξηη-+a14uξηη+-Ci,j+O(maxm|γm|h3).The coefficients{a1,a2,…,a14} of (12) depend on the position of the stencil point relative to the interface:
(13)a1=∑(m∈M-)‍γm;a2=∑(m∈M+)‍γm;a3=∑(m∈M-)‍γmξm;a4=∑(m∈M+)‍γmξm;a5=∑(m∈M-)‍γmηm;a6=∑(m∈M+)‍γmηm;a7=12∑(m∈M-)‍γmξm2;a8=12∑(m∈M+)‍γmξm2;a9=12∑(m∈M-)‍γmηm2;a10=12∑(m∈M+)‍γmηm2;a11=∑(m∈M-)‍γmξmηm;a12=∑(m∈M+)‍γmξmηm;a13=12∑(m∈M-)‍γmξmηm2;a14=12∑(m∈M+)‍γmξmηm2,
where the sets M± are defined as
(14)M±={m:(ξm,ηm)∈Ω±}.
At this end it is important to specify the relations of the solution and its derivatives at both sides of the interface. After some algebraic manipulation, from the contact conditions (3) the following relations are deduced:
(15)u+=u-+αK-uξ-,uξ+=ρuξ-,uη+=(1-αχ′′K-)uη-+αKη-uξ-+αK-uξη-,uξξ+=ρuξξ--[Kη+uη+-Kη-uη-+Kξ+uξ+-Kξ-uξ-+K+uηη+-K-uηη-]K+;uηη+=uηη--χ′′(uξ+-uξ-)+α(χ′′Kξ-+Kηη-)uξ-+αKη-(uξη--χ′′uη-)+αKη-(uξη--χ′′uη-)+αK-(χ′′uξξ-+uξηη-2χ′′uηη-)-α(χ′′)2K-uξ-,uξη+=(Kη-uξ--Kη+uξ+)K++ρ(uξη--χ′′uη-)+χ′′uη+,uξηη+=ρuξηη-+χ′′(ρuξξ--uξξ+)+[uξ-(χ′′Kξ-+Kηη-)-uξ+(χ′′Kξ++Kηη+)]K++2[Kη-uξη--Kη+uξη+-χ′′(Kη-uη--Kη+uη+)]K+-2χ′′(ρuηη--uηη+),
with ρ=K-/K+. The curvature χ′′ of the interface at (x*, y*) is the second order derivative of function ξ=χ(η) with respect to η. We note also that at (x*, y*) we have χ(0)=0 and for a smooth interface (as assumed here) χ′(0)=0.The interface relations (15) are used to recast (12) in a compact form using only quantities from one side of the interface (the Ω- side, e.g.):
(16)κi,j=B1u-+B2uξ-+B3uη-+B4uξξ-+B5uηη-+B6uξη-+B7uξηη-+O(maxm|γm|h3),
where Bi (i=1,7) are expressions of coefficients {a1,a2,…,a14} and other quantities used in (15). Since we are looking for a second order accuracy, the coefficients Bi of κi,j expression must vanish in the same time with the correction term Ci,j. These constraints lead to following system of equations with respect to variables ai:
(17)a1+a2=0,a2αK-+a3+a4ρ+a6αKη-+a8c8(2)+a10c10(2)+a12c12(2)+a14c14(2)=Kξ-,a5+a6(1-αχ′′K-)+a8c8(3)-2αχ′′Kη-a10+a14c14(3)+χ′′(1-ρ-αχ′′K-)a12=Kη-,a7+(ρ-αχ′′K-)a8+αχ′′K-a10+3α(χ′′)2K-a14=K-,a9+(1-2αχ′′K-)a10+(ρ-1+2αχ′′K-)a8+3χ′′(1-ρ-2αχ′′K-)a14=K-,αK-a6-(2αKη--αρKη+)a8+2αKη-a10+a11+(ρ+αχ′′K-)a12+a14c14(6)=0,-αK-a8+αK-a10+a13+(ρ+3αχ′′K-)a14=0,
where
(18)c8(2)=-(αKη-Kη++αK+Kηη--Kξ-+ρKξ+)/K+-(1+αKξ--ρ)χ′′+αK-(χ′′)2,c10(2)=χ′′(1-ρ)+αKηη-+αχ′′(Kξ--χ′′K-),c12(2)=(Kη--ρKη+)K++αχ′′Kη-,c14(2)=(-2Kη-Kη++2ρ(Kη+)2)(K+)2+(αχ′′Kη-Kη++Kηη--ρKηη+)K++3χ′′(αKηη-+χ′′+αχ′′Kξ--ρχ′′(1+αχ′′K+)),c8(3)=(Kη--Kη+)K++αχ′′(2Kη-+ρKη+),c14(3)=-χ′′(3Kη--Kη+-2ρKη+)K+-(χ′′)2(6αKη-+αρKη+),c14(6)=2(Kη--ρKη+)K++αχ′′(6Kη-+ρKη+).
Since all coefficients ai are defined as functions of γm (see (13)), one could write the system (17) in a compact form as
(19)[A]·γ_=b_
with the vector γ_=[γ1,γ2,…,γns]t and b_=[0,Kξ-,Kη-,K-,K-,0,0]t.Thus, the central point to obtain the solution of the problem (2) by using IIM is to construct and resolve the system (19). Note that, at irregular points, the set of (19) could be obtained using at least six points stencil as proposed in the original IIM method [21]. However, as discussed in some recent works [16, 17, 19], nine-point stencil (ns=9) is preferred because the resulting linear system of equations is in this case a block tridiagonal one.However, using of a nine points stencil (ns=9) leads to an underdetermined system (19) whose solution could be obtained by using the maximum principle [16, 17, 19] and an optimization approach. For that, the following constrained quadratic optimization problem is considered:
(20)minγ{12∥γ-g∥22},
with
(21)γm≥0if(im,jm)≠(0,0),γm<0if(im,jm)=(0,0).
In this equation, the vector g_=[g1,g2,…,gns]t has the following components:
(22)gm=Ki+im/2,j+jm/2h2,if(im,jm)∈{(-1,0),(1,0),(0,-1),(0,1)},gm=-Ki,j-1/2+Ki-1/2,j+Ki+1/2,j+Ki,j+1/2h2,if(im,jm)=(0,0),gm=0otherwise.
This vector g_ represents in effect the coefficients of the standard five points stencil to whom the solution of the optimization problem has to coincide in case of continuous properties of constituent phases (i.e., when [K]=0).The existence of the solution of the optimization problem (20) as well as the convergence of the substitution method in that case has been demonstrated in [16, 17] and will not be discussed here. We just mention that, as in a previous work of these authors [18], this optimization problem is solved by using the optimization toolbox available in MATLAB.


## 2.3. Computing Solution, Its Derivatives at Interface, and Effective Properties

As previously outlined, the non linear transfer problem is treated here by using an iterative procedure consisting for each iteration to solve a linearized partial differential equation. The resolution of this later by using the IIM method needs to determine the set of weighting coefficientsγm. At each irregular point, these coefficients have to be evaluated by solving an undetermined system via an optimization approach. The coefficients K±, Kξ±, and Kη± of this system of (17) are in reality defined as functions of the solution u(x*,y*) at interface. For the substitution scheme followed here, these coefficients are calculated at each step from the solution of the previous iteration uk at projection point (x*, y*).The solution value at this pointu(x*,y*) is interpolated from solution at grid points using a least squares scheme. For this interpolation, grid points situated at the same side of the considered point (xi,yj) are used. As discussed in chapter 6 of [17] as well as in the contribution of [19], using this approach by involving in least squares approximation at least six points of the same side avoids to use the interface relations and maintains the second order accuracy. For example, if the irregular point (xi,yj) belongs to Ω-, one can interpolate the solution at its projection point (x*,y*) from the values of points in the same side Ω- as follows:
(23)u-(x*,y*)=∑(xm,yn)∈Ω-,|(xm,yn)-(x*,y*)|≤Rε‍γm,num,n,
where Rε can be chosen between 4.1h and 6.1h.The same method of least squares scheme is also used to calculate the derivatives of the solutionux and uy using uξ by applying relation (9). Note that once the value of the solution u- is known (and so are the values of K-=K(u-)) as well as its derivative uξ- at Ω- side of the interface, it is possible to evaluate the respective values at the other side Ω+ of interface by using the following interface relation:
(24)u+=u-+αK-uξ-;K+=K(u+).The same approach can be used to determine the coefficientsKx- and Ky- at (x*, y*) by involving values of Km,n=K(um,n) determined at grid points situated at the same side Ω- (or Ω+). We obtain finally Kξ± and Kη± from the relation written in (9).Once all coefficients of the finite difference scheme are known, it is possible to solve the system of (5) and obtain the solution (uk+1) of the (k+1)th step at all grid nodes. The iterative calculation is finished when the convergence criterion is reached. For all numerical simulations presented in the next sections, the norm of residue of the iterative scheme is fixed to 10-6. These results numerically confirm the convergence of the iterative procedure which is obtained after only a small number of iteration (~5 iterations).From the solution of the transfer problem, one could use it to calculate the overall conductivity of the heterogeneous media from the well-known relation of volume average quantities:(25)〈q_〉=-Keff·{〈grad(u)_〉+1|Ω|∫Γ‍α[q_(s)·n_(s)]n_(s)ds},
where 〈·〉 denotes the average over all elements of the representative volume.


## 3. Variation of the Effective Thermal Conductivity of Geomaterials with respect to Temperature and Pressure

### 3.1. Temperature Dependence

Beginning from the pioneer work of [3], the temperature-dependent thermal properties of geomaterials have been studied for a long time in relation with various engineering applications such as geothermal reservoirs, underground storage of nuclear waste or petroleum and natural gas geology. This subject is always in continuous progress with a remarkable number of results conducted on different types of soils and rocks in the last decade [2, 4, 6–9]. A general tendency of decreasing of thermal conductivity as a function of temperature can be stated from these contributions. At the same time, a thermal damage is reported as a consequence of the contrast of constituent properties of geomaterials when the temperature increases. This type of damage known as thermal cracking developed principally in contacts of constituent phases of the heterogeneous media like geomaterials.In this part, we use the IIM method to solve the problem of effective thermal conductivity estimation of geomaterials whose constituent thermal conductivities depend on temperature. In addition, the effects of cracking interfaces with temperature on the effective thermal conductivity are considered using imperfect interface. In that case, the nonlinearity of effective thermal property reflects not only the nonlinearity of the intrinsic properties of matrix and inclusions but also the properties of interfaces, which when debonding produce an extra thermal resistance. In order to describe the progressive cracking of interfaces when the temperature evolves, the properties of interface are supposed to vary with coordinates following a known function. In the simplest case a piecewise interfacial resistance function that takes some values for debonding parts of interfaces and some other values for intact ones could be used. As an example, the debonding parts of the interface of a circular inclusion could be given in function of an angleϕ measured from an arbitrary chosen direction (Figure 2). In extreme cases, full interface could be considered as debonded. One could anticipate that, because of the temperature dependent conductivity of local constituents and interface debonding, the overall conductivity would be highly temperature dependent [14, 15, 22–26]. As demonstrated in a previous contribution from the authors [18], the contact resistance leads also to a size-dependent behavior of heterogeneous materials as observed in reality: for the same volumetric fraction of inclusions, a more pronounced effect of the interface is observed for smaller inclusions. If not otherwise stated in following simulations, the size of inclusion is equal to R=1μm.(a) Single partial debonding of circular inclusion. (b) Double partial debonding of circular inclusion.
(a)
(b)For the sake of clarity, the geometry as well as the initial and boundary conditions of the considered problem is depicted in Figure2. More precisely, a constant initial temperature u is considered in whole sample, and the mixed boundary conditions consist in applying the temperatures u and u+Δu at the two opposing faces (to the Oy direction in the examples presented here) of the rectangular domain, whereas the zero heat flux is supposed to be to two other faces in the orthogonal direction (faces normal to the Ox direction) to the first. In what follows, a constant 320×320 grid is used for all calculations. This value obtained from the analysis of the sensitivity of numerical results to the grid density in previous work [18] gives excellent agreement with the analytical solutions in the context of the linear problem.For illustration purposes, we consider a sample of Callovo-Oxfordian argillite (see, e.g., [18] for a brief description) constituted by a clay matrix with quartz inclusions at the mesoscopic level. The constant thermal conductivity of the matrix is taken as K1(u)=1.8(W·m-1·K-1). Otherwise the thermal conductivity of quartz significantly decreases with temperature and following evolution is reported in published works [1, 3]: K2(u)=7.7/(0.0045×u-0.3863)(W·m-1·K-1). It must emphasize however that despite simple expressions of thermal conductivities used for this material, one could use any expression provided by experimental considerations. Concerning the thermal resistance of interface in debonded parts, an arbitrary value is taken equal to 10-6(m2·K·W-1), while perfect interface conditions (α=0) are supposed for bonded parts.For simulations, firstly the sample of heterogeneous material with a single circular inclusions and partially (simple or double) debonded interface is considered. As the illustration in Figure3, the temperature fields are presented on the material containing a single circular inclusion with, respectively, simple and double debonding interface. The debonding, introduced in these simulations as an extra resistance of the interface, is the reason of temperature jumps manifested in the temperature profiles (Figure 4), respectively, in one side and in two sides of inclusion according to debonding interface geometry (single partial debonding in Figure 2(a) or double partial debonding in Figure 2(b)).Temperature field in the heterogeneous media with. (a) Double partial debonding of circular inclusion; (b) single partial debonding of circular inclusion.
(a)
(b)Distribution of temperature following the vertical and horizontal cuts at the center of circular inclusion: (a) double partial debonding; (b) single partial debonding.
(a)
(b)The variation of the overall thermal conductivity of the considered heterogeneous material with temperature as well as with the volume fraction of constituents is shown in Figure5. Firstly, due to the decrease of the thermal conductivity of solid inclusion, the same tendency is obtained for the homogenized thermal property. The reduction between 300∘K and 500∘K of this property is much more pronounced in case of perfect interface with higher volume fraction of minerals inhomogeneity (Figure 5(a)). For the perfect interface case, the effective thermal conductivity is always higher than that of the clay matrix K1=1.8(W·m-1·K-1) because of a higher thermal conductivity of inclusion in all considered range of temperature. On the contrary, in the case of fully debonded interface, under the effect of thermal resistance of the interface, the effective thermal conductivity could be smaller than that of the matrix (Figure 5(b)).Effective thermal conductivityKeffyy as a function of temperature and volume fraction of inclusion: (a) case of perfect interface; (b) case of completely debonded interface.
(a)
(b)During thermal cracking, interfaces debonding is progressive. If this debonding is orientated and happened preferentially at some directions, then the overall properties could become anisotropic. Such situation is simulated by supposing a partially debonding of a circular inclusion. Figure6 illustrates the evolution of effective thermal conductivities in the vertical and horizontal directions in such situation. In both directions the thermal conductivity decreases when the debonding angle (angle ϕ2(b)) increases from ϕ=0 (perfect interface) to ϕ=π/2. However, because of the orientation of the debonded part of the interface, the thermal conductivity in the vertical direction decreases quicker than in horizontal direction introducing an anisotropy on macroscopic thermal conductivity of the homogenized material.Effective thermal conductivity as a function of temperature and debonding angle: (a) vertical directionKeffyy; (b) horizontal directions Keffxx.
(a)
(b)While single inclusion structure model assumes a periodic material, in cases when this hypothesis is not satisfied a more realistic description of the microstructure is needed. The IIM could be successfully used in such conditions where complex structures are studied. This feature is illustrated here by considering a clayey rock with 60 quartz inclusions counting for 20% of the volume. While all inclusions have the same diameter, their positions are randomly obtained by using a Monte Carlo procedure. As described in our previous work [18], the number of inclusion used in simulations is sufficient for having convergent results for homogenized thermal conductivity of the random medium.For illustration purposes, we consider a hypothetic cracking evolution described by a linear evolution of debonding angle with temperatureϕ=a×u+b with a=π/300 and b=-273×π/300, which means that the interfaces are perfect at 273∘K and completely debonded at 573∘K. Note also that the position of debonded part is randomly generated at the boundary of each inclusion.The contours of temperature field in elementary representative volume at some temperature levels are presented in Figure7, in which one could clearly distinguish the progressive debonding of inclusions illustrated by the increase of the debonded arc length of inclusions. Finally, as performed in Figure 8(a), the comparison between two cases of perfectly bonded and partially debonded of random inclusions has elucidated the important increase of the nonlinearity of effective thermal conductivity with temperature in the latter case. Further, due to the random distribution of circular inclusions as well as the debonded part of each inclusion, these numerical results always display the quasi-isotropic of the overall thermal property where Keffyy/Keffxx≈1 as observed in Figure 8(b).Contour of temperature in the heterogeneous media at different state of applied temperature.
(a)
(b)Effective thermal conductivity of the randomly heterogeneous media as a function of temperature: (a) result estimated in the vertical directionKeffyy; (b) anisotropy of the effective thermal conductivity Keffyy/Keffxx.
(a)
(b)


### 3.2. Pressure Dependence

In addition to the temperature dependence, the variation of the thermal conductivity of geomaterials with respect to applied pressure is also an important issue of research. Based on experimental results realized on different types of rock, many authors showed that the increase of pressure augments the thermal conductivity of materials [2, 27–29]. For example, in their work, Abdulagatova et al. [28] measured the thermal conductivity of sandstone at pressures up to 400 MPa. They found that the thermal conductivity of a rock increases rapidly between 0.1 and 100 MPa, and at high pressure (>100 MPa), a weak linear dependence of the thermal conductivity with pressure was observed. The increasing effect of pressure on thermal behavior was explained by the closure of microcrack which caused the approaching of the rock grains to each other at moderate pressure. If pressure is still further increased, the reduction of the rock’s intrinsic porosity may take place. These results confirm the observation presented in the previous work of Clauser and Huenges [1].The influence of pressure on the effective thermal conductivity will be carried out here by supposing that the applied pressure will close continuously the debonded part of inclusions. This simplification aims to illustrate the evolution of the thermal conductivity under pressure but obviously a more sophisticated study needs to be conducted to explicitly account for the mechanism of this closure of debonded inclusions. This case of study can be realized in the context of the thermomechanical coupled behaviour of geomaterials and will be discussed in our near future work. Thus, in the present paper, we will define the variation of debonding angle with respect to pressure as follows:ϕ=ϕ0(1+(P/P0)n)(1-n)/n. This function presents a rapid decrease of debonding angle with moderate pressure P and is inversely proportional with high pressure. Otherwise, to take into account the temperature effect, the initial debonding angle ϕ0 of random inclusions may vary as a function of temperature as considered in previous section. As an example, in Figure 9, we present the variation of the effective thermal conductivity of the heterogeneous rock with respect to pressure as well as temperature with the value P0=10 MPa and n=2. These results show the quick increase of the overall thermal conductivity between the atmospheric pressure and 100MPa, while a quasi linear evolution can be observed for the range of pressure from 100 MPa to 200 MPa.Effective thermal conductivity of the randomly heterogeneous media as a function of pressure and temperature.
(a)
(b)


### 3.3. Comparison with Experimental Results

In this subsection, we compare our numerical prediction with some experimental results published in the literature. As an example, in their contribution, Jobmann and Buntebarth [7] investigated the evolution of the thermal conductivity of the bentonite-quartz mixture with respect to temperature within 20∘C to 150∘C. Their main purpose of using this admixture is to accelerate the heat flow through the geotechnical barrier consisting of bentonite for the deep geological disposal of waste. To increase the thermal conductivity of the barrier, which needs to be similar to ones of the host rock, different quartz contents in the admixture have been tested. The grain size of quartz was about 0.3mm. The results of these measures showed that the thermal conductivity of the bentonite-quartz mixtures slightly increase with respect to the quartz content. Moreover, due to the temperature dependence of the thermal conductivity of quartz while the thermal property of bentonite is almost constant about 0.95(W·m-1·K-1), the overall thermal conductivity of the mixture decreases with respect to temperature.In Figure10, we present the comparison of the numerical and experimental results with different quartz content of the bentonite-quartz mixture. In these numerical simulations the thermal resistance of interface is equal to 2×10-4(m2·K·W-1) while the two coefficients of the debonding angle’s evolution function with temperature (see Section 3.1) are a=π/720 and b=-53×π/720, respectively. These values are obtained from an inverse analysis so that both experiment and simulation agree at best. We can state that the tendency observed in the experiment is well captured by the numerical estimations using the IIM method.Figure 10
Contour of temperature in the heterogeneous media at different state of applied temperature.


## 3.1. Temperature Dependence

Beginning from the pioneer work of [3], the temperature-dependent thermal properties of geomaterials have been studied for a long time in relation with various engineering applications such as geothermal reservoirs, underground storage of nuclear waste or petroleum and natural gas geology. This subject is always in continuous progress with a remarkable number of results conducted on different types of soils and rocks in the last decade [2, 4, 6–9]. A general tendency of decreasing of thermal conductivity as a function of temperature can be stated from these contributions. At the same time, a thermal damage is reported as a consequence of the contrast of constituent properties of geomaterials when the temperature increases. This type of damage known as thermal cracking developed principally in contacts of constituent phases of the heterogeneous media like geomaterials.In this part, we use the IIM method to solve the problem of effective thermal conductivity estimation of geomaterials whose constituent thermal conductivities depend on temperature. In addition, the effects of cracking interfaces with temperature on the effective thermal conductivity are considered using imperfect interface. In that case, the nonlinearity of effective thermal property reflects not only the nonlinearity of the intrinsic properties of matrix and inclusions but also the properties of interfaces, which when debonding produce an extra thermal resistance. In order to describe the progressive cracking of interfaces when the temperature evolves, the properties of interface are supposed to vary with coordinates following a known function. In the simplest case a piecewise interfacial resistance function that takes some values for debonding parts of interfaces and some other values for intact ones could be used. As an example, the debonding parts of the interface of a circular inclusion could be given in function of an angleϕ measured from an arbitrary chosen direction (Figure 2). In extreme cases, full interface could be considered as debonded. One could anticipate that, because of the temperature dependent conductivity of local constituents and interface debonding, the overall conductivity would be highly temperature dependent [14, 15, 22–26]. As demonstrated in a previous contribution from the authors [18], the contact resistance leads also to a size-dependent behavior of heterogeneous materials as observed in reality: for the same volumetric fraction of inclusions, a more pronounced effect of the interface is observed for smaller inclusions. If not otherwise stated in following simulations, the size of inclusion is equal to R=1μm.(a) Single partial debonding of circular inclusion. (b) Double partial debonding of circular inclusion.
(a)
(b)For the sake of clarity, the geometry as well as the initial and boundary conditions of the considered problem is depicted in Figure2. More precisely, a constant initial temperature u is considered in whole sample, and the mixed boundary conditions consist in applying the temperatures u and u+Δu at the two opposing faces (to the Oy direction in the examples presented here) of the rectangular domain, whereas the zero heat flux is supposed to be to two other faces in the orthogonal direction (faces normal to the Ox direction) to the first. In what follows, a constant 320×320 grid is used for all calculations. This value obtained from the analysis of the sensitivity of numerical results to the grid density in previous work [18] gives excellent agreement with the analytical solutions in the context of the linear problem.For illustration purposes, we consider a sample of Callovo-Oxfordian argillite (see, e.g., [18] for a brief description) constituted by a clay matrix with quartz inclusions at the mesoscopic level. The constant thermal conductivity of the matrix is taken as K1(u)=1.8(W·m-1·K-1). Otherwise the thermal conductivity of quartz significantly decreases with temperature and following evolution is reported in published works [1, 3]: K2(u)=7.7/(0.0045×u-0.3863)(W·m-1·K-1). It must emphasize however that despite simple expressions of thermal conductivities used for this material, one could use any expression provided by experimental considerations. Concerning the thermal resistance of interface in debonded parts, an arbitrary value is taken equal to 10-6(m2·K·W-1), while perfect interface conditions (α=0) are supposed for bonded parts.For simulations, firstly the sample of heterogeneous material with a single circular inclusions and partially (simple or double) debonded interface is considered. As the illustration in Figure3, the temperature fields are presented on the material containing a single circular inclusion with, respectively, simple and double debonding interface. The debonding, introduced in these simulations as an extra resistance of the interface, is the reason of temperature jumps manifested in the temperature profiles (Figure 4), respectively, in one side and in two sides of inclusion according to debonding interface geometry (single partial debonding in Figure 2(a) or double partial debonding in Figure 2(b)).Temperature field in the heterogeneous media with. (a) Double partial debonding of circular inclusion; (b) single partial debonding of circular inclusion.
(a)
(b)Distribution of temperature following the vertical and horizontal cuts at the center of circular inclusion: (a) double partial debonding; (b) single partial debonding.
(a)
(b)The variation of the overall thermal conductivity of the considered heterogeneous material with temperature as well as with the volume fraction of constituents is shown in Figure5. Firstly, due to the decrease of the thermal conductivity of solid inclusion, the same tendency is obtained for the homogenized thermal property. The reduction between 300∘K and 500∘K of this property is much more pronounced in case of perfect interface with higher volume fraction of minerals inhomogeneity (Figure 5(a)). For the perfect interface case, the effective thermal conductivity is always higher than that of the clay matrix K1=1.8(W·m-1·K-1) because of a higher thermal conductivity of inclusion in all considered range of temperature. On the contrary, in the case of fully debonded interface, under the effect of thermal resistance of the interface, the effective thermal conductivity could be smaller than that of the matrix (Figure 5(b)).Effective thermal conductivityKeffyy as a function of temperature and volume fraction of inclusion: (a) case of perfect interface; (b) case of completely debonded interface.
(a)
(b)During thermal cracking, interfaces debonding is progressive. If this debonding is orientated and happened preferentially at some directions, then the overall properties could become anisotropic. Such situation is simulated by supposing a partially debonding of a circular inclusion. Figure6 illustrates the evolution of effective thermal conductivities in the vertical and horizontal directions in such situation. In both directions the thermal conductivity decreases when the debonding angle (angle ϕ2(b)) increases from ϕ=0 (perfect interface) to ϕ=π/2. However, because of the orientation of the debonded part of the interface, the thermal conductivity in the vertical direction decreases quicker than in horizontal direction introducing an anisotropy on macroscopic thermal conductivity of the homogenized material.Effective thermal conductivity as a function of temperature and debonding angle: (a) vertical directionKeffyy; (b) horizontal directions Keffxx.
(a)
(b)While single inclusion structure model assumes a periodic material, in cases when this hypothesis is not satisfied a more realistic description of the microstructure is needed. The IIM could be successfully used in such conditions where complex structures are studied. This feature is illustrated here by considering a clayey rock with 60 quartz inclusions counting for 20% of the volume. While all inclusions have the same diameter, their positions are randomly obtained by using a Monte Carlo procedure. As described in our previous work [18], the number of inclusion used in simulations is sufficient for having convergent results for homogenized thermal conductivity of the random medium.For illustration purposes, we consider a hypothetic cracking evolution described by a linear evolution of debonding angle with temperatureϕ=a×u+b with a=π/300 and b=-273×π/300, which means that the interfaces are perfect at 273∘K and completely debonded at 573∘K. Note also that the position of debonded part is randomly generated at the boundary of each inclusion.The contours of temperature field in elementary representative volume at some temperature levels are presented in Figure7, in which one could clearly distinguish the progressive debonding of inclusions illustrated by the increase of the debonded arc length of inclusions. Finally, as performed in Figure 8(a), the comparison between two cases of perfectly bonded and partially debonded of random inclusions has elucidated the important increase of the nonlinearity of effective thermal conductivity with temperature in the latter case. Further, due to the random distribution of circular inclusions as well as the debonded part of each inclusion, these numerical results always display the quasi-isotropic of the overall thermal property where Keffyy/Keffxx≈1 as observed in Figure 8(b).Contour of temperature in the heterogeneous media at different state of applied temperature.
(a)
(b)Effective thermal conductivity of the randomly heterogeneous media as a function of temperature: (a) result estimated in the vertical directionKeffyy; (b) anisotropy of the effective thermal conductivity Keffyy/Keffxx.
(a)
(b)


## 3.2. Pressure Dependence

In addition to the temperature dependence, the variation of the thermal conductivity of geomaterials with respect to applied pressure is also an important issue of research. Based on experimental results realized on different types of rock, many authors showed that the increase of pressure augments the thermal conductivity of materials [2, 27–29]. For example, in their work, Abdulagatova et al. [28] measured the thermal conductivity of sandstone at pressures up to 400 MPa. They found that the thermal conductivity of a rock increases rapidly between 0.1 and 100 MPa, and at high pressure (>100 MPa), a weak linear dependence of the thermal conductivity with pressure was observed. The increasing effect of pressure on thermal behavior was explained by the closure of microcrack which caused the approaching of the rock grains to each other at moderate pressure. If pressure is still further increased, the reduction of the rock’s intrinsic porosity may take place. These results confirm the observation presented in the previous work of Clauser and Huenges [1].The influence of pressure on the effective thermal conductivity will be carried out here by supposing that the applied pressure will close continuously the debonded part of inclusions. This simplification aims to illustrate the evolution of the thermal conductivity under pressure but obviously a more sophisticated study needs to be conducted to explicitly account for the mechanism of this closure of debonded inclusions. This case of study can be realized in the context of the thermomechanical coupled behaviour of geomaterials and will be discussed in our near future work. Thus, in the present paper, we will define the variation of debonding angle with respect to pressure as follows:ϕ=ϕ0(1+(P/P0)n)(1-n)/n. This function presents a rapid decrease of debonding angle with moderate pressure P and is inversely proportional with high pressure. Otherwise, to take into account the temperature effect, the initial debonding angle ϕ0 of random inclusions may vary as a function of temperature as considered in previous section. As an example, in Figure 9, we present the variation of the effective thermal conductivity of the heterogeneous rock with respect to pressure as well as temperature with the value P0=10 MPa and n=2. These results show the quick increase of the overall thermal conductivity between the atmospheric pressure and 100MPa, while a quasi linear evolution can be observed for the range of pressure from 100 MPa to 200 MPa.Effective thermal conductivity of the randomly heterogeneous media as a function of pressure and temperature.
(a)
(b)


## 3.3. Comparison with Experimental Results

In this subsection, we compare our numerical prediction with some experimental results published in the literature. As an example, in their contribution, Jobmann and Buntebarth [7] investigated the evolution of the thermal conductivity of the bentonite-quartz mixture with respect to temperature within 20∘C to 150∘C. Their main purpose of using this admixture is to accelerate the heat flow through the geotechnical barrier consisting of bentonite for the deep geological disposal of waste. To increase the thermal conductivity of the barrier, which needs to be similar to ones of the host rock, different quartz contents in the admixture have been tested. The grain size of quartz was about 0.3mm. The results of these measures showed that the thermal conductivity of the bentonite-quartz mixtures slightly increase with respect to the quartz content. Moreover, due to the temperature dependence of the thermal conductivity of quartz while the thermal property of bentonite is almost constant about 0.95(W·m-1·K-1), the overall thermal conductivity of the mixture decreases with respect to temperature.In Figure10, we present the comparison of the numerical and experimental results with different quartz content of the bentonite-quartz mixture. In these numerical simulations the thermal resistance of interface is equal to 2×10-4(m2·K·W-1) while the two coefficients of the debonding angle’s evolution function with temperature (see Section 3.1) are a=π/720 and b=-53×π/720, respectively. These values are obtained from an inverse analysis so that both experiment and simulation agree at best. We can state that the tendency observed in the experiment is well captured by the numerical estimations using the IIM method.Figure 10
Contour of temperature in the heterogeneous media at different state of applied temperature.


## 4. Conclusions

The immersed interface method (IIM) was adapted in the present work to solve the elliptical transfer equation taking into account the contact resistance of interface in non linear heterogeneous composite-like geomaterials. This numerical tool is then used to derive the temperature- and pressure-dependent effective thermal conductivity of geomaterial with imperfect contact between matrix and inclusion. The imperfection of interface modeled as the partial or full debonding of inclusions is created during a thermophysical phenomenon called thermal cracking under heat load while a continuous closure of this imperfect interface is considered under pressure. The numerical results highlight that the overall conductivity at the macroscale depend strongly on the applied temperature and pressure through the state of debonding inclusions. The efficiency of numerical tool is demonstrated both in periodic structures and random ones allowing us to study in details temperature fields in a heterogeneous material with interfaces and its effective properties counting for the role of the microstructure complexity on the overall properties.



---

*Source: 456931-2013-07-07.xml*