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

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

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

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

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

對于淹沒水射流這種復雜情況的流動來說，數值模擬是一個十分有效的方法,但其不能完全替代實驗，數值計算的結果一般都需要驗證，并用規范的實驗來確認[8]。

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

1 淹沒水射流實驗

1.1 淹沒水射流結構

一般把淹沒水射流分為初始段和主體段。在噴嘴出口處，射流與環境水發生劇烈的動量交換和湍動擴散，其中心線附近部分射流介質仍會保持噴嘴出口時的速度，該部分介質稱為等速核心也稱為勢流核，射流初始段由勢流核與介質交換的混合區組成。該區域射流速度大，其與環境水會造成強烈的剪切作用產生空化現象。勢流核消散后的區域稱為射流主體段，此時射流的湍動性增強，進一步卷吸環境介質，射流軸心速度迅速降低，在射流核心區產生的空化氣泡會在這個區域內發展增大，并隨著射流運動到一定距離后湮滅。空化現象涉及復雜多變的物理現象，如氣體與液體之間的兩相轉變，其對射流速度等特性影響不可忽略。

1.2 水射流實驗

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

圖1 水射流實驗設備Fig.1 The equipment of water jet experiment

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

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

圖2 軸心滯止壓力的實驗值Fig.2 Experimental results of axial stagnationpressure

2 數值模擬

2.1 模型的控制方程

湍流模型選擇Realizablek-ε模型來對流場進行模擬計算。與標準k-ε模型相比，Realizablek-ε采用了新的湍流黏度公式，公式中的系數Cμ不再為一個常量，而是與湍流層流應變和旋度等有關的函數，而且為湍動能耗散率方程增加了新的傳輸公式。這使Realizablek-ε能更精確地預測平板和圓柱射流的發散性，而且對于預測旋轉流動、強逆壓梯度的邊界層流動、流動分離和二次流也有很好的表現[14]。

Realizablek-ε模型的控制方程為

Gk+Gb-ρε-YM+Sk

(1)

(2)

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

(3)

(4)

式(4)中：系數C1與時均應變率Sij聯系在一起，由式(5)求得，即

(5)

式(5)中：

(6)

Gk為湍動能的產生項，表示為

(7)

對式(7)用Boussinesq假設處理，Gk的表達式可改為

Gk=μtS2

(8)

(9)

在標準k-ε模型中系數Cμ為常數，而在Realizablek-ε模型Cμ按式(10)確定，即

(10)

式(10)中：

(11)

(12)

(13)

(14)

(15)

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

2.2 模型與邊界條件

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

圖3 噴嘴結構Fig.3 Structure of nozzle

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

圖4 網格劃分Fig.4 Gird partition

圖5 噴嘴網格加密Fig.5 nozzle mesh refinement

3 模擬結果及其對比

3.1 仿真結果及其分析

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

圖6 無因次動壓與射程關系Fig.6 Relationship between dimensionless dynamic pressure and range

3.2 對模型常數的改良

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

3.3 σk和σε對計算結果的影響

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

表1 σk和σε取值

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

(16)

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

從圖9可以看出，射流壓力值5～15 MPa時，σk的變化較大，15～40 MPa時σk的變化較小。這是因為從低壓力到高壓力時，淹沒射流的空化云體積的增長速率隨著壓力增大是先增大后減小的。這與毛寧等[16]用高速攝影技術觀察的空化現象一致。

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

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

圖9 湍流Prandtl數和射流壓力的關系Fig.9 Relationship between turbulent Prandtl numbers σk and jet pressure

3.4 模擬結果

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

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

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

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

圖13 沿程動壓分布對比Fig.13 Dimensionless dynamic pressure distribution of group

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

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

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

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

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

4 結論

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

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

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

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