本发明涉及飞机起落架理论建模技术领域，具体地，涉及一种滑跃甲板上起落过程中起落架对飞机作用力的确定方法。
背景技术：
飞机在移动平台上采用滑跃甲板起飞，就受载模式而言，它既不同于常规陆基起飞情况，也与弹射起飞大相径庭。由于是依靠滑跃甲板的倾斜曲面强制施加的载荷实现轨迹上升，这种情况的起落架地面载荷较常规陆基起飞滑跑情况更加难以确定。而现有技术中对起落架对飞机作用力的计算方法中至少存在以下问题：
1、现有方法多不考虑飞机的三维运动，运动过程中假设飞机在一个二维平面上运动。
2、现有方法不考虑甲板摇摆对飞机起落架力的影响。
3、在对飞机起落架轮胎摩擦力的计算中，现有方法均只考虑轮胎在滚动方向的摩擦力，使得其不能计算在甲板摇晃时，轮胎侧向力的改变。
技术实现要素：
针对现有技术中的缺陷，本发明的目的是提供一种滑跃甲板上起落过程中起落架对飞机作用力的确定方法。
根据本发明提供的一种滑跃甲板上起落过程中起落架对飞机作用力的确定方法，包括：
步骤1，对飞机、甲板、起落架轮胎分别建立坐标系，根据当前状态下飞机和甲板在地面惯性坐标系中的对地位置和对地速度获得飞机相对于所述甲板的位置和速度；
步骤2，根据所述飞机相对于所述甲板的位置判断所述飞机是否在所述甲板上，若所述飞机不在所述甲板上则令起落架对飞机作用力为0后结束，所述飞机在所述甲板上则进入下一步；
步骤3，确定飞机坐标系相对于地面惯性坐标系的第一旋转矩阵、甲板坐标系相对于地面惯性坐标系的第二旋转矩阵、甲板坐标系相对于飞机坐标系的第三旋转矩阵，根据飞机对地角速度和甲板对地角速度获得所述第一旋转矩阵的导数值、第二旋转矩阵的导数值、第三旋转矩阵的导数值；
步骤4，根据起落架运动约束方程、甲板拟合模型以及步骤3中获得的所述第一旋转矩阵的导数值、第二旋转矩阵的导数值、第三旋转矩阵的导数值，计算起落架压缩量和起落架压缩速度；
步骤5，根据所述起落架压缩量和起落架压缩速度以及起落架减震柱的刚度和阻尼获得所述起落架轮胎的正压力，
根据起落架轮胎的转向角确定所述起落架轮胎坐标系相对于飞机坐标系的第四旋转矩阵，根据飞机坐标系相对于甲板坐标系的旋转矩阵、所述第四旋转矩阵、起落架相对于甲板的速度在甲板坐标系的坐标以及lugre摩擦模型获得起落架轮胎滚动摩擦系数和侧向摩擦系数；
步骤6，根据所述起落架轮胎摩擦系数和所述起落架轮胎的正压力最终获得所述起落架对飞机作用力。
作为一种优化方案，所述步骤1进一步包括：
飞机和甲板的相对于地面惯性坐标系的位置向量分别为
飞机和甲板的相对于地面惯性坐标系的速度向量
作为一种优化方案，所述步骤3进一步包括：
获取所述飞机坐标系相对于地面惯性坐标系的第一旋转矩阵bri、甲板坐标系相对于与地面惯性坐标系的第二旋转矩阵cri、甲板坐标系相对于飞机坐标系的第三旋转矩阵crb；
根据飞机对地角速度获得飞机对地角速度在飞机坐标系中坐标的叉乘矩阵
根据甲板对地角速度获得甲板对地角速度在甲板坐标系中坐标的叉乘矩阵
则根据
作为一种优化方案，所述步骤4之前还包括：
确定所述起落架运动约束方程为pg，c(s)＝pac，c+crb(rg，b+sng，b)，在确定起落架运动约束方程后还包括对起落架运动约束方程pg，c(s)＝pac，c+crb(rg，b+sng，b)求导获得起落架速度约束方程
其中，pg，c(s)为起落架沿减震柱方向的延长线与甲板接交点在甲板坐标系中的坐标，pac，c为飞机坐标系原点相对于甲板坐标系原点的位置向量在甲板坐标系中的坐标，crb为所述第三旋转矩阵，rg，b为起落架减震柱在不受压缩时轮胎接地点位置相对于飞机坐标系原点的位置向量在飞机坐标系中的坐标，ng，b为起落架减震柱方向向量在飞机坐标系中的坐标，vg，c为起落架相对甲板的速度在甲板坐标系的坐标，vac，c为飞机相对于甲板坐标系原点的速度向量在甲板坐标系中的坐标，
作为一种优化方案，所述步骤4之前还包括确定甲板拟合模型为：
所述z＝ax2+bx+c的解为
其中，fc，c(pg，c)为甲板曲面点在甲板坐标系中的表达式，所述甲板包括平直部分和上翘部分，hc为所述上翘部在甲板坐标系中的高，l为所述平直部分在甲板坐标系中的长，l0为所述上翘部分在甲板坐标系中xcyc平面投影沿xc方向的长，l1为所述甲板在甲板坐标系中的宽，α为所述上翘部分的倾角。
作为一种优化方案，所述步骤4进一步包括：
将所述起落架运动约束方程、起落架速度约束方程、甲板拟合模型代入公式
作为一种优化方案，所述步骤5之前还包括根据所述甲板拟合模型确定起落架轮胎与甲板接触点处甲板法向量在甲板坐标系中的方程：
其中，fc，c(rg，c)为甲板曲面点在甲板坐标系中的表达式，所述甲板包括平直部分和上翘部分，hc为所述上翘部在甲板坐标系中的高，l为所述平直部分在甲板坐标系中的长，l0为所述上翘部分在甲板坐标系中xcyc平面投影沿xc方向的长，l1为所述甲板在甲板坐标系中的宽，α为所述上翘部分的倾角。
作为一种优化方案，所述步骤5中所述起落架轮胎的正压力的获得过程进一步包括
根据所述起落架轮胎与甲板接触点处甲板法向量在甲板坐标系中的方程计算
其中，ng，c为起落架减震柱方向向量在飞机坐标系中的坐标，ks为起落架减震柱的刚度、kr为起落架减震柱的阻尼，kvrebounce为回弹阻尼，kvimpact为压缩阻尼，s为起落架压缩量，
作为一种优化方案，所述步骤5中所述第四旋转矩阵为
其中，θg为起落架轮胎的转向角，brc为飞机坐标系相对于甲板坐标系的旋转矩阵，vg，c为起落架相对甲板的速度在甲板坐标系的坐标；
所述起落架轮胎滚动摩擦系数和侧向摩擦系数的获得过程进一步包括：
根据lugre摩擦模型获得所述起落架轮胎滚动摩擦系数和侧向摩擦系数的计算方程为
其中，σ(v)为起落架轮胎的摩擦系数，包括滚动方向分量和侧向分量，σ0、σ1、σd、σs都为lugre摩擦模型的时间常数，z为摩擦力系统内部状态系数，σstatic为静摩擦系数，σdynamic为滚动摩擦系数或滑动摩擦系数，v为起落架轮胎在起落架轮胎坐标系中xg轴方向和yg轴方向的速度。
作为一种优化方案，所述步骤6所述起落架对飞机作用力
其中，所述brg为飞机坐标系相对于起落架轮胎的旋转矩阵，所述σroll为起落架轮胎摩擦系数σ(v)的滚动方向分量，σroll为起落架轮胎摩擦系数σ(v)的侧向分量，n为起落架轮胎的正压力。
与现有技术相比，本发明具有如下的有益效果：
本方法考虑飞机三维运动，更加贴近真实的起飞状态。本方法考虑在甲板摇晃时飞机起落架作用力的变化。本方法使用lugre摩擦力模型，可以同时计算出轮胎的滚动摩擦力和侧向摩擦力。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单的介绍，显而易见，下面描述中的附图仅仅是本发明的一些实施例，对于本领域技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。附图中：
图1是滑跃甲板、飞机、地面坐标系及相关变量定义示意图；
图2是滑跃甲板上起落过程中起落架对飞机作用力的确定方法已知量输入和未知量输出示意图；
图3是滑跃甲板结构示意图；
图4是起落架轮胎坐标系示意图；
图5是起落架轮胎速度向量示意图；
图6是起落架轮胎受力示意图；
图7是可选的一种滑跃甲板上起落过程中起落架对飞机作用力的确定方法流程图。
具体实施方式
下文结合附图以具体实施例的方式对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，还可以使用其他的实施例，或者对本文列举的实施例进行结构和功能上的修改，而不会脱离本发明的范围和实质。
在本发明提供的一种滑跃甲板上起落过程中起落架对飞机作用力的确定方法的实施例中，包括：
步骤1，对飞机、甲板、起落架轮胎分别建立坐标系，根据当前状态下飞机和甲板在地面惯性坐标系中的对地位置和对地速度获得飞机相对于所述甲板的位置和速度；
步骤2，根据所述飞机相对于所述甲板的位置判断所述飞机是否在所述甲板上，若所述飞机不在所述甲板上则令起落架对飞机作用力为0后结束，所述飞机在所述甲板上则进入下一步；
步骤3，确定飞机坐标系相对于地面惯性坐标系的第一旋转矩阵、甲板坐标系相对于地面惯性坐标系的第二旋转矩阵、甲板坐标系相对于飞机坐标系的第三旋转矩阵，根据飞机对地角速度和甲板对地角速度获得所述第一旋转矩阵的导数值、第二旋转矩阵的导数值、第三旋转矩阵的导数值；
步骤4，根据起落架运动约束方程、甲板拟合模型以及步骤3中获得的所述第一旋转矩阵的导数值、第二旋转矩阵的导数值、第三旋转矩阵的导数值，计算起落架压缩量和起落架压缩速度；
步骤5，根据所述起落架压缩量和起落架压缩速度以及起落架减震柱的刚度和阻尼获得所述起落架轮胎的正压力，
根据起落架轮胎的转向角确定所述起落架轮胎坐标系相对于飞机坐标系的第四旋转矩阵，根据飞机坐标系相对于甲板坐标系的旋转矩阵、所述第四旋转矩阵、起落架相对于甲板的速度在甲板坐标系的坐标以及lugre摩擦模型获得起落架轮胎滚动摩擦系数和侧向摩擦系数；
步骤6，根据所述起落架轮胎摩擦系数和所述起落架轮胎的正压力最终获得所述起落架对飞机作用力。
步骤1建立坐标系的过程如图1和图4，包括：
建立i系为地面惯性坐标系，所述i系的原点oi为地面上的一点，以北向为xi轴方向，以东向为yi轴方向；
建立b系为飞机坐标系，所述b系的原点ob设于飞机质心，以飞机机头方向为xb轴方向，以飞机底部下方为zb轴方向；
建立g系为起落架轮胎坐标系，所述g系的原点og设于轮胎中心，以轮胎滚动方向为xg轴方向，沿减振柱向下延伸为zg轴方向；
建立c系为甲板坐标系，所述c系的原点oc设于甲板上，以飞机飞行方向为xc轴方向，沿法向朝下为zc轴方向。
本技术方案采用的矢量坐标表示方法如下：矢量a用符号
在脚标之后再加上一个数表示坐标的元素。如ac，i，2表示ac，i坐标的第2个元素。坐标系之间的旋转矩阵具有以下特性：
其中，e表示单位矩阵。在本技术方案中，矢量和其坐标随时间的导数表示为相应符号上加“.”，即：
对于相对旋转运动的坐标系，技术方案中定义旋转的相对角速度。
其中矩阵
根据式(1)和式(4)，同一个矢量在两个运动坐标系中随时间导数关系为：
本技术方案中一共需要使用4个坐标系：地面惯性坐标系、飞机机体坐标系、轮胎坐标系和甲板坐标系。如图1所示，地面惯性坐标系给出飞机和甲板绝对位置，用i系来表示，其原点oi为地面某一个点，以北东地方向为xiyizi轴方向。飞机机体坐标系和飞机本体固连，原点ob在飞机质心，xbybzb表示机体的前右下方向，表示飞机的位置和姿态，用b系来表示。轮胎坐标系和轮胎固连，原点og在轮胎中心，zg沿减振柱向下，xg和轮胎滚动方向重合，yg由xg和zg经右手螺旋得到，表示轮胎滚动方向和转向方向，用g系来表示。甲板坐标系和飞行甲板固连，表示甲板的运动和姿态，用c系来表示。根据图1，飞机和航母的相对于地面惯性坐标系的位置矢量分别用
对于位置向量
表1符号和坐标系定义
作为一种实施例，所述步骤1进一步包括：
飞机和甲板的相对于地面惯性坐标系的位置向量分别为
飞机和甲板的相对于地面惯性坐标系的速度向量
本方案讨论的方法是飞行仿真计算中起落架作用力计算模块。作为飞行仿真算法的一部分，本方案需要飞行仿真其他模块的数据输入来完成计算。输入量为飞机飞行姿态、速度、位置以及航母甲板运动速度和角速度，如图2所示。
作为一种实施例，所述步骤3进一步包括：
获取所述飞机坐标系相对于地面惯性坐标系的第一旋转矩阵bri、甲板坐标系相对于与地面惯性坐标系的第二旋转矩阵cri、甲板坐标系相对于飞机坐标系的第三旋转矩阵crb；
根据飞机对地角速度获得飞机对地角速度在飞机坐标系中坐标的叉乘矩阵
根据甲板对地角速度获得甲板对地角速度在甲板坐标系中坐标的叉乘矩阵
则根据
作为一种实施例，所述步骤4之前还包括：
确定所述起落架运动约束方程为pg，c(s)＝pac，c+crb(rg，b+sng，b)，在确定起落架运动约束方程后还包括对起落架运动约束方程pg，c(s)＝pac，c+crb(rg，b+sng，b)求导获得起落架速度约束方程
其中，pg，c(s)为起落架沿减震柱方向的延长线与甲板接交点在甲板坐标系中的坐标，pac，c为飞机坐标系原点相对于甲板坐标系原点的位置向量在甲板坐标系中的坐标，crb为所述第三旋转矩阵，rg，b为起落架减震柱在不受压缩时轮胎接地点位置相对于飞机坐标系原点的位置向量在飞机坐标系中的坐标，ng，b为起落架减震柱方向向量在飞机坐标系中的坐标，vg，c为起落架相对甲板的速度在甲板坐标系的坐标，vac，c为飞机相对于甲板坐标系原点的速度向量在甲板坐标系中的坐标，
作为一种实施例，所述步骤4之前还包括确定甲板拟合模型为：
所述z＝ax2+bx+c的解为
其中，fc，c(pg，c)为甲板曲面点在甲板坐标系中的表达式，所述甲板包括平直部分和上翘部分，hc为所述上翘部在甲板坐标系中的高，l为所述平直部分在甲板坐标系中的长，l0为所述上翘部分在甲板坐标系中xcyc平面投影沿xc方向的长，l1为所述甲板在甲板坐标系中的宽，α为所述上翘部分的倾角。
作为一种实施例，所述步骤4进一步包括：
将所述起落架运动约束方程、起落架速度约束方程、甲板拟合模型代入公式
作为一种实施例，所述步骤5之前还包括根据所述甲板拟合模型确定起落架轮胎与甲板接触点处甲板法向量在甲板坐标系中的方程：
其中，fc，c(pg，c)为甲板曲面点在甲板坐标系中的表达式，所述甲板包括平直部分和上翘部分，hc为所述上翘部在甲板坐标系中的高，l为所述平直部分在甲板坐标系中的长，l0为所述上翘部分在甲板坐标系中xcyc平面投影沿xc方向的长，l1为所述甲板在甲板坐标系中的宽，α为所述上翘部分的倾角。
作为一种实施例，所述步骤5中所述起落架轮胎的正压力的获得过程进一步包括
根据所述起落架轮胎与甲板接触点处甲板法向量在甲板坐标系中的方程计算
其中，ng，c为起落架减震柱方向向量在飞机坐标系中的坐标，ks为起落架减震柱的刚度、kv为起落架减震柱的阻尼，kvrebounce为回弹阻尼，kvimpact为压缩阻尼，s为起落架压缩量，
作为一种实施例，所述步骤5中所述第四旋转矩阵为
其中，θg为起落架轮胎的转向角，rc为飞机坐标系相对于甲板坐标系的旋转矩阵，vg，c为起落架相对甲板的速度在甲板坐标系的坐标；
所述起落架轮胎滚动摩擦系数和侧向摩擦系数的获得过程进一步包括：
根据lugre摩擦模型获得所述起落架轮胎滚动摩擦系数和侧向摩擦系数的计算方程为
其中，σ(v)为起落架轮胎的摩擦系数，包括滚动方向分量和侧向分量，σ0、σ1、σd、σs都为lugre摩擦模型的时间常数，z为摩擦力系统内部状态系数，σstatic为静摩擦系数，σdynamic为滚动摩擦系数或滑动摩擦系数，v为起落架轮胎在起落架轮胎坐标系中xg轴方向和yg轴方向的速度。
作为一种实施例，所述步骤6所述起落架对飞机作用力
其中，所述brg为飞机坐标系相对于起落架轮胎的旋转矩阵，所述σroll为起落架轮胎摩擦系数σ(v)的滚动方向分量，σroll为起落架轮胎摩擦系数σ(v)的侧向分量，n为起落架轮胎的正压力。
以下根据图7所示的流程步骤为例实施上述滑跃甲板上起落过程中起落架对飞机作用力的确定方法。根据两物体的相对摇摆获得相应的旋转矩阵是本领域的现有技术，而本发明将旋转矩阵加入滑跃甲板上起落过程中起落架对飞机作用力的确定过程却是对现有方法的一种较大改进，能够将坐标系之间的相对变化计算在内，提高了最终结算结果的准确性。在航天领域中，任何一计算模型的误差都可能会对航空安全产生无法估计的影响，本发明的计算模型中考虑到细微的摇摆因素无疑能够大大提高飞机的安全性与稳定性。
通过步骤3计算出飞机和飞行甲板旋转矩阵随时间的导数的过程为：根据图1和2中坐标系的定义和图3中输入参数，本方案中旋转矩阵的导数由如下公式确定：
其中，根据表1第3和4行的定义有：
步骤4中所述起落架运动约束方程的确定过程为：轮胎和飞机的相对位置用
将式(11)表示在c系中，有：
pg，c＝pac，c+crb(rg，b+sng，b).(12)
将式(12)看成是s的函数，获得所述起落架运动约束方程：
pg，c(s)＝pac，c+crb(rg，b+sng，b).(13)
令
由式(14)可知，交点的移动速度vg，c由飞机的运动速度vac，c，飞机旋转运动
根据图3可以对甲板进行数学拟合获得甲板数学模型。为了求解起落架压缩量，需要建立甲板的数学模型。如图3所示，航母的甲板由上翘部分和平直部分组成。在本技术方案中，上翘部分用二次曲线来拟合甲板的形状，已知平直部分长l，滑跃甲板长l0，高h，滑跃角为α，上翘部分和平直部分的宽度相等为l1。甲板的方程为：
z＝ax2+bx+c.(15)
可以解出参数：
将上翘部分和平直部分组合，定义甲板的参数方程fc，c和甲板法线方程nc，c，其中脚标第一个c代表航母甲板，第二个c表示参数方程使用c系中xyz坐标来表示关系。得到甲板的参数方程为：
根据图1，对于轮胎与甲板接触点位置和相对运动速度计算过程为：当已知飞机和航母的相对运动以及甲板几何形状之后，就可以计算出起落架的压缩量s和压缩速度
fc，c(pg，c)＝0.(19)
根据式(18)甲板曲面法向量在c系中的坐标为
nc，c＝nc，c(pg，c).(20)
由式(8)、(13)和(14)再联合式(17)，期落架的压缩量s和压缩速度
其中
对于轮胎正压力和摩擦力的计算方法，起落架轮胎坐标系和轮胎受力情况的计算如图4、图5、图6所示，得出s和
其中ng，c为起落架减振柱方向向量在c系中坐标，
如图4所示，定义g系来表示轮胎的滚动方向和转向情况。g坐标系的x轴和y轴分别沿轮胎转动方向和垂直于轮胎转动方向。当飞机起落架有转向角θg时，g系和b系之间的旋转矩阵grb为：
如图5、图6所示，由于轮胎在g系的x和y方向有相对甲板的运动速度
vg，g＝grbbrcvg，c.(24)
其中根据图2中输入量的定义，brc为从甲板坐标系变换到机体坐标系的旋转矩阵，是已知输入量。根据式(24)计算出轮胎滚动速度vg，g，1，轮胎侧滑速度vg，g，2和轮胎的压缩速度vg，g，3。在lugre摩擦模型中，再获得vg，g，1和vg，g，2之后，就可以计算出轮胎滚动摩擦力和侧向摩擦力。用σ(v)表示摩擦系数。v表示在g系xy方向上的速度的分量即vg，g，1或vg，g，2，摩擦力系统的内部状态设为z，它们之间有如下关系
式(25)表示当相应方向上有速度v时，摩擦力内部状态要发生变化，用
本方法考虑飞机三维运动，更加贴近真实的起飞状态。本方法考虑在甲板摇晃时飞机起落架作用力的变化。本方法使用lugre摩擦力模型，可以同时计算出轮胎的滚动摩擦力和侧向摩擦力。
以上所述仅为本发明的较佳实施例，本领域技术人员知悉，在不脱离本发明的精神和范围的情况下，可以对这些特征和实施例进行各种改变或等同替换。另外，在本发明的教导下，可以对这些特征和实施例进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本申请的权利要求范围内的实施例都属于本发明的保护范围。