Time Stable Reduced Order Modeling by an Enhanced Reduced Order Basis of the Turbulent and Incompressible 3D Navier–Stokes Equations

In the following paper, we consider the problem of constructing a time stable reduced order model of the 3D turbulent and incompressible Navier–Stokes equations. The lack of stability associated with the order reduction methods of the Navier–Stokes equations is a well-known problem and, in general, it is very difficult to account for different scales of a turbulent flow in the same reduced space. To remedy this problem, we propose a new stabilization technique based on an a priori enrichment of the classical proper orthogonal decomposition (POD) modes with dissipative modes associated with the gradient of the velocity fields. The main idea is to be able to do an a priori analysis of different modes in order to arrange a POD basis in a different way, which is defined by the enforcement of the energetic dissipative modes within the first orders of the reduced order basis. This enables us to model the production and the dissipation of the turbulent kinetic energy (TKE) in a separate fashion within the high ranked new velocity modes, hence to ensure good stability of the reduced order model. We show the importance of this a priori enrichment of the reduced basis, on a typical aeronautical injector with Reynolds number of 45,000. We demonstrate the capacity of this order reduction technique to recover large scale features for very long integration times (25 ms in our case). Moreover, the reduced order modeling (ROM) exhibits periodic fluctuations with a period of 2 . 2 ms corresponding to the time scale of the precessing vortex core (PVC) associated with this test case. We will end this paper by giving some prospects on the use of this stable reduced model in order to perform time extrapolation, that could be a strategy to study the limit cycle of the PVC.


Introduction
Reduced order modeling (ROM) of the complete Navier-Stokes equations by the projection of these equations upon a reduced order space, that is generated by a minimal number of spatial modes, is still an attractive research area especially when we consider turbulent flows such as the ones encountered in aeronautical engines and that feature a large range of scales. The most important issue to be addressed when performing an order reduction of the turbulent and unsteady Navier-Stokes equations is the construction of a minimal reduced order space that could cover properly the large range of scales of a turbulent fluid flow. If we find a minimal subspace that could verify these properties, then we can do further studies concerning the efficient adaptivity of the associated reduced order equations in terms of design and optimization of the components of an aeronautical engine, such as the combustion chamber and the fuel injection system. equations, whereas many computations are performed in CFD codes using the non-descriptor form of these equations. In [12], the authors proposed to use the stability of the reduced-order Galerkin models in incompressible flows in order to study the limit cycle of the hydrodynamic vortex acting on a circular cylinder. Amsallem et al. [13] have proposed a stabilization of the projection-based linear reduced order models because of a convex optimization problem that operates directly on available reduced order operators. This method was tested for computational fluid dynamics-based model of a linearized unsteady supersonic flow, the reduction of a computational structural dynamic system, and the stabilization of the reduction of a coupled computational fluid dynamics-computational structural dynamics model of a linearized aeroelastic system in the transonic flow regime.
Besides these reduced order techniques which are intrusive for the computational fluid dynamics physics, there are some new non-intrusive reduced order ones, as they rely only on the available snapshots data, without taking into account all the equations of physics as constraints, but rather some physical properties as the turbulent kinetic energy conservation by learning the orthogonal projection coefficients of the solutions over a POD basis using a metamodel or a neural network, see the work of Wang et al. [14]. These non-intrusive techniques would take into account also the parameters calibration of the closure terms and the turbulent sub-grid scale modeling in the large eddy simulation models or the reynolds average Navier-Stokes equations ones, see the work of Lapeyre et al. [15]. There are also some research directions towards combining the machine learning non-intrusive reduced order techniques with the physics based reduced order ones in order to improve the quality of these latters, see the work of Xie et al. [16].
In this paper, we are interested in the stabilization of POD-Galerkin reduced order modeling for the turbulent and incompressible Navier-Stokes equations and we propose a new stabilization and purely physical approach for the POD-Galerkin approximation of the turbulent and incompressible Navier-Stokes equations. More precisely, a simplified POD-Galerkin projection of the complete Navier-Stokes equations is performed within an extended and minimal reduced order subspace, which reproduces accurately all the scales contained in the different terms of the Navier-Stokes equations, in order to recover a proper evolution of the fluctuating TKE. The proposed approach is based on the POD representation of the velocity gradient. A solution for the issue due to the combination of POD velocity and POD gradient modes is proposed. Moreover, we point out that the proposed stabilization is based on the a priori analysis of the different velocity modes and gradient velocity modes, in order to enforce the energetic dissipative modes within the first vectors of the new reduced order basis. This a priori enrichment by scale seperation is a key point to our proposed approach, and this will lead to the desired time stability of the reduced order model. We also point out that our proposed strategy is different from the ones that propose a scalar product change while defining the correlations matrix of the singular value decomposition (SVD) step of the POD method, in order to take into account the gradient scales within the scalar product computations, typically as in [17].
The manuscript is organized as follows: the proposed theoretical framework for the construction of the enhanced reduced order basis is detailed in Section 2. In this section, we motivate the use of an a priori enhancement of the classical POD basis by POD modes associated with the gradient velocity. All the numerical results are detailed in Section 3 for a benchmark problem of a typical aeronautical injection system at Re = 45,000 and 14 millions mesh elements. The flow solver of the High Fidelity Navier-Stokes equations is first exposed. Reduced order modeling via the enriched POD is presented and analyzed. In Section 4, we show the dynamic temporal coefficients obtained by running the enhanced fluid dynamics reduced order model for very long time integrations, even longer than the one of the high-fidelity solutions that generated the enhanced reduced order basis. In Section 5, we give some conclusions and prospects to this work. More precisely, we introduce our future work concerning the use of this stabilization technique for the extrapolation of the reduced order model in time so that, we can study efficiently the limit cycle of the PVC in an aeronautical injector [18]. It is well known that the Q-criterion of the smallest vortices is larger than the one of the large coherent structures, then the PVC is masked by small scales surrounding turbulence in the LES simulation. The enhanced ROM enables us to filter the PVC throughout the reduced order simulation even for large values of Q-criterion and for very long integration times.

POD-Galerkin Reduced Order Modeling Applied to the Unsteady and Incompressible Navier-Stokes Equations
We denote by X = [L 2 (Ω)] 3 the functional Hilbert space of the squared integrable functions over a bounded 3D−open set Ω. The corresponding inner product is the kinetic energy-based one associated with the X-functional norm. They will be denoted respectively by (., .) and . . Consider v(t) ∈ X the velocity field of an unsteady incompressible flow. Denotev(t) the filtered field obtained by a given LES model. A reduced order POD subspace is obtained by the snapshots method [19]. More precisely, if we discretize the time interval to M points, then the snapshots set is given as follows: S = {v(t i ) i = 1, ..., M}. The associated POD eigenmodes Φ n , n = 1, ..., M are solutions of the following eigenvalues problems given the temporal correlations matrix: of size M × M. We denote by A n = (A i,n ) 1≤i≤M for n = 1, ..., M, a set of orthonormal eigenvectors of the matrix C. Then, the POD-eigenmodes associated withv, are given by: where (λ n ) n=1,...,M is the decreasing sequence of the positive eigenvalues of the correlations matrix C.
To achieve the POD reduced order modeling of the filtered incompressible Navier-Stokes equations, the approximated velocity field is expressed in the reduced order POD subspace: where N << 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: where div denotes the divergence operator, p(t) is the pressure field, ρ the density, ν denotes the kinematic viscosity, v 0 is the initial condition of the velocity field and H 0 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 continuous high-fidelity incompressible Navier-Stokes equations without any turbulence model taken into account. So, our reduced order modeling formulation is the one associated with the continuous Navier-Stokes equations. However, it is clear that the POD computation is associated with High-Fidelity snapshotsv(t) which are usually obtained from LES of the Navier-Stokes equations.
In general, the first POD mode Φ 1 (x) 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. However in our case the first POD mode Φ 1 (x) is kept within the ROM. This could be very valuable when we are interested in using the reduced order modeling in order to predict the flow with respect to parametric variations, or even for new geometries [20]. This enables the ROM to consider naturally the influence of the velocity fluctuations on the velocity mean.
We point out the following two remarks concerning our formulation of the reduced order modeling: Remark 1. 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.
Remark 2. 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 inlet and outlet.
We propose to tackle these remarks on account of a physical stabilization by satisfying the kinetic energy budget.

Enrichment of the POD-Galerkin ROM with the Flow Rate Driving Forces
If we integrate by part the pressure term in the reduced order model (4), then we get the following equality: where n is the normal vector to the domain boundaries δΩ.
Using the fact that the incompressibility constraint is also verified by the velocity POD modes, the pressure term could be written: We propose to model the pressure difference between the inlet Γ in and the outlet, because of a penalization for the flow rate. This could be written mathematically as follows: where v Γ in is the inlet boundary condition of the corresponding High-Fidelity Navier-Stokes equations. The reduced order model satisfying the flow rate production forces is now written as follows: where τ is the penalization coefficient and which has been chosen equal to 10,000 in the online resolution of the reduced order modeling (8).
The flow rate penalization will enforce the following equality: where V is a steady velocity field, that should represents the mean motion. Denote by v reduced (t) = N ∑ n=2 a n (t)Φ n the fluctuating reduced order velocity, then the evolution of the turbulent kinetic energy within the ROM (8) is given by: The terms (13) in the assessment of the kinetic energy represents the production rate of the kinetic energy.

Enrichment of the POD-Galerkin ROM with the Most Dissipative Scales Based on the Velocity Gradient
To recover the dissipation rate of the fluctuating TKE in (8), we propose the following numerical algorithm, based on the enrichment of velocity-based POD modes by gradient velocity-based POD modes following a new a priori approach.
The proposed enrichment algorithm is the following: .., M and truncate at N << M these POD modes. We note that N is intentionally chosen to be less than the needed number of the POD modes to represent all the features of the coherent energetic scales of the kinetic energy.

•
Compute the fluctuating POD gradient modes Ψ n = 1  9 , i, j = 1, ..., M (W the mean velocity gradient being removed from these correlations), and (β n ) n=1,...,M is the sequence of the eigenvalues of this latter matrix.
• Compute the following velocity basis functions: Perform the Gram-Schmidt orthonormalization process for the enriched set with respect to the energy-based inner product (., .). This step is the key of the enforcement of dissipative energy modes with high singular values (β n ) n=1,...,N in early ranks of the reduced order basis, which is the opposite case when considering only the classical velocity-based POD modes (dissipative energy modes are classified respectively with very small singular values).
We will show that the new reduced order basis features new modes that represent larger ranges of spatial scales than the ones encountered in the energetic classical velocity POD modes. We point out also that the intentional a priori enforcement of these new modes within the first ranks of the reduced order basis has a major role on the quality of the resulting reduced order model. We will show that it ensures the stability of the reduced order model in an efficient fashion as a result of the availability of the driving forces and the dissipative ones within a reduced number of velocity modes.

Flow Solver
For the present simulations, the low-Mach number solver YALES2 [21] 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 [22,23]. It features a central fourth-order scheme for spatial discretization while time integration of convective terms is performed with an explicit fourth-order temporal scheme. 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 [23].

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 in order to compute precisely the aerodynamic field in the primary zone of the combustion chamber. The so-called PRECCINSTA test case [24,25] is presented in Figure 1. This lean-premixed burner has been widely studied in the combustion community to validate large-eddy simulation models [22,[26][27][28][29][30][31]. The 3D turbulent flow in the complex configuration presented in Figure 1 is considered. The kinematic viscosity ν = 10 −5 m 2 /s yields a Reynolds number 45,000 based on the inlet velocity and the length of the duct. The present High Fidelity simulation runs over 512 cores during 5 days in order to obtain a physical simulation time equal to 250 ms. A velocity-based POD-Galerkin reduced order modeling is performed, and an evaluation of its accuracy and efficiency is done before and after applying the a priori enrichment by the dissipative modes. By efficiency, we mean the online time needed to solve the mesh-independent ordinary differential Equation (8). In order to build the reduced basis, 2500 snapshots of the solution and its gradient are taken, extracted at each time step of the original HF simulation. We point out the fact that these 2500 snapshots are taken from 6644 time steps of the high fidelity simulation corresponding to the final 25 ms of its total physical time. We precise that these 25 ms represent two times the flow through time (FTT) of this test case.

POD Modes Computation for the Preccinsta
The velocity-based and gradient velocity-based POD modes were computed through a distributed snapshots POD performed in the YALES2 code. The CPU ressources needed for this computation are 768 cores (24 nodes), to guarantee a memory availability to read the 2500 time snapshots. The computation runs during 6 hours for the velocity-based POD modes and 9 hours for the gradient velocity-based POD modes. These POD modes for the velocity and gradient velocity fields are shown respectively starting from                 We can see that the velocity-based POD modes contain the high-scales of the principal coherent structures of the flow.
Interestingly, compared to the velocity-based POD modes, the velocity gradient-based POD ones feature high-scales in the dissipative regions such as in the wake of the two channels of the swirler and in the wake of the combustion chamber.
Moreover, if we compare the cumulative kinetic energies (Figures 18 and 19) associated respectively with the velocity-based POD modes and the gradient velocity-based ones, we can see that fewer than 10 velocity-based POD modes are sufficient to reproduce 90% of the high-scales TKE, however we need a larger number of velocity gradient-based POD modes in order to reproduce 90% of the small and dissipatives scales of the TKE. Then, it is clear that the dissipative scales are not considered in the velocity-based POD modes and should be added in order to preserve the energy conservation within the ROM.

The Enhanced Reduced Order Basis
We apply our a priori enforcement of the dissipative velocity modes defined previously by our new approach in the following fashion:

•
We choose N = 4 and start the enforcement by the new velocity modes from the 5th rank. This choice is made because we want to limit the number of classical global POD modes which do not exhibit at the end very large features of spatial scales, as we can see on the modes Φ 5 , Φ 6 , Φ 7 , Φ 8 and Φ 9 .

•
We choose N = 50 because, as already discussed, we need a large number of velocity gradient-based POD modes in order to reproduce 90% of the small and dissipatives scales of the TKE as shown on Figure 19. We give some further remarks in what follows: • We recall that the choice N = 4 is done intentionally in order to retrieve some dissipative modes at earlier stages than in the classical POD technique where we can see that even after 12 modes we do not have any modes of large scale's features.

•
The fact that the dissipative energy modes appear at late stages in the classical POD technique with very small singular values is the reason why we are not able to exploit their physical significance even if we increase the dimension of the classical POD reduced order model.

•
We add starting n = 5 velocity-based modes of high singular values and large features of scales.

•
This enrichment by small scale enforcement and separation is the key to multi-scale reproduction within the reduced order modeling. We precise once more that this approach is very different than the ones based on the change of the inner product that defines the matrix of the correlations between the instanteneous snapshots, typically the approach where the H 1 inner product is considered instead of the L 2 inner product. By our approach we enable scale separation, then small scale's enforcement, which is very hard to distinguish when performing a H 1 correlations matrix and then retrieving a complete POD basis: the small scales will remain dominated by the L 2 correlated large scales even if we perform this inner product change. Some authors use mathematical calibration in order to retrieve the small scales [17].
The new velocity-based modes Φ E 5 , Φ E 6 , Φ E 7 , Φ E 8 and Φ E 9 show very large features of spatial scales which was not observed within the classical global POD modes Φ 5 , Φ 6 , Φ 7 , Φ 8 and Φ 9 . The margin of variation of these large features, as we can see on Figures 24-28, ranges from 0 to 380 (see Figures 25,27 and 28). 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 Figures 14-17).
In what follows, we call "dissipative ROM" the reduced order model computed using the proposed enriched basis Φ E n n , whereas "non-dissipative ROM" refers to the ROM using the classical POD basis (Φ n ) n .

The Temporal Coefficients and Kinetic Energy of the Enriched Reduced Order Model and the Comparaison with the Classical POD-Galerkin Reduced Order Model
Figures 32-36 show the time history of the stabilized ROM amplitudes, when the stabilization algorithm is performed by enrichment of N = 4 POD velocity modes with N = 50 dissipative modes. We can see that these temporal coefficients obtained from the resolution of the reduced order model for a time interval equal to 25 ms which corresponds to the total time from which our data set was extracted in order to compute the POD modes for the velocity and the gradient velocity fields, tend to stabilize at the end of this resolution (from time step 2000). This could be explained by the order reduction using a limited number of modes, which means that the ROM needs to retrieve its equilibrium before the conservation of the kinetic energy. The ROM exhibits periodic fluctuations with a period of 2.2 ms which is the time scale of the precessing vortex core (PVC) associated with this test case.
This proves that the dissipative modes play a major role in the evolution of the Turbulent Kinetic Energy. Their introduction in the set of POD modes of the Galerkin ROM enables us to recover a better time evolution of the TKE in the system with fewer modes, see Figure 37. If we compare on this plot the kinetic energy evolutions respectively for the dissipative ROM and the non dissipative ROM, we can see that in the non dissipative case the plot of the kinetic energy is far from stabilization on the same time interval which is 25 ms.      These results are very encouraging to test the dynamic extrapolation of the stabilized reduced order model, so that we could access in real time (without any further offline operations) the evolution of the turbulent and incompressible flow outside the original snapshots data set. The first results of the dynamic extrapolation are shown in Section 4 in what follows.

3D Time Fields Obtained by the ROM and the High-Fidelity Model
In what follows we show in Figure 38 plots of the 3D reduced order velocity fields when the stabilized ROM is applied, compared to the 3D high fidelity velocity fields obtained by LES. Large scale features of the flow are clearly reproduced by the ROM even for very long time integrations.

CPU Time for Offline and Online Computation
In Table 1, we evaluate the efficiency of the stabilized reduced order modeling with respect to the High Fidelity simulation. Furthermore, we evaluate approximately the cost of the offline phase (including the snapshots POD, the Galerkin projection and the stabilization when applied) and the online ROM phase.
We precise that the speed-up is defined by the ratio of the ROM return time and the YALES2 return time. As a consequence of the proposed strategy we are able to enhance the accuracy of the reduced order modeling with a very good efficiency, regarding the online resolution. Furthermore, the offline effort associated with the additional stabilization algorithm scales with the high fidelity YALES2 return time.
It is important to note that the steps which are the most CPU consuming in the offline stage are the velocity-based POD and the gradient velocity-based one computations, followed by the stabilization by Gram-Schmidt. This took 18 h over 768 cores (24 nodes are required), because of the memory cost needed to read all of the 2500 time snapshots. This operation was not well distributed over the 768 cores due to the following issue: a temporal snapshot was not post-processed as one file per subdomain, i.e. the number of solution files per time step was less than the number of mesh partitions. This means that the running cores are working the available memory on the saved nodes in order to read lots of data per process. However, the Galerkin projection of the Navier-Stokes equations' operators took only three minutes over 768 cores. This is a consequence of the fact that we do not need to read any snapshots but, we read only the reduced number of the enhanced reduced order vectors and, we perform distributed scalar product and classical differentiation operations which scale with the mesh complexity but are very well parallelized due to distributed tasks on the mesh parts.

Temporal Extrapolation of the Dissipative ROM
Running the reduced order model for 250 ms (i.e., 10 times longer than 25 ms that is the time interval over which the POD basis has been performed), yields the dynamic weight coefficients (Figures 39-43) and the evolution of the turbulent kinetic energy represented on Figure 44. These coefficients were obtained as a consequence of the run of the stable ROM over 1 core. In this case, we can legitimately state that the speed up of the reduced order modeling is 10 8 , due to the fact that we are accessing physical solutions that were not seen by the offline phase and the learning phase of the ROM.

Conclusions and Prospects
A new methodology is proposed for the stabilisation of Galerkin reduced order models by POD for the turbulent and incompressible 3D Navier-Stokes equations. The method is based on adding the necessary physics in the new reduced order space, so that all the scales modeled in the high-fidelity Navier-Stokes equations are taken into account by the reduced order model. The only ingredient which is not represented by the retained POD modes for the reduction process in the classical methodology, is the small rank scales which are responsible for the dissipation of the turbulent kinetic energy. This ingredient is added as a result of an a priori enrichment strategy and an enforcement to the velocity-based POD modes, by a minimal number of new velocity modes which contain the low dissipative energy in the new reduced order basis. This strategy shows a very good performance when applied to an unsteady turbulent flow of Reynolds 45,000 in a typical aeronautical injection system.
The prospects of this work are the use of the proposed stable reduced model in order to perform time extrapolation, that could be a way to study the limit cycle of the Precessing Vortex Core of an aeronautical injection system.

Conflicts of Interest:
The authors declare no conflict of interest.