超空化流動數值模擬

0空化器特性仿真研究

超空氣泡法是一種高效、低可見水支持技術的方法。其原理是在水與航行體表面之間形成

穩定的氣層,達到減小水的粘性阻力，從而提高水下航行體的運行速度。為了能夠生成穩定

的、高效減阻的、形態可控的超空泡，位于航行體頭部的空化器起著決定性作用。空化器

形狀直接影響超空泡形態及水下航行體減阻效果，因此優化空化器形狀是獲得穩定高效減

阻超空泡的前提,其中圓盤空化器具有控制超空泡形態的特點。因此國內外的學者針對空

化器作了大量實驗研究即理論分析工作Savchenko利用控制空化器參數來獲得穩定的超空

泡,并總結出經大量實驗驗證過的經驗公式。Stinebring等對頭部為圓盤空化器進行了研

究，并給出了空泡長度及厚度與空化器參數關系的公式。國內何友聲等對系列頭體的變化

器作了實驗研究，在速度比較低即空泡數比較大的情況下，給出了一系列的經驗公式。賈力

平、王海斌對空化器與阻力關系進行了研究。袁緒龍、張宇文應用商業軟件FLUENT建立

多相流CFD模型對典型空化器進行了數值模擬。

以往對空化器進行的數值模擬都是針對二維軸對稱體空化器進行研究。但是，當考慮直力

影響以及帶攻角空化器生成的超空泡問題時,二維仿真不能解決此問題，因此本文針對三維

圓盤空化器超空化流動進行數值模擬。首先,對在不同空化數下無攻角圓盤空化器生成超

空泡進行仿真,將仿真結果同Savchenko經驗公式進行對比，總結其變化規律;其次,針對帶

攻角圓盤空化器進行數值模擬,將模擬結果與無攻角圓盤空化器所得結果進行對比,總結其

變化規律。文中主要是針對帶攻角圓盤空化器進行研究。為了解決三維數值模擬中龐大計

算量問題和獲得更好的仿真結果，本文采用構架于深騰1800高性能服務器上的商業軟件

CFX進行仿真。選擇深騰I800高性能服務器是利用其多個高速CPU以及大容量內存的優

點解決龐大計算量問題，而擇商業軟件CFX是因為其高效的并行計算方法，以及可以獲得較

精確仿真結果。

1控制方程和算法

運用N-S方程及b￡湍流模型對均質流場求解，引入Raylcigh-Plcssct方程模擬超空泡

問題。

1.1連續方程

dpdt+V-(PU)=0,(1)

式中：P表示密度;U表示速度;t表示時間。

1.2流體動力sm廣義源項

{a(pu)+v.(PU?U-u(VU+(VU)T))=SM；P=NPEa=lyaPa：

y=NPEa=lyaua,(2)

式中:a是相的序號且lWaWNp(Np是相數)；u為流體動力粘度;SM為廣義源項；Ya為a

相體積分數。

1.3t+u的選取

{d(pk)dt+V(PUk)=V[(y+utok)Vk]+pk-PE；

d(P￡)&+▽(PU￡)-V[(11+ptaE)V￡]+￡k(CEIpk-C￡2P￡).(3)

式中:pk=utvU(VU+VUT)-

23VU(3uIVU+Pk);ut=CuPk2e;Ce1=1.44;Ce2=1.92;Cn=0.09;ok=1.0;0￡=1.3.

1.4密度和流場密度

自然空化狀態是相變過程:為描述均質多相流中此相變過程定義傳輸方程：

(YaPa)+V(yapaUa)=Sa+ra,(4)

式中:Sa為廣義源項；「a為由于相間質量傳遞產生的進入a相的單位體積質量源。其中

ra=NpEb=irab(rab：ib相到a相單位質量流率,本文中a相為水,b相為水蒸汽),

rab=*mabAab,(5)

式中:*mab是b相到a相單位界面密度的單位質量流率;Aab為相間界面密度。

根據空泡動力學知,Rayleigh-Plesset方程描述空泡與水之間的相變過程。方程如下：

RBd2RBdt2+32dRBdt2+2TPLRB=pc-pPL,(6)

式中:RB為氣泡半徑；vL為運動粘度；PL為液體密度；T為表面張力;p為流場參考壓

力;pc為水的蒸汽壓。

在忽略表面張力以及高次項后化簡公式dRBdt=±J23pc-ppL,令蒸汽體積分數

Yb=VBNB=43nRB3NB,其中NB、VB分別是單位體積氣泡數和氣泡體積。

單位體積總質量傳遞率

mlg=NBdmBdt=F3ybPgRB23pc-p|PLsgn(pc-p),(7)

式中:F為通過實驗獲得的經驗系數,Fvap=50,Fcond=0.01.令代Ynuc(l-yb)替yb,則

m'lg=F3ynuc(1-yb)PgRB231pc-p|PLsgn(pc-p),(8)

式中：Ynuc為成核位置的體積率。通常在數值模擬中計算空泡模型時采用人為給定單位體

積氣泡數NB=9X106(單位:個)，而對于RBGmm的氣泡無法得到好的仿真結果。采用

丫nuc(l-yb)代替Yb方法,可以計算小尺寸氣泡半徑如RB=1um,即使氣泡半徑再縮小2

個數量級均可以仿真。這是因為RB=1um,其體積分數yb<10-12,根據方程(8)可知，生成

的空泡由成核位置的體積率Ynuc決定，這樣避免由于計算時產生舍入誤差影響仿真結果。

1.5流域外邊界的確定

基丁有限元的有限體積法:N-S方程中的對流項采用CFX高精度格式,它利用調整混合因子

提高計算精度且平衡收斂精度。以一階迎風為例,對流項

6二樹&VeA「，⑼

式中:巾為變量巾第i個節點值:?為變量力迎風節點值；C為混合因子;V巾為迎風節點的

節點梯度；C.為從迎風節點到i節點的矢量。

針對不同條件下，C會作出相應的調整滿足計算要求獲得精確的解,例如在i節點臨近迎風

節點的節點梯度變化劇烈:則混合因子趨近于0,反之趨近于1.

邊界條件為：上游邊界給定來流速度和各流體體積率；下游邊界給定壓強；流域外邊界夾用

速度為零物面邊界,或者根據文獻采用上游邊界，但是根據模擬研究發現流域外邊界采用上

游邊界時，仿真的收斂效果不如采用物面邊界，即RMS殘差曲線振蕩較劇烈；另一方面,當

v>12m/s,阻塞比(即，模型截面面積與流域截面面積之比)小于0.01時,流域外邊界采用速

度為零的物面邊界和采用上游邊界模擬結果無差別；當阻塞比大于0.01時,流域截面直徑

影響空泡的形態。因此,本研究工作均采用小阻塞比的物面邊界作為流域的外邊界;初始條

件:設定初始來流速度和冬流體體積分數以及流域初始壓力。

三維物體數值仿真具有計算量龐大,精度高等特點,這就使得進行高精度數值仿真時普通計

算機無法完成任務。因此本仿真實驗使用了哈爾濱工業大學能源學院的深騰1800高性能

服務器。

2下覆蓋物的數值模擬

利用商業CFD軟件CFX,針對三維圓盤空化器模型進行數值模擬,所得結果將為水下超空泡

水下高速航行體減阻機理研究提供重要的理論基礎和依據。

2.1空化器網格劃分

運動體幾何形狀如圖1所示，(a)圖表示無攻角圓盤空化器，圓盤空化器尺寸為R=6cm,L=l

cm;(b)圖表示帶a攻角圓盤空化器正視圖。尺寸與無攻角圓盤空化器一致。在仿真時二

者表面均認為是速度為零的光滑物面。模型所在坐標系為三維直角坐標系,XY平面為正視

圖,XZ平面為俯視圖,YZ平面為側視圖。

對于本文問題,若采用非結構化網格,生成的空泡外邊界會出現凸凹不平的毛刺，影響對于

仿真結果的判斷。若想消除毛刺現象,就要對網格進行細分，但會帶來網格質量差和計算量

呈指數倍增長的缺點。因比采用基于六面體結構化網格劃分方法。無攻角圓盤空化器網格

分布圖和流場域網格分布弱如圖2所示。總共六面體網格數為2165124,節點數為2204

421,其中圓盤空化器表面網格數為5378.流域總長為14m,空化器前面流域長度為1口,為

滿足阻塞比要求,流域截面積為圓形,其半徑為L5m.其中流域設定14m是因為在小空化數

下空泡的長度可達十幾米長。如表1,在速度為100m/s時,空泡長度達到11.04m.此外,根

據軟件ICEMCFD評定網格質量方法可知，結構化網格的網格質量大于等于0.2足以滿足仿

真要求。此空化器模型的網格質量為0.23,完全滿足仿真要求。帶攻角圓盤空化器所在流

域的網格數量與其相差無幾,只是空化器帶有攻角，網格質量為0.2.因此,根據本文研究內

容,采用基于六面體結構化網格劃分方法,滿足數值模擬的必要條件。

3初始條件及研究目標的確定

空化數是衡量流體空化程度的一個物理量,計算公式。=(p8-pc)/0.5Pv2,其中

P8,pc,P,v分別為某基準位置上的絕對靜壓力、空泡內絕對壓力、液體密度、相對于物

體的未擾動液體的流速。研究中設定初始條件P8、25°C時的水的pc和P.

雷諾數反應流體粘性的物理量,其公式Re=vooDn/v,在本文研究中其變化范圍為

3.56X106-1.19X107.由于空化數較大情況下無法生成超空泡,反之，則由于情況比較復

雜超出軟件的計算能力，因此設定空化數0.019-0.23為研究目標。由于計算軟件定義的

上游邊界條件和初始條件是以速度或者壓力為參數，所以在仿真的時候需要將空化數轉換

成速度參數，即取速度范圍在30?100m/s范圍內為研究目標,并以每次仿真速度遞增lOm/s

進行采樣。

3.1仿真結果與實驗結果對比

直徑Dn=120mm無攻角圓盤空化器，速度范圍為30?100m/s,仿真結果取空泡形態參數，即空

泡的長度Lc和中截面最大直徑De,見表1,空泡長度與厚度隨著速度地增加而增加，即隨空

化數的減小而變大。表1是以初始速度為30m/s并以10m/s為步長增加至速度為100n/s

的仿真結果.

烏克蘭專家Savchenko依據實驗結果和空泡截面擴大獨立性原理，利用自由流線擴展的漸

近律,得到并用實驗驗證了：在自然空化情形下，空泡的基本幾何形狀是一個近似的橢球體。

基本的超空泡尺寸有中截面直徑De和長度Lc,在2%?4%的精度范圍內可以由以下漸近關

系式決定：

DcDn=CxO(l+o)。，LcDn=l。CxO(l+o)lnl。,(10)

式中：De為中截面直徑即超空泡最大直徑;Lc為超空泡最大長度;Dn為空化器直徑；。為空

化數,CxO空化器阻力系數。

由于本文采用的是圓盤空化器,從文獻中可以得知CxO=(i.82.

為進一步驗證仿真結果的可靠性,將仿真結果與經驗公式(1)進行對比，如圖3、4所示,由

公式所得超空泡長度和中截面最大直徑變化曲線可知,其長度和中截面最大直徑是隨著空

化數的增加而減少,仿真結果的數據點大部分落在該曲線上。即使在較大空化數下，略有不

同，但是均在2%?4%精度范圍內，因此二者具有很好的一致性。其中D'c=Dc/C和

L'c/C分別為經過無量綱化的空泡中截面最大直徑與長度，C為特征長度。此外,對比結

果也說明了使用商業軟件CFX可以進行三維自然超空化的數值模擬研究。它的計算結果可

以為以后的超空泡實驗研究提供參考依據。

以速度為70m/s的仿真結果為例，如圖5所示,空泡最大中截面直徑Dc=O.55m和長度

Lc=4.94m,圖5為三維視圖中的正視圖(XY平面)。選取方法是在25℃水蒸氣的等高圖中

選取最小含水蒸汽率為5佻處定為輪廓邊界，利用軟件自身的測量工具測定De和Ie其中

設定最小含水蒸汽率為5佻是以在上述邊界條件下,速度范圍在30?100m/s中進行仿真結

果與(1)式相比時最匹配為標準而設定的。另一個方面，根據(1)式計算Lc=4.94458m和

Dc=O.55286m,兩者比較之后，誤差在2%?4%精度范圍之內，說明二者結果相吻合,其中自然

超空泡所取正視圖的水蒸汽體積率為50%-100%.圖5中，在超空泡尾部呈現雙渦管閉合方

式,這與文獻中關于超空泡尾部閉合方式相一致。

如圖6所示,是其三維視圖中的三維形狀的超空泡,取含水蒸氣50%的等值曲面圖形。速度

70m/s時3-D超空泡形狀是個橢球體,從其側面圖圖3看到的雙渦管形狀是超空泡尾部形成

高逆壓梯度區域,迫使部分水倒流,形成回注射流。對于回注射流現象在后面會進一步證明。

3.2帶攻角圓空化器的模擬

采用帶攻角圓盤空化器可以為航行體提供必要升力，對于水下航行體航行起控制作用，因此

具有很重要的研究意義。

3.2.1無攻角部分空化器回注射流仿真條件

帶攻角圓盤空化器生成超空泡形態可由圖7(a)?(b)所示，以60°攻角為例,Dn=12cm圓盤

空化器，。=0.219,航行速度v=30m/s,水蒸汽含汽率為53%?】00%.從俯視圖和側視圖對比

可知兩種角度觀測到的空泡厚度是不一致的，而長度一致;空泡尾部均有凹入部分,其含汽

率為0,表明此處的流體為水，因此可以證明產生了十分明顯的回注射流,且形狀不同。

無攻角圓盤空化器生成超空泡形態可由圖7(c)?(d)所示。仿真條件除無攻角外均與上面

一致。與帶攻角圓盤空化器生成超空泡形態比較,從側視圖和俯視圖觀察其空泡厚度相等，

長度也一致,尾部產生的回注射流形狀相同。

如圖