On the need of mode interpolation for data-driven Galerkin models of a transient flow around a sphere

We present a low-dimensional Galerkin model with state-dependent modes capturing linear and nonlinear dynamics. Departure point is a direct numerical simulation of the three-dimensional incompressible flow around a sphere at Reynolds numbers 400. This solution starts near the unstable steady Navier–Stokes solution and converges to a periodic limit cycle. The investigated Galerkin models are based on the dynamic mode decomposition (DMD) and derive the dynamical system from first principles, the Navier–Stokes equations. A DMD model with training data from the initial linear transient fails to predict the limit cycle. Conversely, a model from limit-cycle data underpredicts the initial growth rate roughly by a factor 5. Key enablers for uniform accuracy throughout the transient are a continuous mode interpolation between both oscillatory fluctuations and the addition of a shift mode. This interpolated model is shown to capture both the transient growth of the oscillation and the limit cycle.


Introduction
In this study, we propose enablers for reduced-order models of transient flows. Focus is placed on data-driven Galerkin models based on solutions of computational fluid mechanics (CFD). In contrast to most experimental techniques, CFD provides information about the whole flow field. This information enables the derivation of dynamic Galerkin systems from first principles, i.e., the Navier-Stokes equations.
One driving force behind such reduced-order models (ROM) is the need for hundred of thousands of costly high-fidelity CFD simulations for new aircraft design [37]. The large number of test cases includes variations of the aircraft configuration, mass distribution, flow conditions (e.g., gusts, angle of attack) besides the whole flight envelope from starting to landing. Such a number of simulations is very challenging even in the case of simplified flow models (like URANS) and massive parallel simulations. ROM may provide cheap substitute plants for a range of operating conditions, thus significantly reducing the computational cost [16]. Another driving force is flow control. Simple control-oriented ROM are key enablers for model-based control design [8,40]. The application of linear control design methods to a CFD discretization constitutes, for instance, a immense numerical challenge.
Myriad of successful ROMs have been developed for a large range of fluid flows [36,53]. Most ROMs fall in the category of data-driven Galerkin models [17], the focus of this study. These Galerkin models are based on empirical modes like proper orthogonal decomposition (POD) [17] or dynamic mode decompositions (DMD) [39,42]. The dynamical system is derived from first principles via Galerkin projection [13] or from data with powerful methods of 4D data assimilation [9] or combinations thereof, e.g., the linear term from a fit and the quadratic term from Galerkin projection [14]. With the current powerful arsenal of methods, it is difficult to conceive flows whose essential features cannot be reproduced by a ROM in a limited observation domain.
However, transients and multiple operating conditions still constitute a challenge to ROM. Deane et al. [12] already noted the fragility of their model to small changes in their computational procedure. One cause is that the stability of the Navier-Stokes solution is not preserved under Galerkin projection. In an ingenious example, Rempfer showed how a stable fixed point can become unstable upon Galerkin projection in a lower-dimensional subspace [35]. The authors have identified additional modes for base flow variations as one root cause of fragile ROM [30,32]. The proposed shift modes were inspired by the structure of the pioneering wall turbulence model of Aubry et al. [4] and mean-field theory [49].
On a kinematic side, the coherent structures change during change of operating conditions, as addressed in sequential POD [20], or during transients which can be conceptualized as sequence of operating conditions. In a simple lookup table approach, different stages of the transient can be catalogued as one operating condition with corresponding modal expansion [24]. Double POD [43] addresses the change of modes in similar spirits. A more elegant approach is the construction of low-dimensional manifolds with associated continuously deforming modes [5]. This approach has been applied to the cylinder wake but comes with significant computational cost. In this study, we pursue continuous mode interpolation [27] as simple computationally inexpensive and promising path. As benchmark, we choose the incompressible three-dimensional flow around a sphere as prototype of a three-dimensional bluff-body wake.
The paper is organized as follows. Section 2 contains the description of the flow past a sphere and the details of the employed direct numerical simulation. The empirical Galerkin method is reviewed in Sect. 3. Galerkin models of the flow past a sphere, based on fixed modal bases resulting from limit-cycle oscillations and fixed point dynamics, are presented, respectively, in Sects. 4 and 5. In Sect. 6, a Galerkin model is proposed which interpolates between both bases. Thus, the whole transient is described with good accuracy. The results are discussed in Sect. 7.

Direct numerical simulation of a transient flow past a sphere
The flow past a sphere may be considered as a prototype wake of axisymmetric bluff bodies [22] or even of immersed bluff bodies [19] which are frequently encountered in many engineering applications [41]. The symmetric body may give rise to non-symmetric three-dimensional unsteady wake by symmetry braking, just like for the circular cylinder wake. Complex wake patterns may even be observed at moderate Reynolds numbers and complex vortical interactions. This kaleidoscope of phenomena has made the flow around a sphere a well-studied benchmark.
Experimental investigations [52] indicate that the separation from the rear of the sphere occurs at Re ≈ 24. The resulting axisymmetric vortex ring is observed up to Re = 210. In the range 210 < Re < 270, the flow becomes non-axisymmetric and the ring vortex is shifted off-axis [19,26]. Numerical simulations done by Tomboulides et al. [54] have corroborated these results. The deadwater zone is first observed at Re = 20, and steady axisymmetric flow prevails at Re < 212. At Reynolds number Re ≈ 280 [19,26,29,52], an unsteady double-thread wake appears leading to hairpin vortices. The frequencies of vortex shedding change with Reynolds number from St = 0.126 at Re = 250 to St = 0.167 [54] or St ≈ 0.18 [41] at Re = 500.
In this paper, a viscous, incompressible flow around a sphere is considered at Reynolds number Re = 400 based on diameter. The flow is described in a Cartesian coordinate system x = (x, y, z) with the origin in the center of the sphere and the x-axis aligned with the flow. Similarly, the velocity is presented in analog Cartesian components u = (u, v, w). The time and pressure are denoted by t and p. The fluid has constant density ρ and kinematic viscosity ν. In the following, all quantities are assumed to be non-dimensionalized The flow evolution is described by the Navier-Stokes equations augmented by boundary and initial conditions. The spatial discretization of these equations is performed by a second-order finite element method with Galerkin formulation and quadratic interpolation functions. The computational domain extends 20 diameters downstream of the sphere. The cross section of the domain is a square of size 10 × 10 diameters. The domain is discretized ( Fig. 1) using 353,310 second-order tetrahedral elements and 480,472 nodes. Boundary conditions include constant velocity (inlet) at front, top, bottom and side faces, free-stream condition (outlet) at back face of the domain and zero velocity ("no slip wall") on the sphere. The numerical simulation is performed using an in-house Navier-Stokes solver [28]. A penalty formulation (LBB condition [7]) eliminates the pressure term. The time integration employs an iterative Newton-Raphson approach. For a parallel computation, the mesh is decomposed onto several sub-domains using the Metis tool [21]. The communication between the computational nodes is done using a Message Passing Interface (MPI) library.
The computation starts from steady solution at Re = 400 (Fig. 2, left) and advances until a stable limit cycle is reached. An instantaneous solution for the post-transient periodic flow is depicted in Fig. 3. These oscillations are characterized by counter-rotating vortices stretched in the streamwise direction [10]. The vortex shedding frequency is approximately St = 0.134. Time-averaged flow (Fig. 2, right), calculated for limit-cycle oscillations, is characterized by slightly non-axisymmetric separation, shorter than in the case of steady solution.

Model order reduction based on dynamic mode decomposition (DMD)
Model order reduction techniques aim to derive a dynamical system of a few equations that reproduces the dynamics of the full-scale, high-fidelity model-at least for the training data. A prominent example is the Galerkin method [13] in which state variables (like velocities u) are approximated by a base flow u 0 , for example, a steady solution or time-averaged flow and an expansion in terms of N modes u j : The spatial modes u i (x) allow to express the flow state compactly in terms of N mode amplitudes a i (t). The resulting data compression are typically many orders of magnitudes. The truncation of the expansion to the N most dominant modes leads to a non-vanishing residuum R [N ] in the Navier-Stokes equations. In the "traditional" Galerkin method [13], the residuum R [N ] is projected onto the space spanned by the expansion modes (3): This approach leads to a system of ordinary differential equations for the mode amplitudes a i (t): The coefficients of the viscous (l ν i j ) and convective term (q c i jk ) are given by The Galerkin pressure term representation f p i = −(u i , p) Ω can typically be neglected for a large observation domain of absolutely unstable wake flows [32]. The impact of the truncation error can be reduced by further calibration of the reduced-order model [11,47].
The expansion modes can be classified in terms of mathematical, physical and empirical approaches, and the right choice of the modal basis has a significant impact on the scope of applicability of the model [30]. Up to now, empirical modes, mostly linear combinations of the snapshot data, are the most widely used. The first popular method was proper orthogonal decomposition [17] and its variants like method of snapshots [44], bi-orthogonal decomposition (BOD) [3], double POD (DPOD) [43] and balanced POD (BPOD) [38]. A more recent alternative is the dynamic mode decomposition (DMD) [39,42] and its variants like oscillation pattern decomposition (OPD) [56], multi-resolution DMD [23] and recursive DMD [33].
DMD can approximate stability eigenmodes from transient data [45] and Fourier modes from post-transient data. In this method, the unsteady data are represented by a time-resolved sequence of n snapshots V 0...n = {q 0 , q 1 , q 2 , . . . , q n } sampled with constant frequency f s = 1/Δt. These snapshots are assumed to be linearly dependent: where the linear operator is defined asÃ = e t A . Here, the companion matrix S is given by The parameters c = [c 0 , . . . , c n ] T are obtained from the solution of the overdetermined system of Eq. (6), by the minimization of the residual norm V 0...n C − V n+1 2 .
The eigenvectors b i of matrix S are projected onto the input snapshots in order to obtain the (complex) DMD modes v i , while the eigenvalues λ i determine modal growth ratios (real part) and angular frequencies (imaginary parts): In the following, also f i = ω i /2π are used.
In the following, the DMD modes are assumed to be sorted by magnitude: The Galerkin expansion modes are built from real DMD modes (for real eigenvalues λ i = σ i ) or the real and imaginary parts of complex conjugate DMD mode pairs (for complex λ i,i+1 = σ i ± ıω i ). The employed modes u i in (2) are Gram-Schmidt orthonormalized.

Reduced-order model for the periodic flow
We use the snapshots of six post-transient periods for the dynamic mode decomposition (DMD) of the snapshots taken from 6 cycles of fully periodic flow that has been done (Fig. 4). The resulting DMD modes are sorted by magnitude as discussed in Sect. 3. Figure 4 displays the corresponding eigenvalues and norms.
The eigenfrequencies and growth rates are resolved in Fig. 5. The first mode has the primary flow frequency ( f 1 = St = 0.134). Higher-order modes exhibit multiples of this frequency in a monotonically increasing order. The growth rates almost vanish as predicted by Koopman theory for long-term post-transient data [39]. The first three DMD modes are presented in Figs. 6, 7 and 8. The structures slowly grow with dominant spatial periodicity in streamwise direction. The streamwise wave number increases in proportion to the frequency.
Orthonormalized DMD modes captured from limit cycle obviously resemble POD modes, shown in Fig. 9. While DMD modes yield the same subspace as POD modes, they point in different directions-what is observed as the rotation of the modes around x-axis. More details on the relationship between POD and DMD might be found in a paper by Tu et al. [55], and the comparison of Galerkin models of 2D cylinder wake, based on POD and DMD modes, might be found in [45].
A six-dimensional Galerkin model is build from the depicted DMD modes as discussed in Sect. 3. For easy reference, this ROM will be called Galerkin model A or (GM-A). The mode amplitudes of DNS are well reproduced by the ROM within few percent error, as expected from earlier similar studies starting with Deane et al. [12]. The fluctuation energy is slightly overpredicted as the energy flow to neglected higher-order modes is not accounted by an eddy viscosity model.   Next, the ROM is test for off-design conditions, the transient from steady solution to limit cycle (Sect. 2). The initial condition of the mode amplitudes is obtained from straightforward Galerkin approximation: Evidently, the ROM underpredicts the growth rate. The limit cycle is researched after about 250 units as opposed to 50. This was to be expected as the onset of oscillations is described by the dominant stability eigenmodes which have, by construction, a larger growth rate than DMD modes. One reason of GM-A poor reproduction of the transients is the unresolved base flow change. This deficit is cured by augmenting the model with the shift mode [30], the difference between averaged and steady solution orthonormalized with respect to 6 modes. The shift mode is also known as zeroth mode [34] and depicted in Fig. 11. The resulting seven-dimensional model, called GM-B for brevity, exhibits a larger growth rate at the expense of a significant overshoot (Fig. 12). The overshoot also be observed in cylinder wake flows for the analog ROM and is explained in [51]. Our study corroborates a frequent observation that empirical Galerkin models can well reproduce periodic flows but are challenged for off-design conditions, like transients.

Reduced-order model for the linear dynamics near the fixed point
In this section, we construct a ROM, referred to as Galerkin model C, which is optimized for the linear dynamics near the fixed point. The training data comprise the first 300 snapshots of the transient discussed in Sect. 2.
The linear dynamics is governed by one exponentially growing oscillatory stability mode in accord with the literature. The direct computation of this mode requires the solution of huge, non-Hermitian eigenvalue problems. Global stability analysis is still hardly performable for three-dimensional flows, even for simple geometries without symmetries [53]. DMD offers a convenient data-driven alternative to distill the most dominant stability modes. The corresponding eigenvalues and norms of the DMD modes are depicted in   A Galerkin model based on these DMD modes diverges. Hence, we add the stabilizing shift mode to obtain a seven-dimensional Galerkin model C, the analog of model B for the limit-cycle data. The transient of GM-C (see Fig. 18) has a larger initial growth rate, but it is unable to predict the asymptotic fluctuation energy and

Reduced-order model based on continuous mode interpolation
In this section, we integrate the good performances of Galerkin model B and C in their respective ranges by a continuous mode interpolation (CMI) of the 6 oscillatory expansion modes. This interpolation has been applied to a transient of a cylinder wake [27,48] and to a flow around an airfoil with changing Reynolds number and varying angles of attack [46]. Starting point of CMI is the Fredholm eigenproblem in space domain where the kernel R is the cross-correlation matrix of the velocity components between locations x and y. Let κ be the interpolation with κ = 0 referring to one state, like the beginning of a transient and κ = 1 referring with eigenvectors u κ i . Thus, the eigenvectors at states κ = 0 and κ = 1 can be linearly interpolated with proper tracking of the eigenvalue.
We "purify" the kernel from neglected eigenmodes by reconstructing the initial and final functions from where λ κ i , i = 1, . . . , N represents the eigenvalues of (10). In this study, we use the fluctuation energy K (9) as interpolation parameter K 0 and K 1 are the initial and post-transient levels of fluctuation energy. A good approximation is K 0 ≈ 0 and thus κ = K (t)/K 1 . K measures the magnitude of the Reynolds stress and thus, to first order, also the mean-field deformation from the fixed point or, equivalently, the shift-mode amplitude.
If κ exceeds 1, its value is averaged using the previous value of κ. The resulting Galerkin expansion starts with the unstable steady solution as base flow, contains the interpolated oscillatory modes and includes the shift mode as seventh mode: where u 0 = u s , u κ 7 = u 7 , a 0 ≡ 1 and N = 7. The oscillatory modes i = 1, . . . , 6 are interpolated and are orthonormal by construction: Orthogonality is guaranteed by the symmetry of the kernel, and orthonormality can be enforced by normalization. The steady solution u s and shift mode u 7 span the base flow approximation. The shift mode is exactly orthogonal to antisymmetric modes and approximately orthogonal to the remaining oscillatory modes. The corresponding Galerkin system (4)  The coefficients l ν,κ i j , q c,κ i jk become functions of κ, i.e., the fluctuation energy. Figure 19 displays the transient solution of the resulting Galerkin model D (GM-D). The new model significantly outperforms GM-B and GM-C as expected from construction. Figure 20 shows the corresponding phase portraits in good agreement with the DNS data.

Conclusions and outlook
This study concerns model order reduction for an incompressible, three-dimensional, viscous flow past a sphere at Reynolds number 400. The flow around sphere is a well-investigated textbook prototype of 3D bluff-body wakes. Focus is placed on the transient dynamics from an unstable steady solution to the posttransient limit cycle. This transient flow constitutes a challenge for reduced-order models (ROM) as the base flow, and the associated coherent structures change considerably. The snapshot data are obtained by a direct numerical solution of the Navier-Stokes equations. One motivation for this transient comes from need of a control-oriented ROM to describe both unforced and forced transients [8].
Starting point are low-order Galerkin models based on the most dominant DMD modes. We emphasize that the dynamical system is derived from first principles, while the expansion is extracted from snapshot data. As expected from other wake studies, the DMD model incorporating the first 3 harmonics (6 modes) describes well the limit cycle with corresponding training data but underpredicts the transient time by roughly a factor of 5. Conversely, the DMD model optimized for the initial linear dynamics is far more accurate but diverges in the long term. Even the inclusion of a shift mode does not cure a significant overshoot. These results are well aligned with earlier cylinder wake studies starting with the pioneering POD model of Deane et al. [12] and subsequent later refinements [30]. Here, DMD of the initial transient features stability eigenmodes and DMD of the limit cycle replaces POD modes.
In this study, we address the coherent structure change from linear to nonlinear dynamics with a continuous mode interpolation [27] between the initial and final DMD expansion of sixth order. This interpolated expansion resolves the first, second and third harmonics. In addition, the shift mode is added as seventh mode to resolve changes in the zeroth harmonics or, equivalently, the base flow. The resulting Galerkin system is derived, again, from the Navier-Stokes equation. The effect of mode interpolation on the model accuracy is quite dramatic. Both the initial exponential growth and the limit cycle are well resolved, i.e., we arrive at a uniformly accurate Galerkin model for the whole transient.
Insight for DMD related to change from linear to nonlinear dynamics is provided in [2]. We conjecture that the presented model order reduction with shift mode(s) for the base flow variation and mode interpolations for coherent structure changes is applicable for a large range of transient flow dynamicsin particular, oscillatory coherent structures. There is no algorithmic restriction. For multiple or broadband frequency dynamics, one should keep in mind that the interpolation between two states is far from trivial. An oscillation at one state may slowly deform into another oscillation during a transient, as in our study. This would require a single deformable mode basis. The transient may also display frequency cross talk, e.g., the mitigation of one frequency by the excitation of another [8,15,57]. In this case, the modes of both frequencies need to be included during the transient [25]. Otherwise, the continuous mode interpolation might lead to a sudden jump between both oscillations by mode switching in the corresponding interpolated eigenvalue problem. This would be a very coarse approximation of the transient behavior. Summarizing the general basis interpolation between two operating conditions with complex behavior remains a grand challenge problem. The presented model order reduction constitutes a simple method to start with.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Shift mode
Galerkin models based on empirical modes like POD have proven to be efficient in the resolution of the kinematics of the flow. Unfortunately, they are also sensitive to the changes in flow parameters, operating and boundary conditions. For example, the first two POD modes for the circular cylinder wake capture about 95 % of the fluctuation energy, but the model based on these modes alone is unstable [31]. Stabilization of the low-dimensional model might be done with the use of mean-field theory. Additional mode, representing slow, non-oscillatory changes in the flow field, is included in the Galerkin model in order to shift the center of the attractor of the (changing) system. For the transient wake dynamics, the shift mode is defined as a normalized difference between time-averaged flow computed for stable limit cycle u 0 and unstable steady solution u s [30]: As the steady solution of Navier-Stokes equations above critical Reynolds numbers might be difficult to obtain (although methods like selective frequency damping [1,6] are an option), another ways of the shiftmode computation might be considered: I One possible approach is the modal decomposition of the entire transient, for example, using POD. The most dominant of these modes is a candidate for a global shift mode, II Local mean-field correction might be computed using modal approximation, defined in (I), as the low-pass filtered time derivative of Fourier coefficients u t = i For more details on these methods, we refer the reader to [50].