Phase Field Fracture


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

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

I hope that helps.