基于熵產分析的二維運動體自然超空化計算

湯冬冬，呂續艦，運洪祿，陳少松

(南京理工大學 能源與動力工程學院, 南京 210094)

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

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

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

1 數值模型

1.1 控制方程

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

(1)

(2)

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

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

(3)

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

當pv>p時：

(4)

當pv≤p時：

(5)

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

(6)

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

(7)

對于二維情況：

(8)

(9)

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

1.2 計算方法和條件

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

圖1 二維運動體模型及計算域

圖2 運動體局部流域網格

湍流模型采用k-ε標準模型，對應的壓力和速度的耦合采用SIMPLE算法。SIMPLE算法在滿足Navier-Stokes方程的同時，基于質量守恒方程來確定速度和壓強的關系。壓力項采用PRESTO方法進行處理，該方法采用離散連續平衡法去計算交錯網格上的壓力，常用于兩相流動中?；旌厦芏?、動量方程、湍動能和湍流耗散等均采用二階迎風格式離散處理。

1.3 計算模型的確認和驗證

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

圖3 阻力和熵產隨網格數變化曲線

運動體近壁面存在較大的法向速度梯度，對計算阻力存在影響。y+反映了數值計算對運動體近壁面信息捕捉的準確性，與近壁面第一層網格高度、流體密度、壁面摩擦速度和粘性系數有關。具體定義如下：

(10)

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

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

通過不同速度下典型工況計算發現，本文關注的運動體來流速度在非空化階段(V≤12 m/s)和空化階段(12

2 計算結果和分析

2.1 計算模型驗證

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

圖4 阻力和熵產隨壁面y+變化曲線

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

圖6 空泡形態隨時間變化情況

為了更好地描述空泡的最終形態、展示不同空化數下空泡的尺寸規律，采用空泡相對長度Lc/Dn和相對直徑Dc/Dn來描述空泡形態。其中，Lc為空泡的最大長度，Dc是為空泡的最大直徑，Dn是為空化器的直徑。Savchenko通過實驗給出了空泡直徑和長度的經驗公式如下：

(11)

(12)

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

圖7 空泡形態隨空化數變化曲線

2.2 熵產和阻力分析

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

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

圖8 阻力和流場熵產隨空化數的變化曲線

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

F=kSgen+b

(13)

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

圖9 阻力隨流場熵產的變化曲線

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

圖10 典型熵產(左)與(右)空泡形態分布

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

圖12展示了典型工況下沿運動體表面的熵產分布曲線，可以發現熵產分布峰值集中在運動體頭部，進一步說明了運動體的能耗主要發生在前緣處。表1表示熵產沿運動體長度面積分值。從表1中直觀看出，運動體前1/10段內集中約一半的熵產。

圖12 典型工況下熵產沿運動體表面分布

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為熵產沿前1/10段的面積分；B1熵產沿整段的面積分。

3 結論

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

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