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

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

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

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

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

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

1 數(shù)學計算模型

1.1 控制方程

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

1.1.1 連續(xù)性方程

(1)

1.1.2 動量方程

(2)

1.1.3 體積分數(shù)方程

(3)

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

1.1.4 湍流方程

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

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

(4)

(5)

1.2 驗證計算模型及邊界條件設(shè)置

圖1 空化器及陣風發(fā)生器的布置

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

圖2 計算域及邊界條件

1.3 航行體計算模型及邊界條件設(shè)置

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

圖3 計算模型

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

Fig.4 Computational domain and boundary conditions of vehicle

2 模型驗證

在頻率f=20 Hz，最大擺動角度θ0=±6°的條件下，監(jiān)測點B的垂直速度計算結(jié)果和Lee的實驗結(jié)果對比如圖5所示?？梢钥闯霰O(jiān)測點B的垂直速度在陣風發(fā)生器的作用下呈現(xiàn)周期性變化，數(shù)值計算結(jié)果和實驗結(jié)果有很好的一致性。Lee的研究結(jié)果表明監(jiān)測點的垂直速度變化頻率等于陣風發(fā)生器的擺動頻率f，數(shù)值模擬得到的陣風流作用下通氣超空泡形態(tài)變化和實驗結(jié)果對比如表1所示。兩者空泡形態(tài)變化的規(guī)律一致，充分驗證本文數(shù)值模擬方法的有效性。

圖5 監(jiān)測點的垂直速度

表1 數(shù)值模擬結(jié)果和實驗結(jié)果[8]對比

3 計算結(jié)果及分析

3.1 流域中波的特性

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

圖6 流場監(jiān)測線布置俯視圖

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

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

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

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

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

3.2 流體動力特性

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

表2 穩(wěn)定流條件下，不同流速下的空泡形態(tài)

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

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

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

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

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

圖9 航行體的流體動力系數(shù)

表3 一周期內(nèi)不同時刻的超空泡形態(tài)

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

時刻通氣超空泡形態(tài)及沾濕區(qū)域

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

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

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

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

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)的數(shù)值達到100%以上，這是因為超空泡內(nèi)部的回射流對航行體的前進方向有一定驅(qū)動作用，文中i時刻航行體表面的流線如圖10所示?；厣淞髯饔檬购叫畜w總阻力減少，而空化器本身所受到的阻力是航行體阻力的主要部分，它受到回射流的影響較小，所以空化器所受到的阻力比例超過100%；而時刻d、時刻h的沾濕區(qū)域較大，使得空泡內(nèi)部的回射流作用減弱，另外沾濕區(qū)域也分擔了部分阻力，導(dǎo)致空化器所占的阻力比例有所減小。

圖10 航行體表面回射流

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

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

3.4 沾濕區(qū)域的壓力分布特性

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

(a) 沾濕區(qū)域位置

(b) 空泡內(nèi)部流線圖

(c) 沾濕區(qū)域壓力云圖

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

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

4 結(jié) 論

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

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

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

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