王建儒，何國強(qiáng)，許團(tuán)委，田維平

(1.西北工業(yè)大學(xué)航天學(xué)院，西安 710072;2.中國航天科技集團(tuán)公司四院四十一所，西安 710025)

燃燒室對接狹縫設(shè)計參數(shù)對壓強(qiáng)振蕩的影響研究①

王建儒1，何國強(qiáng)1，許團(tuán)委2，田維平2

(1.西北工業(yè)大學(xué)航天學(xué)院，西安 710072;2.中國航天科技集團(tuán)公司四院四十一所，西安 710025)

為了獲得分段式固體發(fā)動機(jī)藥柱對接部位狹縫寬度對發(fā)動機(jī)壓強(qiáng)振蕩的影響規(guī)律，針對某典型分段式固體發(fā)動機(jī)建立了二維軸對稱模型。采用大渦模擬(LES)方法，以“突擴(kuò)比”為重要設(shè)計參數(shù)，完成了6種不同“突擴(kuò)比”設(shè)計條件下的分段式發(fā)動機(jī)流場穩(wěn)定性數(shù)值模擬，并分別對其壓強(qiáng)-時間曲線進(jìn)行了FFT分析。結(jié)果表明，隨著“突擴(kuò)比”的逐漸增大，燃燒室壓強(qiáng)振蕩頻率基本不變，但壓強(qiáng)振蕩的幅值總體上呈現(xiàn)下降的趨勢。所得分析結(jié)果已得到試驗的初步驗證。

分段式固體發(fā)動機(jī);壓強(qiáng)振蕩;大渦模擬;突擴(kuò)比

0 引言

分段式固體發(fā)動機(jī)是實現(xiàn)大推力的有效技術(shù)途徑之一，國外紛紛采用該技術(shù)實現(xiàn)推力百噸級、甚至千噸級的大型固體發(fā)動機(jī)的應(yīng)用［1］。分段式固體發(fā)動機(jī)燃燒室設(shè)計結(jié)構(gòu)的獨有特點對工作過程中的內(nèi)彈道穩(wěn)定性帶來了一定影響。由于燃燒室分段對接部位不可避免地需要設(shè)計狹縫，這種設(shè)計結(jié)構(gòu)導(dǎo)致燃燒室內(nèi)流道幾何構(gòu)型趨于復(fù)雜［2］，在發(fā)動機(jī)工作過程中會對燃?xì)饬鲃拥姆€(wěn)定性帶來一定影響［3-4］。特別是在點火初期(本文指0.5～2 s內(nèi))，由于燃?xì)馔ǖ涝趯硬课坏耐蝗蛔兓?，?dǎo)致壓力產(chǎn)生一定的變化，再加之對接結(jié)構(gòu)帶來的各種渦脫落的綜合影響，在很大程度上形成了燃燒室壓力在一定范圍內(nèi)的波動。

國外在分段式固體發(fā)動機(jī)流場穩(wěn)定性技術(shù)領(lǐng)域給予了高度關(guān)注。針對多個發(fā)動機(jī)開展了流動不穩(wěn)定性的影響分析研究［5-8］，以法國為主的歐洲啟動兩大計劃研究壓強(qiáng)振蕩問題，分別是分段式固體發(fā)動機(jī)氣體動力學(xué)(ASSM，Aerodynamics of Segmented Solid Motors)和壓強(qiáng)振蕩計劃(POP，Pressure Osillations Program)。ASSM計劃的主要目標(biāo)是對渦脫落深入理解和建模，促進(jìn)數(shù)值模擬技術(shù)的發(fā)展。

POP計劃是利用P230的縮比模型發(fā)動機(jī)開展實驗研究，獲得實驗和數(shù)值數(shù)據(jù)庫。美國對此也投入了大量的人力和經(jīng)費，開展了多學(xué)科大學(xué)研究倡議(MURI，Multi-Disciplinary University Research Initiative)，試圖從基礎(chǔ)化學(xué)、燃燒和流體動力學(xué)的角度深入研究燃燒不穩(wěn)定課題［9］。但大多停留在已有結(jié)構(gòu)的驗證分析方面，對于對接狹縫的結(jié)構(gòu)尺寸變化對流動穩(wěn)定性的影響未見報道。隨著國內(nèi)針對大型分段式固體發(fā)動機(jī)技術(shù)研究的不斷深入，開展分段對接狹縫寬度對壓強(qiáng)振蕩的影響研究至關(guān)重要。

本文針對某典型分段式發(fā)動機(jī)模型，建立了內(nèi)流場數(shù)學(xué)模型，采用LES方法，針對發(fā)動機(jī)不同“突擴(kuò)比”的設(shè)計結(jié)構(gòu)，獲得了關(guān)鍵設(shè)計參數(shù)對于分段式發(fā)動機(jī)壓力振蕩特性的影響規(guī)律，并將研究成果用于試驗發(fā)動機(jī)，獲得了具有一定工程應(yīng)用價值的成果。

1 物理模型及工況參數(shù)

本文重點研究分段對接狹縫設(shè)計參數(shù)對于壓強(qiáng)振蕩頻率與振幅特性的影響規(guī)律。因此，在模型選擇上，選取某典型(φ400 mm)二分段式固體發(fā)動機(jī)開展研究工作。計算區(qū)域如圖1所示。從多種工況進(jìn)行計算比較以及節(jié)省計算資源等方面考慮，進(jìn)行了二維軸對稱簡化處理，如圖2所示。

圖1 直徑400 mm/2分段式發(fā)動機(jī)計算區(qū)域Fig.1 Calculation region for a two segmented solid rocket motor(φ400 mm)

圖2 二維軸對稱計算區(qū)域圖Fig.2 Calculation region for a 2D axisymmetric model

分段式發(fā)動機(jī)分段對接結(jié)構(gòu)主要由狹縫的寬度、深度和燃?xì)馔ǖ赖热笠蛩亟M成。本文重點研究點火初期的的壓力振蕩。因此，對于狹縫深度的影響暫不考慮。在保持裝藥內(nèi)孔半徑不變下，通過改變分段狹縫寬度的值(即提出“突擴(kuò)比”的概念)，便可得到一組同一藥柱內(nèi)孔、不同狹縫寬度設(shè)計的分段式發(fā)動機(jī)模型。分別建立模型，按照同一計算方法，依次開展不同工況下的大渦模擬研究，所取工況如表1所示。

表1 分段對接狹縫寬度取值Table 1 Width values of joint gap in SSRM

表2 計算工況參數(shù)Table 2 Calculation parameters

2 流場數(shù)學(xué)模型建立與數(shù)值算法

2.1 控制方程

針對氣動聲學(xué)、燃燒室中湍流動量輸運(yùn)及分離流或渦流區(qū)流動等問題，LES以其耗散小、精度高等特點比RANS具有更大優(yōu)勢。因此，可利用LES來直接計算壓強(qiáng)振蕩引起的聲波運(yùn)動。本文的研究涉及的正是分離流、旋渦脫落、壓強(qiáng)振蕩及聲學(xué)。LES采用過濾器對N-S方程進(jìn)行過濾，將小尺度脈動過濾掉，得到亞格子應(yīng)力項。然后，直接計算大尺度脈動，而使用亞格子模型作為亞格子應(yīng)力的封閉模式，代表小尺度脈動的貢獻(xiàn)。因此，亞格子模型的選擇至為重要。本文利用空間濾波器G，將流場變量分為可解尺度脈動量與亞格子尺度脈動量:

考慮到氣體的可壓縮性，利用Favre平均對控制方程進(jìn)行簡化:

式中 “-”表示Reynold平均;“～”表示Favre平均。

本文不考慮化學(xué)反應(yīng)，僅計算單組分工質(zhì)，濾波后連續(xù)方程、動量方程與能量方程分別為

式中μ為動力粘性系數(shù);cp為定壓比熱容;普朗特數(shù)Pr=0.75。

2.2 亞格子模型選擇及控制方程離散

在具有后向臺階的發(fā)動機(jī)中，剪切層及流動自身的不穩(wěn)定性引起了表面渦脫落及轉(zhuǎn)角渦脫落，而剪切層與近壁區(qū)的流動有直接關(guān)系。由于Smagorinsky模型無法準(zhǔn)確預(yù)測近壁區(qū)流動及湍流轉(zhuǎn)捩現(xiàn)象，WALE(Wall-Adaping Local Eddy-viscosity)亞格子模型可對近壁區(qū)進(jìn)行修正。因此，本文選擇WALE亞格子模型，對亞格子應(yīng)力張量、亞格子通量張量以及亞格子尺度粘性力變形功進(jìn)行封閉。

本文共邀請6名智慧城市建設(shè)領(lǐng)域相關(guān)專家進(jìn)行打分，采用百分制，分值越大，表示該指標(biāo)的比重越大.具體打分結(jié)果如表1所示.

對于連續(xù)方程與動量方程，為了避免中心差分格式容易產(chǎn)生數(shù)值振蕩，采用 BCD(Bounded central differencing)格式進(jìn)行離散，該格式可發(fā)現(xiàn)求解區(qū)域中波長小于2Δx的波動，并加以抑制;能量方程則采用Power Law格式，以加速收斂。同時采用二階隱格式，計算步長均為10-5，庫朗特數(shù)在計算過程中，根據(jù)收斂穩(wěn)定性與速度進(jìn)行調(diào)整。為提高計算的收斂性，先采用空間二階離散格式的RANS方法得到流場的穩(wěn)態(tài)解;然后，再調(diào)用LES方法求解整個非定常流動。

2.3 計算方法驗證

鑒于國外對于VKI(Von Karman Institute for Fluid Dynamics)冷流實驗?zāi)Ｐ?，從實驗和?shù)值模擬方面已做了大量工作，并公開了模型的幾何尺寸、實驗結(jié)果和計算結(jié)果［10-11］。本節(jié)亦選用VKI模型進(jìn)行數(shù)值計算，計算區(qū)域如圖3所示。模型發(fā)動機(jī)為二維軸對稱結(jié)構(gòu)，含有絕熱環(huán)及潛入式噴管。

圖3 模型發(fā)動機(jī)計算區(qū)域Fig.3 Computational model and domain

圖4為模型發(fā)動機(jī)頭部壓強(qiáng)隨時間變化曲線及其FFT分析結(jié)果。通過統(tǒng)計處理，可得出壓強(qiáng)振蕩平均值為0.190 MPa，第一階和第二階振蕩頻率分別為397 Hz和786 Hz，相對應(yīng)的振幅與平均壓強(qiáng)的比率分別為 0.818%和 0.128%。

圖4 模型發(fā)動機(jī)頭部壓強(qiáng)隨時間變化曲線及其FFT分析結(jié)果Fig.4 Evolution of the pressure on head point in mode motor and its FFT

由表3可看出，振頻數(shù)值計算結(jié)果與實驗結(jié)果較吻合，除了第二階振蕩頻率數(shù)值計算結(jié)果與實驗結(jié)果誤差相對偏大，最大誤差約為8.98%;通過振幅對比，數(shù)值結(jié)果與實驗結(jié)果相差約一個數(shù)量級。分析認(rèn)為，由于實驗采用多孔介質(zhì)進(jìn)氣，聲波可穿透此處繼續(xù)向上游傳播，從而帶來一定的聲能耗散。此外，數(shù)值計算中沒有考慮壁面阻尼、噴管阻尼以及漩渦的三維特性。因此，數(shù)值計算結(jié)果與實驗所測振幅差異較大。但數(shù)值計算結(jié)果的規(guī)律性與實驗結(jié)果較接近，兩者之間的誤差不影響對流動不穩(wěn)定的規(guī)律性探討。

綜上所述，本文采用的數(shù)值計算方法能較準(zhǔn)確地獲取燃燒室內(nèi)的壓強(qiáng)振蕩頻率與振幅特性，為發(fā)動機(jī)壓強(qiáng)振蕩研究提供了一種很好的有效數(shù)值計算方法。

表3 計算結(jié)果與實驗結(jié)果對比分析Table 3 Comparable analysis for calculation and experimental data

3 計算結(jié)果分析

首先開展了穩(wěn)態(tài)流場計算，在此結(jié)果基礎(chǔ)上開展了純氣相的LES計算。在對模擬結(jié)果的處理上，分別從宏觀方面和微觀方面進(jìn)行分析。從宏觀角度，分析了結(jié)合渦旋分布云圖，分析了流場的整體穩(wěn)定情況，總結(jié)了旋渦運(yùn)動特點;從微觀角度，在計算區(qū)域布置了4個監(jiān)測點，依次位于發(fā)動機(jī)頭部、狹縫處、狹縫下游圓柱段中部以及裝藥后錐段中部，對壓強(qiáng)隨時間變化的過程進(jìn)行了監(jiān)測，并分別對其進(jìn)行了FFT分析。

3.1 渦旋強(qiáng)度分布云圖及分析

圖5為某兩分段式發(fā)動機(jī)在6種不同“突擴(kuò)比”設(shè)計條件下在點火初始時刻的內(nèi)流場旋渦強(qiáng)度分布云圖。從圖5中可看出，隨著突擴(kuò)比由0.61逐漸增大到1.63，狹縫空腔所容納的旋渦數(shù)量也由多逐漸變少，旋渦強(qiáng)度逐漸減弱，流場逐漸趨于平穩(wěn)。

分析認(rèn)為，當(dāng)突擴(kuò)比較小時，即對接狹縫寬度h較大，發(fā)動機(jī)第一分段端面處近似于90°的轉(zhuǎn)角所形成的轉(zhuǎn)角渦脫落之后，進(jìn)入到了狹縫空腔里，在凹腔里形成了紊亂的高低壓強(qiáng)交替，使得發(fā)動機(jī)內(nèi)部整體振幅比增大。隨著渦脫落的繼續(xù)生成和填充空腔，迫使先進(jìn)入凹腔的渦流出狹縫，與發(fā)動機(jī)第二分段燃面處所形成的表面渦脫落發(fā)生連續(xù)的碰撞與聚合等過程，形成尺寸較大的渦旋，沿流向運(yùn)動尺寸較大的渦旋對剪切層的分離產(chǎn)生明顯的促進(jìn)作用，使渦旋的尺寸和強(qiáng)度都進(jìn)一步得到了放大，同時渦旋的結(jié)構(gòu)形狀也發(fā)生了明顯的改變;當(dāng)狹縫寬度h逐漸變小時，此時的凹腔結(jié)構(gòu)為狹長型，即寬度小而深度大，分離形成的渦旋很少，或者無法“侵入”凹腔內(nèi)部，越過凹腔繼續(xù)向下游傳播，在運(yùn)動過程中被耗散掉了。

3.2 監(jiān)測點壓強(qiáng)隨時間變化曲線及分析

經(jīng)過對比分析，模擬發(fā)動機(jī)各點處所獲得的壓強(qiáng)隨時間變化曲線幾乎一致。因此，研究中主要以發(fā)動機(jī)頭部的監(jiān)測點為例，對不同工況下該點處的壓強(qiáng)隨時間變化曲線進(jìn)行FFT分析，以獲得壓強(qiáng)振蕩的頻率與振幅特性，見圖6。

圖5 發(fā)動機(jī)內(nèi)流場旋渦強(qiáng)度分布云圖Fig.5 Vorticity distribution in SSRM(No.6)

圖6 頭部壓強(qiáng)隨時間變化曲線及FFT分析結(jié)果Fig.6 Evolution of the pressure on head point in chamber and its FFT

在點火初始時刻不同突擴(kuò)比設(shè)計條件下前兩階壓強(qiáng)振蕩頻率大小對比見圖7(a)，不同突擴(kuò)比設(shè)計條件下前兩階振蕩頻率對應(yīng)的振幅大小對比見圖7(b)。數(shù)值模擬結(jié)果表明，隨著狹縫寬度的減小，前兩階頻率的降低幅度或者增大幅度都很小，一階呈現(xiàn)微弱下降趨勢，二階頻率呈現(xiàn)微弱上升趨勢;前兩階振蕩頻率所對應(yīng)的振幅隨著狹縫寬度的減小均有降低的趨勢。

圖7 不同突擴(kuò)比設(shè)計條件下前兩階壓強(qiáng)振蕩頻率及其對應(yīng)的振幅比Fig.7 Frequency and amplitude ratio of pressure oscillation for different sudden expand ratio

通過分析認(rèn)為，第一階頻率主要是對接結(jié)構(gòu)部位的轉(zhuǎn)角渦脫落頻率占主導(dǎo)。由于在不同突擴(kuò)比(即藥柱內(nèi)孔不變，狹縫尺寸逐漸減小)條件下，渦脫落的產(chǎn)生形式、機(jī)理沒有變化，只是渦與渦之間的相互作用時間、大小發(fā)生了一定的變化，因此其一階頻率變化相對穩(wěn)定;第二階頻率主要是脫落渦撞擊噴管的頻率占主導(dǎo)，其值與主流速度及渦脫落發(fā)生點與在噴管上的撞擊點之間距離有關(guān)，而渦進(jìn)入狹縫無形中增加了渦的運(yùn)動距離，即增加了渦從脫落到撞擊噴管的時間。也正是上述原因，導(dǎo)致隨著對接狹縫間隙的不斷減小，在一定程度上增加了渦脫落和噴管的撞擊頻率，使得二階頻率呈現(xiàn)小范圍的增加趨勢。因此，分段式“突擴(kuò)比”的變化對壓強(qiáng)振蕩的振幅影響較大。在此結(jié)論下，針對分段式發(fā)動機(jī)可提出在進(jìn)行初始裝藥設(shè)計時，分段對接的“突擴(kuò)”應(yīng)設(shè)計盡可能大，具有抑制點火初期壓力波動的可能性。

4 試驗驗證

為了進(jìn)一步驗證分段式發(fā)動機(jī)對接結(jié)構(gòu)中“突擴(kuò)比”設(shè)計參數(shù)對發(fā)動機(jī)點火初期壓力振蕩的影響分析結(jié)果，設(shè)計完成了2發(fā)不同突擴(kuò)比的試驗發(fā)動機(jī)進(jìn)行熱試驗證，2臺發(fā)動機(jī)分別設(shè)計突擴(kuò)比為0.8和1.2，如圖8所示。

突擴(kuò)比為0.8的試驗發(fā)動機(jī)在點火初期3 s以內(nèi)，其壓力振幅較大;而突擴(kuò)比1.2的驗證發(fā)動機(jī)在點火初期，幾乎未發(fā)現(xiàn)有明顯的壓力振蕩。由此可見，通過合理設(shè)計對接部位的突擴(kuò)比，是有效抑制發(fā)動機(jī)工作初期的壓力脈動的一條途徑。

圖8 驗證發(fā)動機(jī)實測壓力Fig.8 Pressure evolution(0 ～3 s)for a test SSRM

5 結(jié)論

(1)本文所采用的大渦模擬(LES)方法能較準(zhǔn)確地獲取固體火箭發(fā)動機(jī)燃燒室內(nèi)的壓強(qiáng)振蕩頻率與振幅特性，是一種合理有效的數(shù)值計算方法。

(2)隨著裝藥初始突擴(kuò)比的增大，前兩階頻率的降低幅度或者增大幅度都很小，一階呈現(xiàn)微弱下降趨勢，二階頻率呈現(xiàn)微弱上升趨勢;而前兩階振蕩頻率所對應(yīng)的振幅隨著突擴(kuò)比的增大均有降低趨勢。

(3)在分段式發(fā)動機(jī)設(shè)計時，通過對裝藥初始突擴(kuò)比進(jìn)行優(yōu)化選擇，有利于發(fā)動機(jī)工作工程初期的穩(wěn)定，這一結(jié)論已得到試驗的初步驗證。

［1］葉定友，高波，等.重型運(yùn)載火箭大型固體助推器技術(shù)研究［J］.載人航天，2011(1).

［2］John R McGuire.Pressure-actuated joint system［P］.U S Patent＆ Trademark Office，US 6711890B2，2004.

［3］Dotson K，Koshingoe S，Pace K.Vortex driven pressure oscillations in the titan IV solid rocket motor upgrade［C］//31st AIAA/ASME/SAE/ASEE/Joint Propulsio Conference and Exhibit.AIAA Paper 1995-2732.

［4］Fred S Blomshield.Lessons learned in solid rocket combustion instability［C］//43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference ＆ Exhibit.AIAA 2007-5803.

［5］Avalon G，Josset Th.Cold gas experimental applied to the understanding of aeroacoustic phenomena inside solid propellant boosters［R］.AIAA 2006-5111.

［6］Stany Gallier，Michel Prevost，Jouke Hijlkema，et al.Effects of cavity on thrust oscillations in subscale solid rocket motors［C］//45stAIAA /ASME/SAE /ASEE Joint Propulsion Conference and Exhibit.AIAA 2009-5253.

［7］Prevost M.Thrust oscillations in reduced scale solid rocket motors，partⅠ:Experimental Investigations［R］.AIAA 2005-4003.

［8］Prévost M，Le A Quellec，Godon J C.Thrust oscillations in reduced scale solid rocket motors，a new configuration for the MPS of ariane5［C］//42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference＆ Exhibit.AIAA 2006-4418.

［9］Fred S Blomshied.Summary of multi-disciplinary university research initiative in solid propellant combustion instability［R］.AIAA 2000-3172.

［10］Stella F，Paglia F.Pressure oscillations in solid rocket motors:numerical study［J］.Aerospace Science and Technology，2011，15(1).

［11］Anthoine J J，Buchlin M，Guery J F.Effect of nozzle cavity on resonance in large SRM:numerical simulations［J］.Journal of Propulsion and Power，2003，19(3).

Effect of internal joint design parameters on pressure oscillation in combustion chamber

WANG Jian-ru1，HE Guo-qiang1，XU Tuan-wei2，TIAN Wei-ping2
(1.College of Astronautics，Northwestern Polytechnical University，Xi＇an 710072，China;2.The 41st Institute of the Fourth Academy of CSAC，Xi＇an 710025，China)

Aiming at revealing the relation between the internal joint width and the pressure oscillation for SSRM(Segmented Solid Rocket Motor)，a two-dimensional axisymmetric model was built.LES(Large Eddy Simulation)was carried out for six different expansion ratio.Finally，pressure evolution and its FFT(Fast Fourier Transform)results were analyzed.Results indicate that the amplitude of pressure oscillation displays a decrease tendency as the expansion ratio increasing，but the frequency of pressure oscillation is unchanged.The result has gained validation through SSRM static test.

segmented solid rocket motor;pressure oscillation;large eddy simulation;expansion ratio

V435

A

1006-2793(2012)04-0474-05

2012-03-14;

2012-04-20。

王建儒(1978—)，男，博士，從事固體運(yùn)載火箭發(fā)動機(jī)總體設(shè)計。E-mail:wjr104zah@sina.com

book=35,ebook=287

(編輯:崔賢彬)