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

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

Info

Publication number
CN114186452A
CN114186452A CN202111394494.2A CN202111394494A CN114186452A CN 114186452 A CN114186452 A CN 114186452A CN 202111394494 A CN202111394494 A CN 202111394494A CN 114186452 A CN114186452 A CN 114186452A
Authority
CN
China
Prior art keywords
smooth
body system
sensitivity
contact
variable
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
Application number
CN202111394494.2A
Other languages
English (en)
Other versions
CN114186452B (zh
Inventor
彭海军
张孟茹
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Dalian University of Technology
Original Assignee
Dalian University of Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Dalian University of Technology filed Critical Dalian University of Technology
Priority to CN202111394494.2A priority Critical patent/CN114186452B/zh
Publication of CN114186452A publication Critical patent/CN114186452A/zh
Application granted granted Critical
Publication of CN114186452B publication Critical patent/CN114186452B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)

Abstract

本发明属于动力学系统优化技术领域，提供一种非光滑多体系统动力学灵敏度分析方法。首先，给出非光滑多体系统动力学优化问题的目标函数、设计变量和约束条件；其次，基于非光滑接触方法，建立非光滑多体系统的动力学模型；再次，采用伴随变量法，建立非光滑多体系统的解析灵敏度分析方法或半解析灵敏度分析方法的计算公式；最后，通过数值积分方法，求解非光滑多体系统动力学灵敏度。本发明采用伴随变量法，建立非光滑多体系统动力学灵敏度分析框架，以解决非光滑多体系统动力学灵敏度分析问题，目的在于提供一套完整、统一的非光滑多体系统动力学灵敏度分析的新策略，以解决涉及接触碰撞和摩擦行为的多体系统的灵敏度分析问题。

Description

一种非光滑多体系统动力学灵敏度分析方法
技术领域
本发明属于多体系统动力学优化技术领域，涉及一种非光滑多体系统动力学的灵敏度分析方法。
背景技术
近年来，随着人们对机械系统的研究越来越趋于轻量化、复杂化和快速化，不可避免地会涉及多体系统的接触碰撞和摩擦问题，例如智能机器人系统与物体之间的接触运动、航天可展空间结构组件之间的接触运动。由于接触状态和接触边界都不能事先确定，使得对涉及接触碰撞和摩擦的多体系统动力学优化问题的研究仍面临挑战。
灵敏度分析是解决多体系统动力学优化问题的关键环节。多体系统动力学灵敏度分析的主要方法为有限差分法、直接微分法和伴随变量法。其中，有限差分法是近似方法，直接微分法和伴随变量法是解析方法。由于多体系统动力学的状态变量(位置、速度、加速度等)与时间有关，动力学灵敏度分析包括两部分：状态变量关于设计变量的导数和目标函数关于设计变量的导数，其中前者称为状态灵敏度，后者称为设计灵敏度。直接微分法是通过直接求导得到状态灵敏度和设计灵敏度；而伴随变量法是通过引进伴随变量消除状态灵敏度，进而得到系统的设计灵敏度，近年来由于其计算速度相对较快而得到了广泛的应用。对于不考虑接触碰撞和摩擦的光滑多体系统，许多学者采用两种解析方法做了大量的研究工作。然而，对于考虑接触碰撞和摩擦的非光滑多体系统，由于互补不等式方程和摩擦约束的复杂性，使得动力学灵敏度分析方法的研究十分困难，很少有研究报道。最近，Corner等基于速度状态变量和灵敏度变量的跳跃条件，分别采用直接微分法和伴随变量法研究了混合系统的动力学灵敏度，但摩擦问题没有被细节的涉及。
发明内容
本发明提出一种非光滑多体系统动力学灵敏度分析方法。该方法采用伴随变量法，建立非光滑多体系统动力学灵敏度分析框架，以解决非光滑多体系统动力学灵敏度分析问题，目的在于提供一套完整、统一的非光滑多体系统动力学灵敏度分析的新策略，以解决涉及接触碰撞和摩擦行为的多体系统的灵敏度分析问题。
为了达到上述目的，本发明采用的技术方案如下：
一种非光滑多体系统动力学灵敏度分析方法，首先，根据待解决的优化问题，给出非光滑多体系统动力学优化问题的目标函数、设计变量和约束条件；其次，基于非光滑接触方法，建立非光滑多体系统的动力学模型；再次，采用伴随变量法，建立非光滑多体系统的解析灵敏度分析方法或半解析灵敏度分析方法的计算公式；最后，通过数值积分方法，求解非光滑多体系统动力学灵敏度。包括以下步骤：
第一步，建立非光滑多体系统动力学优化问题的目标函数
由于非光滑多体系统动力学问题与时间t有关，其优化问题的目标函数一般表示为：
Figure BDA0003369501600000021
式中，b＝[b1,b2,…,bp]T为设计变量，其中，下标p为设计变量的数目；q,
Figure BDA0003369501600000022
和
Figure BDA0003369501600000023
分别为位置矢量、速度矢量、加速度矢量和拉格朗日乘子矢量，其中，n为系统自由度的数目，m为系统约束的数目；上标0和f分别表示相应变量在系统初始和终止时刻的值；G0和Gf分别为与系统初始和终止状态有关的项；H为与系统运动过程有关的项；t0和tf为初始和终止时刻，由下式得到：
Figure BDA0003369501600000024
式中，Ω0和Ωf为系统的初始和终止时刻条件。
系统的初始状态需要满足以下附加相容条件：
Figure BDA0003369501600000025
Figure BDA0003369501600000026
式中，
Figure BDA0003369501600000027
和
Figure BDA0003369501600000028
分别为关于位置和速度的初始附加相容条件。
第二步，建立非光滑多体系统的动力学模型
考虑多体系统动力学中的接触碰撞和摩擦现象，基于非光滑接触方法，本发明在离散时间内将非光滑接触分解为法向碰撞和切向摩擦两部分，分别建立非光滑多体系统的光滑部分、非光滑接触碰撞部分和非光滑接触摩擦部分的动力学模型。
所述的光滑部分的动力学模型表示为：
Figure BDA0003369501600000029
Φ(q,b,t)＝0 (6)
式中，
Figure BDA00033695016000000210
为系统的质量矩阵；
Figure BDA00033695016000000211
为系统的位置约束；
Figure BDA00033695016000000212
为位置约束的Jacobian矩阵；
Figure BDA00033695016000000213
为系统的广义力矢量。
所述的非光滑接触碰撞部分的动力学模型表示为：
Figure BDA00033695016000000214
Figure BDA00033695016000000215
Figure BDA00033695016000000216
式中，
Figure BDA00033695016000000217
为由接触碰撞引起的速度跃变；
Figure BDA00033695016000000218
和
Figure BDA00033695016000000219
分别为接触碰撞的双边约束和单边约束的拉格朗日乘子，其中，mc为激活的单边约束的数目；Λi,n为拉格朗日乘子Λn的第i个分量；
Figure BDA0003369501600000031
为接触碰撞的方向矩阵；Di,n＝[0,...,0,-nT,0,...,0,nT,0,...,0]为方向矩阵Dn的第i个分量，其中，n＝[nx,ny,nz]T为法向量，nx,ny和nz分别为法向量n在三个坐标轴上的方向分量；
Figure BDA0003369501600000032
为发生碰撞前的速度矢量；c为碰撞系数；ε为一个小值。
所述的非光滑接触摩擦部分的动力学模型表示为：
Figure BDA0003369501600000033
Figure BDA0003369501600000034
Figure BDA0003369501600000035
式中，
Figure BDA0003369501600000036
为由接触摩擦引起的速度跃变；
Figure BDA0003369501600000037
为接触摩擦的双边约束的拉格朗日乘子；
Figure BDA0003369501600000038
为切向相对速度；μ为摩擦系数；γ1,
Figure BDA0003369501600000039
分别为摩擦力在两个正交基向量上的分量；γi,1和γi,2分别为与第i个接触对应的摩擦力在两个正交基向量上的分量；
Figure BDA00033695016000000310
Figure BDA00033695016000000311
分别由定义摩擦方向的两个正交基向量组成；
Figure BDA00033695016000000312
和
Figure BDA00033695016000000313
Figure BDA00033695016000000314
分别为与第i个接触对应的两个正交基向量方向，其中，u1和u2分别表示为：
Figure BDA00033695016000000315
Figure BDA00033695016000000316
需要注意的是，式(12)和式(8)均表示为速度级别的双边约束，表示形式是相同的，但分别关于接触碰撞和摩擦的计算结果并不相同。
第三步，建立非光滑多体系统动力学灵敏度计算公式
采用伴随变量法，给出完整统一的非光滑多体系统动力学灵敏度分析的解析框架，并在此基础上，结合有限差分法，给出非光滑多体系统动力学半解析灵敏度分析方法。
(a)非光滑多体系统动力学的解析灵敏度计算公式
一般地，非光滑多体系统动力学的设计灵敏度计算公式表示为：
Figure BDA00033695016000000317
式中，下标表示关于相关变量的导数；qb,
Figure BDA0003369501600000041
和λb是状态灵敏度，为未知量；
Figure BDA0003369501600000042
和
Figure BDA0003369501600000043
分别为初始和终止时刻关于设计变量的导数，由下式得到：
Figure BDA0003369501600000044
式中，变量上方的点表示该变量关于时间的导数。
采用伴随变量法，引进伴随变量
Figure BDA0003369501600000045
Figure BDA0003369501600000046
和
Figure BDA0003369501600000047
消除状态灵敏度，得到下列光滑部分、非光滑接触碰撞和摩擦部分、以及与非光滑多体系统初始和终止状态有关的伴随变量方程：
Figure BDA0003369501600000048
Figure BDA0003369501600000049
Figure BDA00033695016000000410
Φqκ＝0 (20)
Figure BDA00033695016000000411
Figure BDA00033695016000000412
Φqτ＝0 (23)
Figure BDA00033695016000000413
Figure BDA00033695016000000414
Figure BDA00033695016000000415
Figure BDA00033695016000000416
Figure BDA00033695016000000417
Figure BDA00033695016000000418
其中
Figure BDA00033695016000000419
伴随变量可由上述一系列伴随变量方程求解得到，将它们代入消除状态灵敏度后的非光滑多体系统动力学的设计灵敏度计算公式：
Figure BDA00033695016000000420
即可得到非光滑多体系统动力学灵敏度。
(b)非光滑多体系统动力学的半解析灵敏度计算公式
结合有限差分法，将非光滑多体系统动力学解析灵敏度计算公式中关于设计变量的导数项替换为有限差分，得到半解析灵敏度计算公式：
Figure BDA0003369501600000051
式中，
Figure BDA0003369501600000052
为光滑部分的有限差分项；
Figure BDA0003369501600000053
为非光滑接触碰撞和摩擦部分的有限差分项；
Figure BDA0003369501600000054
和
Figure BDA0003369501600000055
为与非光滑多体系统初始和终止状态有关的有限差分项，分别表示为：
Figure BDA0003369501600000056
Figure BDA0003369501600000057
Figure BDA0003369501600000058
Figure BDA0003369501600000059
Figure BDA00033695016000000510
Figure BDA00033695016000000511
Figure BDA00033695016000000512
Figure BDA00033695016000000513
Figure BDA00033695016000000514
Figure BDA00033695016000000515
Figure BDA00033695016000000516
Figure BDA00033695016000000517
Figure BDA00033695016000000518
式中，Δb表示设计变量的扰动。
本发明的有益效果为：
(1)本发明同时考虑多体系统中广泛存在的接触碰撞和接触摩擦现象，采用伴随变量法给出了一个完整、统一的非光滑多体系统动力学灵敏度分析的解析框架，相比于直接微分法，可以准确高效的实现非光滑多体系统动力学灵敏度分析，极大地提高了涉及接触碰撞和摩擦行为的非光滑多体系统动力学优化问题的求解效率。
(2)本发明在给出的非光滑多体系统动力学灵敏度分析的解析框架的基础上，进一步给出了非光滑多体系统动力学的半解析灵敏度分析方法的计算公式。该方法比有限差分法具有更好的计算精度和更高的计算效率，避免了解析公式中关于设计变量导数项的公式推导和计算机程序编码工作，对设计变量类型的处理更为灵活，极大地便利了结构复杂、设计变量种类较多的非光滑多体系统动力学灵敏度的计算。
附图说明
图1是本发明实施流程图。
图2是五连杆系统动力学模型。
图3是五连杆系统的节点C在2s时间范围内的仿真结果；其中，图(a)为y轴方向上的位置；图(b)为y轴方向上的速度。
图4是有限差分法、半解析方法和解析方法的灵敏度计算结果随设计变量扰动值变化曲线；其中，图(a)为目标函数关于第一根弹簧原长的灵敏度；图(b)为目标函数关于第二根弹簧原长的灵敏度。
图5是半解析方法和解析方法的灵敏度计算时间历程曲线；其中，图(a)为目标函数关于第一根弹簧原长的灵敏度；图(b)为目标函数关于第二根弹簧原长的灵敏度。
具体实施方式
以下结合具体实施例对本发明做进一步说明。
结合图1，本发明具体按以下步骤实现：
第一步，根据待解决的优化问题，给出非光滑多体系统动力学优化问题的目标函数、设计变量和约束条件。
第二步，基于非光滑接触方法，建立非光滑多体系统的动力学模型，具体步骤为：
(1)给出该非光滑多体系统的质量矩阵、位置约束、位置约束的Jacobian矩阵和广义力矢量；
(2)给出该非光滑多体系统是否产生接触碰撞和摩擦的判断条件；
(3)给出该非光滑多体系统接触碰撞的方向矩阵和碰撞系数；
(4)给出该非光滑多体系统接触摩擦的两个正交基向量方向和摩擦系数；
(5)在步骤(1)-(4)的基础上，建立非光滑多体系统的动力学方程。
第三步，采用伴随变量法，建立非光滑多体系统动力学灵敏度计算公式，具体步骤为：
(1)建立一般形式的非光滑多体系统动力学的设计灵敏度计算公式Ψb；
(2)采用伴随变量法，引入伴随变量μ,ν,κ,η,χ,τ,ζ,σi,ρi,ξi,π和β(i＝0,f)消除状态灵敏度qb,
Figure BDA0003369501600000071
和λb，得到一系列伴随变量方程；
(3)建立非光滑多体系统的解析灵敏度分析方法或半解析灵敏度分析方法的设计灵敏度计算公式，步骤分别为：
(a)解析灵敏度分析方法
通过步骤(2)消掉步骤(1)公式中的状态灵敏度后，得到的Ψ′b即为非光滑多体系统的解析的设计灵敏度计算公式；
(b)半解析灵敏度分析方法
将解析设计灵敏度计算公式中关于设计变量的导数项替换为有限差分项，得到的Ψ′b′即为非光滑多体系统的半解析的设计灵敏度计算公式。
第四步，通过数值积分方法，求解非光滑多体系统动力学方程和伴随变量方程，将计算结果代入解析灵敏度分析方法或半解析灵敏度分析方法的设计灵敏度计算公式中，即可得到非光滑多体系统动力学灵敏度。
仿真实例：利用本发明方法，针对考虑接触碰撞和摩擦现象的五连杆系统算例，进行数值仿真验证。
图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之间的固定连接为第五根连杆，连杆的长度分别为LAB＝LDE＝1.4142m,LBC＝LCD＝1.8027m，质量分别为mAB＝mDE＝1kg,mBC＝mCD＝1.5kg，两个弹簧的刚度系数为k1＝k2＝100N/m。
五连杆系统的初始构型如图2所示，系统在重力和两根弹簧的弹性力作用下运动。此外，y＝-2.35m处为地面，当系统最下面的节点C与地面发生接触时，将会受到碰撞力的作用而向上反弹；由于节点C有水平方向的速度，将会受到与地面之间的摩擦力的影响，摩擦系数为0.05。图3给出了节点C在2s时间范围内y轴方向上的位置和速度的仿真结果，从中可以看出，当节点C与地面接触时，速度变量会发生突变，这是由接触碰撞和摩擦引起的速度变量跳跃，产生了明显的非光滑现象。
非光滑多体系统动力学优化问题的设计变量为b＝[Lk1,Lk2]T，分别表示为两个弹簧的原长，取设计变量的初值为Lk1＝2.2360m,Lk2＝2.0615m。当五连杆系统运动的起始和终止时刻分别为t0＝0,tf＝2时，给定优化问题的目标函数为
Figure BDA0003369501600000072
其中，(x,y)为节点C的位置坐标，该优化问题被定义以最小化五连杆系统运动过程中节点C与坐标原点之间距离。接下来，分别采用有限差分法、本发明的解析灵敏度分析方法和半解析灵敏度分析方法，进行考虑接触碰撞和摩擦现象的五连杆系统的灵敏度分析计算。
图4为取不同的设计变量扰动值δ＝1,2,…,10时，分别采用有限差分法(FDM)、半解析灵敏度分析方法(SAM)和解析灵敏度分析方法(AM)计算五连杆系统的目标函数关于两个设计变量的灵敏度计算结果，其中，δ＝-log(Δb/b)。从中可以得到：有限差分法的计算结果受设计变量扰动值影响明显，当δ较小时，有限差分法的计算结果无法保持稳定；而半解析灵敏度分析方法基本不受设计变量扰动值的影响，其计算结果始终与解析灵敏度分析方法一致，说明了本发明提出的非光滑多体系统的半解析灵敏度分析方法的数值稳定性。
表1当δ＝7时，不同方法的灵敏度计算结果和计算时间
Figure BDA0003369501600000081
具体地，表1为设计变量扰动值δ＝7时，采用有限差分法(FDM)、本发明提出的半解析灵敏度分析方法(SAM)和解析灵敏度分析方法(AM)，得到的五连杆系统的目标函数关于设计变量的灵敏度计算结果和计算时间。从中可以得到：三种灵敏度分析方法的灵敏度计算结果基本一致，证明了本发明提出的非光滑多体系统动力学灵敏度分析方法的正确性；考虑计算时间，有限差分法的计算时间明显长于解析方法和半解析方法，其中，半解析方法的计算时间比解析方法长，但相差很小，证明了本发明提出的非光滑多体系统动力学灵敏度分析方法具有较高的计算效率。
图5为设计变量扰动值δ＝7时，分别采用半解析灵敏度分析方法(SAM)和解析灵敏度分析方法(AM)计算五连杆系统在2s时间范围内，目标函数关于两个设计变量的灵敏度计算结果随时间变化曲线。从中可以得到：半解析灵敏度分析方法的计算结果与解析灵敏度分析方法几乎完全一致，证明了本发明提出的非光滑多体系统的半解析灵敏度分析方法具有较高的计算精度。
以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制。应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些均属于本发明的保护范围。

Claims (7)

1.一种非光滑多体系统动力学灵敏度分析方法，其特征在于，首先，根据待解决的优化问题，给出非光滑多体系统动力学优化问题的目标函数、设计变量和约束条件；其次，基于非光滑接触方法，建立非光滑多体系统的动力学模型；再次，采用伴随变量法，建立非光滑多体系统的解析灵敏度分析方法或半解析灵敏度分析方法的计算公式；最后，通过数值积分方法，求解非光滑多体系统动力学灵敏度；包括下述步骤：
第一步、建立非光滑多体系统动力学优化问题的目标函数；
非光滑多体系统动力学问题与时间t有关，其优化问题的目标函数表示为：
Figure FDA0003369501590000011
式中，b＝[b1,b2,…,bp]T为设计变量，其中，下标p为设计变量的数目；
Figure FDA0003369501590000012
和
Figure FDA0003369501590000013
分别为位置矢量、速度矢量、加速度矢量和拉格朗日乘子矢量，其中，n为系统自由度的数目，m为系统约束的数目；上标0和f分别表示相应变量在系统初始和终止时刻的值；G0和Gf分别为与系统初始和终止状态有关的项；H为与系统中间过程有关的项；t0和tf为初始和终止时刻；
系统的初始状态需要满足以下附加相容条件：
Figure FDA0003369501590000014
Figure FDA0003369501590000015
式中，
Figure FDA0003369501590000016
和
Figure FDA0003369501590000017
分别为关于位置和速度的初始附加相容条件；
第二步、建立非光滑多体系统的动力学模型；
基于多体系统动力学中的接触碰撞和摩擦现象，基于非光滑接触方法，在离散时间内将非光滑接触分解为法向碰撞和切向摩擦两部分，分别建立非光滑多体系统的光滑部分、非光滑接触碰撞部分和非光滑接触摩擦部分的动力学模型；
第三步、建立非光滑多体系统动力学灵敏度计算公式；
采用伴随变量法给出非光滑多体系统动力学灵敏度分析的解析框架，并在此基础上，结合有限差分法，给出非光滑多体系统动力学半解析灵敏度分析方法；具体如下：
(a)非光滑多体系统动力学的解析灵敏度计算公式；
非光滑多体系统动力学的设计灵敏度计算公式表示为：
Figure FDA0003369501590000018
式中，下标表示关于相关变量的导数；qb,
Figure FDA0003369501590000019
和λb是状态灵敏度，为未知量；
Figure FDA00033695015900000110
和
Figure FDA00033695015900000111
分别为初始和终止时刻关于设计变量的导数，由下式得到：
Figure FDA0003369501590000021
式中，变量上方的点表示该变量关于时间的导数；
采用伴随变量法，引进伴随变量
Figure FDA0003369501590000022
Figure FDA0003369501590000023
和
Figure FDA0003369501590000024
消除状态灵敏度，得到光滑部分、非光滑接触碰撞和摩擦部分、以及与非光滑多体系统初始和终止状态有关的伴随变量方程，其中，光滑部分、与系统初始和终止状态有关的伴随变量方程表示为：
Figure FDA0003369501590000025
Figure FDA0003369501590000026
Figure FDA0003369501590000027
Figure FDA0003369501590000028
Figure FDA0003369501590000029
Figure FDA00033695015900000210
Figure FDA00033695015900000211
Figure FDA00033695015900000212
伴随变量由一系列伴随变量方程求解得到，将它们代入消除状态灵敏度后的非光滑多体系统动力学的设计灵敏度计算公式，最终得到非光滑多体系统动力学灵敏度；
(b)非光滑多体系统动力学的半解析灵敏度计算公式；
结合有限差分法，将非光滑多体系统动力学解析灵敏度计算公式中关于设计变量的导数项替换为有限差分，得到半解析灵敏度计算公式：
Figure FDA00033695015900000213
式中，
Figure FDA00033695015900000214
为光滑部分的有限差分项；
Figure FDA00033695015900000215
为非光滑接触碰撞和摩擦部分的有限差分项；
Figure FDA00033695015900000216
和
Figure FDA00033695015900000217
为与非光滑多体系统初始和终止状态有关的有限差分项，其中，光滑部分、与系统初始和终止状态有关的有限差分项表示为：
Figure FDA00033695015900000218
Figure FDA00033695015900000219
Figure FDA00033695015900000220
Figure FDA00033695015900000221
Figure FDA0003369501590000031
Figure FDA0003369501590000032
Figure FDA0003369501590000033
Figure FDA0003369501590000034
Figure FDA0003369501590000035
式中，Δb表示设计变量的扰动。
2.根据权利要求1所述的一种非光滑多体系统动力学灵敏度分析方法，其特征在于，在步骤二中，所述的光滑部分的动力学模型表示为：
Figure FDA0003369501590000036
Φ(q,b,t)＝0 (25)
式中，
Figure FDA0003369501590000037
为系统的质量矩阵；
Figure FDA0003369501590000038
为系统的位置约束；
Figure FDA0003369501590000039
为位置约束的Jacobian矩阵；
Figure FDA00033695015900000310
为系统的广义力矢量。
3.根据权利要求1所述的一种非光滑多体系统动力学灵敏度分析方法，其特征在于，在步骤二中，所述的非光滑接触碰撞部分的动力学模型表示为：
Figure FDA00033695015900000311
Figure FDA00033695015900000312
Figure FDA00033695015900000313
式中，
Figure FDA00033695015900000314
为由接触碰撞引起的速度跃变；
Figure FDA00033695015900000315
和
Figure FDA00033695015900000316
分别为接触碰撞的双边约束和单边约束的拉格朗日乘子，其中，mc为激活的单边约束的数目；Λi,n为拉格朗日乘子Λn的第i个分量；
Figure FDA00033695015900000317
为接触碰撞的方向矩阵；Di,n＝[0,...,0,-nT,0,…,0,nT,0,...,0]为方向矩阵Dn的第i个分量，其中，n＝[nx,ny,nz]T为法向量，nx,ny和nz分别为法向量n在三个坐标轴上的方向分量；
Figure FDA00033695015900000318
为发生碰撞前的速度矢量；c为碰撞系数；ε为一个小值。
4.根据权利要求1所述的一种非光滑多体系统动力学灵敏度分析方法，其特征在于，在步骤二中，所述的非光滑接触摩擦部分的动力学模型表示为：
Figure FDA00033695015900000319
Figure FDA00033695015900000320
Figure FDA0003369501590000041
式中，
Figure FDA0003369501590000042
为由接触摩擦引起的速度跃变；
Figure FDA0003369501590000043
为接触摩擦的双边约束的拉格朗日乘子；
Figure FDA0003369501590000044
为切向相对速度；μ为摩擦系数；
Figure FDA0003369501590000045
分别为摩擦力在两个正交基向量上的分量；γi,1和γi,2分别为与第i个接触对应的摩擦力在两个正交基向量上的分量；
Figure FDA0003369501590000046
分别由定义摩擦方向的两个正交基向量组成；
Figure FDA0003369501590000047
和
Figure FDA0003369501590000048
分别为与第i个接触对应的两个正交基向量方向。
5.根据权利要求1所述的一种非光滑多体系统动力学灵敏度分析方法，其特征在于，在步骤三(a)中，所述的非光滑接触碰撞和摩擦部分的伴随变量方程表示为：
Figure FDA0003369501590000049
Φqκ＝0 (33)
Figure FDA00033695015900000410
Figure FDA00033695015900000411
Φqτ＝0 (36)
式中，
Figure FDA00033695015900000412
6.根据权利要求1所述的一种非光滑多体系统动力学灵敏度分析方法，其特征在于，在步骤三(a)中，所述的消除状态灵敏度后的非光滑多体系统动力学的设计灵敏度计算公式表示为：
Figure FDA00033695015900000413
7.根据权利要求1所述的一种非光滑多体系统动力学灵敏度分析方法，其特征在于，在步骤三(b)中，所述的非光滑接触碰撞和摩擦部分的有限差分项表示为：
Figure FDA00033695015900000414
Figure FDA00033695015900000415
Figure FDA00033695015900000416
Figure FDA00033695015900000417
CN202111394494.2A 2021-11-23 2021-11-23 一种非光滑多体系统动力学灵敏度分析方法 Active CN114186452B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111394494.2A CN114186452B (zh) 2021-11-23 2021-11-23 一种非光滑多体系统动力学灵敏度分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111394494.2A CN114186452B (zh) 2021-11-23 2021-11-23 一种非光滑多体系统动力学灵敏度分析方法

Publications (2)

Publication Number Publication Date
CN114186452A true CN114186452A (zh) 2022-03-15
CN114186452B CN114186452B (zh) 2024-09-13

Family

ID=80541260

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111394494.2A Active CN114186452B (zh) 2021-11-23 2021-11-23 一种非光滑多体系统动力学灵敏度分析方法

Country Status (1)

Country Link
CN (1) CN114186452B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017000396A1 (zh) * 2015-06-30 2017-01-05 中国空间技术研究院 基于多体分析试验的桁架天线反射器展开动力学建模方法
CN108121865A (zh) * 2017-12-15 2018-06-05 北京航空航天大学 一种基于伴随变量法的含有多不确定参数的结构区间响应传播分析方法
CN110147571A (zh) * 2019-04-15 2019-08-20 浙江吉利控股集团有限公司 一种组件结构的拓扑优化方法及装置
CN111159636A (zh) * 2019-12-04 2020-05-15 大连理工大学 基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法
CN113806886A (zh) * 2021-09-13 2021-12-17 中山大学 一种基于增强响应灵敏度的多体系统参数识别方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017000396A1 (zh) * 2015-06-30 2017-01-05 中国空间技术研究院 基于多体分析试验的桁架天线反射器展开动力学建模方法
CN108121865A (zh) * 2017-12-15 2018-06-05 北京航空航天大学 一种基于伴随变量法的含有多不确定参数的结构区间响应传播分析方法
CN110147571A (zh) * 2019-04-15 2019-08-20 浙江吉利控股集团有限公司 一种组件结构的拓扑优化方法及装置
CN111159636A (zh) * 2019-12-04 2020-05-15 大连理工大学 基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法
CN113806886A (zh) * 2021-09-13 2021-12-17 中山大学 一种基于增强响应灵敏度的多体系统参数识别方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HAIJUN PENG 等: "Semi-Analytical Sensitivity Analysis for Multibody System Dynamics Described by Differential–Algebraic Equations", 《 AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS》, 31 January 2021 (2021-01-31) *
MENGRU ZHANG 等: "Semi-analytical sensitivity analysis approach for fully coupled optimization of flexible multibody systems", 《MECHANISM AND MACHINE THEORY》, 21 January 2021 (2021-01-21) *

Also Published As

Publication number Publication date
CN114186452B (zh) 2024-09-13

Similar Documents

Publication Publication Date Title
Wang et al. Optimal design of lightweight serial robots by integrating topology optimization and parametric system optimization
CN110103220B (zh) 机器人高速高精度运动轨迹规划方法、装置、设备及介质
CN112904728B (zh) 一种基于改进型趋近律的机械臂滑模控制轨迹跟踪方法
CN110539302B (zh) 一种工业机器人整体动力学建模及动力学参数辨识方法
CN111159636B (zh) 一种柔性多体系统动力学半解析灵敏度分析方法
CN110076783B (zh) 一种基于模糊控制和李雅普诺夫函数的平面欠驱动机械臂位置控制方法
CN110414165A (zh) 一种基于全局应力约束的多相材料柔顺机构拓扑优化方法
CN113569359B (zh) 一种基于惯性释放法的双足机器人结构拓扑优化方法
CN111168682B (zh) 并联机器人结合鲁棒精确微分器有限时间收敛滑模控制方法
CN112099516B (zh) 一种基于模糊线性化理论的卫星集群姿态跟踪控制方法
CN109507886A (zh) 针对时滞不确定性系统执行器故障的鲁棒预测容错控制方法
CN111897224A (zh) 基于演员-评论家强化学习和模糊逻辑的多智能体编队控制方法
CN113910238A (zh) 机器人刚度建模、辨识与修正方法及实验系统
Björkenstam et al. Enhancing digital human motion planning of assembly tasks through dynamics and optimal control
CN114578740A (zh) 一种基于改进自抗扰控制的软体驱动器控制方法
CN113377014A (zh) 一种机械臂系统鲁棒镇定控制方法和系统
Wang et al. Modeling and simulation of closed loop multibody systems with bodies-joints composite modules
Xu et al. Adaptive learning control of robot manipulators via incremental hybrid neural network
Xiong et al. Dynamic surface control-backstepping based impedance control for 5-DOF flexible joint robots
CN114186452A (zh) 一种非光滑多体系统动力学灵敏度分析方法
CN113171271A (zh) 一种上肢康复机器人的重力补偿方法
Mohammadi et al. Determining the state of a nonlinear flexible multibody system using an unscented Kalman filter
Redon et al. An efficient, error-bounded approximation algorithm for simulating quasi-statics of complex linkages
Jung et al. The operational space formulation on humanoids considering joint stiffness and bandwidth limit
CN114211478A (zh) 一种模块化机械臂协调操作最优控制方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant