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

Ksenia_Stern · May 15, 2024, 11:54pm

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.

francesco-ballarin · May 16, 2024, 5:11am

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.

Ksenia_Stern · May 16, 2024, 2:39pm

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?

francesco-ballarin · May 16, 2024, 3:30pm

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.

Ksenia_Stern · May 16, 2024, 5:10pm

Thank you.

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

francesco-ballarin · May 16, 2024, 6:32pm

Indeed, you don’t need it.

Ksenia_Stern · May 20, 2024, 4:35pm

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.

francesco-ballarin · May 21, 2024, 5:38am

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

Ksenia_Stern · May 21, 2024, 5:24pm

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):
    return ufl.sym(ufl.nabla_grad(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)
a2 = fem.form(ufl.dot(ufl.nabla_grad(p), ufl.nabla_grad(q)) * ufl.dx)
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])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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()

francesco-ballarin · May 21, 2024, 5:54pm

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.

Ksenia_Stern · May 21, 2024, 6:08pm

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

Topic		Replies	Views
Initial conditions for Navier--Stokes equations General	6	421	March 8, 2023
Errors when using symmetry boundary conditions in 2D stokes channel flow General	2	547	December 20, 2022
Encountering problem in implementing a pressure boundary condition in N-S equations General	0	13	July 27, 2024
Boundary Conditions & Convergence for 3D Navier-Stokes	14	4535	June 26, 2020
Obtaining Incorrect Results for Flow Over a Cylinder - Navier Stokes variational formulation	2	928	July 21, 2021

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

Related topics