Simulations of Two-dimensional High Speed Turbulent Compressible Flow in a Diffuser and a Nozzle Blade Cascade

This research describes the development of a new structured 2D CFD solver for compressible flow. The high-speed turbulent flow in a diffuser and a cascade of nozzle blade are predicted using standard k-ε turbulence model. The new finite volume CFD solver employs secondorder accurate central differencing scheme for spatial discretization and multi-stage Runge-Kutta time integration to solve the set of non-linear governing equations with variables stored at the vertices. Artificial dissipations with pressure sensors are introduced to control solution stability and capture shock discontinuity. In general, the predictions compare well with the experimental measurements.


INTRODUCTION
In general, there are two methods available in solving the compressible Navier-Stokes Equations (NSE). The most prominent method is to solve the system of non-linear equations in segregated manner by employing an iterative procedure in which solutions are alternately obtained for the pressure and velocity fields. The link is then provided by the continuity equation. In this algorithm, pressure is treated as the primary flow variable due to the fact that pressure gradient is always finite regardless of Mach regimes. Therefore, this is the common scheme employed by the modern commercial CFD codes due to the robustness of the numerical procedure. Known as Pressure Based Method (PBM), this algorithm has been applied extensively in incompressible flow field originally and has been extended to compressible flow by Issa and Lockwood [1] , Van Doormaal et al. [2] , McGuirk and Page [3] , Watterson [4] and Jasak [5] . Notwithstanding this, due to the fact that the momentum, continuity and pressure equations are solved in an uncoupled approach, this may result in convergence problems, especially in situations where the gradients of flow variables are relatively large such as stagnation point at the leading edge. The idea of density variation in compressible flow field has led to the emergence of coupled solution technique since density exists as a dependent variable in the system of compressible Navier-Stokes equations. By realizing the capability of time-marching technique to circumvent the numerical difficulties in mixed subsonicsupersonic problem, time marching procedure has been utilized to solve the system of NSE in a coupled manner by many researchers such as Denton [6] , Ni [7] , Dawes [8] , Ollivier [9] and Lassaline [10] .
The vast majority of fluid applications involve turbulence. Cases such as fluid flow in a pipe, flow processes in combustion chamber, flow over an airfoil will exhibit a chaotic complex motion defined as turbulent flow. While the popularity of Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) have become noticeable due to the rapid development of High-Performance Computing (HPC) technology, the general turbulent fluid motions are well described by the Reynolds-Averaged Navier-Stokes (RANS) equations with the inclusion of Reynolds stresses into the original full Navier-Stokes equation, which is computationally cheaper [8] . To resolve the Reynolds stresses, more equations are necessary and these extra equations are classified as turbulence models.
In this study, the new CFD solver will be used to investigate the high-speed compressible flow in a diffuser and nozzle blade cascade with standard k-ε turbulence model for closure. The idea presented is an extension of the original inviscid 2D solver for twophase steam flow [11] .

MATHEMATICAL FORMULATION
The Governing Equations: The two-dimensional Reynolds-averaged Navier-Stokes equations with k-ε turbulence model as closure describing the turbulent flow of a compressible fluid expressed in strong conservation form in the x-, y-Cartesian co-ordinate system may be written as: Where:

NUMERICAL SCHEMES
Solving Procedure: Starting from the flow field variables obtained from the previous time step, the conserved variables in RANS are solved with the appropriate boundary conditions. The updated variables are then substituted in the standard k-ε turbulence model to solve for turbulent kinetic energy and turbulent dissipation rate. Wall functions and turbulent boundary conditions are then imposed. The updated k and ε will be used to calculate the turbulent viscosity and the Reynolds stresses. Subsequently, the new Reynolds stresses will be utilized to solve RANS in the next iteration. The loop continues until convergence is achieved.

Cell-vertex finite-volume spatial discretization:
The flow domain is replaced by a finite number of grid points, which are generated algebraically by the built-in pre-processor [12] . The mesh system is commonly known as H-mesh and divides the physical domain into a set of discrete rectangular control volumes.
A cell-vertex formulation is used in which the flow variables are stored at four vertices of a quadrilateral cell. It has been shown by Martinelli [13] , Dick [14] , Swanson and Radiespiel [15] that cell-vertex formulation offers some advantages over the cell-centred one in which cell-vertex method offers higher accuracy on irregular grid. For a uniform mesh, there would be no difference between the cell-centred and cell-vertex schemes; however, cell-vertex scheme does not require extrapolation to the solid boundary to obtain the wall static pressure, which is necessary in solving the momentum equations for cells adjacent to the solid boundary.
Starting from known values of primitive variables from the previous time-step, the values of F and G are calculated for each node. Then the line integration is performed for each control volume in turn for the six conserved variables (RANS + k-ε Turbulence Model).

The discretized RANS:
(2) The discretized standard k-ε model: - After the spatial integration, the cell residual will take the form: represents the sum residuals from RANS or turbulence model.
The calculated residuals apply to the values of properties within the cell, whereas, the variables are actually stored at the nodes. Consequently, they have to be redistributed to the four surrounding nodes. This is done by sharing the changes equally between the four nodes in the context of central differencing.
Thus, the equivalent discretized equation for a node will be: Artificial dissipations: All second-order centraldifferencing schemes, even with a stable time-step, suffer from certain tendencies to instability due to the odd-even decoupling near a discontinuity. The scheme can be stabilised by introducing a small amount of artificial viscosity, suggested by Jameson et al. [16] . First -order upwind differencing scheme may be used to remedy the stability problem, however, the scheme tends to damp the solution so much and alter the flow physics. Therefore, 2 nd order accurate central differencing scheme, which is consistent to the framework of Navier-Stokes equations, is applied to both RANS and turbulence model in the current work with suitable amount of artificial viscosity.
This artificial viscosity formulation is a blend of second and fourth-order terms with a pressure switch to detect changes in pressure gradient [16] . After the addition of the dissipation terms, Equation (5) becomes: The multi-stage runge-kutta time stepping: Equation (6) is integrated with respect to time by means of a fourstage Runge-Kutta time stepping scheme, as proposed by Jameson et al. [16] : Boundary conditions: The mathematical theory of incompletely parabolic PDEs indicates the number and type of boundary conditions for the unsteady compressible RANS. In the present work, Euler-type of boundary condition is applied except at solid wall where no-slip condition is imposed. At inlet, the total pressure, total temperature, flow angle, turbulent kinetic energy and turbulent dissipation rate are fixed, while the static pressure is extrapolated from the interior if the inflow is subsonic. Otherwise, all variables are fixed at inlet. At exit, if the outflow is subsonic, only the static pressure is fixed, while total pressure, total temperature and flow angle are extrapolated from the interior by zero-order extrapolation. If the exit flow is supersonic, all four variables are extrapolated from the interior. For k-ε transport equations, k and ε are extrapolated for arbitrary exit flow conditions. Periodic boundary condition is essential in simulating flow in turbomachines. The periodicity condition on the bounding streamlines is easily satisfied by treating the calculating points on each of the bounding streamline as if they are interior ones, by assuming that all properties are equal for corresponding points on each of the streamline. Nodes on the periodic boundaries will have contributions from four corresponding cells on both sides of the boundaries.
On solid walls, the values for velocity components as well as k and ε are set to zero. Adiabatic condition is imposed. For nodes adjacent to the wall, wall functions are introduced to calculate k and ε .
Initial conditions: To start the computation, initial flow field variables must be specified at all calculating points. In this approach, a linear variation of pressure between inlet and exit planes is assumed from which the pressures at all calculating points can be obtained. The tangency condition is enforced to obtain the velocity components at each calculating point and no variation of other properties along the pitch is assumed. Using the inlet stagnation temperature and pressure, the assumed static pressure and velocity components, other properties can be calculated using isentropic relations.
For turbulence models, k and ε are set to the values consistent to the inlet k and ε . The fluctuating Reynolds stresses are then calculated using Boussinesq relation.

Blade-to-blade calculations on a turbine nozzle cascade:
In order to validate the current solver, bladeto-blade flow simulation on a turbine nozzle blade cascade will be presented. The blade profile belongs to a stator of a low-pressure steam turbine. The geometry of the blade was generated using the in-house preprocessor of the current solver [12] , as shown in Fig. 1. The experimental surface pressure measurements on the cascade were performed by Mamat [17] . Three cases at overall inlet total to outlet static pressure ratio, P oinlet /P b of 1.49, 1.83 and 2.32 were simulated. The overall pressure ratio of 2.32 corresponded to supersonic outlet, while 1.83 corresponded to transonic outlet. The flow conditions with subsonic outlet were represented by tests at an overall pressure ratio of 1.49.
The mesh consisted of 33 x 230 grids. The mesh resolution near the wall was adjusted to be higher to account for the boundary layer development. A comparison of measured and calculated values of blade surface static pressure for subsonic flow, transonic flow and supersonic flow are illustrated in Fig. 2-4. In general, the numerical results show good agreement with the experimental data, except the pressure ratio at the trailing edge due to mesh distortion. Sajben diffuser: Transonic turbulent flow has been computed in a two-dimensional converging-diverging duct using the standard k-model. Extensive experimental data are available for this geometry, at a variety of flow conditions (Chen, Sajben and Kroutil, [18] , Bogar, Sajben and Kroutil [19] ,  Salmon, Bogar and Sajben [20] , Sajben, Bogar and Kroutil [21] ,Bogar [22] ). The flow fields being modelled were the weak-and strong-shock diffuser cases of Sajben [23] . The geometry of the diffuser is presented in Fig. 5. A 51x81 body-fitted computational mesh was generated algebraically. Adiabatic no-slip conditions were used on both the top and bottom walls. The pressure ratio (P exit /P oinlet ) was 0.82 and 0.72 for strongand weak-shock case, respectively. Figure 6 and 7 compares the pressure distributions along the bottom and top wall of the diffuser. The shock location predicted by the current solver compares well with the experimental data carried out at P exit /P oinlet = 0.82 (weak shock). Next, the strong-shock case was simulated. Again, the illustrated results are comparable with the experimental data, except that the shock is predicted to occur 1 grid point further downstream as shown in Fig. 8 and 9. However, by considering the coarseness of the mesh employed in the current solver, the result is satisfactory.

CONCLUSION
In the present work, a new two-dimensional compressible flow solver has been developed for structured grid. It uses the second-order accurate cell-vertex finite-volume spatial discretization and Runge-Kutta temporal integration. Standard kturbulence model has been successfully adopted to simulate the compressible turbulent flow in a cascade of nozzle blade and Sajben diffuser. Both cases show good comparison with the experimental data, except the excessive smearing of shock wave at the trailing edge of the nozzle blade cascade. Research is still in progress on the existing turbulence model to compute the turbulent viscosity in a more accurate way.
Further work to be done on the solver includes the extension to 3D environment and modification to handle arbitrary meshes.