基于CFD/CSD方法的跨聲速靜氣動(dòng)彈性數(shù)值模擬應(yīng)用研究

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

0 引 言

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

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

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

1 計(jì)算方法

1.1 CFD、CSD控制方程

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

(1)

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

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

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

us=CFs

(2)

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

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

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

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

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

(3)

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

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

2 計(jì)算方法驗(yàn)證

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

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

(a) 外形輪廓

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

2.2 計(jì)算方法驗(yàn)證

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

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

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

圖2 靜氣動(dòng)彈性變形迭代計(jì)算收斂歷程Fig.2 Iteration convergence history of static aeroelastic deformation

(a) CL～α曲線

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

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

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

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

3.1 幾何變形影響

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

(a) 機(jī)翼彎曲變形

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

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

3.2 壓力分布影響

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

圖5 靜氣動(dòng)彈性對(duì)機(jī)翼表面壓力的影響Fig.5 Effects of static aeroelasticity on wing surface pressure

(a) y/b=20%

(b) y/b=95%

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

3.3 氣動(dòng)性能影響

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

圖7 靜氣動(dòng)彈性對(duì)機(jī)翼氣動(dòng)性能的影響Fig.7 Effects of static aeroelasticity on wing’s aerodynamic performances

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

4 結(jié) 論

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

[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]方明祥, 白俊強(qiáng)等. 基于非結(jié)構(gòu)動(dòng)網(wǎng)格的跨音速顫振研究[J]. 機(jī)械設(shè)計(jì)與制造, 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]. 中國(guó)科學(xué)A輯, 2002, 32(8): 716-724

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

[16]趙永輝. 氣動(dòng)彈性力學(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)國(guó). 飛行器結(jié)構(gòu)學(xué)(第2版)[M]. 北京: 北京航空航天大學(xué)出版社, 2010.