李智生

（91550部隊 大連 116023）

1 引言

對于采用水中點火的垂直發(fā)射航行體而言，通過及時對航行體的水中彈道和姿態(tài)進行控制，可提高航行體在水中運動的穩(wěn)定性，并有助于獲得有利的出水姿態(tài)，從而增強航行體出水后的姿態(tài)控制能力［1］。通常，水下發(fā)射航行體在其尾部離開發(fā)射筒口一定距離處點火進行發(fā)射。如果航行體尾部與發(fā)射筒口距離較近處點火，則發(fā)動機可以利用從發(fā)射筒內(nèi)溢出并附著在航行體尾部的燃氣泡，作為發(fā)動機噴流建立初期燃氣的受納空間，減輕發(fā)動機直接在水中點火所造成的沖擊，提高發(fā)動機工作安全性，但是由于水下點火處距離發(fā)射平臺較近，發(fā)動機高速噴流會對筒蓋及周圍發(fā)射平臺結(jié)構(gòu)產(chǎn)生破壞，影響后續(xù)發(fā)射［2］。如果航行體遠離發(fā)射筒口處（簡稱筒口）點火，附著在航行體尾部的燃氣泡體積變小，同時燃氣泡內(nèi)滲人水，發(fā)動機若直接在此環(huán)境中點火，由于水的巨大慣性，則發(fā)動機噴流會受到水的阻礙，導致噴管內(nèi)壓強過高，從而威脅發(fā)動機工作安全性［3～4］。顯然，在決定采用水中點火方案時，必須明確水中點火對航行體載荷的影響，并充分估計水中點火對發(fā)射筒及發(fā)射平臺所造成的威脅［5～6］。

本文在構(gòu)建垂直發(fā)射航行體水中點火數(shù)值仿真模型的基礎(chǔ)上，研究航行體水中點火過程中尾部流場變化情況，分析不同點火深度下發(fā)動機射流流場對發(fā)射平臺壁面壓力和溫度的影響，同時對水中點火燃氣后效進行估計。

2 控制方程

建立水下發(fā)射過程中的航行體動力模型，通過實時求解航行體受力，通過動力學方程求得加速度和角加速度，然后進行數(shù)值積分，獲得航行體當前速度、位移及姿態(tài)等。航行體受力如圖1所示，圖1中oxyz為發(fā)射坐標系［7］，相對于地球靜止，ox1y1z1為航行體坐標系，固聯(lián)于航行體；v為航行體運動速度；F為航行體所受的浮力；G為自身重力；f為航行體軸向阻力；ox1為航行體軸向，ox1與v的夾角為α；oy1與v的夾角為θ，T為發(fā)動機產(chǎn)生的總推力；?為發(fā)動機噴管為糾正航行體出水姿態(tài)而偏轉(zhuǎn)的擺角。

圖1 發(fā)射坐標系

航行體在水下發(fā)射過程中，首先利用高壓燃氣將航行體彈射出筒，航行體尾部離開發(fā)射筒后，筒內(nèi)燃氣進入外界環(huán)境，與周圍的水汽混合并附著在尾部形成燃氣泡。發(fā)動機接收到點火指令后，航行體運動的發(fā)射動力主要由兩部分組成：彈射力和發(fā)動機提供的軸向推力。

航行體在水下發(fā)射過程中，其航行體動力模型為［8］

其中：

式中：t為航行體運動時間，F(xiàn)'為發(fā)射筒高壓燃氣彈射力，T'為發(fā)動機點火后航行體受到的軸向推力，m為航行體質(zhì)量，s為航行體的迎流面積，C為阻力系數(shù)，ρ為海水密度，λ為附加質(zhì)量；α和θ值利用航行體加速度計測量的各個方向的加速度值進行解算求取。式（1）中的航行體運動參數(shù)值為后續(xù)計算燃氣后效參數(shù)提供動態(tài)輸入。

3 模型建立

3.1 計算區(qū)域及網(wǎng)格劃分

CFD計算模型具體采用二維軸對稱的無限長圓柱體模型。CFD求解器計算域如圖2所示。網(wǎng)格劃分的原則是，全部流場區(qū)域采用結(jié)構(gòu)化網(wǎng)格。該方式能夠消除計算結(jié)果對網(wǎng)格的依賴性［9～10］。在劃分噴管時，在軸向方向、噴管及尾部附近處的網(wǎng)格大小接近一致。外場劃分梯度網(wǎng)格，靠近航行體位置網(wǎng)格細密，遠場較為稀疏。按照以上要求劃分網(wǎng)格，得到航行體尾部局部網(wǎng)格如圖3所示。

圖2 CFD求解器計算域

圖3 航行體尾部局部網(wǎng)格

將式（1）和式（2）所示的航行體動力模型輸入到CFD求解器中，利用CFD求解器對航行體水下運動參數(shù)進行求解。在得到航行體運動參數(shù)后對網(wǎng)格節(jié)點參數(shù)進行更新，進而對航行體能量方程及控制方程進行數(shù)值離散求解，最終得出航行體近體流場壓力及溫度分布參數(shù)。由于模型的對稱性，計算域模型采取二維軸對稱模型。在建模時不考慮航行體尾部產(chǎn)生的空化現(xiàn)象，且由于研究的是尾部流場變化，頭部對仿真流場的影響較小，因此將計算模型簡化為尾部只有噴管擴張段的無限長圓柱體［11］。

3.2 計算方法

采用Mixture多相流模型對流場進行求解，其控制方程基本形式如下［12］：

式中，ρm為混合物密度：

式中，αk為第k相的體積分數(shù)。利用Mixture多相流模型進行求解時，分別設(shè)置了主相（Primary phase）和副相（Secondary phase），所有相體積分數(shù)之和為1，副相可為多個，對于副相p的體積分數(shù)求解方法為

動量方程是牛頓第二定律在流體動力學中的體現(xiàn)。對于Mixture多相流模型的動量方程，可通過將所有相各自的動量方程相加的方法獲得，其表示式為

式中，n為相的序號，為體積力，μm為混合物相的粘性，其表達式為

能量方程是能量守恒定律在流體動力學中的體現(xiàn)，Mixture多相流模型中混合物的能量方程為

式中，keff為有效傳導系數(shù)，表達式為

其中，kt為湍動熱傳導系數(shù)，kk為第k相的湍動能，其取值取決于所選擇的湍流模型。式（13）右端項中第一項代表了傳導引起的能量遷移，SE為其余的體積熱源項。

對于可壓縮流動，Ek的表達式為

對于不可壓縮流動，表達式為

其中，hk為第k相的焓。

上述偏微分方程組，包括式（3）～式（16），通過網(wǎng)格節(jié)點在計算域流場內(nèi)進行離散化求解，具體步驟如下：

1）在計算過程的每個時間步，利用CFD方法聯(lián)合求解Mixture模型控制方程、動量方程、能量方程和航行體動力模型，獲得流場和航行體運動參數(shù)，通過更新網(wǎng)格將二者進行耦合，計算完畢后轉(zhuǎn)入下一時間步；

2）當?shù)谝徊降玫降乃俣仍诰植坎粷M足連續(xù)方程時，從連續(xù)性方程和線化動量方程推導出壓力校正的泊松方程，然后解出壓力校正方程，獲取壓力和速度場；

3）用上一個時間步長更新的、除了壓力速度和溫度外的其它變量值，解出湍流、能量和輻射等標量；

4）當存在相間耦合時，用離散相軌跡計算來更新連續(xù)相的源項；

5）根據(jù)計算的殘差曲線檢查設(shè)定方程的收斂性，當所有變量的殘差值都降到10-3時，就認為計算收斂，即完成了對壓力和溫度的離散求解。

3.3 計算設(shè)置

表1給出了一組計算初始參數(shù)設(shè)置。

表1 主要參數(shù)設(shè)置

在水域邊界設(shè)置中，初始壓力輸入數(shù)值通過重力梯度法獲得。重力梯度法就是根據(jù)遠場邊界網(wǎng)格節(jié)點坐標和航行體發(fā)射水深來計算原場邊界壓力輸入的方法。設(shè)標準大氣壓為P0=101325Pa，重力加速度g=9.81m/s2，海水密度ρ=1.02kg/m3，則初始遠場壓力總壓，簡稱為總壓，為P=P0+ρgh，其中，h為壓力邊界網(wǎng)格節(jié)點高度與航行體發(fā)射深度的差值。總溫及回流總溫初始值由發(fā)射筒內(nèi)溫度傳感器讀取。

在CFD模型設(shè)置中，湍流模型采用RNGk-ε，可以達到較好的收斂性；由于流場變量在壁面附近存在很大梯度的流動，因此壁面函數(shù)采用非平衡壁面函數(shù)法；流場初始化過程中，先對整流場按表1進行參數(shù)初始化，而后標記出全部為水的區(qū)域并對此區(qū)域按表1中初始項給出的參數(shù)修正，再迭代計算；在計算步長選擇上，首先保證計算精度，其次要保證單時間步迭代計算的收斂性，最后要考慮計算的經(jīng)濟性。

在對航行體近體流場進行離散求解后，即可得出流場屬性的關(guān)鍵參數(shù)（如壓力、密度和溫度等）。

4 仿真計算與結(jié)果分析

不同位置點火時，發(fā)動機噴管所黏附的燃氣泡體積不同，如果燃氣泡內(nèi)點火，則可以起到緩沖作用，從而減弱對航行體、發(fā)射平臺的沖擊載荷。仿真計算中，通過設(shè)置不同的點火距離，分析燃氣后效中壓力和溫度對發(fā)射平臺壁面的影響。

采用如表1所示的參數(shù)來設(shè)置CFD計算模型。假設(shè)發(fā)動機距筒口點火的距離分別為2.2、4.3、6.5（數(shù)值進行了無量綱化處理），對應(yīng)時間為0.4、0.6、0.8（數(shù)值進行了無量綱化處理）。在點火時刻，F(xiàn)luent設(shè)置中將喉部的邊界wall改為壓力入口（PressureInlet），通過Fluent自定義函數(shù)（UDF）定義航行體質(zhì)量和轉(zhuǎn)動慣量等參數(shù)，通過網(wǎng)格節(jié)點離散求解航行體水動力及壓力等參數(shù)，湍流模型選用增強RNGk-ε模型，計算迭代步長為1E-5。

通過自定義函數(shù)UDF設(shè)置壓力監(jiān)測點位置，得到的發(fā)射平臺壁面的壓力與溫度的變化規(guī)律曲線分別如圖4和圖5所示。從圖4中可以看出，發(fā)射平臺壁面各個測點的壓力脈動變化規(guī)律基本一致，且離筒口中心越遠，脈動幅值越小，當t>0.7后，壓力值基本維持在當?shù)厮顗毫Φ?.55倍；在點火后，測點Pb1的壓力值出現(xiàn)3次壓力峰值，這是由于Pb1是發(fā)射筒與發(fā)射平臺壁面的交界點，發(fā)動機燃氣射流與發(fā)射筒內(nèi)高壓燃氣同時對Pb1產(chǎn)生作用。圖5所示，發(fā)射筒中心處各個測點的溫度變化規(guī)律趨于一致，溫度分布符合熱傳導規(guī)律，有明顯的滯后，從筒口開始自上而下，溫度降低；在發(fā)動機點火發(fā)射開始后，測點P1的溫度瞬間升至3000K，其他測點溫度也逐漸升高；其后又出現(xiàn)兩次溫度波動，其峰值要遠低于第一次，期間出現(xiàn)溫度降低的原因是：取兩個典型時刻t=0.231s和t=0.471s，在密度云圖6中這兩個時刻的溫度云圖和流線圖也可以看出，此刻燃氣泡幾乎處于閉合狀態(tài)，燃氣溫度無法傳遞。

圖4 發(fā)射平臺壁面不同測點壓力變化圖

圖5 發(fā)射平臺壁面不同測點溫度變化圖

5 結(jié)語

本文基于Mixture多相流模型、動網(wǎng)格技術(shù)及三維軸對稱模型，建立了航行體水中點火數(shù)值仿真模型，仿真計算得到了不同點火深度條件下燃氣后效中壓力和溫度的變化特性數(shù)據(jù)，結(jié)果分析表明，發(fā)射平臺壁面各個測點的壓力脈動變化規(guī)律基本一致，且離筒口中心越遠，脈動幅值越?。粶囟确植挤蠠醾鲗б?guī)律，有明顯的滯后，從筒口開始自上而下，溫度降低。研究結(jié)論能夠為航行體水中點火時機的選擇提供決策依據(jù)，可為發(fā)射筒筒口流場分析提供數(shù)據(jù)支撐。

圖6 燃氣泡的非定常流狀態(tài)