Thank you for answering, I will leave more details below.
I am solving the problem for 3 fields: displacement, viscous deformation, and chemical potential.
The equilibrium equations are the following (linear viscoelasticity, I already checked the first two with the tutorial):
- \text{div}[\sigma_{ij}]+f=0
-
\dot{\varepsilon}^v_{ij}-\frac{1}{\eta_v}\sigma^\beta_{ij}=0 (viscous strain of the maxwell element)
-
\dot{\varepsilon}_{kk}-\frac{\kappa}{\eta\Omega} \ \text{div}[\text{grad}(\mu_s)]=0 (as defined above, \eta: viscosity of the polymer network)
and the stress is the sum of the contribution of the linear standard solid model and the diffusion term:
\sigma_{ij}=2G^\alpha \left[ \varepsilon_{ij} + \frac{\nu^\alpha}{1-2\nu^\alpha}\varepsilon_{kk} \right] + 2G^\beta \left[ (\varepsilon_{ij} - \varepsilon^v_{ij})+ \frac{\nu^\beta}{1-2 \nu^\beta}(\varepsilon_{kk}- \varepsilon^v_{kk}) \right] + \frac{\mu_e - \mu_s}{\Omega}\delta_{ij}
where:
-
\mu_s: chemical potential of the solvent
-
\mu_e: chemical potential of the enviroment
My problem is that when simulating creep, as the strain tends to an equilibrium, the chemical potential gradient does not do so, it even remains constant.
The following is an example based on the linear viscoelasticity tutorial. It is a rectangular body with restricted displacement at the bottom and to the left, subjected to an external chemical potential at the top and right. At time zero, the solvents are in equilibrium \mu_s=\mu_e=8.430953e6\ [\text{J} \text{mol}^{-1}].
This is my complete code:
from dolfin import *
import numpy as np
## Mesh
L, H = 0.00125, 0.005
mesh = RectangleMesh(Point(0., 0.), Point(L, H), 5, 10)
## Material parameters
E0 = Constant(70e3)
E1 = Constant(20e3)
eta1 = Constant(1e3)
nu = Constant(0)
#
Ga = E0/2/(1+nu)
Gb = E1/2/(1+nu)
NUa= nu
NUb= nu
# solvent (water)
Omega= 18000
kappa= 6e-12
eta = 1e-10
## Functions
def eps(u):
return sym(grad(u))
def Stress(u,epsv,mu):
# Total stress
I = Identity(2)
# Elastic
EPSij = eps(u)
EPSkk = tr(EPSij)*I
SIGa = Constant(2.*Ga)*(EPSij + Constant(NUa/(1.-2.*NUa))*EPSkk)
# Viscous
EPSvij = epsv
EPSvkk = tr(EPSvij)*I
SIGb = Constant(2.*Gb)*(EPSij-EPSvij + Constant(NUb/(1.-2.*NUb))*(EPSkk-EPSvkk))
# Chemical
SIGc = (mu_e - mu)/Omega * I
return SIGa + SIGb + SIGc
def StressB(u,epsv):
# Stress in the maxwell element
I = Identity(2)
# Elastic strain
EPSij = eps(u)
EPSkk = tr(EPSij)*I
# Viscous
EPSvij = epsv
EPSvkk = tr(EPSvij)*I
SIGb = Constant(2.*Gb)*(EPSij-EPSvij + Constant(NUb/(1.-2.*NUb))*(EPSkk-EPSvkk))
return SIGb
## Functional spaces
Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Qe = TensorElement("DG", mesh.ufl_cell(), 0)
MUe = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([Ve, Qe, MUe]))
w = Function(W, name="Variables at current step")
## Initial values
w.interpolate(Constant((0.,0.,0.,0.,0.,0.,8.430953e6)))
(u, epsv, mu_s) = split(w)
w_ = TestFunction(W)
(u_, epsv_, mu_s_) = split(w_)
dw = TrialFunction(W)
(du, depsv, dmu_s) = split(dw)
w_old = Function(W, name="Variables at previous step")
(u_old, epsv_old, mu_s_old) = split(w_old)
## Dirichlet bc
facets = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
left = AutoSubDomain(lambda x: near(x[0], 0))
right = AutoSubDomain(lambda x: near(x[0], L))
bot = AutoSubDomain(lambda x: near(x[1], 0))
top = AutoSubDomain(lambda x: near(x[1], H))
facets.set_all(0)
left.mark(facets, 1)
right.mark(facets, 2)
bot.mark(facets, 3)
top.mark(facets, 4)
dx = Measure('dx', mesh)
ds = Measure('ds', mesh, subdomain_data=facets)
mu_e = Constant(8.430953e6) # environment
bcs = [DirichletBC(W.sub(0).sub(0), Constant(0), facets, 1),
DirichletBC(W.sub(0).sub(1), Constant(0), facets, 3),
DirichletBC(W.sub(2), mu_e, facets, 4),
DirichletBC(W.sub(2), mu_e, facets, 2)]
## Variational problem
dt = 1e-2
sigc = 50
T = Constant((0.0,0.0))
## mechanical part [J]
mech = inner(Stress(du, depsv, dmu_s),eps(u_))*dx - dot(T,u_)*ds(1)
## Viscoelastic part [1/s]
visc = inner((depsv-epsv_old)/dt - 1/eta1*StressB(du,depsv), epsv_)*dx
## Chemical part [J mm3/mol/s]
chem = tr(grad(du)-grad(u_old))/dt*mu_s_*dx + (kappa/eta/Omega)*dot(grad(dmu_s),grad(mu_s_))*dx
## Total
form = mech + visc + chem
# Time-stepping loop
tTotal = 0.5
Nsteps = int(tTotal/dt+1)
Time = np.linspace(0, tTotal, Nsteps)
Sigyy = np.zeros((Nsteps, 1))
Epsyy = np.zeros((Nsteps, 2))
displacements = np.zeros((Nsteps,1))
vtk_chemical = File("chemical/chemical.pvd")
vtk_disp = File("displacement/disp.pvd")
for i, d in enumerate(Time):
# creep
if i == 1:
T.assign(Constant((0.0,sigc)))
w_old.assign(w)
solve(lhs(form) == rhs(form), w, bcs)
u , epsv, mu_s = w.split()
u_old , epsv_old, mu_s_old = w_old.split()
# get average stress/strain
Sigyy[i] = assemble(Stress(u, epsv, mu_s)[1, 1]*dx)/L/H
Epsyy[i, 0] = assemble(eps(u)[1, 1]*dx)/L/H
Epsyy[i, 1] = assemble(epsv[1, 1]*dx)/L/H
displacements[i] = u(L/2,H)[1]
vtk_chemical << (mu_s,d)
vtk_disp << (u,d)
####
## analytical (linear viscoelasticity)
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
I used paraview to see the evolution of the chemical potential but it remains constant.
Thanks for any help!
