Subsections

# 3.3 Constitutive equation: stress-strain relationships

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:

• (TOP) Hooke's law: the force [N] required to produce an elongation [m] in a spring with mechanical constant [N/m] is

 (3.9)

• (MIDDLE) Consider a prismatic solid bar with cross sectional area and length . The force required to produce an elongation [m] is inversely proportional to , and proportional to proportional to , and (the stiffness modulus of the solid), such that

 (3.10)

which can be expressed as stress and strain

 (3.11)

• (BOTTOM) In general, the stress tensor is proportional to the strain tensor through the stiffness tensor

 (3.12)

## 3.3.1 Linear isotropic elasticity

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)

 (3.13)

The solid will most likely tend to enlarge in the direction perpendicular to the stress applied. The Poisson ratio (greek letter nu) is defined as (-1) times the ratio between the strain perpendicular to the applied stress (or ) and the strain in the direction of the applied stress

 (3.14)

These two coefficients are the two coefficients conventionally used as elasticity constants in continuum mechanics. We will see later that in the subsurface we almost never find conditions of laterally “unconfined” stress loading like the one shown in Figure 3.11.

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:

• A soft mudrock with E = 1 GPa
• A soft sandstone with E = 10 GPa
• A hard limestone with E = 50 GPa

SOLUTION
Let us work in SI units:

psi    MPa    MPa

Then

GPa

GPa

GPa

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.

## 3.3.2 The isotropic solid in Voigt notation

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.

 (3.15)

In addition, shear strains are proportional to the applied shear stress through shear modulus , such that,

 (3.16)

Hence, all six equations permit putting together the shear strain tensor as a function of the stress tensor through compliance fourth-order tensor .

 (3.17)

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:

 (3.18)

For example, let us apply a stress as in example in Figure 3.11. The result of is

which are the same strains we found above in the definition of and .

The inverse of the compliance matrix is the stiffness matrix and let us calculate stress as a function of strain.

 (3.19)

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:

or equivalently

where the Lamé parameters are:

 (3.20)

and

 (3.21)

Notice that , the shear modulus as defined above. Putting equations in all directions together yields the complete set of Lamé's equations:

 (3.22)

EXAMPLE 3.4: Write the Lamé equations (Eq. 3.22) in matrix format using the Voigt notation.

SOLUTION

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.

## 3.3.3 Effective stress and elasticity

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:

 (3.23)

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.

 (3.24)

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.

## 3.3.4 Calculation of horizontal stress according to linear elasticity

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.

 (3.25)

The effective vertical stress will be

 (3.26)

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

 (3.27)

Let us express as a function of , and plug it in the equation for horizontal stresses and . The result is

 (3.28)

or equivalently

 (3.29)

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

 (3.30)

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,

 (3.31)

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,

 (3.32)

Let us now substitute in the equations of and as a function of . The result is:

 (3.33)

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:

 (3.34)

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:

1. Calculate total vertical stress
2. Calculate pore pressure .
1. If hydrostatic, use psi/ft.
2. If non-hydrostatic, use the method based on shale porosity or use directly if is given.
3. Calculate effective vertical stress:
4. Calculate effective horizontal stresses and :
1. Assuming linear isotropic elasticity, use Eq. 3.34 if tectonic strains are not zero or simply Eq. 3.29 if tectonic strains are zero.
2. Assuming subsurface stresses are at yield, use Eqs. in the Chapter about “Stresses on Faults”.
3. Assuming subsurface stresses are affected by visco-elastic response, add stress relaxation component (Fig. 3.26)
5. Calculate total horizontal stress by adding pore pressure to effective horizontal stresses:

and

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 .

SOLUTION
At a depth of 7950 ft

Hence, the effective vertical stress is

Now, we are in conditions of using Eq. 3.34. Let us first calculate the plane strain modulus:

Then,

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 .

## 3.3.5 Calculation of reservoir compressibility with linear elasticity

The rock pore volume compressibility is a critical parameter in the fluid flow mass conservation equation and therefore on the diffusivity equation (1D example):

 (3.35)

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:

 (3.36)

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

 (3.37)

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:

 (3.38)

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,

 . (3.39)

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.

SOLUTION
The constrained modulus is

Hence,

## 3.3.6 Generalized linear elasticity problem

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:

 (3.40)

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.