袁曉龍，何小瀧,汪凱迪

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

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

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

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

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

1 數值模擬方法

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

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

(1)

(2)

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

(3)

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

(4)

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

(5)

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

(6)

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

(7)

為滿足熱力學一致性，本文中外力采用Li等[24]提出的偽勢模型外力項，其表達式為

(8)

(9)

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

(10)

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

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

(11)

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

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

2 模型驗證

2.1 熱力學一致性驗證

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

圖1 汽液黏滯系數比對熱力學一致性的影響

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

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

圖2 計算域示意圖

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

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

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

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

3 汽相和液相黏滯系數變化的影響

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

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

為研究液相黏滯系數變化對空化泡潰滅過程中最大微射流流速、最大壓強、最大潰滅時間的影響，選取液相松弛系數和對應的液相黏滯系數如表2所示。圖5展示了液相黏滯系數對空化泡潰滅過程中最大微射流流速、最大潰滅壓力和最大潰滅時間的影響，可以看出空化泡最大潰滅時間隨著液相黏滯系數的減小而減小，而空化泡潰滅最大微射流流速和最大潰滅壓強均隨著液相黏滯系數的減小而增大，當液相黏滯系數由0.075減小到0.012 5時，空化泡最大壓強由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%，而最大潰滅時間則由721 tu縮短到680 tu。LBM數值模擬結果表明降低液體黏滯系數會在空化時產生更大的微射流流速和潰滅壓強，與Popinet等[27]的試驗所得結論一致。

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

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

4 結 論

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

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

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