李郝林 應(yīng)杏娟 遲玉倫
上海理工大學(xué),上海,200093
在精密機(jī)床的切削加工中,熱源對(duì)加工精度的影響極大,提高工件的加工精度必須對(duì)機(jī)床的熱變形作定量研究,并在加工過程中作合理控制與補(bǔ)償。機(jī)床導(dǎo)軌與工作臺(tái)的相對(duì)運(yùn)動(dòng)會(huì)產(chǎn)生大量的摩擦熱,造成導(dǎo)軌與工作臺(tái)的熱變形,影響刀具與工件之間的相對(duì)位置關(guān)系,最終影響工件的加工精度。為提高刀具的加工精度,需要獲得導(dǎo)軌整個(gè)行程范圍內(nèi)的熱變形動(dòng)態(tài)變化量,以便進(jìn)行實(shí)時(shí)誤差補(bǔ)償。但是,對(duì)導(dǎo)軌全行程范圍內(nèi)進(jìn)行熱誤差實(shí)時(shí)測(cè)量目前還沒有有效的方法[1-4]。
隨著有限元理論的不斷完善和數(shù)值模擬技術(shù)的日趨成熟,數(shù)值模擬技術(shù)在熱特性分析領(lǐng)域得到了廣泛應(yīng)用,成為熱變形分析的重要手段。有限元數(shù)值模擬技術(shù)可以定量地計(jì)算出溫度分布狀態(tài)和由溫度產(chǎn)生的熱位移、應(yīng)力和應(yīng)變等數(shù)據(jù),這些模擬結(jié)果對(duì)機(jī)床的熱特性分析有一定的指導(dǎo)意義。但由于有限元分析邊界條件的不確定性,故無法給出精確的計(jì)算結(jié)果,只能給出定性的分析結(jié)果。要得到接近實(shí)際值的最優(yōu)計(jì)算結(jié)果,必須將有限元數(shù)值模擬技術(shù)與優(yōu)化技術(shù)相結(jié)合[5-6]。文獻(xiàn)[7-8]提出了有限元模型的實(shí)驗(yàn)校準(zhǔn)方法,將計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果相比較,從而對(duì)有限元模型進(jìn)行校準(zhǔn)。本文提出了采用近似模型優(yōu)化方法進(jìn)行機(jī)床導(dǎo)軌熱變形有限元模型邊界條件優(yōu)化選擇的方法。采用試驗(yàn)設(shè)計(jì)方法進(jìn)行數(shù)據(jù)采樣,為了使數(shù)值模擬結(jié)果逼近實(shí)驗(yàn)結(jié)果,以導(dǎo)軌熱變形計(jì)算誤差平方和作為優(yōu)化目標(biāo)函數(shù),以有限元分析的邊界條件為設(shè)計(jì)變量,采用響應(yīng)面法建立設(shè)計(jì)變量和目標(biāo)函數(shù)的近似模型,利用導(dǎo)軌測(cè)量點(diǎn)位置處熱變形的實(shí)際測(cè)量數(shù)據(jù),修正有限元計(jì)算邊界條件,從而得到接近實(shí)際測(cè)量值的全行程范圍內(nèi)的導(dǎo)軌熱變形誤差實(shí)時(shí)動(dòng)態(tài)變化值。
在研究導(dǎo)軌熱變形誤差有限元分析過程中,由于有限元模型的邊界條件很難在理論上精確確定,故使得有限元計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果之間往往存在明顯的誤差,為了縮小這種誤差,本文借助響應(yīng)面方法進(jìn)行有限元模型修正。
響應(yīng)面法以試驗(yàn)設(shè)計(jì)為基礎(chǔ),用于的處理多變量問題的建模和分析。通過近似構(gòu)造一個(gè)目標(biāo)函數(shù)響應(yīng)y和一組設(shè)計(jì)變量x相關(guān)的低階多項(xiàng)式,以顯式的響應(yīng)面近似模型逼近目標(biāo)函數(shù)響應(yīng)y與設(shè)計(jì)變量x之間復(fù)雜的隱式關(guān)系[9-11]。本文使用該方法,根據(jù)導(dǎo)軌熱變形的實(shí)際測(cè)量數(shù)據(jù),修正有限元模型的邊界條件(對(duì)流傳熱系數(shù)和熱源)。
工程中最廣泛采用的響應(yīng)面近似函數(shù)為二階模型[12-13]:
式中,n 為設(shè)計(jì)變量總數(shù);xi、xj為 設(shè)計(jì)變量;α0、αi、αii、αij(i<j)為多項(xiàng)式的待定系數(shù)。
將上式轉(zhuǎn)換成多元線性回歸模型,可得統(tǒng)一的簡單形式:
其中,β0、βi為待定系數(shù),k是待定系數(shù)βi的個(gè)數(shù)。為了確定βi,需要做m次的獨(dú)立實(shí)驗(yàn),m≥k,通過求解得到相應(yīng)的系數(shù)βi。對(duì)于4個(gè)設(shè)計(jì)變量的情況,待定系數(shù)βi的個(gè)數(shù)k=15,至少進(jìn)行15次獨(dú)立實(shí)驗(yàn)。
設(shè)ε為服從正態(tài)分布N(0,σ2)的響應(yīng)面回歸值與實(shí)驗(yàn)值之間的誤差。為了估計(jì)二次多項(xiàng)式的系數(shù)βi,可以用最小二乘法,使誤差平方和S最小,即
根據(jù)變分原理,令
解式(1)即可得到待定系數(shù)βi的無偏差估計(jì)值。以上是響應(yīng)面法的基本原理及響應(yīng)面系數(shù)求解的具體方法。
1.2.1 復(fù)相關(guān)系數(shù)R2
復(fù)相關(guān)系數(shù)R2描述響應(yīng)面函數(shù)對(duì)實(shí)驗(yàn)數(shù)據(jù)的擬合程度,定義如下:
R2是一個(gè)在[0,1]之間變化的值,其值越接近1說明誤差的影響越小,即回歸方程越精確[14]。
1.2.2 響應(yīng)面函數(shù)的顯著性檢驗(yàn)
通過方差分析和F檢驗(yàn)確定回歸模型是否可以作為有意義的近似模型。
服從F分布[14]。對(duì)于給定的顯著水平α,由F分布表查出Fα(k,m-k-1),若F >Fα(k,m-k-1)成立,則認(rèn)為在α水平下該響應(yīng)面模型有意義。
熱誤差的測(cè)量實(shí)驗(yàn)在精密導(dǎo)軌試驗(yàn)臺(tái)上進(jìn)行,試驗(yàn)裝置和傳感器安裝方式如圖1所示??紤]工作臺(tái)的運(yùn)動(dòng)問題,傳感器選用非接觸位移傳感器,選用 HEIDENHAIN公司的capaNCDT6100精密位移傳感器。對(duì)于數(shù)控主軸磨床類機(jī)床,軸向(Z方向)和縱向(Y方向)對(duì)工件加工精度影響不大,所以本文考慮分析具有典型意義的橫向(X方向)熱變形誤差沿導(dǎo)軌運(yùn)動(dòng)方向的分布,坐標(biāo)方向見圖1。
圖1 機(jī)床導(dǎo)軌熱變形測(cè)量試驗(yàn)裝置
沿導(dǎo)軌運(yùn)動(dòng)方向(Z方向)安裝r個(gè)位移傳感器(本文以3個(gè)傳感器為例說明所提出的計(jì)算方法)。其中1個(gè)位移傳感器放置在中間位置,其余2個(gè)分布在導(dǎo)軌兩端,圖1中C1、C2、C3測(cè)點(diǎn)位置即傳感器安裝位置。工作臺(tái)運(yùn)行速度為60m/min,連續(xù)運(yùn)行240min(4h),每隔一定時(shí)間讀取位移傳感器的測(cè)量值,并根據(jù)測(cè)量值計(jì)算出水平方向的熱變形量誤差。
設(shè)傳感器 1、2、3 的讀數(shù)分別為 δ1(i)、δ2(i)、δ3(i),i=0,1,…,N(N為測(cè)量序號(hào)),本例設(shè)定每20min測(cè)量一次,共測(cè)量12次,則N=12。測(cè)量開始時(shí)主軸熱變形為0,設(shè)傳感器1、2、3的讀數(shù)分別為 δ1(0)、δ2(0)、δ3(0),則第 i次測(cè)量 X 方向 1、2、3傳感器位置處的熱變形量分別為
本文研究的對(duì)象(自制精密導(dǎo)軌試驗(yàn)臺(tái))主要由導(dǎo)軌、滑塊、工作臺(tái)、床身等組成,主要熱源是滾珠導(dǎo)軌與滑塊之間以及滾珠絲杠副的摩擦發(fā)熱。
建立響應(yīng)面模型之前,先通過對(duì)給定邊界條件(發(fā)熱量、對(duì)流系數(shù)等)導(dǎo)軌熱變形進(jìn)行有限元數(shù)值熱變形分析,獲得導(dǎo)軌全程范圍內(nèi)的熱變形結(jié)果,并與實(shí)際實(shí)驗(yàn)結(jié)果比較。
計(jì)算中使用有限元分析軟件建立并經(jīng)過簡化的導(dǎo)軌試驗(yàn)臺(tái)有限元模型,如圖2所示。簡化原則為:①忽略對(duì)計(jì)算結(jié)果影響不大的細(xì)小結(jié)構(gòu);②試驗(yàn)臺(tái)床身機(jī)架、絲杠等零件不參與建模;③床身底部用螺釘固定,在 X、Y、Z方向沒有熱位移。
溫度場(chǎng)分析采用SOLID87單元,熱變形分析轉(zhuǎn)換為結(jié)構(gòu)單元SOLID92。計(jì)算時(shí)假定環(huán)境溫度為20℃。
圖2中C1、C2、C3表示與實(shí)驗(yàn)中傳感器位置對(duì)應(yīng)的測(cè)點(diǎn)位置。當(dāng)工作臺(tái)進(jìn)給速度為60m/min時(shí),計(jì)算得到前后導(dǎo)軌和滑塊表面的熱流密度分別為1000W/m2和3000W/m2;導(dǎo)軌和床身靜止表面與空氣間的對(duì)流換熱系數(shù)取10W/(m2?K);當(dāng)工作臺(tái)運(yùn)動(dòng)速度為60m/min時(shí),計(jì)算得到對(duì)流換熱系數(shù)為50W/(m2?K)[15]。
圖2 導(dǎo)軌試驗(yàn)臺(tái)有限元模型
在上述邊界條件下,計(jì)算得到導(dǎo)軌運(yùn)動(dòng)4h的熱誤差分布如圖3所示。計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)比較,最大誤差為 6.4μm,誤差平方和為813μm2??梢钥闯?計(jì)算值和實(shí)驗(yàn)數(shù)據(jù)有較大的差別,這是因?yàn)橛邢拊P蜔嵩春瓦吔鐥l件的不確定性影響了計(jì)算精度。為使計(jì)算結(jié)果逼近實(shí)驗(yàn)結(jié)果,需要對(duì)影響計(jì)算結(jié)果的有限元模型的邊界條件進(jìn)行優(yōu)化,以便對(duì)機(jī)床導(dǎo)軌熱誤差進(jìn)行優(yōu)化計(jì)算。
圖3 熱誤差比較曲線(優(yōu)化前)
本文提出應(yīng)用近似模型方法進(jìn)行優(yōu)化計(jì)算,必須首先建立對(duì)應(yīng)的數(shù)學(xué)模型,包括選取設(shè)計(jì)變量,設(shè)計(jì)目標(biāo)函數(shù),給出約束條件等。
為了使數(shù)值模擬結(jié)果逼近實(shí)驗(yàn)結(jié)果,取導(dǎo)軌上3個(gè)傳感器安裝位置的熱變形計(jì)算誤差平方和作為目標(biāo)函數(shù),用F(Ω)表示,作為設(shè)計(jì)變量Ω的優(yōu)化計(jì)算指標(biāo)。
對(duì)C1、C2、C3測(cè)點(diǎn)熱變形誤差模擬計(jì)算值Δ1i、Δ2i、Δ3i與實(shí)際測(cè)量值δ1i、δ2i、δ3i進(jìn)行比較 ,其
誤差分別為 ω1i= Δ1i-δ1i;ω2i= Δ2i-δ2i;ω3i=Δ3i-δ3i,i為測(cè)量序號(hào),i=0,1,…,N 。
我們的優(yōu)化目標(biāo)是上述誤差越小越好,所以取其誤差平方和構(gòu)造設(shè)計(jì)變量 Ω的優(yōu)化計(jì)算指標(biāo)
本例中工作臺(tái)運(yùn)動(dòng)4h,每20min讀取一次熱誤差值,i=12。
在對(duì)機(jī)床導(dǎo)軌熱變形誤差進(jìn)行有限元分析的過程中,上述計(jì)算的有限元模型的邊界條件(對(duì)流傳熱系數(shù)和熱源)與實(shí)際有較大誤差,很難從理論上精確確定,使得有限元計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果之間往往存在明顯的誤差。因此本文選擇導(dǎo)軌和滑塊表面對(duì)流系數(shù)(ξ1、ξ2),導(dǎo)軌和滑塊發(fā)熱表面的熱流密度(Q1、Q2)作為設(shè)計(jì)變量 Ω。
在理論計(jì)算的基礎(chǔ)上結(jié)合經(jīng)驗(yàn)選取設(shè)計(jì)變量的因子水平范圍,如表1所示。
表1 設(shè)計(jì)變量因子水平范圍表
為使實(shí)驗(yàn)點(diǎn)分布得比較均勻,將4個(gè)設(shè)計(jì)變量在各自范圍內(nèi)均勻取5水平。為便于計(jì)算研究,將設(shè)計(jì)變量進(jìn)行中心編碼[-1,1]轉(zhuǎn)換,形成新的設(shè)計(jì)變量
4因子的二階響應(yīng)面采用完整二次型,其實(shí)驗(yàn)次數(shù)至少需進(jìn)行(4+1)(4+2)/2=15組[14]。均勻設(shè)計(jì)方法是利用數(shù)論方法在設(shè)計(jì)空間中均勻分布實(shí)驗(yàn)點(diǎn)集,是一種有效的部分因子實(shí)驗(yàn)方法,特別適合多因素、多水平的試驗(yàn)設(shè)計(jì)[16]。所以本文采用4因素16水平的均勻試驗(yàn)表(采用擬水平法),另外在零點(diǎn)增加一組實(shí)驗(yàn),對(duì)總共17組實(shí)驗(yàn)進(jìn)行有限元分析。對(duì)于每組實(shí)驗(yàn)(有限元計(jì)算),計(jì)算導(dǎo)軌的X方向變形量,并與采用前述方法獲得的測(cè)量值進(jìn)行比較,計(jì)算出優(yōu)化計(jì)算指標(biāo)F(Ω)值。
采用第1節(jié)所述響應(yīng)面方法公式建立優(yōu)化計(jì)算指標(biāo)F(Ω)和設(shè)計(jì)變量 Ω之間的近似模型。得到響應(yīng)面函數(shù)表達(dá)式
根據(jù)式(2)計(jì)算得到復(fù)相關(guān)系數(shù)R2=0.9825,接近1,說明回歸方程足夠精確。
對(duì)目標(biāo)函數(shù)的二階響應(yīng)面模型方程進(jìn)行F檢驗(yàn),根據(jù)式(3)計(jì)算統(tǒng)計(jì)量F,得到F=120.5>F0.05(14,2)=19.4,說明響應(yīng)面模型是有意義的。
對(duì)于數(shù)控機(jī)床導(dǎo)軌有限元模型,為使數(shù)值模擬結(jié)果逼近實(shí)際測(cè)量結(jié)果,其優(yōu)化目標(biāo)應(yīng)盡可能減少由于有限元模型邊界條件難于精確確定所導(dǎo)致的有限元計(jì)算誤差。
綜上所述,導(dǎo)軌有限元分析模型優(yōu)化的數(shù)學(xué)模型為
這樣,有限元模型優(yōu)化問題就轉(zhuǎn)化為對(duì)上式的求解。
采用上述優(yōu)化數(shù)學(xué)模型并選擇優(yōu)化算法,對(duì)導(dǎo)軌有限元模型邊界條件進(jìn)行優(yōu)化設(shè)計(jì)。得到一組 最 優(yōu) 設(shè) 計(jì) 變 量
根據(jù)優(yōu)化后的邊界條件,利用有限元分析軟件計(jì)算導(dǎo)軌在整個(gè)行程范圍內(nèi)的熱變形量。
優(yōu)化后的導(dǎo)軌橫向(X方向)熱誤差比較曲線如圖4所示。仍選取C1、C2、C3三個(gè)測(cè)點(diǎn),將有限元計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比,計(jì)算值與實(shí)驗(yàn)結(jié)果已非常接近。對(duì)照?qǐng)D3優(yōu)化前的曲線,其最大誤差(實(shí)驗(yàn)值與計(jì)算值之間的誤差)已經(jīng)由6.4μm 降低為0.83μm,誤差平方和F(Ω)由優(yōu)化前的 813μm2降低為 11.2μm2。
在研究導(dǎo)軌熱變形誤差有限元分析的過程中,由于有限元模型的邊界條件很難從理論上精確確定,使得有限元計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果之間往往存在明顯的誤差。響應(yīng)面方法的優(yōu)勢(shì)在于可以通過較少的實(shí)驗(yàn)(有限元計(jì)算)獲得設(shè)計(jì)變量(有限元模型邊界條件)和目標(biāo)函數(shù)(熱變形計(jì)算誤差平方和)之間足夠準(zhǔn)確的相互關(guān)系,并且可以只用簡單代數(shù)表達(dá)式展現(xiàn)出來。因此,為了縮小計(jì)算誤差,本文借助響應(yīng)面方法進(jìn)行有限元模型修正。通過實(shí)驗(yàn)得到3個(gè)測(cè)點(diǎn)的熱誤差變化值,建立了優(yōu)化計(jì)算指標(biāo)(3個(gè)測(cè)點(diǎn)的計(jì)算誤差平方和)和設(shè)計(jì)變量(熱源和對(duì)流系數(shù))之間的近似模型,對(duì)導(dǎo)軌有限元邊界條件進(jìn)行了優(yōu)化選擇,使得計(jì)算誤差明顯降低。獲得了導(dǎo)軌整個(gè)行程范圍內(nèi)的接近實(shí)際測(cè)量值的熱變形動(dòng)態(tài)變化計(jì)算值,以便進(jìn)行實(shí)時(shí)誤差補(bǔ)償。但本文中的優(yōu)化選擇是在一個(gè)運(yùn)行條件下以3個(gè)測(cè)點(diǎn)為例進(jìn)行的,對(duì)于大型機(jī)床導(dǎo)軌,為了提高計(jì)算精度可以設(shè)置3個(gè)以上測(cè)點(diǎn)進(jìn)行優(yōu)化,另外,如何將近似模型應(yīng)用于復(fù)雜運(yùn)行條件(如運(yùn)行時(shí)開啟關(guān)閉機(jī)床)需要進(jìn)一步研究。
圖4 熱誤差比較曲線(優(yōu)化后)
[1] Koevoets A H,Eggink H J,van der Sanden J,et al.Optimal Sensor Configuring Techniques for the Compensation of Thermo-elastic Deformations in High-precision Systems[C]//13th International Workshop on Thermal Investigation of ICs and Systems.Budapest,2007:208-213.
[2] Yang JG,Ren Y Q,Liu G L,et al.Testing Variable Selecting and Modeling of Thermal Errors on a INDEX-G200 Turning Center[J].International Journal of Advanced Manufacturing Technology,2005,26:814-818.
[3] 楊建國,任永強(qiáng),朱衛(wèi)斌,等.數(shù)控機(jī)床熱誤差補(bǔ)償模型在線修正方法研究[J].機(jī)械工程學(xué)報(bào),2003,39(3):81-84.
[4] 閆嘉鈺,張宏韜,劉國良,等.基于灰色綜合關(guān)聯(lián)度的數(shù)控機(jī)床熱誤差測(cè)點(diǎn)優(yōu)化新方法及應(yīng)用[J].四川大學(xué)學(xué)報(bào)(工程科學(xué)版),2008,40(2):160-164.
[5] Zhao Haitao,Yang Jiangguo,Shen Jinhua.Simulation of Behavior of a CNCMachine Tool Spindle[J].International Journal of Machine Tools and Manufacture,2007,47:1003-1010.
[6] 郭勤濤,張令彌,費(fèi)慶國.結(jié)構(gòu)動(dòng)力學(xué)有限元模型修正的發(fā)展——模型確認(rèn)[J].力學(xué)進(jìn)展,2006,36(1):36-42.
[7] Carrion F J,Lozano A,Fabela M J,et al.Simplified Experimental Calibration for a Tridilosa-type Bridge Finite Element Model[J].Experimental Mechanicaics,1999,39(4):324-328.
[8] Choi J K,Lee D G.Thermal Characteristics of the Spindle Bearing System with a Gear Located on the Bearing Span[J].International Journal of Machine Tools&Manufacture,1998,38(9):1017-1030.
[9] Romero V J,Swiler L P,Giunta A A.Construction of Response Surface Based on Progressive-latticesampling Experimental Designs with Application to Uncertainty Propagation[J].Structural Safety,2004,26(2):201-219.
[10] 李郝林,應(yīng)杏娟.數(shù)控機(jī)床主軸系統(tǒng)熱誤差溫度測(cè)量點(diǎn)的最優(yōu)化設(shè)計(jì)方法[J].中國機(jī)械工程,2010,21(7):804-808.
[11] Kim T,Park H,Rhee S.Optimization of Welding Parameters for Resistance Spot Welding of TRIP Steel with Response Surface Methodology[J].International Journal of Production Research.2005,43(21):4643-4657.
[12] 席光,王志恒,王尚錦.葉輪機(jī)械氣動(dòng)優(yōu)化設(shè)計(jì)中的近似模型方法及其應(yīng)用[J].西安交通大學(xué)學(xué)報(bào),2007,41(2):125-135.
[13] 趙祖德,陳學(xué)文,陳軍,等.基于近似模型和數(shù)值模擬的連桿熱鍛成形工藝設(shè)計(jì)優(yōu)化[J].上海交通大學(xué)學(xué)報(bào),2008,42(5):748-751.
[14] 張小云.雙相鋼點(diǎn)焊熔核界面撕裂失效機(jī)理與控制方法研究[D].上海:上海交通大學(xué),2008.
[15] 王斌.機(jī)械設(shè)計(jì)手冊(cè)[M].北京:機(jī)械工業(yè)出版社,2004.
[16] Montgomery D C.Design and Analysis of Experiments[M].北京:中國統(tǒng)計(jì)出版社,1998.