面向实时控制应用的过约束重型并联机床动力学建模方法

摘要

本发明公开了一种面向实时控制应用的过约束重型并联机床动力学建模方法，属于机床制造技术领域；该方法包括，根据机床的结构求出各个主动关节的位置，速度和加速度；以及各构件质心的位置，速度和加速度，建立动平台和各个主动关节的速度关系；采用牛顿-欧拉方法建立过约束重型并联机床中各构件的力平衡方程和力矩方程，根据过约束重型并联机床的结构参数和支链的杆件轴向变形情况，建立其变形与动平台输出误差之间的协调方程，最后，结合力平衡方程、力矩方程以及协调方程，建立动力学模型。本发明考虑重型并联机床中刚度较差连杆的变形，不仅克服了传统的刚体动力学模型精度低的问题，而且建立的动力学模型满足控制系统的实时计算要求。

说明书

技术领域

本发明属于机床制造技术领域，尤其涉及于过约束重型并联机床动力学建模法。

背景技术

并联机床在理论上具刚度质量比大、响应速度快和环境适应性强等优势，然而也存在 工作空间小、奇异位形多等缺点。驱动冗余是解决这一问题的有效方法，并且驱动冗余可 以进行并联机床的自标定。此外，通过引入额外连杆，将并联机床设计为过约束机构，可 以提高机床的刚度。因此，驱动冗余的过约束并联机床具有较好的静、动态性能，是一类 具有广阔应用前景的机床。为了实现高动态性能，在控制系统以及结构设计中，需要高精 度动力学模型，特别是基于动力学模型的高速高精度并联机床的控制系统中对动力学模型 的精度和实时性要求更严格。

一般常用的刚体动力学建模方法，如牛顿－欧拉方法法、达朗伯－拉格朗日方程法 和虚功原理法，可以用来建立过约束和驱动冗余并联机床的刚体动力学模型，但是，驱动 冗余并联机床刚体动力学逆问题的解不唯一。例如，Nahon和Angeles概述了针对驱动冗余 刚体动力学逆问题的四种方法：正交分解法、伪逆法、拉格朗日法、直接法。对于过约束 并联机构动力学模型，其静态方程个数少于未知力的个数，不能直接求解。为了求解所有 未知力，在建模中通常将过约束并联机构简化为非过约束机构，此时，静态力方程的个数 等于未知力的个数。尽管这些方法可以建立过约束及驱动冗余并联机床的动力学模型，但 是建模中没有考虑到刚度较差杆件的变形以及将过约束机构简化为非过约束机构带来的 误差，所建立的动力学模型不准确。尤其是过约束重型并联机床中连杆的变形更大，如果 不考虑连杆变形，所建立的刚体动力学模型精度较低。一种过约束重型并联机床的结构如 图1所示,该机床的并联机构部分由三个运动学支链组成，其中两个支链1为平行四边形结 构，分别由滑块B₁B₂，B₃B₄、转动副（又称为关节）、定长连杆A₁B₁，A₂B₂，A₃B₃，A₄B₄及 动平台（A₁A₂A₃A₄A₅）组成，由交流伺服电机(图中未示出)驱动，并采用滚珠丝杠传动； 第三个支链3由转动副、伸缩连杆包括构件51（B₅C）和构件52（CA₅）及动平台A₁A₂A₃A₄A₅ 组成，由交流伺服电机驱动伸缩连杆，其刀架固联于动平台。虽然重型并联机床弹性动力 学更精确，但是弹性动力学太复杂，逆动力学计算的时间长，不能用于控制系统。控制系 统在一个插补周期内需要大量的复杂计算，计算方法不同，插补的时间也就不同，这种重 型机床的插补周期一般不超过2ms。而且在设计计算机控制系统时，通常要求运动控制计 算时间小于控制周期的10%，从而需要重型并联机床的数控系统在0.2ms的时间内完成插补 计算、运动学逆解计算、动力学逆解计算和各控制器计算。在上述各项计算任务中，动力 学逆解计算最为耗时，应该在0.1ms以内完成动力学逆解的计算。

重型并联机床的运动控制效果主要依赖于模型精度，因此建立重型并联机床高精度动 力学模型，并且能够应用于实时控制，对提高重型并联机床精度，提高其实用化水平具有 重要意义。

发明内容

本发明的目的是公开一种面向实时控制应用的过约束重型并联机床动力学建模方法， 克服目前刚体动力学模型精度低、弹性动力学模型不能应用于实时控制系统的缺陷。

本发明的动力学建模方法，包括以下步骤：

1）根据过约束重型并联机床的结构，在运动学分析的基础上，求出机床各个主动关 节的位置，速度和加速度；以及各构件质心在机床不同坐标系中的位置，速度和加速度， 同时还根据雅可比矩阵建立动平台和各个主动关节的速度关系；

2）根据对过约束重型并联机床各个构件进行受力分析的情况,采用牛顿‐欧拉方法建立 过约束重型并联机床中各个构件的力平衡方程和力矩方程；

3）根据过约束重型并联机床的结构参数和刚度较差支链的杆件轴向变形情况，建立 过约束重型并联机床支链中的轴向变形和动平台的输出误差的协调方程；

4）最后，将力平衡方程、力矩方程和误差协调方程联立建立动力学模型，并求出模 型中的参数。

本发明提供一种在重型并联机床的初始设计阶段，面向实时控制应用的重型并联机床 高精度动力学建模法方法，其特点和有益效果是：

1.与现有的并联机床动力学建模方法相比，本方法所建立的动力学模型考虑了刚度较 差杆件的变形，与刚体动力学模型相比，提高了模型精度。此外，这种方法适合驱动冗余 和过约束并联机床的动力学建模。

2.本专利所建立的动力学模型计算效率高，满足控制系统实时性要求，可以应用于并 联机床基于模型的控制方法中，从而提高机床精度。

附图说明

以下结合附图和具体实施方式对本发明作进一步详细说明。

图1是本发明对驱动冗余过约束并联机床结构的分析图；

图2是本发明的动平台误差示意图。

图3是本发明的动平台的位置误差和角度误差示意图，其中(a)为x和y方向的位置误 差，(b)为转动误差。

具体实施方式

本发明的方法实施方式包括以下步骤如下：

1)根据过约束重型并联机床的结构，通过机床进行的运动学分析，求出机床各个关节 的位置，速度和加速度；以及各构件质心在机床不同坐标系中的位置，速度和加速度，同 时还根据雅可比矩阵建立该机床的动平台和各个主动关节的速度关系；具体包括：

11）建立固定在机床机架上的固定坐标系O-XY，如图1所示，O点为机床两个立柱 与地面的接触点C₁、C₃的中点；X轴水平向右，Y轴垂直向上；建立在动平台上的固定坐 标系O′-xy，动平台的中心位置记为O′，x轴水平向右，y轴垂直向上。

在关节点B_i建立以B_iA_i方向为x′轴的动坐标系B_i-x′y′以及与固定坐标系O-XY平 行的动坐标系B_i-xy。根据图中的几何位置关系得到：

r+a_i＝b_i+q_ie₂+ln_i，i＝1,2,3,4    (1)

其中，r＝[x y]^T为O′点在O-XY坐标系中的位置矢量，a_i是A_i点在O′-xy坐标系 中的位置矢量，a_i＝[a_ix a_iy]^T，b_i表示从点O到C_i的矢量，q_i表示关节点B_i在坐标系 O-XY中的Y坐标，e₂＝[0 1]^T，l和n_i分别表示连杆A_iB_i的长度和单位矢量。

连杆A₅B₅的限制方程可以表示为

r+a₅＝b₅+(l₅₁+l₅₂)n₅    (2)

其中，a₅为A₅在O′-xy坐标系的位置矢量，b₅为从点O到点B₅的矢量，l₅₁和l₅₂分别为 构件51和构件52的长度，n₅为连杆A₅B₅的单位矢量。

12）表达式（1）的两边同时对时间求一阶导数，得出各滑块的速度为：

${\stackrel{·}{q}}_i=\frac{{n_i}^T\stackrel{·}{r}}{n_{\mathrm{iy}}},(i=\mathrm{1,2,3,4}),{\stackrel{·}{q}}_5=0---\left(3\right)$

其中，为O′点在O-XY坐标系中的速度矢量，连杆A_iB_i和A₅B₅的角速度分别为： ${\omega }_i=\frac{(\stackrel{·}{r}-{\stackrel{·}{q}}_i{e}_2)}{l}·\left(\begin{array}{c}{-n}_{\mathrm{iy}}\\ n_{\mathrm{ix}}\end{array}\right),$${\omega }_5=\frac{\stackrel{·}{r}}{l_{51}+l_{52}}·\left(\begin{array}{c}{-n}_{5y}\\ n_{5x}\end{array}\right),$其中，n_ix和n_iy分为n_i在X和Y轴上的分量 (i＝1,2,3,4,5)。

构件51及构件52的质心在O-XY坐标系中的速度分别表示为：

${\stackrel{·}{r}}_{g51}=\frac{l_{51}{\omega }_5}{2}\left(\begin{array}{c}-n_{5y}\\ n_{5x}\end{array}\right)---\left(4\right)$

${\stackrel{·}{r}}_{g52}=\frac{{\stackrel{·}{l}}_{52}}{2}{n}_5+(l_{51}+\frac{l_{52}}{2}){\omega }_5\left(\begin{array}{c}{-n}_{5y}\\ n_{5x}\end{array}\right)---\left(5\right)$

13）表达式（1）的两边同时对时间求二阶导数，得出各滑块的加速度为：

${\stackrel{··}{q}}_i=\frac{n_i\stackrel{··}{r}+{\omega }_i^2{l}_i}{n_{\mathrm{iy}}},(i=\mathrm{1,2,3,4}),{\stackrel{··}{q}}_5=0---\left(6\right)$

其中，为O′点在O-XY坐标系中的加速度矢量，连杆A_iB_i和A₅B₅的角加速度分别 为：${ϵ}_i=\frac{\stackrel{··}{r}-{\stackrel{··}{q}}_i{e}_2}{l}·\left(\begin{array}{c}-n_{\mathrm{iy}}\\ n_{\mathrm{ix}}\end{array}\right),$${ϵ}_5=\frac{\stackrel{··}{r}·{\left(\begin{array}{cc}-n_{5y}& n_{5x}\end{array}\right)}^T-2{\stackrel{·}{l}}_{52}{\omega }_5}{l_{51}+l_{52}}.$构件51和构件52的质心在坐标 系O-XY中的加速度分别为：

${\stackrel{··}{r}}_{g51}=\frac{l_{51}{ϵ}_5}{2}\left(\begin{array}{c}-n_{5y}\\ n_{5x}\end{array}\right)-\frac{l_{51}{\omega }_5^2{n}_5}{2}---\left(7\right)$

${\stackrel{··}{r}}_{g52}=\frac{l_{52}}{2}{n}_5+{\stackrel{·}{l}}_{52}{\omega }_5\left(\begin{array}{c}-n_{5y}\\ n_{5x}\end{array}\right)+(l_{51}+\frac{l_{52}}{2}){ϵ}_5\left(\begin{array}{c}-n_{5y}\\ n_{5x}\end{array}\right)-(l_{51}+\frac{l_{52}}{2}){\omega }_5^2{n}_5---\left(8\right)$

2）根据对过约束重型并联机床各个构件进行受力分析的情况，利用牛顿-欧拉方法分 别写出各个构件的牛顿方程和欧拉方程，具体包括：

21）连杆A_iB_i的力平衡方程如公式（9）所示：

$F_{\mathrm{Bi}}+F_{\mathrm{Bi}}+m_{i}g-m_i{\stackrel{··}{r}}_{\mathrm{gi}},(i=\mathrm{1,2,3,4})---\left(9\right)$

其中，m_i是A_iB_i杆的质量，F_Ai是连杆A_iB_i作用在动平台上的力，且F_Ai＝[F_Aix′ F_Aiy′]^T， 是杆件A_iB_i质心加速度矢量，g＝[0 g]^T，且g为重力加速度,F_B1和F_B2是滑块B₁B₂作 用在连杆A₁B₁和A₂B₂的约束力，F_B3和F_B4是滑块B₃B₄作用在连杆A₃B₃和A₄B₄的约束力， 并且F_Bi＝[F_Bix′ F_Biy′]^T。连杆A_iB_i的力矩方程如公式（10）所示：

$F_{\mathrm{Ai}}.l\left(\begin{array}{c}-n_{\mathrm{iy}}\\ n_{\mathrm{ix}}\end{array}\right)-m_i(g+{\stackrel{··}{q}}_i).r_{\mathrm{cxi}}-m_i{\stackrel{··}{r}}_{\mathrm{ci}}·\left(\begin{array}{c}-r_{\mathrm{cyi}}\\ r_{\mathrm{cxi}}\end{array}\right)-J_i{ϵ}_i=0,(i=\mathrm{1,2,3,4})---\left(10\right)$

其中，J_i是连杆A_iB_i基于B_i-xy坐标系的转动惯量，r_ci是连杆A_iB_i在B_i-xy坐标系下的位 置矢量，r_cxi和r_cyi分别是x轴和y轴的分量。

连杆A₅B₅的力平衡方程和力矩方程分别如公式（11）和（12）所示：

$F_{A5}+F_{B5}+m_{51}(g-{\stackrel{··}{r}}_{g51})+m_{52}(g-{\stackrel{··}{r}}_{g52})=0---\left(11\right)$

$F_{A5}.(l_{51}+l_{52})\left(\begin{array}{c}-n_{5y}\\ n_{5x}\end{array}\right)-m_{51}g.r_{\mathrm{cx}51}-m_{51}{\stackrel{··}{r}}_{g51}\left(\begin{array}{c}-r_{\mathrm{cy}51}\\ r_{\mathrm{cx}51}\end{array}\right)-m_{52}g{.r}_{\mathrm{cx}52}-m_{52}{\stackrel{··}{r}}_{g52}\left(\begin{array}{c}-r_{\mathrm{cy}52}\\ r_{\mathrm{cx}52}\end{array}\right)-(J_{51}+J_{52}){ϵ}_5=0---\left(12\right)$

其中，J₅₁和J₅₂分别是连杆B₅C和CA₅相对于坐标系B₅-xy的惯性矩，m₅₁和m₅₂分别为构 件51和构件52的质量，F_A5和F_B5表示动平台和横梁作用在连杆A₅B₅上作用力，r_cx51和r_cy51分别表示构件B₅C质心在坐标系B₅-xy中位置矢量沿x和y的分量。

21）动平台的力平衡方程公如公式（13）所示：

$F_e+M(g-\stackrel{··}{r})-\underset{i=1}{\overset{5}{\Sigma }}{F}_{\mathrm{Ai}}=0---\left(13\right)$

其中，M为动平台的质量，F_e为施加在动平台上的外力。

设动平台的质心相对于坐标系O′-xy的位置矢量为R_c＝[R_cx R_cy]^T，则动平台的欧拉方程 如公式（14）所示：

$-M·g·R_{\mathrm{cx}}-M\stackrel{··}{r}\left(\begin{array}{c}-R_{\mathrm{cy}}\\ R_{\mathrm{cx}}\end{array}\right)+M_e-\underset{i=1}{\overset{5}{\Sigma }}(F_{\mathrm{Ai}}·\left(\begin{array}{c}-A_{\mathrm{iy}}\\ A_{\mathrm{ix}}\end{array}\right))=0---\left(14\right)$

其中，M_e是施加在平台上的力矩，A_ix和A_iy是点O′到A_i的矢量沿F_Aix'和F_Aiy′方向的分量。

3)根据过约束重型并联机床的结构参数和支链的杆件轴向变形情况，建立过约束重 型并联机床支链中的轴向变形和动平台的输出误差的协调方程如下：

(在实际实施工作时，由于机床中的各个支链都很长，所以施加在它们上的力产生变形。 杆件的变形是较大的，对动力学模型会产生很大的影响。)

连杆A_iB_i的轴向变形如公式（15）：

${\delta }_i={\int }_0^l{F}_u\mathrm{du}=\frac{-n_i}{E{S}_i}(F_{\mathrm{Bi}}+m_i\mathrm{gl}/2-m_i{\stackrel{··}{r}}_{\mathrm{gi}}l/2),(i=\mathrm{1,2,3,4})---\left(15\right)$

其中，E表示连杆A_iB_i的弹性模量，δ_i表示连杆A_iB_i的轴向变形，S_i是连杆A_iB_i的横 截面积，F_u表示连杆A_iB_i上到B_i距离为u点所受的力，是连杆A_iB_i质心的加速度。

同样可以计算出连杆A₅B₅的轴向变形如公式(16)：

${\delta }_5=\frac{1}{2E·S_{52}}(2F_{{A5x}^{\prime }}+m_{52}{l}_{52}{n}_5(g-{\stackrel{··}{r}}_{g52}))+\frac{1}{2E·S_{51}}{m}_{51}{l}_{51}{n}_5(g-{\stackrel{··}{r}}_{g51})---\left(16\right)$

其中，S₅₁和S₅₂是构件51和构件52的横截面积。

理论上，动平台不能转动，但是考虑到杆件轴向变形，动平台具有微小的转动dγ，如 图2所示。连杆A_iB_i的轴向变形和动平台位置误差之间关系如公式（17）：

δ_i＝n_i(a′_i-a_i)＝n_i((T-I)a_i+[dx dy]^T)(i=1,2,3,4,5)    (17)

式中，dx为X方向误差，dy为Y方向误差，$T=\left(\begin{array}{cc}\cos \left(\mathrm{d\gamma }\right)& \sin \left(\mathrm{d\gamma }\right)\\ -\sin \left(\mathrm{d\gamma }\right)& \cos \left(\mathrm{d\gamma }\right)\end{array}\right),$a′_i表示动平台 发生微小转动之后点A_i在坐标系O′-xy中的位置矢量，I为二维的单位方阵。

方程(17)可以写成

$\left(\begin{array}{c}{\delta }_1\\ {\delta }_2\\ {\delta }_3\\ {\delta }_4\\ {\delta }_5\end{array}\right)=Q\left(\begin{array}{c}\mathrm{dx}\\ \mathrm{dy}\\ \mathrm{d\gamma }\end{array}\right)---\left(18\right)$

式中，$Q=\left(\begin{array}{ccc}{n}_{1x}& n_{1y}& (n_{1x}{a}_{1y}-n_{1y}{a}_{1x})\\ n_{2x}& n_{2y}& (n_{2x}{a}_{2y}-n_{2y}{a}_{2x})\\ n_{3x}& n_{3y}& (n_{3x}{a}_{3y}-n_{3y}{a}_{3x})\\ n_{4x}& n_{4y}& (n_{4x}{a}_{4y}-n_{4y}{a}_{4x})\\ n_{5x}& n_{5y}& (n_{5x}{a}_{5y}-n_{5y}{a}_{5x})\end{array}\right).$

联立方程（15）、（16）和（18）建立方程组如公式（19）所示：

$Q\left(\begin{array}{c}\mathrm{dx}\\ \mathrm{dy}\\ \mathrm{d\gamma }\end{array}\right)=\left(\begin{array}{c}\frac{1}{2E{S}_1}(2F_{A1x^{\prime }}+m_{1}l{n}_1(g-{\stackrel{··}{r}}_{g1}))\\ \frac{1}{2E{S}_2}(2F_{A2x^{\prime }}+m_{2}l{n}_2(g-{\stackrel{··}{r}}_{g2}))\\ \frac{1}{2E{S}_3}(2F_{A3x^{\prime }}+m_{3}l{n}_3(g-{\stackrel{··}{r}}_{g3}))\\ \frac{1}{2E{S}_4}(2F_{A4x^{\prime }}+m_{4}l{n}_4(g-{\stackrel{··}{r}}_{g4}))\\ \frac{1}{2E{S}_{52}}(2F_{A5x^{\prime }}{l}_{52}+m_{52}{l}_{52}{n}_5(g-{\stackrel{··}{r}}_{g52}))+\frac{m_{51}{l}_{51}{n}_5(g-{\stackrel{··}{r}}_{g51})}{2E{S}_{51}}\end{array}\right)---\left(19\right)$

4）最后，联立构件的力平衡方程(9)、(11)、(13)和力矩方程(10)、(12)、(14)以及连杆 轴向变形和动平台输出误差的协调方程(19)得到动力学模型，并求出模型中的参数；所述 参数包括：连杆A_iB_i给动平台作用力沿杆长方向分量F_Aix'、动平台的位置误差dx、dy及转 动误差dγ，滑块作用在连杆上的力F_Bi及驱动力，具体为首先联立公式(13)，(14)和(19)， 求出F_Aix'(i＝1,2,3,4,5)、dx、dy和dγ，然后根据公式(10)和(12)求出F_Aiy'，进一步根据公式 (9)和(11)求出F_Bi，最后对滑块B₁B₂和B₃B₄进行受力分析，根据滑块以及构件B₅C的力平 衡方程可以计算出驱动力。

根据驱动力求解方法进行计算机仿真，与不考虑连杆变形的模型相比，得到的X、Y 轴的位置和转角的误差，见图3(a)，其中，圆形曲线和正常曲线分别是X和Y轴位置误差， 图3(b)中曲线代表角度误差，由此可看出本发明建立的模型较准确，提高了动力学模型的 精度，进而提高了重型机床的精度。

应用本发明执行系统用配置为英特尔Duo2 2.2GHZ CPU和1.96GB内存微机计算过 约束的冗余重型并联机床的逆动力学的时间为0.05ms。计算效率能够满足实时控制系统的 要求。