# 2D Navier-Stokes with given inlet velocity - resulting pressure field is not physically right

Hello,
I am working with the Navier-Stokes equation - the flow in the rectangular channel. And I was following the (provided example for channel flow. Thanks Dokken for really helpful example.
In my case, I have to consider the flow with given uniform velocity at the inlet = 2 m/s. Also, I should consider the Neumann boundary condition for pressure, but I just started to use FEniCS and didn’t get how to do this - I will appreciate any help with this.

Everything else is the same as is in example. The issue is that the resulting pressure doesn’t look right.
So here are the boundary conditions I have changed:

``````#define the boundary conditions
gdim = domain.geometry.dim
fdim = domain.topology.dim - 1

#walls
def walls(x):
return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], width))

w_b_facets = mesh.locate_entities_boundary(domain, fdim, walls)

u_noslip = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
bcu_noslip = fem.dirichletbc(u_noslip, fem.locate_dofs_topological(V, fdim, w_b_facets), V)

# Inlet
def inflow(x):
return np.isclose(x[0], 0)

class InletVelocity():
def __init__(self, t):
self.t = t

def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0, :] = 2
return values

in_b_facets = mesh.locate_entities_boundary(domain, fdim, inflow)

u_inlet = fem.Function(V)
inlet_velocity = InletVelocity(t)
u_inlet.interpolate(inlet_velocity)
bcu_inflow = fem.dirichletbc(u_inlet, fem.locate_dofs_topological(V, fdim, in_b_facets))

#outlet
def outflow(x):
return np.isclose(x[0], lenght)

out_b_facets = mesh.locate_entities_boundary(domain, fdim, outflow)

bcp_outlet = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(Q, fdim, out_b_facets), Q)

bcu = [bcu_noslip, bcu_inflow]
bcp = [bcp_outlet]
``````

For the images below I used density (rho) = 1 and viscosity (mu) = 1e-2
And here is a resulting pressure:

And resulting velocity:

Also, one more question, when I changed the density = 1000 and dynamic viscosity = 1e-3 (as right for water) the solution seems really unreasonable.

I will appreciate any help, to solve this issue. Thank you.

Please clarify what pressure you were expecting.

I can guess that you were expecting a pressure like in the Poiseuille flow. Note however that your velocity boundary conditions have jumps at the corners on the left, and this results in a singularity of the pressure, which is likely what you are seeing in your plot.

Furthermore, in FEM we typically do not impose strong boundary conditions on p at the outlet.

Changing density/viscosity means increasing a lot the Reynolds number. You may need some stabilization.

Yes, I am expecting the pressure like in Poiseuille flow (decreasing gradient towards the outlet). As you noticed, I have those jumps in the corners, might be because of the high Re number, will try to play with the variables.

I agree about strong boundary conditions for pressure. I prefer to set up the Neumann BC at the outlet, but I cannot figure out how to do it. Could you please help with this?

then you must impose a parabolic profile at the inlet, not a flat one

you don’t need to do anything to set a homogeneous Neumann BC, because it would result in a boundary integral with zero integrand.

Thank you.

So, you mean that I don’t need to set up the outlet boundary conditions for pressure?

Indeed, you don’t need it.

Hi,
I follow your suggestion and run the code with no outlet boundary conditions.
But the solution is not right. With rho = 1, mu = 0.1, inlet velocity = 1, and lengths of the channel = 2 - the Reynolds number = 20, which mean laminar flow:
But here you can see the jumps in the inlet corners, and the values for pressure -9.11e+11 = too much.

So, I do need some boundary conditions for the pressure and with Dirichlet BC I used, at least the values are not that high.

I can’t help much more unless I see the code.

Here are the code. I really appreciate your help. Thank you.

``````from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista

from dolfinx import fem, mesh, io, plot
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc

import ufl
import numpy as np

#define the constants:
lenght = 2
width = 0.5

mu = 1e-2
rho = 1

#create the rectangular mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([lenght, width])],
[41, 6], cell_type=mesh.CellType.triangle)

#time steps
t = 0
Tend = 5
num_steps = 5000
dt = Tend / num_steps

#create the function space
from basix.ufl import element
v_cg2 = element("Lagrange", domain.topology.cell_name(), 2, shape=(domain.geometry.dim, ))
s_cg1 = element("Lagrange", domain.topology.cell_name(), 1)

V = fem.functionspace(domain, v_cg2)
Q = fem.functionspace(domain, s_cg1)

#create the trial and test functions
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
p, q = ufl.TrialFunction(Q), ufl.TestFunction(Q)

#define the boundary conditions
gdim = domain.geometry.dim
fdim = domain.topology.dim - 1

#walls
def walls(x):
return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], width))

w_b_facets = mesh.locate_entities_boundary(domain, fdim, walls)

u_noslip = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
bcu_noslip = fem.dirichletbc(u_noslip, fem.locate_dofs_topological(V, fdim, w_b_facets), V)

# Inlet
def inflow(x):
return np.isclose(x[0], 0)

class InletVelocity():
def __init__(self, t):
self.t = t

def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0, :] = 2
return values

in_b_facets = mesh.locate_entities_boundary(domain, fdim, inflow)

u_inlet = fem.Function(V)
inlet_velocity = InletVelocity(t)
u_inlet.interpolate(inlet_velocity)
bcu_inflow = fem.dirichletbc(u_inlet, fem.locate_dofs_topological(V, fdim, in_b_facets))

#outlet
def outflow(x):
return np.isclose(x[0], lenght)

out_b_facets = mesh.locate_entities_boundary(domain, fdim, outflow)

bcp_outlet = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(Q, fdim, out_b_facets), Q)

bcu = [bcu_noslip, bcu_inflow]
bcp = [bcp_outlet]

#define the variational problem
u_n = fem.Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = ufl.FacetNormal(domain)
f = fem.Constant(domain, PETSc.ScalarType((0, 0)))
k = fem.Constant(domain, PETSc.ScalarType(dt))
mu = fem.Constant(domain, PETSc.ScalarType(mu))
rho = fem.Constant(domain, PETSc.ScalarType(rho))

# Define strain-rate tensor
def epsilon(u):

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

# Define the variational problem for the first step
p_n = fem.Function(Q)
p_n.name = "p_n"
F1 = rho * ufl.dot((u - u_n) / k, v) * ufl.dx
F1 += rho * ufl.dot(ufl.dot(u_n, ufl.nabla_grad(u_n)), v) * ufl.dx
F1 += ufl.inner(sigma(U, p_n), epsilon(v)) * ufl.dx
F1 += ufl.dot(p_n * n, v) * ufl.ds - ufl.dot(mu * ufl.nabla_grad(U) * n, v) * ufl.ds
F1 -= ufl.dot(f, v) * ufl.dx
a1 = fem.form(ufl.lhs(F1))
L1 = fem.form(ufl.rhs(F1))

A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)

# Define variational problem for step 2
u_ = fem.Function(V)
L2 = fem.form(ufl.dot(ufl.nabla_grad(p_n), ufl.nabla_grad(q)) * ufl.dx - (rho / k) * ufl.div(u_) * q * ufl.dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# Define variational problem for step 3
p_ = fem.Function(Q)
a3 = fem.form(rho * ufl.dot(u, v) * ufl.dx)
L3 = fem.form(rho * ufl.dot(u_, v) * ufl.dx - k * ufl.dot(ufl.nabla_grad(p_ - p_n), v) * ufl.dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

# Solver for step 1
solver1 = PETSc.KSP().create(domain.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")

# Solver for step 2
solver2 = PETSc.KSP().create(domain.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.BCGS)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(domain.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

from pathlib import Path
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = io.VTXWriter(domain.comm, folder / "u.bp", u_n, engine="BP4")
vtx_p = io.VTXWriter(domain.comm, folder / "p.bp", p_n, engine="BP4")
vtx_u.write(t)
vtx_p.write(t)

for i in range(num_steps):
# Update current time step
t += dt

# Step 1: Tentative veolcity step
with b1.localForm() as loc_1:
loc_1.set(0)
assemble_vector(b1, L1)
apply_lifting(b1, [a1], [bcu])
set_bc(b1, bcu)
solver1.solve(b1, u_.vector)
u_.x.scatter_forward()

# Step 2: Pressure corrrection step
with b2.localForm() as loc_2:
loc_2.set(0)
assemble_vector(b2, L2)
apply_lifting(b2, [a2], [bcp])
set_bc(b2, bcp)
solver2.solve(b2, p_.vector)
p_.x.scatter_forward()

# Step 3: Velocity correction step
with b3.localForm() as loc_3:
loc_3.set(0)
assemble_vector(b3, L3)
solver3.solve(b3, u_.vector)
u_.x.scatter_forward()
# Update variable with solution form this time step
u_n.x.array[:] = u_.x.array[:]
p_n.x.array[:] = p_.x.array[:]

# Write solutions to file
vtx_u.write(t)
vtx_p.write(t)

# Close xmdf file
vtx_u.close()
vtx_p.close()
b1.destroy()
b2.destroy()
b3.destroy()
solver1.destroy()
solver2.destroy()
solver3.destroy()
``````

I won’t be able to run it for a couple of days now, but at first glance, the inlet profile is not parabolic, which was suggested above.

1 Like

Yes, I have to stuck with this profile for the class project, unfortunately.