胡艷海 周林飛

摘要：石佛寺水庫是遼河干流上唯一一座控制性水利工程，將其作為研究對象，基于水動力軟件MIKE 21，建立了石佛寺水庫二維非穩(wěn)態(tài)水動力模型。利用質(zhì)量平衡進(jìn)行模型可靠性分析，結(jié)果表明，不合理洪水過程線是有效的，即模型是可靠的;利用2018年水庫實測流量資料進(jìn)行模型主要參數(shù)的率定，采用偏差統(tǒng)計法進(jìn)行誤差分析，結(jié)果表明，絕對值偏差為72 m3/s，相對誤差為7.5%，模型計算誤差在允許范圍內(nèi)，滿足精度要求;實現(xiàn)了水庫流場變化規(guī)律的模擬計算。利用所建立的模型，進(jìn)行石佛寺水庫水動力情況的模擬與分析，進(jìn)一步驗證了模型的可靠性，可為水庫的日常管理、生態(tài)治理、水質(zhì)改善等提供技術(shù)支持。該模型的建立為水庫水質(zhì)模型提供了運行載體，是水質(zhì)參數(shù)運移、消散、衰減過程研究的基礎(chǔ)。

關(guān)鍵詞：石佛寺水庫;水動力模型;數(shù)值模擬;參數(shù);率定

中圖分類號： S271文獻(xiàn)標(biāo)志碼： A文章編號：1002-1302（2020）13-0265-06

收稿日期：2019-08-13

基金項目：遼寧省水利科學(xué)技術(shù)項目（編號：20170147）。

作者簡介：胡艷海（1993—），男，遼寧凌源人，碩士研究生，主要從事水文學(xué)及水資源研究。E-mail：2207567047@qq.com。

通信作者：周林飛，博士，副教授，主要從事水環(huán)境研究。E-mail：zlf924@163.com。目前，在水庫研究中水動力學(xué)過程發(fā)揮著重要的作用，對水體與氣體、水體本身、水體與沉積物之間的能量轉(zhuǎn)化、物質(zhì)傳輸以及能量交換等一系列過程產(chǎn)生一定程度的影響，是水庫變化規(guī)律研究的前提[1]。從20世紀(jì)80年代開始，平面二維水動力模型逐漸變得成熟，研究者不斷對數(shù)值模擬的方法進(jìn)行改善與革新。在大型數(shù)字計算機(jī)以及地理信息系統(tǒng)技術(shù)的不斷快速發(fā)展過程中，平面二維水動力模型（MIKE 21 flow model），以完全圣維南方程為基本前提來求解兩維水流的水力學(xué)方法，可提供有效的參數(shù)和完整的設(shè)計條件。該模型利用有限體積算法靈活完成復(fù)雜地形的處理工作，在短時間內(nèi)能夠反映水庫的水動力流場的變化運移情況。本研究以石佛寺水庫為研究背景，建立其二維非穩(wěn)態(tài)水動力模型，并經(jīng)過模型參數(shù)的率定和水動力模型可靠性的驗證，最終得到合理可行的水庫水動力模型，為進(jìn)一步和水質(zhì)模型的建立和耦合提供前提條件。

遼河河源位于河北省平泉縣七老圖山脈的光頭山，該河流河源段由東遼河與西遼河2段組成，其干流主要流經(jīng)河北省、內(nèi)蒙古自治區(qū)、吉林省、遼寧省4個省份，流域全長和流域面積分別為1 345 km和21.9萬km2。石佛寺水庫根據(jù)地形分類，屬于河道形平原水庫，是遼河干流上唯一一座同時兼顧防洪、生態(tài)和供水三大功能的大型控制性水利工程，其中防洪是水庫的最主要功能，主要體現(xiàn)在對下游人口、基礎(chǔ)設(shè)施和耕地的保護(hù)，水庫行政區(qū)處于遼寧省沈陽市新城子區(qū)境內(nèi)，與沈陽市區(qū)相距47 km。石佛寺水庫主要水利樞紐有泄洪閘、主壩和副壩。樞紐工程等級為Ⅱ等工程，永久性建筑物等級為2級，共布置泄洪閘16孔，閘孔總寬度為248.5 m。水庫開工于2003年5月，完工于2005年10月，并于2009年開始進(jìn)行生態(tài)蓄水。

1MIKE 21水動力計算原理

MIKE 21 flow model（FM，平面二維非穩(wěn)態(tài)數(shù)學(xué)模型）憑借著其自身強(qiáng)大的卡片設(shè)置功能，干、濕節(jié)點和干、濕單元的設(shè)置功能，以及非常優(yōu)秀的前、后處理功能，被廣泛地應(yīng)用于二維水動力學(xué)現(xiàn)象的研究中。該模型常被推薦為水庫、湖泊、河流、河口和海岸水流的二維仿真模擬工具。

1.1計算原理

MIKE 21 flow model（FM）是基于三向不可壓縮和雷諾（Reynolds）值分布的納維耶斯托克斯（Navier-Stokes）方程，同時服從于靜水壓力的假定和布辛涅斯克（Boussinesq）假設(shè)。由于各水力要素沿垂向變化速率遠(yuǎn)不及水平方向的變化顯著，因此，考慮將沿垂向積分的二維非恒定流淺水方程作為模型的控制方程。連續(xù)方程見式（1）;x、y方向的動量方程如式（2）和式（3）所示[2]。

1.2計算方法

二維水動力數(shù)值模擬的計算方法有以下3種形式：有限差分法、有限單元法和有限體積法。基于各計算方法有著較為明顯的優(yōu)缺點，用戶全面綜合的考慮計算速度、要求的計算精度以及相應(yīng)的研究對象等多方面因素進(jìn)行計算方法的選擇和應(yīng)用，本研究針對計算區(qū)域的空間離散選擇有限體積法（finite volume method），將該連續(xù)統(tǒng)一體細(xì)分為互相不重疊的單元，單元可以是任意形狀的多邊形，本研究只考慮非結(jié)構(gòu)化的三角形和四邊形混合形式[3]。

2石佛寺水庫水動力模型的建立

2.1數(shù)據(jù)收集與前處理

2.1.1數(shù)據(jù)收集模型研究區(qū)域起始于遼河干流石佛寺水庫入口斷面，終止于石佛寺水庫出口斷面，收集2個斷面2017年7月1日至2017年8月31日的水位和流量數(shù)據(jù)，以及2018年8月石佛寺水庫入口、出口流量數(shù)據(jù)。采用石佛寺水庫電子版地形（比例尺為1 ∶5 000）作為網(wǎng)格建立的地形資料。

2.1.2數(shù)據(jù)前處理數(shù)據(jù)的提?。簽榱司_對石佛寺水庫庫區(qū)地表形態(tài)的描述和模型所需網(wǎng)格的生成，使用南方cass軟件聯(lián)合Auto CAD軟件，對石佛寺水庫地形圖內(nèi)部高程點進(jìn)行提取，共提取8 621個高程點信息，并利用GPS獲取庫區(qū)水面高程點，隨后利用可折疊式水深測桿在獲得高程點位置測量水深值，計算得到該點庫底高程，完成一個采樣點高程驗證，庫區(qū)共布置采樣點237個;利用Auto CAD對水庫地形圖進(jìn)行邊界高程點提取，并采用同樣的方法使用GPS和測深桿對提取邊界高程數(shù)據(jù)進(jìn)行采樣驗證，邊界共布置采樣點225個，完成建模高程數(shù)據(jù)的提取工作。

數(shù)據(jù)格式轉(zhuǎn)化：將所提取的高程點，采用Excel和Notepad轉(zhuǎn)換成擴(kuò)展名為*.xyz的格式文件，以適合于MIKE軟件模型，供網(wǎng)格生成使用。所使用的高程點投影坐標(biāo)為Beijing1954坐標(biāo)，并以此坐標(biāo)作為工作空間坐標(biāo)和導(dǎo)入庫區(qū)高程點數(shù)據(jù)的坐標(biāo)。

2.2計算網(wǎng)格制作

研究區(qū)計算網(wǎng)格質(zhì)量的好壞很大程度上決定了模型計算是否穩(wěn)定。采用Beijing_1954_GK_Zone_21N投影坐標(biāo)，利用MIKE ZERO網(wǎng)格生成器完成網(wǎng)格制作，計算網(wǎng)格共2 832個，并采用網(wǎng)格編輯器對網(wǎng)格進(jìn)行優(yōu)化處理。工作區(qū)選擇完成后，導(dǎo)入已經(jīng)制作好的*.xyz邊界文件，進(jìn)行邊界命名同時對邊界以200 m距離進(jìn)行均勻化處理，以最大三角形面積、最小三角形角度、最大節(jié)點個數(shù)作為網(wǎng)格生成的控制條件，其值分別為20 000 m2、30°、150 000 個。利用Smoth mesh功能對生成的網(wǎng)格進(jìn)行了10 000次光滑處理，并使用網(wǎng)格分析模塊分別檢查最小計算時間、最小網(wǎng)格角度和最小三角形邊長，即嚴(yán)重畸形的網(wǎng)格，并利用網(wǎng)格編輯菜單反復(fù)修改使地形文件（Bathemetry）可靠，進(jìn)而保證模型的穩(wěn)定性和準(zhǔn)確度[4]。地形文件是模型完成計算所需要的區(qū)域，該文件的后綴名為*.mesh，本研究將其命名為Beijing54SFS.mesh。

研究區(qū)計算網(wǎng)格見圖1，庫區(qū)二維地形見圖2。

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

MIKE 21 FM二維非穩(wěn)態(tài)水動力模型參數(shù)分為計算參數(shù)和物理參數(shù)2類。計算參數(shù)包括計算步長（time step interval）、計算步數(shù)（number of step）、克朗數(shù)（Courant-Friedrich Levy number）和計算精度（solution technique）等影響方程組迭代求解的有關(guān)參數(shù);物理參數(shù)包括干濕水深（flood and dry）、渦黏系數(shù)（eddy viscosity）、底摩擦力（bed resistance）、水邊界條件（boundary conditions）等。

2.3.1計算參數(shù)考慮到數(shù)據(jù)的完整性以及可用性，選擇石佛寺水庫入口及出口斷面2017年主汛期（2017-07-01—2017-08-31）水位、流量資料，歷時62 d。因此，模擬時間為2017年7月1日至2017年8月31日，模擬時段為62 d，計算步長為 10 800 s，計算步數(shù)為488步。

克朗數(shù)（Courant-Friedrich Levy number，簡稱CFL）是指一種能夠保證模型穩(wěn)定計算的關(guān)鍵因素，CFL值正比于水動力模塊計算所需交換數(shù)據(jù)時間、反比于三角形網(wǎng)格的面積[5]。通過對網(wǎng)格反復(fù)修改并試驗，來保證模型的計算穩(wěn)定性，同時為確保模型的精度和準(zhǔn)確性，采用高階算法，最終確定水動力模塊計算所需交換數(shù)據(jù)時間為0.1～30 s，CFL值為0.85。

2.3.2物理參數(shù)初始條件（initial condition）：將石佛寺水庫正常運行情況下的水位作為模型運算的初始水位，即46.2 m;根據(jù)水動力模型所具有的能夠在很短的時間內(nèi)就能夠使計算達(dá)到穩(wěn)定而形成真實流速的特點，將x、y方向的初始流速條件定義為0 m/s。

邊界條件（boundary condition）：水動力模型通常情況下邊界條件分為陸地邊界和開邊界。對于開邊界又被細(xì)分為上邊界（進(jìn)口）和下邊界（出口）。在MIKE 21中一般情況下將進(jìn)口斷面的流量作為上邊界，將出口斷面的水位過程作為下邊界。因此本研究將石佛寺水庫入口斷面2017年7月1日至2017年8月31日實測流量數(shù)據(jù)作為進(jìn)口邊界，將石佛寺水庫出口斷面2017年7月1日至2017年8月31日實測水位數(shù)據(jù)作為出口邊界，并相應(yīng)地將流量和水位數(shù)據(jù)制成符合MIKE軟件建模要求的時間序列文件（即*.dfs0），以此作為模型運算的邊界條件。

渦黏系數(shù)（eddy viscosity）：描述空間上和時間上的不確定性物理過程，在控制方程中即為相應(yīng)的附加應(yīng)力項。選擇Smagorinsky公式的推薦值0.28作為最終的取值。

底摩擦力（bed resistance）：底摩擦力用于防止模型的計算失穩(wěn)，作為模型的率定參數(shù)之一，通常情況下將以曼寧系數(shù)作為相應(yīng)的輸入項，且通常選取模型的缺省值作為最終的取值，即32 m1/3/s。

干濕水深（flood and dry）：干、濕動邊界主要是為了避免模型計算中出現(xiàn)失穩(wěn)問題，即當(dāng)模型在水深小于干水深時對應(yīng)的網(wǎng)格不參與計算。需要設(shè)定3個參數(shù)，干水深（drying depth）、淹沒水深（flooding water depth）、濕水深（wetting depth）。本研究3個值均取模型的默認(rèn)值，分別為hdry=0.005 m、hflood=0.05 m、hwet=0.1 m。

3模型驗證與分析

3.1模型的水量平衡分析

在二維水動力模型結(jié)果輸出時，將文件形式設(shè)置為水量平衡。其中水流可作為該結(jié)果文件的輸出項目，并且在水流項目組件下又可以包括誤差（error）等項目，誤差能夠用于模型的可靠性分析。在進(jìn)行模型可靠性分析時，通過不合理洪量過程線來詮釋，即當(dāng)不合理洪量過程線很低時說明模型輸出結(jié)果的總水量變化符合能量守恒，也就是模型可靠性良好。所建立的研究區(qū)二維非穩(wěn)態(tài)水動力模型的不合理洪量過程線見圖3，僅有7月1日至7月4日誤差線的值為非負(fù)，其他值均小于0，其中最小值為-3.0×10-7 m3/s，據(jù)此數(shù)據(jù)充分說明誤差線是有效的，即模型是可靠的。

3.2模型的率定與誤差分析

利用MIKE 21 FM水動力學(xué)模型對石佛寺水庫庫區(qū)進(jìn)行模擬，在確定模型參數(shù)的基礎(chǔ)上，利用2018年8月石佛寺水庫入口斷面的流量數(shù)據(jù)進(jìn)行

計算，將出口斷面模擬的流量結(jié)果輸出，并使用二次摩擦定律公式τbρ0=cf·ub·|ub|，在經(jīng)驗值范圍內(nèi)調(diào)節(jié)研究區(qū)水域拖拽系數(shù)（cf）與底床流速（ub），ρ0為研究區(qū)水體密度，最終糙率值（τb）調(diào)整為0.022 2時較為合適。在MIKE 21軟件中所使用的曼寧系數(shù)為糙率值τb的倒數(shù)，因此最終曼寧系數(shù)被確定為45 m1/3/s時率定效果最佳。

率定結(jié)果如圖4所示，圖中虛線為2018年8月出口斷面流量實測值，實線為模擬值，從圖中可以看出二者擬合程度良好。采用偏差統(tǒng)計法進(jìn)行誤差計算，結(jié)果表明最大絕對值偏差為72 m3/s，相對誤差為7.5%，計算結(jié)果在誤差允許范圍內(nèi)，說明模擬結(jié)果比較合理，滿足精度要求，表明所構(gòu)建的模型有效，可以為構(gòu)建水質(zhì)模型提供可靠的前提條件，進(jìn)一步驗證了模型的可靠性。

3.3水動力情況模擬與分析

3.3.1水深分析利用建立的石佛寺水庫水動力模型輸出水庫水深分布，來直觀表達(dá)庫區(qū)水面高程分布情況。圖5-a表示水庫正經(jīng)歷重現(xiàn)期為300年一遇的洪水時的水深分布，圖5-b表示水庫恢復(fù)正常運行狀態(tài)時的水深分布。無論是在洪水時段還是正常運行時段庫區(qū)水深分布情況大致趨勢呈現(xiàn)出相似的規(guī)律，在泄洪閘附近水深出現(xiàn)最大值、庫中主槽區(qū)域水深值次之、水深在遠(yuǎn)離主槽區(qū)出現(xiàn)最小值，在洪水時段整個庫區(qū)水深值絕大部分高于水庫正常運行狀態(tài)時的水深值。據(jù)此，表明模型的運行結(jié)果與實際情況一致，說明模型模擬庫區(qū)水深分布趨勢是正確的。在研究水體水質(zhì)改善問題時，水深分布提供的庫區(qū)不同分區(qū)的水深值，可為后續(xù)水質(zhì)模型研究水質(zhì)監(jiān)測點的布置和水質(zhì)參數(shù)擴(kuò)散起始地點的選擇提供直觀的參考依據(jù)。

3.3.2流速與流態(tài)分析根據(jù)水動力模型輸出石佛寺水庫正常運行情況下的流速、流態(tài)分布，來表達(dá)庫區(qū)流速值分布及流態(tài)變化情況，如圖6、圖7所示。根據(jù)流速、流態(tài)分布圖可知庫區(qū)主槽區(qū)域流速最顯著，其平均值大致為0.60 m/s;泄洪閘區(qū)域流速值較水閘區(qū)域小些，平均流速約為 0.36 m/s;由于水生植物的局部阻礙作用，水生植物區(qū)域的水流速度最小，其平均流速約為 0.12 m/s。水動力模型提供的庫區(qū)流場變化情況為污染物質(zhì)的運移消散路徑的初步確定提供了依據(jù)。因此，可根據(jù)石佛寺水庫的流速流態(tài)，合理布置水生植物，以增加局部糙率來改變污染物質(zhì)運移路徑，及利用水生植物凈化水質(zhì)。

3.3.3瞬時流速分析輸出水庫正常運行情況主槽區(qū)、相對密閉區(qū)域模擬時段內(nèi)的瞬時流速，如圖8所示，可以看出主槽區(qū)流速明顯大于相對密閉區(qū)。

瞬時流速作為影響污染物運移、擴(kuò)散的主要影響因素，其大小對污染物濃度變化的影響也相當(dāng)顯著，具體表現(xiàn)為流速大時水質(zhì)濃度相對較小，反之較大，因此，要注意相對密閉區(qū)的水質(zhì)監(jiān)測。輸出水庫正常運行情況下水閘位置模擬時段內(nèi)的瞬時流速，將水閘位置的流速情況與實際情況相比，如圖9所示，可以看出輸出結(jié)果與實測結(jié)果基本吻合，進(jìn)一步驗證了模型的可靠性。

4結(jié)論

本研究采用MIKE 21 flow model（FM）建立石佛寺水庫二維非穩(wěn)態(tài)水動力模型，通過模型穩(wěn)定分析及率定，模型真實可信。利用水動力模型輸出的庫區(qū)模擬時段內(nèi)任意時刻的各水力要素（如水深、流速、流態(tài)等）的流場分布情況，以反應(yīng)各要素在水庫的不同區(qū)域和不同時間段的變化關(guān)系，模型模擬結(jié)果與實測結(jié)果對比后的吻合程度表明，模擬結(jié)果可信，趨勢正確，能夠用于水庫水情變化情況的預(yù)測，可為水庫日常管理提供依據(jù)。

模型輸出的庫區(qū)不同分區(qū)的水深值，可為水質(zhì)模型水質(zhì)監(jiān)測點的布置和水質(zhì)參數(shù)擴(kuò)散起始地點的選擇提供依據(jù);運用水庫水動力模型計算獲得的水流流速可為水質(zhì)濃度提供水動力條件。這為進(jìn)一步加載水質(zhì)模型進(jìn)行水庫水質(zhì)研究提供了基礎(chǔ)保障。

參考文獻(xiàn)：

[1]郭雪蕊. 基于水動力學(xué)的武漢市東湖水質(zhì)模擬[D]. 西安：西安理工大學(xué)，2018.

[2]周紅玉，劉操，吳曉輝. 大型跨流域調(diào)水工程對水源地水動力的影響研究[J]. 中國農(nóng)村水利水電，2017（12）：104-108.

[3]李前君，楊理棟，東愛明. 防城港市漁港工程潮流數(shù)值模擬研究[J]. 中國水運，2016（11）：105-110.

[4]鄭婷婷，徐明德，景勝元，等. 汾河水庫水動力及水質(zhì)數(shù)值模擬[J]. 水利水運工程學(xué)報，2016（3）：105-113.

[5]姜灝. 遼寧省石佛寺水庫2013年汛情分析[J]. 中國防汛抗旱，2014，24（2）：72-73.彭峰生. 核桃仁脫皮機(jī)設(shè)計[J]. 江蘇農(nóng)業(yè)科學(xué)，2020，48（13）：271-275.doi：10.15889/j.issn.1002-1302.2020.13.054