寧夏沙坡頭河段二維水沙運動數(shù)值模擬研究

楊 程 ，李春光 ，景何仿 ，呂歲菊

（1.寧夏大學(xué)土木與水利工程學(xué)院，寧夏 銀川 750021；2.北方民族大學(xué)數(shù)值計算與工程應(yīng)用研究所，寧夏 銀川 750021；3.寧夏節(jié)水灌溉與水資源調(diào)控工程技術(shù)研究中心，寧夏 銀川 750021；4.旱區(qū)現(xiàn)代農(nóng)業(yè)水資源高效利用教育部工程研究中心，寧夏 銀川 750021）

多數(shù)天然河流、海岸、水庫、內(nèi)河等水域的水平尺度遠(yuǎn)大于垂直尺度，因此沿水深平均的平面二維水沙數(shù)學(xué)模型是目前研究河床變形使用較為廣泛的數(shù)學(xué)模型[1]。在運用離散數(shù)學(xué)模型之前必須先對計算區(qū)域進(jìn)行網(wǎng)格劃分，而網(wǎng)格質(zhì)量的好壞將直接影響到數(shù)值計算的可行性、收斂的速度以及計算結(jié)果的精度。目前的網(wǎng)格一般分為結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格。對于復(fù)雜邊界的流動問題，結(jié)構(gòu)網(wǎng)格為了能實現(xiàn)其貼體性，會使邊界處的網(wǎng)格不能正交或減弱其正交性，從而影響網(wǎng)格的質(zhì)量。而非結(jié)構(gòu)網(wǎng)格的幾何特性更加靈活。特別是它能更為精確和有效地表示復(fù)雜的流場邊界，并根據(jù)流體計算的精度需要調(diào)整網(wǎng)格單元的大小和密度分布[2]。故本文以沙坡頭河段為例，選用非結(jié)構(gòu)網(wǎng)格對計算區(qū)域進(jìn)行離散，研究其二維水沙運動，為進(jìn)一步研究水庫淤積規(guī)律奠定基礎(chǔ)，為水庫管理運行提供一定的參考依據(jù)。

1 數(shù)學(xué)模型

沿水深平均的二維水沙數(shù)學(xué)模型包括兩個模塊:水流模塊和泥沙模塊。本文選擇二維RNG紊流水沙數(shù)學(xué)模型。

水流模塊包括[3]:連續(xù)性方程、x方向動量方程、y方向動量方程、k-方程和ε-方程，其通用方程為

式中，通用變量φ，擴(kuò)散系數(shù)Γφ和源項Sφ的含義及各系數(shù)的計算公式參見文獻(xiàn)[3]。

泥沙模塊包括[4-5]懸移質(zhì)不平衡輸沙方程、推移質(zhì)不平衡輸沙方程和河床變形方程。

（1）懸移質(zhì)不平衡輸沙方程

式中，S為垂線平均含沙量；εx，εy分別為x，y方向的泥沙擴(kuò)散系數(shù)，此處取εx=εy=νt；ω為沉速，由張瑞瑾公式計算；α為垂線平均恢復(fù)飽和系數(shù)，由武漢水利電力學(xué)院韋直林公式[6]得；S*為水流挾沙力，由張紅武公式[7]計算。

（2）懸移質(zhì)河床變形方程

式中，γ′s為泥沙干密度；Zs為懸移質(zhì)泥沙引起的河床沖淤厚度。

（3）推移質(zhì)河床變形方程

式中，Zb，gbx，gby分別是推移質(zhì)泥沙引起的河床沖淤厚度、x方向和y方向的單寬輸沙率，其中單寬輸沙率由竇國仁公式[8]計算。

2 數(shù)值計算方法

2.1 模型離散及邊界處理

模型離散。采用格心格式的中心有限體積法，將變量存在單元的中心，單元的邊界為控制體，則

對式（5）在控制單元上求體積分并利用高斯定理將體積分化為面積分

式中，U為網(wǎng)格單元的流速矢量；n為網(wǎng)格界面上的法向量。

通過構(gòu)造輔助點，利用有限體積法可對式 （6）進(jìn)行離散，其中瞬態(tài)項和源項的離散與結(jié)構(gòu)網(wǎng)格一樣，即瞬態(tài)項采用向前差分，源項采用線性化處理即可。具體過程可參見文獻(xiàn)[9]。單元面的各通量項的離散表達(dá)式及壓力校正方程的離散參見文獻(xiàn)[10]。

式（3）、式（4）的離散:式（3）等式左邊的瞬態(tài)項可采用向前差分處理，右邊的源項采用線性化處理；式（4）的第一項為瞬態(tài)項采用向前差分處理，第二、三項與處對流項的方式相同。

2.2 計算時間步長

為了保證計算時模型的穩(wěn)定，時間步長的選取需要滿足CFL線性穩(wěn)定條件[12]

3 計算實例

3.1 工況概述

沙坡頭水利樞紐位于寧夏中衛(wèi)縣境內(nèi)的黃河干流黑山峽峽谷出口處，是以灌溉、發(fā)電為主的大型綜合性水利樞紐，控制流域面積25.34萬km2，多年平均徑流量為336億m3，年平均輸沙量1.6億t，總庫容0.26億m3，灌溉面積為5.847萬hm2，總裝機(jī)容量為120.3 MW，年均總發(fā)電量6.06億kW·h。沙坡頭水庫屬峽谷河道型，全長14.3 km。

項目組在2008年的7月16日～18日和12月5日～7日分別對沙坡頭的SH1-SH11斷面進(jìn)行實測，得到斷面的水深、流速、流量等相關(guān)數(shù)據(jù)。本文以2008年7月實測數(shù)據(jù)作為數(shù)值模擬的初始值，模擬計算時間為2008年7月16日至12月7日。

3.2 模擬區(qū)域及初始條件

模擬區(qū)域。所選模擬河段共設(shè)置了20個斷面，本文選擇16個斷面為研究對象，模擬區(qū)域及斷面位置如圖1所示，其中斷面SH11是進(jìn)水口，斷面SH1是出水口。本文采用前沿推進(jìn)法[13-14]在計算區(qū)域內(nèi)生成三角形網(wǎng)格。網(wǎng)格劃分共5103個節(jié)點，9406個單元。

圖1 模擬區(qū)域

初始條件的選定。由于實測資料有限，根據(jù)兩次實測資料、部分文獻(xiàn)以及沙坡頭水庫附近的下沿河水文站的部分測量數(shù)據(jù)，擬定數(shù)值模擬的進(jìn)出口初始值。兩次實測進(jìn)口斷面SH11的流量分別為:1034.50 m3/s和588.71 m3/s。沙坡頭水庫汛期一般為7月～9月，入庫徑流主要來自蘭州以上,多年平均流量820 m3/s，多年平均含沙量3.47 kg/m3[15]。汛期水量占全年的44%，輸沙量占全年的82.7%[16]。兩次實測出口斷面SH1的水位分別為1240.50 m和1240.51 m，水庫正常運行水位1240.50 m。經(jīng)計算分析，本文選定的進(jìn)出口的流量、含沙量和水位初始值如表1所示。

3.3 數(shù)值模擬結(jié)果及分析

模擬區(qū)域有兩個彎道，根據(jù)水的流向，將第一個彎道記作A彎道，第二個彎道記作B彎道。選擇A彎道中斷面SH9和SHJ5以及B彎道附近的斷面SH5和SH6為例，各斷面位置如圖1所示。將斷面垂向流速和河底高程模擬值與實測值進(jìn)行對比并分析。模擬區(qū)域的流場圖及局部放大圖如圖2所示。

斷面SH5、SH6、SH9和SHJ5斷面垂向平均流速和河底高程的模擬值與實測值的對比如圖3～圖6所示。

由圖3～圖6可知:①模擬值與實測值較為一致，說明本文建立的模型、采用的方法以及初始條件的選定等是合理的；②斷面SH5所處位置可看作是B彎道的出口處，受到彎道離心力的影響，斷面的左岸附近有稍許沖刷，右岸淤積比較嚴(yán)重；③斷面SH6處于彎道中，左岸沖刷、右岸淤積較為明顯；④斷面SHJ5處于A彎道中，右岸沖刷左岸淤積較為明顯；⑤由于斷面SH9垂向平均流速的模擬值稍小于實測值，導(dǎo)致該斷面泥沙淤積量的模擬值稍大于實測值。

圖2 模擬區(qū)域流場

表1 進(jìn)出口初始值

圖3 斷面SH5

圖4 斷面SH6

4 結(jié) 論

（1）本文建立的平面二維水沙模型是合理的，能夠較為精確模擬沙坡頭庫區(qū)二維水沙運動。

（2）研究所選斷面垂線平均流速、庫底變形的模擬值和實測值較為一致，說明非結(jié)構(gòu)網(wǎng)格有限體積法是適用的。

圖5 斷面SHJ5

圖6 斷面SH9

（3）研究結(jié)果為后續(xù)研究、預(yù)測該水庫的泥沙淤積情況及多水庫聯(lián)合調(diào)度奠定了必要的基礎(chǔ)。

［1］景何仿，李春光.黃河上游大柳樹——沙坡頭河段河床變形數(shù)值模擬［J］.中國農(nóng)村水利水電， 2011（11）:5-9， 13.

［2］尹河.在非結(jié)構(gòu)自適應(yīng)網(wǎng)格上對二維Euler方程進(jìn)行數(shù)值模擬［D］.西安:西北工業(yè)大學(xué)博士學(xué)位論文，2001.

［3］吳修廣.曲線坐標(biāo)系下水流和污染物擴(kuò)散輸移的湍流模型［D］.大連:大連理工大學(xué)博士學(xué)位論文，2004.

［4］竇希萍，李來，竇國仁.長江口全沙數(shù)學(xué)模型研究［J］.水利水運科學(xué)研究， 1999（2）:136-145.

［5］李春光，景何仿，呂歲菊，等.大柳樹-沙坡頭河段泥沙運移二維數(shù)值模擬［J］.水利水運工程學(xué)報， 2011（4）:102-107.

［6］錢意穎，曲少軍，曹文洪，等.黃河泥沙沖淤數(shù)學(xué)模型［M］.鄭州:黃河水利出版社，1998.

［7］張俊華，張紅武，王嚴(yán)平，等.多沙水庫準(zhǔn)二維泥沙數(shù)學(xué)模型［J］.水動力學(xué)研究與進(jìn)展， 1993， 14（1）:45-50.

［8］竇國仁，趙世清，黃亦芬.河道二維全沙數(shù)學(xué)模型的研究［J］.水利水運科學(xué)研究，1987（2）:1-12.

［9］LAI X J,WANG D G,CHEN Y.Pressure correction method on unstructured Grids［J］.Hydrodynamics， 2004， 16（3）:316-324.

［10］賴錫軍，汪德貜，傅源方.非結(jié)構(gòu)同位網(wǎng)格SIMPLE類算法收斂性能比較［J］.空氣動力學(xué)報， 2004， 22（3）:289-294.

［11］馬方凱.河口三維水沙輸移過程數(shù)值模擬研究［D］.北京:清華大學(xué)博士學(xué)位論文，2010.

［12］岳志遠(yuǎn)，曹志先，李有為，等.基于非結(jié)構(gòu)網(wǎng)格的非恒定淺水二維有限體積數(shù)學(xué)模型研究［J］.水動力學(xué)研究與進(jìn)展，2011，26 （3）:359-367.

［13］陶文銓.計算傳熱學(xué)的近代進(jìn)展［M］.北京:科學(xué)出版社，2001:67-71.

［14］孫力勝，鄭建靖，陳建軍，等.二維自適應(yīng)前沿推進(jìn)網(wǎng)格生成［J］.計算機(jī)工程與應(yīng)用， 2011， 47（3）:146-148.

［15］丁義斌，姚志霞，謝建勇，等.沙坡頭水庫淤積現(xiàn)狀分析與評價［M］.水利科技與經(jīng)濟(jì)， 2009， 15（10）:915-916.

［16］鄭賀新.黃河沙坡頭水庫泥沙沖淤分析［M］.水利水電工程設(shè)計， 1998（3）:28-29， 34， 56.