Improving Accuracy in α-Models of Turbulence through Approximate Deconvolution

Abstract

In this report, we present several results in the theory of $\alpha$-models of turbulence with improved accuracy that have been developed in recent years. The $\alpha$-models considered herein are the Leray-$\alpha$ model, the zeroth Approximate Deconvolution Model (ADM) turbulence model, the modified Leray-$\alpha$ and the Navier–Stokes-$\alpha$ model. For all of the models from above, the accuracy is limited to ${\alpha }^2$ in smooth flow regions. Better accuracy requires decreasing the filter radius $\alpha$, which, in turn, requires a smaller mesh width that will lead in the end to a higher computational cost. Instead, one can use approximate deconvolution (without decreasing the mesh size) to attain better accuracy. Such deconvolution methods have been considered recently in many studies that show the efficiency of this approach. For smooth flows, periodic boundary conditions and van Cittert deconvolution operator of order N, the expected accuracy is ${\alpha }^{2N+2}$. In a bounded domain, such results are valid only in case special conditions are satisfied. In more general conditions, the author has recently proved that, in the case of the ADM, the expected accuracy of the finite element method with Taylor–Hood elements and Crank–Nicolson time stepping method is $\Delta t^2+h^2+K^N{\alpha }^2$, where the constant $K<1$ depends on the ratio $\alpha /h$, which is assumed constant. In this study, we present the extension of the result to the rest of the models.

1. Introduction

At a high Reynolds number, turbulence is not efficient to simulate using the Navier–Stokes equations because they require a very fine mesh and a high computational cost [1]. Therefore, in practical applications, reduced order turbulence models, such as the $\alpha$-models of turbulence discussed herein, have to be used.

The $\alpha$-models of turbulence are regularizations of the Navier–Stokes equations (NSE) whose goal is to allow stable computations on coarser grids than NSE. They are given by

$$\begin{array}{ccc}\hfill {\mathbf{w}}_t+N\left(\mathbf{w}\right)-\nu \Delta \mathbf{w}+\nabla p& =\mathbf{f}\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ \hfill \nabla ·\mathbf{w}& =0\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ \hfill \mathbf{w}(\mathbf{x},0)& ={\mathbf{w}}_0\left(\mathbf{x}\right)\hfill & \mathrm{in}\mathsf{\Omega }.\hfill \end{array}$$

(1)

Here, the nonlinearity N is $N\left(\mathbf{w}\right)=\overline{\mathbf{w}}·\nabla \mathbf{w}$ in the Leray-$\alpha$ model [2], $N\left(\mathbf{w}\right)=\mathbf{w}·\nabla \overline{\mathbf{w}}$ in the case of the modified Leray-$\alpha$ model [3], $N\left(\mathbf{w}\right)=\overline{\mathbf{w}}·\nabla \overline{\mathbf{w}}$ in the case of the zeroth Approximate Deconvolution Model (ADM) model [4] (here, the ADM is formulated as in [5]), and $N\left(\mathbf{w}\right)=(\nabla \times \mathbf{w})\times \overline{\mathbf{w}}$ in the case of the Navier–Stokes-$\alpha$ model.

In the above formula, the average $\overline{\mathbf{w}}$ is the solution of the BVP [6]

$$\begin{array}{ccc}\hfill \overline{\mathbf{w}}-{\alpha }^2\Delta \overline{\mathbf{w}}+\nabla p& =\hfill & \mathbf{w}\mathrm{in}\mathsf{\Omega },\\ \hfill \overline{\mathbf{w}}& =\hfill & \mathbf{w}\mathrm{on}\partial \mathsf{\Omega }.\end{array}$$

(2)

The study of $\alpha$-model of turbulence started with the work of Leray who introduced the Leray-$\alpha$ model in [7,8] (with the filter defined as the convolution with the Gaussian) and later studied in [2,9,10]. The other three models have been introduced and studied in papers such as [3,4,11,12,13,14,15,16,17,18,19,20].

One property that all four models have in common is that their accuracy $\left|\right|{\mathbf{u}}_{NSE}-{\mathbf{w}}_{\alpha -model}\left|\right|$ is limited to $\mathcal{O}\left({\alpha }^2\right)$ for smooth flows. In practice, this may lead to over-diffusivity and a lack of accuracy that requires smaller filter radius $\alpha$ and, as a consequence, a smaller mesh size that in turn causes higher computational costs. One way to correct this deficiency is to use an approximate deconvolution operator D having the approximation property $\mathbf{u}\approx D\overline{\mathbf{u}}$. One example of such deconvolution operator is the van Cittert deconvolution operator that approximates $\mathbf{u}$ by $D_N\overline{\mathbf{u}}$ defined as

$$D_N:=\sum _{n=0}^N{(I-G)}^n,N=0,1,2,\dots .$$

(3)

For $N=0,1,2$, the velocity field $\mathbf{u}$ will be approximated by

$$\begin{array}{cccc}\mathbf{u}\hfill & \approx \hfill & \overline{\mathbf{u}}\hfill & N=0\hfill \\ \mathbf{u}\hfill & \approx \hfill & 2\overline{\mathbf{u}}-\overline{\overline{\mathbf{u}}}\hfill & N=1\hfill \\ \mathbf{u}\hfill & \approx \hfill & 3\overline{\mathbf{u}}-3\overline{\overline{\mathbf{u}}}+\overline{\overline{\overline{\mathbf{u}}}}\hfill & N=2\hfill \end{array}$$

(4)

The finite element analysis of the van Cittert deconvolution method has been studied in [21]. The $\alpha$-models enhanced using approximate deconvolution become

$$\begin{array}{ccc}\hfill {\mathbf{w}}_t+DN\left(\mathbf{w}\right)-\nu \Delta \mathbf{w}+\nabla p& =\mathbf{f}\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ \hfill \nabla ·\mathbf{w}& =0\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ \hfill \mathbf{w}(\mathbf{x},0)& ={\mathbf{w}}_0\left(\mathbf{x}\right)\hfill & \mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(5)

where the nonlinearity $DN$ is $DN\left(\mathbf{w}\right)=D\overline{\mathbf{w}}·\nabla \mathbf{w}$ in the Leray-deconvolution model, $DN\left(\mathbf{w}\right)=\mathbf{w}·\nabla D\overline{\mathbf{w}}$ in the case of the modified Leray-deconvolution model, $DN\left(\mathbf{w}\right)=D\overline{\mathbf{w}}·\nabla D\overline{\mathbf{w}}$ in the case of the ADM model (formulated as in [21]) and $DN\left(\mathbf{w}\right)=(\nabla \times \mathbf{w})\times D\overline{\mathbf{w}}$ in the case of the enhanced Navier–Stokes-$\alpha$ model.

The Leray-deconvolution model has been introduced in [11] and later studied numerically in [22]. The ADM model (together with its zeroth version) has been introduced and studied in [13,14,23] and later studied in [20,24,25,26,27,28,29] (see also the monograph [30]). The deconvolution-enhanced Navier–Stokes-$\alpha$ model has been introduced and studied in [19] and also investigated numerically in [31,32,33]. To the author’s knowledge, the modified Leray-deconvolution model has not been studied so far.

Most of the analysis for the deconvolution-enhanced $\alpha$-models has been carried out in the context of van Cittert deconvolution. For periodic boundary conditions and assuming a smooth NSE solution, the modeling error is

$$\left|\right|{\mathbf{u}}_{NSE}-{\mathbf{w}\left|\right|}_{L^{\infty }\left(L^2\right)}+\left|\right|{\mathbf{u}}_{NSE}-\mathbf{w}{\left|\right|}_{L^2\left(H^1\right)}=\mathcal{O}\left({\alpha }^{2N+2}\right)$$

(6)

where $\mathbf{w}$ is the solution of the deconvolution enhanced Leray-$\alpha$ or the ADM or the deconvolution enhanced Navier–Stokes-$\alpha$ model [5,11,19,24,34]. Here, N is the order of the van Cittert deconvolution operator $D_N$. In a bounded domain, such estimates are valid only if special boundary conditions are satisfied by the exact NSE solution $\mathbf{u}$ [35].

The finite element analysis has been carried out for the Leray-deconvolution and the deconvolution enhanced Navier–Stokes-$\alpha$ in [19,22] with Crank–Nicolson time discretization and $P2/P1$ elements and the main error estimate is

$$\left|\right|{\mathbf{u}}_{NSE}-{\mathbf{w}}_h{\left|\right|}_{L^{\infty }\left(L^2\right)}+\left|\right|{\mathbf{u}}_{NSE}-{\mathbf{w}}_h{\left|\right|}_{L^2\left(H^1\right)}=\mathcal{O}(\Delta t^2+h^2+{\alpha }^{2N+2})$$

(7)

in the case of periodic boundary conditions. Here, the filter radius $\alpha$ satisfies $\alpha =\mathcal{O}\left(h\right)$. In the above formula, the term ${\alpha }^{2N+2}$ is induced by the deconvolution modeling. The case $N=0$ corresponds to the classical $\alpha$ turbulence models. In the case of a bounded domain, the formula above is no longer valid [35], and one can only prove that, for all N, the error satisfies

$$\left|\right|{\mathbf{u}}_{NSE}-{\mathbf{w}}_h{\left|\right|}_{L^{\infty }\left(L^2\right)}+\left|\right|{\mathbf{u}}_{NSE}-{\mathbf{w}}_h{\left|\right|}_{L^2\left(H^1\right)}=\mathcal{O}(\Delta t^2+h^2+{\alpha }^2)$$

(8)

Recently, in [29], the following estimate has been proved for the ADM:

$$\left|\right|{\mathbf{u}}_{NSE}-{\mathbf{w}}_h{\left|\right|}_{L^{\infty }\left(L^2\right)}+\left|\right|{\mathbf{u}}_{NSE}-{\mathbf{w}}_h{\left|\right|}_{L^2\left(H^1\right)}=\mathcal{O}(\Delta t^2+h^2+K^N{\alpha }^2)$$

(9)

In the above estimate, N is the order of the van Cittert deconvolution operator and $K<1$ is a constant that only depends on the ratio $\alpha /h$ (which is assumed constant) and the sequence of quasiuniform meshes used in the computation. The term $K^N{\alpha }^2$ is due to the deconvolution and it can be made small by increasing N and decreasing the ratio $\alpha /h$ [21]. This estimate supports the claim that high order deconvolution operators indeed improve accuracy, a behavior that has been observed in the numerical tests in [22,29].

We believe that this is the general behavior of deconvolution enhanced $\alpha$-models of turbulence in case the van Cittert deconvolution procedure is used. Using similar techniques such as those in [29], similar estimates can be obtained for the Leray-deconvolution model, the Navier–Stokes $\alpha$ deconvolution model and the modified Leray-$\alpha$ deconvolution model.

2. Mathematical Context

The mathematical notations and concepts are similar to the ones used in [19,22,29]. $\mathsf{\Omega }\subset {\mathbb{R}}^d$, $d=2,3$, is regular, bounded, polyhedral domain and $\parallel ·\parallel$ and $(·,·)$ is the $L^2\left(\mathsf{\Omega }\right)$ norm and inner product. Lebesgue spaces and their norms are denoted by $L^p\left(\mathsf{\Omega }\right)$, ${\parallel ·\parallel }_{L^p}$. We use standard notations for the Lebesque and Sobolev spaces and their norms.

In the variational problem, we use the functional spaces [36]:

$$\begin{array}{c}\mathrm{Velocity}\mathrm{Space}-X:={\left[H_0^1\left(\mathsf{\Omega }\right)\right]}^d\\ \mathrm{Pressure}\mathrm{Space}-Q:=L_0^2\left(\mathsf{\Omega }\right)=\left\{q\in L^2\left(\mathsf{\Omega }\right):{\int }_{\mathsf{\Omega }}qdx=0\right\}\\ \mathrm{Divergence}-\mathrm{free}\mathrm{Space}-V:=\left\{v\in X:{\int }_{\mathsf{\Omega }}q\nabla ·\mathbf{v}dx=0,\forall q\in Q\right\}\end{array}$$

(10)

The finite element analysis is carried out on a family of triangulations ${\left({\mathcal{T}}_h\right)}_h$ on $\mathsf{\Omega }$ that are quasiuniform [37]. We also consider an inf-sup stable pair of finite elements $(X_h,Q_h)$ where $X_h\subset X$, $Q_h\subset Q$, which are assumed to satisfy [37]:

$$\begin{array}{c}\parallel \mathbf{u}-I{\mathbf{u}}_h\parallel \le C{h}^{k+1}{\parallel \mathbf{u}\parallel }_{k+1},\mathbf{u}\in H^{k+1}{\left(\mathsf{\Omega }\right)}^d\\ \parallel \mathbf{u}-I{\mathbf{u}}_h{\parallel }_1\le C{h}^k{\parallel \mathbf{u}\parallel }_{k+1},\mathbf{u}\in H^{k+1}{\left(\mathsf{\Omega }\right)}^d\\ \parallel p-I{p}_h\parallel \le C{h}^k{\parallel p\parallel }_k,p\in H^k\left(\mathsf{\Omega }\right)\end{array}$$

(11)

where $I{\mathbf{u}}_h\in X_h$ is an interpolant of $\mathbf{u}$ and $I{p}_h\in Q_h$ is an interpolant of p.

The discretely divergence-free space $V_h$ is

$$\begin{array}{c}{V}_h:=\left\{\mathbf{v}\in X_h:(q,\nabla ·\mathbf{v})=0,\forall q\in Q_h\right\}\end{array}$$

(12)

We assume that $V_h$ satisfies the first two approximation properties in the inequalities (11) above if $\mathbf{u}\in V$.

We let $t_n=n\Delta t$, $t_{n+1/2}=(n+1/2)\Delta t$, $T:=N_T\Delta t$. Here, $\Delta t$ is the chosen time step, $N_T$ is the number of time steps, and T is the final time. We will also use the notations ${\mathbf{v}}^n=\mathbf{v}(n\Delta t)$ and ${\mathbf{v}}^{n+1/2}=({\mathbf{v}}^n+{\mathbf{v}}^{n+1})/2$ for a function $\mathbf{v}$.

We will assume that the $L^2$ projection operator onto $V_h$ denoted by $P_{V_h}:L^2\left(\mathsf{\Omega }\right)\to V_h\subset L^2\left(\mathsf{\Omega }\right)$ is $H^1$ stable. The analysis also use the inverse inequality [37]:

Lemma 1.

$$\left|\right|\nabla {\mathbf{v}}_h\left|\right|\le \frac{C_{\mathcal{T}}}{h}\left|\right|{\mathbf{v}}_h\left|\right|,\forall {\mathbf{v}}_h\in X^h$$

(13)

The constant $C_{\mathcal{T}}$ does not depend on h.

Definition 1.

[30] The discrete Stokes operator $A_h:V_h\to V_h$ satisfies $(A_{h}v,v_h)=(\nabla v,\nabla v_h)$ for any $v_h\in V_h$.

One may easily show that $A_h$ is a bijective, self-adjoint, positive operator with respect to the inner products $(·,·)$ and $(\nabla ·,\nabla ·)$ on $V_h$.

The eigenvalues of $A_h$ are denoted by ${\lambda }_1<{\lambda }_2<\cdots <{\lambda }_M.$

One may show using the inverse inequality (13) that

$${\lambda }_M\le {\displaystyle \frac{C_{\mathcal{T}}^2}{h^2}},\left|\right|A_h\mathbf{v}\left|\right|\le {\displaystyle \frac{C_{\mathcal{T}}}{h}}\left|\right|\nabla \mathbf{v}\left|\right|,\left|\right|A_h\mathbf{v}\left|\right|\le \frac{C_{\mathcal{T}}^2}{h^2}\left|\right|\mathbf{v}\left|\right|$$

(14)

for $\mathbf{v}\in V_h$.

Definition 2.

(Discrete Stokes average [30]) For given $\alpha >0$, we let $G:V_h\to V_h$,

$$G\mathbf{v}={\overline{\mathbf{v}}}^h={(I+{\alpha }^2{A}_h)}^{-1}\mathbf{v}$$

(15)

When $\mathbf{v}\in X$, we set ${\overline{\mathbf{v}}}^h:=G{P}_{V_h}\mathbf{v}$.

Theorem 1.

[29] Due to the spectral mapping theorem, the eigenvalues of the discrete Stokes filter G are ${\displaystyle \frac{1}{1+{\alpha }^2{\lambda }_i}}$, $i=1,2,\cdots ,M$, and they satisfy

$${\displaystyle \frac{1}{1+C_{\mathcal{T}}^2{\displaystyle \frac{{\alpha }^2}{h^2}}}}\le {\displaystyle \frac{1}{1+{\alpha }^2{\lambda }_M}}<{\displaystyle \frac{1}{1+{\alpha }^2{\lambda }_2}}<\cdots <{\displaystyle \frac{1}{1+{\alpha }^2{\lambda }_1}}<1$$

(16)

Therefore, $\left|\right|{\overline{\mathbf{v}}}^h\left|\right|\le \left|\right|\mathbf{v}\left|\right|$, $\left|\right|\nabla {\overline{\mathbf{v}}}^h\left|\right|\le \left|\right|\nabla \mathbf{v}\left|\right|$ for all $\mathbf{v}\in V_h$.

One may also notice that the largest eigenvalue of $I-G$ is bounded by

$$K={\displaystyle \frac{C_{\mathcal{T}}^2{\displaystyle \frac{{\alpha }^2}{h^2}}}{1+C_{\mathcal{T}}^2{\displaystyle \frac{{\alpha }^2}{h^2}}}}<1$$

(17)

Therefore,

$$\left|\right|I-G\left|\right|\le K$$

(18)

Definition 3.

[30] The N-th order discrete van Cittert deconvolution operator $D_N:V_h\to V_h$ is given by:

$$D_N:=\sum _{n=0}^N{(I-G)}^n,N=0,1,2,\dots .$$

(19)

Several properties of $D_N$ that can be found in [30] are listed below. $D_N$ is self-adjoint, positive with respect to the inner products $(·,·)$ and $(\nabla ·,\nabla ·)$. $D_{N}G=I-{(I-G)}^{N+1};$ therefore, the eigenvalues of $D_{N}G$ are $1-{\left({\displaystyle \frac{{\alpha }^2{\lambda }_i}{1+{\alpha }^2{\lambda }_i}}\right)}^{N+1}$, $i=1\cdots M$, and can be bounded as [29]:

$$1-{\left({\displaystyle \frac{C_{\mathcal{T}}^2{\displaystyle \frac{{\alpha }^2}{h^2}}}{1+C_{\mathcal{T}}^2{\displaystyle \frac{{\alpha }^2}{h^2}}}}\right)}^{N+1}<1-{\left({\displaystyle \frac{{\alpha }^2{\lambda }_M}{1+{\alpha }^2{\lambda }_M}}\right)}^{N+1}<\cdots <1-{\left({\displaystyle \frac{{\alpha }^2{\lambda }_1}{1+{\alpha }^2{\lambda }_1}}\right)}^{N+1}<1$$

(20)

Following [19], we define two norms on $V_h$

$${\left|\right|\mathbf{v}\left|\right|}_E:={(\mathbf{v},D_{N}G\mathbf{v})}^{1/2}{,\left|\right|\mathbf{v}\left|\right|}_{ϵ}:={(\nabla \mathbf{v},\nabla D_{N}G\mathbf{v})}^{1/2}$$

(21)

Using inequalities (20), it follows that, see also for more details [29],

$$\begin{array}{ccccc}\sqrt{1-K^{N+1}}\left|\right|\mathbf{v}\left|\right|\hfill & \le \hfill & {\left|\right|\mathbf{v}\left|\right|}_E\hfill & \le \hfill & \left|\right|\mathbf{v}\left|\right|\hfill \\ \sqrt{1-K^{N+1}}\left|\right|\nabla \mathbf{v}\left|\right|\hfill & \le \hfill & {\left|\right|\mathbf{v}\left|\right|}_{ϵ}\hfill & \le \hfill & \left|\right|\nabla \mathbf{v}\left|\right|\hfill \end{array}$$

(22)

and

$$\begin{array}{c}\left|\right|D_{N}G\mathbf{v}\left|\right|\le {\left|\right|\mathbf{v}\left|\right|}_E,\left|\right|\nabla D_{N}G\mathbf{v}\left|\right|\le {\left|\right|\mathbf{v}\left|\right|}_{ϵ},\forall \mathbf{v}\in V_h\hfill \\ \left|\right|D_{N}G{P}_{V_h}\mathbf{v}\left|\right|\le \left|\right|\mathbf{v}\left|\right|,\left|\right|\nabla D_{N}G{P}_{V_h}\mathbf{v}\left|\right|\le C\left|\right|\nabla \mathbf{v}\left|\right|,\forall \mathbf{v}\in X\hfill \end{array}$$

(23)

We will also use the next estimate of the discrete deconvolution error $\mathbf{v}-D_N{\overline{\mathbf{v}}}^h$ that has been obtained in [21].

Theorem 2.

[21] If $\mathbf{v}\in H^{k+1}\left(\mathsf{\Omega }\right)\cap V$, $k\ge 1$ and $\alpha =\mathcal{O}\left(h\right),$ then

$$\begin{array}{ccc}\left|\right|\mathbf{v}-D_N{\overline{\mathbf{v}}}^h\left|\right|\hfill & \le \hfill & C\left(h^2{K}^N\right||\Delta \mathbf{v}||+h^{k+1}|\mathbf{v}{|}_{k+1})\hfill \\ \left|\right|\nabla \mathbf{v}-\nabla D_N{\overline{\mathbf{v}}}^h\left|\right|\hfill & \le \hfill & C\left(h{K}^N\right||\Delta \mathbf{v}||+h^k|\mathbf{v}{|}_{k+1})\hfill \end{array}$$

(24)

Here, N is the order of the discrete deconvolution operator and $K<1$ has been previously defined in (17).

Proof. 

We recall the main ideas of the proof presented in [21]. First, we split the error $\left|\right|\mathbf{v}-D_N{\overline{\mathbf{v}}}^h\left|\right|$ using the triangle inequality

$$\left|\right|\mathbf{v}-D_N{\overline{\mathbf{v}}}^h\left|\right|\le \left|\right|\mathbf{v}-P_{V_h}\mathbf{v}\left|\right|+\left|\right|P_{V_h}\mathbf{v}-D_N{\overline{\mathbf{v}}}^h\left|\right|$$

(25)

The term $\left|\right|\mathbf{v}-P_{V_h}\mathbf{v}\left|\right|$ is bounded by $C{h}^{k+1}{\left|\mathbf{v}\right|}_{k+1}$ due to the approximation properties of $V_h$. It remains then to estimate $\left|\right|P_{V_h}\mathbf{v}-D_N{\overline{\mathbf{v}}}^h\left|\right|$. Since $D_{N}G=I-{(I-G)}^{N+1},$ this term is written as

$$\left|\right|P_{V_h}\mathbf{v}-D_N{\overline{\mathbf{v}}}^h\left|\right|={\left|\right|(I-G)}^{N+1}{P}_{V_h}\mathbf{v}\left|\right|$$

(26)

It follows that

$$\left[\right||P_{V_h}\mathbf{v}-D_N{\overline{\mathbf{v}}}^h||\le ||I-G|{|}^N\left|\right|(I-G)P_{V_h}\mathbf{v}\left|\right|\le K^N{\alpha }^2\left|\right|A_h{\overline{\mathbf{v}}}^h\left|\right|$$

(27)

The discrete Laplacian $\left|\right|A_h{\overline{\mathbf{v}}}^h\left|\right|$ can be estimated as

$$\left|\right|A_h{\overline{\mathbf{v}}}^h\left|\right|=\left|\right|{\overline{A_h{P}_{V_h}\mathbf{v}}}^h\left|\right|\le C\left|\right|A_h{P}_{V_h}\mathbf{v}\left|\right|\le C\left|\right|\Delta \mathbf{v}\left|\right|$$

(28)

The last inequality from above can be found in [38]. It follows that

$$\left|\right|P_{V_h}\mathbf{v}-D_N{\overline{\mathbf{v}}}^h\left|\right|\le C{K}^N{\alpha }^2\left|\right|\Delta \mathbf{v}\left|\right|$$

(29)

This proves the first inequality in inequalities (24).

The second inequality in inequalities (24) is immediatelly proved using the same arguments combined with inverse inequalities.

 ☐

3. Numerical Scheme for the $\alpha$-Models of Turbulence

In our analysis $b^{*}(·,·,·):X\times X\times X\to \mathrm{I}\mathrm{R}$ will denote the standard skew-symmetrized trilinear form [36]:

$$b^{*}(\mathbf{u},\mathbf{v},\mathbf{w}):=\frac{1}{2}(\mathbf{u}·\nabla \mathbf{v},\mathbf{w})-\frac{1}{2}(\mathbf{u}·\nabla \mathbf{w},\mathbf{v})$$

(30)

One can show that the trilinear form satisfies [39]:

$$|b^{*}(\mathbf{u},\mathbf{v},\mathbf{w})|+|b^{*}{(\mathbf{v},\mathbf{u},\mathbf{w})|\le C|\left|\mathbf{u}\right|\left|\right(\left|\right|\mathbf{v}\left|\right|}_{L^{\infty }}+{\left|\right|\nabla \mathbf{v}\left|\right|}_{L^3}\left)\right||\nabla \mathbf{w}||$$

(31)

It follows then using inequalities (35) and (36) in [38] that

$$|b^{*}(\mathbf{u},\mathbf{v},\mathbf{w})|+|b^{*}(\mathbf{v},\mathbf{u},\mathbf{w})|\le C|\left|\mathbf{u}\right|\left|\right||A_h\mathbf{v}\left|\right|\left|\right|\nabla \mathbf{w}\left|\right|$$

(32)

for $\mathbf{u},\mathbf{w}\in X$, $\mathbf{v}\in V_h$.

In the case of the Navier–Stokes in rotational form, the trilinear form is

$$\tilde{b}(\mathbf{u},\mathbf{v},\mathbf{w})=((\nabla \times \mathbf{v})\times \mathbf{u},\mathbf{w})$$

(33)

For body force $\mathbf{f}\in L^2((0,T],X)$ and initial velocity ${\mathbf{u}}_0\in X$, the discrete solution ${\mathbf{w}}_h^{n+1}\in V_h$ of the accuracy enhanced $\alpha$ model of turbulence at step $n+1$ for $n=0,1,2,\cdots ,N_T-1$ satisfies:

$$\begin{array}{ccc}\hfill ({\displaystyle \frac{{\mathbf{w}}_h^{n+1}-{\mathbf{w}}_h^n}{\Delta t}},\mathbf{v})& +\hfill & B({\mathbf{w}}_h^{n+1/2},\mathbf{v})-(p_h^{n+1/2},\nabla ·\mathbf{v})\hfill \\ \hfill +\nu (\nabla {\mathbf{w}}_h^{n+1/2},\nabla \mathbf{v})& =\hfill & (\mathbf{f}((n+1/2)\Delta t),\mathbf{v}),\forall \mathbf{v}\in X_h\hfill \\ \hfill (q,\nabla ·{\mathbf{w}}_h^{n+1})& =\hfill & 0,\forall q\in P_h\hfill \\ \hfill ({\mathbf{w}}_h^0,\mathbf{v})& =\hfill & ({\mathbf{w}}_0,\mathbf{v}),\forall \mathbf{v}\in X_h\hfill \end{array}$$

(34)

where

$$\begin{array}{c}B({\mathbf{w}}_h^{n+1/2},\mathbf{v})=b^{*}(D_N{\overline{{\mathbf{w}}_h^{n+1/2}}}^h,{\mathbf{w}}_h^{n+1/2},\mathbf{v}),\mathrm{Leray}-\mathrm{deconvolution}\mathrm{model} [22],\hfill \\ B({\mathbf{w}}_h^{n+1/2},\mathbf{v})=b^{*}({\mathbf{w}}_h^{n+1/2},D_N{\overline{{\mathbf{w}}_h^{n+1/2}}}^h,\mathbf{v}),\mathrm{modified}\mathrm{Leray}-\mathrm{deconvolution}\mathrm{model},\hfill \\ B({\mathbf{w}}_h^{n+1/2},\mathbf{v})=b^{*}(D_N{\overline{{\mathbf{w}}_h^{n+1/2}}}^h,D_N{\overline{{\mathbf{w}}_h^{n+1/2}}}^h,\mathbf{v}),\mathrm{ADM} [\mathrm{5,29}],\hfill \\ B({\mathbf{w}}_h^{n+1/2},\mathbf{v})=\tilde{b}(D_N{\overline{{\mathbf{w}}_h^{n+1/2}}}^h,{\mathbf{w}}_h^{n+1/2},\mathbf{v}),\mathrm{NS}-\alpha \mathrm{deconvolution}\mathrm{model} [19].\hfill \end{array}$$

(35)

Thorough the analysis, we assume that $\alpha /h$ is constant.

Lemma 2.

(Stability) A finite element solution ${\mathbf{w}}_h^{n+1}$ of the model (34) exists at each time step $n=0,\dots N_T-1$ and satisfies the stability estimate

$$\begin{array}{c}\underset{0\le n<N_T}{sup}\left|\right|{\mathbf{w}}_h^{n+1}{\left|\right|}^2+\Delta t\nu {\displaystyle {\sum }_{n=0}^{N_T-1}}\left|\right|\nabla {\mathbf{w}}_h^{n+1/2}{\left|\right|}^2\hfill \\ \hfill \le C\left(\frac{\Delta t}{\nu }{\displaystyle \sum _{n=0}^{N_T-1}}\left|\right|f\left(t^{n+1/2}\right){\left|\right|}_{-1}^2+\left|\right|{\mathbf{w}}_0{\left|\right|}^2\right)\end{array}$$

(36)

where C does not depend on $h,\alpha$ (but it depends on $\alpha /h$).

Proof. 

The Leray-deconvolution model case has been studied (for a slightly different filter operator, but the argument is still valid) in [22]. The NS-$\alpha$-deconvolution case in [19]. The modified Leray-deconvolution case can be proved using arguments similar to the NS $\alpha$-deconvolution (since it requires a similar multiplier to cancel the nonlinearity).

The ADM case has been studied in [29] and shares great similarity with the stability proofs for the Leray-deconvolution and the NS-$\alpha$-deconvolution models mentioned previously. We outline the ideas from [29] here.

To prove the existence of a discrete solution ${\mathbf{u}}_h^{n+1}\in V_h$, we consider the operator $L:V_h\to V_h$, $L\phi =\psi$ such that $\phi ,\psi$ satisfy the equality

$$\begin{array}{c}2{(\psi ,\mathbf{v})}_E+\Delta t\nu {(\psi ,\mathbf{v})}_{ϵ}=2{({\mathbf{u}}_h^n,\mathbf{v})}_E-\Delta t{b}^{*}(D_N{\overline{\phi }}^h,D_N{\overline{\phi }}^h,D_{N}G\mathbf{v})\hfill \\ +\Delta t(\mathbf{f}\left(t^{n+1/2}\right),D_{N}G\mathbf{v})\forall \mathbf{v}\in V_h\hfill \end{array}$$

(37)

We notice that, if $\phi$ is a fixed point of L, then ${\mathbf{u}}_h^{n+1}=2\phi -{\mathbf{u}}_h^n$ is a solution of the ADM model (34).

In order to find a fixed point of $\phi$, we apply the Leray–Schauder fixed point theorem. To show that L is a compact operator, we notice that it is the composition of a continuous, linear operator $T:V_h\to V_h$, $T\zeta =\psi$, where

$$\begin{array}{c}2{(\psi ,\mathbf{v})}_E+\Delta t\nu {(\psi ,\mathbf{v})}_{ϵ}=(\zeta ,\mathbf{v})\forall \mathbf{v}\in V_h\hfill \end{array}$$

(38)

and a nonlinear, continuous and bounded operator $N:V_h\to V_h$, $N\phi =\zeta$, where

$$(\zeta ,\mathbf{v})=2{({\mathbf{u}}_h^n,\mathbf{v})}_E-\Delta t{b}^{*}(D_N{\overline{\phi }}^h,D_N{\overline{\phi }}^h,D_{N}G\mathbf{v})+\Delta t(\mathbf{f}\left(t^{n+1/2}\right),D_{N}G\mathbf{v})\forall \mathbf{v}\in V_h$$

(39)

To show that N is bounded and continuous, we use the equivalence of all norms on the finite element space $V_h$, similar to the proof in [19].

It remains only to show that the set

$$\{\phi \in V_h|\phi =\lambda L\phi \mathrm{for}\mathrm{some}\lambda \in [0,1)\}$$

(40)

is bounded in the $L^2$ norm independent of $\lambda$.

Assume $\phi =\lambda L\phi \in V_h$, i.e.,

$$\begin{array}{c}2{(\phi ,\mathbf{v})}_E+\Delta t\nu {(\phi ,\mathbf{v})}_{ϵ}=\lambda (2{({\mathbf{u}}_h^n,\mathbf{v})}_E-\Delta t{b}^{*}(D_N{\overline{\phi }}^h,D_N{\overline{\phi }}^h,D_{N}G\mathbf{v})\hfill \\ +\Delta t(\mathbf{f}\left(t^{n+1/2}\right),D_{N}G\mathbf{v}))\forall \mathbf{v}\in V_h\hfill \end{array}$$

(41)

and set $\mathbf{v}=\phi$, cancel the nonlinearity and get

$${2\left|\right|\phi \left|\right|}_E^2+{\Delta t\nu \left|\right|\phi \left|\right|}_{ϵ}^2=\lambda (2{({\mathbf{u}}_h^n,\phi )}_E+\Delta t(\mathbf{f}\left(t^{n+1/2}\right),D_{N}G\phi ))\forall \mathbf{v}\in V_h$$

(42)

Next, using the standard Cauchy–Schwartz and Yound inequalities on the right terms, we get that

$${\left|\right|\phi \left|\right|}_E^2+{\Delta t\nu \left|\right|\phi \left|\right|}_{ϵ}^2\le C\left(\right||\Delta t{\nu }^{-1}\mathbf{f}\left(t^{n+1/2}\right){\left|\right|}_{*}^2+\left|\right|{\mathbf{u}}_h^n{\left|\right|}_E^2)$$

(43)

Therefore, by the Leray–Schauder Theorem, it follows that the operator L has a fixed point $\phi ={\mathbf{u}}_h^{n+1/2}$ and the discrete solution ${\mathbf{u}}_h^{n+1}$ of the ADM system exists.

To obtain the stability estimate from above, we set $D_N{\overline{{\mathbf{u}}_h^{n+1/2}}}^h$ as a test function in the ADM model (34) and the nonlinear term will vanish:

$${\displaystyle \frac{1}{2\Delta t}}\left|\right|({\mathbf{u}}_h^{n+1}{\left|\right|}_E^2-\left|\right|{\mathbf{u}}_h^n{\left|\right|}_E^2)+\nu ||{\mathbf{u}}_h^{n+1/2}{\left|\right|}_{ϵ}^2=(\mathbf{f}\left(t^{n+1/2}\right),D_{N}G{\mathbf{u}}_h^{n+1/2}))$$

(44)

Using the standard Cauchy–Schwartz and Young inequalities on the right term, we get

$${\displaystyle \frac{1}{\Delta t}}\left|\right|({\mathbf{u}}_h^{n+1}{\left|\right|}_E^2-\left|\right|{\mathbf{u}}_h^n{\left|\right|}_E^2)+\nu ||{\mathbf{u}}_h^{n+1/2}{\left|\right|}_{ϵ}^2\le C{\nu }^{-1}\left|\right|\mathbf{f}\left(t^{n+1/2}\right){\left|\right|}_{*}^2$$

(45)

Summing up from $n=0$ to $n=N_T-1$, we get the required inequality

$$\begin{array}{c}\underset{0\le n<N_T}{sup}\left|\right|{\mathbf{w}}_h^{n+1}{\left|\right|}_E^2+\Delta t\nu {\displaystyle {\sum }_{n=0}^{N_T-1}}\left|\right|\nabla {\mathbf{w}}_h^{n+1/2}{\left|\right|}_{ϵ}^2\hfill \\ \hfill \le C\left(\frac{\Delta t}{\nu }{\displaystyle \sum _{n=0}^{N_T-1}}\left|\right|f\left(t^{n+1/2}\right){\left|\right|}_{-1}^2+\left|\right|{\mathbf{w}}_0{\left|\right|}_E^2\right)\end{array}$$

(46)

Here, $C=C\left(\mathsf{\Omega }\right)$.

The norms $\left|\right|·{\left|\right|}_E,\left|\right|·{\left|\right|}_{ϵ}$ can be replaced by the standard $L^2$, respectively $H^1$ norm via the inequality (22), and this will add a dependence of C on K, (see inequality 17), i.e., on the ratio $\alpha /h$.

 ☐

The convergence of the discrete solution ${\mathbf{w}}^h$ of the enhanced $\alpha$-model (34) to the NSE exact solution $\mathbf{u}$ is stated in the next theorem.

Theorem 3.

We let ($\mathbf{u}$,p) be a smooth strong solution of the NSE and ${\mathbf{w}}_h^n$, $n=1\dots N_T$ be the discrete solution of (34). We assume $\alpha /h$ is constant and ${\mathbf{w}}_0={\mathbf{u}}_0$. Then, if $\Delta t$ is picked small enough, there holds

$$\begin{array}{c}\underset{0\le n<N_T}{sup}\left|\right|{\mathbf{u}}^{n+1}-{\mathbf{w}}_h^{n+1}{\left|\right|}^2+\Delta t\nu {\displaystyle {\sum }_{n=0}^{N_T-1}}\left|\right|\nabla {\mathbf{u}}^{n+1/2}-\nabla {\mathbf{w}}_h^{n+1/2}{\left|\right|}^2\hfill \\ \hfill \le C\left(h^{2k}+\Delta t^4+K^{2N}{\alpha }^4\right)\end{array}$$

(47)

where C is a constant that does not depend on $h,\alpha$ (though it depends on $\alpha /h$) and K is given in inequality (17) .

Remark 1.

The proof uses the discrete Gronwall’s inequality [38], which requires that $\Delta t$ satisfies $\Delta t=\mathcal{O}\left(\right({\nu }^{-3}{\left|\right||\nabla \mathbf{u}|\left|\right|}_{\infty ,0}+1{)}^{-1})$.

Proof. 

The convergence proof for the Leray-deconvolution model is done in [22] and for the NS-$\alpha$ deconvolution model in [19], but the two proofs do not use the improved deconvolution error estimate (24) and the error is proved to be in the bounded domain case on the order of

$$\mathcal{O}\left(h^{2k}+\Delta t^4+{\alpha }^4\right)$$

(48)

In the two proofs in [19,22], one bounds

$$\begin{array}{c}\underset{0\le n<N_T}{sup}\left|\right|{\mathbf{u}}^{n+1}-{\mathbf{w}}_h^{n+1}{\left|\right|}^2+\Delta t\nu {\displaystyle {\sum }_{n=0}^{N_T-1}}\left|\right|\nabla {\mathbf{u}}^{n+1/2}-\nabla {\mathbf{w}}_h^{n+1/2}{\left|\right|}^2\hfill \end{array}$$

(49)

by several terms which can be left unchanged herein and also by the modelling error ( the Intp term in [22], the Interp term in [19]), which itself contains several terms that can be left unchanged but also a term that is bounded by the deconvolution error $\left|\right|{\mathbf{u}}^{n+1/2}-D_N{\mathbf{u}}^{n+1/2}\left|\right|$. Applying the improved deconvolution estimate (24) provides the improved $K^{2N}{\alpha }^4$ instead of ${\alpha }^4$.

The proof of the above result for the ADM is done in [29]. We outline here the main ideas of the proof presented in [29], which are very similar to the ones in [22] except that the nonlinearity is handled differently.

Similar to [19,22], the NSE is written in the form

$$\begin{array}{c}({\displaystyle \frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\Delta t}},\mathbf{v})+b^{*}(D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,\mathbf{v})+\nu (\nabla {\mathbf{u}}^{n+1/2},\nabla \mathbf{v})\hfill \\ -(p^{n+1/2},\nabla ·\mathbf{v})=(\mathbf{f}((n+1/2)\Delta t),\mathbf{v})-Interp(\mathbf{u},p,n,\mathbf{v})\forall \mathbf{v}\in X_h\hfill \end{array}$$

(50)

where

$$\begin{array}{c}Interp(\mathbf{u},p,n,\mathbf{v})=({\mathbf{u}}_t\left(t^{n+1/2}\right)-{\displaystyle \frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\Delta t}},\mathbf{v})+b^{*}(\mathbf{u},\mathbf{u},\mathbf{v})\hfill \\ -b^{*}(D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,\mathbf{v})+\nu (\nabla \mathbf{u}\left(t^{n+1/2}\right)-\nabla {\mathbf{u}}^{n+1/2},\nabla \mathbf{v})\hfill \\ +(p\left(t^{n+1/2}\right)-p^{n+1/2},\nabla ·\mathbf{v})\hfill \end{array}$$

(51)

Next, one derives the error equation for ${\mathbf{e}}^{n+1}={\mathbf{u}}^{n+1}-{\mathbf{u}}_h^{n+1}$ by subtracting (34) from (50). We let ${\mathbf{U}}^{n+1}\in V_h$ be the $L^2$ projection of ${\mathbf{u}}^{n+1}$ onto $V_h$, split the error as ${\mathbf{e}}^{n+1}=({\mathbf{u}}^{n+1}-{\mathbf{U}}^{n+1})-({\mathbf{u}}_h^{n+1}-{\mathbf{U}}^{n+1}):={\eta }^{n+1}-{\phi }_h^{n+1}$ and set $\mathbf{v}=D_N{\overline{{\phi }_h^{n+1/2}}}^h$ in the error equation. Upon rearranging terms, we obtain

$$\begin{array}{c}{\displaystyle \frac{1}{2}}\left(\right||{\phi }_h^{n+1}{\left|\right|}_E^2-\left|\right|{\phi }_h^n{\left|\right|}_E^2+\nu \Delta t\left|\right|{\phi }_h^{n+1/2}{\left|\right|}_{ϵ}^2=({\eta }^{n+1}-{\eta }^n,D_N{\overline{{\phi }_h^{n+1/2}}}^h)+\hfill \\ \nu \Delta t(\nabla {\eta }^{n+1/2},\nabla D_N{\overline{{\phi }_h^{n+1/2}}}^h)-\Delta t{b}^{*}(D_N{\overline{{\eta }^{n+1/2}}}^h,D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\phi }_h^{n+1/2}}}^h)\hfill \\ +\Delta t{b}^{*}(D_N{\overline{{\phi }_h^{n+1/2}}}^h,D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\phi }_h^{n+1/2}}}^h)\hfill \\ -\Delta t{b}^{*}(D_N{\overline{{\mathbf{u}}_h^{n+1/2}}}^h,D_N{\overline{{\eta }^{n+1/2}}}^h,D_N{\overline{{\phi }_h^{n+1/2}}}^h)\hfill \\ +\Delta t(p^{n+1/2},\nabla ·D_N{\overline{{\phi }_h^{n+1/2}}}^h)+\Delta tInterp(\mathbf{u},p,n,D_N{\overline{{\phi }_h^{n+1/2}}}^h)\hfill \end{array}$$

(52)

On the right side, all above terms excepting the two trilinear terms appearing in the interpolating term (51) will be estimated exactly as in [22]. In all of these terms, wherever needed (to get rid of $D_{N}G$ terms and switch to the norms $\left|\right|·{\left|\right|}_E,\left|\right|·{\left|\right|}_{ϵ}$) to keep the arguments similar to the ones in [22], inequalities (23) will be used.

The two trilinear terms in $Interp$ that are treated differently compared to [22] take the form

$$\begin{array}{c}{b}^{*}({\mathbf{u}}^{n+1/2},{\mathbf{u}}^{n+1/2},D_N{\overline{{\phi }_h^{n+1/2}}}^h)-b^{*}(D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\phi }_h^{n+1/2}}}^h)=\hfill \\ b^{*}({\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,{\mathbf{u}}^{n+1/2},D_N{\overline{{\phi }_h^{n+1/2}}}^h)\hfill \\ +b^{*}(D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,{\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\phi }_h^{n+1/2}}}^h)\hfill \end{array}$$

(53)

The first term on the right side above is estimated similar to inequality 3.26 in [22]:

$$\begin{array}{c}|b^{*}({\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,{\mathbf{u}}^{n+1/2},D_N{\overline{{\phi }_h^{n+1/2}}}^h)|\hfill \\ \le {\displaystyle \frac{\nu }{24}}\left|\right|{\phi }_h^{n+1/2}{\left|\right|}_{ϵ}^2+C{\nu }^{-1}\left|\right|{\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h{\left|\right|}^2\hfill \end{array}$$

(54)

For the second term, we first use the inequality (32) to get

$$\begin{array}{c}|b^{*}(D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,{\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\phi }_h^{n+1/2}}}^h)|\hfill \\ \le C\left|\right|{\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h\left|\right|\left|\right|A_h{D}_N{\overline{{\mathbf{u}}^{n+1/2}}}^h\left|\right|\left|\right|\nabla D_N{\overline{{\phi }_h^{n+1/2}}}^h\left|\right|\hfill \end{array}$$

(55)

Using again the argument in inequality (28) to estimate the middle term and Young’s inequality, we get further

$$\begin{array}{c}C\left|\right|{\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h\left|\right|\left|\right|A_h{D}_N{\overline{{\mathbf{u}}^{n+1/2}}}^h\left|\right|\left|\right|\nabla D_N{\overline{{\phi }_h^{n+1/2}}}^h\left|\right|\hfill \\ \le C\left|\right|{\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h\left|\right|\left|\right|\Delta {\mathbf{u}}^{n+1/2}\left|\right|\left|\right|{\phi }_h^{n+1/2}{\left|\right|}_{ϵ}\hfill \\ \le C{\nu }^{-1}\left|\right|{\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h{\left|\right|}^2+{\displaystyle \frac{\nu }{24}}\left|\right|{\phi }_h^{n+1/2}{\left|\right|}_{ϵ}^2\hfill \end{array}$$

(56)

Therefore, both trilinear terms have been bounded by the deconvolution error and some terms that will eventually be hidden on the left side of Equation (52). The proof continues as in [22] by applying the discrete Gronwall’s inequality. The application of the Gronwall’s inequality conditions the time-step to satisfy $\Delta t\le C({\nu }^{-3}{\left|\right||\nabla \mathbf{u}|\left|\right|}_{\infty ,0}{+1)}^{-1}$.

For the modified Leray-deconvolution model, the proof follows the steps outlined above down to equality (53) which, in the case of the ML-$\alpha$, will take the form

$$\begin{array}{c}{b}^{*}({\mathbf{u}}^{n+1/2},{\mathbf{u}}^{n+1/2},D_N{\overline{{\phi }_h^{n+1/2}}}^h)-b^{*}({\mathbf{u}}^{n+1/2},D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\phi }_h^{n+1/2}}}^h)=\hfill \\ b^{*}({\mathbf{u}}^{n+1/2},{\mathbf{u}}^{n+1/2}-D_N{\overline{{\mathbf{u}}^{n+1/2}}}^h,D_N{\overline{{\phi }_h^{n+1/2}}}^h)\hfill \end{array}$$

(57)

which is further estimated as in inequality (54) from above. From here on, the proof proceeds as in the ADM case.

 ☐

4. A Numerical Experiment

In this section, we verify the theoretical convergence rates on a 2d NSE problem from [29] with exact solution

$$\begin{array}{c}\begin{array}{ccc}{u}_1(x,y,t)\hfill & =\hfill & sin\left(2\pi y\right)e^{-4\nu {\pi }^{2}t}\hfill \\ u_2(x,y,t)\hfill & =\hfill & sin\left(\pi x\right)e^{-\nu {\pi }^{2}t}\hfill \\ p(x,y,t)\hfill & =\hfill & 0\hfill \end{array}\end{array}$$

(58)

therefore, $\mathbf{f}=(f_1,f_2)$ is equal to

$$\begin{array}{c}\begin{array}{ccc}{f}_1(x,y,t)\hfill & =\hfill & 2\pi cos\left(2\pi y\right)sin\left(\pi x\right)e^{-5\nu {\pi }^{2}t}\hfill \\ f_2(x,y,t)\hfill & =\hfill & \pi cos\left(\pi x\right)sin\left(2\pi y\right)e^{-5\nu {\pi }^{2}t}\hfill \end{array}\end{array}$$

(59)

Herein, we check the Leray-deconvolution model with $N=0,1$ (thus suplementing the numerical results presented in [22]) and the modified Leray-deconvolution model with $N=0,1$. The ADM with $N=0,1$ has been checked in this problem in [29] on a slightly different mesh and is therefore omitted here. The NS-$\alpha$ deconvolution model has been checked numerically in [19].

The scope of the test is to verify the convergence rates and also to check numerically that, as N increases, the corresponding model enters the predicted convergent regime faster.

The test is carried out with the FreeFEM++ package [40]. The computational domain is the unit square $\mathsf{\Omega }=(0,1)\times (0,1)$. The kinematic viscosity is set equal to $\nu =1.0$ and the final time T is set equal to $T=0.1$. An initial mesh $M_0$ is generated with five even nodes on the edges $x=0,x=1$, 4 even nodes on $y=0$ and 6 even nodes on $y=1$. Then, the edges are succesivelly halved to generate the computational meshes $M_1,M_2,\cdots ,M_5$, see Figure 1. For the space discretization, we use $P2/P1$ Taylor–Hood finite elements. In our computation, we set the filter radius $\alpha$ equal to the mesh size h, i.e., $\alpha =h$.

First, we solve the Leray-deconvolution and modified Leray-deconvolution with $N=0,1$ on $M_1$ with $\Delta t=1/80$ and eight iterations. When solving the next mesh $M_2$, the time step $\Delta t$ is halved and we double the number of time iterations. We keep applying this procedure up to the last mesh $M_5$.

To stabilize pressure, the model and the filter equation are augmented with an $L^2$ pressure stabilization term with parameter ${10}^{-9}$. The resulting algebraic system corresponding to the discrete model is solved using a fixed point iteration until the $L^{\infty }$ norm of two successive iterates is less than ${10}^{-12}$. The rate presented in the

Table 1. $L^{\infty }\left(L^2\right)$ and $L^2\left(H^1\right)$ errors and rates for the the Leray-$\alpha$ model. The final time is $T=0.1$. The predicted rate order is 2.

Level	Nr Iter	$\mathbf{\Delta }\mathit{t}$	h	${\mathit{L}}^{\mathbf{\infty }}\left({\mathit{L}}^2\right)$ Error		${\mathit{L}}^2\left({\mathit{H}}^1\right)$ Error
$M_1$	8	1/80	0.198493	0.00349919	rate	0.00635556	rate
$M_2$	16	1/160	0.0992465	0.000966758	1.85	0.0017491	1.86
$M_3$	32	1/320	0.0496232	0.000267232	1.85	0.000472424	1.88
$M_4$	64	1/640	0.0248116	$6.78601\times {10}^{-5}$	1.97	0.000118784	1.99
$M_5$	128	1/1280	0.0124058	$1.68665\times {10}^{-5}$	2.00	$2.9308\times {10}^{-5}$	2.01

,

Table 2. $L^{\infty }\left(L^2\right)$ and $L^2\left(H^1\right)$ errors and rates for the Leray-deconvolution model with $N=1$. The final time is $T=0.1$. The predicted rate order is 2.

Level	Nr Iter	$\mathbf{\Delta }\mathit{t}$	h	${\mathit{L}}^{\mathbf{\infty }}\left({\mathit{L}}^{\mathbf{2}}\right)$ Error		${\mathit{L}}^{\mathbf{2}}\left({\mathit{H}}^{\mathbf{1}}\right)$ Error
$M_1$	8	1/80	0.198493	0.00321879	rate	0.00599466	rate
$M_2$	16	1/160	0.0992465	0.000735394	2.12	0.00131779	2.18
$M_3$	32	1/320	0.0496232	0.000174567	2.07	0.000291998	2.17
$M_4$	64	1/640	0.0248116	$4.29922\times {10}^{-5}$	2.02	$6.89069\times {10}^{-5}$	2.08
$M_5$	128	1/1280	0.0124058	$1.07004\times {10}^{-5}$	2.00	$1.67512\times {10}^{-5}$	2.04

,

Table 3. $L^{\infty }\left(L^2\right)$ and $L^2\left(H^1\right)$ errors and rates for the modified Leray-$\alpha$ model. The final time is $T=0.1$. The predicted rate order is 2.

Level	Nr Iter	$\mathbf{\Delta }\mathit{t}$	h	${\mathit{L}}^{\mathbf{\infty }}\left({\mathit{L}}^2\right)$ Error		${\mathit{L}}^2\left({\mathit{H}}^1\right)$ Error
$M_1$	8	1/80	0.198493	0.0059773	rate	0.00995501	rate
$M_2$	16	1/160	0.0992465	0.00259845	1.20	0.00434837	1.19
$M_3$	32	1/320	0.0496232	0.000900568	1.52	0.00149832	1.53
$M_4$	64	1/640	0.0248116	0.000257406	1.80	0.000427468	1.80
$M_5$	128	1/1280	0.0124058	$6.76934\times {10}^{-5}$	1.92	0.000112259	1.92

and

Table 4. $L^{\infty }\left(L^2\right)$ and $L^2\left(H^1\right)$ errors and rates for the modified Leray-deconvolution model with $N=1$. The final time is $T=0.1$. The predicted rate order is 2.

Level	Nr Iter	$\mathbf{\Delta }\mathit{t}$	h	${\mathit{L}}^{\mathbf{\infty }}\left({\mathit{L}}^2\right)$ Error		${\mathit{L}}^2\left({\mathit{H}}^1\right)$ Error
$M_1$	8	1/80	0.198493	0.00453629	rate	0.00788281	rate
$M_2$	16	1/160	0.0992465	0.00111931	2.01	0.00198618	1.98
$M_3$	32	1/320	0.0496232	0.000209968	2.41	0.000373953	2.40
$M_4$	64	1/640	0.0248116	$4.41691\times {10}^{-5}$	2.24	$7.65631\times {10}^{-5}$	2.28
$M_5$	128	1/1280	0.0124058	$1.06889\times {10}^{-5}$	2.04	$1.80013\times {10}^{-5}$	2.08

is computed as the $lo{g}_2$ of the quotient of two successive errors. The theoretical rate in the $L^{\infty }\left(L^2\right)$ and $L^2\left(H^1\right)$ norms for the Leray-deconvolution and modified Leray-deconvolution models is 2 and is confirmed by this numerical test. Moreover, the predicted convergent regime is reached faster for $N=1$.

5. Conclusions

This report presents some results on the theory of $\alpha$ models of turbulence that have been obtained in recent years. They show that the accuracy of these models can be improved using higher order approximate deconvolution operators. Some numerical tests are also presented to support the theoretical results. The technique to estimate the deconvolution error can be used in the larger context of models using approximate deconvolution such as the time-relaxation model investigated in [41] or for other models, such as the MHD, that have been recently investigated using approximate deconvolution methods.

References