CN105912776A - 一种基于多刚体模型的跳伞员坠落运动仿真方法 - Google Patents
一种基于多刚体模型的跳伞员坠落运动仿真方法 Download PDFInfo
- Publication number
- CN105912776A CN105912776A CN201610220029.XA CN201610220029A CN105912776A CN 105912776 A CN105912776 A CN 105912776A CN 201610220029 A CN201610220029 A CN 201610220029A CN 105912776 A CN105912776 A CN 105912776A
- Authority
- CN
- China
- Prior art keywords
- parachutist
- axis
- centerdot
- global
- angle
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000004088 simulation Methods 0.000 title claims abstract description 51
- 238000000034 method Methods 0.000 title claims abstract description 45
- 230000005484 gravity Effects 0.000 claims description 31
- 238000006243 chemical reaction Methods 0.000 claims description 16
- 210000000689 upper leg Anatomy 0.000 claims description 13
- 210000004197 pelvis Anatomy 0.000 claims description 11
- 230000001133 acceleration Effects 0.000 claims description 10
- 239000013598 vector Substances 0.000 claims description 10
- 230000009471 action Effects 0.000 claims description 9
- 210000000245 forearm Anatomy 0.000 claims description 9
- 210000001624 hip Anatomy 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 7
- 210000003127 knee Anatomy 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000000087 stabilizing effect Effects 0.000 claims description 6
- 210000003414 extremity Anatomy 0.000 claims description 5
- 230000009466 transformation Effects 0.000 claims description 5
- 230000002146 bilateral effect Effects 0.000 claims description 4
- 244000309466 calf Species 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 4
- 210000002310 elbow joint Anatomy 0.000 claims description 4
- 210000002683 foot Anatomy 0.000 claims description 4
- 210000000629 knee joint Anatomy 0.000 claims description 4
- 210000000544 articulatio talocruralis Anatomy 0.000 claims description 2
- 238000013016 damping Methods 0.000 claims description 2
- 238000013461 design Methods 0.000 claims description 2
- 230000010354 integration Effects 0.000 claims description 2
- 210000003857 wrist joint Anatomy 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 7
- 238000012549 training Methods 0.000 abstract description 2
- 230000036544 posture Effects 0.000 description 19
- 238000004422 calculation algorithm Methods 0.000 description 8
- 210000002414 leg Anatomy 0.000 description 6
- 230000009191 jumping Effects 0.000 description 3
- 238000005094 computer simulation Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 239000003381 stabilizer Substances 0.000 description 2
- 210000001364 upper extremity Anatomy 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 229960001948 caffeine Drugs 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- RYYVLZVUVIJVGH-UHFFFAOYSA-N trimethylxanthine Natural products CN1C(=O)N(C)C(=O)C2=C1N=CN2C RYYVLZVUVIJVGH-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/30—Circuit design
- G06F30/36—Circuit design at the analogue level
- G06F30/367—Design verification, e.g. using simulation, simulation program with integrated circuit emphasis [SPICE], direct methods or relaxation methods
Landscapes
- Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Microelectronics & Electronic Packaging (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
Abstract
本发明公开了一种基于多刚体模型的跳伞员坠落运动仿真方法本发明基于人体多刚体模型首次实现了跳伞员非刚体模型下的离机坠落运动仿真;在坠落运动过程中,跳伞员可根据需要调整头部和四肢相对于躯干的位置和姿态,进而改变整体下落轨迹和姿态;本发明的方法已经跳伞员伞降训练模拟系统中得到应用,相比传统的仿真方法,在真正意义上实现了跳伞员坠落运动轨迹和复杂姿态的仿真。
Description
技术领域
本发明涉及飞行力学计算机仿真领域,尤其涉及一种基于多刚体模型的跳伞员坠落运动仿真方法。
背景技术
跳伞员离机坠落运动仿真有两种方法:一是空气动力学仿真方法;二是飞行力学仿真方法。采用空气动力学方法时,仿真跳伞员的离机坠落运动是一个极其复杂的计算流体力学问题,可以采用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(粗体表示矢量,下同)为跳伞员总的空气动力,Fi为身体各部位空气动力,ΔF为模型简化以及各部位互扰造成的空气动力误差。
在整个自由坠落过程中,ΔF较小,将其忽略。从而,跳伞员受到的总的空气动力为
这是整个仿真的基础。
(2)为了描述跳伞员模型在空间中的位置和姿态,以图2所示的姿态为基准,建立参考坐标系。在图中,跳伞员双臂、双腿并拢,且大小臂和大小腿分别形成一条直线,躯干、头、骨盆形成一条直线,双臂轴线、双腿轴线与躯干、头、骨盆轴线平行,跳伞员身体上下水平、左右对称。从而,建立坐标系如下:
跳伞员体轴坐标系:坐标系原点取在跳伞员重心;x轴为跳伞员身体的左右对称中心线,从脚部指向头部为正;y轴与跳伞员身体的左右对称面垂直,指向跳伞员身体左侧;z轴根据 x轴和y轴通过右手螺旋定则确定;用英文字母b表示,也称全局体轴系。
各刚体体轴坐标系:坐标系原点取在各刚体重心;xi轴为各刚体的左右对称中心线,躯干、骨盆坐标系指向颈部为正,头部坐标系指向头顶为正,大臂坐标系指向肩部为正,小臂坐标系指向肘关节为正,大腿坐标系指向臀部为正,小腿坐标系指向膝关节为正;yi轴与各刚体左右对称面垂直,指向各刚体左侧;zi轴根据xi轴和yi轴通过右手螺旋定则确定;用英文字母bi表示,也称局部体轴系。显然,在跳伞员模型基准姿态下,各刚体体轴坐标系的坐标轴方向与全局体轴系的各坐标轴方向相同,只是将坐标原点移到了各刚体重心。
跳伞员风轴坐标系:坐标系原点取在跳伞员重心;以空速相反方向作为x1轴;z1轴在跳伞员身体左右对称平面内,与x1轴垂直,向上为正;y1轴根据x1轴和z1轴通过右手螺旋定则确定;用英文字母w表示,也称全局风轴系。
各刚体风轴坐标系:坐标系原点取在各刚体重心;以空速相反方向作为xi1轴;zi1轴在各刚体左右对称平面内,与xi1轴垂直,向上为正;yi1轴根据xi1轴和zi1轴通过右手螺旋定则确定;用英文字母wi表示,也称局部风轴系。显然,在跳伞员模型基准姿态下,各刚体风轴坐标系的坐标轴方向与全局风轴系的各坐标轴方向相同,只是将坐标原点移到了各刚体重心。
铰链坐标系:颈、肩坐标系与局部坐标中的躯干坐标系方向相同,肘坐标系与局部坐标中的大臂坐标系方向相同,膝坐标系与局部坐标中的大腿坐标系方向相同;所区别的只是将坐标原点移到各铰接点处;体轴系用英文字母bi′表示,风轴系用英文字母wi′表示。
跳伞员地面坐标系:坐标系原点取在跳伞员离机时刻重心在地面上的投影;x2轴指向此时速度矢量在地面投影的方向;z2轴与地面垂直,向上为正;y2轴根据x2轴和z2轴通过右手螺旋定则确定;用英文字母r表示,该系为基准坐标系,即惯性参考系。
使用下述变量记号:
右下标描述变量的对象或坐标关系,如表示头部风轴系相对于颈部风轴系的滚转 角;
右上标描述变量所在坐标系,如F1 w表示躯干在全局风轴系中的气动力;
左上标描述变量导数所在的坐标系,如表示ω在体轴系的导数。
左下标描述变量的时间,如表示ω在t0时刻的瞬时值。
给出下述坐标转换算法:
A201、跳伞员模型各刚体在各自局部风轴系的局部攻角αi、侧滑角βi和滚转角γi可通过其相对于铰接点的肢体定位以及铰接点攻角αi′、侧滑角βi′和滚转角γi′(根据坐标系的定义,铰接点攻角、侧滑角、滚转角与其所在刚体的对应角相同)的组合来计算。计算过程是一个角量的旋转变换问题,用四元数法计算如下:
定义跳伞员模型各刚体相对铰接点的滚转角俯仰角偏航角则
B0=cos(-αi′/2)·cos(βi′/2)·cos(-γi′/2)-sin(-αi′/2)·sin(βi′/2)·sin(-γi′/2)
B1=cos(-αi′/2)·cos(βi′/2)·sin(-γi′/2)+sin(-αi′/2)·sin(βi′/2)·cos(-γi′/2)
B2=cos(-αi′/2)·sin(βi′/2)·sin(-γi′/2)+sin(-αi′/2)·cos(βi′/2)·cos(-γi′/2)
B3=cos(-αi′/2)·sin(βi′/2)·cos(-γi′/2)-sin(-αi′/2)·cos(βi′/2)·sin(-γi′/2)
C0=A0·B0-A1·B1-A2·B2-A3·B3
C1=A0·B1+A1·B0+A2·B3-A3·B2
C2=A0·B2-A1·B3+A2·B0+A3·B1
C3=A0·B3+A1·B2-A2·B1+A3·B0
从而,各刚体在各自局部风轴系中的攻角αi、侧滑角βi和滚转角γi为:
βi=arcsin(2·(C1·C2+C3·C0))
为引用方便,将上述过程记为函数形式:
特别地,对于躯干和盆骨,w1′=w2′=w,即全局风轴系,且γ=0。
A202、设攻角α、侧滑角β和滚转角γ,则采用下面的方法将全局(或局部)体轴系中的向量[xb,yb,zb]T转换到全局(或局部)风轴系[xw,yw,zw]T:体轴系与风轴系之间的坐标变换无论是在全局坐标系或是局部坐标系之间都是一样的,因此不区分;
A0=cos(-α/2)·cos(β/2)·cos(-γ/2)-sin(-α/2)·sin(β/2)·sin(-γ/2)
A1=cos(-α/2)·cos(β/2)·sin(-γ/2)+sin(-α/2)·sin(β/2)·cos(-γ/2)
A2=cos(-α/2)·sin(β/2)·sin(-γ/2)+sin(-α/2)·cos(β/2)·cos(-γ/2)
A3=cos(-α/2)·sin(β/2)·cos(-γ/2)-sin(-α/2)·cos(β/2)·sin(-γ/2)
B0=-xb·A1-yb·A2-zb·A3
B1=xb·A0+yb·A3-zb·A2
B2=-xb·A3+yb·A0+zb·A1
B3=xb·A2-yb·A1+zb·A0
从而,全局(或局部)风轴系向量[xw,yw,zw]T为:
xw=A0·B1-A1·B0-A2·B3+A3·B2
yw=A0·B2+A1·B3-A2·B0-A3·B1
zw=A0·B3-A1·B2+A2·B1-A3·B0
A203、设攻角α、侧滑角β和滚转角γ,则采用下面的方法将全局风轴系中的向量[xw,yw,zw]T转换到全局体轴系[xb,yb,zb]T:
A0=cos(γ/2)·cos(-β/2)cos(α/2)·+sin(γ/2)·sin(-β/2)·sin(α/2)
A1=-cos(γ/2)·sin(-β/2)·sin(α/2)+sin(γ/2)·cos(-β/2)·cos(α/2)
A2=cos(γ/2)·cos(-β/2)·sin(α/2)-sin(γ/2)·sin(-β/2)·cos(α/2)
A3=cos(γ/2)·sin(-β/2)·cos(α/2)+sin(γ/2)·cos(-β/2)·sin(α/2)
B0=-xw·A1-yw·A2-zw·A3
B1=xw·A0+yw·A3-zw·A2
B2=-xw·A3+yw·A0+zw·A1
B3=xw·A2-yw·A1+zw·A0
从而,全局体轴系向量[xb,yb,zb]T为:
xb=A0·B1-A1·B0-A2·B3+A3·B2
yb=A0·B2+A1·B3-A2·B0-A3·B1
zb=A0·B3-A1·B2+A2·B1-A3·B0
A204、设全局体轴系相对于地面坐标系的姿态角则Ωb/r到四元数q之间的转换关系按以下方法计算:
将q′归一化得到q
简记作
地面坐标系到全局体轴系的转换矩阵Cb/r按以下方法计算:
四元数q到欧拉角Ωb/r的转换关系按以下方法计算:
简记作fΩ/q(q0,q1,q2,q3);
(3)空气密度为ρ、空速大小为VA/r,则来流动压头为q∞=0.5·ρ·VA/r。从而,跳伞员受到的空气动力和力矩可按以下方法计算:
A301、令i=1,2…,11,则跳伞员模型各刚体在局部风轴系中的空气动力为:
其中,分别为无量纲的空气动力系数在局部风轴系各坐标轴上的分量,其值与局部攻角αi和侧滑角βi有关,求αi和βi的算法由A201给出;Si为各刚体的特征面积;
由此,跳伞员模型各刚体所受空气动力在全局风轴系中的阻力、升力和侧力分量为:
其中,γi为跳伞员模型各刚体在其局部风轴系中的滚转角,其算法由A201给出;
由于跳伞员坠落运动是在降落伞的稳定伞的作用下进行的,因此需要考虑稳定伞牵引作用。忽略过渡过程,则稳定伞牵引力的作用线总是与来流风向相同,视为阻力,即其中,为稳定伞的阻力系数,Sp为稳定伞的特征面积。
从而,计算出跳伞员模型在全局风轴系中的受到的气动力:
A302、令i=1,2…,11,则跳伞员模型各刚体在局部风轴系中的空气动力矩为:
其中,分别为无量纲的空气动力矩系数在局部风轴系各坐标轴上的分量,其值与局部攻角αi和侧滑角βi有关,求αi和βi的算法由A201给出;Li为各刚体的特征长度;
由此,跳伞员模型各刚体在全局风轴系中的滚转、俯仰和偏航力矩为:
其中,γi为跳伞员模型各刚体在其局部风轴系中的滚转角,其算法由A201给出;
令i=1,2…,11,为跳伞员模型重心在全局风轴系中的位置坐标, 为各刚体重心在全局风轴系中的位置坐标,为各铰接点在全局风轴系中的位置坐标,Cgw、和的值可由其在全局体轴系的相应的和转换得到,算法由A202给出。
于是,各刚体重力臂在全局风轴系中的分量:
令i=1,2,
令i=3,4,6,8,10,
令i=5,7,9,11,
设稳定伞作用点在全局体轴系的坐标为用A202转换到全局风轴系,则稳定伞作用点的力臂计算如下:
其中,为稳定伞作用点在全局风轴系的坐标;
于是,跳伞员模型在全局风轴系中的滚转、俯仰和偏航力矩:
跳伞员在全局体轴系中受到的阻尼力矩为
其中,为跳伞员角速度在全局体轴系各坐标轴的分量;LW为跳伞员模型特征宽度,L为特征长度,为特征面积;
从而,计算出跳伞员模型在全局体轴系中受到的气动力矩的合力距为:
其中,表示全局体轴系中受到的气动力矩,其各量由 经A203转换而来;
根据算法A201~A204和算法A301、A302,就可以进行坠落运动的仿真了,其具体步骤包括:
步骤一:进行如下简化假设,
将大地当作静止平面,不考虑地球的曲率和旋转;
跳伞员地面坐标系为惯性参考系;
重力加速度为恒值,不随位置和高度而变化;
步骤二:初始化。
S201、建立跳伞员各坐标系:全局体轴系、全局风轴系,各刚体体轴系、各刚体风轴系,铰链坐标系,地面坐标系;
S202、初始化各参数:当地风速、跳伞员在地面坐标系的位置、地速、姿态角、角速度, 采用A204的转换矩阵Cb/r将位置、地速、姿态角和角速度转换到全局体轴系,采用A204的函数将姿态角转化为四元数;
跳伞员模型各刚体相对铰接点的姿态角,采用A202从局部体轴系转换到局部风轴系;
设置仿真步长Δt;
S203、设置仿真终止条件:设坠落持续时间为T、仿真时间为tk,则tk≥T时仿真就可以终止了;
步骤三:根据A204的转换矩阵Cb/r,将当地风速vW/r转化到全局体轴系:
其中,为vW/r在全局体轴系中的表示,为vW/r在地面坐标系中的表示;表示vW/r在全局体轴系的三个分量,表示vW/r在地面坐标系的三个分量;
步骤四:跳伞员在全局体轴系中的地速空速和当地风速三者形成速度三角形关系:
其中,表示S202获得的跳伞员在全局体轴系中的当前时刻地速;
从而,可计算跳伞员在全局体轴系中的空速:
其中,VA/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计算得到的跳伞员受到的全局风轴系中的空气动力Fw,并用A203的方法,将空气动力Fw转换到跳伞员全局体轴系,得Fb;
S602、根据步骤八得到的跳伞员模型各刚体局部攻角αi、侧滑角βi和滚转角γi,采用骤 步A302的方式,计算跳伞员受到的全局体轴系中的空气动力矩Mb;
步骤七:计算当地风速的导数
其中,表示当地风速在跳伞员全局体轴系中的导数,表示当地风速在跳伞员地面坐标系中的导数;Cb/r按A204得出;“×”表示向量叉乘,在具体运算时,可通过下述方法将叉乘转化为点乘(下同):
步骤八:计算跳伞员在全局体轴系内的空速加速度、角加速度以及地速度、角速度,
S801、跳伞员空速加速度和角加速度
其中,Cb/r按A204得出;m表示跳伞员的质量;gr=[0,0,-gE]T,gE为重力加速度;Jb为转动惯量,其值为
S802、跳伞员地速度和角速度
其中,表示步骤四获得的当前时刻跳伞员空速;表示tk+1时刻的风速;表示S202获得的当前时刻跳伞员角速度;
S803、以四元数方程表示的跳伞员角速度:
其中,为tk+1时刻的角速度的三个分量;
步骤九:假设仿真时刻为tk+1,采用欧拉数值积分方法计算跳伞员在全局体轴系内的位置和姿态,
S901、跳伞员在全局体轴系的位置为
其中,表示S202获得的当前时刻跳伞员位置,由步骤八获得;
S902、采用A204计算跳伞员在全局体轴系的姿态为
其中,和是tk+1时刻的四元数的四个分量,且
其中,表示S202中由姿态角转化的当前时刻四元数;由步骤八获得;
对于具体来说
当θ∈(-π/2,π/2)、时,
θb/r=arcsin(-Cb/r(1,3))
当θ∈(-π,-π/2)∪(π/2,π)、时,
θb/r=-π·sgn(Cb/r(1,3))-arcsin(-Cb/r(1,3))
步骤十:检查是否满足跳伞员坠落运动仿真的终止条件,若满足则转步骤十一;否则,利用步骤八获得的地速度和角速度以及步骤九获得的位置和四元数作为初始值,返回步骤三,即重复步骤三至步骤十的过程,直至满足终止条件;
步骤十一:输出跳伞员在每一仿真时刻的位置和姿态,仿真结束。
综上所述,以上仅为本发明的较佳实施例而已,并非用于限定本发明的保护范围。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (2)
1.一种基于多刚体模型的跳伞员坠落运动仿真方法,其特征在于,包括如下步骤:
步骤一:建立跳伞员人体空气动力学模型;
首先,将跳伞员躯干和骨盆简化为椭圆柱,头部简化为椭球体,大臂、小臂、大腿、小腿简化为椭圆台;忽略腕关节和踝关节的相对运动,将小臂和手作为一个整体,将小腿和脚作为一个整体;颈、肩和臀部简化为具有固定回转中心的球铰链,铰接点位于各自的重心,且具有3个自由度;肘关节、膝关节的铰接点分别位于各关节的重心位置;
其次,将跳伞员模型简化为由11块刚体用铰链连接的多刚体系统,包括:躯干1,骨盆2,头部3,左大臂4、右大臂6,左小臂5、右小臂7,左大腿8、右大腿10,左小腿9以及右小腿11;每块刚体的物理特征均包括特征面积、特征长度、重心位置以及姿态角;各刚体之间采用依附的级联方式联系在一起,其中,所述铰链,即为关节,包括颈3′,左肩4′、右肩6′,左肘5′、右肘7′,左臀8′、右臀10′,左膝9′以及右膝11′;
最终,将跳伞员坠落运动过程中受到的空气动力表示为其中,Fi为第i块所述刚体所受到的空气动力,i=1,2,...,11;
步骤二:首先,建立参考坐标系,包括跳伞员体轴坐标系,即全局体轴系;各刚体的体轴坐标系,即局部体轴系;跳伞员风轴坐标系,即全局风轴系;各刚体的风轴坐标系,即局部风轴系;铰链坐标系以及跳伞员地面坐标系;
然后确定参考坐标系之间的转换关系,包括如下具体步骤:
A201、针对跳伞员模型的各个刚体,通过各自相对于铰接点的肢体定位以及铰接点攻角αi′、侧滑角βi′和滚转角γi′,获得各刚体在各自局部风轴系的攻角αi、侧滑角βi和滚转角γi,具体为:
定义跳伞员模型各刚体相对铰接点的姿态角,即滚转角俯仰角偏航角则:
B0=cos(-αi′/2)·cos(βi′/2)·cos(-γi′/2)-sin(-αi′/2)·sin(βi′/2)·sin(-γi′/2)
B1=cos(-αi′/2)·cos(βi′/2)·sin(-γi′/2)+sin(-αi′/2)·sin(βi′/2)·cos(-γi′/2)
B2=cos(-αi′/2)·sin(βi′/2)·sin(-γi′/2)+sin(-αi′/2)·cos(βi′/2)·cos(-γi′/2)
B3=cos(-αi′/2)·sin(βi′/2)·cos(-γi′/2)-sin(-αi′/2)·cos(βi′/2)·sin(-γi′/2)
C0=A0·B0-A1·B1-A2·B2-A3·B3
C1=A0·B1+A1·B0+A2·B3-A3·B2
C2=A0·B2-A1·B3+A2·B0+A3·B1
C3=A0·B3+A1·B2-A2·B1+A3·B0
从而,各刚体在各自局部风轴系中的攻角αi、侧滑角βi和滚转角γi为:
βi=arcsin(2·(C1·C2+C3·C0))
将攻角αi、侧滑角βi和滚转角γi的计算公式记为函数形式:
A202、将全局体轴系和局部体轴系中的向量分别转换到全局风轴系和局部风轴系中;
A203、将全局风轴系中的向量[xw,yw,zw]T转换到全局体轴系[xb,yb,zb]T中;
A204、设全局体轴系相对于地面坐标系的姿态角获得姿态角Ωb/r到四元数q之间的转换关系
获得地面坐标系到全局体轴坐标系的转换矩阵Cb/r;
获得四元数q到欧拉角Ωb/r的转换关系fΩ/q(q0,q1,q2,q3);
步骤三、获得跳伞员受到的空气动力和力矩,具体为:
A301、跳伞员模型各刚体在局部风轴系中的空气动力为:
其中,分别为无量纲的空气动力系数在局部风轴系各坐标轴上的分量,其值根据A201中的局部攻角αi和侧滑角βi获得;Si为各刚体的特征面积;q∞=0.5·ρ·VA/r表示来流动压头,ρ为空气密度,VA/r为空速大小;
则,跳伞员模型各刚体所受空气动力在全局风轴系中的阻力、升力和侧力分量为:
其中,γi为A201中跳伞员模型各刚体在其局部风轴系中的滚转角;
计算出跳伞员模型在全局风轴系中的受到的空气动力:
其中,表示稳定伞牵引力, 为稳定伞的阻力系数,Sp为稳定伞的特征面积;
A302、跳伞员模型各刚体在局部风轴系中的空气动力矩为:
其中,分别为无量纲的空气动力矩系数在局部风轴系各坐标轴上的分量,其值根据A201中局部攻角αi和侧滑角βi计算获得,Li为各刚体的特征长度;
由此,跳伞员模型各刚体在全局风轴系中的空气动力矩为:
其中,γi为A201中跳伞员模型各刚体在其局部风轴系中的滚转角;和分别为滚转、俯仰和偏航力矩;
计算各刚体重力力臂在全局风轴系中的分量:
即躯干1和盆骨2的重力力臂在全局风轴系中的分量为:头部3、左大臂4、右大臂6、左大腿8以及右大腿10的重力力臂在全局风轴系中的分量为:
左小臂5、右小臂7、左小腿9以及右小腿11的重力力臂在全局风轴系中的分量为:
其中,Cgw为跳伞员模型重心在全局风轴系中的位置坐标,为各刚体重心在全局风轴系中的位置坐标,为各铰接点在全局风轴系中的位置坐标,和的值由各自在全局体轴系的相应位置坐标Cgb、和转换得到,由A202给出;
令为稳定伞作用点在全局体轴系的坐标,采用A202的方法将其转换到全局风轴系,则稳定伞作用点的力臂计算如下:
其中,为稳定伞作用点在全局风轴系的坐标;
则,跳伞员模型在全局风轴系中的滚转、俯仰和偏航力矩:
跳伞员在全局体轴系中受到的阻尼力矩为:
其中,为跳伞员角速度在全局体轴系各坐标轴的分量;LW为跳伞员模型特征宽度,L为跳伞员模型特征长度,为跳伞员模型特征面积;
从而,计算出跳伞员模型在全局体轴系中受到的空气动力的合力矩为:
其中,表示跳伞员模型在全局体轴系中受到的空气动力矩,其各分量由根据A203的方法获得;
步骤四,假设:大地为静止平面,不考虑地球的曲率和旋转;跳伞员地面坐标系为惯性参考系;重力加速度为恒值,不随位置和高度而变化;
步骤五:初始化:
S202、初始化各参数:获得当地风速,跳伞员在地面坐标系中的位置、地速、姿态角以及角速度,采用A204的转换矩阵Cb/r,将跳伞员在地面坐标系位置、地速、姿态角以及角速度转换到全局体轴系,并采用A204的函数将姿态角转化为四元数;
采用A202的方法,将跳伞员模型各刚体相对铰接点的姿态角从局部体轴系转换到局部风轴系;
S203、设置仿真终止条件:设坠落持续时间为T、仿真时间为tk,则tk≥T时仿真终止;
步骤六:根据A204的转换矩阵Cb/r,将当地风速vW/r转化到全局体轴系:
其中,为vW/r在全局体轴系中的表示,为vW/r在地面坐标系中的表示;表示vW/r在全局体轴系的三个分量,表示vW/r在地面坐标系的三个分量;
步骤七:根据跳伞员在全局体轴系中的地速空速和当地风速形成的速度三角形关系:
其中,表示S202获得的跳伞员在全局体轴系中的当前时刻地速;
从而,计算跳伞员在全局体轴系中的空速
其中,VA/r表示空速的大小;
据此,计算跳伞员的全局攻角α和侧滑角β:
步骤八:根据步骤七中获得的跳伞员的全局攻角α和侧滑角β,将S202中得到的各刚体相对姿态角代入到A201中跳伞员模型各刚体在各自局部风轴系下的攻角αi、侧滑角βi和滚转角γi的计算公式中,由此得到跳伞员模型各刚体攻角αi、侧滑角βi和滚转角γi;
步骤九:计算跳伞员在全局体轴系中的空气动力和力矩:
S601、根据步骤八得到的跳伞员模型各刚体局部攻角αi、侧滑角βi和滚转角γi,结合A301计算得到的跳伞员受到的全局风轴系中的空气动力Fw,并用A203的方法,将空气动力Fw转换到跳伞员全局体轴系,得Fb;
S602、根据步骤八得到的跳伞员模型各刚体局部攻角αi、侧滑角βi和滚转角γi,采用A302的方法,计算跳伞员受到的全局体轴系中的空气动力矩Mb;
步骤十:计算当地风速的导数
其中,表示当地风速在跳伞员全局体轴系中的导数,表示当地风速在跳伞员地面坐标系中的导数;
步骤十一:首先,获得跳伞员空速加速度和角加速度
其中,m表示跳伞员质量;gr=[0,0,-gE]T,gE为重力加速度;Jb为转动惯量;
其次,获得跳伞员的地速度和角速度
其中,表示步骤七获得的当前时刻跳伞员空速;表示tk+1时刻的风速;表示S202获得的当前时刻跳伞员角速度;Δt表示仿真步长;
最后,以四元数方程表示跳伞员角速度:
其中,和为tk+1时刻的角速度的三个分量;
步骤十二:假设仿真时刻为tk+1,采用欧拉数值积分方法计算跳伞员在全局体轴系内的位置和姿态,具体为:
S901、tk+1时刻跳伞员在全局体轴系的位置为:
其中,表示S202获得的当前时刻跳伞员位置,由步骤十一获得;
S902、采用A204的方法,计算跳伞员在全局体轴系的姿态角为:
其中,和是tk+1时刻的四元数的四个分量,且
其中,表示S202中由姿态角转化的当前时刻四元数;由步骤十一获得;
步骤十三:检查是否满足跳伞员坠落运动仿真的终止条件,若满足则转步骤十四;否则,利用步骤十一获得的地速度和角速度以及步骤十二获得的位置和四元数作为初始值,返回步骤六;
步骤十四:输出跳伞员在每一仿真时刻的位置和姿态,仿真结束。
2.如权利要求1所述的一种基于多刚体模型的跳伞员坠落运动仿真方法,其特征在于,所述跳伞员体轴坐标系的原点取在跳伞员重心;x轴为跳伞员身体的左右对称中心线,从脚部指向头部为正;y轴与跳伞员身体的左右对称面垂直,指向跳伞员身体左侧;z轴根据x轴和y轴通过右手螺旋定则确定;
各刚体体轴坐标系的原点取在各刚体重心;xi轴为各刚体的左右对称中心线,躯干、骨盆坐标系指向颈部为正,头部坐标系指向头顶为正,大臂坐标系指向肩部为正,小臂坐标系指向肘关节为正,大腿坐标系指向臀部为正,小腿坐标系指向膝关节为正;yi轴与各刚体左右对称面垂直,指向各刚体左侧;zi轴根据xi轴和yi轴通过右手螺旋定则确定;
跳伞员风轴坐标系的原点取在跳伞员重心;以空速相反方向作为x1轴;z1轴在跳伞员身体左右对称平面内,与x1轴垂直,向上为正;y1轴根据x1轴和z1轴通过右手螺旋定则确定;
各刚体风轴坐标系的原点取在各刚体重心;以空速相反方向作为xi1轴;zi1轴在各刚体左右对称平面内,与xi1轴垂直,向上为正;yi1轴根据xi1轴和zi1轴通过右手螺旋定则确定;
铰链坐标系中,颈和肩坐标系的各坐标轴与躯干坐标系各坐标轴方向相同,肘坐标系的各坐标轴与大臂坐标系的各坐标轴方向相同,膝坐标系的各坐标轴与大腿坐标系的各坐标轴方向相同;铰链坐标系的坐标原点为对应铰接点处;
跳伞员地面坐标系的原点取在跳伞员离机时刻重心在地面上的投影;x2轴指向此时速度矢量在地面投影的方向;z2轴与地面垂直,向上为正;y2轴根据x2轴和z2轴通过右手螺旋定则确定。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610220029.XA CN105912776B (zh) | 2016-04-07 | 2016-04-07 | 一种基于多刚体模型的跳伞员坠落运动仿真方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610220029.XA CN105912776B (zh) | 2016-04-07 | 2016-04-07 | 一种基于多刚体模型的跳伞员坠落运动仿真方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105912776A true CN105912776A (zh) | 2016-08-31 |
CN105912776B CN105912776B (zh) | 2019-03-08 |
Family
ID=56745804
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610220029.XA Expired - Fee Related CN105912776B (zh) | 2016-04-07 | 2016-04-07 | 一种基于多刚体模型的跳伞员坠落运动仿真方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105912776B (zh) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106446442A (zh) * | 2016-10-12 | 2017-02-22 | 北京理工大学 | 一种火星伞舱组合体的降落伞展开过程稳定性评估方法 |
CN106767781A (zh) * | 2016-11-29 | 2017-05-31 | 中国地质大学(武汉) | 跌落测试预埋传感器的六自由度运动轨迹数据处理方法 |
CN108268427A (zh) * | 2018-01-10 | 2018-07-10 | 中国地质大学(武汉) | 一种任意球进球概率分析方法、设备及存储设备 |
CN114202633A (zh) * | 2021-12-06 | 2022-03-18 | 北京理工大学 | 一种耦合空气动力学的跳台滑雪肌骨动力学模拟方法 |
CN114638444A (zh) * | 2022-05-19 | 2022-06-17 | 珠海翔翼航空技术有限公司 | 针对突发故障下的跳伞时间预测方法、系统、设备 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR20060071302A (ko) * | 2004-12-21 | 2006-06-26 | 한국전자통신연구원 | 사용자 인터페이스 설계 및 평가 시스템 및 손 상호 작용기반의 사용자 인터페이스 설계 및 평가 시스템 |
CN102520726A (zh) * | 2011-12-19 | 2012-06-27 | 南京航空航天大学 | 大攻角飞行状态下的大气攻角及侧滑角估计方法 |
CN102880732A (zh) * | 2011-12-28 | 2013-01-16 | 南京康尼机电股份有限公司 | 一种轨道交通车辆门系统动力学联合仿真分析方法 |
CN103921947A (zh) * | 2014-03-27 | 2014-07-16 | 中国科学院长春光学精密机械与物理研究所 | 跳伞模拟训练半实物仿真系统及其操作方法 |
KR101510058B1 (ko) * | 2014-06-30 | 2015-04-08 | 성결대학교 산학협력단 | 이산 사건 시스템 명세 기반의 임베디드 장치 모의 모델링 및 시뮬레이션 시스템 및 그의 방법 |
-
2016
- 2016-04-07 CN CN201610220029.XA patent/CN105912776B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR20060071302A (ko) * | 2004-12-21 | 2006-06-26 | 한국전자통신연구원 | 사용자 인터페이스 설계 및 평가 시스템 및 손 상호 작용기반의 사용자 인터페이스 설계 및 평가 시스템 |
CN102520726A (zh) * | 2011-12-19 | 2012-06-27 | 南京航空航天大学 | 大攻角飞行状态下的大气攻角及侧滑角估计方法 |
CN102880732A (zh) * | 2011-12-28 | 2013-01-16 | 南京康尼机电股份有限公司 | 一种轨道交通车辆门系统动力学联合仿真分析方法 |
CN103921947A (zh) * | 2014-03-27 | 2014-07-16 | 中国科学院长春光学精密机械与物理研究所 | 跳伞模拟训练半实物仿真系统及其操作方法 |
KR101510058B1 (ko) * | 2014-06-30 | 2015-04-08 | 성결대학교 산학협력단 | 이산 사건 시스템 명세 기반의 임베디드 장치 모의 모델링 및 시뮬레이션 시스템 및 그의 방법 |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106446442A (zh) * | 2016-10-12 | 2017-02-22 | 北京理工大学 | 一种火星伞舱组合体的降落伞展开过程稳定性评估方法 |
CN106446442B (zh) * | 2016-10-12 | 2019-12-13 | 北京理工大学 | 一种火星伞舱组合体的降落伞展开过程稳定性评估方法 |
CN106767781A (zh) * | 2016-11-29 | 2017-05-31 | 中国地质大学(武汉) | 跌落测试预埋传感器的六自由度运动轨迹数据处理方法 |
CN108268427A (zh) * | 2018-01-10 | 2018-07-10 | 中国地质大学(武汉) | 一种任意球进球概率分析方法、设备及存储设备 |
CN114202633A (zh) * | 2021-12-06 | 2022-03-18 | 北京理工大学 | 一种耦合空气动力学的跳台滑雪肌骨动力学模拟方法 |
CN114202633B (zh) * | 2021-12-06 | 2024-05-24 | 北京理工大学 | 一种耦合空气动力学的跳台滑雪肌骨动力学模拟方法 |
CN114638444A (zh) * | 2022-05-19 | 2022-06-17 | 珠海翔翼航空技术有限公司 | 针对突发故障下的跳伞时间预测方法、系统、设备 |
CN114638444B (zh) * | 2022-05-19 | 2022-09-02 | 珠海翔翼航空技术有限公司 | 针对突发故障下的跳伞时间预测方法、系统、设备 |
Also Published As
Publication number | Publication date |
---|---|
CN105912776B (zh) | 2019-03-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105912776B (zh) | 一种基于多刚体模型的跳伞员坠落运动仿真方法 | |
CN102073755B (zh) | 近空间高超声速飞行器运动控制仿真方法 | |
CN105116910B (zh) | 一种对地面点凝视成像的卫星姿态控制方法 | |
CN106446442B (zh) | 一种火星伞舱组合体的降落伞展开过程稳定性评估方法 | |
CN107145071B (zh) | 一种空间系绳系统的拉格朗日动力模型及控制器 | |
CN103942085B (zh) | 一种舰载直升机飞行模拟器的着舰仿真方法 | |
CN109319171B (zh) | 一种空间碎片横向角速度抑制和自旋方向控制方法 | |
CN112906212B (zh) | 一种基于绝对节点坐标法的裸电动力绳系统建模方法 | |
CN116027905A (zh) | 基于惯性传感器的双人皮划艇上肢动作捕捉方法 | |
CN106249616A (zh) | 一种在轨服务机械臂动力学建模方法和系统 | |
HANAVAN JR | A personalized mathematical model of the human body. | |
Li et al. | Modeling and control of tethered undersea kites | |
CN114756955A (zh) | 一种折叠翼飞行器分离模拟方法 | |
CN115270313A (zh) | 一种伞-机组合体的建模方法、装置、服务器和存储介质 | |
CN108733858B (zh) | 应用于高空飞行器系统的建模方法和装置 | |
Bingbing et al. | A new modeling scheme for powered parafoil unmanned aerial vehicle platforms: Theory and experiments | |
Calogero et al. | A dynamic spar numerical model for passive shape change | |
CN105320807B (zh) | 一种大规模空降空投着陆点预评估方法 | |
CN106363646A (zh) | 基于视觉伺服控制的多旋翼和机载机械臂联合位姿控制法 | |
Azouz et al. | Dynamic analysis of airships with small deformations | |
Xiu et al. | Dynamic modelling and simulation of a deployable quadrotor | |
CN110298083B (zh) | 一种系绳式洋流发电机运动建模方法 | |
Di Capua et al. | Body Pose Measurement System: System Validation and Range of Motion/Kinematic Analysis of Three Pressure Suits | |
Calogero et al. | A Dynamic Spar Numerical Model for Passive Shape Change | |
Vance et al. | A universal technique for grappling non-cooperative objects using planar, multi-jointed arms |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190308 Termination date: 20200407 |
|
CF01 | Termination of patent right due to non-payment of annual fee |