A segregated reduced order model of a pressure-based solver for turbulent compressible flows
]Matteo<EMAIL_ADDRESS>]Giovanni<EMAIL_ADDRESS>]Gianluigi<EMAIL_ADDRESS>Mathematics Area, mathLab, SISSA, via Bonomea 265, I-34136 Trieste,

\begin@twocolumnfalse

This article provides a reduced order modelling framework for turbulent compressible flows discretized by the use of finite volume approaches. The basic idea behind this work is the construction of reduced order models capable of providing fully consistent solutions with respect to the high fidelity flow fields. Full order solutions are often obtained through the use of segregated solvers, employing slightly modified conservation laws so that they can be decoupled and then solved one at a time. Classical reduction architectures, on the contrary, rely on the Galerkin projection of a complete Navier-Stokes system to be projected all at once, causing a mild discrepancy with the high order solutions. In this article we rely on segregated reduced order algorithms for the resolution of turbulent and compressible flows in the context of physical and geometrical parameters. At the full order level turbulence is modeled using an eddy viscosity approach. Since there are a variety of different turbulence models for the approximation of this supplementary viscosity, one of the aims of this work is to provide reduced order models which are independent on this selection. This goal is reached by the application of hybrid methods where Navier-Stokes equations are projected in a standard way while the viscosity field is approximated by the use of data-driven interpolation methods or by the evaluation of a properly trained neural network. By exploiting the aforementioned expedients it is possible to reduced the computational cost associated wit fluid flow problems characterized by high Reynolds numbers and elevated Mach numbers.

\end@twocolumnfalse

§ INTRODUCTION

In the last decades fluid flows simulations have progressively enlarged their applicability and their influence in many different research fields (general overviews can be found in [1, 2, 3]). Nowadays Computational Fluid Dynamics (CFD) applications have reached widely spread ambits as, for example, shape optimization for naval/automotive/aerospace engineering ([4, 5]), cardiovascular in real time surgery ([6]), chemistry industrial processes ([7, 8]) or weather forecasts ([9]). While increasing the demand on reliability and usability of CFD methods, the computational capability of the employed hardware architectures are no more sufficient in terms of time consuming. For this reason the search for new efficient methods able to reduce computational time keeps on covering a relevant amount of CFD research efforts.

A popular research field, related to the aforementioned applications, is the analysis of the dynamics relative to parametrized Partial Differential Equations (PDEs). In this case an infinite number of solutions is available for every slightly different value of the parameter. For some specific ambits, many of them have to be evaluated in order to find out the one that is best performing in terms of prefixed requirements fulfilling. The cost of such a procedure can easily become unaffordable, in particular when the dimension of the problem is big enough (see, e.g., [10]). Typical applications where such a difficulty is often arising are, for example, shape optimization problems, uncertainty quantification studies or optimal control frameworks.

Recently many different techniques have been taken into consideration to try to overtake this problem. Galerkin projection has widely been employed to develop new reduction strategies capable of exploiting the information of just a few full order solutions for different parameter values in order to perform efficient, accurate and much cheaper solutions for a different selection of the parameter (for fluid flows applications, [11, 12, 13, 14, 15] represents relevant works). Many different choices are possible to take advantage of the dynamical content contained into high fidelity solutions. The most used ones are the Proper Orthogonal Decomposition (POD) (see [16, 17, 18, 19, 20]), the Proper Generalized Decomposition (PGD) (see [21, 22]), the Dynamic Mode Decomposition (DMD) (see [23, 24]) or the Greedy algorithm (see [25, 26]).

Lately a new research branch has risen in this research sector: Machine Learning. By the use of neural networks the resolution of parametric PDEs has surely become easier and more accessible ([27, 28, 29, 30, 31, 32]).

Both the two mentioned approaches have some valuable aspects together with critical points to be underlined. Projection techniques are strongly connected with physical laws of the problem since they use modal basis functions obtained by real solutions to extract the main dynamics and they employ those modes to project and reconstruct conservation equations solutions manifolds. Unfortunately non-linearity and non-affinity of the parametrized formulation can be difficult issues to be carefully treated. Moreover sometimes it may happen that the equations are not directly available and in that case these methods are not employable. A classical example is constituted by commercial software where a deep description of the employed laws is not provided.  Again, in the literature one may find many examples of very slow decay of the eigenvalues for specific problems, that means many modal basis functions are needed to reach a decent approximation for the solutions, wasting, de facto, the whole possible gain coming from such an architecture ([33]). A very well known problem affected by such an issue is represented by turbulent flows (see e.g. [34, 35, 36, 37, 38]). On the contrary Artificial Intelligence (AI) techniques are very versatile. They just need a set of solutions to be trained and, despite the complexity of the mathematical formulation of the starting problem, they can be constructed and modeled to pursue good approximations in a short time. The dark side of these approaches is the fact that they have a much weaker connection with the real physics of what they are approximating and the actual meaning of every single part constituting their architecture is arguably comprehensible in terms of phenomena representation. For this reason they may give inaccurate results thanks to impossibility in having a deeper check on networks responses.

Taking all the aforementioned examinations under consideration, in this work we provide a new mixed technique for Navier Stokes compressible problems, capable of merging the advantages of projection techniques together with neural networks data driven architectures. In particular, in our approach, classical projection methods are used for the Favre Averaged Navier Stokes (FANS) equations while a neural network gets trained to provide the eddy viscosity solutions in a turbulence modeling approach. The main goal is to propose an architecture proficient in dealing with different types of parametrizations for compressible flows. Moreover one of the most relevant focuses concerning this work is constituted by a coherent approach between full order and reduced order solutions, by developing a new reduced compressible SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm.

This work will be divided into six different sections: <ref>, <ref> and <ref> present the equations we will use for this work and their Finite Volume discretization; <ref> explains the Proper Orthogonal Decomposition (POD) procedure employed to obtain the modal basis functions. In <ref> the core algorithm used for our technique is introduced together with <ref> where the AI architecture for turbulence treatment is shown. Two different test cases, a physically parametrized and a geometrically parametrized ones, are exposed in <ref> and <ref> respectively. Finally into <ref> few considerations on the results and some possible developments for this work are present.

§ THE COMPRESSIBLE NAVIER-STOKES EQUATIONS

In this work we want to deal with parametrized compressible Navier-Stokes equations poblems. To manage the compressibility of the fluid, we selected a common strategy for this kind of applications: the Favre averaging. The equations describing the physics are the following ones:

\begin{equation} \label{eq:FANS}
\begin{cases}
\displaystyle \frac{\partial \overline{\rho}}{\partial t} + \nabla \cdot \left( \overline{\rho} \tilde{\bm{u}}\right)= 0 & \text{in} \ \Omega(\pi),\\
\displaystyle \frac{\partial \overline{\rho} \tilde{\bm{u}}}{\partial t} + \nabla \cdot \left[ \overline{\rho} \tilde{\bm{u}} \otimes \tilde{\bm{u}} - \tilde{\bm{\tau}}_{turb} - \tilde{\bm{\tau}} + \overline{p} \bm{I}\right] = 0 & \text{in} \ \Omega(\pi),\\
\displaystyle \frac{\partial \overline{\rho} \tilde{e}_0}{\partial t} + \nabla \cdot \bigg[ \overline{\rho} \tilde{\bm{u}} \tilde{e}_0 - C_p \frac{\mu}{Pr}\nabla \tilde{T} - C_p \frac{\mu_t}{Pr}_{t}\nabla \tilde{T} \big. \\
\big. + \overline{p} \tilde{\bm{u}} - \tilde{\bm{u}} \cdot \tilde{\bm{\tau}}  - \tilde{\bm{u}} \cdot \tilde{\bm{\tau}}_{turb} \bigg] = 0 & \text{in} \ \Omega(\pi),\\
\tilde{\bm{u}} = \bm{g}_D & \text{in} \ \Gamma_D,\\
\nu \displaystyle \frac{\partial \tilde{\bm{u}}}{\partial \bm{n}} - \overline{p} \bm{n} = g_N &\text{in} \ \Gamma_{N},
\end{cases}
\end{equation}

where $\overline{p} \in \mathbb{Q}$, $\tilde{\bm{u}} \in \mathbb{V}$ and $\tilde{e} \in \mathbb{E}$ are the unknowns of the problem, $\Omega(\pi)$ is the domain of definition and it can be, in geometrical parametrization cases, directly dependent on the parameter $\pi$, $\Gamma_D$ stands for the part of the boundary where the Dirichled condition $\bm{g}_D$ is imposed while $\Gamma_N$ is the part of the boundary where the Neumann condition $\bm{g}_N$ is imposed, $C_p$ means specific heat at constant pressure, $\bm{\tau}$ indicates the viscous stress tensor, $\bm{\tau}_{turb}$ stands for the extra viscosity term due to turbulence, $\mu$ is the dynamic viscosity, $\mu_{turb}$ is the extra viscosity owing to turbulence, $\nu$ is the kinematic viscosity, $Pr$ indicates the Prandtl number and $Pr_t$ its turbulent counterpart which is a constant value.

The set of equations above have to be accompanied by an equation of state defining the relation existing between density $\rho$, pressure $p$ and temperature $T$: $\rho = \rho(p,T)$. In our case we suppose air to be an ideal gas, leading to

$$\rho = \frac{p}{RT},$$

being $R$ the gas constant.

In the Favre Averaged Navier-Stokes (FANS) equations, all the variables (density $\rho$, pressure $p$, velocity $\bm{u}$, total energy $e_0$, temperature $T$ and internal energy $e$) are decomposed into an averaged part and a fluctuating one as follows:

$$\rho = \overline{\rho} + \rho', \hspace{0.5cm} p = \overline{p} + p', \hspace{0.5cm} \bm{u} = \tilde{\bm{u}} + \bm{u}'',$$

$$e_0 = \tilde{e_0} + e_0'', \hspace{0.5cm} T = \tilde{T} + T'', \hspace{0.5cm} e = \tilde{e} + e'' .$$

Superscript $\tilde{\square}$ indicates the Favre averaging which correspond to a density weighted Reynolds averaging $\overline{\square}$. Given a certain variable $\Phi(t)$, we have:

$$\overline{\Phi} = \frac{1}{T} \int_T \Phi(t) dt \Rightarrow \Phi' = \Phi - \overline{\Phi},$$

$$\tilde{\Phi} = \frac{\overline{\rho \phi}}{\overline{\rho}} \Rightarrow \Phi'' = \Phi - \tilde{\Phi} .$$

<ref> is obtained after some approximations and assumptions from an eddy viscosity point of view. The reader interested on the averaging procedure and modeling should refer to [39].

From now on, <ref> will be considered only into its steady-state formulation.
All the averaged variables are dependent on the parameter $\pi$ but they will be indicated without any explicit reference to it for sake of simplicity:

$$\overline{\rho} = \overline{\rho}(t), \hspace{0.5cm} \overline{p} = \overline{p}(t),$$

$$\tilde{\bm{u}} = \tilde{\bm{u}}(t), \hspace{0.5cm} \tilde{T} = \tilde{T}(t), \hspace{0.5cm} \tilde{e} = \tilde{e}(t).$$

Moreover, the energy equation will be be rewritten into the internal energy form, leading to:

\begin{equation*}
\begin{cases}
\nabla \cdot \left( \overline{\rho} \tilde{\bm{u}}\right)= 0 & \text{in} \ \Omega(\pi),\\
\nabla \cdot \left[ \overline{\rho} \tilde{\bm{u}} \otimes \tilde{\bm{u}} - \tilde{\bm{\tau}}_{turb} - \tilde{\bm{\tau}} + \overline{p} \bm{I}\right] = 0 & \text{in} \ \Omega(\pi),\\
\displaystyle \nabla \cdot \Bigg[ \overline{\rho} \tilde{\bm{u}} \left( \tilde{e} + \frac{\tilde{\bm{u}} \cdot \tilde{\bm{u}}}{2} + \underbrace{\widetilde{\frac{\bm{u}'' \cdot \bm{u}''}{2}}}_{\star}\right) - \frac{C_p}{C_v}*\\
\left( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \right) \nabla \tilde{e}
+ \overline{p} \tilde{\bm{u}} - \underbrace{\tilde{\bm{u}} \cdot \tilde{\bm{\tau}}}_{\star}  - \underbrace{\tilde{\bm{u}} \cdot \tilde{\bm{\tau}}_{turb}}_{\star} \Bigg] = 0 & \text{in} \ \Omega(\pi),\\
\tilde{\bm{u}} = \bm{g}_D & \text{in} \ \Gamma_D,\\
\nu \displaystyle \frac{\partial \tilde{\bm{u}}}{\partial \bm{n}} - \overline{p} \bm{n} = g_N &\text{in} \ \Gamma_{N}.
\end{cases}
\end{equation*}

Star terms can be reasonably neglected if compared with the other terms present into the energy equation, resulting in the system:

\begin{equation*}
\begin{cases}
\nabla \cdot \left( \overline{\rho} \tilde{\bm{u}}\right)= 0 & \text{in} \ \Omega(\pi),\\
\nabla \cdot \left[ \overline{\rho} \tilde{\bm{u}} \otimes \tilde{\bm{u}} - \tilde{\bm{\tau}}_{turb} - \tilde{\bm{\tau}} + \overline{p} \bm{I}\right] = 0 & \text{in} \ \Omega(\pi),\\
\displaystyle \nabla \cdot \bigg[ \overline{\rho} \tilde{\bm{u}} \left( \tilde{e} + \frac{\tilde{\bm{u}} \cdot \tilde{\bm{u}}}{2}\right) - \frac{C_p}{C_v}*\\
\left( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \right) \nabla \tilde{e} + \overline{p} \tilde{\bm{u}} \bigg] = 0 & \text{in} \ \Omega(\pi),\\
\tilde{\bm{u}} = \bm{g}_D & \text{in} \ \Gamma_D,\\
\nu \displaystyle \frac{\partial \tilde{\bm{u}}}{\partial \bm{n}} - \overline{p} \bm{n} = g_N &\text{in} \ \Gamma_{N}.
\end{cases}
\end{equation*}

We can write explicitly also the viscous stress terms into the momentum equation:

$$\tilde{\bm{\tau}} = 2\mu \tilde{\bm{S}}, \hspace{0.5cm} \tilde{\bm{\tau}}_{turb} = 2\mu_t \tilde{\bm{S}},$$

where $\tilde{\bm{S}} = \frac{\nabla \tilde{\bm{u}} + \nabla \tilde{\bm{u}}^T}{2} - \frac{1}{3} \nabla \cdot \tilde{\bm{u}} \bm{I}$. Also the momentum equation can then be rewritten exploiting this last definition:

\begin{equation} \label{eq:rhosimple}
\begin{cases}
\nabla \cdot \left( \overline{\rho} \tilde{\bm{u}}\right)= 0 & \text{in} \ \Omega(\pi),\\
\displaystyle \nabla \cdot \bigg[ \overline{\rho} \tilde{\bm{u}} \otimes \tilde{\bm{u}} - \big( \mu + \mu_t \big)\bigg( \nabla \tilde{\bm{u}} + \nabla \tilde{\bm{u}}^T\\
- \frac{2}{3} \nabla \cdot \tilde{\bm{u}} \bm{I}\bigg) + \overline{p} \bm{I}\bigg] = 0 & \text{in} \ \Omega(\pi),\\
\displaystyle \nabla \cdot \bigg[ \overline{\rho} \tilde{\bm{u}} \bigg( \tilde{e} + \frac{\tilde{\bm{u}} \cdot \tilde{\bm{u}}}{2}\bigg) - \frac{C_p}{C_v} *\\
\bigg( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \bigg) \nabla \tilde{e} + \overline{p} \tilde{\bm{u}} \bigg] = 0 & \text{in} \ \Omega(\pi),\\
\tilde{\bm{u}} = \bm{g}_D & \text{in} \ \Gamma_D,\\
\nu \displaystyle \frac{\partial \tilde{\bm{u}}}{\partial \bm{n}} - \overline{p} \bm{n} = g_N &\text{in} \ \Gamma_{N}.
\end{cases}
\end{equation}

It is now clear that all the turbulence-related terms of the equations rely on $\mu_t$ to be calculated. For this reason, since only the eddy viscosity is required, a common 2-equations turbulent model as, e.g., $k-\epsilon$ or $k-\omega$, is sufficient as a closure for the problem.

§ FINITE VOLUME DISCRETIZATION

The first step towards a Finite Volume discretization (for a deeper insight see [40]) of the problem is the division of the domain $\Omega(\pi)$ into a tessellation $\mathcal{T}(\pi)$ composed by a certain number $N_h$ of cells $\Omega_i(\pi)$, so that:

\begin{equation*}
\mathcal{T}(\pi) = \{ \Omega_i(\pi) \}_{i=1}^{N_h}, \hspace{1cm} \bigcup_{i=1}^{N_h} \Omega_i (\pi) = \Omega (\pi),
\end{equation*}

where every cell $\Omega_i$ can be constructed as a non-convex polyhedron.

The Finite Volume variables can be here introduced: $\overline{p}_h \in \mathbb{Q}_h$, $\tilde{\bm{u}}_h \in \mathbb{Q}_h$ and $\tilde{e}_h \in \mathbb{E}_h$. They are not continuous and they are constant in the interior part of each cell assuming everywhere the value at the center of the cell. For sake of simplicity in this section we will keep on referring to this variables without the $\square_h$ subscript to not make the formulas too heavy.

(-1,-0.5) – (0,0) – (-0.5,1) – (-1.5,1.5);
(-0.5,1) – (0.5,1.7) – (0.5,2.5);
(0.5,1.7) – (2,1.5) – (1,0) – (1,-1);
(0,0) – (1,0) – (3,-0.5) – (4,-1);
(3,-0.5) – (4,1.5) – (2,1.5) – (2.5,2.5);
(4,1.5) – (5.5,2);
[dashed]  (0.65,0.8) – (2.5,0.6);
[thick,->] (1.5,0.7) – (2.5,0.25);
[blue!] (0.65,0.8) circle (0.1cm);
[blue!] (2.5,0.6) circle (0.1cm);
at (0.65,0.5) $v_i$;
at (2.5,0.9) $v_j$;
at (2.5,-0.1) $\bm{S}_f$;
at (0.65,1.4) $\Omega_i$;
at (3.5,1.2) $\Omega_j$;
at (1.2,1) $\bm{d}_f$;

Relation between two neighbor cells of the tessellation $\mathcal{T}$ for a certain variable $v$.

We use in this work a segregated approach based on a compressible formulation of the SIMPLE algorithm. This aspect has to be kept in mind for the Finite Volume discretization strategy.

We can now take into consideration the equations one by one, starting from the continuity constrain. A detailed treatment for what concerns fluid flows equations discretization, the interested reader can refer to [41]. The equation can be written into its integral form over each cell as follows:

\begin{equation*}
\int_{\Omega_i} \nabla \cdot \left( \overline{\rho} \tilde{\bm{u}} \right) dV = 0.
\end{equation*}

By exploiting the divergence theorem, the equation above leads to:

\begin{equation*}
\int_{\delta \Omega_i} \overline{\rho} \tilde{\bm{u}} \cdot d\bm{S} \simeq \sum_{f=1}^{Nf_i} \left( \overline{\rho} \tilde{\bm{u}} \right)\rvert_{f} \cdot \bm{S}_f =  \sum_{f=1}^{Nf_i} F_f,
\end{equation*}

where the subscript $\square\rvert_{f}$ indicates that those variables are evaluated at the center of the face $f$ and $\bm{S}_f$ is the oriented surface of the same face while $Nf_i$ if the total number of faces surrounding the $i$-th cell while $F_f$ stands for the mass flux crossing the face $f$ (see <ref>).

The discretized version of the continuity equation then reads:

\begin{equation}\label{eq:contEq}
\sum_{f=1}^{Nf_i} F_f = 0.
\end{equation}

Of course all the variables in a Finite Volume scheme are known only at the cell center so their values at the center of the faces has to be calculated by interpolating the neighbor cell center values.

Let us now take into consideration the momentum equation. It has to be integrated over the volume of every cell and it can then be analysed term by term, starting from the convective one:

\begin{multline*}
\int_{\Omega_i} \bm{\nabla} \cdot (\overline{\rho}\tilde{\bm{u}} \otimes \tilde{\bm{u}}) \, dV = \int_{\delta \Omega_i} \overline{\rho} \tilde{\bm{u}} \otimes \tilde{\bm{u}} \cdot d\bm{S} \simeq\\
\simeq \sum_{f=1}^{Nf_i} \bm{S}_f \cdot \left( \overline{\rho} \tilde{\bm{u}} \otimes \tilde{\bm{u}}\right)\rvert_f = \sum_{f=1}^{Nf_i} \bm{F}_f \tilde{\bm{u}}_f.
\end{multline*}

The first part of the diffusion term is discretized as follows:

\begin{multline*}
\int_{\Omega_i} \bm{\nabla} \cdot \left[ \left(\mu + \mu_t \right) \nabla \tilde{\bm{u}} \right] \, dV = \int_{\delta \Omega_i} \left[ \left(\mu + \mu_t \right) \nabla \tilde{\bm{u}} \right] \cdot d\bm{S} \simeq \\
\simeq \sum_{f=1}^{Nf_i} \left[ \left(\mu + \mu_t \right) \nabla \tilde{\bm{u}} \right] \rvert_f \cdot \bm{S}_f.
\end{multline*}

For orthogonal meshes we can approximate this term as:

\begin{equation*}
\sum_{f=1}^{Nf_i} \left[ \left(\mu + \mu_t \right) \nabla \tilde{\bm{u}} \right] \rvert_f \cdot \bm{S}_f \simeq \left(\mu + \mu_t \right)\rvert_f \ \lvert \bm{S}_f \rvert \frac{\tilde{\bm{u}}_i - \tilde{\bm{u}}_j}{\lvert \bm{d}_f \rvert},
\end{equation*}

being $\bm{d}_f$ the oriented vector bridging the cell centers of two neighbor cells. When this is  not the case, a non-orthogonal correction is added:

\begin{multline*}
\sum_{f=1}^{Nf_i} \left[ \left(\mu + \mu_t \right) \nabla \tilde{\bm{u}} \right] \rvert_f \cdot \bm{S}_f \simeq \sum_{f=1}^{Nf_i} \left(\mu + \mu_t \right)\rvert_f *\\
\left[ \lvert \bm{\mathcal{P}}_f \rvert \frac{\tilde{\bm{u}}_i - \tilde{\bm{u}}_j}{\lvert \bm{d} \rvert} + \bm{\mathcal{O}}_f \cdot \nabla \tilde{\bm{u}}\rvert_f \right],
\end{multline*}

where we have $\bm{\mathcal{P}}_f \parallel \bm{d}_f$, $\bm{\mathcal{O}}_f \bot \bm{d}_f$ and $\bm{\mathcal{P}}_f + \bm{\mathcal{O}}_f = \bm{S}_f$ while $\nabla \tilde{\bm{u}}\rvert_f $ is evaluated starting from its value at the cell centers $\nabla \tilde{\bm{u}}_i$ and $\nabla \tilde{\bm{u}}_j$ by interpolation.

The second part of the diffusion term is treated following the previous steps:

\begin{multline*}
\int_{\Omega_i} \bm{\nabla} \cdot \left[ \left(\mu + \mu_t \right) \nabla \tilde{\bm{u}}^T \right] \, dV = \int_{\delta \Omega_i} \left[ \left(\mu + \mu_t \right) \nabla \tilde{\bm{u}}^T \right] \cdot d\bm{S} \simeq\\
\simeq \sum_{f=1}^{Nf_i} \left[ \left(\mu + \mu_t \right) \nabla \tilde{\bm{u}}^T \right] \rvert_f \cdot \bm{S}_f.
\end{multline*}

In this case, the face center evaluation is treated explicitly so that this term is considered to be a forcing term:

\begin{multline*}
\nabla \tilde{\bm{u}}^T \rvert_f \cdot \bm{S}_f = \begin{bmatrix}
\frac{\partial \tilde{u}_x}{\partial x} & \frac{\partial \tilde{u}_x}{\partial y} & \frac{\partial \tilde{u}_x}{\partial z}\\
\frac{\partial \tilde{u}_y}{\partial x} & \frac{\partial \tilde{u}_y}{\partial y} & \frac{\partial \tilde{u}_y}{\partial z}\\
\frac{\partial \tilde{u}_z}{\partial x} & \frac{\partial \tilde{u}_z}{\partial y} & \frac{\partial \tilde{u}_z}{\partial z}\\
\end{bmatrix}^T_f \cdot \begin{bmatrix}
\end{bmatrix}_f = \\
= \begin{bmatrix}
\frac{\partial \tilde{u}_x}{\partial x} S_x + \frac{\partial \tilde{u}_y}{\partial x} S_y + \frac{\partial \tilde{u}_z}{\partial x} S_z\\
\frac{\partial \tilde{u}_x}{\partial y} S_x + \frac{\partial \tilde{u}_y}{\partial y} S_y + \frac{\partial \tilde{u}_z}{\partial y} S_z\\
\frac{\partial \tilde{u}_x}{\partial z} S_x + \frac{\partial \tilde{u}_y}{\partial z} S_y + \frac{\partial \tilde{u}_z}{\partial z} S_z\\
\end{bmatrix}.
\end{multline*}

The same applies for the last part of the diffusive term:

\begin{multline*}
\int_{\Omega_i} \bm{\nabla} \left[ \left(\mu + \mu_t \right) \frac{2}{3} \nabla \cdot \tilde{\bm{u}} \right] \, dV = \\
= \int_{\delta \Omega_i} \left[ \left(\mu + \mu_t \right) \frac{2}{3} \nabla \cdot \tilde{\bm{u}} \right] d\bm{S} \simeq \sum_{f=1}^{Nf_i} \left[ \left(\mu + \mu_t \right) \frac{2}{3} \nabla \cdot \tilde{\bm{u}} \right] \Big\rvert_f \bm{S}_f,
\end{multline*}

where once again the divergence of the velocity is interpolated to the surface and treated explicitly leading to an additional forcing term.

The last term to be considered is the pressure gradient:

\begin{equation*}
\int_{\Omega_i} \bm{\nabla} \overline{p} \, dV = \int_{\delta \Omega_i} \overline{p} \ d\bm{S} \simeq \sum_{f=1}^{Nf_i} \overline{p}_f \bm{S}_f,
\end{equation*}

In the momentum equation, pressure is interpolated to the faces and then treated explicitly as a source term. The final momentum equation reads:

\begin{equation}\label{eq:momEq}
\begin{aligned}
\sum_{f=1}^{Nf_i} \left[ \bm{F}_f \tilde{\bm{u}}_f - \left(\mu + \mu_t \right)\rvert_f \left( \lvert \bm{\mathcal{P}}_f \rvert \frac{\tilde{\bm{u}}_i - \tilde{\bm{u}}_j}{\lvert \bm{d} \rvert} + \bm{\mathcal{O}}_f \cdot \nabla \tilde{\bm{u}}\rvert_f \right) \right]=\\ =\sum_{f=1}^{Nf_i} \left[ \left(\mu + \mu_t \right)\rvert_f \left( \nabla \tilde{\bm{u}}^T_f \cdot \bm{S}_f - \frac{2}{3} \nabla \cdot \tilde{\bm{u}}_f \bm{S}_f\right) - \overline{p}_f \bm{S}_f \right],
\end{aligned}
\end{equation}

where all the terms composing the right-hand side of the equation are treated explicitly as source terms.

<ref> can be rewritten into its Finite Volume matrix form as follows:

\begin{equation}\label{eq:mMomEq}
\bm{A}_u(\tilde{\bm{u}}) = -\nabla \overline{p} \Rightarrow \bm{A} \tilde{\bm{u}} = \bm{H} (\tilde{\bm{u}}) -\nabla \overline{p},
\end{equation}

where $\bm{A}_u(\tilde{\bm{u}})$ is the Finite Volume discretized form containing all the terms related to velocity of both left-hand and right-hand sides of <ref>, $\bm{A}\tilde{\bm{u}}$ is the diagonal part of $\bm{A}_u(\tilde{\bm{u}})$ while $-\bm{H} (\overline{\bm{u}})$ is its extra diagonal part so that $\bm{A}_u(\tilde{\bm{u}}) = \bm{A}\tilde{\bm{u}} - \bm{H} (\tilde{\bm{u}})$.

The last equation to be analysed regards the energy conservation:

\begin{multline*}
\int_{\Omega_i} \nabla \cdot \left[ \overline{\rho} \tilde{\bm{u}} \tilde{e} \right] \, dV = \int_{\delta \Omega_i} \overline{\rho} \tilde{\bm{u}} \tilde{e} \cdot d\bm{S} \simeq \\
\simeq\sum_{f=1}^{Nf_i} \tilde{e}_f \overline{\rho}_f \tilde{\bm{u}}_f \cdot \bm{S}_f = \sum_{f=1}^{Nf_i} \tilde{e}_f F_f.
\end{multline*}

The kinetic part of the total energy is treated explicitly and leads to:

\begin{multline*}
\int_{\Omega_i} \nabla \cdot \left[ \overline{\rho} \tilde{\bm{u}} \frac{\tilde{\bm{u}} \cdot \tilde{\bm{u}}}{2} \right] \, dV = \int_{\delta \Omega_i} \overline{\rho} \tilde{\bm{u}} \frac{\tilde{\bm{u}} \cdot \tilde{\bm{u}}}{2} \cdot d\bm{S} \simeq \\
\simeq \sum_{f=1}^{Nf_i} \frac{\tilde{\bm{u}}_f \cdot \tilde{\bm{u}}_f}{2} \overline{\rho}_f \tilde{\bm{u}}_f \cdot \bm{S}_f = \sum_{f=1}^{Nf_i} \frac{\tilde{\bm{u}}_f \cdot \tilde{\bm{u}}_f}{2} F_f.
\end{multline*}

The diffusive term reads:

\begin{multline*}
\int_{\Omega_i} \nabla \cdot \left[ \frac{C_p}{C_v} \left( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \right) \nabla \tilde{e} \right] \, dV = \\
= \int_{\delta \Omega_i} \frac{C_p}{C_v} \left( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \right) \nabla \tilde{e} \cdot d\bm{S} \simeq \\
\simeq \sum_{f=1}^{Nf_i} \frac{C_p}{C_v} \left( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \right)\Big \rvert_f \nabla \tilde{e}_f \cdot \bm{S}_f.
\end{multline*}

Once again the energy gradient is not available at the center of the faces but it can be approximated:

\begin{multline*}
\sum_{f=1}^{Nf_i} \frac{C_p}{C_v} \left( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \right)\Big \rvert_f \nabla \tilde{e}_f \cdot \bm{S}_f \simeq \\
\simeq \sum_{f=1}^{Nf_i} \frac{C_p}{C_v} \left( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \right)\Big \rvert_f \left[ \lvert \bm{\mathcal{P}}_f \rvert \frac{\tilde{e}_i - \tilde{e}_j}{\lvert \bm{d} \rvert} + \bm{\mathcal{O}}_f \cdot \nabla \tilde{e}\rvert_f \right].
\end{multline*}

Finally also the pressure term is discretized and treated explicitly:

\begin{equation*}
\int_{\Omega_i} \nabla \cdot \left[ \overline{p} \tilde{\bm{u}} \right] \, dV = \int_{\delta \Omega_i} \overline{p} \tilde{\bm{u}} \cdot d\bm{S} \simeq \sum_{f=1}^{Nf_i} \overline{p}_f \tilde{\bm{u}}_f \cdot \bm{S}_f = \sum_{f=1}^{Nf_i} \frac{\overline{p}_f}{\overline{\rho}_f} F_f.
\end{equation*}

The resulting equation reads:

\begin{multline}\label{eq:enEq}
\sum_{f=1}^{Nf_i} \bigg[ \tilde{e}_f F_f - \frac{C_p}{C_v} \bigg( \frac{\mu}{Pr} + \frac{\mu_t}{Pr}_{t} \bigg)\Big \rvert_f *\\
\bigg( \lvert \bm{\mathcal{P}}_f \rvert \frac{\tilde{e}_i - \tilde{e}_j}{\lvert \bm{d} \rvert} + \bm{\mathcal{O}}_f \cdot \nabla \tilde{e}\rvert_f \bigg) \bigg] = \\
= - \sum_{f=1}^{Nf_i} \bigg(  \frac{\tilde{\bm{u}}_f \cdot \tilde{\bm{u}}_f}{2} + \frac{\overline{p}_f}{\overline{\rho}_f} \bigg) F_f.
\end{multline}

Also <ref> can be written into its matrix form as follow:

\begin{equation}\label{eq:mEnEq}
\bm{E}(\tilde{e}) = \bm{F}(\overline{p},\tilde{\bm{u}}).
\end{equation}

§.§ Pressure equation for compressible flows

By following what has been done in [42], let us localize <ref> at a generic $\Omega_i$ cell center, we get:

\begin{equation*}
\tilde{\bm{u}}_i = \frac{\bm{H}(\tilde{\bm{u}})}{a_i}  - \frac{\nabla \overline{p}_i}{a_i}.
\end{equation*}

Let us rename $\tilde{\bm{u}} = \tilde{\bm{u}}^*$ and $\overline{\rho} = \overline{\rho}^*$ both velocity and density we have at this point, after having solved the momentum equation, for a reason that will be clarified in a moment.

The mass flux, at the generic cell center, can be obtained as:

\begin{equation*}
\overline{\rho}_i^*\tilde{\bm{u}}_i^* = \overline{\rho}_i^* \frac{\bm{H}(\tilde{\bm{u}}^*)}{a_i}  - \overline{\rho}_i^* \frac{\nabla \overline{p}^{n-1}_i}{a_i}.
\end{equation*}

Since the pressure gradient has to be calculated explicitly, we indicate it as $\nabla \overline{p}^{n-1}$ meaning that the pressure field has to be previously calculated.

It is easy to realize that the set $\overline{\rho}^*, \tilde{\bm{u}}^*, \overline{p}^{n-1}$ will not satisfy the mass conservation constrain since velocity field has been evaluated by the use of the pressure gradient at step $n-1$. We can then imagine to introduce some corrections to all the terms so that $\overline{\rho} = \overline{\rho}^* + \overline{\rho}', \tilde{\bm{u}} = \tilde{\bm{u}}^* + \tilde{\bm{u}}', \overline{p} = \overline{p}^{n-1} + \overline{p}'$.

It is now possible to rewrite the mass flux as:

\begin{multline*}
\left(\overline{\rho}_i^* + \overline{\rho}_i'\right)\left(\tilde{\bm{u}}_i^* + \tilde{\bm{u}}_i'\right) = \left(\overline{\rho}_i^* + \overline{\rho}_i'\right) \left[\frac{\bm{H}(\tilde{\bm{u}}^*)}{a_i} + \frac{\bm{H}(\tilde{\bm{u}}')}{a_i} \right]\\
- \left(\overline{\rho}_i^* + \overline{\rho}_i'\right) \left[\frac{\nabla \overline{p}^{n-1}_i}{a_i} + \frac{\nabla \overline{p}_i'}{a_i} \right].
\end{multline*}

By the definition of compressibility $\Psi$, we can write $\overline{\rho} = \Psi \overline{p}$ and then $\rho = \overline{\rho}^* + \overline{\rho}' = \Psi \overline{p}^{n-1} + \Psi \overline{p}' \Rightarrow \overline{\rho}' = \Psi p - \Psi \overline{p}^{n-1} = \Psi \overline{p}' .$

We can then interpolate that expression to obtain the variables evaluations at the faces and finally sum over all the faces surrounding the cell $\Omega_i$ to get the mass conservation equation in its pressure correction shape:

\begin{multline*}
\sum_{f=1}^{Nf_i}\left(\overline{\rho}_i^* + \Psi \overline{p}_i' \right) \left[\frac{\bm{H}(\tilde{\bm{u}}^*)}{a_i} + \frac{\bm{H}(\tilde{\bm{u}}')}{a_i} \right] \Bigg\rvert_f  - \sum_{f=1}^{Nf_i}\left(\overline{\rho}_i^* + \underbrace{\overline{\rho}_i'}_{\star}\right) *\\
\left[\frac{\nabla \overline{p}^{n-1}_i}{a_i} + \frac{\nabla \overline{p}_i'}{a_i} \right]\Bigg\rvert_f = 0.
\end{multline*}

The $\star$ term can be neglected obtaining the correction equation for pressure. The only term that has to be modeled in some way is $\bm{H}(\tilde{\bm{u}}')$. In the SIMPLE-based algorithms the correction extra diagonal velocity term is neglected leading to the following final pressure correction equation:

\begin{multline}\label{eq:mPresEq}
\sum_{f=1}^{Nf_i}\left(\overline{\rho}_i^* + \Psi \overline{p}_i' \right) \left[\frac{\bm{H}(\tilde{\bm{u}}^*)}{a_i} \right] \Bigg\rvert_f  \\
- \sum_{f=1}^{Nf_i}\left(\overline{\rho}_i^* \right) \left[\frac{\nabla \overline{p}^{n-1}_i}{a_i} + \frac{\nabla \overline{p}_i'}{a_i} \right]\Bigg\rvert_f = 0.
\end{multline}

§ REDUCED ORDER MODELING ARCHITECTURE

§.§ Proper Orthogonal Decomposition procedure

The scope of this work is to find an efficient and reliable reduced order model to be able to solve <ref> for many different values of the parameter $\pi$ without solving the Finite Volume discretized equations every time from scratch. For this reason we decided to develop a new procedure based on a POD-Galerkin scheme.

The whole machinery is divided into two main steps: an offline phase which consists on the resolution of a certain number $N_{\pi}$ of full order solutions, trying to extract as much information as possible from this set, and an online phase consisting on the resolution of a dimensionally reduced problem for all the different needed parametric configurations. What is new in this method is to be capable of resulting as general as possible with respect to the selected full order turbulence model and, at the same time, as coherent as possible with respect to high fidelity solutions.

Let us introduce the training parameters set as $\mathbb{P} = \{ \pi_1, \ldots, \pi_{N_{\pi}} \}$. For every parameter $\pi_i \in \mathbb{P}$, the full order problem can be solved to obtain the corresponding solution $\bm{s}_i$. All these offline solutions are then stored into the snapshots matrix:

\begin{equation*}
\bm{S} = \begin{bmatrix}
s_{1_1} & s_{2_1} & \dotsc & s_{{N_{\pi}}_1}\\
\vdots & \vdots & \vdots & \vdots\\
s_{1_{N_h}} & s_{2_{N_h}} & \dotsc & s_{{N_{\pi}}_{N_h}}\\
\end{bmatrix}.
\end{equation*}

In our case we want to construct an online solver able to mimic the offline convergence dynamics. For this reason the use of a monolithic approach for the reduced problem is not a good choice: since the offline solutions are obtained relying on a segregated solver, also at the online level a segregated strategy has to be applied to obtain solutions as consistent as possible. For a discussion on a similar consistent approach in the context of explicit time integration schemes the reader is referred to [43].

To obtain an algorithm able to properly follow the behavior of the high fidelity algorithm, it is not possible to just rely on the converged solutions: starting from an attempted initial solution, there is a high possibility for the algorithm to get lost since it has not enough information on the path to follow to reach the minimum of the residual, i.e. the correct solution. In other words, the intermediate solutions provided by the SIMPLE algorithm cannot be accurately represented by the use of just the dynamics obtained by final converged iterations. Thus, by projecting intermediate equations over the subspaces spanned by those converged basis functions, a stable convergence of the algorithm may not be possible. To improve the convergence performances, we decided to enrich the set of snapshots by saving a certain amount of intermediate solutions $\bm{s}_i^j$ obtained during the offline iterations. The distance between exported intermediate solutions is set to $\Delta$ (see <ref>). By adding some non-physical solutions to the snapshots matrix, which is what is happening by inserting non-converged fields, we are somehow polluting the physical content but the convergence propertied of the algorithm are improved in any case. To reach a balance between convergence and reliability, $\Delta$ can be varied and the total amount $N_{int}$ of selected intermediate solutions can be modified. The new snapshots matrix then reads:

\begin{equation*}
\bm{S} = \left[ \bm{s}_1^1, \bm{s}_1^2, \ldots, \bm{s}_1^{N_{int}}, \bm{s}_1, \dots, \bm{s}_{N_{\pi}}^1, \bm{s}_{N_{\pi}}^2, \cdots, \bm{s}_{N_{\pi}}^{N_{int}}, \bm{s}_{N_{\pi}}\right],
\end{equation*}

where $\bm{s}_i^j$ is the solution obtained at the $(j \, \Delta)$-th iteration for the $i$-th offline parameter.

3[surf, domain=-1:1] (x^2) + (y^2) ;
(-0.8,1.8) – (-0.4,1.7) – (-1,1.5) – (-0.3,1.3) – (-0.7,1.1) – (-0.2,0.9) – (-0.4,0.7) – (-0.1,0.5) – (0,0.1);
[black!] (-0.8,1.8) circle (0.15cm);
at (-0.9,2) $\bm{s}_i^{init}$;
[black!] (-0.4,1.7) circle (0.1cm);
at (-0.25,1.7) $\bm{s}_i^1$;
[blue!] (-1,1.5) circle (0.1cm);
at (-1.15,1.5) $\bm{s}_i^2$;
[black!] (-0.3,1.3) circle (0.1cm);
at (-0.15,1.3) $\bm{s}_i^3$;
[blue!] (-0.7,1.1) circle (0.1cm);
at (-0.85,1.1) $\bm{s}_i^4$;
[black!] (-0.2,0.9) circle (0.1cm);
at (-0.05,0.9) $\bm{s}_i^5$;
[blue!] (-0.4,0.7) circle (0.1cm);
at (-0.55,0.7) $\bm{s}_i^6$;
[black!] (-0.1,0.5) circle (0.1cm);
at (0.05,0.5) $\bm{s}_i^7$;
[red!] (0,0.1) circle (0.2cm);
at (-0.2,0.1) $\bm{s}_i$;

Scheme of the snapshots selection for $\Delta = 2$: black dots are discarded intermediate solutions, blue dots are saved intermediate solutions while the red dot represents the final solution.

In a POD-Galerkin approach, the reduced order solution $\bm{s}^r$ is obtained as a linear combination of some precalculated basis functions $\bm{\xi}$:

\begin{equation*}
\bm{s}^r(\pi) = \sum_{i=1}^{N_r} \beta_i(\pi) \bm{\xi}_i(\bm{x}),
\end{equation*}

where $N_r < N_{\pi}$ is the number of basis functions to be used for the reconstruction and the $\beta_i$ are the coefficients depending only on the parameter representing the reduced solution.

Once provided a certain amount $N_t$ of high fidelity solutions, with $N_t > N_{\pi}$ because of the intermediate snapshots, the best reduced order model we can get is the one able to fully reproduce the training offline solutions with no error with respect to it. Of course this is not achievable but we would like the $L^2$ norm of the error $E_{ROM}$ between all the offline solutions and the respective online ones to be as low as possible:

\begin{equation*}
E_{ROM} = \sum_{i=1}^{N_t} \lvert \lvert \bm{s}^{ROM}_i - \bm{s}_i \rvert \rvert_{L^2} = \sum_{i=1}^{N_t} \Bigg|\Bigg| \sum_{j=1}^{N_r} \beta_j(\pi) \bm{\xi}_j(\bm{x}) - \bm{s}_i \Bigg| \Bigg|_{L^2}.
\end{equation*}

It is well known (see, e.g., [19]) that the basis functions best performing in this sense are the ones obtained through a Proper Orthogonal Decomposition (POD) applied to the snapshots matrix $\bm{S}$. The eigenproblem

\begin{equation*}
\bm{C} \bm{V} = \bm{V} \bm{\lambda},
\end{equation*}

has to be resolved, where $\bm{C} \in \mathbb{R}^{N_t \times N_t}$ is the correlation matrix containing all the $L^2$ inner products products of the snapshots between each others, $\bm{V} \in \mathbb{R}^{N_t \times N_t}$ is the matrix containing its eigevectors while $\bm{\lambda}$ is the diagonal matrix containing the eigenvalues.

The basis functions are then constructed as just a linear combination of the snapshots contained in $\bm{S}$:

\begin{equation*}
\bm{\xi}_i(\bm{x}) = \frac{1}{N_t} \sum_{j=1}^{N_t} \bm{V}_{ji} \bm{s}_j (\bm{x}).
\end{equation*}

The basis functions matrix is then defined as:

\begin{equation*}
\bm{\Xi} = \left[ \bm{\xi}_1, \dotsb, \bm{\xi}_{N_r} \right] \in \mathcal{R}^{N_h \times N_r} .
\end{equation*}

The interested reader may refer to [44, 45, 25] for a detailed explanation regarding POD approaches.

§.§ A segregated reduced approach for compressible flows

In <ref> we introduced the intention of this work: the construction of a new segregated approach for compressible flows. By the procedure explained in the previous section, it is possible to construct three different snapshots matrices for pressure, velocity and energy respectively: $\bm{S}_{p} \in \mathbb{R}^{N_h \times N_t}$, $\bm{S}_{u} \in \mathbb{R}^{d \, N_h \times N_t}$ and $\bm{S}_{e} \in \mathbb{R}^{N_h \times N_t}$ being $d$ the spacial dimension of the problem. They are employed to extract three basis functions sets so that all the different variables of <ref> can be written into their reduced form:

\begin{align*}
&\overline{p}^r = \sum_{i=1}^{N_p} a_i (\pi) \varphi_i (\bm{x}) = \bm{\Phi} \bm{a} \in \mathbb{Q}_r, \\
&\tilde{\bm{u}}^r = \sum_{i=1}^{N_u} b_i (\pi) \bm{\psi}_i (\bm{x}) = \bm{\Psi} \bm{b} \in \mathbb{V}_r, \\
&\tilde{e}^r = \sum_{i=1}^{N_e} c_i (\pi) \theta_i (\bm{x}) = \bm{\Theta} \bm{c} \in \mathbb{E}_r,
\end{align*}

where $\mathbb{Q}_r = \text{span} \{\varphi_i\}_{i=1}^{N_p} \subset \mathbb{Q}_h$, $\mathbb{V}_r = \text{span} \{\psi_i\}_{i=1}^{N_u} \subset \mathbb{V}_h$ and $\mathbb{E}_r = \text{span} \{\theta_i\}_{i=1}^{N_e}  \subset \mathbb{V}_h$ are the reduced spaces, $\bm{a} \in \mathbb{R}^{N_p}$, $\bm{b} \in \mathbb{R}^{N_u}$ and $\bm{c} \in \mathbb{R}^{N_e}$ are the vector containing the coefficients $a_i$, $b_i$ and $c_i$ depending only on the parameter value related to pressure, velocity and energy respectively, $\bm{\Phi} \in \mathbb{R}^{N_h \times N_p}$, $\bm{\Psi} \in \mathbb{R}^{d \, N_h \times N_u}$ and $\bm{\Theta} \in \mathbb{R}^{N_h \times N_e}$ are the matrices containing the modal basis functions $\varphi_i$, $\bm{\psi}_i$ and $\theta_i$ related to pressure, velocity and energy respectively while $N_p < N_t$, $N_u < N_t$ and $N_e < N_t$ are the numbers of modal basis functions selected for pressure, velocity and energy to reconstruct their reduced solutions. This means that the procedure explained in <ref> has to be applied three times to the three different solutions sets $\bm{S}_{p}$, $\bm{S}_{u}$ and $\bm{S}_{e}$.

The Reduced Order SIMPLE algorithm

Input: first attempt reduced pressure and velocity and energy coefficients $\bm{a}^\star$, $\bm{b}^\star$ and $\bm{c}^\star$; modal basis functions matrices for pressure, velocity and energy $\Phi$, $\Psi$ and $\Theta$

Output: reduced pressure, velocity and energy fields $\overline{p}_r$, $\tilde{\bm{u}}_r$ and $\tilde{e}_r$
[1]
From $\bm{a}^\star$, $\bm{b}^\star$ and $\bm{c}^\star$, reconstruct reduced fields $\overline{p}^\star$, $\tilde{\bm{u}}^\star$ and $\tilde{e}^\star$:
$$\overline{p}^\star = \Phi \bm{a}^\star {,} \hspace{1.5cm} \tilde{\bm{u}}^\star = \Psi \bm{b}^\star {,} \hspace{1.5cm} \tilde{e}^\star = \Theta \bm{c}^\star;$$

Evaluate the eddy viscosity field $\mu_t$;
Momentum predictor step : assemble <ref>, relax it employing prescribed under-relaxation factor $\alpha_u$, project it over the velocity basis functions $\psi_i$ and solve it to obtain new reduced velocity coefficients vector $\bm{b}^{\star\star}$;
Reconstruct the new reduced velocity $\tilde{\bm{u}}^{\star\star}$ and calculate the off-diagonal component $\bm{H}(\tilde{\bm{u}}^{\star\star})$;
Energy equation step : assemble <ref>, relax it employing prescribed under-relaxation factor $\alpha_e$, project it over the energy basis functions $\theta_i$ and solve it to obtain new reduced energy coefficients vector $\bm{c}^{\star\star}$;
Reconstruct the new reduced energy $\tilde{e}^{\star\star}$;
Calculate both density $\overline{\rho}^{\star\star}$ and temperature $\tilde{T}^{\star\star}$ fields starting from $\overline{p}^\star$, $\tilde{\bm{u}}^{\star\star}$ and $\tilde{e}^{\star\star}$ by the use of the state equation;
Pressure correction step: assemble <ref>, project it over the pressure basis functions $\varphi_i$ to get new reduced pressure coefficients $\bm{a}^{\star\star}$; then correct the velocity explicitly after having reconstructed the new pressure $\overline{p}^{\star\star}$;
Relax the pressure field with the prescribed under-relaxation factor $\alpha_p$. The under-relaxed field is called $\overline{p}^{ur}$;
$\overline{p}_r = \overline{p}^{ur}$, $\tilde{\bm{u}}_r = \tilde{\bm{u}}^{\star\star}$ and $\tilde{e}_r = \tilde{e}^{\star\star}$
Assemble the conservative face fluxes $F_{f}$:

$$F_{f} = \tilde{\bm{u}}_{f} \cdot \bm S_{f};$$

set $\overline{p}^\star = \overline{p}^{ur}$, $\tilde{\bm{u}}^{\star} = \tilde{\bm{u}}^{\star\star}$ and $\tilde{e}^{\star} = \tilde{e}^{\star\star}$;
iterate from step 1.

§.§ Turbulence treatment

every neuron input/.style=
minimum size=0.5cm,
every neuron input2/.style=
minimum size=0.5cm,
every neuron hidden/.style=
minimum size=0.5cm,
every neuron output/.style=
minimum size=0.5cm,
neuron missing/.style=
text height=0.333cm,
execute at begin node=$\vdots$
layer missing/.style=
text height=0.333cm,
execute at begin node=$\dots$

[x=1.5cm, y=1.5cm, >=stealth, font=, scale=0.6]

ı[count=] in 1, missing, 2, 3, missing, 4
=4 =6[every neuron input2/.try, neuron ı/.try] (input-ı) at (0,2.5-) [every neuron input/.try, neuron ı/.try] (input-ı) at (0,2.5-) ;

[count=] in 1, 2, missing,3
[every neuron hidden/.try, neuron /.try ] (hidden1-) at (2,2-) ;

[count=] in 1, missing, 2, 3
[every neuron hidden/.try, neuron /.try ] (hidden2-) at (4,2-) ;

[count=] in 1,missing,2
[every neuron output/.try, neuron /.try ] (output-) at (6,1.5-) ;

[<-] (input-3) – ++(-1,0)
node [above, midway] $b_1$;
[<-] (input-4) – ++(-1,0)
node [above, midway] $b_{N_{u}}$;
[<-] (input-1) – ++(-1,0)
node [above, midway] $\pi_1$;
[<-] (input-2) – ++(-1,0)
node [above, midway] $\pi_{N_{\pi}}$;

ł[count=ı] in 1,N_μ_t
[->] (output-ı) – ++(1,0)
node [above, midway] $m_{\l}$;

ıin 1,...,4
ȷin 1,...,3
[->] (input-ı) – (hidden1-ȷ);

ıin 1,...,3
ȷin 1,...,3
[->] (hidden1-ı) – (hidden2-ȷ);

ıin 1,...,3
ȷin 1,...,2
[->] (hidden2-ı) – (output-ȷ);

ł[count=from 0] in Input
layer, Hidden
layer  $f_{1}$
[align=center, above] at (*2,2) ł;
[align=center, above] at (4,2) Hidden
layer $f_2$;
[align=center, above] at (6,2) Output
layer $f_{out}$;

Schematic perspective of a fully connected neural network composed by an input layer, two hidden layers and an output layer, linking parameters $\pi_i$ and reduced velocity coefficients $b_i$ to reduced eddy viscosity coefficients $m_i$, being $N_{\pi}$ the number of parameters possibly existing in the problem.

In this work some assumptions were taken in <ref> and <ref> leading to a simplified FANS system, <ref>. Turbulence effects in <ref> are all due to the presence of the eddy viscosity field $\mu_t$. A technique has to be selected to model the eddy viscosity. Within this scope, many different approaches are possible (see [46, 47, 48, 49]).

To make our architecture as independent as possible on the turbulence model used during the offline phase to evaluate the $\mu_t$ field, we decided to combine a classical POD-Galerkin approach for what concerns the physical variables $\overline{p}, \tilde{\bm{u}}$ and $\tilde{e}$ together with a data driven scheme for what concerns the eddy viscosity evaluation in the Boussinesq hypothesis (see [50]).

Let us imagine to approximate the eddy viscosity field similarly to what has been done for all the other variables:

\begin{equation*}
\mu_{t_r} = \sum_{i=1}^{N_{\mu_t}} m_i (\bm{x}) \eta_i (\pi),
\end{equation*}

where $N_{\mu_t}$ is the number of basis functions selected to reconstruct the eddy viscosity field, $m_i$ are the coefficients depending only on the position $\bm{x}$ while $\eta_i$ are the $\mu_t$ basis functions depending only on the parameter. During the offline phase, together with all the other saved solutions, also the eddy viscosity fields are exported and stored. Those snapshots are then collected into the $\bm{S}_{\mu_t}$ matrix and used, as explained in <ref>, to obtain the requested basis functions.
For what concerns the spacial coefficients, they are evaluated through a Neural Network (NN) scheme linking the parameters of the problem $\pi_i$ and the reduced velocity coefficients $b_i$ to the $m_i$. In fact it is well known that, no matter what turbulence model is employed, the eddy viscosity $\mu_t$ depends on the velocity field but, especially for geometrically parametrized problem, it also depends on the parameter itself. The reduced problem is thus completely independent on the choice of the turbulence model and point $2$ into <ref> can be performed in an efficient way. This would not have been the case if turbulence equations were projected: in case there was the necessity of changing the adopted turbulence model, all the architecture had to be modified.

In this work we selected a fully connected Neural Network composed by an input layer, two hidden layers and an output layer. The input vector $\bm{z}$ and output vector $\bm{m}$ are defined as mentioned before:

\begin{equation*}
\bm{z}=
\begin{bmatrix}
\pi_1\\
\vdots\\
\pi_{N_{\pi}}\\
\vdots\\
\end{bmatrix}, \hspace{1.5cm}
\bm{m}=
\begin{bmatrix}
\vdots\\
\end{bmatrix}.
\end{equation*}

It is clear that the Neural Network has to be trained in some way. To this scope the snapshots contained into $\bm{S}_{\mu_t}$ are projected over their own basis functions $\eta_i$ to obtain the set of real coefficients $\{\bm{m}_i\}_{i=1}^{N_t}$. They can be compared with the NN estimated coefficients $\{\tilde{\bm{m}}_i\}_{i=1}^{N_t}$ into a loss function to target the training procedure. The loss function $\ell$ we adopted is a widely used quadratic one:

\begin{equation*}
\ell = ||\bm{m} - \tilde{\bm{m}}||_{L^2}.
\end{equation*}

The quantity $\mathcal{L}$ to be minimized during the training of the network is the sum of the loss function evaluated for all the different snapshots:

\begin{equation*}
\mathcal{L} = \sum_{i=1}^{N_t} ||\bm{m}_i - \tilde{\bm{m}}_i||_{L^2}.
\end{equation*}

The coefficients estimated by the network can be written as:

\begin{equation*}
\tilde{\bm{m}} = \bm{\textit{f}}_{out} \left( \bm{W}_{out} \, \bm{\textit{f}}_2\left( \bm{W_2} \bm{\textit{f}}_1 \left( \bm{W_1} \bm{x} + \bm{b}_1\right) + \bm{b}_2\right) + \bm{b}_{out} \right),
\end{equation*}

where $\bm{\textit{f}}_1$, $\bm{\textit{f}}_2$ and $\bm{\textit{f}}_{out}$ are the activation functions, $\bm{W_1}$, $\bm{W_2}$ and $\bm{W_{out}}$ are the weights while $\bm{b}_1$, $\bm{b}_2$ and $\bm{b}_{out}$ are the biases, related to the first and the second hidden layer and to the output layer respectively. For the hidden layers the best performing activation function appears to be the hyperbolic tangent while the output layer has been simply implemented as a linear combination of the received data. The previous formula can then be simplified as follows:

\begin{equation*}
\tilde{\bm{m}} = \bm{W}_{out} \, \tanh\left( \bm{W_2} \, \tanh \left( \bm{W_1} \bm{x} + \bm{b}_1\right) + \bm{b}_2\right) + \bm{b}_{out},
\end{equation*}

where $\tanh(\bm{y}) = \begin{bmatrix} \tanh(y_1)\\ \vdots\\ \tanh(y_{dim_y})\\ \end{bmatrix}$ being $\bm{y} = [y_1, \hdots, y_{dim_y}]$ is a generic vector quantity.

§ NUMERICAL RESULTS

§.§ Physical parametrization test case

The first test case we present in this work is a physically parametrized external flow: a NACA0012 airfoil is immersed into a fluid with variable viscosity $\mu$. The unperturbed velocity is fixed and is equal to $\tilde{\bm{u}}_{inlet}=\begin{bmatrix} 250, 0, 0 \end{bmatrix}^T$kg.m while the chord of the airfoil is equal to one. As already said, the viscosity can vary so that $\mu \in [10^{-5}, 10^{-2}]$. The speed of sound at the inlet can easily be evaluated by taking into consideration the thermophysical properties of the gas we are working with. We consider perfect gasses. Thus the specific heat transfer at constant pressure is sufficient to evaluate $\gamma = \frac{C_p}{C_v} = \frac{C_p}{C_p-R}$ where $C_p=1005\si{\joule\kilogram\per\kelvin}$ while $R=8,314 \si{\joule\mol\per\kelvin}$ is the constant for perfect gasses. We suppose our airfoil to move into air so that $M=28,9$  where $M$ stands for the molar weight. Temperature is fixed at $T=298\si{\K}$. Collecting all these data together, we end up with

$$C = \sqrt{\frac{\gamma R T}{M}} = 341.17 \si{\m\per\s}.$$

This means that at the inlet the mach number can be calculated as

$$\text{Mach}=\frac{\tilde{\bm{u}}_{inlet}}{C} \simeq 0.73.$$

For this test case, consequently, a compressible treatment for the flow is needed since we are approaching the transonic regime and compressible effects are pretty significant.
At the inlet, pressure is fixed to $10^5$ Pa. Then also the Reynolds number can be evaluated as

\text{Re}=\frac{\rho L \tilde{\bm{u}}_{inlet}}{\mu}=\frac{p L \tilde{\bm{u}}_{inlet} M}{\mu R T}.

The resulting Reynolds number is then $\text{Re}\in 2.92 \times \left[10^4, 10^7 \right]$, which clearly requires a treatment for turbulence since the system is operating in fully turbulent regime.

For the offline phase, $50$ random values have been selected: $\pi_i \in [10^{-5}, 10^{-2}] \ \text{for} \ i=1,\ldots,50$ where $[\pi_1, \ldots, \pi_{50}]=\mathbb{P}$. Full order eddy viscosity is calculated by the resolution of a $k-\omega$ turbulence model.

scaled y ticks=false

xlabel=POD mode no.,
ylabel=Cumulated eigenvalues,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,
y tick label style=
/pgf/number format/.cd,
legend style=at=(0.98,0.4),anchor=north east,
+[no marks]
table[x=xAx, y=eigsU] physicalData/eigs.dat;

+[no marks]
table[x=xAx, y=eigsP] physicalData/eigs.dat;

+[no marks]
table[x=xAx, y=eigsE] physicalData/eigs.dat;

Velocity, Pressure, Energy;

at=(plot1.right of south east),
anchor=left of south west,
xlabel=POD mode no.,
ylabel=Cumulated eigenvalues,
legend style=at=(0.98,0.3),anchor=north east,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,
y tick label style=
/pgf/number format/.cd,

+[no marks]
table[x=xAx, y=eigsNut] physicalData/eigs.dat;

Eddy viscosity;

Cumulated eigenvalues trends.

<ref> shows the trends of the cumulated eigenvalues for velocity, pressure, energy and eddy viscosity. As we may notice, by just considering a few modes for every variable, the amount of discarded information is pretty low. For this reason, just the first $20$ modal basis functions have been selected for velocity, pressure and energy while $30$ modal basis functions are used to reconstruct the eddy viscosity field. This is due to the fact that, analyzing <ref>, it is clear that a higher number of basis functions are needed in order to approach the unity in the cumulated eigenvalues plot.

For what concerns the neural network for the eddy viscosity coefficients, as explained in <ref>, two hidden layer are present, the first one composed by $256$ neurons and the second one composed by $64$ neurons, resulting in a fully connected network where only $\tanh$ activation functions are used. Offline solutions, including the intermediate steps, are retained to train the network.

legend style=at=(0.98,0.98),anchor=north east,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,
x tick label style=
/pgf/number format/.cd,
scaled x ticks=false
+[no marks]
table[x=xAxL, y=trainLoss] physicalData/loss2030.dat;

+[no marks]
table[x=xAxL, y=testLoss] physicalData/loss2030.dat;

Train, Test;

Loss function decay for both train and test sets.

The training procedure is carried out in $2 \times 10^3$ epochs. $20$ new random offline solutions have been performed to obtain a testing set which was not correlated to the solutions used for the training stage. A mean squared error loss function is used to evaluate the reconstruction capability of the network for both training and testing sets. The decay behaviour of both losses are depicted in <ref>. The training stage has been stopped afted $2 \times 10^3$ epochs since test loss was no more decreasing and the distance between test and train losses was starting to increase significantly.

xlabel=Parameter no.,
ylabel=$L^2$ relative error,
legend style=at=(0.98,0.7),anchor=north east,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,

+[only marks]
table[x=parameter, y=errorU] physicalData/errors.dat;
+[only marks]
table[x=parameter, y=errorP] physicalData/errors.dat;
+[only marks]
table[x=parameter, y=errorE] physicalData/errors.dat;

Velocity, Pressure, Energy;

at=(plot1.right of south east),
anchor=left of south west,
xlabel=Parameter no.,
ylabel=$L^2$ relative error,
legend style=at=(0.98,0.3),anchor=north east,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,

+[only marks]
table[x=parameter, y=errorNut] physicalData/errors.dat;

Eddy viscosity;

$L^2$ norm relative errors.

<ref>, left, shows the $L^2$ norm relative errors for all the different parameters in the online set concerning velocity, pressure and internal energy. <ref>, right, shows the $L^2$ norm relative error for the eddy viscosity between full order and reduced order for the whole online parameter set. As we may notice, even if the order of magnitude of the $\nu_t$ error is equal to $10^{-2}$, it is sufficient to ensure a lower error for the quantities of interest, i.e. velocity, pressure and energy. By this observation we are allowed to employ such a small neural network which is not compromising the computational cost, still ensuring good performances.

Comparison between full order (top picture) and reduced order (bottom picture) solutions: velocity magnitude (1st, 2nd pictures), pressure (3rd, 4th pictures) and energy (5th, 6th) on the right. These fields refer to the resolution of the problem for $\pi = \mu = 0.21 \times 10^{-3}$ which has been selected as a random value in the online parameter set.

Comparison between full order (top picture) and reduced order (bottom picture) eddy viscosity solutions. These fields refer to the resolution of the problem for $\pi = \mu = 0.21 \times 10^{-3}$ which has been selected as a random value in the online parameter set.

In <ref> and <ref> is depicted a comparison between full order and reduced order solutions for a random value of the parameter, included in the online set. By analysing the depicted fields, full order and reduced order solutions appear to be very similar and most important areas into the domain, i.e. the zone surrounding the airfoil together with the wake created by the body, are well reconstructed.

§.§ Geometrical parametrization test case

In this section we present a second test case, focused on a geometrically parametrized problem: the shape of the airfoil used into <ref> is modified by the use of a bump function. In particular the foil is divided into a top and a bottom part by the chord. The bump function depicted in <ref> is added to the top and subtracted to the bottom surface, premultiplied by two different amplitude scalar factors: every solution is parametrized uniquely by two different scalar values. We use the same thermophysical properties used for <ref> but the dynamic viscosity is fixed and equal to $1.74 \times 10^{-5} Pa \, s$. Moreover the inlet velocity has been slightly decreased since the random modification of the geometry may lead to high curvature areas where the flow could eventually become supersonic: $\tilde{\bm{u}}_{inlet}=\begin{bmatrix} 170, 0, 0 \end{bmatrix}^T \si{\m\per\s}$. This means that the Mach number at the inlet is now around $0.5$.

For the offline phase, $50$ random values have been selected:$\pi_{top_i}, \pi_{bottom_i} \in [0,0.1]$ for $i= 1,...,50$ where

$$\begin{bmatrix} &\pi_{top_1}, &\pi_{bottom_1} \\ &\vdots &\vdots \\ &\pi_{top_{50}}, &\pi_{bottom_{50}} \end{bmatrix} = \mathbb{P}.$$

Full  order  eddy  viscosity  is  calculated  by  the  resolution  of  a $k-\omega$ turbulence model.

The general POD approach described in <ref> is not directly applicable to a geometrical parametrization problem since the $L^2$ norm used for the inner products is not well defined in case of multiple different domains. The mesh in our case is moved thanks to a Radial Basis Functions (RBF) algorithm where the points on the moving boundaries are displaced by the application of the desired law and their displacements are used as boundary conditions for an interpolation procedure, performed in order to move all the remaining points of the grid. The interested reader may find a deeper explanation of this technique in [51] or some applications in [52] and [53]. By exploiting the aforementioned method, the mesh is modified for each offline solution. To take into account the fact that all the snapshots are defined over a different mesh, the grid is taken back to its undeformed state before starting the POD procedure: the mass matrix we consider to evaluate the norms is then the reference unperturbed one.

legend style=at=(0.98,0.4),anchor=north east,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,
xmin=-0.1, xmax=1.1,
ymin=-0.01, ymax=0.15]

Shape of the employed bump function.

scaled y ticks=false

xlabel=POD mode no.,
ylabel=Cumulated eigenvalues,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,
y tick label style=
/pgf/number format/.cd,
legend style=at=(0.98,0.4),anchor=north east,
+[no marks]
table[x=xAx, y=eigsU] geometricalData/eigs.dat;

+[no marks]
table[x=xAx, y=eigsP] geometricalData/eigs.dat;

+[no marks]
table[x=xAx, y=eigsE] geometricalData/eigs.dat;

Velocity, Pressure, Energy;

at=(plot1.right of south east),
anchor=left of south west,
xlabel=POD mode no.,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,
ylabel=Cumulated eigenvalues,
legend style=at=(0.98,0.3),anchor=north east,

+[no marks]
table[x=xAx, y=eigsNut] geometricalData/eigs.dat;

Eddy viscosity;

Cumulated eigenvalues trends.

To test the online performances, $20$ new scalar amplitude couples have been randomly selected. $30$ modal basis functions have been picked  for the reconstruction of velocity, pressure and internal energy fields while $15$ modal basis functions have been emploied for $\nu_t$. This choice is supported by what is shown in <ref>: the increasing trend of the cumulated eigenvalues is pretty fast and this fact allows the discarding of the modes higher than the fixed quantity. For every new parameter couple the mesh motion has to be performed but the procedure is very efficient since the coefficients for the RBF have to be evaluated and stored just once (see [53]).

legend style=at=(0.98,0.98),anchor=north east,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,
x tick label style=
/pgf/number format/.cd,
scaled x ticks=false
+[no marks]
table[x=xAxL, y=trainLoss] geometricalData/loss3015.dat;

+[no marks]
table[x=xAxL, y=testLoss] geometricalData/loss3015.dat;

Train, Test;

Loss function decay for both train and test sets.

The same neural network used for <ref> is employed here for what concerns the eddy viscosity. Again, looking at <ref>, it can be noticed that the learning of the net seems to stabilize after $2\times10^4$ epochs which is the threshold we fixed for the training procedure.

xlabel=Parameter no.,
ylabel=$L^2$ relative error,
legend style=at=(0.98,0.5),anchor=north east,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,

+[only marks]
table[x=parameter, y=errorU] geometricalData/errors.dat;
+[only marks]
table[x=parameter, y=errorP] geometricalData/errors.dat;
+[only marks]
table[x=parameter, y=errorE] geometricalData/errors.dat;

Velocity, Pressure, Energy;

at=(plot1.right of south east),
anchor=left of south west,
xlabel=Parameter no.,
ylabel=$L^2$ relative error,
legend style=at=(0.98,0.8),anchor=north east,
grid style=line width=.1pt, draw=gray!10,
minor tick num=1,

+[only marks]
table[x=parameter, y=errorNut] geometricalData/errors.dat;

Eddy viscosity;

$L^2$ norm relative errors.

The resulting $L^2$ norm errors for all the parameter couples in the online set are shown in <ref>. Once again a discrepancy of about one order of magnitude can be noticed between the relative errors for the quantities of interest and the one calculated for the eddy viscosity. This is due to the fact that we are using a very simple and small network but it reveals to be reliable enough to make the online algorithm work fine.

Comparison between full order (top picture) and reduced order (bottom picture) solutions: velocity magnitude (1st, 2nd pictures), pressure (3rd, 4th pictures) and energy (5th, 6th pictures). These fields refer to the resolution of the problem for $\pi_{top}\simeq0.004$ and $\pi_{bottom}\simeq0.086$ which has been selected as a random value in the online parameter set.

Comparison between full order (top picture) and reduced order (bottom picture) eddy viscosity solutions. These fields refer to the resolution of the problem for $\pi_{top}\simeq0.004$ and $\pi_{bottom}\simeq0.086$ which has been selected as a random value in the online parameter set.

Comparison between full order (top picture) and reduced order (bottom picture) solutions: velocity magnitude (1st, 2nd pictures), pressure (3rd, 4th pictures) and energy (5th, 6th pictures). These fields refer to the resolution of the problem for $\pi_{top}\simeq0.095$ and $\pi_{bottom}\simeq0.003$ which has been selected as a random value in the online parameter set.

Comparison between full order (left) and reduced order (right) eddy viscosity solutions. These fields refer to the resolution of the problem for $\pi_{top}\simeq0.095$ and $\pi_{bottom}\simeq0.003$ which has been selected as a random value in the online parameter set.

In <ref>, <ref>, <ref> and <ref> a comparison between offline and online solutions is depicted for two different parameter couples selected from the online set. Even if the two solutions are obtained for airfoil geometries that are perturbed in opposite directions, in both cases the method exhibits good reliability properties.

§ CONCLUSIONS AND FUTURE PERSPECTIVES

In this work we focused on compressible flows proposing a new mixed technique, capable of merging the reliability of Galerkin-projection methods together with the versatility of data-driven strategies. The good results obtained for both a physical parametrization test case and a geometrically parametrized problem make this approach very promising. From one side the possibility to freely select the turbulence model avoiding the necessity of changing the whole architecture is attractive, from another point of view the warranty of a strong connection with physical aspects given by the projection of the conservation laws is reassuring.

The segregated compressible algorithm proposed in <ref> also introduces a way to provide accurate reduced solutions without any kind of stabilization: the employment of a decoupled approach for the compressible turbulent Navier Stokes equations relies on the chipping of the saddle point formulation. For this reason no stabilization for pressure is required: as shown in both <ref> and <ref>, pressure field solutions do not exhibit significant instability or inaccuracy issues. This aspect helps the procedure on being more consistent without pollution of the resulting solution due to stabilization.

A natural extension of this work is the application of neural networks to also approximate the functional evaluations required by the online phase to overtake the necessity of reconstruct the full fields at each iteration. This aspect would increase the performances but it has to be carefully calibrated to avoid possible drifting of the algorithm resulting on the loss of the convergence.

A final aspect that can be improved is the neural network itself: a weighted strategy where eigenvalues play a relevant role in the loss function would, in principle, enhance the training stage since the first modal basis functions, represented by the highest eigenvalues, are the  most significant ones on the reconstruction procedure.

§ ACKNOWLEDGMENTS

This research has been supported by the European Union Funding for Research and Innovation -Horizon 2020 Program- in the framework of European Research Council Executive Agency: Consolidator Grant H2020 ERC CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics" (PI Prof. Gianluigi Rozza) and by the H2020 MSCA RISE ARIA (grant 872442) project. We also acknowledge the INDAM-GNCS project “Tecniche Numeriche Avanzate per Applicazioni Industriali", the support by MIUR (Italian Ministry for University and Research) FARE-X-AROMA-CFD project and PRIN “Numerical Analysis for Full and Reduced Order Methods for Partial Differential Equations" (NA-FROM-PDEs). The main computations in this work were carried out by the usage of ITHACA-FV [54], a library maintained at SISSA mathLab, an implementation in OpenFOAM [55] for reduced order modelling techniques. We acknowledge developers and contributors of each of the aforementioned libraries.

[1]
J. D. Anderson and J. F. Wendt, Computational fluid dynamics, vol. 206.
Springer, 1995.

[2]
J. F. Wendt, Computational Fluid Dynamics: An Introduction.
Springer Science & Business Media, 2008.

[3]
J. Blazek, Computational Fluid Dynamics: Principles and Applications.
Butterworth-Heinemann, 2015.

[4]
A. Jameson, “Aerodynamic Shape Optimization Using the Adjoint Method,” Lectures at the Von Karman Institute, Brussels, 2003.

[5]
J. Samareh, “Aerodynamic Shape Optimization Based on Free-Form
Deformation,” in 10th AIAA/ISSMO multidisciplinary analysis and
optimization conference, p. 4630, 2004.

[6]
L. Formaggia, A. Quarteroni, and A. Veneziani, Cardiovascular
Mathematics: Modeling and simulation of the circulatory system, vol. 1.
Springer Science & Business Media, 2010.

[7]
B. Murthy, R. Ghadge, and J. Joshi, “CFD simulations of gas–liquid–solid
stirred reactor: Prediction of critical impeller speed for solid
suspension,” Chemical Engineering Science, vol. 62, no. 24,
pp. 7184–7195, 2007.

[8]
J. Van Baten and R. Krishna, “CFD simulations of mass transfer from Taylor
bubbles rising in circular capillaries,” Chemical Engineering
Science, vol. 59, no. 12, pp. 2535–2545, 2004.

[9]
F. Castellani, M. Burlando, S. Taghizadeh, D. Astolfi, and E. Piccioni, “Wind
energy forecast in complex sites with a hybrid neural network and CFD based
method,” Energy Procedia, vol. 45, pp. 188–197, 2014.

[10]
N. Morozova, F. Trias, R. Capdevila, C. D. Pérez-Segarra, and A. Oliva,
“On the feasibility of affordable high-fidelity CFD simulations for indoor
environment design and control,” Building and Environment, vol. 184,
p. 107144, 2020.

[11]
D. Amsallem and C. Farhat, “Interpolation Method for Adapting Reduced-Order
Models and Application to Aeroelasticity,” AIAA journal, vol. 46,
no. 7, pp. 1803–1813, 2008.

[12]
I. Tezaur, J. Fike, K. Carlberg, M. Barone, D. Maddix, E. Mussoni, and
M. Balajewicz, “Advanced Fluid Reduced Order Models for Compressible
Flow,” in Tech. Rep., Sandia National Lab.(SNL-NM), Albuquerque, NM
(United States), 2017.

[13]
“6 model reduction in computational aerodynamics,” in Applications,
pp. 201–236, De Gruyter, Dec. 2020.

[14]
J. Weller, E. Lombardi, M. Bergmann, and A. Iollo, “Numerical methods for
low-order modeling of fluid flows based on POD,” International Journal
for Numerical Methods in Fluids, vol. 63, no. 2, pp. 249–268, 2010.

[15]
G. Stabile and G. Rozza, “Finite volume POD-Galerkin stabilised reduced order
methods for the parametrised incompressible Navier–Stokes equations,” Computers & Fluids, vol. 173, pp. 273–284, 2018.

[16]
I. Akhtar, A. H. Nayfeh, and C. J. Ribbens, “On the stability and extension
of reduced-order Galerkin models in incompressible flows,” Theoretical
and Computational Fluid Dynamics, vol. 23, no. 3, pp. 213–237, 2009.

[17]
J. Baiges, R. Codina, and S. R. Idelsohn, “Reduced-Order Modelling Strategies
for the Finite Element Approximation of the Incompressible Navier-Stokes
Equations,” in Numerical Simulations of Coupled Problems in
Engineering, pp. 189–216, Springer, 2014.

[18]
J. Burkardt, M. Gunzburger, and H.-C. Lee, “POD and CVT-based reduced-order
modeling of Navier–Stokes flows,” Computer methods in applied
mechanics and engineering, vol. 196, no. 1-3, pp. 337–355, 2006.

[19]
K. Kunisch and S. Volkwein, “Galerkin Proper Orthogonal Decomposition Methods
for a General Equation in Fluid Dynamics,” SIAM Journal on Numerical
analysis, vol. 40, no. 2, pp. 492–515, 2002.

[20]
O. San and T. Iliescu, “Proper orthogonal decomposition closure models for
fluid flows: Burgers equation,” International Journal of Numerical
Analysis & Modeling, 2014.

[21]
A. Dumon, C. Allery, and A. Ammar, “Proper Generalized Decomposition method
for incompressible Navier–Stokes equations with a spectral
discretization,” Applied Mathematics and Computation, vol. 219,
no. 15, pp. 8145–8162, 2013.

[22]
F. Chinesta, P. Ladeveze, and E. Cueto, “A Short Review on Model Order
Reduction Based on Proper Generalized Decomposition,” Archives of
Computational Methods in Engineering, vol. 18, no. 4, pp. 395–404, 2011.

[23]
P. J. Schmid, “Dynamic Mode Decomposition of numerical and experimental
data,” Journal of Fluid Mechanics, vol. 656, pp. 5–28, July 2010.

[24]
N. J. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic Mode
Decomposition: Data-Driven Modeling of Complex Systems.
SIAM, 2016.

[25]
J. S. Hesthaven, G. Rozza, and B. Stamm, Certified Reduced Basis Methods
for Parametrized Partial Differential Equations, vol. 590.
Springer, 2016.

[26]
K. Urban, S. Volkwein, and O. Zeeb, “Greedy Sampling Using Nonlinear
Optimization,” in Reduced Order Methods for modeling and computational
reduction, pp. 137–157, Springer, 2014.

[27]
J. S. Hesthaven and S. Ubbiali, “Non-intrusive reduced order modeling of
nonlinear problems using neural networks,” Journal of Computational
Physics, vol. 363, pp. 55–78, 2018.

[28]
Q. Wang, J. S. Hesthaven, and D. Ray, “Non-intrusive reduced order modeling of
unsteady flows using artificial neural networks with application to a
combustion problem,” Journal of computational physics, vol. 384,
pp. 289–307, 2019.

[29]
A. T. Mohan and D. V. Gaitonde, “A Deep Learning based Approach to Reduced
Order Modeling for Turbulent Flow Control using LSTM Neural Networks,” arXiv preprint arXiv:1804.09269, 2018.

[30]
A. Mannarino and P. Mantegazza, “Nonlinear aeroelastic reduced order modeling
by recurrent neural networks,” Journal of Fluids and Structures,
vol. 48, pp. 103–121, 2014.

[31]
D. Papapicco, N. Demo, M. Girfoglio, G. Stabile, and G. Rozza, “The Neural
Network shifted-Proper Orthogonal Decomposition: a Machine Learning Approach
for Non-linear Reduction of Hyperbolic Equations,” Computer Methods in
Applied Mechanics and Engineering, vol. 392, p. 114687, 2022.

[32]
F. Romor, G. Stabile, and G. Rozza, “Non-linear manifold ROM with
Convolutional Autoencoders and Reduced Over-Collocation method,” Submitted, arXiv:2203.00360, 2022.

[33]
B. Haasdonk, “Convergence Rates of the POD–Greedy Method,” ESAIM:
Mathematical Modelling and Numerical Analysis, vol. 47, pp. 859–873, Apr.

[34]
Y. Wang, I. M. Navon, X. Wang, and Y. Cheng, “2D Burgers equation with large
Reynolds number using POD/DEIM and calibration,” International Journal
for Numerical Methods in Fluids, vol. 82, no. 12, pp. 909–931, 2016.

[35]
M. Sen, K. Bhaganagar, and V. Juttijudata, “Application of proper orthogonal
decomposition (POD) to investigate a turbulent boundary layer in a channel
with rough walls,” Journal of Turbulence, p. N41, 2007.

[36]
S. Hijazi, G. Stabile, A. Mola, and G. Rozza, “Data–Driven POD–Galerkin
reduced order model for turbulent flows,” Journal of Computational
Physics, vol. 416, p. 109513, 2020.

[37]
S. Georgaka, G. Stabile, K. Star, G. Rozza, and M. J. Bluck, “A hybrid
reduced order method for modelling turbulent heat transfer problems,” Computers & Fluids, vol. 208, p. 104615, 2020.

[38]
G. Stabile, F. Ballarin, G. Zuccarino, and G. Rozza, “A reduced order
variational multiscale approach for turbulent flows,” Advances in
Computational Mathematics, vol. 45, no. 5-6, pp. 2349–2368, 2019.

[39]
D. C. Wilcox et al., Turbulence modeling for CFD, vol. 2.
DCW industries La Canada, CA, 1998.

[40]
R. Eymard, T. Gallouët, and R. Herbin, “Finite volume methods,” Handbook of numerical analysis, vol. 7, pp. 713–1018, 2000.

[41]
F. Moukalled, L. Mangani, and M. Darwish, The finite volume method in
computational fluid dynamics, vol. 113.
Springer, 2016.

[42]
H. Jasak, “Error analysis and estimation for the finite volume method with
applications to fluid flows.,” PhD thesis, 1996.

[43]
S. K. Star, B. Sanderse, G. Stabile, G. Rozza, and J. Degroote, “Reduced order
models for the incompressible navier-stokes equations on collocated grids
using a `discretize-then-project' approach,” International Journal for
Numerical Methods in Fluids, vol. 93, no. 8, pp. 2694–2722, 2021.

[44]
A. Quarteroni, A. Manzoni, and F. Negri, Reduced Basis Methods for
Partial Differential Equations: an Introduction, vol. 92.
Springer, 2015.

[45]
P. Benner, M. Ohlberger, A. Patera, G. Rozza, and K. Urban, Model
reduction of parametrized systems.
Springer, 2017.

[46]
M. S. Siddiqui, A. Rasheed, and T. Kvamsdal, “Numerical assessment of RANS
turbulence models for the development of data driven reduced order models,”
Ocean Engineering, vol. 196, p. 106799, Jan. 2020.

[47]
D. C. Wilcox and R. Traci, “A complete model of turbulence,” in 9th
Fluid and PlasmaDynamics Conference, p. 351, 1976.

[48]
G. Alfonsi, “Reynolds-Averaged Navier–Stokes Equations for Turbulence
Modeling,” Applied Mechanics Reviews, vol. 62, no. 4, 2009.

[49]
L. Sirovich, “Turbulence and the dynamics of coherent structures, Parts
I-III,” Quarterly of applied mathematics, pp. 561–590, 1987.

[50]
F. G. Schmitt, “About Boussinesq's turbulent viscosity hypothesis: historical
remarks and a direct evaluation of its validity,” Comptes Rendus
Mécanique, vol. 335, no. 9-10, pp. 617–627, 2007.

[51]
A. De Boer, M. Van der Schoot, and H. Bijl, “Mesh deformation based on radial
basis function interpolation,” Computers & structures, vol. 85,
no. 11-14, pp. 784–795, 2007.

[52]
G. Stabile, M. Zancanaro, and G. Rozza, “Efficient geometrical parametrization
for finite-volume-based reduced order methods,” International Journal
for Numerical Methods in Engineering, vol. 121, no. 12, pp. 2655–2682,

[53]
M. Zancanaro, M. Mrosek, G. Stabile, C. Othmer, and G. Rozza, “Hybrid Neural
Network Reduced Order Modelling for Turbulent Flows with Geometric
Parameters,” Fluids, vol. 6, 2021.

[54]
“ITHACA-FV website.” <https://github.com/mathLab/ITHACA-FV>.

[55]
“OpenFOAM website.” <https://openfoam.com/>.