# Experimental and Numerical Investigation on Ballast Flight from Perspective of Individual Particles
## 2 Materials and Methods
### Computational Fluid Dynamics Simulation
#### 2.3.1 Model Formulation
The continuous equation is as follows:  
\[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{ \partial z}=0 \tag{1}\]  
The momentum equation is as follows:  
\[\frac{\partial u}{\partial t}+(u\frac{\partial u}{\partial x}+v\frac{\partial u }{\partial y}+w\frac{\partial u}{\partial z})=-\frac{1}{\rho}\frac{\partial P }{\partial x}+(v+v_{T})(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2 }u}{\partial y^{2}}+\frac{\partial^{2}u}{\partial z^{2}}) \tag{2}\]  
\[\frac{\partial v}{\partial t}+(u\frac{\partial v}{\partial x}+v\frac{\partial v }{\partial y}+w\frac{\partial v}{\partial z})=-\frac{1}{\rho}\frac{\partial P }{\partial y}+(v+v_{T})(\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^ {2}v}{\partial y^{2}}+\frac{\partial^{2}v}{\partial z^{2}}) \tag{3}\]  
\[\frac{\partial w}{\partial t}+(u\frac{\partial w}{\partial x}+v\frac{\partial w }{\partial y}+w\frac{\partial w}{\partial z})=-\frac{1}{\rho}\frac{\partial P }{\partial z}+(v+v_{T})(\frac{\partial^{2}w}{\partial x^{2}}+\frac{\partial^ {2}w}{\partial y^{2}}+\frac{\partial^{2}w}{\partial z^{2}}) \tag{4}\]  
The turbulent kinetic energy (\(k\)) equation is as follows:  
\[ \begin{array}{ll}&\frac{\partial k}{\partial t}+\frac{\partial(uk)}{ \partial x}+\frac{\partial(vk)}{\partial y}+\frac{\partial(wk)}{\partial z} \\ =&(G-\rho\varepsilon)+\frac{\partial}{\partial x}(\frac{v_{T}}{\partial_{k}} \frac{\partial k}{\partial x})+\frac{\partial}{\partial y}(\frac{v_{T}}{ \partial_{k}}\frac{\partial k}{\partial y})+\frac{\partial}{\partial z}(\frac {v_{T}}{\partial_{k}}\frac{\partial k}{\partial z})\end{array}  \tag{5}\]  
The turbulence dissipation (\(\varepsilon\)) equation is as follows:  
\[ \begin{array}{ll}&\frac{\partial c}{\partial t}+\frac{\partial(u \varepsilon)}{\partial x}+\frac{\partial(v\varepsilon)}{\partial y}+\frac{ \partial(w\varepsilon)}{\partial z}\\ =&\frac{\varepsilon}{k}(C_{1}G-C_{2}\rho\varepsilon)+\frac{\partial}{ \partial x}(\frac{v_{T}}{\partial_{\varepsilon}}\frac{\partial\varepsilon}{ \partial x})+\frac{\partial}{\partial y}(\frac{v_{T}}{\partial_{\varepsilon}} \frac{\partial\varepsilon}{\partial y})+\frac{\partial}{\partial z}(\frac{v_{ T}}{\partial_{\varepsilon}}\frac{\partial\varepsilon}{\partial z})\end{array}  \tag{6}\]  
where \(u\), \(v\) and \(w\) are the local velocity components of the flow field in each coordinate direction; \(\rho\), \(P\) and \(v\) is the density, the pressure and the viscosity coefficient of the airflow respectively; \(v_{T}\) is the effective viscosity coefficient; \(G\) is the production term of the turbulent kinetic energy (\(k\)); \(\sigma_{k},\sigma_{\varepsilon},C_{1}\) and \(C_{2}\) are the empirical constants. The second-order upwind style of the finite volume method (FVM) was used to discretize the governing equations. The algebraic equations with discretized initial conditions and boundary conditions were solved by the semi-implicit method for pressure linked equations (SIMPLE), based on which the numerical solution of airflow was obtained.