Linear Viscoelasticity

General
1 

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)
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 = F
# ----------------------------------------------------------------------------------------------------------------------------------


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()
    
    # 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()
#  ## 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")

Dear All, I tried to implement and run the tutorial linear viscoelasticity code in my fenics but I get the error below:

fenics@DESKTOP-DCV0R8U:~$ /bin/python3 "/home/fenics/linear_viscoelasticity (4).py"
Form has no parts with arity 2.
Traceback (most recent call last):
  File "/home/fenics/linear_viscoelasticity (4).py", line 152, in <module>
    viscoelastic_test("relaxation")
  File "/home/fenics/linear_viscoelasticity (4).py", line 94, in viscoelastic_test
    solve(lhs(form) == rhs(form), w, bc)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 268, in _solve_varproblem
    problem = LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py", line 62, in __init__
    cpp.fem.LinearVariationalProblem.__init__(self, a, L, u._cpp_object, bcs)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 0).
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

What is my mistake and what should I do to get rid of this error ? Thank you in advance.

2

Hi,
A variational linear problem requires a TestFunction for the right-hand side and a pair of TrialFunction/TestFunction for the left-hand side. In your MWE, the form includes a TestFunction w_ and a function w, making the left-hand side improperly defined. You can either pass this form as an argument to a nonlinear solver, which will then solve the problem in a single iteration, or you can replace w with dw by adding:

 form = replace(F, {w: dw}).

I think you just forgot to copy-paste this line from
Linear viscoelasticity — Numerical tours of continuum mechanics using FEniCS master documentation)
Sincerely,
Paul

3

Thank you for your response, I also tried to use it but it also didnt work then I changed my code to the current version. When I used it it gives this error:

Code:

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)
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()
    
    # 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()
#  ## 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")

Error:


fenics@DESKTOP-DCV0R8U:~$ /bin/python3 "/home/fenics/linear_viscoelasticity (4).py"
Traceback (most recent call last):
  File "/home/fenics/linear_viscoelasticity (4).py", line 61, in <module>
    form = replace(F, {w: dw})
  File "/home/fenics/.local/lib/python3.10/site-packages/ufl/algorithms/replace.py", line 63, in replace
    mapping2 = dict((k, as_ufl(v)) for (k, v) in iteritems(mapping))
  File "/home/fenics/.local/lib/python3.10/site-packages/ufl/algorithms/replace.py", line 63, in <genexpr>
    mapping2 = dict((k, as_ufl(v)) for (k, v) in iteritems(mapping))
  File "/home/fenics/.local/lib/python3.10/site-packages/ufl/constantvalue.py", line 416, in as_ufl
    raise UFLValueError("Invalid type conversion: %s can not be converted"
ufl.log.UFLValueError: Invalid type conversion: v_1 can not be converted to any UFL type.
4

See Announcement: ufl_legacy and legacy dolfin

To fix that:

  1. you must uninstall your local installation of ufl, i.e. the one in /home/fenics/.local/lib/python3.10/site-packages/
  2. you must import ufl_legacy and not ufl
5

Hi, thank you for your answer. I removed the ufl files and tried to run the code with importing ufl_legacy but I took the error below:

 fenics@DESKTOP-DCV0R8U:~$ /bin/python3 "/home/fenics/linear_viscoelasticity (4).py"
Traceback (most recent call last):
  File "/home/fenics/linear_viscoelasticity (4).py", line 61, in <module>
    form = replace(F, {w: dw})
NameError: name 'replace' is not defined
6

Did you change this in

from ufl_legacy import replace

?

7

If you want to backwards compatible I would use

try:
    from ufl import replace
except ModuleNotFoundError:
    from ufl_legacy import replace