Help with Diffusion problem (Poroviscoelasticity)

General
1

Hi, I am trying to replicate the results of a solvent (water) diffusion simulation inside a polymeric matrix. The problem is implemented using 3 fields: displacements, chemical potential (water), and viscous deformation (linear viscoelasticity).

Considering the following figure (it is a quarter of the body):

image
image1126×635 98.2 KB

I must define a chemical potential \mu_e that simulates the humidity of the environment on the outer faces 1-2-5, and a solvent flux of zero value on faces 3-4-6. BC for the displacements on 3-4-6. And a traction in the top face of 50\ \text{Mpa} to simulate creep.

My doubts are:

  1. Will it be correct to define Dirichlet bc for the solvent chemical potential field to simulate the environment potential?
  2. To impose the zero-valued flux, \partial\mu/\partial x=0, is it sufficient to define nothing?
  3. To sum variational expressions, in this case: linear momentum + mass + viscosity, must all be in the same unit (energy)? I ask because I have seen examples where Joules + 1/sec are summed.

I am not getting the right results, not even close.

I leave my code if you want to check it out, thank you very much!

###### Functional interpolation spaces
element_u    = VectorElement("CG", mesh.ufl_cell(), 2)     # Displacement
element_mu_s = FiniteElement("CG", mesh.ufl_cell(), 1)     # Solvent
element_epsv = TensorElement("DG", mesh.ufl_cell(), 0)     # Viscous strain
mixed_element = MixedElement([element_u, element_mu_s, element_epsv])
V  = FunctionSpace(mesh, mixed_element)
w                  = Function(V)

# Tensor space
ET2   = TensorFunctionSpace(mesh, "DG",0)
# Vector space
ET1   = VectorFunctionSpace(mesh, "CG",2)
# Scalar space
ET0   = FunctionSpace(mesh, "CG", 1)    

##### Initial values
u0  = Function(ET1)
u00 = Expression(("0.0","0.0","0.0"), degree=2)
u0.interpolate(u00)

eps0  = Function(ET2)
eps00 = Expression((("0.0","0.0","0.0"),("0.0","0.0","0.0"),("0.0","0.0","0.0")), degree=0)
eps0.interpolate(eps00)

mus0  = Function(ET0)
mus00 = Expression("8.430953e6", degree=1)
mus0.interpolate(mus00)

assigner = FunctionAssigner(V, [ET1,ET0,ET2])
assigner.assign(w, [u0,mus0,eps0])


u  , mu_s, epsv   = split(w) 
vu , vmu_s, vepsv = TestFunctions(V)
dw                = TrialFunction(V) 
du , dmu_s, depsv = split(dw) 

wOLD                       = Function(V)
u_old, mu_s_old, epsv_old  = split(wOLD) #, mu_e_old

##### Dirichlet bc
# Displacement
bc1 = DirichletBC(V.sub(0).sub(0), Constant(0.), facets, 3) # ux = 0
bc2 = DirichletBC(V.sub(0).sub(2), Constant(0.), facets, 4) # uz = 0
bc3 = DirichletBC(V.sub(0).sub(1), Constant(0.), facets, 6) # uy = 0
# Enviroment
bc4 = DirichletBC(V.sub(1), Constant(8.430953e6), facets, 1) # μₑ
bc5 = DirichletBC(V.sub(1), Constant(8.430953e6), facets, 2) # μₑ
bc6 = DirichletBC(V.sub(1), Constant(8.430953e6), facets, 5) # μₑ

###### Variational formulation
# Linear momentum
mechanical = inner(Stress(du,depsv,dmu_s),eps(vu))*dx - dot(T,vu)*ds(1)

# Convervation of mass
# volumetric strain is the change of the solvent concentration
# (εₖₖ - εₖₖ₀)/dt + Ωdiv(J) = 0
chemical   = tr(grad(du)-grad(u_old))/Constant(dt)*vmu_s*dx +\
             Constant(Omega)*dot(-D*grad(dmu_s),grad(vmu_s))*dx

# Viscosity 
# (εᵥ - εᵥ₀)/dt - 1/η·σᵥ = 0
viscosity  = inner((depsv-epsv_old)/dt - 1/Gvb*StressB(du,depsv), vepsv)*dx

total_energy = mechanical + chemical + viscosity