Numerical Simulations of Liquid-Solid Flows in A Vertical Pipe by MPS-DEM Coupling Method

In the process of deep-sea mining, the liquid-solid flows in the vertical transportation pipeline are very complex. In the present work, an in-house solver MPSDEM-SJTU based on the improved MPS and DEM is developed for the simulation of hydraulic conveying. Firstly, three examples including the multilayer cylinder collapse, the Poiseuille flow and two-phase dam-break are used to validate the precision of the DEM model, the pipe flow model and MPS-DEM coupling model, respectively. Then, the hydraulic conveying with coarse particles in a vertical pipe is simulated. The solid particle distribution is presented and investigated in detail. Finally, the coupling method is successfully applied for the simulation of the liquid-solid flows in a vertical pipe with rotating blades, which shows the stability of the solver under rotating boundary conditions. This fully Lagrangian model is expected to be a new approach for analyzing hydraulic conveying.


Introduction
With further exploitation of deep-sea resources, researchers have paid more attention to the hydraulic conveying (Ali and Yeung, 2014;Silva et al., 2015). The minerals and the seawater are transported from the seabed to the floating bodies through the pipelines (Zhou et al., 2019). The vertical pipes and the rotating machineries tend to suffer from continuous abrasion due to the impact of the gravels, which may affect the safety of offshore operation. Therefore, it is important to investigate the distribution of solid particles and predict the attrition of the vertical pipe.
The fluid-solid flows' pattern and the distribution of solid particles in the pipelines can be recorded with experimental techniques accurately. However, the experimental equipment such as PIV and high-speed camera is very expensive. In recent years, several numerical techniques have been proposed and applied for fluid-solid two-phase problems, which can provide more detailed information of the flow field. The solid part of two-phase is modelled as the continuum (Khanpour et al., 2016;Tajnesaie et al., 2018) or the discrete bodies (Feng and Yu, 2004;Zhou et al., 2010). The former is so called Two Fluid Model (TFM).
The solid phase is also governed by continuity equation and the momentum equation as the fluid phase, which is suitable for the large-scale simulation, such as the sediment transport (Shi et al., 2019) and landslides (Jandaghian et al., 2021). The solid phase is regarded as a whole, while the movement of individual solid particle and the internal collision between solid particles are ignored. In addition, it is difficult to model the collision between particles and walls, especially moving boundaries (Delacroix et al., 2021). Discrete Element Method (DEM) was firstly proposed by Cundall and Strack (1979) and proven to be effective to simulate the individual behavior of each solid body. DEM has been used to solve many problems, such as landslides (Xu and Dong, 2021), structure deformation (Wang, 2020), and ice breaking (Liu and Ji, 2021;Long et al., 2020).
The DEM method can be coupled with other method to simulate the solid−liquid flows. The liquid phase can be solved by mesh-based methods or particle-based methods. The CFD-DEM as a Eulerian−Lagrangian model has been widely used to investigate the problems of pipeline transportation. Zhou et al. (2019Zhou et al. ( , 2020Zhou et al. ( , 2021 used unresolved CFD-DEM to simulate the hydraulic conveying in the vertical pipe, horizontal pipe and bends. The solid particles' distribution, flow regime and erosion rate are investigated in detail. Filho et al. (2022) conducted a resolved CFD-DEM simulation based on the open-source code CFDEM (Kloss et al., 2012) and investigated the hydraulic conveying of coarse grains through a bending. In the cases with highspeed rotating boundaries, the meshes are prone to be distorted and remeshing is needed to avoid the ill meshes, which affects the accuracy and efficiency of the simulation.
Smoothing Particle Hydrodynamics (SPH) method (Gingold and Monaghan, 1977) and Moving Particle Semi-implicit (MPS) method (Koshizuka and Oka, 1996) based on the lagrangian view are two main particle-based methods. Compared with original weakly compressible SPH, the pressure field and volume conservation of MPS is enhanced. However, the solution process of MPS is time-consuming and is parallelized with more efficiency loss than SPH . There is no fixed topology for particles, the particlebased method shows its superiority in the simulation of moving boundaries (Zhang and Wan, 2018). The particlebased methods developed rapidly in the past decades. The gradient model of momentum conservation (Tanaka and Masunaga, 2010), the mixed source term of Pressure Poisson Equation (PPE) (Tanaka and Masunaga, 2010;Khayyer and Gotch, 2011) and the Particle Shifting Technique (PST) (Khayyer et al., 2017;Duan et al., 2018) have been proposed, which significantly improve the stability and accuracy of incompressible particle-based method. Based on the above improvement, particle-based methods have been applied for the simulation of multiphase flows (Wen et al., 2021a(Wen et al., , 2021b(Wen et al., , 2021cXie et al., 2021b) and Fluid-Structure Interaction (FSI) (Gotoh et al., 2021;Xie et al., 2021a;Zhang et al., 2021). In the actual pipeline system, there are pumps rotating fast to lift the seawater and minerals. In this work, a simplified pump with blades is placed in the vertical pipe and the particle-based method is used to simulate the fluid and moving boundary.
The coupling of MPS/SPH and DEM can be classified as resolved coupling and unresolved coupling. According to the concept of resolved coupling, DEM particles are driven by the pressure and viscous force. They serve as the boundaries for the fluid (Zhang et al., 2009;Zhou et al., 2022). In order to fit the curve boundary, the selected particle spacing of boundary is small. Therefore, the computational cost is heavy if simulation is carried out with uniform resolution. The unresolved coupling is based on the local average technique (Anderson and Jackson, 1969). DEM particles can overlap with fluid particles and it is unnecessary to use extra boundary particles to prevent the penetration of fluid particles. In addition to the gravity and contact force, DEM particles are moved by the empirical drag force and buoyancy from the fluid. Simulations of water entry of a sphere are carried out by SPH−DEM resolved coupling method (Zhou et al., 2022) and the SPH−DEM unresolved coupling method (Xu et al., 2019). The SPH-DEM unresolved coupling method with low resolution even performs better than the resolved coupling method with higher resolution. It is obvious that the unresolved coupling method has higher computational efficiency than the resolved coupling method. The unresolved coupling method has been employed to solve the liquidsolid two-phase problems. Robinson et al. (2014) used an SPH−DEM method to simulate the sediment of single particle and multiple particles. The influence of SPH resolution, drag laws, porosity and fluid physical properties on the numerical results were investigated. Xu and Dong (2021) coupled SPH and DEM to investigate the effect of different type of slides to the generation of tsunamis.  developed an SPH−DEM−FEM model to investigate the interaction between the fluid-particle flows and elastic structures. The simulation of Wenjia gully debris flow was conducted and the destruction of the structures was estimated. Kim et al. (2021) coupled ISPH and Coarse-Grained DEM (CGDEM) to investigate the scouring behaviors. Compared with original DEM, the computational cost of CGDEM was very small. The local volume fraction was not considered in fluid phase, which had little effect on the results of large-scale simulation. Although the MPS−DEM and SPH−DEM methods have been widely used to simulate liquid-solid two-phase flow with free surface, they are rarely used to solve the problem of hydraulic conveying in the pipe.
In this study, a 3-D MPS−DEM coupling method is developed and firstly applied for simulating the fluid−solid flow in the hydraulic conveyer with rotating boundary. The paper is organized as follows. Firstly, an improved MPS method, a DEM method, their coupling scheme and rotating boundary condition are presented briefly. Then, three examples including the multilayer cylinders collapse, the Poiseuille flow and liquid−solid two-phase dam-break are simulated to validate the precision of the DEM model and pipe flow model and MPS−DEM coupling model respectively. Finally, two numerical tests are conducted. The fluidsolid flows in a vertical pipe are simulated. The distribution of solid particles under different velocities of fluid is investigated. In another test, a rotor with blades is added to the vertical pipe. The influence of rotor on the fluid pattern and the distribution of solid particles during the conveying process is discussed.

Liquid phase
The local average technique is used to couple the MPS and DEM, which was proposed by Anderson and Jackson (1969). The governing equations considering the local volume fraction are given by ∂ ∂t where subscript f represents the fluid particle. , , , , , and t are the local volume fraction, the density, the velocity vector of fluid, the pressure, the viscosity, the gravity acceleration vector and the physical time, respectively.
is the reaction force from the solid particles.

W (r)
The interactions among MPS particles are governed by the kernel function , which plays a role of weight function in the discretization process. The kernel function proposed by Zhang et al. (2014) is adopted here, which can avoid numerical singularity.
r r e r e r 0 r 0 where denotes the distance between two particles and denotes the radius of particles interaction. is set to 4.1 , being the initial particle spacing. The number density represents the degree of aggregation of MPS particles, which is directly proportional to the density of the fluid, calculated as: ⟨n⟩ where subscripts i and j represent the interacting MPS particles and r denotes the position vector. For the simulation of single-phase incompressible fluid, the particle number density should keep the same with the initial particle number density . To solve the problem of solid-fluid flows, the particle number density required to be a constant value considering the local volume fraction (Sakai et al., 2012). ⟨n⟩ n 0 n 0 In the particle interaction models, is simply replaced by . The particle interaction models consist of gradient model, divergence model and Laplacian model. ⟨∇ϕ⟩ where indicates the physical quantities of the MPS particles, represents the number of space dimensions, is the parameter used to make the variance increase to be equal to the analytical solution (Koshizuka and Oka, 1996). Through solving the Pressure Poisson equation (PPE), the pressure field can be obtained. A mixed source term method (Khayyer and Gotoh, 2011;Tanaka and Masunaga, 2010) is employed to solve PPE, which is defined by is the pressure of the next step, is the blending parameter and its value varies from 0 to 1, is the time step, and denote the intermediate velocity and particle density, respectively. In this work, the value of is set to 0.01.

Solid phase
DEM method is used to simulate the movement of solid particles, which is based on the Newton's Second Law, given by where the subscripts k and l represent the solid particles contacting with each other; , , and represent the mass, the moment of inertia, the translational and rotational velocities of the DEM particle k, respectively; and are the contact force and moment due to the solid− solid contact, respectively; is the hydrodynamic force applied to the DEM particles.
The DEM particles are regarded as soft spheres (Cundall and Strack, 1979) and the contact force is decided by the size of overlapping parts. The contact model is made up of dashpots, springs and sliders. The contact force can be decomposed into the normal and tangential components. The normal force is given by, where and denote the relative displacement and velocity in normal direction between DEM particle k and particle l. and are the spring stiffness and damping coefficient, respectively. The tangential force is given by where and denote the relative displacement and velocity in tangential direction between DEM particles k and l; is the friction coefficient. The damping coefficient is given by where e is the restitution coefficient.
544 XIE Feng-ze et al. China Ocean Eng., 2022, Vol. 36, No. 4, P. 542-552 2.3 Two-phase interaction The fluid forces applied to DEM particles include the drag force, pressure gradient force, virtual mass force and lubrication force. Because the drag force and pressure gradient force dominate (Sakai et al., 2012), the other forces are not considered in this work. The coupling strategy of twophase interaction is presented in Fig. 1.  The empirical drag force can be calculated as: where and are the volume and the local volume fraction of solid particle k. is the fluid velocity in the center of solid particle k. Two drag models of Ergun (1952) and Wen and Yu (1996) are combined and the coefficient of the interphase momentum exchange is given by where is the diameter of the DEM particle k. is the drag coefficient, which is related to the Reynolds number, written as: , The Reynolds number in the center of DEM particle k is defined as: An additional weight function is employed to build the fluid−particle interaction and it is developed by Sun et al. (2014), given by The local volume fraction of fluid particle i is calculated as: The local volume fraction and fluid velocity in the center of DEM particle k are calculated by neighboring MPS particles with the weight function , given by .
The drag force and the pressure gradient force can be combined and the total hydrodynamic force applied on the DEM particles can be written as: The reaction force exerted on the MPS particles is required to satisfy the momentum conservation law, given by where is the volume of the MPS particle.

Boundary condition
The solid particles may collide with the rotor during the conveying process. Therefore, the collision model should be added to the rotor. The MPS−DEM boundary condition of the pipe is the same as that in the reference (Xie et al., 2021b), which is time-saving. The shape of the rotor is not as regular as the pipe and another approach is developed for the pipe. Fig. 2 presents the sketch of MPS−DEM boundary for the rotor. DEM boundary made of DEM particles completely overlaps the MPS boundary. However, the forces applied to DEM boundary particles are not calculated. The velocity and the displacement of DEM boundary particles are updated according to their corresponding MPS boundary particles. The calculation of boundary force on DEM particle flow follows the common solid-solid interaction. Fig. 3 shows the flow chart of MPS−DEM coupling method. The whole algorithm contains the MPS part for the XIE Feng-ze et al. China Ocean Eng., 2022, Vol. 36, No. 4, P. 542-552 simulation of fluid and the DEM part for the simulation of solid. Firstly, the program reads the properties and the initial position of MPS−DEM particles from the input file. Then, the main loop starts. If liquid-solid flows are simulated, DEM calculation is performed after MPS particle neighbor list is created. Meanwhile, the position and the velocity of the DEM boundary particles are transferred from MPS solver to DEM solver. The DEM program can be divided into two parts: MPS−DEM interaction and solid−solid calculation. In the MPS−DEM part, the fluid-particle neighbor list is created and the fluid−solid interaction, including the voidage and internal forces, are calculated. The solid−solid calculation follows creating solid-solid list, calculating contact force and updating position/velocity of DEM particles. The DEM solver cycle multiple times until the DEM physical time match the MPS physical time. The interaction parameters of voidage and internal forces are transferred from DEM solver to MPS solver. In addition to the creation of MPS particle neighbor list, the computational flow chart of MPS method follows explicitly calculating gravity and viscous force, updating temporal velocity of MPS particles, calculating particle number density, detecting free surface, solving the pressure Poisson equation, calculating pressure gradient, updating position and velocity of MPS fluid particles and updating position/velocity of MPS−DEM boundary particles.

Verification of the DEM model
In this sub-section, the experiment carried out by Zhang et al. (2009) is used to validate the DEM part of the MPS-DEM-SJTU solver. Multi-layer aluminium cylinders are initially arranged as a hexagonal packing scheme at the left of the tank as shown in Fig. 4. The cylinders have a diameter of 0.01 m and a density of 2700 kg/m 3 . The tank is made of PVC and has a length of 0.26 m and a width of 0.1 m. There is a movable plate near the right-most cylinders, which is 0.06 m from the left wall of the tank. At the beginning of the experiment, the plate moves upward immediately allowing cylinders to collapse under the gravity. Other parameters of the simulation are presented in Table 1.   Fig. 5 presents the comparison between the experimental and numerical results. It can be noticed that the distribution of cylinders observed in experiment and simulation is almost the same at different instants. The quantitative comparison is further conducted and the average location of the mass center (Zhang et al., 2009) is recorded as shown in Fig. 6. It can be noted that the numerical results agree well with the experimental data. In summary, the DEM solver can accurately predict the collision among the solids.

Verification of the pipe flow model
Poiseuille flow (Morris et al., 1997) is a Laminar flow in pipes, which is driven only by pressure. In this sub-section, Poiseuille flow is simulated by MPS. The pipe has a radius of 0.1 m and a length of 1 m. The velocity of the inlet flow is 0.002 m/s and its corresponding Reynold number is 400. More detailed parameters of the simulation can be found in Table 2. Fig. 7 shows the distribution of MPS particles at different instants. It can be noted that the distribution of MPS particles takes on the shape of a rotating parabola. The velocity distribution in the middle section of the pipe is shown in Fig. 8. In comparison with the analytical result, the maximum error of MPS results is smaller than 3%, showing that the present solver can be used to simulate the internal flow in the pipe.

Verification of MPS-DEM model
In this sub-section, an experiment (Sun et al., 2013) of liquid-solid two-phase dam-break is used to verify the MPS−DEM model. Fig. 9 presents the two-phase dam-break at t = 0 s. The tank is 0.2 m long and 0.1 m wide. A water column containing solid sediment is restricted to the right side by a plate. The water depth is 0.1 m. At the beginning of the simulation, the plate moves at a constant velocity of 0.68 m/s. Other parameters of the simulation are presented in Table. 3. Fig. 10 shows flow patterns of experiment and simulation at different instants. It can be noticed that both the shape of sediment and the profile of free surface simulated by MPS−DEM are in good agreement with experimental results. The 3-D effect can be observed in the results simulated by MPS−DEM as shown in Fig. 10d. The similar phenomenon can also be observed in the numerical results obtained by SPH−DEM (He et al., 2018). Fig. 11 presents the time his-tories of the leading front of liquid−solid flows. The numerical results are in general agreement with the experimental data. Some small discrepancies of the leading front of the sediment between the numerical results and experimental data can be observed. The main reason may be that the physical parameters of solid phase in simulation is roughly estimated (Sun et al., 2013;Markauskas et al., 2018). In general, liquid−solid two-phase can be simulated by MPS−DEM solver accurately.

Hydraulic conveying in a vertical pipe
The hydraulic conveying in a vertical pipe is simulated. The pipe is 4 m long and its radius is 0.1 m. Several DEM particles are placed in the pipe at the initial time as shown in    Fig. 12. Once the DEM particles rush out of the outlet boundary, they will be moved to the inlet boundary immediately and other physical parameters they carry will not change. Other physical parameters of the simulation are listed in Table 4. Two cases with inlet fluid velocity V f being 3 m/s and 5 m/s are simulated. Fig. 13 shows the distribution of solid particles when two-phase flows are fully developed. Several solid particle clusters can be observed during the hydraulic conveying process with the inlet fluid velocity being 5 m/s. When the inlet velocity is high, it can be noticed that there      are some gaps in the solid particle flow, while when the inlet velocity is low, it seems that the particle flow is relatively continuous.

Hydraulic conveying in a vertical pipe with a rotating rotor
In this sub-section, the MPS−DEM coupling model is applied to the simulation of hydraulic conveying in a vertical pipe with a rotating rotor. The rotor is placed 3 m away from the inlet of the pipe and rotates at 60 rpm. The inlet velocity of MPS−DEM two-phase flows is 2 m/s. Other parameters are similar to those in the simulation of hydraulic conveying in a vertical pipe without rotor. Fig. 14 shows the DEM particles' distribution.
The fluid velocity vectors in different pipe sections at t = 2.1 s is shown in Fig. 15, where colors represent the solid velocities in y direction. It can be noticed that the fluid rotates in the downstream region of the rotor. The movement of the fluid is also affected in the upstream region of the rotor, especially near the rotor. The velocity vectors of the fluid far away from the rotor point out the outlet of the pipe as shown in Fig. 15d. Fig. 16 presents the distribution of solid particles during the hydraulic conveying process. Solid particles collide with the rotor at around t = 1.6 s and their velocities are reduced immediately. Therefore, they gather near the rotating shaft of the rotor, where the tangential velocities of the rotor and fluid are relatively small. Fig. 17 shows the velocity vectors of solid particles. Owing to the existence of rotating rotor, the trajectory of solid particles under the tangential force of fluid and rotor changes. Solid particles bypass the rotor and are accelerated by the fluid again as shown in Fig. 17b. The solid particles are lifted spirally as shown in Fig.17d. Although the fluid has been affected before it comes into contact with the rotor, the motion of solid particles is not changed obviously. The main reason is that the solid particles tend to concentrate    near the axis of the pipe, where the circumferential component of fluid velocity vector is small. Figs. 18 and 19 present time histories of the solid-solid contact force exerted on the rotor and the pipe, respectively. The x component of the force exerted to the rotor increases with relatively small fluctuation as the solid flows come into contact with the rotor. When the solid flows fully contact the rotor, the fluctuation of the curve is enlarged. The force in x direction is obviously larger than the force in other directions. As the trajectory of particle flow is changed by the rotor, the pipe wall is also impacted by particles as shown in Fig. 18. The above analysis of fluid velocity, particle distribution and dynamitic loads show that MPSDEM-SJTU has the ability to solve the hydraulic conveying problems with rotating boundary.

Conclusions
In this paper, a fully Lagrangian method 3-D MPS− DEM is developed and firstly applied to the problems of    hydraulic conveying. Rigorous validation of MPS, DEM and MPS−DEM methods are carried out firstly. The multilayer cylinder collapse, the 3-D Poiseuille flow and liquidsolid two-phase dam-break are simulated. The numerical results agree well with the experimental and analytical results, showing the accuracy of numerical models. Then, the simulation of solid-liquid flows in a vertical pipe is conducted. The distributions of solid particles under different inlet boundary velocities are different. When the inlet fluid velocity is high, solid particle clusters can be observed during the hydraulic conveying process, but when the inlet velocity is low, it seems that the particle flow is relatively continuous. Finally, solid-liquid flows in a vertical pipe with a rotating rotor are simulated. Under the influence of the rotating rotor, the pattern of solid particle flows is changed and the solid particles are lifted spirally. The numerical results of those two cases are the same with the actual situation, which shows that MPS−DEM has the potential to solve the problems of hydraulic conveying with rotating boundaries.
In this work, the shape of the rotor is relatively simple, which is not as complex as the real pump. In order to match the model of the real pump, the MPS particle resolution in the whole computational domain will be very high, which increases the computational cost sharply. In the future, the acceleration techniques, including the GPU parallel technique (Xie et al., 2020) and the adaptive multi-resolution technique (Sun et al., 2017), will be introduced to the MPS−DEM solver and more realistic simulations of the hydraulic conveying process will be conducted.

Right and permission
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.