# Variational formulation in diffusion problem

I am trying to implement a 3-field poroviscoelastic problem. My question is about the diffusion part:

\dot{C}+\text{div}(J_i) =0 ,

with:

• \varepsilon_{kk}=\left(C-C_0\right)\Omega (coupling between volumetric strain and solvent concentration)
• J_i=-\left(\frac{\kappa}{\eta\Omega^2}\right)\frac{\partial \mu_s}{\partial x_i} (Ficks law, \mu: chemical potential)

Then:

is it well defined?:

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)
vw = TestFunction(W)
(vu, vepsv, vmu_s) = split(vw)
dw = TrialFunction(W)
(du, depsv, dmu_s) = split(dw)

chem = tr(eps(du)-eps(u_old))/dt*vmu_s*dx +\



I ask because I can’t get the correct results. I appreciate any help.

Thank you!

It’s difficult to give advice when you’ve omitted a lot of information. At a glance the equation looks fine, but your FE formulation does not have a sufficient number of equations (1) for unknowns (3).

If I were working on your problem I’d carefully build up the problem until you arrive at your full model.

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):

1. \text{div}[\sigma_{ij}]+f=0
2. \dot{\varepsilon}^v_{ij}-\frac{1}{\eta_v}\sigma^\beta_{ij}=0 (viscous strain of the maxwell element)
3. \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):

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))
right = AutoSubDomain(lambda x: near(x, L))
bot = AutoSubDomain(lambda x: near(x, 0))
top = AutoSubDomain(lambda x: near(x, 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]
## 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)
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!