Hey guys,
Im currently working on my bachelor thesis and i started with coding a linear elasticity problem with a phase field to detect a crack on my domain. Unfortunately the phi in the calculation is not evolving the way I want it to. In the picture above you can see that the phi evolves more like a surface than a crack and I do not really know why. My code is based on this paper by Hirshikesh et al. The mesh is made with Gmsh, I will also post the gmsh script below. Maybe some of you have a clue how to fix this code, I would be very grateful for any help.
Paper: A FEniCS implementation of the phase field method for
quasi-static brittle fracture
Regards, Leandro
**GMSH:**
import gmsh
gmsh.initialize()
gmsh.model.add("testmesh3")
lc = 5
width=1
gmsh.model.geo.addPoint(0,0,0, lc)
gmsh.model.geo.addPoint(90,0,0, lc)
gmsh.model.geo.addPoint(100,0,42, lc)
gmsh.model.geo.addPoint(110,0,0, lc)
gmsh.model.geo.addPoint(200,0,0, lc)
gmsh.model.geo.addPoint(200,0,50, lc)
gmsh.model.geo.addPoint(0,0,50, lc)
gmsh.model.geo.addLine(1,2)
gmsh.model.geo.addLine(2,3)
gmsh.model.geo.addLine(3,4)
gmsh.model.geo.addLine(4,5)
gmsh.model.geo.addLine(5,6)
gmsh.model.geo.addLine(6,7)
gmsh.model.geo.addLine(7,1)
gmsh.model.geo.addCurveLoop([1,2,3,4,5,6,7])
field_tag = gmsh.model.mesh.field.add("MathEval")
# middle of the profile, mesh should be more fine
x_feiner = 100
z_feiner = 47
lc=gmsh.model.mesh.field.setString(field_tag, "F", f"0.3 + 0.03 * sqrt((x - {x_feiner})^2 + (z - {z_feiner})^2)")
# activate field
gmsh.model.mesh.field.setAsBackgroundMesh(field_tag)
gmsh.model.geo.synchronize()
gmsh.model.geo.addPlaneSurface([1])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=1)
gmsh.model.setPhysicalName(2, 1, "front")
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.Algorithm", 5) # or 5/7/11
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.recombine()
# Extrude the 2D mesh to create a 3D mesh
vol = gmsh.model.geo.extrude([(2, 1)], 0, breite, 0, #2,1 -> dim, tag und dim=2 weil 0=punkt, 1= linie ...
[1], recombine=True)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(2, [19], tag=2)
gmsh.model.setPhysicalName(2, 2, "bottom_left")
gmsh.model.addPhysicalGroup(2, [31], tag=3)
gmsh.model.setPhysicalName(2, 3, "bottom_right")
gmsh.model.addPhysicalGroup(2, [35], tag=4)
gmsh.model.setPhysicalName(2, 4, "right_side")
gmsh.model.addPhysicalGroup(2, [39], tag=6)
gmsh.model.setPhysicalName(2, 6, "top")
gmsh.model.addPhysicalGroup(2, [43], tag=5)
gmsh.model.setPhysicalName(2, 5, "left_side")
gmsh.model.addPhysicalGroup(2, [44], tag=7)
gmsh.model.setPhysicalName(2, 7, "backside")
gmsh.model.addPhysicalGroup(3, [1], tag=8)
gmsh.model.setPhysicalName(3, 8, "domain")
# Generate the 3D mesh
gmsh.model.mesh.generate(3)
gmsh.write("mesh_linearelastFINALRadial.msh")
gmsh.fltk.run()
gmsh.finalize()
import numpy as np
from dolfinx import log, default_scalar_type,fem,io
from dolfinx import fem, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from dolfinx.nls.petsc import NewtonSolver
import ufl
from numpy import random
from petsc4py import PETSc
from mpi4py import MPI
import pyvista
from dolfinx.io import gmshio
from dolfinx.io import XDMFFile
#Mesh von Gmsh reinladen
domain, cell_tags, facet_tags = gmshio.read_from_msh("mesh_linearelastFINALRadial.msh", MPI.COMM_WORLD)
#Functionspaces
ndegree=1
FS= fem.functionspace(domain, ("CG", ndegree, (domain.geometry.dim, )))
FSPF=fem.functionspace(domain, ("CG", ndegree))
#Phasefield "phi"
phi = fem.Function(FSPF)
phi.x.array[:]= 0
#Boundary Conditions
disp= fem.Constant(domain, PETSc.ScalarType(0.0))
zero=np.zeros(domain.geometry.dim, dtype=default_scalar_type)
zero_scalar=PETSc.ScalarType(0.0)
dofs_bottom_left = fem.locate_dofs_topological(FS.sub(2), domain.topology.dim - 1, facet_tags.find(2))
dofs_bottom_right = fem.locate_dofs_topological(FS.sub(2), domain.topology.dim - 1, facet_tags.find(3))
dofsleft = fem.locate_dofs_topological(FS.sub(0), domain.topology.dim-1, facet_tags.find(5))
dofsright = fem.locate_dofs_topological(FS.sub(0), domain.topology.dim-1, facet_tags.find(4))
dofsfront = fem.locate_dofs_topological(FS.sub(1), domain.topology.dim-1, facet_tags.find(1))
bc_bottom_l = fem.dirichletbc(zero_scalar, dofs_bottom_left, FS.sub(2))
bc_bottom_r = fem.dirichletbc(zero_scalar, dofs_bottom_right, FS.sub(2))
bc_right = fem.dirichletbc(disp, dofsright, FS.sub(0))
bc_front = fem.dirichletbc(zero_scalar, dofsfront, FS.sub(1))
bc_left = fem.dirichletbc(zero_scalar, dofsleft, FS.sub(0))
BC= [bc_bottom_l, bc_bottom_r, bc_right, bc_front, bc_left ]
#Constants
E = 210000
nu = 0.2
mu = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
# lambda_=110.5 # 4000
# mu=82
Gc=0.03
l=0.6
tol= 1e-4
diff= 1
phi_old = fem.Function(FSPF)
phi_old.x.array[:] = phi.x.array
#Testfunktionen und Verschiebung
u=ufl.TrialFunction(FS)
v=ufl.TestFunction(FS)
dphi = ufl.TrialFunction(FSPF) # TrialFunction für Phasefield
theta = ufl.TestFunction(FSPF)
#Funktionen und Konstanten
I= ufl.Identity(len(u))
f= fem.Constant(domain, default_scalar_type((0.0, 0.0,0))) #Bodyforce
T = fem.Constant(domain, default_scalar_type((0, 0, 0))) #Traction
H= fem.Function(FSPF)
#Integrationen
metadata = {"quadrature_degree": 2}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tags, metadata=metadata) # to integrate over the surface
dx = ufl.Measure("dx", domain=domain, metadata=metadata) # to integrate over the domain
def sigma(u):
eps_reg= 1e-5
return (lambda_ * ufl.nabla_div(u)* I + mu*(ufl.grad(u)+ufl.grad(u).T))*((1-phi)**2+ eps_reg)
def epsilon(u):
return 0.5*(ufl.grad(u)+ ufl.grad(u).T)
#L
L = ufl.inner(f,v)*dx+ufl.inner(T, v)*ds
#a(u,v)
a = ufl.inner(sigma(u), epsilon(v))*dx
dispmax= 1 # 20
delta= 0.01
loadfactors = np.arange(delta,1+delta,delta)
Traction_list =[]
Disp_List =[]
normal = ufl.FacetNormal(domain)
xdmf_phi = XDMFFile(MPI.COMM_WORLD, "phi_output2.xdmf", "w")
xdmf_phi.write_mesh(domain)
# xdmf_phi.write_function(phi, 0.0)
xdmf_uh = XDMFFile(MPI.COMM_WORLD, "uh_output2.xdmf", "w")
xdmf_uh.write_mesh(domain)
counter=0
for loadfactor in loadfactors:
diff=1.0
disp.value = loadfactor * dispmax
counter = counter + 1
print(f"Iteration:{counter}")
while diff > tol:
problem = LinearProblem(a, L, bcs=BC, petsc_options={"ksp_type": "cg","pc_type": "hypre","ksp_rtol": 1e-8})
uh = problem.solve()
psi_expression = fem.Expression(0.5*lambda_*(ufl.tr(epsilon(uh)))**2 + mu * (ufl.inner(epsilon(uh), epsilon(uh))),FSPF.element.interpolation_points())
psi= fem.Function(FSPF)
psi.interpolate(psi_expression)
with H.vector.localForm() as H_local, psi.vector.localForm() as psi_local:H_local.array[:] = np.maximum(H_local.array, psi_local.array)
#a fürs phasefield
a_pf= ufl.inner(ufl.grad(theta), Gc * l * ufl.grad(dphi)) * dx + (Gc / l + 2 * H) * theta * dphi * dx
#L fürs phasefield
L_pf=2*H*theta*dx
PFproblem = LinearProblem(a_pf, L_pf, petsc_options={"ksp_type": "cg", "pc_type": "hypre", "ksp_rtol": 1e-8})
phi = PFproblem.solve()
#Konvergenz?
diff = np.linalg.norm(phi.x.array - phi_old.x.array)
phi_old.x.array[:] = phi.x.array
print(f"diff = {diff:.4e}")
Disp_List.append(loadfactor * dispmax)
#Reaktionskraft zwischenspeichern
storedtraction = fem.assemble_scalar(fem.form(ufl.inner(ufl.as_vector([1.0, 0.0, 0.0]), ufl.dot(sigma(uh),normal))*ds(4))) #top surface traction
Traction_list.append(storedtraction)
xdmf_phi.write_function(phi, loadfactor)
xdmf_uh.write_function(uh, loadfactor)
xdmf_phi.close()
xdmf_uh.close()
#pdf für Kraft-Verschiebung
from matplotlib import pyplot as plt
plt.plot(Disp_List, Traction_list, '-')
plt.xlabel("Displacement")
plt.ylabel("Force")
plt.legend(["Kraft-Verschiebungs Kurve"])
plt.grid(True)
plt.show()
please note that i changed my code from
L_pf=2*H*theta*dx + ufl.inner(ufl.grad(phi), ufl.FacetNormal(domain))*theta*ds
to
L_pf=2*H*theta*dx
otherwise phi did not convergate for a small Gc
