Replacement relations for thermal conductivity of a porous rock

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

Contents lists available at ScienceDirect

International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms

Replacement relations for thermal conductivity of a porous rock a

c,⁎

b

b

MARK a

Fengjuan Chen , Yuri A. Popov , Igor Sevostianov , Raissa Romushkevich , Albert Giraud , Dragan Grgica a b c

GeoRessources Laboratory, Université de Lorraine (ENSG), CNRS, CREGU, F-54501 Vandoeuvre-les-Nancy, France Center for Hydrocarbon Recovery, Skolkovo Institute of Science & Technology, Moscow, Russia Department of Mechanical and Aerospace Engineering, New Mexico State University, Las Cruces, NM 88003, USA

A R T I C L E I N F O

A BS T RAC T

Keywords: Replacement relations Thermal conductivity Sandstone Non-ellipsoidal inhomogeneity

We focus on derivation and experimental verification of the replacement relations that link the overall thermal conductivities of heterogeneous materials with the same matrix and microstructure but having inhomogeneities with different properties. First, we derive replacement inequalities based on Hashin-Shtrikman bounds that relate overall thermal conductivity of a composite to thermal conductivity of a.porous material. HashinShtrikman bounds are also used to derive analogy of Gassmann equation for thermal conductivity. Then, we use formalism of property contribution tensors to obtain replacement relations for anisotropic materials containing ellipsoidal inhomogeneities. These relations coincide for different homogenization schemes – non-interaction approximation, Mori-Tanaka scheme, and Maxwell scheme, provided, that both effective conductivity of a composite and porous material are calculated in the framework of the same method. In the case of the overall isotropy, all of these relations coincide with Gassmann equation derived from Hashin-Shtrikman bounds. We check the possibility to apply these relations to 3-D non-ellipsoidal inhomogeneities on example of a supersphere using numerical simulations. The replacement relations are approximate with satisfactory accuracy for convex superspheres, while the error is significant for concave shapes. For this case, we suggest a modification that involves an extra shape factor that can be determined, for example, from comparison of the average Eshelby tensor for conductivity and conductivity contribution tensor of a pore. To verify the approach, thermal conductivity of 85 quartz sandstone specimens – dry and saturated with water and kerosene – of porosity varying from 0.14 to 0.29 is measured using optical scanning technique. The average pore shape and thermal conductivity of the dense quartz matrix are determined from best fitting of the conductivity-porosity curve for dry sandstone. Then these parameters are used in the replacement relation for sandstone saturated with water and kerosene. The comparison of the experimental data with theoretical predictions shows a good accuracy of the proposed approach.

1. Introduction In the present paper we formulate the replacement relations used in the evaluation of overall conductive properties to the case of thermal conductivity. The replacement relations link overall properties of heterogeneous materials that have the same matrix material and microstructure while properties of the inclusions are different. They can be used to predict the change in overall thermal conductivity of a porous material upon the saturation. The replacement relations have been first proposed by Gassmann1 in the context of the effect of saturation on seismic properties of rock in geomechanics. He proposed to express the bulk modulus of fully saturated rock in terms of the elastic properties of dry rock. Brown and Korringa2 generalized Gassmann equation for inhomogeneous ⁎

Corresponding author. E-mail address: [email protected] (I. Sevostianov).

http://dx.doi.org/10.1016/j.ijrmms.2017.06.008 Received 16 February 2017; Received in revised form 12 June 2017; Accepted 22 June 2017 1365-1609/ © 2017 Elsevier Ltd. All rights reserved.

anisotropic material. Further discussion has been done by Han and Batzle3 who mentioned that several physical quantities (porosity, density and speed velocity) should be consistent and constrained and proposed to use Voigt-Reuss bounds and critical porosity limits constraints to get upper and lower bounds for the fluid-saturation effect. Ciz and Shapiro4 obtained the replacement relations for the general case of anisotropic two-phase material with anisotropic constituents and specified it for bulk and shear moduli of the isotropic twophase composite (see also5 for corrected equation). Sevostianov and Kachanov6 formulated the replacement relations for anisotropic materials in terms of the property contribution tensor for elasticity problem. They showed that these relations are exact for ellipsoids and approximately accurate for certain non-ellipsoidal shapes (a cube and various 2D shapes). Chen et al.7 further developed

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

2. Replacement relations based on Hashin-Shtrikman bounds

their results checking the applicability of the relations to 3-D nonellipsoidal inhomogeneities on example of a supersphere (the shape is mathematically described by equation x1 2p + x2 2p + x3 2p = 1). They showed that the replacement relations can be used to estimate with a satisfactory accuracy the effective elastic properties of a material containing inhomogeneities of convex shape, the error of estimation for concave inhomogeneities is significant. Saxena and Mavko8 discussed the replacement relations under the assumption that the elastic fields inside the inhomogeneities are uniform (overall properties and properties of the constituents are isotropic). Note that this assumption is actually equivalent to the statement that the inhomogeneities are ellipsoids subjected to the uniform external field9–11. Saxena and Mavko12,13 discussed the impactful applications of these relations in geophysics, they highlighted the importance of these relations in evaluating the effective properties of a heterogeneous material from those of a porous material and the benefices of calculating effective properties using Eshelby tensor. In this context, the difference between the Eshelby tensor and the property contribution tensor need to be clarified. Eshelby tensor, referred in the first Eshelby problem, is a dimensionless quality that interrelate the resulted elastic field to the eigenstrain that would have been exist inside the inclusion. Property contribution tensor, referred in the second Eshelby problem, is used to give a quantitative description of the contribution of the inhomogeneity into the overall properties. For ellipsoids, the first and the second Eshelby problems are mathematically equivalent, all the tensors (Eshelby tensor, Hill tensor, property contribution tensor etc.) are uniform inside the inclusion/ inhomogeneity (position independent) and they are linearly interrelated. For non-ellipsoids, the interrelation between the Eshelby tensor and the property contribution tensor is not valid, and the Eshelby tensor can’t be used in the calculation of the overall properties. In the context of thermal conductivity problem, Schärli and Rybach14 compared air saturated and water saturated low-porosity granitic rocks and reported that thermal conductivity of water-saturated samples is 30% higher than “dry” conductivities. Zimmerman15 evaluated thermal conductivity of air- and water- saturated rock using Fricke formula16 for electrical conductivity of a material containing randomly oriented spheroidal inhomogeneities (that represents generalization of Maxwell formula17). In particular, Zimmerman15 proposed a methodology to evaluate thermal conductivity of water saturated rock from the dry rock measurements and validated the prediction on experimental data for sedimentary and granitic rock data. He stated, that the procedure requires inversion of the equations and, while to of the required equations can be inverted in closed form, the final results can be obtained only numerically. In the present paper, we derive the replacement relations analytically in closed form. First we use Hashin-Shtrikman bounds18 to obtain replacement inequalities and analogy of Gassmann equation1 for the case of thermal conductivity of isotropic material. Then we use formalism of resistivity/conductivity contribution tensors to derive these relations for anisotropic materials containing ellipsoidal inhomogeneities. The replacement relations constitute an important connection between Eshelby tensor (the first Eshelby problem) and the resistivity/conductivity contribution tensor (the second Eshelby problem). We check the possibility to use replacement relations for materials containing non-ellipsoidal inhomogeneities using example of a superspherical shape7. We show that, similarly to the elastic case, the replacement relations can be used as a good approximation for materials with convex inhomogeneities, while application to concave shapes lead to a serious error. In such cases, we propose to use modified replacement relation that involves a shape factor that can be determined experimentally. The approach is verified by comparison with experimentally measured thermal conductivities of dry and saturated sandstone. To the best of our knowledge, replacement relations have never been derived for anisotropic materials and their applicability has never been checked against the pore shape.

Hashin-Shtrikman bounds for effective thermal conductivity keff of an isotropic statistically homogeneous two-phase material (consisting of two isotropic phases with conductivities k 0 and k1) have the following form18:

k1 +

1−ϕ ϕ 3k 1

+

1 k 0 − k1

≤ keff ≤ k 0 +

ϕ 1−ϕ 3k 0

+

1 k1 − k 0

(2.1)

Where ϕ is volume concentration of the phase with conductivity k1. For a porous material k1 = 0 and the lower bound vanishes while the upper bound can be written as

k0 +

ϕ 1−ϕ 3k 0

+

1 k1 − k 0

⎛ 3ϕ ⎞ = k 0 ⎜1 − ⎟ 2 + ϕ⎠ ⎝

(2.2)

Thus, the Hashin-Shtrikman bounds for overall conductivity kdry of a porous material can be written as

⎛ 3ϕ ⎞ 0 ≤ kdry ≤ k 0 ⎜1 − ⎟ 2 + ϕ⎠ ⎝

(2.3)

It leads to the following bounds for porosity ϕ :

ϕ≤

2(k 0 − kdry ) 2k 0 + kdry

; 1−ϕ≥

3kdry 2(k 0 + kdry )

(2.4)

Substitution of the expression (2.4) into the lower HashinShtrikman bound (2.1) yields the following lower bound for the effective conductivity of a composite material in terms of the conductivity of porous material having the same structure:

⎡ 2k 0 (k 0 + 2k1) + kdry (7k 0 − 4k1) ⎤ ⎥ ≤ k1 + k1 ⎢ ⎣ 2k 0 (k 0 + 2k1) − kdry (2k 0 − 5k1) ⎦

1−ϕ ϕ 3k 1

+

1 k 0 − k1

≤ keff (2.5)

Similarly, the upper bound for the effective conductivity can be written as

keff ≤ k 0 +

ϕ 1−ϕ 3k 0

+

1 k1 − k 0

⎡ 2k 0 k1 + kdry (2k 0 − k1) ⎤ ⎥ ≤ k0 ⎢ ⎢⎣ ⎥⎦ 2k 02 + kdry k1

(2.6)

Combining (2.5) and (2.6), we can write the replacement bounds Hashin-Shtrikman bounds for effective conductivity of a two-phase material in terms of conductivities of two phases and conductivity of a porous material (having porosity equal to the volume concentration of the phase with conductivity k1)

⎡ 2k 0 k1 + kdry (2k 0 − k1) ⎤ ⎡ 2k 0 (k 0 + 2k1) + kdry (7k 0 − 4k1) ⎤ ⎥ ⎥ ≤ keff ≤ k 0 ⎢ k1 ⎢ ⎥⎦ 2k 02 + kdry k1 ⎣ 2k 0 (k 0 + 2k1) − kdry (2k 0 − 5k1) ⎦ ⎣⎢ (2.7) These bounds can be used, for example, to evaluate thermal conductivity of saturated rock if properties of the dry one, as well as properties of two phases are known while the porosity is not known. Fig. 1 compares bounds (2.7) with the Hashin-Shtrikman bounds (2.1) for the case of porous calcite either saturated with crude oil or having pores filled with saturated clay (the properties are given in Table 1). As expected, the wideness of resulting bounds is about the same as those of the original bounds (2.1). In particular, when the thermal conductivity of two constituents, k 0 and k1, differ two times (conductivities of calcite and saturated clay), the upper and lower bounds almost coincide (Fig. 1a and b). Hashin-Shtrikman upper bounds for effective conductivities of two statistically isotropic composites can be used to derive the analogy of Gassmann's equation for conductivity. Indeed, the upper bound

keff = k 0 +

65

ϕ 1−ϕ 3k 0

+

1 k1 − k 0

(2.8)

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

Fig. 1. Comparison of the replacement bounds (2.7) (figures (a) and (c)) and Hashin-Shtrikman bounds (2.1) (figures (b) and (d)) for calcite saturated with crude oil and calcite with pores filled by saturated clay. k 0 is the bulk thermal conductivity of the calcite, k dry is the effective thermal conductivity of the porous calcite; keff is the effective thermal conductivity of the calcite with filled pores and 1 − ϕ is the volume fraction of the calcite.

ductivity of an isotropic two phase material. After some algebra, Eq. (2.12) can also be written for thermal resistivities:

Table 1 Thermal conductivities of constituents used for illustration in Figs. 1–3.

reff =

k 0 [35]

koil [36]

kclay

3.30W /m⋅K

0.12W /m⋅K

1.52W /m⋅K

keff k 0 − keff

2k 0 + k1 − 2ϕ (k 0 − k1) 3(k 0 − k1)

=

(2.9)

For a composite with one insulating phase (porous material, k1 = 0 ) it is reduced to

ϕ

kdry k 0 − kdry

=

2 (1 − ϕ) 3

keff k 0 − keff

−ϕ

kdry k 0 − kdry

3. Property contribution tensors and replacement relations for anisotropic materials containing ellipsoidal inhomogeneities

=

k1 k 0 − k1

(2.11) 3.1. Replacement relations for a single inhomogeneity were derived in 6 for elastic properties. They used the concept property contribution tensors, first introduced as compliance contribution tensors in the context of pores and cracks in 21. For the thermal (or electrical conductivity) problem, the resistivity contribution tensors were introduced in 22. Following this work, we consider a homogeneous material (matrix) of volume V with thermal resistivity tensor r0 containing an inhomogeneity of volume V* with thermal resistivity tensor r1, subjected to a remotely applied heat flux on the infinite boundary. For a twophase material, assuming the temperature gradient ∇T is linearly related to the heat flux vector q per volume (Fourier law) for both the constituents, second rank resistivity contribution tensor R of an inhomogeneity relates the changes in average over V temperature gradient due to the presence of the inhomogeneity:

This expression can be written in the explicit form for keff as a function of kdry and porosity ϕ :

keff = kdry +

(2.13)

(2.10)

Subtraction of (2.10) from (2.8) yields

ϕ

ϕ (r0 − rA) + (r0 − rdry )

Remark.. Eqs. (2.12) and (2.13) are derived using the upper HashinShtrikman bounds for thermal conductivities of dry and saturated materials. However, these bounds coincide with the predictions of effective field homogenization schemes for a material containing spherical pores/ inhomogeneities (Maxwell scheme and Mori-Tanaka scheme, see, for example,20). In this context, Eqs. (2.12) and (2.13) (as well as the original Gassmann equation) are based on the assumptions that the morphology of the material is formed by isolated spherical inhomogeneities. Below, we show that the ellipsoidal shape of the inhomogeneities does not violate the form of the equations.

can be rewritten after some algebra as

ϕ

ϕ rdry (r0 − rA) + r0 (r0 − rdry )

k1 (k 0 − kdry ) 2 (1 − ϕ) k 0 k1 − k1 kdry + ϕ k 02

(2.12)

Note that the pore shape is not entering expressions (2.11) and (2.12) explicitly. It appears, however, if kdry is expressed as a function of porosity. Relation (2.12) repeats the form of Gassmann equation1 for effective bulk modulus of fully saturated rock in terms of the corresponding elastic properties of dry rock (if the bulk moduli of all the materials are replaced by thermal conductivities). The same equation can also be obtained using cross-property connection of Lutz and Zimmerman19 that links bulk modulus and thermal con66

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

V Δ (∇T ) = * R⋅q V

NIA r eff = r0 + ϕ RA

(3.1)

(superscript NIA indicates the non-interaction approximation). In particular, the effective resistivity of a porous material is given by:

Dual to R conductivity contribution tensor K is used to express the extra heat flux q (per reference volume V ) when the material is subjected to a remotely applied temperature gradient ∇T on the remote boundary.

V Δq = * K⋅∇T V

NIA r dry = r0 + ϕ Rpore

NIA NIA r eff = r0 + ϕ (r dry − r0)

For an ellipsoidal inhomogeneity, R and K can be written in explicit form (3.3)

∫V

*

NIA = reff

∂G (x − x′) d x′ ∂x′ j

(3.4)

G (x − x′) =

1 = 4π r

is Green's function for conductivity pro-

−1 −1 − (r − r )−1 R−1 B 0 A − RB = (rA − r0)

(3.6)

K−1 A

(3.7)

−

= (k A −

k 0)−1

− (kB −

k 0)−1

−1

RA = [(rA − r0)−1 + R−1 pore]

−1

MTB r dry = r0 +

1 V

∑ V*n R(n)

NIA (1 − ϕ) k 0 k1 − k1 kdry + ϕ k 02

(3.15)

(3.16)

ϕ Rpore 1−ϕ

(3.17)

Combining (3.16) and (3.17) with the replacement relations (3.8) yields

(3.9)

−1 −1

MTB MTB r eff = r0 + ϕ [(rA − r0)−1 + ϕ (r dry − r0) ]

For materials containing multiple inhomogeneities, relations (3.6)– (3.9) lead to generalization of Gassmann's equation to the general anisotropic case. Below, we consider three homogenization schemes – non-interaction approximation, Mori-Tanaka scheme, and Maxwell scheme. 3.2. Non-interaction approximation, being rigorous at infinitesimal concentration of inhomogeneities, shows good accuracy for concentration of inhomogeneities up to 15%25. Besides, it serves as a basic building block for various homogenization schemes (self-consistent, differential, Mori-Tanaka's etc.) by placing non-interaction inhomogeneities into some sort of “effective environment” – effective field or effective media. In the framework of this approximation, contributions of individual inhomogeneities into overall properties are just summed up. So that the change in the thermal resistivity due to the presence of the inhomogeneities is

Δr =

NIA 2 ) k1 (k 0 − kdry

(superscript MTB indicates Mori-Tanaka-Benveniste method). In particular, for a dry porous material,

(3.8)

]

(3.13)

(3.14)

MT r eff = r0 + ϕ RA⋅[ϕ (rA − r0)−1RA + (1 − ϕ) I]−1

In the opposite case of a superconducting materials (kB → ∞), expression (3.7) takes the form

KA = (kA − k 0)−1 + K−1 pore

−1

3.3. Mori-Tanaka-Benveniste scheme is widely used due to its simplicity and good accuracy predictions. It belongs to the class of effective field homogenization schemes where each inhomogeneity, treated as a single one, is placed into the unaltered matrix material; interactions are accounted for by assuming that the inhomogeneity is subjected to the field that differs from the remotely applied one. The basic idea of the method has roots in works of Mossotti (see chapter 11 in 26). The Mori-Tanaka-Benveniste scheme27,28 is based on the assumption that the effective field acting on each inhomogeneity is equal to the average over the matrix. Then the macroscopic properties may be calculated from the non-interaction approximation with appropriate change of the remotely applied field. Then, the effective resistivity for a material containing inhomogeneities “A” can be expressed as

In particular, for an insulating material (rB → ∞), the replacement relation (3.6) is simplified to:

[

NIA ϕ (r0 − rA) + (r0 − rdry )

NIA NIA keff = kdry +

blem. To the best of our knowledge it was first derived for the general case of anisotropy in 23. Barthélémy24 proposed an alternative method of derivation using transformation of coordinates and also derived explicit expressions for components of tensor Pij in the case of arbitrarily oriented ellipsoid in orthotropic matrix. Note that tensors Q and P depend on the shape of the inhomogeneity (aspect ratios of the ellipsoid) and thermal conductivities of the matrix material. If we consider now two inhomogeneities, A and B, with different properties, but having the same shape and embedded in the same matrix, expression (3.3) applied to these two inhomogeneities leads, after elimination of hill's tensor, to the following replacement relations:

K−1 B

]

NIA − r0)−1 (r dry − r0) + ϕ I

NIA NIA ϕ rdry (r0 − rA) + r0 (r0 − rdry )

(3.5)

det r 0 er ⋅ r0 ⋅ er

A

In terms of conductivities, this expression has the same form as (2.12)

and the second Hill's tensor Q that is related to P by 0 Qij = kim [δmj − Pmn k nj0 ]

[ (r

In particular, for an isotropic material with randomly oriented NIA is identical ellipsoidal inhomogeneities, the effective resistivity reff NIA expressed in terms of the effective resistivity of dry material rdry having the same microstructure as

in terms of the first Hill's tensor

∂ Pij = ∂xi

(3.12)

Using replacement relations (3.8), one obtains

(3.2)

R = [(r1 − r0 )−1 + Q] −1 and K = [ (k1 − k0 )−1 + P] −1

(3.11)

(3.18)

For an isotropic material, (3.18) is reduced to MTB reff = r0 −

MTB ) φ (r0 − rA) (r0 − rdry MTB (r0 − rdry ) + ϕ (r0 − rA)

=

MTB MTB ϕrdry (r0 − rA) + r0 (r0 − rdry ) MTB ϕ (r0 − rA) + (r0 − rdry )

(3.19) or, in terms of thermal conductivities, MTB MTB keff = kdry +

MTB 2 ) k1 (k 0 − kdry MTB (1 − ϕ) k 0 k1 − k1 kdry + ϕ k 02

(3.20)

3.4. Maxwell scheme was originally proposed by Maxwell17 for electrical conductivity of a material containing randomly located spherical inhomogeneities is probably the oldest homogenization scheme. We use its interpretation proposed29,30, where Maxwell scheme is formulated in terms of property contribution tensors. Applications of Maxwell homogenization scheme to calculation of the effective thermal conductivity of isotropic and anisotropic rock have been developed in 15 and 31.

(3.10)

For a material containing identical inhomogeneities “A” of arbitrary conductive properties, the effective resistivity is expressed as: 67

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

The basic idea of the method is that the far field produced by the considered set of inhomogeneities is equated to the far field produced by a fictitious domain Ω of certain shape that possesses the (yet unknown) effective properties. Then, applying expressions (3.3) for this domain, the effective resistivity for a material containing inhomogeneities “A” can be expressed as M r eff = r0 + ϕRA⋅[I − ϕQ Ω⋅RA]−1

4.1. Calculation of Eshelby tensors for a superspherical inclusion in the conductivity problem Eshelby tensor SC for conductivity problem does not have clear physical meaning (in contrast with the elasticity problem). This secondrank tensor relates a prescribed heat flux-free temperature gradient ∇T * inside an inclusion (occupying domain Ω ) with the resulting temperature gradient. The latter is non-uniform if Ω has non-ellipsoidal shape. Then the average value has to be considered:

(3.21)

(subscript M indicates Maxwell scheme). For a dry material, it takes form M r eff

= r0 + ϕRpore⋅[I −

ϕQ Ω⋅Rpore]−1

∇T = SC ⋅∇T *

(3.22)

with

(note that Q Ω in (3.21) and (3.22) is the same if morphologies of two materials coincide). Combining the last two expressions with the relation (3.6) give us −1 −1

M M r eff = r0 + ϕ [(rA − r0)−1 + ϕ (r dry − r0) ]

C

Sij = −

M ϕ (r0 − rA) (r0 − rdry )

(r 0 −

M rdry

) + ϕ (r0 − rA)

=

1 V*

ni (x) dS (x)

(4.3)

denotes thermal where overbar stays for average over the volume, conductivity of the matrix, n is unite normal vector outward the surface of the inclusion, G (x , x′) is Green's function for conductivity problem that gives the temperature at point x produced by a unit point heat sources at point x′. For an infinite isotropic medium,

(3.23)

M M ϕrdry (r0 − rA) + r0 (r0 − rdry ) M ϕ (r0 − rA) + (r0 − rdry )

1 1 1 4π k 0 ρ (x − x′)

(4.4)

(x1 − x′1) 2 + (x2 − x′2 ) 2 + (x3 − x′3 ) 2

(4.5)

G (x − x′) = − (3.24) with

or, equivalently, M M keff = kdry +

⎞

⎛

∫∂ Ω (x) ⎜⎝ ∫ ∂ Ω (x′) kkj0 G (x, x′) nk (x′) dS (x′) ⎟⎠

kkj0

For an isotropic material, it takes the following form: M reff = r0 −

(4.2)

ρ (x − x′) =

M 2 ) k1 (k 0 − kdry

(1 − ϕ) k 0 k1 −

M k1 kdry

+ϕ

k 02

so that

(3.25) 7

Note, that, similarly to the elasticity problems , the replacement relations given by non-interaction approximation, Mori-Tanaka method and Maxwell scheme coincide with each other and with one provided by Gassmann's Eq. (2.12). As seen, from the structure of formulas (3.15), (3.20), and (3.25) different homogenization schemes produce the same replacement relations provided that the properties of materials with dry and filled pores are calculated with the same method. It is related to the fact that the effective properties in the mentioned schemes are expressed in terms of the resistivity contribution tensor for a single inhomogeneity (for a pore, this tensor is just an inverse of the second Hill's tensor Q ). We now illustrate the replacement relation on example of porous calcite matrix filled by two different materials: saturated clay (solid) and crude oil (liquid) (the properties of the materials are given in Table 1). Fig. 2(a) and (c) show dependences of the isotropic effective thermal conductivity of materials with filled pores on corresponding thermal conductivity of the dry material at ϕ = 0.10, 0.25 and 0.50. Dependences on volume fraction are illustrated in Fig. 2(b) and (d) for kdry / k 0 = 0.80, 0.60 and 0.40.

C

Sij =

1 4π V*

⎛

∫∂ Ω (x) ⎜⎝∫∂ Ω (x′)

⎞ dS (x′) ⎟ ni (x) dS (x) ρ (x , x′) ⎠ n j (x′)

Volume of supersphere V* is calculated as: ⎡ ⎛ 1 ⎞ ⎤3 2 ⎢Γ ⎜ 2p ⎟ ⎥ ⎣ ⎝ ⎠⎦ V* = ⎛3⎞ 3p 2 Γ ⎜ 2p ⎟ ⎝ ⎠

(4.7) C S11

C S22

C S33

= = and Due to the cubic symmetry of the supersphere, C C C Sij = 0 if i ≠ j , only S11 need to be calculated. Fig. 4 shows that S11 is almost independent of p. 4.2. Calculation of the resistivity contribution tensor for a superspherical pore and verification of the replacement relations Eshelby and Hill's tensors for conductivity problems are interrelated as C

Pmn = S mi⋅rin0 4. Applicability of the replacement relations to nonellipsoidal inhomogeneity: case of the supersphere

(4.8)

For a pore (r1 → ∞) in isotropic material (kij0 = k 0 δij ), substitution of (4.8) into (3.5) and (3.3) yields the following expression for the resistivity contribution tensor

For a non-ellipsoidal shape, Hill's tensors P and Q are not uniform inside the inhomogeneity,P = P (x), Q = Q (x) and relations (3.3) do not take place. However, in some cases, these expressions can be used approximately with tensors P = P (x), Q = Q (x) being replaced by their volume averages P and Q. It also leads to the approximation of the replacement relations (3.6) and (3.7) since they follow from the representation (3.3). In this next section we check the accuracy of such approach for an inhomogeneity of a superspherical shape:

x1 2p + x2 2p + x3 2p = 1

(4.6)

C −1

k 0 Rij = [δij − Sij ]

(4.9)

Numerical evaluation of this tensor for a superspherical pore is given in 32 (see Appendix); it is showed that the tensor can be approximated with good accuracy by the following expression

k 0 Rij =

(4.1)

3(5p − 1) ⎡ V0 ⎤ ⎥ δij = ηδij ⎢ ⎣ V* ( p ) ⎦ 8

(4.10)

(where V0 is the volume of the unit sphere). Fig. 5 provides comparison of the resistivity contribution tensor calculated using (4.9) and (4.6) with the results of 32. We observe that the agreement is reasonably good for convex supersphere ( p > 0.5), however, for concave supersphere, the relative error is significant. It means that the replacement relations (3.6), (3.7) can be applied for convex shapes only.

where parameter p describes the derivation of the shape from a sphere. The shape is concave when p < 0.5 and convex when p > 0.5; it becomes a sphere when p = 1. Hereafter parameter p is called concavity parameter. Fig. 3 illustrates examples of the superspherical shape for different values of p . 68

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

Fig. 2. Dependence of the isotropic effective thermal conductivity keff of a porous material filled with crude oil or saturated clay on corresponding thermal conductivity k dry of the dry material (figures (a) and (c)) and on volume fraction of the inhomogeneities ϕ (figures (b) and (d)).

Generally, for non-ellipsoidal shapes, effects of the material contrast and shape of the inhomogeneity cannot be separated from each other. To take it into account we formulate the following hypothesis

with numerically calculated one. Eq. (4.11) leads to the following connection between thermal conductivities of dry and saturated materials:

Hypothesis.. Property contribution tensors R and K are related to average Hill's tensors Q and P by the following relations

reff = r0 + ϕ [α (rA − r0)−1 + ϕ (rdry − r0)−1]−1

R = α [(r1 − r0 )−1 + Q ] −1 and K = α [(k1 − k0 )−1 + P ]

−1

(4.12)

(4.11)

where shape-dependent parameter α can be obtained from comparison of resistivity contribution tensors of a pore obtained by this equation

For isotropic material, this hypothesis leads to the following replacement relation

Fig. 3. Superspherical shapes at different values of the concavity parameter p.

69

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

Fig. 6. Dependence of parameter α entering (4.11)–(4.14) on the shape factor of supersphere p.

shown that quartz content in the rock ranges within 90–95%. The rock specimens have fine-grained texture, sometimes with admixture of gravel granules, and indistinct cross-bedded, often massive structure. Microstructure of the specimens is illustrated in Fig. 7. The specimens were represented by core plugs of the dimensions 30×30 mm, drilled from full sized cores (recovered from a vertical well to provide standard petrophysical laboratory measurements). A core plug drilling direction was perpendicular to the full size core axis. An optical scanning technique had been used for the thermal conductivity measurements. An important peculiarity of the optical scanning technique is that it allowed us to perform non-contact nondestructive measurements just on standard core plugs without any mechanical treatment of them to keep them for other petrophysical measurements and investigations. The measurement technique provided the determination of the principal thermal conductivity tensor components from optical scanning in two perpendicular directions33. Two scans in mutually perpendicular directions were performed on every flat surface (two bottoms) of every core plug. Thermal conductivity profiles along every scan line were recorded during the measurements, and average thermal conductivity values were determined from every profile. The measurement results for parallel scans on both flat surfaces were averaged to be processed further for determination of the thermal conductivity tensor components. 2D model of rock anisotropy (when azimuthal anisotropy is neglected) was used to determine two principal thermal conductivity tensor components - along and perpendicular to bedding plane. Total errors (uncertainties) of the thermal conductivity measurements estimated from special metrological control on certified thermal conductivity references were estimated as ± 2,5% for the thermal conductivity tensor component along a bedding plane, ± 5% for the thermal conductivity tensor component perpendicular to the bedding plane, and ± 6% (for confidential level of 0.95 in all the cases) for thermal anisotropy coefficient. The thermal anisotropy coefficient was calculated for every rock sample as a ratio of the principal thermal conductivity tensor components parallel and perpendicular to the bedding plane of a rock sample. The heat penetration and correspondingly the effective thickness of a rock sample layer involved in every measurement depend on the measurement regime parameters and rock thermal diffusivity value33. The thermal conductivity of quartz in aggregate is 7.60 W/(m·K)34, that causes relatively high values of the thermal conductivity established for rock samples studied. At the measurement regime parameters used for the collection study the effective thickness of rock samples ranged within 6–12 mm and from the measurements on both flat surfaces of core plugs resulted in 12–24 mm totally for every core plug. The thermal conductivity measurements were carried out sequentially for core plugs (1) dried in an oven, (2) saturated with water (a brine model), and (3) saturated with kerosene that can be considered

Fig. 4. Dependence of the (isotropic) average Eshelby tensor for conductivity problem C S11 for a superspherical inclusion on concavity parameter p.

Fig. 5. Comparison of the resistivity contribution tensor R of a superspherical pore calculated using (4.9) and (4.6) with the numerical results of 32 in dependence on the concavity parameter p.

reff =

αϕ rdry (r0 − rA) + r0 (r0 − rdry ) αϕ (r0 − rA) + (r0 − rdry )

(4.13)

We will verify the workability of this relation using experimental data obtained on dry and saturated quartz sandstone approximating the pore shape by superspheres. Note, first of all, that, for any shape, trace of the average Eshelby tensor for conductivity problem is equal to one. It means that for any 2 1 isotropic shape Pij = 3 δij , and Qij = 3 δij . Then, for a superspherical pore, dependence of α on the concavity parameter p is given by

α=

3 2η

(4.14)

with η defined in (4.10). This dependence is illustrated in Fig. 6. 5. Thermal conductivity measurements The thermal conductivity measurements have been performed on a set of specimens of quartz sandstones. Mineralogical analysis has 70

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

Fig. 7. Microphotographs of colored thin cross-sections of sand stone specimens (light color corresponds to the grains of quartzite, dark – porous space. The images correspond to the following specimens in Table 2: (a) specimen 42; (b) specimen 48, (c) specimen 25, (d) specimen 84. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 and in Table 2. Porosity values were determined from weighting rock samples in their dry, water- and kerosene-saturated states. Accounting for essential number rock samples under study, the determination coefficients R2 calculated (see Fig. 8) should be estimated as statistically significant for all states of rock samples (dry, water- and kerosene-saturation) as in all cases R2 exceed essentially a corresponding statistical bound for 85 elements. It is seen from Fig. 8, that correlation fields for three states of rock samples (dry, water-saturated and kerosene-saturated) differ in their compactness: the correlation field for water-saturated rock samples is the most compact, and the correlation field for dry rock samples is the less compact. An assumption could be done that a decrease in contrast of porous fluid and rock matrix (from dry to kerosene-saturated and to water-saturated rock samples) causes a decrease in influence of variations in pore space geometry within the rock collection on effective thermal conductivity of rock samples. 6. Analysis of the experimental data and verification of the replacement relation (4.13) Fig. 8. Measured thermal conductivities of dry and saturated sandstone specimens in dependence on the volume fraction of pores ϕ , and best fit exponential curves.

At the first step, we used the data on thermal conductivity of dry rock to evaluate parameters of the model – average shape of the pores and thermal conductivity k 0 of the solid phase. From the analysis of photomicrographs (Fig. 7) we chose a supersphere to approximate the average pore shape. To evaluate the concavity parameter p we used Maxwell homogenization scheme (3.20). A custom MATLAB script was developed to determine this parameter using nonlinear least squares fitting. The Trust Region Reflective method was used to minimize the squared second norm of the vector e denoting the error between the overall conductivity of the dry material calculated using (3.20) and the experimental data:

as a thermal model of oil as the thermal conductivity values for kerosene and oil are very close (0.12–0.14 W/(m K) at atmospheric temperature). For dry core plugs the thermal anisotropy coefficient values were established to be not more than 1.10 for 57 core plugs from 89 core plugs studied. For four core plugs thermal coefficient values were estimated to be larger than 1.20 that could be caused by their crossbedded structure. The thermal anisotropy coefficient of porous/fractured rocks reduces normally after saturation of the rock samples with water and kerosene as thermal conductivity of water and kerosene (0.60 W/(m K) and 0.12 W/(m K) correspondingly) are in smaller contrast with the rock matrix thermal conductivity than the air thermal conductivity (0.025 W/(m K))33. The results of the measurements for 85 core plugs are presented in

min

e

( p)

where 71

2 2

= min

M kdry ( p,

( p)

M kdry ( p , ϕ) − k exp (ϕ)

k exp

2

2

(6.1)

ϕ) is the vector of the overall thermal conductivity of the

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

Table 2 (continued)

Table 2 Experimental data of effective conductivity of sandstones.

Sample Sample

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70

Thermal conductivity, W/ (m K)

Dry

Watersaturated

Oil (kerosene)saturated

1.53 2.51 2.71 2.42 2.23 2.33 1.90 2.44 1.96 2.22 1.90 1.82 3.35 1.98 1.72 1.80 1.54 1.90 2.00 1.78 1.78 1.83 1.86 1.87 1.50 1.75 1.96 1.78 2.34 1.71 1.68 1.78 1.71 1.84 1.71 1.73 1.91 1.91 1.57 1.86 1.98 1.96 1.90 1.62 1.99 1.79 1.85 1.72 1.84 2.07 1.74 2.18 2.60 2.01 1.98 2.11 2.17 2.34 1.73 1.61 1.67 1.82 1.65 1.85 1.72 1.89 1.69 2.13 2.05 1.65

4.11 4.05 4.35 3.78 4.07 3.97 3.85 3.77 4.01 3.96 3.66 3.52 4.65 4.05 3.61 3.78 3.75 3.74 3.92 3.94 3.81 3.99 3.84 3.87 4.05 3.86 3.73 3.56 4.10 3.92 3.73 3.88 3.92 3.86 3.87 3.97 3.82 3.93 3.76 3.77 3.74 3.60 3.77 3.83 3.62 3.95 3.87 3.98 3.90 3.87 3.92 3.90 4.23 3.82 4.06 3.72 4.02 3.97 3.96 3.74 3.76 3.60 3.82 3.73 3.94 3.86 3.89 3.71 4.07 3.83

2.86 3.30 3.30 3.03 3.06 3.09 2.88 3.39 2.83 2.94 2.61 2.55 3.92 2.87 2.65 2.85 2.52 2.75 2.69 2.96 2.65 2.79 2.89 2.57 2.77 2.62 2.80 2.51 3.23 2.71 2.47 2.82 2.88 2.70 2.74 2.77 2.70 2.71 2.60 2.75 2.74 2.85 2.73 2.66 2.84 2.71 2.66 2.77 2.78 2.68 2.77 2.84 3.39 2.95 2.90 2.88 2.95 2.85 2.75 2.61 2.50 2.69 2.63 2.42 2.62 2.70 2.62 2.67 2.98 2.55

Thermal anisotropy coefficienta

Porosity

1.00 1.00 1.15 1.00 1.00 1.00 1.00 1.17 1.00 1.00 1.11 1.00 1.18 1.14 1.11 1.11 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.18 1.00 1.00 1.00 1.16 1.18 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.16 1.00 1.12 1.00 1.14 1.00 1.00 1.15 1.00 1.00 1.00 1.00 1.00 1.11 1.00 1.00 1.00 1.11 1.00 1.00 1.00 1.18 1.00 1.12 1.11 1.00

0.261 0.210 0.203 0.219 0.243 0.202 0.251 0.220 0.233 0.231 0.264 0.259 0.139 0.260 0.279 0.220 0.262 0.256 0.265 0.219 0.231 0.253 0.236 0.275 0.261 0.242 0.270 0.252 0.240 0.233 0.266 0.256 0.253 0.264 0.261 0.256 0.267 0.257 0.273 0.281 0.269 0.244 0.284 0.277 0.247 0.260 0.257 0.258 0.259 0.238 0.249 0.213 0.212 0.237 0.248 0.235 0.224 0.213 0.251 0.241 0.263 0.246 0.258 0.254 0.255 0.240 0.268 0.257 0.247 0.242

71 72 73 74 75 76 77 78 79 80 81 82 83 84 85

Thermal conductivity, W/ (m K)

Dry

Watersaturated

Oil (kerosene)saturated

1.96 1.74 2.00 2.46 2.22 2.75 2.19 2.11 2.22 2.15 1.95 1.50 1.59 1.53 2.75

4.08 3.73 3.95 4.18 3.95 4.48 3.91 3.91 4.01 3.95 3.80 4.05 3.91 4.03 4.60

2.83 2.67 2.95 3.18 2.90 3.82 2.86 2.91 2.79 2.89 2.78 2.77 2.69 2.86 3.74

Thermal anisotropy coefficienta

Porosity

1.12 1.00 1.00 1.00 1.00 1.00 1.15 1.00 1.18 1.00 1.00 1.00 1.19 1.00 1.20

0.253 0.250 0.241 0.228 0.267 0.148 0.265 0.267 0.257 0.272 0.286 0.259 0.242 0.259 0.170

a Thermal anisotropy coefficient Λ = k// /k⊥ . k// and k⊥ are thermal conductivities parallel and perpendicular to rock bedding. respectively.

dry material calculated by (3.20) for the given concavity parameter p and the porosity ϕ , k exp (ϕ) is the vector of experimentally measured conductivities, and k exp is the mean value of k exp (ϕ). The step of variation of the concavity parameter was chosen as 0.01. We observed that independently of the value of k 0 , this parameter varies as 0.27 ≤ p ≤ 0.31 with the best fit provided by p ≈ 0.30 (Fig. 9). Then, we used thus obtained value of p and repeated the procedure to evaluate k 0 . The obtained best fit is provided by k 0 = 6.80 W / m⋅K that is slightly less than the value reported in literature for dense quartz34,35. This difference may be explained by the effect of the thermal resistance of the intergrain contacts. The examples of fitting curves are presented in Fig. 10. Thus obtained values of k 0 and p are used to check the hypothesis formulated at the end of Section 4 and to validate the replacement relation (4.13). The value of the shape parameter α according to (4.14) is 0.434. Using now the thermal conductivities of oil and water (Table 3), we obtain the curves shown in Fig. 11. The regression analysis allow us to state that the agreement is statistically significant for both types of rock saturation (water and kerosene). 7. Conclusions We derived and experimentally validated the replacement relations that link the overall thermal conductivities of heterogeneous materials having the same microstructure and matrix material but filled with different inhomogeneities. First, we derived replacement inequalities based on Hashin-Shtrikman bounds18 and showed how analogy of Gassmann equation1 can be obtained for thermal conductivity. Then we used formalism of property contribution tensors to obtain replacement relations for anisotropic materials containing ellipsoidal inhomogeneities. These relations coincide with Gassmann equation in the case of isotropy. We checked the possibility to apply these relations to 3-D non-ellipsoidal inhomogeneities on example of a supersphere. The replacement relations are approximate with satisfactory accuracy for convex superspheres, while the error is significant for concave shapes. For this case, we suggested a modification (4.13) that involves additional shape dependent parameter that can be determined for nonellipsoidal shapes if Eshelby tensor and compliance contribution tensor of the pore are known. To verify the workability of thus modified replacement relation, we measured thermal conductivity of dry sandstone and sandstone saturated with water and kerosene and compared the model predictions with the experimental data inferred from the measurements of thermal conductivity of sandstone – dry and satu72

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al.

Fig. 9. Best fit for the average concavity parameter p based on the experimental data on dry sandstone.

Fig. 11. Verification of the replacement relation (4.13): comparison of the analytical prediction for keff with experimentally measured thermal conductivities of sandstone saturated with kerosene and water; ϕ is the volume fraction of (fully saturated) pares.

References 1 Gassmann F. Über die elastizität porpöser medien. Vierteljahrsschr der Nat Gesellscaft Zur. 1951;96:1–23. 2 Brown RJS, Korringa J. On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid. Geophysics. 1975;40:608–616. 3 Han D-H, Batzle ML. Gassmann's equation and fluid-saturation effects on seismic velocities. Geophysics. 2004;69:398–405. 4 Ciz R, Shapiro SA. Generalization of Gassmann equations for porous media saturated with a solid material. Geophysics. 2007;72:A75–A79. 5 Shapiro SA. Fluid-Induced Seismicity, Cambridge University Press; 2015. 6 Sevostianov I, Kachanov M. Relations between compliances of inhomogeneities having the same shape but different elastic constants. Int J Eng Sci. 2007;45:797–806. 7 Chen F, Sevostianov I, Giraud A, Grgic D. Accuracy of the replacement relations for materials with non-ellipsoidal inhomogeneities. Int J Solids Struct. 2017;104– 105:73–80. 8 Saxena N, Mavko G. Impact of change in pore-fill material on P-wave velocity. Geophysics. 2014;79:D399–D407. 9 Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc Roy Soc Lond. 1957;A241:376–396. 10 Lubarda VA, Markenscoff X. On the absence of Eshelby property for ellipsoidal inclusions. Int J Solids Struct. 1998;35:3405–3411. 11 Rodin GJ. Eshelby's inclusion problem for polygons and polyhedral. J Mech Phys Solids. 1996;44:1977–1995. 12 Saxena N, Mavko G. Exact equations for fluid and solid substitution. Geophysics. 2014;79:L21–L32. 13 Saxena N, Mavko G. The embedded-bound method for estimating the change in rock moduli under pore fill and mineral phase substitution. Geophysics. 2015;80:L1–L10. 14 Schtirli U, Rybach L. On the thermal conductivity of low-porosity crystalline rocks. Tectonophysics. 1984;103:307–313. 15 Zimmerman RW. Thermal conductivity of fluid-saturated rocks. J Pet Sci Eng. 1989;3:219–227. 16 Fricke H. A mathematical treatment of the electric conductivity and capacity of disperse systems. Phys Rev. 1924;24:575–587. 17 Maxwell JC. A Treatise on Electricity and Magnetism, Oxford: Clarendon Press; 1873. 18 Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behavior of multiphase materials. J Mech Phys Solids. 1962;11:127–140. 19 Lutz MP, Zimmerman RW. Effect of an inhomogeneous interphase zone on the bulk modulus and conductivity of a particulate composite. Int J Solids Struct. 2005;42:429–437. 20 Markov KZ. Elementary micromechanics of heterogeneous mediaMarkov KZ, Preziozi L, eds. Heterogeneous Media: Micromechanics Modeling Methods and Simulations, Boston: Birkhauser; 2000, pp.1–162. 21 Horii H, Nemat-Nesser S. Overall moduli of solids with microcracks: load-induced anisotropy. J Mech Phys Solids. 1983;31:155–171. 22 Sevostianov I, Kachanov M. Explicit cross-property correlations for anisotropic twophase composite materials. J Mech Phys Solids. 2002;50:253–282. 23 Willis JR. Bounds and self-consistent estimates for the overall properties of anisotropic composites. J Mech Phys Solids. 1977;25:185–202. 24 Barthélémy J-F. Effective permeability of media with a dense network of long and micro fractures. Transp Porous Med. 2009;76:153–178.

Fig. 10. Best fit for k 0 based on the experimental data on dry sandstone (p=0.30).

Table 3 Parameters used for validation of the replacement relations.

k0

koil

kwater

p

α

6.80W /m⋅K

0.12W /m⋅K

0.60W /m⋅K

0.30

0.434

rated with water and kerosene. The obtained agreement is quite good and statistically significant for both types of rock saturation (water and kerosene). Acknowledgement Financial support from the FP7 Project TAMER IRSES-GA-2013– 610547 and New Mexico Space Grant Consortium contained in the NASA Cooperative Agreement NNX13AB19A to New Mexico State University are gratefully acknowledged. 73

International Journal of Rock Mechanics & Mining Sciences 97 (2017) 64–74

F. Chen et al. 25 Sevostianov I, Sabina F. Cross-property connections for fiber reinforced piezoelectric materials. Int J Eng Sci. 2007;45:719–735. 26 Feynman RP, Leighton RB, Sands M. The Feynman Lectures on Physics, Reading, MA: Addison-Wesley; 1964. 27 Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall. 1973;21:571–574. 28 Benveniste Y. On the Mori-Tanaka method for cracked solids. Mech Res Commun. 1986;13:193–201. 29 Sevostianov I. On the shape of effective inclusion in the Maxwell homogenization scheme for anisotropic elastic composites. Mech Mater. 2014;75:45–59. 30 Sevostianov I, Giraud A. Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape. Int J Eng Sci. 2013;64:23–36. 31 Giraud A, Sevostianov I, Chen F, Grgic D. Effective thermal conductivity of oolitic rocks using the Maxwell homogenization method. Int J Rock Mech Min Sci. 2015;80:379–387.

32 Chen F, Sevostianov I, Giraud A, Grgic D. Evaluation of the effective elastic and conductive properties of materials containing concave pores. Int J Eng Sci. 2015;97:60–68. 33 Popov Yu, Beardsmore G, Clauser C, Roy S. ISRM suggested methods for determining thermal properties of rocks from laboratory tests at atmospheric pressure. Rock Mech Rock Eng. 2016;49:4179–4207. 34 Popov Yu, Berezin V, Soloviov G, et al. Thermal conductivity of minerals. Izv Phys Solid Earth. 1987;23:245–253. 35 Clauser C, Huenges E. Thermal conductivity of rocks and mineralsAhrens T, ed. . A Handbook of Physical Constants: Rock Physics and Phase Relations.., Washington: The American Geophysical Union; 1995, pp.105–126. 36 Elam SK, Tokura I, Saito K, Altenkirch RA. Thermal conductivity of crude oils. Exp Therm Fluid Sci. 1989;2:1–6.

74