Hello everyone,
I have a code for a three point bend test using phase fiel method. This code returns a Force vs displacement graph.
from fenics import *
from dolfin import *
from mshr import *
parameters['form_compiler']['cpp_optimize'] = True
parameters['form_compiler']['optimize'] = True
mesh = Mesh('mesh.xml')
boundaries = MeshFunction('size_t', mesh,'.boundaries.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 = 0.54
l = 0.03
lmbda =12e3
mu = 8e3
# 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):
str_ele = 0.5*(grad(u) + grad(u).T)
IC = tr(str_ele)
ICC = tr(str_ele * str_ele)
return (0.5*lmbda*IC**2) + mu*ICC
# Boundary conditions
class top(SubDomain):
def inside(self, x, on_boundary):
#return x[0] - 3.8 > 1e-4 and x[0] - 4.2 < 1e-4 and abs(x[1] - 2) < 1e-4 and on_boundary
return near(x[0],4,0.2) and near(x[1],2,0.3)
class bot1(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],0,0.2) and near(x[1],0,0.2)
class bot2(SubDomain):
def inside(self, x,on_boundary):
return near(x[0],8,0.2) and near(x[1],0,0.2)
top = top()
bot1 = bot1()
bot2 = bot2()
load = Expression("-t", t = 0.00001, degree=1)
bcbot1x= DirichletBC(W.sub(0), 0, bot1, method="pointwise")
bcbot1y= DirichletBC(W.sub(1), 0, bot1, method="pointwise")
bcbot2= DirichletBC(W.sub(1),0, bot2, method="pointwise")
bctop = DirichletBC(W.sub(1), load, top, method="pointwise")
bc_u = [bcbot1x, bcbot1y, bcbot2, bctop]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds", domain = mesh, subdomain_data=boundaries)
n = FacetNormal(mesh)
# Variational form
unew, uold = Function(W), Function(W)
pnew, pold, Hold = Function(V), Function(V), Function(V)
phi_conv,u_conv = Function(V),Function(W)
E_du = ((1.0-pold)**2+1e-6)*inner(grad(v),sigma(u))*dx
E_phi = Gc*l*inner(grad(p),grad(q))*dx +\
((Gc/l) + 2.*psi(unew))*inner(p,q)*dx-\
2.*psi(unew)*q*dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi), rhs(E_phi), pnew)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)
# Initialization of the iterative procedure and output requests
min_load = 0.01
max_load = 0.12
ut = 1
t = 1e-2
deltaT = 1e-2
Force = Function(W)
conc_f = File ("./Bend_Test/phi.pvd")
fname = open('ForcevsDisp_Bend_Test.txt', 'w')
# Staggered scheme
while t<=max_load:
load.t=t*ut
if 0.04 <= t <= 0.041:
deltaT = 1e-4
if t >0.041:
deltaT = 1e-2
uold.assign(u_conv)
pold.assign(phi_conv)
iter = 1
tol=1e-2
err = 1
maxiter = 100
while err > tol:
solver_disp.solve()
solver_phi.solve()
err_u = errornorm(unew,uold,"L2")/norm(unew,"L2")
err_phi = errornorm(pnew,pold,"L2")/norm(pnew,"L2")
err = max(err_u,err_phi)
uold.assign(unew)
pold.assign(pnew)
iter = iter + 1
#Hold.assign(project(psi(unew), WW))
if err < tol:
print('--------------------------------------------------------------------------')
print ('Iterations:', iter, ', Total time', t)
u_conv.assign(unew)
phi_conv.assign(pnew)
conc_f << pnew
A = assemble(E_du)
Force.vector()[:] = A*unew.vector()
reaction_forces = Force.split(deepcopy=True)[1]
point = (0.0,0.0)
point1 = (8.0,0.0)
react = reaction_forces(point)
react1 = reaction_forces(point1)
a = react+react1
print("Force:", str(round(a*1e-3,3)),'kN')
print("Displacement:",str(round(t,3)),'mm')
print('--------------------------------------------------------------------------')
fname.write(str(round(t,4)) + "\t")
fname.write(str(round(a,4)) + "\n")
t += deltaT
fname.close()
print ('Simulation completed')
The mesh is:
based on Ambati et al and Hirshikesh et al.
My question is: how I can change the displacement to CMOD?