基于動網格技術的固體燃料沖壓發動機燃面瞬態退移速率研究①

魏 韜，武曉松

(南京理工大學機械工程學院航空宇航系，南京 210094)

基于動網格技術的固體燃料沖壓發動機燃面瞬態退移速率研究①

魏 韜，武曉松

(南京理工大學機械工程學院航空宇航系，南京 210094)

為了研究固體燃料沖壓發動機(SFRJ)燃面退移速率在工作過程中的變化特性，基于發動機工作特點及動網格技術，考慮到燃燒流動及燃料表面的對流、輻射換熱與燃料熱解退移等過程耦合的影響，建立了SFRJ燃面瞬態退移速率預示方法，并對某帶補燃室、以聚乙烯(PE)為燃料的試驗發動機的燃燒室-噴管統一內流場進行數值計算，得到在移動邊界條件下的瞬態流場分布，并分析了內彈道參數云圖及其隨時間的變化規律。結果表明，燃燒主要發生在當量比函數φ在-2～2之間的區域;隨著發動機工作，燃速逐漸降低，且再附點向下游移動，燃料通道出口處流速和溫度有降低趨勢;此外，在小型發動機工作初期，燃料通道尾部出現類似固體火箭發動機的侵蝕燃燒現象。研究表明，該方法能成功求解發動機復雜的非定常工作過程，較好揭示燃面退移過程。所得結論對發動機設計和試驗具有一定指導意義。

固體燃料沖壓發動機;瞬態退移速率;動網格;數值仿真

0 引言

在固體燃料沖壓發動機工作過程中，燃面瞬態退移速率隨時間和軸向位置變化，形成的燃氣內流場幾何邊界十分復雜，流動區域不斷變化，且受結構和工作條件影響，難以用統一公式描述。用試驗方法獲得燃面瞬態退移速率所需系統復雜，成本昂貴，而傳統數值計算方法采用準定常假設，顯然不能準確描述發動機復雜非定常流場的細節。目前，針對SFRJ燃面瞬態退移速率應用動網格技術和瞬態方法，較詳細地反映流場變化規律的研究，國內外很少有報道。

SFRJ燃面退移速率的傳統分析方法［1-4］是基于邊界層假設的傳熱理論，來確定向燃面的換熱量及燃速。研究者通常采用準定常假設［3-6］，忽略熱輻射［1，3-5］，甚 至 假 設 燃 面 溫 度 和 有 效 汽 化 熱 為 常數［2-4］，這樣的簡化不能考慮到各種復雜物理化學過程的相互作用。因此，針對SFRJ在運動邊界條件下的非定常流場，建立更為精確的燃面瞬態退移速率預示方法具有重要意義和必要性。

針對以上問題，基于動網格技術，將燃燒流動及固體燃料表面的對流、輻射換熱與燃料熱解退移等過程耦合計算。采用低Re數k-ε湍流模型，并考慮加質對對流換熱的影響和溫度對燃料有效汽化熱的影響，建立SFRJ燃面瞬態退移速率的預示方法，并對代夫特科技大學試驗發動機［3］的燃燒室、噴管流場進行一體化計算。研究移動邊界條件下的瞬態流場，并分析內彈道參數云圖、燃面瞬態退移速率、再附點、燃料通道出口參數及其隨時間的變化規律。研究結果可為SFRJ的工程設計和試驗提供參考。

1 物理數學模型

1.1 物理模型

圖1為試驗發動機結構簡圖［3］。燃料藥柱通道直徑為45 mm，長度為300 mm，入口直徑為15 mm，突擴臺階高度為15 mm，補燃室長度為100 mm，喉部直徑為20 mm。固體燃料為聚乙烯(PE)。

圖1 試驗固體燃料沖壓發動機簡圖Fig.1 Configuration of experimental SFRJ

1.2 基本假設

固體燃料沖壓發動機工作過程是一個極其復雜的能量轉化過程。其中，燃料熱解退移、燃燒、湍流流動和傳熱等均是復雜的物理化學過程，且相互之間高度耦合。目前，對其進行嚴格的數學描述尚存在諸多困難。因此，本文假設:(1)發動機為二維軸對稱模型;(2)聚乙烯熱解產物只含C2H4，燃氣為純氣相;(3)化學反應為C2H4-空氣兩步反應模型;(4)不考慮壁面燒蝕與傳熱。

1.3 控制方程

1.3.1 流動控制方程

對于具有運動邊界的有限體積控制體，采用任意拉格朗日-歐拉(Arbitrary Lagrange-Euler-ALE)有限體積法描述的積分形式守恒型控制方程為

式中ρ為流體密度;φ為通用變量為燃氣速度為燃面網格運動速度;Γφ為廣義擴散系數;Sφ為廣義源項。對流項和擴散項采用二階迎風格式離散。

1.3.2 氣-固相交界面控制方程

固體燃料和燃氣的耦合，通過交界面質量和能量守恒方程實現，即

式中為燃料熱解產物的質量流率;ρs為固體燃料密度，取926 kg/m3;˙rb為當地燃面退移速率;和分別為對流和輻射換熱項=(Tw-T0)為燃料內部溫升吸熱率，其中cs為固體燃料比熱容，取2 142 J/(kg·K)，Tw為當地燃面溫度，T0為燃料初溫，取300 K;=為燃料熱解、氣化潛熱，其中hv為燃料有效汽化熱。

1.4 數學子模型

除以上控制方程外，還須附加數學子模型來描述化學反應、湍流、燃料熱解、對流和輻射換熱以及燃料有效汽化熱。

1.4.1 氣相化學反應模型

固體燃料沖壓發動機再附點下游的燃料表面上形成湍流邊界層，在邊界層內形成湍流擴散火焰，發動機中的燃燒主要為擴散燃燒，反應速率受氣體擴散過程控制。因此，燃燒模擬采用渦耗散模型(EDM)。用簡化的兩步總包反應模擬氣相化學反應，即

1.4.2 湍流模型

由于燃面存在不均勻的質量和能量輸運，在靠近燃面的區域，湍流脈動動能強烈衰減，而耗散率達到最大值［7］。為使計算能從高Re數區域一直進行到燃面(該處Ret=0)，本文使用考慮了近壁面效應的低Re數k-ε湍流模型(low-Re k-εmodel)［8-9］封閉方程組。湍流輸運引入的源項為

其中，ηt=cμ∣fμ∣，f1=1.0，f2=1-0.22 ×exp［-(Ret/6)2］，fμ=1-exp(-0.011 5y+)，模型常數cμ=0.09，c1=1.35，c2=1.8，σk=1.0，σε=1.3。

1.4.3 固體燃料熱解模型

根據Hadar Ian［10］等人的研究，碳氫型固體燃料的退移速率與燃面溫度的關系符合Arrhenius公式，即

式中A、Ea和R分別為指前因子、活化能和通用氣體常數。

模型常數通過熱分解實驗［11］獲得，對于PE，A=8.25 ×105mm/s，Ea=133 539.35 J/mol。將式(6)代入式(3)，采用牛頓迭代法解得Tw，進而得到當地燃速。

1.4.4 對流換熱模型

燃氣向固體燃料對流換熱熱流密度為

式中T∞為邊界層內火焰面溫度;h為對流換熱系數。

對圓管內壁有質量(即燃料熱解氣體)加入的湍流流動，h的計算式為［12］

式中cp、ρ和Pr分別為燃面附近燃氣定壓比熱容、密度和普朗特數;uin為入口流速;Red為按uin計算的通道雷諾數為平均燃速。

1.4.5 輻射換熱模型

燃氣向燃料輻射換熱凈熱流密度為［13］

式中σ為Stefan-Boltzmann常數;ε＇w=(εw+1)/2為燃料的有效發射率，將燃料壁面考慮為漫射的灰體表面［14］，并考慮到炭黑附著，發射率取εw=0.8;εg和αg分別為燃氣的發射率和吸收率;Tc為燃氣靜溫。

熱力計算表明，火焰區主要氣體為CO2、H2O、N2和O2，其中N2和O2為非極性對稱型雙原子氣體，對熱射線的發射和吸收能力微弱，可認為是透明體，故不考慮其對輻射換熱的貢獻，而CO2和H2O等多原子氣體是主要輻射源。因此，燃氣發射率εg和吸收率αg分別為

查文獻［13］中圖 7-1到 7-3并計算得，εg=0.024 85，αg=0.175 55。

1.4.6 燃料有效汽化熱

燃料的有效汽化熱與溫度有關，并取決于傳熱機理，準確的hv對于精確預測燃速非常關鍵［3］。根據文獻［15］，hv可表達為

式中cp為熱解產物的定壓比熱容;hmel為解鏈反應熱;hvap為熱解產物的汽化熱;hpyr為產生熱解產物所需的反應熱。

對于 PE，hmel=225 kJ/kg，hvap=485 kJ/kg，hpyr=3 335 kJ/kg。

Rihani-Doraiswamy(RD)基團貢獻法在求解SFRJ燃燒室環境中化合物的定壓比熱容時有更高的精度［6］，它將cp表示為

式中ni為第i種基團的數量;ai、bi、ci和di分別表示不同基團對化合物定壓比熱容的貢獻。

應用RD法所得C2H4定壓比熱容的基團貢獻值見表1。C2H4的定壓比熱容為

表1 Rihani-Doraiswamy基團貢獻法C2H4定壓比熱容的基團貢獻值Table 1 Heat capacity of C2H4organic compounds from group contribution by RD method

1.5 初始計算網格

利用Gambit生成計算網格，局部網格見圖2。由于燃面節點退移速率不同，變形后燃面形狀復雜。為了動網格的實現，在燃面附近采用加密的非結構網格，以適應不規則的氣流通道。為提高計算精度和速度，對計算域進行分區，在遠離燃面的區域，采用稍稀疏的四邊形結構網格，并對喉部和近壁區進行加密。初始網格數為48 746。

圖2 局部計算網格(燃燒室頭部)Fig.2 Partial computation mesh

1.6 邊界及初始條件

初始條件采用初場賦值方法:首先在定網格條件下，利用DEFINE_PROFILE宏定義燃速，計算定常流場。待求解收斂后，采用動網格技術模擬燃面退移的非定常流場。

邊界條件:(1)燃面邊界。包括加質和移動邊界條件。燃面節點退移速率由式(6)給出。(2)固體壁面邊界。絕熱、無滑移條件。(3)軸對稱邊界。沿軸線法向速度為零，各物理量在軸線上的法向梯度為零。(4)空氣入口邊界。采用質量流率入口邊界條件=0.15 kg/s，總溫T0=600 K。(5)噴管出口邊界。噴管出口為超音速流，流動參數由上游參數二階外推插值得到。

1.7 動網格技術

為了動態模擬SFRJ工作過程中燃面瞬態退移過程，需使用動網格技術，根據計算得到的當地燃面退移速率實時更新網格。

1.7.1 動網格更新算法

動態網格更新算法通常有動態層鋪法、彈性光順法及局部網格重構法。針對燃面變形復雜和運動幅度較小的特點，本文聯合使用彈性光順法和局部網格重構法。

(1)彈性光順法。該網格更新算法假定相鄰兩節點之間的邊為理想彈簧。變形前的初始位置為平衡狀態，當節點運動時，與其相連的各邊將產生正比于邊長變化的力。將作用于節點的力寫為胡克定律的形式，即

平衡狀態下，作用于節點的合力為零，可得迭代方程為

當得到節點位移后，通過對內部節點的Jacobi迭代來求解方程。當迭代收斂后，則按式(17)更新節點位置

式中n+1和n分別表示下一時刻和當前時刻的節點位置。

(2)局部網格重構法。當邊界位移較小時，僅靠網格變形就能適應邊界運動，則采用彈性光順法;當邊界位移相對當地網格尺寸較大時，網格變形會形成嚴重扭曲的單元，或發生網格退化，下一時刻的求解可能會產生收斂性問題。局部網格重構法將不符合畸變率或尺寸標準的網格聚集起來，通過插值運算，在局部重新生成網格，以保證計算精度。若新網格滿足畸變率和尺寸標準，就對網格進行局部更新，否則放棄新網格。

1.7.2 計算流程

圖3為流場計算程序流程圖。主程序主要有更新流場參數、流場計算和動網格技術更新網格等模塊。其中，更新流場參數模塊在每個時間步開始時，利用上一時刻流場參數的計算結果，對當前流場進行更新;流場計算模塊為Fluent提供的流場計算程序;動網格模塊的作用是在當前時間步下流場計算收斂后，計算并更新燃面在下一時刻的位置，同時添加源項。其中，燃面的運動利用DEFINE_GRID_MOTION宏定義。

圖3 流場計算程序流程圖Fig.3 Flow chart of CFD program

1.8 模型校驗

為驗證所建立計算方法的準確性，對文獻［3］的試驗發動機進行計算，燃面退移速率計算值與試驗結果的對比如圖4所示。可見，兩者基本吻合，表明所建立的預示方法適用于SFRJ燃面退移速率的計算。

圖4 燃面退移速率計算值與試驗數據對比Fig.4 Comparison between experimental data and computational results

2 計算結果與分析

選用Fluent提供的隱式、非定常耦合求解器，采用UDF編程方式處理燃面退移。選取時間步長為0.05 s，迭代步數為600，內迭代最大步數為100。由于發動機長徑比較大，為有效描述內流場特性，將發動機的X軸(軸向)與Y軸(徑向)坐標之比設為1/3。

圖5為14.0 s時刻發動機中速度及流線分布圖。可看出，突擴臺階后及節流板后形成兩個回流區，它們將在一定程度上強化空氣和燃料的摻混，并延長氣體在燃燒室的停留時間，有利于擴散燃燒的充分進行;此外，突擴臺階后的回流區還會加強燃氣對燃料的對流換熱，這有利于提高局部燃面退移速率。發動機軸線附近的流速在突擴臺階和節流板后，受回流渦旋的影響有一個減小過程;流速沿軸向有較大變化，這是由于燃料熱解加質、燃燒及通道結構引起的;同時，燃面附近流速很小，且存在回流區，這正是SFRJ能維持穩定燃燒的原因。

圖5 發動機中速度和流線圖Fig.5 Velocity and path-line in SFRJ

圖6為發動機工作過程中速度沿軸線的變化曲線(噴管加速過程未完全畫出)。可見，軸線上的速度變化趨勢相同;燃氣進入補燃室之前，流速有先減小后增大的過程，減速是由于節流板的阻滯，使流動以較大的徑向速度向軸線匯集，加速是由于節流板通道直徑較小，流動連續性要求使然;此外，在突擴臺階和節流板后流速有不同程度的降低，這也印證了對圖5的速度分析。應當指出，隨著發動機工作，燃料通道截面逐漸增大，通道出口處(x=0.29 m)速度逐漸降低，在6 s時刻，速度為130 m/s，到26 s時刻，速度僅為102 m/s;燃燒室入口處的速度基本維持在180 m/s，這是由于空氣質量流率保持不變。

為直觀地表述燃燒室中燃料和氧化劑的分布狀況，引入組分當量比函數φ=lg(Yoxid/Yfuel/φ0)。其中，Yoxid和Yfuel分別為當地氧化劑和燃料的摩爾分數;φ0為氧化劑和燃料的摩爾恰當比。在φ=0的區域，氧化劑和燃料完全反應;在φ＜0的區域，燃料剩余;在φ＞0的區域，氧化劑剩余。

圖6 不同時刻沿發動機軸線的速度曲線Fig.6 Curves of velocity along axis at different time

圖7為14.0 s時刻發動機中組分當量比函數φ和溫度等值線圖。可見，燃燒主要集中在燃料表面和軸線之間的一定區域內(-2＜φ＜2)，在此區域外，氧燃比偏離反應恰當比φ0較大，幾乎無燃燒發生;軸線前半部及燃面附近氧氣與燃料混合較差，故溫度較低;火焰面最高溫度達2 400 K，沿軸線向后，由于氧氣的消耗，燃燒區逐漸擴展至軸線上;此外，補燃室溫度達2 200 K以上，說明補燃室可提高燃燒效率。

圖7 發動機中當量比函數和溫度等值線圖Fig.7 Contours of species equivalent ratio function and temperature in SFRJ

圖8為不同時刻發動機軸線上溫度變化曲線。可看出，溫度變化規律相同;結合圖7可知，在燃燒室前部，由于氧氣充足，燃料未擴散至軸線處即被消耗，因此溫度基本等于來流溫度;在x=0.2 m處，溫度開始迅速上升，這是由于沿軸線向后，隨著氧氣消耗，燃料與氧氣在越來越靠近軸線的區域混合充分，并發生擴散燃燒;此外，補燃室軸線溫度相對穩定，氣體在噴管中膨脹加速，因而溫度迅速下降。應當指出，隨著發動機工作，燃料通道出口處溫度有降低趨勢，在6 s時刻溫度為2 435 K，到30 s時刻溫度為2 285 K。這是因為燃料通道截面逐漸增大，空氣相對集中在軸線附近，燃料和氧氣混合更困難，使火焰面遠離軸線。

圖8 不同時刻沿發動機軸線的溫度曲線Fig.8 Curves of temperature along axis at different time

圖9為發動機工作過程中燃面瞬態退移速率變化曲線。可發現，燃面退移速率沿軸向有較大變化，且不同時刻燃面退移速率的變化趨勢基本相同，都是先增大后減小，這是由于燃面退移速率與熱解表面溫度之間符合Arrhenius熱解速率公式。這造成燃料通道呈現出兩頭細中間粗的型面，如圖5和圖7所示;緊靠突擴臺階后面存在一個流動滯止點，燃料難以與氧氣充分反應，使局部燃料絕對富余(當量比函數φ為-7左右)，因此燃速很低，約為0.05 mm/s;在此滯止點后，回流區的產生強化了燃料和氧氣的摻混，以及燃氣向燃料表面的對流換熱，使得燃速迅速上升，并在再附點處達到峰值，在6 s時刻再附點(x=104 mm)燃速為0.399 mm/s;峰值后邊界層充分發展，隨著氧氣的消耗，火焰面向軸線靠近，對固體燃料熱解表面的傳熱減弱，因而燃面退移速率逐漸降低。

隨著發動機工作，燃面退移速率逐漸降低，這是由于燃料通道截面逐漸增大，使空氣相對集中在軸線附近，燃料和氧氣在遠離燃面的區域發生擴散燃燒，從而傳向燃面的熱量減少。還可發現，隨著發動機工作，再附點逐漸向下游移動。這是由于燃料通道直徑增大，使得臺階高度增加引起的。此外，在發動機工作初期(t＜10 s)，燃料通道出現類似于固體火箭發動機中的侵蝕燃燒現象，即燃料通道尾部的局部退移速率呈緩慢增大趨勢。這是因為燃燒初期小型發動機的燃氣通道較小，平行于燃料表面的流速大于出現侵蝕效應的臨界速度。

圖9 不同時刻燃面瞬態退移速率曲線Fig.9 Curves of instantaneous regression rate of the solid fuel grain at different times

3 結論

(1)通過對比計算結果與試驗數據，驗證了所建立的固體燃料沖壓發動機燃面瞬態退移速率預示方法的準確性。該方法成功求解了發動機復雜的非定常工作過程，較好地揭示了燃面不規則退移過程。

(2)固體燃料沖壓發動機的流場具有明顯不均勻性，燃燒主要發生在燃料表面和軸線之間的一定區域，該區域的當量比函數φ在-2～2之間。

(3)在燃燒室入口處，燃速很低，之后形成的回流區強化了燃料和氧氣的摻混及對流換熱，使局部燃面退移速率迅速升高，并在再附點處達到峰值，峰值后平緩降低。

(4)隨著發動機工作，燃面退移速率逐漸降低，且再附點逐漸向下游移動;燃料通道出口處的流速和溫度有降低趨勢。

(5)在小型發動機工作初期，燃料通道尾部出現類似固體火箭發動機的侵蝕燃燒現象。

［1］Krishnan S，Philmon G.Solid fuel ramjet combustor design［J］.Progress in Aerospace Sciences，1998，34:219-256.

［2］Korting P A O G，Van der Geld C W M，Wijchers T，et al.Combustion of polymethylmethacrylate in a solid fuel ramjet［J］.Journal of Propulsion and Power，1990，6(3):263-270.

［3］Elands P J M，Korting P A O G，Wijchers T，et al.Comparison of combustion experiments and theory in polyethylene solid fuel ramjets［J］.Journal of Propulsion and Power，1990，6(6):732-739.

［4］向敏.固體燃料沖壓增程炮彈工作過程仿真及性能分析研究［D］.長沙:國防科技大學，2006.

［5］劉巍.固體燃料沖壓發動機燃燒組織技術研究［D］.長沙:國防科技大學，2010.

［6］夏強.固體燃料沖壓發動機工作過程研究［D］.南京:南京理工大學，2011.

［7］So R M C，Zhang H S，Speziale C G.Near-wall modeling of the dissipation rate equation［J］.AIAA Journal，1991，29(12):2069-2076.

［8］陶文銓.數值傳熱學(第2版)［M］.西安:西安交通大學出版社，2001.

［9］Chien K.Predictions of channel and boundary-layer flows with a low-Reynolds-number turbulence model［J］.AIAA Journal，1982，20(1):33-38.

［10］Hadar I，Gany A.Fuel regression mechanism in a solid fuel ramjet［J］.Propellants，Explosives，Pyrotechnics，1992，17(2):70-76.

［11］De Wilde J P.Fuel pyrolysis effects on hybrid rocket and solid fuel ramjet combustor performance［D］.Delft:Delft University of Technology，1991.

［12］王守范.固體火箭發動機燃燒與流動［M］.北京:北京工業學院出版社，1987.

［13］鄭亞，陳軍，鞠玉濤，等.固體火箭發動機傳熱學［M］.北京:北京航空航天大學出版社，2006.

［14］楊世銘，陶文銓.傳熱學(第四版)［M］.北京:高等教育出版社，2006.

［15］De Wilde J P.The heat of gasification of polyethylene and polymethylmethacrylate［M］.Delft:Delft University of Technology，1988.

Study of instantaneous regression rate in solid fuel ramjet based on dynamic mesh

WEI Tao，WU Xiao-song
(Department of Aeronautics and Astronautics，School of Mechanical Engineering，Nanjing University of Science and Technology，Nanjing 210094，China)

In order to study the changes in regression rate of solid grain in solid fuel ramjet(SFRJ)，a method to indicate the instantaneous regression rate of solid grain in SFRJ was developed based on the operational feature，dynamic mesh，as well as coupled simulation of gas-phase combustion，heat transfer and regression of solid grain.Then whole inner flow of an experimental motor with additional chamber and polyethylene(PE)fuel was numerically simulated.The time-dependent flow was obtained with dynamic boundaries，and the internal ballistic parameters were analyzed.The results show that combustion mainly occurs in the area where φ is between-2 and 2;as motor works，the regression rate decreases，the reattachment point moves downstream，and the velocity and temperature reduce at the outlet of solid grain;in addition，there is erosive burning in small SFRJ during initial operation stage.The study shows that this method can simulate the unsteady working process and fuel regression.The conclusions offer instruction for designing and experiment of solid fuel ramjet.

solid fuel ramjet;instantaneous regression rate;dynamic mesh;numerical simulation

V435

A

1006-2793(2012)04-0450-07

2011-06-10;

2012-02-17。

國防預研項目。

魏韜(1986—)，男，碩士，主要從事固體燃料沖壓發動機工作過程研究。E-mail:nustweitao@gmail.com

(編輯:崔賢彬)