Development of Direct Analysis Model of Gas-Liquid-Solid Three Phase Flow
Using Particle Simulation
粒子シミュレーションによる固気液三相の直接流れ解析モデルの開発
Shungo Natsui
Assistant Professor
Division of Material Science, Faculty of Engineering,
Hokkaido University

1. Introduction
Control of multiphase flow in high temperature melt is ubiquitous in several material processes.
In the steel making or copper smelting processes, gas flow which passes between melting metal andslag affects reaction rate and refinement efficiency because of the change of gas-liquid interfacial
area (Bird et al. 2006; Nakamura et al. 1988; Sano et al. 1976; Pehlke et al. 1974). Metal electrolysis
process is also often accompanied with gas-liquid two phase flow such as gas generation around
anode or reduction of gas at cathode in molten salt (Otake et al.; 2013). In ironmaking process,because the gas flow permeability in the melting zone differs remarkably from that in the lumpy
zone, the control of this region has a large influence on gas flow. This means accurate understandingof the characteristics of the melt behavior is extremely important for constructing a low carbon blast
furnace (Hayashi et al., 2014).
In order to analyze multiphase behavior with molten metal or slag at high temperature, in-situ
observation by high temperature X-rays penetration equipment has been performed (Nakamura;
1988, Han; 2003). Visualization of melt and surface measurement are possible by using this method.
However, it is difficult to carry out the X-ray penetration of the inside of molten metal, and since theobservation domain is extremely limited, three-dimensional bubbles behavior cannot be known in
two dimensions. Because a flow and mass transfer of gas-liquid (with solid) interface are greatly
dependent on the physical properties of melt and very complicated, it is difficult to un-experiential
predict these phenomena and optimization of an operational condition.
From the background of progress of the latest computational science and technology, a highaccuracy numerical simulation raises expectation for solution approach to this complicated problem.
The authors have performed the numerical simulation of gas or solid or liquid flow in packed bed
process. (Natsui, 2011, 2012a, 2012b, 2012c, 2013). By using large-scale numerical simulation, we
can understand the unsteady three-dimensional information which is difficult to observe by
experimental method.
Now, in gas-liquid two phase flow, although many computational models have been proposed,－ 35 －

there are still many problems. Since an interface shape changes with time while a physical properties
changes discontinuously in a gas-liquid interface, it is difficult to analyze from microscopic dynamicbalance conditions. Thus the macroscopic two-fluid model which does not calculate gas-liquidinterface shape is still used widely practical. Two-fluid model requires to apply an empirical formulafor the momentum exchange between gas and liquid in compensation for skip of gas-liquid interface
shape calculation, and the calculation accuracy will be greatly dependent on the reliability of theexperimental data.
On the other hand, there is the interface tracking method as un-experiential approach which cancatch gas-liquid interface form with directly. The finite element method (Bonnerot; 1977, Lynch;1982) and a boundary-fitted coordinates system (Ryskin; 1984a, b) are the methods of making acalculation cell following interface modification, and need reconstruction of a calculation lattice.
The Volume of Fluid (VOF) method (Hirt; 1981) which pursues the share of fluid by the Euler
method is known widely. In VOF method, in the case of a complicated interface and advection of
bubbles, in order to prevent accumulation of the error of the interface form by coarse graining, andwhat is called numerical diffusion, a special technic is needed. These are, for example,
Donor-Acccptor method (Hirt, 1981), Level-set method (Osher, 1988), CIP method (Yabe, 1985),
PLIC(Kothe, 1998), MARS method (Kunugi, 1997), and CICSAM (Ubbink and Issa, 1999). Thesemethods have merit and demerit respectively, about reconstruction of calculation cells, volumeconservation, and the interface reconstruction method, etc., thus it is necessary to choose the optimal
method according to calculation condition.
Above-mentioned problems are caused by discretization of cells. Then, particle method is a
Lagrangian approach which does not use a mesh, and discretizes fluid by moving particles, andnumerical diffusion does not exist in essence. In particle method, particles themselves move byadvection with mass and volume. Therefore, an algorithm is simple and each particles keep interface
sharp. The Smoothed Particle Hydrodynamics (SPH) method (Lucy, 1977) is one of particle method
for fluid dynamics. It was originally developed for astrophysical applications and was later expanded
for applications in solid and fluid mechanics (Monaghan, 1994). It works by dividing the fluid into
asset of discrete elements, referred as particles and is based on integral representation of quantities
and spatial derivatives.

SPH method is suitable for the compressive fluid on condition ofnoticeable density change. On the other hand, Koshizuka (1996) proposed the Moving Particle
Semi-implicit (MPS) method which discretized incompressible Navier-Stokes equation by Taylor
series.
In particle method, gas phase movement was ignored and only liquid phase has been tracked likeas single-phase flow. This is in order to avoid the numerical instability by large pressure gradient due
to the large density difference. However, new algorithms for solving multi-density fluid directly is
proposed in recent years. Das and Das (2011) proposed Diffuse Interface (DI) model for SPH－ 36 －

method to solve the high density ratio fluid, and it is also performing to explore physical parametersof such bubbly flows (Das and Das, 2013). Shakibaenia and Jin (2012) proposed Weakly
Compressible MPS (WC-MPS) method, and applied to the Rayleigh-Taylor instability problem.
They compared with VOF method, and have concluded the WC-MPS method have a more sharp
interface. Grenier (2013) and Szewc (2013) applied the Multiphase-SPH model to more practicalbubbly flow. Although the density ratio of their intended phenomena was still not large such as
water-oil system, it was indicated to be able to apply this method to large scale industrial process.
It is required to numerically stabilize the gas-liquid interface more in order to apply tometallurgical process such as molten metal-slag-gas system. And solids form a packed bed, atwo-phase solid-gas system exists in the upper furnace, and a three-phase solid-gas-liquid system
exists in the lower furnace. If the stable gas-liquid model by particle simulation is developed,
seamless coupling with Discrete Element Method (DEM) which is one of solid particle simulation ispossible. In the present research, stability of gas-liquid interface was improved, and we verified thecomputational accuracy of this model based on a comparison of the numerical results withanalytical/experimental results.
2. Method
In the particle method, fluid flow is tracked by discretized particles in Lagrangian system. All
fluids (liquid and gas) are assumed to be incompressible fluid because of the sufficiently-small Mach
number. The fluids are spatially discretized by particles which have each character of density �.

However, in a gas/liquid two–phase flow, because the density at phase boundaries is not a steady, �is newly defined as the density function of each particle i.

(1)�� = 〈�〉�

where 〈 〉 represents the interparticle interaction, i.e. particle i is influenced by surroundings. Thismeans that the density of the fluid varies with the spatial distribution. The large density difference at
a two–phase interface dramatically affects the momentum exchange. To smooth the pressure gradient
between the gas and the liquid, a smoothed interface density function is introduced in boundary
particles.
〈�〉� = � �� ��� �� ��� ��
���

�

(2)

where �� is the mass of the particle j, � is the weight function, and � is the volume of the fluid.

The basic idea of the MPS method is to express incompressible fluid motion by keeping eachparticle ambient density constant. The sum of the weight functions of neighboring particles isdefined as the particle number density � (Koshizuka and Oka, 1996).

－ 37 －

��− 1 ���� � �� �〈�〉� = � ��� � ��� = ����� ���� � �� �
���

(3)

where ��� is distance between particle i and j �= ��� − �� ��, and �� is the influence radius of theweight function. In this study, �� = ��1�� is assumed (Tanaka and Masunaga, 2010). To maintainthe incompressibility of the fluid, a standard particle number density �� is given from the initialarrangement of uniform fluid particles.

�� = � ��� �

(4)

���

All the particles try to maintain a constant density and be close to �� . Therefore, numericalstabilization of pressure distribution is affected by �� , then the initial arrangements of particles will

be important. A simple cubic (SC) arrangement is commonly used, then the filling fraction for thisarrangement is relatively low (� � ����). In the case of the gas/liquid flow, particles move easily
under a pressure gradient (i.e. density difference ratio is approximately 1:1000 between gas andliquid). In the previous report, we introduced the checkerboard (CB) arrangement which is a highfilling fraction, pressure fluctuation was partially decreased. As same this idea, we employedhexagonal close-packed (HCP) in order to calculation more stably. In the HCP arrangement, a denser
state is assumed as the standard for a non–compressed state, and the numerical stability at a gas–liquid interface can be improved.
Flow of each fluid phase is governed by the equation of continuity and incompressible
Navier–Stokes equation.
��+ �� � � = �
��

��
1
1
= − �� + �� � � + � + ��
��
�
�

(5)
(6)

where � is the velocity, � is the time, � is the density, � is the pressure, � is the kinetic
viscosity, � is the gravity, and �� is the surface tension translated into a force per unit fluid
volume. Equation (6) is discretized by MPS–based method, i.e. pressure gradient term is calculated
by implicit method, other terms are calculated by explicit method.

����
�
��
1
1
= �− ���
+ ���� ��� + ���� + � �� �
��
�
�

(7)

An explicit solution to the particle motion equation is obtained by substituting the differential
operators in the governing equations into the MPS discrete model (Tanaka and Masunaga, 2012).
�1�∗ = �� + �� ��� � � + � + �� ��(8)where �∗ is velocity in the prediction step. The pressure term at the next step remains unknown－ 38 －

value. The pressure is solved from the particle number density in the prediction step 〈�∗ 〉� , which is

explicitly solved using the particle velocity and initial particle density. The gas phase as well as theliquid phase is discretized by particles in this model, and the numerical destabilization is often
increased owing to the very different density at the gas/liquid interface. This is caused by keeping
strict compressibility at the interface, thus the limited compressibility is assumed to stabilize thepressure at the interface. Because a smooth density function 〈�〉� is specified at the gas/liquidinterface, it should not be given as a constant. In this model, the spatial differentiation of 〈�〉� is

also considered. The following Poisson’s equation is employed.

〈�∗ 〉�
1γ� � �∗
���������� � �� � �� � � � 1�
〈�〉�
��
��
�

(9)

This formulization is also similar to a stable algorithm for calculating a compressible flow by an
Eulerian approach, and a convergent calculation is applied instead of an explicit calculation, thusinterfacial stability can be obtained. We employed γ � 0�01 considering the numerical stability andvolume conservation (Natsui et al., 2014). Therefore, each particle velocity in next step can be given
using obtained ������ as:

�

����

����
��
� � � �� ��� �
�
∗

(10)

Particle-particle interaction forces �� are added to the momentum conservation equation

to generate surface tension. The continuum surface force (CSF) approach model was localized to a

fluid interface by applying them to fluid elements in the transition region of the interface (Morris,
2000). Surface tension is converted to force per unit volume by:1 ��� � �� � ����� � �� ���� � �� � �
0 ��� � �� � �(11)where � is a delta function for judgment of surface

particles, the neighboring same phase particle number ��

is used. Here, � and � � are the standard neighboring

particle number in an incompressible state and acoefficient for surface judgment, respectively. In thepresent research, � � 0��� was assumed (Tanaka andMasunaga, 2010.). �� is the force per unit area given by

the interparticle potential force model. According to theprevious research, the exact form of this force is not very
important as long as it is anti-symmetric and short-range

Fig. 1. Schematic diagram of 3-phaseinterface in particle base discretization
model.

repulsive, long-range attractive (Kordilla et al. 2013). In this study, we construct �� following

Kondo et al. (2007) who employed this type of interaction force as same as Dissipative Particle
Dynamics (DPD) models.

－ 39 －

��� �� � ����� � �� ����� ��

��

���

� ���� � ������ � �� �

(12)

where � is the surface tension coefficient, and � is the interparticle potential. Natsui et al. (2012)verified that this form can estimate the absolute value of surface tension. In the case of multi-phaseflow, pairwise different phase will be formed with interface particles as shown in Figure 1. In thesurface, two interface particles should maintain a smooth pairwise relation. Here, a different liquid

phase A and B are defined, and let us consider a liquid A particle which are contact with gaseous
particles. We assume that surface tension �� is changed smoothly according to the particle number

ratio between liquid A and B. It is described as:�� � 〈�〉�

then,〈�〉� �

(13)��� �� � ��� ����� � ���(14)

where ��� is number of particle ��� �� around particle i which contact with gaseous particles, ���is number of particle ��� �� around particle i which contact with gaseous particles, �� is surface
tension coefficient of liquid A, �� is surface tension coefficient of liquid B, respectively. Thanks tothis simple assumption, surface tension working at A-gas and B-gas of two different interfaceschange smoothly. We employed this “mutual interface” model for gas/liquid/liquid simulation.
Liquid/liquid interface tension coefficient ��� remains as unknown value. It is derived by applying
the “Antonoff’s rule” (1907) described as:

3. Results and discussion��� � |�� � �� |(15)

In this section, we discuss the calculation verification in the bubble shape in a single liquid
phase. The time change in the maximum liquid heightand diameter of a gas bubble are shown in Figure 2.
Distilled water was employed in this experiment. The
rising of the bubble as the result of the buoyancy
caused by the density difference between the gas andthe liquid was seen. At the same time, a bubble has
become a horizontally long shape. The bubble behavior
was in good agreement with the experimental results.
This reason may be because initial HCP particle
arrangement and mutual interface model stabilized the
bubble behavior. The influence of liquid viscosity on－ 40 －

Fig. 2. Comparison of experimental andcalculated results for the change overtime in the maximum height, width.

bubble behavior was investigated. Figure 3 shows snapshots of the bubble behavior obtained by thesimulation at various times. Two types of silicon oil, 10cs and 100cs, were employed. Although thedata show a two–dimensional picture of transparent object, it can be seen that a dome–shaped bubble
rose in the liquid, changed its interfacial shapes. The difference of liquid viscosity changed bubble
shape; i.e. high viscosity shows the vertically long shape.
Next we discuss application of this model to gas/liquid/liquid system. Figure 4 shows thesnapshots of the bubble rising through the water/oil interface obtained by the simulation and theexperiment. At t = 0.20 s, a rising bubble modified water/oil interface, after that it pulls the waterphase upwards after the passage through the oil/water interface. The height of the water column
grows up and the stretched water column makes a neck around z = 0.09 m. The bubble height shows
good agreement with experimental data, thus interface shape change by convection might becorrectly estimated. On the other hand, the water column height obtained by calculation increased
later than the experiment. Observed liquid surface was little bit thicker than simulated result at t =0.20 s in Figure 4. Rapid rupture of “thin water film” would derive from a finite interface thickness
in this simulation model. The error resulting from thin film is remarkable after interface invasion of
air bubbles. In the experimental data, thin water film are observed in front and back side at t = 0.35
and 0.50 s in Figure 4. This is because the formation of a thin film on the liquid surface isdisregarded in this calculation condition, or liquid/liquid interface tension might not be able topredict enough accurately.
The behavior of even a liquid including simple single bubble changes to a fully
3-dimensional, so we believe that just our 2-dimensional information in the cold model experimental
picture is insufficient to predict the behavior in actual gas/liquid/liquid behavior. The main novelty

Figure 3. Snapshots of a gas bubble rising thorough oil.

－ 41 －

Figure 4. Snapshots of a single bubble
rising thorough a water/oil interface.

of this article is direct prediction of the liquid/liquid surface area �. Here, � can be predicted using
surface area at initial condition �� and surface delta function as:�� �� � � � �� ��� � �� �� �
�
�

(16)

�

Figure 5 shows the time change of simulated liquid/liquid surface area. � was increased withfluctuation until rupturing a thin water film (Fig. 6a). Although � was temporary decreased after
rupturing a thin film, it was significantly increased with water column growing. From this result, it isderived that even a single bubble stirs the liquid/liquid interface. When a bubble reached to surfaceof oil at t = 0.45 s, the water column stopped growth and � indicates the maximum at the point ofFig. 12b. Then coalescences of water droplet in oil decreased surface area. In the stage of
sedimentation, droplets slowly drop and it maintains �⁄�� ~1.20. A bubble affects an interface area,even after passing through the water/oil interface.

Figure 5. Change over time in the surface areaof a water/oil interface; a: breaking of thin film,
b: breaking of bubble.

4. Conclusion
We newly developed a particle method computational fluid dynamics model for molten
metal-slag-gas multi-phase flow, and have investigated the model accuracy based on comparing withcold model experiments. In this model, the gas phase and liquid phase are directly discretized as－ 42 －

particles, therefore it is practical and useful for application to actual metallurgical process.
Numerical stability was improved by using denser initial particle arrangement. The surface tension
model based on inter-particle potential expanded to immiscible three fluids flow by considering
multi-interface force balance.
For the first step, we applied this model to simple liquid/gas system, and compared the results ofsimulation with experimental photographs. Interface modification was recorded using high-speed
video camera, and experimental data was analyzed. The time change of bubble height and shapewere quantitatively in good agreement with the experimental results, thus this model successfully
predicted bubble rising behavior in each viscous fluid.
Next, we investigated an entrainment phenomenon of the water column pulled by a risingbubble through a water/oil interface. The time change of interface shape by convection might becorrectly estimated. However, rupture of thin water film which covered around bubble has not beennecessarily predicted. This is because the liquid/liquid interface tension might not be able to predictenough accurately. Nevertheless, this 3-dimensional unsteady numerical model gives new
information to predict the time change of liquid/liquid interface area based on fluid dynamics.
Prediction of such a continuous interface change has not been obtained by the conventional static
model. Another advantage of this technique is to be able to apply to further complicated system suchas actual metallurgical systems. In the future, further developing of interface accuracy will be aneffective method for predicting the surface area by simulating the position, velocity, pressure, andother parameters of the fluid elemental volume.
Now we are also trying to develop DEM-based
block model for seamless analysis of gas-liquid-solid3-phase flow. The behavior of an individual block is
usually tracked using a rigid body model as shown in
Figure 6. DEM and MPS have common characteristics,
such as being Lagrangian models, the passively
moving solid model shows high affinity with DEM. In
a present stage, this direct 3-phase flow model
indicates poor convergency due to the inter-particle
potential force of surface tension. It needs to develop
the further numerical stability especially at the 3-phase

Figure 6. Sample calculation of rigid
body DEM model. Color scale shows thevelocity of each element.

interfaces.
Acknowledgements
The financial support by the JFE 21st Century Foundation is gratefully acknowledged.

－ 43 －

References
Bird, R. B., W. E. Stewart and E. N. Lightfoot (2006): “Transport Phenomena”, 2nd ed., John Wiley& Sons, New York, 554-562.
Bonnerot, R., and P. Jamet (1977): Numerical computation of the free boundary for the
tow-dimensional Stefan problem by Space-time Finite Element, J. Comput. Phys., 25, 163-181.
Das, A. K. and P. K. Das (2011), Int. J. Numer. Meth. Fluids, 67, 671–699.
Das, A. K. and P. K. Das (2013), “Bubble Evolution and Necking at a Submerged Orice for theComplete Range of Orice Tilt”, AIChE J, 59, 630–642.
Gotoh, H. and T. Sakai (2006), “Key issues in the particle method for computation of wave
breaking”, Coastal Engineering, 53, 171-179.
Grace, J. R.: “Shapes and Velocities of Bubbles Rising in Infinite Liquids” Trans. Inst. Chem. Eng.,
51 (1973),116-120.
Grenier, N., D. Le Touzé, A. Colagrossi, M. Antuono, and G. Colicchio: Viscousbubbly flows
simulation with an interface SPH model, Ocean Eng., 69(2013), 88-102.
Han, Z., and L. Holappa (2003), Mechanism of iron entrainment into slag due to rising gas bubbles,
ISIJ International, 43, 292-297.
M. Hayashi, S. Sukenaga, K. Ohno, S. Ueda, K. Sunahara, and N. Saito, Tetsu-to-Hagané, 100(2014),
211.
Hirt, C. W., and B. D. Nichols (1981): Volume of Fluid(VOF) method for dynamics of freeboundaries, J. Comput. Phys., 39, 201-225.
Hestenes, M. R. and E. Stiefel (1952), “Methods of conjugate gradients for solving linear systems”,
Journal of Research of the National Bureau of Standards, 49, 409–436.
Ikeno, T., and T. Kajishima: “Finite-difference immersed boundary method consistent with wallconditions for incompressible turbulent flow simulations”, J. Comput. Phys., 226 (2007), 1485-1508.
Khayyer, A. and H. Gotoh (2013), “Enhancement of performancw and stability of MPS mesh-freeparticle method for multiphase flows characterized by high density ratio”, Journal of Computational
Physics, 242, 211-233.
Kon T., S. Natsui, S. Ueda, R. Inoue, and T. Ariyama (2012), Numerical Simulation of Dripping
Behavior of Droplet in Packed Bed Using Particle Method, ISIJ International, 52, 1565-1573.

－ 44 －

Kon, T., S. Natsui, S. Ueda, R. Inoue, and T. Ariyama (2013), Influence of Physical Properties ofMelt on Liquid Dripping in Packed Bed Analyzed by MPS Method, ISIJ International, 53, 590-597.
Koshizuka, S. and Y. Oka (1996): Moving-particle semi-implicit method for fragmentation ofincompressible fluid, Nuclear Science and Engineering, 123, 421-434.
Kunugi, T., (1997), “Direct Numerical Algorithm for Multiphase Flow with Free Surface andInterface”, Trans. Japan Soc. Mech. Eng. Ser. B, 63, 1576-1584.
Lucy, L. B. (1977), A numerical approach to the testing of the fission hypothesis, Astronomical
Journal, 82, 1013-1024.
Lynch, D. R. (1982): Unified approach to simulation on deforming elements with application tophase change problems, J. Comput. Phys., 47, 387–411.
Monaghan, J. J. (1989): On the Problem of Penetration in Particle Methods, J. Comput. Phys., 82, 1–
15.
Monaghan, J. J. (1994): Simulating free surface flows with SPH, J. Comput. Phys., 110, 399–406.
Monaghan, J. J. (2012): Smoothed particle hydrodynamics and its diverse applications, Ann. Rev.
Fluid Mech. 44, 323-346.
Nakamura, T., F. Noguchi, Y. Ueda, S. Nakajo (1988), Floating Behaviour of Molten Copper andCu-matte in Cu-slag Study on interfacial phenomena in phase separation of copper smelting (2ndReport), J. MMIJ, 104, 531-536.
Natsui, S., H. Nogami, S. Ueda, J. Kano, R. Inoue and T. Ariyama (2011): Simultaneous
three-dimensional analysis of gas–solid flow in blast furnace by combining discrete element method
and computational fluid dynamics, ISIJ International, 51, 41-50.
Natsui, S., S. Ueda, H. Nogami, J. Kano, R. Inoue and T. Ariyama (2012a), Gas-Solid Flow
Simulation of Fines Clogging a Packed Bed using DEM-CFD, Chemical Engineering Science, 71,
274-282.
Natsui S., R. Soda, T. Kon, S. Ueda, J. Kano, R. Inoue and T. Ariyama (2012b), Wettability Model
Considering Three-phase Interfacial Energetics in Particle Method, Materials Transactions, 53,
662-670.
Natsui S., T. Kon, S. Ueda, J. Kano, R. Inoue, T. Ariyama and H. Nogami (2012c), Analysis of Heatand Mass Transfer in a Packed Bed by Considering Particle Arrangement, Tetsu-to-Hagane,
98(2012), 341-351.

－ 45 －

Natsui S., R. Shibasaki, T. Kon, S. Ueda, R. Inoue and T. Ariyama: Effect of High Reactivity Cokefor Mixed Charge in Ore Layer on Reaction Behavior of Each Particle in Blast Furnace, ISIJ Int., 53(2013), 1770–1778.
Otake, K., H. Kinoshita, T. Kikuchi, and R. O. Suzuki (2013): CO2 gas decomposition to carbon byelectro-reduction in molten salt, Electrochimica Acta, 100, 293-299.
Pehlke, R. D., W. F. Porter, R. F. Urban (1974): BOF steelmaking, vol. 2, New York, USA, 287.
Rider, W. J. and Kothe, D. B. (1998): "Reconstructing Volume Tracking", J. of Comp. Phys., 141,112-152.
Ryskin, G. and L. G. Leal (1984): Numerical solution of free-boundary problems in fluid mechanics,
Part1. The Finite-difference technique, J. Fluid Mech., 148, 1-17.
Ryskin, G. and L. G. Leal (1984): Numerical solution of free-boundary problems in fluid mechanics,
Part2. Buoyancy-driven motion of a gas bubble through a quiescent liquid, J. Fluid Mech., 148,
19-35.
Sano, M., Y. Fujita and K. Mori, (1976): Metal. Trans. B, 7B, 300.
Shakibaenia, A., and Y. C. Jin, (2012): MPS mesh-free particle method for multiphase flows, 13,229-232.
Szewc, K., J. Pozorski, and J. P. Minier, (2013): Simulations of single bubbles rising through viscous
liquids using Smoothed Particle Hydrodynamics, Int. J. Multiphase Flow, 50(2013), 98.
Tanaka, M., and T. Masunaga (2010): Stabilization and smoothing of pressure in MPS method by
Quasi-Compressibility, Journal of Computational Physics, 229, 4279–4290.
Ubbink, O. and R. I. Issa, (1999): A Method for Capturing Sharp Fluid Interfaces on Arbitrary
Meshes, J. Comput. Phys. 153 , 26-50.
Takewaki, H., A. Nishiguchi and T. Yabe: The Cubic-Interpolated Pseudo-Particle(CIP) Method for
Solving Hyperbolic-Type Equations, J. Comput. Phys., 61 (1985), 261-268.
Yabe, T. and P. Y. Wang: Unified Numerical Procedure for Compressible and Incompressible Fluid, J.
Phys. Soc. Jpn., 60 (1991), 2105-2108.

－ 46 －