郭秋亭，孫巖，郭正，劉光遠

1. 國防科技大學 空天科學學院，長沙 410073

2. 中國空氣動力研究與發(fā)展中心 高速空氣動力研究所，綿陽 622762

3. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

一個多世紀以來，風洞試驗伴隨著飛行器的誕生，在飛行器研制中發(fā)揮了推動性的作用，目前仍然是開展飛行器設計、性能評估和流動機理研究的重要手段。高質量的風洞試驗數(shù)據(jù)能夠有效降低設計風險、縮短研制周期和挖掘飛行器設計潛能。但受設備物理條件的限制，風洞試驗難以對全部的流場相似參數(shù)進行模擬。例如，高速風洞試驗中，由于試驗段尺寸和模型縮比的影響，試驗雷諾數(shù)通常要比真實飛行雷諾數(shù)低1～3個量級[1]。

雷諾數(shù)是描述流動尤其是壁面附近邊界層內流動演化的重要相似參數(shù)[2]。早期由于條件限制，飛行器的風洞試驗主要模擬馬赫數(shù)，較少考慮雷諾數(shù)的影響，為此付出過沉重的代價。20世紀50年代，美國空軍研制的C-141運輸機在設計時沒有妥善修正雷諾數(shù)對氣動性能的影響，焦點位置預判失準，險些釀成重大的飛行事故[3]。

現(xiàn)代大型運輸類飛行器為了提高飛行品質、降低油耗和增大航程，廣泛采用跨聲速性能更優(yōu)異的超臨界翼型。但這類翼型對雷諾數(shù)變化十分敏感，真實飛行雷諾數(shù)與試驗雷諾數(shù)下的氣動性能差異明顯[4]。為了探究飛行雷諾數(shù)下的流動機理、推動先進大型運輸類飛行器的研制，美國和歐洲先后建成了具有飛行雷諾數(shù)模擬能力的低溫跨聲速風洞，分別為國家跨聲速風洞設備[5](National Transonic Facility, NTF)和歐洲跨聲速風洞[6](European Transonic Windtunnel, ETW)，能夠準確預測雷諾數(shù)和靜氣動彈性影響。

低溫高雷諾數(shù)風洞主要通過降低試驗段氣流總溫與增大氣流總壓相結合的方式來提高風洞試驗的雷諾數(shù)。但增壓的同時也增大了速壓，增加了風洞模型受到的氣動載荷，導致模型出現(xiàn)較大的結構變形，引入了不可忽視的數(shù)據(jù)誤差。數(shù)值計算和風洞試驗研究結果均表明，模型結構變形與雷諾數(shù)變化引起的氣動載荷變化量相近，且多數(shù)情況下二者的作用相反，使得模型靜彈性變形掩蓋了雷諾數(shù)對氣動特性的影響，造成偽雷諾數(shù)效應[7-9]。低溫高雷諾數(shù)風洞的總溫、總壓可以獨立控制，能夠保持速壓不變，開展不同雷諾數(shù)下的風洞試驗，單獨研究雷諾數(shù)對氣動特性的影響；也可以保持雷諾數(shù)不變，開展不同速壓下的風洞試驗，單獨研究模型靜彈性變形對氣動特性的影響，從而將雷諾數(shù)與靜彈性變形的影響進行分離[10]。

國內尚未建成總溫、總壓可獨立控制的風洞設備，高雷諾數(shù)試驗主要通過增壓方式進行[11]，結構靜彈性變形影響難以從試驗數(shù)據(jù)分離。目前，對于增壓試驗中雷諾數(shù)/靜氣彈效應分離問題，尚無有效的解決手段?；跀?shù)值模擬或變形測量的模型變形影響修正技術在復雜外形上，仍然存在操作上的復雜性和結果的不確定性[12-14]。

為此，通過分析飛行器氣動特性隨雷諾數(shù)和靜氣彈變形的變化規(guī)律，提出一種基于風洞試驗數(shù)據(jù)的雷諾數(shù)/靜氣彈效應快速分離方法，為風洞增壓試驗提供一種數(shù)據(jù)修正技術。

1 數(shù)值計算方法

1.1 NNW-FSI軟件

模型繞流流場和氣動彈性數(shù)值模擬均在NNW-FSI軟件平臺[15]開展。對湍流的模擬采用有限體積方法離散求解雷諾平均奈維爾-斯托克斯方程，該方程中的雷諾應力項通過SST(Shear Stress Transport) 渦粘模型[16]求解。氣動/結構耦合數(shù)值模擬中結構靜變形計算采用柔度矩陣方法，耦合界面數(shù)據(jù)傳遞采用薄板樣條插值方法，網(wǎng)格運動采用徑向基函數(shù)-超限插值復合方法實現(xiàn)。

NNW-FSI是國家數(shù)值風洞多學科耦合應用軟件系統(tǒng)下規(guī)劃的一款流固耦合模擬軟件，在“in-house”求解軟件TRIP[17]近20年研究和開發(fā)工作的基礎上，通過模塊化設計、組件封裝、代碼重構和功能拓展，實現(xiàn)了體系化設計和基礎版本的開發(fā)。NNW-FSI軟件包含前置處理、后置處理、核心解算和運行環(huán)境管理4個主要功能模塊，能夠開展定常/非定常流動、靜態(tài)氣動/結構耦合、顫振、抖振等流動現(xiàn)象的數(shù)值仿真，具有較高的可信度[18-19]。

1.2 計算準確性驗證

以高雷諾數(shù)氣動結構(High Reynolds Number Aero-Structural Dynamics,HIRENASD)機翼對靜氣動彈性耦合計算結果檢驗數(shù)值計算方法的準確性。HIRENASD機翼是氣動彈性預測會議(Aeroelastic Predication Workshop, AePW)選擇的標準外形，采用典型超臨界翼型，氣動外形和幾何參數(shù)見圖1[20]。HIRENASD機翼由3段組成，機翼的前緣后掠角為34°、半翼展b/2為1.286 m、機翼參考面積Aref為0.392 6 m2、平均氣動弦長cref為0.344 5 m。

圖1 HIRENASD機翼幾何信息[20]

HIRENASD機翼計算網(wǎng)格采用粗密度多塊結構對接網(wǎng)格，如圖2所示。流場網(wǎng)格采用H型拓撲結構，機身和機翼周圍分別包裹O型網(wǎng)格用于模擬邊界層內部流動，機身表面第1層網(wǎng)格高度為6.0×10-7m，機翼表面第1層網(wǎng)格高度為4.4×10-7m，整個網(wǎng)格包含319個網(wǎng)格塊、3 158 849 個網(wǎng)格點、3 088 384個網(wǎng)格單元。

圖2 HIRENASD模型流場計算網(wǎng)格

圖3給出了HIRENASD機翼彎曲變形計算與試驗結果[21]的對比曲線。其中,Ma、Re和α分別表示馬赫數(shù)、雷諾數(shù)和迎角,Ma=0.80,Re=1.4×107，α=3.0°;q/E為無量綱載荷因子，定義為速壓q與結構彈性模量E的比值，用于表征氣動載荷作用下結構變形的大小,q/E=4.7×10-7。圖3中TE(Trailing Edge)表示機翼后緣，LE(Leading Edge)表示機翼前緣，Δy表示機翼彎曲位移，η為歸一化的機翼展向坐標。結果表明，計算與試驗得到的機翼變形沿展向分布從趨勢與量值上都比較一致，說明采用的靜氣動彈性耦合計算方法能夠準確預測跨聲速繞流流場和模型變形情況。

圖3 HIRENASD機翼彎曲變形

2 雷諾數(shù)/靜氣動彈性效應

在驗證了數(shù)值模擬方法準確性的基礎上，仍以1.2節(jié)的HIRENASD機翼為研究對象，進一步分析雷諾數(shù)/靜氣動彈性變形對模型氣動系數(shù)的影響規(guī)律，為分離方法的建立提供理論支撐。

2.1 有限元模型與計算條件

圖4(a)給出了HIRENASD機翼的結構有限元模型，共包含193 457個網(wǎng)格點、103 458個四面體網(wǎng)格單元，固支約束添加在與天平連接的端面。為了降低結構靜變形計算中的柔度矩陣規(guī)模，提高耦合數(shù)據(jù)傳遞的效率與魯棒性，選擇有限元模型表面359個網(wǎng)格點構建簡化柔度矩陣，如圖4(b)所示。

圖4 HIRENASD機翼有限元模型

流場與耦合計算條件依據(jù)總溫不可調節(jié)風洞的試驗條件給定，即雷諾數(shù)調節(jié)只能通過總壓變化實現(xiàn)，雷諾數(shù)Re與載荷因子q/E之間為線性變化的關系。

計算條件采用3組不同的雷諾數(shù)Re，對應不同的載荷因子q/E，詳細的計算條件參見表1。

表1 HIRENASD機翼計算條件

2.2 雷諾數(shù)/靜氣動彈性效應

圖5給出了HIRENASD機翼氣動力系數(shù)變化量ΔCL、ΔCD隨迎角的變化曲線，其中的氣動力系數(shù)變化量均是以Re為0.7×107時剛性外形的氣動力系數(shù)為基準。

從圖5可以看出，剛性機翼(Rigid)隨著雷諾數(shù)Re的增加，升力系數(shù)CL增大，阻力系數(shù)CD減小。對于HIRENASD機翼所采用的超臨界翼型，雷諾數(shù)Re增加會減小邊界層厚度，使得機翼的等效彎度增大，從而升力系數(shù)增大；而Re增加會減小機翼受到的摩擦阻力，從而使阻力系數(shù)CD降低。但當迎角大于一定值，機翼上翼面開始出現(xiàn)流動分離，雷諾數(shù)Re對邊界層流動的影響減弱，引起的升力系數(shù)變化量ΔCL也隨著迎角增加開始減小，同時機翼的阻力主要由分離流動引起的壓差阻力貢獻，摩擦阻力在阻力中的占比已經(jīng)很小，表現(xiàn)為迎角較大情況下Re對阻力系數(shù)CD的影響非常小。

圖5 HIRENASD機翼氣動力系數(shù)變化量

而對于彈性機翼(Elastic)，氣動力系數(shù)隨雷諾數(shù)Re的變化與剛性狀態(tài)下表現(xiàn)出完全迥異的情況。靜氣動彈性變形對升力系數(shù)CL的影響與雷諾數(shù)Re的影響相反，使得升力系數(shù)降低，且兩者的影響量級接近，因此靜氣動彈性變形可能完全掩蓋住Re對升力系數(shù)CL的真實影響，造成虛假的雷諾數(shù)效應。而對于阻力系數(shù)CD，靜氣動彈性變形對CD的影響與Re相同，使得阻力系數(shù)減小，但與雷諾數(shù)Re不同的是，靜氣動彈性變形主要影響的是機翼的誘導阻力，且阻力系數(shù)變化量隨著迎角增加急劇增大，機翼靜氣動彈性變形將放大雷諾數(shù)Re對阻力系數(shù)CD的影響。

2.3 氣動力系數(shù)隨參數(shù)的變化規(guī)律

圖6給出了幾個特定迎角下剛性機翼氣動力系數(shù)變化量ΔCL和ΔCD隨lg(Re)的變化曲線，剛性機翼氣動力系數(shù)變化量為相對于Re為0.7×107時剛性機翼氣動力系數(shù)的差量。這里選擇氣動力系數(shù)變化量是因為變化量與絕對量之間只相差一個偏移量，二者隨lg(Re)的變化規(guī)律相同，采用變化量可以將不同迎角下的曲線繪制在同一幅圖內，并不影響對剛性機翼氣動力系數(shù)隨參數(shù)變化規(guī)律的分析。從圖6可以看出，對于HIRENASD機翼，在當前的計算狀態(tài)下，剛性機翼氣動力系數(shù)或氣動力系數(shù)變化量與雷諾數(shù)的對數(shù)lg(Re)之間表現(xiàn)為強線性變化關系。

圖6 剛性機翼氣動力系數(shù)變化量隨lg(Re)的變化

圖7實線部分給出了迎角α=-1.0°,1.0°,3.0°,5.0°時，彈性機翼氣動力系數(shù)變化量ΔCf隨載荷因子q/E的變化曲線，這里的氣動力系數(shù)變化量是指機翼靜氣動彈性變形引起的氣動力系數(shù)改變量，為每一個對應雷諾數(shù)Re和載荷因子q/E下機翼彈性外形與剛性外形氣動力系數(shù)的差值。可以看出，彈性變形引起的氣動力系數(shù)變化量與載荷因子之間也呈現(xiàn)出強烈的線性變化關系。

圖7 彈性機翼氣動力系數(shù)變化量隨載荷因子q/E的變化

3 分離方法

在雷諾數(shù)效應風洞試驗中，不考慮其它參數(shù)的變化，彈性外形下的氣動力系數(shù)Cf,elastic可以表示為雷諾數(shù)Re和速壓q的函數(shù)：

Cf,elastic=f(Re,q)

(1)

速壓會引起機翼結構的靜氣動彈性變形，因此可以將氣動力系數(shù)分成剛性外形氣動力系數(shù)Cf,rigid與機翼靜氣動彈性變形引起氣動力系數(shù)變化量ΔCf,elastic2個部分，其中Cf,rigid僅隨Re變化，ΔCf,elastic隨Re和q變化，關系式為

Cf,elastic(Re,q)=Cf,rigid(Re)+ΔCf,elastic(Re,q)

(2)

為了將雷諾數(shù)Re與速壓q引起的靜氣動彈性變形效應分離，提出2個假設：

假設1模型靜氣動彈性變形引起的氣動力系數(shù)變化量ΔCf,elastic隨載荷因子q/E呈線性變化關系。

假設2雷諾數(shù)Re變化引起ΔCf,elastic的改變量很小，可以忽略不計。

對于假設1，機翼結構在線性小變形范圍內，如果不考慮變形前后氣動載荷分布的變化，靜氣動彈性變形引起的氣動力系數(shù)變化量ΔCf,elastic正比于機翼的結構靜變形量，而變形量正比于結構受到的氣動載荷、反比于結構材料的彈性模量E，而氣動載荷又正比于速度壓力q，因此推斷ΔCf,elastic與q/E表現(xiàn)為線性比例關系，2.3節(jié)中的圖7也印證了這個觀點。而當載荷因子q/E等于零時，機翼不受到氣動載荷的作用，靜氣動彈性變形為零，因此ΔCf,elastic也等于零。

對于假設2，從雷諾數(shù)Re對升力系數(shù)的影響可以看到，Re增大會使升力增加，一般情況下增加量相對于絕對量是小量，而機翼靜氣動彈性變形主要由機翼升力載荷引起，因此Re增大對靜氣動彈性變形的影響很小，由此推斷Re變化對ΔCf,elastic的影響量非常小。

根據(jù)假設1與假設2，式(2)可以變化為

Cf,elastic(Re,λ)=Cf,rigid(Re)+k2λ

(3)

式中：λ=107×q/E；k2為待確定的系數(shù)。由式(3) 可以看出，雷諾數(shù)/靜氣動彈性效應分離的關鍵是確定系數(shù)k2，進而確定ΔCf,elastic。但通過式(3)無法直接確定比例系數(shù)k2，因此提出第3個假設，通過給定Cf,rigid隨雷諾數(shù)Re的函數(shù)關系確定系數(shù)k2。

假設3剛性外形氣動力系數(shù)Cf,rigid隨雷諾數(shù)的對數(shù)lg(Re)呈線性變化關系。

飛行器氣動力特性的雷諾數(shù)效應研究結果表明[22]，對于給定的飛行器外形和繞流條件，存在一個雷諾數(shù)的自準區(qū)，在該范圍內，氣動力特性與雷諾數(shù)的對數(shù)呈現(xiàn)近似線性變化的關系，第2節(jié)的數(shù)值模擬結果(圖6)也十分吻合這種規(guī)律。風洞試驗中，為了外插獲取飛行雷諾數(shù)下的氣動數(shù)據(jù)，試驗條件通常會根據(jù)經(jīng)驗選擇自準區(qū)范圍內的雷諾數(shù)。假設3用于確定式(3)中的系數(shù)k2，但在有些情況下假設3可能并不嚴格滿足，后文將對不滿足情況的影響做進一步詳細的分析。

基于假設3，式(3)可以變換為

Cf,elastic(ξ,λ)=k1ξ+k2λ+k3

(4)

式中：ξ=lg(Re/Re0)；Re0=106；k1、k3為待確定的系數(shù)。利用3組不同雷諾數(shù)Re與載荷因子q/E組合下的氣動力試驗數(shù)據(jù)即可求解線性方程組確定系數(shù)k1、k2、k3為

(5)

當不同雷諾數(shù)Re與載荷因子q/E的氣動力數(shù)據(jù)超過3組，可以利用最小二乘方法計算待定系數(shù)k1、k2、k3。

待定系數(shù)k1、k2、k3確定后，利用式(4)即可得到氣動力系數(shù)隨雷諾數(shù)與載荷因子的變化特性，其中式(4)右端第1項表示隨Re變化的氣動力系數(shù)，右端第2項表示模型靜氣動彈性變形引起的氣動力系數(shù)變化量，第3項k3表示雷諾數(shù)Re等于Re0時的剛性氣動力系數(shù)，即

(6)

4 分離方法測試

采用大展弦比翼/身/平尾組合體風洞模型的試驗數(shù)據(jù)對同步分離方法進行測試，并通過不同載荷因子外插得到的剛性外形風洞試驗數(shù)據(jù)對分離方法獲得的結果進行驗證和評估。

4.1 測試對象及數(shù)據(jù)

測試對象選擇NASA設計的通用研究模型(Common Research Model,CRM)的翼/身/平尾組合體構型CRM-WBT0，外形如圖8所示。CRM模型代表了一類典型的現(xiàn)代寬體客機外形，其設計巡航馬赫數(shù)Ma為0.85，巡航升力系數(shù)為0.50，是第4～6屆阻力預測會議(Drag Prediction Workshop,DPW)選擇的標準外形。CRM模型的詳細幾何參數(shù)見文獻[23]。

圖8 CRM-WBT0模型外形

CRM模型先后在Langley研究中心的NTF風洞、Ames研究中心的11 ft (1 ft=0.305 m)跨聲速風洞以及歐洲跨聲速風洞ETW開展過風洞試驗。其中，CRM模型在ETW開展的風洞試驗得到了歐洲戰(zhàn)略風洞研究潛能改進項目(European Strategic Wind Tunnels Improved Research Potential,ESWIRP)的資助，獲得了豐富的試驗數(shù)據(jù)。因此，采用CRM-WBT0模型在互聯(lián)網(wǎng)上公開發(fā)布的ETW風洞試驗數(shù)據(jù)開展分離方法的測試與驗證。

表2給出了測試使用的3組試驗數(shù)據(jù)條件，可以利用車次號在ETW發(fā)布的CRM模型試驗數(shù)據(jù)中獲取對應的氣動力系數(shù)。

表2 試驗參數(shù)

此外，為了對氣動力系數(shù)分離效果進行評估，還采用了1組與表2中第2組具有不同載荷因子的試驗數(shù)據(jù)，用于修正得到剛性外形的氣動力系數(shù)數(shù)據(jù)。這組試驗數(shù)據(jù)的載荷因子q/E為5.058×10-7，其他參數(shù)條件與第2組相同，對應的ETW試驗的車次號為213。

4.2 氣動力系數(shù)分離結果

圖9給出了CRM-WBT0模型的氣動力系數(shù)分離結果。圖中、ETW Rigid表示利用相同雷諾數(shù)不同載荷因子下試驗數(shù)據(jù)外插得到的剛性外形氣動力數(shù)據(jù);ETW Elastic為ETW發(fā)布的風洞試驗數(shù)據(jù);Modeling為采用分離方法預測的剛性外形氣動力系數(shù)，通過對比Modeling表征的氣動力預測數(shù)據(jù)與ETW Rigid表征的修正剛性外形氣動力數(shù)據(jù)，評估預測結果的準確性。具體的插值思路是：認為在載荷因子q/E為零時，機翼不受氣動載荷的作用，即機翼沒有變形，利用不同載荷因子數(shù)據(jù)外插得到q/E等于零的試驗數(shù)據(jù)，即對應剛性外形下的氣動力。對于其他雷諾數(shù)，前文分析了雷諾數(shù)對靜氣動彈性變形引起的氣動力系數(shù)變化量影響很小，可以利用給定雷諾數(shù)下的彈性與剛性外形氣動力系數(shù)差量修正得到其他雷諾數(shù)下的剛性外形氣動力數(shù)據(jù)。這種基于真實的風洞試驗數(shù)據(jù)修正得到剛性外形氣動力系數(shù)，用于評估與驗證分離方法的結果具有較高的可信度。

可以看出，在模型迎角α≤4.0°時，采用分離方法得到的剛性外形氣動力系數(shù)與試驗修正數(shù)據(jù)吻合較好，僅阻力系數(shù)預測趨勢存在小量偏差。原因主要有2個方面。一是模型的迎角較小，機翼受載荷作用后產(chǎn)生的結構靜彈性變形對阻力系數(shù)的影響非常小，分離方法預測偏差超出了機翼結構靜彈性變形對阻力系數(shù)的影響量，圖9(a)中ETW Rigid和ETW Elastic 2組阻力系數(shù)之間的最大差量僅為1.07×10-4。二是分離方法建立在3個前提假設之上，對不滿足這些條件的預測結果存在一定的偏差，如圖9(b)～圖9(d)所示。此外，在迎角α=5.0°時，流動出現(xiàn)分離，原有的理論假設不再滿足，升力系數(shù)CL與俯仰力矩系數(shù)Cm分離結果與ETW Rigid之間存在明顯的偏差。

圖9 CRM-WBT0模型氣動力系數(shù)分離結果

4.3 分離方法的分析與討論

在3個假設完全成立的前提下，建立的分離方法有著嚴格的數(shù)學推導，能夠從含有模型靜氣動彈性變形試驗數(shù)據(jù)中分離得到理想的剛性外形氣動力結果。假設1與假設2有著理論上的依據(jù)，并不依賴于模型外形與繞流條件。但假設3的滿足條件與繞流條件相關，并不是完全成立的。如圖9(e)中，此時迎角較大流動出現(xiàn)了分離現(xiàn)象，升力系數(shù)CL與俯仰力矩系數(shù)Cm隨雷諾數(shù)的對數(shù)并不是線性變化的，從而使得待定系數(shù)k1、k2、k3的計算出現(xiàn)不準確，容易造成分離方法獲得的結果與剛性試驗數(shù)據(jù)修正結果之間出現(xiàn)較大的偏離。

式(3)中，Cf,rigid隨雷諾數(shù)變化的函數(shù)關系作為系數(shù)k2的定解條件，當Cf,rigid隨lg(Re)滿足線性變化時，利用式(5)能夠得到合理的k2值，從而得到合理的剛性外形氣動力系數(shù)分離結果，如圖9(a)～圖9(d)所示；而當Cf,rigid隨lg(Re)不滿足線性變化時，采用式(5)得到分離結果偏差較大，如圖9(e)所示。

從上述分析可以看出，依據(jù)式(5)的分離方法能夠對Cf,rigid隨lg(Re)呈現(xiàn)強線性變化的情況得到很好的分離結果。而對于其他情況，可以將分離方法進行變換，由式(3)可以得到：

Cf,rigid(Re)=Cf,elastic(Re,λ)-k2λ

(7)

式(7)中，首先利用式(5)得到一個k2的初值，然后在初值附近變換不同的k2值，得到不同的Cf,rigid隨Re的變化數(shù)據(jù)，并以此分析Cf,rigid隨Re的變化規(guī)律。此外，該分離方法也可以與其他技術一起，如CFD/CSD耦合數(shù)值模擬，共同評估靜氣動彈性變形影響的系數(shù)k2，然后分離得到無彈性變形影響的雷諾數(shù)效應試驗數(shù)據(jù)。

5 結 論

利用靜氣動彈性耦合數(shù)值模擬方法研究了雷諾數(shù)與靜氣動彈性變形對氣動力系數(shù)的影響規(guī)律，基于線性理論提出了3個假設，并發(fā)展了一種雷諾數(shù)/靜氣動彈性變形效應分離方法，通過大展弦比風洞模型的試驗數(shù)據(jù)對分離方法的準確性進行了檢驗。結果表明提出的雷諾數(shù)/靜氣動彈性變形效應分離方法能夠為常規(guī)風洞的雷諾數(shù)效應試驗提供一種簡便快捷的數(shù)據(jù)分析技術，預估雷諾數(shù)與靜氣動彈性變形對氣動力特性的影響，且在剛性模型氣動力系數(shù)隨雷諾數(shù)對數(shù)呈現(xiàn)近似線性變化的情況下能夠得到合理的氣動力系數(shù)分離結果。