黃江濤，高正紅，白俊強(qiáng)，周 鑄，趙 軻

（1.西北工業(yè)大學(xué) 翼型葉柵空氣動力學(xué)國防科技重點實驗室，西安 710072；2.中國空氣動力研究與發(fā)展中心，綿陽 621000）

0 引 言

飛機(jī)氣動彈性特性研究傳統(tǒng)上劃為靜氣動彈性和顫振特性研究兩個獨立部分［1－2］。靜氣動彈性特性計算求解過程往往采用靜變形和定常氣動力的交替迭代計算，最后得出發(fā)散速度或特定飛行條件下的機(jī)翼靜變形。傳統(tǒng)的飛行器設(shè)計主要是針對型架外形設(shè)計。然而在飛行器經(jīng)過氣動彈性變形之后，巡航狀態(tài)氣動性能與設(shè)計性能存在一定的區(qū)別。而對于民用客機(jī)來講首先進(jìn)行巡航外形設(shè)計，通過靜氣動彈性計算、修正得到型架外形，從而保證型架外形在巡航狀態(tài)飛行時，在氣動彈性作用下氣動外形與性能能夠恢復(fù)到巡航設(shè)計狀態(tài)。在型架修正及分析計算方面國內(nèi)外采用不同的方法進(jìn)行了不同程度的研究。例如國內(nèi)楊國偉采用稀疏背景網(wǎng)格進(jìn)行動網(wǎng)格變形、RBF進(jìn)行流固耦合，從而實現(xiàn)型架修正以及靜氣彈分析研究［3］，詹浩等采用反設(shè)計以及簡化結(jié)構(gòu)模型方法進(jìn)行了機(jī)翼型架設(shè)計及分析計算［4］，國外C B Allen等采用RBF結(jié)合TFI技術(shù)進(jìn)行了氣動彈性的分析計算［5］，取得了合理的計算結(jié)果。然而對于氣動彈性計算來講，變形網(wǎng)格的工作效率一定程度上決定了氣動彈性分析計算效率，尤其在非定常動氣動彈性計算中需要不斷的調(diào)用變形網(wǎng)格技術(shù)，這時網(wǎng)格變形工作效率更為重要。

對于靜氣動彈性來說，其關(guān)鍵技術(shù)為空間網(wǎng)格變形技術(shù)以及CFD／CSD數(shù)據(jù)交換技術(shù)。其中變形網(wǎng)格技術(shù)效率決定著靜氣動彈性計算效率、精度，其穩(wěn)健性決定了能否進(jìn)行飛行器大運動量下網(wǎng)格變形且能否保持較高網(wǎng)格質(zhì)量，以及所求氣動彈性變形外形的光滑程度；CFD／CSD數(shù)據(jù)交換技術(shù)決定著空氣動力學(xué)模型與結(jié)構(gòu)模型的數(shù)據(jù)能否進(jìn)行高效合理的交換。

本文基于徑向基函數(shù)與圖映射技術(shù)建立了高效的網(wǎng)格變形技術(shù)，與相關(guān)文獻(xiàn)不同的是，本文在合理提取飛行器表面稀疏點序列前提下，有效的利用多塊網(wǎng)格的塊頂點進(jìn)行Delaunay圖映射單元的構(gòu)建（網(wǎng)格分塊盡量均勻），避免了傳統(tǒng)Delaunay圖單元構(gòu)建時，相關(guān)點輸入的復(fù)雜操作，并通過RBF技術(shù)對網(wǎng)格頂點進(jìn)行高效操作，實現(xiàn)映射單元變形；通過建立共用的“映射曲面”，實現(xiàn)了CFD／CSD數(shù)據(jù)高效高精度交換，構(gòu)建了穩(wěn)健性強(qiáng)的靜氣動彈性計算平臺；利用該技術(shù)對某型支線客機(jī)進(jìn)行了型架外形設(shè)計、修正，進(jìn)一步通過靜氣動彈性分析計算，驗證本文所建立方法的正確性。

1 Delaunay圖映射技術(shù)

本文建立了基于Delaunay圖映射的網(wǎng)格變形技術(shù)。Dirichlet圖是由一組連接兩相鄰點直線的垂直平分線組成的連續(xù)多邊形，而Delaunay三角網(wǎng)是由與相鄰Dirichlet多邊形共享一條邊的相關(guān)點連接而成的三角形，換句話說Dirichlet三角網(wǎng)是Delaunay圖的偶圖。

基于Delaunay圖映射的網(wǎng)格變形技術(shù)基本依據(jù)為：在平面或空間給定包括流場邊界點在內(nèi)的一組點，即可以進(jìn)行唯一的Delaunay圖三角化。在Delaunay圖四面體單元構(gòu)建時，依據(jù)球面準(zhǔn)則進(jìn)行交換測試形式的優(yōu)化，根據(jù)插入點坐標(biāo)定位點所在格網(wǎng)，如果該格網(wǎng)記錄有三角形編號，則從該編號開始在相鄰三角形中定位插入點所在三角形；由近及遠(yuǎn)遍歷格網(wǎng)的相鄰網(wǎng)格，若存在三角編號則遍歷結(jié)束；在最新生成三角形中定位插入點。從而可以完成計算區(qū)域的三角化圖覆蓋，進(jìn)一步對空間計算域任意點進(jìn)行四面體映射單元定位。

計算域網(wǎng)格節(jié)點定位完畢后，需建立映射關(guān)系［6］，對于平面Delaunay三角化背景網(wǎng)格，建立如圖1的計算網(wǎng)格與平面Delaunay三角單元關(guān)系，對于任意網(wǎng)格節(jié)點O處于三角形MNQ中，由O和三角形頂點構(gòu)造的三角形MON、QOM以及NOQ的面積分別為S1、S2、S3，令：

則可以建立Delaunay三角單元節(jié)點與計算網(wǎng)格節(jié)點如下關(guān)系式：

飛行器氣動外形變形時，相當(dāng)于Delaunay圖背景網(wǎng)格三角單元變化，即 ［XMXNXQ］更新為［X′MX′NX′Q］，對于三維網(wǎng)格變形方法與二維原理一致，如圖2所示。依據(jù)建立的Delaunay圖映射關(guān)系式，計算網(wǎng)格節(jié)點更新為：

以氣動彈性變形后飛行器表面區(qū)域頂點位移量為輸入量，建立RBF徑向基函數(shù)插值系統(tǒng)，進(jìn)行空間區(qū)域頂點位移量插值，操作空間網(wǎng)格區(qū)域頂點運動，從而進(jìn)行Delaunay圖四面體映射單元的變形，利用式（2）進(jìn)行計算網(wǎng)格的變形，完成網(wǎng)格更新。而Delaunay映射變換之前首先要完成空間區(qū)域頂點的更新，即采用RBF徑向基函數(shù)技術(shù)，其工作示意圖見圖3。

圖1 網(wǎng)格節(jié)點與Delaunay三角形Fig.1 Grid points and Delaunay triangle

圖2 網(wǎng)格節(jié)點與Delaunay四面體Fig.2 Grid points and Delaunay tetrahedron

圖3 RBF工作示意圖Fig.3 RBF work map

2 RBF徑向基函數(shù)

RBF徑向基函數(shù)基本概念由Buhmann和Wendland提出［7－8］，該技術(shù)在精確插值、神經(jīng)網(wǎng)絡(luò)構(gòu)建、數(shù)據(jù)預(yù)測等方面獲取了廣泛的應(yīng)用［9－13］。基于RBF方法的插值數(shù)學(xué)模型可以表示為：

本文對p（x）的選取采用多項式方法：

其中M 為徑向基中心矢量的維數(shù)。p（x）的建立不僅能夠精確的描述物體或網(wǎng)格的平移以及旋轉(zhuǎn)運動方式，且進(jìn)行徑向基插值時能夠保證輸入點系統(tǒng)“力平衡”及“力矩平衡”［14］。其中，xi為徑向基中心，相當(dāng)于圖3中ci；αi為插值系數(shù)，相當(dāng)于圖3中λi。‖x－xi‖為x到中心xi的徑向歐氏距離。且徑向基函數(shù)φ具有以下性質(zhì)：

以輸入數(shù)據(jù)為RBF技術(shù)構(gòu)建基礎(chǔ)，可以寫為如下矩陣形式：

其中：Φ ＝［φij］，φij＝φ（‖xi－xj‖），（i，j＝1，2，…，N ）；Φ 為插值矩陣，ΓT＝Φ－1Y，Γ ＝ ［λ1，λ2，…，λN］。從而實現(xiàn)RBF徑向基函數(shù)框架構(gòu)建，RBF（x）則為通過所有輸入點的連續(xù)可微函數(shù)。

3 CFD／CSD數(shù)據(jù)耦合技術(shù)與靜力學(xué)方程求解

在CFD／CSD數(shù)據(jù)交換時，本文采用構(gòu)建空氣動力模型與結(jié)構(gòu)模型的共用“映射曲面”［15］，根據(jù)機(jī)翼形狀建立一曲面，該曲面由兩個方向邏輯參數(shù)進(jìn)行描述，網(wǎng)格中第i行，第j列處的節(jié)點坐標(biāo)為（x，y），且該點處的參數(shù)化坐標(biāo)定義為（u ，v），u ＝i，v ＝j(luò)。將映射網(wǎng)格控制點參數(shù)化后得到參數(shù)空間。

對于CFD／CSD網(wǎng)格交界面上節(jié)點的參數(shù)化就是要計算該點在參數(shù)空間中的坐標(biāo)，首先利用面積坐標(biāo)方法［15］定位源／目標(biāo)學(xué)科節(jié)點所處映射曲面單元（i，j），進(jìn)一步通過雙線性曲面插值計算該點在其所在的結(jié)構(gòu)化映射網(wǎng)格單元內(nèi)的位置（i0，j0），則該節(jié)點在參數(shù)空間的位置為（i＋i0，j＋j0）。

本文采用雙三次Bezier曲面［15－16］進(jìn)行源學(xué)科節(jié)點與目標(biāo)學(xué)科節(jié)點變量插值：

其中，U、V為目標(biāo)學(xué)科參數(shù)空間邏輯坐標(biāo)，P為源學(xué)科節(jié)點矩陣，M 為Bezier基函數(shù)系數(shù)矩陣。

圖4 “映射曲面”工作示意圖Fig.4 Mapping surface work map

采用柔度系數(shù)矩陣方法［17－18］求解結(jié)構(gòu)靜力學(xué)方程時，僅考慮飛機(jī)機(jī)翼的彈性變形，發(fā)動機(jī)及其掛架運動與機(jī)翼連接處保持一致。

其中，i＝（1，2，…，n），常數(shù)cij為柔度影響系數(shù)矩陣。則結(jié)構(gòu)靜力學(xué)方程為：

qs為結(jié)構(gòu)節(jié)點在三個方向的變形位移矩陣；C為結(jié)構(gòu)柔度矩陣；Fs為作用在節(jié)點上的氣動力矩陣。對于型架構(gòu)型逆計算來講，只需將式（7）中力方向反向。綜上所述，本文建立的靜氣動彈性數(shù)值模擬流程見圖5。

圖5 靜氣動彈性計算流程Fig.5 Static aeroelastic work flow chart

4 某型支線客機(jī)型架外形修正與靜氣彈計算

以某型支線客機(jī)為型架外形為計算算例，進(jìn)行型架設(shè)計后修型，進(jìn)一步進(jìn)行靜氣動彈性計算，驗證所求型架構(gòu)型正確性以及本文相關(guān)技術(shù)的效率、可靠性。

該型客機(jī)設(shè)計巡航狀態(tài)為：

該型客機(jī)結(jié)構(gòu)有限元采用桿梁單元，建立雙梁單塊式結(jié)構(gòu)模型，主要結(jié)構(gòu)元素包括前梁、后梁、翼肋等。利用有限元分析計算了該結(jié)構(gòu)模型的剛度矩陣?？諝鈩恿W(xué)模型采用多區(qū)域點對接結(jié)構(gòu)網(wǎng)格進(jìn)行空間離散，空間流場流場共分為330個區(qū)域，表面劃分為160個區(qū)域，共1300萬個網(wǎng)格單元（半模）。圖6給出了該支線客機(jī)的表面網(wǎng)格分布示意圖。氣動分析主控方程為Navier－Stokes方程；空間離散方法為二階Roe格式；湍流模型采用一方程SA模型；采用LU－SGS隱式時間推進(jìn)。通過線性插值的方法求出設(shè)計升力系數(shù)下所對應(yīng)的巡航攻角α＝1.95°。圖7為巡航狀態(tài)下飛機(jī)表面壓力云圖。

根據(jù)CFD數(shù)值模擬結(jié)果，采用文中建立的CFD／CSD共用的“映射曲面”技術(shù)進(jìn)行氣動載荷向結(jié)構(gòu)載荷的映射，通過求解靜力學(xué)方程，求出結(jié)構(gòu)節(jié)點位移，再一次通過“映射曲面”技術(shù)將結(jié)構(gòu)位移映射到氣動網(wǎng)格上，從而完成型架構(gòu)型的逆向計算，如圖所示，綠色外形為型架構(gòu)型，淺色為巡航設(shè)計外形，可以看出機(jī)翼展向下彎較為明顯，表1給出了展向10個站位型架外形機(jī)翼與巡航外形的相對扭角。

對于生產(chǎn)與裝配工藝來講，要求機(jī)翼前緣為一直線且與巡航設(shè)計外形重合，因此在保持型架構(gòu)型不變的情況下對機(jī)翼做前緣一致性修正，以滿足裝配工藝要求。

圖6 CFD表面網(wǎng)格Fig.6 CFD surface grid

圖7 巡航狀態(tài)表面壓力云圖Fig.7 Pressure contour of cruise shape

表1 不同展向位置的機(jī)翼相對扭角Table 1 Relative torsional angles of the wing at different spanwise sections

利用文中建立的相關(guān)技術(shù)對型架外形進(jìn)行靜氣彈計算，本文主要驗證了型架構(gòu)型在巡航狀態(tài)下氣動特性。在利用稀疏表面網(wǎng)格節(jié)點與空間網(wǎng)格頂點構(gòu)建Delaunay圖時，Delaunay單元數(shù)為6090個，計算時間為30s，之后由稀疏表面網(wǎng)格節(jié)點與空間網(wǎng)格頂點進(jìn)行映射單元變形。本文采用以上方法既避免了Delaunay單元構(gòu)建輸入點的復(fù)雜操作，又從很大程度上降低了RBF大型矩陣的求解量，從而大幅度提高網(wǎng)格變形效率。表2給出了相同網(wǎng)格條件下（以DLR－F6翼身組合體為對象），不同變形網(wǎng)格技術(shù)與本文方法工作效率的對比，可以看出本文提出的新型變形網(wǎng)格方法具備較高的網(wǎng)格更新速度。圖10給出了經(jīng)過靜氣動彈性計算后不同網(wǎng)格區(qū)域頂點位置在RBF徑向基函數(shù)技術(shù)操作下的更新。圖11進(jìn)一步給出了空間頂點更新后對應(yīng)Delaunay圖空間截面。在該計算模型靜氣彈計算過程中，RBF徑向基函數(shù)技術(shù)對330個不同網(wǎng)格區(qū)域1000個頂點的操作僅需要0.05s時間，而對于頂點更新后Delaunay圖單元變形進(jìn)行的空間1300萬網(wǎng)格節(jié)點變形僅需要1.5s，可見本文構(gòu)建的徑向基函數(shù)與Delaunay圖映射網(wǎng)格變形技術(shù)效率較高，比較適用于大型工程靜氣動彈性計算中。整個靜氣動彈性經(jīng)過10次迭代基本收斂。圖12給出了型架構(gòu)型巡航外形與巡航設(shè)計外形的對比，可以看出兩者基本完全重合，驗證了修正方法的正確性與靜氣動彈性分析的可靠性與較高效率。表3給出了兩者在氣動特性方面的對比，兩者氣動特性接近，誤差較小。圖13～圖16給出了展向典型站位壓力分布對比，壓力分布形態(tài)相同基本重合。

表2 不同變形網(wǎng)格技術(shù)效率對比Table 2 Efficient comparison between different grid deformation methods

圖8 型架外形與巡航外形對比Fig.8 Jig shape and cruise shape

圖9 型架外形與巡航外形局部視圖Fig.9 Local view of jig shape and cruise shape

圖10 RBF操作下區(qū)域頂點分布Fig.10 Block vertexs distribution under RBF operation

圖11 Delaunay四面體空間截面Fig.11 Delaunay space section

圖12 巡航設(shè)計外形與型架巡航外形對比Fig.12 Comparision of cruise designed shape and jig cruise shape

表3 巡航設(shè)計外形與型架巡航外形氣動特性對比Table 3 Comparison of aerodynamic performance between cruise design shape and jig shape at cruise condition

圖13 y／b＝0.19展向位置壓力分布Fig.13 Pressure distribution of y／b＝0.19spanwise section

圖14 y／b＝0.3展向位置壓力分布Fig.14 Pressure distribution of y／b＝0.3spanwise section

圖15 y／b＝0.5展向位置壓力分布Fig.15 Pressure distribution of y／b＝0.5spanwise section

圖16 y／b＝0.8展向位置壓力分布Fig.16 Pressure distribution of y／b＝0.8spanwise section

5 結(jié)束語

本文研究了采用RBF徑向基函數(shù)、Delaunay圖映射技術(shù)的網(wǎng)格變形方法。研究建立了基于靜氣動彈性計算的飛行器型架外形修正設(shè)計、驗證方法。對某型支線客機(jī)進(jìn)行了型架外形逆推、修正及靜氣動彈性驗證，得出如下結(jié)論：

（1）基于RBF徑向基函數(shù)、Delaunay圖映射的網(wǎng)格變形技術(shù)充分利用多塊網(wǎng)格頂點，在避免Delaunay圖單元構(gòu)建時輸入點的復(fù)雜操作的同時充分利用RBF對網(wǎng)格頂點的高效操作，與傳統(tǒng)變形網(wǎng)格技術(shù)相比，計算效率較高，適應(yīng)大幅度網(wǎng)格變形，適合于大型工程靜氣動彈性計算。

（2）逆推的飛行器型架構(gòu)型表面光滑，變形合理，說明所采用“映射曲面”CFD／CSD數(shù)據(jù)交換技術(shù)插值精度高。

（3）通過對型架外形的氣動彈性計算，驗證了文中型架構(gòu)型修正的正確性以及靜氣動彈性分析方法的可靠性、高效性。

［1］XIE C C，WU Z G，YANG C.Aeroelastic analysis of flexible large aspect ratio wing［J］.Journal of Beijing University of Aeronautics and Astronautics，2003，29（12）：1087－1090.（in Chinese）謝長川，吳志剛，楊超.大展弦比柔性機(jī)翼的氣動彈性分析［J］.北京航空航天大學(xué)學(xué)報，2003，29（12）：1087－1090.

［2］GUAN D.Aeroelastic experiment［M］.Beijing Aeronautics College press，1986.（in Chinese）管德.氣動彈性試驗［M］.北京航空學(xué)院出版社，1986.

［3］YANG G W，ZHENG G N.Aircraft jig shape correction method based on static aeroelastic analyses［J］.Advances in Aeronautical Science and Engineering，2011，2（2）：113－150.（in Chinese）楊國偉，等.基于靜氣動彈性效應(yīng)的飛機(jī)型架外形修正方法研究［J］.航空工程進(jìn)展，2011，2（2）：113－150.

［4］ZHAN H，CHENG S X，ZHU J.An effective aerodynamic and aeroelastic coupled design method for complicated wing shape［J］.Journal of Northwestern Polytechnical University，2009，27（1）：100－104.（in Chinese）詹浩，等.考慮氣動彈性影響的機(jī)翼復(fù)雜氣動外形設(shè)計研究［J］.西北工業(yè)大學(xué)學(xué)報，2009，27（1）：100－104.

［5］ALLEN C B.Unified approach to CFD－CSD interpolation and mesh motion using radial basis functions［R］.AIAA Paper 2007－3084，2007.

［6］LIU X Q，QIN N.Fast dynamic grid deformation based on Delaunay graph mapping［J］.Journal of Computational Physics，2006，211：405－423.

［7］BUHMANN M.Radial basis functions［M］.Cambridge University Press，2005.

［8］CIZMAS P，GARGOLOFF J.Mesh generation and deformation algorithm for aeroelasticity simulations［R］.AIAA 2007－556.45th Aerospace Sciences Meeting［C］.Reno，NV，2007.

［9］PARK J，SANDBERG I W.Universal approximation using radialbasis－function networks［J］.Neural Computation，1991.

［10］LIGHT W A.Some aspects of radial basis function approximation.approximation theory，spline functions and applications［J］.NATO ASI Series，1992，356：163－190.

［11］KARIM A，ADELI H.Radial basis function neural network for work zone capacity and queue estimation［J］.J.Transp.Eng.，2003，129（5）：494－503.

［12］ZHANG X，SONG K Z，LU M W，et al.Meshless methods based on collocation with radial basis functions［J］.Computational Mechanics，26（4）：333－343.

［13］SARLER B.A radial basis function collocation approach in computational fluid dynamics［J］.CMES，2005，7（2）：185－193.

［14］WENDLAND H.Fast evaluation of radial basis functions：Methods based on partition of Unity［M］／／Approximation Theory X：Wavelets，Splines，and Applications.Nashville，Texas，USA：Vanderbilt University Press，2002：473－483.

［15］LI L Z.Turbine blade temperature transfer using the load surface method［J］.Computer Aided Design，2007，39：494－505.

［16］ZHU X X.Free curves and surfaces modeling technology［M］.Science Press，2000.（in Chinese）朱心雄.自由曲線曲面造型技術(shù)［M］.科學(xué)出版社，2000.

［17］ZHAO Y H.Aeroelastic dynamics and control［M］.Science Press，2007.（in Chinese）趙永輝.氣動彈性力學(xué)與控制［M］.科學(xué)出版社，2007.

［18］SHEN K Y，GUAN D.The principle of elastic mechanics［M］.Shanghai science and Technology Press，1982.（in Chinese）沈克揚，管德.彈性力學(xué)原理［M］.上?？茖W(xué)技術(shù)文獻(xiàn)出版社，1982.