基于旋轉閥的固體火箭發(fā)動機燃燒室壓強振蕩特性

席運志， 王寧飛， 李軍偉， 張智慧

(北京理工大學 宇航學院, 北京 100081)

0 引言

固體火箭發(fā)動機(SRM)出現的燃燒不穩(wěn)定，通常伴隨燃燒室內壓強不規(guī)則振蕩[1]。極端情況下，可能導致SRM在飛行過程熄火或爆炸，直接影響其總體飛行性能。因此，在SRM設計過程有必要對其內彈道特性和固體推進劑燃燒穩(wěn)定性進行評估[2-3]。

對于SRM內彈道特性，可通過流場仿真[4]或內彈道方程[5]計算獲取。然而，對于固體推進劑，因其具有組分復雜、燃燒環(huán)境惡劣等特點，目前尚無完備的理論模型可以精確預估其燃燒穩(wěn)定性。當前主要通過實驗的方式獲取壓強耦合響應函數，對固體推進劑的燃燒穩(wěn)定性進行預估[6-7]。在現有測試方案中[8]，由于T型燃燒器及其改進方案[9-12]具有結構簡單、操作方便等特點而被廣泛應用，但其也存在費用高、理論不完備、測試誤差高達30%～50%、不易開展低頻實驗等缺點[6]。為克服T型燃燒器存在的上述不足，文獻[13-15]提出旋轉閥法，該方法將小型SRM和旋轉閥裝置相結合，具有測試頻域廣、經濟適用性好、測試結果最接近實際SRM等優(yōu)點[6]，是T型燃燒器的一種良好替代或補充方案。

但旋轉閥法對測控系統精度要求極高，需精確測量燃燒室內壓強振蕩與旋轉閥有效排氣面積變化之間的相位差。該相位差直接決定實驗測試結果的準確性。為此，Brown等[13]采用接觸式測量方法直接獲取有效排氣面積變化，但存在磨損以及高頻工況下誤差大等缺點。由于缺少面積測量的相關改進方案，致使旋轉閥法目前仍未被廣泛使用，尤其是國內缺少對旋轉閥法的相關實驗及仿真研究。現階段SRM中出現的燃燒不穩(wěn)定多為低頻[16]，鑒于T型燃燒器對低頻測試的不足，亟需開展旋轉閥法相關研究，了解旋轉閥工作過程，優(yōu)化旋轉閥系統設計并提出準確相位角測控方案。

為進一步了解旋轉閥有效排氣面積變化與燃燒室壓強之間的關系，本文首先搭建一套基于旋轉閥的冷流實驗系統，并推導了有效排氣面積變化及燃燒室壓強理論計算模型。同時，基于動網格技術[17-19]和ANSYS/Fluent軟件平臺建立了旋轉閥工作過程的三維瞬態(tài)流場仿真模型(簡稱仿真模型)，該模型使用自定義函數(UDF)將旋轉閥轉動運動等效為周期性的擺動運動，用于開展內流場仿真與可視化研究。最后，通過實驗、理論模型及仿真模型之間的相互對比，驗證了理論模型和仿真模型的有效性。為研究旋轉閥壓強振蕩及流場特征提供了新的思路和方法，也進一步為熱流實驗中開展燃燒室壓強相對于有效排氣面積的相位角測量方案研究及旋轉閥改進設計工作奠定實驗基礎及理論驗證方法。

1 實驗方法和理論模型

1.1 實驗裝置

基于旋轉閥的冷流實驗裝置如圖1(a)所示，主要組成部分包含高壓氮氣瓶、燃燒室、旋轉閥和數據采集儀。實驗中，將高壓氣體經減壓閥降低至實驗需求壓強，而后經減壓的氣體進入燃燒室和旋轉閥。在減壓閥和燃燒室之間有壓強計和電磁閥，壓強計用于測量實際供氣壓強，電磁閥用于控制氣流的供給。其中，燃燒室內徑為72 mm，體積為320 cm3. 實驗主要流程為，先啟動旋轉閥并達到預定轉速，然后打開電磁閥開關，對燃燒室進行供氣，燃燒室內壓強變化由數據采集儀采集。圖1(b)中RED為轉子排氣通道，SED為定子排氣通道。

圖1 冷流實驗裝置和旋轉閥示意圖Fig.1 Schematic diagrams of cold-flow device and rotary valve

旋轉閥轉速范圍介于0～3 000 r/min，其內部結構剖視如圖1(b)所示，主要由轉子、轉子軸、定子、聯軸器、伺服電機等組成。其中，轉子與定子為石墨材質，轉子軸為中空軸，其材質為鋼。轉子和轉子軸裝配在一起形成轉動部件同步轉動。同時，轉子與伺服電機由聯軸器連接，伺服電機運轉由伺服執(zhí)行器控制，即轉子轉速實際由伺服執(zhí)行器控制。轉子外徑為76 mm，在轉子中心位置，沿其周向等間距開有16個半徑為1.5 mm的RED. 定子與燃燒室裝配在一起，燃燒室內高壓氣體可以通過定子中心位置的SED排氣。其中SED與RED位置中心對齊，為消除諧波組分[13]及安裝精度因素的影響，SED截面形狀設計為矩形，其長為4.5 mm，寬為3.0 mm. 在轉子轉動過程，16個RED與SED形成周期性的排氣通道，使得高壓氣體可從燃燒室經SED與RED周期性的排出，引發(fā)燃燒室內壓強振蕩，壓強振蕩頻率由轉子轉速決定。因此，通過控制供氣壓強和轉子轉速，可引發(fā)燃燒室內產生不同工作壓強和振蕩頻率的壓強振蕩。實驗中，供氣管路的壓強由杭州美控自動化技術有限公司產P300壓強計(量程0～10 MPa)測量，燃燒室內的壓強由西安杰誠傳感器測控技術有限公司產CYG4100高頻壓強傳感器(量程0～4 MPa)獲取。壓強傳感器采集的壓強數據及伺服電機轉動信號分別通過線路1、線路2傳輸給江蘇東華測試技術股份有限公司產DH5922D數據采集儀(16通道)。該采集儀兼具數據采集和分析的能力，最大采樣頻率200 kHz，實驗中使用采樣頻率為20 kHz.

1.2 理論模型

1.2.1 有效排氣區(qū)域面積求解模型

旋轉閥工作過程如圖2所示，轉子旋轉使得RED與SED之間形成周期性的排氣通道，燃燒室內高壓氣體從該通道排出引起壓強周期性振蕩。排氣通道的開- 閉由SED與RED交界面形成的有效排氣區(qū)域表征(見圖3)。當有效排氣區(qū)域面積Se大于0 mm2時，則表示排氣通道打開，反之關閉。為了解旋轉閥工作過程燃燒室內壓強振蕩與Se和轉速r之間的關系，本文對旋轉閥工作過程進行了理論建模。

圖2 旋轉閥轉動示意圖Fig.2 Motion of rotary valve

圖3 有效排氣區(qū)域隨時間變化示意圖Fig.3 Schematic diagram of effective exhaust region

(1)

式中：ω為角速度，ω=2πr，r為轉速(r/min)；Dr為轉子外徑(mm)。

為定量描述非排氣階段與排氣階段行程對燃燒室壓強振蕩的影響，引入行程比：

(2)

本實驗轉子上有16個RED，每個排氣周期對應行程比ΩL為1.485，則一個完整排氣周期時長滿足：

(3)

由圖2和圖3可知，Se隨轉子的轉動呈周期性變化。在一個排氣周期內，非排氣階段(A0→A1)，Se為0 mm2. 對于排氣階段(A1→A5)，Se變化關于A3位置對稱(見圖3)，其由0 mm2(A1)逐漸增加至1/2Se,max(A2，其中Se,max代表Se最大值)、Se,max(A2)而后逐漸減小至1/2Se,max(A4)、0 mm2(A5)，然后再次進入非排氣階段。下一個RED重復上述變化規(guī)律(A0至A5)進行周期性排氣。

圖4 有效排氣面積求解示意圖(排氣階段)Fig.4 Schematic diagram of solution of effective exhaust area (exhaust phase)

在一個完整排氣周期內，Se直接影響排氣量，進而影響燃燒室內壓強波動幅值，其大小由有效排氣區(qū)域的位置決定。Se的求解示意如圖4所示，其中Δx=vr(t-tA1)為RED在排氣階段轉過的弧線距離，tA1及tA5(同圖3)分別為排氣階段開始和結束的時刻。

在Δx介于0～2R0，即tA1

(4)

θcl=arcos(1-Δx/R0).

(5)

因為Se變化關于Δx=2R0位置對稱，因此，在Δx介于2R0～4R0，即tA3

(6)

θcr=arcos(3-Δx/R0).

(7)

本模型中，轉子上有16個RED均勻的分布在轉子圓周上，因此Se呈周期性變化，基于(1)式～(7)式，可以推導旋轉閥轉動過程Se隨時間變化規(guī)律：

(8)

(9)

式中：n為排氣周期數。

1.2.2 壓強振蕩求解模型

根據質量守恒定律，可建立旋轉閥轉動過程Se與壓強振蕩之間的理論計算模型。燃燒室中的氣體動態(tài)質量變化由(10)式定義：

(10)

根據對2017屆高三8個應屆文科班實施了“201010”課堂模式的339人進行問卷調查，調查問卷發(fā)放339份，收回339份，回收率100%，全部有效。

(11)

(12)

Cd為排氣流量修正系數[21-22]，

(13)

γ為氣體比熱比。

結合理想氣體狀態(tài)方程pV=mRgTg，將(11)式～(13)式代入(10)式，得

(14)

式中：V為燃燒室自由容積。

基于(14)式可計算旋轉閥轉動過程中轉子轉速r、行程比ΩL、RED半徑R0及供氣壓強pin等參數對燃燒室壓強振蕩規(guī)律的影響。冷流實驗中主要參數如表1所示。

考慮到通過實驗方式開展不同參數對燃燒室壓強振蕩影響規(guī)律研究，具有成本高、操作復雜及獲取有效數據點有限等缺點，以及理論模型存在無法獲取旋轉閥內部瞬時流動細節(jié)、流場結構特征等不足。為進一步解決實驗和理論模型存在的上述問題，本文同時基于旋轉閥實驗裝置(見圖1)建立了三維瞬態(tài)流場仿真模型，該模型能夠提供更加全面的三維流場可視化信息及流動細節(jié)，便于開展旋轉閥轉動過程壓強振蕩特性研究。

表1 冷流實驗中主要參數

2 仿真模型

2.1 物理模型

由圖2可知，當旋轉閥工作時，在一個排氣周期內，轉子上的16個RED僅有一個進行排氣, 其余處于非排氣狀態(tài)。如果仿真模型對16個RED同時進行計算，則計算量非常大。為了縮短計算時間，降低計算費用，本文將16個RED單向旋轉運動等效為一個RED左右擺動運動，其擺動模型如圖5所示，其中Ⅰ、Ⅲ分別表示RED在擺動過程的左右極限位置，Ⅱ為SED與RED中心對正時的位置。

圖5 RED擺動模型示意圖Fig.5 Schematic diagram of swing model of RED

在擺動模型中，RED擺動周期及線速度大小分別為T和vr與旋轉運動一致，其左右極限位置與SED中心線夾角θs均為vrT/2Dr. 擺動過程，RED依次經過位置Ⅰ→位置Ⅱ(見圖5(a))→位置Ⅲ(見圖5(b))向左擺動完成一個排氣周期，而后向右擺動依次經過位置Ⅲ→位置Ⅱ(見圖5(c))→位置Ⅰ(見圖5(d))，進入下一個排氣周期，循環(huán)往復，實現周期性排氣。因此，該擺動模型可以模擬多孔轉子轉動引起的周期性排氣過程。

2.2 網格生成及邊界條件

使用Solidworks軟件構建旋轉閥內部流體域的三維幾何模型，而后將其導入ANSYS/ICEM軟件對流體域進行網格劃分。該流體域主要包含3部分：燃燒室域(含SED)、環(huán)境域及RED域，仿真計算中前兩者處于靜止狀態(tài)，而RED按照擺動模型進行運動。由于涉及RED擺動域，且RED在擺動過程與SED存在非接觸時段(見圖5(a))。因此，三者網格需獨立劃分，最后通過網格融合技術[17]進行裝配。

對于流場仿真，計算域網格質量是準確獲取計算結果的關鍵因素[17,23]，而網格模型劃分類型決定了網格劃分難度及仿真計算時間。針對燃燒室，由于其幾何結構相對復雜，本文主要采用四面體非結構網格對其進行網格劃分，并在近壁面添加3層邊界層，邊界層第1層網格寬度為0.2 mm，網格總數約50萬。

相對于燃燒室尺寸，其壓強入口Inlet及SED的尺寸較小，為保證計算精度同時減少網格數量，對上述兩區(qū)域附近的網格進行局部加密處理，如圖6(a)所示。由于RED處于擺動狀態(tài)，最大速度和壓強梯度將出現在有效排氣區(qū)域(重疊面1)附近，同時考慮到非結構網格具有無序性，為便于后期對附近流場數據進行分析和處理，將SED下部與RED接觸的部分采用六面體結構網格進行劃分，并對結構網格和非結構網格交界面進行節(jié)點對齊保證數據傳遞的連續(xù)性，如圖6(b)所示，圖中YL為RED的無量綱長度。

圖6 計算域三維網格劃分Fig.6 Topology of 3D mesh

圖7 有效排氣區(qū)域網格示意圖Fig.7 Mesh of effective exhaust region

圖8 RED三維運動(局部)示意圖(從位置Ⅲ到Ⅱ)Fig.8 Mesh position of RED in swing process (from Ⅲ to Ⅱ)

將生成的計算域網格導入ANSYS/Fluent軟件平臺中對圖6所示邊界進行定義，具體如下：

1) 進氣口: 此邊界設置為壓強入口，與實驗裝置的供氣管路入口位置一致。其數值大小等于經減壓閥后的高壓氣體壓強，方向垂直于邊界。

2) 壓強出口：此邊界設置為壓強出口，出口壓強與大氣環(huán)境一致為0.1 MPa，其余流場參數由空間內部迭代計算獲取。

3) 有效排氣區(qū)域(重疊面1，EER)：定義RED上表面及SED下表面邊界類型為干涉面(見圖7)。流場仿真過程采用動網格技術[17,24]和自UDF[18]相結合的方式實現RED左右擺動(見圖5)。當RED處于排氣階段時，RED上表面與SED下表面局部重疊形成有效排氣區(qū)域，重疊部分自動識別為內部邊界，非重疊部分自動識別為壁面[17, 19]，即來自燃燒室內的高壓氣流可從SED經有效排氣區(qū)域流入RED，但不能經非重疊區(qū)域的干涉面進入RED. 同樣的，當RED處于非排氣階段時, RED上表面和SED小表面不發(fā)生重疊，均自動識別為壁面。同理，重疊面2定義同重疊面1.

4)壁面: 除了上述邊界條件外，其余表面均設定為無滑移壁面。

2.3 湍流模型

在旋轉閥工作過程，RED處于周期性的擺動狀態(tài)，來自燃燒室的高壓氣流在RED內部以自由剪切流和受壁面效應強烈影響的壁面邊界流[24]并存的形式流動，且在有效排氣區(qū)域附近的流場壓強及馬赫數變化最劇烈。為準確獲取流場變化，需根據流場雷諾數確定流動類型。考慮到旋轉閥的幾何結構及流體介質的物理參數，可獲取雷諾數Re的變化范圍：從非排氣階段的最小值0至排氣階段的最大值16 000，即部分階段的流動類型為湍流流動。

本文流場仿真基于ANSYS/Fluent軟件平臺展開，其提供了多種可選擇的湍流模型[17]。根據旋轉閥的工作特點以及前人的相關研究[17,19-20,24]，本文采用的湍流模型為RNGk-ε模型，該模型可給出有效排氣區(qū)域兩側的湍流量預估值，并且能夠正確估計高壓自由射流和壁面的邊界區(qū)域[25]，這對確定RED內部流動特性具有重要意義。與此同時，理想的可壓縮氣體模型、能量方程及基于壓力基的求解器也一并被采用。對于控制方程，本文采用有限體積法對其進行離散化處理。對于擴散相，采用一階差分法；對于對流、湍流動能及湍流耗散率則采用一階迎風格式進行處理。基于上述求解方法，同時設定了以下收斂標準：除了能量殘差的絕對值等于10-6外，其余均為10-4. 對于瞬態(tài)流動，時間步長與轉子轉速相關，其變化范圍為10～2 μs.

針對瞬態(tài)流場計算中的RED擺動問題，最合適的策略為UDF和動網格技術相結合的方案[17-19]。RED擺動過程其附近及內部流場參數可由(15)式求解：

(15)

式中：ρ為密度；u為流場速度矢量；ug為動網格運動矢量；Γ為擴散系數；Sφ為源相。

在仿真計算中忽略次要因素的影響，譬如RED與SED之間的漏氣、轉子轉速誤差等。為進一步獲取流場內部數據，在計算域的不同位置分別設置了壓強監(jiān)控點p1、p2、p3及質量流率監(jiān)控面：進氣口、有效排氣區(qū)域,如圖6所示，其中監(jiān)控點p1位置和實驗測壓點一致。

3 仿真模型有效性驗證

為節(jié)省流場仿真計算時間，本文采用擺動模型(見圖5)代替單向旋轉運動。為驗證簡化方法可行性以及仿真模型的準確性，首先采用仿真模型對轉子周期性擺動進行了計算，并將擺動運動中有效排氣面積變化與理論模型相比較。其次，采用冷流實驗和理論計算方法獲得了旋轉閥轉動引起的燃燒室內壓強振蕩，并與流場仿真結果進行對比。

本節(jié)有效性驗證主要開展了3種排氣頻率f分別為20 Hz、100 Hz及300 Hz的驗證實驗、理論及流場仿真計算，具體如表2所示。除表2中參數外，其余條件相同：供氣壓強pin為2.15 MPa，流體介質為氮氣，RED半徑為1.5 mm.

表2 有效性驗證工況

圖9 RED周期性擺動的UDF邏輯框圖Fig.9 UDF logic diagram for periodic swing of RED

圖10 有效排氣面積及質量流率隨時間變化Fig.10 Variations in mass flow rate and effective exhaust area with time

對于仿真模型的有效性驗證，本文針對工況1～工況3分別開展了冷流實驗和理論計算，獲取了燃燒室內p1位置的壓強振蕩數據并與仿真結果對比，如圖11(a)所示，其中pmax和pmin分別為壓強最大值和最小值。由圖11(a)可知，3種方法獲得的壓強曲線同步變化，均由壓強最小值上升至供氣壓強，而后維持不變，最后再減小，呈周期性變化。其中仿真計算中監(jiān)控點p1～p3處的壓強數值及變化規(guī)律一致，即燃燒室內壓強處于整體振蕩狀態(tài)。因此，后續(xù)將采用監(jiān)控點p1處數據開展壓強振蕩規(guī)律研究。

對于工況2和工況3，理論模型和仿真計算獲取的燃燒室內壓強最大值均無法達到供氣壓強2.15 MPa，說明轉速(或排氣頻率)影響壓強振蕩幅值。同時，由于實驗中采用減壓閥手動控制供氣壓強，因而存在誤差，比理論計算和瞬態(tài)流場仿真供氣壓強高約0.01 MPa. 進一步對圖11數據進行處理，獲取燃燒室壓強振蕩表征參數，如表3所示，ε為壓強振蕩幅值比，

圖11 燃燒室內壓強隨時間變化曲線Fig.11 Pressure oscillation in combustor vs. time

(16)

式中：Δp為壓強峰- 峰值，Δp=pmax-pmin.

表3 壓強振蕩表征參數

由表3可知，對于工況1，3種方法對應燃燒室最大壓強pmax均能達到供氣壓強，實驗排氣周期T比理論和仿真計算多0.15 ms，誤差小于1%. 同時，理論和仿真計算的壓強峰- 峰值Δp較接近約0.04 MPa但比實驗值0.058 MPa小。對于工況2和工況3，3種方法對應燃燒室最大壓強pmax均未能達到供氣壓強，兩種工況的實驗排氣周期T比理論和仿真計算分別多0.11 ms和少0.01 ms，誤差約1%. 同時，3種方法獲取的壓強峰- 峰值Δp基本一致，分別約為0.01 MPa和0.005 MPa. 由上分析可知，本文建立的瞬態(tài)流場仿真模型是有效的，可利用其進一步開展不同工作條件下燃燒室壓強振蕩特性研究。

4 旋轉閥工作過程仿真分析

本節(jié)以工況2為例, 對旋轉閥工作過程進行瞬態(tài)流場仿真，獲取了燃燒室內壓強、質量流率變化曲線以及典型時刻的流場馬赫數、壓強分布云圖。

圖12 質量流率及燃燒室壓強曲線Fig.12 Variations in mass flow rate and pressure with time

在旋轉閥工作期間RED位置不斷變化，有效排氣區(qū)域附近的流場參數變化劇烈，進一步獲取圖12中典型時刻t=46.3 ms、47.98 ms以及有效排氣區(qū)域面積為1/2Se,max對應時刻46.97 ms的流場馬赫數Ma及壓強分布云圖，如圖13所示。在t=46.3 ms時，RED剛進入非排氣階段，其與SED下表面部分重疊，高壓氣流在RED內部形成超音速氣流，最大馬赫數約3.2. 隨著有效排氣區(qū)域面積的增加，在t=46.97 ms時，SED內壓強減小，RED中最大馬赫數降為2.6且在RED內部形成旋渦。而到了t=47.98 ms時，有效排氣區(qū)域面積達到最大值，此時凈排氣質量也達到了最大值，SED內部壓強繼續(xù)減小，RED最大馬赫數降為1.2.

圖13 不同時刻的RED和SED內部流場Fig.13 Flow fields inside RED and SED at different times

5 結論

本文通過搭建旋轉閥冷氣實驗系統，推導燃燒室壓強振蕩理論計算模型，構建旋轉閥三維瞬態(tài)內流場計算模型，研究了旋轉閥有效排氣面積變化與燃燒室壓強之間的關系，分析了旋轉閥內部瞬態(tài)流動特征，得到如下結論：

1)燃燒室壓強呈周期性變化，振蕩頻率由旋轉閥排氣頻率決定，且排氣頻率越高壓強振蕩幅值比越小，逐漸由20 Hz的2.68%降至300 Hz的0.23%.

2)燃燒室壓強振蕩相對于有效排氣面積變化存在延遲，延遲時間由燃燒室凈質量流率決定。

3)RED內部流場馬赫數劇烈變化程度隨有效排氣區(qū)域面積增加而減弱，其最大值由初始階段的3.2變?yōu)樽畲竺娣e時的1.2.