Model Order Reduction of Turbulent Mixed Heat Transfer Problems by Projection Method

Computational fluid dynamics modeling has been widely applied to study flow and heat transfer problems in nuclear engineering. However, due to the prohibitively high computation cost, fast or real-time solution of nontrivial problems is still a challenging task under context of many-quires. Model order reduction is an efficient way to achieve significant speedup. For the turbulent mixed heat transfer problems, the Proper Orthogonal Decomposition (POD)-Galerkin projection method is introduced and then applied to solve the problem of non-isothermal mixing in a T-junction pipe with the inlet velocities treated as parameters. It is shown that a speedup of two orders of magnitude has been achieved. For new parameters out of the scope of the samples, the reduced order model still can give satisfactory results.


INTRODUCTION
With tremendous advances in computational power, computational fluid dynamics (CFD) has gained wide applications in the study of flow and heat transfer problems in nuclear engineering. However, CFD simulations often face great computational challenges considering finely resolved grids have to be used for the nontrivial transient applications. As early as 2012, the CASL (Consortium for the Advanced Simulation of Light Water Reactors) project attempted to couple the commercial CFD software, i.e., STAR-CCM+, and the neutronics software to simulate the flow and the heat transfer behavior of a 1/4 reactor core [1]. Due to the large scale of the problem, the computation took 10,000 CPU cores and 1 TB of memory. A recent study applied the CFD simulation to study the thermal shock of a reactor pressure vessel [2]. The simulation of a 300s transient process took about 200 hours (Hardware: 64 cores, 64 GB memory; Software: STAR-CCM+). Obviously, it is impractical to conduct parametric studies or do optimal design by using the costly CFD simulations. Despite the great advances in computation power in recent years, it is still difficult to achieve a speedup of several orders of magnitude just through parallel computation considering the limitation imposed by the communication overhead between nodes [3].
Conducting parametric studies and investigating major effects are the main purpose of numerical modelling. When conducting parametric studies, the underlying governing equation of the problem remains the same but the physical parameters (e.g., viscosity coefficient), boundary conditions (e.g., inlet velocity and temperature) and even the spatial domain (e.g., the diameter of the inlet) may change for each simulation. If the coefficients or the boundary conditions are considered as parameters, such problems can be described by a parameterized governing equation.
The most straightforward means of speeding up the solution of parametric problems is to apply model order reduction (MOR) methods to reduce the number of unknowns to be solved [4]. The basic idea of model order reduction is to project a full order model (FOM) containing large number of degrees of freedom (DOFs) into a low-order subspace, transforming a complex full order model into a reduced order model (ROM) containing only much fewer degrees of freedom. The most widely used model order reduction method for linear problems is the reduced basis finite element method [5]. However, when the nonlinear model described by the Navier-Stokes (N-S) equation is treated by the reduced basis finite element method, the inner product of the transport term and the reduced basis function leads to a trilinear term, which poses great difficulties in deriving the reduced order model [6]. Unlike the finite element method, the finite volume method integrates the governing equation directly within the control volume without defining the test function and the local stabilization term. Therefore, the model order reduction methods based on the finite volume method are more suitable for the flow problems.
Proper orthogonal decomposition (POD) and greedy algorithm (GA) are the two most commonly used methods for constructing the low-order subspaces. The former directly analyzes the typical solution snapshots while the latter requires a posteriori error estimator. The basic principle of the POD method is to use singular value decomposition (SVD) to filter out the most representative modes from the solution snapshots collected at different moments and parameters. Then an orthogonal basis = [ 1 , 2 , … , ] is constructed by using, for example, the Gram-Schmidt process [7]. The full order solution can thus be represented by a linear combination of the small set of basis functions as shown in Eq. (1), The freedom of the model is greatly reduced by using the basis function with global support.
The model order reduction methods built upon the POD can be divided into the nonintrusive methods and the intrusive methods. The non-intrusive methods do not depend on the governing equation and essentially build a surrogate model between the input data and the output data. Within the context of POD, the surrogate model is usually used to obtain the decomposition coefficient (i.e. in Eq. (1)) for the solution of a new parametric configuration. Radial basis function (RBF) and Kriging are two popular methods for construction of surrogate models [8,9]. The non-intrusive reduced order model, although simple to construct and fast to compute, is essentially a fitting of the decomposition coefficients and may produce large errors when solving transient problems. On the other hand, the intrusive method still requires solving the governing equations commonly projected onto the subspace spanned by the POD basis [10]. The intrusive reduced order model retains the physical properties of the original problem. So its mathematical form is more rigorous and the accuracy is usually higher even the new parameters are out of the scope for sampling the solution snapshots. Use of the POD-Galerkin projection method to construct reduced order models for flow problems has received much attention in recent years [11].
Compared to the applications in engineering fields such as the mechanical engineering and the aerospace engineering, application of the model order reduction method in nuclear engineering is rather limited. This paper investigates the model order reduction of a typical problem in nuclear engineering, i.e., turbulent mixed heat transfer, within the framework of finite volume method. The Reynolds averaged Navier-Stokes (RANS) equations are projected on the subspace spanned by the POD basis with the inlet conditions treated as parameters. Finally, the accuracy and the efficiency of the reduced order model of a T-Junction of pipes are studied. The results show a speedup of two orders of magnitude is achieved with acceptable accuracy. Finally, the predictive capability of the reduced order model is evaluated and discussed.

Governing equation
For the problem of turbulent mixed heat transfer of the incompressible fluid (only passive heat transfer is considered in this paper), the governing equations are given in Eq. (2) -(4), where is the time-averaged velocity field, is the temperature, � is the Renault mean strain rate tensor, is the density, and is the unit matrix. is the effective turbulent viscosity defined as = t + 1 of which is the turbulent viscosity and 1 is the molecular turbulent viscosity. * is the pressure field modified by the RANS model. is the effective turbulent thermal diffusion coefficient defined as = + 1 of which is the turbulent thermal diffusion coefficient and 1 is the molecular turbulent thermal diffusion coefficient. The two coefficients are respectively defined as = and 1 = 1 where is the turbulent Prandtl number and is the molecular turbulent Prandtl number. In addition, the boundary conditions and the initial conditions are given in Eq. (5), where = ∪ ∪ 0 , and , and 0 are the inlet, the outlet and the solid wall boundaries, respectively. The function f(x) represents the inhomogeneous boundary condition of the velocity field, g(x) represents the initial condition of the velocity field at moment 0, h(x) represents an inhomogeneous boundary condition of the temperature field, and k(x) represents the initial condition of the temperature field at moment 0.

Introduction of POD
The POD method was first applied to study the pseudo-sequential structure in turbulent flows [12]. After decades of development, the POD method has become a common method for building reduced order model and has found many applications in hydromechanics, heat transfer and other fields [13][14][15].
Before applying the POD technique, it is necessary to prepare sufficient sample data, i.e. the flow and the temperature fields of different parametric configurations at different moments. These field data are called snapshots and can be obtained from full order CFD simulations or experimental measurements. The matrix of snapshots consisting of a series of snapshots is represented as U ∈ ℛ × , where is the number of snapshots (i.e., one snapshot is stored as a column vector) and n is the length of each snapshot. U can be decomposed as shown in Eq. (6), where the columns of the matrices X ∈ ℛ × and Y ∈ ℛ × are the left and the right singular vectors of U, respectively, and Σ∈R n s ×n s =diag(σ 1 ,σ 2 ,…,σ ns ) is the singular values of U where 1 ≥ 2 ≥ ⋯ ≥ .The POD modes = [ 1 , 2 , … , ] are defined as N leftmost singular vectors of U, which correspond to the N largest singular values, respectively. POD provides an efficient low-dimensional representation of the snapshot data and the PODbasis is optimal in an ℓ 2 sense over the parameter space. It can be rigorously demonstrated that the square of the reconstruction error is equal to the first discarded singular value [7]. Therefore, the singular values are the main indicator for selection of the POD modes.

POD-Galerkin projection method
After obtaining the dominant modes of the system, the velocity field, the pressure field, and the temperature field can be expressed as a linear combination of these POD modes as shown in Eq. (7) -(9), ( ), ( ) and ( ) are the POD modes obtained by decomposing the snapshot matrix. α(µ,t), b(µ,t) and c(µ,t) are the decomposition coefficients of the velocity field, the pressure field and the temperature field, respectively. , and are the number of the POD modes of the three fields. µ is the vector of the parameters of the system (e.g., inlet velocity, inlet temperature, etc.).
The momentum equations are then projected onto the subspace spanned by the velocity modes, the continuous equation onto the subspace spanned by the pressure modes and the scalar transport equation onto the subspace spanned by the temperature modes as shown in Eq. (10) -(12), The RANS equations (i.e. Eq. (2)-(4)) can thus be converted into a set of ordinary differential equations with coefficients to be determined as shown in Eq. (13) -(23), . where Considering the above reduced matrix are precomputed low-order matrices, the degrees of freedom of the reduced order model (i.e., Eq. (13)-(15)) are greatly reduced.
For the RANS model, the turbulent viscosity can be decomposed into a linear combination of its POD modes in the same way as shown in Eq. (24), where ξ(x) and l(µ,t) denote the POD modes of turbulent viscosity and the corresponding coefficient. Considering that multiple forms of turbulence models (e.g., Spalart-Allmaras, kε, k-ω, etc.) [16] are available and that turbulence mainly introduces a new term in the momentum equation that approximates the vortex viscosity, it is enough to reduce the turbulent viscosity field only. To achieve a separation of reduced turbulent viscosity field from the specific turbulence model, the interpolation by using radial basis functions (RBFs) instead of the POD-Galerkin projection is applied in the present study. The basic form of the RBF interpolation is given in Eq. (25), where is the coefficient in Eq. (24), is the parameter vector and is the m-th sample parameter, ( ) is the kernel function and is the weight. For each sample parameter , the corresponding coefficient ( ) is known. So, the weight can be determined by solving a simple linear set of equations. For each new parameter vector (e.g., different moments or different boundary conditions), the POD decomposition coefficient for turbulent viscosity can be quickly obtained by Eq. (25). The form of the kernel function selected in this paper is given in Eq. (26), where the parameter 0 controls the region of influence of the kernel function.

Parameterization of boundary conditions
The main purpose of model order reduction is not to reconstruct known solutions, but to quickly solve solutions of new configurations. The homogeneous boundary conditions are included in the governing equations after projection and therefore do not require special treatment. For models that contain parametrized inhomogeneous boundary conditions, a lifting function can be used to homogenize the solution so that the boundary conditions can be separated from the snapshot matrix. For example, the homogeneous velocity field can be expressed as Eq. (27), where is the number of boundary conditions, is the velocity on the boundary, is the lift function at each boundary which takes the value of 1 at the boundary and 0 at other locations. The procedure is shown in the pseudo-code of Algorithm 1.
The POD decomposition of the homogeneous snapshot matrix yields the dominant modes of each field. And then the homogeneous solution of a new parameter can be obtained by solving the reduced order model (i.e., Eqs. (13)- (15)). Finally, the inhomogeneous solution can be recovered by imposing the boundary conditions through the lifting function. For example, for the velocity field as shown in Eq. (28),

TEST AND RESULTS
The above method is verified using the non-isothermal mixing in a T-junction pipe with the inlet velocities treated as parameters (Fig. 1). The diameter of the main pipe is =160mm and the diameter of the branch pipe is =40mm. The total length of the main pipe is =1.2m. The bending radius is =160mm. The branch pipe is located at 0.5 . Γ 1 , Γ 2 , Γ , Γ 0 are respectively the boundaries of the main pipe inlet, the branch pipe inlet, the main pipe outlet and the solid wall. The boundary conditions are listed in Table 1. In order to test the predictive capability of the reduced order model, the snapshot matrix is prepared by collecting the solution snapshots only with the boundary conditions shown Table 1. The reduced order model are built upon the POD modes and then used to obtain the solutions at other inlet velocities.
The full order model is solved using the open source CFD software OpenFOAM and the coupling of the flow and the passive heat transfer models (i.e., Eq (4)) is achieved by modifying the pisoFoam solver. The convection term is discretized by the second-order upwind and central difference scheme.
The diffusion term is discretized by the second-order central difference scheme. And the temporal term is discretized by the first-order backward Euler scheme. The two-equation k-ω model is adopted for the turbulence term. The total simulation time is 50s and the time step is 0.05s. The relative error of the reduced order model with respect to the full order model is defined as shown in Eq. (29), It is found that that an interval of 0.1s is accurate enough for collecting the solution snapshots. Figure 2 shows the cumulative energy errors of the velocity field, the temperature field, the pressure field and the turbulent viscosity field after POD decomposition. The required numbers of modes for the velocity field, the temperature field, the pressure field, and the turbulent viscosity field are respectively 20, 10, 5 and 25 if the truncation error is set to be 10 -4 . The flow and the temperature fields with other different boundary conditions are solved using the reduced order model. The maximum error of the velocity field with is shown in Figure 3. For the boundary conditions with an inlet velocity of 0.8 m/s and a branch pipe inlet velocity of 1.6 m/s, the solutions obtained by the reduced order model are compared with those obtained by the full order model (Table 2). It can be seen that the reduced order model constructed using the POD-Galerkin method shows good stability and satisfactory predictive capability.
In these calculations, the full order model takes about 232.6 seconds (single-core CPU, 32G memory) while the reduced order model takes only 3.3 seconds. A speedup of 2 orders of magnitude has been achieved.
In this paper, the solution snapshots are collected from only one set of parameters as shown in Table 1. It is anticipated when the new parameters are closer to the sampled parameters, the built reduced order model can yield accurate results (Fig.3). However, the accuracy of the reduced order model decreases when the new parameters are far from the sampled parameters. To further improve the accuracy of the reduced order model, it is a standard practice to collect the solution snapshots at multiple sets of parameters.

CONCLUSIONS
The POD-Galerkin projection method is an accurate and efficient method for solving parametric problems under many-queries context. In this paper, a model order reduction technique for a typical problem in nuclear engineering, i.e., the mixed turbulent heat transfer, is investigated. It is shown that a speedup of two orders of magnitude has been achieved.
The classical POD method usually requires a large snapshot matrix and multiple full order solutions have to be prepared. For large-scale industrial problems, however, every solution is very expensive. Further work may combine the greedy algorithm and the surrogate error model to develop an efficient sampling method.

DISCLOSURE STATEMENT
No potential conflict of interest was reported by the authors.
ACKNOWLEDGEMENT This work was supported by the Guangdong Provincial Science and Technology Project (Industry-University-Research Collaboration), grant number 2015B090901051.