Research Article  Open Access
Development of Moving Particle Simulation Method for MultiliquidLayer Sloshing
Abstract
The mixed oil and gas including water and sand are extracted from well to offshore structure. This mixed fluid must be separated for subsequent processes by using wash tanks or separators. To design such a system, a proper numericalprediction tool for multiphase fluids is required. In this regard, a new moving particle simulation (MPS) method is developed to simulate multiliquidlayer sloshing problems. The new MPS method for multifluid system includes extra search methods for interface particles, boundary conditions for interfaces, buoyancycorrection model, and surfacetension model for interface particles. The new particle interaction models are verified through comparisons with published numerical and experimental data. In particular, the multiliquid MPS method is verified against Molin et al’s (2012) experiment with three liquid layers. In case of excitation frequency close to one of the internallayer resonances, the internal interface motions can be much greater than top freesurface motions. The verified multiliquid MPS program is subsequently used for more nonlinear cases including multichromatic multimodal motions with larger amplitudes, from which various nonlinear features, such as internal breaking and more particle detachment, can be observed. For the nonlinear case, the differences between with and without buoyancycorrection and surfacetension models are also demonstrated.
1. Introduction
As the demand of oil/gas increases rapidly, offshore oil/gas production is continuously extended to deeper waters and more floating offshore structures with storage space, such as FPSOs and SPARs, are to be manufactured and installed including multiwell system.
In order to increase the production/processing efficiency for floating platforms, many new ideas have been suggested. For instance, subsea separators are newly introduced for several deepwater wells. Another new idea is placing largescale “wash tanks” inside the hull as preprocessors for conventional separators. The “wash tank” is for initial separation of mixed fluids that consisted of oil and sea water produced through production risers. For maximum processing efficiency, an emulsion fluid can be added in the wash tank to help in initial separation. Since the size of wash tank is large, the innerfluid motions cannot be ignored in predicting platform motions. Therefore, the analysis of the innerfluid dynamics induced by platform motions and vice versa is important for the design of such vessels with liquid storage system.
The sloshing problem by a single liquid has been the subject of many studies (e.g., Chen [1]; Lee et al. [2, 3]; Loysel et al. [4]). Its interaction with vessel motions has also been extensively investigated (e.g., Kim et al. [5, 6]; Gaillarde et al. [7]; Cho et al. [8], Lee and Kim [9]; Kim et al. [10, 11]). However, sloshing problems by multilayerliquid system have rarely been attempted. Twoliquidlayer sloshing in a rectangular tank was theoretically and experimentally studied by La Rocca et al. [12]. Molin et al. [13] investigated threeliquidlayer sloshing and their numerical simulations were compared against their experiments. Both potential theory and VOF simulations were used in their numerical predictions. In this paper, a new MPS method is developed for multiliquidlayer system and the simulation results are compared with the threeliquidlayersloshing experiments by Molin et al. [13]. Since the developed MPS method is suitable for the simulations of very violent flows including breaking, mixing, slamming, splashing, fragmentation, and coalescence, it is also applied to highly nonlinear cases with large motion amplitude and irregularity.
The incompressible viscous Lagrangian CFD method, MPS, has been used by many researchers, such as Koshizuka and Oka [14], Toyota et al. [15], Gotoh [16], Lee et al. [2], and Hwang et al. [17], for many engineering applications. The original MPS of Koshizuka and Oka [14], which includes nonphysical pressure fluctuations and lessoptimal treatments in the Poisson equation, gradient/collision model, and tracing method for freesurface particles, has also been improved by many researchers (e.g., Lee et al. [3] and Gotoh [16]). However, the improved MPS models have been limited for singlephaseliquid problems. In order to simulate multiphasefluid motions with multiple interfaces, Shirakawa et al. [18] proposed a buoyancy model and Nomura et al. [19] suggested a surfacetension model to better simulate interfaces. Khayyer and Gotoh [20] suggested a method to solve the multiphase interfaces of high density ratio. There has been similar development based on the similar Lagrangian approach, called SPH (smoothed particle hydrodynamics) method (e.g., Monaghan and Kocharyan [21] and Hu and Adams [22]).
In the present study, a new MPS method for multiliquid layers including improved buoyancy model is developed (Kim et al. [11]). The buoyancy model has been devised based on relevant physics in a more rigorous form compared to the empirically devised formula of Shirakawa et al. [18] with arbitrary parameters. The inner interfaces are traced by the analogous but more improved method compared to that of finding freesurface particles. A surfacetension model is also developed for interface particles. A hydrostaticcorrection model is also suggested. The simulation results by the newly developed MPS program are compared with the experimental results by Molin et al. [13] for threelayerliquid sloshing. The case studies show that the newly developed program is very useful in simulating various engineering problems including multiple fluids and interfaces.
2. MPS Method for SingleLiquidLayer
In this section, the MPS method for singleliquidlayer is explained. More detailed information about improved MPS method can be found in Lee et al. [3] and Kim et al. [10]. The improvements compared to original MPS by Koshizuka and Oka [14] include optimal kernel function, improved gradient model and source term for Poisson equation, and improved method for freesurface searching and collision model.
2.1. Governing Equations
The MPS method is based on Lagrangian approach including the effect of particle interactions. To solve fluid dynamics by MPS, the continuity and NavierStokes equations are employed as governing equations for incompressible viscous fluid. Consider where is the density, is the particle velocity, is the pressure, is the kinematic viscosity, is the gradient, is the time, and is the external force.
The righthand side of NavierStokes equation, (2), consists of pressure gradient, viscous, and external force terms, respectively. To simulate the incompressible viscous flows, all terms of differential operators should be replaced with the particle interaction models, which include gradient, divergence, and Laplacian models. The particle interaction models are based on the kernel function which represents the effect of neighboring particles with respect to the distance from center particle. In this study, the following kernel function is employed: where is distance between particles and represents the effective range of particle interactions. If particle distance, , exceeds the effective range, , the kernel function becomes zero.
2.2. Gradient, Divergence, and Laplace Models
Gradient model represents a local weighted average of the gradient vectors between center particle and its neighboring particles . The gradient model can be expressed as follows: where denotes particle interaction model, is an arbitrary scalar function, is number of dimensions, is position vector of particle, and is the particle number density at initial arrangement. Substituting pressure instead of , the pressure gradient can be obtained as The particle number density, , can be calculated by following formula: Divergence model can also be obtained from The Laplacian model can be expressed as where The symbol in (8) is the parameter to make the increase of variance measured by distribution of particles equal to an increase of variance from the unsteady diffusion equation. The Laplacian model is used to calculate viscous effects. It is measured by viscosity of center particle, which means that the weighted/averaged viscosity of neighboring particles inside effective radius is used here.
2.3. Incompressibility Model
The algorithm of incompressibility model in MPS method is similar to SMAC (simplified markerandcell) method in grid based CFD system. It consists of two stages, explicit and implicit stages. In the explicit stage, the velocity and position are predicted explicitly. In this stage, particles move with the intermediate velocity, , calculated by viscous and gravitational forces. Since the continuity equation has to be satisfied with the particle number density, the velocity corrector, , which adjusts particle arrangement, needs to be introduced. The velocity corrector can be calculated by the pressure gradient and the pressure gradient can be obtained through PPE (Poisson pressure equation) implicitly. Tanaka and Masunaga [23] suggested the mixed source term of PPE to achieve less nonphysical pressure fluctuation. This mixed source term has been improved by Lee et al. [3] in the following form: where the asterisk denotes temporal quantities obtained by explicit stage, is the pressure at time , is the relaxation parameter, is the particle number density at initial arrangement, is the intermediate velocity, and is time interval. The improved formula (10) can be used with (5) in a very stable manner. In this study, flexible time interval is used, which means that the time interval can be adjusted by the following equation: where is maximum particle velocity, is particle distance at initial stage, and is Courant number and it is set to be 0.2 in this simulation. The program runs with the initial time interval if the calculated by (11) exceeds the initial time interval. Otherwise, the adjusted time interval is used. This way, the time interval is adjusted in every single cycle.
2.4. Boundary Conditions
The kinematic and dynamic boundary conditions are imposed on free surface. Due to tracing moving particles on the free surface, the kinematic boundary condition is automatically satisfied. The dynamic freesurface boundary condition is satisfied by taking the pressure on the freesurface particles to be the same as atmospheric pressure .
To apply the dynamic freesurface condition, the freesurface particles should be identified first. In the freesurface region, the particle density decreases since there are no particles in air region. Thus, the following simple conditions are used to identify the freesurface particles: where and are parameters below 1.0, is the particle number density of particle , is the number of neighboring particles within the effective range of particle interaction , and is the maximum number of neighboring particles fully submerged in the initial distribution. Especially, the freesurface parameters are used to judge whether the particles are on the free surface or not and and are used in this study. The choice is based on the extensive numerical experiments (Lee et al. [3]). Using this freesurface boundary condition, we can also simulate fragmentation and coalescence of freesurface flow.
2.5. Collision Model
When particles get close in the fluid, the local pressure causes repulsive force, which can be adjusted by introducing collision model. Particularly, since the pressure of freesurface particles is set to be zero, the repulsive force might be improperly generated when particles get close enough. Also, the sudden increase of particle number density can occur when particles from outside come into collision with freesurface particles. Since the sudden increase of particle number density influences identifying freesurface particles, the spatial stability of pressure can be reduced. In order to avoid the reduction of stability, a special treatment for collision model is necessary, especially near the free surface.
At the initial stage, particles are arranged uniformly with constant distance of . When the distance between two particles is less than , the collision model is applied. After applying the conservation of momentum, the repulsive velocity can be calculated by using the coefficient related to repulsiveness and defined as the ratio of particle velocities before and after collision. Here, and are set as 0.97 and 0.2, respectively, after extensive numerical optimalvalue search including both static and dynamic problems. Particularly, the value selected reduces nonphysical pressure fluctuations compared to other values lower than that (e.g., Lee et al. [3]). This collision model is applied to all fluid particles including freesurface particles.
2.6. Hydrostatic Pressure Correction Model
The hydrostatic pressure of the toplayer particles is lost due to imposing free surface boundary condition, as shown in Figure 1. In order to adjust this hydrostatic error, the HSC (hydrostatic pressure correction) model is suggested here.
(a)
(b)
The hydrostaticpressure correction is only applied to freesurface particles, but the detached/splashed/fragmented particles are excluded. By adding (13) to the freesurfaceparticlesearching method, the particles for HSC can be identified. Consider where is a proper number obtained by numerical test. In this study, 0.4 is used for . can be calculated by the following equation: When th particle is freesurface and HSCrequired particle, its center pressure is given by, which is equivalent to applying the dynamic freesurface condition at the true free surface (top of freesurface particle):
2.7. Validation of SinglePhase Case
In order to validate the developed singlephase MPS, the sloshing experiment conducted by Kishev et al. [24] is numerically reproduced. The model is shown in Figure 2. The tank is forced to move by horizontal sinusoidal motion with the period of 1.5 seconds. The simulated pressure is compared with the experiment result. As shown in Figure 3, both cases are in good agreement with each other. Not only magnitude but also specific trends are accurately reproduced by the present MPS.
3. MPS Method for MultiliquidLayer
3.1. Interface Searching Method and Interface Boundary Conditions
In the multiliquid problem, there is not only free surface but also interfaces between liquids. Since interface is similar to free surface, the concept of freesurface searching method, (12) and (13), can be used for interface. However, it is necessary to modify the freesurface searching method to be suitable for the interface tracing method.
The interface searching has one more condition compared to the freesurface searching, which is added to eliminate detached particles. According to (12), a particle can be identified as freesurface particle when its particle number density is less than . At interfaces, an additional condition is given for minimum value for interface particles. In this model, the new particle number density is introduced, which can be measured only by the same phase of fluid. Therefore, the interface particle can be identified through the following conditions: where The second subscript, , represents the phase of particle. In the present study, , , and are used after extensive numerical experiments. The difference between freesurface searching and interface searching is that detached particles are not included in applying dynamic interface boundary conditions that is pressure continuity at the interface in the normal direction. The kinematic interface condition is again automatically satisfied by using particle tracing. Figure 4 shows that interface particles are correctly traced by the present approach. By additionally introducing , the detached particles from interface are not recognized as interface particles, which is more robust in applying interface boundary condition at the actual interface.
(a)
(b)
After identifying the interface particles, the dynamic boundary condition at the interface can be expressed as follows: where is the gravitational acceleration and and are the vertical positions for particle and , respectively. Since the interface is a very thin layer and the particle method cannot formulate this very thin layer, the static pressure adjustment needs to be considered as the second term of lefthand side of (19) to improve accuracy.
3.2. BuoyancyCorrection Model
MPS uses particle number density to interact between particles. However, this method can underestimate the selfbuoyancy of center particle for multiliquid cases. Therefore, the correction of buoyancy effect is necessary for multilayerliquid problems.
Since buoyancy force is related to pressure gradient, the buoyancycorrection model can be implemented into the implicit stage of adjusting the velocity corrector.
In (5), the pressure gradient model, , and include both dynamic and hydrostatic pressure. However, the selfbuoyancy of the center particle due to surrounding fluid of different density is not included because the kernel function is omnispread from the center value with axial symmetry. In case of single fluid, the center particle is neutrally buoyant and the selfbuoyancy correction is not necessary.
In order to find the selfbuoyancy of center particle, let us consider the following two cases shown in Figure 5. In both cases, particles of different density are arranged vertically. In the first case, the lower particle is regarded as center particle. In the second case, the upper particle is regarded as center particle. In both cases, it is assumed that the center particle has lighter density than neighboring particles. The selfbuoyancy of the center particle can, for example, be obtained by subtracting the pressure of upper particles from that of lower particles. Therefore the selfbuoyancy can be expressed as follows: In this equation, the particle radius, , can be generalized to half of the vertical distance between and particles; that is, . The first terms of the righthand sides of (20) are already included in the pressure gradient model of MPS method, while the second terms of the righthand sides of (20) can be implemented into the pressure gradient model as buoyancycorrection model. The top and bottom particles are used for illustration and all surrounding particles contribute to the center particle buoyancy in a weighted manner. Therefore, the pressure gradient for multiphase MPS with buoyancycorrection model can be expressed as follow: where is the kernel function for buoyancycorrection model, which is the same as (3), but the effective range, , is 1.5. The symbols and are vertical positions of particles and , and is velocityadjustment parameter similar to that used in Shirakawa et al. [18]. Its default value is 0.5. The wall and dummy particles are not included in this correction. The present buoyancycorrection model is more straightforward and theoretically clear compared to that of Shirakawa et al. [18] that includes several arbitrary empirical parameters. In case of the singlephase fluid, the buoyancycorrection term vanishes.
In order to validate the developed buoyancycorrection model, a simulation case shown in Figure 6 is considered. The rectangular tank is filled with fluids which are water and oil. In this case, total 17020 particles are used including 12959 particles for water and 1 particle for oil. The initial particle distance is 0.005 m and time interval is 0.001 second. The Oil particle is located near the bottom of tank and has less density (950 kg/m^{3}) than surrounding water particles (1000 kg/m^{3}). The oil particle is to freely rise up to the free surface by selfbuoyancy effect. In order to avoid sidewall effect, the length of tank is set to be large enough as 1.08 m.
Because the buoyancy force is constantly applied to the oil particle, it will rise vertically with constant terminal velocity. However, the MPS method without buoyancycorrection model does not show that trend, as shown in Figure 7. Due to the underestimation of the selfbuoyancy effect, the rising velocity is extremely slow and irregular. By applying the present buoyancycorrection model, the oil particle rises with constant vertical velocity induced by the net selfbuoyancy effect. The necessity of the buoyancycorrection model can be proved through this example. The rising velocity can be adjusted by using the velocityadjustment parameter, as shown in the figure. As the parameter increases, the rising velocity is also increased. If the oilparticle ( m) is replaced by an idealized spherical solid particle of the same density, its rising terminal velocity ( where is density, is diameter of particle, is drag coefficient, and subscription and denote object and surrounding fluid, resp.) is estimated to be 0.05 m/s with the drag coefficient of 1.0, which is very similar to the present case with the default value of . The velocityadjustment parameter can be tuned if necessary depending on specific problems. The different size of oil droplet is also tested. When oil droplet diameter is increased twice ( m), the simulated rising velocities with various are also compared with the theoretical terminal velocity. The estimated terminal velocity 0.068 m/s with the default value is again the closest to the above formula. Therefore, it is seen that the developed buoyancy model is not sensitive to the droplet size.
(a) Particle diameter is 0.01 m
(b) Particle diameter is 0.005 m
3.3. Surface Tension Model
The surface tension is also present at the interface, so it needs to be applied to interface particles. In order to simulate surface tension for interface particles, two more particle number densities are newly adopted: where is the effective range for surface tension. It is set as 3.1, while the effective range for other models is 2.1.
The curvature, , can be calculated with the newly adopted particle number density for surface tension: where is particle number density of surface tension model which can be measured by (21). The normal vector at the interface can be calculated from the following equation: Then, the surface tension force, , is additionally applied to only the interface particles, where is surface tension coefficient (Figure 8). In order to validate the surfacetension model, vibration of an ethanol droplet in twodimensional domain without gravitational force is simulated, as shown in Figure 9. Total 900 particles are used for fluid particles. Particles are initially arranged to make a square with the side length of 75 mm and only surface tension is applied as external force. The density of fluid is 798 kg/m^{3}, the kinematic viscosity is m^{2}/s, and the surface tension is N/m. As time grows, the square shape becomes more like circle due to its surface tension force after several steps of vibrations, which can also be observed in both experiment and theory (Lamb, [25]). The vibration period also quantitatively agrees with the theoretical value of Lamb [25].
3.4. ThreePhaseSloshing Example
For further verification, the developed multiphaseliquid MPS program is run to compare with Molin et al.’s [13] threeliquidlayersloshing experiment. The detailed information of the threeliquidlayer experiment can be found in Molin et al. [13]. All the simulation conditions are set equal to the experiment. The rectangular tank is filled with three liquids whose properties are given in Table 1. The heaviest but less viscous liquid, dichloromethane, has 15 cm height from the bottom of tank. In the middle, water is located with the height of 15 cm from dichloromethane. At the top layer, cyclohexane is located with 38 cm height from water (see Figure 10). The air region above the cyclohexane is regarded as vacuum without particles. The resonant sloshing frequencies at the three interfaces are obtained by Molin et al. [13] with linear potential theory. Their natural frequencies are given in Table 2 at corresponding natural sloshing modes. The sloshing modes mean patterns of wave motions with respective wave lengths. For the present MPS simulation, a total of 33436 particles are used 6510 particles for bottom layer, 6510 for middle layer, and 16356 for top layer. 4060 particles are also used to represent solid tank wall. For this arrangement, the initial particle distance is 0.005 m. In order to test convergence with the number of particles, two more initial distances of 0.0075 m and 0.01 m are also considered. In case of 0.0075 m, the total number of particles, 15700, consists of 2740 wall particles, 2859 particles for middle and bottom layers, and 7242 particles for top layer. The initial time interval is 0.0005 seconds for all cases. The width and height of the tank are 1.08 m and 0.9 m, respectively.


The tank has forced harmonic roll motion with 1degree amplitude and 1.83rad/sec frequency, which is close to the lowest natural frequency of the middle layer. The interface elevation history is measured at leftside of tank and middle of tank at every time interval. The convergencetest result is shown in Figure 11. It shows that convergence is achieved for the cases of 0.0075 m and 0.005 m. From this point on, the results of 0.005 m will be presented. The simulated results are also compared with Molin et al.’s [13] experimental results.
Figure 12 shows the comparisons of three interfaces between the present simulation and Molin et al.’s [13] experiment. They generally show good agreement. The big kinks at trough in the experimental result of WD interface are assumed to be a problem associated with measurement. Since the oscillation is made near the lowest natural frequency of the middle layer, the corresponding lowest modal shape is generated for the middle layer (see Figure 14) and subsequently influences the other interfaces. The WD (waterdichloromethane) interface and CW (cyclohexanewater) interface have different amplitudes and they are larger than the actual oscillation of the top free surface. The freesurface oscillation looks more wiggly having both high and given frequencies and its phase is slightly shifted compared to others, which is also observed in the experiment. Figure 13 plots 3interface time histories at left tank wall and midtank. Interestingly, the WD interface at left tank wall has appreciably larger trough than crest and the WD interface at midtank has secondary hump near trough. Figure 14 shows two snapshots by the present MPS at 20 and 50 seconds. It is clearly seen that the internal sloshing motion is larger than that of top free surface.
(a)
(b)
(a)
(b)
3.5. More Violent Cases of ThreePhase Liquid Sloshing
In this section, two more cases are considered. The first case is forced roll harmonic motion with 2degree amplitude at 1.83rad/sec. The second case is more violent and bichromatic with 3degree roll amplitude at 1.83rad/sec and 4 cm sway amplitude at 3.62rad/sec (close to 3rd natural frequency of water layer).
Case 1. In this case, the roll amplitude is just doubled compared to the previous case. The oscillations are doubled at the free surface and WD interface, but the increase rate is less for the CW interface (see Figure 15(a)). The nonlinear phenomenon of secondary hump near trough at midtank is more pronounced for both WD and CW interfaces. An interesting point of this case is that wave breaking and mixture start to appear near midtank, as shown in Figure 16(b). The small internal breaking waves are reconstructed by buoyancy and surfacetension effects subsequently.
(a)
(b)
(a)
(b)
Case 2. In order to simulate a case of more violent interface motions, the tank is forced harmonically by roll amplitude of 3 degrees at 1.83rad/sec and sway amplitude of 4 cm at 3.62rad/sec. Figure 17 shows interface oscillations at left wall and midtank. Despite violent fluid motions, the time histories of interface elevations have repeatability. The increase of amplitude for internal interfaces is more toward trough side than crest side. The snapshots of Figure 18 show strong nonlinearity more clearly due to resonant bichromatic internal interfacial motions. One notable point is that the wave breaking is observed at the middle interface, but there is no breaking at the free surface. The severe internal motions can significantly hamper the effectiveness of separators or wash tanks although the 3deg roll angle is not considered severe.
(a)
(b)
(a)
(b)
Finally in Figure 19, we compare the two cases (a) with and (b) without buoyancy and surfacetension effects. They are snapshots at 23 seconds. Without the buoyancy and surfacetension model, the top layer penetrates much deeper with more particle detachment toward the middle layer. Their differences become smaller as nonlinearity decreases.
(a)
(b)
4. Concluding Remarks
The multiliquidsloshing simulation is investigated by developing a new multifluidlayer MPS. The new MPS method for multifluid system includes extra search methods for interface particles, boundary conditions for interfaces, buoyancycorrection model, and surfacetension model for interface particles. The new particle interaction models are verified through comparisons with published numerical and experimental data. In particular, the multiliquid MPS method is verified by comparing to the experimental results of Molin et al. [13] with three liquid layers. The significant point of multiliquid sloshing observed both in experiment and simulation is that the internal sloshing can be significantly larger than that of free surface when the frequencies of excitation are close to the natural frequencies of internal layers.
The verified multiliquid MPS program was subsequently used for more nonlinear cases including multichromatic multimodal motions with larger amplitudes, from which various nonlinear features, such as internal breaking and more particle detachment, can be observed. For this kind of nonlinear case, the differences between with and without buoyancycorrection and surfacetension models are demonstrated. The lessphysical phenomena of the latter can be fixed by the former. In case of internal resonance, the large internal motions can significantly hamper the effectiveness of separators or wash tanks although the imposed motion amplitude is not considered severe.
Conflict of Interests
The authors declare that they have no conflict of interests.
References
 H. C. Chen, “CFD Simulation of compressible twophase sloshing flow in a LNG tank,” Ocean Systems Engineering, vol. 1, no. 1, pp. 29–55, 2011. View at: Google Scholar
 B. H. Lee, J. C. Park, M. H. Kim, S. J. Jung, M. C. Ryu, and Y. S. Kim, “Numerical simulation of impact loads using a particle method,” Ocean Engineering, vol. 37, no. 23, pp. 164–173, 2010. View at: Publisher Site  Google Scholar
 B. H. Lee, J. C. Park, M. H. Kim, and S. C. Hwang, “Stepbystep improvement of MPS method in simulating violent freesurface motions and impactloads,” Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 9–12, pp. 1113–1125, 2011. View at: Publisher Site  Google Scholar
 T. Loysel, E. Gervaise, S. Moreau, and L. Brosset, “Results of the 20122013 sloshing model test benchmark,” in Proceedings of the 23rd International Offshore and Polar Engineering Conference (ISOPE '13), Anchorage, Alaska, USA, 2013. View at: Google Scholar
 Y. Kim, Y. Shin, W. Kin, and D. Yue, “Study on sloshing problem coupled with ship motion in waves,” in Proceedings of the 8th International Conference on Numerical Ship Hydrodynamics, Busan, Korea, 2003. View at: Google Scholar
 Y. Kim, B. W. Nam, D. W. Kim, and Y. S. Kim, “Study on coupling effects of ship motion and sloshing,” Ocean Engineering, vol. 34, no. 16, pp. 2176–2187, 2007. View at: Publisher Site  Google Scholar
 G. Gaillarde, A. Ledoux, and M. Lynch, “Coupling between liquefied gas and vessel’s motion for partially filled tanks: effect on seakeeping,” in Design & Operation of Gas Carriers, The Royal Institution of Naval Architects, London, UK, 2004. View at: Google Scholar
 S. Cho, S. Y. Hong, J. Kim, and I. Park, “Studies on the coupled dynamics of ship motion and sloshing including multibody interactions,” in Proceedings of the 17th International Offshore and Polar Engineering Conference (ISOPE '07), vol. 3, pp. 1900–1904, Lisbon, Portugal, July 2007. View at: Google Scholar
 S. J. Lee and M. H. Kim, “The effects of innerliquid motion on LNG vessel responses,” Journal of Offshore Mechanics and Arctic Engineering, vol. 132, no. 2, pp. 1–8, 2010. View at: Publisher Site  Google Scholar
 K. S. Kim, B. H. Lee, M. H. Kim, and J. C. Park, “Simulation of sloshing effect on vessel motions by using MPS (moving particle simulation),” CMESComputer Modeling in Engineering and Sciences, vol. 79, no. 34, pp. 201–221, 2011. View at: Google Scholar  MathSciNet
 K. S. Kim, M. H. Kim, and J. C. Park, “Dynamic coupling between ship motion and threelayerliquid separator by using MPS, (Moving Particle Simulation),” in Proceedings of the 23rd International Offshore and Polar Engineering Conference, Anchorage, Alaska, USA, 2013. View at: Google Scholar
 M. La Rocca, G. Sciortino, C. Adduce, and M. A. Boniforti, “Experimental and theoretical investigation on the sloshing of a twoliquid system with free surface,” Physics of Fluids, vol. 17, no. 6, pp. 1–17, 2005. View at: Publisher Site  Google Scholar
 B. Molin, F. Remy, C. Audiffren, and R. Marcer, “Experimental and numerical study of liquid sloshing in a rectangular tank with three fluids,” in Proceedings of the 22nd International Offshore and Polar Engineering Conference (ISOPE '12), Rhodes, Greece, 2012. View at: Google Scholar
 S. Koshizuka and Y. Oka, “Movingparticle semiimplicit method for fragmentation of incompressible fluid,” Nuclear Science and Engineering, vol. 123, no. 3, pp. 421–434, 1996. View at: Google Scholar
 E. Toyota, H. Akimoto, and S. Kubo, “A particle method with variable spatial resolution for incompressible flows,” in Proceedings of the 19th Japan Society of Fluid Mechanics, vol. A92, 2005. View at: Google Scholar
 H. Gotoh, “Lagrangian particle method as advanced technology for numerical wave flume,” International Journal of Offshore and Polar Engineering, vol. 19, no. 3, pp. 161–167, 2009. View at: Google Scholar
 S. C. Hwang, B. H. Lee, J. C. Park, and H. G. Sung, “ParticleBased Simulation for sloshing in a rectangular tank,” Korean Society of Coastal and Ocean Engineers, vol. 24, no. 2, pp. 964–989, 2010. View at: Google Scholar
 N. Shirakawa, H. Horie, Y. Yamamoto, and S. Tsunoyama, “Analysis of the void distribution in a circular tube with the twofluid particle interaction method,” Journal of Nuclear Science and Technology, vol. 38, no. 6, pp. 392–402, 2001. View at: Google Scholar
 K. Nomura, S. Koshizuka, Y. Oka, and H. Obata, “Numerical analysis of droplet breakup behavior using particle method,” Journal of Nuclear Science and Technology, vol. 38, no. 12, pp. 1057–1064, 2001. View at: Google Scholar
 A. Khayyer and H. Gotoh, “Enhancement of performance and stability of MPS meshfree particle method for multiphase flows characterized by high density ratios,” Journal of Computational Physics, vol. 242, no. 1, pp. 211–233, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 J. J. Monaghan and A. Kocharyan, “SPH simulation of multiphase flow,” Computer Physics Communications, vol. 87, no. 12, pp. 225–235, 1995. View at: Google Scholar
 X. Y. Hu and N. A. Adams, “A multiphase SPH method for macroscopic and mesoscopic flows,” Journal of Computational Physics, vol. 213, no. 2, pp. 844–861, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 M. Tanaka and T. Masunaga, “Stabilization and smoothing of pressure in MPS method by QuasiCompressibility,” Journal of Computational Physics, vol. 229, no. 11, pp. 4279–4290, 2010. View at: Publisher Site  Google Scholar
 Z. R. Kishev, C. Hu, and M. Kashiwagi, “Numerical simulation of violent sloshing by a CIPbased method,” Journal of Marine Science and Technology, vol. 11, no. 2, pp. 111–122, 2006. View at: Publisher Site  Google Scholar
 H. Lamb, Hydrodynamics, Cambridge University Press, Cambridge, UK, 1945. View at: MathSciNet
Copyright
Copyright © 2014 Kyung Sung Kim et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.