李智生

（91550部隊(duì) 大連 116023）

1 引言

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

本文在構(gòu)建垂直發(fā)射航行體水中點(diǎn)火數(shù)值仿真模型的基礎(chǔ)上，研究航行體水中點(diǎn)火過程中尾部流場變化情況，分析不同點(diǎn)火深度下發(fā)動(dòng)機(jī)射流流場對(duì)發(fā)射平臺(tái)壁面壓力和溫度的影響，同時(shí)對(duì)水中點(diǎn)火燃?xì)夂笮нM(jìn)行估計(jì)。

2 控制方程

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

圖1 發(fā)射坐標(biāo)系

航行體在水下發(fā)射過程中，首先利用高壓燃?xì)鈱⒑叫畜w彈射出筒，航行體尾部離開發(fā)射筒后，筒內(nèi)燃?xì)膺M(jìn)入外界環(huán)境，與周圍的水汽混合并附著在尾部形成燃?xì)馀荨０l(fā)動(dòng)機(jī)接收到點(diǎn)火指令后，航行體運(yùn)動(dòng)的發(fā)射動(dòng)力主要由兩部分組成：彈射力和發(fā)動(dòng)機(jī)提供的軸向推力。

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

其中：

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

3 模型建立

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

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

圖2 CFD求解器計(jì)算域

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

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

3.2 計(jì)算方法

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

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

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

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

式中，n為相的序號(hào)，為體積力，μm為混合物相的粘性，其表達(dá)式為

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

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

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

對(duì)于可壓縮流動(dòng)，Ek的表達(dá)式為

對(duì)于不可壓縮流動(dòng)，表達(dá)式為

其中，hk為第k相的焓。

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

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

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

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

4）當(dāng)存在相間耦合時(shí)，用離散相軌跡計(jì)算來更新連續(xù)相的源項(xiàng)；

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

3.3 計(jì)算設(shè)置

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

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

在水域邊界設(shè)置中，初始?jí)毫斎霐?shù)值通過重力梯度法獲得。重力梯度法就是根據(jù)遠(yuǎn)場邊界網(wǎng)格節(jié)點(diǎn)坐標(biāo)和航行體發(fā)射水深來計(jì)算原場邊界壓力輸入的方法。設(shè)標(biāo)準(zhǔn)大氣壓為P0=101325Pa，重力加速度g=9.81m/s2，海水密度ρ=1.02kg/m3，則初始遠(yuǎn)場壓力總壓，簡稱為總壓，為P=P0+ρgh，其中，h為壓力邊界網(wǎng)格節(jié)點(diǎn)高度與航行體發(fā)射深度的差值?？倻丶盎亓骺倻爻跏贾涤砂l(fā)射筒內(nèi)溫度傳感器讀取。

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

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

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

不同位置點(diǎn)火時(shí)，發(fā)動(dòng)機(jī)噴管所黏附的燃?xì)馀蒹w積不同，如果燃?xì)馀輧?nèi)點(diǎn)火，則可以起到緩沖作用，從而減弱對(duì)航行體、發(fā)射平臺(tái)的沖擊載荷。仿真計(jì)算中，通過設(shè)置不同的點(diǎn)火距離，分析燃?xì)夂笮е袎毫蜏囟葘?duì)發(fā)射平臺(tái)壁面的影響。

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

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

圖4 發(fā)射平臺(tái)壁面不同測點(diǎn)壓力變化圖

圖5 發(fā)射平臺(tái)壁面不同測點(diǎn)溫度變化圖

5 結(jié)語

本文基于Mixture多相流模型、動(dòng)網(wǎng)格技術(shù)及三維軸對(duì)稱模型，建立了航行體水中點(diǎn)火數(shù)值仿真模型，仿真計(jì)算得到了不同點(diǎn)火深度條件下燃?xì)夂笮е袎毫蜏囟鹊淖兓匦詳?shù)據(jù)，結(jié)果分析表明，發(fā)射平臺(tái)壁面各個(gè)測點(diǎn)的壓力脈動(dòng)變化規(guī)律基本一致，且離筒口中心越遠(yuǎn)，脈動(dòng)幅值越??；溫度分布符合熱傳導(dǎo)規(guī)律，有明顯的滯后，從筒口開始自上而下，溫度降低。研究結(jié)論能夠?yàn)楹叫畜w水中點(diǎn)火時(shí)機(jī)的選擇提供決策依據(jù)，可為發(fā)射筒筒口流場分析提供數(shù)據(jù)支撐。

圖6 燃?xì)馀莸姆嵌ǔＡ鳡顟B(tài)