<figure><img src="image_1.png" /></figure>

HAL
open science

# An updated Gappy-POD to capture non-parameterized geometrical variation in fluid dynamics problems

Nissrine Akkari, Fabien Casenave, David Ryckelynck, Christian Rey

► To cite this version:

Nissrine Akkari, Fabien Casenave, David Ryckelynck, Christian Rey. An updated Gappy-POD to capture non-parameterized geometrical variation in fluid dynamics problems. Advanced Modeling and Simulation in Engineering Sciences, 2022, 9 (1), pp.3. 10.1186/s40323-022-00215-x. http://hal-03697156

HAL Id: hal-03697156

https://hal.science/hal-03697156v1

Submitted on 8 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

Open Access

<figure><img src="image_4.png" /></figure>

Check for
updates

# An updated Gappy-POD to capture non-parameterized geometrical variation in fluid dynamics problems

Nissrine Akkari¹*, Fabien Casenave¹, David Ryckelynck² and Christian Rey¹

## Abstract

In this work, we propose a new method to fill the gap within an incomplete turbulent and incompressible data field in such a way to satisfy the topological and intensity changes of the fluid flow after a non-parameterized geometrical variation in the fluid domain. This work extends the one that has been published as a conference proceeding to the 2018 AIAA Scitech Forum and Exposition (Akkari et al. in Geometrical reduced order modeling (ROM) by proper orthogonal decomposition (POD) for the incompressible navier stokes equations. In: 2018 AIAA Aerospace Sciences Meeting, AIAA SciTech Forum, (AIAA 2018-1827), 2018). A single baseline large eddy simulation (LES) is assumed to be performed prior geometrical variations. The proposed method is an enhancement of the Gappy-POD method proposed by Everson and Sirovich in 1995, in the case where the given set of empirical eigenfunctions is not sufficient and is not interpolant for the recovering of the modal coefficients for each Gappy snapshot by a least squares procedure. This happens when the available data cannot be written as an interpolation of the baseline POD modes. This is typically the case when we introduce non-parameterized geometrical modifications in the fluid domain. Here, after the baseline simulation, additional solutions of the incompressible Navier-Stokes equations are solely performed over a restricted fluid domain, that contains the geometrical modifications. These local LESs that we will call hybrid simulations are performed by using the immersed boundary technique, which uses of a fluid boundary condition and the baseline velocity field. Then, we propose to update the POD modes using a local modification of the baseline POD modes in the restricted fluid domain. Furthermore, we will propose a physical correction of the latter enhanced Gappy-POD modal coefficients thanks to a Galerkin projection of the Navier-Stokes equations upon the new modes of the available data. This enhancement procedure on the global velocity reconstruction by the physical constraint was tested on a 3D semi-industrial test case of a typical aeronautical injection system and, a 2D laminar and unsteady incompressible test case. The speed-up relative to this new technique is equal to 100, which allows us to perform an exploration of two new designs of the aeronautical injection system.

**Keywords:** Proper orthogonal decomposition (POD), Updated Gappy-POD, Gappy-POD, Non-parameterized geometrical variation, Hybrid approach, Locally available data, Galerkin projection, Dynamic extrapolation, Navier-Stokes equations, Large eddy simulation (LES), Aeronautical injection system, Design exploration in industry, Efficiency, Robustness

<figure><img src="image_5.png" /></figure>

<footnote>*Correspondence:
nissrine.akkari@safrangroup.com
¹Safran Tech, Digital Sciences and Technologies Department,
Magny-Les-Hameaux, France
Full list of author information is
available at the end of the article</footnote>
Introduction

A large number of complex simulations of 3D unsteady and incompressible turbulent flows encountered in aeronautical engines, associated with different geometrical configurations, are crucial for designing new technologies. We consider the conception of the injection system in an aeronautical engine. We need multiple 3D incompressible and unsteady simulations of the non-reacting fluid flow in the primary zone of the combustor that occurs before the ignition point. These simulations are associated with different geometries of the injection system, so that the topology of the resulting incompressible and turbulent fluid flow could influence the flame stability in the combustion zone of the gas turbine. This industrial process is challenging because of the size and the complexity of these numerical simulations. We propose to use reduced order modeling technologies to speed computational return times. Recently, we have proposed a new physical methodology to stabilize the classical POD-Galerkin Reduced Order Modeling (ROM) for the turbulent and incompressible Navier-Stokes equations, in order to cover a proper evolution of the Turbulent Kinetic Energy (TKE) spectrum and guarantee a conservation of the kinetic energy within the ROM, see [1,2]. Nevertheless, if the proposed ROM is accurate for a given geometrical configuration, its accuracy is not guaranteed for complex geometrical variations, such as non-parameterized topological ones. The geometrical inaccuracy within a baseline reduced order model associated with a reference reduced order POD basis comes mainly from the first POD mode corresponding to the mean motion. The velocity field needs to be predicted correctly in the entire domain, so that the Galerkin projection is predictable, for the new geometrical configuration. To remedy for this inaccuracy, we propose an approach based on the application of the Gappy-POD algorithm. The Gappy-POD has been first introduced in [3] as a technique for reconstruction of incomplete data field, assuming that the incomplete data vector represents a solution whose behavior can be characterized with an existing snapshots set. In [4], the POD technique was applied for inverse design purpose, in order to determine the optimal airfoil shape as an interpolation of known designs. In this paper, the POD technique was also applied in an iterative procedure in order to determine a POD associated with an incomplete pressure field on an airfoil at a given angle of attack. This iterative scheme was proposed for the first time in [3] for finding empirical eigenfunctions from the gappy data and it was shown numerically that the method yields to eigenfunctions that are close to those obtained from unmarred data. In [5], the Gappy-POD was applied in order to complete velocity data obtained using particle image velocimetry, which is often marred by missing data in various spatial locations due to inconsistent seeding and other factors. The results demonstrate that the Gappy-POD can provide an estimate that is accurate within the experimental uncertainty of the measured data. In [6], the Gappy-POD approach was performed on the inverse design of various airfoil shapes. In [7], the Gappy-POD method was applied to replace erroneous measurements in digital particle velocimetry (DPIV), where a locally adaptive criterion allows for determination of the optimum number of POD modes required for the reconstruction of each replaced measurement. In [8], the Gappy-POD was applied for aerodynamic shape optimization. In [9], the Gappy-POD is used to fuse wind-tunnel measurements and computational fluid dynamics (CFD) data to provide a consistent and more comprehensive output of greater utility. In [10], the Gappy-POD was applied in a sensor-CFD data fusion procedure for airflow field estimation.
Other than Gappy-POD methods have been proposed to deal with flows in variable geometries. Ref. [11] have used the Lagrangian sensitivity with respect to shape parameters of a baseline POD, in order to compute subsequently a reduced order model at perturbed states. However, there is no guarantee these bases will be divergent free once mapped on any other geometry. This technique was applied for the 2D incompressible Navier-Stokes equations for a flow over a square cylinder. The considered mapping is the rotation of the cylinder. Ref. [12] developed a POD method for modeling nonlinear flows with deforming meshes, thanks to dynamic functions that depend on parameters associated with flow unsteadiness. Ref. [13] developed a reduction method for immersed boundary techniques, where a new method for snapshots collection and compression describes the position of moving boundaries. Ref. [14] proposed a parametric and geometric reduced basis method for the Stokes equations using the immersed boundary equations.

In this work, we consider the incompressible and unsteady Navier-Stokes equations and we assume a non-parameterized geometric variation locally in the fluid flow domain. We want to predict correctly and efficiently the velocity field in the scope of this new geometry, from a high-fidelity computational knowledge available in a Restricted Fluid Domain (RFD). The RFD surrounds the geometrical modifications. A Gappy-POD is used as a first prediction step of the flow field in the scope of the new geometry defined by the interpolation of a pre-computed POD basis associated with one baseline high-fidelity simulation. The interpolation modal coefficients are determined via an optimization problem of the squared $L^2$-error between the available high-fidelity data around the non-parameterized geometric modification and the interpolation, evaluated only in the RFD.

We propose a second correction step which requires the local modification of the baseline POD basis vectors in order to fit the coherent structures of the local geometrical available physical data. The new basis vectors are termed updated modes. The residual of the optimization Gappy-POD problem will be equal to zero in this case. We precise also that our proposed algorithm for the correction of the predicted missing data is efficient, because only one local POD performed with the local available data followed by a Gappy-POD is sufficient to correct and refine the topology of the fluid flow in the scope of the geometrical modification. Furthermore, it is important to precise that the local data around the geometrical modification is obtained thanks to the solution of a hybrid model coupling the restricted fluid domain to the reduced projection of the equations outside this domain. The coupling is an application of a penalization approach of the velocity field by setting a reference reduced order solution associated with a reference configuration as a physical fluid immersed boundary that defines unsteady and non-homogeneous boundary conditions on the RFD. This means that at the beginning we do not have access to any time and space information of the fluid flow in the new geometrical configuration, and we will access this information successively thanks to the hybrid approach and the newly proposed Gappy reduced order method. This method is termed updated Gappy-POD (UGPOD).

The paper is organized as follows: in "Proper orthogonal decomposition (POD)" section, we briefly recall the proper orthogonal decomposition technique. In "Gappy-POD" section, we recall the Gappy-POD technique. In "Model order reduction by POD" section, we give the mathematical framework of the POD-Galerkin projection of the 3D Navier-Stokes equations.
In "Proposition of an updated Gappy-POD for CFD applications" section, we propose the UGPOD technique as described above. In "Numerical framework and experiments" section, numerical applications on a 3D semi-industrial aeronautical injection system, are presented. This enables us to perform a robust design exploration of this semi-industrial injection system with respect to non-parameterized geometrical variations in the swirler, such as the opening diameter of the primary zone of the combustor or the opening angle of this zone, with a speed up of the order of 100 with respect to a LES computation.

## Proper orthogonal decomposition (POD)

We denote by $X = [L^2(\Omega)]^3$ the functional Hilbert space of the squared integrable functions over a bounded 3D-open set $\Omega$. The corresponding inner product is the one associated with the kinetic energy of an incompressible fluid flow. They will be denoted respectively by $(\cdot, \cdot)$ and $\|\cdot\|$. Consider $U(t) \in X$ the baseline velocity field of an unsteady incompressible flow, prior any geometrical modification. A reduced order POD subspace is obtained thanks to the snapshots method [15]. More precisely, if we discretize the time interval to $M$ points, then the snapshots set is given as follows: $S = \{U(t_i) \mid i = 1, \dots, M\}$. The associated POD eigenmodes $\Phi_n, n = 1, \dots, M$, computed via the snapshots POD [15] start with the solution of the following eigenvalues problem given the temporal correlations matrix:

$$C_{ij} = (U(t_i), U(t_j)), \qquad (1)$$

of size $M \times M$. We denote by $(A_n)_{n=1,...,M} = (A_{i,n})_{1\le i \le M}$ and $(\lambda_n)_{n=1,...,M}$ for $n = 1, \dots, M$, sets of respectively orthonormal eigenvectors and eigenvalues of the matrix $C$. Then, the POD-eigenmodes associated with $\lambda_n$, are given by:

$$\Phi_n(x) = \frac{1}{\sqrt{\lambda_n}} \sum_{i=1}^{M} A_{i,n} U(t_i, x), \forall x \in \Omega, \forall n = 1, \dots, M. \qquad (2)$$

## Gappy-POD

By following [3], a mask has to be defined. Here, the support of this mask is the RFD where the POD modes of the baseline simulation are not consistent with the new geometry. In what follows, we precise some notations for the remainder of the paper:

* $\Omega$: the fluid domain.

* $U$: the nominal high-fidelity (HF) solution.

* $G^{New}$: the new configuration of the fluid domain.

* $\Omega_R$: the RFD containing the geometrical modification of the initial fluid domain.

* $U(G^{New})$: the HF velocity for the new geometry of the fluid domain.

* $U^{predicted}(G^{New})$: Gappy-POD prediction.

* $U^{ROM}(G^{New})$: New ROM prediction.

These notations are represented on the following Fig. 1, for more clarification.

Now, given a baseline POD $(\Phi_i)_{i=1,...,M}$ with a baseline snapshots set of the nominal solution $U$, then the predicted intermediate velocity field for the new geometry is determined as follows in the Gappy-POD algorithm:

$$\sum_{n=1}^{N} b_n(t, G^{New}) \Phi_n(x), \qquad (3)$$
<figure><img src="image_6.png" /><figcaption>Fig. 1 Notations</figcaption></figure>

where the interpolation POD coefficients $b_n$ are determined by the following minimization of the instantaneous error between the interpolation (3) and the HF solution, on $\Omega_R$ only:

$$ (b_1(t, G^{\text{New}}), \dots, b_N(t, G^{\text{New}})) = \operatorname*{argmin}_{\beta_1, \beta_2, \dots, \beta_N} \left\| \sum_{n=1}^{N} \beta_n(t) \Phi_n - U(t, G^{\text{New}}) \right\|_{[L^2(\Omega_R)]^3}^2. \quad (4) $$

Finally, the complete predicted velocity field $U^{\text{predicted}}$ for the new geometry is obtained by:

* $U^{\text{predicted}}(t, G^{\text{New}}, x) = U(t, G^{\text{New}}, x)$ if $x \in \Omega_R$,

* $U^{\text{predicted}}(t, G^{\text{New}}, x) = \sum_{n=1}^{N} b_n(t, G^{\text{New}}) \Phi_n(x)$ otherwise.

We precise that we supposed that we know $U(t, G^{\text{New}})$ restricted to $\Omega_R$. We will propose in “Hybrid approach: local HF solution/global reduced order solution” section a procedure to derive a fastly computed prediction for $U(t, G^{\text{New}})|_{\Omega_R}$.

## Model order reduction by POD

Let us denote by $U^{\text{ROM}}$ the reduced approximation of the field given by a LES model.

To achieve the POD reduced order modeling, the approximated velocity field is expressed in the reduced order POD subspace:

$$ U^{\text{ROM}}(t, x) = \sum_{n=1}^{N} a_n(t) \Phi_n(x), \quad \forall x \in \Omega, \quad (5) $$

where, $N \ll M$ denotes the number of retained high energetic POD modes, and $a_1(t)$, $a_2(t)$, ..., $a_N(t)$ are the temporal weights which are solutions of the following coupled dynamical system:

$$ \left\{ \begin{array}{l} \frac{da_n}{dt} + (\text{div}(U^{\text{ROM}}(t)) \otimes U^{\text{ROM}}(t)), \Phi_n \right) = \nu (\Delta U^{\text{ROM}}(t), \Phi_n) - \frac{1}{\rho} (\nabla p(t), \Phi_n) \\ (q, \text{div}(U^{\text{ROM}}(t)))_{H^0} = 0, \forall q \in H^0 \\ a_n(0) = (U(0), \Phi_n) \end{array} \right. \quad (6) $$
where *div* denotes the divergence operator, *p*(*t*) is the pressure field, *ρ* the density, *ν* denotes the kinematic viscosity, *U*(0) is the initial condition of the velocity field and *H*⁰ is the subspace of the divergence free *X*-functions.

We point out the fact that the equations upon which we perform the POD-Galerkin projection are the HF incompressible Navier-Stokes equations without any turbulence model and numerical scheme taken into account. However, the POD computation is associated with HF snapshots *U*(*t*) obtained from LES of the Navier-Stokes equations.

In general, the first POD mode which describes the mean topology of the fluid flow is not kept and a ROM of the fluid dynamics equations represents only the fluctuating part. Here, POD modes are not restricted to the fluctuation part of *U*, they also approximate the mean motion. This is necessary because we are interested in using the reduced order modeling in order to predict the flow for new geometries [16]. This enables the ROM to consider naturally the influence of the velocity fluctuations on the velocity mean.

So, we point out the following two remarks concerning our formulation of the reduced order modeling:

* The POD modes contain only the energetic scales of the flow. The dissipative scales at the Taylor macro-scale are not present in the basis.

* The flow rate in the flow domain is not guaranteed except if penalization is added in the pressure term to take into account the pressure difference between the inlet and the outlet.

We proposed in [1,2] to tackle these limits thanks to a physical stabilization by satisfying the kinetic energy budget. It is based on the enrichement of the POD-Galerkin ROM with the flow rate driving forces and with the most dissipative scales based on the velocity gradient. We refer to [2] for more details concerning this enrichment strategy of the POD-Galerkin ROM.

This stabilization step is done because we need to have a reference POD basis which is rich enough to take into account the large features of scales in the case of turbulent and incompressible fluid flows. We will denote by $\Phi^E$ the dissipative basis. This is very important from the following two points of view:

* The reference dissipative POD basis constitutes a good candidate when applying the UGPOD method in order to reconstruct new incomplete snapshots set.

* The reference reduced order modeling obtained by the POD Galerkin projection of the Navier-Stokes equations upon the reference and enriched POD basis is stable in the sense of the kinetic energy conservation physical constraint. So, it reproduces well the non-homogeneous and unsteady inlet conditions of the LES in the RFD with global penalization of the reference reduced order velocity field.

**Proposition of an updated Gappy-POD for CFD applications**

**Hybrid approach: local HF solution/global reduced order solution**

We will begin our fluid flow computation by applying a hybrid approach between $\Omega_R$ and $\Omega \setminus \Omega_R$. The fluid domain outside the RFD, $\Omega \setminus \Omega_R$, is defined by a reference fluid flow forced around the RFD. In the RFD we run the finite volume HF equations of the Navier-Stokes model with unsteady and non-homogeneous boundary conditions. We denote $U^{\text{hybrid}}$ the velocity field obtained by this approach. In this paper, we propose to
solve this coupled problem by applying the classical immersed boundary method [17], widely used in the literature. Therefore, we propose to solve the Cauchy problem defined by the Navier–Stokes equations with unsteady boundary conditions on the RFD thanks to the immersed boundary technique. By construction:

$$ \bullet U^{\text{hybrid}}(G^{\text{New}})_{\Omega_R} = U(G^{\text{New}}) $$

$$ \bullet U^{\text{hybrid}}(G^{\text{New}})_{\Omega \setminus \Omega_R} = \sum_{n=1}^{N} a_n \Phi_n^E. $$

So practically, the immersed boundary solver is run on the whole fluid domain $\Omega$ by setting a level set on $\Omega \setminus \Omega_R$ where the velocity is equal to the baseline reduced order velocity fluid flow field. Hence, indirectly we define the boundary conditions of the local HF solution by the reference reduced order velocity on the interface of the RFD. We precise also that the reference reduced order solution is saved on a coarser grid for $\Omega \setminus \Omega_R$ in order to gain efficiency when computing the pressure field during the hybrid (immersed boundary) simulation along the whole fluid domain.

## Updated Gappy-POD method (UGPOD)

We will apply the UGPOD in order to update the velocity field in the wake of the new geometry ($\Omega \setminus \Omega_R$), as follows:

1. We define updated modes $\Psi_n$ from the local modification of the reference ones in such a way to include the coherent structures of the newly computed HF solution in the RFD by the hybrid approach. These updated modes are obtained thanks to the following steps:

* Compute $U^{\text{predicted}}(t, G^{\text{New}}, x)$ with the classical Gappy-POD approach (4). We precise that the Gappy-POD modes $\Phi_n$ in (4) are here the reference POD modes $\Phi_n^E$ projected on the coarse mesh of $\Omega \setminus \Omega_R$.

* Computation of a new POD basis associated with the previous reconstructed aerodynamic fields $U^{\text{predicted}}(t, G^{\text{New}}, x)$: $\Psi_n(x) = \frac{1}{\sqrt{c_n}} \sum_{i=1}^{M} c_{i,n} U^{\text{predicted}}(t_i, G^{\text{New}}, x)$, where $(c_n)_{n=1,...,M}$ are the eigenvectors of the correlations matrix defined by: $(U^{\text{predicted}}(t_i, G^{\text{New}}), U^{\text{predicted}}(t_j, G^{\text{New}}))_{[L^2(\Omega)]^3}$, and $(d_n)_{n=1,...,M}$ is the sequence of the associated eigenvalues.

2. Galerkin projection step of the Navier–Stokes equations (still with no modelling of the turbulence) over the new POD basis $\Psi_n: \vec{c}(t, G^{\text{New}})$ are the temporal coefficients of the Galerkin projection of the Navier–Stokes equations upon $\Psi_n$.

* $U^{\text{ROM}}(t, G^{\text{New}}, x) = U^{\text{hybrid}}(t, G^{\text{New}}, x)$ if $x \in \Omega_R$.

* $U^{\text{ROM}}(t, G^{\text{New}}, x) = \sum_{n=1}^{N} c_n(t, G^{\text{New}}) \Phi_n^E(x)$ otherwise, where $\Phi_n^E$ are the reference global POD modes defined on the refined reference grid.

It is important to note that the modes $\Psi_n$ are mainly with local support $\Omega_R$ and are the locally updated modes defined in step 1 by the coherent structures obtained by data compression of the snapshots set $\{U^{\text{hybrid}}(t, G^{\text{New}}, x), x \in \Omega_R\}$. Hence, the latter temporal coefficients $c_n(t, G^{\text{New}})$ describe mainly the dynamics of the new hybrid solution defined on $\Omega_R$. Moreover, it is important to notice that the
UGPOD reconstruction is done as usual following $\Phi_n^E$ but with temporal coefficients $c_n(t, G^{\text{New}})$ that on the contrary to the Gappy-POD will describe the dynamics of the new hybrid solution defined on $\Omega_R$.

We recall the notations used so far.

* $U$: the nominal HF solution.

* $U^{\text{hybrid}}(G^{\text{New}})$: obtained by the immersed boundary technique (see "Hybrid approach: local HF solution/global reduced order solution" section).

* $U^{\text{predicted}}(G^{\text{New}})$: Gappy-POD prediction.

* $U^{\text{ROM}}(G^{\text{New}})$: New ROM prediction.

We mention that the choice of the RFD is based only on the a priori knowledge of the place where the geometrical changes are imposed, with the only condition that it is a connected domain around these changes. Furthermore, we show numerically that the length of the RFD did not affect the results obtained by the UGPOD.

## Consistence of the UGPOD

Without geometrical modification, the UGPOD is error free, under the following assumption: the reference POD modes $\Phi_n^E$ are orthogonal with respect to the Gappy inner product $(\cdot, \cdot)_{[L^2(\Omega_R)]^3}$.

**Proposition 1** *The UGPOD proposed in "Updated Gappy-POD method (UGPOD)" section is consistent.*

*Proof* Without geometrical modification $U(G^{\text{New}}) = U$. There exists $N \le M$ such that $\forall i = 1, \dots, M, \|U(t_i) - \sum_{n=1}^N (U(t_i), \Phi_n^E)_{[L^2(\Omega)]^3} \Phi_n^E\|_{[L^2(\Omega)]^3}^2 = M \times 0$. In particular, this squared difference is zero on $\Omega_R \subset \Omega: \forall i = 1, \dots, M,$

$$ \|U(t_i) - \sum_{n=1}^N (U(t_i), \Phi_n^E)_{[L^2(\Omega)]^3} \Phi_n^E\|_{[L^2(\Omega_R)]^3}^2 = 0. $$ Hence, thanks to the fact that the reference POD modes are orthogonal using the Gappy inner product $(\cdot, \cdot)_{[L^2(\Omega_R)]^3}$, the unique solution of (4) for $\Phi = \Phi^E$ is $b_n(t_i, G^{\text{New}} = G^{\text{ref}}) = (U(t_i), \Phi_n^E)_{[L^2(\Omega)]^3}$. Then
$$ U_{\Omega \setminus \Omega_R}^{\text{predicted}} = \sum_{n=1}^N (U(t_i), \Phi_n^E)_{[L^2(\Omega)]^3} (\Phi_n^E)_{\Omega \setminus \Omega_R} = U(t)_{\Omega \setminus \Omega_R}. $$ Besides, $U_{\Omega_R}^{\text{predicted}} = U_{\Omega_R}$. Hence $U^{\text{predicted}} = U$, and the POD basis on these snapshots are identical, namely $\Phi^E$ and $\Psi$ are identical. $\square$

The hypothesis made for this proposition is simply based on the fact that the Gappy-POD on the RFD $\Omega_R$ is accurate because it is performed with the dissipative POD basis $\Phi^E$ associated with the same complete data over $\Omega$ that has been restricted to $\Omega_R$.

## Numerical framework and experiments

### Flow solver

For the presented simulations, the low-Mach number solver YALES2 [18] for unstructured grids is retained. This flow solver has been specifically tailored for the direct numerical simulation and large-eddy simulation of turbulent reacting flows on large meshes counting several billion cells using massively parallel super-computers [19,20]. The Poisson equation that arises from the low-Mach formulation of the Navier-Stokes equations is solved with a highly efficient Deflated Preconditioned Conjugated Gradient method [20].
In what follows we consider a 3D semi-industrial test case of a typical aeronautical injection system and, a 2D laminar and unsteady incompressible test case.

## 2D laminar and unsteady incompressible test case

### Test case presentation

In what follows, we apply the UGPOD for a 2D unsteady, laminar and incompressible fluid flow in a rectangular channel, given a constant inlet velocity on the channel inlet, an outlet boundary condition on the channel outlet and a wall boundary condition on the upper and lower walls of the channel: we have the access to two high fidelity aerodynamic simulations corresponding to two geometric configurations of the rectangular channel, as shown on Figs. 2 and 3. The 2D laminar flow in the rectangular channel presented in Fig. 3 is considered. The kinematic viscosity is $\mu = 10^{-5} \text{ m}^2/\text{s}$. The rectangular channel is discretized into 10,000 triangles. The presented two HF simulations with 100 snapshots each, are concatenated in order to build in what follows the reference POD basis.

<figure><img src="image_7.png" /><figcaption>**Fig. 2** On the left, the level set function with an obstacle on the upper and lower walls. On the right, the level set function with a new obstacle translated with respect to the first one in the horizontal direction, and with a new length</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 3** On the left, a velocity field snapshot with the first channel configuration. On the right, a velocity field snapshot with the second channel configuration, see Fig. 2</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 4** On the left, $\Phi_1^E$ the first POD mode with the set of all snapshots with the two configurations. On the right, $\Phi_2^E$ the second POD mode with this set</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 5** On the left, the level set function with two obstacles on the lower and upper walls situated respectively at the same position of the above ones, but with a new intermediate length and two different widths. On the right, a mask vector describing the restricted fluid domain $\Omega_R$ of the new solution by the red color, and the unknown flow field data zone by the blue color</figcaption></figure>

## Construction of a reference POD basis and limits of the classical Gappy-POD

In this section, we illustrate the limits of the Gappy-POD for our particular application. Let $(\Phi_i^E)_{i=1,...,N}$ be the POD basis (Fig. 4) with the baseline instantaneous snapshots set where all the snapshots are completely known.

Let $U(t, G^{\text{New}})$ (see Figs. 5 and 6) be another solution vector with a new configuration. We have two classical ways to build the fluid velocity field associated with the new technology: (1) either we apply directly an order reduction of the Navier–Stokes equations upon the POD basis associated with the baseline snapshots. Then, it is clear that the baseline POD basis functions will fail to characterize the flow data close to the new geometry. We note from Fig. 4 that the enriched POD modes have a non-zero velocity of the fluid flow even in the obstacles associated with each ones of the baseline test cases shown on Fig. 3. Hence, these POD velocity modes could not be used in a confident way to model correctly a new geometrical fluid flow. (2) Another possibility is to apply a Gappy-POD procedure if we have the access to some high-fidelity informations with the corresponding mask $n$, as shown on Fig. 5.

The complete predicted velocity field $U^{\text{predicted}}(t, G^{\text{New}})$ by Gappy-POD for the new geometry is shown on Fig. 6, and the associated high-fidelity solution in Fig. 6.

In this case when compared to the high-fidelity model solution, the Gappy-POD under estimates the velocity field in the fluid domain. This could be explained because the baseline POD modes are not interpolant within the local zoom in red defined by the mask vector.
<figure><img src="image_11.png" /><figcaption>**Fig. 6** On the left, the final time step predicted flow field $U^{\text{predicted}}(t_{\text{final}}, G^{\text{New}})$ by Gappy-POD. On the right, $U(t_{\text{final}}, G^{\text{New}})$ the final time step 2D laminar incompressible flow in the new channel configuration</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 7** On the left, the final time step corrected flow field $U^{\text{ROM}}(t_{\text{final}}, G^{\text{New}})$ by the UGPOD. On the right, $U(t_{\text{final}}, G^{\text{New}})$ the final time step 2D laminar incompressible flow in a channel with the new configuration</figcaption></figure>

## Application of the UGPOD

Now we apply the UGPOD to the laminar 2D case where we identified some limits of the classical Gappy-POD approach to tackle geometrical variations in “Construction of a reference POD basis and limits of the classical Gappy-POD” section. The new geometrical entire flow data is now illustrated on Fig. 7.

The relative $[L^2(\Omega)]^2$ instantaneous error between the corrected flow field and the high fidelity one, $\frac{\|U(t, G^{\text{New}}) - U^{\text{ROM}}(t, G^{\text{New}})\|^2_{[L^2(\Omega)]^2}}{\|U(t, G^{\text{New}})\|^2_{[L^2(\Omega)]^2}}$, is plotted on Fig. 8.

## A posteriori verification of the ROM quality on a 2D quantity of interest

In what follows we consider a 2D quantity of interest, which is the pointwise relative error on the magnitude of the unsteady velocity fields in the wake of the two obstacles, between the Gappy-POD and the proposed UGPOD. The zone in the wake of these obstacles is chosen so that the high fidelity fluid flow in not equal to zero on each cell of the mesh and outside the RFD. On Fig. 9 we illustrate the high fidelity flow of interest in the wake of the two obstacles.

The pointwise relative errors with respect to the fluid flow in the wake of the new obstacles shown on Fig. 9 associated respectively with the application of the Gappy-POD and the UGPOD are shown on Fig. 10.
<figure><img src="image_13.png" /><figcaption>**Fig. 8** The relative $[L^2(\Omega)]^2$ instantaneous error, $\frac{\|U(t, G^{\text{New}}) - U^{\text{ROM}}(t, G^{\text{New}})\|_{[t/2, \Omega]}^2}{\|U(t, G^{\text{New}})\|_{[t/2, \Omega]}^2}$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 9** The high fidelity flow of interest in the wake of the two obstacles</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 10** On the left, the pointwise errors with the Gapp-POD. On the right, the pointwise errors with the UGPOD</figcaption></figure>

## 3D turbulent and incompressible semi-industrial test case

### Test case presentation

In what follows, we apply our new approach for a 3D unsteady, turbulent and incompressible fluid flow in a fuel injection system. The main objective is to be able to have an efficient strategy for the computation of the aerodynamic field in the primary zone of the combustion chamber. The Preccinsta test case [21,22] is presented in Fig. 11. This lean-premixed burner has been widely studied in the combustion community to validate large-eddy simulation models [19,23–28].
<figure><img src="image_16.png" /><figcaption>**Fig. 11** The 3D unsteady turbulent and incompressible flow in a fuel injection system and in the primary zone of the combustion chamber, given a constant inlet velocity, an outlet boundary condition on the channel outlet and a wall boundary condition on the upper and lower walls of the channel</figcaption></figure>

The 3D turbulent flow in the complex configuration presented in Fig. 11 is considered. The kinematic viscosity $ν = 10^{-5} \text{ m}^2/\text{s}$ yields a Reynolds number 45,000 based on the inlet velocity and the length of the duct. The presented HF simulation with 14 million tetrahedra runs over 512 cores during 5 days in order to obtain a physical simulation time equal to 250 ms. In order to build the dissipative reduced basis, 2500 snapshots of the solution are taken, extracted at each time step of the original HF simulation. We point out the fact that these 2500 snapshots are taken among 6644 time steps of the HF simulation corresponding to the final 25 ms of its total physical time.

## Construction of a reference POD basis

As explained in the previous sections, we first need to construct a reference POD basis which is rich enough to take into account the large features of scales in the case of turbulent and incompressible fluid flows. This reference POD basis is obtained in this semi-industrial case of the Preccinsta burner, in association with snapshots data generated from only one LES HF simulation of the Navier-Stokes equations with the reference geometrical configuration presented in Fig. 11. We illustrate that the UGPOD is robust with only one complete HF simulation, the baseline simulation, and one dissipative large scale reduced order basis in association. All the detailed informations concerning the dissipative reduced order basis construction can be found in [2]. However, we precise that the number of the baseline POD modes used in the UGPOD method is the same as the one that have been chosen in order to ensure an accuracy of the reference reduced order solution with no geometrical change. It is in this case equal to 54 as discussed in [2].

The velocity-based and gradient velocity-based POD modes were computed through a snapshots POD. The CPU resources needed for this computation are 768 cores, to guarantee a memory availability to read the 2500 time snapshots. The computation runs during 6 h for the velocity-based POD modes and 9 h for the gradient velocity-based POD modes. However, these operations were not well distributed over the 768 cores due to the following issue: in YALES2, the post-processing would lead to at least one file per snapshot. In this case, a temporal snapshot was not post-processed as one file per subdomain, i.e.
<figure><img src="image_17.png" /><figcaption>**Fig. 12** Velocity mode $\Phi_1^E = \Phi_1$ [2]</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 13** Velocity mode $\Phi_2^E = \Phi_2$ [2]</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 14** Velocity mode $\Phi_3^E = \Phi_3$ [2]</figcaption></figure>

the number of solution files per time step was less than the number of mesh partitions which is 128 in this case. This is due to the limited number of files that we might save on the super computer, especially when considering 2500 snapshots.

By applying the dissipative POD approach, we get a new velocity-based reduced order basis as shown from Figs. 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22 and 23. The enforcement of the small scales is done starting from the 5th mode in the reduced order basis. This choice is made in order to enforce the highly non-linear scales of the flow within the earliest ranks of the reduced order basis and so to take advantage in the stability this may offer to the reduced dynamical system by the Galerkin projection. The new velocity-based modes $\Phi_5^E, \Phi_6^E, \Phi_7^E, \Phi_8^E, \dots, \Phi_{12}^E$ show very large features of spatial scales which were not observed within the classical POD modes. Moreover, the largest scales exhibit local structures in the fluid domain which are the small vortices carrying out the dissipative energy, by analogy with the gradient velocity-based POD modes (see [2] for more details concerning the dissipative reduced basis construction).

## Introducing new geometric configurations

Two new configurations of the injection system are introduced, see Fig. 24.

We point out the fact that our geometrical modifications are non-parameterized because the level set functions are not necessarily defined by a parametric function but rather by a "if loop" that delimits the spatial location of the immersed solid boundary, please see [17].
<figure><img src="image_20.png" /><figcaption>**Fig. 15** Velocity mode $\Phi_4^E = \Phi_4$ [2]</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 16** Velocity mode $\Phi_5^E$ [2]</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 17** Velocity mode $\Phi_6^E$ [2]</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 18** Velocity mode $\Phi_7^E$ [2]</figcaption></figure>

We will apply our UGPOD in order to determine for each of these two configurations, the new fluid flow topology and the new recirculation zones.

**Hybrid approach: local HF solution/global reduced order solution**

As mentioned in "Proposition of an updated Gappy-POD for CFD applications" section, the first step of our enhancement algorithm is to compute locally the new geometrical velocity field in a restricted fluid domain that is identified by a mask, see Fig. 25.

The result of the hybrid approach is illustrated in Fig. 26, when applied to configuration 1 of Fig. 24. The physical time needed in order to obtain a solution with good statistics is far less than the time we fixed in this case which is 21 ms. This physical time is ten times greater than the flow through time associated with the zoom box in red. Therefore, the maximal CPU time needed in order to perform the LES in the RFD in red is in this
<figure><img src="image_24.png" /><figcaption>Fig. 19 Velocity mode $\Phi_8^E$ [2]</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 20 Velocity mode $\Phi_9^E$ [2]</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 21 Velocity POD mode $\Phi_{10}^E$ [2]</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 22 Velocity mode $\Phi_{11}^E$ [2]</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 23 Velocity mode $\Phi_{12}^E$ [2]</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Fig. 24** Two new configurations are considered by introducing solid obstacles in the swirler's zone. These obstacles are defined using level set functions</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. 25** Mask vector that defines the RFD in red $\Omega_R$ of interest around the geometrical variations. In blue, $\Omega \setminus \Omega_R$</figcaption></figure>

case 3.5 h on 128 cores. We have a very important reduction in the CPU time needed for this local LES with respect to the complete LES, as a consequence of the reduction of the computation domain. $\Omega \setminus \Omega_R$ is defined by the reference reduced order velocity field associated with the reference configuration, which has been saved on a coarser grid for this region. We finally point out the fact that the pressure field is computed all over the fluid domain by the HF solver, because our reference dissipative reduced order model does not contain the pressure field. The pressure computation in the hybrid simulation is efficient as we did a derefinement step by 40% for the reference dissipative velocity modes in $\Omega \setminus \Omega_R$ in blue, see Fig. 25. In other words, the dissipative velocity modes, see Figs. 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22 and 23, are saved on a coarser grid in the blue zone of Fig. 25.
<figure><img src="image_31.png" /><figcaption>**Fig. 26** An unsteady solution field of the hybrid simulation: the inlet unsteady and non-homogeneous boundary conditions on the left of the RFD are well reproduced thanks to the nominal dissipative reduced order solution: the LES in the RFD is stable</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Fig. 27** Locally updated velocity modes $\Psi_1, \Psi_2, \Psi_3$, with respect to the HF solution in the RFD associated with the first configuration. The first two modes are of global support $\Omega$ and starting from the third mode, they all have a local support $\Omega_R$</figcaption></figure>

The flow must be determined in the wake of the new swirler. In what follows, we apply our UGPOD method in order to adapt the fluid flow topology outside the red box of Fig. 25 for it to follow the geometrical variation of the swirler.

### Application of the UGPOD for the two new configurations

We now apply the UGPOD method proposed in “Proposition of an updated Gappy-POD for CFD applications” section to the two new configurations, see Figs. 27 and 28 respectively, then we get the reduced order model coefficients associated respectively with these POD modes with the predicted fields by Gappy-POD (see Figs. 29, 30, 31, 32, 33 and 34) and, the calibrated entire mean flow fields with respect to the geometrical variation are represented in Figs. 35 and 36. We can remark that by the UGPOD, the reconstructed velocity outside the RFD follows the new topology of the HF velocity in the scope of the RFD. A detailed analysis of these results is considered in the following part.

Nevertheless, it is important to note that a discontinuity on the interface between the RFD and the remainder of the fluid domain is very clear on the POD modes (see Figs. 27 and 28) with the predicted fields by the original Gappy-POD, because of the impact of the actual geometrical modification on the mean motion, hence the drawback of the Gappy-
<figure><img src="image_33.png" /><figcaption>**Fig. 28** Locally updated velocity modes $\Psi_1, \Psi_2, \Psi_3$, with respect to the HF solution in the RFD associated with the second configuration. Only the first mode is of global support $\Omega$ and starting from the second mode they have a local support $\Omega_R$</figcaption></figure>

POD approach in this case. This naturally violates the divergence free condition of these POD modes, however thanks to the fact that the fluctuating POD modes within this basis are for the most of them (of local support) null in the remainder of the fluid domain, the space integrals of the terms of the Galerkin projection upon these POD modes will be performed naturally only within the RFD where the divergence free condition on the modes is verified.

The fluctuating POD modes are in their majority of local support because of the drawback of the Gappy-POD to represent the new mean motion and so this follows to all the unsteady scales which could not be taken into account, hence the Gappy-POD gives almost steady reconstructed velocity fields. Now, the reason why we might have a second mode of global support as it can be seen on Fig. 27 was investigated in [29], it is rather a vector potential mode related to the unsteady and non-homogeneous inlet boundary conditions on the RFD. The existence of such mode with a second source of potential has been proven theoretically in [29]. The rank and the number of these velocity potential modes depend only on the potential created by the unsteady and non-homogeneous boundary conditions within the corresponding fluid configuration.

Moreover a comparison has been done between the distributions of the POD modes with the predicted velocity fields by Gappy-POD and the correct POD modes associated with the HF LES solutions after the non-parameterized geometrical modifications of configuration 2. These POD modes of the global HF LES are shown on Fig. 37. It is clear that on the RFD, the POD modes with the predicted fields by the original Gappy-POD are different from the POD modes with the HF LES solutions after the geometrical modification. A POD mode is a linear combination of the snapshots following the eigenelements of an eigenvalue problem with the temporal correlations matrix. These temporal correlations will contain information on the resulting coherent structures in the POD modes of the corresponding solutions. The correlations between the predicted fields by Gappy-POD will be very strong on the RFD where we have the HF unsteady information by the hybrid simulation, however they will be very poor between the RFD and $\Omega \setminus \Omega_R$ as the
<figure><img src="image_34.png" /><figcaption>**Fig. 29** The first temporal weight by the Galerkin projection upon $\Psi: c_1(t, G^1)$</figcaption></figure>

reconstructed velocity fields by Gappy-POD are almost steady in this case, so don't contain fluctuating scales. Hence, the eigenelements with these correlations already defined in “Updated Gappy-POD method (UGPOD)” section by, $(C_n)_{n=1,...,M}$ and $(d_n)_{n=1,...,M}$, are different from the ones with the correlations between the (global unsteady) HF LES solutions with the new configuration. Hence on the RFD, the linear combination of the HF snaphots (by the hybrid approach) following $(C_n)_{n=1,...,M}$ and $(d_n)_{n=1,...,M}$ is different from the one following the singular value decomposition of the HF LES correlations i.e. the resulting POD modes in the two cases are different on the RFD.

We precise that the computation of the updated POD basis $\Psi$ was done in a completely distributed fashion by processing all the predicted snapshots as one file per subdomain which yields 128 files. Indeed, we were able to save the predicted snapshots by Gappy-POD in one HDF5 file as the Gappy-POD is performed outside the HF solver, so we are able to control the data processing. Hence, the distributed Snapshots POD over a multiple of 128 cores will be able to read effciently a large number of snapshots (2500 in this case) per computational subdomain i.e. per CPU process.

The Galerkin projection coefficients respectively with the locally updated modes for the two new configurations are shown on Figs. 29, 30, 31, 32, 33 and 34. These coefficients describe the dynamics of the local HF LES on the RFD during at least 21 ms which is the time duration of this hybrid simulation. Notice that we performed a time extrapolation of the ROM with the updated POD modes of the second configuration, until 50 ms. This is a time extrapolation of a solution with a cyclic behavior under a constant in time loading. This means that the $\bar{c}(t, G^2)$ coefficients defined in “Updated Gappy-POD method (UGPOD)” section of the UGPOD, are obtained by resolving the reduced order equations until 50 ms, while the updated POD modes $\Psi_n$ are associated with the dynamics of the local HF LES on the RFD during 21 ms. The purpose of this extrapolation is to emphasize the stability of this reduced order model (with no modelling of the turbulence) thanks to the stable reduced order strategy developed in [2] applied to this whirlpool flow.

To summarize the steps that have been done, we refer to the flowchart of Fig. 38 that represents an illustration of the applied procedure on the geometry of configuration 1.
<figure><img src="image_35.png" /><figcaption>**Fig. 30** The second temporal weight by the Galerkin projection upon $\Psi: c_2(t, G^1)$</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>**Fig. 31** The third temporal weight by the Galerkin projection upon $\Psi: c_3(t, G^1)$</figcaption></figure>

## Detailed discussion

We present in what follows the instantaneous results obtained with the UGPOD applied to configuration 1 of the burner. We compare these results with the instantaneous ones associated with the nominal (configuration) dissipative reduced order solution (see for instance [2]). Each one of these results is compared to the corresponding time step high fidelity fluid flow velocity field: on Fig. 39, we show the result of the nominal dissipative reduced order model at two different time steps. We can see the accuracy of this latter in reproducing a large scale of spatial features with respect to integration times. On Fig. 40, we show two different time steps of the UGPOD approach applied to configuration 1. We can see that the UGPOD approach allows the reproduction of the unsteady behavior with a satisfactory accuracy. This issue was one of the problem encountered in the application of the Gappy-POD approach, as proved by the structure of the updated fluctuating POD modes $\Psi_n$, being in their majority of local support.
<figure><img src="image_37.png" /><figcaption>**Fig. 32** The first temporal weight by the Galerkin projection upon $\Psi: c_1(t, G^2)$</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>**Fig. 33** The second temporal weight by the Galerkin projection upon $\Psi: c_2(t, G^2)$</figcaption></figure>

For more results concerning the accuracy of the UGPOD applied to configuration 1, we show all the three components of the mean and the root mean square flow fields obtained by this technique. These fields are compared with the corresponding high fidelity mean and root mean square components with configuration 1 and with the nominal configuration. These results are summarized from Figs. 41, 42, 43, 44, 45 and 46. Very comparable results are obtained with a RFD of which the length is smaller than the one considered until now, please see Figs. 47, 48, 49, 50, 51 and 52.
<figure><img src="image_39.png" /><figcaption>**Fig. 34** The third temporal weight by the Galerkin projection upon $\Psi: c_3(t, G^2)$</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>**Fig. 35** The new geometrical mean velocity field for configuration 1</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Fig. 36** The new geometrical mean velocity field for configuration 2</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>**Fig. 37** Velocity modes with respect to the HF solution in the complete fluid domain associated with the second configuration</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Fig. 38** A schematic representation of the UGPOD: this flowchart shows the final time step obtained by the UGPOD for configuration 1</figcaption></figure>

## Specification for the design in combustion: a posteriori verification of the ROM quality on a 1D quantity of interest

In what follows we show the recirculation zones of each one of the two mean fields with the two new configurations, given by the new Gappy approach and the complete LES. These recirculations zones are the same in reacting and non reacting cases. In order to plot these recirculation zones, we consider three axes in the fluid domain along the y-direction for $z = 0$ and for different x-positions outside the RFD, in order to compare the result of the UGPOD strategy with the one of the LES, see Fig. 53. We add also to these validation results, a comparison with the recirculation zones obtained when the classical Gappy-POD approach is applied outside the reduced fluid domain. All these results are summarized on Figs. 54, 55, 56, 57, 58 and Table 1: it is important to point out the fact that the small sized fluctuations observed on Figs. 56 and 58 are related to the fact that the temporal coefficients with the Galerkin projection of the Navier-Stokes equations upon the POD modes with the predicted velocity fields, are obtained over a time duration that corresponds to the physical time (21 ms) of the local HF simulation which is ten times
<figure><img src="image_44.png" /><figcaption>**Fig. 39** On the left, two different time steps of the nominal dissipative reduced order. On the right, the same two time steps LES velocity fields with the nominal configuration</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>**Fig. 40** One the left, two different time steps of the UGPOD applied to configuration 1. On the right, the same two time steps LES velocity fields with configuration 1</figcaption></figure>
<figure><img src="image_46.png" /><figcaption>**Fig. 41** On the left, the first component of the mean flow field by the UGPOD applied to configuration 1. In the middle, the first component of the high fidelity mean flow field with configuration 1. On the right, the first component of the high fidelity mean flow field with the nominal configuration.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>**Fig. 42** On the left, the second component of the mean flow field by the UGPOD applied to configuration 1. In the middle, the second component of the high fidelity mean flow field with configuration 1. On the right, the second component of the high fidelity mean flow field with the nominal configuration.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>**Fig. 43** On the left, the third component of the mean flow field by the UGPOD applied to configuration 1. In the middle, the third component of the high fidelity mean flow field with configuration 1. On the right, the third component of the high fidelity mean flow field with the nominal configuration.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>**Fig. 44** On the left, the first component of the root mean square flow field by the UGPOD applied to configuration 1. In the middle, the first component of the high fidelity root mean square flow field with configuration 1. On the right, the first component of the high fidelity root mean square flow field with the nominal configuration.</figcaption></figure>

smaller than the one of the global HF simulation. These fluctuations are less important for configuration 2 (see Fig. 58) thanks to time extrapolation of the temporal coefficients with this Galerkin projection upon 50 ms, which is the natural statistical behavior in fluid dynamics.

We propose in what follows two different criteria in order to evaluate the quality of the UGPOD to recover the mean velocity field with respect to the complete LES.
<figure><img src="image_50.png" /><figcaption>**Fig. 45** On the left, the second component of the root mean square flow field by the UGPOD applied to configuration 1. In the middle, the second component of the high fidelity root mean square flow field with configuration 1. On the right, the second component of the high fidelity root mean square flow field with the nominal configuration</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>**Fig. 46** On the left, the third component of the root mean square flow field by the UGPOD applied to configuration 1. In the middle, the third component of the high fidelity root mean square flow field with configuration 1. On the right, the third component of the high fidelity root mean square flow field with the nominal configuration</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>**Fig. 47** Smaller RFD. On the left, the first component of the mean flow field by the UGPOD applied to configuration 1. On the right, the first component of the high fidelity mean flow field with configuration 1</figcaption></figure>

1. The position of the recirculation zones on each 1D y-axis with respect to each one of the two walls of the aeronautical injector. If we denote by $v(y)$ a velocity value along the corresponding y-axis, then this position is defined on each wall side by:

$$y_g = \frac{\int_{1D\ y-axis} y(v(y) - \bar{v})dy}{\int_{1D\ y-axis} (v(y) - \bar{v})dy}$$

as we can see on Fig. 59. Then, we compute respectively the two relative errors with respect to the position of the two recirculation zones within the LES.

2. The $L^2$-relative error with respect to the LES on each 1D y-axis.
<figure><img src="image_53.png" /><figcaption>**Fig. 48** Smaller RFD. On the left, the second component of the mean flow field by the UGPOD applied to configuration 1. On the right, the second component of the high fidelity mean flow field with configuration 1</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>**Fig. 49** Smaller RFD. On the left, the third component of the mean flow field by the UGPOD applied to configuration 1. On the right, the third component of the high fidelity mean flow field with configuration 1</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>**Fig. 50** Smaller RFD. On the left, the first component of the root mean square flow field by the UGPOD applied to configuration 1. On the right, the first component of the high fidelity root mean square flow field with configuration 1</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>**Fig. 51** Smaller RFD. On the left, the second component of the root mean square flow field by the UGPOD applied to configuration 1. On the right, the second component of the high fidelity root mean square flow field with configuration 1</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>**Fig. 52** Smaller RFD. On the left, the third component of the root mean square flow field by the UGPOD applied to configuration 1. On the right, the third component of the high fidelity root mean square flow field with configuration 1</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>**Fig. 53** The three y-axes for the recirculation zones: for (x, z) = (0.035 m, 0), (x, z) = (0.05 m, 0) and (x, z) = (0.055 m, 0)</figcaption></figure>

The results are shown in Tables 1 and 2.

We remark first that when no geometrical modification is introduced, see Fig. 54, the ROM velocity field is accurately the one obtained by the complete LES. We remark also that the UGPOD approach allows the reduced order solution to follow the real topology of the fluid flow outside the RFD, after non-parameterized geometrical modifications. We see that the reduced order solution respects the recirculation zones of the fluid flow, see Figs. 56 and 58, and Table 2. These recirculation zones are different from the baseline
<figure><img src="image_59.png" /><figcaption>**Fig. 54** Reference geometry: comparison of the recirculation zones obtained respectively as a consequence of the dissipative Galerkin ROM by $\Phi^E$ in orange and the HF LES in black</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>**Fig. 55** G1: comparison of the recirculation zones obtained respectively as a consequence of the classical Gappy-POD in red and the HF LES in black</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>**Fig. 56** G1: comparison of the recirculation zones obtained respectively as a consequence of the UGPOD in orange and the HF LES in black</figcaption></figure>

configuration, see Fig. 54. We see also that the classical Gappy-POD approach does not provide accurate results neither in the sense of the recirculation zones, nor in the sense of the 1D $L^2$-relative errors all along the three y-axis respectively as shown in Table 1.
<figure><img src="image_62.png" /><figcaption>**Fig. 57** G2: comparison of the recirculation zones obtained respectively as a consequence of the classical Gappy-POD in red and the HF LES in black</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>**Fig. 58** G2: comparison of the recirculation zones obtained respectively as a consequence of the UGPOD in orange and the HF LES in black</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>**Fig. 59** Illustration of the positions of the two recirculation zones on each wall side</figcaption></figure>
<table><caption>**Table 1** Evaluation of the 1D L²-relative errors along the three y-axis respectively</caption><thead><tr><th>Geometry and y-axis</th><th>Gappy-POD (%)</th><th>UGPOD (%)</th></tr></thead><tbody><tr><td>G1: x = 35 mm</td><td>32</td><td>24</td></tr><tr><td>G2: x = 35 mm</td><td>36</td><td>27</td></tr><tr><td>G1: x = 50 mm</td><td>41</td><td>34</td></tr><tr><td>G2: x = 50 mm</td><td>50</td><td>27</td></tr><tr><td>G1: x = 55 mm</td><td>31</td><td>46</td></tr><tr><td>G2: x = 55 mm</td><td>32</td><td>28</td></tr></tbody></table>

<table><caption>**Table 2** Evaluation of the position of the recirculation zones of the two new geometries and their relative errors, Err1 for our method and Err2 for the Gappy-POD, with respect to the LES ones</caption><thead><tr><th>Config.</th><th>HF</th><th>UGPOD</th><th>Gappy-POD</th><th>Err1 (%)</th><th>Err2 (%)</th></tr></thead><tbody><tr><td>G1: x = 35 mm</td><td>y<sub>g1</sub> = -0.0211</td><td>y<sub>g1</sub> = -0.019</td><td>y<sub>g1</sub> = -0.0246</td><td>6</td><td>16</td></tr><tr><td>G1: x = 35 mm</td><td>y<sub>g2</sub> = 0.022</td><td>y<sub>g2</sub> = 0.020</td><td>y<sub>g2</sub> = 0.023</td><td>6</td><td>6</td></tr><tr><td>G1: x = 50 mm</td><td>y<sub>g1</sub> = -0.0256</td><td>y<sub>g1</sub> = -0.024</td><td>y<sub>g1</sub> = -0.031</td><td>5</td><td>22</td></tr><tr><td>G1: x = 50 mm</td><td>y<sub>g2</sub> = 0.03</td><td>y<sub>g2</sub> = 0.028</td><td>y<sub>g2</sub> = 0.032</td><td>6</td><td>9</td></tr><tr><td>G1: x = 55 mm</td><td>y<sub>g1</sub> = -0.027</td><td>y<sub>g1</sub> = -0.024</td><td>y<sub>g1</sub> = -0.031</td><td>10</td><td>14</td></tr><tr><td>G1: x = 55 mm</td><td>y<sub>g2</sub> = 0.028</td><td>y<sub>g2</sub> = 0.031</td><td>y<sub>g2</sub> = 0.032</td><td>10</td><td>13</td></tr><tr><td>G2: x = 35 mm</td><td>y<sub>g1</sub> = -0.021</td><td>y<sub>g1</sub> = -0.020</td><td>y<sub>g1</sub> = -0.024</td><td>5</td><td>15</td></tr><tr><td>G2: x = 35 mm</td><td>y<sub>g2</sub> = 0.0195</td><td>y<sub>g2</sub> = 0.0197</td><td>y<sub>g2</sub> = 0.0234</td><td>0.9</td><td>19</td></tr><tr><td>G2: x = 50 mm</td><td>y<sub>g1</sub> = -0.027</td><td>y<sub>g1</sub> = -0.025</td><td>y<sub>g1</sub> = -0.031</td><td>8</td><td>14</td></tr><tr><td>G2: x = 50 mm</td><td>y<sub>g2</sub> = 0.023</td><td>y<sub>g2</sub> = 0.026</td><td>y<sub>g2</sub> = 0.032</td><td>16</td><td>41</td></tr><tr><td>G2: x = 55 mm</td><td>y<sub>g1</sub> = -0.029</td><td>y<sub>g1</sub> = -0.028</td><td>y<sub>g1</sub> = -0.031</td><td>4</td><td>6</td></tr><tr><td>G2: x = 55 mm</td><td>y<sub>g2</sub> = 0.027</td><td>y<sub>g2</sub> = 0.031</td><td>y<sub>g2</sub> = 0.032</td><td>16</td><td>20</td></tr></tbody></table>

<table><caption>**Table 3** Details of the computational costs</caption><thead><tr><th>Operation</th><th>Wall clock time</th></tr></thead><tbody><tr><td>HF YALES2 solver (512 cores)</td><td>5 days</td></tr><tr><td>HF over the RFD Ω<sub>R</sub> (128 cores)</td><td>3.5 h</td></tr><tr><td>Classical Gappy-POD on Ω<sub>R</sub> (512 cores)</td><td>3 min</td></tr><tr><td>Distributed POD Ψ with the predicted fields <sup>Upredicted</sup> (t, G<sup>New</sup>) (512 cores)</td><td>3 min</td></tr><tr><td>Galerkin projection of the Navier-Stokes equations upon Ψ (512 cores)</td><td>3 min</td></tr><tr><td>Resolution of the reduced equations (1 core)</td><td>3.7 s</td></tr><tr><td>Speed up factor</td><td>100</td></tr></tbody></table>

### CPU time reduction

In Table 3, we evaluate the efficiency of the UGPOD approach with respect to the complete HF simulation.

It is important to note that the most CPU-consuming step in the proposed approach is the hybrid computation by global penalization of the velocity field: the pressure field is still computed all over the fluid domain, so an effort has been done by coarsening the grid in the scope of a new geometry. We reached a maximum CPU time equal to 3.5 h on 128 cores.

### Conclusion and prospects

In this paper we proposed a reduced order method in order to tackle geometrical variabilities within the unsteady and incompressible Navier-Stokes equations. The proposed updated Gappy-POD is based on a hybrid approach in order to compute local LES around
a new design definition and an adaptation of the classical Gappy-POD approach in order to tackle geometrical variations of turbulent fluid flows that might introduce irregularities in the flow topology and for which the classical Gappy-POD is no longer sufficient. The main idea is to update the reference POD basis of the Gappy-POD method by the coherent structures of the local LES of the new design. The optimization problem of the Gappy-POD is replaced in our method by a Galerkin projection of the governing Navier-Stokes equations on global and local POD modes with the new design. Our approach proved good results when applied to a typical aeronautical injection system. The speed-up associated with this technique is equal to 100, by taking into account the hybrid computation of the local large eddy simulation in the Reduced Fluid Domain.

**Acknowledgements**

The authors of the paper would like to thank Vincent Moureau for his technical assistance with Yales2 and fruitful discussions. Authors knowledge the support of the CNRS GdR 3587 AMORE.

**Authors' contributions**

NA derived the methodology, the algorithm and the results, and wrote the manuscript. FC contributed to the a posteriori verification criteria of the methodology. DR and FC contributed to the consistency result and reviewed the manuscript. CR reviewed the manuscript. Every co-author participated to fruitful discussions. All authors read and approved the final manuscript.

**Funding**

The study has been funded by Safran.

**Availability of data and materials**

Similar data can be reproduced using fluid dynamics softwares.

**Declarations**

**Competing interests**

The authors declare that they have no competing interests.

**Author details**

¹Safran Tech, Digital Sciences and Technologies Department, Magny-Les-Hameaux, France, ²Centre des Matériaux, CNRS UMR 7633, Mines ParisTech PSL Research University, Corbeil-Essonnes, France.

Received: 8 December 2020 Accepted: 12 February 2022

Published online: 11 March 2022

**References**

1. Akkari N, Mercier R, Lartigue G, Moureau V. Stable POD-Galerkin reduced order models for unsteady turbulent incompressible flows. In: 55th AIAA Aerospace Sciences Meeting, AIAA Science and Technology Forum and Exposition. Grapevine, Texas, USA, 2017.

2. Akkari N, Casenave F, Moureau V. Time stable reduced order modeling by an enhanced reduced order basis of the turbulent and incompressible 3D navier stokes equations. Math Comput Appl. 2019;24(2):45. https://doi.org/10.3390/mca24020045

3. Everson R, Sirovich L. Karhunen loeve procedure for gappy data. J Opt Soc Am A. 1995;8(12):1657-64.

4. Bui-Thanh T, Damodaran M, Willcox K. Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition. AIAA J. 2004;42(8):1505-16.

5. Murray NE, Seiner JM. The effects of gappy POD on higher-order turbulence quantities. In: 46th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada. 2008.

6. Csi X, Ladeinde F. A comparison of two POD methods for airfoil design optimization. In: 35th AIAA fluid dynamics conference and exhibit, Toronto, Ontario, Canada, 2005.

7. Raben SG, Charonko JJ, Vlachos PP. Adaptive gappy POD for particle image velocimetry data reconstruction. Meas Sci Technol. 2012;23(2):025303.

8. Duan Y, Cai J, Li Y. Gappy POD-based two step optimization for airfoil design. AIAA J. 2012;50(4):968-71.

9. Mifsud M, Vendl A, Hansen L-U, Gortz S. Fusing wind-tunnel measurements and cfd data using constrained gappy proper orthogonal decomposition. Aerosp Sci Technol. 2019;86:312-26. https://doi.org/10.1016/j.ast.2018.12.036

10. Jiang C, Soh YC, Li H. Sensor and cfd data fusion for airflow field estimation. Appl Therm Eng. 2016;92:149-61. https://doi.org/10.1016/j.applthermaleng.2015.09.078

11. Hay A, Borggaard J, Akhtar I, Pelletier D. Reduced-order models for parameter dependent geometries based on shape sensitivity analysis. J Comput Phys. 2010;229(4):1327-52.

12. Freno BA, Cizmas PGA. A proper orthogonal decomposition method for nonlinear flows with deforming meshes. Int J Heat Fluid Flow. 2014;50:145-59.
13. Balajewicz M, Farhat C. Reduction of nonlinear embedded boundary models for problems with evolving interfaces. J Comput Phys. 2014;274:489–504. https://doi.org/10.1016/j.jcp.2014.06.038.

14. Karatzas EN, Stabile G, Nouveau L, Scovazzi G, Rozza G. A reduced basis approach for pdes on parametrized geometries based on the shifted boundary finite element method and application to a stokes flow. Comput Methods Appl Mech Eng. 2019;347:568–87. https://doi.org/10.1016/j.cma.2018.12.040.

15. Sirovich L. Turbulence and the dynamics of coherent structures. Part III: dynamics and scaling. Q Appl Math. 1987;45:583–90.

16. Akkari N, Mercier R, Moureau V. Geometrical reduced order modeling (ROM) by proper orthogonal decomposition (POD) for the incompressible navier stokes equations. In: 2018 AIAA Aerospace Sciences Meeting, AIAA SciTech Forum, (AIAA 2018-1827). 2018.

17. Peskin CS. The immersed boundary method. Acta Numer. 2002;11:479–517.

18. Moureau V, Domingo P, Vervisch L. Design of a massively parallel CFD code for complex geometries. C R Mécanique. 2011;339(2–3):141–8. https://doi.org/10.1016/j.crme.2010.12.001.

19. Moureau V, Domingo P, Vervisch L. From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: filtered laminar flame-PDF modeling. Combust Flame. 2011;158(7):1340–57. https://doi.org/10.1016/j.combustflame.2010.12.004.

20. Malandain M, Maheu N, Moureau V. Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines. J Comput Phys. 2013;238:32–47.

21. Meier W, Weigand P, Duan XR, Giezendanner-Thoben R. Detailed characterization of the dynamics of thermoacoustic pulsations in a lean premixed swirl flame. Combust Flame. 2007;150(1–2):2–26. https://doi.org/10.1016/j.combustflame.2007.04.002.

22. Weigand P, Duan XR, Meier W, Meier U, Aigner M, Bérat C. Experimental investigations of an oscillating lean premixed ch4-air swirl flame in a gas turbine model combustor. In: Proceedings of the European combustion meeting, 2005. p. 235.

23. Lartigue G, Meier U, Berat C. Experimental and numerical investigation of self-excited combustion oscillations in a scaled gas turbine combustor. Appl Therm Eng. 2004;24(11):1583–92.

24. Roux S, Lartigue G, Poinsot T, Meier U, Bérat C. Studies of mean and unsteady flow in a swirled combustor using experiments, acoustic analysis, and large eddy simulations. Combust Flame. 2005;141(1–2):40–54. https://doi.org/10.1016/j.combustflame.2004.12.007.

25. Moureau V, Minot P, Pitsch H, Bérat C. A ghost-fluid method for large-eddy simulations of premixed combustion in complex geometries. J Comput Phys. 2007;221(2):600–14.

26. Fiorina B, Vicquelin R, Auzillon P, Darabiha N, Gicquel O, Veynante D. A filtered tabulated chemistry model for LES of premixed combustion. Combust Flame. 2010;157:465–75.

27. Franzelli B, Riber E, Gicquel LYM, Poinsot T. Large Eddy Simulation of combustion instabilities in a lean partially premixed swirled flame. Combust Flame. 2012;159(2):621–37. https://doi.org/10.1016/j.combustflame.2011.08.004.

28. Lourier JM, Stöhr M, Noll B, Werner S, Fiolitakis A. Scale adaptive simulation of a thermoacoustic instability in a partially premixed lean swirl combustor. Combust Flame. 2017;183:343–57. https://doi.org/10.1016/j.combustflame.2017.02.024.

29. Akkari N. A velocity potential preserving reduced order approach for the incompressible and unsteady Navier–Stokes equations. 2020. https://doi.org/10.2514/6.2020-1573.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Submit your manuscript to a SpringerOpen® journal and benefit from:

▶ Convenient online submission
▶ Rigorous peer review
▶ Open access: articles freely available online
▶ High visibility within the field
▶ Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com