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.

Figure 3.9: Example of deformations imparted by an applied stress.

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 $F$ [N] required to produce an elongation $\Delta x$ [m] in a spring with mechanical constant $k$ [N/m] is

  $$F = k \Delta x$$ (3.9)

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

  $$\frac{F}{A} = E \frac{\Delta x}{L}$$ (3.10)

  which can be expressed as stress $\sigma$ and strain $\varepsilon$

  $$\sigma = E \varepsilon$$ (3.11)

• (BOTTOM) In general, the stress tensor $\uuline{\sigma}$ is proportional to the strain tensor $\uuline{\varepsilon}$ through the stiffness tensor $\uuline{C}$

  $$\uuline{\sigma} = \uuline{C} \: \uuline{\varepsilon}$$ (3.12)

Figure 3.10: From Hooke's law to generalized 3D linear elasticity.

3.3.1 Linear isotropic elasticity

Consider a prismatic solid with length $L$ to which we apply a stress $\sigma_{33}$ 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 $E$ is defined as the ratio between the applied stress $\sigma_{33}$ and the resulting strain (in the direction of the applied stress) $\varepsilon_{33}$

$$E = \frac{\sigma_{33}}{\varepsilon_{33}}$$ (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 $\varepsilon_{11}$ (or $\varepsilon_{22}$) and the strain in the direction of the applied stress $\varepsilon_{33}$

$$\nu = -\frac{\varepsilon_{11}}{\varepsilon_{33}}$$ (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.

Figure 3.11: Unconfined stress loading (compression) of a linear elastic isotropic solid. Because the solid is isotropic, the same equations are valid for compression in any other direction, and also in tension.

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 $\sigma_a$ in the vertical axis and axial strain $\varepsilon_a$ 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 $E$ when lateral strain transducers are available.

Figure 3.12: Schematic of stress-strain curve during rock uniaxial loading in the laboratory.

The Young modulus of sediments and rocks varies widely. Figure 3.13 shows typical values of Young's modulus.

Figure 3.13: Typical Young moduli for sediments and rocks. The values correspond to “quasi-static” loading from stress-strain response. SS: sandstone; SiS: Siltstone; Sh: Shale; Gr: Granite.

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:

$$\sigma_a = 3000$$    psi$$= 3000 \frac{1}{145}$$    MPa$$= 20.7$$    MPa

Then

$$\varepsilon_a(1$$    GPa$$) = \frac{\sigma_a}{E} = \frac{20.7 \text{ MPa}}{1000 \text{ MPa}} = 0.0207 = 2.07 \%$$

$$\varepsilon_a(10$$    GPa$$) = \frac{\sigma_a}{E} = \frac{20.7 \text{ MPa}}{10000 \text{ MPa}} = 0.00207 = 0.207 \%$$

$$\varepsilon_a(50$$    GPa$$) = \frac{\sigma_a}{E} = \frac{20.7 \text{ MPa}}{50000 \text{ MPa}} = 0.00041 = 0.041 \%$$

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. $\: \: \blacksquare$

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.

$$\left\lbrace \begin{array}{rcl} \varepsilon_{11} & = & +\cfra... ...\: \sigma_{22} + \cfrac{1}{E} \: \sigma_{33} \end{array}\right.$$ (3.15)

In addition, shear strains $\varepsilon_{ij}$ are proportional to the applied shear stress through shear modulus $G=E/[2(1+\nu)]$, such that,

$$\left\lbrace \begin{array}{rcl} 2 \varepsilon_{12} & = & \cfr... ...silon_{23} & = & \cfrac{1}{G} \: \sigma_{23} \end{array}\right.$$ (3.16)

Hence, all six equations permit putting together the shear strain tensor $\uuline{\varepsilon}$ as a function of the stress tensor $\uuline{\sigma}$ through compliance fourth-order tensor $\uuline{D}$.

$$\uuline{\varepsilon} = \uuline{D} \: \uuline{\sigma}$$ (3.17)

For ease of calculation, we will express the stress and strain tensors as $6 \times 1$ matrices, such that $\uuline{D}$ will be a $6 \times 6$ matrix. This notation is called Voigt notation. Hence, fourth-order tensor $\uuline{D}$ can be expressed as a $6 \times 6$ matrix:

$$%compliance matrix \left[ \begin{array}{c} \varepsilon_{11} ... ...ma_{13} \cfrac{}{}\\ \sigma_{12} \cfrac{}{} \end{array}\right]$$ (3.18)

For example, let us apply a stress $\uuline{\sigma} = [0,0,\sigma_{33},0,0,0]^T$ as in example in Figure 3.11. The result of $\uuline{D} \: \uuline{\sigma}$ is

$$\uuline{\varepsilon} = \left [ -\cfrac{\nu}{E} \: \sigma_{33},-\cfrac{\nu}{E} \: \sigma_{33},\cfrac{1}{E} \: \sigma_{33},0,0,0 \right]^T$$

which are the same strains we found above in the definition of $E$ and $\nu$.

The inverse of the compliance matrix is the stiffness matrix $\uuline{C} = \uuline{D}^{-1}$ and let us calculate stress as a function of strain.

$$%stiffness matrix \left[ \begin{array}{c} \sigma_{11} \\ \s... ... \cfrac{}{}\\ 2 \varepsilon_{12} \cfrac{}{} \end{array}\right]$$ (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 $\lambda$ and $\mu$ instead of $E$ and $\nu$. For example, let us write the first equation of the product of the stiffness tensor and the strain tensor in Voigt notation:

$$\sigma_{11} = \cfrac{E}{(1+\nu)(1-2\nu)} \left[ (1-\nu) \varepsilon_{11} + \nu \varepsilon_{22} + \nu \varepsilon_{33} \right]$$

or equivalently

$$\sigma_{11} = \cfrac{\nu E}{(1+\nu)(1-2\nu)} (\varepsilon_{11} + ... ...\nu E}{(1+\nu)(1-2\nu)} \left( -1 + \cfrac{1-\nu}{\nu} \right) \varepsilon_{11}$$

where the Lamé parameters are:

$$\lambda = \cfrac{\nu E}{(1+\nu)(1-2\nu)}$$ (3.20)

and

$$2\mu = \cfrac{\nu E}{(1+\nu)(1-2\nu)} \left( -1 + \cfrac{1-\nu}{\nu} \right) = \cfrac{E}{(1+\nu)}$$ (3.21)

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

$$\left\lbrace \begin{array}{rcl} \sigma_{11} & = & (\lambda + ... ... & = & 2 \mu \: \varepsilon_{12}\\ \end{array}\right. \text{.}$$ (3.22)

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

SOLUTION

$$%compliance matrix \left[ \begin{array}{c} \sigma_{11} \\ \... ...3} \\ 2 \varepsilon_{12} \end{array}\right] \: \: \blacksquare$$

Remember that there are only two independent constitutive parameters in linear isotropic elasticity. The usual pair choice is $E$ and $\nu$. However, there are other options depending on the application and equations used, e.g, $K$ and $G$. 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.

Figure 3.14: Table with equivalencies of the most used elastic moduli in linear isotropic elasticity.

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 $\uuline{\varepsilon} = \uuline{D} \: \uuline{S}$ is incorrect. Instead, the stress-strain relationship requires effective stress:

$$\uuline{\varepsilon} = \uuline{D} \: (\uuline{S} - P_p \uuline{I}) = \uuline{D} \: \uuline{\sigma}$$ (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.

Figure 3.15: The effective stress tensor.

Rigorously, the effective stress tensor needs a correction of pore pressure by the Biot coefficient $\alpha$ that accounts for solid grain deformation with changes in pore pressure.

$$\uuline{\varepsilon} = \uuline{D} \: (\uuline{S} - \alpha P_p \uuline{I})$$ (3.24)

For most problems, the assumption of $\alpha \sim 1$ is satisfactory. The rock matrix of tight sandstones and shales may have a Biot coeffiecient as low as $\alpha \sim 0.5$. 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 ( $S_v = S_{33}$) is a function of depth and rock bulk mass density.

$$S_v(z) = \int_0^z \rho_{bulk}(z) g \: dz$$ (3.25)

The effective vertical stress will be

$$\sigma_v(z) = S_v(z) - P_p$$ (3.26)

Let us now assume that the half space did not deform in horizontal directions ( $\varepsilon_{11}=\varepsilon_{22}=0$), 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.

Figure 3.16: Uniaxial strain condition with sedimentation in a tectonically passive environment. As sediment burial progresses the sediment compacts only in vertical direction.

Let us now use Equation 3.23 together with the equilibrium equation. Shear strains are zero. Hence $\uuline{\varepsilon} = [0,0,\varepsilon_{33},0,0,0]^T$. Then, the multiplication of $\uuline{\sigma} = \uuline{D} \: \uuline{\varepsilon}$, results in

$$\left\lbrace \begin{array}{l} \sigma_{11} = \sigma_{22} = \fr... ...(1-\nu) E}{(1+\nu)(1-2\nu)} \varepsilon_{33} \end{array}\right.$$ (3.27)

Let us express $\varepsilon_{33}$ as a function of $\sigma_{33}$, and plug it in the equation for horizontal stresses $\sigma_{11}$ and $\sigma_{22}$. The result is

$$\sigma_{11} = \sigma_{22} = \frac{\nu}{1-\nu} \sigma_{33}$$ (3.28)

or equivalently

$$\sigma_{h} = \frac{\nu}{1-\nu} \sigma_{v}$$ (3.29)

For typical values of $\nu \sim 0.25$, the horizontal stress coefficient is $\nu/(1-\nu) \sim 1/3$ (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 $\nu \rightarrow 0.5$ implies $\nu/(1-\nu) \rightarrow 1$. An “effective” $\nu \sim 0.5$ is applicable for fluids, soft rocks under undrained loading, and salt rocks.

Figure 3.17: Lateral stress coefficient as a function of Poisson's ratio. Solids do not push sideways with all its weight when compacted vertically.

The total horizontal stresses are obtained by adding pore pressure to the effective horizontal stresses: $S_{11} = \sigma_{11} + P_p$ and $S_{22} = \sigma_{22} + P_p$.

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 $S_h$ (assuming zero tectonic strains):

$$S_h = \sigma_{h} + P_p = \frac{\nu}{1-\nu} \sigma_{v} + P_p = \fr... ...nu}{1-\nu} (S_v - P_p) + P_p = \frac{\nu}{1-\nu} S_v + \frac{1-2\nu}{1-\nu} P_p$$ (3.30)

The gradient is the variation of pressure (or stress) with depth, i.e., derivative with respect to depth $z$. Assuming that the material properties are constant, then,

$$\frac{\Delta S_h}{\Delta z} = = \left( \frac{\nu}{1-\nu} \right) ... ..._v}{\Delta z} + \left( \frac{1-2\nu}{1-\nu} \right) \frac{\Delta P_p}{\Delta z}$$ (3.31)

For example, for onshore conditions with typical values $\Delta S_v / \Delta z \sim 1$ psi/ft, $\Delta P_p / \Delta z \sim 0.44$ psi/ft, and $\nu = 0.25$, the fracture gradient is $\Delta S_h / \Delta z = 0.63$ psi/ft. Figure XX shows a schematic example of the calculated fracture gradient.

Figure 3.18: Schematic representation of the prediction of minimum horizontal total stress $S_h$, a lower bound for the fracture gradient, using Eq. 3.31

Let us now relax the assumption of horizontal strains equal to zero, such that they are not zero $\varepsilon_{11} \neq \varepsilon_{22} \neq 0$, but are known quantities. We use the equation $\uuline{\sigma} = \uuline{C} \: \uuline{\varepsilon}$ again. Shear strains are zero. Hence $\uuline{\varepsilon} = [\varepsilon_{11},\varepsilon_{22},\varepsilon_{33},0,0,0]^T$. The resulting equations are,

$$\left\lbrace \begin{array}{l} \sigma_{11} = \cfrac{(1-\nu) E}... ...(1-\nu) E}{(1+\nu)(1-2\nu)} \varepsilon_{33} \end{array}\right.$$ (3.32)

Let us now substitute $\varepsilon_{33}$ in the equations of $\sigma_{11}$ and $\sigma_{22}$ as a function of $\sigma_{33}$. The result is:

$$\left\lbrace \begin{array}{l} \sigma_{11} = \cfrac{\nu}{1-\nu... ...1} + \cfrac{E}{1-\nu^2} \varepsilon_{22} \\ \end{array}\right.$$ (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 $\varepsilon _{Hmax}$ the maximum (compressive) tectonic strain, and $\varepsilon_{hmin}$ the minimum tectonic strain in a given direction. As a result the maximum effective horizontal stress and minimum horizontal stresses are:

$$\left\lbrace \begin{array}{l} \sigma_{Hmax} = \cfrac{\nu}{1-\... ...arepsilon_{Hmax} + E' \varepsilon_{hmin} \\ \end{array}\right.$$ (3.34)

where $E' = \frac{E}{1-\nu^2}$ 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 $S_v$
2. Calculate pore pressure $P_p$.
   1. If hydrostatic, use $\Delta P_p / \Delta z \sim 0.44$ psi/ft.
   2. If non-hydrostatic, use the method based on shale porosity $\phi_{shale}$ or use directly $P_p = \lambda_p S_v$ if $\lambda_p$ is given.
3. Calculate effective vertical stress: $\sigma_v = S_v - P_p$
4. Calculate effective horizontal stresses $\sigma_{hmin}$ and $\sigma_{Hmax}$:
   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:
   $$S_{Hmax} = \sigma_{Hmax} + P_p$$
   and
   $$S_{hmin} = \sigma_{hmin} + P_p.$$

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 $\mathrm{d}S_v/\mathrm{d}z = 23.8$ MPa/km, overpressure parameter $\lambda_p = 0.7$, shale Young's modulus $E = 5 \times 10^6$ psi, Poisson's ratio $\nu = 0.22$, and tectonic strains $\varepsilon_{hmin} = 0$ and $\varepsilon_{Hmax} = 0.0002$.

SOLUTION
At a depth of 7950 ft

$$S_v = 23.8 \frac{\text{MPa}}{\text{km}} \times \frac{1 \frac{\tex... ...}}}{23 \frac{\text{MPa}}{\text{km}}} \times 7950 \text{ ft} = 8227 \text {psi}$$

$$P_p = \lambda_p S_v = 0.7 \times 8227$$$$\text {psi} = 5759 \text{ psi}$$

Hence, the effective vertical stress is

$$\sigma_v = S_v - P_p = 8227$$$$\text {psi} - 5759 \text{ psi} = 2468 \text{ psi}$$

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

$$E' = \frac{E}{1-\nu^2} = \frac{5 \times 10^6 \text{ psi}}{1-0.22^2} = 5.25 \times 10^6 \text{ psi}$$

Then,

$$\left\lbrace \begin{array}{l} \sigma_{Hmax} = \frac{\nu}{1-\n... ...xt{ psi} \times 0.0002 = 927 \text{ psi} \\ \end{array}\right.$$

Finally, let us compute the total horizontal stresses by adding pore pressure:

$$\left\lbrace \begin{array}{l} S_{Hmax} = \sigma_{Hmax} + P_p ... ...{ psi} = 6686 \text{ psi} \end{array}\right. \: \: \blacksquare$$

Horizontal stresses tend to be different in most basins, hence, $S_{Hmax} > S_{hmin}$. In practice, differences $S_{Hmax} - S_{hmin} = \sigma_{Hmax} - \sigma_{hmin} > 0.7$ MPa (100 psi) tend to be enough to force hydraulic fractures to grow mostly perpendicular to $S_{hmin}$ 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 $S_{Hmax} \sim S_{hmin}$.

3.3.5 Calculation of reservoir compressibility with linear elasticity

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

$$\frac{\partial P_p}{\partial t}= \frac{k}{\mu C_t} \frac{\partial^2 P_p}{\partial x^2}$$ (3.35)

Where the total compressibility is $C_t = C_g S_g + C_w S_w + C_o S_o + C_{pp}$. Reservoir simulators usually calculate the fluid compressibility $(C_g S_g + C_w S_w + C_o S_o)$ based in phase behavior, hence, the only required input is $C_{pp}$. For example, compaction drive is proportional to rock compressibility (see https://petrowiki.org/Compaction_drive_reservoirs).

The pore volume compressibility $C_{pp}$ tells us what the change of pore volume $V_p$ is due to a change in pore pressure:

$$C_{pp} = \left. \frac{1}{V_p} \frac{\mathrm{d}V_p}{\mathrm{d}P_p} \right\vert _{S_v,\varepsilon_h}$$ (3.36)

The equation above captures reservoir boundary conditions in which the total vertical stress $S_v$ remains constant (overburden above the reservoir does not change) and there is no change of lateral strain $\varepsilon_h$, 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 $C_{pp}$ are derived from bulk volume measurements. Let us assume that the change of pore volume $\mathrm{d}V_p$ is equal to the change of bulk volume $\mathrm{d}V_b$, which means that all bulk deformation is caused by change of porosity. Hence it is possible to rewrite the definition of $C_{pp}$ as

$$C_{pp} = \frac{1}{\frac{V_p}{V_b}} \left( \frac{1}{V_b} \left. \frac{\mathrm{d}V_b}{\mathrm{d}P_p} \right\vert _{S_v,\varepsilon_h} \right)$$ (3.37)

Porosity is defined as $\phi = V_p/V_b$ and the term between parenthesis is defined as the bulk compressibility under uniaxial condition $C_{bp}$ (notice $\varepsilon_{vol} = \mathrm{d}V_b/V_b$). Hence, the parameter $C_{pp}$ is linked to the bulk rock compressibility $C_{bp}$ through porosity:

$$C_{pp} = \frac{C_{bp}}{\phi}$$ (3.38)

where the bulk compressibility with no lateral strain is approximately equal to the inverse of the bulk constrained modulus $C_{bp} \sim M^{-1}$, where $M = (1-\nu) E / [(1+\nu)(1-2\nu)]$ 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 $E$ and $\nu$ as,

$$C_{pp} = \frac{(1+\nu)(1-2\nu)}{(1-\nu) E \phi}$$ (3.39)

Figure 3.19: Pressure depletion causes reservoir compaction. The reservoir compressibility is inversely proportional to the constrained modulus $M$.

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 $\times 10^{-6}$    psi$^{-1}$, where $\times 10^{-6}$    psi$^{-1} = \mu$   sip. Stiff well cemented rocks have low pore volume compressibility $\sim 2 \: \mu$   sip while uncemented loose sediments tend to have high pore volume compressibility $\sim 30 \: \mu$   sip.

EXAMPLE 3.6: Calculate the pore compressibility of a rock tested in the laboratory with porosity $\phi = 0.20$, Young's modulus $E = 10$ GPa, and $\nu = 0.20$. Provide the solution in [$10^{6}$psi]$^{-1}$ units.

SOLUTION
The constrained modulus is

$$M = \frac{(1-\nu) E}{(1+\nu)(1-2\nu)} = \frac{(1-0.20) 10 \text{ ... ...(1+0.20)(1-2 \times 0.20)} = 11.11 \text{ GPa} = 1.6 \times 10^{6} \text{ psi}$$

Hence,

$$C_{pp} = \frac{1}{M \phi} = \frac{1}{1.6 \times 10^{6} \text{psi}... ...} = 3.1 \: [10^{6} \text{psi}]^{-1} = 3.1 \: \mu \text{sip} \: \: \blacksquare$$

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 $\uline{u}$ as the unknown:

$$(\lambda + G) \nabla (\nabla \cdot \uline{u}) + G \nabla^2 \uline{u} + \rho \uline{b} = 0$$ (3.40)

where $\lambda = (\nu E)/[(1+\nu)(1-2\nu)]$ is the first Lamé parameter, $G$ is the shear modulus, $\rho$ is the rock bulk mass density, and $\uline{b}$ is the body force acceleration vector (usually gravity). A review of the gradient $\nabla ()$, divergence $\nabla \cdot ()$ and Laplacian $\nabla^2 ()$ 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.

Figure 3.20: General continuum mechanics problem.