# Taylor-Hood finite hexahedral elements, pressure diverging

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?

• 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.
– knl
Commented Sep 30, 2018 at 16:53
• @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. Commented Oct 1, 2018 at 10:59
• @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... Commented Jun 19, 2020 at 5:58

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.