Rock Plasticity from Microtomography and Upscaling

0. INTRODUCTION

The advantages of cost effective, reliable, sustainable, and environment friendly, geothermal energy has often been promoted over the recent years. The most important scientific problem relevant with geothermal energy is the permeability of reservoirs, which has attracted a great deal of attention. Another important problem is the strength of rocks, which is particularly crucial for enhanced geothermal system, as permeability must be increased for poor or low permeability reservoirs to improve production and/or injection performance. In this paper, we focus on rock plasticity from the microscopic point of view for the purpose to extend our understanding of rock strength at micro-scale using microtomographic images.

Microtomography makes it possible to detect internal structure of rocks on micro- to nano-scales thus opens a new way to quantify the relationship between the microstructures and their mechanical and transport properties. Numerical computations on microtomographic data have shown good agreement with experimental data for fluid flow and elastic properties (Zhu et al., 2012; Chen et al., 2010; Fredrich et al., 2006; Knackstedt et al., 2006; Arns et al., 2005, 2001). However, the extension of the methodology to plastic properties still is a poorly understood area.

Owing to the finite size of microtomographic images, the trend is that the higher the resolution, the smaller is the length scale of the sample. This leads to the immediate challenge to detect the microstructure of rocks across fine scales and connect the images to the macro-scale essential for describing the petroleum, geothermal or mining reservoir.

In this paper we use the workflow of upscaling of rock properties proposed previously (Regenauer-Lieb et al., 2013; Liu et al., 2012) and focus on rock strength and the upscaling of the yield stress. Two microtomographic datasets of carbonate, which show very different structures, are studied following the workflow. We compare the results with the previous study (Liu et al., 2012) and try to (1) find the relationships of the cohesion and the angle of friction as a function of the microstructures; (2) understand the characteristics of the upscaling of physical properties of heterogeneous rocks.

1. Workflow

Our work starts from the gray scale images of microtomography and involves three main components based on the binary data after segmentation (see Fig. 1). Three main components in the workflow are as following.

Figure  1.  Workflow for upscaling of rock properties from microtomography.

DownLoad: Full-Size Img PowerPoint

(1) Geometrical analyses (left column)—General parameters, such as porosity, connectivity, specific surface area, and the fractal dimension of pores are conducted. The size of geometrical representative volume element (RVE) is determined by using a stochastic analysis of the moving window method (Liu et al., 2009). The geometrical RVE is used in the computations of permeability to obtain reliable results with minimum computational cost.

(2) Mechanical analyses (middle column)—The size of mechanical RVE is determined by the computations of a series of models of different sizes with upper and lower bounds based on the theory of maximum and minimum entropy productions (Regenauer-Lieb et al., 2011, 2010; Terada et al., 2000). Then the mechanical RVE is used for computing elastic and plastic properties at micro-scale.

(3) Extracting critical exponents—A shrinking/expanding algorithm (Liu and Regenauer-Lieb, 2011) is used to create a series of derivative models with different porosity. Critical exponents, which describe the scale-independent exponential change of properties, can be extracted from these derivative models at micro-scale.

With any two critical exponents and/or the fractal dimension, scaling laws are defined (Stauffer and Aharony, 1994). Scaling laws, the percolation threshold, and some other parameters such as crossover length (refer to Liu and Regenauer-Lieb, 2011) are constraints used for the upscaling of properties from micro-scale to large scale.

The workflow is suitable for properties of fluid flow and mechanics. A detailed analysis of fluid flow and the upscaling of permeability is given in Liu et al. (2014). In this paper we focus on mechanical RVE, rock plasticity and critical exponent of yield stress in the workflow.

2. Datasets and characteristics

Two carbonate samples, named as CP10 and CP39, are studied in this paper. The scanned images are 2 000³ voxels, i.e., 2 000 slices and 2 000 by 2 000 pixels in each slice. After removing the blurred top and bottom slices and non-material regions (four corners of the square slice outside the round shape of the sample) of each slice, we use a 1 000³ volume for both of the samples for the following analysis. The resolution of CP10 is 1.85 μm and CP39 is 1.88 μm.

The Indicator Kriging (IK) segmentation method (Oh and Lidquist, 1999) is used to segment these samples. Segmentation thresholds are carefully tested and compared and the best-segmented model is selected based on the visual comparisons of images and the results of experimental porosity. IK segmentation tends to ignore isolated tiny pores and the porosity is generally smaller than that obtained from simple segmentation using one threshold value, and smaller than the experimental measurement as well. We define the results acceptable if the difference is less or equal to 2 units of porosity in the segmentation step. One typical greyscale image and the corresponding segmented binary image of each of the sample are shown in Fig. 2.

Figure  2.  Typical greyscale images and segmented binary images. (a) CP10 greyscale; (b) CP10 binary; (c) CP39 greyscale; (d) CP39 binary. Image size: 1 000×1 000 pixels.

DownLoad: Full-Size Img PowerPoint

Quantitative analyses of these binary images show that: (1) Porosities are 26.13% for CP10 and 17.87% for CP39. (2) There is a percolating pore network in both of the samples. The percolating pore network is isotropic for sample CP10 and anisotropic for CP39. (3) Except the percolating pore-structure, sample CP39 has enormous isolated pores occupying nearly 20% of the porosity; in contrast, CP10 has fewer isolated pores that only occupy 4% of the total porosity of the samples. (4) The fractal dimensions of the two samples are 2.36 (CP10) and 2.31 (CP39), which are close to each other. (5) The correlation length, which represents the average distance between any two voxels belonging to the same pore-structure, are very different, from 231 (CP10) to 407 μm (CP39).

3. mechanical analyses

3.1 Mechanical RVE

We use the thermodynamic upper and lower bound principles to determine the size of mechanical RVE. A series of models of different sizes under displacement and force loads are analysed to quantify whether the mechanical responses of the models converge, which indicates that they deliver homogenized values. The method was proposed by Terada et al. (2000) for elastic deformation and extended to more complex deformations by Regenauer-Lieb et al.(2011, 2010) based on the thermodynamic principle of maximum and minimum entropy production. Theoretically, upper bound (displacement boundary condition) and lower bound (force boundary condition) give different responses for models of small size; the difference diminishes as the model size increases (see Fig. 3). The size of mechanical RVE is determined from the converging solutions. We analyse elastic response only in this stage to determine mechanical RVE.

Figure  3.  Thermodynamic homogenization of upper and lower bounds of microstructures (from Regenauer-Lieb et al., 2010).

DownLoad: Full-Size Img PowerPoint

For each sample, a series of models are cropped from the 1 000³ volume. Models are cropped from a large one to a small one, step by step, and each smaller one is inside the nearby larger one. The cropped models have the porosity close to the biggest model, i.e., ±1 porosity unit.

Finite element computations are performed using Abaqus®, starting from smallest models to larger ones, and stopped when convergence can be seen. For small models, i.e., side-length L=40 and 60, hexahedral elements are used, in which one voxel equals to one element. For models with side-length larger than 60, tetrahedral elements are used to reduce the computing time. Tetrahedral meshes are created by using Avizo® package.

For each of these models, boundary surfaces of x=x_min, y=y_min, z=z_min are normal displacement constrained, boundaries x=x_max, y=y_max are free. Displacement boundary condition (BC) or pressure BC are loaded on the surface of z=z_max. We limit the displacement BC and pressure BC to be a low value to ensure only small deformation occurs and deformation is in elastic regime. In the binary models, pores are void and matrix is solid. We give the input parameters of solid of elastic modulus of 50 GPa and Poisson's ratio of 0.20 referring to the mechanical parameters of rocks (Schön, 2011).

Figure 4 shows the homogenization of elastic modulus and shear modulus of the two samples. Mechanical RVE size can be determined when the curves of upper and lower bounds are convergent. For sample CP10, the mechanical RVE size is 200³. For sample CP39, it is 300³. These mechanical RVEs are used to analyse the plastic response at micro-scale in the next step. The RVE models are illustrated in Fig. 5.

Figure  4.  Elastic modulus and shear modulus of upper and lower bounds dissipation of different volume sizes. (a) Sample CP10; (b) sample CP39.

DownLoad: Full-Size Img PowerPoint

Figure  5.  Microstructure and mesh of mechanical RVEs of CP10 L=200 (a) and CP39 L=300 (b).

DownLoad: Full-Size Img PowerPoint

3.2 Plastic Strength of Microstructures

Drucker-Prager plasticity and only elastic perfectly plastic behaviour are considered as matrix behaviour in our analyses of plastic strength on mechanical RVEs. The linear Drucker- Prager relationship σ_y=σ_ntanβ+c, where σ_y is the yield stress, σ_n is pressure (or mean stress), c is cohesion, and β is the angle of internal friction, is used to derive cohesion and the angle of friction from pressure dependent yield stress. Thus two cases of different pressures are computed. Case 1 is using displacement BC on z_max surface and normal displacement constraints on x_min, y_min, and z_min surfaces. Case 2 is using the same loading on z_max surface and normal displacement constraints on all other surfaces.

Two more parameters are necessary for Drucker-Prager plasticity as input: the angle of internal friction of the solid matrix, which is specified as 40° in this study; and cohesion, which is specified as 30 MPa.

Results show that: (1) For Case 1, stress-strain relationship, which represented by the minimum principal stress and strain averaged over all elements, performs as elastic perfect plasticity; while Case 2 performs as plastic hardening behaviour (not pictured here, refer to Liu et al., 2012). (2) Equivalent plastic strain concentrates in areas of weak links, see Fig. 6 a typical slice of our model. It shows a reasonable result that failure starts from the fragile connections between large pores. (3) For Case 1, von Mises stresses and pressure form a linear relationship before yielding, they are constant after yielding. (4) For Case 2, these two variables roughly form two linear relationships, before and after yielding, respectively. (5) Points of Case 1 and Case 2 after yielding are linearly distributed and from the line we can fit the cohesion and the angle of friction, see Fig. 7. (6) The fitted cohesion is close to each other for sample CP10 and CP39, around 21 MPa which is 70% of the input solid parameter. (7) The angle of friction is 7.6° for sample CP10 and 16.4° for sample CP39, which are both much lower than the input solid parameter of 30°.

Figure  6.  Equivalent plastic strain of two typical slices. Left. Case 1; right. Case 2.

DownLoad: Full-Size Img PowerPoint

Figure  7.  Relationship of Mises stress and pressure, and the fitting of cohesion and friction angle from points after yielding of two cases.

DownLoad: Full-Size Img PowerPoint

3.3 Upscaling of Properties

Only the critical exponent of yield stress is analysed in this paper. A shrinking/expanding algorithm (Liu and Regenauer-Lieb, 2011) is used to create a series of derivative models with different volume fractions. The percolation threshold p_c, is here defined as the lowest volume fraction for an unbroken structure connecting at least two opposite outside boundaries, and it can be detected from these derivative models. Then the derivative models close to the percolation threshold are used to simulate the deformation and yield phenomenon. With a series of results of yield stress of models, it is possible to fit the critical exponent of yield stress T_y in the form of σ_y=(p–p_c)^Ty, when the volume fraction T_y is a scale-independent parameter describing the change of yield stress along with volume fraction for all scales.

The original models have the solid volume fractions of 73.87% (CP10) and 82.23% (CP39), respectively. By deflating 10 steps, the model of CP10 is percolating only in the y direction and we define its volume fraction 7.57% as the percolation threshold. Five derivative models that are close to the percolation threshold and are percolating in 3 directions are analysed to determine their yield stresses. The derivative model after 14 deflating steps of the sample CP39 is only percolating in y direction. Thus the percolation threshold is the volume fraction of this model, which is 3.23%. As the models of 10th to 13th shrinking (deflating) are not percolating in all 3 directions, they cannot be used for the computations for fitting the critical exponent of yield stress. The derivative models used for computations and their results of yield stress are listed in Table 1.

Table  1.  Derivative models and the results of sample CP10 and CP39 for extracting critical exponent of yield stress

CP10 pc=7.57%	Deflation	8	7.5	7	6.5	6
	(p)	0.127 9	0.135 7	0.165 3	0.175 3	0.212 5
	(p–pc)	0.052 2	0.060 0	0.089 6	0.099 6	0.136 8
	Nodes	972 771	1 069 938	1 391 835	1 730 534	2 103 375
	σy	0.035 9	0.042 8	0.071 1	0.116 2	0.369
CP39 pc=3.23%	Deflation	9	8	7	6	5
	(p)	0.104 3	0.134 3	0.173 1	0.223 0	0.286 4
	(p–pc)	0.072	0.102	0.140 8	0.190 7	0.254 1
	Nodes	3 292 071	4 379 419	5 721 366	7 465 946	9 825 960
	σy	0.71	0.88	1.08	1.35	Unable

 | Show Table

DownLoad: CSV

As these derivative models are all close to the percolation threshold, generally links are very weak. We found problems with a tetrahedral element mesh as this may break links and change the connectivity of the models. Therefore we use hexahedral element meshes for all these derivative models. The disadvantage is that the number of elements is very large and this causes high computing cost. We found that the parallel capability of Abaqus® Explicit is limited to eight million element nodes. This was enough for most calculations and we only had one model with more than nine million nodes, which is unable to perform the computation (see the bottom of the last column of Table 1).

By fitting the yield stress with respect to the difference of volume fraction and the percolation threshold |p–p_c| in log-log plot (Fig. 8), we obtained a critical exponent of yield stress of 2.3 from sample CP10, and it is 0.66 from sample CP39. Previously, we have obtained a value of 1.33 (Liu et al., 2012) from another carbonate sample. These three values are covering an extremely wide span. Laboratory experiments with a regular structure (Benguigui et al., 1987) reported a value of ~2.5 and Sieradzki and Li (1986) obtained a value of 1.7. In their experiments, a metal plate was punched to create more and more holes thereby approaching the percolation threshold. The distribution of holes was arbitrary or aligned. Our 3D models are strongly heterogeneous comparing with their 2D laboratory models. Considering the great difference of the models, we find it plausible that the values of critical exponent of the natural rocks are different from the experiment results.

Figure  8.  Stress-strain relationships of five derivative models close to the percolation threshold (left) and the fitting of the critical exponent of yield stress (right).

DownLoad: Full-Size Img PowerPoint

4. DISCUSSION AND Conclusions

We have established and tested a workflow of upscaling of petrophysical properties using microtomography and percolation theory. The workflow is valid for both properties of fluid flow and mechanical response. In this paper we focused on plastic response using two carbonate samples with different structures in order to study the relationship between mechanical strength and microstructures for geothermal and petroleum exploitation.

The size of mechanical RVE is determined by using upper and lower bound asymptotinc homogenisation finite element computations. Our examples show good convergence of upper and lower bounds when the volume is increasing and the mechanical RVEs can be determined. The size of mechanical RVE is not directly related to the porosity but relies on the matrix structure. Different from the feature of permeability, mechanical properties are not directly related to the pore structure. But we see the RVE sizes of these two samples follow the trend of the correlation lengths of the samples.

Two cases of different pressure are analysed for Drucker-Prager plasticity and macrosopic cohesion and the angle of friction are derived. For both cases, we can see high plastic strain in weak links inside the microstructures. Von Mises stress and pressure of the two cases after yielding form a straight line, which is inconsistent with the linear Drucker-Prager relationship. We have found that the cohesion for these two samples is close to each other while the angle of friction is quite different. Sample CP10 with higher porosity has much lower friction angle.

The critical exponents of yield stress of the two samples are extracted by fitting the yield stress of derivative models and the difference of volume fraction and the percolation threshold in log-log plots. Our results span 0.66 to 2.3. This variability is in direct contradiction with classical percolation theory of random media where it should be a constant. Considering, however, that the specific structures are far from being random models this conflict with classical percolation theory can be expected. It is interesting to detect the relationship of the critical exponent of yield stress with respect to specific structures. Progress on this problem will not only significantly contribute to percolation theory for structured media, but also greatly improve the understanding of rock strength at micro-scale and its link to large scale. Our results of the pilot study therefore prompt for additional work in this exciting research area.