固體火箭發動機約束方式對內流場與結構共振影響數值研究①

王 革,周博成,楊銘義,蒲煒強,關 奔,楊海威,楊澤南

(哈爾濱工程大學,哈爾濱 150001)

0 引言

固體火箭發動機廣泛用于火箭、導彈等飛行器,燃燒不穩定一直困擾著研究人員[1-3]。一些發動機即使通過了地面試車試驗,在實際飛行過程中仍可能出現壓力振蕩和結構異常振動等燃燒不穩定問題[4-5],可總結為固體火箭發動機的天地不一致問題。為了在更低成本的地面試車中模擬發動機的飛行工作狀態,誘發飛行過程中出現的燃燒不穩定,國內外學者在地面靜止臺架上對發動機施加脈沖[6],或利用火箭撬[7]、過載臺[8]和振動臺等設備施加載荷,部分再現了燃燒不穩定現象。上述方法中發動機受到的約束狀態與飛行狀態差別較大[9],發動機與試車臺固定連接,提高了發動機的整體剛度和質量,從而改變發動機的固有頻率,可能抑制發動機因流場與結構頻率一致而發生的共振。

在眾多誘發燃燒不穩定的因素中,除了渦脫落主導和燃燒響應主導的兩種典型因素以外[9],發動機結構與燃燒室聲腔的耦合,也可能誘發燃燒不穩定[2,10]。DOTSON等[11]研究證實了大型固體運載火箭發動機結構和流動之間的共振會增加壓力振蕩的振幅,激振力和系統響應之間存在近似恒定的相位差,體現為飛行試驗中壓力振蕩振幅高于相同燃燒時間的地面試驗中測得的振幅。MASON等[12]通過數值模擬評估了由渦流脫落引起的發動機內部壓力振蕩、發動機內部聲學模態和發動機結構振動模態之間的復雜相互作用,研究了在發動機燃燒后期的結構模態與一階聲學模態重合導致的推力振蕩增加。ZHANG等[13]將穩定飛行狀態下的固體火箭簡化為雙層時變截面軸向運動自由梁,從非線性自由振動的角度研究了發動機在不同工況下的非線性固有頻率。加拿大學者GREATRIX[14-19]用數值和實驗方法在發動機殼體加速度對燃燒不穩定的影響問題上進行了長期研究,建立了描述結構振動加速度對推進劑燃燒響應影響的準一維歐拉方程模型[14],并在之后的研究中,陸續探討了兩相流[15]、藥柱形狀[16-17]、尺度效應[18]、加速度方向角[19]等因素對燃燒不穩定的影響,形成了一套相對完備的理論體系。

對于固體火箭發動機天地不一致問題,如今亟需構建一種合適的方法在地面試車條件下對發動機建立更加接近真實自由飛行環境的載荷[9],以用于在更低成本的地面試驗條件下評估發動機飛行狀態的工作穩定性。本文使用商業軟件ANSYS Workbench對某型固體火箭發動機進行模態分析與雙向流固耦合計算,通過設計模擬自由飛行和地面試車外載荷的兩種約束條件,分析因結構約束方式不同導致的流場與結構共振的現象及規律,為固體火箭發動機天地不一致問題的研究提供參考。

1 研究方法

采用模態分析方法求解聲固有頻率與結構固有頻率,預估發動機結構模態與聲模態的共振頻率。通過ANSYS Workbench商業軟件進行模態分析,分別使用Modal Acoustics和Modal模塊計算聲模態與結構模態。采用雙向流固耦合方法求解固體火箭發動機殼體、絕熱層、推進劑等結構和燃燒室內流場耦合作用下的發動機工作不穩定問題。流體和固體區域分別使用Fluent和Transient Structural模塊求解,數據通過System Coupling模塊實現雙向傳遞。

1.1 物理模型

某型固體火箭發動機在工作后期的某時刻出現了壓強異常躍升和極限環振蕩的燃燒不穩定問題,該時刻發動機的結構剖面如圖1所示,此時推進劑的厚度較薄。

圖1 發動機結構示意圖

發動機結構由碳纖維殼體、絕熱層、推進劑和噴管組成,具體材料參數如表1所示。碳纖維殼體根據纖維纏繞角度采用各向異性材料,絕熱層和推進劑假設為線彈性材料,噴管材料較為復雜,這里將其簡化假設為鈦合金材料。在殼體和推進劑外表面總長度的10%、50%、90%處設置3個三向位移監測點xshn和xprn,其中n=1、2、3。在內流場軸線總長度的10%、50%、90%處設置內彈道壓力監測點pn,其中n=1、2、3。

表1 發動機結構材料參數

假設流場中燃氣為化學平衡流的理想氣體,氣體在流動過程中不發生化學反應,燃氣總溫T=3532 K,工作壓力約10 MPa,燃氣分子摩爾質量Mgas=30.3 kg/kmol,燃氣比熱容比γ=1.16[20]。

1.2 數值方法

(1)流場控制方程

在固體火箭發動機工作過程中,高雷諾數可壓縮流滿足RANS(Reynolds Average Navier-Stokes)方程,其質量、動量和能量守恒方程為

(1)

(2)

(3)

利用理想氣體狀態方程,使方程組閉合。

p=ρRT

(4)

式中ρ、u、p分別代表密度、速度矢量、壓力;I是單位張量;En=e+1/2u2為總能量;e為內能(e=h-p/ρ);h為生成焓;T為溫度;氣體性質參數μ和Г分別為粘性系數和熱導率;μt和Prt為湍流渦粘系數和普朗特數;R為氣體常數,R=408.3 J/(kg·K)。

湍流粘性系數μt由湍流動能k和比耗散率ω計算得到,利用SSTk-ω湍流模型[20]對其進行求解。

(2)結構控制方程

發生燃燒不穩定時,發動機內流場的壓力分布隨時間變化,其作為載荷作用于發動機結構時會產生結構響應,而這種響應主要受結構慣性力和阻尼作用的影響。發動機的殼體、絕熱層、藥柱和噴管等結構可視為一個線性的結構動力學系統,通過虛功原理建立的有限元半離散運動方程表示如下:

(5)

式(5)用Newmark方法[22]進行時間積分求解,其中速度和加速度通過以下方程確定:

(6)

(7)

式中 矢量的角標n和n+1表示該矢量在時間步tn和在下一時間步tn+1;δ和α為Newmark積分參數;Δt為時間步長。

(3)雙向流固耦合方程

雙向流固耦合的求解遵循耦合交界面處流體與固體變量相等或守恒,耦合的交界面需要滿足如下耦合邊界條件。

動力耦合邊界條件:

(8)

運動耦合邊界條件:

(9)

式中Q和F為界面上的力;上標f代表流體與結構交界面處的流體界面,s代表交界面處的固體界面,fs代表流固耦合交界面;下標c代表耦合。

1.3 邊界條件及計算流程

(1)內流場邊界條件

使用用戶自定義函數(UDF)在Fluent中定義推進劑與流場接觸的表面為質量輸入面,通過定義質量和能量源項方式加入12.56 kg/(m2·s)的質量通量,以使燃燒室流場壓力穩定在10 MPa。同樣使用UDF定義前封頭絕熱層內表面為脈沖輸入面,假設在燃燒不穩定發生初期頭部產生了一個幅值為1 MPa的壓力脈沖并向下游傳遞,質量通量為788.31 kg/(m2·s)。流場的邊界條件如圖2所示,絕熱層和噴管壁面為無滑移壁面,推進劑表面為質量入口,噴管出口為壓力出口,總壓0.01 MPa,總溫300 K。

圖2 流場邊界條件

(2)結構約束方式

圖3為三種發動機約束方式示意圖,圖中黃色高亮區域為約束施加位置。

圖3 發動機三種約束方式示意圖

圖3(a)通過禁止前封頭環形區域位移和禁止尾部環形區域的徑向位移,模擬地面試車條件下頭部和尾部受試車臺架約束的實際狀況。因其相對于其他兩種約束方式的強度最強,以下簡稱“強約束”狀態。發動機前封頭通過襯裙等結構與飛行器的前艙連接,圖3(b)的約束禁止前封頭環形區域產生軸向位移。當火箭在均勻來流的均速穩定飛行過程中,發動機前封頭環形區域受力與發動機推力大小一致,方向相反,且結構的徑向位移不受限制。該約束方式近似模擬了火箭在均勻來流的穩定飛行過程中,前封頭通過襯裙等結構與飛行器前艙連接的受力狀態。因其約束強度較弱,命名為“弱約束”狀態。圖3(c)為自由狀態,但因模型中不含有彈頭等額外部件,該自由狀態與飛行器實際飛行狀態并不相符。且在雙向流固耦合計算中由于其沒有結構的約束限制,會導致計算不收斂,因此僅作為模態分析的一種對比算例,命名為“無約束”狀態。如果將發動機視為一根等截面梁,可將上述約束條件簡化為圖中右側部分的分析簡圖。約束條件施加在等截面梁的兩端,三種狀態分別代表固支-簡支,簡支-自由和自由-自由條件。

(3)流固耦合計算流程

首先對流場進行穩態計算,收斂結果作為初始流場條件進行流固耦合瞬態計算,計算總時長0.05 s,時間步長2.5×10-5s,在0時刻對流場施加1×10-3s的壓力脈沖后觀察自由振動結果。流體域和固體域的交界面設置為Data Transfer耦合面,每個時間步內信息交換5次保證收斂。

1.4 網格繪制及無關性驗證

流體區域使用ICEM繪制六面體結構網格,使用動網格技術實現流體區域因結構位移造成的形變,每個迭代步進行10次網格重構,用以提升網格質量。

對流體域網格進行無關性驗證,共設置4個不同的網格量:800 000、1 600 000、3 200 000和6 400 000。在頭部脈沖面施加一個1 MPa的壓力躍升,監測點p2的壓力振蕩結果如圖4所示。

圖4 流體域網格無關性驗證

800 000網格的預測峰值壓力明顯低于其他網格。當網格數大于3 200 000時,結果的一致性較好,因此使用該網格量進行后續計算。

對于結構域使用ANSYS的Mesh模塊進行網格繪制。由于物理模型的殼體和絕熱層厚度較薄,在給定整體網格尺寸的同時進行局部加密,保證薄壁結構在厚度方向上至少有4層網格。使用ANSYS的Static Structural模塊對結構網格進行靜載荷分析,以驗證網格無關性,共設置網格量為50 000、100 000、200 000和400 000的四種網格。在發動機內壁面施加10 MPa的壓力,記錄發動機殼體監測點xsh2的位移。400 000網格的結構最大位移為3.512 4 mm,其余網格分別與其相差4.72%、0.24%、0.08%,因此選擇200 000網格用于結構計算。

2 結果與討論

2.1 模態分析

模態分析常用于預測發動機聲腔與結構的固有頻率。圖5為發動機前四階無量綱聲壓分布,由于噴管喉口下游流動速度為超音速,下游信息無法向前傳遞,因此可將計算域從噴管喉口處截去。

圖5 發動機前四階聲模態聲壓分布

聲腔前四階固有頻率為157、322、491、658 Hz,因為發動機內流場為細長體結構,聲模態主要為軸向聲壓分布。對于一階聲振型,燃燒室兩端為聲壓相位相反的聲壓波腹,中間為波節;二階兩端和中間均為波腹;更高階數則出現了更多的波節和波腹。

結構的固有頻率主要受剛度分布、質量分布和結構阻尼的影響,在同一物理模型條件下,約束方式通過改變結構的剛度分布影響固有頻率。使用預應力模態分析方法,首先通過施加內部靜載荷獲得發動機在實際工作壓力下的應力應變場,再根據該結果計算得到模型的固有頻率(如表2所示),主要振型如圖6所示。圖6中,陰影部分為模型初始外形。為了清晰地顯示結構的振型特征,結構形變被適當放大。發動機的振動形式主要分為彎曲、拉伸和呼吸三種,其中彎曲振型又可根據振型的波峰數量分為一階、二階和三階。軸向拉伸振型體現為發動機在沿對稱軸方向的往復拉伸-收縮形變,呼吸振型體現為發動機厚度較薄的中后部區域的往復徑向壓扁-回彈過程。

表2 發動機結構固有頻率

圖6 發動機前五階固有振型

傳統的聲不穩定是推進劑燃燒過程與發動機燃燒室聲學過程相互作用的結果,當壓強振蕩的頻率與燃燒室聲腔的固有頻率相近時發動機是一個自激聲振系統[5]。當考慮結構振動時,需兼顧結構振動與聲場的關系,對比發動機聲模態頻率與結構固有頻率如圖7所示。圖7中,橫線為發動機固有聲頻,數據點的縱坐標為不同約束方式的結構固有頻率。類比自激聲振系統的定義,當發動機聲模態的頻率與結構固有振動頻率相同或相近時,發動機聲場與結構場在一定的內外部激勵條件誘發下可能發生共振。由圖7可知,不同約束方式對一階彎曲頻率和一階呼吸頻率的影響較小。在二階彎曲和三階彎曲形式中,弱約束和強約束的頻率較自由狀態大,且約束強度越高頻率越高。而對于軸向拉伸頻率,約束使該狀態的固有頻率減小。一階聲模態的固有頻率與一階彎曲頻率接近,二階聲模態與弱約束和強約束的軸向拉伸頻率接近,在這兩個頻率下可能發生共振,而其他模態未體現出明顯的共振趨勢。模態分析能夠為聲場與結構的共振提供預測,但是否發生共振還需要進行流固耦合計算或試車試驗。

圖7 不同約束方式固有聲頻與結構固有頻率對比

2.2 不同脈沖條件流場壓力振蕩對比

為了探究不同脈沖條件的燃氣質量注入對內彈道壓力振蕩的影響,規避因脈沖條件不同導致的流場振蕩頻率和幅值差異,共設計了5種脈沖工況進行流場計算,具體脈沖強度和作用時間如圖8標注所示(注:脈沖時間-脈沖強度)。流場計算結果表明,在壓力穩定的流場中施加一定時間和一定強度的質量脈沖會引起壓力振蕩,脈沖強度越大振蕩幅值越大,脈沖時間越長振蕩幅值也越大。

圖8 不同脈沖條件壓力振蕩時域曲線

對不同脈沖工況數據結果進行短時傅里葉變換(STFT),得其頻譜如圖9所示。在同一流場中,內彈道壓力振蕩的各階頻率以基頻的倍數增長,基頻與發動機一階聲振頻率相同,為157 Hz。不同脈沖條件導致的內彈道振蕩頻率基本一致,但幅值有所不同。質量脈沖以對流場輸入質量和能量形式引起壓力振蕩,脈沖的能量與脈沖幅值和作用時間的乘積成正相關,而能量又與振蕩幅值成線性關系。相同能量所激發的一階脈沖幅值基本一致,如1 ms-0.5 MPa和0.5 ms-1 MPa的一階振幅接近。這說明固體火箭發動機工作時,注入額外的燃氣質量會引起壓力振蕩,擾動所激發的內彈道壓力振蕩頻率基本一致,幅值有所不同。

圖9 不同脈沖條件壓力頻譜

2.3 流固耦合分析

(1)弱約束結構振動分析

對弱約束發動機進行1 MPa脈沖激勵條件下的流固耦合計算,獲得結構振動時域曲線并進行頻譜分析。計算結果中,x方向的振幅較低,而徑向y、z方向的振動曲線基本一致。因此,僅展示弱約束發動機各監測點y方向振動的時域曲線,如圖10所示。在脈沖觸發后,結構首先出現高頻高幅值振動,這是由于脈沖觸發導致的結構受迫振動。隨后,其衰減至低頻低幅值振動,該階段主要為發動機內彈道振蕩與結構的耦合振動。對比不同位置的監測結果,發動機殼體尾端監測點xsh3的振動幅值最大,且在低頻段出現更加規律的波動。推進劑的尾部監測點xpr3位于xsh2和xsh3之間,因此主要體現為發動機中后部振動特性,與xsh2接近。中部殼體xsh2和推進劑xpr2的振幅也較大,頭部位移最小。其中,殼體頭部xsh1的振動幅值因量級過小在圖10中近似于直線。

圖10 弱約束發動機振動時域曲線

為了分析發動機的振動頻率,將弱約束發動機振動時域曲線的低頻段進行短時傅里葉變換(見圖11),利用靜載最大位移(3.512 4 mm)對監測點振動數據進行無量綱處理。在低頻段主要出現了四階頻率,分別為151、353、554、1210 Hz。其中,1210 Hz與壓力脈沖有關,為結構脈沖觸發響應頻率。前三階頻率為發動機結構振動與壓力振蕩的耦合作用結果,發動機不同位置的監測點測得振動頻率基本一致,但幅值有所不同。發動機尾部監測點因距離約束位置最遠,其振幅最大。除了發動機殼體尾部監測點xsh3的主頻為554 Hz和推進劑尾部xpr3為151 Hz外,其他位置的主頻均為第二階的353 Hz,且除了距離約束位置最遠的殼體尾部和推進劑尾部監測點處存在一階151 Hz的振動頻率外,其他位置均不存在該階頻率。

圖11 弱約束發動機結構振動低頻段頻譜

(2)強約束結構振動分析

強約束發動機結構y方向振動時域曲線見圖12。振動初始為脈沖激勵引起的高頻強迫振動,隨后為低頻耦合振動。相比于弱約束方式,強約束發動機在其低頻段頭部與尾部監測點的振幅明顯較低,僅發動機中部的xsh2和xpr2有較明顯波動。殼體尾部監測點xsh3位于固定約束作用位置上,因此振動幅值為零。

圖12 強約束發動機振動時域曲線

對強約束振動時域曲線的低頻段進行短時傅里葉變換如圖13所示,結構振動的主頻為314 Hz,不同監測點位置的振動頻率基本一致。相對于弱約束條件,強約束頻譜中不存在151 Hz頻率的振動。由模態分析可知,該振動形式為一階彎曲,發動機強約束方式通過模擬地面試車狀態限制了其頭尾的徑向位移從而抑制了一階彎曲振型。強約束提升了結構剛度,改變了發動機的結構動力學特性。發動機振動的最大幅值出現在殼體中部的xsh2位置上,推進劑中部xpr2次之。發動機中部因遠離兩端約束位置,其結構振幅最大。圖13中,1178 Hz為脈沖觸發響應頻率,與弱約束工況中的1210 Hz基本一致。

圖13 強約束發動機結構振動低頻段頻譜

(3)流場壓力振蕩對比

對比弱約束和強約束發動機流固耦合計算的流場壓力振蕩結果與不考慮耦合的流場(以下簡稱純流場)計算結果見圖14。壓力初始為高頻低幅值振蕩,隨后出現規律的低頻高幅值振蕩。隨著時間的推移,振蕩幅值有所衰減。在初始幾個峰值中弱約束工況和強約束工況比純流場工況的壓力絕對值低。后續發展過程中,對比第4個完整波形,純流場的壓力幅值0.56 MPa低于弱約束和強約束結果的幅值0.59 MPa,幅值增大5.3%,且在其他壓力峰中,弱約束和強約束工況的壓力振幅也均比純流場工況的幅值高,說明結構振動對壓力振蕩的幅值有放大作用,結構通過振動對流場做功。兩種約束方式在前4個峰值未體現出較大差異,而在第5個峰值后均出現了雙峰現象,且波形開始出現不同。

圖14 不同約束方式壓力振蕩

對流場壓力數據進行短時傅里葉變換,如圖15所示。弱約束、強約束和純流場工況的主頻均為第一階167 Hz,幅值0.27 MPa。弱約束方式的內彈道壓力振蕩頻率與純流動頻率一致,為167、334和500 Hz,基本滿足倍頻關系。強約束使流場的第二階頻率略有升高,為350 Hz,且幅值有所降低。這說明結構的約束方式會對內彈道壓力振蕩頻率產生影響,在一定范圍內,結構的剛性越強,壓力振蕩的頻率越高。

圖15 不同約束方式壓力振蕩頻譜

(4)耦合結果分析

綜合模態分析和流固耦合計算結果,數值預測的發動機共振對比如圖16所示。圖16中,橫實線為流固耦合計算得到的結構振動頻率,橫點線代表模態分析得到的結構固有頻率,豎虛線為流固耦合計算得到的流場壓力振蕩頻率,豎點劃線為模態分析得到的聲頻率。橫縱坐標等比例繪制,圖中沿45°角設定一條寬為10 Hz的灰色區域,該區域的橫縱坐標相差10 Hz,此處稱為共振帶。當結構模態頻率與聲模態頻率相交在灰色區域時,模態分析預測在該頻率范圍可能存在共振現象。而當結構振動頻率和流場振蕩頻率也相交在該頻率范圍時,則判定流場與結構發生共振。紅色圓形區域代表共振區域,藍色圓形區域代表模態分析預測固體與聲場頻率接近,但流固耦合結果沒有出現共振。

弱約束工況結果如圖16(a)所示,在紅色區域內4條線相交于共振帶。模態分析結果中聲頻157 Hz和結構固有頻率166 Hz預測該頻率附近可能發生共振,流固耦合結果的流場振蕩頻率167 Hz和弱約束結構一階振動頻率151 Hz與模態預測頻率相近,驗證了共振特征。因此,在150～170 Hz頻率范圍內發動機可能發生流場與結構的共振。同時,圖16(a)藍色區域的結構模態與流場頻率和聲頻率相交,雖然模態分析預測該頻率范圍存在共振現象,但結構振動頻率高于預測頻率,沒有發生共振。強約束工況結果見圖16(b),由模態分析可知,在150～170 Hz附近可能出現共振,但該階數的結構頻率不存在,因此流場與結構在該頻率不發生共振。在310～340 Hz的藍色區域,模態分析預測出的共振頻率不與流場和結構頻率重合,也不會發生共振。因此,數值結果表明,模擬自由飛行條件的弱約束發動機存在150～170 Hz頻率范圍的流場與結構的共振現象,而模擬試車條件的強約束發動機不發生共振。

(a)Comparison of weak constraint resonance (b)Comparison of strong constraint resonance

3 結論

本文針對某型固體火箭發動機在工作后期出現燃燒不穩定的天地不一致問題,模擬自由飛行和地面試車狀態設計了弱約束和強約束兩種結構約束方式。針對兩種約束方式的發動機進行了模態分析與雙向流固耦合計算,探究了約束方式對發動機流場與結構共振的影響,得出以下結論與展望:

(1)不同的約束方式改變了發動機結構的固有頻率,可能導致結構振動頻率與流場振蕩頻率一致而發生共振。

(2)脈沖作用下發動機流場與結構的共振會使結構對流場做功,耦合結果的壓力振幅較純流動的振幅更大。

(3)模擬自由飛行條件的弱約束發動機存在頻率范圍在150～170 Hz的流場與結構共振。而模擬地面試車狀態的強約束發動機沒有出現共振。

(4)在固體發動機工作過程中,地面試車與真實飛行條件的約束方式不同,可能會導致在飛行過程中出現地面試車沒有測量到的流場與結構共振現象。因此,亟需發展一種能夠在地面試車條件下模擬飛行條件的工作穩定性預估模型和實驗手段,使因結構和流場共振的燃燒不穩定現象能夠在地面條件下提前暴露。本文提出的兩種約束方式需要地面試車和飛行數據驗證,以證明該方法的準確性。