A Numerical Model for the Confinement of Oil Spill with Floating Booms

Spill Science & Technology Bulletin, Vol. 7, Nos. 5–6, pp. 249–255, 2002 Ó 2002 Elsevier Science Ltd. All rights reserved Printed in Great Britain 1353-2561/02 $ - see front matter

PII: S1353-2561(02)00042-7

A Numerical Model for the Confinement of Oil Spill with Floating Booms SONG-PING ZHU * & DMITRY STRUNINà  Mechanical Engineering Department, University of Northern Ohio, Ada, OH 45810, USA àDepartment of Mathematics and Computing, University of Southern Queensland, Toowoomba, QLD 4350, Australia

Zhu and Strunin (Appl. Math. Model., 2001) presented a simple mathematical model for the confinement of oil spills, through the use of floating barriers such as booms, in an open ocean. In that paper, solutions with an indirect approach to solve the nonlinear integral equation system were discussed. In this paper, we shall present some recent results of the model obtained with a direct approach. This direct solution procedure enables us to quantitatively discuss the relationship among three important physical parameters, the Froude number, F, the barrier submergence depth, l, and the amount of oil trapped by the barrier, q. Although a further upgrading of our model to somehow incorporate the shear stress on the oil–water interface appears to be necessary, the complexity involved in such upgrading warrants the current model as an initial step in our modelling exercises of this complicated flow phenomenon. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Modelling oil spill booms, confinement, Froude number, submergence, shear stress

Introduction Oil spills pose serious threats to our environmental and ecological system. As our society becomes more and more dependent on petroleum and petroleum products (e.g., gasoline, diesel, etc.), the demand of using supertankers to ship crude oils is undoubtedly going to keep increasing and consequently disasters of oil spills off the surface of open oceans or bays are unavoidable. There have been many oil spill incidents from supertankers in the past ten years already. Examples of these type of accidents range from the most recent one taking place in Tokyo Bay with about 10,000 bbl of crude oils spilled on July 2, 1997 (The Oil

*Corresponding author. Tel.: +1-419-772-2385; fax: +1-419-7722404. E-mail address: [email protected] (S.-P. Zhu).

and Gas Journal, 1997) to the most notorious Exxon Valdez tanker spill off AlaskaÕs Prince William Sound on May 24, 1989 (Times, 1989). As demanded by the devastating nature of oil spills, they have been studied quite intensively in the past decade. Majority of these studies are dedicated to modelling the spread of spilled oil on the surface of an open ocean, understanding the basic physico-chemical mechanisms of the spreading and developing hindcast models for some specific accidents. In the earlier seventies, Hoult (1972) undertook a detailed study of the oil spreading mechanisms and outlined three stages of the spreading––inertial, viscous and surface-tensiondriven. For each stage, a relevant theory was proposed and a comparison with some available experimental data were presented. Later on, further work was done in this direction, for example, by Sebastiao and Soares (1995) who took into account evaporation, emulsification and change in viscosity and density of the oil. 249

S.-P. ZHU & D. STRUNIN

Empirical relations were extensively used, which gave good correlation with field observations for the problem with such a high degree of complexity. However, models heavily depending on empirical formulae generally tend to be less versatile. Reed and Gundlach (1989a,b) developed a numerical oil spill model with a number of those factors taken into account and tested it against data on the Amoco Cadiz event. The data from this accident were also used by Cuesta et al. (1990) to test their numerical model. Based on a twodimensional depth-averaged barotropic model, Noye et al. (1993) presented a nice numerical model for the oil slick movement in the Northern Spencer Gulf of South Australia. Another big accident, the grounding of the tanker Braer, also stimulated considerable modelling efforts, although the results were not quite satisfactory. Turrell (1994) and Proctor et al. (1994) analysed these results and indicated possible ways of improvement of the models. Relatively, there are very few papers in the literature that address the confinement of spilled oil with floating booms, let alone the modelling of the flow near a floating boom. A very practical approach of confining the oil spill on the surface of an open ocean is to use mechanical barriers, such as floating booms. However, the confined oil may leak beneath the floating boom if either the towing speed of the boom or the amount of the confined oil or both exceeds certain critical values. The lack of quantitative understanding of the relationship between the submergence depth of the boom and these critical values requires some mathematical modelling of this practical engineering problem. To tackle this problem of considerable practical importance, Strunin and Zhu (1998) first conducted some simple experiments to observe the flow characters qualitatively. Then, Zhu and Strunin (2001) made an attempt towards building a mathematical model to simulate the flow near a floating barrier. In their model a potential flow was considered with a layer of heavy fluid (water) flowing underneath a finite volume of light fluid (oil) floating near but trapped by a vertical barrier. To find the solution, a similar approach to that of Zhang and Zhu (1996) was adopted to solve a set of nonlinear integral equations numerically. Critical Froude numbers were defined as those, corresponding to which convergent solutions no longer exist. Therefore, when the actual Froude numbers are larger than these critical Froude numbers, oil starts to leak underneath the barrier. It was shown that the oil leakage under the barrier is impossible if the contact point is a stagnation one. For other (nonstagnation) cases, the flow were computed up to a critical situation where the oil starts leaking underneath the barrier. However, since only an indirect approach was adopted to solve the nonlinear integral equations for the purpose of numerical efficiency, the water velocity at the 250

water–oil upstream contact point became (along with the Froude number) one of input parameters. It was not clear how the interface positions vary individually with each one of the three physical parameters (Froude number, submergence depth of the barrier and the total mass of the trapped oil). In the present paper we present some new results on this model with a direct approach, which we have recently incorporated into the model. Some discussions on these newly obtained results lead us to believe that upgrading the model by taking into account the viscous effect, at least that within the boundary layer near the oil–water interface should be our next step of a series of attempts to eventually simulate and understand this quite complicated flow phenomenon with the aid of modern computers. The section ‘‘Model formulation’’ contains the necessary description of the model formulation. In the section ‘‘Numerical implementation and result discussion’’ we present the numerical implementation of the new direct solution approach and make some discussions on some very preliminary numerical results on the full model. Our concluding remarks are given in the last section.

Model Formulation As described in details in Zhu and Strunin (2001), the model is characterised by the assumption that the viscosity in the water is so small that it can be neglected and the water flow is irrotational. As a result of this assumption, potential theory is adopted to describe the water flow and the normal pressures across the interface are matched to determine the unknown interfacial position. An important dimensionless parameter in the model is the Froude number, which is defined as pffiffiffiffiffiffiffi F ¼ V1 = gH ; where V1 is the upstream velocity, H the upstream water depth and g the gravitational acceleration. From now on, all the quantities are nondimensional unless stated otherwise. Although the details of the model formation has been given in Zhu and Strunin (2001), we shall still describe the key equations in the model for the completeness of this paper. We first introduce a complex potential f ¼ / þ iw with / being a potential function and w being a stream function. The complex potential f is an analytical function of the complex variable z ¼ x þ iy with x and y being the horizontal and vertical coordinates of the coordinate system shown in Fig. 1. The total discharge can be nondimensionalised to unity, therefore Spill Science & Technology Bulletin 7(5–6)

MODELLING THE CONFINEMENT OF OIL SPILL WITH FLOATING BOOMS

apply the CauchyÕs formula to the upper half of the fplane: Z 1 1 XðnÞ dn; Imf > 0; ð3Þ XðfÞ ¼ 2pi 1 n  f where the integral is understood to be in the sense of the Cauchy principal value. For real f, Eq. (3) becomes Z 1 1 XðnÞ XðfÞ ¼ dn; Imf ¼ 0; ð4Þ pi 1 n  f which can be separated into real and imaginary parts Z Z 1 1 dðnÞ 1 1 sðnÞ sðfÞ ¼  dn; dðfÞ ¼ dn: p 1 n  f p 1 n  f ð5Þ

Fig. 1 Physical plane and the coordinator system.

the bottom and top boundaries are two streamlines with w ¼ 1 and 0 respectively. The point of zero potential is chosen to be the lowest point E of the barrier (see Fig. 1). After applying the conformal mapping f ¼ epf

ð1Þ

the region occupied by the water is transformed to the upper half of the f-plane as shown in Fig. 2. The bottom and free surfaces are respectively mapped onto the negative and positive real axes. It is convenient to introduce the logarithmic velocity variable X  d þ is defined by u  iv 

On the interface and the free surface, the Bernoulli equation must be satisfied. After a nondimensionisation based on g, V1 , H and the density of the heavy fluid q1 (the density of the light fluid is q2 ), the Bernoulli equation is, in terms of the new variable f, of the form

df ð zðfÞÞ ¼ e½iXðfÞ : dz

ð2Þ

It can be shown that d represents the direction of the velocity and es the magnitude of the velocity. Assuming X ! 0 with jfj ! 1 (physically, this corresponds to the situation where the flow far upstream is parallel to the bottom and velocity is uniform) we can

P ðfÞ þ

 F 2  2sðfÞ  1 þ yðfÞ ¼ 0; e 2

ð6Þ

where P ðfÞ and yðfÞ are respectively the dimensionless pressure on the interface and the interface elevation. The corresponding dimensional variables are P ¼ P = ðq1 gH Þ; y ¼ y =H : With the assumption that the oil is motionless, the pressure in it must be hydrostatic, that is, it has no horizontal gradient and linearly varies vertically. The continuity of the pressure on the interface demands that     P y þ Dy ¼ P y  q2 gDy ;

or

DP ¼ q2 g; Dy

which can be written in nondimensional form as 

dP q2 ¼  a: dy q1

ð7Þ

Now, differentiating (6) with respect to f yields dP dy ds ¼   F 2 e2s : df df df

ð8Þ

The relation (7) can be rewritten as dP df ¼ a: df dy

ð9Þ

Eliminating dP =df from (8) and (9) leads to F 2 e2s Fig. 2 The f-plane. Spill Science & Technology Bulletin 7(5–6)

ds df ¼ a  1; df dy

fA < f < f B :

ð10Þ

251

S.-P. ZHU & D. STRUNIN

It follows from (2) that on the real axis of the fplane dx es ¼ cos d; df pf

dy es ¼ sin d; df pf

ð11Þ

which are kinematic conditions on the free surface and interface. Taking into account that dy=df is represented by the second equation in (11), we integrate (10) to express s in terms of d on the interface:  1 3ð1  aÞ sðfÞ ¼ ln expð3sB Þ  3 pF 2 Z fB sin dðnÞ dn ; fA < f < fB ; ð12Þ n f where sB is the value of s at the point B. On the free surface, instead of (10), the equation F 2 e2s

ds df ¼ 1; df dy

0 < f < 1;

ð13Þ

which is a particular case of (10) with a ¼ 0, must be satisfied. Integrating (13) yields  1 sðfÞ ¼ ln expð3sD Þ 3 Z f 3 sin dðnÞ þ 2 dn ; 0 < f < 1; ð14Þ pF 0 n where sD is the s value at the point D, which denotes far downstream. On the other hand, the first equation in (5) can be written as



Z 1 1 dðnÞ p

fA  f

sðfÞ ¼  dn  ln

p 2 f1

0 nf Z fB dðnÞ þ dn ; f > 0 ð15Þ fA n  f in view of the boundary conditions: d  p=2 for 1 < f < fA (vertical barrier) and d  0 for f > fB (rigid lid) and for f < 0 (bottom). Equating (15) and (12), we obtain an integral equation for the unknown function dðfÞ in the range fA < f < fB :

Z fB Z 1 1 dðnÞ p

fA  f

dðnÞ dn  ln

þ dn  p 2 f1

fA n  f 0 nf  Z fB 1 3ð1  aÞ sin dðnÞ ¼ ln expð3sB Þ  dn : 3 pF 2 n f ð16Þ Further, equating (15) and (14) gives an integral equation for dðfÞ in the range 0 < f < 1: 252



Z fB dðnÞ p

fA  f

dðnÞ dn  ln

dn þ 2 f1

0 nf fA n  f  Z f 1 3 sin dðnÞ ¼ ln expð3sD Þ þ 2 dn : ð17Þ 3 pF 0 n

1  p

Z

1

In addition to these two questions, the total amount of oil, defined by the integral

Z

q ¼

yB

yA

Z



xðyÞdy

¼

fB

fA

expðsðnÞÞ sin dðnÞ

dn

xðnÞ pn ð18Þ

must be a given constant if no leaking takes place. The submergence depth of the barrier, l, is calculated using the second equation in (11):

Z

l ¼

1

fB

expðsðnÞÞ sin dðnÞ

dn : pn

ð19Þ

Equations (16)–(19) form a closed integral equation system. With a given set of parameters a, F, q, l, these equations are to be solved to yield the complete information about the flow, in particular, the values of sB , sD , fA , fB , dðfÞ values for fA < f < fB and 0 < f < 1 and the positions of the free surface and interface. This is what we called the direct approach, which is physically and intuitively more natural. But, it is more expensive computationally. Previously, Zhu and Strunin (2001) adopted an indirect approach in which a, F, sB , fB are taken as the input parameters and the volume of oil q and the submergence depth of the barrier l are a result of the calculation. Less computational effort associated with the indirect approach proposed by Zhu and Strunin (2001) is attributed to the fact that only two out of four integral equations need to be solved simultaneously first and also with given sB and fB , the limits in Eqs. (17) and (18) are fixed. Once the unknown function dðnÞ is found from the solution of Eqs. (17) and (18), the values of q and l can be easily obtained using (18) and (19) respectively. On the other hand, the indirect method has a fatal flaw; it is very difficult to get some physically useful information from the model. For example, one often wishes to examine how the increase of the volume of the trapped oil would change the interfacial position for a given Froude number and submergence depth of the barrier. One may also wish to see how the positions of free surface and interface change for a fixed volume of the trapped oil and submergence depth of the barrier. This motivated us to refine our numerical calculation, implementing the direct approach as we shall present in the next section. Spill Science & Technology Bulletin 7(5–6)

MODELLING THE CONFINEMENT OF OIL SPILL WITH FLOATING BOOMS

Numerical Implementation and Result Discussion On top of Eqs. (16)–(19) we have in our system of equations, there are other equations that arise from various geometrical and dynamical constrains. First of all, as shown by Zhu and Strunin (2001), the Bernoulli equation and the condition of motionlessness of the oil lead to an inequality of VA > VB . Hence, VA > 0 and point A can not be a stagnation point. Consequently, the interface re-attaches the barrier tangentially at point A: dA ¼ p=2:

ð20Þ

At the lowest point E of the barrier, we are looking for a solution in which the flow is supposed to separate smoothly, that is dE ¼ p=2. Far downstream, we confine our current problem to a downstream wavefree solution, which demands that dD ¼ 0:

ð21Þ

For simplicity, the value of a has been fixed at 0.85 throughout this paper. This is a typical value for the oil–water density ratio. The integral Eqs. (16) and (17) were discretised, which resulted in a system of nonlinear algebraic equations. We then called an IMSL subroutine NEQNF (a variation of NewtonÕs method of the nonlinear iteration) to solve this nonlinear algebraic system. In the indirect method, for every set of a, F, fB and sB , we solved the discretised system of equations and tabulated our results against various given fA and sD values. More specifically, instead of solving Eqs. (16), (17), (20) and (21) together, we adjusted fA and sD until Eqs. (20) and (21) were satisfied. As discussed in Zhu and Strunin (2001), there are two type of solutions: the stagnation point solution (with B as a stagnation point) and a nonstagnation point solution. Here, we should focus only on the latter, in which we set dB ¼ 0 since the velocity must be along the horizontal direction. Unlike the indirect approach, the convergence of the NEQNF in solving the nonlinear system of equations resulting from the direct approach turned out to be very sensitive to the change in magnitudes of fA and sD , especially fA . Consequently, sometimes the solution appeared to oscillate for a long time and no convergent solution could be found even though its indirect counterpart has proved to have a convergent solution. To cope with this problem, we designed an iteration loop in which solutions from our indirect approach is utilised to expedite the convergence. In this loop, we first solve Eqs. (16)–(20) with appropriate boundary conditions ðdA ¼ p=2 and dD ¼ 0Þ for a given set of a, F, q and l. The obtained values of sB , fB , Spill Science & Technology Bulletin 7(5–6)

dD and dA are then used as the input for the code we developed based on the indirect approach, which always exhibited very good convergence rate. As a result of calling this subroutine based on the indirect approach, we obtained a better approximation of fA and sD . Although the q and l are no longer the same as those specified, they are not too much different from those we input into the iteration loop in Step 1 either. Now using the newly obtained q and l in conjunction with the given a, F values, which were not altered in Step 2, we can proceed to the next round of iteration. Our numerical experiments showed that after about five to seven times of these two-step iterations, convergent solutions are obtained for a given set of a, F, q and l. Some typical results of the calculations from the direct method are presented in Figs. 3–5. In these plots we show how the flow profiles change against the parameters F, q, l. Fig. 3 shows the positions of the interface and free surface when l varies while q, F are kept constant. Clearly, with smaller submergence

Fig. 3 Computed free and interfacial positions with fixed F ¼ 2 and q ¼ 0:02. l ¼ 0:3(1), 0.4(2), 0.6(3).

Fig. 4 Computed free and interfacial positions with fixed F ¼ 4 and l ¼ 0:4. q ¼ 0:01(1), 0.06(2), 0.17(3). 253

S.-P. ZHU & D. STRUNIN

value of the terms having F as a coefficient turns out to be very small for all considered values of F. To certain extent, this smallness is a consequence of smallness of 1  a 0:15. Despite the value of a being a constant for the given fluids (water and oil) it is interesting to see what happens if a is not close to unity so that 1  a is not that small. For this purpose we conducted the calculations for a ¼ 0:2. Figure 5b shows that the faster the flow the shorter the depth of the contact point A; the block of oil is actually pushed up and backwards rather than down and forwards as one might have expected. This is right in line with the assumption we made; with the shear stress between the water and oil being ignored, a faster flow would exert larger pressure on the front of the interface to push it upwards. The correction of this perhaps replies on the next phase of modelling, which is to properly include the shear stress on the oil–water interface in an elegant way that modelling effort is not substantially increased. Another important observation we can make regarding the basic Eqs. (16) and (17) is that if F ! 1 the integrals containing F tend to 0, and the flow assumes some finite limiting shape.

Conclusions Fig. 5 Computed free and interfacial positions with fixed l ¼ 0:4 and q ¼ 0:02, F ¼ 1(1), 4(2): (a) a ¼ 0:85; (b) a ¼ 0:2.

depth, the detachment point B moves upstream and the reattachment point moves upwards. However, the interfacial position only seems to have changed slightly when l is changed from 0.3 to 0.6. The large difference of the downstream free surfaces is a direct result of different submergence depths of the barrier. The interfacial and free surface profiles for different values of q and fixed F, l are exhibited in Fig. 4. As q is increased the reattachment point A moves downwards as expected. As it approaches the vicinity of the tip of the barrier, E, convergent solutions do not exist any more. Such a situation infers that a critical value of q is reached beyond which the oil starts to leak under the barrier. Note that the points A and E never coincide exactly since otherwise the pressure at this point would be discontinuous (indeed, from the oil side the pressure would equal the weight of the oil column whereas from the free-surface side it would be zero). The results of altering parameter F while q, l are constants are shown in Fig. 5. One can see that the interface shape changes very slightly. Mathematically, this outcome is easy to understand after estimating the value of the terms in the integral equations. A typical 254

Some recent results of the model proposed by Zhu and Strunin (2001) are presented in this paper. The direct solution procedure enables us to quantitatively discuss the relationship among three important physical parameters, the Froude number, the submergence depth, and the amount of oil trapped by the barrier. Although all the results obtained so far appear to be in line with our assumptions, a further relaxation of some assumptions such as no shear stress on the water–oil interface may be necessary to improve the model. The inclusion of the shear stress on the interface would definitely make the model much more complicated than it is now, the research is being undertaken at the moment and the results will be reported in the near future. Acknowledgements—Authors gratefully acknowledge the financial support from the Australian Research Council under the grant A89700141 for this project.

References Cuesta, I., Grau, F.X., Giralt, F., 1990. Numerical simulation of oil spills in a generalized domain. Oil Chem. Pollut. 7, 143–159. Hoult, D.P., 1972. Oil spreading on the sea. Ann. Rev. Fluid Mech. 4, 341–368. Noye, J., Bills, P., Lewis, G., 1993. In: Stewart, D. et al. (Eds.), Prediction of Oil Slick Movement in Northern Spencer Gulf CTAC93, pp. 320–328. Proctor, R., Elliot, A.J., Flather, R.A., 1994. Forecast and hindcast simulations of the Braer oil spill. Marine Pollut. Bull. 28, 219– 229. Spill Science & Technology Bulletin 7(5–6)

MODELLING THE CONFINEMENT OF OIL SPILL WITH FLOATING BOOMS

Reed, M., Gundlach, E., 1989a. A coastal zone oil spill model: development and sensitivity studies. Spill Sci. Tech. Bull. 5, 411–449. Reed, M., Gundlach, E., 1989b. Hindcast of the Amoco Cadiz event with a coastal zone oil spill model. Spill Sci. Tech. Bull. 5, 451–476. Sebastiao, P., Soares, C.G., 1995. Modelling the fate of oil spills at sea. Spill Sci. Tech. Bull. 2, 121–131. Strunin, D.V., Zhu, S.P., 1998. An experimental study of the confinement of spilled oil by floating booms. In: Proc. of the 3rd Int. Conf. on Hydrodynamics, ICHD-98. UIAM Publishers, Seoul, pp. 839–844.

Spill Science & Technology Bulletin 7(5–6)

The Oil and Gas Journal, 1997. Cleanup continues on downgraded Japan Spill. 95 (28), 25. Times, 1989, 134 (9), 57. Turrell, W.R., 1994. Modelling the Braer oil spill––a retrospective view. Oil Chem. Pollut. 28, 211–218. Zhang, Y., Zhu, S.-P., 1996. Open channel flow past a bottom obstruction. J. Eng. Math. 30, 487–499. Zhu, S.P., Strunin, D.V., 2001. Modelling the confinement of spilled oil with floating booms. Appl. Math. Modell. 25, 713– 729.

255