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

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

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

0 引 言

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

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

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

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

1 數(shù)值計(jì)算方法

1.1 控制方程

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

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

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

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

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

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

1.2 數(shù)值算法

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

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

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

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

式（8）和（9）即為邊界元計(jì)算的基本方程，通過(guò)這兩個(gè)方程即可求得物面及空泡壁面的λ分布，進(jìn)而確定流場(chǎng)中速度勢(shì)、速度及壓力。

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

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

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

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

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

1.3 尾空泡閉合模式

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

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

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

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

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

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

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

2.2 瞬態(tài)計(jì)算

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

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

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

圖7 收縮過(guò)程不同時(shí)刻物面及空泡表面動(dòng)壓分布對(duì)比Fig.7 Comparison of dynamic pressure distribution

針對(duì)尾空泡收縮瞬態(tài)變化過(guò)程，本文流場(chǎng)壓力計(jì)算采用非穩(wěn)態(tài)伯努利方程，如式（5）所示。基于此，壓力系數(shù)為Cp=1-q/q∞-ρΦ˙/q∞。式中，ρ為流體密度，Φ˙=?Φ/?t為速度勢(shì)變化率?？梢?jiàn)引起壓力變化的因素為兩部分：1-q/q∞為動(dòng)壓變化所導(dǎo)致的壓力系數(shù)變化，其分布如圖7所示；-ρΦ˙/q∞為速度勢(shì)變化率導(dǎo)致的流場(chǎng)壓力變化。兩者綜合作用表現(xiàn)為物面壓力系數(shù)的分布，如圖8所示。

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

圖8 收縮過(guò)程壓力分布對(duì)比Fig.8 Comparison of static pressure distribution

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

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

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

圖10 不同初始?jí)毫ξ部张菔湛s引起的物面壓力變化Fig.10 Pressure change during the tail cavity shrinking process with different initialpressures

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

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

3 結(jié) 論

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

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

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

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