Solving non-newtonian fluid

1

Dear FEniCS community. Could you help me point out whats wrong with my coding? My first itearation is already nan. Error that appears:

Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** DOLFIN version: 2017.2.0
*** Git changeset: unknown

My coding:

 from __future__ import print_function
 from fenics import *
 import matplotlib.pyplot as plt

    # define dimensionless parameter
    Re = 1000
    Pr = 0.71
    n = 3

    # create mesh
    mesh = RectangleMesh(Point(0,0),Point(2,1),32,32)

    # define function space
    V = VectorElement('P', mesh.ufl_cell() , 2)
    Q = FiniteElement('P', mesh.ufl_cell(), 1)
    element = MixedElement([V, Q])
    W = FunctionSpace(mesh, element)

    # define boundaries
    walls = 'x[1] < DOLFIN_EPS'
    inflow = 'x[0] < DOLFIN_EPS'
    outflow = 'x[0] > 2 - DOLFIN_EPS'
    upper = 'x[1] > 1 - DOLFIN_EPS'

    # define boundary condition

#walls
bcu_walls = DirichletBC(W.sub(0), Constant((0,0)), walls)

#inflow
bcu_inflow = DirichletBC(W.sub(0), Constant((1,0)), inflow) #u_e=1

bcu = [bcu_walls, bcu_inflow]

bcs = bcu 

#define test functions
(v, q) = TestFunctions(W)

# define functions
u = Function(W)
p = Function(W)

# split functions
w = Function(W)
(u, p) = split(w)

# define expression used in variational forms
Re = Constant (Re)
Pr = Constant (Pr)
n = Constant(n)

# define epsilon
def epsilon(w):
    return sym(nabla_grad(w))

# define secondinvariant
def secondinvariant(w):
    return pow(0.5*inner(epsilon(w),epsilon(w)),0.5)

# define powerlaw
def powerlaw(w):
	return pow(secondinvariant(w),n-1)

# define variational problem
F1 = div(u)*q*dx 

F2 = dot(dot(u,nabla_grad(u)),v)*dx + dot(grad(p),v)*dx \
+ 2*pow(Re,-1)*powerlaw(u)*inner(epsilon(u),epsilon(v))*dx 

F = F1 + F2 

# solve nonlinear variational problem
J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['linear_solver'] = 'lu'
prm['newton_solver']['convergence_criterion'] = 'residual'
prm['newton_solver']['absolute_tolerance'] = 1E-10
prm['newton_solver']['relative_tolerance'] = 1E-9
prm['newton_solver']['maximum_iterations'] = 10
solver.solve()

# Create VTK files for visualization output
vtkfile_u = File('steady_powerlaw/u.pvd')
vtkfile_p = File('steady_powerlaw/p.pvd')

# Save solution to file (VTK)
(u, p) = w.split()
vtkfile_u << u
vtkfile_p << p

# Plot solution
plot(u, title='Velocity')
plot(p, title='Pressure')    

# Hold plot
plt.show()

Thank you for your help. It really helps me a lot.

2

You can fix the nan by adding a very small regularization to the square root:

# define secondinvariant
def secondinvariant(w):
   return pow(0.5*inner(epsilon(w),epsilon(w)) + DOLFIN_EPS,0.5)

You’ll also want to integrate the \nabla p\cdot\mathbf{v} term in F2 by parts (and maybe add a desired nonzero pressure boundary condition on the Neumann/outflow boundaries). The missing boundary condition on the subdomain defined by upper may also be a problem; the code didn’t converge for me unless I either set

walls = 'x[1] < DOLFIN_EPS || x[1] > 1-DOLFIN_EPS'

or added some backflow stabilization on the Neumann boundaries.

1 Like
3

Thank you for your help. The reason my boundary condition becomes like that because I am trying to solve the problem past a vertical plate. If the boundary is a problem, what are the proper boundary for me to choose for the case of vertical plate?