Rarefied gas flow simulation is an important field of computational fluid dynamics (CFD). Recent years, with the development of aerospace technology and micro-electromechanical system (MEMS), the rarefied flow simulation has attracted more and more attention. The discrete velocity method (DVM) [1, 2, 3, 4, 5] , or also known as the discrete ordinate method (DOM), is a classical numerical method playing an important role and widely used in the area of rarefied flow simulation. Unlike another important method for rarefied flow simulation called direct simulation Monte Carlo (DSMC)  which is based on the particle statistics, DVM is a deterministic method solving the Boltzmann equation or its model equations. In the conventional DVM, the operator splitting is used, the transportation term and the collision term of the governing equation are handled separately, which endows the scheme with a series of advantages like conciseness, easiness in programming, high efficiency and good accuracy in the simulation of high Knudsen (Kn) number. These assets just ensure the position of the conventional DVM in the area of rarefied flow simulation. However, on the other hand, due to the decoupling of the transportation term and the collision term, during the transportation process of the conventional DVM, particles are always transferring freely without any collision, causing an overly nonequilibrium particle velocity distribution when the numerical scale (i.e. the mesh size and the time step) is much larger than the kinetic scale (i.e. the mean free path and the mean collision time). This means that in the simulation of low Kn number where intensive particle collision occurs, the conventional DVM needs to adopt very small mesh size and time step, which requires huge amount of computation, otherwise the free-transport mechanism of the decoupled transportation term will make the result very dissipating. Thus, the conventional DVM can be hardly used to do fine simulations of flow in continuum regime.
The limitation of the conventional DVM makes it mainly applicable to single scale simulations where the numerical scale is comparable with the kinetic scale. However, nowadays in the field of CFD there are more and more important problems involving the continuum flow and the rarefied flow simultaneously in one flow field. These problems include but are not limited to the flow around the hypersonic vehicle which needs to travel between ground and near space atmosphere, or the flow in the reaction control system (RCS) of the spacecraft. The single-scale conventional DVM is not proper to deal with these multi-regime problems and the multiscale numerical method is required. On this point, Xu and Huang proposed the unified gas-kinetic scheme (UGKS)  for gas flows in all flow regimes. UGKS is presented in the finite volume framework, and is also solving the kinetic equation based on the discrete particle velocity space just as the conventional DVM. What distinguishes UGKS from the conventional DVM is that UGKS couples the particle transportation with the particle collision by the analytical solution of the governing equation in the calculation of the flux at the cell interface, which means that in the transportation process of UGKS the particle is not transferring freely as what happens in the conventional DVM but transferring simultaneously with particle collision. Then when the cell Kn number (the ratio of the mean free path and the mesh size) is very low, namely the numerical scale is much larger than the kinetic scale, the velocity distribution can evolve to near-equilibrium state in a numerical time step and the scheme will yield results consistent with the Navier-Stokes (NS) equation. Hence, by considering the particle collision in the calculation of the interface flux, UGKS turns into a multiscale scheme which can adopt mesh size and time step comparable with that of the traditional NS-equation-based CFD method in the continuum flow simulation, and it is very suitable for multiscale problems  involving both continuum and rarefied flows. The discrete unified gas-kinetic scheme (DUGKS) proposed by Guo et al. [9, 10] is another multiscale scheme based on the similar idea of UGKS. Instead of using the analytical solution, in DUGKS a difference scheme of the governing equation is applied at the cell interface to calculate the multiscale flux which couples the particle transportation and the particle collision. Generally speaking UGKS and DUGKS have similar multiscale mechanisms and can yield similar results in all flow regimes.
Recently, Su et al. 
put forward a general synthetic iteration scheme (GSIS) to find steady-state solutions of the kinetic equation. In their method the discrete velocity framework is adopted and the high-order terms obtained from the moments of the discrete velocity distribution function are extracted into the macroscopic governing equations which contain Newton’s law for stress and Fourier’s law for heat conduction explicitly. The method can quickly get the steady state solutions of gas flow in all flow regimes, and they explained the multiscale mechanism as that the macroscopic equations help to adjust the solution obtained from the microscopic discrete velocity system so that the constraint of the cell Kn number is removed just as UGKS and DUGKS.
Here, inspired by the work of Su et al. , we present another idea to achieve the multiscale property in the scheme adopting DVM framework. The present method takes the finite volume DVM framework, and just similar to UGKS and DUGKS the flux at the interface is specially treated. Unlike UGKS and DUGKS, the flux at the interface is not calculated by the analytical solution nor the difference scheme of the governing equation, but by the reconstruction of the moments obtained from the discrete velocity distribution function at the cell center, which is relatively concise and has some similarity in multiscale mechanism with the method of Su et al.  although macroscopic equations containing the NS constitutive relations are not used in the present work to accelerate the convergence. The method is applicable to both steady and unsteady flows in all flow regimes, and the accuracy has been verified by test cases.
The remainder of the paper is organized as follows. In Section 2, firstly a DVM-based moment method which works well in the continuum regime is introduced, then the multiscale DVM is presented by modifying the DVM-based moment method in a simple way. In Section 3, the accuracy and the multiscale property of the present method are testified by several test cases. In Section 4 a summary about the present work is given.
2 Numerical method
In this work, the monatomic gas is considered and the BGK-type kinetic model equation  is solved, which has the form
where is the gas particle velocity distribution function and is the particle velocity. is the collision time calculated as
where and are the dynamical viscosity coefficient and the pressure respectively. The equilibrium state has a form of Maxwellian distribution,
or if the Shakhov model  is used,
where is the peculiar velocity and is the macroscopic gas velocity, is the heat flux, is a variable related to the temperature by . Macroscopic variables can be obtained by taking moments of ,
is the vector of the conservative variables,is the vector of moments , is the velocity space element,
is the stress tensor. The collision term on the right-hand side of the model equation (1) follows the conservation constraint,
The present method is based on the finite volume framework and the integral form of the governing equation (1) is concentrated,
where is the control volume, is the volume element, is the surface area element and is the outward normal unit vector. Taking moments of to Eq. (9) will yield the corresponding macroscopic governing equation,
where the flux is calculated as
This paper focuses on the numerical scheme for Eq. (9) and Eq. (10). In the following, the paper will first introduce a DVM-based moment method about Eq. (9) and Eq. (10), and then modify it to a multiscale scheme based on the DVM framework. All data reconstructions involved are done by the linear reconstruction based on least square method and the Venkatakrishnan limiter  is used. Second-order accuracy in space and first-order accuracy in time are achieved in the current work. Accordingly, the conventional DVM mentioned below refers to the method with second-order accuracy in space and first-order accuracy in time.
2.1 Scheme I: DVM-based moment method
In this scheme, just like the conventional DVM, at the th time step the distribution function at the cell center is discretized both in physical space (denoted by the subscript ) and velocity space (denoted by the subscript ). The main idea of this scheme is that, solving the microscopic governing equation (9) to obtain the non-equilibrium moments and meanwhile solving the macroscopic governing equation (10) to update the macroscopic variables.
2.1.1 General construction
We here start from the update of the macroscopic variables and the macroscopic governing equation (10) is discretized as
where denotes the neighboring cell of cell and is the set of all of the neighbors of cell . The subscript denotes the variable at the interface between cell and , and is the interface area. The marching time step is restricted by the CFL condition based on the particle velocity , and is calculated as
where is the Heaviside function defined as
The CFL number is set as in the work of the present paper. In Eq. (12) the flux at the interface needs to be carefully treated to avoid excess numerical dissipation in the case of low cell Kn number. For the conventional DVM of first-order temporal accuracy, this flux is calculated as the moments of the distribution function at the interface,
where the interface distribution function is obtained directly from the reconstruction of the initial distribution function data. Such a treatment will lead to more dissipating solution in the low-cell-Kn-number case because the reconstructed distribution function at the interface suffers from a large deviation, which is of the order of
due to the interpolation error for second-order accuracy, from the equilibrium state. In the real physical situation, the deviation should be of the order of the collision timewhen the kinetic scale is much smaller than the flow characteristic scale, just as what is revealed by the Chapman-Enskog expansion . Larger deviation means additional viscous stress and heat flux, making the conventional DVM very dissipating when the cell size is much larger than the mean free path (see Fig. 8(d) in Section 3.3 for how this excess numerical dissipation influences the results). In UGKS [7, 8] and DUGKS [9, 10], the problem is solved by evolving the interface distribution function with a CFL time step from the initial reconstructed data based on the governing equation (1). After the evolution the interface distribution function will be sufficiently close to the equilibrium state in the low-cell-Kn-number case, thus there is no additional dissipation caused by the above problem. Here, inspired by the work of Su et al. , the paper provides another idea. Considering that the distribution function stored at the cell center can be obtained directly and is free from the data reconstruction error, and should be sufficiently close to the equilibrium state due to the collision process of BGK equation (1) in the low-cell-Kn-number case, it is reasonable to think that the flux calculated from the moments of the cell-center distribution function is reliable, i.e.
where is the flux tensor at the cell center. On the other side, directly interpolating the flux at the cell interface from the cell-center flux data may cause stability issue. Hence, here the flux is divided into the equilibrium flux G and the nonequilibrium flux H, i.e.
Then at the interface, the flux is also split as
where and will be calculated in different ways. Such a treatment is similar to the separate handling of convection term and diffusion term in the traditional NS-equation-based scheme.
For the nonequilibrium flux tensor H at the cell center, it is calculated as
where is the equilibrium state determined from the macroscopic variables which are obtained from the cell-center distribution function ,
The second part on the right-hand side of Eq. (20) is directly the Euler flux tensor determined by . After that, through the reconstruction of the nonequilibrium flux tensor H, at the interface the nonequilibrium flux is calculated as
For the equilibrium flux, we don’t calculate the tensor but directly calculate the flux at the cell interface. The equilibrium flux is actually the Euler flux which should be determined from the macroscopic variables at the interface. Here is calculated by the way used in gas-kinetic scheme (GKS) . First, the macroscopic variables are reconstructed and the variables can be interpolated at the two sides of the interface as
Then the macroscopic variables for the interface are calculated as
where and are determined by and respectively. After that, to make the scheme stable, the equilibrium flux at the interface is calculated as a weighting of the Euler flux and the flux of kinetic flux vector splitting (KFVS)  which is of high stability, i.e.
where the first part on the right-hand side is corresponding to the KFVS flux (for the calculation about the moments one can see the appendix of Ref.  for a quick guide) and the second part is the Euler flux determined from . The weight factors and are constructed based on the form of DUGKS [9, 10]. is the collision time determined by . is the physical local time step determined by the CFL condition based on the macroscopic characteristic velocity, i.e.
where is the acoustic velocity and is the maximum cross section area of cell . For , it is set as in the work of the present paper. Decreasing will increase the weight of the KFVS flux, resulting in better stability but more numerical dissipation. can be viewed as the characteristic time scale based on the cell size for the interface . According to the construction of the weight factors in Eq. (27), in the low-cell-Kn-number case (continuum case) the flux can recover the accurate Euler flux while in the high-cell-Kn-number case (rarefied case) it can recover the KFVS flux to stabilize the scheme. For more details about these weight factors one can refer to Ref.  and Ref. .
Now that we have obtained both the nonequilibrium flux and the equilibrium flux , the flux is then calculated by Eq. (19) and the macroscopic variables can be updated to through Eq. (12). The last thing to do is the update of the distribution function at the cell center. The microscopic governing equation (9) is discretized as
where the collision term on the right-hand side adopts an implicit form to ensure the stability, and are both calculated by the updated macroscopic variables . Eq. (30) is further arranged as
The transportation term adopts the same treatment as in the conventional DVM, and the interface distribution function is calculated through the reconstruction of the distribution function as
The distribution function at the cell center can be then updated to by Eq. (31).
2.1.2 Boundary condition
The treatment of the boundary is one notable thing for the scheme. Here we only discuss the implementation of the diffuse reflection boundary condition with full thermal accommodation  on the isothermal solid wall, while for other boundary conditions one can handle them in the conventional manner. On the full-diffuse-wall boundary, just as at the interior face, the microscopic flux (flux for distribution function) and the macroscopic flux are treated separately. For the microscopic flux, the treatment is just the same as in the conventional DVM and the distribution function on the boundary face is
where is the center of the boundary face, is the outward normal unit vector relative to the adjacent cell . is the Maxwellian distribution for the outcoming reflected particles and is calculated as
in which is the wall velocity and is obtained from the wall temperature. can be solved from the no-penetration condition
Then on the boundary face is completely determined and can be used to calculate the microscopic flux in Eq. (31) to update the distribution function. For the macroscopic flux which can be calculated as
the incoming flux to the wall is divided into the equilibrium flux and the nonequilibrium flux just as what we have done at the interior face to avoid excess numerical dissipation in the low-cell-Kn-number case, i.e.
The nonequilibrium flux is calculated directly from the variables at the center of the adjacent cell as
which means the first-order accuracy for this part of flux. The equilibrium flux is calculated by
where is determined by the macroscopic variables which are obtained through reconstruction with second-order accuracy as
Such a construction makes the main part of second-order accurate in the continuum flow case. After is obtained, the outcoming flux is calculated as
where the reflected Maxwellian distribution is
which is different from the discrete obtained by Eq. (34) because here should be constrained by the no-penetration condition for the macroscopic flux
where and are the incoming mass flux and the outcoming mass flux in and respectively. The macroscopic flux is then completely determined and can be used in the update of the macroscopic variables Eq. (12).
2.1.3 Artificial dissipation
When applying the above scheme to the supersonic flow simulation with high Re number, artificial dissipation should be added to suppress the oscillation induced by the shock wave. The artificial dissipation is introduced by the following two measures. First, the weight of the KFVS flux for in Eq. (27) is increased at the discontinuity region by calculating the collision time as
where and are pressures calculated from the reconstructed macroscopic variables and at the two sides of the interface. By applying Eq. (44), the pressure discontinuity will increase the weight of the KFVS flux and the oscillation induced by such discontinuity can be suppressed. Second, an artificial viscous term is introduced explicitly into the interface flux by amplifying the nonequilibrium flux , that is
where the factor is used to switch off the artificial viscosity in the case of rarefied flow.
At the end of this part, we make some remarks about the above scheme. By interpolating the nonequilibrium flux from the cell center to the cell interface, the scheme avoids the excessive dissipation suffered by the conventional DVM in the continuum flow simulation. The scheme adopts the DVM framework, but the microscopic flux (the flux of the distribution function) in Eq. (31) and the macroscopic flux calculated by Eq. (19) or Eq. (45) are not consistent, so the macroscopic variables are not equal to the macroscopic variables calculated by Eq. (21) from the integration of the distribution function at the cell center. The discrete distribution function is only used to obtain the nonequilibrium flux and the method is essentially a moment method. In the rarefied flow simulation with a high cell Kn number, the relaxation (or collision) process of the microscopic governing equation drives the distribution function to the equilibrium state determined by the macroscopic variables very slowly (see Eq. (31) for reference), making the distribution function and the nonequilibrium flux nearly decoupled with the macroscopic variables . Due to this decoupling, in the rarefied flow simulation, the macroscopic variables may present a fake result which is spoiled by the discretization error. Our test cases show that the scheme I works well in the continuum flow cases but works improperly in the rarefied flow cases. One recipe (not applied in the current work) for this problem is to view the macroscopic variables obtained from by Eq. (21) as the real solution of the scheme. This strategy is useful in the single-scale single-regime simulation, but may not work well in the multiscale multi-regime simulation because the unreliable variables in the rarefied region can still spoil the solution of the adjacent (either in time or in space) continuum region.
2.2 Scheme II: Multiscale DVM
In the scheme I, the macroscopic variables and the distribution function get decoupled when the cell Kn number is high and the scheme becomes unreliable in such a case. Here a very concise treatment can be used to fix this problem. Just as mentioned in the end of Section 2.1, the macroscopic variables obtained from the integration of are more reliable. Therefore, after the update of the distribution function through Eq. (31), the macroscopic variables at the next time step are recalculated by integrating as
This means that the macroscopic variable updated through the macroscopic flux by Eq. (12) is only an intermediate to get the equilibrium state in Eq. (31), and to avoid ambiguousness these equations are rewritten in the scheme II as
where are the intermediate macroscopic variables and is the intermediate equilibrium state determined by . By simply applying Eq. (46), the macroscopic variables and the distribution function are closely coupled and the scheme turns into a multiscale scheme for all flow regimes. Further investigation for the mechanism will ensue. Combining Eq. (46), Eq. (47) and Eq. (48) will yield
and it can be arranged as
Eq. (50) shows that the macroscopic variables are intrinsically updated through the weighting of and , where the former is the conventional DVM flux calculated by Eq. (16) and the latter is the flux calculated by Eq. (45) which is of relatively low numerical dissipation in the low-cell-Kn-number case. Such a construction is quite similar to other multiscale methods like UGKS [7, 8] and DUGKS [9, 10]. In those methods the flux at the cell interface couples the particle free transportation with the particle collision through the BGK equation, and can be also approximately viewed as the weighting of the conventional DVM flux and the flux based on NS equation. Furthermore, the weight factors for the fluxes in Eq. (50) are very similar to those in DUGKS [9, 10], which are deduced from the discretization of the BGK equation likewise. Thus, the scheme applying Eq. (46) can get the multiscale property based on the same mechanism as UGKS and DUGKS.
It is worth noting that, the scheme using Eq. (46) works well in the uniform mesh, but in the non-uniform mesh due to the constraint of the global CFL condition the marching time step can be very small and the weight for the conventional DVM flux will be overestimated as shown in Eq. (50), leading to excessive dissipation. In view of this, in the scheme II the macroscopic variables are not directly updated by Eq. (46), but by modifying Eq. (50) as
which can be further simplified as
where is the physical local time step obtained from Eq. (28). In addition, considering that in Eq. (51) (and also in Eq. (50)) there is already a weighting process which introduces the dissipative conventional DVM flux into the flux at the interface, the calculation of the equilibrium flux (i.e. Eq. (27)) is modified in the scheme II as
To make the scheme II more clear, the computation procedure in one step is listed as follows:
- Step 1.
- Step 2.
Calculate the equilibrium flux at the cell interface by Eq. (53) with data reconstruction.
- Step 3.
- Step 4.
Calculate the intermediate macroscopic variables by Eq. (47).
- Step 5.
Update the distribution function by Eq. (48), during which the reconstruction of the distribution function is implemented.
- Step 6.
Update the macroscopic variables by Eq. (52).
Finally, we’d like to point out that in Eq. (51), because the weight factors for the fluxes are calculated by the cell-related variables and , the actual macroscopic flux across the cell interface will have different values for the two cells on the two sides, which will cause conservation issue. Fortunately, the fluxes and satisfy the conservation constraint separately (i.e. the same for cells on both sides), so this conservation issue is not so serious. For simulations which have initial-value dependence, some simple tricks (such as the mass compensation tricks used in  and ) can be applied to fix the problem. Also, one can further modify Eq. (51) and calculate the weight factors using the interface-related variables such as and , which solves this issue completely with increased computational complexity and computational cost. In the present work, we directly apply Eq. (52) without any additional tricks, and the numerical results seem not spoiled by the conservation issue.
3 Numerical results
To investigate the performance of the scheme II (the multiscale DVM) as well as the scheme I, three sets of test cases are carried out in this section. The cases include the Sod’s shock tube, the lid-driven cavity flow and the laminar boundary layer flow over a flat plate, covering a variety of conditions such as unsteady, steady, continuum, rarefied. In all of the test cases the working gas is monatomic and the Shakhov model Eq. (4) with a Prandtl number equal to is used.
3.1 Sod’s shock tube
The 1D test case of Sod’s shock tube is performed to show the capability of the schemes in unsteady flow simulation for various flow regimes. The computational domain is with a length and the initial field is
with a jump at the center of the field. The computational domain is discretized into 100 uniform cells. For the particle velocity space, the discrete domain is set as and it is also discretized uniformly into 100 cells. The hard sphere (HS) model is used with heat index . Three Kn numbers, which are defined by the length of the physical-space domain , are considered, i.e. , and . Because of the initial discontinuity in the field, the Venkatakrishnan limiter  is used. The time step is set as and the result at is investigated.
The results are shown in Fig. 1, Fig. 2 and Fig. 3. The reference results are calculated by UGKS [7, 8] under the same simulation setup. For the case in the continuum regime, the three sets of results agree well, and the results obtained from the scheme I and the scheme II completely overlap. The results all suffer from some oscillations at the position of the initial discontinuity. These oscillations can be suppressed by varieties of measures, such as decreasing in Eq. (28), adding more artificial dissipation, modifying the reconstruction method, but here we don’t pay much effort on this. For the case , the three sets of results are also in good consistence. But comparing with the results of , slight differences exist at the position of the shock wave between the results from scheme I and the results from scheme II. For the case , the results from UGKS and scheme II agree well while the results from scheme I have some observable errors, which implies the invalidation of the scheme I in such rarefied case just as discussed in Section 2.1.4. Thus, it can be seen that the multiscale DVM (the scheme II) works well in this set of test cases for various flow regimes, while the scheme I is only valid for the continuum flow regime.
3.2 Lid-driven cavity flow
The 2D test case of lid-driven cavity flow is conducted to research whether the schemes can accurately simulate the viscous effect in various flow regimes. In this test case the gas flow inside a closed square cavity is driven by the tangentially moving lid. The Mach number defined by the lid velocity and the acoustic velocity is around , where is defined by the fixed wall temperature . Aiming to investigate the property of the schemes from continuum to rarefied regimes, three degrees of viscosity (or degrees of rarefaction), including and , are considered, where the reference length is the cavity width . The HS molecular model is applied and the heat index is . The physical-space computational domain is discretized into structured mesh. For the case a nonuniform mesh is used with a minimum mesh size near the walls (as shown in Fig. 4(a)), and for the cases a uniform mesh is used. The discrete velocity-space domain is set as a square , and is discretized into , , uniform cells for the cases and respectively. On the walls of the cavity, the diffuse reflection boundary condition with full thermal accommodation is imposed and the implementation of this boundary condition is discussed in Section 2.1.2.
The results are shown in Fig. 4, Fig. 5 and Fig. 6. For the continuum case , the results from scheme I and scheme II overlap with each other, and they agree well with the results obtained from GKS , which can yield NS solutions in second-order accuracy, under the same simulation setup. For the case , generally speaking the three sets of results agree well. More precisely, the results from scheme II coincide with the results from UGKS perfectly, while the results from scheme I agree well with other two sets of results except some very minor and insignificant mismatches near the wall boundary. For the case , the results from scheme II and UGKS are in very good consistence while the results of scheme I suffer some small deviations, which are larger than the mismatches in the case , near the wall. In conclusion, for this set of test cases the multiscale DVM (the scheme II) can give results identical to the results of GKS in the continuum case and the results of UGKS in rarefied cases, meanwhile the scheme I can also give acceptable results except some small deviations in the rarefied cases.
3.3 Laminar boundary layer flow over a flat plate
For a multiscale scheme, besides describing the dynamics of particle transportation and collision when the numerical cell size is comparable with the kinetic scale, it should recover the accurate viscous hydrodynamic mechanism when the numerical scale is much larger than the kinetic scale. Here the test case of laminar boundary layer flow over a flat plate is conducted to test if the schemes can recover the accurate constitutive relations in the continuum flow regime just as an NS-equation-based scheme without excess numerical dissipation under a cell size much larger than the kinetic scale. The freestream condition is , where is the length of the flat plate. As shown in Fig. 7, the physical-space computational domain has a length of in and a height of to diminish the influence of the domain boundary, and the leading edge of the plate is placed at . In -direction the computational domain is discretized into cells, where cells are in the range along the flat plate. The increasing rates of the cell size from the leading edge of the plate are around upstream and downstream, with the minimum cell width near the leading edge. In -direction, three different resolutions are adopted and the minimum cell heights are , with the same increasing rate around from the surface of the plate to the far field and the discretization numbers are , , respectively. The setup of the boundary condition is shown in Fig. 7, where the top and left boundaries of the domain are set as the freestream condition, the right boundary is set as the outlet condition, on the bottom boundary during the range ahead of the plate the symmetric condition is applied, and it is notable that the no-slip bounce-back boundary condition (implemented by simply mirroring the variables into the ghost cells), but not the diffuse reflection boundary condition discussed in Section 2.1.2, is imposed on the surface of the plate. For the discrete velocity-space domain, the range is where is the freestream acoustic velocity, and it is discretized into uniform cells.
For purpose of comparison, the simulations under the same conditions have also been performed by GKS , which is a numerical scheme working in the hydrodynamic (Navier-Stokes) scale, and the conventional DVM (more exactly, the conventional DVM of first-order temporal accuracy discussed in Section 2), which is a numerical scheme working in the kinetic (mean-free-path) scale. The results of the velocity profiles at obtained from the four methods are shown in Fig. 8. It is shown that for the scheme I, the scheme II and GKS, they can give roughly accurate velocity profiles even under the coarsest mesh where the minimum mesh height is and there are only cells in the boundary layer. The deviations of the vertical velocity profiles are relatively larger under such a coarse resolution. Besides, the vertical velocity profiles obtained from these three methods under all resolutions fall below the analytic solution in the far field, which may due to the singularity of the flow at the leading edge because the cross section is close to there. Generally speaking the three methods are in the same level of accuracy under all of the mesh resolutions for this test case. In contrast, for the conventional DVM, it fails to give an acceptable result even under the finest mesh, and all of its results deviate largely from the analytical solution, showing a much thicker boundary layer and suggesting excessive numerical dissipation.