I am trying to compute a 3D channel flow with a parabolic inflow-profile, homogeneous Dirichlet BCs at the walls and a homogeneous Neumann BC at the outflow. For now I am discretising in time using IMEX Euler and using Taylor-Hood elements in space. My code is based on a similar 2D implementation (and which works without problems in 2D), however in 3D the pressure explodes after one time step at the point (0,0,0).
I’ve checked the BC names and they are all where they are meant to be.
With the Neumann do-nothing BC at the outflow, the pressure should be well defined.
Using a different sparse solver didn’t help.
Increasing the pressure regularisation term helps but does not remove the problem.
Is there anything different in 3D compared to 2D (in NGSolve) which I haven’t taken into account?
think the problems comes from the missing LBB-stability of Taylor Hood for general unstructured meshes.
Try to add bubbles to the velocity space:
V = VectorH1(mesh, order=k, dirichlet="wall|inflow")
print ("ndofV = ", V.ndof)
My recommendation is to use HDG (or MCS) methods instead
thank you! That seems to have solved the issue!
I wasn’t aware, that Taylor-Hood elements have inf-sup stability problems and thought that besides their relative simplicity, the inf-sup stability of TH was one the major advantage of these elements. Also the meshes from Netgen have always looked nice and regular to me, which makes the lack of stability here even more surprising to me.
Apart from my belief that Taylor-Hood would be inf-sup stable, one of my main motivations to use Taylor-Hood was to try and keep things simple at this point. But I am aware that HDG would bring some advantages!
I think, In principle you can already have problems in 2D with Taylor-Hood if you have Dirichlet boundary conditions and the element in a corner has no interior vertex. Similar statements hold in 3D. You can always make meshes which have an interior vertex for elements on a corner. Netgen meshes however do not guarantee that from the start. You could make subdivisions manually of a corner element to achieve the interior point condition.
A simple trick to avoid elements which have all edges on Dirichlet boundaries is to take a mesh and apply one NGSolve-Refinement. So another quick fix for your script could be the following:
# ----------------------- PARAMETERS -----------------------
h_max = 0.22 # Global mesh size
thank you! Refining also worked The lack of an interior vertex in the original mesh also explains why the pressure instability was in one of the corners.