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):
    return sym(grad(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[1], 0.5) && on_boundary")
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")
def Crack(x):
    return abs(x[1]) < 1e-03 and x[0] <= 0.0
load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(1), load, top)
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)
E_du = ((1.0-pold)**2)*inner(grad(v),sigma(u))*dx
E_phi = (Gc*l*inner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))\
            *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
    load.t=t*u_r
    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[1]*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?

Thanks for the advice, Sven

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:
image

different from the paper, which is like this:

image

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:

This is your mesh:

The perspective of the images is a little off but they are both unit square meshes. Consider the number of the nodes in the y direction to be roughly 11, that means you have a spacing of 0.1 between nodes as opposed to the critical length scale of 0.015. As you displayed above, you can see the mesh resolution in your solution thus your crack looks “odd”. Your density is not sufficient to create a diffuse crack, the value along the center is likely 1 or “damaged” while the next adjacent element is 0 or “undamaged” - hence the zigzag appearance.

Of note, a lot of the author’s meshing is unnecessary. With this starting crack along the middle and a tensile load applied in the y direction - the crack is going to propagate in nearly a straight line. If you make your own mesh (GMSH or other software), you could reduce a lot of the elements and focus that density around where the phenomena is actually occurring.

Additionally there is no real need to physically remove geometry from the mesh to generate the starting crack - this is what Initial Conditions are for in the Phase Field method. That is one of the convenient parts of the method (no remeshing or explicit crack front tracking). In your case it helps increase the mesh density around the crack but that is easily done in a meshing software.

If you are interested in using the Phase Field method I highly recommend reading other pieces of literature that can provide additional perspectives.

1 Like