石佛寺水庫(kù)二維水動(dòng)力數(shù)值模擬

胡艷海 周林飛

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

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

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

收稿日期：2019-08-13

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

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

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

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

1MIKE 21水動(dòng)力計(jì)算原理

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

1.1計(jì)算原理

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

1.2計(jì)算方法

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

2石佛寺水庫(kù)水動(dòng)力模型的建立

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3模型驗(yàn)證與分析

3.1模型的水量平衡分析

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

3.2模型的率定與誤差分析

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

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

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

3.3水動(dòng)力情況模擬與分析

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

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

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

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

4結(jié)論

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

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

參考文獻(xiàn)：

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

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

[3]李前君，楊理?xiàng)潱瑬|愛明. 防城港市漁港工程潮流數(shù)值模擬研究[J]. 中國(guó)水運(yùn)，2016（11）：105-110.

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

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