自主空中加油技术(Autonomous aerial refueling，AAR)兼具战略与战术双重价值，无论在军事还是民用方面都具有较大的应用潜力，受到了诸多国家的重视。软管-锥套式空中加油(软式加油)和伸缩管式空中加油(硬式加油)是两种常用的空中加油技术^[1]，并且各具优缺点。软式加油技术可以同时为多架无人机加油，并且成本低，方式简单。但是，在大气紊流情况下，系统不稳定，且传输速率低。与软式加油技术相比，硬式加油技术具有较快的传输速率并且易于与受油机连接。但是它需要的加油设备复杂且一次只能给一架无人机加油。随着无人机技术的迅速发展，AAR也成为国内外学者的研究热点。我国主要实施了软式加油方式，加油软管为一根长达十几米的柔性体，并受到加油机尾流场、大气紊流等因素的影响，对软管-锥套系统进行动力学建模和性能分析有助于提高与受油机安全对接成功率，同时对我国软式自主空中加油技术的发展有着积极意义^[2]。

目前，国内外学者对软管-锥套动态建模做了大量研究工作，取得了显著的成果。Zhu等人^[3]从材料力学角度提出了弹性横梁有限元模型，该方法克服了经典缆绳理论中处理旋转及大变形问题的困难，但是理论推导复杂，计算量大，难以满足实时性要求。Kamman等人^[4-5]基于多刚体运动学理论，提出一种恒定长度的多级串联“球杆”三维运动模型，并通过计算流体力学和风洞实验测试了锥套的气动特性。Vassberg等人^[6]采用计算流体力学法对KC-10加油机软管-锥套的气动特性进行了较全面的研究。胡孟权^[7-8]等将软管视为质点系，以静力平衡条件计算软管平衡拖拽位置，并研究了大气紊流对软管锥套运动的影响。王海涛等人^[9]在Kamman的基础上，提出了一种变长度的多级理想单摆系软管-锥套模型，并对软管“甩鞭”现象进行了分析。这些方法均采用传统的经典力学法建模，而凯恩方程法对多刚体运动学给出了统一的分析公式，而不必进行动态导数及微分运算，从而减小计算量，提高了运算速度。

本文根据多刚体动力学理论，提出了一种基于凯恩方程法的软管-锥套动态模型。该方法能够较简明地建立多自由度力学系统动力学方程，而不必进行内部约束力的分析。本文定义了描述系统状态的广义坐标及广义速率，导出了软管段位置多级递推公式和系统运动方程，估算了软管在加油机尾流、定常流和大气扰动下的气动载荷。数值仿真验证了所建模型的正确性，并进一步研究了系统的稳态性能和动态性能。

1 坐标系定义 1.1 建模假设

采用集中参数法思想，假设软管-锥套系统为多级串联“球杆”模型。其中软管由有限段等长的圆柱形光滑刚性杆组成，连接处由无摩擦的铰链相连，且各级刚性杆质量及所受外力均集中于连接处。锥套被视为固连于最后一级软管末端的质点，软管另一端与加油机吊舱相连，并随加油机运动。

1.2 坐标系及广义坐标定义

首先定义描述系统位姿所需的坐标系和状态量。

大地坐标系S_g：采用传统的北、东、地方向作为惯性系。

加油机航迹坐标系S_p：坐标系中X_p轴与加油机飞行速度方向重合一致，Z_p轴位于包含飞行速度在内的铅垂面，与X_p轴垂直并指向下方，Y_p轴垂直于O_pX_pZ_p平面，按右手定则确定。

杆件坐标系S_k：坐标原点位于各级杆件质心，X_k轴沿着杆件的对称轴方向指向加油机，Y_k轴垂直于对称面指向右侧，Z_k轴在对称面内指向杆件下方。

由S_p和S_k的关系，可以定义如下状态量描述各级杆件与加油机的相对位置。

相对俯仰角θ_k1：O_kX_k在O_pX_pZ_p平面的投影与航迹系O_pX_p轴的夹角，下倾为正。

相对偏航角θ_k2：O_kX_k与航迹系O_pX_pZ_p平面的夹角，左偏为正。

相对滚转角θ_k3：由于软管材料特性和外形结构约束使其绕管身中轴线的扭转运动可以忽略，即θ_k3≡0。

软管-锥套系统及坐标关系示意图如图 1所示。

由上述定义，可得由S_p系到S_k系的坐标变换矩阵为

${S_p}^k = \left[ {\matrix{ {cos{\theta _{k1}}cos{\theta _{k2}}} & {sin{\theta _{k2}}} & { - sin{\theta _{k1}}cos{\theta _{k2}}} \cr { - cos{\theta _{k1}}sin{\theta _{k2}}} & {cos{\theta _{k2}}} & {sin{\theta _{k1}}sin{\theta _{k2}}} \cr {sin{\theta _{k1}}} & 0 & {cos{\theta _{k1}}} \cr } } \right]$

(1)

根据前面建模假设及坐标系定义，设软管-锥套系统由n级刚性杆串联而成，第k级刚性杆在空间的位置可以用角度θ_k1和θ_k2加以描述。则系统的位姿可用2n个独立广义坐标来描述：(θ₁₁,θ₁₂,…,θ_k1,θ_k2,…,θ_n1,θ_n2)。相应的广义速率定义为：(${\dot \theta }$₁₁,${\dot \theta }$₁₂,… ,${\dot \theta }$_k1,${\dot \theta }$_k2,…,${\dot \theta }$_n1,${\dot \theta }$_n2) $ \buildrel \Delta \over = $(u₁,u₂,…,u_r,…u_2n)。

2 凯恩动力学方程推导

凯恩方程采用广义速率和偏速度描述系统^[10]，表达式为

${F_r} + {F_r}^* = 0$

(2)

式中：F_r和F_r^*分别为系统对第r个广义速率u_r的广义主动力和广义惯性力^[11]。

$\left\{ \matrix{ {F_r} = \mathop \sum \limits_{k = 1}^n \left( {{F_k}\cdot{V_{k,r}}^{} + {M_k}\cdot{\omega _{k,r}}^{}} \right) \hfill \cr {F_r}^* = \mathop \sum \limits_{k = 1}^n ({F_k}^*\cdot{V_{k,r}}^{} + {M_k}^*\cdot{\omega _{k,r}}^{}) \hfill \cr} \right.$

(3)

式中：V_k,r与ω_k,r为第k个刚体对第r个广义速率u_r的偏速度和偏角速度；F_k和M_k为作用于第k个刚体的合外力和合外力矩；F_k^*和M_k^*为第k个刚体的惯性力和惯性力矩。

2.1 系统运动学分析

在软管-锥套坐标系统关系如图 1所示。假设各级刚性杆的长度均为l_k，由第k-1级指向第k级的距离矢量为p_k，第k级刚性杆的空间位置向量在杆件坐标系S_k中可表示为(l_k,0,0) ^T，则p_k在S_p系下的坐标为

${p_k} = {\left( {{S_p}^k} \right)^T}\left[ \matrix{ - {l_k} \hfill \cr 0 \hfill \cr 0 \hfill \cr} \right] = {l_k}\left[ \matrix{ - {C_{k1}}{C_{k2}} \hfill \cr - {S_{k2}} \hfill \cr {S_{k1}}{C_{k2}} \hfill \cr} \right]$

(4)

式中：C_ki=cosθ_ki，S_ki=sinθ_ki(i=1，2) 。因此铰链k在空间的位置矢量r_k在S_p系下可表示为

${r_k} = {r_{k - 1}} + {p_k}$

(5)

第k级刚性杆在S_p系下的角速度和角加速度分别为

${\omega _k} = \left[ \matrix{ {{\dot \theta }_{k2}}{S_{k1}} \hfill \cr {{\dot \theta }_{k1}}{_{k2}} \hfill \cr \dot \theta {C_{k1}} \hfill \cr} \right]$

(6)

${\alpha _k} = {{d{\omega _k}} \over {dt}} = \left[ \matrix{ {{\ddot \theta }_{k2}}{S_{k1}} + {{\dot \theta }_{k1}}{{\dot \theta }_{k2}}{C_{k1}} \hfill \cr {{\ddot \theta }_{k1}} \hfill \cr {{\ddot \theta }_{k2}}{C_{k1}} - {{\dot \theta }_{k1}}{{\dot \theta }_{k2}}{S_{k1}} \hfill \cr} \right]$

(7)

则铰链k在空间的运动速度V_k及加速度a_k在S_p系下为

${V_k} = {V_{k - 1}} + {\omega _k} \times {p_k}$

(8)

${a_k} = {a_{k - 1}} + {\alpha _k} \times {p_k} + {\omega _k} \times ({\omega _k} \times {p_k})$

(9)

2.2 偏速度、偏角速度及惯性力、惯性力矩

将铰链k在空间的速度及角速度在S_p系下表示成广义速率的形式，即

${V_k} = {V_{k - 1}} + {l_k}\left[ \matrix{ {S_{k1}}{C_{k2}}_{k1} + {C_{k1}}{S_{k2}}{{\dot \theta }_{k2}} \hfill \cr - {C_{k2}}_{k2} \hfill \cr {C_{k1}}{C_{k2}}_{k1} - {S_{k1}}{S_{k2}}{{\dot \theta }_{k2}} \hfill \cr} \right] = \mathop \sum \limits_{i = 1,2} {V_{k,ki}}{\dot \theta _{ki}}$

(10)

${\omega _k} = \mathop \sum \limits_{i = 1,2} {\omega _{k,ki}}{{\dot \theta }_{ki}}$

(11)

式中：V_k,ki和ω_k,ki为铰链k对广义速率${\dot \theta }$_ki(i=1,2) 的偏速度和偏角速度向量。

${{V}_{k,ki}}=\frac{\partial {{V}_{k}}}{\partial ~{{{\dot{\theta }}}_{ki}}}=\left\{ \begin{align} & {{l}_{k}}{{\left[ {{S}_{k1}}{{C}_{k2}}0{{C}_{k1}}{{C}_{k2}} \right]}^{T}}~i=1 \\ & {{l}_{k}}{{\left[ {{C}_{k1}}{{S}_{~}}k2-{{C}_{k2}}-{{S}_{k1}}{{S}_{k2}} \right]}^{T}}i=2 \\ \end{align} \right.$

(12)

${{\omega }_{k,ki}}=\frac{\partial {{\omega }_{k}}}{{{\partial }_{~}}\dot{\theta }ki}=\left\{ \begin{align} & {{\left[ 010 \right]}^{T~}}i=1 \\ & {{\left[ {{S}_{k1}}0{{C}_{k1}} \right]}^{T}}~i=2 \\ \end{align} \right.$

(13)

那么，铰链k对广义速率${\dot \theta }$_j1,${\dot \theta }$_j2(j=1,2,…,n)的偏速度及偏角速度的通用表达式可表示为

${V_{k,j1}} = {{\partial {V_k}} \over {\partial {{\dot \theta }_{j1}}}} = \left\{ \matrix{ {l_k}{\left[ {{S_{j1}}{C_{j2}}0{C_{j1}}{C_{j2}}} \right]^T}k \ge j{\rm{ }} \hfill \cr {\left[ {000} \right]^T}k ＜ j \hfill \cr} \right.$

(14)

${V_{k,j2}} = {{\partial {V_k}} \over {\partial {{\dot \theta }_{j2}}}} = \left\{ \matrix{ {l_k}{\left[ {{C_{j1}}{S_{j2}} - {C_{j2}} - {S_{j1}}{S_{j2}}} \right]^T}k \ge j \hfill \cr {\left[ {000} \right]^T}k ＜ j \hfill \cr} \right.$

(15)

${\omega _{k,j1}} = {{\partial {\omega _k}} \over {\partial {{\dot \theta }_{j1}}}} = \left\{ \matrix{ {\left[ {010} \right]^T}k = j \hfill \cr {\left[ {000} \right]^T}k \ne j \hfill \cr} \right.$

(16)

${{\omega }_{k,j2}}=\frac{\partial {{\omega }_{k}}}{\partial {{{\dot{\theta }}}_{j2}}}=\left\{ \begin{align} & {{\left[ {{S}_{j1}}0{{C}_{j1}}^{~} \right]}^{T}}~k=j \\ & {{\left[ 000 \right]}^{T}}~k\ne j \\ \end{align} \right.$

(17)

下面计算在S_p系下铰链k的惯性力及惯性力矩。

${F_k}^* = - {m_k}{a_k} = - \mathop \sum \limits_{j = 1}^i {m_k}({V_{k,j1}}{{\ddot \theta }_{j1}} + {V_{k,j2}}{{\ddot \theta }_{j2}}) - {m_k}{\xi _k}$

(18)

式中

$\begin{align} & {{\xi }_{k}}=\underset{j=1}{\overset{k}{\mathop{\sum }}}\,{{l}_{k}}[-2{{S}_{j1}}{{S}_{j2}}{{_{j1}}_{j2}}+{{C}_{j1}}{{C}_{j2}}\left( {{{\dot{\theta }}}_{j1}}^{2}+{{{\dot{\theta }}}_{j2}}^{2} \right) \\ & {{S}_{j2}}{{{\dot{\theta }}}_{j2}}^{2}~-2{{C}_{j1}}{{S}_{j2}}{{{\dot{\theta }}}_{j1}}{{{\dot{\theta }}}_{j2}}-{{S}_{j1}}{{C}_{j2}}({{{\dot{\theta }}}_{j1}}^{2}+{{{\dot{\theta }}}_{j2}}^{2}){{]}^{T}} \\ \end{align}$

(19)

$\begin{align} & {{M}_{k}}^{*}=-{{J}_{k}}{{\alpha }_{k}}-{{\omega }_{k}}\times \left( {{J}_{k}}{{\omega }_{k}} \right)= \\ & -{{({{S}_{p}}^{k})}_{k}}T{{J}^{k}}{{S}_{p}}^{k}{{\omega }_{k,k1}}{{{\ddot{\theta }}}_{k1}}+{{\omega }_{k,k2}}{{{\ddot{\theta }}}_{k2}}-{{\eta }_{k}} \\ \end{align}$

(20)

$\begin{align} & {{\eta }_{k}}={{J}_{yk}}[2{{C}_{k1}}{{S}_{k2}}^{2}{{{\dot{\theta }}}_{k1}}{{{\dot{\theta }}}_{k2}}+{{S}_{k1}}{{S}_{k2}}{{C}_{k2}}{{{\dot{\theta }}}_{k1}}^{2} \\ & -2{{C}_{k2}}{{S}_{k2}}{{{\dot{\theta }}}_{k1}}{{{\dot{\theta }}}_{k2}}~-2{{S}_{k1}}{{S}_{k2}}^{2}{{{\dot{\theta }}}_{k1}}{{{\dot{\theta }}}_{k2}}+{{C}_{k1}}{{C}_{k2}}{{S}_{k2}}{{{\dot{\theta }}}_{k1}}^{2}] \\ \end{align}$

(21)

式中：J_k^k为第k级刚体件在S_k系下的惯量矩阵，由于各级刚性杆的直径d_k远小于其长度l_k且结构对称，因此设J_k^k=diag(0,J_yk,J_yk)。

2.3 软管锥套主动力分析

由受力分析可知，作用在铰链k上的主动力F_k包括第k级软管的重力m_kg和气动阻力D_k两部分。由于假设各级铰链所受外力均作用于质心位置，则主动力矩M_k=0。则有

${{F}_{k}}={{m}_{k}}g+\frac{{{D}_{k-1}}+{{D}_{k}}}{2}$

(22)

式中：气动阻力D_k由切向气动力D_t,k和法向气动力D_n,k组成^[12]。其中

${{D}_{t,k}}=\frac{1}{2}\rho {{V}_{t,k}}^{2}\pi {{d}_{k}}{{l}_{k}}{{C}_{t,k}}~$

(23)

${{D}_{n,k}}=\frac{1}{2}\rho {{V}_{n,k}}^{2}\pi {{d}_{k}}{{l}_{k}}{{C}_{n,k}}$

(24)

式中：ρ为空气密度；V_t,k=V_k-u_k为第k段软管的相对气流速度切向分量，u_k为铰链k处的定常流、加油机尾流和大气紊流等矢量和；V_n,k=V_k-u_k-V_t,k为第k段软管的相对气流速度法向分量；C_t,k和C_n,k分别为切向和法向摩擦力系数，其大小由当地气流雷诺数决定^[5]。

由于锥套安装在软管末端，同样受到重力和气动阻力的影响。则第n级铰链所受的主动力为

${{F}_{n}}=({{m}_{n}}+{{m}_{dro}})g+\frac{{{D}_{n}}}{2}+{{D}_{dro}}$

(25)

式中：m_dro为锥套质量，D_dro为锥套所受气动阻力。

${{D}_{dro}}~=\frac{1}{2}\rho {{({{V}_{n~}}-{{u}_{n~}})}^{2}}\left( \frac{\pi {{d}_{dro}}{{~}^{2}}}{4} \right){{C}_{dro}}~$

(26)

式中：d_dro为锥套直径，C_dro为锥套阻力系数，与锥套的几何外形有关^[9]。

2.4 系统动力学方程及仿真流程

将式(14~17) ，式(22~26) 所得的偏速度、偏角速度及主动力表达式代入式(3) ，得到对应于第r个广义速率u_r的广义主动力F_r；将式(18~21) 所得惯性力及惯性力矩代入式(3) ，得到对应于第r个广义速率u_r的广义惯性力F_r^*；且F_r和F_r^*均为标量表达式。最后根据凯恩动力学原理，得到2n个二阶非线性系统动力学方程如下

$\begin{align} & \underset{i=1}{\overset{n}{\mathop{\sum }}}\,\left( {{m}_{i}}{{a}_{i}}\cdot {{V}_{i,r}} \right)+\underset{i=1}{\overset{n}{\mathop{\sum }}}\,\left( {{J}_{i}}{{\alpha }_{i}}+{{\omega }_{i}}\times \left( {{J}_{i}}~{{\omega }_{i}} \right) \right)\cdot {{\omega }_{i,r}}+ \\ & {{m}_{dro}}{{a}_{n}}\cdot {{V}_{n,r}}+({{J}_{dro}}{{\alpha }_{n}}+{{\omega }_{n}}\times ({{J}_{dro}}{{\omega }_{n}}))\cdot {{\omega }_{n,r}}= \\ & \underset{i=1}{\overset{n-1}{\mathop{\sum }}}\,{{F}_{i}}\cdot {{V}_{i,r}}+{{F}_{n}}\cdot {{V}_{~}}n,r\text{ }r=1,2,\ldots ,2n \\ \end{align}$

(27)

把方程写成矩阵形式如下

$\left[ {\matrix{ {{H_{11}}} & {{H_{12}}} & \ldots & {{H_{1(2n)}}} \cr {{H_{21}}} & {{H_{22}}} & \ldots & {{H_{2(2n)}}} \cr {} & \vdots & \vdots & {} \cr {{H_{\left( {2n - 1} \right){\rm{ }}1}}} & {{H_{\left( {2n - 1} \right){\rm{ }}2}}} & \ldots & {{H_{\left( {2n - 1} \right){\rm{ }}(2n)}}} \cr {{H_{\left( {2n} \right)1}}} & {{H_{\left( {2n} \right)2}}} & \ldots & {{H_{\left( {2n} \right)(2n)}}} \cr } } \right]\left[ \matrix{ {{\ddot \theta }_{11}} \hfill \cr {{\ddot \theta }_{12}} \hfill \cr \vdots \hfill \cr {{\ddot \theta }_{n1}} \hfill \cr {{\ddot \theta }_{n2}} \hfill \cr} \right] = \left[ \matrix{ {R_1} \hfill \cr {R_2} \hfill \cr \vdots \hfill \cr {R_{2n - 1}} \hfill \cr {R_{2n}} \hfill \cr} \right]$

(28)

即

$\left[ H \right]\left[ {\ddot \theta } \right] = \left[ R \right] \Rightarrow \left[ {\ddot \theta } \right] = {\left[ H \right]^{ - 1}}\left[ R \right]$

(29)

给定加油机的飞行状态及软管的初始状态，采用如图 2所示的系统仿真流程来跟踪软管-锥套的动态过程。

3 数值仿真与算法性能评估

基于上述数学模型，对系统进行了Matlab仿真。将软管模型分为20段，软管拖拽点位于飞机右翼二分之一机翼处。仿真参数如表 1所示。

3.1 模型的正确性

为了验证基于凯恩方程的软管-锥套动态模型的可行性，采用Hallock-Burnham尾涡模型^[12]来近似模拟加油机尾涡流。设加油机在高度3 000 m，空速为100 m/s状态下平飞，观察软管-锥套系统的动态位置如图 3所示。

由于加油机尾涡流包含大量的不稳定风干扰，软管-锥套会持续飘摆，而Hallock-Burnham模型是针对尾涡流主要成分的简化模型，本质上是稳定的有旋风场，因此仿真结果最终趋于稳定。改变加油机的飞行条件，仿真验证系统都能在20 s内趋于稳定。该结论与文献^{[4, 9]}结论一致。

3.2 稳态性能

如图 1所示，定义锥套的平衡位置$D = {{{V_D}} \over {{L_H}}}$。在平稳大气中，对加油机在不同空速和高度条件下，锥套稳态阻力和平衡位置如图 4所示，软管-锥套平衡拖拽尾迹如图 5，6所示。

仿真结果显示，当加油机以100 m/s的空速在2 000～8 000 m的高空飞行时，加油锥套相对加油机拖拽点的稳定位置大概在2.78～4.66 m。即随着飞行高度的增加，由于空气密度减小，锥套气动阻力减小，所以锥套下沉量增加。当加油机在3 000 m高空以100～150 m/s的空速飞行时，加油锥套相对加油机拖拽点的稳定位置大概在3.04～1.77 m。即随着加油机飞行速度的增加，锥套受到的气动阻力增加，所以锥套下沉量减小。

3.3 动态性能

根据飞行试验报道，大气紊流对软管-锥套系统的影响很大。在晴朗的天气，当系统遇到轻微的紊流时，锥套就会发生高频无规则的振荡而无法稳定下来，从而导致加油任务的失败。在软管锥套系统中加入Dryden大气扰动模型^[13]，从t=10 s开始，产生一阵紊流强度为1.5 m/s，紊流尺度为530 m，时长为20 s的各项同性紊流。锥套在加油机不同飞行速度下(高度为3 000 m)受大气紊流影响的侧向和垂向位移以及第10段、第15段和第20段软管段的受扰运动分别如图 7，8所示。

由仿真图可见，锥套在大气扰动下发生不规则的振荡运动，并随着加油机飞行速度的增加，振荡幅度减小。这是因为在同一高度，所取紊流相同，而加油机速度增大时紊流相对自由来流的比值减小，从而对锥套的运动影响减小。同时，离拖拽点越远的软管段摆动越剧烈。如图 9所示，锥套在大气扰动下的运动轨迹杂乱无章而无法达到稳定位置，并且锥套阻力增加，该结论与报道的飞行试验结果一致^[6]。

4 结论

本文采用有限元法建立了基于凯恩方程的软管-锥套动态模型，该模型能够较为真实地反映软管-锥套的动力学特性，可以作为软管-锥套式空中加油研究的通用模型，并通过数值仿真验证了模型的正确性，分析了系统性能，得到如下结论：

(1) 平稳大气中，在不同的加油机飞行条件下，软管锥套系统能够在短时间内达到稳定状态。

(2) 加油机的飞行高度和飞行速度对软管锥套系统的稳定位置有一定的影响，且随着加油机速度增加，锥套阻力增加，下沉量减小；加油机高度增加，锥套阻力减小，下沉量增加。

(3) 大气紊流对软管锥套系统具有明显的干扰作用，同一高度时，加油机飞行速度增大，锥套受扰影响减小。受扰时，锥套发生杂乱无章的运动且阻力增大，将直接影响受油机安全对接操作。

在本文研究结论的基础上，后续可展开针对不确定性环境因素引起的误差进行在线建模参数实时调整研究以及需要对锥套系统进行进一步的增稳控制器研究以抑制紊流干扰。