Hello all,
I am trying to model the thermoviscoelastic material behaviour and I have a FEniCS model for viscoelastic model. I am sharing my model below. With using the prony coefficients and corresponding relaxation time and WLF parameters, how can I change my model to thermoviscoelastic model? Thank you in advance.
from dolfin import *
from ufl import replace
import numpy as np
import matplotlib.pyplot as plt
L, H = 0.1, 0.2
mesh = RectangleMesh(Point(0., 0.), Point(L, H), 5, 10)
E0 = Constant(70e3)
E1 = Constant(20e3)
E2 = Constant(20e3)
E3 = Constant(20e3) # Third Maxwell element
E4 = Constant(20e3) # Fourth Maxwell element
eta1 = Constant(1e3)
eta2 = Constant(2e2)
eta3 = Constant(2e2) # Third Maxwell element
eta4 = Constant(2e2) # Fourth Maxwell element
nu = Constant(0.)
dt = Constant(0.) # time increment
sigc = 100. # imposed creep stress
epsr = 1e-3 # imposed relaxation strain
def left(x, on_boundary):
return near(x[0], 0.) and on_boundary
def bottom(x, on_boundary):
return near(x[1], 0.) and on_boundary
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], H) and on_boundary
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Top().mark(facets, 1)
ds = Measure("ds", subdomain_data=facets)
Ve = VectorElement("CG", mesh.ufl_cell(), 1)
Qe = TensorElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([Ve, Qe, Qe, Qe, Qe])) # Four Qe for four Maxwell elements
w = Function(W, name="Variables at current step")
(u, epsv1, epsv2, epsv3, epsv4) = split(w) # Split into four viscoelastic strains
w_old = Function(W, name="Variables at previous step")
(u_old, epsv1_old, epsv2_old, epsv3_old, epsv4_old) = split(w_old)
w_ = TestFunction(W)
(u_, epsv1_, epsv2_, epsv3_, epsv4_) = split(w_)
dw = TrialFunction(W)
def eps(u):
return sym(grad(u))
def dotC(e):
return nu/(1+nu)/(1-nu)*tr(e)*Identity(2) + 1/(1+nu)*e
def sigma(u, epsv1, epsv2, epsv3, epsv4):
return E0*dotC(eps(u)) + E1*dotC(eps(u) - epsv1) + E2*dotC(eps(u) - epsv2) + E3*dotC(eps(u) - epsv3) + E4*dotC(eps(u) - epsv4)
def strain_energy(u, epsv1, epsv2, epsv3, epsv4):
e = eps(u)
return 0.5*(E0*inner(e,dotC(e)) + E1*inner(e-epsv1, dotC(e-epsv1)) + E2*inner(e-epsv2, dotC(e-epsv2)) + E3*inner(e-epsv3, dotC(e-epsv3)) + E4*inner(e-epsv4, dotC(e-epsv4)))
def dissipation_potential(depsv1, depsv2, depsv3, depsv4):
return 0.5*eta1*inner(depsv1, depsv1) + 0.5*eta2*inner(depsv2, depsv2) + 0.5*eta3*inner(depsv3, depsv3) + 0.5*eta4*inner(depsv4, depsv4)
Traction = Constant(0.)
incremental_potential = strain_energy(u, epsv1, epsv2, epsv3, epsv4)*dx \
+ dt*dissipation_potential((epsv1-epsv1_old)/dt, (epsv2-epsv2_old)/dt, (epsv3-epsv3_old)/dt, (epsv4-epsv4_old)/dt)*dx \
- Traction*u[1]*ds(1)
F = derivative(incremental_potential, w, w_)
J = derivative(F, w, dw)
dimp = Constant(H*epsr) # imposed vertical displacement instead
bcs = [DirichletBC(W.sub(0).sub(0), Constant(0), left),
DirichletBC(W.sub(0).sub(1), Constant(0), bottom),
DirichletBC(W.sub(0).sub(1), dimp, facets, 1)]
def viscoelastic_test(case, Nsteps=50, Tend=1.):
# Solution fields are initialized to zero
w.interpolate(Constant((0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.))) # 18 zeros for 18 DoFs
# Define proper loading and BCs
if case in ["creep", "recovery"]: # imposed traction on top
Traction.assign(Constant(sigc))
bc = bcs[:2] # remove the last boundary conditions from bcs
t0 = Tend/2. # traction goes to zero at t0 for recovery test
elif case == "relaxation":
Traction.assign(Constant(0.)) # no traction on top
bc = bcs
# Time-stepping loop
time = np.linspace(0, Tend, Nsteps+1)
Sigyy = np.zeros((Nsteps+1, ))
Epsyy = np.zeros((Nsteps+1, 5))
for (i, dti) in enumerate(np.diff(time)):
if i>0 and i % (Nsteps/5) == 0:
print("Increment {}/{}".format(i, Nsteps))
dt.assign(dti)
if case == "recovery" and time[i+1] > t0:
Traction.assign(Constant(0.))
w_old.assign(w)
solve(F == 0, w, bc, J=J)
# get average stress/strain
Sigyy[i+1] = assemble(sigma(u, epsv1, epsv2, epsv3, epsv4)[1, 1]*dx)/L/H
Epsyy[i+1, 0] = assemble(eps(u)[1, 1]*dx)/L/H
Epsyy[i+1, 1] = assemble(epsv1[1, 1]*dx)/L/H
Epsyy[i+1, 2] = assemble(epsv2[1, 1]*dx)/L/H
Epsyy[i+1, 3] = assemble(epsv3[1, 1]*dx)/L/H
Epsyy[i+1, 4] = assemble(epsv4[1, 1]*dx)/L/H
# Define analytical solutions
if case == "creep":
if float(E0) == 0.:
eps_an = sigc*(1./float(E1)+time/float(eta1))
else:
Estar1 = float(E0*E1/(E0+E1))
tau1 = float(eta1)/Estar1
Estar2 = float(E0*E2/(E0+E2))
tau2 = float(eta2)/Estar2
Estar3 = float(E0*E3/(E0+E3))
tau3 = float(eta3)/Estar3
Estar4 = float(E0*E4/(E0+E4))
tau4 = float(eta4)/Estar4
eps_an = sigc/float(E0)*(1 - float(Estar1/E0)*np.exp(-time/tau1) - float(Estar2/E0)*np.exp(-time/tau2)
- float(Estar3/E0)*np.exp(-time/tau3) - float(Estar4/E0)*np.exp(-time/tau4))
sig_an = sigc + 0*time
elif case == "relaxation":
if float(E1) == 0.:
sig_an = epsr*float(E0) + 0*time
else:
tau1 = float(eta1/E1)
tau2 = float(eta2/E2)
tau3 = float(eta3/E3)
tau4 = float(eta4/E4)
sig_an = epsr*(float(E0) + float(E1)*np.exp(-time/tau1) + float(E2)*np.exp(-time/tau2)
+ float(E3)*np.exp(-time/tau3) + float(E4)*np.exp(-time/tau4))
eps_an = epsr + 0*time
elif case == "recovery":
Estar1 = float(E0*E1/(E0+E1))
tau1 = float(eta1)/Estar1
Estar2 = float(E0*E2/(E0+E2))
tau2 = float(eta2)/Estar2
Estar3 = float(E0*E3/(E0+E3))
tau3 = float(eta3)/Estar3
Estar4 = float(E0*E4/(E0+E4))
tau4 = float(eta4)/Estar4
eps_an = sigc/float(E0)*(1 - float(E1/(E0+E1))*np.exp(-time/tau1) - float(E2/(E0+E2))*np.exp(-time/tau2)
- float(E3/(E0+E3))*np.exp(-time/tau3) - float(E4/(E0+E4))*np.exp(-time/tau4))
sig_an = sigc + 0*time
time2 = time[time > t0]
sig_an[time > t0] = 0.
eps_an[time > t0] = sigc*(float(E1/E0/(E0+E1))*(np.exp(-(time2-t0)/tau1) - np.exp(-time2/tau1))
+ float(E2/E0/(E0+E2))*(np.exp(-(time2-t0)/tau2) - np.exp(-time2/tau2))
+ float(E3/E0/(E0+E3))*(np.exp(-(time2-t0)/tau3) - np.exp(-time2/tau3))
+ float(E4/E0/(E0+E4))*(np.exp(-(time2-t0)/tau4) - np.exp(-time2/tau4)))
# Plot strain evolution
plt.figure()
plt.plot(time, 100*eps_an, label="analytical solution")
plt.plot(time, 100*Epsyy[:, 0], '.', label="FE solution")
plt.plot(time, 100*Epsyy[:, 1], '--', linewidth=1, label="viscous strain 1")
plt.plot(time, 100*Epsyy[:, 2], '--', linewidth=1, label="viscous strain 2")
plt.plot(time, 100*Epsyy[:, 3], '--', linewidth=1, label="viscous strain 3")
plt.plot(time, 100*Epsyy[:, 4], '--', linewidth=1, label="viscous strain 4")
plt.ylim(0, 1.2*max(Epsyy[:, 0])*100)
plt.xlabel("Time")
plt.ylabel("Vertical strain [\%]")
plt.title(case.capitalize() + " test")
plt.legend()
plt.savefig("Strain")
plt.show()
# Plot stress evolution
plt.figure()
plt.plot(time, sig_an, label="analytical solution")
plt.plot(time, Sigyy, '.', label="FE solution")
plt.ylim(0, 1.2*max(Sigyy))
plt.ylabel("Vertical stress")
plt.xlabel("Time")
plt.title(case.capitalize() + " test")
plt.legend()
plt.savefig("Stress")
plt.show()
viscoelastic_test("relaxation")