text:
\begin{document}

\title{Parametric solutions of turbulent incompressible flows in OpenFOAM via the proper generalised decomposition}

\author{
\renewcommand{\thefootnote}{\arabic{footnote}}
			  Vasileios Tsiolakis\footnotemark[1]\textsuperscript{ \ ,}\footnotemark[2]\textsuperscript{ \ ,}\footnotemark[3] \ ,
			  Matteo Giacomini\footnotemark[2]\textsuperscript{ \ ,}\footnotemark[4]\textsuperscript{ \ ,}*  ,
			  Ruben Sevilla\footnotemark[3], \\
\renewcommand{\thefootnote}{\arabic{footnote}}
             Carsten Othmer\footnotemark[1] \ and
             Antonio Huerta\footnotemark[2]\textsuperscript{ \ ,}\footnotemark[4]
}

\date{\today}

\maketitle

\renewcommand{\thefootnote}{\arabic{footnote}}

\footnotetext[1]{Volkswagen AG, Brieffach 011/1777, D-38436, Wolfsburg, Germany}
\footnotetext[2]{Laboratori de C\`alcul Num\`eric (LaC\`aN), ETS de Ingenieros de Caminos, Canales y Puertos, Universitat Polit\`ecnica de Catalunya, Barcelona, Spain}
\footnotetext[3]{Zienkiewicz Centre for Computational Engineering, College of Engineering, Swansea University, Wales, UK}
\footnotetext[4]{Centre Internacional de M\`etodes Num\`erics en Enginyeria (CIMNE), Barcelona, Spain.
\vspace{5pt}\\
* Corresponding author: Matteo Giacomini. \textit{E-mail:} \texttt{matteo.giacomini@upc.edu}
}

\begin{abstract}
An \emph{a priori} reduced order method based on the proper generalised decomposition (PGD) is proposed to compute parametric solutions involving turbulent incompressible flows of interest in an industrial context, using OpenFOAM.

The PGD framework is applied for the first time to the incompressible Navier-Stokes equations in the turbulent regime, to compute a generalised solution for velocity, pressure and turbulent viscosity, explicitly depending on the design parameters of the problem.

In order to simulate flows of industrial interest, a minimally intrusive implementation based on OpenFOAM SIMPLE algorithm applied to the Reynolds-averaged Navier-Stokes equations with the Spalart-Allmaras turbulence model is devised.

The resulting PGD strategy is applied to parametric flow control problems and achieves both qualitative and quantitative agreement with the full order OpenFOAM solution for convection-dominated fully-developed turbulent incompressible flows, with Reynolds number up to one million.
\end{abstract}

\begin{keywords}
Reduced order models, proper generalised decomposition, turbulent incompressible flows, parametrised flows, OpenFOAM
\end{keywords}

\section{Introduction}
\label{sc:Intro}

Parametric studies involving flows of industrial interest require robust computational fluid dynamics (CFD) solvers and efficient strategies to simulate multiple queries of the same problem.

Finite volume (FV) methods represent the most common approach in industry to perform flow simulations~\cite{LeVeque:02,Toro:09,Sonar-MS:07,Ohlberger-BHO:17,Eymard-EGH:00,RS-SGH:18,RS-VGSH:20,MG-RS-20} and different strategies have been proposed to simulate flows in the turbulent regime~\cite{Deardorff-70,Hughes-HMJ-00,Spalart-TSSS-00}.

A widespread approach is represented by the Reynolds-averaged Navier-Stokes (RANS) equations~\cite{Reynolds-1895} coupled with the one-equation Spalart-Allmaras (SA) turbulence model~\cite{SpalartAllmaras-92}.

This work focuses on such strategy and relies on its cell-centred FV implementation available in OpenFOAM~\cite{OpenFOAM} and validated by the industry.

When the simulation requires testing a large number of different configurations - e.g. for shape optimisation, uncertainty quantification, inverse and control problems - numerical strategies to reduce the cost of the overall computation are critical. Reduced order models (ROM)~\cite{AH-CHRW:17,Gunzburger-PWG-18} construct an approximation of the solution in a lower dimensional space, for which an appropriate basis needs to be devised.

It is known that numerical difficulties arise when reduced basis (RB) and proper orthogonal decomposition (POD) techniques are applied to convection-dominated problems~\cite{Iliescu-IW-13,Pacciarini-PR-14,Iliescu-GIJW-15}.

This is especially critical in the context of flow simulations when the Reynolds number is increased and turbulent phenomena need to be accounted for.

More precisely, the most relevant POD modes are associated with the highest energy scales of the problem under analysis, whereas small scales, which play a critical role in the dissipation of turbulent kinetic energy, are poorly represented by POD-ROM~\cite{Aubry-AHLS-88}.

To remedy this issue, closure models stemming from traditional description of turbulence have been extended to ROMs, leading to Galerkin projection-based POD with dynamic subgrid-scale~\cite{Iliescu-WABI-12}, variational multiscale~\cite{Iliescu-IW-14,Rozza-SBZR-19}, $k-\omega$ SST~\cite{Rozza-LCLR-16} models and to a certified Smagorinsky RB strategy~\cite{Ballarin-RAMBR-17}.

Moreover, strategies to improve efficiency and accuracy of POD-ROM in the context of realistic and turbulent flows have been proposed by coupling the projection-based framework with residual minimisation~\cite{Farhat-CFCA-13}, nonlinear least-squares optimisation~\cite{Zimmermann-ZG-12}, interpolation based on radial basis functions~\cite{Rozza-GSSRB} and a constrained greedy approach~\cite{Taddei-FMPT-18}.

In the context of machine learning-based reduced order models~\cite{Iliescu-XMRI-18,Hesthaven-GH-19,Hesthaven-WHR-19}, a strategy coupling a traditional projection-based POD for velocity and pressure with a data-driven technique for the eddy viscosity has been recently proposed in~\cite{Rozza-HSMR-20}.

All above contributions involve the development of \emph{a posteriori} ROMs, namely RB and POD, in which the basis of the low-dimensional approximation space is computed starting from a set of snapshots. On the contrary, PGD~\cite{Chinesta-Keunings-Leygue,Chinesta-CLBACGAAH:13} constructs a reduced basis of separable functions explicitly depending on space and on user-defined parameters, with no \emph{a priori} knowledge of the solution of the problem.

The resulting PGD \emph{computational vademecum} provides a generalised solution which can thus be efficiently evaluated in the online phase via interpolation in the parametric space, that is, no extra problem needs to be solved in the low-dimensional reduced space as in POD.

In the context of flow problems, PGD was originally utilised to develop efficient solvers for the incompressible Navier-Stokes equations by separating spatial directions~\cite{Allery-DAA-10} and space and time~\cite{Allery-DAA-11,Allery-LA-14}.

In addition, problems involving parametrised geometries have been solved using PGD~\cite{AH-AHCCL:14,SZ-ZDMH:15}, with special emphasis on incompressible flows in geometrically parametrised domains~\cite{PD-DZH:17,RS-SZH:20}.

To foster the application of \emph{a priori} model order reduction techniques to problems of industrial interest, a non-intrusive PGD implementation in the CFD software OpenFOAM has been recently proposed in~\cite{Tsiolakis-TGSOH-20} to solve parametrised incompressible Navier-Stokes equations in the laminar regime.

Following the work on PGD for convection phenomena~\cite{2010-IJMF-GDCCDH,Gonzalez-GCCDH-13} and for viscous incompressible Navier-Stokes flows~\cite{Tsiolakis-TGSOH-20}, the present contribution proposes the first \emph{a priori} ROM for turbulent incompressible flows.

The proposed PGD \emph{computational vademecum} relies on a separated approximation of the RANS equations coupled with the SA turbulence model. The PGD-ROM methodology mimicks the structure of the \texttt{simpleFoam} algorithm with SA turbulence model, resulting in a minimally intrusive approach within OpenFOAM.

The resulting strategy provides a generalised solution of incompressible flows, explicitly depending on user-defined parameters, in convection-dominated regimes.

The remainder of this paper is organised as follows. Section~\ref{sc:RANS-SA} recalls the full order RANS-SA equations and the corresponding cell-centred FV approximation utilised by OpenFOAM. The rationale of the PGD-ROM for the turbulent incompressible Navier-Stokes equations is introduced in section~\ref{sc:Par-RANS-SA}, where the details of the algorithms to devise the separated velocity-pressure approximation of the flow equations (\texttt{PGD-NS}) and the separated form of the eddy (\texttt{PGD-SA}) and turbulent (\texttt{PGD-$\nu_t$}) viscosities via the SA equation are presented. Numerical experiments involving flow control in external aerodynamics, in two and three dimensions, with Reynolds number ranging up to $1,000,000$ are reported in section~\ref{sc:simulations}, whereas section~\ref{sc:Conclusion} summarises the contributions of this work. Three appendices provide the technical details of the derivation of the PGD formulation of the Navier-Stokes and SA equations, as well as the expressions of the coefficients for the corresponding spatial and parametric iterations of the alternating direction scheme.

\section{The Reynolds-averaged Navier-Stokes equations and the Spalart-Allmaras turbulence model}
\label{sc:RANS-SA}

To simulate turbulent incompressible flows using the RANS equations, the velocity-pressure pair $(\bu,p)$ is decomposed into a mean flow component $(\bU,P)$ and a perturbation $(\bu',p')$, that is $\bu {=} \bU {+} \bu'$ and $p {=} P {+} p'$.

Given an open bounded computational domain $\Omega \subset \mathbb{R}^{d}$ in $d$ spatial dimensions, the boundary $\partial\Omega$ is partitioned such that $\partial\Omega {=} \GammaW \cup \GammaI \cup \GammaO$, where the three disjoint portions $\GammaW$, $\GammaI$ and $\GammaO$ denote material walls, inlet and outlet surfaces, respectively.

The steady-state RANS equations for the mean flow variables $(\bU,P)$ are given by

\begin{equation} \label{eq:RANS}
\left\{\begin{aligned}
\grad {\cdot} (\bU {\otimes} \bU) - \grad {\cdot} ((\nu {+} \nu_t)  \grad \bU) + \grad P &= \bm{0}       &&\text{in $\Omega$,}\\
\grad {\cdot} \bU &= 0  &&\text{in $\Omega$,}\\
\bU &= \bm{0}  &&\text{on $\GammaW$,}\\
\bU &= \bUin  &&\text{on $\GammaI$,}\\
(\nu \grad \bU {-} p \Insd ) \bn &= \bm{0}  &&\text{on $\GammaO$,}
\end{aligned}\right.
\end{equation}

where the velocity profile $\bUin$ is imposed on the inlet surface $\GammaI$, whereas it is customary to use homogeneous Dirichlet and Neumann boundary data to model no-slip and free-traction conditions, respectively, on fixed material walls $\GammaW$ and outlet surfaces $\GammaO$, respectively.

In Equation~\eqref{eq:RANS}, $\Insd$ denotes the $d {\times} d$ identity matrix, $\nu$ represents the physical viscosity of the fluid, whereas the turbulent viscosity $\nu_t$ has been introduced in the momentum equation to model the perturbations to the mean flow due to turbulence.

The one-equation SA turbulence model, derived by means of dimensional analysis and empirical observations~\cite{SpalartAllmaras-92}, introduces the relationship

\begin{equation}\label{eq:nuT}
\nu_t = \nuT f_{v1}
\end{equation}

between the turbulent viscosity $\nu_t$ and the eddy viscosity $\nuT$.

Under the assumption of fully turbulent flows, the trip term in the SA turbulence model is neglected and $\nuT$ is solution of

\begin{equation}\label{eq:SA}
\left\{\begin{aligned}
\grad {\cdot} (\bU \nuT) - \frac{1}{\sigma} \grad {\cdot} ( (\nu {+} \nuT) \grad \nuT ) - \frac{c_{b2}}{\sigma} \grad \nuT {\cdot} \grad \nuT &= s       &&\text{in $\Omega$,}\\
\nuT &= \nuT_D  &&\text{on $\GammaI \cup \GammaW$,}\\
\grad \nuT {\cdot} \bn &= 0  &&\text{on $\GammaO$,}
\end{aligned}\right.
\end{equation}

where $\nuT_D$ is the eddy viscosity datum and the source term is given by

\begin{equation}\label{eq:sourceSA}
s := c_{b1} \widetilde{S} \nuT - c_{w1} \frac{f_w}{\widetilde{d}^2} \nuT^2 .
\end{equation}

The terms on the left-hand side of equation~\eqref{eq:SA} account for the convection and diffusion of the eddy viscosity, whereas the right-hand side, see equation~\eqref{eq:sourceSA}, features the production and destruction contributions, respectively. The SA turbulence model is thus closed by introducing the following definitions

\begin{equation}\label{eq:functionsSA}
\hspace{-7pt}
\begin{alignedat}{5}
\bm{\omega} &:= \frac{\grad \bU - \grad \bU^T}{2} , \quad && \widetilde{S} &&:= \left[ 2\bm{\omega} : \bm{\omega} \right]^{1/2} + \frac{\nuT}{\kappa^2 \widetilde{d}^2} f_{v2} , \quad && \chi &&:= \frac{\nuT}{\nu} , \\
f_w &:= g \left[ \frac{1 + c_{w3}^6}{g^6 + c_{w3}^6} \right]^{1/6} ,  \quad && f_{v2} &&:= 1 - \frac{\chi}{1+ \chi f_{v1}} , \quad && f_{v1} &&:= \frac{\chi^3}{\chi^3 + c_{v1}^3} , \\
c_{w1} &:= \frac{c_{b1}}{\kappa^2} + \frac{1 + c_{b2}}{\sigma} , \quad && g &&:= r + c_{w2}(r^6 -r) , \quad && r &&:= \frac{\nuT}{\widetilde{S} \kappa^2 \widetilde{d}^2} ,
\end{alignedat}
\end{equation}

where $\widetilde{d}$ is the minimum distance to the closest wall, $\sigma {=} 2/3$, $\kappa {=} 0.41$, $c_{b1} {=} 0.1355$, $c_{b2} {=} 0.622$, $c_{v1} {=} 7.1$, $c_{w2} {=} 0.3$ and $c_{w3} {=} 2$.

\subsection{A finite volume formulation of the RANS-SA equations}
\label{sc:CCFV-NS-SA}

In order to discretise the turbulent Navier-Stokes equations, OpenFOAM cell-centred finite volume rationale is considered~\cite{OpenFOAM}. The computational domain is subdivided in $N$ cells $V_i, \, i{=}1,\ldots,N$ such that $V_i {\cap} V_j {=} \emptyset, \text{for $i {\neq} j$}$ and $\Omega {=} \bigcup_{i=1}^{N} V_i$.

In each cell $V_i$, the integral form of equation~\eqref{eq:RANS} is defined as

\begin{equation} \label{eq:weak-RANS}
\!\left\{\begin{aligned}
\int_{V_i}{\!\! \grad {\cdot} (\bU {\otimes} \bU) \, dV}
- \int_{V_i}{\!\!  \grad {\cdot} ((\nu {+} \nu_t)  \grad \bU) \, dV}
+ \int_{V_i}{\!\! \grad P \, dV}
&= \bm{0} , \\
\int_{V_i}{\!\! \grad {\cdot} \bU \, dV} &= 0 ,
\end{aligned}\right.
\end{equation}

where $(\bU,P)$ are cell-by-cell constant approximations of the velocity and pressure fields, respectively, and $\bU {=} \bUin \, \text{on} \, \GammaI$ and $\bU {=} \bm{0} \, \text{on} \, \GammaW$.

In a similar fashion, the cell-centred finite volume approximation of the SA equation~\eqref{eq:SA} is: compute $\nuT$ constant in each cell such that $\nuT = \nuT_D \, \text{on} \, \GammaI \cup \GammaW$ and it holds

\begin{equation} \label{eq:weak-SA}
\begin{aligned}
\int_{V_i}{\!\! \grad {\cdot} (\bU \nuT)  \, dV}
- \frac{1}{\sigma} \int_{V_i}{\!\!  \grad {\cdot} \left( (\nu {+} \nuT) \grad \nuT \right) \, dV}
- \frac{c_{b2}}{\sigma} & \int_{V_i}{\!\!  \grad \nuT {\cdot} \grad \nuT \, dV} \\
=
c_{b1} \int_{V_i}{\!\! \widetilde{S} \nuT \, dV}
& - c_{w1}  \int_{V_i}{\!\! \frac{f_w}{\widetilde{d}^2} \nuT^2 \, dV} .
\end{aligned}
\end{equation}

\subsection{A turbulent Navier-Stokes solver in OpenFOAM}
\label{sc:OpenFOAM-NS-SA}

OpenFOAM strategy to solve the RANS equation with SA turbulence model relies on a staggered approach.

First, the flow equations~\eqref{eq:weak-RANS} are solved using a seed value of $\nu_t$. More precisely, the integrals over each cell in~\eqref{eq:weak-RANS} are approximated by means of the corresponding fluxes across the boundaries of the cell~\cite{Ohlberger-BHO:17,Eymard-EGH:00}. In addition, the semi-implicit method for pressure linked equations (SIMPLE) algorithm~\cite{Patankar-PS:72}, that is, a fractional-step Chorin-Temam projection method~\cite[Sect. 6.7]{Donea-Huerta}, is utilised to handle incompressibility and a relaxation approach is employed for the nonlinear convection term.

Second, the velocity field $\bU$ obtained using \texttt{simpleFoam} is employed to compute the quantities in~\eqref{eq:functionsSA} and to solve the SA equation~\eqref{eq:weak-SA}. It is worth noting that equation~\eqref{eq:weak-SA} is highly nonlinear and a relaxation strategy is implemented in OpenFOAM to improve the convergence of the numerical algorithm~\cite{OpenFOAM}.

Finally, the updated value of the turbulent viscosity $\nu_t$ is determined according to equation~\eqref{eq:nuT} and the \texttt{simpleFoam} routine is utilised to recompute the turbulent velocity and pressure fields.

\section{Proper generalised decomposition for parametric turbulent flow problems}
\label{sc:Par-RANS-SA}

In the context of parametric studies, viscosity coefficient, reference velocity or boundary conditions of the problems may depend on a set of $M$ user-defined parameters $\bmu {=} (\mu_1,\ldots,\mu_M)^T$.

The solution of the RANS-SA equations is thus denoted by the velocity-pressure pair $(\bU(\bx,\bmu),P(\bx,\bmu))$ and the eddy viscosity $\nuT(\bx,\bmu)$, which are now functions of the spatial, $\bx \in \Omega \subset \RR^d$, and parametric, $\bmu\in\bI\subset\mathbb{R}^{M}$, variables.

More precisely, $\bU(\bx,\bmu),P(\bx,\bmu)$ and $\nuT(\bx,\bmu)$ fulfil the high-dimensional RANS-SA equations obtained by integrating~\eqref{eq:weak-RANS}-\eqref{eq:weak-SA} in the parametric space $\bI$.

\subsection{Preliminary notions on the proper generalised decomposition}
\label{sc:PGD}

The PGD separated approximation of the high-dimensional RANS-SA equations, denoted by $(\upgd^n, \ppgd^n)$ for the flow field and $\nuTpgd^m$ for the eddy viscosity, is constructed as a sum of $n$ (respectively, $m$) separable modes, each being the product of functions which depend on either the spatial or one of the parametric variables $\mu_j, \ j=1,\ldots,M$. Henceforth, only space, $\bx$, and parameters, $\bmu$, are separated for the sake of simplicity.

Following~\cite{Tsiolakis-TGSOH-20}, a \emph{predictor-corrector} approximation is utilised, namely

\begin{subequations}\label{eq:sep-increment}
\begin{equation}\label{eq:sep-increment-up}
\left\{\begin{aligned}
 \upgd^n(\bx,\bmu)   &=
   \upgd^{n-1}(\bx,\bmu) + \sigmaU^n\fu^n(\bx)\phi^n(\bmu) + \sigmaU^n\de\upgd^n(\bx,\bmu)  \\
 \ppgd^n(\bx,\bmu)   &=
  \ppgd^{n-1}(\bx,\bmu) +\sigmaP^n\fp^n(\bx)\phi^n(\bmu) + \sigmaP^n\de\ppgd^n(\bx,\bmu)
 \end{aligned}\right.
\end{equation}

for the flow variables and

\begin{equation}\label{eq:sep-increment-nu}
 \nuTpgd^m(\bx,\bmu)  =
  \nuTpgd^{m-1}(\bx,\bmu) +\sigmaNU^m\fnu^m(\bx)\psi^m(\bmu) + \sigmaNU^m\de\nuTpgd^m(\bx,\bmu)
\end{equation}
\end{subequations}

for the eddy viscosity.

It is worth noticing that the same scalar parametric function $\phi(\bmu)$ is selected for both velocity and pressure~\cite{PD-DZH:17}, whereas a different scalar function, $\psi(\bmu)$, is considered for the eddy viscosity.

In~\eqref{eq:sep-increment}, $(\upgd^{n-1},\ppgd^{n-1})$ and $\nuTpgd^{m-1}$ feature the contributions of the previous $n{-}1$ and $m{-}1$ PGD modes, respectively, and $\sigmaU^n\fu^n\phi^n$,$\sigmaP^n\fp^n\phi^n$ and $\sigmaNU^m\fnu^m\psi^m$ represent the predictions of the current term in the PGD expansion of velocity, pressure and eddy viscosity.

The terms $\sigmaU^n\de\upgd^n,\sigmaP^n\de\ppgd^n$ and $\sigmaNU^m\de\nuTpgd^m$ account for the corrections utilised in the computation of the current mode and feature variations $\De$ in the spatial and parametric functions, namely

\begin{equation}\label{eq:PGD-corr}
\left\{\begin{aligned}
\de\upgd^n(\bx,\bmu) &:= \De\fu(\bx)\phi^n(\bmu)+\fu^n(\bx)\De\phi(\bmu)+\De\fu(\bx)\De\phi(\bmu) , \\
\de\ppgd^n(\bx,\bmu) &:= \De\fp(\bx)\phi^n(\bmu)+\fp^n(\bx)\De\phi(\bmu)+\De\fp(\bx)\De\phi(\bmu) , \\
\de\nuTpgd^m(\bx,\bmu) &:= \De\fnu(\bx)\psi^m(\bmu)+\fnu^m(\bx)\De\psi(\bmu)+\De\fnu(\bx)\De\psi(\bmu) ,
\end{aligned}\right.
\end{equation}

where the high-order variation introduced by the last term is henceforth neglected.

\begin{remark}\label{rmrk:amplitude}
The coefficients $\sigmaU^n$, $\sigmaP^n$ and $\sigmaNU^m$ denote the amplitudes of the velocity, pressure and eddy viscosity modes, respectively. Given $\| \phi^n {+} \De\phi \| {=} 1$ and $\| \psi^m {+} \De\psi \| {=} 1$, they are defined as

\begin{equation*}
\sigmaU^n := \| \fu^n + \De\fu \| , \,\, \sigmaP^n := \| \fp^n + \De\fp \| \,\,  \text{and} \,\, \sigmaNU^m := \| \fnu^m + \De\fnu \| ,
\end{equation*}

where the dependence of the modes on $\bx$ and $\bmu$ is omitted for the sake of readability. For the simulations in section~\ref{sc:simulations}, $\eltwo$ norms over the spatial and parametric domains have been considered for the normalisation step.
\end{remark}

It is worth noticing that the final number of modes needed for the PGD approximations, denoted by the super-indexes $n$ and $m$ in~\eqref{eq:sep-increment}, is not known \emph{a priori} and is in general different for the flow variables and the eddy viscosity.

More precisely, the number of terms in the PGD expansion is automatically determined by the algorithm which ends the enrichment procedure when a user-defined stopping criterion is fulfilled. Classical definitions of the stopping criterion include the relative amplitude of the last computed mode with respect to the first one or to the sum of all previosuly computed terms~\cite{Tsiolakis-TGSOH-20}.

\subsection{Proper generalised decomposition of the flow equations}
\label{sc:PGD-RANS}

In this section, the spatial and parametric steps of the non-intrusive PGD algorithm applied to the Navier-Stokes equations~\eqref{eq:weak-RANS} are recalled. First, the laminar case is considered, that is the turbulent viscosity $\nu_t$ is set to zero. The derivation of such formulation is briefly reported in appendix~\ref{sc:appPGD-NS}, whereas the interested reader is referred to~\cite{Tsiolakis-TGSOH-20} for a detailed presentation of the method.

Recall that in order to construct a PGD approximation, a separated representation of the user-defined data is required~\cite{SZ-ZDMH:15}, e.g.  $\nu(\bx,\bmu) {=} D(\bx)\zeta(\bmu)$ for the physical viscosity.

In addition, the first mode $(\upgd^0, \ppgd^0)$ is selected to verify the boundary condition on the inlet surface and, more generally, all inhomogeneous Dirichlet-type boundary conditions. The following terms of the PGD expansion of velocity and pressure are computed via an alternating direction scheme~\cite{PGD-CCH:14, Chinesta-Keunings-Leygue}, henceforth named \texttt{PGD-NS}.

At each iteration of the alternating direction approach, a spatial mode $(\sigmaU^n\fu^n,\sigmaP^n\fp^n)$ is computed using \texttt{simpleFoam} and an algebraic problem is solved to determine the corresponding parametric function $\phi^n$, as detailed in the following subsections.

The alternating direction algorithm stops when the relevance of the increments $\De\fu$ and $\De\fp$ is negligible with respect to the amplitudes $\sigmaU^n$ and $\sigmaP^n$ of the corresponding modes.

Similarly, the global procedure stops when the contribution of the current mode in the PGD enrichment, measured by means of its relative amplitude, is negligible.

\subsubsection{\texttt{PGD-NS}: the spatial iteration}

First, the high-dimensional problem, obtained by integrating~\eqref{eq:weak-RANS} in $\bI$, is restricted to the spatial direction fixing the parametric function $\phi^n$. The spatial increments $(\sigmaU^n \De\fu, \sigmaP^n \De\fp)$ are thus computed as the FV solution of the flow equations

\begin{equation} \label{eq:PGD-NS-spatial}
\left\{\begin{aligned}
\alpha_0 \!\! \int_{V_i}\!\! \grad {\cdot} (\sigmaU^n \De\fu {\otimes} \sigmaU^n & \De\fu)  dV \\
- \alpha_1 \!\! \int_{V_i} \grad {\cdot} (D \grad (\sigmaU^n & \De\fu))  dV
 + \alpha_2 \!\! \int_{V_i}{\grad (\sigmaP^n \De\fp)  dV} \\
&= R_U^n - \!\! \int_{V_i}\!\! \grad {\cdot} \Big( \sum_{j=1}^n \! \alpha_3^j \sigmaU^j \fu^j {\otimes} \sigmaU^{k-1} \De\fu^{k-1} \Big)  dV   \\
& \hspace{4em} - \!\! \int_{V_i}\!\! \grad {\cdot} \Big( \sigmaU^{k-1} \De\fu^{k-1} {\otimes}\!\sum_{j=1}^n \!\alpha_3^j \sigmaU^j \fu^j \Big)   dV , \\
\alpha_2 \!\! \int_{V_i}{\grad {\cdot} (\sigmaU^n \De\fu)  dV} &= R_P^n ,
\end{aligned}\right.
\end{equation}

where the coefficients $\alpha_j, \, j{=}0,\ldots,3$, reported in appendix~\ref{sc:appCoeff}, only depend on user-defined data and parametric functions and can thus be efficiently precomputed.

On the right-hand sides of the momentum and continuity equations, the residuals $R_U^n$ and $R_P^n$ are determined using the previously computed terms $(\upgd^{n-1},\ppgd^{n-1})$ and the predictions $(\sigmaU^n\fu^n\phi^n,\sigmaP^n\fp^n\phi^n)$ of the current mode in the PGD expansions of velocity and pressure , namely

\begin{subequations}\label{eq:NSspatialRes-general}
\begin{align}
&
\begin{aligned}
R_U^n :=
&- \int_{\bI} \phi^n \int_{V_i}{\!\! \grad {\cdot} \left( [\upgd^{n-1} + \sigmaU^n\fu^n \phi^n ] {\otimes} [\upgd^{n-1} + \sigmaU^n\fu^n \phi^n ] \right) dV \, d\bI} \\
&+ \int_{\bI} \phi^n \int_{V_i}{\!\! \grad {\cdot} \left( \nu \grad (\upgd^{n-1} + \sigmaU^n\fu^n \phi^n) \right)  dV \, d\bI} \\
&- \int_{\bI} \phi^n \int_{V_i}{\!\! \grad (\ppgd^{n-1} + \sigmaP^n\fp^n \phi^n)  \, dV \, d\bI} ,
\end{aligned}
\label{eq:NSspatialResU-general} \\

&
R_P^n :=
- \int_{\bI} \phi^n \int_{V_i}{\!\! \grad {\cdot} \left(\upgd^{n-1} +  \sigmaU^n\fu^n \phi^n \right) dV \, d\bI} .
\label{eq:NSspatialResP-general}
\end{align}
\end{subequations}

It is worth recalling that the factor $\phi^n$ in the expressions of $R_U^n$ and $R_P^n$, see equation~\eqref{eq:NSspatialRes-general}, follows from the restriction of the residuals defined in the high-dimensional space $\Omega {\times} \bI$ to the tangent manifold associated with the spatial direction~\cite{Tsiolakis-TGSOH-20}.

In order to perform an efficient computation of such residuals, the separated expressions of $(\upgd^{n-1},\ppgd^{n-1})$ as a product of spatial and parametric functions are exploited, leading to

\begin{subequations}\label{eq:NSspatialRes}
\begin{align}
&
\begin{aligned}
R_U^n
= -& \sum_{j=1}^n \sum_{\ell=1}^n \alpha_4^{j\ell} \!\! \int_{V_i}{\grad {\cdot} (\sigmaU^j\fu^j {\otimes} \sigmaU^\ell\fu^\ell) \, dV} \\
&+ \sum_{\ell=1}^n \alpha_5^\ell \int_{V_i}{\!\! \grad {\cdot} ( D \grad (\sigmaU^\ell\fu^\ell ) ) \, dV}
-   \sum_{\ell=1}^n \alpha_6^\ell \int_{V_i}{\!\! \grad (\sigmaP^\ell\fp^\ell )  \, dV} ,
\end{aligned}
\label{eq:NSspatialResU} \\

&
R_P^n
= - \sum_{\ell=1}^n \alpha_6^\ell \int_{V_i}{\!\! \grad {\cdot} (\sigmaU^\ell \fu^\ell ) \, dV} ,
\label{eq:NSspatialResP}
\end{align}
\end{subequations}

where the coefficients $\alpha_j, \, j{=}4,\ldots,6$ encapsulate the information of the previously computed parametric modes and are defined in appendix~\ref{sc:appCoeff}.

In order to devise a strategy for the spatial iteration of the \texttt{PGD-NS} step which is non-intrusive with respect to the OpenFOAM SIMPLE algorithm, two linear contributions arising from the relaxation of the convection term appear on the right-hand side of the momentum equation in~\eqref{eq:PGD-NS-spatial}. More precisely, the last two integrals are evaluated using the last computed increment $\sigmaU^{k-1}\De\fu^{k-1}$ in the SIMPLE iterations.

It is straightforward to observe that the resulting structure of the left-hand side of equations~\eqref{eq:PGD-NS-spatial} mimicks the traditional Navier-Stokes equations~\eqref{eq:weak-RANS}, with $\nu_t {=} 0$.

Hence, the PGD spatial iteration is solved using the \texttt{simpleFoam} algorithm, natively implemented in OpenFOAM, and a non-intrusive model reduction strategy for parametric problems is obtained~\cite{Tsiolakis-TGSOH-20}.

\subsubsection{\texttt{PGD-NS}: the parametric iteration}

To compute the parametric increment $\De\phi$, the value of the spatial functions \\ $(\sigmaU^n\fu^n,\sigmaP^n\fp^n) {\gets} (\sigmaU^n [\fu^n {+} \De\fu], \sigmaP^n [\fp^n {+} \De\fp])$ is updated and fixed.

The high-dimensional problem is thus restricted to the parametric direction $\bI$ and the following algebraic equation is obtained

\begin{subequations}\label{eq:PGD-NS-param}
\begin{equation}\label{eq:NSparamMatrix}
a_0 (\De\phi)^2 + \left( - a_1 \zeta + a_2 + \sum_{j=1}^n a_3^j \phi^j \right) \De\phi
= r_U^n + r_P^n ,
\end{equation}

where $r_U^n$ and $ r_P^n$ are given by

\begin{equation}\label{eq:NSparamRes}
r_U^n := \sum_{\ell=1}^n \left( - \sum_{j=1}^n a_4^{j\ell} \phi^j + a_5^\ell \zeta - a_6^\ell \right) \phi^\ell , \quad
r_P := - \sum_{\ell=1}^n a_7^\ell \phi^\ell
\end{equation}
\end{subequations}

and represent the residuals of the discretised momentum and continuity equations, respectively, in the parametric space.

\begin{remark}\label{rmrk:diffLaminar}
Contrary to the parametric problem in~\cite{Tsiolakis-TGSOH-20}, the second-order variation is maintained in~\eqref{eq:NSparamMatrix}. Although this term was negligible in laminar simulations, it has been verified numerically that its presence improves the stability of the solution of the parametric step of the Navier-Stokes equations in the turbulent regime.
\end{remark}

The algebraic problem~\eqref{eq:PGD-NS-param} is fully determined by the coefficients $a_i, \ i=0,\ldots,7$, which are detailed in appendix~\ref{sc:appCoeff}. Similarly to what observed in the previous section, these coefficients only depend on user-defined data and on previously computed spatial modes and can be efficiently precomputed. The resulting parametric problem is solved through a collocation approach.

\subsection{Proper generalised decomposition of the turbulence model}
\label{sc:PGD-SA}

In the derivation of the PGD formulation of the flow equations in the previous section, the turbulent viscosity $\nu_t$ was neglected.

In order to devise a PGD strategy for parametric turbulent flow problems, a key aspect is represented by the construction of a separated approximation of $\nu_t$. In this section, such a task is performed via the formulation of spatial and parametric PGD iterations to compute the eddy viscosity $\nuT$ according to the SA turbulence model.

Similarly to the procedure developed for the flow equations, the first mode $\nuTpgd^0$ is arbitrarily selected to match the imposed value of the turbulent viscosity on $\GammaI \cup \GammaW$. More precisely, Dirichlet data for $\nuTpgd$ are selected as full order solutions of the SA equation computed using the boundary condition modes of velocity.

Assuming $\nuTpgd^{m-1}$ known, the $m$-th term in the PGD expansion is computed by means of a greedy approach via an alternating direction method, named \texttt{PGD-SA}.

More precisely, spatial, $\sigmaNU^m\fnu^m$, and parametric, $\psi^m$, modes of the eddy viscosity are determined by alternatively solving a PDE in $\Omega$ and an algebraic problem in $\bI$, until the corrections $\De\fnu$ and $\De\psi$ are negligible with respect to the amplitude $\sigmaNU^m$, as detailed in the following subsections.

Finally, the separated form of $\nu_t$ is retrieved according to relationship~\eqref{eq:nuT}, via the \texttt{PGD-$\nu_t$} routine.

\subsubsection{\texttt{PGD-SA}: the spatial iteration}
\label{sc:PGD-SA-spatial}

Following the rationale presented for the flow equations, the parametric function $\psi^m$ is fixed and the high-dimensional SA equation is restricted to the spatial direction. The increment $\sigmaNU^m\De\fnu$ acts as unknown of the spatial iteration of the PGD procedure for the parametric SA equation. More precisely, a cell-by-cell constant approximation $\sigmaNU^m\De\fnu$ is computed solving

\begin{equation} \label{eq:PGD-SA-spatial}
\hspace{-10pt}
\begin{aligned}
 \int_{V_i}\!\! \grad {\cdot} &\Big(\!\!\Big(\sum_{j=1}^n \beta_1^j \sigmaU^j \fu^j \Big) \sigmaNU^m \De\fnu \Big)  dV \\
- \frac{1}{\sigma} & \int_{V_i}{\!\!  \grad {\cdot} \Big( \Big[ \Big( \beta_2 D {+} \sum_{j=1}^m \beta_3^j \sigmaNU^j \fnu^j \Big) {+} \beta_4 \sigmaNU^m \De\fnu \Big] \grad (\sigmaNU^m \De\fnu) \! \Big) dV} \\
-& \frac{\beta_4 c_{b2}}{\sigma} \int_{V_i}{\!\!  \grad (\sigmaNU^m \De\fnu) {\cdot} \grad (\sigmaNU^m \De\fnu) \, dV} \\
&- \beta_5 c_{b1} \int_{V_i}{\!\! \Sx{\widetilde{S}^m} \sigmaNU^m \De\fnu \, dV}
 + \beta_6 c_{w1}  \int_{V_i}{\!\! \frac{\Sx{f_w^m}}{\widetilde{d}^2} \left(\sigmaNU^m \De\fnu \right)^2 \, dV} \\
&\hspace{85pt} = R_{\nu}^m  + \frac{1}{\sigma} \int_{V_i}{\!\!  \grad {\cdot} \Big( \sigmaNU^{k-1} \De\fnu^{k-1} \grad \Big(\sum_{j=1}^m \beta_3^j \sigmaNU^j \fnu^j \Big) \! \Big) dV} \\
&\hspace{95pt} + \frac{2 c_{b2}}{\sigma} \int_{V_i}{\!\!  \grad \Big(\sum_{j=1}^m \beta_3^j \sigmaNU^j \fnu^j \Big)  {\cdot} \grad (\sigmaNU^{k-1} \De\fnu^{k-1}) \, dV} \\
&\hspace{98pt} - 2 c_{w1}  \int_{V_i}{\!\! \frac{\Sx{f_w^m}}{\widetilde{d}^2} \Big( \sum_{j=1}^m \beta_7^j \sigmaNU^j \fnu^j \Big) \sigmaNU^{k-1} \De\fnu^{k-1} \, dV} ,
\end{aligned}
\end{equation}

where $\widetilde{S}^m$ and $f_w^m$ are evaluated using the $m$ previously computed modes of $\nuTpgd^m$ and $\Sx{\circledcirc}$ denotes the spatial modes of the function $\circledcirc$.

\begin{remark}\label{rmrk:nonSeparable}
It is worth emphasising that neither $\widetilde{S}(\bx,\bmu)$ nor $f_w(\bx,\bmu)$ are separable exactly via an analytical procedure. In this context, in order to efficiently compute the coefficients and the integral terms in~\eqref{eq:PGD-SA-spatial}, a high-dimensional reconstruction of these functions in the space $\Omega \times \bI$ to later interpolate in $\Omega$ or a numerical procedure for PGD separation~\cite{PD-DZGH-18,Diez-DZGH-19} need to be performed. The former strategy is employed in the simulations in section~\ref{sc:simulations}.
\end{remark}

As previously observed for the spatial iteration of the flow equations, the right-hand side of equation~\eqref{eq:PGD-SA-spatial} features the residual $R_{\nu}^m$ obtained using the values of the previous terms $\nuTpgd^{m-1}$ in the PGD expansion of the eddy viscosity and the prediction $\sigmaNU^m\fnu^m\psi^m$ of the current mode, namely

\begin{equation}\label{eq:SAspatialRes-general}
\begin{aligned}
\mathcal{R}_\nu :=
&- \int_{\bI} \psi^m \int_{V_i}{\grad {\cdot} \left(\upgd^n (\nuTpgd^{m-1} + \sigmaNU^m\fnu^m\psi^m) \right) \, dV  \, d\bI} \\
&+ \frac{1}{\sigma} \int_{\bI} \psi^m \int_{V_i}{\!\!  \grad {\cdot} \Big( (\nu + \nuTpgd^{m-1} + \sigmaNU^m\fnu^m\psi^m ) \grad (\nuTpgd^{m-1} + \sigmaNU^m\fnu^m\psi^m) \! \Big) dV  \, d\bI} \\
&+ \frac{c_{b2}}{\sigma} \int_{\bI} \psi^m \int_{V_i}{\!\!  \grad (\nuTpgd^{m-1} + \sigmaNU^m\fnu^m\psi^m) {\cdot} \grad (\nuTpgd^{m-1} + \sigmaNU^m\fnu^m\psi^m) \, dV  \, d\bI} \\
&+ c_{b1} \int_{\bI} \psi^m \int_{V_i}{\!\! \widetilde{S} (\nuTpgd^{m-1} + \sigmaNU^m\fnu^m\psi^m) \, dV  \, d\bI} \\
&- c_{w1}  \int_{\bI} \psi^m \int_{V_i}{\!\! \frac{f_w}{\widetilde{d}^2} (\nuTpgd^{m-1} + \sigmaNU^m\fnu^m\psi^m)^2 \, dV  \, d\bI} ,
\end{aligned}
\end{equation}

where the parametric function $\psi^m$ in the integrals above stems from the restriction of the high-dimensional residuals to the tangent manifold in the spatial direction, see appendix~\ref{sc:appPGD-SA}. By eploiting the separated expression of $\nuTpgd^{m-1}$, the residual can be rewritten as

\begin{equation}\label{eq:SAspatialRes}
\begin{aligned}
R_{\nu}^m
= -& \sum_{j=1}^n \sum_{\ell=1}^m \beta_8^{j\ell} \!\! \int_{V_i}{\grad {\cdot} (\sigmaU^j\fu^j \sigmaNU^\ell\fnu^\ell) \, dV} \\
& + \frac{1}{\sigma} \sum_{\ell=1}^m \int_{V_i}{\!\!  \grad {\cdot} \Big(\!\!\Big( \beta_9^\ell D {+} \sum_{j=1}^m \beta_{10}^{j\ell} \sigmaNU^j \fnu^j \Big) \grad (\sigmaNU^\ell \fnu^\ell) \! \Big) dV} \\
& + \frac{c_{b2}}{\sigma} \sum_{j=1}^m \sum_{\ell=1}^m \beta_{10}^{j\ell} \int_{V_i}{\!\!  \grad (\sigmaNU^j \fnu^j) {\cdot} \grad (\sigmaNU^\ell \fnu^\ell) \, dV} \\
& + c_{b1} \sum_{\ell=1}^m \beta_{11}^\ell \int_{V_i}{\!\! \Sx{\widetilde{S}^m} \sigmaNU^\ell \fnu^\ell \, dV} \\
& - c_{w1}  \sum_{j=1}^m \sum_{\ell=1}^m \beta_{12}^{j\ell} \int_{V_i}{\!\! \frac{\Sx{f_w^m}}{\widetilde{d}^2} \sigmaNU^j \fnu^j \sigmaNU^\ell \fnu^\ell \, dV} .
\end{aligned}
\end{equation}

In order for the discussed PGD-ROM implementation to be non-intrusive with respect to the SA solver natively implemented in OpenFOAM, three terms arising from the restriction of the high-dimensional SA equation to the spatial direction, see appendix~\ref{sc:appPGD-SA}, are treated explicitly via a relaxation approach. Hence, in the last three terms on the right-hand side of equation~\eqref{eq:PGD-SA-spatial}, $\sigmaNU^{k-1} \De\fnu^{k-1}$ denotes the last computed increment in the iterative procedure to solve the SA equation.

The left-hand side of the equation~\eqref{eq:PGD-SA-spatial} thus presents the same structure as the original SA equation~\eqref{eq:weak-SA}, where the spatial integrals are now weighted by means of appropriately computed parametric coefficients, and the OpenFOAM strategy for the solution of the turbulence model equation is exploited.

It is worth noticing that since the functions $\widetilde{S}(\bx,\bmu)$ and $f_w(\bx,\bmu)$ are not separable exactly using an analytical procedure, see remark~\ref{rmrk:nonSeparable}, equation~\eqref{eq:PGD-SA-spatial} cannot be solved in a complete non-intrusive way using OpenFOAM. The resulting implementation of the \texttt{PGD-SA} algorithm in OpenFOAM is thus minimally intrusive as the structure of equation~\eqref{eq:PGD-SA-spatial} is the same as equation~\eqref{eq:weak-SA}  but a numerical separation of $\widetilde{S}(\bx,\bmu)$ and $f_w(\bx,\bmu)$ needs to be integrated in the existing OpenFOAM routine for the SA solver.

The spatial problem~\eqref{eq:PGD-SA-spatial}-\eqref{eq:SAspatialRes} is fully closed by the definitions of the coefficients $\beta_j$, $j{=}1,\ldots,12$ in appendix~\ref{sc:appCoeff}. These quantities depend soley on user-defined data and parametric functions and can be precomputed for an efficient implementation of the PGD spatial iteration.

\subsubsection{\texttt{PGD-SA}: the parametric iteration}
\label{sc:PGD-SA-param}

The parametric iteration for the turbulence model is devised by fixing the newly computed eddy viscosity spatial mode $\sigmaNU^m \fnu^m {\gets} \sigmaNU^m [\fnu^m {+} \De\fnu]$ and restricting the high-dimensional SA equation to the parametric direction $\bI$.

It follows the algebraic equation

\begin{equation}\label{eq:PGD-SA-param}
\begin{aligned}
 (& -b_4 + b_6 \Smu{f_w^m} )(\De\psi)^2 \\
 &+ \left( - b_2 \zeta - b_5 \Smu{\widetilde{S}^m} + \sum_{j=1}^n b_1^j \phi^j - \sum_{j=1}^m (b_3^j + b_7^j \Smu{f_w^m}) \psi^j \right) \De\psi
= r_{\nu}^m ,
\end{aligned}
\end{equation}

in which the unknown is the parametric increment $\De\psi$, discretised in the nodes of the interval $\bI$ according to a collocation method.

The right-hand side of equation~\eqref{eq:PGD-SA-param} features the residual of the SA equation in the parametric space, namely

\begin{equation}\label{eq:SAparamRes}
r_{\nu}^m :=  \sum_{\ell=1}^m \left( b_9^\ell \zeta +  b_{11}^\ell \Smu{\widetilde{S}^m} + \sum_{j=1}^n b_8^{j\ell} \phi^j + \sum_{j=1}^m (b_{10}^{j\ell} - b_{12}^{j\ell} \Smu{f_w^m} ) \psi^j \right) \psi^\ell  .
\end{equation}

As previosuly observed, $\widetilde{S}^m$ and $f_w^m$ are evaluated using the $m$ previously computed modes of $\nuTpgd^m$. Moreover, the functions $\widetilde{S}(\bx,\bmu)$ and $f_w(\bx,\bmu)$ are not separable analytically. Following remark~\ref{rmrk:nonSeparable}, $\Smu{\circledcirc}$ denotes the parametric mode of the non-separable function $\circledcirc$, obtained either via a high-order reconstruction in $\Omega \times \bI$ and consequent interpolation in the parametric space or via the PGD separation described in~\cite{PD-DZGH-18,Diez-DZGH-19}.

The expression of the coefficients $b_j, \ j=1,\ldots,12$ is detailed in appendix~\ref{sc:appCoeff}. As for the parametric step of the flow equations, they only depend on data provided by the user and on spatial functions and they can thus be efficiently precomputed.

\subsubsection{\texttt{PGD-$\nu_t$}: devising a separated turbulent viscosity}

The PGD approximation $\turbNUpgd$ of the turbulent viscosity is obtained introducing the separated form of the eddy viscosity $\nuTpgd$ in~\eqref{eq:nuT}.

It is worth noticing that the function $f_{v1}$ is not separable analytically. Hence, as detailed in remark~\ref{rmrk:nonSeparable}, a high-dimensional reconstruction in the space $\Omega \times \bI$ to perform interpolation in $\Omega$ and $\bI$ separately, or a numerical strategy for PGD separation~\cite{PD-DZGH-18,Diez-DZGH-19} is required. The former strategy is employed in the simulations in section~\ref{sc:simulations}.

For the sake of readability, consider $f_{v1}^m$ obtained using the $m$ previously computed modes of $\nuTpgd^m$ and introduce the approximation

\begin{equation}\label{eq:sepFw1}
f_{v1}^m \simeq \Sx{f_{v1}^m} \Smu{f_{v1}^m} ,
\end{equation}

where $\Sx{f_{v1}^m}$ and $\Smu{f_{v1}^m}$ denote the spatial and parametric modes of the function $f_{v1}^m$, respectively, obtained by means of either the interpolation of its high-dimensional reconstruction or the PGD numerical separation.

The resulting PGD approximation of the turbulent viscosity is

\begin{equation}\label{eq:sep-turbNu}
 \turbNUpgd^q(\bx,\bmu)  =  \turbNUpgd^{q-1}(\bx,\bmu) +\sigmaT^q\ft^q(\bx)\xi^q(\bmu) ,
\end{equation}

where $\ft^q$ and $\xi^q$ denote the spatial and parametric modes, respectively.

\begin{remark}
A computationally efficient implementation of~\eqref{eq:nuT} exploiting the separated nature of all the involved variables can be devised in terms of elemental arithmetic operations of separated functions~\cite{Diez-DZGH-19}. More precisely, the spatial modes $\ft^q$ are obtained from the product of the spatial functions $\Sx{f_{v1}^m}$ and $\fnu^m$, whereas the parametric terms $\xi^q$ stem from the product of the parametric modes $\Smu{f_{v1}^m}$ and $\psi^m$.

It is worth noting that the product of separated functions leads to a number of terms in the PGD expansion larger than the number of modes of its factors. Nonetheless, such operation can be efficiently performed via a separated multiplication~\cite{Diez-DZGH-19} and the result can be compressed to eliminate redundant information~\cite{Modesto-MZH:15,PD-DZGH-18,Diez-DZGH-19}.
\end{remark}

\subsection{A minimally intrusive PGD implementation of a parametric solver for turbulent Navier-Stokes flows in OpenFOAM}
\label{sc:PGD-RANS-SA}

In sections~\ref{sc:PGD-RANS} and~\ref{sc:PGD-SA}, the PGD formulations of the Navier-Stokes equations and the SA turbulence model have been separately presented, introducing the so-called \texttt{PGD-NS}, \texttt{PGD-SA} and \texttt{PGD-$\nu_t$} algorithms.

In order to devise a minimally intrusive parametric solver for turbulent incompressible Navier-Stokes flows in OpenFOAM, these routines are integrated in a unique framework. More precisely, the PGD algorithms for the flow equations and the turbulence model are coupled by mimicking the segregated structure of \texttt{simpleFoam} with SA turbulence model implemented in OpenFOAM, leading to an overall PGD-ROM strategy able to exploit the capabilities of the validated CFD library.

First, recall that the strategy of the full order solver involves the following steps~\cite{OpenFOAM}:

\begin{enumerate}[label=(\Alph*)]
\item Compute velocity and pressure using a fractional step projection approach to solve~\eqref{eq:weak-RANS} with a user-prescribed value of the turbulent viscosity (RANS solver via \texttt{simpleFoam}).
\item Use the value of the computed velocity to solve~\eqref{eq:weak-SA} and determine the eddy viscosity (SA solver).
\item Update the turbulent viscosity according to~\eqref{eq:nuT}.
\item Recompute velocity and pressure solving~\eqref{eq:weak-RANS} with the newly determined turbulent viscosity (RANS solver via \texttt{simpleFoam}).
\end{enumerate}

Following this rationale, the corresponding parametric solver is described in algorithm~\ref{alg:PGD-RANS-SA-OF}.

For step (A), the PGD algorithm solves the parametrised flow equations (Algorithm 1 - Step 5) via the non-intrusive \texttt{PGD-NS} strategy. The initial value of the turbulent viscosity utilised in this computation is selected starting from the boundary conditions modes of the velocity.

\begin{remark}\label{rmrk:RANS-modification}
It is worth emphasising that the contribution of the turbulent viscosity $\nu_t$ in the Navier-Stokes momentum equation needs to be accounted for in step (A) and (D), that is, an additional term is introduced in the PGD formulation for the laminar Navier-Stokes equations presented in section~\ref{sc:PGD-RANS}. More precisely, given the PGD approximation $\turbNUpgd^q$ of the turbulent viscosity obtained from the \texttt{PGD-$\nu_t$} routine, the term

\begin{equation}\label{eq:PGD-NS-turb-viscous-term}
-\int_{\bI} \int_{V_i}{\!\! \grad {\cdot} ((\nu + \turbNUpgd^q) \grad \upgd^n )  dV \, d\bI}
\end{equation}

leads to the contribution

\begin{equation}\label{eq:PGD-NS-turb-LHS-spatial}
- \int_{V_i}{\!\! \grad {\cdot} \Big(\!\! \Big(\alpha_1 D + \sum_{j=1}^q\alpha_7^j \sigmaT^j\ft^j \Big)\grad (\sigmaU^n \De\fu) \! \Big)  dV}
\end{equation}

on the left-hand side of the spatial equation~\eqref{eq:PGD-NS-spatial}, whereas the corresponding term in the spatial residual~\eqref{eq:NSspatialResU} is given by

\begin{equation}\label{eq:PGD-NS-turb-RHS-spatial}
\begin{aligned}
\int_{\bI} \phi^n & \int_{V_i}{\!\! \grad {\cdot} \left( (\nu + \turbNUpgd^q) \grad (\upgd^{n-1} + \sigmaU^n\fu^n \phi^n) \right)  dV \, d\bI} \\
&= \sum_{\ell=1}^n \int_{V_i}{\!\! \grad {\cdot} \Big(\!\! \Big(\alpha_5^\ell D + \sum_{j=1}^q\alpha_8^{j\ell} \sigmaT^j\ft^j\Big) \grad (\sigmaU^\ell\fu^\ell ) \! \Big) \, dV} ,
 \end{aligned}
\end{equation}

where the newly introduced coefficients $\alpha_j, \ j{=}7,8$, depending only on parametric functions, are detailed in appendix~\ref{sc:appCoeff}.

The corresponding terms in the parametric problem are

\begin{equation}\label{eq:PGD-NS-turb-LHS-param}
- \left( a_1 \zeta + \sum_{j=1}^q a_8^j \xi^j \right) \De\phi
\end{equation}

on the left-hand side of equation~\eqref{eq:NSparamMatrix} and

\begin{equation}\label{eq:PGD-NS-turb-RHS-param}
\sum_{\ell=1}^n \left( a_5^\ell \zeta + \sum_{j=1}^q a_9^{j\ell} \xi^j \right) \phi^\ell
\end{equation}

in the definition of $r_U^n$ in~\eqref{eq:NSparamRes}, $a_j, \ j{=}8,9$ being two coefficients depending solely on spatial functions, see appendix~\ref{sc:appCoeff}.
\end{remark}

The PGD enrichment procedure for the RANS equations continues by alternatively solving the spatial and parametric problems with the above modifications, until a user-prescribed threshold is achieved by the amplitude of the computed velocity and pressure modes.

Once a \emph{sufficiently accurate} PGD approximation of velocity and pressure is obtained, step (B) and (C) compute separated representations of the eddy and turbulent viscosities, respectively, by means of the minimally intrusive \texttt{PGD-SA} (Algorithm 1 - Step 9) and \texttt{PGD-$\nu_t$} (Algorithm 1 - Step 12) routines.

Finally, the PGD approximation of the flow equations is reset and a new PGD enrichment is performed using the newly computed PGD approximation of the turbulent viscosity (Step (D)) via the non-intrusive \texttt{PGD-NS} algorithm.

\begin{algorithm}
\caption{An OpenFOAM implementation of turbulent \texttt{pgdFoam}}\label{alg:PGD-RANS-SA-OF}
\begin{algorithmic}[1]
\REQUIRE{Stopping criteria $\eta_\diamond^\star$ and $\eta_\nu^\star$ for the PGD enrichment of the flow equations and the turbulence model. $\diamond=U,P$}
\STATE{Compute boundary condition modes: the spatial mode is solution of~\eqref{eq:weak-RANS} using \texttt{simpleFoam} with SA turbulence model and the parametric modes are set equal to $1$.}
\STATE{Set $n \gets 1$, $i \gets 0$ and initialise the amplitudes $\sigma_\diamond^1 \gets 1$.}

\WHILE{$\sigma_\diamond^n > \eta_\diamond^\star\,\sigma_\diamond^1$}
\STATE{Set the enrichment threshold for turbulent viscosity $\eta_t^i$.}
\STATE{Compute the spatial, $(\sigmaU^n\fu^n,\sigmaP^n\fp^n)$, and parametric, $\phi^n$, modes of velocity and pressure using \texttt{PGD-NS}.}

\IF{$\sigma_\diamond^n < \eta_t^i$}
\STATE{Set $m \gets 1$ and initialise the amplitude $\sigma_\nu^1 \gets 1$.}

\WHILE{$\sigma_\nu^m > \eta_\nu^\star\,\sigma_\nu^1$}
\STATE{Compute the spatial, $\sigmaNU^m\fnu^m$, and parametric, $\psi^m$, modes of the eddy viscosity using \texttt{PGD-SA}.}
\STATE{Update the mode counter: $m \gets m+1$.}
\ENDWHILE

\STATE{Compute the spatial, $\sigmaT^q\ft^q$, and parametric, $\xi^q$, modes of the turbulent viscosity using \texttt{PGD-$\nu_t$}.}
\STATE{Increment viscosity update counter: $i \gets i+1$.}
\STATE{Reinitialise the mode counter: $n \gets 0$.}
\ENDIF

\STATE{Update the mode counter: $n \gets n+1$.}

\ENDWHILE

\end{algorithmic}
\end{algorithm}

\begin{remark}\label{rmrk:turbTolerance}
The overall cost of the parametric solver for the turbulent Navier-Stokes equations depends on the number of updates of the turbulent viscosity.

Effective numerical strategies are devised by tuning the accuracy $\eta_\diamond^\star$ required to stop the enrichment of the velocity and pressure approximations and its relationship with the threshold value $\eta_t^i$ determining the number of updates of the turbulent viscosity.

More precisely, the threshold value $\eta_t^i$ is decreased after each update of the turbulent viscosity to improve the accuracy of the velocity and pressure modes computed by means of \texttt{PGD-NS}.

In the simulations in section~\ref{sc:simulations}, for the $i$-th iteration it holds

\begin{equation}\label{eq:criterionTurbUpdate}
\eta_t^i = 10^{-(i+\gamma)} ,
\end{equation}

that is, starting from an initial accuracy of $10^{-\gamma}$, an exponentially decreasing tolerance is defined for each new update of the turbulent viscosity.

An alternative approach to control the accuracy of the separated representation of eddy and turbulent viscosities may be devised modifying step 6 of algorithm 1 and fixing \emph{a priori} the number of modes in the PGD approximation of the velocity field required to run \texttt{PGD-SA} and $\texttt{PGD-$\nu_t$}$ routines.
\end{remark}

\section{Numerical experiments}
\label{sc:simulations}

In this section, numerical simulations of the NASA wall-mounted hump are presented to demonstrate the potential of the proposed methodology. This problem is a quasi-two-dimensional NASA benchmark devised to validate turbulence modelling, starting from an experimental set-up. The domain consists of a Glauert-Goldschmied type body mounted on a splitter plate between two endplates, see figure~\ref{fig:hump_model}. Following~\cite{Seifert-Pack:Hump-FlowControl,Moin-06}, the characteristic length of the problem is set equal to the chord length of the hump $c {=} \SI{0.42}{m}$, whereas its maximum height is $\SI{0.0537}{m}$ and its span is $\SI{0.5842}{m}$. Flow separation along the hump is controlled via a suction jet acting through an opening located at $65\%$ of the chord $c$, as detailed in figure~\ref{fig:hump_slot}. In the experimental set-up the opening is connected to a plenum, on the bottom of which suction pumps are installed; for the numerical simulations in the following sections, the plenum is removed and the suction effect is imposed as a boundary condition on the opening via a mass flow rate of $1.518\times 10^{-2}\SI{}{kg/s}$ for the jet.

\begin{figure}[ht]
	\centering
	\subfigure[Experimental set-up]{\includegraphics[width=0.49\textwidth]{3dBody}\label{fig:hump_model}}
	\subfigure[Jet location]{\includegraphics[width=0.49\textwidth]{slotLocation}\label{fig:hump_slot}}
	\caption{NASA wall-mounted hump: (a) representation of the experimental set-up and (b) 2D section of the hump, location of the jet (blue rectangle) and detail of the jet opening. Source: {\color{blue}\texttt{https://cfdval2004.larc.nasa.gov/case3.html}}}
	\label{fig:hump_model-hump_slot}
\end{figure}

In the analysis of this problem, the quantity of interest is represented by the effect of the suction jet on the flow separation and on the position of the reattachment point. Experimental and numerical studies~\cite{Greenblatt-PYHSW:Hump,Rumsey:Hump-RANS} verified the quasi-two-dimensional nature of the phenomena identifying minor three-dimensional effects located near the endplates.

Henceforth, the PGD results will be compared to the full order OpenFOAM approximation, considered as \emph{reference solution}.

The NASA wall-mounted hump problem being quasi-two-dimensional, in the following sections both the 2D and the 3D cases are studied. A parametric flow control problem is devised by varying the maximum amplitude of the suction jet between $10\%$ and $100\%$ of the module of a peak velocity $\hat{U}$. In two dimensions, a sinusoidal velocity profile is defined as

\begin{equation}\label{eq:Hump-jetProfile}
\bU^{\text{jet}}_{\!\! \hat{y}} = \mu \frac{\hat{U}}{2}(1-\cos(2\pi \hat{x})) ,
\end{equation}

where $\hat{x}$ is the normalised curvilinear abscissa of the jet patch, that is $\hat{x} \in \left[0 , 1\right]$, and the resulting jet is pointing in the direction $\hat{y}$ orthogonal to the boundary.

Similarly, in the 3D case the jet defined on the plane $(\hat{x},\hat{z})$ and pointing in the orthogonal direction $\hat{y}$ is

\begin{equation}\label{eq:Hump-jetProfile-3D}
\bU^{\text{jet}}_{\!\! \hat{y}} = \mu \frac{\hat{U}}{4}(1-\cos(2\pi \hat{x}))(1-\cos(2\pi \hat{z})),
\end{equation}

where the normalised coordinate $\hat{z}$ is

\begin{equation}\label{eq:Hump-hatZ}
\hat{z}=
\begin{cases}
0 													& \text{ for } z < 0.4 c \\
\displaystyle\frac{5z-2c}{c} 								& \text{ for } 0.4 c \leq z \leq 0.6 c \\
0 													& \text{ for } z > 0.6 c .
\end{cases}
\end{equation}

It is worth noting that the module of the peak velocity $\hat{U}$ is selected such that the ratio between the mass flow rate of the jet and of the inlet is $1.5 \times 10^{-3}$, reproducing the value in the experimental set-up of the NASA wall-mounted hump with a plenum~\cite{Seifert-Pack:Hump-FlowControl}. In addition, both in~\eqref{eq:Hump-jetProfile} and~\eqref{eq:Hump-jetProfile-3D}, the interval of the parametric variable is defined as $\I {=} [0.1,\,1]$.

\subsection{Two-dimensional NASA wall-mounted hump with parametrised jet}
\label{sc:simulation-2D}

The computational domain for the two-dimensional NASA wall-mounted hump is a channel of height $c$, extending $6.39 c$ upstream and $5 c$ downstream as displayed in figure~\ref{fig:hump2D}. The resulting mesh consists of $114,000$ cells.

Homogeneous Dirichlet boundary conditions are imposed on both the velocity and the eddy viscosity on the bottom wall and on the hump. A symmetry condition is imposed on the top wall, whereas on the outlet free traction is enforced. At the inlet, a parabolic profile is imposed for both velocity and eddy viscosity. More precisely, the variables range between a null value on the bottom wall and a maximum value at $y {=} \SI{0.02}{m}$. For the velocity, the peak value is $\SI{34.6}{m/s}$, whereas the free-stream eddy viscosity $\nuT {=} 3\nu$ is selected. The kinematic viscosity is $\nu {=} 1.55274\times 10^{-5}\SI{}{m^2/s}$, thus the resulting Reynolds number is approximately $\text{Re} {=} 936,000$. On the jet patch, the velocity profile~\eqref{eq:Hump-jetProfile} with $\hat{U} {=} \SI{23.4}{m/s}$ is imposed and a homogeneous Neumann condition is considered for the eddy viscosity.

\begin{figure}
	\centering
	\includegraphics[width=0.8\textwidth]{domain2D}
	\caption{Computational domain for the two-dimensional NASA wall-mounted hump.}
	\label{fig:hump2D}
\end{figure}

It is worth noting that on the hump the mesh is such that $y^+<1$, whence no wall treatment is required for the turbulent viscosity.

The boundary conditions are imposed using two modes computed via \texttt{simpleFoam} with SA turbulence model. The first mode is a full order solution with the jet acting at $100\%$ of the mass flow rate ($\mu {=} 1$) and the second one is associated with the jet acting at $10\%$, that is $\mu {=} 0.1$. The corresponding parametric modes are

\begin{equation}\label{eq:parModesHump}
\phi^1(\mu)=\frac{10}{9}(1-\mu) \ \text{ and } \ \phi^2(\mu)= 1 - \phi^1(\mu),
\end{equation}

respectively.

The tolerance for the enrichment of the flow variables is set to $\eta_u^\star {=} \eta_p^\star {=} 10^{-4}$, whereas the tolerance for the turbulence model is selected as $\eta_{\nu}^\star {=} 10^{-2}$. The criterion to update the turbulent viscosity is detailed in remark~\ref{rmrk:turbTolerance}, with $\gamma {=} 1$.

The turbulent \texttt{pgdFoam} algorithm achieves convergence with eight velocity-pressure modes computed using \texttt{PGS-NS} and three corrections by means of \texttt{PGD-SA} and \texttt{PGD-$\nu_t$}. Each \texttt{PGD-SA} loop reached the prescribed tolerance within three computed modes. The relative amplitude of the computed modes, as the turbulent viscosity is updated, is reported in figure~\ref{fig:hump-ampl}.

\begin{figure}
	\centering
	\includegraphics[width=0.6\textwidth]{amplitudeHumpJets2D}
	\caption{Relative amplitude of the computed velocity-pressure modes as the turbulent viscosity is updated.}
	\label{fig:hump-ampl}
\end{figure}

It is worth recalling that each time the relative amplitude of the modes drops by one order of magnitude, the separated representation of the turbulent viscosity is updated via the \texttt{PGD-SA} and \texttt{PGD-$\nu_t$} routines (see remark~\ref{rmrk:turbTolerance}) and the PGD approximation for velocity and pressure is recomputed using the updated turbulent viscosity.

The importance of updating the PGD approximation of the turbulent viscosity to correctly compute the turbulent velocity and pressure fields is displayed in figure~\ref{fig:hump-L2-2D-PGD-SA-comparison}. Considering the result of the full order \texttt{simpleFoam} with SA turbulence model for $\mu {=} 0.5$ as a reference solution, figure~\ref{fig:hump-L2-2D-PGD-SA-comparison} compares the relative $\eltwo(\Omega)$ error of the PGD approximation computed via the \texttt{PGD-NS}, \texttt{PGD-SA} and \texttt{PGD-$\nu_t$} strategy described in algorithm~\ref{alg:PGD-RANS-SA-OF} with the one obtained by omitting the turbulent viscosity update. Without the \texttt{PGD-SA} and \texttt{PGD-$\nu_t$} routines, the error of both velocity and pressure approximations stagnates from the first computed mode and the overall value is one order of magnitude larger than the one achieved by the methodology in algorithm~\ref{alg:PGD-RANS-SA-OF}.

\begin{figure}[ht]
	\centering
	\subfigure[Influence of turbulent viscosity update]{\includegraphics[width=0.49\textwidth]{errorHumpPGD_SA_vs_INTERP_2D}\label{fig:hump-L2-2D-PGD-SA-comparison}}
	\subfigure[Relative $\eltwo(\Omega)$ error]{\includegraphics[width=0.49\textwidth]{errorHumpJets2D}\label{fig:hump-L2-2D}}
	\caption{Accuracy of the PGD approximations of velocity and pressure with respect to the full order solutions as a function of the number of modes. (a) Influence of turbulent viscosity update for the case of $\mu {=} 0.5$. (b) Relative $\eltwo(\Omega)$ error for different values of $\mu$. The vertical lines separates the two boundary condition modes and the computed modes.}
	\label{fig:hump-L2}
\end{figure}

Figure~\ref{fig:hump-L2-2D} reports the relative $\eltwo(\Omega)$ error of the PGD approximation with turbulent viscosity update for three configurations, that is $\mu {=} 0.25$, $\mu {=} 0.5$ and $\mu {=} 0.75$. The results clearly display that the PGD approximation achieves comparable accuracy throughout the parametric interval $\I$ using two boundary condition modes and three computed modes. The following modes only introduce minor corrections to the solution as identified by their corresponding amplitudes, see figure~\ref{fig:hump-ampl}.

As mentioned in the problem statement, the quantities of interest in this study are the position of the reattachment point and the effect of the suction jet on the recirculation bubble.

Figure~\ref{fig:hump-velocity-2D} displays the velocity field after the hump and the recirculation bubble for three values of the parameter $\mu$. The influence of the jet in reducing the flow separation and moving the reattachment point towards the hump is well captured by the PGD approximation which is in excellent agreement with the full order solution.

\begin{figure}[ht]
	\centering
	\begin{tabular}[c]{@{}c@{}c@{ }c@{ }c@{ }}
		$\upgd$ &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{U_attach_stream_PGD_0_25}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{U_attach_stream_PGD_0_5}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{U_attach_stream_PGD_0_75}} \\[2em]
		$\uref$ &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{U_attach_stream_bench_0_25}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{U_attach_stream_bench_0_5}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{U_attach_stream_bench_0_75}} \\
		  & $\mu{=}0.25$ & $\mu{=}0.5$ & $\mu{=}0.75$ \\
	\end{tabular}

	\includegraphics[width=0.6\textwidth]{U_legend}
\caption{Comparison of the PGD approximation (top) and the full order solution (bottom) of the recirculation bubble for $\mu{=}0.25$, $\mu{=}0.5$ and $\mu{=}0.75$. The vertical line denotes the position of the reattachment point.}
\label{fig:hump-velocity-2D}
\end{figure}

Qualitative comparisons of the pressure field and the turbulent viscosity for different values of the parameter $\mu$ are presented in figure~\ref{fig:hump-pressure-2D} and~\ref{fig:hump-viscosity-2D}, respectively. Using eight computed modes, the PGD approximation is able to accurately approximate localised variations in the flow pattern, throughout the interval $\I$.

\begin{figure}[ht]
	\centering
	\begin{tabular}[c]{@{}c@{}c@{ }c@{ }c@{ }}
		$\ppgd$ &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{pISOS_PGD_25}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{pISOS_PGD_50}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{pISOS_PGD_75}} \\[2em]
		$\pref$ &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{pISOS_bench_25}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{pISOS_bench_50}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{pISOS_bench_75}} \\
		  & $\mu{=}0.25$ & $\mu{=}0.5$ & $\mu{=}0.75$ \\
	\end{tabular}

	\includegraphics[width=0.6\textwidth]{p_legend}
\caption{Comparison of the PGD approximation (top) and the full order solution (bottom) of the pressure field after the hump for $\mu{=}0.25$, $\mu{=}0.5$ and $\mu{=}0.75$.}
\label{fig:hump-pressure-2D}
\end{figure}

\begin{figure}[ht]
	\centering
	\begin{tabular}[c]{@{}c@{}c@{ }c@{ }c@{ }}
		$\nu_{t_{\texttt{PGD}}}$ &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{nut_PGD_0_25}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{nut_PGD_0_5}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{nut_PGD_0_75}} \\[2em]
		$\nu_{t_{\texttt{REF}}}$ &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{nut_bench_0_25}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{nut_bench_0_5}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{nut_bench_0_75}} \\
		  & $\mu{=}0.25$ & $\mu{=}0.5$ & $\mu{=}0.75$ \\
	\end{tabular}

	\includegraphics[width=0.6\textwidth]{nut_legend}
\caption{Comparison of the PGD approximation (top) and the full order solution (bottom) of the turbulent viscosity after the hump for $\mu{=}0.25$, $\mu{=}0.5$ and $\mu{=}0.75$.}
\label{fig:hump-viscosity-2D}
\end{figure}

In addition, table~\ref{tab:reattach-2D} confirms the accuracy of the PGD approximation by quantitatively comparing the estimated location of the reattachment point with respect to the position computed using the full order \texttt{simpleFoam} solver with SA turbulence model. The reduced model with eight computed modes provides errors below $0.2\%$.

\begin{table}[ht]
\begin{center}
	\begin{tabular}{|c||c|c|c|}
		\hline
		$ \mu $	   &	$0.25$       &  $0.50$       &  $0.75$      \\
		\hline
		  PGD          &	$1.183 c$  &  $1.156 c$   &  $1.129 c$ \\
		  full order &	$1.184 c$  &  $1.154 c$   &  $1.131 c$ \\
		  \hline
		  Relative error         &  $0.84 \times 10^{-3}$ & $0.17 \times 10^{-2}$ & $0.17 \times 10^{-2}$ \\
		\hline
	\end{tabular}
	\caption{2D NASA wall-mounted hump: position of the reattachment point computed using the PGD approximation and the full order solver for different values of $\mu$.}
	\label{tab:reattach-2D}
\end{center}
\end{table}

Finally, figure~\ref{fig:hump-Cf-2D} and~\ref{fig:hump-Cp-2D} report the skin friction coefficient and the pressure coefficient, respectively. Both figures focus on the area after the jet and compare the PGD approximation using different number of modes with the full order solution. More precisely, the PGD approximations based on the boundary condition modes, i.e. $n {=} 2$, provide results qualitatively comparable with the full order solution, whereas perfect agreement is achieved using eight computed modes ($n {=} 10$).

\begin{figure}[ht]
	\centering
	\subfigure[Skin friction coefficient]{\includegraphics[width=0.49\textwidth]{Cf_hump_2D_scaled}\label{fig:hump-Cf-2D}}
	\subfigure[Pressure coefficient]{\includegraphics[width=0.49\textwidth]{Cp_hump_2D_scaled}\label{fig:hump-Cp-2D}}
	\caption{Comparison of (a) the skin friction coefficient and (b) the pressure coefficient of the full order solution and the PGD approximation, for different number of PGD modes and for the three values of the parameter $\mu$.}
	\label{fig:hump-Cf-Cp}
\end{figure}

\subsection{Three-dimensional NASA wall-mounted hump with parametrised jet}
\label{sc:simulation-3D}

The computational domain for the three-dimensional problem, see figure~\ref{fig:hump-domain-3D}, is obtained by extruding the 2D domain described in the previous section in the $z$ direction by $0.8$ chord lengths.

\begin{figure}
	\centering
	\includegraphics[width=0.8\textwidth]{domain3Dtransparent}
	\caption{Computational domain for the three-dimensional NASA wall-mounted hump.}
	\label{fig:hump-domain-3D}
\end{figure}

The problem inherits the set of boundary conditions utilised in the 2D case. On the additional external surfaces, slip boundary conditions are imposed. The peak value of the inlet velocity is set to $\SI{3.46}{m/s}$ and the profile of the jet suction is defined as in~\eqref{eq:Hump-jetProfile-3D} with $\hat{U} {=} \SI{2.34}{m/s}$. The kinematic viscosity being $\nu {=} 1.55274\times 10^{-5}\SI{}{m^2/s}$, the Reynolds number for the 3D case is $\text{Re} {=} 93,600$ and the computational mesh consists of $2.34$ million cells.

Similarly to the two-dimensional case, the boundary conditions are enforced using the two parametric modes in~\eqref{eq:parModesHump} and two spatial modes corresponding to the \texttt{simpleFoam} solutions with SA turbulence model for $\mu {=} 0.1$ and $\mu {=} 1$.

The values $\eta_u^\star {=} \eta_p^\star {=} 0.5 \times 10^{-3}$ and $\eta_{\nu}^\star {=} 10^{-2}$ are considered for the tolerance of the enrichment loops of the flow variables and the turbulent viscosity, respectively. To reduce the overall cost of the \texttt{PGD-NS}, \texttt{PGD-SA} and \texttt{PGD-$\nu_t$} procedure, the number of turbulent viscosity updates is reduced by considering a lower initial tolerance in criterion~\eqref{eq:criterionTurbUpdate}, namely $\gamma {=} 2$.

Algorithm~\ref{alg:PGD-RANS-SA-OF} achieves convergence with four modes computed by the \texttt{PGS-NS} routine and two \texttt{PGD-SA} and \texttt{PGD-$\nu_t$} corrections. Each \texttt{PGD-SA} loop reached the prescribed tolerance within two computed modes. The PGD approximation is then compared with the corresponding full order solution provided by \texttt{simpleFoam} with the SA turbulence model: the relative $\eltwo(\Omega)$ error for $\mu {=} 0.25$, $\mu {=} 0.5$ and $\mu {=} 0.75$ is displayed in figure~\ref{fig:hump-L2-3D}, reporting that the reduced order model is able to provide errors in velocity and pressure below $0.1 \%$ and $0.5 \%$, respectively.

\begin{figure}[ht]
	\centering
	\includegraphics[width=0.6\textwidth]{errorHumpJets3D}
	\caption{Relative $\eltwo(\Omega)$ error of the PGD approximations of velocity and pressure with respect to the full order solutions as a function of the number of modes for different values of $\mu$. The vertical line separates the two boundary condition modes and the computed modes.}
	\label{fig:hump-L2-3D}
\end{figure}

The effect of the suction jet on flow recirculation is displayed in figure~\ref{fig:hump-Cf-3D}, where a top view of the wall shear stress on the bottom wall is reported starting from the jet patch up to $1.6 c$ downstream. A qualitative comparison between the reduced order and the full order solution confirms the ability of the PGD to accurately reproduce the turbulent flow in the entire range of values $\I$ of the parameter.

\begin{figure}[ht]
	\centering
	\begin{tabular}[c]{@{}c@{}c@{ }c@{ }c@{ }}
		$\btau_{\!\! w_{\texttt{PGD}}}$ &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{TOPVIEW_wallShearStress_PGD_25}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{TOPVIEW_wallShearStress_PGD_50}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{TOPVIEW_wallShearStress_PGD_75}} \\[2em]
		$\btau_{\!\! w_{\texttt{REF}}}$ &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{TOPVIEW_wallShearStress_bench_25}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{TOPVIEW_wallShearStress_bench_50}} &
		\parbox[c]{0.3\textwidth}{\includegraphics[width=0.3\textwidth]{TOPVIEW_wallShearStress_bench_75}} \\
		  & $\mu{=}0.25$ & $\mu{=}0.5$ & $\mu{=}0.75$ \\
	\end{tabular}

	\includegraphics[width=0.6\textwidth]{wallShearStressScale}
\caption{Comparison of the PGD approximation (top) and the full order solution (bottom) of the wall shear stress on the bottom wall for $\mu{=}0.25$, $\mu{=}0.5$ and $\mu{=}0.75$. Detail of the region starting from the jet patch up to $1.6 c$ downstream.}
	\label{fig:hump-Cf-3D}
\end{figure}

In addition, figure~\ref{fig:hump-velocity-3D} displays the velocity profile on the hump, computed using the PGD, for different values of the parameter $\mu$, whereas figure~\ref{fig:hump-streamlines-3D} reports the corresponding  streamlines. The results display that the recirculation effects are reduced when increasing the suction jet and the PGD is able to capture the vortex structure with comparable accuracy with respect to the full order solution.

\begin{figure}[ht]
	\centering
	\subfigure[$\mu{=}0.25$]{\includegraphics[width=0.49\textwidth]{velocityPGD_25Z}}
	\subfigure[$\mu{=}0.75$]{\includegraphics[width=0.49\textwidth]{velocityPGD_75Z}}
	\caption{Comparison of the PGD approximation of the velocity profile on the hump for $\mu{=}0.25$ and $\mu{=}0.75$.}
	\label{fig:hump-velocity-3D}
\end{figure}

\begin{figure}[ht]
	\centering
	\begin{tabular}[c]{@{}c@{}c@{ }c@{ }}
		$\upgd$ &
		\parbox[c]{0.45\textwidth}{\includegraphics[width=0.45\textwidth]{Umin3D_streamsZ}} &
		\parbox[c]{0.45\textwidth}{\includegraphics[width=0.45\textwidth]{Umax3D_streamsZ}} \\[2.5em]
		$\uref$ &
		\parbox[c]{0.45\textwidth}{\includegraphics[width=0.45\textwidth]{UbenchMin3D_streamsZ}} &
		\parbox[c]{0.45\textwidth}{\includegraphics[width=0.45\textwidth]{UbenchMax3D_streamsZ}} \\
		  & $\mu{=}0.25$ & $\mu{=}0.75$ \\
	\end{tabular}
\caption{Detail of the vortex structure in the recirculation region computed using the PGD approximation (top) and the full order solution (bottom) for $\mu{=}0.25$ and $\mu{=}0.75$.}
	\label{fig:hump-streamlines-3D}
\end{figure}

Finally, the position of the reattachment point in correspondance of the location of the peak of the jet profile is reported in table~\ref{tab:reattach-3D} for different values of the parameter $\mu$. The PGD approximation with four computed modes is in excellent agreement with the full order solver, with relative errors below $0.5\%$.

\begin{table}[ht]
\begin{center}
	\begin{tabular}{|c||c|c|c|}
		\hline
		$ \mu $	   &	$0.25$       &  $0.50$       &  $0.75$      \\
		\hline
		  PGD          &	$1.102 c$  &  $1.062 c$   &  $1.024 c$ \\
		  full order  &	$1.103 c$  &  $1.059 c$   &  $1.019 c$ \\
		  \hline
		  Relative error         &  $0.91 \times 10^{-3}$ & $0.28 \times 10^{-2}$ & $0.49 \times 10^{-2}$ \\
		\hline
	\end{tabular}
	\caption{3D NASA wall-mounted hump: position of the reattachment point computed using the PGD approximation and the full order solver for different values of $\mu$.}
	\label{tab:reattach-3D}
\end{center}
\end{table}

\section{Conclusion}
\label{sc:Conclusion}

A PGD strategy to compute parametric solutions of turbulent incompressible flow problems in OpenFOAM has been proposed. The methodology is based on the incompressible Reynolds-averaged Navier-Stokes equations with Spalart-Allmaras turbulence model and mimicks the segragated approach implemented in the industrially-validated solver OpenFOAM to devise a minimally intrusive PGD-ROM for convection-dominated flow problems of industrial interest.

First, the velocity and pressure modes are computed using the non-intrusive PGD strategy \texttt{PGD-NS} developed in~\cite{Tsiolakis-TGSOH-20} using a seed value for the turbulent viscosity. The PGD approximation of the velocity is then used to improve the turbulent viscosity representation via the minimally intrusive \texttt{PGD-SA} and \texttt{PGD-$\nu_t$} routines. Finally, the resulting separated turbulent viscosity is utilised to recompute the PGD expansions of velocity and pressure.

The importance of an accurate approximation of the turbulent viscosity has been verified by comparing the solution of the above algorithm with the one computed without solving the SA equation: the latter solution quickly stagnates providing errors of one order of magnitude larger than the proposed methodology.

The developed strategy has been validated in two and three spatial dimensions using a benchmark problem of turbulent external flow, the NASA wall-mounted hump, with $\text{Re} {=} 93,600$ and $\text{Re} {=} 936,000$. A flow control problem of industrial interest has been devised by introducing a suction jet on the hump to reduce the recirculation effects.

The proposed PGD-based reduced order model has proved to be able to compute a reduced basis with no \emph{a priori} knowledge of the solution, for convection-dominated viscous incompressible flows achieving both qualitative and quantitative agreement with the full order solution computed via \texttt{simpleFoam} with SA turbulence model, throughout the interval of the parametric variable.

More precisely, the reduced model provided accurate approximations of the velocity and pressure fields, with relative $\eltwo$ errors below $0.1 \%$ and $1 \%$, respectively. In addition, it proved to be able to capture localised flow features and estimate quantities of engineering interest such the position of the reattachment point, the skin friction coefficient and the pressure coefficient with errors below $0.5 \%$.

\section*{Acknowledgements}

This work was partially supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sk\l odowska-Curie Actions (Grant agreement No. 675919) that financed the Ph.D. fellowship of the first author.

The second, third and last authors were supported by the Spanish Ministry of Economy and Competitiveness (Grant agreement No. DPI2017-85139-C2-2-R).

The second and last authors are also grateful for the financial support provided by the Spanish Ministry of Economy and Competitiveness through the Severo Ochoa programme for centres of excellence in RTD (Grant agreement No. CEX2018-000797-S) and the Generalitat de Catalunya (Grant agreement No. 2017-SGR-1278).

\bibliographystyle{elsarticle-num}
\bibliography{Ref-PGD}

\appendix

\section{PGD spatial and parametric coefficients}
\label{sc:appCoeff}

\subsection*{PGD for the Navier-Stokes equations}

The coefficients for the spatial iteration~\eqref{eq:PGD-NS-spatial} are

\begin{equation}\label{eq:NScoef-spatial}
\begin{alignedat}{5}
  \alpha_0 & := \int_{\bI}{ \left[\phi^n\right]^3 \, d\bI} , \quad
  && \alpha_1 && := \int_{\bI}{ \left[\phi^n\right]^2 \zeta \, d\bI} , \quad
  && \alpha_2 && := \int_{\bI}{ \left[\phi^n\right]^2 \, d\bI} , \\
  \alpha_3^j & := \int_{\bI}{ [\phi^n]^2 \phi^j \, d\bI}, \quad
  && \alpha_4^{j\ell} && := \int_{\bI}{ \phi^n \phi^j \phi^\ell \, d\bI} , \quad
  && \alpha_5^\ell && := \int_{\bI}{ \phi^n \phi^\ell \zeta \, d\bI} , \\
  \alpha_6^\ell & := \int_{\bI}{ \phi^n \phi^\ell \, d\bI} , \quad
  && \alpha_7^j && := \int_{\bI}{ \left[\phi^n\right]^2 \xi^j \, d\bI} , \quad
  && \alpha_8^{j\ell} && := \int_{\bI}{ \phi^n \phi^\ell \xi^j \, d\bI} .
\end{alignedat}
\end{equation}

The coefficients for the parametric iteration~\eqref{eq:PGD-NS-param} are

\begin{equation}\label{eq:NScoef-param}
\begin{aligned}
  a_0 &:= \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \bigl[\grad {\cdot} (\sigmaU^n\fu^n {\otimes} \sigmaU^n\fu^n)\bigr] \, dV} , \\
  a_1 &:=  \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \bigl[\grad {\cdot} (D \grad (\sigmaU^n\fu^n))\bigr] \, dV} , \\
  a_2 &:= \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \grad (\sigmaP^n\fp^n)  \, dV}
			 + \!\! \int_{V_i}{\!\!\! \sigmaP^n\fp^n \grad {\cdot} (\sigmaU^n\fu^n) \, dV} , \\
  a_3^j &:= \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \bigl[\grad {\cdot} (\sigmaU^n\fu^n {\otimes} \sigmaU^j\fu^j) + \grad {\cdot} (\sigmaU^j\fu^j {\otimes} \sigmaU^n\fu^n)\bigr] \, dV} , \\
  a_4^{j\ell} &:= \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \bigl[\grad {\cdot} (\sigmaU^j\fu^j {\otimes} \sigmaU^\ell\fu^\ell)\bigr] \, dV} , \\
  a_5^\ell &:= \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \bigl[\grad {\cdot} (D \grad (\sigmaU^\ell\fu^\ell))\bigr] \, dV} , \\
  a_6^\ell &:= \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \grad (\sigmaP^\ell\fp^\ell)  \, dV} , \\
  a_7^\ell &:=  \int_{V_i}{\!\!\! \sigmaP^n\fp^n \grad {\cdot} (\sigmaU^\ell\fu^\ell) \, dV}  , \\
  a_8^j &:=  \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \bigl[\grad {\cdot} (\sigmaT^j\ft^j \grad (\sigmaU^n\fu^n))\bigr] \, dV} , \\
  a_9^{j\ell} &:= \int_{V_i}{\!\!\! \sigmaU^n\fu^n {\cdot} \bigl[\grad {\cdot} (\sigmaT^j\ft^j \grad (\sigmaU^\ell\fu^\ell))\bigr] \, dV} .
  \end{aligned}
\end{equation}

\subsection*{PGD for the Spalart-Allmaras equations}

The coefficients for the spatial iteration~\eqref{eq:PGD-SA-spatial} are

\begin{equation}\label{eq:SAcoef-spatial}
\begin{alignedat}{3}
  \beta_1^j & := \int_{\bI}{ \left[\psi^m\right]^2 \phi^j \, d\bI} , \quad
  && \beta_2 && := \int_{\bI}{ \left[\psi^m\right]^2 \zeta \, d\bI} , \\
  \beta_3^j & := \int_{\bI}{ \left[\psi^m\right]^2 \psi^j\, d\bI} , \quad
  && \beta_4 && := \int_{\bI}{ \left[\psi^m\right]^3 \, d\bI} , \\
  \beta_5 & := \int_{\bI}{ \left[\psi^m\right]^2 \Smu{\widetilde{S}^m} \, d\bI} , \quad
  && \beta_6 && := \int_{\bI}{ \left[\psi^m\right]^3 \Smu{f_w^m} \, d\bI} , \\
  \beta_7^j & := \int_{\bI}{ [\psi^m]^2 \psi^j \Smu{f_w^m} \, d\bI} , \quad
  && \beta_8^{j\ell} && := \int_{\bI}{ \psi^m \psi^\ell \phi^j \, d\bI} , \\
  \beta_9^\ell & := \int_{\bI}{ \psi^m \psi^\ell \zeta \, d\bI} , \quad
  && \beta_{10}^{j\ell} && := \int_{\bI}{ \psi^m \psi^\ell \psi^j \, d\bI} , \\
 \beta_{11}^\ell & := \int_{\bI}{ \psi^m \psi^\ell \Smu{\widetilde{S}^m} \, d\bI} , \quad
  && \beta_{12}^{j\ell} && := \int_{\bI}{ \psi^m \psi^\ell \psi^j \Smu{f_w^m} \, d\bI} .
\end{alignedat}
\end{equation}

The coefficients for the parametric iteration~\eqref{eq:PGD-SA-param} are

\begin{equation}\label{eq:SAcoef-param}
\begin{aligned}
  b_1^j :=& \, \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad {\cdot} (\sigmaU^j\fu^j \sigmaNU^m\fnu^m) \bigr] \, dV} , \\
  b_2 :=& \, \frac{1}{\sigma} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad {\cdot} (D \grad (\sigmaNU^m \fnu^m)) \bigr]  \, dV} , \\
  b_3^j :=& \, \frac{1}{\sigma} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad {\cdot} (\sigmaNU^j \fnu^j \grad (\sigmaNU^m \fnu^m))
  + \grad {\cdot} (\sigmaNU^m \fnu^m \grad (\sigmaNU^j \fnu^j)) \bigr]  \, dV} \\
  &+ \frac{2 c_{b2}}{\sigma} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad (\sigmaNU^j \fnu^j ) {\cdot} \grad (\sigmaNU^m \fnu^m) \bigr]  \, dV} , \\
  b_4 :=& \, \frac{1}{\sigma} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad {\cdot} (\sigmaNU^m \fnu^m \grad (\sigmaNU^m \fnu^m)) \bigr]  \, dV} \\
  &+ \frac{c_{b2}}{\sigma} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad (\sigmaNU^m \fnu^m ) {\cdot} \grad (\sigmaNU^m \fnu^m) \bigr]  \, dV} , \\
  b_5 :=& \, c_{b1} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\Sx{\widetilde{S}^m} \sigmaNU^m\fnu^m\bigr] \, dV} , \\
  b_6 :=& \, c_{w1} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \Bigl[\frac{\Sx{f_w^m}}{\widetilde{d}^2} (\sigmaNU^m\fnu^m)^2 \Bigr] \, dV} , \\
  b_7^j :=& \, 2 c_{w1} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \Bigl[\frac{\Sx{f_w^m}}{\widetilde{d}^2} \sigmaNU^j\fnu^j\sigmaNU^m\fnu^m \Bigr] \, dV} , \\
  b_8^{j\ell} :=& \, \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad {\cdot} (\sigmaU^j\fu^j \sigmaNU^\ell\fnu^\ell) \bigr] \, dV} , \\
  b_9^\ell :=& \, \frac{1}{\sigma} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad {\cdot} (D \grad (\sigmaNU^\ell \fnu^\ell)) \bigr]  \, dV} , \\
  b_{10}^{j\ell} :=& \, \frac{1}{\sigma} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad {\cdot} (\sigmaNU^j \fnu^j \grad (\sigmaNU^\ell \fnu^\ell)) \bigr]  \, dV} \\
  &+ \frac{c_{b2}}{\sigma} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\grad (\sigmaNU^j \fnu^j ) {\cdot} \grad (\sigmaNU^\ell \fnu^\ell) \bigr]  \, dV} , \\
  b_{11}^\ell :=& \, c_{b1} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \bigl[\Sx{\widetilde{S}^m} \sigmaNU^\ell\fnu^\ell\bigr] \, dV} , \\
  b_{12}^{j\ell} :=& \, c_{w1} \int_{V_i}{\!\!\! \sigmaNU^m\fnu^m \Bigl[\frac{\Sx{f_w^m}}{\widetilde{d}^2} \sigmaNU^j\fnu^j\sigmaNU^\ell\fnu^\ell \Bigr] \, dV} .
  \end{aligned}
\end{equation}

\section{Derivation of the proper generalised decomposition of the Navier-Stokes equations}
\label{sc:appPGD-NS}

Consider the velocity-pressure pair $(\bU(\bx,\bmu),P(\bx,\bmu))$ defined in the high-dimensional space $\Omega \times \bI$ as solution of the steady laminar parametrised Navier-Stokes equations, namely

\begin{equation} \label{eq:weak-NS-PGD}
\left\{\begin{aligned}
\begin{aligned}[b]
\int_{\bI}\int_{V_i}\!\! & \grad {\cdot} (\bU {\otimes} \bU) \, dV \, d\bI \\
&- \int_{\bI}\int_{V_i} \!\! \grad {\cdot} (\nu \grad \bU) \, dV \, d\bI
+\int_{\bI}\int_{V_i}{\!\! \grad P  \, dV \, d\bI}
\end{aligned}
&= \bm{0} , \\
\int_{\bI}\int_{V_i}{\!\! \grad {\cdot} \bU \, dV \, d\bI} &= 0 ,
\end{aligned}\right.
\end{equation}

in each cell $V_i, \, i{=}1,\ldots,N$ of the spatial computational domain.

According to the predictor-corrector approximation rationale introduced in~\eqref{eq:sep-increment}, the PGD approximation is splitted into the unknown increments $(\sigmaU^n\de\upgd^n,\sigmaP^n\de\ppgd^n)$ and the known functions obtained from the $n-1$ previous modes $(\upgd^{n-1},\ppgd^{n-1})$ and the predictions $(\sigmaU^n\fu^n\phi^n,\sigmaP^n\fp^n\phi^n)$ of the $n$-th term.
In addition, define the estimated rank-$n$ PGD approximation for velocity and pressure as the sum of the converged $n{-}1$ terms and the prediction of the current mode, that is, $\uBpgd^n {:=} \upgd^{n-1} + \sigmaU^n\fu^n\phi^n$ and $\pBpgd^n {:=} \ppgd^{n-1} +\sigmaP^n\fp^n\phi^n$.

Thus, from equation~\eqref{eq:weak-NS-PGD}, it follows

\begin{equation} \label{eq:PGD-NS-increment}
\left\{\begin{aligned}
\int_{\bI}\int_{V_i}{\!\! \grad {\cdot} (\sigmaU^n\de\upgd^n {\otimes} \sigmaU^n\de\upgd^n) \, dV \, d\bI}
\hspace{182pt} & \\
+ \int_{\bI}\int_{V_i}{\!\! \grad {\cdot} (\sigmaU^n\de\upgd^n {\otimes} \upgd^n) \, dV \, d\bI}
+ \int_{\bI}\int_{V_i}{\!\! \grad {\cdot} (\upgd^n {\otimes} \sigmaU^n\de\upgd^n) \, dV \, d\bI}
\hspace{18pt} & \\
- \int_{\bI}{ \zeta \int_{V_i} \!\! \grad {\cdot} (D \grad (\sigmaU^n\de\upgd^n))  \, dV \, d\bI}
+ \int_{\bI}\int_{V_i}{\!\! \grad (\sigmaP^n\de\ppgd^n)  \, dV \, d\bI} = \mathcal{R}_u ,
 & \\
\int_{\bI}\int_{V_i}{\!\! \grad {\cdot} (\sigmaU^n\de\upgd^n) \, dV \, d\bI} = \mathcal{R}_p ,
&

\end{aligned}\right.
\end{equation}

where the unknowns $(\sigmaU^n\de\upgd^n,\sigmaP^n\de\ppgd^n)$ have been gathered on the left-hand side, whereas the right-hand side features the following residuals

\begin{equation}\label{eq:resNS}
\begin{aligned}
\mathcal{R}_u &:=
\begin{aligned}[t]
&- \int_{\bI}\int_{V_i}{\!\! \grad {\cdot} (\uBpgd^n {\otimes} \uBpgd^n) \, dV \, d\bI} \\
&+ \int_{\bI} \int_{V_i}{\!\! \grad {\cdot} (\nu \grad \uBpgd^n)  \, dV \, d\bI}
- \int_{\bI}\int_{V_i}{\!\! \grad \pBpgd^n  \, dV \, d\bI} ,
\end{aligned}
\\
\mathcal{R}_p &:= - \int_{\bI}\int_{V_i}{\!\! \grad {\cdot} \uBpgd^n \, dV \, d\bI} .
\end{aligned}
\end{equation}

In order to devise the alternating direction algorithm described in section~\ref{sc:PGD-RANS}, equation~\eqref{eq:PGD-NS-increment} is alternatively restricted to the tangent manifold associated with either the spatial or the parametric coordinate.

It is worth noting that the terms in equation~\eqref{eq:PGD-NS-increment} are multiplied by $\phi^n$ and $(\sigmaU^n\fu^n,\sigmaP^n\fp^n)$, when restricted to the spatial and parametric direction, respectively.

Hence, setting the value of the parametric function $\phi^n$, the increments $(\sigmaU^n \De\fu, \sigmaP^n \De\fp)$ are computed by solving the spatial PDE~\eqref{eq:PGD-NS-spatial}. Similarly, the algebraic parametric problem~\eqref{eq:PGD-NS-param} provides the increment $\De\phi$, when the spatial functions $(\sigmaU^n \fu^n,\sigmaP^n \fp^n)$ are assumed to be known. Interested readers are referred to~\cite{Tsiolakis-TGSOH-20} for a detailed explanation of such derivation.

An efficient implementation of this procedure may be devised by exploiting the affine decomposition of the forms in~\eqref{eq:PGD-NS-increment} and~\eqref{eq:resNS}, see~\cite{Patera-Rozza:07,Rozza:14}, and the separated structure of the unknowns, see~\eqref{eq:sep-increment}.

\begin{remark}\label{rmrk:relaxNS}
It is straightforward to observe that equations~\eqref{eq:PGD-NS-increment} mimick the structure of the full order Navier-Stokes equations~\eqref{eq:weak-NS-PGD}. In order to devise a PGD algorithm non-intrusive with respect to OpenFOAM spatial solver \texttt{simpleFOAM}, in~\cite{Tsiolakis-TGSOH-20} the authors proposed to relax the second and third terms in~\eqref{eq:PGD-NS-increment} during the PGD spatial step, evaluating them during the last computed SIMPLE iteration. Such approach is also followed in section~\ref{sc:PGD-RANS} of the present work.
\end{remark}

\section{Derivation of the proper generalised decomposition of the Spalart-Allmaras equation}
\label{sc:appPGD-SA}

In each cell $V_i, \, i{=}1,\ldots,N$ of the computational domain, the parametric version of the steady Spalart-Allmaras turbulence model is given by

\begin{equation} \label{eq:weak-SA-PGD}
\begin{aligned}
\int_{\bI}\int_{V_i} \!\! \grad {\cdot}(\bU \nuT) & \, dV  \, d\bI
- \frac{1}{\sigma} \int_{\bI}\int_{V_i}{\!\!  \grad {\cdot} \left( (\nu {+} \nuT) \grad \nuT \right) \, dV  \, d\bI} \\
-& \frac{c_{b2}}{\sigma} \int_{\bI}\int_{V_i}{\!\!  \grad \nuT {\cdot} \grad \nuT \, dV  \, d\bI}
- c_{b1} \int_{\bI}\int_{V_i}{\!\! \widetilde{S} \nuT \, dV  \, d\bI} \\
&\hspace{90pt} + c_{w1}  \int_{\bI}\int_{V_i}{\!\! f_w \left( \frac{\nuT}{\widetilde{d}} \right)^2 \, dV  \, d\bI} = 0 ,
\end{aligned}
\end{equation}

where the eddy viscosity $\nuT(\bx,\bmu)$ is defined in the high-dimensional space $\Omega \times \bI$.

First, replace the velocity field $\bU$ by its rank-$n$ PGD approximation $\upgd^n$ constructed using the procedure described in section~\ref{sc:PGD-RANS}.

Following the rationale utilised for the flow equations, the PGD approximation of the eddy viscosity in~\eqref{eq:sep-increment-nu} is introduced into~\eqref{eq:weak-SA-PGD} and the unknown increment $\sigmaNU^m\de\nuTpgd^m$ is gathered on the left-hand side of the equation, namely,

\begin{equation}\label{eq:PGD-SA-increment}
\begin{aligned}
\int_{\bI}\int_{V_i}\!\! \grad {\cdot} \Big( \upgd^n \sigmaNU^m\de\nuTpgd^m \Big)  dV  \, d\bI \hspace{12em} & \\
- \frac{1}{\sigma} \int_{\bI}\int_{V_i}{\!\!  \grad {\cdot} \Big( (D\zeta + \nuBpgd^m + \sigmaNU^m\de\nuTpgd^m) \grad (\sigmaNU^m\de\nuTpgd^m) \! \Big) dV  \, d\bI} & \\
- \frac{1}{\sigma} \int_{\bI}\int_{V_i}{\!\!  \grad {\cdot} \Big( \sigmaNU^m\de\nuTpgd^m \grad \nuBpgd^m \! \Big) dV  \, d\bI} \hspace{7em} & \\
- \frac{c_{b2}}{\sigma} \int_{\bI}\int_{V_i}{\!\!  \grad (\sigmaNU^m\de\nuTpgd^m) {\cdot} \grad (\sigmaNU^m\de\nuTpgd^m) \, dV  \, d\bI} \hspace{3em} & \\
- \frac{2 c_{b2}}{\sigma} \int_{\bI}\int_{V_i}{\!\!  \grad \nuBpgd^m  {\cdot} \grad (\sigmaNU^m\de\nuTpgd^m) \, dV  \, d\bI} \hspace{4em} & \\
- c_{b1} \int_{\bI}\int_{V_i}{\!\! \widetilde{S} \sigmaNU^m\de\nuTpgd^m \, dV  \, d\bI} \hspace{7em} & \\
+ c_{w1}  \int_{\bI}\int_{V_i}{\!\! \frac{f_w}{\widetilde{d}^2} \left(\sigmaNU^m\de\nuTpgd^m \right)^2 \, dV  \, d\bI} \hspace{3em} & \\
+ 2 c_{w1}  \int_{\bI}\int_{V_i}{\!\! \frac{\Sx{f_w^m}}{\widetilde{d}^2} \nuBpgd^m \sigmaNU^m\de\nuTpgd^m \, dV  \, d\bI} &=\mathcal{R}_\nu ,
\end{aligned}
\end{equation}

where $\nuBpgd^m {:=} \nuTpgd^{m-1} +\sigmaNU^m\fnu^m\psi^m$ denotes the estimated rank-$m$ approximation of the eddy viscosity obtained using the $m{-}1$ previously converged modes $\nuTpgd^{m-1}$ and the prediction $\sigmaNU^m\fnu^m\psi^m$ of the $m$-th term.

Moreover, the residual on the right-hand side is defined as

\begin{equation}\label{eq:resSA}
\begin{aligned}
\mathcal{R}_\nu :=
&- \int_{\bI}\int_{V_i}{\grad {\cdot} (\upgd^n \nuBpgd^m) \, dV  \, d\bI} \\
&+ \frac{1}{\sigma} \int_{\bI}\int_{V_i}{\!\!  \grad {\cdot} \Big( (\nu +\nuBpgd^m ) \grad \nuBpgd^m \! \Big) dV  \, d\bI} \\
&+ \frac{c_{b2}}{\sigma} \int_{\bI}\int_{V_i}{\!\!  \grad \nuBpgd^m {\cdot} \grad \nuBpgd^m \, dV  \, d\bI} \\
&+ c_{b1} \int_{\bI}\int_{V_i}{\!\! \widetilde{S} \nuBpgd^m \, dV  \, d\bI}
 - c_{w1}  \int_{\bI}\int_{V_i}{\!\! \frac{f_w}{\widetilde{d}^2} (\nuBpgd^m)^2 \, dV  \, d\bI} .
\end{aligned}
\end{equation}

Problem~\eqref{eq:PGD-SA-increment} is solved by means of a greedy procedure relying on an alternating direction algorithm. More precisely, the equation is alternatively restricted to the tangent manifold associated with either the spatial or the parametric coordinate.

Equation~\eqref{eq:PGD-SA-increment} and the residual~\eqref{eq:resSA} are thus multiplied by the parametric function $\psi^m$ and by the spatial function $\sigmaNU^m\fnu^m$ during the spatial and parametric iteration, respectively.

In the former step, the value of the parametric function $\psi^m$ is fixed and the increment $\sigmaNU^m \De\fnu$ is computed by solving the spatial PDE~\eqref{eq:PGD-SA-spatial}. In the latter, the spatial function $\sigmaNU^m \fnu^m$ is set and the algebraic parametric problem~\eqref{eq:PGD-SA-param} is solved to determine the increment $\De\psi$.

The steps of the PGD algorithm for the parametric SA turbulence model are detailed in section~\ref{sc:PGD-SA}.

\begin{remark}
Similarly to the observation in remark~\ref{rmrk:relaxNS}, the first, second, fourth, sixth and seventh term in equation~\eqref{eq:PGD-SA-increment} mimick the original SA turbulence model~\eqref{eq:weak-SA-PGD}. Since the remaining integrals are not accounted for in the full order model, and in order to minimise intrusiveness with respect to OpenFOAM core routines, these terms are relaxed during the spatial iteration of the PGD algorithm, described in section~\ref{sc:PGD-SA}.
\end{remark}

\end{document}

config: arxiv


TITLE: Deriving Implicit Euler Method Update Rule to be used in iterations
QUESTION [0 upvotes]: Starting with the approximation
$$y'(t) = \frac{y(t) - y(t - h)}h$$
arriving at an update rule of the form:
$$y(t + h) = y(t) + hf(t + h, y(t + h))$$
Derive the implicit Euler update rule for this differential equation:
$$y'(t) = 4e^{-t} y(t) - [y(t)]^2$$  with $$y(0) = 2$$
I really have no idea what to do. Any help would be greatly appreciated.

REPLY [0 votes]: This is known as the Backward Euler method, which is a one-step method for numerically solving the equation $y' = f(t,y).$ In your case, the right hand side is:
$$f(t,y) = 4e^{-t}y - y^2$$
and the initial value is given as $y(t_0) = 2 = y_0,$ where $t_0=2.$ Given step size $h>0,$ one usually denotes $y_n$ as the $n$-th iteration or equivalently $y_n$ is the approximation to the true value $y(t_n).$
For instance: $y_1 = y_0+hf(t_1, y_1)$ and in general:
$$y_{n+1} = y_n + hf(t_{n+1}, y_{n+1}).$$
Now, the rest is just a matter of plug and chug...


TITLE: Finite Group with $n$-automorphism map
QUESTION [4 upvotes]: If $G$ is a finite group and $\phi(x) = x^{p+1}$ is an automorphism of $G$  with $order(\phi) |p$ then $G$ is a $p$-group...?
If the order of $\phi$ is $1$ then $\phi(x) = x = x^{p+1} = x^px \rightarrow x^p =e$ so the order $\forall x\in G$ is $p$ therefore $G$ is a $p$-group
I'm having trouble when considering that the order of $\phi$ is $p$
If the order of $\phi$ is $p$ then
$\phi^p(x) =x= x^{(p+1)^p}$
using the binomial theorem, I get
$\forall x \in G \, \, \, order(x) | \ \displaystyle\sum_{k=1}^p \binom{p}{k}p^k = (p+1)^p-1$
At first, I thought that the only divisors of this was powers $p$ so I got this:
Suppose for contradiction that $G$ is not a $p$-group. Let $|G| = kp^n$ where $k$ is not a multiple of $p$  ($p\nmid k$). If $k>1$ then take a $q$ in $k$'s prime factorization. So we have $q|k$ and by Cauchy's Theorem $\exists y \in G$ with $y^q = e$ i.e $order(y) = q$.
Since  $\forall x \in G \, \, \, order(x) | \ \displaystyle\sum_{k=1}^p \binom{p}{k}p^k$ we must have that $q|\displaystyle\sum_{k=1}^p \binom{p}{k}p^k$. A contradiction, therefore $k=1$ thus $G$ is a $p$-group.
But through examples, the divisors of $(p+1)^p-1$ also include other primes :'(.  How can I  show $G$ is a $p$-group?
Thanks so much :D
Update:
Taking the advice from the comments, Assume every proper subgroup of $G$ is a $p$-group, now take the Sylow $p$-subgroup of $G$, say $P$ then $P$ is normal since it is characteristic i.e $\phi(P) = P$. Since $P$ is a normal Sylow subgroup, it contains all $p$ subgroups of $G$ therefore $P$ contains all subgroups of $G$.
Can I conclude that $P=G$ since $P$ contains all subgroups of $G$?

REPLY [5 votes]: This is false. Let $G$ be cyclic of order 7, and let $p=3$. Then $\phi:G \to G : g \mapsto g^4$ is an automorphism of order 3.
Update: The question has been changed now. You can read the original version here: https://math.stackexchange.com/revisions/823130/6


TITLE: Group action on Grasmannian
QUESTION [3 upvotes]: A book i'm reading introduces the notion of a group $G$ acting on a set $X$ and then lists these two examples. I have two questions:
Example 1) The group $GL(V)$ of linear bijections from a vector space $V$ to itself
Example 2a) The group O($E$) of linear bijections that are also isometries from a Euclidean space $E$ to itself.
Example 2b) The book then defines $G_{E,p}$ to be the set of $p$-dimensional subspaces of $E$. And says that O($E$) also acts on $G_{E,p}$.
Question 1) Does the group $GL(V)$ act on the corresponding set $G_{V,p}$? Or do we somehow need $V$ to carry a Euclidean structure?
Question 2) The book then introduces the notion of a stabilizer of an element under a group action, and says that for $X = G_{E,p}$, we have that the stabilizer of an element of this set is isomorphic to $O(p)xO(n-p)$ where $O(p) = O(\mathbf{R^p})$.
Can someone clarify this?

REPLY [1 votes]: For 1) Yes. An element $\phi \in \mbox{GL}(V)$ must send a basis of $W\in G_{V,p}$ into $p$ linearly independent vectors which generate $\phi(W)$.
Now for 2) fix $W \in  G_{E,p}$. We have that $E = W \oplus W^\perp$ and any $\phi \in O(E)$ satisfying $\phi(W) = W$ must satisfy $\phi(W^\perp) = W^\perp$. Chosing an appropriate basis for $E$ we see that
$$
\phi = \begin{pmatrix}
A & 0 \\
0 & B
\end{pmatrix}
$$
where $(A,B) \in \mbox{GL}(W)\times \mbox{GL}(W^\perp)$. From $\phi$ being orthogonal, it follows that $$(A,B) \in O(W)\times O(W^\perp)\simeq O(p)\times O(n-p).$$


TITLE: Polynomials Question: Proving $a=b=c$.
QUESTION [7 upvotes]: Question:

Let $P_1(x)=ax^2-bx-c, P_2(x)=bx^2-cx-a \text{   and } P_3=cx^2-ax-b$
  , where $a,b,c$ are non zero reals. There exists a real $\alpha$ such
  that $P_1(\alpha)=P_2(\alpha)=P_3(\alpha)$. Prove that $a=b=c$.

The questions seems pretty easy for people who know some kind of calculus. Since, the question is from a contest, no calculus can be used. I have got a solution which is bit long(No calculus involved), I'm looking for a simplest way to solve this, which uses more direct things like $a-b|P(a)-P(b)$ .

REPLY [8 votes]: We have,
$a \alpha^2-b \alpha-c= b \alpha^2-c \alpha-a =c \alpha ^2-a \alpha-b=k$
Eliminating $\alpha^2$ would give me,
$(ca-b^2)\alpha-(bc-a^2)=k(b-a)$
$(ab-c^2) \alpha-(ca-b^2) =k(c-b)$
$(bc-a^2)\alpha -(ab-c^2)=k(a-c)$
Adding all of them would give me,
$(a^2+b^2+c^2-ab-ca-bc)(\alpha-1)=0 \implies \alpha =1$ or the other term is zero.
We know what happens when  $(a^2+b^2+c^2-ab-ca-bc)$ is zero. When $\alpha=1$
$a-b-c=b-c-a=c-a-b \implies a=b=c$.


TITLE: Abel Summation for d(n)/n
QUESTION [1 upvotes]: If $d(n) = \sum_{d|n}1$ is the divisor function, it is a fact that $d$ has summatory behaviour
\begin{equation}
\sum_{n \leq x} d(n) = x \log x + (2 \gamma - 1)x + O(\sqrt{x}).
\end{equation}
I'm interested in the behaviour of $\sum_{n \leq x} \frac{d(n)}{n}$ for large $x$. Now, this answer uses the Dirichlet Hyperbola Method, which I understand, to get the result
\begin{equation}
\sum_{n \leq x} \frac{d(n)}{n} = \frac{1}{2} \log^2 x + 2 \gamma \log x + \gamma^2 -2\gamma_1 + O\left(\frac{1}{\sqrt x}\right).
\end{equation}
However, if I try using Abel summation (motivated by this similar question), I get the bound
\begin{equation}
\sum_{n \leq x} \frac{d(n)}{n} = \frac{1}{x}\left(x \log x + (2 \gamma - 1) x + O(\sqrt x)\right) -1 + \int_1^x\frac{t \log t + (2 \gamma - 1) t + O(\sqrt t)}{t^2}dt
\end{equation}
\begin{equation}
=\frac{1}{2}\log^2 x + 2 \gamma\log x + 2(\gamma - 1) + O\left(\frac{1}{\sqrt x} \right).
\end{equation}
Now the numerical constants $\gamma^2-2\gamma_1$ and $2(\gamma-1)$ are not equal; more worryingly I see no way that my approach could have yielded $\gamma_1$ leading me to believe the idea was fundamentally flawed. Would appreciate if someone could aide my understanding.

REPLY [5 votes]: You probably made the following mistake:
$$\int_1^x \frac{O(\sqrt{t})}{t^2}\,dt = O\left(\frac{1}{\sqrt{x}}\right)$$
That's not right.
Observe that $$\int_1^x \frac{O(\sqrt{t})}{t^2}\,dt = O\left(\int_1^x \frac{\sqrt{t}}{t^2}\,dt\right)=O\left(2-\frac{2}{\sqrt{x}}\right) = O(1)$$
which means that using Abel summation, the best approximation you can get is $$\sum_{n\le x} \frac{d(n)}{n} = \frac{1}{2}\log^2 x + 2 \gamma\log x + O(1)$$