基于径向基函数插值的气动弹性计算方法

林言中 陈 兵 徐 旭

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

气动弹性力学[1]涉及到非定常流场、结构应力/应变场的相互作用，属于典型的流固耦合[2]问题.随着航空航天技术的发展，未来飞行器普遍具有轻质、大翼展的特点，结构在气动载荷作用下容易产生变形，而结构的变形又会反作用于流场，改变流体载荷的分布和大小，气动载荷与弹性结构响应之间的耦合现象十分突出，甚至可能发生操纵面反转、颤振等不稳定现象，使得气动弹性问题成为飞行器设计过程中不得不考虑的关键问题之一[3].

气动弹性问题的计算涉及流体域与固体域各自的求解，耦合作用仅发生在两相交界面上，通过引入平衡和协调关系实现两场的统一.气动力预测与结构动态响应分析方法是气动弹性问题的核心内容.早期对气动弹性问题的研究多在频域内进行，采用线性结构模型和线化气动力模型进行分析.这种方法在高速、大攻角、大变形的条件下精度不高，容易失效[3].随着计算流体力学(CFD，Computational Fluid Dynamics)、计算结构力学(CSD，Computational Structure Dynamics)以及计算机技术的发展，高精度的CFD/CSD耦合求解[4]成为气动弹性问题的有效研究方法:流场与结构场在时域内各自推进，在交界面上通过传递压力、位移等参数实现耦合计算.CFD求解完全的三维非定常流场，由流场结果直接计算任意时刻作用在弹性体表面的气动力分布，精度高、适用范围广;CSD求解也可考虑结构的几何非线性及材料非线性，使计算结果更加真实.

在CFD/CSD耦合求解过程中，由于结构变形破坏了流场的拓扑结构，使得流场的计算网格在每个物理时间步内需要更新;由于流场网格与结构网格在界面上往往不匹配，需要通过一定的插值方法进行参数的传递.可见，除了CFD与CSD求解器外，动网格与界面数据传递方法也是气动弹性计算的关键技术.动网格算法要求高的计算效率以及强的适应大变形能力;界面数据传递方法要求精确插值以及满足运动、动力连续和能量守恒等条件[5]，任何一个环节出现问题都会影响耦合计算的精度和效率.

对于动网格问题，Batina[6]首先提出线弹性模型，之后由 Anderson[7]及 Farhat等[8]进行了改进和补充，极大地提高了其计算能力，使得弹簧比拟方法成为当前最常用的动网格方法之一.然而繁琐的数据结构以及计算效率、适应大变形能力等因素制约了弹簧比拟法在气动弹性问题中的应用.对于界面数据传递问题，无限平板样条法、多重二次曲面-双调和法、常体积转换法等方法被广泛地应用于气动弹性计算[5]，取得了一系列成果，但寻求高精度、高效率的界面数据传递方法仍是气动弹性计算发展的重点问题.

近年来，以径向基函数(RBF，Radial Basis Function)为基础的插值方法引起了数学界的广泛关注，Frank[9]在大量实例中对各种散乱插值方法进行了对比，发现径向基函数插值的效果最佳.径向基函数插值方法已经被广泛地应用于工程实际问题之中，Beckert[10]最早利用 RBF 插值处理机翼变形问题，Allen等[11]进一步将其用于网格运动，取得了理想的效果.

国内，苏波等[12]分析了径向基函数应用于流固耦合界面插值的可行性，并进行了任一时间步上的流固交互作用分析;林言中等[13]将RBF插值用于网格运动，并将其与弹簧比拟方法进行对比，说明了RBF方法在计算效率和适应大变形能力上的优势.

本文将基于径向基函数插值的动网格及界面数据传递方法用于CFD/CSD耦合计算，通过对AGARD 445.6机翼颤振问题的非定常数值模拟验证了该方法的可行性.径向基函数插值用于界面数据传递能够充分地保证运动连续、动力连续及能量守恒条件且兼有整体插值与局部插值方法的特点;用于动网格技术有着数据结构简单、适应大变形能力强等优势，为气动弹性问题的计算提供了非常有力的技术支撑.

1 数值方法

1.1 流动控制方程及数值方法

动网格上，积分形式的可压缩欧拉方程[14]可写成以下形式:

式中，W为守恒变量;F为对流通量，即

式中，V为流体质点的速度矢量;Vw为网格运动速度矢量;Ix，Iy，Iz分别为惯性系下沿坐标轴方向的单位矢量.

对式(1)采用有限体积法以一定的差分格式进行空间离散后，可得如下半离散方程:

本文的空间离散方案采用二阶精度的Van Leer格式.对式(3)可采用双时间步法求解.

1.2 结构动力学方程及数值方法

瞬态动力学平衡方程[15]为

式中，[M]为质量矩阵;[C]为阻尼矩阵;[K]为刚度矩阵;{F(t)}为随时间变化的外载荷;分别为节点加速度、速度及位移向量.采用有限元分析中的Newmark方法对式(4)进行求解，可得到结构各时刻的位移和应力分布.阻尼矩阵采用瑞利阻尼:

1.3 CFD/CSD耦合计算方法

当前针对CFD/CSD耦合已经发展了全耦合、松耦合及紧耦合3种方式，分别对应不同的时间推进方案.全耦合方案具有很高的时间精度，但非线性程度高，难以应用于复杂工程问题;松耦合方法计算复杂度低，有利于实现程序模块化，但仅有一阶时间精度;紧耦合方法可取较大的时间步长，但为满足时间精度需反复内迭代，计算耗时.本文采用二阶时间精度的松耦合格式，具体计算过程参见文献[16].

2 径向基函数插值

2.1 径向基函数及插值方法

径向基函数以空间距离 x为基本变量，通常可分为全局(global)函数、局部(local)函数和紧支(compact)函数[17]3 类，如表 1 所示.紧支函数在其中心距离达到某一特定值sr(称为紧支半径support radius)后函数值恒为0.表1中紧支函数的变量ζ为中心距离与紧支半径之比，即ζ=，当 ζ＞1 时，函数值为0.

表1 常用的RBF基函数

令d维欧氏空间中的一组中心点X={xc1，xc2，…，xcn}对应的标量值为 g={gc1，gc2，…，gcn}，则对于径向基函数存在连续函数:

通过所有中心点，确定插值系数γ的过程即为径向基函数插值.其中，为欧氏距离.对于三维空间而言，可直接用r表示为

2.2 径向基函数插值用于动网格

针对气动弹性计算中的动网格问题，可直接应用式(6)确定RBF插值函数.此时中心点即为耦合界面上流场区域的网格点，而相应的标量值为流体网格点3个方向上的位移分量Dx，Dy，Dz.插值系数可由以下线性系统确定:

式中n阶方阵M为插值矩阵，其元素为

矩阵M对称正定，可用共轭梯度法[18]快速求解.将插值系数 γx，γy，γz代入式(6)可得 3 个坐标分量的插值函数.将流场计算域内每个网格点的坐标代入RBF插值函数即可得到内场网格点的位移，从而实现整个计算域内的网格移动:

在气动弹性计算过程中，每个时间步的网格更新都可从零时刻的网格出发，这样位移插值系数只需计算一次，极大地提高了计算效率.需要指出的是，对于流体区域的网格运动问题，耦合界面的网格运动仅对界面附近区域的网格有较大影响，对于那些远离界面的网格点可认为位移很小甚至不运动.紧支基函数可通过调整紧支半径来控制中心点的影响区域，此特点与动网格的要求相一致，且求解方程组的系数矩阵呈稀疏、带状分布，对于存储量和计算效率都有着极强的适应能力，因此在动网格计算中更倾向于使用紧支型基函数.

2.3 径向基函数插值用于界面数据传递

针对气动弹性计算中的界面数据传递问题，中心点变为耦合界面上结构区域的网格点，而标量值则为结构网格点3个方向上的位移分量.为了使RBF插值函数满足平移及旋转不变性，通常需在式(6)的基础上加上d维低阶多项式p(x)，有

对于小于等于2阶的正定或条件正定基函数，p(x)可以采用如下线性多项式:

式(11)的2个定解条件为

式中q(x)为满足deg(q(x))≤deg(p(x))的所有多项式，由式(12)可知，q为阶数不超过1的多项式，即 q(x)的基函数为{1，x，y，z}.

插值系数 γi及多项式系数 β0，β1，β2，β3可由以下线性系统确定:

式中矩阵M元素同式(9)，矩阵P为

由于式(15)有零块，不宜直接求解线性方程组，可通过矩阵运算求得

由式(11)、式(17)及式(18)可建立流体网格与结构网格在界面上的位移传递关系:

得到位移传递矩阵H后，可根据能量守恒原理传递表面压力:

在气动弹性计算中，若整个过程结构变形不是很大，则可认为传递矩阵H基本不变，在初始时刻计算一次即可.对于基函数的选择，若采用全局函数，则RBF插值是一种整体插值方法;若采用紧支函数，则当紧支半径较小时，RBF插值呈现局部插值的特点.

3 AGARD 445.6机翼颤振算例

AGARD 445.6机翼颤振风洞实验是由美国国家航天局兰利研究中心在跨音速动态风洞中完成的，有着较为完备的风洞实验数据[19].目前，AGARD 445.6机翼软模型已经成为国际上跨音速气动弹性程序考核的标准算例，通过对其颤振问题的计算来考察RBF方法在气动弹性计算中的应用.AGARD 445.6机翼几何参数为:展弦比1.644，梢根比 0.659 2，根弦长 0.558 7 m，1/4 弦线后掠角为45°，沿流向翼型为NACA 65A004.流场计算采用结构化网格，网格总数746 916.机翼几何尺寸及对称面上的网格分布如图1所示.

图1 AGARD 445.6机翼几何参数及对称面网格

进行结构有限元计算时，机翼根部为固支边界条件.首先对机翼进行模态分析，各阶固有频率与实验数据对比如表2所示.计算结果与实验数据吻合较好，本文建立的机翼模型能够正确模拟实验模型的材料特性.

表2 AGARD 445.6机翼模态固有频率 Hz

首先取来流马赫数为0.499，静温297.2 K，攻角为0°的工况进行计算，计算初始时刻在翼尖处给一个微小法向初速度作为初始扰动.对于机翼颤振的气动弹性问题，若来流速度小于临界颤振速度，则振动收敛;若来流速度恰好等于临界颤振速度，则机翼呈现等幅谐振的状态;若来流速度大于临界颤振速度，则振动发散.通过调整无量纲来流速度Vf观察机翼的振动情况，观测对象为翼梢后缘点的法向位移.

无量纲速度定义为

式中，V∞为实际来流速度;bs为参考长度，一般定义为机翼的半根弦长;ωα为参考频率，该问题定义为机翼的一阶扭转频率;μ为质量比，定义为

式中，m为机翼质量;ρ∞为来流密度;V为体积系数.

由图2～图4的计算结果可以看出，随着来流无量纲速度的增加，机翼的动态响应曲线由收敛到发散，在Vf=0.455时响应曲线为等幅振动，此时机翼发生颤振，此速度为该马赫数下机翼的颤振速度，颤振频率为18.2 Hz.文献[19]中实验得到的机翼颤振速度为0.45，颤振频率20.4 Hz，与本文计算结果较为吻合.

图2 Vf=0.432，Ma=0.499时机翼的动态响应曲线

图3 Vf=0.455，Ma=0.499时机翼的动态响应曲线

由实验报告[19]可知，机翼跨声速颤振计算的马赫数范围为0.499～1.141，本文取马赫数为0.499，0.678，0.901，0.960 和 1.072 这 5 个状态进行颤振点的预测，计算结果与实验数据对比如图5所示.

图4 Vf=0.477，Ma=0.499时机翼的动态响应曲线

图5 机翼临界颤振速度边界

由图5可见，CFD/CSD耦合计算捕捉到了跨声速“凹坑”现象，计算的临界颤振速度在各来流马赫数状态下均与实验数据符合较好，对AGARD 445.6机翼的颤振计算结果表明径向基函数方法可以有效地应用于气动弹性计算之中.

4 结论

本文将径向基函数插值应用于气动弹性计算的两大关键技术:动网格及界面数据传递，介绍了径向基函数的基本概念及其用于气动弹性计算的数值方法.采用二阶时间精度的松耦合方案，对AGARD 445.6机翼颤振问题进行了CFD/CSD耦合计算.在较宽泛的马赫数范围内捕捉到了临界颤振速度及跨声速凹坑现象，计算结果与实验数据符合良好.结果表明，径向基函数插值方法可以有效地应用于气动弹性计算，为气动弹性问题研究的发展提供了有力的技术途径.

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，Farhat 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]苏波，石启印，钱若军.径向基函数应用于流固耦合分析初探[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]林言中，陈兵，徐旭.径向基函数插值方法在动网格技术中的应用[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].北京:冶金工业出版社，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