Skip to content
Snippets Groups Projects
Commit 3c3765f8 authored by Matthes's avatar Matthes Committed by Lars Bilke
Browse files

jupyter notebook phasefield beam

parent 7edf5b9f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 662 additions and 44 deletions
ProcessLib/PhaseField/Tests.cmake
+ 5
− 0
View file @ 3c3765f8
......@@ -137,4 +137,9 @@ AddTest(
if(OGS_USE_PETSC)
NotebookTest(NOTEBOOKFILE PhaseField/surfing_jupyter_notebook/surfing_pyvista.ipynb RUNTIME 25 RESOURCE_LOCK PYVISTA)
NotebookTest(NOTEBOOKFILE PhaseField/beam_jupyter_notebook/beam.ipynb RUNTIME 500 RESOURCE_LOCK PYVISTA)
endif()
if(TEST nb-PhaseField/beam_jupyter_notebook/beam-LARGE)
set_tests_properties(nb-PhaseField/beam_jupyter_notebook/beam-LARGE PROPERTIES PROCESSORS 3)
endif()
Tests/Data/PhaseField/beam_jupyter_notebook/beam.ipynb 0 → 100644
+ 373
− 0
View file @ 3c3765f8
%% Cell type:raw id:90c995c2 tags:
notebook = "PhaseField/beam_jupyter_notebook/beam_pyvista.ipynb"
author = "Matthes Kantzenbach, Keita Yoshioka, Mostafa Mollaali"
date = "2022-11-28"
title = "Beam"
web_subsection = "phase-field"
<!--eofm-->
%% Cell type:code id:41d0a5bc tags:
``` python
from ogs6py import ogs
import os
import ogs6py
import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import time
from xml.dom import minidom
from types import MethodType
```
%% Cell type:markdown id:aa987d95 tags:
## Problem description
In this example, it is shown, how phasefield models can be applied to simulate crack initiation.
Consider a bar $\Omega=[0,1]\times [0,0.05]\times [0,0.05]$ where the left end $\partial \Omega_{Left}$ is fixed and a displacement $U_x(t)$ is applied to the right end $\partial \Omega_{Right}$, resulting in a uniaxial state of stress $\boldsymbol{\sigma} = \sigma_{33} \, \mathbf{n}_3 \otimes \mathbf{n}_3$. At a certain point, fracture will initiate resulting in failure of the beam. $\newline$ This will be examined for positive $U_x(t)$ (tensile stress) and for negative $U_x(t)$ (compressive stress). As energy-split models the Isotropic and the Volumetric Deviatoric are used and as phasefield models $\texttt{AT}_1$ and $\texttt{AT}_2$. So in total 8 different cases are tested.
The aim is to estimate the ultimate load before fracture initiates for each case, as this can be veryfied against analytical results.
![Schematic view of beam](./figures/bar_.png#one-half "Schematic view of beam.")
Also considering homogeneous damage development (i.e.$\nabla v = 0$), the analytical solution for the strengths under uniaxial tension and compression are as follows.
| | Isotropic | Vol-Dev |
|:---:|:---:|:---:|
|Tensile $$\texttt{AT}_1$$| $$\sqrt{ \dfrac{3 G_\mathrm{c} E} {8 \ell}}$$| $$\sqrt{\dfrac{3 G_\mathrm{c} E}{8 \ell}}$$ |
|Tensile $$\texttt{AT}_2$$| $$ \dfrac{3}{16} \sqrt{ \dfrac{3 G_\mathrm{c} E } { \ell}}$$ | $$\dfrac{3}{16} \sqrt{\dfrac{3 G_\mathrm{c} E}{\ell}}$$ |
|Compressive $$\texttt{AT}_1$$| $$-\sqrt{ \dfrac{3 G_\mathrm{c} E} {8 \ell}}$$ | $$- \sqrt{ \dfrac{9 G_\mathrm{c} E} {16 \ell(1+\nu)}}$$ |
|Compressive $$\texttt{AT}_2$$| $$ -\dfrac{3}{16} \sqrt{ \dfrac{3 G_\mathrm{c} E } { \ell}}$$ | $$ -\dfrac{9}{32} \sqrt{\dfrac{2 G_\mathrm{c} E}{ \ell (1+\nu)}}$$ |
%% Cell type:markdown id:7ac26dee tags:
## Define some helper functions
%% Cell type:code id:b137ce19 tags:
``` python
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 method to be assigned to model, to replace a specific curve, given by name
# (analogue to replace_parameter method)
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 method to change timstepping 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]])
```
%% Cell type:markdown id:d7d09975 tags:
## Define function generating mesh, modifying project file and running ogs with given parameters
%% Cell type:code id:a34e61bd tags:
``` python
def ogs_beam(phasefield_model, energy_split_model, mesh_size = 0.01, length_scale = 0.02, bc_displacement = 5, ts_coords='0 0.05 1', values ='0 0.25 1', repeat_list=None, delta_t_list=None, hypre = False):
##phasefield_model: 'AT1' or 'AT2'
##energy_split_model: 'VolumetricDeviatoric' or 'Isotropic'
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'
#file's name
prj_name = "beam.prj"
print(f"> Running beam model {phasefield_model} – {energy_split_model} ... <")
logfile = f"{out_dir}/log_{phasefield_model}_{energy_split_model}.txt"
#beam dimensions
beam_height = 0.05
beam_depth = beam_height
beam_length = 1.
#mesh properties
h = mesh_size # distance between nodes
ls = length_scale
#generate prefix from properties
if energy_split_model == 'VolumetricDeviatoric':
prefix = phasefield_model + '_vd'
elif energy_split_model == 'Isotropic':
prefix = phasefield_model + '_iso'
else :
raise ValueError('"'+energy_split_model+'"'+' is no valid input for energy_split_model, choose between "VolumetricDeviatoric" and "Isotropic"')
if bc_displacement > 0:
prefix = prefix + '_tensile'
else :
prefix = prefix + '_compressive'
#generate mesh
! generateStructuredMesh -o {out_dir}/bar_.vtu -e hex --lx {beam_length} --nx {round(beam_length/h)} --ly {beam_height} --ny {round(beam_height/h)} --lz {beam_depth} --nz {round(beam_depth/h)} > {logfile}
! NodeReordering -i {out_dir}/bar_.vtu -o {out_dir}/bar.vtu >> {logfile}
! ExtractSurface -i {out_dir}/bar.vtu -o {out_dir}/bar_left.vtu -x 1 -y 0 -z 0 >> {logfile}
! ExtractSurface -i {out_dir}/bar.vtu -o {out_dir}/bar_right.vtu -x -1 -y 0 -z 0 >> {logfile}
! partmesh -s -o {out_dir} -i {out_dir}/bar.vtu >> {logfile}
! partmesh -m -n 3 -o {out_dir} -i {out_dir}/bar.vtu -- {out_dir}/bar_right.vtu {out_dir}/bar_left.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=ls)
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(prefix, xpath="./time_loop/output/prefix")
model.replace_parameter_value(name="dirichlet_right", value=bc_displacement)
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.write_input()
#run ogs
t0 = time.time()
print(" > OGS started execution ...")
! mpirun -n 3 ogs {out_dir}/{prj_name} -o {output_dir} >> {logfile}
tf = time.time()
print(" > OGS terminated execution. Elapsed time: ", round(tf - t0, 2), " s.")
```
%% Cell type:markdown id:324f380a tags:
## Input data
The values of the material properties are choosen as follows.
| **Name** | **Value** | **Unit** | **Symbol** |
|--------------------------------|--------------------|--------------|------------|
| _Young's modulus_ | 1 | Pa | $E$ |
| _Critical energy release rate_ | 1 | Pa$\cdot$m | $G_{c}$ |
| _Poisson's ratio_ | 0.15 | $-$ | $\nu$ |
| _Regularization parameter_ | 3$h$ | m | $\ell$ |
| _Length_ | $1$ | m | $L$ |
| _Height_ | $0.05$ | m | $H$ |
| _Depth_ | $0.05$ | m | $D$ |
%% Cell type:markdown id:02482e8b tags:
## Run Simulations
%% Cell type:code id:71c6b710 tags:
``` python
pf_ms = ['AT1', 'AT2']
es_ms = ['VolumetricDeviatoric','Isotropic']
displ = [4.,-4.]
'''
for a in pf_ms:
for c in displ:
ogs_beam(a,es_ms[1],bc_displacement = c,mesh_size = 0.01, length_scale = 0.03)
ogs_beam(pf_ms[1],es_ms[0],bc_displacement = 4,mesh_size = 0.01, length_scale = 0.03)
# run AT1_vd_tensile with smaller timesteps in critical time range
ogs_beam(pf_ms[0],es_ms[0],bc_displacement = 5,mesh_size = 0.01, length_scale = 0.03,repeat_list=['62','2','20','1'], delta_t_list=['1e-2','1e-3','1e-4','1e-2'])
# run VolumetricDeviatoric in Compression with Hypre and smaller timesteps in critical time range
ogs_beam(pf_ms[1],es_ms[0],bc_displacement = -4.5,mesh_size = 0.01, length_scale = 0.03, hypre = True, repeat_list=['70','4','30','1'], delta_t_list=['1e-2','1e-3','1e-4','1e-2'])
# loosen relative error tolerance for displacement process in order to get convergence for the AT1 case
prj_path='./'
prj_name = "beam.prj"
model = ogs.OGS(INPUT_FILE=prj_path+prj_name, PROJECT_FILE=prj_path+prj_name, MKL=True)
model.replace_text('1e-6', xpath="./time_loop/processes/process/convergence_criterion/reltol",occurrence=0)
model.write_input()
ogs_beam(pf_ms[0],es_ms[0],bc_displacement = -4.95, mesh_size = 0.01, length_scale = 0.03, hypre= True, repeat_list=['66', '8','3','3','20','1'], delta_t_list=['1e-2','1e-3','1e-4','1e-5','1e-6','1e-2'],ts_coords='0 0.1 1', values ='0 0.5 1')
model = ogs.OGS(INPUT_FILE=prj_path+prj_name, PROJECT_FILE=prj_path+prj_name, MKL=True)
model.replace_text('1e-14', xpath="./time_loop/processes/process/convergence_criterion/reltol",occurrence=0)
model.write_input()
'''
## run only cases easy to handle with coarse timestepping:
for a in pf_ms:
for b in es_ms:
for c in displ:
if a == 'AT1' and b == 'VolumetricDeviatoric':
continue
if a == 'AT2' and b == 'VolumetricDeviatoric' and c < 0:
ogs_beam(a,b,bc_displacement = c,mesh_size = 0.01,length_scale = 0.03, hypre = True, repeat_list=['1'],delta_t_list=['1e-1'])
else:
ogs_beam(a,b,bc_displacement = c,mesh_size = 0.01,length_scale = 0.03, repeat_list=['1'],delta_t_list=['1e-1'])
```
%% Cell type:markdown id:44016656 tags:
## Results
The final timestep of the AT1 iso tensile case is shown exemplary for the phasefield after fracture.
![Phasefield after fracture](./figures/beam_final.png "phasefield after fracture")
%% Cell type:markdown id:d2b7386b tags:
## Post-processing
The force applied to the beam is compared to the change in length of the beam.
%% Cell type:code id:b1357ab9 tags:
``` python
# define function to obtain displacement applied at the right end of the beam from rvtu file
def displ_right(filename):
data = pv.read(filename)
data.point_data["displacement"]
max_x = max(data.points[:,0])
return np.mean(data.point_data["displacement"][:,0], where= np.transpose(data.points[:,0]==max_x))
# define fuction to obtain force acting on the right end of the beam from vtu file
def force_right(filename):
data = pv.read(filename)
data.point_data["NodalForces"]
max_x = max(data.points[:,0])
return np.sum(data.point_data["NodalForces"][:,0], where= np.transpose(data.points[:,0]==max_x))
# define function applying obove functions on all vtu file listet in pvd file, returning force-displacement curve
def force_displ_from_pvd(pvd):
doc = minidom.parse(pvd)
DataSets = doc.getElementsByTagName("DataSet")
vtu_files = [x.getAttribute("file") for x in DataSets]
forces_right = [force_right(f"{out_dir}/{x}") for x in vtu_files]
displs_right = [displ_right(f"{out_dir}/{x}") for x in vtu_files]
return [forces_right,displs_right]
```
%% Cell type:markdown id:5375302a tags:
## Plot force-strain curves
%% Cell type:code id:5ec1718d tags:
``` python
prefixes = ['AT1_vd_tensile','AT1_iso_tensile','AT2_vd_tensile','AT2_iso_tensile', 'AT1_vd_compressive', 'AT1_iso_compressive', 'AT2_vd_compressive', 'AT2_iso_compressive']
labels = [r'$\texttt{AT}_1$ vol-dev tensile',r'$\texttt{AT}_1$ iso tensile',r'$\texttt{AT}_2$ vol-dev tensile',r'$\texttt{AT}_2$ iso tensile', r'$\texttt{AT}_1$ vol-dev compressive', r'$\texttt{AT}_1$ iso compressive', r'$\texttt{AT}_2$ vol-dev compressive', r'$\texttt{AT}_2$ iso compressive']
ls=['-','--']
colors = ['#ffdf4d','#006ddb','#8f4e00','#ff6db6','#920000','#b66dff','#db6d00','#490092']
fig, ax = plt.subplots()
plt.rc('text', usetex=True)
fig.set_size_inches(18.5, 10.5)
for i,pre in enumerate(prefixes):
pvd = f"{output_dir}/{pre}.pvd"
if os.path.isfile(pvd) :
curve = force_displ_from_pvd(pvd)
ax.plot(curve[1],curve[0],ls[i%2], label = labels[i],linewidth=5, color = colors[i], alpha= 1)
plt.rcParams['xtick.labelsize'] = 16
plt.rcParams['ytick.labelsize'] = 16
ax.grid(linestyle='dashed')
ax.set_xlabel('$\Delta [m]$',fontsize =18)
ax.set_ylabel('$F_y [N]$',fontsize =18)
plt.legend(fontsize =18, ncol = 2)
ax.axhline(y = 0, color = 'black',linewidth=1)
ax.axvline(x = 0, color = 'black',linewidth=1)
ax.set_xlim(-4.5, 4.5);
```
%% Cell type:markdown id:9346910f tags:
Running all the cases with finer timesteps yields the following figure. The failure of the beam in the simulation is observed, when the loading reaches the analytical strength for the corresponding case.
![Results](./figures/beam_result_with_analytical.png "Results")
Tests/Data/PhaseField/beam_jupyter_notebook/beam.prj 0 → 100644
+ 284
− 0
View file @ 3c3765f8
<?xml version='1.0' encoding='ISO-8859-1'?>
<OpenGeoSysProject>
<meshes>
<mesh>bar.vtu</mesh>
<mesh>bar_left.vtu</mesh>
<mesh>bar_right.vtu</mesh>
</meshes>
<processes>
<process>
<name>PhaseField</name>
<type>PHASE_FIELD</type>
<coupling_scheme>staggered</coupling_scheme>
<integration_order>2</integration_order>
<phasefield_model>AT1</phasefield_model>
<energy_split_model>Isotropic</energy_split_model>
<irreversible_threshold>0.05</irreversible_threshold>
<constitutive_relation>
<type>LinearElasticIsotropic</type>
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
<phasefield_parameters>
<residual_stiffness>k</residual_stiffness>
<crack_resistance>gc</crack_resistance>
<crack_length_scale>ls</crack_length_scale>
</phasefield_parameters>
<solid_density>rho_sr</solid_density>
<process_variables>
<phasefield>phasefield</phasefield>
<displacement>displacement</displacement>
</process_variables>
<secondary_variables>
<secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
<secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
</secondary_variables>
<specific_body_force>0 0 0</specific_body_force>
</process>
</processes>
<time_loop>
<global_process_coupling>
<max_iter> 1000 </max_iter>
<convergence_criteria>
<!-- convergence criterion for the first process -->
<convergence_criterion>
<type>DeltaX</type>
<norm_type>INFINITY_N</norm_type>
<abstol>1.e-1</abstol>
<reltol>1.e-1</reltol>
</convergence_criterion>
<!-- convergence criterion for the second process -->
<convergence_criterion>
<type>DeltaX</type>
<norm_type>INFINITY_N</norm_type>
<abstol>1.e-4</abstol>
<reltol>1.e-10</reltol>
</convergence_criterion>
</convergence_criteria>
</global_process_coupling>
<processes>
<process ref="PhaseField">
<nonlinear_solver>basic_newton_u</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<reltol>1e-14</reltol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial>0</t_initial>
<t_end>1.0</t_end>
<timesteps>
<pair>
<repeat>1</repeat>
<delta_t>1e-1</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
<process ref="PhaseField">
<nonlinear_solver>petsc_snes</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<reltol>1.e-14</reltol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial>0</t_initial>
<t_end>1.0</t_end>
<timesteps>
<pair>
<repeat>1</repeat>
<delta_t>1e-1</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
</processes>
<output>
<variables>
<variable>displacement</variable>
<variable>phasefield</variable>
<variable>sigma</variable>
<variable>epsilon</variable>
</variables>
<type>VTK</type>
<prefix>AT1_iso_compressive</prefix>
<timesteps>
<pair>
<repeat>1</repeat>
<each_steps>1</each_steps>
</pair>
</timesteps>
</output>
</time_loop>
<parameters>
<!-- Mechanics -->
<parameter>
<name>E</name>
<type>Constant</type>
<value>1.0</value>
</parameter>
<parameter>
<name>nu</name>
<type>Constant</type>
<value>0.15</value>
</parameter>
<parameter>
<name>k</name>
<type>Constant</type>
<value>1e-8</value>
</parameter>
<parameter>
<name>gc</name>
<type>Constant</type>
<value>1.0</value>
</parameter>
<parameter>
<name>ls</name>
<type>Constant</type>
<value>0.03</value>
</parameter>
<parameter>
<name>rho_sr</name>
<type>Constant</type>
<value>0.0</value>
</parameter>
<parameter>
<name>displacement0</name>
<type>Constant</type>
<values>0 0 0</values>
</parameter>
<parameter>
<name>phasefield_ic</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>phasefield_bc</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>dirichlet0</name>
<type>Constant</type>
<value>0</value>
</parameter>
<parameter>
<name>Dirichlet_spatial</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>dirichlet_right_time</name>
<type>CurveScaled</type>
<curve>dirichlet_time</curve>
<parameter>dirichlet_right</parameter>
</parameter>
<parameter>
<name>dirichlet_right</name>
<type>Constant</type>
<value>-4.0</value>
</parameter>
</parameters>
<curves>
<curve>
<name>dirichlet_time</name>
<coords>0 0.05 1</coords>
<values>0 0.25 1</values>
</curve>
</curves>
<process_variables>
<process_variable>
<name>phasefield</name>
<components>1</components>
<order>1</order>
<initial_condition>phasefield_ic</initial_condition>
<boundary_conditions>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>left</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>phasefield_bc</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>phasefield_bc</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>PhaseFieldIrreversibleDamageOracleBoundaryCondition</type>
<component>0</component>
</boundary_condition>
</boundary_conditions>
</process_variable>
<process_variable>
<name>displacement</name>
<components>3</components>
<order>1</order>
<initial_condition>displacement0</initial_condition>
<boundary_conditions>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>left</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>dirichlet0</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>dirichlet_right_time</parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
</process_variables>
<nonlinear_solvers>
<nonlinear_solver>
<name>petsc_snes</name>
<type>PETScSNES</type>
<max_iter>200</max_iter>
<linear_solver>linear_solver_d</linear_solver>
</nonlinear_solver>
<nonlinear_solver>
<name>basic_newton_u</name>
<type>Newton</type>
<max_iter>10000</max_iter>
<linear_solver>linear_solver_u</linear_solver>
</nonlinear_solver>
</nonlinear_solvers>
<linear_solvers>
<linear_solver>
<name>linear_solver_d</name>
<eigen>
<solver_type>BiCGSTAB</solver_type>
<precon_type>ILUT</precon_type>
<max_iteration_step>10000</max_iteration_step>
<error_tolerance>1e-16</error_tolerance>
</eigen>
<petsc>
<parameters>-ksp_type cg -pc_type bjacobi -ksp_atol 1e-10 -ksp_rtol 1e-10 -snes_type vinewtonrsls -snes_linesearch_type l2 -snes_atol 1.e-8 -snes_rtol 1.e-8 -snes_max_it 1000 -snes_monitor</parameters>
</petsc>
</linear_solver>
<linear_solver>
<name>linear_solver_u</name>
<petsc>
<parameters>-ksp_type cg -pc_type bjacobi -ksp_atol 1e-14 -ksp_rtol 1e-14</parameters>
</petsc>
</linear_solver>
</linear_solvers>
</OpenGeoSysProject>
Tests/Data/PhaseField/beam_jupyter_notebook/figures/bar_.png 0 → 100644
+ 0
− 0
View file @ 3c3765f8
Tests/Data/PhaseField/beam_jupyter_notebook/figures/bar_.png

10.5 KiB

Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_final.png 0 → 100644
+ 0
− 0
View file @ 3c3765f8
Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_final.png

47 KiB

Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_result_with_analytical.png 0 → 100644
+ 0
− 0
View file @ 3c3765f8
Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_result_with_analytical.png

144 KiB

web/content/docs/benchmarks/phase-field/phasefield/Miao_Biot2017.pdf deleted 100755 → 0
+ 0
− 0
View file @ 7edf5b9f
File deleted
web/content/docs/benchmarks/phase-field/phasefield/beam.png deleted 100644 → 0
+ 0
− 0
View file @ 7edf5b9f
web/content/docs/benchmarks/phase-field/phasefield/beam.png

11.7 KiB

web/content/docs/benchmarks/phase-field/phasefield/beam_d.png deleted 100644 → 0
+ 0
− 0
View file @ 7edf5b9f
web/content/docs/benchmarks/phase-field/phasefield/beam_d.png

16.3 KiB

web/content/docs/benchmarks/phase-field/phasefield/beam_u.png deleted 100644 → 0
+ 0
− 0
View file @ 7edf5b9f
web/content/docs/benchmarks/phase-field/phasefield/beam_u.png

17.3 KiB

web/content/docs/benchmarks/phase-field/phasefield/error_u.png deleted 100644 → 0
+ 0
− 0
View file @ 7edf5b9f
web/content/docs/benchmarks/phase-field/phasefield/error_u.png

58.8 KiB

web/content/docs/benchmarks/phase-field/phasefield/index.md deleted 100644 → 0
+ 0
− 44
View file @ 7edf5b9f
+++
author = "Xing-Yuan Miao"
date = "2017-05-19T09:10:33+01:00"
title = "Crack beam under tension"
image = "beam_u.png"
+++
{{< data-link >}}
## Problem description
We solve a homogeneous beam model under a given displacement loading. The length of the beam is 2\,mm. Detailed model description can refer [this PDF](Miao_Biot2017.pdf).
## Results and evaluation
Results show crack Phase-Field and displacement field distributions through the length of the beam:
{{< img src="beam.png" >}}
{{< img src="beam_d.png" >}}
{{< img src="beam_u.png" >}}
For highlighting the deviation between the analytical and numerical solution, we provide the absolute error of the analytical solution and numerical simulation as follows:
{{< img src="error_u.png" >}}
The analytical solution is:
$$
\begin{equation}
d (x) = 1 - {\mathrm{e}}^{\frac{- |x|}{2 \varepsilon}}
\end{equation}
$$
$$
\begin{equation}
u (x) = \dfrac{\sigma}{E} \int_0^x \dfrac{1}{d (x)^2 + k} \mathrm{d}x
\end{equation}
$$
with
$$
\begin{equation}
\sigma = \dfrac{E u_0}{I (\varepsilon, k)}
\end{equation}
$$
\begin{equation}
I (\varepsilon, k) = \int_0^1 \dfrac{1}{d (x)^2 + k} \mathrm{d}x
\end{equation}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment