龔艷春，武文遠，黃雁華，戴斌飛，吳成國

（解放軍理工大學(xué)理學(xué)院，江蘇南京 211101）

激光輻照下固體熱輸運的隨機算法及應(yīng)用

龔艷春，武文遠，黃雁華，戴斌飛，吳成國

（解放軍理工大學(xué)理學(xué)院，江蘇南京 211101）

采用隨機算法對激光輻照下固體目標溫度場的分布進行了計算機模擬。該算法將熱輸運過程作為動力學(xué)問題來對待，不考慮具體的傳輸機制，從而詳盡地考慮系統(tǒng)的細節(jié)，相比于傳統(tǒng)的數(shù)值算法，該算法簡單明了，占用計算機資源少，普適性強，尤其在處理復(fù)雜系統(tǒng)問題時具有獨特的優(yōu)點。

激光輻照；熱輸運；統(tǒng)計算法

隨著激光技術(shù)的廣泛應(yīng)用，尤其是激光武器的出現(xiàn)，激光損傷機制和激光防護技術(shù)越來越受到人們的關(guān)注。研究表明激光束入射到目標體上，破壞目標體的直接因素是激光加熱引起的溫升，因此研究目標體的激光破壞，必須考慮目標體在激光輻照下溫度場的時間和空間分布，從而為激光防護提供依據(jù)。激光輻照下目標體的溫升取決于目標材料的性質(zhì)，包括激光吸收系數(shù)、激光反射系數(shù)、熱容率、熱傳導(dǎo)率等，而這些性質(zhì)又與材料的微結(jié)構(gòu)和表面狀況存在密切的聯(lián)系，因此研究激光的輻照效應(yīng)與目標材料特征參數(shù)之間的關(guān)系，有助于人們尋求新的激光防護材料以及采取新的激光防護措施［1－2］。

關(guān)于激光輻照下目標物體溫度場的計算，傳統(tǒng)的方法是給定激光加載條件（通常把吸收的激光能量或光強作為面熱源（表層吸收）或體熱源（深層吸收）處理）、目標邊界或界面條件、初始溫度分布以及物體的熱物理特性，將激光加熱問題歸結(jié)為不定常熱傳導(dǎo)方程的解析描述，通過有限差分（finite-difference）或有限元（finite-unite）的方法進行數(shù)值計算［3］。不過，目前大多數(shù)詳細的計算方法都是針對簡單系統(tǒng)設(shè)計的，只能處理幾何形狀和加熱方式都比較簡單的問題［4－5］。這些方法在處理復(fù)雜系統(tǒng)，比如復(fù)雜的材料結(jié)構(gòu)、多種材料組分、材料特征參數(shù)對溫度具有依賴性、空間分布不規(guī)則的擴展光源等問題時，往往需要進行各種各樣的近似和簡化，這些都將影響模擬計算結(jié)果的準確性。有些物理問題，比如納米尺度上的熱傳輸問題，由于很難進行實驗測量，計算機模擬就成為基本的研究方法，而計算精度的不確定性將帶來很大問題。若將一個物理問題的各種復(fù)雜因素均考慮在內(nèi)，又會使數(shù)值計算變得很難處理。對于微觀機制不明以致不能寫出解析表達式的熱過程的模擬計算，傳統(tǒng)方法更是顯得無能為力。

闡述了一種新的隨機算法（stochastic methods），這種算法將熱輸運過程作為動力學(xué)問題來對待，不必考慮具體的傳輸機制，而是用統(tǒng)計的方法進行模擬。這種算法詳盡地考慮系統(tǒng)的細節(jié)，而計算過程又相當簡化，可以用來處理各種各樣的熱過程，因而具有普適性，尤其在處理復(fù)雜系統(tǒng)問題時具有獨特的優(yōu)點。筆者采用這種方法對固體硅在激光作用下的熱傳輸過程進行了詳細的計算，得出了一些有益的結(jié)論。

1 計算方法

隨機方法最初是由BUNKER和GILLESPIE提出的，用以模擬復(fù)雜的化學(xué)反應(yīng)，他們可以快速準確地獲得一組微分方程的解，比傳統(tǒng)的數(shù)值模擬方法要簡捷得多。HOULE最先將這種算法用于固體材料在激光照射下的溫度變化過程的模擬［6］，其基本思想是將熱傳輸過程看成是一個隨機的動力學(xué)過程，這一點與Monte Carlo方法類似。

人們知道激光輻照目標時，激光輻射能在目標表面處一部分被反射，一部分折射進入目標體內(nèi)，使目標材料被加熱、熔解以及蒸發(fā)，目標上的溫度分布隨著時間的變化可以用熱輸運方程來描述：

式中：T是目標體的溫度，它是時間t和空間位置的函數(shù)；C（T），ρ（T）以及K（T）分別是目標體的熱容、質(zhì)量密度以及熱傳導(dǎo)率，它們與溫度有關(guān)；vp是伴隨蒸發(fā)導(dǎo)致的固體目標表面的退化速率；α（T），R（T）分別是激光吸收系數(shù)和表面反射系數(shù)，它們也與溫度有關(guān)；I0（t）是入射激光強度。式（1）實際給出的是單位空間中能量的時間變化率。

按照隨機算法，將時間劃分成許多小的時間步長，將空間劃分成許多子空間，每一子空間具有一定數(shù)量的“虛粒子”，其數(shù)目取決于該子空間的溫度，或者是該子空間中的能量，一個“虛粒子”就相當于一個能量量子，當固體中發(fā)生能量的入射或能量的轉(zhuǎn)移時，就用能量量子的轉(zhuǎn)移來代表。

考慮具有軸對稱分布的入射激光束，以入射激光束的中心軸線為主軸將目標體按照柱坐標細分，i代表徑向的指標，基數(shù)用r表示，j代表厚度方向（軸向）的指標，基數(shù)用d表示，如圖1所示。指標為（i，j）的子空間的體積用V i，j表示，相應(yīng)的虛粒子數(shù)目為qi，j，qi，j的變化將引起該子空間能量的變化。溫度隨能量量子數(shù)qi，j的變化可以用式（2）來表示。

式中：Tinit是初始溫度，Q0是“虛粒子”的能量單位，Ci，j是該子空間的熱容。這樣式（1）的熱輸運過程就轉(zhuǎn)化成了各單元子空間上能量（或溫度）的局部變化問題。

子空間的大小可以根據(jù)實際問題的需要精細劃分，因而每一子空間的溫度將隨時間變化，但任一子空間上的溫度分布是均勻一致的。同時規(guī)定一個子空間上目標體的質(zhì)量密度相同，只有一

種材料組分甚至是空隙，與溫度有關(guān)的特征參數(shù)如熱容、熱傳導(dǎo)率、光吸收系數(shù)等在整個子空間上各向同性、均勻一致。這樣統(tǒng)計方法即使是處理具有復(fù)雜結(jié)構(gòu)和材料組分的系統(tǒng)，也不會帶來更多復(fù)雜的計算。各子空間上相應(yīng)于激光吸收、蒸發(fā)以及熱傳導(dǎo)的事件概率可表示如下（根據(jù)圖1所示的柱坐標）：激光吸收概率

這里Sr和S h分別是子空間的側(cè)面積和底面積，d r和d h分別是子空間的寬度和高度。

按照Possion統(tǒng)計，兩個動力學(xué)事件之間的時間步長依賴于瞬時發(fā)生的所有事件的可能性，可以表示為

其中R0是介于0與1之間的隨機數(shù)。

2 計算結(jié)果

采用上述算法，本文計算了激光輻照下單晶硅的熱效應(yīng)問題。計算中假定入射激光為高斯型連續(xù)激光，光斑半徑為3.54μm，光功率可調(diào)，環(huán)境溫度和介質(zhì)硅初始溫度為300 K，表1為常溫下硅的熱物理參數(shù)，各熱物理參數(shù)隨溫度變化規(guī)律如下：

式中α，γ，β為可調(diào)參數(shù)，T0是參考溫度。

圖2為從介質(zhì)表面離激光光斑中心不同距離處的溫度分布，同一時刻介質(zhì)表面溫度場的徑向分布與激光場能量分布相似，基本上也呈高斯狀分布特征。當介質(zhì)熱導(dǎo)率較大（k＝457 W／（m·K））時，介質(zhì)溫升幅度較小；當介質(zhì)的熱導(dǎo)率較小（k＝157 W／m K）時，溫升反而較大。這是因為隨著介質(zhì)熱導(dǎo)率的增加，介質(zhì)在軸向和徑向方向的熱擴散概率相應(yīng)增大，熱擴散量相對增加，介質(zhì)光斑中心處從入射激光束吸收的能量能更及時地沿徑向和軸向擴散出去，因此溫度也相對較低。

表1 Si的熱物理參數(shù)

圖3給出了不同激光輻照功率下，光斑中心處介質(zhì)軸向上的溫度分布。由圖3可知，介質(zhì)軸向溫升隨著激光輻照功率的增加而增加，介質(zhì)溫度沿著軸向呈指數(shù)衰減趨勢。

圖4、圖5給出了5 W激光功率下不同反射率及比熱時，介質(zhì)表面溫度的徑向分布和表面中心處溫度隨輻照時間的變化。從圖中可以看出，考慮介質(zhì)熱物理性質(zhì)隨溫度的變化和不考慮熱物理性質(zhì)隨溫度變化，能導(dǎo)致結(jié)果的較大差異。而實際情況下物體熱性質(zhì)一般都與溫度有關(guān)，并且溫度越高熱性質(zhì)的溫度相關(guān)性越不容忽視。激光輻照固體材料，尤其是高功率激光輻照固體材料時，能短時間內(nèi)把物質(zhì)加熱到很高溫度，因此通過近似的方法解析求解激光熱效應(yīng)問題，能極大地影響結(jié)果的精度。

圖4 不同參數(shù)下介質(zhì)徑向溫度的時間變化

圖5 不同參數(shù)下介質(zhì)中心溫度的時間變化

3 結(jié) 語

在熱傳導(dǎo)方程的基礎(chǔ)上，對比分析了隨機模擬方法與傳統(tǒng)計算方法的優(yōu)劣，并在柱坐標系下建立起激光輻照固體材料熱效應(yīng)問題隨機模擬的計算模型。由于該模型中各個物理量均可調(diào)，而且也可以考慮熱物理性質(zhì)隨溫度的非線性變化，因此具有普遍的適用性。該方法的突出特點是計算過程簡單，卻能適用于復(fù)雜情形，并且能很容易通過計算機模擬實現(xiàn)最終結(jié)果。

以單晶硅材料為例進行了隨機模擬計算。結(jié)果表明：溫升不僅與材料熱性質(zhì)（熱容、熱導(dǎo)率、密度）密切相關(guān)，而且還與材料的吸收系數(shù)，激光輻照參數(shù)（功率密度、能量分布、光斑大小、作用時間）有關(guān)。

［1］劉曉鳳，李 笑.SiO2保護膜對高反膜激光損傷特性的改善［J］.強激光與粒子束，2010，22（12）：2 861-2 864.

［2］徐立君，張喜和.強激光輻照硅光電探測器的損傷判別研究［J］.激光與紅外，2009，39（12）：1 263-1 266.

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

［4］鄭艷麗，杜太焦，束慶邦，等.不同氣流環(huán)境下激光輻照金屬材料溫升的數(shù)值模擬［J］.強激光與粒子束，2010，22（11）：2 531-2 534.

［5］焦路光，趙國民，陳敏孫.二維情形下組合激光輻照45＃鋼靶溫升效應(yīng)的模擬［J］.紅外與激光工程，2010，39（1）：42-46.

［6］HOULE F A，HINSBERG W D.Stochastic simulation of heat flow with application to laser-solid interactions［J］.Applied Physics A，1998，66（2）：143-151.

TN24

A

1008-1542（2011）12-0011-04

2011-06-20；責(zé)任編輯：陳書欣

龔艷春（1970-），女，江蘇連云港人，副教授，主要從事激光物理與技術(shù)及電磁環(huán)境效應(yīng)方面的研究工作。