CN107194045A - 一种用于空中加油的受油机前扰波建模方法 - Google Patents

Abstract

本发明公开了一种用于空中加油的受油机前扰波建模方法，属于空中加油技术领域。本发明包括：通过计算流体力学CFD仿真，得到受油机机头附近速度流场数据；以多项式的函数形式用最小二乘法拟合出速度分别沿空间三个轴向上的分布规律，研究一维的速度与位置的函数关系；根据一维函数关系推测速度在机头横截面上的二维表达式，用多项式进行拟合；用平面插值法求解出速度与三维位置坐标的解析表达式，并根据误差曲线来选择插值密度；最后根据标准误差选择合适的插值方向。本发明方法可以有效弥补无人机前扰波建模现有方法过于复杂，模型不精确等问题，计算量小，结果较为精确，对于其他机型的无人机，均可采用本方法进行前扰波建模。

Description

一种用于空中加油的受油机前扰波建模方法

技术领域

本发明涉及一种空中加油自动对接过程中，受油机前扰波的建模方法，属于空中加油技术领域。

技术背景

空中加油对于增加无人机航时和作战半径具有重要作用，近年来得到了越来越多国家和科研机构的关注。随着无人驾驶飞行器的广泛应用，特别是无人作战飞机的兴起，针对于自动空中加油技术的研究成为热点。空中加油目前主要分为硬管式和软管式，软管式空中加油有五个阶段，即会合，编队，对接，加油及分离。其中对接段是最为关键，也是最困难的部分，需要加油机和受油机以很高的精度保持一定的相对位置，并且具有一定的抗扰能力。在对接过程中，扰动主要为大气紊流，加油机尾涡以及受油机前扰波。对于大气紊流和加油机尾涡的建模与分析，国内外已有了相应的研究，发展较为成熟。但是对于前扰波效应，即受油机头部对加油锥套的气动影响，还缺乏足够的重视和充分的研究，飞翼型布局无人机前扰波与常规飞机有较大差异，更是处于空白状态。

在软管式空中加油对接过程中，当受油机接近加油锥套时，受油机前扰波作用所产生的气动影响，会使锥套产生不同程度的飘摆和偏离，严重制约着软式空中加油对接成功率。国外对受油机前扰波的研究主要基于飞行试验、数学模型和计算流体力学等方法。NASA德莱顿飞行研究中心通过F/A-18双机自动空中加油飞行测试，研究了飞行条件、受油机机动等因素对软管锥套组合体飘摆特性影响，并分析了前扰波影响范围和受油机动态对锥套运动影响规律。布里斯托尔大学利用兰金半体模型模拟前扰波影响，并设计了人工补偿器补偿受油机前扰波造成的锥套偏离。目前国内对受油机前扰波的研究较少，戴训华和魏子博分别运用流函数法和拟合法对前扰波进行建模，得到了F16的前扰波模型(参考文献1：Dai Xunhua,Wei Zibo,Quan Quan.Modeling and simulation of bow wave effect inprobe and drogue aerial refueling，2016；参考文献2：Wei Z B,Dai X,Quan Q,etal.Drogue dynamic model under bow wave in probe-and-drogue refueling，2016)。

总的来说，在空中加油领域中，对于受油机前扰波建模的研究相对较少。现有的研究成果存在着建模方法过于复杂，模型不精确等问题。迫切需要新的前扰波研究方法，给出相对准确，简单的前扰波模型。

发明内容

本发明的目的是为了解决在无人机空中加油自动对接过程中，现有的方法不能准确简洁地建立出受油机前扰波模型的问题，提出了一种基于最小二乘的多项式拟合建模方法，该拟合方法根据计算流体动力学(CFD)仿真数据在XYZ三个轴向上进行多项式拟合，基于拟合的结果进而在二维平面上进行拟合，再运用平面插值法得到整个三维区域内气流速度与坐标关系的解析表达式，完成受油机前扰波建模。

本发明提供的用于空中加油的受油机前扰波建模方法，包括步骤：

步骤一：进行受油机前扰波CFD仿真，仿真获取三维空间内大量坐标点(x,y,z)所对应的三维速度值(v_x,v_y,v_z)。

步骤二：用多项式拟合一维的速度与位置的函数关系；

对任一坐标轴上的速度v，用多项式分别拟合该速度v随每一轴上位置变化的关系，具体是：固定z值，选取不同y值，用多项式拟合速度随x值变化的关系；固定y值，选取不同x值，用多项式拟合速度随z值变化的关系；固定x值，选取不同z值，用多项式拟合速度随y值变化的关系；

多项式拟合表达式为：其中c_k为待定系数，n为最高次项，t＝x,y,z。

步骤三：在一维速度与位置的函数关系基础上，进行二维的平面拟合。

对任一坐标轴上的速度v，在tr平面的多项式拟合表达式为：

其中，r＝x,y,z，t和r为不同坐标轴；p,q分别为步骤二中确定的在t轴方向和r轴方向拟合多项式的阶数；c_pq为待定系数。

步骤四：基于步骤三得到的二维平面拟合表达式，选用插值法来得到三维空间中的解析表达式；三维空间的三个轴向对应三个插值方向，分别进行计算；

对任一坐标轴上的速度v，设步骤三计算得到t轴方向上的两个相邻的平行二维平面对应的v取值分别为v_d和v_u，则利用线性插值计算速度v与所述两个相邻的平行二维平面之间的点的关系，表达式如下：

v＝v_d+(v_u-v_d)(t-t_d)/d

式中，t是待计算点在t轴方向上的坐标值，t_d为其中一个二维平面在t轴方向上的坐标值，d为插值密度。

步骤五：选取三个插值方向中平面拟合的标准误差最小的方向作为最合适的插值方向，对速度和坐标进行插值拟合。

本发明的用于空中加油的受油机前扰波建模方法，优点和积极效果在于：

(1)本方法可以直接给出受油机机头附近三维流场速度与位置关系的解析表达式，即给出任意一点的坐标值，可通过解析式得出其速度值。

(2)无论什么样的函数形式，都可以展开成多项式，并且多项式的阶数越高，拟合效果越好。所以采用的多项式拟合精度高，通用性好，可以适用于任何机头形状前扰波的研究。

(3)相比其他函数形式，计算机计算多项式的速度较快，本发明方法实时性较好。

(4)建模简单，占用的存储空间小。

附图说明

图1是受油机机头模型图；

图2是PointWise软件生成的网格划分图；

图3是Fluent软件计算所得速度和静压云图；

图4是z＝-1.2平面速度在X轴上分量v_x的等值线；

图5是z＝-1.2平面X轴向上4阶多项式拟合曲线；

图6是y＝0.8平面速度在X轴上分量v_x等值线；

图7是y＝0.8平面z轴向上5阶多项式拟合曲线；

图8是x＝1.4平面速度在X轴上分量v_x等值线；

图9是x＝1.4平面y轴向上3阶多项式拟合曲线；

图10是三种平面拟合方式示意图；

图11是y＝0.8平面多项式拟合曲面；

图12是三个方向上平面插值法示意图；

图13是Y轴方向上各插值平面标准方差。

具体实施方式

下面结合实例和附图对本发明作进一步详细说明。

现有针对前扰波进行建模在建模复杂，模型不精确等问题。例如位流法，一种较为常用的受油机前扰波建模方法，该方法是利用流体学方法进行理论推导，它将机头作为边界条件，把问题转化为直匀流和多个偶极子叠加的绕流问题。该方法可以直接推导出前扰波模型的解析解，但是公式复杂，计算量大，不适合工程实践。对于其他的拟合方法，也只是采用近似形状的函数，如指数函数幂函数等来表示。

本发明采用多项式进行函数拟合，多项式拟合函数形式简单，拟合精度高，通用性好，可以适用于任何机头形状前扰波的研究中，相比于三角函数、对数函数等其他函数形式，计算机对多项式的计算要快的多，也提高了实时性。

下面结合附图对本发明提供的用于空中加油的受油机前扰波建模方法的各实现步骤进行说明。

步骤一：受油机前扰波CFD仿真，得到受油机机头附近速度流场数据；

采用AutoCAD建立受油机的三维模型，模型为飞翼型布局，为了便于数值计算与仿真，选取的坐标系为CFD坐标系，即坐标原点位于机头，X轴正方向指向机尾，Y轴朝上，Z轴指向机身左侧，机头模型及插头位置如图1所示。本发明实施例中，受油机机头模型的尺寸为长4米，宽8米，高0.85米，加油插头的坐标位置为(1.4m，0.8m，-1.2m)。

选用PointWise软件来生成网格，基于点生成线，线生成面，最后合成体的思想，生成非结构网格，设定的求解器为ANSYS FLUENT。如图2所示，为了能足够表示出前扰波的影响范围，并且使求解的计算量尽量小，包含机头的求解范围为一个立方体，尺寸为6m×6m×10m。每条线上的节点数约为60个，最后生成的网格总数约为300万个。将PointWise生成的网格导入Fluent软件中并且检查网格质量，再返回PointWise中调整网格，直到满足求解要求为止。

结合空中加油的基本情况设置CFD的仿真求解条件。本发明实施例中，飞行高度为1500m，大气密度1.177kg/m³，焓值323608.9J/kg，设定的飞机飞行速度为210m/s，由于飞行速度较快，因此空气被假设为三维可压缩的粘性流体。选择的求解器为隐式耦合求解器，考虑热交换，用Spalart-Allmaras单方程模型进行湍流计算，选择适当的材料密度、粘度和热传导系数，边界条件分别设为压力远场和壁面。设定松弛因子为0.8，最后初始化流场的解，进行流场迭代计算，迭代130步后收敛在较小范围内。仿真结果包含各种物理量，包括在CFD坐标系三个轴向上的速度(v_x,v_y,v_z)，速度的大小v和压力，温度等。速度的大小v和静压力的求解结果如图3所示，其中速度是本次建模所需的物理量，图3中左图为CFD求解出的速度大小云图，右图为压力云图。

步骤二：研究一维的速度与位置的函数关系，即速度分别在X、Y、Z坐标轴上变化的规律；

在无人机空中加油对接过程中，锥套受到前扰波的气流扰动会漂浮不定，为了判断受油机前扰波对加油锥套的影响以及预测锥套的位置变化，需要对机头附近的速度进行分析。但Fluent软件仿真得到的数据量过于庞大，用查表法得到受油机头部附近的速度值计算量太大。因此需要给出受油机头部附近的三维速度与位置关系的解析表达式，根据表达式从而快速得到相应点的位置所对应的速度，便于控制器的设计。另外在对接过程中，加油锥套只会在附近一定范围之内运动，因此只需对局部区域进行重点分析，本发明实施例设置分析的区域范围为X轴上x取值为-0.2m～3.5m，Y轴上y取值0.5m～3.5m(无人机自身有一定厚度)，Z轴上z取值-3m～3m。

为了得到三维空间中速度与位置关系的解析表达式，首先要研究一维的速度与位置的函数关系，下文以X轴向上速度的分量v_x为例来进行说明。图4表示z＝-1.2平面上，即插头所在z平面，v_x的等值线。在此平面上，选取不同的y值，以此来研究v_x与x轴坐标值的对应关系。

受油机插头所在位置的y坐标值为0.8m，选取的y值范围是0.5m～3.5m，每隔0.2m进行一次数据采集，对导出的数据进行最小二乘拟合，选择的拟合函数形式为多项式。在满足精度的前提下，最终选用4次多项式来进行拟合，函数表达式如下，图5展示了其中三条多项式拟合曲线。

v_x＝c₄x⁴+c₃x³+c₂x²+c₁x+c₀

其中v_x为速度在x轴方向上的分量，x为坐标，(c₀,c₁,c₂,c₃,c₄)为待定系数。

用标准误差来决定拟合程度的好坏，标准误差表示为：

其中，ε_i为数据点的CFD计算值与拟合值之间的误差，当多项式的阶次越高，σ越小，考虑到计算精度与计算量之间的关系，选取一个满足计算精度的较为合适的多项式拟合阶次。

本发明实施例中，此处所有曲线中，标准误差最大值为0.11(m/s)，由此可以看出四阶多项式能够在满足精度要求的情况下很好的拟合出速度随x轴上位置变化的关系。在设置的z取值范围-3m～3m内，也可选取不同的z值，可按照上面过程，拟合得到对应z平面上v_x与x轴坐标值的对应关系，结论类似。

与此类似，在y＝0.8平面上拟合出v_x与z坐标的关系，由于无人机关于z＝0平面是对称的，因此只需拟合z＞0的数值，再沿直线z＝0做对称即可，对于z＞0只需5次多项式就能拟合出很好的效果。图6是y＝0.8平面速度v_x等值线，对称后的拟合函数图形如图7所示。同理，在x＝1.4平面拟合出的v_x与y坐标的函数关系如图9所示，图8为x＝1.4平面v_x等值线。

步骤三：在一维速度与位置的函数关系基础上，进行二维平面拟合；

掌握了v_x在xyz三个轴向上分布的规律之后，就可以对在机头横截面上的数据进行拟合，共有三种平面拟合方式，如图10所示。

以xz平面为例，由于在x轴向上v_x是关于x的四次多项式，在z轴向上v_x是关于z的五次多项式，推测出在xz平面上v_x与xz的基本函数关系如下：

v_x＝c₀₀+c₁₀x+c₀₁z+c₂₀x²+c₁₁xz+c₀₂z²+c₃₀x³+c₂₁x²z+

c₁₂xz²+c₀₃z³+c₄₀x⁴+c₃₁x³z+c₂₂x²z²+c₁₃xz³+

c₀₄z⁴+c₄₁x⁴z+c₃₂x³z²+c₂₃x²z³+c₁₄xz⁴+c₀₅z⁵

以上式为函数原型，根据最小二乘拟合出多个Y轴方向平面上的数据点，以平面y＝0.8为例，如图11所示，标准误差为0.14(m/s)，其他方向上的平面拟合与此类似。

步骤四：基于插值法进行三维解析表达式求解；

对于三维空间内的数据点，想要直接得到v_x与xyz坐标之间的解析表达式较为困难，为此选用插值法基于二维平面的解析表达式对三维空间内的数据进行计算，三种平面拟合方式对应着三种插值方向，如图12所示。以Y轴向上插值为例，选取多个xz平面进行二维平面拟合，得到多个平面拟合表达式，对于两个平面之间的点则根据这上下两个平面的计算结果进行线性插值。

当选取的平面过多，则平面拟合的表达式过多，不利于计算，如果选取的平面个数过少，则插值所带来的误差较大。因此就要结合插值平面的误差，合理的选择插值密度，在0.5≤y≤2范围内，若取插值密度为0.1，则各插值平面的误差曲线为图13。对于误差较大的位置插值密度应该较大，而误差较小的部分，插值密度可适当降低。在0.5m～1.0m之内插值密度选为0.05，即相邻平面间隔为0.05m，在1.0m～3.5m内将插值密度扩大为0.10。

速度与位置的基本计算公式如下：

式中v_xd和v_xu是由步骤三计算得到的相邻上下两个插值平面对应的速度值，y是待计算点的坐标值，y_d是下平面对应点的坐标值。

步骤五：对插值方向进行选择；

由上文可知，有三种形式的平面拟合，即xz平面拟合，xy平面拟合以及yz平面的拟合，分别对应三个插值方向，即Y轴，Z轴和X轴。选取插头附近空间区域x∈[0.5m，1.7m]，y∈[0.5m，2.0m]，z∈[-0.5m，-2m]为验证区域，对三种形式的拟合进行标准误差计算。结果如下：

插值方向	X	Y	Z
				标准误差(ms)	0.22	0.18	0.24

在Y轴方向上最小，可见，沿Y轴方向进行插值拟合是最为合适的，这也可以从直观上给出解释，由图9可知，在Y轴方向上，速度变化较为平缓，y＞2时速度基本不变，可看作线性变化，因此在图9的曲线上进行线性插值误差较小。

按照同样的方法，重复步骤二，步骤三和步骤四，对v_y,v_z速度分量进行拟合，可以得到v_y,v_z与三维位置之间关系的表达式。整体的函数关系可以表示为：

本方法与其他方法效果对比可见下表：

本发明对受油机前扰波进行数学建模，可以得到在受油机机头附近区域内任意一点所对应的流场速度，从而可以推算出前扰波对加油锥套的影响。在控制器中加入这一气动影响，对前扰波效应进行补偿，可以大大提高无人机空中加油自动对接的成功率。而且，本发明方法最终得到的是受油机机头附近流场速度与位置的解析表达式，可以快速、方便地给出锥套在任意一点所受的气动影响。

Claims (7)

1.一种用于空中加油的受油机前扰波建模方法，其特征在于，包括步骤：

步骤一：进行受油机前扰波CFD仿真，仿真获取三维空间中坐标点(x,y,z)对应的三维速度值(v_x,v_y,v_z)；

步骤二：用多项式拟合一维的速度与位置的函数关系；

对任一坐标轴上的速度v，用多项式分别拟合该速度v随每一轴上位置变化的关系；

多项式拟合表达式为：其中c_k为待定系数，n为最高次项，t＝x,y,z；

步骤三：在一维的速度与位置的函数关系基础上，进行二维的平面拟合；

对任一坐标轴上的速度v，在tr平面的多项式拟合表达式为：

其中，r＝x,y,z，t≠r；p,q分别为步骤二中确定的在t轴方向和r轴方向拟合多项式的阶数；c_pq为待定系数；

步骤四：基于步骤三得到的二维平面拟合表达式，选用插值法来得到三维空间中的解析表达式；三维空间的三个轴向对应三个插值方向，分别进行计算；

对任一坐标轴上的速度v，设步骤三计算得到t轴方向上的两个相邻的平行二维平面对应的v取值分别为v_d和v_u，则利用线性插值计算速度v与所述两个相邻的平行二维平面之间的点的关系，表达式如下：

v＝v_d+(v_u-v_d)(t-t_d)/d

其中，t是待计算点在t轴方向上的坐标值，t_d为其中一个二维平面在t轴方向上的坐标值，d为插值密度；

步骤五：选取三个插值方向中平面拟合的标准误差最小的方向作为最合适的插值方向，对速度和坐标进行插值拟合。

2.根据权利要求1所述的用于空中加油的受油机前扰波建模方法，其特征在于，所述的步骤一包括：采用AutoCAD建立受油机的三维模型，选取的坐标系为CFD坐标系，坐标原点位于机头，X轴正方向指向机尾，Y轴朝上，Z轴指向机身左侧；采用网格生成软件PointWise读取CAD文件并生成网格，将PointWise生成的网格导入Fluent软件中并检查网格质量，再返回PointWise中调整网格，直到满足求解要求；设置CFD的仿真求解条件，得出仿真结果，仿真结果包括三维空间内坐标点所对应的三维速度值。

3.根据权利要求1所述的用于空中加油的受油机前扰波建模方法，其特征在于，步骤二中所述的用多项式分别拟合速度v随每一轴上位置变化的关系，具体是：固定z值，选取不同y值，用多项式拟合速度v随x值变化的关系；固定y值，选取不同x值，用多项式拟合速度随z值变化的关系；固定x值，选取不同z值，用多项式拟合速度随y值变化的关系。

4.根据权利要求1或3所述的用于空中加油的受油机前扰波建模方法，其特征在于，所述的步骤二中，对于多项式拟合的好坏，用标准误差σ来判定，选取满足计算精度的多项式拟合阶次。

5.根据权利要求1或3所述的用于空中加油的受油机前扰波建模方法，其特征在于，所述的步骤二，根据加油锥套的运动范围，设置三维空间中各坐标轴向上的分析范围，在多项式拟合时，所需要固定的坐标值在该分析范围内选取。

6.根据权利要求1所述的用于空中加油的受油机前扰波建模方法，其特征在于，所述的步骤三中，根据加油锥套的运动范围，设置三维空间中各坐标轴向上的分析范围，在二维平面拟合时，在同一坐标轴向上，根据选择该坐标轴向上位于分析范围内的不同值，拟合得到不同的平行平面。

7.根据权利要求1所述的用于空中加油的受油机前扰波建模方法，其特征在于，所述的步骤四，三维空间内的解析表达式根据二维平面表达式利用插值法得到，并且利用线性插值计算时，插值密度根据各个插值平面拟合的误差曲线设置，拟合误差大的插值密度大，拟合误差小的插值密度小。