Intergranular pressure solution in chalk: a multiscale approach

Computers and Geotechnics 34 (2007) 291–305 www.elsevier.com/locate/compgeo

Intergranular pressure solution in chalk: a multiscale approach D. Lydzba a, S. Pietruszczak b

b,*

, J.F. Shao

c

a Wroclaw University of Technology, Wroclaw, Poland Department of Civil Engineering and Engineering Mechanics, McMaster University, Hamilton, Ont., Canada L8S 4L7 c Laboratoire de Mecanique de Lille, Villeneuve d’Ascq Cedex, France

Available online 6 April 2007

Abstract This paper outlines a multiscale approach for modelling of chemo-mechanical interaction in chalk. The primary focus is on quantitative micromechanical analysis of dissolution/diffusion phenomena that occur within the intergranular interface. First, the influence of the kinetics of the dissolution process on the underlying physical mechanisms is investigated. Subsequently, a macroscopic chemoplasticity formulation is outlined. Extensive numerical studies are carried out. Those include the analysis at micro-level, simulation of macroscopic creep deformation characteristics as well as application in a boundary value problem that involves a water injection test.  2007 Elsevier Ltd. All rights reserved. Keywords: Chalk; Intergranular pressure solution; Chemical degradation

1. Introduction Modelling of short and long-term behaviour of porous chalk is of significant interest in the areas of mining and petroleum engineering. Examples of practical applications involve a stability analysis of an underground cavity, subsidence prediction of oil reservoir due to water injection and infiltration, etc. Over the last few decades, extensive experimental studies have been performed combined with theoretical investigations aimed at describing the main features of the mechanical behaviour. In spite of intense research efforts, some basic issues still remain open. Those pertain to description of physico-chemical mechanisms associated with water saturation as well as with the longterm creep deformation. The mechanical response of chalk is complex and depends on many factors, such as mineralogical composition of the skeleton, chemical composition of the pore fluid, temperature and the stress state. In this work, the consideration is primarily focused on the so-called pure chalk, composed of more than 98% of CaCO3 and with less than 0.8% *

Corresponding author. Tel.: +1 905 525 9140; fax: +1 905 524 2121. E-mail address: [email protected] (S. Pietruszczak).

0266-352X/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2007.02.002

of SiO2. The porosity of a typical pure chalk can range from 40% to 45%. The intrinsic permeability varies from 1 · 1016 m2 to 1 · 1015 m2, while the size of calcite grains ranges from 0.1 lm to 1 lm (Schroeder [1]). Due to high porosity, the basic mechanical behaviour of chalk is strongly affected by the confining pressure (e.g., Halleux et al. [2], Andersen [3], Monjoie et al. [4], Schroeder et al. [5], Homand and Shao [6]). Under a low confinement, the behaviour is brittle-elastic or brittle-elastoplastic. In this case, the failure of the sample is generally accompanied by the formation of localized shear bands, whose orientation depends on the value of confining pressure. With the increase in confinement, a progressive transition from brittle to ductile behaviour is observed. At high confining pressures, no peak stress is observed until a very large value of axial strain is attained (generally, more than 10%). The slope of stress–strain curve is continuously increasing resulting in a concave characteristic, similar to that observed in hydrostatic compression test. The failure of the sample is the result of a complete destruction of the pore structure; the cohesive chalk is transformed, after a significant volumetric deformation, into a compacted powder. Mechanical behaviour of chalk is strongly influenced by the chemical composition of the pore fluid. This aspect has

292

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

been extensively examined through numerous experimental investigations (e.g., Brignoli et al. [7], Andersen [3], Schroeder et al. [5], Risnes and Flaageng [8], Homand and Shao [6], Schroeder [1], Xie [9]). By comparing the response of oil and water-saturated samples, it has been demonstrated that the presence of water triggers a significant degradation in mechanical properties. In general, both the deformation and strength characteristics are affected; the axial strength may be reduced by as much as 50%. In addition to standard triaxial tests, the water injection tests have also been conducted on samples initially saturated with oil (Homand and Shao [6], Schroeder [1]). In this case, an abrupt increase in strain has been observed at the locations adjacent to the propagating water injection front. The effects of water saturation are associated with different time scales. For example, while the degradation of strength is almost instantaneous, the creep deformation is a long-term phenomenon. In order to describe the evolution of characteristics in the presence of water, various constitutive formulations have been employed. Among these, the most common approach is perhaps that based on the concept of suction that originates from the mechanics of unsaturated soils (cf. Papamichos et al. [10], Collin et al. [11]). In this approach, the decrease in mechanical strength is attributed to the reduction in suction that is defined as the differential pressure between fluids, i.e. water and oil. This concept, although simple and computationally attractive, cannot adequately explain the drastic decrease in mechanical strength triggered by the presence of water and/or the long-term process of creep deformation. Furthermore, the duality of the time scale associated with both these phenomena cannot be accounted for. The creep phenomenon itself, in chalks and other rock materials, has been typically handled within the phenomenological framework of viscoplasticity (cf. Cristescu [12], Jin and Cristescu [13], De Gennaro et al. [14]). Again, this approach cannot explain the underlying physical mechanism or the effect of chemical composition of the pore fluid on the process. In recent years there have been several attempts to describe the chemo-mechanical interaction in terms of intergranular pressure solution occurring at the grain contacts. An example here is the work of Hellemann et al. [15], who had studied the kinetics of the creep process and its relation to water chemistry, temperature, etc. The authors concluded that the stationary creep phenomenon in chalk is inherently related to the slow dissolution of interfaces between grains and that the dissolution is enhanced by temperature. More recently, Heggheim et al. [16] have studied the phenomenon of calcium dissolution by seawater with a high sulphate concentration. The mechanical strength of chalk saturated with the modified seawater was 20–25% lower than that of samples saturated with artificial Ekofisk formation brine. The degradation of chalk properties was then discussed in terms of a chemical dissolution/precipitation process and the chemistry of the thin water film close to inter-granular contacts. Further, Risnes et al. [17] studied the weakening of chalk skeleton by water–glycol mix-

tures. The authors demonstrated that the mechanical strength of the material increases with decreasing water activity. They suggested a mechanism that involved an additional pressure exerted on the grains by attraction of water molecules to the chalk surface. This surface adsorption was said to be responsible for an increase in local pore pressure and the associated reduction in the cohesion of the material. By evaluating different mechanisms related to chalk–water interaction, the authors concluded that pressure solution and adsorption pressure may be the two most important factors affecting the time dependent deformation of water saturated chalk. Based on the experimental information as outlined above, Pietruszczak et al. [18] have recently proposed a plasticity formulation describing the chemo-mechanical coupling in chalk. Both macro- and micromechanical frameworks were developed attributing the time-dependent effects to the intergranular pressure solution process (IPS). Two distinct time scales were introduced; one related to a rapid dissolution of contacts leading to degradation of mechanical strength and the other one associated with the subsequent much slower process of diffusion/dissolution along the interface boundaries. The purpose of this work is to enhance the formulation given in [18] by conducting an in-depth quantitative micromechanical analysis of the dissolution/diffusion phenomena. The paper is written in the following sequence. First, the micromechanics of IPS is examined and some numerical examples are given. In particular, the influence of kinetic of the process on the associated physical mechanism is investigated. Subsequently, the macroscopic formulation is outlined with the focus on the long-term creep phenomenon. Again, some numerical examples are provided. The paper is concluded by finite element simulation of a water injection test as conducted by Homand [19]. 2. Micromechanics of IPS 2.1. Mathematical description The actual contact is assumed to take place along isolated islands that are surrounded by a network of interconnected channels containing fluid (e.g. Raj [20], Spiers et al. [21]). The IPS process involves the dissolution of the solid phase, diffusion of the dissolved material along the interface as well as precipitation at the grain boundaries that are in contact with the pore fluid. The mathematical description of the process employs a continuum formulation that is obtained by incorporating the averaging – first, at the molecular level and, subsequently, at the level of representative volume element (RVE). The thickness of the intergranular layer is in the range of nanometers, which is more than 10 times the diameter of a single molecule. In view of this, the transport phenomena may be described by employing the molecular averaging (Murdoch [22]). Let RVEm (yi) be a region of the averaging at the molecular level. Here, yi is a position vector defining

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

the location of the geometric center of RVEm with respect to the origin, Fig. 1a. Assume that the representative volume contains N molecules at a given time t. Each molecule a (a = 1, . . ., N), located at the position zai , has a mass of ma and the velocity vai . The averaging is carried out by employing an ‘envelope’ function wðzai  y i Þ defined as 8 > < 1 if zai 2 RVEm ðy i Þ 3 a wðzi  y i Þ ¼ 4=3rðmÞ ð1Þ > : a 0 if zi 62 RVEm ðy i Þ where r(m) is the radius of RVEm. The above representation allows to assign to a given material point yi a set of averaged physical quantities. In particular, the density and the linear momentum vector, both averaged over RVEm, are defined as qðy i ; tÞ ¼

N X

ma wðzai  y i Þ

ð2Þ

a¼1

pi ðy i ; tÞ ¼

N X

ma vai ðtÞwðzai  y i Þ

ð3Þ

a¼1

Note that N N X oq owðzai  y i Þ X owðzai  y i Þ a ðy i ; tÞ ¼ ¼ ma ma z_ k ot ot ozak a¼1 a¼1

¼

N X a¼1

ma vak

owðzai  y i Þ oy k

ð4Þ

and N N opi ðy; tÞ X oðwvai Þ X ow ¼ ma ¼ ma vai oy i oy oy i i a¼1 a¼1

ð5Þ

In Eq. (7), the terms density and velocity both refer to the fluid filling the intergranular contact space. This fluid is actually a mixture of water and the dissolved substances (in particular, chalk). Assuming that the fluid is stationary, the mass conservation law for the dissolved material may be written in the form analogous to Eq. (7), i.e. oc oJ i þ ¼0 ot oy i

ð8Þ

where c is the concentration of the dissolved substance (in mass per unit volume) and Ji is the flux vector. The intergranular area has an evolving microstructure in which the actual contact between grains occurs along a set of isolated islands that are surrounded by fluid. An idealized geometry of the intergranular interface is shown in Fig. 2. The interface is perceived as a cylindrical slab with thickness d and the radius R. The islands at which the contact between grains is said to take place are idealized as cylinders of the radius j. Furthermore, it is assumed that the microstructure of the interface, prior to the dissolution process, is periodic; i.e. the geometry is obtained by invoking periodicity of the RVE. Thus, the topology of RVE consists of a solid cylinder of radius j bounded by adjacent grain surfaces and submerged in fluid, Fig. 2c. During the IPS process, the dissolution takes place along the phase transition boundaries C. These represent the surfaces of solids (grains/islands) which are in a direct contact with fluid. Along the moving phase transition front, the continuity of flux must be satisfied. The latter requires that the flux of the dissolved substance be equal to the rate of dissolution per unit area. Thus ½qs  ðqf þ cÞv þ J ¼ 0

Thus, defining the macroscopic velocity as

ð9Þ s

p ðy ; tÞ vi ðy i ; tÞ ¼ i i qðy i ; tÞ

ð6Þ

the relations (4) and (5) lead to the standard form of the averaged mass conservation law oqðy i ; tÞ þ divðqvi Þ ¼ 0 ot

ð7Þ

y3

293

f

In the equation above, q and q denote the intrinsic densities of the solid and liquid phase, respectively. Furthermore J = Jini represents flux component in the direction normal to the interphase boundary and v = vini is the speed of the interphase boundary. A continuum formulation of the transport phenomenon, valid over the entire intergranular area, can be obtained by employing the spatial (volume) averaging within RVE. The

z3

x3

y3

vα r

r zα

yi RVE m(yi) yi

r(m)

xi y2

x2

z2

y2

RVE m(0) y1

z1

RVE(xi)

RVE(0) x1

y1

Fig. 1. Upscaling by the molecular and volume averaging methods.

294

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

Γ

R

κ v δ

Fig. 2. Idealized geometry of intergranular interface.

averaging scheme incorporates the following transformation rules (Cushmann [23])   Z oJ i oh J i i 1 þ J dS ð10Þ ¼ oxi V RVE C oy i   Z oc oh c i 1  cvi ni dS ð11Þ ¼ ot ot V RVE C In the expressions above, the symbol Z 1 hi ¼  dV V RVE V f

ð12Þ

represents the volume averaging operator; VRVE is the volume of RVE, while Vf is the volume of the fluid within RVE. Furthermore, xi is a macroscopic space variable and yi is the microscopic variable defined within RVE, Fig. 1b. Introducing the transformation rules (10) and (11) in the mass conservation law (8), and utilizing the boundary condition (9), yields Z ohci ohJ i i 1 þ ¼ ðqs  qf Þv dS ð13Þ ot oxi V RVE C Assuming now that the average value of the flux vector is defined by Fick’s law hJ i i ¼ D

ohci oxi

leads to ohci 1  Dr2 hci ¼ ot V RVE

ð14Þ Z

ðqs  qf Þv dS

ð15Þ

C

where $2 is the Laplace operator and D is the effective coefficient of diffusion. Eq. (15) represents a smoothened form of the mass conservation law for the dissolved substance. Note that the dissolution process along the phase transition boundaries C is described here in terms of mass sources that are smeared over the representative volume. The velocity of the moving interphase boundary is governed by the kinetics of dissolution process. The latter is assumed to be proportional to the ‘distance’ from the state of chemical equilibrium, i.e. qs v ¼ bðy i Þðc  hciÞ

ð16Þ

where c is the saturation concentration (solubility limit). Note that the left-hand-side of Eq. (16) represents the rate of change of the solid mass per unit surface area. Thus, Eq. (16) is analogous to classical kinetics laws that govern the rate of dissolution/precipitation (cf. Jeschke and Dreybrodt [24]). Also note that hci > hci corresponds to dissolution process, while the inverse implies precipitation. The proportionality coefficient b is, in general, stressdependent so that its value is a function of the location at which the dissolution process takes place. The latter stems from the fact that the stress distribution along the interface is not uniform, thus prompting different rates of dissolution. Substituting the constitutive Eq. (16) in the relation (15) yields, after some algebraic transformations: Z ohci ðqs  qf Þ 1  Dr2 hci ¼ ð c  hciÞ bðy i Þ dS ð17Þ ot qs V RVE C For the axisymmetric geometry as depicted in Fig. 2, Eq. (17) can be expressed in terms of polar coordinates as    t  oc 1 oc o2 c qs  qf A ðr;tÞ Ar ðr;tÞ br t D þ þ ¼ s ðc  cÞb 2 ot r or or2 A A bt qd ð18Þ Here, At is the area of the upper/lower surfaces that are in contact with fluid, Ar refers to outer surface of the solid islands, while A is the total cross-sectional area of RVE. Furthermore, bt, br are the coefficients that govern the rate of the dissolution at the upper/lower phase transition boundaries and along the outer surface of cylindrical contacts, respectively. Note that, for clarity of presentation, the averaging symbol Æ æ has been omitted here. As mentioned earlier, an individual contact is represented by a cylinder with the radius j and the height d, so that At ðr; tÞ Ar jðr; tÞ þ ¼1 A 2Ad

ð19Þ

The current value of j(r, t) is defined through Eq. (16), i.e. qs

oj ¼ br ðc  cÞ ot

which leads to

ð20Þ

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

jðr; tÞ ¼

8 <

jðr; 0Þ 

:

br qs

Z

t

ðc  cÞ dt

0

0 6 t 6 td

ð21Þ

0 t P td Here td is referred to as the ‘life-span’ of the contact and its value, at a fixed location r, is defined by the relation Z br td jðr; 0Þ  s ðc  cÞ dt ¼ 0 ð22Þ q 0 where j(r, 0) = j0 is the initial radius at the onset of dissolution which, in view of the periodicity of the structure, is uniform in the entire domain. Introduce now the following geometric measure j2 ðr; tÞ 1ðr; tÞ ¼ 1  ð23Þ j20 Apparently, given the relation (21), the value of 1 is a function of concentration as well as the time, i.e. 8  2 Z t > br < 1 1 ðc  cÞ dt 0 6 t 6 td 1ðr; tÞ ¼ j0 qs 0 > : 1 t > td ð24Þ The above variable is a local measure which describes the progressive reduction in the effective contact area. The associated global measure f, which represents the average value of 1, is defined as Z R 1 f¼ 2 21ðr; tÞr dr ð25Þ R 0 Utilizing now the expressions (19) and (23), the differential equation (18) can be transformed to      oc 1 oc o2 c 2bt qs  qf D þ 2 ¼ ð c  cÞ 1  g0 ð1  1ðr;tÞÞ ot r or or d qs  d br pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  1ðr;tÞ ð26Þ þg0 j0 bt where g0 ¼ pj20 =A. Eq. (26) provides an averaged description of the dissolution and diffusion phenomena within the intergranular interface. In general, the mechanism of the pressure dissolution can be described as follows. At grain-to-grain contacts, the solids dissolve due to high stress concentrations. The dissolved substance is transported into the pore space, where the precipitation takes place at the grain boundaries as a result of oversaturation of solutes in the pore fluid. The description of the transport phenomenon and the precipitation in the pore space invokes: – mass conservation ocp oJ pi þ ¼0 ot oxi

ð27Þ

– continuity of flux along the interphase boundary between the grain and the pore fluid ½qs  ðqf þ ^cp Þvp þ J p ¼ 0

ð28Þ

295

– continuity condition at the boundary between the pore space and the intergranular interface undergoing IPS  oc p ð29Þ J D  ¼0 or r¼R – kinetics of precipitation at the grain boundary in contact with pore fluid qs vp ¼ bp ðcp  ^cp Þ

ð30Þ

In the expressions above, the superscript ‘p’ refers to the pore space; J p ¼ J pi ni is the component of the flux vector normal to the grain boundary and vp is the speed of the grain interphase boundary. Moreover, ^cp is the value of concentration cp at the grain boundaries, while cp is the solubility limit. Finally, bp is the coefficient governing the kinetics of the precipitation/dissolution along the grain boundaries and the  oc term D  denotes the flux of the dissolved substance or r¼R transported out of the intergranular contact. In what follows, the main focus is on the assessment of the evolution of concentration of the dissolved substance at the periphery of the intergranular interface. This is required in order to complete the formulation of the initial boundary value problem, as defined by Eq. (26). Here, a simplified approach is pursued for the description of the transport phenomena within the pore space, viz. Eqs. (27)–(30), which incorporates averaging process together with some explicit approximations. Note first that by averaging Eq. (27) over the representative pore space and utilizing Eqs. (28) and (29), the following expression is obtained:  p ohcp i 2pRd oc p kS k þD p ¼0 ð31Þ  ðq  q Þv s f ot kV k or  kV p k r¼R

Here, iV pi is the volume of the pore space and iSpi is the surface area of grains in contact with the pore fluid. Incorporating now the law of kinetics of the process (30), Eq. (31) can be transformed to  p ohcp i 2pRd oc ðqs  qf Þ p p p kS k ^ ¼ D p ¼0 þ b ð c  c Þ ot kV k or  qs kV p k r¼R

ð32Þ Note that the representation above employs an implicit assumption that the concentration is uniform along the grain boundary. Furthermore, it is assumed that there are no additional flux sources such as, for example, those resulting from the filtration through the pore space. An approximate form of Eq. (32) is obtained by assuming that the average concentration in the pore space may be defined using a linear estimate based on the value at the periphery of the intergranular contact c(R) and that at the grain boundaries ^cp , i.e. ð33Þ hcp i ¼ xcðRÞ þ ð1  xÞ^cp where x is a proportionality constant. Similarly, the flux at the grain boundary can be estimated as

296

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

ocp cðRÞ  ^cp ffi D on Rp p where R is an average radius of the grain.

J p ¼ D

ð34Þ

Apparently, there is kS p k ffi ðRp Þ2 ;

kV p k ffi ðRp Þ3

ð35Þ

Using the approximations above, together with the boundary condition (28), Eq. (32) can be simplified to the form  ocðRÞ 2pRd oc ½x þ cð1  xÞ ¼  p 3 D  ot ðR Þ or r¼R s

þ

f

q q p p b ðc  cðRÞÞ qs Rp

ð36Þ

where qs  qf p b qs ð37Þ c¼1 qs  qf p p D=R þ b qs The above representation defines, in an approximate sense, the evolution of the concentration of the dissolved substance at the periphery of the intergranular contact. In case when the void space is filled with a non-aqueous solvent, such as oil, the intergranular interface is in the state of thermodynamic equilibrium, i.e. the concentration of the dissolved substance is at the solubility limit. The injection of an aqueous solvent triggers the IPS process, i.e. dissolution and diffusion of the dissolved material. During this process the microstructure of the interface evolves, as the individual contacts reach their life-span. The initial boundary-value problem is defined in terms of Eqs. (26) and (36), subject to the following boundary and initial conditions: cðR; tÞ ¼ 0 for t ¼ 0 oc ð38Þ ¼ 0 for r ¼ 0; t > 0 or cðr; tÞ ¼ c for r 2 ½0; RÞ; t ¼ 0 2.2. Numerical analysis Apparently, the general formulation of the IPS phenomenon is complex; the underlying axisymmetric problem, viz. Eq. (26), is strongly non-linear and involves evolving microstructure with vanishing mass sources. In view of this complexity, the solution has been obtained numerically by employing the finite difference scheme. The segment [0, R] has been divided into m equal intervals, each with the length of Dx = R/m. The nodes were numbered sequentially from 1 at r = 0 to m + 1 at r = R. An explicit finite difference approximation of Eq. (26), with the boundary condition (38b), takes the form ckþ1 ¼ ck1 þ 2F ðck2  ck1 Þ þ Bk1 ðc  ck1 ÞDt ð39Þ 1 kþ1 cj ¼ ckj þ F ðckjþ1  2ckj þ ckj1 Þ F ðck  ck Þ þ Bkj ðc  ckj ÞDt; j ¼ 2;...;m þ 2ðj  1Þ jþ1 j1 ð40Þ

where

 ðqs  qf Þbt dbr k k 2 1 þ g ¼ 2 j  g ðj Þ Bkþ1 0 0 j j qs d j0 bt j 8 r < jk  b ðc  ck ÞDt if jkþ1 P 0 kþ1 j j j s jj ¼ q j0 : 0 otherwise DDt F ¼ 2 ðDxÞ

ð41Þ ð42Þ ð43Þ

In the expressions above, the subscript denotes the node number, while the superscript refers to the time step number. Thus, ckj is the concentration of the substance at the node j at the time step k. Note that Eq. (42) is the finite difference approximation of Eq. (21), while Eq. (41) is the result of the reduction in mass sources due to the dissolution of contacts. Furthermore, Dt represents the time increment. Since the explicit scheme is employed here, Dt must be selected such as to ensure the numerical stability. Based on von Neumann stability analysis (see, for example, Gershenfeld [25]), the critical time increment is bounded by Dt 6

ðDxÞ2 =2   s f dbr 2q q t D þ ðDxÞ b 1 þ g  g 0 0 qs j0 bt

ð44Þ

In the numerical simulations described later in this section, the value of Dt was of the order of 103 s. In general, depending on the values of the key parameters, i.e. diffusion coefficient D and the constants governing the kinetics bt, br, bp, the dominant process will be that of diffusion (D), dissolution (bt, br) or precipitation (bp). Unfortunately, there is no accurate information available on the value of these coefficients. For example, the value of bt is within a broad range of 107–1011 m/s. The former value corresponds to the maximum rate of dissolution, which is 109.7 mol/cm2 s (Plummer et al. [26]). The latter one stems from an estimate that the dissolution rate in the direction normal to the surface of a contact undergoing pressure solution creep is of the order of 109 m/day (Dysthe et al. [27]). The same degree of ambiguity applies to the assessment of values of other coefficients. There are virtually no estimates for br, i.e. the coefficient governing the kinetics of dissolution in the cylindrical contacts. It is known that the value is significantly affected by the stress state; there is no reliable quantitative assessment though of this effect. Given the above limitations, the numerical simulations have focused on different qualitative scenarios associated with a range of values of bt and br. In the first set of simulations it was assumed that the concentration at the periphery of the intergranular contact is equal to zero. This corresponds to the case when the radius of the intergranular contact R is much smaller than that of the grain, so that the void space is said to contain pure water, i.e. concentration of the dissolved substance is negligible. The numerical analysis was carried out for some typical values of R = 100 nm, d = 4 nm, j0 = 2 nm, g0 = 0.1, qs = 2700 kg/m3,

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

for bt = 2 · 107 m/s, while the value of the coefficient br was taken as 106 m/s. The results are presented in Figs. 3–6. Fig. 3 shows the distribution of the normalized radius j/j0 of cylindrical contacts at different time intervals in the range up to 300 s. In this case , the dissolution starts at the periphery of the interface and the front gradually propagates towards the center. In general, the dissolution at the peripheries is very fast, while the subsequent propagation of the dissolved front is much slower. This is largely the result of the evolution of the concentration profile along the interface, as depicted in Fig. 4. The concentration of the dissolved substance (calcite) is zero at the periphery of the interface r = R, while that at the center becomes, almost instantaneously, equal to the solubility limit c. Since the kinetics of the process is proportional to the difference c  c, the rate of the dissolution at the peripheries is much faster than that closer to the center. This is evidenced in Fig. 4 which shows the history covering the period of 3000–12,000 s. Note that on the concentration profile c=c (the plot on the right-hand-side), the solid line depicts the

κ/κo 1.0 0.8 300 [s]

3 [s]

270 [s]

0.6

240 [s] 30 [s]

210 [s]

0.4

180 [s] 60 [s]

150 [s]

0.2

120 [s]

r[nm]

90 [s]

0

20

40

60

100

80

Fig. 3. Results of micromechanical analysis corresponding to bt = 2 · 107 m/s and br = 106 m/s; distribution of normalized radius j/j0 of cylindrical contacts at different time intervals.

qf = 1300 kg/m3, c ¼ 2 kg=m3 . For all simulations, the coefficient of diffusion was taken as 3 · 1015 m2/s. The first part of analysis was representative of a very fast kinetics of the dissolution process. It was carried out

a

297

b c/c

κ/κo

1.0

1.0 12000 [s]

0.8

0.8

9000 [s] 0.6

0.6

6000 [s]

0.4 0.2

0.2

r [nm] 20

40

60

t =3000 [s]

0.4

3000 [s]

80

r[nm]

100

0 7

t

20 r

40

60

80

100

6

Fig. 4. Results of micromechanical analysis corresponding to b = 2 · 10 m/s and b = 10 m/s; (a) distribution of j/j0 at different time intervals, (b) concentration profile c=c at t = 3000 s (dots – numerical solution, solid line – distribution corresponding to stationary state, i.e. when all contacts have been dissolved).

t=300 [s]

t = 3 [s]

κ/κo

κ/κo

100

50

t=12000 [s]

κ/κo

1.0 0.8

1.0 0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.0

r [nm] 50

100 100

50

0.0

0.2

r [nm] 50

100 100

50

0.0

r [nm] 50

100

Fig. 5. Results of micromechanical analysis corresponding to bt = 2 · 107 m/s and br = 106 m/s; propagation of the dissolution front at different time intervals.

298

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

ζ

c/c

0.1 [s] 0.2 [s] 0.3 [s] 0.4 [s] 1.0 [s] 10.0 [s] 100.0 [s]

1.0

0.8

0.8

0.6

0.6

0.4

0.4 0.2

0.2

t [s] 0

500

1000

1500

2000

2500

0.0

3000

Fig. 6. Results of micromechanical analysis corresponding to bt = 2 · 107 m/s and br = 106 m/s; evolution of f (dots – numerical solution, solid line – best fit approximation).

distribution corresponding to stationary state, i.e. when all contacts have been dissolved. It is apparent that after 3000 s the actual distribution of c=c (shown with dots) is very close to that corresponding to the stationary conditions, while the dissolution of contacts continues. Further illustration of the dissolution process is provided in Fig. 5. The figure shows a progressive propagation of the dissolution front at different time intervals; the region marked in dark indicates the area that remains unaffected by the dissolution process. Finally, Fig. 6 presents the evolution of the degradation parameter f, Eq. (25), which describes the reduction in the effective contact area. The dots represent the numerical results while the solid line is the best-fit approximation based on an exponential function f ¼ 1  a et=T 1  ð1  aÞ et=T 2

r [nm]

ð45Þ

where a, T1, T2 are constants. It should be noted that the best-fit corresponds to T1/T2 . 102, so that the process is defined by two significantly different time scales. This is consistent with the kinematics of dissolution process as depicted in Fig. 3. For a fixed geometry, T1 and T2 reflect the time scales associated with the rate of dissolution and diffusion, respectively. The next example corresponds to a much slower kinetics of the dissolution process; specifically bt = 2 · 1011 m/s and br = 4 · 109 m/s. All the other geometric/material parameters were the same as before. The numerical results are given in Figs. 7–9. Fig. 7 shows the evolution of the concentration of the dissolved substance along the interface. It is evident that the nature of the distribution is substantially different when compared to the previously considered scenario, Fig. 4. In the present case, the concentration at the center is much lower than the solubility limit. As a result, the rate of dissolution in the region adjacent to the center is similar to that at the peripheries of the interface. This is illustrated in Fig. 8 which shows the distribution of j/j0 at chosen time intervals. There is no rapid dissolution at the peripheries; the process is fairly uniform along the entire intergranular interface. In this case, the

20

40

60

80

100

Fig. 7. Results of micromechanical analysis corresponding to bt = 2 · 1011 m/s and br = 4 · 109 m/s; evolution of the concentration profile c=c along the interface.

κ/κo 5 [s] 30 [s] 60 [s] 90 [s] 120 [s] 150 [s] 180 [s] 210 [s] 240 [s] 270 [s]

1.0 0.8 0.6 0.4 0.2

600 [s]

0.0

20

40

60

80

100 r [nm]

Fig. 8. Results of micromechanical analysis corresponding to bt = 2 · 1011 m/s and br = 4 · 109 m/s: distribution of j/j0 at different time intervals.

gradient of the concentration of dissolved substance is small, which implies that the kinetics of the dissolution is now the key factor governing the phenomenon. This is in contrast to the former case, when the process was largely controlled by the diffusion which was much slower in rela-

ζ 1.0 0.8 0.6 0.4 0.2

t [s] 0.0

200

400

600

800

Fig. 9. Results of micromechanical analysis corresponding to bt = 2 · 1011 m/s and br = 4 · 109 m/s; evolution of f (dots – numerical solution, solid line – best fit approximation).

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

c/c

0.2 [s]

0.8

0.3 [s] 0.4 [s] 1.0 [s] 10.0 [s]

0.6 0.4 0.2 0.0

c/c

0.1 [s]

1.0

299

1.0 0.8

500 [s] 300 [s] 200 [s] 100 [s]

0.6 0.4 0.2

20

40

60

80

r [nm] 100

0.0

20

40

60

80

10 [s] r [nm] 100

Fig. 10. Results of micromechanical analysis accounting for precipitation phenomenon (bt = 2 · 1011 m/s, br = 4 · 109 m/s); evolution of the concentration c=c profile along the interface.

tion to the kinetics of the dissolution. Finally, Fig. 9 shows the evolution of the degradation parameter f. It is apparent that the nature of the distribution is now different and the function (45) is not as effective as before in terms of representing the phenomenon. The next set of simulations was focused on examining the influence of precipitation in the pore space on the kinetics of the dissolution within the intergranular interface. In this case, the representation (36) was employed that provided an explicit boundary condition at the periphery of the interface, i.e. at r = R. The simulations were carried using the same values of parameters b, i.e. bt = 2 · 1011 m/s and br = 4 · 109 m/s. Furthermore, it was assumed that bp ¼ bt , cp ¼ 0:3c and Rp = 2 · 107 m. Note that, for the above set of parameters, the value of the constant c appearing in Eq. (37) is close to unity (c = 0.9993). In this case, Eq. (36) can be simplified to  ocðRÞ 2pRd oc qs  qf ¼  p 3 D  þ s p bp ðcp  cðRÞÞ ð46Þ ot qR ðR Þ or r¼R which, in turn, implies the following incremental form k k k k cp  ckmþ1 ÞDt ckþ1 mþ1 ¼ cmþ1  A1 ð3cmþ1  4cm þ cm1 ÞDt þ A2 ð

ð47Þ where the constants A1, A2 are defined as A1 ¼

pRdD

; 3 2 ðRp Þ ðDxÞ

A2 ¼

ðqs  qf Þbp qs Rp

ð48Þ

It is evident here that, in the present case, the concentration of the dissolved substance at the boundary of the interface, r = R, evolves with time. The results of numerical simulations are provided in Figs. 10–12. Fig. 10 shows the evolution of normalized concentration of the dissolved substance c=c along the interface. At an early stage, soon after the voids are filled with water, there is an abrupt decrease in the concentration that is analogous to that depicted earlier (Fig. 7). Subsequently, as the time elapses, a uniform increase is observed, which stems from a progressive increase of the concentration at the periphery of the contact r = R. These two distinct phases of the process result in the dual time scale of the dissolution phenomenon. The first stage involves a rapid dissolution and is followed by a phase when the pro-

cess becomes much slower. This is clearly illustrated in Fig. 11, which shows the time history of the distribution of j/j0. It needs to be emphasized that the mechanism is similar to that depicted earlier, i.e. the dissolution is uniform along the entire interface. This is different from the first scenario that involved a rapid dissolution at the periphery and a slow propagation of the front towards the center (Figs. 3–5). The duality of the time scale has again implications on the evolution of the degradation parameter f, Eq. (25), whose time history is shown in Fig. 12. Again, the solid line is the best-fit approximation employing representation (45). In this case, the fit is almost perfect. 2.3. Pressure solution creep The last topic briefly addressed here is the pressure solution creep. In general, the dissolution of islands within the intergranular contact does not imply that the IPS process has been terminated. With time passing, there is a continuous dissolution of upper and lower boundaries of the intergranular interface. After a sufficiently long time, the distribution of concentration of the dissolved substance reaches a stationary state. This is the final stage of IPS that is associated with stationary transport of the dissolved material. The distribution of concentration in the stationary regime is described by the following set of equations:

κ/κo 1.0

10 [s] 30 [s] 60 [s] 90 [s] 120 [s] 150 [s] 180 [s] 210 [s]

0.8 0.6

1000 [s] 1150 [s] 1300 [s]

0.4 0.2

r [nm] 0.0

20

40

60

80

100

Fig. 11. Results of micromechanical analysis accounting for precipitation phenomenon (bt = 2 · 1011 m/s, br = 4 · 109 m/s); distribution of j/j0 at different time intervals.

300

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

ζ

hvðrÞi ¼

1.0

2 bt R2 qs

Z

R

rðc  cðrÞÞ dr

ð56Þ

0

Substituting now Eq. (52) and integrating, the following expression is obtained for the average stationary velocity of the propagating dissolution front

0.8 0.6

hvðrÞi ¼ 0.4

c  cp 2bt vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi

s p ffiffiffiffiffiffiffiffiffiffiffi  u Rq u I =D R A 3 uA3 0

t pffiffiffiffiffiffiffiffiffiffiffi  þ A4 A3 D I A =D R 1

0.2

t [s] 0.0

1000

2000

3000

4000

Fig. 12. Results of micromechanical analysis accounting for precipitation phenomenon (bt = 2 · 1011 m/s, br = 4 · 109 m/s); evolution of f (dots – numerical solution, solid line – best fit approximation).

    1 oc o2 c 2bt qs  qf þ D ¼ ðc  cÞ r or or2 d qs  2pRd oc qs  qf  p 3 D  þ s p bp ðcp  cðRÞÞ ¼ 0 qR ðR Þ or r¼R   oc ¼0 or 

ð49Þ ð50Þ ð51Þ

r¼0

Eq. (49) is a simplified form of Eq. (26) that describes the transport phenomenon within the interface. Note that (49) is obtained from (26) by neglecting the term involving the time derivative of c as well as imposing the constraint 1 = 1 (i.e., all islands dissolved) for which the last term within the brackets vanishes. Similarly, Eq. (50) is obtained from Eq. (36) by dropping the first term involving again the time derivative of c. The distribution of concentration of the dissolved substance, under stationary conditions, is described by the following function: rffiffiffiffiffi ! cp  c A3 cðrÞ ¼ c þ I 0 r rffiffiffiffiffi ! rffiffiffiffiffi ! D pffiffiffiffiffiffiffiffiffi A3 A3 I0 R þ A4 A3 D  I 1 R D D

ð57Þ

3

The above value defines directly the rate of convergence of _ i.e. the change in the distance between the centers grains D, of adjacent grains that remain in contact. Thus, c  cp hvðrÞi 2bt ffi D_ ¼ ¼ p s vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

p p ffiffiffiffiffiffiffiffiffiffiffi  R R Rq u u I0 A3 =D R uA3

pffiffiffiffiffiffiffiffiffiffiffi  þ A4 A3 t D I1 A3 =D R

ð58Þ

The variable D_ specifies the magnitude of the so-called pressure solution creep. In macroscopic sense, this is the rate of creep that is triggered by a stationary process of the dissolution along the upper/lower interphase boundaries between the neighboring grains. A quantitative assessment of the average rate of convergence can be performed using expressions (57) and (58). For the set of parameters employed in the last example of the previous section, the following values are obtained hvðrÞi ’ 4  1015 m=s ’ 0:35  109 m=day; D_  2  108 ½1=s  1:7  102 ½%=day Note that the value of the velocity of dissolution front is consistent with that suggested by Dysthe et al. [27]. In the reference quoted, the authors state that, for calcite, dissolution rates normal to the contact surface of a single contact undergoing pressure solution creep are of the order 109 [m/day]. 3. Phenomenological approach 3.1. Description of the kinetics of IPS

ð52Þ where I0 and I1 are hyperbolic Bessel functions and A3, A4 are constants defined by the following relations A3 ¼ 2

bt ðqs  qf Þ ; dqs

A4 ¼ 2pRd

qs p 2

bp ðR Þ ðqs  qf Þ

ð53Þ

The velocity of the dissolution of upper/lower interphase boundary is given by ð54Þ qs vðrÞ ¼ bt ðc  cðrÞÞ At the same time, the average velocity is defined through the relation Z R 2 rvðrÞ dr ð55Þ hvðrÞi ¼ 2 R 0 so that

As demonstrated in the preceding section, the intergranular pressure solution process typically consists of two distinct stages that incorporate, in general, two different time scales. An example here is a mechanism that involves a very fast dissolution of contacts at the periphery of the interface, followed by a much slower dissolution/diffusion in the region adjacent to the center of the interface. Other possible scenario is a rapid uniform dissolution along the horizontal phase transition boundaries, followed by a progressive increase in concentration at the periphery of the interface that is in contact with the pore space. The latter, subsequently, slows down the kinetics of the process. In order to emphasize the duality of the time scale, it is convenient to incorporate two independent scalar parameters f1, f2, which govern the degradation process on the

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

macroscale. Both these parameters can be interpreted as the change in the solid contact area of grains (Ac) in a representative volume of the material, i.e. f1,2  (Aco  Ac)/Aco, where Aco is the initial area at the onset of the reaction. Based on micromechanical considerations in the previous section, the evolution laws may be postulated in an exponential form Z 0

f1 ¼ fð1  eB1 t Þ;

0

f2 ¼ fð1  eB2 t Þ;

t0 ¼

gðcÞ dt

t

ð59Þ In the equations above, t 0 represents the so-called reaction time, whose evolution depends on the solubility c of solid (calcite) in a given fluid. In general, the function gðcÞ satisfies g  Æ0, 1æ, with g ! 0 for most non-aqueous solvents, such as oil (no IPS reaction). The parameters B1, B2 are material constants that define the respective time scales and f represents a stationary state associated with thermodynamic equilibrium. The latter is typically a function of confining pressure and the mathematical representation is chosen in such a way that f h0; 1Þ, which implies f1;2 h0; fi. The macroscopic manifestation of the continuing intergranular pressure solution is a progressive degradation of mechanical properties of the skeleton. Since both stages of the dissolution process contribute to the degradation, this process can be described in terms of variable f h0; fi defined as f ¼ af1 þ ð1  aÞf2 ) f_ ¼ g½B1 aðf  f1 Þ þ B2 ð1  aÞðf  f2 Þ

ð60Þ

where a  Æ0, 1æ is a weighting coefficient. Eq. (60) defines the evolution law that governs the kinetics of the IPS process. It is noted that the strain levels resulting from the degradation process are, in general, much lower than those experienced during the creep [9]. Thus, the creep phenomenon cannot be attributed to pressure solution alone, but is associated with other timedependent processes.

301

dominant. In this case, the description of the creep phenomenon incorporates the evolution parameter n = n(t), which is a measure of the progressive rearrangement of grains. It is postulated that the kinetics of the microstructural arrangement depends, once again, on the ‘distance’ from the stationary state, n. Thus, t ! 1 ) n !  n, which defines the configuration associated with internal equilibrium. Assuming, for convenience, that n  [0, 1], so that  n ¼ 1, the evolution law may be taken in a simple linear form, which is similar to that employed for the IPS process, i.e. n_ ¼ B3 ð1  nÞ ) n ¼ 1  eB3 t

ð61Þ

and B3 is a function of confining pressure, i.e. B3 = B3(rii), which is said to affect the kinetics of the process. In describing the macroscopic deformation, the total creep strain rates e_ cij can be expressed as the sum of the contributions from both phases of the process, i.e. e_ cij ¼ ~e_ ij þ hðt  t0 Þ e_ ij

ð62Þ

Here, ~e_ ij , e_ ij are the strain rates associated with the micro

structural rearrangement and the pressure solution, respectively; h is a Heaviside function and t0 is a threshold value for the onset of stationary conditions. In order to formulate the problem, restrict the considerations to the range of days, so that e_ cij  ~e_ ij , and assume that the macroscopic strain rate ~e_ ij , is co-axial with the stress tensor. Thus, ~e_ ij is an isotropic function of the type _ The general representation has the form ~e_ ij ¼ ~e_ ij ðrij ; nÞ. ~e_ ij ¼ a1 dij þ a2 rij þ a3 rik rkj

ð63Þ

where a’s are functions of n_ and the basic invariants of rij. Assume, for simplicity, that ~e_ ij is linear in rij and that ~e_ ij vanishes whenever rij vanishes. In this case, ~e_ ij ¼ a1 dij þ a2 rij ¼ ^a1 rkk dij þ ^a2 rij

ð64Þ

where a^1 ; ^a2 are functions of n_ alone. Take again a simple _ ^a2 ¼ a2 n, _ so that linear form ^a1 ¼ a1 n;

3.2. Time-dependent creep response

_ 1 rkk dij þ a2 rij Þ ¼ B3 ð1  nÞða1 rkk dij þ a2 rij Þ ~e_ ij ¼ nða

The creep response of chalk is assumed to involve two different physical mechanisms. The first one entails a viscous deformation that is associated with a progressive compaction of the solid skeleton. This mechanism is independent of IPS-induced degradation and is predominant in the early stages of creep, i.e. in the range of days. The second mechanism involves a stationary pressure solution creep as described in the preceding section. The latter is associated with grain convergence triggered by continuing dissolution at the intergranular grain boundaries. The time scale of the pressure solution creep is measured in geological range and while its impact in the early stages is negligible, it results in a significant deformation over a long period of time. In what follows, the focus is on the time interval that is in the range of days, i.e. when the viscous effects are pre-

where a1, a2 are constants and B3 = B3(rkk). Note that the relation above can be expressed in an alternative form ~e_ ij ¼ Cijkl rkl ;

Cijkl ¼ B3 ð1  nÞða1 dij dkl þ a2 dik djl Þ

ð65Þ

ð66Þ

Finally, B3 may be defined by employing, once again, a linear approximation B3 ¼ hhC 1  C 2 rii =fc0 ii

ð67Þ

Here, the symbol ÆÆÆææ stands for the Macauly bracket, C1, C2 are material parameters and fc0 is a normalizing constant identified with uniaxial strength in intact (i.e. oil-saturated) state. The representation (66) governs the evolution of strain rates associated with grain rearrangement. It is noted that this mechanism is present in both water and oil-saturated samples. The latter is consistent with existing experimental

302

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

evidence, which indicates that the creep strain attained in oil-saturated state are similar to those recorded in watersaturated samples [9]. The chemistry of the fluid though will affect the stationary creep rates e_ ij which are associated

with pressure solution. 3.3. Formulation of the problem Given the general description of the IPS process as well as the creep phenomenon, the formulation of the problem follows now the framework of chemo-plasticity analogous to that developed in the original reference [18]. The constitutive relation in the elastic range takes the form e ¼ ½C e r þ ~e ) e_ ¼ ½C e r_ þ

o e _ ½C rf þ ½Cr of

ð68Þ

where [Ce] is the elastic compliance and the operator [C] is defined by Eq. (66). The extension to elastoplastic regime is obtained by invoking the additivity of elastic and plastic strain rates. Thus e_ ¼ ½C e r_ þ

o e _ ½C rf þ ½Cr þ e_ p of

ð69Þ

In order to define the plastic strain rates, the functional form of the yield criterion f = 0 is assumed to be affected by the continuing IPS process, i.e. e_ p ¼ k_

f ðr; ep ; fÞ ¼ 0;

oQ or

ð70Þ

where Q = Q(r, f) represents the plastic potential function. Writing now the consistency condition as  T  T of of oQ of _ þ f¼0 ð71Þ k_ f_ ¼ r_ þ p or oe or of yields 1 k_ ¼ H

(

of or

T

) of _ r_ þ f ; of

 H ¼

of oep

T

oQ or

ð72Þ

Substituting Eq. (72) in Eqs. (70) and (69), the following constitutive relation is obtained after some transformations e_ ¼ ½Cr_ þ bf_ þ ½Cr where 1 ½C ¼ ½C  þ H e

(

 T ) oQ of ; or or

ð73Þ

b¼

o e 1 of oQ ½C r þ of H of or ð74Þ

Note that the first term on the right-hand side of Eq. (73) gives the elastoplastic strain rate due to the change in stress. The second term describes the strain rate resulting from degradation of properties triggered by IPS process, while the last term gives the strain rate due to the creep phenomenon. Clearly, in the absence of aqueous fluids, there is g ! 0 ) t 0 ! 0, Eq. (59) so that the second term, _ does not contribute to the deformation proinvolving f,

cess. Also, note that the elastoplastic response due to the imposed stress rate is assumed to be instantaneous. 3.4. Numerical example In order to provide an illustration of the framework presented above, an example is given here, which involves the numerical analysis of a long-term response of chalk. The simulations pertain to creep deformation of ‘Lixhie chalk’, for which the samples were obtained from a quarry near Liege (Belgium). The tests analysed here were conducted at the confining pressures of 2 MPa and 5 MPa, respectively, and involved water-saturated samples. The samples were tested a few hours after the saturation process was completed. The loading procedure consisted of applying a confining pressure first and, subsequently, increasing the deviatoric stress to a prescribed value using the loading rate typical of that employed in standard triaxial tests. During the creep stage, the stresses were kept constant and the corresponding strain history was recorded. In the numerical analysis, the conditions prior to creep can be established by following the deformation history associated with drained triaxial compression coupled with IPS-induced degradation. This involves the integration of the constitutive relation (73) and the details on the plasticity framework employed, together with specification of material functions, are given in Ref. [18]. Here, the primary focus is on the evolution of creep strain, which has been examined by implementing the constitutive relation (66). The constants appearing in the formulation were selected following the procedure described below. The description of creep employs four material parameters a1, a2, Eq. (65), as well as C1, C2, Eq. (67). For a given stress state, the constants a can be defined by specifying the maximum axial strain intensity ~e1 ! e for n ! 1 and the corresponding lateral strain ~e2 ¼ ~e3 ¼ #~e1 . Choosing, for the identification purpose, the test involving confinement of 5 MPa and the deviatoric stress intensity of 4 MPa, i.e. r = {5, 9, 5, 0, 0, 0} (MPa), and assuming e = 0.038; # = 0.05, yields, based on the representation (65) a1 ¼ 0:002725;

a2 ¼ 0:009975

The constants C1, C2 define the material function B3, Eq. (67), which in turn, controls the kinetics of the process, viz. Eq. (61). Assume that for the test conducted at rii = p0 there is n = n0 at t = t0. Similarly, let for rii = p1; n = n1 at t = t1. Substituting these conditions in representation (61), with B3 defined by Eq. (67), a set of two simultaneous equations is obtained that can be solved for C1, C2. For the simulations conducted here, both these constants were estimated using the conditions for the test at confinement of 5 MPa and deviatoric stress of 4 MPa (n = 0.95 at t = 20 h) as well as the test at confinement of 2 MPa and the deviatoric stress of 3.75 MPa (n = 0.5 at t = 80 h). The resulting values of C’s are: C 1 ¼ 0:140;

C 2 ¼ 0:107

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305

Note that for lateral deformation, only the results for deviatoric stress of 3.75 MPa are presented to make the figure more legible. Fig. 14 shows the simulation for the creep test conducted at the confinement of 5 MPa and the deviatoric stress of 4 MPa. Again, the present formulation can quite adequately depict the time history of both the axial as well as the lateral strain.

-ε1(%)

σ1 − σ3 = 3.6 MPa σ1 − σ3 = 3.75 MPa

3

σ1 − σ3 = 4.05 MPa Experimental

axial/lateral strain

2

4. Numerical example: finite element analysis of a water injection test 1

0

40

80

120

160

0

-ε3(%)

time (hrs)

-1

Fig. 13. Simulations of a series of creep tests at the confining pressure of 2 MPa.

- ε1(%)

The results of numerical simulations are provided in Figs. 13 and 14. Fig. 13 shows the deformation histories for a series of tests conducted at the confining pressure of 2 MPa and the deviatoric stress in the range of 3.6– 4.05 MPa. In general, for tests conducted at deviatoric stress intensities lower than 4 MPa, which is the strength of the fully-degraded material, the stable conditions are reached. On the other hand, for tests carried out at higher intensities, a spontaneous failure of the sample takes place. The results of numerical simulations are quite consistent with the experimental data, not only in qualitative but also the quantitative sense. In particular, both the axial and lateral deformation histories are predicted quite accurately.

5

Experimental

(σ1 − σ3 = 4.0 MPa)

Numerical

creep strain

4

3

2

1

-ε3(%)

303

0

-1

0

40

80

120

160

200

time (hrs)

Fig. 14. Evolution of creep strain for test at confining pressure of 5 MPa.

In this section, the formulation outlined earlier is applied in the context of a boundary value problem. The focus is on the short-term behaviour and the objective is to examine the interaction between the mechanical response of a system and the chemical degradation. The specific example selected here involves the numerical simulation of a water injection test as performed by Homand [19]. In general, a series of tests have been carried out at different initial conditions. The tests involved oil-saturated samples that were first subjected to axial compression at a given confining pressure. When the desired deviatoric stress intensity was reached, the water was injected from the bottom of the sample with a pressure in the range of 1–3 MPa. During the injection process, the axial and lateral strains were monitored at three different positions along the sample height. In all tests, the strain level remained very small before the water front reached the location of the respective strain gage. In the neighbourhood of the propagating front, a significant strain gradient was recorded that can be attributed to degradation triggered by the IPS process. The injection test represents a typical initial-boundaryvalue problem in the scale of the sample. In general, this is a coupled two-phase transient flow problem which involves a chemical degradation of the material. In this work, a simplified finite element analysis has been performed. As the injection pressure is much smaller than the applied stress level, hydromechanical coupling has been neglected and the emphasis was placed on the chemomechanical interaction. The sample was considered as a single-phase medium; the propagation of the water pffi front was described by a simple empirical law H d ¼ a t, where Hd is the current location of the front, measured from the bottom of the sample, and a is a parameter that controls the kinetics of the process. The above formula stems from classical theories of the water infiltration in the porous media (e.g. Kirkham and Powers [28], Dullien [29]). The numerical simulations have been carried out by incorporating the constitutive relation (73). Once again, the details on the plasticity framework employed are provided in Ref. [18]. Here, the formulation has been simplified by introducing a linear, rather than quadratic, form of the failure function. Thus, the evolution of plastic properties involved a decrease in cohesion coupled with the degradation of hardening characteristics. The values of material parameters chosen were consistent with those identified in the original reference. The only difference is

304

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305 5.00 26MPa

axial strain (%)

4.00 17MPa

3.00 Water injection

data

2.00

simulation

1.00 time (seconds)

0.00 0

2000

4000

6000

8000

10000

12000

Fig. 15. Evolution of axial strain at gauge #2 (mid-height of the sample) due to water injection; confining pressure 17 MPa, deviatoric stress of 9 MPa.

the value of the degradation constant G5 (i.e., G5 = 0.58) governing the evolution of hardening properties. The latter was changed (from G5 = 0.67) in view of the fact that the test simulated here involved much higher confining pressure than those examined in Ref. [18]. Typical results of finite element analysis are provided in Fig. 15. The simulations correspond to the test with confining pressure of 17 MPa and the deviatoric stress of 9 MPa. Here, the water was injected with a low pressure of 1 MPa and the strain gauges were positioned at three different levels with an equal distance of 18 mm apart. The comparison between the experimental and numerical results pertains to axial strain history measured at the strain gage #2 located at the mid-height of the sample. In general, there is a fairly good agreement in terms of the magnitude of strain induced by water injection. At the same time, there are some quantitative differences, particularly in terms of the kinetics of the strain evolution (i.e, the strain gradient is underestimated). It is noted though that the axial strain measured in the experiment is an averaged value that depends on the size of the strain gage. In contrast, the strain obtained in numerical simulation corresponds to a local value at a given material point. Thus, some discrepancies differences between these two values are to be expected. Overall, however, the proposed formulation appears to reproduce the chemo-mechanical coupling quite well. 5. Final remarks In this work, the framework developed in Ref. [18] has been enhanced by an extended micromechanical analysis. A multiscale approach has been employed to describe the physico-chemical processes at the level of intergranular contact and, later, to examine the chemo-mechanical coupling on the macroscale. The description of dissolution and diffusion of the dissolved substance has been obtained by incorporating molecular and volume averaging. The intergranular interface has been modelled as an evolving structure in which the contact between grains takes place

along a set of isolated islands surrounded by fluid. The mathematical formulation has been illustrated by extensive qualitative and quantitative studies focussed on examining the process of reduction in effective contact area triggered by dissolution. Two distinct scenarios have been identified which are governed by the kinetics of the dissolution and transport phenomena. In the first case, the islands at the periphery of the interface rapidly vanish and the dissolution front propagates towards the centre. The process is somewhat similar to that of propagation of rarefaction shock waves. The second scenario involves a uniform dissolution along the interface. In this case, the process is initially very fast and later it is slowed down due to the ongoing precipitation in the pore space. Thus, the evolution of contacts involves two distinct time scales. The macroscopic description of chemo-mechanical coupling has been formulated by invoking an evolution law that employed the notion of a dual time scale associated with the degradation process. In addition, the long-term creep response was addressed. Numerical simulations were conducted examining the creep characteristics of Lixhie chalk at different confining pressures and different deviatoric stress intensities. The results, both qualitative and quantitative, were consistent with the experimental evidence. Finally, the framework has been implemented in the finite element code and a numerical analysis of a water injection test was carried out. Again, the trends were similar to those observed experimentally, thus providing an additional degree of confidence in the predictive abilities of the framework. References [1] Schroeder Ch. Du coccolithe au re´servoir pe´trolier; approche phe´nome´nologique du comportement me´canique de la craie en vue de sa mode´lisation a` diffe´rentes e´chelles, Doctoral Thesis (in French), Belgium: University of Lie`ge; 2003. [2] Halleux L, Detiege C, Poot B, Schroeder C, Monjoie A, Debande G, Da Silva F. Mechanical behavior of chalks. In: Proceedings of the 3rd north sea chalk symposium, Copenhagen, Denmark; 1999. [3] Andersen MA. Petroleum research in north sea chalk, RF-Rogaland research, Stavanger, Norway; 1995. p. 179. [4] Monjoie A, Schroeder Ch, Da Silva F. Mechanical behaviour of chalks. In: Proceedings of the 3rd north sea chalk symposium, Stavanger, Norway; 1995. [5] Schroeder C, Bois AP, Maury V, Halle´ G. Water/chalk (or collapsible soil) interaction. Part II. Results of tests performed in laboratory on Lixhe chalk to calibrate water/chalk models. In: Proceedings of Eurock’98 rock mechanics in petroleum engineering, SPE Inc., Trondheim, Norway; 1998. [6] Homand S, Shao JF. Mechanical behaviour of a porous chalk and water/chalk interaction. Part I. Experimental study. Oil Gas Sci Technol 2000;55(6):591–8. [7] Brignoli M, Santarelli FJ, Righetti C. Capillary phenomena in impure chalk. In: Proceedings of Eurock’94 rock mechanics in petroleum engineering, Delft, Netherlans; 1994. [8] Risnes R, Flaageng O. Mechanical properties of chalk with emphasis on chalk–fluid interactions and micromechanical aspects. Oil Gas Sci Technol – Revue IFP 1999;54(6):751–8. [9] Xie SY. Contribution to modelling of mechanical behaviour of a porous rock, Doctoral thesis (in French), USTL, Lille, France; 2005.

D. Lydzba et al. / Computers and Geotechnics 34 (2007) 291–305 [10] Papamichos E, Brignoli M, Santarelli FJ. An experimental and theoretical study of partially saturated collapsible rocks. Mech Cohesive-Frict Mater 1997;2(3):251–78. [11] Collin F, Cui YJ, Schroeder C, Charlier R. Mechanical behaviour of Lixhe chalk partly saturated by oil and water: experiment and modelling. Int J Numer Anal Meth Geomech 2002;26(9):897–924. [12] Cristescu N. Damage and failure of viscoplastic rock-like materials. Int J Plasticity 1986;2(2):189–204. [13] Jin J, Cristescu N. An elastic viscoplastic model for transient creep of rock salt. Int J Plasticity 1998;14(1):85–107. [14] De Gennaro V, Delage P, Cui YJ, Schroeder C, Collin F. Time dependent behaviour of oil reservoir chalk: a multiphase approach. Soils Foundations 2003;43(4):131–48. [15] Hellemann R, Renders PJN, Gratier JP, Guiguet R. Experimental pressure solution compaction of chalk in aqueous solutions. Part 1. Deformation behaviour and chemistry, water–rock interactions. In: Hellmann R, Wood SA, editors. The Geochemical Society, special publication No. 7; 2002. p. 129–52. [16] Heggheim T, Madland MV, Risnes R, Austad T. A chemical induced enhanced weakening of chalk by seawater. J Petrol Sci Eng 2005;46(3):171–84. [17] Risnes R, Madland MV, Hole K, Kwabiah NK. Water weakening of chalk – mechanical effects of water–glycol mixture. J Petrol Sci Eng 2005;48(1):21–36. [18] Pietruszczak S, Lydzba D, Shao JF. Modelling of deformation response and chemo-mechanical coupling in chalk. Int J Numer Anal Meth Geomech 2006;30(10):997–1018. [19] Homand S. Comportement mecanique d’une craie tres poreuse avec prise en compte de l’effet de l’eau: de l’experience a la simulation, Doctoral thesis (in French), USTL, Lille, France; 2000.

305

[20] Raj R. Creep in polycrystalline aggregates by matter transport through a liquid phase. J Geophys Res 1982;87:4731–9. [21] Spiers CJ, Schutjens PM, Brzesowky RH, Peach CJ, Liezenberg JL, Zwart HJ. Experimental determination of constitutive parameters governing creep of rock salt by pressure solution. In: Knipe RJ, Rutter EH, editors. Deformation mechanisms, rheology and tectonics, Geological Society London, Special Publication 1990;54: 215–27. [22] Murdoch I. Fluid flow in porous media: continuum modelling based upon microscopic consideration. In: Modelling coupled phenomena in saturated porous materials. AMAS Lecture Notes 20, Institute of Fundamental Technological Research Polish Academy of Sciences, Warsaw, Poland; 2003, p. 485–540. [23] Cushmann JH. The physics of fluids in hierarchical porous media: Angstroms to Miles. Dordrecht: Kluwer Academic Press; 1997. [24] Jeschke AA, Dreybrodt W. Dissolution rates of minerals and their relation to surface morphology. Geochim Cosmochim Acta 2002;66(17):3055–62. [25] Gershenfeld N. The nature of mathematical modelling. Cambridge (MA): Cambridge University Press; 1999. [26] Plummer LN, Wigley TML, Parkhurst DL. The kinetics of calcite dissolution in CO2–water systems at 5 C to 60 C and 0.0–1.0 atm CO2. Am J Sci 1978;278:179–216. [27] Dysthe DK, Renard F, Feder JG, Jamtveit B, Meakin P, Jøssang T. High-resolution measurements of pressure solution creep. Phys Rev E 2003;68(1):011603. [28] Kirkham D, Powers WL. Advanced soil physics. New York: Wiley– Interscience; 1972. [29] Dullien FAL. Porous media: fluid transport and pore structure. New York: Academic Press; 1979.