# 2D unstationnary plane Poiseuille flow, problem with pressure dirichlet conditions (using DirichletBC method)

This issue is related to Incompressible Navier-Stokes with particular scheme & mixed elements FEniCSx in the sense that the variational formulation is the same. But I am using FEniCS instead of FEniCSx. The case is the 2D unstationnary Poiseuille flow, similar to Test problem 1: Channel flow (Poiseuille flow).

The equations of the problem are:
\frac{\partial\mathbf u}{\partial t}+\mathbf u\overline\otimes\mathbf{\nabla u}=-\frac{1}{\rho}\mathbf\nabla p+\nu\mathbf\Delta\mathbf u\\ \mathbf\nabla\cdot\mathbf u=0

The fluid domain is \Omega=[0, L]\times[-h,h]\subset\mathbb R^2 where h=0.5 is the semi-height of the 2D channel and L=2h is its length. The initial condition is (\mathbf u, p)=0_{\mathbb R^3} inside \Omega. For the boundary conditions we have:

1. a no-slip Dirichlet boundary condition for the velocity at the walls such that \mathbf u(x, y, t)=0_{\mathbb R^2},\;\forall (x,y,t)\in [0,L]\times\lbrace -h,h\rbrace\times\mathbb R_+.
2. a Dirichlet condition for the pressure at the inlet such that p(x,y,t)=p_0>0,\;\forall (x,y,t)\in\lbrace0\rbrace\times[-h,h]\times\mathbb R_+
3. a Dirichlet condition for the pressure at the outlet such that p(x,y,t)=0,\;\forall (x,y,t)\in\lbrace L\rbrace\times[-h,h]\times\mathbb R_+

Theoretically, we impose a constant pressure gradient such that \nabla p=-G\mathbf e_x. Considering some other mathematical constraints, we can derive the steady Poiseuille profile toward which the flow will develop: \mathbf u(x,y,t)=\mathbf u_x(y)\mathbf e_x=\frac{Gh^2}{2\mu}\left(1-\frac{y^2}{h^2}\right)\mathbf e_x, where \mu=\rho\nu.
p(x,y,t)=p(x)=G(L-x), hence p_0=GL.

I use h=0.5, \mu=10^{-2} and G=8\times10^{-2} so that the maximum velocity in the Poiseuille profile is U_0=\frac{Gh^2}{2\mu}=1. Reynolds number is defined as Re=\frac{2h\rho U_0}{\mu}=100.

The variational formulation is the same as in Incompressible Navier-Stokes with particular scheme & mixed elements FEniCSx, so I will not detail it here.

To define the boundary conditions, I use the DirichletBC method, but implementing that in FEniCS leads to a strange solution, far from the expected steady Poiseuille profile. I found a solution on an old FEniCS post that I can’t find anymore, it suggested to weakly impose the pressure. As I understood it, it meant that I had to add some terms in the variational formulation to impose the Dirichlet condition. I have still a lot of things to learn from the finite element method, so I tried something but I am not sure it makes sense. However it seems to work. What I did is add this term inside the weak formulation:

Where \mathbf n is the normal vector to \partial\Omega, \Gamma_\text{inlet}\subset\partial\Omega is the inlet, \mathbf v is the test function for the velocity and \overline\otimes is the simple dot product. Imposing weakly a null pression at the outlet would (as I understand it) be useless because the term \int_{\Gamma_\text{outlet}}\frac{0}{\rho}\mathbf n\overline\otimes\mathbf v\mathrm ds would simply vanish. So for the inlet, I impose the pressure weakly (without using DirichletBC), and for the outlet, I keep using DirichletBC.

I did that because the variational formulation of the incompressible Navier-stokes equations makes the term:
\int_{\partial\Omega}\left(\nu\nabla\mathbf u\overline\otimes\mathbf n-\frac{1}{\rho}p\mathbf n\right)\overline\otimes\mathbf v\mathrm ds naturally appear (by integration by parts), so for me this is how to weakly impose pressure, at least in FEniCS.

My minimum working example is organised as follow:

I have a folder MWE_FEniCS located at the home path (accessed by cd ~/MWE_FEniCS). Inside this folder there is:

• the launch script mwe.py.
• a Mesh folder containing poiseuille_mesh.py
• an empty Results folder

To run the code and see the results you’ll need FEniCS, gmsh, meshio, mpi4py, Paraview (for me, the Xdmf3ReaderS works fine). The run consists of two cases:

• First (problematic) where I use DirichletBC to define all my Dirichlet boundary conditions.
• Second, which seems to work, uses DirichletBC for the velocity at the walls, and the pressure at the outlet. But the Dirichlet condition for the pressure at the inlet is weakly imposed by the adding of (\star) in the variational formulation.

1. Start by going into the Mesh folder and run poiseuille_mesh.py. This will create a .msh file containing the mesh.

2. Then run the script mwe.py two times: first time with the boolean dirichlet set to True (line ~140). This setting will launch the problematic First case. Second time with the same boolean set to False: this setting will launch my “solution” (Second case).

3. You can see the results inside the folder Results. With the command paraview velocity_dirichlet-False.xdmf, you would open the second case velocity solution, which is accurate in my opinion.

Could you please tell me what I am doing wrong in the First case (dirichlet = True) ?

The two scripts needed are written below. If you wish to run my mwe, please make sure you create folders and place the scripts as I detailed above.

Code for mwe.py:

############################################
#--------------- Imports ------------------#
############################################
from dolfin import *
import meshio
from mpi4py import MPI
mesh_comm = MPI.COMM_WORLD
from pathlib import Path
set_log_level(LogLevel.ERROR)

############################################
#----------- Physical parameters ----------#
############################################

# volumic mass
rho = 1

# cinematic viscosity
nu = 0.01

# length of the square domain
L = 1

# Adverse pressure gradient in e_x direction (dp/dx) = -G. The value 0.08 is adopted so that
# the maximum value of the longitudinal velocity is equal to 1
G = 0.08

############################################
#--------- Utility functions --------------#
############################################

# Those to following functions are used to create a mesh & facet xdmf files from the original msh file generated with gmsh
"""
Function that create a meshio.Mesh object from msh mesh using meshio.

Parameters
----------
mesh meshio.Mesh
the mesh object on which the mesh will be written
cell_type : str
see meshio doc
data_name : str
see meshio doc
prune_z: bool
see meshio doc

Returns
-------
meshio.Mesh
the mesh created from the msh mesh file.
"""

# retrieving the cell type (line or triangles)
cells = mesh.get_cells_type(cell_type)

# getting the cell data knowing physical groups of the msh file and the cell type
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)

# retrieving the points corresponding to the cells (the z coordinate is pruned in the 2D case!)
points = mesh.points[:, :2] if prune_z else mesh.points

# creating the meshio.Mesh object from the previously defined parameters
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
data_name: [cell_data]})
return out_mesh

def create_Mesh(mesh_file_path_name):
"""
Method that creates the Mesh object from a msh file. It also creates facet_tags, and updates the
instance parameters attribute 'mesh'.

Returns
-------
the mesh object
"""
mesh = Mesh()

xdmf_file = mesh_file_path_name + ".xdmf"

# xdmf file that will contain the 1D elements mesh (facets)
facet_file = mesh_file_path_name + "_facet.xdmf"

# read mesh from msh file

# create the 1D element mesh
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write(facet_file, line_mesh)

# create the xdmf mesh
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write(xdmf_file, triangle_mesh)

# writing xdmf file into Mesh object
with XDMFFile(xdmf_file) as infile:

# retrieving the facets
mvc1 = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(facet_file) as infile:
facet = cpp.mesh.MeshFunctionSizet(mesh, mvc1)

return mesh, facet

# This function exports the solution into xdmf files (velocity & pressure)
def export_data(xdmf, w, step, first_step=False):

# exporting velocity
xdmf["velocity"].write_checkpoint(w.sub(0), "velocity", step, XDMFFile.Encoding.HDF5, True)

# exporting pressure
xdmf["pressure"].write_checkpoint(w.sub(1), "pressure", step, XDMFFile.Encoding.HDF5, True)

# computing vorticity
u, p = split(w)
vort = project(curl(u), W.sub(1).collapse())

# exporting vorticity
xdmf["vorticity"].write_checkpoint(vort, "vorticity", step, XDMFFile.Encoding.HDF5, True)

############################################
#------------------ Case ------------------#
############################################

"""
dirichlet = True: when the inlet pressure is imposed via DirichletBC method
dirichlet = False: when inlet pressure is weakly imposed
"""

dirichlet = True

############################################
#---------- Problem definition ------------#
############################################

# First we create the mesh objects from the .msh file (The domain is a square (0, -0.5) -- (1, 0.5))

home_path = str(Path.home())

# retrieve the path to the mesh (without its .msh extension)
mesh_file_path_name = home_path + "/MWE_FEniCS/Mesh/poiseuille_mesh"

# creating the xdmf files
mesh, facet = create_Mesh(mesh_file_path_name)

# Now we create the Elements and Spaces for velocity and pressure
velocity_element = VectorElement("Lagrange", mesh.ufl_cell(), 2)
pressure_element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, velocity_element * pressure_element)

# define the relaxation constant (see numerical scheme)
epsilon = Constant(10e-6)

# order of the Backward Differentiation Formula
BDF = 2

# define markers for the boundary conditions
wall_marker, inlet_marker, outlet_marker = 1, 2, 3
markers = {"wall_marker": wall_marker, "inlet_marker": inlet_marker, "outlet_marker": outlet_marker}

# define Dirichlet boundary conditions
bcu_walls = DirichletBC(W.sub(0), Constant((0, 0)), facet, wall_marker) # no slip on the walls
bcp_inlet = DirichletBC(W.sub(1), Constant(G*L), facet, inlet_marker) # pressure at inlet
bcp_outlet = DirichletBC(W.sub(1), Constant(0), facet, outlet_marker) # 0 pressure at outlet

# gather the previously defined boundary conditions
bcs = [bcu_walls, bcp_outlet]
if dirichlet:
bcs += [bcp_inlet]

# define the normal vector on the boundaries
normal = FacetNormal(mesh)

# define the measure to integrate on surfaces
ds = Measure("ds", domain=mesh, subdomain_data=facet)

# Definition of the time step
dt = 1

# Definition of the export path
export_path = home_path + "/MWE_FEniCS/Results/"

# we want to export velocity and pressure
xdmf = dict()
for ele in ["velocity", "pressure", "vorticity"]:
xdmf.update({ele:XDMFFile(mesh_comm, export_path + ele + "_dirichlet-" + str(dirichlet) + ".xdmf")})

# we initialize the solution vector
w = Function(W)

# split solution vector w into velocity (u) and pressure (p)
u, p = split(w)

# initialize finite element test functions, v for velocity and q for pressure
v, q = TestFunctions(W)

# definition of the temporal scheme (here a BDF). A little trick is used to create an adaptative temporal scheme. Meaning that at step 1, BDF1 is used, and after that, BDF2 is used. To better understand, see the numerical scheme detailed in the post.
if BDF >= 1: # BDF1
w_prev1 = Function(W)
u_prev1, p_prev1 = split(w_prev1)
temporal_scheme_list = [dot(u - u_prev1, v)]
assign_list = ['assign(w_prev1, w)']
if BDF >= 2: # BDF2
w_prev2 = Function(W)
u_prev2, p_prev2 = split(w_prev2)
temporal_scheme_list.append(dot(Constant(1.5) * u\
- Constant(2.0) * u_prev1\
+ Constant(0.5) * u_prev2, v))
assign_list.append('assign(w_prev2, w_prev1)')

assign_list = assign_list[::-1]

# function that defines the problem
def define_problem(temporal_scheme):
# definition of the discretization scheme (without BDF)
- Constant(1/rho)* p * div(v)\
- q * div(u)\
- epsilon * p * q) * Constant(dt) * dx
if not dirichlet:
F += dot(Constant(1/rho)*Constant(L*G)*normal, v) * Constant(dt) * ds(markers["inlet_marker"])

# reunion of the BDF (temporal scheme) and F to form the complete variational formulation
F_BDF = temporal_scheme*dx + F

# defining the problem
J = derivative(F_BDF, w)
problem = NonlinearVariationalProblem(F_BDF, w, bcs, J)
solver = NonlinearVariationalSolver(problem)

# specify some solver parameters
solver_settings = {'absolute_tolerance': None,
'convergence_criterion': None,
'error_on_nonconvergence': False,
'linear_solver': "mumps",
'lu_solver': None,
'maximum_iterations': None,
'preconditioner': None,
'relative_tolerance': None,
'relaxation_parameter': None,
'report': True}
for key, value in solver_settings.items():
if value is not None:
solver.parameters["newton_solver"][key] = value
return solver

# Number of steps
n_step = 100

# current step
step = 0

# now that everything is defined, we can start the calculation
while step < n_step:

# if it is the first step that will be calculated, one must apply the boundary conditions to the initial solution
if step == 0:
for bc in bcs:
bc.apply(w.vector())

# write solution on file
export_data(xdmf, w, step, first_step=True)

# define the adatpative temporal scheme
if step < BDF:
if step < BDF - 1:
temporal_scheme = temporal_scheme_list[step - BDF]
else:
temporal_scheme = temporal_scheme_list[-1]

# incrementing the number of step
step += 1

# defining the problem according to the temporal scheme
solver = define_problem(temporal_scheme)

# assign previous current solutions to previous ones (w_prev1 becomes w_prev2 and w becomes w_prev1)
for assignation in assign_list:
eval(assignation)

# solve the problem
try:
n, converged = solver.solve()
except:
raise ExceptionError

print("step number", step, "calculated.")

export_data(xdmf, w, step)

for ele in ["velocity", "pressure", "vorticity"]:
xdmf[ele].close()



Code for poiseuille_mesh.py:

"""
This script takes only care about the gmsh mesh creation
"""

############################################
#------------ Mesh Creation ---------------#
############################################
import gmsh
import time
import sys
import numpy as np
import math

# semi height of the 2D canal
h = 0.5

# length of the canal
L = 2*h

# dimension of the mesh
gdim = 2

# Calculating the time that mesh generation will take
tic = time.perf_counter()

gmsh.initialize()

# target mesh size for any point
lc = 0.04

# test settings
factors = {"walls": 1, "walls_far": 4}

############################################
#---------- Geometry creation -------------#
############################################
## 0 dim : Points definition
# Rectangle delimiting the flow

"""
4----3
|    |
1----2
"""

## 1 dim: Lines and Curves
# outer boundary

gmsh.model.geo.synchronize()

############################################
#-------------- Plane surfaces ------------#
############################################
## 2 dim
# Creation of the plane surface: wake

############################################
#------------- Discretization -------------#
############################################
Nx = 8
Ny = 16
for i in [2, 4]:
gmsh.model.geo.mesh.setTransfiniteCurve(i, Ny)
for i in [1, 3]:
gmsh.model.geo.mesh.setTransfiniteCurve(i, Nx)
gmsh.model.geo.mesh.setTransfiniteSurface(1, "Left", [1, 2, 3, 4])

gmsh.model.geo.synchronize()

############################################
#----- Mesh Saving & Post-processing ------#
############################################
# gmsh.model.geo.mesh.setTransfiniteSurface(1, "Left", [1, 2, 3, 4])

# Generate and write to file
gmsh.model.mesh.generate(gdim)

# ############################################
# #------------ Physical Groups -------------#
# ############################################
# ## Create Physical Groups
gmsh.model.addPhysicalGroup(1, [1, 3], name="walls", tag=1) # walls: up and down sides of the rectangle
gmsh.model.addPhysicalGroup(1, [4], name="inlet", tag=2) # inlet: left side of the rectangle
gmsh.model.addPhysicalGroup(1, [2], name="outlet", tag=3) # outlet: right side of the rectangle
gmsh.model.addPhysicalGroup(2, [1], name="control_volume", tag=4) # Reunion of the two plane surfaces: where fluid will flow

gmsh.model.geo.synchronize()

gmsh.write("poiseuille_mesh.msh")

# gmsh.model.mesh.generate(gdim)
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
gmsh.clear()


Best regards !

Hi @TOP1

Unfortunately, your MWE is too complicated for anyone to quickly test and debug, but I don’t think it makes sense to set DirichletBC for pressure both at the inlet and outlet, as you do in the first case.

There’s only pressure gradient in the Navier-Stokes equations and you can’t set both pressure and pressure gradient to be deterministic in your system. What people typically do is to set pressure at the outlet and velocity at the inlet. Then, as a result of solving for the pressure gradient, you will also get pressure at the inlet. The other option is to give pressure gradient as a body force, which I assume is what you’re doing in the second case.

Well, I fear this is the shortest MWE I can provide. It may seem very long but the code is spaced and commented, also the second script poseuille_mesh.py is there only to generate the mesh and does not need debugging I believe. But I totally understand if it is too long (and also it is using FEniCS and not FEniCSx). However I had to try my luck !

Also you write that you don’t think setting pressure both at inlet & outlet makes sense. I don’t really know about that because it seems that it is exactly what is done in those FEniCS tutorial and FEniCSx tutorial. I don’t see why setting both pressure (inlet & outlet) is problematic in that particular case (starting Poiseuille flow), because it is what will define the pressure gradient and start the (initially static) fluid to flow and goes toward the Poiseuille flow configuration.

For the moment, my “solution” gives good results for the Poiseuille flow configuration, and unless someone can point me what I am doing wrong (in the First or Second cases), I’ll keep the “solution” where I add (\star) in my weak formulation (Second case). For the rest of my work, I have no pressure DirichletBC, hence this “problem” does not appear.