from dolfin import *
from ufl_legacy 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)
eta1 = Constant(1e3)
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]))
w = Function(W, name="Variables at current step")
(u, epsv) = split(w)
w_old = Function(W, name="Variables at previous step")
(u_old, epsv_old) = split(w_old)
w_ = TestFunction(W)
(u_, epsv_) = 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, epsv):
return E0*dotC(eps(u)) + E1*dotC(eps(u) - epsv)
def strain_energy(u, epsv):
e = eps(u)
return 0.5*(E0*inner(e,dotC(e)) + E1*inner(e-epsv, dotC(e-epsv)))
def dissipation_potential(depsv):
return 0.5*eta1*inner(depsv, depsv)
Traction = Constant(0.)
incremental_potential = strain_energy(u, epsv)*dx \
+ dt*dissipation_potential((epsv-epsv_old)/dt)*dx \
- Traction*u[1]*ds(1)
F = derivative(incremental_potential, w, w_)
form = replace(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.,)*6))
# 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, 2))
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(lhs(form) == rhs(form), w, bc)
# get average stress/strain
Sigyy[i+1] = assemble(sigma(u, epsv)[1, 1]*dx)/L/H
Epsyy[i+1, 0] = assemble(eps(u)[1, 1]*dx)/L/H
Epsyy[i+1, 1] = assemble(epsv[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:
Estar = float(E0*E1/(E0+E1))
tau = float(eta1)/Estar
eps_an = sigc/float(E0)*(1-float(Estar/E0)*np.exp(-time/tau))
sig_an = sigc + 0*time
elif case == "relaxation":
if float(E1) == 0.:
sig_an = epsr*float(E0) + 0*time
else:
tau = float(eta1/E1)
sig_an = epsr*(float(E0) + float(E1)*np.exp(-time/tau))
eps_an = epsr + 0*time
elif case == "recovery":
Estar = float(E0*E1/(E0+E1))
tau = float(eta1)/Estar
eps_an = sigc/float(E0)*(1-float(E1/(E0+E1))*np.exp(-time/tau))
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)/tau)
- np.exp(-time2/tau))
# 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")
plt.ylim(0, 1.2*max(Epsyy[:, 0])*100)
plt.xlabel("Time")
plt.ylabel("Vertical strain [\%]")
plt.title(case.capitalize() + " test")
plt.legend()
plt.show()
plt.savefig("strain_evoluation")
# 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.show()
plt.savefig("stress_evoluation")
# ## relaxation test --------------------------------------------------------------------------------------------------------------------------------------
viscoelastic_test("relaxation")
# ### Creep test ----------------------------------------------------------------------------------------------------------------------------
viscoelastic_test("creep")
# ### Recovery test ---------------------------------------------------------------------------------------------------------------------------------
viscoelastic_test("recovery")
# ## Relaxation and creep tests for a Maxwell model
#
# We give here the solutions for a Maxwell model which is obtained from the degenerate case $E_0=0$. We recover that the strain evolves linearly with time for the creep test.
E0.assign(Constant(0.))
viscoelastic_test("relaxation")
viscoelastic_test("creep")
# ## Relaxation and creep tests for a Kelvin-Voigt model
#
E0.assign(Constant(70e3))
E1.assign(Constant(1e10))
viscoelastic_test("relaxation")
viscoelastic_test("creep")
For the viscoelastic material solution above the solution is calculated with using the energy equations but I want to use the weak form of balance of momentum and calculate the same solution in this way. The equation which I want to use is that :
What changes should I do to implement this to my code ?