How to add a outlet conditon?

[HYB777]

In the demo of tutorials/02_navier_stokes, I want to add a outlet condition "p=0" at boundary "l3".

The outlet pressure define are follow.

def p_out_eval(x: np.typing.NDArray[np.float64])-> np.typing.NDArray[petsc4py.PETSc.ScalarType]:
"""Return the zero pressure at the outlet."""
return np.zeros((1, x.shape[1]))

The following code is what I modified.

p_out = dolfinx.fem.Function(Q) # new
p_out.interpolate(p_out_eval) # new
bdofs_V_1 = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundaries_1)
bdofs_V_2 = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundaries_2)
bdofs_Q_3 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_3) # new
inlet_bc = dolfinx.fem.dirichletbc(u_in, bdofs_V_1)
wall_bc = dolfinx.fem.dirichletbc(u_wall, bdofs_V_2)
outlet_bc = dolfinx.fem.dirichletbc(p_out, bdofs_Q_3) # new
bc = [inlet_bc, wall_bc, outlet_bc]

But the results is not convergence and it terminate after only one iteration.

[francesco-ballarin]

In typical finite element formulation for CFD you never impose a Dirichlet condition on the pressure. You instead impose a Neumann condition on p n - \partial u / \partial n = 0.
For the geometry in tutorial 2 the normal n for the outlet section is the unit vector aligned with x, this means that \partial u / \partial n is actually \partial u / \partial x, so the Neumann condition is actually p n - \partial u / \partial x = 0

If you want to use that to impose p = 0, make sure to change the mesh so that the channel is long enough that if you take a look at the velocity at the outlet and at a vertical section just before the outlet, they look indistinguishable. If you manage to achive that, then \partial u / \partial x will be zero, and thus in practice the Neumann condition will enforce p = 0.

[HYB777]

Thanks for your reply!

But I have another problem, for the demo in tutorial 2, I impose a symmetry boundary condtion on "l4", which means u.y=0, the relevant part of the code as follows.

gmsh.model.addPhysicalGroup(1, [l4], 4)
boundaries_4 = boundaries.indices[boundaries.values == 4]
bdofs_V_4 = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundaries_4)
symmetry_bc = dolfinx.fem.dirichletbc(dolfinx.default_scalar_type(0), bdofs_V_4, V.sub(1))
bc = [inlet_bc, wall_bc, symmetry_bc]

But the result seem to be incorrect. The velocity result are shown as follow.

Here is my complete code
stationary_ns.zip

[francesco-ballarin]

I would suggest to first prepare the axysimmetric geometry and verify the solution there, and afterwards move to the axysimmetric case

[HYB777]

Thanks for your message! By the way, I found that if I choose small kinematic viscosity (such as 1e-6), the result seems not convergence and it terminate after only one iteration. And what should be done about this situation?

[francesco-ballarin]

The Reynolds number is very large in that case, you'll need to add some stabilization to your FE scheme. A possible acronym to search for is SUPG, but there is a whole research field of techniques working on this.