An Artificial Compression Reduced Order Model

TL;DR

The paper introduces a novel artificial compression reduced order model (AC-ROM) for simulating viscous incompressible flows, providing accurate velocity and pressure approximations without strict mathematical constraints on the basis functions.

Contribution

It presents a new AC-ROM that relaxes traditional stability conditions and constructs basis functions from non-divergence-free data, with proven error estimates and numerical validation.

Findings

AC-ROM accurately approximates velocity and pressure.

The model does not require inf-sup condition for basis functions.

Numerical experiments confirm the effectiveness of AC-ROM.

Abstract

We propose a novel artificial compression, reduced order model (AC-ROM) for the numerical simulation of viscous incompressible fluid flows. The new AC-ROM provides approximations not only for velocity, but also for pressure, which is needed to calculate forces on bodies in the flow and to connect the simulation parameters with pressure data. The new AC-ROM does not require that the velocity-pressure ROM spaces satisfy the inf-sup (Ladyzhenskaya-Babuska-Brezzi) condition and its basis functions are constructed from data that are not required to be weakly-divergence free. We prove error estimates for the reduced basis discretization of the AC-ROM. We also investigate numerically the new AC-ROM in the simulation of a two-dimensional flow between offset cylinders.

Equations

$$u_t+u\cdot \nabla u+\nabla p-\nu \Delta u=f,\text{and }\nabla \cdot u=0,\text{in}\Omega \times (0,T]$$

$$u=0,\text{on}\partial \Omega \times (0,T],\text{and }u(x,0)=u_0(x),\text{in}\Omega .$$

$$\Big{(}\frac{u^{n+1}_{R}-u^{n}_{R}}{\Delta t},\varphi\Big{)}+b^{\ast}(u_{R}^{n},u^{n+1}_{R},\varphi)+\nu(\nabla u^{n+1}_{R},\nabla\varphi)=(f^{n+1},\varphi),\quad\forall\varphi\in X_{R},$$

$$\underset{q_M\in Q_M}{\inf }\underset{v_R\in X_R}{\sup }\frac{(\nabla \cdot v_R,q_M)}{\Vert \nabla v_R\Vert \Vert q_M\Vert }\ge {\beta }_{is}>0.$$

$$\Delta p_M=-\nabla \cdot ((u_R\cdot \nabla )u_R)\text{mbox}in\Omega ,$$

$$\displaystyle\Big{(}\frac{u^{n+1}_{R}-u^{n}_{R}}{\Delta t},\varphi\Big{)}+b^{\ast}(u_{R}^{n},u^{n+1}_{R},\varphi)+\nu(\nabla u^{n+1}_{R},\nabla\varphi)$$

$$-(p_M^{n+1},\nabla \cdot \phi )=(f^{n+1},\phi )\forall \phi \in X_R$$

$$\epsilon \left(\frac{p_M^{n+1}-p_M^n}{\Delta t},\psi \right)+(\nabla \cdot u_R^{n+1},\psi )=0\forall \psi \in Q_M,$$

$$|||v|||_{2,s}:=\Big{(}\int_{0}^{T}\|v(\cdot,t)\|_{s}^{2}dt\Big{)}^{\frac{1}{2}}\qquad\text{and}\qquad|||v|||_{\infty,s}:=\text{ess\,sup}_{[0,T]}\|v(\cdot,t)\|_{s}.$$

$$\Vert f{\Vert }_{-1}=\underset{0\ne v\in X}{\sup }\frac{(f,v)}{\Vert \nabla v\Vert }\forall f\in H^{-1}(\Omega ).$$

$$X:=$$

$$Q:=$$

$$\{\begin{array}{rlr}(u_t,v)+(u\cdot \nabla u,v)+\nu (\nabla u,\nabla v)-(p,\nabla \cdot v)& =(f,v)& \forall v\in X\\ (\nabla \cdot u,q)& =0& \forall q\in Q\\ u(x,0)& =u^0(x).& \end{array}$$

$$\underset{v_h\in X_h}{\inf }\Vert v-v_h\Vert$$

$$\underset{v_h\in X_h}{\inf }\Vert \nabla (v-v_h)\Vert$$

$$\underset{q_h\in Q_h}{\inf }\Vert q-q_h\Vert$$

$$b(w,u,v)=(w\cdot \nabla u,v)\forall u,v,w\in [H^1(\Omega ){]}^d$$

$$b^{\ast }(w,u,v):=\frac{1}{2}(w\cdot \nabla u,v)-\frac{1}{2}(w\cdot \nabla v,u)\forall u,v,w\in [H^1(\Omega ){]}^d,$$

$$b^{\ast }(w,u,v)\le C_{b^{\ast }}(\Vert w\Vert \Vert \nabla w\Vert {)}^{1/2}\Vert \nabla u\Vert \Vert \nabla v\Vert \forall u,v,w\in X,$$

$$b^{\ast }(w,u,v)\le C_{b^{\ast }}\Vert \nabla w\Vert (\Vert u\Vert \Vert \nabla u\Vert {)}^{1/2}\Vert \nabla v\Vert \forall u,v,w\in X$$

$$\underset{n}{\sup }\Vert u^n-u_h^n{\Vert }^2+h^2\Vert \nabla (u^n-u_h^n){\Vert }^2$$

$$\underset{n}{\sup }\Vert p^n-p_h^n{\Vert }^2$$

$$\displaystyle\Big{(}\frac{u^{n+1}_{h}-u^{n}_{h}}{\Delta t},v_{h}\Big{)}+b^{\ast}(u_{h}^{n},u^{n+1}_{h},v_{h})+\nu(\nabla u^{n+1}_{h},\nabla v_{h})$$

$$-(p_h^{n+1},\nabla \cdot v_h)=(f^{n+1},v_h)\forall v_h\in X_h$$

$$\epsilon \left(\frac{p_h^{n+1}-p_h^n}{\Delta t},q_h\right)+(\nabla \cdot u_h^{n+1},q_h)=0\forall q_h\in Q_h.$$

$$\displaystyle\mathbb{A}=\big{(}{u}_{S}^{1},{u}_{S}^{2},\ldots,{u}_{S}^{N_{V}})\ \ \text{and}\ \ \mathbb{B}=\big{(}{p}_{S}^{1},{p}_{S}^{2},\ldots,{p}_{S}^{N_{P}}),$$

$$X_R:=\text{span}\{{\phi }_i{\}}_{i=1}^R\subset X_h\text{and}{Q}_M:=\text{span}\{{\psi }_i{\}}_{i=1}^M\subset Q_h,$$

$$\displaystyle\min\sum_{n=0}^{N}\Big{\|}u_{h,s}^{n}-\sum_{j=1}^{R}(u_{h,s}^{n},\varphi_{j})\varphi_{j}\Big{\|}^{2}$$

$$\displaystyle\text{subject to }(\varphi_{i},\varphi_{j})=\delta_{ij}\quad\mbox{for $i,j=1,\ldots,R$},$$

$$\displaystyle\min\sum_{n=0}^{N}\Big{\|}p_{h,s}^{n}-\sum_{j=1}^{M}(p_{h,s}^{n},\psi_{j})\psi_{j}\Big{\|}^{2}$$

$$\displaystyle\text{subject to }(\psi_{i},\psi_{j})=\delta_{ij}\quad\mbox{for $i,j=1,\ldots,M$}.$$

$$\mathbb{C}{\vec{a}}_i={\lambda }_i{\vec{a}}_i.$$

$$\mathbb{D}{\vec{b}}_i={\sigma }_i{\vec{b}}_i.$$

$${\vec{\phi }}_i=\frac{1}{\sqrt{{\lambda }_i}}\mathbb{A}{\vec{a}}_i,i=1,\dots ,R.$$

$${\vec{\psi }}_i=\frac{1}{\sqrt{{\sigma }_i}}\mathbb{B}{\vec{b}}_i,i=1,\dots ,M.$$

$$\displaystyle\Big{(}\frac{u^{n+1}_{R}-u^{n}_{R}}{\Delta t},\varphi\Big{)}+b^{\ast}(u_{R}^{n},u^{n+1}_{R},\varphi)+\nu(\nabla u^{n+1}_{R},\nabla\varphi)$$

Full text

{AMS}

65M12, 65M15, 65M60

1 Introduction

We consider the Navier-Stokes equations (NSE) with no-slip boundary conditions:

[TABLE]

Here $u$ is the velocity, $f$ is the known body force, $p$ is the pressure, and $\nu$ is the kinematic viscosity. For the past three decades, reduced order models (ROMs) have been successfully used in the numerical simulation of fluid flows modeled by the NSE (1)  [3, 10, 16, 17, 18, 19, 27, 30, 31, 33]. The ROM construction is similar to the full finite element approximation except we seek a solution in a low dimensional ROM space $X_{R}$ using the basis $\{\varphi_{i}\}_{i=1}^{R}$. These basis functions are often assumed to be weakly divergence-free. This assumption holds true, for example, if the ROM basis functions are constructed from data from a NSE discretization with finite element velocity-pressure pairs that satisfy the inf-sup (Ladyzhenskaya-Babuska-Brezzi (LBB)) condition. In this case, the pressure drops out from the ROM, which yields approximations only for the velocity field: such as for the backward Euler method

[TABLE]

where $b^{\ast}(w,u,v):=\frac{1}{2}(w\cdot\nabla u,v)-\frac{1}{2}(w\cdot\nabla v,u),\ \ \forall u,v,w\in[H^{1}(\Omega)]^{d}\,,$ and the superscript denotes the timestep number.

We emphasize, however, that even when the pressure is not required in the ROM, one may still need a ROM pressure approximation. This happens, for example, in fluid-structure interaction problems, if drag and lift coefficients need to be computed, or if the residual has to be calculated [7]. Another practical issue with velocity only ROMs is that internal (industrial) flows will often have reliable pressure data, but little to no velocity data. A velocity only ROM will be unable to incorporate pressure data to improve accuracy, calibrate the model, or check if a control loop is functioning properly.

When a ROM pressure approximation is required, there are two main approaches that are currently used:

(I) Inf-Sup/LBB Condition:  In the first approach, the velocity and pressure ROM approximations satisfy the inf-sup/LBB condition:

[TABLE]

This approach has been extensively developed in the reduced basis method (RBM) community over the past decade [17, 30]. This approach yields accurate ROM approximations for both velocity and pressure and eliminates the spurious numerical instabilities in the pressure approximation that are often generated by ROMs that do not satisfy the inf-sup condition. Furthermore, rigorous error estimates are proven for the LBB conforming ROM approximations. The RBM has been successfully used in numerous scientific and engineering applications [17, 30]. However, enforcing the inf-sup condition (3) is significantly more challenging for ROMs than for finite elements. Indeed, in the finite element context, the approximation spaces (e.g., piecewise quadratic for the velocity and piecewise linear for the pressure, i.e., the Taylor-Hood element) are specified beforehand and the corresponding discrete inf-sup condition can be investigated a priori. In the ROM context, on the other hand, the approximation spaces are problem-dependent – they are known only after the underlying finite element simulations (or the actual physical experiments) have been carried out. Thus, in a ROM context, the inf-sup condition needs to be enforced for each problem separately. In the RBM context, this is generally achieved by enriching the ROM basis with supremizers, which need to be computed in the offline stage. Thus, in realistic fluid flow applications (e.g., the NSE at high Reynolds numbers), enforcing the inf-sup condition can be prohibitively expensive (see, e.g., Sections 4.2.2 and 4.2.3 in [4]).

(II) Pressure Poisson Equation:  In the second approach to generate ROM approximation for the pressure, the available ROM velocity approximation is used to solve a pressure Poisson equation for the ROM pressure approximation

[TABLE]

which is obtained by taking the divergence of the NSE (1). This approach has been used in, e.g., [1, 7, 28]. We note that this approach faces several significant challenges: We emphasize that the Poisson equation (4) is not valid anymore if the ROM basis functions are not weakly divergence-free. This is the case, for example, if the ROM basis functions are built from data from NSE discretizations with finite element velocity-pressure pairs that do not satisfy the inf-sup/LBB condition, e.g., when the artificial compression, penalty, or projection methods are used [14]. Furthermore, the boundary conditions for (4) are not clear. Finally, the numerical investigation in [7] showed that even when weakly divergence-free snapshots were used, the ROMs that solve the pressure Poisson equation (4) were less competitive in terms of numerical accuracy and computational efficiency.

(III) Novel AC-ROM:  The artificial compression (AC) method and related approaches (e.g., the penalty and projection methods) have found significant success in the CFD community  [8, 11, 14]. The main idea in the AC method is to replace the incompressibility condition in the NSE with an artificial compression condition. Thus, the AC method decouples the velocity and pressure computations, which results in significant savings in execution time and storage. Furthermore, since the velocity and pressure computations are decoupled, the AC method allows the use of finite element pairs that do not satisfy the inf-sup/LBB condition [24]. (We also note that, because the incompressibility condition is not satisfied exactly, the AC method yields velocity fields that are not weakly divergence-free.)

In this paper, we develop a novel AC-ROM that employs AC to decouple the velocity-pressure ROM approximations: The fully discrete algorithm for the AC-ROM algorithm can be written as:

[TABLE]

where $(\{\varphi_{i}\}_{i=1}^{R},\{\psi_{k}\}_{k=1}^{M})$ is the ROM basis for the ROM space $(X_{R},Q_{M})$. The new AC-ROM (5a)–(5b) has several significant advantages over the approaches (I) and (II):

• •

The AC-ROM does not require that the velocity-pressure ROM spaces satisfy the inf-sup/LBB condition, thus avoiding the challenges encountered in approach (I).

• •

The AC-ROM basis functions are constructed from data that do not have to be weakly-divergence free, such as those from NSE discretizations with the artificial compression, penalty, or projection methods. Thus, the AC-ROM avoids the challenges faced by approach (II).

The rest of the paper is organized as follows: In Section 2, we introduce some notation. In Section 3, we describe the proper orthogonal decomposition, which we use to construct the ROM basis. In Section 4 and Section 5, we prove the stability and an error estimate of the AC-ROM (5a)–(5b), respectively. In Section 6, we investigate numerically the new AC-ROM in the simulation of a two-dimensional flow between offset cylinders. Finally, in Section 7, we draw conclusions and outline future research directions.

2 Notation and preliminaries

We denote by $\|\cdot\|$ and $(\cdot,\cdot)$ the $L^{2}(\Omega)$ norm and inner product, respectively, and by $\|\cdot\|_{L^{p}}$ and $\|\cdot\|_{W_{p}^{k}}$ the $L^{p}(\Omega)$ and Sobolev $W^{k}_{p}(\Omega)$ norms, respectively. $H^{k}(\Omega)=W_{2}^{k}(\Omega)$ with norm $\|\cdot\|_{k}$. For a function $v(x,t)$ that is well defined on $\Omega\times[0,T]$, we define the norms

[TABLE]

The space $H^{-1}(\Omega)$ denotes the dual space of bounded linear functionals defined on $H^{1}_{0}(\Omega)=\{v\in H^{1}(\Omega)\,:\,v=0\mbox{ on }\partial\Omega\}$; this space is equipped with the norm

[TABLE]

The solutions spaces $X$ for the velocity and $Q$ for the pressure are respectively defined as

[TABLE]

A weak formulation of the NSE is given as follows: find $u:(0,T]\rightarrow X$ and $p:(0,T]\rightarrow Q$ such that, for almost all $t\in(0,T]$, satisfy

[TABLE]

We denote conforming velocity and pressure finite element spaces based on a regular triangulation of $\Omega$ having maximum triangle diameter $h$ by $X_{h}\subset X$ and $Q_{h}\subset Q.$ We also assume that the finite element spaces satisfy the approximation properties

[TABLE]

where $C$ is a positive constant that is independent of $h$.

We define the trilinear form

[TABLE]

and the explicitly skew-symmetric trilinear form given by

[TABLE]

which satisfies the bound [24]

[TABLE]

To ensure the uniqueness of the NSE solution and ensure that standard finite element error estimates hold, we make the following regularity assumptions on the data and true solution [24]:

Assumption 1**.**

*In (6) we assume that $u^{0}\in X$, $f\in L^{2}(0,T;L^{2}(\Omega))$, $u\in L^{\infty}(0,T;L^{2}(\Omega))\cap L^{4}(0,T;H^{s+1}(\Omega))\cap H^{1}(0,T;H^{s+1}(\Omega))\cap H^{2}(0,T;L^{2}(\Omega))$, and $p\in L^{\infty}(0,T;Q\cap H^{k}(\Omega))$. *

We assume the following error estimate for the finite element solution of (6) used to compute the velocity and pressure snapshots:

Assumption 2**.**

We assume that the finite element errors satisfy the following error estimates

[TABLE]

Remark 3**.**

*Error estimates of this form have been proven for varying amounts of regularity on the continuous solution $u$ and $p$. Some examples include the incremental pressure correction schemes in [13] and chapter 7 of [29]. *

The full space and time model on which we base our method is a backward Euler based artificial compression scheme with a Taylor-Hood spatial discretization, i.e., $P^{s}-P^{s-1}$ with $s\geq 2$. Given $u^{0}_{h}\in X_{h}$, $p^{0}_{h}\in Q_{h}$ for $n=0,1,2,\ldots,N-1$, find $u^{n+1}_{h}\in X_{h}$ and $p_{h}^{n+1}\in Q_{h}$ satisfying

[TABLE]

3 Proper Orthogonal Decomposition

In this section we briefly describe the POD method and apply it to the previously stated artificial compression algorithm. A more detailed description of this method can be found in [23].

Given a positive integer $N$, let $0=t_{0}<t_{1}<\cdots<t_{N}=T$ denote a uniform partition of the time interval $[0,T]$. Denote by $u_{h,S}^{n}(x)\in X_{h}$, $p_{h,S}^{n}(x)\in Q_{h}$, $n=0,\ldots,N$, the finite element solution to (9) evaluated at $t=t_{n}$, $n=1,\ldots,N$.

We denote by ${u}_{S}^{n}$ and ${p}_{S}^{n}$ the vector of coefficients corresponding to the finite element functions $u_{h,S}^{n}(x)$ and $p_{h,S}^{n}$. We then define the velocity snapshot matrix $\mathbb{A}$ and pressure snapshot matrix $\mathbb{B}$ as

[TABLE]

i.e., the columns of $\mathbb{A}$ and $\mathbb{B}$ are the finite element coefficient vectors corresponding to the discrete snapshots. The POD method then seeks a low-dimensional basis

[TABLE]

which can approximate the snapshot data. Let $\delta_{ij}$ denote the Kronecker delta. These bases can be determined by solving the constrained minimization problems

[TABLE]

and

[TABLE]

Defining the correlation matrices $\mathbb{C}=\mathbb{A}^{T}\mathbb{M}\mathbb{A}$ and $\mathbb{D}=\mathbb{B}^{T}\mathbb{M}\mathbb{B}$, where $\mathbb{M}$ denotes the finite element mass matrix, these problems can then be solved by considering the eigenvalue problems

[TABLE]

and

[TABLE]

It can then be shown the POD basis functions will be given by

[TABLE]

and

[TABLE]

Using this POD basis we can now construct the AC-ROM algorithm. The construction is similar to the full finite element approximation except we seek a solution in the POD space $(X_{R},Q_{M})$ using the basis $(\{\varphi_{i}\}_{i=1}^{R},\{\psi_{k}\}_{k=1}^{M})$. The fully discrete algorithm for the AC-ROM algorithm can be written as:

[TABLE]

4 Stability

In this section we prove the unconditional, nonlinear, longtime stability of the AC-ROM algorithm.

Theorem 4** (Unconditional Stability of AC-ROM).**

For any $n$, we have the energy equality

[TABLE]

and energy inequality

[TABLE]

Proof.

Let $\varphi=2\Delta{t}u_{R}^{n+1}$ and $\psi=2\Delta{t}p_{M}^{n+1}$ in (12). By the polarization identity and skew-symmetry of the nonlinearity, we have

[TABLE]

Adding the two equations gives

[TABLE]

Summing from $n=0~{}\text{to}~{}N$ gives the energy equality above. By definition of the dual norm and Young’s inequality, we have the energy inequality

[TABLE]

proving unconditional stability.

5 Error Analysis

Next we provide an error analysis for the AC-ROM scheme. We begin by stating preliminary results.

Let $\mathbb{S}_{R}=(\nabla\varphi_{i},\nabla\varphi_{j})_{L^{2}}$ be the POD stiffness matrix and let $|||\cdot|||_{2}$ denote the matrix $2$-norm. It was shown in [23] that this POD basis satisfies the following inverse inequality.

Lemma 5** (POD inverse estimate).**

[TABLE]

The norm $|||\mathbb{S}_{R}|||_{2}$ on the right hand side of (13) depends on the choice of the POD basis with no universal pattern of growth with R (their number). Since $R$ is small, $|||\mathbb{S}_{R}|||_{2}$ can be precomputed giving a precise number for the right hand side of (13).

We define the $L^{2}$ projection into the velocity space $X_{R}$, and the pressure space $Q_{M}$ as follows.

Definition 1**.**

Let $P_{R}:L^{2}(\Omega)\rightarrow X_{R}$ and $\chi_{M}:L^{2}(\Omega)\rightarrow Q_{M}$ such that

[TABLE]

The following lemmas, proven in [23, 32], provide bounds for the error between the snapshots and their projections onto the POD space.

Lemma 6**.**

[$L^{2}$ POD projection error] With $\lambda_{i}$ the eigenvalues of $\mathbb{C}=\mathbb{A}^{T}\mathbb{M}\mathbb{A}$, we have

[TABLE]

Lemma 7**.**

[$H^{1}$ POD projection error] We have

[TABLE]

The following error estimates then follow easily for the $L^{2}$ projection error into the velocity space $X_{R}$ using the techniques in [23, 32].

Lemma 8**.**

For any $u^{n}\in V$ the $L^{2}$ projection error into $X_{R}$ satisfies the following estimates

[TABLE]

Similarly for the $L^{2}$ projection into the pressure space the following can be proven.

Lemma 9**.**

For any $p^{n}\in Q$ the $L^{2}$ projection error satisfies the following estimates

[TABLE]

Let $e_{u}$ and $e_{p}$ denote the error between the true velocity and pressure solution and their POD approximations respectively, For the error analysis we split the error for the velocity and the pressure using the $L^{2}$ projections into the spaces $X_{R},Q_{M}$

[TABLE]

We will see in Theorem 12 that the convergence rate faces order reduction by a power of $\Delta t^{-1}$ term appearing in the error bound. This occurs due to the term $(\nabla\cdot\eta^{n+1},\pi_{M}^{n+1})$ arising from the continuity equation. Due to the fact the AC-ROM scheme proposed in this paper does not require the ROM velocity-pressure spaces to satisfy the $LBB_{h}$ condition, this order reduction cannot be eliminated via the usual Stokes projection. However, we will show in Theorem 12 that even if the basis does not satisfy the $LBB_{h}$ condition this order reduction in the convergence rate will be improved by a multiplicative constant with size dependent upon the quality of the basis.

To this end, we consider the subspace

[TABLE]

and recall from [9] the strengthened Cauchy-Buniakowskii-Schwarz (CBS) inequality commonly used in the analysis of multilevel methods [2, 9].

Lemma 10**.**

Given a Hilbert space V and two finite dimensional subspaces $V_{1}\subset V$ and $V_{2}\subset V$ with trivial intersection:

[TABLE]

then there exists $0\leq\alpha<1$ such that

[TABLE]

Considering $X^{div}_{R}$ and $Q_{M}$, we are interested in computing the exact constant $\alpha$ between these spaces. This is equivalent to finding the first principal angle defined as

[TABLE]

with $0\leq\theta_{1}\leq\frac{\pi}{2}$.

The problem of computing angles between subspaces was introduced by Jordan in 1875 [21] and studied by Friedrichs in 1937 [12]. Recently, principal angles were used to improve the accuracy of reduced basis schemes for optimization problems in [26]. They can be calculated using either QR factorization or SVD of the orthogonal bases of the spaces in Lemma 10, as outlined in [34]. More efficient and stable schemes for calculating principal angles were also developed in [22]. We note that due to the relatively small size of the pressure and velocity reduced basis, the QR or SVD approach is sufficient in this setting. This procedure will be briefly outlined in section 6.

Using the strengthened CBS inequality, we get the following bound on the error term arising from the continuity equation.

Lemma 11**.**

Let $u^{n+1}=u(t^{n+1},x)$ be the exact solution of the NSE and let $\eta^{n+1}=u^{n+1}-P_{R}(u^{n+1})$ denote the projection error. Defining $\alpha=\cos(\theta_{1})$, where $\theta_{1}$ is given in (20) the following bound holds

[TABLE]

Proof.

Since $u^{n+1}$ is the exact solution to the NSE it follows that $\nabla\cdot u^{n+1}=0$ and therefore $(\nabla\cdot u^{n+1},\psi)=0\ \forall\psi\in Q_{M}$. This gives

[TABLE]

It then follows from the fact that $\nabla\cdot P_{R}u^{n+1}\in X_{R}^{div}$, Lemma 10, and $\nabla\cdot u^{n+1}=0$.

[TABLE]

We are now ready to state the full error estimate.

Theorem 12**.**

Consider AC-ROM (12) and the partition $0=t_{0}<t_{1}<\cdots<t_{N}=T$ used in section 3. Letting $\epsilon=\mathcal{O}(\Delta{t})$ it then holds

[TABLE]

Proof.

The weak solution of the NSE satisfies

[TABLE]

[TABLE]

where

[TABLE]

Now subtracting (12a) from (21) and (12b) from (22) we have

[TABLE]

and

[TABLE]

Setting $\varphi=2\Delta t\xi_{R}^{n+1}$ and $\psi=2\Delta t\pi_{M}^{n+1}$, we use the fact that $\left(\frac{\eta^{n+1}-\eta^{n}}{\Delta t},\xi_{R}^{n+1}\right)=0$ and $\left(\frac{\kappa^{n+1}-\kappa^{n}}{\Delta t},\pi_{M}^{n+1}\right)=0$ by the definition of the $L^{2}$ projection. Adding (24) to (25) and using the polarization identity yields

[TABLE]

By Poincaré and Young’s inequality we bound the first two terms on the right hand side of (26)

[TABLE]

For the third term on the right of (26), adding and subtracting $2\Delta t(\nabla\cdot\eta^{n+1},\pi_{M}^{n})$, applying Young’s inequality, and Lemma 11 yields

[TABLE]

Next, for the nonlinear terms we add and subtract $b^{*}\left(u_{R}^{n},u^{n+1},\xi^{n+1}_{R}\right)$ and

$b^{*}\left(u^{n},u^{n+1},\xi^{n+1}_{R}\right)$. This yields, by skew-symmetry,

[TABLE]

The first three nonlinear terms are now bounded using the Sobolev Imbedding Theorem, Young’s inequality, (7), and (8)

[TABLE]

Dealing with the consistency terms, by Taylor’s Theorem and Young’s inequality we have

[TABLE]

and, by adding and subtracting by $2\Delta{t}\tau_{p}(p^{n+1};\pi^{n}_{M})$, we have

[TABLE]

Letting $\delta_{1}=\delta_{2}=\delta_{5}=\delta_{6}=\delta_{7}=\delta_{8}=\delta_{10}=\frac{1}{14}$, $\delta_{3}=\delta_{11}=\frac{1}{4}$, $\delta_{4}=\delta_{12}=\frac{1}{2}$, $\delta_{9}=1$, and rearranging/combining terms we have

[TABLE]

We note by Theorem 4 and the Cauchy Schwartz inequality that it follows

[TABLE]

By Theorem 4, Assumption 1, (34), the fact that $\epsilon=\mathcal{O}(\Delta{t})$, combining all inequalities, taking a maximum $\widetilde{C}$ over all constants, and summing from $n=0$ to $N$ yields

[TABLE]

Therefore, by a discrete Grönwall inequality,

[TABLE]

By the triangle inequality we have $\|e_{u}^{n+1}\|^{2}\leq 2(\|\eta^{n+1}\|^{2}+\|\xi_{R}^{n+1}\|^{2})$, as well as

$\|e_{p}^{n+1}\|^{2}\leq 2(\|\kappa^{n+1}\|^{2}+\|\pi_{M}^{n+1}\|^{2})$. Applying this and and taking a maximum among constants the result follows.

We see in Theorem 12 that whether or not $\Delta{t}^{-1}$ order reduction occurs depends upon the constant $\alpha$. If $\alpha^{2}<<\Delta{t}$ the order reduction will be alleviated. On the other hand for $\alpha^{2}>>\Delta{t}$ there will be no improvement in the convergence.

In [5, 6] it was determined for a $P^{1}-P^{0}$ finite element pair that $\alpha=\sqrt{\frac{3}{2}}h$. To our knowledge no similar calculations have been done for other finite element pairs such as Taylor-Hood. It is unclear what impact the choice of offline basis functions has on the principal angle for the POD spaces. Additionally the POD basis functions will be problem dependent. It is possible that certain problem settings yield large principal angles independent of the choice of finite element basis functions used in the offline stage.

6 Numerical Experiments

In this section, we perform a numerical investigation of the new AC-ROM algorithm (12). First, we show that the AC-ROM algorithm yields accurate velocity and pressure approximations without enforcing the LBB condition or requiring weakly divergence free snapshots. Then, we illustrate numerically the theoretical scalings proved in Theorem 12. In particular, we show that the AC-ROM algorithm yields first order scalings with respect to the time step, $\Delta t$. All computations are done using the FEniCS software suite [25] and all meshes generated via the built in meshing package mshr.

6.1 Problem Setting

For the numerical experiments we consider the two-dimensional flow between offset cylinders used in [15, 20]. The domain is a disk with a smaller off-center disc inside. Let $r_{1}=1$, $r_{2}=0.1$, $c_{1}=1/2$, and $c_{2}=0$; then, the domain is given by

[TABLE]

The viscosity is $\nu=\frac{1}{100}$ and the body force is given by

[TABLE]

The POD basis is computed from snapshots of the finite element discretization of a backward Euler artificial compression scheme.

Remark 13**.**

*We emphasize that, since the snapshots are generated using an artificial compression scheme they will not be weakly divergence free. This is clearly illustrated in Fig. 6.3, where we plot the magnitude of the divergence of the velocity field. We note that ROMs based on the pressure Poisson equation (i.e., ROMs in approach II in Section 1) cannot be used when the snapshots are not weakly divergence free. *

For the offline calculation, the flow is initialized at rest ($u_{h}^{0}\equiv 0$ and $p_{h}^{0}\equiv 0$). We discretized in space via the $P^{2}$-$P^{1}$ Taylor-Hood element pair. The spaces $X_{h}$ and $Q_{h}$ had 114,224 and 14,421 degrees of freedom respectively. We took $\Delta t=2.5e-4$ and $\varepsilon=1e-6$. The mesh is shown in Fig. 4. The no-slip, no-penetration boundary conditions are imposed on both cylinders. The flow developed into an almost periodic flow after $t=12$. Velocity and pressure snapshots were taken for every $t\in[12,16]$. The resulting singular values are shown in Fig. 5. The POD modes corresponding to the six largest singular values for velocity (resp. pressure) are shown in Fig. 1 (resp. Fig. 2).

Remark 14**.**

*We emphasize that the new AC-ROM uses the same number of velocity and pressure basis functions, i.e., $R=M$ in (12). Thus, we expect that the ROM LBB condition (3) is not satisfied. This shows that the new AC-ROM avoids the ROM LBB condition, which is generally prohibitively expensive for the RBM methods in approach I of Section 1 when those are used in realistic flows (see, e.g., Sections 4.2.2 and 4.2.3 in [4]). *

The force due to drag is the force exerted by the smaller cylinder against the main flow, which is counterclockwise. We calculated this as the line integral of the stress tensor around the smaller cylinder dotted with $(0,-1)$. The force due to lift is the line integral of the stress tensor around the smaller cylinder dotted with $(1,0)$.

With the stress tensor $\tau=\left(\nabla u+(\nabla u)^{T}\right)-pI$, and $\Gamma_{small}$ the boundary restricted to the inner cylinder, these quantities are

[TABLE]

[TABLE]

6.2 Numerical Investigation

We compare the kinetic energy, force due to drag, and force due to lift of the ROM simulations with $R=M=3,5$ and $7$ with the offline simulation in Fig. 6. $R\geq 7$ appears sufficient to capture the kinetic energy, lift and drag accurately. Again, this is in spite of the fact that the LBB condition is not satisfied due to using an equal number of pressure and velocity modes (see Remark 14).

Next, to illustrate numerically the theoretical scalings proved in Theorem 12, we show that the new AC-ROM algorithm yields first order scalings with respect to the time step, $\Delta t$. To test convergence with respect to $\Delta t$, we take $R=50$ using the snapshots for $t\in[12,16]$ with the offline simulation done from $t=12$ to $t=12.24$. $\Delta t$ ranges from $1.6e-3$ to $2.5e-4$, which was the stepsize from the offline simulation. The error is measured by comparing the $u_{R}$ to the corresponding offline solution $u_{h}$. The relative $l^{2}L^{2}$ errors that are shown in Fig. 8 verify the $\mathcal{O}(\Delta t)$ convergence proven in Theorem 12.

We briefly outline the process of computing the first principal angle between the spaces $X_{R}^{div}$ and $Q_{M}$. Let $\{\nabla\cdot\varphi^{orth}_{i}\}_{i=1}^{R}$ denote the orthonormalized basis of $X_{R}^{div}$ (19). We consider the matrices $\mathbb{Q}=[\psi_{1},\psi_{2},\ldots\psi_{M}]$ and $\mathbb{X}=[\nabla\cdot\varphi^{orth}_{2},\nabla\cdot\varphi^{orth}_{2},\ldots\nabla\cdot\varphi^{orth}_{R}]$. Multiplying these two matrices and taking the SVD gives

[TABLE]

The first principal angle will then be given in terms of the first nonzero entry of $\Sigma$, by $\theta_{1}=\arccos(\sigma_{1})$. We measured the influence of the principal angle between the velocity and pressure POD basis using the method outlined above. The results are shown in Fig. 7. $\alpha^{2}$ begins near machine precision and seems to plateau around $1e-2$ when adding more basis functions. This appears to match up with our theoretical results and explains why we do not observe an order reduction in our numerical investigation.

7 Conclusions and Outlook

In this paper, we propose an artificial compression reduced order model (AC-ROM) for the numerical simulation of fluid flows. The new AC-ROM provides approximations for both the velocity and the pressure. Compared to the current ROMs that generate pressure approximations, the new AC-ROM has two main advantages: (i) it does not require the fulfillment of the inf-sup/LBB condition, which can be prohibitively expensive in current ROMs [7]; and (ii) it does not require weakly divergence-free snapshots, which allows it to work with snapshots generated with, e.g., artificial compression, penalty, or projection methods.

In Section 4, we prove the unconditional stability of the finite element discretization of the new AC-ROM. In Section 5, we prove an error estimate for the AC-ROM. In particular, we show that that it is possible to overcome the $\Delta t^{-1}$ order degradation due to lack of inf-sup stability,if the angle between the divergence of the velocity space and pressure space is sufficiently small.

In Section 6, we perform a numerical investigation of the new AC-ROM for a two dimensional flow between two offset cylinders. To generate the snapshots, we use the artificial compression method. Thus, the snapshots used in the AC-ROM construction are not weakly divergence-free, which is illustrated in 3. We also show that the velocity and pressure spaces of the new AC-ROM do not satisfy the LBB condition (see Fig. 7). In the numerical investigation of the new AC-ROM, we first show that the AC-ROM yields accurate velocity and pressure approximations. Specifically, in Fig. 6, we show that it provides energy, drag force, and lift force approximations that are close to the direct numerical simulation results. Next, to illustrate numerically the theoretical scalings proved in Section 5, we show that the new AC-ROM algorithm yields first order scalings with respect to the time step. Finally, in Fig. 7, we show that the constant multiplying the $\Delta t^{-1}$ term in the error estimate is extremely small. This may explain why we do not observe an order reduction in our numerical investigation.

One future research direction will be a further study of the principal angle and its impact on the convergence of the AC-ROM scheme. We will also investigate whether it plays a role in other popular schemes such as penalty methods. Another research direction that we plan to pursue is numerical stabilization of ROMs whose velocity-pressure ROM spaces do not satisfy the inf-sup/LBB condition. To our knowledge, numerical stabilization to account for the violation of the LBB condition in ROMs has been investigated only in [7]. In Section 3.2.3 in [7], it was shown that adding a pressure stabilizing/Petrov-Galerkin (PSPG) term to the ROM formulation yields better results than those produced by the other two velocity-pressure ROMs that were investigated.