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:
- 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_+.
- 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_+
- 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:
\int_{\Gamma_\text{inlet}}\frac{GL}{\rho}\mathbf n\overline\otimes\mathbf v\mathrm ds\qquad(\star)
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
Meshfolder containingpoiseuille_mesh.py - an empty
Resultsfolder
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
DirichletBCto define all my Dirichlet boundary conditions. - Second, which seems to work, uses
DirichletBCfor 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.
Please follow these 3 steps to run my mwe:
-
Start by going into the
Meshfolder and runpoiseuille_mesh.py. This will create a.mshfile containing the mesh. -
Then run the script
mwe.pytwo times: first time with the booleandirichletset toTrue(line ~140). This setting will launch the problematic First case. Second time with the same boolean set toFalse: this setting will launch my “solution” (Second case). -
You can see the results inside the folder
Results. With the commandparaview 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
def create_mesh(mesh, cell_type, data_name="name_to_read", prune_z=False):
"""
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
msh = meshio.read(mesh_file_path_name + ".msh")
# 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:
infile.read(mesh)
# retrieving the facets
mvc1 = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(facet_file) as infile:
infile.read(mvc1, "name_to_read")
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))
# retrieve your home path.
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)
F = (dot(dot(grad(u), u_prev1), v) + dot(dot(grad(u_prev1), u), v) - dot(dot(grad(u_prev1), u_prev1), v)\
+ Constant(nu) * inner(grad(u), grad(v))\
- 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()
gmsh.model.add("poiseuille_plane_2D")
# 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
gmsh.model.geo.addPoint(0, -h, 0, tag=1)
gmsh.model.geo.addPoint(L, -h, 0, tag=2)
gmsh.model.geo.addPoint(L, h, 0, tag=3)
gmsh.model.geo.addPoint(0, h, 0, tag=4)
"""
4----3
| |
1----2
"""
## 1 dim: Lines and Curves
# outer boundary
gmsh.model.geo.addLine(1, 2, tag=1)
gmsh.model.geo.addLine(2, 3, tag=2)
gmsh.model.geo.addLine(3, 4, tag=3)
gmsh.model.geo.addLine(4, 1, tag=4)
gmsh.model.geo.synchronize()
############################################
#-------------- Plane surfaces ------------#
############################################
## 2 dim
# Creation of the plane surface: wake
gmsh.model.geo.addCurveLoop([1, 2, 3, 4], tag=1)
gmsh.model.geo.addPlaneSurface([1], tag=1)
############################################
#------------- 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 !