I would like to solve 3D transient incompressible Navier-Stokes with FEM, Newton method, Schur-based preconditioner, Lagrangean P2/P1 elements (no stabilization), in a rigid pipe discretized with tetrahedrons.
The fluid is initially at rest, I impose constant inlet and outlet pressure boundary conditions (and therefore constant pressure gradient), and I would like to reach Poiseulle steady-state solution after a while. Reynolds number is within the laminar regime, and CFL holds! After couple hundred of timesteps, pressure becomes unstable and everything breaks.
I impose pressure b.c. by adding -$\int_\Gamma p \eta N_i dx$ at the moment equation, where $\eta$ is the normal outward unitary vector, and $N_i$ is the $P_2$ basis function associated to the boundary nodes. Nothing out of the standard, I think.
However, when I visualize the result of this integral in Paraview, the boundary looks like a checkerboard: there are nonzero values at the d.o.f. associated to the edges, but nearly ZERO values at the d.o.f associated to the vertices! The explanation I used to convince myself is that the shape functions associated to vertices have zero integral when they are restricted to triangles (actually, unless I have made a bizarre mistake, this also holds for $P_2$ shape functions associated to the vertices of the triangles have zero integral).
P.S.: the velocity profile around the last successful time-steps before stoping has also a kind of chessboard profile.
Few issues puzzle me:
1) everything works fine if I keep the same pressure gradient, mesh, Reynolds number etc. but impose parabolic velocity profile (given by Poiseulle formulae), or constant velocity (and then I have the Poiseulle profile after the entrance length).
2) the code which evaluates the pressure integral above is the exactly the same used by another (solid mechanics) solver, validated with inhomogeneous & transient traction boundary conditions, which uses P1 elements instead!
3) if I refine the mesh, I still have the checkerboard pattern, but the pressure imposition setup somehow works! I did not expect mesh dependency. To be clear, this checkerboard pattern is not the pressure, but the result of the integral above.
My quite vague question is: what am I missing?
Let me know if the problem statement is clear.
Btw., do you know any other way to impose pressure? I've already heard about introducing 1's at the (nonexistent) boundary pressure entries of the continuity equation ($\nabla \cdot v = 0$), but I am uncomfortable with this idea, if I don't use any stabilization.
Any hint is highly appreciated!