Rock Plasticity from Microtomography and Upscaling - Springer Link

Journal of Earth Science, Vol. 26, No. 1, p. 053–059, February 2015 Printed in China DOI: 10.1007/s12583-015-0520-4

ISSN 1674-487X

Rock Plasticity from Microtomography and Upscaling Jie Liu*1, 2, Reem Freij-Ayoub3, Klaus Regenauer-Lieb4, 3, 2 1. School of Earth Science and Geological Engineering, Sun Yat-Sen University, Guangzhou 510275, China 2. School of Earth and Environment, the University of Western Australia, Perth WA 6000, Australia 3. Earth Science and Resource Engineering, CSIRO, Kensington WA 6151, Australia 4. School of Petroleum Engineering, University of New South Wales, Sydney NSW 2052, Australia ABSTRACT: We present a workflow for upscaling of rock properties using microtomography and percolation theory. In this paper we focus on a pilot study for assessing the plastic strength of rocks from a digital rock image. Firstly, we determine the size of mechanical representative volume element (RVE) by using upper/lower bound dissipation computations in accordance with thermodynamics. Then the mechanical RVE is used to simulate the rock failure at micro-scale using FEM. Two cases of different pressures of linear Drucker-Prager plasticity of rocks are computed to compute the macroscopic cohesion and the angle of internal friction of the rock. We also detect the critical exponents of yield stress for scaling laws from a series of derivative models that are created by a shrinking/expanding algorithm. We use microtomographic data sets of two carbonate samples and compare the results with previous results. The results show that natural rock samples with irregular structures may have the critical exponent of yield stress different from random models. This unexpected result could have significant ramifications for assessing the stability of solid materials with internal structure. Therefore our pilot study needs to be extended to investigate the scaling laws of strength of many more natural rocks with irregular microstructure. KEY WORDS: plastic strength, rock, microtomography, percolation theory, 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 *Corresponding author: [email protected] © China University of Geosciences and Springer-Verlag Berlin Heidelberg 2015 Manuscript received August 19, 2014. Manuscript accepted October 27, 2014.

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. (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

Liu, J., Freij-Ayoub, R., Regenauer-Lieb, K., 2015. Rock Plasticity from Microtomography and Upscaling. Journal of Earth Science, 26(1): 53–59. doi:10.1007/s12583-015-0520-4

Jie Liu, Reem Freij-Ayoub and Klaus Regenauer-Lieb

54 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 0003 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 0003 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 bestsegmented 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. 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.

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

Rock Plasticity from Microtomography and Upscaling

55

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.

Figure 3. Thermodynamic homogenization of upper and lower bounds of microstructures (from Regenauer-Lieb et al., 2010). 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. For each sample, a series of models are cropped from the 1 0003 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.,


56

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 2003. For sample CP39, it is 3003. 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.

side-length L=40 and 60, hexahedral elements are used, in which one voxel equals to one element. For models with sidelength 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=xmin, y=ymin, z=zmin are normal displacement constrained, boundaries x=xmax, y=ymax are free. Displacement boundary condition (BC) or pressure BC are loaded on the surface of z=zmax. 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).

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 DruckerPrager 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 20

Shear modulus

Elastic modulus

40 30 20 10

Displacement BC Pressure BC 0

100 200 Model size ( L in voxel)

(a)



300

20

Shear modulus

Elastic modulus

10 5

300

45 35 25 15

15



300

(b)

15 10 5



300

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

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

Rock Plasticity from Microtomography and Upscaling

57 which are both much lower than the input solid parameter of 30°.

friction from pressure dependent yield stress. Thus two cases of different pressures are computed. Case 1 is using displacement BC on zmax surface and normal displacement constraints on xmin, ymin, and zmin surfaces. Case 2 is using the same loading on zmax 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,

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 pc, 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 Ty in the form of σy=(p–pc)Ty, when the volume fraction Ty is a scaleindependent 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

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

σ y=0.133 σ n+21.04

20

0

0

10

20 30 Pressure (MPa)

(b)

30

Case 1-not yield Case 1-yield Case 2-not yield Case 2-yield

10

Mises stress (MPa)

Mises stress (MPa)

40 (a)

30

40

50

σ y=0.294 9 σ n+21.3

20

Case 1-not yield Case 1-yield Case 2-not yield Case 2-yield

10 0

0

10

20 30 Pressure (MPa)

40

50

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


58

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.

(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. 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–pc| 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

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

Table 1 Derivative models and the results of sample CP10 and CP39 for extracting critical exponent of yield stress CP10 pc=7.57%

CP39 pc=3.23%

Deflation (p) (p–pc) Nodes σy Deflation (p) (p–pc) Nodes σy

8 0.127 9 0.052 2 972 771 0.035 9 9 0.104 3 0.072 3 292 071 0.71

7.5 0.135 7 0.060 0 1 069 938 0.042 8 8 0.134 3 0.102 4 379 419 0.88

5.0E-06 p - p c=0.072 0 p - p c=0.102 0 p - p c=0.140 8 p - p c=0.190 7

1.0E-05 Strain

1.5E-05

Yeild stress

0.1

1.2 Stress (MPa)

p - p c=0.052 2 p - p c=0.060 0 p - p c=0.080 0 p - p c=0.099 6 p - p c=0.136 8

0.2

1.5

0.2

(a) 2.0E-05

0.02 0.03

| p - p c|

0.3

3 σ y=3.962 6( p-p c) 0.656 5

0.9 0.6 0.3

6 0.212 5 0.136 8 2 103 375 0.369 5 0.286 4 0.254 1 9 825 960 Unable

σ y=26.629( p-p c) 2.299 4

0.3

0.0 0.0E+00

6.5 0.175 3 0.099 6 1 730 534 0.116 2 6 0.223 0 0.190 7 7 465 946 1.35

Yeild stress

Stress (MPa)

0.4

7 0.165 3 0.089 6 1 391 835 0.071 1 7 0.173 1 0.140 8 5 721 366 1.08

(b)

0.0 0.3 0.0E+00 1.0E-05 2.0E-05 3.0E-05 4.0E-05 5.0E-05 6.0E-05 0.04 Strain

| p-p c|

0.4

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

Rock Plasticity from Microtomography and Upscaling 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 DruckerPrager 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 loglog 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. ACKNOWLEDGMENTS We are grateful to Petrobras’ financial support of this research and the permission of this publication. We are also thankful to iVEC for the technical support and access to its visualisation facilities and supercomputers. Claudio Delle Piane and Valeriya Shulakova in CSIRO ESRE contributed to data acquisition and image pre-processing of microtomography. We are also thankful to Ben Clennell for the support and suggestive discussions while implementing the project. REFERENCES CITED Arns, C. H., Knackstedt, M. A., Martys, N. S., 2005. CrossProperty Correlation and Permeability Estimation in Sandstone. Physical Review E, 72(4): 046304. doi:10.1103/PhysRevE.72.046304 Arns, C. H., Knackstedt, M. A., Pinczewski, W. V., et al., 2001. Accurate Estimation of Transport Properties from Microtomographic Images. Geophysical Research Letters, 28(17): 3361–3364 Benguigui, L., Ron, P., Bergman, D. J., 1987. Strain and Stress at the Fracture of Percolative Media. J. Physique, 48(9): 1547– 1551 Chen, C., Packman, A. I., Zhang, D. X., et al., 2010. A Multi-Scale Investigation of Interfacial Transport, Pore Fluid Flow, and Fine Particle Deposition in a Sediment Bed. Water Resources Research, 46: W11560. doi:10.1029/2009WR009018 Fredrich, J. T., DiGiovanni, A. A., Noble, D. R., 2006. Predicting

59 Macroscopic Transport Properties Using Microscopic Image Data. Journal of Geophysical Research, 111(B3): B03201. doi:10.1029/2005JB003774 Knackstedt, M. A., Arns, C. H., Saadatfar, M., et al., 2006. Elastic and Transport Properties of Cellular Solids Derived from Three-Dimensional Tomographic Images. Proc. R. Soc. A, 462(2073): 2833–2862 Liu, J., Pereira, G. G., Regenauer-Lieb, K., 2014. From Characterisation of Pore-Structures to Simulations of Pore-Scale Fluid Flow and the Upscaling of Permeability Using Microtomography: A Case Study of Heterogeneous Carbonates. Journal of Geochemical Exploration, 144(Part A): 84–96. doi:10.1016/j.gexplo.2014.01.021 Liu, J., Regenauer-Lieb, K., 2011. Application of Percolation Theory to Microtomography of Structured Media: Percolation Threshold, Critical Exponents and Upscaling. Physical Review E, 83(1): 016106. doi: 10.1103/ PhysRevE.83.016106 Liu, J., Freij-Ayoub, R., Regenauer-Lieb, K., 2012. Determination of the Plastic Strength of Carbonates from Microtomography and Upscaling Using Percolation Theory. In: Eberhardsteiner, J., et al., eds., Proceeding of ECCOMAS 2012. September 10–14, 2012, Vienna Liu, J., Regenauer-Lieb, K., Hines, C., et al., 2009. Improved Estimates of Percolation and Anisotropic Permeability from 3-D X-Ray Microtomographic Model Using Stochastic Analyses and Visualization. Geochem., Geophys., Geosyst.,10: Q05010. doi:10.1029/2008GC002358 Oh, W., Lidquist, W. B., 1999. Image Thresholding by Indicator Kriging. IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(7): 509–602 Regenauer-Lieb, K., Karrech, A., Chua, H. T., et al., 2011. NonEquilibrium Thermodynamics for Multi-Scale THMC Coupling. GeoProc 2011, 6–9 July, Perth Regenauer-Lieb, K., Karrech, A., Chua, H., et al., 2010. TimeDependent, Irreversible Entropy Production and Geodynamics. Philosophical Transactions of the Royal Society London, A, 368(1910): 285–300 Regenauer-Lieb, K., Veveakis, M., Poulet, T., et al., 2013. Multiscale Coupling and Multiphysics Approaches in Earth Science: Theory. Journal of Coupled Systems and Multiscale Dynamics, 1(1): 49–73 Schön, J. H., 2011. Physical Properties of Rocks: A Workbook. Handbook of Petroleum Exploration and Production, 8: 1– 494 Sieradzki, K., Li, R., 1986. Fracture Behavior of a Solid with Random Porosity. Physical Review Letters, 56(23): 2509–2512 Stauffer, D., Aharony, A., 1994. Introduction to Percolation Theory (Second Ed.). Taylor & Francis Ltd., London Terada, K., Hori, M., Kyoya, T., et al., 2000. Simulation of the Multi-Scale Convergence in Computational Homogenization Approaches. International Journal of Solid and Structures, 37: 2285–2311 Zhu, B. J., Cheng, H. H., Qiao, Y. C., et al., 2012. Porosity and Permeability Evolution and Evaluation in Anisotropic Porosity Multiscale-Multiphase-Multicomponent Structure. Chin. Sci. Bull., 57(4): 320327. doi:10.1007/s11434011-4874-4