Abstract
The smoothed particle hydrodynamics (SPH) method is applied to study the oil film diffusion in the water. By modifying the SPH equations of fluid dynamics, the multiphase flow SPH equations are obtained to establish the computational oil film diffusion model. By discussing three kinds of particle pairing schemes in the calculation of oil particle density, the redistribution mode of particle density is determined. The diffusion process of oil film is simulated, the effects of oil viscosity coefficient and particle density on oil film diffusion are analyzed, and the distribution of local pressure near oil particles in the process of oil film spreading is calculated. Finally, the calculated value of the oil film expansion diameter is compared with two other numerical models, and the calculated result shows a high coherence with the others.
1. Introduction
Smoothed particle hydrodynamics (SPH) is a meshless Lagrange algorithm, which was first applied to astrophysics and then was widely used in solid mechanics and fluid mechanics. It was put forward by Lucy and Gingold and Monaghan [1, 2]. This method was initially applied to solve the problem of arbitrary flow in three-dimensional space under the condition of no boundary in astrophysics. In the SPH method, the system state is described by a series of particles with properties. The motion of the particles is in accordance with the law of the conservation governing equation. The SPH method does not need mesh but uses an integral kernel called “kernel function” to approximate the “kernel function estimation.” In the particle support domain, the field function is presented by the integral approximation to smooth function, that is, the kernel approximation. The integrals of function and its derivative are replaced by superposition of the values of adjacent particles in the local region, that is, the particle approximation. This is the particle approximation process of hydrodynamic equations. The SPH method is adaptive; the arbitrary distribution of particles in an arbitrary shape region has no effect on the construction of the equation, so it is suitable for dealing with problems with large deformation. The particularity of the pure Lagrangian particle method is that the fluid in the computational area is replaced by particles, and the particles have properties of the fluid, so the SPH method is more flexible than other methods.
In the field of hydrodynamics, based on continuous medium hypothesis, the continuum is regarded as a series of particles, which are used to simulate the fluid motion. The SPH method simulates the fluid motion with a series of particles, and each particle has its own influence area and interpolation range. The position, velocity, and pressure of the particles and the gradient distribution of physical quantities are obtained by interpolation of the kernel function.
In recent years, the SPH method has been applied to a wide range of problems, which include fluid flows [3–5], wave-structure interaction problems [6, 7], water wave dynamics [8], geophysical flows [9, 10], and industrial applications [11–13]. Some scholars at home and abroad apply the SPH method to many different fields: Groenenboom et al. presented a coupled SPH-FE software tool useful for a large variety of dynamic fluid flow simulations and fluid-structure interaction [14]. Rahmat and Yildiz presented a multiphase ISPH method to simulate the droplet coalescence and electrocoalescence [15]. Ozbulut et al. investigated the wave characteristics in oscillatory motion of partially filled rectangular tanks by the SPH method [16]. Duan et al. presented an incompressible SPH algorithm by relaxing the density-invariant condition for incompressible flow [17]. Hopp-Hirschler et al. presented a smoothed particle hydrodynamics model for thermocapillary flow driven by gradients of the surface tension [18]. Qiu et al. presented a three-dimensional two-way coupled method to simulate moving solids in viscous free surface flows [19]. Liu and his collaborators used the SPH method to simulate a series of underwater explosions [20, 21]. Aristodemo et al. investigated the diffusive term in the continuity equation in order to smooth out the high-frequency numerical noise in the pressure field and simulated wave pressures acting on vertical and slotted coastal structures [22]. Ren presented a periodic density reinitialization smoothed particle hydrodynamics (PDRI-SPH) method to treat the generalized Newtonian free surface flows [23]. Shadloo et al. presented a 2D Lagrangian two-phase numerical model to study the deformation of droplet suspended in a quiescent fluid subjected to the combined effects of viscosity, surface tension, and electric field forces [24]. De Padova et al. simulated the process of transition from supercritical to subcritical flow at an abrupt drop by the SPH method [25]. Vacondio et al. [26–28] improved the closed boundary conditions of SPH with the use of virtual particles and enhanced the resolution of the small depth problems by using the particle splitting process with the shock capturing method.
The conventional SPH method is associated with low accuracy, so some different approaches have been proposed to improve the SPH approximation accuracy, including the finite particle method (FPM) by Liu et al. [29], the kernel gradient free SPH (KGF-SPH) by Huang et al. [30], and the decoupled finite particle method (DFPM) by Zhang et al. [31]. Zhang and Liu [32] have discussed these approaches and found that the DFPM inherits advantages from the conventional SPH and is flexible, effective, and easy in coding and no matrix inversion is required in DFPM, which is different from FPM, CSPM, and KGC.
In this paper, the SPH method is applied to the study of multiphase flow, and the extension of oil film on water is simulated and analyzed. It is of great significance for the analysis and treatment of oil spill accidents in the sea. Density is an important parameter of the model in the process of calculating oil diffusion, and three kinds of particle pairing schemes for calculating oil particle density were discussed; in oil-oil particles pairing scheme, the oil particle diffused better and the expanding diameter increased gradually with time. The calculated diameter was compared with two other numerical models. The results show that the viscosity coefficient of oil has little effect on oil film extension diameter, and the diffusion rate of light oil particles is faster than that of heavy oil. With the oil film spreading, the pressure distribution is gradually stable, and the pressure increases from the top of oil film to the bottom of water.
2. Materials and Methods
The smoothed particle hydrodynamics (SPH) method is a kind of Lagrangian particle method for simulation of fluid motion. It does not need grids to partition the computing domain; instead, the computational domain is discretized into a series of interaction particles, which carry a variety of physical quantities and interact with each other by the kernel function. In this paper, the SPH method is applied to study the oil film diffusion in water, and the calculation flowchart of SPH methodology is shown in Figure 1.
2.1. Theoretical Basis of SPH Method
The construction of SPH function is divided into two steps: the first step is the integral representation, that is, kernel interpolation; the second step is particle approximation. Kernel function interpolation implements interpolation of field function or the gradient of field function, while particle approximation realizes particle discretization of integral expression of kernel function estimation.
2.1.1. Integral Interpolation
Interpolation is the core of the SPH method; as for a continuous smooth function , the value of function in the domain can be expressed aswhere the parameter h, called smoothing length, is used to describe the size of supported domain; the kernel function W is an even function that satisfies the regularization condition and the compact support condition, and it is a Dirac delta function when h ⟶ 0; is coordinate vector; and represents surface differential and volume differential, respectively, in two and three dimensions.
The kernel approximation equation (2) of the first derivative of can be obtained by bring as a function into equation (1):
When the support region is located inside the computational domain or with periodic boundary conditions, equation (2) can be written in equation (3) which is the integral representation of functional derivative.
2.1.2. Particle Approximation
The SPH method is represented by finite particles with independent mass and independent influence domain in the whole computational region. The particle approximation is based on the integral approximation equation by applying the corresponding values of adjacent particles in the local region to superimpose and replace the integral form of function and its derivative. Approximation of particle “a” is as shown in Figure 2. The approximation of particle “a” needs all particles in a circular supporting domain with a radius of “kh.”
Equation (1) can be expressed in the form of all particle summation in the support domain; the transformation process is as follows:
If , , the particle approximation equation of field function is expressed as
2.1.3. Kernel Function
The SPH method is expressed by integration of kernel function. The kernel function transforms complex particle space relations into simple interpolation weight number and transforms the gradient distribution of particle field function into the derivative of the kernel function for calculation. In this paper, the cubic spline kernel function is used as kernel function.where q is the dimensionless distance (relative distance) between two particles , r is the distance between two particles, and h is the smooth distance; the values of are , , and in one, two, and three dimensions, respectively.
Figure 3 shows the cubic spline kernel function distribution; it can be seen that the function is continuous and smooth. Figure 4 shows the contour distribution of cubic spline; the value of decreases slowly from 1 to 0.25, and it decreases quickly from 0.25 to 0; 0.25 is the function value at piecewise point.
2.2. Governing Equation
2.2.1. Continuity Density Equation
The governing equation of the model, Navier–Stokes equation, is taken as the basic equation to describe the fluid motion. The idea of the SPH method is to discretize the partial differential equations of N-S equation in space and get the SPH governing equation suitable for generalized fluid dynamics in Cartesian coordinates.
The continuity equation based on conservation of mass is
The density is included in the gradient operator by using the law of multiplication of partial differential equations:
According to the concept of SPH approximation, the particle approximation of velocity and divergence in governing equation is carried out, and the continuity density equation is obtained [33]:where and subscript means that both particles a and b are from the same fluid phase, while indicates the other way.
2.2.2. Momentum Equation
In the SPH method, the general form of the momentum equation in the Navier–Stokes equation in Cartesian coordinate system iswhere is dissipation, and different forms of momentum equations can be established based on different dissipative formulas. In this paper, the subparticle scale, which is a turbulence simulation method proposed by Gotoh in 2001, is used to describe the turbulence in the MPS model. The scale of turbulence simulated by the SPS method is smaller than that by the SPH method; it is a method based on large eddy simulation (LES). In this method, the N-S equation is filtrated from the particle distribution, and the subparticle stress term is added. Therefore, the emphasis of coupling large eddy simulation with particle method is to establish the subparticle stress model and express it by the particle method.
The momentum conservation equation is as follows:
The stress term is calculated by the implicit formula of the interaction force between particles with a fixed length of density, and the initial term is calculated by the explicit formula. is laminar flow term, and it can be simplified to Lo and Shao [34]:where is the kinematic viscosity coefficient, which usually takes .
In equation (11), is the SPS stress tensor with elements given by and satisfies the Boussinesq vortex-viscosity hypothesis and is calculated using the Favre averaging method of compressible fluids:where is the subparticle shear stress; is the turbulent eddy viscosity, , in which is the distance between particles, is the SPS strain tensor at the position (i, j), is the Smagorinsky constant, with the value of 0.12, and the local strain rate is defined by , ; and is turbulent kinetic energy in which is the average value of the product of the turbulent velocity of two particles.
The particle approximation of momentum equation can be expressed as
2.3. Multiphase Flow SPH Equations and Correlated Particle Equations
2.3.1. Multiphase Flow SPH Equations and State Equation
The traditional SPH method updates the physical quantity by the particle density in the smoothing domain, which leads to the traditional SPH method only suitable for solving the problem of a single medium. In the solution of multiphase flow problems, the computational instability is easy to occur, so it is necessary to improve the traditional SPH method.
The modified SPH equations for multiphase flow proposed by Ott can be used to simulate the multiphase flow effectively [35]. The equations are as follows:
In the SPH approach, the fluid is assumed to be weakly compressible and the pressure is determined via the state equation as follows [36]:where B is the numerical parameter, , in which is reference density and is the maximum velocity of liquid; and is a constant, usually .
2.3.2. Particles Displacement Equation
Considering the stability of the calculation, the integral terms such as displacement, velocity, and density are usually interpolated in the process of numerical integration. The X-SPH technology proposed by Monaghan was adopted to calculate displacement integral [37].where and is an adjustable parameter.
2.3.3. Particles Energy Equation
The energy equation of each particle can be expressed aswhere P is the pressure; is the velocity; is the thermal energy; and is the interpolating kernel. The problem of nonconservation of quality can be solved by using the modified kernel function to build noncompressible solution.
In SPH simulations, there may be large fluctuation in the pressure field of particles because the density error accumulates when time marches. In order to reduce the fluctuation, the density of each particle is modified by the zero-order Shepard particle density filter, which was introduced by Colagrossi and Landrini [38].
2.3.4. Surface Tension
In the modified SPH model, the surface tension is mainly determined by the particles interaction forces. Given the interaction force to the particles in a certain range, the resultant force of the particles in the fluid is zero and the resultant force of the particles at the interface is not zero; thus, the surface tension is calculated. For the two adjacent particles a and b, when , there are two modes of interaction between particles, that is, the repulsive force between particles from different materials [21] and the molecular force between particles from the same material [39].
The repulsive force between particles is calculated as follows.
A highly peaked repulsive force is applied to particles from different materials near the interface, which are approaching and tending to penetrate each other. The penetration pe is detected if
The repulsive force is similar to the molecular force of Lennard-Jones form and is applied pairwise on the two approaching particles along the center of the two particles.where is the distance between particles a and b; h is the smooth length; is the coordinate vector, ; D is the coefficient characterizing the wettability of the interface, D = 105 N·m; and the parameters and are taken as 6 and 4, respectively.
The molecular force between particles is calculated as follows:where is the strength of the force acting between particles a and b and the resultant force of surface tension acting on a is generally calculated by the formula . However, when dealing with the problem of fluid motion with large Weber number, the calculation of surface tension is limited to the range of , which may lead to intersurface fracture or small surface tension due to the lack of particles. Therefore, the surface tension per unit mass of particle a is modified by interpolation:where the calculation range of particle a surface tension can be adjusted to as needed.
The momentum equation after considering the effect of surface tension is as follows:
2.4. Boundary Treatment
In the SPH method, the free surface boundary is automatically captured, and the boundary condition problem is mainly focused on the treatment of the solid boundary. The boundary normal force method proposed by Monaghan is applied in this paper [40]. It is considered that the solid wall particles only provide normal force to water particles, but not tangential force, as shown in Figure 5.
A series of points are arranged on the solid boundary, and the boundary points are assumed to impose an appropriate central repulsive force on the fluid particles nearby. The purpose of this boundary condition is to prevent nonphysical penetration of the solid boundary by the fluid particles. The unit mass forces of fluid particles and boundary particles are in accordance with the Lennard-Jones potential energy model. By interpolating the solid particles to calculate the pressure, the effect of the particle spacing on the wall of the repulsive force is reduced.
The unit mass force on the fluid particle from the boundary particle iswhere is a local unit normal vector that points from the boundary into the fluid; y is the normal distance of the fluid particles to the solid boundary; x is the tangential distance; is the normal projection of fluid velocity; ; the subscripts wp and bp represent fluid particle and boundary particle, respectively; and z is the water level of hydrostatic water depth . is a repulsive function, which can be expressed with the normalized distance of fluid particles from solid boundary aswhere in which is the initial distance between two particles and in which is the acoustic velocity of particle i.
The pressure function is defined by the rule:
The pressure function can ensure that the fluid particles are subjected to a constant repulsion force when they moved parallel to the solid boundary.
is the correction function, and it adjusts the force according to the local water depth and the normal velocity of the fluid particle; the expression is as follows:where is the reference acoustic velocity at reference density, .
In the two-dimensional calculation, the normal vector at the boundary particle j is calculated by the position of boundary particles j − 1 and j + 1 on both sides of j and the relationship between unit tangent vector and normal vector is , where .
3. Results and Discussion
3.1. Oil Particles Density Redistribution
Density is an important parameter in the calculation of oil film spreading motion by the SPH method. The density estimation depends on the distribution of particles in the supporting domain. The particle pressure is calculated by density, and then the position and velocity of particles are calculated. If the average density is calculated for different particles in each oil particle support domain, different motion principles can be reflected.
In this paper, the processing methods of oil particle density are divided into the following three cases:(a)Considering only the nearby oil particles density, and the effect of water particles density is not considered(b)Considering only the nearby water particles density, and the effect of oil particles density is not considered(c)Considering both oil and water particles density at the same time
The experimental tank is 7.5 m × 0.18 m, a 0.5 m × 0.05 m oil block is placed in the middle of the tank, and the distance between the center of the oil block and the left wall of the water tank is 3.8 m. The total number of particles in the model is 14351, of which 14097 are water particles and 254 are oil particles, and the particles distance is 0.01 m. The time step is .
Figure 6 shows the oil film state at 0.1 s of three calculation modes. Figure 6(a) shows an oil-oil particles pairing model; it can be seen from the figure that the distance between the oil particles and water particles is large, the distribution of the oil particles is dispersed, and the diffusion of the oil particles is good, but the arrangement of the water particles is not in good order. Figure 6(b) shows an oil-water particles pairing model; the oil and water particles are arranged in order, especially the oil particles are evenly distributed on the oil and water interface. Figure 6(c) shows the oil particles paired with the nearby oil and water particles, and the simulation results are similar to that of Figure 6(b); the oil and water particles are arranged orderly. It can be seen from the figure that the difference of simulation results between searching for oil and water particles together and searching water particles only is not obvious, and the oil film expands faster for only searching for oil particles than the first two cases. Figure 7 shows the simulation results of oil extension in different modes at 4 s. In Figure 7(a), the oil film spread on the water surface and the oil particles are distributed well; in Figure 7(b) and Figure 7(c), the oil film spread on water surface is not obvious, and the oil particles are arranged closely.
(a)
(b)
(c)
(a)
(b)
(c)
When considering the density of water and oil particles, oil particles form a similar molecular force effect, which inhibits the extension of oil film. In model (a), the oil film expands continuously, the oil particles at the center of the oil film are concentrated slightly, the thickness of the oil film decreases, and the extension diameter increases with time. In the later period of diffusion of the three models, the particles in the interface are well distributed, and the arrangement of particles in model (a) is no longer disordered. Therefore, the simulation of oil film diffusion on water is based on the calculation model of oil-oil particles pairing.
3.2. Oil Film Diffusion Process
The diffusion model of oil spill film in water by the SPH method is composed of oil and water particles. For the simulation of two-phase flow with different densities, oil particles and water particles need to be calculated separately in the equation of state. The experimental tank is 7.5 m × 0.18 m, a 0.5 m × 0.05 m oil block is placed in the middle of the tank, and the distance between the center of the oil block and the left wall of the water tank is 3.8 m. Because the volume of the oil film is small relative to the model volume of the water tank, the number of particles is small; therefore, the zero-order filter correction is performed on the density every 30 steps. The particles distance is 0.01 m, and time step is . In order to avoid calculation instability, the side wall particles are taken as the same density as water particles.
The simulation result is shown in Figure 8. Initial field is square oil block; during 0.1 s and 0.15 s, the oil block sinks briefly due to the action of gravity. At 0.75 s, the buoyancy of the oil block gradually balances with the gravity, and the oil block floats slowly on the surface. After that, the oil film diffused around due to the pressure gradient caused by the difference of the gravity centers between the oil film and the water.
When the oil film is gently placed on the water surface, the oil film will sink for a short time due to the gravity and inertia force of the oil film, but because the oil density is less than the water density, the oil film will not sink into the water but float on the water surface. The gravity center of the oil film is higher than that of the water, and the pressure gradient produced by the difference of the gravity center leads to the diffusion of oil film.
3.2.1. Effect of Viscosity Coefficient on Oil Film Expansion
There may be many different kinds of oil when oil film accidents happened, such as gasoline, diesel, and petroleum. Different kinds of oil have different kinematical viscosity coefficients. In this paper, gasoline and petroleum are selected to simulate. The viscosity coefficients of gasoline and petroleum are 0.65 times and 25 times that of water. Figure 9 shows the expansion diameter of oil with different kinematic viscosity coefficients.
In Figure 9, the red and blue lines represent the expansion diameters of gasoline and petroleum. Before 5 s, the relationship between oil expansion diameter and gasoline expansion diameter was not obvious, while after 5 s, the expansion diameter of gasoline was larger than that of petroleum. The correlation of the two curves was good. It can be seen that the motion viscosity coefficient is inversely proportional to the spread diameter, but the value of the motion viscosity coefficient has little effect on the oil film expansion.
3.2.2. Comparison of Oil Particles Diffusion with Different Densities
Violeau pointed that when densities of oil and water are close, there are no additional problems studying oil film diffusion by the SPH method [41]. When oil spills occur, there are many kinds of oil products. In this paper, two kinds of oil spill with different densities are simulated and compared with each other. The light oil density is 830 kg/m3, and the heavy oil density is 995 kg/m3. Figures 10 and 11 show the diffusion processes of light oil and heavy oil.
It can be seen from Figures 10 and 11 that gravity plays a major role in the movement of oil particles from 0 s to 0.25 s. The oil particles are in the stage of gravity falling and move in the negative direction of the y axis. The light oil density is smaller than that of water, and there is obvious stress instability on the interface between oil and water. The particle distributions between oil-oil, water-water, and oil-water particles are inhomogeneous, especially at 0.1 s, and the oil-water particles contact closely because the oil particles sink to the limit of water compression, while due to the large difference of the densities between the two kinds of particles and the obvious pressure difference calculated by state equation, the shear stress produced by the two particles increases the vortex at the interface, and there is a slight bouncing between the two kinds of particles. At this stage, the gravity of heavy oil is greater than buoyancy, and then it sinks into the water. Because the density of heavy oil is closer to water, the density of the model is well distributed, the stress instability phenomenon on the oil-water interface is not obvious, and the oil and water particles distributed compactly. From 0.25 s to 5 s, the light oil particles mainly expand in the x direction, and the particles are gradually well distributed; the spacing between oil particles is larger than that of water particles; heavy oil particles are mainly diffused horizontally, the gravity center of heavy oil particles is lower than that of light oil particles, and the gravity center of heavy oil particles is always above the water surface. The diffusion velocity and diameter of light oil particles are larger than those of heavy oil particles.
In the process of oil diffusion, the particles distance between oil and water surface is longer compared with other particles. The main reason is there is a repulsive force at the interface between the two fluids in order to prevent the penetration between the two phases, and this may lead to the big distance between oil and water particles at the interface.
3.2.3. Oil Spill under Water
The experimental tank is 7.5 m × 0.18 m, a 0.5 m × 0.05 m oil film is placed in the middle of the tank and under the water surface, the distance between the canter of the oil block and the left wall of the water tank is 3.8 m, and the particles distance is 0.01 m. The time step is . Figure 12 shows the diffusion process of oil film under water on static water. From 0 s to 0.15 s, because the density of oil is smaller than that of water, the oil film is mainly floating under the effect of buoyancy, and the horizontal expansion is not obvious. After 0.15 s, the oil film spreads on the water and the thickness of oil block gradually decreases. The oil film diameter extends from 0.6 m to 2.5 m between 0.15 s and 4 s. The expansion process of oil film diameter is similar to that of oil spill on water.
3.3. Local Pressure Distribution
Figure 13 shows the pressure distribution and variation process during the oil film expansion with a density of 830 kg/m3. At the beginning of calculation, the oil block has a transient subsidence trend, and it is deformed slightly; the contact surface between the oil block and water oscillates. The oil block is slightly placed on the water surface at 0 s, and the pressure distribution between the oil and water is in accordance with the distribution of hydrostatic pressure. The surface pressure is 0 pa, and the pressure increases with the depth. At 0.05 s, when the oil block is in contact with water, the oil block sinks momentarily and vibrates slightly under the action of gravity and the inertial force, and the particle density changes greatly. There is a large negative pressure at the interface between oil and water, especially near the center of the lower surface of the oil block, and there is a small negative pressure in the water surface near the oil block. At 0.1 s, the oil block has a slight rebound, the stress at the bottom of the block is released along the water toward the two ends of the block, and the negative pressure gradually decreases. At 0.25 s, the oil block floats slightly to a more balanced state, the negative pressure was released to both sides and decreased, and the interface between oil and water appears a high-pressure region. From 0.25 s to 0.45 s, the oil block slightly oscillates again. The negative pressure on the lower surface gradually decreases due to the spreading of the oil film, and the upper surface of the oil film has a small negative pressure due to the decreasing number of particles.
From 0 s to 0.45 s, the oil film extension is mainly affected by gravity and inertia force and oscillates continuously in the balance of the two forces. In this process, the diameter of oil film increases continuously, and the oil-water interface oscillates continuously, which is caused by the diameter enlargement and results in negative pressure at the oil-water interface and high pressure of particles at the interface. The hydrodynamic pressure transmission is more obvious among the oil particles but not transferred into the water, and the pressure of particles at the oil-water interface is high. The transmission of dynamic pressure to the underwater is relatively small, and the local pressure of the underwater is disordered. After 0.45 s, the oil film extends on the water surface. In this process, the diameter of oil film increases with the decrease of thickness. Oil particles increase laterally and decrease longitudinally. The pressure distribution is gradually stable, showing the trend of increasing pressure from the upper surface to the lower surface of the oil film and water.
3.4. Modeling Verification
In this paper, the oil film expansion diameter of this model is compared with the Fay model (1971) [42] and Xiaokong Liu model (1981) [43].
The simulated and two theoretical curves of oil film expansion diameters with time are shown in Figure 14; the red line is the calculated value of oil expansion diameter of the oil film diffusion model, the blue line is the theoretical value of oil expansion diameter of the Fay model [42], and the green line is the theoretical value of oil expansion diameter of the Xiaokong Liu model [43].
Comparing the calculated value with the Fay model, it can be seen from the figure that in the first two seconds, the simulation value of the expanded diameter is slightly smaller than the theoretical value, the absolute error is between 0.1 m and 0.6 m, and the relative error is between 4% and 50%. As the expansion diameter of the oil film becomes larger, the relative difference between the calculated value and the theoretical value decreases. After 2 s, the relative difference between the theoretical value and the simulated value is relatively small, and the relative error is between 1% and 12%. The reason for large relative error in the early stage of simulation may be that the initial diameters of simulated and theoretical values are different and the initial diameter of the oil block is 0.5 m while the theoretical value is 0 m.
Comparing the calculated value with the Xiaokong Liu model, it can be seen from the figure that in the first second, the calculated value is closer to this model compared with the Fay model and the maximum relative error is 34.3%, while from 1 s to 23 s, the absolute error becomes larger and then becomes smaller, the maximum absolute error is 0.36 m, and the maximum relative error is 10.3%. By comparing the oil film diameter of the two models, the oil film diameter expand trend of the three models is the same, the simulation value of the model has a good correlation with the theoretical value, and the model can be used to simulate the oil film diffusion.
4. Conclusions
(1)Due to the good adaptability of the SPH method, the random distribution of particles in any shape calculation region will not have much influence on the construction of the equation, especially for the numerical simulation problem with large deformation. In this paper, the SPH method is applied to the study of multiphase flow, and the expansion of oil film on water is simulated and analyzed. It is of great significance for the analysis and treatment of oil spill accidents in the sea.(2)By interpolation and particle approximation of the smooth kernel function, the hydrodynamic SPH formula is derived and modified to obtain the multiphase flow SPH equations, which are applied to the simulation of oil film diffusion in water. Three kinds of particle pairing schemes for calculating oil particle density are discussed. In the pairing scheme of oil-oil particles pairing, the oil particles diffuse better, the oil film thickness decreases gradually to both sides, and the expanding diameter increases gradually with time; this scheme is adopted in the model.(3)The model simulates the inertial sinking, floating, and spreading processes of oil film. The effects of oil viscosity coefficient and oil particle density on oil film diffusion are studied. The results show that the viscosity coefficient of oil has little effect on oil film expansion diameter; the stress distribution of heavy oil is more stable than that of light oil, and the density of heavy oil is more uniform. The diffusion rate of light oil particles is faster than that of heavy oil.(4)The local pressure near the oil particles in the process of oil film diffusion is calculated. With the oil film spreading, the pressure distribution is gradually stable, and the pressure increases from the top of oil film to the bottom of water.(5)The calculated values of the oil film expansion diameter are compared and verified with the Fay model and Xiaokong Liu model, the simulation value of the model has a good correlation with the theoretical values, and the model can be used to simulate the oil film diffusion.
Data Availability
The data cannot be shared because of the cooperation unit confidentiality requirements.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
Daming Li and Zhu Zhen contributed to the study design. Zhu Zhen contributed to the data collection, performed the data analysis, and wrote the manuscript. Hongqiang Zhang, Yanqing Li, and Xingchen Tang contributed to the discussion and interpretation of the results. All authors participated in reading and finalizing the manuscript.
Acknowledgments
This research was financially supported by the National Natural Science Foundation of China Innovative Research Group Science Fund Project, “Basic Theoretical Research on the Safety of Major Hydraulic Engineering” (51321065). This study was also supported by the State Key Laboratory of Hydraulic Engineering Simulation and Safety of Tianjin University.