Hi I have a rookie question for shear loading.

This code is adapted from :
Hirshikesh, S. Natarajan, R. K. Annabattula, E. Martinez-Paneda.
Phase field modelling of crack propagation in functionally graded materials.
Composites Part B: Engineering 169, pp. 239-248 (2019)
doi: 10.1016/j.compositesb.2019.04.003

Emilio Martinez-Paneda (mail@empaneda.com)
University of Cambridge

So I have a working code for tension loading (vertical loading) 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.

Mesh: unit square mesh

## Question:

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?

2. What is the proper way to calculate shear reaction force?
I tried:
fx = Traction[0]*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…

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

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

Here is the code I have:

#### 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 = l
lmbda = 121.1538e3
mu = 80.7692e3

#### Constituive functions

def epsilon(u):
def sigma(u):
return 2.0muepsilon(u)+lmbdatr(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)

#### y = 1

top = CompiledSubDomain(“near(x[1], 1) && on_boundary”)

#### y = 0

bot = CompiledSubDomain(“near(x[1], 0) && on_boundary”)

#### Crack near y=0.5 and x<=0.3

def Crack(x):
return ((0.5-1e-03)< x[1]<(0.5+1e-03)) and (x[0] <= 0.5)

load = Expression(“t”, t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)

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

elif method == ‘shear’:
bctop = DirichletBC(W.sub(0), 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)

#### define tangent directions on mesh

tang = as_vector([n[1], -n[0]])

#### Variational form

unew, uold = Function(W), Function(W)
pnew, pold, Hold = Function(V), Function(V), Function(V)
E_phi = (Gc
l
inner(p,q)-2.0H(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 # mm
deltaT = 0.1
tol = 1e-3

print(‘Creating Results file…’)

conc_f = File(os.path.join(outputpath, ‘ResultsDir/phi.pvd’))
fname = open(os.path.join(outputpath, ‘ForcevsDisp.txt’), ‘w’)

print(‘Starting the iterative procedure’)

#### Staggered scheme

start = timeit.default_timer()

if method == ‘tension’:
while t<=1.0:
t += deltaT
if t >=0.7:
deltaT = 0.01

``````        #### fast test
# deltaT = 0.01

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

fy = ufl.inner(sigma(unew) * n, n)*ds(1)
# R = dolfinx.fem.assemble_scalar(t)

# Traction = dot(sigma(unew),n)
# fy = Traction[1]*ds(1)
# fy = sigma(unew)[0, 1]*ds(1)
# print('I am writing to file')
fname.write(str(t*u_r) + "\t")
fname.write(str(assemble(fy)) + "\n")

fname.close()
``````

elif method == ‘shear’:
while t<=1.0:
t += deltaT
if t >=0.7:
deltaT = 0.01

``````        #! fast test
# deltaT = 0.01

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)
fx = Traction[0]*ds(1)
# print('I am writing to file')
fname.write(str(t*u_r) + "\t")
fname.write(str(assemble(fx)) + "\n")

fname.close()
``````

end = timeit.default_timer()
print(‘Start time’, start)
print(‘End time’, end)
sim_t = end - start
sim_t = round(sim_t/60, 4)

print(‘Finished the iterative procedure’)
print(‘Simulation time:’, sim_t)

Please encapsulate your code with 3x `
ti make sure that it gets formatted properly.

Oh sorry I was trying to figure out how to do that too. Fixed it!

I really appreciate any input/leads… I’ve been staring at this for too long…