 Phase Field Fracture

Dear community,
I am using the following code for a research, but I did not understant how read and check the results in Paraview. Can someone help me?
The code was made by Hirshikesh, S. Natarajan, R. K. Annabattula, E. Martinez-Paneda.

from dolfin import *
mesh = Mesh('mesh.xml')
# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
p, q = TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)

# Introduce manually the material parameters
Gc =  2.7
l = 0.015
lmbda = 121.1538e3
mu = 80.7692e3

# Constituive functions
def epsilon(u):
def sigma(u):
return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))
def psi(u):
return 0.5*(lmbda+mu)*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2+mu*inner(dev(epsilon(u)),dev(epsilon(u)))
def H(uold,unew,Hold):
return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)

# Boundary conditions
top = CompiledSubDomain("near(x, 0.5) && on_boundary")
bot = CompiledSubDomain("near(x, -0.5) && on_boundary")
def Crack(x):
return abs(x) < 1e-03 and x <= 0.0
load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bc_u = [bcbot, bctop]
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)

# Variational form
unew, uold = Function(W), Function(W)
pnew, pold, Hold = Function(V), Function(V), Function(V)
*inner(p,q)-2.0*H(uold,unew,Hold)*q)*dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi), rhs(E_phi), pnew, bc_phi)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)

# Initialization of the iterative procedure and output requests
t = 0
u_r = 0.007
deltaT  = 0.1
tol = 1e-3
conc_f = File ("./ResultsDir/phi.pvd")
fname = open('ForcevsDisp.txt', 'w')

# Staggered scheme
while t<=1.0:
t += deltaT
if t >=0.7:
deltaT = 0.0001
iter = 0
err = 1
while err > tol:
iter += 1
solver_disp.solve()
solver_phi.solve()
err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
err_phi = errornorm(pnew,pold,norm_type = 'l2',mesh = None)
err = max(err_u,err_phi)
uold.assign(unew)
pold.assign(pnew)
Hold.assign(project(psi(unew), WW))
if err < tol:
print ('Iterations:', iter, ', Total time', t)
if round(t*1e4) % 10 == 0:
conc_f << pnew
Traction = dot(sigma(unew),n)
fy = Traction*ds(1)
fname.write(str(t*u_r) + "\t")
fname.write(str(assemble(fy)) + "\n")

fname.close()
print ('Simulation completed')

Caio, I have worked with this exact paper before. On that note I think you have posted the code twice in your example above. I replaced the mesh.xml with RectangleMesh() and made some other minor modifications, otherwise our codes are the same.

Could you clarify your question? This specific implementation only has two outputs, a .pvd of the crack phase field and .txt of the Load vs. Displacement. I assume you want to visualize the crack movement?

Sven, I just ran the code again and it seems it was working, but it was taking a while. So I remove this part

if t >=0.7:
deltaT = 0.0001

However, from the 0.8 second, the number of iterations stopped appearing and these messages were repeating infinitely:

Solving linear variational problem. Warning: Degree of exact solution may be inadequate for accurate result in errornorm.

Does the code really take long? Should I not remove that part of the code?

Caio,

Short answer is yes, it takes a long time and don’t remove that line.

For this post’s reference the paper can be found here. As is common with many coupled phase field approaches, a staggered approach is taken here. A simpler problem like this could be solved with a coupled solution - but it would still take a long time.

The best and worst part of the phase field method is that it requires very small steps when a phenomena is active. I.e. in this paper the authors choose inputs (Gc, l, ...) such that the necessary energy to propagate the crack is not reach until a load in the vicinity of t=0.7. This implementation uses load driven displacement load.t=t*u_r so the load is tied directly to the current time. There is no reason to take tiny steps prior to the initiation of crack movement hence:

if t >=0.7:
deltaT = 0.0001

Before this time the problem is essentially linear elastic in nature. Once crack propagation starts the solution can drastically change with just small changes in the applied load. I recommend that you go through each line and make sure you are aware of the author’s choices. The reason you stop seeing updates with large time steps is:

err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
err_phi = errornorm(pnew,pold,norm_type = 'l2',mesh = None)

If you take large steps while the crack is propagating, the method will not converge and in this specific code you will be stuck in an infinite loop of those messages (no max iterations). In my experience with staggered approaches you take large time steps until the solution fails to converge and then refine the timestep at that time. If you changed any of the inputs from this paper, you would likely need to change where the step refinement occurs. Phase field methods can be very stable but at the cost of small timesteps (and no lack of fidgeting). Consider how quickly a crack can initiate and propagate in a real material - and now we want to actually capture it moving in small steps. Hope this helps.

1 Like

I got it, the code is running now. Thanks a lot, Sven.

Hello, please tell me, I also have the same code and I want to use it for a problem in dynamic and mixed mode. Can you tell me what modifications you have made. Thanks

Dear Sven,

I wanted to run the code faster, so I modified the mesh to have fewer elements, but the result was this: different from the paper, which is like this: I made the mesh with tihs code:

from dolfin import *
from mshr import *
domain=Rectangle(Point(-0.5,-0.5),Point(0.5,0.5))-Rectangle(Point(-0.5,-0.0001),Point(0,0.0001))
mesh=generate_mesh(domain,10)

The code used was the same of the paper, I just modified the mesh.
Can you show me how to solve the problem?

Caio,

Again I advise you to consider the outcomes you want in the context of the physics at play. The last word that comes to my mind when Phase Field is discussed is “fast”. If you have used the same code then the following should be true:

l = 0.015

This is the critical length of the crack, the driving factor of the “phase” or smeared crack approach that the paper is based on. Generally your spatial discretization should be on the order of the phase field parameter, otherwise how can you actually diffuse a crack over multiple nodes? Consider the following:

This is the Author’s mesh: