一种基于多刚体模型的跳伞员坠落运动仿真方法

摘要

本发明公开了一种基于多刚体模型的跳伞员坠落运动仿真方法本发明基于人体多刚体模型首次实现了跳伞员非刚体模型下的离机坠落运动仿真；在坠落运动过程中，跳伞员可根据需要调整头部和四肢相对于躯干的位置和姿态，进而改变整体下落轨迹和姿态；本发明的方法已经跳伞员伞降训练模拟系统中得到应用，相比传统的仿真方法，在真正意义上实现了跳伞员坠落运动轨迹和复杂姿态的仿真。

说明书

技术领域

本发明涉及飞行力学计算机仿真领域，尤其涉及一种基于多刚体模型的跳伞员坠落运动仿真方法。

背景技术

跳伞员离机坠落运动仿真有两种方法：一是空气动力学仿真方法；二是飞行力学仿真方法。采用空气动力学方法时，仿真跳伞员的离机坠落运动是一个极其复杂的计算流体力学问题，可以采用ANSYS FLUENT等专业软件，通过建立跳伞员的三维实体模型，基于实体模型进行动态仿真计算。这种方法的困难在于：就跳伞员实体模型来说，由于人体本身并不具有刚性结构，身体的每个部位如头、四肢等都是可以自由活动的，从而造成三维实体模型的不规则性和多样性，若要进行完全仿真，就必须建立各种姿态下的跳伞员模型并进行大量计算，需要花费大量的计算时间和人力、物力、财力。因此，这类方法只能进行基于跳伞员刚体模型的单一姿态下的仿真，对于基于跳伞员非刚体模型的复杂姿态下的仿真，当前还没有见到实例。采用飞行力学方法时，同样由于人体本身的非刚性，难以建立复杂姿态下的空气动力参数，仿真也只能在单一姿态下进行。

对于坠落运动的仿真，文献“人-伞系统自由坠落运动规律分析”(见《空军航空大学学报》2013年第3期)进行了人-伞系统的自由坠落运动仿真，但不涉及对跳伞员姿态的仿真；文献“单兵伞降运动轨迹计算问题研究”(见《桂林空军学院学报》2010年第5期)进行了单兵伞的自由坠落运动仿真，但也没涉及对跳伞员姿态的仿真。在运动仿真中，通常可将人体视为多刚体系统，例如Hanavan 15刚体模型、Yeadon 11刚体模型等。文献“飞机弹射救生过程人椅运动仿真”(见《航空计算技术》2009年第2期)建立了飞行员3刚体上肢模型，进行了人椅系统弹射后的飞行员上肢运动仿真。

由于在离机坠落运动过程中，跳伞员需要根据不同情况调整自身姿态，这种调整主要是依靠头部和四肢相对于躯干的位置和姿态变化来实现的。在这种情形下，由于跳伞员模型的复杂性和多变性，而难以利用传统的方法进行完整的运动仿真，必须寻找新的方法。

发明内容

有鉴于此，本发明提供了一种基于多刚体模型的跳伞员坠落运动仿真方法，实现了跳伞员非刚体模型下的离机坠落运动仿真。

本发明具有如下有益效果：

本发明基于人体多刚体模型首次实现了跳伞员非刚体模型下的离机坠落运动仿真；在坠落运动过程中，跳伞员可根据需要调整头部和四肢相对于躯干的位置和姿态，进而改变整体下落轨迹和姿态；本方法已经跳伞员伞降训练模拟系统中得到应用，相比传统的仿真方法，在真正意义上实现了跳伞员坠落运动轨迹和复杂姿态的仿真。

附图说明

图1为本发明中跳伞员人体的11块刚体模型；

图2为本发明中跳伞员坐标系参考姿态；

其中，1-躯干、2-骨盆、3-头部、4-左大臂、6-右大臂，5-左小臂、7-右小臂，8-左大腿、10-右大腿、9-左小腿、11-右小腿、3′-颈，4′-左肩、6′-右肩、5′-左肘、7′-右肘、8′-左臀、10′-右臀、9′-左膝、11′-右膝。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

(1)跳伞员人体空气动力学模型(简称跳伞员模型)简化：

首先，进行如下假设，

跳伞员躯干和骨盆简化为椭圆柱，头部简化为椭球体，大臂、小臂、大腿、小腿简化为 椭圆台；

忽略腕关节和踝关节的相对运动，将小臂和手作为一个整体，将小腿和脚作为一个整体；

颈、肩和臀部简化为具有固定回转中心的球铰链，铰接点位于各自的重心，且具有3个自由度；

肘关节、膝关节的铰接点分别位于各关节的重心位置。

其次，在假设条件下，将跳伞员模型简化为由11块刚体用铰链连接的多刚体系统(图1)，包括：躯干1，骨盆2，头部3，左4、右6大臂，左5、右7小臂，左8、右10大腿，左9、右11小腿等。每块刚体的物理特征均包括特征面积、特征长度以及重心位置、姿态角等。各刚体之间采用依附的级联方式联系在一起，每个铰链，即每个关节(颈3′，左4′、右6′肩，左5′、右7′肘，左8′、右10′臀，左9′、右11′膝)均可通过3轴旋转(滚转、俯仰、偏航)设置成包括坠落姿态在内的任意合理姿态。

如此，跳伞员坠落运动过程中受到的空气动力可表示为其中，F(粗体表示矢量，下同)为跳伞员总的空气动力，F_i为身体各部位空气动力，ΔF为模型简化以及各部位互扰造成的空气动力误差。

在整个自由坠落过程中，ΔF较小，将其忽略。从而，跳伞员受到的总的空气动力为>

这是整个仿真的基础。

(2)为了描述跳伞员模型在空间中的位置和姿态，以图2所示的姿态为基准，建立参考坐标系。在图中，跳伞员双臂、双腿并拢，且大小臂和大小腿分别形成一条直线，躯干、头、骨盆形成一条直线，双臂轴线、双腿轴线与躯干、头、骨盆轴线平行，跳伞员身体上下水平、左右对称。从而，建立坐标系如下：

跳伞员体轴坐标系：坐标系原点取在跳伞员重心；x轴为跳伞员身体的左右对称中心线，从脚部指向头部为正；y轴与跳伞员身体的左右对称面垂直，指向跳伞员身体左侧；z轴根据 x轴和y轴通过右手螺旋定则确定；用英文字母b表示，也称全局体轴系。

各刚体体轴坐标系：坐标系原点取在各刚体重心；x_i轴为各刚体的左右对称中心线，躯干、骨盆坐标系指向颈部为正，头部坐标系指向头顶为正，大臂坐标系指向肩部为正，小臂坐标系指向肘关节为正，大腿坐标系指向臀部为正，小腿坐标系指向膝关节为正；y_i轴与各刚体左右对称面垂直，指向各刚体左侧；z_i轴根据x_i轴和y_i轴通过右手螺旋定则确定；用英文字母b_i表示，也称局部体轴系。显然，在跳伞员模型基准姿态下，各刚体体轴坐标系的坐标轴方向与全局体轴系的各坐标轴方向相同，只是将坐标原点移到了各刚体重心。

跳伞员风轴坐标系：坐标系原点取在跳伞员重心；以空速相反方向作为x1轴；z1轴在跳伞员身体左右对称平面内，与x1轴垂直，向上为正；y1轴根据x1轴和z1轴通过右手螺旋定则确定；用英文字母w表示，也称全局风轴系。

各刚体风轴坐标系：坐标系原点取在各刚体重心；以空速相反方向作为x_i1轴；z_i1轴在各刚体左右对称平面内，与x_i1轴垂直，向上为正；y_i1轴根据x_i1轴和z_i1轴通过右手螺旋定则确定；用英文字母w_i表示，也称局部风轴系。显然，在跳伞员模型基准姿态下，各刚体风轴坐标系的坐标轴方向与全局风轴系的各坐标轴方向相同，只是将坐标原点移到了各刚体重心。

铰链坐标系：颈、肩坐标系与局部坐标中的躯干坐标系方向相同，肘坐标系与局部坐标中的大臂坐标系方向相同，膝坐标系与局部坐标中的大腿坐标系方向相同；所区别的只是将坐标原点移到各铰接点处；体轴系用英文字母b_i′表示，风轴系用英文字母w_i′表示。

跳伞员地面坐标系：坐标系原点取在跳伞员离机时刻重心在地面上的投影；x2轴指向此时速度矢量在地面投影的方向；z2轴与地面垂直，向上为正；y2轴根据x2轴和z2轴通过右手螺旋定则确定；用英文字母r表示，该系为基准坐标系，即惯性参考系。

使用下述变量记号：

右下标描述变量的对象或坐标关系，如表示头部风轴系相对于颈部风轴系的滚转 角；

右上标描述变量所在坐标系，如F₁^w表示躯干在全局风轴系中的气动力；

左上标描述变量导数所在的坐标系，如表示ω在体轴系的导数。

左下标描述变量的时间，如表示ω在t₀时刻的瞬时值。

给出下述坐标转换算法：

A201、跳伞员模型各刚体在各自局部风轴系的局部攻角α_i、侧滑角β_i和滚转角γ_i可通过其相对于铰接点的肢体定位以及铰接点攻角α_i′、侧滑角β_i′和滚转角γ_i′(根据坐标系的定义，铰接点攻角、侧滑角、滚转角与其所在刚体的对应角相同)的组合来计算。计算过程是一个角量的旋转变换问题，用四元数法计算如下：

定义跳伞员模型各刚体相对铰接点的滚转角俯仰角偏航角则

B₀＝cos(-α_i′/2)·cos(β_i′/2)·cos(-γ_i′/2)-sin(-α_i′/2)·sin(β_i′/2)·sin(-γ_i′/2)

B₁＝cos(-α_i′/2)·cos(β_i′/2)·sin(-γ_i′/2)+sin(-α_i′/2)·sin(β_i′/2)·cos(-γ_i′/2)

B₂＝cos(-α_i′/2)·sin(β_i′/2)·sin(-γ_i′/2)+sin(-α_i′/2)·cos(β_i′/2)·cos(-γ_i′/2)

B₃＝cos(-α_i′/2)·sin(β_i′/2)·cos(-γ_i′/2)-sin(-α_i′/2)·cos(β_i′/2)·sin(-γ_i′/2)

C₀＝A₀·B₀-A₁·B₁-A₂·B₂-A₃·B₃

C₁＝A₀·B₁+A₁·B₀+A₂·B₃-A₃·B₂

C₂＝A₀·B₂-A₁·B₃+A₂·B₀+A₃·B₁

C₃＝A₀·B₃+A₁·B₂-A₂·B₁+A₃·B₀

从而，各刚体在各自局部风轴系中的攻角α_i、侧滑角β_i和滚转角γ_i为：

β_i＝arcsin(2·(C₁·C₂+C₃·C₀))

为引用方便，将上述过程记为函数形式：

特别地，对于躯干和盆骨，w_1′＝w_2′＝w，即全局风轴系，且γ＝0。

A202、设攻角α、侧滑角β和滚转角γ，则采用下面的方法将全局(或局部)体轴系中的向量[x^b,y^b,z^b]^T转换到全局(或局部)风轴系[x^w,y^w,z^w]^T：体轴系与风轴系之间的坐标变换无论是在全局坐标系或是局部坐标系之间都是一样的，因此不区分；

A₀＝cos(-α/2)·cos(β/2)·cos(-γ/2)-sin(-α/2)·sin(β/2)·sin(-γ/2)

A₁＝cos(-α/2)·cos(β/2)·sin(-γ/2)+sin(-α/2)·sin(β/2)·cos(-γ/2)

A₂＝cos(-α/2)·sin(β/2)·sin(-γ/2)+sin(-α/2)·cos(β/2)·cos(-γ/2)

A₃＝cos(-α/2)·sin(β/2)·cos(-γ/2)-sin(-α/2)·cos(β/2)·sin(-γ/2)

B₀＝-x^b·A₁-y^b·A₂-z^b·A₃

B₁＝x^b·A₀+y^b·A₃-z^b·A₂

B₂＝-x^b·A₃+y^b·A₀+z^b·A₁

B₃＝x^b·A₂-y^b·A₁+z^b·A₀

从而，全局(或局部)风轴系向量[x^w,y^w,z^w]^T为：

x^w＝A₀·B₁-A₁·B₀-A₂·B₃+A₃·B₂

y^w＝A₀·B₂+A₁·B₃-A₂·B₀-A₃·B₁

z^w＝A₀·B₃-A₁·B₂+A₂·B₁-A₃·B₀

A203、设攻角α、侧滑角β和滚转角γ，则采用下面的方法将全局风轴系中的向量[x^w,y^w,z^w]^T转换到全局体轴系[x^b,y^b,z^b]^T：

A₀＝cos(γ/2)·cos(-β/2)cos(α/2)·+sin(γ/2)·sin(-β/2)·sin(α/2)

A₁＝-cos(γ/2)·sin(-β/2)·sin(α/2)+sin(γ/2)·cos(-β/2)·cos(α/2)

A₂＝cos(γ/2)·cos(-β/2)·sin(α/2)-sin(γ/2)·sin(-β/2)·cos(α/2)

A₃＝cos(γ/2)·sin(-β/2)·cos(α/2)+sin(γ/2)·cos(-β/2)·sin(α/2)

B₀＝-x^w·A₁-y^w·A₂-z^w·A₃

B₁＝x^w·A₀+y^w·A₃-z^w·A₂

B₂＝-x^w·A₃+y^w·A₀+z^w·A₁

B₃＝x^w·A₂-y^w·A₁+z^w·A₀

从而，全局体轴系向量[x^b,y^b,z^b]^T为：

x^b＝A₀·B₁-A₁·B₀-A₂·B₃+A₃·B₂

y^b＝A₀·B₂+A₁·B₃-A₂·B₀-A₃·B₁

z^b＝A₀·B₃-A₁·B₂+A₂·B₁-A₃·B₀

A204、设全局体轴系相对于地面坐标系的姿态角则Ω_b/r到四元数q之间的转换关系按以下方法计算：

将q′归一化得到q

简记作

地面坐标系到全局体轴系的转换矩阵C_b/r按以下方法计算：

四元数q到欧拉角Ω_b/r的转换关系按以下方法计算：

简记作f_Ω/q(q₀,q₁,q₂,q₃)；

(3)空气密度为ρ、空速大小为V_A/r，则来流动压头为q_∞＝0.5·ρ·V_A/r。从而，跳伞员受到的空气动力和力矩可按以下方法计算：

A301、令i＝1,2…,11，则跳伞员模型各刚体在局部风轴系中的空气动力为：

其中，分别为无量纲的空气动力系数在局部风轴系各坐标轴上的分量，其值与局部攻角α_i和侧滑角β_i有关，求α_i和β_i的算法由A201给出；S_i为各刚体的特征面积；

由此，跳伞员模型各刚体所受空气动力在全局风轴系中的阻力、升力和侧力分量为：

其中，γ_i为跳伞员模型各刚体在其局部风轴系中的滚转角，其算法由A201给出；

由于跳伞员坠落运动是在降落伞的稳定伞的作用下进行的，因此需要考虑稳定伞牵引作用。忽略过渡过程，则稳定伞牵引力的作用线总是与来流风向相同，视为阻力，即其中，为稳定伞的阻力系数，S_p为稳定伞的特征面积。

从而，计算出跳伞员模型在全局风轴系中的受到的气动力：

A302、令i＝1,2…,11，则跳伞员模型各刚体在局部风轴系中的空气动力矩为：

其中，分别为无量纲的空气动力矩系数在局部风轴系各坐标轴上的分量，其值与局部攻角α_i和侧滑角β_i有关，求α_i和β_i的算法由A201给出；L_i为各刚体的特征长度；

由此，跳伞员模型各刚体在全局风轴系中的滚转、俯仰和偏航力矩为：

其中，γ_i为跳伞员模型各刚体在其局部风轴系中的滚转角，其算法由A201给出；

令i＝1,2…,11，为跳伞员模型重心在全局风轴系中的位置坐标， 为各刚体重心在全局风轴系中的位置坐标，为各铰接点在全局风轴系中的位置坐标，Cg^w、和的值可由其在全局体轴系的相应的和转换得到，算法由A202给出。

于是，各刚体重力臂在全局风轴系中的分量：

令i＝1,2，

令i＝3,4,6,8,10，

令i＝5,7,9,11，

设稳定伞作用点在全局体轴系的坐标为用A202转换到全局风轴系，则稳定伞作用点的力臂计算如下：

其中，为稳定伞作用点在全局风轴系的坐标；

于是，跳伞员模型在全局风轴系中的滚转、俯仰和偏航力矩：

跳伞员在全局体轴系中受到的阻尼力矩为

其中，为跳伞员角速度在全局体轴系各坐标轴的分量；L_W为跳伞员模型特征宽度，L为特征长度，为特征面积；

从而，计算出跳伞员模型在全局体轴系中受到的气动力矩的合力距为：

其中，表示全局体轴系中受到的气动力矩，其各量由经A203转换而来；

根据算法A201～A204和算法A301、A302，就可以进行坠落运动的仿真了，其具体步骤包括：

步骤一：进行如下简化假设，

将大地当作静止平面，不考虑地球的曲率和旋转；

跳伞员地面坐标系为惯性参考系；

重力加速度为恒值，不随位置和高度而变化；

步骤二：初始化。

S201、建立跳伞员各坐标系：全局体轴系、全局风轴系，各刚体体轴系、各刚体风轴系，铰链坐标系，地面坐标系；

S202、初始化各参数：当地风速、跳伞员在地面坐标系的位置、地速、姿态角、角速度， 采用A204的转换矩阵C_b/r将位置、地速、姿态角和角速度转换到全局体轴系，采用A204的函数将姿态角转化为四元数；

跳伞员模型各刚体相对铰接点的姿态角，采用A202从局部体轴系转换到局部风轴系；

设置仿真步长Δt；

S203、设置仿真终止条件：设坠落持续时间为T、仿真时间为t_k，则t_k≥T时仿真就可以终止了；

步骤三：根据A204的转换矩阵C_b/r，将当地风速v_W/r转化到全局体轴系：

其中，为v_W/r在全局体轴系中的表示，为v_W/r在地面坐标系中的表示；>W/r在全局体轴系的三个分量，表示v_W/r在地面坐标系的三个分量；

步骤四：跳伞员在全局体轴系中的地速空速和当地风速三者形成速度三角形关系：

其中，表示S202获得的跳伞员在全局体轴系中的当前时刻地速；

从而，可计算跳伞员在全局体轴系中的空速：

其中，V_A/r表示空速的大小；

据此，计算跳伞员的全局攻角和侧滑角：

步骤五：根据步骤七中获得的跳伞员的全局攻角α和侧滑角β，将S202中得到的各刚体 相对姿态角代入到A201中跳伞员模型各刚体在各自局部风轴系下的攻角α_i、侧滑角β_i和滚转角γ_i的计算公式中，由此得到跳伞员模型各刚体攻角α_i、侧滑角β_i和滚转角γ_i，具体为：

S501、令i＝1,2，计算躯干、骨盆各角

S502、令i＝3,4,6，计算头、左右大臂各角

令i＝8,10，计算左右大腿各角

S503、令i＝5,7,9,11，计算左右小臂、左右小腿各角

步骤六：计算跳伞员在全局体轴系中的空气动力和力矩：

S601、根据步骤八得到的跳伞员模型各刚体局部攻角α_i、侧滑角β_i和滚转角γ_i，结合A301计算得到的跳伞员受到的全局风轴系中的空气动力F^w，并用A203的方法，将空气动力F^w转换到跳伞员全局体轴系，得F^b；

S602、根据步骤八得到的跳伞员模型各刚体局部攻角α_i、侧滑角β_i和滚转角γ_i，采用骤>b；

步骤七：计算当地风速的导数

其中，表示当地风速在跳伞员全局体轴系中的导数，表示当地风速在跳伞员地面坐标系中的导数；C_b/r按A204得出；“×”表示向量叉乘，在具体运算时，可通过下述方法将叉乘转化为点乘(下同)：

步骤八：计算跳伞员在全局体轴系内的空速加速度、角加速度以及地速度、角速度，

S801、跳伞员空速加速度和角加速度

其中，C_b/r按A204得出；m表示跳伞员的质量；g^r＝[0,0,-g_E]^T，g_E为重力加速度；J^b为转动惯量，其值为

S802、跳伞员地速度和角速度

其中，表示步骤四获得的当前时刻跳伞员空速；表示t_k+1时刻的风速；表示S202获得的当前时刻跳伞员角速度；

S803、以四元数方程表示的跳伞员角速度：

其中，为t_k+1时刻的角速度的三个分量；

步骤九：假设仿真时刻为t_k+1，采用欧拉数值积分方法计算跳伞员在全局体轴系内的位置和姿态，

S901、跳伞员在全局体轴系的位置为

其中，表示S202获得的当前时刻跳伞员位置，由步骤八获得；

S902、采用A204计算跳伞员在全局体轴系的姿态为

其中，和是t_k+1时刻的四元数的四个分量，且>

其中，表示S202中由姿态角转化的当前时刻四元数；由步骤八获得；

对于具体来说

当θ∈(-π/2,π/2)、时，

θ_b/r＝arcsin(-C_b/r(1,3))

当θ∈(-π,-π/2)∪(π/2,π)、时，

θ_b/r＝-π·sgn(C_b/r(1,3))-arcsin(-C_b/r(1,3))

步骤十：检查是否满足跳伞员坠落运动仿真的终止条件，若满足则转步骤十一；否则，利用步骤八获得的地速度和角速度以及步骤九获得的位置和四元数作为初始值，返回步骤三，即重复步骤三至步骤十的过程，直至满足终止条件；

步骤十一：输出跳伞员在每一仿真时刻的位置和姿态，仿真结束。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。