郭洪濤, 陳德華, 張昌榮, 呂彬彬, 王曉冰(. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室, 四川 綿陽 6000；. 中國空氣動力研究與發(fā)展中心 高速空氣動力研究所, 四川 北川 6763)

0 引 言

現(xiàn)代先進大型飛機的巡航速域通常在跨聲速范圍。飛機升力線斜率在跨聲速時最大且變化最劇烈，此時流固耦合作用也最敏感，結(jié)構(gòu)上即使一個細微變化也可能對飛機的氣動特性造成很大影響[1]。試驗表明：在相同速壓情況下，與亞、超聲速比較，跨聲速范圍內(nèi)靜氣動彈性造成的飛機升力損失最大，使焦點前移最甚[2]。隨著氣動彈性計算技術(shù)的進步，近年來國內(nèi)外在這方面已經(jīng)開展了大量的研究工作。

70年代，基于線化位流理論發(fā)展起來的面元法開始廣泛應(yīng)用于飛機的氣動彈性計算，現(xiàn)在采用該方法的NASTRAN等商業(yè)軟件仍被許多工程單位廣泛使用[3]。1974年，Ballhaus與Steger提出了一種二維全隱式有限差分格式構(gòu)造方法，這個方法導(dǎo)致了第一個專業(yè)氣動彈性計算軟件LTRANS的誕生。80年代初期，跨聲速小擾動方程開始被用于跨聲速氣動彈性計算問題，波音公司開發(fā)了XTRAN3S，實際上是把LTRANS從二維擴展到了三維。整個80年代，全速勢方程加邊界層修正方法成為民用飛機設(shè)計的主流CFD方法。1985年，為了滿足未來空軍戰(zhàn)斗機的需求，AFWAL與NASA合作開發(fā)了ATRAN3S，主要是考慮到小展弦比戰(zhàn)斗機機翼的氣動彈性計算[4]。上述方法對氣動彈性計算的發(fā)展做出了許多開創(chuàng)性的貢獻，但不足的是它們都是基于位流理論假設(shè)。當(dāng)遭遇比較復(fù)雜的流動現(xiàn)象，比如當(dāng)有較強的激波通過翼面或者有流動分離時，這種假設(shè)就不太適合了。80年代到90年代，基于Euler方程的各種激波捕捉格式得到發(fā)展，在被應(yīng)用到飛行器氣動力計算的同時，也被運用到飛行器氣動彈性計算中來[5-6]。隨著90年代基于雷諾平均N-S方程(簡稱RANS方程)求解方法的成熟，既可以考慮跨聲速流動中的激波間斷又可以考慮氣動力湍流黏性效應(yīng)的氣動彈性數(shù)值計算方法也得到了有效的發(fā)展[7-8]。Robinson與Bantina在CFL3D中加入了氣動彈性計算的相關(guān)模塊[9]，Ames研究中心也開發(fā)了專門的氣動彈性計算軟件ENSAERO，這些軟件采用Euler/N-S方程求解氣動力，采用隱式B-W格式，B-L、k-ε等湍流模型，采用模態(tài)疊加法求解結(jié)構(gòu)運動方程，可用于解決復(fù)雜外形的氣動力計算和氣動彈性分析研究[10]。但是，這類方法基于模態(tài)運動方程求解結(jié)構(gòu)變形，對于靜氣動彈性計算而言效率太低。近年來，具有更高精度的大渦模擬方法和直接數(shù)值模擬方法在國內(nèi)外也被廣泛研究[11-15]，但這兩種方法因為計算量巨大，目前僅適合于一些基礎(chǔ)理論研究，并不適合于飛機等具有復(fù)雜邊界條件物體的靜氣動彈性數(shù)值計算。

綜上所述，對于跨聲速流動且物面邊界條件復(fù)雜的情況，目前還沒有一套權(quán)威、統(tǒng)一的計算分析方法和代碼，還值得深入細致地研究。本文通過耦合求解RANS方程與靜氣動彈性平衡方程，采用結(jié)構(gòu)化動網(wǎng)格技術(shù)和多物理場數(shù)據(jù)插值技術(shù)，實現(xiàn)了跨聲速復(fù)雜外形飛行器的靜氣動彈性數(shù)值預(yù)測分析，并基于典型風(fēng)洞試驗結(jié)果驗證了計算方法、代碼的有效性。此外，基于數(shù)值模擬結(jié)果，分析了靜氣動彈性對典型大型飛機機翼幾何變形特性、表面壓力分布以及氣動性能的影響。

1 計算方法

1.1 CFD、CSD控制方程

控制方程采用時間相關(guān)的三維守恒型可壓縮RANS方程，在一般曲線坐標(biāo)系(ξ、η、ζ)下，其無量綱形式為：

(1)

式中：t為時間，Q為守恒變量，F(xiàn)、G、H為無黏矢通量，F(xiàn)v、Gv、Hv為黏性矢通量。

湍流模擬選用一方程的S-A模型。物面為絕熱無滑移邊界條件，遠場為壓力遠場無反射邊界條件。為了加速CFD計算收斂，應(yīng)用了多重網(wǎng)格技術(shù)。

物面的結(jié)構(gòu)彈性變形采用靜氣動彈性平衡方程來求解。該方程為常系數(shù)齊次線性方程組，實際上相當(dāng)于模態(tài)運動微分方程組沒有了時間導(dǎo)數(shù)項，相對而言，計算量大大減小。即：

us=CFs

(2)

式中：us為結(jié)構(gòu)點的變形位移矢量，C為結(jié)構(gòu)點柔度影響系數(shù)矩陣，F(xiàn)s為作用在結(jié)構(gòu)點上的氣動力、質(zhì)量力、發(fā)動機推力等合力矢量。

1.2 流固耦合數(shù)據(jù)插值方法

大型飛機靜氣動彈性屬于三維空間的小變形情況，綜合內(nèi)存占用率、計算效率以及插值精度考慮，采用三維薄板樣條插值方法(TPS)是比較適宜的，TPS的詳細計算方法可見參考文獻[16]。

1.3 動網(wǎng)格生成方法

針對結(jié)構(gòu)化流場計算網(wǎng)格，本文采用RBF結(jié)合TFI的方法來進行動網(wǎng)格生成。RBF可以視為面樣條函數(shù)插值方法的三維擴充。其插值公式為：

(3)

其中，ri=(xi,yi,zi)為位移已知點，其數(shù)目為n，φ為關(guān)于空間距離‖r-ri‖的基函數(shù)，本文取φ(‖r-ri‖)=‖r-ri‖3、ψ=b0+b1x+b2y+b3z，插值公式的系數(shù)可通過已知點ri的位移di和平衡條件得到。詳細計算方法見參考文獻[17]。

基于RBF方法插值得到各網(wǎng)格塊邊內(nèi)點的位移后，各網(wǎng)格塊面內(nèi)點的位移可基于網(wǎng)格塊邊通過TFI方法插值得到。同理，各網(wǎng)格塊內(nèi)部點的位移也可基于網(wǎng)格塊面通過TFI方法插值得到[18]。

2 計算方法驗證

2.1 計算模型與CFD網(wǎng)格

本文對某大型飛機翼身組合體模型進行了靜氣動彈性數(shù)值計算，并與風(fēng)洞試驗結(jié)果進行了對比驗證。圖1給出了計算模型的外形輪廓及CFD物面網(wǎng)格。該模型采用了超臨界擴張后緣翼型以及大展弦比后掠下單翼布局，忽略了短艙/掛架、襟翼滑軌整流罩等部件，計算網(wǎng)格單元規(guī)模為500萬。

(a) 外形輪廓

(b) 初始物面網(wǎng)格

2.2 計算方法驗證

驗證計算狀態(tài)為：Ma=0.78，q=35 kPa，雷諾數(shù)Re=7.3×106，迎角α為-4°～8°。

靜氣動彈性計算是一個CFD/CSD的反復(fù)耦合迭代計算過程。對于工程應(yīng)用來說，需要迭代計算效率高，收斂速度快?；诒疚姆椒ǎ粋€計算狀態(tài)只需要6～7個變形迭代步就基本上收斂了，如圖2所示。在計算過程中，由于流場計算網(wǎng)格的拓撲結(jié)構(gòu)與網(wǎng)格數(shù)量不會改變，于是每個變形步后的CFD計算都可以利用上個變形步的CFD結(jié)果作為初值進行續(xù)算，因此，中間變形迭代步的CFD計算只需要較少的計算步數(shù)。一般情況下，本文CFD/CSD方法的計算量是對應(yīng)的剛性模型CFD計算量的1.5倍～2倍。

圖3給出了本文計算結(jié)果與相應(yīng)風(fēng)洞試驗結(jié)果的機翼升力系數(shù)CL和彎曲變形Δz的對比結(jié)果。

圖2 靜氣動彈性變形迭代計算收斂歷程Fig.2 Iteration convergence history of static aeroelastic deformation

(a) CL～α曲線

(b) 不同迎角下Δz～y/b曲線

圖3中給出的CL進行了計算與試驗的相關(guān)性修正。主要是修正了邊界層轉(zhuǎn)捩模擬不一致、風(fēng)洞試驗洞壁干擾、計算模型與試驗?zāi)Ｐ屯庑尾町惖龋芟抻诋?dāng)前的技術(shù)手段，還無法對CFD中的湍流模擬不足與分離特性模擬不準(zhǔn)以及實際機翼變形可能存在的非線性因素等影響進行修正?？傮w上看，本文計算結(jié)果與風(fēng)洞試驗結(jié)果的一致性較好。CL的差異主要出現(xiàn)在翼面氣流出現(xiàn)分離之后的CFD模擬誤差所致；Δz的差異主要存在于小變形范圍，根據(jù)機翼變形應(yīng)該是連續(xù)一致的規(guī)律性推測，應(yīng)該是小變形時測量設(shè)備的測量誤差相對較大。

3 跨聲速靜氣動彈性數(shù)值模擬分析

大型飛機為了提高氣動效率，其巡航Ma通常位于跨聲速范圍。在跨聲速時，即使機翼迎角很小，由于激波與邊界層的強相互作用，黏性影響也不應(yīng)該被忽略。因此，必須利用N-S方程計算氣動力。不失一般性，本文重點考察了Ma=0.78時算例模型的靜氣動彈性特性。

3.1 幾何變形影響

圖4給出了Ma=0.78、q=35 kPa、α∈[-4°～8°]時機翼的Δz與Δε(彈性扭轉(zhuǎn)角)沿展向變化曲線。

(a) 機翼彎曲變形

(b) 機翼彈性扭轉(zhuǎn)變形

可以看出，當(dāng)迎角為正時，靜氣動彈性使機翼產(chǎn)生向上的彎曲撓度變形，并在機翼的順氣流剖面產(chǎn)生負的彈性扭轉(zhuǎn)角。從圖中還可以看出，內(nèi)翼段的變形量較小，在機翼拐折處以外機翼的彈性變形量越來越大，到翼尖處變形量達到最大值。這種變化規(guī)律與機翼剛度沿展向逐步減小的分布特性是相符的。根據(jù)參考文獻[19]可知，對于后掠機翼來說，這種變形特征將減小機翼沿展向各剖面的當(dāng)?shù)赜?，從而影響載荷分布，改變機翼氣動特性，產(chǎn)生所謂的靜氣動彈性效應(yīng)。例如，依據(jù)本文算例，在Ma=0.78、q=35 kPa、α=2°時，翼尖的彎曲變形可以達到機翼半展長的4%，翼尖剖面的順流向彈性扭轉(zhuǎn)角與順流向迎角變化量分別可達-2°與-4°。

3.2 壓力分布影響

圖5與圖6分別給出了Ma=0.78、q=35 kPa、α=2°時剛性與彈性機翼的表面壓力等值線分布圖及不同展向位置的壓力系數(shù)分布圖。從圖6可以看出，與機翼的幾何變形特性相對應(yīng)，靜氣動彈性對翼根的壓力分布特性影響較小，對翼尖的壓力分布特性影響較大，比較明顯的影響是減少了靠近外側(cè)機翼上表面的負壓范圍。

圖5 靜氣動彈性對機翼表面壓力的影響Fig.5 Effects of static aeroelasticity on wing surface pressure

(a) y/b=20%

(b) y/b=95%

從圖6中機翼不同剖面變形前后的壓力系數(shù)對比可知，越靠近翼尖時壓力系數(shù)分布變化越大。這是由于在機翼小迎角(表面未出現(xiàn)大面積的氣流分離)情況下，后掠彈性機翼的氣動載荷使得機翼向上彎曲，同時由于負的彈性扭轉(zhuǎn)角，導(dǎo)致順氣流方向的當(dāng)?shù)赜菧p小，使得彈性變形后機翼的前緣吸力峰值降低，同時使得彈性機翼的升力系數(shù)小于剛性機翼的升力系數(shù)。依據(jù)本文算例，在Ma=0.78、q=35 kPa、α=2°時，彈性相對于剛性機翼的CL降低了15%以上。因此，在跨聲速時靜氣動彈性對大展弦比后掠翼的氣動載荷影響不應(yīng)該被忽略。

3.3 氣動性能影響

圖7給出了Ma=0.78時彈性機翼的典型靜氣動彈性影響因子隨速壓變化曲線?？梢钥闯觯賶簩椥詸C翼升阻特性及縱向靜穩(wěn)定性都有不同程度的影響。其中，K1(升力線斜率影響因子)隨速壓基本上呈線性減小趨勢，在q=35 kPa(對應(yīng)巡航飛行高度)時影響因子大約為0.94，但這種線性規(guī)律性將在速壓達到某一較大值后發(fā)生變化，當(dāng)q>65 kPa時，K1～q曲線開始出現(xiàn)拐折。K2(縱向靜穩(wěn)定性裕度影響因子)隨速壓增大總體上呈現(xiàn)出非線性減小趨勢，也即是說，速壓越大，彈性機翼的靜穩(wěn)定性裕度越小。但是，在q>65 kPa后，K2～q曲線開始出現(xiàn)上翹，即影響特性發(fā)生了逆向變化。對于K3(巡航升阻比影響因子)來說，盡管靜氣動彈性會降低機翼的升力，但是由于降低升力的同時也降低了升致阻力，因此實際上會增大升阻比。從圖中可以看出，K3在小速壓情況下呈線性變化規(guī)律，但在q>45 kPa以后非線性特征逐漸顯現(xiàn)。

圖7 靜氣動彈性對機翼氣動性能的影響Fig.7 Effects of static aeroelasticity on wing’s aerodynamic performances

綜合跨聲速靜氣動彈性對大型飛機縱向氣動特性的影響來看，當(dāng)速壓較小時，靜氣動彈性影響量隨速壓呈線性變化，當(dāng)速壓增大到一定程度后，變化規(guī)律將出現(xiàn)非線性特征。因此，大速壓條件下的靜氣動彈性預(yù)測應(yīng)當(dāng)加密計算工況點，且不宜用小速壓條件下獲得的變化規(guī)律外推至大速壓范圍。

4 結(jié) 論

本文基于RANS方程與靜氣動彈性平衡方程發(fā)展了一種CFD/CSD流固耦合靜氣動彈性計算方法。為了提升計算精準(zhǔn)度、魯棒性及效率，CFD采用了多塊對接結(jié)構(gòu)化網(wǎng)格進行分區(qū)并行計算，并利用多重網(wǎng)格技術(shù)來加快計算收斂速度，利用RBF結(jié)合TFI技術(shù)來生成結(jié)構(gòu)化動網(wǎng)格，通過TPS方法來進行氣動/結(jié)構(gòu)的數(shù)據(jù)插值，并開發(fā)了相應(yīng)代碼給予實現(xiàn)。通過典型風(fēng)洞試驗結(jié)果驗證了方法和代碼的有效性。基于某大型飛機翼身組合體模型，數(shù)值模擬研究了跨聲速時某典型大展弦比機翼的靜氣動彈性特性，獲得了一些具有代表性的結(jié)論，可以提供給大型飛機氣動/結(jié)構(gòu)設(shè)計人員參考。

[1]Thomas G I. Unique testing capabilities of the NASA langley transonic dynamics tunnel, an exercise in aeroelastic scaling[R]. AIAA 2013-2625

[2]Mawilkinson D G, Blackerby W T, et al. Correlation of full-scale drag predictions with flight measurements on the C-141A aircraft-phase II, wind tunnel test, analysis, and prediction techniques[R]. NASA Technical Paper 2333, Feb. 1974

[3]Rodden W P, Johnson E H. MSC/Nastran Aeroelastic analysis user’s guide V68[M]. Los Angeles: MSC Software Corporation, 1994: 44-65

[4]Tian B Y. Computational aeroelastic analysis of airplane wings including geometry nonlinearity[D]. Graduate Education and Research University of Cincinnati, 2003, Paper for the degree of Doctor

[5]Zhang Z C, Liu F. Calculations of unsteady flow and flutter by an Euler and integral boundary-layer method on Cartesian grids[R]. AIAA 2004-5203

[6]方明祥, 白俊強等. 基于非結(jié)構(gòu)動網(wǎng)格的跨音速顫振研究[J]. 機械設(shè)計與制造, 2007, 11: 11-19

[7]Livne E, Weisshaar T A. Aeroelasticity of nonconventional configurations: past and future[J]. Journal of Airplane, 2003, 40(6): 1047-1065

[8]Ramji Kamakoti, Wei S, et al. Time dependent RANS computation for an aeroelastic wing[R]. AIAA 2004-0886

[9]Robinson B A, Batina J T, et al. Aeroelastic analysis of wings using the Euler equation with a deforming mesh[J]. Journal of Airplane, 1991, 28(11): 781-788

[10]Guruswamy G P. Unsteady aerodynamic and aeroelastic calculations for wings using Euler equations[J]. AIAA Journal, 1990, 28(3): 461-469

[11]Rizzetta D P, Visbal M R, Gaitonde D V. Direct numerical and large-eddy simulation of supersonic flows by a high-order method[R]. AIAA 2000-2408

[12]Guarini S E, Moser R D, et al. Direct numerical simulation of a supersonic turbulent boundary layer at Mach 2.5[J]. Fluid Mech, 2000, 414: 1-33

[13]Weber C, Ducros F. A large eddy simulation of complex turbulent flows[R]. AIAA 98-2651

[14]李新亮, 傅德薰, 馬延文. 可壓縮均勻各向同性湍流直接數(shù)值模擬[J]. 中國科學(xué)A輯, 2002, 32(8): 716-724

[15]馬漢東, 潘宏祿, 王強. 超聲速平板邊界層斜波失穩(wěn)轉(zhuǎn)捩過程研究[J]. 力學(xué)學(xué)報, 2007, 39(2): 153-157

[16]趙永輝. 氣動彈性力學(xué)與控制[M]. 北京: 科學(xué)出版社, 2007

[17]Holger W. Computational aspects of radial basis function approximation[J]. Studies in Computational Mathematics, 2006, 2: 231-256

[18]Soni B K. Two- and three-dimensional grid generation for internal flow applications of computational fluid dynamics[C]//7th Computational Physics Conference, Cincinnati, OH, AIAA 1985-1526

[19]酈正能, 程小全, 方衛(wèi)國. 飛行器結(jié)構(gòu)學(xué)(第2版)[M]. 北京: 北京航空航天大學(xué)出版社, 2010.