CN110287553B - 一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法 - Google Patents

一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法 Download PDF

Info

Publication number
CN110287553B
CN110287553B CN201910498260.9A CN201910498260A CN110287553B CN 110287553 B CN110287553 B CN 110287553B CN 201910498260 A CN201910498260 A CN 201910498260A CN 110287553 B CN110287553 B CN 110287553B
Authority
CN
China
Prior art keywords
machine tool
error
numerical control
matrix
tij
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.)
Active
Application number
CN201910498260.9A
Other languages
English (en)
Other versions
CN110287553A (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.)
Beijing University of Technology
Original Assignee
Beijing 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 Beijing University of Technology filed Critical Beijing University of Technology
Priority to CN201910498260.9A priority Critical patent/CN110287553B/zh
Publication of CN110287553A publication Critical patent/CN110287553A/zh
Application granted granted Critical
Publication of CN110287553B publication Critical patent/CN110287553B/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
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Geometry (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Evolutionary Computation (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Numerical Control (AREA)

Abstract

本发明公开了一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析新方法，属于机床精度设计领域，具体涉及到多轴数控机床的空间误差建模方法以及基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法。本发明运用多体***理论建立数控机床空间误差模型，根据蒙特卡洛模拟采样机理，对机床加工误差模型进行全局灵敏度分析，获得影响机床加工误差的关键几何误差参数，在机床设计的初期阶段，提出新的机床设计理念，为提升数控机床的加工精度以及关键几何误差参数补偿奠定了理论基础。

Description

一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析 方法
技术领域
本发明涉及一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，属于机床精度设计技术领域。
背景技术
近年来，数控机床作为制造业的工业母机，其功能在不断的完善，为了满足日益复杂形状工件精密加工的要求，五轴机床得到越来越广泛地应用。五轴机床结构复杂，涉及误差因素较多，加工成型机理也比三轴机床复杂。相比于三轴机床，五轴机床具有更好的通用性和灵活性、更高的效率和精度等优点。由于五轴数控机床各项技术的不断成熟，大大提高了复杂曲面零件的可加工性和加工精度。
数控机床的精度指标主要有加工精度、定位精度和重复定位精度，其中加工精度是数控机床追求的最终精度，反映出机械制造业的制造能力和发展水平。影响机床加工精度的因素主要包括机床各个零部件的几何误差、热误差、切削力引入误差和刀具磨损等，其中几何误差对加工精度的影响最大。机床的几何误差包括直线度误差，滚摆误差，偏摆误差，定位误差，颠摆误差和垂直度误差等。数控机床各项几何误差参数的作用相互耦合，对加工精度的影响权重各不相同。如何有效的识别出对数控机床加工精度影响较大的关键几何误差项，并且在机床设计的初期阶段有效的控制它们是有效提高机床加工精度的关键问题。
这一关键问题的解决方法分为两个步骤：
第一、基于多体***运动学理论，建立机床的空间误差模型；
目前国内外学者已经开展了许多关于机床精度建模方法的研究，先后出现了二次关系模型法、几何建模法、误差矩阵法、刚体运动学法和多体***理论法。基于多体***运动学理论，将五轴机床抽象为多体***，用拓扑结构图以及低序体阵列表来描述机床的结构和各个体之间的关联关系，分析数控机床的几何误差，建立广义坐标系，用相邻体间的特征矩阵表示位置关系，用齐次变换矩阵表示多体***间的相互关系，最终建立机床的空间误差模型；
第二、结合空间误差模型，提出基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法；
由于拟蒙特卡洛方法的计算准确度高，收敛速度快，所以选用拟蒙特卡洛方法进行计算。将机床误差模型作为黑盒处理，所有几何误差参数在定义域内的随机采样值作为输入变量，机床加工误差为输出值，根据各输入变量的方差对输出值方差的影响，评估各项几何误差参数对机床加工误差模型的影响程度。灵敏度系数越大，说明该属性对模型输出的影响越大，根据灵敏度系数的大小对各项几何误差参数进行排序，最后筛选出灵敏度系数较大的几何误差项。本发明采用基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法来分析五轴数控机床关键几何误差参数。
本发明基于多体***运动学理论，建立了机床的空间误差模型，然后结合空间误差模型提出了一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，最后有效的辨识出了五轴数控机床的关键几何误差参数。
发明内容
本发明的目的是提供一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法。通过建立机床的空间误差模型和误差敏感度分析模型，辨识出关键几何误差参数，为提高加工精度奠定基础。
为了实现上述目的，本发明采用的技术方案为一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，本发明基于多体***运动学理论，建立了机床的空间误差模型，然后结合空间误差模型提出了一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，最后有效的辨识出了数控机床的关键几何误差。
本方法具体包括如下步骤：
步骤一：建立数控机床的空间误差模型；
基于多体***运动学理论，用多体***示意图以及低序体阵列表对机床的结构进行简化，如图2和表1所示。分析数控机床的几何误差参数，建立广义坐标系，用相邻体间的特征矩阵表示各零部件之间的位置关系，用齐次变换矩阵表示多体***间的相互关系；
步骤1.1建立数控机床的拓扑结构；
数控机床是一个多分支的复杂***，从B1处分为两个分支，除了B1体外每个物体都有一个相邻的较低序体，用Ln(j)表示，称为低序体阵列表，如表1所示，j表示物体的序号，j＝1,2,3…n，n表示机床所包含典型体的个数；
表1：数控机床低序体阵列
L<sup>0</sup>(j) 1 2 3 4 5 6
L<sup>1</sup>(j) 0 1 1 3 4 5
L<sup>2</sup>(j) 0 0 0 1 3 4
L<sup>3</sup>(j) 0 0 0 0 1 3
L<sup>4</sup>(j) 0 0 0 0 0 1
L<sup>5</sup>(j) 0 0 0 0 0 0
典型体的编号规则如下：
首先任选一典型体为B1，然后沿远离B1体的方向，依自然增长的数列依次标定每个物体的序号；
步骤1.2数控机床的几何误差分析
在空间坐标系中任意物体均有6个自由度，在运动过程中必然产出6项与位置有关的误差，包括3项线位移误差和3项角位移误差，X、Y、Z三条导轨间存在3项不垂直度误差，C轴与X、Y轴，A轴与Y、Z轴之间共存在4项垂直度误差，因此共37项误差如表2所示；
表2：数控机床几何误差参数
Figure BDA0002089337570000021
Figure BDA0002089337570000031
步骤1.3建立数控机床的特征矩阵；
根据数控机床各部件之间的运动关系，建立各相邻体之间的变换矩阵如表3所示；
表3：相邻体间的变换矩阵
Figure BDA0002089337570000032
Figure BDA0002089337570000041
其中：[Tij]p表示Bj体相对于Bi体的相对位置变换矩阵；
[Tij]pe表示Bj体相对于Bi体的相对位置误差变换矩阵；
[Tij]s表示Bj体相对于Bi体的相对运动变换矩阵；
[Tij]se表示Bj体相对于Bi体的相对运动误差变换矩阵；
x表示X轴平移的距离；
y表示Y轴平移的距离；
z表示Z轴平移的距离；
a表示A轴转动的角度；
c表示C轴转动的角度；
几何误差的敏感度分析方法使用过程中，忽略除几何误差之外的所有误差因素；
步骤1.4建立机床的空间误差模型
理想情况下相邻体运动关系模型的建立；
设P点为Bj体上任意一点，P在Bi体坐标系Oi-XiYiZi中的位置矩阵表达式为；
Pji＝[Tij]p[Tij]srj(1)
式中：Pji为P点在坐标系Oi-XiYiZi中的位置矩阵表达式；
rj为P点在坐标系Oj-XjYjZj中的位置矩阵表达式；
[Tij]p表示Bj体相对于Bi体的相对位置变换矩阵；
[Tij]s表示Bj体相对于Bi体的相对运动变换矩阵；
有误差情况下相邻体运动关系模型的建立；
设P点为Bj体上任意一点，P在Bi体坐标系Oi-XiYiZi中的位置矩阵表达式为；
Pji＝[Tij]p[Tij]pe[Tij]s[Tij]serj(2)
式中：Pji为P点在坐标系Oi-XiYiZi中的位置矩阵表达式；
rj为P点在坐标系Oj-XjYjZj中的位置矩阵表达式；
[Tij]p表示Bj体相对于Bi体的相对位置变换矩阵；
[Tij]pe表示Bj体相对于Bi体的相对位置误差变换矩阵；
[Tij]s表示Bj体相对于Bi体的相对运动变换矩阵；
[Tij]se表示Bj体相对于Bi体的相对运动误差变换矩阵；
刀具中心点在刀具坐标系中的坐标为：
rt＝[0,0,l,1]T (3)
l表示刀具长度；
下标t表示刀具
理想情况下刀具中心点P按“数控机床-工件”分支到惯性坐标系中的位置矩阵表达式：
Figure BDA0002089337570000051
理想情况下刀具中心点P按“数控机床-刀具”分支到惯性坐标系中的位置矩阵表达式：
Figure BDA0002089337570000052
数控指令精密加工方程：
Pw I＝Pt I (6)
理想情况下，数控指令到工件坐标系中的位置矩阵表达式：
Figure BDA0002089337570000053
实际情况下刀具中心点P按“机床-工件”分支到惯性坐标系中的位置矩阵表达式：
Figure BDA0002089337570000054
实际情况下刀具中心点P按“机床-刀具”分支到惯性坐标系中的位置矩阵表达式：
Figure BDA0002089337570000055
实际情况下，数控指令到工件坐标系中的位置矩阵表达式：
Figure BDA0002089337570000061
则数控机床的空间误差模型表示为：
E＝rw-rw I (11)
步骤二：基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析；
全局灵敏度分析是一种针对***模型的研究方法，而不是针对模型的某些特定解的分析。对机床关键几何误差因素进行灵敏度分析的本质，是要求解机床各项几何误差对机床加工误差的固有影响程度。应用基于方差的灵敏度分析方法，可以同时考虑所有几何误差因素对几何误差模型的影响。Sobol提出的基于蒙特卡洛的灵敏度分析法，是目前最为常用的灵敏度分析方法。
令I为单位向量，In为n维单位立方空间，x∈In，以下每一项变量的积分区间均为[0,1]。设***方程为y＝f(x)，其中y为模型输出，x＝(x1,x2,...xn)为模型的n个输入变量。f(x)的高维模型分解表示(analysis of variance，ANOVA)为公式(12)。
Figure BDA0002089337570000062
其中f0＝E(y)，fi＝E(y|xi)-E(y)，fij＝E(y|xi,xj)-fi-fj-E(y)
公式(12)中，f(x)被分解为2n项。当各变量相互独立且正交时，这种分解方式唯一。
对公式(12)两边同时求方差，得：
Figure BDA0002089337570000063
其中Vi＝V(fi(xi))＝V[E(y|xi)]，Vij＝V(fij(xi,xj))＝V(E(y|xi,xj))-Vi-Vj，
Vijk＝V(fijk(xi,xj,xk))＝V(E(y|xi,xj,xk))-Vij-Vik-Vjk-Vi-Vj-Vk
令Si＝Vi/V(y)，Sij＝Vij/V(y)，方程两边同时除以V(y)，得：
Figure BDA0002089337570000064
其中Si为1阶灵敏度指标，表示每一项输入对输出方差的影响程度，为主灵敏度指标；Sij为2阶灵敏度指标，为xi和xj对输出方差的联合影响程度减去各自的主灵敏度指标，表示xi和xj的2阶交叉灵敏度指标；更高阶的灵敏度指标的定义以此类推。Si越大，xi对输出方差的影响程度越大。
根据公式13，计算***的一阶灵敏度需要计算两项参数V(y)和V[E(y|xi)]。
设y为n个输入变量的函数
y＝f(x1,x2,...xn) (15)
假设各自独立变量的联合概率密度函数为
Figure BDA0002089337570000065
由此可得y的期望和方差表达如下
Figure BDA0002089337570000071
Figure BDA0002089337570000072
令xj(j＝1,2,...n)取固定值
Figure BDA0002089337570000073
则
Figure BDA0002089337570000074
其中
Figure BDA0002089337570000075
和
Figure BDA0002089337570000076
分别为输入变量
Figure BDA0002089337570000077
时***输出的方差和期望。
Figure BDA0002089337570000078
通过xj的概率密度函数计算
Figure BDA0002089337570000079
的期望，可以消除其对
Figure BDA00020893375700000710
数值的依赖。
Figure BDA00020893375700000711
Figure BDA00020893375700000712
V(y)＝E[V(y|xj)]+V[E(y|xj)] (23)
由此可以得出以下关系
Figure BDA00020893375700000713
令
Figure BDA00020893375700000714
Figure BDA00020893375700000715
以上Uj的方程可以用下式表达
Figure BDA00020893375700000716
F(X)由2n-1个独立变量决定。对每一项变量进行N次采样后，可以估算f和f*输出值的数学期望。f的输出值由N×n维的输入变量采样矩阵计算。将该矩阵的第j列固定，其他数据进行重采样，可以计算f*的输出值。根据已知的X的分布函数，构造两个N×n的随机矩阵A、B。
Figure BDA0002089337570000081
Figure BDA0002089337570000082
将矩阵B的第j列用矩阵A的第j列替代，得矩阵Cj。
Figure BDA0002089337570000083
将以上样本矩阵A、Cj作为输入，带入***方程，得到输出响应
yiA＝f(xi1,xi2,...xin) (27)
Figure BDA0002089337570000084
对于离散变量x，Uj可由下式估计
Figure BDA0002089337570000085
Figure BDA0002089337570000086
Figure BDA0002089337570000087
由式计算输入变量xi的灵敏度指标为
Figure BDA0002089337570000088
根据灵敏度系数的大小确定几何误差参数对机床空间误差影响程度；灵敏度系数小说明该项几何误差对机床空间误差影响较小，反之，灵敏度系数大说明该项几何误差对机床空间误差影响较大。然后根据灵敏度系数列出误差源参数敏感度队列；根据敏感度分析结果，在机床设计初期，对相应的主要误差进行严格的限制，从而提高机床的加工精度。
与现有技术相比，本发明具有如下有益效果。
1、在机床设计的初期阶段，可为工程师在机床的设计，装配和加工方面提供指导，从根本上提高机床的加工精度。
2、本文证明了全局灵敏度分析方法可以有效地分析出机床误差模型的敏感参数，找出机床误差参数与加工误差模型存在的固有联系，并且该方法可以有效的避免敏感度分析结果失真的情况。针对某一特定型号的机床，该方法可在机床设计阶段分析出机床的敏感几何误差参数，为机床的设计和制造提供指导。
附图说明
图1为本发明方法的实施流程图
图2为多体***示意图
图3为五轴机床的拓扑结构图
图4为五轴机床的结构示意图
图5为位置点有关的误差参数示意图
图6为位置点无关的误差参数示意图
图7为各项几何误差对机床加工精度的敏感度排序图
具体实施方式
本发明以五轴高架横梁移动龙门数控铣床为例，对上述五轴数控铣床关键性几何误差的辨识方法进行验证。
具体包括如下步骤：
步骤一：以五轴数控机床为例，建立机床的空间误差模型；
基于多体***运动学理论，用拓扑结构图以及低序体阵列表来描述机床的结构和各个体之间的关联关系，如图3和表1所示。分析数控机床的几何误差，建立广义坐标系，用相邻体间的特征矩阵表达位置关系，用齐次变换矩阵表示多体***间的相互关系；
步骤1.1建立五轴数控机床的拓扑结构；
该机床的结构如图4所示。包括床身、工作台、刀具、工件、X轴、Y轴、Z轴、B轴、C轴、主轴；
五轴数控机床是一个多分支的复杂***，从B1处分为两个分支，除了B1体外每个物体都有一个相邻的较低序体，用Ln(j)表示，称为低序体阵列表，如表1所示，j表示物体的序号(j＝1,2,3…n)，n表示机床所包含典型体的个数；
表1：数控机床低序体阵列
L<sup>0</sup>(j) 1 2 3 4 5 6
L<sup>1</sup>(j) 0 1 1 3 4 5
L<sup>2</sup>(j) 0 0 0 1 3 4
L<sup>3</sup>(j) 0 0 0 0 1 3
L<sup>4</sup>(j) 0 0 0 0 0 1
L<sup>5</sup>(j) 0 0 0 0 0 0
典型体的编号规则如下：
首先任选一典型体为B1，然后沿远离B1体的方向，依自然增长的数列依次标定每个物体的序号，从***的一个分支到另一个分支，直到全部物体都标定完毕；
步骤1.2分析五轴数控机床的几何误差；
在空间坐标系中任意物体均有6个自由度，在运动过程中必然产出6项误差，3项线位移误差和3项角位移误差，这些都是与位置点有关的误差，如图5所示。X、Y、Z三条导轨间存在3项不垂直度误差，C轴与X、Y轴，A轴与Y、Z轴之间共存在4项垂直度误差，如图6所示。因此共37项误差如表2所示；
表2：五轴数控机床几何误差参数
Figure BDA0002089337570000101
步骤1.3建立五轴数控机床的特征矩阵；
根据数控机床各部件之间的运动关系，可建立各相邻体之间的变换矩阵如表3所示；
表3：相邻体间的变换矩阵
Figure BDA0002089337570000102
Figure BDA0002089337570000111
本发明是一种几何误差的敏感度分析方法，使用过程中忽略除几何误差之外的所有误差因素；
步骤1.4建立机床的空间误差模型；
刀具中心点在刀具坐标系中的坐标为：
rt＝[0,0,l,1]T(1)
l表示刀具长度；
下标t表示刀具
理想情况下刀具中心点P按“机床-工件”分支到惯性坐标系中的位置矩阵表达式：
Pw I＝[T12]p[T12]srw I(2)
理想情况下刀具中心点P按“机床-刀具”分支到惯性坐标系中的位置矩阵表达式：
Pt I＝[T13]p[T13]s[T34]p[T34]s[T45]p[T45]s[T56]p[T56]srt I(3)
数控指令精密加工方程：
Pw I＝Pt I (4)
理想情况下，数控指令到工件坐标系中的位置矩阵表达式：
rw I＝([T12]p[T12]s)-1[T13]p[T13]s[T34]p[T34]s[T45]p[T45]s[T56]p[T56]srt I(5)
实际情况下刀具中心点P按“机床-工件”分支到惯性坐标系中的位置矩阵表达式：
Pw＝[T12]p[T12]pe[T12]s[T12]serw (6)
实际情况下刀具中心点P按“机床-刀具”分支到惯性坐标系中的位置矩阵表达式：
Figure BDA0002089337570000121
实际情况下，数控指令到工件坐标系中的位置矩阵表达式：
Figure BDA0002089337570000122
则机床的空间误差模型表示为：
E＝rw-rw I (9)
步骤二：基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析；
全局灵敏度分析是一种针对***模型的研究方法，而不是针对模型的某些特定解的分析。对机床关键几何误差因素进行灵敏度分析的本质，是要求解机床各项几何误差对机床加工误差的固有影响程度。应用基于方差的灵敏度分析方法，可以同时考虑所有几何误差因素对几何误差模型的影响。Sobol提出的基于蒙特卡洛的灵敏度分析法，是目前最为常用的灵敏度分析方法。
令I为单位向量，In为n维单位立方空间，x∈In，以下每一项变量的积分区间均为[0,1]。设***方程为y＝f(x)，其中y为模型输出，x＝(x1,x2,...xn)为模型的n个输入变量。f(x)的高维模型分解表示(analysis ofvariance，ANOVA)为公式(10)。
Figure BDA0002089337570000123
其中f0＝E(y)，fi＝E(yxi)-E(y)，fij＝E(yxi,xj)-fi-fj-E(y)
公式(10)中，f(x)被分解为2n项。当各变量相互独立且正交时，这种分解方式唯一。
对公式(10)两边同时求方差，得：
Figure BDA0002089337570000124
其中Vi＝V(fi(xi))＝V[E(y|xi)]，Vij＝V(fij(xi,xj))＝V(E(y|xi,xj))-Vi-Vj，
Vijk＝V(fijk(xi,xj,xk))＝V(E(y|xi,xj,xk))-Vij-Vik-Vjk-Vi-Vj-Vk。
令Si＝Vi/V(y)，Sij＝Vij/V(y)...，方程两边同时除以V(y)，得：
Figure BDA0002089337570000125
其中Si为1阶灵敏度指标，表示每一项输入对输出方差的影响程度，为主灵敏度指标；Sij为2阶灵敏度指标，为xi和xj对输出方差的联合影响程度减去各自的主灵敏度指标，表示xi和xj的2阶交叉灵敏度指标；更高阶的灵敏度指标的定义以此类推。Si越大，xi对输出方差的影响程度越大。
根据公式11，计算***的一阶灵敏度需要计算两项参数V(y)和V[E(y|xi)]。
设y为n个输入变量的函数
y＝f(x1,x2,...xn) (13)
假设各自独立变量的联合概率密度函数为
Figure BDA0002089337570000131
由此可得y的期望和方差表达如下
Figure BDA0002089337570000132
Figure BDA0002089337570000133
令xj(j＝1,2,...n)取固定值
Figure BDA0002089337570000134
则
Figure BDA0002089337570000135
其中
Figure BDA0002089337570000136
和
Figure BDA0002089337570000137
分别为输入变量
Figure BDA0002089337570000138
时***输出的方差和期望。
Figure BDA0002089337570000139
通过xj的概率密度函数计算
Figure BDA00020893375700001310
的期望，可以消除其对
Figure BDA00020893375700001311
数值的依赖。
Figure BDA00020893375700001312
Figure BDA00020893375700001313
V(y)＝E[V(y|xj)]+V[E(y|xj)] (21)
由此可以得出以下关系
Figure BDA00020893375700001314
令
Figure BDA00020893375700001315
Figure BDA00020893375700001316
以上Uj的方程可以用下式表达
Figure BDA0002089337570000141
F(X)由2n-1个独立变量决定。对每一项变量进行N次采样后，可以估算f和f*输出值的数学期望。f的输出值由N×n维的输入变量采样矩阵计算。将该矩阵的第j列固定，其他数据进行重采样，可以计算f*的输出值。根据已知的X的分布函数，构造两个N×n的随机矩阵A、B。
Figure BDA0002089337570000142
Figure BDA0002089337570000143
将矩阵B的第j列用矩阵A的第j列替代，得矩阵Cj。
Figure BDA0002089337570000144
将以上样本矩阵A、Cj作为输入，带入***方程，得到输出响应
yiA＝f(xi1,xi2,...xin) (25)
Figure BDA0002089337570000149
对于离散变量x，Uj可由下式估计
Figure BDA0002089337570000145
Figure BDA0002089337570000146
Figure BDA0002089337570000147
由式计算输入变量xi的灵敏度指标为
Figure BDA0002089337570000148
根据灵敏度系数的大小确定几何误差参数对机床空间误差影响程度；灵敏度系数小说明该项几何误差对机床空间误差影响较小，反之，灵敏度系数大说明该项几何误差对机床空间误差影响较大。然后根据灵敏度系数列出误差源参数敏感度队列；为了便于分析，将计算结果用柱状图表示，如图7所示。
敏感度分析结果表明：
εy(x),εy(y),εx(x),εy(z),εy(b),εy(c),εx(z)这7项误差源参数对应的误差敏感度系数之和为0.93，其他误差源参数的敏感度系数之和仅为0.07，因此基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法可以有效的识别出对数控机床加工精度影响较大的几何误差项。
综上所述，为了提高机床的加工精度，在机床设计初期阶段，应该对εy(x),εy(y),εx(x),εy(z),εy(b),εy(c),εx(z)7项几何误差进行严格的控制。

Claims (3)

1.一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，其特征在于：基于多体***运动学理论，建立机床的空间误差模型，然后结合空间误差模型，最后辨识出数控机床的关键几何误差；
本方法具体包括如下步骤：
步骤一：建立数控机床的空间误差模型；
基于多体***运动学理论，用多体***示意图以及低序体阵列表对机床的结构进行简化；分析数控机床的几何误差参数，建立广义坐标系，用相邻体间的特征矩阵表示各零部件之间的位置关系，用齐次变换矩阵表示多体***间的相互关系；
步骤1.1建立数控机床的拓扑结构；
数控机床是一个多分支的复杂***，从B1处分为两个分支，除了B1体外每个物体都有一个相邻的较低序体，用Ln(j)表示，称为低序体阵列表，如表1所示，j表示物体的序号，j＝1,2,3…m，m表示机床所包含典型体的个数；
表1：数控机床低序体阵列
L<sup>0</sup>(j) 1 2 3 4 5 6 L<sup>1</sup>(j) 0 1 1 3 4 5 L<sup>2</sup>(j) 0 0 0 1 3 4 L<sup>3</sup>(j) 0 0 0 0 1 3 L<sup>4</sup>(j) 0 0 0 0 0 1 L<sup>5</sup>(j) 0 0 0 0 0 0
典型体的编号规则如下：
首先任选一典型体为B1，然后沿远离B1体的方向，依自然增长的数列依次标定每个物体的序号；
步骤1.2数控机床的几何误差分析
在空间坐标系中任意物体均有6个自由度，在运动过程中必然产出6项与位置有关的误差，包括3项线位移误差和3项角位移误差，X、Y、Z三条导轨间存在3项不垂直度误差，C轴与X、Y轴，A轴与Y、Z轴之间共存在4项垂直度误差，因此共37项误差如表2所示；
表2：数控机床几何误差参数
Figure FDA0003823098390000021
步骤1.3建立数控机床的特征矩阵；
根据数控机床各部件之间的运动关系，建立各相邻体之间的变换矩阵如表3所示；
表3：相邻体间的变换矩阵
Figure FDA0003823098390000022
Figure FDA0003823098390000031
其中：[Tij]p表示Bj体相对于Bi体的相对位置变换矩阵；
[Tij]pe表示Bj体相对于Bi体的相对位置误差变换矩阵；
[Tij]s表示Bj体相对于Bi体的相对运动变换矩阵；
[Tij]se表示Bj体相对于Bi体的相对运动误差变换矩阵；
x表示X轴平移的距离；
y表示Y轴平移的距离；
z表示Z轴平移的距离；
a表示A轴转动的角度；
c表示C轴转动的角度；
几何误差的敏感度分析方法使用过程中，忽略除几何误差之外的所有误差因素；
步骤1.4建立机床的空间误差模型
理想情况下相邻体运动关系模型的建立；
设P点为Bj体上任意一点，P在Bi体坐标系Oi-XiYiZi中的位置矩阵表达式为；
Pji＝[Tij]p[Tij]srj (1)
式中：Pji为P点在坐标系Oi-XiYiZi中的位置矩阵表达式；
rj为P点在坐标系Oj-XjYjZj中的位置矩阵表达式；
[Tij]p表示Bj体相对于Bi体的相对位置变换矩阵；
[Tij]s表示Bj体相对于Bi体的相对运动变换矩阵；
有误差情况下相邻体运动关系模型的建立；
设P点为Bj体上任意一点，P在Bi体坐标系Oi-XiYiZi中的位置矩阵表达式为；
Pji＝[Tij]p[Tij]pe[Tij]s[Tij]serj (2)
式中：Pji为P点在坐标系Oi-XiYiZi中的位置矩阵表达式；
rj为P点在坐标系Oj-XjYjZj中的位置矩阵表达式；
[Tij]p表示Bj体相对于Bi体的相对位置变换矩阵；
[Tij]pe表示Bj体相对于Bi体的相对位置误差变换矩阵；
[Tij]s表示Bj体相对于Bi体的相对运动变换矩阵；
[Tij]se表示Bj体相对于Bi体的相对运动误差变换矩阵；
刀具中心点在刀具坐标系中的坐标为：
rt＝[0,0,l,1]T (3)
l表示刀具长度；
下标t表示刀具；
理想情况下刀具中心点P按“数控机床-工件”分支到惯性坐标系中的位置矩阵表达式：
Figure FDA0003823098390000051
理想情况下刀具中心点P按“数控机床-刀具”分支到惯性坐标系中的位置矩阵表达式：
Figure FDA0003823098390000052
数控指令精密加工方程：
Figure FDA0003823098390000053
理想情况下，数控指令到工件坐标系中的位置矩阵表达式：
Figure FDA0003823098390000054
实际情况下刀具中心点P按“机床-工件”分支到惯性坐标系中的位置矩阵表达式：
Figure FDA0003823098390000055
实际情况下刀具中心点P按“机床-刀具”分支到惯性坐标系中的位置矩阵表达式：
Figure FDA0003823098390000056
实际情况下，数控指令到工件坐标系中的位置矩阵表达式：
Figure FDA0003823098390000061
则数控机床的空间误差模型E表示为：
E＝rw-rw I (11)
步骤二：基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析；
令I为单位向量，In为n维单位立方空间，x∈In，以下每一项变量的积分区间均为[0,1]；设***方程为y＝f(x)，其中y为模型输出，x＝(x1,x2,...xn)为模型的n个输入变量；f(x)的高维模型分解表示为公式(12)；
Figure FDA0003823098390000062
其中f0＝E(y)，fi＝E(y|xi)-E(y)，fij＝E(y|xi,xj)-fi-fj-E(y)
公式(12)中，f(x)被分解为2n项；当各变量相互独立且正交时，这种分解方式唯一；
对公式(12)两边同时求方差，得：
Figure FDA0003823098390000063
其中Vi＝V(fi(xi))＝V[E(y|xi)]，Vij＝V(fij(xi,xj))＝V(E(y|xi,xj))-Vi-Vj，Vijk＝V(fijk(xi,xj,xk))＝V(E(y|xi,xj,xk))-Vij-Vik-Vjk-Vi-Vj-Vk
令Si＝Vi/V(y)，Sij＝Vij/V(y)，方程两边同时除以V(y)，得：
Figure FDA0003823098390000064
其中Si为1阶灵敏度指标，表示每一项输入对输出方差的影响程度，为主灵敏度指标；Sij为2阶灵敏度指标，为xi和xj对输出方差的联合影响程度减去各自的主灵敏度指标，表示xi和xj的2阶交叉灵敏度指标；更高阶的灵敏度指标的定义以此类推；Si越大，xi对输出方差的影响程度越大；
根据公式13，计算***的一阶灵敏度需要计算两项参数V(y)和V[E(y|xi)]；
设y为n个输入变量的函数
y＝f(x1,x2,...xn) (15)
假设各自独立变量的联合概率密度函数为
Figure FDA0003823098390000071
由此可得y的期望和方差表达如下
Figure FDA0003823098390000072
Figure FDA0003823098390000073
令xj(j＝1,2,...n)取固定值
Figure FDA0003823098390000074
则
Figure FDA0003823098390000075
其中
Figure FDA0003823098390000076
和
Figure FDA0003823098390000077
分别为输入变量
Figure FDA0003823098390000078
时***输出的方差和期望；
Figure FDA0003823098390000079
通过xj的概率密度函数计算
Figure FDA00038230983900000710
的期望，可以消除其对
Figure FDA00038230983900000711
数值的依赖；
Figure FDA00038230983900000712
Figure FDA00038230983900000713
V(y)＝E[V(y|xj)]+V[E(y|xj)] (23)
由此得出以下关系
Figure FDA00038230983900000714
令
Figure FDA00038230983900000715
Figure FDA0003823098390000081
以上Uj的方程可以用下式表达
Figure FDA0003823098390000082
F(X)由2n-1个独立变量决定；对每一项变量进行N次采样后，可以估算f和f*输出值的数学期望；f的输出值由N×n维的输入变量采样矩阵计算；将该矩阵的第j列固定，其他数据进行重采样，可以计算f*的输出值；根据已知的X的分布函数，构造两个N×n的随机矩阵A、B；
Figure FDA0003823098390000083
Figure FDA0003823098390000084
将矩阵B的第j列用矩阵A的第j列替代，得矩阵Cj；
Figure FDA0003823098390000085
将以上样本矩阵A、Cj作为输入，带入***方程，得到输出响应
yiA＝f(xi1,xi2,...xin) (27)
Figure FDA0003823098390000086
对于离散变量x，Uj可由下式估计
Figure FDA0003823098390000087
Figure FDA0003823098390000091
Figure FDA0003823098390000092
由式计算输入变量xi的灵敏度指标为
Figure FDA0003823098390000093
根据灵敏度系数的大小确定几何误差参数对机床空间误差影响程度。
2.根据权利要求1所述的一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，其特征在于：灵敏度系数小说明该项几何误差对机床空间误差影响小，反之，灵敏度系数大说明该项几何误差对机床空间误差影响大；然后根据灵敏度系数列出误差源参数敏感度队列。
3.根据权利要求1所述的一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，其特征在于：根据敏感度分析结果，在机床设计初期，对相应的主要误差进行限制，从而提高机床的加工精度。
CN201910498260.9A 2019-06-10 2019-06-10 一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法 Active CN110287553B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910498260.9A CN110287553B (zh) 2019-06-10 2019-06-10 一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910498260.9A CN110287553B (zh) 2019-06-10 2019-06-10 一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法

Publications (2)

Publication Number Publication Date
CN110287553A CN110287553A (zh) 2019-09-27
CN110287553B true CN110287553B (zh) 2022-11-04

Family

ID=68003734

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910498260.9A Active CN110287553B (zh) 2019-06-10 2019-06-10 一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法

Country Status (1)

Country Link
CN (1) CN110287553B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111079252A (zh) * 2019-11-12 2020-04-28 山东交通职业学院 一种用于提高非线性颤振可靠度振动计算精度与计算效率的数据处理方法
CN110837246A (zh) * 2019-11-25 2020-02-25 天津工业大学 五轴数控机床双旋转轴几何误差敏感度分析方法
CN110955979B (zh) * 2019-12-11 2023-11-24 北京工业大学 一种考虑几何误差偏相关性的机床加工精度可靠性灵敏度分析方法
CN111580469B (zh) * 2020-04-24 2021-07-02 成都飞机工业（集团）有限责任公司 一种基于大数据数控机床精度指标数据特征的挖掘方法
CN112001052B (zh) * 2020-08-31 2023-09-29 郑州轻工业大学 一种高精密多轴数控机床精度设计的定量分析方法
CN112257252B (zh) * 2020-10-19 2024-03-01 北京精雕科技集团有限公司 一种仿真分析机床空间误差对工件加工精度影响的方法
CN112989520B (zh) * 2021-04-12 2022-11-22 重庆大学 一种蜗轮母机关键误差项辨识方法
CN113359609B (zh) * 2021-07-06 2022-05-31 宁波大学 五轴数控机床关键几何误差优化配比补偿方法
CN113656902A (zh) * 2021-07-19 2021-11-16 湖南科技大学 一种多轴运动平台工作末端位姿的误差敏感性分析方法
CN113960959B (zh) * 2021-10-25 2023-05-19 安徽科技学院 一种多轴数控机床关键几何误差溯源方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104050316A (zh) * 2014-03-25 2014-09-17 北京工业大学 一种数控机床空间加工误差分布特征分析方法
CN104375460A (zh) * 2014-11-17 2015-02-25 北京工业大学 一种数控机床加工精度可靠性敏感度分析方法
CN108445839A (zh) * 2018-05-06 2018-08-24 北京工业大学 一种基于误差增量的机床精度敏感度分析方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101126808B1 (ko) * 2010-03-02 2012-03-23 경북대학교 산학협력단 다축 제어 기계의 오차 평가 방법 및 장치

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104050316A (zh) * 2014-03-25 2014-09-17 北京工业大学 一种数控机床空间加工误差分布特征分析方法
CN104375460A (zh) * 2014-11-17 2015-02-25 北京工业大学 一种数控机床加工精度可靠性敏感度分析方法
CN108445839A (zh) * 2018-05-06 2018-08-24 北京工业大学 一种基于误差增量的机床精度敏感度分析方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于蒙特卡洛模拟的机床关键几何误差溯源方法;范晋伟等;《北京工业大学学报》;20171110;全文 *
机床误差敏感度分析方法;范晋伟等;《北京工业大学学报》;20190319;全文 *

Also Published As

Publication number Publication date
CN110287553A (zh) 2019-09-27

Similar Documents

Publication Publication Date Title
CN110287553B (zh) 一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法
CN108445839B (zh) 一种基于误差增量的机床精度敏感度分析方法
Li et al. Dual sliding mode contouring control with high accuracy contour error estimation for five-axis CNC machine tools
CN104375460B (zh) 一种数控机床加工精度可靠性敏感度分析方法
CN104965483B (zh) 一种基于稳健设计的多轴数控机床加工精度保持性优化方法
CN103390082B (zh) 一种多轴机床几何精度稳健优配方法
Fu et al. Product-of-exponential formulas for precision enhancement of five-axis machine tools via geometric error modeling and compensation
CN109765848B (zh) 五轴机床运动轴几何误差贡献值计算及影响评估方法
CN112558547B (zh) 一种五轴数控机床平动轴几何误差补偿数据快速优化方法
CN104007700A (zh) 一种基于全局敏感度分析的三轴数控机床的关键性几何误差辨识方法
CN108052747B (zh) 一种基于价值分析法的机床几何精度优化方法
CN103034166A (zh) 一种机床关键性几何误差源识别方法
Zou et al. Sensitivity analysis using a variance-based method for a three-axis diamond turning machine
CN104156519A (zh) 一种面向加工精度可靠度提升的多轴数控机床几何精度设计方法
CN107368637A (zh) 基于区间理论的精密卧式加工中心几何精度优化分配方法
CN113536488B (zh) 基于配准算法的毛坯质量包容性分析和余量优化方法
Jiang et al. A minimal-error-model based two-step kinematic calibration methodology for redundantly actuated parallel manipulators: An application to a 3-DOF spindle head
CN104200063B (zh) 机床空间加工误差的非确定性描述及预测方法
Zhang et al. Kinematic chain optimization design based on deformation sensitivity analysis of a five-axis machine tool
CN116974241B (zh) 面向绿色低碳制造的数控机床几何优化方法及装置
CN113359609A (zh) 五轴数控机床关键几何误差优化配比补偿方法
CN113656902A (zh) 一种多轴运动平台工作末端位姿的误差敏感性分析方法
Song et al. A task-oriented calibration method for a 5-DOF hybrid machining robot based on the equivalent system
CN112989520B (zh) 一种蜗轮母机关键误差项辨识方法
CN114048557A (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