趙 頤, 田曉耕

(西安交通大學(xué) 機械結(jié)構(gòu)強度與振動國家重點實驗室, 西安 710000)

0 引 言

近年來,隨著激光技術(shù)的不斷發(fā)展,其脈沖寬度已經(jīng)可以短至皮秒、飛秒量級,屬于超短脈沖的范疇．超短脈沖激光可以實現(xiàn)對硬脆材料的精密加工,在打孔、切割、劃線、沉積成型等高精度微加工領(lǐng)域有著廣大的應(yīng)用前景[1-3]．目前,對于激光加工的優(yōu)化,多通過重復(fù)實驗尋找合適的加工參數(shù)[4],存在成本高、效率低、適用范圍小的問題．通過研究硬脆材料在超短脈沖激光作用下的熱力響應(yīng),探究不同物理量在多脈沖作用下的演化規(guī)律,獲得相關(guān)參數(shù)對響應(yīng)的影響,預(yù)測加工時工件內(nèi)部的動態(tài)響應(yīng),對于優(yōu)化加工參數(shù)具有重要的指導(dǎo)意義．

對于超短脈沖激光加工熱力響應(yīng)的研究,大多在簡化真實物理過程的基礎(chǔ)上建立數(shù)學(xué)模型,利用有限元方法對物理場進行求解,這導(dǎo)致所得到的結(jié)果不能準確地描述材料的熱力響應(yīng)．Shen等[5]研究了YSZ發(fā)生固體相變過程中的動態(tài)響應(yīng)．秦淵等[6]建立了毫秒激光打孔的軸對稱模型,并對材料從發(fā)生固體相變到熔融去除過程中的孔形狀變化進行了研究．由于聚焦后較強的超短脈沖激光能量在空間上服從Gauss分布[7-8],發(fā)生錐形輻射使研究激光作用材料的溫度分布變得復(fù)雜,Rahaman等[9-10]利用三角波代替Gauss分布,并闡述了其合理性,得到了瞬態(tài)二維熱傳導(dǎo)的解析解．值得注意的是,上述研究均采用經(jīng)典Fourier傳熱定律(擴散方程描述變量以無限速度傳播),而超短時等極端條件下熱以有限速度傳播,經(jīng)典Fourier定律無法給出準確的響應(yīng)預(yù)測[11]．對于這類極端條件下的熱力耦合問題,應(yīng)采用以非Fourier傳熱為基礎(chǔ)的廣義熱彈性理論進行求解．包立平等[12]討論了激光等離子體產(chǎn)生的波模型,應(yīng)用攝動法給出了金屬板中非Fourier和Fourier溫度場的關(guān)系．譚勝等[13]基于雙相延遲模型[14]對飛秒激光燒蝕金屬的現(xiàn)象進行了數(shù)值模擬,獲得了與實驗結(jié)果一致的預(yù)測．也有學(xué)者采用分子動力學(xué)方法對固體在超短脈沖作用下的熱力響應(yīng)進行了研究,如Xiong等[15]研究了銅薄膜在超短脈沖下的熱力響應(yīng);Li等[16]研究了氬晶體在皮秒激光燒蝕過程中的相變．

現(xiàn)有研究中,多采用Fourier熱傳導(dǎo)定律,或是未考慮多次脈沖帶來的影響,對材料的熱力響應(yīng)不能準確描述．本文驗證有限元解法的可靠性后采用以C-V傳熱模型[17-18]為基礎(chǔ)的L-S廣義熱彈性理論[19]對YSZ材料在多脈沖作用下的瞬態(tài)熱力響應(yīng)開展了研究．討論了不同脈沖間隔、材料比熱容隨溫度變化等對熱力響應(yīng)的影響．最后,討論了激光脈沖作用下機械波經(jīng)邊界反射后在材料內(nèi)部的演化．

1 基 本 方 程

本文基于L-S廣義熱彈性理論,其控制方程如下．幾何方程為

(1)

本構(gòu)方程為

σij=λεkkδij+2Gεij-γΘδij;

(2)

各向同性材料不計體力的運動方程為

(3)

溫度控制方程為[20]

(4)

式中,εij為應(yīng)變分量,ui為位移分量,σij為應(yīng)力分量,ρ為材料密度,λ和G為Lamé常數(shù),δij為Kronecker指標,Θ=T-T0,T為絕對溫度,T0為參考溫度,κ為熱傳導(dǎo)系數(shù),γ=(3λ+2G)αt,αt為線性熱膨脹系數(shù),

?2為Laplace算子,τ為熱松弛時間,cE為材料比熱容,Q為熱源．

考慮超短脈沖作用時可簡化為二維軸對稱模型,在柱坐標系下改寫方程(1)、(3)、(4)可得幾何方程為

(5)

運動方程為

(6)

溫度控制方程為

(7)

式中,u和w分別為r方向和z方向的位移．對于時間與空間尺度都很小的超短脈沖問題,為了便于分析,引入以下無量綱量：

(8)

將式(8)代入式(1)—(4),可得軸對稱模型下L-S理論的無量綱控制方程．運動方程為

(9)

本構(gòu)方程為

(10)

溫度控制方程為

(11)

2 YSZ在超短脈沖激光作用下的熱力響應(yīng)

目前,求解廣義熱彈性問題常用的方法有積分變換[21-23],即借助于Laplace變換和Fourier變換將控制方程轉(zhuǎn)化至頻域內(nèi)求解后再進行積分反變換得到時域內(nèi)的解．此方法需要進行兩次積分變換,而且數(shù)值積分反變換會引入離散和截斷誤差,使得熱的波動性不能被精確捕捉[24]．還有學(xué)者采用正則模態(tài)法對一些問題進行了解析求解[25-28],以達到迅速解耦的目的．上述方法求解廣義熱彈性問題時,多限于相對簡單的物理模型與熱載荷情況,對于熱載荷形式相對復(fù)雜的問題,積分變換及解析求解均面臨著不同程度的數(shù)學(xué)困難．因此,采用有限元(FlexPDE)直接在時域?qū)刂品匠踢M行求解是一種有效的方法．

考慮超短脈沖激光輪廓迂回法打孔的工況,將打孔過程簡化為圓柱受到多次脈沖載荷作用的過程．部分材料參數(shù)如下[5,29]：

ρ=6 100 kg/m3,G=48.4 GPa,λ=48.4 GPa,

κ=2.09 W/(m·K),αt=9.4×10-6K-1,T0=294 K．

根據(jù)Vedavarz等[30]的研究,YSZ的熱松弛時間為

(12)

式中,v為固體中的聲子速度,為了簡化計算,可取固體中的聲速近似代替,本文取v=5 842 m/s．將式(12)代入式(8),可得無量綱松弛時間為

(13)

激光脈沖的能量函數(shù)為

(14)

有限元計算時,模型(如圖1)四周取絕熱邊界,其力學(xué)邊界條件上下表面應(yīng)力自由、周邊固定,可表示為

(15)

初始條件為

(16)

為簡化表述,在不引起混淆的情況下,下文將去除無量綱量的上標“～”．

圖1 輪廓迂回法打孔示意圖Fig. 1 Schematic diagram of contour roundabout punching

2.1 有限元求解準確性驗證

為驗證有限元求解的準確性,首先利用有限元方法對不計體力兩端固定且絕熱的均勻各向同性桿在移動熱源作用下的熱力響應(yīng)問題進行求解,并與文獻[20]結(jié)果進行對比．邊界條件、材料參數(shù)均與文獻中相同．由于文獻中熱源項包含Dirac函數(shù),為使有限元模擬描述的問題與文獻一致,根據(jù)Dirac函數(shù)的性質(zhì),移動熱源可表示為

(17)

式中,A為積分系數(shù),d為移動熱源在x方向的寬度,滿足A·d=0.5．本文取d=0.1,因此有A=5．f(a1,a2)為脈沖函數(shù),滿足

(18)

有限元求解時,采用自適應(yīng)網(wǎng)格,不滿足精度要求時,網(wǎng)格加密重新計算．圖2、3分別給出了熱源移動速度為2、熱松弛時間為0.05時的溫度和應(yīng)力結(jié)果．從圖中可知,有限元所得結(jié)果與解析解吻合較好,說明有限元解法可以得到準確的解答．

2.2 比熱容為常數(shù)

由于深度對熱源作用附近的應(yīng)力分布并無影響[31],本文中在豎直方向施加與z無關(guān)的熱源．首先分析比熱容為常數(shù)時(cE=453)多脈沖激光作用下YSZ的熱力響應(yīng)以及周期為6、不同脈寬(無量綱脈寬取2,3,4)對熱力響應(yīng)的影響,如圖4所示．

圖4 激光脈沖隨時間的變化Fig. 4 The variations of laser pulses with time

脈寬為2時,z=5處響應(yīng)的徑向分布如圖5(a)、5(b)所示．脈寬為2時,t=14是熱源剛消失的時刻,t=18是下一次熱源剛出現(xiàn)的時刻．由圖5(a)可知,熱松弛時間的存在導(dǎo)致熱源剛作用后系統(tǒng)的溫度并未達到一個周期內(nèi)的最大值,而是在弛豫之后達到溫度的頂峰．由圖5(b)知,在熱源激勵下,熱源內(nèi)外兩側(cè)材料均受壓,且位移隨時間不斷增大．從應(yīng)力分布曲線中可以看出,隨著脈沖次數(shù)的增多,σrr出現(xiàn)多個壓應(yīng)力峰值,而σθθ在應(yīng)力波波前出現(xiàn)一個“應(yīng)力平臺”,經(jīng)過平臺后再下降至0,整個應(yīng)力曲線的變化規(guī)律與溫度相似．

圖5 脈寬為2時的熱力響應(yīng)Fig. 5 Thermal-mechanical responses with a pulse width of 2

脈寬為4時,z=5處的溫度和應(yīng)力分布曲線如圖6所示．與脈寬為2時相似,t=16為熱源剛消失的時刻,t=18為下一次熱源剛出現(xiàn)的時刻．由圖6(a)知,由于加熱時間間隔減小,在兩次加熱之間的溫度響應(yīng)沒有明顯降低;而圖6(b)、6(c)壓應(yīng)力峰值(r=100處)均隨著時間的推移(熱源的消失及出現(xiàn))呈先增大后減小的變化．由于熱機耦合,溫度變化引起的彈性波會影響變形,導(dǎo)致應(yīng)力發(fā)生變化,應(yīng)力響應(yīng)對熱源的消失重現(xiàn)更為敏感．

在保證輸入熱量和加熱時長不變(時長為10)的條件下,相同周期不同脈寬對z=5處物理量分布的影響如圖7(a)—7(d)所示．由圖7(a)可知,輸入相同熱量的情況下由于所需時間較長,脈寬越短的熱影響區(qū)內(nèi)可達到的溫度峰值越低;由圖7(b)可知,脈寬越長位移峰值越大．同時,由于脈寬為2(較短)時脈沖次數(shù)增大(總加熱時長不變,非加熱時間變長),會導(dǎo)致位移分布曲線在下降時出現(xiàn)“臺階”．研究參數(shù)變化對應(yīng)力分布的影響,可以發(fā)現(xiàn)脈寬變短會導(dǎo)致徑向應(yīng)力σrr和周向應(yīng)力σθθ顯著減小,且由于脈沖間隔變大導(dǎo)致應(yīng)力出現(xiàn)更多的峰值．

圖7 不同參數(shù)下的熱力響應(yīng)Fig. 7 Thermal-mechanical responses under different parameters

2.3 比熱容隨溫度變化

已有研究發(fā)現(xiàn),由于非金屬主要通過晶格振動實現(xiàn)導(dǎo)熱,其導(dǎo)熱形式分為聲子導(dǎo)熱和光子導(dǎo)熱兩種．溫度較低時聲子導(dǎo)熱占主導(dǎo),溫度較高時光子導(dǎo)熱占主導(dǎo),比熱容會隨溫度變化[32-33]．考慮材料比熱容隨溫度變化,可將比熱容寫作無量綱溫度的函數(shù)：

cE=c1(a+φ(Θ))．

(19)

修改式(8)中的無量綱參數(shù)η0為

(20)

重新改寫無量綱溫度控制方程可得

(21)

若式(19)中a+φ(Θ)=1,則式(21)退化為比熱容為常數(shù)時的無量綱控制方程．

根據(jù)單水維[34]測量的結(jié)果,YSZ的比熱容隨溫度的變化可表示為如下擬合形式：

(22)

對熱松弛時間重新進行無量綱化處理可得

(23)

當(dāng)cE隨溫度變化時,考慮到a+φ(Θ)變化很小,對熱松弛時間的影響較小,可取a+φ(Θ)=1.328,得到τ=1.574．

因此,為獲得脈沖激光作用下YSZ響應(yīng)的準確預(yù)測,需考慮比熱容隨溫度變化對熱力響應(yīng)的影響,其余參數(shù)同上、仍保持為常數(shù)．為了更顯著地體現(xiàn)比熱容隨溫度的變化,取脈寬為3時不同物理量的對比結(jié)果如圖8(a)—8(d)所示．

圖8 比熱容是否變化對熱力響應(yīng)的影響Fig. 8 Effects of the specific heat capacity change on the thermal-mechanical response

脈寬為3時,t=8和11分別為加熱及不加熱的時刻．由圖8(a)可知,當(dāng)比熱容隨溫度變化時,相同的加熱條件下溫度峰值更低;從圖8(b)中可知,比熱容的變化會使得位移曲線整體降低(位移減小)．由應(yīng)力分布曲線可知比熱容隨溫度變化會導(dǎo)致徑向正應(yīng)力σrr和周向正應(yīng)力σθθ分布曲線的絕對值低于比熱容為常數(shù)時的情況．

2.4 波的反射

機械波會傳播至材料底端并發(fā)生反射,進而對已經(jīng)受到脈沖作用區(qū)域的熱力響應(yīng)造成影響．為觀察具體的傳播情況,計算域大小保持不變,縱向上熱源僅在z>5的區(qū)域作用,單次作用的無量綱時間為10,其他條件與2.2小節(jié)中保持一致．計算得到不同時刻應(yīng)力的云圖分布與z=10上的應(yīng)力分布曲線,如圖9—11所示．

圖9 應(yīng)力σrr隨時間的變化Fig. 9 Variations of stress σrrwith time

由σrr云圖隨時間的變化可知,在徑向的壓應(yīng)力會隨著時間的推移逐漸分離成兩個獨立的部分分別向圓心和遠離圓心的方向傳播．同時,σrr在縱向上經(jīng)反射后形成的拉應(yīng)力波反作用于熱源位置處,并最終形成拉壓交替的應(yīng)力波在徑向上傳播．由σθθ云圖隨時間的變化及表面處σθθ的徑向分布可知,由壓應(yīng)力波反射后產(chǎn)生的拉應(yīng)力波會沿壓應(yīng)力波的兩側(cè)傳播,直到反作用于計算域的上端,在熱源的周圍形成拉應(yīng)力,且熱源外側(cè)的拉應(yīng)力更為顯著．

圖10 應(yīng)力σθθ隨時間的變化Fig. 10 Variations of stress σθθ with time

圖11 z=10處的應(yīng)力分布Fig. 11 Stress distributions at z=10

3 結(jié) 論

本文基于L-S廣義熱彈性理論、計及材料比熱容隨溫度變化建立了激光脈沖作用下的熱彈控制方程,利用有限元方法研究了YSZ在受到多次熱沖擊下的熱力響應(yīng)．通過分析得到了以下結(jié)論：

1) 隨著脈沖作用次數(shù)的增加,σrr,σθθ均出現(xiàn)多個峰值;除位移值隨著時間一直增加外,溫度、應(yīng)力的響應(yīng)均隨著脈沖產(chǎn)生周期性變化,而應(yīng)力響應(yīng)對脈沖的消失重現(xiàn)更加敏感．

2) 在輸入能量保持不變的情況下,脈寬越短,由于熱量耗散導(dǎo)致溫度、位移、應(yīng)力的響應(yīng)峰值越?。?/p>

3) 考慮材料比熱容隨溫度變化,溫度、位移和應(yīng)力響應(yīng)均降低,表明考慮比熱容與溫度的相關(guān)性對響應(yīng)的準確預(yù)測至關(guān)重要．

4) 周向應(yīng)力σθθ在邊界處的反射會在熱源兩側(cè)產(chǎn)生拉應(yīng)力,且外側(cè)更為顯著;徑向應(yīng)力σrr在傳播中沿徑向上發(fā)生分離,受反射波的作用最終形成拉壓交替的應(yīng)力波．隨應(yīng)力波的傳播,材料則受到拉壓交替的疲勞載荷作用．