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

摘要

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

主权项

一种基于多刚体模型的跳伞员坠落运动仿真方法，其特征在于，包括如下步骤：步骤一：建立跳伞员人体空气动力学模型；首先，将跳伞员躯干和骨盆简化为椭圆柱，头部简化为椭球体，大臂、小臂、大腿、小腿简化为椭圆台；忽略腕关节和踝关节的相对运动，将小臂和手作为一个整体，将小腿和脚作为一个整体；颈、肩和臀部简化为具有固定回转中心的球铰链，铰接点位于各自的重心，且具有3个自由度；肘关节、膝关节的铰接点分别位于各关节的重心位置；其次，将跳伞员模型简化为由11块刚体用铰链连接的多刚体系统，包括：躯干1，骨盆2，头部3，左大臂4、右大臂6，左小臂5、右小臂7，左大腿8、右大腿10，左小腿9以及右小腿11；每块刚体的物理特征均包括特征面积、特征长度、重心位置以及姿态角；各刚体之间采用依附的级联方式联系在一起，其中，所述铰链，即为关节，包括颈3′，左肩4′、右肩6′，左肘5′、右肘7′，左臀8′、右臀10′，左膝9′以及右膝11′；最终，将跳伞员坠落运动过程中受到的空气动力表示为其中，F_i为第i块所述刚体所受到的空气动力，i＝1,2，...,11；步骤二：首先，建立参考坐标系，包括跳伞员体轴坐标系，即全局体轴系；各刚体的体轴坐标系，即局部体轴系；跳伞员风轴坐标系，即全局风轴系；各刚体的风轴坐标系，即局部风轴系；铰链坐标系以及跳伞员地面坐标系；然后确定参考坐标系之间的转换关系，包括如下具体步骤：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为：${\alpha }_i=arctan(2·(C_2·C_0-C_1·C_3)/(C_0^2+C_1^2-C_2^2-C_3^2))$β_i＝arcsin(2·(C₁·C₂+C₃·C₀))${\gamma }_i=-arctan(2·(C_1·C_0-C_3·C_2)/(C_0^2-C_1^2+C_2^2-C_3^2))$将攻角α_i、侧滑角β_i和滚转角γ_i的计算公式记为函数形式：A202、将全局体轴系和局部体轴系中的向量分别转换到全局风轴系和局部风轴系中；A203、将全局风轴系中的向量[x^w,y^w,z^w]^T转换到全局体轴系[x^b,y^b,z^b]^T中；A204、设全局体轴系相对于地面坐标系的姿态角获得姿态角Ω_b/r到四元数q之间的转换关系获得地面坐标系到全局体轴坐标系的转换矩阵C_b/r；获得四元数q到欧拉角Ω_b/r的转换关系f_Ω/q(q₀,q₁,q₂,q₃)；步骤三、获得跳伞员受到的空气动力和力矩，具体为：A301、跳伞员模型各刚体在局部风轴系中的空气动力为：$F_i^{w_i}=\left[\begin{array}{c}{F}_{ix}^{w_i}\\ F_{iy}^{w_i}\\ F_{iz}^{w_i}\end{array}\right]=\left[\begin{array}{c}{c}_{ix}^{w_i}({\alpha }_i,{\beta }_i)\\ c_{iy}^{w_i}({\alpha }_i,\beta )\\ c_{iz}^{w_i}({\alpha }_i,\beta )\end{array}\right]·q_{\infty }·S_i$其中，分别为无量纲的空气动力系数在局部风轴系各坐标轴上的分量，其值根据A201中的局部攻角α_i和侧滑角β_i获得；S_i为各刚体的特征面积；q_∞＝0.5·ρ·V_A/r表示来流动压头，ρ为空气密度，V_A/r为空速大小；则，跳伞员模型各刚体所受空气动力在全局风轴系中的阻力、升力和侧力分量为：$F_i^w=\left[\begin{array}{c}{F}_{ix}^w\\ F_{iy}^w\\ F_{iz}^w\end{array}\right]=\left[\begin{array}{ccc}-1& 0& 0\\ 0& -{\mathrm{sin\gamma }}_i& -{\mathrm{cos\gamma }}_i\\ 0& {\mathrm{cos\gamma }}_i& -{\mathrm{sin\gamma }}_i\end{array}\right]·F_i^{w_i}$其中，γ_i为A201中跳伞员模型各刚体在其局部风轴系中的滚转角；计算出跳伞员模型在全局风轴系中的受到的空气动力：$F^w={[F_x^w,F_y^w,F_z^w]}^T=\underset{i=1}{\overset{11}{\Sigma }}{F}_i^w+F_p^w$其中，表示稳定伞牵引力，为稳定伞的阻力系数，S_p为稳定伞的特征面积；A302、跳伞员模型各刚体在局部风轴系中的空气动力矩为：$M_i^{w_i}=\left[\begin{array}{c}{M}_{ix}^{w_i}\\ M_{iy}^{w_i}\\ M_{iz}^{w_i}\end{array}\right]=\left[\begin{array}{c}{m}_{ix}^{w_i}({\alpha }_i,{\beta }_i)\\ m_{iy}^{w_i}({\alpha }_i,{\beta }_i)\\ m_{iz}^{w_i}({\alpha }_i,{\beta }_i)\end{array}\right]·q_{\infty }·S_i·L_i$其中，分别为无量纲的空气动力矩系数在局部风轴系各坐标轴上的分量，其值根据A201中局部攻角α_i和侧滑角β_i计算获得，L_i为各刚体的特征长度；由此，跳伞员模型各刚体在全局风轴系中的空气动力矩为：$M_i^w=\left[\begin{array}{c}{M}_{ix}^w\\ M_{iy}^w\\ M_{iz}^w\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\mathrm{cos\gamma }}_i& -{\mathrm{sin\gamma }}_i\\ & {\mathrm{sin\gamma }}_i& {\mathrm{cos\gamma }}_i\end{array}\right]·M_i^{w_i}$其中，γ_i为A201中跳伞员模型各刚体在其局部风轴系中的滚转角；和分别为滚转、俯仰和偏航力矩；计算各刚体重力力臂在全局风轴系中的分量：即躯干1和盆骨2的重力力臂在全局风轴系中的分量为：头部3、左大臂4、右大臂6、左大腿8以及右大腿10的重力力臂在全局风轴系中的分量为：左小臂5、右小臂7、左小腿9以及右小腿11的重力力臂在全局风轴系中的分量为：其中，Cg^w为跳伞员模型重心在全局风轴系中的位置坐标，为各刚体重心在全局风轴系中的位置坐标，为各铰接点在全局风轴系中的位置坐标，和的值由各自在全局体轴系的相应位置坐标Cg^b、和转换得到，由A202给出；令为稳定伞作用点在全局体轴系的坐标，采用A202的方法将其转换到全局风轴系，则稳定伞作用点的力臂计算如下：$L_p^w={[L_{px}^w,L_{py}^w,L_{pz}^w]}^T={\mathrm{Cp}}^w-{\mathrm{Cg}}^w$其中，为稳定伞作用点在全局风轴系的坐标；则，跳伞员模型在全局风轴系中的滚转、俯仰和偏航力矩：$M_x^{\prime w}=\underset{i=1}{\overset{11}{\Sigma }}(-F_{ix}^w·L_{iy}^w-F_{iy}^w·L_{iz}^w+M_{ix}^w)$$M_y^{\prime w}=\underset{i=1}{\overset{11}{\Sigma }}(F_{ix}^w·L_{ix}^w-F_{iz}^w·L_{iz}^w+M_{iy}^w)-F_{pz}^w·L_{pz}^w$$M_z^{\prime w}=\underset{i=1}{\overset{11}{\Sigma }}(F_{iz}^w·L_{iy}^w+F_{iy}^w·L_{ix}^w+M_{iz}^w)+F_{pz}^w·L_{py}^w$跳伞员在全局体轴系中受到的阻尼力矩为：$M_d^b=-0.05·q_{\infty }·S·{[L_W·{\left({\omega }_{b/r}^b\right)}_x,L·{\left({\omega }_{b/r}^b\right)}_y,L·{\left({\omega }_{b/r}^b\right)}_z]}^T$其中，为跳伞员角速度在全局体轴系各坐标轴的分量；L_W为跳伞员模型特征宽度，L为跳伞员模型特征长度，为跳伞员模型特征面积；从而，计算出跳伞员模型在全局体轴系中受到的空气动力的合力矩为：$M^b=M^{\prime b}+M_d^b$其中，表示跳伞员模型在全局体轴系中受到的空气动力矩，其各分量由根据A203的方法获得；步骤四，假设：大地为静止平面，不考虑地球的曲率和旋转；跳伞员地面坐标系为惯性参考系；重力加速度为恒值，不随位置和高度而变化；步骤五：初始化：S202、初始化各参数：获得当地风速，跳伞员在地面坐标系中的位置、地速、姿态角以及角速度，采用A204的转换矩阵C_b/r，将跳伞员在地面坐标系位置、地速、姿态角以及角速度转换到全局体轴系，并采用A204的函数将姿态角转化为四元数；采用A202的方法，将跳伞员模型各刚体相对铰接点的姿态角从局部体轴系转换到局部风轴系；S203、设置仿真终止条件：设坠落持续时间为T、仿真时间为t_k，则t_k≥T时仿真终止；步骤六：根据A204的转换矩阵C_b/r，将当地风速v_W/r转化到全局体轴系：$v_{W/r}^b={[{\left(v_{W/r}^b\right)}_x,{\left(v_{W/r}^b\right)}_y,{\left(v_{W/r}^b\right)}_z]}^T=C_{b/r}·v_{W/r}^r=C_{b/r}·{[{\left(v_{W/r}^r\right)}_x,{\left(v_{W/r}^r\right)}_y,{\left(v_{W/r}^r\right)}_z]}^T$其中，为v_W/r在全局体轴系中的表示，为v_W/r在地面坐标系中的表示；表示v_W/r在全局体轴系的三个分量，表示v_W/r在地面坐标系的三个分量；步骤七：根据跳伞员在全局体轴系中的地速空速和当地风速形成的速度三角形关系：$v_{E/r}^b=v_{A/r}^b+v_{W/r}^b$其中，表示S202获得的跳伞员在全局体轴系中的当前时刻地速；从而，计算跳伞员在全局体轴系中的空速$v_{A/r}^b=\left[\begin{array}{c}{\left(v_{A/r}^b\right)}_x\\ {\left(v_{A/r}^b\right)}_y\\ {\left(v_{A/r}^b\right)}_z\end{array}\right]=\left[\begin{array}{c}{\left(v_{E/r}^b\right)}_x\\ {\left(v_{E/r}^b\right)}_y\\ {\left(v_{E/r}^b\right)}_z\end{array}\right]-C_{b/r}·\left[\begin{array}{c}{\left(v_{W/r}^r\right)}_x\\ {\left(v_{W/r}^r\right)}_y\\ {\left(v_{W/r}^r\right)}_z\end{array}\right],V_{A/r}={({\left(v_{A/r}^b\right)}_x^2+{\left(v_{A/r}^b\right)}_y^2+{\left(v_{A/r}^b\right)}_z^2)}^{1/2}$其中，V_A/r表示空速的大小；据此，计算跳伞员的全局攻角α和侧滑角β：$\alpha =arctan({\left(v_{A/r}^b\right)}_z/{\left(v_{A/r}^b\right)}_x),\beta =arcsin({\left(v_{A/r}^b\right)}_y/V_{A/r})$步骤八：根据步骤七中获得的跳伞员的全局攻角α和侧滑角β，将S202中得到的各刚体相对姿态角代入到A201中跳伞员模型各刚体在各自局部风轴系下的攻角α_i、侧滑角β_i和滚转角γ_i的计算公式中，由此得到跳伞员模型各刚体攻角α_i、侧滑角β_i和滚转角γ_i；步骤九：计算跳伞员在全局体轴系中的空气动力和力矩：S601、根据步骤八得到的跳伞员模型各刚体局部攻角α_i、侧滑角β_i和滚转角γ_i，结合A301计算得到的跳伞员受到的全局风轴系中的空气动力F^w，并用A203的方法，将空气动力F^w转换到跳伞员全局体轴系，得F^b；S602、根据步骤八得到的跳伞员模型各刚体局部攻角α_i、侧滑角β_i和滚转角γ_i，采用A302的方法，计算跳伞员受到的全局体轴系中的空气动力矩M^b；步骤十：计算当地风速的导数${}^b\stackrel{·}{v}_{W/r}^b=\frac{d}{dt}(C_{b/r}·v_{W/r}^r)=C_{b/r}^r·{}^r\stackrel{·}{v}_{W/r}^r-{\omega }_{b/r}^b\times C_{b/r}·v_{W/r}^r$其中，表示当地风速在跳伞员全局体轴系中的导数，表示当地风速在跳伞员地面坐标系中的导数；步骤十一：首先，获得跳伞员空速加速度和角加速度${}^b\stackrel{·}{v}_{A/r}^b=\frac{F^b}{m}{}^b\stackrel{·}{v}_{W/r}^b-{\omega }_{b/r}^b\times (v_{A/r}^b+v_{W/r}^b)+C_{b/r}·g^r$${}^b\stackrel{·}{\omega }_{b/r}^b={\left(J^b\right)}^{-1}·(M^b-{\omega }_{b/r}^b\times (J^b·{\omega }_{b/r}^b)$其中，m表示跳伞员质量；g^r＝[0,0,‑g_E]^T，g_E为重力加速度；J^b为转动惯量；其次，获得跳伞员的地速度和角速度其中，表示步骤七获得的当前时刻跳伞员空速；表示t_k+1时刻的风速；表示S202获得的当前时刻跳伞员角速度；Δt表示仿真步长；最后，以四元数方程表示跳伞员角速度：$\stackrel{·}{q}=\left[\begin{array}{c}{\stackrel{·}{q}}_0\\ {\stackrel{·}{q}}_1\\ {\stackrel{·}{q}}_2\\ {\stackrel{·}{q}}_3\end{array}\right]=\frac{1}{2}·\left[\begin{array}{cccc}0& -{\left({\omega }_{b/r}^b\right)}_x& -{\left({\omega }_{b/r}^b\right)}_y& -{\left({\omega }_{b/r}^b\right)}_z\\ {\left({\omega }_{b/r}^b\right)}_x& 0& {\left({\omega }_{b/r}^b\right)}_y& -{\left({\omega }_{b/r}^b\right)}_z\\ {\left({\omega }_{b/r}^b\right)}_y& -{\left({\omega }_{b/r}^b\right)}_z& 0& {\left({\omega }_{b/r}^b\right)}_x\\ {\left({\omega }_{b/r}^b\right)}_z& {\left({\omega }_{b/r}^b\right)}_y& -{\left({\omega }_{b/r}^b\right)}_x& 0\end{array}\right]·\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]$其中，和为t_k+1时刻的角速度的三个分量；步骤十二：假设仿真时刻为t_k+1，采用欧拉数值积分方法计算跳伞员在全局体轴系内的位置和姿态，具体为：S901、t_k+1时刻跳伞员在全局体轴系的位置为：其中，表示S202获得的当前时刻跳伞员位置，由步骤十一获得；S902、采用A204的方法，计算跳伞员在全局体轴系的姿态角为：其中，和是t_k+1时刻的四元数的四个分量，且其中，表示S202中由姿态角转化的当前时刻四元数；由步骤十一获得；步骤十三：检查是否满足跳伞员坠落运动仿真的终止条件，若满足则转步骤十四；否则,利用步骤十一获得的地速度和角速度以及步骤十二获得的位置和四元数作为初始值，返回步骤六；步骤十四：输出跳伞员在每一仿真时刻的位置和姿态，仿真结束。