Abstract
This study is devoted to investigating a flow around a stationary or moving sphere by using direct numerical simulation with immersed boundary method (IBM) for the three-dimensional compressible Navier-Stokes equations. A hybrid scheme developed to solve both shocks and turbulent flows is employed to solve the flow around a sphere in the equally spaced Cartesian mesh. Drag coefficients of the spheres are compared with reliable values obtained from highly accurate boundary-fitted coordinate (BFC) flow solver to clarify the applicability of the present method. As a result, good agreement was obtained between the present results and those from the BFC flow solver. Moreover, the effectiveness of the hybrid scheme was demonstrated to capture the wake structure of a sphere. Both advantages and disadvantages of the simple IBM were investigated in detail.
1. Introduction
Acoustic wave from rocket nozzle on the rocket launch may critically damage satellites inside the payloads of the rocket since the wave sometimes becomes quite strong. The prediction of the crucial acoustic field around rocket launching sites has been conducted by the past enormous studies and testing data [1]. However, more accurate and reliable prediction is required to construct a new rocket launch site more safely. Various numerical studies have been conducted to realize qualitative prediction [2–4] until now. In real large-scale liquid launch vehicles, generated acoustic waves are suppressed by water injection. For high accurate prediction of the acoustic wave, the interaction between the water droplets and the flows should be clarified from the physical aspect.
Fukuda et al. showed that acoustic waves are primarily attenuated by the interactions between particles and turbulence [5] from theoretical analysis. These particles that interfered with turbulence in their study are alumina particles exhausted from the rocket nozzle or water droplets from water injection at the rocket launch. However, the mechanism of the phenomenon for the attenuation has not been well understood.
We develop the Cartesian flow solver with an immersed boundary method (IBM) [6] to investigate the interactions between particles and shocks by the direct numerical simulation (DNS). The objective of this study is to develop the flow solver to investigate the supersonic flow as rocket plumes around micro-meter-sized particles like alumina particles or water droplets. Accordingly, the flow becomes high-Mach- and low-Reynolds-number flow. To resolve real interaction between flows and particle in this study, an Euler-Euler approach is employed in the equally spaced Cartesian mesh with sharp interface IBM techniques. A hybrid scheme is also utilized to deal with flows of both shocks and vortical wake around a fixed or moving sphere. A drag coefficient is compared with numerical results obtained from the high accurate boundary-fitted coordinate (BFC) grid solver [7] and theoretical values [8]. We extended two-dimensional numerical methods [6] to three-dimensional one to solve flows around multiple particles passing through shock waves. The aim of this study is to investigate gas-particle flow involving real interaction between flows and particles by large-scale flow simulation. One of the most significant points for this purpose is to keep the algorithm simple for optimization and fast computation. Both simplicity and applicability of the code should be preserved in this study.
2. Numerical Method
2.1. Governing Equations
The governing equations of the present flow simulations are three-dimensional compressible Navier-Stokes equations:where , , and are inviscid fluxes in the , , and directions, , , and are viscous fluxes, and contains conservative variables. The pressure is related to the total energy per unit mass by an equation of state:Here, viscous stress tensor components are given as
All variables are normalized by the freestream values of the density, the sound speed, and the reference length. The above equations are discretized on an equally spaced Cartesian mesh with 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 third-order monotone upstream-centered scheme for conservation laws- (MUSCL-) Roe scheme [9] and the second-order pseudo skew-symmetric scheme [10].
In our previous studies, we proposed the hybrid scheme estimated by the Ducros sensor [11] , pseudo skew-symmetric scheme , and MUSCL-Roe scheme as follows:Here, the subscript denotes cell interface to evaluate the numerical flux along with -axis in this study. The sensor in (4) for switching the scheme is estimated by
Here, = 10−15 is a small value that prevents division by zero, and = 0.6 is the switching threshold. The vector is a velocity vector . The numerical fluxes for and in (4) are calculated by the formulas of (6) and (7):
Here a value is set to be 0.05 to suppress the spurious numerical oscillationwhere is the flux Jacobian and the subscript Roe denotes the Roe-averaged quantity of the left and right states. Herewhere and are the right and left eigenmatrices of , respectively, and is a diagonal matrix.
In the previous study [6], the symmetric central difference parts of (7) that are first and second terms of right hand side were replaced by the ones of the pseudo skew-symmetric form without losing the formal order of accuracy of the equation. However, in this study, the skew-symmetric form is not applied to (7) from the preliminary investigation. Detailed information for the investigation is described in Appendix A.
2.2. Immersed Boundary Method
A simple IBM based on the level set function is used to represent object boundaries in this study. The level set function is determined as a signed distance of minimal distance from the object surface. In the case of multiple objects, multiple level set functions are calculated based on simple minimum distance approach [12]. Whole cells are classified into three categories of fluid cells (FC), ghost cells (GC), and object cells (OC) by using the level set function expressed as (9). The ghost cells behave as guard cells between the fluid cells and object cells and are assigned for two layers under the present definition:
Figure 1 shows distributions of the classified cells around a cut sphere in the present IBM. The blue, brown, and green regions are assigned to the fluid, ghost, and object cells. The blue, brown, and green hemispheres are surfaces of the image points (), zero level set (), and the threshold of the object cells (). Here, the circle symbols of GC, IP, and IB represent a cell centroid of ghost cell, an image point, and an immersed boundary point. Developed numerical method in the previous study [6] in 2D is extended to 3D directly. Flow quantities of ghost cells are determined by using those of the image points. The flow variables of the image points at the edge of a “probe” are defined by trilinear interpolation from surrounding fluid cells. The probe is a vector from GC through IB to IP in the direction normal to the zero level set surface. The probe length expressed as in (10) is the fixed value of 1.75 times of the mesh size to avoid recursive reference in this study. The flow variables of the ghost cell are determined by assumed distributions on the probe. In this study, three velocity components of the ghost cell are calculated by the linear extrapolation while the density and the pressure are determined by Neumann condition along with the probe expressed as
2.3. Force Evaluation
Aerodynamic forces acting on the object surface as pressure and friction can be evaluated on the cell face between the fluid and the ghost cells. In the present method, a simple algorithm where the surface polygons are not necessary to estimate the forces by using staircase representation is employed. The blue and white cells correspond to the fluid and ghost cells in Figure 2. In algorithm A, an object boundary is considered as an assemblage of piecewise linear surfaces. Accordingly, cut faces should be determined to calculate the force components between the fluid and ghost cells in algorithm A. In algorithm B, on the other hand, the object is treated as a simple staircase agglomeration. Therefore, the cell faces become simple quadrangles because of the staircase representation in algorithm B. Consequently, the latter algorithm is quite simpler than the former algorithm because there is no need to calculate the precise surface area for the force evaluation. This becomes considerable especially in the case of moving boundary problems. The validity and applicability of the methodologies are described in Appendix B.
2.4. Sponge Region
The reflection of sound waves from outer boundaries may affect the flow in the subsonic case. Therefore, a numerical approach of the sponge layer proposed by Mani [13] and Freund [14] is employed to suppress the reflection. The method has several control parameters of and as the strength of the sponge layer depending on the flow field as follows:where , , , and are the damping ratio of the reflection, number of sponge layer cells, and parameters to control the strength of the sponge. The sponge strength should be adjusted based on the damping ratio.
2.5. Flowchart of the Flow Solver
Figure 3 shows a flowchart of the current computational processes. The level set function around a sphere is reconstructed due to the simple geometry. To solve a flow around 3D arbitrary geometry, convection and reinitialization of the level set function should be implemented to accomplish rapid computation.
3. Comparison of Drag Coefficient
3.1. Computational Condition
A drag coefficient around a sphere is compared between the current Cartesian and BFC solvers. All computational conditions are summarized in Table 1. Two kinds of numerical schemes for evaluating the inviscid flux are as follows: (A) the developed hybrid scheme and (B) the MUSCL-Roe scheme. Reynolds number based on the sphere diameter and freestream values is fixed at 300, while the freestream Mach numbers are set to be 0.3 and 1.2 as subsonic and supersonic flows. Two kinds of mesh resolutions 0.1D and 0.05D are prepared for the comparison. Total numbers of computational cells are 2 million (200 × 100 × 100) for 0.1D and 16 million (400 × 200 × 200) for 0.05D. The schematic of a computational domain is shown in Figure 4. All flow variables are specified by Dirichlet boundary conditions at the inflow boundary and on density alone at the outflow boundary only for the case of subsonic flow. Neumann conditions are used at other boundaries for all variables. The sponge regions are arranged next to the outer boundaries with ten layers especially in the subsonic cases. The sphere geometry is represented by the Cartesian mesh as shown in Figure 5. The black solid lines are grid lines connecting cell centers. The bold red line is the boundary of a sphere immersed in the Cartesian mesh. The blue, white, and orange regions show fluid, ghost, and object cells, respectively. Table 2 shows the boundary layer thickness . From the estimation, just one cell may be allocated in the boundary layer even in the case with the fine mesh.
(a) 0.1D
(b) 0.05D
The BFC solver is based on WENOCU6 method developed by Nonomura et al. [7, 15]. The number of grid points of the BFC simulation is 107 × 48 × 177 as shown in Figure 20 and minimum grid resolution is 6.5 × 10−3D.
3.2. Computational Result at Mach 0.3
Instantaneous Mach number distributions at Mach 0.3 are obtained as shown in Figure 6. Even though asymmetric flow appears like in Figure 6(a) in the BFC result, a symmetrical flow is observed in the results of the MUSCL scheme with the coarse mesh as shown in Figure 6(c). When the fine mesh is used, weak asymmetric flow is observed even in the MUSCL scheme in Figure 6(e). With regard to the proposed scheme, however, asymmetric flows can be observed even in the coarse mesh in Figure 6(b). In the case of the fine mesh with the proposed scheme, the flow feature is quite similar to one obtained by the BFC solver, as shown in Figure 6(d). As a result, qualitative agreements are confirmed from the visualization.
(a) BFC
(b) M03D010A
(c) M03D010B
(d) M03D005A
(e) M03D005B
In the previous 2D study [6], we confirmed that the hybrid scheme worked well by visualizing the spatial distributions of the switching parameter . Figure 7 shows instantaneous distributions at center planes and isosurface of . In the previous study, we could observe that the MUSCL scheme corresponding to the white region was wrongly selected due to the acoustic wave reflection from outer boundaries. In this study, however, the MUSCL scheme is appropriately utilized for the specified region since the unphysical reflection is suitably suppressed owing to the sponge region. In the case of the fine mesh, the regions solved by the MUSCL and pseudo skew-symmetric schemes are clearly separated. The white region for the MUSCL scheme can be spread according to acoustic dipole from the sphere. From both 2D images of Figure 7, the wake of the sphere is covered by the pseudo skew-symmetric scheme, so that the asymmetric flow can be produced by the low dissipation nature of the scheme.
(a) M03D010A
(b) M03D005A
The wake structure of a sphere at Mach 0.3 is visualized by second invariant of velocity tensors (-criterion) in Figure 8. Hairpin type vortices are obtained by the BFC solver as shown in Figure 8(a) and the proposed hybrid scheme as shown in Figures 8(b) and 8(d). In the case of the MUSCL scheme, however, the wake structure was not captured due to large numerical dissipation even in the fine mesh.
(a) BFC
(b) M03D010A
(c) M03D010B
(d) M03D005A
(e) M03D005B
Figure 9 summarizes time histories of drag and lift coefficients. Results from only the MUSCL scheme show constant values in contrast with the ones from the hybrid scheme. From the results, the hybrid scheme has an advantage in capturing the unsteadiness that cannot be captured by the MUSCL scheme:
Figure 10 shows time-averaged drag coefficients, where black lines, red squares, and blue circles are results from the BFC solver, hybrid scheme (A), and only MUSCL scheme (B), respectively. In Figure 10(a), theoretically predicted value from Carlson’s formula [8] in (12) is presented for the comparison. Both the total drag coefficients and the friction drag coefficients of the current results were underestimated while the pressure drag coefficients were overestimated. This could be caused from lack of mesh resolution based on minimum mesh sizes to resolve the boundary layer thickness described in Table 2. The drag coefficients of (A) and (B) showed almost the same values even though the wake structures were quite different as denoted above. Consequently, the proposed hybrid scheme is appropriate to capture both drag coefficient and wake simultaneously.
(a) Total drag coefficient
(b) Pressure drag coefficient
(c) Friction drag coefficient
Flow fields around a sphere are solved at freestream Mach number 1.2 to confirm the validity of the current numerical method for supersonic conditions. In the visualization of Figure 11, the distortion of the density contours is confirmed in the present results owing to the smaller computational domain compared with the one for the BFC solver. With regard to the wake structure, there is a little difference between the hybrid and the MUSCL schemes. The wake solved by the MUSCL scheme demonstrates more dissipative feature than the one of the hybrid scheme. From the visualization of the switching parameter in Figure 12, we can observe that the wake was solved by the pseudo skew-symmetric scheme. The low dissipation characteristics of the hybrid scheme work well to keep the flow structure. Even in the coarse mesh computation, the distribution of the switching parameter does not show big change. This indicates that the flow structure can be captured by using the mesh size at this flow condition.
(a) BFC
(b) M12D010A
(c) M12D010B
(d) M12D005A
(e) M12D005B
(a) M12D010A
(b) M12D005A
Time-averaged drag coefficients at Mach 1.2 in Figure 13 show the same trends as those at the Mach number of 0.3. The pressure drag coefficients are overestimated while the friction drag coefficients are underestimated.
(a) Total drag coefficient
(b) Pressure drag coefficient
(c) Friction drag coefficient
Consequently, the accuracy of the present flow solver was confirmed by the discrepancy with the drag coefficient from the high accurate BFC flow solver. All results could be categorized into “A” (less than 5 percent discrepancy) and “B” (5–10 percent discrepancy) in Table 3. In spite of the low mesh resolution for the boundary layer thickness, this flow solver demonstrates good applicability to predict the total drag coefficient around a sphere in the flow condition.
4. Flow Simulation around a Moving Sphere
4.1. Computational Condition
To solve a real gas-particle flow at the rocket launch, we have to confirm the validity of a flow around moving particles. Therefore, next validation is conducted to compare drag coefficient around a fixed sphere and a moving sphere. In applying the sharp interface IBM to the moving boundary problem, the so-called “fresh cell” problem [16, 17] can occur due to sudden appearance/disappearance of the fluid/ghost cell. This problem can be suppressed by applying cut-cell like approach [17]. In this simulation, however, no special treatment is implemented to keep simple IBM to solve the flow around moving multiple particles. Therefore, the accuracy of the present flow simulation should be investigated by the flow around a moving object.
Computational conditions are summarized in Table 4. The Reynolds number based on the diameter of a sphere is fixed to be 300 while the Mach numbers based on the moving velocities are set to be 0.3 and 1.2. Mesh resolutions are 0.1D and 0.05D. The number of grid points is 200 × 100 × 100 for 0.1D and 400 × 200 × 200 for 0.05D cases. The schematic of a computational domain is drawn in Figure 14. The boundary conditions are the same as the previous fixed sphere cases. The first capital characters of F and M declare the case of a fixed sphere in uniform freestream and a moving sphere in static flow, respectively.
4.2. Computational Result
Figure 15 demonstrates density contours around a moving sphere at Mach 0.3 and Mach 1.2 at the beginning of the motion. The flow visualization around a moving sphere shows qualitative agreements with the one of the fixed sphere as shown in Figures 6 and 11. The flow structure observed around a fixed sphere can be captured around a moving sphere.
(a) = 1.6
(b) = 5.0
(c) = 8.3
A developed Cartesian solver underestimated the friction drag and overestimated the pressure drag for the fixed sphere analysis as shown in Figures 10 and 13. In the case of the moving sphere at Mach 0.3 and Mach 1.2, however, opposite trends to the fixed cases were observed in both the pressure and the friction drag coefficients in Figures 16 and 17. This trend was also confirmed for the 2D cylinder as shown in Appendix B. Furthermore, the pressure and friction drag coefficients at Mach 0.3 were scattered as shown in Figure 16. In the case of the 2D cylinder, on the other hand, the values of the friction drag coefficients were almost identical at both Mach numbers. Therefore, this scattering could occur from the small computational domain along with the -axis direction. Moreover, high frequency oscillations of the drag coefficients were observed in Figures 16 and 17 especially in the pressure drag coefficients of the coarse mesh resolution. Although this pressure oscillation from the fresh cell problem is a drawback of the present simple IBM which contains no special treatments for the conservation law, it could be suppressed by using finer mesh resolution of 0.05D. In a flow simulation of gas-particle flow involving multiple particles and shocks, the numerical oscillation may be relatively weakened due to the high-speed flow around particles. Therefore, it is considered that the developed numerical method is applicable for a flow simulation of the gas-particle flow with moving particles.
5. Conclusion
A simple IBM for compressible flows was developed and applied to flows around a stationary and a moving sphere. To solve gas-particle flows consisting of shocks and turbulence, a hybrid scheme was employed with the simple IBM. Moreover, a simple and rapid algorithm for force evaluation was implemented and validated. Obtained drag coefficients around a sphere were compared with the ones from high accurate BFC solver. As a result, we could confirm several features of the present numerical method as follows:
Unsteady flow characteristics such as the wake structure and time variance of drag coefficient at Mach 0.3 could be captured by the hybrid scheme. On the other hand, the MUSCL-Roe scheme was not suitable to resolve the flow unsteadiness due to the large numerical dissipation.
The total drag coefficient of a fixed sphere showed good agreement with the one from the BFC solver. The pressure and friction drag components, however, were over- and underestimated in the present IBM by the current mesh resolutions in both subsonic and supersonic flows.
The total drag coefficient of a moving sphere could show grid convergence as same as the case of a fixed sphere. The pressure and the friction drag components, however, showed reversal trends with the results of a fixed sphere.
The numerical oscillation of the drag coefficients caused from the fresh cell problem could be suppressed by the fine mesh resolution. Consequently, we could confirm the applicability of the developed numerical method to a flow simulation of gas-particle flow with moving particles.
Next issue for the gas-particle flow simulation is to treat the momentum and heat exchange between particles and a flow by coupled simulation.
Appendices
A. Symmetric Central Difference Parts in the Hybrid Scheme
We have two options for the symmetric central difference parts of (7) as original terms or modified terms of used in the previous study [6]. Although the latter showed good performance in the previous study, the former is more straightforward expression than the latter. Therefore, drag coefficient of a fixed sphere was compared between these schemes at Mach 0.3 and Mach 1.2. The results obtained from former (A), latter (C), and only MUSCL scheme (B) are summarized in Table 5. The values show discrepancies between the evaluated drag coefficients and the reference value obtained by the highly accurate BFC solver. The original MUSCL-Roe (A) showed smaller discrepancies than the previous scheme (C) as shown in Table 5. Consequently, scheme (A) consisting of usual MUSCL-Roe fluxes was suitable for the symmetric central difference parts in this flow regime.
B. Aerodynamic Force Evaluation without Surface Polygons
To realize large-scale computations by high performance computers, simple and easy programs are preferable from the aspect of extensibility and rapid modification. The aerodynamic force evaluation in IBM with Cartesian mesh is one of the concerns due to the implicit boundary representation. The object boundary is approximated as assemblage of planes in 3D or lines in 2D in computational cells in the present IBM. Although it is straightforward procedure that the aerodynamic force is estimated by using the cut planes, it may become time consuming work especially for a moving boundary problem. Nonomura et al. proposed simple force measurement technique for arbitrary geometry based on the staircase representation. We applied the technique to a fixed and moving object and numerical tests were conducted as shown in Table 6. Algorithms “A” and “B” indicate schematics of Figure 2. Algorithm A is a force calculation method by using cut planes between the ghost and fluid cells. Algorithm B, on the other hand, is based on the staircase object representation proposed by Nonomura et al. The object geometry is a 2D cylinder in these cases. The cylinder moves along with -axis. The computational domain sizes are 70D and 10D along with - and -axes based on the diameter of the cylinder D.
From the time-averaged drag coefficient at Mach 0.3 (Figure 18) and Mach 1.2 (Figure 19), trends of the grid convergence could be observed for both force evaluation methods even though opposite tendencies were also clarified between the cases of a fixed and a moving cylinder. From the comparison of algorithms A and B, reasonable prediction of the aerodynamic forces could be accomplished even in using simpler algorithm B. Therefore, we employed algorithm B to estimate drag coefficient for all computations in this paper.
(a) Total drag coefficient
(b) Pressure drag coefficient
(c) Friction drag coefficient
(a) Total drag coefficient
(b) Pressure drag coefficient
(c) Friction drag coefficient
C. Computational Mesh for Boundary-Fitted Coordinate Solver
See Figure 20.
Conflict of Interests
The authors declare no conflict of interests.
Acknowledgment
Part of the computations in this study was conducted by NEC SX-9 supported by HPCI (High Performance Computing Infrastructure) hp140138 of the Research Organization for Information Sciences and Technology.