劉 鋒,邱秀云,周 著

(1.新疆農(nóng)業(yè)大學(xué)發(fā)展規(guī)劃處,新疆烏魯木齊 830052;2.新疆農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院,新疆烏魯木齊 830052)

1995年,畢慈芬、李桂芬和于倬德針對(duì)我國(guó)砒砂巖地區(qū)水土流失對(duì)溝道徑流產(chǎn)粗沙以及暴雨洪水侵蝕單一的特點(diǎn),提出采用沙棘植物“柔性壩”來(lái)攔截溝壑粗沙[1]。歷經(jīng)數(shù)年,已將其部分粗沙固于溝壑,減少了進(jìn)入黃河的粗沙,同時(shí)也使這些地區(qū)的生態(tài)環(huán)境得以改善。所謂沙棘植物“柔性壩”(以下簡(jiǎn)稱植物壩)就是在溝道或河渠主槽內(nèi),按一定的株距和行距種植一定數(shù)量的沙棘植物。經(jīng)調(diào)查,砒砂巖在新疆普遍存在,而且因其生存環(huán)境惡劣,產(chǎn)生的粗沙相對(duì)較多,因此更迫切需要解決。

植物壩的研究為西部地區(qū),尤其為新疆地區(qū)的防沙、治沙和水土流失治理工作提供了新思路。新疆農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院水工實(shí)驗(yàn)室于2002年和2003年進(jìn)行了植物壩研究的室內(nèi)水槽清水模型試驗(yàn),分別探討了植物壩在行列對(duì)齊及交錯(cuò)布置方式下,不同種植密度對(duì)河渠水流的影響[1-2]。為了進(jìn)一步研究植物壩的固沙機(jī)理,筆者所在的課題組進(jìn)行了室內(nèi)水槽渾水模型試驗(yàn),并利用Fluent軟件對(duì)植物壩覆蓋區(qū)域的流場(chǎng)進(jìn)行了數(shù)值模擬。

1 物理試驗(yàn)簡(jiǎn)介

1.1 試驗(yàn)設(shè)備

試驗(yàn)在新疆農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院水工實(shí)驗(yàn)室的活動(dòng)玻璃水槽內(nèi)進(jìn)行(圖1)。水槽長(zhǎng)22m,寬0.4m,高0.5m,槽首設(shè)有進(jìn)水閥門和靜水池,靜水池中設(shè)消能柵板,用于穩(wěn)定水流。水槽末端設(shè)有可調(diào)節(jié)尾門,水流經(jīng)尾門后跌入水箱。水箱出口為一矩形薄壁堰,用于量測(cè)流量。水槽上設(shè)有活動(dòng)測(cè)針架,可測(cè)讀水位、床面高程等。水槽頂部有帶刻度的軌道,可供測(cè)針架移動(dòng)和讀取所測(cè)斷面位置。將水槽軌道零點(diǎn)刻度所在的斷面記為SC0+00斷面,1m刻度所在的斷面記為SC1+00,其余斷面根據(jù)所處位置的刻度依此類推。水槽中端底部設(shè)一固定鉸,首、尾端底部下分別設(shè)活動(dòng)鉸各1座,用于調(diào)節(jié)水槽底坡。

試驗(yàn)中,水位測(cè)量采用精度為0.1mm的水位測(cè)針量測(cè),流速測(cè)定選用南京水利科學(xué)研究院生產(chǎn)的LGY-Ⅲ型多功能智能流速儀。

圖1 水槽布置示意圖

1.2 模型比尺

該研究的主要目的是為了解決新疆河流的多粗沙問(wèn)題。通過(guò)野外調(diào)查,并且考慮便于研究,主要依據(jù)模型樹的高度來(lái)確定長(zhǎng)度比尺λL=10,選用1∶10的正態(tài)模型,采用中值粒徑d50=0.34mm的非均勻天然沙為推移質(zhì)模型沙。

1.3 模型樹的選擇

為了較真實(shí)地模擬天然沙棘對(duì)水沙的影響,同時(shí)考慮樹木的阻力和變形相似,該試驗(yàn)采用裝飾用的綠色塑料樹作為模型樹。塑料樹全高約10cm,主干直徑約2mm,上部有樹枝和樹葉,樹頂最大橫向?qū)挾葹?cm,并具有一定的柔韌性。模型樹的根部有可用于固定的圓形孔,試驗(yàn)時(shí)將模型樹固定在有插座的網(wǎng)格上。模型樹排成清水試驗(yàn)中優(yōu)選出的阻水效果最好的6cm×9cm×9cm行列交錯(cuò)布置方式,如圖2所示,種植棵數(shù)為63棵,種植密度為248棵/m2。

圖2 模型樹布置方式示意圖(單位:cm)

2 流場(chǎng)的數(shù)值模擬

水流在植物壩內(nèi)部的流動(dòng)可以看作是一種多孔介質(zhì)內(nèi)部的流動(dòng)[3]。無(wú)論是從物理的角度還是從數(shù)學(xué)的角度看,多孔介質(zhì)內(nèi)部的孔隙流都是一種十分復(fù)雜的流動(dòng)現(xiàn)象,本文做出如下假定:①多孔介質(zhì)為各向同性、均勻的剛性介質(zhì);②多孔介質(zhì)上部邊界為剛性可滲邊界。

很多問(wèn)題中包含多孔介質(zhì)的計(jì)算,比如流場(chǎng)中包括過(guò)濾紙、分流器、多孔板等邊界時(shí)就需要使用多孔介質(zhì)條件。在計(jì)算中可以定義某個(gè)區(qū)域或邊界為多孔介質(zhì),并通過(guò)參數(shù)輸入,定義通過(guò)多孔介質(zhì)后流體的壓力降。

2.1 多孔介質(zhì)模型的假設(shè)和限制條件

多孔介質(zhì)模型采用經(jīng)驗(yàn)公式定義多孔介質(zhì)上的流動(dòng)阻力。從本質(zhì)上說(shuō),多孔介質(zhì)模型就是在動(dòng)量方程中增加了一個(gè)代表動(dòng)量消耗的源項(xiàng)。因此,多孔介質(zhì)模型需要滿足下面的限制條件:①因?yàn)槎嗫捉橘|(zhì)的體積在模型中沒(méi)有體現(xiàn),在缺省情況下,Fluent軟件在多孔介質(zhì)內(nèi)部使用基于體積流量的名義速度來(lái)保證速度矢量在通過(guò)多孔介質(zhì)時(shí)的連續(xù)性。如果希望更精確地進(jìn)行計(jì)算,也可以讓Fluent軟件在多孔介質(zhì)內(nèi)部使用真實(shí)速度。②多孔介質(zhì)對(duì)湍流的影響僅僅是近似。③在移動(dòng)坐標(biāo)系中使用多孔介質(zhì)模型時(shí),應(yīng)該使用相對(duì)坐標(biāo)系,而不是絕對(duì)坐標(biāo)系,以保證獲得正確的源項(xiàng)解。

2.2 多孔介質(zhì)模型對(duì)湍流的處理

Fluent軟件在多孔介質(zhì)計(jì)算中通過(guò)求解標(biāo)準(zhǔn)守恒型方程計(jì)算湍流變量。在計(jì)算過(guò)程中,通常假設(shè)固體介質(zhì)對(duì)湍流的生成和耗散沒(méi)有影響。在多孔介質(zhì)滲透率很大,而介質(zhì)的幾何尺度對(duì)湍流渦結(jié)構(gòu)沒(méi)有影響時(shí),這個(gè)假設(shè)是合理的。在另外一些情況中,可能不需要考慮流動(dòng)中的湍流,即假定流動(dòng)為層流。如果計(jì)算中使用的湍流模型是或Spalart-Allmaras模型中的一種,則可以通過(guò)將湍流黏度 μ設(shè)為零的方式忽略湍流的影響。如果 μ被設(shè)置為零,則在計(jì)算過(guò)程中仍然會(huì)將湍流變量輸運(yùn)到介質(zhì)的另一面,但是湍流對(duì)動(dòng)量輸運(yùn)過(guò)程的影響則完全被消除。在流體面板中將多孔介質(zhì)區(qū)設(shè)為層流區(qū)選項(xiàng)就可以將 μ設(shè)為零。

有限厚度的多孔介質(zhì)壓強(qiáng)變化是用Darcy定律和一個(gè)附加的慣性損失項(xiàng)來(lái)定義的:式中:μ為湍流黏度;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ù)為零,這樣可以保證流過(guò)多孔介質(zhì)的層流中,壓力降正比于速度,同時(shí)可以通過(guò)試算來(lái)確定a值。

2.3 網(wǎng)格剖分

網(wǎng)格生成是計(jì)算流體力學(xué)的一個(gè)重要分支,數(shù)值計(jì)算中用離散的網(wǎng)格來(lái)代替原物理問(wèn)題中的連續(xù)空間,網(wǎng)格中的節(jié)點(diǎn)則是所求解物理量的幾何位置[4]。需要先對(duì)計(jì)算區(qū)域進(jìn)行網(wǎng)格劃分,而后在劃分的網(wǎng)格上對(duì)控制方程進(jìn)行離散、求解。

網(wǎng)格質(zhì)量對(duì)計(jì)算精度和穩(wěn)定性有很大的影響。網(wǎng)格質(zhì)量包括節(jié)點(diǎn)分布、光滑性以及網(wǎng)格變形程度。

本文采用較常用的Gambit軟件來(lái)生成網(wǎng)格,選用結(jié)構(gòu)網(wǎng)格對(duì)模型進(jìn)行分區(qū)網(wǎng)格劃分,這樣保證了網(wǎng)格質(zhì)量,減少了整體網(wǎng)格單元數(shù)目,同時(shí)可以較好地驗(yàn)證物理試驗(yàn)中以某個(gè)斷面為對(duì)象進(jìn)行測(cè)量的結(jié)果。植物壩覆蓋區(qū)段共生成網(wǎng)格單元約2800個(gè),網(wǎng)格節(jié)點(diǎn)約2911個(gè)。其他區(qū)域網(wǎng)格單元約4 641個(gè),網(wǎng)格節(jié)點(diǎn)約7200個(gè)。

2.4 邊界條件的設(shè)置

在數(shù)值模擬計(jì)算域內(nèi)的流動(dòng)是由邊界條件驅(qū)動(dòng)的,從某種意義上說(shuō),求解實(shí)際問(wèn)題的過(guò)程,就是將邊界線或邊界面上的數(shù)據(jù)根據(jù)控制方程外推擴(kuò)展到計(jì)算域內(nèi)部的過(guò)程。因此,提供符合物理實(shí)際且適當(dāng)?shù)倪吔鐥l件極其重要。邊界條件設(shè)置時(shí)必須保證在合適的位置選擇合適的邊界條件,同時(shí)讓邊界條件不要過(guò)約束,也不要欠約束。本文水氣兩相流場(chǎng)計(jì)算所涉及的邊界條件有:進(jìn)口邊界條件、出口邊界條件、固壁邊界條件和自由水面。

2.4.1進(jìn)口邊界條件

進(jìn)口邊界包括:進(jìn)水口邊界(圖3);模型頂部為空氣進(jìn)口邊界。

圖3 邊界條件示意圖(縱向剖面)

進(jìn)水口全部被水充滿,可以看做是均勻來(lái)流,采用速度進(jìn)口邊界條件,在2個(gè)物理模型中均以進(jìn)口水流的方向?yàn)?x軸正向,以豎直向上為y軸正向。邊界條件表達(dá)式為式中:ux,uy,uz分別為x,y,z方向的分速度;U0為來(lái)流速度,根據(jù)實(shí)測(cè)入流量而定。

空氣進(jìn)口邊界采用壓力進(jìn)口邊界條件,壓強(qiáng)設(shè)定為大氣壓,即p=0。全部為氣體。

2.4.2出口邊界條件

出口邊界為壓力出口和自由出流。出口由水和空氣兩部分組成,壓力出口壓強(qiáng)設(shè)定為大氣壓,即p=0;自由出流完全發(fā)展的區(qū)域:

式中:n為出口截面外法線方向;k為紊動(dòng)能;ε為紊動(dòng)耗散率。

2.4.3固壁邊界條件

在固壁邊界上采用無(wú)滑移邊界條件,即并利用實(shí)測(cè)的數(shù)據(jù)確定固壁糙率。

2.4.4自由水面

采用VOF方法捕捉水氣交界面,得到自由水面。

對(duì)于邊界上的紊流參數(shù)k和ε,其初始值如果不能從試驗(yàn)中取得,只能采用一定的經(jīng)驗(yàn)公式來(lái)估算。由式(5)(6)來(lái)估計(jì)進(jìn)口的k和ε。式中:I為底坡;u0為入口處的速度;Cμ為躍升系數(shù),可取0.09;l為湍流長(zhǎng)度;L為關(guān)聯(lián)尺寸,對(duì)于充分發(fā)展的湍流,可取L等于水力半徑。

2.5 數(shù)值計(jì)算過(guò)程

本文模型的數(shù)值計(jì)算過(guò)程分2步,先是模擬水槽內(nèi)無(wú)植物壩覆蓋時(shí)水流定常流流動(dòng),直至收斂,然后在這個(gè)收斂流場(chǎng)基礎(chǔ)上模擬水槽內(nèi)有植物壩覆蓋時(shí)的水氣兩相非定常流流動(dòng)。否則如果直接模擬水氣兩相非定常流,迭代計(jì)算容易發(fā)散[5]。在此基礎(chǔ)上應(yīng)用VOF方法捕捉自由水面,變成水氣兩相非定常流模擬。數(shù)值模型的邊界條件設(shè)置:水流入口采用速度入口,水流入口水位根據(jù)試驗(yàn)時(shí)的實(shí)際工況設(shè)定;入口水位以上至下游出口斷面設(shè)置為氣體壓力進(jìn)口,靜壓設(shè)為零;下游出口分為空氣出口和液態(tài)水流出口兩部分,空氣出口部分采用壓力出口,靜壓為零;各邊界條件的紊流參數(shù)均按照經(jīng)驗(yàn)公式確定湍動(dòng)能和紊動(dòng)耗散速率,固體邊壁當(dāng)量粗糙度以試驗(yàn)取0.05m。壓力插值格式選用body force weighted格式,動(dòng)量、紊動(dòng)能、紊動(dòng)耗散速率的離散均采用一階迎風(fēng)格式[6]。時(shí)間步長(zhǎng)為0.001s,需求時(shí)間步設(shè)定為100000,每個(gè)時(shí)間步內(nèi)的最大迭代計(jì)算次數(shù)為20次。

圖4為模擬有植物壩覆蓋水槽的3個(gè)典型斷面上的流速值與試驗(yàn)值的對(duì)比。

圖4 模擬流速值與試驗(yàn)值的對(duì)比

水流通過(guò)植物壩前端時(shí),隨著植物壩阻水所產(chǎn)生的壅高水流流速急劇減小;水流流至植物壩中段時(shí),流速趨于穩(wěn)定,沿水深縱剖面上流速值變化不大;水流通過(guò)植物壩末端時(shí),流速受植物枝葉對(duì)水流反作用力的影響,隨著水深的增加流速值再次減小。

由圖4中3組典型斷面數(shù)值模型計(jì)算流速值與物理模型試驗(yàn)值的比較可以看出,3個(gè)斷面上的模擬流速沿水深分布趨勢(shì)與試驗(yàn)測(cè)得的結(jié)果基本一致。由此表明,本文利用Fluent軟件計(jì)算和VOF方法耦合的數(shù)學(xué)模型是正確的,用多孔介質(zhì)域模擬植物壩覆蓋區(qū)域是可行的。

3 結(jié) 語(yǔ)

a.水槽內(nèi)放置植物壩后,槽內(nèi)水流的流場(chǎng)與植樹前有了明顯的不同,其中上游壅水區(qū)是植物壩對(duì)上游水流阻擋作用的直接結(jié)果。

b.將數(shù)值計(jì)算得到的流速值與模型試驗(yàn)流速值實(shí)測(cè)結(jié)果進(jìn)行了比對(duì),兩者吻合較好。表明本文的數(shù)學(xué)模型是正確的,用多孔介質(zhì)域模擬植物壩覆蓋區(qū)域是可行的。

由數(shù)值模型計(jì)算流速值與物理模型試驗(yàn)值的比較表明,Fluent軟件計(jì)算和VOF方法耦合的數(shù)學(xué)模型是正確的,用多孔介質(zhì)域模擬植物壩覆蓋區(qū)域是可行的[7]。通過(guò)較系統(tǒng)的有植物壩河道的試驗(yàn)研究分析,得到了一些有價(jià)值的數(shù)據(jù)和結(jié)論,但這僅僅是初步的研究工作,還有待于深入的理論分析和采用更多的數(shù)學(xué)工具進(jìn)行更多的物理模型試驗(yàn)和數(shù)值模擬方面的研究。

[1]程艷,李森,邱秀云,等.河渠種樹水流特性試驗(yàn)研究[J].新疆農(nóng)業(yè)大學(xué)學(xué)報(bào),2003,26(2):59-64.

[2]閻潔,周著,邱秀云,等.植物“柔性壩”阻水試驗(yàn)及固沙機(jī)理分析[J].新疆農(nóng)業(yè)大學(xué)學(xué)報(bào),2004,27(3):40-45.

[3]袁夢(mèng),黃本勝,邱秀云,等.水葫蘆覆蓋區(qū)水流阻力效應(yīng)試驗(yàn)研究[J].廣東水利水電.2008(2):29-34.

[4]王福軍.計(jì)算流體動(dòng)力學(xué)分析:CFD軟件原理與應(yīng)用[M].北京:清華大學(xué)出版社,2004.

[5]陳界仁,朱曉丹.含柔性沉水植物明渠水流運(yùn)動(dòng)特性試驗(yàn)研究[J].水利水電科技進(jìn)展,2008,28(3):16-19.

[6]呂升齊,唐洪武,閆靜.有無(wú)植物條件下明渠水流紊動(dòng)特性對(duì)比[J].水利水電科技進(jìn)展,2007,27(6):57-60.

[7]吳航,拾兵,陳志娟.淹沒(méi)狀態(tài)下河道植物偏斜高度與糙率的關(guān)系[J].水利水電科技進(jìn)展,2008,28(1):20-22.