基于格子玻爾茲曼方法的近壁區(qū)空化泡生長潰滅過程模擬

冉彬君,汪 磊,袁 浩,何小瀧

(1.重慶交通大學(xué)航運與船舶工程學(xué)院,重慶 400074; 2.百色樞紐通航投資有限公司,廣西 百色 533000;3.重慶交通大學(xué)西南水利水運工程科學(xué)研究院,重慶 400074)

空化泡初生是指當(dāng)液體內(nèi)部的局部壓力降低或局部溫度升高時,氣核附近壓力低于飽和蒸汽壓,氣核生長為空化泡。隨著周邊壓力環(huán)境恢復(fù),空化泡發(fā)生潰滅,在潰滅瞬間,由于短時相變引發(fā)高溫高壓環(huán)境,伴有強烈的沖擊波和高速微射流[1-2]。這種現(xiàn)象常常伴隨著噪聲、振動和空化侵蝕,對水利水運機械的運行產(chǎn)生重大影響[3-4],導(dǎo)致設(shè)備性能下降。盡管研究者們開展了大量且系統(tǒng)的研究,但受限于試驗觀測手段、近壁區(qū)空化泡的非對稱演化、宏觀數(shù)值模擬中熱流耦合技術(shù)的不成熟等,仍有復(fù)雜的水動力學(xué)和熱動力學(xué)機理尚未被揭示。

空化泡試驗研究主要基于高速攝像技術(shù)或者電火花、激光等[1,5],但由于測量技術(shù)的局限性,難以獲得空化泡周邊流場與溫度場的連續(xù)演化過程。受輸入能量、輸入尺寸、邊界壓力等諸多因素的影響,不同研究者獲得的空化泡尺寸、潰滅時間、最大流速和壁面壓力均存在差異,可重復(fù)性不高。隨著計算機技術(shù)的發(fā)展,計算流體力學(xué)(CFD)方法已經(jīng)成為空化研究的重要工具,與試驗相比,數(shù)值模擬可以實現(xiàn)不同條件組合甚至極端條件下空化泡的演變過程。基于玻爾茲曼動理學(xué)的格子玻爾茲曼方法(LBM)作為一種介觀研究方法,是模擬多相現(xiàn)象的一種穩(wěn)定有效的工具。與傳統(tǒng)CFD方法相比,LBM具有計算效率高、能夠處理復(fù)雜邊界、可通過狀態(tài)方程直接求解壓力等優(yōu)點[6-8]。

目前基于偽勢LBM的空化模擬大多僅考慮了空化泡的潰滅過程,但初始流場調(diào)整時間往往占據(jù)著整個演化過程的1/3[9],因此流場調(diào)整帶來的誤差往往難以簡單消除。本文利用初始輸入高溫激發(fā)氣核,采用He等[10]構(gòu)建的LBM熱流耦合模型模擬了激光/電火花空化泡的初生、生長和潰滅的全過程,并著重探討了泡壁距離和初始激發(fā)溫度對流場、最大半徑演化的影響,以期為深入揭示空化泡潰滅過程中的復(fù)雜規(guī)律打下基礎(chǔ)。

1 數(shù)值模擬方法

1.1 LBM水動力模型

考慮大密度比水氣多相流現(xiàn)象的數(shù)值穩(wěn)定性,水動力部分模擬采用基于多松弛碰撞算子的LBM偽勢模型,其粒子分布函數(shù)為

fi,eq(x,t))-(I-0.5Λ)SΔt

(1)

其中

(2)

(3)

當(dāng)引入非理想流體狀態(tài)方程后,ψ(x)可表述為

(4)

其中

式中:cs為格子聲速;pEOS為壓力,本文選取Carnahan-Starling(C-S)狀態(tài)方程計算:

(5)

式中:a、b、R為狀態(tài)方程控制參數(shù),a=1,b=4,R=1;T為節(jié)點對應(yīng)溫度。

1.2 LBM熱動力模型

本文采用Li等[12]提出的LBM溫度場演化計算方法以計算空化泡生長潰滅時的熱動力過程,其粒子分布方程為

gi(x+eiΔt,t+Δt)=gi(x,t)-

(6)

式中:gi為用于計算溫度場的粒子分布函數(shù);gi,eq為其對應(yīng)的粒子平衡態(tài)分布函數(shù);O為源項,可表述為

O=(O0,0,0,0,0,0,0,0,0)T

(7)

其中

式中:φ為宏觀方程中對應(yīng)的修正項,其一階和二階偏導(dǎo)數(shù)采用二階中心差分格式求解;k為熱傳導(dǎo)系數(shù);α為熱擴散系數(shù),α=k/ρcv;cv為流體比熱。本文液相熱擴散系數(shù)αl取0.15,而氣相熱擴散系數(shù)αg取0.45,各計算點上熱擴散系數(shù)取值按照線性插值求得:

(8)

式中ρg、ρl分別為氣相密度和液相密度。本研究中液相比熱cv_l為9.0,氣相比熱cv_g為3.0,各計算點上比熱取值按照線性插值求得:

(9)

利用矩陣M,可以將式(6)投影到矩空間:

(10)

式中:n為gi在矩空間上的投影;neq為gi,eq在矩空間上的投影:

neq=T(1,-2,2,ux,-ux,uy,-uy,0,0)T

(11)

(12)

(13)

式中:τj、τq為Λ中的松弛系數(shù);n4,neq、n6,neq為n4和n6的非平衡部分,可以通過ni,neq=(ni-ni,eq)/Δt求得。

1.3 熱流耦合過程

前人研究空化泡潰滅時的熱動力過程多采用被動標(biāo)量法[11],并未考慮熱動力過程對流場的影響。研究[8]表明,空化泡的演化伴隨著溫度的改變,其中表面張力隨著溫度的升高而降低,將直接影響空化泡的演化形態(tài)。因此本文考慮了溫度場與流場的耦合,采用t時刻的宏觀密度ρ(x,t)和流速u(x,t)求解t時刻溫度場中的平衡態(tài)分布函數(shù)neq(x,t),進而得到t時刻的溫度T(x,t);而下一時刻的壓力pEOS(x,t+Δt)和偽勢ψ(x,t+Δt)則采用求解T(x,t),進而得到下一步的宏觀密度ρ(x,t+Δt)與流場u(x,t+Δt)。值得注意的是,本文分析均采用格子單位,其中長度單位為lu,質(zhì)量單位為mu,時間單位為ts,溫度單位為tu。

2 模型驗證

模型驗證采用無限域內(nèi)單空化泡的初生和潰滅過程,并與模型試驗和Rayleigh-Plesset(R-P)方程的理論解進行定性與定量對比。考慮熱效應(yīng)的二維R-P方程[13]可描述為

(14)

其中

式(14)可通過四階Runge-Kutta法求解,其中每一步的氣相壓力為對應(yīng)的LBM計算獲得的空化泡內(nèi)平均壓力。計算域初始激發(fā)溫度場為

(15)

式中:Tl為初始液相溫度;Tv為初始氣相溫度;(x0,y0)為初始高溫中心;r0為初始圓形高溫區(qū)域半徑;s為高溫到低溫區(qū)域過渡距離。

計算域中心初始化一個半徑r0=5lu、溫度Tv=Tc的高溫區(qū)域。通過高溫區(qū)域瞬間汽化產(chǎn)生氣核。初始邊界壓力和計算域內(nèi)液相壓力設(shè)置為pl=-2.18×10-2mu/(lu·ts2),邊界壓力則設(shè)置為時間函數(shù),如圖1所示。當(dāng)邊界壓力低于平衡態(tài)壓力時,空化泡處于生長過程,而當(dāng)邊界壓力高于空化泡內(nèi)壓力時,空化泡則處于潰滅過程。此外采用無量綱時間t*=t/tc對比LBM模擬結(jié)果與試驗結(jié)果,其中tc為空化泡潰滅時間。

圖1 邊界壓力隨時間演化過程

LBM模擬結(jié)果與試驗結(jié)果[14]中對應(yīng)的無量綱時間的空化泡半徑與最大半徑比如表1所示,模擬結(jié)果表明空化泡在生長和演化過程中均保持圓形,兩者演化過程中的無量綱半徑比例接近。此外,模擬結(jié)果顯示在空化泡生長過程中,主要發(fā)生液相向氣相方向的相變,伴隨著吸熱過程,導(dǎo)致空化泡內(nèi)溫度低于周邊液體;而在潰滅過程中,隨著液體周邊壓力大于氣泡內(nèi)壓力,主要發(fā)生氣相向液相方向的相變,并伴隨著放熱,空化泡內(nèi)溫度逐漸升高。

圖2為LBM模擬結(jié)果與R-P理論方程獲得的空化泡半徑演化過程的對比,結(jié)果表明兩者吻合良好,僅在空化泡初生階段出現(xiàn)較小誤差,這可能是由于在LBM模擬中初始階段流場和溫度場的調(diào)整導(dǎo)致。

圖2 LBM數(shù)值模擬結(jié)果與R-P理論解獲得的空化泡半徑演化過程對比

3 結(jié)果與討論

為研究近壁區(qū)空化泡的初生和潰滅過程,選取500lu×500lu矩形計算域,溫度初始化如式(15)所示,在氣核激發(fā)點初始化一個半徑為r0=5lu的高溫區(qū)域。只有當(dāng)初始激發(fā)溫度高于流體當(dāng)前壓力下的臨界溫度時,才能實現(xiàn)液相轉(zhuǎn)化為氣相,形成氣核。空化泡的生長通過邊界及計算域內(nèi)的流體低壓實現(xiàn),而潰滅則通過邊界高壓傳遞到計算域內(nèi)實現(xiàn)。選取5種初始激發(fā)溫度探究其對空化泡的影響,分別為1.0Tc、1.2Tc、1.4Tc、1.6Tc和1.8Tc。計算域邊界為非平衡壓力邊界,底部邊界為無滑移邊界。計算域溫度邊界為非平衡外推常溫邊界,其頂部邊界和左右兩側(cè)邊界設(shè)置為T∞=0.5Tc,壁面溫度同樣設(shè)置為Tw=0.5Tc。空化泡內(nèi)的初始壓力和邊界壓力與第二節(jié)設(shè)置相同。

3.1 泡壁距離的影響

設(shè)置初始泡壁距離d=10～180lu,其對應(yīng)的泡壁無量綱距離為γ=d/rmax=0.19～2.94。圖3顯示空化泡最大半徑rmax隨著無量綱泡壁距離的增大而增大。當(dāng)γ<1.6時,空化泡生長過程中底部受壁面邊界限制,當(dāng)空化泡達到最大半徑時無法保持圓形,空化泡最大半徑與泡壁無量綱距離之間關(guān)系可表示為rmax∝γθ,在本研究中θ的取值范圍為0.028～0.038。而當(dāng)γ≥1.6時,空化泡生長始終保持圓形,空化泡最大半徑rmax則與γ之間存在線性關(guān)系。圖4為壁面最大流速uwmax、最大壓力pwmax和最高無量綱溫度Twmax/Tc隨泡壁無量綱距離的變化。由于γ=0.19時空化泡在潰滅點即發(fā)生數(shù)值發(fā)散,因此無法捕捉到壁面最高溫度,這里僅用空心點代替。由于泡壁無量綱距離增加,當(dāng)γ<1時,空化泡潰滅后產(chǎn)生的微射流、壓力波和高溫將直接作用于壁面上,而隨著泡壁無量綱距離的增加,壁面受到的沖擊壓力和沖擊流速迅速減小。當(dāng)γ≥2.2時,由于空化泡處于泡壁影響區(qū)外,壁面受到的沖擊壓力和最大流速逐漸趨于一致。

圖3 泡壁無量綱距離和空化泡最大半徑的關(guān)系

圖4 Tini=1.2Tc時壁面最大流速、最大壓力和最高無量綱溫度隨泡壁無量綱距離的變化

3.2 初始激發(fā)溫度的影響

初始輸入能量可以表示為

(16)

式中A為計算域面積。初始輸入能量與空化泡最大半徑之間的關(guān)系如圖5所示,對于附壁型空化泡(0.18≤γ≤0.19)和近壁區(qū)游離型空化泡(1.55≤γ≤1.59,2.86≤γ≤2.92),與試驗研究結(jié)論相同,空化泡最大半徑與輸入能量之間均為線性關(guān)系。這意味著可通過調(diào)節(jié)輸入能量來獲得可控大小的空化泡。圖6展示了雙對數(shù)坐標(biāo)下無量綱輸入溫度Tini/Tc與空化泡最大半徑rmax之間的關(guān)系,兩者大致為冪函數(shù)關(guān)系,即rmax∝(Tini/Tc)θ,不同γ條件下的θ非常接近(對應(yīng)圖6的斜率),其值約為0.08。

圖5 初始輸入能量與空化泡最大半徑的關(guān)系

圖6 空化泡最大半徑與初始輸入溫度的關(guān)系

圖7為d=90lu時不同初始激發(fā)溫度條件下空化泡潰滅的形態(tài)變化,其中左側(cè)為空化泡初生過程,右側(cè)為空化泡潰滅過程。在空化泡潰滅的最終階段,由于泡內(nèi)溫度隨著初始激發(fā)溫度的升高而增加,導(dǎo)致高溫處的表面張力更小,空化泡抵抗非對稱形變的能力變?nèi)酢Ｗ罱K導(dǎo)致空化泡在更大的半徑時頂部即產(chǎn)生凹陷,形成更為集中的微射流。

圖7 d=90lu時不同初始激發(fā)溫度條件下空化泡形態(tài)演化

空化泡演化數(shù)值模擬研究中,當(dāng)流體工作溫度遠低于其臨界溫度時,通常將空化過程看作絕熱過程,不考慮空化泡的熱動力過程對水動力過程的影響,在這種情況下,可采用被動標(biāo)量法對空化問題進行研究[15-17]。但對于工作溫度處于臨界溫度附近的低溫流體,其熱動力過程尤為重要,空化泡演化過程中的溫度變化往往導(dǎo)致界面上的表面張力處于持續(xù)變化的狀態(tài),因此有必要考慮熱動力過程對水動力過程的影響。因此本文采用被動標(biāo)量模型與熱流耦合模型進一步對比了He等[10]的模擬結(jié)果,如圖8所示,其中左側(cè)為生長階段,右側(cè)為潰滅階段。對于被動標(biāo)量模型,在模擬中僅考慮流場對溫度場的影響,而并不考慮溫度場對流場的影響。對于被動標(biāo)量模型,其計算壓力場時式(5)中溫度本文取為常溫T=0.5Tc,因此其界面上的表面張力則保持常數(shù)σ=0.0102。模擬結(jié)果表明,當(dāng)空化泡達到最大半徑時,空化泡內(nèi)與邊界上的壓力差極小,但對于熱流耦合模型,在空化泡半徑生長過程中,邊界壓力持續(xù)向計算域內(nèi)做功,導(dǎo)致部分機械能轉(zhuǎn)變?yōu)閮?nèi)能,空化泡周邊環(huán)境溫度高于0.5Tc,導(dǎo)致其表面張力較被動標(biāo)量模型中更小。根據(jù)r=σ/Δp,表面張力越大,對應(yīng)最大半徑也越大,因此被動標(biāo)量模型獲得的最大半徑更大。同時在潰滅過程中由于泡內(nèi)氣體發(fā)生相變轉(zhuǎn)化為液體,并伴隨著放熱,其溫度遠高于0.5Tc,最終導(dǎo)致其射流體積更大,空化泡最終被擊穿。

圖8 不同模型獲得的空化泡形態(tài)演化對比(d=90lu,Tini/Tc=1.4)

不同模型獲得空化泡的潰滅時間、最大半徑和空化強度的對比如表2所示。由于被動標(biāo)量模型能獲得更大的空化泡體積,且其表面張力更大,最終導(dǎo)致其潰滅時間比熱流耦合模型獲得的潰滅時間大9.8%。空化泡潰滅強度隨著表面張力增大而增大[15],因此被動標(biāo)量模型獲得的空化泡潰滅時最大流速和最大壓力比熱流耦合模型分別大3.61%和6.52%,而潰滅溫度則高9.3%。

表2 不同模型計算獲得的空化泡潰滅時間、最大半徑和潰滅強度(d=90lu,Tini/Tc=1.4)

4 結(jié) 論

a.受壁面限制,當(dāng)泡壁無量綱距離γ<1.6時,空化泡生長過程中無法保持圓形,其最大半徑與泡壁無量綱距離之間為冪函數(shù)關(guān)系。當(dāng)γ≥1.6,壁面影響減小,空化泡最大半徑與γ之間為線性關(guān)系。

b.空化泡最大半徑與輸入能量之間為線性關(guān)系,與初始輸入無量綱溫度之間為冪函數(shù)關(guān)系,且θ值并不隨泡壁無量綱距離的改變而改變。

c.熱流耦合模型考慮了潰滅時高溫對水氣界面上表面張力帶來的影響,與被動標(biāo)量模型相比,在潰滅階段獲得的射流體積更大,微射流更集中,但潰滅強度小于被動標(biāo)量模型。