Abstract
A computational code adopting immersed boundary methods for compressible gas-particle multiphase turbulent flows is developed and validated through two-dimensional numerical experiments. The turbulent flow region is modeled by a second-order pseudo skew-symmetric form with minimum dissipation, while the monotone upstream-centered scheme for conservation laws (MUSCL) scheme is employed in the shock region. The present scheme is applied to the flow around a two-dimensional cylinder under various freestream Mach numbers. Compared with the original MUSCL scheme, the minimum dissipation enabled by the pseudo skew-symmetric form significantly improves the resolution of the vortex generated in the wake while retaining the shock capturing ability. In addition, the resulting aerodynamic force is significantly improved. Also, the present scheme is successfully applied to moving two-cylinder problems.
1. Introduction
The acoustic waves from rocket plumes are sufficiently strong to damage the satellites inside the fairing of a rocket. These waves are widely assessed by empirical prediction methods [1], but the low accuracy of these methods renders them unsuitable for new rocket launch sites. To improve the prediction accuracy of acoustic waves from rocket plumes, a sophisticated model based on the underlying physics is required. Numerical simulations are an essential component of new model development [2–5]. The behavior of acoustic waves from rocket plumes is difficult to predict, because actual plumes are very hot, with very high speed, and of the multiple phase conditions. In real rocket systems, acoustic waves are suppressed by a water injection system installed at the rocket launch site. Fukuda et al. [6] showed that rarefaction or absorption of acoustic waves by particles exerts no significant effect and that acoustic waves might be primarily attenuated by interactions between particles and turbulence. However, this scenario is not well modeled in their study. To more accurately evaluate acoustic wave suppression by particle-turbulence interactions, further fundamental analyses are necessary. Therefore, we consider a multiscale analysis of gas-particle multiphase high-speed compressible flows. Because the target is a rocket plume, we propose a simultaneous treatment of the turbulence and the shock waves. Figure 1 shows an overview of the proposed numerical approach. The simultaneous high-resolution simulation of the particles and turbulence is conducted on the microscale, the large eddy simulation (LES) modeling only the particle behavior is conducted on the intermediate scale, and the complex flow fields are modeled by the Reynolds-averaged Navier-Stokes (RANS) simulation.
Speeds of these scattering particles cover a wide range of Mach numbers from subsonic to supersonic, while Reynolds numbers are quite low since the sizes of the particles are small. Therefore, flowfields around the small particles can be solved by flow simulation without any turbulence models though the flowfield is macroscopically turbulent. There are several kinds of numerical methods to solve the problem treating a number of moving objects, such as dynamic mesh method [7] or overset method [8]. However, simple implementation and rapid computation are difficult to achieve in using these methods because additional numerical processes are included in the flow simulation like mesh regeneration or interpolation between computational grids. We select level set method [9] to track a number of moving particles in this study. There are some other options to trace many moving objects with high accuracy like phase field method or front tracking method. However, our main focus is investigation of the acoustic wave characteristics under interference between turbulence and particles. Therefore, level set method is selected based on computational efficiency. The represented boundaries by level set method are imposed by immersed boundary method [10] in equally spaced Cartesian mesh in this study. Immersed boundary method (IBM) is widely used in the community of Cartesian mesh method from the simplicity and applicability. The methodology of IBM can be classified into several categories such as continuous forcing method [11], direct forcing method, and ghost-cell method [10]. Although IBM was originally proposed and used for incompressible flow simulations, it is also applied to compressible flow simulations [12–15], recently. Takahashi et al. have been developing several Cartesian flow solvers [15–21] and investigating the performance of this kind of numerical method. Based on the background, we adopted ghost-cell method [10, 15] with equally spaced Cartesian mesh by the points of simple implementation, robustness, and extensibility.
Here, we should recall our flow features that consist of both turbulence and shocks. In general, an upwind scheme is often employed to evaluate inviscid fluxes in IBM flow solver to stabilize a flow with numerical dissipation. The dissipation is not preferable to solve our flows of turbulence part clearly, while it should be added appropriately to capture the shocks in a part of our flows. In other words, we should minimize numerical dissipation in the turbulent region, while the dissipation must be added to prevent spurious oscillations around shocks. In resolving both the turbulence and shock waves, dissipation switching scheme can play a major role. In this study, a switching procedure for numerical dissipation is introduced and examined through the two-dimensional test cases. This study overviews the computational code and demonstrates its efficacy in simulations of two-dimensional static cylinder flows under various Mach number conditions. The present numerical method is developed to three-dimensional problem of interference among turbulence, shocks, and particles with high-performance parallel computing based on the previous studies [16, 17, 20, 21].
2. Computational Methods
2.1. Governing Equations and Numerical Method
In the present study, flows are governed by two-dimensional compressible Navier-Stokes equations. No averaging or filtering process is involved and the flows are solved without any turbulence models: where and are the inviscid fluxes in the - and -directions, respectively, and are the corresponding viscous fluxes, and contains the conservative variables. Here the stress tensor components are given as The pressure is related to the total energy per unit mass by the equation of state: The heat flux term in the equation of total energy is computed by where the equation is transformed by the ideal gas equation and Prandtl number as follows: All variables are nondimensionalized by the freestream conditions of density, sound speed, and unit length. The above equations are discretized on an equally spaced Cartesian mesh with a cell-centered arrangement. To eliminate additional numerical dissipation everywhere, except in the vicinities of shock waves and potential flows, the inviscid terms are computed by a hybrid scheme that combines the pseudo skew-symmetric central difference scheme [22] with the monotone upstream-centered scheme for conservation laws (MUSCL)-Roe scheme [23, 24]. The flux in the turbulent region is modeled by a pseudo skew-symmetric central difference scheme with a minimum dissipation term as follows: Here the subscript denotes the quantity on the th grid point and subscripts and denote the left- and right-side states, respectively, interpolated by the third-order MUSCL scheme [24] with van Albada’s limiter [25].
On the other hand, the flux in the shock and potential flow region, computed by the second-order MUSCL-Roe scheme, is written as follows: where is the flux Jacobian and the subscript Roe denotes the Roe-averaged quantity of the left and right states. Here where and are the right and left eigenmatrices of , respectively, and is a diagonal matrix.
The symmetric central difference part of (8) can be replaced by that of the pseudo skew-symmetric form without losing the formal order of accuracy of the equation. The proposed scheme adopts the following new form of (7): Combining this with the digital switching function, we obtain the following hybrid scheme: Excessive dissipation is added to the shock or potential flow region when beta is unity, whereas dissipation is minimized when beta is zero. is defined in terms of the Ducros-type sensor [26] as follows: Here is a small value that prevents division by zero and is the switching threshold. The divergence and rotation in (11) are computed by a second-order finite difference scheme. In the present study, the Ducros-type sensor [26] alone is used in both the shock and potential flow regions, although previous studies have combined this sensor with the Jameson sensor [27] in the shock region [26, 28]. Furthermore, in one of our proposed schemes, is set to unity in cells close to the body. Finally, the flux derivative is approximated as follows: The derivative is obtained similarly.
The diffusive terms are treated by a second-order, central difference scheme using the mid-point flux. The time marching is conducted by the three-stage total variation diminishing Runge-Kutta scheme [29]. In this study, time increment is determined as follows:
2.2. Boundary Representation
The boundary is defined by the level set method [9, 14]. The level set function is determined in whole cells as a signed distance from the object boundary. A schematic of the cells around a boundary is shown in Figure 2. The level set method effectively computes the normal vector from the object surface on the basis of a gradient operation. In the present study, flows around multiple moving objects are solved by extending the level set method to multiple level set functions based on simple minimum distance approach [8].
On the basis of the level set function (14), all cells are classified into three categories: fluid cell, ghost cell, and object cell, as shown in Figure 2. The ghost cells behave as guard cells between the fluid and object regions and are assigned in two layers under the present definition as follows:
Ghost cells are used for imposing boundary condition in the present method [10, 15]. An image point set in the region of fluid cells is used to collect flow information for a ghost cell. A primary advantage of the present ghost-cell method adopting the image point approach is its robustness. A schematic of the present immersed boundary method is shown in Figure 3. The image point is located at the edge of a probe that extends from a ghost cell through the object boundary in the direction normal to the surface. The length of the probe, denoted as in Figure 3, is an important parameter that eliminates recursive interpolation. Here we fix the length of as 1.75 times the mesh size, considering the extension to the three-dimensional problem. Because the probe is longer than times the mesh size, the nodes enclosed by the orange dotted square in Figure 3 are classified only as fluid cells. The primitive variables on the image point are interpolated by the bilinear function based on the surrounding cells. Finally, the value of the ghost cell is defined by the value of the image point. The flow velocity, computed by (15), is assumed to be linearly distributed along the probe from the boundary point to image point. For the pressure and density, a Neumann condition is assumed on the boundary and the ghost cell is assigned the value of its corresponding image point:
2.3. Estimation of Friction Force
The fluid force acting on an object is estimated by surface integration. Fluid forces must be evaluated at both ghost and fluid cells because the object boundary can straddle both cells, as shown in Figure 4. For the ghost cells, the image point for the immersed boundary method is directly available for calculating the fluid force on the surface. In the case of fluid cells, on the other hand, the small distance between the surface and fluid cell can magnify the friction force wrongly. Therefore, the friction force on the fluid cell is estimated by using the image point method with changed probe length to times the mesh size. This value is determined to adopt only fluid cells as interpolation cells for the image point. In this case, the viscous force is estimated by (16). The velocity component in the following equation corresponds to the component parallel to the boundary:
Our code supports a simpler but accurate force estimation method without the information of a level set function [30], but this option is not considered here.
2.4. Outer Boundary Conditions
While supersonic flows can be solved by applying Neumann conditions at the outer boundary, this approach can induce instabilities in subsonic flow simulations. Therefore, we apply a sponge region [31] near the outer boundary to suppress sound reflections with imposing a Dirichlet condition on the density at the outflow boundary. This boundary condition ensures numerical stability of flows with vortical wakes.
2.5. Coding Flowchart
Figure 5 shows a flowchart of the computational procedures. The loop of the three-stage Runge-Kutta time integration is iterated three times in the blue box. Immediately before flux computation, the ghost cells that specify the inner boundary condition of the fluid-cell neighboring boundary are updated. If the flows contain moving objects, the level set function, cell classification, and image point identification are reexecuted before the Runge-Kutta loop as shown in the red box. Prior to computing the viscous forces using the ghost cell values, the ghost cells are again updated.
3. Computational Result of Flow around a Two-Dimensional Circular Cylinder
3.1. Computational Condition
The proposed scheme is evaluated through a series of numerical tests. The proposed schemes are compared with the conventional MUSCL-Roe scheme at different grid resolutions. Three methods for estimating the inviscid flux as shown in Table 1 are compared: (A) the present scheme (10) over the whole region, (B) enforcing the MUSCL-Roe scheme for nearby objects and the present scheme for other regions, and (C) the MUSCL-Roe scheme over the whole region. In the present implementation, there is no difference about the computational cost to calculate the fluxes based on (A), (B), and (C). The characteristics of the schemes are investigated on subsonic and supersonic flows around a two-dimensional circular cylinder. The Reynolds number, based on the cylinder diameter and freestream values (including viscosity), is fixed at 300, while the freestream Mach numbers are varied as Mach 0.3, 1.2, and 2.0. Four mesh sizes are compared: 0.200, 0.100, 0.050, and 0.025, where the diameter of the circular cylinder is fixed at 1. The parameters in all trials are summarized in Table 2. The computational domain is shown in Figure 6. We investigated the validity to employ the size of computational domain with comparing square domain of 40. The flow simulations are conducted with a Courant number of 0.4. It was confirmed that there was no significant effect about the Courant number with comparing results from Courant number 0.2. Dirichlet conditions are imposed on all flow variables at the inflow boundary and on density alone at the outflow boundary. Other variables are assigned Neumann conditions at the outflow boundary. Neumann conditions are imposed at the top and bottom boundaries for all variables. The circular cylinders represented by the four mesh sizes are illustrated in Figure 7. The black solid lines are grid lines connecting cell centers. The bold red line is the boundary of the circular cylinder immersed in the Cartesian mesh. The blue, white, and red regions show fluid, ghost, and object cells, respectively. The boundary layer thickness is roughly estimated as Re, that is, 0.058. Consequently, the boundary layer is discretized by two or three cells in using 0.025 mesh size at Reynolds number 300.
In addition, the present results are compared with the recent “state-of-the-art” WENO body-fitted coordinate computational code. This code is based on the sixth-order WENO-CU developed by Hu et al. [32] and extended to the body-fitted coordinate system by Nonomura et al. [33]. The number of grid points is 208 × 177. The result obtained by this code is used as a reference solution. We confirmed that the 208 × 177 mesh is sufficiently fine to generate a reference solution by comparing the results with those obtained on a finer grid (410 × 268). The results output by this code are presented in the Appendix.
3.2. Comparison of Computational Results at M 0.3
Figure 8 shows the density distribution at Mach 0.3. The top, central, and bottom columns show the results from schemes A, B, and C, respectively. The left and right panels were computed on mesh sizes of 0.200 and 0.050, respectively. Clearly, the vortices are collapsed when the MUSCL-Roe scheme is implemented on the coarse mesh (M03-C-0200D; bottom left of Figure 8). On the other hand, the present scheme A preserves the flow features at all mesh resolutions. All schemes yield similar results on the fine mesh.
(a) M03-A-0200D
(b) M03-A-0050D
(c) M03-B-0200D
(d) M03-B-0050D
(e) M03-C-0200D
(f) M03-C-0050D
The differences among the three schemes are confirmed from pressure distributions and the time history of aerodynamic coefficients at mesh size 0.050. The results are shown in Figure 9. In scheme A, although high vortex resolution is achieved in the pressure contours, pressure oscillation is also observed near the cylinder. The oscillation is not preferable for stable computation. On the other hand, scheme C yields smooth pressure distribution and history of aerodynamic coefficients. The vortices in scheme C are more dissipative than those of scheme A. The strong features of schemes A and C, vortex conservation and pressure oscillation suppression, are realized in scheme B.
(a) M03-A-0050D
(b) M03-A-0050D
(c) M03-B-0050D
(d) M03-B-0050D
(e) M03-C-0050D
(f) M03-C-0050D
Figure 10 plots the distributions of the switching parameter in (11). Schemes A and B are adopted at mesh sizes 0.100 and 0.025. The black and white regions are solved by the present scheme (10) and MUSCL-Roe [22, 23] scheme, respectively. On the coarse mesh, the switching function is rendered ineffective by numerical error and small perturbation. In scheme B, however, where the MUSCL-Roe scheme is applied only near the object, the performance of the switching function is superior to that of scheme A. While the switching scheme improves at finer mesh resolution in both cases, scheme B shows good performance at all mesh resolutions.
(a) M03-A-0100D
(b) M03-A-0025D
(c) M03-B-0100D
(d) M03-B-0025D
Figure 11 summarizes the results at Mach 0.3. The reference values were obtained from a boundary-fitted mesh simulation shown in the Appendix. Although grid convergence is almost obtained for the lift coefficient amplitudes, Strouhal number, and average drag coefficient, the convergence is not monotonic. One of the reasons of the inflected convergence that can occur is interference between the present switching scheme of different spatial accuracy (third-order MUSCL-Roe and second-order pseudo skew-symmetric) and immersed boundary method. The lift and drag coefficients are overestimated and underestimated, respectively, in schemes A and C. In the case of scheme B, the lift and drag coefficients are intermediate between schemes A and C. The friction drag coefficient follows the same trend in all schemes and never completely converges. The boundary layer, which is discretized by several grid points even in the finest grid resolution 0.025, is not sufficiently resolved to show grid convergence since the thickness is roughly estimated by as 0.058.
(a) Average drag coefficient
(b) Average pressure drag coefficient
(c) Average friction drag coefficient
(d) Amplitude of lift coefficient
(e) Strouhal number
3.3. Comparison of Computational Results at M 1.2
Figure 12 plots the density contours over the whole computational region and in the near field of the object at mesh size 0.100. Obtained flowfield becomes almost symmetrically different from the previous subsonic case. The trend of the flow resolution is similar to the previously discussed subsonic case. Scheme A captures a sharper distribution than scheme C but develops weak numerical oscillation. Reasonable results are obtained by scheme B, in which regions far and near the object are evolved under schemes A and C, respectively.
(a) M12-A-0100D
(b) M12-A-0100D
(c) M12-B-0100D
(d) M12-B-0100D
(e) M12-C-0100D
(f) M12-C-0100D
The distribution of at mesh size 0.100 is similar among the three schemes (Figure 13). At Mach 1.2, the white region (solved by scheme C) enlarges relative to that at Mach 0.3. The dependency of drag coefficient on grid resolution is similar to that of Mach 0.3 (Figure 14), although the drag coefficients are more similar among the three schemes. A likely reason for this trend is that the region solved by scheme C becomes wider in Mach 1.2 than in Mach 0.3.
(a) M12-A-0100D (far)
(b) M12-A-0100D (close)
(c) M12-B-0100D (far)
(d) M12-B-0100D (close)
(a) Drag coefficient
(b) Pressure drag coefficient
(c) Friction drag coefficient
3.4. Comparison of Computational Results at M 2.0
At Mach 2.0, the flows computed by scheme A destabilize even under various restart conditions probably because the present scheme consists of little numerical dissipation. Figure 15 plots the density distribution calculated by schemes B and C at mesh size 0.100. While the distributions yielded by both schemes are very similar, those of scheme B are slightly sharper than those of scheme C.
(a) M20-B-0100D
(b) M20-B-0100D
(c) M20-C-0100D
(d) M20-C-0100D
The distributions of the switching parameter at mesh sizes 0.100 and 0.025 are displayed in Figure 16. Both distributions are obtained from scheme B because computations by scheme A blow up. The black regions solved by the central difference scheme appear both upstream and downstream of the cylinder. Thus, the present switching scheme adequately captures nondissipative flows. However, the drag coefficients calculated by schemes B and C are very similar, a likely consequence of the wide MUSCL-Roe region (see Figure 17). Thus, while the wake resolution is unambiguously clarified, the drag coefficients are not affected.
(a) M20-B-0100D (far)
(b) M20-B-0100D (close)
(c) M20-B-0025D (far)
(d) M20-B-0025D (close)
(a) Drag coefficient
(b) Pressure drag coefficient
(c) Friction drag coefficient
3.5. Comparison of Surface Pressure Coefficient at M 1.2 and M 2.0
Here, surface pressure coefficients in supersonic cases are compared with BFC results to investigate resolution near the boundary. Figure 18 shows the pressure coefficient distributions obtained from all schemes with fine (0.025) and coarse (0.100) mesh resolutions. The present results from fine and coarse meshes are drawn by circles and other symbols, respectively. Discrepancies are observed between BFC and results from coarse meshes, while almost agreements are obtained in the cases of fine mesh resolution. In the case of 1.2, however, the stagnation pressure is overestimated rather than BFC even in all cases of fine mesh. It can be affected from the size of upstream region from the cylinder and the location of the shock wave. Although scheme A shows oscillatory distributions due to the characteristic of the central difference scheme, schemes B and C show almost same distributions without the feature. In the case of coarse meshes, the beginning of the adverse pressure gradient is not captured clearly; that is, the resolution of the present IBM should be investigated precisely in addition to the effect of the mesh resolution.
| (a) 1.2 |
| (b) 2.0 |
4. Flow Simulation around Relatively Moving Cylinders
The present study aims to solve flows around multiple moving objects. To this end, we simulate the flows induced by two moving cylinders at Mach 1.2. Schematics of the applications are shown in Figure 19. All cylinders in this section are forced to move with fixed velocities and directions. All computations are performed by scheme B, and the diameters of all cylinders are 1. The Reynolds number, based on the moving velocity, flow viscosity, and cylinder diameter, is fixed at 300. For comparison, the simulations are conducted on two mesh sizes: 0.050 and 0.025. In flow simulations with moving objects, the flow simulations are conducted with a Courant number of 0.4. Cartesian meshes can cause the so-called “fresh cell” problem because the cell properties alter with the moving boundary [34]. In this study, although we do not employ any special treatment for the fresh cell problem, the computations are accurate and stable.
(a) Application 1
(b) Application 2
4.1. Application 1: Passing Cylinders over Each Other
In the first application, flow is simulated around cylinders passing each other along the -axis. Figure 20 plots the distributions of density and vorticity magnitude at nondimensional times , 13.6, 19.1, and 23.2. The shock wave generated by a moving cylinder interacts with the wake of the partner cylinder. Vortices are generated from the interaction between the shock wave and shear layer of the wake.
(a) Density
(b) Vorticity
Figure 21 displays the distribution of the switching parameter at . The distribution of around a moving cylinder is similar to that around a fixed cylinder. The black region, solved by the central difference scheme, spreads in the wake and vortical regions. The nondissipative scheme encourages instabilities to develop in the shear layer.
The fluid force, estimated by surface integration on each object, is normalized identically to the usual aerodynamic coefficient based on the velocity of moving objects. Figure 22 plots the axial force coefficient calculated in this manner, as a function of nondimensional time. The red and blue lines represent the coefficients of cylinders and , respectively. The solid and dotted lines represent coarse and fine meshes, respectively. The coefficients are almost independent of grid resolution. Apart from the initial impulse, most of the variation is caused by the shock wave intercepting from the partner cylinder. Along the -axis, the force coefficient jumps at as the shock wave strikes and then decreases nonlinearly under interference between the shock wave from the partner cylinder and shear layer in the cylinder’s own wake. Finally, the coefficient returns to its initial value, having been reduced by half following the nonlinear variation. Along the -axis, on the other hand, the force coefficient peaks around 0.7 as the shock wave arrives. Thus, the numerical method allows quantitative evaluation of the fluid force around the moving cylinders.
| (a) -axis |
| (b) -axis |
4.2. Application 2: Crossing Cylinders
The second application is flow simulation around crossing cylinders. In this flowfield, the shock wave, wake, and objects mutually interact. Figure 23 shows the distributions of density and vorticity magnitude at , 12.9, 18.1, and 22.0. The shock waves propagated from the cylinders diagonally interfere ahead of the cylinders at . The shear flow in the wake is disordered following the crossing at . The flows formed by the interaction of wake and shocks are resolved well by the present switching scheme (see Figure 24).
(a) Density
(b) Vorticity
As in application 1, we now evaluate the axial force coefficients around a pair of crossing cylinders. The temporal changes are plotted in Figure 25. Initially, the axial force on is enlarged by the shock wave and shear flow from , similar to application 1. However, the force coefficient falls to zero, negated by the shear flow of . Along the -axis, the mesh size introduces a 10% variation in the peak coefficient of (occurring around ). At the peak, the flows are highly compressed by the shock waves from and . Moreover, the shock wave is damped by the shear flow around the cylinder. While the grid resolution affects the sharpness of both shock wave and shear flow, the decreased peak value at the higher grid resolution may be attributable to excitation of the shear flow around . Along the -axis, the force coefficient of enhances around , as the circulation around is diminished by intercepting the shear flow.
| (a) -axis |
| (b) -axis |
5. Conclusions
To enable simulation of high-speed gas-particle multiphase flows, we developed a high-resolution computational code that captures shock behavior and applied it to the compressible Navier-Stokes equations on an equally spaced Cartesian mesh. The second-order pseudo skew-symmetric and MUSCL-Roe schemes, together with the immersed boundary method, were combined into a hybrid scheme. The hybrid scheme yielded much higher vortex resolution than the MUSCL-Roe scheme while capturing the shock waves with the same effectiveness. The scheme was evaluated on Mach 0.3 subsonic flows and Mach 1.2 and 2.0 supersonic flows around a two-dimensional circular cylinder at Re = 300. The high resolution enabled accurate estimation of the aerodynamic forces. The modified hybrid scheme, which enforces the MUSCL-Roe scheme only around nearby objects, showed more accurate and stable characteristics than the original hybrid scheme because it dampens oscillations near the body. These oscillations are induced by the insufficient grid resolution near the object. The effectiveness of the scheme was emphasized in the flow simulations on coarse meshes.
Moreover, flows were simulated around two moving cylinders at Mach 1.2. The results of these simulations verified the applicability and robustness of the proposed method. Grid convergence results were obtained and the flow features, including interference among shock waves, shear flows, and objects, were well captured by the code developed in this study.
Our main focus is characteristics of acoustic wave in flow involving turbulence, shocks, and particles. Now we develop the present numerical method to three-dimensional simulation with investigating the applicability to turbulence and computational performance carefully to achieve high performance parallel computing. In near future, detailed phenomenon of this problem is going to be clarified from large-scale gas-particle flow simulation by using over billion computational cells.
Appendix
The flowfields obtained by the body-fitted coordinate code are shown in Figure 26.
(a) Computational mesh
(b) Density distribution at Mach 0.3
(c) Density distribution at Mach 1.2
(d) Density distribution at Mach 2.0
Nomenclature
| Speed of sound | |
| : | Courant number |
| : | Specific heat at constant pressure |
| : | Distance from object boundary (= value of level set function) |
| : | Cylinder diameter |
| : | Total energy per unit mass |
| : | Conservative variable vector |
| : | Inviscid flux vector along the -axis |
| : | Viscous flux vector along the -axis |
| : | Inviscid flux vector along the -axis |
| : | Viscous flux vector along the -axis |
| : | Pressure |
| : | Prandtl number |
| : | Conservative variable vector (used as a solution vector in this study) |
| : | Reynolds number |
| : | Time |
| : | Horizontal (-axis) velocity component |
| : | Vertical (-axis) velocity component |
| : | Velocity vector |
| : | Horizontal Cartesian coordinate |
| : | Vertical Cartesian coordinate |
| : | Density |
| : | Specific heat ratio |
| : | Switching function of inviscid scheme |
| : | Viscosity |
| : | Viscous stress tensor |
| : | Thermal conductance. |
| : | Freestream value |
| : | Image point |
| : | Fluid cell |
| : | Ghost cell |
| : | Object cell |
| : | Object value. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Authors’ Contribution
Shun Takahashi and Taku Nonomura equally contributed to this paper.
Acknowledgments
This research is partially supported by KAKENHI (24656522). The computations are carried out by NEC SX-9 at the Cyber Science Center of Tohoku University.