## 3 Results  
### Validation around a Fixed Particle  
The accuracy of the flow solver was determined by comparing the calculated flow structure and the drag coefficient around a fixed sphere with those of previous studies and drag models [21; 22; 23; 24]. The convection term was evaluated using two schemes: U, the first-order upwind finite-difference scheme and S, the second-order skew-symmetric scheme. The Reynolds numbers was set to \(Re\) = 300 and 400 depending on the freestream and particle diameter values. A comparison was made using two mesh resolutions 0.05\(D\) (fine) and 0.1\(D\) (coarse), where \(D\) is the particle diameter. The computational conditions are summarized in Table 1. The computational domain (see Figure 3) is 20\(D\)\(\times\) 10\(D\)\(\times\) 10\(D\). The boundary values of the velocity and pressure were defined by Dirichlet boundary conditions at the inlet and outlet boundaries, and all variables were subject to Neumann conditions at the other boundaries. All variables were normalized by the freestream values of the density, velocity and the reference length (here taken as \(D\)). Figure 4 plots the mean drag coefficient versus Reynolds number. The drag models of [PERSON] et al. [21] and [PERSON] et al. [22] are respectively:  
\[C_{D}=\frac{24}{Re}\bigg{(}1+0.15 Re^{0.687}+\frac{0.0175 Re}{1+4.25\times 10^{4} Re^{-1.16}}\bigg{)}, \tag{12}\] \begin{table}
\begin{tabular}{c c c c} \hline \hline
**Reynolds Number** & **Mesh Size** & **Scheme** & **Case** \\ \hline \multirow{3}{*}{300} & 0.10\(D\) & & Re300D010U \\  & 0.05\(D\) & upwind. & Re300D05U \\  & 0.10\(D\) & & Re300D010S \\ \hline \multirow{3}{*}{400} & 0.05\(D\) & & Re300D05S \\  & 0.10\(D\) & skew-sym. & Re400D010S \\ \cline{1-1}  & 0.05\(D\) & & Re400D005S \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Test cases in the proposed model. Figure 3: Computational domain. The red sphere is the object. \[C_{D}=\frac{24}{Re}\Big{(}1+0.15 Re^{0.687}\Big{)}. \tag{13}\]  
To show the adaptability of the scheme, the drag coefficients in schemes U and S were compared at \(Re=300\). The difference between these two schemes became small at high grid resolutions. Moreover, the difference between scheme S and the previous models was smaller than that of scheme U. The present results captured the trends produced by previous drag models and the drag coefficients obtained agreed well with those of previous numerical results. Figures 5 and 6 summarize the nondimensional time histories of the drag and the lift coefficients at \(Re=300\), respectively. The coefficients oscillated in all cases except Re300D010A. The amplitudes of the oscillations were higher for a fine-mesh than for a coarse-mesh and were almost twice as high in scheme S than in scheme U. From the distributions of the time-averaged pressure coefficients on the center plane (Figure 7); the time-average was processed in five cycles of the steady wake vortex. We observed the same non-axisymmetric flows as in a previous study [23]. We saw a weak vortex in the wake of scheme U instead of scheme S because of numerical dissipation. Furthermore, a smooth contour distribution was established as the grid resolution increased. The instantaneous wake structure was visualized using the second invariant of the velocity tensors \(Q\) (Q-criterion) in Figure 8. At \(Re=300\), [PERSON] et al. [23] observed a hairpin vortex in the wake of the particle. In the present study, all cases except Re300D010A showed a hairpin vortex in the wake structure. In scheme U on fine mesh (e.g., see Re300D005A in Figure 8b), the hairpin vortex was distorted. In scheme S, however, the vortex was highly resolved (see Figure 8c,d). The wake structure was missing in scheme U because of the convection term by the first-order upwind finite-difference scheme, which cannot properly resolve the unsteady phenomena because of the large numerical dissipation. The same trend was observed at \(Re=400\) (see Figure 8e,f).