石佛寺水庫二維水動力數值模擬

胡艷海 周林飛

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

關鍵詞：石佛寺水庫;水動力模型;數值模擬;參數;率定

中圖分類號： S271文獻標志碼： A文章編號：1002-1302（2020）13-0265-06

收稿日期：2019-08-13

基金項目：遼寧省水利科學技術項目（編號：20170147）。

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

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

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

1MIKE 21水動力計算原理

MIKE 21 flow model（FM，平面二維非穩態數學模型）憑借著其自身強大的卡片設置功能，干、濕節點和干、濕單元的設置功能，以及非常優秀的前、后處理功能，被廣泛地應用于二維水動力學現象的研究中。該模型常被推薦為水庫、湖泊、河流、河口和海岸水流的二維仿真模擬工具。

1.1計算原理

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

1.2計算方法

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

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

2.1數據收集與前處理

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

2.1.2數據前處理數據的提取：為了精確對石佛寺水庫庫區地表形態的描述和模型所需網格的生成，使用南方cass軟件聯合Auto CAD軟件，對石佛寺水庫地形圖內部高程點進行提取，共提取8 621個高程點信息，并利用GPS獲取庫區水面高程點，隨后利用可折疊式水深測桿在獲得高程點位置測量水深值，計算得到該點庫底高程，完成一個采樣點高程驗證，庫區共布置采樣點237個;利用Auto CAD對水庫地形圖進行邊界高程點提取，并采用同樣的方法使用GPS和測深桿對提取邊界高程數據進行采樣驗證，邊界共布置采樣點225個，完成建模高程數據的提取工作。

數據格式轉化：將所提取的高程點，采用Excel和Notepad轉換成擴展名為*.xyz的格式文件，以適合于MIKE軟件模型，供網格生成使用。所使用的高程點投影坐標為Beijing1954坐標，并以此坐標作為工作空間坐標和導入庫區高程點數據的坐標。

2.2計算網格制作

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

研究區計算網格見圖1，庫區二維地形見圖2。

2.3模型參數設置

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

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

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

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

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

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

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

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

3模型驗證與分析

3.1模型的水量平衡分析

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

3.2模型的率定與誤差分析

利用MIKE 21 FM水動力學模型對石佛寺水庫庫區進行模擬，在確定模型參數的基礎上，利用2018年8月石佛寺水庫入口斷面的流量數據進行

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

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

3.3水動力情況模擬與分析

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

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

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

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

4結論

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

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

參考文獻：

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

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

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

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

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