單股射流傾斜撞壁影響因素的數(shù)值模擬

戴青雯,楊建文,楊 斌,王 瑩

(1.上海理工大學 能源與動力工程學院 上海市動力工程多相流動與傳熱重點實驗室,上海 200093;2.西安航天動力研究所,陜西 西安 710100)

0 引言

傾斜射流撞擊固體表面設計簡單,并且生產(chǎn)成本低,廣泛應用于燃料霧化[1]、液體火箭發(fā)動機[2]的液膜冷卻[3-4]、射流清洗[5-7]等領域。在液體火箭發(fā)動機工作過程中,核心區(qū)域溫度高達3 000 ℃,為避免燃燒室壁面被高溫熱氣流燒蝕,常采用液膜冷卻的方式來降低傳遞到壁面的熱量,達到保護壁面的作用[8]。在液膜冷卻設計中,發(fā)動機推力室內冷卻液膜的覆蓋面積是影響冷卻性能的重要指標,因此有必要對液膜形態(tài)的影響因素開展深入的研究。

文獻[9]的研究表明,射流角度、射流速度、流體物性及壁面接觸角等因素皆會影響射流撞壁后的液膜形態(tài),例如液膜呈現(xiàn)圓形水力跳躍、辮狀流、撞壁后反彈及先鋪展后破碎等形態(tài)。針對傾斜射流撞壁后的液膜流動這一典型氣液兩相流問題,國內外很多學者通過開展一系列的實驗來研究液膜的形態(tài)變化規(guī)律。文獻[10]針對疏水和超疏水壁面開展了射流撞壁的實驗研究,發(fā)現(xiàn)接觸角增大,而液膜鋪展面積減小,同時還對比分析了壁面接觸角、射流傾角及韋伯數(shù)對射流反射角的影響。Good等探究了不同射流傾角對撞擊后親水表面液膜分布的影響,揭示了射流速度與液膜最大寬度和壁面液膜濺射量之間的關系[11]。唐亮等在開展單股圓柱射流撞壁的實驗研究的基礎上,得出了射流撞壁后液膜最大寬度、液膜長度等幾何參數(shù)的表達式[12]。林慶國等開展了射流撞擊曲壁的實驗研究,采用探針法測量液膜厚度,探究了射流角度和壁面曲率半徑對液膜厚度和液膜形態(tài)的影響[13]。然而,射流撞壁后液膜的鋪展過程復雜且實驗獲取液膜幾何參數(shù)的成本較高。

在實驗研究的基礎上,國內外很多學者通過建立半經(jīng)驗模型開展液膜變化規(guī)律的理論研究。Wang等建立了半經(jīng)驗模型來分析射流速度、射流角度、流體黏度和表面張力等因素對液膜關鍵特征的影響,并對射流撞壁開展了一系列實驗,得出的半經(jīng)驗模型計算結果與實驗結果相符合[14]。Hasson等提出兩股射流對撞的流動滯止點位于對撞后形成的橢圓形液體薄片的一個焦點上,并在不考慮邊界層理論和黏性作用的情況下建立了液膜厚度分布的理論模型[15]。Watson考慮邊界層理論,建立了射流垂直撞擊壁面后層流和湍流兩種情況下的液膜厚度分布模型[16]。更進一步地,Inamura等[17]考慮到層流邊界層的發(fā)展,對Hasson等[15]的模型進行修正后建立了適用于射流撞擊曲壁的液膜厚度分布模型。文獻[18]在考慮黏性作用的基礎上建立了單股圓柱射流撞擊壁面后的液膜鋪展模型,在實驗驗證模型可靠性的基礎上研究了表面親疏水性對液膜形態(tài)的影響。為進一步深入研究液膜的形態(tài)變化,文獻[19]還建立了兩股流體射流的理論分析模型,研究了射流角度、韋伯數(shù)和雷諾數(shù)對液膜形態(tài)的影響。然而,傳統(tǒng)經(jīng)驗模型目前還無法開展下游流動復雜處的壁面壓力和液膜變化規(guī)律的研究。

相比于實驗和理論模型預測,傾斜射流撞壁在數(shù)值仿真方面的研究相對較少。對于垂直射流撞擊壁面,Gradeck等開展了不同射流速度下撞擊運動表面的實驗和數(shù)值仿真研究,通過對比數(shù)值仿真結果與實驗結果得到的水躍位置,驗證了數(shù)值模型的準確性[20]。在傾斜射流撞擊壁面方面,Fard等運用三維數(shù)值模擬研究了工業(yè)噴嘴射流撞壁后液膜的形成及破碎成液滴過程,并將仿真得到的液膜厚度與實驗結果進行對比,驗證了數(shù)值模型的可靠性[21]。林慶國通過數(shù)值仿真得到壁面上的液膜分布,研究了不同壁面曲率和不同壓降對液膜表面形態(tài)的影響,且只定性研究了液膜的基本形態(tài),并未定量對比數(shù)值仿真結果與實驗結果[22]。邱添開展了射流沖擊平板的全三維和準三維數(shù)值模擬研究,由于其全三維模擬研究的網(wǎng)格不夠細密,使得模擬得到的液膜形態(tài)存在失真現(xiàn)象[23]。上述文獻并未對不同影響因素,例如不同射流速度下壁面的壓力變化和不同孔徑下液膜表面速度變化開展相關研究。

對于單股射流撞壁后形成的液膜,可以用液膜長度、寬度和厚度等關鍵幾何參數(shù)來描述。但當前在液膜幾何參數(shù)方面的定量研究較少,一是實驗獲取液膜幾何參數(shù)的耗時長、成本高,二是傳統(tǒng)經(jīng)驗模型會忽略黏性、重力等因素對液膜的影響導致不能精準捕捉液膜的基本形態(tài)。因此,本文采用數(shù)值模擬方法定量研究射流撞壁后的流場區(qū)域。首先,比較實驗和數(shù)值模擬所得的液膜寬度,驗證了本文采用的VOF多相流模型及數(shù)值模擬方法的可靠性;其次,分析液膜形成過程中的關鍵特征;最后,定量研究射流傾角、射流速度及射流孔徑對液膜形態(tài)的影響。

1 數(shù)值求解方法

1.1 控制方程

在本文射流撞壁數(shù)值仿真的研究中,研究對象為水和空氣。在計算中,流體的物性設為常數(shù),因此連續(xù)性方程和動量方程為

(1)

(2)

式中:ρ為流體密度;v為流體速度;t為時間;p為壓強;μ為動力黏度;g為重力加速度;f為動量源項。

在數(shù)值模擬過程中,氣液兩相的體積分數(shù)相加和為1,引入流體體積比函數(shù)φ的概念,φ表示流體體積與網(wǎng)格單元體積的比值[24]。φ=1時,表征區(qū)域內只有液態(tài)水;φ=0時,表征區(qū)域內只有空氣;0<φ<1時,表征區(qū)域內既有液態(tài)水也有空氣。

因此,密度和動力黏度可用流體體積比函數(shù)φ進行插值計算,即

(3)

(4)

式中:ρl和ρg分別為液體和氣體的密度;μl和μg分別為液體和氣體的動力黏度。

1.2 幾何模型

在傾斜射流撞壁的過程中,液態(tài)水從圓形噴孔以一定的角度射出,到達壁面后液態(tài)水先在平板表面向四周鋪展并繼續(xù)向出口處流動,然后在表面張力的作用下收縮,形成一個類似橢圓形的液相區(qū)域,最后匯聚成一股流動從出口流出,射流撞壁后形成的液膜如圖1所示[14]。

圖1 液膜圖像Fig.1 Picture of the liquid sheet

由于液態(tài)水射流撞壁鋪展成液膜的全過程是關于噴孔中沿軸線對稱的,計算域模型采用半模型,可以減小計算量;同時液膜厚度很小,在不影響液膜形態(tài)的前提下可以減小計算域模型在液膜厚度方向上的尺寸。為了充分展示流體從射流入口以一定的速度v沿著射流傾角θ方向流出后液膜的鋪展情況,故計算域大小為9 mm×51 mm×6 mm,如圖2所示。圖2中,y軸正方向為重力方向,射流傾角θ為速度與y軸之間的夾角,同時液膜主流方向與y軸正方向保持一致。

圖2 計算域三維示意圖(單位:mm)Fig.2 Three-dimensional schematic chart of computational domain(unit:mm)

1.3 計算方法

選擇基于有限體積法的Fluent求解器計算,采用的VOF模型是一種基于歐拉法的多相流表面跟蹤方法[25],具有易守恒和占用內存小的優(yōu)點,適用于相界面變化復雜的液體流動計算。使用VOF模型計算傾斜射流撞壁時,在多相流模型設置中將空氣和液態(tài)水分別設置為第一和第二相,初始計算域內均為空氣。

將射流入口設置為速度入口,壁面設置為無滑移壁面,與壁面的接觸角設為73.3°,其他各個面都設置為壓力出口。物性參數(shù)設置參考水的物性參數(shù),表面張力系數(shù)設為72.8 mN/m,液體的密度設為998.2 kg/m3,黏性系數(shù)設為0.001 Pa·s。湍流模型選擇RNGk-ε模型,采用Menter-Lechner壁面函數(shù)處理近壁面流動,壓力離散選用PRESTO算法,壓力與速度耦合關系選用SIMPLE算法,時間步長為10 μs,每時間步長最大迭代步數(shù)為40。

1.4 網(wǎng)格劃分和無關性驗證

采用多面體網(wǎng)格劃分方法,并且為了更好地捕捉液膜形態(tài),在壁面附近設置邊界層,如圖3所示。通過4套不同尺寸的網(wǎng)格進行網(wǎng)格無關性驗證,網(wǎng)格數(shù)量分別為54萬、99萬、115.5萬和175萬。參考文獻[26]取射流孔徑為0.8 mm,設置進口速度為5 m/s,射流傾角為20°,以y=6 mm、y=12 mm和y=20 mm這3個截面處的液膜寬度為判斷依據(jù)。

圖3 計算域網(wǎng)格及壁面邊界層網(wǎng)格局部放大圖Fig.3 Computational domain mesh and partial magnified drawing of wall boundary layer mesh

圖4為在射流孔徑為0.8 mm、速度為20 m/s和傾角為20°的工況下,射流撞擊固體表面后形成的液膜形態(tài)及液膜關鍵幾何參數(shù)示意圖。射流撞壁后可分為自由射流區(qū)、薄膜區(qū)、邊緣水躍區(qū)及液膜匯集區(qū)。

圖4 射流撞壁后的液膜形態(tài)及液膜關鍵幾何參數(shù)示意圖Fig.4 Diagram of liquid film shape after jet impingement and key geometric parameters of liquid film

圖4還展示了液膜的關鍵幾何參數(shù)液膜長度l和液膜寬度w的提取方式。液膜長度為沿著流動方向從壁面上開始出現(xiàn)液膜到液膜第一次匯集成單股流動的位置;液膜寬度為液膜兩側邊緣在某一截面上的鋪展距離。

網(wǎng)格無關性驗證結果如圖5所示。當網(wǎng)格數(shù)量為54萬時,這3個截面上的液膜寬度與其他3套網(wǎng)格液膜寬度的相對誤差超過5%,而99萬、115.5萬和175萬3套網(wǎng)格之間計算出來的液膜寬度相對誤差均低于5%。為兼顧計算精度和成本,最終應用115.5萬的網(wǎng)格數(shù)量進行后續(xù)的仿真計算。

圖5 網(wǎng)格無關性驗證結果Fig.5 The result of mesh independence verification

1.5 模型驗證

當液態(tài)水以射流的形式撞擊固體壁面,由于慣性力的作用,液態(tài)水以撞擊點為中心向周圍擴散,并形成兩個突出的邊緣,即邊緣水躍區(qū);在邊緣水躍區(qū)之間形成的一層薄液體層稱為薄膜區(qū);在液膜達到最大寬度后,由于表面張力克服慣性力的作用,邊緣水躍區(qū)開始向內收縮流動,最終在液膜匯集點變成一股流體繼續(xù)向下游流動;在液膜匯集點之后的區(qū)域稱為液膜匯集區(qū),至此,液膜匯集區(qū)上游的液膜穩(wěn)定發(fā)展,液膜鋪展的關鍵特征不再發(fā)生變化。圖4中觀察到的現(xiàn)象與文獻[14]實驗結果相似,證明了本文模型設置的準確性。

為了進一步驗證數(shù)值仿真結果的可靠性,將數(shù)值仿真結果與實驗結果進行定量對比。當液膜穩(wěn)定后,在y=6 mm、y=12 mm和y=20 mm這3個截面上提取液膜寬度,并與文獻[26]實驗得到的液膜寬度進行對比,如表1所示。水射流撞擊平板后形成液膜,液膜與壁面之間存在一定的相互作用,而在本文數(shù)值仿真過程中簡化了液膜與壁面之間溫度變化和壁面接觸角變化的相互作用關系,這使得仿真得到的液膜寬度略小于實驗值。但數(shù)值模擬結果與實驗結果之間的誤差均在20%以內,低于文獻[26]仿真結果的誤差。因此認為本文的數(shù)值仿真結果可靠,可采用本文的數(shù)值模擬方法進行射流撞壁的相關研究。

表1 不同截面的液膜寬度對比Tab.1 Width comparison of liquid film in different cross sections

2 結果與討論

2.1 射流傾角對液膜的影響

圖6為傾斜射流撞壁原理的局部示意圖。

圖6 傾斜射流撞壁原理的局部示意圖Fig.6 Partial schematic diagram of oblique jet impingement principle

Kate等通過理論預測和實驗驗證分析了射流撞壁后的流場變化,推導得出射流撞壁形成液膜的滯止點位置[27]。液態(tài)水從噴孔射出后,在與壁面平行的自由射流區(qū)截面A-A上呈橢圓形,橢圓的焦點S1在分離流線上,分離流線與射流中心線平行,且在分離流線右側射流向液膜匯集區(qū)方向鋪展,左側射流向反方向鋪展,將分離流線與壁面的交點稱為滯止點,如圖6中S所示。

在不同的射流傾角下,液膜的鋪展方式相同。首先液膜在鋪展過程中,沿著液膜鋪展方向,液膜寬度先增大后減小,但由于不同射流傾角下的液膜在各個方向上的流動分速度不同,使得液膜寬度和液膜厚度有所差異。

取滯止點所在截面為y=0 mm截面,圖7是射流傾角θ分別為20°、15°和10°時沿y方向不同截面位置的液相體積分數(shù)分布云圖。由圖7可知,射流傾角越大,液膜的寬度越大。雖然液膜傾角越大,在各個截面上的鋪展面積明顯增大,但是傾角越小,液膜的能量和受力越集中,在下游越不容易發(fā)生破碎。

圖7 沿y軸不同截面的液相體積分數(shù)分布云圖Fig.7 Liquid phase volume fraction in the y-direction at different cross sections

圖8對比分析了不同射流傾角下y=4 mm、y=10 mm、y=16 mm及y=22 mm這4個截面上的液膜厚度變化。射流傾角越大,不同橫截面上薄膜區(qū)與水躍區(qū)之間的液膜厚度變化越劇烈。這是由于射流傾角越大,液膜軸向分速度越小,流動越慢,橫向分速度越大,橫向擴展面積越大。此外,由于在不改變射流流量情況下改變傾角,液膜在薄膜區(qū)與邊緣水躍區(qū)厚度變化越劇烈,液膜穩(wěn)定性越差,且過大的射流傾角會使邊緣水躍區(qū)在橫向不斷向外擴展,在軸向不斷向內收縮,從而可能會導致后續(xù)液膜兩側破碎現(xiàn)象的產(chǎn)生。

圖8 不同射流傾角下的液膜厚度變化Fig.8 Thickness variation of liquid film at different inclination angles

結合圖7和圖8分析可知:當液膜完全鋪展并達到穩(wěn)定后,在薄膜區(qū),液膜厚度從中心軸線向兩側減小;當液膜鋪展到邊緣水躍區(qū)時,厚度突然增大。因此,可認為在液膜匯集區(qū)上游各個截面上液膜整體呈中間薄兩側厚的分布。

射流傾角越大,薄膜區(qū)鋪展面積越大,邊緣水躍區(qū)的厚度也越大;同時液膜的寬度越大,撞擊點與液膜匯集點之間的距離越長。文獻[23]中也觀察到同樣的現(xiàn)象,這是因為不同射流角度下流量在各個方向上分配不一致導致的,射流傾角越大,y方向流量越大,這也解釋了不同入射條件下撞壁后液膜輪廓不同的原因。

2.2 射流速度對液膜的影響

射流孔徑和射流傾角一定時,不同射流速度會影響液膜在流動方向上的延伸距離和液膜橫向鋪展方向上寬度的變化。表2是不同射流速度下液膜的長度和最大寬度的對比。

表2 不同射流速度下液膜最大寬度和液膜長度Tab.2 Maximum width and length of liquid film at different inclination velocity

圖9是不同射流速度下的液膜形態(tài),雖然射流速度不同,但描述液膜關鍵特征的自由射流區(qū)、薄膜區(qū)、邊緣水躍區(qū)及液膜匯集區(qū)的輪廓和邊界線都清晰可見。結合表2和圖9分析可得,射流速度越大,液膜的最大寬度越大,長度也越大。這是因為射流速度越大會導致射流總流量和動量越大,液膜的鋪展面積越大。

圖9 不同射流速度下的液膜形態(tài)Fig.9 Liquid film shape at different inclination velocities

液態(tài)水撞擊壁面后鋪展時,壁面會受到水的壓力。忽略壁面粗糙度、黏性等參數(shù)的影響,由于力的作用是相互的,該壁面壓力值也能反映液膜的壓力分布情況。圖10是不同射流速度下在0～150 Pa范圍內壁面的壓力分布。紅色區(qū)域表示此處壁面的壓力高于150 Pa,藍色區(qū)域表示此處壁面的壓力低于0 Pa。

圖10 不同射流速度下壁面的壓力分布Fig.10 Wall pressure distribution at different inclination velocities

從圖10中可以看出,液態(tài)水在剛撞擊壁面時的壓力較高,之后在薄膜區(qū)壁面壓力迅速降低,趨近于環(huán)境壓力;邊緣水躍區(qū)則是外邊界壓力最大,從水躍區(qū)外邊界到內邊界壓力遞減,外側壓力大能抑制薄膜區(qū)的持續(xù)鋪展,推動液膜向內側收縮;而在薄膜區(qū)和邊緣水躍區(qū)內邊界之間的分界處,壓力值小于0 Pa,認為該處壓力低于環(huán)境壓力。同時,在液膜匯集區(qū)之后的壁面會出現(xiàn)高壓,文獻[9]在進行超疏水壁面傾斜射流時的研究也發(fā)現(xiàn)液膜匯集區(qū)之后的壁面存在高壓區(qū)。因此,推測本文數(shù)值模擬的親水表面出現(xiàn)這一現(xiàn)象也是因為傾斜射流是由固體壁面支撐的,壁面壓力與射流的切向動量密切相關,液膜在匯集區(qū)發(fā)生堆積,之后在重力作用下回落到壁面并繼續(xù)向下游流動,這就使得匯集區(qū)之后的切向動量增大,從而匯集區(qū)之后的壁面出現(xiàn)高壓區(qū)。

圖10對比分析了不同射流速度時壁面的壓力分布,速度越大,撞壁后壓力高于150 Pa的面積越大,可認為射流撞壁區(qū)域越大;且速度越大,撞壁處液膜的慣性力越大,使得液膜在薄膜區(qū)鋪展面積越大,液膜寬度越大;當速度越小時,靠近液膜匯集區(qū)的水躍區(qū)壓力越大,且速度越小慣性力越小,表面張力越大,較高的表面張力推動著液膜收縮,使得液膜長度越短,但液膜的穩(wěn)定性會有所提高。

2.3 射流孔徑對液膜的影響

當保持射流傾角θ=20°和射流速度v=5 m/s時,射流孔徑越大,液態(tài)水的質量流量也越大,質量流量的差異導致液膜的長度、寬度和厚度發(fā)生改變。圖11是不同射流孔徑下液膜表面的速度分布,在撞擊點附近流速最大,使得液膜向周圍鋪展,而隨著液膜向下游運動,速度不斷減小,在水躍區(qū)匯集點附近速度達到最小。此時由于重力和氣體回流的影響,液膜下游初步呈現(xiàn)向內側聚攏的現(xiàn)象。

圖11 不同射流孔徑下液膜表面速度分布Fig.11 The velocity distribution of liquid film surface at different inclination diameter

不同射流孔徑下液膜寬度和長度變化如表3所示。由表3數(shù)據(jù)分析可知,射流孔徑越大,質量流量越大,液膜鋪展面積越廣,液膜的寬度和長度也越大。

表3 不同射流孔徑下液膜最大寬度和液膜長度Tab.3 Maximum width and length of liquid film at different inclination diameters 單位:mm

圖12是不同射流孔徑下中心軸線上的液膜厚度分布。從圖12中可以看出,從撞擊點處開始隨著液態(tài)水不斷地向下游流動,液膜厚度先逐漸減小再逐漸增大。發(fā)生這一現(xiàn)象的主要原因是,一開始射流撞擊壁面,較大的速度使得液膜迅速在壁面鋪展,液膜在撞擊點附近由厚變薄,之后由于壁面粗糙度及回流等因素的影響,流動速度變得緩慢,液膜逐漸開始收縮堆積,最后形成液膜匯集區(qū)。在這一過程中,中心軸線處的液膜厚度不斷增加,甚至在液膜匯集區(qū)的最大厚度與薄膜區(qū)的最小厚度相差一個量級。結合表3和圖12分析可知,射流孔徑越大,液膜最大寬度越大,液膜長度越長,液膜厚度越大。這與文獻[28]中理論模型計算得到液膜輪廓結論一致。表3中射流孔徑從0.6 mm增大至0.8 mm時,液膜最大寬度和長度的增幅分別為37.6%和37.7%,而孔徑變化對射流速度在各個方向的分量影響幾乎可忽略,可推測孔徑變化對液膜穩(wěn)定性的影響較小。

圖12 不同射流孔徑下中心軸線上液膜厚度變化Fig.12 Thickness variation of liquid film on central axis at different inclination diameters

3 結論

本文使用VOF多相流模型并結合多面體混合網(wǎng)格,對多工況下傾斜射流撞壁進行了數(shù)值模擬,并得出以下結論。

1) 本文數(shù)值仿真得到的液膜寬度與實驗值變化趨勢保持一致,且兩者之間誤差均在20%以內,表明本文構建的數(shù)值模型適用于傾斜射流撞壁的數(shù)值模擬研究。

2)在液膜寬度方向,液膜厚度整體呈中間薄兩側厚的分布;在液膜主流方向,從撞擊點處開始,中心軸線上液膜厚度先逐漸減小再迅速增大;射流孔徑越大,液膜的中心軸線上最大厚度也越大,當射流孔徑為0.6 mm時,液膜最大厚度約為1.3 mm,當射流孔徑為0.8 mm時,液膜的最大厚度超過2 mm。隨著射流傾角增大,液膜的寬度增大而長度減小;隨著射流速度和孔徑的增大,液膜的寬度和長度都增大。

3) 射流撞壁后,壁面壓力從液膜的邊緣水躍區(qū)到中心軸線先逐漸減小后增大再減小,且在薄膜區(qū)和邊緣水躍區(qū)內邊界之間的分界處,壓力值小于0 Pa。在撞擊點附近,射流速度越大,液膜的慣性力越大,使得壁面受到的壓力越大,而在液膜匯集區(qū),表面張力起主導作用,表面張力增大時慣性力減小,因而壁面受到的壓力越小。

4) 在液膜匯集點之前,沿著液膜長度方向,液膜的流動速度先增大后減小;液膜在撞擊點附近速度最大,近似等于射流速度,在水躍區(qū)匯集點附近速度最小;當射流孔徑從0.6 mm增大到0.8 mm時,液膜的最大質量流量增大,體現(xiàn)在液膜最大寬度從5.24 mm增大到8.4 mm,液膜長度從19.55 mm增大到31.55 mm,液膜鋪展面積越廣,液膜最大速度和最小速度區(qū)域的面積越大。