基于L-S廣義熱彈性理論YSZ在超短脈沖下的熱力響應

趙 頤, 田曉耕

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

0 引 言

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

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

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

1 基 本 方 程

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

(1)

本構方程為

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

(2)

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

(3)

溫度控制方程為[20]

(4)

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

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

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

(5)

運動方程為

(6)

溫度控制方程為

(7)

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

(8)

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

(9)

本構方程為

(10)

溫度控制方程為

(11)

2 YSZ在超短脈沖激光作用下的熱力響應

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

考慮超短脈沖激光輪廓迂回法打孔的工況,將打孔過程簡化為圓柱受到多次脈沖載荷作用的過程．部分材料參數如下[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．

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

(12)

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

(13)

激光脈沖的能量函數為

(14)

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

(15)

初始條件為

(16)

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

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

2.1 有限元求解準確性驗證

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

(17)

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

(18)

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

2.2 比熱容為常數

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

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

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

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

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

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

圖7 不同參數下的熱力響應Fig. 7 Thermal-mechanical responses under different parameters

2.3 比熱容隨溫度變化

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

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

(19)

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

(20)

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

(21)

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

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

(22)

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

(23)

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

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

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

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

2.4 波的反射

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

圖9 應力σrr隨時間的變化Fig. 9 Variations of stress σrrwith time

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

圖10 應力σθθ隨時間的變化Fig. 10 Variations of stress σθθ with time

圖11 z=10處的應力分布Fig. 11 Stress distributions at z=10

3 結 論

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

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

2) 在輸入能量保持不變的情況下,脈寬越短,由于熱量耗散導致溫度、位移、應力的響應峰值越小．

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

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