第38卷第4期________________________________计算机仿真___________________________________2021年4月文章编号：1006 -9348(2021 )04 -0047-06
火箭结构气动载荷样条插值中的数据预处理
张洋洋、于哲峰^ ，毛玉明2，王吉飞2
(1.上海交通大学航空航天学院，上海200240;
2.上海宇航系统工程研究所，上海201109)
摘要：在运载火箭结构载荷精细化设计中，需要把三维流场压力载荷向结构网格等效转换，采用了薄板样条法和无限板样条法进行具有复杂表面火箭结构的载荷转换。针对应用中因线性方程组条件数巨大而导致解的误差较大、薄板样条法与无限
板样条法插值的选择和火箭表面某些外形突变处插值误差较大的三个问题，提出了对插值点和待插值点的坐标进行预处
理、使用平面度判断方法和智能判断算法三个预处理方法。运用上述方法后线性方程组的条件数大大降低，压力突变处插
值误差大大减少，最终使转换前后合力和合力矩偏差在可接受范围内。以上处理方法的应用有效地防止计算出错，提高了 计算稳定性和载荷转换精度，增强了火箭结构三维气动载荷计算的工程可实施性。
关键词：火箭结构；气动载荷转换；样条；条件数；线性方程组
中图分类号:V448. 15 +3 文献标识码：B文竹的养殖方法
Data Preprocessing of Spline Interpolation in Aerodynamic
Load of Rocket Structure
ZHANG Yang - yang1，YU Zhe - feng丨，*，MAO Yu - ming2，WANG Ji - fei2
(1. School of Aeronautics and Astronautics,Shanghai Jiao Tong University,Shanghai 200240, China ；
2. Shanghai I n s t i t u t e of Aerospace System Engineering, Shanghai 201109,China)
A B S T R A C T：In the rearch and development of launch vehicle structure load,i t i s necessary t o transform the three
-dimensional flow f i e l d pressure load t o the structure network.Therefore,thin - plate spline method and i n f i n i t e plate spline method were adopted t o transform the load of rocket structure with complex surface.The large condition number of linear equations not only leads t o a larger error of solution but also the interpolation choice of thin plate spline method and i n f i n i t e plate spline method and the interpolation error of abrupt change of rocket surface shape.In t h i s regard,this paper ud the f latness judgment method and i n t e l l i g e n t judgment algorithm t o preprocess the coordinates of inteq^olation points and points t o be interpolated.Data preprocessing reduces the condition number of the l i near equations and the interpolation error a t the pressure mutation.The resu l t s show t h a t t h i s method can e f f e c t i v e l y prevent calculation errors,improve the s t a b i l i t y of calculation and the accuracy of load conversion,and enhance the engineering f e a s i b i l i t y of three - dimensional aerodynamic load calculation of rocket structure.
K E Y W O R D S：Rocket structure；Aerodynamic loads transformation；Spline；Condition number；Linear equations
i引言
运载火箭结构设计时需要计算其在各种外力作用下各 横截面上的轴力、剪力和弯矩等内力。其所受载
荷n)u]主要 有飞行时的气动力、发动机推力、竖立在发射台上时的风载 荷等。准确计算气动载荷、风载荷[3]的动力响应对箭体结构
基金项目：民用航天项目（18-g f b-k g-021 -53)
收稿日期:2019-08-05修回日期:2019 -09-18设计、各部段强度计算有重要的意义。目前计算火箭动力学 响应时，一般把火箭简化成一维梁模型进行计算[4U5]。更为 精细化的做法是将计算流体力学方法得出的分布气动载荷 转换到三维的结构有限元模型上。这就涉及到气动网格载 荷向结构网格节点转换的问题，也就是流固载荷传递问题。
国内外学者在近半个世纪的研究中提出了多种插值方 法来解决不匹配网格之间的数据传递问题:6]，如反距离插值 法[7]、Kriging方法™9]、常体积转换法[1°]和径向基函数法 (radial basis function,R B S)[n][12]等。Frank」.* 进行了多种插
47
值方法的比较，发现径向基函数法的综合性能较好。吴宗 分析这些方法对插值精度的影响。敏[14]、1?此11131111;15]从数学的角度分析总结了不同基函数的
径向基函数特性，其中Harder(1971)[16]、Duchon(1979)[17]2面样条插值原理和使用方法
等人提出的无限板样条（I n f i n i t e Plate Spline,IPS)、薄板样条 (Thin Plate Spline,TPS)径向基函数法作为航空宇航工程中 最受欢迎的一种插值方法，特别适用于三维表面的插值。
为了进行流固载荷传递，结构网格要尽量与流体网格匹 配，因此会保留较多的表面细节，这增加了载荷传递的复杂 程度。对于火箭这样的大型结构，在直接使用T P S法和IPS 法进行数据插值时也有些独特的问题。如线性方程组条件 数过大、流体网格点接近于平面时应使用T P S法还是IPS 法、在火箭表面某些外形突变处插值误差较大。为了解决上 述问题，本文分别提出了对应的预处理插值算法，基于算例2.1 I P S法原理
IPS法是一种二维插值方法，适用于结构网格点与气动 网格点位于相同平面的情况。由无限薄板受分布载荷q引起弯曲小变形的模型推导得出，具体可参见参考文献[16] [18]。其中需要求解方程组
N
w(x,y)=a0 + a xx+a2y+X F/M r,(l)
i=1
式中W为插值的物理量，r f = (* +(y-yi)2,yv为已知插值点的个数。此时，方程（1)的/V+ 3个未知量（a(),ai，a2, 6，&,•••，,）可由如下的线性方程组求出
『丨_
■1尤17i0r n l n r n• • r i v l n r Lv a〇
1欠2r2r l\l n r210• • d v l n r^a i
1无3y3r31l n r31r32l n r32• • r3,\^n r3N o2 ^41丨474☆ l n r^r42l n r42• r4A l n r4A r
1X N y N r m l n r vl r A2^n r.\2• • 0乙-3 000011• • 1f,-2
0000
尤i^2X N
0--000
r i J l* * y N」-F n.
(2)
求解式(2)，可求出a。',，％，&,F2，…，乙，同时把待插 T P S法是一种三维插值方法，适用于曲面
插值，其推导值的坐标一起代入式（1)可得待插值点的值，本文中也称为 计算过程与IPS法相同，只是比IPS法多了一维，此时对应的结构网格点的压力值，也称为插值点的值。N +4个未知数的线性方程组为
2.2 TPS法原理
w,-"1x l r i0r12l n r12 •r lN\n r lN
a〇1X2j i22r21l n r210•• r2N'n r2N a l 1X3r323r31l n r3.r32 l n r32 •• r3N'n r3N a2 1X4Z4r】丨
l n r^r42l n r42 *• r4A,l n r4A.a3 1X5r5Z5r51l n r51r52l n r52 *• r5N^n r5N
1x N y N Z5r2m l n r^,,r^l n r^• 0
000001 1 •• 1^N-3
00000
尤i
x2•• x N心-2
00000
r i r2 •• y N F n-\
0-
-0000h Z2•Z N --P's-
(3)
求出系数a。,a W p M F v，心后，同时把待插值 点的坐标U,y,z)代人
N
u>{x,y,z)= a0 + a^x+a2y+a}z+ ^F^lnr] (4)
i = i
可得待插值点的值。
由方程组（3)可知，要正确解出系数a，和厂，需要保证 线性方程组（3)系数矩阵非奇异。矩阵非奇
异出现在两种情况下[1S]：第一，如果两个气动网格点有相同的和z坐标 (或者离得非常近），此时，只选取其中一个代人计算，舍去另 一个；第二，所有的气动网格点位于同一平面内。火箭模型 中大部分含有平面的整流罩会出现第二种情况，可改用/PS 法。对于火箭的弧形表面，当气动网点距离较近时，也会近 似地在一个平面内，这时也要将数据点投影到一个平面内，采用IPS法进行插值。
48
2.3 IPS法的具体使用方法
使用IPS法前需要程序自动判断/V个气动点与结构点 是否都在同一平面。判断共面的思路为先从W个气动点中 选出3个点4、fl、C确定一个平面，然后判断第四个点II是 否在此平面内。计算出向量，并单位化，若
x A^-A S=0(5)则认为四点共面。以此方法判断其它所有点是否在/1BMC 所在的平面上。
然而，实际使用式(5)过程中，由于存储误差或者计算精 度误差，或者多个气动网格点近似位于同一平面的情况，等 式右端通常不可能等于0。故需要将判断条件改为小于一个 参考值
AB•AD<=a(6)根据式(6)判断好最近的/V个气动网格点全部在同一平面后，采用/PS法求解的下一步需要将原气动点坐标进行 转换。由于采用样条插值的一个前提是结构外表面外形需 要与气动外形完全
戚风小蛋糕
一致，而实际上两者一般不一致，在实际 处理时默认结构网格点和与之最近的/V个气动网格点在同一个曲面或平面。将气动点的坐标原点移动到结构点，将距 结构点最远的一个气动点连线为新坐标系x轴，将气动点共 面的平面的法向量作为新坐标系的z轴，最后按右手规则确 定y轴，建立新坐标系。
/V个气动点原点移到结构点后，此时可直接使用坐标转 换式(7)[19]，计算新坐标系下气动点的坐标值
y\=
■
式中:A，.K’*,
(7) .cos(*>,x)cos(x' ,y)cos(x,,z) xh
cos(y,,x)cos(y,,y) cos(y,,i) yk
■cosCz" ,i) cos(z' ,y)cos(z' ,z)L Z),.
匕，匕为新坐标系下的气动点坐标值，万，万为 原坐标系下的坐标值表示原坐标系•轴与新坐标系x轴夹角的cos函数值。
当~个气动网格点严格在同一平面上，经过转换后的坐 标值中Z坐标都应该为0，如果近似在一个平面上，则Z坐标 为气动网格点距此平面的距离，接近0。然后将新坐标系下 的坐标值代人式（1)即可最终可求出该结构点的压力值。风水口诀
3数据处理方法
3.1结构网格点与气动网格点坐标值量级的预处理
由于流体网格点非常密，相邻流体网格中心点平均距离 一般为几毫米，但火箭长度一般为几十米，直径为几米，而结 构网格点附近最近的W个气动网格点覆盖的范围很小。因此方程(2)、（3)中，会出现*,7,2很大，『很小的情况，在采用 主元消去法求解时，会多次出现特别小的主元，此时矩阵的 条件数非常大。Boyc^2°、Diederichs'〜等人认为样条径向基 函数方程组条件数与插值点个数N、插值点之间的距离和分 布是否规律等因素有关。为了降低方程系的条件数，需要使系数矩阵元素间数量级尽量接近[22]，因此对结构网格点坐
标与气动网格点坐标进行预处理^
观察方程（2)、（3)矩阵中有较多固定元素0和1，因此
可使矩阵中其它元素的绝对值大小尽量均布于〇至1之间，
同时综合考虑T P S法和IPS法的基函数f(r) =r2hr2的图
像，如图1所示，应使最近的N个气动点坐标值与其之间距
离的大小在0到1.3左右之间。通过以下两个步骤，可实现
这个条件：
1) 在每个结构点压力的计算过程中，均以结构点为原 点，即将最近的N个气动点的坐标均减去该结构点的坐标。
在使用IPS法进行坐标转换时，已经包含了这种处理，所以61年属什么
这一个处理主要针对T P S法；
2) 将气动网格点与结构网格点的坐标值放大适当系数
k。 k可取为经过第一步处理后最近N个气动点的x,y，z坐
标值中的最大值的倒数。然后将最近的N个气动点坐标值
均乘以k。此时,N个气动点x，y，z坐标绝对值均小于等于
l。 结构网格点与气动网格点的坐标的单位对系数矩阵的条
件数无影响。求解每一个结构网格点压力时k均取为不
同值。
实践证明，通过上述两步预处理后可提高样条插值的精
度和鲁棒性。
图1函数f(r> s H/n r2的图像
3.2 TPS法和IPS法的选取（平面度的选取）
当N个流体网格点接近于平面时，接近平面的程度可用
参数平面度a来衡量，若式（8)中平面度a取很小，即把很
接近共面的气动点判断为非共面，从而采用T P S法，由于气
动点接近共面，线性方程组条件数通常会非常大，方程x的
系数矩阵为病态矩阵，解的误差会很大。图2中举了某个例
子说明了条件数与平面度a的关系，该例中点数为8。若a
取较大时，此时程序会把不严格共面的气动点判断为共面，
从而采用IPS法，程序实际求得的结果为把不共面的气动点
和结构点投影到该平面，最终实际求得的压力值为结构点在
在该平面的投影点位置的压力值。若气动点和结构点离该
平面越大，则计算误差也会相对越大。综合考虑以上两点因
49苹果手机电脑
素，需要选取一个合适大小的a 值，保证综合误差最小。
图2
条件数与平面度的关系
实践证明，当£»取〇.〇1数量级及以上时，不会出现极端 情况，包括使用单精度（float型），采用主元消去法解方程组 出现大数除小数而导致无法求解的情况，且方程条件数通常 在1〇7以下。
3.3压力突变区域气动数据点的取舍
在火箭局部表面存在外形突变的区域（如整流罩等处）， 故分布在不同曲面的气动点的压力也存在明显的突变。图3 U )中类似凸台左半边压力明显大于右半边，存在一个明显 的压力突变截面。此时，若待插值点在此截面附近，如图（b ) 中的24442号单元，则最近的N 个插值点容易在该截面两端 均有分布，从而给插值运算带来较大误差。例如待插值点在 此截面左边，理论上只应选取截面左边的插值点，若N 个插 值点中有部分点在截面右端，则实际求得的待插值点压力值 会偏小。
为了解决此问题，同时针对整个火箭多种不同外形突 变，本文提出了一种数据筛选算法，使最终选取的插值点能 与待插值点的处于突变区域的同一侧表面上。该智能判断 算法可分两步完成：
1)
师训平台利用压力的突变，区分出该截面两端各自包含的插
值点，即将N 个插值点分为两类。首先根据相邻气动网格点 压力差值平均大小，人为设定一个压力差判断阈值，阈值大 小可取为离结构单元最近的一个气动插值点A ，的压力值的 20% ，然后分别计算出剩下9个插值点与A ，的压力差值，若 差值小于压力差判断阈值，则认为该插值点与A ，属于同一 表面(第一类），若差值大于阈值，则判断为该插值点属于另 外的表面（第二类）。
2)
判断出结构点属于哪一类。分别计算出待插值点距
上述两类点中心距离。若待插值点距第一类中心点距离近， 则判断为待插值点也属于第一类，此时剔除距离较大的第二 类插值点，只按第一类插值点进行插值。
该智能判断算法不需要人为指出表面突变区域，可对全
人中旁边有痣
L
(b )结构网格局部图
图3压力突变面示意图
箭压力突变区域自动识别并分类挑选插值点，使用方便，且 大大减小了压力突变带来的误差。
4算例
本文的算例为某型运载火箭，火箭外形及气动载荷压力
分布如图4所示，表面气动网格点共614209个，有限元模型 中表面结构网格点共40780个，整体网格如图5所示，图6为 部分整流罩及变截面处局部细节图。气动模型与结构模型 的外表面一致，计算在某自由飞行工况下的气动压力向结构 有限元网格的插值转换。
图4
火箭外形及气动载荷压力分布
4.1结构网格点坐标与气动网格点坐标的预处理
为方便说明，令N = 10，取第44385号结构单元为例，其
50
含有反义词图5结构有限元网格示意图号结构单元中心点及其最近的10个气动网格点相对位置如 图7所示，左边7个插值点压力值范围在4890 P a至5195 Pa 之间，右边3个插值点压力值在2916 P a至3397 P a之间。离结构单元最近的一个气动插值点A，压力值为51%，阈值大 小为1039，显然左边7个插值点与A，的压力差值小于阈值，认为该插值点与A，属于同一表面（第一类）。右边3个插值 点与A，的压力差值大于阈值，属于另外的表面（第二类）。
2)通过待插值点距上述两类点中心距离判断出结构点 属于哪一类。本例中结构点离第一类中心点较近，故采用第 一类的气动网格点进行插值。
图7和图8分别表示使用和未使用该智能判断算法进 行的选点及插值结果，表1为数值的对比，可见，筛选前后两 者结果相差较大，前者结果更为合理。
表1使用和未使用智能判断算法的插值结果对比是否使用智能
判断算法
插值结果/P a插值点个数使用5092 7
未使用4411 10
原始系数矩阵A未处理情况下T P S法的系数矩阵的条件数 为470618。通过3.1节所述的两个预处理步骤后，第44385 号结构单元的条件数降为61936。对于本文算例中40780个 火箭结构单元，即40780个方程的平均条件数由处理前的9.32e +8降为处理后的3.89e +4，平均下降了约23598倍，即下降了 5个数量级。
4.2 TPS法和IPS法的选取
本文算例中取a =0. 01，共有6781个结构网格点的气动 点被判断为在同一平面，采用了 IPS法进行了插值计算。此 时，方程组(3)、(6)条件数大于10e + 5的个数约占比整体待 插值点总数的17% ，这使求解方程时的误差较小。
4.3压力突变区域插值点的取舍
结合本文算例，以图3中第24442号结构单元的插值计 算为例，说明数据筛选算法的具体步骤：
I)利用压力的突变，将N个插值点分为两类。第24442
图7未使用智能判断算法筛选插值点的插值结果
图8压力突变区域数据筛选后的数据示意图
表2和表3分别为使用和未使用该智能判断算法进行 气动压力载荷向结构节点载荷转换后，两者整体合力、合力 距的误差对比。表中数据表明该算法可提高整体插值精度。
5
1
