天然連續(xù)彎道中水流運動的三維數(shù)值模擬

呂歲菊，馮民權(quán)，李春光

(1.西安理工大學(xué) 陜西省西北旱區(qū)生態(tài)水利工程重點實驗室，陜西 西安 710048;2.北方民族大學(xué) 土木建筑工程學(xué)院，寧夏 銀川 750021)

對彎道水流運動規(guī)律的研究，在水利工程的許多領(lǐng)域占有重要的位置，研究成果在河流治理、港口興建、引水防沙以及改善河道航運等方面都得到廣泛的應(yīng)用［1］.連續(xù)彎道的水流特性影響著彎道泥沙運動和河床演變的基本規(guī)律，是河流動力學(xué)研究的主要問題之一.因此，認(rèn)識河流彎道的水流結(jié)構(gòu)，準(zhǔn)確地模擬彎道處水流特性尤其重要.

河流彎道處水流會出現(xiàn)表面流指向凹岸，底流指向凸岸的螺旋流(環(huán)流)結(jié)構(gòu)，環(huán)流運動引起的橫向輸沙，使彎道呈現(xiàn)出凸岸淤長，凹岸沖蝕的現(xiàn)象，泥沙在橫向上的輸移作用對河床沖淤非常明顯，而平面二維數(shù)學(xué)模型不能處理次生環(huán)流問題.次生環(huán)流對污染物、泥沙運動及河床變形起著關(guān)鍵作用，因此，隨著計算機速度的提高和存儲量的增大，三維數(shù)學(xué)模型得到迅速發(fā)展［2－7］.在網(wǎng)格劃分、紊流模型選取、壁面邊界條件處理、計算程序優(yōu)化等方面也取得了進步，尤其是以SIMPLE算法結(jié)合k－ε紊流模型的數(shù)值模擬發(fā)展已比較成熟.文獻［8］建立了非正交曲線坐標(biāo)下彎曲河流的標(biāo)準(zhǔn)k－ε紊流模型對彎曲河流進行了三維模擬.文獻［9］采用帶旋流修正的紊流模型對連續(xù)畸彎的水流運動進行了數(shù)值模擬研究.文獻［10］利用標(biāo)準(zhǔn)k－ε紊流模型對室內(nèi)90°彎道水流運動進行了三維數(shù)值模擬.但對天然連續(xù)彎道而言，由于受諸多因素的影響，三維水流的模擬還處于起步階段.本文在對標(biāo)準(zhǔn)k－ε模型分析的基礎(chǔ)上，考慮時均應(yīng)變率對Reynolds正應(yīng)力及紊動動能對ε方程有較大影響，基于上述原因，為了能使模擬結(jié)果更符合天然連續(xù)彎道的實際情況，對標(biāo)準(zhǔn)k－ε模型進行了改進.

本文模擬區(qū)域為天然河道，形狀極不規(guī)則，采用適體坐標(biāo)變換將笛卡爾坐標(biāo)系下三維Reynolds應(yīng)力方程變換到計算坐標(biāo)ξ－η－ζ下.針對天然彎曲河流，由于彎道螺旋流和水面比降的影響，水面沿程變化較大.首先通過實測水位插值得到計算節(jié)點上的水位，用這一水位作為計算區(qū)域的表面，采用“剛蓋”假定近似自由水面.利用有限體積法離散方程，變量布置采用同位網(wǎng)格，將u，v，w，k，ε布置在同一個網(wǎng)格節(jié)點上，并運用SIMPLEC算法進行求解.

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

1.1 三維紊流水流控制方程

Navier-Stokes方程是描述流體運動最基本的方程組，考慮水流紊流的隨機脈動，則不可壓水流運動的連續(xù)性方程和動量方程為:

式中:下標(biāo) i，j=1，2，3，遵循求和約定;u1，u2，u3分別為 x，y，z方向的 Reynolds時均流速.

1.2 適體坐標(biāo)下紊流水流模型

由于天然河流邊界曲折，地形復(fù)雜，對式(1)～(5)采用適體坐標(biāo)變換，從物理區(qū)域(x，y，z)變換到計算域(ξ，η，ζ)，其統(tǒng)一形式為:

表1 適體坐標(biāo)系下通用控制方程中的變量、系數(shù)和源項Tab.1 Variable，coefficients and source item of governing equations in body-fitted coordinate

2 數(shù)值模型

2.1 控制方程的離散

采用有限體積法對式(6)進行離散，非恒定項采用隱式差分格式，對流項采用迎風(fēng)格式，擴散項采用中心差分格式，源項采用線性化處理的方法，其離散形式為:

2.2 定解條件

固壁邊界采用無滑移條件，并由壁面函數(shù)法處理［14－15］;進口流速由實測值給定，紊動動能及紊動動能耗散率由k=0.003 75，ε =0.09k3/2給出［6］，Uin為進口斷面平均流速，取0.85 m/s;出口斷面利用充分發(fā)展條件;自由水面根據(jù)實測水位通過插值計算節(jié)點處的水位.

3 計算區(qū)域及網(wǎng)格剖分

采用上述模型對黃河寧夏沙坡頭河段連續(xù)彎道的三維水流進行編程模擬計算，計算河段長約10 km，平均河寬為230.65 m，計算流量為 739.37 m3/s，上游水位 1 240.85 m，下游水位為1 240.65 m.計算區(qū)域如圖1所示，共布設(shè)15個計算斷面.研究區(qū)域采用結(jié)構(gòu)網(wǎng)格進行剖分，劃分為124×31×15個網(wǎng)格，共計57 660個節(jié)點，時間步長取為12 s.

4 實測和計算結(jié)果對比分析

圖1 河道示意圖Fig.1 Sketch of river bends

4.1 沙坡頭河段實測數(shù)據(jù)

為了檢驗數(shù)值模型和確定初始條件，2009年7月，項目組成員利用聲學(xué)多普勒流速剖面儀、回聲測深儀等先進設(shè)備赴黃河寧夏沙坡頭河段對其中12個典型斷面的河床高程、三維流速、河寬和水位等進行了實測.表2給出了典型斷面的部分實測數(shù)據(jù).

從實測數(shù)據(jù)中可以看出，水位從上游到下游呈下降趨勢，斷面SH1較窄，水流湍急;斷面SH7較寬，水流相對較緩.

表2 典型斷面部分實測數(shù)據(jù)Tab.2 Experimental data of typical sections

4.2 平面流場分布

本文采用MATLAB語言編寫程序，對黃河寧夏沙坡頭河段天然連續(xù)彎道進行了模擬計算.

圖2給出了計算流速場的分布，從計算結(jié)果來看，在連續(xù)彎道中，一個彎頂處的平面流速分布和下一彎頂處的流速分布呈橫向?qū)ΨQ，且在兩彎頂之間開始過渡.下一彎道入口尚殘存有上一彎道未完全衰減的剩余環(huán)流，它將表層的高速水流移向凸岸，隨著下一彎道中心角的增大，反向環(huán)流的影響沿程衰減，使縱向流速逐漸過渡到均勻分布，并使反向彎道環(huán)流得到充分發(fā)展.圖3給出了兩彎頂附近處縱向垂線平均流速沿河寬的分布，可以看出，模型計算值和實測值吻合較好.

圖2 計算流場及平面流場Fig.2 Calculation velocity field and surface velocity field

圖3 計算流速和實測流速比較Fig.3 Comparison between simulated and measured velocity distributions

4.3 縱向流速沿垂線的分布

對連續(xù)彎道中典型斷面的主流流速的計算值和實測值進行對比分析，斷面SH1～SH5為第1個彎道，斷面SH6～SH9為第2個彎道，y表示距左岸的距離.

圖4 典型斷面縱向流速的垂線分布Fig.4 Vertical distribution of longitudinal velocities of the typical sections

圖4((a)～(h))為彎道斷面上5條垂線處的縱向流速分布，可以看出，模型計算結(jié)果與實測結(jié)果較為接近，均體現(xiàn)了水流在連續(xù)彎道中的分布規(guī)律.水流進入第一個彎道后，主流開始逐漸向凹岸轉(zhuǎn)移，且在彎道以下相當(dāng)長的距離內(nèi)，主流緊貼凹岸(圖4(a)～(d));當(dāng)水流進入下一個彎道時，開始過渡，過渡段的流速分布比較均勻(圖4(e)～(f));越過彎頂后，又逐漸偏向凹岸(圖4(g)～(h)).

4.4 橫向流速分布

在彎道中，由于慣性離心力的存在，水流沿軸向運動的同時，所有的水質(zhì)點又受到同樣的徑向壓力梯度，使得近水面的水質(zhì)點向凹岸運動，而近底的水質(zhì)點向凸岸運動(見圖5).由圖5可見，水流進入第一個彎段后，最大橫向環(huán)流強度出現(xiàn)在彎頂后半段(斷面SH4～SH5).中間過渡段(SH6)既是反向彎道的過渡段也是環(huán)流轉(zhuǎn)化方向的過渡段，此處環(huán)流方向已經(jīng)轉(zhuǎn)為反向，只是環(huán)流強度比較弱.下一彎道中環(huán)流的發(fā)展由于受到上一彎道的抑制作用，環(huán)流在彎頂后半部分出現(xiàn)且較明顯，斷面(SH8～SH10)橫向環(huán)流發(fā)展并出現(xiàn)加速的過程.總體上來看，本模型的計算結(jié)果是合理的，平面流速、斷面縱向垂線流速和二次流都較好地反映了實際流場的分布.

圖5 斷面橫向流速矢量場Fig.5 Transverse velocity distribution of sections

5 結(jié)語

本文對笛卡爾坐標(biāo)系下的Reynolds應(yīng)力方程進行了三維適體坐標(biāo)變換，建立了天然連續(xù)彎道河流的修正k－ε湍流模型.針對天然彎曲河流，由于環(huán)流和水面比降的影響，水面沿程變化較大.首先通過實測水位插值得到計算節(jié)點上的水位，用這一水位作為計算區(qū)域的表面，采用“剛蓋”假定近似自由水面.利用有限體積法離散控制方程，基于SIMPLEC算法編程求解，并利用黃河寧夏沙坡頭河段的實測資料進行了驗證.結(jié)果表明:模型計算的流速分布與實測資料吻合較好，說明該模型可以模擬天然河流的三維水流特性，為下一步研究該河段的泥沙沖淤和河床演變情況奠定了基礎(chǔ).計算結(jié)果還表明，自然界中河彎多以非對稱形式出現(xiàn)，非對稱河彎在兩彎頂中間存在較多低流速區(qū)，在這些區(qū)域?qū)纬赡嗌秤俾?，另外對岸流速一般較大，從而推動彎道整體以較快的速度向下游遷徙.

［1］宋志堯.羅索夫斯基彎道環(huán)流公式的統(tǒng)一［J］.水科學(xué)進展，2003，14(2):218-221.(SONG Zhi-yao.Generality of the Rozovskii's formulas on circulation at river bends［J］.Advances in Water Science，2003，14(2):218-221.(in Chinese))

［2］許棟，劉召平，乾愛國，等.彎曲河道中水流運動的三維數(shù)值模擬［J］.水利學(xué)報，2010，41(12):1423-1431.(XU Dong，LIU Zhao-ping，QIAN Ai-guo，et al.Three-dimensional numerical simulation of flow in river bends［J］.Journal of Hydraulic Engineering，2010，41(12):1423-1431.(in Chinese))

［3］ZHANG Xiao-feng，CUI Zhan-feng，LU Xin-hua.Simulation of sediment deposition in a cavity with free surface［J］.Journal of Hydrodynamics，2010，22(50):626-633.

［4］姚仕明，王興奎，張超，等.曲線同位網(wǎng)格的三維水流數(shù)學(xué)模型［J］.清華大學(xué)學(xué)報:自然科學(xué)版，2006，46(3):336-340.(YAO Shi-ming，WANG Xing-kui，ZHANG Chao，et al.3-D numerical flow model in curvilinear coordinates with collocated grid［J］.Journal of Tsinghua University(Science and Technology)，2006，46(3):336-340.(in Chinese))

［5］SHEN Yong-ming，LIU Cheng.A three-dimensional k-ε-kp model in curvilinear coordinates for sediment movement and bed evolution［J］.Sci China Ser E-Tech Sci，2009，52(4):1090-1100.

［6］陳青毅，吳衛(wèi)，劉樺.河口水庫庫內(nèi)三維流動數(shù)值模擬［J］.水動力學(xué)研究與進展:A輯，2011，26(2):140-149.(CHEN Qing-yi，WU Wei，LIU Hua.Numerical simulation of three dimensional flows in estuarine reservoir［J］.Chinese Journal of Hydrodynamics，2011，26(2):140-149.(in Chinese))

［7］易雨君，王兆印，張尚弘.考慮彎道環(huán)流影響的平面二維水沙數(shù)學(xué)模型［J］.水力發(fā)電學(xué)報，2010，29(1):126-136.(YI Yu-jun，WANG Zhao-yin，ZHANG Shang-hong.Two-dimensional sedimentation model of channel bend(Part 2)-Verification and application［J］.Journal of Hydroelectric Engineering，2010，29(1):126-136.(in Chinese))

［8］吳修廣，沈永明，潘存鴻.天然彎曲河流的三維數(shù)值模擬［J］.力學(xué)學(xué)報，2005，37(6):679-696.(WU Xiu-guang，SHEN Yong-ming，PAN Cun-hong.Three dimensional numerical simulation of natural curved river［J］.Acta Mechanica Sinica，2005，37(6):679-696.(in Chinese))

［9］李琳琳，余錫平.裁彎對上下游流態(tài)影響的三維數(shù)學(xué)模型研究［J］.水力發(fā)電學(xué)報，2010，29(5):183-189.(LI Lin-lin，YU Xi-ping.Study on effects of river cutoffs by a 3-D numerical model［J］.Journal of Hydroelectric Engineering，2010，29(5):183-189.(in Chinese))

［10］ABHARI M N，GHODSIAN M，VAGHEFI M，et al.Experimental and numerical simulation of flow in a 90°bend［J］.Flow Measurement and Instrumentation，2010，21:292-298.

［11］SHIH T H，LIOU W W，SHABBIR A，et al.A new k－ε eddy viscosity model for high Reynolds number turbulent flows［J］.Computers and Fluids，1995，24(3):227-238.

［12］ZHANG Jian-min，CHEN Jian-gang，XU Wei-lin，et al.Three-dimensional numerical simulation of aerated flows downstream sudden fall aerator expansion-in a tunnel［J］.Journal of Hydrodynamics，2011，23(1):71-80.

［13］YE Fei，GAO Xue-ping.Numerical simulations of the hydraulic characteristics of side inlet/outlets［J］.Journal of Hydrodynamics，2011，23(1):48-54.

［14］陶文銓.數(shù)值傳熱學(xué)［M］.西安:西安交通大學(xué)出版社，1995.(TAO Wen-quan.Numerical heat transfer［M］.Xi'an:Xi'an Jiaotong University Press，1995.(in Chinese))

［15］景何仿，李春光，呂歲菊，等.黃河沙坡頭連續(xù)彎道水流運動三維數(shù)值模擬［J］.寧夏大學(xué)學(xué)報:自然科學(xué)版，2010，31(3):198-202.(JING He-fang，LI Chun-guang，LV Sui-ju，et al.Three-dimensional numerical simulation of the Yellow River flow in Shapotou consective curves［J］.Journal of Ningxia University(Natural Science Edition)，2010，31(3):198-202.(in Chinese))