Single Fracture Under Shear Compression

Context

In this example, a single fracture is simulated using a Lagrange contact model in a 2D infinite domain and subjected to a constant uniaxial compressive remote stress (Franceschini et al., 2020). An analytical solution (Phan et al., 2003) is available for verifying the accuracy of the numerical results, providing an analytical form for the normal traction and slip on the fracture surface due to frictional contact. In this example, the TimeHistory function and a Python script are used to output and postprocess multi-dimensional data (traction and displacement on the fracture surface).

Input file

Everything required is contained within these GEOS input files and one mesh file located at:

inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml

inputFiles/lagrangianContactMechanics/SingleFracCompression_benchmark.xml

inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_benchmark.xml

inputFiles/lagrangianContactMechanics/crackInPlane_benchmark.vtu

Description of the case

We simulate an inclined fracture under a compressive horizontal stress ( $\sigma$ ), as shown below. This fracture is placed in an infinite, homogeneous, isotropic, and elastic medium. Uniaxial compression and frictional contact on the fracture surface cause mechanical deformation to the surrounding rock and sliding along the fracture plane. For verification purposes, plane strain deformation and Coulomb failure criterion are considered in this numerical model.

../../../../../../_images/sketch2.png — Fig. 16 Sketch of the problem

To simulate this phenomenon, we use a Lagrange contact model. Displacement and stress fields on the fracture plane are calculated numerically. Predictions of the normal traction ( $t_N$ ) and slip ( $g_T$ ) on the fracture surface are compared with the corresponding analytical solution (Phan et al., 2003).

$t_N = - \sigma {( \text{sin} \left( {\psi} \right) )}^{ 2 }$

$g_T = \frac{ 4 ( 1- {\nu}^{ 2 }) }{ E } {( \sigma \text{sin} \left( {\psi} \right) { ( \text{cos} \left( {\psi} \right) - \text{sin} \left( {\psi} \right) \text{tan} \left( {\theta} \right) )} )} \, \sqrt{ { b }^{ 2 } - { (b - \xi) }^{ 2 } }$

where $\psi$ is the inclination angle, $\nu$ is Poisson’s ratio, $E$ is Young’s modulus, $\theta$ is the friction angle, $b$ is the fracture half-length, $\xi$ is a local coordinate on the fracture varying in the range [ $0$ , $2b$ ].

In this example, we focus our attention on the Mesh tags, the Constitutive tags, and the FieldSpecifications tags.

Mesh

The following figure shows the mesh used in this problem.

Here, we load the mesh with VTKMesh (see Importing the Mesh). The syntax to import external meshes is simple: in the XML file, the mesh file crackInPlane_benchmark.vtu is included with its relative or absolute path to the location of the GEOS XML file and a user-specified label (here CubeHex) is given to the mesh object. This unstructured mesh contains quadrilaterals elements and interface elements. Refinement is performed to conform with the fracture geometry specified in the Geometry section.

  <Mesh>
    <VTKMesh
      name="CubeHex"
      file="crackInPlane_benchmark.vtu"/>
  </Mesh>

  <Geometry>
    <Rectangle
      name="fracture"
      normal="{-0.342020143325669, 0.939692620785908, 0.0}"
      origin="{0.0, 0.0, 0.0}"
      lengthVector="{0.939692620785908, 0.342020143325669, 0.0}"
      widthVector="{0.0, 0.0, 1.0}"
      dimensions="{ 2, 10 }"/>

    <Rectangle
      name="core"
      normal="{-0.342020143325669, 0.939692620785908, 0.0}"
      origin="{0.0, 0.0, 0.0}"
      lengthVector="{0.939692620785908, 0.342020143325669, 0.0}"
      widthVector="{0.0, 0.0, 1.0}"
      dimensions="{ 2, 10 }"/>

    <Box
      name="rightPoint"
      xMin="{ 39.9, -40.1, -0.001}"
      xMax="{ 40.1,  40.1,  0.051}"/>

    <Box
      name="leftPoint"
      xMin="{-40.1, -40.1, -0.001}"
      xMax="{-39.9,  40.1,  0.051}"/>

    <Box
      name="topPoint"
      xMin="{-40.1, 39.9, -0.001}"
      xMax="{ 40.1, 40.1,  0.051}"/>

    <Box
      name="bottomPoint"
      xMin="{-40.1, -40.1, -0.001}"
      xMax="{ 40.1, -39.9,  0.051}"/>

    <Box
      name="front"
      xMin="{-40.1, -40.1, -0.001}"
      xMax="{ 40.1,  40.1,  0.001}"/>

    <Box
      name="rear"
      xMin="{-40.1, -40.1, 0.049}"
      xMax="{ 40.1,  40.1, 0.051}"/>

    <Box
      name="xmin"
      xMin="{-40.1, -40.1, -0.001}"
      xMax="{-39.9,  40.1,  0.051}"/>

    <Box
      name="xmax"
      xMin="{39.9, -40.1, -0.001}"
      xMax="{40.1,  40.1, 0.051}"/>
  </Geometry>

Solid mechanics solver

GEOS is a multi-physics platform. Different combinations of physics solvers available in the code can be applied in different regions of the domain and be functional at different stages of the simulation. The Solvers tag in the XML file is used to list and parameterize these solvers.

To specify a coupling between two different solvers, we define and characterize each single-physics solver separately. Then, we customize a coupling solver between these single-physics solvers as an additional solver. This approach allows for generality and flexibility in constructing multi-physics solvers. Each single-physics solver should be given a meaningful and distinct name because GEOS recognizes these single-physics solvers based on their given names to create the coupling.

To setup a coupling between rock and fracture deformations, we define three different solvers:

For solving the frictional contact, we define a Lagrangian contact solver, called here lagrangiancontact. In this solver, we specify targetRegions that includes both the continuum region Region and the discontinuum region Fracture where the solver is applied to couple rock and fracture deformation. The contact constitutive law used for the fracture elements is named fractureMaterial, and defined later in the Constitutive section.
Rock deformations are handled by a solid mechanics solver SolidMechanics_LagrangianFEM. This solid mechanics solver (see Solid Mechanics Solver) is based on the Lagrangian finite element formulation. The problem is run as QuasiStatic without considering inertial effects. The computational domain is discretized by FE1, which is defined in the NumericalMethods section. The solid material is named rock, and its mechanical properties are specified later in the Constitutive section.
The solver SurfaceGenerator defines the fracture region and rock toughness.

  <Solvers
    gravityVector="{0.0, 0.0, 0.0}">
    <SolidMechanicsLagrangeContact
      name="lagrangiancontact"
      timeIntegrationOption="QuasiStatic"
      stabilizationName="TPFAstabilization"
      logLevel="1"
      discretization="FE1"
      targetRegions="{ Region, Fracture }">
      <NonlinearSolverParameters
        newtonTol="1.0e-8"
        logLevel="2"
        newtonMaxIter="10"
        maxNumConfigurationAttempts="10"
        lineSearchAction="Require"
        lineSearchMaxCuts="2"
        maxTimeStepCuts="2"/>
      <LinearSolverParameters
        solverType="direct"
        directParallel="0"
        logLevel="0"/>
    </SolidMechanicsLagrangeContact>
  </Solvers>

Constitutive laws

For this specific problem, we simulate the elastic deformation and fracture slippage caused by uniaxial compression. A homogeneous and isotropic domain with one solid material is assumed, with mechanical properties specified in the Constitutive section.

Fracture surface slippage is assumed to be governed by the Coulomb failure criterion. The contact constitutive behavior is named fractureMaterial in the Coulomb block, where cohesion cohesion="0.0" and friction coefficient frictionCoefficient="0.577350269" are specified.

  <Constitutive>
    <ElasticIsotropic
      name="rock"
      defaultDensity="2700"
      defaultBulkModulus="16.66666666666666e9"
      defaultShearModulus="1.0e10"/>

    <Coulomb
      name="frictionLaw"
      cohesion="0.0"
      frictionCoefficient="0.577350269"/>
  </Constitutive>

Recall that in the SolidMechanics_LagrangianFEM section, rock is the material of the computational domain. Here, the isotropic elastic model ElasticIsotropic is used to simulate the mechanical behavior of rock.

All constitutive parameters such as density, bulk modulus, and shear modulus are specified in the International System of Units.

Time history function

In the Tasks section, PackCollection tasks are defined to collect time history information from fields. Either the entire field or specified named sets of indices in the field can be collected. In this example, tractionCollection and displacementJumpCollection tasks are specified to output the local traction fieldName="traction" and relative displacement fieldName="displacementJump" on the fracture surface.

  <Tasks>
    <PackCollection
      name="tractionCollection"
      objectPath="ElementRegions/Fracture/faceElementSubRegion"
      fieldName="traction"/>

    <PackCollection
      name="displacementJumpCollection"
      objectPath="ElementRegions/Fracture/faceElementSubRegion"
      fieldName="displacementJump"/>   
    </Tasks>

These two tasks are triggered using the Event management, with PeriodicEvent defined for these recurring tasks. GEOS writes two files named after the string defined in the filename keyword and formatted as HDF5 files (displacementJump_history.hdf5 and traction_history.hdf5). The TimeHistory file contains the collected time history information from each specified time history collector. This information includes datasets for the simulation time, element center defined in the local coordinate system, and the time history information. Then, a Python script is used to access and plot any specified subset of the time history data for verification and visualization.

Initial and boundary conditions

The next step is to specify fields, including:

The initial value (the remote compressive stress needs to be initialized),

The boundary conditions (the constraints of the outer boundaries have to be set).

In this tutorial, we specify an uniaxial horizontal stress ( $\sigma_x$ = -1.0e8 Pa). The remaining parts of the outer boundaries are subjected to roller constraints. These boundary conditions are set up through the FieldSpecifications section.

  <FieldSpecifications>
    <FieldSpecification
      name="frac"
      initialCondition="1"
      setNames="{ fracture }"
      objectPath="faceManager"
      fieldName="ruptureState"
      scale="1"/>

    <FieldSpecification
      name="separableFace"
      initialCondition="1"
      setNames="{ core }"
      objectPath="faceManager"
      fieldName="isFaceSeparable"
      scale="1"/>

    <FieldSpecification
      name="xconstraint"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="0"
      scale="0.0"
      setNames="{ leftPoint, rightPoint }"/>

    <FieldSpecification
      name="yconstraint"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="1"
      scale="0.0"
      setNames="{ bottomPoint, topPoint }"/>

    <FieldSpecification
      name="zconstraint"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="2"
      scale="0.0"
      setNames="{ front, rear }"/>

    <FieldSpecification 
      name="Sigmax"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions/Region"
      fieldName="rock_stress"
      component="0"
      scale="-1.0e8"/>
  </FieldSpecifications>

Note that the remote stress has a negative value, due to the negative sign convention for compressive stresses in GEOS.

The parameters used in the simulation are summarized in the following table.

Symbol	Parameter	Unit	Value
$K$	Bulk Modulus	[GPa]	16.67
$G$	Shear Modulus	[GPa]	10.0
$\sigma$	Compressive Stress	[MPa]	-100.0
$\theta$	Friction Angle	[Degree]	30.0
$\psi$	Inclination Angle	[Degree]	20.0
$b$	Fracture Half Length	[m]	1.0

Inspecting results

We request VTK-format output files and use Paraview to visualize the results. The following figure shows the distribution of $u_{y}$ in the computational domain.

../../../../../../_images/displacement_yy.png — Fig. 18 Simulation result of $u_{y}$ 

The next figure shows the distribution of relative shear displacement values on the fracture surface.

../../../../../../_images/slip1.png — Fig. 19 Simulation result of fracture slip

The figure below shows a comparison between the numerical predictions (marks) and the corresponding analytical solutions (solid curves) for the normal traction ( $t_N$ ) and slip ( $g_T$ ) distributions on the fracture surface. One can observe that the numerical results obtained by GEOS and the analytical solutions are nearly identical.

(Source code)

../../../../../../_images/singleFracCompressionFigure.png

To go further

Feedback on this example

For any feedback on this example, please submit a GitHub issue on the project’s GitHub page.