一種火箭發射場坪流場的環形線匯描述方法

謝 建，權 輝，謝 政

(火箭軍工程大學兵器發射理論與技術軍隊重點實驗室，西安 710025)

0 引 言

火箭發射時，發射場坪上的流場發生劇烈變化[1]，會對設施設備產生巨大影響：一是發射過程中燃氣流對場坪上的建筑設施及設備產生沖擊；二是場坪上的小型設備、顆粒雜物等被氣流卷吸進入發射井，對火箭及井內設施設備產生危害。所以，發射場坪上流場的流動規律是值得重點關注的問題之一。

對于火箭發射時的流場，王飛采用流體計算軟件對W型地下井發射過程進行了數值模擬，分析了壓力和溫度的變化規律[1]。李軍采用數值模擬方法研究了固體火箭燃氣流場形成和發展的機理[2]。馬溢清采用氣液兩相流模型對火箭發動機尾焰注水降噪效果進行了研究并做了縮比試驗[3]。傅德彬應用動網格技術對火箭發射過程中的燃氣流場進行了數值模擬，分析了燃氣流場對火箭氣動力的影響[4-5]。姜毅采用一維等熵流方法和有限體積法分析了固體火箭發動機尾焰流場的流場結構和燃氣組分分布[6-7]。楊宏偉對火箭垂直發射時的流場進行了數值模擬，得到了流場變化規律[8]。謝政等采用數值方法研究燃氣流的二次燃燒現象以及點火壓力脈沖形成機理和影響因子[9-10]。單時卓采用數值方法對艦載導彈垂直發射過程中的燃氣流場進行了分析[11]。周笑飛對井下發射燃氣射流的數值模擬方法進行了深入研究[12]。這些研究多數針對某一具體問題進行研究，其結果對于后續相關研究具有重要指導意義。

隨著火箭技術的不斷進步和火箭運用方式的不斷發展，發射場坪流場的研究變得尤為迫切，必須進行深入分析，建立一套行之有效的分析方法。本文基于勢流法提出環形線匯模型對發射場坪上的流場進行描述，分析發射場坪流場的速度變化規律，并與FLUENT仿真結果進行對比，驗證方法的可行性，為后續發射井相關問題的研究提供新的思路。

1 環形線匯方法

1.1 工程背景介紹

某型火箭發射井構型如圖1所示，火箭由發射臺支撐豎立在發射井中，點火后燃氣尾焰由排焰道排除井外，火箭在燃氣射流的反作用推力下起飛出井。從火箭點火到出井過程中，發射場區的流場經歷爆震波階段、引流階段和失穩階段等三個明顯的階段變化[1][9]。

在引流階段，火箭燃氣射流在井下狹窄空間內產生引射作用，整個發射井相當于大型的引射器，在引射器作用下發射場坪上形成強烈的引流流場。引流流場中流體介質主要為場坪上的空氣，溫度為常溫，流動過程中流速低于聲速，相對于高溫高速燃氣射流，可以作為不可壓縮流體進行處理。由于空氣邊界層很薄，而邊界層外流體黏性作用可以忽略[13]，為簡化分析，本文采用理想流體進行計算。

為方便研究，對發射場坪流場進行適當簡化。如圖2和圖3所示，將火箭與井壁間隙(以下簡稱箭井間隙，并用d表示)簡化為環形線匯，半徑取為火箭半徑與發射井半徑的平均值，于是發射井出口處的流場可以用環形線匯流場來描述。

1.2 環形線匯流場速度勢推導

設環形線匯強度為Q1，半徑為R，如圖3所示，以發射場坪中心為原點，場坪所在平面為XOY平面，沿發射井軸向指向火箭頭部方向為Z軸正向，建立右手O-XYZ笛卡爾坐標系，場坪上空中一點M(x,y,z)受到井口一點N1(Rcosθ,Rsinθ, 0)作用的速度勢[14]為

(1)

沿井口圓環對φ1積分可得M1受到整個環形線匯作用的速度勢為

(2)

式中：γ為OM1在XOY平面的投影線OM1＇與X軸正向夾角。

(3)

式中：K(ζ)為第一類橢圓積分[15-16]，ζ為模數。

設井壁與火箭間隙流場為均勻流動，速度為v，則Q1可以表示為

(4)

于是m1式可以化為

m1=v(R2-R1)

(5)

令d=R2-R1，將(5)式代入(3)式可得

(6)

1.3 環形線匯流場速度的推導

根據勢流理論

(7)

式中：u1、v1、w1分別代表發射場坪X、Y、Z三個方向的流速(X、Y、Z方向如圖3所示)。

將(6)式代入(7)式，結合第一類橢圓積分的求導公式[18]可得

(8)

式中：E(ζ)為第二類橢圓積分，ζ為其模數。

同理可得

(9)

(10)

1.4 環形線匯流場理論計算

根據公式(8)(9)(10)，取v=-100，d=0.5，y=0，R=1.55，z=0.1，x∈(-12, 12)得到X方向速度vx如圖4所示。從圖上可以看出，發射場坪地面附近水平速度在靠近井口時急劇增大，且在井口半徑以內沿半徑增大方向，在井口半徑以外沿半徑減小方向，這顯示出流體向井的環形出口處聚集。

取v=-100，d=0.5，y=0，x=R=1.55，z∈(0.1, 12)得到Z方向速度vz如圖5所示。從圖上可以看出，在箭井間隙附近流體豎直方向速度沿Z軸負向急劇增大，隨著高度增加，速度逐漸趨近于0。

由ζ表達式可知，環形線匯模型的奇點位置為場坪地面上半徑為R的圓，在流場奇點附近速度趨近于∞，與實際流場有很大差距，在圖4、圖5上可以直觀看到這一現象，所以奇點附近的取值不宜作為實際流場的估算數據。

2 FLUENT仿真模型的建立

2.1 計算模型

本文以FLUENT仿真結果對比分析勢流理論計算結果的偏差。為簡化分析計算，采用圓柱軸對稱模型，選取的幾何模型如圖6所示，軸對稱處理的平面模型如圖7所示。為估計火箭半徑(r)以及箭井間隙(d)大小對環形線匯模擬發射場坪流場的準確度的影響，r值分別取為1.0 m、1.5 m、2.0 m、2.5 m、3.0 m，d值分別取為0.1 m、0.3 m、0.5 m、0.7 m、0.9 m，建立相應的幾何模型。為幾何模型劃分四邊形網格，以r=1.0 m，d=0.5 m為例得到的網格模型如圖8所示。

邊界條件和監測線的設置如圖7所示。在軟件界面中選擇圓柱軸對稱、穩態模型，k-ε湍流模式[17]，以常溫常壓全局初始化，入口速度v分別選擇-60 m/s、-100 m/s、-200 m/s、-300 m/s進行計算，得到監測線上的速度數據。

2.2 網格無關性驗證

在使用大型仿真軟件進行流場計算中，先將流場區域劃分成微小有限體積，再將控制微分方程組離散化為差分方程組，使用迭代方法求解，在此過程中會產生截斷誤差和舍入誤差，截斷誤差隨著網格劃分的加密逐漸縮小，而舍入誤差會隨著計算迭代次數的增多而增大。在實際計算機仿真中，計算結果的精度會隨著網格密度的增加而提高，而且提高量會相對越來越小，但是所需要的計算資源卻會急劇增加。為使用有限的計算資源達到可靠的計算精度，需要進行網格無關性驗證，即確定最小的網格數以取得所要求的計算精度。

本文中以r=1.0 m，d=0.5 m時幾何模型進行網格無關性驗證，得到適用于本項研究的網格密度。分別對網格數為19641、31376、42455、56204的四套網格進行計算，得到徑向監測線上的速度幅值如圖9所示。從圖上可以看出，網格數為19641的網格模型其監測線上速度幅值的計算結果與其他三種網格模型差距明顯，最大值為88.4347，而其余三種分別為81.2527、81.5286、82.6366。為保證計算精度，同時有效減少計算量，本文選取網格數為31376的網格模型的網格密度作為網格劃分的方案。

3 理論計算結果與仿真結果的對比

皮爾遜相關系數[18]廣泛用于度量兩個變量之間的線性相關程度，其值介于1和-1之間，趨近于1表明正向線性相關性趨強，趨近于-1表明負向線性相關趨強，趨近于0表明線性相關性趨弱。本文中用兩種方法計算同一個物理量——速度，數據結果之間沒有明顯平移和比例差距(見圖10～12)，所以僅用相關系數表明兩種結果的接近程度。

3.1 不同速度入口條件下環形線匯模型與仿真結果的差異

取火箭半徑r=1.0 m，分別取d=0.1 m、0.3 m，將入口速度v分別設為-60 m/s、-100 m/s、-200 m/s、-300 m/s，得到軟件仿真結果。分別提取徑向監測線上的徑向速度數據和軸向監測線上的軸向速度數據與環形線匯理論計算結果進行對比，計算兩種方法得到曲線的相關系數如表1和表2所示，表中ρr表示徑向速度相關系數，ρa表示軸向速度相關系數。從表中可以看出，d=0.1 m時，理論計算結果與軟件仿真結果吻合得非常好；d=0.3 m時，理論計算結果與軟件仿真結果的差距增大；不同入口速度對理論計算和軟件仿真結果的相關系數沒有顯著影響。

為直觀地比較兩種方法計算結果的差距，取d=0.3 m時兩種方法得到的徑向監測線上徑向速度數據畫出曲線如圖10所示，圖中水平坐標軸是徑向監測線上位置，單位是m，縱坐標是徑向速度，單位是m/s，實線是FLUENT仿真結果，虛線是公式計算結果。從圖上可以看出，雖然速度絕對值變大，但是兩種方法得到的曲線差值變化趨勢基本一致。

表1 d=0.1 m時不同入口速度條件下兩種方法計算結果的相關系數Table 1 Correlations of the results from the two methods under different inlet velocities when d=0.1 m

表2 d=0.3 m時不同入口速度條件下兩種方法計算結果的相關系數Table 2 Correlations of the results from the two methods under different inlet velocities when d=0.3 m

3.2 不同r值、d值條件下環形線匯模型計算結果與仿真結果的差異

取入口速度v為-100 m/s，用軟件計算不同r值、d值的流場模型得到25組計算結果。從計算結果中分別提取軸向監測線上的軸向速度數據(記為vaa)，計算其與環形線匯理論計算結果的相關系數如表3所示，分別提取徑向監測線上的徑向速度數據(記為vrr)，計算其與環形線匯理論計算結果的相關系數如表4所示，表中ρ0.1、ρ0.3、ρ0.5、ρ0.7、ρ0.9分別表示d=0.1 m、0.3 m、0.5 m、0.7 m、0.9 m時的相關系數。從表中數據可以看出：r值的變化對兩種方法計算結果的相關系數沒有顯著影響；軸向監測線上的軸向速度數據的相關性明顯好于徑向；d=0.1 m時，理論計算結果和軟件仿真結果符合得較好；隨著d值的增大，理論計算與仿真結果的差距逐漸增大。

表3 兩種方法計算得到的“vaa”的相關系數Table 3 Correlations of “vaa” from the two methods

表4 兩種方法計算得到的“vrr”的相關系數Table 4 Correlations of “vrr” from the two methods

為了直觀地比較兩種方法計算結果的差距，取d=0.5 m時不同r值下兩種方法得到的徑向監測線上徑向速度數據畫出曲線如圖11所示(圖例意義同圖10)，取r=1.0 m時不同d值下兩種方法得到的徑向監測線上徑向速度數據畫出曲線如圖12所示(圖例意義同圖10)。從圖11可以看出，理論計算值隨著r值增大而減小，但是兩種方法得到的曲線差值變化趨勢基本一致。從圖12可以看出，隨著d值增大，速度絕對值增大，同時兩種曲線差距也越來越大。

4 結 論

本文采用勢流法對發射場坪上的流場進行研究，推導了環形線匯流場速度勢和速度的計算公式，在FLUENT軟件中建立發射場坪流場計算模型并進行求解，在不同r值、d值和不同速度入口條件下分別對環形線匯與軟件仿真兩種方法流場計算結果的差異進行了分析，得到了以下結論：

(1)不同火箭半徑和流場引流入口速度對理論計算和軟件仿真結果的誤差沒有顯著影響；

(2)兩種方法得到的軸向監測線上軸向速度數據的相關性明顯好于徑向監測線上的徑向速度；

(3)環形線匯模型在d值較小時與軟件仿真得到的流場速度較為一致，隨著d值得增大，兩種方法得到的流場速度差距逐漸增大，當d>0.5 m時，徑向監測線上的徑向速度數據相關系數在85%以下；

(4)環形線匯模型能夠反映發射場坪流場速度的變化趨勢，在d值較小且誤差允許時可以作為發射場坪流場的簡化描述。