亚声速条件下五孔探针非线性插值计算方法研究*

2020-09-28史远鹏王佳璐陆华伟孔晓治

风机技术 2020年4期

史远鹏王佳璐陆华伟孔晓治王宇

（1.大连海事大学船舶与海洋工程学院；2.大连海事大学信息科学技术学院）

0 引言

风洞测试离不开流场测控，目前测量流场的方法主要有多孔气动探针[1]、激光多普勒测速仪[2]、热线风速仪等[3]。多普勒测速仪可以测得流场的速度大小以及速度方向，但是设备较昂贵，应用场景范围小，在工程应用上受到限制[4]，主要应用于实验室中；热线风速仪的探头容易被空气中的尘埃颗粒损坏，且测量的并非是一个点上的流场数据，而是一个小平面内的平均值。

多孔气动探针测量法作为测量流场最经典的方法，不仅可以方便测得流场数据，而且具有测量精度高、成本较低、使用方便和测量原理简单等优点，所以多孔气动探针在实验室和工程上受到了广泛的应用。气动压力探针的相关理论及应用已经十分成熟，早在1971年Bryer和Pankhurst已经对气动压力探针进行了研究[5]，其中锥形五孔探针在Ma=0.1～1.2范围内有较好的测量效果[6]，本文研究的对象就是应用较为广泛的五孔探针。根据Treaster和Yocum[7]的研究，将五孔探针测量三维流场分为三种方法：对向测量法、半对向测量法、非对向测量法。

1）对向测量法：比较直观，但是需要较多时间寻找1,3孔，4,5孔的压力平衡，实际测量较为复杂；2）半对向测量法：操作简单，只需要寻找一对压力孔的压力平衡，所需数据处理量较少，适合均匀流场测量；3）非对向测量法通过采集五个孔的压力值，对照该探针的校准数据，可以插值计算出测量点的俯仰角、偏航角、总压、静压，这种方法的使用最为普遍，数据处理量较大，需要较为准确的插值方法，实际测试系统借助计算机软件编程处理数据，有较快的测试速度。

在探针的制造过程中，每支探针的头部五孔并不严格对称，这就导致了每支探针的气动特性都不相同，工程中往往要求五孔探针在多个马赫数下使用，需要对各个马赫数分别进行标定，且标定结果只能在一定的马赫数变化范围内使用。

国内陈浮[8]等人提出了五孔探针的线性插值方法，该方法对凸四边形有良好的插值效果，并实现了对原始校准数据奇异点的修正。在实际测量过程中，当来流速度与标定速度不一致时，传统的线性插值法选择就近标定马赫数所对应的特性曲线作为校准数据，利用测得的五孔压力值计算出Kα，Kβ，寻找出点（Kα，Kβ）在校准曲线中的位置后计算点（Kα，Kβ）距所在四边形各边的距离，用线性插值的方法得到探针感受到的气流角度α，β，再根据α，β数值在总压、静压校准系数曲线中得到测点处的CPt，CPs，最后由相应公式求出总压、静压及空间3个方向的速度分量，并以此作为最终所测量点的流场参数，此过程流程图如图1。

图1 线性插值计算流程Fig.1 Linear interpolation calculation flow

此外，岳国强[9]等人对原有的线性插值法进行了改进。尤其是对于特性较差的探针，该方法能够在不改变探针特性曲线的基础上保证插值精度，提高了线性插值法的实用性。在实际的流场测量中，线性插值方法对于变化的流场没有很好的适应性，当来流马赫数超出标定速度一定范围时，计算结果存在较大误差。为进一步提高测量精度，本文对线性插值法进行了理论研究，并在此基础上提出了可以减小误差的三维非线性插值法。

1 五孔探针结构与标定方法

目前实际应用的五孔探针形式众多，例如有：圆锥形、球形、棱形和蛇形等[10]。其中圆锥形五孔探针在测量高速流场时效果较好，圆锥形五孔探针的头部角度多为60°，常见的设计为“L”型，即直角弯曲圆锥形五孔探针，探针头部与后面的探针杆夹角90°。五孔探针的直径尽量小能减少探针对流场的扰动，但制造难度也相应增加、孔径太小空气中的灰尘等杂质容易堵塞测量孔。国内的探针多将头部设计成外径2～5mm,孔径1mm[11]，余世策等人利用三维雕刻机实现了小尺寸多孔探针的快速加工[12]。

本文以L型圆锥形头束状五孔探针为研究对象，该型五孔探针孔序及各角度定义如图2和图3,定义P1～P5为1～5孔的压力值，Pt为总压，Ps为静压。

图2 五孔探针的孔序定义Fig.2 The definition of pore sequence of five-hole probe

图3 五孔探针的角度定义Fig.3 The definition of angles of five-hole probe

非对向测量法，标定时已知来流角度、总压、静压，利用这五个压力值，通过公式(1)～(5)[7]获得对应来流角度下的Kα，Kβ，CPt，CPs。其中Pˉ为1，3，4，5孔平均压力；Kα为α方向校准系数；Kβ为β方向校准系数；CPt为总压校准系数；CPs为静压校准系数[13]。

标定的任务就是找出α，β，Kα，Kβ，CPt，CPs之间的关系，绘出探针特性曲线，使用时根据五孔探针所采集到的各孔压力值计算各个标定系数，在特性曲线中进行插值得出流场参数。

在实际的流场测量中，探针前的来流速度往往与标定马赫数存在一定差值，偏离标定马赫数时，传统方法是选择就近标定马赫数对应的特性曲线近似代替，这会导致插值结果偏离真实值，结果的误差随来流速度偏离标定马赫数的变大而变大，因此在测量来流较复杂的流场时，就需要增加大量的探针标定马赫数，费用高，耗时长，为解决这个问题本文提出来三维非线性插值法。

2 三维非线性插值法

2.1 校准系数拟合阶数的选择原则

拟合函数表达式的阶数选择对拟合效果影响较大，因此拟合阶数的选择是进行三维非线性插值之前需要探究的一个重要问题。为反映不同标定马赫数下，特性曲线间的变化规律，找到最佳的拟合阶数，使得插值得到的特性曲线更准确，本文利用最小二乘法对特性曲线进行了拟合，得出相应的函数表达式[14]，总结出了进行三维非线性拟合时的具体方法。

函数原型如上式，其中a0，a1，a2，a3，…，ak为一元的多项式系数；k为拟合阶数；n为标定马赫数；Y为Kα，Kβ，CPt，CPs校准系数。

以一组五孔探针标定的Kα校准系数实验值为例，马赫数范围从0.2～0.8马赫，共7个值，由最小二乘法的性质可知，七个实验值最高可拟合出六阶多项式。下面对数据分别进行一阶、二阶、三阶、四阶、五阶、六阶拟合，将插值结果与标定实验值对比，讨论各阶的拟合效果。

根据多项式本质，多项式是无边界的振荡函数，因此，它们不太适合外插有界数据或单调数据（递增或递减）的数据，从图4中同样可以得到印证，当超出实验值的范围，拟合函数快速偏离实际值，高阶六阶、五阶尤其明显；高阶多项式在数据点之间可能是振荡的，导致对数据的拟合较差，即并不是阶次越高，拟合效果越好，在图4中可以看到，0.2～0.3马赫段的高阶六阶、五阶函数产生了较大幅度的振荡。

图4 Kα实验值的不同阶拟合曲线Fig.4 Different order fitting curves ofKαexperimental values

用上述7个实验值中的6个进行1～5阶基于最小二乘法的多项式拟合，插值得到另外一个的计算值并与实验值比较，讨论阶数的选择造成误差的不同。例如用0.3～0.8马赫的6个实验值分别进行1～5阶拟合，插值得到0.2马赫下的Kα校准系数计算值。特别说明的是，插值0.2马赫、0.8马赫属于上面所述的外插情况，计算结果偏离了正常范围。

观察图5数据，插值0.2马赫、0.8马赫相当于外插的情况，误差远大于其他马赫数下的值，插值时应避免这种情况，当拟合阶数k为2阶、3阶或4阶时，各个数据点误差值的绝对值变化较小，折线图波动较小。同时从图6可以看出，拟合阶数k为2阶、3或4阶时的误差均值明显小于1阶和5阶时的平均误差值，即校准系数拟合阶数k为2阶、3阶或4阶时插值结果与各马赫下的实际标定实验值之间的误差最小。

图5 误差值折线图Fig.5 Error value line chart

图6 平均误差值柱状图Fig.6 Bar chart of average error value

通过对多只不同特性的探针的Kα,Kβ,CPt,CPs进行上述拟合插值，均能得到类似的结果。从而，在选择拟合阶数时，提出以下四点参考原则：

1）不选高阶，高阶多项式虽然更贴合数据点，但在数据点之间存在较大振荡，数据的拟合效果并不好，插值结果不准确；

2）不选低阶，阶数太低，数据点离曲线较远，插值结果与实际值误差也较大；

3）在亚声速范围内进行标定，数据点足够时，拟合阶数尽量选取2阶、3阶或4阶；

4）在数据点区间外，拟合函数很快偏离实际值，该方法不适合进行外插。

2.2 三维非线性插值法原理

三维非线性插值法的具体方法——先对一个范围内若干个等间距马赫数进行标定，通过对所标定各个马赫数特性曲线相同α，β位置的特性值（Kα,Kβ,CPt,CPs）与马赫数进行拟合，得到马赫数与Kα,Kβ,CPt,CPs的关系，即若干一元多次函数，通过插值便可得到该马赫区间内所有马赫数下的特性曲线。图7是0.2马赫、0.5马赫、0.8马赫的特性曲线（黑色）进行三维非线性拟合，插值得到0.6马赫特性曲线（红色）的示意图，图中作出了四条（Kα，Kβ-Ma）拟合曲线。

图7 三维非线性插值示意图Fig.7 3-D nonlinear interpolation graph

2.3 寻找最佳标定马赫数个数

本文对一五孔探针在0.2，0.3，0.4，0.5，0.6，0.7，0.8马赫下进行了标定，实验过程中通过改变风洞鼓风机频率大小调节标定风洞出口马赫数，自动测试系统对步进电机发送指令，实现对坐标架的精确控制。在各个马赫数下，探针扫过的角度范围为±30°，角度间隔为5°，采集点为169个，即采取的角度分别为α=±30°，±25°，±20°，±15°，±10°，±5°，0°，β=±30°，±25°，±20°，±15°，±10°，±5°，0°。先将β角转至-30°，α角由-30°至30°，以5°为间隔逐一扫过，之后β角增大5°，α角重复上述步骤，直至13×13个点的数据全部采集完毕。

根据所采集到的数据，计算若干来流角度下的四个校准系数。理论上标定马赫数间隔越小，插值结果越准确，但由于时间、设备、人员等成本的限制，本次标定在0.2～0.8上等间距选择了七个马赫数。如图8在0.2～0.8马赫范围内，逐渐增加参与拟合的特性曲线个数，同时作出各种情况下插值得到的0.6马赫特性曲线，以0.6马赫下实际标定的特性曲线为对照，计算并绘制出平均绝对误差云图，寻找该区间最佳的标定个数。在选择标定马赫数时尽可能的均匀，同时依照前面所述的拟合阶数选择原则，使各个对比实验均达到了较好的拟合效果。

图8 三维非线性插值结果和误差云图Fig.8 3-D Nonlinear Interpolation Results and Error nephogram

由以上结果，得知两个马赫数一阶拟合得到特性曲线误差较大，小角度下误差在1.5%左右，在大角度下的误差达到4%以上，拟合精度较差；三个马赫数二阶拟合得到的特性曲线仅在大角度时存在较大误差，最大误差为3.5%，在80%角度范围内误差在2%以下，小角度时误差在1%以下，与两个马赫数拟合结果相比，拟合精度大幅增加；四个马赫数三阶拟合得到的特性曲线虽然误差整体误差均在2.5%以下，但有50%的角度区域误差大于等于2%，有80%的区域误差在1%以上，较三个马赫数拟合精度相比有所下降；五个马赫数三阶拟合得到的特性曲线整体拟合精度较高，最大拟合误差不到2%，有50%的角度区域误差在1%以下，与三个马赫数拟合相比，拟合精度有所提高；六个马赫数四阶拟合，最大误差小于2%，误差小于1%的区域可达70%，小角度区域误差均小于0.5%，拟合精度是所有平行实验中最高的。

表1 平行实验的误差比较Tab.1 Deviation comparison of parallel experiments

通过上述对比实验，我们讨论了不同数量的特性曲线对于非线性插值精度的影响。对大多数工程我们将特性曲线的平均相对误差控制在1%以内，即三个马赫数拟合、五个马赫数拟合、六个马赫数拟合均可满足工程要求。通常对五孔探针在一个马赫数下进行标定需要20～40分钟，权衡考虑节约时间与保证精度时，在亚声速流场测量区间，对三个马赫数下的特性曲线进行三维非线性拟合便可完成整个区间的标定，即只需标定测量来流范围的边界与该范围内中间或者较中间的马赫数所对应的特性曲线，通过拟合插值便可得到此范围内所有马赫数的特性曲线。

2.4 三维非线性插值法与传统线性插值法校准风洞下的误差分析

由于上述叶栅出口截面测试实验，叶栅后流场参数不能准确得知，无法计算三维非线性插值法与线性插值法的具体误差，作者在大连海事大学凌海校区的跨声速校准风洞进行了实验。测试时风洞出口马赫数为0.6，将该校准风洞已知的流场参数作为参考值，三维非线性插值法使用0.2马赫，0.5马赫，0.8马赫校准数据拟合得到的Kα，Kβ，CPt，CPs三维校准曲线为校准数据；线性插值法以0.5马赫校准曲线为校准文件。两种方法对各参数计算结果与参考值比较，得到三维非线性插值法具体的误差如表2所示，表3所示为线性插值法的各参数误差。

表2 三维非线性插值法的误差Tab.2 Error of 3D Nonlinear Interpolation

表3 线性插值法的误差Tab.3 Error of Linear Interpolation

上述误差云图，颜色从蓝到红表示误差由小到大，三维非线性插值法总压误差有80%左右区域在0.2%以下，而线性插值法40%左右的区域总压误差在0.2%以下；同样三维非线性插值法其它各参数误差计算结果均显著优于线性插值法。

另外，实验时由于静压值较小，压力传感器读值时的跳动相对较大，所以静压误差整体相对较大；在测量0.6Ma真值与验证实验过程中，两次测量时五孔探针位置不能保证绝对相同，也可能引入α角与β角的绝对误差。

2.5 实际流场测试实验对比

本文采用五孔探针对某一矩形叶栅出口截面进行了实际测量，风洞出口流场主流区马赫数为0.7，分别利用三维非线性插值法和传统处理方法计算流场参数。对比组采用传统处理方法，即根据0.7马赫的标定数据线性插值法计算，三维非线性插值法使用0.2马赫，0.5马赫，0.8马赫校准数据拟合得到的Kα，Kβ，CPt，CPs三维校准曲线为校准数据。

实验时，根据五孔探针各个孔的压力值由公式（1）、（2）、（3）计算出Kα,Kβ，以马赫（0）=0.7（可任意）的特性曲线作为初次计算的校准文件，利用线性插值法由公式（4）、（5）可以得到来流各方向角、CPt与CPs，根据CPt，CPs得到的总压、静压计算出马赫（1），判断马赫（1）与马赫（0）的差值，若差值大于万分之一，则以马赫（1）所对应的特性曲线再一次进行线性插值计算，经过上述相同步骤即可计算出马赫（2），再比较马赫（2）与马赫（1）差值，以此类推，直到马赫（n）与马赫（n-1）差值小于万分之一，输出最终计算出的来流各方向角、总压、静压以及马赫数等流场参数，即通过不断振荡来逼近真实的来流马赫数，用计算机实现自动采集和校准[15]，上述过程能在极短时间内完成，计算流程如图10所示。

图9 误差云图Fig.9 Error Nephogram

图10 三维非线性插值方法的计算流程Fig.10 Calculation flow of three-dimensional nonlinear interpolation method

上图中，a1,a2为出口截面马赫数云图，其中a1为三维非线性插值法的处理结果，a2为传统方法的处理结果。b1,b2为β角云图，其中b1为三维非线性插值法的处理结果，b2为传统方法的处理结果。c1,c2为总压云图，其中c1为三维非线性插值法的处理结果，c2为传统方法的处理结果。

通过图11的对比，从整体定性上来看两种方法处理的结果规律几乎完全一致，尤其是叶栅流道区域两种方法处理的结果几乎相同。但在高损失区域由于三维非线性插值法对变马赫数流场的较强适应性得到的结果较传统处理方法更为全面。在高损失区域定量分析可得两种方法处理的马赫数、总压最大相差分别是0.369%和39.02%。因此在实际流场测量过程中，尤其是在测量不稳定来流时，三维非线性插值法能适应流场变化，减小来流速度与标定速度不一致所造成的测量误差，得到更综合、细腻的测量结果。