尤天慶，權(quán)曉波，劉元清，鮑文春，王占瑩

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

0 引 言

航行體水下垂直干發(fā)射情況下，彈射氣體會附著在尾部低壓區(qū)，形成附著尾空泡[1]?？紤]不可壓縮流體下游的擾動會對上游流動產(chǎn)生影響，尾空泡狀態(tài)及演化過程必然會對上游航行體表面壓力產(chǎn)生不同程度的干擾[2]，進而對航行體整體流體動力特性產(chǎn)生影響[3]。

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

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

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

1 數(shù)值計算方法

1.1 控制方程

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

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

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

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

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

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

1.2 數(shù)值算法

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

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

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

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

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

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

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

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

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

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

1.3 尾空泡閉合模式

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

2 數(shù)值計算結(jié)果分析

2.1 穩(wěn)態(tài)計算

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

圖4 穩(wěn)態(tài)計算所形成的附著空泡Fig.4 Tail bubble profile of steady calculation

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

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

2.2 瞬態(tài)計算

為分析空泡瞬態(tài)變化對物面壓力的影響，采用與穩(wěn)態(tài)計算相同的航行體外形，在穩(wěn)態(tài)計算結(jié)果上，降低空泡內(nèi)部壓力，計算空泡收縮過程流場變化?？张菔湛s過程形態(tài)變化如圖6所示。由于空泡內(nèi)部壓力降低，空泡收縮長度縮短，空泡輪廓曲線變化趨于劇烈?？张葑兓鸬牧鲌鏊俣燃凹铀俣雀淖儯M一步導(dǎo)致了上游物面壓力的變化。下面主要針對空泡收縮過程的t0+dt和t0+4dt兩時刻流場參數(shù)變化進行分析。

圖6 空泡收縮過程形態(tài)變化Fig.6 Cavity profile change while shrinking

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

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

針對尾空泡收縮瞬態(tài)變化過程，本文流場壓力計算采用非穩(wěn)態(tài)伯努利方程，如式（5）所示?；诖耍瑝毫ο禂?shù)為Cp=1-q/q∞-ρΦ˙/q∞。式中，ρ為流體密度，Φ˙=?Φ/?t為速度勢變化率?？梢娨饓毫ψ兓囊蛩貫閮刹糠郑?-q/q∞為動壓變化所導(dǎo)致的壓力系數(shù)變化，其分布如圖7所示；-ρΦ˙/q∞為速度勢變化率導(dǎo)致的流場壓力變化。兩者綜合作用表現(xiàn)為物面壓力系數(shù)的分布，如圖8所示。

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

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

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

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

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

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

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

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

3 結(jié) 論

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

（1）不同狀態(tài)尾空泡流場穩(wěn)態(tài)數(shù)值計算表明，隨著空化數(shù)的提高，尾空泡壓力系數(shù)逐漸降低，導(dǎo)致臨近尾空泡的物面壓力系數(shù)下降，尾空泡對上游物面壓力影響范圍也相應(yīng)擴大；

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

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