王靜靜，王金磊

（1.天津市水利勘測設(shè)計(jì)院，天津 300204；2.天津普澤工程咨詢有限責(zé)任公司，天津 300204）

1 MIKE21FM水動(dòng)力模型簡介

數(shù)值模擬使用的工具是由丹麥DHI公司開發(fā)的一個(gè)用于數(shù)值模擬各種流場問題（如海域、港灣、河流等）和基于流場下的環(huán)境問題（如污染物平流擴(kuò)散、水質(zhì)、重金屬、泥沙輸移）等工程問題的軟件包。本文采用MIKE21系列軟件模擬二維平面自由表面流，軟件模塊包括水動(dòng)力、水質(zhì)、泥沙、粒子追蹤等。不同模塊的模擬可以通過控制輸入文件的不同來進(jìn)行，但其中水動(dòng)力模塊是必須的，是其他模塊的基礎(chǔ)。

選用MIKE21FM模塊進(jìn)行流場及漂污物運(yùn)移規(guī)律模擬，使用以三角形為基礎(chǔ)的非結(jié)構(gòu)網(wǎng)格。在平面二維自由表面數(shù)值模擬方面具有強(qiáng)大功能。

MIKE21FM包括Hydrodynamic Module（水動(dòng)力模塊）、Transport （輸移模塊）、ECOLab （水質(zhì)模塊）、Particle Tracking（粒子追蹤模塊）、Mud Transport（黏性泥沙模塊）、Sand Transport（非黏性泥沙模塊）。MIKE21FM水動(dòng)力模型主要適合應(yīng)用在近岸沿海、海洋地區(qū)和河口、河道水流等運(yùn)移現(xiàn)象問題研究，由于這些地區(qū)的地形邊界復(fù)雜性，十分需要利用非結(jié)構(gòu)化網(wǎng)格來進(jìn)行模擬。

2 水動(dòng)力數(shù)值模型的建立

2.1 確定研究區(qū)域

該水利樞紐工程近壩河段上起廟河下至樂天溪，大壩以上壩軸線至廟河河段長17km，大壩以下壩軸線至樂天溪河段長10km，全長27km，河寬在500～1000m之間。本文研究區(qū)域選為某水利樞紐工程上游河段從美人沱至壩前，共計(jì)10km。其中美人沱至百歲溪，流程2.5km，河道微彎；百歲溪至太平溪，流程1.5km，河道彎曲明顯，南凸北凹；太平溪至壩前，流程6km，河道順直。三峽工程近壩區(qū)河段河勢如圖1。通過工具讀取AUTOCAD地形圖的地理位置坐標(biāo)和相應(yīng)的河底高程，獲得所需要的研究區(qū)域地形信息資料。研究區(qū)域河底高程如圖2。

圖1 研究區(qū)域河段河勢

圖2 研究區(qū)域河底高程

2.2 研究區(qū)域地形采集

通過工具讀取AUTOCAD地形圖的地理位置坐標(biāo)和相應(yīng)的河底高程，獲得所需要的研究區(qū)域地形信息資料。

2.3 計(jì)算網(wǎng)格劃分

網(wǎng)格構(gòu)建應(yīng)滿足以下條件[5-6]。

（1）根據(jù)水力特征、流速大小和過水能力的大小使網(wǎng)格疏密有致；

（2）構(gòu)建模型與所估計(jì)的流線平行；

（3）邊界、流場復(fù)雜區(qū)域用三角形網(wǎng)格；流速均勻、航道、湖區(qū)等地使用四邊形網(wǎng)格。

本文主要研究某水利樞紐近壩區(qū)域流場及漂污物運(yùn)移規(guī)律，采用平面二維水動(dòng)力數(shù)值模型，研究區(qū)域邊界、地形條件較復(fù)雜，綜合考慮，擬采用非結(jié)構(gòu)化三角形網(wǎng)格[7]。圖3為插值地形后的模型網(wǎng)格圖，圖4為研究區(qū)域河底三維效果圖。由圖可以看出，所生成網(wǎng)格較好地體現(xiàn)了計(jì)算區(qū)域內(nèi)的地形特征，同時(shí)，對(duì)計(jì)算區(qū)域岸線擬合也較好。

圖3 插值河底地形散點(diǎn)后網(wǎng)格劃分結(jié)果

圖4 河底三維效果

2.4 計(jì)算模型

模型控制方程為二維圣維南方程，采用有限體積法計(jì)算。

2.4.1 控制方程組

二維淺水流動(dòng)連續(xù)性方程式 （1）、x方向動(dòng)量方程式（2）和y方向動(dòng)量方程式（3）。

式中p，q分別為x，y方向上的流通通量；H為水深，H=h+ζ，其中ζ，h分別為水位和水深；w為風(fēng)速；wx，wy為風(fēng)速在x和y方向上的分量；fw為風(fēng)阻力系數(shù)；τxx，τyy，τxy為有效剪切力分量；C為謝才系數(shù)；ρ為水的密度；g為重力加速度；f為科氏力系數(shù)。

2.4.2 計(jì)算模型離散

2.4.2.1 計(jì)算區(qū)域離散

數(shù)值計(jì)算前，首先要將計(jì)算區(qū)域進(jìn)行離散，把所要計(jì)算的區(qū)域劃分成互不重疊的多個(gè)子區(qū)域，確定每個(gè)子區(qū)域中節(jié)點(diǎn)所代表的控制體積及該節(jié)點(diǎn)的位置?？紤]到果河蜿蜒，地形條件復(fù)雜，計(jì)算采有有限體積法非結(jié)構(gòu)化網(wǎng)格。

2.4.2.2 控制方程離散

有限體積法控制方程可寫成如下通用形式：

將φ取值為不同的變量，并取擴(kuò)散系數(shù)Γ和源項(xiàng)為適當(dāng)?shù)谋磉_(dá)式，可得到連續(xù)性方程、動(dòng)量方程、能量方程、紊動(dòng)能方程和紊動(dòng)耗散率方程。

對(duì)二維問題，離散方程可表示為：

其中，aw=Dw+max（0，F(xiàn)w）；

式中ΔV是控制體積的體積；Sp為隨時(shí)間和物理量φ變化的項(xiàng)；Sc是常數(shù)；F為通過界面上單位面積的對(duì)流質(zhì)量通量，簡稱對(duì)流質(zhì)量通量；D為界面的擴(kuò)散傳導(dǎo)性；φ為溫度、速度、濃度等一些待求的物理量，是廣義變量。

3 計(jì)算結(jié)果及驗(yàn)證

模型計(jì)算采用某工程右岸地下電廠建成并投入運(yùn)行后工況條件，為了與物理試驗(yàn)計(jì)算結(jié)果進(jìn)行對(duì)比，設(shè)計(jì)4種計(jì)算工況，工況S1與S2水位相同，為145m，流量分別為3.5萬，4.5萬m3/s。2種工況下深孔泄洪，排漂孔、洞關(guān)閉。選工況S1與工況S2的實(shí)測流速進(jìn)行模擬結(jié)果驗(yàn)證。

取模型試驗(yàn)壩前500m斷面cs1和壩前1500m斷面cs2。cs1斷面與cs2斷面實(shí)測流速測點(diǎn)間距100m。cs1斷面和cs2斷面實(shí)測流速與數(shù)值模擬流速值擬合如圖5～圖8。

圖5 工況S1斷面cs1擬合

圖6 工況S1斷面cs2擬合

圖7 工況S2斷面cs1擬合

圖8 工況S2斷面cs2擬合

由于物理模型所測為河流表面的流速大小，本文計(jì)算的流速是沿水深方向的垂線平均流速，因此用于驗(yàn)證的實(shí)測數(shù)據(jù)是物理模型數(shù)據(jù)乘以一個(gè)沿水深進(jìn)行積分的系數(shù)，此系數(shù)取經(jīng)驗(yàn)值0.85。分析對(duì)比兩種工況下cs1斷面與cs2斷面實(shí)測值與模擬值對(duì)比結(jié)果可知：數(shù)值模擬結(jié)果多數(shù)測點(diǎn)兩者流速趨勢較為一致，與實(shí)測結(jié)果吻合較好，誤差可視為在合理解釋范圍內(nèi)。

通過數(shù)值模擬所得到的某水利樞紐近壩區(qū)域流場計(jì)算結(jié)果如圖9～圖12。分別給出工況S1至工況S2條件下的水位流場分布圖、流速分布圖。圖中可以看出，工況S1條件下，壩前區(qū)最大流速集中分布在中間偏右側(cè)部分，流速集中在0.4~0.9m/s之間，水流在壩前左右兩側(cè)形成兩個(gè)回流區(qū)；右岸回流區(qū)由于右電廠與地下電廠之間隔有條形山體，故回流范圍較左岸更大，壩前流態(tài)較不穩(wěn)定。工況S2條件下，壩前最大流速仍集中在中間偏右部分，流速集中在0.6~1.2m/s之間，壩前流態(tài)與工況S1大致相似。

圖9 工況S1流場分布

圖10 工況S1流速分布

圖11 工況S2流場分布

圖12 工況S2流速分布

4 結(jié)語

（1）對(duì)某水利樞紐近壩區(qū)域在給定工況下的流態(tài)及流速進(jìn)行了驗(yàn)證，計(jì)算所得流場與實(shí)際觀測流場基本一致，計(jì)算所得流速與實(shí)際測量流速數(shù)據(jù)基本吻合。

（2）模擬流場顯示在兩種計(jì)算工況下，近壩區(qū)流場變化不大：①美人沱至百歲溪河段，右岸微凸，左岸微凹，水流較順暢，在兩岸山坳有小范圍緩流、回流產(chǎn)生；②百歲溪到太平溪段，河段彎曲明顯，主流居中央；③太平溪至引航道口門處，河道逐漸順直變寬，主流略偏左，至九嶺山段由于左岸凹陷，此處出現(xiàn)緩流、回流區(qū)，并使主流右偏，流態(tài)較不穩(wěn)定。

（3）工況S1條件下，壩前區(qū)域最大流速區(qū)集中分布在中間偏右側(cè)部分，流速值在0.4~0.9m/s之間，水流在壩前左右兩側(cè)形成兩個(gè)回流區(qū)；右岸回流區(qū)由于右電廠與地下電廠之間隔有條形山體，故回流范圍較左岸更大，壩前流態(tài)較不穩(wěn)定。工況S2條件下，壩前最大流速仍集中在中間偏右部分，流速值集中在0.6~1.2m/s之間，壩前流態(tài)與工況S1大致相似。模擬流場與物理模型試驗(yàn)流態(tài)對(duì)比結(jié)果較一致，流場模擬較好。

[1]黃鵬，江帆.Fluent高級(jí)應(yīng)用與實(shí)例分析[M].北京： 清華大學(xué)出版社，2008.

[2]謝漢龍，朱紅鈞，林元華.流體分析及仿真實(shí)用教程[M].北京：人民郵電出版社，2010.

[3]馮民權(quán)，鄭邦民，周孝德.河流及水庫流場與水質(zhì)的數(shù)值模擬[M].北京： 科學(xué)出版社，2007.

[4]鄧義祥，張愛軍.環(huán)境水動(dòng)力數(shù)值模擬的研究現(xiàn)狀與展望[J].資源開發(fā)與市場，1998（6）：276-278.

[5]文岑.平面二維水沙運(yùn)動(dòng)與船舶運(yùn)動(dòng)的數(shù)值模擬研究及其應(yīng)用[D].重慶：重慶大學(xué)，2012.

[6]Haliva-Cohen，A.Sources and transport routes of fine detritusmaterial to the Late Quaternary Dead Sea basin[J].Quaternary Science Reviews.2012，50（4）.

[7]施奇.SMS二維水動(dòng)力數(shù)學(xué)模型在淮干方邱湖河段整治中的應(yīng)用[J].江淮水利科技，2011（3）：11-13.