Numerical analysis of a projection-based stabilized POD-ROM for incompressible flows

TL;DR

This paper introduces a new stabilized projection-based POD reduced order model for incompressible flows that improves computational efficiency and accuracy by relaxing traditional stability conditions and incorporating pressure modes.

Contribution

The paper presents a novel Local Projection Stabilization ROM that avoids the discrete inf-sup condition and allows use of non-divergence-free snapshots, enhancing CFD simulations.

Findings

The new LPS-ROM provides accurate velocity and pressure approximations.

Numerical results demonstrate improved stability and efficiency.

The method effectively models flow past a circular obstacle.

Abstract

In this paper, we propose a new stabilized projection-based POD-ROM for the numerical simulation of incompressible flows. The new method draws inspiration from successful numerical stabilization techniques used in the context of Finite Element (FE) methods, such as Local Projection Stabilization (LPS). In particular, the new LPS-ROM is a velocity-pressure ROM that uses pressure modes as well to compute the reduced order pressure, needed for instance in the computation of relevant quantities, such as drag and lift forces on bodies in the flow. The new LPS-ROM circumvents the standard discrete inf-sup condition for the POD velocity-pressure spaces, whose fulfillment can be rather expensive in realistic applications in Computational Fluid Dynamics (CFD). Also, the velocity modes does not have to be neither strongly nor weakly divergence-free, which allows to use snapshots generated for instance with penalty or projection-based stabilized methods. The numerical analysis of the fully Navier-Stokes discretization for the new LPS-ROM is presented, by mainly deriving the corresponding error estimates. Numerical studies are performed to discuss the accuracy and performance of the new LPS-ROM on a two-dimensional laminar unsteady flow past a circular obstacle.

Tables

Table 1. Table 1: Maximum drag coefficient c D m ​ a ​ x superscript subscript 𝑐 𝐷 𝑚 𝑎 𝑥 c_{D}^{max} , maximum lift coefficient c L m ​ a ​ x superscript subscript 𝑐 𝐿 𝑚 𝑎 𝑥 c_{L}^{max} , and Strouhal number for the FOM solution with do nothing BC (first row) and Dirichlet BC (second row) at the outlet, compared with reference intervals from [ 37 ] (third row).

$c_D^{max}$ $c_L^{max}$ $St$ FOM (Outlet do nothing BC) $3.22$ $0.96$ $0.303$ FOM (Outlet Dirichlet BC) $3.16$ $0.92$ $0.303$ Reference results from [37] $[3.22,3.24]$ $[0.99,1.01]$ $[0.295,0.305]$

Equations

$$\left\{\begin{array}[]{rcll}\partial_{t}{\boldsymbol{u}}+({\boldsymbol{u}}\cdot\nabla){\boldsymbol{u}}-\nu\Delta{\boldsymbol{u}}+\nabla p&=&{\boldsymbol{f}}&\,\,\,\mbox{in}\,\,\,\Omega\times(0,T),\\ \nabla\cdot{\boldsymbol{u}}&=&0&\,\,\,\mbox{in}\,\,\,\Omega\times(0,T),\\ {\boldsymbol{u}}&=&\bf{0}&\,\,\,\mbox{on}\,\,\,\Gamma\times(0,T),\\ {\boldsymbol{u}}({\boldsymbol{x}},0)&=&{\boldsymbol{u}}_{0}({\boldsymbol{x}})&\,\,\,\mbox{in}\,\,\,\Omega.\end{array}\right.$$

$$\mathbf{X}={\mathbf{H}}_0^1=[H_0^1(\Omega ){]}^d=\left\{\mathbf{v}\in [H^1(\Omega ){]}^d:\mathbf{v}=0\text{ on }\Gamma \right\}.$$

$$Q=L_0^2(\Omega )=\left\{q\in L^2(\Omega ):{\int }_{\Omega }qd\mathbf{x}=0\right\}.$$

$$\left\{\begin{array}[]{rcll}\displaystyle\frac{d}{dt}({\boldsymbol{u}},{\boldsymbol{v}})+b({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{v}})+\nu(\nabla{\boldsymbol{u}},\nabla{\boldsymbol{v}})-(p,\nabla\cdot{\boldsymbol{v}})&=&\langle{\boldsymbol{f}},{\boldsymbol{v}}\rangle&\forall{\boldsymbol{v}}\in{\boldsymbol{X}},\quad\text{in }\mathcal{D}^{\prime}(0,T),\\ (\nabla\cdot{\boldsymbol{u}},q)+\sigma(p,q)&=&0&\forall q\in Q,\quad\text{a.e. in }(0,T),\\ {\boldsymbol{u}}(0)&=&{\boldsymbol{u}}_{0},\end{array}\right.$$

$$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].$$

$$\left\{\begin{array}[]{rcll}\displaystyle\frac{d}{dt}({\boldsymbol{u}}_{h},{\boldsymbol{v}}_{h})+b({\boldsymbol{u}}_{h},{\boldsymbol{u}}_{h},{\boldsymbol{v}}_{h})+\nu(\nabla{\boldsymbol{u}}_{h},\nabla{\boldsymbol{v}}_{h})-(p_{h},\nabla\cdot{\boldsymbol{v}}_{h})&=&\langle{\boldsymbol{f}},{\boldsymbol{v}}_{h}\rangle&\text{in }\mathcal{D}^{\prime}(0,T),\\ (\nabla\cdot{\boldsymbol{u}}_{h},q_{h})+\sigma(p_{h},q_{h})&=&0&\text{a.e. in }(0,T),\\ {\boldsymbol{u}}_{h}(0)&=&{\boldsymbol{u}}_{0h},\end{array}\right.$$

$$\left\{\begin{array}[]{lll}&Y_{h}^{l}=V_{h}^{l}(\Omega)=\{v_{h}\in C^{0}(\overline{\Omega}):v_{h}|_{K}\in\mathbb{R}^{l}(K),\,\forall K\in{\cal T}_{h}\},&\\ \\ &{\mathbf{Y}}_{h}^{l}=[Y_{h}^{l}]^{d}=\{{\boldsymbol{v}}_{h}\in[C^{0}(\overline{\Omega})]^{d}:{\boldsymbol{v}}_{h}|_{K}\in[\mathbb{R}^{l}(K)]^{d},\,\forall K\in{\cal T}_{h}\},&\\ \\ &{\mathbf{X}}_{h}={\mathbf{Y}}_{h}^{l}\cap{\bf H}_{0}^{1}.&\end{array}\right.$$

$$(\cdot ,\cdot {)}_{\tau }:{\mathbf{L}}^2(\Omega )\times {\mathbf{L}}^2(\Omega )\to \mathbb{R},(\mathbf{v},\mathbf{w}{)}_{\tau }=\sum _{K\in {\mathcal{T}}_h}{\tau }_K(\mathbf{v},\mathbf{w}{)}_K,$$

$$\Vert \mathbf{v}{\Vert }_{\tau }=(\mathbf{v},\mathbf{v}{)}_{\tau }^{1/2},$$

$$\left\{\begin{array}[]{rcll}\displaystyle\frac{d}{dt}({\boldsymbol{u}}_{h},{\boldsymbol{v}}_{h})+b({\boldsymbol{u}}_{h},{\boldsymbol{u}}_{h},{\boldsymbol{v}}_{h})+\nu(\nabla{\boldsymbol{u}}_{h},\nabla{\boldsymbol{v}}_{h})-(p_{h},\nabla\cdot{\boldsymbol{v}}_{h})&=&\langle{\boldsymbol{f}},{\boldsymbol{v}}_{h}\rangle&\text{in }\mathcal{D}^{\prime}(0,T),\\ (\nabla\cdot{\boldsymbol{u}}_{h},q_{h})+(\Pi_{h}^{*}(\nabla p_{h}),\Pi_{h}^{*}(\nabla q_{h}))_{\tau}+\sigma(p_{h},q_{h})&=&0&\text{a.e. in }(0,T),\\ {\boldsymbol{u}}_{h}(0)&=&{\boldsymbol{u}}_{0h},\end{array}\right.$$

$$\left\{\begin{array}[]{rcl}\left(\displaystyle\frac{{\boldsymbol{u}}_{h}^{n+1}-{\boldsymbol{u}}_{h}^{n}}{\Delta t},{\boldsymbol{v}}_{h}\right)+b({\boldsymbol{u}}_{h}^{n+1},{\boldsymbol{u}}_{h}^{n+1},{\boldsymbol{v}}_{h})+\nu(\nabla{\boldsymbol{u}}_{h}^{n+1},\nabla{\boldsymbol{v}}_{h})&&\\ -(p_{h}^{n+1},\nabla\cdot{\boldsymbol{v}}_{h})&=&\langle{\boldsymbol{f}}^{n+1},{\boldsymbol{v}}_{h}\rangle,\\ \\ (\nabla\cdot{\boldsymbol{u}}_{h}^{n+1},q_{h})+(\Pi_{h}^{*}(\nabla p_{h}^{n+1}),\Pi_{h}^{*}(\nabla q_{h}))_{\tau}+\sigma(p_{h}^{n+1},q_{h})&=&0,\end{array}\right.$$

$$\min \Delta t\sum _{n=1}^N{\Vert {\mathbf{u}}_h^n-\sum _{i=1}^{r_v}{\left({\mathbf{u}}_h^n,{\mathbf{\phi }}_i\right)}_{{\mathcal{H}}_v}{\mathbf{\phi }}_i\Vert }_{{\mathcal{H}}_v}^2,$$

$$\min \Delta t\sum _{n=1}^N{\Vert p_h^n-\sum _{i=1}^{r_p}{\left(p_h^n,{\psi }_i\right)}_{{\mathcal{H}}_p}{\psi }_i\Vert }_{{\mathcal{H}}_p}^2,$$

$$K^v{\mathbf{a}}_i={\lambda }_i{\mathbf{a}}_i,\text{ for }1,\dots ,r_v,$$

$$K^p{\mathbf{b}}_i={\gamma }_i{\mathbf{b}}_i,\text{ for }1,\dots ,r_p,$$

$$K_{mn}^v=\Delta t{\left({\mathbf{u}}_h^n,{\mathbf{u}}_h^m\right)}_{{\mathcal{H}}_v},\text{ for }m,n=1,\dots ,N,$$

$$K_{mn}^p=\Delta t{\left(p_h^n,p_h^m\right)}_{{\mathcal{H}}_p},\text{ for }m,n=1,\dots ,N,$$

$${\mathbf{\phi }}_i(\cdot )=\frac{1}{\sqrt{{\lambda }_i}}\sqrt{\Delta t}\sum _{n=1}^N({\mathbf{a}}_i{)}_n{\mathbf{u}}_h^n,1\le i\le r_v,$$

$${\psi }_i(\cdot )=\frac{1}{\sqrt{{\gamma }_i}}\sqrt{\Delta t}\sum _{n=1}^N({\mathbf{b}}_i{)}_n{p}_h^n,1\le i\le r_p,$$

$$\Delta t\sum _{n=1}^N{\Vert {\mathbf{u}}_h^n-\sum _{i=1}^{r_v}{\left({\mathbf{u}}_h^n,{\mathbf{\phi }}_i\right)}_{{\mathcal{H}}_v}{\mathbf{\phi }}_i\Vert }_{{\mathcal{H}}_v}^2=\sum _{i=r_v+1}^{M_v}{\lambda }_i,$$

$$\Delta t\sum _{n=1}^N{\Vert p_h^n-\sum _{i=1}^{r_p}{\left(p_h^n,{\psi }_i\right)}_{{\mathcal{H}}_p}{\psi }_i\Vert }_{{\mathcal{H}}_p}^2=\sum _{i=r_p+1}^{M_p}{\gamma }_i,$$

$${\mathbf{X}}_r=\text{span}\left\{{\mathbf{\phi }}_1,\dots ,{\mathbf{\phi }}_r\right\}\subset {\mathbf{X}}_h,$$

$$Q_r=\text{span}\left\{{\psi }_1,\dots ,{\psi }_r\right\}\subset Q_h.$$

$$\mathbf{u}(\mathbf{x},t)\approx {\mathbf{u}}_r(\mathbf{x},t)=\sum _{i=1}^r{a}_i(t){\mathbf{\phi }}_i(\mathbf{x}),$$

$$p(\mathbf{x},t)\approx p_r(\mathbf{x},t)=\sum _{i=1}^r{b}_i(t){\psi }_i(\mathbf{x}),$$

$$\left\{\begin{array}[]{rcl}\left(\displaystyle\frac{{\boldsymbol{u}}_{r}^{n+1}-{\boldsymbol{u}}_{r}^{n}}{\Delta t},\boldsymbol{\varphi}\right)+b({\boldsymbol{u}}_{r}^{n+1},{\boldsymbol{u}}_{r}^{n+1},\boldsymbol{\varphi})+\nu(\nabla{\boldsymbol{u}}_{r}^{n+1},\nabla\boldsymbol{\varphi})&&\\ -(p_{r}^{n+1},\nabla\cdot\boldsymbol{\varphi})&=&\langle{\boldsymbol{f}}^{n+1},\boldsymbol{\varphi}\rangle,\\ \\ (\nabla\cdot{\boldsymbol{u}}_{r}^{n+1},\psi)+(\Pi_{h}^{*}(\nabla p_{r}^{n+1}),\Pi_{h}^{*}(\nabla\psi))_{\tau}+\sigma(p_{r}^{n+1},\psi)&=&0,\end{array}\right.$$

$$b(\mathbf{u},\mathbf{v},\mathbf{v})=0,$$

$$b(\mathbf{u},\mathbf{v},\mathbf{w})\le C\Vert \nabla \mathbf{u}{\Vert }_{{\mathbf{L}}^2}\Vert \nabla \mathbf{v}{\Vert }_{{\mathbf{L}}^2}\Vert \nabla \mathbf{w}{\Vert }_{{\mathbf{L}}^2}.$$

$$\Vert y{\Vert }_X+\Vert z{\Vert }_X\le C\Vert y+z{\Vert }_X\forall y\in Y,z\in Z.$$

$$|(y,z{)}_X|\le \alpha \Vert y{\Vert }_X\Vert z{\Vert }_X\forall y\in Y,z\in Z.$$

Full text

2010 Mathematics Subject Classification: Primary 65M12, 65M15, 65M60;

Secondary 76D03, 76D05.

Keywords: Navier–Stokes Equations, Projection Stabilization, Proper Orthogonal Decomposition, Reduced Order Models, Incompressible Flows, Numerical Analysis.

1 Introduction

Reduced Order Models (ROM) have been applied to numerical design in modern engineering as a tool that is wide-spreading in the scientific community in the recent years in order to solve complex realistic multi-parameters, multi-physics and multi-scale problems. Among the most popular ROM approaches, Proper Orthogonal Decomposition (POD) strategy provides optimal (from the energetic point of view) bases or modes to represent the dynamics from a given database (snapshots) obtained by a full order system. Onto these reduced bases, a Galerkin projection of the governing equations can be employed to obtain a low-order dynamical system for the bases coefficients. This has led researchers to apply POD-ROM to a variety of physical and engineering problems, including Computational Fluid Dynamics (CFD) problems in order to model the Navier–Stokes Equations (NSE), see e.g. [6, 10, 21, 28, 31, 40].

In this context, POD velocity modes are usually assumed to be at least weakly divergence-free. To be this assumption true, the POD velocity modes should be generated, for instance, by a Full Order Model (FOM) which consists in a NSE space discretization using inf-sup stable Finite Element (FE) for the velocity-pressure pair. In this way, the contribution of the pressure formally drops out from the ROM, which thus only approximates the velocity field through the POD velocity modes. Despite the appealing computational efficiency of only velocity ROM, there exist however important settings in which the pressure should be somehow considered. Indeed, the pressure is needed in many CFD applications, e.g. in the computation of relevant physical quantities, such as drag and lift forces on bodies in the flow, and for incompressible shear flows, as the mixing layer or the wake flow [34], where neglecting it may lead to large amplitude errors. On the other hand, note also that the weakly divergence-free property does not hold for many popular discretizations of the NSE. This is the case, for instance, of using equal order FE for the velocity-pressure pair, for which a suitable numerical stabilization becomes essential to circumvent the violation of the standard discrete inf-sup condition, as considered in this work. Furthermore, the pressure approximation allows the computation of the residual associated to the strong form of the NSE, often needed in stabilized discretizations (cf. [8]). Altogether, all these reasons pushed us to propose and fully analyze a new robust and stable ROM that directly incorporates an approximation of the pressure, driven by numerical stabilization motivations, and recovers it correctly, avoiding spurious pressure oscillations.

The new method draws inspiration from successful numerical stabilization techniques used in the context of FE methods, such as Local Projection Stabilization (LPS) methods (cf. [1, 2]). In particular, the new LPS-ROM is a coupled velocity-pressure ROM that uses pressure modes as well to compute the reduced order pressure. In order to avoid pressure instabilities sources, the new LPS-ROM circumvents the standard discrete inf-sup condition for the POD velocity-pressure spaces, whose fulfillment can be rather expensive in realistic applications in CFD, see for instance [7, 39], where an offline strategy based on the supremizer enrichment of the reduced velocity space has been proposed and applied in the POD context, adapted from the Reduced Basis (RB) method framework. Also, with respect to other proposals existing in the current ROM literature that provides velocity-pressure approximations, the velocity modes does not have to be neither strongly nor weakly divergence-free for the new LPS-ROM, which allows to use snapshots generated for instance with penalty or projection-based stabilized methods. This is not the case, for instance, of ROM based on a pressure Poisson equation approach (see, for instance, the first two methods investigated in [11] and also [39]), for which the velocity snapshots, and hence the POD velocity modes must be at least weakly divergence-free. This requirement also holds for the last method proposed and investigated in [11], which uses a residual-based stabilization mechanism in order to overcome a possible violation of the discrete inf-sup condition in the ROM framework by considering a decoupled approach for the reduced velocity-pressure pair.

The main contribution of the present paper has been to perform a stability and convergence analysis of the arising fully discrete LPS-ROM applied to the unsteady incompressible NSE, by mainly deriving the proof of a rigorous error estimate that considers all contributions: the spatial discretization error (due to the FE discretization), the temporal discretization error (due to the backward Euler method), and the POD truncation error. To the best of our knowledge, the LPS-ROM introduced is novel, and the numerical analysis for a stabilization-motivated ROM to take into account the violation of the discrete inf-sup condition cannot be found in the literature so far. Indeed, on the one hand a thorough numerical analysis has been recently performed for a stabilization-motivated ROM accounting for instabilities due to convection-dominated phenomena [4], which only uses velocity POD modes. On the other hand, only few numerical investigations of stabilization-motivated ROM can be found in the literature (cf. [8] and last method in [11]). In particular, in [11] the authors showed that adding a Pressure-Stabilizing Petrov–Galerkin (PSPG) term to account for stabilizing the violated discrete inf-sup condition and recover the reduced pressure provided more efficient and accurate results with respect to ROM based on a pressure Poisson equation approach, which would also require an ad-hoc treatment of the pressure boundary conditions. Parallel and independently to the current paper, a velocity-pressure ROM has been very recently proposed and analyzed in [20], which however relies on an artificial compression method to compute the reduced velocity-pressure approximations for which the pressure must be initialized.

Note that the detailed numerical analysis corroborated with numerical studies makes apparent an interesting link between the number of POD velocity-pressure modes used in the LPS-ROM and the angle $\theta$ between the space spanned by the divergence of the POD velocity modes and the POD pressure space. Indeed, for the numerical example proposed, where the same numerical stabilization technique used for the ROM is initially applied also to the FOM (LPS-FOM) to generate the snapshots, so that these latter are not weakly divergence-free, we have found that for small values of $r$, which is common in practice, the saturation constant $\alpha=cos(\theta)$ [13] is rather small, and this allows to ease the convergence order reduction due to the violation of the discrete inf-sup stability condition.

Numerical studies performed on a two-dimensional laminar unsteady flow past a circular obstacle have also been used to assess the accuracy and efficiency of the new LPS-ROM. Despite the fact that the discrete inf-sup condition is not fulfilled by the new LPS-ROM, using a small equal number of POD velocity-pressure modes already provides accurate approximations, close to the LPS-FOM results, and theoretical considerations suggested by the numerical analysis are recovered in practice.

The outline of the paper is as follows: In Section 2, we introduce the model problem and its continuous variational formulation for time-dependent NSE. In Section 3, we consider the FE-LPS full order discretization used to generate the snapshots for the online phase. Section 4 briefly describes the POD methodology and introduce the formulation of the new LPS-ROM for the NSE. The stability and error analysis for the full discretization (FE in space and backward Euler in time) of the new model is presented in Section 5. The numerical investigation of the new method is proposed in Section 6 for the simulation of a two-dimensional flow past a circular obstacle, in order to test on the one hand some theoretical predictions of the performed numerical analysis and to show on the other hand the accuracy and efficiency of the proposed method in preventing spurious pressure instabilities due to the violation of the discrete inf-sup condition for the reduced system. Finally, Section 7 presents the main conclusions of this work and ongoing research directions.

2 Time-dependent NSE: model problem and variational formulation

We introduce an Initial-Boundary Value Problem (IBVP) for the incompressible evolution NSE. For the sake of simplicity, we just impose the homogeneous Dirichlet boundary condition on the whole boundary.

Let $[0,T]$ be the time interval and $\Omega$ a bounded polyhedral domain in $\mathbb{R}^{d}$, $d=2$ or $3$, with a Lipschitz-continuous boundary $\Gamma=\partial\Omega$. The transient NSE for an incompressible fluid are given by:

Find ${\boldsymbol{u}}:\Omega\times(0,T)\longrightarrow\mathbb{R}^{d}$ and $p:\Omega\times(0,T)\longrightarrow\mathbb{R}$ such that:

[TABLE]

The unknowns are the velocity ${\boldsymbol{u}}({\boldsymbol{x}},t)$ and the pressure $p({\boldsymbol{x}},t)$ of the incompressible fluid. The data are the source term ${\boldsymbol{f}}({\boldsymbol{x}},t)$, which represents a body force per mass unit (typically the gravity), the kinematic viscosity $\nu$ of the fluid, which is a positive constant, and the initial velocity ${\boldsymbol{u}}_{0}({\boldsymbol{x}})$.

To define the weak formulation of problem (2.1), we need to introduce some useful notations for functional spaces [9]. We consider the Sobolev spaces $H^{s}(\Omega)$, $s\in\mathbb{R}$, $L^{p}(\Omega)$ and $W^{m,p}(\Omega)$, $m\in\mathbb{N}$, $1\leq p\leq\infty$. We shall use the following notation for vector-valued Sobolev spaces: ${\bf H}^{s}$, ${\bf L}^{p}$ and ${\bf W}^{m,p}$ respectively shall denote $[H^{s}(\Omega)]^{d}$, $[L^{p}(\Omega)]^{d}$ and $[W^{m,p}(\Omega)]^{d}$ (similarly for tensor spaces of dimension $d\times d$). Also, the parabolic Bochner function spaces $L^{p}(0,T;X)$ and $L^{p}(0,T;{\bf X})$, where $X$ (${\bf X}$) stands for a scalar (vector-valued) Sobolev space, shall be denoted by $L^{p}(X)$ and $L^{p}({\bf X})$, respectively. In order to give a variational formulation of problem (2.1), let us consider the velocity space:

[TABLE]

This is a closed linear subspace of ${\bf H}^{1}$ and thus a Hilbert space endowed with the ${\bf H}^{1}$-norm. Thanks to Poincaré inequality, the ${\bf H}^{1}$-norm is equivalent on ${\bf H}_{0}^{1}$ to the norm $\|{{\boldsymbol{v}}}\|_{{\bf H}_{0}^{1}}=\|{\nabla{\boldsymbol{v}}}\|_{{\bf L}^{2}}$. Also, let us consider the pressure space:

[TABLE]

Note that the null mean condition for the pressure is usually introduced in order to fix the constant the pressure is determined up through the formulation. However, this condition could be relaxed by adding a suitable penalty term to the variational formulation of (2.1) that allows to simply take $Q=L^{2}(\Omega)$ as pressure space.

We shall thus consider the following variational formulation of (2.1):

Given ${\boldsymbol{f}}\in L^{2}({\bf H}^{-1})$, find ${\boldsymbol{u}}:(0,T)\longrightarrow{\boldsymbol{X}}$, $p:(0,T)\longrightarrow Q$ such that

[TABLE]

where $(\cdot,\cdot)$ stands for the $L^{2}$-inner product in $\Omega$, $\langle\cdot,\cdot\rangle$ stands for the duality pairing between ${\boldsymbol{X}}$ and its dual ${\boldsymbol{X}}^{{}^{\prime}}={\bf H}^{-1}$, and $\mathcal{D}^{\prime}(0,T)$ is the space of distributions in $(0,T)$. The trilinear form $b$ is given by: for ${\boldsymbol{u}},\,{\boldsymbol{v}},\,{\boldsymbol{w}}\in{\boldsymbol{X}}$

[TABLE]

The term with factor $\sigma$ denotes the penalty term that permits to fix the constant the pressure is determined up through the formulation, for a small positive value of $\sigma$ (e.g., $\sigma=\mathcal{O}(10^{-6})$), and thus $Q=L^{2}(\Omega)$ in (2.2) and hereafter. Note that in the continuous model problem (2.1) the incompressibility condition is no more satisfied exactly, and thus we are going to search for a velocity field approximation that does not have to be neither strongly nor weakly divergence-free.

3 Finite element full order model

In order to give a FE approximation of (2.2), let $\{{\cal T}_{h}\}_{h>0}$ be a family of affine-equivalent, conforming (i.e., without hanging nodes) and regular triangulations of $\overline{\Omega}$, formed by triangles or quadrilaterals ($d=2$), tetrahedra or hexahedra ($d=3$). For any mesh cell $K\in{\cal T}_{h}$, its diameter will be denoted by $h_{K}$ and $h=\max_{K\in{\cal T}_{h}}h_{K}$. We consider ${\boldsymbol{X}}_{h}\subset{\boldsymbol{X}}$, $Q_{h}\subset Q$ being suitable FE spaces for velocity and pressure, respectively. The FE approximation of (2.2) can be written as follows:

Find $({\boldsymbol{u}}_{h},p_{h}):(0,T)\longrightarrow{\boldsymbol{X}}_{h}\times Q_{h}$ such that

[TABLE]

for any $({\boldsymbol{v}}_{h},q_{h})\in{\boldsymbol{X}}_{h}\times Q_{h}$, and the initial condition ${\boldsymbol{u}}_{0h}$ is some stable approximation to ${\boldsymbol{u}}_{0}$ in $L^{2}$-norm belonging to ${\boldsymbol{X}}_{h}$.

In order to circumvent the standard discrete inf-sup condition and thus use equal order interpolation for velocity and pressure, and also provide an extra-control on the high frequencies components of the pressure gradient that could lead to unstable discretizations, we introduce a filtered pressure stabilizing term of high-order. In this way, the considered FE method falls into the class of Local Projection Stabilization (LPS) methods (cf. [1, 2]). The stabilization effect is achieved by adding a least-square term that give a weighted control on the fluctuations of the pressure gradient, based upon a specific locally stable projection or interpolation operator on a continuous buffer space. This provides an efficient discretization with a reduced computational cost that keeps the same high-order accuracy with respect to standard projection-stabilized methods.

To describe this approach, we define hereafter the specific choice of FE spaces done both for the numerical analysis and practical computations in the present work. Given an integer $l\geq 2$ and a mesh cell $K\in{\cal T}_{h}$, denote by $\mathbb{R}^{l}(K)$ either $\mathbb{P}^{l}(K)$ (i.e., the space of Lagrange polynomials of degree $\leq l$, defined on $K$), if the grids are formed by triangles ($d=2$) or tetrahedra ($d=3$), or $\mathbb{Q}^{l}(K)$ (i.e., the space of Lagrange polynomials of degree $\leq l$ on each variable, defined on $K$), if the family of triangulations is formed by quadrilaterals ($d=2$) or hexahedra ($d=3$). We consider the following FE spaces for the velocity:

[TABLE]

Hereafter, ${\mathbf{Y}}_{h}^{l}$ (resp., $Y_{h}^{l}$) will constitute the discrete foreground vector-valued (resp., scalar) spaces in which we will work on. Also, $Q_{h}=Y_{h}^{l}\subset L^{2}$, since we use equal order FE. We define the scalar product:

[TABLE]

and its associated norm:

[TABLE]

where for any $K\in{\cal T}_{h}$, $\tau_{K}$ is in general a positive local stabilization parameter.

The considered FE-FOM is given by:

Find $({\boldsymbol{u}}_{h},p_{h}):(0,T)\longrightarrow{\boldsymbol{X}}_{h}\times Q_{h}$ such that

[TABLE]

for any $({\boldsymbol{v}}_{h},q_{h})\in{\boldsymbol{X}}_{h}\times Q_{h}$, where $\Pi_{h}^{*}=Id-\Pi_{h}$ is the “fluctuation operator”, being $Id$ the identity operator, and $\Pi_{h}$ some locally stable (in $L^{2}$-norm) projection or interpolation operator from ${\bf L}^{2}(\Omega)$ on the foreground vector-valued space ${\mathbf{Y}}_{h}^{l-1}$ (also called “buffer space” in this context), satisfying optimal error estimates (cf. [17]). In practical implementations, we choose $\Pi_{h}$ as a Scott–Zhang-like [38] linear interpolation operator in the space ${\mathbf{Y}}_{h}^{l-1}$ (see [14], Sect. 4 for its construction), implemented by the software FreeFEM [23].

To state the full space-time discretization of the unsteady LPS-FOM (3.3), consider a positive integer number $N$ and define $\Delta t=T/N$, $t_{n}=n\Delta t$, $n=0,1,\ldots,N$. We compute the approximations ${\boldsymbol{u}}_{h}^{n}$, $p_{h}^{n}$ to ${\boldsymbol{u}}^{n}={\boldsymbol{u}}(\cdot,t_{n})$ and $p^{n}=p(\cdot,t_{n})$ by using, for simplicity of the analysis, a backward Euler scheme:

• •

Initialization. Set: ${\boldsymbol{u}}_{h}^{0}={\boldsymbol{u}}_{0h}.$

• •

Iteration. For $n=0,1,\ldots,N-1$: Given ${\boldsymbol{u}}_{h}^{n}\in{\mathbf{X}}_{h}$, find $({\boldsymbol{u}}_{h}^{n+1},p_{h}^{n+1})\in{\mathbf{X}}_{h}\times Q_{h}$ such that:

[TABLE]

for any $({\boldsymbol{v}}_{h},q_{h})\in{\mathbf{X}}_{h}\times Q_{h}$.

4 Proper orthogonal decomposition reduced order model

We briefly describe the POD method, following [29], and apply it to the projection-based stabilized FOM (3.4).

Let us consider the ensembles of velocity snapshots $\chi^{v}=\text{span}\left\{{\boldsymbol{u}}_{h}^{1},\ldots,{\boldsymbol{u}}_{h}^{N}\right\}$ and pressure snapshots $\chi^{p}=\text{span}\left\{p_{h}^{1},\ldots,p_{h}^{N}\right\}$, given by the FE solutions to (3.4) at time $t_{n}$, $n=1,\ldots,N$. The POD method seeks low-dimensional bases $\left\{\boldsymbol{\varphi}_{1},\ldots,\boldsymbol{\varphi}_{r_{v}}\right\}$ and $\left\{\psi_{1},\ldots,\psi_{r_{p}}\right\}$ in real Hilbert spaces $\mathcal{H}_{v}$, $\mathcal{H}_{p}$ that optimally approximate the velocity and pressure snapshots in the following sense:

[TABLE]

and

[TABLE]

subject to the conditions $\left(\boldsymbol{\varphi}_{j},\boldsymbol{\varphi}_{i}\right)_{\mathcal{H}_{v}}=\delta_{ij}$, $1\leq i,j\leq r_{v}$ and $\left(\psi_{j},\psi_{i}\right)_{\mathcal{H}_{p}}=\delta_{ij}$, $1\leq i,j\leq r_{p}$, where $\delta_{ij}$ is the Kronecker delta. To solve the optimization problems (4.1)-(4.2), one can respectively consider the eigenvalue problems:

[TABLE]

and

[TABLE]

where $K^{v},K^{p}\in\mathbb{R}^{N\times N}$ are the velocity, pressure snapshots correlation matrices, respectively with entries:

[TABLE]

and

[TABLE]

${\bf a}_{i},{\bf b}_{i}$ are the $i$-th eigenvector, and $\lambda_{i},\gamma_{i}$ are the associated eigenvalues. The eigenvalues are positive and sorted in descending order: $\lambda_{1}\geq\ldots\geq\lambda_{r_{v}}>0$ and $\gamma_{1}\geq\ldots\geq\gamma_{r_{p}}>0$. It can be shown that the solutions of (4.1)-(4.2), i.e. the POD velocity-pressure bases functions, are respectively given by:

[TABLE]

and

[TABLE]

where $({\bf a}_{i})_{n},({\bf b}_{i})_{n}$ are the $n$-th components of the eigenvectors ${\bf a}_{i},{\bf b}_{i}$, respectively. It can also be shown that the following POD projection error formulas hold [24, 29]:

[TABLE]

and

[TABLE]

where $M_{v},M_{p}$ are the rank of $\chi^{v}$ and $\chi^{p}$, respectively. Although $\mathcal{H}_{v},\mathcal{H}_{p}$ can be any real Hilbert spaces, in what follows we consider $\mathcal{H}_{v}={\bf L}^{2}$ and $\mathcal{H}_{p}=L^{2}$. Also, we are going to take the same number of velocity and pressure POD bases functions, i.e. $r_{v}=r_{p}=r$ in what follows. Thus, we expect that also the POD velocity-pressure spaces will not satisfy the standard discrete inf-sup condition, and the POD-Reduced Order Model (POD-ROM) we are going to consider must circumvent it. To overcome this restriction, we draw inspiration from the FOM (3.4) in order to construct the new projection-based stabilized POD-ROM.

We respectively consider the following velocity and pressure spaces for the POD setting:

[TABLE]

and

[TABLE]

Remark 4.1**.**

Since, as shown in (4.5), the POD velocity modes are linear combinations of the velocity snapshots, the POD velocity modes satisfy the boundary conditions in (2.1). This is because of the particular choice we have made at the beginning to work with homogeneous Dirichlet boundary conditions. In general, one has to manipulate the velocity snapshots set. This is the case, for instance, of steady-state non-homogeneous Dirichlet boundary conditions, for which is preferable to consider a proper lift in order to generate POD velocity modes for the lifted velocity snapshots, satisfying homogeneous Dirichlet boundary conditions. This would lead to work with centered-trajectory method in the POD-ROM setting [22].

The standard Galerkin projection-based POD-ROM uses both Galerkin truncation and Galerkin projection. The former yields an approximation of the velocity and pressure fields by a linear combination of the corresponding truncated POD basis:

[TABLE]

and

[TABLE]

where $\left\{a_{i}(t)\right\}_{i=1}^{r}$ and $\left\{b_{i}(t)\right\}_{i=1}^{r}$ are the sought time-varying coefficients representing the POD-Galerkin velocity and pressure trajectories. Note that $r<<\mathcal{N}$, where $\mathcal{N}$ denotes the number of degrees of freedom (d.o.f.) of the equal order FE velocities-pressure in FOM (3.4). Replacing the velocity-pressure FE pair $({\boldsymbol{u}}_{h},p_{h})$ with $({\boldsymbol{u}}_{r},p_{r})$ in the FE approximation (3.4) and projecting the resulted equations onto the POD product space $({\boldsymbol{X}}_{r},Q_{r})$ using the POD basis $\left(\left\{\boldsymbol{\varphi}_{i}\right\}_{i=1}^{r},\left\{\psi_{i}\right\}_{i=1}^{r}\right)$, the full space-time discretization of the new projection-based stabilized POD-ROM reads as:

• •

Initialization. Set: ${\boldsymbol{u}}_{r}^{0}=\displaystyle\sum_{i=1}^{r}({\boldsymbol{u}}_{0h},\boldsymbol{\varphi}_{i})\boldsymbol{\varphi}_{i}.$

• •

Iteration. For $n=0,1,\ldots,N-1$: Given ${\boldsymbol{u}}_{r}^{n}\in{\boldsymbol{X}}_{r}$, find $({\boldsymbol{u}}_{r}^{n+1},p_{r}^{n+1})\in{\boldsymbol{X}}_{r}\times Q_{r}$ such that:

[TABLE]

for any $(\boldsymbol{\varphi},\psi)\in{\boldsymbol{X}}_{r}\times Q_{r}$.

An alternative time discretization could be given by the semi-implicit Euler method, where the trilinear form in (4.11) is discretized by $b({\boldsymbol{u}}_{r}^{n},{\boldsymbol{u}}_{r}^{n+1},\boldsymbol{\varphi})$. Note that considering a semi-implicit time discretization of the new projection-based stabilized POD-ROM is less costly from the computational point of view with respect to a fully implicit one, which yields a nonlinear algebraic system of equations to be solved. However, the numerical analysis will be performed in detail for the more technical case of the fully implicit time discretization given by (4.11).

5 Analysis of the projection-based stabilized POD-ROM

In this section, we perform the numerical analysis of the proposed unsteady POD-ROM (4.11), which we will call in the sequel LPS-ROM.

5.1 Technical background

This section provides some technical results that are required for the numerical analysis. Throughout the paper, we shall denote by $C$ a positive constant that may vary from a line to another, but which is always independent of the FE mesh size $h$, the FE velocity-pressure equal interpolation order $l$, the time step $\Delta t$, and the velocity, pressure eigenvalues $\lambda_{i}$, $\gamma_{i}$.

Lemma 5.1** (See Lemma 13 in [30]).**

For any function ${\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}}\in{\boldsymbol{X}}$, the skew-symmetric trilinear form $b(\cdot,\cdot,\cdot)$ satisfies:

[TABLE]

[TABLE]

Definition 5.2**.**

Let $X$ be a Hilbert space and $Y$, $Z$ two finite-dimensional subspaces of $X$ with intersection reduced to the zero function. The pair of finite-dimensional spaces $(Y,Z)$ is called to satisfy the saturation property if there exists a positive constant $C$ such that:

[TABLE]

Thus, the saturation property can be viewed as an inverse triangular inequality.

Lemma 5.3**.**

The saturation property is equivalent to the existence of a constant $\alpha<1$ such that:

[TABLE]

Actually, we may take $\alpha=1-2/C^{2}$ (see Remark 2 in [13]), and in the sequel we will call $\alpha$ the saturation constant. Then, we can interpret the saturation property in the sense that the angle between spaces $Y$ and $Z$, defined by:

[TABLE]

is uniformly bounded from below by a positive angle, and $\alpha=\cos(\theta)$.

Remark 5.4**.**

Lemma 5.3 will be essential in Theorem 5.14 to bound the error term coming from the continuity equation, which cannot be removed by using the standard Stokes projection since the reduced velocity-pressure spaces in the LPS-ROM (4.11) proposed in this work violate the standard discrete inf-sup condition.

In the numerical studies performed in Section 6 we will observe that the saturation constant $\alpha$, for the chosen numerical setup, starts with small values for small $r$ and seems to experience a flattening effect around $10^{-2}$ when adding more POD modes.

Note that the argument of saturation property has been used in [13] to develop a stabilized post-processing of the Galerkin FE solution of convection-dominated flows and very recently extended [5] to POD-ROM approximations to propose a cure for instabilities due to advection-dominance in POD solution to advection-diffusion-reaction equations [36].

To ensure error estimates in Theorem 5.14 (main result of the present paper), we make the following regularity assumption on the continuous solution:

Hypothesis 5.5**.**

We assume that the continuous solution $({\boldsymbol{u}},p)$ of the unsteady NSE (2.2) has augmented regularity, i.e. $({\boldsymbol{u}},p)\in L^{\infty}({\bf H}^{s+1})\times L^{\infty}(H^{s})$, $s\geq 1$, such that $\partial_{t}^{2}{\boldsymbol{u}}\in L^{2}({\bf L}^{2})$.

For the subsequent numerical analysis, we need the following technical hypothesis on the stabilization parameters $\tau_{K}$:

Hypothesis 5.6**.**

The stabilization parameters $\tau_{K}$ satisfy the following condition:

[TABLE]

Remark 5.7**.**

The question whether the stabilization parameters should depend on the spatial resolution of the underlying FE space, or on the number of POD basis functions used has been addressed in [22], by means of numerical analysis arguments. In that work, numerical investigations using both definitions suggested that the one based on estimates from the underlying FE discretization provides a better suppression of numerical oscillations, and thus guarantees a more effective numerical stabilization. For this reason, we make here assumption 5.6 on the stabilization parameters, which is also essential for the subsequent numerical analysis.

Definition 5.8**.**

Let ${\boldsymbol{v}}_{r}$ be the $L^{2}$-orthogonal projection of ${\boldsymbol{u}}$ on the reduced velocity space ${\boldsymbol{X}}_{r}$:

[TABLE]

and $z_{r}$ be the $L^{2}$-orthogonal projection of $p$ on the pressure space $Q_{r}$:

[TABLE]

We have the following error estimates for ${\boldsymbol{v}}_{r}$ and $z_{r}$ (see [25], Lemma 3.3):

Lemma 5.9** ($L^{2}$-projection error estimates).**

[TABLE]

[TABLE]

Lemma 5.10** ($H^{1}$-projection error estimates).**

[TABLE]

with $\|{S^{v}_{r}}\|_{2}$ denoting the $2$-norm of the stiffness velocity matrix with entries $[S^{v}_{r}]_{ij}=(\nabla\boldsymbol{\varphi}_{j},\nabla\boldsymbol{\varphi}_{i})$, $i,j=1,\ldots,r$.

[TABLE]

with $\|{S^{p}_{r}}\|_{2}$ denoting the $2$-norm of the stiffness pressure matrix with entries $[S^{p}_{r}]_{ij}=(\nabla\psi_{j},\nabla\psi_{i})$, $i,j=1,\ldots,r$.

The appearance of $\|{S^{v}_{r}}\|_{2}$, $\|{S^{p}_{r}}\|_{2}$ in Lemma 5.10 comes from the use of the POD inverse estimates (see [29], Lemma 2):

[TABLE]

5.2 Existence and stability results for LPS-ROM

We have the following existence and unconditional stability result for the LPS-ROM (4.11):

Theorem 5.11**.**

Problem (4.11) admits a solution that satisfies the following bound:

[TABLE]

for $k=0,\ldots,N$.

Proof. Problem (4.11) can be written as:

[TABLE]

for any $(\boldsymbol{\varphi},\psi)\in{\boldsymbol{X}}_{r}\times Q_{r}$, where $\widetilde{a}({\boldsymbol{u}}_{r}^{n+1},\boldsymbol{\varphi})=\Delta t^{-1}({\boldsymbol{u}}_{r}^{n+1},\boldsymbol{\varphi})+\nu(\nabla{\boldsymbol{u}}_{r}^{n+1},\nabla\boldsymbol{\varphi})$ and $\widetilde{{\boldsymbol{f}}}^{n+1}={\boldsymbol{f}}^{n+1}+\Delta t^{-1}({\boldsymbol{u}}_{r}^{n},\boldsymbol{\varphi})$. This problem fits into the same functional framework as for implicit discretizations of the steady NSE (see [15], for instance), since $\widetilde{a}$ is an inner product on space ${\boldsymbol{X}}$ that generates a norm equivalent to the ${\bf H}^{1}$-norm. Then, the existence of a solution follows from Brouwer’s fixed point theorem [9] as for the steady case.

• •

Velocity-pressure estimate.

To prove estimate (5.15), we set $\boldsymbol{\varphi}=2\Delta t{\boldsymbol{u}}_{r}^{n+1}$, $\psi=2\Delta tp_{r}^{n+1}$ in (4.11) and add both equations. Using the polarization identity:

[TABLE]

and noting that $b({\boldsymbol{u}}_{r}^{n+1},{\boldsymbol{u}}_{r}^{n+1},{\boldsymbol{u}}_{r}^{n+1})=0$ by (5.1), we obtain:

[TABLE]

By definition of the dual norm and Young’s inequality, from (5.17) we get:

[TABLE]

Then, the stability estimate (5.15) follows by summing (5.18) from $n=0$ to $k\leq N-1$. $\Box$

Remark 5.12**.**

The stability estimate (5.15), which makes apparent the estimate of the pressure stabilization term, guarantees an extra-control on the high frequencies of the pressure gradient that could lead to unstable discretization.

Remark 5.13**.**

An alternative unconditional pressure stability estimate for LPS-ROM (4.11) could be obtained if we do not add the penalty term with factor $\sigma$ to the variational formulation. In this case, when considering for instance Dirichlet boundary conditions on the whole boundary, the constant the pressure is determined up through the formulation could be fixed prescribing its value at a point of the boundary. When considering do nothing boundary conditions on a part of the boundary (outflow), the constant the pressure is determined up through the formulation is already fixed by the prescribed boundary conditions.

To do so, one sets $\psi=0$ in (4.11). For any $\boldsymbol{\varphi}\in{\boldsymbol{X}}_{r}$, this yields:

[TABLE]

Let $P_{r}^{n+1}=\displaystyle\sum_{k=0}^{n}\Delta tp_{r}^{k+1}$, then summation over the discrete times gives:

[TABLE]

Thus, we get:

[TABLE]

where we have applied triangle inequality, Cauchy–Schwarz inequality, the definition of the dual norm, and the standard estimate (5.2) for the convective term. So, if we define the norm:

[TABLE]

by Cauchy–Schwarz and triangle inequalities, from (5.19) we have:

[TABLE]

Using estimate (5.15) to bound the terms on the r.h.s. of (5.21), we obtain:

[TABLE]

which proves the unconditional stability of the time primitive of the pressure in the space induced by the norm (5.20), which we will call $\widetilde{Q}_{r}$ in the sequel. Indeed, due to the stability of ${\boldsymbol{u}}_{0h}$ in $L^{2}$-norm and the regularity of ${\boldsymbol{f}}$, we have:

[TABLE]

and we have denoted $P_{r}(t)=\displaystyle\int_{0}^{t}\widetilde{p}_{r}(s)\,ds$, being $\widetilde{p}_{r}$ the piecewise constant in time function that takes the value $p_{r}^{n+1}$ on $(t_{n},t_{n+1})$. Note that for the ROM pressure we can only obtain:

[TABLE]

which is a not uniform bound (with respect to $\Delta t$) in $L^{1}(\widetilde{Q}_{r})$ space of space-time functions (cf. [16], Remark 10.2), thus the technical trick of considering the time primitive of the pressure, as originally introduced in [16] for a FE-FOM, to prove its unconditional stability in this case. On the other side, when considering the penalty term with factor $\sigma>0$ to fix the constant the pressure is determined up through the formulation (for Dirichlet boundary conditions on the whole boundary), as in Theorem 5.11, thus we have:

[TABLE]

5.3 Error estimates for LPS-ROM

We are now in position to prove the following error estimate result for the LPS-ROM defined by (4.11):

Theorem 5.14**.**

Under the regularity assumption on the continuous solution (Hypothesis 5.5), the assumption on the stabilization parameters (Hypothesis 5.6), and supposing that $\|{{\boldsymbol{u}}^{0}-{\boldsymbol{u}}_{r}^{0}}\|_{{\bf L}^{2}}=\mathcal{O}(h^{s})$, the solution of the LPS-ROM (4.11) satisfies the following error estimate:

[TABLE]

*where $C^{*}$ is a positive constant that will be determined throughout the proof. *

Proof. We start deriving the error bounds by splitting the error for the velocity and the pressure into two terms:

[TABLE]

[TABLE]

In (5.27), the first term, $\boldsymbol{\eta}^{n+1}={\boldsymbol{u}}^{n+1}-{\boldsymbol{v}}_{r}^{n+1}$, represents the difference between ${\boldsymbol{u}}^{n+1}$ and its $L^{2}$-orthogonal projection on ${\boldsymbol{X}}_{r}$ (5.7).The second term, $\boldsymbol{\phi}_{r}^{n+1}={\boldsymbol{u}}_{r}^{n+1}-{\boldsymbol{v}}_{r}^{n+1}$, is the remainder. Similarly, in (5.28), the first term, $\rho^{n+1}=p^{n+1}-z_{r}^{n+1}$, represents the difference between $p^{n+1}$ and its $L^{2}$-orthogonal projection on $Q_{r}$ (5.8).The second term, $s_{r}^{n+1}=p_{r}^{n+1}-z_{r}^{n+1}$, is the remainder.

Next, we construct the error equation. We first evaluate the weak formulation of the continuous NSE (2.2) at $t=t_{n+1}$ and let ${\boldsymbol{v}}=\boldsymbol{\varphi}$, then subtract the LPS-ROM (4.11) from it. For any $\boldsymbol{\varphi}\in{\boldsymbol{X}}_{r}$, we obtain:

[TABLE]

By adding and subtracting the difference quotient term $\displaystyle\frac{1}{\Delta t}({\boldsymbol{u}}^{n+1}-{\boldsymbol{u}}^{n},\boldsymbol{\varphi})$ in (5.29), and applying the decompositions (5.27)-(5.28), we get, for any $(\boldsymbol{\varphi},\psi)\in{\boldsymbol{X}}_{r}\times Q_{r}$:

[TABLE]

Note that $\displaystyle\frac{1}{\Delta t}(\boldsymbol{\eta}^{n+1}-\boldsymbol{\eta}^{n},\boldsymbol{\varphi})=0$, since ${\boldsymbol{v}}_{r}^{n+1}$ is the $L^{2}$-orthogonal projection of ${\boldsymbol{u}}^{n+1}$ on ${\boldsymbol{X}}_{r}$, and $\sigma(\rho^{n+1},\psi)=0$, since $z_{r}^{n+1}$ is the $L^{2}$-orthogonal projection of $p^{n+1}$ on $Q_{r}$. Choosing $\boldsymbol{\varphi}=2\Delta t\boldsymbol{\phi}_{r}^{n+1}$, $\psi=2\Delta ts_{r}^{n+1}$ in (5.30) and letting $\boldsymbol{c}^{n}=\partial_{t}{\boldsymbol{u}}^{n+1}-\displaystyle\frac{{\boldsymbol{u}}^{n+1}-{\boldsymbol{u}}^{n}}{\Delta t}$, we get:

[TABLE]

Using the polarization identity:

[TABLE]

from (5.31) we obtain:

[TABLE]

We estimate the terms on the r.h.s. of (5.32) one by one. By definition of the dual norm and Young’s inequality, we get for the first term on the r.h.s. of (5.32):

[TABLE]

for some small positive constant $\varepsilon_{1}$ (to be determined later).

The nonlinear convective terms in (5.32) can be written as follows:

[TABLE]

where we have used $b({\boldsymbol{u}}_{r}^{n+1},\boldsymbol{\phi}_{r}^{n+1},\boldsymbol{\phi}_{r}^{n+1})=0$, which follows from (5.1). Next, we estimate each term on the r.h.s. of (5.34). Since ${\boldsymbol{u}}_{r}^{n+1},\boldsymbol{\eta}^{n+1},\boldsymbol{\phi}_{r}^{n+1},{\boldsymbol{u}}^{n+1}\in{\boldsymbol{X}}$, we can apply the standard bound (5.2) for the trilinear form $b(\cdot,\cdot,\cdot)$ and use Young’s inequality to get:

[TABLE]

[TABLE]

For the last nonlinear convective term, applying Hölder’s inequality, Sobolev embedding theorem, and Young’s inequality yields:

[TABLE]

By Cauchy–Schwarz and Young’s inequalities, we bound the fourth and fifth terms on the r.h.s. of (5.32):

[TABLE]

For the sixth term on the r.h.s of (5.32), applying triangle, Cauchy–Schwarz and Young’s inequalities plus Lemma 5.3, we have:

[TABLE]

for some small positive constant $\varepsilon_{2}$ (to be determined later). In particular, we have applied Lemma 5.3 with $Y=\text{span}\{\nabla\cdot\boldsymbol{\varphi}_{1},\ldots,\nabla\cdot\boldsymbol{\varphi}_{r}\}$ and $Z=Q_{r}$.

Finally, for the last term on the r.h.s. of (5.32) we have:

[TABLE]

where we have used Cauchy–Schwarz and Young’s inequalities.

Substituting inequalities (5.33)-(5.41) in (5.32) and taking $\varepsilon_{1}=\nu/12$, $\varepsilon_{2}=\sigma/8$, we obtain:

[TABLE]

where in the last inequality we have used the identity $z_{r}^{n+1}=p^{n+1}-\rho^{n+1}$ from (5.28), triangle inequality, and the inequality $\|{\nabla\cdot{\boldsymbol{u}}^{n+1}}\|_{L^{2}}^{2}\leq\sigma^{2}\|{p^{n+1}}\|_{L^{2}}^{2}$ from (2.2).

Summing (5.42) from $n=0$ to $k\leq N-1$, we have:

[TABLE]

Next, we estimate each term on the r.h.s. of (5.43).

The first term on the r.h.s. of (5.43) can be estimated as follows:

[TABLE]

where the last inequality follows from the fact that ${\boldsymbol{v}}_{r}^{0}$ is the $L^{2}$-orthogonal projection of ${\boldsymbol{u}}^{0}$ on ${\boldsymbol{X}}_{r}\subset{\boldsymbol{X}}_{h}$, so that it satisfies optimal approximation properties similar to standard FE interpolations (cf. [17]), and we have supposed $\|{{\boldsymbol{u}}^{0}-{\boldsymbol{u}}_{r}^{0}}\|=\mathcal{O}(h^{s})$.

By Taylor’s theorem, the second term on the r.h.s. of (5.43) can be estimated as follows:

[TABLE]

where we have used the regularity assumption 5.5 on the continuous velocity.

To estimate the third term on the r.h.s. of (5.43), we use Theorem 5.11 and the fact that ${\boldsymbol{v}}_{r}^{n+1}$ is the $L^{2}$-orthogonal projection of ${\boldsymbol{u}}^{n+1}$ on ${\boldsymbol{X}}_{r}\subset{\boldsymbol{X}}_{h}$, so that it satisfies optimal approximation properties as standard FE interpolations (cf. [17]):

[TABLE]

By using the regularity assumption 5.5 on the continuous velocity and $H^{1}$ velocity projection error estimate (5.11), the fourth term on the r.h.s. of (5.43) can be estimated as follows:

[TABLE]

Similarly, we have:

[TABLE]

and, by using $L^{2}$ pressure projection error estimate (5.10):

[TABLE]

The ninth term on the r.h.s. of (5.43) makes apparent the convergence order reduction linked to the diffusive nature of the penalty term. Indeed, by using the regularity assumption 5.5 on the continuous pressure, we get:

[TABLE]

For the last two terms on the r.h.s. of (5.43), we obtain:

[TABLE]

and:

[TABLE]

where we have used the assumption on stabilization parameter 5.6, the regularity assumption 5.5 on the continuous pressure, stability and optimal error estimates of $\Pi_{h}$ (cf. [17]), and $H^{1}$ pressure projection error estimate (5.12).

Collecting (5.44)-(5.52), estimate (5.43) becomes:

[TABLE]

The discrete Grönwall’s lemma (cf. Lemma 27 in [30]) implies the following inequality:

[TABLE]

where $C^{*}=Ce^{C_{1}\Delta t\sum_{n=0}^{N-1}\|{\nabla{\boldsymbol{u}}^{n+1}}\|_{{\bf H}^{1}}^{2}}$, and $C_{1}$ is a positive constant depending on $\nu^{-1}$. Finally, using in (5.54) the inequality:

[TABLE]

triangle inequality, $L^{2}-H^{1}$ projection error estimates for velocity (5.9)-(5.11) and $L^{2}$ projection error estimate for pressure (5.10), we get:

[TABLE]

This concludes the proof. $\Box$

Remark 5.15**.**

Note that in error estimate (5.26) the convergence order is limited by the penalty constant $\sigma$. To remove this order limitation, one could simply replace bound (5.40) in Theorem 5.14 by:

[TABLE]

applying Cauchy–Schwarz and Young’s inequalities. Taking $\varepsilon_{2}=\sigma/2$, this leads to the error estimate:

[TABLE]

So, we have removed the order limitation due to the penalty constant $\sigma$, but changed the constant $(1/\nu+\alpha^{2}/\sigma)$ by $(1/\nu+1/\sigma)$. In any case, note that to prove error estimates, we have to assume to work with sufficiently regular flows, which is a common approach in the derivation of error estimates for POD-ROM (cf. [25, 32]). Thus, the convergence orders obtained are generally valid in laminar flow settings or for sufficiently regular flows, but are usually not valid in realistic turbulent flow settings, since the convergence order decreases with the regularity. In this context, the kinematic viscosity is usually of the order $\nu=\mathcal{O}(10^{-3})$. Also, we have observed in the numerical studies performed in Section 6 that the saturation constant $\alpha$ starts with small values for small $r$ and seems to experience a flattening effect around $10^{-2}$ when adding more POD modes. Thus, for practical relatively coarse FE mesh size $h$ and time step $\Delta t$, we have that estimate (5.26), even if contains explicitly the penalty constant $\sigma$, improves estimate (5.57), since the presence of the saturation constant $\alpha$ allows to ease the convergence order reduction due to $\sigma^{-1}$.

Remark 5.16**.**

Following Remark 5.13, an alternative unconditional pressure error estimate for LPS-ROM (4.11) could be obtained if we set $\psi=0$ in the error equation (5.30). Thus, we obtain for any $\boldsymbol{\varphi}\in{\boldsymbol{X}}_{r}$:

[TABLE]

with $\boldsymbol{\varepsilon}_{r}^{n+1}$ denoting the consistency error, defined as:

[TABLE]

Let $S_{r}^{n+1}=\displaystyle\sum_{k=0}^{n}\Delta ts_{r}^{k+1}=\displaystyle\sum_{k=0}^{n}\Delta t(p_{r}^{k+1}-z_{r}^{k+1})$, then summation over the discrete times gives:

[TABLE]

Thus, we get:

[TABLE]

where we have applied triangle inequality, Cauchy–Schwarz inequality, the definition of the dual norm, and the standard estimate (5.2) for the convective term. Then, using the norm defined in (5.20), by Cauchy–Schwarz inequality, the stability result (5.15) for the reduced order velocity and the regularity assumption 5.5 on the continuous velocity, from (5.58) we have:

[TABLE]

Using estimates (5.45)-(5.48) to bound the third term on (5.59), and estimate (5.54) to bound the rest of terms in (5.59), we obtain:

[TABLE]

where $E$ denotes the sum of all terms within brackets in (5.26). From triangle inequality and the $L^{2}$ projection error estimates for pressure (5.10), it follows:

[TABLE]

and we have denoted $P(\cdot,t)=\displaystyle\int_{0}^{t}\widetilde{p}(\cdot,s)\,ds$, being $\widetilde{p}$ the piecewise constant in time function that takes the value $p^{n+1}$ on $(t_{n},t_{n+1})$. Thus, we have derived a new error estimate on the time-average of the reduced order pressure with respect to the one obtained in Theorem 5.14:

[TABLE]

6 Numerical studies

In this section, we present numerical results for the LPS-ROM introduced and analyzed in the previous sections, for which the standard discrete inf-sup condition is circumvented and neither strongly nor weakly divergence-free POD modes are required. The numerical experiments are performed on the benchmark problem of the 2D laminar unsteady flow around a cylinder with circular cross-section [37]. The open-source FE software FreeFEM [23] has been used to run the numerical experiments.

6.1 Setup for numerical simulations

Following [37], the computational domain is given by a rectangular channel with a circular hole (see Figure 1 for the computational grid used):

[TABLE]

No slip boundary conditions are prescribed on the horizontal walls and on the cylinder, and a parabolic inflow profile is provided at the inlet:

[TABLE]

with $U_{m}={\boldsymbol{u}}(0,H/2,t)=1.5\,\rm{m/s}$, and $H=0.41\,\rm{m}$ the channel height. At the outlet, we perform a comparison using on one side outflow (do nothing) boundary conditions $(\nu\nabla{\boldsymbol{u}}-p\,Id){\boldsymbol{n}}={\bf 0}$, with ${\boldsymbol{n}}$ the outward normal to the domain, for which we can remove the penalty term with factor $\sigma$ to the variational formulation, since the constant the pressure is determined up through the formulation is already fixed by the prescribed boundary conditions. On the other side, we impose the same parabolic inflow profile for the outflow velocity, for which we use the penalty term with factor $\sigma$ to fix the constant the pressure is determined up through the formulation. Although imposing Dirichlet boundary conditions at the outlet is unphysical, we have seen that it simplifies the theoretical analysis performed. Also, we have observed in the numerical studies that this do not influence too much quantities of interest such as the upstream drag and lift coefficients, and it is often used [26, 35], also in the ROM framework [33, 41].

The kinematic viscosity of the fluid is $\nu=10^{-3}\,\rm{m^{2}/s}$, and there is no external (gravity) forcing, i.e. ${\boldsymbol{f}}={\bf 0}\,\rm{m/s^{2}}$. Based on the mean inflow velocity $\overline{U}=2U_{m}/3=1\,\rm{m/s}$, the cylinder diameter $D=0.1\,\rm{m}$ and the kinematic viscosity of the fluid $\nu=10^{-3}\,\rm{m^{2}/s}$, the Reynolds number considered is $Re=\overline{U}D/\nu=100$. In the fully developed periodic regime, a vortex shedding can be observed behind the obstacle, resulting in the well-known von Kármán vortex street (see Figures 2-3).

For the evaluation of computational results, we are interested in studying the temporal evolution of the following quantities of interest. The kinetic energy of the flow is the most frequently monitored quantity, given by:

[TABLE]

Other relevant quantities of interest are the drag and lift coefficients. In order to reduce the boundary approximation influences, in the present work these quantities are computed as volume integrals:

[TABLE]

[TABLE]

for arbitrary test functions ${\boldsymbol{v}}_{D},{\boldsymbol{v}}_{L}\in{\bf H}^{1}$ such that ${\boldsymbol{v}}_{D}=(1,0)^{T}$ on the boundary of the cylinder and vanishes on the other boundaries, ${\boldsymbol{v}}_{L}=(0,1)^{T}$ on the boundary of the cylinder and vanishes on the other boundaries. In the actual computations, we have used the approach described in [27] to fix the test functions ${\boldsymbol{v}}_{D},{\boldsymbol{v}}_{L}$ and evaluate the drag and lift coefficients $c_{D},c_{L}$. Reference intervals for these coefficients were given in [37] (see last row of Table 1), together with the Strouhal number $St=Df/\overline{U}$, where $f$ is the frequency of the vortex shedding.

6.2 LPS-FOM for snapshots generation

The numerical method used to compute the snapshots is the LPS-FOM described in Section 3, with a spatial discretization using EO ${\bf P}^{2}-\mathbb{P}^{2}$ FE for the pair velocity-pressure on a relatively coarse computational grid (see Figure 1, $h=2.76\times 10^{-2}\,\rm{m}$), resulting in $32\,488$ d.o.f. for velocities and $16\,244$ d.o.f. for pressure. We perform a comparison using do nothing Boundary Conditions (BC) at the outlet, for which $\sigma=0$, and Dirichlet BC at the outlet, for which $\sigma=10^{-6}$. The following expression of the stabilization parameters is used in the computations:

[TABLE]

by adapting the form proposed in [18, 19], designed by a specific Fourier analysis applied in the framework of stabilized methods. For the time discretization, a semi-implicit Backward Differentiation Formula of order 2 (BDF2) has been applied, which guarantees a good balance between numerical accuracy and computational complexity (cf. [3]). In particular, the discrete time derivative has been approximated by the operator $D_{t}^{2}$ defined as:

[TABLE]

and we have considered the following extrapolation for the convection velocity by means of Newton–Gregory backward polynomials (cf. [12]): $\widehat{{\boldsymbol{u}}}_{h}^{n}=2{\boldsymbol{u}}_{h}^{n}-{\boldsymbol{u}}_{h}^{n-1}$, $n\geq 1$, in order to achieve a second-order accuracy in time. For the initialization $(n=0)$, we have considered ${\boldsymbol{u}}_{h}^{-1}={\boldsymbol{u}}_{h}^{0}={\boldsymbol{u}}_{0h}$, being ${\boldsymbol{u}}_{0h}$ the initial condition, so that the time scheme reduces to the semi-implicit Euler method for the first time step $(\Delta t)^{0}=(2/3)\Delta t$. In the FOM simulations, an impulsive start is performed, i.e. the initial condition is a zero velocity field, and the time step is $\Delta t=2\times 10^{-3}\,\rm{s}$. Time integration is performed till a final time $T=7\,\rm{s}$. In the time period $[0,5]\,\rm s$, after an initial spin-up, the flow is expected to develop to full extent, including a subsequent relaxation time. Afterwards, it reaches a periodic-in-time (statistically- or quasi-steady) state (see Figure 4, where we plot in particular drag and lift coefficients, and kinetic energy temporal evolution for the FOM solutions).

Observe from Figures 2-3 that results are qualitatively very close using do-nothing and Dirichlet BC at the outlet. The only very little difference can be noticed in the vorticity plot, where some numerical oscillations of the vorticity field appear when using Dirichlet BC at the outlet. Nevertheless, they have low magnitude and just slightly influence drag and lift coefficients reported in Table 1. Indeed, results with do nothing BC for all quantities of interest agree quite well with reference results from [37] and other numerical studies [11, 33, 41], and are only slightly less accurate when using Dirichlet BC at the outlet. This behavior is reflected also in Figure 4, where in addition we have also plotted the “weak-strong” divergence temporal evolution for the FOM solution. In parrticular, the curve of the “strong” divergence has been obtained by plotting $\|{\nabla\cdot{\boldsymbol{u}}_{h}^{n}}\|_{L^{2}}$, while the curve of the “weak” divergence has been obtained by plotting $\max_{q_{h}\in Q_{h}}|(\nabla\cdot{\boldsymbol{u}}_{h}^{n},q_{h})|$. From this figure, it is evident that the computed snapshots are neither strongly nor weakly divergence-free.

The POD modes are generated in $L^{2}$ by the method of snapshots with velocity centered-trajectories [22] by storing every fifth FOM solution (with do nothing and Dirichlet BC at the outlet) in the stable response time interval $[5,7]\,\rm{s}$, so that $200$ snapshots were used both for velocity and pressure. Figures 5-6 display the Euclidean norm of the first POD velocity modes, obtained using do nothing BC and Dirichlet BC at the outlet, respectively. Similarly, Figures 7-8 display the first POD pressure modes, obtained using do nothing BC and Dirichlet BC at the outlet, respectively.

Only slight differences can be noticed between the POD modes using do nothing BC and Dirichlet BC at the outlet, being the most noticeable ones given by the numerical noise of the second POD mode using Dirichlet BC at the outlet. In Figure 9, we show the decay of POD velocity ($\lambda_{i}$) and pressure ($\gamma_{i}$) eigenvalues (top), together with the corresponding captured system’s energy (bottom), computed respectively as $100\sum_{i=1}^{r}\lambda_{i}/\sum_{i=1}^{M_{v}}\lambda_{i}$ and $100\sum_{i=1}^{r}\gamma_{i}/\sum_{i=1}^{M_{p}}\gamma_{i}$.

Note that the first $r=5$ POD modes already capture more than $99\%$ of the system’s velocity-pressure energy.

We recall that the new LPS-ROM uses the same number $r$ of velocity and pressure POD modes. Thus, we expect that also the POD velocity-pressure spaces do not satisfy the standard discrete inf-sup condition and the LPS-ROM try to circumvent it trough numerical stabilization. Effectively, we have checked numerically that, varying $r$, the discrete inf-sup constant $\beta$ for the POD velocity-pressure spaces remains very close to zero (below $10^{-7}$), as we can observe from Figure 10 (top). In Figure 10 (bottom), we also display the saturation constant $\alpha$ (see Lemma 5.3) between the spaces $Y=\text{span}\{\nabla\cdot\boldsymbol{\varphi}_{i},\ldots,\nabla\cdot\boldsymbol{\varphi}_{r}\}$ and $Z=Q_{r}$ in terms of $r$:

We can observe that the saturation constant $\alpha$, for the chosen numerical setup, starts with small values for small $r$ (around $10^{-5}$, $10^{-6}$) and experiences a flattening effect with values within $(10^{-2},10^{-1})$ when adding more POD modes ($r>10$). This could explain why we do not experience in practice so much improvement when adding more POD modes and we compare the ROM and the FOM solutions. In terms of the theoretical analysis performed, we have that $\alpha^{2}<\sigma$ for small $r$, which allows to significantly ease the convergence order reduction due to $\sigma^{-1}$, but then $\alpha^{2}$ becomes predominant and experiences a flattening effect, which implies no much improvement in the convergence order when adding already more than $7$ POD modes. This could give a sort of criterion (or at least an idea) on how many POD modes are needed to reach, using the new LPS-ROM, a reasonable accuracy with respect to the FOM solution, and thus well catch physical quantities of interest at a very reduced computational cost.

6.3 Numerical results for LPS-ROM

With POD velocity-pressure modes generated, the fully discrete LPS-ROM is constructed as discussed in Section 4 using the semi-implicit BDF2 time scheme as for the LPS-FOM, and run in the stable response time interval $[5,7]\,\rm{s}$ with $\Delta t=2\times 10^{-3}\,\rm{s}$.

To assess the behavior of the new LPS-ROM, the temporal evolution of the drag and lift coefficients, and kinetic energy are monitored and compared to the FOM solutions in the stable response time interval $[5,7]\,\rm{s}$, corresponding to six periods for the lift coefficient. The corresponding numerical results with do nothing and Dirichlet BC at the outlet are reported in Figures 11-12. The LPS-ROM has been computed with $r=3,5,7$ POD modes for velocity and pressure.

We observe that the new LPS-ROM is able to replicate the temporal evolution of the lift coefficient computed with the LPS-FOM already with $r=3$. As for the drag coefficient, $r=5$ POD modes are already sufficient to capture its corresponding FOM values. More sensitive seems to be the kinetic energy, for which however $r=7$ POD modes are already able to catch up the FOM solution with a reasonable accuracy. Apart from the temporal evolution of the drag coefficient with $r=3$, results with do nothing BC and Dirichlet BC at the outlet are almost similar for the new LPS-ROM. Qualitatively, it is interesting to observe that, as the LPS-FOM, the new LPS-ROM shows a fully periodic time evolution for all monitored quantities of interest at all levels ($r=3,5,7$).

To better assess on the one hand the behavior of the new LPS-ROM and partially illustrate on the other hand the theoretical convergence order predicted by the numerical analysis and stated in Theorem 5.14, we plot relative errors with respect to the LPS-FOM solution. In particular, in Figure 13, we first plot the temporal evolution of the discrete $L^{2}$ relative error (in semilogarithmic scale) of the reduced order velocity and pressure with respect to the full order ones: $\|{{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{r}}\|_{{\bf L}^{2}}/\|{{\boldsymbol{u}}_{h}}\|_{{\bf L}^{2}}$ and $\|{p_{h}-p_{r}}\|_{L^{2}}/\|{p_{h}}\|_{L^{2}}$, for $r=3,5,7$, obtained with do nothing BC. Similarly, Figure 14 shows the temporal evolution of the discrete $L^{2}$ relative error of the reduced order velocity and pressure with respect to the full order ones for $r=3,5,7$, obtained with Dirichlet BC at the outlet.

As expected, the errors decrease when increasing $r$, but then they almost stabilize for $r\geq 7$. It can be observed that the use of do nothing and Dirichlet BC at the outlet provides almost similar results. However, in order to limit the influence of the spatial error due to the FE discretization and the temporal error due to the time-stepping scheme, we have decided to test the $\ell^{2}(L^{2})$ relative errors of the reduced velocity and pressure with respect to the $L^{2}$ projection of the full order ones on the respective POD spaces, and plot them in terms of $\Lambda_{r}=\sum_{i=r+1}^{M_{v}}\lambda_{i}$ and $Z_{r}=\sum_{i=r+1}^{M_{p}}\gamma_{i}$, respectively. In this way, the analysis is focused on the ROM error due to the POD truncation and one would expect a linear regression behavior (in logarithmic scale), as suggested by the performed numerical analysis. Effectively, this behavior is numerically recovered in Figure 15 for both reduced velocity and pressure.

7 Summary and conclusions

In this paper, we have proposed a new stabilized projection-based reduced order method (LPS-ROM) for the numerical simulation of incompressible flows. In particular, the new LPS-ROM is a velocity-pressure ROM that uses pressure modes as well to compute the reduced order pressure, needed for instance in the computation of relevant quantities, such as drag and lift forces on bodies in the flow. With respect to other approaches existing in the current ROM literature that provides velocity-pressure approximations, the new LPS-ROM circumvents the standard discrete inf-sup condition for the POD velocity-pressure spaces, whose fulfillment can be rather expensive and inefficient in realistic applications in CFD, see for instance [7, 39], where an offline strategy based on the supremizer enrichment of the reduced velocity space has been proposed and applied in the POD context, adapted from the RB method framework. Also, the velocity modes for the new LPS-ROM does not have to be neither strongly nor weakly divergence-free, which allows to use snapshots generated for instance with penalty or projection-based stabilized methods. This is not the case, for instance, of ROM based on a pressure Poisson equation approach (see, for instance, the first two methods investigated in [11] and also [39]), for which the velocity snapshots, and hence the POD velocity modes must be at least weakly divergence-free. This requirement also holds for the last method proposed and investigated in [11], which uses a residual-based stabilization mechanism in order to overcome a possible violation of the discrete inf-sup condition in the ROM framework by considering a decoupled approach for the reduced velocity-pressure pair.

The main contribution of the present paper has been to perform a stability and convergence analysis of the arising fully discrete LPS-ROM applied to the unsteady incompressible NSE, by mainly deriving the proof of a rigorous error estimate that considers all contributions: the spatial discretization error (due to the FE discretization), the temporal discretization error (due to the backward Euler method), and the POD truncation error. In particular, the numerical analysis corroborated with numerical studies makes apparent an interesting link between the number of POD velocity-pressure modes used in the ROM and the angle $\theta$ between the space spanned by the divergence of the POD velocity modes and the POD pressure space. Indeed, for the numerical setting proposed, where the same numerical stabilization technique used for the ROM is initially applied also to the FOM (LPS-FOM) to generate the snapshots, so that these latter are not weakly divergence-free, we have found that for small values of $r$, which is common in practice, the saturation constant $\alpha=cos(\theta)$ is rather small, and this allows to ease the convergence order reduction due to the violation of the discrete inf-sup stability condition.

Numerical studies performed on a two-dimensional laminar unsteady flow past a circular obstacle have also been used to assess the accuracy and efficiency of the new LPS-ROM. Despite the fact that the discrete inf-sup condition is not fulfilled by the new LPS-ROM, using a small equal number of POD velocity-pressure modes already provides accurate approximations, close to the LPS-FOM results, and theoretical scalings suggested by the numerical analysis are recovered in practice.

We plan to extend the theoretical work of this paper to the numerical analysis of the last method proposed and investigated in [11] but for not exactly weakly divergence-free POD velocity modes, by also performing a numerical investigation that both supports the analytical results and illustrates the comparison of that method with the new LPS-ROM here proposed as stabilization-motivated ROM. This theoretical and computational study is today in progress and shall appear in a forthcoming paper.

Acknowledgments: This work has been partially supported by the Spanish Government-EU Feder grant MTM2015-64577-C2-1-R. The research of the author has been also funded by the Spanish State Research Agency through the national programme Juan de la Cierva-Incorporación 2017. The author acknowledges Prof. M. Azaïez (Laboratoire I2M CNRS UMR5295, Institut Polytechnique de Bordeaux) for some advices and fruitful suggestions, and Prof. T. Chacón (Departamento EDAN & IMUS, Universidad de Sevilla) for some helpful discussions, especially on Remark 5.13.