Streamline upwind Petrov–Galerkin pressure-stabilizing Petrov–Galerkin formulation for incompressible Navier–Stokes equations explained

The streamline upwind Petrov–Galerkin pressure-stabilizing Petrov–Galerkin formulation for incompressible Navier–Stokes equations can be used for finite element computations of high Reynolds number incompressible flow using equal order of finite element space (i.e.

P _{k-P} _{k} ) by introducing additional stabilization terms in the Navier–Stokes Galerkin formulation.^{[1] [2]}

The finite element (FE) numerical computation of incompressible Navier–Stokes equations (NS) suffers from two main sources of numerical instabilities arising from the associated Galerkin problem.^[1] Equal order finite elements for pressure and velocity, (for example,

P _{k-P} _{k,} \forall k \ge 0 ), do not satisfy the inf-sup condition and leads to instability on the discrete pressure (also called spurious pressure).^[2] Moreover, the advection term in the Navier–Stokes equations can produce oscillations in the velocity field (also called spurious velocity).^[2] Such spurious velocity oscillations become more evident for advection-dominated (i.e., high Reynolds number Re ) flows.^[2] To control instabilities arising from inf-sup condition and convection dominated problem, pressure-stabilizing Petrov–Galerkin (PSPG) stabilization along with Streamline-Upwind Petrov-Galerkin (SUPG) stabilization can be added to the NS Galerkin formulation.^{[1] [2]}

The incompressible Navier–Stokes equations for a Newtonian fluid

Let

\Omega \subset R ^{3} be the spatial fluid domain with a smooth boundary \partial \Omega \equiv \Gamma , where \Gamma = \Gamma _{N} \cup \Gamma _{D} with \Gamma _{D} the subset of \Gamma in which the essential (Dirichlet) boundary conditions are set, while \Gamma _{N} the portion of the boundary where natural (Neumann) boundary conditions have been considered. Moreover, \Gamma _{N} = \Gamma \setminus \Gamma _{D} , and \Gamma _{N} \cap \Gamma _{D=\emptyset} . Introducing an unknown velocity field u ( x ,t):\Omega x [ 0 ,T] → R ^{3} and an unknown pressure field p( x ,t):\Omega x [ 0 ,T] → R , in absence of body forces, the incompressible Navier–Stokes (NS) equations read^[3] where \hat{n } is the outward directed unit normal vector to \Gamma _{N} , \boldsymbol{\sigma} is the Cauchy stress tensor, \rho is the fluid density, and \nabla and \nabla ⋅ are the usual gradient and divergence operators.The functions g and h indicate suitable Dirichlet and Neumann data, respectively, while u _{0} is the known initial field solution at time t =0 .

For a Newtonian fluid, the Cauchy stress tensor

\boldsymbol{\sigma} depends linearly on the components of the strain rate tensor:^[3] $$\backslash boldsymbol\; (\backslash mathbf\; u,p)=-p\; \backslash mathbf\; +2\backslash mu\; \backslash mathbf\; S(\backslash mathbf\; u),$$where \mu is the dynamic viscosity of the fluid (taken to be a known constant) and I is the second order identity tensor, while S( u) is the strain rate tensor$$\backslash mathbf\; S(\backslash mathbf\; u)=\backslash frac\; \backslash big[\backslash nabla\; \backslash mathbf\; u\; +\; (\backslash nabla\; \backslash mathbf\; u)^T\; \backslash big].$$u _{0} , g , and h are assigned.

Hence, the strong formulation of the incompressible Navier–Stokes equations for a constant density, Newtonian and homogeneous fluid can be written as:^[3]

Find,

\forall t \in ( 0 ,T] , velocity u( x ,t) and pressure p( x ,t) such that:where, \nu = (\mu)/(\rho) is the kinematic viscosity, and \hat p = (p)/(\rho) is the pressure rescaled by density (however, for the sake of clearness, the hat on pressure variable will be neglect in what follows).

In the NS equations, the Reynolds number shows how important is the non linear term,

( u ⋅ \nabla ) u , compared to the dissipative term, \nu \nabla ⋅ S( u): ^[4] $$\backslash frac\backslash approx\; \backslash frac=\backslash frac\; =\; \backslash mathrm.$$

The Reynolds number is a measure of the ratio between the advection convection terms, generated by inertial forces in the flow velocity, and the diffusion term specific of fluid viscous forces.^[4] Thus,

Re can be used to discriminate between advection-convection dominated flow and diffusion dominated one.^[4] Namely:

• for "low"

Re , viscous forces dominate and we are in the viscous fluid situation (also named Laminar Flow),^[4]

• for "high"

Re , inertial forces prevail and a slightly viscous fluid with high velocity emerges (also named Turbulent Flow).^[4]

The weak formulation of the Navier–Stokes equations

v and q , respectively, belonging to suitable function spaces, and integrating these equation all over the fluid domain \Omega .^[3] As a consequence:^[3]

By summing up the two equations and performing integration by parts for pressure (

\nabla p ) and viscous (\nabla ⋅ S ( u) ) term:^[3] $$\backslash int\_\backslash Omega\; \backslash frac\backslash cdot\; \backslash mathbf\; v\backslash ,d\backslash Omega+\; \backslash int\_(\backslash mathbf\; u\; \backslash cdot\; \backslash nabla)\; \backslash mathbf\; u\; \backslash cdot\; \backslash mathbf\; v\; \backslash ,d\backslash Omega\backslash ,+\backslash int\_\backslash Omega\; \backslash nabla\; \backslash cdot\; \backslash mathbf\; u\; \backslash ,\; q\; \backslash ,d\backslash Omega-\; \backslash int\_p\; \backslash nabla\; \backslash cdot\; \backslash mathbf\; v\; \backslash ,d\backslash Omega+\backslash int\_p\; \backslash mathbf\; v\; \backslash cdot\; \backslash mathbf\; \backslash ,d\backslash Gamma\; \backslash ,+\; \backslash int\_\backslash Omega\; 2\backslash nu\; \backslash mathbf\; S(\backslash mathbf\; u)\; :\; \backslash nabla\; \backslash mathbf\; v\; \backslash ,d\backslash Omega-\backslash int\_2\backslash nu\; \backslash mathbf\; S(\backslash mathbf\; u)\; \backslash cdot\; \backslash mathbf\; v\; \backslash cdot\; \backslash mathbf\; \backslash ,d\backslash Gamma\; \backslash ,\; =0.$$

Regarding the choice of the function spaces, it's enough that

p and q , u and v , and their derivative, \nabla u and \nabla v are square-integrable functions in order to have sense in the integrals that appear in the above formulation. Hence,

Having specified the function spaces

l{V} , l{V} _{0} and l{Q} , and by applying the boundary conditions, the boundary terms can be rewritten as$$\backslash int\_p\; \backslash mathbf\; v\; \backslash cdot\; \backslash mathbf\; \backslash ,d\backslash Gamma+\; \backslash int\_\; -2\backslash nu\; S(\backslash mathbf\; u)\; \backslash cdot\; \backslash mathbf\; v\; \backslash cdot\; \backslash mathbf\; \backslash ,d\backslash Gamma,$$where \partial \Omega = \Gamma _{D} \cup \Gamma _{N} . The integral terms with \Gamma _{D} vanish because v| _{\Gamma _{D}} = 0 , while the term on \Gamma _{N} become$$\backslash int\_\; [p\; \backslash mathbf\; I\; -2\backslash nu\; S(\backslash mathbf\; u)]\; \backslash cdot\; \backslash mathbf\; v\; \backslash cdot\; \backslash mathbf\; \backslash ,d\backslash Gamma\; =\; -\backslash int\_\; \backslash mathbf\; h\; \backslash cdot\; \backslash mathbf\; v\; \backslash ,d\backslash Gamma,$$

The weak formulation of Navier–Stokes equations reads:^[3]

Find, for all

t \in ( 0 ,T] , ( u,p)\in \{ l{V} x l{Q}\} , such that$$\backslash left(\backslash frac,\backslash mathbf\; v\; \backslash right)+c(\backslash mathbf\; u,\backslash mathbf\; u,\backslash mathbf\; v)+b(\backslash mathbf\; u,q)-b(\backslash mathbf\; v,p)\; +\; a(\backslash mathbf\; u,\backslash mathbf\; v)\; =\; f(\backslash mathbf\; v)$$

with

u | _{t=0} = u _{0} , where

Finite element Galerkin formulation of Navier–Stokes equations

\Omega _{h} , composed by tetrahedra l{T} _{i} , with i = 1 , \ldots, N _{l{T} } (where N _{l{T} } is the total number of tetrahedra), of the domain \Omega and h is the characteristic length of the element of the triangulation.^[3] l{V} _{h} and l{Q} _{h} , approximations of l{V} and l{Q} respectively, and depending on a discretization parameter h , with \dim l{V} _{h} = N _{V} and \dim l{Q} _{h} = N _{Q} ,^[3] $$\backslash mathcal\_h\; \backslash subset\; \backslash mathcal\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; \backslash mathcal\_h\; \backslash subset\; \backslash mathcal,$$the discretized-in-space Galerkin problem of the weak NS equation reads:^[3]

Find, for all

t \in ( 0 ,T] , ( u _{h,p} _{h)\in} \{ l{V} _{h} x l{Q} _{h\}} , such thatwith u _{h|} _{t=0} = u _{h,0} , where g _{h} is the approximation (for example, its interpolant) of g , and$$\backslash mathcal\_=\backslash .$$[ 0 ,T] into N _{t} time step of size \delta t ^[3] $$t\_n=n\backslash delta\; t,\; \backslash ;\backslash ;\backslash ;\; n=0,1,2,\backslash ldots,N\_t\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; N\_t=\backslash frac.$$

For a general function

z , denoted by z ^{n} as the approximation of z(t _{n)} . Thus, the BDF2 approximation of the time derivative reads as follows:^[3] $$\backslash left(\backslash frac\; \backslash right)^\; \backslash simeq\; \backslash frac\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; \backslash text\; n\; \backslash geq\; 1.$$

So, the fully discretized in time and space NS Galerkin problem is:^[3]

Find, for

n = 0 , 1 , \ldots, N _{t-1} , ( u ^{n +1} _{h)\in} \{ l{V} _{h} x l{Q} _{h\}} , such thatwith u ^{0} _{h} = u _{h,0} , and u ^{*} _{h} is a quantity that will be detailed later in this section.

The main issue of a fully implicit method for the NS Galerkin formulation is that the resulting problem is still non linear, due to the convective term,

c( u ^{n +1} _{h, v} _{h)} .^[3] Indeed, if u ^{n +1} _{h} is put, this choice leads to solve a non-linear system (for example, by means of the Newton or Fixed point algorithm) with a huge computational cost.^[3] In order to reduce this cost, it is possible to use a semi-implicit approach with a second order extrapolation for the velocity, u ^{*} _{h} , in the convective term:^[3] $$\backslash mathbf\; u^*\_h=2\backslash mathbf\; u^\_h-\backslash mathbf\; u^\_h.$$

Finite element formulation and the INF-SUP condition

Let's define the finite element (FE) spaces of continuous functions,

X ^{r} _{h} (polynomials of degree r on each element l{T} _{i} of the triangulation) as^[3] $$X\_h^r=\; \backslash left\backslash \; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; r=0,1,2,\backslash ldots,$$where, P _{r} is the space of polynomials of degree less than or equal to r .

Introduce the finite element formulation, as a specific Galerkin problem, and choose

l{V} _{h} and l{Q} _{h} as$$\backslash mathcal\_h\backslash equiv\; [X\_h^r]^3\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; \backslash mathcal\_h\backslash equiv\; X\_h^s\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; r,s\; \backslash in\; \backslash mathbb.$$

The FE spaces

l{V} _{h} and l{Q} _{h} need to satisfy the inf-sup condition(or LBB):^[6] $$\backslash exists\; \backslash beta\_h\; >0\; \backslash ;\backslash text\; \backslash ;\; \backslash inf\_\backslash sup\_\; \backslash frac\; \backslash geq\; \backslash beta\_h\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; \backslash forall\; h>0,$$

with

\beta _{h} > 0 , and independent of the mesh size h. This property is necessary for the well posedness of the discrete problem and the optimal convergence of the method. Examples of FE spaces satisfying the inf-sup condition are the so named Taylor-Hood pair P _{k+1} - P _{k} (with k \geq 1 ), where it can be noticed that the velocity space l{V} _{h} has to be, in some sense, "richer" in comparison to the pressure space l{Q} _{h.} Indeed, the inf-sup condition couples the space l{V} _{h} and l{Q} _{h} , and it is a sort of compatibility condition between the velocity and pressure spaces.

The equal order finite elements,

P _{k} - P _{k} (\forall k ), do not satisfy the inf-sup condition and leads to instability on the discrete pressure (also called spurious pressure). However, P _{k} - P _{k} can still be used with additional stabilization terms such as Streamline Upwind Petrov-Galerkin with a Pressure-Stabilizing Petrov-Galerkin term (SUPG-PSPG).

In order to derive the FE algebraic formulation of the fully discretized Galerkin NS problem, it is necessary to introduce two basis for the discrete spaces

l{V} _{h} and l{Q} _{h} $$\backslash \_^\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; \backslash \_^,$$in order to expand our variables as^[3] $$\backslash mathbf\; u^n\_h\; =\; \backslash sum\_^\; U^n\_j\; \backslash boldsymbol\_j(\backslash mathbf\; x),\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; q^n\_h=\backslash sum\_^\; P^n\_l\; \backslash psi\_l(\backslash mathbf\; x).$$

The coefficients,

U ^{n} _{j} (j =1 ,\ldots,N _{V} ) and P ^{n} _{l} (l =1 ,\ldots,N _{Q} ) are called degrees of freedom (d.o.f.) of the finite element for the velocity and pressure field, respectively. The dimension of the FE spaces, N _{V} and N _{Q} , is the number of d.o.f, of the velocity and pressure field, respectively. Hence, the total number of d.o.f N _{d.o.f} is N _{d.o.f} = N _{V+N} _{Q} .^[3]

Since the fully discretized Galerkin problem holds for all elements of the space

l{V} _{h} and l{Q} _{h} , then it is valid also for the basis.^[3] Hence, choosing these basis functions as test functions in the fully discretized NS Galerkin problem, and using bilinearity of a( ⋅ , ⋅ ) and b( ⋅ , ⋅ ) , and trilinearity of c( ⋅ , ⋅ , ⋅ ) , the following linear system is obtained:^[3] $$\backslash begin\backslash displaystyle\; M\; \backslash frac\; +\; A\backslash mathbf\; U^\; +C(\backslash mathbf\; U^*)\backslash mathbf\; U^+\backslash displaystyle\backslash \backslash \backslash displaystyle\backslash end$$where M \in R ^{N _{V} x N _{V}} , A \in R ^{N _{V} x N _{V}} , C( U ^{*)} \in R ^{N _{V} x N _{V}} , B \in R ^{N _{Q} x N _{V}} , and F \in R ^{N _{V}} are given by^[3] and U and P are the unknown vectors^[3] $$\backslash mathbf\; U^n=\backslash Big(U^n\_1,\backslash ldots,U^n\_\; \backslash Big)^T,\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; \backslash mathbf\; P^n=\backslash Big(P^n\_1,\backslash ldots,P^n\_\; \backslash Big)^T.$$

Problem is completed by an initial condition on the velocity

U(0)= U _{0} . Moreover, using the semi-implicit treatment U ^{*} =2 U ^{n} - U ^{n-1} , the trilinear term c( ⋅ , ⋅ , ⋅ ) becomes bilinear, and the corresponding matrix is^[3] $$C\_=c(\backslash mathbf\; u^*,\backslash boldsymbol\_j,\backslash boldsymbol\_i)=\backslash int\_(\backslash mathbf\; u^*\; \backslash cdot\; \backslash nabla)\; \backslash boldsymbol\_j\; \backslash cdot\; \backslash boldsymbol\_i\; \backslash ,d\backslash Omega,$$

Hence, the linear system can be written in a single monolithic matrix (

\Sigma , also called monolithic NS matrix) of the form^[3] where $K=\backslash fracM+A+C(U^*)$.

Streamline upwind Petrov–Galerkin formulation for incompressible Navier–Stokes equations

NS equations with finite element formulation suffer from two source of numerical instability, due to the fact that:

• NS is a convection dominated problem, which means "large"

Re , where numerical oscillations in the velocity field can occur (spurious velocity);

• FE spaces

P _{k-P} _{k} (\forall k) are unstable combinations of velocity and pressure finite element spaces, that do not satisfy the inf-sup condition, and generates numerical oscillations in the pressure field (spurious pressure).

To control instabilities arising from inf-sup condition and convection dominated problem, Pressure-Stabilizing Petrov–Galerkin(PSPG) stabilization along with Streamline-Upwind Petrov–Galerkin (SUPG) stabilization can be added to the NS Galerkin formulation.^[1]

$$s(\backslash mathbf\; u^\_h,\; p^\_h\; ;\backslash mathbf\; v\_h,\; q\_h)=\backslash gamma\; \backslash sum\_\; \backslash tau\_\; \backslash int\_\backslash left[\backslash mathcal\{L\}(\backslash mathbf\; u^\{n+1\}\_h,\; p^\{n+1\})\; \backslash right]^T\; \backslash mathcal\_(\backslash mathbf\; v\_h,\; q\_h)d\backslash mathcal,$$where

\gamma> 0 is a positive constant, \tau _{l{T} } is a stabilization parameter, l{T} is a generic tetrahedron belonging to the finite elements partitioned domain \Omega _{h} , l{L}( u, p) is the residual of the NS equations.^[1]

$$\backslash mathcal(\backslash mathbf\; u,\; p)\; =\; \backslash begin\backslash frac+\; (\backslash mathbf\; u\; \backslash cdot\; \backslash nabla)\; \backslash mathbf\; u\; +\; \backslash nabla\; p\; -2\backslash nu\; \backslash nabla\; \backslash cdot\; \backslash mathbf\; S(\backslash mathbf\; u)\; \backslash \backslash \; \backslash nabla\; \backslash cdot\; \backslash mathbf\; u\; \backslash end,$$and

l{L} _{ss} ( u, p) is the skew-symmetric part of the NS equations^[1] $$\backslash mathcal\_(\backslash mathbf\; u,\; p)\; =\; \backslash begin(\backslash mathbf\; u\; \backslash cdot\; \backslash nabla)\; \backslash mathbf\; u\; +\; \backslash nabla\; p\; \backslash \backslash \; \backslash mathbf\; 0\; \backslash end\; .$$

The skew-symmetric part of a generic operator

l{L}( u, p) is the one for which l ( l{L}( u, p),( v, q) r ) = - l ( ( v, q), l{L}( u, p) r ). ^[5]

Since it is based on the residual of the NS equations, the SUPG-PSPG is a strongly consistent stabilization method.^[1]

The discretized finite element Galerkin formulation with SUPG-PSPG stabilization can be written as:^[1]

Find, for all

t = 0 , 1 , \ldots, N _{t-1,} ( u ^{n +1} _{h)\in} \{ l{V} _{h} x l{Q} _{h\}} , such thatwith u ^{0} _{h = u} _{h,0} , where^[1]

and

\tau _{M,l{T} } , and \tau _{C,l{T} } are two stabilization parameters for the momentum and the continuity NS equations, respectively. In addition, the notation $\backslash left(a\; \backslash boldsymbol\; \backslash ;\; b\; \backslash right)\_=\backslash int\_ab\; \backslash ;\; d\backslash mathcal$ has been introduced, and u ^{*} _{h} was defined in agreement with the semi-implicit treatment of the convective term.^[1]

In the previous expression of

s \left( ⋅ ; ⋅ \right) , the term $\backslash sum\_\; \backslash tau\_\; \backslash left(\backslash nabla\; p^\_\; \backslash boldsymbol\; \backslash ;\; \backslash frac\; \backslash right)\_,$ is the Brezzi-Pitkaranta stabilization for the inf-sup, while the term$\backslash sum\_\; \backslash tau\_\; \backslash left(u\_h^*\; \backslash cdot\; \backslash nabla\; \backslash mathbf\; u^\_h\; \backslash boldsymbol\; \backslash ;\; u\_h^*\; \backslash cdot\; \backslash nabla\; \backslash mathbf\; v\_h\; \backslash right)\_,$corresponds to the streamline diffusion term stabilization for large Re .^[1] The other terms occur to obtain a strongly consistent stabilization.^[1]

Regarding the choice of the stabilization parameters

\tau _{M,l{T} } , and \tau _{C,l{T} } :^[2] $$\backslash tau\_=\backslash left(\backslash frac\; +\; \backslash frac\; +\; C\_k\backslash frac\; \backslash right)^,\; \backslash ;\backslash ;\backslash ;\backslash ;\backslash ;\; \backslash tau\_=\backslash frac,$$

where:

C _{k=60} ⋅ 2 ^{k-2} is a constant obtained by an inverse inequality relation (and k is the order of the chosen pair P _{k-P} _{k} ); \sigma _{BDF} is a constant equal to the order of the time discretization; \delta t is the time step; h _{l{T} } is the "element length" (e.g. the element diameter) of a generic tetrahedra belonging to the partitioned domain \Omega _{h} .^[7] The parameters \tau _{M,l{T} } and \tau _{C,l{T} } can be obtained by a multidimensional generalization of the optimal value introduced in^[8] for the one-dimensional case.^[9]

Notice that the terms added by the SUPG-PSPG stabilization can be explicitly written as follows^[2]

where, for the sake of clearness, the sum over the tetrahedra was omitted: all the terms to be intended as $s^\_\; =\; \backslash sum\_\; \backslash tau\_\backslash left(\backslash cdot\; \backslash ,,\; \backslash cdot\; \backslash right)\_\backslash mathcal$; moreover, the indices

I,J in s ^{(n)} _{(I,J)} refer to the position of the corresponding term in the monolithic NS matrix, \Sigma , and n distinguishes the different terms inside each block^[2]

Hence, the NS monolithic system with the SUPG-PSPG stabilization becomes^[2] where $\backslash tilde=K+\backslash sum\backslash limits\_^4\; S^\_$, and $\backslash tilde=B+\backslash sum\backslash limits\_^3S^\_$.

It is well known that SUPG-PSPG stabilization does not exhibit excessive numerical diffusion if at least second-order velocity elements and first-order pressure elements (

P _{2-P} _{1} ) are used.^[8]

References

1. Tayfun Tezduyar . Tezduyar . T. E. . Stabilized Finite Element Formulations for Incompressible Flow Computations††This research was sponsored by NASA-Johnson Space Center (under grant NAG 9-449), NSF (under grant MSM-8796352), U.S. Army (under contract DAAL03-89-C-0038), and the University of Paris VI. . Advances in Applied Mechanics . 1 January 1991 . 28 . 1–44 . Elsevier. 10.1016/S0065-2156(08)70153-4 .
2. Tobiska . Lutz . Lube . Gert . A modified streamline diffusion method for solving the stationary Navier–Stokes equation . Numerische Mathematik . 1 December 1991 . 59 . 1 . 13–29 . 10.1007/BF01385768 . 123397636 . en . 0945-3245.
3. Book: Alfio Quarteroni . Quarteroni . Alfio . Numerical Models for Differential Problems . 2014 . Springer-Verlag . 9788847058835 . 2 . en.
4. Book: Pope . Stephen B. . Turbulent Flows by Stephen B. Pope . 2000 . Cambridge University Press . 9780521598866 . en.
5. Book: Alfio Quarteroni . Quarteroni . Alfio . Sacco . Riccardo . Saleri . Fausto . Numerical Mathematics . 2007 . Springer-Verlag . 9783540346586 . 2 . en.
6. Book: Mixed and Hybrid Finite Element Methods . 15 . 1991 . en-gb. 10.1007/978-1-4612-3172-1 . Springer Series in Computational Mathematics . 978-1-4612-7824-5 . Brezzi . Franco . Fortin . Michel .
7. Forti . Davide . Dedè . Luca . Semi-implicit BDF time discretization of the Navier–Stokes equations with VMS-LES modeling in a High Performance Computing framework . Computers & Fluids . August 2015 . 117 . 168–182 . 10.1016/j.compfluid.2015.05.011.
8. Shih . Rompin . Ray . S. E. . Mittal . Sanjay . Tezduyar . T. E. . Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements . Computer Methods in Applied Mechanics and Engineering . 95 . 2 . 221 . 1992 . en. 1992CMAME..95..221T . 10.1016/0045-7825(92)90141-6 . 31236394 .
9. Kler . Pablo A. . Dalcin . Lisandro D. . Paz . Rodrigo R. . Tezduyar . Tayfun E. . SUPG and discontinuity-capturing methods for coupled fluid mechanics and electrochemical transport problems . Computational Mechanics . 1 February 2013 . 51 . 2 . 171–185 . 10.1007/s00466-012-0712-z . en . 1432-0924. 2013CompM..51..171K . 123650035 . 11336/1065 . free .