# Nonlinear variational problem with user specified function

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_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 = (
- 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 = (
- 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. )

# 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_invariant = sqrt( 0.5 * inner(gamma, gamma))
return gamma_invariant

F = (
- 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...")



rease format all code snippets with 3x syntax, ie

python

`