A Modular Regularized Variational Multiscale Proper Orthogonal Decomposition for Incompressible Flows

TL;DR

This paper introduces a modular, post-processing VMS-POD method for incompressible flows, enabling easier implementation and independent stabilization tuning, supported by theoretical analysis and numerical validation.

Contribution

It presents a novel decoupled, post-processing VMS stabilization integrated with POD for Navier-Stokes equations, enhancing flexibility and effectiveness.

Findings

The method effectively stabilizes simulations of incompressible flows.

Numerical tests confirm the theoretical properties and improved accuracy.

The approach is easily incorporated into existing computational frameworks.

Abstract

In this paper, we propose, analyze and test a post-processing implementation of a projection-based variational multiscale (VMS) method with proper orthogonal decomposition (POD) for the incompressible Navier-Stokes equations. The projection-based VMS stabilization is added as a separate post-processing step to the standard POD approximation, and since the stabilization step is completely decoupled, the method can easily be incorporated into existing codes, and stabilization parameters can be tuned independent from the time evolution step. We present a theoretical analysis of the method, and give results for several numerical tests on benchmark problems which both illustrate the theory and show the proposed method's effectiveness.

Figures

Tables

Table 1. Table 1 : Convergence of the VMS-POD for varying R 𝑅 R .

$r$	$R$	$\epsilon$	$\Delta t$	${\Vert 𝐮-{𝐮}_r\Vert }_{L^{\infty }\left(L^2\right)}$	rate	${\Vert \nabla \left(𝐮-{𝐮}_r\right)\Vert }_{L^2\left(H^1\right)}$	rate
12	3	18.1675	0.002	0.0294	-	1.4060	-
12	5	10.4747	0.002	0.0158	1.13	1.0639	0.51
12	7	4.7609	0.002	0.0076	0.93	0.7039	0.52
12	9	2.6534	0.002	0.0034	1.37	0.6131	0.24
12	11	1.2920	0.002	0.0090	-1.35	0.7580	-0.29

Table 2. Table 2 : Convergence of the VMS-POD for varying Δ ​ t Δ 𝑡 \Delta t .

$r$	$R$	$\Delta t$	${\Vert 𝐮-{𝐮}_r\Vert }_{L^{\infty }\left(L^2\right)}$	rate	${\Vert \nabla \left(𝐮-{𝐮}_r\right)\Vert }_{L^2\left(H^1\right)}$	rate
16	14	0.032	0.4343	-	15.6394	-
16	14	0.016	0.3587	0.2759	10.8125	0.5325
16	14	0.008	0.1296	1.4687	4.5266	1.2562
16	14	0.004	0.0267	2.2792	1.0657	2.0866

Equations

$$\begin{array}[]{rcll}{\bf u}_{t}-\nu\Delta{\bf u}+({\bf u}\cdot\nabla){\bf u}+\nabla p&=&{\bf f}&\mathrm{in}\ (0,T]\times\Omega,\\ \nabla\cdot{\bf u}&=&0&\mathrm{in}\ [0,T]\times\Omega,\\ {\bf u}&=&\mathbf{0}&\mathrm{in}\ [0,T]\times\partial\Omega,\\ {\bf u}(0,{\bf x})&=&{\bf u}_{0}&\mathrm{in}\ \Omega,\\ \displaystyle\int_{\Omega}p\ d{\bf x}&=&0,&\mathrm{in}\ (0,T].\end{array}$$

$${\Vert \mathbf{f}\Vert }_{-1}=\underset{\mathbf{v}\in X}{\sup }\frac{|(\mathbf{f},\mathbf{v})|}{\Vert \nabla \mathbf{v}\Vert }.$$

$$|||{\bf v}|||_{\infty,p}:=\max_{0\leq n\leq M}\|{\bf v}^{n}\|_{p}\quad|||{\bf v}|||_{m,p}:=\bigg{(}\Delta t\sum_{n=0}^{M}\|{\bf v}^{n}\|_{p}^{m}\bigg{)}^{1/m}.$$

$$\begin{array}[]{rcll}({\bf u}_{t},{\bf v})+(\nu\nabla{{\bf u}},\nabla{{\bf v}})+b({\bf u},{\bf u},{\bf v})-(p,\nabla\cdot{\bf v})&=&({\bf f},{\bf v})&\forall\ {\bf v}\in{\bf X},\\ (q,\nabla\cdot{\bf u})&=&0&\forall\ q\in Q,\end{array}$$

$$b(\mathbf{u},\mathbf{v},\mathbf{w})=\frac{1}{2}\left(((\mathbf{u}\cdot \nabla )\mathbf{v},\mathbf{w})-((\mathbf{u}\cdot \nabla )\mathbf{w},\mathbf{v})\right)$$

$$b(\mathbf{u},\mathbf{v},\mathbf{w})$$

$$\underset{q_h\in Q^h}{\inf }\underset{{\mathbf{v}}_h\in {\mathbf{X}}^h}{\sup }\frac{(q_h,\nabla \cdot {\mathbf{v}}_h)}{||\nabla {\mathbf{v}}_h||||q_h||}\ge \beta >0.$$

$$\underset{{\mathbf{v}}_h\in X^h}{\inf }\left(\Vert (\mathbf{u}-{\mathbf{v}}_h)\Vert +h\Vert \nabla (\mathbf{u}-{\mathbf{v}}_h)\Vert \right)$$

$$\underset{q_h\in Q^h}{\inf }\Vert p-q_h\Vert$$

$${\mathbf{V}}^h=\{{\mathbf{v}}_h\in {\mathbf{X}}^h:(\nabla \cdot {\mathbf{v}}_h,q_h)=0,\forall q_h\in Q^h\},$$

$$\begin{array}[]{rcll}({\bf u}_{h,t},{\bf v}_{h})+(\nu\nabla{{\bf u}_{h}},\nabla{{\bf v}_{h}})+b({\bf u}_{h},{\bf u}_{h},{\bf v}_{h})=({\bf f},{\bf v}_{h}),\quad\forall\ {\bf v}_{h}\in{\bf V}^{h}.\end{array}$$

$$\mathcal{R}=span\{\mathbf{u}(\cdot ,t_1),\dots ,\mathbf{u}(\cdot ,t_M)\}$$

$$\min \frac{1}{M}\sum _{k=1}^M\Vert \mathbf{u}(\cdot ,t_k)-\sum _{i=1}^r(\mathbf{u}(\cdot ,t_k),{\mathbf{\psi }}_i(\cdot )){\mathbf{\psi }}_i(\cdot ){\Vert }^2,$$

$${\mathbf{\psi }}_l(\cdot )=\frac{1}{\sqrt{{\lambda }_l}}\sum _{i=1}^M({\mathbf{v}}_l{)}_i\mathbf{u}(\cdot ,t_i),1\le l\le r,$$

$$\frac{1}{M}\sum _{k=1}^M\Vert \mathbf{u}(\cdot ,t_k)-\sum _{i=1}^r(\mathbf{u}(\cdot ,t_k),{\mathbf{\psi }}_i(\cdot )){\mathbf{\psi }}_i(\cdot ){\Vert }^2=\sum _{i=r+1}^d{\lambda }_i.$$

$$\begin{array}[]{rcll}({\bf u}_{r,t},{\boldsymbol{\psi}})+(\nu\nabla{{\bf u}_{r}},\nabla{\boldsymbol{\psi}})+b({\bf u}_{r},{\bf u}_{r},{\boldsymbol{\psi}})=({\bf f},{\boldsymbol{\psi}}),\quad\forall\ {\boldsymbol{\psi}}\in{\bf X}^{r}.\end{array}$$

$${\mathbf{u}}_r(\mathbf{x},t):=\sum _{i=1}^r{a}_j(t){\mathbf{\psi }}_j(\mathbf{x}),$$

$$(M_r{)}_{i,j}$$

$$(S_r{)}_{i,j}$$

$$\Vert \nabla {\mathbf{u}}_r\Vert$$

$$\frac{1}{M}\sum _{k=1}^M\Vert {\mathbf{u}}_h(\cdot ,t_k)-\sum _{m=1}^r({\mathbf{u}}_h(\cdot ,t_k),{\mathbf{\psi }}_m(\cdot )){\mathbf{\psi }}_m(\cdot ){\Vert }_1^2=\sum _{m=r+1}^d\Vert {\mathbf{\psi }}_m{\Vert }_1^2{\lambda }_m.$$

$$(\mathbf{u}-P_r\mathbf{u},{\mathbf{\psi }}_r)=0,\forall {\mathbf{\psi }}_r\in {\mathbf{X}}^r.$$

$$\frac{1}{M}\sum _{n=1}^M\Vert {\mathbf{u}}^n-P_r{\mathbf{u}}^n{\Vert }^2$$

$$\frac{1}{M}\sum _{n=1}^M\Vert \nabla ({\mathbf{u}}^n-P_r{\mathbf{u}}^n){\Vert }^2$$

$$\Vert {\mathbf{u}}^n-P_r{\mathbf{u}}^n{\Vert }^2$$

$$\Vert \nabla ({\mathbf{u}}^n-P_r{\mathbf{u}}^n){\Vert }^2$$

$${\mathbf{X}}^R$$

$${\mathbf{L}}^R$$

$$(\mathbf{u}-P_R\mathbf{u},{\mathbf{v}}_R)=0,\forall {\mathbf{v}}_R\in {\mathbf{L}}^R.$$

$$\displaystyle\qquad\bigg{(}\frac{{\bf w}_{r}^{n+1}-{\bf u}_{r}^{n}}{\Delta t},{\boldsymbol{\psi}}\bigg{)}+b({\bf w}_{r}^{n+1},{\bf w}_{r}^{n+1},{\boldsymbol{\psi}})+\nu(\nabla{\bf w}_{r}^{n+1},\nabla{\boldsymbol{\psi}})=({\bf f}^{n+1},{\boldsymbol{\psi}}).$$

$$\displaystyle\qquad\bigg{(}\frac{{\bf w}_{r}^{n+1}-{\bf u}_{r}^{n+1}}{\Delta t},{\boldsymbol{\psi}}\bigg{)}=(\nu_{T}(I-P_{R})\nabla\frac{({\bf w}_{r}^{n+1}+{\bf u}_{r}^{n+1})}{2},(I-P_{R})\nabla{\boldsymbol{\psi}}),$$

$$\Vert {\mathbf{w}}_r^{n+1}{\Vert }^2=\Vert {\mathbf{u}}_r^{n+1}{\Vert }^2+2{\nu }_T\Delta t{\Vert (I-P_R)\nabla \frac{({\mathbf{w}}_r^{n+1}+{\mathbf{u}}_r^{n+1})}{2}\Vert }^2.$$

Full text

Keywords: proper orthogonal decomposition, projection-based variational multiscale, reduced order models, post-processing

1 Introduction

We consider the incompressible Navier-Stokes equations (NSE) on a polyhedral domain $\Omega\subset\mathbb{R}^{d},d\in\{2,3\}$ with boundary $\partial\Omega$:

[TABLE]

Here, ${\bf u}(t,{\bf x})$ is the fluid velocity and $p(t,{\bf x})$ the fluid pressure. The parameters in (1.1) are the kinematic viscosity $\nu>0$, the prescribed body forces ${\bf f}(t,{\bf x})$ and the initial velocity field ${\bf u}_{0}({\bf x})$.

It is known that due to the wide range of scales in many complex fluid flows, simulating these flows by a direct numerical simulation (DNS) can be very expensive, and sometimes is even infeasible. In particular, in the engineering design process, flow simulations must be run many times, e.g. to perform parameter studies or for system control purposes; this multiplies the DNS cost by at least several times. The concept of reduced order models was introduced as a way of lowering the computational complexity in such settings. The proper orthogonal decomposition (POD) approach has proven to be quite successful for generating reduced order models [28, 30, 29] that capture many (sometimes most) of the dominant flow features, and has been applied in many areas such as image processing, pattern recognition, unsteady aerodynamic applications etc., see e.g. [2, 5, 6]. POD reduces the complexity of systems, often by orders of magnitude, by representing it with only its most energetic structures.

However, it is well known that for certain problems, and incompressible NSE in particular, a POD by itself can perform quite poorly without some sort of numerical stabilization [1, 10]. One type of stabilization that has been successful in this endeavor is the combining of POD with the variational multiscale (VMS) method. VMS was introduced in [9] in a variational setting, and it (and its variants) have been studied extensively in finite element frameworks [8, 18, 12, 13, 14]. Using VMS in POD was pioneered in [10, 11, 26], and their studies showed this could increased numerical accuracy for convection-dominated convection-diffusion equations [11] and for NSE [10, 26]. Furthermore, in [11], an analysis was performed that showed optimal error bounds could be obtained (in terms of mesh width, time step size, and eigenvalues and eigenvectors removed from the system) for convection-diffusion-reaction VMS-POD systems.

The main objective of this report is to extend the novel ideas of [21] to the POD setting, in particular to create a VMS-POD where stabilization is added as a completely decoupled second step, in a time stepping scheme for incompressible flow simulation. That is, at each step there is a two step procedure at each time step: step 1 evolves with a standard POD (i.e. unstabilized Galerkin POD), and then the second step is a weighted POD projection that adds (in a sense) extra viscosity to the lower POD modes (hence the method can be easily incorporated into a standard Galerkin POD code). We formally introduce the method in section 3 (after giving necessary notation and preliminaries in section 2), and prove it is stable and optimally convergent (in terms of the mesh width, time step size, and removed POD modes/eigenvalues). In section 4 we provide extensive numerical results that illustrate the effectiveness of the method. Finally, conclusions are given in section 5.

2 Preliminaries and Notations

We assume that $\Omega\subset\mathbb{R}^{d},\,d=2,3$ is a polygonal or polyhedral domain, with boundary $\partial\Omega$. Throughout the paper standard notations for Sobolev spaces and their norms will be used. The norm in $(H^{k}(\Omega))^{d}$ is denoted by $\|\cdot\|_{k}$ and the norms in Lebesgue spaces $(L^{p}(\Omega))^{d}$, $1\leq p<\infty$, $p\neq 2$ by $\|\cdot\|_{L^{p}}$. The space $L^{2}(\Omega)$ is equipped with the norm and inner product $\|\cdot\|$ and $(\cdot,\cdot)$, respectively.

The continuous velocity and pressure spaces are denoted by ${\bf X}=(H^{1}_{0}(\Omega))^{d}$ and $Q=L_{0}^{2}(\Omega)$, and we denote the dual space of $H_{0}^{1}(\Omega)$ by $H^{-1}$ with norm

[TABLE]

We use the following notation for discrete norms, for ${\bf v}^{n}\in{\bf H}^{p}(\Omega)$, n=0,1,2,…,M:

[TABLE]

The variational formulation of (1.1) reads as follows: Find ${\bf u}\ :\ (0,T]\to{\bf X},\ p\ :\ (0,T]\to Q$ satisfying

[TABLE]

with ${\bf u}(0,{\bf x})={\bf u}_{0}({\bf x})\in{\bf X}$ and

[TABLE]

is the skew-symmetric form of the convective term.

The following properties for the skew symmetric form are well known [7, 19].

Lemma 2.1**.**

The trilinear skew-symmetric form $b({\bf u},{\bf v},{\bf w})$ satisfies

[TABLE]

This paper considers a conforming finite element method for (2.1), with spaces ${\bf X}^{h}\subset{\bf X}$ and $Q^{h}\subset Q$ satisfying the inf-sup condition: there is a constant $\beta$ independent of the mesh size $h$ such that

[TABLE]

It will also be assumed for simplicity that the finite element spaces ${\bf X}^{h}$ and $Q^{h}$ are composed of piecewise polynomials of degrees at most $m$ and $m-1$, respectively (the analysis can be extended without significant difficulty to any inf-sup stable pair). In addition, we assume that the spaces satisfy the following approximation properties

[TABLE]

for ${\bf u}\in{\bf X}\cap{\bf H}^{m+1}(\Omega)$ and $p\in Q\cap H^{m}(\Omega)$.

We denote the discretely divergence free space by

[TABLE]

The inf-sup condition (2.2) implies that the space ${\bf V}^{h}$ is closed subspace of ${\bf X}^{h}$ and the formulation in ${\bf X}^{h}$ is equivalent to ${\bf V}^{h}$. Thus, the Galerkin finite element approximation of (2.1) in ${\bf V}^{h}$ has the following form: Find ${\bf u}_{h}\in{\bf V}^{h}$ satisfying

[TABLE]

@fb@secFB

2.1 POD preliminaries

In this section the essentials of the POD are described. For a more detailed description of the method, see e.g. [23]. With the same notation of [10], consider the finite number of the time instances,

[TABLE]

at time $t_{i}=i\Delta t$, $i=1,\dots,M$ and let $\Delta t=\frac{T}{M}$, where $rank(\mathcal{R})=d$. In what follows, these time instances will be assumed to come from a DNS computed with a finite element spatial discretization. That is, we replace ${\bf u}(\cdot,t_{k})$ with DNS computed solution ${\bf u}_{h}(\cdot,t_{k})$.

Let $\{{\boldsymbol{\psi}}_{1},{\boldsymbol{\psi}}_{2},\dots,{\boldsymbol{\psi}}_{r}\}$ be the low-dimensional ordered POD basis functions to approximate these time instances. The low-dimensional space is obtained by solving the minimization problem of

[TABLE]

such that $({\boldsymbol{\psi}}_{i},{\boldsymbol{\psi}}_{j})=\delta_{ij}$, $1\leq i,j\leq r$ and $r<<d$. The solution of the problem (2.7) yields

[TABLE]

where $({\bf v}_{l})_{i}$ is the $i^{th}$ component of the eigenvector ${\bf v}_{l}$ corresponding to $\lambda_{l}$ which is the eigenvalue of the snapshots correlation matrix. We note that all eigenvalues are sorted in descending orders. Thus, the basis functions $\{{\boldsymbol{\psi}}_{1},{\boldsymbol{\psi}}_{2},\dots,{\boldsymbol{\psi}}_{r}\}$ correspond to the first $r$ largest eigenvalues. Moreover, the error estimation satisfies (see [16]):

[TABLE]

Let ${\bf X}^{r}=span\{{\boldsymbol{\psi}}_{1},{\boldsymbol{\psi}}_{2},\dots,{\boldsymbol{\psi}}_{r}\}$ be the POD-space, then POD-Galerkin (POD-G) formulation of the NSE is : Find ${\bf u}_{r}\in{\bf X}^{r}$ satisfying

[TABLE]

Note that the POD-G solution of the NSE is constructed by writing,

[TABLE]

where $a_{j}(t)$ are time varying coefficients.

To carry out the error analysis, we state the following error estimations without their proofs. Let $M_{r}$ and $S_{r}$ denote the POD mass matrix and stiffness matrix, respectively, with

[TABLE]

Note that, since we use $L^{2}(\Omega)$ to generate snapshots, $M_{r}=I_{r\times r}$.

Our analysis will utilize the following POD inequalities, proven in [16]: for all ${\bf u}_{r}\in{\bf X}^{r}$,

[TABLE]

where $\|\cdot\|_{2}$ denotes the matrix $2$-norm.

The results of the following lemmas will be applied to bound the POD projection error. The proofs can be found in [10].

Lemma 2.2**.**

The error in POD projection for the snapshots ${\bf u}_{h}(\cdot,t_{k})$, $k=1,\dots,M$ satisfies

[TABLE]

The proof of the finite element error estimate consists of the splitting the error into an approximation term and a finite element remainder term. To decompose the error term we use the $L^{2}$ projection of ${\bf u}$, which fulfills certain interpolation estimates. Lemma 2.3 estimates the error between the snapshots and their $L^{2}$ projection into ${\bf X}^{r}$.

Let $P_{r}$ denote a projection operator $P_{r}:{\bf L}^{2}\rightarrow{\bf X}^{r}$ that satisfies

[TABLE]

Lemma 2.3**.**

For ${\bf u}^{n}$ the true NSE solution at time $t^{n}$, the difference ${\bf u}^{n}-P_{r}{\bf u}^{n}$ satisfies

[TABLE]

Proof.

This follows from the projection error estimate from [10] along with the triangle inequality, after adding and subtracting the FEM-DNS snapshots to the true solution (at corresponding times). ∎

Assumption: We make the assumption that

[TABLE]

This assumption essentially says that the projection error the true solution is of the same order of magnitude at most of the time steps. It is possible the assumption may not hold (in which case the $C$ on the right hand side with become $C\Delta t^{-1}$) in diabolical cases where the true NSE solution lives in the ROM space at most of the times $t^{n}$. The assumption is common in error analysis for POD type methods [10, 31].

In the projection-based VMS method, besides the standard finite element spaces representing all resolved scales an additional large resolved scale is needed. For VMS-POD setting, the following spaces are used for $R<r$:

[TABLE]

Note that from the construction, we have ${\bf X}^{R}\subset{\bf X}^{r}\subset{\bf X}^{h}\subset{\bf X}$.

The $L^{2}$ orthogonal projection $P_{R}:{\bf L}^{2}\rightarrow{\bf L}^{R}$ will be needed in the VMS formulation, and is defined by

[TABLE]

3 Post-Processed VMS-POD Schemes

This section proposes fully discrete VMS-POD methods for solving (2.10). For simplicity, in section 3.1, we analyze the backward Euler temporal discretization. In section 3.2, extension to BDF2 time stepping is considered. In our analysis, we assume that the eddy viscosity coefficient $\nu_{T}$ is known bounded, positive and element-wise constant. The results of this report can be extended in the case $\nu_{T}$ is nonconstant even nonlinear. The consideration of a nonlinear $\nu_{T}$ requires more complex mathematical theory due to the strong monocity, see [17].

In that analysis that follows, we denote variables at time $t^{n}=n\Delta t,n=0,1,2,\dots,M,T:=M\Delta t$ using superscripts, e.g. ${\bf f}^{n}:={\bf f}(t^{n})$.

@fb@secFB

3.1 Backward Euler

The two step VMS-POD scheme equipped with backward Euler time stepping reads as follows:

Algorithm 3.1**.**

Let ${\bf f}\in L^{2}(0,T;{\bf H}^{-1}(\Omega))$ and ${\bf u}_{r}^{0}={\bf w}_{r}^{0}$ be given with $L^{2}$ projection of ${\bf u}_{0}$ in ${\bf X}^{r}$. Given ${\bf u}_{r}^{n}\in{\bf X}_{r}$ compute ${\bf u}_{r}^{n+1}$ by applying the following two steps:

Step 1.* Calculate ${\bf w}_{r}^{n+1}\in{\bf X}^{r}$ satisfying $\forall{{\boldsymbol{\psi}}\in{\bf X}^{r}}$,*

[TABLE]

Step 2.* Post-process ${\bf w}_{r}^{n+1}$ by applying projection $P_{R}$ to obtain ${\bf u}_{r}^{n+1}\in{\bf X}_{r}$, $\forall{{\boldsymbol{\psi}}\in{\bf X}^{r}}$:*

[TABLE]

We note that Step 1 is the standard Galerkin POD method, and step 2 is completely decoupled VMS stabilization step. The projection in Step 2 is not a filter but constructed to recover VMS eddy viscosity term as in [21].

Note that if we let ${\boldsymbol{\psi}}=\frac{({\bf w}_{r}^{n+1}+{\bf u}_{r}^{n+1})}{2}$ in (3.2), the numerical dissipation induced from Step 2 is immediately seen to be

[TABLE]

3.1.1 Stability of Algorithm 3.1

We now prove stability of Algorithm 3.1.

Lemma 3.1**.**

The post-processed-VMS-POD approximation is unconditionally stable in the following sense: for any $\Delta t>0$,

[TABLE]

Proof.

Letting ${\boldsymbol{\psi}}={\bf w}^{n+1}_{r}$ in (3.1) and using the polarization identity yields

[TABLE]

Substitute (3.3) in (3.4) and multiply both sides by $2\Delta t$, which provides

[TABLE]

Bounding the forcing term in the usual way, and then summing over the time steps gives the stated result. ∎

The result of Lemma 3.1 also establishes the stability of ${\bf w}_{r}^{M}$.

Corollary 3.1**.**

(Stability of ${\bf w}_{r}^{M}$)

[TABLE]

Proof.

Expand the summation in Lemma 3.1 for $n=M-1$ and use (3.3). ∎

3.1.2 A Priori Error Estimation

In this section, we present the error analysis of the true solution of Navier Stokes equations and POD approximation (3.1)-(3.2). The optimal asymptotic error estimation requires the following regularity assumptions for the true solution:

[TABLE]

Theorem 3.1**.**

Suppose (3.7) holds and ${\bf u}_{r}^{n}$ and ${\bf w}_{r}^{n}$ given by Algorithm (3.1)-(3.2). For sufficiently small $\Delta t$, i.e. $\Delta t\leq[C\nu^{-3}\|\nabla{\bf u}\|^{4}_{\infty,0}]^{-1}$ we have the following estimation:

[TABLE]

where $C$ is independent from $\Delta t,h,\nu$ and $\nu_{T}$.

Remark 3.1**.**

Under the assumptions of Theorem 3.1 and the finite element spaces $({\bf X}^{h},Q^{h})$ with piecewise polynomials of degree $m$ and $m-1$, respectively. We obtain the following asymptotic error estimation:

[TABLE]

Proof.

We begin the proof by deriving error equations. From (2.1), we have that true solution $({\bf u},p)$ at time level $t=t^{n+1}$ satisfies, for ${\boldsymbol{\psi}}_{r}\in X_{r}$,

[TABLE]

where

[TABLE]

We define the following notations:

[TABLE]

where ${\bf\mathcal{U}^{n}}$ is $L^{2}$ projection of ${\bf u}^{n}$ in ${\bf X}_{r}$.

Subtracting (3.1) from (3.8) yields

[TABLE]

Substitute ${\boldsymbol{\varepsilon}}_{r}^{n}={\boldsymbol{\eta}}^{n}-{\boldsymbol{\phi}}_{r}^{n}$ and ${\bf e}_{r}^{n}={\boldsymbol{\eta}}^{n}-{\boldsymbol{\theta}}_{r}^{n}$ with ${\boldsymbol{\psi}}={\boldsymbol{\phi}}_{r}^{n+1}$ in the last equation we obtain

[TABLE]

From the definition of $L^{2}$ projection (2.13), we note that $({\boldsymbol{\eta}}^{n},{\boldsymbol{\phi}}_{r}^{n+1})=0$ and $({\boldsymbol{\eta}}^{n+1},{\boldsymbol{\phi}}_{r}^{n+1})=0$. Using this along with the polarization identity and that ${\boldsymbol{\phi}}_{r}^{n+1}\in{\bf X}^{r}\subset{\bf V}^{h}$, we obtain the bound

[TABLE]

The first term on the right hand side of (3.10), and the pressure term, can be bounded using Cauchy-Schwarz and Young’s inequalities,

[TABLE]

For the nonlinear terms, first add and subtract terms to get

[TABLE]

Now using Lemma 2.1, Young’s and Poincaré’s inequalities, the terms in (3.13) are estimated as follows:

[TABLE]

The consistency error in (3.10) is estimated by

[TABLE]

Collecting all bounds for the right hand side terms of (3.10) and multiplying both sides by $2\Delta t$ gives

[TABLE]

We next get a bound for $\|{\boldsymbol{\phi}}_{r}^{n+1}\|^{2}$. Write (3.2) by adding and subtracting the true solution projection ${\bf\mathcal{U}}^{n+1}$ on both sides to get

[TABLE]

and then choose ${\boldsymbol{\psi}}=\frac{({\boldsymbol{\phi}}_{r}^{n+1}+{\boldsymbol{\theta}}_{r}^{n+1})}{2}$ in (3.16) to obtain

[TABLE]

Noting ${\bf\mathcal{U}}^{n+1}={\bf u}^{n+1}-{\boldsymbol{\eta}}^{n+1}$ and inserting (3.17) into (3.15) results into

[TABLE]

Assume now \Delta t\leq\dfrac{1}{8C}\Big{[}\frac{\|\nabla{\bf u}\|^{4}}{\nu^{3}}\Big{]}^{-1}. Then, we can bound the remaining right hand terms of (3.18) as follows :

[TABLE]

Similarly, using the preceding bound and Young’s inequality we get

[TABLE]

and

[TABLE]

Substitute the bounds (3.1.2)-(3.20) into (3.18) and sum from $n=0$ to $M-1$. This gives

[TABLE]

Note that since ${\bf u}_{r}^{0}={\bf\mathcal{U}}^{0}$, the first right hand side term in (3.21) vanishes. The second right hand side term is majorized using Lemma 2.3, as

[TABLE]

Third term is bounded by using (2.17) as follows,

[TABLE]

To bound the fourth term, one can proceed as follows. Using $\|(I-P_{R})\nabla{\boldsymbol{\eta}}^{n+1}\|\leq\|\nabla{\boldsymbol{\eta}}^{n+1}\|$ along with Lemma 2.3 leads to

[TABLE]

For the first term on the right-hand side of (3.24), we use (2.19) to find that

[TABLE]

where ${\bf\mathcal{U}}_{R}$ is the large scale representation of the projection. Now using Lemma 2.3, the final estimation for (3.24) becomes

[TABLE]

Corollary 3.1, the stability of ${\bf w}_{r}$ and (2.17) provide an estimation for the fifth term in the right-hand side of (3.21):

[TABLE]

The estimation of sixth term in the right hand side of (3.21) uses approximation property (2.4) to find

[TABLE]

Finally, for the last term in the right hand side of (3.21), Taylor series expansion with remainder in integral form is used along with Cauchy Schwarz and the triangle inequality to obtain

[TABLE]

Collecting all the bounds (3.22)-(3.29) for (3.21) yields

[TABLE]

Again using the assumption that \Delta t\leq\frac{1}{8C}\bigg{[}\nu^{-3}{|||\nabla{\bf u}|||_{\infty,0}^{4}}\bigg{]}^{-1} allows us to apply the discrete Gronwall inequality, which yields

[TABLE]

Finally, the triangle inequality is applied to produce the stated result. ∎

@fb@secFB

3.2 Extension to second order time stepping

We consider now an extension of Algorithm 3.1 to BDF2 time stepping.

Algorithm 3.2**.**

*Let ${\bf f}\in L^{2}(0,T;{\bf H}^{-1}(\Omega))$ and initial conditions ${\bf u}_{r}^{0}$ and ${\bf u}_{r}^{-1}$ be given in ${\bf X}^{r}$. Then for n=0,1,2,…

Step 1. Calculate ${\bf w}_{r}^{n+1}\in{\bf X}^{r}$ satisfying $\forall{{\boldsymbol{\psi}}\in{\bf X}^{r}}$,*

[TABLE]

Step 2.* Post-process ${\bf w}_{r}^{n+1}$ to obtain ${\bf u}_{r}^{n+1}\in{\bf X}^{r}$ satisfying $\forall{{\boldsymbol{\psi}}\in{\bf X}^{r}}$,*

[TABLE]

We note the post-processing step is exactly the same as in the backward Euler case. Also as in the case of the backward Euler method above, without Step 2, Algorithm 3.2 reduces to the classical Galerkin POD formulation for the NSE, although now using BDF2 time stepping.

We now prove stability of Algorithm 3.2. A convergence result can be obtained by combining the ideas of Algorithm 3.1’s convergence proof with the stability proof below. Such a proof is thus long and technical, but produces the expected result (i.e. same convergence as Algorithm 3.1, but second order in $\Delta t$ instead of first order). This expected second order temporal convergence is illustrated in the numerical experiments section below.

Lemma 3.2**.**

The post-processed VMS-POD approximation is stable for the eddy viscosity term $\nu_{T}<4\nu$ in the following sense:

[TABLE]

Proof.

Note that if we let ${\boldsymbol{\psi}}=\frac{({\bf w}_{r}^{n+1}+{\bf u}_{r}^{n+1})}{2}$ in Step 2, the numerical dissipation induced from Step 2 is given by

[TABLE]

Letting ${\boldsymbol{\psi}}={\bf w}^{n+1}_{r}$ in (3.30) and using the identity

[TABLE]

yields

[TABLE]

Substitute (3.32) in (3.33), multiply both sides by $4\Delta t$, add $\|2{\bf u}_{r}^{n+1}-{\bf u}_{r}^{n}\|^{2}$ to both sides, and then apply the Cauchy-Schwarz inequality to get

[TABLE]

We now consider the term $\|2{\bf w}_{r}^{n+1}-{\bf u}_{r}^{n}\|^{2}-\|2{\bf u}_{r}^{n+1}-{\bf u}_{r}^{n}\|^{2}$ in (3.34) . By using the properties of $L^{2}$ inner product, the equality (3.32) and rearranging terms gives

[TABLE]

Inserting (3.35) in (3.34) and applying Cauchy-Schwarz and Young’s inequalities gives

[TABLE]

Adding and subtracting terms for the last term in the previous inequality with the use of $\|I-P_{R}\|\leq 1$, we get

[TABLE]

Summing over the time step $i=1,\cdots,M$ gives

[TABLE]

With the assumption ${\bf w}_{r}^{0}={\bf w}_{r}^{1}=0$ and $\|I-P_{R}\|<1$, we obtain the stated result. ∎

4 Numerical Studies

This section gives results for three numerical experiments. In all cases we use Algorithm 3.2, i.e. the scheme with second order time stepping. Our first test considers the predicted convergence rates of the previous section, with respect to varying $R$ and $\Delta t$. For the second test, we compare accuracy of the proposed VMS-POD scheme compared with the usual Galerkin POD method (i.e. unstabilized POD, computed by eliminating the post-processing step of the VMS-POD) in 2D channel flow past a cylinder. Finally we consider the VMS-POD for a 3D turbulent channel flow simulation.

Our first two numerical tests use the test problem of 2D channel flow past a cylinder. Here, the domain is a $2.2\times 0.41$ rectangle with a circle radius $=0.05$ centered at $(0.2,0.2)$ (see figure 1). The test problem uses no slip boundary conditions for the walls and cylinder, and the time dependent inflow and outflow profiles are given by

[TABLE]

The kinematic viscosity $\nu=10^{-3}$, and there is no forcing (${\bf f}={\bf 0}$). The POD is created as described in section 2, by taking snapshots from a BDF2-finite element DNS simulation using Taylor-Hood elements, after a periodic-in-time solution is reached (a more detailed description of the setup is given in [4]). The tests with the cylinder problem use the projection of the T=7 DNS solution as the initial condition, and run for $10$ time units.

@fb@secFB

4.1 Numerical test 1: Convergence in $R$ and $\Delta t$

We now test the predicted convergence rates of the previous section. This is a particularly difficult task for this problem, as the parameters are the spatial mesh width $h$, time step $\Delta t$, POD cutoff $r$, and VMS cutoff $R$. We assume the spatial mesh width is sufficiently small so that this error source is negligible. Our particular interest here is the scaling of the error with the VMS cutoff $R$ and with the time step size $\Delta t$. Due to having several parameters, to see convergence with respect to a particular parameter, it must be part of the dominant error source. Hence in our tests, different parameter choices are made to see the various scalings.

To test the scaling of the error with $R$, we fix $r=12$, $\Delta t=0.002$, $\nu_{T}=0.0003$, and compute with varying $R$. The error estimates are depend on $R$ via the term $\varepsilon=\sqrt{\sum\limits_{j=R+1}^{d}\|\varphi_{j}\|_{1}^{2}\lambda_{j}}$, as well as $\|S_{R}\|_{2}$, and our interest is the scaling with $\varepsilon$. We note that for sufficiently large $R$, the term $\|S_{R}\|_{2}$ can become sufficiently large so that it becomes a dominant error source. Results for this test are shown in table 1, and we observe the expected convergence rate of approximately $1/2$ in the $L^{2}(0,T;H^{1}(\Omega))$ norm, and seemingly a higher rate in $L^{\infty}(0,T;L^{2}(\Omega))$ . However, for larger $R$, we see a deterioration of the rate, and by $R=11$, the error increases (which is expected, due to increase in $\|S_{R}\|_{2}$ with $R$).

To test the scaling with respect to $\Delta t$, we fix $r=16$, $R=14$, and vary $\Delta t$. Errors and rates are given in table 2, and we observe rates consistent with second order.

@fb@secFB

4.2 Numerical test 2: Error comparison of VMS-POD versus POD-G for 2D channel flow past a cylinder

We now consider error comparison of the VMS-POD again the standard POD-G method. The statistics of interest are the maximal drag and the maximal lift coefficients at the cylinder. The reference intervals are given for $c_{dmax}$ and $c_{lmax}$ in [27],

[TABLE]

and recent computations of [4], as well as our DNS that created the snapshots, are in agreement with these numbers.

We compute with $r=8$ modes, and find that the POD-G method is not good, and in particular we observe in figure 2 that the energy is growing (seemingly) linearly with time. But $T=10$, a significant increase in energy has occurred, leading to very inaccurate lift and drag predictions as well.

The VMS-POD does a much better job with prediction of energy, lift and drag, however. Shown in figure 2 is the energy, lift and drag with $R=5$ (and $\nu_{T}=0.0003$ was found to be optimal to within $\pm$ $0.000005$, in terms of matching the DNS energy at $T=10$). We note that with POD, such a parameter optimization is quite cheap. Hence this is an excellent example of how the post-processing step of VMS-POD can remove the nonphysical energy growth that occurs with POD-G, and provide good reduced order solutions.

@fb@secFB

4.3 Numerical test 3: VMS-POD for $Re_{\tau}$=395 turbulent channel flow

For our last test, we consider the proposed VMS-POD for a benchmark turbulent channel flow problem at $Re_{\tau}=395$ [22, 15]. The domain $\Omega=(-2\pi,2\pi)\times(0,2)\times\left(-2\pi/3,2\pi/3\right)$, and for boundary conditions, no slip is enforced on the top and bottom walls $y=2$ and $y=0$, and periodic boundary conditions are enforced the remaining sides. The kinematic viscosity is taken to be $\nu=\frac{1}{Re_{\tau}}$, and the flow is forced in the $x$-direction via ${\bf f}=\langle 1,0,0\rangle$, and is also dynamically adjusted at each time step as in [15] to maintain the bulk velocity. We compute snapshots using the rNS-$\alpha$ model and scheme from [24, 25], by saving snapshots at each time step ($\Delta t=0.002$) of the simulation from $T=60$ to $T=70$, and then create the POD basis as described in section 2. Plots of some of the POD basis functions are displayed in figure 3. We construct the reduced order model using $r=10$ modes, and will test the ability of the VMS-POD to give stable solutions. We note that even though the snapshots were created with a turbulence model, we use the VMS-POD in the same form as above (i.e. step 1 is NSE, not the turbulence model).

Turbulent channel flow is well known to be a challenging numerical problem which causes many codes to fail and/or give poor solutions. Underresolved simulations suffer from a nonphysical accumulation of energy near the smallest resolved scales. In time, this accumulation of energy creates spurious oscillations, inaccuracies, and instability, often causing simulations to fail [20, 3]. Thus, we test here the ability of the VMS-POD to produce a stable solution. We use the projections of the $T=60$ and $T=60.002$ turbulence model solutions as the initial conditions, and run the reduced order models for 5 time units.

We first test the usual POD-G, without any stabilization (i.e. $R=r=10$). Energy versus time is shown for this simulation in figure 4, and we observe that there is an energy blowup. The instability occurs around $t=1$, and then the energy grows exponentially.

We now test the VMS-POD with the $r=10$, and using $R=3$, $2$, and $1$ ($R=3$ is the largest $R$ for which we could get stable energy up to $t=5$). Plots for energy versus time and mean streamwise velocity (average taken in both periodic directions and in time) are shown in figures 5-7 for these choices of $R$, and for each $R$ we test with several values of $\nu_{T}$. The main observation we make from these results is that with sufficient stabilization using $R$ and $\nu_{T}$, the VMS-POD is able to produce energy stable solutions whose mean streamwise velocity profile matches the DNS values of Moser et al [22] (we note that the first mode ${\boldsymbol{\psi}}_{1}$ also matches this mean streamwise profile as well). For smaller values of $R$, smaller values of $\nu_{T}$ can be sufficient to obtain energy stable solutions.

5 Conclusion

We proposed, analyzed and tested a VMS-POD method for incompressible NSE simulation, where the stabilization is completely decoupled into the second step of a two step implementation at each time step. Decoupling of the stabilization has the advantage of easily being incorporated into existing POD-G codes, and also that stabilization parameters can be adjusted only in the stabilization step (and not as part of the evolution equation). We rigorously prove an error estimate for the model, in terms of the number of POD modes $r$, the stabilization parameters $R$ (number of modes not to add stabilization to) and $\nu_{T}$, as well as the time step size $\Delta t$ and the mesh width $h$ of the underlying FEM simulation that produced the POD modes.

Results from several numerical experiments are provided that show how effective the method can be. In particular, we show for 2D channel flow past a step, POD-G has an energy growth that causes poor lift and drag prediction, especially for longer times. The proposed VMS-POD is able to fix this by stabilizing so that the energy matches the DNS energy, which in turn leads to excellent lift and drag prediction, even up to $t=10$ (and from the plots, it appear the accurate predictions can continue for even longer time). Finally, we tested the ability of the VMS-POD to produce stable solutions for turbulent channel flow. With $r=10$ modes, POD-G quickly becomes unstable, however the VMS-POD is able to produce simulations with stable energy and that yield accurate mean streamwise velocity profiles.