Finite element stabilization schemes for incompressible flow

Question

I am looking for an easy to implement stabilization scheme that can be used with equal order ($P_1-P_1$ or $Q_1-Q_1$) finite elements for fluid flow. Is there something like this or should I stick to Taylor Hood elements? It appears that there are a number of methods but the ease of implementation is not quite clear.

Answer 1

Warning

First, be warned that stabilization is generally much more fragile than use of stable elements. If you solve a narrow class of problems, you might be able to choose reasonable parameters and get on with your work, but if you solve a wide range of problems or need to solve demanding under-resolved scenarios, especially with less regular boundaries, coefficient structure, and constitutive relations, stabilization will generally require constant vigilance to detect and quantify the impact of numerical artifacts, and to re-tune parameters for.

For these reasons, I strongly recommend stable elements (especially those with discontinuous pressure space, thus locally conservative) if you want a robust method that does not need to be frequently revisited. In my opinion, reliability and clean analysis is worth the modest increase in data structure complexity.

Stabilization methods

I'll outline the two most popular choices in the context of a stabilization term $S_h(\mathbf u_h, p_h, \mathbf f; \mathbf v_h, q_h)$ to be added to the weak form.

• Pressure-stabilized Petrov-Galerkin (PSPG) This stabilization has the form $$ S_h^{\text{pspg}}(\mathbf u_h,p_h;q_h) = \sum_{T \in \mathcal{T}_h} \alpha_T\Big(-\nabla\cdot(\eta D \mathbf u_h) + \nabla p_h - \mathbf f, \nabla q_h\Big)_T$$ where the inner product $(\cdot,\cdot)_T$ is evaluated over element $T$. This method was introduced by Hughes, Franca, and Balestra in 1986 and is a residual-based method. It requires second derivatives on elements (if applicable; they are only zero for lowest order spaces with affine elements). When the solution is not very well resolved, the element integral can be very different from the actual cell residual (because the influence from neighbor cells is discarded). The stabilization produces a non-physical boundary condition for the pressure, leading to local inaccuracies impacting functionals of interest (such as drag). There are more robust and accurate edge-stabilized modifications, but these are more complicated and grow the stencil. Residual stabilization adds a large number of couplings for more complicated problems (like reaction and turbulence models) and significant implementation complexity. The formulation is cumbersome for time-dependent problems because the "affine" term enters in an inner product, thus typically precluding a method of lines approach. PSPG also breaks symmetry (relevant for adjoint-consistency) and increases regularity requirements for an adjoint, making it less attractive for optimization and sensitivity analysis.

• Local pressure projection This conceptually simpler method uses a stabilization of the form $$ S_h^{\text{lps}} = (\alpha \tilde \pi p_h, \tilde \pi q_h)$$ where $\tilde \pi$ is a local fluctuation operator and $\alpha$ is a scaling factor. Dohrmann and Bochev (2004) studied what is perhaps the simplest choice, $\tilde \pi = 1 - \pi_{k-1}$, where $\pi_{k-1}$ is cell-wise projection into a polynomial space of degree one less than the continuous space. There are many patch-based and edge-stabilized choices, see Braack and Lube's 2009 review. These methods are simpler to implement and preserve symmetry, but are often less robust.