修正Realizable k-ε模型在高壓淹沒水射流中的應(yīng)用

余健翔， 劉 劍， 王觀石*， 羅嗣海， 汪 杰

(1.江西理工大學(xué)土木與測繪工程學(xué)院， 贛州 341000； 2.江西理工大學(xué)資源與環(huán)境工程學(xué)院， 贛州 341000；3.攀枝花學(xué)院土木與建筑工程學(xué)院， 攀枝花 617000)

高壓水射流是通過增壓泵等設(shè)備和特定形狀的噴嘴產(chǎn)生的具有極高的能級(jí)密度的高速射流束。高壓水射流按工作介質(zhì)與環(huán)境介質(zhì)分類，可分為淹沒水射流和非淹沒水射流，淹沒水射流在石油開采、煤炭開采、航道疏浚、地基處理等方面應(yīng)用廣泛且具有安全、降低設(shè)備損耗、適應(yīng)性強(qiáng)等優(yōu)點(diǎn)[1]。離子型稀土礦原地浸礦的高效收液系統(tǒng)中有大量的水平填砂導(dǎo)流孔[2]，需要通過高壓淹沒水射流的沖蝕能力，在風(fēng)化層礦土中快速鉆進(jìn)導(dǎo)流孔；水射流的滯止壓力直接影響射流的沖蝕能力；研究水射流滯止壓力沿程變化是確定導(dǎo)流孔施工技術(shù)參數(shù)的關(guān)鍵，進(jìn)而也影響到該工藝的生產(chǎn)成本和生產(chǎn)周期。

目前對(duì)高壓淹沒水射流的研究方法主要有實(shí)驗(yàn)和數(shù)值模擬。楊樹人等[3]對(duì)空化射流噴嘴的內(nèi)部流場進(jìn)行了數(shù)值模擬研究；劉海霞等[4]、Nobel等[5]通過實(shí)驗(yàn)研究淹沒射流的動(dòng)壓沖擊作用；張偉等[6]、張欣瑋等[7]采用數(shù)值模擬與粒子圖像測速(particle image velocimetry, PIV)實(shí)驗(yàn)相結(jié)合的方法研究流場特性。從現(xiàn)有的研究來看，高壓淹沒水射流研究的困難性主要在于兩方面：一方面是缺少準(zhǔn)確有效的實(shí)驗(yàn)方法對(duì)流場各部分的速度與壓力等數(shù)據(jù)進(jìn)行測量，現(xiàn)有的實(shí)驗(yàn)研究大多是對(duì)中低壓力的淹沒水射流進(jìn)行研究，對(duì)于高壓的淹沒水射流實(shí)驗(yàn)研究較少且大多是基于實(shí)驗(yàn)現(xiàn)象給出定性分析，少有給出定量的表達(dá)；另一方面在于對(duì)空化機(jī)理認(rèn)識(shí)不充分，含空化氣泡的射流速度計(jì)算沒有普遍適用的方法。

對(duì)于淹沒水射流這種復(fù)雜情況的流動(dòng)來說，數(shù)值模擬是一個(gè)十分有效的方法,但其不能完全替代實(shí)驗(yàn)，數(shù)值計(jì)算的結(jié)果一般都需要驗(yàn)證，并用規(guī)范的實(shí)驗(yàn)來確認(rèn)[8]。

綜上，現(xiàn)采用實(shí)驗(yàn)和數(shù)值模擬相結(jié)合的方法對(duì)淹沒水射流進(jìn)行研究。實(shí)驗(yàn)測定現(xiàn)有實(shí)驗(yàn)較少涉及的0 ～ 40 MPa射流壓力范圍下的沿程滯止壓力。根據(jù)實(shí)驗(yàn)工況建立模型對(duì)淹沒水射流進(jìn)行數(shù)值模擬，考慮到現(xiàn)有的空化模型不完善[9]，采用改良模型常數(shù)的方法，通過與實(shí)驗(yàn)結(jié)果對(duì)比驗(yàn)證得出各射流壓力對(duì)應(yīng)的模型常數(shù)，提高模型預(yù)測淹沒射流的準(zhǔn)確性。該方法對(duì)比使用空化模型來說，更簡便，占用的計(jì)算機(jī)資源較少，對(duì)于工程應(yīng)用具有較大的實(shí)際意義。

1 淹沒水射流實(shí)驗(yàn)

1.1 淹沒水射流結(jié)構(gòu)

一般把淹沒水射流分為初始段和主體段。在噴嘴出口處，射流與環(huán)境水發(fā)生劇烈的動(dòng)量交換和湍動(dòng)擴(kuò)散，其中心線附近部分射流介質(zhì)仍會(huì)保持噴嘴出口時(shí)的速度，該部分介質(zhì)稱為等速核心也稱為勢流核，射流初始段由勢流核與介質(zhì)交換的混合區(qū)組成。該區(qū)域射流速度大，其與環(huán)境水會(huì)造成強(qiáng)烈的剪切作用產(chǎn)生空化現(xiàn)象。勢流核消散后的區(qū)域稱為射流主體段，此時(shí)射流的湍動(dòng)性增強(qiáng)，進(jìn)一步卷吸環(huán)境介質(zhì)，射流軸心速度迅速降低，在射流核心區(qū)產(chǎn)生的空化氣泡會(huì)在這個(gè)區(qū)域內(nèi)發(fā)展增大，并隨著射流運(yùn)動(dòng)到一定距離后湮滅?？栈F(xiàn)象涉及復(fù)雜多變的物理現(xiàn)象，如氣體與液體之間的兩相轉(zhuǎn)變，其對(duì)射流速度等特性影響不可忽略。

1.2 水射流實(shí)驗(yàn)

近年來對(duì)淹沒水射流實(shí)驗(yàn)研究主要采用粒子圖像測速法對(duì)射流流場進(jìn)行測量并采用高速攝影技術(shù)對(duì)淹沒射流中空化現(xiàn)象進(jìn)行觀測[10-13]。上述方法主要適用于較低射流壓力及流速的情況，對(duì)于高射流壓力及流速的情況，PIV采用的示蹤粒子對(duì)介質(zhì)的跟隨性差且分布不均勻；高濃度的空化氣泡在圖像上重疊，難以辨認(rèn)射流內(nèi)部結(jié)構(gòu)。對(duì)于本文實(shí)驗(yàn)所采用的高壓力工況，采用經(jīng)典的皮托管法測量沿程滯止壓力，實(shí)驗(yàn)裝置如圖1所示。

圖1 水射流實(shí)驗(yàn)設(shè)備Fig.1 The equipment of water jet experiment

實(shí)驗(yàn)裝置主要包括水箱、鋼板箱、高壓泵、高壓閥、高壓噴嘴、滑臺(tái)、高壓管、壓力變送器、高速采集儀。鋼板箱使用8塊4 mm厚鋼板組成，尺寸為60 cm× 40 cm× 80 cm，底部裝有固定壓力變送器的支架，壓力變送器和皮托管通過螺栓剛性連接固定在支架上，皮托管管直徑為8 mm，測壓孔的內(nèi)徑為1.0 mm，外徑為1.5 mm，測壓孔距容器底面20 cm，底部壁面對(duì)射流影響較小故射流可以作為自由射流。高壓泵為柱塞泵可輸出的最大壓力為50 MPa，最大流量為23 L/min。高壓閥為節(jié)流調(diào)壓閥?；_(tái)固定在鋼板箱側(cè)邊，正面有兩個(gè)滑塊固定高壓管，滑臺(tái)可以上下移動(dòng)，并在任意位置鎖定。高壓管為長約1.2 m的無縫鋼管，內(nèi)徑 7 mm，固定在支架上，下部連接可更換高壓噴嘴，上部接一個(gè)三通管，三通管一端接水管，另一端接一個(gè)量程為50 MPa的壓力變送器。實(shí)驗(yàn)采用的壓力變送器為裝有阻尼器的擴(kuò)散硅式壓力變送器，輸出信號(hào)為4～20 mA電信號(hào)，測得滯止壓力在高速采集儀上實(shí)時(shí)顯示，采集頻率為10 Hz。

對(duì)于噴嘴出口處的射流壓力，由于射流速度的波動(dòng)和空化爆鳴劇烈的影響，近距離測試易損傷傳感器芯體，故本實(shí)驗(yàn)射程為6～16 cm處。實(shí)驗(yàn)測定了射流壓力5 ～ 40 MPa時(shí)的沿程軸心滯止壓力，其結(jié)果如圖2所示，模坐標(biāo)為距噴嘴軸心距離x與噴嘴直徑d的比值(x/d)。

圖2 軸心滯止壓力的實(shí)驗(yàn)值Fig.2 Experimental results of axial stagnationpressure

2 數(shù)值模擬

2.1 模型的控制方程

湍流模型選擇Realizablek-ε模型來對(duì)流場進(jìn)行模擬計(jì)算。與標(biāo)準(zhǔn)k-ε模型相比，Realizablek-ε采用了新的湍流黏度公式，公式中的系數(shù)Cμ不再為一個(gè)常量，而是與湍流層流應(yīng)變和旋度等有關(guān)的函數(shù)，而且為湍動(dòng)能耗散率方程增加了新的傳輸公式。這使Realizablek-ε能更精確地預(yù)測平板和圓柱射流的發(fā)散性，而且對(duì)于預(yù)測旋轉(zhuǎn)流動(dòng)、強(qiáng)逆壓梯度的邊界層流動(dòng)、流動(dòng)分離和二次流也有很好的表現(xiàn)[14]。

Realizablek-ε模型的控制方程為

Gk+Gb-ρε-YM+Sk

(1)

(2)

式中：ρ為液體密度，kg/m3；μ為液體黏度，Pa·s；μt為液體渦黏系數(shù)，Pa·s；k為湍流動(dòng)能，m2/s2;ε為湍流耗散率，m2/s3;Gk為由于平均速度梯度而引起的湍動(dòng)能k的產(chǎn)生項(xiàng)；Gb為由于浮力而引起的湍流動(dòng)能k的產(chǎn)生項(xiàng)；YM代表可壓湍流中脈動(dòng)擴(kuò)張的貢獻(xiàn)；C2、C1ε為常量；σk、σε為k方程和ε方程的Prandtl數(shù)；Sk和Sε為用戶自定義的源項(xiàng)。由于把淹沒水射流考慮為單一介質(zhì)的不可壓縮流動(dòng)，且恒定射流不考慮用戶定義的源項(xiàng)，有Gb=0，YM=0，Sk=0，Sε=0。Realizablek-ε模型的控制方程變?yōu)?/p>

(3)

(4)

式(4)中：系數(shù)C1與時(shí)均應(yīng)變率Sij聯(lián)系在一起，由式(5)求得，即

(5)

式(5)中：

(6)

Gk為湍動(dòng)能的產(chǎn)生項(xiàng)，表示為

(7)

對(duì)式(7)用Boussinesq假設(shè)處理，Gk的表達(dá)式可改為

Gk=μtS2

(8)

(9)

在標(biāo)準(zhǔn)k-ε模型中系數(shù)Cμ為常數(shù)，而在Realizablek-ε模型Cμ按式(10)確定，即

(10)

式(10)中：

(11)

(12)

(13)

(14)

(15)

模型常數(shù)：C1ε=1.44,C2=1.9,σk=1.0,σε=1.2。

2.2 模型與邊界條件

對(duì)流場進(jìn)行三維建模，建立一個(gè)半徑為5 cm，高為25 cm的圓柱形計(jì)算域，噴嘴位于其中央頂部。噴嘴結(jié)構(gòu)為錐直型噴嘴，結(jié)構(gòu)如圖3所示，其中進(jìn)口直徑D=7 mm，出口直徑d=1.5 mm，L1=20 mm,收縮段長度L2=10 mm,出口段長度L3=5 mm,收縮角度θ=30°。

圖3 噴嘴結(jié)構(gòu)Fig.3 Structure of nozzle

進(jìn)口和出口邊界條件設(shè)置為壓力入口和壓力出口，噴嘴末端7 mm管道為進(jìn)口邊界,射流壓力值與實(shí)驗(yàn)一致。圓柱頂部除噴嘴進(jìn)口外為出口邊界，設(shè)置靜壓值為0，其他邊界全部設(shè)置為壁面，壁面采用無滑移條件，近壁面采用標(biāo)準(zhǔn)壁面函數(shù)處理，y+值取為11～30。網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格，網(wǎng)格總數(shù)約為150萬。網(wǎng)格劃分如圖4所示，沿射流方向采用等間距結(jié)構(gòu)化網(wǎng)格，節(jié)點(diǎn)數(shù)共有303個(gè)；對(duì)于垂直射流方向的網(wǎng)格采用變間距網(wǎng)格，靠近噴嘴壁面的網(wǎng)格間距越小，網(wǎng)格間距增長率為1.2，共有91個(gè)節(jié)點(diǎn)。在噴嘴處對(duì)網(wǎng)格進(jìn)行加密，如圖5所示，主要對(duì)沿射流方向網(wǎng)格進(jìn)行加密，噴嘴出口處圓柱及收縮段的網(wǎng)格間距約為0.3 mm，而噴嘴進(jìn)口處的圓柱段網(wǎng)格間距約為0.8 mm。為模擬淹沒環(huán)境，流場為單一介質(zhì)水，其密度為998.2 kg/m3，動(dòng)力黏度為0.001 Pa·s。壓力-速度耦合算法采用SIMPLE算法，為滿足計(jì)算精度，各控制方程離散格式均采用二階迎風(fēng)格式，其他設(shè)置保持默認(rèn)。

圖4 網(wǎng)格劃分Fig.4 Gird partition

圖5 噴嘴網(wǎng)格加密Fig.5 nozzle mesh refinement

3 模擬結(jié)果及其對(duì)比

3.1 仿真結(jié)果及其分析

采用上述設(shè)置，在Fluent中對(duì)各射流壓力的淹沒水射流進(jìn)行數(shù)值模擬，以5、20、40 MPa計(jì)算結(jié)果為例。圖6顯示了無因次動(dòng)壓與無因次距離的關(guān)系，橫坐標(biāo)為距噴嘴軸心距離x與噴嘴直徑d的比值，縱坐標(biāo)為壓力值P與射流壓力Pjet的比值。數(shù)值模擬的結(jié)果具有相似性，即無因次動(dòng)壓只是與無因次軸心距離的函數(shù)，而與射流壓力和環(huán)境壓力無關(guān)。而實(shí)驗(yàn)值不具有相似性，隨著壓力值的增大，數(shù)值模擬與實(shí)驗(yàn)值之間的差值越大，這是由于淹沒水射流在射流過程中產(chǎn)生空化而造成的。文獻(xiàn)[5]指出滯止壓力的增大是由于空化氣泡減少了射流與周圍靜水之間的動(dòng)量交換，產(chǎn)生了減阻效果。只考慮單一介質(zhì)的淹沒射流數(shù)值模擬無法預(yù)測這種減阻效果，故在噴嘴附近的空化區(qū)內(nèi)模擬結(jié)果與實(shí)驗(yàn)值差較大，離噴嘴距離較遠(yuǎn)的位置比較接近實(shí)驗(yàn)值。

圖6 無因次動(dòng)壓與射程關(guān)系Fig.6 Relationship between dimensionless dynamic pressure and range

3.2 對(duì)模型常數(shù)的改良

由于現(xiàn)有的空化模型和多相流模型對(duì)于模擬空化射流仍不完善。將含空化氣泡的水射流簡化為混合流，以單相流來對(duì)淹沒水射流進(jìn)行數(shù)值模擬?？紤]k-ε模型常數(shù)對(duì)數(shù)值模擬結(jié)果的影響，該模型常數(shù)主要根據(jù)一些特定條件下的實(shí)驗(yàn)確定，其適用性有一定的局限性[15]，對(duì)模擬高壓淹沒水射流的情況需要對(duì)其默認(rèn)常數(shù)進(jìn)行修改。Realizablek-ε模型的模型常數(shù)有C2、σk、σε。為降低模型常數(shù)研究的復(fù)雜性，只對(duì)k方程和ε方程的湍流Prandtl數(shù)σk和σε進(jìn)行修改。

3.3 σk和σε對(duì)計(jì)算結(jié)果的影響

采用幾組不同的模型常數(shù)σk和σε進(jìn)行計(jì)算，根據(jù)其計(jì)算所得的結(jié)果值來研究不同模型常數(shù)對(duì)模型預(yù)測性能的影響。σk和σε的取值如表1所示。

表1 σk和σε取值

以射流壓力20 MPa為例，對(duì)照組A的各軸心速度分布如圖7所示，可以看出，σk主要影響射流主體段，隨著σk減少，射流初始段勢流核長度減小，射流主體段軸心速度增大，且σk從1變?yōu)?.8時(shí)其速度增長率較小，從0.8變至0.3時(shí)，速度增長率較大。而σk增大時(shí)，射流初始段勢流核長度增加，而射流主體段軸心速度幾乎不變，與默認(rèn)常數(shù)的計(jì)算結(jié)果一致。圖8表示了對(duì)照組B的軸心速度分布，增大σε對(duì)射流的軸心速度影響很小，其計(jì)算結(jié)果基本不變。而減小σε會(huì)減小射流在軸心距離10～75 mm的軸心速度，且變化幅度較小。通過比較和分析，可見σk對(duì)射流軸心速度的影響遠(yuǎn)大于σε。

(16)

式(16)中：m=1.658 MPa；n=0.171。

從圖9可以看出，射流壓力值5～15 MPa時(shí)，σk的變化較大，15～40 MPa時(shí)σk的變化較小。這是因?yàn)閺牡蛪毫Φ礁邏毫r(shí)，淹沒射流的空化云體積的增長速率隨著壓力增大是先增大后減小的。這與毛寧等[16]用高速攝影技術(shù)觀察的空化現(xiàn)象一致。

圖7 對(duì)照組A軸心速度沿程分布Fig.7 Velocity distribution of control group A along the axial way

圖8 對(duì)照組B軸心速度沿程分布Fig.8 Velocity distribution of control group B along the way

圖9 湍流Prandtl數(shù)和射流壓力的關(guān)系Fig.9 Relationship between turbulent Prandtl numbers σk and jet pressure

3.4 模擬結(jié)果

改良參數(shù)后對(duì)5～40 MPa淹沒射流進(jìn)行模擬，以20 MPa計(jì)算結(jié)果為例如圖10所示，改良常數(shù)后其模型預(yù)測準(zhǔn)確性在(40～70)d(6～11 cm)處大大提高，與實(shí)驗(yàn)值的誤差在10%左右，而在70d之后的非空化區(qū)之后，由于空化氣泡的消失，導(dǎo)致其滯止壓力迅速減小。只考慮單向流體的數(shù)值模擬無法預(yù)測這種變化導(dǎo)致其在12 cm后位置預(yù)測結(jié)果誤差較大。對(duì)于其他壓力依然有這個(gè)現(xiàn)象，5～20 MPa射流壓力的計(jì)算結(jié)果如圖11所示，為了便于對(duì)比，其局部放大圖如圖12所示。

圖10 20 MPa無因次動(dòng)壓沿程分布Fig.10 Distribution of dimensionless dynamic pressure at a jet pressure of 20 MPa

圖11 5～20 MPa無因次動(dòng)壓沿程分布Fig.11 Distribution of dimensionless dynamic pressure along the axial way at the jet pressure of 5～20 MPa

圖12 5～20 MPa局部無因次動(dòng)壓沿程分布Fig.12 Local distribution of dimensionless dynamic pressure along the axial way at the jet pressure of 5～20 MPa

圖13 沿程動(dòng)壓分布對(duì)比Fig.13 Dimensionless dynamic pressure distribution of group

對(duì)于近噴嘴區(qū)域，參考其他學(xué)者的實(shí)驗(yàn)值。Nobel等[5]測定了淹沒條件3 mm直徑錐形噴嘴在壓力為5、10、14.5、19.5 MPa的情況下沿程軸心動(dòng)壓分布。Shen等[17]用皮托管測定了錐形噴嘴2.8 mm不同射流速度下的沿程軸心動(dòng)壓分布。將本文數(shù)值模擬結(jié)果與Nobel等[5]和Shen等[17]測量的實(shí)驗(yàn)值進(jìn)行對(duì)比如圖13所示。本文結(jié)果與Shen等[17]的測量值吻合較好，對(duì)于Nobel等[5]的測量值，在射流壓力為5 MPa時(shí)，數(shù)值模擬與其實(shí)驗(yàn)值較接近，而隨著壓力增加其實(shí)驗(yàn)值偏小。這是由于不同實(shí)驗(yàn)的研究工況不同，Nobel等[5]實(shí)驗(yàn)壓力值是在圍壓為0.13 MPa的情況下測定的，導(dǎo)致其空化程度減弱，在空化程度較明顯的高壓力下與數(shù)值模擬結(jié)果差異較大。而本文與Shen等[17]實(shí)驗(yàn)的工況較類似，故數(shù)值模擬結(jié)果與其實(shí)驗(yàn)值吻合較好。

對(duì)于25～40 MPa壓力下的淹沒水射流，為了便于對(duì)比只給出(40～106)d射程范圍內(nèi)的無因次軸心動(dòng)壓分布圖與局部軸心動(dòng)壓分布圖，分別如圖14、圖15所示。25～40 MPa下數(shù)值模擬的預(yù)測準(zhǔn)確性與20 MPa時(shí)一致，在(40～70)d(6～11 cm)處預(yù)測結(jié)果較好，而在11 cm后的位置預(yù)測值略高于實(shí)測值，其原因也與20 MPa相同，都是由空化氣泡的減少所導(dǎo)致的。因?yàn)檩^遠(yuǎn)的射程滯止壓力低于500 kPa，破壞力很弱，對(duì)射流工作能力的影響不大。故對(duì)于較遠(yuǎn)距離的射流壓力來說，預(yù)測是否準(zhǔn)確對(duì)實(shí)際工程的應(yīng)用影響不大。

圖14 25～40 MPa局部動(dòng)壓沿程分布 Fig.14 Local distribution of dynamic pressure along the axial way at the jet pressure of 25～40 MPa

圖15 25～40 MPa局部無因次軸心動(dòng)壓分布Fig.15 Local distribution of dimensionless dynamic pressure along the axial way at the jet pressure of 25～40 MPa

通過比較0～40 MPa射流壓力下的軸心動(dòng)壓模擬值與實(shí)驗(yàn)值，可以發(fā)現(xiàn)修正模型常數(shù)可以有效提高數(shù)值模擬在(0～70)d射程范圍的準(zhǔn)確性，使誤差降低至10%以內(nèi)。且對(duì)比使用多相流和空化模型進(jìn)行模擬，該方法占用計(jì)算機(jī)資源更少，較為簡便。

4 結(jié)論

(1)由于空化作用，高壓淹沒水射流的無因次速度/滯止壓力不再具有自相似性，模型常數(shù)不變的Realizablek-ε模型無法準(zhǔn)確預(yù)測射流的沿程滯止壓力。

(2)Realizablek-ε模型模擬高壓淹沒水射流，修正模型常數(shù)可以顯著提高(0～70)d射程范圍內(nèi)的動(dòng)壓預(yù)測能力。

(3)Realizablek-ε模型模擬高壓淹沒水射流，模型常數(shù)σk對(duì)模型預(yù)測軸心動(dòng)壓的影響較大，而模型常數(shù)σε的影響較小。

(4)應(yīng)用修正模型常數(shù)的Realizablek-ε模型模擬高壓淹沒水射流，σk與射流壓力的關(guān)系滿足反比例型函數(shù)的關(guān)系，射流壓力越大，σk越小。