郝林南，楊小宸

（1.遼寧省防汛抗旱指揮部辦公室，遼寧沈陽110003；2.沈陽農(nóng)業(yè)大學(xué)，遼寧沈陽110044；3.遼寧省水利水電科學(xué)研究院，遼寧沈陽110003）

遼河口三維水動力數(shù)值模擬

郝林南1,2，楊小宸3

（1.遼寧省防汛抗旱指揮部辦公室，遼寧沈陽110003；2.沈陽農(nóng)業(yè)大學(xué)，遼寧沈陽110044；3.遼寧省水利水電科學(xué)研究院，遼寧沈陽110003）

文中基于有限體積海岸海洋模型FVCOM，建立了遼河口三維水動力數(shù)值模型，利用潮位、流速實測資料對模型進行驗證并分析遼河口三維水動力運動特性，可為進一步研究遼河口泥沙運動及鹽度、污染物擴散等提供基礎(chǔ)。

模型驗證；遼河口；三維水動力；數(shù)值模型

遼河口地區(qū)位于遼寧省的南部，是我國海岸線的最北端，面臨渤海；區(qū)域內(nèi)自東向西有大遼河口、遼河口（雙臺子河口）、大凌河口三個獨立入?？陂T。遼河流域?qū)贉貛Ъ撅L(fēng)氣候，年徑流量為89億m3；年降水量約為350～1 000 mm，并呈自東南向西北遞減的特點。流域夏季多暴雨，強度大、頻率高、集流快，常使水位陡漲猛落，造成下游地區(qū)洪澇。此外，遼河的含沙量較高，僅次于黃河、海河，為中國第三位，年輸沙量達2 098萬t。遼河口地區(qū)為平原型淤泥質(zhì)海岸，有著豐富的海陸相來沙，近百年來發(fā)展的結(jié)果已成為攔門沙發(fā)達的堆積型河口。不斷淤積的河口使中水河槽過水面積相對減少，使納潮量減少和泄洪能力降低。此外，遼河口地區(qū)還存在淡水資源嚴重不足，水環(huán)境污染嚴重。

遼河口近海水體中的泥沙運動及水環(huán)境污染等問題是水動力作用過程的直接表現(xiàn)，對遼河口水動力特點的合理描述是研究上述問題的基礎(chǔ)。因此，下文建立遼河口三維水運動數(shù)學(xué)模型，分析河口區(qū)域水流運動規(guī)律，可為整治工程提供指導(dǎo)依據(jù)，可緩解河口地區(qū)淤積情況，同時可為河口附近的港口選址、油氣開發(fā)平臺、航道疏通工程和海灘旅游等建設(shè)帶來的相關(guān)環(huán)境問題的分析提供技術(shù)支撐，從長遠來看，對于遼河口的資源開發(fā)利用、河口治理、環(huán)境變化和經(jīng)濟發(fā)展都具有重要的現(xiàn)實指導(dǎo)意義。

1 模型建立

基于有限體積海岸海洋模型FVCOM建立遼河口三維水動力數(shù)值模型。FVCOM模型是由美國馬薩諸塞大學(xué)海洋科技研究院與伍茲霍爾海洋研究所共同開發(fā)的基于無結(jié)構(gòu)網(wǎng)格、有限體積法、三維原始方程的海洋模型。模型包括動量方程、連續(xù)性方程、溫鹽輸運方程和狀態(tài)方程，通過二階紊流閉合模型封閉方程。水平方向使用非結(jié)構(gòu)化三角形網(wǎng)格，能更好地描述近岸地區(qū)的復(fù)雜岸線和建筑物，垂向采用坐標變換。數(shù)值方法采用有限體積法，能夠很好地將有限差分法離散結(jié)構(gòu)簡單、計算效率高的優(yōu)點和有限元方法處理復(fù)雜邊界的幾何靈活性結(jié)合起來，可以更好地保證具有復(fù)雜地形和岸線的近岸、河口區(qū)域的質(zhì)量、動量、鹽度以及熱量的守恒。

1.1 控制方程

FVCOM三維水動力模型基于自由表面的三維原始控制方程，由動量方程、連續(xù)性方程和狀態(tài)方程構(gòu)成，基于笛卡爾坐標系的三維原始方程可表示為：

其中：x，y和z為笛卡爾坐標系下的三維坐標，分別為東西方向、南北方向以及垂向的坐標，m；u，v和w分別為x，y和z軸的速度分量，m/s；g為重力加速度，m/s2；t為時間，s；ρ為密度，kg/m3；p為壓力，Pa；f為科氏力參數(shù)；Km為垂向渦粘系數(shù)，m2/s；Fu和Fv為水平動量擴散項，m/s；T為溫度，℃；S為鹽度，‰。

垂向采用坐標變換：

其中，ζ為自由表面，m；H為海床相對于基準面的距離，m。因此，D=H+ζ為總水深，σ坐標的變化范圍為[-1，0]。

1.2 數(shù)值離散

數(shù)值離散采用有限體積法，將計算區(qū)域劃分為若干互不重疊且相互鄰接的三角形網(wǎng)格單元。模型采用二維外模和三維內(nèi)模交替計算的方法，由外模態(tài)計算垂向平均流速和水位，由內(nèi)模態(tài)計算三維流速、泥沙運動和溫鹽輸運等物理量。二維外模將動量方程和連續(xù)性方程在控制體內(nèi)積分，通過改進的四階Runge-Kutta時間積分求解。三維內(nèi)模采用隱式與顯式結(jié)合的方法，對流項采用顯式求解，垂向擴散項采用隱式求解。具體的數(shù)值離散過程見文獻。

2 遼河口三維水動力數(shù)值模擬

2.1 模型設(shè)置

模型采用非結(jié)構(gòu)化的三角形網(wǎng)格，通過大、小模型嵌套計算，由大模型為小模型提供邊界潮位條件。大模型計算范圍為122.5°E以西海域，包括整個渤海及渤海海峽。大模型區(qū)域共包含12 856個節(jié)點，劃分為24 969個三角形單元，最小網(wǎng)格步長約為1 800 m。垂向平均分為5層，大模型開邊界水位由中國海域潮汐預(yù)報軟件Chinatide提供，外模計算時間步長為1.0 s。小模型計算范圍包括西至葫蘆島、東至鲅魚圈連線北側(cè)海域，同時小模型北端向遼河口上游延伸，其中雙臺子河上游延伸至盤錦市盤錦縣的盤山閘，大遼河延伸至鞍山海城市的三岔河，見圖1。小模型區(qū)域共包含18 118個節(jié)點，劃分為34 152個三角形單元，在河口及河道內(nèi)加密，最小網(wǎng)格步長約為30 m。垂向平均分為10層，河口地形取自1∶500 00海圖，河道地形資料采用2011年大斷面資料。小模型開邊界潮位由大模型計算結(jié)果提供，河流流量邊界根據(jù)水文站實測資料確定，雙臺子河近似采用六間房測站數(shù)據(jù)，大遼河近似采用邢家窩棚和唐馬寨測站之和，大凌河根據(jù)凌海測站資料顯示流量相對較小，暫不考慮。模型外模時間步長取為0.1 s，內(nèi)模時間步長為1.0 s。

圖1 小模型計算范圍及地形圖

2.2 模型驗證

收集到2013年9月20日10:00至 2013年9月21日16:00大潮期間在計算海域的水文全潮觀測結(jié)果，用于模型驗證。測站位置見圖1，其中包括潮位測站T1，T2及分層流速、流向測站C1～C4。計算時間選為2013年9月18日10:00至2013年9月21日16：00，給水動力計算充分的穩(wěn)定時間。通過大、小模型嵌套計算，由大模型為小模型提供開邊界潮位條件，河流流量根據(jù)水文站數(shù)據(jù)確定為雙臺子河191 m3/s，大遼河111 m3/s。計算結(jié)果與實測數(shù)據(jù)吻合較好，能夠較好地反映遼河口海域潮流變化情況。

計算結(jié)果表明：落潮流速略大于漲潮流速，表層和底層流向基本一致，表層流速大于底層流速。圖2給出了漲急時刻和落急時刻計算區(qū)域表層的流場圖，由圖2中可以看出，遼河口海域潮流基本呈往復(fù)流運動。漲潮期間，外海潮流自西南向東北方向近岸傳播，蓋州灘東側(cè)深槽（即測站附近）為主流區(qū)，流速大、流勢強，水流主要通過這股漲潮主流分別流向雙臺子河口和大遼河口。在雙臺子河口外約12 km處，河流流量和外海潮流形成頂托，并向西北延伸至近岸區(qū)域，該現(xiàn)象可能會導(dǎo)致泥沙在此處滯留并沉積。落潮期間，流向與河流流向一致，雙臺子河口水流在蓋州灘處分為兩股流向外海，蓋州灘東側(cè)深槽依然是落潮流的主流區(qū)。

3 結(jié) 語

基于FVCOM模型建立的遼河口三維水動力數(shù)值模型，采用非結(jié)構(gòu)化網(wǎng)格建立大、小模型并模擬了大潮期間遼河口海域的水動力運動，研究結(jié)論主要有：

1）通過與實測資料對比，模型計算結(jié)果與實測數(shù)據(jù)吻合較好，此文所建模型能夠合理描述遼河口三維水動力特性，可為進一步研究遼河口泥沙運動及鹽度、污染物擴散等提供基礎(chǔ)。

2）遼河口海域潮流運動呈往復(fù)流，蓋州灘東側(cè)深槽為漲落潮的主流區(qū)。落潮流速略大于漲潮流速，表層流速大于底層流速。漲潮期間，雙臺子河口外潮水與河流的頂托作用明顯，可能對泥沙產(chǎn)生滯留沉積作用。

［1］Chen C，Liu H，Beardsley R C.An unstructured,finitevolume,threedimensional，primitive equation ocean model：application to coastal ocean and estuaries［J］.Journal of Atmospheric and Oceanic Technology,2003，20（1）：159-186.

TV131

A

1002－0624（2016）12－0044－02