米博宇，張小峰，任 實(shí)

(1.武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072; 2.中國(guó)長(zhǎng)江三峽集團(tuán)公司，湖北 宜昌 443133)

斜坡上的密度流數(shù)值模擬研究

米博宇1，張小峰1，任 實(shí)2

(1.武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072; 2.中國(guó)長(zhǎng)江三峽集團(tuán)公司，湖北 宜昌 443133)

為了探討均勻密度環(huán)境水體中斜坡密度流的運(yùn)動(dòng)規(guī)律，建立了立面二維RNGk-ε紊流數(shù)學(xué)模型，通過與已有試驗(yàn)資料對(duì)比，驗(yàn)證了該模型的合理性與準(zhǔn)確性。利用該模型模擬不同坡角和流量下的斜坡密度流，研究結(jié)果表明:密度流頭部流速與坡角具有一定的函數(shù)關(guān)系，存在一個(gè)最優(yōu)坡角使得相同條件下的頭部流速最大；頭部流速與浮力通量的三次方根之間并非嚴(yán)格的正比例函數(shù)關(guān)系，在流量較小或較大時(shí)將發(fā)生偏離；密度流運(yùn)動(dòng)過程中，頭部形態(tài)不斷擴(kuò)大，頭部的厚長(zhǎng)比逐漸減??；頭部擴(kuò)大的速率隨著坡角和流量的增大而增大，最后逐漸趨于一個(gè)穩(wěn)定速率。研究結(jié)果能夠幫助進(jìn)一步了解斜坡密度流的運(yùn)動(dòng)規(guī)律。

密度流；斜坡；頭部流速；數(shù)值模擬；RNGk-ε紊流模型

1 研究背景

本文將在前人研究成果的基礎(chǔ)之上，采用數(shù)值模擬的方法，與已有試驗(yàn)資料進(jìn)行對(duì)比分析，進(jìn)一步探討均勻密度環(huán)境水體中斜坡密度流的運(yùn)動(dòng)規(guī)律，重點(diǎn)分析斜坡坡角與進(jìn)口流量對(duì)密度流頭部流速與形態(tài)的影響。

2 數(shù)學(xué)模型與驗(yàn)證

Choi等[5]采用標(biāo)準(zhǔn)k-ε紊流模型模擬斜坡上的密度流，驗(yàn)證了其可行性，本文在此基礎(chǔ)上改用RNGk-ε紊流模型，根據(jù)Patankar[6]求解流場(chǎng)的思想，將方程組中各方程以統(tǒng)一形式表達(dá)為

(1)

控制方程包括Reynolds平均Navier-Stokes方程(簡(jiǎn)稱RANS)和RNGk-ε紊流模型方程(具體參考文獻(xiàn)[7])。計(jì)算方法采用Patankar提出的SIMPLER算法。壓力采用SIMPLER算法中的壓力方程直接求解，與傳統(tǒng)的采用Boussinesq假設(shè)將密度的變化轉(zhuǎn)換成浮力項(xiàng)這一做法相比，本文的數(shù)學(xué)模型直接將壓力差所導(dǎo)致的浮力還原到壓力場(chǎng)中，因?yàn)楂@取壓力場(chǎng)時(shí)并沒有采用常密度假設(shè)，而是根據(jù)實(shí)時(shí)的密度場(chǎng)與速度場(chǎng)求解得來。

邊界條件參考陶文銓[7]在其著作中給出的邊界提法，其中壓力與壓力修正方程在各類邊界上采用二類邊界條件，動(dòng)量方程與k-ε方程在邊壁上采用壁面函數(shù)法，k-ε模型方程中的參數(shù)取值與文獻(xiàn)[7]中的推薦取值一致。

采用典型算例對(duì)數(shù)學(xué)模型進(jìn)行驗(yàn)證，算例取自雙層密度流體中的侵入密度流試驗(yàn)，這類問題已經(jīng)被Britter等[8]、Lowe等[9]多位學(xué)者廣泛研究過，并取得了基本一致的結(jié)論。算例設(shè)置如圖1(a)所示，豎直平面內(nèi)有一長(zhǎng)(x方向)1.0 m、高度(y方向)>0.2 m的矩形水槽，3種不同密度的水體分別位于A，B，C 3個(gè)區(qū)域，其中ρA=(ρB+ρC)/2，ρB<ρC，A區(qū)域與B，C區(qū)域之間有一塊擋板隔離，初始處于靜止?fàn)顟B(tài)。水槽頂部為空氣邊界，四周及底部為固壁邊界。突然撤去中間的擋板，由于ρB<ρA<ρC，重力作用下，A與B，C的交界面上會(huì)產(chǎn)生壓強(qiáng)差，驅(qū)使原來靜止的流體開始流動(dòng)，從而形成侵入密度流，如圖1(b)所示。

圖1 驗(yàn)證算例計(jì)算區(qū)域與侵入密度流形態(tài)Fig.1 Calculation areas and intrusive density current shape of verification example

表1 本文數(shù)值模擬結(jié)果與Lowe等[9]試驗(yàn)結(jié)果對(duì)比

圖2 斜坡計(jì)算區(qū)域示意圖Fig.2 Schematic diagram of calculation area of slope

3 計(jì)算工況

對(duì)于斜坡上的密度流，本文的計(jì)算區(qū)域如圖2所示，斜坡頂端有一入口，高密度水體從此處流入環(huán)境水體中，斜坡角度為θ，豎直高度為H2，斜坡頂端距水面高度為H1，整個(gè)計(jì)算區(qū)長(zhǎng)度根據(jù)θ的不同而略有變化。所有計(jì)算工況中，進(jìn)口高度均為1 cm，為了盡量減小水面對(duì)密度流的影響，取H1=15 cm，H2=15 cm，坡角θ范圍為5°～40°，環(huán)境水體密度ρa(bǔ)=1 000 kg/m3，入流水體密度ρi=1 015 kg/m3。

為了分析坡角和流量對(duì)密度流的影響，本文的計(jì)算工況分為2大類(工況編號(hào)以“I-”和“II-”標(biāo)志)。第1類工況共13個(gè)，進(jìn)口流量(單寬)統(tǒng)一為3×10-4m3/(m·s)，坡角標(biāo)在工況編號(hào)中，如“I-5”表示第1類中坡角為5°的工況，以此類推；第2類工況共11個(gè)，坡角統(tǒng)一為10°，進(jìn)口流量標(biāo)在工況編號(hào)中，如“II-3”表示第2類中進(jìn)口流量為3×10-4m3/(m·s)的工況，以此類推。

4 模擬結(jié)果

斜坡上每個(gè)工況中的密度流都有一個(gè)明顯的頭部，后面緊跟較薄的尾部，部分坡角下的密度流在同一時(shí)刻的形態(tài)如圖3所示。

圖3 密度流形態(tài)Fig.3 Shape of density current

從圖3中可以看出，頭部最前端都有一定的抬升，都是離開邊壁的，坡角越大，在同一時(shí)刻形成的密度流頭部也越大，尾部形態(tài)則沒有明顯區(qū)別。這與Britter等[1]的試驗(yàn)現(xiàn)象相吻合。

取各個(gè)工況中密度流頭部的位置和形態(tài)與時(shí)間的對(duì)應(yīng)關(guān)系進(jìn)行分析，可以發(fā)現(xiàn)密度流在沿斜坡向下運(yùn)動(dòng)的過程中頭部的速度Uf經(jīng)過一定的調(diào)整后保持恒定，頭部大小(包括長(zhǎng)度L和厚度h)則因卷吸作用不斷增大，且與時(shí)間近似呈線性關(guān)系。具體數(shù)值列于表2中，其中Uf是恒定以后的頭部速度，αL是L隨時(shí)間的增長(zhǎng)率，αh是h隨時(shí)間的增長(zhǎng)率，頭部雷諾數(shù)Re=(g′0q)1/3h/ν，g′0=2g(ρi-ρa(bǔ))/(ρi+ρa(bǔ))為入流的相對(duì)重力加速度。由于頭部厚度h是不斷變化的，本文選取分析時(shí)段內(nèi)的最小Re作為代表值，即表2中的Remin。帶星號(hào)(*)的2個(gè)工況中，L和h與時(shí)間t的關(guān)系更接近于一條曲線，而非其他工況中的線性關(guān)系，說明增長(zhǎng)率αL和αh未能穩(wěn)定，因此其相應(yīng)結(jié)果不參與本文數(shù)據(jù)分析，但頭部位置與時(shí)間的關(guān)系明確，且與其他工況規(guī)律一致，故將其Uf值納入分析范圍。

表2 各工況特性

4.1 頭部流速的規(guī)律

圖4 頭部流速與坡角的關(guān)系Fig.4 Relations of head velocity vs. slope gradient and

密度流自進(jìn)口流入環(huán)境水體以后，初始流速發(fā)展到最后的恒定流速，會(huì)經(jīng)歷一個(gè)調(diào)整階段。本文第1類工況的進(jìn)口流速均為0.03 m/s，從圖4(a)可以看出，所有工況中的密度流頭部都會(huì)經(jīng)歷一個(gè)加速過程到它的恒定流速，這一過程歷時(shí)為5～12 s不等，且歷時(shí)長(zhǎng)短與圖4(a)曲線走勢(shì)大致吻合。同樣在第2類工況中也有類似情況。調(diào)整完成以后進(jìn)入穩(wěn)定階段，頭部流速不再變化。

影響密度流頭部流速的主要力學(xué)因素有密度流在環(huán)境水體中的凈重力G、斜坡上的阻力fs與支持力N、以及環(huán)境水體的阻力fa，這些力之間的平衡決定了頭部流速的恒定。在沿斜坡方向，G的分力是密度流流動(dòng)的驅(qū)動(dòng)力，流動(dòng)過程中形成的阻力fs和fa是與流速呈正相關(guān)的。當(dāng)初始流速過小時(shí)，驅(qū)動(dòng)力大于阻力，密度流會(huì)加速流動(dòng)，反之則會(huì)減速，直到某一流速恰好使得阻力等于驅(qū)動(dòng)力以后，密度流處于受力平衡狀態(tài)，流速不再變化。同樣條件下，坡角越大，驅(qū)動(dòng)力也越大，達(dá)到平衡時(shí)需要的流速也越大。但fa不僅與流速有關(guān)，還與密度流頭部的形態(tài)有關(guān)，頭部越大，受到環(huán)境水體的阻力自然越大。從圖3可以看出，同一時(shí)刻，坡角越大，密度流頭部也越大，由此帶來的阻力減小了對(duì)于流速所引起的阻力的需求，因而密度流頭部流速并不是隨著坡角增大而一直增大的。Britter等[1]的試驗(yàn)數(shù)據(jù)與本文的數(shù)值結(jié)果都顯示在坡角為20°左右時(shí)，密度流流速最大，至于這一最優(yōu)坡角的具體值以及跟其他因素間的關(guān)系還有待進(jìn)一步研究。

圖5 Uf與(工況II-3至II-11)，q的關(guān)系Fig.5 Relations of Uf vs.(in conditions from II-3 to II-11) and q

密度流頭部流速與驅(qū)動(dòng)力大小密切相關(guān)，當(dāng)坡角一定時(shí)，驅(qū)動(dòng)力大小取決于頭部的相對(duì)重力加速度g′0h，而g′0h由頭部與環(huán)境水體密度決定。取同一時(shí)刻各個(gè)流量下的密度流頭部中心位置的密度分布(見圖6)進(jìn)行分析，分布曲線從內(nèi)到外流量依次增加。

圖6 同一時(shí)刻各個(gè)流量下頭部密度分布Fig.6 Distribution of head density under various flow rates at the same time

4.2 頭部形態(tài)的規(guī)律

密度流頭部在沿斜坡向下運(yùn)動(dòng)的過程中，由于尾部高密度水體的匯入，以及與環(huán)境水體間的摻混作用，頭部是不斷增大的[1]。以頭部厚度h與長(zhǎng)度L作為參考指標(biāo)，本文數(shù)值模擬的結(jié)果顯示h和L與時(shí)間t有著較好的線性關(guān)系。從表2可以看出，h和L隨時(shí)間的增長(zhǎng)率有明顯的差別，L增長(zhǎng)的速率要快于h，說明密度流頭部在運(yùn)動(dòng)過程中是不斷扁平化發(fā)展的。各工況的計(jì)算結(jié)果表明，隨著坡角的增大，同一時(shí)刻形成的密度流頭部的厚長(zhǎng)比h/L是逐漸增大的，在θ=40°工況中最大能到0.72，圖3也直觀地反映了這一現(xiàn)象。隨著密度流頭部的運(yùn)動(dòng)， h/L是逐漸減小，在θ=5°工況中最小能到0.44。這些現(xiàn)象均與試驗(yàn)觀測(cè)相符。分析h和L的增長(zhǎng)率與坡角和流量的關(guān)系，如圖7所示。

圖7 h和L的增長(zhǎng)率與坡角、流量的關(guān)系Fig.7 Relations of growth rates of h and L vs. slope gradient and discharge

從圖7中可以看出，隨著坡角或流量的增大，h和L增長(zhǎng)率逐漸趨于一個(gè)穩(wěn)定水平，并非線性增長(zhǎng)。

5 結(jié) 論

(1) RNGk-ε紊流模型可以準(zhǔn)確地模擬密度流的流動(dòng)，對(duì)于頭部流速等參數(shù)的模擬結(jié)果與試驗(yàn)觀測(cè)一致，因此可以用來進(jìn)行密度流的相關(guān)研究。

(3) 坡角對(duì)密度流頭部流速的影響較小，數(shù)值模擬的結(jié)果與試驗(yàn)數(shù)據(jù)均表明存在一個(gè)最優(yōu)坡角使得相同條件下密度流的頭部流速最大。造成這一現(xiàn)象的原因是坡角的變化在改變密度流驅(qū)動(dòng)力的同時(shí)也改變了頭部的形態(tài)，使得阻力也發(fā)生了變化，二者的相互制約形成了頭部流速與坡角間的關(guān)系。

(4) 密度流在運(yùn)動(dòng)過程中頭部形態(tài)不斷擴(kuò)大，坡角越大，流量越大，擴(kuò)大的速率越快，但這一速率并不是隨坡角或流量線性增長(zhǎng)的，而是逐漸趨于一個(gè)穩(wěn)定值。頭部在擴(kuò)大的同時(shí)，形狀也在發(fā)生改變，隨著時(shí)間推移，頭部越來越扁平。

[1]BRITTERRE,LINDENPF.TheMotionoftheFrontofaGravityCurrentTravellingdownanIncline[J].JournalofFluidMechanics, 1980, 99(3): 531-543.

[2]ALAVIANV.BehaviorofDensityCurrentsonanIncline[J].JournalofHydraulicEngineering, 1986, 112(1): 27-42.

[3]BAINESPG.MixinginFlowsdownGentleSlopesintoStratifiedEnvironments[J].JournalofFluidMechanics, 2001, 443: 237-270.

[4]CENEDESEC,ADDUCEC.MixinginaDensity-drivenCurrentFlowingdownaSlopeinaRotatingFluid[J].JournalofFluidMechanics, 2008, 604: 369-388.

[5]CHOISU,GARCíAMH. k-εTurbulenceModelingofDensityCurrentsDevelopingTwoDimensionallyonaSlope[J].JournalofHydraulicEngineering, 2002, 128(1): 55-63.

[6]PATANKARSV.NumericalHeatTransferandFluidFlow[M].USA:HemispherePublishingCorporation, 1980: 1-154.

[7] 陶文銓. 數(shù)值傳熱學(xué)(第2版) [M]. 西安: 西安交通大學(xué)出版社, 2001: 347-360.

[8]BRITTERRE,SIMPSONJE.ANoteontheStructureoftheHeadofanIntrusiveGravityCurrent[J].JournalofFluidMechanics, 1981, 112: 459-466.

[9]LOWERJ,LINDENPF,ROTTMANJW.ALaboratoryStudyoftheVelocityStructureinanIntrusiveGravityCurrent[J].JournalofFluidMechanics, 2002, 456: 33-48.

[10]BENJAMINTB.GravityCurrentsandRelatedPhenomena[J].JournalofFluidMechanics, 1968, 31(2): 209-248.

(編輯：羅 娟)

Numerical Simulation of Density Current on a Slope

MI Bo-yu1, ZHANG Xiao-feng1, REN Shi2

(1.State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China；2.China Three Gorges Corporation, Yichang 443133，China)

A vertical two-dimensional RNGk-εturbulent model is established and its reasonability and accuracy are verified by comparison with existing experimental data. Density current on a slope in the presence of different slope gradients and discharges is simulated, and results reveal that 1) the head velocity of density current has a function relationship with the slope gradients, and there is an optimal slope gradients which maximizes the head velocity under the same condition; 2)the relation between head velocity and the cubic root of buoyance flux is not a strict proportional function, deviating under small or large discharge; 3)the head shape is enlarged in the motion process of density current, and the ratio of thickness to length of the head decreases gradually; 4)the growth rate of the head increases with the increase of slope gradients and discharges, and finally tends to a steady rate. These results could help further understand the motion pattern of density current on slope.

density current; slope; head velocity; numerical simulation; RNGk-εturbulent model

2016-04-12；

2016-06-03

米博宇(1992-)，男，湖北安陸人，碩士研究生，主要從事水力學(xué)數(shù)值模擬方面的研究，(電話)15527820486(電子信箱)miboyu@whu.edu.cn。

10.11988/ckyyb.20160346

2017,34(7):60-64

TV133.2

A

1001-5485(2017)07-0060-05