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

Abstract

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.

1. 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, 2011), free surface fluid flow through porous media (Sun et al., 2019), 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 (Sakai et al., 2010, 2014; 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.

2. 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 boundary method (IBM). For a better understanding of the combinational technique, the SDF is first introduced briefly, and then the combined SDF–IBM is described.

2.1 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

  

$$\phi (\mathit{x})=s(\mathit{x})d(\mathit{x}),$$(1)

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 δ_pwn are given by

  

$${\mathit{n}}_{\text{pw}}=\frac{\mathit{\text{∇}}\phi }{\left|\mathit{\text{∇}}\phi \right|}$$(2)

and

  

$${\mathit{\delta }}_{{\text{pw}}_{\text{n}}}=\left(\phi -r_{\text{p}}\right){\mathit{n}}_{\text{pw}},$$(3)

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

  

$${\mathit{F}}_{{\text{C}}_{\text{n}}}=\left(-k{\mathit{\delta }}_{{\text{pw}}_{\text{n}}}\left|\mathit{\text{∇}}\phi \right|-\eta {\mathit{\nu }}_{{\text{r}}_{\text{n}}}\right).$$(4)

The tangential component of the contact force is obtained in the same way as the traditional DEM as

  

$${\mathit{F}}_{{\text{C}}_{\text{t}}}=\{\begin{array}{ll}-k{\mathit{\delta }}_{{\text{pw}}_{\text{t}}}-\eta {\mathit{\nu }}_{{\text{r}}_{\text{t}}}\hfill & \left|{\mathit{F}}_{{\text{C}}_{\text{t}}}\right|<-\mu \left|{\mathit{F}}_{{\text{C}}_{\text{n}}}\right|\hfill \\ -{\mu }_{\text{S}}\left|{\mathit{F}}_{{\text{C}}_{\text{n}}}\right|\frac{{\mathit{v}}_{{\text{r}}_{\text{t}}}}{\left|{\mathit{v}}_{{\text{r}}_{\text{t}}}\right|}\hfill & \left|{\mathit{F}}_{{\text{C}}_{\text{t}}}\right|\ge -\mu \left|{\mathit{F}}_{{\text{C}}_{\text{n}}}\right|\hfill \end{array},$$(5)

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 (Sakai et al., 2015), and powder die-filling (Tsunazawa et al., 2015).

2.2 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.

Fig. 1

Overview of typical computational grid systems.

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; X. Sun and Sakai, 2016; 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 multi-phase 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:

  

$$\frac{\partial \varepsilon }{\partial t}+\mathit{\text{∇}}·(\varepsilon \mathit{u})=0,$$(6)

and

  

$$\begin{array}{l}\frac{\partial (\varepsilon \rho \mathit{u})}{\partial t}+\mathit{\text{∇}}·(\varepsilon \rho \mathit{uu})\\ =-\varepsilon \mathit{\text{∇}}p+{\mathit{f}}_{\mathbf{s}}+\mathit{\text{∇}}·(\varepsilon \mathit{\tau })+\varepsilon \rho \mathit{g},\end{array}$$(7)

where ɛ, u, ρ, p, f_s, τ, and g are, respectively, the void fraction, fluid velocity, fluid density, fluid pressure, solid–fluid 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

  

$$\mathit{U}=(1-\alpha )\hspace{0.17em}\mathit{u}+\alpha {\mathit{u}}_{\text{obj}}.$$(8)

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 volume-weighted average velocity U, Eq. (7) can be updated as

  

$$\begin{array}{l}\frac{\partial (\varepsilon \rho \mathit{U})}{\partial t}+\mathit{\text{∇}}·(\varepsilon \rho \mathit{UU})\\ =-\varepsilon \mathit{\text{∇}}p+{\mathit{f}}_{\mathbf{s}}+\mathit{\text{∇}}·(\varepsilon \mathit{\tau })+\varepsilon \rho \mathit{g}+{\mathit{F}}_{\text{IB}},\end{array}$$(9)

where F_IB is the correction term in the IBM. The force term F_IB is given by

  

$${\mathit{F}}_{\text{IB}}=\frac{\alpha \rho \left({\mathit{u}}_{\text{obj}}-\mathit{U}\right)}{\mathrm{\Delta }t}.$$(10)

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

  

$${\rho }_{\text{scale}}=\{\begin{array}{ll}\gamma \rho \hfill & \hfill (\phi <0)\\ \rho \hfill & \hfill (otherwise)\end{array},$$(11)

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

  

$$\begin{array}{l}\mathit{\text{∇}}·\left(\frac{{\varepsilon }^{n+1}}{{\rho }_{\text{scale}}}\mathit{\text{∇}}{p}^{n+1}\right)\\ =\frac{1}{\mathrm{\Delta }t}\left(\frac{{\varepsilon }^{n+1}-{\varepsilon }^n}{\mathrm{\Delta }t}+\mathit{\text{∇}}·{\varepsilon }^{n+1}{\mathit{U}}^{*}\right),\end{array}$$(12)

with the density scaling factor being introduced. In Eq. (12), superscript n + 1 and n are updated iteration number and present iteration number.

2.3 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.

Fig. 2

Schematic diagram of a die-filling system (Yao et al., 2018) Copyright: (2018) Elsevier B.V.

The target powder was NONPAREIL-108^® (FREUND CORPORATION). The particle size and the particle density were, respectively, 200 μm and 1,379 kg/m³, 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.

Fig. 3

Wall boundary modeling by the SDF and IBM.

Fig. 4

Effect of air flow on the powder filling.

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.

2.4 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; Takabatake et al., 2018), 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).

Fig. 5

Three-dimensional fluidized bed system (Mori, Wu and Sakai, 2019) Copyright: (2019) Elsevier B.V.

Fig. 6

Wall boundary model in a gas-solid flow simulation.

Fig. 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.

Fig. 7

Validation test result in a fluidized bed (Mori, Wu and Sakai, 2019) Copyright: (2019) Elsevier B.V.

3. 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 (Wang et al., 2018), 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.

3.1 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 liquid state is described in the subsequent sections.

The governing equations for phase change simulation by the MPS method include a Lagrangian description of the continuity equation, the Navier–Stokes equation, and energy conservation equation,

  

$$\frac{D\rho }{Dt}+\rho \mathit{\text{∇}}·\mathit{u}=0,$$(13)

  

$$\frac{D\mathit{u}}{Dt}=-\frac{1}{\rho }\mathit{\text{∇}}p+\frac{{\mu }_f}{\rho }{\mathit{\text{∇}}}^2\mathit{u}+\mathit{g}+\frac{1}{\rho }\sigma \kappa \delta \mathit{n},$$(14)

and

  

$$\frac{Dh}{Dt}=\mathit{\text{∇}}·(\lambda \mathit{\text{∇}}T)+Q,$$(15)

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,

  

$$Q=-\mathit{\text{∇}}·\mathit{q},$$(16)

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

  

$${⟨\mathit{\text{∇}}\mathrm{\Phi }⟩}_{\text{i}}=\frac{D}{n^0}\sum _{j\ne i}\left[\frac{{\mathrm{\Phi }}_{\text{j}}-{\mathrm{\Phi }}_{\text{i}}}{{\left|{\mathit{r}}_{\text{j}}-{\mathit{r}}_{\text{i}}\right|}^2}\left({\mathit{r}}_{\text{j}}-{\mathit{r}}_{\text{i}}\right)w\left({\mathit{r}}_{\text{j}}-{\mathit{r}}_{\text{i}}\right)\right],$$(17)

and

  

$${⟨{\mathit{\text{∇}}}^2\mathrm{\Phi }⟩}_{\text{i}}=\frac{2D}{\mathrm{\Omega }{n}^0}\sum _{j\ne i}\left[\left({\mathrm{\Phi }}_{\text{j}}-{\mathrm{\Phi }}_{\text{i}}\right)w\left({\mathit{r}}_{\text{j}}-{\mathit{r}}_{\text{i}}\right)\right],$$(18)

where D, n⁰, 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,

  

$$h_{\text{s}}=\rho C_{\text{p}}{T}_{\text{m}},$$(19)

  

$$h_1=h_{\text{s}}+\rho L,$$(20)

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,

  

$$\frac{{\mathit{u}}^{*}-{\mathit{u}}^k}{\mathrm{\Delta }t}=\mathit{g}-\frac{1}{\rho }\mathit{\text{∇}}p,$$(21)

  

$$\frac{{\mathit{u}}^{k+1}-{\mathit{u}}^{*}}{\mathrm{\Delta }t}=\frac{{\mu }_{\text{f}}}{\rho }{\mathit{\text{∇}}}^2{\mathit{u}}^{k+1}.$$(22)

A detailed discussion on the MPS method is presented in the original papers (Koshizuka and Oka, 1996).

3.2 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²/s, 2500 kg/m³, 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⁻⁴, which implies that the computation could be accelerated up to 10⁴ 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.

Fig. 8

Highly viscous fluid flow in a glass melter (Sun et al., 2012) Copyright: (2012) Elsevier B.V.

3.3 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).

Fig. 9

Verification test system for the heat flux model in the MPS method (Takabatake et al., 2016) Copyright: (2016) Elsevier B.V.

As boundary conditions, heat flux was 25,000 W/m² 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³, 1.0 × 10⁻³ m²/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.

Fig. 10

Temperature distribution in the verification test (Takabatake et al., 2016) Copyright: (2016) Elsevier B.V.

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³, 1.0 × 10⁻³ m²/s, 100 J/kg K, 50 W/m K, 10,000 K, and 1.0 × 10⁴ 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² 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.

Fig. 11

Snapshots obtained from a phase change simulation (Takabatake et al., 2016) Copyright: (2016) Elsevier B.V.

4. 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.

References

•  Adam  A.,  Pavlidis  D.,  Percival  J.R.,  Salinas  P.,  Xie  Z.,  Fang  F.,  Pain  C.C.,  Muggeridge  A.H.,  Jackson  M.D., Higher-order conservative interpolation between control-volume meshes: Application to advection and multiphase flow problems with dynamic mesh adaptivity, Journal of Computational Physics., 321 (2016) 512–531. https://doi.org/10.1016/j.jcp.2016.05.058
•  Anderson  T.B.,  Jackson  R., Fluid mechanical description of fluidized beds. Equations of motion, Industrial and Engineering Chemistry Research, 6 (1967) 527–539. https://doi.org/10.1021/i160024a007
•  Baserinia  R.,  Sinka  I.C., Powder die filling under gravity and suction fill mechanisms, International Journal of Pharmaceutics, 563 (2019) 135–155. https://doi.org/10.1016/j.ijpharm.2019.01.068
•  Basinskas  G.,  Sakai  M., Numerical study of the mixing efficiency of a ribbon mixer using the discrete element method, Powder Technology, 287 (2016a) 380–394. https://doi.org/10.1016/j.powtec.2015.10.017
•  Basinskas  G.,  Sakai  M., Numerical study of the mixing efficiency of a batch mixer using the discrete element method, Powder Technology, 301 (2016b) 815–829. https://doi.org/10.1016/j.powtec.2016.07.017
•  Chan  E.L.,  Washino  K., Coarse grain model for DEM simulation of dense and dynamic particle flow with liquid bridge forces, Chemical Engineering Research and Design, 132 (2018) 1060–1069. https://doi.org/10.1016/j.cherd.2017.12.033
•  Chen  R.,  Oka  Y.,  Li  G.,  Matsuura  T., Numerical investigation on melt freezing behavior in a tube by MPS method, Nuclear Engineering and Design, 273 (2014) 440–448. https://doi.org/10.1016/j.nucengdes.2014.03.049
•  Cundall  P.A.,  Strack  O.D.L., A discrete numerical model for granular assemblies, Geotechnique, 29 (1979) 47–65. https://doi.org/10.1680/geot.1979.29.1.47
•  Dhenge  R.M.,  Washino  K.,  Cartwright  J.J.,  Hounslow  M.J.,  Salman  A.D., Twin screw granulation using conveying screws: Effects of viscosity of granulation liquids and flow of powders, Powder Technology, 238 (2013) 77–90. https://doi.org/10.1016/j.powtec.2012.05.045
•  Duan  G.,  Yamaji  A.,  Koshizuka  S., A novel multiphase MPS algorithm for modeling crust formation by highly viscous fluid for simulating corium spreading, Nuclear Engineering and Design, 343 (2019) 218–231. https://doi.org/10.1016/j.nucengdes.2019.01.005
•  Ergun  S., Fluid flow through packed columns. Chemical Engineering Progress Symposium Series, 48 (1952) 89–94.
•  Ge  W., Chang  Q.,  Li  C.,  Wang  J., Multiscale structures in particle – fluid systems : Characterization, modeling, and simulation, Chemical Engineering Science, 198 (2019) 198–223. https://doi.org/10.1016/j.ces.2018.12.037
•  Govender  N.,  Wilke  D.N.,  Wu  C.,  Tuzun  U.,  Kureck  H., A numerical investigation into the effect of angular particle shape on blast furnace burden topography and percolation using a GPU solved discrete element model, Chemical Engineering Science, 204 (2019) 9–26. https://doi.org/10.1016/j.ces.2019.03.077
•  Govender  N.,  Wilke  D.N.,  Wu  C.Y.,  Rajamani  R.,  Khinast  J.,  Glasser  B.J., Large-scale GPU based DEM modeling of mixing using irregularly shaped particles, Advanced Powder Technology, 29 (2018) 2476–2490. https://doi.org/10.1016/j.apt.2018.06.028
•  Hu  C.,  Luo  K.,  Wang  S.,  Sun  L.,  Fan  J., Influences of operating parameters on the fluidized bed coal gasification process : A coarse-grained CFD-DEM study, Chemical Engineering Science, 195 (2019) 693–706. https://doi.org/10.1016/j.ces.2018.10.015
•  Hu  R., Fang  F.,  Salinas  P.,  Pain  C.C.,  Domingo  N.D.S.,  Mark  O., Advances in water resources numerical simulation of floods from multiple sources using an adaptive anisotropic unstructured mesh method, Advances in Water Resources, 123 (2019) 173–188. https://doi.org/10.1016/j.advwatres.2018.11.011
•  Kafashan  J.,  Wiącek  J.,  Abd  N.,  Jieqing  R., Two-dimensional particle shapes modelling for DEM simulations in engineering : A review, Granular Matter, (2019) 21:80. https://doi.org/10.1007/s10035-019-0935-1
•  Kajishima  T.,  Takiguchi  S.,  Hamasaki  H.,  Miyake  Y., Turbulence structure of particle-laden flow in a vertical plane channel due to vortex shedding, JSME International Journal Series B, 44 (2001) 526–535. https://doi.org/10.1299/jsmeb.44.526
•  Kanaoka  C., Fine particle filtration technology using fiber as dust collection medium, KONA Powder and Particle Journal, 36 (2019) 88–113. https://doi.org/10.14356/kona.2019006
•  Klinzing  G.E., Historical review of pneumatic conveying, KONA Powder and Particle Journal, 35 (2018) 150–159. https://doi.org/10.14356/kona.2018010
•  Koshizuka  S.,  Oka  Y., Moving-particle semi-implicit method for fragmentation of incompressible fluid, Nuclear Science and Engineering, 123 (1996) 421–434. https://doi.org/10.13182/NSE96-A24205
•  Krok  A.,  García-tri  P.,  Peciar  M.,  Wu  C.-Y., Finite element analysis of thermomechanical behaviour of powders during tableting, Chemical Engineering Research and Design, 110 (2016) 141–151. https://doi.org/10.1016/j.cherd.2016.03.019
•  Lu  L.,  Xu  J.,  Ge  W.,  Gao  G.,  Jiang  Y.,  Zhao  M.,  Liu  X.,  Li  J., Computer virtual experiment on fluidized beds using a coarse-grained discrete particle method—EMMS-DPM, Chemical Engineering Science, 155 (2016) 314–337. https://doi.org/10.1016/j.ces.2016.08.013
•  Monaghan  J.J., Simulating free surface flows with SPH, Journal of Computational Physics, 110 (1994) 399–406. https://doi.org/10.1006/jcph.1994.1034
•  Mori  Y.,  Wu  C.-Y.,  Sakai  M., Validation study on a scaling law model of the DEM in industrial gas-solid flows, Powder Technology, 343 (2019) 101–112. https://doi.org/10.1016/j.powtec.2018.11.015
•  Nie  Z.,  Zhu  Y.,  Wang  X.,  Gong  J., Investigating the effects of Fourier-based particle shape on the shear behaviors of rockfill material via DEM, Granular Matter, 21 (2019) 22. https://doi.org/10.1007/s10035-019-0875-9
•  Osher  S.,  Fedkiw  R., Level Set Methods and Dynamic Implicit Surfaces, Springer, 2002. ISBN 978-0-387-22746-7.
•  Peskin  C.S., The immersed boundary method, Acta Numerica, 11 (2002) 479–517. https://doi.org/10.1017/S0962492902000077
•  Russell  M.A.,  Souto-iglesias  A.,  Zohdi  T.I., Numerical simulation of Laser Fusion Additive Manufacturing processes using the SPH method, Computer Methods in Applied Mechanics and Engineering, 341 (2018) 163–187. https://doi.org/10.1016/j.cma.2018.06.033
•  Sakai  M., How should the discrete element method be applied in industrial systems ?: A review, KONA Powder and Particle Journal, 33 (2016) 169–178. https://doi.org/10.14356/kona.2016023
•  Sakai  M.,  Abe  M.,  Shigeto  Y.,  Mizutani  S.,  Takahashi  H.,  Viré  A.,  Percival  J.R.,  Xiang  J.,  Pain  C.C., Verification and validation of a coarse grain model of the DEM in a bubbling fluidized bed, Chemical Engineering Journal, 244 (2014) 33–43. https://doi.org/10.1016/j.cej.2014.01.029
•  Sakai  M.,  Kazuya  T.,  Kotaro  T.,  Hatano  K.,  Minakuchi  A., Why do wet-particles adhere to a high-speed roll in a three-roll mill ? Physics of Fluids, 31 (2019) 033302. https://doi.org/10.1063/1.5085693
•  Sakai  M.,  Koshizuka  S., Large-scale discrete element modeling in pneumatic conveying, Chemical Engineering Science, 64 (2009) 533–539. https://doi.org/10.1016/j.ces.2008.10.003
•  Sakai  M.,  Shigeto  Y.,  Basinskas  G.,  Hosokawa  A.,  Fuji  M., Discrete element simulation for the evaluation of solid mixing in an industrial blender, Chemical Engineering Journal, 279 (2015) 821–839. https://doi.org/10.1016/j.cej.2015.04.130
•  Sakai  M.,  Shigeto  Y.,  Sun  X.,  Aoki  T.,  Saito  T.,  Xiong  J.,  Koshizuka  S., Lagrangian–Lagrangian modeling for a solid–liquid flow in a cylindrical tank. Chemical Engineering Journal, 200–202 (2012a) 663–672. https://doi.org/10.1016/j.cej.2012.06.080
•  Sakai  M.,  Takahashi  H.,  Pain  C.C.,  Latham  J.-P.,  Xiang  J., Study on a large-scale discrete element model for fine particles in a fluidized bed, Advanced Powder Technology, 23 (2012b) 673–681. https://doi.org/10.1016/j.apt.2011.08.006
•  Sakai  M.,  Yamada  Y.,  Shigeto  Y.,  Shibata  K.,  Kawasaki  V.M.,  Koshizuka  S., Large-scale discrete element modeling in a fluidized bed, International Journal for Numerical Methods in Fluids, 64 (2010) 1319–1335. https://doi.org/10.1016/j.ces.2008.10.003
•  Schmelzle  S.,  Nirschl  H., DEM simulations : Mixing of dry and wet granular material with different contact angles, Granular Matter, 20 (2018) 19. https://doi.org/10.1007/s10035-018-0792-3
•  Shibata  K.,  Koshizuka  S.,  Sakai  M.,  Tanizawa  K., Transparent boundary condition for simulating nonlinear water waves by a particle method, Ocean Engineering, 38 (2011) 1839–1848. https://doi.org/10.1016/j.oceaneng.2011.09.012
•  Shigeto  Y.,  Sakai  M., Arbitrary-shaped wall boundary modeling based on signed distance functions for granular flow simulations, Chemical Engineering Journal, 231 (2013) 464–476. https://doi.org/10.1016/j.cej.2013.07.073
•  Shigeto  Y.,  Sakai  M., Parallel computing of discrete element method on multi-core processors, Particuology, 9 (2011) 398–405. https://doi.org/10.1016/j.partic.2011.04.002
•  Shrestha  S.,  Kuang  S.,  Yu  A.,  Zhou  Z., Particle shape effect on bubble dynamics in central air jet pseudo-2D fluidized beds : A CFD-DEM study, Chemical Engineering Science, 201 (2019) 448–466. https://doi.org/10.1016/j.ces.2019.02.030
•  Steuben  J.C.,  Iliopoulos  A.P.,  Michopoulos  J.G., Discrete element modeling of particle-based additive manufacturing processes, Computer Methods in Applied Mechanics and Engineering, 305 (2016) 537–561. https://doi.org/10.1016/j.cma.2016.02.023
•  Stroh  A.,  Daikeler  A.,  Nikku  M.,  May  J.,  Alobaid  F.,  Bohnstein  M. Von,  Ströhle  J.,  Epple  B., Coarse grain 3D CFD-DEM simulation and validation with capacitance probe measurements in a circulating fluidized bed, Chemical Engineering Science, 196 (2019) 37–53. https://doi.org/10.1016/j.ces.2018.11.052
•  Sun  X.,  Sakai  M., A liquid bridge model for spherical particles applicable to asymmetric configurations, Chemical Engineering Science, 182 (2018) 28–43. https://doi.org/10.1016/j.ces.2018.02.034
•  Sun  X.,  Sakai  M., Immersed boundary method with artificial density in pressure equation for modeling flows confined by wall boundaries, Journal of Chemical Engineering of Japan, 50 (2017) 161–169. https://doi.org/10.1252/jcej.16we115
•  Sun  X.,  Sakai  M., Direct numerical simulation of gas-solid-liquid flows with capillary effects: An application to liquid bridge forces between spherical particles, Physical Review E, 94 (2016) 063301. https://doi.org/10.1103/PhysRevE.94.063301
•  Sun  Xiaosong,  Sakai  M., Numerical simulation of two-phase flows in complex geometries by using the volume-of-fluid /immersed-boundary method, Chemical Engineering Science, 139 (2016) 221–240. https://doi.org/10.1016/j.ces.2015.09.031
•  Sun  X.,  Sakai  M.,  Sakai  M.-T.,  Yamada  Y., A Lagrangian–Lagrangian coupled method for three-dimensional solid– liquid flows involving free surfaces in a rotating cylindrical tank, Chemical Engineering Journal, 246 (2014) 122–141. https://doi.org/10.1016/j.cej.2014.02.049
•  Sun  X.,  Sakai  M.,  Shibata  K.,  Tochigi  Y.,  Fujiwara  H., Numerical modeling on the discharged fluid flow from a glass melter by a Lagrangian approach, Nuclear Engineering and Design, 248 (2012) 14–21. https://doi.org/10.1016/j.nucengdes.2012.04.004
•  Sun  X.,  Sakai  M.,  Sun  M.,  Takabatake  K.,  Pain  C.C., Numerical simulation of free surface fluid flows through porous media by using the explicit MPS method, Transport in Porous Media, 127 (2019) 7–33. https://doi.org/10.1007/s11242-018-1178-z
•  Takabatake  K.,  Mori  Y.,  Khinast  J.G.,  Sakai  M., Numerical investigation of a coarse-grain discrete element method in solid mixing in a spouted bed, Chemical Engineering Journal, 346 (2018) 416–426. https://doi.org/10.1016/j.cej.2018.04.015
•  Takabatake  K.,  Sun  X.,  Sakai  M.,  Pavlidis  D.,  Xiang  J.,  Pain  C.C., Numerical study on a heat transfer model in a Lagrangian fluid dynamics simulation, International Journal of Heat and Mass Transfer, 103 (2016) 635–645. https://doi.org/10.1016/j.ijheatmasstransfer.2016.07.073
•  Tsuji  Y.,  Kawaguchi  T.,  Tanaka  T., Discrete particle simulation of two-dimensional fluidized bed, Powder Technology, 77 (1993) 79–87. https://doi.org/10.1016/0032-5910(93)85010-7
•  Tsunazawa  Y.,  Fujihashi  D.,  Fukui  S.,  Sakai  M.,  Tokoro  C., Contact force model including the liquid-bridge force for wet-particle simulation using the discrete element method, Advanced Powder Technology, 27 (2016) 652–660. https://doi.org/10.1016/j.apt.2016.02.021
•  Tsunazawa  Y.,  Shigeto  Y.,  Tokoro  C.,  Sakai  M., Numerical simulation of industrial die filling using the discrete element method, Chemical Engineering Science, 138 (2015) 791–809. https://doi.org/10.1016/j.ces.2015.09.014
•  Udono  H.,  Sakai  M., A numerical study on dynamic inertial focusing of microparticles in a confined flow, Granular Matter, 19 (2017) 79. https://doi.org/10.1007/s10035-017-0758-x
•  Wang  X.,  Shao  M.,  Gao  S.,  Gau  J.,  Tang  H.,  Jin  H.,  Liu  H., Numerical simulation of laser impact spot welding, Journal of Manufacturing Processes, 35 (2018) 396–406. https://doi.org/10.1016/j.jmapro.2018.08.028
•  Wen  C.Y.,  Yu  Y.H., Mechanics of fluidization, Chem. Eng. Progr. Symposium Series, 62 (1966) 100–111.
•  Xiong  J.,  Koshizuka  S.,  Sakai  M., Investigation of droplet impingement onto wet walls based on simulation using particle method, Journal of Nuclear Science and Technology, 48 (2011) 145–153. https://doi.org/10.1080/18811248.2011.9711689
•  Xiong  J.,  Koshizuka  S.,  Sakai  M., Numerical analysis of droplet impingement using the moving particle semi-implicit method, Journal of Nuclear Science and Technology, 47 (2010) 314–321. https://doi.org/10.3327/jnst.47.314
•  Yamada  Y.,  Sakai  M., Lagrangian–Lagrangian simulations of solid–liquid flows in a bead mill, Powder Technology, 239 (2013) 105–114. https://doi.org/10.1016/j.powtec.2013.01.030
•  Yang  J.,  Guo  Y.,  Buettner  K.E.,  Curtis  J.S., DEM investigation of shear flows of binary mixtures of non-spherical particles, Chemical Engineering Science, 202 (2019) 383–391. https://doi.org/10.1016/j.ces.2019.03.027
•  Yao  H.,  Mori  Y.,  Takabatake  K.,  Sun  X.,  Sakai  M., Numerical investigation on the influence of air flow in a die filling process, Journal of the Taiwan Institute of Chemical Engineers, 90 (2018) 9–17. https://doi.org/10.1016/j.jtice.2017.11.031
•  Yokoi  K., Numerical method for interaction between multiparticle and complex structures, Physical Review E, 72 (2005) 046713. https://doi.org/10.1103/PhysRevE.72.046713
•  Yue  Y.,  Wang  T.,  Sakai  M.,  Shen  Y., Particle-scale study of spout deflection in a flat-bottomed spout fluidized bed, Chemical Engineering Science, 205 (2019) 121–133. https://doi.org/10.1016/j.ces.2019.04.031
•  Zakhvatayeva  A.,  Zhong  W.,  Makroo  H.A.,  Hare  C.,  Wu  C.Y., An experimental study of die filling of pharmaceutical powders using a rotary die filling system, International Journal of Pharmaceutics, 553 (2018) 84–96. https://doi.org/10.1016/j.ijpharm.2018.09.067
•  Zhou  Z.Y.,  Pinson  D.,  Zou  R.P.,  Yu  A.B., Discrete particle simulation of gas fluidization of ellipsoidal particles, Chemical Engineering Science, 66 (2011) 6128–6145. https://doi.org/10.1016/j.ces.2011.08.041

Nomenclature

DEM

discrete element method

CFD

computational fluid dynamics

MPS

moving particle semi-implicit

SDF

signed distance function

IBM

immersed boundary method

CPU

central process unit

GPU

graphical processing unit

D

dimension number (−)

F_Cn

normal component of contact force (N)

F_Ct

tangential component of contact force (N)

F_IB

correction term in the IBM (N/m³)

f_s

solid–fluid momentum exchange term (N/m³)

g

gravitational acceleration (m/s²)

h_s

enthalpy of solid phase (J/m³)

h_l

enthalpy of liquid phase (J/m³)

k

spring constant (N/m)

L

latent heat (J/kg)

n

normal vector (−)

n_pw

normal vector in the particle-wall interaction (−)

n⁰

initial particle number density (−)

p

fluid pressure (Pa)

Q

heat source (W/m³)

r

distance between two fluid particles (m)

r_p

particle radius (m)

Δt

time step (s)

U

volume-weighted average velocity (m/s)

u_obj

object velocity (m/s)

v_r

relative velocity of the interacting particles (m/s)

α

local volume fraction of a object (−)

δ_pw

displacement in particle-wall interaction (m)

δ

delta function (−)

ɛ

void fraction (−)

φ

signed distance function (m)

γ

density scaling factor (kg/s)

η

damping coefficient (−)

κ

curvature (m⁻¹)

λ

heat conductivity (W/m K)

μ_f

viscosity (Pa s)

μ_s

friction coefficient of a solid particle (−)

Ω

Lapracian constant in the MPS method (−)

ρ

fluid density (kg/m³)

σ

coefficient of surface tension (N/m)

τ

viscous stress (N/m²)

Authors’ Short Biographies

Mikio Sakai

Dr. Mikio Sakai is currently Associate Professor in Resilience Engineering Research Center in The University of Tokyo. He earned his Ph.D. degree from The University of Tokyo in 2006. Then, he became Assistant Professor in 2007 and Associate Professor in 2008 in the same university. He has been Visiting Reader at Imperial College London since 2016. He extensively studies modeling of granular flows, multi-phase flows and the heat transfer, besides the parallel computation techniques. He is a world-leading professor in computational granular dynamics, and hence, has delivered lots of invited lectures in conferences. He holds important posts in powder technology community, such as Director of Society of Powder Technology of Japan and Head of Simulation & Modeling Division in Association of Powder Process Industry and Engineering, JAPAN. At present, he is also Associate Editor of Chemical Engineering Science and Editor of Granular Matter.

Yuki Mori

Mr. Yuki Mori is currently a freshman PhD student in Department of Nuclear Engineering and Management, School of Engineering, The University of Tokyo, and is JSPS Research Fellow (DC1). He joined Prof Sakai’s group since 2016 and won Dean’s Award (Annual Best Student Award) upon completion of his masteral course in The University of Tokyo.

Xiaosong Sun

Dr. Xiaosong Sun earned his Ph.D. degree from The University of Tokyo in 2016. Then, he became Project Researcher in Resilience Engineering Research Center, School of Engineering, The University of Tokyo. He joined Prof Sakai’s group since 2011 and won Dean’s Award (Annual Best Student Award) upon completion of his masteral and doctorate courses in The University of Tokyo.

Kazuya Takabatake

Mr. Kazuya Takabatake is currently a third-year PhD student in Department of Nuclear Engineering and Management, School of Engineering, The University of Tokyo, and is JSPS Research Fellow (DC2) since 2018. He joined Prof Sakai’s group since 2015 and won Dean’s Award (Annual Best Student Award) upon completion of his masteral course in The University of Tokyo.