常用的外部气动力一般包括3种：高网格密度的CFD气动网格压力、中等网格密度的高阶面元法气动网格压力和低密度无网格的风洞试验气动压力数据点。同时，现有成熟的快速弹性飞行载荷分析方法包括高阶面元法和低阶面元法。因而在实际工程中，飞行载荷映射过程中包括3种网格尺度的对应关系：密到疏、疏密相当和疏到密。而且还存在网格交界面不完全匹配的情况，如将三维气动力等效映射至二维平板面元的情况。单纯使用传统的方法难以在上述复杂情况下实现气动力的等效映射。
飞行载荷分析中外部气动力的等效映射与流体分析中混合重叠网格中网格之间的流体信息传递类似，其传递方法一般包括两大类[4]：①守恒型方法，此类方法精度高，能保证气动力传递过程中力守恒，对网格相对尺度依赖小，但通常涉及网格求交运算，处理过程复杂；②非守恒型方法，此类方法一般使用插值方法，实现简单，但传递过程仅能保证气动力的分布一致，无法保证物理量的守恒。
常用的守恒型方法包括嵌入网格方法[5]、分段交界面方法[6]、守恒重映方法[7]、总体守恒隐式插值方法[8]，及之后Farrell、Maddison[9]和徐春光等[10]发展的超网格等方法。但这些方法都需要对原有网格进行重构或求交，虽然能严格保证插值守恒，但计算较为复杂。
非守恒型插值方法中，田书玲[11]将单元标准化，实现了基于非结构动态网格的Lagrange插值方法；Allen和Rendall[12]将径向基插值方法同时应用于流固耦合数据传递和动网格变形的计算中，并与多种方法的结果进行了对比；林言中等[13]将径向基插值方法应用于AGARD445.6机翼的流固耦合计算，并与实验结果进行了对比；孙岩等[14]基于径向基函数与混合背景网格实现了大变形时位移的精确插值。将这些非守恒插值应用于飞行载荷分析的气动压力时，虽然能获得较好的气动压力分布，但不能保证总体载荷的严格守恒。
综上所述，传统守恒型方法虽然能严格保证插值守恒，但计算复杂；非守恒型方法计算效率高，但不能保证载荷严格守恒。针对这一情况，本文利用等式约束保证积分载荷相等，利用基于加入线性基的径向基函数(RPIM)的最小二乘保证气动压力的分布特征，将气动载荷等效传递问题转化为一个带等式约束的二次规划问题。从而在非守恒方法的基础上构造出新的守恒型方法，兼顾了2种方法的优点。该方法在牺牲了一定气动压力分布精度的基础上，严格保证了总积分载荷的相等，同时也继承了非守恒型插值方法的高效性。并以某大展弦比机翼作为算例，将三维机翼上的气动力等效地映射至二维平板面元，验证了本文方法的有效性。
1 径向基插值方法 添加多项式的径向基插值函数可表示为[15]
 | (1) |
式中：x为插值面上任意一点的坐标，u(x)为该点的函数值；Ri(x)为径向基函数(RBF)；n为已知点的个数；pj(x)为多项式基函数；m为多项式基函数的个数，如m=0，则为单纯的RBF插值，否则为添加了m个多项式基函数的RPIM插值；ai和bj为待定系数；R(x)、a、P(x)和b分别为Ri(x)、ai、pj(x)和bi的矩阵和向量形式。
在实际应用径向基插值方法时，一般会使用添加线性基的RPIM插值方法，即P(x)=[1x]，以改善矩阵的条件数，并保证插值结果的线性再生性。
设有n个已知点的广义坐标x0T，这些点两两之间构成径向基函数矩阵R0和多项式基函数矩阵Pm，U0为已知点的函数值。再根据正交约束条件添加m个约束条件(左侧)，从而推出RPIM插值矩阵(右侧)为
 | (2) |
式中：H为已知点函数值U0的增广矩阵; G为R0和Pm构成的径向基插值矩阵。
对线性方程组(2)求解得到待定系数a0，再代入式(1)，即可得到待插值点的函数值，而这一函数值在本文的飞行载荷外部气动力等效映射问题中代表气动压力系数。
常用的径向基插值函数[16]Ri(x)包括高斯型、对数型(又称薄板样条)、复合2次型、逆复合2次型，以及Wendland[17]提出的一系列紧支型，如表 1所示，α和c为形状参数。
表 1 径向基函数 Table 1 Radial basis function
| 类型 | R(x)定义式 |
| 高斯型(G) | e－α||x||2 |
| 对数型(TPS) | ||x||2ln||x|| |
| 复合2次型 | (c2+||x||2)1/2 |
| 逆复合2次型 | 1/(c2+||x||2)1/2 |
| Wendland’s C0 | (1－||x||)+2 |
| Wendland’s C2 | (1－||x||)+4(4||x||+1) |
| Wendland’s C4 | (1－||x||)+6(35||x||2+18||x||+3) |
表选项
径向基函数中, x表示待求点x到已知点x0之间的欧氏距离，设其广义坐标共有q维，
 | (3) |
设||x||=||x||/δ为紧支距离，其中：δ为形状参数，用于定义紧支半径。则紧支含义如下：
 |
Smith等[18]在流固耦合分析中对各种类型的插值方法以及各种形式径向基函数进行了测试，结果表明对数型(TPS)的径向基函数在流体问题中的插值精度表现最好。因而本文在RPIM插值中选用对数型径向基函数，即
 | (4) |
2 等效映射方法 基于RPIM的插值型气动力映射方法虽然能获得更为精确的压力系数分布，但在飞行载荷分析中，力和力矩的守恒也十分关键。要保证力和力矩的守恒，势必要对原插值方法得到的结果进行调整, 而此问题正好契合数值优化中的等式约束二次规划问题。
2.1 等式约束二次规划问题 等式约束二次规划问题可以描述为[19](本文对原文献的变量符号和正负号有所调整)
 | (5) |
式中：G∈Rn×n为对称阵；h∈Rn；A∈Rn×m，m≤n；f∈Rm，一般要求A列满秩。
采用Lagrange方法求解等式约束二次规划问题(5)，可得该问题的Lagrange函数为
 | (6) |
式中：λ为Lagrange乘子。则式(6)的KKT(Karush-Kuhn-Tucker)条件为
 | (7) |
写成矩阵形式可表示为
 | (8) |
求解式(8)所形成的线性方程组，即能得到等式约束二次规划问题的最优解x0。
2.2 气动力等效映射问题的转化 设外部气动力等效中的气动压力点个数为n，已知的气动压力系数列向量为Cp0；待映射计算气动面网格数为l，坐标值为xT=[x1 x2…xl]q×l，q为坐标值的维度，待映射的气动压力分布为Cp。本文采用添加线性多项式的径向基函数进行插值，其中线性基为
 | (9) |
根据式(1)和式(2)可以得到基于RPIM的待映射气动压力分布Cp和基于外部气动力数据的RPIM插值矩阵G定义为
 | (10) |
 | (11) |
以二维情况为例，此时积分载荷系数F可用基于RPIM插值矩阵的系数a0表示。此时q=2，xT=[xi, yi]，对应积分载荷系数F包含1个力系数和2个力矩系数。则由面积系数[s1 s2 … sl]对气动网格右手坐标系原点积分的力和力矩系数积分矩阵
 | (12) |
则对应的积分载荷系数F可表示成
 | (13) |
原问题欲保证积分载荷等效和保持原气动压力分布，即寻求系数矩阵aL满足式(13)，并满足式(11)的最小二乘要求。而式(11)的最小二乘可表示为
 | (14) |
式中：W为加权系数对角阵，本文暂时取为单位阵。又RPIM形成的插值矩阵G为对称阵，则有G=GT，及aLTGTH=(aLTGTH)T=HTGaL，再展开式(14)等号右边可得
 | (15) |
另外，由式(13)可以得到积分载荷系数守恒的等式约束条件为
 | (16) |
由2.1节可知，式(15)和式(16)构成的等式约束二次规划问题的KKT条件为
 | (17) |
求解式(17)的线性方程组，可得满足原问题要求的插值矩阵系数aL，再将其代入式(10)，即可得到能够保证积分载荷等效并保持原气动压力分布的计算气动网格气动压力分布Cp。
3 算例与计算结果 本文选用某大展弦比机翼的上表面作为原气动网格(三维)，并以NASTRAN飞行载荷分析中的低阶平板面元网格(二维)作为待映射计算气动网格(见图 1和图 2)。将Ma=0.7，迎角1°下风洞测压数据处理至原气动网格后的气动压力分布作为外部气动力数据，分别用RPIM插值方法和本文提出的基于RPIM的二次规划方法进行气动力映射。
 |
| 图 1 原三维气动网格 Fig. 1 Original three-dimensional aerodynamic mesh |
| 图选项 |
 |
| 图 2 待映射二维平板面元网格 Fig. 2 Two-dimensional panel mesh for mapping |
| 图选项 |
原气动网格和经2种不同气动力映射方法得到的压力系数Cp分布云图如图 3所示。
 |
| 图 3 原气动网格、RPIM插值方法及二次规划方法得到的压力系数分布 Fig. 3 Pressure coefficient distribution obtained by original aerodynamic mesh, RPIM interpolation method and quadratic programming method |
| 图选项 |
图 4给出了从翼根至翼尖20%、40%、60%和80%展向位置3种方法的截面压力系数结果，xc为弦向位置。从气动压力分布结果来，看基于RPIM插值的二次规划方法要略逊于RPIM插值方法，但保持了原有结果的分布特征。这与图 3(a)原气动网格的压力系数分布的结果是一致的。
 |
| 图 4 20%、40%、60%和80%展向位置截面压力系数 Fig. 4 Section pressure coefficient in 20%, 40%, 60% and 80% spanwise position |
| 图选项 |
通过图 3(a)原气动网格和图 4气动压力分布和各展向位置截面压力系数结果的对比，以及表 2映射前后总积分载荷系数可知，虽然非守恒的RPIM插值方法得到的气动压力分布能够更好地拟合原气动网格，但是得到的总积分载荷系数与原气动网格存在误差；而本文提出的基于RPIM的二次规划方法不仅能够保证映射后的总积分载荷系数与原模型相等，同时还能保持原模型的气动压力分布特征。
表 2 映射前后总积分载荷系数 Table 2 Total integral load coefficient before and after mapping
| 系数 | 原气动 网格 | RPIM插值方法 (误差) | 二次规划方法 (误差) |
| Cz | 0.149 08 | 0.147 0(1.4%) | 0.149 08(0%) |
| CMx | 0.027 657 5 | 0.027 081(2.08%) | 0.027 657 5(0%) |
| CMy | -0.060 124 | -0.059 259(1.44%) | -0.060 124(0%) |
表选项
表 2的结果基于翼根处的坐标系积分得到，其中Cz为升力系数，CMy为扭矩系数，CMx为弯矩系数。为了对比映射前后的局部积分载荷情况，图 5给出了由翼尖(100%)积分至翼根(0)各展向位置yb的积分载荷系数。积分参考坐标系见图 6。
 |
| 图 5 升力系数、扭矩系数和弯矩系数沿展向积分结果 Fig. 5 Integral lift coefficient, torque coefficient and bending moment coefficient along spanwise direction |
| 图选项 |
 |
| 图 6 积分载荷坐标系 Fig. 6 Coordinate system for integrating loads |
| 图选项 |
由升力系数Cz的展向积分结果可知，除了在20%展向位置处二次规划方法的结果优于RPIM插值方法之外，其他位置的结果RPIM插值方法更接近于原气动网格。弯矩系数CMx的结果二次规划方法更优，而扭矩系数的结果RPIM插值方法更优。这说明二次规划方法为了保证全局总积分载荷系数的守恒，虽然一定程度上损失了压力系数分布的精度，但是在飞行载荷分析更为关
注的某些局部积分载荷分布结果，相比RPIM插值方法却更接近原气动网格。
另外，在计算速度上，基于RPIM的二次规划方法的插值矩阵相比RPIM插值方法只多出了与等式约束个数相等的阶数，因而效率相当。而常用的守恒型插值方法由于存在复杂的搜索和网格重构过程，在基于相同结点数量的情况下，为达到相同的计算精度，其效率只有非守恒型插值方法的四分之一[20]。另一方面，RPIM插值方法还可以通过合理选点进一步提高效率。
4 结论 1) 提出了基于RPIM的守恒型二次规划方法，在牺牲一定的气动压力分布精度的情况下，兼顾了守恒型和非守恒型两者的优点。既能够严格保证气动力映射前后总积分载荷的相等，也保持了原气动压力的分布特征，同时也继承了非守恒型插值的高效性。
2) 所提方法能够适用于其他各类复杂情况的气动力等效映射，同时还可拓展至其他学科类似的数据映射问题中。
参考文献
| [1] | 《飞机设计手册》总编委会. 飞机设计手册第6册:气动设计[M]. 北京: 航空工业出版社, 2001: 49-186. General Editorial Board of Aircraft Design Manual. Aircraft design manual(6th volume):Aerodynamic design[M]. Beijing: Aviation Industry Press, 2001: 49-186. (in Chinese) |
| [2] | 万志强, 邓立东, 杨超, 等. 基于非线性试验气动力的飞机静气动弹性响应分析[J]. 航空学报, 2005, 26(4): 439-445. WAN Z Q, DENG L D, YANG C, et al. Aircraft static aeroelastic response analysis based on nonlinear experimental aerodynamic data[J]. Acta Aeronautica et Astronautica Sinica, 2005, 26(4): 439-445. DOI:10.3321/j.issn:1000-6893.2005.04.012 (in Chinese) |
| [3] | 安伟刚, 梁生云, 陈殿宇. 一种局部动态数据交换方法在流固耦合分析中的应用[J]. 航空学报, 2013, 34(3): 541-546. AN W G, LIANG S Y, CHEN D Y. Local dynamic data exchange in fluid structure interaction analysis[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(3): 541-546. (in Chinese) |
| [4] | 黄宇, 阎超, 王文, 等. 混合重叠网格插值方法的改进及应用[J]. 北京航空航天大学学报, 2017, 43(2): 285-292. HUANG Y, YAN C, WANG W, et al. An improved interpolation method for hybrid overset grid and its application[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(2): 285-292. (in Chinese) |
| [5] | YAO Z, LIOU M S.Progress in the three-dimensional dragon grid scheme: AIAA-2001-2540[R].Reston: AIAA, 2001. |
| [6] | XU K L, SUN G, CAI J Y.On interface conservative attributions for computations of the complex flows of high-lift system based on chimera technique[C]//International Conference on Computer and Automation Engineering.Piscataway, NJ: IEEE Press, 2010: 279-283. |
| [7] | 张宇飞, 陈海昕, 符松. 基于高阶守恒重映对窗口嵌入技术的改进[J]. 计算物理, 2011, 28(2): 167-173. ZHANG Y F, CHEN H X, FU S. Improvement on window embedment technology with high order conservative remapping[J]. Chinese Journal of Computational Physics, 2011, 28(2): 167-173. DOI:10.3969/j.issn.1001-246X.2011.02.002 (in Chinese) |
| [8] | XIANG Z.An implicit and globally conservative unstructured chimera grid method: AIAA-2011-777[R].Reston: AIAA, 2011. |
| [9] | FARRELL P E, MADDISON J R. Conservative interpolation between volume meshes by local galerkin projection[J]. Computer Methods in Applied Mechanic and Engineering, 2011, 200(1): 89-100. |
| [10] | 徐春光, 董海波, 刘君. 基于单元相交的混合网格精确守恒插值方法[J]. 爆炸与冲击, 2016, 36(3): 305-312. XU C G, DONG H B, LIU J. An accurate conservative interpolation method for the mixed grid based on the intersection of grid cells[J]. Explosion and Shock Waves, 2016, 36(3): 305-312. (in Chinese) |
| [11] | 田书玲.基于非结构网格方法的重叠网格算法研究[D].南京: 南京航空航天大学, 2008: 63. TIAN S L.Investigation of overset unstructured grids algorithm[D].Nanjing: Nanjing University of Aeronautics and Astronautics, 2008: 63(in Chinese). |
| [12] | ALLEN C B, RENDALL T C S. Unified fluid-structure interpolation and mesh motion using radial basis functions[J]. International Journal for Numerical Methods in Engineering, 2008, 74: 1519-1559. DOI:10.1002/nme.2219 |
| [13] | 林言中, 陈兵, 徐旭. 基于径向基函数插值的气动弹性计算方法[J]. 北京航空航天大学学报, 2014, 40(7): 953-958. LIN Y Z, CHEN B, XU X. Numerical method of aeroelasticity based on radial basis function interpolation[J]. Journal of Beijing University of Aeronautics and Astronautics, 2014, 40(7): 953-958. (in Chinese) |
| [14] | 孙岩, 孟德虹, 王运涛, 等. 基于径向基函数与混合背景网格的动态网格变形方法[J]. 航空学报, 2016, 37(5): 1462-1472. SUN Y, MENG D H, WANG Y T, et al. Dynamic grid deformation method based on radial basis function and hybrid background grid[J]. Acta Aeronautica et Astronautica Sinic, 2016, 37(5): 1462-1472. (in Chinese) |
| [15] | LIU G R, GU Y T.无网格法理论及程序设计[M].王建明, 周学军, 译.济南: 山东大学出版社, 2007: 45-48. LIU G R, GU Y T.An Introduction to meshfree methods and their programming[M].WANG J M, ZHOU X J, translated.Jinan: Shandong University Press, 2007: 45-48(in Chinese). |
| [16] | 陈文, 傅卓佳, 魏星. 科学与工程计算中的径向基函数方法[M]. 北京: 科学出版社, 2014: 11-13. CHEN W, FU Z J, WEI X. Radial basis function method in scientific and engineering computation[M]. Beijing: Science Press, 2014: 11-13. (in Chinese) |
| [17] | WENDLAND H. Piecewise polynomial, positive definite and compactly supported radial basis functions of minimal degree[J]. Advances Computational Mathematics, 1995, 4(1): 389-396. DOI:10.1007/BF02123482 |
| [18] | SMITH M J, HODGES D H, CESNIK C E S. Evaluation of computational algorithms suitable for fluid-structure interactions[J]. Journal of Aircraft, 2000, 37(2): 282-294. DOI:10.2514/2.2592 |
| [19] | 高立. 数值最优化方法[M]. 北京: 北京大学出版社, 2014: 217-226. GAO L. Numerical optimization method[M]. Beijing: Peking University Press, 2014: 217-226. (in Chinese) |
| [20] | BOER A D, ZUIJLEN A H V, BIJL H. Review of coupling methods for non-matching meshes[J]. Computer Methods in Applied Mechanics and Engineering, 2007, 196(8): 1515-1525. DOI:10.1016/j.cma.2006.03.017 |
