# An Implicit Gas-Kinetic Scheme for Turbulent Flow on Unstructured Hybrid Mesh

###### Abstract

In this study, an implicit scheme for the gas-kinetic scheme (GKS) on the unstructured hybrid mesh is proposed. The Spalart-Allmaras (SA) one equation turbulence model is incorporated into the implicit gas-kinetic scheme (IGKS) to predict the effects of turbulence. The implicit macroscopic governing equations are constructed and solved by the matrix-free lower-upper symmetric-Gauss-Seidel (LU-SGS) method. To reduce the number of cells and computational cost, the hybrid mesh is applied. A modified non-manifold hybrid mesh data(NHMD) is used for both unstructured hybrid mesh and uniform grid. Numerical investigations are performed on different 2D laminar and turbulent flows. The convergence property and the computational efficiency of the present IGKS method are investigated. Much better performance are obtained compared with the standard explicit gas-kinetic scheme. Also, our numerical results are found to be in good agreement with experiment data and other numerical solutions, demonstrating the good applicability and high efficiency of the present IGKS for the simulations of laminar and turbulent flows.

###### pacs:

47.45.Ab, 02.70.-c, 47.11.Df, 47.27.E-## I Introduction

Over the past half-century, the gas-kinetic theory has been extensively studied for incompressible and compressible fluid flows. As an alternative method based on the Bhatnagar-Gross-Krook (BGK) collision model, the gas kinetic scheme (GKS) proposed by Xu and Prendergast xu1994numerical () has been rapidly developed in the last two decades and has attracted an increasing amount of attention from the CFD community. Through new modifications as reported in the literature xu2001gas (), the revised GKS scheme is of a certain advantage of resolving dissipative structure, especially in capturing shock waves. In the traditional Navier-Stokes (NS) solvers, the gas is assumed to stay in two equilibrium states on both sides of the shock wave and the shock wave appears as a discontinuity; that is to say, the physical dissipation in a cell size is replaced by numerical one. In the flux reconstruction approach, an additional artificial dissipation is also introduced by NS solver due to application of upwind scheme and/or central difference method. To remove the spurious dissipation, a general non-equilibrium state is considered for constructing the initial gas distribution function at each time step and the equilibrium state at a cell interface utilizing the Chapman-Enskog expansion xu2001gas (). In the complex flow simulations, the GKS can give us more real description from the aspect of physical evolution than traditional NS solvers.

In practice, as reported by May et al. may2007improved (), the flux evaluation in GKS is slightly more complicated compared with that generally used in the other finite volume method (FVM) solvers, which means that GKS takes more computational cost than other FVM solvers under the same situation. So, in practical applications, acceleration methods for GKS must be developed. To reduce the computational time, the parallel implementation is a good strategy. Ilgaz and Tuncer ilgaz2006parallel () investigated the parallel implementation of the gas-kinetic BGK scheme on unstructured grids by using the domain decomposition method. Kumar et al. kumar2013weno () proposed a parallelized WENO-enhanced GKS for the three-dimensional (3D) direct simulations of compressible transitional and turbulent flows.

On the other hand, for the steady flow simulations concerned in this study, using implicit scheme is an efficient way to accelerate convergency and reduce computational cost. Chit et al. chit2004implicit () applied an approximate factorization and alternating direction-implicit (AF-ADI) method for the GKS simulations of the inviscid compressible flow on structured grids, in which an improved algorithm was achieved with fast convergence and the capability for using a large Courant-Friedrichs-Lewy (CFL) number, and this implicit method exhibited better results than the traditional explicit method. Recently, Li et al. li2014implicit () presented an implicit gas-kinetic method with the matrix-free lower-upper Symmetric Gauss-Seidel (LU-SGS) time marching implicit scheme to simulate the hypersonic inviscid flows on unstructured mesh, in which the implicit CFL number was determined by decreasing residual and the good robustness of this kind of implicit gas-kinetic scheme (IGKS) was validated in their works. However, the convergent behavior was not shown with the residual curve in their works. As mentioned above, several kinds of IGKS have been applied for compressible flow simulations chit2004implicit (); li2014implicit (). For unsteady flows, a dual time-stepping strategy of gas-kinetic scheme is proposed by Li et al PhysRevE.95.053307 (). In their works, fluid flows from laminar to turbulent and from incompressible to compressible are accurately simulated with better efficiency than previous method. However, there were few studies on implementing IGKS on unstructured hybrid grids for the compressible turbulent flows around/in complex geometries. This study focus on presenting an IGKS with the LU-SGS scheme on unstructured hybrid grids for complex turbulent flows, and the convergent behavior of IGKS is also investigated in detail.

For the simulations of compressible flow around/in complex geometries, the unstructured hybrid grid becomes a promising choice, due to its advantage in balancing the accuracy and the computational cost. For example, for flow around complex realistic configuration, the body fitted grid can be used to resolve the boundary layer region, while the unstructured grid with a suitable growth rate can be applied to fill all other computational domain. Besides, the grid refinement and coarsening are quite easy to be implemented mavriplis1997unstructured (), compared to the structured or block-structured grid technique which relies on regular connectivity of quadrilateral or hexahedral cells. In problems with complex configurations, high-order finite volume method under unstructured grids can also obtain more elaborate and precise results liu2016high (). In compressible cases, limiter for unstructured grid is easy to implement liu2017accuracy (). As we known, the complexity of hybrid grids, including elements, edges, nodes, and connectivity, needs an efficient mesh data structure to reduce the extra computational costs. In this study, the extended non-manifold hybrid mesh data (NHMD) ebeida2009mesh (); FLD4239 () is employed to build an accessible library for mesh.

Most flows encountered in engineering applications are of turbulent nature. For the compressible turbulent flow simulations, compared to the direct numerical simulation (DNS) and the large eddy simulation (LES), the turbulence models require coarse grid resolution and have been validated to be suitable for the turbulent flow simulations in engineering. In particular, the one-equation Splalart-Allmaras (SA) model spalart1992one () became quite popular because of its satisfactory results for a wide range of flow problems and its reliable numerical properties. On the basis of the explicit GKS, we have successfully carried out the simulation of compressible turbulent flows with shock waves on unstructured meshes, in which the SA turbulence model FLD4239 () and SST turbulence model PhysRevE.95.053307 () are incorporated to include the effect of turbulence. In this study, the SA turbulence model is also selected on account of its good accurate and stability under relatively coarse grid near the wall, as well as its low computational cost. Other choices, such as the two-equation models, might also be good candidates, which will be developed in the future.

In this paper, an IGKS coupled with SA turbulence model is proposed for the incompressible and compressible flow simulations on unstructured hybrid mesh, and several flow simulations including the lid-driven cavity flow, the laminar and the turbulent flow around a flat plate, as well as the turbulent flow around a multi-element airfoil are performed here.

The rest of the paper is organized as follows. The basic GKS proposed by Xu and Prendergast xu1994numerical (), the matrix-free LU-SGS scheme and the implementation of IGKS, the expanded NHMD structure and the coupled SA turbulence model, as well as four different boundary conditions are described in Section II. Then, the numerical validations and several turbulent flow cases are carried out to show the accuracy and reliability of the present IGKS-NHMD method in Section III. Finally, some remarks concluded from this study are grouped in Section IV.

## Ii IMPLICIT METHOD for GAS-KINETIC SCHEME

### ii.1 Gas-Kinetic Scheme

The two-dimensional (2D) gas-kinetic scheme based on BGK model is written as xu1994numerical ()

(1) |

where and are the gas distribution functions, which are the functions of space , particle velocity , time , and internal variable . is the particle collision time. is the Maxwell distribution function which has the following form

(2) |

where is the internal freedom degree with for the 2D diatomic molecule gas flows. The variable , is the molecular mass, is the Boltzmann constant, and is the temperature. is the density, and are the and components of the macroscopic velocity in 2D, respectively. Note that, for the gas system with K freedom degree, the square of internal variable can be taken as

(3) |

On the cell interface, the general solution of Eq. (1) at time can be given as xu1994numerical (); xu2001gas ():

(4) |

where the particle trajectory . For simplicity, setting t the cell interface in the following text.

To include the physical nature of non-equilibrium phenomena in the complex flows, the deviation of a distribution function away from a Maxwellian velocity distribution should be applied to construct the gas distribution function. In Eq. (4) the initial gas distribution function at the beginning of each time step is written as

(5) |

where and represent the Maxwellian distribution functions on the left and right sides of interface, , , , , and are related to the derivatives of the Maxwellian velocity distribution in space and time. The additional time and the spatial derivative terms in Eq. (5), and , denote the non-equilibrium effects, which have no direct contribution to the macroscopic conservative variables.

All the derivative terms in Eq.( 5), including , , , , , and have the form of Taylor expansion

(6) | ||||

These terms can be obtained from the relations between the gas distribution functions and the macroscopic variables.

(7) | ||||

where , is the vector of collisional invariants, and the macroscopic variables .

To reach second-order accuracy, the equilibrium distribution should contain the first-order spatial and time derivative terms,

(8) |

where , , , and are spatial derivatives and time derivative of distribution function for equilibrium state. In Eq. (5) and Eq. (8), , , and represent left, right of the interface and on the interface, respectively.

In the present simulations on the unstructured mesh, all spatial derivative terms in Eq. (5) can be obtained by the least squares method. Taking the simple schematic grid shown in Fig. 1 for example, firstly, the fitting formulation is applied

(9) |

To determine the derivative terms and , the least-squares regression equations can be reconstructed as

(10) |

where denotes the total number of cells adjacent to the cell . Then, the normal and tangential derivatives of the macroscopic variables can be computed from the known derivative terms and . Finally, the concerned derivative terms (e.g. , ) can be computed by following the detailed process reported in Ref. xu2001gas ().

For compressible flow simulations, the limiter is often introduced to compute the gradients to avoid a huge value and/or to eliminate the possible spurious oscillations. In this paper, the minmod limiter sweby1984high (); de1993quadtree () is used for updating the derivatives of macroscopic variables,

(11) |

where the and are the minimum and maximum value of around the cell, which are used to calculate the gradients of at the cell. The and are the minimum and maximum value of at the node of the cell, which are reconstructed by the gradients. represents the macroscopic variables at the center of the cell. Note that, all the macroscopic variables have their own minmod variety coefficients , and Eq. (9) can be rewritten as

(12) |

To determine the time derivative terms, , and , the flow conservation constraints must be utilized: the non-equilibrium part of the gas distribution function should have no direct contribution to the macroscopic conservative variable

(13) |

and the macroscopic variables mass, momentum, and energy should be conserved during particle collisions,

(14) |

Thus, the systems of linear equations about those derivative terms can be built, and can be solved by the Gaussian elimination.

Now, the gas distribution function at the cell interface can be expressed as

(15) | ||||

where the Heaviside function reads

Finally, we can obtain the flux term of the macroscopic variables across the cell interface from the gas distribution function at the cell interface, e.g., the flux in x-direction can be computed as

(16) |

### ii.2 Matrix-Free Lu-Sgs Scheme

Generally, the explicit time marching of finite volume method for the conserved variables are written as

(17) |

where is the volume of the control volume. represent the area of the interfaces of this control volume.

In this study, considering an implicit scheme for the time marching. Firstly, Eq. (17) can be rewritten as

(18) |

Following the implicit algorithm presented in Ref. li2014implicit (), for simplicity, the right hand term of Eq. (18) is referred to as a residual term , and Eq. (18) can be expressed as the differential form

(19) |

Then, above equation can be further rewritten as

(20) |

where the term is a Jacobian matrix, which takes the form in first-order Roe’s scheme in which the Roe’s flux can be linearized chen2000fast (). According to the idea of LU-SGS method, we first split this Jacobian matrix into three parts: lower triangular matrix, upper triangular matrix and diagonal terms. Thus, Eq. (20) can be rewritten as

(21) |

with

(22) |

where reads

(23) |

with represents velocity vector, denotes the sound speed at the cell , and is area of the interface .

Then, the LU-SGS scheme is applied to solve Eq. (21):

(24) |

When computing the RHS term at time n-level, the n-level’s explicit time step should be known. Here, this explicit time step is computed by

(25) |

To accelerate convergency (time marching), an implicit time step is introduced with the similar form,

(26) |

The implicit CFL number can be adjusted with the decreasing residual. We choose infinity norm of residual to be the criterion that how to adjust the implicit CFL number tidriri1997preconditioning ().

(27) |

### ii.3 Turbulence Model

In order to solve the turbulent flow problems, the SA one-equation turbulence model is employed, which solves a transport equation for a transformed eddy viscosity . In this study, the governing transport equation followed in the work by Moro et al. moro2011navier () reads

(28) |

where the three terms at the right hand side represent turbulence production, turbulence destruction due to wall and turbulence diffusion/propagation, respectively. For the computation/definitions of the eddy viscosity and all other parameters and constants involved in the SA model, please find the details in Ref. moro2011navier (). Here, we focus on introducing how to solve this SA transport equation on unstructured hybrid mesh by FVM.

Following the approach proposed in Ref. shenderobust (), a semi-discrete form in cell for the discretization of SA model equation reads

(29) |

where is volume of element , denotes the source terms,

(30) |

and denotes the flux term which including three parts: inviscid flux , viscous flux and anti-diffusion flux ,

(31) |

where represents the area of the th interface of cell and we define

(32) |

with

(33) |

Note that, in above two equations, represents the distance vector from the center of cell point to the center of cell which shares the th interface and denotes unit normal vector of the cell interface.

Now, Eq. (29) can be rewritten as the following implicit form

(34) |

For the derivative terms, taking cell and its neighbors into consideration, Eq. (34) becomes

(35) |

Here, we can obtain a system of implicit linear equations which can be solved by using the LU-SGS scheme.

In BGK model bhatnagar1954model (), the particle average collision time equal to the ratio of the dynamic viscosity to the static pressure,

(36) |

In order to couple the turbulence effect into the present GKS model, the computation of relaxation time should be modified by incorporating the eddy viscosity into the total viscosity xiong2011numerical ()

(37) |

where and denote the laminar dynamic viscosity and the eddy dynamic viscosity, respectively, and is computed from the updated values of the density and the eddy viscosity , .

### ii.4 Non-Manifold Hybrid Mesh Data Structure

For the numerical simulation on unstructured hybrid mesh, an efficient mesh data structure is required to access all mesh data without causing high cost. In the present study, the non-manifold hybrid mesh data (NHMD) ebeida2009mesh () is applied, in which six kinds of elements are considered: nodes, lines, edges, faces, cell and entity. These elements form a library and can be indexed each other ebeida2009mesh (); FLD4239 (). Thus, we do not need to repeat the process of searching for geometrical information.

For example, as shown in Fig. 2, NHMD is categorized into three parts: 1) the number of nodes and every kind of cells; 2) the locations of all nodes denoted by , , ; 3) the connectivity of every cell represented by the index number of all vertex of this cell. By using this method, we can easily convert arbitrary hybrid mesh into this kind of mesh data for the present IGKS solver since there is no requirement for the geometric features of every cell.

### ii.5 Boundary Conditions

Following the “ghost cell” method shown in Xu’s studies xu2001gas (), this subsection is mainly about how to implement this idea numerically for the present IGKS. Here, four kinds of boundary conditions including free-stream, symmetric, non-slip wall and outflow boundary conditions are introduced, and both physical values and their increments are considered due to the used implicit scheme.

#### ii.5.1 Free-stream boundary condition

For the physical variables including four conservative variables and the transformed eddy viscosity , simply set as

For increments of these physical variables,

#### ii.5.2 Symmetric boundary condition

For the physical variables,

Similarly, for increments of the physical variables,

#### ii.5.3 Non-slip boundary condition

For the physical variables,

For increments of the physical variables,

#### ii.5.4 Outflow boundary condition

For the outflow boundary conditions, the Riemann invariants presented by Carlson carlson2011inflow () are introduced in the present treatment of outflow boundary condition to reduce the disturbance propagating into flow field. As reported in Ref. carlson2011inflow (), if the physical state (outside the domain) is on the right side of boundary and the solution space (inside the domain) is on its left, we define two Riemann invariants correspond to the incoming and outgoing characteristic waves, respectively

where the two subscripts and represents the normal values/components outside and inside the flow field, respectively. For the subsonic outflow boundary conditions, both incoming and outgoing characteristic waves exist, so the velocity and the speed of sound at the boundary are computed as the sum and difference of the invariants, respectively

Then, the density and energy on the boundary can be calculated as follows

where represents inside the flow field, and the transformed eddy viscosity on the boundary just set as .

For the increments of macroscopic variables, all are simply set at zero. This strategy does not produce negative effect in convergence.

## Iii Numerical Results

In this section, we perform some numerical simulations to validate the present IGKS under unstructured hybrid mesh in different types of fluid flow; from laminar to turbulent flows. In all simulations of steady flows, the convergent criterion towards the steady state is set as follows

(38) |

where is the number of mesh cells in fluid field.

### iii.1 Lid-driven cavity flow on uniform and unstructured hybrid meshes

The 2D lid-driven cavity flow has been studied by many researchers and employed as a benchmark case of incompressible viscous flow. In this part, the numerical simulations for the 2D lid-driven cavity flow on uniform and unstructured hybrid meshes are performed.

The flow configuration of this case is consists of a 2D square cavity with its top plate moving with a uniform horizontal velocity , while other three walls are set to be static non-slip boundary. In this series of simulations, the Reynolds number is defined as and set at unless otherwise notified, in which the height/width of the cavity is defined as the reference length, , and the Mach number is set to be . Besides, the uniform grid and an unstructured hybrid grid with a total number of mesh cells, are respectively used for the present two simulations. Here, the unstructured hybrid grid is composed of three kinds of cells, quadrilaterals, pentagons, and hexagons, as shown in Fig. 3. In the present 2D lid-driven cavity flow simulation with the unstructured hybrid grid, the mesh shown in Fig. 3 is the final adaptive mesh which are generated and controlled by performing refinement in the region where the gradients and vorticities are larger than their mean level.

Fig. 4 shows the two computed streamlines (stream functions ). It is clear that the same flow structures are obtained from the two simulations with different meshes, which are in good agreement with those results shown in previous literatures ghia1982high (); zhuo2012filter (); hou1995simulation (). Fig. 5 presents the comparisons of the computed velocity profiles along the mid-height () of the cavity against the numerical results of multi-grid method reported by Ghia et al. ghia1982high (). It can be found from this plot that both results of the present simulations with uniform and unstructured hybrid grids agree well with the benchmark data of Ghia et al. ghia1982high ().

In order to demonstrate the accuracy of the present IGKS, the locations of the primary and secondary vortices, as well as the values of vorticity at these locations are grouped in Table 1. The values of locations and vorticity are normalized by and , respectively. For purpose of comparison, the numerical results reported by Ghia et al. ghia1982high () and two lattice Boltzmann method (LBM) results obtained with uniform grid of zhuo2012filter (); hou1995simulation () are also grouped into this table. The present comparisons show that, both IGKS results obtained with uniform and unstructured hybrid grids are in good agreement with those results given in the previous literatures ghia1982high (); zhuo2012filter (); hou1995simulation (): the relative differences between the present results and those benchmark data of Ghia et al. ghia1982high () are less than , and for the values of locations, vorticity and stream function , respectively.

To further compare the present results obtained with two different kinds of meshes, Fig. 6 shows the computed pressure contours of two simulations at the same contour levels. It is clear that, the two pressure contour lines agree well in the most flow regions, even near the top-right corner in which there is singular point on the top-moving wall. Also, there are a slight differences between both pressure contour lines at the primary vortex, near the wall and the top-left corner. Overall, this comparison shows that the present IGKS with unstructured hybrid grids can give us the satisfactory results.

For the present two simulations with different kinds of meshes, the variations of the residual value determined by Eq. (38) and the CFL number over the iteration step are grouped in Fig. 7 and Fig. 8, respectively. It can be observed from the two plots that, the two simulations with uniform and unstructured hybrid grids reach the steady state quickly and only take about and iteration steps, respectively,, and both variable CFL numbers grow very fast with setting the initial implicit CFL number and their final values increase up to the order of one million.

When compared to the explicit GKS, it is no doubt that implicit iteration in present IGKS needs more computational time (takes more time than the explicit GKS per iteration step, which can be calculated from the data grouped in Table 2) to complete one iteration step, but IGKS can achieve convergence much faster in the simulations of steady flows since the quite large CFL number can be used in the IGKS simulations. To show the advantages of the present IGKS and assess its computational efficiency, another simulation is also performed by using the standard explicit GKS with the uniform mesh for the purpose of comparisons. Fig. 9 shows the comparison of convergency curves for the present IGKS and explicit GKS simulations. It is clear that, the IGKS simulation only take iteration steps to converge while about iteration steps are need for the explicit GKS simulation. This result is consistent with the work of Zhu et al (2016) which is based on the unified gas kinetic scheme (zhu2016implicit, ). Considering the computational efficiency, the computational time of the present IGKS and the explicit GKS simulations on uniform and unstructured hybrid grid have been computed and grouped in Table 2. It is noted that, the total numbers of mesh cells and the cost of computational time for running the corresponding iteration steps are included here for all simulation cases performed in this paper. Table 2 shows that, for the present steady cavity flow cases, the IGKS simulations are about two orders of magnitude faster than the explicit GKS simulations. Meanwhile, due to few mesh cells used for the unstructured hybrid grid, the computational time of the present simulations of cavity flows on unstructured hybrid grid is just about of those taken by the corresponding simulations on uniform grid.

In all, it can be concluded from these comparisons that, the present IGKS with unstructured hybrid grid not only has capability for accurately predicting the steady incompressible viscous flows, but also has good computational efficiency compared with the explicit GKS.

### iii.2 Laminar flows around a zero-pressure-gradient flat plate

In this part, we consider a laminar boundary flow with Mach number and over a flat plate. In Fig. 10 the schematic view of the computational domain is given along with the corresponding physical boundary conditions. It is noted that, the length of flat plate is set to be the characteristic length, , and the leading edge point locates at ; the inflow at the left boundary is set as uniform with a constant horizontal velocity, and the Riemann invariants are considered for the upper boundary and the outflow boundary; for bottom boundary, the symmetric and non-slip boundary conditions are considered.

As we known, the mesh resolution for boundary layer is crucial, thus using hybrid mesh is a wise strategy to enhance mesh efficiency and reduce computational cost as low as possible. As shown in Fig. 11, a hybrid mesh with cells is used in this simulation, and the local view near wall is plotted in Fig. 11. In Fig. 12, the computed U and V velocity profiles at three locations and are plotted, where is the dimensionless distance from the flat plate and the exact Blasius solutions are also grouped into these plots. Fig. 13 shows the computed skin friction coefficient compared with the exact Blasius solution. In those comparisons, the present IGKS results of the velocity profiles and the skin friction coefficient profile fit the exact solutions very well.

Besides, for this case, the IGKS simulation converges at th iteration steps and the corresponding computational time is about seconds, and which also shows about two order of magnitude improvement in the computational times when compared to that of explicit GKS, as shown in Table 2. The history curves of the CFL number and residual value are plotted in Fig. 14.

### iii.3 Turbulent flow around a zero-pressure-gradient flat plate

In this part, to validate the present IGKS-SA method in simulation of turbulent flows, we continue to investigate the turbulent boundary layer case with Mach number and over a flat plate. The delicate flow structures in boundary layer, such as velocity profile, skin friction coefficient and eddy viscosity, are computed for comparisons.

In this test, the length of the flat plate is 2.0, the height of the whole flow field is 1.0, and the leading edge of the flat plate is at . The final refined mesh reported in our previous simulations with the explicit GKS FLD4239 () is used in the present simulation, where the total number of cells is and most of the cells gather within the boundary layer, as shown in Fig. 15. Near the wall boundary, the high resolution of mesh is guaranteed where the first grid point off the wall is located at , the local view of mesh near the wall is shown in Fig. 15.

For this case, all boundary conditions used in above laminar boundary layer case are still applied. For the SA turbulence model, we set for the left inflow boundary and the initial flow field allmaras2012modifications (), while set for the top, outflow and wall boundaries.

Fig. 16 and 16 show the present computed velocity profiles at two positions and , respectively. The benchmark solutions of CFL3D reported by Wilcox LangleyResearchCenter () are grouped into the two plots for the purpose of comparisons. It is clear that, both profiles computed by the present IGKS are in good agreement with the benchmark data of CFL3D. Note that, in Wilcox’s study LangleyResearchCenter (), the number of cells used in the simulation is which is about eight times the present computational grid. In Fig. 17 and Fig. 18, we illustrate the present results of the skin friction coefficient and the dimensionless eddy viscosity profile at , respectively, which also agree well with those benchmark solutions of CFL3D LangleyResearchCenter ().

In Fig. 19, we plot the variations of CFL number and residual value over the iteration. The present IGKS simulation can reach convergency at about iteration steps, while for the previous explicit GKS simulation, it takes up to about iteration steps (see Table 2). To investigate the computational efficiency of the present IGKS for the turbulent flows, the computational time of the present IGKS simulation and the previous simulation of the explicit GKS are also included into Table 2, which shows that the present IGKS simulation only takes of the computational time of the explicit GKS simulation. It should be pointed out that, compared to the value of observed from the laminar flow cases, the ratio of computational time () between explicit GKS and IGKS in the turbulent flow simulations is a little bit smaller since the implicit scheme is also considered for the SA turbulence model.

### iii.4 Turbulent flow around NHLP multi-element airfoil

In order to validate the applicability of the present IGKS-SA solver for predicting the complex turbulent flow with complex boundaries, a turbulent flow around a two-dimensional supercritical airfoil with high-lift devices, NHLP-2D three-element airfoil, is investigated in this study.

The configuration of NHLP-2D includes a leading-edge slat and a single slotted-flap, where is the chord length of the nested configuration. Fig. 20 shows this geometry with the computational grid near airfoil. It is clear that, to reduce the total number of mesh cells, the cell scale increased fast within the region out of boundary layer and the outer of wake region. However, in boundary layer and near-wall wake regions, the fine grids are used in this study; as shown in Figs. 21 and 21, the body-fitted quadrilaterals are applied to capture the thin boundary layer and the fine triangles are used to fill the two gaps between slat and flap with the main element. Here, the first node point off the airfoil surface for the finest grid is located such that , and the total number of mesh cells used in this simulation is . The flow conditions are Mach number , Reynolds number and an angle of attack of .

The contours of density, pressure, Mach number and eddy viscosity computed from the present IGKS-SA solver are shown in Figs. 22, 23, 24 and 25, respectively. It can be found from those plots that, the slat wake merges with the boundary layer of main element and seems to be disappearing near the trailing edge of main element. Fig. 26 presents the surface pressure coefficients calculated with present IGKS-SA solver, which are in good agreement with the experiment data reported by Morrison morrison1998numerical (). It is important to point out that the present SA turbulence model is not calibrated to predict transition, and actually, the fully turbulent simulation is performed here.

In Fig. 27, the variations of CFL number and residual value over the iteration are illustrated. The present IGKS simulation can reach convergence at iteration steps, while for the previous explicit GKS simulation, it needs about iteration steps (see Table II). To investigate the computational efficiency of the present IGKS for the turbulent flows, the computational time of the present IGKS simulation and the previous simulation of the explicit GKS are also included into Table II; the present IGKS also exhibits times improvement than the classic explicit GKS.

In all, the present IGKS solver is capable of the laminar flow simulations on both regular structured grid and unstructured hybrid mesh, and when coupled with the turbulence model, e.g. SA model, the present IGKS solver shows good convergence behavior and accuracy for simulating the complex turbulent flow with complex boundaries which contains arbitrary complex configurations.

## Iv Conclusions

This paper presents an implicit GKS coupled with Splalart-Allmaras turbulence model based on LU-SGS scheme under unstructured hybrid mesh. As a kind of efficient mesh data structure for unstructured hybrid mesh, an extended NHMD is used in this study to reduce the extra computational costs. In this paper, to validate the capability of the present IGKS solver and show the advantage of the present IGKS compared with the explicit GKS, four steady viscous flow cases including lid-driven cavity flows, laminar and turbulent boundary flows over a flat plate, and the turbulent flow around NHLP-2D airfoil are performed.

For the lid-driven cavity flows, both simulations on uniform and unstructured hybrid mesh are considered, and almost the same results have been obtained by the present IGKS solver which also agree well with the benchmark data reported by other researchers. In the present simulation for the laminar boundary flow over a flat plate, the velocity profiles at three different locations and the skin friction coefficient profiles along the flat plate computed from the present IGKS solver agree quite well with the Blasius analytical solutions; also, for the turbulent boundary flow case, the delicate flow structures in boundary layer are predicted by the present IGKS coupled with SA model, show us the satisfactory results compared with the available data of CFL3D. At last, the steady turbulent flow around the NHLP-2D airfoil is considered to validate the capability of the present IGKS-SA solver for predicting the complex flow with complex configurations. The merging of the slat wake with the main element boundary layer are clearly observed from the present simulation and the present results of surface pressure coefficients are also in good agreement with the available experimental data. Besides, the comparisons of computational time for the steady flow simulations, between the present IGKS and explicit GKS on uniform and/or unstructured hybrid mesh, have been carried out, indicating that the present IGKS with unstructured hybrid mesh have the best computational efficiency for steady state solutions to get the convergent flow field.

Through this study, the IGKS coupled with SA turbulence model has demonstrated to be fully capable of simulating the turbulent flows on unstructured hybrid mesh, as indicated by various comparison tests. The present IGKS solver can be extended to the 3-D case by following a parallel algorithm of LU-SGS in 3D problems reported by Gong et al (2016) (gong2016an, ) and considering the neighborhood connectivity of cells carefully, and the relevant work will be reported in the future. As such, the present IGKS solver developed in this study looks promising in its extensive applications.

## Acknowledgements

The project has been financially supported by the National Natural Science Foundation of China (Grant No. 11472219), the Natural Science Basic Research Plan in Shaanxi Province of China (Program No. 2015JM1002), the 111 Project of China (B17037), the ATCFD Project (2015-F-016), as well as the National Pre-Research Foundation of China.

## References

- (1) K. Xu, K. H. Prendergast, Numerical Navier-Stokes solutions from gas kinetic theory, Journal of Computational Physics 114 (1) (1994) 9–17.
- (2) K. Xu, A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method, Journal of Computational Physics 171 (1) (2001) 289–335.
- (3) G. May, B. Srinivasan, A. Jameson, An improved gas-kinetic BGK finite-volume method for three-dimensional transonic flow, Journal of Computational Physics 220 (2) (2007) 856–878.
- (4) M. Ilgaz, I. Tuncer, Parallel implementation of gas-kinetic BGK scheme on unstructured hybrid grids, AIAA Paper 2006-3919.
- (5) G. Kumar, S. S. Girimaji, J. Kerimo, WENO-enhanced gas-kinetic scheme for direct simulations of compressible transition and turbulence, Journal of Computational Physics 234 (2013) 499–523.
- (6) O. J. Chit, A. A. Omar, W. Asrar, M. M. Hamdan, Implicit gas-kinetic Bhatnagar-Gross-Krook scheme for compressible flows, AIAA Journal 42 (7) (2004) 1293–1301.
- (7) W. Li, M. Kaneda, K. Suga, An implicit gas kinetic BGK scheme for high temperature equilibrium gas flows on unstructured meshes, Computers & Fluids 93 (2014) 100–106.
- (8) J. Li, C. Zhong, Y. Wang, C. Zhuo, Implementation of dual time-stepping strategy of the gas-kinetic scheme for unsteady flow simulations, Physical Review E 95 (2017) 053307.
- (9) D. Mavriplis, Unstructured grid techniques, Annual Review of Fluid Mechanics 29 (1) (1997) 473–514.
- (10) Y. Liu, W. Zhang, Y. Jiang, Z. Ye, A high-order finite volume method on unstructured grids using RBF reconstruction, Computers & Mathematics with Applications 72 (4) (2016) 1096–1117.
- (11) Y. Liu, W. Zhang, Accuracy preserving limiter for the high-order finite volume method on unstructured grids, Computers & Fluids 149 (2017) 88–99.
- (12) M. S. Ebeida, E. Mestreau, Y. Zhang, S. Dey, Mesh insertion of hybrid meshes, in: Proceedings of the 18th International Meshing Roundtable, Springer, 2009, pp. 359–375.
- (13) D. Pan, C. Zhong, J. Li, C. Zhuo, A gas-kinetic scheme for the simulation of turbulent flows on unstructured meshes, International Journal for Numerical Methods in Fluids 82 (2016) 748–769.
- (14) P. Spalart, S. Allmaras, A one equation turbulence model for aerodynamic flows, AIAA Journal 94 (1) (1992) 5–21.
- (15) P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM Journal on Numerical Analysis 21 (5) (1984) 995–1011.
- (16) D. L. De Zeeuw, A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations, Ph.D. thesis, Citeseer (1993).
- (17) R. Chen, Z. Wang, Fast, block lower-upper symmetric Gauss-Seidel scheme for arbitrary grids, AIAA Journal 38 (12) (2000) 2238–2245.
- (18) M. Tidriri, Preconditioning techniques for the Newton-Krylov solution of compressible flows, Journal of Computational Physics 132 (1) (1997) 51–61.
- (19) D. Moro, N. Nguyen, J. Peraire, Navier–Stokes solution using hybridizable discontinuous Galerkin methods, AIAA paper 3407 (2011) 2011.
- (20) N. V. Shende, Y. Mor-Yossef, Robust implementation of the Spalart-Allmaras turbulence model for unstructured grid, in: 5Th European Conference on Computational Fluid Dynamics ECCOMAS CFD, 2010, pp. 1–14.
- (21) P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Physical Review 94 (3) (1954) 511.
- (22) S. Xiong, C. Zhong, C. Zhuo, K. Li, X. Chen, J. Cao, Numerical simulation of compressible turbulent flow via improved gas-kinetic BGK scheme, International Journal for Numerical Methods in Fluids 67 (12) (2011) 1833–1847.
- (23) J.-R. Carlson, Inflow/outflow boundary conditions with application to FUN3D, National Aeronautics and Space Administration, Langley Research Center, 2011.
- (24) U. Ghia, K. N. Ghia, C. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, Journal of Computational Physics 48 (3) (1982) 387–411.
- (25) C. Zhuo, C. Zhong, J. Cao, Filter-matrix lattice Boltzmann model for incompressible thermal flows, Physical Review E 85 (4) (2012) 046703.
- (26) S. Hou, Q. Zou, S. Chen, G. Doolen, A. C. Cogley, Simulation of cavity flow by the lattice Boltzmann method, Journal of Computational Physics 118 (2) (1995) 329–347.
- (27) Y. Zhu, C. Zhong, K. Xu, Implicit unified gas-kinetic scheme for steady state solutions in all flow regimes, Journal of Computational Physics 315 (2016) 16–38.
- (28) S. R. Allmaras, F. T. Johnson, Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model, in: Seventh International Conference on Computational Fluid Dynamics (ICCFD7), 2012, pp. 1–11.
- (29) C. Rumsey, Wilcox2006 expected results - 2D zero pressure gradient flat plate, http://turbmodels.larc.nasa.gov/flatplate_w06.html.
- (30) J. H. Morrison, Numerical study of turbulence model predictions for the MD 30P/30N and NHLP-2D three-element highlift configurations, National Aeronautics and Space Administration, Langley Research Center, 1998.
- (31) C. Gong, W. Bao, J. Liu, G. Tang, Y. Jiang, An efficient wavefront parallel algorithm for structured three dimensional LU-SGS, Computers & Fluids 134 (2016) 23–30.

Primary vortex | Left secondary vortex | Right secondary vortex | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

Structured grid | 2.0861 | -0.1178 | 0.5335 | 0.5681 | 2.2461 | 0.0820 | 0.0716 | 1.7123 | 0.8720 | 0.1132 |

Error (structured grid) | 1.7% | 0.1% | 0.4% | 1.0% | 2.8% | 4.5% | 8.3% | 2.2% | 1.5% | 3.5% |

Unstructured mesh | 2.0902 | -0.1179 | 0.5405 | 0.5721 | 2.2870 | 0.0899 | 0.0862 | 1.7289 | 0.8752 | 0.1145 |

Error (Unstructured grid) | 2.0% | 0.0% | 1.7% | 1.7% | 1.1% | 4.7% | 10.4% | 1.3% | 1.8% | 4.7% |

Hou et al hou1995simulation () | 2.0760 | -0.1178 | 0.5333 | 0.5647 | 2.2200 | 0.0902 | 0.0784 | 1.6900 | 0.8667 | 0.1137 |

Zhuo et al zhuo2012filter () | 2.0570 | -0.1179 | 0.5311 | 0.5662 | 2.2667 | 0.0828 | 0.0770 | 1.7066 | 0.8645 | 0.1120 |

Ghia et al ghia1982high () | 2.0497 | -0.1179 | 0.5313 | 0.5625 | 2.3113 | 0.0859 | 0.0781 | 1.7510 | 0.8594 | 0.1094 |

Number of cells | Explicit method | Implicit method | Speedup | |||
---|---|---|---|---|---|---|

Iteration steps | Time(s) | Iteration steps | Time(s) | |||

Cavity (structured grid) | 10000 | 608000 | 34820.2 | 4430 | 543.2 | 64.1 |

Cavity (unstructured mesh) | 8933 | 590000 | 30184.0 | 4270 | 484.4 | 62.3 |

Laminar flat plate | 13441 | 605000 | 64685.6 | 4320 | 717.4 | 90.2 |

Turbulent flat plate | 26386 | 810000 | 78599.5 | 5340 | 1087.7 | 72.3 |

NHLP multi-element airfoil | 64344 | 520000 | 123614.2 | 3650 | 1308.4 | 94.5 |

(a) Hybrid mesh for incompressible laminar flow, (b) Mesh near wall.

(a) CFL number, (b) Residual.

(a) Hybrid mesh, (b) Mesh near wall.

(a) , (b) .

(a) CFL number, (b) Residual.

(a) Between slat and airfoil, (b) Between airfoil and flap.

(a) CFL number, (b) Residual.