# Computational upscaling of Drucker-Prager plasticity from micro-CT images of synthetic porous rock

Computational upscaling of Drucker-Prager plasticity from micro-CT images of synthetic porous rock Abstract Quantifying rock physical properties is essential for the mining and petroleum industry. Microtomography provides a new way to quantify the relationship between the microstructure and the mechanical and transport properties of a rock. Studies reporting the use microtomographic images to derive permeability and elastic moduli of rocks are common; only rare studies were devoted to yield and failure parameters using this technique. In this study, we simulate the macroscale plastic properties of a synthetic sandstone sample made of calcite-cemented quartz grains using the microscale information obtained from microtomography. The computations rely on the concept of representative volume elements (RVEs). The mechanical RVE is determined using the upper and lower bounds of finite-element computations for elasticity. We present computational upscaling methods from microphysical processes to extract the plasticity parameters of the RVE and compare results to experimental data. The yield stress, cohesion and internal friction angle of the matrix (solid part) of the rock were obtained with reasonable accuracy. Computations of plasticity of a series of models of different volume-sizes showed almost overlapping stress-strain curves, suggesting that the mechanical RVE determined by elastic computations is also valid for plastic yielding. Furthermore, a series of models were created by self-similarly inflating/deflating the porous models, that is keeping a similar structure while achieving different porosity values. The analysis of these models showed that yield stress, cohesion and internal friction angle linearly decrease with increasing porosity in the porosity range between 8 and 28 per cent. The internal friction angle decreases the most significantly, while cohesion remains stable. Microstructure, Plasticity, diffusion and creep, Numerical modelling 1 INTRODUCTION Understanding upscaling of microscale physical processes to the scale of a core is a standard problem for the hydrocarbon, geothermal, water and mining industries. Elasticity, plasticity, permeability and other properties of rocks are related to the structure of pores and cracks as well as their connectivity. Microtomography enables the determination of the internal structure of rocks at the micro- to nanoscales and allows quantification of the relationship between microstructure and mechanical and transport properties through computational forward simulations of the underlying processes (Spanne et al.1994; Stock 2008; Blunt et al.2013). Numerical computations on microtomographic data have shown good agreement with experimental data for fluid flow and elastic properties (Arns et al.2001, 2005; Fredrich et al.2006; Knackstedt et al.2006; Chen et al.2010). However, the extension of the methodology to plastic properties still represents a poorly understood area. One aspect is the difficulty of defining a representative volume element (RVE) for dissipative processes such as plastic deformation from microtomographic data. This paper aims to address this problem by suggesting an energy-based up-scaling procedure (Veveakis & Regenauer-Lieb 2015). The study of rock properties based on microtomographic images has been developing for two decades since the middle of 1990s. Permeability is the first rock property that has been assessed through computational homogenisation method based on microtomographic structure (Spanne et al.1994). Arns et al. (2001) followed on from this early pioneering work and performed a benchmark analysis of transport properties derived from microtomographic images of Fontainebleau sandstone. Porosity and permeability, and their scale dependency have also been analysed for other porous rocks (Fredrich et al.2006; Chen et al.2010). The power of the new method led to many other works. Elastic properties derived from microtomographic images were, for instance, studied by Arns et al. (2002). Roberts & Garboczi (2002) computed the linear elastic properties of random porous materials with a wide variety of microstructures. Almqvist et al. (2011) applied differential effective medium and X-ray microtomography to calculate the elastic properties of porous and anisotropic rock aggregates; Shulakova et al. (2013) also computed the elastic properties of natural sandstones from microtomographic images; Knackstedt et al. (2006) analysed the correlation between the up-scaled diffusivity, elasticity, permeability and conductivity of industrial cellular foams. Rock deformation in nature occurs far beyond elasticity and no robust method for upscaling post-elastic rock behaviour exists to date. In this work, we discuss the derivation of plastic properties of rocks as a natural extension of the methods introduced for computational homogenisation of elastic properties. In doing so, we first define the yield envelope and the plastic potential of the solid skeleton of the porous rock. Laboratory work on porous silicate rocks show that Drucker-Prager and Mohr-Coulomb adequately describe plastic yield. Because of its smooth yield surface, the Drucker-Prager criterion appears better suited than Mohr-Coulomb criterion for numerical simulations. Liu et al. (2012, 2015) simulated the Drucker-Prager plastic response of carbonate samples at the microscale using microtomography. However, due to the extremely strong heterogeneity of carbonates involving at least two different scales of pores, those pilot studies remain inconclusive. The available micro-CT scans, were too small to convincingly test the validity of RVE and hence the results on the relationship between porosity and plastic properties could not be fully understood. In this study, we therefore analyse a well-defined synthetic sandstone sample and derive the yield stress, cohesion and friction angle from 3-D X-ray images of the microstructure. The computations are based on the RVE defined from elastic computations, and a subsequent test of RVE for computational homogenisation of plastic properties. In a second step, the computed properties are compared with experimental measurements. Finally, the dependence of plastic properties on porosity is deduced from a series of models with similar structure but different porosity values. 2 ROCK MATERIAL AND EXPERIMENTAL DATA 2.1 Microtomographic data In this study, we used a CIPS (Calcite In situ Precipitation System, a patent technique) synthetic sandstone sample composed of pure quartz grains (diameter 0.15–0.35 mm) and cemented by calcite. The data collected for this CIPS is freely available from CSIRO data repository (Yang et al.2012). The data includes 3-D X-ray microtomography images and mineral phase distributions. A cylindrical plug, 38 mm in diameter, was drilled from the man-made block. The gas permeability of the sample was measured and compared well with the numerically simulated permeability from the point of view of percolation theory (Stauffer & Aharony 1994; Liu & Regenauer-Lieb 2011). A small sample, 5 mm in diameter, was prepared and X-ray CT-scanned with a resolution of 2.5 μm. After segmentation, a porosity of 24.4 per cent was derived from the images, which is lower than the experimental measurement of Helium porosity of 28.5 per cent. This discrepancy is most likely due resolution limits preventing imaging smaller size pores below 2.5 μm. Detailed information about image processing and geometrical characteristics of this sample can be found in Liu et al. (2009). The processed 3-D image of a subvolume of the scanned sample is shown in Fig. 1. For statistical analysis we use a Cartesian reference frame and investigate percolation in three orthogonal directions thus show cropped parallelepipeds out of the cylindrical sample. Figure 1. View largeDownload slide Visualization of microtomographic data of the synthetic sandstone sample analysed in this paper. (a) 3-D microstructure, in the top slice the binary image shows the solid in black and the pores in white (volume size is 600 × 600 × 400 voxels); (b) the coordinate frame used. Figure 1. View largeDownload slide Visualization of microtomographic data of the synthetic sandstone sample analysed in this paper. (a) 3-D microstructure, in the top slice the binary image shows the solid in black and the pores in white (volume size is 600 × 600 × 400 voxels); (b) the coordinate frame used. 2.2 Rock deformation data A cylindrical CIPS plug, 25.4 mm in diameter and 47 mm in length, was subjected to an axial loading while maintaining uniaxial strain conditions. The experiment is conducted in a conventional (axisymmetric) triaxial stress vessel, using oil as the confining medium acting mechanically on the rock through a sealing flexible membrane. It consists in applying an increasing axial displacement on the specimen's end face at a constant rate (1.8 × 10−5 s−1 in this instance) while the oil pressure is servo-controlled to cancel the radial strain induced by Poisson's effect. Therefore, the rock specimen is nominally subjected to a uniaxial strain field, similar to a conventional oedometric tests on soils. Therefore, the strain field within the specimen is nominally uniaxial at mid-height of the specimen, that is away from the frictional boundary conditions induced by the loading platens. The specimen is tested in dry conditions (no pore fluid) so that there is no internal pore pressure within the rock, and therefore the external stresses applied to the specimen are equal to the effective stresses. Fig. 2 depicts the evolution of the differential, radial and mean effective stresses as a function of the applied axial strain during the experiment. For sake of simplicity the nominally null radial strain is not represented in this figure. Figure 2. View largeDownload slide Evolution of the differential stress (black curve), effective mean stress (dotted curve), and radial stress (dashed curve) with increasing axial strain during the uniaxial strain experiment. The radial strain remains nominally zero. Figure 2. View largeDownload slide Evolution of the differential stress (black curve), effective mean stress (dotted curve), and radial stress (dashed curve) with increasing axial strain during the uniaxial strain experiment. The radial strain remains nominally zero. Note the stresses in Fig. 2 are not starting from zero. This is because of a technical limitation of the equipment used to conduct this uniaxial strain test in a conventional triaxial stress vessel: we have to apply a small confining pressure (1 MPa) prior to applying the displacement-controlled axial stress and servo-controlling the confining pressure to cancel the radial strain. If the loading is started at zero confining pressure, the servo-control of the oil pressure is inaccurate, limiting our ability to drive properly the experiment in the range 0–1 MPa confining pressure. Below a differential stress of about 15 MPa and an axial strain of about 0.25 per cent (point A), the behaviour of the rock is essentially linear (see the dashed grey line). Beyond this threshold point A, and after a non-linear transition (between points A and B), the behaviour assumes another linear trend but with a different slope (between points B and C). Beyond point C the behaviour of the rock behaves essentially non-linearly and the differential stress reaches a plateau while the specimen is deforming. This overall behaviour is typical of a uniaxial strain experiment. At low stress, before point A the rock first exhibits a linear elastic behaviour. Between points A and B, at higher stress, the non-linear response of the rock suggests that energy-dissipative (irreversible) deformation might be at play. Following this non-linear stage, the specimen exhibits a linear response between points B and C, although with a different slope. Beyond point C, the specimen behaves essentially non-linearly, which is again attributed to irreversible deformation, and the non-zero slope of the curve in this region suggests the occurrence of strain hardening. Initiation of non-linear deformation might not always correspond to inception of irreversible deformation (damage/plasticity) as noted by Wong et al. (1997). However, non-linear reversible (elastic) deformation under increasing shear stress is not likely to occur in cemented granular materials at room temperature. The weakly cemented synthetic sandstone tested here would be more prone to grains re-arrangement and possible cement breakage. In the following, the stress strain curves obtained during this experiment will be numerically simulated using a reconstructed mesh from the microtomography data reported in the previous section. 3 NUMERICAL MODELLING 3.1 RVE Rocks are heterogeneous at all scales. There exists, however, a range of scales where the continuum assumption holds and material heterogeneity at smaller scale can be averaged into a RVE within a prescribed scale range. Therefore, for the estimation of macroscopic properties derived from microscale structural data, different positions and volume sizes within the imaged region are selected to conduct the numerical simulations and test the continuum assumption. The main source heterogeneity for the rock studied here is the presence of the pore space within the mineral phase. Upscaling methods allow us to derive a homogenised expression of a physical property of the medium, provided that a statistical RVE can be defined for that property (Charalambakis 2010; Shulakova et al.2016). The RVE is defined as the smallest subvolume of a heterogeneous medium that is of sufficient size for providing all structural information necessary for describing the homogenized property of interest (Nemat-Nasser & Hori 1993). It represents the effective (up-scaled) response of the microstructure and minimises the computational effort. Terada et al. (2000) proposed that asymptotic computational homogenisation techniques, using upper and lower bounds of finite element computations for converging solutions (with increasing volume size), can provide a reliable material parameter for the structure under investigation. The proposed method was demonstrated on 2-D models with randomly created microstructures for elastic properties. Regenauer-Lieb et al. (2010, 2013) extended the postulate to more complex deformation mechanisms and coupled problems based on the thermodynamic principle of maximum and minimum entropy production (Fig. 3). The upper bound corresponds to thermodynamic flux boundary condition (i.e. a constant velocity/displacement boundary condition for mechanical properties) and the lower bound corresponds to a thermodynamic force (i.e. a constant stress boundary condition for mechanical properties). The size of a RVE can be the volume size when a convergent result is obtained, for the physical process at hand (see also Veveakis & Regenaeur-Lieb 2015). Figure 3. View largeDownload slide Thermodynamic homogenisation procedure considering combinations of constant thermodynamic force and constant thermodynamic flux (from Regenauer-Lieb et al. 2013). Figure 3. View largeDownload slide Thermodynamic homogenisation procedure considering combinations of constant thermodynamic force and constant thermodynamic flux (from Regenauer-Lieb et al. 2013). Accordingly, a series of models were cropped step by step from the volume of microtomographic image as shown in Fig. 1. Each model is selected using two criteria: (1) the porosity should be close to that of the whole model, that is in the range of 24.0 ± 0.5 per cent in our case and (2) a smaller model is cropped from the nearby larger model. For each model, we analyse the response in terms of displacement and pressure/stress loads, respectively. The following boundary conditions are prescribed (coordinate frame refers to Fig. 1): (1) a zero normal displacement constraint is prescribed on the surfaces x = xmin, y = ymin and z = zmin; (2) surfaces x = xmax and y = ymax are free and (3) a displacement or stress loading is applied on the surface of z = zmax. Mathematically, this constraint is equivalent to the unconfined compression of a sample with mirror images of the volume in the x-, y- and z-directions. In the analysis of the mechanical response of the porous medium, only the solid skeleton is considered and the pores are treated as filled with a very compressible fluid (air), which is consistent with the experimental testing conditions reported in the previous section. Within the solid phase, quartz grains and calcite cement may play different roles. However, for sake of simplification in this study, we consider them as a solid mixture of unknown properties to focus on the influence of the pore-structure. Only the elastic response was analysed for these models in the determination of the RVE. Thus, only the elastic parameters of the solid skeleton are required as input. The magnitude of the loading is small enough to ensure only elastic deformation occurs. The RVE determined by convergence of the asymptotic homogenization is an elastic RVE. Based on the results reported by Regenauer-Lieb et al. (2010, 2013), we assume that the size of the RVE may also be a good first guess for plastic yielding. We test the validity of this assumption a posteriori in the next section. We emphasize that in this contribution only the onset of plastic yielding is considered. The post-yield hardening law will be the subject of future studies. For the evaluation of asymptotic homogenization trends eighteen cubes are cropped with a side-length L = 400, 320, 280, 240, 200, 180, 160, 140, 120, 100, 90, 80, 70, 60, 50, 40, 30 and 20. Finite element meshes are created for these cubes and two types of element are used—hexahedra or tetrahedra. For small volumes (100-cube and smaller), hexahedral elements are used, where one voxel of solid in the image data is converted to one hexahedral element in the finite element mesh. For all other volumes, tetrahedral elements are used in order to reduce computation time. The meshing is implemented by using a large number of triangles to represent the iso-surfaces separating the pore space from the solid. The tetrahedra are then generated from these triangles inward the solid phase. The Avizo® software package (https://www.fei.com/software/avizo3d/) is used to create these tetrahedral meshes. Elements created in this way have similar size and the size is controlled by the size of the iso-surface triangles specified in Avizo. The input elastic parameters of the solid phase (calcite-cemented quartz grains) are first taken to be: Young's modulus E = 100 GPa, Poisson's ratio $$\upsilon$$ = 0.2. At this stage, we do not claim that these values are representative of the actual elastic parameters of the minerals present in this sample. They are only used as reference values for the minerals in order to compare the effective (up-scaled) elastic properties of the porous medium for a series of models. The validity of these values of Young's modulus and Poisson's ratio of the solid phase will be discussed further in the next section. Abaqus® Explicit is used to perform the required finite element computations. The motivation underlying the use of an explicit solver is given in Liu et al. (2016) and is recalled in the discussion (Section 4). Fig. 4(a) shows the evolution with the model size of the computed elastic parameters using the upper and lower bounds. We see that the up-scaled Young's and shear moduli are lower for the stress loading boundary condition (lower bound) than for displacement boundary condition (upper bound). For models with L ≥ 240, stress and displacement loads give reasonably stable and convergent results. Overall, we notice that the stress boundary condition yields less stable predictions (larger oscillations) with increasing volume size than the displacement boundary condition. This instability (oscillations) seems to settle for L > 320. Figure 4. View largeDownload slide Asymptotic homogenization through upper and lower bounds for models of different sizes. (a) Elastic modulus and shear modulus; (b) Poisson's ratio. Figure 4. View largeDownload slide Asymptotic homogenization through upper and lower bounds for models of different sizes. (a) Elastic modulus and shear modulus; (b) Poisson's ratio. The average values of upper and lower bounds for L = 400 can be used as the asymptotic, convergent predictions, that is Young's modulus of 61.5 GPa and shear modulus of 28.2 GPa. The relative errors in the computed moduli for the various volumes and for both boundary conditions are listed in Table 1. The relative errors for all volumes lower than L = 400 are calculated by assuming that the average values for L = 400 (E = 61.5 GPa and μ = 28.2 GPa) are the true predictions. In view of the relative errors obtained, the RVE size can reasonably be taken to be L = 240. Table 1. Computed Young's and shear moduli for models with L ≥ 200, and their relative error compared to the asymptotic and convergent values obtained at L = 400. Size L  Upper bound E  Lower bound E  Upper bound μ  Lower bound μ    Value  Error (per cent)  Value  Error (per cent)  Value  Error (per cent)  Value  Error (per cent)  200  59.22  –3.70  52.48  −14.66  26.93  –4.58  23.42  –17.01  240  62.89  2.27  55.63  –9.53  28.96  2.59  25.26  –10.52  280  64.51  4.90  59.32  –3.53  29.60  4.88  26.77  –5.17  320  64.49  4.87  52.76  −14.20  29.86  5.79  24.18  –14.34  400  62.93  2.34  60.05  –2.34  28.93  2.49  27.52  –2.49  Size L  Upper bound E  Lower bound E  Upper bound μ  Lower bound μ    Value  Error (per cent)  Value  Error (per cent)  Value  Error (per cent)  Value  Error (per cent)  200  59.22  –3.70  52.48  −14.66  26.93  –4.58  23.42  –17.01  240  62.89  2.27  55.63  –9.53  28.96  2.59  25.26  –10.52  280  64.51  4.90  59.32  –3.53  29.60  4.88  26.77  –5.17  320  64.49  4.87  52.76  −14.20  29.86  5.79  24.18  –14.34  400  62.93  2.34  60.05  –2.34  28.93  2.49  27.52  –2.49  View Large 3.2 Drucker-Prager plasticity 3.2.1 Computational scheme We now proceed with the verification of this REV at yield. To model rock yield, we select the pressure-dependent Drucker-Prager criterion in which cohesion and friction angle can be used as intrinsic rock characteristics. A linearized version of the Drucker-Prager criterion is used  $$\ {\sigma _y} = {\sigma _n}{\rm{\ tan}}\beta + c,$$ (1)where σy is the yield stress, σn is the mean stress, c is the cohesion and β is the friction angle. Using eq. (1) c and β can be determined with at least two couples of σy and σn. Two cases are computed for the mechanical RVE in order to determine the cohesion and the friction angle of the rock. For both cases, the following boundary conditions are prescribed (coordinate frame refers to Fig. 1): (1) a zero normal displacement constraint on surfaces x = xmin, y = ymin and z = zmin; (2) the same displacement loading on z = zmax. For Case 1, the lateral surfaces x = xmax and y = ymax are free (similar to the elastic case described in Section 3.1); and for Case 2, a zero normal displacement is prescribed on these lateral surfaces x = xmax and y = ymax. The latter case corresponds to the uniaxial strain loading condition, which also corresponds to the experimental testing conditions reported in Section 2.2. For each case, the stress-strain relationship is computed and the yield stress σy can be obtained. The loading on z = zmax is numerically simulated up to a total strain of 3.5 per cent in the z (axial) direction. 3.2.2 Input parameters: Properties of the solid phase Similar to the computations of elasticity in the determination of the RVE, we need to specify the input parameters of the solid phase (minerals) in order to compute the overall plasticity parameters of the mechanical RVE. In addition to the elastic modulus and Poisson's ratio, two more parameters, yield stress and friction angle, are required as input for the Drucker-Prager plasticity in Abaqus computations. The solid phase (without pores) in our model is a mixture of quartz and calcite for which no common values are found in the literature. Often, values previously reported in the literature relate porous samples rather than to the pure solid mixture (e.g. Airey 1993). Thus, we investigate the parameters by fitting experimental data. To this end, we use the numerical model of Case 2 along with the corresponding experimental data reported in Fig. 2, that is uniaxial strain conditions. For the computed stress-strain relationships, stress is taken to be the differential stress averaged over all elements in the model; the strain is taken to be the relative shortening in the z-direction. The fitting procedure is done by trial-and-error and gradually narrowing the value ranges of the parameters. The Young's modulus is the first fitted as it is independent from the other two parameters, yield stress and friction angle are iteratively narrowed down. The detailed fitting results are given in the Appendix. We conclude that the parameters of the solid phase of the synthetic sandstone, which is composed of calcite-cemented quartz grains, are: Young's modulus between 9.1 and 9.3 GPa; yield stress between 35 and 38 MPa; and friction angle between 52° and 55°. Young's modulus appears to be much lower than values generally accepted for crystalline quartz, that is 70–100 GPa, or crystalline calcite, that is ∼86 GPa (Chen et al.2001). The value of Young's modulus derived for our synthetic sandstone is closer to that of the grey calcite marble (e.g. Martineze-Martineze et al.2012), or that of limestones (e.g. Yasar & Erdogen 2004). This suggests that the stiffness of the calcite-cemented sandstone is mostly controlled by the cement. The friction angle for the pure solid phase (no porosity) appears to be higher than the values often reported for porous sandstones or confined sand packs (Handin 1969; Jaeger & Cook 1976; Twiss & Moores 1992). It is however close to values reported for quartzite with a low to very low porosity (e.g. 48° for the friction angle of Sioux quartzite after Goodman 1989). As no unique and simple value for a given parameter is derived from fitting procedures (see the Appendix), for convenience, when we discuss these results we refer to a representative value for the corresponding range. In the following computations, a Young's modulus of 9100 MPa, a yield stress of 38 MPa and a friction angle of 52° are used as the properties of the solid phase (mineral assemblage composing the synthetic sandstone). 3.2.3 Output parameters: Properties of the porous rock The properties of the porous rock can be investigated with the fitted solid skeleton parameters. Fig. 5(a) represents the stress-strain relationship obtained for Case 1 and Case 2. The curves are plotted using the maximum compressive principal stress (average of all elements) and overall strain (relative shortening in the z-direction). For Case 1 the rock behaves in a typical elastic perfectly plastic way, whereas for Case 2 it exhibits a plastic hardening behaviour. Note that Case 1 represents a uniaxial stress condition with zero radial stress, which typically leads to a brittle failure of the actual rock beyond a peak stress point. Only a small amount of plastic deformation is expected before failure in this case. Because the model is not designed to simulate the brittle failure of rocks, only the maximum value of the stress reached during this simulated test is relevant for our analysis, that is ∼26 MPa. The plateau observed beyond a strain of ∼0.006 is not expected to replicate accurately the actual post-peak response of the actual rock in this brittle case. Compared to the yield stress of the solid skeleton (i.e. 38 MPa), the overall porous structure exhibits a yield stress about 30 per cent lower in Case 1. Figure 5. View largeDownload slide Stress-strain relationships (a) and mean stress—Mises stress relationships (b) of plastic response for the volume of L = 240. Figure 5. View largeDownload slide Stress-strain relationships (a) and mean stress—Mises stress relationships (b) of plastic response for the volume of L = 240. Fig. 5(b) shows the mean stress and the von Mises stress averaged over all elements in the model. Note that the von Mises stress and the mean stress after yielding correspond to σy and σn in eq. (1), respectively. For Case 1, the relationship between the von Mises stress and the mean stress is linear before yielding; both the von Mises stress and the mean stress appear constant after yielding. For Case 2, the von Mises stress and the mean stress show two different linear segments: one at the beginning of the deformation and another after yielding. The trend line of the points after yielding is connected to the clustered points obtained after yielding for Case 1. The slope and the intercept of this trend line are 0.4533 and 20.8, respectively. Fitting eq. (1) to this trend line, where Sy represents σy and Sn represents σn (see Fig. 5b) gives a cohesion of 20.8 MPa and a friction angle of 24.4°. These values are typical for porous sandstones (e.g. Goodman 1989) 3.3 Effect of the volume size An assumption used in this paper is that the size of the mechanical RVE determined from the asymptotic computational homogenization of elastic parameters is also a suitable first guess RVE for plastic yield. A more rigorous way to define the size of RVE would be to perform upper bound and lower bound plastic computations of all the series of models. Owing to the enormous computational cost this is not a very practical suggestion. We therefore developed the above alternative methodology for simplifying the workflow for future work. In order to prove the validity of the approach we performed plastic computations for very large volumes (L ≥ 240) and compare the so-derived plastic parameters. We postulate that our alternative method is robust if we obtain similar results. This would imply that mechanical RVEs can be derived by running simpler elastic analyses. Figs 6(a) and (b) plot the stress-strain relationships of models L = 240, 280, 320 and 400 for Case 1 and Case 2, respectively. Case 1 (Fig. 6a) shows an almost elastic perfectly plastic response for L < 320, and an emerging peak behaviour for L = 320 and 400, with a non-zero residual strength. Experimentally, in unconfined tests (Case 1) at room temperature, porous sandstones usually fail in a brittle way beyond the peak stress, while the stress drops to nearly zero. In this process, strain localisation occurs in the form of shear fractures. Damage prior to failure could be the dominant dissipative mechanism. Therefore, the present model is not expected to be adequate for this case as soon as strain localisation occurs: the RVE determined from elastic considerations becomes essentially irrelevant, plastic deformation localises within the localisation feature leaving the rest of the sample deforming essentially elastically. Figure 6. View largeDownload slide Stress-strain relationships of models L = 240, 280, 320 and 400 for Case 1 (a) and Case 2 (b). Figure 6. View largeDownload slide Stress-strain relationships of models L = 240, 280, 320 and 400 for Case 1 (a) and Case 2 (b). Case 2 (Fig. 6b) exhibits a plastic hardening behaviour for L = 240–400, and all curves are very similar. This confirms a posteriori the validity of the volume size L = 240 used in this case. Similar to Fig. 5(b), we plot the post-yield linear trends for Case 1 and Case 2 for the volumes L = 240, 280, 320 and 400 and fit the linearized Drucker-Prager criterion. We obtain a very narrow range of intercepts ranging from 20.8 to 21.8 MPa. The range of friction angles obtained is also quite narrow, that is 24.4°, 25.3°, 25.5° and 23.6° for L = 240, 280, 320 and 400, respectively. Based on the above results, we conclude that the mechanical RVE determined for elastic properties is also suitable for determining plastic yield parameters, provided that uniaxial strain experimental data (Case 2), as well as high-resolution CT-scan images of the sample, are available to calibrate the model. With this result, we anticipate that the experimental requirements, computational effort and hence the numerical costs can be reduced for similar studies of material characterisation in the future. 3.4 Effects of the porosity In order to investigate the influence of porosity on the yielding behaviour of the rock, we need several models with different porosity. The microstructures of these derivative models should be self-similar to ensure that the structure effect can be ruled out. To this end, we digitally create derivative models from the imaged sandstone sample by self-similarly inflating/deflating the pore structure. The shrinking procedure consists of moving the boundaries separating a pore and the surrounding solid towards the pore by one voxel in the x-, y- and z-directions. In contrast, the expansion procedure consists in moving this boundary towards the solid. This technique was initially developed to detect the percolation threshold of the imaged microstructure of a rock (Liu & Regenauer-Lieb 2011). As these derivative models have very similar microstructures but different porosity, they are ideally suited for investigating the influence of porosity on mechanical properties. The volume L = 240 (the elastic-plastic RVE derived earlier) is used in the following. In addition to the original model, four shrunk derivative models and one expanded model are derived and analysed (see Table 2). The porosity of these six models ranges from 8 to 28 per cent, which covers the most common range of porosity for typical reservoir rocks. The structure of all six models can be seen in Fig. 7. Figure 7. View largeDownload slide Structure of the six models for investigating the effect of porosity on Drucker-Prager plasticity. Pore space is reduced starting from (a) Md_P + 1, to (b) Md_P0, to (c) Md_P-1, to (d) Md_P-2, to (e) Md_P-3, to (f) Md_P-4. Figure 7. View largeDownload slide Structure of the six models for investigating the effect of porosity on Drucker-Prager plasticity. Pore space is reduced starting from (a) Md_P + 1, to (b) Md_P0, to (c) Md_P-1, to (d) Md_P-2, to (e) Md_P-3, to (f) Md_P-4. Table 2. Models used for investigating the effect of porosity on plasticity. Name  Porosity (per cent)  Description  Md_P + 1  28.14  One step of expanding of pores from original model  Md_P0  23.57  Original model  Md_P-1  19.14  One step of shrinking of pores from original model  Md_P-2  15.06  Two steps of shrinking of pores from original model  Md_P-3  11.39  Three steps of shrinking of pores from original model  Md_P-4  8.23  Four steps of shrinking of pores from original model  Name  Porosity (per cent)  Description  Md_P + 1  28.14  One step of expanding of pores from original model  Md_P0  23.57  Original model  Md_P-1  19.14  One step of shrinking of pores from original model  Md_P-2  15.06  Two steps of shrinking of pores from original model  Md_P-3  11.39  Three steps of shrinking of pores from original model  Md_P-4  8.23  Four steps of shrinking of pores from original model  View Large We have already presented the Drucker-Prager plastic computations of the original model of L = 240 in the previous section. Similar computations are conducted for Case 1 and Case 2 for the five derivative models. The stress-strain relationship for the six models are similar to the original model shown in Fig. 5(a)—an apparent elastic perfectly plastic behaviour is observed for Case 1, and a plastic hardening behaviour is observed for Case 2. The only difference for the models is: for a larger porosity, the plateau of the elastic perfectly plastic response is lower, and the plastic hardening is decreasing. Fig. 8 illustrates a few results obtained by fitting the linearized Drucker-Prager yield criterion to the linear trends obtained for Case 1 and 2. Fig. 8(a) corresponds to the model Md_P-4 with the lowest porosity (8.23 per cent). The value of von Mises stress post-yield is around 30 MPa for Case 1 and reaches over 150 MPa for Case 2. The slope of the fitting line is 0.9717 that corresponds to a friction angle of 44.2°. The cohesion for this model is 22.8 MPa corresponding to the intercept. Fig 8(b) shows the model Md_P-2 with a porosity of 15.06 per cent, the fitted friction angle is 35.6° and cohesion is 22.4 MPa; and Fig. 8(c) shows the model Md_P + 1 with the largest porosity of 28.1 per cent, the fitted friction angle is 19.6° and cohesion is 19.4 MPa. Fig. 8. View largeDownload slide Mean stress versus von Mises stress showing the fit of the linear Drucker-Prager relationship for describing the effects of porosity of models, (a) Md_P-4 (ϕ = 8.23 per cent); (b) Md_P-2 (ϕ = 15.06 per cent); (c) Md_P + 1 (ϕ = 28.14 per cent). Fig. 8. View largeDownload slide Mean stress versus von Mises stress showing the fit of the linear Drucker-Prager relationship for describing the effects of porosity of models, (a) Md_P-4 (ϕ = 8.23 per cent); (b) Md_P-2 (ϕ = 15.06 per cent); (c) Md_P + 1 (ϕ = 28.14 per cent). Using the data of all six models, we plot the relationships between the plasticity parameters and the porosity in Fig. 9. The yield stress, the cohesion and the friction angle decrease in an essentially linear way with increasing porosity. The yield stress and the friction angle for a porosity of 0 are 39 MPa and 54.4°, matching the input parameters of the solid phase. We notice that the friction angle decreases more significantly than cohesion with increasing porosity. Figure 9. View largeDownload slide Plastic parameters versus porosity. Figure 9. View largeDownload slide Plastic parameters versus porosity. 4 DISCUSSION We have documented in this paper a new approach for assessing the effects of microstructure on key mechanical properties of rocks. The approach differs from conventional micromechanical analyses of plasticity in that it is not restricted to randomly created mathematical composite models with specific geometries and volume fractions. We have instead studied the plastic response of a synthetic calcite-cemented porous sandstone (CIPS) based on the 3-D images of its microstructures by X-ray microtomography. Although the CIPS is a relatively simple synthetic rock simulating an idealized rock, its inherent complex topology and connectivity have significant impact on the macroscopic properties. It is a step further toward complexity compared to a digitally generated rock microstructure. The computational and meshing aspects and additional complexities of using X-ray microtomographic images of rocks are discussed in detail in Liu et al. (2016). Computational efficiency and convergence issues were solved by selecting a parallel explicit finite element technique (Abaqus® Explicit) which provided numerically acceptable results. The use of an implicit solver such as Abaqus® Standard on such complex microstructures was unsuccessful as the computations did not converge. Numerical issues were encountered through mesh sensitivity affecting the final results of the computations. For instance, the difference in the resulting yield stress for different discretization (hexahedral versus tetrahedral meshes at different resolutions) may be 10 per cent. This error is much higher than for the estimation elastic modulus as reported by Liu et al.2016 (Table 2). Meshing error mainly occurs at the initial step of the meshing process because stepped voxel surfaces must be replaced by triangular iso-surfaces. Doing so, critical connections between grains may break or be reinforced, and these connections significantly affect the mechanical response due to possible stress/strain concentration in these areas (see Fig. 10). With such additional complexity, it is difficult to estimate the extent of the error caused by the meshing. In principle, the finer mesh, the more accurate results. In this study, we did not fully investigate this specific source of error. However, we checked the model with the coarsest mesh, then re-meshed more finely and recomputed all the models listed in Figs 5–9 (and in the Appendix) using the exact same finer mesh, that is the maximum size of the triangles was the same as in Avizo® to generate the required tetrahedral elements. Using the same mesh resolution is crucial to obtain overlapping curves in Fig. 6 and clearly linear relationships in Fig. 9. Furthermore, the intercepts at zero porosity of the fitted yield stress and of the angle of friction lines match the input skeleton parameter reasonably well. This suggests that the discrete error in this study did not systematically overestimate or underestimate the trend of plastic parameters. Figure 10. View largeDownload slide Equivalent plastic strain on two typical slices of the model L = 240 in Case 2. Plastic strain concentrates at grain boundaries. Figure 10. View largeDownload slide Equivalent plastic strain on two typical slices of the model L = 240 in Case 2. Plastic strain concentrates at grain boundaries. For the purpose of the current pilot study models with <10 per cent error were assumed acceptable. Future work should focus more on minimizing numerical artifacts through comparison of different element types (e.g. homogeneous or adaptive cubes, homogeneous or adaptive tetrahedra). Since the calculations are performed on several millions 3-D finite elements a comparison of the resulting accuracy in light of computational costs is necessary. The RVE size was determined by asymptotic computational homogenization deriving the upper and lower bounds of finite element analyses. Owing to the high computation cost this technique was only practical for deriving the size of the elastic RVE. This size of the mechanical RVE is strictly only valid for elastic predictions. We proposed that the elastic RVE is also a good first guess of the RVE for the plasticity calculations. The analysis of the effect of volume size in Section 4.1 showed a posteriori that volume sizes from 240-cube to 400-cube (voxels) have almost exactly the same plastic response. This implies that the RVE size determined by the upper and lower bounds elastic computations is acceptable for plastic yielding. The error in this assumption is less than error introduced by the numerical artefacts of meshing and explicit solution technique (<10 per cent). For obtaining more accurate results future work should consider a full asymptotic derivation of the plastic RVE. This could become possible when the above described numerical resolution and convergence issues are fully resolved in an implicit solver. Since plastic deformation can propagate to larger scales and post-yield bifurcations can introduce new scaling lengths through introduction of the width of shear bands we only discuss initial plastic yield in this paper. We emphasize that the RVE size derived from elastic analysis is inappropriate for the later stages of plastic deformation after yielding. This statement applies specifically if strain localization effects are resolved and shear bands or compaction bands are propagating through the sample. While investigating the bulk plastic response of a porous rock, pore space is treated as mechanically inert in this study. The choice of the input parameters of the solid phase is crucial for the accuracy of the resulting bulk properties. For the specific sandstone sample studied here, the skeleton is made of calcite-cemented quartz grains. There is not experimental data reporting the mechanical parameters of this specific combination of minerals. By fitting the experimental curve of the bulk response of the rock to uniaxial strain loading we were able to estimate the parameters required for the solid phase: Young's modulus was estimated between 9100 and 9300 MPa and the yield stress between 35 and 38 MPa. The friction angle was bracketed between 52° and 55°. As expected the value of Young's modulus obtained is lower than the commonly accepted value for pure quartz. This implies that composite quartz-calcite skeleton appears to have a much lower elastic modulus than the pure quartz or calcite constituents. The friction angle of the solid is higher than that usually reported for sandstones, but close to that reported for (low porosity) quartzites. The output bulk parameters obtained in this study are all within reasonable ranges. Quartz and calcite have different mechanical parameters (Ahrens 1995). We modelled this synthetic sandstone as two-phase medium comprising the pore phase and the solid phase. The latter is considered as a mixture of two solids (quartz and calcite), the properties of which are a priori unknown. We therefore focused on the impact of the pore space, as derived from the 3-D X-ray microtomographic images, on the macroscopic mechanical behaviour of the rock, especially the post-elastic domain of deformation. Previous studies of the plastic behaviour of rocks in relation to their microstructure have been reported by (Liu et al.2012, 2015). In this study, we focus specifically on the impact of the size of the RVE and of the porosity on the macroscopic yield parameters of the rock assuming that the solid phase is a mixture of quartz and calcite. This is meant to be a first step for benchmarking future studies in which the arrangement and respective properties of the constituting minerals (quartz and calcite) could then be added. However, this study already suggests (see Fig. 10) that the macroscopic deformation of the porous rock is strongly controlled by the cemented grain contacts where high strain values occur. Interesting results of this paper are the simple linear relationships providing quantitative insights into the effect of microstructure on macroscopic rock properties. A good example is the reduction of the friction angle to with increasing porosity as shown in Fig. 9. The results show that the yield stress, the cohesion and the friction angle of the rock linearly decrease with increasing porosity, for a range of porosity comprised between 8 and 28 per cent. This range covers the porosity values encountered for most conventional reservoir rocks. Because the macroscopic behaviour of the rock is found to be strongly controlled by its microstructure, and since the models studied here are geometrically self-similar, we do not recommend to extrapolate the relationships derived here to models with very different structures. More work is required in that case. It is not recommended either to extrapolate the linear relationships obtained here further to an untested range. According to percolation theory, parameters such as yield stress change exponentially when the solid volume fraction approaches the percolation threshold (Stauffer & Aharony 1994; Sahimi 1998). Based on the linear trends obtained in Fig. 9, it appears that the volume fractions of the solid in the models analysed here are far from the percolation threshold. According to the basic rule of mixture (Voigt 1889; Reuss 1929; Alger 1997), the property Ec of a composite material made of two phases should be comprised between an upper bound Ec = fEf + (1 − f)Em and a lower bound $${E_c} = {( {\frac{f}{{{E_f}}} + \frac{{1 - f}}{{{E_m}}}} )^{ - 1}}\$$, where f is the volume fraction of the inclusions, and Em and Ef are the properties of the matrix and the inclusions, respectively. In our two-phase model (solid and pores), f corresponds to porosity, and pores are voids with Ef = 0. We obtained a linear relation between the parameters of the solid phase and the macroscopic properties. Note that the expression of the upper bound leads to a similar linear relation with a slope of (1 − f), which differs from the slope we found in Fig. 9. Because the upper and lower bounds have been derived for the linear elastic deformation regime, this difference is not surprising. 5 CONCLUSIONS In this pilot study, we have clearly shown a promising new path for evaluating rock mechanical properties from micro-CT imaging. The analysis has been deliberately chosen to be close to the ideal mathematical rock models by the choice of imaging a synthetic rock composite. We have shown that the suggested computational homogenization workflow is suitably robust and reliable within a 10 per cent error of the mechanical properties. While the described computational homogenization method is well established (Terada et al.2000), we have introduced three new elements: (1) generalization of the classical computational homogenization methods in terms of thermodynamic bounds; (2) proposed equivalence of the elastic RVE to a plastic RVE at initial yield and (3) calculation of quantitative microstructural trends derived from derivative models (inflated/deflated microstructures). The successful test of these new improvements on a synthetic porous sandstone opens promising avenues for future assessments of the effect of microstructure on mechanical properties for complex microstructures encountered in rock mechanics. Future development includes the consideration of multiphase involved in a sample, such as from quartz, calcite and pores to quartz, calcite and saturated/unsaturated fluid and the consideration of other constitutions, such as viscoplasticity and elasto-viscoplasticity at different temperature and pressure. The approach developed here and its upgraded version will ultimately lead to new ways of characterizing rock properties while relying less on time consuming and costly laboratory measurements, with potential application in crustal strength and deformation in geophysics, borehole stability, hydraulic fracturing, or reservoir compaction studies in the oil and gas, geothermal or CO2 geo-sequestration industries. Acknowledgements We thank Keyu Liu, China University of Petroleum (East China), formerly CSIRO Energy, provided the synthetic sandstone sample for experiments and the micro CT images. Jie Liu would like to acknowledge the support of National Natural Science Foundation of China grants no. 41574087. Part of the computing was performed on Tianhe-2 supercomputer and the computing time was supported by the Special Program for Supercomputing of the NSFC-Guangdong Joint Fund. Klaus Regenauer-Lieb would also like to acknowledge support of the Australian Research Council through ARC Discovery grants no. DP140103015, DP170104550, DP170104557. REFERENCES Ahrens T.J., 1995. Mineral Physics and Crystallography: A Handbook of Physical Constants , American Geophysical Union. Google Scholar CrossRef Search ADS   Airey D.W., 1993. Triaxial testing of naturally cemented carbonate soil. J. Geotech. Eng. , 119( 9), 1379– 1398. Google Scholar CrossRef Search ADS   Alger M.S.M., 1997. Polymer Science Dictionary.  2nd edn. Springer Publishing. ISBN 0412608707. Almqvist B.S.G., Mainprice D., Madonna C., Burlini L., Hirt A.M., 2011. Application of differential effective medium, magnetic pore fabric analysis, and X-ray microtomography to calculate elastic properties of porous and anisotropic rock aggregates, J. geophys. Res.: Solid Earth , 116, 1– 17. Google Scholar CrossRef Search ADS   Arns C.H., Knackstedt, M.A., Pinczewski W.V., Lindquist W.B., 2001. Accurate estimation of transport properties from microtomographic images, Geophys. Res. Lett. , 28, 3361– 3364. Google Scholar CrossRef Search ADS   Arns C.H., Knackstedt M.A., Pinczewski W.V., Garboczi E.J., 2002. Computation of linear elastic properties from microtomographic images: methodolgy and agreement between theory and experiment, Geophysics , 67( 5), 1396– 1405. Google Scholar CrossRef Search ADS   Arns C.H., Knackstedt M.A., Martys N.S., 2005. Cross-property correlation and permeability estimation in sandstone, Phys. Rev. E , 72, 046304, doi:10.1103/PhysRevE.72.046304. Google Scholar CrossRef Search ADS   Blunt M.J.,et al.  , 2013. Pore-scale imaging and modelling, Adv. Water Resour. , 51, 197– 216. Google Scholar CrossRef Search ADS   Charalambakis N., 2010. Homogenization techniques and micromechanics. A survey and perspectives, Appl. Mech. Rev. , 63( 3), 030803. Google Scholar CrossRef Search ADS   Chen C.C., Lin C.C., Liu L.G., Sinogeijin S.V, Bass J.D., 2001. Elasticity of single-crystal calcite and rhodochrosite by Brillouin spectroscopy, Am. Mineral. , 86, 1525– 1529. Google Scholar CrossRef Search ADS   Chen C., Packman A.I, Zhang D., Gaillard J.-F., 2010. A multi-scale investigation of interfacial transport, pore fluid flow, and fine particle deposition in a sediment bed, Water Resour. Res. , 46, W11560, doi:10.1029/2009WR009018. Fredrich J.T., DiGiovanni A.A., Noble D.R. 2006. Predicting macroscopic transport properties using microscopic image data, J. geophys. Res. , 111, B03201, doi:10.1029/2005JB003774. Google Scholar CrossRef Search ADS   Goodman R.E., 1989. Introduction to Rock Mechanics , 2. Wiley. Handin J., 1969. On the Coulomb–Mohr failure criterion, J. geophys. Res. , 74, 5343– 5348. Google Scholar CrossRef Search ADS   Jaeger J.C., Cook N.G.W., 1976. Fundamentals of Rock Mechanics , Chapman & Hall, Wiley, 585 pp. Knackstedt M.A.et al.  , 2006. Elastic and transport properties of cellular solids derived from three-dimensional tomographic images, Proc. R. Soc. A , 462, 2833– 2862. CrossRef Search ADS   Liu J., Regenauer-Lieb K., 2011. Application of percolation theory to microtomography of structured media: percolation threshold, critical exponents, and upscaling, Phys. Rev. E , 83, 1– 13. Liu J., Regenauer-Lieb K., Hines C., Liu K., Gaede O., Squelch A., 2009. Improved estimates of percolation and anisotropic permeability from 3-D X-ray microtomography using stochastic analyses and visualization, Geochem., Geophys., Geosyst. , 10, doi:10.1029/2008GC002358. Liu J., Freij-Ayoub R., Regenauer-Lieb K., 2012. Determination of the plastic strength of carbonates from microtomography and upscaling using percolation theory, in Proceedings of the ECCOMAS 2012 , ed. Eberhardsteiner J.et al.  , Vienna, Austria, September 10–14, 2012. Liu J., Freij-Ayoub R., Regenauer-Lieb K., 2015. Rock plasticity from microtomography and upscaling, J. Earth Sci. , 26( 1), 053– 059. Google Scholar CrossRef Search ADS   Liu J., Pereira G.G., Liu Q., Regenauer-Lieb K., 2016. Computational challenges in the analyses of petrophysics using microtomography and upscaling: a review, Comput. Geosci. , 89, 107– 117. Google Scholar CrossRef Search ADS   Martinez-martinez J., Benavente D., Angeles G.-C.M., 2012. Comparison of the static and dynamic elastic modulus in carbonate rocks, Bull. Eng. Geol. Environ. , 71, doi:10.1007/s10064-011-0399-y. Nemat-Nasser S., Hori M., 1993. Micromechanics: Overall Properties of Heterogeneous Materials, Elsevier. Regenauer-Lieb K.et al.  , 2010. Time-dependent, irreversible entropy production and geodynamics, Phil. Trans. R. Soc. Lond., A , 368( 1910), 285– 300. Google Scholar CrossRef Search ADS   Regenauer-Lieb K.et al.  , 2013. Multiscale coupling and multiphysics approaches in earth science: theory, J. Coupled Syst. Multiscale Dyn. , 1( 1), 49– 73. Google Scholar CrossRef Search ADS   Reuss A., 1929. Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle, ZAMM - J. Appl. Math. Mech./ Zeitschrift für Angewandte Mathematik und Mechanik , 9, 49– 58. Google Scholar CrossRef Search ADS   Roberts A.P., Garboczi E.J., 2002. Elastic properties of model random three-dimensional open-cell solids, J. Mech. Phys. Solids , 50( 1), 33– 55. Google Scholar CrossRef Search ADS   Sahimi M., 1998. Non-linear and non-local transport processes in heterogeneous media: from long-range correlated percolation to fracture and materials breakdown, Phys. Rep. , 306, 213– 395. Google Scholar CrossRef Search ADS   Shulakova V.et al.  , 2013. Computational elastic up-scaling of sandstone on the basis of X-ray micro-tomographic images, Geophys. Prospect. , 61( 2), 287– 301. Google Scholar CrossRef Search ADS   Shulakova V., Sarout J., Pimienta L., Lebedev M., Mayo S., Clennell M.B., Pervukhina M., 2016. Effect of supercritical CO2 on carbonates: savonnieres sample case study, Geophys. Prospect. , 61, 287– 301. Google Scholar CrossRef Search ADS   Stauffer D., Aharony A., 1994. Introduction to Percolation Theory , Taylor & Francis. Spanne P., Thovert J.F., Jacquin C.J., Lindquist W.B., Jones K.W., Adler P.M., 1994. Synchrotron computed microtomography of porous media: topology and transports, Phys. Rev. Lett. , 73, 2001– 2004. Google Scholar CrossRef Search ADS PubMed  Stock S.R., 2008. Recent advances in X-ray microtomography applied to materials, 53( 3), 129– 181. Terada K., Hori M., Kyoya T., Kikuchi N., 2000. Simulation of the multi-scale convergence in computational homogenization approaches, Int. J. Solid Struct. , 37, 2285– 2311. Google Scholar CrossRef Search ADS   Twiss R.J., Moores E.M., 1992. Structural Geology , W.H. Freeman and Company, 532 pp. Veveakis E., Regenauer-Lieb K., 2015. Review of extremum postulates, Curr. Opin. Chem. Eng. , 7, 40– 46. Google Scholar CrossRef Search ADS   Voigt W., 1889. Ueber die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper, Annalen der Physik , 274, 573– 587. Google Scholar CrossRef Search ADS   Wong T.-F., David C., Zhu W., 1997. The transition from brittle faulting to cataclastic flow in porous sandstones: mechanical deformation, J. geophys. Res. , 102( B2), 3009– 3025. Google Scholar CrossRef Search ADS   Yang S., Liu K., Mayo S., Tulloh A., 2012. CIPS sandstone microstructure, v2, CSIRO Data Collection , doi:org/10.4225/08/5476787A1A50F. Yasar E., Erdogan Y., 2004. Correlating sound velocity with the density, compressive strength and Young's modulus of carbonate rocks, Int. J. Rock Mech. Mining Sci. , 41, 871– 875. Google Scholar CrossRef Search ADS   APPENDIX Here we show the comparisons of experimental curve with computed curves that the ranges of elastic modulus, yield stress and the angle of friction have been preliminarily obtained. Fig. A1 Shows the stress-strain relationships derived for Case 2, at small strain, for three values of Young's modulus, and the comparison with the experimental curve. We see both 9100 and 9300 MPa are almost parallel to the experimental curve—the experimental curve is not starting from zero but at a small value of stress, the computed curve with E = 9100 and 9300 MPa can be seen as the best fit of the experimental curve when strain is less than 0.005. In these computations, the yield stress and friction angle are taken to be 43 MPa and 47°, respectively. These two yield parameters do not affect the deformation response in the elastic domain. Therefore, the elastic modulus of the solid skeleton in this sample is estimated to be between 9100 and 9300 MPa. Figure A1. View largeDownload slide Three computed stress-strain curves of different elastic moduli comparing with experimental measurement to investigate the elastic modulus of the solid. Figure A1. View largeDownload slide Three computed stress-strain curves of different elastic moduli comparing with experimental measurement to investigate the elastic modulus of the solid. Fig. A2 compares the stress-strain relationships using three different yield stress values as input solid parameter with the experimental curve. In these computations, Young's modulus and the friction angle are 9100 MPa and 47°, respectively. We see that when the yield stress of the solid phase is taken to be 35 MPa, the computed stress-strain curve is lower than the experimental curve. However, the trends of the computed and experimental curves are very similar (parallel curves). Overall, a better fit over the wider strain range of 0–0.01 is obtained with a yield stress value of 38–41 MPa. Figure A2. View largeDownload slide Three computed stress-strain curves of different yield stress comparing with experimental measurement to investigate the yield stress of the solid. Figure A2. View largeDownload slide Three computed stress-strain curves of different yield stress comparing with experimental measurement to investigate the yield stress of the solid. Fig. A3 compares the stress-strain relationships from different values of the friction angle. In these computations, we used a Young's modulus of 9100 MPa and a yield stress value of 38 MPa. From the figure, we see that 50° is too low; 52° is slightly lower than the experimental curve in the strain range 0.01–0.03, but replicates quite accurately the yield point; 55° matches the experimental curve very well except for the high values of strain. Therefore, the friction angle of the solid phase is estimated to be constrained between 52° and 55°. Figure A3. View largeDownload slide Three computed stress-strain curves of different friction angle comparing with experimental measurement to investigate the angle of friction of the solid. Figure A3. View largeDownload slide Three computed stress-strain curves of different friction angle comparing with experimental measurement to investigate the angle of friction of the solid. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Computational upscaling of Drucker-Prager plasticity from micro-CT images of synthetic porous rock

, Volume 212 (1) – Jan 1, 2018
13 pages

/lp/ou_press/computational-upscaling-of-drucker-prager-plasticity-from-micro-ct-ARJf84jhpk
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx409
Publisher site
See Article on Publisher Site

### Abstract

Abstract Quantifying rock physical properties is essential for the mining and petroleum industry. Microtomography provides a new way to quantify the relationship between the microstructure and the mechanical and transport properties of a rock. Studies reporting the use microtomographic images to derive permeability and elastic moduli of rocks are common; only rare studies were devoted to yield and failure parameters using this technique. In this study, we simulate the macroscale plastic properties of a synthetic sandstone sample made of calcite-cemented quartz grains using the microscale information obtained from microtomography. The computations rely on the concept of representative volume elements (RVEs). The mechanical RVE is determined using the upper and lower bounds of finite-element computations for elasticity. We present computational upscaling methods from microphysical processes to extract the plasticity parameters of the RVE and compare results to experimental data. The yield stress, cohesion and internal friction angle of the matrix (solid part) of the rock were obtained with reasonable accuracy. Computations of plasticity of a series of models of different volume-sizes showed almost overlapping stress-strain curves, suggesting that the mechanical RVE determined by elastic computations is also valid for plastic yielding. Furthermore, a series of models were created by self-similarly inflating/deflating the porous models, that is keeping a similar structure while achieving different porosity values. The analysis of these models showed that yield stress, cohesion and internal friction angle linearly decrease with increasing porosity in the porosity range between 8 and 28 per cent. The internal friction angle decreases the most significantly, while cohesion remains stable. Microstructure, Plasticity, diffusion and creep, Numerical modelling 1 INTRODUCTION Understanding upscaling of microscale physical processes to the scale of a core is a standard problem for the hydrocarbon, geothermal, water and mining industries. Elasticity, plasticity, permeability and other properties of rocks are related to the structure of pores and cracks as well as their connectivity. Microtomography enables the determination of the internal structure of rocks at the micro- to nanoscales and allows quantification of the relationship between microstructure and mechanical and transport properties through computational forward simulations of the underlying processes (Spanne et al.1994; Stock 2008; Blunt et al.2013). Numerical computations on microtomographic data have shown good agreement with experimental data for fluid flow and elastic properties (Arns et al.2001, 2005; Fredrich et al.2006; Knackstedt et al.2006; Chen et al.2010). However, the extension of the methodology to plastic properties still represents a poorly understood area. One aspect is the difficulty of defining a representative volume element (RVE) for dissipative processes such as plastic deformation from microtomographic data. This paper aims to address this problem by suggesting an energy-based up-scaling procedure (Veveakis & Regenauer-Lieb 2015). The study of rock properties based on microtomographic images has been developing for two decades since the middle of 1990s. Permeability is the first rock property that has been assessed through computational homogenisation method based on microtomographic structure (Spanne et al.1994). Arns et al. (2001) followed on from this early pioneering work and performed a benchmark analysis of transport properties derived from microtomographic images of Fontainebleau sandstone. Porosity and permeability, and their scale dependency have also been analysed for other porous rocks (Fredrich et al.2006; Chen et al.2010). The power of the new method led to many other works. Elastic properties derived from microtomographic images were, for instance, studied by Arns et al. (2002). Roberts & Garboczi (2002) computed the linear elastic properties of random porous materials with a wide variety of microstructures. Almqvist et al. (2011) applied differential effective medium and X-ray microtomography to calculate the elastic properties of porous and anisotropic rock aggregates; Shulakova et al. (2013) also computed the elastic properties of natural sandstones from microtomographic images; Knackstedt et al. (2006) analysed the correlation between the up-scaled diffusivity, elasticity, permeability and conductivity of industrial cellular foams. Rock deformation in nature occurs far beyond elasticity and no robust method for upscaling post-elastic rock behaviour exists to date. In this work, we discuss the derivation of plastic properties of rocks as a natural extension of the methods introduced for computational homogenisation of elastic properties. In doing so, we first define the yield envelope and the plastic potential of the solid skeleton of the porous rock. Laboratory work on porous silicate rocks show that Drucker-Prager and Mohr-Coulomb adequately describe plastic yield. Because of its smooth yield surface, the Drucker-Prager criterion appears better suited than Mohr-Coulomb criterion for numerical simulations. Liu et al. (2012, 2015) simulated the Drucker-Prager plastic response of carbonate samples at the microscale using microtomography. However, due to the extremely strong heterogeneity of carbonates involving at least two different scales of pores, those pilot studies remain inconclusive. The available micro-CT scans, were too small to convincingly test the validity of RVE and hence the results on the relationship between porosity and plastic properties could not be fully understood. In this study, we therefore analyse a well-defined synthetic sandstone sample and derive the yield stress, cohesion and friction angle from 3-D X-ray images of the microstructure. The computations are based on the RVE defined from elastic computations, and a subsequent test of RVE for computational homogenisation of plastic properties. In a second step, the computed properties are compared with experimental measurements. Finally, the dependence of plastic properties on porosity is deduced from a series of models with similar structure but different porosity values. 2 ROCK MATERIAL AND EXPERIMENTAL DATA 2.1 Microtomographic data In this study, we used a CIPS (Calcite In situ Precipitation System, a patent technique) synthetic sandstone sample composed of pure quartz grains (diameter 0.15–0.35 mm) and cemented by calcite. The data collected for this CIPS is freely available from CSIRO data repository (Yang et al.2012). The data includes 3-D X-ray microtomography images and mineral phase distributions. A cylindrical plug, 38 mm in diameter, was drilled from the man-made block. The gas permeability of the sample was measured and compared well with the numerically simulated permeability from the point of view of percolation theory (Stauffer & Aharony 1994; Liu & Regenauer-Lieb 2011). A small sample, 5 mm in diameter, was prepared and X-ray CT-scanned with a resolution of 2.5 μm. After segmentation, a porosity of 24.4 per cent was derived from the images, which is lower than the experimental measurement of Helium porosity of 28.5 per cent. This discrepancy is most likely due resolution limits preventing imaging smaller size pores below 2.5 μm. Detailed information about image processing and geometrical characteristics of this sample can be found in Liu et al. (2009). The processed 3-D image of a subvolume of the scanned sample is shown in Fig. 1. For statistical analysis we use a Cartesian reference frame and investigate percolation in three orthogonal directions thus show cropped parallelepipeds out of the cylindrical sample. Figure 1. View largeDownload slide Visualization of microtomographic data of the synthetic sandstone sample analysed in this paper. (a) 3-D microstructure, in the top slice the binary image shows the solid in black and the pores in white (volume size is 600 × 600 × 400 voxels); (b) the coordinate frame used. Figure 1. View largeDownload slide Visualization of microtomographic data of the synthetic sandstone sample analysed in this paper. (a) 3-D microstructure, in the top slice the binary image shows the solid in black and the pores in white (volume size is 600 × 600 × 400 voxels); (b) the coordinate frame used. 2.2 Rock deformation data A cylindrical CIPS plug, 25.4 mm in diameter and 47 mm in length, was subjected to an axial loading while maintaining uniaxial strain conditions. The experiment is conducted in a conventional (axisymmetric) triaxial stress vessel, using oil as the confining medium acting mechanically on the rock through a sealing flexible membrane. It consists in applying an increasing axial displacement on the specimen's end face at a constant rate (1.8 × 10−5 s−1 in this instance) while the oil pressure is servo-controlled to cancel the radial strain induced by Poisson's effect. Therefore, the rock specimen is nominally subjected to a uniaxial strain field, similar to a conventional oedometric tests on soils. Therefore, the strain field within the specimen is nominally uniaxial at mid-height of the specimen, that is away from the frictional boundary conditions induced by the loading platens. The specimen is tested in dry conditions (no pore fluid) so that there is no internal pore pressure within the rock, and therefore the external stresses applied to the specimen are equal to the effective stresses. Fig. 2 depicts the evolution of the differential, radial and mean effective stresses as a function of the applied axial strain during the experiment. For sake of simplicity the nominally null radial strain is not represented in this figure. Figure 2. View largeDownload slide Evolution of the differential stress (black curve), effective mean stress (dotted curve), and radial stress (dashed curve) with increasing axial strain during the uniaxial strain experiment. The radial strain remains nominally zero. Figure 2. View largeDownload slide Evolution of the differential stress (black curve), effective mean stress (dotted curve), and radial stress (dashed curve) with increasing axial strain during the uniaxial strain experiment. The radial strain remains nominally zero. Note the stresses in Fig. 2 are not starting from zero. This is because of a technical limitation of the equipment used to conduct this uniaxial strain test in a conventional triaxial stress vessel: we have to apply a small confining pressure (1 MPa) prior to applying the displacement-controlled axial stress and servo-controlling the confining pressure to cancel the radial strain. If the loading is started at zero confining pressure, the servo-control of the oil pressure is inaccurate, limiting our ability to drive properly the experiment in the range 0–1 MPa confining pressure. Below a differential stress of about 15 MPa and an axial strain of about 0.25 per cent (point A), the behaviour of the rock is essentially linear (see the dashed grey line). Beyond this threshold point A, and after a non-linear transition (between points A and B), the behaviour assumes another linear trend but with a different slope (between points B and C). Beyond point C the behaviour of the rock behaves essentially non-linearly and the differential stress reaches a plateau while the specimen is deforming. This overall behaviour is typical of a uniaxial strain experiment. At low stress, before point A the rock first exhibits a linear elastic behaviour. Between points A and B, at higher stress, the non-linear response of the rock suggests that energy-dissipative (irreversible) deformation might be at play. Following this non-linear stage, the specimen exhibits a linear response between points B and C, although with a different slope. Beyond point C, the specimen behaves essentially non-linearly, which is again attributed to irreversible deformation, and the non-zero slope of the curve in this region suggests the occurrence of strain hardening. Initiation of non-linear deformation might not always correspond to inception of irreversible deformation (damage/plasticity) as noted by Wong et al. (1997). However, non-linear reversible (elastic) deformation under increasing shear stress is not likely to occur in cemented granular materials at room temperature. The weakly cemented synthetic sandstone tested here would be more prone to grains re-arrangement and possible cement breakage. In the following, the stress strain curves obtained during this experiment will be numerically simulated using a reconstructed mesh from the microtomography data reported in the previous section. 3 NUMERICAL MODELLING 3.1 RVE Rocks are heterogeneous at all scales. There exists, however, a range of scales where the continuum assumption holds and material heterogeneity at smaller scale can be averaged into a RVE within a prescribed scale range. Therefore, for the estimation of macroscopic properties derived from microscale structural data, different positions and volume sizes within the imaged region are selected to conduct the numerical simulations and test the continuum assumption. The main source heterogeneity for the rock studied here is the presence of the pore space within the mineral phase. Upscaling methods allow us to derive a homogenised expression of a physical property of the medium, provided that a statistical RVE can be defined for that property (Charalambakis 2010; Shulakova et al.2016). The RVE is defined as the smallest subvolume of a heterogeneous medium that is of sufficient size for providing all structural information necessary for describing the homogenized property of interest (Nemat-Nasser & Hori 1993). It represents the effective (up-scaled) response of the microstructure and minimises the computational effort. Terada et al. (2000) proposed that asymptotic computational homogenisation techniques, using upper and lower bounds of finite element computations for converging solutions (with increasing volume size), can provide a reliable material parameter for the structure under investigation. The proposed method was demonstrated on 2-D models with randomly created microstructures for elastic properties. Regenauer-Lieb et al. (2010, 2013) extended the postulate to more complex deformation mechanisms and coupled problems based on the thermodynamic principle of maximum and minimum entropy production (Fig. 3). The upper bound corresponds to thermodynamic flux boundary condition (i.e. a constant velocity/displacement boundary condition for mechanical properties) and the lower bound corresponds to a thermodynamic force (i.e. a constant stress boundary condition for mechanical properties). The size of a RVE can be the volume size when a convergent result is obtained, for the physical process at hand (see also Veveakis & Regenaeur-Lieb 2015). Figure 3. View largeDownload slide Thermodynamic homogenisation procedure considering combinations of constant thermodynamic force and constant thermodynamic flux (from Regenauer-Lieb et al. 2013). Figure 3. View largeDownload slide Thermodynamic homogenisation procedure considering combinations of constant thermodynamic force and constant thermodynamic flux (from Regenauer-Lieb et al. 2013). Accordingly, a series of models were cropped step by step from the volume of microtomographic image as shown in Fig. 1. Each model is selected using two criteria: (1) the porosity should be close to that of the whole model, that is in the range of 24.0 ± 0.5 per cent in our case and (2) a smaller model is cropped from the nearby larger model. For each model, we analyse the response in terms of displacement and pressure/stress loads, respectively. The following boundary conditions are prescribed (coordinate frame refers to Fig. 1): (1) a zero normal displacement constraint is prescribed on the surfaces x = xmin, y = ymin and z = zmin; (2) surfaces x = xmax and y = ymax are free and (3) a displacement or stress loading is applied on the surface of z = zmax. Mathematically, this constraint is equivalent to the unconfined compression of a sample with mirror images of the volume in the x-, y- and z-directions. In the analysis of the mechanical response of the porous medium, only the solid skeleton is considered and the pores are treated as filled with a very compressible fluid (air), which is consistent with the experimental testing conditions reported in the previous section. Within the solid phase, quartz grains and calcite cement may play different roles. However, for sake of simplification in this study, we consider them as a solid mixture of unknown properties to focus on the influence of the pore-structure. Only the elastic response was analysed for these models in the determination of the RVE. Thus, only the elastic parameters of the solid skeleton are required as input. The magnitude of the loading is small enough to ensure only elastic deformation occurs. The RVE determined by convergence of the asymptotic homogenization is an elastic RVE. Based on the results reported by Regenauer-Lieb et al. (2010, 2013), we assume that the size of the RVE may also be a good first guess for plastic yielding. We test the validity of this assumption a posteriori in the next section. We emphasize that in this contribution only the onset of plastic yielding is considered. The post-yield hardening law will be the subject of future studies. For the evaluation of asymptotic homogenization trends eighteen cubes are cropped with a side-length L = 400, 320, 280, 240, 200, 180, 160, 140, 120, 100, 90, 80, 70, 60, 50, 40, 30 and 20. Finite element meshes are created for these cubes and two types of element are used—hexahedra or tetrahedra. For small volumes (100-cube and smaller), hexahedral elements are used, where one voxel of solid in the image data is converted to one hexahedral element in the finite element mesh. For all other volumes, tetrahedral elements are used in order to reduce computation time. The meshing is implemented by using a large number of triangles to represent the iso-surfaces separating the pore space from the solid. The tetrahedra are then generated from these triangles inward the solid phase. The Avizo® software package (https://www.fei.com/software/avizo3d/) is used to create these tetrahedral meshes. Elements created in this way have similar size and the size is controlled by the size of the iso-surface triangles specified in Avizo. The input elastic parameters of the solid phase (calcite-cemented quartz grains) are first taken to be: Young's modulus E = 100 GPa, Poisson's ratio $$\upsilon$$ = 0.2. At this stage, we do not claim that these values are representative of the actual elastic parameters of the minerals present in this sample. They are only used as reference values for the minerals in order to compare the effective (up-scaled) elastic properties of the porous medium for a series of models. The validity of these values of Young's modulus and Poisson's ratio of the solid phase will be discussed further in the next section. Abaqus® Explicit is used to perform the required finite element computations. The motivation underlying the use of an explicit solver is given in Liu et al. (2016) and is recalled in the discussion (Section 4). Fig. 4(a) shows the evolution with the model size of the computed elastic parameters using the upper and lower bounds. We see that the up-scaled Young's and shear moduli are lower for the stress loading boundary condition (lower bound) than for displacement boundary condition (upper bound). For models with L ≥ 240, stress and displacement loads give reasonably stable and convergent results. Overall, we notice that the stress boundary condition yields less stable predictions (larger oscillations) with increasing volume size than the displacement boundary condition. This instability (oscillations) seems to settle for L > 320. Figure 4. View largeDownload slide Asymptotic homogenization through upper and lower bounds for models of different sizes. (a) Elastic modulus and shear modulus; (b) Poisson's ratio. Figure 4. View largeDownload slide Asymptotic homogenization through upper and lower bounds for models of different sizes. (a) Elastic modulus and shear modulus; (b) Poisson's ratio. The average values of upper and lower bounds for L = 400 can be used as the asymptotic, convergent predictions, that is Young's modulus of 61.5 GPa and shear modulus of 28.2 GPa. The relative errors in the computed moduli for the various volumes and for both boundary conditions are listed in Table 1. The relative errors for all volumes lower than L = 400 are calculated by assuming that the average values for L = 400 (E = 61.5 GPa and μ = 28.2 GPa) are the true predictions. In view of the relative errors obtained, the RVE size can reasonably be taken to be L = 240. Table 1. Computed Young's and shear moduli for models with L ≥ 200, and their relative error compared to the asymptotic and convergent values obtained at L = 400. Size L  Upper bound E  Lower bound E  Upper bound μ  Lower bound μ    Value  Error (per cent)  Value  Error (per cent)  Value  Error (per cent)  Value  Error (per cent)  200  59.22  –3.70  52.48  −14.66  26.93  –4.58  23.42  –17.01  240  62.89  2.27  55.63  –9.53  28.96  2.59  25.26  –10.52  280  64.51  4.90  59.32  –3.53  29.60  4.88  26.77  –5.17  320  64.49  4.87  52.76  −14.20  29.86  5.79  24.18  –14.34  400  62.93  2.34  60.05  –2.34  28.93  2.49  27.52  –2.49  Size L  Upper bound E  Lower bound E  Upper bound μ  Lower bound μ    Value  Error (per cent)  Value  Error (per cent)  Value  Error (per cent)  Value  Error (per cent)  200  59.22  –3.70  52.48  −14.66  26.93  –4.58  23.42  –17.01  240  62.89  2.27  55.63  –9.53  28.96  2.59  25.26  –10.52  280  64.51  4.90  59.32  –3.53  29.60  4.88  26.77  –5.17  320  64.49  4.87  52.76  −14.20  29.86  5.79  24.18  –14.34  400  62.93  2.34  60.05  –2.34  28.93  2.49  27.52  –2.49  View Large 3.2 Drucker-Prager plasticity 3.2.1 Computational scheme We now proceed with the verification of this REV at yield. To model rock yield, we select the pressure-dependent Drucker-Prager criterion in which cohesion and friction angle can be used as intrinsic rock characteristics. A linearized version of the Drucker-Prager criterion is used  $$\ {\sigma _y} = {\sigma _n}{\rm{\ tan}}\beta + c,$$ (1)where σy is the yield stress, σn is the mean stress, c is the cohesion and β is the friction angle. Using eq. (1) c and β can be determined with at least two couples of σy and σn. Two cases are computed for the mechanical RVE in order to determine the cohesion and the friction angle of the rock. For both cases, the following boundary conditions are prescribed (coordinate frame refers to Fig. 1): (1) a zero normal displacement constraint on surfaces x = xmin, y = ymin and z = zmin; (2) the same displacement loading on z = zmax. For Case 1, the lateral surfaces x = xmax and y = ymax are free (similar to the elastic case described in Section 3.1); and for Case 2, a zero normal displacement is prescribed on these lateral surfaces x = xmax and y = ymax. The latter case corresponds to the uniaxial strain loading condition, which also corresponds to the experimental testing conditions reported in Section 2.2. For each case, the stress-strain relationship is computed and the yield stress σy can be obtained. The loading on z = zmax is numerically simulated up to a total strain of 3.5 per cent in the z (axial) direction. 3.2.2 Input parameters: Properties of the solid phase Similar to the computations of elasticity in the determination of the RVE, we need to specify the input parameters of the solid phase (minerals) in order to compute the overall plasticity parameters of the mechanical RVE. In addition to the elastic modulus and Poisson's ratio, two more parameters, yield stress and friction angle, are required as input for the Drucker-Prager plasticity in Abaqus computations. The solid phase (without pores) in our model is a mixture of quartz and calcite for which no common values are found in the literature. Often, values previously reported in the literature relate porous samples rather than to the pure solid mixture (e.g. Airey 1993). Thus, we investigate the parameters by fitting experimental data. To this end, we use the numerical model of Case 2 along with the corresponding experimental data reported in Fig. 2, that is uniaxial strain conditions. For the computed stress-strain relationships, stress is taken to be the differential stress averaged over all elements in the model; the strain is taken to be the relative shortening in the z-direction. The fitting procedure is done by trial-and-error and gradually narrowing the value ranges of the parameters. The Young's modulus is the first fitted as it is independent from the other two parameters, yield stress and friction angle are iteratively narrowed down. The detailed fitting results are given in the Appendix. We conclude that the parameters of the solid phase of the synthetic sandstone, which is composed of calcite-cemented quartz grains, are: Young's modulus between 9.1 and 9.3 GPa; yield stress between 35 and 38 MPa; and friction angle between 52° and 55°. Young's modulus appears to be much lower than values generally accepted for crystalline quartz, that is 70–100 GPa, or crystalline calcite, that is ∼86 GPa (Chen et al.2001). The value of Young's modulus derived for our synthetic sandstone is closer to that of the grey calcite marble (e.g. Martineze-Martineze et al.2012), or that of limestones (e.g. Yasar & Erdogen 2004). This suggests that the stiffness of the calcite-cemented sandstone is mostly controlled by the cement. The friction angle for the pure solid phase (no porosity) appears to be higher than the values often reported for porous sandstones or confined sand packs (Handin 1969; Jaeger & Cook 1976; Twiss & Moores 1992). It is however close to values reported for quartzite with a low to very low porosity (e.g. 48° for the friction angle of Sioux quartzite after Goodman 1989). As no unique and simple value for a given parameter is derived from fitting procedures (see the Appendix), for convenience, when we discuss these results we refer to a representative value for the corresponding range. In the following computations, a Young's modulus of 9100 MPa, a yield stress of 38 MPa and a friction angle of 52° are used as the properties of the solid phase (mineral assemblage composing the synthetic sandstone). 3.2.3 Output parameters: Properties of the porous rock The properties of the porous rock can be investigated with the fitted solid skeleton parameters. Fig. 5(a) represents the stress-strain relationship obtained for Case 1 and Case 2. The curves are plotted using the maximum compressive principal stress (average of all elements) and overall strain (relative shortening in the z-direction). For Case 1 the rock behaves in a typical elastic perfectly plastic way, whereas for Case 2 it exhibits a plastic hardening behaviour. Note that Case 1 represents a uniaxial stress condition with zero radial stress, which typically leads to a brittle failure of the actual rock beyond a peak stress point. Only a small amount of plastic deformation is expected before failure in this case. Because the model is not designed to simulate the brittle failure of rocks, only the maximum value of the stress reached during this simulated test is relevant for our analysis, that is ∼26 MPa. The plateau observed beyond a strain of ∼0.006 is not expected to replicate accurately the actual post-peak response of the actual rock in this brittle case. Compared to the yield stress of the solid skeleton (i.e. 38 MPa), the overall porous structure exhibits a yield stress about 30 per cent lower in Case 1. Figure 5. View largeDownload slide Stress-strain relationships (a) and mean stress—Mises stress relationships (b) of plastic response for the volume of L = 240. Figure 5. View largeDownload slide Stress-strain relationships (a) and mean stress—Mises stress relationships (b) of plastic response for the volume of L = 240. Fig. 5(b) shows the mean stress and the von Mises stress averaged over all elements in the model. Note that the von Mises stress and the mean stress after yielding correspond to σy and σn in eq. (1), respectively. For Case 1, the relationship between the von Mises stress and the mean stress is linear before yielding; both the von Mises stress and the mean stress appear constant after yielding. For Case 2, the von Mises stress and the mean stress show two different linear segments: one at the beginning of the deformation and another after yielding. The trend line of the points after yielding is connected to the clustered points obtained after yielding for Case 1. The slope and the intercept of this trend line are 0.4533 and 20.8, respectively. Fitting eq. (1) to this trend line, where Sy represents σy and Sn represents σn (see Fig. 5b) gives a cohesion of 20.8 MPa and a friction angle of 24.4°. These values are typical for porous sandstones (e.g. Goodman 1989) 3.3 Effect of the volume size An assumption used in this paper is that the size of the mechanical RVE determined from the asymptotic computational homogenization of elastic parameters is also a suitable first guess RVE for plastic yield. A more rigorous way to define the size of RVE would be to perform upper bound and lower bound plastic computations of all the series of models. Owing to the enormous computational cost this is not a very practical suggestion. We therefore developed the above alternative methodology for simplifying the workflow for future work. In order to prove the validity of the approach we performed plastic computations for very large volumes (L ≥ 240) and compare the so-derived plastic parameters. We postulate that our alternative method is robust if we obtain similar results. This would imply that mechanical RVEs can be derived by running simpler elastic analyses. Figs 6(a) and (b) plot the stress-strain relationships of models L = 240, 280, 320 and 400 for Case 1 and Case 2, respectively. Case 1 (Fig. 6a) shows an almost elastic perfectly plastic response for L < 320, and an emerging peak behaviour for L = 320 and 400, with a non-zero residual strength. Experimentally, in unconfined tests (Case 1) at room temperature, porous sandstones usually fail in a brittle way beyond the peak stress, while the stress drops to nearly zero. In this process, strain localisation occurs in the form of shear fractures. Damage prior to failure could be the dominant dissipative mechanism. Therefore, the present model is not expected to be adequate for this case as soon as strain localisation occurs: the RVE determined from elastic considerations becomes essentially irrelevant, plastic deformation localises within the localisation feature leaving the rest of the sample deforming essentially elastically. Figure 6. View largeDownload slide Stress-strain relationships of models L = 240, 280, 320 and 400 for Case 1 (a) and Case 2 (b). Figure 6. View largeDownload slide Stress-strain relationships of models L = 240, 280, 320 and 400 for Case 1 (a) and Case 2 (b). Case 2 (Fig. 6b) exhibits a plastic hardening behaviour for L = 240–400, and all curves are very similar. This confirms a posteriori the validity of the volume size L = 240 used in this case. Similar to Fig. 5(b), we plot the post-yield linear trends for Case 1 and Case 2 for the volumes L = 240, 280, 320 and 400 and fit the linearized Drucker-Prager criterion. We obtain a very narrow range of intercepts ranging from 20.8 to 21.8 MPa. The range of friction angles obtained is also quite narrow, that is 24.4°, 25.3°, 25.5° and 23.6° for L = 240, 280, 320 and 400, respectively. Based on the above results, we conclude that the mechanical RVE determined for elastic properties is also suitable for determining plastic yield parameters, provided that uniaxial strain experimental data (Case 2), as well as high-resolution CT-scan images of the sample, are available to calibrate the model. With this result, we anticipate that the experimental requirements, computational effort and hence the numerical costs can be reduced for similar studies of material characterisation in the future. 3.4 Effects of the porosity In order to investigate the influence of porosity on the yielding behaviour of the rock, we need several models with different porosity. The microstructures of these derivative models should be self-similar to ensure that the structure effect can be ruled out. To this end, we digitally create derivative models from the imaged sandstone sample by self-similarly inflating/deflating the pore structure. The shrinking procedure consists of moving the boundaries separating a pore and the surrounding solid towards the pore by one voxel in the x-, y- and z-directions. In contrast, the expansion procedure consists in moving this boundary towards the solid. This technique was initially developed to detect the percolation threshold of the imaged microstructure of a rock (Liu & Regenauer-Lieb 2011). As these derivative models have very similar microstructures but different porosity, they are ideally suited for investigating the influence of porosity on mechanical properties. The volume L = 240 (the elastic-plastic RVE derived earlier) is used in the following. In addition to the original model, four shrunk derivative models and one expanded model are derived and analysed (see Table 2). The porosity of these six models ranges from 8 to 28 per cent, which covers the most common range of porosity for typical reservoir rocks. The structure of all six models can be seen in Fig. 7. Figure 7. View largeDownload slide Structure of the six models for investigating the effect of porosity on Drucker-Prager plasticity. Pore space is reduced starting from (a) Md_P + 1, to (b) Md_P0, to (c) Md_P-1, to (d) Md_P-2, to (e) Md_P-3, to (f) Md_P-4. Figure 7. View largeDownload slide Structure of the six models for investigating the effect of porosity on Drucker-Prager plasticity. Pore space is reduced starting from (a) Md_P + 1, to (b) Md_P0, to (c) Md_P-1, to (d) Md_P-2, to (e) Md_P-3, to (f) Md_P-4. Table 2. Models used for investigating the effect of porosity on plasticity. Name  Porosity (per cent)  Description  Md_P + 1  28.14  One step of expanding of pores from original model  Md_P0  23.57  Original model  Md_P-1  19.14  One step of shrinking of pores from original model  Md_P-2  15.06  Two steps of shrinking of pores from original model  Md_P-3  11.39  Three steps of shrinking of pores from original model  Md_P-4  8.23  Four steps of shrinking of pores from original model  Name  Porosity (per cent)  Description  Md_P + 1  28.14  One step of expanding of pores from original model  Md_P0  23.57  Original model  Md_P-1  19.14  One step of shrinking of pores from original model  Md_P-2  15.06  Two steps of shrinking of pores from original model  Md_P-3  11.39  Three steps of shrinking of pores from original model  Md_P-4  8.23  Four steps of shrinking of pores from original model  View Large We have already presented the Drucker-Prager plastic computations of the original model of L = 240 in the previous section. Similar computations are conducted for Case 1 and Case 2 for the five derivative models. The stress-strain relationship for the six models are similar to the original model shown in Fig. 5(a)—an apparent elastic perfectly plastic behaviour is observed for Case 1, and a plastic hardening behaviour is observed for Case 2. The only difference for the models is: for a larger porosity, the plateau of the elastic perfectly plastic response is lower, and the plastic hardening is decreasing. Fig. 8 illustrates a few results obtained by fitting the linearized Drucker-Prager yield criterion to the linear trends obtained for Case 1 and 2. Fig. 8(a) corresponds to the model Md_P-4 with the lowest porosity (8.23 per cent). The value of von Mises stress post-yield is around 30 MPa for Case 1 and reaches over 150 MPa for Case 2. The slope of the fitting line is 0.9717 that corresponds to a friction angle of 44.2°. The cohesion for this model is 22.8 MPa corresponding to the intercept. Fig 8(b) shows the model Md_P-2 with a porosity of 15.06 per cent, the fitted friction angle is 35.6° and cohesion is 22.4 MPa; and Fig. 8(c) shows the model Md_P + 1 with the largest porosity of 28.1 per cent, the fitted friction angle is 19.6° and cohesion is 19.4 MPa. Fig. 8. View largeDownload slide Mean stress versus von Mises stress showing the fit of the linear Drucker-Prager relationship for describing the effects of porosity of models, (a) Md_P-4 (ϕ = 8.23 per cent); (b) Md_P-2 (ϕ = 15.06 per cent); (c) Md_P + 1 (ϕ = 28.14 per cent). Fig. 8. View largeDownload slide Mean stress versus von Mises stress showing the fit of the linear Drucker-Prager relationship for describing the effects of porosity of models, (a) Md_P-4 (ϕ = 8.23 per cent); (b) Md_P-2 (ϕ = 15.06 per cent); (c) Md_P + 1 (ϕ = 28.14 per cent). Using the data of all six models, we plot the relationships between the plasticity parameters and the porosity in Fig. 9. The yield stress, the cohesion and the friction angle decrease in an essentially linear way with increasing porosity. The yield stress and the friction angle for a porosity of 0 are 39 MPa and 54.4°, matching the input parameters of the solid phase. We notice that the friction angle decreases more significantly than cohesion with increasing porosity. Figure 9. View largeDownload slide Plastic parameters versus porosity. Figure 9. View largeDownload slide Plastic parameters versus porosity. 4 DISCUSSION We have documented in this paper a new approach for assessing the effects of microstructure on key mechanical properties of rocks. The approach differs from conventional micromechanical analyses of plasticity in that it is not restricted to randomly created mathematical composite models with specific geometries and volume fractions. We have instead studied the plastic response of a synthetic calcite-cemented porous sandstone (CIPS) based on the 3-D images of its microstructures by X-ray microtomography. Although the CIPS is a relatively simple synthetic rock simulating an idealized rock, its inherent complex topology and connectivity have significant impact on the macroscopic properties. It is a step further toward complexity compared to a digitally generated rock microstructure. The computational and meshing aspects and additional complexities of using X-ray microtomographic images of rocks are discussed in detail in Liu et al. (2016). Computational efficiency and convergence issues were solved by selecting a parallel explicit finite element technique (Abaqus® Explicit) which provided numerically acceptable results. The use of an implicit solver such as Abaqus® Standard on such complex microstructures was unsuccessful as the computations did not converge. Numerical issues were encountered through mesh sensitivity affecting the final results of the computations. For instance, the difference in the resulting yield stress for different discretization (hexahedral versus tetrahedral meshes at different resolutions) may be 10 per cent. This error is much higher than for the estimation elastic modulus as reported by Liu et al.2016 (Table 2). Meshing error mainly occurs at the initial step of the meshing process because stepped voxel surfaces must be replaced by triangular iso-surfaces. Doing so, critical connections between grains may break or be reinforced, and these connections significantly affect the mechanical response due to possible stress/strain concentration in these areas (see Fig. 10). With such additional complexity, it is difficult to estimate the extent of the error caused by the meshing. In principle, the finer mesh, the more accurate results. In this study, we did not fully investigate this specific source of error. However, we checked the model with the coarsest mesh, then re-meshed more finely and recomputed all the models listed in Figs 5–9 (and in the Appendix) using the exact same finer mesh, that is the maximum size of the triangles was the same as in Avizo® to generate the required tetrahedral elements. Using the same mesh resolution is crucial to obtain overlapping curves in Fig. 6 and clearly linear relationships in Fig. 9. Furthermore, the intercepts at zero porosity of the fitted yield stress and of the angle of friction lines match the input skeleton parameter reasonably well. This suggests that the discrete error in this study did not systematically overestimate or underestimate the trend of plastic parameters. Figure 10. View largeDownload slide Equivalent plastic strain on two typical slices of the model L = 240 in Case 2. Plastic strain concentrates at grain boundaries. Figure 10. View largeDownload slide Equivalent plastic strain on two typical slices of the model L = 240 in Case 2. Plastic strain concentrates at grain boundaries. For the purpose of the current pilot study models with <10 per cent error were assumed acceptable. Future work should focus more on minimizing numerical artifacts through comparison of different element types (e.g. homogeneous or adaptive cubes, homogeneous or adaptive tetrahedra). Since the calculations are performed on several millions 3-D finite elements a comparison of the resulting accuracy in light of computational costs is necessary. The RVE size was determined by asymptotic computational homogenization deriving the upper and lower bounds of finite element analyses. Owing to the high computation cost this technique was only practical for deriving the size of the elastic RVE. This size of the mechanical RVE is strictly only valid for elastic predictions. We proposed that the elastic RVE is also a good first guess of the RVE for the plasticity calculations. The analysis of the effect of volume size in Section 4.1 showed a posteriori that volume sizes from 240-cube to 400-cube (voxels) have almost exactly the same plastic response. This implies that the RVE size determined by the upper and lower bounds elastic computations is acceptable for plastic yielding. The error in this assumption is less than error introduced by the numerical artefacts of meshing and explicit solution technique (<10 per cent). For obtaining more accurate results future work should consider a full asymptotic derivation of the plastic RVE. This could become possible when the above described numerical resolution and convergence issues are fully resolved in an implicit solver. Since plastic deformation can propagate to larger scales and post-yield bifurcations can introduce new scaling lengths through introduction of the width of shear bands we only discuss initial plastic yield in this paper. We emphasize that the RVE size derived from elastic analysis is inappropriate for the later stages of plastic deformation after yielding. This statement applies specifically if strain localization effects are resolved and shear bands or compaction bands are propagating through the sample. While investigating the bulk plastic response of a porous rock, pore space is treated as mechanically inert in this study. The choice of the input parameters of the solid phase is crucial for the accuracy of the resulting bulk properties. For the specific sandstone sample studied here, the skeleton is made of calcite-cemented quartz grains. There is not experimental data reporting the mechanical parameters of this specific combination of minerals. By fitting the experimental curve of the bulk response of the rock to uniaxial strain loading we were able to estimate the parameters required for the solid phase: Young's modulus was estimated between 9100 and 9300 MPa and the yield stress between 35 and 38 MPa. The friction angle was bracketed between 52° and 55°. As expected the value of Young's modulus obtained is lower than the commonly accepted value for pure quartz. This implies that composite quartz-calcite skeleton appears to have a much lower elastic modulus than the pure quartz or calcite constituents. The friction angle of the solid is higher than that usually reported for sandstones, but close to that reported for (low porosity) quartzites. The output bulk parameters obtained in this study are all within reasonable ranges. Quartz and calcite have different mechanical parameters (Ahrens 1995). We modelled this synthetic sandstone as two-phase medium comprising the pore phase and the solid phase. The latter is considered as a mixture of two solids (quartz and calcite), the properties of which are a priori unknown. We therefore focused on the impact of the pore space, as derived from the 3-D X-ray microtomographic images, on the macroscopic mechanical behaviour of the rock, especially the post-elastic domain of deformation. Previous studies of the plastic behaviour of rocks in relation to their microstructure have been reported by (Liu et al.2012, 2015). In this study, we focus specifically on the impact of the size of the RVE and of the porosity on the macroscopic yield parameters of the rock assuming that the solid phase is a mixture of quartz and calcite. This is meant to be a first step for benchmarking future studies in which the arrangement and respective properties of the constituting minerals (quartz and calcite) could then be added. However, this study already suggests (see Fig. 10) that the macroscopic deformation of the porous rock is strongly controlled by the cemented grain contacts where high strain values occur. Interesting results of this paper are the simple linear relationships providing quantitative insights into the effect of microstructure on macroscopic rock properties. A good example is the reduction of the friction angle to with increasing porosity as shown in Fig. 9. The results show that the yield stress, the cohesion and the friction angle of the rock linearly decrease with increasing porosity, for a range of porosity comprised between 8 and 28 per cent. This range covers the porosity values encountered for most conventional reservoir rocks. Because the macroscopic behaviour of the rock is found to be strongly controlled by its microstructure, and since the models studied here are geometrically self-similar, we do not recommend to extrapolate the relationships derived here to models with very different structures. More work is required in that case. It is not recommended either to extrapolate the linear relationships obtained here further to an untested range. According to percolation theory, parameters such as yield stress change exponentially when the solid volume fraction approaches the percolation threshold (Stauffer & Aharony 1994; Sahimi 1998). Based on the linear trends obtained in Fig. 9, it appears that the volume fractions of the solid in the models analysed here are far from the percolation threshold. According to the basic rule of mixture (Voigt 1889; Reuss 1929; Alger 1997), the property Ec of a composite material made of two phases should be comprised between an upper bound Ec = fEf + (1 − f)Em and a lower bound $${E_c} = {( {\frac{f}{{{E_f}}} + \frac{{1 - f}}{{{E_m}}}} )^{ - 1}}\$$, where f is the volume fraction of the inclusions, and Em and Ef are the properties of the matrix and the inclusions, respectively. In our two-phase model (solid and pores), f corresponds to porosity, and pores are voids with Ef = 0. We obtained a linear relation between the parameters of the solid phase and the macroscopic properties. Note that the expression of the upper bound leads to a similar linear relation with a slope of (1 − f), which differs from the slope we found in Fig. 9. Because the upper and lower bounds have been derived for the linear elastic deformation regime, this difference is not surprising. 5 CONCLUSIONS In this pilot study, we have clearly shown a promising new path for evaluating rock mechanical properties from micro-CT imaging. The analysis has been deliberately chosen to be close to the ideal mathematical rock models by the choice of imaging a synthetic rock composite. We have shown that the suggested computational homogenization workflow is suitably robust and reliable within a 10 per cent error of the mechanical properties. While the described computational homogenization method is well established (Terada et al.2000), we have introduced three new elements: (1) generalization of the classical computational homogenization methods in terms of thermodynamic bounds; (2) proposed equivalence of the elastic RVE to a plastic RVE at initial yield and (3) calculation of quantitative microstructural trends derived from derivative models (inflated/deflated microstructures). The successful test of these new improvements on a synthetic porous sandstone opens promising avenues for future assessments of the effect of microstructure on mechanical properties for complex microstructures encountered in rock mechanics. Future development includes the consideration of multiphase involved in a sample, such as from quartz, calcite and pores to quartz, calcite and saturated/unsaturated fluid and the consideration of other constitutions, such as viscoplasticity and elasto-viscoplasticity at different temperature and pressure. The approach developed here and its upgraded version will ultimately lead to new ways of characterizing rock properties while relying less on time consuming and costly laboratory measurements, with potential application in crustal strength and deformation in geophysics, borehole stability, hydraulic fracturing, or reservoir compaction studies in the oil and gas, geothermal or CO2 geo-sequestration industries. Acknowledgements We thank Keyu Liu, China University of Petroleum (East China), formerly CSIRO Energy, provided the synthetic sandstone sample for experiments and the micro CT images. Jie Liu would like to acknowledge the support of National Natural Science Foundation of China grants no. 41574087. Part of the computing was performed on Tianhe-2 supercomputer and the computing time was supported by the Special Program for Supercomputing of the NSFC-Guangdong Joint Fund. Klaus Regenauer-Lieb would also like to acknowledge support of the Australian Research Council through ARC Discovery grants no. DP140103015, DP170104550, DP170104557. REFERENCES Ahrens T.J., 1995. Mineral Physics and Crystallography: A Handbook of Physical Constants , American Geophysical Union. Google Scholar CrossRef Search ADS   Airey D.W., 1993. Triaxial testing of naturally cemented carbonate soil. J. Geotech. Eng. , 119( 9), 1379– 1398. Google Scholar CrossRef Search ADS   Alger M.S.M., 1997. Polymer Science Dictionary.  2nd edn. Springer Publishing. ISBN 0412608707. Almqvist B.S.G., Mainprice D., Madonna C., Burlini L., Hirt A.M., 2011. Application of differential effective medium, magnetic pore fabric analysis, and X-ray microtomography to calculate elastic properties of porous and anisotropic rock aggregates, J. geophys. Res.: Solid Earth , 116, 1– 17. Google Scholar CrossRef Search ADS   Arns C.H., Knackstedt, M.A., Pinczewski W.V., Lindquist W.B., 2001. Accurate estimation of transport properties from microtomographic images, Geophys. Res. Lett. , 28, 3361– 3364. Google Scholar CrossRef Search ADS   Arns C.H., Knackstedt M.A., Pinczewski W.V., Garboczi E.J., 2002. Computation of linear elastic properties from microtomographic images: methodolgy and agreement between theory and experiment, Geophysics , 67( 5), 1396– 1405. Google Scholar CrossRef Search ADS   Arns C.H., Knackstedt M.A., Martys N.S., 2005. Cross-property correlation and permeability estimation in sandstone, Phys. Rev. E , 72, 046304, doi:10.1103/PhysRevE.72.046304. Google Scholar CrossRef Search ADS   Blunt M.J.,et al.  , 2013. Pore-scale imaging and modelling, Adv. Water Resour. , 51, 197– 216. Google Scholar CrossRef Search ADS   Charalambakis N., 2010. Homogenization techniques and micromechanics. A survey and perspectives, Appl. Mech. Rev. , 63( 3), 030803. Google Scholar CrossRef Search ADS   Chen C.C., Lin C.C., Liu L.G., Sinogeijin S.V, Bass J.D., 2001. Elasticity of single-crystal calcite and rhodochrosite by Brillouin spectroscopy, Am. Mineral. , 86, 1525– 1529. Google Scholar CrossRef Search ADS   Chen C., Packman A.I, Zhang D., Gaillard J.-F., 2010. A multi-scale investigation of interfacial transport, pore fluid flow, and fine particle deposition in a sediment bed, Water Resour. Res. , 46, W11560, doi:10.1029/2009WR009018. Fredrich J.T., DiGiovanni A.A., Noble D.R. 2006. Predicting macroscopic transport properties using microscopic image data, J. geophys. Res. , 111, B03201, doi:10.1029/2005JB003774. Google Scholar CrossRef Search ADS   Goodman R.E., 1989. Introduction to Rock Mechanics , 2. Wiley. Handin J., 1969. On the Coulomb–Mohr failure criterion, J. geophys. Res. , 74, 5343– 5348. Google Scholar CrossRef Search ADS   Jaeger J.C., Cook N.G.W., 1976. Fundamentals of Rock Mechanics , Chapman & Hall, Wiley, 585 pp. Knackstedt M.A.et al.  , 2006. Elastic and transport properties of cellular solids derived from three-dimensional tomographic images, Proc. R. Soc. A , 462, 2833– 2862. CrossRef Search ADS   Liu J., Regenauer-Lieb K., 2011. Application of percolation theory to microtomography of structured media: percolation threshold, critical exponents, and upscaling, Phys. Rev. E , 83, 1– 13. Liu J., Regenauer-Lieb K., Hines C., Liu K., Gaede O., Squelch A., 2009. Improved estimates of percolation and anisotropic permeability from 3-D X-ray microtomography using stochastic analyses and visualization, Geochem., Geophys., Geosyst. , 10, doi:10.1029/2008GC002358. Liu J., Freij-Ayoub R., Regenauer-Lieb K., 2012. Determination of the plastic strength of carbonates from microtomography and upscaling using percolation theory, in Proceedings of the ECCOMAS 2012 , ed. Eberhardsteiner J.et al.  , Vienna, Austria, September 10–14, 2012. Liu J., Freij-Ayoub R., Regenauer-Lieb K., 2015. Rock plasticity from microtomography and upscaling, J. Earth Sci. , 26( 1), 053– 059. Google Scholar CrossRef Search ADS   Liu J., Pereira G.G., Liu Q., Regenauer-Lieb K., 2016. Computational challenges in the analyses of petrophysics using microtomography and upscaling: a review, Comput. Geosci. , 89, 107– 117. Google Scholar CrossRef Search ADS   Martinez-martinez J., Benavente D., Angeles G.-C.M., 2012. Comparison of the static and dynamic elastic modulus in carbonate rocks, Bull. Eng. Geol. Environ. , 71, doi:10.1007/s10064-011-0399-y. Nemat-Nasser S., Hori M., 1993. Micromechanics: Overall Properties of Heterogeneous Materials, Elsevier. Regenauer-Lieb K.et al.  , 2010. Time-dependent, irreversible entropy production and geodynamics, Phil. Trans. R. Soc. Lond., A , 368( 1910), 285– 300. Google Scholar CrossRef Search ADS   Regenauer-Lieb K.et al.  , 2013. Multiscale coupling and multiphysics approaches in earth science: theory, J. Coupled Syst. Multiscale Dyn. , 1( 1), 49– 73. Google Scholar CrossRef Search ADS   Reuss A., 1929. Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle, ZAMM - J. Appl. Math. Mech./ Zeitschrift für Angewandte Mathematik und Mechanik , 9, 49– 58. Google Scholar CrossRef Search ADS   Roberts A.P., Garboczi E.J., 2002. Elastic properties of model random three-dimensional open-cell solids, J. Mech. Phys. Solids , 50( 1), 33– 55. Google Scholar CrossRef Search ADS   Sahimi M., 1998. Non-linear and non-local transport processes in heterogeneous media: from long-range correlated percolation to fracture and materials breakdown, Phys. Rep. , 306, 213– 395. Google Scholar CrossRef Search ADS   Shulakova V.et al.  , 2013. Computational elastic up-scaling of sandstone on the basis of X-ray micro-tomographic images, Geophys. Prospect. , 61( 2), 287– 301. Google Scholar CrossRef Search ADS   Shulakova V., Sarout J., Pimienta L., Lebedev M., Mayo S., Clennell M.B., Pervukhina M., 2016. Effect of supercritical CO2 on carbonates: savonnieres sample case study, Geophys. Prospect. , 61, 287– 301. Google Scholar CrossRef Search ADS   Stauffer D., Aharony A., 1994. Introduction to Percolation Theory , Taylor & Francis. Spanne P., Thovert J.F., Jacquin C.J., Lindquist W.B., Jones K.W., Adler P.M., 1994. Synchrotron computed microtomography of porous media: topology and transports, Phys. Rev. Lett. , 73, 2001– 2004. Google Scholar CrossRef Search ADS PubMed  Stock S.R., 2008. Recent advances in X-ray microtomography applied to materials, 53( 3), 129– 181. Terada K., Hori M., Kyoya T., Kikuchi N., 2000. Simulation of the multi-scale convergence in computational homogenization approaches, Int. J. Solid Struct. , 37, 2285– 2311. Google Scholar CrossRef Search ADS   Twiss R.J., Moores E.M., 1992. Structural Geology , W.H. Freeman and Company, 532 pp. Veveakis E., Regenauer-Lieb K., 2015. Review of extremum postulates, Curr. Opin. Chem. Eng. , 7, 40– 46. Google Scholar CrossRef Search ADS   Voigt W., 1889. Ueber die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper, Annalen der Physik , 274, 573– 587. Google Scholar CrossRef Search ADS   Wong T.-F., David C., Zhu W., 1997. The transition from brittle faulting to cataclastic flow in porous sandstones: mechanical deformation, J. geophys. Res. , 102( B2), 3009– 3025. Google Scholar CrossRef Search ADS   Yang S., Liu K., Mayo S., Tulloh A., 2012. CIPS sandstone microstructure, v2, CSIRO Data Collection , doi:org/10.4225/08/5476787A1A50F. Yasar E., Erdogan Y., 2004. Correlating sound velocity with the density, compressive strength and Young's modulus of carbonate rocks, Int. J. Rock Mech. Mining Sci. , 41, 871– 875. Google Scholar CrossRef Search ADS   APPENDIX Here we show the comparisons of experimental curve with computed curves that the ranges of elastic modulus, yield stress and the angle of friction have been preliminarily obtained. Fig. A1 Shows the stress-strain relationships derived for Case 2, at small strain, for three values of Young's modulus, and the comparison with the experimental curve. We see both 9100 and 9300 MPa are almost parallel to the experimental curve—the experimental curve is not starting from zero but at a small value of stress, the computed curve with E = 9100 and 9300 MPa can be seen as the best fit of the experimental curve when strain is less than 0.005. In these computations, the yield stress and friction angle are taken to be 43 MPa and 47°, respectively. These two yield parameters do not affect the deformation response in the elastic domain. Therefore, the elastic modulus of the solid skeleton in this sample is estimated to be between 9100 and 9300 MPa. Figure A1. View largeDownload slide Three computed stress-strain curves of different elastic moduli comparing with experimental measurement to investigate the elastic modulus of the solid. Figure A1. View largeDownload slide Three computed stress-strain curves of different elastic moduli comparing with experimental measurement to investigate the elastic modulus of the solid. Fig. A2 compares the stress-strain relationships using three different yield stress values as input solid parameter with the experimental curve. In these computations, Young's modulus and the friction angle are 9100 MPa and 47°, respectively. We see that when the yield stress of the solid phase is taken to be 35 MPa, the computed stress-strain curve is lower than the experimental curve. However, the trends of the computed and experimental curves are very similar (parallel curves). Overall, a better fit over the wider strain range of 0–0.01 is obtained with a yield stress value of 38–41 MPa. Figure A2. View largeDownload slide Three computed stress-strain curves of different yield stress comparing with experimental measurement to investigate the yield stress of the solid. Figure A2. View largeDownload slide Three computed stress-strain curves of different yield stress comparing with experimental measurement to investigate the yield stress of the solid. Fig. A3 compares the stress-strain relationships from different values of the friction angle. In these computations, we used a Young's modulus of 9100 MPa and a yield stress value of 38 MPa. From the figure, we see that 50° is too low; 52° is slightly lower than the experimental curve in the strain range 0.01–0.03, but replicates quite accurately the yield point; 55° matches the experimental curve very well except for the high values of strain. Therefore, the friction angle of the solid phase is estimated to be constrained between 52° and 55°. Figure A3. View largeDownload slide Three computed stress-strain curves of different friction angle comparing with experimental measurement to investigate the angle of friction of the solid. Figure A3. View largeDownload slide Three computed stress-strain curves of different friction angle comparing with experimental measurement to investigate the angle of friction of the solid. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations