Recent Progress on Mesh-free Particle Methods for Simulations of Multi-phase Flows: A Review †

The discrete element method (DEM) and the moving particle semi-implicit (MPS) method are the most popular mesh-free particle methods in the discontinuum and continuum. This paper describes a state-of-the-art modeling on multi-phase flows using these mesh-free particle methods. Herein, a combinational model of the signed distance function (SDF) and immersed boundary method (IBM) is introduced for an arbitrary-shaped wall boundary in the DEM simulation. Practically, this model uses a simple operation to create the wall boundary. Although the SDF is a scalar field for the wall boundary of the DEM, it is useful for the wall boundary of the CFD through combination with the IBM. Validation tests are carried out to demonstrate the adequacy of the SDF/IBM wall boundary model. Regarding the mesh-free particle method for continuum, the phase change problem is one of the challenging topics, as the solid state is usually modeled by extremely high viscous fluid in the phase change simulation. The phase change simulation is shown to be efficiently performed through an implicit algorithm and a heat flux model in the MPS method. The adequacy of these models is verified by the numerical examples.


Introduction
The discrete element method (DEM) (Cundall and Strack, 1979), the smoothed particle hydrodynamics (SPH) (Monaghan, 1994) and the moving particle semi-implicit (MPS) method (Koshizuka and Oka, 1996) are well known as a mesh-free particle methods for simulations of multi-phase flows. Although they require different targets, these methods are capable of easily simulating the violent motion of the media. The DEM, for instance, has been employed for simulating discontinuum as in the case of granular media and is known to inevitably and easily introduce inter-particle force. Consequently, the DEM can introduce cohesive forces, such as van der Waals force (Sakai et al., 2012) and liquid bridge force (Dhenge et al., 2013;Tsunazawa et al., 2016;Chan and Washino, 2018;Schmelzle and Nirschl, 2018;Sun and Sakai, 2018;Sakai et al., 2019). Very recently, modeling of the particle shape has been extensively studied in the DEM (Zhou et al., 2011;Govender et al., 2019;Kafashan et al., 2019;Nie et al., 2019;Shrestha et al., 2019;Yang et al., 2019). In contrast, the continuum media have been simulated using the SPH and the MPS methods, which have been consistently employed in the computation of free surface fluid flows, such as droplet impingement (Xiong et al., 2010(Xiong et al., , 2011, free surface fluid flow through porous media , discharge of highly viscous fluid flow (Sun et al., 2012), and wave migration (Shibata et al., 2011). Hence, the mesh-free particle methods make it possible to simulate complex multi-phase flow systems. From the background, this article highlights a discussion on the recent progress of these mesh-free particle methods from the numerical simulation viewpoint of multi-phase flows.
Regarding the multi-phase flows involving solid particles, the DEM-CFD method (Tsuji et al., 1993) has been widely employed, e.g., a fluidized bed Yue et al., 2019), a pneumatic conveyor (Klinzing, 2018), and a particle filtration technology (Kanaoka, 2019). With this method, the local volume average technique (Anderson and Jackson, 1967) has been introduced into the governing equations for the fluid phase. Moreover, staggered grids are usually employed in the fluid flow calculation, where the grid size is sufficiently larger than the particle size. Usage of the staggered grid makes it feasible to simulate the fluid flow precisely and stably. On the other hand, in the existing DEM-CFD method, modeling of an arbitrary-shaped wall and moving wall is substantially impossible, because of the rectangular shape of the staggered grid and the fixed grid points. A curved wall boundary cannot be created by the staggered grid, and no useful model exists for the moving wall boundary. Apart from the problem regarding the wall boundary, industrial-scale DEM-CFD simulation has not been performed on a single PC, though latest multi-core processors such as central processing unit (CPU) and graphical processing unit (GPU) have been developed (Shigeto and Sakai, 2011;Govender et al., 2018). This is because the DEM-CFD simulations cannot be finished within practical calculation time when the number of the computational particles approaches to the actual one. From this fact, the number of calculated particles is substantially restricted. Performing a calculation on a single PC is quite important for industrial engineers, because they can perform the calculation by themselves without any expertise regarding parallel computation and without introducing supercomputers. In summary, the existing DEM-CFD method has setbacks from the viewpoint of creating arbitrary shape wall boundary and due to a substantial limit of the number of calculated particles.
Phase change is one of the significantly important and challenging problems in numerical modeling of particle technology, e.g., additive manufacturing (Steuben et al., 2016;Russell et al., 2018) and thermo-mechanical behavior of powder (Krok et al., 2016). The MPS method can simulate phase change at ease, and hence has been applied to various phase change problems (Chen et al., 2014;Duan et al., 2019). In the modeling, enthalpy has been introduced to simulate the phase change, and the solid state is simulated by extremely highly viscous fluid. On the contrary, heat flux has not been given on free surface in the MPS method; instead, temperature supplied in the previous studies. Eventually, heat balance has not been assured during the heat transfer simulation using the MPS method.
In this manuscript, state-of-the-art technologies for the mesh-free particle methods, especially the DEM-CFD and MPS methods, are introduced, where these have been developed to address solutions to the above problems.

Arbitrary shape wall boundary model in DEM-CFD method
A useful wall boundary model for the DEM-CFD method is described herein, utilizing mainly a combination of the signed distance function (SDF) and the immersed bou nd a r y method (I BM ). For a bet ter understanding of the combinational technique, the SDF is first introduced briefly, and then the combined SDF-IBM is described.

Wall boundary model for the DEM
The SDF has been originally used in the level-set method (Osher and Fedkiw, 2002) in the CFD, and represents an implicit surface in space. Likewise, it has been shown being utilized in the DEM (Yokoi, 2005). The SDF consists of distance and sign; particularly, the distance from the nearest wall is calculated and saved for each point in a grid enclosing the boundary shape. When the computational particle is inside the calculation domain, the sign becomes positive, and vice versa.
Utilizing the SDF makes the computation of interaction forces between solid particles and wall boundary an easy task, as it does not necessitate particle contact detection for the surfaces, edges, and vertices comprising the wall, as compared to the conventional mesh-based wall boundary models. Once the distance from the particle center to the wall is obtained, the overlap between them can be easily estimated from the neighboring reference points.
Nevertheless, from a view point of energy conservation, the existing SDF model might lead to instability in non-dissipative systems. The new SDF model (Shigeto and Sakai, 2013) has been proposed to solve this problem. In this model, the normal component of the contact force is derived from the elastic potential energy; as such, this model enables the particle-wall interaction to retain the energy conservation in a non-dissipative system.
The SDF, which is denoted by ϕ, is given by where d and s are, respectively, the distance from the nearest wall boundary, and the sign that indicates the inside or outside of the calculation domain. The normal vector in the particle-wall interaction n pw and normal compact of overlap δ pw n are given by where r p is the particle radius. Finally, in the new SDF, the normal component of the contact force between the particle and the wall is given by The tangential component of the contact force is obtained in the same way as the traditional DEM as where v rt and μ are the tangential components of relative velocity of the interacting particles and friction coefficient, respectively. Note that ϕ is not necessary in the tangential component of the elastic force.
A number of validation tests have proven the adequacy of the SDF model, specifically, for a ribbon mixer (Basinskas and Sakai, 2016b), a pot blender (Basinskas and Sakai, 2016a), a twin-screw kneader , and powder die-filling (Tsunazawa et al., 2015).

Wall boundary model for the CFD
Some CFD techniques have been proposed for modeling arbitrary-shaped wall boundaries, one of which adopts an unstructured grid (Adam et al., 2016;R. Hu et al., 2019), where meshes are generated along the boundary. Although this technique is well established, it poses a problem on the mesh size, namely, the mesh size becomes extremely small when highly complicated boundaries are modeled. In addition, treatment of a moving wall is extremely difficult.
The IBM can be employed to overcome these problems. With the IBM, the boundary profiles are projected to the non-conforming computational grids. As such, Cartesian grids can be used even with the existence of complicated walls. The IBM eliminates the necessity of mesh generation along the boundary. Thus, the IBM is more efficient than boundary-fitted meshes. Fig. 1 shows schematically the boundary-fitted coordination system and the IBM.
Incidentally, the IBM falls under two categories (Lu et al., 2016). One is the so-called "continuous forcing" (Peskin, 2002) approach, in which control points are distributed over the boundary surface to modify the velocity field. The other approach is called "direct forcing" (Kajishima et al., 2001;Udono and Sakai, 2017), in which the boundary profile is projected directly to the fluid grid. Moreover, transfer of wall information to the computational grid plays a dominant role, for which an efficient method of projection is required.
The numerical technique combining the IBM with the SDF geometrical representation has been proposed (Xiaosong Sun and Sakai, 2016), with an objective of obtaining the local boundary profile easily. In the multiphase flows, the SDF provides a unified description of wall boundary for both the solid and fluid phases, defining superior computational efficiency. The numerical algorithm is shown briefly as follows.
Governing equations for the fluid phase include the continuity and Navier-Stokes equations, where the local volume average technique is introduced. These equations are expressed as follows: where ε, u, ρ, p, f s , τ, and ɡ are, respectively, the void fraction, fluid velocity, fluid density, fluid pressure, solidfluid momentum exchange term, viscous stress, and gravitational acceleration. Combined Ergun (Ergun, 1952) and Wen-Yu (Wen and Yu, 1966) equation is often used in the evaluation of the drag force acting on a computational particle.
The IBM can simulate the solid-fluid interaction and even utilizes structured grids, for which the object size should be sufficiently larger than the mesh size. With the IBM, the volume-weighted average velocity is used to analyze the interaction between the fluid and the solid object. The volume-weighted average velocity U is calculated by In Eq.(8) α and u obj indicate the local volume fraction (LVF) of the object and the object velocity. Here α can be determined efficiently by counting the saved points of the SDF in a CFD grid. Note that this volume-weighted velocity is used instead of fluid velocity, and the whole calculation domain is solved as a fluid, even inside the solid object.
When fluid velocity u is replaced by the volumeweighted average velocity U, Eq. (7) can be updated as where F IB is the correction term in the IBM. The force term F IB is given by Very recently, the authors have shown that fluid flow might not be calculated suitably when the IBM is used as the wall boundary. This problem can be solved by considering the density scaling coefficient (Sun and Sakai, 2017), which is given by where γ is the density scaling factor. Actually, the density is scaled in the pressure Poisson equation, where the grids are filled by the object. Eventually, the pressure Poisson equation is be given by with the density scaling factor being introduced. In Eq.
(12), superscript n + 1 and n are updated iteration number and present iteration number.

Powder flow involving airflow in die-filling
As a typical gas-solid flow system with a moving wall boundary, the powder die-filling is a preferable example. To date, fundamental studies of die-filling (Zakhvatayeva et al., 2018;Baserinia and Sinka, 2019) have been extensively focused in pharmaceutics. As far as the pharmaceutical particles are concerned, not only the particle size is small but also their density. Hence, the drag force acting on these particles is exceptionally high. This fact implies that solid particle behavior can be precisely simulated if solid-fluid interaction, gas flow, and moving wall boundary are precisely modeled.
In the authors' group, the Advanced DEM-CFD method (Yao et al., 2018) was proposed for investigating the effect of air flow on the particle filling state in die-filling. In the Advanced DEM-CFD method, combination of the SDF and the IBM was introduced into the DEM-CFD method. The schematic diagram of the die-filling system is shown in Fig. 2; here, the shoe was moved by 0.10 m/s toward the left side. In this study, calculations and experiments were performed under fair conditions.
The target powder was NONPAREIL-108 ® (FREUND CORPORATION). The particle size and the particle density were, respectively, 200 μm and 1,379 kg/m 3 , and the number of the calculated particles was 150,000. A combination of the SDF and IBM was used for the wall boundary model for the DEM and the CFD. Fig. 3 displays the schematic diagram of the die-filling system. Fig. 3(a) is the scalar field of the SDF, where the sign of the calculation domain and inside of the objects respectively became positive and negative. The SDF was created independently to each object. The interval of the SDF was usually one-quarter of the particle size. Fig. 3(b) designates the LVF for the IBM. When the object was fulfilled in a CFD grid, the value of the LVF of the gas phase, namely, 1 -α, became zero. By calculating the values of the LVF, interactions between the object and fluid could be simulated through the IBM. In order to show the influence of airflow on the flow of the solid particles, two kinds of simulations were performed: a DEM simulation for a single granular flow and a DEM-CFD simulation for a gas-solid flow. The calculation results are presented in Fig. 4. Die-filling simulation by the DEM-CFD method, where air flow was taken into consideration, demonstrated the appearance of a clearance on top of the powder bed, as illustrated in Fig. 4(a). On the contrary, the powder was fulfilled when the DEM was simulated without air flow.
Actually, the powder bed shape obtained from the DEM-CFD simulation exhibited an excellent agreement with that of the experiment. Hence, consideration of the air flow was shown as an indispensable factor for the simulation of the die-filling in pharmaceutical powders.

Industrial fluidized bed
A fluidized bed is extensively investigated by numerical simulation, where the DEM-CFD method is widely employed. At present, application of the DEM-CFD method to industrial fluidized beds (Sakai, 2016;C. Hu et al., 2019;Ge et al., 2019;Stroh et al., 2019) is among the important research topics. Not the shape of industrial fluidized beds is always simple, though rectangular systems have been employed in the existing fluidized bed simulations. Nonetheless, the structure of actual fluidized beds might be extremely complex, as tubes might have been inserted, a slope might have been set at the bottom, or they must have been equipped with cyclone. As a matter of course, number of the solid particles becomes huge in the industrial systems. In this sense, besides modeling of the calculation domain shape, the number of the computational particles becomes a critical problem as well. The authors' group attempted to solve the problem with the limit of the calculated particles by developing a coarse graining DEM (Sakai and Koshizuka, 2009;Sakai et al., 2010;Sakai et al., 2012b;Sun et al., 2014;, in which largely modeled particles represent the group of the original particles where total energy between them is assumed to agree. The adequacy of such approach has been shown through verification and validation tests, in which the combined SDF-IBM was employed in the numerical simulations for the fluidized bed. Very recently, the compatibility of the coarse graining DEM and the wall boundary model of SDF/IBM has been proven through validation tests (Mori, Wu and Sakai, 2019), whose results are addressed in this paper.
A schematic diagram of the fluidized bed system is illustrated in Fig. 5. In the fluidized bed, four tubes were inserted. Gas was injected at 0.2 m/s from the bottom side; glass beads were used for the powder and gas for air. The coarse graining DEM was employed in this study; hence, the number of computational particles could be drastically reduced than the actual one. Namely, the coarse grain ratio was 5.0 for 250,000 calculated particles. Actually, the number of solid particles reached the order of seven digits, and hence, calculation could not be performed without the coarse graining DEM on a single PC. Hence, application of the coarse graining DEM makes it possible to perform large-scale gas-solid systems with a smaller number of the calculated particles, as compared to using the actual number. The SDF and the IBM were employed as the wall boundary model, as depicted in Fig. 6. In Fig. 6(a), the calculation domain was expressed as positive in the SDF. LVF of the gas phase (1 -α) was zero inside of the object in the IBM, as illustrated in Fig. 6(b).    7 designates snapshots obtained from the validation tests. Note the excellent agreement between the calculated and experimental results. Specifically, bubble behavior and spatial location of the particles in the simulation and in the experiment were in good agreement. Besides, the pressure drop quantitatively agreed in the computation and the experiment. Thus, the coarse graining DEM was shown capable of simulating the macroscopic behavior of the original gas-solid flow system, although the particle size was larger than the actual one. This is because the total energy was modeled to agree between the coarse grain particle and a group of the original particles. In fact, the calculation has been done prior to the experimentation. Surprisingly, the agreement could be obtained without any necessary adjustment in the calculation.

Mesh-free particle method for a phase change problem
Modeling the phase change between the solid and liquid phases is one of the challenging problems in the CFD. The mesh-free particle method is a promising approach to simulate the phase change, and has been applied into various systems, such as additive manufacturing (Russell et al., 2018), welding , and nuclear severe accidents (Duan et al., 2019). During a phase change, the solid and liquid states are mixed, thus, such state is regarded as a solid-liquid flow in a broad sense. In the simulation of the phase change by the mesh-free methods, the efficient calculation for highly viscous fluid flow, as well as supplying heat flux on the free surface, is essential. The authors' group has developed certain models for the highly viscous fluid flow and the heat flux on the fluid surface, which are described herein.

The MPS method for phase change simulation
The original MPS method (Koshizuka and Oka, 1996) has been developed for an incompressible fluid flow, where the pressure Poisson equation is calculated implicitly. Recently, weakly incompressible algorithm has been introduced into the MPS method by calculating the pressure explicitly. Such method is referred to as E-MPS method. Here, the kernel functions are the same as the original ones. In general understanding, the accuracy of the E-MPS method is quite similar to that of the original MPS method, as both have been employed even in the solid-liquid flow states in chemical engineering (Sakai et al., 2012;Yamada and Sakai, 2013;Sun et al., 2014). Modeling of the phase change between the solid state and and   , where μ f , σ, κ, δ and n, h, λ and Q are, respectively, fluid viscosity, coefficient of surface tension, curvature, delta function, normal vector, enthalpy, heat conductivity and heat source. In Eq. (15), the heat flux can be given as the volumetric heat source, with q as the heat flux. Moreover, in the MPS method, the pressure and viscous terms in the governing equations are modeled by a gradient model and a Lapracian model. Regarding fluid particle i, these models are given by where D, n 0 , r and Ω are, respectively, dimension, initial particle number density, position vector, and Lapracian constant in the MPS method. In these models, weight functions are used to give smoothed interpolation for quantities through the interaction between fluid particles. In order to perform a phase change, the enthalpy is introduced as, where h s , C p , T m , h l , and L are, respectively, enthalpy of the solid phase, heat capacity, melting temperature, enthalpy of the solid phase, and latent heat. Incidentally, solid state is simulated as exceptionally highly fluid expressed by Arrhenius viscosity formula. When the viscous term is calculated implicitly, updated fluid velocity u k+1 can be obtained by the expressions, * 1 , A detailed discussion on the MPS method is presented in the original papers (Koshizuka and Oka, 1996).

Investigation on the calculation efficiency of a highly viscous fluid
First, implicit calculation of the viscous term is shown to be efficiently performed in the simulation of a discharged fluid flow in a glass melter.
In this system, kinematic viscosity, fluid density, and surface tension coefficient were, respectively, 0.03 m 2 /s, 2500 kg/m 3 , and 0.3 N/m. The highly viscous fluid was flown by only the gravitational force. Fig. 8 illustrates the snapshots obtained from a simulation of a discharged flow of the highly viscous fluid in a glass melter (elapsed time: 5.0 s, 10 s, 15 s and 20 s). Here, the volume flow rate in the simulation agreed well with the experiment, in the steady state. Subsequently, the acceleration due to the implicit model was investigated based on the maximum value of the time step. Specifically, the time step ratio of the Courant condition to diffusion number was roughly 10 -4 , which implies that the computation could be accelerated up to 10 4 times faster than the original viscous calculation model in this system. Thus, the efficiency of the implicit calculation algorithm was proven in the highly viscous fluid flow system.

Adequacy of the heat flux model in a meshfree particle method
Subsequently, the adequacy of the heat flux model is mentioned in the MPS simulation. In the MPS method, instead of the heat flux, the temperature is usually given (the Dirichlet boundary condition) at the free surface due to ease. On the contrary, heat balance cannot be kept if the heat flux (the Neumann boundary condition) is not given. As shown in Eq. (16), the heat flux can be introduced into the mesh-free particle method as heat source. This attempt has been done very recently by the authors (Takabatake et al., 2016).
The adequacy of the heat flux model was examined via heat transfer simulation in a simple system (Fig. 9).
As boundary conditions, heat flux was 25,000 W/m 2 at the top of the rectangle, while the temperature at the bottom was 300 K. Regarding the physical properties, the density, kinetic viscosity, specific heat, and thermal conductivity were, respectively, 1,000 kg/m 3 , 1.0 × 10 -3 m 2 /s, 10 J/(kg K), and 50 W/(m K). Temperature distribution in the calculation and theoretical results were later compared.
The calculation result is shown in Fig. 10. Temperature increased gradually from the initial value (300 K). The temperature distribution achieved a steady state in 20 s and was found to quantitatively agree with the theoretical results. This result implies that the heat flux can be given even in the mesh-free particle method, where the free surface is violently moved.
This aside, a phase change simulation was performed through the supplied heat flux. Two types of materials, namely, Material-l and Material-2, were used in the calculations. The density, kinetic viscosity, specific heat, thermal conductivity, latent heat and melting temperature of Material-1 were, respectively, 1,000 kg/m 3 , 1.0 × 10 -3 m 2 /s, 100 J/kg K, 50 W/m K, 10,000 K, and 1.0 × 10 4 J/kg. Material-2 had the same density, kinetic viscosity, and latent heat, while the specific heat, thermal conductivity, and melting temperature were, respectively, 10 J/kg K, 500 W/m K, and 320 K. As far as the boundary conditions were concerned, the initial heat flux was 100,000 W/m 2 on the upper side and the temperature was 300 K at the bottom side.
Snapshots obtained from the phase change calculation are shown in Fig. 11. Due to the heat flux, the upper side was molten at the beginning. Solidification occurred as the molten fluid approached the floor, as the floor temperature was kept at 300 K. Hence, simulation of the phase change was demonstrated by the MPS method.

Conclusions
This paper described a cutting-edge modeling on simulations for multi-phase flows using the DEM and the MPS methods.
Relative to the DEM, a brand-new wall boundary model was introduced, namely, a combinational model of the SDF with the IBM, which was shown to be useful in generating an arbitrary shape wall boundary in the DEM-CFD simulation. Such new model is advantageous in the sense that an arbitrary shape wall boundary could be created by a simple operation, and thus, it might be standardized in the future DEM-CFD simulations.
Regarding the mesh-free particle method for continuum such as the MPS method, the phase change problem has generally been a challenging topic. An implicit calculation method, along with the heat flux model, was developed for efficiently simulating the phase change. The usefulness of such heat transfer modeling was demonstrated in this paper. Hereafter, it is expected to contribute to the numerical simulation of a thermo-mechanical problem.

Acknowledgments
Parts of this study were financially supported by JSPS KAKENHI Grant Numbers 17H02825 and 17KK0110.