The Application of Discontinuous Galerkin Methods in Conjugate Heat Transfer Simulations of Gas Turbines

The performance of modern heavy-duty gas turbines is greatly determined by the accurate numerical predictions of thermal loading on the hot-end components. The purpose of this paper is: (1) to present an approach applying a novel numerical technique—the discontinuous Galerkin (DG) method—to conjugate heat transfer (CHT) simulations, develop the engineering-oriented numerical platform, and validate the feasibility of the methodology and tool preliminarily; and (2) to utilize the constructed platform to investigate the aerothermodynamic features of a typical transonic turbine vane with convection cooling. Fluid dynamic and solid heat conductive equations are discretized into explicit DG formulations. A centroid-expanded Taylor basis is adopted for various types of elements. The Bassi-Rebay method is used in the computation of gradients. A coupled strategy based on a data exchange process via numerical flux on interface quadrature points is simply devised. Additionally, various turbulence Reynolds-Averaged-Navier-Stokes (RANS) models and the local-variable-based transition model γ-Reθ are assimilated into the integral framework, combining sophisticated modelling with the innovative algorithm. Numerical tests exhibit good consistency between computational and analytical or experimental results, demonstrating that the presented approach and tool can handle well general CHT simulations. Application and analysis in the turbine vane, focusing on features around where OPEN ACCESS Energies 2014, 7 7858 there in cluster exist shock, separation and transition, illustrate the effects of Bradshaw’s shear stress limitation and separation-induced-transition modelling. The general overestimation of heat transfer intensity behind shock is conjectured to be associated with compressibility effects on transition modeling. This work presents an unconventional formulation in CHT problems and achieves its engineering applications in gas turbines.


Introduction
The continuous increase of the combustor outlet temperature has been producing incessant breakthroughs in the performance of modern aero engines and heavy-duty gas turbines.Accordingly, the thermal load of hot-end components has been steadily increasing, which necessitates the application of more complex cooling structures, such as various types of film cooling holes, impingement cooling holes, fins and serpentine-shaped channels, etc.The use of these structures increases the conjugate effects of heat transfer between hot flow, cooling flow and solid.It challenges engineers to precisely evaluate the lifetime of components or even the safety of whole engines, so the development of accurate and efficient conjugate heat transfer (CHT) simulation methods has become an important research area.In 1999, Martin et al. [1] developed a sequential conjugate strategy, by which a code for gas flow field computation and another one for solving heat conduction in solid blade was combined through an inner iteration process.This method taking advantage of available single-filed solvers had then been promptly adopted by some investigations [2][3][4][5][6], most of which used mature parabolic or elliptic equation solvers based on finite element (FE) methods in the solid domain.Meanwhile, partially aiming at enhancing the stability of the algorithm, the direct or full conjugate strategy [7][8][9][10][11] that integrates the governing equations in both fluid and solid domains using a single code also began to be developed and applied.
On or around the interface of two adjacent domains, the numerical accuracy together with the data exchange process has great effects on the computational performance in CHT problems.For improvement in this respect, more accurate algorithms can be enforced in both domains.Furthermore, more direct or localized process of data exchange is also rewarding.On these grounds, the discontinuous Galerkin (DG) methods can be considered to be applied in CHT simulations.
The DG methods that were first put forward by Reed and Hill [12] in 1973 and then further generalized by Cockburn and Shu [13][14][15][16], and Bassi and Rebay [17,18] have been successfully used in the simulations of inviscid, laminar and turbulent flows, to this day.The methods adopt the cell-wise approximate solutions in polynomial space based on the variational principle to obtain arbitrary high-order accuracy, whilst the eigen-decomposition (Riemann Solver usually) is used on cell-boundaries to account for the physics of wave propagation, which ensures good resolution of discontinuities.Comparing with the traditional finite volume (FV) framework for space discretization in computational fluid dynamics (CFD), the DG methods achieve high-order accuracy without the high-order interpolation process that comes at the expense of stencil expansions.As a consequence, it has better compactness and its solution accuracy has a lower dependence on the geometry and mesh properties.For the above reasons, together with the fact that it has been showing superior performance in CFD during last two decades, the introduction of DG methods to the whole domain of CHT problems can be beneficial for improving the accuracy in the region around the interface and localizing the data exchange process, and further enhancing the convergency and stability of the entire computation to some degree, especially in the situation of an interface with complex geometry or meshes of poor quality.Besides, one of the potential advantages of the idea is that the DG methods present a more natural way than the classic FE methods to represent the thermal discontinuities inside the solid, caused by thermal contact resistance [19], and so on.
In this paper, we present a primary framework for CHT simulations applying the DG methods on unstructured grids and a corresponding code fit for general engineering problems in gas turbines is developed.To this end, the compressible fluid dynamic equations and solid heat conduction equations are discretized in space and time under the unified algorithm.Taylor basis expandion at the centroid of each element, regardless of the element geometry type, is adopted.The Bassi-Rebay method [17], the ideology of which is borrowed from the mixed finite element method, is used to compute the gradients, including the temperature gradients of solids.A fully coupled strategy based on the data exchange process via the numerical flux of interface quadrature points is adopted.Moreover, considering the great effects of turbulence and transition behavior of boundary layer on the heat transfer, the code also assimilates various turbulence models and a local-variable-based two-equation transition model, γ-Reθ developed by Menter and Langtry [20,21] in recent years.All the modeling equations are solved using the DG methods as the basis equations, consequently forming a coupling equation system under a unified algorithm.Some test cases are presented to validate the feasibility of the method.Then we apply it in the simulation and analysis of a convention cooled transonic turbine vane.

Fluid Domain
Here is a description of the DG methods we used in viscous fluid dynamics.The compressible Navier-Stokes or Reynolds-averaged Navier-Stokes (RANS) equations in the fluid domain ΩF are given by: F , in where: (2) The total stress tensor T is calculated by viscous constitutive relation (in which the molecular viscosity coefficient obeys the Sutherlands formula for air) and Boussinesq assumption, and the heat flux q is calculated by Fourier's law (in which the thermal conductivity is λ = cp(μ/Pr + μt/Prt).The turbulence and transition variables are calculated through the modeling equations which will be presented below.By means of the inner product of Equation ( 1) and a test function v in ΩF and the Green's formula, we obtain the semi-discrete weak formulation: where   where: we can translate Equation (3) into the form: Uu xx (6) For the calculation of gradients, the auxiliary variants Pj such that: F 0, in are introduced.The weak form can be derived from the inner product of Equation ( 7) and a test function v in ΩF: Assuming the numerical solution eq h, h j  P , we obtain the gradients calculation: where , Pp xx (10) In this paper, the All-Speed-Roe type scheme devised by Li and Gu [22] is selected as the inviscid numerical flux ˆj F in Equation ( 5) to guarantee its reasonable numerical dissipation in the situation of low Mach number, and consequently contributes improvements to the accuracy and stability in the boundary layer regions.The numerical flux ˆj D and Û use the central scheme following [17].

Solid Domain
We use an analogous numerical method for the solid domain with the fluid one.The Fourier Equation in the solid domain ΩS are given by: S p,S S ρ 0, in where: Analogizing with the DG manner Equations ( 3)-( 6), we can derive the formulation: where Equation ( 12) can be also transformed into the equation below as Equation ( 8) into Equation ( 10): where , hj E q such that: Again the central scheme is used for the numerical flux ˆj q and T .

Fluid-Solid Interface
The continuity in physics of the temperature and normal heat flux on the fluid-solid interface should be ensured in any computational process.For this purpose, the precept is adopted that in each time step, the value of temperature computed from the side of solid domain is specified to the side of fluid meanwhile that of normal heat flux computed from the side of fluid domain is specified on the solid side, and the updated data will be treated as new boundary conditions of the respective domains in the next time step.In the DG framework, the solution distribution in each cell is naturally known at a certain time step.If the interface mesh on one side is consistent with that on the other (the situation in which the cases of this paper are), the above data exchange process can be directly achieved via the numerical flux calculations at quadrature points of the interface, without any interpolation.For any monotonicity preserving flux scheme represented by local variables on two sides of cell-boundaries, the numerical flux on the fluid side of the interface can be written as: and the flux on the solid side can be: The formulations of Equations ( 17) and ( 18) can implicitly guarantee the thermal coupling workable.

The Taylor Basis
In the framework, we choose the Taylor series expansion at the centroid of each element as the basis functions to represent the numerical polynomial solutions.The basis can be constructed naturally while generating a mass matrix with relatively small stiffness, and also enable the formulation adaptive to any types of element geometry.Considering the two-dimensional situation, for example, and aiming at rewriting Equation ( 6) (similarly for the Equations ( 10), ( 14) and ( 16)) as the series formula:  (19) where the subscript c represents the values at the centroid of E k , we can derive a set of appropriate basis functions as: and: in which k E is the volume of the element.The element scales, ∆x and ∆y can be chosen as:

Numerical Treatments for Turbulence and Transition Models
Various turbulence models, including k-ω [23], Baseline k-ω (BSL) and Shear Stress Transport (SST) [24], are tested in our framework.For these models, the conservative variables and source terms are: The logarithm of turbulence frequency ῶ = lnω is chosen to be in the conservative variable instead of ω, because of the smoother distribution pattern of the former in the near wall region.While the expressions of the production terms Pk, Pῶ, dissipation terms Dk, Dῶ and cross diffusion term CDῶ can be trivially derived from those in relevant references, the extra term Yῶ coming from the diffusion term in the original transport equation of ω has the formula: To decrease the adverse impact rooted in the stiffness of turbulence source terms on the computations to some extent, we briefly derive an approximate implicit treatment of the source terms to make the models more compatible with the entire explicit code.Following the approximate linearization advised by [24]: and ignoring Yῶ, we can deduce the approximation: Thus, the source-Jacobian matrix can be written as: The transition model used in the paper, γ-Reθ, consists the transport equations of turbulence intermittency factor and transition onset Reynolds number based on the momentum thickness of boundary layer.The conservative variables are: The model is strictly based on local variables, so it has good compatibility with CFD methods on unstructured grids and makes it possible to predict transitions in complex three-dimensional flow fields.More detailed descriptions of the model can be found in the original literatures [20,21], and the modification of separation-induced transition in [21] is also used in this paper because of the shock-induced separation in the case of turbine.

Method Verification
For the test cases in the present work, a vertex-based hierarchical slope limiter devised by Kuzmin [25] is used for fluid elements to restrain the pseudo-oscillation near discontinuities and consequently improve the numerical stability.We take second-order accuracy in space, i.e., first-order polynomial base functions for all the numerical examples below.Explicit three-stage Runge-Kutta scheme is used for the time integration of Equations ( 5) and ( 13) synchronously to achieve the steady solution ultimately.Unified domain decomposition process of the combined domain, including the fluid and solid, is employed through the software METIS and then the parallel computation can be implemented in each part between which communications are based on the Message Passing Interface (MPI) technology.

Heat Convection along a Flat Plate
The first case to validate the method presented in this paper is a 2-D simulation of the heat convection of laminar air flow with a total temperature of Ttin = 328.97K and total pressure of ptin = 107.853kPa, along a flat plate with the length of L = 66 mm and thickness of d = 3 mm, as showed in Figure 1.The isothermal condition Tb = 288.15K is enforced on the bottom surface.The Reynolds and Mach numbers are ReL = 4.4 × 10 5 and Ma = 0.30.The mesh is devised to ensure more than 20 grid lines located inside the boundary layer region.Various kinds of solid materials are selected such that the thermal conductivity ratios of solid to fluid λs/λF range from 1 × 10 1 to 1 × 10 4 .Figure 2 illustrates that under the method in this paper, the residual of temperature in the fluid domain can decrease nearly five orders of magnitude and that in the solid it can decrease over seven in this case.The chokes of decreasing residuals are conjectured to mainly originate from the lack of p-multigrid strategy in the available program, which may typically improve the convergence performance of DG formulations with natural or Taylor basis, although the residual minor oscillations in the results do not affect their reliability.Figure 3 shows that the temperature inside the solid has a tendency to the homogeneity with the increase of λf/λs, and the heat transfer surface also tends to the isothermal state.Figure 4 shows that when λf/λs ≈ 9 × 10 3 , the computed distribution of Nusselt number Nu = hL/λF along the surface is identical with the analytical correlation derived by Eckert and Weise [26].

Heat Transfer of a Poiseuille Pipe Flow
The second case is the conjugate heat transfer of a Poiseuille pipe flow, of which the analytical solution was derived by Montenay et al. [2].The geometry of the case is sketched in Figure 5. Air flows into a pipe with the length of L = 3 mm, inner diameter of 2a = 0.2 mm and outer diameter of 4a = 0.4 mm.The outer surface and both ends of the pipe are imposed heat flux q1 = 138.5 W/m 2 , q2 = −2272 W/m 2 and q3 = −138.5W/m 2 .The inlet static pressure is pin = 102454 Pa and the static temperature is set to have the profile represented by the formula: where r is the distance from a certain point to the pipe central line.The outlet condition is pout = 100.00kPa.The thermal conductivity is set to be λs = 0.1 W).According to our mesh independent study, the fluid and solid domains are partitioned into 8470 and 6160 elements respectively, as shown in Figure 6.The computational temperature distributions at three different axial locations in Figure 7 appear as parabolic profiles in the fluid region and logarithmic profiles in the solid region, and also show good agreement with the analytical results in [2].

and Analysis in a Turbine Vane
Here we take an application of the constructed platform in the simulation of a convection cooled high-pressure turbine nozzle guide vane named Mark II and analyze the numerical results.2-D simulations of the condition #5411 (Table 1) in the experiment in [27] are carried out.The average temperatures of the cooling air and heat transfer coefficients of the inner cooling surfaces are given according to the experimental results (Table 2).ASTM 310 stainless steel is chosen to be the material of the vane in the simulation.According to the preparative mesh dependence study, we ultimately partition the fluid and solid regions into grids with the number of 18,597 and 6284, respectively, and with the type of mixed rectangular-triangle elements filled (Figure 8).The y + value of the near wall gird is about 1.0.Figure 8 also shows the domain decomposition for the parallel computation.Different models including k-ω, BSL, SST and SST + γ-Reθ are used to make a comparison of the results.The reference pressure is pref = ptin = 337.080kPa and the reference temperature is Tref = 811 K.

Aerodynamic Features
All the numerical models predict very similar aerodynamic performance illustrated by the good agreements between computed and experimental surface pressure distributions shown in Figure 9. On the suction side, the flow accelerates around the leading edge and forms a strong shock at the location at about a half of the axial chord.It accelerates again to supersonic downstream of the strong shock, and forms the second, relatively weak shock near the trailing edge.On the pressure side, the flow gently accelerates in general.Figure 10 shows the pressure contour and isolines around the shock regions, giving a visualized observation of the computed aerodynamic characteristics.The velocity vectors in Figure 11 illustrate that the shock-induced separation occurs behind the strong shock near suction surface in the results of SST and SST + γ-Reθ models but doesn't in those of k-ω and BSL models.The distinction mainly stems from the fact that the SST model adopts the Bradshaw assumption [28] to decrease the computed Reynolds shear stress inside boundary layer and consequently enhance the trend of its separation under adverse pressure gradient.

Heat Transfer Features
The comparison of computed and experimental results of vane surface temperature distributions is illustrated in Figure 12.All the results indicate that on the suction side, the surface temperature reaches a local maximum at the stagnation point around the leading edge and then decreases due to the flow acceleration and the cooling effects of the first three channels.It rises sharply when traversing the strong shock.Behind the shock, the surface temperature appears a rising tendency with fluctuations due to the decrease of the vane thickness and the discrete cooling channels.On the pressure side, the minimum of surface temperature exists at the inflection point of the vane profile and then appears a similar tendency with the suction side.Figure 13 shows the temperature contour in the flow passage and inside the vane.From Figure 12 we can observe the result of the SST + γ-Reθ model has the best agreement with the experimental one in general, especially in the region from the stagnation point to the strong shock on the suction surface.A more detailed observation is shown in Figure 14 which illustrates that in this region the boundary layer appears to have much lower turbulence intensity predicted by SST + γ-Reθ model than the other three (k-ω, BSL, and SST), and the transition to turbulence occurs near the end of this region due to strong deceleration induced by the shock.The laminar region existing in the result of SST + γ-Reθ model makes the computed surface thermal resistance larger than the ones of the other models, and leads to lower (so closer to the real) temperature on the surface and also inside the vane (as in Figure 15).The comparison above implies that numerical models containing no reasonable modes to predict transitions may do not have the potential to precisely simulate this kind of heat transfer problem.  Figure 12 also indicates that all the models overestimate the surface temperature just behind the strong shock.The deviation is conjectured to result from the overestimation of turbulence intensity inside the boundary layer behind the shock.In fact, most two-equation turbulence models without compressibility corrections predict excessive turbulence kinetic energy behind shocks [29].Moreover, since the transition process is triggered around the shock region concurrently, an alternative reason of the temperature deviation for SST + γ-Reθ model may be that the start of full turbulent boundary layer (i.e., the end of the transition process) predicted by the γ-Reθ model is located upstream of the real one, partially due to the model's overlook of the compressible effects on the transition onset and length.The fact that the compressibility can defer the transition onset and increase its length has been realized by some researches (e.g., [30,31]), although complete experimental correlations as a function of Mach number haven't been established entirely.So both turbulence and transition models considering compressible effects are worth exploring in the future.

Conclusions
In this paper, we have presented a primary framework for the simulation of CHT using DG methods on unstructured grids, and a corresponding simulation platform fit for general engineering problems in gas turbines is developed.The centroid-expanded Taylor basis is adopted for various element types.The Bassi-Rebay method is used in the gradient computation inside both fluid and solid domains.A strategy for data exchange through the fluid-solid interface based on the numerical flux of quadrature points is devised to enforce the fully coupled implementation.
The compressible RANS equations combined with turbulence and transition modeling equations, and the heat conduction Fourier's equation, are solved under the unified DG algorithm to realize the applications of the advanced numerical models and algorithm in CHT simulations.An approximate implicit treatment of the source terms in k-ω/BSL/SST models, with the turbulence frequency transportation represented by its logarithm, is derived to decrease the adverse impact rooted in the source term stiffness and make them more compatible with the entire explicit framework.Some cases were studied to validate the feasibility of the method as well as the numerical platform developed.In its application on Mark II, the result of the code combined with the SST + γ-Reθ model showed very good agreement with the experiment for both aerodynamic and heat transfer features.The numerical reasons of the general overestimation of surface temperature just behind the shock were conjectured to mainly result from the ignorance of the compressibility for both the turbulence models and the transition model based on experiment calibration with low Mach number.
The expected advantages of DG methods in the CHT simulations compared with the traditional FV methods, as briefly analyzed in the first section in this paper, need to be further investigated and confirmed, especially when elements with poor quality or non-conformable elements are used near the interface.In addition, more efficient time integration methods will be considered to alleviate the huge computational cost of the DG methods (the computational time consumed by the available explicit DG program was usually about one order of magnitude higher than the FV-based commercial software for these cases).The present work can be viewed as a first attempt towards further investigations of DG and other high-accuracy methods' applications in the engineering CHT problems where their superiority can be seen.
element partition of ΩF, ˆj F and ˆj D are inviscid and viscous numerical fluxes respectively.Assuming that the numerical solution eq hh  U

Figure 1 .
Figure 1.Computational domain and mesh for heat convection along a flat plate.

Figure 7 .
Figure 7.Comparison of temperature distributions of the pipe flow.

Figure 8 .
Figure 8. Mesh partition and succedent domain decomposition (into 32 subdomains) processes for the case of the turbine vane Mark II.

Figure 9 .
Figure 9. Pressure distributions on the vane surface.

Figure 11 .
Figure 11.Near-wall velocity vectors just behind the strong shock wave.

Figure 12 .
Figure 12.Temperature distributions on the vane surface.

Figure 14 .
Figure 14.Turbulence intensity contours in the boundary layers above the suction surface around leading edge of the vane.

Figure 15 .
Figure 15.Temperature contours around leading edge of the vane.

Table 1 .
Test condition parameters of main flow for the case of the turbine vane Mark II.

Table 2 .
Test condition parameters of inner holes boundaries for the case of the turbine vane Mark II.

index Average temperature of cooing air (K) Average heat transfer coefficient (W/m 2 •K)
Pj gradient component of the conservative variable U P m space of all the polynomials with the degree at most m Pk production term of turbulence kinetic energy Pω production term of turbulence frequency Pῶ production term of logarithm of turbulence frequency Pr Prandtl number (set to be 0.71 for air) Prt Turbulent Prandtl number (set to be 0.90 for air) PS pressure surface of vane or blade ReL Reynolds number based on the length of the flat plate Rex Reynolds number based on the distance away from the leading edge of the flat plate u the s-th freedom of the numerical solution of conservative variable Uh in the k-th element , η k js the s-th freedom of the numerical solution of heat flux component qh,j in the k-th element k s  the s-th freedom of the numerical solution of temperature Th in the k-th