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:
-
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) 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… -
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.7692e3Constituive functions
def epsilon(u):
return sym(grad(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)Boundary conditions
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)Define shear/tension loading
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_du = ((1.0-pold)**2)inner(grad(v),sigma(u))dx
E_phi = (Gclinner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))
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-3print(‘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 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 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 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) 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)