Phase Field Fracture

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.


  1. 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?

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

fx = Traction[0]*ds(1)
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…

  1. 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…