Mixed incompressible Navier-Stokes

Hello there,

I am trying to implement a time-dependent incompressible Navier-Stokes flow, and it ‘sorta works’, but it gets unstable at the outlet when the velocity gets high.

I have tried to run the code on a finer mesh and with smaller time steps, but that does not change anything about the instability - it still happens at the same time.

Therefore I believe something is wrong with my formulation of Naiver-Stokes. Can someone help me find my mistake? I have attached a minimal working code where the problem can be seen.

#region: Load libraries----------------------------
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div, max_value
from progress.bar import Bar
import numpy as np
import mshr
import random
import sys
#import time

#region: define square mesh ------------------------------
domain = mshr.Rectangle(Point(0.0), Point(1,1))
mesh = mshr.generate_mesh(domain, 30)

d = mesh.topology().dim()
mf = MeshFunction("size_t",mesh,d,mesh.domains())
#endregion ----------------------------------------

#region: define boundaries ------------------------
class Top_wall(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[1], 1.0))

class Bottom_wall(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[1], 0.0))

class leftWall(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 0.0))

class rightWall(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 1.0))

#call boundaries
top = Top_wall()
bottom = Bottom_wall()

#mark boundaries
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1,0)
#names for easyness

top.mark(boundaries, Top)
bottom.mark(boundaries, Bottom)

leftwall.mark(boundaries, inlet)
rightwall.mark(boundaries, outlet)

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

boundary_file = File("boundaries.pvd")
boundary_file << boundaries
#endregion ----------------------------------------

#region: define parameters ------------------------
t = 0.0; dt = 0.0256; T = 3
force = Constant((0., 0.))
nu = Constant(0.01)
n = FacetNormal(mesh)

#region: define functionspaces and functions-------
#define function spaces
V = FunctionSpace(mesh, 'RT', 1) # velocity space, Raviart–Thomas space
Q = FunctionSpace(mesh, 'DG', 0) # pressure space, Discontinuous Lagrange space
u1 = Function(V) #for storing velocity solution
p1 = Function(Q) #for current pressure

#define mixed functionspace
P2 = VectorElement("CG", mesh.ufl_cell(), 2)
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
element = MixedElement([P2, P1])
X = FunctionSpace(mesh, element)

# #define functions
xh = Function(X)

tentative = Function(X)
u0, p0 = split(tentative)

#define trial and test functions
(u, p) = TrialFunctions(X) #for velocity, pressure
(v,q) = TestFunctions(X)

#region: define boundary conditions----------------
bcp1 = DirichletBC(X.sub(1), pd, boundaries, inlet)
bcp2 = DirichletBC(X.sub(1), Constant(0), boundaries, outlet)

bcu1 = DirichletBC(X.sub(0), (0,0), boundaries, Top)
bcu2 = DirichletBC(X.sub(0), (0,0), boundaries, Bottom)

bc = [bcp1, bcp2, bcu1, bcu2]

#region: define weak formulation Navier-Stokes------------
# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*nu*epsilon(u) - p*Identity(len(u))

U = 0.5*(u+u0)

#IPCS method, ds over inlet and outlet
F1_ns = dot((u - u0) / dt, v)*dx \
    + dot(dot(u, nabla_grad(u0)), v)*dx  \
    + inner(sigma(U, p), epsilon(v))*dx \
    - dot(force, v)*dx \
    + dot(nabla_grad(p), nabla_grad(q))*dx\
    + (dot(p*n, v)- dot(nu*nabla_grad(U)*n, v))*ds(inlet) \
    + (dot(p*n, v)- dot(nu*nabla_grad(U)*n, v))*ds(outlet)

a1 = lhs(F1_ns)
L1 = rhs(F1_ns)

#region: create save files--------------------------
#get scripts name
name_of_self = sys.argv[0].split('/')[-1].split('.')[0]

#for storing solutions

xdmffile_velocity = XDMFFile(directory_xdmf+'velocity.xdmf') 
xdmffile_pressure = XDMFFile(directory_xdmf+'pressure.xdmf')
xdmffile_temperature = XDMFFile(directory_xdmf+'temperature.xdmf')
xdmffile_velocity.parameters["flush_output"] = True
xdmffile_pressure.parameters["flush_output"] = True
xdmffile_temperature.parameters["flush_output"] = True

#region: solve problem---------------------------
bar = Bar('Processing', max=T/dt)
while t +dt < T + DOLFIN_EPS:
    tentative.vector()[:] = xh.vector()
    u0, p0 = split(tentative)
    dt = min(T-t, dt)

    solve(a1 == L1, xh, bc)

    u1.assign(xh.sub(0, deepcopy=True))
    p1.assign(xh.sub(1, deepcopy=True))

    xdmffile_velocity.write(u1, t); xdmffile_pressure.write(p1, t)

    t += dt

Thank you!

I would suggest starting with verifying your solver using the Taylor-Green vortex problem or the DFG-2D benchmark, which are both covered in: Boundary Conditions & Convergence for 3D Navier-Stokes

After further investigation it has proved, that the instability is a problem with small viscosities. With nu=0.1 there are no problems. @dokken do you know how to work around that?

As the Reynolds number directly relates to the viscosity, this is not suprising. As written in the fenics tutorial.

First of all, it shows that all fluids behave in the same way: it does not matter whether we have oil, gas, or water flowing between two plates, and it does not matter how fast the flow is (up to some criticial value of the Reynolds number where the flow becomes unstable and transitions to a complicated turbulent flow of totally different nature). This means that one simulation is enough to cover all types of channel flow! In other applications, scaling shows that it might be necessary to set just the fraction of some parameters (dimensionless numbers) rather than the parameters themselves. This simplifies exploring the input parameter space which is often the purpose of simulation. Frequently, the scaled problem is run by setting some of the input parameters with dimension to fixed values (often unity).

For turbulence modelling, see for instance:


Hi dokken,

Does that mean that all the equations in Fenics (incl. tutorials) should be (are) given in dimensional values? And if so, what are these dimensions?

Reading from the section you cited, I am confused on why we cannot formulate a weak formulation and a variational problem for the dimensionless Navier-Stokes equation, which is given in this section (3.4.3) too.

Many thanks,

The dimensions in FEniCS is whatever you set them to. If you input a mesh of length 1, it is you as the user who decides if it is given in mm, m or km. They key is that one must make sure every input parameter has the same dimension.

If you read through the rest of section 3.4.3, it states that it is non-trivial to implement the dimensionless form. This is due to the fact if Re is dominating the formulation, it becomes tricky to solve this accurately numerically (Then you should look at turbulence modelling).
As you can see in the implementation of the Navier-Stokes solver, they have not implemented the dimensionless formulation.
If you for instance have a look at the membrane problem (chapter 2.4) in the tutorial, the solve a scaled problem

1 Like

Thank you for that!

I tried to solve the channel problem with H,L=1mm and aiming for a maximum velocity of 0.01m/s.
That results (according to the book) to a Dp of 8e+4! That can be correct, right? The numerical solution blows up almost immediately…

So what is your visocity and density? The Reynolds number less than 100 should be fine.

I have tried both the following case:

  1. rho = 1, mu = 1
  2. rho = 1000 (kg/m3), mu = 0.001 (Ns/m2)

Aiming for velocity of 0.01m/s for a channel of height 10-3 (so 0 < y < H).

According to the book, the value β for the inflow pressure is 8*umax/H^2 (just by changing the assumption U=H=1). But that gives that very high Dp.

The second case will Give an extremely high reynolds number, and is definitely not feasible.

In general, your parameters gives a very high dP as you have made the channel very short, and the pressure has to drop from 8 to 0 in 1e-3 meters.

The second case gives a Reynolds number of 10, why is that high?

I want to find the “correct” Dp in order to get that desired max velocity.

Could a potential longer channel fix the issue?

Thank you :slight_smile:

I would think so yes. Have you tried increasing the size of the channel to for instance 1?

I still do not understand why for the second case you consider the Re extremely high. I computed as 10.

Yes, it works better. However, I notice that the maximum velocity decreases significantly, compared to the unit square case with max.vel=1 for the 2nd case of density and viscosity.

Right, the Reynolds number is 10 in that case, my mistake.
I would suggest you start by varying one parameter in your simulation from the scaled problem to your case, and see how that affects the flow field, instead of varying them all at once.
For instance, start with the standard case L=H=\rho=\mu=U=1 and start varying the height/length of the channel, and consider the solution.
You also have to satisfy the CFL condition: https://en.wikipedia.org/wiki/Courant–Friedrichs–Lewy_condition

1 Like

I am using the tutorial for channel flow in which I think the variational problem is not in scaled form.
Would I need to rewrite that?

Noted! :slight_smile:

No, just simply set the unscaled implementation to use these parameters as a starting point.

Oooh okay, I asked because you mentioned the scaled problem to my case.
I have a more specific case study I want to look at, hence the jump to these values.

I will definitely try your suggestion! Thanks again :slight_smile:


Is there any way I can estimate the Dp that corresponds to the maximum velocity I want?

Many thanks,

See: Section 3.4.3 of the FEniCS tutorial, which goes through how to estimate the pressure difference and maximum velocity for the channel flow.

Hi dokken,

Unfortunately, I find this paragraph (and section) very confusing because it starts assuming a non-dimensional formulation and then returns to dimensions.It presents a characteristic velocity U, which it returns only as 1 (as being non-dim).

Assuming that my box is square and H=L=0.001, that gives me a Dp=4x10^5, which I don’t think it can be right. I used the following formula for that:

> u=(1/2)*beta*y*(H-y)

where beta=Dp (just for consistency with the tutorial).

What am I missing here…?

Many thanks,