發(fā)射井口引射流場理論計算與仿真①

謝 建，權(quán) 輝，謝 政，韓小霞

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

0 引言

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

對于火箭發(fā)射時的流場，以往的研究主要關(guān)注點在火箭發(fā)動機(jī)燃?xì)馍淞髁鲌?。從單純針對火?導(dǎo)彈燃?xì)饬鲌鲎杂缮淞饕?guī)律[4-5]，慢慢擴(kuò)展到關(guān)注某一具體發(fā)射裝置中燃?xì)馍淞鞯牧鲃右?guī)律，例如地下發(fā)射井流場規(guī)律的研究[6]，艦載導(dǎo)彈垂直發(fā)射裝置燃?xì)饬鲌龅难芯縖7-8]，同心筒中燃?xì)饬鲌龅难芯縖9]，雙噴管發(fā)動機(jī)燃?xì)饬鲌龅难芯縖10]等。近年來，又有學(xué)者關(guān)注到燃?xì)馍淞鞯囊恍┘?xì)節(jié)現(xiàn)象，如固體火箭發(fā)動機(jī)燃?xì)馍淞鹘鼒鲂纬珊桶l(fā)展的機(jī)理[11]，燃?xì)饬鲌龅亩稳紵F(xiàn)象以及點火壓力脈沖形成機(jī)理和影響因子[12-13]。這些研究均采用了CFD方法，部分還進(jìn)行了大量實驗驗證，具有較高的準(zhǔn)確度，對于火箭發(fā)射技術(shù)的發(fā)展具有重要指導(dǎo)意義，但當(dāng)前專門針對發(fā)射場坪流場的研究卻鮮見報道。

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

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

1 改進(jìn)的環(huán)形線匯方法

1.1 工程背景介紹

引射流場中流體介質(zhì)主要為場坪上的空氣，溫度為常溫，流動過程中流速低于聲速，相對于高溫高速燃?xì)馍淞鳎勺鳛椴豢蓧嚎s流體進(jìn)行處理。忽略空氣黏性和重力，引射流場可用勢流法[15]進(jìn)行計算。某型發(fā)射井結(jié)構(gòu)圖如圖1所示。

(b) Top view

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

圖2 環(huán)形線匯流場中某點的速度勢

1.2 環(huán)形線匯流場速度勢的推導(dǎo)

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

(1)

其中

在箭井間隙區(qū)域?qū)Ζ辗e分，可得M受到整個環(huán)形線匯作用的速度勢為

(2)

(3)

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

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

Q=v·π(R22-R12)

(4)

于是

m=v

(5)

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

(6)

1.3 環(huán)形線匯流場速度的推導(dǎo)

設(shè)ur、uθ、uz分別代表發(fā)射場坪徑向、周向、軸向的速度，根據(jù)勢流理論：

(7)

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

(8)

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

同理可得

(9)

(10)

1.4 環(huán)形線匯流場理論計算

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

取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)較遠(yuǎn)時，速度幅值等值線形成以(0, 0)點為中心的圓，在距離(0, 0)較近時，速度幅值等值線分別形成以左右箭井間隙為中心的封閉區(qū)域；速度幅值等值線越靠近箭井間隙越密集，顯示出速度在箭井間隙附近迅速增大。

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

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

2 FLUENT仿真模型

2.1 仿真模型的建立

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

圖6 流場幾何模型

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

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

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

2.2 網(wǎng)格無關(guān)性驗證

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

圖8 網(wǎng)格模型

圖9 不同網(wǎng)格模型計算速度比較

3 理論計算結(jié)果與仿真結(jié)果的對比

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

3.1 不同入口速度下結(jié)果對比

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

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

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

3.2 不同火箭半徑下結(jié)果對比

取箭井間隙d=0.5，入口速度為100，火箭半徑分別為1、2、3，采用兩種方法計算得到的軸向監(jiān)測線上的軸向速度如圖12所示，徑向監(jiān)測線上的徑向速度如圖13所示?？煽闯觯碚撚嬎愫虵LUENT仿真結(jié)果吻合很好，但火箭半徑對軸向速度和徑向速度的最大幅值沒有影響。觀察徑向速度曲線可發(fā)現(xiàn)，隨著火箭半徑的增大，箭井間隙靠近火箭一側(cè)的速度峰值也就是圖13中徑向速度的正向最大值緩慢增大，這說明引流速度的最大值和流量密切相關(guān)。

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

3.3 不同箭井間隙下結(jié)果對比

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

(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 發(fā)射場坪流場特征曲線

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

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

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

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

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

曲線2是由流場的軸對稱性質(zhì)導(dǎo)致，可直接得到其徑向坐標(biāo)為r=0。

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

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

(11)

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

當(dāng)被積函數(shù)F為0時，可得

(12)

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

(13)

取二階近似

(14)

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

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

(15)

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

(16)

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

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

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

(17)

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

(18)

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

(19)

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

當(dāng)r>0時，R取正值：

(20)

當(dāng)r<0時，R取負(fù)值：

(21)

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

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

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

(22)

5 結(jié)論

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

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

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

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

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