林言中 陳 兵 徐 旭

(北京航空航天大學 宇航學院，北京100191)

氣動彈性力學[1]涉及到非定常流場、結構應力/應變場的相互作用，屬于典型的流固耦合[2]問題.隨著航空航天技術的發(fā)展，未來飛行器普遍具有輕質、大翼展的特點，結構在氣動載荷作用下容易產生變形，而結構的變形又會反作用于流場，改變流體載荷的分布和大小，氣動載荷與彈性結構響應之間的耦合現(xiàn)象十分突出，甚至可能發(fā)生操縱面反轉、顫振等不穩(wěn)定現(xiàn)象，使得氣動彈性問題成為飛行器設計過程中不得不考慮的關鍵問題之一[3].

氣動彈性問題的計算涉及流體域與固體域各自的求解，耦合作用僅發(fā)生在兩相交界面上，通過引入平衡和協(xié)調關系實現(xiàn)兩場的統(tǒng)一.氣動力預測與結構動態(tài)響應分析方法是氣動彈性問題的核心內容.早期對氣動彈性問題的研究多在頻域內進行，采用線性結構模型和線化氣動力模型進行分析.這種方法在高速、大攻角、大變形的條件下精度不高，容易失效[3].隨著計算流體力學(CFD，Computational Fluid Dynamics)、計算結構力學(CSD，Computational Structure Dynamics)以及計算機技術的發(fā)展，高精度的CFD/CSD耦合求解[4]成為氣動彈性問題的有效研究方法:流場與結構場在時域內各自推進，在交界面上通過傳遞壓力、位移等參數(shù)實現(xiàn)耦合計算.CFD求解完全的三維非定常流場，由流場結果直接計算任意時刻作用在彈性體表面的氣動力分布，精度高、適用范圍廣;CSD求解也可考慮結構的幾何非線性及材料非線性，使計算結果更加真實.

在CFD/CSD耦合求解過程中，由于結構變形破壞了流場的拓撲結構，使得流場的計算網格在每個物理時間步內需要更新;由于流場網格與結構網格在界面上往往不匹配，需要通過一定的插值方法進行參數(shù)的傳遞.可見，除了CFD與CSD求解器外，動網格與界面數(shù)據(jù)傳遞方法也是氣動彈性計算的關鍵技術.動網格算法要求高的計算效率以及強的適應大變形能力;界面數(shù)據(jù)傳遞方法要求精確插值以及滿足運動、動力連續(xù)和能量守恒等條件[5]，任何一個環(huán)節(jié)出現(xiàn)問題都會影響耦合計算的精度和效率.

對于動網格問題，Batina[6]首先提出線彈性模型，之后由 Anderson[7]及 Farhat等[8]進行了改進和補充，極大地提高了其計算能力，使得彈簧比擬方法成為當前最常用的動網格方法之一.然而繁瑣的數(shù)據(jù)結構以及計算效率、適應大變形能力等因素制約了彈簧比擬法在氣動彈性問題中的應用.對于界面數(shù)據(jù)傳遞問題，無限平板樣條法、多重二次曲面-雙調和法、常體積轉換法等方法被廣泛地應用于氣動彈性計算[5]，取得了一系列成果，但尋求高精度、高效率的界面數(shù)據(jù)傳遞方法仍是氣動彈性計算發(fā)展的重點問題.

近年來，以徑向基函數(shù)(RBF，Radial Basis Function)為基礎的插值方法引起了數(shù)學界的廣泛關注，F(xiàn)rank[9]在大量實例中對各種散亂插值方法進行了對比，發(fā)現(xiàn)徑向基函數(shù)插值的效果最佳.徑向基函數(shù)插值方法已經被廣泛地應用于工程實際問題之中，Beckert[10]最早利用 RBF 插值處理機翼變形問題，Allen等[11]進一步將其用于網格運動，取得了理想的效果.

國內，蘇波等[12]分析了徑向基函數(shù)應用于流固耦合界面插值的可行性，并進行了任一時間步上的流固交互作用分析;林言中等[13]將RBF插值用于網格運動，并將其與彈簧比擬方法進行對比，說明了RBF方法在計算效率和適應大變形能力上的優(yōu)勢.

本文將基于徑向基函數(shù)插值的動網格及界面數(shù)據(jù)傳遞方法用于CFD/CSD耦合計算，通過對AGARD 445.6機翼顫振問題的非定常數(shù)值模擬驗證了該方法的可行性.徑向基函數(shù)插值用于界面數(shù)據(jù)傳遞能夠充分地保證運動連續(xù)、動力連續(xù)及能量守恒條件且兼有整體插值與局部插值方法的特點;用于動網格技術有著數(shù)據(jù)結構簡單、適應大變形能力強等優(yōu)勢，為氣動彈性問題的計算提供了非常有力的技術支撐.

1 數(shù)值方法

1.1 流動控制方程及數(shù)值方法

動網格上，積分形式的可壓縮歐拉方程[14]可寫成以下形式:

式中，W為守恒變量;F為對流通量，即

式中，V為流體質點的速度矢量;Vw為網格運動速度矢量;Ix，Iy，Iz分別為慣性系下沿坐標軸方向的單位矢量.

對式(1)采用有限體積法以一定的差分格式進行空間離散后，可得如下半離散方程:

本文的空間離散方案采用二階精度的Van Leer格式.對式(3)可采用雙時間步法求解.

1.2 結構動力學方程及數(shù)值方法

瞬態(tài)動力學平衡方程[15]為

式中，[M]為質量矩陣;[C]為阻尼矩陣;[K]為剛度矩陣;{F(t)}為隨時間變化的外載荷;分別為節(jié)點加速度、速度及位移向量.采用有限元分析中的Newmark方法對式(4)進行求解，可得到結構各時刻的位移和應力分布.阻尼矩陣采用瑞利阻尼:

1.3 CFD/CSD耦合計算方法

當前針對CFD/CSD耦合已經發(fā)展了全耦合、松耦合及緊耦合3種方式，分別對應不同的時間推進方案.全耦合方案具有很高的時間精度，但非線性程度高，難以應用于復雜工程問題;松耦合方法計算復雜度低，有利于實現(xiàn)程序模塊化，但僅有一階時間精度;緊耦合方法可取較大的時間步長，但為滿足時間精度需反復內迭代，計算耗時.本文采用二階時間精度的松耦合格式，具體計算過程參見文獻[16].

2 徑向基函數(shù)插值

2.1 徑向基函數(shù)及插值方法

徑向基函數(shù)以空間距離 x為基本變量，通常可分為全局(global)函數(shù)、局部(local)函數(shù)和緊支(compact)函數(shù)[17]3 類，如表 1 所示.緊支函數(shù)在其中心距離達到某一特定值sr(稱為緊支半徑support radius)后函數(shù)值恒為0.表1中緊支函數(shù)的變量ζ為中心距離與緊支半徑之比，即ζ=，當 ζ＞1 時，函數(shù)值為0.

表1 常用的RBF基函數(shù)

令d維歐氏空間中的一組中心點X={xc1，xc2，…，xcn}對應的標量值為 g={gc1，gc2，…，gcn}，則對于徑向基函數(shù)存在連續(xù)函數(shù):

通過所有中心點，確定插值系數(shù)γ的過程即為徑向基函數(shù)插值.其中，為歐氏距離.對于三維空間而言，可直接用r表示為

2.2 徑向基函數(shù)插值用于動網格

針對氣動彈性計算中的動網格問題，可直接應用式(6)確定RBF插值函數(shù).此時中心點即為耦合界面上流場區(qū)域的網格點，而相應的標量值為流體網格點3個方向上的位移分量Dx，Dy，Dz.插值系數(shù)可由以下線性系統(tǒng)確定:

式中n階方陣M為插值矩陣，其元素為

矩陣M對稱正定，可用共軛梯度法[18]快速求解.將插值系數(shù) γx，γy，γz代入式(6)可得 3 個坐標分量的插值函數(shù).將流場計算域內每個網格點的坐標代入RBF插值函數(shù)即可得到內場網格點的位移，從而實現(xiàn)整個計算域內的網格移動:

在氣動彈性計算過程中，每個時間步的網格更新都可從零時刻的網格出發(fā)，這樣位移插值系數(shù)只需計算一次，極大地提高了計算效率.需要指出的是，對于流體區(qū)域的網格運動問題，耦合界面的網格運動僅對界面附近區(qū)域的網格有較大影響，對于那些遠離界面的網格點可認為位移很小甚至不運動.緊支基函數(shù)可通過調整緊支半徑來控制中心點的影響區(qū)域，此特點與動網格的要求相一致，且求解方程組的系數(shù)矩陣呈稀疏、帶狀分布，對于存儲量和計算效率都有著極強的適應能力，因此在動網格計算中更傾向于使用緊支型基函數(shù).

2.3 徑向基函數(shù)插值用于界面數(shù)據(jù)傳遞

針對氣動彈性計算中的界面數(shù)據(jù)傳遞問題，中心點變?yōu)轳詈辖缑嫔辖Y構區(qū)域的網格點，而標量值則為結構網格點3個方向上的位移分量.為了使RBF插值函數(shù)滿足平移及旋轉不變性，通常需在式(6)的基礎上加上d維低階多項式p(x)，有

對于小于等于2階的正定或條件正定基函數(shù)，p(x)可以采用如下線性多項式:

式(11)的2個定解條件為

式中q(x)為滿足deg(q(x))≤deg(p(x))的所有多項式，由式(12)可知，q為階數(shù)不超過1的多項式，即 q(x)的基函數(shù)為{1，x，y，z}.

插值系數(shù) γi及多項式系數(shù) β0，β1，β2，β3可由以下線性系統(tǒng)確定:

式中矩陣M元素同式(9)，矩陣P為

由于式(15)有零塊，不宜直接求解線性方程組，可通過矩陣運算求得

由式(11)、式(17)及式(18)可建立流體網格與結構網格在界面上的位移傳遞關系:

得到位移傳遞矩陣H后，可根據(jù)能量守恒原理傳遞表面壓力:

在氣動彈性計算中，若整個過程結構變形不是很大，則可認為傳遞矩陣H基本不變，在初始時刻計算一次即可.對于基函數(shù)的選擇，若采用全局函數(shù)，則RBF插值是一種整體插值方法;若采用緊支函數(shù)，則當緊支半徑較小時，RBF插值呈現(xiàn)局部插值的特點.

3 AGARD 445.6機翼顫振算例

AGARD 445.6機翼顫振風洞實驗是由美國國家航天局蘭利研究中心在跨音速動態(tài)風洞中完成的，有著較為完備的風洞實驗數(shù)據(jù)[19].目前，AGARD 445.6機翼軟模型已經成為國際上跨音速氣動彈性程序考核的標準算例，通過對其顫振問題的計算來考察RBF方法在氣動彈性計算中的應用.AGARD 445.6機翼幾何參數(shù)為:展弦比1.644，梢根比 0.659 2，根弦長 0.558 7 m，1/4 弦線后掠角為45°，沿流向翼型為NACA 65A004.流場計算采用結構化網格，網格總數(shù)746 916.機翼幾何尺寸及對稱面上的網格分布如圖1所示.

圖1 AGARD 445.6機翼幾何參數(shù)及對稱面網格

進行結構有限元計算時，機翼根部為固支邊界條件.首先對機翼進行模態(tài)分析，各階固有頻率與實驗數(shù)據(jù)對比如表2所示.計算結果與實驗數(shù)據(jù)吻合較好，本文建立的機翼模型能夠正確模擬實驗模型的材料特性.

表2 AGARD 445.6機翼模態(tài)固有頻率 Hz

首先取來流馬赫數(shù)為0.499，靜溫297.2 K，攻角為0°的工況進行計算，計算初始時刻在翼尖處給一個微小法向初速度作為初始擾動.對于機翼顫振的氣動彈性問題，若來流速度小于臨界顫振速度，則振動收斂;若來流速度恰好等于臨界顫振速度，則機翼呈現(xiàn)等幅諧振的狀態(tài);若來流速度大于臨界顫振速度，則振動發(fā)散.通過調整無量綱來流速度Vf觀察機翼的振動情況，觀測對象為翼梢后緣點的法向位移.

無量綱速度定義為

式中，V∞為實際來流速度;bs為參考長度，一般定義為機翼的半根弦長;ωα為參考頻率，該問題定義為機翼的一階扭轉頻率;μ為質量比，定義為

式中，m為機翼質量;ρ∞為來流密度;V為體積系數(shù).

由圖2～圖4的計算結果可以看出，隨著來流無量綱速度的增加，機翼的動態(tài)響應曲線由收斂到發(fā)散，在Vf=0.455時響應曲線為等幅振動，此時機翼發(fā)生顫振，此速度為該馬赫數(shù)下機翼的顫振速度，顫振頻率為18.2 Hz.文獻[19]中實驗得到的機翼顫振速度為0.45，顫振頻率20.4 Hz，與本文計算結果較為吻合.

圖2 Vf=0.432，Ma=0.499時機翼的動態(tài)響應曲線

圖3 Vf=0.455，Ma=0.499時機翼的動態(tài)響應曲線

由實驗報告[19]可知，機翼跨聲速顫振計算的馬赫數(shù)范圍為0.499～1.141，本文取馬赫數(shù)為0.499，0.678，0.901，0.960 和 1.072 這 5 個狀態(tài)進行顫振點的預測，計算結果與實驗數(shù)據(jù)對比如圖5所示.

圖4 Vf=0.477，Ma=0.499時機翼的動態(tài)響應曲線

圖5 機翼臨界顫振速度邊界

由圖5可見，CFD/CSD耦合計算捕捉到了跨聲速“凹坑”現(xiàn)象，計算的臨界顫振速度在各來流馬赫數(shù)狀態(tài)下均與實驗數(shù)據(jù)符合較好，對AGARD 445.6機翼的顫振計算結果表明徑向基函數(shù)方法可以有效地應用于氣動彈性計算之中.

4 結論

本文將徑向基函數(shù)插值應用于氣動彈性計算的兩大關鍵技術:動網格及界面數(shù)據(jù)傳遞，介紹了徑向基函數(shù)的基本概念及其用于氣動彈性計算的數(shù)值方法.采用二階時間精度的松耦合方案，對AGARD 445.6機翼顫振問題進行了CFD/CSD耦合計算.在較寬泛的馬赫數(shù)范圍內捕捉到了臨界顫振速度及跨聲速凹坑現(xiàn)象，計算結果與實驗數(shù)據(jù)符合良好.結果表明，徑向基函數(shù)插值方法可以有效地應用于氣動彈性計算，為氣動彈性問題研究的發(fā)展提供了有力的技術途徑.

References)

[1]陳桂彬，鄒叢青，楊超.氣動彈性設計基礎[M].北京:北京航空航天大學出版社，2004

Chen Guibin，Zou Congqing，Yang Chao.Aeroelastic design[M].Beijing:Beihang University Press，2004(in Chinese)

[2]邢景棠，周盛，崔爾杰.流固耦合力學概述[J].力學進展，1997，27(1):19 －38

Xing Jingtang，Zhou Sheng，Cui Erjie.A survey on the fluid-solid interaction mechanics[J].Advances in Mechanics，1997，27(1):19－38(in Chinese)

[3]安效民，徐敏，陳士櫓.多場耦合求解非線性氣動彈性的研究綜述[J].力學進展，2009，39(3):284 －298

An Xiaomin，Xu Min，Chen Shilu.An overview of CFD/CSD coupled solution for nonlinear aeroelasticity[J].Advances in Mechanics，2009，39(3):284 － 298(in Chinese)

[4]Baum J D，Luo H，Mestreau E L，et al.Recent developments of a coupled CFD/CSD methodology[C]//15th AIAA Computational Fluid Dynamics Conference 2001.Reston，VA:AIAA，2001:1087－1097

[5]蘇波，錢若軍，袁行飛.流固耦合界面信息傳遞理論和方法研究進展[J].空間結構，2010，16(1):3 －10

Su Bo，Qian Ruojun，Yuan Xingfei.Advances in research on theory and method of data exchange on coupling interface for FSI analysis[J].Spatial Structures，2010，16(1):3 －10(in Chinese)

[6]Banita J T.Unsteady Euler airfoil solutions using unstructured dynamic meshes[J].AIAA Journal，1990，28(8):1381 － 1388

[7]Anderson W K，Venkatakrishnan V.Aerodynamic design optimization on unstructured grids with a continuous adjoint formulation[J].Computers ＆ Fluids，1999:28(4/5):443 －480

[8]Degand C，F(xiàn)arhat C.A three-dimensional spring analogy method for unstructured dynamic meshes[J].Computers and Structures，2002，80(3):305 －316

[9]Frank R.Scattered data interpolation:tests of some methods[J].Math Comp，1982，38(157):181 －200

[10]Beckert A，Wendland H.Multivariate interpolation for fluid-structure-interaction problems using radial basis functions[J].Aerospace Science and Technology，2001，5(2):125 －134

[11]Rendall T C S，Allen C B.Efficient mesh motion using radial basis functions with data reduction algorithms[J].Journal of Computational Physics，2009，228:6231 －6249

[12]蘇波，石啟印，錢若軍.徑向基函數(shù)應用于流固耦合分析初探[J].工程力學，2013，30(1):59 －63

Su Bo，Shi Qiyin，Qian Ruojun.Preliminary study on the use of radial basis function in fluid-structure interaction analysis[J].Engineering Mechanics，2013，30(1):59 －63(in Chinese)

[13]林言中，陳兵，徐旭.徑向基函數(shù)插值方法在動網格技術中的應用[J].計算物理，2012，29(2):191 －197

Lin Yanzhong，Chen Bing，Xu Xu.Radial basis function interpolation in moving mesh technique[J].Chinese Journal of Computational Physics，2012，29(2):191 －197(in Chinese)

[14]閻超.計算流體力學方法及應用[M].北京:北京航空航天大學出版社，2006

Yan Chao.Computational fluid dynamics method and its application[M].Beijing:Beihang University Press，2006(in Chinese)

[15]陳道禮，饒剛，魏國前.結構分析有限元法的基本原理及工程應用[M].北京:冶金工業(yè)出版社，2012

Chen Daoli，Rao Gang，Wei Guoqian.Finite element analysis and engineering application[M].Beijing:Metallurgical Industry Press，2012(in Chinese)

[16]安效民，徐敏，陳士櫓.二階時間精度的CFD/CSD耦合算法研究[J].空氣動力學學報，2009，27(5):547 －552

An Xiaomin，Xu Min，Chen Shilu.Analysis for second order time accurate CFD/CSD coupled algorithms[J].Acta Aerodynamica Sinica，2009，27(5):547 －552(in Chinese)

[17]Wendland H.On the smoothness of positive definite and radial functions[J].Computational and Applied Mathematics，1999，101(1):177－188

[18]Buhmann M D.Radial basis functions[J].Acta Numerica，2000，9:1 －38

[19]Yates E C.AGARD standard aeroelastic configurations for dynamic response I-wing 445.6[R].AGARD Report No.765，1988