劉 凱, 苑凱華, 田海濤
(北京機電工程研究所, 北京 100074)
氣動彈性問題一直是飛行器設計中備受關注的重要問題, 在未來的飛行器設計中仍將占有重要地位。現代飛行器設計日益追求高速度、 高機動性, 使得飛行器越來越呈現出輕結構特點, 氣動彈性問題也日漸突出[1]。隨著高速飛行器研究工作的逐漸深入, 飛行器的性能不斷提升, 高速飛行器的結構質量占比需要大幅降低, 使得結構設計人員將面臨巨大的壓力和挑戰(zhàn)[2-8]。顯然, 結構設計的余量越來越小, 結構越來越“纖細”, 導致結構的剛性降低, 飛行器的氣動彈性影響日益突出。
發(fā)散是氣動彈性力學中的一個重要問題, 屬于靜氣動彈性穩(wěn)定性問題。它的現象是: 在一個完全確定的臨界風速下, 彈性升力系統(tǒng)受定常升力的作用, 使其扭轉變形不斷擴大直至破壞[9]。在進行模型臨界發(fā)散風洞試驗時, 可以通過觀察模型在受到干擾后是否能夠保持穩(wěn)定來判斷發(fā)散現象, 即在不斷增加風速的過程中, 分別在各個固定風速下給模型一個擾動, 如果在某個風速下受擾動后的模型變形不再回復到初始平衡位置, 并出現變形擴大, 那么這個動壓即為發(fā)散臨界動壓[10]。
氣動彈性發(fā)散邊界預測方面, 主要集中為前掠翼風洞試驗驗證[11-12], 細長體火箭彈采用工程氣動力結構近似方法獲取彈體的發(fā)散邊界[13], 工程應用方面往往取足夠的安全裕度來保證飛行器不發(fā)生靜氣彈發(fā)散[14]。采用CFD/CSD耦合的方法獲取結構靜發(fā)散邊界研究很少, 本文提出的氣動力/結構耦合的方法可準確預測結構發(fā)散邊界, 能有效地應用于飛行器氣動外形和設計中。
靜氣動彈性計算中, 在氣動力作用下, 結構的平衡方程可寫作[15]
Ku=F
其中,K為結構剛度矩陣,u為結構自由度位移矢量,F為氣動載荷矢量。
依據初始構型獲得的氣動力, 根據公式可以計算出相應的結構變形, 結構形狀改變影響原有的氣動外形, 從而又使氣動力發(fā)生改變, 即氣動力F隨著結構外形u改變而變化。由結構變形引起的氣動力變化須通過重新進行流場計算才能獲得, 所以要求解結構的真實變形, 須采用循環(huán)迭代方法求解,如圖1所示。
圖1 靜氣動彈性計算流程圖Fig. 1 Static aeroelastic calculation flowchart
對于靜氣彈發(fā)散邊界預測方面, 通過計算不同動壓下的靜氣彈收斂歷程。假設動壓Q當前步計算得到最大位移為un, 下一步最大變形位移為un+1, 下兩步最大變形位移為un+2。通過式(1)求得耦合迭代過程位移增長比qs為
(1)
繪制各個動壓Q~qs的曲線, 并利用1階最小二乘法對數據進行處理, 求得qs=1對應的動壓為發(fā)散邊界。
本文計算的模型為矩形安定面, 長度為300 mm, 高度125 mm, 厚度為2.5 mm, 外形如圖2所示。
圖2 安定面外形Fig. 2 Configuration of the stabilizer
安定面結構由兩部分組成, 安定面和安裝底座, 材料均為鋁合金, 底座與安定面采用焊接形式連接, 底座前緣與安定面前緣距離為163 mm, 底座打螺紋孔, 通過螺紋安裝固定。
本次計算網格采用ICEM劃分, 共186個模塊, 網格規(guī)模為3.5×106, 15層附面層, 附面層首層高度為2×10-5m, 附面層總厚度為1.5×10-3m。圖3給出了物面網格和空間拓撲結構。
圖3 物面網格Fig. 3 Surface mesh
數值模擬采用Fluent軟件來求解流場, 控制方程為定常三維N-S方程, 空氣狀態(tài)方程按理想氣體考慮, 空氣黏性系數采用Sutherland公式修正, 湍流模型為標準k-ε模型, 壁面速度分布采用壁面函數修正, 計算時選用基于密度的顯式格式求解。來流邊界給定遠場壓力條件, 固定壁面給定無滑移、 無滲透邊界條件。
結構采用六面體純實體模擬, 網格規(guī)模為8 000, 底座采用SPC約束六自由度。結構有限元求解采用Nastran靜力分析模塊求解, 圖4是結構有限元模型。
圖4 結構有限元模型Fig. 4 Finite element model of the structure
分別計算了Ma=1.5狀態(tài), 攻角0°, 60, 65, 70, 80, 85, 90 kPa這6個動壓下靜氣彈特性, 這里挑選4種列出。在初始數值擾動下, 靜氣彈收斂情況如圖5所示。
(a) 60 kPa
通過式(1)求得不同動壓下的位移增長比, 擬合動壓和位移增長比曲線見圖6, 獲得當前結構的靜氣彈發(fā)散動壓為83.3 kPa。
圖6 不同動壓與位移增長比擬合曲線Fig. 6 Fitted curve of displacement growth ratio at different dynamic pressures
分別計算了Ma=1.5狀態(tài), 攻角-2°, 65 kPa和70 kPa兩個動壓下的靜氣彈特性。在-2°攻角氣動力作用下, 靜氣彈耦合迭代步收斂情況如圖7所示。
圖7 -2°攻角位移收斂曲線Fig. 7 Convergence curve of displacement at -2° angle of attack
65 kPa和70 kPa狀態(tài)安定面位移都呈收斂狀態(tài), 其中65 kPa狀態(tài)位移收斂值為13.2 mm, 70 kPa狀態(tài)收斂值為25.7 mm, 如表1所示。
表1 兩個狀態(tài)位移對比
表2 發(fā)散動壓對比
對比65 kPa狀態(tài)和70 kPa狀態(tài), 動壓增加了7.69%, 位移增長了94.7%。位移增長比例比動壓增長比例大了10倍多。在實際計算狀態(tài), 70 kPa狀態(tài)在耦合迭代收斂后的位移, 安定面大面區(qū)域已經超出鋁合金許用應力, 見圖8。
圖8 70 kPa狀態(tài)位移和應力云圖Fig. 8 Displacement and stress contour map at 70 kPa state
本文的試驗在0.6 m風洞開展, 側壁支撐, 為了避免沖擊, 采用側壁插入機構, 如圖9所示, 待流場穩(wěn)定后插入均勻流場。
圖9 氣彈風洞試驗插入機構Fig. 9 Insertion mechanism for aeroelastic wind tunnel test
試驗在安定面部分粘貼應變片, 測量風洞試驗過程中的應變特性, 如圖10所示。
圖10 風洞試驗應變測量Fig. 10 Strain measurement in wind tunnel test
試驗采用階梯變動壓, 0°攻角, 選取了4個動壓(64.5, 77.3, 80.1, 84.2 kPa), 其中在80.1~84.2 kPa動壓變化過程中安定面出現靜發(fā)散破壞。試驗現場照片如圖11所示。
圖11 0°攻角風洞試驗高速攝像Fig. 11 High-speed camera recording of wind tunnel test at 0° angle of attack
試驗過程中安定面有小幅度振動, 動壓由80.1 kPa增壓至84.2 kPa安定面突然出現破壞, 如圖12所示, 破壞時刻動壓為81.7 kPa。應變信號表現為隨動壓增加應變幅值增加, 在安定面破壞時刻, 應變突增。
圖12 安定面破壞后試驗照片Fig. 12 Experimental photo after destruction of the stabilizer
安定面預置-2°攻角, 如圖13所示, 采用2個階梯變動壓進行試驗, 實際動壓由64.5 kPa到71.5 kPa安定面出現塑性變形, 如圖14所示。
圖13 -2°攻角風洞試驗高速攝像Fig. 13 High-speed camera recording of wind tunnel test at -2° angle of attack
圖14 安定面塑性變形Fig. 14 Plastic deformation of the stabilizer
通過計算值和風洞試驗結果對比, 發(fā)現采用CFD/CSD耦合進行靜氣彈發(fā)散分析具有良好的數值精度。在0°攻角靜發(fā)散對比方面, 偏差小于2%。
同樣在-2°攻角狀態(tài), 風洞試驗和計算均出現結構材料失效的現象, 表明計算與試驗吻合良好。
本文采用氣動力/結構耦合的方法對一種安定面結構靜氣動彈性發(fā)散進行分析, 并開展了相關風洞試驗研究, 得出以下結論:
1) CFD/CSD耦合可準確預測靜氣彈發(fā)散邊界, 本文計算的安定面發(fā)散邊界與試驗值偏差小于2%;
2) 本文選用的安定面模型為靜氣彈不穩(wěn)定結構, 計算和試驗過程中表現為應變與動壓關系呈非線性;
3) 靜氣動彈性不穩(wěn)定結構, 彈性變形造成氣動力進一步增加, 有攻角狀態(tài)首先出現彈性載荷超出結構強度極限導致的破壞, 破壞時動壓離發(fā)散邊界仍有一定余量。