基于Matlab的加工钻头螺旋槽砂轮截形的干涉误差分析

贠 良，董黎敏，王泽巍，王 瑞，董 良

（1．天津理工大学 机械工程学院 天津市先进机电系统设计与智能控制重点实验室，天津300384；2．天津市宏维刀具有限公司，天津300110）

标准砂轮在工件成型中应用广泛，选定标准砂轮进行螺旋槽成型加工时虽然可以磨削出满足常规参数的螺旋槽，但是有些特定钻头螺旋槽无法用标准砂轮精确磨削，需要根据已知的螺旋槽进行非标砂轮的设计.在砂轮反求设计中，通常需确定端截形的摆放位置，砂轮与螺旋槽的相对位置，而反求中，砂轮位置的不确定性和其他因素，必然导致反求误差[1-2].干涉误差分析就是对反求得到的砂轮回转体与当时螺旋面的位置关系进行定量分析，从而准确的判断砂轮位置参数设置的合理性，间接预测出所求砂轮磨削出的螺旋槽与已知螺旋槽的差异程度[3-4].

本文首先建立了加工螺旋槽的砂轮数学模型，为砂轮截形的反求提供了理论基础，在此模型的基础上对砂轮位置参数进行干涉误差分析判断，提出了计算干涉误差ΔR的数学模型，基于所求砂轮截形上相应点对干涉误差值数组ΔR进行了计算.然后在砂轮的计算实例中，用Matlab编制截形绘制程序和数组ΔR的计算程序，运行程序得到砂轮截形离散图、拟合图以及ΔR数组.最后根据得到的数组ΔR，可以了解砂轮截形上相应点对应截面的干涉值，并求出最大干涉值.改变砂轮位置参数值，观察所求砂轮截形和最大干涉值的变化.

1 加工螺旋槽的砂轮数学模型

1.1 坐标系的建立及坐标变换

按照磨削加工时钻头、砂轮和机床的相对位置关系建立坐标系（如图1所示）.在机床上建立坐标系O3-x3y3z3，在砂轮上建立坐标系O1-x1y1z1，以上两坐标系固定，位置关系由O1点在机床坐标系O3-x3y3z3上的坐标（a，b，c）和砂轮轴z1与钻头轴线的夹角λ确定.在钻头上建立坐标系O2-x2y2z2.砂轮坐标系先绕x1轴旋转λ角，然后平移（a，b，c），得到机床坐标系.由此可反求得到机床坐标系到砂轮坐标系的转换式为

图1三个坐标系的相对位置Fig.1 Relative position of three coordinate systems

1.2 钻头螺旋面数学模型

钻头坐标系O2-x2y2z2起始位置与机床坐标系O3-x3y3z3重合，此时在O2-x2y2平面上建立的端截形，随着钻头坐标系或钻头做螺旋运动，即绕钻头轴线旋转同时沿其轴向移动，便可得到机床坐标系下的螺旋曲面.则在初始O2-x2y2平面上摆放端截形的位置不同决定了螺旋槽摆放位置的不同，其摆放位置对应的不同方程可以由两个坐标的旋转变换实现，如图2所示.

端截形的方程式为

图2旋转坐标变换Fig.2 Rotation coordinate transformation

式中u2为参变数.

由式（2）可以得到螺旋曲面的方程[5]

式中v2为角度参数，表示端截形从开始位置绕z3轴转过的角度.顺着z3轴看，以顺时针方向转动为正；ρ为螺旋参数，ρ=导程T/2π，表示端截形绕着z3轴转过单位角度时，沿z3轴线方向移动的距离.

1.3 砂轮回转表面数学模型

在砂轮坐标系O1-x1y1z1上，砂轮回转表面可以由平面O1-x1z1上的母线绕z1轴旋转得到.砂轮轴向截形的方程式表示为

式中，u1为参变数.由式（4）可以得到砂轮坐标系下的砂轮回转表面方程式[6]

式中，v1为角度参数，表示砂轮轴向截形所在的平面绕砂轮轴线旋转的角度.

由砂轮回转曲面方程式（5）得到砂轮坐标系下垂直于砂轮轴线的截面方程，表示为

把螺旋面方程式（3）代入机床坐标系到砂轮坐标系的转换式（1）中得到砂轮坐标系下的螺旋面方程式.

2 干涉误差分析

2.1 干涉分析理论

由上述砂轮截形的反求理论可知，在端截形摆放位置和砂轮与螺旋槽相对位置确定的前提下，便可反求出砂轮截形，干涉误差分析就是对反求得到的砂轮回转体与当时螺旋面的位置关系进行判断[7].

为说明干涉分析，假想将求出的砂轮实体模型和已知的螺旋槽，按照反求砂轮时端截形摆放位置和砂轮与螺旋槽的相对位置放置，用垂直于砂轮轴线的平面所截砂轮和螺旋槽，判断平面所截砂轮形成的圆端面与平面所截螺旋槽形成的端截面是否有交集.若在砂轮厚度范围内所有垂直于砂轮轴线所截的两个端截面存在交集即为干涉，若所有平面所截的两个端截面不存在交集即为不干涉.

由剖平面所截砂轮得到截圆，所截螺旋槽得到一截面，如图3所示.通过剖截面所截螺旋槽形成的截面外表曲线上的点到截圆中心的距离，与在此剖面下截圆的半径求差来定量分析干涉误差.如图4所示，A点为外表曲线上任意一点，L为A到圆心O的距离，圆的半径为R，ΔR=L-R，若ΔR大于等于零，说明A点在圆的外面及边缘，不干涉；反之，ΔR小于零，说明A点在圆的内部，说明所求的砂轮会干涉螺旋槽，此时的ΔR值表明了干涉程度.

图3螺旋槽和砂轮剖面图Fig.3 Sectional view of spiral groove and grinding wheel

本文是在垂直于砂轮轴线的剖面中研究的，所以需要求出螺旋槽在此剖面的方程式.根据在同一剖面内砂轮回转曲面和螺旋面的z1相等，由（6）和（7）得到等式

图4圆与外表曲线的位置关系Fig.4 Positional relationship between circle and outer surface curve

根据式（8），给定一个g2（u1）=Z1，即选定一截面，若给定一个u2，便会得到对应的一个v2.其含义是在确定一垂直于砂轮轴线的剖面上，选取之前代表端截形上的点，通过式（8）求解出端截形上的点，成为与此点共螺旋线的所截螺旋面剖面上外表曲线上的点时旋转的角度.将代表端截形点的u2与求出的v2代入到式（7）中，得到砂轮坐标系下与端截形上各点共螺旋线的剖面上外表曲线上点的坐标.在定剖面内，由外表曲线上点的坐标得到此点距所截砂轮形成圆中心的距离，然后与圆的半径相减，得到ΔR的表达式如下.

2.2 基于所求砂轮截形上相应点的干涉分析判断

基于砂轮截形上的相应点来判断干涉情况，每取一个所求砂轮截形上的相应点，即得到一个点坐标（x1，z1），其中z1对应于一个z1剖平面，x1表示z1平面所截砂轮形成截圆圆半径.将z1=g2（u1）代入（8）中，得到一个关于u2与v2的等式，根据等式由端截形上的m个u2分别得到相应的m个v2.将这m组u2，v2和x1=g1（u1）代入式（9）中，得到z1剖面下m个ΔR.由m个砂轮截形上的相应点，得到m×m个ΔR.判断这m×m个ΔR值，若存在负数，就说明存在干涉.负数的大小确定了干涉程度.同样干涉分析的精确度也和反求时端截形上取点量有关.

3 基于Matlab的砂轮截形干涉误差分析求解实例

以半径R=5.347 mm，半芯厚r=1.266 mm，α=30°（钻头螺旋角）的硬质合金钻头为例，编制Matlab计算程序[8].首先定义用到的常量和变量，然后通过数组之间的运算来实现接触点的运算和坐标变换；最后通过Plot命令绘制砂轮轴向截形的散点图及拟合图，再通过所求砂轮截形上的相应点得到关于ΔR的数组.通过对ΔR数组中的每行元素寻找最小值，若最小值大于等于零，说明这行对应的截面不发生干涉，记为零；若最小值小于零，即发生干涉，记为最小值，也为这个截面的最大干涉；得到一个一维数组，再寻找最小值，便会得到整个ΔR数组中的最小值，即最大干涉误差.

3.1 基于Matlab的砂轮轴向截形的绘制

图5砂轮轴向截形Fig.5 Axial section of grinding wheel

由钻头端截形反求砂轮截形时，为简化计算，将端截形型线划分为四段圆弧：圆弧AB为倒角1，大圆弧段BC，小圆弧段CD，圆弧DE为倒角2，如图5所示.每段圆弧对应的方程式和圆心角都不同，求得的砂轮轴向截形也不同，故而在Matlab中使用linspace命令，采用ployfit函数进行曲线拟合，分别在每段圆弧的圆心角范围内均分得到10个u2，很显然端截形由40个点代替进行计算.如图6所示.

图6砂轮轴向截形上的40个相应点的散点图及拟合图Fig.6 Scatter diagram and fitting diagram of 40corresponding points on axial section of grinding wheel

本文以a=80 mm，λ=60°的砂轮截形绘制程序为例，求得砂轮轴向截形上相应点的坐标数组W3，也就是得到了砂轮截形上的40个相应点坐标，见表1.

3.2 所求砂轮截形上相应点的ΔR数组计算程序

由数组W3可知其中有3个坐标值相同，去除掉相同点得到数组W4，即求出轮截形上的37个相应点，同时对应于倒角1圆弧上对应10个离散点，大圆弧上对应9个离散点，小圆弧上对应9个离散点，倒角2圆弧上对应9个离散点.由于端截形由四段圆弧组成，所以求解ΔR数组时，仍需要每段圆弧分别求解最后整合.针对倒角1对应圆弧的ΔR数组DDD求解，由所求砂轮截形上37个相应点集的z1坐标，得到砂轮厚度范围内的垂直于砂轮轴线的37个截面，每取一个截面，就需要求出与此段圆弧10个离散点共螺旋线的此截面所截螺旋面外表线上的10个相应点对应的ΔR，从而求出此段圆弧37个截面对应的ΔR数组DDD，其计算程序如图7所示.

表1 a=80 mm，λ=60°时反求得到砂轮截形上40个相应点的坐标Tab.1 a=80 mm andλ=60，the coordinates of 40 corresponding points on the grinding wheel section are obtained by reverse calculation

图7倒角1圆弧对应的ΔR数组计算程序Fig.7ΔR array calculation program corresponding to chamfer 1 arc

由以上程序同理可得到大圆弧段、小圆弧段和倒角2所对圆弧对应的ΔR数组DDD1，DDD2和DDD3，最后合并四个数组得到总的ΔR数组P（37，37）.数组P中37行代表37个截面，每行的37列元素代表此截面所截螺旋槽外表线上的37个点的ΔR值.

求出的ΔR数组P（37，37）有1 369个元素，通过判断每个元素的大小寻找到每行数据对应的最小值，判断最小值的正负，若ΔRmin≥0，说明此行元素都大于等于0，肯定不会有干涉，用0来表示此行不干涉；若ΔRmin＜0，说明此截面会发生干涉同时最小值也是最大干涉程度，用最小值来表示此行干涉和干涉程度.最终得到数组H（1，37），经过在数组H中寻找最小值得到最大干涉值MAX.数组P的元素搜寻流程图如图8所示.

图8数组P的元素搜寻流程图Fig.8 Element searching flow chart of array P

同样以a=80 mm，λ=60°的砂轮截形反求为例，运行以上程序，可以得到数组H.

H=[-2.06e-15，-1.2e-14，0，0，-1.06e-14，-4.8e-15，0，-3.18e-15，-5.09e-15，0，-1.05e-15，0，0，-4.34e-15，-3.92e-15，-1.45e-14，-3.68e-15，0，0，0，0，0，-5.33e-16，0，-7.41e-16，0，-5.40e-15，0，0，0，0，0，0，0，0，0，0]，MAX=-1.45e-14

由得到的数组H，可以发现反求砂轮截形上37个相应点对应的37个截面上有14个截面会产生干涉误差，其误差值大约都在-e-14或-e-15左右，误差值极小，间接地说明所求砂轮磨削出的螺旋槽与已知螺旋槽差异极小.同时由MAX值，可以得到最大干涉值，可以说明：干涉误差小于-1.45e-14 mm（高精度数控刀具干涉误差应小于-1.5e-10 mm），砂轮位置参数的选取是合理的，由此反求的砂轮是可以加工出已知的螺旋槽的[9-10].

3.3 改变参数的砂轮截形和MAX值

1λ=60°时，a改变，砂轮截形和MAX值的变化

编制完以上程序后，根据一般砂轮直径大小，分别取a=40、50、60、70、80、90（单位为mm），通过更改a的大小得到不同的M文件，运行程序可以得到相应的砂轮截形和MAX值，并将其图形放置于同一坐标系进行比较，可以发现砂轮直径的变化对截形的变化影响较小，最大干涉值变化也不太明显.a变化对截形的影响如图9、10所示.

图9λ=60°a变化时，砂轮截形和最大干涉值Fig.9 Whenλ=60 a changes，the grinding wheel profile and maximum interference value

图10λ=60°，a变化时，所求砂轮截形和最大干涉值及比较Fig.10λ=60，when a changes，the profile of grinding wheel and the maximum interference value are obtained and compared

2a=80 mm，λ改变，砂轮截形和MAX值的变化

a取80 mm时，λ在60°附近取值，观察砂轮截形和MAX值的变化.λ分别取57°、58°、59°、60°、61°、62°，λ变化的砂轮截形和最大干涉值变化如图11所示.

由图11可以发现，在λ=57°、58°、59°、60°、61°时最大干涉值在-e-14或-e-15范围内，取λ=62°时截形直接出现畸形.

图11 a=80mm，λ变化时，砂轮截形和最大干涉值Fig.11 a=80 mm，whenλchanges，grinding wheel profile and maximum interference value

4 结论

由非标整体硬质合金钻头，通过分析砂轮和螺旋槽的基本啮合条件，建立了砂轮坐标系下砂轮截圆和钻头端截形的数学模型.通过剖截面所截螺旋槽形成的截面外表曲线上的点到截圆中心的距离，与在此剖面下截圆的半径求差来定量分析干涉误差，提出了干涉误差数组ΔR的计算方法.在实例计算中，使用Matlab软件编制了绘制砂轮截形和计算数组ΔR的程序.运行程序得到砂轮的轴向截形以及ΔR数组.由ΔR数组得到了最大干涉误差值为-1.45e-14 mm.根据高精度数控刀具干涉误差应小于-1.5e-10 mm，故砂轮位置的选取是合理的.进行干涉误差数组ΔR的计算，能定量地研究砂轮位置参数选取的合理性，同时也为砂轮反求时位置的确定提供了一定的支持.