Calibration of projection-based reduced-order models for unsteady compressible flows
CCalibration of projection-based reduced-order modelsfor unsteady compressible flows
Victor Zucatti a , William Wolf a , Michel Bergmann b a University of Campinas, Campinas, Brazil, 13086-860 b Institut de Mathmatiques de Bordeaux, Talence, France, 33405
Abstract
An analysis of calibration for reduced-order models (ROMs) is presented in this work. TheGalerkin and least-squares Petrov-Galerkin (LSPG) methods are tested on compressible flowsinvolving a disparity of temporal scales. A novel calibration strategy is proposed for the LSPGmethod and two test cases are analyzed. The first consists of a subsonic airfoil flow whereboundary layer instabilities are responsible for trailing-edge noise generation and the secondcomprises a supersonic airfoil flow with a transient period where a detached shock wave propa-gates upstream at the same time that shock-vortex interaction occurs at the trailing edge. Resultsshow that calibration produces stable and long-time accurate Galerkin and LSPG ROMs for bothcases investigated. In order to reduce the computational costs of the LSPG models, an accel-erated greedy missing point estimation (MPE) algorithm is employed for hyper-reduction. Forthe first case investigated, LSPG solutions obtained with hyper-reduction show good comparisonwith those obtained by the full order model. However, for the supersonic case the transient fea-tures of the flow need to be properly captured by the sampled points of the accelerated greedyMPE method. Otherwise, the dynamics of the moving shock wave are not fully recovered. Theimpact of di ff erent time-marching schemes is also assessed and, di ff erently than reported in lit-erature, Galerkin models are shown to be more accurate than those computed by LSPG whenthe non-conservative form of the Navier-Stokes equations are solved. For the supersonic case,the Galerkin and LSPG models (without hyper-reduction) capture the overall dynamics of thedetached and oblique shock waves along the airfoil. However, when shock-vortex interactionoccurs at the trailing-edge, the Galerkin ROM is able to capture the high-frequency fluctuationsfrom vortex shedding while the LSPG presents a more dissipative solution, not being able torecover the flow dynamics. Keywords:
Reduced-order modeling, model calibration, Galerkin projection, least-squaresPetrov-Galerkin, proper orthogonal decomposition
1. Introduction
The higher computational power achieved in the last fifty years allowed the application oftime-accurate numerical simulations of complex engineering problems. However, despite the
Email addresses: [email protected] (Victor Zucatti), [email protected] (William Wolf), [email protected] (Michel Bergmann)
Preprint submitted to Elsevier August 17, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] A ug mprovement in computer performance, accurate numerical solutions of unsteady flows are stillcostly. In such problems, high resolution numerical schemes are typically employed to resolvethe broad range of spatial and temporal scales. On one hand, small time steps are required tocapture the higher frequencies of the flow. On the other hand, simulations need to be carried outfor long periods to obtain meaningful statistics related to the low frequencies.Direct numerical simulation (DNS), large eddy simulation (LES) and detached eddy simu-lation (DES) are the typical high-fidelity methodologies applied in the study of unsteady flows.While DNS resolves all spatial and temporal scales associated with the flow, LES is able to re-solve the larger, more energetic scales, modeling the smaller, more isotropic and universal scalesof turbulence. The computational cost is reduced in DES since it combines LES with a Reynolds-averaged Navier-Stokes (RANS) approach to solve flows over more complex configurations. Allthese methodologies are associated with high computational costs, especially in applications oflow Mach number flows involving propagation of acoustic waves [1, 2] or supersonic flows withshock-turbulence interaction [3].It is in this context that reduced-order models (ROMs) stand out when compared to con-ventional methods used for computational fluid dynamics. The application of ROMs allows theconstruction of simpler models by dimensionality reduction of the problem, what leads to lowercost simulations [4]. It should be clear that ROMs will not replace traditional CFD methods, butimprove the arsenal of tools available for solving complex engineering problems. Reduced-ordermodels must be stable and accurate for long-time integration and they should be able to recoverthe main physical features of the unsteady flows investigated. Models with quick turnaroundsolution find application in preliminary design, optimization and flow control, to name a few. Inrecent years, ROMs have attracted the attention from mathematicians, physicists and engineersinterested in the solution of complex non-linear dynamical systems such as those found in fluidflows [5, 6, 7, 8].Projection of the flow governing equations into a low-dimensional subspace is probably themost widely found class of ROMs in the literature. Among these methods, we mention Galerkin[9] and least-squares Petrov-Galerkin (LSPG) projections [10]. Succinctly, order reduction ismade possible by first extracting a low-dimensional basis from data previously collected using afull order model (FOM), for instance, DNS, LES or experiments. This step is important to feedthe ROM and obtain physical insight of the problem. The low-dimensional subspace is usuallyobtained via proper orthogonal decomposition (POD) [11, 12]. Finally, a system of non-linearordinary di ff erential equations with fewer degrees of freedom is obtained in the projection stepwhile tying the problem to physical grounds.Recently, data-driven ROMs based on regression have also been gaining attention [13, 14, 15]and unsteady flow problems have been successfully modeled with techniques based on sparseregression [16] and deep neural networks [17]. Most of these methods do not impose any con-straints from the governing equations, di ff erently from Galerkin-type methods. It is shown in[18] that projection-based ROMs overcome regression methods for unsteady flows involvingnon-periodic patterns. In such cases, the regression schemes su ff er from overfitting while theprojection-based ROMs are able to represent the physical features from chaotic flows.Stability and accuracy issues are common problems when designing reduced-order models,being well documented in literature [9, 19, 20]. Typically, POD-basis truncation is pointed outas the culprit since it acts as a filter that eliminates high-frequency modes responsible for energydissipation. However, the wavenumber interactions appear due to the non-linear convection op-erators which are responsible for the creation of high frequencies and the growth of instabilities.This phenomenon is also a problem in full order models such as LES, which require subgrid2cale modeling to represent filtered small-scale flow structures. A large number of methods havebeen proposed in the last 20 years to overcome such issues and render ROMs more robust tocomplex engineering problems.Artificial viscosity models are probably the most popular methods to overcome unstableROMs since they can be directly adapted from the CFD community. This class of closure modelsis based in adding the lacking dissipation e ff ects of truncated POD modes using corrective termsin the dynamical system [21]. The simplest model considers a change in the viscosity coe ffi cientthat impacts equally on all POD modes. This idea can be easily adapted so the viscosity canimpact di ff erently on each individual mode; for example, Smagorinsky models introduce an arti-ficial viscosity that is variable in both space and time [22]. Overall, little overhead is added to theROM although the dissipation intensity is a free parameter to be not always easily determined.Adopting di ff erent inner products may also be a solution. According to [23], the impor-tance of the smaller scales can be strengthened adopting an H Sobolev norm. Nonetheless, thisapproach also has a free parameter to be defined. Studies comparing di ff erent physics-based clo-sure models are not common, but can be found for the simulation of the one-dimensional Burgersequation [22] and the 3D flow past a cylinder at Re = ff erent way to tackle the problem would be to model ROM inconsistencies using a globalblack box approach instead of trying to account for error sources individually. The ROM shouldbe capable of recovering the temporal modes in the training window, and hopefully beyond, butthis can be a challenging task even for simple flow configurations. This can be imposed byminimizing the error between the POD and ROM temporal modes which leads to a non-linearminimization problem [25]. A cheaper linear approximation can be considered [26] and it isthe method used in the present study. This methodology was previously used in reduced-ordermodeling of unsteady flows around airfoils at low Reynolds numbers [27, 28]. This class oftechniques is normally referred to as calibration (as opposed to closure) and further technicaldetails will be discussed in section 2.4.1.The performance of di ff erent calibration methods is studied in [29, 30] for a cylinder wakeflow at Re =
200 using six POD modes. For this case, linear least squares was shown to be3000 times cheaper than non-linear optimization which was, however, slightly more accurate.In the same reference, a second problem was investigated consisting of a separated flow aroundan ONERA-D airfoil at Re = where a 60-mode POD basis was obtained from PIV snap-shots. The initially unstable ROM was stabilized using the linear least-squares methodology, butthe model results were not accurate. Moreover, the iterative non-linear optimization could notprovide a stable solution for this case.An adaptation of the energy conserving sampling and weighting (ECSW) method [31] wasrecently presented for Galerkin and LSPG ROMs. The method imposes conservation by elementweighting in order to compensate for missing energy contributions and it was applied to solve3onvection-dominated problems [8, 32]. In [32], instead of using the POD temporal modes,data calibration was accounted for by directly using a subset of the FOM snapshots also usedin the POD-basis construction. Another alternative, grounded in di ff erential geometry, is able tomodel the e ff ects of truncated modes by a minimal rotation of the projection subspace [33]. Thismethod avoids adding supplementary terms to the systems of di ff erential equations. Althoughthis approach was applied to a non-conservative compressible Navier-Stokes framework, show-ing both stable and accurate results, the optimal POD basis representation is lost. Furthermore,minimization of the error between the POD and ROM temporal modes is also enforced.Imposing the error minimization of temporal modes as a constraint may be undesirable for alarge number of applications, such as flow control and many-query problems, because in theseapplications one may desire to have the existing features of the POD spatial modes with a di ff er-ent temporal evolution. In general, ROMs are built from snapshots of a single parameter spacerealization. Therefore, physical structures of di ff erent parameters may (unsurprisingly) be ill-represented or absent. Fragility to parameter variation is well-documented [21, 34] and possiblesolutions usually involve POD-basis interpolation techniques [34]. This being said, calibrationtailored for a specific spatial and temporal evolution would most probably be inadequate whenapplied in a di ff erent context and, for this reason, we focus on calibration for accurate long-timeprediction.In this work, calibration of ROMs is assessed for projection-based methods and we presenta new calibration approach for the LSPG scheme. In order to reduce the computational costs ofthis method, hyper-reduction is applied. Calibration of Galerkin and LSPG ROMs is tested forthe compressible flow past an airfoil where boundary layer hydrodynamic instabilities lead tosecondary tones in the acoustic radiation. Then, the methods are tested for a supersonic transientflow past an airfoil where the performance of calibrated ROMs is assessed. The comparisonof calibrated Galerkin and LSPG ROMs for unsteady compressible flow problems involving adisparity of temporal scales is one of the contributions of this work together with new insightson the use of hyper-reduction for transient problems.
2. Reduced-order modeling
In the ROMs studied in this work, proper orthogonal decomposition (POD) is applied tocompute low-dimensional subspaces of specific volume, velocity and pressure fields given by ζ , u , v and p , respectively. These variables are chosen because the flows of interest are compressibleand more details are provided in Section 2.3. The unsteady flow fields can be decomposed asfollows q ( x , t ) = ¯q ( x ) + N (cid:88) i = Φ i ( x ) a i ( t ) , (1)where q = { ζ, u , v , p } (cid:62) , ¯q ( x ) is the mean flow, Φ i are the orthonormal spatial eigenfunctions, a i represent the temporal modes and {·} (cid:62) is the transpose of {·} . The parameter N is the number ofdata sets extracted from the numerical simulation and i represents the mode index.The POD consists of looking for the deterministic functions Φ i that are most similar in anaveraged sense to the realizations q ( x , t ). An alternative technique introduced in [11] is adopted4ere and the resulting constrained optimization problem reduces to the following Fredholm inte-gral eigenvalue problem (cid:90) T C i j a i ( t (cid:48) ) dt (cid:48) = λ i a i ( t ) , (2)where the temporal covariance matrix C is defined by C i j = T (cid:90) Ω q ( x , t i ) q ( x , t j ) d x ≈ T (cid:104) q ( x , t i ) , q ( x , t j ) (cid:105) Π . (3)In the previous equation, Π is a symmetric positive definite matrix defining the inner-product (cid:104) q ( x , t i ) , q ( x , t j ) (cid:105) Π = (cid:104) q i , q j (cid:105) Π = q Ti Π q j . Here, matrix Π is diagonal with non-zero elements de-fined as Π ii = A i , where A i is the area associated to the i -th vector element. This choice ofinner-product [9, 12] is equivalent to the quadrature rule (cid:82) Ω u v dx ≈ (cid:80) N g i = u i v i A i , where N g isthe total number of grid points. The covariance matrix C is symmetric positive semidefinite and,therefore, allows the use of singular value decomposition to compute the eigenvalues and eigen-vectors (modes) of the POD reconstruction. Such modes are calculated so that the reconstructionis optimal in the sense of truncated mean quadratic error. The idea of writing a temporal covari-ance matrix comes from the fact that solution cost grows rapidly for large computational grids.This is an issue especially in multidimensional problems.In Eq. 2, a i are the i − th time-dependent POD eigenfunctions, also called POD temporalmodes [12], that form an orthogonal set satisfying the condition1 T (cid:90) T a i ( t ) a j ( t ) dt = λ i δ i j . (4)The associated eigenvectors Φ i , also called empirical eigenfuctions or POD spatial modes, forma complete orthogonal set and are normalized so that they can verify (cid:104) Φ i , Φ j (cid:105) Π = δ i j . Theseeigenvectors will be used in the Galerkin and LSPG projections to reconstruct the system ofordinary di ff erential equations that, in turn, will determine the evolution of temporal modes. Thespatial basis functions Φ i can then be calculated from the realizations q i and the temporal modes a i with Φ i ( x ) = T λ i (cid:90) T q ( x , t ) a i ( t ) dt . (5)Finally, reconstruction of the fluctuation fields of specific volume, velocity and pressure canbe approximated by q (cid:48) ( x , t ) ≈ M (cid:88) i = Φ i ( x ) a i ( t ) , (6)where M is the number of modes used in the reconstruction of fluctuation fields. In practicalROM applications, one seeks M (cid:28) N . The POD method is widely used in literature and we referto the following references for further information and applications [4, 9, 11, 12, 20]. Let us consider the system of non-linear partial di ff erential equations F ( q ) defined in a con-nected open region Ω ⊂ R N whose boundary Γ is well defined F ( q ) = d q dt − G ( q ) = in Ω q ( t = t ) = q q = g on Γ . (7)5n the system above, q is a function of space and time, and the non-linear operator G ( q ) is givenby the convective and di ff usive terms appearing in the mass, momentum and energy equations,herein referred to as Navier-Stokes equations. Let Φ define an orthonormal basis obtained byPOD. The state variable q is then approximated as the linear combination of this basis vector as q ≈ ˆq = ¯q + M (cid:88) i = Φ i a i , (8)where the explicit dependencies on space and time are omitted for simplicity.In general, after discretization and approximation, F ( q ) ≈ R ( ˆq ) (cid:44) for the physical problembeing solved. A solution is sought by enforcing the residual R ( ˆq ) orthogonality as M (cid:88) i = (cid:104) Ψ i , R ( ˆq ) (cid:105) Π = , (9)where Ψ i is the test space. A projection method is generally called Galerkin (Petrov-Galerkin)when the test and solution bases are equal (di ff erent), i.e., Ψ = Φ ( Ψ (cid:44) Φ ) for the followingprojection M (cid:88) i = (cid:104) Ψ i , R ( ¯q + M (cid:88) j = Φ j a j ) (cid:105) Π = . (10)Boundary conditions must be implicitly satisfied by the POD solution basis, otherwise theproblem may lead to an ill-conditioned or ill-posed reduced-order model. Homogeneous Dirich-let or Neumann boundary conditions can be inherited by the spatial modes Φ from the snapshotcollection [4]. Galerkin projection is the most popular alternative for reduced-order modeling of time de-pendent problems. This can be attributed to its implementation simplicity and solid mathematicalfoundation. Applying the Galerkin projection method ( Ψ = Φ ) to Eq. 7 we obtain M (cid:88) i = (cid:104) Φ i , M (cid:88) j = Φ j ˙a j (cid:105) Π = M (cid:88) i = (cid:104) Φ i , G ( ¯q + M (cid:88) j = Φ j a j ) (cid:105) Π . (11)This equation can be further simplified since the functions Φ i are orthonormal. Hence, a systemof ordinary di ff erential equations arises for the temporal modes as M (cid:88) i = ˙a i = M (cid:88) i = (cid:104) Φ i , G ( ¯q + M (cid:88) j = Φ j a j ) (cid:105) Π , (12)with initial conditions obtained by projection of a single snapshot in the vector basis M (cid:88) i = a i ( t ) = M (cid:88) i = (cid:104) Φ i ( x ) , q (cid:105) Π . (13)The previous system of ordinary di ff erential equations represents the ROM associated to theFOM and can be solved using a time-marching method. The right-hand side of Eq. 12 should6ot scale with the full-order model so as to achieve reduced computational cost. Following thePOD-Galerkin approach, the ROM obtained for the non-conservative form of the compressibleNavier-Stokes equations (see details in Section 2.3) can be written as M (cid:88) i = ˙a i = M (cid:88) i = e i + M (cid:88) i , j = A ij a j + M (cid:88) i , j , k = N ijk a j a k , (14)where the ODE coe ffi cients e , A and N can be found in Appendix A. It is worth mentioningthat these coe ffi cients are time-independent and, thus, need to be calculated only once, in a pre-processing step. Test and trial bases are di ff erent when a Petrov-Galerkin approach is used ( Ψ (cid:44) Φ ). In thiswork, we employ the least-squares Petrov-Galerkin (LSPG) method that has shown promisingresults for reduced order modeling [10, 35]. According to these references, the Petrov-Galerkinapproach improves robustness compared to the Galerkin technique, despite the absence of a prioristability guarantees for general non-linear problems [6]. The least-squares variant is obtained bysolving the fully discrete residual (i.e., the residual after temporal and spatial discretization)minimization problem at each n -th time-step asminimize ˆq f n ( R ( ˆq )) . (15)The objective function f n ( R ( ˆq )) is defined in the following special form f n ( ˆq ) = (cid:107) R n ( ˆq ) (cid:107) Π = (cid:104) R n ( ˆq ) , R n ( ˆq ) (cid:105) Π , (16)or equivalently M (cid:88) i = a ni = arg min (cid:107) R ( ¯q + M (cid:88) i = Φ i a ni ) (cid:107) Π , (17)where the initial conditions are also given by Eq. 13.Optimality conditions are derived from Taylor’s theorem and can be determined by exam-ining the gradient ∇ f n ( ˆq ) and Hessian ∆ f n ( ˆq ) matrices [36]. The derivatives of f n ( ˆq ) can beexpressed in terms of the Jacobian J ( ˆq ) = ∂ R ∂ ˆq . Applying the first-order necessary condition ∇ f n ( ˆq ) = ∇ f ( ˆq ) = (cid:104) J ( ˆq ) , R ( ˆq ) (cid:105) Π = (cid:42) ∂ R ( ˆq ) ∂ ˆq , R ( ˆq ) (cid:43) Π = , (18)and applying the chain-rule to the previous equation we have (cid:42) ∂ R ( ˆq ) ∂ a Φ , R ( ˆq ) (cid:43) Π = . (19)Therefore, in the LSPG method, the discrete test basis Ψ is given by Ψ = J = ∂ R ( ˆq ) ∂ a Φ . (20)7roof of the equivalence of Eqs. 15 and 19 is available in [8]. Algorithms that follow theline-search framework are commonly used for solving problems such as Eq. 15. The main ideain the line-search approach is to chose a direction ( k ) p leading to a decrease of the objectivefunction f n when moving from the current iterate ( k ) a ni to a new one ( k + ) a ni . These algorithmshalt when the accuracy criteria has been satisfied or when further progress is impossible.The steepest descent algorithm is a first-order line-search method where the steps taken areproportional to the negative gradient ( k ) p = −∇ ( k ) f . The slow first-order convergence of thismethod is counterbalanced by requiring only gradients. A second-order alternative is Newton’smethod which requires the calculation of second derivatives. When applied in the solution of Eq.15, it results in the following iterations ∆ ( k ) f ( k ) p = −∇ ( k ) f , (21a) ( k + ) a ni = ( k ) a ni + ( k ) α ( k ) p , (21b)where k = , . . . , K with K satisfying the convergence criterion. The step length ( k ) α ∈ R ∗ + cansimply be set to unity, ( k ) α =
1, or computed using a line search in the direction ( k ) p satisfying,for example, Wolfe or Goldstein conditions [36]. The systematic choice of the step length iscrucial, since it should be cheap and significantly reduce the objective function f . Also, it shouldbe noted that problems can emerge when using second derivatives, especially when the Hessianis not a symmetric positive definite matrix.The Gauss-Newton algorithm, an alternative to Newton’s method, is capable of avoidingsome, but not all, of the issues potentially emerging when using the Hessian. In this method,second derivatives are neglected leading to an approximated Hessian using only the Jacobian ∆ f ≈ (cid:104) J , J (cid:105) Π . LSPG reduced-order models combined with a Gauss-Newton solver have beenwidely and successfully used by [6, 5, 10, 37]. Applying the modified Hessian to Eq. 21 resultsin the following iterations for the search direction ( k ) p (cid:104) ( k ) J , ( k ) J (cid:105) Π ( k ) p = −(cid:104) ( k ) J , ( k ) R (cid:105) Π . (22)However, the Gauss-Newton method does not address one of the main problems of Newton’smethod and iterations may fail when the Hessian is near, or exactly, rank-deficient. The Levenberg-Marquardt method is employed in this work and it considers ∆ f ≈ (cid:104) J , J (cid:105) Π + λ I to ensure full rank,where I is the identity matrix and λ ≥ (cid:104) J , J (cid:105) Π (cid:29) λ I the search direction is similarto the direction given by the Gauss-Newton algorithm and, when (cid:104) J , J (cid:105) Π (cid:28) λ I , the method issimilar to the steepest descent.A detailed theoretical and computational comparison of Galerkin and LSPG projection meth-ods is provided by [6]. A number of interesting outcomes from this previous reference areoutlined in the following. The LSPG solution may converge to that obtained from Galerkinprojection in some situations, for example, as the time step converges to zero. Also, the regres-sion problem from Eq. 15 can be linear or non-linear depending on the time integration schemeemployed (explicit or implicit) and the set of equations being solved. Non-linear least-squaresincreases the ROM cost considerably and may become an issue. Finally, finding that the errordoes not decrease as the time step approaches zero, but is optimal for an intermediate value, isat the same time a surprising and inconvenient finding. Despite these issues, the LSPG methodhas demonstrated in many situations to be a better alternative to the Galerkin method for generalnon-linear dynamical systems in spite of the higher level of complexity. As also shown in [6], the8ethod is more robust (without hyper-reduction) than Galerkin projection and it has desirablestability properties if implemented together with implicit time-marching schemes. Projection-based reduced-order modeling may fail to produce significant computational timegains. Problems containing strong non-linearities (i.e. non-polynomial) or non-a ffi ne parameterdependence impede pre-computation of Galerkin coe ffi cients. Consequently, full order scal-ing inner-products must be systematically and consistently calculated during time integration.Similarly, the reduced POD space is not su ffi cient to provide low-cost to LSPG because theresidual minimization problem from Eq. 15 is always grid-size dependent. Hyper-reductiontechniques are capable of providing the supplementary approximation needed to obtain reason-able computational savings in the so-called “reduce-then-project” approach. Succinctly, spatialmodes are reduced using a non-random sampling algorithm before projection. In some cases,pre-computation of strong non-linearities can be made possible by the use of lifting transforma-tions [38]. On one hand, the invasiveness of the additional layer of approximation introduced byhyper-reduction is avoided. On the other hand, a non-uniquely defined lifted system has to bederived and supplementary variables are introduced.The gappy POD method [39], which was originally used in facial image reconstruction withincomplete data sets, provides the framework to the reduce-then-project procedure also usedin other studies [40, 10]. Within this approach, given a subset of indices J = { j , . . . , j s } ⊂{ , . . . , N g } , a solution vector q is approximated by shrinking the spatial modes using a maskprojection matrix P = [ e j , . . . , e j s ] T ∈ R N g × s to construct the estimated solution ˜ q ≈ P T Φ a .Here, s is the number of indices retained from the original vector of size N g and e j k denotes thevector with a 1 in the j k -th coordinate and 0’s elsewhere. Temporal modes a are then calculatedby minimizing the error (cid:15) between the gapless solution q and ˜ q . Hence, the error is defined as (cid:15) = (cid:107) P T q − P T Φ a (cid:107) Π , (23)which is equivalent to the minimization problem (cid:107) I − (cid:104) ˜ Φ i , ˜ Φ j (cid:105) Π (cid:107) , (24)where ˜ Φ = P T Φ . Random sampling is used when no dynamical insight of the problem beingsolved is previously available. This approach usually leads to considerably bigger mask matricesand randomness makes reproduction of results impossible. Fortunately, this is not the case whenusing a POD basis. Hyper-reduction has been a very active research topic in the past two decadesand, as a result, several methods have been developed and tested. The optimal solution of Eq.24 for a set of given size is an intractable combinatorial optimization problem even for relativelysmall problems. However, the missing point estimation (MPE) method [41] eliminates points byevaluating the condition number c (cid:16) (cid:104) ˜ Φ i , ˜ Φ j (cid:105) (cid:17) ≡ λ max (cid:16) (cid:104) ˜ Φ i , ˜ Φ j (cid:105) (cid:17) λ min (cid:16) (cid:104) ˜ Φ i , ˜ Φ j (cid:105) (cid:17) (25)of the approximated identity matrix (cid:104) ˜ Φ i , ˜ Φ j (cid:105) ≈ I . This procedure is performed looping over allcomponents until a user specified near-optimal set of points is left. Computational cost of theMPE method may quickly become prohibitive, although the method is cheaper than the optimalsolution. 9he discrete empirical interpolation method (DEIM) [42] is a popular substitute to MPE [43,44, 45, 46]. This method was derived from the empirical interpolation method (EIM) [47] andit first appeared as a hyper-reduction approach for equations containing strong non-linear terms.However, even though it is simple and cost e ff ective, the DEIM has an important shortcoming:the number of interpolation points has to be at most the same as the number of POD modes usedin the reconstruction. This is an obvious limiting element, especially for non-trivial evolutionequations or problems using few POD modes. Hence, one may require alternative samplingmethods capable of circumventing this limitation in order to obtain stable and accurate ROMs.The accelerated greedy MPE procedure [48] introduces a massive improvement to the classicMPE and, therefore, considerable time complexity reductions are possible while still using theMPE principle. Briefly, the greedy point selection is equivalent to picking the spatial modeindex responsible for the largest growth in the condition number of the modified eigenvalueproblem. The costly modified symmetric eigenvalue problem is not actually solved, but theindex selection occurs by analyzing properties of the candidate indices. In the present work, weadopt the accelerated greedy MPE technique for hyper-reduction. Generally, the conservative form of the compressible Navier-Stokes equations is preferred inCFD applications, especially in applications involving discontinuities. However, the non-linearrational functions of the variables impede speed-ups when applying a Galerkin projection frame-work because pre-computation of the ROM coe ffi cients is impossible in the conservative form.This issue can be avoided either by applying hyper-reduction methods to cut down complexityor choosing a non-conservative formulation as shown by [23]. The conservative form of thecompressible Navier-Stokes equations was recently applied with hyper-reduction in the contextof projection-based reduced-order models [8]. Hyper-reduction is unattractive and to be avoidedwhen possible since conservation properties are lost. Moreover, the second degree polynomialnon-linearity of the non-conservative equations allows for pre-computation and cost reductionwithout any further approximations, thus being the approach used in this study.In this work, we solve the 2 D non-conservative compressible Navier-Stokes equations pre-sented by [23]. They can be written as ζ t = ζ ( u x + v y ) − u ζ x − v ζ y , (26a) u t = − uu x − vu y − ζ p x + MaRe ζ (cid:20)(cid:18) u x − v y (cid:19) x + ( v x + u y ) y (cid:21) , (26b) v t = − uv x − vv y − ζ p y + MaRe ζ (cid:20)(cid:18) v y − u x (cid:19) y + ( v x + u y ) x (cid:21) , and (26c) p t = − up x − vp y − γ p ( u x + v y ) + γ MaRePr [( p ζ ) xx + ( p ζ ) yy ] + ( γ − MaRe (cid:20) u x (cid:18) u x − v y (cid:19) + v y (cid:18) v y − u x (cid:19) + ( v x + u y ) (cid:21) , (26d)where ζ = /ρ is the specific volume, u and v are the x and y velocity components, respectively,and p is the pressure. In the set of Eqs. 26, Pr , Re and Ma denote the reference Prandtl, Reynoldsand Mach numbers, respectively. Subscripts denote partial derivatives and γ is the specific heatratio. 10 .4. Calibration of projection-based reduced-order models Calibration methods are mostly developed and applied to Galerkin-type ROMs. They provideadditional terms to the generally non-linear system of ordinary di ff erential equations arising fromprojecting the POD modes in the governing equations ˙a = f ( a ) , (27)where f ( a ) is the right-hand side of Eq. 14. Calibration could be, in principle, applied to anyROM originating from an initial value problem. For example, calibration of data-driven ROMs[16, 17] is a work in progress. Usually, only linear terms are adopted in the calibration process,but non-linear terms can also be considered and have been used to some extent [49]. Here,constant e c and linear A c calibration terms will be added as ˙a = f ( a ) + e c + A c a , (28)and details of the numerical procedure are provided in the following sections. In a linear least-squares calibration, the goal may be to minimize the error E between tem-poral modes obtained by solving the system of non-linear equations of the reduced-order model a ROM and the original POD temporal modes a POD in a user specified training window 0 ≤ t ≤ T as E = M (cid:88) i = (cid:90) T (cid:16) a PODi ( t ) − a ROMi ( t ) (cid:17) dt . (29)Similarly, temporal modes should also satisfy ˙a = f ( a ) and, hence, a second error norm E canbe defined by the ODE system and this is the preferred method used in this work E = M (cid:88) i = (cid:90) T (cid:16) ˙a PODi ( t ) − f i ( a ROM ) (cid:17) dt . (30)Linearization is obtained by enforcing f ( a ROM ) = f ( a POD ) in the error norm E [26], where ˙a PODi ( t ) is computed using a high-order finite di ff erence scheme. Suppression of the non-linearconstraint to obtain an a ffi ne operator may seem extreme at first, but this can be understood asa measure of how the non-calibrated ROM alters the solution in each time step. If successful,calibration should systematically compensate for any deviations from the POD temporal modes E = M (cid:88) i = (cid:90) T (cid:16) ˙a PODi ( t ) − f i ( a POD ) (cid:17) dt . (31)Here, E will be used to indicate the non-calibrated a ffi ne approximation error while E c willindicate the same error with constant e c and linear A c calibration terms included E c ( e c , A c ) = M (cid:88) i = (cid:90) T (cid:16) ˙a PODi ( t ) − f i ( e c , A c , a POD ) (cid:17) dt . (32)Calibrating the ODE system by directly minimizing the error E c inside the training windowcould fail to generalize the knowledge to longer time periods. Overfitting can be overcome by11ontrolling the relation between the calibrated and original ROM terms. In this case, a functional J ( e c , A c ,θ ) can be defined as J ( e c , A c ,θ ) = θ E c ( e c , A c ) E + (1 − θ ) (cid:107) e c (cid:107) + (cid:107) A c (cid:107) (cid:107) e (cid:107) + (cid:107) A (cid:107) , (33)where the first term in the right hand side corresponds to the normalized prediction error and thesecond one provides the relative weight of the calibration coe ffi cients compared to the coe ffi cientsof the original ROM. The parameter θ ∈ [0 , θ close to 0 add more importance tothe original ROM while θ close to 1 add a higher weight to the prediction quality of the calibratedmodel along the training window. Other functional choices, favoring for example constant terms,could be explored. In this work, we use the following functional form˜ J ( e c , A c , θ ) = E c ( e c , A c ) + ˜ θ ( (cid:107) e c (cid:107) + (cid:107) A c (cid:107) ) (34)with ˜ θ = − θθ E (cid:107) e (cid:107) + (cid:107) A (cid:107) . (35)Minimizing the functional ˜ J gives rise to a linear system. For clarity, the di ff erent calibrationcoe ffi cients are grouped such that K c = [ e c A c ] and the enriched temporal mode vector a ∗ = [1 a . . . a M ] = [1 a ], so that the functional to be minimized is written as˜ J ( K c , θ ) = M (cid:88) i = (cid:90) T ˙ a i ( t ) − f i ( a POD ) − M + (cid:88) j = K ci j a ∗ j ( t ) dt + ˜ θ (cid:107) K c (cid:107) , (36)where f i is the general non-linear system of ordinary di ff erential equations defining the reduced-order model being calibrated.For a given parameter θ , the optimality condition ∂ ˜ J /∂ K ci j = M linearsystems of size M + i = , . . . , M D T K ci = b i , (37)where K ci is the i th row of K c , D i j = (cid:90) T a ∗ i ( t ) a ∗ j ( t ) dt + ˜ θδ i j (38)and b ji = (cid:90) T (cid:16) ˙ a j ( t ) − f j ( t ) (cid:17) a ∗ i ( t ) dt . (39)Matrix D is calculated only once while vector b i must be evaluated for each mode during thecalibration phase. Further details and explanations can be found in [50]. The LSPG method should provide better stability properties compared to Galerkin projec-tion. However, the method still lacks a priori stability and accuracy guarantees, especially whencombined with hyper-reduction. Thus, the technique could possibly benefit from calibration butit has to be tailored to fit the approach developed in Section 2.4.1 because it solves the fully12iscrete residual minimization problem. In other words, the problem must be applicable in theerror norm given by Eq. 32.This inconvenience can be worked around by extracting a function f ( a ) equivalent to theright-hand side of Eq. 27 from the LSPG formulation. Instead of solving the fully discreteresidual minimization problem, the trick is to solve the spatial-discrete residual minimizationversion of the problem for the time derivative of the temporal mode ˙a = f ( a ) at each n -th time-step using the POD trial basis projection with hyper-reduction f ( a nPOD ) = arg min (cid:107) R ( a nPOD ) (cid:107) Π . (40)Although hyper-reduction is non-compulsory in Eq. 40 when evaluating f ( a nPOD ), it is stronglyencouraged for a cheaper calibration procedure. After all values of f ( a nPOD ) are determined, cali-bration terms can be easily calculated according to Section 2.4.1. Finally, constant e c and linear A c calibration terms are integrated systematically afterwards to the model at each time-step in apredictor-corrector fashion as ˜a = arg min (cid:107) R ( a n ) (cid:107) Π , (41a) a n + = ˜a + (cid:90) ∆ t ( e c + A c a n ) dt . (41b)Calibration could be particularly beneficial to systematically accounting for the additionalerrors introduced by hyper-reduction. Moreover, smaller mask matrices (i.e., a more aggressivehyper-reduction) could be envisioned when combined with calibration, what would reduce boththe computational costs and model errors.
3. Results and Discussion
In this section, we analyze the performance of Galerkin and LSPG methods with and withoutcalibration for generating reduced-order models of compressible flows. Both ROM approacheswere previously tested in [18] for the solution of incompressible flows involving convective heattransfer. In this previous reference, the LSPG method was able to present accurate solutionswith an aggressive hyper-reduction that used only 0 .
05% of the grid nodes. In this work, ROMcalibration is tested on compressible flows involving wave propagation. Firstly, we study theflow past a NACA0012 airfoil at freestream Mach number M ∞ = . Re c = , M ∞ = . Re c = , ff erence scheme [51] for spatial discretiza-tion. The method employs a staggered grid formulation and is able to capture shock wavesfor low Mach number supersonic flows ( M ∞ < .
3) without the explicit addition of a shockcapturing scheme. A hybrid implicit-explicit framework is employed for time marching of the13 a) -1 Frequency -8 -7 -6 -5 -4 (b)Figure 1: Sketch of 2D compressible flow at moderate Reynolds number over NACA0012 airfoil (left). Flow instabilitiesdeveloping over suction side boundary layer are advected along the trailing edge leading to acoustic scattering. Theseflow instabilities are modulated by a low-frequency motion of the separation bubble, what induces the appearance ofequidistant secondary tones in the noise radiation as shown in [55]. Pressure spectrum computed one chord above thetrailing edge showing a main tone plus secondary tones (right). equations through a combination of a third-order Runge Kutta scheme with a modified Beam-Warming implementation [52]. All flow simulations are conducted with O-type grids with N x × N y grid points in the streamwise and wall-normal directions, respectively. The present grids have (cid:16) N x × N y (cid:17) = (800 × × In the present analysis, a NACA0012 airfoil with a rounded trailing edge is immersed in a M ∞ = . Re c = , ff ect the acoustic scatteringmechanism, leading to multiple equidistant secondary tones in the acoustic spectrum as shownFig. 1(b). Therefore, despite the simple geometrical configuration, this airfoil flow at moderateReynolds number is a suitable candidate to evaluate the performance of ROMs since it o ff ers richdynamics. More details about this flow can be found in [55].Results obtained from the full order model are recorded for 10 ,
000 snapshots with a constantnon-dimensional time step ∆ t snapshot = × − . Half of the snapshots are used to construct a14
50 100
POD index mode R I C (a) POD index mode -4 -3 -2 -1 S i ngu l a r v a l ue (b)Figure 2: Relative information content (left). Spetrum of singular values (right) reduced-order basis by the snapshot-POD method introduced in section 2.1. This is equivalent to15 non-dimensional time units. The length of the training window and the number of samplingsnapshots need to be chosen wisely to avoid an ill-resolved POD basis. In this case, the low-frequency motion of the separation bubble requires a wide training window. On the other hand,high sampling frequency rates are essential for resolving the finer flow scales appearing along theboundary layer. This combination of features is crucial for stability and accuracy of the presentnon-linear ROM.The optimality property of the POD method is expected to produce a basis where only asmall number of modes should be necessary to reconstruct the input data and, thus, the benefit ofincluding additional modes is expected to rapidly decay. Usually, the number of POD modes usedin the ROM is chosen according to the relative information content RIC( M ) = (cid:80) Mi = λ i / (cid:80) Ni = λ i and should satisfy a predefined threshold. The evolution of the RIC for this case is presented inFig. 2(a) for a basis composed of the first 100 modes (out of 5,000). This basis represents 99 . ffi cient to lead to an accurate flow representation.Additionally, the corresponding spectrum of singular values is presented in Fig. 2(b), where it ispossible to see the fast magnitude decay of the first modes. As can be also seen in this figure, themagnitudes are similar for mode pairs, what indicates that such modes contain similar frequencyinformation, di ff ering only with respect to phase.It is expected that the first mode pairs represent the most energetic coherent structures re-sponsible for the flow dynamics. In order to have a better understanding of the present flow, thePOD spatial eigenfunctions are shown in Fig. 3 for modes 1, 5, 9, 13 and 17. These modesare presented for u and v velocity components, and pressure fluctuations ( ζ (cid:48) has the same spatialdistribution as p (cid:48) ). They are chosen according to their dynamical content, depicted in Fig. 4. Forall modes, it is possible to notice that the flow structures appear along the suction side boundarylayer and wake regions. Analyzing both figures together, it is clear that mode 1 contains the mostenergetic large-scale structures on the boundary layer, which are associated with the main tone ofthe spectrum shown in Fig. 1(b). Mode 5 is associated with finer scales and displays a modulated15 igure 3: Contours of POD eigenfunctions for modes 1, 5, 9, 13 and 17 (top to bottom) for u -velocity (left), v -velocity(center) and pressure (right). higher frequency content. Modes 9, 13 and 17 share similarities in terms of flow structures as canbe observed from Fig 3. The temporal dynamics of the former two modes are similar and includemultiple frequencies strongly modulated. However, the main frequency of mode 9 is lower thanthat of mode 13. On the other hand, the temporal dynamics of mode 17 appears to have a moreclear pattern composed of sinusoids with weaker modulation.For this case, the ROMs are constructed using the first 100 POD modes. For the Galerkinscheme, the constant, linear and non-linear coe ffi cients are computed and stored in a pre-processingstage following the equations shown in Appendix A. The spatial derivatives of the POD modesare computed using a 10 th -order accuracy compact finite di ff erence scheme for both Galerkinand LSPG techniques. Zucatti et al. [18] show that the error evolution of projection-basedROMs is sensitive to the computation of spatial derivatives, being considerably reduced whenhigh-resolution schemes are applied.The present Galerkin model is computed using a maximum calibration parameter ( θ = ff erent values. For this particular flow, the solution error is notvery sensitive with respect to the calibration parameter, as can be observed in Fig. 5. This figureshows that for θ = (cid:107) · (cid:107) F of the constant and linear calibratedterms is an order of magnitude smaller than that computed for the original Galerkin operators.16 t -0.02-0.0100.010.02 a t -0.010.01 a t -0.010.01 a t -0.010.01 a t -0.010.01 a Figure 4: Temporal dynamics of POD modes 1, 5, 9, 13 and 17. -4 -5 -4 -3 -2 -1 Figure 5: Ratio of Frobenius norms computed for calibrated and Galerkin coe ffi cients as a function of approximationerror for di ff erent values of θ (see Eq. 34 for details). Thus, relatively low-intrusive terms are obtained in the calibration procedure even though theprediction quality along the training window is favored during calibration. However, this is notalways the case and, for other flows, it is important to choose θ aiming to improve the modelquality and, at the same time, avoid overfitting. The calibrated LSPG models do not dependexplicitly on θ since they are computed using Eq. 41.As previously discussed, the LSPG method should provide better stability properties than theGalerkin projection. However, as stated by [10], its computational cost can scale with that of theFOM if all degrees of freedom are used on the model reconstruction. In order to reduce the modelassembly cost, we use the hyper-reduction technique previously described. Before a particularnumber of grid points is employed in the hyper-reduction, it is important to assess the conditionnumber behavior of the system. This is a good way to have an a priori estimate of the number ofpoints to be used in the ROM. On one hand, we desire to use as few as possible grid points in themodel reconstruction to reduce its computational cost. On the other hand, the condition numbershould be kept as small as possible to reduce the errors associated with the procedure and to beable to represent the relevant flow dynamics. Figure 6 shows the condition number as a functionof the number of grid points used in the hyper-reduction.One can see that the condition number does not fall monotonically, what is expected sincethe present hyper-reduction algorithm is locally optimal but not globally optimal, as previouslydiscussed. Eventually, as the number of grid points is increased, the condition number willbecome unity once the approximation mode matrix recovers the original full POD spatial modematrix. As can be observed from the figure, there is a first considerable drop in the conditionnumber when approximately 50 grid points are used in the mask matrix. However, the LSPGROMs built for this range of condition number were either unstable or inaccurate. Increasingthe number of grid points resulted in a second drop in the condition number which providedstable models. For the present problem, hyper-reduction is then applied using 29 ,
000 grid points(6 .
04% of the total) which resulted in a condition number c ( (cid:104) ˜ Φ i , ˜ Φ j (cid:105) ) = .
65. Figure 7 showsthe sample points computed by the accelerated greedy-MPE algorithm. In Fig. 7(a), one can seethat the hyper-reduction approach picks the grid points along the wake and suction-side boundarylayer since these regions contain the most relevant flow dynamics. A detail view of the trailing-18 Number of points C ond i t i on nu m be r Figure 6: Condition number of approximation mode matrix (cid:104) ˜ Φ i , ˜ Φ j (cid:105) obtained from hyper-reduction. edge region is presented in Fig. 7(b) and it is possible to notice that some grid points are alsoselected in the near acoustic field region, where cylindrical sound waves radiate from the trailingedge.Reduced-order models obtained either by Galerkin or LSPG schemes consist in a set of cou-pled non-linear ordinary di ff erential equations. In order to describe the dynamics of the system athand, these equations need to be temporally integrated using a time-marching scheme. A suit ofexplicit and implicit techniques is available and, here, some are tested to assess the accuracy andstability properties of the ROMs. All non-linear least squares problems emerging from implicittime integration are solved using the LevenbergMarquardt algorithm previously described. Themodel time step is selected the same as that between snapshots ∆ t ROM = × − and this valueis 16 times larger compared to the time step used in the FOM simulation. For this time step thehighest frequencies observed in the POD temporal modes are still resolved with at least 20 pointsper wavelength. Carlberg et al. present a detailed theoretical assessment of time integrators forGalerkin and LSPG schemes in [6]. Here, we analyze the performance of implicit and explicitmethods together with calibration techniques. Results are computed for a probe located on theboundary layer at ( x , y ) = (0 . , . u -velocity fluctuations computed by the FOM andthe Galerkin and LSPG models. Solutions are presented for the training region (0 ≤ t ≤
15) andfor an extrapolation period for which FOM solutions are available (15 < t ≤ ff erent time marching schemes and, in Fig. 8, the implicitEuler method is evaluated. This 1st-order scheme is tested due to its stability properties. Ascan be observed in the figure, both Galerkin and LSPG models are stable along the training andextrapolation periods even without calibration, but they show a large amplitude error. Calibrationleads to an excessive damping for the Galerkin method. On the other hand, the LSPG solutionshows a reasonable agreement with the FOM, displaying a small error in amplitude and phase.19 a)(b)Figure 7: Sample points (in blue) chosen by the accelerated greedy-MPE hyper-reduction algorithm. In Fig. 9, the ROMs are integrated using the trapezoidal method. One should expect animprovement in terms of amplitude and phase errors for this 2nd-order implicit scheme, whileretaining stability. For this case, the Galerkin scheme becomes unstable without calibration whilethe LSPG maintains stability, but presents an inaccurate solution. When calibration is employed,both ROMs become stable and relatively accurate, with small discrepancies in the capture ofhigh-frequency oscillations. It is also worth mentioning that calibration render implicit mod-els much cheaper because of the faster convergence rate towards the solution despite imposingadditional terms.Solutions obtained for the explicit fourth-order Runge Kutta (RK4) scheme are presented inFig. 10. Both Galerkin and LSPG solutions become unstable when this explicit time integrationscheme is applied without calibration. However, calibration provides stable Galerkin and LSPGmodels, with the Galerkin scheme being visually more accurate in terms of amplitude.Table 1 shows the computed absolute and root mean-squared errors for the di ff erent timeintegration schemes. The mean-squared error is computed for the fluctuation field which is inte-grated along the entire mesh from 0 ≤ t ≤
30 asError = (cid:107) u (cid:48) FOM ( x , t ) − u (cid:48) ROM ( x , t ) (cid:107) L (cid:107) u (cid:48) FOM ( x , t ) (cid:107) L (42)and, therefore, is prone to any small phase and amplitude di ff erences with respect to the FOMsolution. Results obtained for diagonally-implicit second and third-order Runge Kutta schemes(DIRK2 and DIRK3) are also included in the table. As can be noticed, calibration not only sta-bilizes the solutions but leads to a considerable error reduction. Except for the implicit Eulerscheme employed with the Galerkin model, the mean-squared error is always reduced when cali-bration is applied. However, LSPG combined with implicit Euler performed reasonably well andwould be the preferred method if calibration was not to be used. The trapezoidal method obtainsthe most accurate solutions for both the calibrated Galerkin and LSPG models. For the present20 t -0.3-0.2-0.100.10.2 u FOMG Implicit EulerG Implicit Euler calibrated (a) Galerkin model t -0.3-0.2-0.100.10.2 u FOMLSPG Implicit EulerLSPG Implicit Euler calibrated (b) LSPG modelFigure 8: Time histories of u -velocity fluctuations computed by implicit Euler time marching scheme at ( x , y ) = (0 . , . t -0.3-0.2-0.100.10.2 u FOMG TrapezoidalG Trapezoidal calibrated (a) Galerkin model t -0.3-0.2-0.100.10.2 u FOMLSPG TrapezoidalLSPG Trapezoidal calibrated (b) LSPG modelFigure 9: Time histories of u -velocity fluctuations computed by trapezoidal time marching scheme at ( x , y ) = (0 . , . t -0.3-0.2-0.100.10.2 u FOMG RK4G RK4 calibrated (a) Galerkin model t -0.3-0.2-0.100.10.2 u FOMLSPG RK4LSPG RK4 calibrated (b) LSPG modelFigure 10: Time histories of u -velocity fluctuations computed by RK4 time marching scheme at ( x , y ) = (0 . , . ff ects, a calibrated LSPG model using the RK4 withouthyper-reduction was built. Although better than its low-cost counterpart, the gapless model wasnot only incapable of outperforming the Galerkin results but it was also very expensive due tothe large number of degrees of freedom in the optimization problem. The root mean-squared andabsolute errors for the gapless LSPG-RK4 model are 0 . . Table 1: Mean-square and absolute errors for the di ff erent time integration schemes. Snapshots of u -velocity and pressure fluctuations obtained at t =
24 are presented in Figs.11 and 12, respectively. These figures allow a comparison of results between the FOM and cal-ibrated ROMs using the trapezoidal method for time integration. Although small discrepanciesbetween the ROM and FOM solutions can be observed, especially for small-scale flow struc-tures, the main features of the flow are recovered. In the FOM, vortex merging taking place atthe mid-chord location leads to flow instabilities that develop along the boundary layer. Figure 11shows the velocity fluctuations resulting from these instabilities along the trailing-edge region.As one can see, both the Galerkin and LSPG models are able to capture the relevant flow features.Acoustic scattering occurs due to the advection of the flow instabilities along the trailing edge,what leads to sound waves that propagate upstream closing a feedback loop mechanism [55].Figure 12 shows that the ROMs accurately capture the acoustic waves generated at the trailingedge. In this figure, one can also observe a hydrodynamic wavepacket that is generated at theairfoil mid-chord being transported along the boundary layer towards the trailing edge.
The second test case investigated is a supersonic flow over a NACA0012 airfoil. The freestreamMach number is set as M ∞ = .
2, the Reynolds number is Re c = u -velocity component in Fig. 13. The first mode shown in this figure is solely24 .6 0.7 0.8 0.9 1.0 1.1 1.20.20.10.00.10.2 Full Order Model, t = 24 (a) FOM
Galerkin Calibrated Trapezoidal, t = 24 (b) Galerkin
LSPG Calibrated Trapezoidal, t = 24 (c) LSPGFigure 11: Contours of u -velocity computed at t =
24 using calibrated models with trapezoidal integration.
Full Order Model, t = 24 (a) FOM
Galerkin Calibrated Trapezoidal, t = 24 (b) Galerkin
LSPG Calibrated Trapezoidal, t = 24 (c) LSPGFigure 12: Pressure contours computed at t =
24 using calibrated models with trapezoidal integration. a) Mode 1. (b) Mode 5.(c) Mode 10. (d) Mode 15.Figure 13: Contours of POD eigenfunctions for u -velocity. related to the shock waves and is similar to a mean flow where the detached shock appears sta-tionary upstream the airfoil. Modes 5 and 10 exhibit an oscillatory behavior upstream the airfoilto represent the motion of the bow shock. In these modes, oblique shock waves are observedon the trailing edge and they also have an oscillatory pattern to model the initial transient ofthe fish-tail shock. Downstream the airfoil, one can also see some oscillations that represent theadvection of the starting vortex. Mode 15 shows the oscillatory pattern of vortex shedding whichis formed at the trailing edge after the initial transient.Results obtained by the FOM are sampled for 10 ,
000 snapshots with a constant non-dimensionaltime step ∆ t snapshot = .
01. The first 7 ,
000 snapshots, which represent 70 non-dimensional timeunits, are used to construct a reduced-order basis by the snapshot-POD method. This period issu ffi cient to capture the transient motion of the detached shock wave. The ROMs are built us-ing the first 30 POD modes which contain 99 .
91% of the model RIC. The Galerkin models arebuilt using the maximum calibration parameter ( θ =
1) and this choice is made after an analysisof Fig. 14 that shows the low intrusiveness of the calibration terms relatively to the originalGalerkin operators. For θ =
1, the error E c is minimized and the ratio of the Frobenius normsfrom the calibrated and Galerkin terms is still lower than unity.The application of hyper-reduction for the present transient flow is a challenge since the entirepathway of motion of the detached shock wave has to be captured by the sampled points. More-over, the initial transient also includes the starting vortex that is advected along the airfoil wake.Figure 15 shows the sampled points chosen by di ff erent levels of the accelerated greedy-MPE26 E -5 -4 -3 -2 -1 Figure 14: Ratio of Frobenius norms computed for calibrated and Galerkin coe ffi cients as a function of approximationerror for di ff erent values of θ (see Eq. 34 for details).(a) 20,000 points. (b) 50,000 points. (c) 100,000 points.Figure 15: Sample points (in blue) chosen by the accelerated greedy-MPE hyper-reduction algorithm. hyper-reduction. As shown in Fig. 15(a), only the initial pathway of the bow shock appears infirst 20,000 points selected by the hyper-reduction. At the same time, the stronger oblique shockon the airfoil suction side is better represented by the sampled points. Figure 15(b) shows thata less aggressive hyper-reduction with 50,000 points captures a wider pathway of the detachedshock and adds sample points to both oblique shocks at the trailing edge. Finally, when 100,000points (out of ≈ ff ective as for the previous subsonic flow since thedetached shock wave moves upstream. Hence, in order to capture the entire motion of the shock,the hyper-reduction needs to include its entire pathway ahead of the airfoil, increasing the overallcost of the LSPG ROM.In Fig. 16, results are presented for calibrated ROMs computed with the explicit fourth-orderRunge-Kutta (RK4) time integrator for the LSPG method with and without hyper-reduction andfor the Galerkin method. Without calibration, all ROMs obtained by the RK4 became unstable.The LSPG method is presented for a hyper-reduction including 100,000 sample points. For moreaggressive applications of the accelerated greedy-MPE algorithm, the LSPG method could notprovide stable ROMs. In the figure, temporal modes 1, 5, 10 and 15 are shown and, as observed27efore, the first three modes are associated with the shock motion, being well recovered by bothGalerkin and LSPG methods without hyper-reduction. When the LSPG method is tested withhyper-reduction, a small phase and amplitude distortion is observed for modes 5 and 10. Mode1 appears as a moving average since the present flow is transient and a true mean flow is notobtained for the entire period of simulation. Modes 5 and 10 are composed of low frequencieswhich have some modulation while mode 15 is governed by the vortex shedding dynamics. Thislatter mode is composed of much higher frequencies (compared to modes 5 and 10) with aninitial growth and modulation. For the initial time period (0 < t < t >
10, the calibrated Galerin model is able recover the periodical dynamics of thismode while the LSPG model has a more dissipative behavior and cannot reproduce the high-frequency oscillations. In the Galerkin model, the delay in the vortex shedding developmentmay occur because the initial flow period coincides with the motion of the shock waves which areassociated with more important singular values of the POD decomposition. Once the shocks arereaching their steady-state, the calibrated model is able to capture the less energetic fluctuationsassociated with vortex shedding.The RK4, implicit Euler and trapezoidal time-marching methods are tested for the presentcase. We observed that all non-calibrated Galerkin and LSPG ROMs were unstable except byone, the LSPG ROM computed with the implicit Euler method, which was stable but inaccurate.Figure 17 presents the temporal histories of two probes that capture the u -velocity fluctuations ofthe detached and trailing-edge shock waves. Resuts are provided for the FOM and the Galerkinand LSPG models computed with the RK4 scheme. Solutions obtained by the LSPG with dif-ferent levels of hyper-reduction are also provided for the trapezoidal scheme. All results arepresented for the training region (0 ≤ t ≤
70) and an extrapolation period for which the FOMsolutions are available (70 < t < ff ectedby the hyper-reduction. When 20,000 points are used in the hyper-reduction, the probe missesthe detached shock because the points sampled by the accelerated greedy MPE algorithm do notcontain the entire pathway of the shock wave. The solution obtained without hyper-reduction isconsiderably improved but is still less accurate than that computed by the RK4 LSPG.The u -velocity time histories of vortex shedding are exhibited in Fig. 18 for the probe locationindicated by the yellow symbol in Fig. 19. This probe captures both the motion of the separationbubble on the suction side of the airfoil and the trailing-edge shock motion induced by vortexshedding. In this case, solutions are presented in detail views for the first 15 and last 5 time units.Here, both the Galerkin and LSPG ROMs are not able to reproduce the initial dynamics of theFOM, as can be seen in Fig. 18(a). However, such dynamics is captured by the Galerkin modelfor t >
10 leading to vortex shedding with a small phase and amplitude error as shown in Fig.18(b). In contrast, the LSPG ROM is not able to reproduce the higher frequency vortex sheddingdynamics.Figure 19 presents contours of divergence of velocity, in gray scale, and vorticity, in color, at t =
10 20 30 40 50 60 70 t -3-2-1012 a t -1-0.500.51 a t -0.4-0.200.20.4 a t -0.1-0.0500.050.1 a Figure 16: Temporal dynamics of POD modes 1, 5, 10 and 15. The LSPG is tested with (10 sample points) and withouthyper-reduction.
20 40 60 80 100 t -0.100.10.20.30.4 u (a) Probe location indicated by the purple symbol in Fig.19. t -0.8-0.6-0.4-0.200.2 u (b) Probe location indicated by the green symbol in Fig.19. t -0.100.10.20.30.4 u (c) Probe location indicated by the purple symbol in Fig.19. t -0.8-0.6-0.4-0.200.2 u (d) Probe location indicated by the green symbol in Fig.19.Figure 17: Time histories of u -velocity fluctuations for probes placed at the shock wave locations. t -1.5-1-0.500.511.5 u FOMG RK4LSPG RK4 (a)
95 95.5 96 96.5 97 97.5 98 98.5 99 99.5 100 t -1-0.500.51 u FOMG RK4LSPG RK4 (b)Figure 18: Time histories of u -velocity fluctuations computed for the probe indicated by the yellow symbol in Fig. 19. a) FOM (b) Galerkin RK4 model(c) LSPG RK4 model without hyper-reduction (d) LSPG trapezoidal model with hyper-reductionFigure 19: Contours of dilatation in gray scale and vorticity in color. The purple, green and yellow symbols represent theprobe locations from Figs. 17 and 18. fluctuations at the trailing-edge shock and vortex shedding are dissipated by the LSPG model butare captured by the Galerkin ROM as can be seen in the figure. The LSPG ROM obtained by thetrapezoidal time-marching scheme is also shown and, for this case, hyper-reduction with 20,000points is employed in the model reconstruction. Although the oblique shocks are well capturedby the model, the detached shock wave is positioned downstream compared to the FOM solution.This occurs because the hyper-reduction sampled points do not contain the entire pathway of theshock wave as shown in Fig. 15 and, hence, the model is not able to recover upstream informationof the shock motion.
4. Conclusions
Galerkin and least-squares Petrov-Galerkin reduced-order models are applied to compress-ible flows involving a disparity in temporal scales. The first case investigated comprises a sub-sonic flow past a NACA0012 airfoil for which boundary layer hydrodynamic instabilities lead totrailing-edge noise generation. For this case, a separation bubble induces a low-frequency mod-ulation of the higher-frequency boundary layer instabilities that, in turn, leads to the appearanceof multiple secondary tones in the acoustic radiation. The second case consists of a supersonicflow over a NACA0012 for which shock-vortex shedding interaction appears at the trailing edge.The dataset employed in the model construction includes a transient period of the flow wherea detached shock wave propagates upstream the airfoil leading edge. The initial transient alsoincludes the formation and advection of a starting vortex.32alibration is employed to construct stable and accurate models and, here, we propose anew form of calibration for the LSPG method. Moreover, di ff erently from other references,we employ the Levemberg-Marquadt method to solve the minimization problem appearing inthe LSPG models since it solves the problem of Newton’s method when the Hessian is rank-deficient. For both the Galerkin and LSPG methods, calibration is applied only on constant andlinear terms appearing in the set of non-linear ODEs resulting from the ROMs. While the presentLSPG calibration does not depend on input parameters, the calibration of Galerkin ROMs hasa free parameter that balances the model error and the intrusiveness of the calibration terms.For both cases analyzed, an assessment of this free parameter shows that even with maximumintrusiveness the weight of the calibrated coe ffi cients, measured in terms of the Frobenius norm,is still one order of magnitude smaller than that of the original Galerkin coe ffi cients.Di ff erent time-marching schemes are employed with the Galerkin and LSPG ROMs and, forthe subsonic airfoil flow, the non-calibrated LSPG method was able to build stable methods usingthe implicit Euler and trapezoidal schemes. For the supersonic flow, a stable method was alsoobtained by the non-calibrated LSPG method using the implicit Euler scheme. However, for allthese cases, the models were not accurate in comparison with the FOM solutions. Except for theimplicit Euler applied to the subsonic flow, the non-calibrated Galerkin models were unstableindependently of the time-marching scheme. When calibration was applied to the Galerkin andLSPG methods, solutions became stable and accurate. For the first case investigated, di ff erentimplicit and explicit time-marching schemes provided accurate solutions while, for the secondcase, the fourth-order Runge Kutta presented the most accurate results.In order to avoid rational functions of the model variables, the non-conservative compress-ible Navier-Stokes equations are solved in this work. Therefore, the second degree polynomialnonlinearities of the formulation allow pre-computation and storage of constant, linear and non-linear coe ffi cients and, hence, model cost reduction. Since the overall cost of the LSPG models isconsiderably higher than that of Galerkin models, an accelerated greedy missing point estimationhyper-reduction technique is applied to select the most dynamically relevant points in the flowsfor model reconstruction. In the first problem analyzed, the combination of hyper-reduction andcalibration allows accurate and stable LSPG models to be reconstructed using only 6% of thetotal flow information. However, for the second problem, hyper-reduction is not as e ffi cient dueto the transient nature of the flow. In this latter case, the sampled points need to capture theentire pathway of the detached shock wave in order to provide an accurate model. Otherwise,its motion is compromised as shown in the results. In order to capture all transient aspects ofthe flow, at least 30% of the flow information had to be used in the hyper-reduction procedure.For the supersonic flow, calibration was able to account for some of the damage inflicted byhyper-reduction.This work shows that calibration is important to construct long-time stable and accurateGalerkin and LSPG methods. In the present studies, the calibrated Galerkin ROMs providecheaper and more accurate solutions than the LSPG ROMs. Moreover, they do not require hyper-reduction when implemented with the non-conservative form of the Navier-Stokes equations.Both the Galerkin and LSPG models could recover the small hydrodynamic scales and acousticwaves generated in the subsonic airfoil flow. A comparison between ROM and FOM shows thatresults remain accurate during and beyond the training window. Without hyper-reduction, bothmodels also accurately capture the shock waves present in the airfoil supersonic flow. However,the LSPG models present a more dissipative behavior and could not capture the vortex sheddingmechanism along the airfoil wake. On the other hand, the Galerkin models could reproduce thedynamics of the shedding after an initial transient.33 cknowledgments The authors of this work would like to acknowledge Fundac¸ ˜ao de Amparo `a Pesquisa do Es-tado de S˜ao Paulo, FAPESP, for supporting the present work under research grants No. 2013 / / / Appendix A. Galerkin Coe ffi cients Consider the Galerkin coe ffi cients given by the following tensors e , A and N from Eq. 14 andcomputational domain Ω . These coe ffi cients are functions of the spatial basis Φ = [ Φ ζ Φ u Φ v Φ p ] (cid:62) obtained by POD, mean flow field ¯q = [ ¯ ζ ¯u ¯v ¯p ] (cid:62) , initial conditions [ ζ u v p ] (cid:62) , specific heatratio γ , and reference Prandtl ( Pr ), Reynolds ( Re ) and Mach ( Ma ) numbers. It is convenient todecompose e = e i = e ζ i + e ui + e vi + e pi , which leads to the following terms e ζ i = (cid:90) (cid:90) Ω Φ ζ i (cid:32) ¯ ζ ∂ ¯u ∂ x + ¯ ζ ∂ ¯v ∂ y − ¯u ∂ ¯ ζ∂ x − ¯v ∂ ¯ ζ∂ y (cid:33) dx dy (A.1) e ui = (cid:90) (cid:90) Ω Φ u i (cid:34) − ¯u ∂ ¯u ∂ x − ¯v ∂ ¯u ∂ y − ¯ ζ ∂ ¯p ∂ x + MaRe (cid:32) ¯ ζ ∂ ¯u ∂ x − ¯ ζ ∂ ¯v ∂ x ∂ y + ¯ ζ ∂ ¯v ∂ y ∂ x + ¯ ζ ∂ ¯u ∂ y (cid:33)(cid:35) dx dy (A.2) e vi = (cid:90) (cid:90) Ω Φ v i (cid:34) − ¯u ∂ ¯v ∂ x − ¯v ∂ ¯v ∂ y − ¯ ζ ∂ ¯p ∂ y + MaRe (cid:32) ¯ ζ ∂ ¯v ∂ y − ¯ ζ ∂ ¯u ∂ y ∂ x + ¯ ζ ∂ ¯u ∂ x ∂ y + ¯ ζ ∂ ¯v ∂ x (cid:33)(cid:35) dx dy (A.3) e pi = (cid:90) (cid:90) Ω Φ p i (cid:34) − ¯u ∂ ¯p ∂ x − ¯v ∂ ¯p ∂ y − γ (cid:32) ¯p ∂ ¯u ∂ x + ¯p ∂ ¯v ∂ y (cid:33) + γ MaRePr (cid:32) ∂ ¯p ∂ x ¯ ζ + ∂ ¯p ∂ x ∂ ¯ ζ∂ x + ¯p ∂ ¯ ζ∂ x + ∂ ¯p ∂ y ¯ ζ + ∂ ¯p ∂ y ∂ ¯ ζ∂ y + ¯p ∂ ¯ ζ∂ y (cid:33) + ( γ − MaRe (cid:32) ∂ ¯u ∂ x ∂ ¯u ∂ x − ∂ ¯u ∂ x ∂ ¯v ∂ y + ∂ ¯v ∂ y ∂ ¯v ∂ y − ∂ ¯v ∂ y ∂ ¯u ∂ x + ∂ ¯v ∂ x ∂ ¯v ∂ x + ∂ ¯v ∂ x ∂ ¯u ∂ y + ∂ ¯u ∂ y ∂ ¯u ∂ y (cid:33)(cid:35) dx dy (A.4)The term A is also conveniently decomposed as A = A i j = A ζ i j + A ui j + A vi j + A pi j A ζ i j = (cid:90) (cid:90) Ω Φ ζ i (cid:32) ¯ ζ ∂ Φ u j ∂ x + ¯ ζ ∂ Φ v j ∂ y − ¯u ∂ Φ ζ j ∂ x − ¯v ∂ Φ ζ j ∂ y + ∂ ¯ u ∂ x Φ ζ j + ∂ ¯ v ∂ y Φ ζ j − ∂ ¯ ζ∂ x Φ u j − ∂ ¯ ζ∂ y Φ v j (cid:33) dx dy (A.5) A ui j = (cid:90) (cid:90) Ω Φ u i (cid:40) − ¯u ∂ Φ u j ∂ x − ¯v ∂ Φ u j ∂ y − ¯ ζ ∂ Φ p j ∂ x − ∂ ¯ u ∂ x Φ u j − ∂ ¯ u ∂ y Φ v j − ∂ ¯ p ∂ x Φ ζ j + MaRe (cid:34) (cid:32) ¯ ζ ∂ Φ u j ∂ x + ∂ ¯ u ∂ x Φ ζ j (cid:33) − (cid:32) ¯ ζ ∂ Φ v j ∂ x ∂ y + ∂ ¯ v ∂ x ∂ y Φ ζ j (cid:33) + ¯ ζ ∂ Φ v j ∂ y ∂ x + ∂ ¯ v ∂ y ∂ x Φ ζ j + ¯ ζ ∂ Φ u j ∂ y + ∂ ¯ u ∂ y Φ ζ j (cid:35)(cid:41) dx dy (A.6)34 vi j = (cid:90) (cid:90) Ω Φ v i (cid:40) − ¯u ∂ Φ v j ∂ x − ¯v ∂ Φ v j ∂ y − ¯ ζ ∂ Φ p j ∂ y − ∂ ¯ v ∂ x Φ u j − ∂ ¯ v ∂ y Φ v j − ∂ ¯ p ∂ y Φ ζ j + MaRe (cid:34) (cid:32) ¯ ζ ∂ Φ v j ∂ y + ∂ ¯ v ∂ y Φ ζ j (cid:33) − (cid:32) ¯ ζ ∂ Φ u j ∂ y ∂ x + ∂ ¯ u ∂ y ∂ x Φ ζ j (cid:33) + ¯ ζ ∂ Φ u j ∂ x ∂ y + ∂ ¯ u ∂ x ∂ y Φ ζ j + ¯ ζ ∂ Φ v j ∂ x + ∂ ¯ v ∂ x Φ ζ j (cid:35)(cid:41) dx dy (A.7) A pi j = (cid:90) (cid:90) Ω Φ p i (cid:40) − ¯u ∂ Φ p j ∂ x − ¯v ∂ Φ p j ∂ y − ∂ ¯ p ∂ x Φ u j − ∂ ¯ p ∂ y Φ v j − γ (cid:32) ¯p ∂ Φ u j ∂ x + ¯p ∂ Φ v j ∂ y + ∂ ¯ u ∂ x Φ p j + ∂ ¯ v ∂ y Φ p j (cid:33) + γ MaRePr (cid:34) ∂ ¯p ∂ x Φ ζ j + ¯ ζ ∂ Φ p j ∂ x + (cid:32) ∂ ¯ p ∂ x ∂ Φ ζ j ∂ x + ∂ ¯ ζ∂ x ∂ Φ p j ∂ x (cid:33) + ¯p ∂ Φ ζ j ∂ x + ∂ ¯ ζ∂ x Φ p j + ∂ ¯p ∂ y Φ ζ j + ¯ ζ ∂ Φ p j ∂ y + (cid:32) ∂ ¯ p ∂ y ∂ Φ ζ j ∂ y + ∂ ¯ ζ∂ y ∂ Φ p j ∂ y (cid:33) + ¯p ∂ Φ ζ j ∂ y + ∂ ¯ ζ∂ y Φ p j (cid:35) + ( γ − MaRe (cid:34) (cid:32) ∂ ¯ u ∂ x ∂ Φ u j ∂ x + ∂ ¯ u ∂ x ∂ Φ u j ∂ x (cid:33) − (cid:32) ∂ ¯ u ∂ x ∂ Φ v j ∂ y + ∂ ¯ v ∂ y ∂ Φ u j ∂ x (cid:33) + (cid:32) ∂ ¯ v ∂ y ∂ Φ v j ∂ y + ∂ ¯ v ∂ y ∂ Φ v j ∂ y (cid:33) − (cid:32) ∂ ¯ v ∂ y ∂ Φ u j ∂ x + ∂ ¯ u ∂ x ∂ Φ v j ∂ y (cid:33) + (cid:32) ∂ ¯ v ∂ x ∂ Φ v j ∂ x + ∂ ¯ v ∂ x ∂ Φ u j ∂ y + ∂ ¯ u ∂ y ∂ Φ v j ∂ x + ∂ ¯ u ∂ y ∂ Φ u j ∂ y (cid:33)(cid:35)(cid:41) dx dy (A.8)and, finally, the third order tensor is written as N = N i jk = N ζ i jk + N ui jk + N vi jk + N pi jk N ζ i jk = (cid:90) (cid:90) Ω Φ ζ i (cid:32) Φ ζ j ∂ Φ u k ∂ x + Φ ζ j ∂ Φ v k ∂ y − Φ u j ∂ Φ ζ k ∂ x − Φ v j ∂ Φ ζ k ∂ y (cid:33) dx dy (A.9) N ui jk = (cid:90) (cid:90) Ω Φ u i (cid:34) − Φ u j ∂ Φ u k ∂ x − Φ v j ∂ Φ u k ∂ y − Φ ζ j ∂ Φ p k ∂ x + MaRe (cid:32) Φ ζ j ∂ Φ u k ∂ x − Φ ζ j ∂ Φ v k ∂ x ∂ y + Φ ζ j ∂ Φ v k ∂ y ∂ x + Φ ζ j ∂ Φ u k ∂ y (cid:33)(cid:35) dx dy (A.10) N vi jk = (cid:90) (cid:90) Ω Φ v i (cid:34) − Φ u j ∂ Φ v k ∂ x − Φ v j ∂ Φ v k ∂ y − Φ ζ j ∂ Φ p k ∂ y + MaRe (cid:32) Φ ζ j ∂ Φ v k ∂ y − Φ ζ j ∂ Φ u k ∂ y ∂ x + Φ ζ j ∂ Φ u k ∂ x ∂ y + Φ ζ j ∂ Φ v k ∂ x (cid:33)(cid:35) dx dy (A.11)35 pi jk = (cid:90) (cid:90) Ω Φ p i (cid:34) − Φ u j ∂ Φ p k ∂ x − Φ v j ∂ Φ p k ∂ y − γ (cid:32) Φ p j ∂ Φ u k ∂ x + Φ p j ∂ Φ v k ∂ y (cid:33) + γ MaRePr (cid:32) Φ ζ j ∂ Φ p k ∂ x + ∂ Φ p j ∂ x ∂ Φ ζ k ∂ x + Φ p j ∂ Φ ζ k ∂ x + ∂ Φ p j ∂ y Φ ζ k + ∂ Φ p j ∂ y ∂ Φ ζ k ∂ y + Φ p j ∂ Φ ζ k ∂ y (cid:33) + ( γ − MaRe (cid:32) ∂ Φ u j ∂ x ∂ Φ u k ∂ x − ∂ Φ u j ∂ x ∂ Φ v k ∂ y + ∂ Φ v j ∂ y ∂ Φ v k ∂ y − ∂ Φ v j ∂ y ∂ Φ u k ∂ x + ∂ Φ v j ∂ x ∂ Φ v k ∂ x + ∂ Φ v j ∂ x ∂ Φ u k ∂ y + ∂ Φ u j ∂ y ∂ Φ v k ∂ x + ∂ Φ u j ∂ y ∂ Φ u k ∂ y (cid:33)(cid:35) dx dy (A.12)The initial condition a of Eq. 14 is obtained by projection of the first snapshot on the basisvector as a i = (cid:90) (cid:90) Ω (cid:16) ( ζ − ¯ ζ ) Φ ζ i + ( u − ¯u ) Φ u i + ( v − ¯v ) Φ v i + ( p − ¯p ) Φ p i (cid:17) dx dy . (A.13) References [1] W. R. Wolf, J. L. F. Azevedo, S. K. Lele, Convective e ff ects and the role of quadrupole sources for aerofoilaeroacoustics, Journal of Fluid Mechanics 708 (2012) 502538.[2] T. R. Ricciardi, W. R. Wolf, R. Speth, Acoustic prediction of LAGOON landing gear: Cavity noise and coherentstructures, AIAA Journal 56 (2018) 4379–4399.[3] B. J. Olson, S. K. Lele, A mechanism for unsteady separation in over-expanded nozzle flow, Physics of Fluids 25(2013) 110809.[4] M. Bergmann, L. Cordier, J.-P. Brancher, Optimal rotary control of the cylinder wake using proper orthogonaldecomposition reduced-order model, Physics of Fluids 17 (2005) 097101.[5] D. Amsallem, M. J. Zahr, C. Farhat, Nonlinear Model Order Reduction Based on Local Reduced-Order Bases,International Journal for Numerical Methods in Engineering 92 (2012) 891–916.[6] K. Carlberg, M. Barone, H. Antil, Galerkin v. least-squares Petrov-Galerkin projection in nonlinear model reduc-tion, Journal of Computational Physics 330 (2017) 693–734.[7] M. Bergmann, A. Ferrero, A. Iollo, E. Lombardi, A. Scardigli, H. Telib, A zonal Galerkin-free POD model forincompressible flows, Journal of Computational Physics 352 (2018) 301–325.[8] S. J. Grimberg, C. Farhat, N. Youkilis, On the stability of projection-based model order reduction for convection-dominated laminar and turbulent flows, Journal of Computational Physics 419 (2020) 1–28.[9] C. W. Rowley, T. Colonius, R. M. Murray, Model reduction for compressible flows using POD and Galerkinprojection, Physica D: Nonlinear Phenomena 189 (2004) 115 – 129.[10] K. Carlberg, C. Farhat, D. Amsallem, The GNAT method for nonlinear model reduction: E ff ective implementationand application to computational fluid dynamics and tubulent flows, Journal of Computational Physics 243 (2013)623–647.[11] L. Sirovich, Turbulence and the Dynamics of Coherent Structures. Part 1: Coherent Structures, Quarterly ofApplied Mathematics 45 (1987) 561–571.[12] L. Cordier, M. Bergmann, Proper Orthogonal Decomposition: an overview, in: Lecture series 2002-04, 2003-03and 2008-01 on post-processing of experimental and numerical data, Von Karman Institute for Fluid Dynamics,2008., VKI, 2008, p. 46 pages.[13] S. Pan, K. Duraisamy, Long-time predictive modeling of nonlinear dynamical systems using neural networks,Complexity 2018 (2018) 1–26.[14] Z. Y. Wan, P. Vlachas, P. Koumoutsakos, T. Sapsis, Data-assisted reduced-order modeling of extreme events incomplex dynamical systems, PLOS ONE 13 (2018) 1–22.
15] P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, P. Koumoutsakos, Data-driven forecasting of high-dimensionalchaotic systems with long short-term memory networks, Proceedings of the Royal Society of London A: Mathe-matical, Physical and Engineering Sciences 474 (2018).[16] S. L. Brunton, J. L. Proctor, N. J. Kutz, Discovering governing equations from data by sparse identification ofnonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (2016) 3932–3937.[17] H. F. S. Lui, W. R. Wolf, Construction of reduced-order models for fluid flows using deep feedforward neuralnetworks, Journal of Fluid Mechanics 872 (2019) 963994.[18] V. Zucatti, H. F. S. Lui, D. B. Pitz, W. R. Wolf, Assessment of reduced-order modeling strategies for convectiveheat transfer, Numerical Heat Transfer, Part A: Applications 77 (2020) 702–729.[19] W. Cazemier, R. W. C. P. Verstappen, A. E. P. Veldman, Proper orthogonal decomposition and low-dimensionalmodels for driven cavity flows, Physics of Fluids 10 (1998) 1685–1699.[20] B. R. Noack, P. Papas, P. A. Monkewitz, The need for a pressure-term representation in empirical Galerkin modelsof incompressible shear flows, Journal of Fluid Mechanics 523 (2005) 339–365.[21] M. Bergmann, C. H. Bruneau, A. Iollo, Enablers for robust POD models, Journal of Computational Physics 228(2009) 516–538.[22] O. San, T. Iliescu, Proper orthogonal decomposition closure models for fluid flows: Burgers equation, InternationalJournal of Numerical Analysis and Modeling, Series B 5 (2014) 217–237.[23] A. Iollo, S. Lanteri, J.-A. Dsidri, Stability Properties of POD-Galerkin Approximations for the CompressibleNavier-Stokes Equations, Theoretical and Computational Fluid Dynamics 13 (2000) 377–396.[24] Z. Wang, I. Akhtar, J. Borggaard, T. Iliescu, Proper orthogonal decomposition closure models for turbulent flows:A numerical comparison, Comput. Methods Appl. Mech. Eng. 237 (2012) 10–26.[25] B. Galletti, A. Bottaro, C. H. Bruneau, A. Iollo, Accurate model reduction of transient and forced wakes, EuropeanJournal of Mechanics B / Fluids 26 (2006) 354–366.[26] M. Couplet, C. Basdevant, P. Sagaut, Calibrated reduced-order POD-Galerkin system for fluid flow modelling,Journal of Computational Physics 207 (2005) 192–220.[27] R. Bourguet, M. Braza, A. Dervieux, Reduced-order modeling for unsteady transonic flows around an airfoil,Physics of Fluids 19 (2007) 111701.[28] R. Bourguet, M. Braza, A. Svrain, A. Bouhadji, Capturing transition features around a wing by reduced-ordermodeling based on compressible NavierStokes equations, Physics of Fluids 21 (2009) 094104.[29] J. Favier, Contrle d’coulements : approche exprimentale et modlisation de dimension rduite, Ph.D. thesis, 2007.[30] J. Favier, L. Cordier, A. Kourta, A. Iollo, Calibrated POD Reduced-Order Models of Massively Separated Flowsin the Perspective of Their Control, Fluids Engineering Division Summer Meeting, pp. 743–748.[31] C. Farhat, P. Avery, T. Chapman, J. Cortial, Dimensional reduction of nonlinear finite element dynamic modelswith finite rotations and energy-based mesh sampling and weighting for computational e ffi ciency, InternationalJournal for Numerical Methods in Engineering 98 (2014) 625–662.[32] S. J. Grimberg, C. Farhat, Hyperreduction of CFD Models of Turbulent Flows using a Machine Learning Approach.[33] M. Balajewicz, I. Tezaur, E. Dowell, Minimal subspace rotation on the Stiefel manifold for stabilization andenhancement of projection-based reduced order models for the compressible Navier-Stokes equations, Journal ofComputational Physics 321 (2016) 224–241.[34] D. Amsallem, C. Farhat, Interpolation Method for Adapting Reduced-Order Models and Application to Aeroelas-ticity, AIAA Journal 46 (2008).[35] A. M. A. Quarteroni, F. Negri, Reduced Basis Methods for Partial Di ff erential Equations: An Introduction,Springer, Switzerland, 1st edition edition, 2016.[36] J. Nocedal, S. J. Wright, Numerical Optimization, Springer, 2006.[37] K. Carlberg, C. Bou-Mosleh, C. Farhat, E ffi cient non-linear model reduction via a least-squares petrovgalerkinprojection and compressive tensor approximations, International Journal for Numerical Methods in Engineering86 (2011) 155–181.[38] B. Kramer, K. E. Willcox, Nonlinear model order reduction via lifting transformations and proper orthogonaldecomposition, AIAA Journal 57 (2019) 2297–2307.[39] R. Everson, L. Sirovich, Karhunun-Love procedure for gappy data, Optical Society of America 12 (1995) 1657–1664.[40] K. Willcox, Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition, Computers &Fluids 35 (2006) 208–226.[41] P. Astrid, S. Weiland, K. Willcox, T. Backx, Missing point estimation in models described by proper orthogonaldecomposition, IEEE Transactions on Automatic Control 53 (2008) 2237–2251.[42] S. Chaturantabut, D. C. Sorensen, Nonlinear Model Reduction via Discrete Empirical Interpolation, SIAM Journalon Scientific Computing 32 (2010) 2737–3764.[43] B. Peherstorfer, K. Willcox, Online adaptive model reduction for nonlinear systems via low-rank updates, SIAMJ. Sci. Comput. 37 (2015) A2123–A2150.