基于靈敏度分析的主軸熱特性有限元模型修正

葉緒康，李郝林

(上海理工大學(xué) 機(jī)械工程學(xué)院，上海 200093)

?

基于靈敏度分析的主軸熱特性有限元模型修正

葉緒康，李郝林

(上海理工大學(xué) 機(jī)械工程學(xué)院，上海 200093)

提出了一種基于靈敏度分析的模型修正方法。以某機(jī)床廠高精度磨床主軸系統(tǒng)為研究對(duì)象，通過(guò)有限元法仿真主軸系統(tǒng)溫度場(chǎng)，結(jié)合試驗(yàn)測(cè)得幾個(gè)關(guān)鍵測(cè)量點(diǎn)的實(shí)測(cè)溫度，推導(dǎo)出修正邊界條件的靈敏度矩陣，并反求出待修正參數(shù)的優(yōu)化值，修正該主軸系統(tǒng)的熱特性有限元模型。試驗(yàn)結(jié)果表明，基于靈敏度分析的模型修正方法經(jīng)過(guò)4次迭代就能收斂到最優(yōu)值，且經(jīng)修正后的模型參數(shù)能使仿真的溫度場(chǎng)最大誤差從11.99%減小到2.88%。

靈敏度分析；有限元模型修正；熱特性

機(jī)床主軸系統(tǒng)熱變形誤差是影響精密機(jī)床加工精度的主要原因之一[1-3]。近年來(lái)，國(guó)內(nèi)外諸多學(xué)者提出了機(jī)床主軸系統(tǒng)熱特性研究的有限元分析法，其關(guān)鍵所在是建立精確的有限元熱分析模型[4-5]。

由于初始建模的多種簡(jiǎn)化假設(shè)以及邊界條件模擬的不精確等因素的影響，導(dǎo)致了有限元模型計(jì)算得到的結(jié)果與實(shí)驗(yàn)測(cè)試結(jié)果之間存在較大誤差[6]。需要評(píng)估證明其可信度，由此有限元模型修正技術(shù)便應(yīng)運(yùn)而生[7]。近年來(lái)，將有限元模型修正問(wèn)題轉(zhuǎn)化為參數(shù)優(yōu)化問(wèn)題進(jìn)行求解是模型修正領(lǐng)域的發(fā)展趨勢(shì)。其中，基于靈敏度分析的模型修正方法具有物理意義明確、計(jì)算過(guò)程直觀等特點(diǎn)，是模型修正方法中最具應(yīng)用前景的方法之一[8-10]。但由于靈敏度矩陣需要求解目標(biāo)函數(shù)的偏導(dǎo)數(shù)，通常的數(shù)值解析法計(jì)算量大，且當(dāng)特征值與反求參數(shù)間是隱函數(shù)關(guān)系時(shí)，參數(shù)的偏導(dǎo)數(shù)可能無(wú)法求出。故研究靈敏度的數(shù)值計(jì)算方法具有重大意義。

本文以某機(jī)床廠高精度磨床主軸系統(tǒng)為研究對(duì)象，引入Lyness和Moler[11]提出的復(fù)變量求導(dǎo)法，計(jì)算靈敏度矩陣?；谠摲椒▽?duì)主軸有限元模型進(jìn)行修正，以提高模型參數(shù)修正的精度與效率。

1　靈敏度分析

1.1靈敏度分析

熱誤差有限元模型修正過(guò)程實(shí)際上是參數(shù)反求的過(guò)程，將該反問(wèn)題描述為對(duì)應(yīng)的目標(biāo)方程即

(1)

(2)

{Δt}=S{Δp}

(3)

其中，{Δt}={t(p+Δp)}-{t(p)}為溫度殘差矩陣；{Δp}=[Δp1,Δp2,…,Δpn]T為邊界條件參數(shù)的變化量；S為t對(duì)P的一階靈敏度矩陣。即

由于環(huán)境和試驗(yàn)條件限制，儀器本身精度、讀數(shù)技巧、安裝位置和方法等因素的影響，特征測(cè)量值難免存在誤差。為避免模型修正過(guò)程中解的收斂困難，有必要根據(jù)試驗(yàn)中各測(cè)點(diǎn)的測(cè)量精度引入加權(quán)矩陣[10]。由式(3)可看出，第i個(gè)溫度測(cè)點(diǎn)的誤差僅對(duì){Δt}的第i行有影響；而第j個(gè)參數(shù)加載誤差僅對(duì){Δt}的第j列有影響。因此，通過(guò)對(duì)殘差矩陣{Δt}左乘測(cè)點(diǎn)對(duì)角加權(quán)矩陣[Wt]和右乘參數(shù)對(duì)角矩陣[Wp]來(lái)處理。即

(4)

對(duì)n個(gè)設(shè)計(jì)參數(shù)，m個(gè)特征量測(cè)試點(diǎn)，進(jìn)行N次線性無(wú)關(guān)加載，則[Wt]為一對(duì)角陣

(5)

同理，[Wp]也為對(duì)角矩陣

(6)

1.2靈敏度矩陣求解

由式(3)可看出，靈敏度分析法的主要步驟在于靈敏度矩陣的計(jì)算。而求解靈敏度矩陣的關(guān)鍵在于求解特征量對(duì)設(shè)計(jì)參數(shù)的偏導(dǎo)數(shù)，對(duì)于溫度場(chǎng)與邊界參數(shù)間是復(fù)雜的隱式關(guān)系，一般的解析法無(wú)法實(shí)現(xiàn)。本文應(yīng)用復(fù)變量求導(dǎo)法來(lái)確定靈敏度矩陣中的各系數(shù)。復(fù)變量求導(dǎo)法由Lyness和Moler提出，用于計(jì)算實(shí)函數(shù)的偏導(dǎo)數(shù)，被應(yīng)用于輻射反問(wèn)題[11]，取得較滿(mǎn)意的計(jì)算結(jié)果。復(fù)變量求導(dǎo)法的數(shù)學(xué)描述如下：將t(p)轉(zhuǎn)化為以p+ih為自變量的復(fù)變函數(shù)t(p+ih)。當(dāng)h較小時(shí)將其泰勒展開(kāi)為

(7)

對(duì)比式(7)等式兩邊的實(shí)部和虛部可得

(8)

(9)

由式(8)可看出，求解函數(shù)的一階偏導(dǎo)數(shù)只要將原函數(shù)轉(zhuǎn)化為復(fù)變函數(shù)，取其虛部除以一個(gè)小量h(h取10-20)即可。則靈敏度矩陣S轉(zhuǎn)化成

(10)

由式(10)可看出，參數(shù)的偏導(dǎo)數(shù)求解轉(zhuǎn)化為復(fù)域函數(shù)值的計(jì)算，求解n個(gè)帶修正參數(shù)的靈敏度矩陣只需進(jìn)行n次正計(jì)算，且無(wú)需計(jì)算函數(shù)差值，避免了迭代過(guò)程中的急劇振蕩問(wèn)題。

1.3修正參數(shù)反求

工程中常用Matlab編寫(xiě)迭代程序來(lái)反求待修正參數(shù)。取原始模型計(jì)算得來(lái)的設(shè)計(jì)參數(shù)值P0作為初值。由于式(3)是一個(gè)變態(tài)矩陣方程，因此應(yīng)用最小二乘法得到{Δp}穩(wěn)定的解，于是有

(11)

利用式(10)和式(11)實(shí)現(xiàn)參數(shù)迭代反求。對(duì)于第k+1次迭代，Gauss法有下述迭代關(guān)系式

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

當(dāng)?shù)玫綒埐罹仃嚨亩稊?shù)<ε時(shí)迭代終止，即max{Δpk-1}<ε1;max{Δpk-1/pk-1}<ε2時(shí)。其中，ε1和ε2是較小的值。本文計(jì)算時(shí)取ε1=10-6,ε2=10-4。

2　有限元模型修正試驗(yàn)研究

2.1主軸系統(tǒng)溫度場(chǎng)試驗(yàn)概況

本文以某高精度磨床的主軸系統(tǒng)為例，證明所提出模型修正方法的有效性。該主軸系統(tǒng)由主軸、主軸箱、主軸軸承以及連接件等組成，主要熱源是軸承的摩擦發(fā)熱[12-13]。計(jì)算中使用有限元分析軟件建立了經(jīng)過(guò)簡(jiǎn)化的有限元模型，如圖1所示。設(shè)定環(huán)境溫度為23 ℃，當(dāng)主軸轉(zhuǎn)速為2 700r/min時(shí)，計(jì)算得到主軸軸承的發(fā)熱量分別為125W、110W和86W，將軸承發(fā)熱量作為熱源邊界條件以生熱率的形式施加到主軸有限元模型上。主軸與空氣間的對(duì)流換熱系數(shù)可按努謝爾特準(zhǔn)則方程計(jì)算得到對(duì)流傳熱系數(shù)為395W/(m·K)，箱壁自由換熱表面的對(duì)流傳熱系數(shù)取10W/(m·K)。

圖1　機(jī)床主軸系統(tǒng)熱分析有限元模型

根據(jù)有關(guān)機(jī)床溫度傳感器布置的優(yōu)化實(shí)驗(yàn)和理論分析結(jié)果[12]，結(jié)合本文所用磨床主軸的結(jié)構(gòu)特點(diǎn)，確定關(guān)鍵溫度傳感器分布在前后軸承對(duì)應(yīng)的箱體表面，如圖2所示。

圖2　溫度傳感器的布置

本文將利用溫度巡檢儀提取圖2所示的9個(gè)關(guān)鍵測(cè)點(diǎn)的實(shí)驗(yàn)數(shù)據(jù)，并導(dǎo)出25組測(cè)點(diǎn)對(duì)于坐標(biāo)的有限元仿真數(shù)據(jù)，結(jié)合上文所提的靈敏度法對(duì)機(jī)床主軸系統(tǒng)熱分析有限元模型進(jìn)行修正。

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

由于基于靈敏度分析的模型修正的效果以及效率受到待修正參數(shù)的影響較大。因此，正確的選定待修正參數(shù)成為模型修正一項(xiàng)重要任務(wù)[14]。通過(guò)對(duì)有限元模型的邊界條件進(jìn)行相關(guān)性分析，可知熱源對(duì)溫度場(chǎng)的影響因子最大。因此，本文以生熱率pj為待修正參數(shù)，以溫度場(chǎng)ti(p1,p2,…,pn)為特征量對(duì)所選用的高精度磨床主軸系統(tǒng)有限元模型進(jìn)行修正。

由經(jīng)驗(yàn)公式計(jì)算所得生熱率p1=0.50,p2=0.40,p3=0.20(單位:106W/m3)為初始邊界條件進(jìn)行加載，并對(duì)初始模型數(shù)據(jù)進(jìn)行分析，可獲得表1所示的對(duì)比結(jié)果。通過(guò)對(duì)比初始生熱率加載情況下，溫度測(cè)試點(diǎn)的計(jì)算值與實(shí)測(cè)值可發(fā)現(xiàn)，關(guān)鍵測(cè)點(diǎn)處有限元模型的分析偏差均在5%以上，測(cè)點(diǎn)4的偏差更是達(dá)到11.99%。故需要對(duì)初始有限元模型進(jìn)行修正。

表1　修正前模型計(jì)算結(jié)果與測(cè)試結(jié)果比較

基于靈敏度分析的模型修正法關(guān)鍵在于計(jì)算參數(shù)的靈敏度矩陣。本文根據(jù)靈敏度分析方法建立該磨床主軸系統(tǒng)的熱誤差有限元模型修正問(wèn)題的目標(biāo)方程(1)，通過(guò)Matlab對(duì)關(guān)鍵測(cè)點(diǎn)的25組分析值及實(shí)驗(yàn)值進(jìn)行擬合分析。并結(jié)合上文所提的靈敏度矩陣計(jì)算方法，由式(10)計(jì)算靈敏度矩陣。通過(guò)Matlab軟件編寫(xiě)Gauss法參數(shù)反求的迭代程序，可得到如圖3所示的收斂過(guò)程。

結(jié)果表明，經(jīng)過(guò)4次迭代各個(gè)待修正參數(shù)的估計(jì)值均趨于穩(wěn)定，3個(gè)待修正生熱率p1、p2、p3均能收斂到最優(yōu)值。經(jīng)過(guò)靈敏度迭代反求出修正后的生熱率分別為：p1=0.41,p2=0.32,p3=0.18(單位:106W/m3)。

圖3　參數(shù)修正迭代曲線

將反求得到的修正參數(shù)重新加載到有限元模型中，提取測(cè)點(diǎn)的溫度數(shù)據(jù)。表2給出了修正后模型計(jì)算值與測(cè)試值比較?？梢钥闯?，經(jīng)過(guò)靈敏度分析法修正后的模型與實(shí)際主軸的偏差大幅減小，測(cè)點(diǎn)的最大偏差為2.88%。因此，說(shuō)明基于靈敏度分析的模型修正方法可有效提高有限元模型修正精度，經(jīng)修正后的模型能預(yù)測(cè)主軸的熱特性。

表2　修正后模型計(jì)算結(jié)果與測(cè)試結(jié)果比較

(12)

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

(13)

3　結(jié)束語(yǔ)

(1)本文通過(guò)對(duì)某機(jī)床廠的高精度磨床主軸系統(tǒng)進(jìn)行有限元熱分析，結(jié)合幾個(gè)關(guān)鍵測(cè)點(diǎn)的實(shí)際溫度測(cè)量值，運(yùn)用基于靈敏度分析的模型修正法對(duì)主軸有限元模型進(jìn)行修正，經(jīng)修正后的模型溫度場(chǎng)與測(cè)試結(jié)果之間的最大偏差從11.99%減小到2.88%，證明與實(shí)驗(yàn)相結(jié)合的機(jī)床有限元仿真能使所建立的模型更加精確；

(2)靈敏度分析模型修正方法的修正對(duì)象是生熱率、材料、幾何參數(shù)等邊界條件，修正后的模型物理意義明確，修正的結(jié)果可方便地應(yīng)用到有限元分析軟件，便于工程應(yīng)用；

(3)靈敏度分析以泰勒展開(kāi)為基礎(chǔ)，受到計(jì)算效率及計(jì)算難度的限制，由于存在著泰勒展開(kāi)的截?cái)嗾`差。因此，靈敏度分析方法是一種近似方法。盡管靈敏度分析技術(shù)較為成熟，但仍存在不少問(wèn)題：如重根問(wèn)題、迭代收斂性問(wèn)題、病態(tài)矩陣問(wèn)題等，且目標(biāo)函數(shù)選取問(wèn)題；當(dāng)待修正的結(jié)構(gòu)參數(shù)較多時(shí)，靈敏度矩陣計(jì)算量過(guò)大，阻礙了物理參數(shù)型修正方法的實(shí)際應(yīng)用。

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

[2]楊建國(guó),任永強(qiáng),朱衛(wèi)斌,等.數(shù)控機(jī)床熱誤差補(bǔ)償模型在線修正方法研究[J].機(jī)械工程學(xué)報(bào),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]蔣興奇.主軸軸承熱特性及對(duì)速度和動(dòng)力學(xué)性能影響的研究[D].杭州:浙江大學(xué),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]秦仙蓉.基于靈敏度分析的結(jié)構(gòu)計(jì)算模型修正技術(shù)及相關(guān)問(wèn)題研究[D].南京:南京航天航空大學(xué),2001.

[10]袁愛(ài)民.基于靈敏度分析的有限元模型修正技術(shù)若干關(guān)鍵問(wèn)題研究[D].南京:東南大學(xué),2006.

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

[12]李郝林,應(yīng)杏娟.數(shù)控機(jī)床主軸系統(tǒng)熱誤差溫度測(cè)量點(diǎn)的最優(yōu)化設(shè)計(jì)方法[J].中國(guó)機(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-)，男，碩士研究生。研究方向：數(shù)控技術(shù)。李郝林(1961-)，男，博士，教授，博士生導(dǎo)師。研究方向：精密測(cè)量與智能控制等。

10.16180/j.cnki.issn1007-7820.2016.09.022

TG502

A

1007-7820(2016)09-079-04