王璟，婁安剛，曹振東

（中國海洋大學(xué)環(huán)境工程與科學(xué)學(xué)院，山東青島266100）

崖城油氣田附近海域溢油漂移軌跡數(shù)值模擬預(yù)測

王璟，婁安剛，曹振東

（中國海洋大學(xué)環(huán)境工程與科學(xué)學(xué)院，山東青島266100）

通過采用不規(guī)則的三角網(wǎng)格和有限體積法的FVCOM模式，建立三維潮流數(shù)值模型。利用大海域計(jì)算得到的調(diào)和常數(shù)值作為開邊界的輸入值，模擬出崖城附近海域的潮流和潮位變化情況。在潮流、潮位驗(yàn)證正確的前提下，利用歐拉—拉格朗日追蹤方法，建立了溢油軌跡預(yù)測模型，進(jìn)行崖城油氣田附近海域溢油中心軌跡的預(yù)測，同時(shí)預(yù)測了溢油漂移的平均速率和油膜抵達(dá)敏感區(qū)的時(shí)間，為油氣田實(shí)施應(yīng)急措施提供技術(shù)支持。

FVCOM模型；崖城；溢油；數(shù)值模擬

1 引言

崖城位于海南島最南端，距離三亞市西40 km處，其附近海域擁有大量的油氣田，在油氣田開發(fā)和運(yùn)營的過程中可能會對周圍的海域和一些敏感資源造成一定的污染[1]，因此崖城附近海域水動力場的改變直接影響到這些工程造成的污染物輸送。本文以崖城13-4平臺附近發(fā)生突發(fā)性溢油為例，通過建立三維潮流模型模擬出崖城附近海域的潮流潮位情況，并以此為基礎(chǔ)分別在平均風(fēng)和最強(qiáng)風(fēng)下預(yù)測出溢油的漂移軌跡、平均速率、抵岸時(shí)間等。

在之前對南海的研究中，多以北部灣[2-3]為主要研究對象，而對海南島南部海域的研究較少。本文重點(diǎn)針對海南島南部海域，對其潮流潮位隨時(shí)間變化情況加以研究。

2 潮流數(shù)值模型

2.1 流體動力學(xué)基本方程組

運(yùn)用FVCOM三維模型進(jìn)行數(shù)值模擬，采用sigma坐標(biāo)，三維非線性潮波運(yùn)動方程為：

式中：g為重力加速度；f=2Ωsinφ為柯氏參數(shù)；Ω為地轉(zhuǎn)角速度；φ為地理緯度；ζ為自靜止水面算起的水位高度；h為自靜止水面算起的水深；H=h+ζ；u、v、w為流速，p為壓強(qiáng)，ρ為海水密度，τ為切應(yīng)力。

2.2 邊界條件

海面邊界條件：

固體側(cè)邊界條件：

在開邊界上水位邊界條件：

這里H和g分別是調(diào)和常數(shù)的振幅和遲角，下標(biāo)C代表著M2、S2、K1、O1四個(gè)分潮中的一個(gè)。ω代表著分潮頻率，f為交點(diǎn)因子，u為交點(diǎn)訂正角，V0是天文潮的初位相。

2.3 計(jì)算海域和網(wǎng)格設(shè)置

由于計(jì)算域開邊界上調(diào)和常數(shù)實(shí)際中難以得到，本文首先建立了大海域海流模型（見圖1），由大海域計(jì)算得出的調(diào)和常數(shù)值[4]，作為計(jì)算海域開邊界的輸入值。

圖1 大海域網(wǎng)格及地形圖

本文所建立的海域數(shù)學(xué)計(jì)算域范圍及網(wǎng)格見圖2，從圖中可以看出，三角形網(wǎng)格較好的反映出計(jì)算域內(nèi)的海岸線和地形特征。其中坐標(biāo)范圍為16.5°—19°N，108°—110.5°E，水深及岸線根據(jù)海圖確定。崖城13-4氣田周圍海域網(wǎng)格加密，最小網(wǎng)格邊長800 m。垂向δ分20層。

圖2 計(jì)算海域網(wǎng)格及地形圖

水界強(qiáng)迫水位：利用大海域的計(jì)算結(jié)果作為工程海域水界的輸入資料。水界可提供4個(gè)分潮(M2、S2、K1、O1)的調(diào)和常數(shù)。由于是不正規(guī)全日潮海域[3]，水界可輸入下列調(diào)和方程：

式中，右方第1項(xiàng)表示平均日分潮，第2項(xiàng)表示平均半日分潮。

3 模擬結(jié)果驗(yàn)證分析

3.1 潮流與潮位的驗(yàn)證

在計(jì)算海域中，選取A、B、C三點(diǎn)（坐標(biāo)位置見表1）分別作為潮流驗(yàn)證點(diǎn)和潮位驗(yàn)證點(diǎn)。在這些點(diǎn)分別將數(shù)值計(jì)算的結(jié)果與實(shí)測的資料進(jìn)行了驗(yàn)證，其驗(yàn)證結(jié)果見圖3。

表1 驗(yàn)證點(diǎn)坐標(biāo)位置

圖3 流向驗(yàn)證曲線

從圖3c—3h可以看出，流速的大小以及方向，轉(zhuǎn)流發(fā)生時(shí)刻的計(jì)算值與實(shí)測值基本一致；圖3a和3b顯示A、B兩站位的潮位振幅和位相計(jì)算值亦與實(shí)測值基本一致。

潮位和潮流的驗(yàn)證結(jié)果表明建立的潮流模型是可行的，適合本海區(qū)。

3.2 評價(jià)海域潮流計(jì)算結(jié)果

圖4、5給出漲、落潮中間時(shí)的表層計(jì)算潮流場的潮流場圖。由圖中可以看出崖城13-4氣田附近潮流分布表現(xiàn)為往復(fù)流：漲潮為西北向流，落潮為東南向流，落潮潮流流速大于漲潮潮流流速。通過計(jì)算得知，崖城13-4氣田附近區(qū)域漲潮最大流速為48 cm/s，平均流速為31cm/s，落潮最大流速為60 cm/s，平均流速為38 cm/s。海南島東南部海域流速遠(yuǎn)低于其他海域流速，這是由于東南部有一由100 m變深至1000 m的海底陡坡，與文獻(xiàn)[5]的結(jié)論一致。

圖4 漲潮中間時(shí)計(jì)算潮流場

圖5 落潮中間時(shí)計(jì)算潮流場

圖6和圖7為崖城13-4氣田附近區(qū)域漲潮時(shí)和落潮時(shí)垂向最大流速分布圖，由圖可以看出，最大流速隨著深度的增加而減小，與文獻(xiàn)[6]的結(jié)論一致。

圖6 漲潮中間時(shí)最大流速垂向分布

圖7 落潮中間時(shí)最大流速垂向分布

4 溢油預(yù)測模型

4.1 油膜漂移軌跡

在環(huán)境動力模型提供的環(huán)境動力參數(shù)的基礎(chǔ)上，采用歐拉-拉格朗日追蹤方法，進(jìn)行油膜中心軌跡的預(yù)測[7]。油膜中心漂移速度，取決于海面風(fēng)速與表層流，是空間和時(shí)間的函數(shù)，其值用油膜中心點(diǎn)所在網(wǎng)格點(diǎn)上的速度內(nèi)插而得?？臻g每個(gè)網(wǎng)格節(jié)點(diǎn)上的x、y方向上的速度在某時(shí)刻為：

式中，Vrx、Vry為網(wǎng)格點(diǎn)上表層潮流流速加上風(fēng)海流的x、y方向分量，表層潮流流速由環(huán)境動力學(xué)模型求出，風(fēng)海流的模擬根據(jù)埃克曼公式計(jì)算：

式中，w為風(fēng)速，Ф為緯度。Vwind為網(wǎng)格點(diǎn)上的風(fēng)速,α為風(fēng)因子，θ0為風(fēng)向，θ為油粒子受風(fēng)影響的漂移偏角。θ的取值與風(fēng)速的大小有關(guān)，公式為：

油粒子漂移軌跡計(jì)算公式為：

式中：S0為初始時(shí)刻，S為油膜中心點(diǎn)所在位置，Vl(x(t),y(t),t)為拉格朗日追蹤速度，

由于空間和時(shí)間不同，流況不同，有時(shí)風(fēng)速、風(fēng)向也不同，所以在不同地點(diǎn)、不同時(shí)刻發(fā)生溢油后所追蹤到的油膜中心運(yùn)移軌跡就不同。

4.2 油膜擴(kuò)展輸移預(yù)測

剪流和湍流引起的擴(kuò)散過程屬于隨機(jī)運(yùn)動，可用隨機(jī)走動法實(shí)現(xiàn)模擬。由于每個(gè)粒子的隨機(jī)運(yùn)動而導(dǎo)致整個(gè)粒子云團(tuán)在水體中的擴(kuò)散過程。對于水體表面隨機(jī)擴(kuò)散過程可用下式描述[8]：

式中：ra為α=(x,y)方向上的湍動擴(kuò)散距離；R為[-1，1]間均勻分布隨機(jī)數(shù)。kα為α方向上的湍流擴(kuò)散系數(shù)，Δt為時(shí)間步長。

第i個(gè)粒子在Δt時(shí)段內(nèi)的位移可表示為：

式中：rx，ry為在x，y方向上的隨機(jī)移動距離；ui，vi為拉格朗日速度在x，y方向上的分量。

4.3 油的化學(xué)過程

溢油在輸移和擴(kuò)展的過程中，也同時(shí)發(fā)生各種生物和化學(xué)變化過程，使得海上溢油的量不斷減少。

4.3.1 溢油的揮發(fā)

溢油揮發(fā)過程受油性質(zhì)、油厚度、風(fēng)及油組分控制。采用Stiver和Mackay提出的一個(gè)暴露模式來計(jì)算油的揮發(fā)：

式中，B為系數(shù)，常取10.3；TG為揮發(fā)曲線梯度；T為油的表面溫度，通常與大氣溫度相近；T0為初始時(shí)油揮發(fā)溫度；Pa為大氣壓；V為油分子體積；R為大氣常數(shù)；θ為揮發(fā)系數(shù)，常取2.5×10-3，Uw為風(fēng)速；T0、TG的數(shù)值常參考如下常數(shù)：

由YC13-1氣田混合凝析油的試驗(yàn)得知API度(15.56℃)為43.70。

4.3.2 溢油的乳化

乳化過程受風(fēng)速、波浪、油的厚度、環(huán)境溫度、油風(fēng)化程度等因素的影響，一般用含水率來表示乳化程度。

式中：

YWi為第i個(gè)油粒子含水率，Uw為風(fēng)速，WAx為油的含蠟量，Asph為油的瀝青質(zhì)量含量，η0為油的無水動力粘性系數(shù)，YWsat為穩(wěn)定含水量，K1、K2為常數(shù)，分別為5.0×10-7和1.2×10-5。ηi為乳化后油的運(yùn)動粘性系數(shù)。

式中，ηoil為乳化前油的運(yùn)動粘性系數(shù)。

4.4 溢油量及溢出方式

4.4.1 溢油點(diǎn)

崖城13-4氣田主體區(qū)，坐標(biāo)17°37'35.697"N，109°09'44.123"E。

4.4.2 溢油方式及溢油量

點(diǎn)源一次性連續(xù)排放2 h，溢油量為50 m3，平均溢油速率為25 m3/h。

4.5 風(fēng)場

根據(jù)提供的風(fēng)速風(fēng)向分布圖以及風(fēng)向玫瑰圖，本文選取表2中所示年平均最大風(fēng)頻風(fēng)場（風(fēng)速：m/s），同時(shí)選取表3中所示5個(gè)不利極值風(fēng)，考慮溢油對近海海域可能造成的最大影響。

4.6 預(yù)測結(jié)果

（1）油膜漂移軌跡

圖8、9為崖城13-4氣田主體區(qū)發(fā)生溢油后，分別在平均風(fēng)和極值風(fēng)（見表2、3所示風(fēng)場）下所模擬出油膜的漂移軌跡示意圖。

圖8 平均風(fēng)下油膜漂移軌跡示意圖

圖9 極值風(fēng)下油膜漂移軌跡示意圖

表2 年平均風(fēng)場

表3 年極值風(fēng)場

（2）油膜漂移的平均速率、抵岸時(shí)間及抵達(dá)功能區(qū)的時(shí)間

通過計(jì)算可知，無論在什么時(shí)候發(fā)生溢油，油膜經(jīng)過一定時(shí)間會溢出計(jì)算域或者抵達(dá)岸邊。如果溢油量偏少并且油膜的漂移時(shí)間較長，油膜會在溢出計(jì)算域前或抵岸前消失。表4、5給出平均風(fēng)和極值風(fēng)況與流場耦合的條件下，在崖城13-4氣田主體區(qū)發(fā)生溢油后，油膜漂移的距離，平均速率、抵岸時(shí)間和抵達(dá)重點(diǎn)海域功能區(qū)的時(shí)間等。由表5可以看出，發(fā)生溢油之后，在出現(xiàn)極值不利風(fēng)向下，油膜會在較短時(shí)間內(nèi)抵達(dá)需要重點(diǎn)保護(hù)的海南島南部重點(diǎn)海域功能區(qū)（最短為偏西南風(fēng)向下，21.4 h即可到達(dá)敏感海域，對其造成影響）。因此需要各方在溢油事故發(fā)生后及時(shí)采取有效措施對溢油進(jìn)行攔截。

表4 平均風(fēng)下溢油不同條件下油膜漂移時(shí)間等

表5 極值風(fēng)下溢油不同條件下油膜漂移時(shí)間等

4 結(jié)語

本文采用FVCOM三維潮流模式模擬海南島南部海域的潮流場，看出崖城油氣田附近海域的潮流表現(xiàn)為往復(fù)流，漲潮為西北向流，落潮為東南向流，并且落潮潮流大于漲潮潮流。在潮流模擬正確的情況下進(jìn)行了崖城附近海域溢油軌跡的模擬，并預(yù)測了油膜抵達(dá)油氣田附近重點(diǎn)海域功能區(qū)的時(shí)間。在出現(xiàn)極值不利風(fēng)向下，油膜會在較短時(shí)間內(nèi)抵達(dá)需要重點(diǎn)保護(hù)的海南島南部重點(diǎn)海域功能區(qū)（最短為偏西南風(fēng)向下，21.4 h即可到達(dá)敏感海域，對其造成影響），油氣田附近一旦發(fā)生溢油，應(yīng)立即采取行之有效的應(yīng)急措施，保護(hù)溢油可能到達(dá)的重點(diǎn)保護(hù)海域，防止其遭受污染。

[1]萬修全,鮑獻(xiàn)文,吳德興,等.膠州灣及其鄰近海域潮流和污染物擴(kuò)散的數(shù)值模擬[J].海洋科學(xué),2003,27(5)：31-36.

[2]徐振華,雷方輝,婁安剛,等.北部灣潮汐潮流的數(shù)值模擬[J].海洋科學(xué),2010,34(2)：10-14.

[3]趙昌,呂新剛,喬方利,等.北部灣潮波數(shù)值研究[J].海洋學(xué)報(bào), 2010,32(4)：1-11.

[4]曹德明,方國洪.北部灣潮汐和潮流的數(shù)值模擬[J].海洋與湖沼, 1990,21(2)：105-113.

[5]孫洪亮,黃衛(wèi)民.北部灣潮汐潮流的三維數(shù)值模擬[J].海洋學(xué)報(bào), 2001,23(2)：1-8.

[6]葉安樂.分潮流最大流速發(fā)生時(shí)刻隨深度變化的特征[J].海洋通報(bào),1984,3(2)：1-8.

[7]婁安剛,王學(xué)昌,孫長青,等.膠州灣海面溢油軌跡的數(shù)值模擬[J].黃渤海海洋,2001,19(1)：1-8.

[8]婁安剛,吳德星,王學(xué)昌,等.三維海洋溢油預(yù)測模型的建立[J].青島海洋大學(xué)學(xué)報(bào)：自然科學(xué)版,2001,31(4)：473-479.

The numerical simulation and forecast of spilled oil trajectory in Yacheng oil and gas field sea area

WANG Jing,LOUAn-gang,CAO Zhen-dong
(College of Environmental Science&Engineering,Ocean University of China,Qingdao,266100 China)

Based on the unstructured grid finite-volume coastal ocean model(FVCOM),the complete three-dimensional fine-solution equation group of shallow-sea tide is used to compute the tidal current and sea-level variations in Yacheng sea area.The open boundary conditions are provided by the open oceanic model.The simulation of the tidal current and sea-level are verified with the observation data.By Euler-Lagarangian tracing method,a trajectory forecasting model for spilled oil is established to predict the trajectory of the oil-film centre,the averaging moving rate and the reaching time of the oil to the sensitive area in Yacheng oil and gas field sea area. It offers some effective technologies for the emergency response to the oil and gas field.

FVCOM model;Yacheng;oil spill;numerical simulation

book=19,ebook=19

X55

：A

：1003-0239（2012）03-0073-08

2011-05-09

王璟（1987-），男，碩士研究生，主要從事海洋環(huán)境動力和數(shù)值模擬研究。E-mail：qiaokeli16520245@163.com