鄒建鋒，何小瀧，袁 浩，汪凱迪，周昔東

(1.中國(guó)電建集團(tuán)貴陽勘測(cè)設(shè)計(jì)研究院有限公司，貴陽 550081； 2.重慶交通大學(xué) 重慶西南水運(yùn)工程科學(xué)研究所，重慶 400074； 3.中國(guó)三峽建工(集團(tuán))有限公司，成都 610000； 4.重慶交通大學(xué) 河海學(xué)院，重慶 400074)

1 研究背景

水利樞紐中，高水頭和高流速往往伴隨著負(fù)壓的產(chǎn)生，導(dǎo)致空化空蝕破壞，影響水利樞紐運(yùn)行的安全，因此空化空蝕在水利研究中一直倍受重視。前人研究表明，空化泡外流體與泡內(nèi)氣體壓差在短時(shí)間內(nèi)迅速變化，導(dǎo)致空化泡潰滅，潰滅過程中產(chǎn)生的微射流、沖擊波、高壓、高溫是空蝕破壞的主要原因[1-3]。

受限于復(fù)雜的邊界和有限的觀測(cè)手段，基于理論分析和試驗(yàn)的空化泡潰滅過程研究多著眼于空化泡形態(tài)的演化，難以獲得空化泡及周邊流體流場(chǎng)的連續(xù)變化過程。與之相比，計(jì)算流體力學(xué)可以直觀地反映出空化泡潰滅過程中流場(chǎng)、壓力場(chǎng)和溫度場(chǎng)的連續(xù)變化過程，有助于對(duì)空化泡潰滅中豐富的力學(xué)機(jī)理展開系統(tǒng)研究。傳統(tǒng)宏觀空化泡潰滅模型主要基于Navier-Stokes(N-S)方程和Rayleigh-Plesset(R-P)方程，此類方法需要引入界面捕捉模型(Volume of fluid，Level-set)對(duì)氣液交界面進(jìn)行捕捉，通常需要較小的步長(zhǎng)和較密的網(wǎng)格尺度，往往會(huì)帶來極大的計(jì)算量[4-5]。近30 a來，格子玻爾茲曼方法(Lattice Boltzmann Method，LBM)迅速發(fā)展為多相流研究的重要方法之一，該方法基于玻爾茲曼動(dòng)理學(xué)理論，求解介觀粒子分布函數(shù)，統(tǒng)計(jì)獲得宏觀尺度的密度場(chǎng)與流場(chǎng)。與傳統(tǒng)宏觀數(shù)值模擬方法相比，LBM無需求解宏觀方程中的對(duì)流項(xiàng)和泊松方程，適用于復(fù)雜邊界，且具有高并行性，已被用于模擬沸騰、蒸發(fā)、空化等復(fù)雜多相流現(xiàn)象[6-8]。

基于LBM的空化模型主要依托于偽勢(shì)模型，該類模型利用粒子間的作用勢(shì)，自動(dòng)形成氣液交界面[8]。Sukop等[9]首先模擬了各向異性條件下空化泡的生長(zhǎng)過程，結(jié)果與R-P理論解一致，但其液氣密度比僅為6.5。Chen等[10]引入精確差分法(exact-difference-method，EDM)外力格式后，將原始偽勢(shì)空化模型密度比提高到60。Shan等[11]引入多松弛(multi-relaxation-time，MRT)碰撞算子，進(jìn)一步將LBM偽勢(shì)空化模型密度比提高到750，液氣密度比接近于真實(shí)水與水蒸氣密度比。隨著熱力學(xué)方程被引入，LBM空化模型對(duì)空化現(xiàn)象的研究進(jìn)一步拓展到空化潰滅時(shí)溫度演化規(guī)律的研究[12-13]。

近壁區(qū)空化泡潰滅過程中，常伴隨著氣液交界面的劇烈變形，局部氣液交界面具有極大的曲率，在這些地方，表面張力往往對(duì)空化泡演化形態(tài)、潰滅強(qiáng)度、微射流分布具有重要的影響。而目前基于LBM偽勢(shì)模型的空化泡潰滅過程的模擬極少考慮表面張力對(duì)空化泡演化的影響。因此本研究引入表面張力調(diào)節(jié)外力項(xiàng)，建立可調(diào)節(jié)表面張力的LBM偽勢(shì)空化模型，分析不同表面張力對(duì)空化泡潰滅過程的影響，揭示表面張力對(duì)空化泡微射流流速分布、最大壓力變化、空化泡形態(tài)演化的影響。

2 數(shù)值模擬方法

2.1 LBM偽勢(shì)模型

本研究中采用多松弛系數(shù)碰撞算子、Li修正外力格式[14]、Carnahan-Starling(C-S)狀態(tài)方程等手段獲得大密度比、大黏滯系數(shù)比條件下，數(shù)值穩(wěn)定性更好的LBM偽勢(shì)模型，其對(duì)應(yīng)粒子分布方程為[14]

(1)

(2)

式中：ωi為平衡態(tài)函數(shù)對(duì)應(yīng)權(quán)重；ρ為密度；ueq為平衡態(tài)速度；cs為格子聲速。對(duì)于D2Q9格子模型，ω0=4/9，ω1～4=1/9，ω5～8=1/36。離散速度ei為

(3)

式中：c為格子常數(shù)，c=Δx/Δt；松弛矩陣Λ由不同松弛參數(shù)構(gòu)成，即

(4)

(5)

宏觀密度ρ和流場(chǎng)實(shí)際速度u可通過以下計(jì)算方式求得：

(6)

外力項(xiàng)Si采用Li等[14]提出的改進(jìn)外力格式，滿足大密度比條件下熱力學(xué)一致性，即

meq=Mfeq=

(7)

式中：ε為調(diào)節(jié)熱力學(xué)一致性的參數(shù)；ψ為粒子間的偽勢(shì)；τe為能量松馳系數(shù)；F為作用在粒子上的合力；Fm為粒子間的相互作用力，可以通過偽勢(shì)強(qiáng)度進(jìn)行求解，即

(8)

式中：G為粒子間的相互作用強(qiáng)度；w為不同方向的權(quán)重，w0=0，w1～4=1/3，w5～8=1/12[15]。引入非理想流體狀態(tài)方程后，可獲得求解偽勢(shì)ψ[16]的表達(dá)式為

(9)

式中P為壓力，本研究可用C-S狀態(tài)方程求解[16]，即

其中：

式中：Tc為臨界溫度；Pc為臨界壓力。

本模型采用線性插值獲得不同密度條件下的松弛系數(shù)，以獲得不同黏滯系數(shù)比[8]，即

(11)

式中：τυ(x)為混合流體黏滯系數(shù)；τg為氣相黏滯系數(shù)；τl為液相黏滯系數(shù)；ρg為初始?xì)庀嗝芏?；ρl為初始液相密度。

為探究表面張力對(duì)空化泡潰滅過程的影響，本模型引入源項(xiàng)C，獲得可單獨(dú)調(diào)節(jié)表面張力的偽勢(shì)模型[17]，即

(12)

式中Q可以通過下列方程求解，即

(13)

式中κ為表面張力調(diào)節(jié)參數(shù)。

2.2 邊界條件及初始參數(shù)

近壁區(qū)空化泡潰滅計(jì)算示意圖如圖1所示，計(jì)算區(qū)域?yàn)?00 lu×200 lu矩形區(qū)域，r0為空化泡初始半徑，d為空化泡中心與壁面之間的間距，λ=d/r0為空化泡與壁面之間的無量綱距離。計(jì)算域頂部邊界為Zou-He壓力邊界[8]，底部邊界為無滑移邊界，左右邊界為非平衡外推邊界。

圖1 近壁區(qū)空化泡潰滅計(jì)算域示意圖Fig.1 Computational domain for cavitation bubble collapse in near-wall region

為促使模擬過程盡快達(dá)到平衡態(tài)，不影響空化泡演化過程，空化泡及周邊計(jì)算域流場(chǎng)按照式(14)進(jìn)行初始化。

(14)

式中：(x0，y0)為液滴中心坐標(biāo)；w為界面厚度，本研究中選取界面厚度w=5 lu。

3 模型驗(yàn)證

本研究采用Laplace定律和無限域空化泡潰滅過程驗(yàn)證空化模型的可靠性。

3.1 Laplace定律

首先利用Laplace定律驗(yàn)證并獲得不同κ值條件下的表面張力。Laplace定義為，當(dāng)空化泡或液滴達(dá)到平衡狀態(tài)時(shí)，空化泡內(nèi)外壓差與半徑的倒數(shù)成正比，即其斜率即為該條件下的表面張力，Δp=σ/r。式中Δp為空化泡內(nèi)外壓差，σ為表面張力，req為空化泡平衡態(tài)條件下半徑。

本文選取5種初始半徑用于Laplace定律的研究，分別為r0=35、40、45、50、55 lu。不同值條件下空化泡內(nèi)外壓差Δp和平衡態(tài)半徑r之間關(guān)系如圖2所示。本模型獲得的空化泡內(nèi)外壓差與其平衡時(shí)刻的半徑的倒數(shù)具有較好的線性關(guān)系。隨著κ值的減小，表面張力由0.001 9增大到0.018 3，表明當(dāng)前模型可以有效調(diào)節(jié)表面張力大小。

圖2 不同κ值下空化泡內(nèi)外壓差和平衡態(tài)半徑的關(guān)系Fig.2 Relations of pressure difference inside and outside of the bubble against bubble radius under different κ values

3.2 無限域空化泡潰滅過程

為進(jìn)一步驗(yàn)證模型的準(zhǔn)確性，本研究還對(duì)無限域內(nèi)空化泡潰滅過程進(jìn)行模擬，并將LBM模擬結(jié)果與Rayleigh-Plesset方程數(shù)值解進(jìn)行對(duì)比驗(yàn)證?？紤]表面張力，對(duì)應(yīng)的Rayleigh-plesset空化泡演化方程為[18]

(15)

模擬中設(shè)置計(jì)算域四周均為Zou-He壓力邊界，初始液相壓力p∞=0.006 79 mu/(lu·tu2)，空化泡內(nèi)初始?jí)毫v=7×10-7mu/(lu·tu2)，空化泡半徑r=50 lu，表面張力σ=0.010 2 mu/tu2，其他參數(shù)與2.2節(jié)一致。無限域內(nèi)空化泡潰滅中半徑演化過程的LBM模擬結(jié)果與Rayleigh-Plesset方程對(duì)比如圖3所示，圖中r*為無量綱空化泡半徑，t*為無量綱空化泡潰滅時(shí)間，LBM模擬結(jié)果與R-P方程理論解一致，盡管在空化泡潰滅階段初期存在部分誤差，其可能是由于LBM模擬中對(duì)初始化流場(chǎng)調(diào)整而引入的。

圖3 R-P方程理論解與LBM模擬結(jié)果對(duì)比Fig.3 Comparison between the result of Rayleigh-Plesset equation and LBM simulation

4 表面張力對(duì)空化過程的影響

本節(jié)選取空化泡內(nèi)外初始?jí)翰瞀=0.003 8、0.006 8、0.010 2 mu/(lu·tu2)三種壓力，無量綱距離λ=1.6條件下，不同表面張力系數(shù)對(duì)空化泡潰滅過程中潰滅形態(tài)、流場(chǎng)、最大微射流速度和最大潰滅壓力的影響。

圖4展示了Δp=0.006 8 mu/(lu·tu2)不同表面張力條件下空化泡潰滅最終形態(tài)及流場(chǎng)。對(duì)比不同表面張力條件下壓力分布，隨著表面張力減小，空化泡上方高壓區(qū)面積增大，同時(shí)空化泡與壁面之間的低壓區(qū)面積減小。在較小表面張力條件下，空化泡更易發(fā)生變形，界面曲率半徑更小，空化泡潰滅最終時(shí)刻形變更大，導(dǎo)致潰滅形成的微射流流速更為集中，這與Zhang等[19]的研究結(jié)論一致。

圖4 不同表面張力條件下空化泡潰滅最終時(shí)刻密度、速度和壓力分布密度Fig.4 Distribution of density，velocity，and pressure at the final collapse stage under different surface tensions

圖5展示了Δp=0.006 8 mu/(lu·tu2)時(shí)不同表面張力條件下空化泡潰滅過程中最大流速和最大壓力演化過程。潰滅過程中更大的表面張力條件下空化泡積蓄了更多的表面能，并在潰滅最后階段迅速釋放，導(dǎo)致了更大的空化泡潰滅速度。隨著表面張力增加，空化泡潰滅時(shí)間增加，同時(shí)空化泡潰滅最大流速也相應(yīng)增加。當(dāng)表面張力由0.001 9 mu/tu2增加到0.018 3 mu/tu2時(shí)，空化泡潰滅流速由0.404 lu/tu增加到0.512 lu/tu，速度增加了26.8%。由圖5(b)可知，在空化泡潰滅最后階段，由于空化泡內(nèi)水蒸氣在短時(shí)間內(nèi)發(fā)生相變，潰滅壓力迅速增加。且更大的表面張力阻礙了空化泡的形變和潰滅，導(dǎo)致空化泡潰滅壓力減小。在Δp=0.006 8 mu/(lu·tu2)壓力條件下，當(dāng)表面張力由0.001 9 mu/tu2增加到0.018 3 mu/tu2時(shí)，空化泡潰滅過程中最大壓力由0.049 mu/(lu·tu2)增加到0.064 mu/(lu·tu2)。

圖5 不同表面張力條件下空化泡潰滅過程中最大速度和最大壓力演化過程Fig.5 Evolution of maximum velocity and maximum pressure in the collapse process under different surface tensions

不同表面張力和初始空化泡內(nèi)外壓差條件下空化泡最大微射流流速、壓力和潰滅時(shí)間變化規(guī)律如圖6所示。隨著初始內(nèi)外壓差增大，不同表面張力條件下空化泡潰滅產(chǎn)生的最大微射流流速、最大潰滅壓力均隨之增大，但潰滅時(shí)間則隨之減小。初始空化泡內(nèi)外壓差為0.003 8 mu/(lu·tu2)時(shí)，當(dāng)表面張力由0.001 9 mu/tu2增加到0.018 3 mu/tu2，最大微射流流速增加了55%，最大潰滅壓力增加了74%。初始空化泡內(nèi)外壓差增加到0.010 2 mu/(lu·tu2)時(shí)，空化泡潰滅最大微射流流速僅增加了22%，潰滅壓力僅增加了6%，最大微射流流速和最大潰滅壓力增加幅度均隨著空化泡初始內(nèi)外壓差的增加而減小，說明空化泡初始內(nèi)外壓差的增加會(huì)減小表面張力對(duì)潰滅過程的影響。

圖6 不同壓力條件下空化泡潰滅過程中各參數(shù)隨表面張力的變化規(guī)律Fig.6 Vaiation law of parameters with surface tension under different pressures during the process of civitation bubble collapse

而在相同壓力條件下，空化泡潰滅速度、潰滅壓力和潰滅時(shí)間均隨著表面張力的增加而增加。根據(jù)Bjerknes力的定義[20]，有

(16)

空化泡潰滅時(shí)間隨著表面張力的減小而減小，而時(shí)間質(zhì)量變化率則隨之增大，導(dǎo)致Bjerknes力增大，說明較小的表面張力會(huì)促使空化泡朝向壁面發(fā)生潰滅。

依托于試驗(yàn)研究和宏觀數(shù)值模擬方法，研究者們針對(duì)表面張力對(duì)空化的影響開展了系統(tǒng)性研究[21-22]。LBM模擬結(jié)果表明表面張力對(duì)空化泡潰滅強(qiáng)度具有重要影響，其Bjerknes力隨表面張力減小而增大，但潰滅強(qiáng)度卻隨之減小，模擬結(jié)果定性上與前人試驗(yàn)結(jié)果和數(shù)值模擬結(jié)果一致。摻氣減蝕是有效減小空化空蝕的重要手段，研究表明摻氣濃度將減小流體表面張力，進(jìn)而減小空蝕強(qiáng)度。

5 結(jié) 論

本文基于可調(diào)節(jié)表面張力的LBM偽勢(shì)空化模型，研究了表面張力對(duì)近壁區(qū)空化泡潰滅的影響，得到結(jié)論如下：

(1)表面張力直接影響空化泡附近區(qū)域壓力演化，并對(duì)空化強(qiáng)度產(chǎn)生影響。較小的表面張力在增加空化泡上部高壓范圍，同時(shí)減小壁面與空化泡之間的低壓區(qū)范圍。空化泡水汽交界面抗形變能力隨之減弱，導(dǎo)致微射流范圍更集中。

(2)空化泡潰滅時(shí)產(chǎn)生最大微射流和潰滅壓力均隨表面張力增加而增加。較大的表面張力條件下空化泡演化過程中積蓄了更多的表面能，并在潰滅階段釋放，加劇潰滅強(qiáng)度。

(3)空化泡潰滅時(shí)間隨表面張力減小而減短，導(dǎo)致Bjerknes力增強(qiáng)，加劇了空化泡朝壁面潰滅的趨勢(shì)。