沙坡頭庫區三維水沙運移數值模擬研究

馮淑萍，楊 程

（1.寧夏大學研究生院，寧夏 銀川 750021；2.北方民族大學土木工程學院，寧夏 銀川 750021）

沙坡頭水利樞紐位于寧夏回族自治區中衛市境內，在黃河干流黑山峽峽谷的出口處，距離上游擬建大柳樹水利樞紐壩址12.1 km，距下游已建成的青銅峽水利樞紐122 km。沙坡頭水利樞紐是以灌溉、發電為主的大型綜合性工程。其控制流域面積25.34萬km2，多年平均徑流量為336億m3，年平均輸沙量1.6億t，總庫容0.26億m3，灌溉面積為5.85萬hm2，總裝機容量為120.3 MW，年均總發電量6.06億kW·h。沙坡頭水庫庫區屬峽谷河道型，全長14.3 km。

1 數學模型

三維水沙數學模型包括水流模塊和泥沙模塊。其中，水流模塊采用三維非恒定流可實現的 紊流模型，泥沙模塊采用全沙模型。

1.1 水流模塊

水流模塊在直角坐標系下的下形式為[1-2]

式中，Φ為通用變量，擴散系數ΓΦ和源項SΦ所代表的物理量如表1所示。

表1 直角坐標系下三維水沙模型通用方程中各變量在不同控制方程中的含義

表1中，u，v，w分別為x，y，z方向的Reynolds時均流速；p為時均動水壓強；ρ為水的密度；t為時間；k為紊動動能；ε為紊動動能的耗散率；Gk為紊動動能產生項；v為水流的分子粘性系數；vt為渦粘性系數。參數的計算見文獻[2]。

1.2 泥沙模塊

泥沙模塊控制方程在直角坐標系下的形式[3-5]:

（1）懸移質輸沙方程

式中，Cs為泥沙濃度；ω 為泥沙沉速；εx，εy，εz分別為x，y，z向的泥沙質量擴散系數，此處取εx=εy=εz=εs，εs=v+vt/σs，σs為 Schmidt數，σs=1.0。

（2）推移質不平衡輸沙方程

式中，qb、qb*分別為實際推移質輸沙率和平衡推移質輸沙率；Db、Eb分別為懸移質與推移質運動交界面上的泥沙下沉通量和上浮通量；qbx、qby分別為沿x和y方向的推移質輸沙率；Ls為推移質不平衡輸沙距離。平衡推移質輸沙率由長江科學院提出的推移質輸沙經驗公式求得。

（3）河床變形方程

式中，γs′為泥沙干容重；Z河床高程；qsx、qsy分別為通過沿水深積分得到的沿x和y方向的懸沙通量。

2 數值計算方法

2.1 模型離散

采用格心格式的非結構網格有限體積法，對方程進行離散，采用SIMPLEC算法和TDMA算法求解。水流模塊、懸移質輸沙方程可統一寫作

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

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

構造輔助點對式 （6）進行離散，瞬態項采用一階格式離散，源項采用線性化處理，具體過程可參見文獻 [7]。單元面的各通量項的離散表達式及壓力修正方程的離散參見文獻[8]。

2.2 邊界條件

兩岸邊界:按照無滑移處理，假定沒有泥沙交換。

出口邊界:根據實測數據，出口給定水位，各物理量在出口按照充分發展。

自由表面邊界:對于水面邊界采取剛蓋假定，由水位變化引起的動邊界，采取 “干濕邊界”法來處理，假定沒有泥沙交換。

3 計算實例

3.1 模擬區域及網格劃分

模擬區域——沙坡頭水庫屬于河道型水庫，庫區內共設置了20個斷面，本文選擇斷面SH1～SH11共16個斷面為研究對象，模擬區域、斷面位置及網格劃分如圖1所示。本文采用Delaunay三角化法[10-11]生成三角形網格，斷面SHJ4、SH6和SH5所在彎道處的網格進行了加密處理。

3.2 平面流場模擬結果分析

模擬區域的網格在縱向分為15層，選取表層（第15層）、中層 （第7層）和底層 （第1層）為例，對模擬的流場結果進行分析 （見圖2）。

圖1 沙坡頭庫區模擬區域及網格剖分

圖2 沙坡頭庫區三層流場

由圖2可知:①從底層到表層，流速逐漸增大。底層到中層的變化較為明顯，中層到表層流速的增長較為緩慢。②水深位置較水淺處的流速大，中間的流速較岸邊的流速大。以上兩點符合水流運動的一般規律。放大圖顯示，稀疏的網格不能精確模擬出靠近河岸邊界的流場，而網格密集的區域能精細的反映出邊界附近的流場。這也說明了數值模擬中網格對模擬結果的影響很大。

3.3 斷面流場模擬結果分析

以彎道處的橫向斷面SH6和SH5為例進行分析，3個斷面流場的等速線圖如圖3所示。由于模擬計算的初始值2008年7月的實測值，故模擬結果中各斷面是以7月份的左岸水邊點作為原點，而實測結果是2008年12月實測值，兩次實測時同一斷面的水面寬度稍有不同，故圖3中部分斷面模擬結果與實測結果的橫坐標稍有差異。

以斷面SHJ4為例，將模擬計算的斷面4等分，即得3條垂直于河面的垂線，每條線上取5個點，將斷面的流速的模擬值與實測值進行對比，如表2所示。其中，y1，y2，y3為從左岸起三條垂線，表2中對應的為垂線上點的高程值，y1，y2，y3距左岸分別為:65、130、195 m；u1為模擬值；u2為實測值。

圖3 沙坡頭庫區部分橫向斷面等流線

由圖3和表2可知:①模擬結果與實測結果較為一致，表明本文建立的三維水沙模型中的水流模塊是適合對庫區進行模擬計算的；②從河底到河面流速逐漸增大，符合水流運動的基本規律；③結合圖1中深泓線的位置，幾個斷面的最大流速與深泓線均勻不同程度的偏移，其中SH6和SH5偏移較多。由此可知，一般而言，最大水深處流速一般最大，但是在彎道的拐點、進出口、突擴或突縮等流場變化劇烈處主流與深泓線分離。

3.4 斷面垂向平均流速、主流與二次流的比較

通過上述3.2、3.3兩個部分的分析，說明流場的模擬結果符合水流運動的一般規律，并與實測結果較為一致。為了進一步驗證模擬結果與實測結果的一致性，下面以斷面SH5的主流和二次流的模擬值與實測值進行比對、分析。

將模擬計算中的斷面5等分，即得4條垂直于河面的垂線，然后將垂線上的模擬值與對應位置的實測值進行比較、分析 （見圖4）。斷面SH5的4條垂線距左岸的距離分別為47、94、141、188 m。

由圖4可知，第一，模擬值與實測值較為接近，且主流流速沿垂線的分布均符合對數律；第二，離河床床面越遠，主流流速在垂向的變化逐漸減小；第三，斷面的都是靠近水面處的橫向流速為負，靠近床面的橫向流速為正。

圖4 斷面主流、二次流比較 （定義從左岸指向右岸為正）

3.5 河床變形模擬結果分析

在二維水沙運移模擬中已對部分斷面的河床變形進行了比對。這里選擇斷面SH7和SHJ4進行簡要分析。4個斷面的河床高程對比如圖5所示。

圖5 沙坡頭庫區部分斷面河床高程對比

由圖5可知，模擬值與實測值較為一致，表明本文建立的三維水沙模型中的泥沙模塊能應用于沙坡頭庫區河床變形的模擬中。

4 結 語

以沙坡頭庫區斷面SH1-SH11段為研究對象，利用三維水沙紊流模型，采用非結構網格有限體積法、壓力修正算法，對庫區三維水沙運移進行數值模擬。對庫區三個典型層面的流場以及斷面流速分布、主流、二次流及河床高程的模擬值與實測值進行的對比、分析結果表明，兩者較為一致。得到以下結論:①從庫區底部到表層流速逐漸增大；②水深處流速一般較大；③彎道的拐點、突擴及進出口等流速變化距離的地方主流與深泓線出現分離現象；④主流流速沿垂線的分布均符合對數律；⑤離河床床面越遠，主流流速在垂向的變化逐漸減小。這些驗證了模型的準確性，及非結構網格有限體積法的適用性。

［1］陶文銓.數值傳熱學[M].2版.西安:西安交通大學出版社,2001.

［2］景何仿,李春光.黃河上游連續彎道水流運動及泥沙運移數值模擬研究[M].鄭州:黃河水利出版社,2012.

［3］假冬冬,邵學軍,王虹,等.荊江典型河灣河勢變化三維數值模型[J].水利學報,2010,41（12）:1451-1460.

［4］CANICO L,NEVES R.Three-dimensional model system for baroclinic estuarine dynamics and suspended sediment transport in a mesotidal estuary[C]//2nd International conference on Computer Modelling of Seas and Coastal Regions,Computational Mechanics Publications,1995,353-360.

［5］陳國祥,陳界仁,沙撈·巴里.三維泥沙數學模型的研究進展[J].水利水電科技進展，1998，18（1）: 13-19.

［6］唐日長.泥沙研究[M].北京:水利水電出版社，1990.

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

［8］賴錫軍，汪德貜，傅源方.非結構同位網格SIMPLE類算法收斂性能比較[J].空氣動力學報，2004，22（3）: 289-294.

［9］KANG H,CHOI S U.Reynolds stress modeling of rectangular open-channel flow[J].International Journal for Numerical Methods in Fluids,2006,51（11）: 1319-1334.

［10］陶文銓.傳熱與流動問題的多尺度數值模擬:方法與應用[M].北京:科學出版社,2009.

［11］陶文銓.計算傳熱學的近代進展[M].北京:科學出版社,2001.