固體火箭發(fā)動機大渦模擬技術(shù)研究范文

本站小編為你精心準(zhǔn)備了固體火箭發(fā)動機大渦模擬技術(shù)研究參考范文，愿這些范文能點燃您思維的火花，激發(fā)您的寫作靈感。歡迎深入閱讀并收藏。

《現(xiàn)代防御技術(shù)雜志》2016年第3期

摘要:

研究固體火箭發(fā)動機工作中發(fā)生的聲渦耦合導(dǎo)致的不穩(wěn)定燃燒，以及燃燒室空腔變化對壓力振蕩的影響。對VKI縮比發(fā)動機流場采用大渦模擬(LES)方法進行數(shù)值模擬，確定數(shù)值方法可行。采用有限元和大渦模擬結(jié)合的方法，對某固體火箭發(fā)動機聲渦耦合進行仿真分析，對聲場分別采用理論方法和有限元數(shù)值方法對聲模態(tài)和聲學(xué)特性計算，確定有限元方法精度可靠及大渦模擬方法對流場旋渦流動不穩(wěn)定捕捉準(zhǔn)確，計算結(jié)果和試驗中壓力振蕩結(jié)果相符。表明研究的該固體火箭發(fā)動機發(fā)生了聲渦耦合引起的縱向1階聲不穩(wěn)定，且獲得了自由容積對壓力振蕩頻率及幅值的影響規(guī)律。

關(guān)鍵詞:

聲渦耦合;不穩(wěn)定燃燒;壓力振蕩;數(shù)值模擬;聲模態(tài);有限元;大渦模擬

0引言

發(fā)動機燃燒工作過程中，由于某種隨機擾動影響燃燒速率，燃燒室內(nèi)流場和燃燒過程、空腔的反射、聲場等相互作用，引起的不規(guī)則的燃燒室壓力變化，稱為不穩(wěn)定燃燒，也常稱為振蕩燃燒［1］。壓力振蕩會導(dǎo)致較嚴(yán)重的推力振蕩，輕則發(fā)動機振動，中斷燃燒，嚴(yán)重的情況下，和飛行器耦合共振，導(dǎo)致系統(tǒng)可靠性降低、失效乃至災(zāi)難性的爆炸事故發(fā)生［2－3］。近年來，隨著固體運載要求的提出以及導(dǎo)彈技術(shù)發(fā)展，對固體火箭發(fā)動機提出了高裝填、大推力比、初始大推力、大長細比等要求。國內(nèi)的發(fā)動機設(shè)計者為了滿足高性能指標(biāo)，提高設(shè)計標(biāo)準(zhǔn)，導(dǎo)致戰(zhàn)術(shù)地空導(dǎo)彈用固體火箭發(fā)動機屢次出現(xiàn)燃燒不穩(wěn)定現(xiàn)象。除此之外，國外研制的大型助推器，如美國的SpaceShuttle和歐洲的Ariane5［4－5］采用的大型分段式固體發(fā)動機相繼出現(xiàn)了較為嚴(yán)重的不穩(wěn)定燃燒，不穩(wěn)定燃燒研究也越來越重要［6］。渦脫落產(chǎn)生的聲學(xué)振蕩被認為是不穩(wěn)定燃燒問題中一個關(guān)鍵問題，最初Flandro和Jacobs［7］提出由于渦脫落激發(fā)的聲模態(tài)是引起不穩(wěn)定燃燒的一個重要因素，早期證據(jù)表明復(fù)雜幾何裝藥的戰(zhàn)術(shù)導(dǎo)彈上面級發(fā)動機的不穩(wěn)定燃燒由旋渦脫落驅(qū)動。后來有研究表明聲渦耦合導(dǎo)致的聲加強剪切層不穩(wěn)定的現(xiàn)象會放大了壓力震蕩，當(dāng)渦脫落頻率與聲振頻率相同時會形成了不穩(wěn)定燃燒［8］。

1問題描述

本文同時采用理論計算和有限元結(jié)合的方法，針對在工作末期發(fā)生了不穩(wěn)定燃燒的某固體火箭發(fā)動機不同燃燒時刻燃燒室聲腔進行分析和數(shù)值計算。并對該發(fā)動機翼柱燒盡、振蕩前、振蕩時、振蕩后等關(guān)鍵工作時刻建模，對不同時刻燃燒室流場采用大渦模擬方法數(shù)值模擬，初步分析認為由于流動不穩(wěn)定性產(chǎn)生了渦脫落，導(dǎo)致聲渦耦合，從而產(chǎn)生不穩(wěn)定燃燒，研究該不穩(wěn)定燃燒發(fā)生時刻是否發(fā)生聲渦耦合，對耦合情況下振蕩特性分析，并研究工作過程中自由容積變化對振蕩特性的影響。 

2數(shù)值方法與算例校驗

2．1聲學(xué)共振頻率理論方法燃燒室固有聲學(xué)頻率計算公式［9］:f縱=nc2L，(1)f切=0．293ncR，(2)f徑=0．61ncR，(3)式中:n為模態(tài)階數(shù);c為當(dāng)?shù)芈曀佟?/p>

2．2聲學(xué)共振頻率有限元方法采用基于有限元的方法，直接離散三維波動方程，可得到復(fù)雜裝藥聲腔的聲學(xué)特性。使用迦遼金法對簡化的亥姆霍茲方程進行離散，聲腔模態(tài)方程的單元矩陣［10］形式為Kf－ω2a(M)fp=0，(4)式中:Kf為聲剛度矩陣;Mf為聲質(zhì)量矩陣;p為聲特征向量;ωa為特征根。模型表面定義0位約束，使用有限元法求得圓周頻率ωa進而求得聲振頻率fa。

2．3大渦模擬方法通過在傅里葉空間或構(gòu)型空間將隨時間變化的N-S方程進行濾波可得到控制方程。考慮到氣體的可壓縮性，利用Favre平均對控制方程按式(5)簡化。f～=(ρf)ρ．(5)本文不考慮化學(xué)反應(yīng)，僅計算單組分工質(zhì)，濾波后連續(xù)方程、動量方程與能量方程分別為ρt+ρu～ixi=0，(6)t(ρu～i)+xi(ρu～iu～j)=－pxi+xj(σ～ij－τsgsij)，(7)t(ρe～)+xi(ρu～ie～)=xi(－pu～i－qi+u～iσ～ij－Hsgsi－Θsgsi)，(8)式中:“－”表示Reynolds平均;“～”表示Favre平均。濾波后分子粘性應(yīng)力σij與熱通量q～i分別為σ～ij=μ(u～ixj+u～jxi)－23μu～kxkδij，(9)qi=－KT～xi，式中:μ為動力粘性系數(shù);導(dǎo)熱系數(shù)K=CpμPr。將氣體工質(zhì)作理想氣體處理。采用WALE［11］亞格子模型對亞格子應(yīng)力張量τsgsij進行封閉。該模型對壁面區(qū)域進行了修正，改進了傳統(tǒng)的Smagorinksy亞格子模型耗散過大，不能用于湍流轉(zhuǎn)捩預(yù)測的不足。亞格子熱通量張量Hsgsij以及亞格子尺度粘性力變形功Θsgsij分別為Hsgsi=－μtprHixi，(11)Θsgs=Cερ(ksgs)3/2/V1/3，(12)式中:常數(shù)Cε為1．1。

2．4大渦模擬算例校驗為研究障礙物導(dǎo)致渦脫落，VKI實驗室設(shè)計了含有絕熱環(huán)和潛入式噴管空腔縮比冷流實驗發(fā)動機，Anthoine［12］還對該發(fā)動機進行試驗研究。流場計算區(qū)域尺寸結(jié)構(gòu)如圖1所示。VKI發(fā)動機尺寸已知，且經(jīng)過冷流實驗和數(shù)值模擬驗證，可以對本文采用的大渦模擬方法進行校驗，為此選用此發(fā)動機進行方法校驗。旋渦具有三維特性，產(chǎn)生、拉伸過程均為三維空間過程。二維模型計算對旋渦的捕捉和模擬總有欠缺，但三維模擬計算量太大，Mason［13］和張嶠［14］等人采用二維和三維大渦模擬對比后表明二維模擬可以成功的對ETM－03，RSRM，VKI等發(fā)動機進行數(shù)值模擬，且對振蕩頻率等特性很好的預(yù)估，為此本文綜合計算精度和效率采用二維模型。采用結(jié)構(gòu)化網(wǎng)格，結(jié)構(gòu)化網(wǎng)格質(zhì)量好，易于計算收斂;方便控制局部區(qū)域網(wǎng)格稀疏，易于對流動復(fù)雜區(qū)域和關(guān)注區(qū)域網(wǎng)格加密;減少其他區(qū)域網(wǎng)格，從而減少計算時間。此處對障礙區(qū)域、產(chǎn)生旋渦流動區(qū)域、旋渦和噴管頭部碰撞區(qū)域進行網(wǎng)格加密，局部加密網(wǎng)格如圖2所示。在發(fā)動機流場內(nèi)選取關(guān)鍵點監(jiān)測壓力波動，入口為質(zhì)量流率入口，溫度為285K，質(zhì)量流率為66.3kg/m2的理想氣體，采用軸向進氣［15］，比熱比為1.4，Pr數(shù)取0.71，出口為壓力出口，壓力為101325Pa，溫度為285K。流場初始計算先用穩(wěn)態(tài)計算得到一個初始流場，待流場穩(wěn)定后，改湍流模型為大渦模擬，為了避免連續(xù)方程和動量方程中心差分格式產(chǎn)生數(shù)值振蕩，采用BCD格式進行離散，并同時改為時間步長為4e－7的非穩(wěn)態(tài)計算，待流場穩(wěn)定開始監(jiān)測頭部壓力。圖3為部分時間段內(nèi)發(fā)動機頭部監(jiān)測點壓力曲線，由圖可以看出，發(fā)動機內(nèi)流場存在明顯的壓力振蕩，壓力振蕩區(qū)間為(170000～210000Pa)，對壓力曲線進行快速傅里葉變換后頻率分布圖如圖4，在0～2000Hz頻率內(nèi)，其中振蕩主頻為424Hz，壓力振蕩主要振蕩頻率為424，848，1272，1697Hz，呈明顯倍頻，主頻壓力振幅是其他振蕩頻率下振幅的7倍。Anthoine等人對此發(fā)動機分別進行過數(shù)值模擬和試驗，本文得到的振蕩頻率與其試驗和數(shù)值模擬結(jié)果對比如圖5和表1所示，四階主頻與Anthoine實驗和數(shù)值結(jié)果吻合，其中一階模態(tài)誤差最大，最大誤差僅為3.4%，表明本文所采用的大渦模擬方法對研究旋渦脫落導(dǎo)致流場不穩(wěn)定可行。

3數(shù)值方法與算例校驗

3．1計算模型本文研究的發(fā)動機為大長徑比固體火箭發(fā)動機，裝藥為前內(nèi)孔后翼，為此選取以下4個關(guān)鍵工作時刻進行研究，發(fā)動機在6s時后翼燒盡，燃燒室空腔變化趨于均勻，在15.5s左右發(fā)生振蕩燃燒，后翼燒盡后很長一段時間均會產(chǎn)生壓力振蕩，壓力振蕩頻率為173Hz。在6～15.5s工作時間段內(nèi)燃燒室內(nèi)自由容積隨著燃面的退移產(chǎn)生顯著的變化，會導(dǎo)致燃通比等參數(shù)明顯變化。為此，本節(jié)主要以后翼燃燒結(jié)束時刻6s，振蕩發(fā)生前的時刻14s，壓力振蕩發(fā)生開始時刻15.5s，振蕩后期時刻16s，4個時刻的發(fā)動機建模，計算模型如圖6。壓力監(jiān)測點在流場區(qū)域面積最小的16s時刻模型中選取，監(jiān)測點位置分別為:燃燒室頭部，突擴前，突擴后，突擴段中間，突擴段尾部，如圖7所示。

3．2邊界條件與計算參數(shù)燃面在不同時刻略有變化，燃燒室工作壓強在不同燃燒時刻也不同，這些因素在本文研究中不是主要因素，為此忽略燃面、壓強的細微不同，所有時刻燃面入口統(tǒng)一設(shè)置為質(zhì)量流率入口，質(zhì)量流量均為15．5s時刻的預(yù)估質(zhì)量流率，出口為壓力出口101325Pa，出口溫度300K，燃氣參數(shù)見表2。噴管出口為超聲速，出口界面參數(shù)由外推得到，且由于氣體在噴管加速至超聲速，下游對上游無影響，壁面采用無滑移壁面。

3．3燃燒室聲振頻率理論計算應(yīng)用式(1)計算燃燒室聲振頻率如表3所示。

3．4燃燒室聲學(xué)模態(tài)有限元仿真計算理論計算聲腔聲模態(tài)和振型時要求聲腔幾何構(gòu)型簡單，且邊界條件明確，鑒于許多燃燒室初始聲腔極復(fù)雜，本文中燃燒室相對簡單，但依然存在階梯凸臺等特征，為此對燃燒室三維聲腔建立有限元模型，噴管部分超聲速，相當(dāng)于燃燒室兩端封閉，只選取燃燒室主體部分，流體單元選取fluid30，模型表明定義零位移約束，取聲介質(zhì)密度為4.0kg/m3，平均聲速為1061m/s。圖8中為某固體火箭發(fā)動機發(fā)生不穩(wěn)定燃燒時前4階縱向振型聲壓分布云圖，前4階聲振頻率如表4所示，一階縱向頻率165Hz與試驗發(fā)動機壓力振蕩頻率相近，聲壓分布中，發(fā)動機頭尾部為聲壓波腹，表明發(fā)生了與一階縱向聲頻耦合的聲渦耦合。圖9為有限元方法和理論計算方法下聲振頻率對比，對比結(jié)果表明隨著模態(tài)階數(shù)增大，兩者計算誤差也有所增大，最大誤差約為5%。

3．5燃燒室各時刻壓力振蕩分析在發(fā)動機工作6s后翼柱燒盡，此后發(fā)動機內(nèi)腔變化趨于均勻，不同時刻計算入口條件一致的情況下，待監(jiān)測壓力波動穩(wěn)定后，觀測得知從6s到最后16s過程中均會產(chǎn)生一定的壓力振蕩。6s時壓力在波動中下降，14，16s時刻壓力振蕩在波動中上升，分別對以下4個時刻壓力波動周期平均計算后發(fā)現(xiàn)壓力波動周期處于5.86～6.26ms之間，壓力波動頻率均處于159～161Hz之間。若對壓力波動濾去波動量后，14，15.5，16s時壓力波動類似正弦波動，15.5s時的壓力波動頻率如圖10，周期性明顯。盡管都存在壓力振蕩，可以看出14，15.5s時刻壓力振蕩幅值較大，15.5s壓力振蕩相比其他時刻明顯高一個數(shù)量級，6s時由于后階凸臺明顯，會產(chǎn)生明顯的轉(zhuǎn)角渦脫落，可能會導(dǎo)致頭部壓力波動較大，14s時的壓力波動幅值是前2個時刻的10倍左右。15.5s壓力振蕩幅值遠大于其他時刻，是前兩時刻壓力波動20倍左右。模擬結(jié)果得到的15.5s的壓力振蕩最大，和發(fā)動機實際工作情況相符。監(jiān)測點head處的平均壓力隨發(fā)動機工作時間變化如圖11所示，隨著發(fā)動機工作，自由容積的變大，平均壓力先逐漸降低，在振蕩發(fā)生時刻附近上升，15.5s時最大，然后下降，平均壓力的最大波動范圍在5%左右。對壓力時間曲線，進行開窗快速傅里葉變換后，得到頻率、振幅曲線，均出現(xiàn)壓力振蕩，如圖12所示，15.5s時振幅最大，波峰對應(yīng)的前2階頻率分別為167，332Hz，其167Hz與一階縱向聲頻161Hz和165Hz接近。隨著時間變化，燃燒室自由容積增大，1階壓力振蕩頻率幾乎不變，其他頻率峰值逐漸平移，待發(fā)生壓力和聲振耦合共振后，隨著時間推移，產(chǎn)生共振的壓力振蕩頻率消失，到燃燒室工作末期，壓力振蕩幾乎消失。在燃面藥形設(shè)計中避開14～15.5s時刻的變化壓力振蕩明顯時的空腔，即可只發(fā)生低頻、小振幅微弱振蕩，而避免發(fā)生共振。

4結(jié)束語

大渦模擬和有限元方法結(jié)合的方法可以對某火箭發(fā)動機數(shù)值仿真與分析，表明該發(fā)動機發(fā)生了縱向1階頻率主導(dǎo)的聲渦耦合聲不穩(wěn)定燃燒，壓力振蕩的時刻特性與發(fā)動機試車工作情況相符。隨著發(fā)動機工作，自由容積增大的過程中，平均壓力先增大最然后下降，在振蕩發(fā)生時達到峰值，平均壓力的最大波動范圍在5%左右。表明本文方法可以對聲渦耦合引起的不穩(wěn)定燃燒情況進行預(yù)計，可以為發(fā)動機設(shè)計和抑振提供參考校驗。

參考文獻:

［1］孫維申．固體火箭發(fā)動機不穩(wěn)定燃燒［M］．北京:北京工業(yè)學(xué)院出版社，1987．

［2］劉佩進，金秉寧，李強．戰(zhàn)術(shù)導(dǎo)彈固體發(fā)動機燃燒不穩(wěn)定研究概述［J］．固體火箭技術(shù)，2012，35(4):446－456．

［3］王寧飛，張嶠，李軍偉，等．固體火箭發(fā)動機不穩(wěn)定燃燒研究進展［J］．航空動力學(xué)報，2011，26(6):1406－1412．

［5］胡大寧，何國強，劉佩進，等．翼柱型藥柱固體火箭發(fā)動機不穩(wěn)定燃燒研究［J］．固體火箭技術(shù)，2010，33(5):502－506．

［8］熊永亮，楊丹，郜冶，等．固體火箭發(fā)動機燃燒室內(nèi)聲波振蕩數(shù)值研究［J］．推進技術(shù)，2006，27(6):484－487．

［9］唐金蘭．固體火箭發(fā)動機原理［M］．北京:國防科技出版社，2013．

［10］楊向明，劉佩進，陳曉龍．翼柱型裝藥固體火箭發(fā)動機燃燒室聲場分析［J］．宇航學(xué)報，2008，29(5):1593－1597．

［14］張嶠，李軍偉，王偉臣，等．固體火箭發(fā)動機渦聲耦合特性數(shù)值研究［J］．推進技術(shù)，2011，32(3):349－354．

作者：韓磊 卞云龍 張衛(wèi)平 孫再庸 單位： 內(nèi)蒙動力機械研究所 中國航天科工集團第二研究院