Modelling and Simulation in Engineering

Volume 2012 (2012), Article ID 398387, 7 pages

doi:10.1155/2012/398387

Research Article

Eulerian Oil Spills Model Using Finite-Volume Method with Moving Boundary and Wet-Dry Fronts

Ehsan Sarhadi Zadeh¹ and Kourosh Hejazi²

¹ Hydraulics Group, Department of Civil Engineering, K. N. Toosi University of Technology, Tehran 1996715433, Iran

²Department of Civil Engineering, K. N. Toosi University of Technology, Tehran 1996715433, Iran

Received 3 April 2012; Revised 27 July 2012; Accepted 5 August 2012

Academic Editor: Farouk Yalaoui

Copyright © 2012 Ehsan Sarhadi Zadeh and Kourosh Hejazi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The world production of crude oil is about 3 billion tons per year. The overall objective of the model in present study is supporting the decision makers in planning and conducting preventive and emergency interventions. The conservative equation for the slick dynamics was derived from layer-averaged Navier-Stokes (LNS) equations, averaged over the slick thickness. Eulerian approach is applied across the model, based on nonlinear shallow water Reynolds-averaged Navier-Stokes (RANS) equations. Depth-integrated standard k-ε turbulence schemes have been included in the model. Wetting and drying fronts of intertidal zone and moving boundary are treated within the numerical model. A highly accurate algorithm based on a fourth-degree accurate shape function has been used through an alternating-direction implicit (ADI) scheme which separates the operators into locally one-dimensional (LOD) components. The solution has been achieved by the application of KPENTA algorithm for the set of the flow equations which constitutes a pentadiagonal matrix. Hydrodynamic model was validated for a channel with a sudden expansion in width. For validation of oil spill model, predicted results are compared with experimental data from a physical modeling of oil spill in a laboratory wave basin under controlled conditions.

1. Introduction

Oil transport, exploration, and storage facilities are all possible sources of spills. The 1991 Persian Gulf (Kuwait War) oil spill was estimated at 143 million liters or 38 million gallons [1]. The fate and transport of spilled oil is governed by the advection due to current, wave, and wind; horizontal spreading of the surface slick due to turbulent diffusion, gravitational, inertia, viscous, and surface tension forces; emulsification; weathering processes such as evaporation, dispersion, and dissolution; interaction of oil with shoreline; photochemical reaction and biodegradation.

The chemical and biological processes generally occur a long time after the oil spill [1]. Tkalich [2] applied a consistent Eulerian model; the slick thickness is computed using layer-averaged Navier-Stokes (LNS) equations, and the advection-diffusion equation is employed to simulate oil dynamics in the water column. A high-order accuracy numerical scheme is developed, to match the observed balance between advection, diffusion, and spreading phenomena.

Wang et al. [3] showed that the amount of oil released at sea is distributed among a large number of particles tracked individually. Oil particles are driven by a combination of current-, wave- and wind-induced speed and move in a 3D space. Horizontal and vertical diffusion are taken into account using a random walk technique. Periáñez [4] developed a numerical model that simulated the dispersion of contaminants, including oil spills, in the Sea.

The hydrodynamic equations were solved in advance and include a barotropic model for calculating tides and a reduced gravity model for the average circulation. Dispersion was calculated using particle-tracking methods. Turbulent diffusion and specific processes for contaminants (for instance, decay, biodegradation, or oil evaporation) were simulated by Monte Carlo techniques.

Guo et al. [5] developed a hybrid particle tracking by an Eulerian-Lagrangian approach for the simulation of spilled oil in coastal areas. To acquire accurate environment information, the model was coupled with the 3D free surface hydrodynamics model, Princeton ocean model (POM), and the third generation wave model, simulated waves Nearshore (SWAN). By simulating the oil processes, it has the ability to predict the horizontal movement of surface oil slick, vertical distribution of oil particles, oil concentration in the water column, and mass balance of spilled oil.

2. Oil Spill Processes

The fate and transport of oil spilled in water is dominated by complex physicochemical processes that depend on oil properties and hydrodynamic and environmental conditions [6]. Figure 1 is a schematic diagram showing these processes. In this section, some of the available algorithms describing physical and chemical weathering processes are described.

Figure 1: Transport and fate of oil slick in seas.

2.1. Spreading on Water Surface

The spreading of an oil slick would pass through three distinct phases. In the beginning phase, only gravity and inertia forces are important. In the intermediate phase, the gravity and viscous forces dominate. The final phase is governed by the balance between surface tension and viscous forces [1].

Lehr et al. [7] proposed an elongated ellipse model along the direction of the wind to account for the observed nonsymmetrical spreading of oil slicks

where and are the lengths of the minor and major axes, respectively (m); is the slick area (m²); : and are densities of water and oil, respectively; is total volume of the spilled oil (barrel); is wind speed (knot); is time (min).

2.2. Evaporation

Evaporation causes a very significant mass loss in many kinds of oil and has a profound effect on density, viscosity, and other properties of oil. In the current work, the analytical method proposed by Mackay [6] has been adopted:

where is the volume fraction evaporated; is evaporation exposure; is mass transfer coefficient for evaporation (m/s); is the slick area (m²); is time (min); is initial volume of spilled oil (m³); is wind speed (m/s); is oil slick diameter (m); is the Schmidt number which represents the surface roughness and is taken as 2.7; is environmental temperature (Kelvin = °C + 273); , , , and are derived constants from distillation data [8]. initial boiling point at equal to zero (K); are gradient of the boiling point (K); environmental temperature (K); , constants taken 6.3 and 10.3 for undefined crudes, respectively (Table 1).

Table 1: Data from the Bobra evaporation experiments.

Oil
(K)
(K)

Density (g/mL) 15°C
Viscosity (cp) 15°C

Adgo
551.0
195.4
24.278
20.535
0.9530
61.6

ASMB
397.0
539.1
8.133
11.594
0.8390
9.2

Amauligak
471.0
370.0
12.074
14.485
0.8896
14.0

Bent Horn
406.0
484.0
11.408
13.834
0.8181
24.0

Diesel
517.0
139.8
20.274
18.052
0.8270
2.7

Endicott
453.7
1400.0
−0.815
7.107
0.9149
84.0

North Slope
430.6
722.0
4.452
10.137
0.8936
23.0

Panuke
268.1
367.8
7.174
11.437
0.7757
1.1

2.3. Dissolution

The most soluble oil components are usually the most toxic. In the present study, the mass of soluble is negligible compared to the dispersed oil droplet near the surface but of the same order of magnitude in the deeper water.

Mackay [6] developed a multicomponent theory to compute the rate of oil dissolution

where is the total dissolution rate of the oil slick (gr/s), is dissolution mass transfer coefficient ( m/s), is the slick area (m²), and is the oil solubility in water.

Huang et al. [9] suggested the solubility for typical oil:

where is the solubility for fresh crude oil; is decay constant; is time (hr).

2.4. Emulsification

Emulsions formation changes the properties and characteristics of oil drastically. Mackay [6] proposed the following expression for the incorporation rate of water into an oil slick:

where is the fractional water content; is curve fitting constant that varies with wind speed (); is mousse viscosity constant (0.7 for crude oils and heavy fuel oil, and 0.25 for home heating oil); is wind speed (m/s).

2.5. Weathering Module

The behavior of oil spill depends not only on the prevailing conditions but also on the physicochemical properties of the oil itself. While the former are site and time dependent, the latter change while interacting with each other with oil movements.

The evaporation process, together with dissolution and the mousse formation, leads to an increase in volume and density, respectively [5]:

where is the oil volume; is initial volume of spilled oil (m³); is the fraction dissolved; is the density of the remaining oil.

Evaporation and emulsification also cause an increase in viscosity [5] in the subsequent equation:

where is the Mousse viscosity; is oil-dependent constant; is parent oil viscosity; is percentage Asphaltene.

3. Numerical Modeling

Computations proceed with an Eulerian model on the oil dynamics equations that describe slick advection spreading on water surface [2] which are given as

where is the oil slick thickness; is slick drifting velocity; is shear stress due to wind; is slick spreading-diffusion function; is “oil film-water surface” friction coefficient (0.02 kg/m²s); is physical-chemical kinetic terms; is gravity acceleration; .

The oil spill model couples with a 2D hydrodynamic model. A staggered grid system is used to solve the equations [10]. The oil slick thickness is placed at the center of the grid, while velocity components and are placed at the midpoints of the interfaces (Figure 2).

Figure 2: 2D computational staggered grid.

A 4th-degree shape function [11], for computing h, u, and v, has been adopted (Figure 3)

Figure 3: 4th-degree shape function.

3.1. Governing Equations

In matrix form, a conservation law of the 2D nonlinear shallow water equations and continuity equation, including the effects of bed friction stresses, vegetation shear stresses, wind shear, viscous terms, and the earth’s rotation (Coriolis effects) [12], may be written as

where denotes time, and are Cartesian coordinates, and are vectors representing the conserved variables, fluxes in - and y-direction, and source terms, respectively; represents free surface elevation above , which is still water depth, and so the total water depth equals ; and are depth-averaged velocity components in the and directions, respectively; is Coriolis parameter due to earth’s rotation ; is angular rotational speed of earth; is geographical angle of latitude; are viscous shear stress; , bed friction stresses; are vegetation friction stresses; are wind friction stresses.

For simulating of the flooding-drying process over the coastal regions and in flooding problems, a threshold value is needed, which has been set to 1 mm for all of the computations in the current study. For a dry grid point, the velocity components are set to zero to avoid numerical error, while the water level remains unchanged. In this method, the numerical computational domain covers the entire intertidal area. The moving boundary between the land and water (where water flux is equal to zero) is reconfigured for each time step.

The alternative direction implicit (ADI) scheme has been deployed for solving depth-integrated Navier-Stokes equations [12]. ADI technique involves the subdivision of each time step into two half time steps; therefore, a 2D implicit scheme can be applied by considering only one dimension implicitly for each half time step (Figure 4).

Figure 4: Solving the N-S equations by sweeping.

In the 1st half time step, the water elevation, velocity component, and oil slick thickness are solved implicitly in the x-direction, whilst the other variables are represented explicitly. Similarly in the 2nd half time step, the water elevation, velocity component, and oil slick thickness are solved implicitly in the y-direction, whilst the other variables are represented explicitly.

The equations of continuity and momentum in x-direction for the 1st half time step, for example, are coupled in the form

where and are discharges per unit width; and are known quantities.

By sweeping the coupled equations for each row and column of the whole domain, the set of the equations for flow which constitutes a pentadiagonal matrix is obtained. The solution has been achieved by the application of KPENTA algorithm [13]

4. Model Validation and Verification

4.1. Hydrodynamic Model

To verify the hydrodynamic model, a uniform flow in a straight open channel has been simulated. Rodi [14] reported measured values of longitudinal velocities across a rectangular channel cross-section, having a width-to-depth ratio of 30 and a Manning roughness factor of .

Figure 5 shows comparisons between the computed and measured values. It can be observed to agree well with laboratory measurements for the boundary-layer region near the bank and shows an excellent agreement within the domain. The measured square of the Pearson product moment correlation coefficient for 4th- and 1st-order numerical model is and , respectively.

Figure 5: Horizontal velocity distribution.

4.2. Oil Spill Model

For validation of oil spill model, predicted results are compared with experimental data from a physical model of an oil spill in a laboratory wave basin [15] (Figure 6).

Figure 6: Model bathymetry.

Regular waves were produced at one end of the basin by a piston-type paddle. These waves were generated with a period of 1.29 s and offshore deep water wave height of 0.098 m. At the other end of the basin, the wave energy was dissipated on an upper plane beach.

Figure 7 shows the measured and predicted steady state wave height as a function of distance from the origin. In general, there is reasonable correlation between the measurements and numerical predictions , although the numerical simulation has smoothed out some of the observed behavior.

Figure 7: Offshore wave height distribution.

Figure 8 shows the temporal growth in area of the oil slick with thickness greater than 6 μm. It is clear that the areal dimensions of the physical slick increase more rapidly than predicted values by the numerical model, .

Figure 8: Growth of oil slick area.

Figures 9 and 10 show the movement of the slick over 40 second (30 wave periods) after the initial spill in numerical and experimental model, respectively. Allowing for the disparity between the centers of the nearshore circulation cells, it may be judged that the two advection paths are in reasonable agreement. The large differences between the physical and numerical slick shapes underline the importance of scale effects, for example, oil density, surface tension, turbulent eddy sizes, and so on. Besides, these differences may be explained partly by the discrepancies in modeling wave current interaction. To achieve better results in the small oil spills, Fay’s equations for surface tension-viscous spreading phase should be considered which are not represented here [16].

(a)

(b)

(c)

(d)

Figure 9: Development of 6 μm oil thickness contours in numerical model in 4 steps.

Figure 10: Development of 6 μm oil thickness contours in experimental model.

5. Conclusion

A 2DH numerical model has been developed based on nonlinear shallow water Reynolds-averaged Navier-Stokes (RANS) equations. The model deploys an Eulerian domain description, with a staggered grid system. Classical empirical formulations for weathering processes of oil slick such as spreading, evaporation, dissolution, and emulsification have been included. A 4th-order alternating direction implicit (ADI) finite volume scheme has been used to reduce the truncation error. The set of the equations which constitute a pentadiagonal matrix has been solved by the application of KPENTA fast computing algorithm which greatly reduces the computation time. In similar models, many simplification techniques have been used that lead to a tridiagonal matrix solved by the Gaussian elimination time-consuming method which is not suitable for operational conditions.

Results from the numerical model are compared with experimental data. Reasonable agreement is obtained between predicted and experimental nearshore circulation patterns. Discrepancies occurred are due mainly to reflections from the smooth impermeable beach area. These may reduce for sandy or pebble beach areas. It is recommended that the wave ray formulation be used by a more comprehensive approach, such as the one proposed by Copeland [17]. The major discrepancies in oil slick modeling were caused by underestimation of the wave height and the no inclusion of surface tension-viscous spreading.

References

ASCE Task Committee, “State-of-the-art review of modeling transport and fate of oil spills,” Journal of Hydraulic Engineering, vol. 122, no. 11, pp. 594–609, 1996.

P. Tkalich, “A CFD solution of oil spill problems,” Environmental Modelling and Software, vol. 21, no. 2, pp. 271–282, 2006.

S. D. Wang, Y. M. Shen, Y. K. Guo, and J. Tang, “Three-dimensional numerical simulation for transport of oil spills in seas,” Ocean Engineering, vol. 35, no. 5-6, pp. 503–510, 2008.

R. Periáñez, “Chemical and oil spill rapid response modelling in the Strait of Gibraltar-Alborán Sea,” Ecological Modelling, vol. 207, no. 2-4, pp. 210–222, 2007.

W. J. Guo and Y. X. Wang, “A numerical oil spill model based on a hybrid method,” Marine Pollution Bulletin, vol. 58, no. 5, pp. 726–734, 2009.

D. Mackay, Oil Spill Processes and Models, Environmental Protection Service, Canada, 1980.

W. J. Lehr, H. M. Cekirge, R. J. Fraga, and M. S. Belen, “Empirical studies of the spreading of oil spills,” Oil and Petrochemical Pollution, vol. 2, no. 1, pp. 7–11, 1984.

M. Bobra, “A study of the evaporation of petroleum oils,” Environment Canada EE-135, Ontario, Canada, 1992.

J. C. Huang and F. C. Monastero, Review of the State-of-the-Art of Oil Spill Simulation Models, American Petroleum Institute, 1982.

S. V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, NY, USA, 1980.

J. H. Ferziger and M. Peric, Computational Method for Fluid Dynamic, Springer, New York, NY, USA, 3rd edition, 2002.

Hirsch and C., Numerical Computation of Internal and External Flows, Elsevier, Oxford, UK, 2nd edition, 2007.

A. A. Karawia, “A computational algorithm for solving periodic penta-diagonal linear systems,” Applied Mathematics and Computation, vol. 174, no. 1, pp. 613–618, 2006.

W. Rodi, “Turbulence models and their application in hydraulics,” 1984.

A. G. L. Borthwick and S. A. Joynes, “Laboratory study of oil slick subjected to nearshore circulation,” Journal of Environmental Engineering, vol. 118, no. 6, pp. 905–922, 1992.

J. A. Fay, “The spread of oil slick on a calm sea,” in Oil on the Sea, D. P. Hoult, Ed., pp. 53–63, Plenum Press, 1969.

G. J. M. Copeland, “Practical radiation stress calculations connected with equations of wave propagation,” Coastal Engineering, vol. 9, no. 3, pp. 195–219, 1985.