Code for bend test using phase field fracture

variational formulation
1

I want to make a code bend test using phase field method for brittle fracture. This is the code I am making, following the idea of the paper (Hirshikesh et al):

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import math
# Mesh
vertices=[Point(0,0.4),Point(0.1,0),Point(4,0),Point(4,2),Point(-4,2),Point(-4,0),Point(-0.1,0)]
domain=Polygon(vertices)
mesh=generate_mesh(domain, 100)

# 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 = 0.54
l = 2*mesh.hmin()
lmbda =12e3 
mu = 8e3


# 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
class top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0,1e-2) and near(x[1],2,1e-2)
    
def bot1(x,on_boundary):
    return on_boundary and near(x[0],-4) and near(x[1],0)
    
def bot2(x,on_boundary):
    return on_boundary and near(x[0],4) and near(x[1],0)
    
load = Expression("t", t = 0.0, degree=1)
bcbot1= DirichletBC(W, Constant((0,0)), bot1)
bcbot2= DirichletBC(W, Constant((0,0)), bot2)
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot1, bcbot2, bctop]
bc_phi = []
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,2)
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.001
deltaT  = 0.1
tol = 1e-3
conc_f = File ("./Bend_Test/phi.pvd")
fname = open('ForcevsDisp_Bend_Test.txt', 'w')
x=[]
y=[]
# Staggered scheme
while t<=5:
    t += deltaT
    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")
                x.append(str(t*u_r))
                y.append(str(round(assemble(fy)*1e-3,1)))
                print("Force:", str(round(assemble(fy)*1e-3,1)),'kN')
                print("Displacement:",str(t*u_r),'mm')
               
fname.close()
print ('Simulation completed')

And this is the mesh, with different restrictions for the bottom:

image

So I got an error with the bctop:

RuntimeError: Invalid argument

Can anyone help me?

Paper

2

Caio,

Since you declare top as:

class top(SubDomain):

You either need to declare an instance of this class to use in the DirichletBC or use

bctop = DirichletBC(W.sub(1), load, top())

You will run into a similar issue with marking the boundaries

top.mark(boundaries,2)
top().mark(boundaries,2)
3

I changed this, but this message appers:

Warning: Found no facets matching domain for boundary condition

Why?

4

Caio,

The message says it all. I imagine you receive the message three times, none of your BCs return true. If you intend to apply the bot1 and bot2 to a single point you could specify as DirichletBC:

bcbot1= DirichletBC(W, Constant((0,0)), bot1, method="PointWise)

Otherwise on_boundary returns false for both cases. I highly recommend plotting your MeshFunction to verify your boundary is being applied correctly. If you do so you will see that top() doesn’t return any facets due to your mesh density and the tolerance your provided the function near().

5

Thanks for the reply, Sven. I plotted the MeshFunction and verified. It works, but now the problem is about the load.

For the bot it worked with “pointwise”. So I changed the top, because was giving the same message, and it works now.

 > class top(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and  abs(x[0]-0) < 0.2 and abs(x[1]-2) < 0.2

However, the load does not increase. I think it’s wrong in that line:

load = Expression("-t", t = 0.0, degree=1)

How can I apply a top load in the vertical direction to downward, which increases with time?

6

Caio,

Can you specify what you mean by:

The expression is fine, consider looking at:

  • Your material inputs, is the simulation prematurely fracturing?
  • Do your phase field results look correct visually?
  • Does your displacement field look correct?
  • Etc.

I apologize if I wasn’t specific before but you should remove on_boundary from your bot1 and bot2 expressions if you are using method=pointwise. That may be causing your issues.

1 Like
7

The problem is that it seems that the load is not being applied to the mesh, so I can’t even see results, because it is always zero.

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import math
# Mesh
vertices=[Point(0,0.4),Point(0.1,0),Point(4,0),Point(4,2),Point(-4,2),Point(-4,0),Point(-0.1,0)]
domain=Polygon(vertices)
mesh=generate_mesh(domain, 50)
#plot(mesh)
#plt.show()

# 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 = 0.54
l = 2*mesh.hmin()
lmbda =12e3 
mu = 8e3


# 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
class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and  abs(x[0]-0) < 0.2 and abs(x[1]-2) < 0.2
    
def bot1(x,on_boundary):
    return on_boundary and near(x[0],-3.9) and near(x[1],0)
    
def bot2(x,on_boundary):
    return on_boundary and near(x[0],3.9) and near(x[1],0)
    
load = Expression("-t", t = 0.0, degree=1)
bcbot1= DirichletBC(W, Constant((0,0)), bot1, method="pointwise")
bcbot2= DirichletBC(W, Constant((0,0)), bot2, method="pointwise")
bctop = DirichletBC(W.sub(1), load, top())
bc_u = [bcbot1, bcbot2, bctop]
bc_phi = []
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.001
deltaT  = 0.1
tol = 1e-3
conc_f = File ("./Bend_Test/phi.pvd")
fname = open('ForcevsDisp_Bend_Test.txt', 'w')
x=[]
y=[]
# Staggered scheme
while t<=5:
    t += deltaT
    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")
                x.append(str(t*u_r))
                y.append(str(round(assemble(fy)*1e-3,1)))
                print("Force:", str(round(assemble(fy),1)),'N')
                print("Displacement:",str(t*u_r),'mm')
               
fname.close()
print ('Simulation completed')

You can check this running the code, the load remains zero. that’s why I asked how can I apply a top load in the vertical direction to downward, which increases with time, because I think the mistake is in this.

8

Caio,

When you say that the load is remaining zero I assume you are referring to the output from your print functions and your output file. In my opinion, for any new code I develop I always make sure to constantly plot all of my functions to make sure they make numerical and physical sense. If you plot your displacement, u, say in ParaView and warp it to a vector - you will see your mesh fly away. This comes to a couple different points:

Why did you change

to

If you don’t provide near() a third argument it uses DOLFIN_EPS=3e-16. Again you are using method="pointwise" but have left on_boundary in your function. Your displacement boundary condition is functioning as intended, your other BCs are the issue.

9

Ok, Sven. Thanks

The fact that the mesh goes up when applying the force means that the load is up?. So I wanted to know how I apply the load down.

I’m sorry, I didn’t remove on_boundary before, but I did it. Do you think I need to do this for top() too?

Regarding the near() function, what tolerance do you indicate for my mesh? Considering that I am using, again:

def bot1(x,on_boundary):
    return on_boundary and near(x[0],-4) and near(x[1],0)
10

Caio,

The load expression is fine, in 2D your W.sub(1) is assigning the displacement in the y-direction where “up” is positive thus load = Expression("-t", t = 0.0, degree=1) is appropriate. I did not say that the mesh went up. If you warp the magnitude of the vector in ParaView it would go up of course, but plotting just the y-component of unew will show the load is correctly applied.

That’s up to you, is your load applied at a single point or across an area? Is it only applied on the surface (on_boundary) or on interior facets? I don’t know. I assume you will want to leave your expression as is.

I don’t know, that’s your call. You are using mshr to generate your meshes so you have less control over the mesh than GMSH would give you. If you leave it as is, it would work just fine because you have defined a point at that coordinate (the corners are clearly defined). Otherwise you would need to review your mesh coordinates to determine where your points of interest are (hence why return on_boundary and near(x[0],-3.9) and near(x[1],0) is not guaranteed to return any points). Hope that helps.

11

Sven,

I want to apply the load at a single point (0,2), on the surface. So I think I will leave my expression as is.

I have changed some things in the code, regarding the support points and the way the boundary conditions are applied. In paper the rightmost point of the mesh has vertical and horizontal restriction, and the leftmost has only vertical restriction. I changed the mesh to be equal of the paper.

class top(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0]-0) < 0.2 and abs(x[1]-2) < 0.2

def bot1(x):
    return near(x[0],4,0.2) and near(x[1],0,0.2)
    
def bot2(x):
    return near(x[0],4,0.2) and near(x[1],0,0.2)
    
load = Expression("-t", t = 0.0, degree=1)
bcbot1= DirichletBC(W, Constant((0,0)), bot1, method="pointwise")
bcbot2= DirichletBC(W.sub(1),0, bot2, method="pointwise")
bctop = DirichletBC(W.sub(1), load, top(), method="pointwise")

I ran the code (with a low displacement) just for viewing in Paraview, and noticed that the boundary conditions looked good.

image
image723Ă—229 2.46 KB

So I changed the remote displacement to 0.01 mm, because the crack propagation starts at 0.04, and this happened:

image

I know the mesh may be more refined. However, the maximum load in my code is higher than in the paper. For example, in the paper, there is the point (0.02,0.02) and in my code it is (0.016,0.02). Also, the fracture occurs near 0.05 mm, which does not occur in my model.

image

Whats is the problem?

1 Like
12

Caio,

Without reading through the paper and having your complete code that is a tough question to answer. Your code is running but it seems there might be some materiel inputs or implementation differences that are causing this. I’d check through each line and confirm you’ve replicated the literature as closely as possible.

Sven

13

Sven,

This is the complete code

from fenics import *
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import math
# Mesh
#vertices=[Point(0,0.4),Point(0.1,0),Point(4,0),Point(4,2),Point(-4,2),Point(-4,0),Point(-0.1,0)]
#vertices=[Point(0,1.6),Point(-0.1,2),Point(-4,2),Point(-4,0),Point(4,0),Point(4,2),Point(0.1,2)]
#domain=Polygon(vertices)
#mesh=generate_mesh(domain, 200)
mesh = Mesh('./Murilo/geom_julia.xml')
boundaries = MeshFunction('size_t', mesh,'./Murilo/geom_julia_facet_region.xml')
#plot(mesh)
#plt.show()

# 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 = 0.54
l = 0.03
lmbda =12e3 
mu = 8e3

# 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):
    str_ele = 0.5*(grad(u) + grad(u).T)
    IC = tr(str_ele)
    ICC = tr(str_ele * str_ele)
    return (0.5*lmbda*IC**2) + mu*ICC
    #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
class top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],4,3e-1) and near(x[1],2,3e-1)
    
def bot1(x):
    return near(x[0],0,1e-5) and near(x[1],0,1e-5)
    
def bot2(x):
    return near(x[0],8,1e-5) and near(x[1],0,1e-5)

load = Expression("-t", t = 0.0, degree=1)
bcbot1= DirichletBC(W, Constant((0,0)), bot1, method="pointwise")
bcbot2= DirichletBC(W.sub(1),0, bot2, method="pointwise")
bctop = DirichletBC(W.sub(1), load, top())
bc_u = [bcbot1, bcbot2, bctop]
bc_phi = []
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.01
deltaT  = 1
tol = 1e-3
conc_f = File ("./Bend_Test/phi.pvd")
fname = open('ForcevsDisp_Bend_Test.txt', 'w')
x=[]
y=[]
# Staggered scheme
while t<=11:
    t += deltaT
    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('-----------------')
            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(round(-assemble(fy)*1e-3,3)) + "\n")
                x.append(str(t*u_r))
                y.append(str(round(-assemble(fy)*1e-3,3)))
                print("Force:", str(-1e-3*assemble(fy)),'kN')
                print("Displacement:",str(t*u_r),'mm')
                print('-----------------')
               
fname.close()
print ('Simulation completed') 

#Deformação
eps=TensorFunctionSpace(mesh,"CG", degree=1)
deformação=Function(eps)
deformação.assign(project(epsilon(unew),eps))
File('Bend_Test/deformação.pvd')<<deformação

#TensĂŁo
Vsig = TensorFunctionSpace(mesh, "CG", degree=1)
tensões = Function(Vsig)
tensões.assign(project(sigma(unew),Vsig))
File('Bend_Test/tensões.pvd')<<tensões

#Deslocamento
File('Bend_Test/deslocamento.pvd')<<unew

#Gráfico 

plt.title('FORÇA X DESLOCAMENTO')
plt.xlabel('Displacement (mm)')
plt.ylabel('Force (kN)')
plt.plot(x,y)
plt.show()

I changed the mesh for a more refined one in the middle, to change the phase field lenght to be equal to the paper.

If you run the code, you will see that at the beginning (elastic zone) the results are close. But the problem is that there is no drop-off in my graph.

I would check the paper and confirm I have replicated the literature as closely as possible. One thing I noticed is that according to the change of the third parameter of near() function in

class top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],4,3e-1) and near(x[1],2,3e-1)

there is a clear change in the applied load.

14

Caio,

By applied load I assume you are referring to the loadfile output from:

fy = Traction[1]*ds(1)
fname.write(str(t*u_r) + "\t")
fname.write(str(round(-assemble(fy)*1e-3,3)) + "\n")

In your last post you were using method='pointwise', but you have removed that:

This coupled with fy = Traction[1]*ds(1) means your load scales with the size of ds(1) as you are prescribing the applied displacement through DirichletBC. If you aren’t going to be using point loading consider normalizing fy by the size of ds(1):

fy = Traction[1]*ds(1)
loadArea = Constant(1)*ds(1)
fname.write(str(round(-(assemble(fy)/assemble(loadArea))*1e-3,3)) + "\n")

Sven

15

Sven,

I add this in my code:

but this message appers:

ufl.log.UFLException: This integral is missing an integration domain.

I tried to use method='pointwise', to be a point loading

class top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],4,0.2) and near(x[1],2,0.2)
bctop = DirichletBC(W.sub(1), load, top(), method="pointwise")

but the same problem occurs. And if I use method='pointwise' in bctop and loadArea = Constant(1)*ds(1), the same message appers.

16

Caio,

Sorry, I usually use the following command for defining the boundaries:

ds = Measure("ds", domain = mesh, subdomain_data=boundaries)

Changing this will fix the issue. Again, I would defer to the literature to determine how the load is being applied and where they are measuring the force. Doing what I suggested will get you the average traction, not average force.

Sven

17

Sven,

Maybe I’m doing something very wrong. The paper is here. Then, here is the domain and boundary conditions:
image

I think the boundary conditions are being applied correctly (top,bc1 and bc2). The load-displacement curve is that:

image

If the average traction, not average force, how can I have average force equal to this curve?