动量环式主动增稳加油锥套建模与飘摆抑制

2024-05-09朱喆黄江涛章胜李飞杜昕单恩光唐骥罡王春阳

航空学报 2024年7期

朱喆，黄江涛，章胜，李飞，杜昕，单恩光，唐骥罡，王春阳

中国空气动力研究与发展中心空天技术研究所，绵阳 621000

在大气紊流、加油机尾流、翼尖涡流、受油机头波等干扰下，空中加油软管-锥套组合体产生非定常、非线性、快时变、强耦合飘摆运动，对空中加油高精度对接和飞行安全带来不利影响［1-2］。因此，准确建模模拟软管-锥套动态特性，并发展相应的锥套飘摆抑制方法对提高空中加油效能具有重要意义。

为模拟软管-锥套动态特性，国内外提出了多种建模方法。Ro等［3-4］利用多体动力学理论来对加油软管锥套系统建模，基于有限元分析法将软管锥套近似为多级串联的球杆模型，建模方法直观简单。文献［5-7］通过扩展Hamition 原理建立了变长度软管-锥套模型并开展仿真分析。Dai等［8］通过系统辨识得到了变长度软管-锥套动力学模型。李明哲等［9］基于绝对节点坐标法建立了软管-锥套多体动力学模型，进一步提高了计算效率。王海涛等［10-11］基于集中参数法建立了多级串联单摆系软管-锥套运动模型，并对软管收放、外部扰动和甩鞭现象等特性进行了研究。针对本文研究提出的自旋动量环式主动增稳加油锥套方案，上述建模方法难以对动量环转动时产生的动量矩进行表征，需要针对性开展建模研究。

软管-锥套空中飘摆运动抑制方面，目前对软管-锥套空中飘摆运动抑制的研究主要集中于翼舵式主动增稳控制方案。García-Fogeda等［12］通过在软管与锥管连接处添加翼舵控制面，改善了软管锥管系统的动态特性。张进等［13-14］基于十字舵面式锥套方案进行飘摆抑制控制器设计，有效减小了软管-锥套在大气紊流中的波动幅度。孟中杰等［15］设计了一种舵面式自主空中加油试验锥套，用来进行软式空中加油风洞试验验证。吴玲等［16-17］开展了基于比例积分微分和线性二次调节器的翼舵式软管锥套主动控制策略研究。翼舵式主动增稳锥套方案存在控制逻辑直观、飘摆抑制明显等优势，但其在回收过程中存在与加油机吊舱口撞击的风险，对于不同尺寸的加油吊舱适用性较差，同时对隐身性能产生一定程度的破坏。

本文提出一种自旋动量环式主动增稳锥套方案，通过在锥套锥杯处集成2 个具有偏置角度的主动控制动量环实现飘摆抑制，通过动量环转动产生动量矩进行增稳控制，克服传统翼舵式主动增稳锥套舵面与吊舱口撞击风险，以及吊舱/锥套锥杯尺寸约束下可折叠式翼舵结构复杂等问题。对动量环转动时产生的动量矩进行表征，基于拉格朗日方程开展针对性建模研究，建立考虑动量环旋转的多级串联单摆系空中加油软管-锥套飘摆运动模型，对比分析加油机洗流场内模型的准确性。分析阵风及大气紊流干扰下软管-锥套的飘摆运动特性。提出基于动量环控制的锥套飘摆抑制控制方法并开展仿真研究。

1 自旋动量环式主动增稳锥套设计

设计一种自旋动量环式主动增稳空中加油锥套方案，其原理是在锥套锥杯处集成主动控制动量环，通过直流电机控制动量环的转速及旋转方向，利用动量环产生的动量矩及进动效应对锥套飘摆进行抑制。虽然单动量环方案结构更为精简，但是前期研究表明由于该种方案的欠驱动特性影响，难以同时对锥套纵向及侧向飘摆运动进行有效抑制。因此，本文进一步提出双偏置动量环方案，即在锥套锥杯处集成2 个具有偏置角度的主动控制动量环实现飘摆抑制。2 个偏置动量环均与锥杯横截平面成一定平面角布置，从而增强了对锥套飘摆运动的抑制能力。锥杯内布置D-GPS（Differential Global Position System，差分全球定位系统）传感器，实时测定锥套运动位置及速度，如图1 所示。自旋动量环式方案相对于传统翼舵式主动增稳锥套（见图2），结构更加紧凑，释放回收过程安全性更高，隐身性能更好。

图1 自旋动量环式主动增稳锥套Fig.1 Actively stabilized drogue with momentum rings

图2 传统翼舵式主动增稳锥套Fig.2 Traditional actively stabilized drogue with rudders

2 软管-锥套组合体建模

2.1 建模假设及坐标系定义

基于拉格朗日方程，建立考虑动量环转动的多级串联单摆系软管-锥套运动模型，并作出如下假设：

1）加油机在地面系下做定直匀速平飞。

2）模型抽象为逐级串联单摆系，各段软管与相应摆杆长度相同，质量和外力均集中于球形铰链。

3）模型不计软管扭转与弹性，软管长度不变。

4）软管-锥套系统承受重力、自由来流、大气扰动形成的空气阻力。

5）不考虑加油机姿态变化影响，拖曳点系相对于地平系无转动。

6）锥套轴向平行于自由来流方向。

坐标系O0-X0Y0Z0为加油机拖曳点坐标系，其中，X0为加油机地速方向，Z0为加油机重力方向，Y0、Z0由右手定则得到。坐标系O0-XYZ 为重力与无风条件下自由来流阻力合力场下的等效惯性坐标系，其中，X 为重力与无风条件下自由来流阻力的合力负方向，与加油机拖曳点坐标系X0夹角为γ，Y 与Y0一致，Z 方向由右手定则得到，与Z0夹角为γ。锥套与加油机拖曳点之间通过软管连接，软管长度为l，如果软管由n 级摆杆组成，第k 级摆杆质量为ml，锥套位于末级摆杆铰链，质量为md。第k 级摆杆铰链在等效惯性系O0-XYZ 中的空间坐标可通过其相对O0XZ、O0XY 平面的偏转角θk1、θk2加以描述，如图3 所示。锥套即第n 级摆杆铰链在等效惯性系O0-XYZ 中的空间坐标，可通过其相对XO0Z、XO0Y平面的偏转角θn1、θn2加以描述，如图4 所示。锥套锥杯处集成了2 个偏置动量环W1、W2，与YZ 平面的安装角分别为χ、-χ，与XY 平面安装角均为0°，旋转角速率分别为，如图5所示。

图3 软管-锥套系统模型Fig.3 Hose-drogue model

图4 锥套摆动模型Fig.4 Drogue pendulous motion model

图5 动量环布置Fig.5 Momentum rings arrangement

2.2 软管-锥套系统建模

令铰链k 的空间位置矢量为rk，则rk在惯性系O0-XYZ 中可表示为

式中：pk为铰链k-1 指向铰链k 的距离矢量；pk在惯性系O0-XYZ 中的坐标为

对式（1）求导，得到铰链k 的运动速度vk为

由式（3）可得

不考虑加油机姿态变化影响，对式（2）求一次、二次导数得到分别为

式中：pk，θki为pk关于θki的偏导数为θki的一阶导数为θki的二阶导数。

令铰链k 在惯性系O0-XYZ 中绕O0转动的角速度为ωk，则

则锥套和动量环在惯性系O0-XYZ 中绕O0转动的角速度ωD、ωW1、ωW2分别为

建立拉格朗日方程：

系统总动能由系统各质点平动动能TE，锥套转动动能TD，动量环转动动能TW1、TW2组成，即

考虑锥套在其体系下的转动惯量为IDB，两动量环在各自体系下的转动惯量分别为IWB1、IWB2，锥套体系至等效惯性系的旋转矩阵为SIB，动量环1、2 各自体系至等效惯性系的旋转矩阵分别为SIB1、SIB2，则锥套和动量轮在惯性系O0-XYZ 中下的转动惯量ID、IW1、IW2分别为

考虑将重力、自由来流以及大气扰动形成的空气阻力作为广义力，则系统具有的势能V=0。

根据拉格朗日函数得

分别取各级摆杆的摆角θk1、θk2，以及动量环转动角度β1、β2作为系统的广义坐标，坐标维数共计2n+2。

1）对于坐标qk1=θk1，θk1对应的广义力为

式中：Fi为重力Fgi及空气阻力Fdi在惯性系O0-XYZ 中的合力。

由式（9）可得

当k=1，2，…，n-1时，拉格朗日方程中TE关于θk1的偏导数及相关导数项具体形式为

拉格朗日方程中TD关于θk1的偏导数及相关导数项具体形式为

拉格朗日方程中TW关于θk1的偏导数及相关导数项具体形式为

当k=n时，拉格朗日方程中TE关于θn1的偏导数及相关导数项具体形式为

拉格朗日方程中TD关于θn1的偏导数及相关导数项具体形式为

拉格朗日方程中TW关于θn1的偏导数及相关导数项具体形式为

根据式（15）、式（17）～式（22）可得，当k=1，2，…，n-1时，铰链k关于θk1坐标的动力学方程为

当k=n时，加油锥套关于θn1坐标的动力学方程为

2）对于坐标qk2=θk2，θk2对应广义力为

根据式（9）～式（11）、式（14）可得，当k=1，2，…，n-1时，铰链k关于坐标θk2的动力学方程为

当k=n时，加油锥套关于坐标θn2的动力学方程为

3）对于坐标qj=βj，令动量轮转矩为τ，则β对应的广义力为

参考文献［18-19］，考虑直流电机传动比为a∶1，令电机的转矩常数为Km，电机的电动势常数为Ke，电机的电阻为R，则β 对应的广义力为

根据式（9），拉格朗日方程中L 关于βj的偏导数及相关导数项具体形式为

根据式（9）～式（11）、式（31）、式（32）可得，加油锥套关于坐标βj的动力学方程为

根据式（23）、式（24）、式（27）、式（28）、式（33）可得到锥套-软管组合体在球坐标系下飘摆运动方程，并可将方程表示为

3 干扰条件下锥套飘摆抑制

加油锥套受前机尾流、阵风、大气紊流等扰动影响产生的随机飘摆，扰动气动力作用于锥套绕拖曳点O0转动所产生的动量矩为MD，其可沿等效惯性系O0-XYZ 下沿O0Y 轴和O0Z 轴分解为MDY、MDZ。两动量环分别绕锥套旋转，产生的动量矩为MW1、MW2，可在XY、XZ 平面内分解为MWXY1、MWXY2、MWXZ1、MWXZ2。在XY 平面内，调整锥套转速及方向，可使MWXY1与MWXY2合动量矩具有抵抗MDY的分量MWY；在XZ 平面内，调整锥套转速及方向，可使MWXZ1与MWXZ2合动量矩产生抵抗MDZ的分量MWZ，如图6、图7 所示。

图6 XY 平面内锥套飘摆抑制Fig.6 Pendulous suppression of drogue in XY plane

图7 XZ 平面内锥套飘摆抑制Fig.7 Pendulous suppression of drogue in XZ plane

综合考量系统响应特性及工程应用实际，拟采用位于加油锥套锥杯内部的D-GPS 传感器实时测定锥套位置及速度，通过锥杯内直流电机控制动量环的转速及方向，调整动量轮产生的动量矩大小及方向，实现锥套飘摆运动抑制，控制框图如图8 所示。

图8 锥套飘摆抑制控制框图Fig.8 Pendulous suppression control block diagram

控制系统输入状态量为锥套在等效惯性系下横向飘摆位置XDY、纵向飘摆位置XDZ、横向飘摆速度VDY、纵向飘摆速度VDZ，以及动量环旋转速率。XDY、XDZ、VDY、VDZ由D-GPS 数据及坐标系转换矩阵SBI确定由旋转速率传感器得到，锥套稳定状态下位置XDYcmd、XDZcmd由无干扰条件下软管-锥套模型仿真及飞行测试实际确定。

采用PID 方法进行控制律设计，控制律包括位置误差比例项、速度阻尼项，侧向通道控制力δY和纵向通道控制力δZ可表示为

式中：kpxy、kpxz分别为侧向、纵向位置比例控制增益；kpvy、kpvz分别为侧向、纵向速度阻尼增益。

为实现侧向、纵向控制力δY、δZ，需要进一步将其转化为电机电压控指令u1、u2，根据工程实际限制，电压指令限幅为±60 V。由图6、图7 可知，当两动量环同向旋转时，将产生抑制锥套侧向飘摆的动量矩MWZ；当两动量环反向旋转时，将产生抑制锥套纵向飘摆的动量矩MWY。因此，实现锥套纵向飘摆抑制时，两动量环伺服机构产生的控制力矩相反，实现锥套侧向飘摆抑制时，两动量环伺服机构产生的控制力矩相同，故在电压生成指令中纵向控制力反馈采用不同符号，侧向控制力反馈采用相同符号。同时，控制指令中考虑引入动量环角速度阻尼项，以达到改善控制效果的目的。综上，控制系统输出量即电机电压指令u1、u2为

4 仿真验证及分析

4.1 仿真环境建立及模型准确性验证

令加油软管长度10 m、加油锥套直径0.45 m，加油锥套质量7 kg，两动量环在体轴系下的转动惯量均为IWBxx=0.035 kg ⋅m2、IWByy=0.009 kg ⋅m2、IWBzz=0.009 kg ⋅m2，动量环安装角χ=20°。加油机飞行高度1 km、飞行速度85 m/s，拖拽点位于加油机腹正下方，基于文中软管-锥套运动模型建立仿真环境，如图9所示。

图9 加油锥套飘摆运动仿真Fig.9 Pendulous motion simulation of drogue

根据文献［20］，采用基于毕奥-萨伐尔定律、修正后的亥姆霍兹马蹄涡模型建立前机洗流场，与漩涡中心线距离为r 的风速大小表示为

式中：Γ 为涡旋强度；rc为涡核半径；v 为黏度参数；τ 为涡龄。

分别开展无控状态下基于拉格朗日方程的加油锥套飘摆运动仿真、基于集中参数法［10］的加油锥套飘摆运动仿真，YZ 平面内锥套飘摆轨迹如图10 所示，锥套飘摆直角坐标、极坐标对比如图11、图12 所示。可以看出，基于2 种方法建立的软管锥套飘摆运动轨迹一致，进一步对比验证了方法的准确性。

图10 锥套飘摆轨迹对比Fig.10 Comparison of drogue motion track

图11 锥套直角坐标对比Fig.11 Comparison of drogue rectangular coordinate

图12 锥套极坐标对比Fig.12 Comparison of drogue polar coordinate

4.2 阵风干扰下飘摆抑制仿真

空中加油对接过程中，阵风扰动将对锥套稳态位置带来极大影响。为对阵风扰动条件下，锥套飘摆及增稳情况进行分析，在Y、Z 方向上分别添加阵风干扰，阵风强度为6 m/s，阵风干扰时间为0.5 s。开展无控状态下和增稳状态下锥套飘摆运动仿真，锥套飘摆轨迹如图13 所示。

图13 阵风干扰下锥套Y、Z 方向无控/增稳运动轨迹对比Fig.13 Comparison of drogue Y，Z motion track with control/non-control under gust interference

可以看出，在阵风干扰下，加油锥套将产生周期飘摆运动。无控状态下，Y 方向上摆动幅度最大为0.68 m，Z 方向上摆动幅度最大为0.67 m，锥套受到干扰后16 s 收敛至平衡状态；增稳状态下，Y方向上摆动幅度最大为0.47 m，Z 方向上摆动幅度最大为0.40 m，干扰后5 s 收敛至平衡状态，锥套飘摆运动得到快速抑制。

4.3 大气紊流干扰下飘摆抑制仿真

大气紊流将使锥套飘摆运动呈现出非定常、非线性、快时变、强耦合的系统特征。大气紊流形成原因复杂，如地形诱导、风切边、热交换等，本文采用Dryden 模型［21］仿真不同速度下大气紊流对软管-锥套运动的影响，其纵向与横向频谱函数表达式为

式中：Ω 为空间频率；Lu、Lv、Lw为大气紊流尺度；σu、σv、σw为大气紊流强度。

添加Dryden 模型模拟风场干扰，分别开展无控和增稳状态加油锥套飘摆运动仿真，Y、Z 方向上锥套飘摆轨迹如图14 所示，动量环转速如图15所示，动量环电机电压如图16 所示，根据工程实际限制，电机电压限幅为±60 V。

图14 紊流干扰下锥套Y、Z 方向无控/增稳运动轨迹对比Fig.14 Comparison of drogue Y，Z motion track with control/non-control under turbulence interference

图15 动量环转速Fig.15 Rotation rate of momentum rings

图16 动量环电机电压Fig.16 Motor voltage of momentum rings

可以看出，在大气紊流的干扰下加油锥套将产生随机飘摆运动，无控状态下，Y 方向上最大摆动幅度大小为0.89 m，平均飘摆幅度大小为0.23 m；Z 方向上最大摆动幅度大小为0.73 m，平均飘摆幅度大小为0.18 m。开启增稳模式后，Y方向上最大摆动幅度大小为0.27 m，平均飘摆幅度大小为0.07 m；Z 方向上最大摆动幅度大小为0.18 m，平均飘摆幅度大小为0.06 m。增稳状态下，Y 方向上平均摆动幅度降低了69.57%，Z 方向上平均摆动幅度降低了66.67%，锥套飘摆幅度得到明显抑制。

锥套在紊流干扰下的高频抖振给空中加油精确对接带来了慢动态追踪快动态问题，进一步对无控/增稳状态下锥套飘摆运动数据开展频谱分析，进行快速傅里叶变换，结果如图17 所示。

图17 Y、Z 方向锥套飘摆频谱分析Fig.17 Spectrum analysis of drogue Y，Z motion

可以看出，无控状态下，锥套抖振能量主要集中于0.47～1.10 Hz频率范围内，该范围内Y 方向最大幅值强度达到0.16 m，Z 方向最大幅值强度达到0.11 m；开启增稳模式后，锥套在大于0.48 Hz范围的高频抖振得到明显削弱，整体幅值强度降低至0.05 m 以下，有效削弱了锥套高频动态响应，为空中加油精确对接控制实现创造了良好条件。