# Three point bending test

## Problem description

The notched beam under three-point bending is simulated to show how the phase field models can capture crack propagation and the structural strength. The geometry and boundary conditions of the model are depicted below. The bottom left point is fixed on both $x$ and $y$ directions and bottom right point is fixed on the $y$ direction. A displacement boundary condition $u^{\ast}$ is applied on the midpoint of the upper edge. The crack will initiate at the notch tip and propagate upward to the upper edge.

• Three phase field models are implemented in OGS, i.e., $\texttt{COHESIVE}$, $\texttt{AT}_1$ and $\texttt{AT}_2$.
• The energy-split models include EffectiveStress, Isotropic and the VolumetricDeviatoric.
• The softening rules implemented for the $\texttt{COHESIVE}$ model are Linear and Exponential.

The force-displacement curves ($F^{\ast}$-$u^{\ast}$) for various sets of the above models will be recorded and compared with the experimental data.

## Define some helper functions

import os
import shutil
import time
from pathlib import Path
from types import MethodType
from xml.dom import minidom

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyvista as pv
from ogs6py import ogs
data_dir = os.environ.get("OGS_DATA_DIR", "../../..")

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

output_dir = out_dir

# define function to replace a specific curve, given by name

def replace_curve(
self,
name=None,
value=None,
coords=None,
parametertype=None,
valuetag="values",
coordstag="coords",
):
root = self._get_root()
parameterpath = "./curves/curve"
parameterpointer = self._get_parameter_pointer(root, name, parameterpath)
self._set_type_value(parameterpointer, value, parametertype, valuetag=valuetag)
self._set_type_value(parameterpointer, coords, parametertype, valuetag=coordstag)

# define function to change time_stepping in project file

def set_timestepping(model, repeat_list, delta_t_list):
model.remove_element(
xpath="./time_loop/processes/process/time_stepping/timesteps/pair"
)
for i in range(len(repeat_list)):
blocktag="pair",
parent_xpath="./time_loop/processes/process/time_stepping/timesteps",
taglist=["repeat", "delta_t"],
textlist=[repeat_list[i], delta_t_list[i]],
)

## Define functions generating mesh, modifying project file and running ogs with given parameters

def ogs_TPB(
phasefield_model,
energy_split_model,
softening_curve="Linear",
length_scale=5.0,
bc_displacement=-1.0,
ts_coords="0 1.0",
values="0 1.0",
repeat_list=None,
delta_t_list=None,
hypre=False,
MPI=True,
ncores=4,
):
## define input file

without_hypre = "-ksp_type cg -pc_type bjacobi -ksp_atol 1e-14 -ksp_rtol 1e-14"
with_hypre = "-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_strong_threshold 0.7 -ksp_atol 1e-8 -ksp_rtol 1e-8"

prj_name = "TPB.prj"
print(
f"> Running three point bending test {phasefield_model} - {energy_split_model} - {softening_curve} ... <"
)
logfile = f"{out_dir}/log_{phasefield_model}_{energy_split_model}.txt"  # noqa: F841
model = ogs.OGS(INPUT_FILE=prj_name, PROJECT_FILE=f"{out_dir}/{prj_name}", MKL=True)
# generate prefix from properties
prefix = f"{phasefield_model}" + f"_{energy_split_model}"
if phasefield_model == "COHESIVE":
prefix = (
f"{phasefield_model}" + f"_{energy_split_model}" + f"_{softening_curve}"
)

if MPI:
# partition mesh
! NodeReordering -i TPB.vtu -o {out_dir}/TPB.vtu >> {logfile}
! constructMeshesFromGeometry -m {out_dir}/TPB.vtu -g TPB.gml >> {logfile}
shutil.move("TPB_left.vtu", f"{out_dir}/TPB_left.vtu")
shutil.move("TPB_right.vtu", f"{out_dir}/TPB_right.vtu")
shutil.move("TPB_top.vtu", f"{out_dir}/TPB_top.vtu")
! partmesh -s -o {out_dir} -i {out_dir}/TPB.vtu >> {logfile}
! partmesh -m -n {ncores} -o {out_dir} -i {out_dir}/TPB.vtu -- {out_dir}/TPB_right.vtu {out_dir}/TPB_left.vtu {out_dir}/TPB_top.vtu >> {logfile}
else:
! NodeReordering -i TPB.vtu -o {out_dir}/TPB.vtu >> {logfile}
# change properties in prj file
model = ogs.OGS(INPUT_FILE=prj_name, PROJECT_FILE=f"{out_dir}/{prj_name}", MKL=True)
model.replace_parameter_value(name="ls", value=length_scale)
model.replace_text(phasefield_model, xpath="./processes/process/phasefield_model")
model.replace_text(
energy_split_model, xpath="./processes/process/energy_split_model"
)
model.replace_text(softening_curve, xpath="./processes/process/softening_curve")
model.replace_text(prefix, xpath="./time_loop/output/prefix")
model.replace_curve = MethodType(replace_curve, model)
model.replace_curve(name="dirichlet_time", value=values, coords=ts_coords)
if repeat_list is not None and delta_t_list is not None:
set_timestepping(model, repeat_list, delta_t_list)
else:
set_timestepping(model, ["1"], ["1e-2"])
if hypre is True:
model.replace_text(
with_hypre,
xpath="./linear_solvers/linear_solver/petsc/parameters",
occurrence=1,
)
else:
model.replace_text(
without_hypre,
xpath="./linear_solvers/linear_solver/petsc/parameters",
occurrence=1,
)
model.replace_text("./TPB.gml", xpath="./geometry")
model.write_input()
# run ogs
t0 = time.time()
if MPI:
print(f"  > OGS started execution with MPI - {ncores} cores...")
! mpirun --bind-to none -np {ncores} ogs {out_dir}/{prj_name} -o {output_dir} >> {logfile}
else:
print("  > OGS started execution ...")
! ogs {out_dir}/{prj_name} -o {output_dir} >> {logfile}
tf = time.time()
print("  > OGS terminated execution. Elapsed time: ", round(tf - t0, 2), " s.")

## Input data

The general parameters are chosen as follows.

Name Value Unit Symbol
Young’s modulus 20000 MPa $E$
Poisson’s ratio 0.2 $-$ $v$
Critical energy release rate 0.113 N/mm $G_{c}$
Regularization parameter 5 mm $\ell$
Minimum element size 1 mm $h$

The material parameters for the COHESIVE model are as follows.

Name Value Unit Symbol
Tensile strength 2.4 MPa $f_t$
Irwin’s length 392.4 mm $\ell_{ch}$

## Run Simulations

In the following, coarse time steps are chosen for the sake of less computational cost.

# phasefield_model = ["COHESIVE",'AT1', 'AT2']
# energy_split_model = ["EffectiveStress",'VolumetricDeviatoric','Isotropic']
# softening_curve = ['Line','Exponential']

disp = -1.0
ls = 5.0
mpi_cores = 4  # MPI cores
## only run selected cases
# For the COHESIVE model
for a in ["Linear", "Exponential"]:
ogs_TPB(
"COHESIVE",
"EffectiveStress",
a,
length_scale=ls,
bc_displacement=disp,
repeat_list=["1"],
delta_t_list=["1e-1"],
ncores=mpi_cores,
)

# For AT1 and AT2 models with isotropic split
for b in ["AT1", "AT2"]:
ogs_TPB(
b,
"Isotropic",
length_scale=ls,
bc_displacement=disp,
repeat_list=["1"],
delta_t_list=["1e-1"],
ncores=mpi_cores,
)
> Running three point bending test COHESIVE - EffectiveStress - Linear ... <

  > OGS started execution with MPI - 4 cores...

--------------------------------------------------------------------------
prterun detected that one or more processes exited with non-zero status,
thus causing the job to be terminated. The first process to do so was:

Process name: [prterun-envinf2-1712074@1,1] Exit code:    1
--------------------------------------------------------------------------

  > OGS terminated execution. Elapsed time:  1.72  s.
> Running three point bending test COHESIVE - EffectiveStress - Exponential ... <

  > OGS started execution with MPI - 4 cores...

--------------------------------------------------------------------------
prterun detected that one or more processes exited with non-zero status,
thus causing the job to be terminated. The first process to do so was:

Process name: [prterun-envinf2-1712635@1,2] Exit code:    1
--------------------------------------------------------------------------

  > OGS terminated execution. Elapsed time:  1.69  s.
> Running three point bending test AT1 - Isotropic - Linear ... <

  > OGS started execution with MPI - 4 cores...

--------------------------------------------------------------------------
prterun detected that one or more processes exited with non-zero status,
thus causing the job to be terminated. The first process to do so was:

Process name: [prterun-envinf2-1713073@1,1] Exit code:    1
--------------------------------------------------------------------------

  > OGS terminated execution. Elapsed time:  1.75  s.
> Running three point bending test AT2 - Isotropic - Linear ... <

  > OGS started execution with MPI - 4 cores...

--------------------------------------------------------------------------
prterun detected that one or more processes exited with non-zero status,
thus causing the job to be terminated. The first process to do so was:

Process name: [prterun-envinf2-1713472@1,2] Exit code:    1
--------------------------------------------------------------------------

  > OGS terminated execution. Elapsed time:  1.73  s.


## Results

The phase field profile in the final time step of the COHESIVE EffectiveStress case with Exponential softening is shown below.

## Post-processing

The force-displacement curves ($F^{\ast}$-$u^{\ast}$) applied to the beam are compared.

# define function to obtain displacement applied at the loading point from vtu file
def displ_midpoint(filename):
return np.sum(
data.point_data["displacement"][:, 1],
where=(data.points[:, 0] == 225) * (data.points[:, 1] == 100),
)

# define function to obtain force at the loading point from vtu file

def force_midpoint(filename):
return (
np.sum(
data.point_data["NodalForces"][:, 1],
where=(data.points[:, 0] == 225) * (data.points[:, 1] == 100),
)
/ 10.0
)

# define function to apply the above functions on all vtu files listed in pvd file, returning force-displacement curves

def force_displ_from_pvd(pvd):
doc = minidom.parse(pvd)
DataSets = doc.getElementsByTagName("DataSet")
vtu_files = [x.getAttribute("file") for x in DataSets]
forces_midpoint = [force_midpoint(f"{out_dir}/{x}") for x in vtu_files]
displs_midpoint = [displ_midpoint(f"{out_dir}/{x}") for x in vtu_files]
return [forces_midpoint, displs_midpoint]

## Plot force-displacement curves

# Load experimental data
data_upper = pd.read_csv("figures/experiment_data_upper_limit.csv")
prefixes = [
"AT1_Isotropic",
"AT2_Isotropic",
"COHESIVE_EffectiveStress_Linear",
"COHESIVE_EffectiveStress_Exponential",
]
labels = [
r"${AT}_1$ Isotropic",
r"${AT}_2$ Isotropic",
r"${COHESIVE}$ EffectiveStress Linear",
r"${COHESIVE}$ EffectiveStress Exponential",
]
ls = ["-.", "--", ".", "-"]
colors = ["#ffdf4d", "#006ddb", "#8f4e00", "#ff6db6"]

fig, ax = plt.subplots()
for i, pre in enumerate(prefixes):
pvd = f"{output_dir}/{pre}.pvd"
if Path(pvd).is_file():
curve = force_displ_from_pvd(pvd)
plt.plot(
curve[1],
curve[0],
ls[i % 2],
label=labels[i],
linewidth=2,
color=colors[i],
alpha=1,
)

plt.rcParams["xtick.labelsize"] = 12
plt.rcParams["ytick.labelsize"] = 12
ax.grid(linestyle="dashed")
ax.set_xlabel("$u^{\\ast}$ [mm]", fontsize=12)
ax.set_ylabel("$F^{\\ast}$ [kN]", fontsize=12)
plt.xlim(plt.xlim()[::-1])
plt.ylim(plt.ylim()[::-1])
plt.legend(ncol=1)
ax.axhline(y=0, color="black", linewidth=1)
ax.axvline(x=0, color="black", linewidth=1)
plt.fill_between(
data_upper.iloc[:, 0], 0, data_upper.iloc[:, 1], facecolor="green", alpha=0.3
)
plt.fill_between(
data_lower.iloc[:, 0], 0, data_lower.iloc[:, 1], facecolor="white", alpha=1
)
/tmp/ipykernel_1710661/59277585.py:38: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
plt.legend(ncol=1)

<matplotlib.collections.PolyCollection at 0x7e12733e1520>


Running all the selected cases with fine time steps (100 steps) yields the following figure in which the green patch bounds the experimental data. The cohesive models match the experimental results better than $\texttt{AT}_1$ and $\texttt{AT}_2$.

This article was written by Tao You, Keita Yoshioka, Mostafa Mollaali. 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: November 28, 2022