Siphon pipe parameter optimization of the toilet using CFD-DEM coupling method

<!--==Siphon pipe Parameter Optimization of the Toilet Using CFD-DEM Coupling Method ==

Zhanfu Li<sup>1,3</sup>,Xiaole Ge<sup>2</sup>,Youji Zhan<sup>1</sup>,Xin Tong<sup>1*</sup>,Zhihong Wu<sup>3</sup>, Feng Zhang<sup>3</sup>

(1.Mechanical Engineering training base, School of Mechanical and Automobile Engineering, Fujian University of Technology, Fuzhou, China,350118

2. College of Mechanical and Electrical Engineering, Huangshan University, Huangshan,China,245041

<span id='_GoBack'></span>
3. CSCEC Strait Construction and Development Co., Ltd, Fuzhou, China,350118

*Corresponding author: Email [mailto:addresses:Tongxin_hqu_fjut@163.com addresses:Tongxin_hqu_fjut@163.com], [mailto:19871985@fjut.edu.cn 19871985@fjut.edu.cn])
-->
==Abstract==

The computational fluid dynamics (CFD) and discrete element method (DEM) coupling method is used to simulate the flushing process of the toilet and this paper analyzes the influence of different structural parameters of siphon pipe on the flushing performance of the toilet. And through the adjusted-parameters toilet,   tlushing experments were carried out to verify the simulations. Meanwhile the orthogonal test of different structural parameters of siphon pipes were conducted to study the flushing Performance. The research results show that the CFD-DEM coupling method can be used to study the regularity of the toilet flushing performance. The toilet can get the better flushing performance when the tilt angle of the angle of inclination is 50°, the curvature width and length are 50 mm and 220 mm, the width and height of the secondary water seal are 100 mm and 25 mm and pipe diameter is 53 mm. The method in this paper can provide a new idea for the study and design of the flushing performance of the toilet.

'''Keyword''': CFD-DEM, flushing performance, structural parameters, coupling method, adjusted-parameters toilet

==1. Introduction==

Toilets are essential necessities in our daily life and water consumption of toilets is large among a variety of sanitary equipment [1]. The toilet flushing performance is an important index of quality assessment. And thus, constantly improving the flushing performance of the toilet has been the goal for researchers and workers in sanitary ware industry.

So far, there are two main research methods for the toilet flushing performance: simulation research and experimental research. In simulation research, CFD method was widely used by most scholars to complete the simulation of the toilet flushing process, which can only consider water flow rule [2,3]. In the experimental research, it is convenient for us to get a quantitative flushing effect, but it is not easy to get the complete flow field data [4-6]. In fact, the complex flushing process of toilet contains the interactions of flow-solids and solids-walls besides the movement process of water flow. All these factors have a certain impact on the flushing performance of the toilet. Therefore, the exploration of a new method considering the movement process of solids in the flushing process will be an important part in the study of the toilet flushing performance.

In order to incorporate the movement process of solids into the flushing process of the toilet, it is necessary for us to seek a computational model to simulate the comprehensive fluid and solid systems. And the CFD-DEM coupling method right provides a computational method for the complex model. The basic idea is that the computational fluid dynamics is used to calculate the flow of fluid, and DEM is used to calculate the force and movement of particles. The advantage of the CFD-DEM coupling method is that it can not only really simulate the flow rule of flow field, but also can take the particles’ shapes, collision, movement path and adhesion of particles into account. More importantly, it has been applied widely in many fields such as geology, sands flow, drainage machinery and so on [7-9]. In this paper, the CFD-DEM coupling method is used to study and simulate the flushing process of the toilet. Further, we verify the simulation rule by physical simulation experiments. At the same time, the orthogonal test is used to study different structural parameters of siphon pipes. All these research contents can provide references for the flushing performance analysis of the toilet and the design of siphon pipes.

==2. Theory==

===2.1 CFD theory===

The flushing process of the toilet is a typical process of fluid flow. The process of fluid flow must observe the continuity equation and the law of conservation of momentum, they can be expressed as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \rho }{\partial t}+\frac{\partial \rho u_i}{\partial x_i}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \rho u_i}{\partial t}+\frac{\partial \rho u_iu_j}{\partial x_j}=\rho f_i-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}(\mu \frac{\partial u_i}{\partial x_j})</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

The flushing process of the toilet contains the flow of free liquid surface. So far, the Volume of Fluid (VOF) method is widely used to solve the free liquid surface. VOF method is a computational method proposed by Hirt and Nichols for dealing with complex problems of free liquid surface [10]. The advantage of this method is that it can represent the complex change process of free surface by only a function. The mode of VOF based on the fluid volume method is widely used in dealing with stratified flow, free-surface flow and sloshing, and it also gets effective practice in the study for the flushing process of the toilet [2].

In the mode of VOF, different fluid components all observe the same momentum equation. In the calculation, the volume fraction of the components of different fluids in each computational element of flow field can all be stored, and the sum of the volume fractions of all phases in each control body is one. For ease of expression, assuming the volume fraction of the <math>j</math> phase in <math>i</math> element is <math>a_j</math>, three different cases are as follows:

:(1) When <math>a_j= 0 </math>, it means that the j phase in i element is empty;

:(2) When <math>a_j= 1 </math>, it means that the j phase in i element is completely full;

:(3) When <math>0 < a_j <1</math>, it means that there are j phase and other phases in i element.

The model of VOF can capture the interface information between different phases by the continuity equation of the volume fraction for one or more phases, such as take i phase as an example and it can be expressed as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial a_i}{\partial t}+\overset{\rightarrow}{v}\cdot \nabla a_i=</math><math>\frac{S_{a_i}}{{\rho }_i}+\frac{1}{{\rho }_i}\sum_{j=1}^n\left({\dot{m}}_{ji}-\right. </math><math>\left. {\dot{m}}_{ij}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where  <math>{\dot{m}}_{ji}</math> means the mass transport from <math>j</math> phase to i phase and  <math>{\dot{m}}_{ij}</math> means the mass transport from <math>i</math> phase to <math>j</math> phase.

Generally speaking, the value of the item on the right is 0. When the equation of the volume fraction is not based on the principal phase, the volume fraction of the principal phase can be calculated by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\sum_{i=1}^na_i=1</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

when there are only two phases in the whole system, the density of each element can be calculated by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\rho =a_2{\rho }_2+(1-a_2){\rho }_1</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}


By the calculation of a single momentum equation in all regions, the results are shared by all phases. The volume fraction of the phase expressed by momentum equation can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial }{\partial t}(\rho \overset{\rightarrow}{v})+</math><math>\nabla \cdot (\rho \overset{\rightarrow}{v}\overset{\rightarrow}{v})=</math><math>-\nabla p+\nabla \cdot [\mu (\nabla \overset{\rightarrow}{v}+</math><math>\nabla {\overset{\rightarrow}{v}}^T)]+\rho \overset{\rightarrow}{g}+</math><math>\overset{\rightarrow}{F}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

===2.2 DEM Theory===

DEM is a new computational method proposed by Cundall and Track [11], which can be used to solve the laws of motion and mechanical properties of complex discrete systems, and it is first applied effectively in geomechanics. DEM has been playing an important role in many fields such as granular engineering since it was proposed [12,13]. Its basic idea is to divide a discontinuity into a series of individuals with many rigid elements, and make all rigid elements satisfy the equations of motion. Then, the equations of motion of all rigid elements are calculated by time-step iteration method to get the global motion state of the discontinuity. The calculation principle of discrete element method can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\ddot{x}}_j(T)=\partial {\dot{x}}_j(T)/\partial T={\left(\sum F(T)\right)}_j/m_j</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\ddot{\theta }}_j(T)=\partial {\dot{\theta }}_j(T)/\partial T={\left(\sum M(T)\right)}_j/I_j</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

where <math>x_j</math> is the translational shift of <math>j</math> unit, <math>\sum F (T)</math> is the resultant force of j unit, <math>M_j</math> is the sum of <math>j</math> unit’s quality, <math>\theta_j</math> is the angular shift of <math>j</math> unit, <math>\sum M (T)</math> is the resultant moment of <math>j</math> unit.

A new shift of the unit can be obtained by the difference method or Euler method, and the new contact force and torque can be generated. Then it returns to formulas (7) and formulas (8) to update the calculation to get the motion information of all units. In the flushing process of the toilet, the solid particles in the flow field will be affected by many effect for the gravity of the particles themselves, the contact-impact force between the particles and the inner wall of siphon pipes, the fluid drag force, the fluid buoyancy, the pressure-gradient force, the false mass force and the lift force. The related results show that the drag force is the main interaction force between fluid and particles, and the effect of the false mass force and the lift force on particles can be ignored [14].

The motion of particles in the flow field can be expressed in two forms: translational motion and rotational motion. For any particle <math>i</math>, the equation of translational motion can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>m_i\frac{d\mbox{v}_i}{dt}=m_i\mbox{g}+\mbox{F}_b+\mbox{F}_c+\mbox{F}_d+\mbox{F}_p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

where <math>m_i</math> is the mass of particle <math>i</math>; <math>v_i</math> is the speed of particle <math>I, \, F_b</math>  is the buoyancy of particle <math>I, \, F_c</math> is the contact-impact resultant force of particle <math>I, \, F_d</math> is the drag force of particle <math>I, \, F_P</math> is the pressure-gradient force of particle <math>i</math>.

The buoyancy force of particles in the fluid <math>F_b</math> can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{F}_b=\rho gV_p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}


The contact-impact force of particles <math>F_c</math> can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{F}_c=\mbox{F}_{cni}+\mbox{F}_{cti}+\mbox{F}_{dni}+</math><math>\mbox{F}_{dti}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

where <math>F_{cni}</math> is the normal contact force of particles, <math>F_{cti}</math> is the tangential contact force of particles, <math>F_{dni}</math> is the normal damping force of particles, <math>F_{dti}</math> is the tangential damping force of particles.

The magnitude of contact-impact force depends on the contact model when calculating, the contact-impact force is from between particles and between particles and the inner walls of siphon pipes. Hertz-Mindlin’s no slip model is used to describe this contact-impact behavior. The basic idea for Hertz-Mindlin’s no slip model is to use Hertz contact theory to calculate the normal force and Mindlin-Deresiewicz contact model to calculate the tangential contact force[15]. In this paper, spherical particles are used to study the flushing process of the toilet, the free resistance model is used to study the drag force model, and this model can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{F}_d=0.5C_D\rho A\vert v\vert v</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where <math>C_D</math> is the drag coefficient, <math>\rho</math> is the fluid density, <math display="inline"> A </math> is the projected area of particles, <math display="inline">  |v|  </math>  is the variable size of speed, <math display="inline"> v  </math> is the relative speed.

The drag coefficient <math display="inline">C_D  </math> can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>C_D=\lbrace \begin{array}{c}
\mbox{                 24/}Re\mbox{                          }Re\leq 0.5\\
24(1.0+0.25Re^{0.687})/Re\mbox{    }0.5<Re\leq 1000\\
\mbox{                 }0.44\mbox{                            }Re>1000
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

here <math display="inline"> Re  </math> is the Reynolds number. The pressure-gradient force can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{F}_p=-V_i\nabla p_s</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">  V_i </math> is the volume of particle <math display="inline">  i </math>,  <math>\nabla p_s</math> is the static pressure gradient.

For any particle i, the equation of rotational motion in the flow field can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>I_i\frac{d{\mbox{ω}}_i}{dt}=\sum_{i=1}^n\left(\mbox{T}_i\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline"> I_i  </math> is the moment of inertia of particles, <math display="inline"> \omega_i  </math> isthe angular speed of particles, <math display="inline">T_i  </math> is the torque of particles.

===2.2 CFD-DEM coupled method===

CFD-DEM coupling method is a numerical analysis method for multiphase flow problems of particles and fluids. It has been widely used in many fields. Its basic idea is to model the fluid as a continuous medium and calculate it in the CFD method，and DEM and Newton's second law are combined to calculate the motion of particles which are considered to be a discrete phase.

The flushing process of the toilet consists of two parts: one is the flow of water flow, the other is the movement of particles under the effect of water flow. In the process of the flow of water flow, water flow is considered to be a continuous phase, it can be calculated by solving the continuity equation and N-S equation. In the process of particle motion, particles are considered to be a discrete phase, the force and motion of particles can be calculated by Newton's second law. The data for the process of the flow of water flow and particle motion can be transmitted through the interface card.The coupling process for the flushing process of the toilet in FLUENT and EDEM is shown in [[#img-1|Figure 1]]  [16-18].

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Li_210280261-image19.jpeg|600px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 1'''. The calculation process of CFD-DEM coupled method
|}


According to [[#img-1|Figure 1]], in the calculation, FLUENT first calculates the flow field of an interpolation time period, and transmits the flow field information of the interpolation point to EDEM at an interpolation time interval. EDEM reads the flow field information from FLUENT, then adds various field forces to the motion equation of particles by the interface plug-in (API) and start to calculate the forces and motion of particles in the corresponding interpolation time interval. After EDEM completes the calculation in the interpolation time period, FLUENT starts to update the flow field in the next interpolation time period in circulation. The interpolation time period is the time interval of data transmission, and the value is taken the smaller and the calculation error is the smaller. The minimum interpolation time interval is a FLUENT calculation time step, that is the time step is the smaller, the calculation accuracy is the higher.

The flow field information transmitted from FLUENT to EDEM includes velocity field, volume fraction field (or density field) and pressure-gradient field. The main reason why the volume fraction field or density field need to be transmitted is that the calculation domain contains two kinds of mediums: water flow and air in the flushing process of the toilet, the field forces (buoyancy force and drag force) of particles in different mediums are different. Therefore, it is necessary to distinguish the properties of mediums in various locations in fluid domain. The interface plug-in between FLUENT and EDEM is compiled by Visual Studio development tool through API. The contents of the interface plug-in include reading the flow field information, distinguishing the properties of mediums in calculation domain, calculating field forces and setting the interpolation time interval for a transient calculation. In this paper, the effect of water flow on particles is mainly considered, and it is a one-way coupling method.

==3. Simulation and experimental conditions==

===3.1 Simulation parameter settings===

====3.1.1 CFD parameter settings====

In this paper, the flushing process of a siphon jet toilet with different structural parameters of siphon pipes is studied by CFD-DEM coupling method. According to the structural characteristics of the toilet model and considering the calculation time and accuracy, the hybrid grid is used to mesh, that is, the water tank part with regular structure is divided into structured grids to reduce the number of grids, and the unstructured grids are used to mesh in other areas to reduce the difficulty of meshing. The number of grids is about 2.3 million ([[#img-2|Figure 2]]a).

In the CFD calculation of the flushing process, the pressure base solver is used for a transient solution. The VOF model is used to calculate the multiphase flow, and the air is setted as the first phase, water is setted as the second phase. The Realizable <math>k</math>-epsilon model is chosen for the turbulence model and the standard wall function is applied to deal with the near-wall. Pressure-velocity coupling method adopts SIMPLE algorithm, the momentum equation is discretized by second-order upwind, and turbulent kinetic energy, turbulent kinetic energy dissipation rate and transient equation are discretized by first-order upwind. According to the actual situation of the flushing process, a pressure inlet boundary condition and two pressure outlet boundary conditions are setted up, as shown in [[#img-2|Figure 2]]. In [[#img-2|Figure 2]]b, the red part is water, the blue part is air, and the water in the water tank is 5.5 L. The computation time step of FLUENT is 0.0002 s, the computation time is 6.6 s, and the maximum iteration step is 20 steps.

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Li_210280261-image20.png|264px]]
| style="padding:10px;"|[[Image:Draft_Li_210280261-image21.png|center|318px]]
|-
|  style="text-align: center;font-size: 75%;"|(a) 
|  style="text-align: center;font-size: 75%;"|(b) 
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding:10px;"| '''Figure 2'''. (a) The assembled grid model. (b) The setting of boundary conditions 
|}

===3.1.2 DEM parameter settings===

The solid medium in the flushing process is polypropylene spherical particles with a radius of 8.75 mm, a density of 890 kg/m<math>^3</math>, a flushing quantity of 130, and the wall material of the toilet is ceramic. The formation of particles is set up by the secondary development of the granular factory, that is all water tested balls are generated at the water seal position of the toilet at the beginning of the calculation, and this ensures that the initial state of the water tested balls is the same in each test in order to avoid the influence of the randomness of particle formation on the flushing results. Considering that the contact force between particles and between particles and the inner walls of the toilet, it is necessary to determine the contact parameters between particles and between particles and the walls. In this paper, a high-speed camera is used to measure the contact parameters. Detailed parameter settings are shown in [[#tab-1|Table 1]]. The computation time step of''' '''EDEM is setted as <math>2\times 10^{-6}</math> s.

<div class="center" style="font-size: 75%;">
'''Table 1'''. The setting of EDEM parameters</div>

<div id='tab-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! Material properties !! Poisson ratio !! Density（kg/m<math>^3</math>） !! Shear modulus（Pa）
|-
|  style="text-align: center;"|Particles
|  style="text-align: center;"|0.41
|  style="text-align: center;"|890
|  style="text-align: center;"|3.2e+08
|-
|  style="text-align: center;"|Pipe wall
|  style="text-align: center;"|0.22
|  style="text-align: center;"|2300
|  style="text-align: center;"|9.06e+10Pa
|-
|  style="text-align: center;"|Collision attributes
|  style="text-align: center;"|Restitution coefficie
|  style="text-align: center;"|Static friction coefficient
|  style="text-align: center;"|Rolling friction coefficient
|-
|  style="text-align: center;"|particles-particles
|  style="text-align: center;"|0.6282
|  style="text-align: center;"|0.12
|  style="text-align: center;"|0.0486
|-
|  style="text-align: center;"|particles-walls 
|  style="text-align: center;"|0.7823
|  style="text-align: center;"|0.2
|  style="text-align: center;"|0.0306
|}

===3.1.3 Structure parameter setting of siphon pipes===

In this paper, the influence of structural parameters of siphon pipes on the flushing performance of the toilet is studied. In order to better analyze, the structures of siphon pipes is divided, which are shown in [[#img-3|Figure 3]].

<div id='img-3'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Li_210280261-image22.png|246px]] 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figure 3'''. The division of structural parameters of siphon pipes 
|}


In [[#img-3|Figure 3]], the structure size of siphon pipes is characterized by six parameters: <math>R</math>, <math>L_1</math>, <math>H_1</math>, <math>L_2</math>, <math>H_2</math> and <math>D</math>. <math>R</math> is the inclination angle of A-B section of siphon pipes, <math>H_0</math> is the water seal height of of siphon pipes, <math>L_1</math> is the length of radian change, <math>H_1</math> is the height of radian change, <math>L_2</math> is the length of E-F section of secondary water seal, <math>H_2</math> is the height of secondary water seal, and <math>D</math> is the diameter of siphon pipes. <math>R</math>, <math>L_1</math>, <math>H_1</math>, <math>L_2</math>, <math>H_2</math> and D are named as inclination angle, curvature width, curvature length, secondary water seal width, secondary water seal height and pipe diameter, respectively. All these parameters are both connected and independent, when one of them is changed, the other parameters remain unchanged. In order to explore the effect of six parameters on the flushing performance, six tests are carried out, and the detailed parameters are shown in [[#tab-2|Table 2]].

<div class="center" style="font-size: 75%;">
'''Table 2'''. The structure parameters of siphon pipes</div>

<div id='tab-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! NO. !! ''R''（°） !! <math>L_1</math> （mm） !! <math>H_1</math> （mm） !! <math>L_2</math> （mm） !! <math>H_2</math> （mm） !! <math>D</math> （mm）
|-
|  style="text-align: center;vertical-align: top;"|1
|  style="text-align: center;vertical-align: top;"|36,40,44,48,52,56
|  style="text-align: center;vertical-align: top;"|55
|  style="text-align: center;vertical-align: top;"|210
|  style="text-align: center;vertical-align: top;"|90
|  style="text-align: center;vertical-align: top;"|30
|  style="text-align: center;vertical-align: top;"|50
|-
|  style="text-align: center;vertical-align: top;"|2
|  style="text-align: center;vertical-align: top;"|48
|  style="text-align: center;vertical-align: top;"|45,50,55,60,65,70
|  style="text-align: center;vertical-align: top;"|210
|  style="text-align: center;vertical-align: top;"|90
|  style="text-align: center;vertical-align: top;"|30
|  style="text-align: center;vertical-align: top;"|50
|-
|  style="text-align: center;vertical-align: top;"|3
|  style="text-align: center;vertical-align: top;"|48
|  style="text-align: center;vertical-align: top;"|65
|  style="text-align: center;vertical-align: top;"|170,190,210,230,250
|  style="text-align: center;vertical-align: top;"|90
|  style="text-align: center;vertical-align: top;"|30
|  style="text-align: center;vertical-align: top;"|50
|-
|  style="text-align: center;vertical-align: top;"|4
|  style="text-align: center;vertical-align: top;"|48
|  style="text-align: center;vertical-align: top;"|65
|  style="text-align: center;vertical-align: top;"|210
|  style="text-align: center;vertical-align: top;"|60,70,80,90,100,110
|  style="text-align: center;vertical-align: top;"|30
|  style="text-align: center;vertical-align: top;"|50
|-
|  style="text-align: center;vertical-align: top;"|5
|  style="text-align: center;vertical-align: top;"|48
|  style="text-align: center;vertical-align: top;"|65
|  style="text-align: center;vertical-align: top;"|210
|  style="text-align: center;vertical-align: top;"|90
|  style="text-align: center;vertical-align: top;"|10,20,30,40,50,60
|  style="text-align: center;vertical-align: top;"|50
|-
|  style="text-align: center;vertical-align: top;"|6
|  colstyle="text-align: center;vertical-align: top;"|48
|  style="text-align: center;vertical-align: top;"|65
|  style="text-align: center;vertical-align: top;"|210
|  style="text-align: center;vertical-align: top;"|90
|  style="text-align: center;vertical-align: top;"|30
|  style="text-align: center;vertical-align: top;"|44,47,50,53,56,59
|}

===3.2 Experimental platform===

In order to study the influence of structural parameters of siphon pipes on the flushing performance, a set of experimental platform with adjustable structural parameters of siphon pipes need to be builded. The simulation model used in the numerical calculation is based on a siphon jet toilet. Therefore, the building of the physical experimental platform is also based on it, and the improved toilet model is shown in [[#img-4|Figure 4]]. In this experimental platform, the shape of siphon pipes can be adjusted arbitrarily within a certain range, and the water volume in the water tank can be controlled automatically through a flowmeter. The structural parameters of the siphon pipes are consistent with the parameters in the simulation experiment, and the water volume in the water tank is 5.5 L. It should be noted that the influence of siphon pipe diameter on the flushing performance has not been studied in the experiment, this is because the siphon pipe diameter is related to the basin structure of the toilet, when the basin structure is determined, the pipe diameter is also determined, and the siphon pipe diameter is 50 mm in this experiment. In addition, because of the size limitation of the toilet model, there are only three values of the curvature length parameters are studied. In the experiment, the material properties and quantity of solid granular medium are the same as those in the simulation, and each experiment is carried out six times, then the average value is taken as the quantity of flushing balls.

<div id='img-4'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Li_210280261-image24.png|228px]]
| style="padding:10px;"| [[Image:Draft_Li_210280261-image25.jpeg|center|264px]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding-bottom:10px;"| '''Figure 4'''. Experimental platform of toilet 
|}

===3.3 Definition of flushing capacity===

In order to analyze quantitatively the flushing performance of the toilet and consider synthetically the relationship between the quantity of flushing balls and the amount of water in the flushing process, the flushing capacity is defined as the evaluation index of the flushing performance of the toilet, and the flushing capacity can be expressed by the following formula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>P=\frac{V_B}{V_W}</math> 
|}
|  style="text-align: right;vertical-align: top;width: 5px;text-align: right;white-space: nowrap;"|(16)
|}

where <math>V_B</math> is the volume of all tested water balls that are washed away, <math>V_W</math> is the volume of water used in the flushing process. According to the formula of the flushing capacity, the flushing capacity can be understood as the volume of tested water balls that are washed away by a unit volume of water.

==4 The discussion of flushing performance==

===4.1 Verification of simulation law===

The relationship between the flushing capacity of the toilet and structural parameters of siphon pipes can be obtained by studying the flushing process of different siphon pipe structural parameters using CFD-DEM coupling method. In order to verify the reliability of the simulation law, the experiments under the same structural parameters are carried out, and the simulation law is compared with the experimental law, the results are shown from [[#img-5|Figure 5]]  to [[#img-9|Figure 9]].

<div id='img-5'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Li_210280261-image27-c.png|294px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 5'''. Comparison of inclination angle
|}

<div id='img-6'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Li_210280261-image28-c.png|center|288px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 6'''. Comparison of curvature width
|}

<div id='img-7'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Li_210280261-image29-c.png|288px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 7'''. Comparison of curvature length
|}

<div id='img-8'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Li_210280261-image30-c.png|center|300px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 8'''. Comparison of secondary water seal width
|}

<div id='img-9'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Li_210280261-image31-c.png|306px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 9'''. Comparison of secondary water seal height 
|}


As are shown from [[#img-5|Figure 5]]  to [[#img-9|Figure 9]], the experimental results and the simulation results have the same change rule, and only the absolute values are different. And the difference are mainly reflected in the value of the flushing capacity in the simulation, which is larger than that in the experiment. On the one hand, it is caused by the inconsistency of physical experiment model and simulation calculation model, on the other hand, it is also affected by the accuracy of coupling algorithm in simulation calculation. Although the simulation law and the experimental law are different in the absolute value, the factors that cause the difference all have the same effect on the absolute value of flushing capacity under all parameters, and as a whole, it can reflect the range of structural parameters with better flushing performance. Therefore, it is practicable to study the regularity of the flushing performance of structural parameters of different siphon pipes. The numerical simulation method is not limited by the fixed size of model, it can study the flushing process in a wider range of structural parameters. And it is also convenient to obtain data which is not easy to get in physical experiments, and it is helpful to further study the flushing mechanism and reduce the development cycle and cost.

Because of the size limitation of experimental platform, there is no study about the effect of diameter change of siphon pipes on the flushing performance. Therefore, we can use CFD-DEM method to simulate the flushing process under different diameters. The simulation results are shown in [[#img-10|Figure 10]].

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Li_210280261-image32-c.png|336px]] 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figure 10'''. The influence of siphon pipe diameter on the flushing performance 
|}


According to [[#img-10|Figure 10]]，the flushing capacity increases gradually and the flushing performance of the toilet increases gradually with the pipe diameter increases from 44 mm to 53 mm. When the pipe diameter increases from 53 mm to 56 mm, the flushing performance of the toilet remains unchanged, it is mainly because all the tested water balls are washed out under the diameters of these two pipes. When the pipe diameter is increased to 59 mm, the flushing capacity begins to decline. It is mainly because the pipe diameter is the larger, the siphon strength is the lower and the flushing force is the smaller under the same water consumption condition, which also shows that the bigger is not the better under the same flushing water consumption condition, and the reasonable pipe diameter is helpful to get better flushing and water-saving comprehensive properties. When the diameters of the pipe are 53mm and 56mm, the flushing performance of the toilet is better. In addition, the diameter of the siphon pipe is the smaller, the velocity of water flow is the higher in the pipe, and the siphon negative pressure is the lower, but the flushing capacity is lower, which also shows that the flushing performance can not be judged directly by the velocity and the negative pressure in the pipe under the condition of the pipe diameter changes, and it is necessary to combine the punching balls test.

==4.2 The study of orthogonal experiment==

'''   '''In order to explore the influence of the structural parameters of each siphon pipe on the flushing performance, and to find a better combination of the structural parameters of siphon pipes, an orthogonal test was carried out in this paper. According to the experimental part, the pipe diameter of the toilet is related to the structure of the basin wall. When the diameter of siphon pipes changes, the model of the basin wall also changes. If the orthogonal experiment is carried out, many different models of the toilet need to be made, which greatly increases the development cycle and cost. According to the above research, CFD-DEM coupling method can better reflect the change rule of the flushing capacity of the toilet, and we can use the simulation method to carry out the numerical analysis when it is not easy to carry out an experimental research. Therefore, we can use the simulation method to carry out the orthogonal test.

There are six parameters to be considered in the structure of siphon pipes, each parameter is taken at three levels, and the L<sub>18</sub>（3<sup>7</sup>）orthogonal table is selected for the test, a total of 18 experiments are required. The parameter settings and results of orthogonal test are shown in Table 3.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Table 3 The parameter settings of orthogonal test'''</div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  rowspan='3' style="border-top: 2pt solid black;border-bottom: 1pt solid black;text-align: center;"|Test number
|  colspan='6'  style="border-top: 2pt solid black;border-bottom: 1pt solid black;text-align: center;"|Factors
|  rowspan='3' style="border-top: 2pt solid black;text-align: center;vertical-align: bottom;"|Flushing capacity
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|A
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|B
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|C
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|D
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|E
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|F
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''R''（°）
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''L<sub>1</sub>''（mm）
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''H<sub>1</sub>''（mm）
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''L<sub>2</sub>''（mm）
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''H<sub>2</sub>''（mm）
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|''D''（mm）
|-
|  style="border-top: 1pt solid black;text-align: center;"|1
|  style="border-top: 1pt solid black;text-align: center;"|48
|  style="border-top: 1pt solid black;text-align: center;"|45
|  style="border-top: 1pt solid black;text-align: center;"|200
|  style="border-top: 1pt solid black;text-align: center;"|90
|  style="border-top: 1pt solid black;text-align: center;"|20
|  style="border-top: 1pt solid black;text-align: center;"|47
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|0.0561
|-
|  style="text-align: center;"|2
|  style="text-align: center;"|48
|  style="text-align: center;"|50
|  style="text-align: center;"|210
|  style="text-align: center;"|100
|  style="text-align: center;"|25
|  style="text-align: center;"|50
|  style="text-align: center;vertical-align: top;"|0.0648
|-
|  style="text-align: center;"|3
|  style="text-align: center;"|48
|  style="text-align: center;"|65
|  style="text-align: center;"|220
|  style="text-align: center;"|110
|  style="text-align: center;"|30
|  style="text-align: center;"|53
|  style="text-align: center;vertical-align: top;"|0.0628
|-
|  style="text-align: center;"|4
|  style="text-align: center;"|50
|  style="text-align: center;"|45
|  style="text-align: center;"|200
|  style="text-align: center;"|100
|  style="text-align: center;"|25
|  style="text-align: center;"|53
|  style="text-align: center;vertical-align: top;"|0.0663
|-
|  style="text-align: center;"|5
|  style="text-align: center;"|50
|  style="text-align: center;"|50
|  style="text-align: center;"|210
|  style="text-align: center;"|110
|  style="text-align: center;"|30
|  style="text-align: center;"|47
|  style="text-align: center;vertical-align: top;"|0.0434
|-
|  style="text-align: center;"|6
|  style="text-align: center;"|50
|  style="text-align: center;"|65
|  style="text-align: center;"|220
|  style="text-align: center;"|90
|  style="text-align: center;"|20
|  style="text-align: center;"|50
|  style="text-align: center;vertical-align: top;"|0.0526
|-
|  style="text-align: center;"|7
|  style="text-align: center;"|52
|  style="text-align: center;"|45
|  style="text-align: center;"|210
|  style="text-align: center;"|90
|  style="text-align: center;"|30
|  style="text-align: center;"|50
|  style="text-align: center;vertical-align: top;"|0.0413
|-
|  style="text-align: center;"|8
|  style="text-align: center;"|52
|  style="text-align: center;"|50
|  style="text-align: center;"|220
|  style="text-align: center;"|100
|  style="text-align: center;"|20
|  style="text-align: center;"|53
|  style="text-align: center;vertical-align: top;"|0.0663
|-
|  style="text-align: center;"|9
|  style="text-align: center;"|52
|  style="text-align: center;"|65
|  style="text-align: center;"|200
|  style="text-align: center;"|110
|  style="text-align: center;"|25
|  style="text-align: center;"|47
|  style="text-align: center;vertical-align: top;"|0.0439
|-
|  style="text-align: center;"|10
|  style="text-align: center;"|48
|  style="text-align: center;"|45
|  style="text-align: center;"|220
|  style="text-align: center;"|110
|  style="text-align: center;"|25
|  style="text-align: center;"|50
|  style="text-align: center;vertical-align: top;"|0.0622
|-
|  style="text-align: center;"|11
|  style="text-align: center;"|48
|  style="text-align: center;"|50
|  style="text-align: center;"|200
|  style="text-align: center;"|90
|  style="text-align: center;"|30
|  style="text-align: center;"|53
|  style="text-align: center;vertical-align: top;"|0.0510
|-
|  style="text-align: center;"|12
|  style="text-align: center;"|48
|  style="text-align: center;"|65
|  style="text-align: center;"|210
|  style="text-align: center;"|100
|  style="text-align: center;"|20
|  style="text-align: center;"|47
|  style="text-align: center;vertical-align: top;"|0.0434
|-
|  style="text-align: center;"|13
|  style="text-align: center;"|50
|  style="text-align: center;"|45
|  style="text-align: center;"|210
|  style="text-align: center;"|110
|  style="text-align: center;"|20
|  style="text-align: center;"|53
|  style="text-align: center;vertical-align: top;"|0.0653
|-
|  style="text-align: center;"|14
|  style="text-align: center;"|50
|  style="text-align: center;"|50
|  style="text-align: center;"|220
|  style="text-align: center;"|90
|  style="text-align: center;"|25
|  style="text-align: center;"|47
|  style="text-align: center;vertical-align: top;"|0.0582
|-
|  style="text-align: center;"|15
|  style="text-align: center;"|50
|  style="text-align: center;"|65
|  style="text-align: center;"|200
|  style="text-align: center;"|100
|  style="text-align: center;"|30
|  style="text-align: center;"|50
|  style="text-align: center;vertical-align: top;"|0.0633
|-
|  style="text-align: center;"|16
|  style="text-align: center;"|52
|  style="text-align: center;"|45
|  style="text-align: center;"|220
|  style="text-align: center;"|100
|  style="text-align: center;"|30
|  style="text-align: center;"|47
|  style="text-align: center;vertical-align: top;"|0.0515
|-
|  style="text-align: center;"|17
|  style="text-align: center;"|52
|  style="text-align: center;"|50
|  style="text-align: center;"|200
|  style="text-align: center;"|110
|  style="text-align: center;"|20
|  style="text-align: center;"|50
|  style="text-align: center;vertical-align: top;"|0.0663
|-
|  style="text-align: center;"|18
|  style="text-align: center;"|52
|  style="text-align: center;"|65
|  style="text-align: center;"|210
|  style="text-align: center;"|90
|  style="text-align: center;"|25
|  style="text-align: center;"|53
|  style="text-align: center;vertical-align: top;"|0.0663
|-
|  style="text-align: center;"|Average value 1
|  style="text-align: center;"|0.0567
|  style="text-align: center;"|0.0571
|  style="text-align: center;"|0.0578
|  style="text-align: center;"|0.0543
|  style="text-align: center;"|0.0583
|  style="text-align: center;"|0.0494
|  style="text-align: center;vertical-align: top;"|
|-
|  style="text-align: center;"|Average value 2
|  style="text-align: center;"|0.0582
|  style="text-align: center;"|0.0583
|  style="text-align: center;"|0.0541
|  style="text-align: center;"|0.0593
|  style="text-align: center;"|0.0603
|  style="text-align: center;"|0.0584
|  style="text-align: center;vertical-align: top;"|
|-
|  style="text-align: center;"|Average value 3
|  style="text-align: center;"|0.0560
|  style="text-align: center;"|0.0554
|  style="text-align: center;"|0.0589
|  style="text-align: center;"|0.0573
|  style="text-align: center;"|0.0522
|  style="text-align: center;"|0.0630
|  style="text-align: center;vertical-align: top;"|
|-
|  style="border-bottom: 2pt solid black;text-align: center;"|Range
|  style="border-bottom: 2pt solid black;text-align: center;"|0.0022
|  style="border-bottom: 2pt solid black;text-align: center;"|0.0030
|  style="border-bottom: 2pt solid black;text-align: center;"|0.0048
|  style="border-bottom: 2pt solid black;text-align: center;"|0.0050
|  style="border-bottom: 2pt solid black;text-align: center;"|0.0081
|  style="border-bottom: 2pt solid black;text-align: center;"|0.0136
|  style="border-bottom: 2pt solid black;text-align: center;vertical-align: top;"|
|}


The range reflects the influence of each structural parameter change of siphon pipes on the flushing capacity, and the bigger range indicates that the structural parameters will cause a bigger change on the flushing capacity in the experiment, so the maximum row of range indicates that the level of its structural parameters is the most important factor affecting the flushing performance of the toilet. According to Table 3, the size relation of range is F>E>D>C>B>A, therefore, the order of the influence of the structural parameters of siphon pipes on the flushing capacity is from primary to secondary is pipe diameter, secondary water seal height, secondary water seal width, curvature length, curvature width and inclination angle. The optimized combination parameter of orthogonal test is A<sub>2</sub>B<sub>2</sub>C<sub>3</sub>D<sub>2</sub>E<sub>2</sub>F<sub>3, </sub>but the parameter combination don’t contain that one in the orthogonal test table. Therefore, it needs to be verified by experiments. The flushing capacity of the optimized combination scheme obtained by simulation method is 0.0663. Compared with the results of the orthogonal test, the flushing capacity of the fourth, eighth, seventeenth and eighteenth groups are all 0.0663, which is the same as that one in the optimized scheme. That is, all the solid particles are flushed out. Therefore, it is not yet able to determine whether the flushing capacity under the optimized parameters of the scheme is the optimal. In order to further determine the flushing performance between the optimized scheme and the fourth, eighth, seventeenth and eighteenth groups, the quantity of flushing balls is increased from 130 to 150, then the above five combination parameters are re-tested, and the results of the experiments are shown in Figure 11.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Li_210280261-image33-c.png|336px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 11 The flushing capacity under different groups</div>

'''  '''According to Figure 11, the flushing capacity is the biggest under the optimized scheme, this indicates that the best flushing capacity can be obtained when the inclination angle is 50°, the curvature width is 50 mm, the curvature length is 220 mm, the secondary water seal width is 100 mm, the secondary water seal height is 25 mm and the pipe diameter is 53 mm.

=5 Conclusions=

According to the above analysis, the following conclusions can be obtained:

1. CFD-DEM coupling method can be used to reflect the overall changing tendency of the flushing performance of the toilet, but the absolute value is different. Therefore, CFD-DEM coupling method can be used to study the regularity of the flushing performance of the toilet. When the experimental method is not easy, the simulation method can be used to analyze qualitatively the flushing performance of the toilet.

2. Better flushing performance of the toilet can be obtained when the inclination angle is 50°, the curvature width is 50 mm, the curvature length is 220 mm, the secondary water seal width is 100 mm, the secondary water seal height is 25 mm and the pipe diameter is 53 mm.

3. It is worth to further study about how to use CFD-DEM coupling method for quantitative analysis of the toilet flushing performance.

==Acknowledgement ==

The authors gratefully acknowledged the support from National Natural Science Foundation of China (51775113). Fujian university-industry cooperation project (2017H6002).Research start-up funds of FuJian University of Technology (GY-Z160139).Fuzhou Science and technology project (2020-GX-16).No part of this paper has published or submitted elsewhere. The authors declared that they have no conflicts of interest to this work. All authors have seen the manuscript and approved to submit to your journal.

==References==

1. Z. Shi-Yi, L. Zi-Jian, and P. Zhi-Wei, "Optimized Design of Toilet Siphon Pipeline Based on Fluent and Its Validated Test," Journal of System Simulation '''20''', 4412 (2008).

2. Y. Wang, G. Xiu, and H. Tan, "Cad and Cae Analysis for Siphon Jet Toilet," in ''International Conference on Optics in Precision Engineering and Nanotechnology, ICOPEN 2011, March 23, 2011 - March 25, 2011'', Ed. by (Elsevier, Singapore, Singapore, 2011), p. 472.

3. J. Mahecha, O. López, and O. Ardila, "Computational Study of the Dynamics of the Flow in a Gravity-Driven Toilet," in ''ASME 2012 Fluids Engineering Division Summer Meeting Collocated with the ASME 2012 Heat Transfer Summer Conference and the ASME 2012 International Conference on Nanochannels, Microchannels, and Minichannels'', Ed. by 2012), p. 1319.

4. I.Y. An, Y.L. Lee, and J.-H. Kim, "A Study of the Characteristics of a Super Water-Saving Toilet with Flexible Trapway by Measuring Accumulated Flow Rate," Journal of Mechanical Science and Technology '''28''', 3067 (2014).

5. I.Y. An, Y.L. Lee, E.D. Jung, and W.S. Cho, "A Study on Performance Optimization of a Toilet by Measuring Accumulated Flow Rate of a Trapway," Journal of Mechanical Science and Technology '''28''', 1319 (2014).

6. W.S. Cheng, R.T. Lee, C.H. Liu, and C.W. Hsia, "A Study on Evacuation Performance of Siphon‐Type Water Closets," AIP Conference Proceedings '''1225''', 75 (2010).

7. J. Chen, Y. Wang, X. Li, R. He, S. Han, and Y. Chen, "Reprint of "Erosion Prediction of Liquid-Particle Two-Phase Flow in Pipeline Elbows Via Cfd-Dem Coupling Method"," Powder Technology '''282''', 25 (2015).

8. K.M. Park, H.S. Yoon, and M.I. Kim, "Cfd-Dem Based Numerical Simulation of Liquid-Gas-Particle Mixture Flow in Dam Break," Communications in Nonlinear Science and Numerical Simulation '''59''', 105 (2018).

9. S. Akhshik, M. Behzad, and M. Rajabi, "Cfd-Dem Simulation of the Hole Cleaning Process in a Deviated Well Drilling: The Effects of Particle Shape," Particuology '''25''', 72 (2016).

10. C.W. Hirt, and B.D. Nichols, "Volume of Fluid (Vof) Method for the Dynamics of Free Boundaries," Journal of computational physics '''39''', 201 (1981).

11. P.A. Cundall, and O.D.L. Strack, "Discrete Numerical Model for Granular Assemblies," Geotechnique '''29''', 47 (1979).

12. J.P. De Bono, and G.R. Mcdowell, "An Insight into the Yielding and Normal Compression of Sand with Irregularly-Shaped Particles Using Dem," Powder Technology '''271''', 270 (2015).

13. Z. Li, X. Tong, and Y. Qiu, "The Research of Particle Sieving under a Creative Mode of Vibration," Journal of Vibroengineering '''19''', 4172 (2017).

14. M.A. Habib, H.M. Badr, R. Ben-Mansour, and S.a.M. Said, "Numerical Calculations of Erosion in an Abrupt Pipe Contraction of Different Contraction Ratios," International Journal for Numerical Methods in Fluids '''46''', 19 (2004).

15. H.D. Mindlin R D, "Elastic Spheres in Contact under Varying Oblique Forces," Journal of Applied Mechanics '''20''', 327 (1953).

16. K.W. Chu, B. Wang, A.B. Yu, and A. Vince, "Cfd-Dem Modelling of Multiphase Flow in Dense Medium Cyclones," Powder Technology '''193''', 235 (2009).

17. K.W. Chu, S.B. Kuang, A.B. Yu, and A. Vince, "Particle Scale Modelling of the Multiphase Flow in a Dense Medium Cyclone: Effect of Fluctuation of Solids Flowrate," in Ed. by (Elsevier Ltd, 2012), p. 34.

18. F. Chaumeil, and M. Crapper, "Using the Dem-Cfd Method to Predict Brownian Particle Deposition in a Constricted Tube," Particuology '''15''', 94 (2014).