刘智侃，刘深深,2,*，刘骁，曾磊，代光月

1.中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000 2.中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000

临近空间类高超声速飞行器由于其独特的飞行环境及特点，具有十分重要的应用价值，已成为当前高超声速飞行器研究的重要发展方向。该类飞行器需要在稠密大气中进行长时间的高超声速飞行，因而面临严酷的气动加热问题。长时间气动加热带来的温升导致结构受热后产生静态变形，而结构变形进一步导致流场气动加热的部分改变，造成表面热流、温度场、结构应力/应变场之间发生复杂的耦合作用，从而引发综合热效应[1-3]。独立分散的热环境、热响应和热应力分析方法人为割裂了实际存在的多物理场共同作用的复杂过程，对存在力/热/结构耦合严重的情况并不适用[4-5]，因此在当前针对高超声速飞行器的热防护设计和优化中必须综合考虑气动力/热/结构的多场耦合计算问题。

当前多物理场耦合数值计算问题的求解主要是松耦合方法[6]，它的基本思想是将紧耦合联立求解的问题简化为递推求解的过程：在极小的时间步内各个模块分别采用各自适用的求解器进行求解，通过不断地在流体域和固体域交界面上进行数据传递来实现实时耦合。松耦合方法可以复用现有的发展比较成熟的气动力、气动热、结构热响应等求解器，同时可以通过控制交互时间步的大小来简化计算开销，因而对复杂的外形具有很好的适用性，应用较为广泛[7-8]。在松耦合方案中，不同物理场之间的数据传递方法是驱动和实现耦合计算的关键桥梁，其精度和效率共同决定了耦合计算能否在可承受的时间代价下获得真实反映物理耦合过程的结果，因而是实现高效高精度耦合数值模拟的关键环节。

径向基函数(Radial Basis Function, RBF)因具有形式简单，不依赖于求解器离散格式、网格拓扑结构，易于编程且精度较高，能保证所传递的力载荷的守恒性[9]等特点在多场耦合计算问题中得到了广泛应用。Rendall等[10-12]开展了一系列研究，将其应用于高超声速机翼气弹问题分析，表明其对CFD/CSD数据插值及大变形位移具有良好的适用性；郭中州等[13]发展的高效K近邻-径向基函数(KNN-RBF)动网格方法提高了网格变形的计算效率和对复杂外形的适用性；刘君等[14]利用RBFs-MSA混合算法实现任意几何形状的计算网格变形，能够有效满足黏性流动模拟需要的异向加密需求；王刚等[15]提出了基于贪心法的点选取方法以消除邻近点人为选择对插值精度的影响；为了减少RBF插值过程中支撑点数目，魏其等[16]提出峰值选择法以减少迭代次数提高选点效率。但采用径向基函数插值时存在的问题是径向基本身是各向同性的，更适用于描述各向同性的物理量分布。而具备复杂外形的高超声速飞行器物理量分布本身与局部流动特征高度相关且梯度变化剧烈，例如飞行器头部、前缘、翼舵干扰区及机身大面积区域压力、热流等物理量分布差异巨大；除去上述因素外，由于复杂外形网格的划分会在外形曲率或物理量变化剧烈区域进行加密处理因而网格点分布也是各向异性的，人为引入了不同网格点分布的各向异性因素，由此采用各向同性基函数进行高度各向异性物理量空间插值时存在缺陷，造成了精度及守恒性提升瓶颈。

为了解决上述问题，需要对径向基函数进行各向异性修正。考虑减少人为网格划分的影响，刘深深等[17]提出的一种基于网格划分几何尺度归一化方法的各向异性修正方法取得了较好的效果，但该项工作未考虑物理量自身分布特征的影响。为了解决物理量各向异性分布对插值精度带来的影响，本文考虑将物理量梯度的分布表征因子引入径向基计算，提出了一种基于物理量梯度进行径向基各方向自适应因子修正的新方法。采用高超声速舵面、翼身组合体算例对该方法在单次数据传递中的可行性及效果进行了验证分析，同时基于该方法对压缩拐角外形开展了气动力/热/结构多场耦合数值模拟，并将耦合计算结果与风洞试验结果进行了对比分析，对该方法在实际耦合问题计算长时间高频次数据传递中的适用性及可靠性开展了分析验证研究。

1 建立基于物理量梯度修正的改进模型

1.1 径向基函数数据传递基本原理

在d维欧几里得空间上给定函数φ:R+→R，对于空间中n个不同的点{xi,gi|i=1,2,…,n}∈Rd⊗R，在这些固定点上的标量值为{g1,g2,…,gn}。根据这些点需要确定一个函数：

(1)

使其满足插值条件：

(2)

(3)

(4)

其中，添加的多项式为

p(x)=λ1+λ2x+λ3y+λ4z

(5)

需要添加的定解条件[18]为

(6)

依据基函数的定义和定解条件便可以推导传递矩阵，完成数据的传递。按照基函数作用域的范围可以将RBF分为全局域RBF和紧支RBF，全局域RBF指每个RBF影响区域为整个定义域，如薄板样条插值(Thin Plate Spline, TPS)、多重二次函数插值(Multiquadric, MQ)等基函数，而紧支RBF的影响区域为指定的影响半径，如Wenland定义的C2紧支基函数[17]。这3种基函数是目前多场耦合计算中最为常用的3种基函数，本文研究主要针对上述基函数进行，具体表达式如表1所示。

表1 3种径向基函数

1.2 基于物理量梯度修正RBF方法

吴宗敏[18]指出理论上依据式(3)计算的欧式距离r是各向同性的，故所述的RBF插值方法更适用于物理上各向同性的问题。而实际流动物理量在各个方向的分布并不是一致的，如图1所示,在机翼前缘附近，由于流动特征及网格划分的影响，热流变化在x和y方向变化十分剧烈，在z方向变化较为缓慢，在固定空间半径范围内对于x和y方向的插值所选取的点在待插值点附近越集中插值精度越高，而物理量特征相差较大的D点等的引入会带来更多空间分布拟合误差，因此需要通过增大x和y方向的距离权重来排除物理量相差过大的点，同时相应的在物理量分布变化缓慢的z方向，需要通过减小该方向的距离权重来增加相应插值点E，以防止图中所示的插值点过少的情况。根据上述思想，可以基于物理量梯度对径向基计算中的x、y、z3方向进行系数缩比调控，以求在相同的插值半径选点范围内用来插值的点能够更具有代表性和聚集性，从而使插值的结果更能表征物理量实际分布特征，进而提高精度。

图1 机翼前缘热流分布及固定半径内选择的点

针对径向基函数中欧式距离r的计算进行如下变换[17]：

(7)

其基本思路是对外形基于物理量梯度分布进行适当的缩放使得网格点在不同方向上的分布更加均匀，其中c1、c2、c3由物理梯度决定，计算方法如下：

步骤1求解各个已知网格点上的物理量分布梯度。对于离散网格点的梯度求解本文采用的方法主要为最小二乘法。其基本原理是基于网格节点上物理量值的一阶Taylor级数展开近似得到

(8)

式中:U为物理量的值。将式(8)应用到与目标节点i相连接的所有节点，得到超定约束线性方程组：

(9)

式中：Δ(·)ij=(·)j-(·)i;M表示与目标节点i相连的网格节点数目；θj表示权系数，如面积、体积或距离等，本文的θj值均为1。通过求解上述超定方程组，可得物理量的梯度分布。

对于网格点物理量梯度的求解，需要用到网格点周围的连接点信息，由附近连接点的物理量值进行梯度的求解，对于网格节点间间距较大的节点该方法求得的梯度值精度较低。对于结构/非结构网格在初始插值时均需进行一次搜寻并存储相关连接信息，耦合计算中多次插值时不需要进行重复搜索，因而不会带来过多额外的计算开销。

步骤2求得所有已知节点在各个方向的梯度总和Gx、Gy、Gz，计算公式为

(10)

式中：N为原始网格点的个数。

(11)

式中：i=1,2,3分别代表x、y、z3个方向。在此基础上同时考虑几何外形的尺度归一，可以进一步提高插值算法的鲁棒性：

(12)

步骤4对所有已知网格点和待插值网格点的坐标进行基于物理量梯度的修正：

步骤5将缩放后的值代入式(3)中，得到更新后的rnew的表达式为

(13)

对比式(7)与式(13)可知：

(14)

式中:c1、c2、c3为基于物理量梯度分布得到的自适应修正各向异性因子。

通过分析式(3)与式(7)可知，对于给定的空间半径网格点的搜索是在一个球内，而梯度修正后的网格点的搜索则是在一个椭球内。对网格缩放后外形上某一点来说，距离r相同的所有点组成了一个球面，这等同于在原外形中距离为r的数据点组成了一个椭球面，即相同的r所包围的数据点坐标范围根据不同方向的物理量梯度分布进行了拉伸，变化剧烈的方向被拉长，变化平缓的方向被压缩，进而能够降低物理量梯度和几何外形的影响。

2 改进方法的基础算例验证

高超声速飞行器具有明显的气动力/热/结构的多场耦合现象，涉及到气动力、热流、温度、位移等物理场的传递。对于气动力、热流等载荷传递，不仅需要点对点的精确插值，还需要整个交界面的通量守恒。一般而言，气动力数据传递可以通过虚功原理来自动保证力学能量守恒，而对于热流的传递来说，只能通过更精确的插值来保证交界面上的通量守恒[19]。因而热流的传递是高超声速飞行器耦合问题插值中的难点，选择热流插值作为考核变量。根据前文对基于物理量梯度修正RBF方法的介绍，采用文献[17]中高超声速飞行器中的舵面结构和某一翼身组合体飞行器作为考核算例，对热流传递的准确性进行对比分析。根据以往的研究[20]，样条函数中的TPS方法、MQ方法和紧支C2基函数方法较为常用，本文将在此基础上分析物理量梯度修正前后热流插值效果。

2.1 高超声速飞行器舵面模型

对于高超声速飞行器来说，翼舵结构是其一大典型外形特征，且在飞行过程中气动力/热/结构耦合较为显著。因此，本文针对高超声速飞行器某一舵面进行了热流的插值传递计算与分析。舵面翼根长2.67 m，翼展长1.27 m，前缘后掠角为30°，前后缘半径均为7 mm。对该舵面进行网格划分，得到流体域网格和固体域网格，分别如图2、图3所示。其中流场表面网格节点数4 073个，单元数4 016个；结构表面网格节点数5 317个，单元数10 520个。

图2 舵面流体域网格

图3 舵面固体域网格

针对该翼面的外形特征，假设舵面上任意坐标数值分别为(x,y)，则表面热流分布表达式为

(15)

根据该热流表达式可以得到舵面的热流分布如图4所示。

图4 舵面流体域网格表面热流分布云图

热流的传递由TPS、MQ和紧支C2基函数6种插值方法(包含原始和梯度修正插值方案)从流体域网格到固体域网格进行传递，其中MQ的形态参数c取0.01，紧支基函数半径固定为算例外形3个方向最大尺度的0.3倍。同时依据表达式计算固体域网格点的表面热流作为标准值来衡量插值的精度，误差的计算公式为

(16)

图5给出了插值后固体域网格的热流分布云图。从图中可以看出，6种插值方法均能较好地捕捉到流场表面的热流分布特征，但是从舵面前缘等流动变化剧烈位置来看，原始的3种插值方法会形成间断性的不规则波纹；与之对比梯度修正后的3种插值方法在舵面前缘部位均未出现波纹，对于流场局部细节的捕捉要优于原始插值方法。其原因主要在于在舵面前缘位置，几何模型曲率较大，且物理量沿z方向梯度较大，导致舵面前缘流场网格呈现长条型，即沿y方向网格节点要比z方向稀疏得多，而固体域网格在两个方向的节点分布则相对均匀得多，因此对于这种从稀疏网格到密网格的插值容易产生数值振荡。

图5 舵面固体域网格热流分布云图

图6为6种插值方法的误差分布云图，从结果比较来看原始方法在前缘部位误差均较大，相对误差在10-2左右；梯度修正的插值方法整体插值效果要优于原始方法，在前缘部位的改进尤为明显。而从舵面的上表面来看其改进效果也相当明显，原始方法的上表面相对误差在10-3左右，梯度修正后的相对误差则主要集中于10-4左右，下降了约一个量级。

图6 舵面表面热流误差分布云图

图7给出了各方法热流插值过程中的最大误差和平均误差统计结果，很显然梯度修正的插值方法的平均误差、最大误差均比原始的插值误差要低，其中TPS方法取对数后的平均误差从-3.56 降低到-4.63，紧支C2基函数方法取对数后的平均误差从-3.90降低到-4.80。综合来看本文采用的改进方法要优于原始的插值方法。

图7 舵面6种插值方法平均误差和最大误差

MQ插值方法和紧支C2基函数插值方法中均存在形态参数，不同的形态参数对插值结果的影响是巨大的。对于不同算例中不合理的参数甚至可能导致插值结果误差过大不能进行下一步的计算。本文以MQ中的形态参数为例分析了原始插值方法和梯度修正后插值方法中形态参数对结果的影响。图8给出了形态参数c分别为0.001、0.01、0.1时插值结果误差对比云图。从图中可以发现原始的MQ方法在形态参数较小时误差较小，随着形态参数的增大，误差越来越大；在相同的形态参数取值下，采用梯度修正能够很好地降低插值误差，并且采用梯度修正后的MQ插值方法可以在较大的形态参数变化范围[0.001，0.05]内保持很好的精度，因而提高了在实际工程问题中的适用性。

图8 不同形态参数插值误差对比云图

2.2 高超声速翼身组合体模型

本节采用类X-37运载器[21]翼身组合体外形考察在复杂外形和流场分布情况下TPS、MQ和紧支C2基函数3种径向基函数在工程实用外形中的实用性和改进插值方法带来的效果。分别采用结构网格、非结构网格生成如图9、图10两套不同的网格，其中流体域网格节点数13 216,网格单元数11 702，固体域网格节点数5 704，网格单元数11 088。采用图9的网格得到了马赫数8.84、攻角0°、壁面温度300 K条件下的数值气动热结果，如图11所示。利用3种不同径向基函数分别采用梯度修正方法和原始方法完成热流的传递，其中MQ的形态参数c取0.001 1，紧支基函数半径固定为算例外形3个方向最大尺度的0.25倍。

图9 翼身组合体流体域网格

图10 翼身组合体固体域网格

图11 翼身组合体流体域网格数值计算热流分布

图12给出了采用原始插值方法和梯度修正插值方法得到的热流插值结果云图，从图12(a)～图12(d)中可以发现TPS和MQ方法能比较准确地捕捉到模型整体的热流分布特征，但是在机翼前缘部位原始MQ方法会出现一些不规则的条纹，而梯度修正的插值方法可以有效解决该类问题。梯度修正后的紧支C2基函数插值方法可以有效改善机身头部和机身中段处的插值精度，如图12(e)和图12(f)所示。从图中可以发现，头部区域的等值线不太光滑，部分区域毛刺较多，与原始热流等值线相比有所偏离，而梯度修正后的紧支C2基函数插值结果较为光滑，更接近原有的热流分布。经过分析发现，头部区域及机翼前缘是热流梯度较大区域，从图9和图10对比可知流体域的网格在此处的长细比较大，但是固体域的网格则相对均匀，因此基于物理量梯度修正的插值方法选择的点更加均匀插值效果也更好。上述结果表明基于物理量梯度修正的方法可以在头部及机翼前缘等流动变化剧烈的位置更精确捕捉到流动细节，实现对复杂物理量分布的精细刻画。

图12 翼身组合体固体域网格表面热流分布云图

由于固体域网格并没有真实的流场结果，为了对比每一个数据点的插值误差，将固体域上插值得到的热流再次采用相同的方法插值回流体域。误差计算公式与2.1节的公式一致，其中qinterpolation为从固体域插回流场的数据，qreal为原始热流分布。为了观察改进前后插值结果的误差对比，本文对误差分布取以10为底的对数，图13为原始插值方法和梯度修正插值方法的误差分布云图，从图中可以直观地看出，原始TPS方法的插值误差主要在机身中段和机翼前缘，紧支C2基函数方法的误差在机身头部区域也较大，而原始MQ方法在整个机身表面误差均较大。经过梯度修正后的插值方法不仅改善了前缘部位的插值，还大大改善了机身中段的插值效果。尤其是MQ方法和紧支C2基函数方法的改进效果均很明显，其中梯度修正后的紧支C2基函数插值效果甚至比梯度修正后的TPS、MQ方法效果更优，即用更少的点达到了全域基函数的效果。此外，经过测试发现当梯度修正方法的插值半径为算例外形3个方向最大尺度的0.03倍时，误差与原始插值方法的0.25倍得到的结果相近。经过统计分析发现对于原始网格点数为5 704，待插值网格点数为13 216的插值算例，梯度修正方法的平均选取点数约为181个，而插值半径为0.25倍的原始方法平均选取点数为1 859个，选点个数下降90%以上，插值时间由254.8 s下降到26.4 s。因此本文提出的梯度修正方法实现了精度和效率的双重兼顾。

图13 翼身组合体流体域网格表面热流误差分布云图

本文进一步对比了在复杂算例中基于物理量梯度修正的紧支C2基函数方法与文献[17]中尺度归一化的紧支C2基函数方法插值效果。其表达式为

(17)

其中：

(18)

式中:Lx、Ly、Lz分别代表3个方向的长度大小。其基本原理是减少飞行器本身3个方向尺度差异过大给插值带来的影响，主要思想是对飞行器外形在插值前进行尺度归一化处理。

图14和图15分别为二者的误差分布云图，从图中可知基于梯度修正的紧支C2基函数方法在头部区域的相对误差集中于-3.5，尺度归一化的误差集中于-2.5，说明前者的插值效果要优于后者。定量分析二者的相对误差可知尺度归一化的插值方法平均相对误差为-2.756，最大相对误差为0.181，而本文基于物理量梯度修正的插值方法平均相对误差为-3.046，最大相对误差为0.173，二者均比尺度归一化的数值要小。因此可知在复杂外形算例中梯度修正的插值方法要优于文献[17]中提出的尺度归一化方法。

图14 梯度修正方法误差云图

图15 尺度归一化方法误差云图

图16给出了改进前后3种插值方法的最大插值误差和平均插值误差。整体来看每种梯度修正后的方法均比原始的方法最大误差要小。此外，梯度修正后方法的平均误差均要低于原始的插值方法，其中紧支C2基函数的平均误差下降幅度最大，改进的插值效果最好。

图16 翼身组合体6种插值方法平均误差和最大误差

为了进一步定量分析各种插值方法的误差变化，图17给出了改进前后各种插值方法的误差累计百分比曲线图。横坐标-8代表相对误差区间为[-8,-7)，纵坐标给出了误差区间内点数之和占总点数的百分比，其他依次类推，如原始紧支C2基函数的第6个点坐标(-3,31.794)表示相对误差小于-2的点数占总点数的31.794%。从统计图来看，梯度修正后的方法曲线上升的速度要大于原始的方法，说明梯度修正后方法的误差均集中在误差较小的区域，如梯度修正后的紧支C2基函数方法误差在[-8,-1]区间内的点数占比达到了98%。而原始的紧支C2基函数方法在该区间内的点数只占到了86.6%，仍有14.4%的点的误差在区间[-1, 0]内。综合来看，对于复杂外形经过物理量梯度修正后的插值方法插值精度均有所提升，紧支C2基函数方法的改进效果最好，原始紧支C2基函数方法误差高于TPS、MQ方法，但梯度修正后的紧支C2基函数方法要优于TPS、MQ方法，其结果与图13云图中分析的结果一致。

图17 翼身组合体6种插值方法的热流误差统计

3 基于风洞试验数据的改进方法耦合算例验证

基于所提出的改进紧支C2基函数数据传递方法，开展了针对压缩拐角外形气动力/热/结构多场耦合风洞试验状态的耦合计算数值模拟，以验证该方法在长时间耦合计算中的适用性和精度。图18给出了该压缩拐角风洞试验模型尺寸[22]，其前缘半径3 mm，长633.40 mm，宽140 mm，第一压缩面角度为7.34°，第二压缩面角度为25.34°，形成18°压缩拐角。

图18 试验模型结构尺寸及示意图

该试验在中国空气动力研究与发展中心∅0.6 m 风洞进行。风洞试验状态参数等效为在高度23.5 km，马赫数5.5，攻角α=10°巡航飞行60 s[22]。试验模型材料选取高温合金IC6，具体参数见《中国航空材料手册》[23]。该试验重点关注并获取了气动力/热/结构耦合状态下的飞行外形变化位移数据。图19为通过特制LED光源照射模型并通过相机获得的初始状态(t=0 s)、3 s时刻、7 s时刻模型变形特性图。从图中可以看到，随着时间的推进，压缩面前体位置出现了明显的位移，时间越长，前体位移量越明显。

图19 典型时刻前景照明图

本文耦合计算采用FL-CAPTER软件平台内置的双向面/体自适应耦合计算策略[24]，结构应力/应变场计算时，在尾部螺栓孔洞位置施加固定支撑边界条件。图20、图21分别为试验模型流场网格示意图和结构网格示意图。流场计算湍流模型采用SST模型，得到该状态下的流场，图22 所示为表面热流分布云图。

图20 试验模型流场网格

图21 试验模型结构网格

图22 试验模型表面热流分布云图

采用梯度修正后的紧支C2基函数将气动力/热/结构耦合算例中的热流、温度进行插值以驱动耦合计算的进行，并截取z=0 mm平面内的热流分布曲线得到图23的热流曲线分布对比图，从图中可以发现即使在前缘和压缩拐角等热流高梯度区，流场表面热流与插值后得到的结构表面热流吻合得很好，因此可以将插值后得到的热流进行热响应分析。由于4 s之前流场还未稳定，给模型带来额外的抖动，热流和变形数据均存在较大波动，有效的试验数据是4 s后的数据，因此本文通过耦合计算分析将4 s后的计算结果与试验结果对比，得到沿y方向模型的变形量对比图，如图24所示。结果表明多场耦合仿真结果与试验结果吻合良好，验证了本文采取的数据传递策略在长时间多频次耦合计算状态下的有效性，因此能够将其用于实际工程中的多场耦合分析。

图23 热流曲线分布对比

图24 变形量计算和试验对比

4 结 论

本文提出了一种基于物理量梯度修正的径向基函数数据传递方法，并采用高超声速舵面、高超声速翼身组合体模型和气动力/热/结构多场耦合算例测试了该方法对TPS、MQ以及紧支C2基函数插值效果的影响，得到以下结论：

1)除模型外形特点和网格分布特征外，目标物理量的分布特性同样会对插值结果产生影响。舵面前缘附近目标物理量变化剧烈，会形成较大的梯度，易造成全局方法的数值振荡，而基于物理量梯度修正的TPS、MQ和紧支C2基函数插值方法能很好地消除这一现象。

2)对于存在形态参数的MQ插值方法，在相同的形态参数取值下，采用梯度修正能够很好地降低插值误差，并且采用梯度修正后的MQ插值方法可以在较大的形态参数变化范围内保持很好的精度，因而提高了在实际工程问题中的适用性。

3)基于梯度修正的紧支C2基函数插值效果要高于原始的TPS、MQ等全域基函数，达到甚至超越了基于梯度修正的全域基函数，尤其在翼身组合体等复杂外形中，改进效果更加明显。而与原始的紧支C2基函数插值方法相比，插值时平均选取点数和插值时间下降幅度较大，实现了精度和效率的兼顾。

4)多场耦合算例的计算结果表明采用基于梯度修正的插值方法获得的耦合计算位移变形数据与试验数据吻合度很好，这说明基于梯度修正的径向基函数数据传递方法能够很好地满足气动力/热/结构长时间高频次耦合计算中的数据传递需求。

综上，基于物理量梯度修正的径向基函数多场耦合数据传递方法能够很好地解决实际复杂外形物理量分布各项异性特征所带来的问题，是一种具备应用前景的多场耦合数据传递方法。