基于CFD/CSD耦合的火箭跨聲速氣動阻尼特性分析

李泳德, 郭 力, 季 辰

(中國航天空氣動力技術研究院, 北京 100074)

引 言

通常情況下人們認為氣動力對火箭的振動起到阻尼作用, 即氣動阻尼為正值。然而隨著大推力火箭發(fā)展, 火箭的長細比逐漸加大, 導致彎曲剛度越來越小, 同時為了滿足有效載荷的外形要求, 火箭頭部整流罩尺寸不斷加大, 后續(xù)箱體的直徑卻保持不變, 形成了典型的錘頭體外形。國內外大量的火箭研制經驗表明[1-9], 對于此類錘頭體外形火箭的氣動設計, 必須要進行動態(tài)氣動載荷與動態(tài)氣彈穩(wěn)定性分析, 否則設計的疏忽可能會導致火箭結構出現(xiàn)毀滅性的破壞進而導致發(fā)射失敗。

目前常用的衡量氣彈穩(wěn)定性的方法是通過風洞試驗來獲取氣動阻尼系數(shù)。早在1963年, 美國國家航空航天局Ames研究中心(NASA Ames Research Center)采用半剛性模型開展試驗研究[10], 獲取火箭頭部的氣動阻尼來評估其穩(wěn)定性, 但這只能用來模擬火箭彎曲振型前節(jié)點之前部分的結構動力學特性。直到蘭利研究中心(NASA Langley Research Center)開發(fā)了全彈性模型氣動阻尼試驗技術, 其可以模擬整體的結構動力學特性以及氣動外形, 并應用于多款運載火箭研制[11-15]。國內, 中國航天空氣動力技術研究院對氣動阻尼問題開展過較多的研究[16-20], 從模型設計方法、 模型制作工藝、 試驗機構設計和數(shù)據(jù)處理等諸多方面, 逐步改進實現(xiàn)了從半剛性模型到全彈性模型的過渡, 并在多個型號上得到驗證。

然而通過風洞試驗研究氣動彈性問題, 技術難度大, 試驗成本高, 同時幾乎不可能開展全尺寸試驗。因此通過數(shù)值計算的方法開展相關研究是另一種重要的手段。劉子強等[21]實現(xiàn)了通過數(shù)值計算確定氣動阻尼系數(shù)的技術和方法, 并與試驗結果進行對比, 證實了該方法的可靠性。冉景洪等[22]通過模態(tài)數(shù)據(jù)結合準定常理論的方法分析了減阻桿加后體這一彈性結構的氣動阻尼, 結果表明減阻桿造成的分離流會對后體的氣動阻尼系數(shù)產生影響。朱劍等[23]針對新一代捆綁式運載火箭發(fā)展了非結構網(wǎng)格下的氣動阻尼計算方法, 并分析了攻角、 Mach數(shù)等參數(shù)對氣動阻尼的影響。

本文在之前的計算方法[23]的基礎上采用IDDES模型, 考慮脈動壓力的影響, 通過強迫振動的方式, 針對捆綁式運載火箭的某一特定模態(tài)進行數(shù)值計算仿真, 研究前節(jié)點位置, 振動振幅, 脈動壓力等參數(shù)對氣動阻尼的影響規(guī)律。

1 計算方法

圖1為本文所用的捆綁式運載火箭的計算模型, 是典型的錘頭體結構。在跨聲速階段, 其頭部會產生激波造成激波邊界層干擾, 而在錘頭體外形的過渡段會出現(xiàn)氣流分離。為探究各部分氣動阻尼的變化, 將整個箭體分為頭部、 過渡段、 彈身3個部分。

圖1 表面網(wǎng)格及區(qū)域劃分Fig. 1 Surface grid and region division

1.1 流場仿真模型

本文分別用Reynolds平均法(Reynolds-averaged Navier-Stokes, RANS)和改進的延遲分離渦模擬(improved delayed detached-eddy simulation, IDDES)[24-25]進行計算, 在RANS方程中, 將變量分為平均值和波動值兩部分, 對于速度分量有

(1)

(2)

其中,k和ω分別代表湍流動能和湍流耗散率,Γk和Γω分別代表k和ω的有效擴散系數(shù),Gk和Gω分別代表k和ω的生成率,Yk和Yω分別代表k和ω的耗散率。因此RANS方法只能計算大尺度的平均流動, 本文采用IDDES方法計算脈動壓力對氣動阻尼的影響。

IDDES方法是由分離渦模擬(detached-eddy simulation, DES)方法改進而來, 其本質思想與DES方法相同, 是想以網(wǎng)格尺度和模型中的特征尺度隱式劃分RANS和大渦模擬(large-eddy simulation, LES)區(qū)域, 使其既能處理RANS方法無法得到的脈動場, 也能降低LES方法在模擬高Reynolds數(shù)流動時所需的計算資源。區(qū)別在于當邊界層較厚或者分離區(qū)域較窄時, DES方法會出現(xiàn)如模型應力損耗(modeled stress depletion, MSD), 網(wǎng)格誘導分離(grid-induced separation, GIS)以及對數(shù)層不匹配(logarithmic-layer mismatch, LLM)問題[24], 而IDDES模型通過改良計算區(qū)域劃分, 結合延遲分離渦模擬(delayed detached-eddy simulation, DDES)和壁面模型大渦模擬(wall-modeled large-eddy simulation, WMLES), 定義新的長度尺度解決了這些問題, 具體公式詳見文獻[25]。

流場網(wǎng)格如圖2、 圖3所示, 邊界層采用棱柱層結構, 并調整第1層網(wǎng)格高度使得y+小于1, 遠場部分采用六面體結構網(wǎng)格, 與邊界層的過渡層采用非結構網(wǎng)格。整體網(wǎng)格單元數(shù)量為4.2×106。

圖2 y方向截面網(wǎng)格示意圖Fig. 2 Schematic diagram of cross-sectional grid in the y-direction

圖3 x方向截面網(wǎng)格示意圖Fig. 3 Schematic diagram of cross-sectional grid in the x-direction

物面邊界條件為無滑移壁面條件, 遠場采用壓力遠場邊界條件, 湍流模型采用SSTk-ω模型, 采用密度基求解, 氣體黏性采用Sutherland定律, 空間離散采用2階迎風格式, 對流通量采用Roe格式。

1.2 結構分析模型

結構與流場耦合分析過程中, 結構部分可以采用模態(tài)方法描述。結構模態(tài)可以通過有限元方法與結構模態(tài)試驗方法獲得。本文采用有限元分析結果獲得的模態(tài), 圖4所示為結構的前3階模態(tài), 本文只分析計算結果中氣動阻尼最小的第2階模態(tài)。

(a) f=1.200 Hz

由于火箭結構外形簡單, 一般不考慮其扭轉影響, 因此可以將其簡化為簡單的梁模型, 這樣就可以給出其模態(tài)振動方程

(3)

式中,qi為第i階模態(tài)的廣義位移,bi為第i階模態(tài)的結構阻尼系數(shù),ωi為第i階模態(tài)的固有頻率,fi為第i階模態(tài)下質量歸一化的廣義氣動力。若將fi按照Taylor展開并略去高階項, 可以將其轉化為氣動阻尼項與氣動剛度項的形式, 則式(3)可寫為

(4)

式中,Bi為氣動阻尼系數(shù),Ki為氣動剛度系數(shù), 研究表明[26], 氣動剛度相對于結構剛度為小量可以忽略不計, 而在計算中結構阻尼往往設置為0, 因此氣動阻尼可以直接反映其氣彈穩(wěn)定性。

1.3 氣動阻尼分析原理

氣動阻尼的分析可以采用強迫振動或者自由振動的方式進行, 這兩種方法獲得的時域數(shù)據(jù)不同, 提取氣動阻尼的方式也不同。強迫振動方法初始演化過程較短, 因此計算量較小, 同時能夠分析某一種振動形式的氣動阻尼, 明確該振動形式是收斂還是發(fā)散。分析過程中能夠獲得不同部位與部件的氣動阻尼。但是對于多模態(tài)相互作用引起的發(fā)散(例如顫振)較難預測。自由振動方法需要一定的自由演化時間才能夠對時域數(shù)據(jù)進行分析, 不過自由振動方法能夠獲得最能夠吸收能量的模態(tài)及其振動頻率。

對于本研究所關注的問題, 氣動載荷對結構振動的過程中氣動阻尼的影響較大, 而對氣動剛度與氣動質量影響較小, 即結構的固有振動頻率受到來流的影響較小, 其穩(wěn)定性問題主要由氣動阻尼的正、 負引起, 所以采用強迫振動方法分析。

強迫振動下結構做簡諧模態(tài)振動

qi(t)=Asin(ωit)

式中,A表示振動的振幅, 將其代入計算氣動力的公式中[21]并做正交積分可得

(5)

式中,Mi為第i階模態(tài)的模態(tài)質量,T為整數(shù)倍周期,G為廣義氣動力。根據(jù)式(5)便可以得到局部或分區(qū)域的氣動阻尼。

1.4 耦合計算流程

首先進行模態(tài)分析, 以確定結構的模態(tài)頻率與振型, 用以設計強迫振動的頻率和振幅。非定常流場計算前先進行定常流場計算, 來加快非定常計算的演化速度并增強收斂性, 結構節(jié)點位移通過徑向基函數(shù)(RBF)插值方法[27]映射到氣動網(wǎng)格節(jié)點上, 來進行網(wǎng)格的變形, 這里徑向基函數(shù)選用Wendland C2, 如下所示

最后將計算出來的廣義力提取出來, 截取演化完畢的整數(shù)倍周期, 進行氣動阻尼計算。耦合計算流程圖如圖5所示。

圖5 耦合計算流程圖Fig. 5 Flow chart of coupled calculation

2 結果分析與討論

2.1 流場分析結果

計算的來流Mach數(shù)范圍為0.7～1.2。其中中截面的壓力分布如圖6所示?？梢钥闯鲈陬^部出現(xiàn)了膨脹波以及跨聲速激波, 在過渡段存在流動分離, 隨著Mach數(shù)的增大, 頭部低壓區(qū)域逐漸擴張, 并且能明顯看到, 在流動再附的位置產生了再附激波。

(a) Ma=0.70

2.2 氣動阻尼分布

通過上述流場分析, 可以看出火箭不同部位流動結構并不相同, 在頭部與箭身上, 流動主要為附著流動, 而在過渡段會出現(xiàn)較為復雜的波系結構以及流動分離。針對不同的流動結構隨流向站位x的變化, 設該位置上廣義力與廣義位移的相位差為φ(x), 并且簡諧振動沒有引入其他模態(tài)的廣義力, 則廣義力的表達式為

G(x,t)=Fgen·sin[ωt+φ(x)]+F0

(6)

其中,Fgen為廣義力的振動幅度,F0為廣義力的常數(shù)偏移量。將式(6)代入到式(5)中得到

其中, 廣義力的常數(shù)偏移量F0的積分為0, 因此省略。通過將等式中的正弦函數(shù)部分進行和差化積得到

(7)

式(7)中第1部分在整個周期中的積分為0, 只有第2部分保留, 因此得到

(8)

式(8)中積分部分恒為正值, 決定整個氣動阻尼的部分只有相位角φ(x)的正弦值sin[φ(x)], 為了能夠更加直觀地獲得相位角與氣動阻尼B之間的關系, 須將符號轉化為對應的正弦函數(shù)轉角, 根據(jù)正弦關系, 此轉角為π, 因此得到

(9)

圖7為氣動阻尼變化曲線, 可以看出隨著Mach數(shù)的增大, 整體氣動阻尼先增大后減少, 在Mach數(shù)為0.98時達到最大值, 過渡段與箭體的氣動阻尼變化趨勢與整體基本相同, 而頭部區(qū)域則不同, 是隨著Mach數(shù)的增大一直增大, 只是增長速率變緩。

圖7 有助推時氣動阻尼變化曲線Fig. 7 Aerodynamic damping change curve with boost

根據(jù)式(9), 得到相位角與氣動阻尼B之間的關系為: 當φ(x)∈(-π, 0)時, 相位角滯后, 氣動阻尼B為負值; 當φ(x)∈(0, π), 相位角提前, 氣動阻尼B為正值; 為當φ(x)=0時, 無相位角差別, 氣動阻尼B為0。

在過渡段上, 復雜的波系結構以及流動分離, 使得氣動力與結構位移之間會出現(xiàn)較為明顯的遲滯現(xiàn)象, 從而導致相位角φ(x)∈(-π, 0), 由此在過渡段上產生了負的氣動阻尼。

計算過程中的廣義力與廣義位移隨時間變化曲線如圖8所示, 可以看出所有工況計算結果都表現(xiàn)良好, 需要注意的是在非定常計算初期, 演化的不完全導致廣義力存在一些突變異常的結果, 計算氣動阻尼時須剔除, 選擇后面演化完全的周期。本文計算了9個周期, 剔除了第1個周期出現(xiàn)的錯誤結果, 采用后8個周期進行氣動阻尼分析。強迫運動振幅為芯級直徑的0.5%。

(a) Ma=0.70

2.3 氣動阻尼影響因素

2.3.1 有無助推對氣動阻尼的影響

捆綁式運載火箭相比于傳統(tǒng)的運載火箭, 最大的區(qū)別就是在尾部四周捆綁了助推器, 使得其流場特性變得復雜, 因此須分析其對氣動阻尼的影響。

圖7、 圖9分別為有無助推時氣動阻尼變化曲線, 可以看出隨著Mach數(shù)的增大整體氣動阻尼先增大后減少, 在Mach數(shù)為0.98時達到最大值, 過渡段與箭體的氣動阻尼變化趨勢與整體基本相同, 而頭部區(qū)域則不同, 是隨著Mach數(shù)的增大一直增大, 只是增長速率變緩。對比兩個圖可知, 助推主要起增大氣動阻尼的作用。還可以看出有無助推情況下頭部的氣動阻尼變化很小, 意味著在箭體尾部施加控制很難影響到頭部的氣動阻尼, 特別是在超聲速流場中。

圖9 無助推時氣動阻尼變化曲線Fig. 9 Aerodynamic damping change curve without boost

2.3.2 前節(jié)點位置影響

為了考察前節(jié)點位置變化對氣動阻尼的影響, 在保持振動頻率不變、 頭部最大振型位置與振幅不變的條件下移動前節(jié)點, 變化后的振型如圖10所示。

(a) Front node after the transition region

根據(jù)對計算結果的分析分別獲得了不同前節(jié)點位置的整體氣動阻尼對比與過渡段氣動阻尼對比, 如圖11、 圖12所示, 可以看出前節(jié)點位置的改變并沒有影響整體氣動阻尼隨Mach數(shù)增大而增大的趨勢, 且前節(jié)點在過渡段上與過渡段前的整體氣動阻尼相差不大, 而前節(jié)點在過渡段后的整體氣動阻尼要高于另兩種情況, 因此過渡段與頭部放在同一側有助于提高氣動阻尼。過渡段的氣動阻尼會隨著前節(jié)點的變化發(fā)生劇烈改變, 前節(jié)點在過渡段前后隨Mach數(shù)增大的變化規(guī)律相反, 節(jié)點前后的振動相位變化導致不同節(jié)點位置過渡段的振動相位不同, 進而導致氣動阻尼發(fā)生變化。

圖11 不同節(jié)點位置的整體氣動阻尼Fig. 11 Overall aerodynamic damping at different node positions

圖12 不同節(jié)點位置的過渡段氣動阻尼Fig. 12 Aerodynamic damping of the transition region at different node positions

2.3.3 強迫振動振幅大小對氣動阻尼的影響

為了考察強迫振動振幅大小對氣動阻尼的影響, 在保證流場結構不發(fā)生改變的前提下, 振動振幅分別為原來的一半和兩倍, 根據(jù)工程經驗, 如果振幅超過芯級直徑的5%, 則須考慮流場結構改變所造成的影響。圖13、 圖14分別為不同振幅下的整體與頭部氣動阻尼。

圖13 不同振幅下整體氣動阻尼Fig. 13 Overall aerodynamic damping at different amplitudes

圖14 不同振幅下頭部氣動阻尼Fig. 14 Aerodynamic damping of the head region at different amplitudes

可以發(fā)現(xiàn)改變振幅無論是對整體氣動阻尼還是頭部氣動阻尼來說變化都很小, 這意味著氣動阻尼的大小主要取決于氣動力與結構振動的相位差, 不依賴于振動幅度的大小。

2.3.4 脈動壓力對氣動阻尼的影響

為了模擬出脈動壓力的影響, 采用IDDES方法對火箭氣動阻尼進行計算, 計算來流Mach數(shù)為0.92, 計算過程中的廣義力與廣義位移如圖15所示, 相較于圖8可以看出廣義力隨時間變化曲線并不光滑, 脈動壓力的存在導致廣義力由多個頻率疊加而成。

圖15 基于IDDES的廣義力與廣義位移變化曲線Fig. 15 Variation cures of generalized force and generalized displacement based on IDDES

由于第2階模態(tài)的頻率為2.46 Hz, 而由分離流、 激波振蕩等引起的脈動壓力頻率往往遠大于此頻率, 因此這里選擇3.5 Hz為分界, 將高于3.5 Hz的部分視為由抖振脈動壓力引起的廣義力, 低于3.5 Hz的部分視為強迫振動引起的廣義力, 通過低通濾波把高于3.5 Hz的廣義力濾掉, 可以獲得由強迫振動引起的廣義力與廣義位移變化曲線, 如圖16所示, 通過此廣義力計算的氣動阻尼為2.08‰。同樣地, 進行高通濾波將低于3.5 Hz的廣義力濾掉, 可以獲得由抖振脈動壓力引起的氣動阻尼為(2.94×10-3)‰, 由此得到脈動壓力引起的氣動阻尼變化為0.14%, 可以忽略不計。同時使用RANS方法計算的氣動阻尼為2.07‰, 與IDDES的計算結果相比誤差約為(2.94×10-3+2.08-2.07)/2.07≈0.48%, 這說明針對氣動阻尼的模擬, 抖振引起的脈動壓力對氣動阻尼的計算結果影響很小, 起主要作用的還是廣義力的變化, 該變化由強迫振動引起的結構邊界變化所導致。

圖16 濾波后的廣義力與廣義位移變化曲線Fig. 16 Variation cures of generalized force and generalized displacement variation curve after filtering

3 結論

本文通過數(shù)值計算方法研究了火箭的氣動阻尼特性。根據(jù)流動特征分析與理論推導, 發(fā)現(xiàn)火箭過渡段幾何外形的收縮導致該區(qū)域出現(xiàn)復雜的分離與激波結構, 從而造成了氣動力相對于結構振動相位的滯后, 導致了該區(qū)域為氣動負阻尼, 即氣動不穩(wěn)定性的主要來源。

在此機理的基礎上, 分析了前節(jié)點位置、 振動振幅、 脈動壓力等因素對氣動阻尼的影響規(guī)律?？梢缘贸鲆韵陆Y論:

1) 助推增加了正阻尼區(qū)域的面積, 從而相對于沒有助推的構型起到了增加氣動阻尼的作用。

2) 前節(jié)點位置的改變對過渡段氣動阻尼影響很大, 節(jié)點前后的振動方向相反, 導致節(jié)點在過渡段前后的氣動阻尼變化規(guī)律也截然相反, 將過渡段與頭部區(qū)域放在節(jié)點的同一側有助于增加氣動阻尼。

3) 在不改變流場結構的前提下, 改變振動的振幅, 氣動力也會產生相應幅度的變化, 因此結構振幅對氣動阻尼的影響可忽略不計。

4) 高頻部分的廣義力對氣動阻尼的貢獻很小, 即結構振動引起的廣義力變化對氣動阻尼起主要作用, 而脈動壓力對計算氣動阻尼影響不大, 可忽略不計。