## Three point bending test

This benchmark is available as a Jupyter notebook: PhaseField/TPB_jupyter_notebook/TPB.ipynb.

from ogs6py import ogs
import os
import shutil
import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import time
import pandas as pd
from xml.dom import minidom
from types import MethodType

## 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

data_dir = os.environ.get('OGS_DATA_DIR', '../../..')
out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')
if not os.path.exists(out_dir):
os.makedirs(out_dir)

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)):
model.add_block(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., bc_displacement = -1., 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"
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 != None and delta_t_list != None:
set_timestepping(model,repeat_list, delta_t_list)
else:
set_timestepping(model,['1'], ['1e-2'])
if hypre == 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("  > OGS started execution with MPI - " f"{ncores} cores...")
! mpirun -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.
ls = 5.
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...
> OGS terminated execution. Elapsed time:  23.92  s.
> Running three point bending test COHESIVE - EffectiveStress - Exponential ... <
> OGS started execution with MPI - 4 cores...
> OGS terminated execution. Elapsed time:  33.36  s.
> Running three point bending test AT1 - Isotropic - Linear ... <
> OGS started execution with MPI - 4 cores...
> OGS terminated execution. Elapsed time:  19.18  s.
> Running three point bending test AT2 - Isotropic - Linear ... <
> OGS started execution with MPI - 4 cores...
> OGS terminated execution. Elapsed time:  18.16  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(f"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 os.path.isfile(pvd) :
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)
<matplotlib.collections.PolyCollection at 0x7f87bdb07690>


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.101.0 in CI job 284016 | Last revision: November 28, 2022