A Numerical Scheme Based on an Immersed Boundary Method for Compressible Turbulent Flows with Shocks: Application to Two-Dimensional Flows around Cylinders

A computational code adopting immersed boundarymethods for compressible gas-particlemultiphase turbulent flows is developed and validated through two-dimensional numerical experiments. The turbulent flow region is modeled by a second-order pseudo skew-symmetric formwithminimum 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 skewsymmetric form significantly improves the resolution of the vortex generated in thewakewhile 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.


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][3][4][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 particleturbulence interactions, further fundamental analyses are necessary.Therefore, we consider a multiscale analysis of gasparticle 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.Target of this study: 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][13][14][15], recently.Takahashi et al. have been developing several Cartesian flow solvers [15][16][17][18][19][20][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 twodimensional 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].

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,  V and  V 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. +1/2 is defined in terms of the Ducrostype sensor [26] as follows: Here  = 10 −15 is a small value that prevents division by zero and  = 0.6 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,  +1/2 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:

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  IP in Figure 3, is an important parameter that eliminates recursive interpolation.Here we fix the length of  IP as 1.75 times the mesh size, considering the extension to the three-dimensional problem.Because the probe is longer than √3 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 √2/2 times the mesh size.This value √2/2 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.

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.

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. The blue, white, and red regions show fluid, ghost, and object cells, respectively.The boundary layer thickness is roughly estimated as 1/√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.

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.
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.
While the switching scheme improves at finer mesh resolution in both cases, scheme B shows good performance at all mesh resolutions.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 /√ Re as 0.058.

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.
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.

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.  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.

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 M 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.
Journal of Applied Mathematics

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.Figure 21 displays the distribution of the switching parameter  at  = 23.2.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  1 and  2 , 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  = 11.1 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.

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  = 7.8, 12.9, 18.1, and 22.0.The shock waves propagated from the cylinders diagonally interfere ahead of the cylinders at  = 7.8.The shear flow in the wake is disordered following the crossing at  = 18.1.The flows formed by the interaction of wake and shocks are resolved well by the present switching scheme (see Figure 24).
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  1 is enlarged by the shock wave and shear flow from  2 , similar to application 1.However, the force coefficient falls to zero, negated by the shear flow of  2 .Along the axis, the mesh size introduces a 10% variation in the peak coefficient of  2 (occurring around  = 10).At the peak, the flows are highly compressed by the shock waves from  1 and  2 .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  2 .Along the -axis, the force coefficient of  2 enhances around = 10, as the circulation around  2 is diminished by  1 intercepting the shear flow.

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,

2 Journal
analysis for the small region) Modeled, cheap (leading to analysis for the wider region)

Figure 1 :
Figure 1: Overview of the proposed gas-particle simulation.

Figure 2 :
Figure 2: Cell construction and classification in the present level set method.

Figure 3 :
Figure 3: Schematic of the present ghost cell approach with image point.

Figure 4 :
Figure 4: Image point projected from a fluid cell to compute the fluid force on the surface.

Figure 5 :
Figure 5: Flowchart of the present numerical method.
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 fluidcell neighboring boundary are updated.If the flows contain

Figure 6 :
Figure 6: Schematic of the computational domain.

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.

Figure 14 :
Figure 14: Effect of mesh size on aerodynamic coefficients at Mach 1.2.The red lines are reference values.

Figure 17 :Figure 18 :
Figure 17: Effect of mesh size on drag coefficients at Mach 2.0.The red lines are reference values.

Figure 19 :
Figure 19: Computational configurations of the applications discussed in the text.

Table 1 :
Numerical schemes for estimating inviscid flux.