The Reduced NS-α Model for Incompressible Flow: A Review of Recent Progress

Abstract

This paper gives a review of recent results for the reduced Navier–Stokes-$\alpha$ (rNS-$\alpha$) model of incompressible flow. The model was recently developed as a numerical approximation to the well known Navier–Stokes-$\alpha$ model, for the purpose of more efficiently computations in the $C^0$ finite element setting. Its performance in initial numerical tests was remarkable, which led to analytical studies and further numerical tests, all of which provided excellent results. This paper reviews the main results established thus far for rNS-$\alpha$, and presents some open problems for future work.

1. Introduction

This paper presents a review of results for the recently introduced reduced Navier–Stokes-$\alpha$ (rNS-$\alpha$) model of incompressible, viscous flow, given by

$$\begin{array}{cc}\hfill -{\alpha }^2\Delta w_t+w_t+(\nabla \times Dw)\times w+\nabla p-\nu \Delta Dw& =f,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \nabla ·w& =0,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill w\left(0\right)& =w_0,\hfill \end{array}$$

(3)

where w represents velocity, p a Bernoulli-type pressure, f a body force, and D a deconvolution operator intended to approximate the inverse of the Helmholtz (also called $\alpha$) differential filter $F={(-{\alpha }^2\Delta +I)}^{-1}$. The parameter $\alpha$ > 0 represents the filtering radius.

For simplicity, only van Cittert approximate deconvolution,

$$\begin{array}{c}\hfill D:=D_N={\displaystyle \sum _{n=0}^N}{(I-F)}^n,\end{array}$$

is considered herein. This is the most common type used in fluid flow modeling [1,2,3]. Other types of approximate deconvolution operators, such as Tikhonov–Lavrentiev, may yield similar results, but all numerical tests done to date for rNS-$\alpha$ have used van Cittert. Some advantages of van Cittert, which likely make it the most common, are that it is simple to use, easy to analyze, is formally very high accuracy, and has performed very well in many computational tests for several models [1,2,3]. Note that the groundbreaking idea of using approximate deconvolution in fluid models is due to Stolz, Adams and Kleiser [4,5,6,7]. Typically N is chosen small, and for all numerical simulations run so far with rNS-$\alpha$, $N\le$ 2.

The rNS-$\alpha$ model was proposed in [8] as a numerical approximation to the well-known NS-$\alpha$ model (see e.g., [9,10,11,12]) that is more efficiently computable in a $C^0$ finite element framework. NS-$\alpha$ is given by

$$\begin{array}{ccc}\hfill v_t+(\nabla \times v)\times \overline{v}+\nabla p-\nu \Delta v& =& f,\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \nabla ·v=\nabla ·\overline{v}& =& 0,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill -{\alpha }^2\Delta \overline{v}+\overline{v}-v& =& 0.\hfill \end{array}$$

(6)

This model (also known as Camassa–Holm equations [11,12,13]) has extensive and attractive theoretical properties (see [9,10] and references therein), including well-posedness [10], frame invariance [14], adherence to Kelvin’s circulation theorem [9], conserving a model energy and helicity [9], requiring significantly fewer degrees of freedom than direct numerical simulation (DNS) of Navier–Stokes [9], and accurately predicting scalings of the turbulent boundary layer [15]. Most other fluid models do not have all (or even most) of these properties, and thus NS-$\alpha$ is widely considered to be a very ‘physically accurate’ model. However, how to construct stable, efficient algorithms for it in the $C^0$ finite element setting is an open problem. The essential issue is that the filtering seemingly cannot be decoupled from the momentum-mass system in a stable way, leaving the practitioner to solve the nonlinear problem at each time step. Not only does this create the need for many system solves at each time step, but if one wishes to use a Newton iteration, then the linear systems must solve simultaneously for the velocity, filtered velocity, and pressure. However, rNS-$\alpha$ can be computed efficiently in this setting, in the sense that it can stably decouple the filtering from the mass/momentum system, solve only one mass/momentum system at each time step, and obtain optimal accuracy. Initial testing of rNS-$\alpha$ in [8,16,17] revealed outstanding results for turbulent channel flow simulation and flow past a cylinder.

The purpose of this paper is to review the recent analytical and numerical results for rNS-$\alpha$. In Section 2, the derivation of the model is shown and how it is related to other models and stabilizations. Section 3 reviews the analytical results established for the model, including a priori bounds, energy trapping, global well-posedness, energy dissipation rate, and energy spectra. Results for sensitivity of the model with respect to the filtering radius are discussed in Section 4. Numerical schemes and their analysis are given in Section 5, both for rNS-$\alpha$ and its associated sensitivity equations from Section 4. Finally, in Section 6, results of numerical tests for channel flow past a cylinder and turbulent channel flow are discussed. Section 7 presents some open problems for possible future work.

2. Derivation of the Model and Connections to Other Models

It will be helpful for the later discussion to present the derivation found in [8] of the rNS-$\alpha$ model, which starts from the NS-$\alpha$ model. Recall from above that the NS-$\alpha$ model is defined as follows, after writing $Fv=\overline{v}$:

$$\begin{array}{ccc}\hfill v_t+(\nabla \times v)\times Fv+\nabla q-\nu \Delta v& =& f,\hfill \\ \hfill \nabla ·Fv& =& 0.\hfill \end{array}$$

It is discussed in detail in [8] why this model is difficult to efficiently compute with in a $C^0$ finite element framework.

The rNS-$\alpha$ model is created through a series of transformations to NS-$\alpha$, as follows. By using $I=F^{-1}F=(-{\alpha }^2\Delta +I)F$ in NS-$\alpha$, one gets

$$\begin{array}{ccc}\hfill -{\alpha }^2\Delta F{v}_t+F{v}_t+(\nabla \times F^{-1}Fv)\times Fv+\nabla q-\nu \Delta F^{-1}Fv& =& f,\hfill \\ \hfill \nabla ·Fv& =& 0.\hfill \end{array}$$

Denote $\tilde{w}=Fv$, which provides

$$\begin{array}{ccc}\hfill -{\alpha }^2\Delta {\tilde{w}}_t+w_t+(\nabla \times F^{-1}\tilde{w})\times \tilde{w}+\nabla q-\nu \Delta F^{-1}\tilde{w}& =& f,\hfill \\ \hfill \nabla ·\tilde{w}& =& 0.\hfill \end{array}$$

Next, use the approximation of the deconvolution operator, i.e., that $F^{-1}\tilde{w}\approx D\tilde{w}$, to get the rNS-$\alpha$ model, written in terms of w and p after renaming variables due to the approximation:

$$\begin{array}{ccc}\hfill -{\alpha }^2\Delta w_t+w_t+(\nabla \times Dw)\times w+\nabla p-\nu \Delta Dw& =& f,\hfill \\ \hfill \nabla ·w& =& 0.\hfill \end{array}$$

2.1. Connections to Other Models

There are several models closely related to rNS-$\alpha$ that have been studied recently, including NS-$\alpha$ [9,11,12,13,18,19], NS–Voigt [20,21,22], and various other regularization models. It is important to make connections between the different models when possible, so existing results of various models can be fully utilized. It is clear from Equations (1)–(3) with $D=D_N$, that for any N, choosing $\alpha =0$ recovers the Navier–Stokes equations. Thus, if a mesh sufficient for a DNS of Navier–Stokes is used with rNS-$\alpha$ and $\alpha$ is chosen on the order of the mesh width, then the solution is expected to match the Navier–Stokes solution closely.

Deconvolution theory from [2] shows that $N=\infty$ formally recovers NS-$\alpha$, and the choice $N=0$ recovers the NS–Voigt model [20,21,22]. Hence, one can interpret the rNS-$\alpha$ model as being ‘in between’ NS-$\alpha$ and NS–Voigt, in the sense of deconvolution order. These connections are summarized in the following

Table 1. Connections summary.

$\mathit{\alpha }$	N	rNS-$\mathit{\alpha }$ Model Reduces to
0	$[0,\infty ]$	Navier–Stokes
>0	0	NS–Voigt
>0	∞	NS-$\alpha$

.

A second important connection for rNS-$\alpha$ is to the generalized $\alpha$-models studied in [23]. By replacing w with $Fv$ and using the notation of [23], rNS-$\alpha$ fits the general form with $A=-\nu \Delta DF$, $M=F$, $N=DF$, and $\chi =1$. Work done for these models reveals their strengths and weaknesses, and thus identifying connections between models is critical, as they can be extended to rNS-$\alpha$. For example, the main arguments in [8] for well-posedness of rNS-$\alpha$ are similar to those made for NS–Voigt in [22,24].

When compared to other $\alpha$-models, rNS-$\alpha$ has displayed better results on wall bounded turbulence, as demonstrated in [8,16,25]. A key difference between these models lies in the viscous term, which has the form $-\nu \Delta D_{N}w$ instead of the standard $-\nu \Delta w$. At the continuous level, the rNS-$\alpha$ viscosity provides no extra smoothing, but it does provide additional dispersion when compared to the standard viscous term. This is due to the fact that $\parallel D_N\parallel \ge 1$, and this can be seen at high wave numbers since here $\parallel D_N\parallel \approx (N+1)$ [2].

Further insight is gained by expanding the deconvolution operator of the viscous term. Following [26], one can write

$$\begin{array}{ccc}\hfill -\nu \Delta D_{N}w& =& -\nu \Delta w-\nu \Delta D_{N-1}(w-\overline{w})\hfill \end{array}$$

(7)

$$\begin{array}{ccc}& =& -\nu \Delta w-\nu \Delta D_{N-1}{w}^{\prime }\hfill \end{array}$$

(8)

$$\begin{array}{ccc}& =& -\nu \Delta w+{\alpha }^2\nu {\Delta }^2{D}_{N-1}\overline{w},\hfill \end{array}$$

(9)

where $w^{\prime }:=w-\overline{w}$ can be considered as fluctuations about the (filtered) mean velocity [1]. In Equation (8), the ‘extra’ term $-\nu \Delta D_{N-1}{w}^{\prime }$ reveals a connection to variational multiscale models (VMMs). It takes the form of a viscosity for velocity fluctuations, which is in the same spirit as VMMs, but differs in that here fluctuations are defined with filtered quantities as means. Interestingly, similar dissipation/penalty/stabilizations for the advection equation [27] and in turbulence simulations [28] were recently considered. These were also constructed by using filtering to define means and fluctuations, and have proven to be very successful at increasing accuracy of simulations on coarse discretizations.

One can also observe the effect of the extra term on the energy balance. Assuming periodic or no-slip boundary conditions for both the velocity and filtered velocities and using Equation (9), test rNS-$\alpha$ with w to obtain

$$\frac{1}{2}\frac{d}{dt}\left({\alpha }^2{\parallel \nabla w\parallel }^2+{\parallel w\parallel }^2\right)+\nu {\parallel \nabla w\parallel }^2+\nu {\alpha }^2(\Delta D_{N-1}^{1/2}\overline{w},\Delta D_{N-1}^{1/2}w)=(f,w),$$

(10)

and using the filter definition yields the energy equality (from regularity results in [8], all operations are justified)

$$\frac{1}{2}\frac{d}{dt}\left({\alpha }^2{\parallel \nabla w\parallel }^2+{\parallel w\parallel }^2\right)+{\nu \parallel \nabla w\parallel }^2+\nu {\alpha }^2{\parallel \Delta D_{N-1}^{1/2}\overline{w}\parallel }^2+\nu {\alpha }^4{\parallel \nabla \Delta D_{N-1}^{1/2}\overline{w}\parallel }^2=(f,w).$$

(11)

This demonstrates the effect on the energy dissipation, which serves to dissipate higher order norms of $D_{N-1}^{1/2}\overline{w}$. Applying the filter definition once again, one can write

$$\frac{1}{2}\frac{d}{dt}\left({\alpha }^2{\parallel \nabla w\parallel }^2+{\parallel w\parallel }^2\right)+\nu {\parallel \nabla w\parallel }^2+\nu {\alpha }^2{\parallel D_{N-1}^{1/2}(w-\overline{w})\parallel }^2+\nu {\alpha }^4{\parallel \nabla D_{N-1}^{1/2}(w-\overline{w})\parallel }^2=(f,w),$$

(12)

which more clearly shows that this extra viscous term acts to dissipate velocity fluctuations.

3. Analytical Results

The first analysis of rNS-$\alpha$ was performed in [8], where it was shown that the model is well-posed for a fixed end time, and that regularity of solutions depended on the regularity of the data. Global in time energy and regularity results were proven in [17], and are stated below, after some initial preliminaries. The treatment of energy by the model was studied in [16]. Results for the energy conservation of the model, energy spectra, and energy dissipation were all very good: an appropriate energy quantity is conserved, the model exhibits a $k^{-5/3}$ energy cascade on the large scales in the inertial range and at $k^{-3}$ on smaller scales. Furthermore, energy is dissipated at the rate $\frac{U^3}{L}$, independent of the Reynolds number, which is consistent with true fluid flow.

Consider in this section a domain $\mathsf{\Omega }\subset {\mathbb{R}}^d$, d = 2 or 3, which is a box. The notation $\parallel ·\parallel$ and $(·,·)$ denotes the $L^2\left(\mathsf{\Omega }\right)$ norm and inner product, with all other norms being labeled clearly. Assume periodic boundary conditions on the box, which is the common setting for such energy studies, as it is typically the only case where well-posedness of solutions is known. For the energy balance and dissipation results, the results extend to the wall-bounded case if the well-posedness results also extend to the wall-bounded case, which is an open problem. The energy spectra results, however, rely heavily on the periodic setting, as they explicitly use Fourier decompositions. Most of the filtering and deconvolution results can be extended to the case of homogeneous Dirichlet boundary conditions, if filtered velocities satisfy no-slip boundary conditions. This is not widely accepted as a correct boundary condition for the filter (on the other hand, this is the boundary condition use in all successful turbulent flow simulations with rNS-$\alpha$ [8,16], and seems to work quite well).

Denote $(X,Q):=(H_{\#}^1\left(\mathsf{\Omega }\right),L_{\#}^2\left(\mathsf{\Omega }\right))\subset (H^1\left(\mathsf{\Omega }\right),L^2\left(\mathsf{\Omega }\right))$ the periodic, zero-mean velocity and pressure spaces. The space $H^{-1}\left(\mathsf{\Omega }\right)$ is the dual space of X.

3.1. Energy Bounds and Well-Posedness

Presented now are several long time energy bounds. The following was proven in [17], and shows the energy is trapped for all times if $f\in L^{\infty }((0,\infty );H^{-1}\left(\mathsf{\Omega }\right))$ and $w_0\in X$.

Lemma 1 (Energy trapping)

Suppose $f\in L^{\infty }((0,\infty );H^{-1}\left(\mathsf{\Omega }\right))$ and $w_0\in X$. Then, $\forall t>0$,

$${\alpha }^2{\parallel \nabla w\left(t\right)\parallel }^2+{\parallel w\left(t\right)\parallel }^2\le max\left\{\frac{1}{{\nu }^2\beta }{\parallel f\parallel }_{L^{\infty }((0,\infty ),H^{-1})}^2,{\alpha }^2{\parallel \nabla w_0\parallel }^2+{\parallel w_0\parallel }^2\right\}:=C_E,$$

where $\beta =min\left\{\frac{1}{2{\alpha }^2},\frac{C_{\mathsf{\Omega }}^2}{2}\right\}$. Thus, the kinetic energy is contained in the ball $B\left(0,C_E\right)$ for all time.

Under the additional restriction that $f\in L^2((0,\infty );H^{-1}\left(\mathsf{\Omega }\right))$, then energy must decay to 0. This is proven in [17], and stated in the following lemma.

Lemma 2 (Energy decay)

Suppose $f\in L^2((0,\infty );H^{-1}\left(\mathsf{\Omega }\right))$ and $w_0\in X$. Then, the energy decays to zero as $t\to \infty$:

$$\underset{t\to \infty }{lim}\left({\alpha }^2{\parallel \nabla w\left(t\right)\parallel }^2+{\parallel w\left(t\right)\parallel }^2\right)=0.$$

Global well-posedness of the model was able to be proven, using the long time a priori energy bounds above [17].

Lemma 3 (Existence and uniqueness of weak solutions)

Suppose $f\in L^2((0,\infty );H^{-1}\left(\mathsf{\Omega }\right))$ or $f\in L^{\infty }((0,\infty );H^{-1}\left(\mathsf{\Omega }\right))$ and $w_0\in X$. Then, there exists a global unique weak solution to Equations (1)–(3).

Additional regularity of the data leads to global in time higher order regularity of solutions, and is proven in [17].

Lemma 4 (Higher order estimates)

Assume a solenoidal initial condition $w_0\in H^2\left(\mathsf{\Omega }\right)\cap X$. If $f\in L^{\infty }((0,\infty ),L^2\left(\mathsf{\Omega }\right))$, then $w\in L^{\infty }((0,\infty ),H^2\left(\mathsf{\Omega }\right))$. If $f\in L^2((0,\infty ),L^2\left(\mathsf{\Omega }\right))$, then ${\parallel w\parallel }_{H^2}\to 0$ as $t\to \infty$.

Remark 1.

Assuming higher order regularity of the data, higher order regularity of solutions can be established using bootstrapping and the techniques used in [17].

3.2. Energy and Helicity Balances

Integral invariants such as energy and helicity are fundamental for a priori estimates used in existence theorems. They can also provide physical insight into a model’s behavior, as well as its physical relevance. It has long been established that these quantities are invariant in true fluid flow (i.e., in the Navier–Stokes equations), and as they are believed to be fundamental to the organization of flow structures, a good model should capture these quantities accurately.

The conserved model energy and helicity take the form:

$$\begin{array}{ccc}\hfill E_{rNS\alpha }\left(t\right)& :=& \frac{1}{2\left|\mathsf{\Omega }\right|}\left({\alpha }^2{\parallel \nabla w\left(t\right)\parallel }^2+{\parallel w\left(t\right)\parallel }^2\right),\hfill \\ \hfill H_{rNS\alpha }\left(t\right)& :=& \frac{1}{\left|\mathsf{\Omega }\right|}\left({\alpha }^2(\nabla \times D_N^{1/2}w\left(t\right),\nabla \times (\nabla \times D_N^{1/2}w\left(t\right)))+(D_N^{1/2}w\left(t\right),\nabla \times D_N^{1/2}w\left(t\right))\right).\hfill \end{array}$$

The energy of this model is the same as for NS–Voigt [22]. As w represents an approximation of $\overline{v}$, with v being the NS-$\alpha$ velocity, one can observe that the rNS-$\alpha$ energy is analogous to the NS-$\alpha$ energy $E_{NS\alpha }:=\frac{1}{2\left|\mathsf{\Omega }\right|}\left({\alpha }^2\parallel \nabla \overline{v}{\left(t\right)\parallel }^2+{\parallel \overline{v}\left(t\right)\parallel }^2\right)$ ([9,10]). By similar reasoning, the rNS-$\alpha$ helicity is connected to the NS–Voigt helicity and NS-$\alpha$ helicity.

It is stated below, and proven in [16], that the energy and helicity of the model are preserved. Furthermore, if the forcing is only spatially dependent, then the energy balance leads to a global in time bound on energy.

Theorem 1.

Let w be the solution to Equations (1)–(3) under periodic boundary conditions, with sufficiently smooth data. Then, the following model energy and helicity balances hold:

$$E_{rNS\alpha }\left(T\right)+\frac{1}{\left|\mathsf{\Omega }\right|}\nu {\int }_0^T{\parallel \nabla D_N^{1/2}w\left(t\right)\parallel }^{2}dt=E_{rNS\alpha }\left(0\right)+\frac{1}{\left|\mathsf{\Omega }\right|}{\int }_0^T(f\left(t\right),w\left(t\right))dt,$$

$$H\left(T\right)+\frac{2}{\left|\mathsf{\Omega }\right|}\nu {\int }_0^T(\nabla \times D_{N}w\left(t\right),\nabla \times (\nabla \times D_{N}w\left(t\right)))dt=H\left(0\right)+\frac{2}{\left|\mathsf{\Omega }\right|}{\int }_0^T(f\left(t\right),\nabla \times D_{N}w\left(t\right))dt.$$

In the case with vanishing viscosity and no external forcing, the model energy and helicity are exactly conserved. That is, for $v=0$ and $f=0$,

$$E_{rNS\alpha }\left(t\right)=E_{rNS\alpha }\left(0\right),H_{rNS\alpha }\left(t\right)=H_{rNS\alpha }\left(0\right).$$

If the forcing is only spatially dependent, the energy balance from Theorem 1 can be further analyzed, revealing that ${ϵ}_{rNS\alpha }:=\frac{\nu }{\left|\mathsf{\Omega }\right|}{\parallel \nabla D_N^{1/2}w\parallel }^2$ is the energy dissipation of the model. The following result is proven in [16].

Suppose that the forcing is solenoidal and only dependent on space, $f(x,t)=f\left(x\right)\in L^2\left(\mathsf{\Omega }\right)$, and the initial condition $w_0\in H^1\left(\mathsf{\Omega }\right)$. Then,

$$\begin{array}{ccc}\hfill \underset{t\in (0,\infty )}{sup}{E}_{rNS\alpha }\left(t\right)& \le & C\left(data\right)<\infty ,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\frac{1}{T}{\int }_0^T{ϵ}_{rNS\alpha }\left(t\right)dt=\frac{1}{T}{\int }_0^T\frac{\nu }{\left|\mathsf{\Omega }\right|}{\parallel \nabla D_N^{1/2}w\parallel }^{2}dt& \le & \frac{1}{T}{E}_{rNS\alpha }\left(0\right)+\frac{1}{{\left|\mathsf{\Omega }\right|}^{1/2}}\parallel f\parallel {\left(\frac{1}{T}{\int }_0^T\frac{1}{\left|\mathsf{\Omega }\right|}{\parallel w\parallel }^2\right)}^{1/2}.\end{array}$$

(14)

3.3. Energy Dissipation Rate

Consider now the rate of energy dissipation by the model, and assume a constant solenoidal forcing, i.e., $f=f\left(x\right)\in H^1\left(\mathsf{\Omega }\right)$, and $\nabla ·f=0$. The scale of the body force, large scale velocity, and length are defined as:

$$\begin{array}{ccc}\hfill F& :=& {\left(\frac{1}{\left|\mathsf{\Omega }\right|}{\parallel f\left(x\right)\parallel }^2\right)}^{1/2},\hfill \\ \hfill U& :=& {\langle \frac{1}{\left|\mathsf{\Omega }\right|}{\parallel w(x,t)\parallel }^2\rangle }^{1/2},\hfill \\ \hfill L& :=& min\left(\frac{F}{{\parallel \nabla f\parallel }_{\infty }},\frac{F}{{\left(\frac{1}{\left|\mathsf{\Omega }\right|}{|\parallel \nabla f\parallel }^2\right)}^{1/2}}\right),\hfill \end{array}$$

where $\langle \varphi \rangle :=lim{sup}_{T\to \infty }\frac{1}{T}{\int }_0^T\varphi \left(t\right)dt$ denotes the time average. It can be shown that each component of the definition for L has units of length.

The following theorem for the scaling of the time averaged energy dissipation is proven in [16]. Recall that N is assumed small, i.e., $N\le 5$.

Theorem 2.

The time averaged energy dissipation of the rNS-α model is bounded as

$$\langle {ϵ}_{rNS\alpha }\rangle =\langle \frac{\nu }{\left|\mathsf{\Omega }\right|}\parallel \nabla D_N^{1/2}w{\parallel }^2\rangle \le \frac{C_N}{2}\left(3+R{e}^{-1}\right)\frac{U^3}{L},$$

where $C_N$ is a constant dependent only on N (e.g., $C_N\le 6$ when $N=2$).

Remark 2.

The above result is consistent with K41 phenomenology, which predicts time averaged energy dissipation to scale with $\frac{U^3}{L}$ [29], with the scalings derived from the Navier–Stokes equations directly [30,31,32], and for the usual NS-α model [33].

3.4. Energy Spectra

It is a generally accepted theory in turbulence that energy is input at large scales, cascaded (but preserved) by the nonlinearity through the ‘inertial range’ of intermediate scales, and then dissipated exponentially fast by viscosity in the small scale range. For accurate computations, it is critical to resolve the flow to where viscosity takes over. For the Navier–Stokes equations, resolving all of these scales is a computationally intractable problem, and it is necessary to model. Studying the energy spectra of a model allows us to gain insight into its computability, as a successful model must have an inertial range that is shorter than that of the Navier–Stokes. The analysis and notation below follows studies done in [9,34,35,36].

Decomposing velocity into its Fourier modes yields the balance of energy

$$\frac{1}{2}\frac{d}{dt}\left((w_k,w_k)+{\alpha }^2(-\Delta w_k,w_k)\right)+\nu (-\Delta {\widehat{D}}_N\left(k\right)w_k,w_k)=T_k-T_{2k},$$

(15)

where ${\widehat{D}}_N\left(k\right)$ is the Fourier transform of $D_N$, and

$$T_k=-((\nabla \times {\widehat{D}}_N\left(k\right)w_k^{<})\times w_k,w_k)+((\nabla \times {\widehat{D}}_N\left(k\right)(w_k+w_k^{>}))\times (w_k+w_k^{>}),w_k^{<}),$$

with

$$w_k^{<}=\sum _{\left|j\right|<k}{w}_j,w_k^{>}=\sum _{\left|j\right|\ge 2k}{w}_j.$$

The quantity $T_k-T_{2k}$ represents net energy transferred into the wave numbers between $[k,2k)$. Time averaging the balance Equation (15) gives the energy transfer equation

$$\langle \nu (-\Delta {\widehat{D}}_N\left(k\right)w_k,w_k)\rangle =\langle T_k\rangle -\langle T_{2k}\rangle .$$

(16)

Define the model energy of a size $k^{-1}$ eddy, based on the rNS-$\alpha$ model energy $E_{rNS\alpha }$, as

$$E_{rNS\alpha }\left(k\right)=(1+{\alpha }^2{k}^2)\sum _{\left|j\right|=k}{\left|{\widehat{w}}_j\right|}^2,$$

and combining with Equation (16) gives

$$\langle T_k\rangle -\langle T_{2k}\rangle \sim \frac{\nu k^3{\widehat{D}}_N\left(k\right)E_{rNS\alpha }\left(k\right)}{1+{\alpha }^2{k}^2}\sim \frac{\nu k^3(N+1)E_{rNS\alpha }\left(k\right)}{1+{\alpha }^2{k}^2},$$

with the last relation coming from [2]: ${\widehat{D}}_N\left(k\right)\sim (N+1)$ for sufficiently large k. This is a key estimate for determining the inertial range, since here one expects no leakage of energy through dissipation, and hence $\langle T_k\rangle \approx \langle T_{2k}\rangle$. This suggests that increasing N leads to a shorter inertial range, which is consistent with the energy spectrum computations in Section 4.

The inertial range’s kinetic energy distribution is now considered. Begin by defining the average velocity of a size $k^{-1}$ eddy, and then relate it to model energy:

$$U_k={\langle (w_k,w_k)\rangle }^{1/2}\sim {\left({\int }_k^{2k}\frac{E_{rNS\alpha }\left(k\right)}{1+{\alpha }^2{k}^2}\right)}^{1/2}\sim {\left(\frac{k{E}_{rNS\alpha }\left(k\right)}{(1+{\alpha }^2{k}^2)}\right)}^{1/2}.$$

From the Kraichnan energy cascade theory [37], one obtains that the eddy turnover time is

$${\tau }_k:=\frac{1}{k{U}_k}=\frac{{(1+k^2{\alpha }^2)}^{1/2}}{k^{3/2}{E}_{rNS\alpha }{\left(k\right)}^{1/2}},$$

and thus the model energy dissipation rate scales like

$${ϵ}_{rNS\alpha }\sim \frac{1}{{\tau }_k}{\int }_k^{2k}{E}_{rNS\alpha }\left(k\right)\sim \frac{k^{5/2}{E}_{rNS\alpha }{\left(k\right)}^{3/2}}{{(1+k^2{\alpha }^2)}^{1/2}}.$$

Solving for $E_{rNS\alpha }\left(k\right)$, one obtains the relation

$$E_{rNS\alpha }\left(k\right)\sim \frac{{ϵ}_{rNS\alpha }^{2/3}{(1+{\alpha }^2{k}^2)}^{1/3}}{k^{5/3}},$$

and this yields a spectrum for kinetic energy ($E=\frac{1}{2}{∥w∥}^2$):

$$E\left(k\right)\sim \frac{E_{rNS\alpha }\left(k\right)}{(1+{\alpha }^2{k}^2)}\sim \frac{{ϵ}_{rNS\alpha }^{2/3}}{k^{5/3}{(1+{\alpha }^2{k}^2)}^{2/3}}.$$

Split the spectrum into two parts: if $k\alpha >O\left(1\right)$, then

$$E\left(k\right)\sim \frac{{ϵ}_{rNS\alpha }^{2/3}}{k^{5/3}{\left({\alpha }^2{k}^2\right)}^{2/3}}\sim \frac{{ϵ}_{rNS\alpha }^{2/3}}{k^3{\alpha }^{4/3}},$$

but if $k\alpha <O\left(1\right)$,

$$E\left(k\right)\sim \frac{{ϵ}_{rNS\alpha }^{2/3}}{k^{5/3}}.$$

These scalings show that on larger inertial range scales, a $k^{-5/3}$ rolloff of energy is expected (which agrees with true fluid flow). However, for wave numbers bigger than $O\left({\alpha }^{-1}\right)$, a rolloff of $k^{-3}$ is expected. These rolloffs are identical to usual NS-$\alpha$ on both the larger and smaller scales in the inertial range [9]. This implies that significantly less energy is held in higher wave numbers compared to a Navier–Stokes DNS, which makes the model more computable than the Navier-Stokes Equations (NSE). Numerical tests in [16] illustrate these scalings. Another important takeaway is that since the rolloffs mirror those of NS-$\alpha$, one should expect computational cost (in terms of mesh width resolution) to be very similar NS-$\alpha$, at least in this periodic setting. This means that total cost in the $C^0$ finite element setting should be much less for rNS-$\alpha$ compared to NS-$\alpha$, since the former is more efficiently computable, as discussed above.

4. Sensitivity of the Model with Respect to the Filtering Radius $\mathbf{\alpha }$

The parameter $\alpha$ in rNS-$\alpha$ represents a filtering radius, and changing $\alpha$ in computations can lead to very different solutions. Error analysis of the model’s discretization above suggests $\alpha =O\left(h\right)$ is a good choice, and choosing $\alpha$ much smaller than this would lead to no regularization, since such a discrete filter would effectively have 0 filtering radius. Throughout the literature, $\alpha$ is often chosen slightly larger than h, such as 2 h or 6 h (our experience makes us favor the choice of 2 h). Due to these various choices, it makes sense to study sensitivity of the model to changes in $\alpha$, and in this section, a sensitivity equation for $\frac{\partial }{\partial \alpha }w$ is derived and analyzed. In later sections, efficient algorithms are proposed for computing it, and finally computations are performed. This section follows work done in [17].

The sensitivity equation for $\frac{\partial }{\partial \alpha }w$ is derived by applying $\frac{\partial }{\partial \alpha }$ to rNS-$\alpha$, then setting $s=\partial w/\partial \alpha$ and $r=\partial p/\partial \alpha$. This gives

$$\begin{array}{cc}\hfill -{\alpha }^2\Delta s_t+s_t+\left(\nabla \times \frac{\partial }{\partial \alpha }{D}_{N}w\right)\times w+(\nabla \times \frac{\partial }{\partial \alpha }{D}_{N}w)\times s+\nabla r-\nu \Delta \frac{\partial }{\partial \alpha }{D}_{N}w& =2\alpha \Delta w_t,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \nabla ·s& =0,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill s\left(0\right)& =0.\hfill \end{array}$$

(19)

This system is reduced by eliminating variables for filtered velocities’ sensitivities after evaluating $\frac{\partial }{\partial \alpha }{D}_{N}w$, and writing the system in terms of $s=\partial w/\partial \alpha$. The system depends on N, and for the first few values of N, one gets

N	${\mathit{D}}_{\mathit{N}}\mathit{w}$	$\frac{\mathit{\partial }\mathit{Dw}}{\mathit{\partial }\mathit{\alpha }}$
0	I	$Ds$
1	$2w-\overline{w}$	$Ds-2\alpha \Delta \overline{\overline{w}}$
2	$3w-3\overline{w}+\overline{\overline{w}}$	$Ds-6\alpha \Delta \overline{\overline{w}}+4\alpha \Delta \overline{\overline{\overline{w}}}$
3	$4w-6\overline{w}+4\overline{\overline{w}}-3\overline{\overline{\overline{w}}}$	$Ds-12\alpha \Delta \overline{\overline{w}}+16\alpha \Delta \overline{\overline{\overline{w}}}-18\alpha \Delta \overline{\overline{\overline{\overline{w}}}}$

Note that for general N, one can write $\frac{\partial D_{N}w}{\partial \alpha }$ can be written as $D_{N}s$ plus a linear combination of filtered w’s,

$$\frac{\partial D_{N}w}{\partial \alpha }=D_{N}s+\alpha \Delta {\displaystyle \sum _{i=1}^N}{\beta }_i^{\left(N\right)}{F}^{i+1}w,$$

where the ${\beta }_i^{\left(N\right)}$’s are integers depending on N. Thus, the general system takes the form

$$\begin{array}{cc}\hfill -{\alpha }^2\Delta s_t+s_t+\left(\nabla \times D_{N}s\right)\times w+(\nabla \times D_{N}w)\times s+\nabla r& -\nu \Delta D_{N}s=\hfill \\ \hfill 2\alpha \Delta w_t+\nu \alpha {\Delta }^2\left({\displaystyle \sum _{i=1}^N}{\beta }_i^{\left(N\right)}{F}^{i+1}w\right)-& \alpha \left(\nabla \times \Delta {\displaystyle \sum _{i=1}^N}{\beta }_i^{\left(N\right)}{F}^{i+1}w\right)\times w,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \nabla ·s& =0,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill s\left(0\right)& =0.\hfill \end{array}$$

(22)

This system is proven in [17] to be well-posed in the periodic setting, and the result reads as follows:

Theorem 3 (Existence and uniqueness of weak sensitivity solutions)

Let $f\in L^{\infty }((0,\infty );L^2\left(\mathsf{\Omega }\right))$ and $w_0\in X$. Then, weak solutions to Equations (20)–(22) exist uniquely, and satisfy for $T<\infty$,

$$\begin{array}{cc}\hfill {\alpha }^2{\left|\right|\nabla s\left(T\right)\left|\right|}^2+& {\left|\right|s\left(T\right)\left|\right|}^2+\nu {\int }_0^T{\left|\right|\nabla s\left(t\right)\left|\right|}^{2}dt\le C\left(data\right).\hfill \end{array}$$

In the numerical tests section below, efficient numerical schemes for computing this sensitivity system are proposed and analyzed.

5. Numerical Schemes and Analysis

An efficient numerical scheme for approximating solutions to rNS-$\alpha$ using a $C^0$ finite element spatial discretization, together with an implicit-explicit BDF2 timestepping, is given. This time discretization linearizes the system at each time step and decouples the filtering/deconvolution equations from the momentum/mass system. After providing some preliminaries, the scheme and results for its stability, well-posedness, and convergence, are presented.

5.1. Problem Setting and Preliminaries

Assume a regular, conforming triangulation/tetrahedralization ${\tau }_h\left(\mathsf{\Omega }\right)$, and associated inf-sup stable discrete velocity-pressure spaces $(X_h,Q_h)\subset (X,Q)$ that are piecewise polynomials on each element. Define the discretely divergence free subspace by

$$V_h:=\{v_h\in X_h,(\nabla ·v_h,q_h)=0\forall q_h\in Q_h\}.$$

Discrete filtering is defined by the standard finite element discretization of the $\alpha$ filter: Given $\varphi \in L^2\left(\mathsf{\Omega }\right)$, find $F_h\varphi ={\overline{\varphi }}^h\in X_h$ satisfying

$${\alpha }^2(\nabla {\overline{\varphi }}^h,\nabla {\chi }_h)+({\overline{\varphi }}^h,{\chi }_h)=(\varphi ,{\chi }_h)\forall {\chi }_h\in X_h.$$

Discrete van Cittert deconvolution is analogous to the continuous case, but defined with the discrete filter:

$$D_N^h=\sum _{n=0}^N{(I-F_h)}^n.$$

An explicit filter operator is now given that filters known quantities exactly, but for unknown quantities, it instead filters a known second order approximation:

$$\widetilde{F_h}{w}_h^j=F_h{w}_h^{j}forj=1,2,\dots ,n,and\widetilde{F_h}{w}_h^{n+1}:=F_h(2w_h^n-w_h^{n-1}).$$

This operator defines an explicit deconvolution operator

$$\widetilde{D_N^h}{w}_h^{n+1}=\sum _{n=0}^N{(I-\widetilde{F_h})}^n.$$

Finally, the implicit-explicit (IMEX) BDF2, $C^0$ finite element scheme for the rNS-$\alpha$ model can be stated.

Algorithm 1: [BDF2] Given endtime $T>0$, timestep $\Delta t>0$, forcing $f\in L^{\infty }(0,T;H^{-1}\left(\mathsf{\Omega }\right))$, and initial velocities $w_h^{-1},w_h^0\in V_h$, set $M=\frac{T}{\Delta t}$ and for $n=1,2,\dots ,M-1$, find $(w_h^{n+1},q_h^{n+1})\in (X_h,Q_h)$ satisfying for all $(v_h,r_h)\in (X_h,Q_h)$,

$$\begin{array}{c}\hfill \frac{{\alpha }^2}{2\Delta t}\left(\nabla (3w_h^{n+1}-4w_h^n+w_h^{n-1}),\nabla v_h\right)+\frac{1}{2\Delta t}\left(3w_h^{n+1}-4w_h^n+w_h^{n-1},v_h\right)\\ \hfill +\left((\nabla \times D_N^h(2w_h^n-w_h^{n-1}))\times w_h^{n+1},v_h\right)-(q_h^{n+1},\nabla ·v_h)\end{array}$$ $$+\nu (\nabla \widetilde{D_N^h}{w}_h^{n+1},\nabla v_h)=(f\left(t^{n+1}\right),v_h),$$ (23) $$(\nabla ·w_h^{n+1},r_h)=0.$$ (24)

5.2. Analysis of the Scheme

Results of the stability and convergence of Algorithm 1 are now stated, and come from [8].

Lemma 5 (Stability)

Suppose that $\Delta t<\frac{{\alpha }^2}{4C\left(N\right)\nu }$. Then solutions of Algorithm 1 satisfy

$${\alpha }^2\parallel \nabla w_h^M{\parallel }^2+\parallel w_h^M{\parallel }^2+\sum _{n=0}^{M-1}\parallel w_h^{n+1}-2w_h^n+w_h^{n-1}{\parallel }^2+\nu \Delta t\sum _{n=1}^M{\parallel \nabla w_h^n\parallel }^2\le C\left(data\right).$$

(25)

Moreover, solutions to Algorithm 1 exist uniquely since the scheme is linear at each timestep. Thus, the algorithm is well-posed.

Remark 3.

In most settings where large eddy simulation (LES) models such as rNS-α is used, it generally holds that $\nu \ll \alpha$. Moreover, typically $N\le$ 3, giving $C\left(N\right)\le 6$. Thus, in practice, the timestep restriction will typically be mild.

Convergence of the scheme is proven in [8] and is stated in the following theorem.

Theorem 4.

Suppose $(X_h,Q_h)$ are chosen to be either $(P_k,P_{k-1})$ Taylor–Hood elements or $(P_k,P_{k-1}^{disc})$ Scott–Vogelius elements and suppose further that $(w,p)$ is a solution of Equations (1) and (2) for a given $\alpha >0$, $\nu >0$, with forcing $f\in L^{\infty }(0,T;H^{-1}\left(\mathsf{\Omega }\right))$ and initial velocity $w_0\in V_h$, satisfying regularity criteria $w\in L^{\infty }(0,T;H^{k+1}\left(\mathsf{\Omega }\right)\cap H^3\left(\mathsf{\Omega }\right))$, $p\in L^{\infty }(0,T;H^k\left(\mathsf{\Omega }\right))$, $w_t,w_{tt},w_{ttt}\in L^{\infty }(0,T;H^{k+1}\left(\mathsf{\Omega }\right)\cap H^3\left(\mathsf{\Omega }\right))$. Suppose further that the stability criterion for the timestep size is satisfied. Then the error of Algorithm 1 in approximating solutions to the rNS-α model satisfies

$$\begin{array}{c}∥w\left(T\right)-w_h^M∥+\sum _{n=0}^{M-1}∥(w\left(t^{n+1}\right)-2w\left(t^n\right)+w\left(t^{n-1}\right))-(w_h^{n+1}-2w_h^n+w_h^{n-1})∥\hfill \\ +{\alpha }^2∥\nabla (w\left(T\right)-w_h^M)∥+\nu {\left(\Delta t\sum _{n=1}^M{\parallel \nabla (w\left(t^n\right)-w_h^n)\parallel }^2\right)}^{1/2}\le C\left(\Delta t^2+h^k\right).\hfill \end{array}$$

(26)

5.3. Numerical Scheme for rNS-$\alpha$ Sensitivity

An efficient numerical scheme is now detailed for the rNS-$\alpha$ sensitivity system to be used together with the IMEX BDF2 scheme above. The key feature for efficiency is that the discrete velocity/pressure and sensitivity systems will have exactly the same coefficient matrices. Thus, preconditioners for the velocity solve can be reused, and also one can take advantage of block solvers that work with multiple right-hand sides, since the proposed algorithms allows solving for $s^n$ and $w^{n+1}$ simultaneously. The discrete scheme is now presented. Step 1 is exactly the rNS-$\alpha$ scheme from above. Step 2 solves for velocity and pressure sensitivity, and is decoupled from Step 1.

Algorithm 2: [rNS-$\alpha$ with sensitivity] Given $w_h^0,w_h^1\in V_h$, $s_h^0=s_h^1=0\in V_h$, $f\in L^{\infty }(0,T;L^2\left(\mathsf{\Omega }\right))$, and timestep $\Delta t=\frac{T}{M}$, for $n=2,3,\dots ,M-1$,

Step 1: Find $w_h^{n+1}\in V_h$ satisfying: $$\begin{array}{c}\frac{1}{2\Delta t}\left(3w_h^{n+1}-4w_h^n+w_h^{n-1},v_h\right)+\frac{{\alpha }^2}{2\Delta t}\left(\nabla (3w_h^{n+1}-4w_h^n+w_h^{n-1}),\nabla v_h\right)\hfill \\ +\left((\nabla \times D_N^h(2w_h^n-w_h^{n-1}))\times w_h^{n+1},v_h\right)+\nu \left(\nabla \widetilde{D_N^h}{w}_h^{n+1},\nabla v_h\right)=\left(f(t^{n+1}),v_h\right)\hfill \end{array}$$ (27)

Step 2: Find $s_h^{n+1}\in V_h$ satisfying: $$\begin{array}{c}\frac{1}{2\Delta t}\left(3s_h^{n+1}-4s_h^n+s_h^{n-1},v_h\right)+\frac{{\alpha }^2}{2\Delta t}\left(\nabla (3s_h^{n+1}-4s_h^n+s_h^{n-1}),\nabla v_h\right)\hfill \\ +\left((\nabla \times D_N^h(2w_h^n-w_h^{n-1}))\times s_h^{n+1},v_h\right)+\nu \left(\nabla \widetilde{D_N^h}{s}_h^{n+1},\nabla v_h\right)\\ =-\frac{\alpha }{\Delta t}\left(\nabla (3w_h^{n+1}-4w_h^n+w_h^{n-1}),\nabla v_h\right)-\left((\nabla \times D_N^h(s_h^n-s_h^{n-1}))\times w_h^{n+1},v_h\right)\\ \hfill -\alpha \left(\left(\nabla \times \Delta \sum _{i=1}^N{\beta }_i^{\left(N\right)}{F}_h^{i+1}{w}_h^{n+1}\right)\times w_h^{n+1},v_h\right)-\nu \alpha \left(\nabla \Delta \sum _{i=1}^N{\beta }_i^{\left(N\right)}{F}_h^{i+1}{w}_h^{n+1},\nabla v_h\right)\end{array}$$ (28)

Stability and convergence of the Step 1 is stated above, and it is proven in [17] that Step 2 is stable and convergent. However, there is a slightly more restrictive condition on the time step size.

Theorem 5.

Let $w_0\in H^1\left(\mathsf{\Omega }\right)$, $f\in L^{\infty }([0,T],L^2\left(\mathsf{\Omega }\right))$, and choose $\Delta t$ such that $\Delta t<min\{\frac{{\alpha }^2}{4C\left(N\right)\nu },O\left(h^{2/3}\right)\}$. Then,

$$\begin{array}{c}\hfill {∥s_h^M∥}^2+{\alpha }^2{∥\nabla s_h^M∥}^2+\nu \Delta t\sum _{n=0}^{M-1}{∥\nabla s_h^{n+1}∥}^2\le C\left(data\right).\end{array}$$

Furthermore, assuming sufficiently smooth data, the convergence estimate holds:

$${∥s\left(T\right)-s_h^M∥}^2+\nu \Delta t\sum _{n=0}^M{∥\nabla (s\left(t^n\right)-s_h^n)∥}^2\le C(\Delta t^2+h^k).$$

6. Numerical Results

Results of computations for rNS-$\alpha$ are presented in this section. First, results of convergence rate tests are given, to illustrate the convergence results above. Next, two application problems are considered, 2D flow past a cylinder and 3D turbulent channel flow. For the two application problems, comparisons of the model solutions to the literature are made, and sensitivities are also computed and discussed.

6.1. Verification of Convergence Rates

In [8], convergence rate verification tests were performed for the IMEX scheme above for rNS-$\alpha$. The test problem used was an exact rNS-$\alpha$ solution (i.e., it solves (1) and (2) with $f=0$) given by

$$w=\left(\begin{array}{c}-cos\left(n\pi x\right)sin\left(n\pi y\right)\\ sin\left(n\pi x\right)cos\left(n\pi y\right)\end{array}\right)e^{\frac{-2C_{D_N}{n}^2{\pi }^2\nu t}{1+2n^2{\pi }^2{\alpha }^2}},$$

(29)

where $C_D$ is defined for each N by $D_{N}w=C_{D_N}w$ (justification of the existence of such a constant is given in [8]).

Approximations to this exact solution were calculated in [8] on $\mathsf{\Omega }={(0,1)}^2$ and $0\le t\le T=0.1$, using Algorithm 1, with $f=0$, $N=1$, $n=1$, $\nu =1$, and $\alpha =\frac{1}{16}$ (with $C_{D_1}=2-\frac{1}{1+2{\pi }^2{\alpha }^2}$). Errors and rates for the scheme were computed, using successively refined uniform meshes and timesteps, with $(P_2,P_1)$ Taylor–Hood elements. Based on Theorem 4, convergence is expected to satisfy

$$\begin{array}{c}\hfill \parallel w-w_h{\parallel }_{2,1}:={\parallel w-w_h\parallel }_{L^2(0,T;H^1\left(\mathsf{\Omega }\right))}=O(h^2+\Delta t^2).\end{array}$$

Thus, a convergence rate of 2 is predicted if the mesh width and the time step are reduced together, by a factor of 2 at each refinement, which is observed to be the case from

Table 2. Velocity errors and rates for the convergence rate test are given in the table, and are consistent with the optimal theoretical rate of 2.

h	$\Delta \mathit{t}$	$\parallel \mathit{w}-{\mathit{w}}_{\mathit{h}}{\parallel }_{2,1}$	Rate
1/4	1/2	1.005 $\times {10}^{-1}$	-
1/8	1/4	6.553 $\times {10}^{-2}$	0.617
1/16	1/8	1.678 $\times {10}^{-2}$	1.965
1/32	1/16	4.089 $\times {10}^{-3}$	2.037
1/64	1/32	1.141 $\times {10}^{-3}$	1.841

.

6.2. 2D Channel Flow Past a Cylinder

Results of tests of the rNS-$\alpha$ model/scheme for 2D flow past a cylinder are now presented. This is a benchmark problem from [38,39,40], with the domain being a 2.2 × 0.41 rectangular channel with a radius = 0.05 cylinder, centered at $(0.2,0.2)$, as shown in Figure 1. No slip boundary conditions are prescribed on the walls and cylinder, and the inflow and outflow profiles are given by

$$\begin{array}{c}{u}_1(0,y,t)=u_1(2.2,y,t)=\frac{6}{0.{41}^2}sin(\pi t/8)y(0.41-y),\hfill \\ u_2(0,y,t)=u_2(2.2,y,t)=0.\hfill \end{array}$$

The kinematic viscosity is set as $\nu ={10}^{-3}$, the forcing $f=0$, and the flow starts from rest. The correct behavior is for a vortex street to form by t = 4 and continue to t = 8. The flow is driven by an interaction of the fluid with the cylinder, which is an important scenario for industrial flows. This flow is not turbulent, but is still very challenging for most numerical models and methods, especially on coarser meshes.

Important statistics for this flow are the predicted lift and drag forces. For a velocity field u and pressure p, lift and drag coefficients are defined by

$$\begin{array}{ccc}\hfill c_d\left(t\right)& =& \frac{2}{\rho L{U}_{max}^2}{\int }_S\left(\rho \nu \frac{\partial u_{t_S}\left(t\right)}{\partial n}{n}_y-p\left(t\right)n_x\right)dS,\hfill \\ \hfill c_l\left(t\right)& =& \frac{2}{\rho L{U}_{max}^2}{\int }_S\left(\rho \nu \frac{\partial u_{t_S}\left(t\right)}{\partial n}{n}_x-p\left(t\right)n_y\right)dS,\hfill \end{array}$$

where $u_{t_S}$ is tangential velocity, $U_{max}=1$ is the max velocity at the inlet, $L=0.1$, $\rho =1$, S is the cylinder, and $n=\langle n_x,n_y\rangle$ is the outward unit normal.

Solutions were computed using Algorithm 1 with $\Delta t=0.002$ and $({\left(P_2\right)}^2,P_1^{disc})$ elements on a very coarse barycenter mesh that provided 5104 velocity degrees of freedom (dof). This choice of elements is made because the model is implemented using a Bernoulli-type pressure, which is significantly more complex than usual pressure and is known to affect the overall error if standard element choices are used [41], but not if divergence-free elements are used. The deconvolution parameter $N=2$ was chosen, and tests were run with $\alpha$ varying from $0.010$ to $0.016$. To evaluate these solutions, values for the maximum drag $c_{d,max}$ and lift $c_{l,max}$ coefficients were calculated, and compared to those found in the resolved benchmark tests of [42]:

$$c_{d,max}^{ref}=2.950918381,c_{l,max}^{ref}=0.47787543.$$

For comparison, the max lift and drag coefficient values were calculated on the same mesh, elements, and time step by the NSE discretized by the following extrapolated Crank–Nicolson finite element method:

$$\begin{array}{ccc}\hfill \frac{1}{\Delta t}\left(u_h^{n+1}-u_h^n,v_h\right)+\left(\left(\nabla \times (\frac{3}{2}{u}_h^n-\frac{1}{2}{u}_h^{n-1})\right)\times u_h^{n+1/2},v_h\right)& & \\ \hfill -\left(p_h^{n+1/2},\nabla ·v_h\right)+\nu \left(\nabla u_h^{n+1/2},\nabla v_h\right)& =& (f(t^{n+1/2}),v_h),\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill (\nabla ·u_h^{n+1},q_h)& =& 0.\hfill \end{array}$$

(31)

One observes from

Table 3. Shown in the table are max lift and drag coefficients for various models and parameters, for 2D flow past a cylinder.

Model	$\mathit{\alpha }$	${\mathit{c}}_{\mathit{d},\mathit{max}}$	% Error	${\mathit{c}}_{\mathit{l},\mathit{max}}$	% Error
Navier–Stokes (NS)	-	3.393	15.0%	0.7499	56.9%
rNS-$\alpha$	0.011	3.074	4.2%	0.511	6.9%
rNS-$\alpha$	0.012	3.107	5.2%	0.506	5.9%
rNS-$\alpha$	0.013	3.142	6.4%	0.501	4.8%
rNS-$\alpha$	0.014	3.180	7.8%	0.497	4.0%
rNS-$\alpha$	0.015	3.220	9.2%	0.492	3.0%
rNS-$\alpha$	0.016	3.264	10.6%	0.487	1.9%

that rNS-$\alpha$ provides much better max lift and drag coefficients than the coarse mesh NSE. The reference values come from upwards of ${10}^7$ dof simulations, and it is remarkable the rNS-$\alpha$ gets this level of accuracy on such a coarse discretization.

6.2.1. Sensitivity to $\alpha$ in Flow Past a Cylinder

For this test problem, sensitivity of the solution to $\alpha$ was considered, and Algorithm 2 was computed with $\alpha$ = 0.001, $\Delta t$ = 0.001, and $(P_2,P_1^{\mathrm{disc}})$ Scott–Vogelius elements on a barycenter refined Delaunay triangulation, which provided 35,948 total dof. Plots of the sensitivity magnitude are shown in Figure 2 at various times. At early times, the sensitivity is created at the cylinder, and then as time progresses, it gets convected into the flow. By t = 6, values of $\left|s\right|$ are near 600. Creation of sensitivity at the solid boundary is not unexpected, since the sensitivity equation resembles a vorticity equation.

Lift and drag sensitivities are also calculated, defined by differentiating the lift and drag formulas above with respect to $\alpha$. For the tests above, these quantities are also calculated:

$$\begin{array}{ccc}\hfill s_d\left(t\right)& =& \frac{2}{\rho L{U}_{max}^2}{\int }_S\left(\rho \nu \frac{\partial s_{t_S}\left(t\right)}{\partial n}{n}_y-r\left(t\right)n_x\right)dS,\hfill \\ \hfill s_l\left(t\right)& =& \frac{2}{\rho L{U}_{max}^2}{\int }_S\left(\rho \nu \frac{\partial s_{t_S}\left(t\right)}{\partial n}{n}_x-r\left(t\right)n_y\right)dS.\hfill \end{array}$$

Here, $s_{t_S}$ denotes the tangential sensitivity. Shown in Figure 3 are plots of $s_d$ and $s_l$ versus time that came from this test. One observes very large sensitivity with the lift, and less so in the drag, but still significant. One also observes that the drag sensitivity is positive for most of the simulation, indicating that increasing $\alpha$ will increase drag.

6.3. Turbulent Channel Flow at $R{e}_{\tau }$ = 590

Results from rNS-$\alpha$ tests on a benchmark turbulent channel flow problem with friction Reynolds number $R{e}_{\tau }$ = 590 from [16] are now discussed. The results below show that rNS-$\alpha$ does an outstanding job at matching the DNS mean velocity profile on a very coarse mesh. To our knowledge, no other model produces results this good on such a coarse mesh. For the widely cited DNS results [43], the number of degrees of freedom used was 113 million (our tests use just 82,821). For comparison, results on the same discretization for Leray-$\alpha$, NS-$\alpha$, and NS–Voigt models were computed on the same coarse mesh in [16] and shown to be much worse than rNS-$\alpha$.

6.3.1. Problem Description and Results

The problem setup is given in detail in [8,16,44], and the interested reader is referred to those works for additional details. The domain is the box $\mathsf{\Omega }=(-2\pi ,2\pi )\times (0,2)\times \left(-2\pi /3,2\pi /3\right)$, with homogeneous Dirichlet boundary conditions enforced on the top and bottom walls $y=0$ and $y=2$, and periodic boundary conditions enforced on the remaining sides. The kinematic viscosity is set at $\nu =\frac{1}{R{e}_{\tau }}$. The initial condition is created by randomly perturbing the DNS average velocity profile pointwise. These perturbations create the turbulence as the simulation is run. The forcing is defined by $f={\langle 1,0,0\rangle }^T$, and is dynamically adjusted at each time step to maintain the bulk velocity. The simulations are run to t = 60, and the time and space average of the velocity is done from t = 20 to t = 60.

$(P_3,P_2^{disc})$ divergence-free Scott–Vogelius elements are used with a barycentric tetrahedral mesh that is weighted towards the walls, and provides 82,821 velocity dof. The time step size $\Delta t=0.002$, deconvolution N = 2, $\alpha$ was varied. The value of $\alpha =0.012$ gives an optimal (to within 0.0005) solution in the $L^2$ sense of best matching the DNS mean velocity profile; note that this value is close to the mesh resolution near the wall. Mean streamwise velocity profiles for rNS-$\alpha$ are shown in Figure 4, together with DNS data. Results this good are seemingly unique for any regularization model on such a coarse mesh, using a finite element discretization.

The $L^2(0,1)$ difference between the DNS mean streamwise velocity, and various model’s mean streamwise velocity, is shown in

Table 4. The $L^2(0,1)$ difference between the model’s predicted velocity profile and the (cubic spline of) Moser, Kim, and Moin velocity profile for $R{e}_{\tau }$ = 590, using results from the model with optimal $\alpha$ from 0.001 to 0.070, with increments of 0.001. No NS-$\alpha$ simulation with the implicit-explicit (IMEX) algorithm from the appendix remained stable for the entire simulation.

${\mathit{Re}}_{\mathit{\tau }}$ = 590
Model	Optimal $\mathbf{\alpha }$	$\parallel \langle {\mathit{u}}_{\mathit{DNS}}\rangle \left(\mathit{y}\right)-\langle {\mathit{u}}_{\mathit{model}}\rangle \left(\mathit{y}\right){\parallel }_{{\mathit{L}}^{\mathbf{2}}(\mathbf{0},\mathbf{1})}$
Navier-Stokes Equations (NSE) coarse mesh	-	1.409
rNS-$\alpha$	0.012	0.184
NS-$\alpha$	DNE	DNE
Leray-$\alpha$	0.040	1.213
NS–Voigt	0.030	1.238

(the values of $\alpha$ are optimal to within 0.0005 in the interval [0,0.1] One immediately observes the rNS-$\alpha$ is much more accurate than Leray-$\alpha$ and NS–Voigt (which both barely improve on the ‘no model’, i.e., NSE on the coarse mesh). NS-$\alpha$ blows up for each choice of $\alpha$ (note the nonlinear problem was not solved at each time step; instead, a linearization on the filtered term in the nonlinearity is done, to make its efficiency the same as the other models (see [16])).

6.4. Model Sensitivity for Turbulent Channel Flow with $R{e}_{\tau }$ = 590

rNS-$\alpha$ sensitivity was tested in [17] for this benchmark turbulent channel flow problem. Since the quantity of interest for this problem is the long time averaged velocity, a sensitivity study of rNS-$\alpha$ for this problem is now performed, to determine the long time averaged velocity sensitivity.

The same discretization as in [16] is used for this test: decoupled IMEX-BDF2-finite element discretization (step 1 of Algorithm 2), and $(P_3,P_2^{disc})$ Scott–Vogelius elements on a tetrahedral mesh that provided 82,281 velocity (and sensitivity) degrees of freedom. Both $N=2$ and $N=0$ (which is the NS–Voigt model) are tested, and $\alpha$ = 0.010, 0.015 and 0.020. Starting from the T = 60, $N=2$, $\alpha =0.015$ (turbulent) solution as the initial condition (created using the setup above), the simulation is run to T = 70 s using a time step size of $\Delta t=0.002$ (1000 total time steps). Algorithm 2 step 2 is used to calculate sensitivities at each time step, and time and space (x and z directions) averages are taken to get the mean streamwise velocity and sensitivity profiles.

From Figure 5, observe that rNS-$\alpha$ ($N=2$) gives an excellent prediction of the mean streamwise velocity profile for each choice of $\alpha$, significantly better than the NS–Voigt ($N=0$) case. For the varying $\alpha$, there is no visible change in the NS–Voigt solution, and a slight visible change in the rNS-$\alpha$ solution.

Figure 6 shows the mean streamwise sensitivity and relative sensitivity magnitudes, for both rNS-$\alpha$ and NS–Voigt. One observes that rNS-$\alpha$ is more sensitive than NS–Voigt to changes in $\alpha$. More interesting is the observation that the sensitivity is largest near the boundary. Although this is perhaps expected due to (1) the boundary layer is critical for turbulent channel flow; and (2) the ‘correct’ boundary conditions for filtered velocities is an open problem for LES/$\alpha$ models (note the N = 0 model does not use filtered velocity boundary conditions).

7. Conclusions

The main results to date for rNS-$\alpha$ have been reviewed in this paper. This includes global well-posedness, energy spectrum analysis, and energy dissipation analysis. An efficient numerical discretization is reviewed, and corresponding stability and convergence results are given. Additionally, sensitivity equations for rNS-$\alpha$ are derived, and a numerical scheme for them is given. Finally, results of several numerical tests are given, which show excellent results on several benchmark problems.

There are several directions for future work. Further testing, further comparison to other LES models, and exploration of the model outside of the finite element framework should all be done. One could also immediately move this model into the multiphysics arena, for non-isothermal flow, or doubly diffusive flow, or even to magnetohydrodynamics. Additional numerical testing on difficult benchmark problems is also important for any newer models. An interesting study could likely be done for rNS-$\alpha$’s near-wall behavior; it remains unclear what makes rNS-$\alpha$ so much more accurate for turbulent channel flow, but one hypothesis is that it could be more accurate near the wall. Another interesting future direction of this work is to consider the rNS-$\alpha$ model, but remove all modeling except in the viscous term. That is, one could consider a simplified model to be

$$\begin{array}{cc}\hfill w_t+(\nabla \times w)\times w+\nabla p-\nu \Delta D_{N}w& =f,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \nabla ·w& =0,\hfill \end{array}$$

(33)

or possibly keep the Voigt term $-{\alpha }^2\Delta w_t$. Proceeding as in Section 2, one can write the momentum equation of this reduced model as

$$w_t+(\nabla \times w)\times w+\nabla p-\nu \Delta w-\nu \Delta D_{N-1}(w-\overline{w})=f.$$

Note that this system can be considered as Navier–Stokes with an additional term meant to dissipate velocity fluctuations. Our discussion in Section 2 shows that the altered viscous term is strongly related to the VMM approach to turbulence modeling, so it is possible that it is only this term that is to thank for the model’s successes on the turbulent channel flow benchmark tests. One also notes that implementing the term would require very little extra work in a simulation, if the velocity term is treated implicitly and the filtered velocity is treated explicitly (by extrapolating from previous time steps).

Another possible direction for future work is choosing $\alpha$ locally and dynamically. Work in [45,46,47] has shown that carefully chosen local filtering radii can give large increases in accuracy in certain settings, and so extending these ideas to rNS-$\alpha$ seems like a logical direction.