沙坡頭庫區(qū)三維水沙運(yùn)移數(shù)值模擬研究

馮淑萍，楊 程

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

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

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

三維水沙數(shù)學(xué)模型包括水流模塊和泥沙模塊。其中，水流模塊采用三維非恒定流可實(shí)現(xiàn)的 紊流模型，泥沙模塊采用全沙模型。

1.1 水流模塊

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

式中，Φ為通用變量，擴(kuò)散系數(shù)ΓΦ和源項(xiàng)SΦ所代表的物理量如表1所示。

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

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

1.2 泥沙模塊

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

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

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

（2）推移質(zhì)不平衡輸沙方程

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

（3）河床變形方程

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

2 數(shù)值計算方法

2.1 模型離散

采用格心格式的非結(jié)構(gòu)網(wǎng)格有限體積法，對方程進(jìn)行離散，采用SIMPLEC算法和TDMA算法求解。水流模塊、懸移質(zhì)輸沙方程可統(tǒng)一寫作

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

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

構(gòu)造輔助點(diǎn)對式 （6）進(jìn)行離散，瞬態(tài)項(xiàng)采用一階格式離散，源項(xiàng)采用線性化處理，具體過程可參見文獻(xiàn) [7]。單元面的各通量項(xiàng)的離散表達(dá)式及壓力修正方程的離散參見文獻(xiàn)[8]。

2.2 邊界條件

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

出口邊界:根據(jù)實(shí)測數(shù)據(jù)，出口給定水位，各物理量在出口按照充分發(fā)展。

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

3 計算實(shí)例

3.1 模擬區(qū)域及網(wǎng)格劃分

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

3.2 平面流場模擬結(jié)果分析

模擬區(qū)域的網(wǎng)格在縱向分為15層，選取表層（第15層）、中層 （第7層）和底層 （第1層）為例，對模擬的流場結(jié)果進(jìn)行分析 （見圖2）。

圖1 沙坡頭庫區(qū)模擬區(qū)域及網(wǎng)格剖分

圖2 沙坡頭庫區(qū)三層流場

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

3.3 斷面流場模擬結(jié)果分析

以彎道處的橫向斷面SH6和SH5為例進(jìn)行分析，3個斷面流場的等速線圖如圖3所示。由于模擬計算的初始值2008年7月的實(shí)測值，故模擬結(jié)果中各斷面是以7月份的左岸水邊點(diǎn)作為原點(diǎn)，而實(shí)測結(jié)果是2008年12月實(shí)測值，兩次實(shí)測時同一斷面的水面寬度稍有不同，故圖3中部分?jǐn)嗝婺M結(jié)果與實(shí)測結(jié)果的橫坐標(biāo)稍有差異。

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

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

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

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

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

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

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

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

3.5 河床變形模擬結(jié)果分析

在二維水沙運(yùn)移模擬中已對部分?jǐn)嗝娴暮哟沧冃芜M(jìn)行了比對。這里選擇斷面SH7和SHJ4進(jìn)行簡要分析。4個斷面的河床高程對比如圖5所示。

圖5 沙坡頭庫區(qū)部分?jǐn)嗝婧哟哺叱虒Ρ?/p>

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

4 結(jié) 語

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

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

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

［3］假冬冬,邵學(xué)軍,王虹,等.荊江典型河灣河勢變化三維數(shù)值模型[J].水利學(xué)報,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］陳國祥,陳界仁,沙撈·巴里.三維泥沙數(shù)學(xué)模型的研究進(jìn)展[J].水利水電科技進(jìn)展，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］賴錫軍，汪德貜，傅源方.非結(jié)構(gòu)同位網(wǎng)格SIMPLE類算法收斂性能比較[J].空氣動力學(xué)報，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］陶文銓.傳熱與流動問題的多尺度數(shù)值模擬:方法與應(yīng)用[M].北京:科學(xué)出版社,2009.

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