權輝,謝建,李良,張力
(火箭軍工程大學 導彈工程學院,陜西 西安 710025)
在火箭地下井式發(fā)射中,短時間內由于高溫高速火箭燃氣射流的作用,在發(fā)射井內部和發(fā)射場坪上會產生劇烈變化的流場[1],致使其中設施設備和雜物承受巨大的流場作用力,可能導致設備損壞或者雜物被氣流卷起、撞擊到火箭箭體和其他設備,對發(fā)射安全造成嚴重影響,因此有必要對火箭發(fā)射時井內及場坪的流場展開研究。但是由于武器裝備的保密性和敏感性,當前對于火箭發(fā)射井的研究缺乏廣泛性和針對性,少有公開報道,即便有部分公開報道也語焉不詳、疏于細節(jié),更缺乏針對某一問題的前赴后繼的深入理論研究。正是由于這些原因,致使火箭發(fā)射井問題的理論滯后于實踐,面對技術迭代的需求,在發(fā)射井安全以及優(yōu)化設計相關問題的討論中,存在大量缺乏理論支撐的情況。
在以往的火箭發(fā)射井流場方面的研究中,王飛等[2]、周笑飛[3]、謝政等[4]、謝建等[5]采用計算流體力學(CFD)方法對現有發(fā)射井內流場規(guī)律進行了深入分析,得出火箭發(fā)射時流場的三段變化規(guī)律,并從歷次發(fā)射試驗結果得到了印證,對復燃、注水降噪等問題進行了探討,提出了發(fā)射井改進方案。但是這些有限的研究并未考慮發(fā)射井流場的各個細節(jié),無法系統(tǒng)解決火箭發(fā)射井安全以及優(yōu)化設計問題。從以往的研究結果看,對火箭發(fā)射井內流場分析得較多,對發(fā)射場坪關注較少,且已有研究均采用CFD方法為主、少量試驗數據驗證為輔的方法。當面對發(fā)射井優(yōu)化設計問題時,初始設計階段采用CFD方法會存在同一參數變化時反復建模、消耗大量計算資源的情況,嚴重影響了設計效率。如果建立發(fā)射井內及場坪流場的近似計算公式,則不但可以為CFD結果提供參照,更可為發(fā)射井分析以及優(yōu)化提供簡便可操作性強的方法。
火箭發(fā)射時,由于引射作用,在井口附近形成引流流場,在排焰道口附近形成亞聲速射流流場[3]。對于井口附近引流流場,文獻[6]已經采用勢流法進行了部分研究,而對于排焰道口射流,目前尚未見針對性理論研究。
當前,對于射流問題的研究主要集中在射流穩(wěn)定性和擬序結構[7-9]、沖擊射流強化換熱[10-12]、浮射流[13]等方面,對于彎曲射流,主要集中于實驗和數值方法研究橫向流體作用下射流的流場特性[14-15]。排焰道口射流是處于井口引流流場中的彎曲射流,引流流場對于射流相當于非均勻風場,對于非均勻風場中的彎曲射流尚未有完善理論。
為了準確深刻地把握發(fā)射場坪射流的變化發(fā)展規(guī)律,本文基于前人關于直線射流和彎曲射流的理論研究成果[16-17],推導準確度更高的橫風作用下射流軸線方程和軸線速度計算公式,并與Fluent軟件仿真結果進行比較,驗證方法的準確性。
近似理論公式的驗證本應以試驗為可靠可信的方法,但是進行試驗驗證的成本較高、難度較大,而CFD方法成本較低,實現較容易,且準確度有一定保證,故本文采取理論和CFD兩種方法進行計算和結果對照分析,以表明理論公式和CFD方法的可靠性。
恒定風速風場作用下射流的彎曲變化規(guī)律已有潘衍強[16]、平浚[17]推導得到的近似計算公式描述,但是準確度并不高,且非均勻風場下射流的彎曲變化規(guī)律沒有成熟的理論描述。本文首先證明射流的剩余動量流率沿初始軸向位置(用x軸表示)守恒,然后對射流半寬變化規(guī)律和質量流量變化規(guī)律作出假設,結合向心力和曲率半徑計算公式,求得橫風作用下射流軸線近似方程和軸線速度。
為了方便理論建模和簡化分析,作出以下設定:
1)火箭發(fā)射場坪排焰道口射流為亞聲速矩形彎曲射流且彎曲方向固定,研究表明[18-19]矩形射流在經過充分發(fā)展后截面形狀會演變?yōu)閳A形,因此為了便于研究,將排焰道口射流考慮為圓形二維彎曲射流。
2)不同噴嘴形狀下,射流以不同入射角度射入環(huán)境流場中,產生的作用效果并不相同。但是本文所關注的問題主要在于射流在橫風場中的分布模式,故均假定射流垂直于噴嘴平面射入環(huán)境流場。
3)在射流的研究中,所研究的射流橫截面方向的尺寸相對于射流長度非常小,從而可以將動量、質量等物理量看作是集中于軸線上。另外,為了便于分析,假定射流只有軸向速度、沒有徑向速度。
4)橫風作用下彎曲射流半寬與軸線長度位置的冪函數呈正比。
5)橫截面對稱分布且速度滿足自相似性。
6)在射流主體段,質量流量的變化服從邏輯斯蒂增長率,即增長有上限。
7)根據假設1,射流橫截面為圓形,然而實際上彎曲射流的橫截面為腎臟形[17], 公式中半寬實際上應為當量半寬,相當于把彎曲射流斷面等效為圓形截面時的射流半徑。本文中公式與CFD彎曲射流解的對照中使用的是對稱面半寬,它與當量半寬并不相同,假設對稱面半寬和當量半寬呈正比關系。
有一射流以初始半徑為b0、初始速度為v0射入無限流體區(qū)域中,并在風場作用下向某一個方向發(fā)生彎曲,即設定射流軸線為二維曲線而非三維曲線。以射流噴口中心O為原點,初始速度方向為x軸方向,射流軸線所在平面上與x軸方向相垂直的方向為y軸,建立坐標系Oyx如圖1所示(下文所稱z軸與x軸、y軸呈右手關系)。圖1中,θ為射流軸線上某點處切線與射流入口平面Oyz的夾角,s表示射流中心軸線。沿射流軸線任意位置取dx厚度的射流微元為研究對象,以射流微元中心O′為原點,建立局部圓柱坐標系O′λx如圖1所示,其中圓柱坐標系半徑方向為平行于y軸的射流微元截面半徑方向(用λ表示),圓柱坐標系軸向平行于x軸,以x軸正向為正方向。
圖1 彎曲射流x軸方向微元示意圖Fig.1 Schematic diagram of infinitesimal element of winding jet in x-axis direction
若設射流沿x軸方向的速度為vx,沿λ軸方向的速度為vλ,由湍流射流的微分方程組[20]可得
(1)
式中:ρ為射流密度;τtλ為湍流切向應力。
對(1)式中的第2個公式兩端沿半徑方向積分,得
(2)
式中:ξ為沿λ軸方向的射流擴展半徑。
對(2)式左端第2項積分[20],得
(3)
(1)式中的第1個公式代入(3)式,得
(4)
式中:ux為風場x軸方向的速度。此處近似認為vx在射流邊界處的值取到ux.
假設湍流切應力在射流邊界處為0[20],有
(5)
(1)式、(4)式、(5)式代入(2)式,得
(6)
即
(7)
由于彎曲射流橫截面并非對稱,中心軸線是指橫截面速度最大位置連接而成的曲線,對應于直線射流中心軸線。以射流中心軸線為s軸,沿流動方向為s軸正向,軸線某點法線為r軸,指向射流彎曲方向為正向,建立隨軸線位置變化的圓柱坐標系O″rs,如圖2所示。圖2中,R為射流軸線上的曲率半徑。設射流軸向速度為vs,環(huán)境流體沿s軸流速為us,沿r軸流速為ur,則沿射流軸向的動量流率可以表示為
(8)
式中:b為沿r軸方向的射流半寬。
圖2 彎曲射流s軸方向微元示意圖Fig.2 Schematic diagram of infinitesimal element of winding jet element in s-axis direction
射流沿軸向的動量流率和沿x軸方向的動量流率有如下關系:
(9)
為了簡化分析,將射流質量看作是集中于軸線上,則射流軸線上某點處的質量流量可以表示為
(10)
(9)式、(10)式代入(7)式,結合射流初始位置參數,可得
(11)
式中:ρ0為射流初始密度;ux0為射流初始位置環(huán)境流體沿射流軸向流速。
設風場y軸方向速度為uy,由圖2可知
(12)
若設射流軸線方程為
x=x(y),
(13)
(14)
x′=tanθ,
(15)
式中:x′=dx/dy.
射流軸線上的曲率半徑可以表示為
(16)
式中:x″=dx′/dy.
微元段受到的法向作用力可以表示為
(17)
式中:k1為繞流阻力系數;ρu為環(huán)境流體密度;2πbds為迎風面積。微元段受到的向心力可以表示為
(18)
對確定的微元段,R可以近似看成常量,于是(18)式可以改寫為
(19)
由于射流彎曲運動的向心力由法向作用力提供,由(17)式和(19)式得
(20)
(20)式代入(11)式,得彎曲射流軸線方程為
(21)
在環(huán)境流場參數已知的情況下,要求解(21)式,需要將射流半寬、密度和速度表達為位置坐標y、x的函數,即對變量分布規(guī)律作出假設。
圓形斷面直線射流的半寬沿射流軸向呈線性變化[20-21],而彎曲射流的半寬變化較為復雜。橫向風作用下和斜向風作用下半寬變化規(guī)律差距較大,不同方向斜向風作用下也存在一定差距[17]。按照橫風作用下彎曲射流半寬與軸線長度位置的冪函數呈正比的假設,半寬的變化規(guī)律可以表示為
(22)
式中:C1為常數;n為半寬指數,與速度比等因素相關,用來表征射流截面積增長的快慢,能夠在一定程度上反映射流形狀。
直線亞聲速射流橫截面內的速度分布具有自相似性,而風場中的彎曲射流在內外兩側的半寬并不對稱,在初始階段尚能保持圓形截面,隨著射流發(fā)展,截面會逐漸變?yōu)槟I臟形[17,21],因此彎曲射流并不滿足橫截面速度的自相似性。根據假設7,對稱面半寬和當量半寬呈正比關系,有
bs=C2b,
(23)
式中:bs為對稱面半寬;C2為常數,取值范圍為0~1,C2反映了彎曲射流截面扁窄的程度。
考慮到橫截面形狀變化的連續(xù)性,隨著軸線長度的增加,不同位置處橫截面歸一化速度分布函數必定形成同一函數空間內連續(xù)變化的一族函數。此時速度分布可以記為
(24)
式中:vm為軸線速度;fs表示函數空間內隨著軸線位置變化的速度分布函數;η為歸一化半寬。
為了便于計算,一般仍假設橫截面對稱分布且速度滿足自相似性,從而對射流進行預估。對于復合射流,在主體段射流時均速度服從如下分布[17]:
(25)
(25)式在直線射流時相當有效,但對于彎曲射流卻無法有效模擬。為此本文提出如下多項式分布律:
(26)
式中:n′為指數,取值為正數,本文為了簡化計算,取固定值2. 根據(26)式,可在已知軸線位置和軸線速度的情況下對射流半寬內的速度分布進行預估。
軸線速度是排焰道口射流問題應重點關注的參數,直線射流的軸線速度與軸線位置呈反比。彎曲射流的軸線速度變化比較復雜,一般假設軸線速度與坐標呈多項式關系[21]。為了保持公式推導的連貫性,在已經對半寬和橫截面內速度分布作出假設的情況下,本文嘗試推導軸線速度的變化規(guī)律。
直線射流時,根據半寬隨著軸線長度線性變化以及軸線速度與軸線長度呈反比[20],可得質量流量隨著軸線長度線性變化。然而彎曲射流的軸線延伸長度遠小于直線射流,不能采用線性假設。
由于橫風的作用,射流流體不斷被吹走,在初始段變化尤其劇烈。因此,在計算主體段射流質量流量時必須考慮初始段的質量流量損失。
為了描述初始段的流量損失,引入初始段流量損失系數k2∈(0,1),用k2m0表示經過初始段損失后的質量流量(其中m0表示初始質量流量),即主體段質量流量的初始值。
隨著軸線長度的增加,彎曲射流的卷吸能力先增加后減小,最后射流消弭于環(huán)境流場。這種變化規(guī)律正好與邏輯斯蒂曲線相符。設定瞬時增長率為C3,最大流量為主體段初始流量的k3倍,則質量流量的變化規(guī)律可以設為
(27)
對(27)式積分,可得
(28)
(24)式、(28)式代入(10)式,可得
(29)
式中:
(30)
將半寬的表達式(22)式代入(29)式,就可以根據軸線位置求解軸線速度。
排焰道口射流為亞聲速射流,為了簡化分析,將其視為不可壓流體。火箭發(fā)射時,由于發(fā)射井的引射作用,在井口附近形成強烈的引射流場,發(fā)射場坪上的雜物在引射作用下有可能被吸入發(fā)射井。排焰道口射流也受到強烈的引射作用,雖然設計中設置了折流角,但仍然會向井口一側發(fā)生彎曲。發(fā)射安全問題重點關注地面附近引射流場的作用,在這種情況下將引射產生的非均勻風場在排焰道口近似為橫向風,既能在一定程度上說明射流彎曲情況也能簡化計算。
設流體不可壓,則射流密度可作為常量。另設環(huán)境流場為
(31)
式中:u為橫風速度。則
(32)
此時軸線方程(21)式可以化為
(33)
式中:
(34)
在求得軸線方程的情況下,可以根據(22)式、(25)式和(29)式求得半寬、軸線速度和速度分布。匯集以上方程,形成如下求解恒定風速橫風作用下不可壓縮射流速度分布的方程組:
(35)
以上方程組反映了射流主體段的速度分布規(guī)律。然而,射流初始段和主體段的流動規(guī)律并不相同,必須對初始段結果進行修正。
初始段高度是求解初始段射流流場的重要參數,一般由經驗公式獲得[20]。本文為了簡化計算和保證計算結果的連貫性,令vm=v0,得到射流核心區(qū)高度x=lc.
根據前述假設,射流速度在初始段混合層的分布也表現為多項式函數,具體表達式可寫為
(36)
式中:bk為射流核心區(qū)半寬,其定義為垂直于射流軸線的核心區(qū)半徑,如圖3所示。圖3中,Ql為對稱平面上射流初始位置左端邊界點,Qr為對稱平面上射流初始位置右端邊界點,Qlc為對稱平面上r軸方向射流左邊界點,Qrc為對稱平面上r軸方向射流右邊界點。
圖3 射流核心區(qū)半寬Fig.3 Half-width of core area of winding jet
假設初始段混合層內外邊界均為直線[20],再假設射流核心區(qū)高度在環(huán)境流場作用下相對于沒有環(huán)境流場時保持不變,則bk可由(37)式計算:
(37)
當x=lc時,可以根據方程組(35)式結果求得射流半寬位置Qlc、Qrc點坐標,進而求得直線QlcQl和直線QrcQr的方程,結合軸線方程可以求得射流初始段半寬。
在求得初始段半寬之后,就可以根據(36)式、(37)式求出初始段速度分布。
取b0=0.01 m,ρ0=ρu,v0=100 m/s,u=20 m/s,θ0=0.5π rad,k1=0.77,C1=0.058,k2=0.64,C2=0.55,k3=4.8,C3=6.1,I1=0.36,計算得到射流速度幅值云圖如圖4所示。為了方便計算,未考慮射流背風一側的渦結構,射流外流場區(qū)域與實際情況并不相符,但是本文關心的主要是射流內部情況,因此圖4中未展示環(huán)境流場。從射流云圖的走向、形狀以及速度變化趨勢來看,近似公式能夠預測射流的大致規(guī)律。
圖4 理論計算射流云圖Fig.4 Velocity contour of calculated winding jet
本文以Fluent軟件仿真結果對比分析彎曲射流方程組計算結果的偏差。幾何模型如圖5所示,射流為圓形入口(圖5中射流入口),流場空間取為棱長為10倍入口直徑的正方體區(qū)域,射流入口位于流場區(qū)域底面中心。為流場區(qū)域劃分網格如圖6所示。射流入口和左側面均設為速度入口,流場區(qū)域底部設為不可滑移壁面,其余面均設為有回流的速度出口,溫度設為300 K. 考慮到所研究流場區(qū)域為亞聲速燃氣湍流流場且已簡化為定常流動,選擇穩(wěn)態(tài)模型、k-ε(k表示湍流脈動動能,ε表示湍流耗散率)湍流模型[14,22-23],以2階迎風格式[24]進行離散,計算得到彎曲射流對稱面上的速度云圖如圖6所示。
圖5 幾何模型Fig.5 Geometric model
圖6 網格模型Fig.6 Mesh model
為了驗證網格無關性,選取射流入口半徑(即射流初始半寬)0.01 m,入口速度100 m/s,流體區(qū)域棱長0.2 m,左側速度入口20 m/s. 分別取網格數為28萬、56萬、86萬和115萬共4組網格進行計算,比較彎曲射流對稱面上射流軸線的位置如圖7所示。從圖7中可以看出,4組網格計算得到的射流軸線位置差異很小,最大縱坐標差距小于3.1%,其中28萬網格的計算結果與其他3組計算結果差距稍大。本文選取56萬的網格密度作為基本依據,根據參數變化建立不同的幾何模型劃分網格并求解。
圖7 網格無關性驗證Fig.7 Comparison of jet axes from different mesh models
以b0=0.01 m、v0=100 m/s、u=20 m/s、θ0=0.5π rad建模仿真,得到對稱平面Oyx上的速度幅值云圖如圖8所示。Oyz平面切片(x-Plane)速度幅值云圖如圖9所示,速度幅值等值面圖如圖10所示。
圖8 Fluent軟件仿真速度云圖Fig.8 Velocity contour from Fluent
圖9 彎曲射流x軸方向斷面形狀變化示意圖Fig.9 Schematic diagram of shape change of winding jet in x-axis direction
比較圖4和圖8可知,圖8完整展現了射流迎風面均勻橫風場和背風面的低速渦結構。在射流初始段橫風將一部分射流流體吹走脫離開射流主流,使主流流量減少。相比圖4,圖8更能準確反映實際情況,但是二者在射流的走向和速度變化趨勢上基本一致。
從圖9中可以看出,隨著沿射流軸向位置的變化,射流截面形狀逐漸由圓形過渡為月牙形,對稱面上的半寬先減小、后增大,而Oxz平面上的半寬逐漸變大,射流中心速度(參照圖8速度標尺)不斷減小。
圖10反映了射流在橫風作用下形狀以及初始段質量流量的損失現象。
圖10 彎曲射流速度幅值等值面圖Fig.10 Iso-surface of winding jet velocity
在本文研究的具體情形下,可供改變的參數為射流初始半寬b0、速度比v0/u. 下面在不同初始半寬和速度比條件下,用Fluent軟件和近似公式分別計算橫風作用下射流速度,調節(jié)近似公式參數使其結果和Fluent軟件仿真計算結果相符合,通過比較參數的變化來分析初始半寬和速度比對射流形狀的影響。
速度比v0/u取為5,初始半寬b0分別取為0.01 m、0.02 m、0.03 m,分別用Fluent軟件仿真和近似公式計算得到橫風作用下射流流場結果。在對稱平面上取得射流軸線、邊界位置和軸線速度幅值,調節(jié)近似公式的參數,使近似公式計算結果與Fluent軟件計算結果相符合。在此過程中得到的軸線位置比較如圖11所示,射流邊界比較如圖12所示,軸線速度比較如圖13所示,近似公式參數如表1所示(為了簡化計算,I1取為定值)。從圖11~圖13中可以看出,近似公式能準確模擬射流軸線位置、邊界和軸線速度,表明本文所提出的近似公式在不同射流初始半寬條件下具有很高的準確度。在表1中,除C1和C3外其余參數變化不大。從C1的變化可以看出,在彎曲射流速度比不變時,隨著初始半寬的增大,C1逐漸增大,表明射流邊界范圍隨著初始半寬的增大而增大,符合實際情況。從C3的變化可以看出,隨著初始半寬的增大,射流質量流量增速隨軸線位置的增大而減小,表明彎曲射流卷吸能力隨著初始半寬的增大而減小。
圖11 不同初始半寬下彎曲射流軸線位置Fig.11 Jet axes under different initial half-widths
圖12 不同初始半寬下彎曲射流邊界Fig.12 Jet boundaries under different initial half-widths
圖13 不同初始半寬下彎曲射流軸線速度Fig.13 Velocities of jet axes under different initial half-widths
b0取為0.01,速度比v0/u分別取為5.0(100/20)、7.5(75/10)、10.0(100/10),分別用Fluent軟件仿真和近似公式計算得到橫風作用下射流流場結果。在對稱平面上取得射流軸線、邊界位置和軸線速度幅值,調節(jié)近似公式的參數,使近似公式計算結果與Fluent軟件計算結果相符合。在此過程中得到的軸線位置比較如圖14所示,射流邊界比較如圖15所示,軸線速度比較如圖16所示,近似公式參數如表2所示。從圖14~圖16中可以看出,近似公式能準確模擬射流軸線位置、邊界和軸線速度,表明本文所提出的近似公式在不同速度比條件下具有很高的準確度。在表2中:除了C1參數不變,其余參數均有變化;n增大表明隨著速度比的增大,射流整體形狀變胖;k1增大表明隨著射流速度比的增大,軸線上法向阻力系數增大;k2增大表明隨著速度比的增大,射流在初始段質量流量損失變小,即被風吹走的流體減少;C2有微小增加,表明隨著速度比增大,射流截面扁平程度減小,與實際情況相符,但是也表明在本文所取的速度比條件下,射流截面變圓程度有限;k3微小變大,表明隨著速度比增大,橫向風對射流最終流量的影響變小,射流最終流量增大,趨近于直線射流;C3減小表明隨著速度比的增大射流卷吸能力逐漸變小;I1減小表明隨著速度比的增加,射流橫截面內速度分布逐漸趨于均勻。
表1 不同初始半寬條件下近似公式可調參數的取值Tab.1 Values of adjustable parameters under different initial half-widths
圖14 不同速度比下彎曲射流軸線位置Fig.14 Jet axes under different velocity ratios
圖15 不同速度比下彎曲射流邊界Fig.15 Jet boundaries under different velocity ratios
圖16 不同速度比下彎曲射流軸線速度Fig.16 Velocities of jet axes under different velocity ratios
表2 不同速度比條件下近似公式可調參數的取值Tab.2 Values of adjustable parameters under different velocity ratios
通過以上兩組比較可以看出:C1取值只與射流初始半寬有關,而與速度比無關;n、k1、k2、k3的變化只與速度比有關,而與射流初始半寬無關;C2受到速度比的影響但影響較??;C3、I1受到初始射流半寬和速度比的雙重影響,其中速度比的影響更為明顯。綜合以上分析,射流初始半寬的變化對射流形狀影響較小,速度比的變化對射流形狀影響較為顯著。
本文推導證明了射流剩余動量流率沿射流初始軸向位置守恒,結合風場中射流向心力計算公式,推導了非均勻風場中彎曲射流的軸線方程。提出了半寬和質量流量的近似公式,推導了均勻橫風作用下不可壓縮射流的速度方程組。用Fluent軟件對均勻橫風作用下的射流進行了模擬。調節(jié)參數使近似公式計算結果與Fluent軟件計算結果相符合,根據參數變化分析了初始半寬和速度比對射流的影響。得到以下主要結論:
1)本文所提近似公式能夠較為準確地模擬恒定風速橫風場中不可壓縮射流軸線位置、對稱面射流邊界、軸線速度以及橫截面速度。
2)在速度比一定的情況下,隨著射流初始半寬的增大,C1逐漸增大,表明射流邊界范圍逐漸增大,C3逐漸減小,表明彎曲射流卷吸能力減小。
3)當射流初始半寬一定的情況下,隨著速度比的增大,射流形狀逐漸變胖,阻力系數增大,初始段質量損失減小,射流最終流量變大,卷吸能力減小,橫截面內速度峰峰比減小。
4)C1取值只與射流初始半寬有關,而與速度比無關;n、k1、k2、k3的變化只與速度比有關,而與射流初始半寬無關;C2受到速度比的影響但影響較?。籆3、I1受到初始射流半寬和速度比的雙重影響,其中速度比的影響更為明顯。射流初始半寬的變化對射流形狀影響較小,速度比的變化對射流形狀影響較為顯著。
5)發(fā)射場坪排焰道口射流為非均勻風場中亞聲速彎曲射流。但是在地面附近引射流場風速變化不大,可以近似為均勻橫風,而在計算精度要求不高時,排焰道口射流也可以近似為不可壓縮射流,由此可見本文提出的恒定風速橫風場中不可壓縮射流的近似公式在分析排焰道口射流時具有可行性。