基于灵敏度分析的主轴热特性有限元模型修正

叶绪康，李郝林

(上海理工大学 机械工程学院，上海 200093)



基于灵敏度分析的主轴热特性有限元模型修正

叶绪康，李郝林

(上海理工大学 机械工程学院，上海 200093)

提出了一种基于灵敏度分析的模型修正方法。以某机床厂高精度磨床主轴系统为研究对象，通过有限元法仿真主轴系统温度场，结合试验测得几个关键测量点的实测温度，推导出修正边界条件的灵敏度矩阵，并反求出待修正参数的优化值，修正该主轴系统的热特性有限元模型。试验结果表明，基于灵敏度分析的模型修正方法经过4次迭代就能收敛到最优值，且经修正后的模型参数能使仿真的温度场最大误差从11.99%减小到2.88%。

灵敏度分析；有限元模型修正；热特性

机床主轴系统热变形误差是影响精密机床加工精度的主要原因之一[1-3]。近年来，国内外诸多学者提出了机床主轴系统热特性研究的有限元分析法，其关键所在是建立精确的有限元热分析模型[4-5]。

由于初始建模的多种简化假设以及边界条件模拟的不精确等因素的影响，导致了有限元模型计算得到的结果与实验测试结果之间存在较大误差[6]。需要评估证明其可信度，由此有限元模型修正技术便应运而生[7]。近年来，将有限元模型修正问题转化为参数优化问题进行求解是模型修正领域的发展趋势。其中，基于灵敏度分析的模型修正方法具有物理意义明确、计算过程直观等特点，是模型修正方法中最具应用前景的方法之一[8-10]。但由于灵敏度矩阵需要求解目标函数的偏导数，通常的数值解析法计算量大，且当特征值与反求参数间是隐函数关系时，参数的偏导数可能无法求出。故研究灵敏度的数值计算方法具有重大意义。

本文以某机床厂高精度磨床主轴系统为研究对象，引入Lyness和Moler[11]提出的复变量求导法，计算灵敏度矩阵。基于该方法对主轴有限元模型进行修正，以提高模型参数修正的精度与效率。

1　灵敏度分析

1.1灵敏度分析

热误差有限元模型修正过程实际上是参数反求的过程，将该反问题描述为对应的目标方程即

(1)

(2)

{Δt}=S{Δp}

(3)

其中，{Δt}={t(p+Δp)}-{t(p)}为温度残差矩阵；{Δp}=[Δp1,Δp2,…,Δpn]T为边界条件参数的变化量；S为t对P的一阶灵敏度矩阵。即

由于环境和试验条件限制，仪器本身精度、读数技巧、安装位置和方法等因素的影响，特征测量值难免存在误差。为避免模型修正过程中解的收敛困难，有必要根据试验中各测点的测量精度引入加权矩阵[10]。由式(3)可看出，第i个温度测点的误差仅对{Δt}的第i行有影响；而第j个参数加载误差仅对{Δt}的第j列有影响。因此，通过对残差矩阵{Δt}左乘测点对角加权矩阵[Wt]和右乘参数对角矩阵[Wp]来处理。即

(4)

对n个设计参数，m个特征量测试点，进行N次线性无关加载，则[Wt]为一对角阵

(5)

同理，[Wp]也为对角矩阵

(6)

1.2灵敏度矩阵求解

由式(3)可看出，灵敏度分析法的主要步骤在于灵敏度矩阵的计算。而求解灵敏度矩阵的关键在于求解特征量对设计参数的偏导数，对于温度场与边界参数间是复杂的隐式关系，一般的解析法无法实现。本文应用复变量求导法来确定灵敏度矩阵中的各系数。复变量求导法由Lyness和Moler提出，用于计算实函数的偏导数，被应用于辐射反问题[11]，取得较满意的计算结果。复变量求导法的数学描述如下：将t(p)转化为以p+ih为自变量的复变函数t(p+ih)。当h较小时将其泰勒展开为

(7)

对比式(7)等式两边的实部和虚部可得

(8)

(9)

由式(8)可看出，求解函数的一阶偏导数只要将原函数转化为复变函数，取其虚部除以一个小量h(h取10-20)即可。则灵敏度矩阵S转化成

(10)

由式(10)可看出，参数的偏导数求解转化为复域函数值的计算，求解n个带修正参数的灵敏度矩阵只需进行n次正计算，且无需计算函数差值，避免了迭代过程中的急剧振荡问题。

1.3修正参数反求

工程中常用Matlab编写迭代程序来反求待修正参数。取原始模型计算得来的设计参数值P0作为初值。由于式(3)是一个变态矩阵方程，因此应用最小二乘法得到{Δp}稳定的解，于是有

(11)

利用式(10)和式(11)实现参数迭代反求。对于第k+1次迭代，Gauss法有下述迭代关系式

{p}k+1={p}k+{Δpk}

当迭代得到残差矩阵的二范数<ε时迭代终止，即max{Δpk-1}<ε1;max{Δpk-1/pk-1}<ε2时。其中，ε1和ε2是较小的值。本文计算时取ε1=10-6,ε2=10-4。

2　有限元模型修正试验研究

2.1主轴系统温度场试验概况

本文以某高精度磨床的主轴系统为例，证明所提出模型修正方法的有效性。该主轴系统由主轴、主轴箱、主轴轴承以及连接件等组成，主要热源是轴承的摩擦发热[12-13]。计算中使用有限元分析软件建立了经过简化的有限元模型，如图1所示。设定环境温度为23 ℃，当主轴转速为2 700r/min时，计算得到主轴轴承的发热量分别为125W、110W和86W，将轴承发热量作为热源边界条件以生热率的形式施加到主轴有限元模型上。主轴与空气间的对流换热系数可按努谢尔特准则方程计算得到对流传热系数为395W/(m·K)，箱壁自由换热表面的对流传热系数取10W/(m·K)。

图1　机床主轴系统热分析有限元模型

根据有关机床温度传感器布置的优化实验和理论分析结果[12]，结合本文所用磨床主轴的结构特点，确定关键温度传感器分布在前后轴承对应的箱体表面，如图2所示。

图2　温度传感器的布置

本文将利用温度巡检仪提取图2所示的9个关键测点的实验数据，并导出25组测点对于坐标的有限元仿真数据，结合上文所提的灵敏度法对机床主轴系统热分析有限元模型进行修正。

2.2基于灵敏度分析的有限元模型修正

由于基于灵敏度分析的模型修正的效果以及效率受到待修正参数的影响较大。因此，正确的选定待修正参数成为模型修正一项重要任务[14]。通过对有限元模型的边界条件进行相关性分析，可知热源对温度场的影响因子最大。因此，本文以生热率pj为待修正参数，以温度场ti(p1,p2,…,pn)为特征量对所选用的高精度磨床主轴系统有限元模型进行修正。

由经验公式计算所得生热率p1=0.50,p2=0.40,p3=0.20(单位:106W/m3)为初始边界条件进行加载，并对初始模型数据进行分析，可获得表1所示的对比结果。通过对比初始生热率加载情况下，温度测试点的计算值与实测值可发现，关键测点处有限元模型的分析偏差均在5%以上，测点4的偏差更是达到11.99%。故需要对初始有限元模型进行修正。

表1　修正前模型计算结果与测试结果比较

基于灵敏度分析的模型修正法关键在于计算参数的灵敏度矩阵。本文根据灵敏度分析方法建立该磨床主轴系统的热误差有限元模型修正问题的目标方程(1)，通过Matlab对关键测点的25组分析值及实验值进行拟合分析。并结合上文所提的灵敏度矩阵计算方法，由式(10)计算灵敏度矩阵。通过Matlab软件编写Gauss法参数反求的迭代程序，可得到如图3所示的收敛过程。

结果表明，经过4次迭代各个待修正参数的估计值均趋于稳定，3个待修正生热率p1、p2、p3均能收敛到最优值。经过灵敏度迭代反求出修正后的生热率分别为：p1=0.41,p2=0.32,p3=0.18(单位:106W/m3)。

图3　参数修正迭代曲线

将反求得到的修正参数重新加载到有限元模型中，提取测点的温度数据。表2给出了修正后模型计算值与测试值比较。可以看出，经过灵敏度分析法修正后的模型与实际主轴的偏差大幅减小，测点的最大偏差为2.88%。因此，说明基于灵敏度分析的模型修正方法可有效提高有限元模型修正精度，经修正后的模型能预测主轴的热特性。

表2　修正后模型计算结果与测试结果比较

(12)

Δfmax=max(|Δfi|),i=1,2,…,n

(13)

3　结束语

(1)本文通过对某机床厂的高精度磨床主轴系统进行有限元热分析，结合几个关键测点的实际温度测量值，运用基于灵敏度分析的模型修正法对主轴有限元模型进行修正，经修正后的模型温度场与测试结果之间的最大偏差从11.99%减小到2.88%，证明与实验相结合的机床有限元仿真能使所建立的模型更加精确；

(2)灵敏度分析模型修正方法的修正对象是生热率、材料、几何参数等边界条件，修正后的模型物理意义明确，修正的结果可方便地应用到有限元分析软件，便于工程应用；

(3)灵敏度分析以泰勒展开为基础，受到计算效率及计算难度的限制，由于存在着泰勒展开的截断误差。因此，灵敏度分析方法是一种近似方法。尽管灵敏度分析技术较为成熟，但仍存在不少问题：如重根问题、迭代收敛性问题、病态矩阵问题等，且目标函数选取问题；当待修正的结构参数较多时，灵敏度矩阵计算量过大，阻碍了物理参数型修正方法的实际应用。

[1]BryanJ.Internationalstatusoftherma1errorresearch[J].KeynotePaperAnnalsoftheCIRP,1990,39(2):645-656.

[2]杨建国,任永强,朱卫斌,等.数控机床热误差补偿模型在线修正方法研究[J].机械工程学报,2003, 39(3):81-84.

[3]QiuJian,LiuChunshi.ThermalerrorsofplanertypeNCmachinetoolsanditsimprovementmeasures[J].JournalofMechanicalEngineering,2010,48(21):149-157.

[4]MianNS,FletcherS,LongstaffAP,etal.EfficientestimationbyFEAofmachinetooldistortionduetoenvironmentaltemperatureperturbations[J].PrecisionEngineering, 2013,37(2):372-379.

[5]蒋兴奇.主轴轴承热特性及对速度和动力学性能影响的研究[D].杭州:浙江大学,2001.

[6]BermanA,NagyEJ.Improvementofalargeanalyticalmodelusingtestdata[J].AIAAJournal,1983,21(8):1168-1173.

[7]NaeemSMian.EfficientestimationbyFEAofmachinetooldistortionduetoenvironmentaltemperatureperturbations[J].PrecisionEngineering,2013(37):372-379.

[8]FoxRL,KapoorMP.Ratesofchangeofeigenvaluesandeigenvectors[J].AIAAJournal,1968,6(2):2426-2429.

[9]秦仙蓉.基于灵敏度分析的结构计算模型修正技术及相关问题研究[D].南京:南京航天航空大学,2001.

[10]袁爱民.基于灵敏度分析的有限元模型修正技术若干关键问题研究[D].南京:东南大学,2006.

[11]LynessJN,MolerCB.Numericaldifferentiationofanalyticfunctions[J].SocietyforIndustrial&AppliedMathematics, 1981,7(4):512-526.

[12]李郝林,应杏娟.数控机床主轴系统热误差温度测量点的最优化设计方法[J].中国机械工程,2010,21(7):804-808.

[13]JinKyungChoi,DaiGilLee.Thermalcharacteristicsofthespindlebearingsystemwithagearlocatedonthebearingspan[J].InternationalJournalofMachineTools&Manufacture, 1998，38(9):1017- 1030.

[14]FriswellM,MottersheadJE,FriswellMichael.Finiteelementmodelupdatinginstructuraldynamics[M].NewYork:WillyPress,2007.

Updating Finite Element Model of Spindle Thermal Properties Based on the Sensitivity Analysis

YE Xukang, LI Haolin

(School of Mechanical Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China)

A theory of finite element (FE) model updating based on the sensitivity analysis is proposed. A finite element model of spindle systems for a high precision grinder is established to analyze the thermal properties. However, the outcomes of the theoretical and experimental modal analyses do not match. Therefore, the analytical models of the structures need to be updated according to the experimental test results. This study addresses an updating algorithm to modify the numerical models by deducing the sensitivity analysis for unknown properties. The result shows that the FE model updating based on the sensitivity analysis can reduce the error between the experimental and analytical results from11.99% to 2.88% by only four iterations.

sensitivity analysis; finite element model updating; thermal properties

2015- 12- 24

叶绪康(1991-)，男，硕士研究生。研究方向：数控技术。李郝林(1961-)，男，博士，教授，博士生导师。研究方向：精密测量与智能控制等。

10.16180/j.cnki.issn1007-7820.2016.09.022

TG502

A

1007-7820(2016)09-079-04