Navier-Stokes equations

$\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \newcommand{\bvec}[1]{\mathbf{#1}} \newcommand{\bx}{\mathbf{x}} \newcommand{\bmv}{\mathbf{v}} \newcommand{\bmw}{\mathbf{w}} \newcommand{\Nh}[1][]{N_{h,#1}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bp}{\bvec{p}} \newcommand{\VV}{\mathbf{V}} \newcommand{\WW}{\mathbf{W}} \newcommand{\QQ}{\mathbf{Q}} \renewcommand{\AA}{\mathbf{A}} \newcommand{\BB}{\mathbf{B}} \newcommand{\FF}{\mathbf{F}} \newcommand{\GG}{\mathbf{G}} \newcommand{\CC}{\mathbf{C}} \newcommand{\MM}{\mathbf{M}} \newcommand{\XX}{\mathbf{X}} \newcommand{\KK}{\mathbf{K}} \newcommand{\UU}{\mathbf{U}} \newcommand{\II}{\mathbf{I}} \renewcommand{\SS}{\mathbf{S}} \newcommand{\HH}{\mathbf{H}} \newcommand{\DD}{\mathbf{D}} \newcommand{\NN}{\mathbf{N}} \newcommand{\bmg}{\mathbf{g}} \newcommand{\R}{\mathbb{R}} \newcommand{\JJ}{\mathbf{J}} \newcommand{\RR}{\mathbf{R}}$

Unsteady Navier-Stokes equations
Inf-Sup stable finite element approximation
- Fully implicit BDF time discretization
  - Algebraic formulation and Newton iterations
- Semi-implicit BDF time discretization
  - Algebraic formulation
SUPG-stabilized finite element approximation
- Semi-implicit BDF time discretization
  - Algebraic formulation
- Fully implicit BDF time discretization
  - Algebraic formulation and Newton iterations
References

Unsteady Navier-Stokes equations

Let $\Omega \subset \mathbb{R}^{d}$ ($d=2,3$) be an open bounded domain with piecewise smooth boundary $\Gamma = \partial \Omega$. The latter is decomposed into Dirichlet and Neumann components such that $\Gamma = \Gamma_D \cup \Gamma_N$. The Navier-Stokes equations for an incompressible, homogeneous, Newtonian fluid read:

$\begin{equation} \left\{ \begin{aligned} \rho \frac{\partial \bmv}{\partial t} + \rho \bmv \cdot \nabla \bmv - \nabla \cdot \bm{\sigma}(\bmv, p) &= \bm{0} & \qquad \quad \text{in } & \Omega \times (0,T) \\ \nabla \cdot \bmv &= 0 & \qquad \quad \text{in } &\Omega\times (0,T)\\ \bmv &= \bm{g} & \text{on } &\Gamma_{D} \times (0,T)\\ \bm{\sigma}(\bmv, p) \bm{n} &= \bvec{h} & \text{on } &\Gamma_{N}\times (0,T) \\ \bmv(0) &= \bmv_0 & \text{in } &\Omega, \end{aligned} \right. \label{eq:NSt_cont_form} \end{equation}$

where $\bmv=\bmv(\mathbf{x},t)$ is the fluid velocity, $p=p(\bvec{x},t)$ the pressure, $\bm{n}$ the (outward directed) normal unit vector to $\Gamma_N$, $\rho$ the fluid density and $\bm{\sigma}$ is the stress tensor defined as \begin{equation} \label{eq:NSt_stress_tensor} \bm{\sigma}(\bmv, p) = -p \mathbf{I} + 2 \mu \bm{\epsilon}(\bmv). \end{equation} Here $\mu$ denotes the dynamic viscosity of the fluid (the kinematic viscosity is given by $\nu = \mu / \rho$), while \begin{equation} \label{eq:NSt_strain_tensor} \bm{\epsilon}(\bmv) = \frac{1}{2} \big( \nabla \bmv + \nabla \bmv ^T \big) \end{equation} is the strain tensor. The functions $\bm{g} = \bm{g}(\bx, t)$, $\bm{h} = \bm{h}(\bx, t)$ indicate the Dirichlet and Neumann data, while $\bmv_0 = \bmv_0(\bx)$ is the initial condition.

Let us introduce the following functional spaces: $V_D = \{ \bmw \in [H^1(\Omega)]^d : \bmw|_{\Gamma_D} = \bm{h}\}$ $V = \{ \bmw \in [H^1(\Omega)]^d : \bmw|_{\Gamma_D} = \mathbf{0} \}$ and $M = L^2(\Omega)$. The weak formulation of \eqref{eq:NSt_cont_form} reads: for all $t \in (0,T]$, find $(\bmv(t), p(t)) \in V_D \times M$ such that \begin{equation} \label{eq:NSt_weak_form} \left( \rho \frac{\partial \bmv}{\partial t}, \bmw \right) + (\rho \bmv \cdot \nabla \bmv ,\bmw) + (\mu (\nabla \bmv + \nabla \bmv^T),\nabla \bmw) - (p, \nabla \cdot \bmw) + (\nabla \cdot \bmv,q) = (\bm{h}, \bmw)_{\Gamma_N} \end{equation} for all $ (\bmw,q) \in V \times M$, with $\bmv(0) = \bmv_0$.

Inf-Sup stable finite element approximation

Let us introduce a FE partition $\mathcal{T}_h$ of the domain $\Omega$ from which we construct conforming finite element spaces $V_h \subset V$ and $M_h \subset M$. For the discrete version of problem \eqref{eq:NSt_weak_form} to be well-posed, it is well known that the velocity and pressure spaces $V_h$ and $M_h$ need to fulfill an inf-sup condition \begin{equation} \label{eq:NSt_infsup} \inf_{q_h \in M_h} \sup_{\bmw_h \in V_h} \frac{(q_h, \nabla \cdot \bmw_h)}{|\bmw_h|_V |q_h|_M} \geq \bar{\beta} > 0. \end{equation}

Two different choices of stable finite element spaces are implemented in redbKIT: $\mathbb{P}_2$-$\mathbb{P}_1$ and $\mathbb{P}_1^{\text{bubble}}$-$\mathbb{P}_1$ velocity-pressure spaces.

The semi-discrete formulation of the Navier-Stokes equations reads: for all $t \in (0,T]$, find $(\bmv_h, p_h) \in V_{D,h} \times M_h$ such that

$\begin{multline} \label{eq:NSt_weak_STABLEh} \left( \rho \frac{\partial \bmv_h}{\partial t}, \bmw_h \right) + ( \rho \bmv_h \cdot \nabla \bmv_h ,\bmw_h) + (\mu (\nabla \bmv_h + \nabla \bmv^T_h),\nabla \bmw_h) - (p_h, \nabla \cdot \bmw_h) \\ + (\nabla \cdot \bmv_h,q_h) = (\bm{h}, \bmw)_{\Gamma_N} \end{multline}$

for all $(\bmw_h,q_h) \in V_h \times M_h$, with $\bmv_h(0) = \bmv_0$.

Fully implicit BDF time discretization

For the time discretization of \eqref{eq:NSt_weak_STABLEh} we consider the Backward Differentiation Formula (BDF) scheme. To begin with, we partition the time interval $[0,T]$ into $N_t$ subintervals of equal size $\Delta t = T / N_t$ and denote by $t_n = n \Delta t$, for $n = 0 , \dots, N_t$ the discrete time instances. Moreover, we denote by $\bmv_h^n$ and $p_h^n$ the approximations of $\bmv_h$ and $p_h$ at time $t_n$, respectively. We approximate the time derivative of the velocity as

$\frac{\partial \bmv_h }{\partial t} \approx \frac{\alpha \bmv_h^{n+1} - \bmv_h^{n,\sigma}}{\Delta t}.$

where (limiting ourselves to BDF schemes of order $\sigma = 1,2$)

$\bmv_h^{n,\sigma} = \begin{cases} \bmv_h^n, & \text{ if } n \geq 0, \quad \text{for } \sigma =1 \text{ (BDF1)}, \\ 2\bmv_h^n - \tfrac{1}{2}\bmv_h^{n-1}, & \text{ if } n \geq 1, \quad \text{for } \sigma =2 \text{ (BDF2)}, \end{cases}$

and

$\alpha = \begin{cases} 1, & \text{ for } \sigma =1 \text{ (BDF1)}, \\ \tfrac{3}{2}, & \text{ for } \sigma =2 \text{ (BDF2)}. \end{cases}$

The (fully discrete, implicit) BDF approximation of the Navier-Stokes equations reads: given $\bmv_h^{n}, \dots, \bmv_h^{n+1-\sigma}$, for $n \geq \sigma-1$ find $(\bmv_h^{n+1}, p_h^{n+1}) \in V_{h} \times M_h$ such that

$\begin{multline} \label{eq:NSt_BDF_STABLE} \left( \rho \frac{ \alpha \bmv_h^{n+1} - \bmv_h^{n,\sigma}}{\Delta t}, \bmw_h \right) + \left(\rho \bmv_h^{n+1} \cdot \nabla \bmv_h^{n+1} ,\bmw_h\right) + \left(\mu (\nabla \bmv_h^{n+1} + (\nabla \bmv^{n+1}_h)^T),\nabla \bmw_h\right) \\ - \left(p_h^{n+1}, \nabla \cdot \bmw_h\right) + \left(\nabla \cdot \bmv_h^{n+1},q_h\right) = g( (\bmw_h, q_h); t_{n+1} ) + (\bm{h}^{n+1}, \bmw_h)_{\Gamma_N} \end{multline}$ for all $(\bmw_h,q_h) \in V_h \times M_h$. Here, the functional $g(\cdot;t) : V_h \times M_h \rightarrow \mathbb{R}$ encodes the action of the nonhomogeneous Dirichlet condition $\bmv_h|_{\Gamma_D} = \bm{h}(t)$.

Algebraic formulation and Newton iterations

We denote by $\{\bm{\varphi}_i\}_{i=1}^{\Nh[v]}$ and $\{\eta_k\}_{k=1}^{\Nh[p]}$ Lagrangian FE bases for $V_h$ and $M_h$ respectively. We also denote by $\bmv^n_h \in \mathbb{R}^{\Nh[v]}$ and $\bp^n_h \in \mathbb{R}^{\Nh[p]}$ the vectors of coefficients in the expansions of $\bm{v}^n_h$ and $p_h^n$ with respect to the FE bases. Finally, we set $\bU_h^n = (\bmv_h^n, \bp^n_h) \in \mathbb{R}^{\Nh}$ . The algebraic formulation of \eqref{eq:NSt_BDF_STABLE} reads: given $\bU_h^{n}, \dots, \bU_h^{n+1-\sigma}$ , for $n \geq \sigma-1$ find $\bU^{n+1}_h \in \mathbb{R}^{\Nh}$ such that

$\begin{equation} \frac{ \rho \alpha }{\Delta t}\MM \bU^{n+1}_h + \AA \bU^{n+1}_h + \mathbf{c}( \bU^{n+1}_h) = \frac{\rho}{\Delta t} \MM\bU_h^{n,\sigma} + \bm{f}(t_{n+1}). \label{eq:NSt_BDF_STABLE_ALGtierI} \end{equation}$

The right-hand side term $\bm{f}(t_{n+1})$ encodes the action of Dirichlet and Neumann boundary conditions (as well as body forces, if any). The $\Nh \times \Nh$ block-partitioned matrices $\MM$ and $\AA$ are defined as

$\begin{equation} \MM = \begin{pmatrix} \MM_v & 0 \\ 0 & 0 \end{pmatrix}, \qquad \AA = \begin{pmatrix} \mu \KK & \BB^T \\ -\BB & 0 \end{pmatrix}, \end{equation}$

where $\MM_v \in \R^{\Nh[v] \times \Nh[v]}$ is the velocity mass matrix, $\KK \in \R^{\Nh[v] \times \Nh[v]}$ is the velocity stiffness matrix, $\BB \in \R^{\Nh[v] \times \Nh[v]}$ is the divergence matrix. The nonlinear convective term $\mathbf{c}(\bU^{n+1}_h)$ is defined as

where $\widetilde{\mathbf{c}}(\bmv_h^{n+1})$ is given by

$\begin{equation} (\widetilde{\mathbf{c}}(\bmv_h^{n+1}))_{i} = \left( \rho \bmv_h^{n+1} \cdot \nabla \bmv_h^{n+1}, \bm{\varphi}_i\right) \qquad i = 1, \dots, \Nh[v]. \end{equation}$

At each timestep, the nonlinear system \eqref{eq:NSt_BDF_STABLE_ALGtierI} is solved by means of Newton’s method: given $\bU_h^{n+1,0} = \bU_h^{n}$, for $k=0,1,\dots$ until convergence we seek $\delta \bU$ such that

$\begin{equation} \JJ(\bU_h^{n+1,k}) \, \delta \bU = - \RR(\bU_h^{n+1,k}), \end{equation}$

and then set $\bU_h^{n+1,k+1} = \bU_h^{n+1,k} + \delta \bU$. Here,

$\begin{equation} \JJ = \frac{ \rho \alpha }{\Delta t}\MM + \AA + \CC_1( \bU^{n+1,k}_h) + \CC_2( \bU^{n+1,k}_h) \end{equation}$

and

$\begin{equation} \RR = \frac{ \rho }{\Delta t}\MM \left( \alpha \bU^{n+1,k}_h - \bU^{n,\sigma}_h \right) + \AA \bU^{n+1,k}_h + \mathbf{c}( \bU^{n+1,k}_h) - \bm{f}(t_{n+1}). \end{equation}$

The convective matrices $\CC_1(\bU_h)$ and $\CC_2(\bU_h)$ are defined as

$\begin{equation} \CC_1(\bU_h) = \begin{pmatrix} \widetilde{\CC}_1(\bmv_h) & 0 \\ 0 & 0 \end{pmatrix}, \qquad \CC_2(\bU_h) = \begin{pmatrix} \widetilde{\CC}_2(\bmv_h) & 0 \\ 0 & 0 \end{pmatrix} \end{equation}$

with

$\begin{equation} (\widetilde{\CC}_1(\bmv_h))_{ij} = \left( \rho \bmv_h \cdot \nabla \bm{\varphi}_j, \bm{\varphi}_i\right), \qquad (\widetilde{\CC}_2(\bmv_h))_{ij} = \left( \rho \bm{\varphi}_j \cdot \nabla \bmv_h, \bm{\varphi}_i\right) \end{equation}$

for $\qquad i,j = 1, \dots, \Nh[v]$.

Semi-implicit BDF time discretization

We now consider the BDF scheme with semi-implicit treatment of the convective term. This approach allows to mitigate the computational burden associated to the use of a fully implicit BDF scheme by linearizing the nonlinear convective terms. The linearization is done by extrapolating the convective velocity via an extrapolation formula of the same order of the BDF used.

We approximate the convective velocity at time $t^{n+1}$ with the Lagrange polynomial interpolating $\bmv^n_h, \dots, \bmv^{n-\sigma+1}_h$ evaluated at time $t_{n+1}$. We thus obtain the following expression for the extrapolated velocity at time $t^{n+1}$

$\begin{equation} \label{eq:NSt_extrV} \bmv_h^{n,*} = \begin{cases} \bmv_h^n, & \text{ if } n \geq 0, \quad \text{for } \sigma =1 \text{ (BDF1)}, \\ 2\bmv_h^n - \bmv_h^{n-1}, & \text{ if } n \geq 1, \quad \text{for } \sigma =2 \text{ (BDF2)}. \end{cases} \end{equation}$

The (fully discrete) semi-implicit BDF approximation of the Navier-Stokes equations reads: given $\bmv_h^{n}, \dots, \bmv_h^{n+1-\sigma}$, for $n \geq \sigma-1$ find $(\bmv_h^{n+1}, p_h^{n+1}) \in V_{h} \times M_h$ such that

$\begin{multline} \label{eq:NSt_BDF_STABLEsi} \left( \rho \frac{ \alpha \bmv_h^{n+1} - \bmv_h^{n,\sigma}}{\Delta t}, \bmw_h \right) + \left(\rho \bmv_h^{n,*} \cdot \nabla \bmv_h^{n+1} ,\bmw_h\right) + \left(\mu (\nabla \bmv_h^{n+1} + (\nabla \bmv^{n+1}_h)^T),\nabla \bmw_h\right) \\ - \left(p_h^{n+1}, \nabla \cdot \bmw_h\right) + \left(\nabla \cdot \bmv_h^{n+1},q_h\right) = g( (\bmw_h, q_h); t_{n+1} ) + (\bm{h}^{n+1}, \bmw_h)_{\Gamma_N} \end{multline}$

for all $(\bmw_h,q_h) \in V_h \times M_h$.

Algebraic formulation

At the algebraic level, we obtain the following linear system to be solved at each timestep:

$\begin{equation} \label{eq:NSt_BDF_STABLEsi_ALGtierI} \left(\frac{ \rho \alpha }{\Delta t}\MM + \AA + \CC_1( \bU^{n,*}_h) \right) \bU^{n+1}_h = \frac{\rho}{\Delta t} \MM\bU_h^{n,\sigma} + \bm{f}(t_{n+1}) . \end{equation}$

SUPG-stabilized finite element approximation

Low-order approximation spaces (such as $\mathbb{P}_1$-$\mathbb{P}_1$ spaces) represent an attractive option for the apprimation of Navier-Stokes equations, since they mitigate the required computational effort; however, they do not satisfy the inf-sup condition \eqref{eq:NSt_infsup}. As a remedy, one can resorts to suitable pressure stabilizations, which allow to circumvent \eqref{eq:NSt_infsup}. Nevertheless, pressure stabilizations turn out to be inappropriate when dealing with advection dominated flows, since in this case additional terms are required to enhance the stability with respect to the convective terms and to control the incompressibility constraint.

For these reasons, we resort to the streamline upwind Petrov-Galerkin (SUPG) stabilization – formulated as in the Variational Multiscale framework – which satisfies all these requisites. To this end, let us first introduce the finite element space

$X_h^r = \{ w_h \in C^0(\overline{\Omega}) : w_h|_{K} \in \mathbb{P}^r \ \forall K \in \mathcal{T}_h \}.$

Then, we define $V_h = V \cap [X_h^r]^d$, $V_{D,h} = V_D \cap [X_h^r]^d$ and $M_h = M \cap X_h^r$; note that we use equal-order FE spaces for the velocity and pressure variables. We also introduce the strong residuals $\bm{r}_M(\bmv_h,p_h)$ and $r_C(\bmv_h)$ of the momentum and continuity equations, respectively: \begin{equation} \label{eq:NSt_ResM} \bm{r}_M(\bmv_h,p_h) = \rho \frac{\partial \bmv_h}{ \partial t} + \rho \bmv_h \cdot \nabla \bmv_h + \nabla p - \mu \Delta \bmv_h, \end{equation} \begin{equation} \label{eq:NSt_ResC} r_C(\bmv_h) = \nabla \cdot \bmv_h. \end{equation}

The semi-discrete SUPG formulation of the Navier-Stokes equations reads: for all $t \in (0,T]$, find $(\bmv_h, p_h) \in V_{D,h} \times M_h$ such that

$\begin{multline} \label{eq:NSt_weak_SUPGh} \left( \rho \frac{\partial \bmv_h}{\partial t}, \bmw_h \right) + ( \rho \bmv_h \cdot \nabla \bmv_h ,\bmw_h) + (\mu (\nabla \bmv_h + \nabla \bmv^T_h),\nabla \bmw_h) - (p_h, \nabla \cdot \bmw_h) + (\nabla \cdot \bmv_h,q_h) \\[3pt] + \sum_{K \in \mathcal{T}_h} (\tau_M \bm{r}_M(\bmv_h,p_h), \rho \bmv_h \cdot \nabla \bmw_h + \nabla q_h)_{K} + \sum_{K \in \mathcal{T}_h} (\tau_C r_C(\bmv_h), \nabla \cdot \bmw_h)_{K} = 0 \end{multline}$

for all $(\bmw_h,q_h) \in V_h \times M_h$, with $\bmv_h(0) = \bmv_0$.

The stabilization parameters $\tau_M = \tau_M(\bmv_h)$ and $\tau_C = \tau_C(\bmv_h)$ are defined (element-wise) as:

$\begin{equation} \label{eq:NSt_tauMC} \tau_M = \bigg( \frac{\rho^2\sigma^2}{(\Delta t)^2} + \rho^2 \bmv_h \cdot \bvec{G}_K \bmv_h + C_I \mu^2 \bvec{G}_K : \bvec{G}_K \bigg)^{-1/2}, \qquad \tau_C = \big(\tau_M \bvec{g}_K \cdot \bvec{g}_K\big)^{-1}, \end{equation}$

where $C_I = 60 \cdot 2^{r-2}$, $\sigma$ is a constant equal to the order of the time discretization and $\Delta t$ is the time step that will be chosen for the time discretization. Moreover, $\bvec{G}_K$ and $\bvec{g}_K$ are metric tensors of the computational domain, which can be derived from the inverse Jacobian of the mapping between the reference and physical elements as

$(\bvec{g}_K)_{i} = \sum_{j=1}^d \frac{\partial \xi_j}{\partial x_i}, \qquad (\bvec{G}_K)_{ij} = \sum_{l=1}^d \frac{\partial \xi_l}{\partial x_i} \frac{\partial \xi_l}{\partial x_j}, \qquad i,j = 1,\dots, d,$

being $\xi_i$ and $x_i$ the reference and physical coordinates, respectively.

Semi-implicit BDF time discretization

For the time discretization of \eqref{eq:NSt_weak_SUPGh} we first consider the BDF scheme with semi-implicit treatment of the convective terms proposed in [Forti, Dede; 2015].

The (fully discrete) semi-implicit BDF-SUPG approximation of the Navier-Stokes equations reads: given $\bmv_h^{n}, \dots, \bmv_h^{n+1-\sigma}$, for $n \geq \sigma-1$ find $(\bmv_h^{n+1}, p_h^{n+1}) \in V_{h} \times M_h$ such that

$\begin{multline} \label{eq:NSt_BDF_SUPG} \left( \rho \frac{ \alpha \bmv_h^{n+1} - \bmv_h^{n,\sigma}}{\Delta t}, \bmw_h \right) + \left(\rho \bmv_h^{n,*} \cdot \nabla \bmv_h^{n+1} ,\bmw_h\right) + \left(\mu (\nabla \bmv_h^{n+1} + (\nabla \bmv^{n+1}_h)^T),\nabla \bmw_h\right) \\ - \left(p_h^{n+1}, \nabla \cdot \bmw_h\right) + \left(\nabla \cdot \bmv_h^{n+1},q_h\right) + \sum_{K} (\tau_M ^* \bm{r}_M^*(\bmv_h^{n+1},p_h^{n+1}), \rho \bmv_h^{n,*} \cdot \nabla \bmw_h + \nabla q_h)_{K} \\ + \sum_{K} (\tau_C^* r_C(\bmv_h^{n+1}), \nabla \cdot \bmw_h)_{K} = g( (\bmw_h, q_h); t_{n+1} ) \qquad \forall (\bmw_h,q_h) \in V_h \times M_h, \end{multline}$

where

$\begin{equation} \label{eq:NSt_ResMd} \bm{r}_M^*(\bmv_h^{n+1},p_h^{n+1}) = \rho \frac{\alpha \bmv_h^{n+1} - \bmv_h^{n,\sigma}}{ \Delta t} + \rho \bmv_h^{n,*} \cdot \nabla \bmv_h^{n+1} + \nabla p^{n+1} - \mu \Delta \bmv_h^{n+1} \end{equation}$

is the residual of the momentum equation, and

$\begin{equation} \label{eq:NSt_tauMCd} \tau_M^* = \bigg( \frac{\rho^2 \sigma^2}{(\Delta t)^2} + \rho^2 \bmv_h^{n,*} \cdot \bvec{G}_K \bmv_h^{n,*} + C_I \mu^2 \bvec{G}_K : \bvec{G}_K \bigg)^{-1/2}, \qquad \tau_C^* = \big(\tau_M^* \bvec{g}_K \cdot \bvec{g}_K \big)^{-1} \end{equation}$

are the stabilization parameters.

Algebraic formulation

The algebraic formulation of \eqref{eq:NSt_BDF_SUPG} reads: given $\bU_h^{n}, \dots, \bU_h^{n+1-\sigma}$ , for $n \geq \sigma-1$ find $\bU^{n+1}_h \in \mathbb{R}^{\Nh}$ such that

$\begin{multline} \label{eq:NSt_BDF_SUPG_ALGtierI} \left(\frac{ \rho \alpha }{\Delta t}\MM + \AA + \CC(\bU_h^{*}) + \SS(\bU_h^{n,*}) \right) \bU^{n+1}_h \\ = \frac{\rho}{\Delta t} \MM\bU_h^{n,\sigma} + \mathbf{f}(t_{n+1}) + \bvec{f}^S(\bU^{n,*}_h, \bU_h^{n,\sigma}) . \end{multline}$

The matrix $\SS(\bU^{n,*}_h) \in \mathbb{R}^{\Nh \times \Nh}$ encoding the SUPG operator is defined as

$\SS(\bU_h^{n,*}) = \begin{pmatrix} \SS_{vv}(\bmv_h^{n,*}) & \SS_{vp}(\bmv_h^{n,*}) \\[2pt] \SS_{pv}(\bmv_h^{n,*}) & \SS_{pp}(\bmv_h^{n,*}) \end{pmatrix},$

with

$\begin{gather*} \left(\SS_{vv}(\bmv_h^{n,*})\right)_{ij} = \left( \tau_M^* \left(\tfrac{\rho \alpha}{ \Delta t}\bm{\varphi}_j + \rho \bmv_h^{n,*} \cdot \nabla \bm{\varphi}_j - \mu \Delta \bm{\varphi}_j\right), \rho \bmv_h^{n,*} \cdot \nabla \bm{\varphi}_i\right) + \left( \tau_C^* \nabla \cdot \bm{\varphi}_j, \nabla \cdot \bm{\varphi}_i\right), \\[5pt] \left(\SS_{pv}(\bmv_h^{n,*})\right)_{kj} = \left( \tau_M^*\left(\tfrac{\rho \alpha}{ \Delta t}\bm{\varphi}_j + \rho \bmv_h^{n,*} \cdot \nabla \bm{\varphi}_j - \mu \Delta \bm{\varphi}_j\right), \nabla \eta_k \right),\\[5pt] \left(\SS_{vp}(\bmv_h^{n,*})\right)_{ik} = \left( \tau_M^* \nabla \eta_k , \rho \bmv_h^{n,*} \cdot \nabla \bm{\varphi}_i\right), \qquad \left(\SS_{pp}(\bmv_h^{n,*}) \right)_{kl} = \left( \tau_M^* \nabla \eta_l, \nabla \eta_k \right), \end{gather*}$

for $i,j=1, \dots, \Nh[v]$ and $k,l=1,\dots, \Nh[p]$. The SUPG contribute to the right-hand side is instead given by $\mathbf{f}^S(\bU_h^{n,*}, \bU_h^{n,\sigma}) = (\mathbf{f}^S_v , \mathbf{f}^S_p)^T$ , where

$(\mathbf{f}^S_v)_i = \left( \tau_M^* \tfrac{\rho}{\Delta t} \bmv_h^{n,\sigma}, \rho \bmv_h^{n,*} \cdot \nabla \bm{\varphi}_i \right), \qquad (\mathbf{f}^S_p)_k = \left( \tau_M^* \tfrac{\rho}{\Delta t} \bmv_h^{n,\sigma}, \nabla \eta_k \right).$

Thanks to the semi-implicit treatment of the nonlinear terms, the fully-discrete system \eqref{eq:NSt_BDF_SUPG} yields a linear problem – rather than a nonlinear one – to be solved at each time $t_n$. Moreover, the matrices $\MM$ and $\AA$ are constant in time and independent of $\bU_h^{n,*}$ , so that they can be assembled once and for all at $t=t_0$. On the other hand, the matrices $\CC$ and $\SS$ depend on $\bU_h^{n,*}$ and need to be assembled at each time step.

Fully implicit BDF time discretization

For the time discretization of \eqref{eq:NSt_weak_SUPGh}, we now consider the fully implicit BDF scheme.

The (fully discrete) implicit BDF-SUPG approximation of the Navier-Stokes equations reads: given $\bmv_h^{n}, \dots, \bmv_h^{n+1-\sigma}$, for $n \geq \sigma-1$ find $(\bmv_h^{n+1}, p_h^{n+1}) \in V_{h} \times M_h$ such that

$\begin{multline} \label{eq:NSt_BDF_SUPGfi} \left( \rho \frac{ \alpha \bmv_h^{n+1} - \bmv_h^{n,\sigma}}{\Delta t}, \bmw_h \right) + \left(\rho \bmv_h^{n+1} \cdot \nabla \bmv_h^{n+1} ,\bmw_h\right) + \left(\mu (\nabla \bmv_h^{n+1} + (\nabla \bmv^{n+1}_h)^T),\nabla \bmw_h\right) \\ - \left(p_h^{n+1}, \nabla \cdot \bmw_h\right) + \left(\nabla \cdot \bmv_h^{n+1},q_h\right) + \sum_{K} (\tau_M \bm{r}_M(\bmv_h^{n+1},p_h^{n+1}), \rho \bmv_h^{n+1} \cdot \nabla \bmw_h + \nabla q_h)_{K} \\ + \sum_{K} (\tau_C r_C(\bmv_h^{n+1}), \nabla \cdot \bmw_h)_{K} = g( (\bmw_h, q_h); t_{n+1} ) + (\bm{h}^{n+1}, \bmw_h) \qquad \forall (\bmw_h,q_h) \in V_h \times M_h, \end{multline}$

where

$\begin{equation} \bm{r}_M(\bmv_h^{n+1},p_h^{n+1}) = \rho \frac{\alpha \bmv_h^{n+1} - \bmv_h^{n,\sigma}}{ \Delta t} + \rho \bmv_h^{n+1} \cdot \nabla \bmv_h^{n+1} + \nabla p^{n+1} - \mu \Delta \bmv_h^{n+1} \end{equation}$

is the residual of the momentum equation.

Algebraic formulation and Newton iterations

The algebraic formulation of \eqref{eq:NSt_BDF_SUPGfi} reads: given $\bU_h^{n}, \dots, \bU_h^{n+1-\sigma}$ , for $n \geq \sigma-1$ find $\bU^{n+1}_h \in \mathbb{R}^{\Nh}$ such that

$\begin{equation} \RR(\bU^{n+1}_h) = \frac{ \rho }{\Delta t}\MM \left(\alpha\bU^{n+1}_h - \bU_h^{n,\sigma} \right) + \AA \bU^{n+1}_h + \mathbf{c}( \bU^{n+1}_h) + \mathbf{s}( \bU^{n+1}_h) - \bm{f}(t_{n+1}) = \mathbf{0}. \label{eq:NSt_BDF_SUPGfi_ALGtierI} \end{equation}$

The nonlinear term $\mathbf{s}(\bU^{n+1}_h)$ accounts for the SUPG stabilization.

At each timestep, this nonlinear system of equations is solved by means of Newton’s method: given $\bU_h^{n+1,0} = \bU_h^{n}$, for $k=0,1,\dots$ until convergence we seek $\delta \bU$ such that

$\begin{equation} \JJ(\bU_h^{n+1,k}) \, \delta \bU = - \RR(\bU_h^{n+1,k}), \end{equation}$

and then set $\bU_h^{n+1,k+1} = \bU_h^{n+1,k} + \delta \bU$. Here,

$\begin{equation} \JJ = \frac{ \rho \alpha }{\Delta t}\MM + \AA + \CC_1( \bU^{n+1,k}_h) + \CC_2( \bU^{n+1,k}_h) + \frac{\partial \mathbf{s}}{\partial \bU} \left( {\bU^{n+1,k}_h} \right). \end{equation}$

References

Navier-Stokes equations:

A. Quarteroni and A. Valli. Numerical Approximation of Partial Differential Equations. Springer-Verlag, Berlin-Heidelberg, 1994.
H.C. Elman, D.J. Silvester, and A.J. Wathen. Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics. Oxford University Press, New York, 2004.

SUPG stabilization:

T.J.R. Hughes, G. Scovazzi, and L.P. Franca. Multiscale and stabilized methods. Wiley Online Library, 2004.
Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali, and G. Scovazzi. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Comput. Methods Appl. Mech. Engrg., 197(1-4):173 - 201, 2007.

Semi-implicit treatment of the convective term:

P. Gervasio, F. Saleri, and A. Veneziani. Algebraic fractional-step schemes with spectral methods for the incompressible Navier–Stokes equations. J. Comput. Phys., 214(1):347–365, 2006.
D. Forti and L. Dede. Semi-implicit BDF time discretization of the Navier-Stokes equations with VMS-LES modeling in a high performance computing framework. Comput. Fluids, 117(0):168–182, 2015.

BDF time discretization:

K. Eleda Brenan, S.L. Campbell, and L. R. Petzold. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. SIAM, Philadelphia, 1996.
A. Quarteroni, R. Sacco, and F. Saleri. Numerical Mathematics. Springer, New York, second edition, 2007.