基于應(yīng)力-位移混合有限元法的激光超聲數(shù)值模擬

許伯強(qiáng)，劉洪凱，徐桂東，徐晨光，李俊敏

基于應(yīng)力-位移混合有限元法的激光超聲數(shù)值模擬

許伯強(qiáng)1，劉洪凱2，徐桂東1，徐晨光1，李俊敏2

（1.江蘇大學(xué)理學(xué)院，鎮(zhèn)江212013；2.江蘇大學(xué)機(jī)械工程學(xué)院，鎮(zhèn)江212013）

為了研究激光輻照各向同性半無限大鋁材料內(nèi)超聲波的激發(fā)和傳播特征，采用理想匹配層和應(yīng)力-位移混合有限元方法建立了半無限大介質(zhì)中激光激發(fā)超聲波的有限元數(shù)值模型。模擬研究了各向同性半無限大鋁材料內(nèi)激發(fā)產(chǎn)生的瞬態(tài)波場圖和垂直表面位移，并與相同幾何模型下采用有限元方法得到的結(jié)果進(jìn)行了對比分析。結(jié)果表明，應(yīng)力-位移混合有限元方法能夠有效地消除模型截?cái)噙吔缣幍姆瓷洳ǎ_地模擬出無限大固體材料內(nèi)激光激發(fā)超聲波的產(chǎn)生和傳播特性。數(shù)值模擬結(jié)果為進(jìn)一步研究微納米薄膜材料中皮秒或飛秒激光激發(fā)超聲波提供了有效的方法。

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

引 言

激光超聲具有非接觸激發(fā)、頻帶寬和多模態(tài)等優(yōu)點(diǎn)［1-3］，研究激光激發(fā)的超聲波在材料內(nèi)的產(chǎn)生和傳播特性，能夠?yàn)椴牧系奈锢硖匦院徒】当O(jiān)測提供有效的分析手段，因而在無損檢測等領(lǐng)域得到廣泛的應(yīng)用［4-5］。

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

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

1 理論分析

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

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

為了有效地衰減傳播到PML區(qū)域的超聲波，引入以下復(fù)坐標(biāo)拉伸函數(shù)S～，將原空間坐標(biāo)映射到復(fù)空間坐標(biāo)系統(tǒng)內(nèi)。

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

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

基于平面應(yīng)變理論，熱彈機(jī)制下激光激發(fā)的超聲波在線彈性材料內(nèi)的控制方程為：

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

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

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

邊界條件可表示為：

初始條件為：

1.2網(wǎng)格和時(shí)間步長

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

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

結(jié)合時(shí)間步長，為了使各個計(jì)算的波形達(dá)到合適的精度要求，網(wǎng)格單元大小l的選擇應(yīng)該滿足以下表達(dá)式：

式中，l為網(wǎng)格單元尺寸，cij表示材料的剛度系數(shù)，是材料內(nèi)激發(fā)的超聲波最小波速。

采用變網(wǎng)格技術(shù)，在激光輻照區(qū)域設(shè)定網(wǎng)格大小為5μm，在遠(yuǎn)離激光輻照區(qū)域設(shè)定網(wǎng)格大小為100μm。初始時(shí)間步長設(shè)定為0.1ns，根據(jù)迭代精度在計(jì)算過程中調(diào)整時(shí)間步長，最大時(shí)間步長不超過10ns，給定計(jì)算誤差為δ＝10－4。選擇瞬態(tài)線性直接求解器進(jìn)行求解。

1.3材料、激光和PM L參量設(shè)置

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

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

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

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

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

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

Table 1 Thermal and mechanical parameters of aluminum

2 數(shù)值模擬結(jié)果及討論

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

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時(shí)達(dá)到最大值，之后由于熱傳導(dǎo)作用隨時(shí)間緩慢下降。在激光輻照中心表面下方50μm（x＝0μm，y＝－50μm）和100μm（x＝0μm，y＝－50μm）由于超過了光趨膚深度的范圍，主要受熱傳導(dǎo)的作用而呈現(xiàn)出溫度緩慢上升的趨勢。這些激光中心表面下不同深度處的溫度隨時(shí)間的變化曲線特征很好地與作者之前研究［8］的有限元數(shù)值模擬結(jié)果相吻合，充分說明了應(yīng)力-位移混合有限元模型中溫度場計(jì)算的準(zhǔn)確性。這種非均勻溫度場等效于熱彈體力源，是激發(fā)產(chǎn)生超聲波的直接驅(qū)動力。

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

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

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

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

［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，F(xiàn)lorida，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

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

許伯強(qiáng)（1963-），男，博士，博士生導(dǎo)師，現(xiàn)主要從事超聲無損檢測的研究。

E-mail：bqxu＠ujs.edu.cn

2013-04-11；

2013-04-23