王司文，韓景龍，員?，|

（1.中國直升機(jī)設(shè)計(jì)研究所直升機(jī)動力學(xué)全國重點(diǎn)實(shí)驗(yàn)室，江西 景德鎮(zhèn) 333001；2.南京航空航天大學(xué)航空學(xué)院航空航天結(jié)構(gòu)力學(xué)及控制全國重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210016）

引言

第一架ABC 共軸驗(yàn)證機(jī)XH-59A［1］由西科斯基公司設(shè)計(jì)。該型直升機(jī)的飛行測試結(jié)果表明，ABC共軸直升機(jī)具有出色的穩(wěn)定性和操縱性，但存在嚴(yán)重的機(jī)體振動問題。從后續(xù)X2 型號［2］和S-97 掠奪者等型號研制過程開始，有效減小直升機(jī)振動水平成為這類直升機(jī)研制的核心技術(shù)之一。

目前已有眾多學(xué)者對ABC 共軸直升機(jī)動力學(xué)進(jìn)行研究，其研究對象多為旋翼［3-6］，但直升機(jī)旋翼與機(jī)身作為一個整體存在強(qiáng)烈的彈性與慣性耦合，僅以旋翼為分析對象不足以滿足全機(jī)振動分析的需求，所以需要綜合考慮機(jī)身與旋翼的耦合關(guān)系，進(jìn)行旋翼/機(jī)身耦合振動分析。目前，建立ABC 共軸直升機(jī)旋翼與機(jī)身的整體模型主要有以下方法。Anusonti-Inthra［7］建立了ABC共軸直升機(jī)旋翼/機(jī)身耦合結(jié)構(gòu)模型，并采用高置信度的CFD/CSD 耦合分析方法對ABC 共軸直升機(jī)旋翼氣彈響應(yīng)、全機(jī)配平計(jì)算進(jìn)行分析。但該文的機(jī)身模型為剛體模型，無法考慮機(jī)身彈性運(yùn)動對旋翼的影響，同時無法計(jì)算機(jī)身特定位置的彈性振動響應(yīng)。Ye 等［8］為計(jì)算XH-59A 機(jī)身振動響應(yīng)，在MSC 軟件中建立了機(jī)身三維有限元模型，在CAMRADⅡ中進(jìn)行旋翼氣動彈性分析，并計(jì)算旋翼槳轂載荷，然后將槳轂載荷通過MSC 軟件以外載荷的形式加載到機(jī)身旋翼軸處，分析機(jī)身振動響應(yīng)。該文給出了ABC 直升機(jī)機(jī)身振動響應(yīng)計(jì)算方法，但該方法并未考慮機(jī)身與旋翼的耦合效應(yīng)，計(jì)算精度有待提高。為解決上述問題，需要建立ABC 共軸直升機(jī)旋翼/機(jī)身耦合氣動彈性分析方法。在目前文獻(xiàn)中還未發(fā)現(xiàn)相關(guān)研究成果，因此ABC 共軸直升機(jī)旋翼/機(jī)身耦合動力學(xué)研究具有重要意義。

為得到可靠的直升機(jī)振動響應(yīng)分析結(jié)果，需要精確的直升機(jī)振動響應(yīng)分析方法。目前，CFD/CSD 耦合是精度最高的直升機(jī)振動響應(yīng)分析方法之一。但目前的直升機(jī)CFD/CSD 耦合方法還存在一些問題亟待解決。傳統(tǒng)CFD/CSD 耦合方法［9-10］的CFD 計(jì)算采用2 階精度和較少的氣動網(wǎng)格，該方法在槳葉附近計(jì)算精度較高，但遠(yuǎn)尾跡的計(jì)算存在數(shù)值耗散問題，這會導(dǎo)致直升機(jī)的槳渦干擾和旋翼/機(jī)身氣動干擾計(jì)算精度不足，最終影響氣動彈性響應(yīng)分析精度。目前的解決辦法多為采用高階精度格式的CFD 方法［11］和數(shù)量巨大的氣動網(wǎng)格［12］。這些改進(jìn)措施對計(jì)算資源與時間的消耗是巨大的。所以，提出一種既可以提高尾跡計(jì)算精度，又可以保證計(jì)算效率的CFD/CSD 耦合計(jì)算方法十分有意義。

為解決上述問題，本文建立了ABC 共軸直升機(jī)旋翼/機(jī)身耦合動力學(xué)模型，并提出了CFD/CSD/自由尾跡耦合分析方法，采用自由尾跡方法提高尾跡計(jì)算精度。綜合上述方法進(jìn)行ABC 共軸直升機(jī)旋翼/機(jī)身耦合振動響應(yīng)分析。本文采用中等變形梁模型對旋翼槳葉進(jìn)行結(jié)構(gòu)建模；利用商用軟件MSC.Patran 創(chuàng)建機(jī)身精細(xì)結(jié)構(gòu)有限元模型，并在MSC.Nastran 軟件中進(jìn)行模態(tài)計(jì)算；基于Euler 控制方程的CFD 方法分析全機(jī)流場；采用自由尾跡方法分析旋翼尾跡，計(jì)算由尾跡引起的旋翼誘導(dǎo)入流，進(jìn)而計(jì)算旋翼準(zhǔn)定常氣動力；最后將上述建模和計(jì)算方法耦合建立旋翼/機(jī)身耦合模型，求解耦合系統(tǒng)配平和振動響應(yīng)，研究了飛行及結(jié)構(gòu)參數(shù)對直升機(jī)重點(diǎn)位置振動響應(yīng)的影響。

1 理論與方法

1.1 共軸直升機(jī)旋翼/機(jī)身耦合氣動彈性模型

本文旋翼槳葉采用15 自由度的中等變形梁單元建模［13］。低精度氣動力采用準(zhǔn)定常氣動力模型和自由尾跡模型。旋翼運(yùn)動方程與機(jī)身方程通過槳轂中心點(diǎn)耦合在一起。

系統(tǒng)的動力學(xué)方程通過Hamilton變分原理得到：

式中δΠ為總變分項(xiàng)，δU為彈性能變分，δT為動能變分，δW為外力虛功，這些能量是所有槳葉和機(jī)身的總和，各能量的變分可表示為：

式中 下標(biāo)“b”和“F”分別表示槳葉和機(jī)身；Nb表示上、下旋翼所有槳葉片數(shù)。上、下旋翼的虛應(yīng)變能求解公式與常規(guī)旋翼相同，由于篇幅原因這里不再給出，可參考文獻(xiàn)［13］，本文僅給出虛動能與外力虛功表達(dá)式的推導(dǎo)過程。

（1）動能

為得到槳葉動能項(xiàng)，需要得到在慣性坐標(biāo)系下上、下旋翼槳葉上任意點(diǎn)的位置。由于ABC 共軸直升機(jī)的特殊構(gòu)型，其槳葉上點(diǎn)的位置表達(dá)式也不相同，分別為：

式中 下標(biāo)“u”代表上旋翼，下標(biāo)“l(fā)”代表下旋翼；TFI，TUI和TDI分別為機(jī)身坐標(biāo)系、未變形坐標(biāo)系和變形坐標(biāo)系到慣性系的轉(zhuǎn)換矩陣，具體公式可參考文獻(xiàn)［12］。{xr，yr，zr}為慣性坐標(biāo)系中下槳轂中心點(diǎn)的剛體位移；{xeu，yeu，zeu+hu}與{xel，yel，zel+hd}分別為上、下旋翼槳轂中心點(diǎn)的彈性位移；{x+ue，v，w}為槳葉未變形坐標(biāo)系上的位置坐標(biāo)，其中，ue，v和w分別為徑向、弦向和垂向位移，η和ζ分別為槳葉截面內(nèi)弦向和垂向位置坐標(biāo)。

根據(jù)文獻(xiàn)［12］公式，求Ru和Rl的變分，得到δRu和δRl；Ru和Rl對時間求二階導(dǎo)，得到au和al。根據(jù)虛動能求解公式，得到槳葉虛動能為：

式中m為質(zhì)量。

（2）氣動力虛功

對上、下旋翼槳葉截面的位置向量表達(dá)式求時間導(dǎo)數(shù)，得到慣性系下上、下旋翼槳葉運(yùn)動速度為：

式中 Ii，Ji和Ki為慣性坐標(biāo)系單位向量。

變形坐標(biāo)系下，槳葉截面相對速度計(jì)算公式為：

式中Vw代表來流速度，μ為前進(jìn)比，{Vixu，Viyu，Vizu}和{Vixl，Viyl，Vizl}為由尾隨渦，遠(yuǎn)尾跡和槳葉附著渦在槳葉的誘導(dǎo)速度。本文采用自由尾跡模型計(jì)算式（8a）～（8b）中的誘導(dǎo)速度。iξu，jηu，kζu和iξl，jηl，kζl分別為上、下旋翼變形坐標(biāo)系單位向量。基于文獻(xiàn)［14］提出的方法計(jì)算得到上、下旋翼尾跡形狀及尾跡渦系環(huán)量。在得到槳葉截面相對速度后，可根據(jù)外力虛功 公式，求解上、下旋翼槳葉外力虛功，具體表達(dá)式可參考文獻(xiàn)［13］。

基于Hamilton 變分原理，將上、下旋翼各自的虛應(yīng)變能、虛動能和外力虛功代入式（1）中，并采用空間有限元離散方法對槳葉能量表達(dá)式進(jìn)行離散，整理得到離散化的槳葉氣動彈性動力學(xué)方程。機(jī)身剛性運(yùn)動方程采用由配平方程得到，機(jī)身彈性運(yùn)動方程由旋翼槳轂力方程計(jì)算得到，綜合上述運(yùn)動方程并對機(jī)身彈性運(yùn)動方程進(jìn)行模態(tài)縮聚，得到共軸直升機(jī)旋翼/機(jī)身耦合運(yùn)動方程為：

式中M，C和K分別為系統(tǒng)的質(zhì)量、阻尼和剛度矩陣；q為物理自由度；p為模態(tài)縮聚后的廣義坐標(biāo)；下標(biāo)“t”，“r”和“e”分別表示與旋翼、機(jī)身剛體運(yùn)動和機(jī)身彈性運(yùn)動自由度；Ft，F(xiàn)r，F(xiàn)e為外力向量。最后，本文采用時間有限元方法對運(yùn)動方程進(jìn)行求解，得到振動響應(yīng)。

1.2 CFD 計(jì)算模型

本文CFD 方法的控制方程基于非定常Euler 方程，采用2階精度的Roe格式，采用重疊網(wǎng)格和動網(wǎng)格方法實(shí)現(xiàn)槳葉運(yùn)動與變形。圖1給出了本文研究的樣例共軸雙旋翼直升機(jī)流場網(wǎng)格劃分，由于采用重疊網(wǎng)格技術(shù)，本文CFD 網(wǎng)格分為氣動網(wǎng)格和旋翼網(wǎng)格，不包含機(jī)身網(wǎng)格。其中氣動網(wǎng)格如圖1 所示，中心網(wǎng)格塊與旋翼網(wǎng)格相重合。本文網(wǎng)格總數(shù)為5831913。

1.3 振動響應(yīng)計(jì)算方法

帶有配平計(jì)算的ABC 共軸直升機(jī)旋翼/機(jī)身耦合氣動彈性響應(yīng)分析的具體步驟如下：

（1）在初始步中，旋翼/機(jī)身耦合氣動彈性方程中的氣動力載荷，采用低置信度氣動力，具體計(jì)算方法見1.1 節(jié)，氣彈方程的氣動載荷記為F0=。對氣彈方程進(jìn)行求解，得到振動響應(yīng)R0與自由尾跡模型，根據(jù)低置信度氣動力計(jì)算方法，僅考慮槳葉附近渦系，計(jì)算得到氣動力，并采用優(yōu)化方法計(jì)算配平參數(shù)T0。

（2）根據(jù)前一步的振動響應(yīng)R0與配平參數(shù)T0，利用CFD 求解器得到氣動載荷氣動力修正量ΔF0=

（3）在第k（k=1，2，3，…）步中，氣動彈性方程中的氣動載荷Fk除了包含全部尾跡的低置信度氣動載荷FLk外，還加入氣動力修正量ΔFk-1，即Fk=+ΔFk-1，求解氣動彈性方程，得到了新的響應(yīng)Rk、氣動力和配平參數(shù)Tk。

（4）將旋翼槳葉振動響應(yīng)Rk和配平參數(shù)Tk帶入CFD 求解器中，得到了新的載荷。并求解新的修正量ΔFk=

（5）若配平結(jié)果不收斂，則繼續(xù)進(jìn)行第（3）～（4）迭代步，直至計(jì)算收斂。當(dāng)配平結(jié)果收斂時，得到了旋翼/機(jī)身耦合系統(tǒng)最終的氣動彈性響應(yīng)和配平結(jié)果。具體分析流程圖如圖2 所示。

圖2 CFD/CSD/自由尾跡耦合分析流程圖Fig.2 CFD/CSD/Free Wake coupled analysis flow chart

2 算例與驗(yàn)證

2.1 共軸雙旋翼自由尾跡計(jì)算模型計(jì)算驗(yàn)證

為驗(yàn)證本文共軸雙旋翼自由尾跡計(jì)算程序的正確性，這里以文獻(xiàn)［15］的試驗(yàn)?zāi)Ｐ蜑樗憷?，模型的具體參數(shù)可參考文獻(xiàn)［15］。

采用本文計(jì)算方法的共軸雙旋翼自由尾跡形狀計(jì)算結(jié)果如圖3 所示，圖中R為旋翼半徑。從圖3 中可以看出，上旋翼的尾跡收縮較嚴(yán)重，被包含在下旋翼尾跡之中；且從圖3（b）中可以看出，上旋翼尾跡下移速度明顯大于下旋翼尾跡下移速度，計(jì)算結(jié)果符合物理意義。圖4給出了上旋翼下方特定位置的誘導(dǎo)速度分布對比，圖中Ω為旋翼轉(zhuǎn)速。從圖4中可以看出，上旋翼下方0.3R和0.4R處的誘導(dǎo)速度計(jì)算值與試驗(yàn)值基本一致，驗(yàn)證了本文的計(jì)算程序的正確性。

圖3 共軸雙旋翼尾跡形狀圖Fig.3 Coaxial twin rotor wake profile

圖4 共軸雙旋翼誘導(dǎo)速度分布圖Fig.4 Induction velocity profile of coaxial rotors

2.2 樣例直升機(jī)模型介紹

本文分析模型為樣例ABC 共軸直升機(jī)，具體模型參數(shù)如表1 所示。

表1 ABC 共軸直升機(jī)模型參數(shù)Tab.1 Model parameters of ABC coaxial helicopter

2.3 槳葉動特性分析

為驗(yàn)證本文旋翼槳葉結(jié)構(gòu)建模的正確性，本文以樣例ABC 共軸直升機(jī)槳葉為算例，計(jì)算了槳葉的固有頻率和模態(tài)振型。表2 給出了算例槳葉截面參數(shù)，表3 給出了槳葉的模態(tài)計(jì)算結(jié)果與試驗(yàn)值，圖5給出了槳葉各階的模態(tài)振型。表3 的對比結(jié)果證明了本文槳葉結(jié)構(gòu)建模方法的正確性。

表2 ABC 共軸直升機(jī)槳葉截面參數(shù)Tab.2 Cross section parameters of ABC coaxial helicopter blades

表3 ABC 共軸直升機(jī)單獨(dú)旋翼槳葉模態(tài)頻率Tab.3 Modal frequencies of single rotor blade of ABC coaxial helicopter

圖5 槳葉各階振型Fig.5 Different modes of blade

2.4 機(jī)身結(jié)構(gòu)模型

本文樣例ABC 共軸直升機(jī)機(jī)身精細(xì)有限元模型在工程軟件MSC.Patran 中進(jìn)行建模。由于缺乏機(jī)身詳細(xì)結(jié)構(gòu)模型，本文基于機(jī)身模態(tài)試驗(yàn)結(jié)果進(jìn)行機(jī)身結(jié)構(gòu)設(shè)計(jì)和結(jié)構(gòu)建模。本文建立的機(jī)身有限元模型如圖6 所示。

圖6 機(jī)身有限元模型Fig.6 Finite element model of fuselage

在MSC.Nastran 軟件中，采用SOL103 模塊進(jìn)行機(jī)身模態(tài)計(jì)算。前5 階固有頻率計(jì)算值與試驗(yàn)值結(jié)果如表4 所示。從對比結(jié)果可以看出本文設(shè)計(jì)的機(jī)身動特性與試驗(yàn)結(jié)果相近，機(jī)身設(shè)計(jì)符合要求。

表4 機(jī)身模態(tài)頻率計(jì)算結(jié)果與試驗(yàn)值Tab.4 Calculation results and test values of airframe modal frequency

2.5 耦合分析方法驗(yàn)證

由于缺乏真實(shí)ABC 共軸直升機(jī)結(jié)構(gòu)模型，本文以S300C 直升機(jī)為算例，對本文所提出的CFD/CSD/自由尾跡分析方法的正確性進(jìn)行驗(yàn)證。分別采用傳統(tǒng)CFD/CSD 耦合方法與本文提出的CFD/CSD/自由尾跡耦合方法對直升機(jī)駕駛員座椅處的3 階諧波振動響應(yīng)進(jìn)行對比分析。S300C 直升機(jī)的具體參數(shù)可參考文獻(xiàn)［16］。圖7 給出了直升機(jī)從懸停到0.3 前進(jìn)比的配平參數(shù)。

圖7 不同前飛速度下的配平參數(shù)Fig.7 Trim parameters at different forward speeds

本文耦合分析方法與傳統(tǒng)分析方法計(jì)算時間對比如表5 所示。從表5 中可以看出，不同于采用高精度CFD 的分析方法會增加數(shù)倍的時間成本，本文方法的計(jì)算時間增加較少。

表5 不同耦合方法時間成本對比Tab.5 Comparison of time cost of different coupling methods

圖8 給出了不同前進(jìn)比的飛行員座椅處不同方向的加速度響應(yīng)。通過采用本文耦合分析方法與傳統(tǒng)CFD/CSD 方法［16］和試驗(yàn)結(jié)果［16］的比較，證明了本文耦合方法的正確性，同時也證明了本文分析方法可顯著提高計(jì)算的準(zhǔn)確性。

圖8 飛行員座椅處加速度響應(yīng)Fig.8 Acceleration response at the pilot seat

3 ABC 共軸直升機(jī)振動響應(yīng)參數(shù)計(jì)算分析

3.1 配平計(jì)算

ABC 共軸直升機(jī)因其特殊的構(gòu)型帶來了更多的操縱配平量，相較于常規(guī)直升機(jī)，還存在差動操縱和輔助推進(jìn)力，這也導(dǎo)致ABC 共軸直升機(jī)的配平參數(shù)個數(shù)大于配平方程數(shù)，常規(guī)配平計(jì)算方法無法使用，所以本文采用優(yōu)化方法進(jìn)行求解，優(yōu)化目標(biāo)是配平量變化最小。

本文以樣例直升機(jī)為算例，計(jì)算在以40 m/s 速度定常前飛狀態(tài)下，帶有輔助推進(jìn)力的直升機(jī)配平參數(shù)。假設(shè)機(jī)身俯仰角為0°，提前操縱角為40°，不考慮垂平尾升力。首先給出配平參數(shù)初值：槳葉總距θ0為7.27°，差動總距Δθ0為0°，橫向周期變距θA1為-3.5°，差動橫向周期變距ΔθA1為0°，縱向周期變距θB1為2.26°，差動縱向周期變距ΔθB1為0°，輔助推進(jìn)力Tpro為926.725 N。

圖9 給出了配平參數(shù)收斂曲線，操縱角的收斂閾值為0.05°，輔助推力的收斂閾值為10 N。圖中θ0表示總距；θA1和θB1分別表示橫/縱向周期變距；dθ0表示差動總距；dθA1和dθB1分別表示差動橫/縱向周期變距；Tpro表示尾槳推力。如圖9 所示，通過7 次迭代計(jì)算，配平參數(shù)收斂，說明本文計(jì)算方法具有較好的收斂性。

圖9 配平變量隨迭代次數(shù)的變化Fig.9 Variation of trim variables with iteration times

3.2 機(jī)身振動響應(yīng)參數(shù)分析

由于本文ABC 共軸直升機(jī)為2×3 片槳葉，所以3 階諧波振動響應(yīng)是振動響應(yīng)的主要組成部分。本文分別計(jì)算了旋翼軸根部及機(jī)身前部的各方向加速度響應(yīng)隨前飛速度和交叉角的變化。

圖10～11 分別給出了不同旋翼交叉角的旋翼軸根部和機(jī)身前部的加速度響應(yīng)隨前飛速度的變化。

圖10 旋翼軸根部加速度3 階諧波量隨前飛速度的變化Fig.10 Variation of rotor shaft root acceleration third-order harmonic with forward speed

圖11 機(jī)身前部加速度3 階諧波量隨前飛速度的變化Fig.11 Variation of 3rd-order harmonic of fuselage front acceleration with forward flight velocity

從圖10～11 中可以看出，旋翼軸根部與機(jī)身前部的振動響應(yīng)隨前飛速度及上、下旋翼交叉角的變化趨勢基本一致。兩部位的加速度響應(yīng)隨著前飛速度的增加呈現(xiàn)先增大再減小而后增大的變化趨勢，這種趨勢與同為上下旋翼2×3 片槳葉的XH-59［9］相 近；兩部位 的加速 度響應(yīng) 在0°交叉角下，縱向、垂向響應(yīng)較60°交叉角下的響應(yīng)要高很多，而不同交叉角下的橫向加速度響應(yīng)大小情況正好相反，該現(xiàn)象與文獻(xiàn)［17］中得到的試驗(yàn)結(jié)論相同。其產(chǎn)生的原因?yàn)椋涸?°交叉角下，上、下旋翼的槳轂振動載荷無相位差，上、下旋翼的縱向和垂向的正方向相同，橫向方向的正方向相反，使得其縱向力和垂向力相互疊加，橫向力相互抵消；在60°交叉角下，上、下旋翼的槳轂振動載荷存在相位差，使得其縱向力和垂向力相互抵消，橫向力相互疊加。

4 結(jié)論

本文提出了CFD/CSD/自由尾跡耦合分析方法，建立了ABC 共軸直升機(jī)旋翼/機(jī)身耦合模型，結(jié)合二者得到ABC 共軸直升機(jī)旋翼/機(jī)身耦合氣動彈性響應(yīng)分析方法。本文首先對分析方法進(jìn)行驗(yàn)證，而后采用該方法進(jìn)行旋翼軸根部及機(jī)身前部隨前飛速度和旋翼交叉角變化的振動響應(yīng)計(jì)算，并對結(jié)果進(jìn)行分析，闡明問題發(fā)生機(jī)理，并得到一些有意義的結(jié)論。具體結(jié)論如下：

（1）通過算例驗(yàn)證了本文CFD/CSD/自由尾跡耦合旋翼/機(jī)身耦合振動響應(yīng)分析方法的正確性；同時，與傳統(tǒng)CFD/CSD 方法和試驗(yàn)結(jié)果進(jìn)行比較，證明了本文分析方法可提高計(jì)算準(zhǔn)確性。

（2）本文分析了不同旋翼交叉角的機(jī)身前部和旋翼軸根部的3 階諧波振動響應(yīng)。計(jì)算結(jié)果表明，在0°交叉角下，3 階諧波振動響應(yīng)的縱向、垂向分量明顯大于橫向；在60°交叉角下，機(jī)身重點(diǎn)部位3 階諧波振動響應(yīng)的縱向、垂向分量明顯小于橫向。

（3）本文分析了不同前飛速度下的機(jī)身重點(diǎn)部位的振動響應(yīng)。計(jì)算結(jié)果表明，振動響應(yīng)值隨前飛速度的增加呈現(xiàn)先增大再變小而后增大的變化規(guī)律。