水下垂直發(fā)射航行體尾空泡振蕩演化特性

權曉波， 尤天慶， 張晨星， 王凡瑜， 孔德才

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

0 引言

航行體水下垂直彈射出筒后，發(fā)射筒內燃氣會附著在航行體尾部形成尾空泡[1]。在水下垂直向上運動過程中，尾空泡呈現(xiàn)顯著的瞬態(tài)振蕩演化特性，對周圍流動產生明顯影響，是航行體水下垂直發(fā)射過程的主要多相流動現(xiàn)象[2-3]。

尾空泡發(fā)展演化會對水下發(fā)射流體動力產生重要影響，對其研究目前主要在壓力振蕩、演化形態(tài)以及氣體與水(簡稱氣水)摻混流態(tài)3方面進行。在壓力振蕩研究方面，鑒于尾空泡對航行體軸向彈道的嚴重影響，研究者基于氣泡動力方程對尾空泡壓力振蕩進行模擬并耦合航行體運動求解，獲得與試驗吻合較好的壓力振蕩規(guī)律[4-5]。在空泡演化形態(tài)方面，研究者基于彈射試驗獲得的尾空泡演化物理特征認識，基于空泡獨立膨脹原理或勢流理論，建立尾空泡形態(tài)計算模型，獲得了出筒速度等因素對空泡演化的影響規(guī)律[6-8]。此外基于Navier-Stokes方程的多相流模擬技術發(fā)展，使尾空泡形態(tài)演化模擬考慮的影響因素更加豐富[2,9]。在尾空泡氣水多相流態(tài)方面，研究多數聚焦于多相流態(tài)對水下噴氣射流影響，以及多相流的膨脹斷裂等復雜流動現(xiàn)象[10-11]。

水下發(fā)射水介質不可壓縮流動特性決定了尾空泡演化對流體動力的影響并不局限于航行體尾部。除上述尾空泡3方面影響外，航行體上游物面附近流動也應考慮尾空泡擾動。針對此問題研究，勢流算法可抓住不可壓流動主要特征，在假設無黏、無旋條件下，以較小計算資源消耗進行影響因素研究。同時文獻[12-15]表明，勢流算法在進行空泡流計算時，體現(xiàn)了較好的精度。

針對水下發(fā)射尾空泡演化對上游流動的影響問題，本文主要基于二維軸對稱勢流邊界元方法，開展了針對水下垂直發(fā)射航行體尾部附著空泡的數值計算研究。重點關注航行體物面壓力計算，通過與水洞試驗數據對比，驗證了數值計算方法的有效性。在此基礎上，針對尾空泡瞬態(tài)振蕩影響，分別從上游物面壓力、沾濕區(qū)阻力以及空泡發(fā)展3方面進行了數值計算研究。

1 數值計算方法

1.1 控制方程

假設空泡附著在物面擾流低壓區(qū)，在無限大流體域內振蕩。并假設流體無黏、無旋且不可壓，則Navier-Stokes方程可簡化為拉普拉斯方程。

(1)

(2)

式中:Φ為總速度勢，Φ=ux+φ，u為航行體運動速度(沿x軸負向)，x為位置坐標，φ為擾動速度勢；v為速度向量。

在物面，流體需滿足不可穿透條件。即壁面法向速度為0 m/s，

(3)

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

在空泡壁面，速度勢Φ需滿足伯努利方程：

(4)

式中：p0為水面壓力；ρ為水的密度；g為重力加速度；pc為空泡內壓力?？张荼诿孀鴺藊需滿足：

(5)

航行體表面壓力變化是本文研究的重點關注參數，流場中任意一點壓力p，可由伯努利方程計算得到：

(6)

式中：?Φ/?t直接采用速度勢時間差分進行求解是比較困難的。本文求解中利用?Φ/?t也滿足拉普拉斯方程的特點，采用與φ求解相同的數值方法。

1.2 數值算法

本文在數值方法上采用間接邊界元法求解，以求解擾動速度勢φ為例，介紹求解方法。設空間場點為p，物面和空泡壁面上分布的點源位置為q，其強度為λ(q)，則其在p點產生的速度勢和速度分別為

(7)

(8)

式中：G(p,q)為格林函數；S為面積；np為p點速度方向矢量；dSq為q點附近微元面積。對于拉普拉斯方程，三維情況下格林函數G(p,q)=-1/4πr，其中r=|p-q|.

數值計算中，航行體與空泡表面被劃分為N個單元，對航行體表面的第i個單元，考慮在邊界j(j=1，2，…，N)上源匯自身對速度的影響，擾動勢所滿足的邊界元積分控制方程如(9)式所示：

(9)

式中：Sj為第j個單元面積；u為航行體運動速度向量；nn為單元法向向量；λi為第i個單元的強度；λj為第j個單元的強度；dSqj為qj點附近微元面積。

空泡分離點為位于第d個單元，其擾動勢所滿足的邊界元積分控制方程如(10)式所示：

(10)

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

為節(jié)省計算資源提高計算效率，本文采用二維軸對稱計算模型，在航行體及空泡表面的邊界單元為圓環(huán)狀。與速度勢及速度相關的源匯均勻沿圓環(huán)分布。對于子午線長度為l圓環(huán)邊界單元，其對空間點p(xp,yp,0)的影響速度勢和速度由圓環(huán)上源匯q(xq，yqcosθ，yqsinθ)積分獲得，θ為圓柱坐標系下的方位角，影響速度勢和速度表達式分別為

(11)

(12)

(13)

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

為驗證數值計算代碼的正確性，在無黏、無旋理想流體條件下，計算了A和B兩種軸對稱旋成體零攻角物面壓力分布，并與水洞測壓試驗數據對比，如圖1所示。圖1中：Cp為壓力系數，Cp=(p-p0)/0.5ρu2；X為航行體壓力監(jiān)測點的位置坐標，D為航行體直徑；A為流線性外形，旋成體表面為全沾濕流動，對比表明數值計算與試驗數據吻合較好；B為錐柱組合外形，旋成體錐肩部附著有低壓空泡，數值計算結果較準確地描述了空泡低壓區(qū)范圍以及空泡末端回射壓力變化趨勢。

圖1 數值計算與試驗數據對比Fig.1 Comparison between calculated and experimental results

需要說明的是，空泡末端實際流動中存在較強的湍流現(xiàn)象，本文計算對此進行了簡化處理，導致數值計算的回射壓力量值與試驗結果存在一定差異，但這并不妨礙后續(xù)肩空泡回射壓力變化趨勢分析。

2 尾空泡振蕩影響分析

水下垂直發(fā)射航行體向水面運動過程中環(huán)境壓力逐漸減小，尾空泡發(fā)展演化以及航行體物面流動參數不斷變化。此過程流場的計算即涉及航行體相對固定的物面流場計算，又涉及空泡動態(tài)發(fā)展過程的計算?；谏鲜龌旌线吔鐥l件的邊界元計算方法，可求得航行體物面和空泡壁面的流場速度、加速度以及壓力等參數。在一個時間步內求得各物理量后，需通過邊界條件和空泡形態(tài)，對求解變量進行更新，本文采用4階龍格- 庫塔法進行數值求解。在時間步推進過程中，由于求解誤差積累，空泡壁面會產生鋸齒狀形態(tài)最終導致計算發(fā)散并終止，因此本文在空泡壁面動網格計算中進行了數值光順處理。

本文將從尾空泡振蕩過程模擬、對上游物面壓力影響、對沾濕區(qū)阻力影響以及對肩空泡影響這四方面進行論述，參見圖2.圖2中：坐標X/D≤0部分為尾空泡上游航行體物面沾濕區(qū)；坐標X/D>0部分為尾空泡輪廓，其內部壓力為pct；當航行體肩部物面壓力較低時，會附著有肩空泡，其內部壓力為pcn.

任何脫離了所在環(huán)境文化的能力都是不成熟的，新生代企業(yè)家的優(yōu)秀的勝任力，一定是要能在現(xiàn)有的中國情境下，仍然能夠出色發(fā)揮，而不是被動地等待環(huán)境來適應新生代企業(yè)家管理理念的勝任力。要融入中國情境，就必須先從多個角度來了解中國文化，符合大眾的價值取向。同時，筆者也試著從中西方科學文化的共同、共通之處為起點，對相異之處進行對比論證，為存在的問題做出合理的解釋。

圖2 尾空泡振蕩影響示意圖Fig.2 Sketch map of oscillation influence of tail cavity

2.1 尾空泡振蕩的模擬

為模擬航行體由一定深度向水面運動過程中尾空泡振蕩過程，首先進行一定深度條件下空泡流場穩(wěn)態(tài)計算，獲得σ0=(p0+ρgh-pct)/0.5ρu2條件下空泡流場，其中σ0為此穩(wěn)態(tài)計算獲得流場的初始空泡數，h為深度。以此穩(wěn)態(tài)計算結果作為后續(xù)瞬態(tài)計算初始條件，改變空泡內初始壓力，使空泡數為σ′0，引入尾空泡振蕩擾動，尾空泡壓力振蕩如圖3所示，其中t1、t2、t3、t4、t5為空泡收縮過程中等時間間隔點。由于在計算過程中航行體垂直向上運動，周圍環(huán)境壓力逐漸降低，因此尾空泡壓力在振蕩中逐漸降低。此外，本文數值算法對空泡壁面運動速度和空泡形態(tài)進行了平滑處理，在計算過程中引入了數值黏性，因此壓力振蕩體現(xiàn)出了衰減特征。

圖3 尾空泡內壓力振蕩數值計算結果Fig.3 Numerically calculated results of pressure oscillation of tail cavity

尾空泡振蕩過程典型時刻空泡形態(tài)如圖4所示。對應空泡瞬態(tài)收縮泡內壓力升高的過程，t1～t5等時間間隔，具體時刻如圖3中標注。由空泡振蕩形態(tài)變化可見，空泡壁面收縮或膨脹為變速過程。例如針對收縮過程，t1～t2時間段空泡壁面具有較大加速度，但收縮速度較小，整體表現(xiàn)t1和t2時刻空泡形態(tài)變化不大。隨著空泡壁面加速度逐漸累加，空泡壁面加速度收縮，空泡形態(tài)差異逐漸明顯。

圖4 空泡收縮過程形態(tài)變化Fig.4 Profile change of cavity during shrinking

2.2 尾空泡振蕩對上游物面壓力影響

尾空泡瞬態(tài)振蕩使周圍水介質運動加速度和速度產生顯著振蕩，進而對空泡上游航行體物面壓力產生明顯影響，如圖5所示。圖5中：T為壓力波動周期，ΔCp，max為各位置壓力系數變化量最大值。隨著航行體向水面運動，環(huán)境壓力逐漸降低，物面壓力逐漸下降，尾空泡振蕩使上游航行體物面壓力下降過程疊加了一定幅值的壓力波動。尾空泡起始位置為X/D=0，相應位置壓力系數變化曲線表明，測點位置越接近尾空泡壓力波動幅度越大。

圖5 尾空泡振蕩引起的沾濕區(qū)物面壓力波動Fig.5 Pressure oscillation on wet surface caused by tail cavity

為分析物面壓力波動規(guī)律，濾去壓力變化時程曲線中上升運動引起的靜壓變化，僅關注流體運動速度和加速度變化所引起的壓力系數變化量ΔCp=Cp-ρgh/0.5ρu2，采用ΔCp，max和T對壓力系數曲線進行無量綱化。取水下垂直發(fā)射彈射試驗航行體物面典型測點壓力，與數值計算結果對比，如圖6所示。二者在壓力變化規(guī)律上吻合較好，一方面體現(xiàn)在壓力隨時時間變化趨勢上，另一方面體現(xiàn)在壓力波動在空間上的衰減趨勢上。

圖6 尾空泡壓力振蕩計算與試驗數據對比Fig.6 Comparison between calculated and experimental results

為進一步分析尾空泡振蕩對航行體上游物面壓力影響在空間上的分布規(guī)律，提取了典型時刻物面壓力波動量值沿空間分布情況，如圖7所示。數值計算結果表明，尾空泡振蕩引起的物面壓力擾動以近似指數函數在空間分布，因此基于上述分析認為，當航行體長徑比較小時，尾空泡振蕩對航行體整體受力及上游流動會有較為明顯的影響。

圖7 不同初始壓力尾空泡收縮引起的物面壓力變化Fig.7 Pressure change of tail cavity during shrinking

2.3 尾空泡振蕩對航行體沾濕區(qū)阻力影響

計算長徑比L/D為6情況下，空泡振蕩條件下航行體垂直向水面運動過程受力情況，如圖8所示。其中，航行體阻力系數Cd由物面壓力積分求得

圖8 尾空泡振蕩對沾濕區(qū)阻力影響Fig.8 Effect of tail cavity oscillation on resistance in wet area

(14)

式中：Sw為航行體沾濕區(qū)域；Sr為參考面積；nx為微元法向向量x軸分量。

尾空泡振蕩引起的航行體上游物面壓力波動明顯地表現(xiàn)在阻力系數上，如圖8所示。其中無尾空泡振蕩影響情況，即在向水面運動過程中空泡形態(tài)保持不變的計算結果。由于向水面運動過程中，隨著環(huán)境壓力逐漸降低，尾空泡振蕩過程中形態(tài)逐漸增長，導致周圍流場加速度和速度變化區(qū)域增加。因此，相比尾空泡無變化工況，尾空泡振蕩對阻力系數影響，體現(xiàn)為振蕩中的緩慢變化。

2.4 尾空泡振蕩對肩空泡流場影響

尾空泡振蕩引起的上游物面壓力波動，不僅對航行體受力產生明顯影響，在上游存在空泡情況下，也會對上游空泡流場演化產生影響。針對此問題，本文數值計算獲得了一定泡內壓力條件下肩尾空泡穩(wěn)態(tài)流場，并以此穩(wěn)態(tài)計算結果為初始條件，通過降低尾空泡泡內壓力，引入流場擾動開始瞬態(tài)計算。計算一定時間段后，獲得肩尾空泡流場，如圖9所示。為方便對比分析，在計算過程中忽略尾空泡壁面速度，使尾空泡形態(tài)保持不變，進而獲得忽略尾空泡振蕩影響的計算結果。

由于在瞬態(tài)計算初始條件給定中降低了尾空泡內壓力，進而導致航行體周圍壓力顯著下降，肩空泡逐漸開始膨脹。由于尾空泡收縮使流場上游壓力提高，因此肩空泡發(fā)展相對緩慢。這種影響體現(xiàn)在受尾空泡影響下肩空泡尺寸相對較小，如圖9所示，同時也體現(xiàn)在肩空泡泡內壓力變化上，如圖10所示。在肩空泡內同樣氣體含量條件下，尾空泡收縮過程明顯使肩空泡泡內壓力提升。

圖9 尾空泡收縮對肩空泡發(fā)展的影響Fig.9 Effect of tail cavity shrinking on evolution of shoulder cavity

圖10 尾空泡收縮對肩空泡內部壓力影響Fig.10 Effect of tail cavity shrinking on pressure of shoulder cavity

在肩空泡末端，空泡壁面水流凹向壁面，使其損失徑向速度，形成回射壓力峰。尾空泡振蕩在改變肩空泡發(fā)展速度情況下，使航行體沾濕區(qū)壓力顯著提升，進而使肩空泡末端回射壓力峰值增大，如圖11所示。

圖11 尾空泡收縮對肩空泡流場壓力系數影響Fig.11 Effect of tail cavity shrinking on pressure distribution of shoulder cavity

3 結論

本文通過求解二維軸對稱航行體水下垂直運動過程空泡流場瞬態(tài)演化，分析了尾空泡振蕩對航行體上游物面壓力、受力以及空泡演化的影響。得到主要結論如下：

1) 通過與水洞測壓數據比對，驗證了本文數值方法對航行體上游物面壓力計算的準確性，同時航行體附著尾空泡的混合邊界元計算可模擬尾空泡泡內壓力以及空泡形態(tài)的瞬態(tài)振蕩過程。

2) 在尾空泡振蕩對上游物面壓力擾動計算方面，計算結果取得與彈射試驗數據相一致的變化規(guī)律；同時計算分析表明，隨著關注位置遠離尾空泡，壓力擾動呈指數函數衰減。

3) 針對航行體上游物面整體受力，隨著航行體垂直向水面運動環(huán)境壓力逐漸降低，尾空泡振蕩過程中長度逐漸增加，其振蕩對阻力系數影響，體現(xiàn)為振蕩中的緩慢變化。

4) 針對航行體上游肩空泡流場演化，肩尾空泡耦合計算表明，尾空泡收縮使流場上游肩部空泡發(fā)展速度減緩，同時使上游肩空泡泡內壓力和末端回射壓力峰值偏大。