一种非光滑多体系统动力学灵敏度分析方法

1.本发明属于多体系统动力学优化技术领域，涉及一种非光滑多体系统动力学的灵敏度分析方法。


背景技术：

2.近年来，随着人们对机械系统的研究越来越趋于轻量化、复杂化和快速化，不可避免地会涉及多体系统的接触碰撞和摩擦问题，例如智能机器人系统与物体之间的接触运动、航天可展空间结构组件之间的接触运动。由于接触状态和接触边界都不能事先确定，使得对涉及接触碰撞和摩擦的多体系统动力学优化问题的研究仍面临挑战。
3.灵敏度分析是解决多体系统动力学优化问题的关键环节。多体系统动力学灵敏度分析的主要方法为有限差分法、直接微分法和伴随变量法。其中，有限差分法是近似方法，直接微分法和伴随变量法是解析方法。由于多体系统动力学的状态变量(位置、速度、加速度等)与时间有关，动力学灵敏度分析包括两部分：状态变量关于设计变量的导数和目标函数关于设计变量的导数，其中前者称为状态灵敏度，后者称为设计灵敏度。直接微分法是通过直接求导得到状态灵敏度和设计灵敏度；而伴随变量法是通过引进伴随变量消除状态灵敏度，进而得到系统的设计灵敏度，近年来由于其计算速度相对较快而得到了广泛的应用。对于不考虑接触碰撞和摩擦的光滑多体系统，许多学者采用两种解析方法做了大量的研究工作。然而，对于考虑接触碰撞和摩擦的非光滑多体系统，由于互补不等式方程和摩擦约束的复杂性，使得动力学灵敏度分析方法的研究十分困难，很少有研究报道。最近，corner等基于速度状态变量和灵敏度变量的跳跃条件，分别采用直接微分法和伴随变量法研究了混合系统的动力学灵敏度，但摩擦问题没有被细节的涉及。


技术实现要素：

4.本发明提出一种非光滑多体系统动力学灵敏度分析方法。该方法采用伴随变量法，建立非光滑多体系统动力学灵敏度分析框架，以解决非光滑多体系统动力学灵敏度分析问题，目的在于提供一套完整、统一的非光滑多体系统动力学灵敏度分析的新策略，以解决涉及接触碰撞和摩擦行为的多体系统的灵敏度分析问题。
5.为了达到上述目的，本发明采用的技术方案如下：
6.一种非光滑多体系统动力学灵敏度分析方法，首先，根据待解决的优化问题，给出非光滑多体系统动力学优化问题的目标函数、设计变量和约束条件；其次，基于非光滑接触方法，建立非光滑多体系统的动力学模型；再次，采用伴随变量法，建立非光滑多体系统的解析灵敏度分析方法或半解析灵敏度分析方法的计算公式；最后，通过数值积分方法，求解非光滑多体系统动力学灵敏度。包括以下步骤：
7.第一步，建立非光滑多体系统动力学优化问题的目标函数
8.由于非光滑多体系统动力学问题与时间t有关，其优化问题的目标函数一般表示为：
9.10.式中，b＝[b1,b2,
…
,b
p
]
t
为设计变量，其中，下标p为设计变量的数目；q,和分别为位置矢量、速度矢量、加速度矢量和拉格朗日乘子矢量，其中，n为系统自由度的数目，m为系统约束的数目；上标0和f分别表示相应变量在系统初始和终止时刻的值；g0和gf分别为与系统初始和终止状态有关的项；h为与系统运动过程有关的项；t0和tf为初始和终止时刻，由下式得到：
[0011][0012]
式中，ω0和ωf为系统的初始和终止时刻条件。
[0013]
系统的初始状态需要满足以下附加相容条件：
[0014][0015][0016]
式中，和分别为关于位置和速度的初始附加相容条件。
[0017]
第二步，建立非光滑多体系统的动力学模型
[0018]
考虑多体系统动力学中的接触碰撞和摩擦现象，基于非光滑接触方法，本发明在离散时间内将非光滑接触分解为法向碰撞和切向摩擦两部分，分别建立非光滑多体系统的光滑部分、非光滑接触碰撞部分和非光滑接触摩擦部分的动力学模型。
[0019]
所述的光滑部分的动力学模型表示为：
[0020][0021]
φ(q,b,t)＝0
ꢀꢀ
(6)
[0022]
式中，为系统的质量矩阵；为系统的位置约束；为位置约束的jacobian矩阵；为系统的广义力矢量。
[0023]
所述的非光滑接触碰撞部分的动力学模型表示为：
[0024][0025][0026][0027]
式中，为由接触碰撞引起的速度跃变；和分别为接触碰撞的双边约束和单边约束的拉格朗日乘子，其中，mc为激活的单边约束的数目；λ
i,n
为拉格朗日乘子λn的第i个分量；为接触碰撞的方向矩阵；d
i,n
＝[0,...,0,-n
t
,0,...,0,n
t
,0,...,0]为方向矩阵dn的第i个分量，其中，n＝[n
x
,ny,nz]
t
为法向量，n
x
,ny和nz分别为法向量n在三个坐标轴上的方向分量；为发生碰撞前的速度矢量；c为碰撞系数；ε为一个小值。
[0028]
所述的非光滑接触摩擦部分的动力学模型表示为：
[0029][0030][0031][0032]
式中，为由接触摩擦引起的速度跃变；为接触摩擦的双边约束的拉
格朗日乘子；为切向相对速度；μ为摩擦系数；γ1,分别为摩擦力在两个正交基向量上的分量；γ
i,1
和γ
i,2
分别为与第i个接触对应的摩擦力在两个正交基向量上的分量；分别为与第i个接触对应的摩擦力在两个正交基向量上的分量；分别由定义摩擦方向的两个正交基向量组成；和和分别为与第i个接触对应的两个正交基向量方向，其中，u1和u2分别表示为：
[0033][0034][0035]
需要注意的是，式(12)和式(8)均表示为速度级别的双边约束，表示形式是相同的，但分别关于接触碰撞和摩擦的计算结果并不相同。
[0036]
第三步，建立非光滑多体系统动力学灵敏度计算公式
[0037]
采用伴随变量法，给出完整统一的非光滑多体系统动力学灵敏度分析的解析框架，并在此基础上，结合有限差分法，给出非光滑多体系统动力学半解析灵敏度分析方法。
[0038]
(a)非光滑多体系统动力学的解析灵敏度计算公式
[0039]
一般地，非光滑多体系统动力学的设计灵敏度计算公式表示为：
[0040][0041]
式中，下标表示关于相关变量的导数；qb,和λb是状态灵敏度，为未知量；和分别为初始和终止时刻关于设计变量的导数，由下式得到：
[0042][0043]
式中，变量上方的点表示该变量关于时间的导数。
[0044]
采用伴随变量法，引进伴随变量采用伴随变量法，引进伴随变量和消除状态灵敏度，得到下列光滑部分、非光滑接触碰撞和摩擦部分、以及与非光滑多体系统初始和终止状态有关的伴随变量方程：
[0045][0046][0047][0048]
φqκ＝0
ꢀꢀ
(20)
[0049][0050][0051]
φqτ＝0
ꢀꢀ
(23)
[0052][0053][0054][0055][0056][0057][0058]
其中
[0059][0060]
伴随变量可由上述一系列伴随变量方程求解得到，将它们代入消除状态灵敏度后的非光滑多体系统动力学的设计灵敏度计算公式：
[0061][0062]
即可得到非光滑多体系统动力学灵敏度。
[0063]
(b)非光滑多体系统动力学的半解析灵敏度计算公式
[0064]
结合有限差分法，将非光滑多体系统动力学解析灵敏度计算公式中关于设计变量的导数项替换为有限差分，得到半解析灵敏度计算公式：
[0065][0066]
式中，为光滑部分的有限差分项；为非光滑接触碰撞和摩擦部分的有限差分项；和为与非光滑多体系统初始和终止状态有关的有限差分项，分别表示为：
[0067][0068][0069][0070][0071]
[0072][0073][0074][0075][0076][0077][0078][0079][0080]
式中，δb表示设计变量的扰动。
[0081]
本发明的有益效果为：
[0082]
(1)本发明同时考虑多体系统中广泛存在的接触碰撞和接触摩擦现象，采用伴随变量法给出了一个完整、统一的非光滑多体系统动力学灵敏度分析的解析框架，相比于直接微分法，可以准确高效的实现非光滑多体系统动力学灵敏度分析，极大地提高了涉及接触碰撞和摩擦行为的非光滑多体系统动力学优化问题的求解效率。
[0083]
(2)本发明在给出的非光滑多体系统动力学灵敏度分析的解析框架的基础上，进一步给出了非光滑多体系统动力学的半解析灵敏度分析方法的计算公式。该方法比有限差分法具有更好的计算精度和更高的计算效率，避免了解析公式中关于设计变量导数项的公式推导和计算机程序编码工作，对设计变量类型的处理更为灵活，极大地便利了结构复杂、设计变量种类较多的非光滑多体系统动力学灵敏度的计算。
附图说明
[0084]
图1是本发明实施流程图。
[0085]
图2是五连杆系统动力学模型。
[0086]
图3是五连杆系统的节点c在2s时间范围内的仿真结果；其中，图(a)为y轴方向上的位置；图(b)为y轴方向上的速度。
[0087]
图4是有限差分法、半解析方法和解析方法的灵敏度计算结果随设计变量扰动值变化曲线；其中，图(a)为目标函数关于第一根弹簧原长的灵敏度；图(b)为目标函数关于第二根弹簧原长的灵敏度。
[0088]
图5是半解析方法和解析方法的灵敏度计算时间历程曲线；其中，图(a)为目标函数关于第一根弹簧原长的灵敏度；图(b)为目标函数关于第二根弹簧原长的灵敏度。
具体实施方式
[0089]
以下结合具体实施例对本发明做进一步说明。
[0090]
结合图1，本发明具体按以下步骤实现：
[0091]
第一步，根据待解决的优化问题，给出非光滑多体系统动力学优化问题的目标函数、设计变量和约束条件。
[0092]
第二步，基于非光滑接触方法，建立非光滑多体系统的动力学模型，具体步骤为：
[0093]
(1)给出该非光滑多体系统的质量矩阵、位置约束、位置约束的jacobian矩阵和广义力矢量；
[0094]
(2)给出该非光滑多体系统是否产生接触碰撞和摩擦的判断条件；
[0095]
(3)给出该非光滑多体系统接触碰撞的方向矩阵和碰撞系数；
[0096]
(4)给出该非光滑多体系统接触摩擦的两个正交基向量方向和摩擦系数；
[0097]
(5)在步骤(1)-(4)的基础上，建立非光滑多体系统的动力学方程。
[0098]
第三步，采用伴随变量法，建立非光滑多体系统动力学灵敏度计算公式，具体步骤为：
[0099]
(1)建立一般形式的非光滑多体系统动力学的设计灵敏度计算公式ψb；
[0100]
(2)采用伴随变量法，引入伴随变量μ,ν,κ,η,χ,τ,ζ,σi,ρi,ξi,π和β(i＝0,f)消除状态灵敏度qb,和λb，得到一系列伴随变量方程；
[0101]
(3)建立非光滑多体系统的解析灵敏度分析方法或半解析灵敏度分析方法的设计灵敏度计算公式，步骤分别为：
[0102]
(a)解析灵敏度分析方法
[0103]
通过步骤(2)消掉步骤(1)公式中的状态灵敏度后，得到的ψ
′b即为非光滑多体系统的解析的设计灵敏度计算公式；
[0104]
(b)半解析灵敏度分析方法
[0105]
将解析设计灵敏度计算公式中关于设计变量的导数项替换为有限差分项，得到的ψ
′b′
即为非光滑多体系统的半解析的设计灵敏度计算公式。
[0106]
第四步，通过数值积分方法，求解非光滑多体系统动力学方程和伴随变量方程，将计算结果代入解析灵敏度分析方法或半解析灵敏度分析方法的设计灵敏度计算公式中，即可得到非光滑多体系统动力学灵敏度。
[0107]
仿真实例：利用本发明方法，针对考虑接触碰撞和摩擦现象的五连杆系统算例，进行数值仿真验证。
[0108]
图2为五连杆系统的动力学模型，采用绝对坐标建模，各节点的坐标分别为a(-0.5,0)m,b(-1.5,-1)m,c(0,-2)m,d(1.5,-1)m和e(0.5,0)m。该系统的四根连杆均为刚体，其中固定节点a、e之间的固定连接为第五根连杆，连杆的长度分别为l
ab
＝l
de
＝1.4142m,l
bc
＝l
cd
＝1.8027m，质量分别为m
ab
＝m
de
＝1kg,m
bc
＝m
cd
＝1.5kg，两个弹簧的刚度系数为k1＝k2＝100n/m。
[0109]
五连杆系统的初始构型如图2所示，系统在重力和两根弹簧的弹性力作用下运动。此外，y＝-2.35m处为地面，当系统最下面的节点c与地面发生接触时，将会受到碰撞力的作用而向上反弹；由于节点c有水平方向的速度，将会受到与地面之间的摩擦力的影响，摩擦系数为0.05。图3给出了节点c在2s时间范围内y轴方向上的位置和速度的仿真结果，从中可
以看出，当节点c与地面接触时，速度变量会发生突变，这是由接触碰撞和摩擦引起的速度变量跳跃，产生了明显的非光滑现象。
[0110]
非光滑多体系统动力学优化问题的设计变量为b＝[l
k1
,l
k2
]
t
，分别表示为两个弹簧的原长，取设计变量的初值为l
k1
＝2.2360m,l
k2
＝2.0615m。当五连杆系统运动的起始和终止时刻分别为t0＝0,tf＝2时，给定优化问题的目标函数为其中，(x,y)为节点c的位置坐标，该优化问题被定义以最小化五连杆系统运动过程中节点c与坐标原点之间距离。接下来，分别采用有限差分法、本发明的解析灵敏度分析方法和半解析灵敏度分析方法，进行考虑接触碰撞和摩擦现象的五连杆系统的灵敏度分析计算。
[0111]
图4为取不同的设计变量扰动值δ＝1,2,
…
,10时，分别采用有限差分法(fdm)、半解析灵敏度分析方法(sam)和解析灵敏度分析方法(am)计算五连杆系统的目标函数关于两个设计变量的灵敏度计算结果，其中，δ＝-log(δb/b)。从中可以得到：有限差分法的计算结果受设计变量扰动值影响明显，当δ较小时，有限差分法的计算结果无法保持稳定；而半解析灵敏度分析方法基本不受设计变量扰动值的影响，其计算结果始终与解析灵敏度分析方法一致，说明了本发明提出的非光滑多体系统的半解析灵敏度分析方法的数值稳定性。
[0112]
表1当δ＝7时，不同方法的灵敏度计算结果和计算时间
[0113][0114]
具体地，表1为设计变量扰动值δ＝7时，采用有限差分法(fdm)、本发明提出的半解析灵敏度分析方法(sam)和解析灵敏度分析方法(am)，得到的五连杆系统的目标函数关于设计变量的灵敏度计算结果和计算时间。从中可以得到：三种灵敏度分析方法的灵敏度计算结果基本一致，证明了本发明提出的非光滑多体系统动力学灵敏度分析方法的正确性；考虑计算时间，有限差分法的计算时间明显长于解析方法和半解析方法，其中，半解析方法的计算时间比解析方法长，但相差很小，证明了本发明提出的非光滑多体系统动力学灵敏度分析方法具有较高的计算效率。
[0115]
图5为设计变量扰动值δ＝7时，分别采用半解析灵敏度分析方法(sam)和解析灵敏度分析方法(am)计算五连杆系统在2s时间范围内，目标函数关于两个设计变量的灵敏度计算结果随时间变化曲线。从中可以得到：半解析灵敏度分析方法的计算结果与解析灵敏度分析方法几乎完全一致，证明了本发明提出的非光滑多体系统的半解析灵敏度分析方法具有较高的计算精度。
[0116]
以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制。应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些均属于本发明的保护范围。