2. Continuum mechanics

This document is a draft. Find hand written notes here: https://github.com/dnicolasespinoza/GeomechanicsJupyter/tree/master/ClassNotes.

2.1 Equations needed to solve a continuum mechanics problem

Figure 2.1: General solution for a continuum mechanics problem.

2.2 Kinematic equations

2.2.1 One-dimensional case

Figure 2.2: Kinematic equations 1D and 3D.

2.2.2 Three-dimensional case

Figure 2.3: Displacement field jacobian matrix.

Figure 2.4: Strain tensor.

Figure 2.5: Volumetric strain.

2.3 Constitutive equations

2.3.1 Linear elasticity

Figure 2.6: Linear elasticity.

Figure 2.7: Voigt notation.

2.3.2 Stiffness matrix reduction to orthotropic elasticity

Figure 2.8: Orthotropic elasticity.

2.3.3 Stiffness matrix for vertical transverse isotropic elasticity

Figure 2.9: Vertical transverse isotropic elasticity.

2.3.4 Stiffness matrix for isotropic elasticity

Figure 2.10: Isotropic elasticity.

Figure 2.11: Compliance matrix for isotropic elasticity.

Figure 2.12: Stiffness matrix for isotropic elasticity.

2.4 1D subsurface mechanical models

2.4.1 Uniaxial strain case

Figure 2.13: Stress solution for uniaxial strain condition in isotropic elastic media.

2.4.2 Triaxial strain case

Figure 2.14: Stress solution for tri-axial strain condition in isotropic elastic media.

2.4.3 Parameters for Mechanical Earth Models

Figure 2.15: Elastic wave velocities and dynamic elastic coefficients.

Figure 2.16: Degradation of elastic modulus (G) with strain magnitude and strain rate.

Figure 2.17: Correlations between static and dynamic parameters.

Figure 2.18: MEM application for hydraulic fracture geometry.

2.5 Parameters for VTI MEM

VTI compliance matrix (3 - vertical direction perpendicular to bedding):

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

where and $G_h = \cfrac{E_h}{2(1+\nu_h)}$ and $G_v$ is not related to the other parameters.

In terms of stiffness coefficients:

$$E_h = \frac{(C_{11}-C_{12}) \left[ C_{33}(C_{11}+C_{12})-2\: C_{13}^2 \right]}{C_{11}C_{33}–C_{13}^2}$$

$$E_v = C_{33} - \frac{2\:C_{13}^2}{C_{11}+C_{12}}$$

$$\nu_h = \frac{C_{12}C_{33}-C_{13}^2}{C_{11}C_{33}-C_{13}^2}$$

$$\nu_v = \frac{C_{13}}{C_{11}+C_{12}}$$

$$G_v = C_{44}$$

$$G_h = C_{66} = \frac{C_{11}-C_{12}}{2}$$

VTI stiffness matrix (3 - vertical direction perpendicular to bedding):

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

or in terms of Young moduli and Poisson ratios

$$C_{11} = \left[ \frac{1}{(1-\nu_h) E_v - 2 \nu_v^2 E_h} \right] \left( \frac{E_h E_v - \nu_v^2 E_h^2}{1+\nu_h} \right)$$

$$C_{33} = \left[ \frac{1}{(1-\nu_h) E_v - 2 \nu_v^2 E_h} \right] (E_v^2 - \nu_h E_v^2)$$

$$C_{12} = \left[ \frac{1}{(1-\nu_h) E_v - 2 \nu_v^2 E_h} \right] \left( \frac{\nu_v^2 E_h^2 + \nu_h E_h E_v}{1+\nu_h} \right)$$

$$C_{13} = \left[ \frac{1}{(1-\nu_h) E_v - 2 \nu_v^2 E_h} \right] (\nu_v E_h E_v)$$

$$C_{66} = \frac{C_{11}-C_{12}}{2} = G_h = \cfrac{E_h}{2(1+\nu_h)}$$

The parameter $C_{44}$ is independent of all other parameters.

2.5.1 Laboratory measurements

Figure 2.19: Measurement of VTI static properties.

Figure 2.20: Measurement of VTI dynamic properties.

Figure 2.21: Correlations between dynamic and static VTI stiffness coefficients.

2.5.2 Anisotropy evaluation and models

Figure 2.22: Dimensionless parameters for anisotropic media.

Figure 2.23: Parallel and series models for lamianted media.

Figure 2.24: Bulk volumetric modulus of a VTI material.

2.5.2.1 VTI 1D MEM

TBD

2.6 WP3: Horizontal Stresses Computed with Linear Elasticity

2.6.1 Exercise 1: One-dimensional mechanical Earth model (MEM)

Download the file LostHills.xls. We would like to know the state of stress in the subsurface and its influence on a hydraulic fracture completion. At every depth (and data-point) along the vertical well:

1. Compute (and plot) total vertical stress as a function of depth (you may assume homogeneous rock above 1,750 ft), and overpressure parameter.
2. Compute dynamic Poisson’s ratio and dynamic Young’s modulus from compressive and shear slowness (be careful with unit conversion).
3. Compute static Young's modulus using a coefficient $E_{static} = 0.65 \cdot E_{dynamic}$.
4. Compute (and plot) static plane strain modulus $E'_{static} = E_{static} / (1-\nu^2)$ (Poisson ratio remains the same).
5. Compute (and plot) horizontal stress assuming theory of elasticity and no tectonic strains.
6. Compute (and plot) total maximum and minimum horizontal stress assuming theory of elasticity and $\varepsilon_{Hmax}=0.0015$ and $\varepsilon_{hmin}=0$.
7. The pay-zone is between 2,100 ft and 2,450 ft. A hydraulic fracture is planned to be executed with a vertical well at a depth between 2,130 ft and 2,160 ft. What will be the height of this fracture? Will it reach out to the entire pay zone?

Figure 2.25: Bulk rock mass density, P-wave slowness, and S-wave slowness along a vertical well in the Lost Hills field.

2.7 Navier's equation

Figure 2.26: General solution for a continuum mechanics problem.

2.7.1 Derivation

Figure 2.27: Combination of equilibrium, constitutive and kinematic equations to obtain Navier's equations.

2.7.2 Analytical and numerical solutions

Figure 2.28: Solution strategies to Navier's equation.

2.7.3 Weak formulation

Figure 2.29: Solution of static problems through virtual work.

Figure 2.30: Variational form of equilibrium equations.

Figure 2.31: Strain energy and inclusion of linear elasticity.

2.7.4 Using FreeFem++ to solve equation of linear elasticity - by Igor Shovkun

2.7.4.1 Introduction

Equilibrium equation:

$$-\nabla \cdot \sigma = f \ \Omega$$ (2.3)

Generic constitutive model:

$$\sigma_{ij} = C_{ijkl}\varepsilon_{ij}$$ (2.4)

Linear elasticity:

$$\sigma_{ij} = \lambda \delta_{ij}\nabla u + 2G\varepsilon_{ij}$$ (2.5)

Definition of strain:

$$\varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)$$ (2.6)

The corresponding variationnal form is:

$$\int_{\Omega}\sigma(u):\varepsilon(v) dV + \int_{\partial \Omega_n}\tau \cdot n \, dS = \int_{\Omega}{f \cdot v \ dV}$$ (2.7)

where : denotes the tensor scalar product, i.e. $a : b = \sum_{i, j} a_{ij}b_{ij}$, $\partial \Omega_n$ is a Neumann boundary, $\tau$ is the boundary force vector, and $n$ is the outward normal vector (perpendicular to the infinitesimal surface $\vec{dS}$).

In the case of plane strain linear elasticity, the variational form look like this:

$$\int_{\Omega} \left[ 2G\varepsilon_{ij}(u)\varepsilon_{ij}(v) + \... ...t] dV + \int_{\partial \Omega_n}\tau \cdot n \, dS = \int_{\Omega}f\cdot v \ dV$$ (2.8)

2.7.4.2 Wellbore problem in FreeFem++

One problem of interest is to find stress component around a wellbore. Mathematically, this problem implicates solving Eq. 2.3 in an infinite domain with a pressurized circular opening with radius $R$ and can be written as follows:

$$\begin{cases} -\nabla \sigma = f \\ \sigma_{xx}(x=\pm\infty, y) = S_1 \\ \sigma_{yy}(x, y=\pm\infty) = S_2 \\ \end{cases}$$ (2.9)

The analytical solution for this problem reads:

$$\begin{aligned} \sigma_{rr} = & \frac{1}{2}(S_1 + S_2 - 2P_w) \... ...eft( \frac{R}{r} \right)^4 \right]\sin(2\theta) \\ \end{aligned}$$

Define problem parameters:

  % Start your code-block

// Dimensions
real ySize = 10.;  // y-size of the domain
real xSize = 10.;  // x-size of the domain
real R = 0.1; // wellbore radius
//Elastic constants
real E = 1e10 ; // young's modulus
real nu = 0.3;  // poisson's ratio

Then we calculate Lame constant and shear modulus that are used in the formulation:

  % Start your code-block

real G = E/(2*(1+nu)); // shear modulus
real lambda = E*nu/((1+nu)*(1-2*nu)); // Lame constant

Finally, we specify values for stress boundary conditions:

  % Start your code-block

// Stresses
real Sx = 10e6;
real Sy = 10e6;
real Pwell = 0.0;

Next we proceed to the definition of the domain. We first define all the domain boundaries and then create mesh:

  % Start your code-block

// First define boundaries
border Right(t=-ySize/2,ySize/2) {x = xSize/2; y = t;}
border Top(t=xSize/2,-xSize/2) {x = t; y = ySize/2;}
border Left(t=ySize/2,-ySize/2) {x = -xSize/2; y = t;}
border Bottom(t=-xSize/2,xSize/2) {x = t; y = -ySize/2;}
border Well(t = 0, -2*pi) {x = R*cos(t); y = R*sin(t);}
// Create mesh
int n = 20; // number of mesh nodes on the outer borders
int nwell = 50; // number of mesh nodes on wellobre
mesh Omega = buildmesh(Right(n)+Top(n)+Left(n)+Bottom(n)+Well(nwell));

The next step is to define Finite Element spaces. These object designate the type of shape functions used to solve the problem. In this example we will be using “P1” elements for displacement: continuous Galerkin polynomials of the first degree defined on triangles. To approximate stresses we will use “P0” elements - piecewise constants.

  % Start your code-block

// FE spaces
fespace Displacement(Omega, P1); // linear shape functions
fespace Stress(Omega, P0); // piecewise constants

Next, we will tell the interpreter to allocate several objects to store displacement components (objects of type Displacement) and stresse components (objects of type Stress). $u1$ and $u2$ denote x and y components of the displacement vector. $v1$ and $v2$ are the components of the virtual displacement vector; we do not need those vectors per se, but we use these objects to code up the system of equations. We also define objects $sigmaxx$, $sigmaxy$, and $sigmayy$ to store stress components. Note that we only need 3 components of the stress tensor because it is symmetric.

  % Start your code-block

Displacement u1, u2, v1, v2;
Stress sigmaxx, sigmayy, sigmaxy;

Next we define macros. Macros are special compiler instructions that are inserted into the part of code where they are called during compilation (as opposed to runtime). In this code we use macros in order to make the system of equation look neat later on. The following code defines two macros for strain and stress “vectors”.

  % Start your code-block

// definition of 2 macros:
// macro for strain
macro e(u1,u2)
  [
    dx(u1),
    (dy(u1)+dx(u2))/2,
    (dx(u2)+dy(u1))/2,
    dy(u2)
  ]
  // eps_xx, eps_xy, eps_yx, eps_yy

// macro for stress
macro sigma(u1,u2)
  [
    (lambda+2.*G)*e(u1,u2)[0] + lambda*e(u1,u2)[3],
    2.*G*e(u1,u2)[1],
    2.*G*e(u1,u2)[2],
    lambda*e(u1,u2)[0] + (lambda+2.*G)*e(u1,u2)[3]
  ] // stress sxx,sxy,syx,syy

Now, we can finally define the system to be solved corresponding to the Formulation (ref).

// Define system of equations
problem Elasticity([u1,u2], [v1,v2]) =
    int2d(Omega) (  sigma(u1,u2)'*e(v1,v2)  )
  // Boundary conditions
  + int1d(Omega, Right) (Sx*v1)
  - int1d(Omega, Left) (Sx*v1)
  + int1d(Omega, Top) (Sy*v2)
  - int1d(Omega, Bottom) (Sy*v2)
  + int1d(Omega, Well) (Pwell*(N.x*v1 + N.y*v2))
  ;

In order to solve the problem, we simply need to call this subroutine:

  % Start your code-block

// Solve system
Elasticity;

Finally, in order to calculate stresses, we use previously define macros again:

// Stresses
sigmaxx = sigma(u1, u2)[0];
sigmayy = sigma(u1, u2)[3];
sigmaxy = sigma(u1, u2)[1]; // we could use [2] as well

2.7.4.3 Stresses around a fracture

The changes in the code are minor. Instead of a wellbore we need to define a fracture:

  % Start your code-block

real xf = 40; // fracture half-length
real fw = 0.1; // fracture half-width
border Frac(t = 0, -2*pi) {x = fw*cos(t); y = xf*sin(t);}

The bondary condition on the fracture boundary are also a little bit different:

// Define system of equations
problem Elasticity([u1,u2], [v1,v2]) =
    int2d(Omega) (  sigma(u1,u2)'*e(v1,v2)  )
  // Boundary conditions
  + int1d(Omega, Right) (Sx*v1)
  - int1d(Omega, Left) (Sx*v1)
  + int1d(Omega, Top) (Sy*v2)
  - int1d(Omega, Bottom) (Sy*v2)
  // condition only on one component
  + int1d(Omega, Frac) (Pfrac*(N.x*v1))
  ;

2.7.5 Exercise 2: Vertical transverse isotropic elastic properties

[To be developed]

2.7.6 Exercise 3: Displacement Field and Strain Field

[To be developed]

2.8 WP4: Solution of Navier's Equation for the Stress Field

2.8.1 Exercise 1: Stresses around a wellbore

Consider a 2D problem of a circular cavity subjected to far field effective stresses $\sigma_{xx}$ = 12 MPa and $\sigma_{yy}$ = 3 MPa. The diameter of the cavity is 0.2 m. Rock properties: $E$ = 10 GPa, $\nu$ = 0.20, unconfined compression strength $UCS$ = 30 MPa, tensile strength $T_s$ = 2 MPa.

1. Using Kirsch equations compute (and plot) $\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$ for a domain $x$ = [-1m, +1m], and $y$ = [-1m, +1m]. You may define a polar grid for $(r,\theta)$. How far does the presence of the wellbore influence stresses?
2. Using Kirsch equations compute (and plot) stresses in a line ($x$ = [0.1m, 1m], $y$ = 0 m) and ($x$ = 0 m, $y$ = [0.1 m, 1 m]). Equations in Ch. 6.2 (https://dnicolasespinoza.github.io/)
3. Using Kirsch equations compute (and plot) $\sigma_{rr}$ and $\sigma_{\theta\theta}$ for $r$ = 0.1 m. Is there any section of the rock in shear or tensile failure? Where?
4. Use FreeFEM++ (http://www3.freefem.org/) or FEniCS (https://fenicsproject.org/) to solve the same problem ( $\sigma_{xx}$, $\sigma_{yy}$ and $\sigma_{xy}$) assuming a domain size 2 m by 2 m. Compute $\sigma_{xx}$ and $\sigma_{yy}$ for the same lines as in point (b), and compare with Kirsch's analytical solution. Repeat the process for a domain size 0.5 m by 0.5 m. Are there any differences? Why?
5. Plot the displacement field.
6. EXTRA: compute principal stresses within FreeFEM++ and plot $\sigma_{rr}$ and $\sigma_{\theta\theta}$.

Hint: An example code for 2D elasticity in FreeFEM++ and the corresponding explanation are available at https://github.com/dnicolasespinoza/GeomechanicsJupyter/: Kirsch_Shovkun.edp and FreeFEM_Tutorial_Shovkun.pdf. You can also try FreeFEM++ online here: https://freefem.org/tryit.

2.8.2 Exercise 2: Stresses around a planar fracture

Consider a 2D problem of an elliptical fracture (half-length $c$ = 10 m). Solve the problem using just half of the domain. Set the fracture along the left boundary of a domain: $x$ = [0 m, 100 m] and $y$ = [-50 m, 50 m], with fracture center at $(x,y)=$ (0,0) m. This boundary will have a pressure boundary condition. All other boundaries will have zero displacement. Rock properties: $E$ = 30 GPa, $\nu$ = 0.20.

Figure 2.32: Planar fracture model.

1. Use FreeFEM++ (http://www3.freefem.org/) or FEniCS (https://fenicsproject.org/) to solve for $\sigma_{xx}$, $\sigma_{yy}$ and $\sigma_{xy}$ imposing a fracture pressure $p$ = 10 MPa. Plot results.
2. Export and plot stress perpendicular to the fracture direction $\sigma_{xx}$ at the middle of the fracture (L1 = ($x$ = [0, 100 m], $y$ = 0 m), Figure 2.32). How far does the influence of the fracture extend?
3. Plot $x$-displacements at the face of the fracture. Compare with analytical equation. Equations in Ch. 7.3.2 (https://dnicolasespinoza.github.io/).
4. Plot $\sigma_{xx}$ along fracture length and beyond fracture tips (line L2 = ($x$ = 0 m, $y=$ [-50, 50]) m, Figure 2.32) and compare with analytical Griffith solution.
5. EXTRA: compare FreeFEM++ solution to analytical solution by Sneddon and Elliot, 1946. (Eq. 20 in https://doi.org/10.1016/j.fuel.2017.01.057).

Figure 2.33: Fracture width and stress beyond tip.

Figure 2.34: Stress intensity factor.