飛行器操縱面嗡鳴的非結構網格并行計算方法

幺虹， 郭承鵬， 張穎

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

?

飛行器操縱面嗡鳴的非結構網格并行計算方法

幺虹， 郭承鵬， 張穎

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

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

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

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

1 計算方法

1.1 非定常氣動力計算方法

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

(1)

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

1.2 操縱面結構運動方程求解

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

(2)

1.3 非結構動網格方法

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

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

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

步驟2 建立計算網格節點和背景網格間的對應關系，計算插值權重系數.

步驟3 背景網格的運動和變形.

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

(a) 運動前 (b) 運動后圖1 非結構網格運動效果圖Fig.1 Schematic motion effect of unstructured grids

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

2 嗡鳴程序并行化設計

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

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

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

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

步驟1 前處理.使用Metis進行網格分區，Delaunay背景網格生成及網格點Delaunay背景圖映射關系建立及權值的計算.

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

步驟3 結構運動方程求解.僅在主進程進行結構運動方程求解，然后將求解得到的操縱面偏轉等信息發送到所有進程，并進行時間層的切換.

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

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

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

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

3 NASP標模嗡鳴計算

3.1 計算模型及網格

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

圖3 NASP嗡鳴計算標準模型 圖4 NASP模型表面及對稱面計算網格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 計算結果分析

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

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

4 結束語

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

[1] 葉正寅，張偉偉，史愛明.流固耦合力學基礎及其應用[M].哈爾濱:哈爾濱工業大學出版社，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 分叉分析及結構參數對嗡鳴特性影響的研究[J].航空學報，1999，20(6):527-532.

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

[7] 張偉偉，葉正寅，史愛明，等.基于Euler 方程的B 型和C 型嗡鳴特性數值研究[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] 周璇，李水鄉，孫樹立，等.非結構網格變形方法研究進展[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，FACULTY 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