Constitutive equations tell us how a solid deforms (in time) as a response to stresses, to changes of temperature and to changes of pore pressure among others. How to choose a constitutive equation depends on the material properties, the magnitude of strain changes, the magnitude of stresses, and the loading rate among other factors.
The simplest constitutive relationship for solids is linear elasticity, in which stresses and strains are linearly related by constant coefficients. The examples in Figure 3.10 correspond to applications of linear elasticity in various dimensions:
Consider a prismatic solid with length to which we apply a stress on top face 3 (Figure 3.11). The bottom face is not allowed to move in direction 3 but it can slide sideways. The four other faces are free to move in all directions. Notice that the top face can also deform in directions 1 and 2. The Young modulus is defined as the ratio between the applied stress and the resulting strain (in the direction of the applied stress)(or ) and the strain in the direction of the applied stress
The real behavior of rocks differs from the linear elastic assumption. Figure 3.12 shows a schematic representation of a typical unconfined loading test. The figure plots axial stress in the vertical axis and axial strain in the horizontal axis. Often, rock plugs are not perfectly parallel or may have some microcracks. Both features make the initial loading stress-strain behavior look less stiff than the actual rock stiffness. After the initial loading, the rock may show a linear response -where the Young modulus is measured- followed by softening approaching rock failure and the peak stress. When the test is performed under unconfined conditions, the peak stress is termed the “unconfined compressive strength (UCS)” of the rock (further explained in Section 4). The Poisson ratio can be measured in the same range of the measurement of when lateral strain transducers are available.
The Young modulus of sediments and rocks varies widely. Figure 3.13 shows typical values of Young's modulus.
EXAMPLE 3.3: Compute the (axial) strain expected for a rock subjected to 3,000 psi of (axial) stress under unconfined axial loading for:
Let us work in SI units:
Notice that rocks can be quite stiff and even for an effective stress as large as 3,000 psi (equivalent to a depth onshore of 5,000 ft under hydrostatic pore pressure), the deformation is in the order of 1% to 0.1% or less.
A generalization of the Young's modulus and Poisson's ratio equations (Eq. 3.13 and 3.14) in all directions leads to the 3 independent equations.
In addition, shear strains are proportional to the applied shear stress through shear modulus , such that,
Hence, all six equations permit putting together the shear strain tensor as a function of the stress tensor through compliance fourth-order tensor .
For ease of calculation, we will express the stress and strain tensors as matrices, such that will be a matrix. This notation is called Voigt notation. Hence, fourth-order tensor can be expressed as a matrix:
For example, let us apply a stress as in example in Figure 3.11. The result of is
The inverse of the compliance matrix is the stiffness matrix and let us calculate stress as a function of strain.
Voigt notation is easier to code in computer codes that work with matrices.
The Lamé equations are the same equations shown above but use the Lamé parameters and instead of and . For example, let us write the first equation of the product of the stiffness tensor and the strain tensor in Voigt notation:
where the Lamé parameters are:
Notice that , the shear modulus as defined above. Putting equations in all directions together yields the complete set of Lamé's equations:
EXAMPLE 3.4: Write the Lamé equations (Eq. 3.22) in matrix format using the Voigt notation.
Remember that there are only two independent constitutive parameters in linear isotropic elasticity. The usual pair choice is and . However, there are other options depending on the application and equations used, e.g, and . A complete list of parameter pairs is available scrolling to the bottom in https://en.wikipedia.org/wiki/Young's_modulus. Figure 3.14 list the most common equivalencies.
Porous solids deform and fail due to the application of effective stresses rather than total stress. Hence, Hooke's law requires to use the effective stress tensor rather than the total stress tensor. The equation is incorrect. Instead, the stress-strain relationship requires effective stress:
Pore pressure has an effect on normal stresses only (fluid pressure would not be able to cause solid shear strains). Hence, only pore pressure is subtracted from the diagonal terms of the total stress tensor. The subtracted value is the same in all directions because pore pressure is the same in all directions at a given point location.
Rigorously, the effective stress tensor needs a correction of pore pressure by the Biot coefficient that accounts for solid grain deformation with changes in pore pressure.
For most problems, the assumption of is satisfactory. The rock matrix of tight sandstones and shales may have a Biot coeffiecient as low as . The theory of poroelasticity is covered in the “Advanced Geomechanics” course with a brief introduction in Section 3.7.1.
Let us revisit the problem of stress calculation in a half-space, such as the Earth's shallow subsurface. We already know that the vertical total stress ( ) is a function of depth and rock bulk mass density.
The effective vertical stress will be
Let us now assume that the half space did not deform in horizontal directions ( ), usually known as a “tectonically passive environment”. This means that the solid is laterally contained at “repose” and no additional horizontal strains have been added either compressive or tensile. Such is the case of a sedimentary basin with no additional tectonic strains.
Let us now use Equation 3.23 together with the equilibrium equation. Shear strains are zero. Hence . Then, the multiplication of , results in
Let us express as a function of , and plug it in the equation for horizontal stresses and . The result is
For typical values of , the horizontal stress coefficient is (Figure 3.17). Thus, the effective horizontal stress is approximately one third of the effective vertical stress. Contrary to a fluid, the solid does not push sideways with all its weight. It pushes sideways with just a fraction of its weight proportionally to its tendency to deform sideways, i.e., the Poisson ratio. Notice that implies . An “effective” is applicable for fluids, soft rocks under undrained loading, and salt rocks.
The total horizontal stresses are obtained by adding pore pressure to the effective horizontal stresses: and .
Equation 3.29 allows us to approximate a lower bound for the fracture gradient, that is, the pressure required to open a hydraulic fracture. Such pressure will be equal or greater than the minimum horizontal total stress (assuming zero tectonic strains):
The gradient is the variation of pressure (or stress) with depth, i.e., derivative with respect to depth . Assuming that the material properties are constant, then,
For example, for onshore conditions with typical values psi/ft, psi/ft, and , the fracture gradient is psi/ft. Figure XX shows a schematic example of the calculated fracture gradient.
Let us now relax the assumption of horizontal strains equal to zero, such that they are not zero , but are known quantities. We use the equation again. Shear strains are zero. Hence . The resulting equations are,
Let us now substitute in the equations of and as a function of . The result is:
Horizontal strains are usually caused by tectonic plate movement. Hence, we can call them “tectonic strains” and say that this is a “tectonically active” environment, particularly for large compressive horizontal strains that may increase horizontal stress over the vertical stress. Let us call the maximum (compressive) tectonic strain, and the minimum tectonic strain in a given direction. As a result the maximum effective horizontal stress and minimum horizontal stresses are:
where is called the plane strain modulus. We will see later that the plane strain modulus, rather than the Young's modulus, appears in many of the equations of interest to subsurface applications. These equations have been coded in a Jupyter notebook available at https://mybinder.org/v2/gh/dnicolasespinoza/GeomechanicsJupyter/master?filepath=HorizontalStresses_Widget.ipynb. The above algorithm further assumes a linear increase of strain with depth.
The following workflow is valid to calculate horizontal total stress with any constitutive (rock property) model:
EXAMPLE 3.5: Calculate the total horizontal stresses in a section of the Barnett Shale located at 7,950 ft (TVD) using the theory of linear elasticity. Assume a constant vertical stress gradient MPa/km, overpressure parameter , shale Young's modulus psi, Poisson's ratio , and tectonic strains and .
At a depth of 7950 ft
Finally, let us compute the total horizontal stresses by adding pore pressure:
Horizontal stresses tend to be different in most basins, hence, . In practice, differences MPa (100 psi) tend to be enough to force hydraulic fractures to grow mostly perpendicular to in places subjected to normal faulting or strike slip stress regime. Locations with small horizontal stress anisotropy (mostly normal faulting sites) impose less control on orientations of hydraulic fractures (Chapter 7) and on orientation of faults (Chapter 5). Polygonal faults (Chapter 5 - [update link]) are an example of normal fault growth with strike in all directions because .
The rock pore volume compressibility is a critical parameter in the fluid flow mass conservation equation and therefore on the diffusivity equation (1D example):
Where the total compressibility is . Reservoir simulators usually calculate the fluid compressibility based in phase behavior, hence, the only required input is . For example, compaction drive is proportional to rock compressibility (see https://petrowiki.org/Compaction_drive_reservoirs).
The pore volume compressibility tells us what the change of pore volume is due to a change in pore pressure:
The equation above captures reservoir boundary conditions in which the total vertical stress remains constant (overburden above the reservoir does not change) and there is no change of lateral strain , a condition also termed as “uniaxial strain” deformation. Such condition is appropriate in long and thin reservoirs with a compliant caprock (Figure 3.19).
The measurements of are derived from bulk volume measurements. Let us assume that the change of pore volume is equal to the change of bulk volume , which means that all bulk deformation is caused by change of porosity. Hence it is possible to rewrite the definition of as
Porosity is defined as and the term between parenthesis is defined as the bulk compressibility under uniaxial condition (notice ). Hence, the parameter is linked to the bulk rock compressibility through porosity:
where the bulk compressibility with no lateral strain is approximately equal to the inverse of the bulk constrained modulus , where for an isotropic elastic solid. The approximation is due to a correction needed to account for grain compressibility. Finally, we can calculate the uniaxial strain pore compressibility using the typical mechanical parameters and as,
Unfortunately, the theory of linear elasticity is quite limited to capture the visco-elasto-plastic behavior of rocks upon depletion during long times and with large strains. Hence, Eq. 3.39 is just a first order approximation.
Typical values of pore volume compressibility vary from 2 to 30 psi, where psi sip. Stiff well cemented rocks have low pore volume compressibility sip while uncemented loose sediments tend to have high pore volume compressibility sip.
EXAMPLE 3.6: Calculate the pore compressibility of a rock tested in the laboratory with porosity , Young's modulus GPa, and . Provide the solution in [psi] units.
The constrained modulus is
The general solution of a linear elasticity problem requires combining the equilibrium, kinematic, and constitutive equations. The result is a differential equation with displacement as the unknown:
where is the first Lamé parameter, is the shear modulus, is the rock bulk mass density, and is the body force acceleration vector (usually gravity). A review of the gradient , divergence and Laplacian operators is available at https://en.wikipedia.org/wiki/Vector_calculus_identities. In summary, these are all derivatives that quantify changes of displacement in space. The full derivation of this equation can be found here: https://youtu.be/1PnQ10H2vV0.
The solution requires knowledge of the domain geometry, boundary conditions and initial conditions. For example, a hydraulic fracture simulator solves numerically these same equations (Figure 3.20). In this class, we will see analytical solutions of this equation for 1) displacements and stresses around wellbores, and 2) displacements and stresses around planar fractures.