時域間斷Galerkin有限元法在激光熱加工過程中應(yīng)用

吳 志 剛， 郭 攀， 武 文 華＊

（1.大連理工大學(xué) 航空航天學(xué)院，遼寧 大連 116024；2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116024）

0 引 言

在工程實際中，由于激光的定向性、單色性好，發(fā)散角小的特點，激光技術(shù)具有很好的局部性以及能量的精確性.目前，激光熱加工、熱處理技術(shù)已經(jīng)得到了廣泛應(yīng)用.當(dāng)利用激光束加熱固體介質(zhì)時，在介質(zhì)表面和內(nèi)部出現(xiàn)光波吸收、反射、折射的效應(yīng)，但其中大部分自由電子的能量通過電子與晶格或離子的相互作用轉(zhuǎn)化為介質(zhì)熱能并進行熱擴散，引起固體介質(zhì)內(nèi)部不同程度的溫度上升[1].常規(guī)條件下，介質(zhì)內(nèi)部的溫度場可以通過經(jīng)典的傅里葉定律理論進行描述.但隨著超短脈沖激光、激光作用下金屬快速凝固等技術(shù)的應(yīng)用，經(jīng)典傅里葉定律存在著明顯的缺陷，熱傳導(dǎo)中非傅里葉效應(yīng)的影響得到人們的重視[2].傅里葉定律理論認(rèn)為熱流與溫度梯度成正比，溫度場可以由拋物線型熱傳導(dǎo)方程來描述.但當(dāng)熱作用瞬時時間達(dá)到皮秒級、飛秒級，熱擴散的速度就有了奇異性.就必須要考慮介質(zhì)熱傳導(dǎo)過程中的松弛行為，即熱流矢量和溫度梯度間存在了時間延遲.現(xiàn)在，已經(jīng)有大量數(shù)學(xué)模型對這種現(xiàn)象進行描述.如單相延遲模型、雙向延遲模型、微觀兩步模型、純聲子散射模型等[3、4].其中單相延遲雙曲型熱傳導(dǎo)模型在工程中應(yīng)用得最為廣泛.

目前主要使用差分法、有限元法等數(shù)值方法進行考慮非傅里葉效應(yīng)的數(shù)值模擬.在時域處理上往往利用差分法對問題進行求解.為了避免數(shù)值振蕩，時域差分法對時間步長的要求非常嚴(yán)格[5].本文應(yīng)用時域間斷Galerkin有限元法針對受激光熱源作用的半無限體和薄膜溫度場分布問題進行求解，并使用算例對結(jié)果進行驗證.

1 單相延遲雙曲型熱傳導(dǎo)模型與時域間斷Galerkin方法

單相延遲雙曲型熱傳導(dǎo)模型中與時間相關(guān)的函數(shù)為

式中：q為熱流矢量，tk為松弛時間，k是導(dǎo)熱系數(shù)，T是溫度[6].

考慮瞬態(tài)能量守恒方程：

式中：ρ是介質(zhì)密度，cp是定壓比熱容，g是所施加的激光熱源.其函數(shù)形式為

式中：I（t）是激光熱源強度，R為介質(zhì)表面反射率，μ是介質(zhì)吸收系數(shù).將式（2）、（3）代入式（1）中，得到雙曲型激光熱傳導(dǎo)方程：

式中：a是熱擴散系數(shù)是熱波動速度.

對上式進行量綱一化處理，表達(dá)式如下：

式中：θ是量綱一溫度，τ是量綱一時間，X是量綱一長度，ψ0是量綱一內(nèi)熱源系數(shù)，η是量綱一熱源，β是量綱一熱吸收系數(shù).

半無限問題和薄膜問題相應(yīng)的量綱一化后初始條件表示如下：

時域間斷 Galerkin方法[7、8]是在傳統(tǒng)的空間域有限元方法分析基礎(chǔ)上，區(qū)別于習(xí)慣對時域（0，T）利用連續(xù)差分法求解，對時域進行有限元間斷插值離散.

雙曲型熱傳導(dǎo)方程空間域離散得

其中

時域離散可表示為

時域離散時允許基本未知函數(shù)θ與其導(dǎo)數(shù)ν在離散點處間斷.即在時刻τn的未知函數(shù)值表示為

在任意時間步In＝ （τn，τn＋1）內(nèi)，對時域基本未知函數(shù)θ采用三次Hermite插值，對時域基本未知函數(shù)導(dǎo)數(shù)ν和熱源Q采用線性插值.

對溫度函數(shù)θ與其導(dǎo)數(shù)ν作為獨立的變量進行變分，利用式（6）、（7）并選取控制條件θ·－ν＝0構(gòu)造典型的間斷Galerkin有限元法的弱形式，表示為

把式（8）～（11）代入上式，可以進一步轉(zhuǎn)化為解耦的式（12），即時域間斷Galerkin有限元法的基本求解公式.從式中不難發(fā)現(xiàn)溫度向量在時域內(nèi)間斷點處不再存在間斷，僅保留溫度時間導(dǎo)數(shù)在時間間斷點處存在間斷，大大減少了基本未知數(shù)的個數(shù)，節(jié)省了求解時間.

第n＋1時刻的溫度為

2 算例和結(jié)果

本文選取工程中常用的幾類激光熱源η形式（式（13）～（15）），利用時域間斷 Galerkin有限元法分別對半無限體一端和薄膜兩端受激光熱源作用進行數(shù)值分析.

3類熱源表達(dá)式如下：

算例1 考慮一個半無限體的一維問題，如圖1所示[4、9].半無限體左端分別受式（13）～（15）的3類激光熱源作用，初始條件為

假定ψ0＝1，β＝1.

圖1 半無限體左側(cè)受激光熱源作用模型Fig.1 Model of semi－infinite body subjected to laser heat source in its left side

在計算中，選取一維單元進行分析，單元長度為0.05，單元數(shù)為200，節(jié)點數(shù)為201，時間步長為0.000 75.

首先考慮熱源為式（13）時的數(shù)值模擬，可求得沿軸方向上在1、3、6、9s時溫度分布，如圖2所示.同時給出數(shù)值解與解析解對比，其良好的吻合程度顯示出計算結(jié)果的可靠性.

然后考慮第二類激光熱源（式（14）），選取矩形四節(jié)點單元作為分析對象，單元的尺寸為0.05×1，單元數(shù)為600，節(jié)點數(shù)為804，時間步長為0.05.與第一類激光熱源作用選取的時間步長相比，第二類激光熱源作用所選取的時間步長明顯增大.當(dāng)熱源滿足式（14）時x、y長度方向上溫度分布如圖3所示.

為了驗證方法的適用性，在式（15）作為激光熱源的工況下，選取3種不同尺寸的矩形四節(jié)點單元進行計算.三類單元尺寸分別為0.05×1、0.02×1、0.01×1，如圖4所示.單元數(shù)為450，節(jié)點數(shù)為604，時間步長為0.05.可以求得x、y長度方向上溫度分布如圖5所示.需要說明的是，在第二和第三類激光熱源的工況下，也進行了一維單元的計算，其計算結(jié)果與二維計算結(jié)果相同.

圖2 算例1第一類激光熱源在1、3、6、9s時溫度分布Fig.2 Temperature distributions at 1，3，6，9 second under 1st type of laser heat source in the first example

圖3 算例1第二類激光熱源1、3、6、9s時溫度分布Fig.3 Temperature distributions at 1，3，6，9 second under 2nd type of laser heat source in the first example

圖4 算例1第三類激光熱源計算構(gòu)型的網(wǎng)格劃分Fig.4 Computational configuration in 3rd laser heat source in the first example

圖5 算例1第三類激光熱源在1、3、6、9s時溫度分布Fig.5 Temperature distributions at 1，3，6，9 second under 3rdtype of laser heat source in the first example

算例2 激光加熱薄膜過程中，薄膜內(nèi)將產(chǎn)生能量沉積，光能轉(zhuǎn)化為熱能，從而導(dǎo)致薄膜破壞.本算例考慮薄膜兩端都受有激光熱源作用的熱傳導(dǎo)問題[1、9].模型如圖6所示，與算例1類似，薄膜兩端分別考慮3類不同的激光熱源作用，同時假定薄膜兩端的初始條件和邊界條件對稱，即初始條件分別如下：

并考慮ψ0＝1，β＝1.

圖6 薄膜激光熱源模型Fig.6 Model of thin film with laser heat source

在計算中，全部選取一維單元，單元長度為0.05，單元數(shù)為200，節(jié)點數(shù)為201，時間步長為0.000 75，薄膜厚度為10個單位.

當(dāng)熱源滿足式（13）時求得薄膜長度方向的溫度分布如圖7所示.

當(dāng)熱源滿足式（14）時同樣可以求得薄膜長度方向的溫度分布如圖8所示.

當(dāng)熱源滿足式（15）時可以求得薄膜長度方向的溫度分布如圖9所示.

圖7 算例2第一類激光熱源下1、3、6、9 s時溫度分布Fig.7 Temperature distributions at 1，3，6，9 second under 1st type of laser heat source in the second example

圖8 算例2第二類激光熱源下1、3、6、9s時溫度分布Fig.8 Temperature distributions at 1，3，6，9 second under 2ndtype of laser heat source in the second example

圖9 算例2第三類激光熱源下1、3、6、10s時溫度分布Fig.9 Temperature distributions at 1，3，6，10 second under 3rd type of laser heat source in the second example

將計算結(jié)果與文獻(xiàn)[9、10]相比，可以看出利用時域間斷Galerkin有限元法所得到的計算結(jié)果與解析解十分吻合.

3 結(jié) 語

本文利用時域間斷Galerkin有限元法對單相延遲的雙曲型激光熱源作用的熱傳導(dǎo)過程進行了數(shù)值模擬，分別考慮了半無限體和薄膜兩種情況下的熱傳導(dǎo)過程.

計算結(jié)果顯示出，由于時域間斷Galerkin有限元法在時域上對時間函數(shù)及其導(dǎo)數(shù)進行了一階到三階的插值，并且允許時間函數(shù)和其導(dǎo)數(shù)在時間節(jié)點上間斷，保證了時間函數(shù)連續(xù)的同時，有效避免了在計算激光熱傳導(dǎo)時對時間函數(shù)進行差分而帶來的數(shù)值振蕩現(xiàn)象.時域間斷Galerkin有限元法是分析瞬態(tài)激光非傅里葉熱傳導(dǎo)問題的有效方法.

[1] 孫承偉.激光輻照效應(yīng)[M].北京：國防工業(yè)出版社，2002

[2] 劉 靜.微米/納米尺度傳熱學(xué)[M].北京：科學(xué)出版社，2001

[3] 蔣方明，劉登瀛.非傅立葉導(dǎo)熱的最新研究進展[J].力學(xué)進展，2002，32（1）：128－140

[4] M.VON奧爾曼.激光束與材料相互作用的物理原理及應(yīng)用[M].北京：科學(xué)出版社，1994

[5] 孔祥謙.熱應(yīng)力有限元法分析[M].上海：上海交通大學(xué)出版社，1999

[6] TAMMA K K，NAMBURU R R.Computational approaches with applications to non－classical and classical thermo mechanical problems [J].Applied Mechanics Reviews，1997，50（9）：514－551

[7] 李錫夔，姚冬梅.彈塑性體中波傳播問題的間斷Galerkin有限元方法[J].固體力學(xué)學(xué)報，2003，24（24）：399－409

[8] 武文華，李錫夔.固體非傅里葉溫度場的間斷Galerkin有限元方法[J].計算力學(xué)學(xué)報，2007，24（2）：219－223

[9] LEWANDOWSKA M.Hyperbolic heat conduction in the semi－infinite body with a time dependent laser heat source [J].Heat and Mass Transfer，2001，37（4－5）：333－342

[10] SHUICHI Torii，YANG Wen－jie. Heat transfer mechanisms in thin film with laser heat source[J].Heat and Mass Transfer，2005，48（3－4）：537－544