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

2013-10-20 08:06:06李春光景何仿呂歲菊

水力發電 2013年5期

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

（1.寧夏大學土木與水利工程學院，寧夏銀川 750021；2.北方民族大學數值計算與工程應用研究所，寧夏銀川 750021；3.寧夏節水灌溉與水資源調控工程技術研究中心，寧夏銀川 750021；4.旱區現代農業水資源高效利用教育部工程研究中心，寧夏銀川 750021）

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

1 數學模型

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

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

式中，通用變量φ，擴散系數Γφ和源項Sφ的含義及各系數的計算公式參見文獻[3]。

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

（1）懸移質不平衡輸沙方程

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

（2）懸移質河床變形方程

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

（3）推移質河床變形方程

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

2 數值計算方法

2.1 模型離散及邊界處理

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

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

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

通過構造輔助點，利用有限體積法可對式（6）進行離散，其中瞬態項和源項的離散與結構網格一樣，即瞬態項采用向前差分，源項采用線性化處理即可。具體過程可參見文獻[9]。單元面的各通量項的離散表達式及壓力校正方程的離散參見文獻[10]。

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

2.2 計算時間步長

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

3 計算實例

3.1 工況概述

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

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

3.2 模擬區域及初始條件

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

圖1 模擬區域

初始條件的選定。由于實測資料有限，根據兩次實測資料、部分文獻以及沙坡頭水庫附近的下沿河水文站的部分測量數據，擬定數值模擬的進出口初始值。兩次實測進口斷面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。經計算分析，本文選定的進出口的流量、含沙量和水位初始值如表1所示。

3.3 數值模擬結果及分析

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

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

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