李徐钰，魏朝阳，徐文东，邵建达

(中国科学院上海光学精密机械研究所，上海 201800)

随动压力分布下的非球面抛光去除函数

李徐钰，魏朝阳*，徐文东，邵建达

(中国科学院上海光学精密机械研究所，上海 201800)

考虑去除函数对数控小工具抛光光学元件精度的影响，提出了如何根据需要加工的非球面参数以及抛光盘参数得到最优去除函数的方法。由于计算非球面上去除函数的核心是准确获得抛光盘与镜面间的动态压力分布，本文提出利用有限元法仿真抛光盘与非球面间的压力分布，并结合经典 Preston 方程与行星运动模型来得到非球面不同位置处的去除函数。 基于随动压力分布模型，分析了沥青盘抛光非球面时在不同抛光位置处去除函数的变化。在曲率半径为1 000 mm的球面上进行了去除函数验证实验。结果表明：基于本文理论得到的去除函线型更接近实际情况，皮尔逊系数达到了 0.977。本文提出的方法可以方便地调整加工位置来得到相应的压力分布，实现去除函数的优化，对提高加工效率与精度有实际指导意义。

非球面抛光；随动压力分布；变曲率去除函数；有限元法

1 引 言

随着现代天文系统、高精度光刻系统等前沿领域的飞速发展，对于大口径、高陡度非球面光学元件的需求也大大增加。由于表面曲率存在变化，非球面元件的加工相比于平面元件以及球面元件的加工要困难的多，如何提高其加工效率与精度是如今超精密光学加工的一个关键课题。

目前国内外非球面光学元件的加工主要依靠计算机控制表面成型(CCOS)，常用的精加工方法有数控小磨头抛光[1]，离子束抛光[2]，磁流变抛光[3]，应力盘抛光[4]，激光抛光等[5]。不管是哪种技术，高精度的面形测量技术、符合实际的去除函数以及精确快速的驻留时间求解算法是提高加工精度的关键。仅从加工的角度来考虑，提高去除函数的准确性是核心，也是限制该方法加工能力的瓶颈。国防科技大学的宋辞[6]针对磁流变加工非球面的去除函数进行了理论建模，并提出了线性与非线性补偿工艺；C.Wang[7]利用有限元分析仿真了气囊抛光在平面上的静态去除函数并进行了实验验证；国防科技大学的聂徐庆[8]利用有限单元法取代桥模型进行分析的光顺工艺的结果与实验结果吻合。对于小磨头抛光工艺，以往分析去除函数时，均按压力均匀分布来处理，北京理工大学的董志超[9]利用有限元分析仿真了平面沥青盘的去除函数，发现整个抛光盘的压力分布呈现中间高，两侧低的形状，依据此理论得到的去除函数经过反复实验修正后与实际去除函数在峰值去除率和体积去除率上的误差均保持在5%以内。Kim[10]提出了在抛光盘夹层中附加非牛顿流体的材料，使抛光盘兼顾了非球面加工对于低频面形的适应性以及高频误差的平滑作用。由于非球面元件相对于平面元件的自由度很高，而且其表面具有一定的曲率，因此在不同的非球面上或同一非球面的不同位置处的压力分布都会发生变化。常用的有限元分析软件虽然可以分析受力分布，但难以总结出非球面上不同位置的分布规律，因此对分析动态的压力分布以及将其应用到实际加工有诸多不变。

本文利用MATLAB软件并结合有限单元法建立数学模型，并进行仿真计算，将得到的压力分布以及去除函数变化数学化，可以便于根据需要加工的非球面参数以及抛光盘参数仿真得到动态的去除函数。加工平面或球面元件时，可以自始至终都用一个去除函数就可完成加工，而加工非球面时第一步需选择合适的加工参数。本文介绍了如何根据具体镜面，调整公转半径参数来得到最优的去除函数。针对沥青盘抛光非球面，仿真了去除函数随径向加工位置的变化趋势。针对在沥青盘上附着聚氨酯抛光垫加工非球面这一工况，在曲率半径为1 000 mm的球面上进行了去除函数验证实验。本文所做工作均是分析去除函数的轮廓特征，去除效率受外界环境影响较大，因此去除效率需要根据加工实验反演校正。

2 理论背景

数控小磨头抛光的去除函数理论研究中，最经典的方程是Preston方程：

(1)

其中p(x,y)为小磨头与元件表面接触的压力分布，v(x,y)为小磨头上一参考点与元件表面的相对速度的分布，K则是取决于材料性质，温度，抛光液浓度等参数的环境系数。

前人在研究去除函数性质时，通常近似认为压力是均匀分布的，但实际上由于非球面加工中，与抛光盘接触的面并不是平面，所以这种近似会带来不小的误差，也会大大降低小磨头抛光这一环节的加工精度，为之后的磁流变，离子束的精修带来更大的工作量。本文提出在进行压力分布仿真计算时，将抛光层近似当作是弹性体，利用有限单元法进行仿真。

3 去除函数分析的有限单元法

有限单元法建模分为单元分割，整体刚度矩阵建立，边界条件添加，最终计算几个步骤。立体结构的分割单元通常有六面体单元、四面体单元以及等参数单元这几类，为方便建模，本文选取正六面体单元进行建模。六面体单元及局部坐标系如图1所示。

图1 正六面体示意图Fig.1 Diagram of cube element

在六面体单元上建立局部坐标系。令(ξ,η,ζ)为归一化局部坐标系，它们满足式(2)关系：

(2)

其中d为所划分立方体的边长，根据图中各个顶点的标号定义形状函数为:

(i=1,2,…,8).

(3)

该单元的位移函数可以表示成形状函数的线性叠加，叠加系数是对应的顶点位移量:

(4)

位移函数的意义是可以将整个单元的位移近似地用顶点位移线性表示，而相应的形状函数值则表示为对应顶点的“权重”，这也是有限单元法可以转化为矩阵运算的基础。单元应变向量可以写成式(5)：

(5)

其中B为几何矩阵，δe为单元顶点位移向量：

B=[B1B2B3B4B5B6]

(6)

而根据著名的胡克定律[11]，单元的应力与应变之间应满足式(7)的关系：

{σ}={σxσyσzτxyτyzτxz}=Dε,

(7)

其中D为单元的弹性矩阵，其值仅取决于材料本身，E为弹性模量，μ为泊松比。根据虚功方程，在一定的虚位移向量Δδ下，有式(8)：

ΔδT·F=∭Δεσdxdydz=ΔδT∭BTDBdxdydz,

(8)

很容易由式(8)得到如下关系:

F=Kδ=∭BTDBdxdydz·δ,

(9)

其中K称为单元刚度矩阵，有限单元法的核心思想就是先构建单元刚度矩阵，再将其整合成整体刚度矩阵，分析整体结构的位移与受力。对于正六面体单元，根据式(9)推导得K：

(10)

其中对于任意的Kij都有如下的形式：

(11)

k33、k22与k11的形式类似，而其余项与k12的形式类似。构建整体刚度矩阵时对所有的节点进行编号，然后写成类似式(9)的形式，将所有单元中对应的分量都加入整体矩阵中，求和就可以得到整体刚度矩阵。本文选择对圆柱形抛光盘进行分析，采用图2所示的划分方式，从抛光盘中心向外根据网格参数进行计数，只有当一个单元的所有顶点都在圆内时，才将此单元列入计算。

图2 抛光盘单元划分Fig.2 Element division of polishing pad

接下来便根据边界条件修改式(9),此时的边界条件可近似表示如下：

(1) 抛光盘最下面一层的节点的z方向位移分布与和其接触的非球面的局部面形相吻合；

(2) 抛光盘最上面一层所有节点的z方向位移相同；

(3) 根据整体的受力平衡有，上表面的节点受到总体向下的抛光力F，下表面节点受到总体向上的抛光力F。

边界条件有多种添加方法，本文采用的是乘大数法以及罚函数法，条件(1)、(2)采用乘大数法，而条件(3)采用罚函数法进行添加。放大系数的选取要经过反复考量，以免使计算精度不够。

4 结果与分析

4.1 受力分布仿真结果

由于抛光非球面时需要考虑表面的适应能力，因此抛光盘不宜选取的太大或硬度太高，此处选取聚氨酯作为抛光材料，抛光盘半径为15 mm，厚度为8 mm。将顶点曲率半径为2 000 mm，偏心率e2为1的抛物非球面作为分析对象，分析位置定在顶点处，由于本文侧重点在于分析去除函数的相对变化，因此Preston方程中的常系数设为1，具体的值需要根据实际实验去除量来调整。得到的抛光盘下表面的压力分布如图 3所示。可以看到由于选取的是凸非球面，因此在中心位置的压力较大，而边缘位置的压力较小。在加工过程中去除函数往往是对一个平面元件经过实验得到，相对于平面元件上抛光盘相对均匀的压力分布，非球面上由非均匀压力所产生的去除函数肯定与去除函数点实验得到的结果有较大差别。

图3 顶点处压力分布Fig.3 Pressure distribution in vertex

4.2 去除函数仿真优化

小磨头抛光通常采用的是行星模式[12]，磨头的运动由公转和自转复合组成。此处根据通常的加工参数，设置公转速度为πr/s,自转速度设置为2πr/s，公转半径为8 mm，得到的去除函数如图4(a)所示，再将自转速度设置为-2πr/s，公转半径保持8 mm不变，得到的去除函数如图4(b)所示。不难发现采用公自转方向相反的加工模式并将公自转转速比控制在-1/2较为合适，此时可以获得更高的去除率，得到的去除函数轮廓也更接近高斯分布。为了分析在小磨头加工非球面时公转半径对去除函数的影响，将公自转转速比设置为-1/2，偏心距从0 mm到14 mm每隔2 mm取值，再将得到的去除函数沿x方向的数据取出，得到图5。

图4 公转半径为8 mm去除函数Fig.4 TIF with revolution radius of 8 mm

图5 去除函数轮廓随公转半径演变Fig.5 Evolution of TIF profile varies with radius of revolution

从图中可以看出当偏心距接近于0 mm时，去除函数呈现近似M型分布，而当偏心距接近8 mm 时，去除函数慢慢向高斯型分布靠近，当偏心距大于8 mm时，去除函数又开始变回M型分布，因此6～10 mm是偏心距最优的选择范围。另外，可以很直观地发现，在此范围内的去除体积率也是最高的，兼具了加工精确性与高效率。利用本文方法可以方便地根据具体要加工的非球面面形调整偏心距与公自转比等参数，从而得到理想的近高斯型的去除函数。

4.3 非球面上去除函数变化

由于大曲率、高陡度的非球面与抛光盘接触部分的曲率随加工位置变化明显，因此分析出去除函数随加工位置的变化趋势是提高修形精度的关键。由于采用沥青抛光的去除函数较采用聚氨酯抛光更加稳定，且容易去除高频误差，得到更好的表面粗糙度结果，因此利用沥青直接抛光非球面也是实际加工中经常用的工序。通常在抛光非球面时会将沥青盘制作成与顶点部分相吻合的类球面曲面形状，这样可以减小由于表面弯曲导致的沥青断裂的风险，但是由于非球面不同位置的曲率不同，因此在顶点处贴合的很好的沥青盘，到了其他位置将会与局部面形有微小出入。由于沥青的弹性模量达到GPa级别，故微米级别的面形差也会对压力分布产生明显的影响。选取顶点曲率半径为1 500 mm的抛物面作为分析对象，抛光层选取沥青，厚度为6 mm，公转半径设置为8 mm，公自转比为-1/2，在径向0～400 mm间每隔50 mm选取一个加工位置。将去除函数的径向轮廓进行比对，就可以得到图6。

图6 去除函数轮廓随径向距离演变曲线Fig.6 Evolution curves of TIF profile varies with radial distances

从图6可以发现，当加工位置从顶点处向400 mm径向位置移动时，去除函数的径向轮廓由高斯型转向接近M型，主要原因是随着加工位置向外移动，局部曲率半径在变大，因此沥青表面与非球面表面的面形差趋于凹形，导致靠近公转中心的压力变小，去除函数形成了中间低，两侧高的形状。

4.4 实验结果比对

为了进一步验证本文理论的适用性，在4 mm厚的沥青层下附加一层2 mm厚的聚氨酯抛光垫，抛光盘半径为15 mm，附加0.15 MPa的气缸压强。由于本文侧重于分析局部面形变化对抛光压力分布的影响，考虑到实验成本与效率，将抛光盘制作成平面形状，实验工件选取曲率半径为1 000 mm的凹球面，面形数据利用轮廓仪测量得到，扫描方向沿球面的径向方向，去除函数实验中选取的公转速率为60 r/min，自转速率为100 r/min，偏心距为8 mm，经过3 min的去除。将打点之后的球面面形数据记录下来，与实验前的面形数据作差，并进行归一化处理，同时根据本文提出的有限单元法得到该抛光盘在球面上的压力分布，进而得到仿真的去除函数，同时仿真得到将压力均匀分布情况下的去除函数，将归一化处理之后的结果与实验数据进行比对，就得到了图7。

图7 去除函数实验结果Fig.7 Experiment results of TIF

可以发现在偏心距为8 mm时，实验得到的去除函数在中心位置有凹陷，呈现M型分布，而根据均匀分布得到的去除函数则呈现高斯型分布，利用本文理论计算得到的去除函数也呈现了M型分布，显示了在线型分布上的一致性。

为了定量描述结果的相似性，引入了皮尔逊相关系数来描述数据的相关性。皮尔逊相关系数

定义如下：

(12)

经过计算得到均匀分布下的皮尔逊系数为0.965，而基于本文理论仿真得到的去除函数的皮尔逊系数为0.977，优于均匀分布的结果。说明了本文理论对由于压力分布改变导致的局部面形变化的去除函数优化方面有一定效果。

5 结 论

基于有限单元法分析了非球面的小磨头抛光中抛光压力的分布情况，并基于双转子行星模型模拟了相应的去除函数，通过更改参数发现公转半径是影响去除函数性能的最主要参数，经过一系列运算得出在抛光盘半径为15 mm时，公转半径最佳为6～10 mm，此时去除函数具有类高斯型式。分析了沥青抛光盘在非球面上不同加工位置处的去除函数。针对常用的在沥青抛光盘上附着聚氨酯抛光垫的抛光方法，在曲率半径为1 000 mm的球面镜片上进行了去除函数实验。结果表明基于本文理论仿真得到的去除函数在线型与相关性上更接近实际结果。此处未考虑抛光材料在抛光过程中的流变性，后续工作以此为切入点展开。本文为小磨头非球面抛光提供了一种方便、快捷的优化去除函数方法，对于提高加工效率与精度有实际的指导意义。

[1] 王毅，倪颖,余景池.小型非球面数控抛光技术的研究[J].光学 精密工程，2007，15(6):1527-1533. WANG Y, NI Y,YU J CH. Computer- controlled polishing technology for small aspheric lens[J].Opt.PrecisionEng., 2007, 15(6):1527-1533.(in Chinese)

[2] 马占龙，谷勇,彭利荣，等.应用五轴离子束修正超高陡度镜面[J].光学学报，2015,35(6):0622001. MA ZH L, GU Y Q,PENG L R,etal.. Ultra-gradient mirror figuring by 5-Axis ion beam[J].ActaOpticaSinica,2015, 35(6):0622001.(in Chinese)

[3] 杨航，何建国,黄文，等.磁流变抛光去除函数形状动态预测方法[J]. 强激光与粒子束，2015,27(9):092011. YANG H，HE J G，HUANG W，etal..Dynamic prediction of the shape of MRF removal function[J].HighPowerLaserandParticleBeams,2015,27(9):092011.(in Chinese)

[4] BURGE J H. Mirror technologies for giant telescope[Z].TheUniversityofArizonaTucson.

[5] 徐刚，戴玉堂，肖翔.超短激光抛光微加工的数学建模与实验研究[J].中国激光，2013，40(6):0603003. XU G, DAI Y T, XIAO X.Study on mathematical modeling and experiment of ultrashort laser polishing micro-processing [J].ChineseJournalOfLasers, 2013,40(6):0603003.(in Chinese)

[6] 宋辞.离轴非球面光学零件磁流变抛光关键技术研究[D].长沙：国防科学技术大学, 2012. SONG C.Studyonthekeytechnologiesofmanetorheologicalfinishingforoff-axisasphereopticalelements[D]. Changsha:National University of Defence Technology,2012.(in Chinese)

[7] WANG C,WANG Z,YANG X,etal..Modeling of the static tool influence function of bonnet polishing based on FEA[J].IntJAdvManufTechnol.,2014, 74(1-4):341-349.

[8] NIE X Q, LI SH Y,SHI F,etal..Generalized numerical pressure distribution model for smoothing polishing of irregular midspatial frequency errors[J].AppliedOptics,2014,53(6):1020-1027.

[9] DONG ZH CH,CHEN H B, TAM H Y. Modified subaperture tool influence functions of a flat-pitch polisher with reverse-calculated material removal rate[J].AppliedOptics,2014,53(11):2455-2464.

[10] Dae Wook Kim, James H.BurgeRigid conformal polishing tool using non-linear viscro-elastic effect[J].OpticsExpress.2010,Vol.18,No.3.

[11] 薛强,马士进,童智能. 弹性力学[M].北京:北京大学出版社,2006,37-38. XUE Q,MA SH J,TONG ZH N.TheoryofElasticity[M].Beijing:Peking University Press,2006,37-38.(in Chinese)

[12] 李圣怡，戴一帆.大中型光学非球面镜制造与测量新技术[M].北京：国防工业出版社，2011，99-100. LI SH Y,DAI Y F.NewTechnologyforManufacturingandMeasurementofLargeandMiddlescaleAsphericSurface[M]. Beijing:National Defence Industy Press,2011,99-100.(in Chinese)

Tool influence function in aspheric polishing under dynamic pressure distribution

LI Xu-yu, WEI Chao-yang*, XU Wen-dong, SHAO Jian-da

(ShanghaiInstituteofOpticsandFineMechanics,ChineseAcademyofSciences,Shanghai201800,China) *Correspondingauthor,E-mail:siomwei@siom.ac.cn

By taking the effect of Tool Influence Function(TIF) used in the NC program on the polishing accuracy for optical elements into account, the method how to obtain the optimized TIF based on the aspheric parameters to be polished and polishing pad parameters is proposed. As the key to calculate the TIF of aspheric elements is to obtain accurately the dynamic pressure distribution between polishing pad and mirror, this paper analyzes the pressure distribution between polishing pad and mirror by finite element method and obtains the TIFs in different positions on the aspheric surface by classic Preston equation and planet motion. The variation tendency of the pitch polishing pad in polishing an aspheric element is analyzed in different polishing positions based on the dynamic pressure distribution model. The TIF spot experiment is carried out on an aspheric element with a curvature of 1 000 mm. The result shows that the shape of TIF simulated by the theory introduced in the paper is more similar with that of the actual situation, and the Pearson correlation coefficient reaches 0.977. It concludes that the method can obtain the pressure distribution conveniently by changing the parameter of polishing position and can optimize the TIFs, which provides a theoretical basis for improving polishing efficiency and polishing accuracy.

aspheric polishing; dynamic pressure distribution; curvature-changing tool influence function; Finite element method

2016-10-21；

2016-11-17.

中国科学院青年创新促进会科研基金资助项目；国家重大科技专项资助项目

1004-924X(2016)12-3061-07

TH703；TN305.2

:Adoi:10.3788/OPE.20162412.3061

李徐钰(1992-)，男，江苏常州人，博士研究生，2014年于东南大学获得学士学位,主要研究先进非球面制造与检测技术。E-mail: jsczlixuyu@sina.cn

魏朝阳(1980-)，男，副研究员，硕士生导师，2002年于郑州大学光电子技术专业获得学士学位,2007年于中国科学院上海光学精密机械研究所获得博士学位，现为中科院上海光学精密机械研究所在读博士研究生，主要研究方向包括先进光学制造、检测技术和激光破坏等。E-mail: siomwei@siom.ac.cn