劉 鋒,邱秀云,周 著
(1.新疆農(nóng)業(yè)大學(xué)發(fā)展規(guī)劃處,新疆烏魯木齊 830052;2.新疆農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院,新疆烏魯木齊 830052)
1995年,畢慈芬、李桂芬和于倬德針對我國砒砂巖地區(qū)水土流失對溝道徑流產(chǎn)粗沙以及暴雨洪水侵蝕單一的特點,提出采用沙棘植物“柔性壩”來攔截溝壑粗沙[1]。歷經(jīng)數(shù)年,已將其部分粗沙固于溝壑,減少了進(jìn)入黃河的粗沙,同時也使這些地區(qū)的生態(tài)環(huán)境得以改善。所謂沙棘植物“柔性壩”(以下簡稱植物壩)就是在溝道或河渠主槽內(nèi),按一定的株距和行距種植一定數(shù)量的沙棘植物。經(jīng)調(diào)查,砒砂巖在新疆普遍存在,而且因其生存環(huán)境惡劣,產(chǎn)生的粗沙相對較多,因此更迫切需要解決。
植物壩的研究為西部地區(qū),尤其為新疆地區(qū)的防沙、治沙和水土流失治理工作提供了新思路。新疆農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院水工實驗室于2002年和2003年進(jìn)行了植物壩研究的室內(nèi)水槽清水模型試驗,分別探討了植物壩在行列對齊及交錯布置方式下,不同種植密度對河渠水流的影響[1-2]。為了進(jìn)一步研究植物壩的固沙機(jī)理,筆者所在的課題組進(jìn)行了室內(nèi)水槽渾水模型試驗,并利用Fluent軟件對植物壩覆蓋區(qū)域的流場進(jìn)行了數(shù)值模擬。
試驗在新疆農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院水工實驗室的活動玻璃水槽內(nèi)進(jìn)行(圖1)。水槽長22m,寬0.4m,高0.5m,槽首設(shè)有進(jìn)水閥門和靜水池,靜水池中設(shè)消能柵板,用于穩(wěn)定水流。水槽末端設(shè)有可調(diào)節(jié)尾門,水流經(jīng)尾門后跌入水箱。水箱出口為一矩形薄壁堰,用于量測流量。水槽上設(shè)有活動測針架,可測讀水位、床面高程等。水槽頂部有帶刻度的軌道,可供測針架移動和讀取所測斷面位置。將水槽軌道零點刻度所在的斷面記為SC0+00斷面,1m刻度所在的斷面記為SC1+00,其余斷面根據(jù)所處位置的刻度依此類推。水槽中端底部設(shè)一固定鉸,首、尾端底部下分別設(shè)活動鉸各1座,用于調(diào)節(jié)水槽底坡。
試驗中,水位測量采用精度為0.1mm的水位測針量測,流速測定選用南京水利科學(xué)研究院生產(chǎn)的LGY-Ⅲ型多功能智能流速儀。
圖1 水槽布置示意圖
該研究的主要目的是為了解決新疆河流的多粗沙問題。通過野外調(diào)查,并且考慮便于研究,主要依據(jù)模型樹的高度來確定長度比尺λL=10,選用1∶10的正態(tài)模型,采用中值粒徑d50=0.34mm的非均勻天然沙為推移質(zhì)模型沙。
為了較真實地模擬天然沙棘對水沙的影響,同時考慮樹木的阻力和變形相似,該試驗采用裝飾用的綠色塑料樹作為模型樹。塑料樹全高約10cm,主干直徑約2mm,上部有樹枝和樹葉,樹頂最大橫向?qū)挾葹?cm,并具有一定的柔韌性。模型樹的根部有可用于固定的圓形孔,試驗時將模型樹固定在有插座的網(wǎng)格上。模型樹排成清水試驗中優(yōu)選出的阻水效果最好的6cm×9cm×9cm行列交錯布置方式,如圖2所示,種植棵數(shù)為63棵,種植密度為248棵/m2。
圖2 模型樹布置方式示意圖(單位:cm)
水流在植物壩內(nèi)部的流動可以看作是一種多孔介質(zhì)內(nèi)部的流動[3]。無論是從物理的角度還是從數(shù)學(xué)的角度看,多孔介質(zhì)內(nèi)部的孔隙流都是一種十分復(fù)雜的流動現(xiàn)象,本文做出如下假定:①多孔介質(zhì)為各向同性、均勻的剛性介質(zhì);②多孔介質(zhì)上部邊界為剛性可滲邊界。
很多問題中包含多孔介質(zhì)的計算,比如流場中包括過濾紙、分流器、多孔板等邊界時就需要使用多孔介質(zhì)條件。在計算中可以定義某個區(qū)域或邊界為多孔介質(zhì),并通過參數(shù)輸入,定義通過多孔介質(zhì)后流體的壓力降。
多孔介質(zhì)模型采用經(jīng)驗公式定義多孔介質(zhì)上的流動阻力。從本質(zhì)上說,多孔介質(zhì)模型就是在動量方程中增加了一個代表動量消耗的源項。因此,多孔介質(zhì)模型需要滿足下面的限制條件:①因為多孔介質(zhì)的體積在模型中沒有體現(xiàn),在缺省情況下,Fluent軟件在多孔介質(zhì)內(nèi)部使用基于體積流量的名義速度來保證速度矢量在通過多孔介質(zhì)時的連續(xù)性。如果希望更精確地進(jìn)行計算,也可以讓Fluent軟件在多孔介質(zhì)內(nèi)部使用真實速度。②多孔介質(zhì)對湍流的影響僅僅是近似。③在移動坐標(biāo)系中使用多孔介質(zhì)模型時,應(yīng)該使用相對坐標(biāo)系,而不是絕對坐標(biāo)系,以保證獲得正確的源項解。
Fluent軟件在多孔介質(zhì)計算中通過求解標(biāo)準(zhǔn)守恒型方程計算湍流變量。在計算過程中,通常假設(shè)固體介質(zhì)對湍流的生成和耗散沒有影響。在多孔介質(zhì)滲透率很大,而介質(zhì)的幾何尺度對湍流渦結(jié)構(gòu)沒有影響時,這個假設(shè)是合理的。在另外一些情況中,可能不需要考慮流動中的湍流,即假定流動為層流。如果計算中使用的湍流模型是或Spalart-Allmaras模型中的一種,則可以通過將湍流黏度 μ設(shè)為零的方式忽略湍流的影響。如果 μ被設(shè)置為零,則在計算過程中仍然會將湍流變量輸運(yùn)到介質(zhì)的另一面,但是湍流對動量輸運(yùn)過程的影響則完全被消除。在流體面板中將多孔介質(zhì)區(qū)設(shè)為層流區(qū)選項就可以將 μ設(shè)為零。
有限厚度的多孔介質(zhì)壓強(qiáng)變化是用Darcy定律和一個附加的慣性損失項來定義的:式中:μ為湍流黏度;a為滲透率;ρ為水的密度;ν為法向速度;C2為壓強(qiáng)躍升系數(shù);Δm為介質(zhì)厚度。
多孔躍升模型的相關(guān)參數(shù)在多孔躍升面板中輸入。需要輸入的參數(shù)如下:多孔躍升區(qū)域,a,Δm,C2。如果需要考慮彌散相,在多孔躍升區(qū)域定義彌散相邊界條件。
本文在Gambit里將多孔介質(zhì)設(shè)置為fluid,在fluent里設(shè)置多孔介質(zhì)參數(shù),壓強(qiáng)躍升系數(shù)為零,這樣可以保證流過多孔介質(zhì)的層流中,壓力降正比于速度,同時可以通過試算來確定a值。
網(wǎng)格生成是計算流體力學(xué)的一個重要分支,數(shù)值計算中用離散的網(wǎng)格來代替原物理問題中的連續(xù)空間,網(wǎng)格中的節(jié)點則是所求解物理量的幾何位置[4]。需要先對計算區(qū)域進(jìn)行網(wǎng)格劃分,而后在劃分的網(wǎng)格上對控制方程進(jìn)行離散、求解。
網(wǎng)格質(zhì)量對計算精度和穩(wěn)定性有很大的影響。網(wǎng)格質(zhì)量包括節(jié)點分布、光滑性以及網(wǎng)格變形程度。
本文采用較常用的Gambit軟件來生成網(wǎng)格,選用結(jié)構(gòu)網(wǎng)格對模型進(jìn)行分區(qū)網(wǎng)格劃分,這樣保證了網(wǎng)格質(zhì)量,減少了整體網(wǎng)格單元數(shù)目,同時可以較好地驗證物理試驗中以某個斷面為對象進(jìn)行測量的結(jié)果。植物壩覆蓋區(qū)段共生成網(wǎng)格單元約2800個,網(wǎng)格節(jié)點約2911個。其他區(qū)域網(wǎng)格單元約4 641個,網(wǎng)格節(jié)點約7200個。
在數(shù)值模擬計算域內(nèi)的流動是由邊界條件驅(qū)動的,從某種意義上說,求解實際問題的過程,就是將邊界線或邊界面上的數(shù)據(jù)根據(jù)控制方程外推擴(kuò)展到計算域內(nèi)部的過程。因此,提供符合物理實際且適當(dāng)?shù)倪吔鐥l件極其重要。邊界條件設(shè)置時必須保證在合適的位置選擇合適的邊界條件,同時讓邊界條件不要過約束,也不要欠約束。本文水氣兩相流場計算所涉及的邊界條件有:進(jìn)口邊界條件、出口邊界條件、固壁邊界條件和自由水面。
2.4.1進(jìn)口邊界條件
進(jìn)口邊界包括:進(jìn)水口邊界(圖3);模型頂部為空氣進(jìn)口邊界。
圖3 邊界條件示意圖(縱向剖面)
進(jìn)水口全部被水充滿,可以看做是均勻來流,采用速度進(jìn)口邊界條件,在2個物理模型中均以進(jìn)口水流的方向為 x軸正向,以豎直向上為y軸正向。邊界條件表達(dá)式為式中:ux,uy,uz分別為x,y,z方向的分速度;U0為來流速度,根據(jù)實測入流量而定。
空氣進(jìn)口邊界采用壓力進(jìn)口邊界條件,壓強(qiáng)設(shè)定為大氣壓,即p=0。全部為氣體。
2.4.2出口邊界條件
出口邊界為壓力出口和自由出流。出口由水和空氣兩部分組成,壓力出口壓強(qiáng)設(shè)定為大氣壓,即p=0;自由出流完全發(fā)展的區(qū)域:
式中:n為出口截面外法線方向;k為紊動能;ε為紊動耗散率。
2.4.3固壁邊界條件
在固壁邊界上采用無滑移邊界條件,即并利用實測的數(shù)據(jù)確定固壁糙率。
2.4.4自由水面
采用VOF方法捕捉水氣交界面,得到自由水面。
對于邊界上的紊流參數(shù)k和ε,其初始值如果不能從試驗中取得,只能采用一定的經(jīng)驗公式來估算。由式(5)(6)來估計進(jìn)口的k和ε。式中:I為底坡;u0為入口處的速度;Cμ為躍升系數(shù),可取0.09;l為湍流長度;L為關(guān)聯(lián)尺寸,對于充分發(fā)展的湍流,可取L等于水力半徑。
本文模型的數(shù)值計算過程分2步,先是模擬水槽內(nèi)無植物壩覆蓋時水流定常流流動,直至收斂,然后在這個收斂流場基礎(chǔ)上模擬水槽內(nèi)有植物壩覆蓋時的水氣兩相非定常流流動。否則如果直接模擬水氣兩相非定常流,迭代計算容易發(fā)散[5]。在此基礎(chǔ)上應(yīng)用VOF方法捕捉自由水面,變成水氣兩相非定常流模擬。數(shù)值模型的邊界條件設(shè)置:水流入口采用速度入口,水流入口水位根據(jù)試驗時的實際工況設(shè)定;入口水位以上至下游出口斷面設(shè)置為氣體壓力進(jìn)口,靜壓設(shè)為零;下游出口分為空氣出口和液態(tài)水流出口兩部分,空氣出口部分采用壓力出口,靜壓為零;各邊界條件的紊流參數(shù)均按照經(jīng)驗公式確定湍動能和紊動耗散速率,固體邊壁當(dāng)量粗糙度以試驗取0.05m。壓力插值格式選用body force weighted格式,動量、紊動能、紊動耗散速率的離散均采用一階迎風(fēng)格式[6]。時間步長為0.001s,需求時間步設(shè)定為100000,每個時間步內(nèi)的最大迭代計算次數(shù)為20次。
圖4為模擬有植物壩覆蓋水槽的3個典型斷面上的流速值與試驗值的對比。
圖4 模擬流速值與試驗值的對比
水流通過植物壩前端時,隨著植物壩阻水所產(chǎn)生的壅高水流流速急劇減小;水流流至植物壩中段時,流速趨于穩(wěn)定,沿水深縱剖面上流速值變化不大;水流通過植物壩末端時,流速受植物枝葉對水流反作用力的影響,隨著水深的增加流速值再次減小。
由圖4中3組典型斷面數(shù)值模型計算流速值與物理模型試驗值的比較可以看出,3個斷面上的模擬流速沿水深分布趨勢與試驗測得的結(jié)果基本一致。由此表明,本文利用Fluent軟件計算和VOF方法耦合的數(shù)學(xué)模型是正確的,用多孔介質(zhì)域模擬植物壩覆蓋區(qū)域是可行的。
a.水槽內(nèi)放置植物壩后,槽內(nèi)水流的流場與植樹前有了明顯的不同,其中上游壅水區(qū)是植物壩對上游水流阻擋作用的直接結(jié)果。
b.將數(shù)值計算得到的流速值與模型試驗流速值實測結(jié)果進(jìn)行了比對,兩者吻合較好。表明本文的數(shù)學(xué)模型是正確的,用多孔介質(zhì)域模擬植物壩覆蓋區(qū)域是可行的。
由數(shù)值模型計算流速值與物理模型試驗值的比較表明,Fluent軟件計算和VOF方法耦合的數(shù)學(xué)模型是正確的,用多孔介質(zhì)域模擬植物壩覆蓋區(qū)域是可行的[7]。通過較系統(tǒng)的有植物壩河道的試驗研究分析,得到了一些有價值的數(shù)據(jù)和結(jié)論,但這僅僅是初步的研究工作,還有待于深入的理論分析和采用更多的數(shù)學(xué)工具進(jìn)行更多的物理模型試驗和數(shù)值模擬方面的研究。
[1]程艷,李森,邱秀云,等.河渠種樹水流特性試驗研究[J].新疆農(nóng)業(yè)大學(xué)學(xué)報,2003,26(2):59-64.
[2]閻潔,周著,邱秀云,等.植物“柔性壩”阻水試驗及固沙機(jī)理分析[J].新疆農(nóng)業(yè)大學(xué)學(xué)報,2004,27(3):40-45.
[3]袁夢,黃本勝,邱秀云,等.水葫蘆覆蓋區(qū)水流阻力效應(yīng)試驗研究[J].廣東水利水電.2008(2):29-34.
[4]王福軍.計算流體動力學(xué)分析:CFD軟件原理與應(yīng)用[M].北京:清華大學(xué)出版社,2004.
[5]陳界仁,朱曉丹.含柔性沉水植物明渠水流運(yùn)動特性試驗研究[J].水利水電科技進(jìn)展,2008,28(3):16-19.
[6]呂升齊,唐洪武,閆靜.有無植物條件下明渠水流紊動特性對比[J].水利水電科技進(jìn)展,2007,27(6):57-60.
[7]吳航,拾兵,陳志娟.淹沒狀態(tài)下河道植物偏斜高度與糙率的關(guān)系[J].水利水電科技進(jìn)展,2008,28(1):20-22.