I was trying to solve thixotropic fluid flow. Thixotropy is a fluid with a varying viscosity according to the internal molecular structure that takes a scalar value in between 0 to 1. The set of governing equation consists of momentum and mass conservation + equation that defines the balance of molecular structure.
The equation of the last effect is defined by a shear breakdown of molecular structure (\xi) and the brownian build up.
where \gamma_s is the second invariant of strain-rate tensor. For example, in Fenics code, it can be defined as
def dgamma(u):
gamma = 2*sym(grad(u))
gamma_invariant = sqrt( 0.5 * inner(gamma, gamma))
return gamma_invariant
However, I found that when the above user expression is included in the variational form, the the iterative solver speaks out ‘nan’ from the first iteration.
I am very suspicious on the above dgamma(u) term since while this is not working
choice 1 --not working
F = (
rho_0*inner(grad(u)*u, u_)*dx + viscosity(xi)*inner(grad(u), grad(u_))*dx
- inner(p, div(u_))*dx + inner(p_, div(u))*dx
+ inner(f, u_)*dx
)
F += inner(-kd*xi * **dgamma(u)** + ka*(1-xi),xi_)* dx
However, The following code is working, , (i.e. it converges to some solution)
choice 2 working
----------------------------------
F = (
rho_0*inner(grad(u)*u, u_)*dx + viscosity(xi)*inner(grad(u), grad(u_))*dx
- inner(p, div(u_))*dx + inner(p_, div(u))*dx
+ inner(f, u_)*dx
)
F += inner( -kd*xi + ka*(1-xi),xi_)* dx
I note that the only difference is whether I included “dgamma(u)” or not.
I am attaching a minimum working & not working code example in the follow.
(It considers the case lid-driven cavity flow. )
I appreciate your help in advance.
Full code
from fenics import *
#hyper parameters
rho_0 = 1
mu = 0.01
Ux = 1
kd=0.01
ka=1.0
mesh = UnitSquareMesh(40, 40, diagonal="crossed")
V_ele = VectorElement("CG", mesh.ufl_cell(), 2)
Q_ele = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([V_ele, Q_ele, Q_ele]))
lid_location = "near(x[1], 1.)"
velocity_x = Constant((Ux, 0.0))
reg_vel_x = Expression( ("16*x[0]*(1-x[0])", "0.0") , degree=1)
top_velocity = DirichletBC(W.sub(0), reg_vel_x , lid_location)
fixed_wall_locations = "near(x[0], 0.) | near(x[0], 1.) | near(x[1], 0.)"
noslip = DirichletBC(W.sub(0), (0, 0), fixed_wall_locations)
bcu = [top_velocity, noslip]
u_, p_ , xi_ = TestFunctions(W)
upxi = Function(W)
u, p, xi = split(upxi)
f = Constant((0, 0)) # force
#Initial Newtonian Solution
def viscosity(xi):
return 1.0 + 0.001 * xi
def dgamma(u):
gamma = 2*sym(grad(u))
gamma_invariant = sqrt( 0.5 * inner(gamma, gamma))
return gamma_invariant
F = (
rho_0*inner(grad(u)*u, u_)*dx + viscosity(xi)*inner(grad(u), grad(u_))*dx
- inner(p, div(u_))*dx + inner(p_, div(u))*dx
+ inner(f, u_)*dx
)
#choice 1 --not working
F += inner(-kd*xi * dgamma(u) + ka*(1-xi),xi_)* dx
#chocie 2 --working
#F += inner( -kd*xi + ka*(1-xi),xi_)* dx ### wokring...
J = derivative(F, upxi)
solve(F == 0, upxi, bcu, solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}})
(u,p,xi) = upxi.split(True);
ofile_u = File("thixo_u.pvd")
ofile_u << u
print("Code is sucessfully finished...")