大凌河口水動(dòng)力數(shù)值模擬分析

郭常來(lái)，崔 健，孫秀波，馬玉祥，艾叢芳

1.中國(guó)地質(zhì)調(diào)查局 沈陽(yáng)地質(zhì)調(diào)查中心，遼寧 沈陽(yáng) 110034，2.大連理工大學(xué)/海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024

0 引言

大凌河發(fā)源于遼寧省凌源市，于錦州凌海市注入渤海遼東灣.大凌河全長(zhǎng)447 km，流域面積達(dá)2.33×104 km2，是遼西地區(qū)最大的河流［1］，其下游及河口三角洲地區(qū)的大凌河口濕地，是我國(guó)重要的蘆葦濱海濕地，具有不可替代的社會(huì)經(jīng)濟(jì)價(jià)值和自然生態(tài)價(jià)值［2］.近年來(lái)，由于大凌河流域內(nèi)人口增多和工農(nóng)業(yè)飛速發(fā)展，以及氣候變化等原因，大凌河徑流有所減少［3-4］.而在河口地區(qū)，由于河流斷面擴(kuò)大，水流速度驟減，因徑流減小而產(chǎn)生的問(wèn)題更為明顯［5］.目前，對(duì)于河道徑流和潮汐的數(shù)學(xué)模型主要是基于一維、二維水動(dòng)力數(shù)學(xué)模型；周躍華等［6］基于GIS與一維水動(dòng)力模型對(duì)陶樂(lè)防洪保護(hù)區(qū)漫溢洪水進(jìn)行了風(fēng)險(xiǎn)分析；孫玲玲［7］通過(guò)建立二維非恒定流數(shù)學(xué)模型模擬了黃壁莊水庫(kù)洪水演進(jìn)過(guò)程；鎖曉南等［8］基于二維水動(dòng)力數(shù)學(xué)模型對(duì)黃河四排口河段進(jìn)行洪水演進(jìn)模擬分析水動(dòng)力場(chǎng)變化狀況.基于二維水動(dòng)力數(shù)學(xué)模型目前廣泛應(yīng)用于河床變形與泥沙運(yùn)移中，對(duì)復(fù)雜地形的計(jì)算精度高、可靠性強(qiáng).本研究基于MIKE 21中的水動(dòng)力模塊，考慮大凌河的河道徑流和潮汐的影響，對(duì)遼東灣北部海域潮汐進(jìn)行模擬，重點(diǎn)分析大凌河口附近海域的水動(dòng)力特征，為進(jìn)一步研究大凌河口海域的污染物輸運(yùn)及泥沙運(yùn)動(dòng)提供準(zhǔn)確的水動(dòng)力基礎(chǔ)條件（圖1）.

圖1 大模型計(jì)算域地形分布、測(cè)站位置及網(wǎng)格劃分Fig.1 Terrain distribution，station location and grid division in large-sized model calculation domain

1 研究方法

MIKE 21是一款用于模擬河流、湖泊、河口、海灣、海岸及海洋的水流、波浪、泥沙及環(huán)境的工程軟件包.本研究采用MIKE 21的水動(dòng)力模塊建立水動(dòng)力數(shù)值模型.該模塊基于Boussinesq［9］假定的二維不可壓縮的Navier-Stokes方程［10］，其控制方程組如下：

1）連續(xù)性方程

2）運(yùn)動(dòng)方程

式中：t表示時(shí)間；η表示自由水面水位（m）；h表示總水深（m）；ρ表示水體密度；ρ0表示水體相對(duì)密度；g為重力加速度；S表示源匯項(xiàng)的流量值；pa表示大氣壓力；us、vs表示源匯項(xiàng)的流速值；f＝2ωsinφ，其中ω是地球自轉(zhuǎn)角速度，φ是地理緯度；（τsx，τbx）、（τsy，τby）是x、y方向表面風(fēng)和海底剪切應(yīng)力的分量；Txx、Txy、Tyx、Tyy是橫向應(yīng)力，包括黏性應(yīng)力、湍流摩擦、平流摩擦；分別表示x、y方向垂線平均速度，采用（4）式計(jì)算.

在該模塊中采用基于單元中心的有限體積法進(jìn)行離散［11］，空間上采用不可重疊的非結(jié)構(gòu)網(wǎng)格劃分，能夠較好地?cái)M合河口地區(qū)復(fù)雜的岸線邊界，滿足對(duì)復(fù)雜地形模擬計(jì)算的需要［12］.

2 研究區(qū)域

2.1 模型建立及設(shè)置

為了更準(zhǔn)確地得到遼東灣北部海域潮流開(kāi)邊界條件，本研究采用大、小嵌套模擬的方法，建立水動(dòng)力數(shù)值模型［13-14］.其中大模型計(jì)算域?yàn)椴澈：Ｓ蚣安糠贮S海海域，范圍大致為117.3—122.4° E、36.5—41.2° N；小模型計(jì)算域?yàn)檫|東灣北部部分海域，范圍大致為121.1—122.3° E、40.2—41° N.采用非結(jié)構(gòu)網(wǎng)格對(duì)計(jì)算域劃分，對(duì)近岸及地形復(fù)雜處進(jìn)行了網(wǎng)格加密，以更好地?cái)M合復(fù)雜的地形條件.在本文中，水深和岸線根據(jù)海圖確定，同時(shí)參考研究地區(qū)最新岸線情況進(jìn)行修正，大、小模型計(jì)算域地形分布及網(wǎng)格劃分見(jiàn)圖1、圖2.

圖2 小模型計(jì)算域地形分布及網(wǎng)格劃分Fig.2 Terrain distribution and grid division in small-sized model calculation domain

本模型中設(shè)置閉邊界處法向速度為零，大模型開(kāi)邊界采用煙臺(tái)至大連之間的實(shí)測(cè)連續(xù)潮位數(shù)據(jù).在大尺度水動(dòng)力模型驗(yàn)證良好的基礎(chǔ)上，導(dǎo)出小模型所需要的潮位邊界資料.徑流邊界主要考慮遼東灣北部沿岸的大凌河、小凌河、遼河和大遼河4條河流的淡水輸入，其中小凌河、遼河和大遼河徑流采用月平均數(shù)據(jù)，大凌河徑流采用實(shí)測(cè)日均徑流數(shù)據(jù)［15］.

2.2 模型驗(yàn)證

利用模型對(duì)研究海域進(jìn)行潮位的數(shù)值模擬，選取2005年10月大潮期間的模擬值與環(huán)渤海8個(gè)代表性潮位測(cè)站的潮汐表值進(jìn)行對(duì)比，以驗(yàn)證模型的有效性，潮位站點(diǎn)見(jiàn)表1，對(duì)比結(jié)果如圖3所示.可以看出，本模型潮位模擬結(jié)果與實(shí)測(cè)潮汐資料吻合較好，在一個(gè)周期內(nèi)潮位變化基本一致.總體來(lái)看，模型能準(zhǔn)確地反映渤海海域的水動(dòng)力變化特性［16-18］.

圖3 環(huán)渤海代表測(cè)點(diǎn)潮位對(duì)比圖Fig.3 Comparison of tide levels at representative stations around Bohai Sea

表1 環(huán)渤海代表潮位站一覽表Table 1 Representative tide stations around Bohai Sea

3 討論與分析

3.1 整體流場(chǎng)特征分析

模型由“冷態(tài)”啟動(dòng)，即在數(shù)值模型進(jìn)行模擬計(jì)算時(shí)，初始時(shí)刻計(jì)算域范圍內(nèi)水位均設(shè)置為0［19-20］，如圖4所示.經(jīng)模擬發(fā)現(xiàn)，t＝2 h時(shí)水位降低的落潮過(guò)程已傳播至渤海范圍，t＝6 h時(shí)落潮過(guò)程傳播至大凌河河口地區(qū)，t＝12 h時(shí)落潮后的漲潮過(guò)程完全傳播至大凌河河口地區(qū).潮位過(guò)程每隔6 h漲潮與落潮交替變化，約12 h為一個(gè)完整潮位變化周期，每天兩次潮起潮落，與實(shí)際潮位變化規(guī)律一致.

圖4 不同時(shí)刻流場(chǎng)狀態(tài)Fig.4 Flow field states at different times

圖5給出了模擬期間遼東灣北部海域最低潮位時(shí)刻和最高潮位時(shí)刻的流場(chǎng)分布結(jié)果.可知，遼東灣北部海區(qū)最低潮位為－0.82 m，最高潮位為2.13 m.經(jīng)模擬計(jì)算，遼東灣平均潮差為2.4 m.該海區(qū)漲潮歷時(shí)要略短于落潮歷時(shí)，并且由于河口河道變窄，由外海傳入的潮波由于摩擦等作用潮能發(fā)生衰減，潮差沿徑流向上游逐漸減小.

圖5 遼東灣北部海域流場(chǎng)分布Fig.5 Distribution of flow fields in northern Liaodong Bay

3.2 大凌河口潮流分布

通過(guò)校驗(yàn)較好的水動(dòng)力數(shù)值模型，模擬了2018年10月遼東灣北部海域的潮流場(chǎng)的分布.模擬結(jié)果顯示，小潮期海域潮流場(chǎng)分布與大潮期基本一致，小潮期流速較大潮期流速稍小.圖6給出了大潮期模擬海域的低潮時(shí)、漲急時(shí)、高潮時(shí)、落急時(shí)潮流場(chǎng)分布.

圖6 大凌河口海域不同時(shí)刻流場(chǎng)分布狀態(tài)Fig.6 Flow field distribution at different times in the sea area of Dalinghe River estuary

大潮期間低潮時(shí)潮流場(chǎng)數(shù)值模擬結(jié)果表明，大凌河河口附近海域潮流整體流向?yàn)槟衔飨?，流速普遍低?.3 m/s，由外海至岸邊潮流流速逐漸減小.近岸處表現(xiàn)為沿岸流，流速普遍低于0.1 m/s.此時(shí)，在潮流的作用下，大凌河河道徑流朝外海加速涌出.

大潮期間漲急時(shí)潮流場(chǎng)數(shù)值模擬結(jié)果表明，大凌河河口附近海域潮流整體由南西向北東流動(dòng)，速度普遍介于0.3～0.5 m/s之間；小凌河河口至大凌河河口近岸潮流表現(xiàn)為由西至東向的沿岸流，速度普遍低于0.2 m/s.在漲潮的作用下，潮流涌向大凌河河道，大凌河口口門處流速較小.

大潮期間高潮時(shí)潮流場(chǎng)數(shù)值模擬結(jié)果表明，大凌河河口附近海域潮流整體發(fā)生倒轉(zhuǎn)，潮流由東向西流動(dòng)，流速約0.3 m/s.近岸處表現(xiàn)為沿岸流，流速普遍低于0.1 m/s.

大潮期間落急時(shí)潮流場(chǎng)數(shù)值模擬結(jié)果表明，大凌河河口附近海域潮流整體發(fā)生轉(zhuǎn)流，整體潮流流向?yàn)槟衔飨?，流速普遍介?.1～0.4 m/s之間；在落潮的作用下，大凌河徑流加速涌出河道，大凌河口口門處流速約為0.25 m/s.

3.3 潮流動(dòng)力特性分析

（1）漲潮動(dòng)力

遼東灣北部海域在漲潮期間潮流流向基本與遼東灣東、西岸線平行，呈現(xiàn)為北東向.由數(shù)值模擬的結(jié)果可知，在漲潮初期，由于受到大凌河口淺灘地形的影響，漲潮流主要在大凌河口淺灘外側(cè)及西側(cè)近海海域活動(dòng).在漲急時(shí)刻及漲潮后期階段，隨著漲潮動(dòng)力的增強(qiáng)，潮流流向主要是北東和北北東向，受河口淺灘地形的影響，潮流流速較弱，形成大范圍的緩流區(qū).

（2）落潮動(dòng)力

在落潮階段，河口水動(dòng)力過(guò)程主要來(lái)自潮水的進(jìn)退.在落潮初期，河口內(nèi)淺灘被水淹沒(méi)，但受地形影響，水深較淺，流速較小，在河口附近形成了大范圍的緩流區(qū).到了落潮后期，落潮動(dòng)力增強(qiáng)，潮間淺灘露出，水流主要經(jīng)潮溝流出，從而加大了潮溝的落潮水流動(dòng)力.大凌河的徑流量較小，徑流流出口門后，河道斷面突然展寬，水動(dòng)力快速消減，難以形成明顯的出海主流.

4 結(jié)論

基于大凌河口水動(dòng)力數(shù)值模型模擬結(jié)果，分析了大凌河口附近海域的水動(dòng)力變化特征.由于大凌河口緊鄰遼河口、大遼河口，該海域?qū)儆谌牒：恿鳑_淡水的交匯混合區(qū)，潮流同沖淡水水流相互作用，加之受潮溝密布和潮灘復(fù)雜地形影響，形成了該地區(qū)獨(dú)特的流場(chǎng)結(jié)構(gòu).

大凌河河口附近海域潮流屬于遼東灣海流系統(tǒng)的一部分，其主流方向呈北東－南西向，與遼東灣軸向基本一致.漲潮流（北東向）流速略大于落潮流（南西向）流速，最大漲潮流約為0.52 m/s，最大落潮流速約為0.4 m/s，潮流漲落平均潮流強(qiáng)度的分布大體與該海域等深線相適應(yīng).