風洞試驗雷諾數/靜氣動彈性效應分離方法

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

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

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

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

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

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

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

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

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

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

1 數值計算方法

1.1 NNW-FSI軟件

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

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

1.2 計算準確性驗證

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

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

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

圖2 HIRENASD模型流場計算網格

圖3給出了HIRENASD機翼彎曲變形計算與試驗結果[21]的對比曲線。其中,Ma、Re和α分別表示馬赫數、雷諾數和迎角,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 雷諾數/靜氣動彈性效應

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

2.1 有限元模型與計算條件

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

圖4 HIRENASD機翼有限元模型

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

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

表1 HIRENASD機翼計算條件

2.2 雷諾數/靜氣動彈性效應

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

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

圖5 HIRENASD機翼氣動力系數變化量

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

2.3 氣動力系數隨參數的變化規律

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

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

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

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

3 分離方法

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

Cf,elastic=f(Re,q)

(1)

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

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

(2)

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

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

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

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

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

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

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

(3)

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

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

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

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

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

(4)

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

(5)

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

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

(6)

4 分離方法測試

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

4.1 測試對象及數據

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

圖8 CRM-WBT0模型外形

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

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

表2 試驗參數

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

4.2 氣動力系數分離結果

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

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

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

4.3 分離方法的分析與討論

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

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

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

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

(7)

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

5 結 論

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