周期性陣風流對通氣航行體超空泡形態及流體動力特性影響

王 威， 王 聰， 李聰慧， 杜嚴峰

(哈爾濱工業大學 航天學院，哈爾濱 150001)

水下航行體借助通氣超空泡的包裹，使得其表面的黏性阻力大幅度減小，航速得到突破性的提高[1]。采用主動通氣的方式形成超空泡的關鍵是要向水下航行體的空化器后部注入空氣或者其它不可冷凝氣體，增大空泡內部的壓力值，使得空泡尺度增加到足以包裹整個航行體表面。主動通氣的方式由Reichard[2]最先提出，與單純的依靠提高航速而生成的自然超空泡相比，通氣超空泡對航行體的運動速度沒有那樣苛刻的要求，更加易于實現，而且通氣量具有可控性，可以調節空泡對航行體表面的包裹范圍，從而改變航行體的沾濕面積。Savchenko[3]提出了航行體在超空泡內部運動的四種模式，指出了沾濕區域可以為航行體提供升力，維持其縱平面運動平衡的作用。

目前對于航行體與空泡相互作用的研究多數是在穩定來流中改變航行體的攻角實現的，比如Lee等[4]在水洞實驗段中安裝了弧形的模型支撐機構，通過此機構連接的支撐桿控制航行體模型在水洞中做俯仰運動，實現航行體和超空泡的相互作用。Pan等[5]根據實驗結果，設計了可以實現通氣航行體俯仰運動的動態計算網格，通過數值計算的方式研究了通氣超空泡內部壓力變化特點。Yu等[6]利用數值計算方法，進一步研究了航行體俯仰運動時流體動力變化的特點。若是運動中的通氣超空泡航行體遇到前方不穩定的來流，或是受到海面附近的波浪影響，其空泡形態往往會產生變化，影響到航行體表面的沾濕區域，進而改變航行體的流體動力，過大的來流擾動甚至可能造成空泡的潰滅。

針對這樣的問題，明尼蘇達大學的Arndt等[7]在水洞實驗段中安裝了由兩片水翼組成的陣風發生器(Gust generator)，依靠上游水翼的往復擺動使下游流域產生周期性的陣風流動，研究了穩定流和陣風流作用下通氣超空泡的生成和保持機制，并給出了通氣率范圍。Lee等[8-9]實驗發現了陣風流作用下的通氣超空泡形態變化特性，發現在陣風發生器的頻率增大到一定程度時，通氣空泡的長度將會縮短；Lee比較了前支撐和后支撐的方式對通氣超空泡的影響，前支撐方式的支撐結構會在空泡上游造成擾動，使空泡壁面不光滑，而后支撐方式會抑制空泡變形，并且影響空泡尾部閉合。Sanabria等[10]實驗研究了超空泡航行體的動力學建模與控制，結果表明陣風流作用下，航行體空化器受到的擾動可以忽略，但是尾翼因擾動產生的流體動力不可忽略。Karn等[11-12]實驗研究了超空泡尾部的幾種不同形式的閉合模式，研究了弗勞德數Fr和通氣率系數Cq對這些閉合模式的影響，并分析了空泡內外壓力變化和閉合轉變的關系。Karn等[13]研究了陣風流作用下的通氣率、流場速度、空化器尺度、流動不穩定性對通氣超空泡產生和潰滅的影響，指出陣風流動會降低氣泡的凝結效率，陣風振幅的增加導致空泡形成和潰滅的通氣率略微單調增加。

目前研究陣風流對通氣超空泡影響大多是利用水洞實驗方式對單獨空化器產生的通氣超空泡進行的，主要是圍繞空泡尺度及空泡尾部閉合機理展開的，受到實驗條件的限制，陣風流場中的垂直速度變化采用單點測量的方式，難以捕捉整個流場垂直速度的變化特性，實驗中的模型需要支撐結構固定，而支撐結構對空泡形態存在一定影響，并且陣風流作用下通氣超空泡航行體的流體動力變化特性還沒有相關研究涉及到。基于以上原因，采用數值模擬的方法，對陣風流作用下通氣超空泡航行體的流體動力進行研究，不僅可以捕捉到整個陣風流場的垂直速度分布，而且可以避免支撐結構對空泡形態的影響，了解航行體的空泡形態變化細節，掌握航行體沾濕區域及流體動力的變化規律。

1 數學計算模型

1.1 控制方程

由于來流速度較低，流域的壓力都遠大于水的飽和蒸汽壓，自然空化現象較弱，所以忽略自然空化的影響[14-15]，主動通氣形成的超空泡內部基本被空氣填滿[16]。參照文獻[17]對通氣空泡流動的計算經驗，本文采用VOF(Volume of fluid)兩相流模型對通氣超空泡進行數值模擬研究。

1.1.1 連續性方程

(1)

1.1.2 動量方程

(2)

1.1.3 體積分數方程

(3)

式中:αg為氣體體積分數，αg+αl=1,αl為水體積分數，如果網格單元充滿水，則αg=0；如果充滿氣體，則αg=1；如果網格單元中包含氣體和水的分界面，則0<αg<1。

1.1.4 湍流方程

根據文獻[18]對通氣空化流動的計算方法，本文選擇RNGk-ε湍流模型進行數值模擬研究，此模型是由Yakhot和Orzag提出的，它可以很好的處理高應變率和流線彎曲程度較大的流動，k和ε的輸運方程如下

Gk+Gb-ρmε-YM+Sk,

(4)

(5)

1.2 驗證計算模型及邊界條件設置

圖1 空化器及陣風發生器的布置

整個計算流域如圖2所示。水翼、空化器模型和支撐桿，采用壁面邊界條件；流域左側及流域四周邊界均采用速度入口邊界條件；流域右側采用壓力出口邊界條件；通氣口采用質量流量入口邊界條件，通氣率Cq=Q/(V∞D2)=0.15，其中Q是通氣率，V∞為來流速度，D為模型最大截面直徑(單獨空化器時為空化器直徑Dn)；陣風發生器外圍的柱面將整個流域分為兩個區域，這兩個區域采用一對重合的圓柱交界面分隔，當陣風發生器擺動時，這對圓柱交界面之間相互滑動，而兩個區域的網格并不發生變化。整個計算域采用六面體網格離散，數值計算開始前，采用不同網格尺度和時間步長進行網格獨立性檢測，確保模擬結果對網格密度和時間步長不再敏感。

圖2 計算域及邊界條件

1.3 航行體計算模型及邊界條件設置

如圖3所示，航行體模型由圓盤空化器、錐柱段彈身、圓柱形尾噴管組成。航行體模型的柱段直徑D=40 mm，全長L=12D，其中圓盤空化器直徑Dn=0.3D，航行體質心距離空化器6.5D,通氣口距離空化器D/10，監測點的布置和圖1的方式相同。航行體模型的計算域如圖4所示，陣風發生器的水翼可以獨自轉動，模型后部沒有支撐桿，邊界條件設置與圖2中的條件相同。

圖3 計算模型

圖4 航行體模型的計算域及邊界條件

Fig.4 Computational domain and boundary conditions of vehicle

2 模型驗證

在頻率f=20 Hz，最大擺動角度θ0=±6°的條件下，監測點B的垂直速度計算結果和Lee的實驗結果對比如圖5所示。可以看出監測點B的垂直速度在陣風發生器的作用下呈現周期性變化，數值計算結果和實驗結果有很好的一致性。Lee的研究結果表明監測點的垂直速度變化頻率等于陣風發生器的擺動頻率f，數值模擬得到的陣風流作用下通氣超空泡形態變化和實驗結果對比如表1所示。兩者空泡形態變化的規律一致，充分驗證本文數值模擬方法的有效性。

圖5 監測點的垂直速度

表1 數值模擬結果和實驗結果[8]對比

3 計算結果及分析

3.1 流域中波的特性

上游陣風發生器的擺動會在下游流場中將產生陣風流動，Lee通過水洞實驗研究陣風發生器的頻率和最大擺動角度對流域下游固定監測點B的最大垂直速度v的影響。受到實驗條件的限制，實驗中以單點測量不能夠顯示流場中不同位置的垂直速度分布情況，利用數值模擬方法可以很好地解決這一問題。在下游自空化器端面沿流動方向布置一條直線，用于監測流域中垂直速度的變化，具體布置如圖6所示。

圖6 流場監測線布置俯視圖

Fig.6 The top view of the arrangement of the flow field monitoring line

在陣風發生器擺動幾個周期以后，當水翼再次回到平衡位置時，流域中垂直速度的分布如圖7和圖8所示。

Lee的研究表明，流域中的波長λ=V∞/f，從圖7可知，來流速度V∞相同時，陣風流的頻率f越大，流域中的垂直速度的波峰之間距離越小，即波長變短，而垂直速度的峰值越高，即波幅越大。

從圖8可知，陣風流的頻率f相同時，來流速度V∞越大，流域的波長越大，且垂直速度的峰值也較高，即波幅越大。因為較大的來流速度可以把陣風發生器產生的擾動能量迅速的傳遞到下游流域，從而產生較大的波幅，而流速較小的情況，陣風發生器擾動的能量還沒有及時傳遞到下游就已在水中因摩擦阻力的作用而快速衰減。

圖7 不同陣風流頻率下的垂直速度分布(V∞=20 m/s)

Fig.7 Vertical velocity distribution at different frequency of gust flow (V∞=20 m/s)

圖8 不同來流速度下的垂直速度分布(f=20 Hz)

Fig.8 Vertical velocity distribution at different flow velocities (f=20 Hz)

Lee的研究發現短波對空泡長度的影響較大，從表1中也可以看出短波會加劇空泡的變形效果，為了研究陣風流作用下航行體流體動力變化，選擇短波流動影響可以讓通氣超空泡變形，使航行體表面接觸到流動介質，產生更大的流體動力變化。波長由陣風發生器擺動頻率f和流速V∞的影響,從圖8的討論中可以發現給定的流速不宜過低;在低流速情況下，陣風發生器下游的波幅沿流動方向衰減較快，不適合研究陣風流動對通氣超空泡航行體的影響。比較穩定流條件下不同流速的通氣超空泡形態如表2所示。從表2可知,在來流速度較低的情況，通氣航行體尾部已經形成了沾濕區域，這種條件不利于比較陣風流作用下形成沾濕區域的異同，而流速較高的情況，航行體模型沒有接觸到流動介質，因此在后面的計算中我們選擇來流速度V∞=20 m/s，陣風發生器的頻率f=20 Hz的工況進行研究。

3.2 流體動力特性

文中的阻力系數Cd、升力系數Cl、力矩系數Cm分別表示為

表2 穩定流條件下，不同流速下的空泡形態

Tab.2 The shape of cavities at different flow rates under stable flow conditions

流速V∞/(m·s-1)空泡形態

其中，Fx和Fy為x方向和y方向的阻力和升力；Mz為作用于航行體質心位置z方向的俯仰力矩；S為航行體柱段截面面積，L為航行體的長度。

陣風發生器擺動幾個周期后，水翼再次回到平衡位置時，開始監測航行體的流體動力系數，無量綱時間T=t/T0，t為陣風發生器的擺動時間，T0為陣風流的周期。

三個周期內的流體動力變化結果如圖9所示。陣風流動的作用下，航行體的流體動力系數出現了周期性變化，峰值左右的曲線并不對稱，因為流場中重力影響使得陣風發生器向上和向下擺動時對航行體造成的擾動效果不同，從而影響到通氣空泡形態及流體動力變化。另外的一個原因是文獻[19]所提到的航行體附加質量力的作用，當陣風發生器高頻率擺動時，流體對航行體作用的強度變大，航行體的附加質量力對流體動力的影響會增大，使航行體的流體動力曲線峰值向前偏移。為了觀察一個周期內空泡形態的變化，比較圖9中的9個特征時刻(a～i)的超空泡形態如表3所示。

圖9 航行體的流體動力系數

表3 一周期內不同時刻的超空泡形態

Tab.3 The supercavitation shape at different time in a period

時刻通氣超空泡形態及沾濕區域

從表3可知:①a時刻處于t=1T，航行體尾部的下表面出現沾濕區域，阻力系數不是最小值，升力系數為正值，升力作用點位于質心后方，所以力矩系數為正(力矩以使航行體逆時針轉動為正)；②b時刻處于t=1.1T，隨著超空泡形態的變化，航行體尾部的沾濕區域逐漸的消失，此時航行體的阻力系數最小，升力和力矩系數接近0，并在b時刻附近保持一段時間；③c時刻處于t=1.15T，此時航行體的上表面出現了一小部分沾濕區域，因此阻力系數有所回升，而升力系數開始向負值方向發展，由于升力的作用點位于航行體質心前部，所以力矩系數會向正方向有一個小波動，然后隨著空泡形態的變化，沾濕區域逐漸擴大并沿航行體上表面向尾部位置移動，那么力矩系數向負值方向增大；④d時刻位于t=1.3T，此時航行體上表面的沾濕區域達到最大，阻力系數達到峰值，而升力和力矩系數達到負向最大值；⑤e時刻位于t=1.45T，此時航行體上表面的沾濕區域逐漸減小，因此流體動力系數也隨之減小；⑥f時刻位于t=1.55T，此時航行體沾濕區域再次消失，阻力系數隨之達到最小值，升力和力矩系數接近0，并在f時刻附近保持一段時間；⑦g時刻位于t=1.65T，此時航行體下表面開始出現沾濕區域，阻力系數逐漸回升，升力系數向正值方向增加，此時的升力作用點位于航行體質心的前部，所以力矩系數會向負方向有一個小波動，然后隨著超空泡的形態變化，沾濕區域逐漸擴大并沿航行體下表面向尾部移動，因此升力系數向著正向增大；⑧h時刻位于t=1.75T，空泡形態因重力作用會上漂，航行體下表面的沾濕區域比下表面的沾濕區域大，所以此時的沾濕區域面積比d時刻大，阻力系數的峰值也比d時刻的情形要高，升力系數和力矩系數接近正向最大值，并且這個峰值要比航行體上表面沾濕的d時刻要大；⑨i時刻位于t=2T，此時剛好過了一個完整的周期，航行體尾部下表面形成沾濕區域，所以流體動力系數與a時刻相當。

3.3 航行體尾部沾濕區域的流體動力特性

通氣超空泡航行體受到空泡的包裹，流體動力系數會大幅度減小，文獻[3]指出了沾濕區域可以為航行體提供升力的作用，但是陣風流作用下，沾濕區域為航行體提供的流體動力的效果并不清楚，定義空化器阻力占總阻力的百分比為fcd；沾濕區域阻力占總阻力的百分比為fwd、升力占總升力的百分比為fwl；沾濕區域面積占航行體總面積的百分比為fA，時刻a～時刻i的這些比值情況如表4所示。

表4a～i時刻的空化器及沾濕區域的流體動力比例/%

Tab.4 The hydrodynamic ratio of cavitator and the wetted area in the working timea-i/%

時刻fcdfwdfwlfAa101.43.499.24.8b000c99.70.34.80.4d90.19.999.812.7e102.22.099.32.5f103.9000g98.02.022.54.3h86.113.999.919.3i101.03.299.14.0

從表4可知，空化器阻力比值fcd中，時刻(a、b、e、f、i)的數值達到100%以上，這是因為超空泡內部的回射流對航行體的前進方向有一定驅動作用，文中i時刻航行體表面的流線如圖10所示。回射流作用使航行體總阻力減少，而空化器本身所受到的阻力是航行體阻力的主要部分，它受到回射流的影響較小，所以空化器所受到的阻力比例超過100%；而時刻d、時刻h的沾濕區域較大，使得空泡內部的回射流作用減弱，另外沾濕區域也分擔了部分阻力，導致空化器所占的阻力比例有所減小。

圖10 航行體表面回射流

沾濕區域阻力比值fwd中，時刻d、時刻h的沾濕區域阻力比例接近或超過10%，它們的沾濕區域面積占航行體表面積的12.7%和19.3%，可見沾濕區域面積對阻力的影響很大，沾濕區域的阻力比例隨著沾濕區域面積的增大而增大，而其它時刻的沾濕區域面積較小，所以阻力比例也小。

沾濕區域升力比值fwl中，時刻(a、d、e、h、i)所占有的比例接近100%，例如i的沾濕區域面積只占有航行體總面積的4%，但卻提供了99.1%的升力，而c和g的沾濕區域占航行體總面積的比例雖小，但升力卻占有總升力的比例較大，可見沾濕區域對航行體尾部升力的影響很大。

3.4 沾濕區域的壓力分布特性

從上面的分析可以發現，沾濕區域的升力占有總升力的比率較大，以g時刻為例，其沾濕區域位于航行體的下表面，分析該時刻g的沾濕區域的特性如圖11所示。

(a) 沾濕區域位置

(b) 空泡內部流線圖

(c) 沾濕區域壓力云圖

圖11(a)中標注了沾濕區域的具體位置并用網格線表示，從圖11(b)可知,空泡內部的流線分布特點，在沾濕區域附近的流速很高，接近遠場來流的速度V∞，高速區域基本處在沾濕區域的附近，因沾濕區域脫離了空泡的包裹直接與外部的流動介質接觸，所以在沾濕區域的邊緣形成了速度很高的回射流，而遠離沾濕區域的位置回射流減弱；從圖11(c)可知,模型表面沾濕區域的最大壓力出現在沾濕區域的前部位置，處于空泡閉合線左側，流場中的回射流使該位置的壓力增大，而且該壓力峰值比沾濕區域的壓力還大。

如果航行體被整個空泡包裹，表面沒有出現沾濕區域的工況f中航行體表面的壓力處于35～36 kPa，而時刻g出現了沾濕區域，其壓力最大值達到了61 kPa，正是這樣的高壓使得沾濕區域的面積雖小，但卻可以給航行體的表面帶來很大的壓力，進而影響航行體的流體動力變化。

4 結 論

本文使用動態網格技術數值模擬研究了周期性陣風流作用下通氣超空泡航行體的空泡形態及流體動力系數的變化過程，獲得的主要結論如下：

(1) 周期性陣風流作用下，通氣超空泡航行體的空泡形態及流體動力呈現周期性變化，阻力和升力系數隨著沾濕區域的面積的增加而增大；一個周期內，沾濕區域相對于航行體質心的位置變化會使力矩系數出現兩次小幅度波動。

(2) 周期性陣風流作用下，空化器阻力占航行體總阻力的比例較大；沾濕區域阻力占航行體總阻力的比例隨著沾濕區域面積增大而增加，沾濕區域升力占航行體總升力的比例較大。

(3) 沾濕區域前端的空泡閉合線附近，受到回射流影響出現了高壓區，使得沾濕區域的面積雖小，卻可以為航行體提供很大的升力。