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

幺虹， 郭承鵬， 張穎

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

?

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

幺虹， 郭承鵬， 張穎

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

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

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

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

1 計(jì)算方法

1.1 非定常氣動(dòng)力計(jì)算方法

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

(1)

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

1.2 操縱面結(jié)構(gòu)運(yùn)動(dòng)方程求解

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

(2)

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

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

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

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

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

步驟3 背景網(wǎng)格的運(yùn)動(dòng)和變形.

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

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

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

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

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

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

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

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

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

步驟2 氣動(dòng)力矩計(jì)算.每個(gè)物理時(shí)間步開始之前并行化計(jì)算氣動(dòng)力矩.

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

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

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

步驟6 進(jìn)行下一時(shí)間步的流場(chǎng)計(jì)算.

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

3 NASP標(biāo)模嗡鳴計(jì)算

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

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

圖3 NASP嗡鳴計(jì)算標(biāo)準(zhǔn)模型 圖4 NASP模型表面及對(duì)稱面計(jì)算網(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 不同風(fēng)洞試驗(yàn)條件下操縱面偏角隨時(shí)間的變化歷程Fig.5 Change history deflection angle of control surface with time under experimental conditions of different wind tunnels

3.2 計(jì)算結(jié)果分析

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

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

4 結(jié)束語

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

[1] 葉正寅，張偉偉，史愛明.流固耦合力學(xué)基礎(chǔ)及其應(yīng)用[M].哈爾濱:哈爾濱工業(yè)大學(xué)出版社，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] 劉千剛，代捷，白俊強(qiáng).跨音速操縱面嗡鳴Hopf 分叉分析及結(jié)構(gòu)參數(shù)對(duì)嗡鳴特性影響的研究[J].航空學(xué)報(bào)，1999，20(6):527-532.

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

[7] 張偉偉，葉正寅，史愛明，等.基于Euler 方程的B 型和C 型嗡鳴特性數(shù)值研究[J].振動(dòng)工程學(xué)報(bào)，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ìn)展[J].力學(xué)進(jìn)展，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.

(責(zé)任編輯： 黃仲一 英文審校： 崔長(zhǎng)彩)

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-)，女，高級(jí)工程師，主要從事飛行器非定?？諝鈩?dòng)力學(xué)的研究.E-mail:617162950@qq.com.

航空科學(xué)基金資助項(xiàng)目(2014ZA27003)

V 211.3

A

1000-5013(2016)06-0676-05