湯冬冬，呂續(xù)艦，運(yùn)洪祿，陳少松

(南京理工大學(xué) 能源與動(dòng)力工程學(xué)院, 南京 210094)

超空泡現(xiàn)象因其良好的減阻特性得到了廣泛的關(guān)注。在超空泡流動(dòng)中，流經(jīng)運(yùn)動(dòng)體的水達(dá)到飽和蒸汽壓后變成粘度更小的水蒸汽，導(dǎo)致運(yùn)動(dòng)體受到的摩擦阻力急劇下降。在這一現(xiàn)象中，空化數(shù)是一個(gè)關(guān)鍵參數(shù)，其表達(dá)式為σ=2(p∞-pv)/ρV2。其中，p∞和pv分別是環(huán)境壓力和飽和蒸汽壓力，V是來流的速度，ρ是來流水的密度。不同的空化數(shù)，對(duì)應(yīng)著不同的空泡形態(tài)。當(dāng)形成的空泡完全包裹運(yùn)動(dòng)體，空泡將運(yùn)動(dòng)體與水隔離，能夠降低阻力，提高航速的維持性和航行距離。這種流動(dòng)現(xiàn)象稱為自然超空化。

針對(duì)空化特性的研究，國內(nèi)外許多學(xué)者通過數(shù)值模擬或者模型實(shí)驗(yàn)開展了大量卓有成效的研究工作。金大橋等[1]通過數(shù)值模擬研究了空化數(shù)對(duì)水下射彈空泡閉合部位及阻力系數(shù)的影響。王復(fù)峰等[2]采用實(shí)驗(yàn)與數(shù)值模擬相結(jié)合方法，得出了在不同空化數(shù)和相同的通氣速率下，繞回轉(zhuǎn)體通氣空化三相流動(dòng)的空泡形態(tài)不同于兩相流動(dòng)的空泡形態(tài)。陳瑋琪[3]通過模擬細(xì)長體的空泡流動(dòng)現(xiàn)象，建立了自由面影響下的氣團(tuán)和空泡相互耦合動(dòng)力學(xué)模型，提出了非定常流場中帶氣體泄漏現(xiàn)象的含氣空泡壓力估算方法。周景軍[4]等發(fā)現(xiàn)當(dāng)速度達(dá)到50 m/s以上時(shí)，重力效應(yīng)對(duì)空泡生成速度影響可以忽略。張曉東[5]求解Paryshev空泡模型的基礎(chǔ)上，計(jì)算了重力場中垂直空泡的長度變化。張紀(jì)華等[6]研究發(fā)現(xiàn)，在保證生成穩(wěn)定超空泡前提下，空泡生成相對(duì)速度和相對(duì)時(shí)間隨通氣流量的增大出現(xiàn)相反的變化趨勢(shì)。白濤等[7]基于水洞實(shí)驗(yàn)，采用強(qiáng)跟蹤濾波算法，得到了超空泡形態(tài)的動(dòng)態(tài)估計(jì)算法。Fu等[8]利用FLUENT軟件模擬超空泡減阻，并發(fā)現(xiàn)運(yùn)動(dòng)體形狀、長徑比和空化數(shù)等的最佳組合可以達(dá)到最優(yōu)的減阻效果。Ahn等[9]通過研究空化器楔形角對(duì)空泡形狀的影響時(shí)發(fā)現(xiàn)，楔形角越大形成的空泡越長、直徑越大。烏克蘭的Savchenko[10]較早研究了不同空化數(shù)下的空泡長度和大小，通過實(shí)驗(yàn)研究，擬合出空泡長度、直徑、阻力系數(shù)與空化數(shù)的關(guān)系。Jiang等[11]計(jì)算了二維自然超空泡運(yùn)動(dòng)體阻力及空泡尺寸等參數(shù)，并與Savchenko的實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。Kadivar等[12]對(duì)全尺寸空化器進(jìn)行了數(shù)值模擬，分析了空泡直徑和長度等參數(shù)，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好。

空化流動(dòng)受到的阻力主要由流動(dòng)過程中不可逆因素引起，具體表現(xiàn)為湍流引起的耗散和溫差產(chǎn)生的耗散，即熵產(chǎn)。熱力學(xué)理論常用熵產(chǎn)來評(píng)估能量傳遞過程的不可逆程度，許多學(xué)者亦將其應(yīng)用于流動(dòng)特性研究。王威等[13]利用熵產(chǎn)最小理論優(yōu)化二維跨音速翼型設(shè)計(jì)，使得翼型在來流馬赫數(shù)為0.73、攻角為2.54°時(shí)升阻比達(dá)到最大；他們還利用熵產(chǎn)分析得出流場熵產(chǎn)和翼型阻力系數(shù)存在線性關(guān)系[14]。Li等[15]更早地使用同樣的方法設(shè)計(jì)翼型，由二維拓展到三維模擬中，并對(duì)比通過流域熵產(chǎn)分析得到的阻力系數(shù)與直接沿物體表面積分得到阻力系數(shù)的偏差，發(fā)現(xiàn)該方法在三維流動(dòng)計(jì)算中的應(yīng)用有待提高。吉宇等[16]利用熵產(chǎn)分析蒸汽射流直接接觸冷凝過程，將兩相流與傳熱傳質(zhì)結(jié)合。Jin等[17]利用熵產(chǎn)分布分析肋片熱量的傳遞以及如何減小肋片阻力。對(duì)于單相流動(dòng)，研究的結(jié)果表明阻力與熵產(chǎn)存在一定聯(lián)系；對(duì)于多相流動(dòng)的熵產(chǎn)，研究其分布特性較多，而分析熵產(chǎn)與運(yùn)動(dòng)體阻力之間關(guān)系的研究非常少。本文擬在常規(guī)超空泡特性研究的基礎(chǔ)上，探索兩相流動(dòng)中運(yùn)動(dòng)體干擾流場引起的熵產(chǎn)與運(yùn)動(dòng)體阻力之間的聯(lián)系，為多相流動(dòng)中阻力的評(píng)估和優(yōu)化探索一種可能的新思路。

1 數(shù)值模型

1.1 控制方程

本文基于均質(zhì)多相流理論，采用基于有限體積法的商業(yè)軟件ANSYS FLUENT開展二維運(yùn)動(dòng)體自然超空化模擬。采用多相流模型，假設(shè)流體由水、水蒸汽和非凝結(jié)氣體組成，水的汽化速度與水蒸汽的冷凝速度相同。在混合模型中，連續(xù)性方程和動(dòng)量矩方程如下：

(1)

(2)

其中，u、p、g和Fi分別為速度矢量、壓力、重力加速度和體積力。ρm表示為混合相的密度，有ρm=(1-α)ρl+αρv；μm表示為混合相的動(dòng)力黏度，有μm=(1-α)μl+αμv；ρl、ρv、μl、μv和α分別為液體密度、蒸汽密度、液體動(dòng)力粘度、蒸汽動(dòng)力粘度和蒸汽體積率。

空化模型采用Schner-Sauer模型，汽化的體積率滿足如下方程：

(3)

其中，Re和Rc分別為蒸汽蒸發(fā)速率和冷凝速率。

當(dāng)pv>p時(shí)：

(4)

當(dāng)pv≤p時(shí)：

(5)

采用的能量守恒方程如下：

(6)

從熱力學(xué)角度看，空化過程是一個(gè)不可逆熵產(chǎn)過程，不可逆因素主要由湍流和溫差引起。湍流因素包括物質(zhì)的粘性耗散和流動(dòng)中的壓力差，而溫差與溫度場分布有關(guān)。兩相流中，首先計(jì)算單個(gè)網(wǎng)格熵產(chǎn)S?gen，然后在流域內(nèi)積分即可得到總熵產(chǎn)Sgen。單個(gè)網(wǎng)格熵產(chǎn)計(jì)算如下[15]：

(7)

對(duì)于二維情況：

(8)

(9)

其中，εq為湍流耗散率；μeff,q為第q相的有效粘性系數(shù)，包括流體本身的粘性系數(shù)和流體脈動(dòng)引起的湍流粘性系數(shù)；ρq為第q相的密度。

1.2 計(jì)算方法和條件

考慮一個(gè)二維運(yùn)動(dòng)體，其完整尺寸為30 mm×7 mm，采用軸對(duì)稱模型來模擬當(dāng)前空化現(xiàn)象，對(duì)應(yīng)的計(jì)算流域須足夠大以減小遠(yuǎn)場影響。取計(jì)算域?yàn)?1 000 mm×100 mm，來流邊界距離運(yùn)動(dòng)體前緣270 mm(圖1)。計(jì)算邊界條件設(shè)置如下：來流邊界為速度入口；出口邊界設(shè)置為壓力出口，壓力值取標(biāo)準(zhǔn)大氣壓 101 325 Pa；流域上側(cè)面為與來流一致的速度入口，下側(cè)面為無滑移壁面。采用ANSYS ICEM進(jìn)行建模和網(wǎng)格劃分，為了提高計(jì)算精度和更好地捕捉到流場信息，對(duì)運(yùn)動(dòng)體周邊敏感區(qū)域進(jìn)行局部加密處理，如圖2所示。

圖1 二維運(yùn)動(dòng)體模型及計(jì)算域

圖2 運(yùn)動(dòng)體局部流域網(wǎng)格

湍流模型采用k-ε標(biāo)準(zhǔn)模型，對(duì)應(yīng)的壓力和速度的耦合采用SIMPLE算法。SIMPLE算法在滿足Navier-Stokes方程的同時(shí)，基于質(zhì)量守恒方程來確定速度和壓強(qiáng)的關(guān)系。壓力項(xiàng)采用PRESTO方法進(jìn)行處理，該方法采用離散連續(xù)平衡法去計(jì)算交錯(cuò)網(wǎng)格上的壓力，常用于兩相流動(dòng)中?；旌厦芏?、動(dòng)量方程、湍動(dòng)能和湍流耗散等均采用二階迎風(fēng)格式離散處理。

1.3 計(jì)算模型的確認(rèn)和驗(yàn)證

在開展計(jì)算研究之前，須對(duì)計(jì)算模型進(jìn)行確認(rèn)，主要考慮網(wǎng)格數(shù)獨(dú)立性和運(yùn)動(dòng)體壁面第一層網(wǎng)格y+獨(dú)立性分析。數(shù)值計(jì)算所耗費(fèi)的時(shí)間與網(wǎng)格數(shù)量直接相關(guān)，網(wǎng)格數(shù)量越多，計(jì)算所需時(shí)間越長。確定合適的網(wǎng)格數(shù)和網(wǎng)格分布形式，在滿足所需計(jì)算精度的同時(shí)，盡量節(jié)省計(jì)算機(jī)時(shí)。圖3所示為速度40 m/s時(shí)阻力F和熵產(chǎn)Sgen隨著網(wǎng)格數(shù)N的變化曲線。網(wǎng)格數(shù)從3.3萬依次增加到6.7萬和10.6萬，阻力和熵產(chǎn)增加顯著，變化率超過35%；當(dāng)網(wǎng)格數(shù)從10.6萬增加到21萬后，阻力變化約1%，熵產(chǎn)變化約3%，所以計(jì)算網(wǎng)格數(shù)取10.6萬(圖3所示的收斂點(diǎn))比較合理。

圖3 阻力和熵產(chǎn)隨網(wǎng)格數(shù)變化曲線

運(yùn)動(dòng)體近壁面存在較大的法向速度梯度，對(duì)計(jì)算阻力存在影響。y+反映了數(shù)值計(jì)算對(duì)運(yùn)動(dòng)體近壁面信息捕捉的準(zhǔn)確性，與近壁面第一層網(wǎng)格高度、流體密度、壁面摩擦速度和粘性系數(shù)有關(guān)。具體定義如下：

(10)

式(10)中：Δy表示為網(wǎng)格第一層的高度；uτ表示為壁面摩擦速度；μ表示為流體粘性系數(shù)。

運(yùn)動(dòng)體近壁面不同位置y+值不同。便于比較近壁面情況，統(tǒng)一取運(yùn)動(dòng)體壁面y+的平均值作為比較標(biāo)準(zhǔn)。圖4展示了運(yùn)動(dòng)體在40 m/s(σ=0.123 9)下y+值的獨(dú)立性計(jì)算結(jié)果，當(dāng)y+從95依次減小到83和39.6時(shí)，阻力和熵產(chǎn)變化顯著，分別增加7%～25%和22%～23%；當(dāng)y+從39.6減小到5.47時(shí)，阻力變化約1.1%，熵產(chǎn)變化約2.7%，所以壁面y+取39.6(收斂點(diǎn))比較合理。

通過不同速度下典型工況計(jì)算發(fā)現(xiàn)，本文關(guān)注的運(yùn)動(dòng)體來流速度在非空化階段(V≤12 m/s)和空化階段(12

2 計(jì)算結(jié)果和分析

2.1 計(jì)算模型驗(yàn)證

超空化運(yùn)動(dòng)體的阻力與其周邊流動(dòng)速度場、壓力場和粘性系數(shù)等密切相關(guān)，能夠間接反映超空化流動(dòng)特性。圖5所示為不同空化數(shù)下運(yùn)動(dòng)體的阻力系數(shù)Cd(Cd=2F/ρV2S，S表示為運(yùn)動(dòng)體最大橫截面積)變化曲線。隨著空化數(shù)的增加，運(yùn)動(dòng)體所受阻力系數(shù)增加，兩者幾乎成線性關(guān)系。與實(shí)驗(yàn)結(jié)果相比，本文數(shù)值計(jì)算結(jié)果誤差為1.2%～2.2%，且總體比實(shí)驗(yàn)結(jié)果偏高。相較于Jiang等[11]的計(jì)算值(總體誤差3%～7%)，本文的計(jì)算結(jié)果更接近實(shí)驗(yàn)值。

圖4 阻力和熵產(chǎn)隨壁面y+變化曲線

自然空泡的演變是一個(gè)不穩(wěn)定過程，選取較小的時(shí)間步長(10-4s)，通過分析空泡動(dòng)態(tài)演變過程，給出空泡沿運(yùn)動(dòng)體的初生、發(fā)展變化規(guī)律。圖6展示了V=50 m/s(σ=0.079 3)的空泡初生和發(fā)展變化情況。由于低壓區(qū)存在，運(yùn)動(dòng)體的頭部兩側(cè)以及尾部后端會(huì)首先產(chǎn)生空泡，產(chǎn)生空泡的范圍和大小隨時(shí)間推移向后發(fā)展，漸漸地包裹住整個(gè)運(yùn)動(dòng)體，最終形成將運(yùn)動(dòng)體與水隔絕的超空泡。

圖6 空泡形態(tài)隨時(shí)間變化情況

為了更好地描述空泡的最終形態(tài)、展示不同空化數(shù)下空泡的尺寸規(guī)律，采用空泡相對(duì)長度Lc/Dn和相對(duì)直徑Dc/Dn來描述空泡形態(tài)。其中，Lc為空泡的最大長度，Dc是為空泡的最大直徑，Dn是為空化器的直徑。Savchenko通過實(shí)驗(yàn)給出了空泡直徑和長度的經(jīng)驗(yàn)公式如下：

(11)

(12)

圖7所示為空泡相對(duì)長度和相對(duì)直徑隨空化數(shù)的變化曲線，可以發(fā)現(xiàn)空泡的相對(duì)長度和相對(duì)直徑隨著空化數(shù)增大而減小。與實(shí)驗(yàn)結(jié)果比較，本文計(jì)算的空泡相對(duì)長度和相對(duì)直徑與Jiang等[11]的計(jì)算結(jié)果誤差相當(dāng)，基本在20%左右，但本文計(jì)算結(jié)果總體趨勢(shì)更接近實(shí)驗(yàn)值，特別是在低空化數(shù)范圍內(nèi)(σ≤0.065 6)?？梢哉J(rèn)為，本文的計(jì)算模型是可行的。

圖7 空泡形態(tài)隨空化數(shù)變化曲線

2.2 熵產(chǎn)和阻力分析

采用ANSYS FLUENT研究空化問題時(shí)，為保證計(jì)算的穩(wěn)定性，一般考慮定溫環(huán)境下的相變情況，即可假設(shè)式(7)中?T/?x=0，?T/?y=0。在考慮流場熵產(chǎn)時(shí)，本文取常溫下(25 ℃)的水及其飽和蒸汽的物理參數(shù)為參考數(shù)據(jù)。

圖8所示為二維運(yùn)動(dòng)體阻力和流場熵產(chǎn)隨空化數(shù)的變化曲線?？梢钥闯觯枇挽禺a(chǎn)均隨空化數(shù)的增大而減小，且變化趨勢(shì)基本一致，即在超空化區(qū)域都隨空化數(shù)的增加迅速降低、在局部空化區(qū)域都隨空化數(shù)增加而有明顯降低但降低趨勢(shì)大幅減緩，而在非空化區(qū)域二者變化都非常小。

圖8 阻力和流場熵產(chǎn)隨空化數(shù)的變化曲線

運(yùn)動(dòng)體阻力的計(jì)算通常采用物面壓力積分、剪切應(yīng)力積分、邊界積分和尾跡動(dòng)量損失積分等方法進(jìn)行考慮。對(duì)于不可壓縮流動(dòng)，一般認(rèn)為阻力與流場熵產(chǎn)存在線性關(guān)系[14-15]，而從熱力學(xué)第二定律能量不可逆的角度可以對(duì)流場的熵產(chǎn)進(jìn)行計(jì)算，那么通過熵產(chǎn)分析求解運(yùn)動(dòng)體的阻力不失為一個(gè)可行的思路。鑒于此，圖9展示了運(yùn)動(dòng)體受到的阻力和流場熵產(chǎn)之間的關(guān)系曲線。根據(jù)是否空化以及空化程度，將該關(guān)系曲線分為非空化區(qū)(AB段)、局部空化區(qū)(BC段)和超空化區(qū)(CD段)。其中點(diǎn)B和C分別表示空化臨界點(diǎn)和超空化臨界點(diǎn)。在上述三個(gè)分區(qū)內(nèi)，熵產(chǎn)和阻力均分別呈線性關(guān)系，可統(tǒng)一表示為

F=kSgen+b

(13)

式(13)中：k和b分別為斜率系數(shù)和截距系數(shù)，對(duì)于非空化區(qū)、局部空化區(qū)和超空化區(qū)，可分別取k1=19.239、b1=0.441；k2=7.278、b2=2.600；k3=3.457、b3=9.113。顯然，上述斜率之間存在如下關(guān)系：非空化區(qū)>局部空化區(qū)>超空化區(qū)。

圖9 阻力隨流場熵產(chǎn)的變化曲線

圖10所示為典型工況下空泡的形態(tài)與流場熵產(chǎn)分布情況。流動(dòng)處于非空化狀態(tài)時(shí)(圖10(a))，熵產(chǎn)主要分布在運(yùn)動(dòng)體頭部；在局部空化時(shí)(圖10(b))，熵產(chǎn)主要分布在運(yùn)動(dòng)體的頭部、汽液交界面處和運(yùn)動(dòng)體尾部；而在超空化情況下(圖10(c))，流場熵產(chǎn)沿著汽液交界面分布，主要集中于運(yùn)動(dòng)體頭部附近，其次為空泡尾部附近。式(7)表明，影響流場熵產(chǎn)的主要因素是速度梯度和有效粘性系數(shù)，如圖11(a)～圖11(d)所示，其中u，v分別表示x，y方向的速度。而湍流情況下湍流粘性系數(shù)μtur是有效粘性系數(shù)μeff的主要組成部分，其分布情況如圖11(e)所示，與熵產(chǎn)分布趨勢(shì)基本一致。運(yùn)動(dòng)體前緣和運(yùn)動(dòng)體尾部的熵產(chǎn)主要由水流折轉(zhuǎn)的速度梯度產(chǎn)生的阻力造成；空泡內(nèi)的熵產(chǎn)幾乎為零，其原因是超空泡內(nèi)水蒸汽的速度梯度小、水蒸汽的粘性系數(shù)遠(yuǎn)低于水且其湍流脈動(dòng)極低；超空泡尾部流場由于尾流汽液交界面上的流動(dòng)不穩(wěn)定性導(dǎo)致較大的湍流粘性系數(shù)和速度梯度。

圖10 典型熵產(chǎn)(左)與(右)空泡形態(tài)分布

圖11 典型速度梯度與湍流粘性系數(shù)分布(V=40 m/s)

圖12展示了典型工況下沿運(yùn)動(dòng)體表面的熵產(chǎn)分布曲線，可以發(fā)現(xiàn)熵產(chǎn)分布峰值集中在運(yùn)動(dòng)體頭部，進(jìn)一步說明了運(yùn)動(dòng)體的能耗主要發(fā)生在前緣處。表1表示熵產(chǎn)沿運(yùn)動(dòng)體長度面積分值。從表1中直觀看出，運(yùn)動(dòng)體前1/10段內(nèi)集中約一半的熵產(chǎn)。

圖12 典型工況下熵產(chǎn)沿運(yùn)動(dòng)體表面分布

V/(m·s-1)σA1(W/K)B1(W/K)A1/B1(%)40.000.1239157.41325.7748.3245.000.0979196.33409.0348.0050.000.0793237.10445.8953.1755.000.0656258.63472.8954.6960.000.0551403.24811.3449.70

注：A1為熵產(chǎn)沿前1/10段的面積分；B1熵產(chǎn)沿整段的面積分。

3 結(jié)論

1) 本文的計(jì)算模型較為合理，得到的阻力系數(shù)和空泡尺寸與文獻(xiàn)實(shí)驗(yàn)值和計(jì)算值吻合較好；運(yùn)動(dòng)體所受阻力隨速度增大而增大，阻力系數(shù)與空化數(shù)幾乎成線性變化，空泡尺寸與空化數(shù)密切相關(guān)；超空泡初生過程中，空泡首先產(chǎn)生于頭部與尾部后端，并逐漸向后發(fā)展包裹住運(yùn)動(dòng)體。

2) 當(dāng)運(yùn)動(dòng)體以不同速度運(yùn)動(dòng)，在非空化、局部空化和超空化情況下，其阻力與流動(dòng)熵產(chǎn)存在分段線性關(guān)系。流場熵產(chǎn)主要與速度梯度和湍流脈動(dòng)引起的粘性系數(shù)密切相關(guān)。非空化時(shí)，熵產(chǎn)主要集中在運(yùn)動(dòng)體的前緣，少量位于運(yùn)動(dòng)體尾部；局部空化時(shí)，熵產(chǎn)主要分布在汽液交界面上和運(yùn)動(dòng)體尾部；超空化時(shí)，熵產(chǎn)沿著汽液交界面上分布，主要集中在運(yùn)動(dòng)體頭部和空泡尾部。