米博宇，張小峰，任 實

(1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072; 2.中國長江三峽集團公司，湖北 宜昌 443133)

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

米博宇1，張小峰1，任 實2

(1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072; 2.中國長江三峽集團公司，湖北 宜昌 443133)

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

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

1 研究背景

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

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

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

(1)

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

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

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

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

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

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

3 計算工況

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

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

4 模擬結(jié)果

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

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

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

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

表2 各工況特性

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

[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