安永寧，楊 鯤，王 瑩，李 晶

(交通運(yùn)輸部 天津水運(yùn)工程科學(xué)研究所，天津 300456)

20世紀(jì)90年代以來(lái)，國(guó)際上出現(xiàn)了很多的二維水力學(xué)模型。丹麥水力學(xué)研究所開(kāi)發(fā)的MIKE21二維數(shù)學(xué)模型，是其中應(yīng)用較為廣泛的商業(yè)模型之一。該模型可用于模擬河流、湖泊、河口、海灣、海岸及海洋的水流、波浪、泥沙及環(huán)境場(chǎng)，可為工程應(yīng)用、海岸及規(guī)劃提供完備、有效的設(shè)計(jì)條件和參數(shù)。該軟件的高級(jí)圖形用戶(hù)界面與高效計(jì)算引擎的結(jié)合，使其在世界范圍內(nèi)成為很多專(zhuān)業(yè)河口海岸工程技術(shù)人員不可缺少的工具，并曾在丹麥、埃及、澳洲、泰國(guó)等許多國(guó)家和地區(qū)得到成功應(yīng)用。如今該軟件已在國(guó)內(nèi)的一些大型工程中得到廣泛應(yīng)用。如：長(zhǎng)江口綜合治理工程、杭州灣數(shù)值模擬、南水北調(diào)工程、重慶市城市排污評(píng)價(jià)、太湖富營(yíng)養(yǎng)模型、香港新機(jī)場(chǎng)工程建設(shè)、臺(tái)灣桃園工業(yè)港興建工程等。

1 MIKE21模型簡(jiǎn)介

MIKE21是丹麥水力研究所(簡(jiǎn)稱(chēng)DHI)開(kāi)發(fā)的系列水動(dòng)力學(xué)軟件之一，屬于平面二維自由表面流模型。DHI以MIKE21作為研究手段，不斷在應(yīng)用中發(fā)展和改進(jìn)該軟件。20多a來(lái)，MIKE21在世界范圍內(nèi)大量工程應(yīng)用經(jīng)驗(yàn)的基礎(chǔ)上持續(xù)發(fā)展，其在平面二維數(shù)值模擬方面具有強(qiáng)大的功能。

1.1 MIKE21模型功能

1)用戶(hù)界面友好，屬于集成Windows圖形界面。

2)具有強(qiáng)大的前、后處理功能。在前處理方面，具有根據(jù)地形資料進(jìn)行網(wǎng)格劃分功能；在后處理方面，具有流場(chǎng)動(dòng)態(tài)演示及動(dòng)畫(huà)制作、計(jì)算斷面流量、實(shí)測(cè)與計(jì)算過(guò)程驗(yàn)證、不同方案比較等分析功能。

3)熱啟動(dòng)設(shè)置功能。用戶(hù)因各種原因暫時(shí)中斷了MIKE21模型，當(dāng)再次開(kāi)始計(jì)算時(shí)，只需將上次計(jì)算設(shè)置的熱啟動(dòng)文件調(diào)入便可繼續(xù)計(jì)算，極大地方便了計(jì)算時(shí)間。

4)干、濕節(jié)點(diǎn)和干、濕單元的設(shè)置功能，能較方便地進(jìn)行灘地水流的模擬。

5)強(qiáng)大的卡片設(shè)置功能，可以進(jìn)行多種控制性結(jié)構(gòu)的設(shè)置。

6)多種類(lèi)型水邊界條件的定義功能，如流量、水位或流速等。

7)廣泛應(yīng)用于二維水力學(xué)現(xiàn)象的研究功能。

1.2 MIKE21模型的2類(lèi)計(jì)算參數(shù)

1)數(shù)值參數(shù)：主要是方程組迭代求解時(shí)的有關(guān)參數(shù)。

2)物理參數(shù)：主要有床面阻力系數(shù)、動(dòng)邊界計(jì)算參數(shù)以及渦動(dòng)粘性系數(shù)等，取值或計(jì)算方法詳見(jiàn)DHI編撰的MIKE21用戶(hù)手冊(cè)(MIKE21 User Guide)。

2 MIKE21 HD計(jì)算原理

MIKE21軟件中的水動(dòng)力學(xué)模塊是其最核心的基礎(chǔ)模塊，可以模擬因各種作用力作用而產(chǎn)生的水位和水流變化以及模擬任何忽略分層的二維表面流[1-5]。

2.1 基本控制方程組

(1)

(2)

(3)

(4)

拖曳力系數(shù)可以根據(jù)Chezy系數(shù)C或Manning系數(shù)M計(jì)算

(5)

(6)

Manning系數(shù)可以根據(jù)底部糙率計(jì)算。

(7)

2.2 定解條件

1)邊界條件

(1)開(kāi)邊界

ζr=ζr(t)或ur=ur(t)或vr=vr(t)

(8)

式中，ζr或ur或vr為開(kāi)邊界γ上已知潮位、流速過(guò)程。

(2)閉邊界

u=0或v=0

(9)

2)初始條件

(10)

式中，ζ0為計(jì)算初始時(shí)刻潮位空間分布函數(shù)。

2.3 計(jì)算方法

模型求解采用非結(jié)構(gòu)網(wǎng)格中心網(wǎng)格有限法求解，其優(yōu)點(diǎn)為計(jì)算速度較快，非結(jié)構(gòu)網(wǎng)格可以擬合復(fù)雜地形。模型對(duì)計(jì)算區(qū)域內(nèi)灘地干濕過(guò)程，采用水位判別法處理，即當(dāng)某點(diǎn)水深小于一淺水深εdry時(shí)，令該處流速為零，灘地干出，當(dāng)該處水深大于εflood時(shí)，參與計(jì)算。

3 應(yīng)用實(shí)例

3.1 研究區(qū)概況

龍口灣位于渤海萊州灣的東北，是指屺姆島頭與界河口連線(xiàn)以東、呈對(duì)數(shù)螺線(xiàn)型半敞開(kāi)的海灣，為萊州灣的一個(gè)附屬海灣。灣內(nèi)水深不足10 m，灣外水深一般為10～20 m。研究區(qū)潮汐屬不規(guī)則半日潮，日不等現(xiàn)象明顯；潮流屬不規(guī)則半日潮流，漲潮方向由北向南，落潮時(shí)由南向北[6]，我們?cè)跀?shù)值模擬中用到的實(shí)測(cè)資料的平面位置分布見(jiàn)圖1。

圖1 研究區(qū)人工島群平面布置及調(diào)查站位分布圖Fig.1 Artificial islands cluster layout and investigation sites

3.2 模型網(wǎng)格設(shè)置

模擬大海域?yàn)檫|寧登沙河、山東雞鳴島兩點(diǎn)以及岸線(xiàn)圍成的北黃海及渤海海域。采用非結(jié)構(gòu)三角網(wǎng)格剖分計(jì)算域，三角網(wǎng)格能較好的擬合陸邊界，網(wǎng)格設(shè)計(jì)靈活且可隨意控制網(wǎng)格疏密，為了獲得更為清楚的潮流特征，對(duì)研究區(qū)海域網(wǎng)格設(shè)置進(jìn)行了加密處理。研究區(qū)海域網(wǎng)格劃分見(jiàn)圖2。

3.3 模型參數(shù)設(shè)置

3.3.1 水深和岸界

水深和岸界選取中國(guó)人民解放軍海軍航海保證部制作的1∶100萬(wàn)海圖(10011號(hào))，1∶15萬(wàn)(11370號(hào)、11570號(hào)、11710號(hào)、11770、1840號(hào)、1910號(hào))海圖及龍口港附近水深地形測(cè)量資料。

3.3.2 大海域模型水邊界輸入

1)開(kāi)邊界

引用遼寧登沙河和山東雞鳴島多年的潮位觀測(cè)資料，通過(guò)調(diào)和分析求得M2、S2、K1和O1四個(gè)主要分潮調(diào)和常數(shù)值輸入下式計(jì)算。

(11)

式中，σi是第i個(gè)分潮(這里共取4分潮：M2，S2，O1和K1)的角速度；Hi和Gi是調(diào)和常數(shù)，分別為分潮的振幅和遲角。

2)閉邊界

以大海域和工程周邊岸線(xiàn)作為閉邊界。

3.3.3 計(jì)算時(shí)間步長(zhǎng)和底床糙率

模型計(jì)算時(shí)間步長(zhǎng)根據(jù)CFL條件進(jìn)行動(dòng)態(tài)調(diào)整，確保模型計(jì)算穩(wěn)定進(jìn)行，最小時(shí)間步長(zhǎng)0.3 s。底床糙率通過(guò)Manning系數(shù)M進(jìn)行控制，M取0.3～0.45 m1/3/s。

3.3.4 水平渦動(dòng)粘滯系數(shù)

采用考慮亞尺度網(wǎng)格效應(yīng)的Smagorinsky[8]公式計(jì)算水平渦粘系數(shù)，表達(dá)式如下：

(12)

式中，cs=0.15，l為特征混合長(zhǎng)度。

(13)

3.4 潮位、潮流驗(yàn)證

將潮位、潮流實(shí)測(cè)數(shù)據(jù)采用傅氏分析方法[8]進(jìn)行調(diào)和分析后，作為數(shù)值模擬結(jié)果的驗(yàn)證依據(jù)。經(jīng)與渤海沿岸多個(gè)驗(yàn)潮點(diǎn)的潮位及龍口灣內(nèi)C1～C4測(cè)流站位的流速流向進(jìn)行驗(yàn)證后，結(jié)果表明，對(duì)應(yīng)觀測(cè)點(diǎn)上潮流模擬結(jié)果與實(shí)測(cè)潮流資料基本吻合，能夠較好地反映工程周邊海域潮流狀況。僅將龍口港的潮位驗(yàn)證曲線(xiàn)和C3和C4站位的流速流向驗(yàn)證曲線(xiàn)列于圖3～圖5。

圖3 潮位驗(yàn)證曲線(xiàn)(龍口港)Fig.3 Tide validation curve(Longkou Harbor)

圖4 流速流向驗(yàn)證曲線(xiàn)(C3站位)Fig.4 Current velocity and direction validation curve(C3 site)

圖5 流速流向驗(yàn)證曲線(xiàn)(C4站位)Fig.5 Current velocity and direction validation curve(C4 site)

3.5 潮流場(chǎng)數(shù)值模擬結(jié)果

3.5.1 萊州灣海域潮流場(chǎng)模擬結(jié)果分析

圖6是漲潮潮流場(chǎng)，計(jì)算域內(nèi)的潮流由萊州灣外向?yàn)硟?nèi)流動(dòng)，潮流流速一般為30～50 cm/s。北部海域潮流流向?yàn)槲髂舷?，至萊州淺灘處流速達(dá)到最大值，為96 cm/s。繞過(guò)萊州淺灘后，潮流轉(zhuǎn)向南流動(dòng)，進(jìn)入太平灣后，向偏東方向流動(dòng)。萊州淺灘以西至黃河口之間海域，潮流由東北往西南流，至灣頂附近逐漸轉(zhuǎn)為西向流，在黃河口南側(cè)近岸海域潮流由東南往西北流。

圖7為落潮潮流場(chǎng)，潮流分布情況與漲急時(shí)潮流場(chǎng)相似，只是流向與漲急時(shí)潮流場(chǎng)剛好相反。潮流整體由萊州灣內(nèi)向?yàn)惩饬鲃?dòng)，潮流流速一般為30～50 cm/s，在萊州淺灘處流速最大值可達(dá)82 cm/s。

圖6 萊州灣現(xiàn)狀漲潮時(shí)刻流場(chǎng)圖Fig.6 Current field in spring tide time in LaiZhou Bay

圖7 萊州灣現(xiàn)狀落潮時(shí)刻流場(chǎng)圖Fig.7 Current field in ebb tide time in LaiZhou Bay

3.5.2 人工島群海域潮流場(chǎng)模擬結(jié)果分析

圖8和圖9分別為人工島群建成后的漲潮和落潮預(yù)測(cè)潮流場(chǎng)。漲潮時(shí)模擬區(qū)內(nèi)潮流繞過(guò)屺姆島高角后，先向南流，后受到人工島群的阻擋，一部分以西南向繞過(guò)人工島群西部的防波堤后改向南流，并在防波堤南側(cè)的港池內(nèi)形成一個(gè)逆時(shí)針旋轉(zhuǎn)的小環(huán)流；一部分以東南向流入人工島群內(nèi)部，后沿人工水道流出。落潮時(shí)模擬區(qū)內(nèi)潮流一部分由南向北流，流經(jīng)人工島群西部的防波堤時(shí)，在防波堤?hào)|北側(cè)形成一個(gè)不太明顯的順時(shí)針旋轉(zhuǎn)的小環(huán)流；一部分沿人工島群內(nèi)部水道流出后，轉(zhuǎn)為西北向流向屺姆島。

圖8 人工島群建設(shè)后漲潮時(shí)刻流場(chǎng)圖Fig.8 Current field in spring tide time after construction of artificial islands

圖9 人工島群建設(shè)后落潮時(shí)刻流場(chǎng)圖Fig.9 Current field in ebb tide time after construction of artificial islands

4 結(jié) 語(yǔ)

根據(jù)MIKE21 HD模塊的計(jì)算原理，以萊州灣大范圍潮流場(chǎng)和龍口人工島群工程小范圍潮流場(chǎng)數(shù)值模擬的實(shí)例，展示了該模型的模擬效果。結(jié)果表明，該模型可較好地反應(yīng)實(shí)際情況，并進(jìn)一步驗(yàn)證了MIKE21 HD模塊模擬大范圍潮流場(chǎng)和小范圍局部潮流場(chǎng)的功能。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)學(xué)模型和數(shù)值試驗(yàn)在科學(xué)研究、工程設(shè)計(jì)中的地位和作用不斷提高。將水動(dòng)力數(shù)值模擬全過(guò)程設(shè)計(jì)成軟件包，配合良好的前、后處理工具，降低技術(shù)門(mén)檻，可使其通用性大大加強(qiáng)，并且通過(guò)圖形界面操作，可大大減輕工作量，減少不必要的重復(fù)工作。今后，水動(dòng)力數(shù)值模擬技術(shù)將朝著高效、高精度、可視化、軟件化等方向迅速發(fā)展。

參考文獻(xiàn)：

[1] 王領(lǐng)元.應(yīng)用MIKE 對(duì)河流一、二維的數(shù)值模擬[D].大連:大連理工大學(xué),2007.

[2] JTS/T 231-2-2010 海岸與河口潮流泥沙模擬技術(shù)規(guī)程[S].北京:人民交通出版社,2010.

[3] 謝亞力.錢(qián)塘江河口區(qū)治江圍涂對(duì)杭州灣風(fēng)暴潮影響研究[D].青島:中國(guó)海洋大學(xué),2006.

[4] 姜衛(wèi)星.黃浦江溢油事故的數(shù)值模擬研究[D].上海:同濟(jì)大學(xué),2007.

[5] 王萬(wàn)戰(zhàn),董利瑾.渤海流場(chǎng)基本特性的MIKE21 模擬研究[J].人民黃河,2007,29(10):32-33.

[6] 中國(guó)海灣志編纂委員會(huì).中國(guó)海灣志第三分冊(cè)(山東半島北部和東部海灣)[M].北京:海洋出版社,1991.

[7] SMAGORINSKY J. General Circulation Experiments with the Primitive Equations. I. The Basic Experiment[J]. Monthly Weather Review, 1963, 91(3): 99-164.

[8] 曲綿旭,王文海,豐鑒章,等.龍口灣自然環(huán)境[M].北京:海洋出版社,1995.