Hi all I am trying to interpret my results for some phase field simulations, but I am confused about the results.
Goal: An input stochastic field, in this case fracture toughness - which outputs a corresponding damage field (solved using phasefield)
Most of the results are consistent, but a few were not (and are probably incorrect?)
Here are my configurations for my phase field simulation:
h == mesh size
l_phasefield = 0.03
l_phasefield/h = 3
boundary condition = fixed at the bottom, tension pull on the top edge
This is an overlayed plot of the damage for 70 samples, and you can see that some simulated damage are at the top and bottom edge, which is probably not right?
An example of the output
Code is at the bottom.
I appreciate any advice!
Best,
Ariana
“”"
Code
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)
realizationpath = os.path.join(outputpath, 'realization_' + str(i))
# stochastic
# Introduce manually the material parameters
# write to xdmf to view in paraview
Gc = Function(WW)
xs = WW.tabulate_dof_coordinates()
xs_x = xs[:, 0]
xs_y = xs[:, 1]
# interpolate using rbf interpolator
# rbf_i = Rbf(rfcentroid_x, rfcentroid_y, z_samples[i], function='gaussian')
rbf_i = Rbf(X, Y, z_samples[i], function='gaussian')
z_centroids = rbf_i(xs_x, xs_y)
Gc.vector()[:] = z_centroids
print('z_centroids.shape: ', z_centroids.shape)
# save Gc to pvd
Gc_file = File(os.path.join(realizationpath, "Gc.pvd"))
Gc_file << Gc
# plot and save Gc
# plot z_sample i
plt.figure()
# c = plt.contourf(rfcentroid_x, rfcentroid_y, z_samples[i].reshape(N, M), cmap=cm.RdBu)
c = plt.contourf(X, Y, z_samples[i].reshape(N_um, M_um), cmap=cm.RdBu)
plt.colorbar(c)
# plot with grid
# plt.grid()
# save figure
plt.savefig(realizationpath + '/z_sample_' + str(i) + '.png')
# plot Gc on mesh
plt.figure(3)
fig = plt.figure(figsize=(7,7))
ax = fig.gca()
plot(mesh)
gc=plot(Gc, title='Gc', cmap=cm.RdBu)
# color bar with 2 values
cbar = plt.colorbar(gc)
# save plot
plt.savefig(os.path.join(realizationpath, 'Gc.png'))
lmbda = 121.1538e3
mu = 80.7692e3
# 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):
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)
# tension loading
if method == 'tension':
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot, bctop]
elif method == 'shear':
bctopy= DirichletBC(W.sub(1), Constant(0.0), bot)
bctop = DirichletBC(W.sub(0), load, top)
bc_u = [bcbot, bctop, bctopy]
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]])
tang = as_vector([n[0], -n[1]])
# 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 = (Gc*l*inner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))\
*inner(p,q)-2.0*H(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(realizationpath, 'ResultsDir/phi.pvd'))
# save z centroids
np.savetxt(os.path.join(realizationpath, 'ResultsDir/z_centroids.txt'), z_centroids)
fname = open(os.path.join(realizationpath, '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.001
#! 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)
Tn = assemble(inner(sigma(unew)*n,n)*ds(1))
# 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(Tn) + "\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)
Tn = assemble(inner(sigma(unew)*n,n)*ds(1))
Tt = assemble((sigma(unew)*n - Tn*n)[0]*ds(1)) # Tangential, x[1] is ~0
# fx = ufl.inner(sigma(unew) * n, tang)*ds(1)
# print('I am writing to file')
fname.write(str(t*u_r) + "\t")
fname.write(str(Tt) + "\n")
fname.close()
“”"