火箭结构气动载荷样条插值中的数据预处理

2021-11-17张洋洋于哲峰毛玉明王吉飞

计算机仿真 2021年4期

张洋洋，于哲峰，*，毛玉明，王吉飞

(1. 上海交通大学航空航天学院，上海 200240；2. 上海宇航系统工程研究所，上海 201109)

1 引言

运载火箭结构设计时需要计算其在各种外力作用下各横截面上的轴力、剪力和弯矩等内力。其所受载荷[1][2]主要有飞行时的气动力、发动机推力、竖立在发射台上时的风载荷等。准确计算气动载荷、风载荷[3]的动力响应对箭体结构设计、各部段强度计算有重要的意义。目前计算火箭动力学响应时，一般把火箭简化成一维梁模型进行计算[4][5]。更为精细化的做法是将计算流体力学方法得出的分布气动载荷转换到三维的结构有限元模型上。这就涉及到气动网格载荷向结构网格节点转换的问题，也就是流固载荷传递问题。

国内外学者在近半个世纪的研究中提出了多种插值方法来解决不匹配网格之间的数据传递问题[6]，如反距离插值法[7]、Kriging方法[8][9]、常体积转换法[10]和径向基函数法(radial basis function，RBS)[11][12]等。Frank[13]进行了多种插值方法的比较，发现径向基函数法的综合性能较好。吴宗敏[14]、Buhmann[15]从数学的角度分析总结了不同基函数的径向基函数特性，其中Harder(1971)[16]、Duchon(1979)[17]等人提出的无限板样条(Infinite Plate Spline，IPS)、薄板样条(Thin Plate Spline，TPS)径向基函数法作为航空宇航工程中最受欢迎的一种插值方法，特别适用于三维表面的插值。

为了进行流固载荷传递，结构网格要尽量与流体网格匹配，因此会保留较多的表面细节，这增加了载荷传递的复杂程度。对于火箭这样的大型结构，在直接使用TPS法和IPS法进行数据插值时也有些独特的问题。如线性方程组条件数过大、流体网格点接近于平面时应使用TPS法还是IPS法、在火箭表面某些外形突变处插值误差较大。为了解决上述问题，本文分别提出了对应的预处理插值算法，基于算例分析这些方法对插值精度的影响。

2 面样条插值原理和使用方法

2.1 IPS法原理

IPS法是一种二维插值方法，适用于结构网格点与气动网格点位于相同平面的情况。由无限薄板受分布载荷q引起弯曲小变形的模型推导得出，具体可参见参考文献[16][18]。其中需要求解方程组

(1)

(2)

求解式(2)，可求出a0，a1，a2，F1，F2，…，FN，同时把待插值的坐标一起代入式(1)可得待插值点的值，本文中也称为结构网格点的压力值，也称为插值点的值。

2.2 TPS法原理

TPS法是一种三维插值方法，适用于曲面插值，其推导计算过程与IPS法相同，只是比IPS法多了一维，此时对应的N+4个未知数的线性方程组为

(3)

求出系数a0，a1，a2，a3，F1，F2，…，FN后，同时把待插值点的坐标(x，y，z)代入

(4)

可得待插值点的值。

由方程组(3)可知，要正确解出系数ai和Fj，需要保证线性方程组(3)系数矩阵非奇异。矩阵非奇异出现在两种情况下[18]：第一，如果两个气动网格点有相同的x，y和z坐标(或者离得非常近)，此时，只选取其中一个代入计算，舍去另一个；第二，所有的气动网格点位于同一平面内。火箭模型中大部分含有平面的整流罩会出现第二种情况，可改用IPS法。对于火箭的弧形表面，当气动网点距离较近时，也会近似地在一个平面内，这时也要将数据点投影到一个平面内，采用IPS法进行插值。

2.3 IPS法的具体使用方法

使用IPS法前需要程序自动判断N个气动点与结构点是否都在同一平面。判断共面的思路为先从N个气动点中选出3个点A、B、C确定一个平面，然后判断第四个点D是否在此平面内。计算出AB、AC、AD向量，并单位化，若

(5)

则认为四点共面。以此方法判断其它所有点是否在AB，AC所在的平面上。

然而，实际使用式(5)过程中，由于存储误差或者计算精度误差，或者多个气动网格点近似位于同一平面的情况，等式右端通常不可能等于0。故需要将判断条件改为小于一个参考值

(6)

根据式(6)判断好最近的N个气动网格点全部在同一平面后，采用IPS法求解的下一步需要将原气动点坐标进行转换。由于采用样条插值的一个前提是结构外表面外形需要与气动外形完全一致，而实际上两者一般不一致，在实际处理时默认结构网格点和与之最近的N个气动网格点在同一个曲面或平面。将气动点的坐标原点移动到结构点，将距结构点最远的一个气动点连线为新坐标系x轴，将气动点共面的平面的法向量作为新坐标系的z轴，最后按右手规则确定y轴，建立新坐标系。

N个气动点原点移到结构点后，此时可直接使用坐标转换式(7)[19]，计算新坐标系下气动点的坐标值

(7)

当N个气动网格点严格在同一平面上，经过转换后的坐标值中z坐标都应该为0，如果近似在一个平面上，则z坐标为气动网格点距此平面的距离，接近0。然后将新坐标系下的x,y坐标值代入式(1)即可最终可求出该结构点的压力值。

3 数据处理方法

3.1 结构网格点与气动网格点坐标值量级的预处理

由于流体网格点非常密，相邻流体网格中心点平均距离一般为几毫米，但火箭长度一般为几十米，直径为几米，而结构网格点附近最近的N个气动网格点覆盖的范围很小。因此方程(2)、(3)中，会出现x，y，z很大，r很小的情况，在采用主元消去法求解时，会多次出现特别小的主元，此时矩阵的条件数非常大。Boyd[20]、Diederichs[21]等人认为样条径向基函数方程组条件数与插值点个数N、插值点之间的距离和分布是否规律等因素有关。为了降低方程系的条件数，需要使系数矩阵元素间数量级尽量接近[22]，因此对结构网格点坐标与气动网格点坐标进行预处理。

观察方程(2)、(3)矩阵中有较多固定元素0和1，因此可使矩阵中其它元素的绝对值大小尽量均布于0至1之间，同时综合考虑TPS法和IPS法的基函数f(r)=r2lnr2的图像，如图1所示，应使最近的N个气动点坐标值与其之间距离的大小在0到1.3左右之间。通过以下两个步骤，可实现这个条件：

图1 函数f(r)=r2lnr2的图像

1)在每个结构点压力的计算过程中，均以结构点为原点，即将最近的N个气动点的坐标均减去该结构点的坐标。在使用IPS法进行坐标转换时，已经包含了这种处理，所以这一个处理主要针对TPS法；

2)将气动网格点与结构网格点的坐标值放大适当系数k。k可取为经过第一步处理后最近N个气动点的x，y，z坐标值中的最大值的倒数。然后将最近的N个气动点坐标值均乘以k。此时，N个气动点x，y，z坐标绝对值均小于等于1。结构网格点与气动网格点的坐标的单位对系数矩阵的条件数无影响。求解每一个结构网格点压力时k均取为不同值。

实践证明，通过上述两步预处理后可提高样条插值的精度和鲁棒性。

3.2 TPS法和IPS法的选取(平面度的选取)

当N个流体网格点接近于平面时，接近平面的程度可用参数平面度α来衡量，若式(8)中平面度α取很小，即把很接近共面的气动点判断为非共面，从而采用TPS法，由于气动点接近共面，线性方程组条件数通常会非常大，方程x的系数矩阵为病态矩阵，解的误差会很大。图2中举了某个例子说明了条件数与平面度α的关系，该例中点数为8。若α取较大时，此时程序会把不严格共面的气动点判断为共面，从而采用IPS法，程序实际求得的结果为把不共面的气动点和结构点投影到该平面，最终实际求得的压力值为结构点在在该平面的投影点位置的压力值。若气动点和结构点离该平面越大，则计算误差也会相对越大。综合考虑以上两点因素，需要选取一个合适大小的α值，保证综合误差最小。

图2 条件数与平面度的关系

实践证明，当α取0.01数量级及以上时，不会出现极端情况，包括使用单精度(float型)，采用主元消去法解方程组出现大数除小数而导致无法求解的情况，且方程条件数通常在107以下。

3.3 压力突变区域气动数据点的取舍

在火箭局部表面存在外形突变的区域(如整流罩等处)，故分布在不同曲面的气动点的压力也存在明显的突变。图3(a)中类似凸台左半边压力明显大于右半边，存在一个明显的压力突变截面。此时，若待插值点在此截面附近，如图(b)中的24442号单元，则最近的N个插值点容易在该截面两端均有分布，从而给插值运算带来较大误差。例如待插值点在此截面左边，理论上只应选取截面左边的插值点，若N个插值点中有部分点在截面右端，则实际求得的待插值点压力值会偏小。

图3 压力突变面示意图

为了解决此问题，同时针对整个火箭多种不同外形突变，本文提出了一种数据筛选算法，使最终选取的插值点能与待插值点的处于突变区域的同一侧表面上。该智能判断算法可分两步完成：

1) 利用压力的突变，区分出该截面两端各自包含的插值点，即将N个插值点分为两类。首先根据相邻气动网格点压力差值平均大小，人为设定一个压力差判断阈值，阈值大小可取为离结构单元最近的一个气动插值点A1的压力值的20%，然后分别计算出剩下9个插值点与A1的压力差值，若差值小于压力差判断阈值，则认为该插值点与A1属于同一表面(第一类)，若差值大于阈值，则判断为该插值点属于另外的表面(第二类)。

2) 判断出结构点属于哪一类。分别计算出待插值点距上述两类点中心距离。若待插值点距第一类中心点距离近，则判断为待插值点也属于第一类，此时剔除距离较大的第二类插值点，只按第一类插值点进行插值。

该智能判断算法不需要人为指出表面突变区域，可对全箭压力突变区域自动识别并分类挑选插值点，使用方便，且大大减小了压力突变带来的误差。

4 算例

本文的算例为某型运载火箭，火箭外形及气动载荷压力分布如图4所示，表面气动网格点共614209个，有限元模型中表面结构网格点共40780个，整体网格如图5所示，图6为部分整流罩及变截面处局部细节图。气动模型与结构模型的外表面一致，计算在某自由飞行工况下的气动压力向结构有限元网格的插值转换。

图4 火箭外形及气动载荷压力分布

图5 结构有限元网格示意图

图6 火箭整流罩及变直径处的局部网格示意图

4.1 结构网格点坐标与气动网格点坐标的预处理

为方便说明，令N=10，取第44385号结构单元为例，其原始系数矩阵A未处理情况下TPS法的系数矩阵的条件数为470618。通过3.1节所述的两个预处理步骤后，第44385号结构单元的条件数降为61936。对于本文算例中40780个火箭结构单元，即40780个方程的平均条件数由处理前的9.32e+8降为处理后的3.89e+4，平均下降了约23598倍，即下降了5个数量级。

4.2 TPS法和IPS法的选取

本文算例中取α=0.01，共有6781个结构网格点的气动点被判断为在同一平面，采用了IPS法进行了插值计算。此时，方程组(3)、(6)条件数大于10e+5的个数约占比整体待插值点总数的17%，这使求解方程时的误差较小。

4.3 压力突变区域插值点的取舍

结合本文算例，以图3中第24442号结构单元的插值计算为例，说明数据筛选算法的具体步骤：

1)利用压力的突变，将N个插值点分为两类。第24442号结构单元中心点及其最近的10个气动网格点相对位置如图7所示，左边7个插值点压力值范围在4890 Pa至5195 Pa之间，右边3个插值点压力值在2916 Pa至3397 Pa之间。离结构单元最近的一个气动插值点A1压力值为5196，阈值大小为1039，显然左边7个插值点与A1的压力差值小于阈值，认为该插值点与A1属于同一表面(第一类)。右边3个插值点与A1的压力差值大于阈值，属于另外的表面(第二类)。

2) 通过待插值点距上述两类点中心距离判断出结构点属于哪一类。本例中结构点离第一类中心点较近，故采用第一类的气动网格点进行插值。

图7和图8分别表示使用和未使用该智能判断算法进行的选点及插值结果，表1为数值的对比，可见，筛选前后两者结果相差较大，前者结果更为合理。

表1 使用和未使用智能判断算法的插值结果对比

图7 未使用智能判断算法筛选插值点的插值结果

图8 压力突变区域数据筛选后的数据示意图

表2和表3分别为使用和未使用该智能判断算法进行气动压力载荷向结构节点载荷转换后，两者整体合力、合力距的误差对比。表中数据表明该算法可提高整体插值精度。

表2 没有进行压力突变区域数据筛选的合力、合力距误差

表3 压力突变区域数据筛选后的合力、合力距误差

上述表中合力、合力距误差差别较小，主要因为有两个。第一，压力突变区域总面积占比较小，40780个结构网格点中共有614个结构点采用了数据删选算法，占比不到2%。第二，采用数据删选算法后计算得出的压力值通常与不采用此方法计算出的压力值相差通常在0-200Pa之间，然而压力值通常在5000-6000Pa，误差占比本身较小，大部分在0-4%之间。若压力值也较小(如也在200Pa左右)，此时是否采用数据删选算法的区别则会大大增加。

使用本文所述方法插值后的结构压力云图如图9所示。