汽液黏滯系數(shù)對空化泡潰滅過程影響的數(shù)值模擬

袁曉龍，何小瀧,汪凱迪

(1.雅礱江流域水電開發(fā)有限公司，四川 成都 610051；2.交通運(yùn)輸部天津水運(yùn)工程科學(xué)研究所工程泥沙交通運(yùn)輸行業(yè)重點(diǎn)實(shí)驗(yàn)室，天津 300456；3.中國三峽建設(shè)管理有限公司，四川 成都 610023)

空蝕空化廣泛存在于水利、航海、液壓系統(tǒng)、人類體液循環(huán)系統(tǒng)這一系列工業(yè)與自然現(xiàn)象中[1-3]。空化泡潰滅過程中，由于空化泡內(nèi)部壓力與周邊液體壓力之間的壓差劇烈變化，泡內(nèi)汽體與泡外液體之間迅速轉(zhuǎn)化，因此，空化泡的潰滅過程是帶有復(fù)雜傳熱傳質(zhì)相變過程的流體力學(xué)問題。前人研究表明空化泡潰滅過程中產(chǎn)生的高溫、微射流和沖擊波的共同作用是空蝕破壞的主要原因。

常用的空化泡潰滅的研究方法有理論分析、試驗(yàn)研究和數(shù)值模擬。理論分析通常需要對問題作出許多簡化，如Plesset等[4]采用Raleigh方程模擬空化泡潰滅時(shí)作出了以下假設(shè)：忽略表面張力、液體可壓縮性、液體黏滯性等，這樣的假設(shè)導(dǎo)致理論分析并不能直接用于實(shí)際的空化泡潰滅的研究中。而試驗(yàn)受限于測量技術(shù)，難以獲得空化泡周邊的速度和壓力變化過程。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，計(jì)算流體力學(xué)因其結(jié)果可以有效展示流場的演化過程，逐漸成為研究流體運(yùn)動的重要手段。與試驗(yàn)研究相比，數(shù)值模擬可以有效地展示出空化泡在潰滅過程中微射流的形成。傳統(tǒng)的空化泡潰滅研究的數(shù)值模擬方法主要分為基于Navier-Stokes方程和非理想流體狀態(tài)方程的可壓縮空化模型、基于Navier-Stokes方程和傳質(zhì)方程的不可壓空化模型及可以模擬多空泡群的離散汽泡法。格子玻爾茲曼方法(LBM)作為一種介觀數(shù)值模擬方法，從玻爾茲曼粒子動理學(xué)出發(fā)，能夠有效模擬結(jié)晶、沸騰等帶有傳質(zhì)相變過程的多相流現(xiàn)象。相較于傳統(tǒng)數(shù)值模擬方法，LBM優(yōu)勢在于僅需要求解一系列線性方程組而避免了對流相的求解，同時(shí)壓力與狀態(tài)方程相關(guān)而避免求解泊松方程，且對于復(fù)雜邊界處理十分簡單[5]。

經(jīng)過近30年的發(fā)展，現(xiàn)有 LBM多相流模型可以歸納為偽勢模型[6-7]、顏色模型[8]、自由能模型[9-10]和相場模型[11]。其中偽勢模型由Shan等[6]提出，該模型利用粒子間的相互作用力自動形成交界面，簡化了汽液交界面的處理。Sukop等[12]利用偽勢模型對各相異性條件下空化泡的生成開展了研究。Chen等[13]結(jié)合偽勢模型與精確差分法(exact difference method，EDM)外力格式，模擬了液相與汽相的密度之比為67.5條件下空化泡在剪切流中的生長過程，其所得結(jié)果與Rayleigh-Plesset方程理論解吻合。Yang等[14-16]利用Bhantagar-Gross-Krood(BGK)模型與EDM外力格式相結(jié)合，模擬了液相與汽相的密度之比為30～50時(shí)單個(gè)空化泡的潰滅過程，獲得了潰滅過程中壓力場和速度場的變化。Shan等[17-18]進(jìn)一步采用多松弛時(shí)間(multi-relaxation-time， MRT)碰撞算子模擬了液相與汽相的密度之比為750時(shí)空化泡的潰滅過程，模擬結(jié)果與Lauterborn等[19]的理論分析結(jié)果吻合良好。以上研究均證明LBM可以有效模擬空化泡的生長和潰滅過程。

黏滯系數(shù)對空化泡潰滅過程中形成的微射流和沖擊波具有較大的影響，但當(dāng)前利用LBM對空化泡潰滅過程的研究中并未考慮黏滯系數(shù)對空化泡潰滅過程的影響，而是將汽液黏滯系數(shù)比通常設(shè)為1。本文采用LBM模擬近壁區(qū)空化泡潰滅過程，分析不同汽相黏滯系數(shù)和液相黏滯系數(shù)對空化泡潰滅過程的影響，以及潰滅過程中產(chǎn)生的最大微射流流速、最大壓力、最大潰滅時(shí)間的變化規(guī)律，為深入揭示空化泡潰滅機(jī)理打下基礎(chǔ)。

1 數(shù)值模擬方法

基于BGK碰撞算子[20]，LBM偽勢模型中粒子分布函數(shù)可以表示為

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

(1)

(2)

粒子平衡態(tài)分布函數(shù)可以表示為

(3)

式中：ρ為宏觀密度；ueq為平衡態(tài)分布速度；ωi為i方向的粒子平衡態(tài)分布函數(shù)的權(quán)重系數(shù)，對于D2Q9模型，當(dāng)i=0時(shí)ωi=4/9，當(dāng)i=1～4時(shí)ωi=1/9，當(dāng)i=5～8時(shí)ωi=1/36。流場實(shí)際速度可以通過下式求得：

(4)

式中：u為宏觀實(shí)際速度；FA為作用在粒子上的合力，包含粒子間相互作用力F和重力、向心力等體積力，F(xiàn)可采用Shan等[7]提出的公式計(jì)算：

(5)

式中：G為粒子之間作用強(qiáng)度；ψ為粒子之間的偽勢；wi為權(quán)重系數(shù)，對于D2Q9模型，當(dāng)i=0時(shí)wi=0，當(dāng)i=1～4時(shí)wi=1/3，當(dāng)i=5～8時(shí)wi=1/12[21]。根據(jù)Yuan等[22]的研究，引入非理想氣體狀態(tài)方程后，粒子間的相互作用勢可以表示為

(6)

式中：pEOS為壓強(qiáng)，可通過流體狀態(tài)方程求得，本文采用Carnahan-Starling(C-S)狀態(tài)方程[23]求解pEOS：

(7)

為滿足熱力學(xué)一致性，本文中外力采用Li等[24]提出的偽勢模型外力項(xiàng)，其表達(dá)式為

(8)

(9)

式中：ε為調(diào)節(jié)參數(shù)，用于滿足模型的熱力學(xué)一致性。密度場初始化如下：

(10)

式中：ρl為初始液相密度；ρg為初始汽相密度；R0為初始半徑；(x0,y0)為液滴中心；初始界面厚度w=5。

為研究流體黏滯系數(shù)對空化的影響，數(shù)值模擬中選取了不同松弛系數(shù)τ，如公式(11)所示[18]：

(11)

式中：τg為汽相松弛系數(shù)；τl為液相松弛系數(shù)；ρg、ρl分別為汽相密度和液相密度。

若無特別申明，文中所討論質(zhì)量單位為mu,長度單位為lu，時(shí)間單位為tu,速度單位為lu/tu，壓力單位為mu·lu-1·tu-2。

2 模型驗(yàn)證

2.1 熱力學(xué)一致性驗(yàn)證

為研究汽液黏滯系數(shù)比對模型的熱力學(xué)一致性的影響，采用單松弛等溫模型分別模擬了初始溫度T0= 0.7Tc、0.8Tc、0.9Tc、0.95Tc條件下的汽液分離過程，計(jì)算結(jié)果如圖1所示，圖中νg為汽相黏滯系數(shù)，νl為液相黏滯系數(shù)。由圖1可見，數(shù)值結(jié)果與理論解一致，在T≥0.7Tc條件下，改變汽液黏滯系數(shù)比不會影響數(shù)值模型熱力學(xué)一致性。

圖1 汽液黏滯系數(shù)比對熱力學(xué)一致性的影響

2.2 近壁區(qū)空化泡潰滅過程

為進(jìn)一步驗(yàn)證模型的準(zhǔn)確性，模擬了單個(gè)空化泡潰滅的演化過程，并與Philipp等[25]的試驗(yàn)現(xiàn)象進(jìn)行對比。計(jì)算域?yàn)?01 lu×401 lu正方形，空化泡初始半徑R0=50 lu，計(jì)算域底部固壁采用標(biāo)準(zhǔn)反彈邊界，左右均采用周期邊界，頂部采用Zou-He壓力邊界[26]，計(jì)算區(qū)域如圖2所示，圖中P為液相壓力，Pv為汽相壓力。

圖2 計(jì)算域示意圖

空化泡初始位置與試驗(yàn)相同，無量綱距離λ=d/R0，其中d為空化泡中心與壁面之間距離。汽相密度設(shè)置為0.009 mu·tu-3，液相密度設(shè)置為0.39 mu·tu-3，對應(yīng)空化泡內(nèi)外壓差Δp=0.003 89 mu·lu-1·tu-2，初始溫度T0=0.7Tc，同時(shí)假設(shè)初始空化泡位于靜止水體中，因此計(jì)算域內(nèi)初始流速均為0。

Philipp等[25]認(rèn)為無量綱距離λ>2.2時(shí)，壁面對空化泡潰滅的影響可以忽略。本文選取兩個(gè)無量綱距離λ=1.6和λ=2.5，與Philipp試驗(yàn)中使用的無量綱距離相對應(yīng)。λ=1.6和λ=2.5時(shí)計(jì)算結(jié)果如圖3和圖4所示，兩工況汽相松弛系數(shù)τg=1.062 5，液相松弛系數(shù)τl=0.187 5，對應(yīng)的汽液黏滯系數(shù)比νg/νl=15。當(dāng)λ=1.6時(shí)，受到壁面影響，近壁區(qū)在空化泡潰滅初始階段，橫向壓縮速率遠(yuǎn)大于縱向壓縮速率，在此階段汽泡主要為橢圓形。隨著空化泡潰滅的發(fā)展，空化泡頂部出現(xiàn)高壓區(qū)，而空化泡與壁面之間壓力較小，形成相對低壓區(qū)?？栈萆舷轮g形成壓力差，即Bjerknes力。在Bjerknes力的作用下，汽泡頂部出現(xiàn)凹陷，并伴隨著微射流的形成，最終空化泡發(fā)生潰滅。而當(dāng)λ=2.5時(shí)，空化泡處于壁面影響區(qū)外，由于空化泡與壁面之間距離過大，空化泡與壁面之間并不會出現(xiàn)相對低壓區(qū)，空化泡上下形成的壓差不足以使空化泡頂部出現(xiàn)凹陷。不同無量綱距離條件下，LBM數(shù)值模擬結(jié)果與Philipp等[25]的試驗(yàn)結(jié)果基本一致，驗(yàn)證了LBM模型的有效性。

圖3 λ=1.6時(shí)LBM模擬的空化泡形態(tài)演化過程

圖4 λ=2.5時(shí)LBM模擬的空化泡形態(tài)演化過程

3 汽相和液相黏滯系數(shù)變化的影響

不同汽相和液相黏滯系數(shù)可以通過調(diào)節(jié)松弛系數(shù)τ獲得。為研究汽相黏滯系數(shù)νg變化對空化泡潰滅過程中最大微射流流速umax、最大壓強(qiáng)pmax、最大潰滅時(shí)間tmax的影響，選取不同汽相松弛系數(shù)τg和相同的液相黏滯系數(shù)νl如表1所示。不同汽相黏滯系數(shù)對最大微射流流速、最大潰滅壓強(qiáng)及最大潰滅時(shí)間的影響如表1所示，可以看到改變汽相黏滯系數(shù)對空化泡潰滅時(shí)產(chǎn)生的最大微射流流速、最大壓強(qiáng)、最大潰滅時(shí)間毫無影響，這是由于汽體質(zhì)量較小，具有較小的慣性力。不同汽相黏滯系數(shù)條件下，空化泡潰滅時(shí)最大微射流流速為0.255 lu·tu-1，最大潰滅壓強(qiáng)為0.007 18 mu·lu-1·tu-2，最大潰滅時(shí)間為721 tu。

表1 汽相黏滯系數(shù)對空化泡潰滅的影響

為研究液相黏滯系數(shù)變化對空化泡潰滅過程中最大微射流流速、最大壓強(qiáng)、最大潰滅時(shí)間的影響，選取液相松弛系數(shù)和對應(yīng)的液相黏滯系數(shù)如表2所示。圖5展示了液相黏滯系數(shù)對空化泡潰滅過程中最大微射流流速、最大潰滅壓力和最大潰滅時(shí)間的影響，可以看出空化泡最大潰滅時(shí)間隨著液相黏滯系數(shù)的減小而減小，而空化泡潰滅最大微射流流速和最大潰滅壓強(qiáng)均隨著液相黏滯系數(shù)的減小而增大，當(dāng)液相黏滯系數(shù)由0.075減小到0.012 5時(shí)，空化泡最大壓強(qiáng)由0.007 18 mu·lu-1·tu-2增加到0.010 03 mu·lu-1·tu-2，增加了39%，最大微射流流速由0.255 4 lu·tu-1增加到0.319 9 lu·tu-1，增加了25%，而最大潰滅時(shí)間則由721 tu縮短到680 tu。LBM數(shù)值模擬結(jié)果表明降低液體黏滯系數(shù)會在空化時(shí)產(chǎn)生更大的微射流流速和潰滅壓強(qiáng)，與Popinet等[27]的試驗(yàn)所得結(jié)論一致。

表2 液相黏滯系數(shù)對空化泡潰滅影響工況

圖5 液相黏滯系數(shù)對空化泡潰滅的影響

4 結(jié) 論

a. LBM可以有效模擬近壁區(qū)空化過程，模擬結(jié)果可以有效觀察到Bjerknes力的產(chǎn)生，并伴隨著微射流的形成。

b. 由于汽相質(zhì)量較小，具有較小的慣性力。相同液相黏滯系數(shù)條件下，改變汽相黏滯系數(shù)，空化泡潰滅時(shí)產(chǎn)生的最大微射流流速、最大潰滅壓力和最大潰滅時(shí)間幾乎不變。

c. 保持汽相黏滯系數(shù)不變，空化泡最大潰滅壓強(qiáng)、最大微射流流速均隨著液相黏滯系數(shù)的增大而減小，而最大潰滅時(shí)間隨著液相黏滯系數(shù)的增大而增大。