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 numerical-prediction tool for multiphase fluids is required. In this regard, a new moving particle simulation (MPS) method is developed to simulate multiliquid-layer sloshing problems. The new MPS method for multifluid system includes extra search methods for interface particles, boundary conditions for interfaces, buoyancy-correction model, and surface-tension 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 internal-layer resonances, the internal interface motions can be much greater than top free-surface 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 buoyancy-correction and surface-tension 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 multi-well 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 large-scale “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 inner-fluid motions cannot be ignored in predicting platform motions. Therefore, the analysis of the inner-fluid 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 multilayer-liquid system have rarely been attempted. Two-liquid-layer sloshing in a rectangular tank was theoretically and experimentally studied by La Rocca et al. [12]. Molin et al. [13] investigated three-liquid-layer 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 multiliquid-layer system and the simulation results are compared with the three-liquid-layer-sloshing 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 less-optimal treatments in the Poisson equation, gradient/collision model, and tracing method for free-surface 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 single-phase-liquid problems. In order to simulate multiphase-fluid motions with multiple interfaces, Shirakawa et al. [18] proposed a buoyancy model and Nomura et al. [19] suggested a surface-tension 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 free-surface particles. A surface-tension model is also developed for interface particles. A hydrostatic-correction 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 three-layer-liquid 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 Single-Liquid-Layer
In this section, the MPS method for single-liquid-layer 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 free-surface 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 Navier-Stokes 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 right-hand side of Navier-Stokes 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 marker-and-cell) 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 free-surface boundary condition is satisfied by taking the pressure on the free-surface particles to be the same as atmospheric pressure .
To apply the dynamic free-surface condition, the free-surface particles should be identified first. In the free-surface region, the particle density decreases since there are no particles in air region. Thus, the following simple conditions are used to identify the free-surface 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 free-surface 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 free-surface boundary condition, we can also simulate fragmentation and coalescence of free-surface 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 free-surface 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 free-surface particles. Since the sudden increase of particle number density influences identifying free-surface 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 optimal-value 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 free-surface particles.
2.6. Hydrostatic Pressure Correction Model
The hydrostatic pressure of the top-layer 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 hydrostatic-pressure correction is only applied to free-surface particles, but the detached/splashed/fragmented particles are excluded. By adding (13) to the free-surface-particle-searching 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 free-surface and HSC-required particle, its center pressure is given by, which is equivalent to applying the dynamic free-surface condition at the true free surface (top of free-surface particle):
2.7. Validation of Single-Phase Case
In order to validate the developed single-phase 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 Multiliquid-Layer
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 free-surface searching method, (12) and (13), can be used for interface. However, it is necessary to modify the free-surface searching method to be suitable for the interface tracing method.
The interface searching has one more condition compared to the free-surface searching, which is added to eliminate detached particles. According to (12), a particle can be identified as free-surface 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 free-surface 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 left-hand side of (19) to improve accuracy.
3.2. Buoyancy-Correction Model
MPS uses particle number density to interact between particles. However, this method can underestimate the self-buoyancy of center particle for multiliquid cases. Therefore, the correction of buoyancy effect is necessary for multilayer-liquid problems.
Since buoyancy force is related to pressure gradient, the buoyancy-correction 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 self-buoyancy of the center particle due to surrounding fluid of different density is not included because the kernel function is omni-spread from the center value with axial symmetry. In case of single fluid, the center particle is neutrally buoyant and the self-buoyancy correction is not necessary.
In order to find the self-buoyancy 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 self-buoyancy of the center particle can, for example, be obtained by subtracting the pressure of upper particles from that of lower particles. Therefore the self-buoyancy 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 right-hand sides of (20) are already included in the pressure gradient model of MPS method, while the second terms of the right-hand sides of (20) can be implemented into the pressure gradient model as buoyancy-correction 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 buoyancy-correction model can be expressed as follow: where is the kernel function for buoyancy-correction 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 velocity-adjustment 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 buoyancy-correction 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 single-phase fluid, the buoyancy-correction term vanishes.
In order to validate the developed buoyancy-correction 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/m3) than surrounding water particles (1000 kg/m3). The oil particle is to freely rise up to the free surface by self-buoyancy effect. In order to avoid side-wall 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 buoyancy-correction model does not show that trend, as shown in Figure 7. Due to the underestimation of the self-buoyancy effect, the rising velocity is extremely slow and irregular. By applying the present buoyancy-correction model, the oil particle rises with constant vertical velocity induced by the net self-buoyancy effect. The necessity of the buoyancy-correction model can be proved through this example. The rising velocity can be adjusted by using the velocity-adjustment parameter, as shown in the figure. As the parameter increases, the rising velocity is also increased. If the oil-particle ( 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 velocity-adjustment 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 surface-tension model, vibration of an ethanol droplet in two-dimensional 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/m3, the kinematic viscosity is m2/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. Three-Phase-Sloshing Example
For further verification, the developed multiphase-liquid MPS program is run to compare with Molin et al.’s [13] three-liquid-layer-sloshing experiment. The detailed information of the three-liquid-layer 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 1-degree amplitude and 1.83-rad/sec frequency, which is close to the lowest natural frequency of the middle layer. The interface elevation history is measured at left-side of tank and middle of tank at every time interval. The convergence-test 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 (water-dichloromethane) interface and CW (cyclohexane-water) interface have different amplitudes and they are larger than the actual oscillation of the top free surface. The free-surface 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 3-interface time histories at left tank wall and mid-tank. Interestingly, the WD interface at left tank wall has appreciably larger trough than crest and the WD interface at mid-tank 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 Three-Phase Liquid Sloshing
In this section, two more cases are considered. The first case is forced roll harmonic motion with 2-degree amplitude at 1.83-rad/sec. The second case is more violent and bichromatic with 3-degree roll amplitude at 1.83-rad/sec and 4 cm sway amplitude at 3.62-rad/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 mid-tank 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 mid-tank, as shown in Figure 16(b). The small internal breaking waves are reconstructed by buoyancy and surface-tension 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.83-rad/sec and sway amplitude of 4 cm at 3.62-rad/sec. Figure 17 shows interface oscillations at left wall and mid-tank. 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 3-deg 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 surface-tension effects. They are snapshots at 23 seconds. Without the buoyancy and surface-tension 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 multiliquid-sloshing simulation is investigated by developing a new multifluid-layer MPS. The new MPS method for multifluid system includes extra search methods for interface particles, boundary conditions for interfaces, buoyancy-correction model, and surface-tension 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 buoyancy-correction and surface-tension models are demonstrated. The less-physical 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.