孫 巖, 黃 勇, 王運濤,*, 孟德虹, 王 昊

(1. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室, 四川 綿陽 621000; 2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000)

TRIP軟件的靜氣動彈性計算模塊開發(fā)及精度驗證

孫 巖1,2, 黃 勇2, 王運濤2,*, 孟德虹2, 王 昊2

(1. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室, 四川 綿陽 621000; 2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000)

在結(jié)構(gòu)網(wǎng)格流場求解軟件TRIP基礎(chǔ)上，發(fā)展了一套靜氣動彈性計算模塊。首先簡要介紹了靜氣動彈性計算模塊總體架構(gòu)、構(gòu)成單元和耦合策略，然后詳細介紹了構(gòu)成單元中采用的一些數(shù)值計算方法。最后利用大展弦比機翼、DLR-F6翼身組合體模型和HIRENASD機翼模型對靜氣動彈性計算模塊進行了測試和精度驗證。大展弦比機翼計算結(jié)果表明柔度矩陣、模態(tài)疊加和有限元數(shù)值求解在模態(tài)選取合適的情況下能夠得到一致的變形計算結(jié)果。DLR-F6和HIRENASD模型不同計算軟件的結(jié)果一致性表明本文靜氣動彈性計算模塊采用了正確的算法流程，而計算與試驗的結(jié)果一致性表明靜氣動彈性計算模塊具有很好的預(yù)測精度。

TRIP；靜氣動彈性；柔度矩陣；弱耦合；HIRENASD；DLR-F6

0 引 言

為了減輕結(jié)構(gòu)重量、提高飛行機動性，飛行器設(shè)計中傾向采用柔性的結(jié)構(gòu)與材料，使得飛行器氣動彈性效應(yīng)日益明顯，需要更加精確的氣動彈性數(shù)值模擬工具用于飛行器的精細氣動設(shè)計[1]。

此外，風(fēng)洞模型的完全剛性假設(shè)在新型先進飛機設(shè)計中受到了挑戰(zhàn)。越來越多的研究表明機翼在風(fēng)洞氣動載荷作用下的變形對模型氣動特性存在明顯的影響，而這種影響在增壓試驗中表現(xiàn)的更為顯著[2]。例如波音的超聲速運輸機模型在美國蘭利的NTF風(fēng)洞開展高雷諾數(shù)風(fēng)洞試驗，翼梢扭轉(zhuǎn)角可以達到-2.2°，由此造成不同速壓下的巡航阻力系數(shù)差異最大接近50個阻力單位[3]。風(fēng)洞模型靜彈性變形對氣動特性的影響已經(jīng)引起了全世界研究者的關(guān)注與重視[4]，高效高精度的靜氣動彈性計算程序逐漸成為研究這一問題的重要手段。

為了提高計算結(jié)果的可信度水平，國外廣泛開展了氣動彈性計算軟件的驗證與確認研究，其中影響范圍較大的是AIAA發(fā)起的氣動彈性預(yù)測研討會AePW(Aeroelastic Prediction Workshop)[5]，得到了世界范圍內(nèi)大量學(xué)者的積極響應(yīng)。

近些年，隨著計算機硬件技術(shù)的快速發(fā)展，國內(nèi)高校與科研機構(gòu)基于各自求解器開發(fā)了多種氣動彈性數(shù)值計算平臺[6-8]，并廣泛應(yīng)用于解決工程實際問題，但靜氣動彈性數(shù)值計算程序精度驗證方面的系統(tǒng)研究還比較少。

本文在課題組開發(fā)的流場計算軟件TRIP上發(fā)展了一套高效的靜氣動彈性計算模塊，下文將簡要介紹計算模塊的架構(gòu)、組成、流程和采用的數(shù)值方法。然后通過2個算例對計算模塊的精度進行驗證。

1 靜氣動彈性計算模塊架構(gòu)

圖1給出了TRIP軟件靜氣動彈性計算模塊的總體架構(gòu)，該模塊包含四個主要功能單元：氣動數(shù)值計算、結(jié)構(gòu)靜變形計算、耦合界面數(shù)據(jù)插值以及動態(tài)網(wǎng)格方法，四個主要功能單元采用弱耦合方式串聯(lián)在一起，如圖1所示。

圖1的左側(cè)給出了氣動載荷和結(jié)構(gòu)位移的傳遞路線。首先CFD計算出氣動載荷，利用界面數(shù)據(jù)插值將氣動載荷傳遞給CSM，CSM利用氣動輸入載荷計算得到結(jié)構(gòu)的位移，通過界面數(shù)據(jù)插值再將結(jié)構(gòu)位移傳遞給流場，并利用動態(tài)網(wǎng)格技術(shù)生成流場物面變形后的新網(wǎng)格，CFD再計算新網(wǎng)格下的流場。如此反復(fù)迭代，直至載荷和位移均達到收斂。圖1的右側(cè)給出了每個功能單元采用的主要計算方法，下節(jié)將對這些計算方法進行更加詳細的介紹。

2 數(shù)值計算方法

2.1流場計算方法

流場解算器采用課題組開發(fā)多年的亞跨超流場解算軟件TRIP。TRIP軟件通過有限體積法離散求解歐拉(Euler)和雷諾平均方程(RANS)，方程對流項離散采用二階MUSCL型的Roe通量差分分裂格式，粘性項離散采用二階中心差分格式，湍流模型采用工程常用的單方程SA模型[9]與兩方程SST模型[10]，計算加速方面采用多重網(wǎng)格和基于MPI信息傳遞的并行計算。課題組針對TRIP已經(jīng)開展了大量的算例驗證與測試，軟件的精度和可信度得到了廣泛的認可[11-12]。

2.2結(jié)構(gòu)靜變形計算方法

結(jié)構(gòu)靜變形計算采用柔度矩陣、模態(tài)疊加和有限元數(shù)值求解三種方法。

柔度矩陣方法中，結(jié)構(gòu)網(wǎng)格點的位移通過下式計算得到:

式中，us為結(jié)構(gòu)網(wǎng)格點的位移矩陣，C為結(jié)構(gòu)的柔度矩陣，fs為結(jié)構(gòu)網(wǎng)格點上的輸入載荷。本文中結(jié)構(gòu)柔度矩陣C通過有限元模型施加單位載荷的位移響應(yīng)組裝得到。

模態(tài)疊加法中，結(jié)構(gòu)網(wǎng)格點的位移通過下式計算得到[13]:

式中，us含義同式(1)，Φs為結(jié)構(gòu)的廣義模態(tài)矩陣，通過結(jié)構(gòu)模態(tài)分析得到，q為系統(tǒng)的廣義位移矢量。對于結(jié)構(gòu)靜變形有:

式中，K為結(jié)構(gòu)的剛度矩陣，fs含義同式(1)，式(2)代入式(3)，可以得到:

式中:

式中，Ω為結(jié)構(gòu)模態(tài)剛度矩陣，是一個對角陣，對角線各元素為廣義模態(tài)圓頻率的平方。式(4)代入式(2)即可求解得到結(jié)構(gòu)網(wǎng)格點的位移:

有限元數(shù)值求解采用了自主編寫的一套結(jié)構(gòu)靜變形求解程序，在方法上借鑒了商業(yè)軟件Nastran的數(shù)據(jù)格式和眾多開源程序的設(shè)計架構(gòu)。

2.3耦合界面數(shù)據(jù)傳遞

氣動/結(jié)構(gòu)之間的數(shù)據(jù)傳遞可以通過一個矩陣進行描述，比如對于位移變量可以表示為：

式中，u為網(wǎng)格點位移矢量，下標a表示氣動，下標s表示結(jié)構(gòu)，H為轉(zhuǎn)換矩陣。為了滿足傳遞過程中的能量守恒，有虛功原理可以得到：

式中，f為網(wǎng)格點上的載荷。

本文靜氣動彈性計算模塊采用RBF、TPS和IPS三種方法構(gòu)建傳遞矩陣H。一般形式的徑向基函數(shù)插值可以表示為[14]:

式中，f(x)為被插函數(shù)，x為插值點位置坐標，xi為插值基點的位置坐標，xj(j=1,…,m)為位置坐標的分量，‖x-xi‖表示兩點之間的距離，N為插值基點的數(shù)量，m為描述問題的維數(shù)，二維m為2，三維m為3，βi(i=0,…,m)、γi(i=1,…,N)為插值系數(shù)，φ為徑向基函數(shù)。

利用RBF方法得到的傳遞矩陣H可以表示為：

其中:

式(11)中的M為氣動網(wǎng)格點的數(shù)目，式(12)中右側(cè)的0為(m+1)×N的零矩陣，I為N×N的單位矩陣。

不同的基函數(shù)φ可以得到不同的傳遞矩陣H，TPS和IPS是RBF插值的一種特殊形式，對應(yīng)不同的基函數(shù)類型，基函數(shù)的選取可以參考文獻[14]。

2.4動態(tài)網(wǎng)格生成

本文動態(tài)網(wǎng)格生成采用網(wǎng)格變形的思路，可以避免引起新的數(shù)值離散誤差。目前計算模塊主要采用面向結(jié)構(gòu)網(wǎng)格的RBF+TFI方法，針對一些特殊問題，也采用RBF+Delaunay和RBF+HBGM的網(wǎng)格變形方法，有關(guān)這些方法可以參考文獻[15-16]。

3 靜氣動彈性計算模塊精度驗證

本節(jié)將首先通過一個大展弦比機翼的變形結(jié)果分析三種不同結(jié)構(gòu)靜變形計算方法的差異，然后利用兩個風(fēng)洞模型的試驗數(shù)據(jù)結(jié)果對靜氣動彈性計算模塊精度進行驗證。下文如不作特殊說明，流場計算中湍流模擬采用兩方程SST模型，耦合數(shù)據(jù)傳遞采用TPS方法，動態(tài)網(wǎng)格變形采用RBF+TFI方法。

3.1大展弦比機翼變形結(jié)果分析

本小節(jié)將分析三種不同結(jié)構(gòu)靜變形求解方法結(jié)果的差異，圖2給出了采用的大展弦比機翼的外形、網(wǎng)格和有限元模型。其中有限元模型采用不同厚度的殼單元構(gòu)建，因此耦合數(shù)據(jù)傳遞采用IPS方法。

圖3給出了三種不同方法計算得到的大展弦比機翼變形結(jié)果，圖中Moden表示采用結(jié)構(gòu)前n階模態(tài)位移的變形計算結(jié)果。

圖3可以看出：柔度矩陣方法與有限元分析結(jié)構(gòu)一致性很好，模態(tài)疊加方法的變形結(jié)果與選取的模態(tài)有關(guān)。但隨著模態(tài)數(shù)量的增加，模態(tài)方法的結(jié)果與前兩種方法的結(jié)果也很快趨向一致。另一方面，不同模態(tài)對彎曲變形和扭轉(zhuǎn)變形的貢獻不同，對于彎曲變形，前2階的結(jié)構(gòu)模態(tài)貢獻了主要的變形，而對于扭轉(zhuǎn)變形，后續(xù)模態(tài)仍然有著明顯的貢獻。因此，在選擇結(jié)構(gòu)模態(tài)計算結(jié)構(gòu)靜變形時，應(yīng)當首先開展結(jié)構(gòu)模態(tài)對靜變形計算結(jié)果的影響分析。如不作特殊說明，本文后續(xù)的算例中均采用柔度矩陣方法計算結(jié)構(gòu)靜變形。

3.2DLR-F6模型靜氣動彈性模擬

DLR-F6是德國宇航院設(shè)計的一類現(xiàn)代民用運輸機標模，在多座風(fēng)洞開展了不同類型的風(fēng)洞試驗，具有豐富的試驗數(shù)據(jù)，有關(guān)DLR-F6模型的詳細可以參考文獻[17]。本文采用DLR-F6的翼身組合體構(gòu)型開展精度驗證，圖4給出了該構(gòu)型的外形、網(wǎng)格拓撲和有限元模型。

流場計算Ma數(shù)為0.75，基于平均氣動弦長的雷諾數(shù)Re為3.0×106，模型迎角為-3°～1.5°，來流速壓q與模型材料彈性模量E的比值q/E為2.0×10-7。圖5給出了該計算參數(shù)條件下DLR-F6翼身組合體模型機翼變形展向分布結(jié)果。圖中Test表示在NASA的NTF風(fēng)洞利用光學(xué)方法測量的機翼變形，TAU為DLR的流場解算程序。

可以看出：不同升力系數(shù)下，試驗Test、TAU和本文Present得到的彎曲變形結(jié)果在變形分布和量值上都比較一致，表明本文靜氣動彈性計算模塊采用了正確的算法流程且具有較好的彎曲變形預(yù)測精度。由于DLR-F6是實體金屬模型，來流的速壓并不大，因此由氣動載荷引起的機翼最大彎曲變形并不明顯，約為6.5 mm(半翼展為585.6 mm)。

3.3HIRENASD模型靜氣動彈性模擬

HIRENASD是亞琛工大設(shè)計的一類典型運輸機機翼外形，在歐洲跨聲速風(fēng)洞ETW開展了大量試驗，有關(guān)HIRENASD模型的詳細參數(shù)可以參考文獻[18]。為了消除壁面邊界層對機翼流場的影響，風(fēng)洞中給HIRENASD機翼模型安裝了一個假機身，本文采用包含假機身的外形開展靜氣動彈性模擬。圖6給出了HIRENASD機翼模型的外形、網(wǎng)格拓撲和有限元模型。其中流場計算網(wǎng)格從氣動彈性預(yù)測會議AePW的官網(wǎng)上下載得到，為ICEM軟件生成的多塊對接網(wǎng)格，如圖6(a)；有限元的固支約束添加在機翼與天平連接的端面上，如圖6(b)。

流場計算的Ma數(shù)為0.80，基于平均氣動弦長的雷諾數(shù)Re為1.4×107，模型迎角α為3°，來流速壓q與模型材料彈性模量E的比值q/E為4.7×10-7。

圖7給出了該參數(shù)條件下的機翼彎曲變形展向分布結(jié)果，可以看出，試驗Test、TAU和本文Present得到的機翼變形結(jié)果十分吻合。

圖8給出了HIRENASD機翼展向相對位置η等于0.95時的翼剖面壓力系數(shù)分布，可以看出，采用剛性外形，計算的壓力系數(shù)分布與試驗結(jié)果之間存在顯著的差異，而采用靜氣動彈性耦合計算得到的壓力系數(shù)分布與試驗結(jié)果十分吻合，進一步驗證了本文靜氣動彈性計算模塊的精度。

4 結(jié) 論

本文在TRIP軟件基礎(chǔ)發(fā)展了靜氣動彈性計算模塊，并通過算例對計算結(jié)果的精度進行了驗證，說明該計算模塊具有較好的靜氣動彈性預(yù)測精度。同時，計算結(jié)果表明風(fēng)洞模型結(jié)構(gòu)靜變形對壓力系數(shù)分布有著明顯影響，應(yīng)當?shù)玫皆囼灁?shù)據(jù)使用單位的重視。

致謝:本文得到了課題組張玉倫、洪俊武、王光學(xué)、張書俊、李偉和楊小川等人的幫助，在此表示感謝。

[1]Mouton S, Sant Y, Lyonnet M. Predication of the aerodynamic effect of model deformation during transonic wind tunnel tests[J]. International Journal of Engineering Systems Modelling and Simulation, 2013, 5(1-3): 44-56.

[2]孫巖, 張征宇, 鄧小剛, 等. 風(fēng)洞模型靜彈性變形對氣動力影響研究[J]. 空氣動力學(xué)學(xué)報, 2013, 31(3): 294-300.

[3]Owens L R, Wahls R A, Rivers S M. Off-design Reynolds number effects for a supersonic transport[J]. Journal of Aricraft, 2005, 42(6): 1427- 441.

[4]Levy D W, Laflin K R, Tinoco E N, et al. Summary of data from the fifth computational fluid dynamics drag prediction workshop[J]. Journal of Aircraft, 2014, 51(4): 1194-1213.

[5]Heeg J, Chwalowski P, Florance J P, et al. Overview of the aeroelastic predication workshop[R]. AIAA 2013-0783, 2013.

[6]楊超, 王立波, 謝長川, 等. 大變形飛機配平與飛行載荷分析方法[J]. 中國科學(xué)(E輯: 技術(shù)科學(xué)), 2012, 42(10): 1137-1147.

[7]徐敏, 邱滋華, 裴曦. 基于CFD/CSD/CAA耦合的聲氣動彈性仿真技術(shù)在壁板顫振中的應(yīng)用[J]. 強度與環(huán)境, 2013, 40(5): 16-24.

[8]黃煒, 陸志良, 唐迪, 等. 基于多點約束的大展弦比機翼靜氣動彈性計算[J]. 北京航空航天大學(xué)學(xué)報, 2014, 40(12): 1666-1671.

[9]Spalart P R, Allmaras S R. A one-equation turbulence model for aerodynamic flows[R]. AIAA-92-0439, 1992.

[10]Menter F R. Two-equation eddy-viscosity turbulence models for engineering application[J]. AIAA Journal, 1994, 32(8): 1598-1605.

[11]王運濤, 王光學(xué), 張玉倫. DPWⅢ機翼和翼身組合體構(gòu)型數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報, 2011, 29(3): 264-269.

[12]王運濤, 張書俊, 孟德虹. DPW4翼/身/平尾組合體的數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報, 2013, 31(6): 739-744.

[13]Ritter M. Static and forced motion aeroelastic simulations of the HIRENASD wind tunnel model[R]. AIAA 2012-1633, 2012.

[14]Rendall T C S, Allen C B. Unified fluid-structure interpolation and mesh motion using radial basis functions[J]. International Journal for Numerical Methods in Engineering, 2008, 74: 1519-1559.

[15]孫巖, 鄧小剛, 王運濤, 等. RBF_TFI結(jié)構(gòu)動網(wǎng)格技術(shù)在風(fēng)洞靜氣動彈性修正中的應(yīng)用[J]. 工程力學(xué), 2014, 31(10): 228-233.

[16]孫巖, 鄧小剛, 王光學(xué), 等. 基于徑向基函數(shù)改進的Delaunay圖映射動網(wǎng)格方法[J]. 航空學(xué)報, 2014, 35(3): 727-735.

[17]Burner A W, Goad W K, Massey E A, et al. Wind deformation measurements of the DLR-F6 transport configuration in the national transonic facility[R]. AIAA 2008-6921, 2008.

[18]Ritter M, Hassan D. Assesment of the ONERA/DLR numerical aeroelastics predication capabilities on the HIRENASD configuration[R]. IFASD-2011-109, 2011.

DevelopmentandprecisionvalidationofstaticaeroelasticcomputationalmoduleonflowsolverTRIP

(1.StateKeyLaboratoryofAerodynamics,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China; 2.ComputationalAerodynamicsInstituteofChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China)

A static aeroelastic computational module is developed on the flow solver TRIP. Firstly, the frame, main component units and coupling strategy of the static aeroelastic computational module are introduced briefly. Then numerical methods used in the component units are described in detail. Finally, the static aeroelastic computational module is tested and validated through three cases: high respect-ratio wing, DLR-F6 wing-body model and HIRENASD wing model. The computational results of high respect-ratio wing show that flexible matrix method, modal superposition method and finite element method can have same deformation results when the structure modes are selected properly. The result consistency of DLR-F6 or HIRENASD model among different solvers demonstrates that the static aeroelastic computational module has adopted a correct algorithm procedure. The results consistency of DLR-F6 or HIRENASD model between computations and tests shows that the present static aeroelastic computational module has a good prediction precision.

TRIP; static aeroelasticity; flexibility matrix; weak coupling; HIRENASD; DLR-F6

V211.3

A

10.7638/kqdlxxb-2015.0154

0258-1825(2017)05-0620-05

2015-08-26;

2015-10-14

國家重點研發(fā)計劃(2016YFB0200700)

孫巖(1986-)，男，安徽鳳陽人，博士，助理研究員，研究方向：計算氣動彈性力學(xué). E-mail：supersunyan@163.com

王運濤*(1967-)，男，博士，研究員，博士生導(dǎo)師，主要研究方向：計算空氣動力學(xué). E-mail: ytwang@skla.cardc.cn

孫巖, 黃勇, 王運濤, 等. TRIP軟件的靜氣動彈性計算模塊開發(fā)及精度驗證[J]. 空氣動力學(xué)學(xué)報, 2017, 35(5): 620-624.

10.7638/kqdlxxb-2015.0154 SUN Y, HUANG Y, WANG Y T, et al. Development and precision validation of static aeroelastic computational module on flow solver TRIP[J]. Acta Aerodynamica Sinica, 2017, 35(5): 620-624.

SUN Yan1,2, HUANG Yong2, WANG Yuntao2,*, MENG Dehong2, WANG Hao2