# A Mesh-Free Particle Method for Simulation of Mobile-Bed Behavior Induced by Dam Break
## Moving Particle Semi-Implicit Method
### Particle Interaction Model
These simultaneous linear symmetric matrix equations can be solved by the incomplete Cholesky conjugate gradient method. By substituting \(\phi\) by \(P^{n+1}\) in Equation (5), which is the gradient model, the pressure gradient can be obtained. Equation (12) is used as a source term for PPE. However, an exaggerated fluctuation of pressure took place. Reference [14] introduced the multi-source term to solve this problem. By adding a divergence of velocity in Equation (12), the pressure fluctuation was suppressed. For this treatment, [PERSON] et al. applied a relaxation parameter to distribute the weight of each terms. The source term of PPE used in this study is expressed below:  
\[\left\langle\
abla^{2}P_{i}^{n+1}\right\rangle_{i}=\gamma_{0}\frac{\rho_{0}}{ \Delta t^{2}}\frac{n^{0}-\langle n^{n}\rangle_{i}}{n^{0}}+(1-\gamma_{0})\frac {\rho_{0}}{\Delta t}\
abla\cdot\overset{\rightarrow}{u}^{*} \tag{13}\]  
where \(\overset{\rightarrow}{u}^{*}\) is the intermediate velocity, which is the temporal quantities obtained in the explicit stage, \(P^{n+1}\) is pressure at time \(n+1\), and \(\gamma_{0}\) is the relaxation parameter. The kernel function is a key factor of MPS, but it has an omni-spread function. It may occur under the distinct pressure head, and the buoyancy force can be underestimated at multi-phase flows. It can be solved by applying a self-buoyancy force to the center particle. Self-buoyancy is an adjustment of underestimated buoyant values. This term was applied to the velocity corrector term:  
\[\left\langle\
abla\cdot\phi\right\rangle_{i}=\frac{d}{n_{0}}\!\!\sum_{j\
eq i }\!\left[\frac{(\phi_{ij})\cdot(r_{ij})}{\big{(}|r_{ij}|\big{)}^{2}}w\big{(}| r_{ij}|\big{)}\right] \tag{5}\]0  
where \(w_{B}\) is the kernel function for the buoyancy-correction model, \(z\) is a vertical position of a particle, and \(C_{B}\) is a velocity-adjustment parameter. \(w_{B}\) is the same equation as Equation (3), but its effective range is set to 1.5, which was obtained from a numerical test. The reason for using \(C_{B}\) is due to the self-buoyancy term remedy being an underestimated buoyancy force. Therefore, the buoyancy correction term should be adjusted by a proper parameter. The value of \(C_{B}\) can be chosen by a numerical experiment. The buoyancy correction model was suggested by [PERSON] et al. [23], and it has been modified by the author using a physical formulation. In general, the surface tension is not a strong element to simulate the flow. However, it is not an ignored component for the interface of multi-phase flows. The fourth term on the left side of Equation (2) is the surface tension term. The particle method for tracking each particle is advantageous for a very large deformation of the free surface and interface. Yet, the weakness of the specialty lead particle method makes it difficult to indicate the interface layer as well. Therefore, the special treatment to define a curvature of the interface is introduced by [PERSON] et al. [21]. For the treatment, two additional particle number densities are employed. \[ \begin{array}{l}n_{i}^{st1}=\sum\limits_{j\
eq i}w^{st1}\big{(}\big{|}r_{ij }\big{|}\big{)}\\ n_{i}^{st2}=\sum\limits_{j\
eq i}w^{st2}\big{(}\big{|}r_{ij}\big{|}\big{)}\\ \text{where}\\ w^{st1}=\left\{\begin{array}{ll}1&\big{(}0<r_{ij}\leq r_{e}^{st}\big{)}\\ 0&\big{(}r_{ij}>r_{e}^{st}\big{)}\\ 1&\big{(}0<r_{ij}\leq r_{e}^{st}\text{and}r_{j}^{st1}>n_{i}^{st1}\big{)}\\ 0&\big{(}\text{otherwise}\big{)}\end{array}\right.\end{array}  \tag{15}\]  
where \(r_{e}^{st}\) is the effective range of the surface tension, and it is set to 1.5. With these newly added particle number densities, the curvature of the interface can be calculated by using Equation (16). \[\kappa=\frac{2}{r_{e}^{st}}\cos\left(0.5\pi\frac{n_{i}^{st2}}{n_{i}^{st1}}\right) \tag{16}\]  
After calculating the interface curvature, a normal vector of the interface needs to be estimated. [PERSON] et al. suggested calculating the normal vector of interface by using particle number density [21].