柔性多体系统

柔性多体系统（精选八篇）

柔性多体系统 篇1

随着科技的进步, 对柔性多体系统的研究越来越广泛[1,2,3,4]。目前柔性多体系统的建模方法主要有以牛顿-欧拉方法为代表的矢量力学方法、以拉格朗日力学为代表的分析力学方法和以凯恩方法为基础的兼顾矢量力学与分析力学优点的建模方法[5,6]。Carrera[7]利用Lagrange方程和有限元方法研究了柔性机器人的反向动力学, 解决了刚体运动与柔性变形的耦合问题, 但算法为O (N3) 量级;Znamenácek等[8]利用递推算法推导了柔性多体系统动力学方程, 比较了“S”算法和“L”算法的计算效率, 得到O (N2) 量级的柔性多体系统动力学方程;Hwang[9]利用广义牛顿-欧拉递推算法建立了链式柔性多体系统动力学方程。

Jain等[10]发展了基于空间算子代数 (spatial operator algebra, SOA) 的多体系统动力学方法, 在算法上提高了计算效率, 该方法为O (N) 量级, 有效地解决了航天飞行器的计算问题。吴洪涛等[11]、熊启家[12]和方喜峰等[13]研究了空间算子代数理论, 并且进行了程序实现。但是在实际工程应用中, 对多体系统的建模不会是单纯的正向动力学和反向动力学, 而是在整个系统中既有关节力矩又有关节加速度, 这种系统常被称为欠驱动多体系统, 在大型航天系统中比较常见, 利用递推形式的正向和反向建模方法无法进行计算。Jain等[14]和Sohl等[15]分别利用利用空间算子代数和伴随算子研究了多刚体系统的混合递推动力学建模, 而尚未推广到多柔体系统建模中。

目前传统的求解微分方程有两种数值积分方法[16]:隐式法和显式法。但一般的积分算法对所建立的微分方程求解效率较低。翟婉明[16]所提出的线性多步快速积分算法已经成功地应用于解决列车大规模多体系统动力学问题的研究, 取得了很好的效果, 对求解大型微分代数方程具有很好的借鉴和启发作用, 为柔性多体系统动力学实时仿真打下基础。

本文以建立柔性多体系统高效率模型为目的, 在空间算子代数理论的基础上, 建立了一套广义柔性多体系统高效率的O (N) 量级混合递推动力学算法, 并且借鉴和移植翟婉明教授的线性多步积分算法, 对广义柔性多体系统微分代数方程进行研究, 编制了可以计算柔性多体系统反向、正向以及混合动力学建模的Mathematica软件包HybridDynamics.m。

1 柔性多体系统空间算子代数理论体系

SOA是研究多体系统的新型数学工具, 其关键创新点在于综合了看上去与力学无关的Kalman滤波方法和基于算子的递推计算方法。SOA基本算子有若干个, 每个算子采用空间递推方式实现, 其运算的数量与系统的自由度N同阶, 故称为O (N) 阶算法。

1.1 柔性多体系统模型

本文研究的对象是具有由n个铰链将n个柔性机械臂连接成的柔性多体系统, 如图1所示。铰链和操作臂的序号从系统的基座向顶端增大。k号铰链表示的是k号机械臂与k+1号机械臂的连接铰链。当k号铰链是主动关节时记为Ia, 否则记为Ip。n+1号铰链可固定到n号体的任意位置。其中的k号柔性体是一个由有限个节点 (在空间位置中定义i) 组成的有限元模型。空间位置可以由连接到物体上的坐标系表示出来。k号柔性体有限元模型对应一个质量矩阵Mm和一个刚度矩阵Km, 假定柔性体的质量矩阵和刚度矩阵是不随时间变化, 可预先计算得到的独立变量。

1.2 基于空间算子代数反向动力学递推[17]

采用空间算子代数理论, 对链式柔性多体系统而言, 可获得反向动力学向外和向内高效递推算法。在本研究中, 为提高计算效率所求解的速度和加速度都是在本体坐标系下定义的。

基于空间算子代数柔性体的速度递推和加速度递推模型如下:

基于空间算子代数的广义力递推和力矩递推模型如下:

式中, v (k) 为k体在本体坐标系下的刚体速度;Φ为柔体k-1到柔体k的模态移位算子;H为k体的模态坐标;q (k) 为k体的广义坐标矢量;α (k) 为k体的刚性加速度项;αm (k) 、am (k) 分别为k体的离心加速度项、科氏加速度项;ϕ (tk, k+1) 为k体经有限元离散后与k+1铰链连接的元素到k+1体的移位算子;bm (k) 为k体离心力项;fm (k) 为k体受到的广义力。

1.3 基于空间算子代数正向动力学递推[17]

按照牛顿定律, 多体系统动力学线性方程可以表示为

式中, undefined为惯性力项;T为主动力项;z为系统的速度的矢量积项, 冗余力项。

根据空间算子代数理论[17], 柔性多体系统广义质量按照LDLT分解可以分解为M=[I+HΦK]D[I+HΦK]T, 进而可以得到系统的广义质量的逆M-1=[I-HΨK]TD-1[I-HΨK]。因此式 (3) 用算子表达的形式得到正向动力学表达式为

undefined

2 柔性多体系统混合动力学建模理论

本节基于空间算子代数理论主要研究混合动力学高效率建模的方法。按照图1表示, 柔性系统是含有主动关节和被动关节的欠驱动系统, 按照传统的正向、反向动力学建模方法不能进行动力学计算, 需要研究适合求解欠驱动系统问题的高效率混合动力学建模方法。本文研究的混合递推动力学建模方法同正向动力学建模相似, 可以包括如下步骤:

(1) 从系统顶端n体到基座1体顺序按照铰接体惯量方法递推计算系统的惯量算子:

(2) 按从顶端n体到基体1体的顺序递推计算速度一次项的冗余力项:

(3) 按从基体1体到顶端n体的顺序递推计算系统各个体的加速度项:

以上便为欠驱动柔性多体系统混合递推动力学建模算法。本算法可看作是一般的柔性多体系统广义动力学建模算法, 其主要特征是同时包含正向动力学和反向动力学两种算法。当系统中全是被动关节时, 该算法变成正向动力学算法;全是主动关节时, 该算法变为反向动力学算法;其余为混合动力学建模方法。

3 数值分析

上文介绍了基于空间算子代数的柔性多体系统混合递推动力学建模方法, 显而易见, 如果在柔性系统中既含有主动关节又含有被动关节, 则根据上述算法建立的动力学方程是微分-代数方程。写成状态方程表示为

根据上述建立的动力学数学模型, 本文在翟婉明[16]所提出的线性多步快速积分算法的基础上研究了求解此微分-代数的数值积分算法。翟婉明所提出的线性多步快速积分算法已经成功地应用于解决列车大规模柔性多体系统动力学问题的研究, 并取得了很好的效果。本文借鉴和移植这类先进的计算方法, 对欠驱动柔性多体系统微分代数方程系统进行研究。研究表明, 根据线性多步快速积分算法求解微分-代数方程在积分精度和积分效率上同样是正确的、有效率的。

在起步阶段, 确定微分方程在t=0时刻的初始条件为q (k, 0) 、undefined, 然后代入到式 (8) 中即变成线性代数方程进而求得undefined、F (k, 0) 。

根据线性多步快速积分方法, 确定如下递推算法:

式中, Δt为时间积分步长;ψ、η为控制方法特性的独立参数, 具体参数确定参考文献[16]。

需要注意, 在起步阶段 (n′=0, 1) 时, 可使ψ=η=0从而按照递推方式进行计算。

将上述值代入到系统的动力学方程中, 即可得到系统n′+1状态较精确的广义加速度值以及广义主动力值:

因此按照可积分递推式 (9) 和式 (10) 逐次将欠驱动柔性系统的大型动力学微分-代数方程转换为多步线性代数方程进行计算, 进而得出对应各步的位移、速度、加速度以及主动力的离散值, 并且求解过程简单, 快捷, 不占用计算机内存。

4 混合动力学软件编制与实例仿真

4.1 软件编制

为验证上述理论的正确性和高效性, 笔者的理论推导过程是采用符号推导软件Mathematica 5.2得到的, 并且编制了求解一般多柔体系统的广义混合符号推导软件包HybridDynamics.m。核心的计算推导模块包括四个部分:第一部分定义柔性体以及铰链的数据结构;第二部分为根据柔性体和铰链信息定义关键空间算子;第三部分为计算所定义柔性系统的速度以及计算每个柔体的科氏力加速度项以及科氏力;第四部分为混合递推动力学的核心算法。

根据定义的柔体和铰链的信息, 计算混合递推过程需要的关键空间算子ϕ (x, y) 、Φ (k+1, k) 、HT (k) 等。随后按照文献[17]计算系统的速度、科氏加速度和科氏力。

柔性多体系统混合递推动力学核心算法程序流程图见图2。整个软件计算过程简洁、明了, 可以计算广义柔性多体系统动力建模。

4.2 实例分析

为验证上述混合动力学建模以及线性多步积分算法的准确性和有效性, 本文基于Mathematica软件平台, 以PUMA560机器人为研究对象, 对欠驱动柔性多体系统进行了O (N) 阶动力学建模和实时仿真进行了研究。柔性PUMA560机器人的柔性操作臂预先通过有限元进行离散来插值其弹性变形。PUMA560柔性机器人参数如表1所示。

在本仿真中, 时间积分步长Δt选择0.005s, ψ=η=0.5, 求解结果结果如图3～图8所示。

图3和图4所示为空间机器人本体的质心位置变化曲线;图5和图6所示为机器人第一关节的速度变化曲线;图7和图8所示为机器人第二关节的加速度变化曲线。通过快速积分算法和NDSolve算法的比较, 快速积分算法的精度是稳定的。

4.3 准确度度与高效率验证

针对上述提出的混合递推动力学建模方法, 本文以PUMA560机器人为例对其进行了正确性验证。验证主要分成三个步骤:①将所研究的柔性多体系统中的所有关节铰设定为主动铰, 系统动力学建模变为反向动力学建模过程, 用本算法与传统的反向动力学建模进行比较;②与步骤①类似, 将系统中的所有关节铰变为被动铰, 系统动力学建模变为正向动力学建模过程, 用本算法与传统的反向动力学建模进行比较;③将柔性多体系统某些铰设定为主动铰, 其余设定为被动铰, 进行动力学建模, 然后将系统中的主动铰变为被动铰, 被动铰变为主动铰再进行仿真, 将两种动力学建模结果进行对比。经过上述三步的正确性验证, 本文提出的基于SOA的高效率混合递推动力学建模方法是正确的、高效的。

在保证精度的情况下, 本文对两种算法的计算效率进行了比较研究。在同等计算条件下 (计算机配置为:Intel (R) Core (TM) 6400 2.31GHz的CPU;0.99GB的内存) , 用Mathematica 中的NDsolve命令求解微分方程所用的时间为11min32s, 并且计算过程占用大量的内存空间;用本文提出的快速积分算法求解微分方程耗时6min28s, 在整个计算过程中很少占用计算机内存空间, 其计算效率提高将近一倍。

5 结束语

(1) 在空间算子代数理论的基础上研究了广义的柔性多体系统高效率O (N) 阶建模的算法。该算法从实际工程出发, 突破了传统的纯粹柔性多体系统正向和反向动力学建模方法。该算法自动地判断系统中铰链的运动信息, 快速、高效、准确地建立柔性多体系统混合动力学模型。本算法可以应用于一般多柔体系统、欠驱动系统, 为大型柔性航天器动力学建模打下坚实的基础。

(2) 以大型柔性多体系统实时动力学仿真为目的, 在翟婉明教授提出的求解大型微分方程的线性多步积分算法的基础上, 研究了求解由混合递推算法建立的大型微分-代数方程的线性多步快速积分算法。

摘要：进一步发展了空间算子代数理论体系, 采用空间算子描述了广义柔性多体系统动力学高效率建模以及实时仿真问题。根据不同类型的铰链特征 (主动关节、被动关节) 描述了广义柔性多体系统特征, 按照两次从系统顶端到基座和一次从基座到顶端的顺序分别计算系统的广义铰接体惯量算子、系统的冗余力算子以及广义加速度和广义主动力矩, 进而建立广义柔性系统O (N) 阶动力学模型。采用线性多步积分算法理论解决了大型微分代数方程的数值积分算法, 实现了实时动力学仿真的目的。最后通过实例结果对比验证了研究内容的正确性和高效性。

柔性多体系统 篇2

机械多体系统动力学模型数值算法与违约修正

首先简要介绍了多体系统动力学的建模方法,总结了几种常用的动力学方程;其次,重点阐述了动力学模型数值求解算法及违约修正的最新研究进展;最后展望了多体系统动力学今后的发展方向.

作 者：高海涛 张志胜 史金飞 GAO Hai-tao ZHANG Zhi-sheng SHI Jin-fei 作者单位：东南大学,机械工程学院,江苏,南京,211189刊 名：中国制造业信息化 ISTIC英文刊名：MANUFACTURING INFORMATION ENGINEERING OF CHINA年，卷(期)：200938(11)分类号：O313关键词：多体系统 数值算法 违约修正 进展

刚柔耦合多体系统动力学模型降阶 篇3

(南京航空航天大学机械结构力学及控制国家重点实验室， 江苏 南京 210016)

引 言

随着空间技术的发展，航天器柔性构件的尺寸越来越大，柔性构件的变形对航天器的动力学行为产生很大的影响，传统多刚体系统动力学理论已经不能满足人们对设备精度的要求。最近几十年，考虑构件柔性的柔性多体系统动力学，已经成为国内外众多学者研究的热点并取得了大量的研究成果。考虑构件弹性变形与大范围刚体运动的耦合，Likins提出了浮动坐标方法[1]，该方法将构件的位形认为是浮动坐标系大范围刚体运动与相对于该局部坐标系的变形的叠加。1987年，Kane在研究梁的高速旋转运动时第一次发现了动力刚化现象[2]。为解决动力刚化问题，Haering等在多体系统动力学建模过程中考虑了高阶项[3,4]。

然而，浮动节点坐标法是基于小变形假设，对于存在大变形的多体系统已经不再适用。Shabana提出了目前广泛应用于分析大变形柔性多体系统动力学的绝对节点坐标方法(ANCF)[5]。该方法中单元节点的坐标定义在全局坐标系下, 采用斜率矢量代替传统有限单元中的节点转角坐标，推导建立的多体系统微分-代数方程的质量矩阵为常数矩阵，且具有不存在科氏力和离心力项的优点。Berzeri等提出了几种基于不同假设的一维梁的弹性力简化模型[6]，并做了比较研究。Escalona等首先将该单元应用于柔性大变形多体系统动力学的研究[7]。Omar等放宽梁中线的切线矢量与梁截面的法线方向重合的的假设[8]，将梁的剪切变形考虑到梁单元中,首次提出了一种平面应变剪变梁单元。该单元由于弯曲应变与轴向应变不一致而带来剪变闭锁问题。为了解决剪变闭锁问题，Kerkkänen等通过改变单元的运动学描述[9,10],提出了一些可有效减轻剪变闭锁问题的线性剪变梁单元。考虑到绝对节点坐标体系下刚度矩阵的强非线性，导致采用绝对坐标方法建立的微分-代数方程的计算效率比较低, García-Vallejo提出不变矩阵法[11], 该方法将非线性刚度矩阵分解为常数刚度矩阵与广义坐标相关的刚度矩阵之和。

绝对节点坐标法虽然能够精确描述多体系统的刚性运动，但是即使是刚体也要划分单元[12]，这必然导致该方法较难处理刚柔耦合多体系统动力学问题。Sugiyama等结合浮动节点坐标法能有效处理刚体运动和绝对节点坐标法能描述柔性体大变形的特点[13～15]，使存在柔性体大变形的刚柔耦合多体系统得到了解决。然而，浮动节点坐标法建立的多体系统动力学方程的质量矩阵与广义坐标相关，因而得到的刚柔耦合多体系统的微分-代数方程的质量矩阵也与广义坐标相关，每次计算都需要对质量阵进行计算，会大大降低计算效率。自然坐标法以其具有建立的多体系统微分-代数方程的质量矩阵为常数矩阵的优点而成为刚柔耦合多体系统中替代浮动节点坐标法用于描述刚体构件的方法。García-Vallejo用自然坐标法与绝对节点坐标法的混合方法对平面刚柔耦合多体系统动力学进行了研究[16]，并且对构件的各种连接情况进行了讨论。在此基础上，García-Vallejo进一步对空间刚柔耦合多体系统动力学问题进行了研究[17]。

上述研究随着柔性体单元数量的增加，刚柔耦合多体系统动力学方程的计算效率将会比较低。为了提高计算效率，需要对绝对节点坐标法建立的柔性多体系统进行模型降阶。传统的模型降阶方法主要有子结构方法[18～21]，Krylov子空间方法等[22,23]。Hurty首先提出了模态综合法的概念[18]，并且应用于对大规模线性系统进行模型降阶。R R Craig和C C Bampton对此方法作了部分修正[19]，形成了现在的固定界面模态综合法。随后，Kobayashi成功将Craig-Bampton固定界面法应用于基于绝对节点坐标法的柔性多体系统的模型降阶[24]。本文采用Craig-Bampton固定界面法对基于绝对节点坐标法和自然坐标法建立的刚柔耦合多体系统进行模型降价。

1 绝对节点坐标法

1.1 基于绝对节点坐标法的梁单元

绝对节点坐标法中，柔性体k的一维梁单元j，如图1所示，该单元上任意一点全局位置为

rkj=Skjekj

(1)

图1 平面梁单元

式中Skj为单元形函数，ekj为单元节点的广义坐标矢量，可表示为

基于上述描述，梁单元的动能可表示为

(2)

(3)

式中εkj和κkj分别为单元j中线上对应点的应变和曲率。

基于虚功原理，可得到该单元的动力学方程为

(4)

(5)

(6)

(7)

1.2 模型降阶

(8)

用Craig-Bampton方法进行减缩时，约束模态仅取前nc阶，则广义坐标ek可减缩为

ek=Tkqk

(9)

(10)

将式(9)代入式(8)，等式两边左乘TkT，则有

(11)

2 自然坐标法

自然坐标法中，用两个基点的绝对坐标矢量描述刚体i的位置和方向，如

(12)

式中di为包含基点C和D的刚体i的坐标矢量(如图2所示)，该坐标矢量是非独立的，C和D之间存在距离约束，有

(13)

图2 两基点刚体

基于上述刚体描述，自然坐标法建立的刚体系统动力学方程为

(14)

固结在刚体上的动坐标系的旋转矩阵可表示为

(15)

3 刚柔耦合多体系统

柔性体之间、刚体之间、柔性体与刚体之间存在各种约束。本文仅描述刚体与柔性体之间的旋转副约束和固结约束(如图3所示)，其他相关约束可参考文献[16]。

图3 柔性体与刚体之间的约束

3.1 旋转副约束

假设刚体i的P点与柔性体k的节点n存在旋转副约束(如图3(a)所示)，约束方程可表示为

(16)

(17)

式中

3.2 固结约束

假设刚体i的P点与柔性体k的节点n存在固结约束(如图3(b)所示)，则约束点除了存在位置约束外还存在方向约束，即

(18)

(19)

式中

3.3 消除线性约束方程

将约束节点作为柔性体的边界，则柔性体的边界节点广义坐标可以用连接刚体的自然坐标和柔性体边界的减缩坐标表示。假设刚体i的P点与柔性体k的节点n存在旋转副约束，则柔性体k经Craig-Bampton方法减缩后的广义坐标可表示为

(20)

由式(11)和(14)得刚柔耦合多体系统动力学方程为

(21)

由式(20)可以对广义坐标P进行减缩

(22)

将式(22)代入式(21)，等式两边左乘TT，则

(23)

假设刚体i的P点与柔性体k的节点n存在固结约束，同理可得消除边界约束后的刚柔耦合多体系统动力学方程为

(24)

4 数值算例

图4为受重力的平面双摆。相关参数参考文献[25],A点与B点均为旋转铰，梁AB与梁BC的长度均为1.8 m，截面积为2.5×10-4m2，密度为2.766 67×103kg/m3。梁AB与梁BC都可以作为刚体或柔性体。本文首先基于3种情况对该双摆系统进行动力学分析：(1)梁AB与梁BC均为柔性体，用绝对节点坐标法对其进行研究；(2)梁AB与梁BC均为刚体，用自然坐标法对其进行研究(NCF)；(3)梁AB为柔性体，梁BC为刚体，用文献[16]的平面刚柔耦合多体系统动力学方法对其进行研究(NCF-ANCF)。如果梁为柔性体，将梁等分为10个单元，杨氏模量为6.895×109Pa，截面惯量矩为1.302×10-9m4，取重力加速度为9.81 m/s2，仿真时间为2.5 s。在一台具有Intel Pentium 3.2 GHz处理器及3GB RAM的PC机上运行。3种情况下，C点Y方向的绝对位移如图5所示

图4 双摆

图5 C点Y方向绝对位移

由图5可知，3种情况下C点Y方向的位移存在差异，这说明弹性变形对双摆端点C的运动会产生影响，但3种情况下C点位移相差不大，因此，在精度要求不高的情况下，可以把梁作为刚体考虑。

多体系统按照构件在运行过程中的变形，可以将构件分为刚体构件和柔性体构件，其中，随着柔性体构件划分单元数量的增加，刚柔耦合多体系统动力学方程会存在计算时间长，计算效率低的问题。因此，本文将模型降阶的方法引入刚柔耦合多体系统动力学，提出基于模态综合法的刚柔耦合多体系统减缩方法，解决了传统存在大变形的刚柔耦合多体系统动力学方程计算效率低的问题。为了验证该方法的正确性和有效性，本文以图4的双摆为例，将梁AB作为柔性体，取其杨氏模量为6.895×108Pa，梁BC作为刚体，用该方法对此系统进行研究。在选取主模态集时，分别取前5阶模态(C-N-A(5))，前7阶模态(C-N-A(7))和前9阶模态(C-N-A(9))，将其与原模型(N-A)进行对比。减缩模型和原模型C点Y方向位移和柔性梁端点横向变形分别如图6，7所示,计算所耗CPU时间如表1所示。

表1 计算所耗CPU时间

由图6可知，取较少的主模态就能使C点Y方向绝对位移误差很小，这是因为柔性梁的振动主要为低频振动，因此，只用选取较少的模态就能粗略描述柔性梁的弹性变形。图7中，原模型的柔性梁端点横向变形曲线与取前9阶模态时的柔性梁端点横向变形曲线几乎重合，而与取前5阶、前7阶模态时的柔性梁端点横向变形曲线存在明显差异。由此可见，只有适当选取更多模态才能更好地描述柔性体的弹性变形。由表1可知，减缩模型的计算时间都比原模型的计算时间要少，而且选择的模态数量越少越节省计算时间。结合图7和表1可以得到：减缩模型仅选取前9阶模态时，能够在满足计算精度的情况下使计算时间仅为原模型计算时间的34.9%。由上述分析可知，适当选取模态就可以在保证计算精度的情况下，减少刚柔耦合多体系统动力学方程的计算时间，从而提高计算效率。

图6 C点Y方向绝对位移

图7 柔性梁端点横向变形

5 结 论

本文考虑到存在大变形的刚柔耦合多体系统动力学方程计算效率比较低，提出了基于模态综合法的刚柔耦合多体系统模型降阶方法。通过双摆系统对该方法的正确性和有效性进行了研究。由数值仿真可知，减缩模型随着选择模态数量的减少，同原模型相比越节省计算时间，而且只要适当选择减缩模态就可以保证计算精度。这说明该方法能够提高刚柔耦合多体系统动力学方程的计算效率。

参考文献：

[1] Likins P W. Finite element appendage equations for hybrid coordinate dynamic analysis[J]. Journal of Solids and Structures, 1972, 8: 709—731.

[2] Kane T R, Ryan R R, Banerjee A K. Dynamics of a cantilever beam attached to a moving base[J]. Journal of Guidance, Control, and Dynamics, 1987, 10(2): 139—151.

[3] Haering W J, Ryan R R, Scott R A. A new flexible body dynamic formulation for beam structures undergoing large overall motion[A]. 33rd Structural Dynamics and Materials Conference[C]. Dallas: AIAA, 19921415-1423.

[4] Mayo J M, García-Vallejo D, Domínguez J. Study of the geometric stiffening effect: comparison of different formulations[J]. Multibody System Dynamics, 2004, 11: 321—341.

[5] Shabana A A. Definition of the slopes and the finite element absolute nodal coordinate formulation[J]. Multibody System Dynamics, 1997, 1: 339—348.

[6] Berzeri M, Shabana A A. Development of simple models for the elastic forces in the absolute nodal coordinate formulation[J]. Journal of Sound and Vibration, 2000, 235(4): 539—565.

[7] Escalona J L, Hussien H A, Shabana A A. Application of the absolute nodal coordinate formulation to multibody system dynamics[J]. Journal of Sound and Vibration, 1998, 214(5): 833—851.

[8] Omar M A, Shabana A A. A two-dimensional shear deformable beam for large rotation and deformation problems[J]. Journal of Sound and Vibration, 2001, 243(3): 565—576.

[9] Kerkkänen K S, Sopanen J T, Mikkola A M. A linear beam finite element based on the absolute nodal coordinate formulation[J]. ASME Journal of Mechanical Design, 2005, 127: 621—630.

[10] García-Vallejo D, Mikkola A M, Escalona J L. A new locking-free shear deformable finite element based on absolute nodal coordinates[J]. Nonlinear Dynamics, 2007, 50: 249—264.

[11] García-Vallejo D, Mayo J, Escalona J L, et al. Efficient evaluation of the elastic forces and the Jacobian in the absolute nodal coordinate formulation[J]. Nonlinear Dynamics, 2004, 35: 313—329.

[12] 田强，张青云，陈立平，等. 柔性多体系统动力学绝对节点坐标方法研究进展[J]. 力学进展, 2010, 40(2): 189—202.TIAN Qiang, ZHANG Yunqing, CHEN Liping, et al. Advances in the absolute nodal coordinate method for the flexible multibody dynamics[J]. Advances in Mechanics, 2010, 40(2): 189—202.

[13] Lee S, Park T, Seo J, et al. The development of a sliding joint for very flexible multibody dynamics using absolute nodal coordinate formulation[J]. Multibody System Dynamics, 2008, 20: 223—237.

[14] Hussein B A, Weed D, Shabana A A. Clamped end conditions and cross section deformation in the finite element absolute nodal coordinate formulation[J]. Multibody System Dynamics, 2009, 21: 375—393.

[15] Sugiyama H, Yamashita H. Spatial joint constraints for the absolute nodal coordinate formulation using the non-generalized intermediate coordinates[J]. Multibody System Dynamics, 2011, 26: 15—36.

[16] García-Vallejo D, Escalona J L, Mayo J, et al. Describing rigid-flexible multibody systems using absolute coordinates[J]. Nonlinear Dynamics, 2003, 34: 75—94.

[17] García-Vallejo D, Mayo J, Escalona J L, et al. Three-dimensional formulation of rigid-flexible multibody systems with flexible beam elements [J]. Multibody System Dynamics, 2008, 20: 1—28.

[18] Hurty W C. Dynamic analysis of structural systems using component modes[J]. AIAA Journal, 1965, 3(4): 678—685.

[19] Craig R R, Bampton M C C. Coupling of substructures for dynamic analysis[J]. AIAA Journal, 1968, 6(7): 1 313—1 319.

[20] Hou S N. Review of modal synthesis techniques and a new approach[J]. Shock and Vibration Bulletin, 1969, 40(4): 25—30.

[21] 王文亮，杜作润，陈康元. 模态综合技术短评和一种新的改进[J]. 航空学报, 1979, 3(2): 32—51.

[22] Antoulas A C. Approximation of Large-scale Dynamical Systems[M]. Philadelphia: SIAM, 2004.

[23] Bai Z. Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems[J]. Applied Numerical Mathematics, 2002, 43: 9—44.

[24] Kobayashi N, Wago T, Sugawara Y. Reduction of system matrices of planar beam in ANCF by component mode synthesis method[J]. Multibody System Dynamics, 2011, 26: 265—281.

柔性多体系统 篇4

考虑最一般的情况,含接触/碰撞的多体系统动力学过程可以分为3个阶段:未接触段、接触/碰撞阶段和分离后的阶段(如图1所示).在接触域确定后,将接触作用视为接触力元[4,5,6]或接触约束[1,7,8,9,10,11,12,13],而这3个阶段的不同之处仅在于多体系统的拓扑结构,故含接触/碰撞的多体系统又称为变拓扑系统[1,2,13,14].对于未接触阶段、接触/碰撞阶段和分离后的阶段,变拓扑系统的动力学方程是统一的,3个阶段的区别仅在于接触力元或接触约束的增加或删除.

对应于上述3个阶段,可以总结出柔性多体系统接触/碰撞动力学中的4个关键问题.接触/碰撞阶段涉及到的最主要的问题为:接触区域的处理、接触作用的力学模型.3个阶段都涉及到的问题为:柔性体的运动学描述、约束问题的数值方法.本文将综合概括柔性多体领域处理上述4个问题的现有方法,一并考虑了有限元中对相应问题采用的方法.

1 接触区域的处理

如图2所示,两个发生接触的物体记为Bk,k=1,2.记接触区域的特征尺寸为d,物体的特征尺寸(曲率半径)为Lk,由接触造成的压缩量为δk,参数

若对计算精度要求较低,且满足ξk<<1,可以用简化的方法处理.在简化方法中,接触域的全部信息完全由两个重合点和接触局部区域的其他一些几何量(公共法线矢量、曲率半径、扭率)表达[1,10,11,15].若接触域尺寸不能满足ξk<<1,或者想得到更高精度的计算结果,则要采用精细的方法处理接触域.精细方法对接触域进行细致的空间离散化,对接触/碰撞进行细致的动力学分析.

1.1 简化方法

如果对物体的变形和应力不关心,可以对接触物体采用刚体模型,则两物体的接触关系可以由点接触约束方程描述[2,15,16].刚体的点接触约束可以视为一随时间变化的“虚铰”,应用此思想,Cai等[10]和Montana[11,17]均给出了两个任意形状物体的滑动点接触和非滑动点接触的约束方程.接触约束方程与动力学方程联立,可求得接触力.对于碰撞问题,由于采用刚体模型,碰撞过程变为一非连续问题或非光滑问题,需引入恢复系数来使方程封闭[12,13,14].

若对物体的变形和应力分布比较关心,可以应用半空间近似,由Boussinesq解在接触域上积分,得到接触力合力F与压缩量δ的解析关系式F(δ,Lk)[18],即将接触作为一个接触力元处理,也就是常用到的Herz方法.显然,半空间近似忽略了许多因素,只能得到近似的解.

简化方法极大地简化了接触/碰撞问题的求解,提高了计算效率,而且许多工程实际问题满足ξk<<1,因而简化方法得到了广泛的应用.多体系统动力学侧重于系统大范围运动和长期动力学过程的仿真,注重计算效率,且多体系统中许多接触/碰撞问题满足ξk<<1,因而长期以来对接触域的处理一直采用简化方法.

简化方法的缺点为:只能用于ξk<<1的问题,精度较低,不能反映出接触域的时间演化历程,很难处理摩擦问题,特别是在接触/碰撞过程中摩擦状态发生改变的问题.近些年来,随着计算机硬件水平的提升和对仿真精度要求的提高,柔性多体系统领域内的研究者借鉴有限元领域的方法,对接触区域采用精细的处理方法,对接触碰撞过程进行细致的动力学分析[4,8,9,19].

1.2 精细的方法

如图2所示,为了提高计算精度,精细方法将接触域细分为多个曲面,在细分后的接触域上,两物体(分别称为主物体和从物体)上发生接触的部分构成许多个接触对,接触对间形成接触力元或接触约束.目前常用的接触对有点/面接触对、面/面接触对等[16,19,20].

点/面或面/面接触对能细致地反映接触域的时间演化历程.在接触/碰撞发生时刻,接触域一般较小,少量的接触对开始起作用.在压缩阶段,接触区域逐渐增大,越来越多的接触对开始起作用.在恢复阶段,接触域又逐渐变小,接触对间依次分离.考虑摩擦时,点/面和面/面接触对能反映出接触域间的滑移或粘连的演化过程.

精细方法对ξk无限制,适用于任意复杂界面的接触问题,且计算精度高,能给出物体内部(尤其是接触区域附近)高精度的应力、应变响应[16,19,20].目前在有限元、边界元或者影响函数[18]方法处理接触/碰撞问题时,均采用精细方法对接触域进行处理,以得到高精度计算结果.

2 接触作用的力学模型

接触/碰撞过程中,两物体通过接触发生相互作用.从多体系统动力学的角度来看,系统中物体的相互作用有两种力学模型,力元和约束(多体系统中一般为铰约束).同样地,接触/碰撞过程的力学模型有两种:接触力元和接触约束.相应的,目前处理接触/碰撞过程的方法有两种:接触力元方法[6,16,21]和接触约束方法[8,16,22].接触力元方法将接触作用视为弹簧阻尼器力元,力元根据接触约束违约量直接由相应模型给出接触力,又称为罚函数法或接触力模型方法.接触约束方法将接触作用视为接触约束,接触约束方程与带Lagrange乘子的动力学方程形成封闭的方程组,该方法又称为Lagrange乘子法或附加接触约束法.

在接触/碰撞问题中,接触虚功项的表达式为式中,Γc为接触域,为法向接触项,为切向接触项.

两种力学模型的不同在δWc的数学表达式中体现.本节讨论限于法向约束.

2.1 接触力元方法

接触力元方法中,检测到接触发生时,在接触区域添加相应的接触力元,直至分离.若简单地将接触作用取为弹簧力元,εN为弹簧刚度(罚因子),接触项表达为

其中gN表示接触对的法向穿透量.

接触力元方法在柔性多体系统动力学和有限元方法中均有广泛的应用.早期在多体系统动力学中,对接触域的处理采用简化方法,接触力模型的系数是通过接触局部区域的几何量和刚度、阻尼比等值确定的,约束违约量是通过两物体接触区域整体的嵌入量和嵌入率确定的[6,21].目前在柔性多体系统动力学和有限元领域中,接触域的处理采用精细方法,接触力元方法是应用于离散后的各个接触对水平[16,23,24].因此,早期多体系统动力学中的接触力模型方法随模型结构和系数不同而有很大差异,而目前的接触力元方法却有较高精度.

接触力元方法的优点在于,不需要求解接触约束方程,且动力学方程的维数是固定的,不随接触状态发生改变.接触力元方法的不足之处在于,理论上只有当弹簧刚度(罚因子)取得足够大时,结果才可靠.但在隐式方法中,随着弹簧刚度(罚因子)的增大,方程的数值性态变差.在显式方法中,增大弹簧刚度(罚因子)需减小积分步长以获得平滑的接触力.

从物理本质上来看,接触力元方法依赖于弹簧刚度(罚因子)是因为接触力完全由弹簧力元提供.为了克服接触力元方法依赖于弹簧刚度(罚因子)的困难,目前广泛采用增广Lagrange乘子法[20].增广Lagrange乘子法中同样将碰撞过程视为力元的作用,但不再是简单的弹簧阻尼力元,而是一个常力力元加上一个弹簧阻尼器力元.通过迭代的方法,逐步确定常力力元的值,从而减小了对弹簧阻尼力元的依赖.

2.2 接触约束方法

接触约束方法中,检测到接触发生时,在接触区域施加相应的接触约束方程.接触约束方程与带Lagrange乘子的动力学方程联立可求得运动学变量和接触力.当计算得到的法向接触力符号改变时,接触分离.接触约束方法中,取

正如多体系统中的铰约束方程一样,位置约束方程与动力学方程构成指数3微分代数方程,带来数值上的困难.另外,从未接触到接触时,速度发生突变,数值方法上也遇到困难.接触约束的数值问题将在第4节中详细讨论.接触约束方法的另一问题是要求解的方程的维数是随接触状态而改变的,这给编制计算程序带来了困难.接触约束方法中有两种比较特殊的方法:冲量/动量方法,应用互补算法的方法.

2.2.1 冲量/动量方法

冲量/动量方法假设碰撞过程持续时间无穷短,碰撞过程中物体(或有限元节点)位置不变,由Newton恢复系数提供速度约束方程.或者将碰撞过程分为压缩阶段和恢复阶段,在压缩阶段结束时碰撞点速度相等,求解Lagrange乘子形式的动量方程可得到压缩段冲量,在恢复阶段由Poisson恢复系数或Stronge恢复系数可确定恢复阶段冲量,代入Lagrange乘子形式的动量方程可得碰撞结束时刻物体的速度.通常又把冲量/动量方法称为不连续的方法,冲量/动量方法在柔性多体系统动力学和有限元法中均有应用[12,14,25].

恢复系数在刚体碰撞问题中反映了碰撞过程的能量损失,用于柔性体碰撞问题则缺乏物理依据,直接用于斜碰撞问题可能带来错误,且该方法不能得到碰撞力.为了解决刚体斜碰撞问题,有些学者提出了步进冲量法,该方法将碰撞过程细分为若干个子过程,在每个子过程可以判断切向冲量的大小和方向使其满足摩擦定律,碰撞过程何时结束可由恢复系数确定[26].

2.2.2 应用互补算法的方法

接触力和接触对的法向相对距离(或速率)、切向摩擦力和接触对相对滑移量(或滑移率)满足互补关系.两对互补关系与带Lagrange乘子的动力学方程构成接触/碰撞问题的互补形式的控制方程.一般假设主体与从体间的滑移量很小,应用事先确定的点/点接触对,则平面碰撞问题简化为线性互补问题,空间碰撞问题简化为非线性互补问题,可以应用互补算法方便地求解.互补方法在多体系统和有限元领域均有应用[15,27,28].互补方法的主要问题在于,很难用于一般的大滑移接触问题.

3 柔性体的运动学描述

柔性体的运动学描述是至关重要的,因为运动学变量的选择在很大程度上决定了动力学方程的性态和仿真计算的效率.与非接触问题相比,在接触/碰撞问题中,为了兼顾计算精度(尤其是接触局部区域的精度)与计算效率,柔性体的运动学描述有一些不同之处.以下将分别考虑两类问题的运动学描述.

3.1 非接触问题的柔性体运动学描述

依运动学变量的选取不同,描述柔性体运动主要有两类方法:浮动参考系法[29,30,31,32]和绝对坐标方法[31,33,34,35,36].浮动参考系法以浮动基位形坐标和柔性体相对浮动基的变形量为变量.在小变形条件下,浮动参考系方法的一个显著优点为可用模态坐标减缩变量维数,大大提高计算效率.

绝对坐标的选取与研究对象相关.对于柔性实体,绝对坐标方法以柔性体各节点在惯性参考系下的位置坐标为变量[37].对于梁、壳等柔性结构,引入适当的运动学假设后,柔性体上任意点的运动可以由参考线或参考面上点的位置坐标和固连于截面的一组基相对惯性基的转动来描述[38],通常称这样的梁壳模型为“几何精确”的梁壳模型.“几何精确”的梁、壳单元中,每个节点分别有6个和5个自由度.转动的运动学描述比较复杂,为了避免引入转动变量,一些研究者选取参考线或参考面上点的位置坐标和参考线或参考面在该点处的斜率矢量来描述梁壳的运动,如Shabana提出的绝对节点坐标方法.

“几何精确”的梁壳模型没有考虑截面的翘曲,而且广义应力应变关系一般很难得到.为了解决上述困难,Hodges[39,40]发展了由Berdichevs提出的变分渐进截面分析方法,通过对梁壳三维应变能进行变分渐进分析,得到翘曲位移与广义应变间的关系,将翘曲位移用广义应变表示,进而得到各阶近似下的广义应力和广义应变关系,该方法在直升机桨叶动力学分析中得到了广泛应用.

3.2 接触/碰撞问题的柔性体运动学描述

对于多体系统接触/碰撞问题,无论是浮动参考系法还是绝对坐标方法,若要提高计算精度,就需要采用更多的变量来描述柔性体的变形,这无疑会降低计算效率.为了平衡精度与效率的矛盾,主要有以下两种方法:分阶段描述与分区描述.

分阶段描述即在非接触/碰撞阶段和接触/碰撞阶段对柔性体采用不同的运动学描述.Ogot等[41]与Eberhard[19]提出了多刚体/有限元算法,对于需要长时间仿真的未接触阶段采用多刚体模型计算,在持续时间较短而精度要求较高的接触/碰撞阶段转化为有限元模型,接触分离后又转化到多刚体模型.进一步地,在未接触阶段可以考虑柔性体的变形,采用模态坐标形式的浮动参考系方法.分阶段描述的不足之处在于,接触/碰撞开始和分离后需要对两组运动学变量进行的相互转化,会带来误差.

分区方法对发生接触物体的非接触区域和接触局部区域,或者是对系统中发生接触的物体与非接触物体,采用不同的描述方法,该思想最早见Benson和Wu的文章[42,43].Benson提出了部分刚体/部分柔性体的方法,对接触/碰撞局部区域采用非线性有限元方法,远离接触/碰撞区域采用刚体模型.部分刚体/部分柔性体的方法同样见于Ambrósio和Lim等的文章[4,44],Ambrósio将部分刚体/部分柔体的方法用于多体系统[4],Lim等研究了部分刚体/部分柔体的计算方法[44].部分刚体/部分柔性体方法的不足之处在于,没有考虑非接触/碰撞区域的变形,计算精度低.Wu将柔性体碰撞过程分为局部区域的碰撞和随后的接触两个阶段,应用子结构方法求解多体系统无摩擦碰撞问题.子结构方法的不足之处在于,只能用于无滑移、无摩擦的简单碰撞问题.也可以采用分物体描述的方法,降低计算量,如Kim等[8]只对多体系统中直接发生接触的物体应用有限元方法.

4 约束问题的数值方法

若采用接触约束的方法,则柔性多体系统接触/碰撞问题涉及到两类约束:铰约束和接触约束.从细部结构来看,铰也是一种接触约束.但在多体系统动力学中,为使问题简化,一般认为铰位于物体上的一点.两种约束的区别在于:铰约束始终加在系统上,为双侧约束;接触约束为单侧约束,它是否起作用取决于两物体的接触状态,且接触约束作用的区域一般是不断变化的.

约束的存在给数值仿真带来困难.理论上,铰约束问题和接触约束问题均属于指数3微分代数方程问题.对于指数3微分代数方程,传统上采用多步法中的后向差分算法和单步法中的隐式荣格/库塔方法[45],近些年来,许多研究者采用有限元中常用的二阶算法如能量守恒的算法或能量耗散的算法(HHT算法、广义α算法、非连续的Galerkin方法)等[35,45,46,47,48].无论采用何种算法,求解指数3微分代数方程的主要困难有以下几点:Jaccobi阵的条件数过大、变步长算法中局部误差的阶减缩以及稳定性问题[45].一般在求解指数3微分代数方程之前要先对其进行处理,如降低指数[45,46,47]、对约束方程或Jaccobi阵进行缩放[45,46,49]或者通过选取独立变量将微分代数方程转变为常微分方程[32,50].

接触约束问题的另一困难为从未接触到接触时(即发生碰撞时),速度发生突变.不同于刚体碰撞问题,柔性体从未接触到接触时刻,只有接触局部区域速度发生突变.在多体系统动力学领域,Wu等[43]用碰撞局部区域的动量守恒方程解决发生碰撞时的速度突变问题.在有限元领域,不考虑速度突变而导致的一个问题就是,在未接触到接触的时间步内,施加某一水平(位置、速度、加速度)的约束,另外两个水平的约束不能同时满足[51].在施加位置接触约束的有限元隐式算法中,若不考虑碰撞发生时刻的速度突变,则计算得到的接触力振荡异常显著[20,51],与实际情况不符.如何在发生碰撞时同时满足(或近似满足)3个水平的约束,成为有限元处理碰撞的一个关键问题.目前的方法主要有弹性波理论校正法[52]、附加乘子法[51]与能量守恒算法[20].在有限元显式方法中(如预测/校正方法),在未接触到接触的时间步内,若需要对计算结果进行校正[22,53].

5 总结

含接触/碰撞的多体系统动力学过程可以分为3个阶段:未接触阶段、接触/碰撞阶段、分离后的阶段.3个阶段中主要涉及到4个关键的问题:接触区域的处理、接触作用的力学模型、柔性体的运动学描述和约束问题的数值方法.

本文针对上述4个问题进行了深入地分析,综合概括了柔性多体领域和有限元领域处理上述4个问题的现有方法,指出了这些方法的优点和不足.柔性多体系统碰撞动力学中,较为重要的问题还有摩擦问题、接触搜索等.

从上述论述,可以概括出柔性多体系统接触/碰撞问题的发展趋势为:

(1)借鉴有限元领域对接触域进行精细处理方法,提高接触/碰撞问题的计算精度,并能处理任意复杂界面的接触/碰撞问题;

(2)通过进行分阶段和分区的运动学描述,兼顾多体系统全过程仿真计算的效率与接触/碰撞阶段接触局部区域的计算精度;

柔性多体系统 篇5

数控机床作为工作母机在机械加工制造业中得到了广泛的应用[1]。随着数控加工技术的迅猛发展,传统的数控加工方法已经无法满足现代产品多样化、个性化的需求,现代数控机床向着高速、高精、高效、复合和环保的方向发展,以满足加工行业对零件加工精度不断提高的要求和对零件加工高速高效的不断追求[2]。复合加工机床(Complex Machine Tools)也称之为完全加工机(Complete Machining Machine Tools),其基本含义就是要在单台复合加工机床上实现零件的大部分或全部工序的加工[3,4,5]。随着机械加工市场不断增加的对复合加工机床的需求,国际上复合加工机床将进入激烈的竞争时代。使复合加工技术[6]成为推动机床结构和制造工艺发展的一个新热点,成为数控加工中心发展的重要方向之一。经过国内外学者多年的努力,多体系统理论得到了充分的发展,其通用性、系统性和方便性都有显著的提高,将其应用于数控机床几何误差建模中会大大简化研究过程[7]。本文以双主轴五轴高速加工中心为研究对象,从影响机床加工精度的几何误差着手,对数控机床几何误差建模理论进行分析和研究。

1 五轴数控机床结构分析与描述

以德国巨浪的DZ08FX双主轴高速五轴加工中心为例,可建立如图所示的机床结构图和拓扑结构图,从图中可以看出该机床总共有13个运动部件,将其分为工件分支和刀具分支两个部分,由四条运动链组成,1-2-3-4-6-7链为床身-刀具分支1,1-2-3-4-5-8链为床身-刀具分支2,1-9-10-11链为床身-工件分支1,1-9-13-12链为床身-工件分支2。

1.床身;2.X导轨;3.Z导轨;4.Y导轨;5.第二主轴;6.第一主轴;7.刀具1;8.刀具2;9.A转台;10.C转台;11.工件1;12.工件2;13.C’转台

2 数控机床运动学模型的建立

根据刀具与工件的不同组合,该机床总共有3种加工模式。分别为刀具与第一主轴组合(B-T2-T3-T4-T6-T7)和(B-W9-W10-W11),称为加工模式一,刀具与第二主轴组合(B-T2-T3-T4-T5-T8)和(B-W9-W13-W12)称为加工模式二。同时加工两个相同零件,同时工作。刀具与工件组合,称为加工模式三。由于加工模式一和加工模式二的建模方法类似,因此本文仅以加工模式一为例进行误差建模分析。

设工件坐标系为Ow-xwywzw,刀具坐标系为Otxtytzt,机床坐标系为O0-x0y0z0,P为刀具成形点,通过以上分析可得,P点在床身坐标系O0-x0y0z0中的位置可表示为:

典型体W(工件)体参考坐标系上被加工点P按“机床—工件”分支在与床身固连的惯性坐标系中的位置Pw的表达式:

式中:{rw}w为P点在典型体W(工件)坐标系Owxwywzw中的位置矩阵表达式;

{Pw}0为P点通过工件分支描述到与床身固连的惯性坐标系O0-x0y0z0中的矩阵表达式。

同理可得典型体T(刀具)上被加工点P分别按“机床—车削”、“机床—铣削”分支在与床身相连的坐标系中位置PTc、PTx的表达式:

式中(rTc)Tc为P点在典型体Tc(车削刀具)坐标系OTc-xTcyTczTc中的位置矩阵表达式;

{PTc}0为P点通过车削分支描述在与床身固连的坐标系O0-x0y0z0中的矩阵表达式。

式中(rTx)Tx为P点在典型体Tx(铣削刀具)坐标系OTx-xTxyTxzTx中的位置矩阵表达式;

{PTc}0为P点通过铣削分支描述在与床身固连的坐标系O0-x0y0z0中的矩阵表达式。

机床作成形运动时由于误差存在使得刀具成形点实际位置和理论位置出现偏差,这种偏差称为空间误差,对空间误差进行建模,是机床误差溯源的基础。刀具成形点空间误差为:

{Epc}、{Epx}分别表示车削、铣削模式下,机床的刀具加工点在惯性系中实际位置与被加工点的理论位置的偏离大小情况。

设与刀轴重合的方向矢量在铣削刀具坐标系中的表达式为:(rTx)Tx=[vTxx,vTxy,vTxz,0]

则其通过铣削分支在与床身固连的惯性坐标系O0-x0y0z0中的矩阵表达式为:

在工件坐标系中被加工点的理论刀具方向矢量为:

则其通过工件分支在与床身固连的惯性坐标系O0-x0y0z0中的矩阵表达式为:

可得铣削模式下刀具姿态误差为:

3 数控机床相对运动约束方程的建立

车削模式下,刀具中心位置相对运动约束方程为:

铣削模式下,刀具中心位置相对运动约束方程为:

铣削模式下,刀具方位相对运动约束方程为:

4 五轴数控机床几何误差分析与建模实例

以德国巨浪的DZ08FX双主轴高速五轴加工中心为例,以多体系统理论为基础,对机床的几何误差进行分析,建立机床各运动体的坐标系。利用各运动体的相对运动坐标变换矩阵分别推导出刀具位置和刀具姿态的理想成形约束方程和实际成形约束方程。

4.1 德国巨浪DZ08FX双主轴五轴数控机床几何误差分析

数控机床的误差源有很多,主要的误差源有几何误差、热误差、载荷误差和伺服系统误差等[8]。其中几何误差,受环境因素影响小,易于测量,且无论哪种误差,其最终表现形式都可以用前面叙述的几何误差分析与运动建模方法来表述。因此几何误差作为基本误差源,是本文研究的重点。五轴数控机床拥有多个运动部件,各个部件的误差可分为与位置点无关误差(静态误差)和与位置点相关误差(运动误差)两个部分[9]。

综上所述,五轴数控机床几何误差参数一共37项,如表1所示。

4.2 机床部件坐标系的建立

为了了解机床的运动,需要利用机床的各个运动部件相对运动坐标变换矩阵来描述其运动,将复杂的运动简化为数学模型。首先需要建立各部件的坐标系,部件坐标系有体坐标系和相对运动参考坐标系。

由前述可知,双主轴五轴数控机床有13个运动部件,首先设定初始条件下各运动体体坐标系和相对运动参考坐标系重合,即初始条件下确定运动体的参考坐标系方位就相当于确定了该运动体体坐标系方位。选取机床坐标系为基准,令床身和X导轨的运动参考坐标系方位和机床坐标系一致。令基准坐标系绕Y轴转过垂直度εx(z)后得到Z导轨运动参考坐标系。令Z导轨分别系绕Y轴、X轴转过垂直度εx(y)、εy(z)后得到Y导轨运动参考坐标系,令基准坐标系分别绕Y轴、Z轴旋转垂直度εAY、εAZ后得到A转台运动参考坐标系。令A转台运动参考坐标系分别绕X轴、Y轴旋转垂直度εCX、εCY后得到C转台的运动参考坐标系。这样就确定了初始条件下各坐标系的方向。各坐标系的位置由如下方法确定:令机床的各个运动部件返回到数控机床的绝对零点,令床身、X导轨、Y导轨、Z导轨和主轴体坐标系及其运动参考坐标系零点和机床坐标系零点重合,令A转台和C转台体坐标系以及运动参考坐标系原点位于A转台旋转中心。由上述方法可以确定双主轴五轴数控机床各运动部件的坐标系位置与方法,为后续的建模工作提供基础。

4.3 特征矩阵

根据上述设定方法便可得到双主轴数控机床各相邻体的理想特征矩阵以及误差特征矩阵。设工件坐标系相对于C轴体参考坐标系的位置阵列为(qwxqwyqwz1),A转台旋转中心在机床坐标系的齐次坐标为(q9xq9yq9z1)T,由此可得各运动部件特征矩阵:

4.4 加工模式一理想条件下机床运动学模型的建立

设刀具成形点在刀具坐标系内的齐次坐标为{rTx}Tx=(0 0–d 1)T,其中d为刀长,可以得到刀具成形点在刀具坐标系的理想成形运动方程:

将变换矩阵带入上述方程,经整理消去高阶无穷小量,可得:

实际刀轴矢量在刀具坐标系中的矢量表达式为{vTx}Tx=(0 0–1 0)T,可以得到刀轴矢量在刀具坐标系的理想运动表达式为:

将变换矩阵带入上述方程,经整理消去高阶无穷小量,可得:

4.5 加工模式一实际条件下机床运动学模型的建立

机床实际加工运动时候由于误差的存在使得实际位置偏移理想位置,因此要将实际的刀具成形点描述到工件坐标系中,和工件坐标系中理想的刀具成形点位置作比较。设工件坐标系下理论刀具成形点的齐次坐标为{rwx}wx=(xwx,ywy,zwz,1)T,实际情况下刀具成形点在工件坐标系运动方程为:

将变换矩阵带入上述方程,经整理消去高阶无穷小量,可得:

设工件坐标系下理论刀轴方向矢量的齐次坐标为{vwx}wx=(xwx,ywy,zwz,0)T,可以得到刀轴矢量在工件坐标系的实际运动表达式为:

将变换矩阵带入上述方程,经整理消去高阶无穷小量,可得:

在生产过程中,实际刀具中心点与理论刀具中心点要重合,即{PTx}0={Pwx}0刀具方向矢量也须满足{VTx}0={Vwx}0,即:

5 结论

1)针对双主轴五轴高速加工中心,对其进行结构分析,利用拓扑结构和低序体阵列对机床结构进行描述,将机床分为多个运动部件,并对每个运动部件建立体坐标系和运动参考坐标系。

2)对双主轴五轴高速加工中心的几何误差进行分析,得出机床37项几何误差。

3)利用多体系统运动学理论推导出机床各运动部件之间的相对运动坐标变换矩阵,进一步推导出机床理想条件下和实际条件下的运动学模型方程,并推出空间误差模型,为今后的误差辨识工作打下基础。

参考文献

[1]刘桂芝,周永良.影响车铣复合机床双主轴系统精度的因素分析[J].机械设计与制造工程,2014,(12):59-62.

[2]孙锡娜,韩秋实.车磨复合机床的发展现状及关键技术[J].精密制造与自动化,2008,(1):4-5.

[3]李德,立显凯.五轴车铣复合加工技术的现状与发展趋势[J].航空制造技术,2009,(12):47-50.

[4]吴宝海,严亚南,罗明,等.车铣复合加工的关键技术与应用前景[J].航空制造技术,2010,(19):42-45.

[5]李彦光.复合加工机床的发展现状[J].机械工程师,2011,(3):5-11.

[6]杨建国,潘志宏,薛秉源.数控机床几何和热误差综合的运动学建模[J].机械设计与制造,1998,(5):31-32.

[7]王晓峰.复合数控机床几何误差补偿及误差影响溯源分析[D].北京工业大学,2014.

[8]范晋伟.基于多体系统运动学的数控机床运动建模及软件误差补偿技术的研究[D].天津大学,1996.

柔性多体系统 篇6

针对网络通信拓扑切换,Belykh等提出耦合根据一定概率进行切换的闪烁网络[5],但其设定用于控制拓扑切换的切换信号通常过于迅速。对于现实中大部分多智能体系统,Liu等提出通信网络切换过于迅速易缩短机器的使用寿命[6]。Liberson等提出切换网络可能在速度较快的切换信号下无法保持稳定[7]。因此,针对一阶多智能体系统,研究人员开始研究通过限制切换信号转换的速度来减慢系统网络拓扑的切换速度。其中,Liu等研究拓扑切换下时延网络的局部指数同步[8],其通信网络拓扑根据ADT信号进行切换,但系统针对的是恒定耦合时延,在实际应用中具有较大的局限性。

基于以上论述的启发,本文主要考察在以下两个约束条件下一阶多智能体系统的一致性问题:(1)ADT转换规则下的受限有向切换信息交互网络;(2)时变的信息交互时延。

1 动力学模型

假定所研究的多智能体系统中包含智能体的个智能体,本文只考虑一阶积分器型时间连续系统,建立动力学模型如下[9]

其中,I={1,2,…,N}表示有限集;xi∈R表示第i个智能体的位置状态;xi·表示该位置状态的一阶导数;

表示一致性控制器[10];Niσ(t)表示在时间t上智能体i的邻接点集。同时,aijσ(t)取值可作如下规定

另外,d(t)表示时变时延,可定义为一种满足以下条件的时变可微函数

其中,h>0且μ为常数。同时,令x(t)=[x1(t),x2(t),…,xN(t)]T,则多智能系统可写成矩阵形式如下式所示

其中,Lσ为有限切换拓扑的拉普拉斯矩阵。在建立一阶带自时延和有限切换拓扑的多智能体系统模型后,本文的目标是分析多智能系统在具有ADT特性的切换拓扑下,达到指数一致性的充分条件。因此,需在此引入ADT的定义。

定义1对于切换信号σ(t)和任意T2≥T1≥0,令Nσ(T1,T2)表示σ(t)在区间(T1,T2)内切换的次数。若存在τa>0和整数N0≥0使得

则称τa为ADT。通常规定N0=0。

2 一致性分析

2.1 状态变换

为分析上述多智能体系统的指数一致性,本文对多智能体系统进行状态变换如下

其中,矩阵E=[-1N-1,IN-1]∈R(N-1)×N。对于i=1,2,…,N-1,有yi(t)=xi+1(t)-x1(t),则矩阵形式的多智能体系统可转化为如下降阶形式

其中,Πσ=-ELσF,F=[0N-1,IN-1]T∈RN×(N-1)。通过状态变换(5),向量y(t)∈R(N-1)可用于描述系统(1)中智能体间的差异。因此,经转换后的降阶系统(6)可用来描述原系统的动态差异,称其为差异系统。

2.2 指数一致性分析

本文采用分段Lyapunov-Krasovskii稳定性理论,来分析上述差异系统(6)的指数稳定性问题[11]。首先,考虑变换后系统中的第i个子系统如下所示

针对子系统(7),构造Lyapunov-Krasovskii泛函如下

其中

其中,Pi>0、Qi>0、Ri>0、Zi>0均为待定的正定矩阵,α>0、h>0均为给定标量。然后,分别对V1i(t)、V2i(t)进行关于时间t的一阶求导。

根据前文对时变时延d(t)所作的规定(2),对于Vi3(t)有

另外,对于Vi4(t),有

综合式(9)~式(12),并结合詹森不等式,可得

其中

假设Φi<0,则有

对其积分可得

针对整个系统(6),构造如下的分段LyapunovKrasovskii泛函

其中

假设存在β≥1,对于j∈M有

则根据式(16)进一步可得,对于切换时刻tj有

假设在时间区间(t0,t)内的切换次数其中t∈[tk,tk+1],结合式(14)和式(17)可得i,j∈M

同时,由式(15)可得

其中,同理可得

其中

结合式(18)~式(20)可得

即

其中,因此,多智能体系统(3)达到了指数一致性。给出总结性定理如下

定理1考虑以σ(t)作为切换信号的受限切换拓扑下的多智能体系统(1),其中σ(t)满足式(4),且时变时延d(t)满足式(2)。令α>0、h>0及μ为给定标量,若存在矩阵Pi>0、Qi>0、Ri>0、Zi>0,使得

其中,Φi11=Qi+Ri+αPi-e-αhZi,Φi12=PiΠi+e-αhZi,Φ22=(μ-e-αh)Ri+h2ΠiTZiΠi-2eαhZi,则当σ(t)带有满足时(其中,β≥1满足式(16)),控制器Ui能使系统达到指数一致性。

定理1给出在线性矩阵不等式框架下,检验受限切换拓扑下一阶多智能体系统是否达到指数一致性的充分条件。其中,LMI形式下的检验标准有以下优势:(1)检验过程中无需改变条件中的矩阵和(或)参数值;(2)通过LMI计算机数值模拟能使已执行问题得到有效验证。同时,在ADT满足式时,可进一步确定切换信号σ(t)使得多智能体系统达到指数一致性。

柔性多体系统 篇7

多智能体系统的分布式协调控制在许多领域有着广泛的应用, 如分布式传感器网络、无人飞行器、蜂拥、多自主体的群集、通信网络的拥堵控制、卫星群的队形控制等。这里分布式协调指的是仅通过智能体之间的局部信息交流以实现多智能体系统全局协调行为的涌现。因此, 该领域的一个基本问题在于设计合适的控制协议, 使得多智能体能够在局部信息通信下达到一致。

近年来, 有关多智能体系统一致性问题的研究引起来自各领域学者的广泛关注, 他们提出了众多新颖的低阶多智能体的控制协议。Olfati-Saber和Murray于2004年首先建立一阶多智能体系统一致性问题的理论框架[1]。Ren和Beard进一步考虑了基于有向切换拓扑的的信息一致性的问题, 同时给出了一些较宽松条件[2]。对于有向固定拓扑通信拓扑和有向切换通信拓扑的多智能体系统, 文献[3]研究了它们在离散时间描述状态下的一致性算法。文献[4]则探讨了一种离散时间的多智能体的新的集合算法。通过选取不同的一致性冗余, 提出了满足线性、周期的、正指数动力学等要求的二阶多智能体系统的协议。最近, 一致性问题的研究热点正向高阶或者非线性动力学的方向发展。

由于多智能体系统的通信能力有限, 比如多智能体系统经常受到信息传输速度、信息传输拥堵以及信息不对称性等因素的影响, 在实际多智能体系统应用中时延是不可避免的。一般来说, 时延可分为输入时延和通信时延, 输入时延指多智能体系统对外部作用和信号的连接或处理时延, 通信时延指多智能体之间的信息传输通过传感器或通信网络交流的传输时延。有关时延多智能体系统的一致性的研究成果, 感兴趣的读者可参考文献[5]。

在本文中, 我们主要研究在有向信息通讯模式下, 二阶时延多智能体系统的一致性问题。假设所考虑的多智能体系统中输入时延较小, 因此我们对其忽略不计而只考虑通信时延。此外, 领导者的状态可以是时变的, 并且仅有部分跟随者可以直接接收来自领导者的信息。通过构造一个新颖的Lyapunov-Krasovskii泛函, 我们获得了一个充分条件使得所有智能体最终趋于领导者状态, 即达到多智能体系统的一致性。值得指出的是, 所获得的充分条件可以表示成线性矩阵不等式的框架, 能方便借助Matlab的LMI工具箱直接算出结果。

2问题描述

2.1 图论

多智能体系统的通信网络拓扑通常用有向图来描述。设G= (V, E, A) 为一个含n个节点的加权有向图, 其中V={v1, v2, …, vn}是一个有限非空的节点集, E⊆V×V为边集, A={aij}∈Rn×n (aij≥0) 为加权邻接矩阵, 节点i代表第i个智能体。在G中从节点i到节点j的有向边eij= (i, j) , eij∈E当且仅当智能体vj能接收到vi的信息。节点vi的邻居节点集表示为Ni={vj∈V|eji∈E}。

邻接矩阵A中的邻接元素aij是正的, 即aij>0⇔eji∈E, aij=0⇔eji∉E。此外, 我们假定aii=0, 也就是说图G不含自环环。undefined称作节点vi的入度, undefined称作节点vi的出度。有向图G的拉普拉斯矩阵L=[lij]∈Rn×n定义为

undefined

这里L不一定是对称的。此外, 拉普拉斯矩阵L的一个重要性质是1n=[1, 1, …, 1]T∈Rn是L的一个右特征向量, 对应的特征值为λ=0。

如果有向图中存在有向路径使得每个节点可达其他任意节点, 则称这个有向图是强连通的。若只有一个特殊的节点可达有向图中的任意节点, 则称该节点为根节点。一个有向图的生成树是图中连接所有节点的边构成的有向树。如果至少存在一个节点, 即根节点, 可达所有其他节点, 我们就称该图有一个生成树。

2.2 一致性算法

假定每个智能体都满足如下二阶动态系统

这里xi∈Rn和vi∈Rn表示智能体i的位置和速度, ui∈Rn表示智能体的控制输入 (或者加速度) 。当系统 (1) 的智能体之间信息传输存在定常通信时延d时, 采用如下一致性算法:

undefined

其中Ni表示节点vi的邻居节点集, τ表示通信时延常数, vd表示领导者的速度, ki表示控制参数 (如果领导者的信息可达智能体vi, 则ki>0, 否则ki=0) 。

注1:注意到在 (2) 中ui是智能体vi的本地控制器, 只依赖于其邻居及领导者的信息。对于多智能体系统 (2) , 若undefined, 则称多智能体系统 (2) 渐近达到一致。

注2:为了简单起见, 在 (1) 中每个智能体只考虑了一维的情况。实际上, 通过引入Kronecker积, 我们能够很容易得到高维的结果。

3主要结果

在固定拓扑G下考虑带领导者的多智能体系统 (2) , 我们预先设定通信时延, 构造合适的Lyapunov-Krasovskii泛函。下面的定理以线性矩阵不等式的方法说明了相关结果。

定理1:假定通信时延τ定常, 如果存在矩阵P>0和Q>0使得以下线性矩阵不等式成立:

则二阶多智能体系统 (1) 在协议 (2) 下渐近达到一致性。其中:

Ω11=PW+WTP+τWTQW-τ-1Q, Ω12=PH+τWTQH+τ-1Q,

Ω21=HTP+τHTQW+τ-1Q, Ω22=τHTOH-τ-1Q;

undefined

证明:根据一致性协议 (2) , 多智能体系统 (1) 的闭环形式为

undefined

令undefined系统误差向量, 则我们可以得到系统 (1) 的误差动力学如下:

undefined

其中

undefined

{k1, k2, …, kn}。

对系统 (5) 构造Lyapunov-Krasovskii泛函:

V (t) =V1+V2 (6)

其中undefined。

记

undefined

, 于是

undefined

另外, 根据Schur Complement引理[6]可得

undefined

此时, 记ζ (t) =[δT (t) , δT (t-τ) ]T, 则我们有

undefined

4结论

本文研究了含有通讯时延的二阶多智能体系统的领导者-跟随一致性问题。尽管领导者的状态是时变的, 我们提出的基于邻居局部信息的通讯协议使得每个智能体都能跟随领导者状态, 最终达到渐近一致。目前我们仅考虑了定常时滞的情形, 接下来我们将进一步深入研究多变时滞、异质时滞以及混合时滞情形下的多智能体系统的一致性问题。

摘要：本文研究了在有向信息通讯模式下, 时延二阶多智能体系统的领导者-跟随一致性问题。我们仅要求二阶多智能体系统的拓扑结构图中含有一个有向生成树, 而且仅有部分跟随者能够直接接收领导者的信息。通过设计一个新颖的Laypunov-Krasovskii泛函, 我们获得了在定常通信时延下使得智能体系统达到一致性的充分条件, 该条件可以表达为一个线性矩阵不等式结构, 很方便通过计算机数值求解。

柔性多体系统 篇8

关键词：多体系统理论,凸轮轴磨削,几何误差建模,几何误差参数辨识

0 引言

20世纪90年代以来,国内汽车工业、内燃机工业迅速发展,作为汽车、摩托车、内燃机的关键零件,凸轮轴的需求量越来越大。由于整个进排气系统是由凸轮轴驱动的,故凸轮的形线精度对发动机性能的影响是极其关键的,凸轮轴的加工质量和加工效率将直接影响到我国汽车产品的质量和汽车工业的发展。

磨削加工是保证产品形位精度的最后一道工序,其加工精度将直接影响产品的质量。然而,由于数控凸轮轴磨床自身存在各种误差因素的影响,导致凸轮廓形的加工精度极难保证。在这些误差因素中几何误差对凸轮轮廓精度的影响最大。因此,如何克服磨床本身的几何误差,提高其加工精度是各国学者和企业普遍关注的问题。误差补偿法[1,2]克服了传统误差防止法造价昂贵、适用性差的缺点,作为提高加工精度的有效手段,得以迅速发展。误差补偿技术是通过分析数控凸轮轴磨床影响加工精度的几何误差来源,建立几何误差数学模型,经过对磨床的几何误差进行辨识和补偿,从而提高加工精度的方法。本文重点研究的是凸轮轴磨削几何误差的建模和辨识方法。

1 建立数控凸轮轴磨床精密加工运动约束条件方程

近年发展起来的多体系统理论的误差建模方法[3,4,5]全面考虑了影响机床精度的各项因素,以特有的低序体阵列来描述复杂系统,具有建模过程程式化、约束条件少、易于解决多体系统运动问题的特点。本文采用这一理论,提出了一种广泛适用的、易于实现计算机自动编程的数控凸轮轴磨床几何误差模型的建模方法,使得建模过程程式化,排除了人为因素对模型推导过程的影响。

应用误差补偿技术提高数控凸轮轴磨床的几何加工精度,关键在于数控系统对实际机床发出经过修正的精密数控指令,即在考虑数控机床运动误差的情况下,得出修正的数控指令来控制刀具精确运动。为建立理想条件和实际条件下的数控凸轮轴磨床运动约束条件方程[6,7],本文根据MKS8332A数控凸轮轴磨床的结构特点,在机床床身导轨、砂轮中心及凸轮轴工件上分别建立了多个动静参考坐标系。图1所示为MKS8332A数控凸轮轴磨床结构及各部件坐标参量,PT为实际砂轮中心在机床床身坐标系内的位置矢量,PW为某个凸轮中心坐标系内给定的理想砂轮中心在机床床身坐标系内的位置矢量。

1.砂轮滑台 2.工件滑台 3.头架顶尖 O—机床床身 T—砂轮 W—凸轮轴

要使数控凸轮轴磨床能实现精密磨削,就必须保证在任意瞬时实际砂轮中心位置与理论砂轮路线的对应点完全一致,即PW=PT,这便是数控机床精密加工条件约束方程。即

式中,rW为理论砂轮中心在工件坐标系中位置列阵;rT为刀具中心在刀具坐标系内位置列阵;[AOW]为误差情况下凸轮轴体坐标系相对于床身坐标系的实际变换矩阵,A表示矩阵,O表示床身,W表示凸轮轴,其他依此类推;带下标p、pe、s、se的矩阵分别表示体间理想静止特征矩阵、体间静止误差特征矩阵、体间理想运动特征矩阵、体间运动误差特征矩阵。

将参数具体化,便可得到该精密加工条件约束方程。表1为数控凸轮轴磨床各相邻体间的变换矩阵列表。其中,x为砂轮滑台位移,δx(x)、δy(x)、δz(x)、εx(x)、εy(x)、εz(x)分别为砂轮滑台在运动过程中三个线位移和三个角位移误差;εxz为X轴与Z轴垂直度误差;θ为头架顶尖回转角;Δαyc、Δβxc为C轴与X轴、Y轴垂直度误差;Δα(θ)、Δβ(θ)、Δγ(θ)、Δx(θ)、Δy(θ)、Δz(θ)为主轴回转误差和主轴径向跳动误差;Z为凸轮坐标系与凸轮左端面工件坐标系间距离;ϕ在下文中会有详细介绍。在实际加工过程中,对于凸轮轴上每一个凸轮来说,工件滑台部件的位置一旦调整确定,就不再做任何运动,对单个凸轮廓形的加工精度几乎没有影响,因而不考虑工件滑台部件的运动误差。

从精密加工约束条件方程可以看出,只要通过砂轮滑台沿X方向的适当移动和头架顶尖沿C轴的适当转动,使得工件在磨床床身坐标系中位置和刀具在磨床床身坐标系中位置始终严格保持一致,就可使凸轮廓形的几何加工误差趋于零。也就是说提高凸轮加工精度的问题就转化为在考虑磨床几何误差因素影响下,给定什么样的数控驱动指令值,使得凸轮廓形的几何加工误差趋于零的问题。

2 凸轮轴上各凸轮坐标系描述方法

加工凸轮轴工件时,以工件左端面上中心孔和定位孔的连线作为工件的X轴。为了方便描述凸轮轴上每个凸轮轮廓,需要在每个凸轮上都建立一个凸轮坐标系。如图2所示,设头架顶尖坐标系绕Y轴顺时针转180°后的坐标系与凸轮左端面工件坐标系的夹角为φ,凸轮左端面工件坐标系与第一个凸轮桃尖夹角为θ1,则第一个凸轮坐标系与凸轮左端面工件坐标系夹角φ1=φ±(90°-θ1),即表1中[AW3W]p中的ϕ可表示为:ϕ=φ±(90°-θ1),当90°-θ1>0时取正号,否则取负号。若第一个凸轮桃尖与第二个凸轮桃尖夹角为θ2,则第二个凸轮坐标系与凸轮左端面工件坐标系的夹角为φ2=φ±[90°-(θ1+θ2)],即表1中[AW3W]p中的ϕ可表示为:ϕ=φ±[90°-(θ1+θ2)]。依次类推我们可以得到凸轮轴上不同方位的凸轮坐标系,即不同方位凸轮精密加工条件约束方程。

3 数控凸轮轴磨床几何误差参数辨识

要想通过精密加工约束方程准确计算该磨床各运动部件的运动情况,以达到误差补偿的目的,还必须建立一套准确的数控机床误差参数辨识方法,以便获得准确的运动误差参数。对于机床直线导轨6项误差,一般采用激光干涉仪用9线法[8]进行辨识;对于回转轴的6项误差,通常采用标准芯棒和五传感器的方法[9]进行测量,单纯依靠某种辨识方法是很难对磨床所有的误差进行辨识的,若将几种辨识方法结合起来则需要的检测设备又太多。本文利用多体系统理论的几何误差参数辨识模型[10,11,12],结合球杆仪测量原理[13,14,15]提出了可以辨识所有几何误差的方法。

从该磨床运动模型的建立过程可以看出,影响凸轮廓形的几何误差有15项,这些误差参数的辨识方法在下文进行详细描述。

3.1 C轴运动5项几何误差参数辨识方法

在头架顶尖坐标系O3中取一点P(xp0,yp0,zp0)。设ΔL=(Δx(θ),Δy(θ),Δz(θ),0)T表示点P在头架绕C轴转动θ角后,由于运动误差的存在而产生的实际运动点和理想运动点的差值。则:

$\begin{array}{l}\text{Δ}\mathit{L}=\left[\begin{array}{c}{x}^{\prime }{}_p\\ y^{\prime }{}_p\\ z^{\prime }{}_p\\ 1\end{array}\right]-\left[\begin{array}{c}x{}_p\\ y{}_p\\ z{}_p\\ 1\end{array}\right]=\left[\begin{array}{cccc}\cos \theta & -\sin \theta & 0& 0\\ \sin \theta & \cos \theta & 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\cdot \\ \left[\begin{array}{cccc}1& -\text{Δ}\gamma {}_{(\theta )}& \text{Δ}\beta {}_{(\theta )}& \text{Δ}x{}_{(\theta )}\\ \text{Δ}\gamma {}_{(\theta )}& 1& -\text{Δ}\alpha {}_{(\theta )}& \text{Δ}y{}_{(\theta )}\\ -\beta {}_{(\theta )}& \text{Δ}\alpha {}_{(\theta )}& 1& \text{Δ}z{}_{(\theta )}\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}x{}_{p0}\\ y{}_{p0}\\ z{}_{p0}\\ 1\end{array}\right]-\\ \left[\begin{array}{cccc}\cos \theta & -\sin \theta & 0& 0\\ \sin \theta & \cos \theta & 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}x{}_{p0}\\ y{}_{p0}\\ z{}_{p0}\\ 1\end{array}\right](1)\end{array}$

如图3所示,球杆仪底座固定在C轴头架上,另一端连于尾座顶尖处。建立磨床误差模型时,我们只考虑了头架安装误差和旋转运动时产生的运动误差,忽略了尾座顶尖与头架间的安装误差,而假设两者位于一条直线上,因为一旦两者轴心不在一条直线上,即使在无误差情况下也会产生杆长变化。

球杆仪的两球设为球A、球B,其中A球固接在C轴上随C轴一起转动,A球球心在头架顶尖坐标系O3中的坐标为(x,y,z),C轴转动θ角后,A球球心在头架顶尖坐标系中的坐标为(x(θ),y(θ),z(θ)),将B球固定在尾座顶尖处,B球的球心坐标为(0,0,z1),设杆长为L,则有

L2=x2(θ)+y2(θ)+[z(θ)-z1]2

求导得

LΔL(θ)=x(θ)Δx(θ)+y(θ)Δy(θ)+

(z(θ)-z1)Δz(θ) (2)

因此可在C轴上选取6个测量点,测量出这6个点从0°～360°旋转不同的角度所对应的杆长变化量,然后联立式(1)、式(2)便可求出C轴旋转时产生的6项旋转运动误差。通过分析,我们知道Δγ(θ)并不引起球杆仪杆长的变化,所以上述方程组得到的Δγ(θ)数值为不可信值,但在C-X轴联动误差辨识模型中我们可以将其求出。

但在实际情况下,头架顶尖和尾座的轴心是不在一条直线上的。如图3所示,设B球在X、Y方向的偏心量为m、n,即B球球心的坐标为(m,n,z1)。同理可得

L2=(x(θ)-m)2+(y(θ)-n)2+(z(θ)-z1)2

求导得

LΔL(θ)=(x(θ)-m)Δx(θ)+

(y(θ)-n)Δy(θ)+(z(θ)-z1)Δz(θ) (3)

由于引入了两个未知量,只要增加两个测点,即再增加两组方程,仍可求得C轴旋转时产生的旋转运动误差。此时,式(3)中的ΔL(θ)不再是球杆仪的杆长变化量,它必须减去头架和尾座两者轴心偏心所引起的变化d:

$\begin{array}{l}d=\sqrt{(x(\theta )-m){}^2+(y(\theta )-n){}^2+(z(\theta )-z{}_1){}^2}-\\ \sqrt{(x-m){}^2+(y-n){}^2+(z-z{}_1){}^2}(4)\end{array}$

此时有ΔL(θ)=ΔL-d。其中,ΔL为球杆仪的杆长变化量。

3.2 X轴运动几何误差参数辨识方法

在砂轮滑台上取一点P(xp0,yp0,zp0)。设ΔL=(Δx(x),Δy(x),Δz(x),0)T表示点P在砂轮滑台移动x距离后,由于运动误差的存在,而产生的实际运动点和理想运动点的差值。则:

$\text{Δ}\mathit{L}=\left[\begin{array}{c}{x}^{\prime }{}_p\\ y^{\prime }{}_p\\ z^{\prime }{}_p\\ 1\end{array}\right]-\left[\begin{array}{c}x{}_p\\ y{}_p\\ z{}_p\\ 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& x\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$

·

$\begin{array}{l}\left[\begin{array}{cccc}1& -\epsilon {}_z(x)& \epsilon {}_y(x)& -\delta {}_x(x)\\ \epsilon {}_z(x)& 1& -\epsilon {}_x(x)& -\delta {}_y(x)\\ -\epsilon {}_y(x)& \epsilon {}_x(x)& 1& -\delta {}_z(x)\\ 0& 0& 0& 1\end{array}\right]\cdot \\ \left[\begin{array}{c}x{}_{p0}\\ y{}_{p0}\\ z{}_{p0}\\ 1\end{array}\right]-\left[\begin{array}{cccc}1& 0& 0& x\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}x{}_{p0}\\ y{}_{p0}\\ z{}_{p0}\\ 1\end{array}\right](5)\end{array}$

如图4所示,球杆仪两球分别设为A球、B球,将B球固定在磨床床身某点处,其在砂轮滑台中坐标系设为(x0,y0,z0),把A球放在砂轮滑台上随砂轮滑台一起移动(砂轮滑台垂直纸面方向运动),A球放置时初始位置与B球的连线和X轴垂直,A球球心在砂轮滑台坐标系O1中的坐标为(x0,y1,z1),砂轮滑台移动距离x之后,A球球心在砂轮滑台坐标系中坐标为(x(x),y(x),z(x))。图5为X轴运动示意图,C点为A球理想位置,D点为A球实际位置。设砂轮滑台移动距离x后杆长为L,则有

L2=(x(x)-x0)2+(y(x)-y0)2+(z(x)-z0)2

求导得

LΔL(x)=(x(x)-x0)Δx(x)+

(y(x)-y0)Δy(x)+(z(x)-z0)Δz(x) (6)

在没有运动误差情况下:

$L{}_0^2=x{}^2+[\sqrt{(y{}_1-y{}_0){}^2+(z{}_1-z{}_0){}^2}]{}^2$(7)

所以

ΔL(x)=L-L0 (8)

因此,可在X轴上选取6个不同测量点,测出这6个点在X向导轨运动不同位置时对应的杆长变化量,然后联立式(5)～式(8)便可求出X轴运动时产生的6项运动误差。

3.3 C-X运动轴联动时部分几何误差参数辨识方法

当凸轮进入加工区域时,其Z向工件滑台导轨不再移动。数控凸轮轴磨床的砂轮滑台沿X向作往复进给运动的同时,磨床主轴沿C轴作变速的旋转运动,这便是C-X运动轴联动情况。如图6所示,球杆仪底座固定在C轴头架上,另一端固定在砂轮滑台上。设BK体为头架,BH体为砂轮滑台,BJ体为机床床身,令rK表示pK点相对于BK体坐标系位置列阵,PKh表示pK点相对于BJ体坐标系位置列阵,rH表示pH点相对于BH体坐标系位置列阵,

PH θ表示pH点相对于BJ体坐标系位置列阵。

BK体上pK点位置方程为

$\left[\begin{array}{c}\mathit{{\rm P}}{}_{\text{Κ}h}\\ 1\end{array}\right]=[A{\rm O}{\rm T}{}_1]\left[\begin{array}{c}\mathit{r}{}_{\text{Κ}}\\ 1\end{array}\right]=$

$[A{\rm O}{\rm T}{}_1]{}_{\text{p}}[A{\rm O}{\rm T}{}_1]{}_{\text{p}\text{e}}[A{\rm O}{\rm T}{}_1]{}_{\text{s}}[A{\rm O}{\rm T}{}_1]{}_{\text{s}\text{e}}\left[\begin{array}{c}\mathit{r}{}_{\text{Κ}}\\ 1\end{array}\right]$

BH体上pH点位置方程为

$\left[\begin{array}{c}\mathit{{\rm P}}{}_{\text{Η}\theta }\\ 1\end{array}\right]=[A{\rm O}W{}_2][AW{}_{2}W{}_3]\left[\begin{array}{c}\mathit{r}{}_{\text{Η}}\\ 1\end{array}\right]=[A{\rm O}W{}_2]{}_{\text{p}}$

·

[AOW2]pe[AOW2]s[AOW2]se[AW2W3]p·

$[AW{}_{2}W{}_3]{}_{\text{p}\text{e}}[AW{}_{2}W{}_3]{}_{\text{s}}[AW{}_{2}W{}_3]{}_{\text{s}\text{e}}\left[\begin{array}{c}\mathit{r}{}_{\text{Η}}\\ 1\end{array}\right]$

假设在初始位置(砂轮滑台在机床床身上的位移量h=0,θ=0)时,各运动部件运动参数为零,且各运动件的所有几何误差参数也为零。当C-X轴联动时,可得在BJ体坐标系内机床砂轮滑台上给定点P1i相对于机床头架顶尖上给定点P3i的相对位移方程:

$\left[\begin{array}{c}{\rm P}{}_{1ih}\\ 1\end{array}\right]-\left[\begin{array}{c}{\rm P}{}_{3i\theta }\\ 1\end{array}\right]-(\left[\begin{array}{c}{\rm P}{}_{1i0}\\ 1\end{array}\right]-\left[\begin{array}{c}{\rm P}{}_{3i0}\\ 1\end{array}\right])=$

$\left[\begin{array}{l}x{}_h+\delta {}_x(x{}_h)-\epsilon {}_z(x{}_h)r{}_{1iy}+\epsilon {}_y(x{}_h)r{}_{1iz}+r{}_{3ix}+\\ \epsilon {}_{xz}r{}_{1iz}+\cos \theta (\text{Δ}\gamma {}_{(\theta )}r{}_{3iy}-r{}_{3ix}-\text{Δ}\beta {}_{(\theta )}r{}_{3iz}-\\ \text{Δ}x{}_{(\theta )})+\sin \theta (\text{Δ}\gamma {}_{(\theta )}r{}_{3ix}+r{}_{3iy}-\text{Δ}\alpha {}_{(\theta )}r{}_{3iz}+\text{Δ}y{}_{(\theta )})\\ \delta {}_y(x{}_h)-\epsilon {}_x(x{}_h)r{}_{1iz}+\epsilon {}_z(x{}_h)r{}_{1ix}+r{}_{3iy}+\\ \cos \theta (-\text{Δ}\gamma {}_{(\theta )}r{}_{3ix}-r{}_{3iy}+\text{Δ}\alpha {}_{(\theta )}r{}_{3iz}-\text{Δ}y{}_{(\theta )})+\\ \sin \theta (\text{Δ}\gamma {}_{(\theta )}r{}_{3iy}-r{}_{3ix}-\text{Δ}\beta {}_{(\theta )}r{}_{3iz}-\text{Δ}x{}_{(\theta )})\\ \delta {}_z(x{}_h)+\epsilon {}_x(x{}_h)r{}_{1iy}-\epsilon {}_{xz}x{}_h-\epsilon {}_y(x{}_h)r{}_{1ix}-\\ \epsilon {}_{xz}r{}_{1ix}+\text{Δ}\beta {}_{(\theta )}r{}_{3ix}-\text{Δ}\beta {}_{xc}r{}_{3ix}-\text{Δ}\alpha {}_{(\theta )}r{}_{3iy}+\\ \text{Δ}\alpha {}_{yc}r{}_{3iy}-\text{Δ}z{}_{(\theta )}+\cos \theta (\text{Δ}\beta {}_{xc}r{}_{3ix}-\text{Δ}\alpha {}_{yc}r{}_{3iy})+\\ \sin \theta (-\text{Δ}\alpha {}_{yc}r{}_{3ix}-\text{Δ}\beta {}_{xc}r{}_{3iy})\end{array}\right]$

若要转角θ及进给量x使球杆仪杆长L保持不变,则上式各项都为零。从上述方程可以看出,第二个方程只有Δγ(θ)是未知量,依照上述方法在砂轮滑台导轨上和头架顶尖上分别确定一点即可求出Δγ(θ);第一个方程只有εxz和Δγ(θ)是未知量,由于Δγ(θ)已经解出,按照上述方法可以很容易地解出εxz;第三个方程中Δβxc、Δαyc是未知量,在此可将其近似看作相等,即可求出这两项误差。然而,在C-X运动轴联动情况下,当C轴转动一圈时,X轴来回作两次往复直线运动,X轴的位移和C轴转动90°所对应的弧长是相等的,此时X轴两端有极小部分无法运动到,由于X轴在加工过程中往复运动的位移很小,X轴所运动的范围完全可以满足该轴几何误差辨识范围的需要。

4 精密数控加工指令值的迭代逆解方法

4.1 构造迭代求解方程

根据精密加工条件约束方程,无法直接对精密加工数控指令x、θ进行求解。需将方程中与磨床运动量x、θ有关的误差参数用某一个已知的运动量即xi、θi所属的误差参数来替代,这样便可解得更加精密的运动量xi+1、θi+1。可推得:

$\sin {}^2\theta {}_{i+1}=(-F+\sqrt{F{}^2-4DE})/(2D)$ (9)

xi+1=z(-Δβxc-cosθΔβ(θ)-sinθΔα(θ))+cosθΔx(θ)-

sinθΔy(θ)-δx(xi)+rW x[cosϕ(-cosθ+sinθΔγ(θ))-

sinϕ(sinθ+cosθΔγ(θ))]+rW y[-cosϕ(sinθ+

cosθΔγ(θ))-sinϕ(sinθΔγ(θ)-cosθ)] (10)

D=(B2+C2)2E=(A2-C2)2

F=2(A2C2-A2B2-B2C2-C4)

A=[δy(xi)-zΔαyc]

B=sinθi[rW x(-cosϕ-sinϕΔγ(θi))+

rW x(sinϕ-cosϕΔγ(θi))+Δx(θi)-zΔβ(θi)]

C=cosθi[rW x(-cosϕΔγ(θi)+sinϕ)+

rW y(sinϕΔγ(θi)+cosϕ)+Δy(θi)+zΔα(θi)]

式(9)和式(10)即为精密数控加工指令值的迭代求解方程,其中θi+1的具体取值根据砂轮路线列表点所在象限来决定。

4.2 确定迭代初值

求解迭代方程,必须先要给定一个迭代初值。本文以理想数控指令作为迭代初值,即假设该数控凸轮轴磨床不存在运动误差,要加工凸轮廓形上某一点,求解对应的理想数控指令值,只要令精密加工约束条件方程中所有的运动误差参数为零,就可得理想加工约束条件方程,即可确定理想数控指令迭代初值:

$\left[\begin{array}{cccc}\cos \theta & -\sin \theta & 0& 0\\ \sin \theta & \cos \theta & 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}-\cos \text{ϕ}& \sin \text{ϕ}& 0& 0\\ \sin \text{ϕ}& \cos \text{ϕ}& 0& 0\\ 0& 0& -1& z\\ 0& 0& 0& 1\end{array}\right]$

·

$\left[\begin{array}{c}r{}_{\text{W}x}\\ r{}_{\text{W}y}\\ 0\\ 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& x\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right]$

由上式可得理想数控指令值:

x0=(-rWxsinϕ-rWycosϕ)sinθ+

(-rWxcosϕ+rWysinϕ)cosθ

$\tan \theta {}_0=\frac{r{}_{\text{W}x}\sin \text{ϕ}+r{}_{\text{W}y}\cos \text{ϕ}}{r{}_{\text{W}x}\cos \text{ϕ}-r{}_{\text{W}y}\sin \text{ϕ}}$

4.3 迭代求解的中止判别条件

为使精密数控加工指令值的迭代求解过程能够在达到满意的程度后停止继续执行,本文提出的迭代求解中止条件是:如果用本次迭代求解所得的精密加工指令值驱动数控磨床后,实际砂轮中心位置rT与理论砂轮中心位置rW间的误差e′小于该数控磨床的运动分辨率,则迭代求解过程停止。

根据数控凸轮轴磨床精密加工条件方程,实际砂轮刀具中心位置rT为

$\left[\begin{array}{c}r{}_{\text{Τ}x}\\ r{}_{\text{Τ}y}\\ 0\\ 1\end{array}\right]=\mathit{D}{}^{-1}\mathit{C}{}^{-1}\mathit{B}{}^{-1}\mathit{A}{}^{-1}\left[\begin{array}{c}x+\delta {}_x(x)+\epsilon {}_{xz}\delta {}_z(x)\\ \delta {}_y(x)\\ \delta {}_z(x)-\epsilon {}_{xz}\delta {}_x(x)-\epsilon {}_{xz}x\\ 1\end{array}\right]$

A=[AW2W3]peB=[AW2W3]s

C=[AW2W3]seD=[AW3W]p

因此该磨床精密数控指令迭代求解中止判别条件就可表示为

$\left[\begin{array}{c}{e}^{\prime }{}_x\\ e^{\prime }{}_y\\ 0\\ 1\end{array}\right]=\left[\begin{array}{c}r{}_{\text{Τ}x}\\ r{}_{\text{Τ}y}\\ 0\\ 1\end{array}\right]-\left[\begin{array}{c}r{}_{\text{W}x}\\ r{}_{\text{W}y}\\ 0\\ 1\end{array}\right]\le \left[\begin{array}{c}\text{Δ}x\\ \text{Δ}y\\ 0\\ 1\end{array}\right]$

式中,Δx、Δy由数控凸轮轴磨床在X轴和C轴上的运动分辨率换算得出。

5 刀具路线的计算方法与曲线的拟合方法

本研究中,凸轮轮廓线偏置的计算方法即为理想刀具路线的计算方法。首先将各类凸轮廓形统一用点列表加以描述,再用三次样条曲线来拟合所有的离散点,然后针对每个凸轮廓形列表点,计算该列表点处的圆弧法线方程,判别法线方向,凹/凸圆弧沿内外法线方向偏移一个砂轮刀具半径,获得该点所对应的砂轮中心路线上的一个列表点。以此方法,依次获得磨削凸轮所需的砂轮刀具路线上的全部列表点。对于新的列表点我们可以再次用三次样条曲线来拟合得到凸轮的等距包络线方程,即砂轮磨削路线。通过砂轮磨削路线,我们可以获得理想的数控加工指令代码,这些数据是误差补偿所必需的基础数据。

本文分别建立了数控凸轮轴磨床的几何误差模型、精密加工约束条件方程和各项误差参数的辨识方法,为精密加工指令值的求解奠定了理论基础。对于加工一个特定的凸轮所需的理想砂轮中心路线上的每一个列表点,按照上述方法,依次求得其所对应的精密数控指令值之后,我们就可以把这一系列的精密数控指令值写成圆弧格式的G代码形式,输入数控系统便可以直接驱动数控凸轮轴磨床,完成精密凸轮磨削加工过程。因此,提高凸轮加工精度的问题就转化为在考虑磨床运动误差因素影响下,给定什么样的数控驱动指令值,能使凸轮廓形的几何加工误差趋于零的问题。

6 结论

(1)基于多体系统理论,能够准确、方便地建立数控凸轮轴磨床几何误差模型。利用这一模型可以推导出该磨床精密加工条件约束方程。

(2)计及磨床自身结构误差和运动误差的精密加工条件方程,隐含了数控指令、刀具路线和刀具轨迹之间的相互映射关系。该方程既可以用于求解精密的数控加工指令值,又可以用于求解真实的刀具轨迹,两者是同一问题的正解和逆解。

(3)利用多体系统理论的几何误差参数辨识模型结合球杆仪测量原理所提出的辨识方法,能够很好地解决该磨床的几何误差参数辨识问题,但在提高方程解的精度方面还需做进一步研究。