旋翼/機身耦合問題的CFD/CSD分析方法

陳全龍， 韓景龍， 員海瑋

(南京航空航天大學機械結構力學與控制國家重點實驗室， 江蘇 南京 210016)

引 言

直升機乘坐舒適性要求的提高，對減小直升機振動水平提出了更高的要求。直升機的振動主要來自于旋翼系統及其與機身之間的相互耦合。旋翼/機身耦合系統的振動響應分析，是進行減振處理，降低全機振動水平的基礎。

已有許多國外學者對直升機振動響應分析進行了研究。早期的工作大多是分別求出旋翼和機身的阻抗，令兩者的阻抗在槳轂處匹配，進而求得機身的振動響應[1～3]。Warmbrodt和Friedmann最先建立了旋翼/機身耦合系統的非線性氣動彈性響應分析模型，考慮了機身的6個剛體自由度[4]。Fledel通過有限元方法，對旋翼/機身耦合系統的振動響應進行了研究，考慮了旋翼與機體的運動和慣性載荷之間的相互耦合[5]。之后，Friedmann和Chopra等分別建立了旋翼/彈性機身的非線性氣動彈性方程，采用梁單元建立了機身的三維有限元模型，對機身典型位置的振動響應進行了研究[6～9]。

準確的機身結構模型和旋翼流場分析，是直升機振動響應分析的關鍵。上述研究均采用準定常或半經驗的二維氣動力公式，機身則采用剛體模型或基于梁單元的簡化模型，其分析精度還有待進一步提高。

近年來，隨著計算機技術和計算流體力學(CFD)的迅速發展，基于非定常Euler/N-S方程的CFD方法，成為目前最有效的旋翼流場分析方法[10～12]。這為旋翼/機身耦合系統分析精度的提高，提供了新的思路。目前CFD方法已經用于求解單獨旋翼的氣動彈性問題，但對于旋翼/機身耦合系統的響應分析，由于問題難度更大，還未見到相關文獻。

本文將CFD方法引入到旋翼/機身耦合系統的振動響應分析中，通過動態嵌套網格與動網格技術來實現槳葉網格的剛體運動和彈性變形，利用CSD軟件建立精細的機身三維有限元模型，然后通過帶配平的松耦合迭代方法求解系統響應。該方法能準確模擬機身結構和旋翼流場的特性，具有較高的分析精度和工程適用性。以某型兩片槳葉蹺蹺板式直升機為算例，對懸停及定常前飛狀態下，機身典型位置的振動響應進行了分析，將計算結果與實驗值進行了比較，并考察了機身彈性對系統響應分析的影響。

1 計算方法及理論

1.1 旋翼/機身耦合動力學方程

旋翼為兩片槳葉蹺蹺板式構型，槳葉采用15自由度的中等變形梁單元[13]。理論氣動力采用準定常氣動力模型和Drees線性入流模型。旋翼運動方程與機身方程通過槳轂中心點耦合在一起，旋翼載荷通過槳轂傳到機身上，而機身的運動也會影響到旋翼載荷。

在常規的旋翼氣動彈性響應分析中，僅考慮槳轂力常值部分對全機配平方程的影響，在定常前飛時，則假定機身固定不動。但槳轂力的諧波部分會使機身產生剛體運動，進而與旋翼運動發生耦合。因此，在機身振動水平分析中，還必須考慮機身的剛體運動自由度。

系統的動力學方程通過哈密頓變分原理得到

(1)

式中 δU為彈性能變分，δT為動能變分，δW為外力虛功。對這些能量有貢獻的是槳葉和機身，各能量的變分可表示為：

(2)

(3)

(4)

下標b，F分別表示槳葉和機身，Nb表示槳葉片數。

此次分析中，通過MSC/Nastran軟件建立機身的有限元模型，并根據模態試驗結果，對有限元模型進行修正，以保證模型有足夠的分析精度。然后將機身方程與旋翼方程組裝起來，得到旋翼/機身耦合系統的動力學方程。

機身有限元模型如圖1所示，機身主承力結構采用梁單元和殼單元，其他設備質量則采用集中質量單元分布到對應位置處。機身前10階彈性模態的計算值和實驗結果見表1，可以看出，有限元模型的固有模態計算值與實驗結果吻合良好。

圖1 機身有限元模型

表1 機身彈性模態

為減小計算量，對槳葉運動方程和機身彈性運動方程進行了模態縮聚。為避免出現矩陣奇異，將機身剛體運動自由度分開來，放到方程右端。慣性坐標系下，得到旋翼-剛體/彈性機身耦合系統的非線性動力學方程為：

(5)

式中M，C，K為系統的質量、阻尼和剛度矩陣；FL和FNL表示線性和非線性力向量；p是模態縮聚后的廣義坐標。下標b，G，r和e分別表示與槳葉、槳轂萬向鉸、機身剛體運動和機身彈性模態相對應的自由度。旋翼/機身耦合系統方程是一個非線性，周期系數常微分方程組，本文采用時間有限元方法進行求解。機身剛體運動則由配平方程來求解，并通過與式(5)迭代，得到耦合方程的解。

1.2 CFD分析模型

旋翼流場采用非定常Euler/N-S方程進行求解，并通過動網格與動態嵌套網格方法來實現槳葉的彈性變形和剛體運動。為節省計算時間，本文控制方程采用非定常Euler方程。計算中采用雙時間步進行時間推進，每個旋轉周期進行360個物理時間步的計算，每個物理時間步內進行20個偽時間步的迭代計算，對流通量采用Roe格式計算。

圖2(a)是旋翼系統的嵌套網格裝配示意圖，其中背景網格由7個網格塊組成，總網格單元數為1 749 000。單片槳葉網格如圖2(b)所示，槳葉采用C-H型結構網格，由8個子網格塊組成，共171 900個網格單元。兩組槳葉網格和背景網格組成整個嵌套網格系統，網格總單元數為2 092 800。

圖2 旋翼系統網格圖

1.3 分析流程

直升機的配平計算是旋翼氣動彈性分析中不可缺少地重要環節，本文采用自由飛行條件下的6自由度平衡方程，來計算直升機的操縱輸入(包括總距θ75和周期變距θ1c,θ1s)、飛行姿態角(俯仰角αs和滾轉角φs)和尾槳總距θT。

旋翼-剛體/彈性機身耦合系統的響應，采用帶配平的松耦合迭代方法進行求解，CFD和CSD數據每旋轉一周交換一次，具體分析步驟如下：

5．重復步驟3和4直到系統的響應收斂。此時機身的振動響應即為旋翼/機身耦合系統的非線性振動響應。

2 結果與討論

本文以某兩片槳葉蹺蹺板式直升機為算例，直升機的具體參數如表2所示。每片槳葉劃分為6個具有15自由度的中等變形梁單元，各個單元的截面參數見表3。其中，li/R是無因次化的單元長度，mi/m0是無因次化的槳葉線密度，該算例中取參考線密度m0=2.63 kg/m。EIz/m0Ω2R4，EIy/m0Ω2R4和GJ/m0Ω2R4分別為無因次化的擺振、揮舞和扭轉剛度。

表2 直升機參數

表3 槳葉單元參數

2.1 懸停狀態

在懸停狀態下，采用本文所提出的旋翼-剛體/彈性機身耦合系統分析方法，對旋翼軸根部和機身前端點(如圖1所示)的振動響應進行分析，并將計算結果與實驗值進行比較。

圖3是旋翼軸根部振動加速度功率譜密度曲線的計算值和實驗值的對比。可以看出，對于2倍通過頻率(2/rev)和4倍通過頻率(4/rev)的振動響應，計算值與實驗值吻合良好。試驗中有一定的通過頻率響應成分，而計算中沒有，這主要是由旋翼各片槳葉的氣動和質量不平衡造成的，而在計算中沒有考慮這些因素。此外，實驗中還有較大的52～54 Hz頻率響應，這是由發動機(額定轉速3 200 r/min)振動造成的。旋翼軸根部靠近發動機，因此其受發動機振動的影響較大，從圖3(c)中可以看出，當主旋翼所引起的響應較小時，發動機振動和噪聲的影響將明顯增大。

圖4是機身前端點處的加速度功率譜密度曲線，由于此點遠離發動機，因此其振動受發動機的影響明顯減小。

圖3 旋翼軸根部振動自功率譜密度曲線

圖4 機身前端點振動自功率譜密度曲線

圖5 機體縱向振動隨前飛速度的變化

2.2 前飛狀態

在定常前飛狀態下，對機身前端點和旋翼軸根部的振動響應進行計算，研究前飛速度對機身振動水平的影響，并比較考慮和不考慮機身彈性運動對分析結果的影響。

圖5～7是機體的縱向、橫向和垂向加速度響應幅值隨前飛速度的變化曲線。從圖中可以看出，當直升機由懸停轉入前飛狀態時，機體振動水平會隨著前飛速度的增加而迅速增大，并在前進比μ=0.05左右時到達一個局部峰值；當前進比μ在0.05～0.2范圍內時，機體振動水平則基本保持不變；當μ>0.2后，機體的振動水平又會隨著前飛速度的增加而迅速增大。

圖6 機體橫向振動隨前飛速度的變化

圖7 機體垂向振動隨前飛速度的變化

對比圖5～7可以看出，機體的縱橫向振動響應明顯高于垂向振動響應，這說明旋翼的橫向力擾動相對較大。在前進比μ<0.2時，各個自由度均以2/rev的頻率響應分量為主，4/rev頻率的響應相對較小。而當μ>0.2時，4/rev頻率的振動響應所占比重則逐漸增大。

對比僅考慮機身剛體運動和綜合考慮剛體/彈性運動兩種分析模型，可以看出，對旋翼軸根部位置，兩者的分析結果誤差較小；而對機身前端點，兩者的誤差則較大。這是因為旋翼軸彈性相對較小，因此旋翼軸根部振動主要以剛體運動為主；而機身前端點遠離旋翼，且結構彈性較大，因而其振動中彈性運動所占的分量比較大。由此可見，在進行直升機振動水平分析時，機身的彈性影響是不能忽略的。

3 結 論

本文提出了一種直升機旋翼-剛體/彈性機身耦合系統振動響應分析方法，綜合考慮了槳葉的幾何非線性，機身的彈性和剛體運動等。槳葉采用15自由度中等變形梁單元，通過基于動網格和動態嵌套網格技術的非定常Euler/N-S方程，與準定常氣動力模型相結合的方法來求解彈性旋翼流場，利用CSD軟件建立精細的機身三維有限元模型，最后用松耦合迭代方法對系統的響應進行分析。該方法能準確模擬機身結構動力學和旋翼流場的特性，具有很好的分析精度和工程適用性，可廣泛應用于直升機振動水平分析和減振設計。

以某小型直升機為例，通過MSC/Patran軟件建立精細的機身三維有限元模型，并根據模態試驗對模型進行修正。計算了懸停狀態下機身典型位置的振動響應，計算結果與實驗值吻合良好。研究了前飛速度對機身振動水平的影響，結果表明，在小速度和大速度前飛時，機身振動水平隨前飛速度的增加而增大，在中等速度段，機體振動水平基本保持不變。

分析中采用了僅考慮機身剛體運動和綜合考慮剛體/彈性運動兩種模型，結果表明，在旋翼軸根部，機身振動以剛體運動為主，而在機身前端處的機體振動響應中，結構彈性的影響非常大。因此，在進行直升機振動水平分析時，必須考慮機身結構彈性的影響。

參考文獻：

[1] Gerstenberger W, Wood E R. Analysis of helicopter aeroelastic characteristics in high speed flight[J]. AIAA Journal, 1963, 1(10): 2 366—2 381.

[2] Staley J A, Sciarra J J. Coupled rotor/airframe vibration prediction methods[R]. NASA SP-352, 1974, 81—90.

[3] Hsu T K, Peters D A. Coupled rotor/airframe vibration analysis by a combined harmonic-balance impedance-matching method [J]. Journal of the American Helicopter Society, 1982, 27(1): 25—34.

[4] Warmbrodt Friedmann. Formulation of coupled rotor/fuselage equations of motion [J]. Vertica, 1979, 3: 245—271.

[5] Fledel. Coupled rotor/airframe vibration analysis[D]. Aerospace Engineering, University of Maryland, 1989.

[6] Thiem Chiu, Peretz P Friedmann. A coupled helicopter rotor/fuselage aeroelastic response model for ACSR[A]. AIAA Conference[C]. New Orleans,USA.1995, 574—600.

[7] Richard C Cribbs, Peretz P Friedmann, Thiem Chiu. Coupled helicopter rotor/flexible fuselage aeroelastic model for control of structural response[J]. AIAA Journal, 2000, 38(10):1 777—1 788.

[8] Senthilvel Vellaichamy, Inderjit Chopra. Aeroelastic response of helicopters with flexible fuselage modeling[A]. AIAA Conference[C]. Dallas,USA.1992： 2 015—2 026.

[9] Hyeonsoo Yeo, Inderjit Chopra. Coupled rotor/fuselage vibration analysis for teetering rotor and test data comparison[J]. Journal of Aircraft, 2001, 38(1):111—121.

[10] Klausdieter Pahlke, Berend G， van der Wall. Chimera simulations of multibladed rotors in high-speed forward flight with weak fluid-structure-coupling [J]. Aerospace Science and Technology, 2005, 9: 379—389.

[11] Mark Potsdam, Hyeonsoo Yeo, Wayne Johnson. Rotor airloads prediction using loose aerodynamic/structural coupling [J]. Journal of Aircraft, 2006，43(3)：732—742.

[12] Jasim Ahmad, Robert T Biedron. Code-to-code comparison of CFD/CSD simulation for a helicopter rotor in forward flight[A]. AIAA Applied Aerodynamics Conference[C]. Honolulu, USA. 2011.