Nonlinear Non-Newtonian Navier Stokes

Question

Hello!

I want to solve the Navier Stokes equations with the Nonlinear solver!
I have written the code below for a pipe (0.25*0.05) ,but i can't get a solution...
Can you help me?

Thanks in advance!

from dolfin import *
# Define mesh and geometry
mesh = RectangleMesh(Point(0, 0), Point(0.25, 0.05),40,10)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
V = VectorFunctionSpace(mesh, "CG" , 2)
Q = FunctionSpace(mesh , "CG", 1)
W = MixedFunctionSpace([V, Q])

# Define Function and TestFunction(s)
w = Function(W); (u, p) = split(w)
(v, q) = TestFunctions(W)
u_prev=Function(V)
u_prev=interpolate(Constant((0.0,0.0)),V)

#----------------------------- Define boundaries------------------------------
inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 0.25)'
walls   = "on_boundary && !(near(x[0], 0.0) || near(x[0], 0.05))"

# Define inflow profile
inflow_profile = Expression(('4.0*0.01*x[1]*(0.05 - x[1]) / pow(0.05, 2)', '0'),degree=2)

# Define boundary conditions
bcu_noslip = DirichletBC(V, Constant((0, 0)), walls)
bcu_inflow = DirichletBC(V, inflow_profile, inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip, bcu_inflow]
bcp = [bcp_outflow]
bcs=[bcu_noslip,bcu_inflow,bcp_outflow]

#--------------------------------Setup the problem ---------------------------
f=Constant((0.0,0.0))
dt=1E-3
k=Constant(dt)

# Define variational form
def epsilon(input):
    return sym(grad(input))

def gamma(w):
    return pow(2*inner(epsilon(w),epsilon(w)),0.5)

def viscocity(w):
    C1=Constant('0.00148')
    C2=Constant('5.12')
    C3=Constant('0.499')
    H=Constant('0.4')
    kappa=C1*exp(C2*H)
    n_exp=1-C3*H

    return kappa*pow(gamma(w),n_exp-1)

F = 1./k*inner(u-u_prev,v)*dx+(2*viscocity(u)*inner(epsilon(u), epsilon(v)) - div(u)*q - div(v)*p-inner(f,v))*dx#+(p0-2*nu*epsilon(u0))*dot(v,n)*ds

# Solve problem
solve(F == 0, w, bcs)

Answer 1

Your problem lies in the choice of the initial guess for the nonlinear problem. Please see below for an updated code (look for comments starting with XXX). I am not sure whether the obtained solution is physical, but this should give you at least a starting point to understand the source of your problem.

from dolfin import *
# Define mesh and geometry
mesh = RectangleMesh(Point(0, 0), Point(0.25, 0.05),40,10)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
# XXX updated to 2016.1.0
V_element = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q_element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W_element = MixedElement(V_element, Q_element)
W = FunctionSpace(mesh, W_element)

# Define Function and TestFunction(s)
w = Function(W); (u, p) = split(w)
(v, q) = TestFunctions(W)
# XXX define a mixed prev function
w_prev=Function(W); (u_prev, p_prev) = split(w_prev)
# XXX there is no need to project the 0. constant, Function is already zero by default

#----------------------------- Define boundaries------------------------------
inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 0.25)'
walls   = "on_boundary && !(near(x[0], 0.0) || near(x[0], 0.05))"

# Define inflow profile
inflow_profile = Expression(('4.0*0.01*x[1]*(0.05 - x[1]) / pow(0.05, 2)', '0'),degree=2)

# Define boundary conditions
# XXX using sub(0) and sub(1) in BCs
bcu_noslip = DirichletBC(W.sub(0), Constant((0, 0)), walls)
bcu_inflow = DirichletBC(W.sub(0), inflow_profile, inflow)
bcp_outflow = DirichletBC(W.sub(1), Constant(0), outflow)
bcu = [bcu_noslip, bcu_inflow]
bcp = [bcp_outflow]
bcs=[bcu_noslip,bcu_inflow,bcp_outflow]

#--------------------------------Setup the problem ---------------------------
f=Constant((0.0,0.0))
dt=1E-3
k=Constant(dt)

# Define variational form
def epsilon(input):
    return sym(grad(input))

def gamma(w):
    return pow(2*inner(epsilon(w),epsilon(w)),0.5)

def viscosity(w):
    C1=Constant('0.00148')
    C2=Constant('5.12')
    C3=Constant('0.499')
    H=Constant('0.4')
    kappa=C1*exp(C2*H)
    n_exp=1-C3*H

    return kappa*pow(gamma(w), n_exp-1)

# XXX you cannot start the nonlinear iteration from w=0, because gamma(0) = 0 and n_exp-1 < 0
# XXX solve e.g. a Stokes problem first
(u_stokes, p_stokes) = TrialFunctions(W)
F_stokes = (2*inner(epsilon(u_stokes), epsilon(v)) - div(u_stokes)*q - div(v)*p_stokes-inner(f,v))*dx
solve(lhs(F_stokes) == rhs(F_stokes), w, bcs)

# XXX now solve the nonlinear problem
F = 1./k*inner(u-u_prev,v)*dx+(2*viscosity(u)*inner(epsilon(u), epsilon(v)) - div(u)*q - div(v)*p-inner(f,v))*dx
solve(F == 0, w, bcs)

# XXX assign current solution for the next time iteration
assign(w_prev, w)

Comments

Great!!
Thank you very much!
I have found another way and i will compare the 2 solutions!

Thank you again . You were very helpful!

Answer 2

The code is not complete and you do not include the error messages. Hence, it is difficult to determine the cause of your problem.

Comments

Hello! Thank you for your answer ! I have edited the initial code.
I have made some changes in the code and simplified my problem !
Basically i took the Stokes problem and adapted it in order to get the non-newtonian Navier stokes.
The problem is that i get a nan answer in the residuals.
(In the boundaries i have the integral of pnv*ds which i thought was equal to zero)

Answer 3

Is there anyone who can help me with the problem?

I think the problem is the PDE form .
I think i have missed the stress vector (-p+T)*n on the boundary but i don't know what is its value.
And the second mistake is the form of the viscocity .

I just need some help with the PDE and i think i can solve all the other problems.

Thanks in advance!