Time-Marching Method for Computations of High-Speed Compressible Flow on Structured and Unstructured Grid

The development of Computational Fluid Dynamics (CFD) compressible codes for 2D structured quadrilateral grid and 3D unstructured hexahedral grid is described. The high-speed flow in a nozzle blade cascade is predicted numerically by solving the 2D/3D Euler Equations in a coupled manner. The new finite volume CFD solvers employ second-order accurate central differencing scheme for spatial discretization and multi-stage Runge-Kutta technique for temporal integration with flow variables stored at the vertices. Artificial dissipations with pressure sensors are introduced to control solution stability and capture shock discontinuity. The predictions have been compared with experimental measurements and good agreement has been found. Keyword: CFD, compressible flow, cell-vertex, time marching


INTRODUCTION
The computation of high-speed fluid flow has gained considerable interest among the CFD research community since the past few decades, which is likely to be pioneered by the aerodynamics community in the early sixties such as the early work published by [1] in computing high speed flow using the panel method. More recently [2] has extended the flow models to fully 3D Navier-Stokes equations by solving the set of nonlinear governing equations in decoupled approach, which seems to be the contemporary solution technique for modern CFD. Nevertheless, it is well known that due to the stiffness of the flow governing equations at high Reynolds number, the solution algorithm is found to be destabilizing at stagnation flow regions when the governing equations are solved in decoupled manner as reported by [3] and coupled solution technique is thus more desirable in computing high-speed flow. In addition, the co-existence of subsonic and supersonic region in most of the high-speed flow cases such as flow over a converging-diverging nozzle and flow over a blunt body moving at supersonic velocity further complicates the solution process. It is mathematically known that the steady subsonic flow is governed by elliptic differential equations whereas the steady supersonic flow is described by hyperbolic differential equations. Hence, numerical difficulties exist when the differential equations are solved from elliptic to hyperbolic. The steady state solutions of subsonic and supersonic regions are therefore treated differently and then patched in the transonic region near the intersection of the two regions. In the mid-1960s, a breakthrough happens when [4] developed the time marching procedure. Instead of solving the steady flow governing equations, the equations are advanced in time until steady state is achieved. The idea is that unsteady Euler equations are always hyperbolic and hence no patching is required near the transonic region. Thus, this technique can be applied to mixed subsonic, supersonic and shocked flow cases without a prior knowledge of the flow regimes and the existence of shocks. Because of this advantage, time-marching procedure set together with coupled solution technique is now widely used for high-speed compressible flow as well as incompressible flow, with some minor modification.
Computational grids serve as the essential element in numerical computations. The process in generating the computational grid in a solution domain is known as meshing where the region is divided into many small distinct cells and the flow governing equations are solved over a period of time for each cell. There are two basic types of computational grid, known as structured grid and unstructured grid. A structured grid looks most generally like a twisted coordinate system with each node of a cell represented by an (x, y, z) notation, as denoted by [5][6][7][8][9][10][11][12] . It is a common practice to use structured mesh in a relatively simple geometry because its generation process is easy, both in practicing and programming. The addressing issue in structured environment is straightforward due to the fact that nodes are arranged in an orderly manner. It does have some disadvantages nevertheless. For example, one is restricted to use curved rectangles and this deteriorates the quality of the rectangles particularly in the vicinity of corners and sharp edges. The widely used multiblock structured grid approach solves this problem by simply decomposing the domain into several logical rectangular blocks as recommended by [13,14] . However, the procedure of blocking and structured meshing is a difficult task, as the flow domain becomes complicated enough where domain decomposition into blocks is not possible. It is for this reason that unstructured meshes are used in the current study.
Contrary to structured meshes, unstructured meshes can be in arbitrary shape and become increasingly difficult to work with. In putting the governing equations in discretized form, the trivial but essential information is the neighbouring cell of a computational mesh. This is easily addressed in structured environment, one can just move one element to the left (i-1) or right (i+1) in a typical 1-D computational problem. This becomes more complicated with unstructured grids, particularly the addressing issue; however, there is a payoff here. Unstructured grids can resolve corners and sharp edges easily and hence it becomes prominent to simulate practical flow, particularly as unstructured grid can be generated automatically by using the modern finite element generator in domain of arbitrary complexity as employed by [15] . While unstructured grid has achieved notable success in solving practical flow problems, it also provides flexibility for adapting to flow features such as shock waves and boundary layers, where flow gradient is significant, e.g., [16,17] .
In this study, the CFD solvers for structured quadrilateral (2D) and unstructured hexahedral (3D) will be developed and used to compute the high-speed compressible flow in a nozzle blade cascade using timemarching method.
Governing Equations: In this study, the 3D Euler equations will be given; its application to 2D problem is straightforward. The three-dimensional continuity, x-, y-and z-momentum and energy equations describing the flow of a compressible fluid expressed in strong conservation form may be written as: Where: w is known as the conserved variables, F , G and H are the overall fluxes in x-, y-and z-direction, respectively.

Numerical schemes:
The flow domain is replaced by a finite number of grid points on a mesh system commonly known as quadrilateral mesh (2D) and unstructured hexahedral mesh (3D). The governing equations are solved simultaneously (coupled solution technique) in their integral form for each compact stencil of finite volumes, which can be expressed as: Or in its discretized form: The spatial integration is performed using second order accurate central discretization. A blend of second and fourth order artificial dissipations with pressure switch is added to the residuals prior to the time integration to remove wiggles from the solution. The temporal integration is done using the second order accurate, m-stage (m = 4 in the current study) Runge-Kutta time stepping method proposed by [18] , which can be written as follows: α is the stage-coefficient taken to be [0.250, 0.333, 0.500, 1.000]. To speed up the convergence, 3 types of convergence acceleration schemes are employed: local time stepping, enthalpy damping and implicit residual averaging. At inlet boundary, the total pressure, total temperature and flow angle are fixed while the static pressure is extrapolated from the interior if the inflow is subsonic. Otherwise, all the variables will be specified. At exit, if the exit flow is subsonic, only the static pressure is fixed, while total pressure, total temperature and flow angle are extrapolated from the interior. If the exit flow is supersonic, all four variables are extrapolated from the interior. At the solid boundary, slip condition is used whereby the velocity fluxes normal to the solid boundary will be eliminated.

RESULTS AND DISCUSSION
2D flow simulation in a nozzle blade cascade: In this study, the blade-to-blade flow simulation in a nozzle blade cascade belonged to a stator of a low-pressure steam turbine will be presented. The geometry of the blade was generated using the in-house pre-processor of the current solver developed by [19] . The experimental measurements on the cascade were performed by [20] , which include the surface pressure measurement, wake transverse and flow visualization by Mach-Zhender photography technique.
Three flow cases at overall inlet total to outlet static, P oinlet /P outlet , pressure ratios of 1.49, 1.83 and 2.32 will be examined. The overall pressure ratio of 2.32 corresponds to supersonic outlet, while 1.83 corresponds to transonic outlet. The flow conditions with subsonic outlet are represented by tests at an overall pressure ratio of 1.49.
The mesh illustrated in Fig. 1 consists of 33×230 grids. For inviscid simulation, uniform mesh spacing is employed to maximize the allowable time-step. A comparison of measured and calculated values of blade surface static pressure for subsonic, transonic and supersonic outflows are presented in Fig. 2, 3 and 4, respectively. In general, the turbulent solver [11] performs slightly well in predicting the static pressure along the blade surfaces, except for the pressure values near the trailing edge, as compared to the current inviscid solver. However, by comparing to the experimental data, the effective strategy in capturing the shock wave seems to be the inviscid computation, which accurately represents the impinged shock waves at the suction side of the blade for transonic and supersonic outflow conditions.

3D flow simulation in a nozzle blade cascade:
The subsequent validation case will be the extension of the previous 2D blade-to-blade flow case to 3D. In this study, the calculation is repeated by using the newly developed 3D unstructured hexahedral solver as the objective here is to test the ability of the current numerical code in predicting 3D flow in a linear turbine cascade. The original 2D flow domain was extruded in the z-direction to form an unstructured block consisting of 30000 elements and 34606 nodes with the lowest and highest z-plane indicating the hub and casing, respectively. The extrusion was performed in layering basis to form ten layers of hexahedral cells in the spanwise direction with no mesh refinement near the hub and casing. The hub and casing were set as inviscid walls. The mesh is as shown in Figure 5.
Similar to the 2D test case, three flow cases have been attempted at overall inlet total to outlet static, P oinlet /P outlet , pressure ratios of 1.49, 1.83 and 2.32 which  [20] for subsonic, transonic and supersonic outflow conditions, respectively. Again, considerably good agreement has been obtained. Under subsonic outflow condition, the first pressure rise immediately after the throat was not captured. However, the second pressure rise close to the trailing edge was well captured. In the case of transonic flow, the agreement was better with the sharp rise in pressure due to the normal shock wave, which was well represented by the current prediction. Similarly for the Fig. 8: Comparisons of measured and predicted static pressure distribution for the 3D nozzle cascade at supersonic outflow condition supersonic outflow condition, where the rise in pressure due to the impingement of the trailing edge shock wave originated from the adjacent blade was predicted very well. However, the strength of this shock has been slightly under predicted.

CONCLUSIONS
In the present research, two compressible flow solvers have been developed for 2D structured and 3D unstructured grid by solving the Euler Equations in time-marching coupled manner, which is suitable for computing high-speed compressible flow. It uses the second-order accurate cell-vertex finite-volume spatial discretization and Runge-Kutta temporal integration. The results have been compared with experimental data and good agreements have been achieved. However, the solution is susceptible to wiggles in certain flow region due to the embedded character of the current differencing scheme. In order to overcome this problem, high-resolution differencing scheme is currently being implemented on the flow solvers. Further development of the flow solvers will be the incorporation of implicit time marching scheme, multigrid convergence accelerator and pseudocompressibility factor to simulate incompressible flow using the present time-marching coupled solution technique.