基于有限元的脈沖激光輻照材料溫度場研究

馬 健，趙 揚，周鳳艷，孫繼華，劉 帥，郭 銳，宋江峰，陳建偉

(1．山東省科學(xué)院激光研究所，山東 濟南250014;2．山東省汽車工業(yè)集團有限公司，山東 濟南250011)

1 引言

作為激光超聲的激勵源，脈沖激光能夠以較大的峰值功率輻照于材料表面，一般認為激光超聲波的產(chǎn)生基于熱彈和燒蝕兩種機理，熱彈與燒蝕的主要區(qū)別為瞬間加載于物體表面上的激光能量是否引起表面熔化、氣化。當激光的功率密度達到106～108W/cm2時，材料表面會發(fā)生氣化。上述數(shù)值僅是一個范圍，對于給定的激光功率密度，材料的溫升還與吸收率、材料參數(shù)密切相關(guān)，能否引起材料的燒蝕卻不好判斷，這就需要建立一套有效的數(shù)值計算方法，用于求解脈沖激光輻照材料的溫度場。

采用有限元分析材料受激光輻照的溫度場是常用的熱分析方法，鑒于目前有限元熱分析文獻中只是定性說明在激光輻照區(qū)域中網(wǎng)格劃分要密集一些，或者直接定量指明網(wǎng)格的大小，以便得到精確的結(jié)果;而在其他區(qū)域單元網(wǎng)格尺寸取大一些，以減少計算量節(jié)省計算機資源［1－2］，但是對于具體的選擇依據(jù)沒有做詳細明確的表述。本文討論了激光輻照材料表層的熱分析時，單元網(wǎng)格的大小對有限元分析結(jié)果的影響，給出了以熱擴散深度為依據(jù)的網(wǎng)格大小選取規(guī)則。最后，基于有限元熱分析的溫度數(shù)據(jù)，給出了應(yīng)力場的計算結(jié)果。

2 理論基礎(chǔ)與數(shù)值計算

2．1 理論模型

圖1所示為材料受激光輻照的模型，x，y，z分別表示直角坐標系中的三個不同方向上的坐標，單位m。受激光輻照的物體為圓柱體，其對稱軸為y軸，輻照中心與坐標系原點重合。激光器發(fā)射出重復(fù)頻率5 Hz直徑為6 mm的脈沖激光，脈沖寬度為10 ns，波長為1064 nm。輻照區(qū)域內(nèi)脈沖激光功率密度的表達式為:

式中，r，τ分別為點距光斑中心的距離與時間，單位相應(yīng)為m，s，此時為光斑中心的功率密度;A為材料表面的吸收率;f(r)和g(τ)分別為激光功率密度的空間分布和時間分布。事實上，激光輻照材料表層的過程極其復(fù)雜，用有限元法完全模擬出激光輻照的實際情況是不可能的，為了盡可能的反應(yīng)實際情況并有效地簡化問題，現(xiàn)提出以下假設(shè):①環(huán)境溫度為300 K。②材料的導(dǎo)熱性能是均勻的而且各向同性。③忽略固相之間的轉(zhuǎn)換對溫升規(guī)律的影響。④將材料與工作臺接觸面上的散熱量忽略，其余各面均存在對流換熱與輻射換熱，引入總換熱系數(shù)h進行分析［3］。

圖1 材料受激光輻照的理論模型Fig．1 Model for material surface irradiated by laser

2．2 導(dǎo)熱微分方程及邊界條件

獲得材料受激光輻照的溫度場分布，需要根據(jù)能量守恒定律和傅里葉定律來建立溫度場的導(dǎo)熱微分方程，無內(nèi)熱源的三維非穩(wěn)態(tài)導(dǎo)熱微分方程［4］如下:

式中，T為溫度，K;λ為導(dǎo)熱系數(shù)，W/(m·K);ρ為材料的密度，kg/m3;c為材料的比熱，J/(kg·K)?？紤]材料的熱物性時，λ、ρ、c均為變量，否則λ、ρ、c均為常量。其余參數(shù)含義同上。

對于式(2)，為獲得方程的唯一解，需要附加一定的邊界條件和初始條件，根據(jù)前提條件及假設(shè)，在脈沖激光輻照時(τ0≤τ≤τ0+τp)，材料y=0表面上激光輻照區(qū)域中的單元存在著熱流載荷［5］:

除了材料與工作臺相接觸的面，其余各表面存在著對流與輻射換熱，換熱邊界條件為:

式(3)、(4)中，τ0為每個脈沖激光的起始時刻;τp為脈沖寬度，本文中τp=10 ns;n表示表面的法向矢量;Γ是材料的邊界;T0是環(huán)境溫度。

初始條件:T0=300 K。

2．3 激光輻照材料溫度場的有限元分析

脈沖激光的時間分布g(τ)取為矩形，空間分布f(r)為高斯平頂分布，經(jīng)由點聚焦鏡透射并通過遮擋去除邊緣部分之后，在材料表面形成直徑D=1．8 mm、能量分布均勻的圓形光斑，圓形光斑具有的單脈沖能量為10．2 mJ，經(jīng)過計算I(r，τ)=1×107W/cm2。首先以航空航天中常用的熱障涂層系統(tǒng)為研究對象，對單次脈沖輻照時的情況進行有限元分析和數(shù)學(xué)解析。其中表面陶瓷層的厚度為200μm，直徑為1．8 mm。為了便于有限元求解結(jié)果與理論解析值的對比，暫不考慮材料的熱物性，忽略熱障涂層各成分界面上的熱交換，并將A、h視為恒定值，大小分別為0．25、0。表1中為表面陶瓷層的材料參數(shù)［6］。通過對比有限元求解結(jié)果和溫度場理論值，驗證輻照材料表層溫度場有限元求解方法的正確性，并得出輻照區(qū)域網(wǎng)格尺寸的劃分依據(jù)。

表1 陶瓷層的參數(shù)Tab．1 Parameters of Ceramic Top Coat

在脈沖激光作用完畢的時刻，熱擴散深度為［7］:

式(5)中，D為熱擴散率，m2/s，其計算方法為:

根據(jù)式(5)求出熱擴散深度為1．6×10－7m。

首先利用有限元分析軟件對陶瓷層進行瞬態(tài)熱分析:①建立2．1節(jié)所述的理論模型并分配材料參數(shù)，根據(jù)模型的軸對稱性選擇四節(jié)點軸對稱單元，將三維熱分析問題轉(zhuǎn)換為軸對稱問題進行處理。為獲得不同網(wǎng)格劃分下的溫度場計算結(jié)果，分別按照脈沖激光結(jié)束時熱擴散深度的1倍、1/2倍、1/4倍、1/8倍、1/10倍離散模型。②選擇分析類型為瞬態(tài)熱分析，設(shè)置環(huán)境溫度為300 K，總的分析時間長度為1μs。為了便于處理，時間軸零點取為單個脈沖的起始時刻。0～10 ns對材料受輻照區(qū)域施加熱流密度載荷，為了獲得精確數(shù)據(jù)，時間步長取為1 ns;隨后卸載載荷，在11 ns～0．2μs區(qū)域中，時間步長取為2 ns。而在其余時間內(nèi)，為了提高計算效率，時間步長取為20 ns。③直接求解獲得溫度場，求解完畢后進行時間歷程分析，提取節(jié)點溫度隨時間的變化規(guī)律。該有限元分析所用計算機的主要配置為8核I7－3770 CPU處理器，16GB內(nèi)存。將各網(wǎng)格劃分下獲得的節(jié)點溫度最大值列于表2中。圖2為h/4網(wǎng)格劃分下，有限元求解出的脈沖激光輻照陶瓷層的溫升曲線。

其次對陶瓷層的溫度場進行理論解析，以驗證有限元熱分析結(jié)果的正確性。當垂直入射的激光束覆蓋了整個y=0表面，以及熱擴散深度遠小于光斑大小以及陶瓷層的厚度時，輻照表面中心下方區(qū)域的溫度計算可以簡化為一維處理，則激光加熱和冷卻過程中的瞬態(tài)溫度場解析公式為:

式中各參數(shù)的含義均同上。

根據(jù)式(7)和已知參數(shù)，利用MATLAB編寫程序，計算輻照表面中心及其正下方位置的溫度場，得到表層圖3中的曲線和溫度最大值為1557 K。對比圖2和圖3，溫度場的有限元分析結(jié)果和理論解析結(jié)果均呈現(xiàn)以下規(guī)律:表面中心點的溫度在脈沖結(jié)束時刻達到最大值，然后溫度緩慢下降;而表面下方各處的溫度升至最高點的時刻逐漸滯后，溫度最大值也隨深度的增加而減小。通過對比可知，溫度場有限元結(jié)果與理論解析結(jié)果有較好的吻合。

表2 有限元計算結(jié)果Tab．2 Data of finite element calculation

表2同時表明隨著網(wǎng)格尺寸的減小，輻照中心節(jié)點的溫度最大值不斷增大，有限元分析結(jié)果逐漸逼近理論值。但所需的計算時間和產(chǎn)生的數(shù)據(jù)量顯著增加，而相對誤差減小的幅度卻逐漸降低。當網(wǎng)格尺寸取熱擴散深度的1/10時，有限元結(jié)果和理論值之間的相對誤差為1．24%，單純追求有限元求解結(jié)果與理論值的逼近性，而大量增加計算機資源的消耗并不可取。一般情況下，可以取網(wǎng)格尺寸為脈沖激光結(jié)束時熱擴散深度的1/4，此時，有限元結(jié)果和理論值之間的相對誤差為2．38%。

圖2 輻照中心及下方節(jié)點的溫升曲線的有限元計算結(jié)果Fig．2 Temperature curve of irradiated center and below nodes for finite element results

圖3 輻照中心及下方節(jié)點的理論溫升曲線Fig．3 Temperature curve of irradiated center and below nodes for theory results

3 有限元網(wǎng)格劃分依據(jù)及應(yīng)力場分析

由上可知，有限元求解脈沖激光輻照材料的溫度場時，應(yīng)以脈沖激光結(jié)束時的熱擴散深度作為依據(jù)，綜合考慮計算精度和所需消耗的計算機資源，選取受輻照區(qū)域的網(wǎng)格大小。對于熱物性比較明顯的材料，在脈沖激光結(jié)束時刻輻照區(qū)域的溫度場未知，與溫度相關(guān)的材料參數(shù)值無法具體選取，相應(yīng)的熱擴散深度亦無法具體求出。此種情況下，則根據(jù)各溫度點求出脈沖激光結(jié)束時熱擴散深度的最小值作為網(wǎng)格尺寸選取的依據(jù)。

相對于理論解析，利用有限元可方便求解復(fù)雜的問題，能夠同時考慮材料熱物性、對流與輻射換熱、材料的相變潛熱等多種因素?，F(xiàn)給出一個具體分析實例:脈沖激光的參數(shù)同上，選取幾何尺寸為直徑20 mm，厚度10 mm鋁質(zhì)工件進行研究，鋁的熱物性參數(shù)［8］如表3所示。

表3 鋁的熱物性參數(shù)Tab．3 Thermal parameters of aluminum(Melting point:933．4 K;boiling point:2740 K)

根據(jù)表3和式(5)求出熱擴散深度的最小值為1．52×10－6m，取有限元模型網(wǎng)格尺寸為h/4(3．8×10－7m)，依據(jù)建立的溫度場有限元求解方法，計算鋁質(zhì)試塊受脈沖激光輻照的溫度場。圖4為10 ns時刻鋁質(zhì)試塊的軸截面右側(cè)局部溫度場，圖5為輻照區(qū)域中心節(jié)點及其下方節(jié)點的溫升曲線。

圖4 10 ns時刻鋁質(zhì)試塊的右側(cè)局部溫度場Fig．4 Partial temperature field of the right side of aluminum test block at 10 ns

圖5 鋁質(zhì)試塊的輻照中心節(jié)點及下方節(jié)點的溫升曲線Fig．5 Temperature curve of aluminum test block for finite element results

由求解的溫度場數(shù)據(jù)可以看出，輻照區(qū)域的最高溫度為415．8 K，低于鋁的熔點933．4 K。10μs時刻輻照中心點的溫度已降至300．1 K，而5 Hz重復(fù)頻率對應(yīng)的脈沖間隔為200 ms。相比之下，鋁質(zhì)工件輻照區(qū)域降溫所消耗的時間非常短，因此鋁質(zhì)工件輻照區(qū)域的溫度數(shù)值可視為無累積效應(yīng)。所以激光超聲檢測鋁質(zhì)工件時，加載脈沖寬度為10 ns，單脈沖能量為10．2 mJ，重復(fù)頻率5 Hz的激光是允許的。

工件的局部的溫度場發(fā)生變化，必然引起應(yīng)力場的變化，進而導(dǎo)致超聲波的產(chǎn)生。根據(jù)熱障涂層和鋁質(zhì)工件溫升的數(shù)值計算結(jié)果，二者輻照區(qū)域的最高溫度均低于相應(yīng)的熔點，所以激光超聲波的產(chǎn)生符合熱彈機理。忽略位移場對溫度場的影響，采用耦合場分析中的間接法，將節(jié)點溫度分析結(jié)果作為體載荷施加到結(jié)構(gòu)應(yīng)力分析中，可獲得節(jié)點豎直方向上位移的時間歷程曲線，圖6和圖7分別為200μm厚、直徑20 mm熱障涂層和10 mm厚、直徑20 mm鋁質(zhì)試塊的上表面節(jié)點位移隨時間的變化，與學(xué)者沈中華［9］和王紀俊［10］得到的結(jié)果一致。

圖6 熱障涂層上表面節(jié)點豎直方向上的位移隨時間的變化Fig．6 Displacement of the node for thermal barrier coating system

圖7 鋁質(zhì)試塊上表面節(jié)點豎直方向上的位移隨時間的變化Fig．7 Displacement of the node for aluminum test block

4 結(jié)論

建立了材料受脈沖激光輻照的理論模型，針對輻照區(qū)域的溫度場提出了有效的有限元求解方法，討論了有限元計算中單元網(wǎng)格尺寸對計算結(jié)果的影響。在利用有限元求解脈沖激光輻照區(qū)域的溫度場時，應(yīng)以脈沖激光結(jié)束時的熱擴散深度為依據(jù)，綜合考慮計算精度和所需消耗的計算機資源等因素，確定網(wǎng)格尺寸。一般情況下，將受輻照區(qū)域中的網(wǎng)格大小取為脈沖激光結(jié)束時熱擴散深度的1/4，這樣既能夠獲得可靠的數(shù)據(jù)，又能降低對計算機資源的消耗。材料溫升的有限元數(shù)值計算結(jié)果能夠為激光能量的加載提供參考，并能夠為后續(xù)應(yīng)力場的瞬態(tài)分析提供體載荷。

［1］ YUAN Liguo，ZHAO Guomin，JIAO Luguang．Preliminary study on the temperature history of metal/liquid structure under pulsed laser irradiation［J］．Laser＆Infrared，2010(1):18－21．(in Chinese)袁立國，趙國民，焦路光．脈沖激光輻照下金屬/液體結(jié)構(gòu)溫度變化初探［J］．激光與紅外，2010(1):18－21．

［2］ AN Fuping，CHEN Zhijun，YAO Jianhua，et al．Numerical simulation of temperature field on precipitation hardening SSafter laser solid solution and alloying synchronous hybrid strengthening［J］．Applied Laser，2013，33(002):119－124．(in Chinese)安福平，陳智君，姚建華，等．沉淀硬化不銹鋼激光固溶合金化同步復(fù)合強化溫度場的數(shù)值模擬［J］．應(yīng)用激光，2013，33(002):119－124．

［3］ LIU Zhuang，WU Zhaoji，WU Jingzhi，et al．Numerical simulation of the heat treatment process［M］．Beijing:Science Press，1996:7－9(in Chinese)．劉莊，吳肇基，吳景之，等．熱處理過程的數(shù)值模擬［M］．北京:科學(xué)出版社，1996:7－9．

［4］ YANG Shiming，TAO Wenquan．Heat transfer［M］．4th ed．Beijing:High Education Press，2006:42．(in Chinese)楊世銘，陶文銓．傳熱學(xué)［M］．第4版．北京:高等教育出版社，2006:42．

［5］ WU Sijie，ZHAO Xiaobei，YANG Dongsheng，et al．Laser damage in IR detector［J］．Infrared and Laser Engineering，2013，42(5):1184－1188．(in Chinese)吳思捷，趙曉蓓，楊東升，等．激光輻照對紅外探測器的損 傷［J］．紅 外 與 激 光 工 程，2013，42(5):1184－1188．

［6］ BAI Yu-mei，XU Ying-qiang，LAI Ming-rong，et al．Analysis of Residual Stresses in Thermal Barrier Coating due to Thermal Mismatch［J］．Science Technology and Engineering，2011，11(14):3126－3129．(in Chinese)白玉梅，徐穎強，賴明榮，等．熱障涂層熱不匹配殘余應(yīng)力的分析研究［J］．科學(xué)技術(shù)與工程，2011，11(14):3126－3129．

［7］ SUN Cheng-wei．Laser Irradiation Effect［M］．Beijing:National Defense Industry Press，2002．1:4－66．(in Chinese)孫承緯．激光輻照效應(yīng)［M］．北京:國防工業(yè)出版社，2002．1:4－66．

［8］ HE Yuejuan，CHEN Guoqing，SHEN Zhonghua，et al．The influence of the line-source laser's width on laser-induced temperature field in aluminum pipe［J］．Laser Journal，2008，29(1):68－69．(in Chinese)何躍娟，陳國慶，沈中華，等．激光線源寬對鋁管中溫度場的影響［J］．激光雜志，2008，29(1):68－69．

［9］ SHEN Zhonghua，XU Baiqiang，NI Xiaowu，et al．Numerical simulation of laser-generated ultrasonic waves in layered plates［J］．Journal of Physics D:Applied Physics，2004，37(17):2364－2370．

［10］WANG Jijun，SHEN Zhonghua，NI Xiaowu，et al．Numerical simulation of laser-generated surface acoustic waves in the transparent coating on a substrate by the finite element method［J］．Optics＆Laser Technology，2007，39(1):21－28．