馬礫， 招啟軍， 趙蒙蒙, 王博

南京航空航天大學(xué) 直升機(jī)旋翼動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室， 南京 210016

基于CFD/CSD耦合方法的旋翼氣動(dòng)彈性載荷計(jì)算分析

馬礫， 招啟軍*， 趙蒙蒙, 王博

南京航空航天大學(xué) 直升機(jī)旋翼動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室， 南京 210016

為提高直升機(jī)前飛狀態(tài)下旋翼非定常氣動(dòng)彈性載荷的預(yù)估精度，在旋翼氣動(dòng)彈性綜合分析方法中引入旋翼CFD模塊，建立了一套基于 CFD/CSD松耦合分析的計(jì)算方法和程序。為高效解決流固耦合方法中由于槳葉揮舞、扭轉(zhuǎn)等彈性變形帶來的旋翼貼體網(wǎng)格變形問題，采用基于代數(shù)變換方法的網(wǎng)格變形技術(shù)，槳葉運(yùn)動(dòng)變形量和旋翼氣動(dòng)力信息通過流固交接面?zhèn)鬟f。旋翼流場(chǎng)分析方法中，主控方程采用耦合S-A湍流模型的Navier-Stokes方程，圍繞旋翼流場(chǎng)的網(wǎng)格采用結(jié)構(gòu)嵌套網(wǎng)格方法生成，無黏通量計(jì)算采用Roe格式，時(shí)間推進(jìn)采用雙時(shí)間法。旋翼結(jié)構(gòu)分析中，考慮旋翼配平，基于Hamilton變分原理和20自由度Timoshenko梁模型求解彈性旋翼非線性運(yùn)動(dòng)方程。分別對(duì)CSD和CFD方法進(jìn)行驗(yàn)證，在此基礎(chǔ)上，計(jì)算了SA349/2旋翼槳葉在前飛狀態(tài)下的非定常氣動(dòng)力、揮舞彎矩、擺振彎矩和扭轉(zhuǎn)力矩，并與飛行測(cè)試數(shù)據(jù)進(jìn)行了對(duì)比。計(jì)算表明： CFD/CSD耦合方法可以顯著提高旋翼非定常氣動(dòng)彈性載荷的分析精度，精確捕捉槳葉表面壓強(qiáng)峰值、激波位置等，表明本文發(fā)展的旋翼CFD/CSD耦合方法可以有效地運(yùn)用到旋翼氣動(dòng)彈性載荷的預(yù)測(cè)分析中。

旋翼； 氣動(dòng)彈性載荷； CFD/CSD耦合方法； Navier-Stokes方程； 網(wǎng)格變形技術(shù)； 前飛狀態(tài)

旋翼是直升機(jī)主要的氣動(dòng)力來源，也是其主要振源之一，準(zhǔn)確獲取旋翼非定常氣動(dòng)彈性載荷是進(jìn)行直升機(jī)減振優(yōu)化、強(qiáng)度校核及氣彈穩(wěn)定性研究工作的前提[1]。對(duì)旋翼氣動(dòng)彈性載荷的準(zhǔn)確預(yù)估也是旋翼氣動(dòng)彈性分析領(lǐng)域中最具挑戰(zhàn)性的難題之一。一方面，直升機(jī)旋翼工作在復(fù)雜三維非定常氣動(dòng)環(huán)境中，非定常氣動(dòng)載荷特征明顯；另一方面，由于直升機(jī)旋翼槳葉是細(xì)長(zhǎng)柔性體，彈性槳葉在旋轉(zhuǎn)過程中同時(shí)存在揮舞、擺振、扭轉(zhuǎn)及軸向拉伸等多自由度的彈性變形運(yùn)動(dòng)，受非定常氣動(dòng)載荷的影響，這些運(yùn)動(dòng)之間存在非線性的耦合關(guān)系，從而產(chǎn)生復(fù)雜的結(jié)構(gòu)(彈性)載荷，槳葉產(chǎn)生結(jié)構(gòu)變形，繼而影響氣動(dòng)力的重新分布。因此，旋翼槳葉的氣動(dòng)載荷和結(jié)構(gòu)運(yùn)動(dòng)響應(yīng)是相互依賴和相互耦合的。要準(zhǔn)確預(yù)估直升機(jī)旋翼槳葉的氣動(dòng)彈性載荷，需要綜合空氣動(dòng)力學(xué)、結(jié)構(gòu)動(dòng)力學(xué)和數(shù)值計(jì)算方法等多學(xué)科知識(shí)，構(gòu)建合理的旋翼氣動(dòng)彈性綜合分析方法[2]。

在早期旋翼氣動(dòng)彈性綜合分析中，旋翼氣動(dòng)力的計(jì)算通常采用滑流理論、葉素理論等簡(jiǎn)單模型，導(dǎo)致氣動(dòng)性能計(jì)算精度較差，后來的自由尾跡方法雖然可以有效獲得旋翼流場(chǎng)的誘導(dǎo)速度分布，提高了旋翼氣動(dòng)特性的計(jì)算精度，但在模擬由于黏性引起的流動(dòng)分離現(xiàn)象及氣動(dòng)阻力時(shí)存在局限性。隨著計(jì)算流體力學(xué)(CFD)技術(shù)的發(fā)展和計(jì)算機(jī)性能的大幅提高，以考慮黏性影響的Navier-Stokes方程為主控方程的CFD方法為準(zhǔn)確分析旋翼流動(dòng)細(xì)節(jié)及其氣動(dòng)特性提供了一個(gè)新途徑。因此，將先進(jìn)的CFD方法引入旋翼氣動(dòng)彈性綜合分析中，構(gòu)建旋翼CFD/計(jì)算結(jié)構(gòu)動(dòng)力學(xué)(CSD)耦合方法，提高旋翼氣動(dòng)彈性載荷的模擬精度已成為當(dāng)前旋翼氣動(dòng)彈性綜合分析領(lǐng)域的最新熱點(diǎn)[3]。

在旋翼CFD/CSD耦合分析方面，1986年，Tung和Caradonna等[4]首次提出了旋翼CFD/CSD松耦合迭代方法，該方法得到了廣泛認(rèn)可。2006年，Potsdam等[5]利用松耦合方法將OVERFLOW-D 旋翼CFD求解器和CAMRAD-II旋翼CSD求解器耦合，對(duì)UH-60A旋翼在高速、低速、動(dòng)態(tài)失速3種飛行狀態(tài)下的載荷進(jìn)行了預(yù)估，計(jì)算結(jié)果表明這種耦合方法對(duì)旋翼氣動(dòng)力和力矩有較高的預(yù)測(cè)精度。2013年，Zhao和He[6]采用基于Lagrange描述的黏性渦粒子(VVPM)與CFD/CSD程序結(jié)合，對(duì)UH-60A旋翼槳葉結(jié)構(gòu)特性進(jìn)行了模擬，盡管提高了旋翼在BVI狀態(tài)下結(jié)構(gòu)載荷的計(jì)算精度，但與試驗(yàn)值仍有一定誤差。2015年，Lim[7]耦合了基于Navier-Stokes方程的OVERFLOW2軟件和CAMRAD-II軟件，考慮多段槳尖影響，對(duì)旋翼結(jié)構(gòu)特性和氣動(dòng)特性嘗試了優(yōu)化設(shè)計(jì)并取得了初步結(jié)果。

國(guó)內(nèi)在運(yùn)用CFD/CSD耦合方法進(jìn)行旋翼氣動(dòng)彈性分析研究方面起步較晚，仍處于發(fā)展階段。2010年，王海[8]嘗試進(jìn)行了基于Euler方程的CFD/CSD耦合分析，并進(jìn)行了簡(jiǎn)單的驗(yàn)證。2015年，王俊毅等[9]發(fā)展了基于Navier-Stokes方程的直升機(jī)前飛狀態(tài)旋翼氣動(dòng)特性分析的CFD/CSD耦合方法，湍流模型采用了零方程的B-L模型，提高了旋翼氣動(dòng)載荷的預(yù)估精度，但針對(duì)旋翼彈性載荷及氣動(dòng)性能的研究尚未開展。綜上所述，旋翼CFD/CSD分析方法中槳葉彈性變形對(duì)旋翼流場(chǎng)和非定常氣動(dòng)彈性載荷帶來的影響研究還不夠深入，為此需要提高旋翼非定常氣動(dòng)特性及彈性載荷的模擬精度，這對(duì)CFD方法的數(shù)值離散格式、湍流模型、氣動(dòng)彈性耦合策略和網(wǎng)格變形方法等都提出了很高的要求。

鑒于此，基于先進(jìn)的CSD和CFD方法，建立了適合前飛非定常狀態(tài)下直升機(jī)旋翼CFD/CSD耦合數(shù)值模擬方法。為了準(zhǔn)確模擬黏性大分離流動(dòng)，流場(chǎng)控制方程采用雷諾平均Navier-Stokes(RANS)方程，湍流模型采用一方程的S-A湍流模型，針對(duì)耦合過程中由于槳葉彈性變形帶來的網(wǎng)格變形問題，發(fā)展了一套高效快速的網(wǎng)格變形方法、彈性網(wǎng)格挖洞方法和貢獻(xiàn)單元搜索方法。應(yīng)用所建立的方法，分別對(duì)旋翼結(jié)構(gòu)動(dòng)力學(xué)(CSD)模塊和剛性旋翼流場(chǎng)分析(CFD)模塊進(jìn)行了驗(yàn)證。采用CFD/CSD耦合方法，通過對(duì)不同飛行狀態(tài)下SA349/2旋翼槳葉的氣動(dòng)載荷、結(jié)構(gòu)載荷以及流場(chǎng)的數(shù)值模擬，對(duì)前飛狀態(tài)下彈性槳葉旋翼流場(chǎng)氣流分離、激波失速、槳-渦干擾等現(xiàn)象進(jìn)行了細(xì)致分析，揭示了旋翼非定常氣動(dòng)彈性載荷變化特性，得到了有意義的結(jié)論。

1 數(shù)值方法

1.1 CSD方法

1.1.1 結(jié)構(gòu)動(dòng)力學(xué)計(jì)算方法

采用Hamilton變分原理建立旋翼槳葉動(dòng)力學(xué)方程[10]，在非保守力學(xué)系統(tǒng)中，從時(shí)刻t1到t2的運(yùn)動(dòng)中，系統(tǒng)的勢(shì)能和動(dòng)能與外力虛功的積分為零，其數(shù)學(xué)表達(dá)式為

(1)

式中：δU、δT和δW分別表示某一時(shí)刻槳葉的應(yīng)變能的變分、動(dòng)能的變分和外力所做的虛功的變分。

圖1 槳葉單元節(jié)點(diǎn)自由度排列 Fig.1 Arrangement of degree of freedom of blade element node

將用插值型函數(shù)離散得到的相應(yīng)位移代入變分表達(dá)式，可以得到

(2)

式中：q為單元自由度向量；KL和KNL分別為線性和非線性剛度矩陣；M為對(duì)稱質(zhì)量矩陣；C為非對(duì)稱的科氏阻尼矩陣；KCF為離心力剛度矩陣；FCF為離心力向量；KI為分布力矩所導(dǎo)致的剛度型矩陣；FI為外力向量。

采用有限元法對(duì)式(2)單元陣進(jìn)行組裝，最后得到整個(gè)槳葉的運(yùn)動(dòng)方程：

(3)

本文采用改進(jìn)Newmark-Beta方法[11]，對(duì)旋翼槳葉運(yùn)動(dòng)方程求解。

1.1.2 槳葉模態(tài)分析

為了確定槳葉的頻率和振型，必須對(duì)槳葉進(jìn)行固有特性分析，這也是為了驗(yàn)證所建立的槳葉結(jié)構(gòu)數(shù)值分析模型能否較為準(zhǔn)確地模擬真實(shí)旋翼的動(dòng)力學(xué)特性，以便用于槳葉氣彈耦合的分析中。模態(tài)分析中，假定槳葉不受氣動(dòng)外力的干擾，即認(rèn)為是真空狀態(tài)下的模態(tài)分析，則動(dòng)力學(xué)方程可簡(jiǎn)化為

(4)

1.1.3 槳葉網(wǎng)格彈性變形策略

為兼顧計(jì)算效率，本文采用高效、高質(zhì)量的代數(shù)變換方法來實(shí)現(xiàn)槳葉貼體網(wǎng)格變形[12]。通過CSD方法計(jì)算得到槳葉每一個(gè)方位角和剖面的彈性位移，便可以得到坐標(biāo)轉(zhuǎn)換矩陣Tdc和線性變形陣dlin，再針對(duì)槳葉網(wǎng)格所有坐標(biāo)(x,y,z)進(jìn)行代數(shù)變換，通過式(5)計(jì)算得到槳葉固連坐標(biāo)系中的變形網(wǎng)格坐標(biāo)：

(5)

式中：(xd,yd,zd)為槳葉網(wǎng)格變換后坐標(biāo)。

圖2為SA349/2旋翼槳葉貼體網(wǎng)格變形示意圖，其中下方為變形前網(wǎng)格，上方為變形后網(wǎng)格。

圖2 槳葉網(wǎng)格變形示意圖Fig.2 Schematic of blade grid deformation

1.2 CFD方法

1.2.1 網(wǎng)格生成方法

旋翼運(yùn)動(dòng)復(fù)雜，在旋轉(zhuǎn)過程中伴隨有揮舞、擺振和周期變距等運(yùn)動(dòng)。因此，本文使用了運(yùn)動(dòng)嵌套網(wǎng)格方法對(duì)旋翼流場(chǎng)進(jìn)行數(shù)值模擬。本文的槳葉貼體網(wǎng)格和背景網(wǎng)格均采用結(jié)構(gòu)網(wǎng)格方法。如圖3所示，槳葉貼體結(jié)構(gòu)網(wǎng)格采用C-O型網(wǎng)格。首先，采用求解Poisson方程[13]生成圍繞二維翼型剖面的C型網(wǎng)格，為提高黏性計(jì)算的精度，采用Higenstock源項(xiàng)修正策略調(diào)整網(wǎng)格疏密和正交程度，為了準(zhǔn)確捕捉附面層內(nèi)流動(dòng)，第1層翼型網(wǎng)格的壁面距離取為5×10-6c(y+≈1)，c為翼型弦長(zhǎng)。然后，根據(jù)槳葉外形的翼型分布和扭轉(zhuǎn)分布規(guī)律，沿展向?qū)型網(wǎng)格進(jìn)行插值來獲得槳葉段間的網(wǎng)格，槳根和槳尖處網(wǎng)格采用繞翼型中弧線翻折的方法生成。

圖3 SA349/2旋翼槳葉網(wǎng)格示意圖Fig.3 Schematic of SA349/2 rotor blade grid

前飛計(jì)算時(shí)，旋翼向前運(yùn)動(dòng)，渦向后脫出，則采用了長(zhǎng)方體笛卡兒網(wǎng)格。對(duì)于不同的計(jì)算狀態(tài)，為準(zhǔn)確模擬旋翼槳尖渦的發(fā)展軌跡，對(duì)包含旋翼尾跡區(qū)域進(jìn)行了加密(見圖2)，利于有效計(jì)算流場(chǎng)速度空間分布等信息。在本文的旋翼CFD計(jì)算中，槳葉貼體網(wǎng)格大小為237×45×61(弦向×法向×展向)，背景網(wǎng)格大小為160×164×152(弦向×法向×展向)。

1.2.2 運(yùn)動(dòng)嵌套網(wǎng)格方法

在彈性槳葉流場(chǎng)計(jì)算過程中，由于槳葉每個(gè)時(shí)間步都存在彈性變形，旋翼旋轉(zhuǎn)時(shí)的每個(gè)方位角都必須進(jìn)行背景網(wǎng)格洞單元識(shí)別和槳葉網(wǎng)格貢獻(xiàn)單元搜索。本文采用高效的“擾動(dòng)衍射”[14]法來快速確定槳葉在背景網(wǎng)格中的洞邊界。采用“最小距離方法”[15]進(jìn)行背景網(wǎng)格人工內(nèi)邊界的貢獻(xiàn)單元搜索，進(jìn)而完成旋翼網(wǎng)格與背景網(wǎng)格之間的流場(chǎng)信息傳遞。圖4給出了本文所采用的運(yùn)動(dòng)嵌套網(wǎng)格系統(tǒng)、洞邊界單元以及搜索到的貢獻(xiàn)單元示意圖。

圖4 運(yùn)動(dòng)嵌套網(wǎng)格系統(tǒng)示意圖Fig.4 Schematic of moving-embedded grid system

1.2.3 流場(chǎng)計(jì)算方法

為精確模擬旋翼非定常流動(dòng)特性，捕捉槳葉表面激波位置、流動(dòng)分離等現(xiàn)象，采用積分形式的RANS方程作為旋翼貼體網(wǎng)格流場(chǎng)求解的控制方程[16]：

(6)

式中：V和S分別為控制體體積和表面積;n為控制體表面法矢；守恒變量、無黏通量和黏性通量分別為

其中：Vr為相對(duì)速度；Vt為牽連速度；p為壓強(qiáng)；ρ為流體單元密度；E為單元總內(nèi)能；[u,v,w]為氣流速度；H為單元總焓能；τij和θi分別為黏性應(yīng)力張量項(xiàng)和熱通量項(xiàng)，與溫度導(dǎo)數(shù)相關(guān)。

采用S-A[17]湍流模型，其計(jì)算量小，對(duì)計(jì)算復(fù)雜流動(dòng)有很強(qiáng)的魯棒性，且湍流的渦黏度場(chǎng)是連續(xù)的。在旋翼非定常流場(chǎng)計(jì)算中，本文采用雙時(shí)間法進(jìn)行求解。為了提高流場(chǎng)求解的效率，偽時(shí)間推進(jìn)方法采用LU-SGS隱式格式[18]。首先對(duì)控制方程中的時(shí)間導(dǎo)數(shù)進(jìn)行離散后得

(7)

式中：m為推進(jìn)層數(shù)；Res為殘值；Q=[ρ,u,v,w,p]為原始變量。然后采用LU-SGS隱式格式對(duì)式(7)進(jìn)行求解。

空間離散采用有限體積法，采用Roe格式[19]計(jì)算無黏通量，具體為

|ARoe|i+1/2(WR-WL)]

(8)

式中：ARoe為平均雅可比矩陣，該矩陣考慮了近似Riemann問題的波的大小和傳播方向，下標(biāo)i+1/2表示單元交界面；WL和WR分別為交界面左右兩側(cè)的守恒量。

1.3 CFD/CSD耦合方法

本文采用松耦合[20-21]方式實(shí)現(xiàn)CFD程序與CSD程序的耦合，主要分為2個(gè)步驟：① 氣動(dòng)信息向結(jié)構(gòu)計(jì)算的傳遞，主要通過傳遞CFD計(jì)算得到的氣動(dòng)力來實(shí)現(xiàn)。CFD程序在計(jì)算至少一周期以后獲得一周期的升力、阻力和俯仰力矩，再代入結(jié)構(gòu)計(jì)算程序(CSD)進(jìn)行計(jì)算。② CFD氣動(dòng)計(jì)算中計(jì)入槳葉彈性變形作用，需要對(duì)網(wǎng)格進(jìn)行彈性變形，由于CFD計(jì)算采用時(shí)間推進(jìn)求解，因而需要在每步計(jì)算前對(duì)當(dāng)前站點(diǎn)上的網(wǎng)格進(jìn)行變形。本文松耦合計(jì)算的流程如圖5所示。

步驟1分別對(duì)旋翼CFD模塊和CSD模塊進(jìn)行初始化。

步驟2對(duì)采用基于簡(jiǎn)單升力線模型的旋翼CSD程序進(jìn)行計(jì)算，得到操縱量以及一個(gè)槳葉旋轉(zhuǎn)周期上每個(gè)方位角收斂的槳葉彈性軸處的變形量(位移和扭轉(zhuǎn)角)。

步驟3將一個(gè)周期的槳葉彈性變形量導(dǎo)入網(wǎng)格變形模塊，并完成槳葉貼體網(wǎng)格變形。

步驟4將變形后的槳葉貼體網(wǎng)格代入流場(chǎng)CFD程序，進(jìn)行下一個(gè)周期的時(shí)間推進(jìn)計(jì)算。

步驟5將一個(gè)周期槳葉(表面網(wǎng)格)上的非定常氣動(dòng)力和力矩通過插值傳遞到槳葉結(jié)構(gòu)結(jié)點(diǎn)上，進(jìn)行下一個(gè)周期的配平和槳葉結(jié)構(gòu)響應(yīng)(CSD)計(jì)算。

步驟6重復(fù)步驟2～步驟5的過程，直到旋翼配平量、槳葉結(jié)構(gòu)響應(yīng)和氣動(dòng)載荷收斂為止，輸出流場(chǎng)及結(jié)構(gòu)載荷等結(jié)果。

圖5 旋翼配平及氣動(dòng)彈性耦合流程圖Fig.5 Flowchart of rotor trim and aeroelastic coupling

2 模型驗(yàn)證

2.1 CSD模型驗(yàn)證

為驗(yàn)證本文發(fā)展的旋翼CSD方法的有效性，對(duì)文獻(xiàn)[10]給出的樣例旋翼模型進(jìn)行驗(yàn)證。圖6給出了本文計(jì)算得到的旋翼共振圖。圖中：1F、1L和1T分別表示一階揮舞、一階擺振和一階扭轉(zhuǎn)模態(tài)；Ω為轉(zhuǎn)速。可以看出，CSD模塊可以準(zhǔn)確地計(jì)算旋翼槳葉的固有特性，各階頻率與試驗(yàn)值符合得都較好。

圖7為模型旋翼在前進(jìn)比μ為0.3時(shí)的響應(yīng)分析與文獻(xiàn)[10]計(jì)算比。圖中：φ為方位角?？梢钥闯?，本文發(fā)展的旋翼CSD方法可以有效地計(jì)算前飛狀態(tài)下旋翼的氣動(dòng)彈性響應(yīng)，表明旋翼CSD模塊為CFD/CSD耦合計(jì)算提供了一個(gè)良好的動(dòng)力學(xué)分析模型。

圖6 模型旋翼槳葉共振Fig.6 Resonance of model rotor blade

圖7 前飛狀態(tài)旋翼不同方位角槳尖響應(yīng)Fig.7 Response of blade-tip with different azimuth angles in forward flight

2.2 CFD模型驗(yàn)證

本文計(jì)算了前飛無升力狀態(tài)下剛性Caradonna-Tung旋翼[22]的槳葉表面壓強(qiáng)系數(shù)Cp來驗(yàn)證旋翼CFD模型的有效性，計(jì)算狀態(tài)參數(shù)為：前進(jìn)比為0.2，槳尖馬赫數(shù)為0.8。計(jì)算結(jié)果和試驗(yàn)值的對(duì)比如圖8所示，R為旋翼半徑?？梢姡?jì)算結(jié)果與試驗(yàn)值吻合較好，表明本文的旋翼非定常流場(chǎng)計(jì)算方法能夠有效地用于CFD/CSD耦合分析中的流場(chǎng)計(jì)算部分。

圖8 前飛狀態(tài)Caradonna-Tung旋翼槳葉表面的壓強(qiáng)系數(shù)分布Fig.8 Distribution of pressure coefficient on blade surface of Caradonna-Tung rotor in forward flight

3 旋翼氣動(dòng)彈性載荷分析

3.1 旋翼槳葉模態(tài)分析

分別采用剛性旋翼CFD方法、自由尾跡方法[23]和CFD/CSD耦合方法對(duì)旋翼進(jìn)行氣動(dòng)彈性載荷和旋翼流場(chǎng)分析。采用具有詳細(xì)結(jié)構(gòu)參數(shù)的SA349/2旋翼[24]作為驗(yàn)證算例。SA349/2旋翼基本參數(shù)如表1所示，其中的結(jié)構(gòu)參數(shù)為UMARC的輸入值。

表2為本文計(jì)算得到的SA349/2旋翼固有頻率和CAMRAD軟件[23]計(jì)算結(jié)果的對(duì)比?？梢钥闯觯?種方法吻合較好。

圖9為本文計(jì)算得到的SA349/2旋翼共振與CAMRAD軟件計(jì)算結(jié)果對(duì)比。可以看出，槳葉的固有頻率隨著轉(zhuǎn)速的增大而增加，在工作轉(zhuǎn)速下槳葉固有頻率與CAMRAD軟件計(jì)算值吻合較好。

表1 SA349/2旋翼基本參數(shù)[24]Table 1 Basic parameters of SA349/2 rotor[24]

表2SA349/2旋翼固有頻率計(jì)算結(jié)果

Table2CalculationresultsofnaturalfrequencyofSA349/2rotor

ModeNaturalfrequency/HzCAMRAD[23]CalculatedRelativeerror/%1L0．590．58131．471F1．021．02110．1082F2．782．77580．151T4．164．17170．281

圖9 SA349/2旋翼槳葉共振Fig.9 Resonance of SA349/2 rotor blade

圖10為SA349/2旋翼槳葉的振型圖，r為剖面沿槳葉展向位置；R為槳葉半徑?？梢钥闯觯疚牟捎玫腃SD模型可以有效捕捉彈性旋翼槳葉的模態(tài)。

表3給出了本文計(jì)算的高速和低速2種不同飛行狀態(tài)時(shí)的參數(shù)。表中：CT為旋翼拉力系數(shù)，σ為旋翼實(shí)度。

圖10 SA349/2旋翼槳葉振型Fig.10 Vibration model of SA349/2 rotor blade

表3 SA349/2旋翼飛行試驗(yàn)狀態(tài)基本參數(shù)Table 3 Basic parameters of SA349/2 rotor flight test

3.2 高速前飛狀態(tài)槳葉氣動(dòng)彈性載荷及流場(chǎng)分析

采用CFD/CSD松耦合方法對(duì)高速飛行狀態(tài)下的彈性SA349/2旋翼槳葉進(jìn)行了載荷預(yù)測(cè)，并與文獻(xiàn)[23]中的試驗(yàn)值進(jìn)行對(duì)比。飛行狀態(tài)1：前進(jìn)比為0.378，CT/σ=0.071。圖11為采用不同方法計(jì)算的在0.97R位置槳葉剖面壓強(qiáng)分布和試驗(yàn)值的對(duì)比?？梢钥闯?，相比于傳統(tǒng)剛性旋翼CFD方法，采用計(jì)入槳葉彈性變形的CFD/CSD耦合方法可以更加精確地捕捉槳葉表面壓力峰值以及激波位置和強(qiáng)度，其中在150°方位角，由于流場(chǎng)計(jì)算中計(jì)入了旋翼在氣動(dòng)載荷作用下的彈性變形影響，槳葉表面未產(chǎn)生氣流分離和激波，槳葉剖面壓強(qiáng)分布得到緩和，這與試驗(yàn)值相一致，而在該方位角下剛性CFD計(jì)算結(jié)果仍有激波出現(xiàn)。

圖11 飛行狀態(tài)1下SA349/2旋翼槳葉表面的壓強(qiáng)系數(shù)分布Fig.11 Distribution of pressure coefficient on blade surface of SA349/2 rotor in Case 1

圖12為計(jì)算獲得的SA349/2旋翼不同槳葉剖面的等效法向力系數(shù)Cn與試驗(yàn)值、剛性CFD以及自由尾跡方法計(jì)算結(jié)果的對(duì)比。圖12(c)中4個(gè)小圖是對(duì)應(yīng)大圖不同位置的剖面壓強(qiáng)系數(shù)?？梢姡捎肅FD/CSD耦合方法的旋翼氣動(dòng)載荷預(yù)測(cè)總體更貼合試驗(yàn)數(shù)據(jù)，由于引入CSD槳葉變形模塊，扭轉(zhuǎn)力矩對(duì)槳葉產(chǎn)生了一定的彈性變形，改變了槳葉剖面的氣動(dòng)分布，與實(shí)際情況相符，升力的幅值和相位相對(duì)于剛性CFD與試驗(yàn)值更加吻合，說明本文旋翼CFD/CSD耦合方法可以有效地預(yù)測(cè)直升機(jī)旋翼氣動(dòng)載荷，為結(jié)構(gòu)載荷的預(yù)測(cè)提供一個(gè)較準(zhǔn)確的外部載荷輸入。

圖13給出了旋翼槳葉剖面法向力系數(shù)在槳盤平面內(nèi)的分布情況。可以看出，在槳葉前行側(cè)，由于槳葉剖面相對(duì)來流速度較大，法向力系數(shù)相對(duì)較大的區(qū)域集中在槳尖區(qū)域，而在270° 后行槳葉一側(cè)，靠近槳葉旋轉(zhuǎn)中心處轉(zhuǎn)速小于μΩR的一段槳葉上，槳葉剖面相對(duì)氣流周向分量自翼型后緣吹向前緣，從而出現(xiàn)“反流區(qū)”，導(dǎo)致法向力出現(xiàn)負(fù)值。

圖14給出了CFD/CSD耦合方法計(jì)算得到的不同方位角槳葉上表面壓強(qiáng)分布，p為無量綱的槳葉表面壓強(qiáng)大小?？梢钥闯觯谛砬靶袀?cè)，由于槳葉的槳尖區(qū)域處于跨聲速區(qū)，該區(qū)域?qū)⒂屑げóa(chǎn)生，從而出現(xiàn)氣流分離現(xiàn)象，導(dǎo)致壓強(qiáng)梯度沿槳葉弦向較大，而后行槳葉槳尖區(qū)域表面展向流動(dòng)比重較大，壓強(qiáng)梯度沿槳葉弦向得到緩和。

圖15給出了CFD/CSD耦合方法計(jì)算得到的旋翼槳葉在各個(gè)方位角處的表面流線及90° 方位角時(shí)槳尖附近表面流動(dòng)細(xì)節(jié)。可以看出，旋翼后行槳葉沿展向的三維流動(dòng)現(xiàn)象非常明顯，而在槳葉前行側(cè)，槳葉表面展向流動(dòng)明顯弱于旋翼后行側(cè)。由于計(jì)算狀態(tài)的前進(jìn)比較大，屬于高速前飛狀態(tài)，在30°～90° 方位角處，槳尖附近會(huì)有氣流分離現(xiàn)象出現(xiàn)，前行槳葉處于跨聲速區(qū)，因而在該區(qū)域會(huì)有激波產(chǎn)生，這與圖11計(jì)算的槳葉剖面壓強(qiáng)分布規(guī)律相對(duì)應(yīng)。

圖12 飛行狀態(tài)1下槳葉不同剖面法向力系數(shù)分布Fig.12 Distribution of normal force coefficient of blade on different cross-sections in Case 1

圖13 法向力系數(shù)沿槳盤分布Fig.13 Distribution of normal force coefficient in the plane disk

圖14 槳葉不同方位角上表面壓強(qiáng)分布 Fig.14 Distribution of surface pressure of blade at different azimuth angles

圖16給出了前進(jìn)比為0.378時(shí)CFD/CSD耦合方法計(jì)算得到SA349/2旋翼槳葉展向不同剖面的非定常揮舞彎矩系數(shù)CFM與試驗(yàn)值及自由尾跡方法計(jì)算結(jié)果的對(duì)比。可以看出，揮舞彎矩相位與幅值捕捉較好，在高速前飛狀態(tài)，槳葉的揮舞結(jié)構(gòu)載荷在90°～220° 方位角內(nèi)較小，但是在30° 附近和270° 附近載荷峰值得到了加強(qiáng)。

圖17給出了非定常擺振彎矩系數(shù)CIM的計(jì)算結(jié)果對(duì)比?？梢钥闯?，擺振彎矩幅值基本一致，但是在槳葉后行側(cè)相位捕捉有一定偏差，可能與缺少擺振阻尼器的模擬有關(guān)。

圖18給出了扭轉(zhuǎn)力矩系數(shù)CTM的計(jì)算結(jié)果對(duì)比。扭轉(zhuǎn)力矩與試驗(yàn)值吻合較好，高速狀態(tài)下，扭轉(zhuǎn)載荷高頻響應(yīng)成分微弱，整體趨勢(shì)較為平緩，90° 方位角附近出現(xiàn)抬頭力矩峰值，而在180° 方位角附近出現(xiàn)低頭力矩峰值。

圖15 槳葉在不同方位角處的表面流線Fig.15 Surface streamlines of blade at different azimuth angles

圖16 飛行狀態(tài)1下?lián)]舞彎矩系數(shù)對(duì)比Fig.16 Comparison of flap bending moment coefficient in Case 1

圖17 飛行狀態(tài)1下擺振彎矩系數(shù)對(duì)比Fig.17 Comparison of lag bending moment coefficient in Case 1

圖18 飛行狀態(tài)1下扭轉(zhuǎn)力矩系數(shù)對(duì)比Fig.18 Comparison of torsion moment coefficient in Case 1

3.3 低速前飛狀態(tài)槳葉氣動(dòng)彈性載荷分析驗(yàn)證

飛行狀態(tài)2：前進(jìn)比為0.140，CT/σ=0.065。圖19為CFD/CSD耦合方法計(jì)算的旋翼槳葉氣動(dòng)載荷，并與文獻(xiàn)[23]中的試驗(yàn)值、剛性CFD方法以及自由尾跡方法結(jié)果進(jìn)行對(duì)比?？梢?，CFD/CSD耦合方法和剛性CFD方法在一定程度上都能較好地反映旋翼非定常氣動(dòng)載荷變化，而自由尾跡方法在低速狀態(tài)下計(jì)算結(jié)果和試驗(yàn)結(jié)果在幅值和相位上都存在一定差距。由于前進(jìn)比較低，該飛行狀態(tài)下，旋翼槳-渦干擾嚴(yán)重，氣動(dòng)載荷的非定常效應(yīng)較飛行狀態(tài)1的高速前飛更加明顯，在圖19(b)和圖19(c)中，測(cè)試剖面靠近槳尖，90° 和270° 處氣動(dòng)載荷振蕩劇烈，相對(duì)于剛性CFD方法，本文發(fā)展的旋翼CFD/CSD耦合方法在預(yù)測(cè)槳-渦干擾時(shí)的氣動(dòng)載荷幅值方面更加精確。

圖19 飛行狀態(tài)2下槳葉不同剖面法向力系數(shù)分布Fig.19 Distribution of normal force coefficient of blade on different cross-sections in Case 2

圖20 飛行狀態(tài)2下?lián)]舞彎矩系數(shù)對(duì)比Fig.20 Comparison of flap bending moment coefficient in Case 2

圖20為飛行狀態(tài)2下的槳葉剖面揮舞彎矩系數(shù)CFD/CSD耦合方法與試驗(yàn)結(jié)果及自由尾跡方法對(duì)比。CSD/CFD耦合方法計(jì)算得到的揮舞彎矩整體分布與試驗(yàn)值的對(duì)比呈現(xiàn)一致的變化趨勢(shì)，從圖20中可以看出，在270° 方位角處出現(xiàn)了一個(gè)明顯的載荷峰值，這主要由于槳-渦干擾造成的沖量作用使槳尖處氣動(dòng)力載荷達(dá)峰值而引起的結(jié)果，而自由尾跡方法由于低速狀態(tài)下非定常氣動(dòng)力計(jì)算存在一定誤差，導(dǎo)致計(jì)算得到的非定常揮舞彎矩與試驗(yàn)值相比誤差較大。

4 結(jié) 論

基于Navier-Stokes方程、彈性槳葉運(yùn)動(dòng)嵌套網(wǎng)格方法、網(wǎng)格變形方法和Hamilton變分原理，建立了一套用于旋翼非定常氣動(dòng)特性和氣動(dòng)彈性載荷分析的旋翼CFD/CSD松耦合方法，并進(jìn)行了較系統(tǒng)的算例驗(yàn)證。

1) 流固耦合策略及槳葉彈性網(wǎng)格變形方法能夠有效地進(jìn)行CFD與CSD模塊的氣動(dòng)力和彈性變形之間的信息交換。

2) 相比于剛性旋翼CFD方法，旋翼CFD/CSD耦合方法在高速前飛狀態(tài)計(jì)算得到的槳尖區(qū)域槳葉剖面的壓強(qiáng)分布相比于剛性槳葉明顯改善，可以更準(zhǔn)確地捕捉槳葉剖面壓力峰值以及激波位置和強(qiáng)度，與實(shí)際情況相符；在低速前飛狀態(tài)，可以更加準(zhǔn)確地模擬由于槳-渦干擾引起的氣動(dòng)載荷非定常效應(yīng)。

3) 相比于傳統(tǒng)氣彈耦合方法中的自由尾跡模型，旋翼CFD/CSD耦合方法可以更加準(zhǔn)確地捕捉揮舞彎矩、擺振彎矩和扭轉(zhuǎn)力矩的相位和幅值，結(jié)果反映氣動(dòng)載荷的精度對(duì)結(jié)構(gòu)載荷的計(jì)算有重要意義。

相比于剛性旋翼CFD方法和自由尾跡方法，CFD/CSD耦合方法普遍適用于直升機(jī)高速前飛狀態(tài)和低速前飛狀態(tài)的氣動(dòng)彈性載荷預(yù)估，為旋翼氣動(dòng)彈性綜合分析中的減振優(yōu)化、強(qiáng)度校核及氣動(dòng)彈性穩(wěn)定性研究工作提供基礎(chǔ)。

[1] CONLISK A T. Modern helicopter rotor aerodynamics[J]. Progress in Aerospace Science, 2001, 37(5): 419-476.

[2] DATTA A. Fundamental understanding, prediction and validation of rotor vibratory loads in steady-level flight[D]. Park: University of Maryland, 2004.

[3] DATTA A, CHOPRA I. Prediction of the UH-60A main rotor structural loads using computational fluid dynamics/comprehensive analysis coupling[J]. Journal of the American Helicopter Society, 2008, 53(4): 351-365.

[4] TUNG C, CARADONNA F X, JOHNSON W. The prediction of transonic flows on an advancing rotor[J]. Journal of the American Helicopter Society, 1986, 32(7): 4-9.

[5] POTSDAM M, YEO H, JOHNSON W. Rotor airloads prediction using loose aerodynamic/structural coupling[J]. Journal of Aircraft, 2006, 43(3): 732-742.

[6] ZHAO J G, HE C J. Rotor blade structural loads analysis using coupled CSD/CFD/VVPM [C]//Proceedings of the 69th Annual Forum of the American Helicopter Society, Phoenix, 2013： 1452-1468.

[7] LIM J W. Consideration of structural constraints in passive rotor blade design for improved performance[J]. Aeronautical Journal, 2015, 119(1222): 1513-1539.

[8] 王海. 計(jì)入槳葉結(jié)構(gòu)彈性的新型槳尖旋翼流場(chǎng)數(shù)值模擬研究[D]. 南京: 南京航空航天大學(xué), 2010.

WANG H. Numerical simulation for the flowfield of new-tip rotors with effect of blade elasticity[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2010 (in Chinese).

[9] 王俊毅, 招啟軍, 馬礫. 直升機(jī)旋翼槳-渦干擾狀態(tài)非定常氣彈載荷高精度預(yù)估[J]. 航空動(dòng)力學(xué)報(bào), 2015, 30(5): 1267-1274.

WANG J Y, ZHAO Q J, MA L. High-precision prediction on unsteady aeroelastic loads of helicopter rotors in BVI condition[J]. Journal of Aerospace Power, 2015, 30(5): 1267-1274 (in Chinese).

[10] YUAN K, FRIEDMANN P. Aeroelasticity and structural optimization of composite helicopter rotor blades with swept tips:NASA CR-4665[R]. Washington, D.C.: NASA, 1995.

[11] YAMAKAWA H, OHNISHI T. Dynamic response analysis with many degrees of freedom using step-by-step transfer matrix method[J]. Transactions of the Japan Society of Mechanical Engineers C, 1982, 48(429): 672-681.

[12] DATTA A, SITARMAN J, CHOPRA I, et al. CFD/CSD prediction of rotor vibratory loads in high-speed flight[J]. Journal of Aircraft, 2006, 43(6): 1698-1709.

[13] HSU K, LEE S L. A numerical technique for two-dimensional grid generation with grid control at all of the boundaries[J]. Journal of Computational Physics, 1991, 96(2): 451-469.

[14] 趙國(guó)慶, 招啟軍, 吳琪. 旋翼非定常氣動(dòng)特性CFD模擬的通用運(yùn)動(dòng)嵌套網(wǎng)格方法[J]. 航空動(dòng)力學(xué)報(bào), 2015, 30(3): 546-554.

ZHAO G Q, ZHAO Q J, WU Q. A universal moving-embedded grid method for CFD simulation of unsteady aerodynamic characteristics of rotor[J]. Journal of Aerospace Power, 2015, 30(3): 546-554 (in Chinese).

[15] 趙國(guó)慶, 招啟軍, 王清. 旋翼翼型非定常動(dòng)態(tài)失速特性的CFD 模擬及參數(shù)分析[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2015, 33(1): 72-81.

ZHAO G Q, ZHAO Q J, WANG Q. Simulations and parametric analyses on unsteady dynamic stall characteristics of rotor airfoil based on CFD method[J]. Acta Aerodynamica Sinica, 2015, 33(1): 72-81 (in Chinese).

[16] BLAZEK J. Computational fluid dynamics: Principles and applications, second edition[M]. Amsterdam: Elsevier, 2005: 16-18.

[17] SPALART P R, ALLMARAS S R. A one-equation turbulence model for aerodynamic flows: AIAA-1991-0439[R]. Reston: AIAA, 1991.

[18] LUO H, BAUM J D, LOHNER R. A fast, matrix-free implicit method for compressible flows on unstructured grids[J]. Journal of Computational Physics, 1998, 146(2): 664-690.

[19] ROE P L. Approximate Riemann solvers, parameter vectors and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372.

[20] JOHNSON W. Milestones in rotorcraft aeromechanics: NASA/TP-2011-215971[R]. Washington, D.C.: NASA, 2011.

[21] POTSDAM M, YEO H, JOHNSON W. Rotor airloads prediction using loose aerodynamic/structural coupling[J]. Journal of Aircraft, 2006, 43(5): 732-742.

[22] CARADONNA F X, LAUB G H, TUNG C. An experimental investigation of the parallel blade-vortex interaction: NASA TM-86005[R]. Washington, D.C.: NASA, 1984.

[23] 趙景根, 徐國(guó)華, 招啟軍. 基于自由尾跡分析的直升機(jī)旋翼下洗流場(chǎng)計(jì)算方法[J]. 兵工學(xué)報(bào), 2006, 27(1): 63-68.

ZHAO J G, XU G H, ZHAO Q J. A calculating method of helicopter rotor downwash flowfield based on free wake analysis[J]. Acta Armamentarii, 2006, 27(1): 63-68 (in Chinese).

[24] HEFFERMAN R M, GAUBERT M. Structural and aerodynamic loads and performance measurements of an SA349/2 helicopter with an advanced geometry rotor: NASA TM-88370[R]. Washington, D.C.: NASA,1986.

(責(zé)任編輯： 鮑亞平)

Computation analyses of aeroelastic loads of rotor based onCFD/CSD coupling method

MALi,ZHAOQijun*,ZHAOMengmeng,WANGBo

NationalKeyLaboratoryofScienceandTechnologyonRotorcraftAeromechanics,NanjingUniversityofAeronauticsandAstronautics,Nanjing210016,China

In order to improve the prediction accuracy of the unsteady aeroelastic load of the rotor in forward flight of the helicopter, the rotor CFD module is introduced into the comprehensive analysis method for rotor aeroelastic, and then an analytical model and corresponding code based on CFD/CSD loose coupling method are established. In order to solve the deformation problem of blade body-fitted grid caused by aeroelastic deflections such as flapping and torsional motion in rotor CFD/CSD coupling process, the blade grids are deformed using algebraic method, and blade motion deformation and rotor aerodynamic force information are transferred through the interface between fluid and solid. In the analysis method for rotor flowfield, the Navier-Stokes equations coupled with S-A turbulence model are used as governing equations, and then the rotor moving-embedded grids are constructed. The Roe scheme is employed in spatial discretization, and a dual-time algorithm is adopted for temporal discretization. In the analysis of the rotor structure, considering the rotor trim, the nonlinear equations for motion of elastic rotor are solved based on Hamilton variational principle and 20 degree of freedom Timoshenko beam model. The CSD and CFD methods are validated . Based on aforementioned work, the aeroelastic loads on the SA349/2 rotor blade, mainly including unsteady aerodynamic loads and moment in flapping, lagging and torsion direction, are calculated and compared with flight test data. Comparisons of the CFD/CSD calculated results with flight test data demonstrate that the CFD/CSD coupling method has high accuracy in calculating the rotor unsteady aerodynamic load and structural load, and can capture the peak of pressure and shock position accurately. The method proposed could be effectively applied to prediction analysis of the aeroelastic loads of the rotor.

rotor; aeroelastic load; CFD/CSD coupling method; Navier-Stokes equations; grid deformation technology; forward flight

2016-09-08;Revised2016-10-06;Accepted2016-11-09;Publishedonline2016-11-280922

URL：www.cnki.net/kcms/detail/11.1929.V.20161128.0922.002.html

s：NationalNaturalScienceFoundationofChina(11572156);AProjectFundedbythePriorityAcademicProgramDevelopmentofJiangsuHigherEducationInstitutions

2016-09-08；退修日期2016-10-06；錄用日期2016-11-09； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-11-280922

www.cnki.net/kcms/detail/11.1929.V.20161128.0922.002.html

國(guó)家自然科學(xué)基金 (11572156)； 江蘇省高校優(yōu)勢(shì)學(xué)科建設(shè)工程資助項(xiàng)目

*

.E-mailzhaoqijun@nuaa.edu.cn

馬礫， 招啟軍， 趙蒙蒙, 等． 基于CFD/CSD耦合方法的旋翼氣動(dòng)彈性載荷計(jì)算分析J． 航空學(xué)報(bào),2017,38(6):120762.MAL，ZHAOQJ,ZHAOMM,etal.ComputationanalysesofaeroelasticloadsofrotorbasedonCFD/CSDcouplingmethodJ.ActaAeronauticaetAstronauticaSinica,2017,38(6):120762.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0292

V211.47

A

1000-6893(2017)06-120762-14

*Correspondingauthor.E-mailzhaoqijun@nuaa.edu.cn