基于CFD／CSD松耦合的直升機(jī)配平分析方法

余瑾，王松，劉勇，楊衛(wèi)東

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

在直升機(jī)旋翼動(dòng)力學(xué)分析中，配平計(jì)算不僅是飛行力學(xué)和操縱品質(zhì)分析的基礎(chǔ)，旋翼的配平解對(duì)旋翼的氣彈穩(wěn)定性和氣彈載荷分析也有著至關(guān)重要的影響［1］。直升機(jī)旋翼氣彈分析是一個(gè)具有挑戰(zhàn)性的多學(xué)科問(wèn)題，精確的槳葉氣動(dòng)建模涉及到三維非定常流場(chǎng)、跨聲速流、動(dòng)態(tài)失速、槳葉剛體運(yùn)動(dòng)和彈性變形等，這就必須與CSD分析相結(jié)合進(jìn)行綜合分析。在完全耦合的氣彈分析中，由于槳葉的剛體運(yùn)動(dòng)、彈性變形、氣動(dòng)載荷和旋翼配平狀態(tài)相互影響，槳葉的氣動(dòng)力計(jì)算和結(jié)構(gòu)動(dòng)力學(xué)分析也是相互依賴(lài)、相互作用的。為了解決這一復(fù)雜問(wèn)題，旋翼綜合分析代碼往往采用基于升力線理論和查表法等簡(jiǎn)化的空氣動(dòng)力學(xué)模型，然而使用這些快速方法進(jìn)行載荷預(yù)測(cè)，往往會(huì)顯示出明顯的不足，難以描述細(xì)致的流場(chǎng)細(xì)節(jié)，無(wú)法滿足氣彈載荷分析的精度要求［2］。本文的研究目的就是將CFD求解器耦合進(jìn)氣彈綜合分析代碼中，以提供高保真的氣動(dòng)模型，克服氣彈綜合分析中氣動(dòng)模型精度不足的缺陷，而氣彈分析模塊則依然進(jìn)行結(jié)構(gòu)動(dòng)力學(xué)和配平分析。

CFD模塊與氣彈綜合分析模塊的耦合有2種思路：松耦合和緊耦合。在松耦合策略中，CFD模塊和CSD模塊分別在時(shí)域內(nèi)推進(jìn)，每隔一段時(shí)間（如旋翼旋轉(zhuǎn)一圈）交互一次數(shù)據(jù)；在緊耦合策略中，CFD模塊和CSD模塊集成在一起，每個(gè)時(shí)間步長(zhǎng)上都要進(jìn)行數(shù)據(jù)交換。盡管緊耦合策略更為嚴(yán)密，但必須要在CFD／CSD兩個(gè)模塊間的時(shí)間精度控制和數(shù)據(jù)交換效率上付出更高的代價(jià)。另一方面，配平分析至少包含2層循環(huán)：外層配平迭代循環(huán)和內(nèi)層氣彈分析響應(yīng)求解迭代循環(huán)。采用緊耦合策略進(jìn)行配平計(jì)算意味著將CFD求解器放進(jìn)了內(nèi)層的響應(yīng)求解迭代運(yùn)算，為此將付出更大的計(jì)算代價(jià)。Johnson等［2］指出，就當(dāng)前的計(jì)算能力而言，是禁止將CFD技術(shù)直接引入到配平循環(huán)內(nèi)的。因此，將緊耦合方法用于配平分析是不妥的。而松耦合策略允許采用模塊化方法，各學(xué)科模塊根據(jù)需要處理其時(shí)間精度，配平解則是綜合分析的自然結(jié)果。

Tung等［3］于1986年提出了CFD／CSD之間的松耦合方法，其中CFD求解器運(yùn)用小擾動(dòng)方程，CSD求解器基于CAMRAD平臺(tái)，通過(guò)升力線法修正CFD計(jì)算的氣動(dòng)力，在該模型中，通過(guò)CSD模塊更新入流角從而修正物面無(wú)穿透條件代替驅(qū)動(dòng)網(wǎng)格變形，該方法中CFD計(jì)算精度不足，故耦合存在收斂的問(wèn)題。CFD／CSD松耦合的概念提出之后，被廣泛認(rèn)可和應(yīng)用。2001年，Costes等［4］基于松耦合方法，運(yùn)用WAVES的CFD求解器和HOST的CSD求解器，進(jìn)行彈性旋翼槳葉前飛狀態(tài)的配平計(jì)算，研究表明該方法對(duì)變距力矩和扭轉(zhuǎn)載荷的預(yù)估具有較大的提高。2006年，Johnson等［2］基于CFD／CSD松耦合的方法，CSD求解器使用CAMRAD平臺(tái)，CFD模塊則基于OVERFLOW 求解器，以UH-60A黑鷹直升機(jī)為算例，對(duì)旋翼不同飛行狀態(tài)下進(jìn)行振動(dòng)載荷的預(yù)估，取得較好的成果，并成功捕捉到槳尖渦現(xiàn)象，該方法具有很好的穩(wěn)定性、收斂性及魯棒性，并且對(duì)法向力和變距力矩的幅值和相位的預(yù)測(cè)與飛行測(cè)試數(shù)據(jù)具有很好的吻合，且沒(méi)有使用修正物面無(wú)穿透條件代替驅(qū)動(dòng)網(wǎng)格變形，而是通過(guò)OVERFLOW的動(dòng)網(wǎng)格和“網(wǎng)格挖洞”技術(shù)實(shí)現(xiàn)。2006年，Chopra等［5］建立基于UMARC的CSD綜合分析平臺(tái)和基于雷諾平均的Navier-Stokes方程的CFD模型，采用松耦合的方法對(duì)UH-60A直升機(jī)旋翼高速前飛狀態(tài)的振動(dòng)載荷進(jìn)行了預(yù)估，分析表明二階精確非線性梁理論對(duì)高速前飛狀態(tài)的結(jié)構(gòu)動(dòng)載荷具有較好的預(yù)估能力，并且準(zhǔn)確預(yù)估到1～3階扭轉(zhuǎn)諧波載荷。國(guó)內(nèi)基于CFD／CSD耦合的方法對(duì)飛行器研究相對(duì)較晚。2010年，王海［6］采用以嵌套網(wǎng)格為基礎(chǔ)的CFD模塊，提出直升機(jī)旋翼CFD／CSD松耦合分析方法。2013年，王俊毅、招啟軍［7］發(fā)展了旋翼穩(wěn)態(tài)前飛狀態(tài)的CFD／CSD耦合分析方法，研究表明，改變剖面的扭轉(zhuǎn)剛度旋翼氣動(dòng)載荷變化較大。2014年，肖宇、徐國(guó)華等［8］基于非慣性系下的Navier-Stokes方程，改進(jìn)了懸停狀態(tài)下的CFD／CSD耦合研究方法，UH-60A直升機(jī)的算例表明耦合策略具有更好的載荷預(yù)估能力。2015年，李建東［9］發(fā)展了二維翼型的CFD／CSD緊耦合策略。2016年，黃道博［10］基于CFD／CSD耦合方法準(zhǔn)確預(yù)估了直升機(jī)穩(wěn)態(tài)前飛狀態(tài)旋翼氣動(dòng)載荷和振動(dòng)載荷。CFD模塊以FLUENT商用軟件為平臺(tái)，CSD模塊基于中等變形梁理論，通過(guò)UDF進(jìn)行槳葉運(yùn)動(dòng)和網(wǎng)格變形信息的交換。

在本文采用的CFD／CSD松耦合策略中，CFD模塊和CSD模塊分別在時(shí)域內(nèi)推進(jìn)，以槳葉彈性軸和變距軸線為媒介，通過(guò)線性插值方法交換氣動(dòng)載荷和響應(yīng)數(shù)據(jù)，旋翼每旋轉(zhuǎn)一圈交互一次。在耦合迭代過(guò)程中，隨著氣動(dòng)載荷的改變，改變旋翼的操縱量和飛行姿態(tài)以滿足配平方程組。本文構(gòu)造了內(nèi)外兩層迭代：內(nèi)層迭代為配平求解過(guò)程，隨著氣動(dòng)載荷的改變，以牛頓迭代法改變旋翼的操縱量和飛行姿態(tài)，滿足配平方程組；外層循環(huán)即以CFD模塊計(jì)算的氣動(dòng)力來(lái)修正配平計(jì)算中氣彈分析的氣動(dòng)力輸入，直到配平量和CFD氣動(dòng)力在外層迭代中不再變化，達(dá)到同時(shí)收斂，即得到了耦合配平解。

1 旋翼流場(chǎng)CFD建模

直升機(jī)旋翼CFD技術(shù)一直以來(lái)是計(jì)算流體力學(xué)領(lǐng)域的熱點(diǎn)和較為困難的問(wèn)題之一，這與直升機(jī)旋翼的復(fù)雜工作環(huán)境及旋翼特有的運(yùn)動(dòng)和操縱方式有關(guān)。網(wǎng)格是CFD計(jì)算的基礎(chǔ)，就直升機(jī)旋翼而言，計(jì)算空間生成單一的計(jì)算網(wǎng)格是相當(dāng)困難的，這會(huì)使得網(wǎng)格局部存在嚴(yán)重的扭曲或畸變，故目前多采用多個(gè)分區(qū)對(duì)接及重疊嵌套網(wǎng)格避免這一問(wèn)題。本文基于OversetGrid網(wǎng)格前處理軟件繪制直升機(jī)旋翼嵌套網(wǎng)格，基于OverCFD軟件，綜合考慮旋翼槳葉運(yùn)動(dòng)與變形、網(wǎng)格自適應(yīng)、“網(wǎng)格挖洞”與插值、湍流模型等計(jì)算旋翼流場(chǎng)信息。

1.1 適合CFD使用的控制方程

流體運(yùn)動(dòng)的基本規(guī)律需滿足物理學(xué)三大定律，即質(zhì)量守恒定律、動(dòng)量守恒定律及能量守恒定律。

一般局部坐標(biāo)系下的Navier-Stokes方程具有如下形式［11］：

式中：ξ、η、ζ分別為微元體的一般局部坐標(biāo)系3個(gè)方向；守恒變量為其中：ρ為流體密度；u、v、w分別為笛卡兒坐標(biāo)xy-z三個(gè)方向上的絕對(duì)速度分量；e0和ei分別為總能量和內(nèi)能。

無(wú)黏通量為

式中：p為壓強(qiáng)。

黏性項(xiàng)為

其中：τxx和bx分別為黏性力的切應(yīng)力項(xiàng)和正應(yīng)力項(xiàng)，下標(biāo)標(biāo)明各自不同方向上的分量。

J為笛卡兒坐標(biāo)與一般坐標(biāo)系間轉(zhuǎn)換的Jacobian矩陣，即

1.2 旋翼槳葉貼體網(wǎng)格

圖1為SA349／2旋翼槳葉貼體網(wǎng)格，其旋翼采用3片OA209翼型的矩形高速槳葉，并具有一定的負(fù)扭轉(zhuǎn)。本文對(duì)其槳根位置進(jìn)行簡(jiǎn)化處理，在離散具有負(fù)扭轉(zhuǎn)區(qū)域操作如下：提取三角面元網(wǎng)格的輪廓線得到槳葉葉素（剖面）線網(wǎng)格，對(duì)剖面線網(wǎng)格重布后，基于OversetGrid軟件平移旋轉(zhuǎn)操作，在線性負(fù)扭轉(zhuǎn)槳葉段生成結(jié)構(gòu)化面網(wǎng)格，若是非線性負(fù)扭轉(zhuǎn)，取多段進(jìn)行以上操作，對(duì)旋翼槳葉結(jié)構(gòu)面網(wǎng)格進(jìn)行重布，在槳尖和槳根位置加密。因此，槳尖和槳根位置應(yīng)具有更多的重疊區(qū)域。

圖1 SA349／2旋翼槳葉貼體網(wǎng)格Fig.1 Body-fitted grid of SA349／2 rotor blade

主體網(wǎng)格生成方法步驟如下：

步驟1 幾何建模軟件CATIA生成旋翼槳葉的幾何模型，進(jìn)而使用ICEM 或OversetGrid離散得到三角元槳葉表面網(wǎng)格。

步驟2 提取三角元網(wǎng)格的輪廓線，生成槳葉線網(wǎng)格，并重布線網(wǎng)格（重布線網(wǎng)格時(shí)，槳根槳尖的展向周向都需加密，以便生成合適的貼體帽子網(wǎng)格）。

步驟3 通過(guò)線網(wǎng)格生成槳葉表面結(jié)構(gòu)化面網(wǎng)格，并在2個(gè)端面生成重疊嵌套的端面網(wǎng)格，稱(chēng)之為“帽子”網(wǎng)格；在表面網(wǎng)格的基礎(chǔ)上采用Hyp.tangent體網(wǎng)格推進(jìn)方法生成近體網(wǎng)格。

兩側(cè)帽子網(wǎng)格方法步驟如下：

步驟1 抽取槳葉主體面網(wǎng)格兩端線網(wǎng)格，從而生成2個(gè)端面網(wǎng)格。

步驟2 提取主體面網(wǎng)格前緣后緣網(wǎng)格（提取的前、后緣網(wǎng)格節(jié)點(diǎn)數(shù)相同）。

步驟3 將步驟1、步驟2提取出的網(wǎng)格進(jìn)行連接，并提取新生成的端面網(wǎng)格上下的線網(wǎng)格。

步驟4 指定主面網(wǎng)格為翼型結(jié)構(gòu)化網(wǎng)格的參考面，以該參考面生成緊貼主面網(wǎng)格的新的面網(wǎng)格，并將其與新生成的端面面網(wǎng)格連接，從而形成帽子網(wǎng)格。

表1為網(wǎng)格規(guī)模，其貼體網(wǎng)格系統(tǒng)總數(shù)為1 150萬(wàn)。

表1 SA349／2旋翼槳葉貼體網(wǎng)格規(guī)模統(tǒng)計(jì)Tab le 1 Num bers of body-fitted grid of SA349／2 rotor blade

1.3 旋翼槳葉背景網(wǎng)格

對(duì)于前飛狀態(tài)，采用笛卡兒自動(dòng)背景網(wǎng)格結(jié)合網(wǎng)格自適應(yīng)，應(yīng)用于動(dòng)態(tài)計(jì)算，采用網(wǎng)格自適應(yīng)功能時(shí)，網(wǎng)格數(shù)量會(huì)隨著計(jì)算過(guò)程增加，故流場(chǎng)和壓力系數(shù)曲線會(huì)更加光順，對(duì)槳尖渦具有很好的捕捉能力，結(jié)果更加準(zhǔn)確。

笛卡兒自動(dòng)背景網(wǎng)格系統(tǒng)如圖2所示，其部分剖面圖如圖3和圖4所示。此背景網(wǎng)格共有9級(jí)，以應(yīng)用于計(jì)算懸停狀態(tài)。第1級(jí)背景網(wǎng)格（L1）由均勻網(wǎng)格間隔（Δ＝10%ctip）構(gòu)成，其他“磚網(wǎng)格”（L2，L3，…）通過(guò)逐層增加到L1網(wǎng)格外圍，迅速擴(kuò)展到計(jì)算外場(chǎng)。計(jì)算外場(chǎng)尺寸為17倍旋翼半徑，每一級(jí)別背景網(wǎng)格是上一級(jí)背景網(wǎng)格間隔的2倍，其中粗網(wǎng)格具有Δ2＝2Δ1，Δ3＝4Δ1，Δ4＝8Δ1，…；在網(wǎng)格自適應(yīng)中，網(wǎng)格間隔具有類(lèi)似關(guān)系式，Δ-1＝Δ1／2，Δ-2＝Δ1／4，…。對(duì)于L1網(wǎng)格，其尺寸和間隔由手動(dòng)定義，旋翼的展向擴(kuò)展為1.25R（R為旋翼半徑），旋翼上方和下方擴(kuò)展為0.3R，其第1級(jí)背景網(wǎng)格（L1）間隔仍取Δ＝10%ctip。

圖2 SA349／2旋翼流場(chǎng)笛卡兒自動(dòng)背景網(wǎng)格Fig.2 Cartesian automatic background grid of SA349／2 rotor flow field

圖3 適合前飛狀態(tài)求解的笛卡兒自動(dòng)背景網(wǎng)格剖面Fig.3 Cartesian automatic background grid section for solving in forward flight status

圖4 旋翼前飛狀態(tài)笛卡兒背景網(wǎng)格自適應(yīng)Fig.4 Cartesian adaptive background grid of rotor for forward flight

笛卡兒第9級(jí)背景網(wǎng)格之間的間隔關(guān)系可表示為

1.4 旋翼網(wǎng)格運(yùn)動(dòng)與變形

對(duì)于穩(wěn)態(tài)飛行狀態(tài)，旋翼槳葉的運(yùn)動(dòng)是在配平位置處的穩(wěn)態(tài)周期響應(yīng)，包括剛體運(yùn)動(dòng)和彈性變形。本文在對(duì)旋翼網(wǎng)格運(yùn)動(dòng)中分別施加旋翼槳葉的剛體運(yùn)動(dòng)和彈性變形，兩者疊加實(shí)現(xiàn)槳葉真實(shí)運(yùn)動(dòng)。旋翼剛體運(yùn)動(dòng)通過(guò)GMP定義，彈性變形通過(guò)Motion（即旋翼旋轉(zhuǎn)一周彈性變形關(guān)于方位角、槳葉展向位置的函數(shù)）文件來(lái)定義。驅(qū)動(dòng)網(wǎng)格運(yùn)動(dòng)，更新網(wǎng)格位置，挖洞插值，循環(huán)直到計(jì)算模擬收斂，流程如圖5所示。

旋翼槳葉剛體運(yùn)動(dòng)是指在工作過(guò)程中沒(méi)有發(fā)生彈性變形。旋翼在繞旋翼軸旋轉(zhuǎn)運(yùn)動(dòng)的同時(shí)，具有的揮舞、擺振、扭轉(zhuǎn)運(yùn)動(dòng)可以設(shè)定為槳葉整體繞不同軸的單自由度轉(zhuǎn)動(dòng)，即繞揮舞鉸的揮舞運(yùn)動(dòng)βb、繞擺振鉸的擺振運(yùn)動(dòng)ζb及繞變距軸線的扭轉(zhuǎn)運(yùn)動(dòng)θb。

為描述旋翼姿態(tài)及槳葉的剛體揮舞、擺振及扭轉(zhuǎn)運(yùn)動(dòng)：

圖5 旋翼網(wǎng)格運(yùn)動(dòng)與變形Fig.5 Motion and deformation of rotor grid

1）槳盤(pán)平面的姿態(tài)角。針對(duì)孤立旋翼情況，旋翼軸的幾何前傾角與配平得到的機(jī)身姿態(tài)角可通過(guò)前處理軟件OversetGrid調(diào)整網(wǎng)格位置實(shí)現(xiàn)，亦可在求解流場(chǎng)時(shí)在來(lái)流中定義。

2）單片槳葉組件（包括槳葉表面和2個(gè)端面網(wǎng)格）。網(wǎng)格剛體運(yùn)動(dòng)過(guò)程中，通過(guò)驅(qū)動(dòng)組件運(yùn)動(dòng)實(shí)現(xiàn)網(wǎng)格運(yùn)動(dòng)，每片槳葉定義為一個(gè)組件，實(shí)現(xiàn)槳葉繞軸的轉(zhuǎn)動(dòng)和揮-擺-扭／操縱的運(yùn)動(dòng)，先定義鉸偏置量，后定義揮-擺-扭／操縱的速率。

3）旋翼組件。定義旋轉(zhuǎn)的轉(zhuǎn)向和速率。同一旋翼的槳葉組件共同定義在一個(gè)“父”組件下。旋翼組件和每個(gè)單片槳葉之間存在相對(duì)運(yùn)動(dòng)，即槳葉組件隨“父”組件旋轉(zhuǎn)運(yùn)動(dòng)外，還均有自身的運(yùn)動(dòng)。

每個(gè)組件生成一組重疊網(wǎng)格，根據(jù)組件定義的運(yùn)動(dòng)更新其所對(duì)應(yīng)的網(wǎng)格位置，再重新進(jìn)行網(wǎng)格的挖洞和插值，并計(jì)算求解新的流場(chǎng)。槳葉的彈性變形處理為一維梁的運(yùn)動(dòng)，忽略剖面的翹曲，通過(guò)本文CSD模塊計(jì)算槳葉的彈性軸運(yùn)動(dòng)，提取每個(gè)展向站點(diǎn)處3個(gè)方向的平動(dòng)和轉(zhuǎn)動(dòng)，分別為拉伸位移、擺振位移、揮舞位移、扭轉(zhuǎn)角、揮舞角、擺振角信息，以此表達(dá)槳葉剖面的運(yùn)動(dòng)；并表示成Motion文件，傳遞給CFD模塊。

Motion文件中定義有第一片槳葉變距軸線信息，在CFD計(jì)算時(shí)，通過(guò)與網(wǎng)格信息的比較找到第一片槳葉位置，利用這些旋翼槳葉數(shù)據(jù)表，采用樣條插值將數(shù)據(jù)傳遞給網(wǎng)格剖面，驅(qū)動(dòng)剖面的變形運(yùn)動(dòng)，最終更新槳葉近體網(wǎng)格的位置，實(shí)現(xiàn)槳葉的彈性變形。

圖6為通過(guò)Motion文件施加旋翼?yè)]舞方向的變形，并與未變形槳葉對(duì)比。在網(wǎng)格變形后，仍通過(guò)挖洞、貢獻(xiàn)單元搜索和插值計(jì)算旋翼流場(chǎng)。

圖6 槳葉彈性變形示意圖Fig.6 Schematic of blade elastic deformation

2 旋翼氣彈綜合分析模型和CFD／CSD松耦合策略

本文基于Hodges和Dowell［12］建立的中等變形梁理論，運(yùn)用Ham ilton原理推導(dǎo)旋翼槳葉的運(yùn)動(dòng)方程，以空間有限元方法［13］離散旋翼槳葉，運(yùn)用模態(tài)疊加法和直接數(shù)值積分法建立直升機(jī)旋翼槳葉氣彈綜合分析模型［14］。

引入了CFD氣動(dòng)力修正的CSD模塊計(jì)算流程如圖7所示。

圖7 CSD模塊與CFD氣動(dòng)力耦合計(jì)算流程Fig.7 Flowchart of CSD／CFD aerodynamic force coupling calculation

圖8為旋翼配平及耦合策略?xún)?nèi)外兩層循環(huán)流程，耦合計(jì)算以CSD模塊平臺(tái)作為主控平臺(tái)。啟動(dòng)計(jì)算，氣動(dòng)力借助查表法獲得，計(jì)算旋翼力

對(duì)于第1（i＝1）次配平，運(yùn)用LB非定常氣動(dòng)模型計(jì)算氣動(dòng)力，并采用上一圈CFD計(jì)算的氣動(dòng)力修正查表法。

圖8 旋翼配平及耦合策略?xún)?nèi)外兩層循環(huán)流程Fig.8 Circulation flowchart of inner and outer loop for rotor trim and coupling strategy

式中：上標(biāo)lb表示運(yùn)用LB非定常氣動(dòng)模型計(jì)算得到的旋翼氣動(dòng)力，CFD表示采用CFD技術(shù)得到的旋翼氣動(dòng)力。

對(duì)于第i（i＞1）次配平，使用上一圈CFD計(jì)算的氣動(dòng)力修正LB非定常氣動(dòng)模型。

式（14）也可理解為第i-1次CFD氣動(dòng)載荷結(jié)果無(wú)法滿足配平平衡方程組時(shí)，通過(guò)LB非定常氣動(dòng)模型結(jié)合查表法修正CFD計(jì)算結(jié)果使其滿足配平方程組并求解得到配平解，即

耦合計(jì)算通過(guò)使用OverCFD求解器及MATLAB編程的CSD綜合分析程序進(jìn)行，基于Linux平臺(tái)，采用MPI并行計(jì)算協(xié)議。數(shù)據(jù)交換通過(guò)CFD求解器計(jì)算的剖面氣動(dòng)力文件及CSD模塊計(jì)算的配平量和響應(yīng)進(jìn)行。為更好地理解整個(gè)耦合配平過(guò)程，下面將分條目闡述各個(gè)計(jì)算流程，并對(duì)環(huán)節(jié)中關(guān)鍵點(diǎn)進(jìn)行闡述說(shuō)明。詳細(xì)配平耦合過(guò)程如下：

步驟1 啟動(dòng)計(jì)算，使用飛行力學(xué)解析配平的配平量為基礎(chǔ)，作為CSD綜合分析模型的初值來(lái)使用，基于動(dòng)態(tài)入流模型、旋翼配平量和運(yùn)動(dòng)誘導(dǎo)，通過(guò)LB非定常氣動(dòng)模型結(jié)合查表法計(jì)算旋翼氣動(dòng)力F／M′0，結(jié)合梁模型和葉素理論，并采用有限元方法求解動(dòng)力學(xué)方程解算結(jié)構(gòu)響應(yīng)并進(jìn)一步計(jì)算旋翼力，將此旋翼力代入配平平衡方程組，以牛頓迭代法更新配平量，求解動(dòng)力學(xué)方程組。輸出旋翼操縱量、姿態(tài)角、剛體穩(wěn)態(tài)響應(yīng)及槳葉彈性變形，并以文件形式輸出剖面馬赫數(shù)分布和剖面力系數(shù)，供下次迭代使用。

步驟2 將旋翼總距、橫／縱向周期變距、俯仰角、側(cè)傾角及揮舞擺振的剛體運(yùn)動(dòng)通過(guò)GMP傳遞給OverCFD模塊（俯仰角、側(cè)傾角亦可通過(guò)來(lái)流設(shè)定），彈性變形通過(guò)Motion文件傳遞給Over-CFD模塊。在CFD計(jì)算中，數(shù)值格式采用2階中心差分格式，湍流模型采用SA模型，時(shí)間推進(jìn)選擇ARC3D Diag，限制器選用Koren，耗散項(xiàng)選擇TLNS3D Diss，每步進(jìn)行30次牛頓內(nèi)迭代，初次CFD計(jì)算旋翼槳葉旋轉(zhuǎn)1圈，輸出旋翼槳葉各剖面的垂直力系數(shù)Cn、繞變距軸的力矩系數(shù)Cm、弦向力系數(shù)Cc在槳盤(pán)平面和槳葉展向的分布（F／，供CSD綜合分析模塊調(diào)用。

也可進(jìn)行槳葉旋轉(zhuǎn)1／Nb圈（Nb為槳葉片數(shù)）后計(jì)算所得各剖面氣動(dòng)力，無(wú)需保證每次迭代CFD計(jì)算氣動(dòng)力收斂，只需保證整個(gè)耦合過(guò)程是收斂即可。

步驟3 在CSD綜合分析程序中，通過(guò)讀取上次耦合迭代的CFD求解器計(jì)算氣動(dòng)力和CSD模塊輸出保存的氣動(dòng)力，并將兩者做差，引入旋翼沿展向和周向各剖面氣動(dòng)力增量，即ΔF／M1。

讀取CSD綜合分析程序上次迭代計(jì)算輸出的配平量，作為程序啟動(dòng)值，計(jì)算新的配平量并輸出下次迭代所需要的全部信息。

步驟5 重復(fù)步驟3中CSD綜合分析模塊，即可得到

步驟6 循環(huán)步驟2～步驟5，當(dāng)2次相鄰迭代氣動(dòng)力收斂，即ΔF／Mi≈0，配平操縱和姿態(tài)角收斂，即視為耦合過(guò)程收斂。此時(shí)CFD模塊、配平模塊與氣彈模塊同時(shí)收斂，且

以上耦合流程可歸納為：在每次耦合迭代過(guò)程中通過(guò)ΔF／Mi來(lái)修正查表法得到的氣動(dòng)力，歸納式（12）～式（21）可得

通過(guò)CFD模塊和CSD綜合分析程序求得氣動(dòng)載荷之差而引入增量ΔF／Mi，其在每次迭代過(guò)程中都做了更新。

3 算例驗(yàn)證

3.1 SA349／2直升機(jī)參數(shù)

本文選取SA349／2直升機(jī)穩(wěn)態(tài)前飛的小速度狀態(tài)（前進(jìn)比為0.14）為算例，此狀態(tài)具有明顯的槳渦干擾現(xiàn)象。通過(guò)CFD／CSD耦合方法計(jì)算旋翼氣動(dòng)力，并觀測(cè)流場(chǎng)現(xiàn)象，計(jì)算旋翼槳葉穩(wěn)態(tài)響應(yīng)和振動(dòng)載荷，并將氣動(dòng)載荷和振動(dòng)載荷計(jì)算結(jié)果與飛行實(shí)測(cè)進(jìn)行對(duì)比，驗(yàn)證本文方法的有效性。

SA349／2直升機(jī)為鉸接式旋翼構(gòu)型，翼型采用OA209，圖9給出文獻(xiàn)［15］及本文使用的槳葉網(wǎng)格負(fù)扭轉(zhuǎn)ψ分布。圖中：r／R表示歸一化半徑。

圖9 SA349／2直升機(jī)槳葉負(fù)扭轉(zhuǎn)對(duì)比Fig.9 Comparison of blade twist of SA349／2 helicopter

3.2 CFD／CSD松耦合計(jì)算結(jié)果

采用第2節(jié)闡述的耦合策略進(jìn)行的配平計(jì)算具有較好的收斂速度，CFD模塊與CSD綜合分析模塊耦合4～5次，所得氣動(dòng)力即可收斂，得到的垂直力系數(shù)與實(shí)測(cè)值［15］具有較好的一致性；且輸出的氣動(dòng)載荷同樣捕捉到了槳渦干擾現(xiàn)象，進(jìn)一步反映出本文方法的有效性。表2為此狀態(tài)的配平解。

圖10為槳葉剖面（r／R＝0.75，0.97）處的垂直力系數(shù)與實(shí)測(cè)值隨配平耦合歷程的對(duì)比情況。圖中：實(shí)線為配平迭代過(guò)程每次循環(huán)CFD仿真計(jì)算的氣動(dòng)載荷，虛線為通過(guò)配平程序在采用第2節(jié)中闡述的耦合策略輸出的氣動(dòng)載荷。

表2 SA349／2直升機(jī)飛行狀態(tài)2［15］的配平解Table 2 Trim solution of SA349／2 helicop ter in flight status 2［15］ （°）

圖10 槳葉剖面配平歷程Fig.10 Blade section trimming process

圖11給出了CFD／CSD耦合計(jì)算收斂后剖面垂直力系數(shù)的本文結(jié)果、飛行實(shí)測(cè)值及文獻(xiàn)分析結(jié)果的對(duì)比。圖中：基于文獻(xiàn)配平解的CFD計(jì)算是指使用文獻(xiàn)［15］給出的配平操縱量及響應(yīng)計(jì)算的剛體模型氣動(dòng)載荷，CFDc815指第5次耦合計(jì)算CFD輸出的氣動(dòng)載荷，配平c815指第5次耦合計(jì)算配平模塊輸出的氣動(dòng)載荷，均勻入流模型、預(yù)定尾跡模型及自由尾跡模型數(shù)據(jù)取自CAMRAD［16］計(jì)算值。

從圖11（a）可以分析得出，對(duì)于r／R＝0.75特征剖面，CFD／CSD耦合計(jì)算值與飛行實(shí)測(cè)值具有很好的一致性。CFD／CSD耦合計(jì)算方法計(jì)算的垂直力較CAMRAD具有更好的精度，在本文方法中，每次配平計(jì)算通過(guò)查表法及引入ΔF／Mi增量的方法得到新的配平氣動(dòng)載荷，同樣很好地捕捉到槳渦干擾現(xiàn)象，且與飛行實(shí)測(cè)值相比具有較好的計(jì)算精度。而基于文獻(xiàn)［15］配平解得到的CFD氣動(dòng)力，配平位置與飛行實(shí)測(cè)值相比過(guò)高。

圖11 低速飛行狀態(tài)槳葉剖面垂直力系數(shù)對(duì)比Fig.11 Comparison of vertical force coefficient of blade section in low-speed forward flight

從圖11（b）可以分析得出，CFD／CSD松耦合計(jì)算的垂直力系數(shù)與飛行實(shí)測(cè)值有一定的差距，而基于文獻(xiàn)配平解的計(jì)算與飛行實(shí)測(cè)值具有較好的吻合度，可能是本文的CFD仿真計(jì)算未充分考慮機(jī)身影響，耦合過(guò)程機(jī)身和安定面氣動(dòng)載荷取自文獻(xiàn)［15］的擬合值，可能存在一定誤差。

在低速飛行狀態(tài)下，旋翼的槳渦干擾效應(yīng)對(duì)氣動(dòng)載荷具有極大的影響。從圖10和圖11可以看出，槳尖渦誘導(dǎo)產(chǎn)生的不均勻入流可導(dǎo)致較為突出的氣動(dòng)載荷，其中槳尖位置受槳渦干擾尤為顯著。槳尖渦向槳葉內(nèi)測(cè)移動(dòng)，在前行槳葉90°方位附近造成槳尖氣動(dòng)載荷突增，出現(xiàn)“下-上脈沖”；渦向槳尖移動(dòng)，在后行槳葉270°方位附近造成槳尖氣動(dòng)載荷突降，出現(xiàn)“上-下脈沖”。對(duì)于槳渦干擾現(xiàn)象，低速飛行狀態(tài)尾跡渦計(jì)算的渦量圖可以得到很直觀的反映。

圖12 低速飛行狀態(tài)槳葉剖面垂直力、俯仰力矩和弦向力分布Fig.12 Distribution of blade normal force、pitching moment and in-plane force coefficient in low-speed forward flight

圖12給出了SA 349／2直升機(jī)小速度穩(wěn)態(tài)前飛狀態(tài)旋翼槳葉剖面垂直力Fn、俯仰力矩Fm和弦向力Fc隨展向和周向的分布情況。從圖12（a）中可以直觀看出，在方位角Ψ為90°和270°附近、槳尖段的槳渦干擾現(xiàn)象在槳葉展向r／R＝0.75處垂直力達(dá)到峰值。

圖13 低速飛行狀態(tài)槳葉剖面壓力系數(shù)對(duì)比Fig.13 Comparison of blade surface pressure coefficient in low-speed forward flight

圖13為小速度穩(wěn)態(tài)前飛狀態(tài)，通過(guò)CFD／CSD耦合方法、飛行實(shí)測(cè)值及基于飛行實(shí)測(cè)操縱和響應(yīng)的CFD計(jì)算值，槳葉剖面處壓力系數(shù)的對(duì)比。圖中：-Cp、x／c分別為壓力系數(shù)負(fù)值、歸一化弦向位置；基于文獻(xiàn)配平解的CFD計(jì)算指用文獻(xiàn)［15］中給出的配平量來(lái)進(jìn)行CFD計(jì)算所獲的結(jié)果；CFD／CSD耦合計(jì)算指本文計(jì)算值；上／下表面實(shí)測(cè)值指文獻(xiàn)［15］中給出的槳葉表面壓力系數(shù)實(shí)測(cè)值。通過(guò)曲線對(duì)比不難發(fā)現(xiàn)，對(duì)于下表面壓力系數(shù)，本文方法的計(jì)算值明顯與飛行實(shí)測(cè)值吻合更好。對(duì)于上表面壓力系數(shù)，本文方法的計(jì)算值略低于飛行實(shí)測(cè)值，尤其r／R＝0.97剖面在90°和135°方位角的計(jì)算值，可能是槳渦干擾的影響，造成部分位置計(jì)算結(jié)果誤差。總體來(lái)講，通過(guò)本文方法計(jì)算的壓力系數(shù)分布與實(shí)測(cè)值具有較好的一致性，驗(yàn)證了本文方法的有效性，且對(duì)計(jì)算旋翼穩(wěn)態(tài)前飛狀態(tài)氣動(dòng)力及流場(chǎng)具有很好的精度。對(duì)于計(jì)算結(jié)果的誤差，因?yàn)楸疚慕⒌木W(wǎng)格是對(duì)真實(shí)槳葉的簡(jiǎn)化，槳葉網(wǎng)格不準(zhǔn)確可能帶來(lái)誤差；此外，CFD求解器未考慮低速預(yù)處理也可能導(dǎo)致壓力系數(shù)不精確。

圖14 低速飛行狀態(tài)笛卡兒自動(dòng)背景網(wǎng)格自適應(yīng)計(jì)算旋翼渦量等軸視圖Fig.14 Rotor vorticity isometric view of Cartesian automatic background grid in low-speed forward flight

圖14為通過(guò)本文方法求得的旋翼渦量圖?？梢钥闯?，本文方法對(duì)旋翼槳尖渦及槳渦干擾現(xiàn)象具有很好的捕捉能力。

4 結(jié) 論

1）本文建立的CFD／CSD松耦合策略能有效將基于CFD方法的高保真氣動(dòng)力引入到氣彈綜合分析模塊中，有效提高氣動(dòng)力模型精度的同時(shí)兼顧時(shí)間效率。

2）本文建立的CFD／CSD松耦合配平計(jì)算方法可有效地用于直升機(jī)前飛狀態(tài)下的全機(jī)配平分析，并具有良好的精度、收斂性和穩(wěn)定性。

3）本文改進(jìn)的耦合分析策略對(duì)計(jì)算旋翼穩(wěn)態(tài)前飛狀態(tài)氣動(dòng)力及流場(chǎng)的計(jì)算具有很好的精度，能有效捕捉到槳渦干擾現(xiàn)象等現(xiàn)象，為高精度的氣彈響應(yīng)和載荷分析奠定了良好的基礎(chǔ)。