2
$\begingroup$

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!

Thanks!

|cite|improve this question
$\endgroup$
  • $\begingroup$ Are you applying this mean pressure gradient only on the boundary? How are you enforcing $\partial_i u_i=0$? To be clear, this checker boarding disappears and your solution converges if your mesh is fine enough? $\endgroup$ – Charles Jun 18 '16 at 5:49
  • $\begingroup$ I apply only on the boundary, but it is not a pressure gradient, it is the pressure itself instead. I use mixed formulation, and my system of equations looks like $[A B; B^T 0]*(v,p) = (f,0)$. The divergence-free is imposed by the block $B^T * v = 0$. The checkerboard never disappears, however if the mesh is fine enough, the solver simply works (I will edit and make it clear in the question). $\endgroup$ – Frederico Teixeira Jun 18 '16 at 8:20
  • $\begingroup$ I suggest writing your BCs for pressure and velocity in an equation. Simply saying pressure BCs at the inlet and outlet is vague. Also, because of this, I'm not exactly sure what is driving your flow. $\endgroup$ – Charles Jun 18 '16 at 16:52
  • $\begingroup$ Thanks for your suggestion. I will edit it soon. Meanwhile, just to clarify that my fluid is driven by pressure difference between inlet and outlet. E.g., set $p_{in}=10 Pa$ and $p_{out} = 5 Pa$ in a pipe of length $L=1m$ aligned with the z-axis. Then, $\dfrac{\partial p}{\partial z} = 5 Pa/m$. This value is already enough to define the velocity of the Poiseulle solution. $\endgroup$ – Frederico Teixeira Jun 18 '16 at 20:03
2
$\begingroup$

As I know, a natural setup for FEM solution of incompressible flow is : velocity specified at inlet and homogeneous condition at outlet for $$\vec{n}\cdot(-\nu\nabla u+pI)|_{\Gamma_N}=0$$.

Then the latter can be applied for facets of outlet as so-called natural BCs.

Hope this helps you.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ That is a completely different physical problem. $\endgroup$ – Vladimir F Apr 13 '17 at 8:17
  • $\begingroup$ This answer is correct. One cannot set Dirichlet conditions for $p$ in the standard finite element formulation for Navier-Stokes. Since $p$ is $L^2$ function, pointwise values are not well-defined. $\endgroup$ – knl May 8 '17 at 13:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.