Prediction Capability of Cartesian Cut-Cell Method with a Wall-Stress Model Applied to High Reynolds Number Flows

The Cartesian cut-cell method is one of the most promising methods for computational fluid dynamics due to its sharp interface treatment. However, the Cartesian cut-cell method and other Cartesian mesh solvers have difficulty with concentrating grid to boundary layers. The wall-modelling of shear stress is one of the most effective methods to reduce computational grids in boundary layers. This study investigated the applicability of a wall-stress model to the Cartesian cut-cell method. In the numerical simulations of the flow around a triangular column, Cartesian cut-cell simulation with the wall-stress model adequately predicted the drag coefficient. In the numerical simulations of the flow around a 30P30N high-lift airfoil configuration, the Cartesian cut-cell simulation with the wall-stress model adequately predicts the lift coefficient. The intermittent vortex structure of the outer layer of the turbulent boundary layer was observed on the suction side of the main element and the flap. The Cartesian cut-cell method with a wall-stress model is useful for predicting high Reynolds number flows at Re ∼ 106.


Introduction
High Reynolds number flows Re 10 6 frequently appear in the aerodynamic design of aircrafts and spacecrafts. In the computational fluid dynamics (CFD), Reynolds-averaged Navier-Stokes (RANS) simulation achieves remarkable success in steady flows. However, RANS simulation is difficult to apply to unknown unsteady flows. Large eddy simulations (LES) are rather reliable in simulations of unknown unsteady flows.
In the turbulent boundary layer of Re 10 6 flows, LES struggles to satisfy the sufficient grid resolution requirement due to the vortex scale decreasing in the inner layer of the turbulent boundary layer, which is about 10% of the turbulent boundary layer thickness. Choi and Moin [1] estimated that the total number of the computational grid required in a wall-resolved LES is proportional to Re 13/7 .
The near-wall treatment such as RANS-LES hybrid simulation [2]; detached eddy simulation (DES) [3]; and wall-stress model [4][5][6][7][8] have success with the prediction of the unsteady high Reynolds number flows. The wall-stress model using one-dimensional equations in wall-normal coordinates for turbulence boundary layer was introduced by Kawai and Larson [5,7,8]. This model yields accurate wall shear stress and heat flux using instantaneous values from the LES. Shear stress and heat flux are provided to the LES in the form of flux across the wall in this model.
The Cartesian cut-cell method [9][10][11][12][13] subdivides the object surface into small flat planes inside the cubical computational cells. In other words, these cells are divided by small cut planes into the fluid side and the object side. The interaction between the object and the fluid is given by wall flux vectors across the cut plane. Due to this sharp interface treatment and the robust grid generation of the Cartesian mesh, the Cartesian cut-cell method is one of the most promising methods in CFD. However, the Cartesian cut-cell method and other Cartesian mesh solvers [14][15][16] struggle to satisfy the required grid resolution for boundary layers because rectangular cells with an extremely large aspect ratio used in boundary fitted mesh are not available with these methods. To overcome this difficulty, an appropriate near-wall treatment is required. The wall-stress model by Kawai and Larson is expected to be one of the most compatible wall treatment techniques with the Cartesian cut-cell method.
The present study performed numerical simulations of the flow past a triangular column and 30P30N high-lift airfoil configuration using the Cartesian cut-cell method with the wall-stress model. The numerical results of the present study revealed the prediction capability of the wall-stress model to the Cartesian cut-cell simulation of unsteady flows.

Governing Equation
The governing equation is the compressible Navier-Stokes equation. Conservation laws over a cuboid control volume (x 1 ≤ x ≤ x 2 , y 1 ≤ y ≤ y 2 , z 1 ≤ z ≤ z 2 ) in the fluid can be written in the following integral form: . The conserved quantities Q; the advection fluxes E a , F a , and G a ; and the viscous and conductive fluxes E d , F d , and G d in Equation (1) are defined as follows: where ρ is the density; u x , u y , and u z are the components of the velocity; p is the pressure; and e is the total energy per unit volume. The variables τ and q are the stress tensor and heat flux, respectively: (3) where a is the speed of sound; µ is the viscosity obtained by Sutherland's law; and the Prandtl number Pr is 0.72. The conservation equation in integral form over a cuboid control volume that extends over the two regions of the air and the object is written in the following form using the Heaviside function H(x, y, z): where the Heaviside function H(x, y, z) gives H = 1 in the air and H = 0 in the object; and the vector n represents the unit normal vector directed outward form the object into the air. Considering the boundary condition u·n = 0 at the object surface, Equation (5) is reduced to the following form: where the variable σ on the right-hand side of Equation (6) represents the wall flux vector (the interaction between the object and the air), which is directed from the object into the air: The vector b is the stream-wise unit vector parallel to the wall. The scalar value τ w of the wall shear stress is obtained by the wall-stress model. Wall flux of mentum through the pressure and the shear stress are included, while wall flux of mass and energy do not exist in the problems considered here. The governing Equation (6) is solved using a cell-centered finite volume method. Cell-averaged quantities of the finite volume method are defined as follows: Hdxdydz, Hdydz, Hdxdy, where α is the volume fraction of the air in the cell; β x , β y , β z are the area fractions of the air in each cell interface; V is the volume of the cell; and S x , S y , and S z are the area of each cell interfaces. Using these cell-averaged values, governing Equation (6) is reduced to the following form: The projection of the cut plane to each cell interface is approximated in the following form: These approximation tend to be exact equations at the limit of small curvature of the cut plane. Thus in the case of the sufficiently small control volume, S cut-plane and n are evaluated as follows: The advection flux at the cell interface is calculated using the Simple Low-dissipation AUSM (SLAU) scheme [17]. The primitive variables on the cell interface are reconstructed using Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) of the fifth-order accuracy [18]. The velocity components are obtained using the modification of the MUSCL scheme proposed by Thornber [19]. The time integration to second-order accuracy is calculated using the Total Variation Diminishing (TVD) Runge-Kutta method [20].

Wall-Stress Model
The wall-stress model by Kawai and Larson [5,7,8] for the turbulent boundary layer was used in the present study. Stream-wise momentum and the total energy conservation equation of the wall-stress model are written in the following ordinary differential equations in the local wall-normal coordinate η: where U is the wall-parallel velocity component. Constants Pr = 0.72 and Pr t,wm = 0.90 are the Prandtl number and turbulent Prandtl number, respectively. The solution of the wall-stress model equations shows the averaged velocity and the averaged temperature distributions in the inner-layer part of the turbulent boundary layer. The turbulent eddy viscosity in the wall-stress model µ t,wm is obtained using the following formula: The von Kármán constant is κ = 0.41, and A + is a constant that equals 17. The velocity scale u τ = τ w /ρ is estimated by the instantaneous wall shear stress and density. The dimensionless wall distance y + in Equation (16) is defined as: The boundary value problem of the wall-stress model Equations (13) and (14) are discretized using the one-dimensional finite volume method and solved by the Thomas algorithm. Boundary values at η = η max = d wm are determined from the instantaneous value of the CFD analysis (Figure 1a). A schematic of the application of the wall-stress model to the Cartesian cut-cell method is shown in Figure 1. An overview of the implementation is provided as follows: The upper bound of the wall-stress model mesh is set on the wall-normal line passing through a CFD cell center (x c , y c , z c ) that includes the cut-plane by the object. The upper bound is placed at distance d wm from the wall. The distance d wm corresponds to y + ∼ 300. The location of the upper bound (x ub , y ub , z ub ) of the wall-stress model in the CFD domain is decided as follows: where ϕ c denotes the distance from the wall to the CFD cell center (x c , y c , z c ). b.
The one-dimensional non-uniform mesh for the wall-stress model is generated from the wall to the upper bound. The wall-stress model mesh is generated using the following formula: where j wm denotes the cell number from the wall; r wm = 1.1 denotes the rate of the increment of the grid distancing between neighboring grid points. The grid distancing increases by 10% between neighboring grid points. Constants c 1 and c 2 is decided to satisfy the following conditions: c. The upper boundary values U and T are decided by inverse distance weighted interpolation using instantaneous values of the neighbor cell-center of the Cartesian cut-cell simulation as follows: The wall boundary condition of the wall-stress model equations is the non-slip and adiabatic wall in the present study. d.
Wall-stress model Equations (13) and (14) are solved using the finite volume method using the one-dimensional wall-model mesh. e.
The obtained wall-stress τ w is provided to the CFD cell at (x c , y c , z c ), which includes the cut-plane by the object.

Flow around a Triangular Column
The flow around a triangular column that contains only sharp-edge separation and hardly depends on the Reynolds number is calculated. The results of the wall-modeled Navier-Stokes simulation is compared using a experiment and an inviscid simulation by the Euler equation.

Simulation Conditions
The flow condition is shown in Table 1.  [21]. Boundary values at y min and y max are fixed to the free-stream condition. The grid system in the z-direction consists of only a uniform grid, and the periodic boundary condition is applied.

Simulation Results
Isosurfaces of the second invariant of the velocity gradient tensor (Q-criterion) of the fine grid simulations are shown in Figure 2. The left figure shows the result of the Euler simulation and the right figure shows the wall-modeled Navier-Stokes simulation. In both figures, flow separation occurs at the downstream-side sharp edges, and a complex vortex structure that contains longitudinal vortices is formed in the wake. The outline of the flow structure by the wall-modelled Navier-Stokes simulation is the same as that of the Euler simulation.
Instantaneous distribution of pressure and velocity in an xy-plane of the fine grid simulation is shown in Figure 3.  Free shear layers form in the wake of the sharp edges, and the vortices form by the rollup of the free shear layers. Pressure in the wake of the wall-modelled Navier-Stokes simulation is similar to that of the Euler simulation.
Time-averaged distribution of pressure and velocity in an xy-plane of the fine grid simulation is shown in Figure 4.  The length of the wake of the wall-modelled Navier-Stokes simulation is almost equal to that of the Euler simulation. Pressure distribution in the wake of the wall-modelled Navier-Stokes simulation is slightly different from that of the Euler simulation. However, the surface pressure distributions of the triangular column in both cases do not have noticeable differences. This suggests that the flow field hardly depends on the Reynolds number, and the boundary layer on the triangular column does not play an important role.
The drag coefficient is shown in Table 2. The drag coefficient of the fine grid of the Euler simulation is 1.282 and the difference from the experimental value [22] is 3.2%. The difference decreases with grid resolution. The pressure drag coefficients of the fine grid simulation of the wall-modelled Navier-Stokes simulation is 1.284. This value is extremely close to that of the Euler simulation. The difference in the drag coefficient of the wall-modelled Navier-Stokes simulation from the experiment is slightly less than that of the Euler simulation. This improvement is mainly provided by the additional surface friction of the wall-modelled Navier-Stokes simulation. Considering effect of viscosity of the fluid in the wall-modelled Navier-Stokes simulation, this relationship is reasonable.
The computational cost of the wall-modelled Navier-Stokes simulation is 3/2 or more of that of the Euler simulation. Both the wall-modelled Navier-Stokes simulation and the Euler simulation are valid when the separation lines are fixed to the sharp edges of the object. Consequently, it is possible to select one of these simulation methods according to accuracy and computational cost requirements when the separation lines are fixed to the sharp edges and the flow field hardly depends on the Reynolds number.

Flow around the 30P30N Three-Element High-Lift Airfoil Configuration
Numerical Euler and wall-modelled Navier-Stokes simulations of the flow around a 30P30N high-lift airfoil configuration were performed. This airfoil is a model of the high-lift device of an aircraft that consists of a slat, main element, and flap. The simulation results were compared with experiments and CFD results using the boundary-fitted mesh method. Based on these comparisons, the prediction capability of the Cartesian cut-cell method with the wall-stress model by Kawai and Larson [5] was discussed.

Simulation Conditions
The flow conditions are provided in Table 3. The chord length c of the 30P30N high-lift airfoil configuration is 0.457 m following the experiment by Murayama et al. [23]. The grid system consists of a uniform grid domain and a non-uniform grid domain in the xy direction. The range of the uniform grid domain that consists of cubic cells with constant grid width ∆x is −0.1 m ≤ x ≤ 0.6 m, −0.1 m ≤ y ≤ 0.1 m, and 0 ≤ z ≤ 0.052 m. Outside the uniform grid domain, the width of the cuboid cells stretches by 10% between neighbors in the x and y directions. The boundary condition at the outer edges of the xy-plane is determined by Riemann invariants [21]. The computational grid in xy-plane is shown in Figure 5. The grid system in the z direction only consists of a uniform grid, and the periodic boundary condition is applied.

Simulation Result by Wall-Modelled Navier-Stokes Simulation
The isosurface of the second invariant of the velocity gradient tensor (Q-criterion) around the 30P30N high-lift airfoil configuration at the angle of attack θ = 5.5 • is shown in Figure 6. Vortices are shed from the trailing edge of the slat (Figure 6a), the main element, and the flap. The intermittent vortices of the outer layer of the turbulent boundary layer are generated on the suction side of the main element ( Figure 6b) and flap. Figure 7 shows the distribution of the velocity component u x at the nearest cell from the airfoil surface and an xy-plane cross-section of the fine grid simulation, providing a close-up view at x/c ∼ 0.15. The dark blue area in Figure 7 shows the negative value of the velocity component, which indicates that the leading-edge separation bubble forms at x/c ∼ 0.15. The location of the separation bubble is consistent with the experiment [24]. This separation bubble formed in the fine grid simulation but not in the medium and coarse grid simulations. Flow separation occurs on the suction side of the upstream side of the trailing edge of the flap in the instantaneous flow field as shown in Figure 8a. It does not reattach before the trailing edge of the flap. In contrast with Figure 8b, the flow does not separate in the time-averaged flow field (Figure 9). This indicates that the separation near the trailing edge of the flap occurs intermittently.
The lift coefficients in this study and the θ-sweep of that in the experiment in Reference [23] are shown in Figure 10.  The Cartesian cut-cell simulation with the wall-stress model accurately predicts the lift coefficient. The lift coefficient of the Cartesian cut-cell simulation with the wall-stress model, the boundary-fitted mesh simulation [25], the hybrid mesh simulation [26], and the experiment [23] are shown in Table 4. The hybrid mesh in Reference [26] is composed of the boundary fit mesh in the near-field region and the quadtree-based Cartesian mesh in the far-field region.
Compared with the experiment, the Cartesian cut-cell simulation with the wall-stress model underestimates the lift coefficient. The difference from the experiment decreases monotonically with increasing grid resolution. However, the difference decreases slower than expected accuracy in the formulation.  [25] 70,445,430 2.77 (1.8%) Hybrid mesh (SA DDES) [26] 242,000,000 2.95 (4.6%) Experiment [23] 2.82 The pressure coefficient in the case of θ = 5.5 • is shown in Figure 11. The surface pressure distribution at the pressure side of this study coincides with the experiment and referenced simulations. It indicates that the surface pressure distribution at the pressure side is independent of the grid type and the near-wall treatment.
The Cartesian cut-cell simulation of the present study well estimates the surface pressure distribution as well as the boundary-fitted mesh simulation [25] and hybrid mesh simulation [26], except for the suction side of the slat and the main element suction side of x/c ≤ 0.15. The curvature of the pressure distribution changes around x/c = 0.15, and this x-coordinate coincides with the location from which the vortices of the outer layer of the turbulent boundary layer develop on the suction side of the main element.
The underestimation of the lift coefficient of the present study in Table 4 is consistent with the surface pressure distribution in Figure 11. The positive error of the surface pressure on the suction side of the main element at x/c ≤ 0.15 and the suction side of the slat results in the underestimation of the lift. This error is caused by the mismatched estimation of the wall shear stress of the laminar parts of the boundary layers in these areas. The Cartesian cut-cell method experiences difficulty with resolving the laminar boundary layer due to its incompatibility with the grid concentration in the wall-normal direction. This problem is expected to be solved by the future development of the wall-stress model for laminar boundary layers.
The results of the present simulations showed that Cartesian cut-cell simulation with the wall-stress model by Kawai and Larson [5] is one of the most useful schemes for the external flows of Re ∼ 10 6 . The difference in the lift coefficient in the experiment in the medium grid simulation in Table 4 is the same as the hybrid mesh simulation [26]. Considering the total number of computational cells, the prediction capability of the present Cartesian cut-cell simulation with the wall-stress model is the same as the hybrid simulation [26]. Future improvements in the wall-stress model for the laminar part of boundary layer is expected to produce more accurate prediction within the same degree of accuracy as the boundary-fitted mesh simulation by Sakai [25]. The lift coefficient of the Euler simulation in the medium grid at θ = 5.5 • is 3.25 and the difference from the experimental result is 13.2%. The Euler simulation overestimates the lift coefficient, and the difference in the lift coefficient in the edium grid of the Euler simulation is obviously larger than that of the coarse grid of the wall-modelled Navier-Stokes simulation. The Euler simulation is valid only if separation lines are fixed at the sharp edges of the body as shown in Figure 2a.

Conclusions
The prediction capability of the Cartesian cut-cell method with a wall-stress model was discussed in this paper.
In the numerical analysis of the flow around a triangular column, the wall-modelled Navier-Stokes simulation adequately predicted the drag coefficient. In this case of sharp-edge separation flow, the drag coefficients of both the wall-modelled Navier-Stokes simulation and the Euler simulation agreed with the experimental value within a small error. The value of the drag coefficient of the wall-modelled Navier-Stokes simulation was closer to the experimental value compared to that of the Euler simulation.
The difference of the lift coefficient between the Cartesian cut-cell simulation and the experiment monotonically decreases with increasing grid resolution, in the numerical analysis of the flow around the 30P30N three-element high-lift airfoil configuration. The intermittent vortex structure of the outer layer of the turbulent boundary layer was observed on the suction side of the main element and the flap. The intermittent flow separation from the suction side of the airfoil surface of the flap was predicted.
The prediction result of the flows at Re ∼ 10 6 by the Cartesian cut-cell method with the wall-stress model were equivalent to that by the boundary-fitted mesh methods, except for the laminar parts of the boundary layers. The Cartesian cut-cell method with the wall-stress model is one of the useful methods for high-Reynolds-number flows at Re ∼ 10 6 .