# Static fracture opening under a constant pressure – Sneddon solution

This page is based on a Jupyter notebook.

## Problem description

### Static fracture opening under a constant pressure

Consider a line fracture $[-a_0, a_0] \times \{0\}$ ($a_0$ = 0.1) with no external loading and an internal fluid pressure of $p=1$ was applied to the fracture surfaces. The results are compared to the analytical solution (Sneddon et al., 1969) for the fracture half-opening: \begin{eqnarray} u(x,0) = \frac{2 p a_0}{ E’} \sqrt{1-(x/a_0)^2} , \label{eq:sneddon_2D_static} \end{eqnarray}

where $u$ is the displacement $E'$ is the plane strain Young’s modulus ($E' = E/(1-\nu^2)$) with $\nu$ is Poisson’s ratio, $p$ is the fluid pressure inside the fracture. To account for the infinite boundaries in the closed-form solution, we considered a large finite domain $[-10a_o,10a_o] \times [-10a_o,10a_o]$. The effective element size, $h$, is $1\times10^{-2}$.

• In order to have a static pressurized fracture, add, $\texttt{static}$ in the project file. It’s important to note that static pressurized fracture implementation assumes $p=1$. Yoshioka et al., 2019 discussed using real material properties and rescaling the phase-field energy functional.

# Input data

The material and geometrical properties are listed in the Table below. It’s worth noting that the properties are dimensionless and scaled.

Name Value Symbol
Young’s modulus 1 $E$
Poisson’s ratio 0.15 $\nu$
Fracture toughness 1 $G_{c}$
Regularization parameter 2$h$ $\ell_s$
Pressure 1 $p$
Length 2 $L$
Height 2 $H$
Initial crack length 0.2 $2a_0$
import math
import os
import time

import gmsh
import matplotlib.pyplot as plt
import numpy as np
from ogs6py import ogs

pi = math.pi
plt.rcParams["text.usetex"] = True
E = 1.0
nu = 0.15
Gc = 1.0
P = 1.0
h = 0.01

Orientation = 0
a0 = 0.1  # half of the initial crack length
n_slices = (
2 * (round(3.0 * a0 / h)) + 1
)  # number of slices for calcute width of fracture

phasefield_model = "AT1"
h_list = [0.01]  # list of mesh sizes (h)
# h_list =[0.01, 0.005, 0.0025]  # list of mesh sizes (h), for mesh sensitivity

# Output directory and project file

# file's name
prj_name = "Kregime_Static.prj"
meshname = "mesh_full_pf"

from pathlib import Path

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
out_dir.mkdir(parents=True)

# Mesh generation

def mesh_generation(lc, lc_fine):
"""
lc ... characteristic length for coarse meshing
lc_fine ... characteristic length for fine meshing
"""
L = 4.0  # Length
H = 4.0  # Height
b = 0.4  # Length/Height of subdomain with fine mesh

# Before using any functions in the Python API, Gmsh must be initialized
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)

# Dimensions
dim1 = 1
dim2 = 2

# Points
gmsh.model.geo.addPoint(-L / 2, -H / 2, 0, lc, 1)
gmsh.model.geo.addPoint(L / 2, -H / 2, 0, lc, 2)
gmsh.model.geo.addPoint(L / 2, H / 2, 0, lc, 3)
gmsh.model.geo.addPoint(-L / 2, H / 2, 0, lc, 4)
gmsh.model.geo.addPoint(-b, -b - lc_fine / 2, 0, lc_fine, 5)
gmsh.model.geo.addPoint(b, -b - lc_fine / 2, 0, lc_fine, 6)
gmsh.model.geo.addPoint(b, b + lc_fine / 2, 0, lc_fine, 7)
gmsh.model.geo.addPoint(-b, b + lc_fine / 2, 0, lc_fine, 8)

# Lines

# Line loops
gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)
gmsh.model.geo.addCurveLoop([5, 6, 7, 8], 2)

# Add plane surfaces defined by one or more curve loops.

gmsh.model.geo.synchronize()

# Prepare structured grid
gmsh.model.geo.mesh.setTransfiniteCurve(
6, math.ceil(2 * b / lc_fine + 2), "Progression", 1
)
gmsh.model.geo.mesh.setTransfiniteCurve(
8, math.ceil(2 * b / lc_fine + 2), "Progression", 1
)
gmsh.model.geo.mesh.setTransfiniteSurface(2, "Alternate")

gmsh.model.geo.mesh.setRecombine(dim2, 1)
gmsh.model.geo.mesh.setRecombine(dim2, 2)

gmsh.model.geo.synchronize()

# Physical groups
Bottom = gmsh.model.addPhysicalGroup(dim1, [1])
gmsh.model.setPhysicalName(dim1, Bottom, "Bottom")

Right = gmsh.model.addPhysicalGroup(dim1, [2])
gmsh.model.setPhysicalName(dim1, Right, "Right")

Top = gmsh.model.addPhysicalGroup(dim1, [3])
gmsh.model.setPhysicalName(dim1, Top, "Top")

Left = gmsh.model.addPhysicalGroup(dim1, [4])
gmsh.model.setPhysicalName(dim1, Left, "Left")

Computational_domain = gmsh.model.addPhysicalGroup(dim2, [1, 2])
gmsh.model.setPhysicalName(dim2, Computational_domain, "Computational_domain")
gmsh.model.geo.synchronize()

output_file = f"{out_dir}/" + meshname + ".msh"
gmsh.model.mesh.generate(dim2)
gmsh.write(output_file)
gmsh.finalize()

# Pre process

def pre_processing(h, a0):
phase_field = np.ones((len(mesh.points), 1))
pv.set_plot_theme("document")

for node_id, x in enumerate(mesh.points):
if (
(mesh.center[0] - x[0]) <= a0 + 0.001 * h
and (mesh.center[0] - x[0]) >= -a0 - 0.001 * h
and (mesh.center[1] - x[1]) < h / 2 + 0.001 * h
and (mesh.center[1] - x[1]) > -h / 2 - 0.001 * h
):
phase_field[node_id] = 0.0

mesh.point_data["phase-field"] = phase_field
mesh.save(f"{out_dir}/mesh_full_pf_OGS_pf_ic.vtu")

# Run the simulation

import pyvista as pv

pv.set_plot_theme("document")
if "PYVISTA_HEADLESS" in os.environ:
pv.start_xvfb()
pv.set_jupyter_backend("static")

def sneddon_numerical(h):
# mesh properties
ls = 2 * h
# generate prefix from properties
filename = f"results_h_{h:0.4f}"
mesh_generation(0.1, h)
# Convert GMSH (.msh) meshes to VTU meshes appropriate for OGS simulation.
input_file = f"{out_dir}/" + meshname + ".msh"  # noqa: F841
! msh2vtu --keep_ids -o {out_dir} {input_file}
# As a preprocessing step, define the initial phase-field (crack).
pre_processing(h, a0)
# change properties in prj file #For more information visit: https://github.com/joergbuchwald/ogs6py
model = ogs.OGS(
INPUT_FILE=prj_name,
PROJECT_FILE=f"{out_dir}/{prj_name}",
MKL=True,
args=f"-o {out_dir}",
)

gml_file = Path("./Kregime_Static.gml").resolve()

model.replace_parameter_value(name="ls", value=ls)
model.replace_text(gml_file, xpath="./geometry")
model.replace_text(filename, xpath="./time_loop/output/prefix")
model.write_input()
# run simulation with ogs
t0 = time.time()
print(">>> OGS started execution ... <<<")

!ogs {out_dir}/{prj_name} -o {out_dir} -m {out_dir} > {out_dir}/ogs-out.txt
assert _exit_code == 0  # noqa: F821
tf = time.time()
print(">>> OGS terminated execution  <<< Elapsed time: ", round(tf - t0, 2), " s.")
for h_j in h_list:
sneddon_numerical(h=h_j)
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 30%] Meshing curve 3 (Line)
Info    : [ 40%] Meshing curve 4 (Line)
Info    : [ 60%] Meshing curve 5 (Line)
Info    : [ 70%] Meshing curve 6 (Line)
Info    : [ 80%] Meshing curve 7 (Line)
Info    : [ 90%] Meshing curve 8 (Line)
Info    : Done meshing 1D (Wall 0.000670768s, CPU 0.001s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)

Info    : [  0%] Blossom: 24032 internal 482 closed

Info    : [  0%] Blossom recombination completed (Wall 0.479341s, CPU 0.472833s): 7878 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 0.793024, min Q = 0.371315
Info    : [ 60%] Meshing surface 2 (Transfinite)
Info    : Done meshing 2D (Wall 0.72184s, CPU 0.709041s)
Info    : 14439 nodes 14848 elements
Info    : Writing '/var/lib/gitlab-runner/builds/e3EQ9HiK/0/ogs/build/release-petsc/Tests/Data/PhaseField/kregime_jupyter_notebook/Kregime_Static_jupyter/mesh_full_pf.msh'...
Info    : Done writing '/var/lib/gitlab-runner/builds/e3EQ9HiK/0/ogs/build/release-petsc/Tests/Data/PhaseField/kregime_jupyter_notebook/Kregime_Static_jupyter/mesh_full_pf.msh'

INFO:root:Output: mesh_full_pf

INFO:root:Original mesh (read)
INFO:root:<meshio mesh object>
Number of points: 14439
Number of cells:
line: 40
line: 40
line: 40
line: 40
Cell sets: Bottom, Right, Top, Left, Computational_domain, gmsh:bounding_entities
Point data: gmsh:dim_tags
Cell data: gmsh:physical, gmsh:geometrical
Field data: Bottom, Right, Top, Left, Computational_domain
INFO:root:Detected mesh dimension: 2
INFO:root:##

INFO:root:Domain mesh (written)
INFO:root:14439 points in 3 dimensions
INFO:root:point_data=['gmsh:dim_tags']
INFO:root:cell_data=['gmsh:physical']
INFO:root:cell_sets=[]
INFO:root:##

INFO:root:Boundary mesh (written)
INFO:root:160 points in 3 dimensions
INFO:root:cells: 160 line
INFO:root:point_data=['gmsh:dim_tags']
INFO:root:cell_data=['gmsh:physical']
INFO:root:cell_sets=[]
INFO:root:##
INFO:root:Submesh Bottom (written)
INFO:root:41 points in 3 dimensions
INFO:root:cells: 40 line
INFO:root:point_data=['gmsh:dim_tags']
INFO:root:cell_data=['gmsh:physical']
INFO:root:cell_sets=[]
INFO:root:##
INFO:root:Submesh Right (written)
INFO:root:41 points in 3 dimensions
INFO:root:cells: 40 line
INFO:root:point_data=['gmsh:dim_tags']
INFO:root:cell_data=['gmsh:physical']
INFO:root:cell_sets=[]
INFO:root:##
INFO:root:Submesh Top (written)
INFO:root:41 points in 3 dimensions
INFO:root:cells: 40 line
INFO:root:point_data=['gmsh:dim_tags']
INFO:root:cell_data=['gmsh:physical']
INFO:root:cell_sets=[]
INFO:root:##
INFO:root:Submesh Left (written)
INFO:root:41 points in 3 dimensions
INFO:root:cells: 40 line
INFO:root:point_data=['gmsh:dim_tags']
INFO:root:cell_data=['gmsh:physical']
INFO:root:cell_sets=[]
INFO:root:##

INFO:root:Submesh Computational_domain (written)
INFO:root:14439 points in 3 dimensions
INFO:root:point_data=['gmsh:dim_tags']
INFO:root:cell_data=['gmsh:physical']
INFO:root:cell_sets=[]
INFO:root:##

>>> OGS started execution ... <<<

>>> OGS terminated execution  <<< Elapsed time:  6.98  s.


# Post processing

# As a post-process, we calculate the fracture opening.
def width_calculation(filename):
# --------------------------------------------------------------------------------
# Active the results of last time step
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------
# define grad d
# --------------------------------------------------------------------------------
mesh_d = mesh.compute_derivative(scalars="phasefield")

"""A helper method to label the gradients into a dictionary."""
keys = np.array(
)
keys = keys.reshape((2, 3))[:, : arr.shape[1]].ravel()

# --------------------------------------------------------------------------------
# define width at nodes
# --------------------------------------------------------------------------------
disp = mesh.point_data["displacement"]
num_points = disp.shape
Wnode = np.zeros(num_points[0])
for i, _x in enumerate(mesh.points):
u_x = disp[i][0]
u_y = disp[i][1]

Wnode[i] = 0.5 * (u_x * gd_x + u_y * gd_y)
mesh.point_data["Wnode"] = Wnode
# --------------------------------------------------------------------------------
# define width at nodes
# --------------------------------------------------------------------------------
# Normalize the vector
normal = [np.cos(Orientation), np.sin(Orientation), 0]
# Make points along that vector for the extent of your slices
point_a = mesh.center + np.array([1.5 * a0, 0, 0])
point_b = mesh.center + np.array([-1.5 * a0, 0, 0])
dist_a_b = ((point_b[0] - point_a[0]) ** 2 + (point_b[1] - point_a[1]) ** 2) ** 0.5
# Define the line/points for the slices
line = pv.Line(point_a, point_b, n_slices)
width_line = np.zeros(len(line.points))
r_i = np.zeros(len(line.points))
# Generate all of the slices
slices = pv.MultiBlock()

for i, point in enumerate(line.points):
slices.append(mesh.slice(normal=normal, origin=point))
slice_mesh = mesh.slice(normal=normal, origin=point)
y_slice = slice_mesh.points[:, 1]
Wnode_slice = slice_mesh.point_data["Wnode"]
width_i = np.trapz(Wnode_slice, x=y_slice)  # noqa: NPY201
if width_i >= 0:
width_line[i] = width_i
r_i[i] = (
(point[0] - point_a[0]) ** 2 + (point[1] - point_a[1]) ** 2
) ** 0.5 - dist_a_b / 2

return r_i, width_line

# Sneddon

def sneddon(h, ls, phasefield_model):
# Effective a for AT1/A2
if phasefield_model == "AT1":
a_eff = a0 * (1 + pi * ls / (4.0 * a0 * (3 * h / 8.0 / ls + 1.0)))
elif phasefield_model == "AT2":
a_eff = a0 * (1 + pi * ls / (4.0 * a0 * (h / (2.0 * ls) + 1.0)))

x = np.linspace(-1.0 * a_eff, a_eff, 40)
uy = []
for i in range(len(x)):
uy.append(
2 * a_eff * (1 - nu**2) * P / E * math.sqrt(1.0 - ((x[i]) / (a_eff)) ** 2)
)

return x, uy, a_eff

## Opening profile

color = ["-.k", "ko", "-.r", "ro", "-.b", "bo", "-.g", "go"]
Label = ["Closed form solution", "VPF-FEM"]
lineWIDTH = [1.5, 1.5, 1.5]

for j, h_j in enumerate(h_list):
r_i_num = []
width_line_num = []
filename = f"results_h_{h_j:0.4f}"
print(filename)
width_calculation(filename)
r_i_num = width_calculation(filename)[0]
width_line_num = width_calculation(filename)[1]
ls = 2 * h_j

x_Sneddon = sneddon(h_j, ls, phasefield_model)[0]
uy_Sneddon = sneddon(h_j, ls, phasefield_model)[1]
a_eff_Sneddon = sneddon(h_j, ls, phasefield_model)[2]
plt.plot(
np.array(x_Sneddon[:]),
np.array(uy_Sneddon[:]),
color[2 * j],
fillstyle="none",
markersize=0,
linewidth=lineWIDTH[0],
label=f"Closed form solution - $h=${h_j:0.4f}",
)
plt.plot(
np.array(r_i_num[:]),
np.array(width_line_num[:]),
color[2 * j + 1],
fillstyle="none",
markersize=6,
linewidth=lineWIDTH[1],
label=f"VPF - $h =${h_j:0.4f}",
)
# ------------------------------------------------------------------------
plt.rcParams["figure.figsize"] = [40, 20]
plt.rcParams["figure.dpi"] = 1600
plt.ylabel("$w_n$ [m]", fontsize=14)
plt.xlabel("$r$ [m]", fontsize=14)
plt.grid(linestyle="dashed")
plt.title(f"{phasefield_model}")

legend = plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
plt.show()
results_h_0.0100


To decrease computing time, we perform the simulation with an one coarse mesh; Mesh sensitivity can be investigated in order to assess the convergence of the opening profile for various mesh discretizations. Following are the mesh sensivity results for for Models $\texttt{AT}_1$ and $\texttt{AT}_2$ ($h=0.01,~0.005,~\text{and} ~0.0025$ m.)

Our mesh size sensitivity study shows even with coarse mesh. The results away from crack tips agree well with the Sneddon solution. Using finer mesh demonstrates that the opening along the whole crack length is accurate compared with the closed-form solution.

It’s worth noting that we estimate width as part of post-processing (using the line integral and a given crack normal vector). However, near crack tips the crack normal vector is different to the given normal vector, so the width values are inaccurate near crack tips. For the sake of simplicity, we neglect determining the normal crack vector here; for more information on the line integral approach, see Yoshioka et al., 2020.

## References

[1] Yoshioka, Keita, Francesco Parisio, Dmitri Naumov, Renchao Lu, Olaf Kolditz, and Thomas Nagel. Comparative verification of discrete and smeared numerical approaches for the simulation of hydraulic fracturing. GEM-International Journal on Geomathematics 10, no. 1 (2019): 1-35.

[2] Yoshioka, Keita, Dmitri Naumov, and Olaf Kolditz. On crack opening computation in variational phase-field models for fracture. Computer Methods in Applied Mechanics and Engineering 369 (2020): 113210.

[3] Sneddon, Ian Naismith, and Morton Lowengrub. Crack problems in the classical theory of elasticity. 1969, 221 P (1969).

This article was written by Mostafa Mollaali, Keita Yoshioka. If you are missing something or you find an error please let us know.
Generated with Hugo 0.122.0 in CI job 475374 | Last revision: December 6, 2022