二维海面溢油预报模型的研究及应用

第11卷第7期中国水运V ol .11

N o.72011年7月Chi na W at er Trans port J ul y 2011

收稿日期：作者简介：臧士文，大连理工大学海岸及近海工程国家重点实验室。

二维海面溢油预报模型的研究及应用

臧士文，梁书秀，孙昭晨

（大连理工大学海岸及近海工程国家重点实验室，辽宁大连116023）

摘

要：文中通过对海面溢油行为与归宿的研究，结合水动力模型，建立了二维海面溢油预报模型。基于“油粒子”

思想，将油膜平均分成具有一定体积的小油斑，每个小油斑在对流、湍流作用下扩展、漂移、扩散。根据Lagran ge 追踪法及粒子云团随机走动原理，考虑风的影响，模拟了溢油在海面漂移、扩散过程。将溢油模型应用到大连新港输油管线爆炸的实际溢油事故中，并将模拟结果与实测资料进行了对比。关键词：溢油；La gra n ge 追踪；模拟中图分类号：X 55文献标识码：A 文章编号：1006-7973（2011）07-0081-04一、前言

随着石油工业和海上运输业的快速发展，海洋石油开采及水上运输量逐年增加，溢油事故频繁发生，大量溢油泄入海洋，海洋生态环境遭到严重破坏。溢油污染严重危害人们

的健康，制约沿海经济的发展，严重的溢油事故会长时间影

响海洋生态环境，并会导致一些海洋生物灭绝。为了应对突发性溢油事故，各国政府部门相应建立了海洋溢油灾害应急响应系统，而这一系统的一个主要部分即是应用溢油数学模型预测海上溢油的漂移、扩散、蒸发乳化等行为归宿。

海上溢油行为是一个复杂的物理化学过程，不仅受风、海流、等水动力因素影响，而且和溢油本身的化学性质有关。目前的研究将溢油行为分为动力学过程和非动力学过程，前者主要包括溢油的漂移、扩展、等物理过程，后者则包括溢油的蒸发、乳化、等生化反应。由于受风应力及波浪力的作用，油膜在扩展过程中逐渐变薄并破碎成油滴。对于溢油轨迹的运动过程，目前大多数学者将“油粒子”思想运用到溢油模型中。“油粒子”思想最早由J oh an s en 和Elliot 等人提出，与传统的用对流扩散方法模拟溢油过程不同，“油粒子”思想将溢油分成大量的小油滴，基于lagran ge 追踪法，模拟每个小油滴在海水中的漂移轨迹，直观的反应了油膜在风、波浪等作用下破碎、分离的实际情况。

本文建立了二维海上溢油预报归宿模型，模拟了海上溢油的漂移、扩散过程。并将模型应用到大连新港输油管线爆炸溢油事故中，与实际观测资料进行了对比。

二、海上溢油过程1．溢油扩展过程

溢油扩展过程受重力、惯性力、粘性力和表面张力影响。Fay [1]（1971）把溢油扩展过程分为三个阶段，并提出了计算各个阶段油膜直径的公式。由于假定油膜为圆形而且没有考虑风及湍流的影响，用Fay 方法计算油膜的扩展范围要小于实际观测的油膜范围。Leh r et al.[2]（1984）考虑了风的影响对Fay 方程做了修改：

2/32/31/21/31/34/3

2270()40()wind o o

A V t V U t

ρρρρ=+（1）

式中，A 为油膜扩展面积（2m ），w o ρρρ=，w

ρ为水的密度，o ρ为油的密度；V 为溢油总体积；wind

U 为风速（m /m in ）；t 为时间（min ）。

根据Leh r 公式，每个时间步的溢油扩展面积便可得到。2．溢油漂移、扩散过程

油粒子的定义是将油滴看做很小的圆球，直径为1~100mm 。而在这个变化范围油膜将被分成数量巨大的小油滴，在计算中将耗费大量时间，并且计算机的容量有限，所以实际模拟溢油过程无法考虑的如此精细。本文将初始油膜平均分为一定数量的小油斑，对每个小油斑规定随时间而变化的相对参考坐标系。用潮流场的流速及风生流速计算油斑的水平漂移，根据粒子云团的随机走动原理[3]计算油斑的扩散过程。这样，便可计算出每个油斑在各个时间步长的坐标，从而可以模拟油膜的轨迹和形状。

（1）水平流速

每个油斑的水平速度可由如下公式计算：d t W

W U U K U =+JG J G J G （2）

式中，t U 为深度平均二维潮流流速，可由二维水动力模

型计算；W U 为海面以上10m 风速；W K 为风漂流因子，反应风应力对水流大小的影响，取值0.03~0.035。

（2）湍流扩散过程

通常可由粒子云团的随机走动原理来计算紊动扩散。由湍流速度引起的粒子位移可由随机走动距离来代替。Al-Rab eh et al.[4]（1989）研究指出，每个油斑的随机走动距离为：

[]1

12h S R D t

=（3）

式中，[]1

0R 为在0~1之间变化的随机数；h D 为水平扩散系数，本文取7m 2/s ；

3．油斑的位移

每个油斑的位移可由下式计算：()cos ()sin x dx y dy L t U t S L t U t S θθ

=+=+（4）

式中，()x L t ，()y L t 分别是油斑,x y 方向的水平位移；

2011-04-20

82中国水运第11卷

dx U ，dy U 分别是,x y 方向的水平流速；

[]1

2R θπ=（5）

三、二维水动力模型的建立

本文采用FVCOM 模式[5]计算二维潮流场。1．控制方程

海岸附近的大尺度潮流运动，水平尺度通常远大于垂向尺度，当忽略垂直方向的变化时，将σ坐标下的三维N-S 方程沿垂向积分，得到深度平均的潮流基本方程组：

()()0uD vD t x

y

ζ

++

=（6）200

'0111

[()]

sx bx x x

uD u D uvD

fvD

t x y

gD D gD D d d d x x x DF G

ζρσσσρσρττρ++=++

++∫∫∫（7）

20

'

0111

[()]

sy by

y y

vD uvD v D fuD

t x y

gD D gD D d d d y y y DF G ζρσσσρσρττρ++=++

++∫∫∫（8）

其中

2

2

(

)

x x

x u D uv D u D uv D G DF DF x

y

x

y

=

++

（9）2

2

(

)

y x

y uvD v D

uvD v D

G DF DF x y

x y

=

++（10）

水平扩散项可近似的给定为：

(2)[(

)]

x m m u u v DF A H

A H x x

y y

x

≈

+

+

（11）[()](2)

y m m u v v

DF A H A H x y x y y

≈++（12）2()x m m u u v

DF A H

A H x x y y x

≈

++（13）(

)2y m m u v v DF A H A H x

y x y y

≈

++（14）

这里的符号“

”表示由底至海表的垂直积分。例如，

对一个给定的变量ψ，则有

1

d ψψσ

=∫2．数值求解

FVCOM 在网格设计上采用非结构化三角形网格，数值求解上采用有限体积法。在给定的控制体内，对连续性方程（12）求面积积分，得到

'

'()()[]n s uD vD dxdy dxdy v Dds t x y

ξ=+=∫∫∫∫∫（15）其中，n v 是通过节点相连的相邻三角形中心点所围成的切面的法向水平速度；'

s 为绕此封闭切面的轨迹。方程（15）通过修正过的四届Ru n ge-Ku t ta 时间积分法进行求解运算。这个被修正后的多步时间积分方法具有二阶计算精度，具体的积分过程由下列各式给出：

0()0[(

)()]

n

j j NT j n n n n

n

n n

R R x y D x y D ζζζζ====

+∑（16）1014,2

k k

k

n j j

j j

j

tR ζ

ζ

ζζαζζ+==（17）

其中，1,2,3,4k =和1234

(,,,)(1/4,1/3,1/2,1)αααα=。上

标n 代表第n 个时间步长。

j ζ

是由通过与j ζ所在节点相

连的三角形中心点和相邻三角形每条边的中心点所围成的闭合线的面积。n

m u 和n

m v 分别定义为

(()),(())n n n

n

m m u u NT m v v NT m ==（18）t 为外膜时间步长。

21221221

2;m m

m m m m

x x x x x x +==（19）

212212212;m m m m m m

y y y y y y +==（20）

相似的，每一个三角形面积，对动量方程组（7）和（8）求积分，具体求解方法同上，在此不再赘述，具体计算方法可参见文献[5]。

四、模型在大连湾海域的应用

1．水动力模型的验证

2010年7月16日18时50分左右，大连新港输油管线爆炸发生溢油事故，持续溢油到17号10时左右。故根据实际情况，将溢油类型定为静止点源连续溢油，应用本模型对此次溢油事故进行了为期9天的数值模拟。将模拟潮位与

大连湾内A 、B 、C 、D 4点实测潮位进行了对比。图1为模拟海域，图2为4个验潮站站位。

大连新港

旅顺军港

大小长山岛

大连湾

图1计算海域示意图

D

B

C

A

图2大连湾潮位监测站图

2121212221

m m m m m m m

m m m m v u v u