Phase Field Fracture

Caio_Cesar · January 14, 2022, 1:28pm

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')

Sven · January 14, 2022, 1:47pm

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?

Caio_Cesar · January 14, 2022, 2:19pm

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

Sven · January 14, 2022, 3:14pm

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.

Caio_Cesar · January 14, 2022, 3:30pm

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

gilles_bidi · January 18, 2022, 11:39am

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

Caio_Cesar · January 20, 2022, 4:13pm

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?

Sven · January 20, 2022, 5:45pm

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.

Caio_Cesar · February 11, 2022, 11:02am

Sven,

The paper code worked well, thanks a lot for it.

Now, I’m working with an experimental test of fracture mecanichs. The only different thing is the mesh, but I made it using dolfin and mshr just to test if the code works. The boundary conditions are the same. This is the modified code for the test:

from fenics import *
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
vertices=[Point(8.5,0),Point(5.1,-0.9),Point(2.9,-0.9),Point(1.5,-3),
Point(0,-2),Point(0,-140),Point(22,-140),Point(22,140),
Point(0,140),Point(0,2),Point(1.5,3),Point(2.9,0.9),Point(5.1,0.9)]

domain=Polygon(vertices)
mesh=generate_mesh(domain, 200)

# 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 = 1
E=Constant(200e3)
nu=Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)


# 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],140) && on_boundary")
bot = CompiledSubDomain("near(x[1],-140) && 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 ("./Results/phi.pvd")
fname = open('ForcevsDisp.txt', 'w')

# Staggered scheme
while t<=1:
    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')

However, the code fails, and gives the folowing message during the iterations:

*** Warning: Found no facets matching domain for boundary condition.

I hope the phase fiel length scale is correct, because the distance between nodes is 1. The mesh is in accordance with a standard for the test. If the code works after this question I will make a mesh more fine around where the phenomena is actually occurring.

Sven · February 11, 2022, 1:30pm

The following is likely the issue.

Your updated mesh doesn’t contain any points that meet this criteria. Consider generating pvd files of your edges to check whether the boundaries are correctly being applied.

Caio_Cesar · February 15, 2022, 12:45pm

I changed the mesh of the standard to make it easier. Now it is a triangle, like this:

But I can’t make the fracture start at the tip of the triangle, as it should. The boundary conditions are the same, what changes is the pre-crack that I made.

How can I solve it?

Sven · February 15, 2022, 1:27pm

Can you elaborate? How are you assigning the crack boundary condition? Please provide more information such as your Crack(x). Have you tried visualising the starter crack conditions to confirm the boundary condition is being applied?

Caio_Cesar · February 15, 2022, 10:22pm

To make it easier, this is my code:

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

# Mesh
vertices=[Point(8.5,0),Point(5.1,-0.9),Point(0,-2),
Point(0,-140),Point(22,-140),Point(22,140),Point(0,140),Point(0,2),Point(5.1,0.9)]
domain=Polygon(vertices)
mesh=generate_mesh(domain, 200)

# 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 = 2*mesh.hmin()
E=Constant(210e3)
nu=Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)


# 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],140) && on_boundary")
bot = CompiledSubDomain("near(x[1],-140) && on_boundary")
def Crack(x):
    return abs(x[1])<1e-3 and x[0]<=8.5
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 ("./Results/phi.pvd")
fname = open('ForcevsDisp.txt', 'w')

# Staggered scheme
while t<=1:
    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')

The phase field length scale is calculated from the mesh.

I have tried visualising the starter crack conditions to confirm the boundary condition is being applied using Paraview, but fails. ´

Sven · February 15, 2022, 11:44pm

You have intentionally created a mesh that excludes the starter crack, therefore there are no points that match your Crack(x) still. The following visually and numerically exemplifies this.

class CrackClass(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1])<1e-3 and x[0]<=8.5
cracks = MeshFunction("size_t", mesh, mesh.topology().dim())
cracks.set_all(0)
crackClass = CrackClass()
crackClass.mark(cracks, 1)
crackFile = File("cracks.pvd")
crackFile << cracks
dz = Measure("dx", mesh, subdomain_data=cracks)
print("Total surface volume is {0}".format(assemble(Constant(1)*dz(1))))

I’m not sure what standard you want to replicate (something like ASTM E1820), but you should consider the expected crack profile when deciding how to apply the starter crack with the geometry you’ve chosen. You are no longer using the author’s mesh so I’d advise you to rework the Crack(x) to match your test’s standard.

arianaqyp · February 16, 2022, 2:40am

Hey Sven, I am using the same code here and I have a question, if I may ask… If I were to change the loading to shear, do I just change these line of code:

bctop = DirichletBC(W.sub(1), load, top)
…
Traction = dot(sigma(unew),n)
fx = Traction[0]*ds(1)
fname.write(str(assemble(fx)) + “\n”)

I am trying to get the reaction force from a shear loading, but the graph I got doesn’t seem right…

Appreciate any input!

This is what I got for tension loading:

And this is for shear loading from changing those lines of code:

I am expecting it to be the same range of reaction force, but I can’t seem to figure out how…

Caio_Cesar · February 16, 2022, 12:16pm

The correct mesh is that:

vertices=[Point(8.5,0),Point(5.1,-0.9),Point(2.9,-0.9),Point(1.5,-3),
Point(0,-2),Point(0,-140),Point(22,-140),Point(22,140),
Point(0,140),Point(0,2),Point(1.5,3),Point(2.9,0.9),Point(5.1,0.9)]

but I can’t rework the Crack(x) to match the test’s standard. It is expected to start at the “tip” (x=8.5)

Sven · February 16, 2022, 1:19pm

If you provide the test standard that could help us provide specific advice. In lieu of that I can suggest the following.

Consider if crack tip plasticity is relevant to your standard. Potentially refine the mesh along the x-axis and create a elastic-plastic zone around x=8.5
As per ASTM E1820, place a fatigue pre-crack along the direction of the crack starting from the tip at x=8.5

Even if you :

The phase field model will diffuse the crack based on your critical crack length l.

Sven · February 16, 2022, 1:25pm

@arianaqyp I see you have already opened a post related to this question : shear loading horizontal loading on a unit mesh. I recommend that you close that post or move this content there to avoid duplication.

arianaqyp · February 16, 2022, 1:37pm

Hey @Sven ! Yes I’ll move the content here.

My question: How do I modify the code for shear loading?

Specific question:

I have a working code for tension loading (vertical loading) (same mesh and code as the authors) and I can compute the reaction force by using Traction = sigma*n where n is the normal to the top boundary surface. And fy = Traction[1]*ds(1) where ds(1) is top boundary.

Questions:

How do I do this for shear loading where the loading is in the x direction on the top boundary.
I did:

if method == ‘tension’:
bctop = DirichletBC(W.sub(1), load, top)
elif method == ‘shear’:
bctop = DirichletBC(W.sub(0), load, top)

Is this right?

What is the proper way to calculate shear reaction force?
I tried:

fx = Traction[0]*ds(1)
and
tang = as_vector([n[1], -n[0]])
fx = ufl.inner(sigma(unew) * n, tang)*ds(1)

but I know for sure it is wrong because the Reaction force do not match the one I had for tension. I’ve seen many papers, and they all have the same range of reaction force but plots… Also, when I did the visualization on paraview, the crack propagated the same way as tension loading and did not curve downwards, which is what I was expecting…

Do I have to change the traction calculation and define a tangential vector? If yes, how?

It might be the way the boundary condition is defined at the top surface, but I can’t seem to figure it out…

I appreciate any help I can get, I have been looking at this for days, but can’t seem to find a solution…

Sven · February 16, 2022, 4:48pm

arianaqyp,

It would help if you included the code from the past post for reference.

Your previous comment

is to be expected. You are applying a Displacement driven loading condition to the top surface, not force/traction driven. You should not always expect similar forces. Additionally when in shear you are loading in the direction of the crack whereas in tension you are loading perpendicular. You are getting shear loading by forcing the top surface to translate to the right through a DirichletBC.

2/3. Since you are using the unit square mesh your loading is along a coordinate direction. If you want the shear force on the top surface acting in the y direction consider this post. I would think any of the following would work.

# Referenced Post // 
Tn = assemble(inner(sigma(unew)*n,n)*ds(1)) # Normal
Tt = assemble((sigma(unew)*n - Tn*n)[0]*ds(1))  # Tangential, x[1] is ~0

# Your approach //
fx = assemble(dot(sigma(unew),n)[0]*ds(1)) # Tangential

# Direct //
fx = assemble(Constant(0.5)*(sigma(unew)[0,1]+sigma(unew)[1,0])*ds(1)) # Tangential

Consider this with the following MWE:

from fenics import *
mesh = UnitSquareMesh(20, 20)
n = FacetNormal(mesh)
Q = TensorFunctionSpace(mesh, 'CG', 1)
exp = Expression(((("1"),("2")),(("2"),("4"))), degree = 2)
sigma = project(exp,Q)
tol = 1e-14

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1, tol) and on_boundary

edges = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
edges.set_all(0)
top = Top()
top.mark(edges, 1)
ds = Measure('ds', domain = mesh, subdomain_data = edges)

Tn = assemble(inner(sigma*n,n)*ds(1))
Tt = assemble((sigma*n-Tn*n)[0]*ds(1))
S21 = assemble(Constant(0.5)*(sigma[0,1]+sigma[1,0])*ds(1))
T3 = assemble(dot(sigma,n)[0]*ds(1))
print(Tt,S21,T3)

I hope that helps.

Topic		Replies	Views
Code for bend test using phase field fracture variational formulation	16	835	April 29, 2022
Phase field modelling for a notched plate with hole variational formulation	0	373	April 20, 2022
Doubt in brittle fracture General	4	162	March 22, 2024
Phase field simulation - damage at the boundary	3	808	April 18, 2022
ArityMismatch error General	9	552	July 19, 2023

Phase Field Fracture

Questions:

Related topics