升力式面對稱運載器大攻角顫振分析方法研究

蔡巧言 呂計男 郭力 王昕江

(1 中國運載火箭技術研究院,北京 100076；2 北京宇航系統工程研究所,北京 100076；3 中國航天空氣動力技術研究院,北京 100074）

0 引言

顫振作為復雜的氣動彈性現象，一旦發生往往在幾個周期內就會使飛行器結構發生破壞。國家軍用標準與民航適航條例中對飛行器顫振分析均有明確的要求。

升力式布局面對稱運載器外形復雜，往往由異型截面機身、大翼展機翼、多組氣動舵面組成，飛行剖面具有大空域、寬速域的特點，因此運載器面臨的氣動力環境非常復雜。特別是運載器返回減速需求使得此類飛行器返回飛行段多采用大攻角飛行形式，當運載器大攻角飛行時運載器機身對舵面產生較大干擾，使得舵面前方來流不再表現為自由來流流動結構，特別是運載器后部的V 尾幾乎完全處于前體尾流干擾區內，進一步增加了其繞流流場的復雜性。同時，運載器結構廣泛使用先進復合材料，V 尾彈性頻率低，模態振型復雜。在復雜的流動環境與結構動力學特性共同作用下，V 尾結構易與氣動力耦合發生顫振，因此需要在方案設計初期即對其顫振特性進行準確分析，減小飛行風險。傳統顫振工程計算方法對小攻角條件下氣動控制面的亞聲速/超聲速顫振分析結果可滿足工程初步設計階段的精度要求。對于跨聲速域及復雜擾流環境，工程氣動力模型難以實現非定常氣動力的準確描述。楊超[1]對高超聲速氣動彈性力學研究中的不同工程氣動力模型及CFD/CSD 時域計算方法進行了綜述，給出了不同方法的適用范圍與計算優劣。楊炳淵等[2]針對翼面大攻角顫振給出了一種當地活塞理論計算方法，使用算例證明了該方法在有限厚度翼面大攻角顫振中的可行性。金偉[3]對先進戰斗機全動 V 尾抖振強度進行了研究，使用RANS/LES 混合方法對大迎角飛行狀態下的V 尾擾流流場進行計算，并使用CFD/CSD 方法對其抖振特性進行了分析。對于升力式布局面對稱運載器超聲速大攻角返回飛行而言，受大攻角飛行下大范圍分離非定常流場復雜流動特性影響，傳統工程氣動力方法難以實現流場準確計算。設計中需要采用基于計算流體力學/計算結構動力學（CFD/CSD）耦合求解的數值方法，在精確模擬流場的基礎上研究運載器顫振特性。

本文分別采用基于CFD 和工程氣動力模型的顫振分析方法，對升力式運載器大攻角飛行狀態下V 尾的顫振特性進行計算，對兩種方法計算結果進行對比，分析了工程氣動力方法在此類飛行器大攻角顫振分析中的不足，提出類似問題需要在工程研制中引起關注。

1 顫振分析方法

常用的顫振仿真分析方法有工程氣動力[4,5]方法和CFD/CSD 方法[6-12]。

1.1 工程氣動力方法

時域顫振計算的控制方程可以表示為式（1）所示

式中，M為質量矩陣，D為阻尼矩陣，K為剛度矩陣，u為結構節點位移向量，Ma∞為來流馬赫數，k為減縮頻率，Fa為作用在結構上的氣動力。

根據工程氣動力模型，每個面元上的壓力值與面元下洗的關系可以表示成影響系數的形式，此時Fa可以表示為

式中，Ma∞為來流馬赫數，k為減縮頻率，ρ∞為來流密度，U∞為來流速度。A為氣動力影響系數矩陣，是減縮頻率與來流馬赫數的函數。w為面元下洗向量。實際分析中，一般通過模態坐標投影，對模態坐標下的氣動彈性控制方程進行求解

式中，M,D,K分別為模態質量矩陣、模態阻尼矩陣與模態剛度矩陣。q為廣義坐標向量，Q=φTFa為廣義氣動力向量。φ為模態振型向量組成的振型矩陣。分析中采用基于統一升力面理論（ZONA7U）的工程氣動力模型。通過求解式(3)得到飛行器顫振特性。

1.2 CFD/CSD 耦合分析方法

CFD/CSD 耦合分析方法中CFD 計算模型對V 尾顫振特性的準確預測至關重要。對于沒有氣動干擾的氣動控制面，只需考慮控制面運動產生的非定常氣動力，采用基于雷諾平均（RANS）模式精確模擬近壁面的流場結構即可滿足分析要求。但對于大攻角飛行時處在飛行器前體干擾下的V 尾，僅模擬壁面附近的流場不能準確反映前機身脫體渦對V 尾的影響，需采用更精確的流場模擬數值方法[3,7]。大渦模擬方法（LES）方法能精確求解大范圍的氣流分離和不穩定流動，但由于計算效率低下的問題，不適用于實際工程問題分析。本文采用RANS/LES 的DES（Detached Eddy Simulation）混合算法，分析V 尾的顫振模式。DES 方法在氣動控制面壁面區域具有RANS方法的性能，在流動分離區表現為LES 方程的性能[3,10,12]。V 尾顫振具有模態頻率低、多模態耦合和響應位移大的特點。本文在CFD/CSD 耦合求解方法引入預估矯正迭代方法進一步提高分析精度，采用基于RBF 方法的網格變形技術保證求解魯棒性[8]。CFD/CSD 計算流程示意圖如圖1 所示。

圖1 CFD/CSD 計算流程示意圖Fig.1 CFD/CSD calculation process diagram

2 結構動力學模型

V 尾顫振計算需要考慮所有和V 尾相關的彈性模態。使用Patran 軟件的Cquad4 和Cquad8 殼單元模擬蒙皮，使用殼單元與Cbar 梁單元模擬主要承力結構，使用Conm2 集中質量單元模擬主要質量部件，建立V 尾有限元動力學模型，如圖2所示，設置自由-自由邊界條件。

圖2 V 尾有限元模型Fig.2 FEM model of V-tail

采用Lanczos 方法對有限元動力學模型進行動力學分析，獲得所有和V 尾相關的彈性模態頻率和振型，其中前4 階為V 尾組合模態，第5 階為V 尾扭轉模態。CFD 顫振計算所用三維振型如表1 第3 行所示。使用薄板樣條TPS 方法（Tin Plate Spline Method）將所有和V 尾相關的彈性模態振型向面元網格插值，工程氣動力方法顫振計算所用振型如表1 第4 行所示。

表1 V 尾彈性模態振型Table 1 Modal shapes of V-tail

3 顫振分析

選取大攻角返回特征飛行狀態，馬赫數2.3攻角18°。分別采用工程氣動力/CSD 方法和CFD/CSD 方法開展顫振分析。

3.1 基于工程方法的顫振分析

采用ZONA7U 方法建立超聲速非定常氣動力模型，采用表1 中V 尾彈性模態，在頻域進行顫振分析。分析得到各階模態“動壓-阻尼”曲線與“動壓-頻率”曲線如圖3 所示。

圖3 各階模態“動壓-阻尼/頻率”曲線Fig.3 Dynamic pressure-damping/frequency curve

分析“動壓-阻尼”曲線與“動壓-頻率”曲線可見，工程方法計算得到的顫振動壓為11.3kPa。其中，第4 階模態動壓-阻尼曲線過零。通過模態剔除分析，確定第4 階、第5 階V 尾模態為顫振的主要參與模態。

3.2 基于CFD 方法的顫振分析

相同來流狀態馬赫數2.3、攻角18°，采用CFD 與V 尾彈性模態耦合的分析方法開展V 尾顫振分析。CFD 計算采用DES 模型，結構模態阻尼取0，采用變動壓分析方法開展CFD/CSD 計算。通過各階廣義模態位移的時域響應，判斷顫振動壓與發散模態。來流動壓11.3kPa 時，CFD/CSD 計算結果顯示（圖4）V 尾第4 階、第5 階彈性模態均未發生顫振發散。不斷提高來流動壓進行顫振計算，直到來流動壓達到22.2kPa 時，V尾第4 階、第5 階彈性模態的廣義坐標才出現小幅近似等幅振動。來流動壓22.2kPa 時，V 尾5階廣義模態位移時域響應曲線如圖5 所示。其中，第1～3 階廣義模態位移時域響應均為收斂狀態，第4 階模態呈現緩慢發散趨勢，第5 階模態呈現極限環等幅震蕩趨勢。

圖4 廣義模態位移時域響應（11.3kPa）Fig.4 Generalized displacements in time domain(11.3kPa)

圖5 廣義模態位移時域響應（22.2kPa）Fig.5 Generalized displacements in time domain (22.2kPa)

不同動壓第4 階廣義位移對比如圖6 所示，可見，11.3kPa 來流動壓時，第4 階廣義位移呈現周期性緩慢衰減狀態，阻尼比約0.2%；22.2kPa來流動壓時第4 階廣義位移呈線性小幅增長狀態，阻尼比約-0.4%。在顫振點附近，阻尼呈現非線性特征。

圖6 不同動壓下V 尾廣義位移Fig.6 Generalized displacements in time domain of V-tail

3.3 顫振特性差異分析

工程氣動力方法與CFD 方法預示顫振動壓分別為11.3kPa 和22.2kPa，兩者相差將近一倍，V 尾顫振特性差別較大。對流場進行精確分析，V 尾三個不同站位（Y=0.68/1.29/1.69m）截面流場當地動壓如圖7 所示，其中，Y為飛行器全局坐標系法向。自由來流經運載器頭部脫體激波后，流場壓力上升，速度下降。經過機身產生的膨脹波后流動加速，壓力下降。V 尾附近流動壓縮產生激波，使得V尾附近流場處于局部亞聲速狀態。由于大攻角逆壓梯度存在、激波附面層干擾，在V 尾前部產生大范圍分離，進而影響V 尾附近流場分布，流動結構非常復雜。在左側 V 尾Y=0.68/1.29/1.69m（圖7）的展向站位上，分別取靠近前緣上表面、后緣上表面、前緣下表面與后緣下表面處的4 個特征點，對特征點處的來流動壓進行測量。隨著展向站位Y 值增加，前緣上表面特征點處的當地動壓分別為1.13kPa、17.75kPa與19.32kPa；后緣上表面特征點處的當地動壓分別為6.61kPa、8.86kPa 與7.29kPa；前緣下表面特征點處的當地動壓分別為22.67kPa、18.12kPa 與20.50kPa；后緣下表面當地動壓分別為12.95kPa、14.09kPa 與16.88kPa。三個站位截面當地動壓與自由來流動壓22.2kPa 存在較大差異；此外，隨著展向站位變化，各特征點處當地動壓值的大小也發生變化。

圖7 V 尾典型截面流場云圖Fig.7 Dynamic pressure contour of typical section

對比工程氣動力方法與CFD 方法得到的壓力系數（圖8）。根據CFD 計算結果，左側V 尾背風面的壓力系數與其迎風面相比可以近似視為均勻分布。此時，截取左側V 尾不同展向站位截面，并提取截面前緣、中部與后緣特征點處的表面壓力系數，與背風面對應位置壓力系數相減，可以得到不同特征點處CFD 計算得到的表面壓差系數。從根部特征截面到梢部特征截面，CFD計算左側V尾前緣特征點處的表面壓差系數分別為0.96、0.14 和0.03；中部特征點處的表面壓差系數分別為0.61、0.23 和0.07；后緣特征點處的表面壓差系數分別為0.51、0.28 和0.11。

不同特征點處CFD 計算的壓力系數顯示，左側V尾前緣高壓力系數區域隨著展向站位的增加不斷向后緣擴展，這與工程氣動力方法計算得到的表面壓差系數分布規律存在較大差異。

圖8 為兩種方法得到的V 尾表面壓力云圖。經分析，在計算狀態下，運載器頭部激波、機身壓縮作用，在V 尾前方已出現激波（自身產生激波），使其表面局部出現亞聲速區，并出現分離區，工程氣動力方法已不能反映該流動結構。CFD方法比較準確的預示了當地的流場結構，對展向和弦向流動引起的壓力變化預示準確，較好的描述了前緣高壓區和尾緣及翼尖的低壓區的位置。

不同馬赫數、飛行攻角前體對V 尾的干擾特性不同，V 尾的當地動壓/壓力分布變化規律也不同。采用固定馬赫數（Ma2.3）和固定來流動壓（22.2kPa）變換攻角（8°/18°/28°）的方式研究攻角的影響（圖9）。整體規律上看，在所計算攻角范圍內，V 尾都會發生小阻尼發散。不同的是，隨著攻角增大，V 尾更早進入等幅震蕩狀態，這主要是不同迎角下，前體對V 尾的干擾不同，V尾局部亞聲速區域范圍不同導致。

圖9 攻角對顫振的影響Fig.9 Influence of angle of attack

總之，通過對流場進行精確分析，發現工程氣動力方法與CFD 方法預示顫振動壓產生較大偏差的直接原因是兩種分析方法得到的壓力分布和當地動壓不一致導致，工程方法沒有反映出機身波系及大攻角分離對舵面前方來流流場的干擾影響。

4 結論

本文通過分析升力式布局面對稱飛行器V尾大攻角飛行狀態下的顫振特性，發現工程氣動力頻域方法和CFD/CSD 時域方法在預示中存在較大差異，工程氣動力方法分析結果過于保守，易使總體設計發生誤判，采取額外的剛度補償，增加結構設計代價。通過對升力式布局面對稱運載器大攻角擾流流場開展分析，認為當地動壓和壓力分布的差異是導致工程氣動力和CFD 方法存在差異的主要原因。工程設計中需要辨識飛行剖面內氣動干擾嚴重的區域，在此區域需采用基于CFD 的數值方法開展精細的流動特性分析才能對顫振特性開展準確的預示與評估。