幺虹， 郭承鵬， 張穎

(中國航空工業(yè)空氣動力研究院 高速高雷諾數(shù)氣動力航空科技重點實驗室， 遼寧 沈陽 110034)

?

飛行器操縱面嗡鳴的非結(jié)構(gòu)網(wǎng)格并行計算方法

幺虹， 郭承鵬， 張穎

(中國航空工業(yè)空氣動力研究院 高速高雷諾數(shù)氣動力航空科技重點實驗室， 遼寧 沈陽 110034)

為了數(shù)值模擬飛行器操縱面的嗡鳴現(xiàn)象，在集群計算機MPI并行計算環(huán)境下建立基于三維非定常歐拉方程耦合結(jié)構(gòu)運動方程的嗡鳴計算方法.氣動流場求解采用基于非結(jié)構(gòu)網(wǎng)格的中心有限體積法進行空間離散，時間推進采用雙時間方法，結(jié)構(gòu)運動方程采用Adams預估校正方法求解.針對翼面與操縱面縫隙間存在的網(wǎng)格運動問題，在非結(jié)構(gòu)網(wǎng)格系統(tǒng)上采用Delaunay圖映射方法實現(xiàn)網(wǎng)格的運動變形.最后，使用飛行器操縱面標準嗡鳴計算模型對計算方法進行驗證，結(jié)果表明:所建立的并行計算方法正確，程序具有很好的計算效率，能夠?qū)︼w行器操縱面嗡鳴進行高效的數(shù)值分析.

嗡鳴； 飛行器； 操縱面； 并行計算； 動網(wǎng)格； Adams預估校正

飛行器操縱面嗡鳴是飛行器在跨音速飛行階段發(fā)生的一種操縱面單自由度振動[1]，其本質(zhì)特征與多模態(tài)耦合的經(jīng)典彎/扭型顫振是完全不同的.Parker等[2]對美國航空航天局(NASA)的三維機翼操縱面標準模型(NASP)開展了嗡鳴響應(yīng)特性的分析研究，發(fā)現(xiàn)操縱面上激波的運動和操縱面振動之間的相位差造成了嗡鳴現(xiàn)象.David[3]建立了操縱面偏轉(zhuǎn)狀態(tài)和鉸鏈力矩之間的數(shù)學模型.Oddvar[4]使用Euler方程數(shù)值計算方法研究了跨音速副翼嗡鳴現(xiàn)象.劉千剛等[5]使用CFD手段開展了跨音速操縱面嗡鳴Hopf分叉分析及結(jié)構(gòu)參數(shù)對嗡鳴特性影響的研究.史愛明等[6]基于Euler 方程，使用單自由度方法對機翼嗡鳴進行了數(shù)值模擬研究.張偉偉等[7]基于Euler 方程分析了二維翼型模型的B 型和C 型嗡鳴特性.楊國偉等[8]基于多塊結(jié)構(gòu)網(wǎng)格研究了飛機副翼嗡鳴現(xiàn)象.對于飛行器操縱面嗡鳴現(xiàn)象，激波及流動分離在其中起著非常重要的作用.由于非定常氣動力和激波與操縱面相互干擾的特點，數(shù)值模擬計算量非常巨大，所以準確、高效的數(shù)值預測飛行器操縱面嗡鳴，對計算方法有較高的要求.本文研究在集群MPI并行計算環(huán)境下的飛行器操縱面嗡鳴計算方法和程序，采用MPI進行并行化處理以進一步提高計算效率，最后使用飛行器嗡鳴標準計算模型對計算方法和程序進行驗證計算.

1 計算方法

1.1 非定常氣動力計算方法

非定常氣動力計算的控制方程采用三維可壓縮非定常積分形式的歐拉方程，其守恒形式為

(1)

對方程(1)在空間方向上采用有限體積法進行離散，在時間方向上采用二階精度的雙時間方法進行離散，在偽時間上采用龍格-庫塔方法進行推進.

1.2 操縱面結(jié)構(gòu)運動方程求解

飛行器操縱面的運動可以用單自由度的結(jié)構(gòu)運動方程描述[3]，忽略阻尼項，其具體形式為

(2)

1.3 非結(jié)構(gòu)動網(wǎng)格方法

當前針對含運動邊界的網(wǎng)格方法主要有網(wǎng)格變形方法、網(wǎng)格重構(gòu)法[9]、浸入邊界法[10-11]和嵌套網(wǎng)格法[12].除網(wǎng)格變形方法外，其他的動網(wǎng)格方法均要對流場進行插值，因而在非定常計算過程中會損失精度，文中采用不增加或減少網(wǎng)格節(jié)點并保持原網(wǎng)格拓撲結(jié)構(gòu)的網(wǎng)格變形方法.近年來研究較多的徑向基函數(shù)(RBF)動網(wǎng)格方法是Boer等[13]于2007年首次提出，也是一種需要全局信息的動網(wǎng)格方法，其計算精度高，但效率低.

考慮到翼面-操縱面縫隙間網(wǎng)格運動效率與魯棒性要求，采用非結(jié)構(gòu)網(wǎng)格進行計算，使用Delaunay方法實現(xiàn)網(wǎng)格的運動變形，具體有如下4個主要步驟.

步驟1 收集所有邊界點，生成Delaunay背景網(wǎng)格圖.

步驟2 建立計算網(wǎng)格節(jié)點和背景網(wǎng)格間的對應(yīng)關(guān)系，計算插值權(quán)重系數(shù).

步驟3 背景網(wǎng)格的運動和變形.

考慮到在整個計算過程中計算網(wǎng)格節(jié)點和背景網(wǎng)格單元之間的對應(yīng)關(guān)系保持不變，即背景網(wǎng)格只需生成一次，同時為了便于實現(xiàn)網(wǎng)格變形模塊的并行化，因此在程序設(shè)計上將前兩步放到前處理程序進行處理；然后，將處理好的背景網(wǎng)格相關(guān)信息放到各個分區(qū)網(wǎng)格中，減少了實際并行計算過程中各線程之間的數(shù)據(jù)傳遞量.

(a) 運動前 (b) 運動后圖1 非結(jié)構(gòu)網(wǎng)格運動效果圖Fig.1 Schematic motion effect of unstructured grids

使用文中開發(fā)出的基于Delaunay背景網(wǎng)格動網(wǎng)格方法對NACA0012翼型進行剛體旋轉(zhuǎn)，動網(wǎng)格效果如圖1所示.通過實際測試表明，使用基于Delaunay背景網(wǎng)格的動網(wǎng)格方法能夠快速實現(xiàn)網(wǎng)格運動.

2 嗡鳴程序并行化設(shè)計

在中國航空工業(yè)空氣動力研究院研發(fā)的非結(jié)構(gòu)混合網(wǎng)格流場計算軟件UNSMB的基礎(chǔ)上進行嗡鳴程序設(shè)計.圖2為嗡鳴程序的整體計算流程.在圖2中，除了非定常流場并行求解以外，還需對操縱面的氣動力矩、CSD計算模塊和更新網(wǎng)格幾何數(shù)據(jù)進行并行化處理.

圖2 嗡鳴程序流程Fig.2 Flow chart of buzz program

操縱面氣動力矩計算并行化設(shè)計的思路是，在每個物理時間開始之前(n-1時刻)，各進程分別計算本分區(qū)的操縱面氣動力矩，使用MPI_Allgather將各區(qū)計算的氣動力矩收集起來，在主進程中求和得到整個操縱面的氣動力矩.在主進程進行結(jié)構(gòu)運動方程的求解，可得到操縱面偏轉(zhuǎn)角β，并輸出該時刻的時間值及偏轉(zhuǎn)角.然后，使用MPI_Bcast將偏轉(zhuǎn)角β值發(fā)送到各個進程中，而各進程根據(jù)β值分別對該分區(qū)操縱面(僅操縱面物面表面網(wǎng)格)進行繞空間軸線旋轉(zhuǎn)操作.使用MPI_Send和MPI_Recv將各進程上各分區(qū)的操縱面網(wǎng)格坐標收集到主進程當中，更新物面網(wǎng)格.

由物面網(wǎng)格、對稱面網(wǎng)格、遠場網(wǎng)格和足夠大的一個背景網(wǎng)格區(qū)域，使用Delaunay方法生成背景網(wǎng)格.空間網(wǎng)格點就被Delaunay背景網(wǎng)格分區(qū)了，可以使用4個體積比表達的權(quán)值來表示空間網(wǎng)格點在Delaunay背景網(wǎng)格中的具體位置.物面邊界更新之后，Delaunay背景網(wǎng)格各點的位置更新，空間網(wǎng)格隨之更新.在此過程中，可以將Delaunay背景網(wǎng)格生成放到前處理部分.如此各分區(qū)更新操縱面表面網(wǎng)格之后，在本區(qū)就可以完成本區(qū)空間網(wǎng)格的更新，然后各進程分別對本分區(qū)的新網(wǎng)格點重新計算控制體幾何信息、網(wǎng)格速度等.總體流程有如下6個主要步驟.

步驟1 前處理.使用Metis進行網(wǎng)格分區(qū)，Delaunay背景網(wǎng)格生成及網(wǎng)格點Delaunay背景圖映射關(guān)系建立及權(quán)值的計算.

步驟2 氣動力矩計算.每個物理時間步開始之前并行化計算氣動力矩.

步驟3 結(jié)構(gòu)運動方程求解.僅在主進程進行結(jié)構(gòu)運動方程求解，然后將求解得到的操縱面偏轉(zhuǎn)等信息發(fā)送到所有進程，并進行時間層的切換.

步驟4 網(wǎng)格運動.在各區(qū)上分別進行操縱面偏轉(zhuǎn)操作，操縱面網(wǎng)格發(fā)生改變，然后使用Delaunay圖映射法對空間網(wǎng)格進行運動.空間網(wǎng)格運動充分利用前處理分區(qū)信息及Delaunay圖映射的映射關(guān)系不變的特點，方便進行并行化.

步驟5 處理網(wǎng)格幾何信息.對更新后的網(wǎng)格計算法矢、面積、體積等幾何信息，計算網(wǎng)格速度.

步驟6 進行下一時間步的流場計算.

嗡鳴計算程序并行化處理的關(guān)鍵是背景網(wǎng)格的處理,它充分利用了Delaunay圖映射方法并行化程序設(shè)計方便、計算效率高的優(yōu)點.經(jīng)算例測試，該嗡鳴算例使用歐拉方程串行計算一個狀態(tài)推進0.5 s(時間步長0.000 1 s)約需7 h，而采用并行計算，其總時間約4 h，并行效率基本呈線性關(guān)系.

3 NASP標模嗡鳴計算

3.1 計算模型及網(wǎng)格

NASP系列機翼是用于進行跨音速操縱面嗡鳴實驗的一組不同平面形狀的全翼展操縱面機翼[2,14].計算所采用的模型為NASA三維操縱面標準模型NASP，如圖3所示.模型相關(guān)參數(shù)如下：操縱面質(zhì)量m=0.267 624 kg，操縱面弦長c=0.101 6 m，轉(zhuǎn)動頻率f=16.5 Hz，轉(zhuǎn)動慣量Iβ=88.17 Mg·m2，轉(zhuǎn)動剛度Kβ=9.484 3 N·m·rad-1，其他幾何參數(shù)詳見文獻[2].計算采用非結(jié)構(gòu)純四面體網(wǎng)格，整體計算域和計算網(wǎng)格示意圖，如圖4所示.圖4中：網(wǎng)格節(jié)點數(shù)約8萬個.

圖3 NASP嗡鳴計算標準模型 圖4 NASP模型表面及對稱面計算網(wǎng)格Fig.3 NASP buzz computation standard mode Fig.4 Computation grids on surface and symmetrical plane for NASP model

(a) M=0.98，Q=1.532 2 kPa

(b) M=0.96，Q=2.082 8 kPa (c) M=1.30，Q=1.819 4 kPa圖5 不同風洞試驗條件下操縱面偏角隨時間的變化歷程Fig.5 Change history deflection angle of control surface with time under experimental conditions of different wind tunnels

3.2 計算結(jié)果分析

對NASP標模M=0.98，Q=1.532 2 kPa工況的嗡鳴特性進行了計算.經(jīng)過測試研究表明，非定常計算的內(nèi)迭代步數(shù)和時間步長對計算結(jié)果(如嗡鳴的邊界)有極大的影響.為此，選取時間步長為5×10-5s，內(nèi)迭代步數(shù)為20步進行計算.共有80萬個迭代，計算使用集群計算機環(huán)境，共使用CPU核心24個，機時4 h.在不同風洞試驗條件下計算得到操縱面偏角隨著時間的變化曲線，如圖5所示.

從圖5(a)可知：在0.5 s之前，操縱面偏角呈發(fā)散狀態(tài)；到0.5 s之后，操縱面偏角在振幅9.3°做極限環(huán)振蕩，計算結(jié)果與試驗結(jié)果一致[2].假使操縱面轉(zhuǎn)軸強度足夠大，則操縱面在此狀態(tài)下將以振幅9.3°做極限環(huán)振蕩，振蕩頻率25.3 Hz，若結(jié)構(gòu)強度不夠，則在此狀態(tài)，操縱面結(jié)構(gòu)發(fā)生破壞.從圖5(b),(c)可知：當把M值降低到0.96或提高到1.3時，操縱面偏角隨時間變化呈收斂趨勢.計算結(jié)果表明，在M<0.96或M>1.3時操縱面不會發(fā)生嗡鳴(等幅的極限環(huán)振蕩)現(xiàn)象.

4 結(jié)束語

在中國航空工業(yè)空氣動力研究院UNSMB軟件的基礎(chǔ)上，基于集群計算機MPI并行計算環(huán)境，研究了飛行器操縱面嗡鳴計算方法和程序模塊，并對NASP操縱面嗡鳴標模進行了數(shù)值計算.結(jié)果表明，所提出的嗡鳴并行計算方法和程序模塊能夠高效、高精度地模擬操縱面嗡鳴現(xiàn)象，具有較高的工程應(yīng)用價值.

[1] 葉正寅，張偉偉，史愛明.流固耦合力學基礎(chǔ)及其應(yīng)用[M].哈爾濱:哈爾濱工業(yè)大學出版社，2010:1-294.

[2] PARKER E， SPAN C，SOISTMANN D. Aileron buzz investigated on several generic NASP wing configurations[C]∥32nd Structures, Structural Dynamics, and Materials Conference.Baltimore:AIAA,1991.DOI:10.2514/6.1991-936.

[3] NIXON D.An analytic model for control surface buzz[C]∥36th AIAA Aerospace Sciences Meeting and Exhibit.Reno:AIAA,1998.DOI: 10.2514/6.1998-417.

[4] ODDVAR O B.Nonclassical aileron buzz in transonic flow[C]∥34th Structures, Structural Dynamics and Materials Conference.La Jolla:AIAA,1993.DOI: 10.2514/6.1993-1479.

[5] 劉千剛，代捷，白俊強.跨音速操縱面嗡鳴Hopf 分叉分析及結(jié)構(gòu)參數(shù)對嗡鳴特性影響的研究[J].航空學報，1999，20(6):527-532.

[6] 史愛明，楊永年，葉正寅.跨音速單自由度非線性顫振: 嗡鳴的數(shù)值分析[J].西北工業(yè)大學學報，2004，22(4):525-528.

[7] 張偉偉，葉正寅，史愛明，等.基于Euler 方程的B 型和C 型嗡鳴特性數(shù)值研究[J].振動工程學報，2005，18(4):458-464.

[8] YANG Guowei，OBAYASHI S,NAKAMICHI J.Aileron buzz simulation using an implicit multiblock aeroelastic solver[J].Journal of Aircraft，2003，40(3):580-589.

[9] 周璇，李水鄉(xiāng)，孫樹立，等.非結(jié)構(gòu)網(wǎng)格變形方法研究進展[J].力學進展，2011，41(5):547-561.

[10] SUDHAKAR Y，VENGADESAN S.Flight force production by flapping insect wings in inclined stroke plane kinematics[J].Computers & Fluids，2010，39(4):683-695.

[11] SHYY W，AONO H，CHIMAKURTHI S K.Recent progress in flapping wing aerodynamics and aeroelasticity[J].Progress in Aerospace Sciences，2010，46(7):284-327.

[12] HEATHCOTE J.Flexible flapping airfoil propulsion at zero freestream velocity [J].AIAA Journal，2004，42(11):2196-2204.

[13] BOER A，SCHOOT M S，F(xiàn)ACULTY H B.Mesh deformation based on radial basis function interpolation[J].Computers and Structures，2007，85(11):784-795.

[14] SPAIN C，SOISTMANN D，PARKER E,et al.An overview of selected NASP aeroelastic studies at the NASA langley research center[C]∥2nd International Aerospace Planes Conference.Orlando:AIAA,1990.DOI:10.2514/6.1990-5218.

(責任編輯： 黃仲一 英文審校： 崔長彩)

Parallel Numerical Simulations of Aircraft Control Surface Buzz on Unstructured Grids

YAO Hong， GUO Chengpeng， ZHANG Ying

(AVIC Aerodynamics Research Institute， Aviation Key Laboratory of Science and Technology on Aerodynamics of High Speed and High Reynolds Number， Shenyang 110034， China)

In order to perform the numerical simulation of the buzz of aircraft control surface, a MPI parallel computation based on 3D unsteady Euler equations coupling with structural motion equations was established for cluster computers. In the solution process of pneumatic flow field, the spatial discretization was carried out using the centered finite volume method based on the unstructured grids. In addition, the time stepping was solved with the dual time method, and the structural motion equations were solved with the Adams predictor correction method. Aiming at the grid motion problem existing in the gap between the airfoil and control surface, the motion and deformation of grids were realized with Delaunay image mapping method in the unstructured grid system. Finally, the computation method was verified with the standard buzz computation model for the aircraft control surface. The results validated the established parallel computation method and indicated the excellent computation efficiency of the proposed model.

buzz； aircraft； control surface； parallel computation； dynamic mesh； Adams predictor-correction method

10.11830/ISSN.1000-5013.201606004

2016-10-10

幺虹(1963-)，女，高級工程師，主要從事飛行器非定?？諝鈩恿W的研究.E-mail:617162950@qq.com.

航空科學基金資助項目(2014ZA27003)

V 211.3

A

1000-5013(2016)06-0676-05