## Static fracture opening under a constant pressure – Sneddon solution

This benchmark is available as a Jupyter notebook: PhaseField/kregime_jupyter_notebook/Kregime_Static_jupyter.ipynb.

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

pi = math.pi
plt.rcParams["text.usetex"] = True

## 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{<pressurized_crack_scheme>static</pressurized_crack_scheme>}$$ 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$$
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"

out_dir = os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out")
os.makedirs(out_dir, exist_ok=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

# 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
gmsh.model.setPhysicalName(dim1, Bottom, "Bottom")

gmsh.model.setPhysicalName(dim1, Right, "Right")

gmsh.model.setPhysicalName(dim1, Top, "Top")

gmsh.model.setPhysicalName(dim1, Left, "Left")

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")
pv.start_xvfb()
pv.set_jupyter_backend("static")

def sneddon_numerical(h):
#mesh properties
ls = 2*h
#generate prefix from properties
filename =  'results_h_%0.4f'%h
mesh_generation(0.1, h)
# Convert GMSH (.msh) meshes to VTU meshes appropriate for OGS simulation.
input_file = f"{out_dir}/"+meshname+".msh"
! msh2vtu --ogs -o {out_dir}/{meshname} {input_file}
# As a preprocessing step, define the initial phase-field (crack).
pre_processing(h,a0)
model = ogs.OGS(INPUT_FILE=prj_name, PROJECT_FILE=f"{out_dir}/{prj_name}", MKL=True, args=f"-o {out_dir}")
model.replace_parameter_value(name="ls", value=ls)
model.replace_text("./Kregime_Static.gml", 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} > {out_dir}/log.txt
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    : [ 50%] 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.000554099s, CPU 0.000808s)
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.3314s, CPU 0.331645s): 7878 quads, 0 triangles, 0 invalid quads, 0 quads with Q < 0.1, avg Q = 0.793024, min Q = 0.371315
Info    : [ 50%] Meshing surface 2 (Transfinite)
Info    : Done meshing 2D (Wall 0.486623s, CPU 0.484048s)
Info    : 14439 nodes 14848 elements
Info    : Writing '/var/lib/gitlab-runner/builds/vZ6vnZiU/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/vZ6vnZiU/0/ogs/build/release-petsc/Tests/Data/PhaseField/kregime_jupyter_notebook/Kregime_Static_jupyter/mesh_full_pf.msh'
Output: /var/lib/gitlab-runner/builds/vZ6vnZiU/0/ogs/build/release-petsc/Tests/Data/PhaseField/kregime_jupyter_notebook/Kregime_Static_jupyter/mesh_full_pf

14439 points in 3 dimensions; cells: 160 line, 14358 quad; point_data=['gmsh:dim_tags']; cell_data=['gmsh:physical', 'gmsh:geometrical']; cell_sets=['Bottom', 'Right', 'Top', 'Left', 'Computational_domain', 'gmsh:bounding_entities']
##
Detected mesh dimension: 2
##
Domain mesh (written)
14439 points in 3 dimensions; cells: 14358 quad; point_data=['original_node_number']; cell_data=['MaterialIDs']; cell_sets=[]
##
Boundary mesh (written)
160 points in 3 dimensions; cells: 160 line; point_data=['bulk_node_ids']; cell_data=['bulk_elem_ids']; cell_sets=[]
##
Submesh Bottom (written)
41 points in 3 dimensions; cells: 40 line; point_data=['bulk_node_ids']; cell_data=['bulk_elem_ids']; cell_sets=[]
##
Submesh Right (written)
41 points in 3 dimensions; cells: 40 line; point_data=['bulk_node_ids']; cell_data=['bulk_elem_ids']; cell_sets=[]
##
Submesh Top (written)
41 points in 3 dimensions; cells: 40 line; point_data=['bulk_node_ids']; cell_data=['bulk_elem_ids']; cell_sets=[]
##
Submesh Left (written)
41 points in 3 dimensions; cells: 40 line; point_data=['bulk_node_ids']; cell_data=['bulk_elem_ids']; cell_sets=[]
##
Submesh Computational_domain (written)
14439 points in 3 dimensions; cells: 14358 quad; point_data=['bulk_node_ids']; cell_data=['MaterialIDs']; cell_sets=[]
##
msh2vtu successfully finished
>>> OGS started execution ... <<<
>>> OGS terminated execution  <<< Elapsed time:  3.43  s.


# Post processing

# As a post-process, we calculate the fracture opening.
def width_calculation(filename):
# --------------------------------------------------------------------------------
# Active the results of last time step
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------
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)
x_slice = slice_mesh.points[:, 0]
y_slice = slice_mesh.points[:, 1]
Wnode_slice = slice_mesh.point_data["Wnode"]
width_i = np.trapz(Wnode_slice, x=y_slice)
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 = "results_h_%0.4f" % h_j
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="Closed form solution - $h=$%0.4f" % h_j,
)
plt.plot(
np.array(r_i_num[:]),
np.array(width_line_num[:]),
color[2 * j + 1],
fillstyle="none",
markersize=6,
linewidth=lineWIDTH[1],
label="VPF - $h =$%0.4f" % h_j,
)
# ------------------------------------------------------------------------
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("%s" % 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.101.0 in CI job 284016 | Last revision: December 6, 2022