# A Simple Immersed Boundary Method for Compressible Flow Simulation around a Stationary and Moving Sphere

Authors: Yusuke Mizuno; Shun Takahashi; Taku Nonomura; Takayuki Nagata; Kota Fukuda
Journal: Mathematical Problems in Engineering
 (2015)

---

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


---

## Body

## 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:(1)
∂
Q
∂
t
+
∂
E
∂
x
+
∂
F
∂
y
+
∂
G
∂
z
=
∂
E
v
∂
x
+
∂
F
v
∂
y
+
∂
G
v
∂
z
,
Q
=
ρ
ρ
u
ρ
v
ρ
w
ρ
e
,
E
=
ρ
u
ρ
u
2
+
p
ρ
v
u
ρ
w
u
ρ
e
+
p
u
,
F
=
ρ
v
ρ
u
v
ρ
v
2
+
p
ρ
w
v
ρ
e
+
p
v
,
G
=
ρ
w
ρ
u
w
ρ
v
w
ρ
w
2
+
p
ρ
e
+
p
w
,
E
v
=
0
τ
x
x
τ
x
y
τ
x
z
τ
x
x
u
+
τ
x
y
v
+
τ
x
z
w
+
κ
T
x
,
F
v
=
0
τ
y
x
τ
y
y
τ
y
z
τ
y
x
u
+
τ
y
y
v
+
τ
y
z
w
+
κ
T
y
,
G
v
=
0
τ
z
x
τ
z
y
τ
z
z
τ
z
x
u
+
τ
z
y
v
+
τ
z
z
w
+
κ
T
z
,where E, F, and G are inviscid fluxes in the x
-, y
-, and z
-directions, E
v, F
v, and G
v are viscous fluxes, and Q contains conservative variables. The pressure p is related to the total energy per unit mass e by an equation of state:(2)
ρ
e
=
p
γ
-
1
+
1
2
ρ
u
2
+
v
2
+
w
2
.Here, viscous stress tensor components are given as(3)
τ
x
x
=
2
3
μ
2
u
x
-
v
y
-
w
z
,
τ
y
y
=
2
3
μ
2
v
y
-
w
z
-
u
x
,
τ
z
z
=
2
3
μ
2
w
z
-
u
x
-
v
y
,
τ
x
y
=
τ
y
x
=
μ
u
y
+
v
x
,
τ
y
z
=
τ
z
y
=
μ
v
z
+
w
y
,
τ
z
x
=
τ
x
z
=
μ
w
x
+
u
z
.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 E
t
u
r
b
u
l
e
n
t, and MUSCL-Roe scheme E
s
h
o
c
k as follows:(4)
E
i
+
1
/
2
=
1
-
β
i
+
1
/
2
E
i
+
1
/
2
,
t
u
r
b
u
l
e
n
t
+
β
i
+
1
/
2
E
i
+
1
/
2
,
s
h
o
c
k
.Here, the subscript i
+
1
/
2 denotes cell interface to evaluate the numerical flux along with x-axis in this study. The sensor β in (4) for switching the scheme is estimated by(5)
β
i
+
1
/
2
=
min
⁡
1
,
ϕ
i
+
ϕ
i
+
1
,
ϕ
=
0
,
∇
·
u
2
∇
·
u
2
+
∇
×
u
2
+
ε
<
θ
,
1
,
∇
·
u
2
∇
·
u
2
+
∇
×
u
2
+
ε
≥
θ
.Here,ε = 10−15 is a small value that prevents division by zero, and θ = 0.6 is the switching threshold. The vector u is a velocity vector u
=
(
u
,
v
,
w
). The numerical fluxes for E
t
u
r
b
u
l
e
n
t and E
s
h
o
c
k in (4) are calculated by the formulas of (6) and (7):(6)
E
i
+
1
/
2
,
t
u
r
b
u
l
e
n
t
=
E
i
+
1
/
2
,
c
e
n
t
-
σ
i
+
1
/
2
,
c
e
n
t
,
E
i
+
1
/
2
,
c
e
n
t
=
1
2
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
1
4
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
u
i
+
u
i
+
1
+
1
2
p
i
+
p
i
+
1
1
4
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
v
i
+
v
i
+
1
1
4
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
w
i
+
w
i
+
1
1
4
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
e
i
+
p
i
ρ
i
+
e
i
+
1
+
p
i
+
1
ρ
i
+
1
,
σ
i
+
1
/
2
,
c
e
n
t
=
1
4
χ
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
0
u
i
+
1
/
2
,
R
-
u
i
+
1
/
2
,
L
v
i
+
1
/
2
,
R
-
v
i
+
1
/
2
,
L
w
i
+
1
/
2
,
R
-
w
i
+
1
/
2
,
L
e
i
+
1
/
2
,
R
-
e
i
+
1
/
2
,
L
.Here a valueχ is set to be 0.05 to suppress the spurious numerical oscillation(7)
E
i
+
1
/
2
,
s
h
o
c
k
=
1
2
E
i
+
E
i
+
1
-
A
i
+
1
/
2
,
R
o
e
Q
i
+
1
/
2
,
R
-
Q
i
+
1
/
2
,
L
,where A is the flux Jacobian and the subscript Roe denotes the Roe-averaged quantity of the left and right states. Here(8)
A
=
R
Λ
L
,where R and L are the right and left eigenmatrices of A, respectively, and Λ
=
L
A
R 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:(9)
φ
O
C
<
-
2
√
3
Δ
x
,
-
2
√
3
Δ
x
≤
φ
G
C
≤
0
,
0
<
φ
F
C
.Figure1 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 (φ
=
1.75
Δ
x), zero level set (φ
=
0), and the threshold of the object cells (φ
=
-
2
√
3
Δ
x). 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 d
I
P 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 (10)
ρ
G
C
=
ρ
I
P
,
u
G
C
=
u
I
P
-
d
I
P
+
φ
G
C
d
I
P
u
I
P
-
u
I
B
,
v
G
C
=
v
I
P
-
d
I
P
+
φ
G
C
d
I
P
v
I
P
-
v
I
B
,
w
G
C
=
w
I
P
-
d
I
P
+
φ
G
C
d
I
P
w
I
P
-
w
I
B
,
p
G
C
=
p
I
P
.Figure 1
Schematic of the present IBM for a sphere.


### 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 Figure2. 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.Figure 2
Representation of sphere for force evaluation.


### 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 x
s
p as the strength of the sponge layer depending on the flow field as follows:(11)
η
t
a
r
g
e
t
=
-
20
2
log
⁡
e
1
-
M
2
∫
0
x
s
p
σ
x
d
x
,
σ
x
=
σ
s
p
x
s
p
-
x
x
s
p
β
,
x
s
p
=
Δ
x
·
L
s
p
,
η
t
a
r
g
e
t
=
-
20
2
log
⁡
e
1
-
M
2
∫
0
x
s
p
σ
s
p
x
s
p
-
x
x
s
p
β
d
x
=
-
20
2
log
⁡
e
1
-
M
2
σ
s
p
x
s
p
β
+
1
,where η
t
a
r
g
e
t, L
s
p, σ
s
p, 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

Figure3 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.Figure 3
Flowchart.


## 2.1. Governing Equations

The governing equations of the present flow simulations are three-dimensional compressible Navier-Stokes equations:(1)
∂
Q
∂
t
+
∂
E
∂
x
+
∂
F
∂
y
+
∂
G
∂
z
=
∂
E
v
∂
x
+
∂
F
v
∂
y
+
∂
G
v
∂
z
,
Q
=
ρ
ρ
u
ρ
v
ρ
w
ρ
e
,
E
=
ρ
u
ρ
u
2
+
p
ρ
v
u
ρ
w
u
ρ
e
+
p
u
,
F
=
ρ
v
ρ
u
v
ρ
v
2
+
p
ρ
w
v
ρ
e
+
p
v
,
G
=
ρ
w
ρ
u
w
ρ
v
w
ρ
w
2
+
p
ρ
e
+
p
w
,
E
v
=
0
τ
x
x
τ
x
y
τ
x
z
τ
x
x
u
+
τ
x
y
v
+
τ
x
z
w
+
κ
T
x
,
F
v
=
0
τ
y
x
τ
y
y
τ
y
z
τ
y
x
u
+
τ
y
y
v
+
τ
y
z
w
+
κ
T
y
,
G
v
=
0
τ
z
x
τ
z
y
τ
z
z
τ
z
x
u
+
τ
z
y
v
+
τ
z
z
w
+
κ
T
z
,where E, F, and G are inviscid fluxes in the x
-, y
-, and z
-directions, E
v, F
v, and G
v are viscous fluxes, and Q contains conservative variables. The pressure p is related to the total energy per unit mass e by an equation of state:(2)
ρ
e
=
p
γ
-
1
+
1
2
ρ
u
2
+
v
2
+
w
2
.Here, viscous stress tensor components are given as(3)
τ
x
x
=
2
3
μ
2
u
x
-
v
y
-
w
z
,
τ
y
y
=
2
3
μ
2
v
y
-
w
z
-
u
x
,
τ
z
z
=
2
3
μ
2
w
z
-
u
x
-
v
y
,
τ
x
y
=
τ
y
x
=
μ
u
y
+
v
x
,
τ
y
z
=
τ
z
y
=
μ
v
z
+
w
y
,
τ
z
x
=
τ
x
z
=
μ
w
x
+
u
z
.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 E
t
u
r
b
u
l
e
n
t, and MUSCL-Roe scheme E
s
h
o
c
k as follows:(4)
E
i
+
1
/
2
=
1
-
β
i
+
1
/
2
E
i
+
1
/
2
,
t
u
r
b
u
l
e
n
t
+
β
i
+
1
/
2
E
i
+
1
/
2
,
s
h
o
c
k
.Here, the subscript i
+
1
/
2 denotes cell interface to evaluate the numerical flux along with x-axis in this study. The sensor β in (4) for switching the scheme is estimated by(5)
β
i
+
1
/
2
=
min
⁡
1
,
ϕ
i
+
ϕ
i
+
1
,
ϕ
=
0
,
∇
·
u
2
∇
·
u
2
+
∇
×
u
2
+
ε
<
θ
,
1
,
∇
·
u
2
∇
·
u
2
+
∇
×
u
2
+
ε
≥
θ
.Here,ε = 10−15 is a small value that prevents division by zero, and θ = 0.6 is the switching threshold. The vector u is a velocity vector u
=
(
u
,
v
,
w
). The numerical fluxes for E
t
u
r
b
u
l
e
n
t and E
s
h
o
c
k in (4) are calculated by the formulas of (6) and (7):(6)
E
i
+
1
/
2
,
t
u
r
b
u
l
e
n
t
=
E
i
+
1
/
2
,
c
e
n
t
-
σ
i
+
1
/
2
,
c
e
n
t
,
E
i
+
1
/
2
,
c
e
n
t
=
1
2
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
1
4
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
u
i
+
u
i
+
1
+
1
2
p
i
+
p
i
+
1
1
4
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
v
i
+
v
i
+
1
1
4
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
w
i
+
w
i
+
1
1
4
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
e
i
+
p
i
ρ
i
+
e
i
+
1
+
p
i
+
1
ρ
i
+
1
,
σ
i
+
1
/
2
,
c
e
n
t
=
1
4
χ
ρ
i
u
i
+
ρ
i
+
1
u
i
+
1
0
u
i
+
1
/
2
,
R
-
u
i
+
1
/
2
,
L
v
i
+
1
/
2
,
R
-
v
i
+
1
/
2
,
L
w
i
+
1
/
2
,
R
-
w
i
+
1
/
2
,
L
e
i
+
1
/
2
,
R
-
e
i
+
1
/
2
,
L
.Here a valueχ is set to be 0.05 to suppress the spurious numerical oscillation(7)
E
i
+
1
/
2
,
s
h
o
c
k
=
1
2
E
i
+
E
i
+
1
-
A
i
+
1
/
2
,
R
o
e
Q
i
+
1
/
2
,
R
-
Q
i
+
1
/
2
,
L
,where A is the flux Jacobian and the subscript Roe denotes the Roe-averaged quantity of the left and right states. Here(8)
A
=
R
Λ
L
,where R and L are the right and left eigenmatrices of A, respectively, and Λ
=
L
A
R 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:(9)
φ
O
C
<
-
2
√
3
Δ
x
,
-
2
√
3
Δ
x
≤
φ
G
C
≤
0
,
0
<
φ
F
C
.Figure1 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 (φ
=
1.75
Δ
x), zero level set (φ
=
0), and the threshold of the object cells (φ
=
-
2
√
3
Δ
x). 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 d
I
P 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 (10)
ρ
G
C
=
ρ
I
P
,
u
G
C
=
u
I
P
-
d
I
P
+
φ
G
C
d
I
P
u
I
P
-
u
I
B
,
v
G
C
=
v
I
P
-
d
I
P
+
φ
G
C
d
I
P
v
I
P
-
v
I
B
,
w
G
C
=
w
I
P
-
d
I
P
+
φ
G
C
d
I
P
w
I
P
-
w
I
B
,
p
G
C
=
p
I
P
.Figure 1
Schematic of the present IBM for a sphere.


## 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 Figure2. 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.Figure 2
Representation of sphere for force evaluation.


## 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 x
s
p as the strength of the sponge layer depending on the flow field as follows:(11)
η
t
a
r
g
e
t
=
-
20
2
log
⁡
e
1
-
M
2
∫
0
x
s
p
σ
x
d
x
,
σ
x
=
σ
s
p
x
s
p
-
x
x
s
p
β
,
x
s
p
=
Δ
x
·
L
s
p
,
η
t
a
r
g
e
t
=
-
20
2
log
⁡
e
1
-
M
2
∫
0
x
s
p
σ
s
p
x
s
p
-
x
x
s
p
β
d
x
=
-
20
2
log
⁡
e
1
-
M
2
σ
s
p
x
s
p
β
+
1
,where η
t
a
r
g
e
t, L
s
p, σ
s
p, 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

Figure3 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.Figure 3
Flowchart.


## 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 Table1. 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.Table 1
Test cases.
Re
Uniform flow
Mesh size
Scheme
Case
300
0.3
0.10D
A
M03D010A
0.05D
M03D005A
0.10D
B
M03D010B
0.05D
M03D005B
1.2
0.10D
A
M12D010A
0.05D
M12D005A
0.10D
B
M12D010B
0.05D
M12D005BTable 2
Boundary layer thickness.
Re
δ
Mesh size
300
0.058D
0.10D
0.05DFigure 4
Computational domain (M03D010).Figure 5
Sphere geometry in Cartesian mesh with IBM.
(a)
0.1D
(b)
0.05DThe 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−3
D.


### 3.2. Computational Result at Mach 0.3

Instantaneous Mach number distributions at Mach 0.3 are obtained as shown in Figure6. 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.Figure 6
Mach number distributions at Mach 0.3.
(a)
BFC
(b)
M03D010A
(c)
M03D010B
(d)
M03D005A
(e)
M03D005BIn 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 β
=
1. 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.Figure 7
Instantaneous distributions of switching parameterβ (black: β
=
0, white: β
=
1) at Mach 0.3 (left: 2D centerplane, right: 3D isosurface at β
=
1).
(a)
M03D010A
(b)
M03D005AThe wake structure of a sphere at Mach 0.3 is visualized by second invariant of velocity tensors (Q-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.Figure 8
Isosurfaces of second invariant of velocity gradient tensor (Q = 1.0 × 10−5) colored by Mach number.
(a)
BFC
(b)
M03D010A
(c)
M03D010B
(d)
M03D005A
(e)
M03D005BFigure9 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:(12)
C
D
=
24
R
e
1
+
0.15
R
e
0.687
1
+
exp
⁡
-
0.427
/
M
4.63
-
3.0
/
R
e
0.88
1
+
M
/
R
e
3.82
+
1.28
e
x
p
-
1.25
R
e
/
M
.Figure 9
Time histories of drag and lift coefficients at Mach 0.3.Figure10 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.Figure 10
Time-averaged total, pressure, and friction drag coefficients at Mach 0.3.
(a)
Total drag coefficient
(b)
Pressure drag coefficient
(c)
Friction drag coefficientFlow 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 Figure11, 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.Figure 11
Density distributions at Mach 1.2.
(a)
BFC
(b)
M12D010A
(c)
M12D010B
(d)
M12D005A
(e)
M12D005BFigure 12
Instantaneous distributions of switching parameterβ (black: β
=
0, white: β
=
1) at Mach 1.2 (left: 2D centerplane, right: 3D isosurface at β
=
1).
(a)
M12D010A
(b)
M12D005ATime-averaged drag coefficients at Mach 1.2 in Figure13 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.Figure 13
Time-averaged total, pressure, and friction drag coefficients at Mach 1.2.
(a)
Total drag coefficient
(b)
Pressure drag coefficient
(c)
Friction drag coefficientConsequently, 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 Table3. 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.Table 3
Performance of drag prediction.
Case
Performance
M03D010A
B
M03D010B
B
M03D005A
A
M03D005B
A
M12D010A
B
M12D010B
B
M12D005A
A
M12D005B
A


## 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 Table1. 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.Table 1
Test cases.
Re
Uniform flow
Mesh size
Scheme
Case
300
0.3
0.10D
A
M03D010A
0.05D
M03D005A
0.10D
B
M03D010B
0.05D
M03D005B
1.2
0.10D
A
M12D010A
0.05D
M12D005A
0.10D
B
M12D010B
0.05D
M12D005BTable 2
Boundary layer thickness.
Re
δ
Mesh size
300
0.058D
0.10D
0.05DFigure 4
Computational domain (M03D010).Figure 5
Sphere geometry in Cartesian mesh with IBM.
(a)
0.1D
(b)
0.05DThe 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−3
D.


## 3.2. Computational Result at Mach 0.3

Instantaneous Mach number distributions at Mach 0.3 are obtained as shown in Figure6. 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.Figure 6
Mach number distributions at Mach 0.3.
(a)
BFC
(b)
M03D010A
(c)
M03D010B
(d)
M03D005A
(e)
M03D005BIn 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 β
=
1. 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.Figure 7
Instantaneous distributions of switching parameterβ (black: β
=
0, white: β
=
1) at Mach 0.3 (left: 2D centerplane, right: 3D isosurface at β
=
1).
(a)
M03D010A
(b)
M03D005AThe wake structure of a sphere at Mach 0.3 is visualized by second invariant of velocity tensors (Q-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.Figure 8
Isosurfaces of second invariant of velocity gradient tensor (Q = 1.0 × 10−5) colored by Mach number.
(a)
BFC
(b)
M03D010A
(c)
M03D010B
(d)
M03D005A
(e)
M03D005BFigure9 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:(12)
C
D
=
24
R
e
1
+
0.15
R
e
0.687
1
+
exp
⁡
-
0.427
/
M
4.63
-
3.0
/
R
e
0.88
1
+
M
/
R
e
3.82
+
1.28
e
x
p
-
1.25
R
e
/
M
.Figure 9
Time histories of drag and lift coefficients at Mach 0.3.Figure10 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.Figure 10
Time-averaged total, pressure, and friction drag coefficients at Mach 0.3.
(a)
Total drag coefficient
(b)
Pressure drag coefficient
(c)
Friction drag coefficientFlow 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 Figure11, 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.Figure 11
Density distributions at Mach 1.2.
(a)
BFC
(b)
M12D010A
(c)
M12D010B
(d)
M12D005A
(e)
M12D005BFigure 12
Instantaneous distributions of switching parameterβ (black: β
=
0, white: β
=
1) at Mach 1.2 (left: 2D centerplane, right: 3D isosurface at β
=
1).
(a)
M12D010A
(b)
M12D005ATime-averaged drag coefficients at Mach 1.2 in Figure13 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.Figure 13
Time-averaged total, pressure, and friction drag coefficients at Mach 1.2.
(a)
Total drag coefficient
(b)
Pressure drag coefficient
(c)
Friction drag coefficientConsequently, 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 Table3. 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.Table 3
Performance of drag prediction.
Case
Performance
M03D010A
B
M03D010B
B
M03D005A
A
M03D005B
A
M12D010A
B
M12D010B
B
M12D005A
A
M12D005B
A


## 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 Table4. 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.Table 4
Test cases.
Re
Velocity
Mesh size
State
Case
300
0.3
0.10D
Fix
M03D010F
0.05D
M03D005F
0.10D
Move
M03D010M
0.05D
M03D005M
1.2
0.10D
Fix
M12D010F
0.05D
M12D005F
0.10D
Move
M12D010M
0.05D
M12D005MFigure 14
Computational domain for moving sphere.


### 4.2. Computational Result

Figure15 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.Figure 15
Initial density contours around a moving sphere (left: Mach 0.3, right: Mach 1.2).
(a)
t = 1.6
(b)
t = 5.0
(c)
t = 8.3A developed Cartesian solver underestimated the friction drag and overestimated the pressure drag for the fixed sphere analysis as shown in Figures10 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 x-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.Figure 16
Time history of drag coefficient for fixed/moving sphere at Mach 0.3 Relative Velocity.Figure 17
Time history of drag coefficient for fixed/moving sphere at Mach 1.2 Relative Velocity.


## 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 Table4. 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.Table 4
Test cases.
Re
Velocity
Mesh size
State
Case
300
0.3
0.10D
Fix
M03D010F
0.05D
M03D005F
0.10D
Move
M03D010M
0.05D
M03D005M
1.2
0.10D
Fix
M12D010F
0.05D
M12D005F
0.10D
Move
M12D010M
0.05D
M12D005MFigure 14
Computational domain for moving sphere.


## 4.2. Computational Result

Figure15 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.Figure 15
Initial density contours around a moving sphere (left: Mach 0.3, right: Mach 1.2).
(a)
t = 1.6
(b)
t = 5.0
(c)
t = 8.3A developed Cartesian solver underestimated the friction drag and overestimated the pressure drag for the fixed sphere analysis as shown in Figures10 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 x-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.Figure 16
Time history of drag coefficient for fixed/moving sphere at Mach 0.3 Relative Velocity.Figure 17
Time history of drag coefficient for fixed/moving sphere at Mach 1.2 Relative Velocity.


## 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:(
1
) 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.(
2
) 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.(
3
) 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.(
4
) 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.



---

*Source: 438086-2015-10-22.xml*