基于應力-位移混合有限元法的激光超聲數值模擬

許伯強，劉洪凱，徐桂東，徐晨光，李俊敏

基于應力-位移混合有限元法的激光超聲數值模擬

許伯強1，劉洪凱2，徐桂東1，徐晨光1，李俊敏2

（1.江蘇大學理學院，鎮江212013；2.江蘇大學機械工程學院，鎮江212013）

為了研究激光輻照各向同性半無限大鋁材料內超聲波的激發和傳播特征，采用理想匹配層和應力-位移混合有限元方法建立了半無限大介質中激光激發超聲波的有限元數值模型。模擬研究了各向同性半無限大鋁材料內激發產生的瞬態波場圖和垂直表面位移，并與相同幾何模型下采用有限元方法得到的結果進行了對比分析。結果表明，應力-位移混合有限元方法能夠有效地消除模型截斷邊界處的反射波，精確地模擬出無限大固體材料內激光激發超聲波的產生和傳播特性。數值模擬結果為進一步研究微納米薄膜材料中皮秒或飛秒激光激發超聲波提供了有效的方法。

激光技術；激光超聲；理想匹配層；混合有限元；數值模擬

引 言

激光超聲具有非接觸激發、頻帶寬和多模態等優點［1-3］，研究激光激發的超聲波在材料內的產生和傳播特性，能夠為材料的物理特性和健康監測提供有效的分析手段，因而在無損檢測等領域得到廣泛的應用［4-5］。

自1963年WHITE首次證實了脈沖激光輻照金屬表面可以激發出高頻超聲波以來，國內外學者從實驗和理論上對激光超聲展開了廣泛的研究［6］。在數值模擬研究方面，通常采用有限元（finite element method，FEM）數值計算方法［7-8］來分析，有限元數值方法建立在嚴密的數學理論基礎上，能夠靈活處理復雜的幾何模型，并且能夠得到精確的數值解。然而用有限元數值模擬超聲波，必須考慮足夠大的幾何尺寸以避免從邊界處超聲波反射對計算的影響，這將消耗大量的計算資源和計算時間，降低了計算效率。為了消除或降低有限元計算中截斷邊界處的反射波的影響，國內外學者做了大量的研究，主要采取兩種方法：無反射邊界條件法［9］和理想匹配層（perfectly matched layer，PML）［10］等。無反射邊界條件方法分析準確度較高，但它需要修改有限元求解程序，無法借助于現有的有限元計算平臺執行；1994年，BERENGER［11］在電磁學領域首次提出理想匹配層方法，基于分割場變量方法，將波動方程進行傅里葉變換到頻域，再逆變換到時域計算，有效地避免了將頻域結果變換到時域結果的卷積運算。這一方法很快被應用到超聲波和彈性力學領域［12］。KUCUKCOBAN等人［12］基于PML方法，針對結構力學控制方程，提出將PML的控制函數進行復延拓，經傅里葉變換和逆變換最后得到PML條件下位移-應力場的混合有限元時域方程。這種方法能夠有效地消除截斷邊界處聲波的反射，直接得到時域內的波形及瞬態波場圖，且求解精度較高。

作者基于應力-位移混合有限元方法，建立激光激發超聲波的數值模型，模擬研究了各向同性無限大鋁材料中激光超聲的產生和傳播特征，并與有限元方法獲得的結果進行對比分析，取得較為理想的結果。

1 理論分析

考慮均勻的各向同性線彈性材料，線源激光垂直入射到材料表面，計算域的長和寬分別為15mm，如圖1所示。為了有效消除激光激發的超聲波在邊界處的反射所帶來的影響，在研究區域Ωr周圍添加一定厚度的理想匹配層區域ΩPML來吸收超聲波。圖中，LPML為PML的厚度，Γ1為對稱邊界，ΓPML為PML的固定外邊界，Γf為自由表面。

Fig.1 Schematic diagram of x，y plane after laser irradiating material surface

為了有效地衰減傳播到PML區域的超聲波，引入以下復坐標拉伸函數S～，將原空間坐標映射到復空間坐標系統內。

式中，ω為角頻率，εs是沿坐標s方向的拉伸函數，s′表示原空間坐標s在積分過程中的積分變量，其表達式為：

1.1PM L混合有限元控制方程

基于平面應變理論，熱彈機制下激光激發的超聲波在線彈性材料內的控制方程為：

式中，s取x或y。對于超聲波而言，實部的作用是對原坐標s的尺度縮放，虛部的作用是衰減進入PML區域的超聲波振幅。并且須滿足以下條件：

式中，σ，ε，C，α分別為應力、應變、材料的剛度張量和材料的線性熱膨脹系數，u為位移矢量；ρ和T分別表示材料的密度和溫度；“▽·”，“∶”，上標“T”和變量上標點號“··”分別為梯度算子、內積算符、矩陣轉置和對應變量的時間2階偏導，ΔT為材料受激光輻照引起的溫度變化量。

為了在控制方程中引入PML吸收層，時域內的控制方程組（4）式先經過傅里葉變換到頻域，然后根據拉伸函數（見（1）式）將原坐標系的微分關系變換到頻域方程，即可得到包含PML的頻域控制方程；為避免將頻域計算結果逆變換到時域的瞬態波的卷積運算，以便在時域內對激光激發超聲波進行瞬態分析，將PML形式的頻域控制方程組經過傅里葉逆變換回時域內求解，得到以下應力-位移混合有限元的時域控制方程：式中，D為材料的柔度矩陣（剛度矩陣C的逆）；系數a，b，c和矩陣，，Λe，Λp為推導過程中得到的包含模型縮放因子和衰減因子的量，它們的表達式見參考文獻［13］；“·”和“··”分別表示對應變量的時間1階偏導和2階偏導。上述變換中，引入了應力的時間積分形式S（x，y，t）＝（x，y，τ）dτ和應變的時間積分形式E（x，y，t）＝∫t0ε（x，y，t）dτ作為新的未知變量來降低方程的階數。

邊界條件可表示為：

初始條件為：

1.2網格和時間步長

在混合有限元數值模擬計算中，時間步長的選取和模型網格大小的劃分直接影響計算的精確性、穩定性和計算效率。一般情況下，選擇越小的時間步長得到的計算結果越精確，但是卻會讓求解過程過于漫長，降低了計算效率；但是過長的時間步長又會導致求解的高頻部分出現較大的誤差而降低計算精度。因此，在數值計算中，假設在計算區域內時間離散誤差都不超過δ，那么就要求時間步長必須滿足以下條件［7］：

式中，δ為時間離散誤差，C0常數，U（n－1）是對應某一時刻t給定的解向量，并將其作為下一個計算時間步長的初始條件，通過計算求得t＋Δt時刻的解向量U（n）。

結合時間步長，為了使各個計算的波形達到合適的精度要求，網格單元大小l的選擇應該滿足以下表達式：

式中，l為網格單元尺寸，cij表示材料的剛度系數，是材料內激發的超聲波最小波速。

采用變網格技術，在激光輻照區域設定網格大小為5μm，在遠離激光輻照區域設定網格大小為100μm。初始時間步長設定為0.1ns，根據迭代精度在計算過程中調整時間步長，最大時間步長不超過10ns，給定計算誤差為δ＝10－4。選擇瞬態線性直接求解器進行求解。

1.3材料、激光和PM L參量設置

作者采用傅里葉熱傳導理論來描述激光輻照材料表面所產生的溫度場分布，計算中采用輸出波長為1064nm的高斯型激光，考慮材料的表面的吸收系數和光穿透效應，材料吸收激光能量表達式為：

式中，Rf為材料表面對光的反射率，那么在材料表面對激光能量的吸收系數為（1－Rf）；β為激光能量在材料深度方向的衰減系數（其倒數1／β為材料對激光的光趨膚深度）；I0為激光輻照區域中心的能量密度；g（t）和f（x）分別為激光能量在時域和空間的分布函數：

式中，τ和r分別為脈沖激光的脈寬和輻照線源半寬。輻照激光能量的時域和空間分布示意圖如圖2所示。

Fig.2 Spatial and temporal profile of the heat source due to line-focused laser illumination

忽略材料同外界的熱交換，傅里葉熱傳導方程的邊界條件為n·（－k▽T）＝0，其中k，n分別為熱傳導率和單位法向量。設定材料的初始溫度為T（x，y，0）＝300K。

本文中采用脈寬為10ns、輻照線源半寬為300μm和能量密度為4.0MW／cm2的入射激光光源，采用鋁板作為模型材料，其光反射率Rf＝0.85，光趨膚深度為40μm。計算中材料的熱學和力學參量如表1所示。當材料參量為已知量時，PML的尺度縮放和衰減因子［13］的參量設置只由單一的變量R來控制，設定反射率為R＝10－8，PML的厚度為LPML＝3mm。

Table 1 Thermal and mechanical parameters of aluminum

2 數值模擬結果及討論

為了精確的模擬激光輻照鋁材料表面時超聲波的產生和傳播特性，首先需要得到精確的激光輻照所產生的溫度場分布。數值模擬過程中激光能量輸入可看作體熱源。圖3中給出了激光與鋁材相互作用后激光輻照中心（x＝0μm）不同深度處溫度隨時間的變化曲線。

Fig.3 Temperature distributions at different depths at laser irradiation center

圖3顯示，在激光輻照中心表面（x＝0μm，y＝0μm）和表面下方10μm（x＝0μm，y＝－10μm）處，材料的溫度在大約70ns時達到最大值，之后由于熱傳導作用隨時間緩慢下降。在激光輻照中心表面下方50μm（x＝0μm，y＝－50μm）和100μm（x＝0μm，y＝－50μm）由于超過了光趨膚深度的范圍，主要受熱傳導的作用而呈現出溫度緩慢上升的趨勢。這些激光中心表面下不同深度處的溫度隨時間的變化曲線特征很好地與作者之前研究［8］的有限元數值模擬結果相吻合，充分說明了應力-位移混合有限元模型中溫度場計算的準確性。這種非均勻溫度場等效于熱彈體力源，是激發產生超聲波的直接驅動力。

Fig.4 Wave snapshots of total displacement at in semi-infinite aluminum

圖4中給出了利用應力-位移混合有限元方法在半無限大鋁材料內不同時刻的總位移場。圖4a顯示在半無限大鋁材料內激光激發產生的超聲波包括縱波P、橫波S、頭波H和Rayleigh表面波R。從圖中可以看出P波的傳播速率最快，在2.0μs時已經到達PML與研究區域的分界面處，并且在圖所示的右上角和左下角處已經有比較明顯的衰減；圖4b中的虛線1和虛線2顯示在4.2μs時，部分H波和S波已經傳播到PML區域并且有很明顯的衰減，R波也已經到達PML區域與研究區域的邊界，此時P波已經被PML區域完全吸收；對比圖4c與圖4b可以明顯地看出，在4.6μs時H波和S波的衰減更加明顯，R波也進入到PML區域中，并且有明顯的衰減；隨著時間的推移，5.2μs時幾乎所有的波都被PML區域完全吸收，只殘留部分橫波在材料內繼續傳播。由此可見，這種PML混合有限元方法能夠有效地衰減進入PML區域的超聲波，從而消除或降低超聲波在PML區域外邊界處的反射，以此精確模擬半無限大介質內激光激發超聲波的數值模型，并且利用這種方法能夠得到頻域分析方法所無法得到的全場位移波場圖，實現結果分析的直觀化、可視化。

圖5是分別采用應力-位移混合有限元方法和有限元方法所獲得垂直表面位移波形圖。圖5a～圖5d分別是在距離激光輻照中心2mm，6mm，8mm和12mm處的波形圖，其中實線表示采用應力-位移混合有限元方法（mixed finite element method，MFEM），虛線表示采用有限元方法（FEM）。兩種方法采用同樣的幾何模型：15mm寬、15mm高，所不同的是，應力-位移混合有限元方法包含3mm厚的PML區域。由圖可知，激光激發超聲波在鋁材表面產生的垂直表面位移波形包含3種模態：掠面縱波Sp，表面橫波Ss和單極性Rayleigh波R。由圖可見掠面縱波的傳播速率最快，最先到達探測點，其次是表面橫波，傳播速率最慢的是Rayleigh波，由于表面橫波和Rayleigh波的波速相差不大，接收到波形顯示Rayleigh波和表面橫波基本疊加在一起（如圖5a所示）；掠面縱波的振幅比較小，并且隨著傳播距離的增加有明顯的衰減，Rayleigh波的振幅比較大，且其幅度幾乎不隨傳播距離的增大而衰減；根據Sp波和R波的波峰到達探測點2mm，6mm，8mm和12mm時的時間可以計算出它們的波速分別為5797.1m／s和2857.1m／s；這與它們的理論計算得出的相速率非常接近。由于Sp波和R波的相速度相差很大，隨著傳播距離的增大，兩個波的波形彼此逐漸分離。然而由圖5a～圖5d中虛線部分顯示的波形圖可見，波形已經受反射波的影響而出現形態變化，距離邊界越近的探測點越早被反射波影響，并且影響的程度越大；由于Rayleigh波的振幅比較大，在波形中占主導地位，所以明顯地看出反射的Rayleigh波Rr是反射波的主要成分。從圖5a～圖5d中實線部分的波形與利用大尺寸幾何模型情況下的FEM［14］計算結果取得較好的一致，且來自邊界處的反射波形明顯消除。并且應力-位移混合有限元方法的模擬結果與1996年DOYLE和SCALA［15］在實驗上得到的Rayleigh波波形吻合得非常好。

Fig.5 Surface normal displacement waveforms at different distances of 2mm，6mm，8mm and 12mm in aluminum by using finite element method and mixed finite element method，respectively

3 結 論

基于應力-位移混合有限元方法，建立了各向同性半無限大線彈性材料中激光激發超聲波的數值模型，得到了材料表面的垂直表面位移波形及全場波場圖。計算結果說明，激光在各向同性半無限大材料內同時激發出縱波、橫波、頭波和Rayleigh波。與有限元方法相比，這種應力-位移混合有限元方法有效地消除邊界處的反射波，能夠有效地降低模型的尺寸和節省計算資源，并能精確地模擬激光超聲波在材料內的產生和傳播特征。數值模擬結果為進一步研究微納米薄膜材料中皮秒或飛秒激光激發超聲波提供了有效的方法。

［1］ ROSSIGNOL C，PERRIN B，LABORDE S，et al.Nondestructive evaluation of micrometric diamond films with an interferometric picosecond ultrasonics technique［J］.Journal of Applied Physics，2004，95（8）：4157-4162.

［2］ KUNDU T.Ultrasonic nondestructive evaluation：engineering and baiological material characterization［M］.Boca Raton，Florida，USA：CRC Press，2004：435-494.

［3］ DAIY，XU B Q，LOU Y，et al.Finite element modeling of the interaction of laser-generated ultrasound with a surface-breaking notch in an elastic plate［J］.Optics＆Laser Technology，2010，42（8）：693-697.

［4］ SOHN Y，KRISHNASWAMY S.Scaning laser line source technique using monopolar Rayleigh waves［J］.American Institute of Physics，2004，700（1）：278-285.

［5］ LIU JS，XU ZH，GU G Q.Numerical study on improvement of signal-to-noise ratio of surface acoustic waves based on laser array［J］.Laser Technology，2011，35（3）：403-406（in Chinese）.

［6］ SUN H X，XU B Q.Numerical analysis of laser-generated lamb wave by finite element method in time and frequency domain［J］.Chinese Journal of Lasers，2010，37（2）：537-542（in chinese）.

［7］ JOHNSON C.Numerical solution of partial differential equations by the finite element methods［M］.New York，USA：Cambridge University Press，1987：14-48.

［8］ XU B Q，SHEN Z H，WANG JJ，et al.Thermoelastic finite element modeling of laser generation ultrasound［J］.Journal of Applied Physics，2006，99（3）：033508.

［9］ GIVOLID，KELLER JB.Non-reflecting boundary conditions for elastic waves［J］.Wave Motion，1990，12（3）：261-279.

［10］ BASU U，CHOPRA A K.Perfectly matched layers for transient elastodynamics of unbounded domains［J］.International Journal for Numerical Methods in Engineering，2004，59（8）：1039-1074.

［11］ BERENGER JP.A perfectly matched layer for the absorption of electromagnetic waves［J］.Journal of Computational Physics，1994，114（2）：185-200.

［12］ FESTAG，NIELSEN S.PML absorbing boundaries［J］.Bulletin of the Seismological Society of America Definition，2003，93（2）：891-903.

［13］ KUCUKCOBAN S，KALLIVOKASL F.Mixed perfectly matched layers for direct transient analysis in 2-D elastic heterogeneous media［J］.Computer Methods in Applied Mechanics and Engineering，2011，200（1／4）：57-76.

［14］ XU B Q，SHEN Z H，NIX W，et al.Numerical simulation of laser-generated ultrasound by the finite element method［J］.Journal of Applied Physics，2004，95（4）：2116-2122.

［15］ DOYLE PA，SCALA CM.Near-field ultrasonic Rayleigh waves from a laser line source［J］.Ultrasonics，1996，34（1）：1-8.

Mixed stress-displacement finite element method for laser-generated ultrasound

XUBaiqiang1，LIUHongkai2，XUGuidong1，XUChenguang1，LIJunmin2
（1.Faculty of Science，Jiangsu University，Zhenjiang212013，China；2.School of Mechanical Engineering，Jiangsu University，Zhenjiang 212013，China）

In order to study the generation and propagation of laser-generated ultrasound in isotropic semi-infinite aluminum material，a laser-generated ultrasound in an arbitrary elastic semi-infinite medium model was established by using mixed stress-displacement finite element method and perfectly matched layer（PML）.The transient wave snapshots and surface normal displacement waveforms in semi-infinite aluminum materials were obtained.The surface normal displacement waveforms were compared with the same geometrical finite element model.The results show that the mixed stressdisplacement finite element method can effectively eliminate reflection waves from truncated boundary，and simulate the generation and propagation of ultrasound in semi-infinite solid material accurately.The simulation results provide an effective method for research of the laser-generated ultrasound waves in micro-nanostructure by picosecond or femtosecond laser irradiation.

laser technique；laser-generated ultrasonic；perfectly matched layer；mixed finite element method； numerical simulation

O426.1；O411.3

A

10.7510／jgjs.issn.1001-3806.2014.02.018

1001-3806（2014）02-0230-06

國家自然科學基金資助項目（11172114）；江蘇省高校自然科學研究資助項目（10KJA140006）；江蘇省六大人才高峰資助項目（2012-ZBZZ-027）

許伯強（1963-），男，博士，博士生導師，現主要從事超聲無損檢測的研究。

E-mail：bqxu＠ujs.edu.cn

2013-04-11；

2013-04-23