Stabilized POD reduced order models for convection-dominated incompressible flows

Abstract

We present a comparative computational study of two stabilized Reduced Order Models (ROMs) for the simulation of convection-dominated incompressible flow (Reynolds number of the order of a few thousands). Representative solutions in the parameter space, which includes either time only or time and Reynolds number, are computed with a Finite Volume method and used to generate a reduced basis via Proper Orthogonal Decomposition (POD). Galerkin projection of the Navier–Stokes equations onto the reduced space is used to compute the ROM solution. To ensure computational efficiency, the number of POD modes is truncated and ROM solution accuracy is recovered through two stabilization methods: (i) adding a global constant artificial viscosity to the reduced dimensional model, and (ii) adding a different value of artificial viscosity for the different POD modes. We test the stabilized ROMs for fluid flow in an idealized medical device consisting of a conical convergent, a narrow throat, and a sudden expansion. Both stabilization methods significantly improve the ROM solution accuracy over a standard (non-stabilized) POD-Galerkin model.

1 Introduction

For a few decades, Reduced Order Models (ROMs) have become the methodology of choice to reduce the computational cost when traditional computational techniques, e.g., Finite Element methods and Finite Volume methods, have to be carried out for several parameter values, as is the case in uncertainty quantification, optimal control, and inverse problems. See, e.g., (Peter et al. 2021; Benner et al. 2020a, b; Hesthaven et al. 2016; Rozza et al. 2008) for reviews on reduced order modeling. ROMs rely on traditional computational methods, also called Full Order Models (FOMs), to build a lower-dimensional approximation of the solution space that captures the essential behavior of the solutions in a given parameter space. In this paper, the parameters will be time and the Reynolds number, which is a parameter related to the physics of the problem. The construction of the lower-dimensional approximation is done in a preliminary phase, called offline, during which one generates a database of several FOM solutions associated to given times and physical parameter values. The database of FOM solutions is used to generate a reduced basis, which is (hopefully much) smaller than the high-dimensional FOM basis but still preserves the essential features of the system. One of the most used techniques to generate the reduced bases space is Proper Orthogonal Decomposition (POD), which extracts the dominant modes from the FOM database. In a second phase, called online, one uses this reduced basis to quickly compute the solution for newly specified times and parameter values. Note that, while the offline phase is performed only once, the online phase is performed as many times as needed.

ROMs are very efficient surrogate models when the number of reduced basis functions is small (i.e., ${\mathcal {O}}$(10)), as is typically the case for diffusion-dominated flows. Unfortunately, though, a large number of basis functions is needed to capture the essential features of convection-dominated flows. If one retains such a large number of modes for the ROM to be accurate, then the computational efficiency suffers. If the number of modes is otherwise kept low, a severe loss of information hinders the accurate reconstruction of the solution. In fact, projection-based ROMs of turbulent flows are affected by energy stability problems related to the fact that POD retains the modes biased toward large, high-energy scales, while the turbulent kinetic energy is dissipated at the level of the small scales. A possible way to tackle this challenging problem is to introduce dissipation via a closure model (Wang et al. 2012; Aubry et al. 1988; Ahmed et al. 2021) or stabilization (Carlberg et al. 2011; Grimberg et al. 2020; Baiges et al. 2013; Reyes and Codina 2020; Parish et al. 2020).

We consider two stabilization techniques that introduce artificial viscosity in the online phase of a POD-Galerkin ROM of the incompressible Navier–Stokes equations. By POD-Galerkin ROM, it is meant that that projection of the system equations onto the reduced space is used to compute the ROM solution during the online phase. The first stabilization method we consider is Heisenberg stabilization (Bergmann et al. 2009; San and Iliescu 2015), where a global constant artificial viscosity is added in the reduced dimensional model. The second method is a variation of the first, where a different value of artificial viscosity is used for the different modes (Cazemier 1997; San and Iliescu 2015). We call this method coefficient-dependent viscosity where coefficient refers to the POD mode index. Specifically, a small amount of dissipation is introduced to the lowest index mode and it increases proportionally to the mode index for every other POD mode. So far, these approaches have been tested on academic problems: 2D square cylinder wake flow in laminar regime, 2D lid-driven cavity flow in turbulent regime and mid-latitude simplified oceanic flow (a 2D problem too).

We test the proposed stabilized POD-Galerkin ROMs on a realistic 3D problem: flow at different Reynolds numbers through a nozzle that contains all the features commonly encountered in medical devices, such as prosthetic heart valves, stents, and heart pumps (flow contraction and expansion, recirculation zones, etc.). See Fig. 1. This benchmark test was originally proposed by the U.S Food and Drug Administration (FDA), which asked three independent laboratories to perform flow visualization experiments on fabricated nozzles for different flow rates up to Reynolds number 6500 (Hariharan et al. 2011). This led to the creation of measurement datasets that can be used to validate Computational Fluid Dynamics (CFD) simulations. The focus of the FDA benchmark is blood flow simulations in larger arteries (hence the cap of the Reynolds number to 6500) to design and assess biomedical devices. However, the numerical challenges this kind of simulations poses are broadly applicable to other fields, e.g., aeronautics or environmental flows. As FOM, we use a Finite Volume method discussed in Girfoglio et al. (2019) and implemented in (Weller et al. 1998; Jasak 1996). The stabilized POD-Galerkin ROMs are implemented in ITHACA-FV (Stabile et al. 2017; Stabile and Rozza 2018), an open-source finite volume C++ library,¹ and compared for Reynolds numbers 2000 (transitional regime), 3500, 5000, and 6500 (progressively more developed turbulent regimes).

The rest of the paper is organized as follows. Section 2 states the Navier–Stokes equations and briefly describes the Finite Volume method used as FOM. Section 3 presents the proposed stabilized POD-Galerkin ROMs. Section 4 presents numerical results obtained for the FDA benchmark. Finally, Sect. 5 summarises the conclusions.

2 Full order model

2.1 Problem formulation

The flow of an incompressible, viscous, and Newtonian fluid in domain $\Omega$ and over time interval $(t_0, T)$ is modeled with the parametrized time-dependent incompressible Navier–Stokes equations:

where $\varvec{x}$ is the space variable, t is time, $\varvec{u}(\varvec{x}, t, \varvec{\mu })$ is the velocity vector, $p(\varvec{x}, t, \varvec{\mu })$ is the pressure, $\nu$ is the kinematic viscosity, and $\varvec{\mu }$ is the parameter vector, i.e., a vector containing all the physical parameters, other than time and viscosity, the problem depends upon. In the rest of the paper, for a lighter notation, we will omit the spatial-temporal dependence of the variables and the dependence on the parameters $\varvec{\mu }$ as well. To characterize the flow regime, we introduce the Reynolds number, which is defined as:

$$\begin{aligned} Re = \frac{U L}{\nu }, \end{aligned}$$

(3)

where U and L are characteristic macroscopic speed and length for the flow. For an internal flow in a cylindrical pipe, U is the mean sectional velocity and L is the pipe diameter. For large Reynolds numbers, inertial forces are dominant over viscous forces and vice versa for small Reynolds numbers.

Let $\Gamma _i$ be the inlet boundary, $\Gamma _w$ the wall, and $\Gamma _o$ the outlet, such that $\Gamma _i \cup \Gamma _w \cup \Gamma _o = \partial \Omega$ and $\Gamma _i \cap \Gamma _w \cap \Gamma _o = \emptyset$. We impose a non-homogeneous Dirichlet boundary conditions at the inlet, no-slip conditions on the wall and a homogeneous Neumann condition at the outlet:

$$\begin{aligned} {\left\{ \begin{array}{ll} \varvec{u} = \varvec{u}_D{(\varvec{x})} & \quad \hbox { in}\ \Gamma _i\times (t_0, T), \\ \varvec{u} = \varvec{0} & \quad \hbox { in}\ \Gamma _w\times (t_0, T), \\ (2\mu \nabla \varvec{u} - p I)\varvec{n} = \varvec{0} & \quad \text {in}\ \Gamma _o\times (t_0, T). \end{array}\right. } \end{aligned}$$

(4)

In (4), $\varvec{u}_D{(\varvec{x})}$ is a given inlet velocity profile and $\varvec{n}$ is the unit normal outward vector to the boundary. For the results in Sect. 4, $\varvec{u}_D = \varvec{u}_D (\varvec{x})$, i.e., the velocity profile is space dependent.

2.1.1 Finite volume approximation

To discretize problem (1)-(2) in space, we adopt a Finite Volume (FV) method. For this purpose, we divide domain $\Omega$ into control volumes $\Omega _i$, for $i = 1,\dots ,N_h$, such that $\bigcup _{i=1}^{N_h} \Omega _i=\Omega$ and $\Omega _i\cap \Omega _j=\emptyset$ for $i,j = 1,\dots ,N_h$ and $i\ne j$.

By integrating equation (1) over control volume $\Omega _i$ and applying the Gauss-divergence theorem, we obtain:

$$\begin{aligned} \int _{\Omega _i} \frac{\partial \varvec{u}}{\partial t} d\Omega + \int _{\partial \Omega _i} (\varvec{u} \otimes \varvec{u}) \cdot d\varvec{A} - \nu \int _{\partial \Omega _i} \nabla \varvec{u} \cdot d\varvec{A} + \int _{\partial \Omega _i} p d\varvec{A} = 0, \end{aligned}$$

(5)

where $\partial \Omega _i$ is the boundary of the control volume $\Omega _i$, $d\Omega$ is an infinitesimal volume element, and $d\varvec{A}$ is an infinitesimal surface element. Boundary terms are approximated with the sum on each face j of cell $\Omega _i$ as follows:

$$\begin{aligned}&\int _{\partial \Omega _i} p d\varvec{A} \approx \sum _{j} p_{i,j} \varvec{A}_{ij}, \quad&\text{(Gradient } \text{ term) } \end{aligned}$$

(6)

$$\begin{aligned}&\int _{\partial \Omega _i} (\varvec{u} \otimes \varvec{u}) \cdot d\varvec{A} \approx \sum _{j} (\varvec{u}_{i,j} \otimes \varvec{u}_{i,j}) \cdot \varvec{A}_{ij} = \sum _{j} (\varvec{u}_{i,j} \cdot \varvec{A}_{ij}) \varvec{u}_{i,j}, \quad&\text{(Convective } \text{ term) } \end{aligned}$$

(7)

$$\begin{aligned}&\int _{\partial \Omega _i} \varvec{\nabla u} \cdot d\varvec{A} \approx \sum _{j} (\varvec{\nabla u_i})_{j} \cdot \varvec{A}_{ij}, \quad&\text{(Diffusion } \text{ term) } \end{aligned}$$

(8)

where $p_{i,j}$, $\varvec{u}_{i,j}$, and $(\varvec{\nabla u}_i)_j$ denote the pressure, the velocity and the gradient of the velocity at the centroid of the face j and $\varvec{A}_{ij}$ is the outward surface vector of the face j in control volume $\Omega _i$. The values at the face centroids are obtained by linear interpolation of values from cell centers to face centers. Centered and second order accurate schemes are employed to balance numerical stability and solution accuracy for our problem. In particular, the Gauss linear scheme available in OpenFOAM is used for the convective and for the gradient term, while the Gauss linear corrected scheme is exploited for the laplacian term, incorporating a correction term for non-orthogonality in the mesh.

For the time discretization, we adopt an equi-spaced grid with time step $\Delta t = {(T-t_0)}/{N_T}$, where $N_T$ is the number of time intervals in the grid. At each time $t^n=t_0+n\Delta t$ for $n=1,\dots ,N_T$, the time derivative of $\varvec{u}$ is approximated with a backward second order scheme:

$$\begin{aligned} \frac{\partial \varvec{u}}{\partial t} \approx \frac{1}{\Delta t}\left( \frac{3}{2}\varvec{u}^{n+1} -2 \varvec{u}^{n} + \frac{1}{2}\varvec{u}^{n-1}\right) , \end{aligned}$$

(9)

where $\varvec{u}^n$ is the approximation of the velocity at the time step $t^n$.

Let ${\textbf{b}}^{n+1}={(4{\textbf{u}}^n-{\textbf{u}}^{n-1})}/{(2\Delta t)}$. Then, the discretized form of Eq. (1) is:

$$\begin{aligned} \frac{3}{2\Delta t} \varvec{u}_i^{n+1}+\sum _j \phi _{i,j} \varvec{u}_{i,j}^{n+1}-\nu \sum _j(\nabla \varvec{u}_i^{n+1})_j\cdot \varvec{A}_{ij}+\sum _j p_{i,j}^{n+1}\varvec{A}_{ij}=\varvec{b}_i^{n+1}, \quad \phi _{i,j} = \varvec{u}_{i,j}^{n} \cdot \varvec{A}_{ij}, \nonumber \\ \end{aligned}$$

(10)

where $\varvec{u}_i^{n+1}$ and $\varvec{b}_i^{n+1}$ are the average velocity and the source term in $\Omega _i$, and $\varvec{u}_{i,j}^{n+1}$ and $p_{i,j}^{n+1}$ the velocity and pressure at the centroid of face j normalized by the volume of $\Omega _i$. Note that we adopt a first order extrapolation for the convective velocity although we use a backward differentiation formula of order 2 for the time discretization of problem (1)-(2) because this is what OpenFOAM solvers do (Weller et al. 1998).

The discrete version of Eq. (2) is derived from the semi-discretized form of (10), where the pressure term is in continuous form while all other terms are in discrete form. The application of the divergence free constraint and the Gauss-divergence theorem leads to the following equation:

$$\begin{aligned} \sum _k (\nabla p_i^{n+1})_k\cdot \varvec{A}_{ik} = \sum _k \Big ( -\sum _j \phi _{i,j} \varvec{u}_{i,j}^{n+1} + \nu \sum _j(\nabla \varvec{u}_i^{n+1})_j\cdot \varvec{A}_{ij} + \varvec{b}_i^{n+1} \Big )_k \cdot \varvec{A}_{ik}, \nonumber \\ \end{aligned}$$

(11)

where $\phi _{i,j}$ is the convective flux associated to the velocity $\varvec{u}^{n}$ through face j of the control volume $\Omega _i$ defined in (10). More details about the derivation of the equations can be found in Girfoglio et al. (2019), Siena et al. (2023).

The above discrete scheme is implemented in the C++ library . The Pressure Implicit with Splitting of Operators (PISO) algorithm (Issa 1986) is used to decouple the computation of the velocity from the computation of the pressure in coupled problem (10)-(11). For this reason, we need to modify boundary conditions (4) as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} \varvec{u} = \varvec{u}_D & \quad \hbox { in}\ \Gamma _i\times (t_0, T), \\ \varvec{u} = \varvec{0} & \quad \hbox { in}\ \Gamma _w\times (t_0, T), \\ (\nabla \varvec{u})\cdot \varvec{n} = \varvec{0} & \quad \hbox { in}\ \Gamma _o\times (t_0, T), \\ (\nabla p) \cdot \varvec{n}= \varvec{0} & \quad \hbox { in}\ (\Gamma _i\cup \Gamma _w)\times (t_0, T), \\ p = 0 & \quad \hbox { in}\ \Gamma _o\times (t_0, T). \end{array}\right. } \end{aligned}$$

(12)

The method described in this subsection assumes that the mesh for space discretization is fine enough to capture all the flow structures, i.e., eddies and vortices. When such a refined mesh cannot be afforded with the given computing facilities and one is forced to use a coarser mesh, then an alternative approach, like Large Eddy Simulation (LES), should be used. For the specific benchmark problem considered in Sect. 4, one easy to implement LES method is presented in Bertagna et al. (2016), Girfoglio et al. (2019).

3 Reduced order model

The ROM adopted for this work is the POD-Galerkin approach (Quarteroni et al. 2011, 2015; Hesthaven et al. 2016; Rowley et al. 2004; Iollo et al. 2000) with stabilization techniques to recover the effect of the neglected modes (Cazemier 1997; Bergmann et al. 2009; San and Iliescu 2015). Here, we briefly recall the main ideas of this methodology.

We assume that the solution $(\varvec{u}, p)$ to problem (1)-(2) can be approximated as a linear combination of basis functions $\left( \varvec{\phi }_i(\varvec{x}), \psi _i(\varvec{x})\right)$ that depend on space only and coefficients $\left( a_i(t, \varvec{\mu }), b_i(t, \varvec{\mu }) \right)$ that depend on time and $\varvec{\mu }$ only:

$$\begin{aligned} \varvec{u} (\varvec{x}, t, \varvec{\mu }) \approx \varvec{u}_{\text {r}}(\varvec{x}, t, \varvec{\mu }) = \sum _{i=1}^{N_{\varvec{u}}} a_i(t, \varvec{\mu }) \varvec{\phi }_i(\varvec{x}), \end{aligned}$$

(13)

$$\begin{aligned} p (\varvec{x}, t, \varvec{\mu }) \approx p_{\text {r}}(\varvec{x}, t, \varvec{\mu }) = \sum _{i=1}^{N_p} b_i(t, \varvec{\mu }) \psi _i(\varvec{x}). \end{aligned}$$

(14)

This assumption allows us to decouple the computations into an expensive phase (offline) to be performed only once and a cheap phase (online) to be performed for every new time and parameter of interest. We recall that $\varvec{\mu }$ is a vector containing all the physical parameters the problem depends upon (see Sect. 2.1).

The specific operations performed during each phase are:

• Offline: given a set of time instances and parameters, the full order solutions (also called snapshots) are computed and collected into a matrix. The lifting function method is adopted to treat non-homogeneous boundary conditions (see Sect. 3.3). The POD algorithm is used to extract the reduced basis space (see Sect. 3.1). The Pressure Poisson Equation (PPE) approach is adopted to satisfy the inf-sup condition.

• Online: given a set of time and parameter values, the corresponding modal coefficients are computed by solving the reduced system. For stabily, an artificial viscosity is added in that system with different approaches (see Sect. 3.2). The approximated solution (13)-(14) is a linear combination between these coefficients and the POD reduced basis functions computed offline.

3.1 Proper orthogonal decomposition

To generate a set of orthonormal basis functions in a least-square setting, we use the POD algorithm (Eckart and Young 1936; Hawkins and Ben-Israel 1973; Bang-Jensen et al. 2004). Below, we summarize how the algorithm works.

Given a discrete set of time instances $\{ t^1, \dots , t^{{N_T}} \} \subset (t_0, T)$ and parameters $\{ \varvec{\mu }^1, \dots , \varvec{\mu }^{{N_{\mu }}} \}$, the corresponding FOM solutions $\Phi =\{p,\varvec{u}\}$ of problem (1)-(2) are stored as the columns of a snapshot matrix $S_{\Phi }$:

$$\begin{aligned} {\mathcal {S}}_{\Phi } = \begin{Bmatrix} \Phi (t^1, \varvec{\mu }^1),&\cdots&, \Phi (t^{N_T}, \varvec{\mu }^1), \Phi (t^1, \varvec{\mu }^2),&\cdots&, \Phi (t^{N_T}, \varvec{\mu }^2),&\cdots&, \Phi (t^{N_T}, \varvec{\mu }^{N_{\mu }}) \end{Bmatrix}. \end{aligned}$$

We are describing the algorithm with generic variable $\Phi$ to avoid repetition since POD is performed separately for the velocity and pressure fields.

Let $N_{\Phi } \ll \text {min}(N_{h}, {M})$ be the size of the reduced basis, with ${M = N_T\cdot N_{\mu }}$. We will explain in Eq. (17) how $N_{\Phi }$ is chosen. The reduced basis is the solution ${\mathcal {V}} = \{\varvec{\phi }_1,\dots ,\varvec{\phi }_{N_{\varvec{\Phi }}}\}$ of the optimization problem (Kunisch and Volkwein 2002; Quarteroni et al. 2015):

$$\begin{aligned} \min _{{\mathcal {V}}} \Vert {\mathcal {S}}_{\Phi }- {\mathcal {V}}{\mathcal {V}}^T{\mathcal {S}}_{\Phi } \Vert \quad s.t. \quad {\mathcal {V}}^T{\mathcal {V}}={\mathcal {I}}. \end{aligned}$$

(15)

Problem (15) seeks to minimize the distance between the snapshots and their projection onto the reduced (POD) space. Note that it can be equivalently formulated as an eigenvalue problem (Kunisch and Volkwein 2002), which is easier to handle:

$$\begin{aligned} C_{\Phi } \varvec{c}^{\Phi }_s = \lambda _s^{\Phi } \varvec{c}_s^{\Phi }, \quad s = 1, \dots , {M}. \end{aligned}$$

(16)

In (16), $C_{\Phi }=\frac{1}{N_{T}}S_{\Phi }^TS_{\Phi }\in {\mathbb {R}}^{{M}\times {M}}$ is the correlation matrix related to the snapshots, and $\varvec{c}^{\Phi }_s$ and $\lambda _s^{\Phi }$ are the corresponding eigenvectors and eigenvalues. Then, the POD basis function are computed as follows:

$$\begin{aligned} \varvec{\phi }_i=\frac{1}{\sqrt{\lambda ^{\Phi }_i}} S_{\Phi }\varvec{c}^{\Phi }_i, \quad i = 1,\dots , N_{\Phi }. \end{aligned}$$

The dimension $N_{\Phi }$ is chosen by comparing the cumulative energy of the eigenvalues with a threshold $\delta$ selected by the user:

$$\begin{aligned} \frac{\sum _{i=1}^{N_{\Phi }}\lambda _i^{\Phi }}{\sum _{i=1}^{{M}}\lambda _i^{\Phi }} \ge \delta . \end{aligned}$$

(17)

The error made by approximating the snapshots matrix $S_{\Phi }$ with the POD basis is given by the sum of the neglected singular values (Quarteroni et al. 2015):

$$\begin{aligned} \Vert {\mathcal {S}}_{\Phi }- {\mathcal {V}}{\mathcal {V}}^T{\mathcal {S}}_{\Phi } \Vert ^2 = \sum _{i=N_{\Phi }+1}^{\text {min}(N_{h},{M})} \lambda _i^{\Phi }. \end{aligned}$$

(18)

Therefore, by varying $N_{\Phi }$, the ROM solution can approximate the FOM solution with arbitrary accuracy. However, we note that there is no guarantee that the $N_{\Phi }$ needed to achieve the desired accuracy is small when compared to $N_{h}$ or M.

3.2 Galerkin projection

The $L^2$ orthogonal projection of the momentum equation (1) onto $\text {span}\{\varvec{\phi }_1,\dots ,\varvec{\phi }_{N_{\varvec{u}}}\}$ leads to Akhtar et al. (2009), Bergmann et al. (2009), Lorenzi et al. (2016):

$$\begin{aligned} \big (\varvec{\phi }_i, \partial _t \varvec{u} + \nabla \cdot (\varvec{u} \otimes \varvec{u}) - \nabla \cdot (\nu \nabla \varvec{u})+\nabla p \big )_{L^2(\Omega )} = 0, \quad \text {for }i = 1, \dots , N_{\varvec{u}}. \end{aligned}$$

(19)

If we plug (13) and (14) into (19) and take into account the orthonormality of the reduced basis, we get the following ordinary differential equation in matrix form:

$$\begin{aligned} \dot{\varvec{a}} = \nu \varvec{B} \varvec{a} - \varvec{a}^T \varvec{C} \varvec{a} - \varvec{K} \varvec{b}, \end{aligned}$$

(20)

where $\varvec{a} = \{ a_i(t)\}_{i=1}^{N_u}$ and $\varvec{b} = \{ b_i(t)\}_{i=1}^{N_p}$ are the ROM coefficients for velocity and pressure and:

$$\begin{aligned} B_{ij}&= \big ( \varvec{\phi }_i, \Delta \varvec{\phi }_j \big )_{L^2(\Omega )}, \quad \text {for } i, j = 1, \dots N_{\varvec{u}}, \end{aligned}$$

(21)

$$\begin{aligned} C_{ijk}&= \big ( \varvec{\phi }_i, \nabla \cdot (\varvec{\phi }_j \otimes \varvec{\phi }_k) \big )_{L^2(\Omega )}, \quad \text {for } i, j, k = 1, \dots N_{\varvec{u}}\end{aligned}$$

(22)

$$\begin{aligned} K_{ij}&= \big ( \varvec{\phi }_i, \nabla \psi _j \big )_{L^2(\Omega )}, \quad \text {for } i= 1, \dots N_{\varvec{u}} \text { and } j =1, \dots , N_{p}. \end{aligned}$$

(23)

To handle the divergence free constraint (2), we adopt the PPE approach (Hesthaven and Ubbiali 2018; Guermond and Quartapelle 1998; Stabile and Rozza 2018). Therefore, we replace Eq. (2) with

$$\begin{aligned}&\Delta p = - \nabla \cdot ( \nabla \cdot (\varvec{u} \otimes \varvec{u})) \quad \hbox { in}\ \Omega , \end{aligned}$$

(24)

$$\begin{aligned}&\frac{\partial p}{\partial \varvec{n}} = -\nu \varvec{n} \cdot (\nabla \times \nabla \times \varvec{u}) \quad \hbox { on}\ \partial \Omega . \end{aligned}$$

(25)

Formulation (24)-(25) is obtained by applying the divergence operator to the momentum equation and by exploiting the incompressibility constraint, under suitable hypotheses of regularity. See, e.g., Stabile et al. (2017) for more details.

The projection of Eq. (24) onto $\text {span}\{\psi _1, \dots , \psi _{N_p}\}$ leads to:

$$\begin{aligned} \big (\nabla \psi _i,\nabla p \big )_{L^2(\Omega )} - \big (\psi _i,\nabla p \varvec{n} \big )_{L^2(\partial \Omega )}= \big ( \psi _i, \nabla \cdot (\nabla \cdot (\varvec{u} \otimes \varvec{u})) \big )_{L^2(\Omega )}, \quad \text {for }i = 1, \dots , N_{p}, \nonumber \\ \end{aligned}$$

(26)

where the left-hand side is obtained integrating by parts and accounting for boundary condition (25). If we plug (13) and (14) into (26), we obtain the matrix form:

$$\begin{aligned} \varvec{D} \varvec{b} - \varvec{N} \varvec{b} = \varvec{a}^T \varvec{G} \varvec{a}, \end{aligned}$$

(27)

where

$$\begin{aligned} D_{ij}&=\big ( \nabla \psi _i, \nabla \psi _j \big )_{L^2(\Omega )}, \quad \text {for } i, j = 1, \dots N_p,\end{aligned}$$

(28)

$$\begin{aligned} N_{ij}&=\big ( \psi _i, \nabla \psi _j \varvec{n} \big )_{L^2(\partial \Omega )}, \quad \text {for } i, j = 1, \dots N_p,\end{aligned}$$

(29)

$$\begin{aligned} G_{ijk}&=\big ( \psi _i, \nabla \cdot (\nabla \cdot (\varvec{\phi }_j\otimes \varvec{\phi }_k)) \big )_{L^2(\Omega )}, \quad \text {for } i = 1, \dots N_p \text { and } j, k = 1, \dots N_{\varvec{u}}. \end{aligned}$$

(30)

We note that all the matrices in (20) and (27) are computed once and for all during the offline stage. Computational challenges associated with nonlinear terms ($\varvec{C}$ and $\varvec{G}$) can be addressed by limiting the number of POD basis functions. In fact, notice that the dimension of tensor $\varvec{C}$ grows as $N_{\varvec{u}}^3$, while dimension of tensor $\varvec{G}$ grows as $N_{\varvec{u}}^2 N_p$. Equations (20) and (27) are solved with Powell’s dog leg method (Powell 1970), an iterative optimization algorithm for the solution of non-linear problems.

The ROM framework introduced so far is able to capture the dynamics of the system if the first $N_{\varvec{u}}$ and $N_p$ modes, respectively for the velocity and the pressure, is large enough to capture the FOM solutions with accuracy. For a convection dominated flow, $N_{\varvec{u}}$ is typically very large, which leads to high computational cost. However, if one keeps $N_{\varvec{u}}$ e $N_p$ low to reduce the computational cost, the accuracy in the reconstruction of the flow at the reduced order level is severely compromised. One possible way to recover information lost to mode truncation is to introduce a stabilization method.

The easiest way to achieve stabilization is by adding a global constant viscosity (Bergmann et al. 2009; Kalb and Deane 2007; Wang et al. 2012), i.e., the viscosity in the projected momentum equation (20) is modified as follows:

$$\begin{aligned} \nu \quad \Rightarrow \quad \nu (1 + \nu _a), \end{aligned}$$

(31)

where $\nu _a$ is a scaling factor. Although trial and error is not the best way to proceed to set $\nu _a$, we would like to remark that (31) is performed only at the reduced order level, i.e., on a system whose size is contained. Hence, though cumbersome, the procedure is not computationally expensive: it takes only a total of 3 s.

A possible easy improvement over (31) is to modify the amount of added viscosity for each mode (San and Iliescu 2015). The idea consists in replacing the viscosity in the projected momentum equation (20) as follows:

$$\begin{aligned} \nu \quad \Rightarrow \quad \nu \left( 1 + \nu _a\frac{k}{N_{\varvec{u}}}\right) , \end{aligned}$$

(32)

where k is the modal index. This means that more diffusivity is added for the modes with less energy content in the system.

3.3 Lifting function method

Non-homogeneous boundary conditions are introduced at reduced level with the lifting function (also called control function) method (Stabile et al. 2017; Fick et al. 2018; Graham et al. 1999; Gunzburger et al. 2007). This entails building a homogeneous reduced basis space by removing non-homogeneous boundary value from the collected snapshots. Then, the non-homogeneous boundary value is added again during the online phase, when the ROM solution is computed.

Let $\varvec{\chi }(\varvec{x})$ be a lifting function, i.e., a divergence free field (like the modes) that takes the value $\varvec{u}_D (\varvec{x})$ where the non-homogeneous boundary condition is enforced. The snapshots are scaled as follows:

$$\begin{aligned} \varvec{u}' (t^n, {\varvec{\mu }^m}) = \varvec{u} (t^n, {\varvec{\mu }^m}) - \varvec{\chi }(\varvec{x}), \end{aligned}$$

(33)

where $\varvec{u} (t^n, {\varvec{\mu }^m})$ is the collected snapshot and $\varvec{u}' (t^n, {\varvec{\mu }^m})$ is the “homogenized” snapshot. A common approach to compute the lifting function $\varvec{\chi }$ is to solve a potential flow problem. This means that if $\Gamma _i$ is the boundary where the non-homogeneous Dirichlet condition is imposed, then $\varvec{\chi }$ is found by solving the following system:

$$\begin{aligned} {\left\{ \begin{array}{ll} \nabla \cdot \varvec{u} = 0 & \quad \hbox { in}\ \Omega , \\ \varvec{u} = \nabla \varvec{\chi } & \quad \hbox { in}\ \Omega , \\ \varvec{\chi } = \varvec{u}_D (\varvec{x}) & \quad \hbox { in}\ \Gamma _i,\\ \varvec{\chi } = \varvec{0} & \quad \hbox { in}\ \Gamma _w,\\ \nabla \varvec{\chi } \cdot \varvec{n} = 0 & \quad \hbox { in}\ \Gamma _o. \end{array}\right. } \end{aligned}$$

Snapshot $\varvec{u}' (t^n, {\varvec{\mu }^m})$, which satisfies homogeneous boundary conditions, is then used to generate the reduced basis as explain in Sect. 3.1. Then, (13) is updated with:

$$\begin{aligned} \varvec{u} (\varvec{x}, t^n, {\varvec{\mu }^m})\approx \varvec{u}_{\text {r}}(\varvec{x}, t^n, {\varvec{\mu }^m}) =\varvec{\chi }(\varvec{x})+ \sum _{i=1}^{N_{\varvec{u}}} a_i(t^n, {\varvec{\mu }^m}) \varvec{\phi }_i(\varvec{x}). \end{aligned}$$

(34)

4 Numerical results

To test the ROM framework described in Sect. 3, we choose the FDA Benchmark 1, i.e., 3D flow in a convergent-divergent nozzle for low-to-moderate Reynolds numbers. The interest in the nozzle geometry, a cross-section of which is shown in Fig. 1, lies in the fact that its features are commonly encountered in medical devices. The advantage of this benchmark is the availability of experimental measurements² against which one can compare the numerical results. The FDA collected the experimental data from three independent laboratories, some of which performed more than one trial, for a total of five data sets.

Fig. 1

Top: Schematic representation of the nozzle geometry used for the FDA benchmark, where $D_i=0.012$ m and $D_t=0.004$ m. Bottom: Identification of the different parts of the geometry

The density and viscosity of the fluid are set to $\rho = 1056$ kg/$\hbox {m}^3$ and $\mu = 0.0035$ Pa/s, with $\nu =\mu /\rho =3.31\cdot 10^{-6}$ (Pa $\cdot$ m) / (s$\cdot$Kg), to match the fluid used in the experiments. A Poiseuille velocity profile is imposed at the inflow boundary (i.e., at $z=0$ m):

$$\begin{aligned} \varvec{u}_D (r, \theta , z) = \left( 0, 0, 2 V_{\text {mean}}\left( 1-\frac{r^2}{R_i^2}\right) \right) , \end{aligned}$$

(35)

where r, $\theta$, and z are the radial, the polar, and axial coordinates, $R_i = 0.012$ m is the inlet radius and $V_{\text {mean}}$ is the mean inlet velocity magnitude to obtained the desired Reynolds number in the throat. We will consider Reynolds numbers in the throat of the nozzle equal to 2000, 3500, 5000 and 6500, i.e., transitional to turbulent regimes. We compute the Reynolds number as follows:

$$\begin{aligned} Re = \frac{V_{\text {mean}}D_t}{\nu }, \end{aligned}$$

(36)

where $D_t=0.004$ m is the diameter of the throat and $\nu$ is the kinematic viscosity of the fluid. We will show that the FOM solutions compare well with the experimental data and then compare the ROM solutions with the FOM solutions.

Table 1 Reynolds number Re in the throat, Kolmogorov length scale $\eta$ and number of cells for a DNS $N_D$ for the flow regimes under consideration

Table 1 reports the Kolmogorov scale $\eta$:

$$\begin{aligned} \eta = Re^{-3/4}D_i, \end{aligned}$$

(37)

and the number of cells $N_D$ required for a DNS for the three flow regimes under consideration. The value of $N_D$ is obtained by assuming a uniform mean grid width equal to $\eta$. As Table 1 shows, a DNS for the Reynolds numbers under consideration requires meshes with millions of cells. See, e.g., (Bertagna et al. 2016; Zmijanovic et al. 2017). In this paper, we consider a much coarser mesh whose features are reported in Table 2. In general, a simulation with this mesh would lead to non-physical oscillations because the mesh size is too large to capture the small scale eddies. However, our FV method gives numerical results that compare well with the experimental measurements even with a mesh much coarser than a DNS requires. A possible reason for the better performance of our FV method with respect to, e.g., the FE method in Bertagna et al. (2016), could be the following: it yields exact conservation (see discussion in Sect. 4.1) and thus provides acceptable results despite the use of a coarse mesh. Furthermore, it might be that the diffusive effect of FV contributes to eliminate numerical oscillations in the FOM solution. However, by employing second-order spatial discretization schemes and second-order backward time discretization with a small time step ($\simeq 10^{-4}$), the diffusion introduced by the FV method is minimized.

Table 2 Features of the mesh used for all the simulations: number of cells, minimum cell size $h_{\text {min}}$, maximum cell size $h_{\text {max}}$, maximum degree of non-orthogonality, and maximum skewness

Figure 2 compares the results obtained with our coarse mesh and the experimental data for all the Reynolds number under consideration. Let us comment first on the comparison for the pressure on the left column in Fig. 2. We observe that for $Re=2000, 3500, 5000$ the full order pressure is in excellent agreement with most experimental data in the entrance channel, conical convergent, and throat. Note that no pressure data in the expansion channel are provided by the FDA for any Reynolds number. For $Re = 6500$, the computed pressure overestimates the measured pressure in the conical convergent and throat (bottom left panel in Fig. 2). Now, let us look at the comparison for the axial velocity on the right column in Fig. 2. For all three Reynolds numbers, we see a good agreement with the experimental data in the entrance channel and conical convergent. However, the computed axial velocity tends to be lower than most data in the throat for all Re. This becomes more evident as the Reynolds number increases. Instead, the length of the jet in the expansion channel is better predicted as the Reynolds number increases. One could improve these results with LES (Girfoglio et al. 2019; Zmijanovic et al. 2017; Bertagna et al. 2016; Janiga 2014; Delorme et al. 2013; Fehn et al. 2019) or RANS (Stewart et al. 2012) models. Since the purpose of this paper is to show how to reconstruct a convection-dominated flow at the reduced level, the results given by our FV approach on a coarse mesh with no LES or RANS model are an equally good starting point. As a further proof of the accuracy of the data used to train the ROM, Fig. 3 compares the measure and computed time-averaged radial velocity for $Re = 3500$ at different radial sections. We observe a great match for the radial velocity at the sections in the entrance channel, narrow throat, and past the jet breakdown point. Note that, despite the mesh coarseness, we can capture well the plug-like velocity profile in the throat. We observe some differences near the sudden expansion ($z=0.126685$ m) in the form of larger recirculations near the walls and a more smoothed-out profile for the computed time-averaged radial velocity. Overall, the results in Fig. 3 confirm that our FOM solutions are good enough to train our ROMs. A final note on the data in Fig. 3, top right panel: the throat is narrower than the entrance channel, hence any value outside the actual geometry should be zero. That is not the case for one of the datasets, which must be an error in the dataset itself.

Fig. 2

Comparison between experimental data (dashed lines) and the FOM solution (solid line) for the time-averaged axial pressure (left) and velocity (right) for different Reynolds numbers

Fig. 3

Comparison between experimental data (dashed lines) and FOM solution (solid line) for the time-averaged radial velocity for $Re = 3500$ at different radial sections

In Sect. 4.1, we present a thorough analysis of the results obtained with the methods presented in Sect. 3.2 for $Re = 3500$, with time as the only parameter. Then, in Sect. 4.2 we show that the conclusions drawn for the $Re = 3500$ case apply also for the transitional case ($Re = 2000$) and the highest Reynolds number case ($Re = 6500$). In Sect. 4.3, the Reynolds number is treated as a parameter, in addition to time.

The time interval of interest for all the simulations is $(t_0, T) = (0, 1.2)$ s, with the fluid at rest at $t_0$. For the FOM, we set $\Delta t = 0.0001$ s and the maximum CFL computed during the simulation is 0.6, as in Girfoglio et al. (2019), Zmijanovic et al. (2017). Starting from the given value of $\Delta t$, OpenFOAM automatically adjusts its during the simulation to ensure that the maximum Courant number does not exceed 0.6.

4.1 Case $Re=3500$

In order to obtain the desired Reynolds number, we set $V_{\text {mean}} = 0.32225$ m/s. For the POD algorithm, we use with 120 equispaced snapshots in time interval (0, 1.2) s.

The cumulative energy of the eigenvalues (17) is shown in Fig. 4a for velocity and pressure. While the pressure reaches $99.92\%$ of the energy with a single mode, the velocity needs 50 modes to achieve only $83.73\%$. The slow growth of velocity eigenvalues is reflected also in the mean relative error (18) for velocity and pressure shown in Fig. 4b. We see that for both variables the use of more than 9 modes does not significantly improve the reconstruction.

Fig. 4

(a) Cumulative energy of the eigenvalues (17) and (b) mean relative error (18) for pressure and velocity

Figure 5 shows the time-averaged axial velocity and pressure along the z-axis by using 9 modes for both variables and no stabilization. The large distance between the FOM and ROM velocity curves reflects the large error in Fig. 4b, while the FOM and ROM pressure curves are practically superimposed. The qualitative comparison of the solutions is reported in Figs. 6 and 7, which confirm that the ROM reconstruction of the pressure is accurate, whereas the ROM velocity shows an evident mismatch with the FOM velocity.

Fig. 5

Comparison of FOM and ROM time-averaged axial velocity (m/s) and pressure ($\hbox {m}^2$/$\hbox {s}^2$) using 9 modes for both variables and no stabilization technique

Fig. 6

Qualitative comparison of the ROM time-averaged pressure without stabilization with the corresponding FOM solution

Fig. 7

Qualitative comparison of the ROM time-averaged axial velocity without stabilization with the corresponding FOM solution

Fig. 8

Comparison of FOM and ROM time-averaged axial velocity (m/s) and pressure ($\hbox {m}^2$/$\hbox {s}^2$) with constant artificial viscosity for (a) $\nu _a=0.8$, (b) $\nu _a=1.5$, (c) $\nu _a=20$, and (d) $\nu _a=500$

Next, we will show how the ROM solution improves with the introduction of artificial viscosity only at the reduced order level. First, we will present the results for strategy (31) and then proceed with strategy (32).

Constant artificial viscosity. Figure 8 shows the time-averaged axial velocity and pressure along the z-axis as the parameter $\nu _a$ increases. We see that the ROM reconstruction of the velocity exhibits oscillations for small values of $\nu _a$ (see Fig. 8a, b). In addition, for $\nu _a=0.8$ the velocity is quite underestimated in the expansion channel. This is overcome with $\nu _a=1.5$: the ROM velocity curve gets closer to the FOM curve, however oscillations are still noticeable. Increasing $\nu _a$ to 20 removes all oscillations. However, for this value of $\nu _a$ we observe an unexpected increase in the time-averaged axial velocity in the expansion channel and the time-averaged ROM pressure moves away from the time-averaged FOM pressure in the entrance channel. See Fig. 8c. The best ROM solutions is offered by $\nu _a=500$, with the exception of the ROM reconstruction of the time-averaged pressure in the entrance chamber. See Fig. 8d. We remark that this trial and error process to select the value of $\nu _a$ is performed during the online phase: the ROM simulation for a given value of $\nu _a$ takes only about 3 s. For all the test cases analyzed in the paper, we began by manually testing a broad range of values $\nu _a$. After identifying a small value of $\nu _a$ that produces oscillations in the ROM predictions but maintains reasonable agreement with the FOM solution, we refined our exploration within a range centered around this value. Since typically the oscillations could not be eliminated for small $\nu _a$, we increased $\nu _a$ until achieving a sufficiently smooth velocity profile while keeping small the velocity and pressure errors between the FOM and the ROM solutions. The “best” $\nu _a$ was therefore selected by minimizing the mean velocity error, while maintaining an accurate approximation of the pressure.

An important observation is in order: even though $\nu _a=500$ gives a good approximation of the time-averaged FOM solution, the ROM solution becomes steady. Table 3 shows the percentage of time dependency (number of time dependent snapshots over the total number of snapshots) as $\nu _a$ increases. We see that for $\nu _a > 20$, the ROM solution can be considered steady. This conclusion is based on a visual evaluation of the reconstructed solution, which shows no apparent temporal variations. This observation does not influence the validity of our results because we are interested in the reconstruction of a time-averaged value for pressure and velocity. However, this approach cannot be used to recover the time-dependent solution.

Table 3 Percentage of time dependency (number of time dependent snapshots over the total number of snapshots) of the ROM solutions, when constant artificial viscosity strategy is adopted

In order to investigate quantitatively the effect of $\nu _a$, we present in Fig. 9 the absolute and relative errors of the time-averaged axial velocity and pressure along the z-axis. The errors of the velocity in Fig. 9a, c confirm that $\nu _a=500$ provides a very accurate ROM solution. The relative errors of the pressure in Fig. 9b show a significant peak around $z=0.17$, also for $\nu _a=500$. Nevertheless, the absolute errors in Fig. 9d show that the difference between FOM and ROM pressure is only about $10^{-2}$, therefore the peak of the relative error is purely due to the low pressure value near $z=0.17$. The curves in Fig. 9b, d show a degeneration in the ROM pressure in the entrance channel as $\nu _a$ is increased. However, we note that the relative error for $\nu _a=20, 500$ is only about $10^{-1}$, which is a rather good compromise to obtain satisfactory approximations for both variables.

Fig. 9

Relative (top row) and absolute (bottom row) errors of the ROM axial velocity (left) and pressure (right) with respect to the FOM solution for constant artificial viscosity $\nu _a = 0, 20, 500$

For a visual comparison, the time-averaged pressure and velocity fields for $\nu _a=500$ are shown in Figs. 10 and 11, respectively.

Fig. 10

Qualitative comparison of the ROM time-averaged pressure obtained with constant artificial viscosity $\nu _{a}=500$ with the corresponding FOM solution

Fig. 11

Qualitative comparison of the ROM time-averaged axial velocity obtained with constant artificial viscosity $\nu _{a}=500$ with the corresponding FOM solution

Since the pressure can be seen as a Langrange multiplier for the enforcement of the incompressibility constraint for the velocity and given the discrepancies between FOM and ROM pressure in the entrance channel, we checked if the ROM solution conserves mass using the following error (Passerini et al. 2013):

$$\begin{aligned} E^\text {mass} = \frac{Q_{\text {FOM}} (z)-Q_{\text {ROM}} (z)}{Q_{\text {FOM}}(z)}, \end{aligned}$$

(38)

where Q is the flow rate. Table 4 reports error (38) for $\nu _a=500$ and in the case of no stabilization ($\nu _a=0$) for given values of the axial coordinate z. We see that the flow rate computed by the ROM with $\nu _a=500$ is within about 1% of the flow rate computed by the FOM, indicating that mass is conserved reasonably well. Instead, when $\nu _a=0$ (no stabilization) error (38) goes up to about 10% in the expansion channel.

Finally, to understand if there is an “optimal” value of $\nu _a$, we plot in Fig. 12 the time- and space-averages relative errors for pressure and axial velocity as $\nu _a$ increases. We observe that both errors peak for $\nu _a=4$ and then reach a plateau for $\nu _a>100$. For $\nu _a>4$ the error of the axial velocity decreases monotonically, whereas the error of the pressure rises slightly after attaining a minimum for $\nu _a=20$, This is consistent with the previous observation of increasing distance between the FOM and ROM pressure curves in the entrance channel in Fig. 8c, d.

Fig. 12

Time- and space-averaged relative errors of the ROM axial velocity and pressure with respect to the FOM solution as constant artificial viscosity $\nu _a$ varies

Modal coefficient-dependent artificial viscosity. Let us present the results obtained with the varying artificial viscosity (32). Our aim is to determine whether it represents an improvement over the constant artificial viscosity approach.

Figure 13 shows the time-averaged axial velocity and pressure along the z-axis for different values of $\nu _a$. Based on the results obtained with the constant artificial viscosity, we have analyzed the ROM solutions for $\nu _a=1.5, 4, 20, 100, 500, 1000$. While $\nu _a=1.5$ provides an accurate ROM solution in the case of constant artificial viscosity (see Fig. 8b), in the case of coefficient-dependent artificial viscosity it leads to an underestimation of the axial velocity in the expansion channel (see Fig. 13a). By increasing the artificial viscosity to $\nu _a=4$, the ROM axial velocity profile gets closer to the FOM counterpart (see Fig. 13b), however it displays rather large oscillations. In particular, we notice larger oscillations with the coefficient-dependent artificial viscosity for $\nu _a=4$ (see Fig. 13b) than with a constant artificial viscosity for $\nu _a=1.5$ (see Fig. 8b). Furthermore, when $\nu _a$ is increased to 20, the coefficient-dependent viscosity approach gives a worse prediction than the constant viscosity approach: more oscillations are present in Fig. 13c than Fig. 8c and the distance between FOM and ROM solutions is significantly higher. In order to obtain similar results to the constant viscosity case for $\nu _a=20$, one needs to increase $\nu _a$ to 100 with the coefficient-dependent artificial viscosity approach. We see that in this case too the ROM axial velocity is too high in the expansion channel and the ROM pressure underestimates the FOM pressure in the entrance channel (see Fig. 13d). In order to improve the accuracy of the ROM solution, we need to further increase the artificial viscosity. A good reconstruction of the velocity along the whole domain can be obtained with $\nu _a=500, 1000$ (see Fig. 13e, f). We observe no apparent distinctions in the ROM solution for these two values.

Fig. 13

Comparison of FOM and ROM time-averaged axial velocity (m/s) and pressure ($\hbox {m}^2$/$\hbox {s}^2$) with coefficient-dependent viscosity for (a) $\nu _a = 1.5$, (b) $\nu _a = 4$, (c) $\nu _a = 20$, (d) $\nu _a = 100$, (e) $\nu _a = 500$, and (f) $\nu _a = 1000$

We deduce that once a certain threshold is exceeded, $\nu _a$ ceases to affect the ROM solution. This result is confirmed in Fig. 14, which shows the time- and space-averaged relative errors for axial velocity and pressure as $\nu _a$ is increased. Similarly to the constant artificial viscosity approach, we see that both time- and space average- relative errors reach a plateau for $\nu _a>500$. Another similarity is the increase in the error of the pressure after reaching a minimum for $\nu _a>50$. A difference with respect to Fig. 12 is that in Fig. 14 the peaks in the errors for axial velocity and pressure do not occur for the same value of $\nu _a$.

Fig. 14

Time- and space-averaged relative errors of the ROM axial velocity and pressure with respect to the FOM solution as $\nu _a$ varies in the coefficient-dependent artificial viscosity approach

Figure 15 reports the relative and absolute errors for axial velocity and pressure as $\nu _a$ increases. In line with the previous considerations, the errors are similar to the constant artificial viscosity case when $\nu _a$ is larger. Compare the curves in Fig. 9 with the curves in Fig. 15. The curves for the pressure in Fig. 15b, d show similar features to Fig. 9b , d: the error rises in the entrance chamber as $\nu _a$ increases, reaching an error of $10^{-1}$ for $\nu _a = 1000$, which is a good compromise to have a satisfactory ROM prediction for both variables. The peak in the relative error for the pressure at $z=0.17$ is still present, but the corresponding absolute error is small ($\approx 10^{-4}$), so the large value of the relative error is due to the pressure being close to zero.

Fig. 15

Relative (top row) and absolute (bottom row) errors of the ROM velocity (left) and pressure (right) with respect to the FOM solution for $\nu _a = 0, 100, 1000$ in the coefficient-dependent viscosity approach

Given that we get a mismatch between the ROM pressure and the FOM pressure in the entrance channel in the case of coefficient-dependent viscosity too, we checked mass conservation. In Table 5, we report error $E^\text {mass}$ (38) at given values of the axial coordinate z for the coefficient-dependent approach at $\nu _a = 500$. Notice that the values in Table 5 are very similar to the values on the second row of Table 4.

The coefficient-dependent viscosity approach does not fix the problem of losing time dependency in the ROM solution as the value of $\nu _a$ increases. Table 6 reports the percentage of time dependency (number of time dependent snapshots over the total number of snapshots) as $\nu _a$ varies. We see that time evolution persists for higher value of $\nu _a$ when compared with the constant artificial viscosity case. However, for $\nu _a > 50$ the ROM reconstruction is close to being steady. Therefore, the coefficient-dependent viscosity approach cannot be used to predict the ROM solution in time as well.

Table 4 Flow rate computed from the FOM solution and error (38) for constant artificial viscosity $\nu _a = 500, 0$ and given values of the axial coordinate z

Table 5 Error (38) at given values of the axial coordinate z for $\nu _a = 500$ with coefficient-dependent artificial viscosity

Table 6 Percentage of time dependency (number of time dependent snapshots / total number of snapshots) of the ROM solutions when then coefficient-dependent viscosity strategy is adopted

The analysis of the results in this subsection clarifies that it is possible to use both a constant artificial viscosity (31) and a coefficient-dependent artificial viscosity (32) to obtain the time average velocity and pressure in a turbulent flow. These two approaches are equivalent, in the sense that neither features a particular advantage over the other. Both procedure, in combination with POD-Galerkin and PPE, are not able to recover a time dependent solution for the turbulent flow and both provide an accurate reconstruction of the time-averaged solution with a sufficient large value of $\nu _a$. For completeness, we tried to use both methods also in combination with the supremizer enrichment (Ballarin et al. 2015; Gerner and Veroy 2012; Rozza and Veroy 2007). However, the distance between FOM and ROM mean pressure increases significantly in the entrance channel when using supremizer enrichment, with a consequent higher error (38). For this reason, the results obtained with the supremizer enrichment are not shown.

Finally, we comment on the computational cost using an Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz 2:30 16GB RAM. Each FOM simulation takes 5 hours and 17 seconds are required for the POD algorithm to complete, while only 1 s is needed for the evaluation phase. Therefore, the speed-up is significant. The reader interested in learning more about the computational savings allows by ROMs in the context of different FOMs is referred to, e.g., (Peter et al. 2021; Benner et al. 2020a, b; Hesthaven et al. 2016; Rozza et al. 2008).

4.2 Cases $Re=2000, 6500$

The purpose of this section is to show that the conclusions on the artificial viscosity and coefficient-dependent artificial viscosity methods drawn for the flow at $Re = 3500$ can be extended to the transitional case ($Re=2000$) and the highest Reynolds number case considered by the FDA ($Re=6500$). Following the previous section, we will present both qualitative and quantitative results.

In the simulations, we set $V_{\text {mean}} = 0.3679616$ m/s to obtain $Re=2000$ and $V_{\text {mean}} = 1.1969$ m/s to reach $Re=6500$. Figure 16 shows the cumulative energy of the eigenvalues (17) and the average relative error (18) for velocity and pressure for $Re=2000, 6500$. We truncate the number of mode to 9, which is clearly not sufficient to provide an accurate reconstruction of the flow velocity and pressure in either case if one uses a standard POD-Galerkin ROM.

Fig. 16

Cumulative energy of the eigenvalues (17) at (a) $Re=2000$, (c) $Re=6500$ and mean relative error at (b) $Re=2000$, (d) $Re=6500$ for pressure and velocity

Like in the previous section, first we will present the results for strategy (31) and then proceed with strategy (32).

Constant artificial viscosity. Figure 17 shows the relative errors for the axial velocity for $\nu _a=0, 20, 50, 500, 1000$ in both the $Re=2000$ and $Re=6500$ case. Similarly to the $Re=3500$ case, the error decreases as $\nu _a$ increases. We observe that, upon reaching a certain value ($\nu _a=500$), higher values of artificial viscosity do not impact the ROM solution in both cases, with the errors becoming larger in the sudden expansion. We observe that the mean (over z) error in Fig. 17a is around $10^{-2}$ for the transitional case, whereas the mean error in Fig. 17b is one order of magnitude larger, roughly $10^{-1}$. It is somewhat expected that the effectiveness of a simple stabilization approach like (31) decreases as the Reynolds number is increased past a certain value.

Fig. 17

Relative errors of the ROM axial velocity with respect to the FOM solution for constant artificial viscosity $\nu _a = 0, 20, 50, 500, 1000$ in the case of (a) $Re=2000$ and (b) $Re=6500$

Figure 18 shows relative and absolute errors of the pressure for the same values of $\nu _a$ as in Fig. 17. As expected, the error of the pressure does not varies significantly as $\nu _a$ is increased. We note that the mean (over z) relative errors for the pressure are around $10^{-1}$ for both $Re=2000$ and $Re=6500$, with isolated peaks that are due to pressure values near zero. In fact, the absolute errors in correspondence to those peaks are small. See Fig. 18c and 18d. Therefore, we can conclude that overall the ROM solution has a satisfactory level of accuracy. This is confirmed by Fig. 19, which reports the time-averaged axial velocity and pressure along the z-axis for $\nu _a=1000$.

Fig. 18

Relative (top) and absolute (bottom) errors of the ROM pressure with respect to the FOM solution for constant artificial viscosity $\nu _a = 0, 20, 50, 500, 1000$

Fig. 19

Comparison of FOM and ROM time-averaged axial velocity (m/s) and pressure ($\hbox {m}^2$/$\hbox {s}^2$) with constant viscosity $\nu _a = 1000$ for (a) $Re=2000$ and (b) $Re=6500$

We conclude with a qualitative comparison of the FOM and ROM time average pressure and velocity in Figs. 20 and 21, respectively. We have used the same color bar for all Reynolds numbers ($Re = 2000, 3500, 6500$) to facilitate the comparisons. We see a good match between ROM and FOM solutions for both $Re = 2000$ and $Re = 6500$ when $\nu _a=1000$.

Fig. 20

Qualitative comparison of the ROM time-averaged axial pressure obtained with constant artificial viscosity $\nu _{a}=1000$ with the corresponding FOM solution for $Re=2000$ (left) and $Re=6500$ (right)

Fig. 21

Qualitative comparison of the ROM time-averaged axial velocity obtained with constant artificial viscosity $\nu _{a}=1000$ with the corresponding FOM solution for $Re=2000$ (left) and $Re=6500$ (right)

Modal coefficient-dependent artificial viscosity. Figure 22 shows the relative errors for the axial velocity for $\nu _a=0, 500, 1000$ when the coefficient-dependent stabilization approach (32) is used for the $Re=2000$ and $Re=6500$. The results closely resemble those of the constant viscosity approach in the sense that: i) the reduced solution that remains unchanged once a given value of artificial viscosity is attained ($\nu _a\ge 500$) and ii) the mean (over z) error in Fig. 22a is around $10^{-2}$ for the transitional case, whereas the mean error in Fig. 22b is one order of magnitude larger, roughly $10^{-1}$.

Fig. 22

Relative errors of the ROM axial velocity with respect to the FOM solution for the coefficient dependent viscosity $\nu _a = 0, 500, 1000$ in the case of (a) $Re=2000$ and (b) $Re=6500$

The absolute and relative errors of the pressure for $Re=2000$ and $Re=6500$ are shown in Fig. 23. A peak appears for $Re=2000$ at $z=0.17$ (see Fig. 23a) and two peaks appear for $Re=6500$ at $z=0.08, 0.17$ (see Fig. 23b). For $Re=6500$, we observed the same peaks when using the constant viscosity strategy: compare Fig. 23b with Fig. 18. As explained previously, these peaks are due to low values of the pressure. This is shown also by the plots of the absolute errors in Fig. 23c, d.

Fig. 23

Relative (top) and absolute (bottom) errors of the ROM pressure with respect to the FOM solution for coefficient-dependent viscosity with $\nu _a = 0, 500, 1000$

To conclude, Fig. 24 displays the FOM and ROM time-averaged axial velocity and pressure distributions along the z-axis when we set $\nu _a=1000$. Once again, we see a great match for the pressure. For $Re=2000$ there is also an excellent match for the axial velocity, while for $Re=6500$ the ROM axial velocity deviates from the FOM axial velocity in the expansion channel. This phenomenon gets worse with the coefficient dependent viscosity (compare Fig. 24b with Fig. 19b) and as Reynolds increases (compare Fig. 24b with Fig. 13f). For ease of comparison, we have reported the curves to be compared in Fig. 25: Fig. 25a compares the ROM time-averaged axial velocity using the linear and the coefficient-dependent approach for Re = 6500 with $\nu _a=1000$, while Fig. 25b illustrates the ROM time-averaged axial velocity for Re = 3500, 6500 with $\nu _a=1000$. A further increase of $\nu _a$ does not improve the match between FOM and ROM axial velocity for $Re=6500$.

Fig. 24

Comparison of FOM and ROM time-averaged axial velocity (m/s) and pressure ($\hbox {m}^2$/$\hbox {s}^2$) for (a) $Re=2000$ and (b) $Re=6500$ when using the coefficient-dependent viscosity with $\nu _a = 1000$

Fig. 25

Comparison of ROM time-averaged axial velocity (m/s) for (a) between linear (L) and coefficient-dependent (CD) viscosity at $Re=6500$ with $\nu _a=1000$ and (b) $Re=3500,6500$ with $\nu _a=1000$

4.3 Physical parameterization

In this section, we will consider the Reynolds number in the throat as a physical parameter varying in [3500, 6500]. Thus, our parameter domain in the Re-t plane is $[3500,6500] \times (0, 1.2)$ s. To sample this parameter domain, we have considered 9 equi-spaced values of Re ($\Delta Re$ = 375) and 60 equispaced snapshots in time (every 0.02 s) for each Re. The set $Re=\{3500, 3875, 4250, 4625, 5375, 5750, 6125, 6500\}$, with the all the corresponding snapshots in time, is used to generate the reduced basis, while $Re=5000$ is used for testing. So, the total number of snapshots employed in the POD is $8\times 60=480$. We remark that this number is a compromise choice between the accuracy of the ROM solution and the computational resources. In fact, it was not possible to perform a POD with 120 snapshots in time (used in Secs. 4.1 and 4.2) for each Re on the machine available to us (see last paragraph in Sect. 4.1).

The cumulative energy of the eigenvalues (17) is shown in Fig. 26a for velocity and pressure. While the pressure reaches $99.91\%$ of the energy with a single mode, the velocity needs 50 modes to achieve only $75.64\%$. The mean relative error (18) for pressure and velocity is shown in Fig. 26b as the number of modes increase. As expected, the error decreases very slowly with the number of modes: for example, when using 30 modes for each variable, we get around $80\%$ error for the velocity error and $40\%$ for the pressure.

Fig. 26

(a) Cumulative energy of the eigenvalues (17) and (b) mean relative error (18) for pressure and velocity when Re is a varying parameter

To contain the computational cost, we cut the number of modes to 12 for each variable. Figure 27 displays the reconstructed time-averaged axial velocity and pressure profiles along the z-axis at $Re = 5000$ if we use 12 modes and no stabilization. The large difference between the FOM and ROM curves is expected given the significant errors depicted in Fig. 26b.

Fig. 27

Comparison of FOM and ROM time-averaged axial velocity (m/s) and pressure ($\hbox {m}^2$/$\hbox {s}^2$) at $Re = 5000$ (test value) using 12 modes for each variable and no stabilization technique

Like in the previous sections, we will present first the results obtained with method (31) and then the results of approach (32).

Constant artificial viscosity. Figure 28 displays the time-averaged axial velocity and pressure along the z-axis at $Re = 5000$ for different values of $\nu _a$. We recall that the difference with the previous subsections is that the FOM solutions for $Re = 5000$ were not used to generate the reduced basis. Nonetheless, the observations are similar to the ones made for the other Reynolds numbers. For low $\nu _a$ (i.e., $\nu _a=1.5$), the velocity is underestimated in the sudden expansion and it exhibits oscillations. The ROM reconstruction improves considerably for both pressure and velocity when $\nu _a=3$, but the fluctuations of the velocity curve in the sudden expansion are still present. For higher values of $\nu _a$, oscillations disappear but this method is not able to recover properly the velocity at the sudden expansion and the pressure in the entrance channel. Increasing the artificial viscosity above $\nu _a=1000$ does not influence the results.

Fig. 28

Comparison of FOM and ROM time-averaged axial velocity (m/s) and pressure ($\hbox {m}^2$/$\hbox {s}^2$) for $Re = 5000$ (test value) with constant artificial viscosity for (a) $\nu _a=1.5$, (b) $\nu _a=3$, (c) $\nu _a=20$, and (d) $\nu _a=1000$

Relative and absolute errors along the z-axis for time-averaged pressure and velocity are shown in Fig. 29. Overall, the axial velocity error decreases as $\nu _a$ increases, with the relative error remaining below $10^{-1}$ along most of the z-axis. The relative error for the pressure oscillates around $10^{-1}$.

Fig. 29

Relative (top row) and absolute (bottom row) errors of the ROM velocity (left) and pressure (right) with respect to the FOM solution for $Re = 5000$ (test value) when using the constant viscosity approach and $\nu _a = 0, 100, 1000$

A qualitative comparison of the FOM and ROM time average pressure and velocity for the test point $Re=5000$ is shown in Fig. 30. Overall, we see a good reconstruction of axial velocity and pressure in the entire domain.

Fig. 30

Qualitative comparison of the ROM time-averaged axial pressure (left) and velocity (right) obtained with constant artificial viscosity $\nu _{a}=1000$ with the corresponding FOM solution for $Re = 5000$ (test value)

Modal coefficient-dependent artificial viscosity. The time-averaged axial velocity and pressure along the z-axis obtained with stabilization technique (32) are depicted in Fig. 31, for $\nu _a=3, 10, 50, 1000$. For low values of $\nu _a$ ($\nu _a=3,10$) we see oscillations in the sudden expansion, like with the constant viscosity approach, but unlike it we have an overestimation of the velocity in the throat. With $\nu _a=50$, we note a significantly higher ROM velocity compared to the FOM solution in the sudden expansion. Finally, with $\nu _a=1000$ we obtain the same reconstruction as with the constant viscosity strategy. Increasing $\nu _a$ does not lead to any improvement in the results.

Fig. 31

Comparison of FOM and ROM time-averaged axial velocity (m/s) and pressure ($\hbox {m}^2$/$\hbox {s}^2$) at $Re = 5000$ (test value) with coefficient-dependent artificial viscosity for (a) $\nu _a=3$, (b) $\nu _a=10$, (c) $\nu _a=50$, and (d) $\nu _a=1000$

Figure 32 shows the relative and absolute error along the z-axis for time-averaged pressure and velocity for increased values of $\nu _a$. The trend and magnitude of the errors are similar to the case with constant artificial viscosity, although we can observe some differences. In fact, while with the constant artificial viscosity method changing $\nu _a$ from 100 to 1000 made little difference in the errors, with the coefficient-dependent artificial viscosity the difference in the errors is more significant.

Fig. 32

Relative (top row) and absolute (bottom row) errors of the ROM velocity (left) and pressure (right) with respect to the FOM solution at $Re = 5000$ (test value) for $\nu _a = 0, 100, 1000$ in the coefficient-dependent viscosity approach

Finally, we would like to comment on the overestimation of the ROM axial velocity in the sudden expansion. While for $Re = 2000, 3500$ it can be minimized by increasing $\nu _a$, we saw in Sect. 4.2 and in this section that increasing $\nu _a$ does not work for $Re=5000, 6500$. We experimented with augmenting the training dataset by including $Re=4813,5188$, i.e., we refined the sampling near the test point, but that did not help. We suspect that a much finer sampling of the parameter domain is needed to lower the ROM axial velocity in the the sudden expansion. However, that comes with considerable computational challenges.

5 Conclusions

We presented two POD-Galerkin reduced order models that stabilize the ROM solution with artificial diffusion. The first ROM introduces a constant artificial viscosity, while the second employs a mode-dependent artificial viscosity. We evaluated the accuracy of both ROMs in the simulation of flow in a nozzle spanning the transitional and turbulent regimes (from $Re=2000$ to $Re=6500$).

First, we considered time as the only parameter. Our numerical results showed that both stabilized ROMs accurately reproduce time-averaged FOM solutions. We obtained a high-accuracy ROM reconstruction for $Re=2000,3500$ when the added artificial viscosity is sufficiently large. In addition, we showed that for such large values, the constant and coefficient-dependent viscosity ROMs become equivalent and both approaches can be used without a notable impact on the accuracy of the solution. The relative error for the time-averaged axial velocity is around $10^{-2}$, while the relative error for the pressure is approximately $10^{-1}$. For $Re=6500$, the ROMs solution becomes less accurate in the sudden expansion region and the relative error for the time-averaged axial velocity increases to $10^{-1}$.

Next, we added the Reynolds number to the parameter space. We observed that the accuracy of the ROM solution depends on the Reynolds number: as the Reynolds number gets higher, the ROM solution becomes progressively less accurate in the sudden expansion, also seen when time is the only parameter. However, for all Re, the relative error for axial velocity ranges between $10^{-1}$ and $10^{-2}$, while the relative error for pressure is around $10^{-1}$.

The main limitation of these stabilized ROMs appears when the added artificial viscosity is large enough to ensure good accuracy: the ROM solution becomes steady. We plan to overcome the limitation through the use of more sophisticated stabilization techniques based on elliptic filters.