I am developing a FEM fluid solver using the Taylor-Hood algorithm, i.e. quadratic interpolation for velocity, and linear for pressure.

I have developed the code for 2-D quadrilaterals and triangles, and everything worked fine. I am trying to do it in 3D, using a combination of 20-noded and 8-noded hexahedrals, and I have the following problem. The velocity solution is correct and converging, but the pressure solution is very unstable. Pressure is not converging. It doesn't diverge either. I tried with under-relaxation and the problem remained.

The code looks alright (I checked it several times), but there's always a possibility to have a problem somewhere. Am I missing something in the physics? The pressure is zeroed at one point, and the problem I am trying to solve is the simple 3D lid driven cavity flow, where velocity is prescribed on all boundaries.

Can you suggest me any relative document?

  • $\begingroup$ It might not be the problem, but I'm always suspicious about fixing point values of $L^2$-functions. Mathematically speaking $L^2$ functions - as the pressure is in the typical weak formulation used for FE implementation - do not have well defined point values. Hence, strictly speaking, you should not be able to specify pressure at "one point". It might work in some cases but I don't know if it works in general. $\endgroup$
    – knl
    Sep 30, 2018 at 16:53
  • $\begingroup$ @knl: That's definitely a point of theoretical contention. One would expect that as the mesh size goes to zero, this leads to an ill-posed problem. But it works well enough in practice, even on very fine meshes. $\endgroup$ Oct 1, 2018 at 10:59
  • $\begingroup$ @knl it sounds like you would be then calling into question the entire finite element method, but the pressure solutions are at least in $H^1$ which means they are a very special subclass of $L^2$ function (dense in the space of continuous functions from memory). I thought the usual argument was that the solution belongs to an $L^2$ equivalence class, and there is always a continuous representative from which you can take point values without excessively cheating. Not sure what numerical integration techniques don't use point values anyway... $\endgroup$ Jun 19, 2020 at 5:58

2 Answers 2


It's of course hard to tell without knowing the code, the testcase, and your results. So I can't really help you out here.

But it is worth pointing out that you are making yourself too much work. Any of the big, open source finite element libraries have these elements already built in and exceedingly well tested. They also have adaptive meshes, excellent solvers, and parallelism available. There is really no reason for you to implement all of this yourself -- save your time and energy to implement the things that are specific to your problem, not generic.

(Disclaimer: I am one of the authors of one of these large libraries -- deal.II. You may be interested in looking at the Stokes equation tutorial program. You can make your way from there to the Navier-Stokes solvers.)


I resolved my problem. I added a small pressure coefficient in the continuity equation, let's say $10^{-10}$. In the weak form, this quantity is multiplied by the integral {$N^T$ $N$} (or mass matrix), where $N$ are the shape functions of the 8 noded hexahedral. The code now works fine.

  • $\begingroup$ You can self-accept your answer. I think it is valuable. $\endgroup$
    – Anton Menshov
    Oct 1, 2018 at 18:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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