Abstract
Improving the physical realism of oil spill scenes in marine simulators can further enhance the emergency response capabilities of officials in charge and crew members and help reduce losses caused by oil spill disasters. In order to uniformly simulate the spreading, drift, breakup, and merging of oil spills at sea, we propose an improved divergence-free position-based fluid (DFPBF) framework based on the particle number density model. In our DFPBF framework, the governing equations for oil spills and ocean are discretized by Lagrangian particles, and the incompressibility of oil spills and ocean is ensured by solving the divergence-free velocity constraint solver and constant density constraint solver. In order to stably simulate the fate and transport of oil spills with higher viscosity, we introduce an implicit viscosity solution scheme for our DFPBF framework. The implicit solver uses a splitting concept to decouple viscosity solve and adopts an implicit scheme to discretize the integration of viscous force. Moreover, our DFPBF framework can ensure a divergence-free velocity field before applying the implicit viscosity scheme, which avoids the undesired bulk viscosity effects. The simulation results show that our DFPBF framework can stably simulate oil spills of various viscosities, especially high-viscosity crude oils.
1. Introduction
Oil spills at sea are one of the most dangerous maritime disasters. Oil pollution has a serious hazard to the marine ecological environment and marine biological resources. Oil spill pollution accidents can cause vast stretches of ocean to become hypoxic, leading to the smothering of sedentary species, such as seaweeds and shellfish [1, 2]. The fishing and mariculture industries will also be seriously affected by oil spills. If the spilled oil is not cleaned up in time, it may also cause explosions and fires, causing more serious economic losses and casualties. Oil spills include releases of crude oil from tankers, offshore platforms, and drilling rigs. In this paper, we mainly focus on the ship-source oil spills. When an oil spill occurs, if the competent authority and ship’s officers can respond effectively and successfully, the damage caused by oil spill will be significantly reduced. However, effective and successful emergency response capability is usually obtained from a lot of practical experience or simulation training [3]. It is unrealistic to accumulate practical experience from actual oil spills because it will often bring huge losses and even dangers, and no mistakes are allowed. The advantage of marine simulator is that it can build virtual oil spill scenes and can simulate oil spill scenes under different weather and sea conditions according to user’s requirements. There is no risk during the training process, and the emergency response capability of competent authority and officers can be trained through a large number of repeated simulation training.
The main study purpose of this paper is to physically model the oil spill scenes at sea and improve the realism of the oil spill scene in the marine simulator so that the marine simulator can play a more important role in the oil spill emergency response training. The fate and transport of oil spills in marine environment include spreading, drift, evaporation, dispersion, emulsification, dissolution, photo-oxidation, sedimentation and sinking, shoreline interaction, and biodegradation. During the early stages of a spill, spreading, drift, evaporation, dispersion, emulsification, and dissolution are most important [4]. In this paper, we mainly focus on the physical movement process of oil spills in the early stage, that is, spreading and drift. The main theme of this paper is as follows:(1)In this paper, we introduce an improved DFPBF (divergence-free position-based fluid) framework to model oil spill scenarios. The DFPBF framework is the application of the concept of position-based dynamics in fluid simulation. The DFPBF framework is a meshless method that discretizes the problem domain into a series of fluid particles carrying physical quantities, which has great advantages in simulating the spreading and drift of oil spills. In order to simulate the coupled movement of the oil spills and the ocean surface, we propose an improved particle number density model. All particles in the support domain, including oil particles and ocean particles, have the same rest density, which can effectively handle the discontinuities of physical quantities at the interfaces.(2)The viscosity of the oil has a great effect on the spreading and drift. Usually, low-viscosity oils spread much faster than those with high viscosity. Moreover, high-viscosity oils are more resistant to deformation during spreading and drift. The viscosity of oil will be very high in low-temperature environment. Although the explicit viscosity solvers are usually more efficient, it requires a smaller time step to ensure stability when simulating higher viscosity oils, which will decrease the performance. In addition, when the viscosity of the oil exceeds a certain range, it is difficult to maintain stability even if the time step is small. Therefore, in order to stably simulate oils with different viscosities, we introduced an implicit viscosity solution scheme for the DFPBF framework. Compared with the explicit solvers, the implicit viscosity solver can stably simulate oil spills with very high viscosity without being limited by the time step. Before applying the implicit viscosity solver, our DFPBF framework ensures a divergence-free velocity field through the velocity-based constraints, which avoids the undesired bulk viscosity effects. The implicit viscosity solver uses the matrix-free conjugate gradient method and the Jacobi preconditioner to accelerate convergence, which significantly improves the computation efficiency.
2. Related Work
2.1. Spreading Models of Oil Spills
Once an oil spill occurs, the oil will start to spread. The spreading speed of the spilled oil is related to the viscosity. Refined products or light crudes with lower viscosity spread faster, and their thickness will gradually reduce as it spreads. In contrast, heavy crudes with high viscosity spread slowly and fragment into patches easily [4]. Oil spill spreading models mainly include regression models and numerical models. Early oil spill spreading models are mainly regression models, which usually need to fit the experimental or observed data. Fay [5, 6] proposed a three-regime model for the spreading of oil spills. The oil spill spreads under the action of gravity, viscosity, and surface tension. On the basis of Fay’s work, Mackay et al. [7] considered the effect of oil thickness. However, both Fay and Mackay thought that the oil spills were circular in shape, which is not consistent with the actual observation. Lehr et al. [8, 9] considered the effect of wind and Langmuir circulation on oil spreading and thought that the oil did not spread uniformly. The modified model proposed by Lehr alters the shape of oil from circle to ellipse with major axis in the direction of wind. Many researchers have improved the spreading of oil spills based on the Fay model, such as providing a surface area correction for the spreading formulas [10]. These regression models have been used to simulate oil spill scenarios in the marine simulator [11]. However, the regression models are difficult to simulate complex fluid-solid interaction scenes, such as the interaction scene between oil and booms.
Numerical models mainly refer to modeling oil spreading based on physical models, such as hydrodynamic methods [12, 13] and advection-diffusion equation [14, 15]. According to the different ways of model discretization, numerical methods can be divided into the Eulerian method [16, 17] and the Lagrangian method [14, 15]. The Eulerian method discretizes the problem domain into meshes, and the physical quantities at the vertices of the meshes change over time. The Eulerian method using volume fraction can easily trace the surface of the oil spill, but mesh generation for the problem domain is a prerequisite for the Eulerian method, and the accuracy is usually affected by the mesh resolution. The Lagrangian method is also called the meshless method, which discretizes the problem into a series of particles that carry physical quantities. The oil particles move according to the governing equations, and there is no need to generate meshes in the problem domain. The Lagrangian methods commonly used in oil spill simulation include the SPH (smoothed particle hydrodynamics) method [18–20] and the MPS (moving particle semi-implicit) method [21, 22]. The meshless method is suitable for simulating the spreading and dispersion of the oil, but it is difficult to track the free surface. Compared with regression models, numerical models have more advantages in simulating the effect of wind and waves on oil spreading and the interaction between oil and boom. On this basis, some researchers have proposed some hybrid methods that include both mesh and particles, such as the NBFLIP (narrow-band fluid implicit particle) method [23]. And, the NBFLIP method has been introduced into the marine simulator for oil spill simulation [24]. In recent years, with the development of artificial intelligence, some neural network models have gradually been used in prediction of oil transport, such as ANN (artificial neural network) [22]. Although these neural network models rely on the data generated by traditional models, the trained models have great application value in the oil spill scenarios of marine simulators.
2.2. Drift Models of Oil Spills
The drift of oil refers to the horizontal transport of oil spilled at sea under the action of ocean currents, waves, and winds. Hoult [25] derived the oil film drift equation, considering the effect of wind, waves, and tidal currents, which can be used to compute the drift velocity of any point on oil edge at different times and in different directions. Coppini et al. [26] adopted the Medsilk [27] drifting model for oil, and the spill is subdivided into a large number of Lagrangian particles. Coppini et al. considered the effect of wind field and turbulence on oil drift, and the information about oil slick features from remote-sensing observations was used to validate the oil drift model. The MIKE 21 SA model was used by De Padova et al. [28] to simulate the spatial evolution of floating oil. The Fokker–Planck equation for suspended oil particles is solved by the Lagrangian discrete parcel method. Dagestad et al. [29] developed an open-source Python-based framework for Lagrangian particle modeling. The framework is designed to be used for drift calculations in the ocean or atmosphere. OpenOil is a subclass of OpenDrift, and it is an oil drift model library including specific descriptions of processes relevant for oil spills. Based on the OpenOil, Röhrs et al. [30] simulated the advection of oil particles under the effect of currents, wind, and waves. It can be seen from the literatures that the Lagrangian particle method is one of the most widely used methods in oil spill drift modeling. With the development of technology and theory, some frameworks for simulating oil spill drift have been developed. However, most of the frameworks or models are used in oil spill forecasting, and there are relatively few applications in virtual reality oil spill simulation.
2.3. Viscosity Solvers for Meshless Frameworks
Viscosity is one of the important properties that affect the fate and transport of oil spills [31]. Considering the wide application of the meshless method in oil spill spreading and drift, we introduced a meshless DFPBF framework for unified modeling of oil spill spreading and drift. Therefore, in this section, we mainly investigate the viscosity solvers that can be used in meshless frameworks. Common viscosity solvers mainly include explicit viscosity solvers and implicit viscosity solvers. It is an intuitively explicit method to discretize the viscosity term by the Laplacian operator of the kernel function, but the second derivative changes the sign inside the kernel domain so that the normalization conditions are very poorly satisfied [32]. Müller et al. [33] designed a viscosity kernel function, which does not change the sign within the range of support. Brookshaw [34] derived the approximate form of the viscosity term based on the finite difference method and the gradient operator of the kernel function, which was later developed into a classic artificial viscosity by Monaghan [35, 36]. The artificial viscosity is Galilean invariant, vanishes for rigid body rotation, and conserves total linear and angular momenta. At the same time, Monaghan [36] proposed a concept of XSPH (“X” denotes the unknown factor). The velocity of the current particle is corrected by interpolating the velocity of the neighboring particles, which is a kind of visual viscosity. The XSPH method has also been used for viscosity calculations in some works [37, 38]. The explicit viscosity solvers have fast calculating speed and are only suitable for simulating fluids with low viscosity.
In order to simulate high-viscosity fluids, some implicit solvers have been proposed. Takahashi et al. [39] derived the first derivative of strain rate tensor to model the implicit viscosity. However, the complex neighborhood search of this implicit scheme decreases the performance. Peer et al. [40] proposed a decomposition form of velocity gradient. They simulated the viscous behavior by modifying the traceless shear rate tensor obtained by decomposition. Bender and Koschier [41] defined a velocity constraint function on the basis of the strain rate tensor. This method adjusts the target strain rate tensor by changing the viscosity parameter to simulate fluids with different viscosities. Although the above implicit viscosity solvers can stably simulate fluids with high viscosity, they are computationally expensive. Therefore, these solvers cannot be directly introduced into the oil spill simulation in the marine simulator for the time being. On the basis of artificial viscosity in [35, 36], Weiler et al. [42] proposed an implicit scheme that does not depend on the strain rate tensor, which is currently the fastest implicit viscosity solver.
3. DFPBF Framework
In the DFPBF framework, the following Navier–Stokes equations are usually used as the governing equations of the ocean waves and oil spills:where is the density, is the velocity, is the external body forces, and is the stress tensor. Equation (1) is the continuity equation, and equation (2) is the momentum equation of viscous flow. For incompressible and isotropic Newtonian fluid, the stress tensor can be defined aswhere is the pressure, is the identity tensor, is the dynamic viscosity, and is the strain rate tensor which is determined aswhere denotes the velocity gradient.
Inserting equations (3) and (4) into equation (2), we can obtain
According to the product rules of , and . For incompressible fluids, the velocity divergence should be zero, that is, . Therefore, the momentum equation can finally be expressed aswhere denotes the kinematic viscosity, is the pressure term, and is the viscosity term.
3.1. Divergence-free Velocity Constraint
The divergence-free velocity field is one of the important conditions to ensure incompressibility and also an important prerequisite of momentum equation (6). The divergence-free velocity constraint solver is a velocity-based solver for continuity equation (1). The DFPBF framework enforces divergence-free velocity field by the following velocity divergence constraint:where is the velocity divergence constraint, denotes the volume of neighboring particle , , is the index of the phase, and is the number of neighboring particles in phase .
The velocity divergence constraint aims to find a velocity correction that can map the velocity field to the divergence-free field, that is, . By applying Taylor expansion, we can obtainwherewhere is the Lagrange multiplier. Based on the concept of the gradient descent (GD) method, we enforce the velocity correction along the direction of the gradient of the constraint. According to equation (7), the gradient of the velocity constraint can be solved as
Similarly, according to Newton's third law, the gradient of the velocity constraint at the neighborhood particle can be defined as
Substituting equations (9)–(11) into equation (8), we can solvewhere and are equal in size and opposite in direction, andwhere is a term introduced to avoid singularities of the denominator.
The velocity correction of the particle under the divergence-free velocity constraints can be derived aswhere and are the rest densities of the fluid phase containing particle and .
Algorithm 1 is the calculation flow of the divergence-free velocity constraint solver (Algorithm 1). is the threshold of the density change rate, and iteriter is the number of iterations. In each iteration, firstly, the current velocity divergence constraint is computed according to equation (7). The Lagrange multiplier for each particle is then computed according to equation (12). We can obtain the velocity correction according to equation (14), and then update the particle velocity using until the average density change rate is less than the threshold , that is, the constraint of divergence-free velocity field is satisfied.
|
3.2. Constant Density Constraint
In the meshless particle frameworks, there are two main ways to solve the pressure term in the momentum equation: EOS-based (equation of state) solvers and PPE-based (pressure Poisson equation) solvers. However, EOS-based solvers usually show weakly compressibility, while PPE is computationally expensive. The DFPBF framework adopts a more controllable way to solve the pressure term, that is, the concept of position-based dynamics (PBD) is introduced into the meshless framework. The PBD solver we proposed works at the position layer and directly solves a position correction according to the nonlinear density constraint so that the density remains constant. It has better incompressibility than the commonly used PPE solvers, and after our optimization, its calculation speed has been significantly improved.
In our optimized PBD solver, we use EOS-like nonlinear density constraint:where is the constraint function of the position of particle and its neighbors and denotes the density of particle , and it is defined as
Equation (16) is an improved form of the particle number density model proposed by Solenthaler and Pajarola in [43]. Since we use a nonuniform sampling model in boundary handling, the volume of boundary particles will be adaptively adjusted with the shape of the rigid body, so we cannot directly compute the density of oil particles and ocean particles based on the particle number density [44]. Therefore, we modified the particle number density model and assumed that all particles in the neighborhood of have only the rest density of particle , as shown in Figure 1.
Density constraint solver aims to find a position correction for each particle, which is used to adjust the position distribution of fluid particles so that the density constraint . Applying Taylor expansion, we can obtainwhere the position correction for particle can be defined aswhere is the Lagrange multiplier for the density constraint solver. The particle position correction is also enforced along the direction of the constraint gradient . It can ensure momentum conservation of the system. The constraint gradient at is
The constraint gradient at is computed analogously to (11):
Substituting equations (18)–(20) into equation (14), we can obtainwhere is defined in equation (13). In principle, the parameter depends on the position distribution of the particles. When the position of the particle updates, the value of also needs to be updated according to the latest particle position. However, repeated experiments show that when simulating oil spill scenarios, the value of does not need to be updated in each iteration, which does not have a significant impact on the stability and convergence of the solver, but can significantly speed up the solution speed.
The position correction of the particle under the constant density constraints can be derived as
In fact, in the process of implementation, the density constraint usually uses inequality constraint, i.e., . This is because we have not considered the influence of air particles. Therefore, ocean particles, especially oil particles, at the free surface, usually have a lower density estimate due to the lack of neighboring particles. If the equality constraint in equation (15) is strictly implemented, it will cause particle clustering or clumping. The inequality constraint only adjusts the region where the particle density estimation is larger than rest density, which also improves the efficiency of the solver.
Algorithm 2 is the calculation flow of the constant density constraint solver (Algorithm 2). is the threshold of density change and iteriter is the number of iterations. Since we adopted a pressure projection scheme to decouple the relationship between pressure term and velocity in momentum equation (2), at the beginning of each iteration, we need to first compute the intermediate density based on the current particle position distribution. Then, we can compute the density constraint and Lagrangian multiplier according to equations (15) and (21), respectively. Finally, we can compute a position correction for each particle according to equation (22) and update the position of the particle using . Continue iterating until the average density change is less than the threshold .
|
When iteratively correct the particle position in Step 7, we introduce a relaxation schemewhere is a user-defined relaxation factor and . The value of can be adjusted within a reasonable range, but should not exceed 1.0. In principle, the larger the value of , the faster the algorithm will converge, but it will also lead to incompressibility and stability issues. Since the parameter is reused in this paper, a relatively smaller value should be selected to keep the algorithm stable. When the velocity is updated in the 8th step of the divergence-free velocity constraint solver shown in Algorithm 1, the relaxation factor can also be introduced to improve the convergence. However, experiments show that the velocity divergence solver only needs a few iterations to project the velocity field to a divergence-free field. Therefore, we only apply the relaxation scheme in the constant density constraint solver.
3.3. Implicit Viscosity Solver
Viscosity is one of the important properties that impact oil spreading and dispersion. Usually, low-viscosity oils flow more easily than those of higher viscosity. The viscosity of oils changes with their temperature and density. Generally, the viscosity of petroleum products is close to that of seawater, but the viscosity of crude oils is relatively high, and the viscosity will increase with the decrease of temperature. Therefore, the treatment of the viscosity term in the momentum equation (6) is very critical to the simulation of oil spill scenarios. At present, the most popular viscosity term solution schemes mainly include explicit schemes and implicit schemes.
3.3.1. Explicit Scheme
In the meshless method, we can directly discretize the viscosity term using particles:where denotes the Laplacian of kernel function. When computing the viscosity term, we only consider the contribution of the neighboring particles in the same phase. However, in the common Gaussian kernel functions, the Laplacian operator will change the sign in the support domain, and the value is generally large, which will cause some nonphysical phenomena.
In order to solve the above problem, the first-order derivative of the kernel function can be combined with the finite difference method to discretize the viscosity term:where is the spatial dimension, is the length of , and is introduced to avoid singularities. Equation (25) has the advantages of conservation of linear and angular momentum and is Galilean invariant. Since it is an explicit scheme, its calculation speed is faster. However, the stability of the explicit scheme is relatively poor and limited by the time step, and it is more suitable for simulating seawater with lower viscosity.
3.3.2. Implicit Scheme
In order to simulate higher viscosity oil spills, we introduced a stable implicit scheme. The second-order derivative of velocity can be obtained by computing the first-order derivative of the strain rate tensor in equation (4) [39]. Then, the first-order derivative can be discretized by gradient operator according the concept of DFPBF framework. This method can stably simulate high-viscosity oil spills and is not limited by the time step size. However, the calculation of two first-order derivatives will lead to a sharp increase in the workload of neighborhood search, which greatly affects the performance. In order to ensure stable simulation and calculation efficiency at the same time, we introduce an implicit scheme based directly on equation (25) instead of using the strain rate. In our DFPBF framework, the viscosity term is solved after enforcing the divergence-free velocity field constraint:where denotes the divergence-free velocity field. Equation (26) is an implicit integration scheme and can be expressed as the following linear system:where is the identity matrix and the matrix consists of 3×3 blocks:where is the average mass to ensure a symmetric and positive definite matrix of coefficients. In order to improve efficiency and save memory consumption, the matrix-free conjugate gradient method is used to solve the linear system. In order to reduce the condition number of the coefficient matrix and improve the convergence, we adopt a preconditioning method. The block version of Jacobi preconditioner consists of the 3 × 3 blocks .
When we couple the implicit viscosity solver into our DFPBF framework, we find that although implicit viscosity solver performs well in the simulation of high-viscosity fluid, it becomes unstable when the fluid particle velocity is high. In fact, under the same initial conditions, the velocity of low-viscosity fluid particles is higher than that of high-viscosity fluid particles. Moreover, the viscosity coefficient in the implicit viscosity solver is not consistent with the actual viscosity of fluid. In meshless frameworks, one of the purposes of introducing artificial viscosity is for the convergence and stability of numerical computation, and there is no strict correspondence with the actual physical meaning. For example, in the implicit viscosity solver, the motion of a fluid with a viscosity coefficient of is almost the same as a fluid with an actual viscosity of 1.0 × 10−3 Pa s. In order to uniformly simulate fluids of different viscosities, we introduce a viscosity correction coefficient:where is a user-defined correction coefficient, and in this paper.
3.4. Adaptive Time Step Size
In the simulation process, in order to speed up the physical calculation, it is usually preferred to choose a larger time step, but considering the stability, the time step size must also be small enough. Based on the above requirements, in meshless methods, the time step size can be estimated adaptively based on the CFL (Courant–Friedrichs–Lewy) condition. The adaptive scheme can maximize the time step while maintaining the stable convergence of the algorithm. Specifically, the time step should satisfy the following condition:where is the diameter of particle, is the maximum velocity of particle, and is a scaling factor, . When , the time step constraint in equation (29) means that the displacement of all particles in each time step should not exceed the particle diameter. On this basis, in order to prevent some particularly high-speed particles from affecting the performance of the algorithm, we specify the maximum and minimum time steps globally:where and denote the minimum and maximum thresholds of the time step, respectively.
Algorithm 3 shows the calculation flow of the DFPBF framework (Algorithm 3). We adopt a splitting scheme to decouple the solution of pressure, viscosity, and external forces. In each time step, we first compute the external forces and predict an intermediate velocity and intermediate position according to the external forces. After that, the adaptive time step is computed to keep the simulation stable. Since the particle position has been updated, it is necessary to perform the neighborhood search again in Step 11. Since we have optimized the simulation loop, the parameter is computed only before running the constraint solver. Adjust the particle distribution according to the constant density constraint solver in Algorithm 2. After applying the density constraint, a new position is obtained, and then, the intermediate velocity can be updated according to . We can get a divergence-free velocity field after applying the divergence-free velocity constraint solver in Algorithm 1. After obtaining the divergence-free velocity field, the viscosity solver is used to solve the viscosity term in the momentum equation (6).
|
4. Results and Discussion
In this section, we simulated some oil spill scenarios based on the oil spill model proposed in this paper. We use C++ to implement the physical calculation part of the oil spill model in Visual Studio 2019 and use OpenGL to visualize the oil spill scenarios. Our project runs on a workstation with an AMD Ryzen 7 3700X CPU (8-Core) and 32 GB RAM. The operating system is windows 10(x64). We adopt the adaptive single layer nonuniform boundary handling method, in which the velocity of static rigid boundary particles is 0.
Figure 2 is the schematic diagram of oil spill test scenario. In the experiment, the volume of ocean particles and oil particles are equal, and the particle radius is 0.015 L. The size of the bounding box is 4.0 L×4.0 L×2.0 L and the number of boundary particle is 175.1k. The size of ocean block is 4.0 L×0.3 L×2.0 L and the number of ocean particle is 78.5k. The size of oil block is 0.5 L×1.0 L×0.5 L and the number of oil particle is 8.2k. The oil block is located above the ocean block and gradually interacts with the ocean under the action of gravity.
(a)
(b)
We conducted the first set of experiments on calm ocean. Figure 3 shows the simulation results. We simulated the movement of three oil spills with different viscosities under the same initial conditions (2.6 × 10−3 Pa s, 9.0 × 10−3 Pa s, and 1.0 × 10−1 Pa s). The density of oil (dark particles) is 840 kg/m3, and the density of ocean (blue particles) is 1027 kg/m3. It can be seen from the simulation results that the low viscosity oil has a larger spreading area and is easily broken into smaller oil films. The simulation results are basically consistent with Fay's assumption [6] that the oil film spreads in a circular shape on calm ocean. Compared with traditional methods, our DFPBF framework can simulate the natural breakup of the oil film.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
We added a stochastic fluctuating wind field to the scene, as shown in Figure 3. The direction of the wind field is the negative direction along the x-axis. It can be seen from Figure 4 that the oil film drifts with the waves. Under the same wind wave conditions, the low viscosity oil spills are more affected by the wind waves. The shape of low-viscosity oil spreading on the ocean is more like an ellipse, with the long axis along the wind direction, which is consistent with Lehr's conclusion [8, 9]. In contrast, high-viscosity oil spills are more resistant to shear stress, so they are less affected by wind waves than low-viscosity oil spills.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Figure 5 is the coupling scene between oil spill and rock. The spilled oil flows out of the pipeline and the number of oil particles is 20k. The dynamic viscosity of oil in Figure 5 is 9.0 × 10−3 Pa s. Due to the viscosity, the spilled oil will stick to the rocks during the fall. When the oil spills drift to the nearshore zone, the animals, plants, and reefs may be coated by the dark crude oil with high viscosity. Our DFPBF framework with the implicit viscosity solver can realistically simulate these scenes, as shown in Figure 5.
(a)
(b)
(c)
(d)
We name the DFPBF framework which is using the relaxation scheme in Section 3.2 as rDFPBF (relaxed DFPBF). Figure 6 shows the running time statistics of the oil spill scenarios shown in Figure 2. The time step size in this experiment is 0.01 s. It can be seen from Figure 6 that the rDFPBF framework shows attractive acceleration performance in both test scenarios. In Figure 6(a), the average time cost per time step of the DFPBF framework is 5.18 s and that of the rDFPBF framework is 4.38 s. In Figure 6(b), the average time cost per time step of DFPBF framework is 5.39 s and that of the rDFPBF framework is 4.59 s. As the viscosity of oil spills increases, the time cost also increases.
(a)
(b)
5. Conclusions
In this article, we proposed a meshless DFPBF oil spill simulation framework. The problem domain is discrete by Lagrangian particles, which is very suitable for simulating the breaking and merging motion of the oil film. The particle number density model was used to calculate the physical quantities of Lagrangian particles, which can accurately estimate the physical quantities and handle density discontinuities at the interface between multiple fluids with varying rest densities. We also introduced an implicit viscosity solver for our DFPBF framework, which can stably simulate oil spills with various viscosities. We conducted a series of simulation experiments based on our DFPBF framework, and the following conclusions can be drawn based on the experimental results:(1)On a calm ocean, the oil spills will spread approximately into a circular shape, and as the viscosity increases, the spreading area gradually decreases. In this respect, our conclusion is consistent with Fay's conclusion [5, 6].(2)Under wind wave conditions, the oil spill will drift with the wind waves, and the final position will deviate from the initial oil spill position. The spreading shape of the oil spill is more like an ellipse, and the long axis of the ellipse is in the direction of the wind waves, which is consistent with Lehr’s conclusion [8, 9].(3)Our oil spill model can simulate the scene of oil spills covering reefs, which helps improve the trainees’ emergency response capability to deal with such oil spills.
We have made efforts to improve the performance of the DFPBF framework, such as adaptive time step size, matrix-free conjugate gradient method, and Jacobi preconditioner. When applied to marine simulators, parallel computing and distributed computing can be used to further improve the performance. Since our DFPBF framework is only designed for virtual scenes in marine simulators and only used to train trainees' emergency response ability, we are more inclined to improve its performance than the accuracy of the algorithm. Therefore, our DFPBF framework cannot be directly introduced to oil spill forecasting. If researchers try to extend our method in oil spill forecasting, it is necessary to compare the simulated results with the measured data. Moreover, some parts of our DFPBF framework need to adopt higher-precision scheme, such as changing Euler integration to higher-order Runge–Kutta and changing the kernel function from cubic spline to a more precise form.
In future work, we will try to introduce surface tension [45] and turbulence effects [46] to our DFPBF framework in order to simulate oil spill spreading and drift transport with Langmuir circulation. In terms of interactive simulation of boom and oil spill, we try to model the boom based on the PBD method to further enhance the realism.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
The authors thank the financial support provided by National Natural Science Foundation of China (Grant no. 52071312), Liaoning Provincial Natural Science Foundation shipping joint fund project (Grant no. 2020-HYLH-29), and Provincial postgraduate education and teaching achievement award cultivation project (Grant no. YJG2020601).