Reduced order models for thermally coupled low Mach flows

In this paper we present a collection of techniques used to formulate a projection-based reduced order model (ROM) for zero Mach limit thermally coupled Navier–Stokes equations. The formulation derives from a standard proper orthogonal decomposition (POD) model reduction, and includes modifications to improve the drawbacks caused by the inherent non-linearity of the used Navier–Stokes equations: a hyper-ROM technique based on mesh coarsening; an implicit ROM subscales formulation based on a variational multi-scale (VMS) framework; and a Petrov–Galerkin projection necessary in the case of non-symmetric terms. At the end of the article, we test the proposed ROM formulation using 2D and 3D versions of the same example: a differentially heated cavity.


Introduction
The main purpose of this paper is to develop a model reduction formulation suitable for thermally coupled flows, by expanding the techniques in projection-based model reduction developed for several applications on fluid dynamics-mostly for incompressible Navier-Stokes equations-to the zero Mach limit Navier-Stokes equations developed in [1,2].
Following the model reduction developments in [3] we choose a POD model reduction approach. The POD, as any projection-based model reduction, aims to describe any phenomena that otherwise would be represented by a 'computationally expensive' numerical method with a surrogate lower dimensional model. This surrogate model is obtained by projecting the computational expensive numerical approximation onto a previously computed reduced space. Thus, the model reduction is arranged in two stages: an offline part, where the solution obtained from a 'high fidelity' full order model (FOM) is used to build the desired reduced order subspace; and an online part, where by projecting the original model onto the reduced subspace the ROM is built and subsequently solved. The traditional POD model reduction approach presents certain drawbacks when considering non-linear complex problems: added computational cost of representing the non-linearities over a linearized computational model, inherent numerical instabilities caused by the Navier-Stokes formulation, and non-optimal projection over the reduced space caused by the asymmetrical nature of the formulation.
To overcome the first issue, a wide variety of methods inspired by the work of Everson and Sirovich [4] have been introduced [5][6][7][8][9][10], with the term 'hyper-ROM' coined by Ryckelink [11]. This family of methods consists in using a sample of the geometrical domain. In this paper, we propose an idea of hyper-ROM that differs from the sampling way: to set the ROM on a geometrical space smaller than the original one. This, applied to mesh based methods, implies the interpolation of the developed ROM-basis includedonto a coarser mesh.
To stabilize both the FOM and ROM, we follow the same approach: a VMS framework. For the FOM case, we follow the stabilized formulation of the thermally coupled zero Mach limit Navier-Stokes equations developed in [12], using the Orthogonal Subgrid-Scales (OSS) as defined in [13]. For the ROM case, we propose an analogous method where the OSS are defined orthogonal to the ROM sub-space.
Lastly, to solve the non-optimality in the reduced space projection, we use a Petrov-Galerkin projection instead of the standard Galerkin projection as proposed in [7]. In the specific case of mesh based methods, the POD-ROM method can be seen as a solution of the original problem using a projection method as in [14], where the Petrov-Galerkin projection is the one that satisfies the non-singularity of the system of equations when the FOM is not symmetric.
This article is organized as follows. In the first section, we present a brief description of the zero Mach number limit Navier-Stokes model. In the second section, we describe briefly the idea of model reduction along with an explanation of the POD compression technique for the construction of the basis. In the third section, we present a modified finite element (FE)-ROM implementation, including the stabilization using VMS, the new proposed hyper-ROM method, and the Petrov-Galerkin projection. In the fourth section, we present examples consisting in a differential heated cavity, with 2D and 3D cases. Finally, we close the paper with some conclusions.

Thermally coupled flow problems
Let us start by describing a general thermally coupled Navier-Stokes formulation as a system of nonlinear convection-diffusion-reaction equations posed in a time interval (0, t f ) and in a domain ∈ R d with a boundary , as: with the nonlinear operator L defined as: and the boundary and initial conditions: where Y = (u, p, T ) is the vector of unknowns (velocity u, pressure p and temperature T ), n is the geometric unit outward normal vector on , F is a known vector of d + 2 components and M, A i , K ij and S are (d + 2 × d + 2) mass, convection, diffusion, and reaction matrices respectively. The usual summation convention is implied, with indices running from 1 to the number of space dimensions d, and the convection matrix i is split in order to define the appropriate Dirichlet ( D ) and Neumann ( N ) boundary conditions.

Low Mach number model
Following the description of the model in [1], the zero Mach limit Navier-Stokes model is based on the splitting of the pressure into two relevant parts: mechanical p, which behaves as in the incompressible and Boussinesq approximations; and thermodynamic pth, which is constant over the domain but dependent on the energy added in such domain and therefore time dependent. Using the previous notation for the thermally coupled flows, this formulation can be written as: where 0 is the vector of R d with zero in all its components, e i is the unit vector in the i-th cartesian direction, I is the d × d identity matrix, and ρ is the fluid density. Viscosity μ, thermal conductivity λ, specific heat at constant pressure c p , and heat capacity ratio γ are used as constant physical properties; and g and Q represent a body force vector and a given external heat source respectively. Considering that a state equation is necessary to close the system, for simplicity, in this work we will focus exclusively on ideal gases, with p = ρRT as the state equation and R the specific gas constant.
Finally, the formulation is complete with the thermodynamic pressure equations: for open flows ( u N = ∅), where the thermodynamic pressure is determined through global conservation equations over (Eq. 5a); and for closed flows ( u N = ∅), where the thermodynamic pressure is determined by the boundary conditions (Eq. 5b):

Model reduction
Projection-based ROMs rely on the existence of a reduced dimensional sub-space that approximates the solution space. Let us define a high dimensional space Y h of dimension M, with ϕ its orthonormal basis. Then, the i-th component of any element Y h ∈ Y h can be written as the linear combination Y h,i = M k=1 (Y h,i , ϕ k i )ϕ k i , with (:, :) a L 2 -inner product. Note that Y h may be an approximation to a continuous-infinite dimensional-space. Since for most cases the exact basis ϕ is unknown, we can define a lower-dimensional space Y r ⊂ Y h of dimension m, which approximates Y h as m → M, with a basis φ. Using this test basis, we can approximate the i-th component of any element where the accuracy of the approximation will depend on how accurate is the basis φ compared to the exact basis ϕ.

Construction of the basis
Following previous works in model reduction in FE [3,7,9], we use the POD statistical procedure as a way to build the reduced order sub-space basis φ. The objective of the POD method is finding a basis for a collection of high-fidelity 'snapshots' to use it as a the basis of the desired reduced sub-space.
Let us organize a set of data that has been previously obtained from a solution of the problem in Eq. 1, as a N -collection of 'snapshots' Y s h,j = Y h,j −Ȳ h in an ensemble: withȲ h the mean value of the snapshots. Now, let us denote as {φ k } m k=1 an orthonormal basis of Y h , where any member Y h,j of the ensemble in Eq. 6 can be approximated as: The POD consists in finding the orthonormal basis {φ k } m k=1 , such that, for every k ∈ {1, . . . , m} the mean square error between the elements Y h,j , 1 ≤ j ≤ N , and the corresponding j-th partial sum of Eq. 7 is minimized on average: To solve the problem in Eq. 8 we follow a singular value decomposition (SVD) of the 'snapshots' collection in Eq. 6, where the resulting left singular-vectors of the decomposition represent the basis {φ k } M k=1 . Here, based on the Eckart-Young-Mirsky theorem we can define a reduced basis simply by truncating the left singular-vectors at the m column as: {φ k } m M k=1 . As a criterion for the truncation, we use the retained energy η defined in [3] as: where λ k are the SVD non-zero singular values.
Remark Note that this is the way we approximate φ, but there are several other ways to find a basis for the reduced order sub-space (as the ones in [15,16]). The projectionbased model reduction formulation below should be valid for any basis regardless of the technique used to obtain it.
Remark In this paper, the 'snapshots' data is obtained using a FE approximation and therefore, the basis φ approximates the FE space Y h -of dimension M-which in turn approximates the space of the continuous problem.

FE-ROM formulation
Having stated the equations that represent the physical problem and the standard model reduction approach, we now describe the FE-ROM approximation of the problem. In this section, let us denote Y as the functional continuous space where Y exist. Instead of following a standard Galerkin approximation of the variational problem-where the FE space is denoted as Y h -we construct the approximation space Y r ⊂ Y for the ROM. The FE space is assumed to be built from a FE partition T h = K of the domain . The order of the FE interpolation is irrelevant for our discussion.
Introducing the notation (:, :) = (:, :) and (:, :) , for the L 2 -inner product on and respectively, or for the integral of the product of two functions-if they are not square integrable but their product can be integrated, we can define the variational problem as finding Y r ∈ Y r , such that: where ϒ r is the test function and the problem is written using the forms B and L defined as: Note that the terms where Y is not approximated by Y r in Eq. 11 involve non-linearities.
In a further section we discuss the linearization scheme and two ways of approximating these terms.

VMS for FE-ROM
Given the well-known lack of stability in the Galerkin standard formulation-present in convective dominated regimes-a stabilization technique is necessary. Inspired in previous works that acknowledge instability issues [17,18] and using a VMS framework as done for several problems in FE approximations, we develop what we call FE-ROM-Subgrid-Scales (SGS), which resembles FE-SGS. We start by splitting the unknown of the continuous problem Y in a part Y r , which can be resolved by the FE-ROM standard approximation; and a remainderY , the FE-ROM-SGS. In this way we can re-define the approximation of the continuous space as Y = Y r ⊕Y , where the SGS spaceY is any space that completes Y r in Y . Then, the variational form of the problem in Eq. 10 expands into two: finding Y r ∈ Y r andY ∈Y , such that: Since it is desirable to avoid derivatives over the subscaleY , we can re-write the bilinear form B(Y ;Y , ϒ r ) in Eq. 13a by integrating by parts within each element and assuming that the exact tractions are continuous across inter-element boundaries, yielding: where (:, :) K is the L 2 inner product over element K , and L * is the adjoint of the lineal operator L(Y ; :). Following the formulation in [12], this adjoint operator L * is:

Subscales approximation
At this point, we have described two sets of equations: a FE-ROM one, in which all the terms have been described; and a FE-ROM-SGS one. In order to solve the second system, we follow a standard VMS formulation. After doing integration by parts in some terms, neglecting the tractions across inter-element boundaries, and replacing u·∇ρ by − ρu T ·∇T ; we can re-write Eq. 13b as: with the residual defined as We can re-write Eq. 16 in the following way: where ϒ ⊥ r is a term that ensures that Eq. 18 belongs toY ⊥ and its definition depends on the choice ofY . Now, using the algebraic approximationL(Y ;Y ) ≈ τ −1 (Y )Y , we can re-write Eq. 18 as: where τ is the matrix of stabilization parameters that depends on K and the coefficients of the operatorL. It is important to notice that so far the matrix τ does not depend on the choice of the approximation space Y r , and therefore we can use the same definition of the stabilization parameters as in the FE approximation of the problem. In this case, we use the definition of the algebraic operator τ given by [2,13]: with the stabilization parameters defined as: where h is the element size and c 1 and c 2 are algorithmic constants.
To complete the approximation of the subscales, it still remains to define the term ϒ ⊥ r , that is, to choose the appropriate subscales spaceY . Based on the orthogonality property of the basis φ and on the work developed in [13,19], we define the subscales spaceY = Y ∩ Y ⊥ r , which can be thought as being orthogonal to the approximation space Y r , that implies having: To numerically compute ϒ ⊥ r we follow the supposition made in [20], where the orthogonality between Y r andY is defined respect to the weighted inner product (Y , ϒ) M = (Y , M(Y )ϒ). Following this condition any subscaleY must satisfy: Now replacing the orthogonality condition (Eq. 23) in Eq. 19, we obtain: from where it follows that ϒ ⊥ r is the following projection onto the space Y with respect to L 2 -inner product, denoted by r : Replacing ϒ ⊥ r into Eq. 19, and using the approximation [20], we find an expression for the FE-ROM-SGS: with ⊥ r . .= I − r defined as the orthogonal projection on Y r , where I is now the identity in Y r . Finally, by introducing the orthogonality condition (Eq. 23) in Eq. 14, neglecting the term M(Y )∂ t Y in the residual R of Eq. 26-due to its orthogonality toY given by M(Y )∂ t Y ,Υ = 0, and rearranging Eqs. 14 and 26, we can state the problem in Eq. 13 as finding Y r ∈ Y r andY ∈Y , such that: with the adjoint operator L * and the residual R redefined as: Remark By following the same analysis performed when deriving the orthogonal SGS in [13,19], we have come to a rather similar definition of the FE-ROM-SGS, where the most important difference lies in the definition of the orthogonal projection ⊥ in Eq. 26. In the FE case the projection is done onto the space Y ⊥ h , while in the FE-ROM approximation it is done onto the space Y ⊥ r .
Remark The choice of the stabilization parameters τ is done following the Fourier analysis done in [13]. Since the information represented by the reduced basis corresponds to the resolved scales from the FOM, the subscales for both the FOM and the ROM are part of the continuous solution which cannot be approximated by the FOM.
Remark The previous definition of the FE-ROM-SGS is equivalent to the dynamic orthogonal SGS model in [20]; it is important to acknowledge that subscales could also be implemented without the temporal term (quasi-static), or not orthogonal to Y r . An extensive analysis of the FE equivalent models is depicted in [12,20].

Time discretization
Any time integration scheme could now be applied to discretize in time the FE-ROM and FE-ROM-SGS equations together (Eq. 27). Considering the results in [21], where it is shown that the time integration for the subscales could be one order less accurate than for the finite element equations without affecting the accuracy of the scheme, we choose for this work two Backward Differentiation Formula (BDFs): of second order for the FE-ROM equations and of first order for the FE-ROM-SGS ones. Setting a uniform partition of the time interval of analysis [0, t f ],with δt the time step and n the superscript that denotes the current time step, we can approximate the temporal derivatives in Eq. 27 using: Replacing the time integration scheme (Eq. 28b) in Eq. 27b we get an equation for the subscales: with the matrix of effective stabilization parameters is defined as τ t (Y n+1 ) . Now, replacing Eq. 29 and the integration scheme (Eq. 28a) in Eq. 27a, we get:

Linearization
To solve the non-linearity present in the terms involving Y , we implement a linearization scheme based on Picard's method. Using the terminology used in [20]

Discrete approximation
We can describe the discrete representation of the FE-ROM problem as a composition of the FE and ROM approximations. In FE the space Y h is defined as made of continuous piece-wise polynomial functions in the domain , where we can write the discrete approximation of the unknown as with N (x i ) the shape function of node i. In contrast, in ROM we approximate the unknown Y as Y (t) ≈Ȳ + m k=1 φ k Y k (t). Using these two approximations, we can describe the space Y r in two ways: as a FE space Y h , represented using orthogonal basis φ, or as a ROM approximation of the problem, discretized in using continuous piece-wise polynomial functions. In that way, we can write the discrete representation of Y r as: (31)

Petrov-Galerkin projection
Since the formulation presented above introduces the non-symmetrical linear operators L and L * , it is necessary to find an optimal projection in a way that it gives us a feasible solution [14]. To solve this lack of optimality in the projection, we replace the traditional Galerkin projection by the Petrov-Galerkin projection defined in [7]. Let us re-write the linearized version of Eq. 30 as the linear system: where is the discrete basis matrix, L and R are the discrete left and right hand sides of Eq. 30, and Y r the discrete unknown.
To apply in a natural way the Petrov-Galerkin projection we define W = L as a matrix whose column vectors form a basis to the projection space that allows us to transform the left hand side of Eq. 32 into a positive semi-define matrix. Projecting Eq. 32 onto such space we get: Using the orthogonality property of the basis = I, Eq. 33 becomes one that resembles the Petrov-Galerkin ROM formulations in [7,9]:

Hyper-ROM
Lastly, in order to reduce the computational cost of evaluating nonlinear terms, we propose a mesh-based hyper-ROM as an alternative to the sampling-based domain reduction algorithms [5][6][7][8][9][10].
The mesh-based hyper-ROM consists in the solution of the described ROM problem using a coarser mesh than the one of the FOM. The implementation of this technique is done straightforwardly by writing the discrete approximation (Eq. 31) in function of the new coarser mesh.
Ideally, the coarsening should be performed as a function of the 'less important' areas of the geometry, which can be achieved using already existing mesh refinement algorithms. In the subsequent examples we test this technique using a uniform coarsening of the mesh.
Remark When the POD basis is obtained by sampling a mesh-based solution-a FE one for example-the coarsening of the mesh implies an interpolation of such basis.

Numerical examples
In this section we present two examples consisting of 2D and 3D versions of the initial transient part of a differentially heated cavity of aspect ratio 1-similar to the one presented in [12,22] (T h +T c ) = 3.55 · 10 6 . For the FOM solution (as a reference solution and for the construction of the snapshots), we follow the VMS formulation presented in [20], using the dynamic orthogonal SGS model. In all cases, we use a constant time step size of δt = 0.01 s. For the finite element meshes used (described below) this time step is slightly higher than the critical time step of an explicit scheme due to advection, which in turn is of the same order as that due to viscous effects. Since we use an implicit scheme, we are not restricted to a critical time step size, but this observation serves to justify that our choice is adequate.
Additionally to the following examples, we have tested a differentially heated cavity of aspect ratio 8 analyzed in [23], getting inconclusive results.

Two dimensional case
In the 2D problem, we use 2 uniform structured meshes composed of quadrilateral elements: one with 10,000 elements and a mesh size h = 0.01, used for the solving the FOM and the ROM; and one with 2500 elements and a mesh size h = 0.02, for testing the hyper-ROM formulation. To construct the basis, we collect 500 snapshots for velocity, pressure and temperature at every 4 time steps in a 20 s interval. Figures 1, 2, and 3, show temperature contours of the solutions obtained using the FOM, ROM, and hyper-ROM formulations respectively.
To evaluate the accuracy of the ROM and hyper-ROM formulations, we perform a series of numerical tests using 12 different sets of basis, varying the retained energy from η 1 = 0.99 to η 12 = 0.3. Figure 4 depicts the number of vectors in the basis in function of the retained energy.
As expected, we observe a more diffusive behaviour, in both the mean and the fluctuation of the Nusselt number, when fewer basis vectors are included. The hyper-ROM results appear to have the same behaviour of the ROM with lower amplitude. Additionally, we and the root squared mean of i -for N time steps-is defined as: (37) Figure 7 shows the convergence of rms over the hot and cold walls, and the mean value over both of them¯ rms . We also include the convergence error for the FOM using the coarse mesh.
Although the convergence error does not have a clear slope, it behaves as expected, with the error decreasing with the addition of basis vectors. It is important to notice that appears to be an optimal value for η = 1-the maximum number of basis vectors-where the error reaches the minimum; considering that this results occur near η = 1, we believe it can attributed to an overfitting phenomenon. The same idea is explored in [24,25], where it is attributed to a lack of smoothness or noisiness in low energy basis vectors. The overall results for the hyper-ROM seem adequate, given that the error is the same order of the ROM and the coarse mesh FOM solutions.
To evaluate the temporal evolution we perform a discrete Fourier transform of the Nusselt number. In Figs. 8 and 9 we compare the Nusselt number spectra for the hot and cold walls respectively.
Again we see how reducing the amount of basis vectors leads to a more diffusive solution. But in contrast to Fig. 7, in Nusselt number spectra we observe that the ROM and hyper-ROM spectra tend to the FOM spectra as we approach η = 1.
Regarding computational time, Fig. 10 shows speedup for all the cases under the same initial conditions, time step, and computer configuration. As expected the computational time is reduced when the coarser mesh is used. It is important to note that the computational time for the ROM and hyper-ROM is decreased not only by the solution of a smaller linear system but as well for having fewer non-linear iterations.

Three dimensional case
For the 3D problem, we use 2 uniform structured meshes composed of regular hexahedral elements: one with 64,000 elements and a mesh size h = 0.025, used for the solving the    To test the 3D problem we perform similar tests as the ones in the 2D problem using 6 different sets of bases, varying the retained energy from η 1 = 0.99 to η 6 = 0.5 Fig. 12 depicts the number of vectors in the basis as a function of the retained energy.
As in the 2D problem, for the 3D case we compute the Nusselt number, the root mean square of the Nusselt error and the discrete Fourier transform of the Nusselt number for the hot and cold walls, getting similar results.
In Figs. 13 and 14 we compare the Nusselt number spectra for the hot and cold walls. Although the results behave in a similar way as in the 2D case, the ROM and hyper-ROM  seems to be more dissipative, in spite of the ratio of the meshes being smaller. We believe this could be attributed to the quality of basis, which directly depends on the quality of the snapshots used to construct it. Figure 15 shows the root mean square of the Nusselt error for ROM and hyper-ROM, showing similar results as the 2D case, including the overfitting phenomena when the retained energy η approaches 1. Figure 16 shows speedup for ROM and hyper-ROM solutions. The reduction in the computational time for the 3D problem is larger than the one achieved for the 2D problem, using a similar amount of basis vectors.

Conclusions
In this work we have developed a formulation that allowed us to perform ROMs on the zero Mach number Navier-Stokes approximation. For that purpose we have described a set of tools that allow us to tackle the main problems that arise: lack of stability, added computational cost due to the non-linearities and non-optimal projection over the reduced space.
To solve the lack of stability, we have proposed a VMS-FE-ROM formulation with the following characteristics: • It is an implicit formulation, reducing the offline costs and easing the implementation. • Since it is built over a VMS framework, it incorporates non-linear terms in the adjoint operator (Eq. 15), which may become relevant in complex flows. • It maintains the orthogonality definition, proposed in [18,26], between the sub-spaces Y r andY . • It is residual-based, as shown in Eq. 18. • It works with the same stabilization coefficients τ as the ones developed for the VMS-FE formulation in [2,12].
Although in this work we focused solely in the use of dynamic subscales, it is possible to construct the same method for quasi-static subscales or other stabilization techniques.
To solve the added computational cost, we have proposed a mesh-based hyper-ROM that contrary to the traditional sample-based methods, does not require any algorithm to select these sampling points. The numerical experiments performed show that this hyper-ROM method behaves appropriately, only inducing an expected added diffusion to the solution. A natural extension of the mesh-based hyper-ROM method is the use of mesh refinement techniques to improve both the computational time and the accuracy of the method.
Finally, to tackle the non-optimal projection of the non-linear formulation, we give a different interpretation to the Petrov-Galerkin projection-originally described in [7]-in order to include it in our formulation.