基于勢流理論的尾空泡對航行體表面壓力影響研究

尤天慶，權曉波，劉元清，鮑文春，王占瑩

（1.北京宇航系統工程研究所，北京 100076；2.中國運載火箭技術研究院，北京 100076）

0 引 言

航行體水下垂直干發射情況下，彈射氣體會附著在尾部低壓區，形成附著尾空泡[1]。考慮不可壓縮流體下游的擾動會對上游流動產生影響，尾空泡狀態及演化過程必然會對上游航行體表面壓力產生不同程度的干擾[2]，進而對航行體整體流體動力特性產生影響[3]。

目前對尾空泡的研究主要集中于空泡形態和泡內壓力的預示，其計算分析方法主要分為簡化工程計算[4]以及基于雷諾平均N-S方程的流場數值模擬兩類[5]。針對流線型航行體及尾空泡組合體，在其流動分離不明顯的情況下，基于勢流的邊界元計算可在有限計算資源下分析流場特性，從計算效率和流場信息獲取等方面比較，是一種介于工程簡化計算和基于雷諾平均N-S方程仿真計算的方法。

運用勢流理論對空泡發展演化的研究，一方面集中于水中孤立氣泡的生長和潰滅[6]，最早見于Blake和Gibson的工作[7]，另一方面集中于航行體肩部附著空泡研究[8-9]。相關研究表明基于勢流的邊界元方法，在空泡形態、流場結構以及物面整體受力等方面的研究分析工作中發揮了重要作用。

本文基于勢流理論，利用邊界元數值方法，建立了二維軸對稱航行體及附著尾空泡的流場計算方法，并形成了計算程序。在此基礎上計算分析了不同空泡數條件下尾空泡對上游航行體物面壓力分布的影響，并以此為基礎進一步分析了瞬態情況下尾空泡收縮過程對上游壓力的擾動，基于數值計算結果結合實驗數據說明了壓力擾動的衰減規律。

1 數值計算方法

1.1 控制方程

假設空泡附著在航行體后部擾流低壓區，在無限大流體域內發展演化。并假設流體無粘、無旋且不可壓，則N-S方程可簡化為拉普拉斯方程。

式中，V為速度向量，U∞為來流速度（設來流為x軸正向），Φ為總速度勢，x為位置坐標，φ為擾動速度勢。

在物面，流體需滿足不可穿透條件，即壁面法向速度為零：

式中，n為物面法向方向向量。

在空泡壁面，流體壓力在流線方向上為恒定，在計算中主要通過伯努利方程計算空泡切向速度，進而通過延空泡壁面流線方向，進行積分獲得空泡壁面上的速度勢分布。t為物面切向方向向量，lc為空泡壁面流線方向長度。

對于空泡周圍流場的瞬態變化，速度勢Φ、速度V及壓力p需滿足伯努利方程：

1.2 數值算法

設空間場點為p，物面和空泡壁面上分布的點源位置為q，其強度為λ(q)，則其在p點產生的速度勢和速度分別為

式中，G(p,q)為格林函數。對于拉普拉斯方程，三維情況下格林函數為G(p,q)=-1/4πr，其中r=|p-q|。

從拉普拉斯方程邊值問題角度，本文為第三類邊界條件，即混合邊界條件。航行體表面為流動不可穿透條件，即第二類邊界條件（Newmann條件）。數值計算中，航行體與空泡表面被劃分為N個單元，對航行體表面的第i個單元，考慮在邊界上源匯自身對速度的影響，擾動勢所滿足的邊界元積分控制方程如式（8）所示，式中Sj為單元面積，ni為單元法向向量。

空泡壁面為第一類邊界條件（Dirichlet條件），考慮沿航行體子午線方向，尾空泡分離點位于第d個單元，其擾動勢所滿足的邊界元積分控制方程如式（9）所示，式中ti為單元切向向量。

式（8）和（9）即為邊界元計算的基本方程，通過這兩個方程即可求得物面及空泡壁面的λ分布，進而確定流場中速度勢、速度及壓力。

為節省計算資源提高計算效率，本文采用二維軸對稱計算模型，在航行體及空泡表面的邊界單元為圓環狀。與速度勢及速度相關的源匯均勻沿圓環分布，如圖1所示。對于子午線長度為l圓環邊界單元，其影響速度勢和速度由圓環上源匯積分獲得，其表達式如下：

圖1 軸對稱單元積分示意圖Fig.1 Configuration of axisymmetric element integration

由于第一類完整橢圓積分函數的奇異性，當R趨近于1時，式（9）～（11）所示的邊界單元速度勢和速度影響函數值趨近于無窮大，但其積分是有限的。本文利用偶數積分點的Gauss積分進行邊界單元積分計算。

為驗證數值計算代碼的正確性，計算一無粘無旋理想流體條件下軸對稱旋成體零攻角物面的壓力分布，并與水洞測壓數據對比。對比表明數值計算與數據吻合較好，如圖2所示。

圖2 數值計算與試驗數據對比Fig.2 Comparison between calculation and experiment results

1.3 尾空泡閉合模式

針對物面階梯后部低壓區附著空泡尾部呈現出很明顯的湍流，如果利用勢流理論對該問題進行計算，就必須將這一區域進行理想簡化處理。目前針對不同類型的空泡，已有多種空泡閉合模式[10]。本文根據附著空泡的形態特征，考慮數值計算的穩定性，采用了圓球壁面條件對空泡尾部流場進行了簡化處理，如圖3所示。空泡壁面邊界條件利用式（8）計算，圓球壁面利用式（9）計算。

2 數值計算結果分析

2.1 穩態計算

采用上述數值計算方法，對不同來流條件下軸對稱半球頭型航行體后部附著尾空泡的二維軸對稱流場進行了迭代計算，獲得收斂后的空泡形態如圖4所示。σ=(p∞-pc)/0.5ρU2∞為空化數，p∞為環境壓力，pc為尾空泡內部壓力，ρ為流體密度，U∞為來流速度大小，X/L和Y/L分別為無量綱的軸向和徑向坐標，L為航行體長度。X/L≤0部分為航行體物面，X/L＞0部分為不同空化數下的附著尾空泡形態。計算結果表明穩態空泡形態與空化數緊密對應，即針對不同的環境壓力、泡內壓力或來流速度，只要空泡數一致，勢流條件下穩態空泡形態即確定。不同空泡數計算結果表明，隨著空泡數的降低，尾空泡長度增加，尾空泡輪廓曲線曲率變化趨緩。

圖4 穩態計算所形成的附著空泡Fig.4 Tail bubble profile of steady calculation

不同空化數條件下航行體和空泡表面的壓力系數分布如圖5所示。與穩態空泡形態和空泡數的對應關系相同，穩態計算的壓力系數分布也與空泡數緊密相關。不同空泡數的數值計算表明，隨著空泡數降低，尾空泡尺寸增大，壓力系數平直段也相應地延長。不同空化數條件下，遠離尾空泡的航行體沾濕區壓力保持一致。尾空泡主要對上游臨近尾空泡的物面壓力產生較大影響，如圖5中放大圖所示。數值計算表明，隨著空化數提高，泡內壓力系數逐漸降低，導致臨近尾空泡的物面壓力越低，其影響范圍也相應擴大。

圖5 沿航行體及空泡表面壓力分布Fig.5 Pressure distribution along vehicle and bubble meridian

2.2 瞬態計算

為分析空泡瞬態變化對物面壓力的影響，采用與穩態計算相同的航行體外形，在穩態計算結果上，降低空泡內部壓力，計算空泡收縮過程流場變化。空泡收縮過程形態變化如圖6所示。由于空泡內部壓力降低，空泡收縮長度縮短，空泡輪廓曲線變化趨于劇烈。空泡變化所引起的流場速度及加速度改變，進一步導致了上游物面壓力的變化。下面主要針對空泡收縮過程的t0+dt和t0+4dt兩時刻流場參數變化進行分析。

圖6 空泡收縮過程形態變化Fig.6 Cavity profile change while shrinking

由于空泡收縮過程中，整個航行體物面和空泡組合體形態發生變化，其長度更短，流場下游空泡輪廓線曲率變化更加劇烈，因此在接近尾空泡的航行體表面流體流速增大。由此動壓系數1-q/q∞在X/L=0附近明顯降低，如圖7所示。其中，q為流體當地處動壓，q∞為無窮遠處動壓。

圖7 收縮過程不同時刻物面及空泡表面動壓分布對比Fig.7 Comparison of dynamic pressure distribution

針對尾空泡收縮瞬態變化過程，本文流場壓力計算采用非穩態伯努利方程，如式（5）所示。基于此，壓力系數為Cp=1-q/q∞-ρΦ˙/q∞。式中，ρ為流體密度，Φ˙=?Φ/?t為速度勢變化率。可見引起壓力變化的因素為兩部分：1-q/q∞為動壓變化所導致的壓力系數變化，其分布如圖7所示；-ρΦ˙/q∞為速度勢變化率導致的流場壓力變化。兩者綜合作用表現為物面壓力系數的分布，如圖8所示。

對比圖7和圖8，在靠近尾空泡的航行體物面區域X/L=0附近，動壓導致的壓力明顯降低，而實際此段壓力在空泡收縮過程中壓力明顯升高，兩者變化趨勢截然相反。導致這種差異的原因是速度勢變化率導致的流場壓力變化，即壓力計算中的-ρΦ˙/q∞，其分布如圖9所示。圖9表明，在空泡收縮過程中，其所引起的流場中速度勢變化率對空泡上游物面壓力有著較大影響，其影響超過流場中速度變化的影響，導致靠近空泡部分物面壓力的上升。速度勢變化率Φ˙=?Φ/?t為歐拉方程中流體加速度項沿流線方向的積分。計算中速度勢變化率導致的流場壓力變化，本質上是由流體加速度所引起。

圖8 收縮過程壓力分布對比Fig.8 Comparison of static pressure distribution

圖9 收縮過程速度勢變化率引起的壓力分布對比Fig.9 Time derivative of the velocity potential

對于圖10顯示的空泡形態變化，初始時刻泡內壓力越低，尾空泡附近流體質點運動加速度就越大。在相同時間內速度變化量也越大，尾空泡收縮程度越劇烈。

在圖10顯示的壓力變化方面，由于尾空泡內壓力在空間上是均勻的，在其收縮過程中空泡壁面壓力整體升高，即X/L＞0部分壓力變化ΔCp在空間上也均勻分布。尾空泡收縮過程引起流場上游航行體壁面壓力變化，這種對流場的擾動在接近尾空泡處最大，在遠離尾空泡過程中逐漸衰減，衰減趨勢近似呈雙指數函數。

圖10 不同初始壓力尾空泡收縮引起的物面壓力變化Fig.10 Pressure change during the tail cavity shrinking process with different initialpressures

為驗證以上壓力擾動空間分布規律的正確性，提取了航行體水下附著尾空泡收縮過程上游物面壓力測點變化情況，如圖10中的實驗數據點。對比結果表明，數值計算壓力空間衰減趨勢與試驗吻合較好。

不同尾空泡收縮初始泡內壓力會導致流體運動加速度產生明顯差異，進而影響尾空泡瞬態收縮過程中的物面壓力，如圖10所示。尾空泡收縮初始泡內壓力越小，流場中流體運動加速度就越大，對上游航行體物面壓力的擾動也越大。

3 結 論

針對水下垂直干發射尾空泡對上游物面壓力影響問題，本文基于勢流理論，考慮航行體物面及空泡壁面混合邊界條件，建立了適用于尾空泡流場發展演化的二維軸對稱邊界元數值計算方法，對航行體及尾空泡流場進行了計算分析，得到的結論如下：

（1）不同狀態尾空泡流場穩態數值計算表明，隨著空化數的提高，尾空泡壓力系數逐漸降低，導致臨近尾空泡的物面壓力系數下降，尾空泡對上游物面壓力影響范圍也相應擴大；

（2）尾空泡收縮瞬態流場演化計算表明，收縮引起的流場加速度對空泡上游物面壓力有著較大影響，導致流場中靠近尾空泡物面壓力的上升；

（3）不同尾空泡收縮過程初始壓力計算表明，尾空泡收縮初始泡內壓力越小，流場中流體運動加速度就越大，對上游航行體物面壓力的擾動也越大，這種壓力擾動在遠離尾空泡過程中呈雙指數函數逐漸衰減。