發射井口引射流場理論計算與仿真①

謝 建，權 輝，謝 政，韓小霞

(1.火箭軍工程大學 兵器發射理論與技術軍隊重點實驗室，西安 710025；2.火箭軍工程設計研究院，北京 100111)

0 引言

火箭發射時，發動機在井下狹窄空間產生高溫高速燃氣射流，經排焰道排出井外，在某個時段內會形成引射現象[1]。引射作用在發射場坪上井口附近形成引射匯聚流場(以下簡稱引射流場)。場坪上小型設備、顆粒雜物等在引射流場作用下，可能被氣流卷吸進入發射井，對火箭及井內設施設備產生危害[2-3]。因此，發射井口引射流場的流動規律是值得重點關注的問題。

對于火箭發射時的流場，以往的研究主要關注點在火箭發動機燃氣射流流場。從單純針對火箭/導彈燃氣流場自由射流規律[4-5]，慢慢擴展到關注某一具體發射裝置中燃氣射流的流動規律，例如地下發射井流場規律的研究[6]，艦載導彈垂直發射裝置燃氣流場的研究[7-8]，同心筒中燃氣流場的研究[9]，雙噴管發動機燃氣流場的研究[10]等。近年來，又有學者關注到燃氣射流的一些細節現象，如固體火箭發動機燃氣射流近場形成和發展的機理[11]，燃氣流場的二次燃燒現象以及點火壓力脈沖形成機理和影響因子[12-13]。這些研究均采用了CFD方法，部分還進行了大量實驗驗證，具有較高的準確度，對于火箭發射技術的發展具有重要指導意義，但當前專門針對發射場坪流場的研究卻鮮見報道。

隨著火箭發射技術的不斷發展，發射井優化設計成為一個重要課題，發射安全以及發射井與周邊環境的相互作用是發射井優化所應考慮的重要方面，發射場坪流場規律是其中的重要環節。以往以CFD為基礎的方法雖然能獲得較高的計算準確度，但在進行發射井設計時，卻要反復建模，無法全面分析某個因素變化范圍內的所有情況，而且消耗大量計算資源和時間，達不到最優設計的效果。因此，在最初設計階段應用近似理論分析方法，在優化設計階段確定具體參數之后，引入CFD方法進行參數修正，最后進行試驗驗證，是一種合理的設計思路。

本文前期針對發射井口引射流場做了部分理論研究，提出了環形線匯方法，以勢流法為基礎，根據點匯模型積分得到發射井口引射流場的描述公式[14]。該方法在距離井口較遠時準確度較高，但由于公式的奇點效應，在井口附近的計算結果與實際偏差較大。本文針對前期簡單環形線匯模型的不足，提出改進措施，克服奇點的影響，提高計算的準確度，為發射場坪流場分析提供簡便可行且準確度較高的方法，進而為發射井優化設計提供一定的理論支撐。

1 改進的環形線匯方法

1.1 工程背景介紹

引射流場中流體介質主要為場坪上的空氣，溫度為常溫，流動過程中流速低于聲速，相對于高溫高速燃氣射流，可作為不可壓縮流體進行處理。忽略空氣黏性和重力，引射流場可用勢流法[15]進行計算。某型發射井結構圖如圖1所示。

(b) Top view

如圖2所示，將火箭與井壁間隙(以下簡稱箭井間隙，并用d表示)簡化為環形線匯(文獻[14]中環形線匯模型僅考慮一個圓環線，此處將箭井間隙考慮為圓環面)，發射井出口處流場可用環形線匯流場描述。

圖2 環形線匯流場中某點的速度勢

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

設引射流場強度為Q，如圖2所示，以發射場坪中心為原點，原點與排焰道口連線方向(x軸)為方向角起始位置，場坪上一點到原點的距離r為徑向距離，向徑與初始方向夾角θ為方向角，垂直發射場坪向上為z軸正向建立圓柱坐標系。另設火箭半徑為R1，發射井半徑為R2，取箭井間隙位置半徑為R的圓為參考圓，場坪上空中一點M(r,ψ,z)受到參考圓上一點N(R,θ,0)作用的速度勢[16]為

(1)

其中

在箭井間隙區域對φ積分，可得M受到整個環形線匯作用的速度勢為

(2)

(3)

式中K(ζ)為第一類橢圓積分[17]；ζ為模數。

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

Q=v·π(R22-R12)

(4)

于是

m=v

(5)

將式(5)代入式(3)，可得

(6)

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

設ur、uθ、uz分別代表發射場坪徑向、周向、軸向的速度，根據勢流理論：

(7)

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

(8)

式中E(ζ)為第二類橢圓積分。

同理可得

(9)

(10)

1.4 環形線匯流場理論計算

根據式(8)～式(10)，取v=100(文中單位均取為國際標準單位)，R∈(1, 1.1)，z∈(0.05, 12)，r∈(-12, 12)，采用數值積分法計算得到發射場坪速度幅值的分布云圖及流線圖如圖3所示，橫軸表示r坐標，縱軸為z坐標，黑色曲線為速度幅值等值線，紅色曲線為流線。

取v=100，R∈(1, 1.1)，z=0.1，r∈(-12, 12)，計算徑向速度如圖4所示。取v=100，R∈(1, 1.1)，z∈(0.05, 12)，r=1.05，計算軸向速度如圖5所示。

圖3 流場速度云圖及流線圖

圖4 理論計算徑向速度

圖5 理論計算軸向速度

由圖3可看出，流體以直線r=0為對稱軸向箭井間隙匯聚；在距離(0, 0)較遠時，速度幅值等值線形成以(0, 0)點為中心的圓，在距離(0, 0)較近時，速度幅值等值線分別形成以左右箭井間隙為中心的封閉區域；速度幅值等值線越靠近箭井間隙越密集，顯示出速度在箭井間隙附近迅速增大。

由圖4可見，發射場坪地面附近徑向速度在靠近箭井間隙時急劇增大，這顯示出流體向井的環形出口處匯聚。由圖5可見，在箭井間隙附近流體軸向速度沿z軸負向迅速增大，隨高度增加，速度逐漸趨近于0。

綜上可見，在靠近箭井間隙時，圓環面的積分結果為有限值，規避了奇點的影響。由于采用數值積分方法，在z=0時，積分結果為不定值，與實際情況不符。而實際流場由于粘性的影響，地面附近速度為0，所以計算中，將z=0位置發射場坪地面速度置為0，箭井間隙處速度置為100。

2 FLUENT仿真模型

2.1 仿真模型的建立

近似理論公式的驗證本應以試驗為可靠可信的方法，但進行試驗驗證的成本較高，難度較大，而CFD方法成本較低，實現較為容易，且準確度有一定保證，故在此采取理論和CFD兩種方法進行計算和結果對照分析，說明理論公式和CFD方法的可靠性。本文以FLUENT作為CFD計算分析的工具。為簡化分析計算，采用圓柱軸對稱模型，選取的幾何模型如圖6所示，軸對稱處理后的平面模型如圖7所示。

圖6 流場幾何模型

圖7 流場圓柱軸對稱幾何模型

為估計火箭半徑R1以及箭井間隙d大小對環形線匯模擬發射場坪流場的準確度的影響，R1值分別取為1.0、2.0、3.0，d值分別取為0.1、0.5、0.9，建立相應的幾何模型。為幾何模型劃分四邊形網格，以R1=1.0，d=0.5為例得到的網格模型如圖8所示。

邊界條件和監測線的設置如圖7所示，其中徑向監測線位置為z=0.1，r∈(-12，12)，軸向監測線位置為z∈(0.05，12)，r=1.05。在軟件界面中選擇圓柱軸對稱、穩態模型，κ-ε湍流模式[19]，以常溫常壓全局初始化，入口速度分別選擇-100、-200、-300進行計算，得到監測線上的速度數據。

2.2 網格無關性驗證

文中以R1=1.0，d=0.5時的幾何模型進行網格無關性驗證，得到適用于本項研究的網格密度。分別對網格數為19 641、31 376、42 455、56 204四套網格進行計算，得到徑向監測線上的速度幅值如圖9所示。可見，網格數為19 641的網格模型其監測線上速度幅值的計算結果與其他三套網格模型差距明顯，最大值為88.434 7，而其他三套分別為81.252 7、81.528 6、82.636 6。為保證計算精度，同時有效減少計算量，本文選取網格數為31 376的網格尺寸作為網格劃分的方案。

圖8 網格模型

圖9 不同網格模型計算速度比較

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

為了驗證本文所采用的環形線匯方法和仿真方法的準確性，分別在不同入口速度、不同火箭半徑和不同箭井間隙條件下對兩種方法的計算結果進行對比，實現兩種方法的相互驗證。

3.1 不同入口速度下結果對比

取火箭半徑R1=3，箭井間隙d=0.9，入口速度分別為100、200、300，采用兩種方法計算得到的軸向監測線上的軸向速度如圖10所示，徑向監測線上的徑向速度如圖11所示。從圖10、圖11可看出，理論計算和FLUENT仿真結果吻合很好，同時也可看到軸向速度和徑向速度的最大幅值隨著引射入口速度的增大而增大，而且近似呈線性比例。

(a)v=100 (b)v=200 (c)v=300

(a)v=100 (b)v=200 (c)v=300

3.2 不同火箭半徑下結果對比

取箭井間隙d=0.5，入口速度為100，火箭半徑分別為1、2、3，采用兩種方法計算得到的軸向監測線上的軸向速度如圖12所示，徑向監測線上的徑向速度如圖13所示。可看出，理論計算和FLUENT仿真結果吻合很好，但火箭半徑對軸向速度和徑向速度的最大幅值沒有影響。觀察徑向速度曲線可發現，隨著火箭半徑的增大，箭井間隙靠近火箭一側的速度峰值也就是圖13中徑向速度的正向最大值緩慢增大，這說明引流速度的最大值和流量密切相關。

(a)R1=1 (b)R1=2 (c)R1=3

3.3 不同箭井間隙下結果對比

取火箭半徑R1=1，入口速度為100，箭井間隙分別取0.1、0.5、0.9，采用兩種方法計算得到的軸向監測線上的軸向速度如圖14所示，徑向監測線上的徑向速度如圖15所示。可見，理論計算和FLUENT仿真結果吻合很好，同時可得到速度的最大幅值隨著箭井間隙的增大而逐漸增大。

(a)R1=1 (b)R1=2 (c)R1=3

(a)d=0.1 (b)d=0.5 (c)d=0.9

(a)d=0.1 (b)d=0.5 (c)d=0.9

4 發射場坪流場特征曲線

取v=100，R∈(1, 1.1)，z∈(0.05, 12)，r∈(-12,12)，采用數值積分法計算得到發射場坪速度幅值局部放大云圖如圖16所示，徑向、軸向速度云圖如圖17所示。從圖16可看出，流體從遠場向箭井間隙匯聚，在靠近發射場坪中心位置，流線發生轉彎，從場坪中心折向箭井間隙。設有與z軸平行的直線z∥與某一流線切于A點，則在A點處，徑向速度(如圖16中水平方向速度)方向發生變化，A點的徑向速度為0。在A附近的很小范圍內，由于速度方向變為豎直，速度的水平分量減小，豎直分量相對出現小范圍內的極大值。這樣數目眾多的切點匯聚起來就形成一條連續變化的曲線，如圖17曲線1所示。

圖16 速度幅值局部放大云圖

在圖17(a)中，云圖等值線除過匯聚形成曲線1之外，還在對稱軸上形成了曲線2。在曲線1和曲線2上，徑向速度為0。

在圖17(b)中，設有z軸的平行線z′、z″與軸向速度的等值線相切，則切點處的軸線速度幅值為該切線上最大值。連接這些切點形成曲線1和曲線3。

下面分別求取曲線1、2、3的方程。

曲線2是由流場的軸對稱性質導致，可直接得到其徑向坐標為r=0。

曲線1有兩種求法，一是根據徑向速度為0求得，二是根據軸向速度幅值出現極值來求解。先采取第一種方法。

在曲線1兩邊，徑向速度的方向發生了改變，曲線1上的徑向速度為0。令ur=0，可得

(11)

于是有，被積函數(設為F)在(R1，R2)上正負積分面積相等，或者被積函數為0。實際上，被積函數的取值隨著v、r、z、R變化，積分變量為R，取v=100、r=15、z=15，可得F隨R的變化曲線如圖18所示。從圖18可知，被積函數F在(R1，R2)上正負積分面積相等無法在大范圍內處處成立，所以被積函數為0。

當被積函數F為0時，可得

(12)

上式涉及到橢圓積分，難于求解，可對橢圓積分級數展開取近似。第一類、第二類橢圓積分的級數展開式為[18]

(13)

取二階近似

(14)

將式(14)代入式(12)，可得

r(r2+rR-2R2+z2)=0

(15)

解得r=0(即再次圖17(a)中曲線2)，或者

(16)

式(16)中，當r>0時，R取正值；當r<0時，R取負值。取R為1.05和-1.05，畫出徑向速度云圖特征曲線，經過與云圖比較，發現式(16)的計算結果偏大，在將橢圓積分取三項近似后仍不夠準確。因此，不能將式(16)作為曲線1的近似方程。

下面采取第二種方法求取曲線1。在圖17(b)中，曲線1和曲線3均為軸線方向軸向速度幅值的極值點匯集而成，可一并求出。

從軸向速度方程(10)可得曲線1和3的準確方程為

(17)

同徑向速度云圖特征曲線的分析，可知被積函數為0，即

(18)

將式(14)代入式(18)，可得

(19)

式(19)中第二式顯然不為負值，所以無實數解。由式(19)中第一式取非負實數解得到

當r>0時，R取正值：

(20)

當r<0時，R取負值：

(21)

取R為1.05和-1.05，畫出曲線1和曲線3如圖19所示。設曲線3與直線r=0的交點為S′，提取圖17(a)中S′點的縱坐標值為0.75，提取圖19中S′點縱坐標值為0.742 5，二者相差很小。

圖19 軸向速度云圖的近似特征曲線

綜合以上分析，在圓柱軸對稱坐標系下僅考慮半徑方向正半軸得到發射場坪流場的特征曲線方程：

(22)

5 結論

本文采用改進的環形線匯方法對發射井口引射流場進行研究，推導了速度勢和速度的計算公式，在FLUENT軟件中建立引射流場計算模型并進行了求解，在不同R1值、d值和v值條件下對環形線匯與軟件仿真兩種方法流場計算結果進行了對照驗證，根據環形線匯速度計算公式采用數值積分法求取了速度云圖，提取了流場特征曲線并推導了近似表達式，得到了以下結論：

(1)不同R1值、d值和v值條件下理論計算和軟件仿真速度吻合很好，表明改進的環形線匯理論速度計算公式可較好地描述引射流場流速。

(2)引射流場的流體在離箭井間隙較遠時速度較小，在接近箭井間隙時速度急劇增大；流場徑向流速在箭井間隙范圍以內和以外方向相反，這和實際情況相符。

(3)隨著入口速度和箭井間隙的增大，軸向速度和徑向速度的最大幅值顯著增大；隨著火箭半徑的增大，軸向速度和徑向速度的最大幅值并無明顯變化，但徑向速度的正向峰值緩慢增大。

(4)特征曲線反映了發射場坪流場的軸對稱性質和鏡像對稱性質，其近似求解方法可較準確地描述特征曲線的位置特征，為分析發射場坪流場速度分布提供了有力工具。