A discrete unified gas kinetic scheme on unstructured grids for viscid compressible flows and its parallel algorithm

1.   Introduction

Compressible flows exist in many practical engineering, such as aerophysics, astrophysics, aeroacoustics, combustion, explosion physics, etc, which is generally illustrated by the Navier-Stokes equations (NSEs) ^[1,2,3,4]. In the conventional computational fluid dynamics (CFD), the discrete method, for example, finite difference (FD) ^[5], finite volume (FV) ^[6], finite element methods ^[7,8] are employed to discretize the NSEs. Instead of solving compressible NSEs, the kinetic schemes, such as the lattice Boltzmann method (LBM) ^[9], the gas kinetic scheme (GKS) ^[10], the lattice Boltzmann flux solver ^[11], the discrete velocity or discrete ordinate methods ^[12] and the (discrete) unified gas kinetic scheme ^[13,14], are derived from the Boltzmann equation to evolve the gas distribution function and calculate the conservative variables by the gas distribution function.

The conventional LBM is a special FD discretization of the Boltzmann equation, in which the fluid domain is discretized into lattices coinciding with the uniform Cartesian grids, and the motion of particles is described by the evolution of particle distribution functions named as "collision" and "stream". The macroscopic variables are calculated by accumulating the particle distribution functions. It has been an alternative method to simulate low-speed incompressible flows. Owing to the restriction of uniform Cartesian grid, the conventional LBM is difficult for simulating the flow problems with strong local gradients or defined on geometrically complex domains. In recent years, a few off-lattice Boltzmann methods were proposed to avoid the Cartesian grid restriction, for example, finite difference LBMs ^[15,16,17], finite element LBMs ^[18,19], FV-LBMs ^[20,21,22], and so on. In the conventional LBM, the equilibrium distribution function (EDF) is obtained by expanded the Maxwellian function by Taylor-series in terms of the Mach number, which can be only used for low Mach number incompressible flows. To simulate high Mach number compressible flows, an alternative method is to replace the Maxwellian function with a simple function avoiding the taylor-series expansion. Sun et al. proposed a locally adaptive LBM, in whic a simplified EDF ($\delta$ function) is used to substitute the Maxwellian function ^[23]. A D2Q13L2 lattice model for inviscid high-speed compressible flows is presented by Qu et al. Instead of the Maxwellian function, a circular function is employed. The EDF is constructred by a polynomial and a FV method with second-order total variation diminishing formular is employed to solve the lattice Boltzmann equation ^[24]. A LBM for high-speed viscid compressible flows based on D2Q17 circular function is proposed by Qu et al. ^[25,26], in which the first to fourth order constraints of the density equilibrium function in the heat flux and the energy dynamics are considered ^[27]. Based on previous works, we presented a FV-LBM on unstructured grids for high Mach number compressible flows ^[22,28].

Recently, the discrete unified gas kinetic scheme (DUGKS), a FV-based mesoscopic CFD method, has attracted more and more attention. The numerical flux in the DUGKS is constructed by the distribution function at a half time-step, in which a numerical characteristic solution of the Boltzmann equation is used. It has been widely studied in multiscale flows ^[29,30], thermal incompressible flows ^[31,32], thermal compressible flows ^[33], moving boundaries ^[34], gas mixture flows ^[35], and nonlinear convection diffusion equations ^[36]. Lately, a simplified DUGKS for incompressible flows is proposed, in which a whole time-step is applied to constructed the numerical flux at the cell interface ^[37], which is applied in this paper. Subsequently, they devised a simplified DUGKS for compressible flows ^[38], in which, the D2Q13 circular distribution function is used to substitute the Maxwellian function. In this paper, the D2Q17 circular function is employed instead of D2Q13 circular function in order to consider the first to fourth order constraints of distribution function.

In this paper, in order to finish the simulation in a short time, a parallel version of the proposed DUGKS on unstructured grids is implemented. In this paper, an open source package ParMETIS is used to partition the unstructured grid into irregular subdomains ^[39], which ensure the minimum of communication among subdomains and promises the load balance among processor cores. The parallel algorithm implemented by PETSc (portable, extensible toolkit for scientific computation) ^[40] is presented, and parallel performance is investigated on a supercomputer.

The remainder of the paper is organized as follows. DUGKS for viscid compressible flows as well as the parallel algorithm is presented in Section 2. Next, numerical results and the parallel performance are given in Section 3. Finally, some conclusions are made in Section 4.

2.   Parallel DUGKS for viscid compressible flows

2.1.   Kinetic model equations

The Boltzmann-BGK model for viscid compressible flows includes the density and the total energy distribution function. For the density distribution function, the discrete Boltzmann-BGK equation is written as

where $f$ represents the density distribution function with velocity $\boldsymbol \xi$ at the position $\mathbf{x}$ and time $t$, $\mu$ the viscosity, $p$ the pressure, $\tau_f = \mu/p$ the particle relaxation time, and $f^{eq}$ the local Maxwellian function expressed as ^[27,41]

where $\rho$ denotes the density, $T$ the temperature, $\mathbf{u}$ the fluid velocity, $R$ the gas constant and $D$ the spatial dimension.

To recover the compressible Navier-Stokes equations, for the total energy distribution function, the governing equation is introduced in the following form

where $h$ represents the total EDF, $\tau_h$ the total energy relaxation time, $\tau_{hf} = \tau_f\tau_h/(\tau_f-\tau_h)$, the Prandtl number $Pr$ is computed by $\tau_f/\tau_h$, the total energy EDF $h^{eq}$ is written as

where $b$ is associated with the specific-heat ratio $\gamma$ by $\gamma = (b+2)/b$. The macroscopic variables are the moments of the distribution function,

where $E = \mathbf{u}^2/2+e$ denotes the total energy, e the internal energy. $f^{eq}$ and $h^{eq}$ have to meet the following velocity moment condition to recover the compressible NSEs:

where $p = (\gamma-1)\rho e$ represents the pressure, $\delta_{ij}$, $\delta_{ik}$ and $\delta_{jk}$ the Kronecker delta functions, $i$, $j$, and $k$ the components along the $x$, $y$ and $z$ directions.

To simulate the high Mach number compressible flows, the local Maxwellian EDF $f^{eq}$ is substituted by a circular function. All the needed statistical relations in (2.3) has to be satisfied with the circular function. A D2Q13 circular function is introduced into DUGKS for high Mach number compressible flows by Zhong et al. ^[38]. In fact, $f^{eq}$ deduced from the D2Q13 circular function satisfies all the 7 equations in (2.3) except the 5th one. Note that A D2Q17 circular function is also proposed by Qu et al. in ^[26], which fulfills all the equations in (2.3). The D2Q17 circular function can be written as

where $c = 2\sqrt{RT}$ denotes the radius of the circle. To discrete the continuous EDF, a fourth-order 17-node Lagrangian interpolating polynomial called assigning function is constructed, which has the following form

Here $p = 0, 1, 2, \cdots, 16$. The corresponding discrete velocity shown in Figure 1 is expressed by

Owing to the $\delta$ property of the Lagrangian plynomial, we get

where $q = 0, 1, \cdots, 16$, $\xi_{x, q}$ and $\xi_{y, q}$ represent the $x$ and $y$ components of the particle velocity in the $q$th direction. The coefficient $a_{p, q}$ can be determined by substituting the discrete velocity into (2.6). The density EDF is then calculated by

Finally, we have

The total energy EDF $h_q^{eq}$ can be computed by

where $T_c$ denotes the characteristic temperature, $\tilde{c} = \sqrt{RT_c}$. $\varpi_0 = 0$, $\varpi_{1-4} = -17/56$, $\varpi_{5-12} = -1/8$, and $\varpi_{13-16} = -3/56$.

2.2.   DUGKS for the velocity and total energy distribution functions

we first discrete the Boltzmann-BGK equation for the density distribution function by DUGKS. Integrating (2.1) over a control volume (CV) $V_j$ located at the cell center $\mathbf{x}_j$ from time $t_n$ to $t_{n+1}$, we get the following equation ^[37]

where $F_{j, meso}^{n+1}$ represents the flux of the density distribution function across the cell interface defined as

$\partial V_j$ and $|V_j|$ are the surface and volume of CV $V_j$, $f_j$ the cell-averaged value of the density distribution function, $\mathbf{n}$ the outward unit normal of the cell interface. (2.8) is simplified as

Note that $F_{j, meso}^{n+1}$ is implicit. To compute $F_{j, meso}^{n+1}$, $f$ at the cell interface $\mathbf{x}_b$ and time $t^{n+1}$ are needed. We discrete and integrate the kinetic equations (2.1) along the particle velocity characteristic lines in one time step $\Delta t$, and employ the trapezoidal rule, we obtain

where $\bar{f}^{+} = f+\frac{\Delta t}{2}\Omega_f$, $\bar{f} = f-\frac{\Delta t}{2}\Omega_f$, $\Omega_f = -\frac{1}{\tau_f}(f-f^{eq})$. $\bar{f}^{+}(\mathbf{x}_b-\boldsymbol \xi\Delta t, \boldsymbol \xi, t_n)$ is evaluated by

where $\mathbf{x}_i$ is the cell center of $\mathbf{x}_j$'s neighbor cell. $\mathcal{L}$ the Venkatakrishnan limiter function ^[42].

DUGKS for (2.1) can also be applied similarly for (2.2). Integrating over CV $V_j$ from $t_n$ to $t_{n+1}$, one can get the following updating equation for $h_j$ for the total energy distribution function

where

is also implicit. To obtain $h(\mathbf{x}_b, \boldsymbol \xi, t)$, we integrate (2.2) within one time step in a similar way and get

where $\bar{h} = f-\frac{\Delta t}{2}\Omega_h$, $\bar{h}^{+} = h+\frac{\Delta t}{2}\Omega_h$, $\Omega_h = -\frac{1}{\tau_h}(h-h^{eq})+ \frac{1}{\tau_{hf}}(\boldsymbol \xi\cdot\mathbf{u})(f-f^{eq})$. $\bar{h}^{+}(\mathbf{x}_b-\boldsymbol \xi\Delta t, \boldsymbol \xi, t_n)$ is also calculated by

Next, the macroscopic variables at $\mathbf{x}_b$ are calculated with $\bar{f}(\mathbf{x}_b, \boldsymbol \xi, t_n+\Delta t)$ and $\bar{h}(\mathbf{x}_b, \boldsymbol \xi, t_n+\Delta t)$ in the following form

Thus, we can obtain the EDFs at the cell interface $\mathbf{x}_b$ based on the computed macroscopic variables. The calculation of $F_{j, meso}^{n+1}$ and $H_{j, meso}^{n+1}$ needs $f(\mathbf{x}_b, \boldsymbol \xi, t_{n+1})$ and $h(\mathbf{x}_b, \boldsymbol \xi, t_{n+1})$, which are got from $\bar{f}(\mathbf{x}_b, \boldsymbol \xi, t_{n+1})$ and $\bar{h}_b(\mathbf{x}, \boldsymbol \xi, t_{n+1})$ by

In (2.8) and (2.11), $f_j^{eq, n+1}$ and $h_j^{eq, n+1}$ are unknown. To solve this, the macroscopic variables at the cell center $V_j$ has to be computed. The macroscopic variables can be deduced from the moments of the Boltzmann equation, which can be written as

where

Now, the update of $f$ and $h$ can be executed in accordance with (2.8) and (2.11).

The CFL condition determines the time step size $\Delta t$ in the following form

where $C$ is the maximal discrete velocity, $\Delta x_{min}$ the minimum grid spacing, and $\eta$ the CFL number. The treatment of the boundary condition is implemented by the ghost cell method and the non-equilibrium extrapolation scheme ^[43]. More details can be found in ^[22].

2.3.   Parallel implementation of DUGKS

To imrpove the computational time, parallel algorithm is required. To do the parallel computing, the first thing is to decompose the unstructured grid into irregular subdomains. In our work, we partition the unstructured grid into irregular subdomains by an open source package ParMETIS ^[39], which ensures the minimum of the parallel communication among subdomains and promises the load balance among processor cores. An example of an unstructured grid and its partition into 16 subdomains using ParMETIS is shown in Figure 2. Each subdomain is distributed to one processor core after decomposition, and the calculations on subdomains are done concurrently. In the DUGKS. most operations are local. Only the density distribution functions and total energy distribution functions of the outmost cells in a subdomain are required by neighboring processor cores, where the parallel communication occurs. We implement the parallel algorithm based on PETSc from Argonne National Laboratory ^[40]. Further reading about the parallel strategy for the DUGKS can be found ^[20].

3.   Numerical simulations

In this section, we validate our method by two numerical siulations, including flows over the NACA0012 airfoil and flows past a circular cylinder. All computations are executed on the Tianhe-2A supercomputer at the National Supercomputer Center in Guangzhou, China. Each node consists of two 12-core Intel Ivy Bridge Xeon CPUs and 64 GB memory. The specific heat ratio $\gamma$ is $1.4$, and the Prandtl number $Pr$ $0.71$ in the simulations.

3.1.   Flows around the NACA0012 airfoil

The first benchmark test is the flow around the NACA0012 airfoil with two different configurations. For the first configuration, the Mach number at farfield $M_\infty = 0.5$, the angle of attack (AoA) $\alpha = 0^{\circ}$. The Reynolds number is set to $Re = 5000$. The Mach number at farfield for the second configuration is$M_\infty = 2.0$, AoA $\alpha = 10^{\circ}$, and $Re = 1000$. We adopt an unstructured grid with $2.87\times10^5$ quadrangular cells in these two numerial simulations. Figure 3 illustrates the computational domain and boundary configurations. A nonslip boundary condition is imposed on the airfoil surface, the non-equilibrium extrapolation method is employed at the inlet and outlet.

The Mach number contours around the airfoil for the first simulation is shown in Figure 4(a), which is inaccrodance with the Mach number contours in ^[44]. Figure 5(a) gives the comparision of the calculated pressure coefficient with the one in ^[44], which shows that the results agree well with the reference data. The Mach number contours and pressure coefficients for the second simulation are given in Figures 4(b) and 5(b), respectively, which indicate that the simulated results are in good accordance with the results in ^[45].

3.2.   Flows past a circular cylinder

In this subsection, two simulations are performed to compute supersonic flows past a circular cylinder ^[46] including (1) $M_\infty = 1.7$ and (2) $M_\infty = 2.5$. The Reynolds number is $2\times 10^5$. The computational domain is a square $[0, 1]\times[0, 1]$ on the $\xi$-$\eta$ plane. The boundary on the $x$-$y$ plane is produced by (as shown in Figure 6)

Here $R_x=4$, $R_y=10$, $a= \rm{tanh}(d)$, and $d=2.3$. We employ an unstructured grid with $32,188$ quadrangular cells in these two simulations. The inflow boundary is imposed at $\eta=0$, where the distribution functions are fixed at their equilibrium states, the outflow boundary and the symmetric boundary are used at $\xi=1$ and $\xi=0$, repsectively, where the nonequilibrium extrapolation method is applied ^[43], the nonslip boundary is imposed at the cylinder surface $\eta=1$.

The Mach number and pressure contours for $M_\infty = 1.7$ and $M_\infty = 2.5$ are illustrated in Figures 7 and 8, respectively. A bow shock located at $x = -3$ and $x = -1.9$ is formed at the upstream. The pressure coefficients $Cp$ for the two cases along the surface of the cylinder compared with the ones in ^[38,46,47] are shown in Figure 9. In ^[47], an immersed boundary method is adopted on Cartesian-grid. The results in ^[46] are computed by a FD LBM. A GKS is used to compute the pressure coefficient in ^[38], which are different from our DUGKS. It can be seen that our results agree well with the ones in the literature.

3.3.   Parallel scalability

The parallel performance of the proposed scheme is evaluated on the Tianhe-2A supercomputer ^[48]. In the test of parallel performance, a flow around the NACA0012 airfoil with two different grid sizes is used. The compute time and parallel efficiency of the scheme is shown in Table 1, in which the first 100 iterations are used to count the compute time. In the table, "Nodes" represent the number of compute nodes, "Cores" the number of processor cores, "Time" the compute time, and "Efficiency" the parallel efficiency, respectively. The speedup and parallel efficiency is illustrated in Figure 10. The larger the grid size is, the higher the speedup and parallel efficiency achieved.

4.   Conclusions

In this paper, we presented a discrete unified gas kinetic scheme on unstructured grids and its parallel algorithm to simulate high Mach number viscous compressible flows, which is validated by two benchmark fluid problems. The validation indicates that the shock waves of high Mach number viscous compressible flows can be accurately captured. The parallel performance of the scheme investigated on the Tianhe-2A supercomputer shows that the parallel algorithm has good scalability on thousands of processor cores, which is suitable to perform high-resolution simulations of high Mach number viscid compressible flows. Based on our presented method, other discrete velocity models and equilibrium distribution functions for compressible flows are also be able to incorporate with our solver, and a similar scalability can be expected. It is convenient to move to three-dimensional compressble flows, we just replace the discrete velocity model and its corresponding equilibrium distribution functions. We expect that our numerical methods can be used for solving some related inverse problems of practical importance ^[1,2]; see also ^{[49,50,51,52,53]} for more related background discussion on inverse problems.

Acknowledgments

This work was supported by the National Key R & D Program of China (Grant No. 2018YFE0198400), the Shenzhen grant (Grant Nos. RCYX20200714114735074, JCYJ20220531100611025 and JCYJ20200109115422828), and the NSFC (Grant Nos. 62161160312, 12071461 and 12101588).

Conflict of interest

The authors declare there is no conflict of interest.