基于有限元的材料紅外無損檢測研究

曹 丹，屈惠明

(南京理工大學(xué)電子工程與光電技術(shù)系，江蘇南京210094)

1 引言

紅外熱成像是目前眾多用于預(yù)防和預(yù)測性維護的無損檢測方法之一。由于其具有簡單、快捷、無接觸、檢測面積大及無有害輻射等諸多優(yōu)勢，在過去的幾十年，紅外熱成像在無損檢測領(lǐng)域獲得了廣泛的應(yīng)用。紅外無損檢測借助于熱傳導(dǎo)模型來映射物體的內(nèi)部特征，由于大部分物體都具有導(dǎo)熱性，所以紅外無損檢測廣泛應(yīng)用于金屬、半導(dǎo)體等多種材料的缺陷探測。

文獻(xiàn)［1］～［3］通過脈沖加熱的方法來獲得物體內(nèi)部缺陷信息，該技術(shù)需要高峰值的能量熱源，同時它對于被檢試樣表面的發(fā)射率和非均勻加熱反饋敏感，而且使用該方法探測到的缺陷深度及表面溫差都較小。利用鎖相加熱法［4］進(jìn)行紅外熱成像無損檢測的方法雖然可以克服脈沖加熱法的某些缺點，但是鎖相加熱法局限性在于:特定頻率的熱源只能探測到特定深度的缺陷，因此，為了探測材料內(nèi)部不同深度的缺陷，需要反復(fù)使用不同頻率熱波鎖相加熱。有限差分法簡單、方便，但僅限于求解形狀規(guī)則和材質(zhì)均勻的物體，它的適應(yīng)性較差，求解導(dǎo)熱方程得到的近似解不準(zhǔn)確，不適用于非均勻網(wǎng)格的問題的求解。邊界元法雖然對問題求解的精度高，但是溫度的測量位置位于物體內(nèi)部，且易受缺陷信息的影響［5］。在過去的二十多年，通過求解一維熱傳導(dǎo)反問題來檢測內(nèi)部缺陷的研究較多，然而對于多維熱傳導(dǎo)的反問題的研究才剛剛起步，而且對于多維熱傳導(dǎo)問題的研究也多是瞬態(tài)問題，介紹穩(wěn)態(tài)多維熱傳導(dǎo)問題文章則較少［6］。

本文主要研究材料內(nèi)部缺陷檢測問題，基于紅外成像無損檢測的原理，通過穩(wěn)態(tài)二維傳熱模型來反演內(nèi)部缺陷的信息。為了克服有限差分和邊界元的不足，提出了有限單元法來求解二維導(dǎo)熱方程，利用ANSYS仿真的結(jié)果作為測點的理論值，借助于共軛梯度法求解二維熱傳導(dǎo)反問題。

2 模型的建立及熱傳導(dǎo)問題的求解

2．1 模型的建立

為了獲取材料內(nèi)部缺陷信息，我們選用材料試樣可用圖1中的模型表示，其中陰影部分表示該材料模型的缺陷部分。

圖1 試樣缺陷檢測模型圖Fig．1 the sample defect testingmodel diagram

對于圖1中的材料模型圖，我們做了以下假定:

(1)材料的內(nèi)部熱流是平穩(wěn)的;

(2)材料x=0，x=l邊分別絕熱;

(3)材料y=0邊施加強度均勻的熱流;

(4)材料y=m邊和空氣形成自熱對流;

(5)紅外熱像儀位于邊y=m上。

綜合以上因素，我們建立材料內(nèi)部的傳熱方程為:

邊界條件為:

其中，λx為材料x方向熱傳導(dǎo)系數(shù);λy為材料y方向熱傳導(dǎo);Ts材料表面的溫度;Tf為空氣的溫度;h為材料和空氣對流換熱系數(shù);qs材料y=0邊所施加的熱流。

2．2 模型的求解

利用有限元分析穩(wěn)態(tài)熱傳導(dǎo)問題時，不需要考慮物體的初始溫度分布對最后穩(wěn)定溫度場的影響，因此不必考慮溫度場的初始條件，而只需考慮換熱邊界條件。穩(wěn)態(tài)溫度場的計算可轉(zhuǎn)化為求解偏微分方程的邊值問題。溫度場是標(biāo)量場，將物體離散成有限單元后，根據(jù)每個單元節(jié)點溫度和形函數(shù)的關(guān)系得到二維模型溫度分布及缺陷的信息。

文中采用Galerkin法來建立單元穩(wěn)態(tài)溫度場分析的一般有限元列式，根據(jù)材料的模型方程(1)和邊界條件方程(2)～(4)得到單元的加權(quán)積分公式為:

其中，N代表有限元形狀函數(shù)矩陣;Ni為形函數(shù)。

按照有限元的格式可以將方程(5)表示為:

其中，矩陣［K］e為單元的導(dǎo)熱矩陣;{T}e為單元的結(jié)點溫度向量;{P}e稱為單元的溫度載荷向量。

整個物體上的加權(quán)積分方程為各個單元積分方程之和，根據(jù)單元結(jié)點的局部編號和整體編號的關(guān)系，直接求和得到整體剛度矩陣，整體方程組為:

其中，矩陣［K］為導(dǎo)熱矩陣;{T}為結(jié)點溫度向量;{F}稱為溫度載荷向量。

對于材料的二維熱傳導(dǎo)問題，將整個矩形區(qū)域離散為m×n網(wǎng)格，利用有限元法求解模型方程(1)，根據(jù)方程(7)即可得到各個結(jié)點處的溫度Ti，進(jìn)而得到內(nèi)部溫度的分布。

3 材料缺陷信息獲取算法的研究

在全空間熱傳導(dǎo)的反問題求解方法中，為了估算含有多個參量的內(nèi)部缺陷問題，可以通過最小化誤差平方和的函數(shù)的方法來實現(xiàn)［7］。在所有優(yōu)化技術(shù)中，共軛梯度法是一種整體收斂速度快而且被廣泛使用一種方法，它可以用來解決整個空間的傳導(dǎo)反問題。這里使用共軛梯度法來最小化誤差平方和的函數(shù)，求解熱傳導(dǎo)反問題，獲取材料內(nèi)部的缺陷信息。

3．1 內(nèi)部缺陷信息的優(yōu)化目標(biāo)函數(shù)

為了獲得材料內(nèi)部的缺陷的深度和大小信息，需要求解材料無損檢測中熱傳導(dǎo)的反問題，為此，把有限元求解方程(1)可得到的材料表面測量點溫度作為Tj，利用ANSYS仿真后獲得的表面測點溫度作為Yj(j=1，2，…，m)，這兩個溫度值的誤差平方和作為目標(biāo)函數(shù)，其中m為表面溫度測量點的個數(shù)。如果以圖1中的矩形缺陷求解，需要確定材料內(nèi)部參數(shù)x1，x2，…，xq，其中q為表示缺陷信息的參數(shù)的個數(shù)，求解材料內(nèi)部缺陷信息的優(yōu)化目標(biāo)函數(shù)為J(r):

式中，r為待反演材料缺陷信息的參數(shù)向量，即r={x1，x2，…，xq}，Tj(r)為測點 j處的計算溫度，可根據(jù)r的猜測值借助于有限元法求解傳熱方程(1)得到。

獲取材料的缺陷信息等價于使目標(biāo)函數(shù)J(r)達(dá)到最小值時的缺陷信息的參數(shù)向量r。

3．2 共軛梯度法(CGM)

共軛梯度法結(jié)合了共軛性及最速下降法的特點，它利用已知點的梯度構(gòu)造共軛方向作為搜索方向，求出誤差平方和函數(shù)的極小值，既克服了最速下降法的鋸齒現(xiàn)象，又避免了牛頓法計算量大和局部收斂的缺點。

缺陷信息參數(shù)r的迭代式為:

其中，α為搜索步長;d為共軛搜索方向，有:

其中，βn為共軛系數(shù)，有:

其中， J為缺陷信息目標(biāo)函數(shù)的梯度:

搜索步長αn可通過缺陷信息的優(yōu)化目標(biāo)函數(shù)J(rn+1)獲得:

式中， Tj為行向量， Tj=(Tj/x1，Tj/x2，…，Tj/xm)。

3．3 缺陷信息優(yōu)化目標(biāo)函數(shù)的收斂條件

缺陷信息優(yōu)化目標(biāo)函數(shù)的收斂條件為:

式中，ε為無限小的一個正數(shù)，如果把測量誤差考慮在內(nèi)，ε可采用的計算公式［8］為:

式中，k表面溫度測量點的數(shù)目;σ為測溫的標(biāo)準(zhǔn)差。

4 仿真結(jié)果與討論

文中通過ANSYS仿真和有限元求解來檢驗上述缺陷檢測算法的有效性。為了利用ANSYS仿真得到材料表面檢測點的溫度 Yj(j=1，2，…，m)，這里首先根據(jù)已知缺陷的材料試樣進(jìn)行ANSYS仿真。使用ANSYS仿真時所選用的參數(shù)如下所示:

材料的規(guī)格:長20 cm，寬10 cm，材料的導(dǎo)熱系數(shù)為 0．077W/(m2·K)，密度為 232 kg/cm3，比熱容為0．88 kJ/(kg·K)，材料的初始溫度及周圍環(huán)境的溫度均為20℃，缺陷部分的導(dǎo)熱系數(shù)為0．027W/(m2·K)，缺陷大小2 cm ×2 cm。

4．1 缺陷深度對缺陷檢測的影響

缺陷和材料表面之間的距離用d表示，當(dāng)d為1 cm，2 cm，3 cm時，利用ANSYS仿真得到試樣內(nèi)部溫度分布云圖如圖2所示。

圖2 不同缺陷深度對應(yīng)的溫度分布云圖Fig．2 temperature distribution contourswith different defect depth

為了對比材料缺陷深度不同表面溫度分布情況，用圖3表示缺陷深度不同表面溫度分布曲線。

圖3 不同缺陷深度對應(yīng)的表面溫度分布Fig．3 the surface temperature distribution curve with different defect depth

由圖2和圖3結(jié)果可以看出，當(dāng)缺陷距離表面1 cm時，表面溫度相差約為1℃，當(dāng)缺陷距離表面2 cm時，表面溫度相差約為0．6℃，當(dāng)缺陷距離表面3 cm處，表面溫度相差約為0．4℃，雖然溫差都大于熱像儀的溫度分辨率，可是由于使用有限元求解方程(1)本身產(chǎn)生的誤差，在缺陷距離為3 cm，誤差也較大，此時的缺陷深度可認(rèn)為已經(jīng)分辨不出了。

在程序中，缺陷的信息使用參量x1，x2，x3表示，x1和x2分別表示缺陷方向兩個邊界的坐標(biāo)，x3表示缺陷y方向上的坐標(biāo)，使用共軛梯度算法對于不同深度缺陷進(jìn)行反演，從程序運行的結(jié)果可以看出，初始值的選取對于最終的結(jié)果有一定的影響，初始值和真實越接近，最終的結(jié)果也越準(zhǔn)確，因此根據(jù)有限元仿真的結(jié)果，需要首先估算x1，x2的值，這樣缺陷的深度值x3和最終結(jié)果也越接近，這是由于有限元計算中網(wǎng)格劃分對于熱傳導(dǎo)求解近似中的誤差所產(chǎn)生的。

缺陷深度不同情況下對應(yīng)的結(jié)果如表1所示。

表1 不同缺陷深度下缺陷檢測的結(jié)果Tab．1 defect testing results under different defect depth

從表1中可以看出，缺陷距離表面的深度越小，最終運行的結(jié)果和真實值的誤差越小，對于深度為1 cm的缺陷，程序迭代的準(zhǔn)確性已相對較高，但是對于2 cm的缺陷還是可以粗略地計算，雖然準(zhǔn)確性不高，但是相較之前的結(jié)果已取得一定的進(jìn)步。目前對于缺陷的檢測也多局限于距離表面1 cm處的缺陷。對于缺陷深度不同的材料，缺陷距離表面深度越小，材料表面的溫差越大，結(jié)果越準(zhǔn)確。

4．2 材料的熱導(dǎo)率［9］對缺陷檢測的影響

對于不同熱導(dǎo)率的材料，當(dāng)缺陷距離表面的深度d都為1 cm時，并且在其他條件相同的情況下，根據(jù)ANSYS仿真結(jié)果整理得到不同熱導(dǎo)率材料對應(yīng)的表面溫度分布曲線如圖4所示。

圖4 不同熱導(dǎo)率材料對應(yīng)的表面溫度分布Fig．4 surface temperature distribution curve with different thermal conductivity ofmaterial

由圖4的分布曲線可以看出，隨著材料導(dǎo)熱系數(shù)的增大，表面溫差呈現(xiàn)為先增大后減小的趨勢。當(dāng)材料的熱導(dǎo)率為20 W/(m2·K)的情況下，材料表面的溫差小于熱像儀的分辨率0．2℃，此時已經(jīng)無法分辨材料是否含有缺陷。

根據(jù)ANSYS仿真得到的各熱導(dǎo)率材料的表面溫度，選取其表面測點溫度作為Yj(j=1，2，…，m)，檢驗共軛梯度算法的可行性。在這里僅列出了缺陷的深度的信息下x3，計算結(jié)果如表2所示。

表2 不同熱導(dǎo)率材料的缺陷檢測結(jié)果Tab．2 defect testing results under different thermal conductivitymaterial

表2的結(jié)果表明，缺陷深度的準(zhǔn)確性和材料的熱導(dǎo)率沒有直接關(guān)系，而直接取決于表面溫差的大小。表面溫差越大，缺陷深度的相對誤差越小，材料缺陷檢測的準(zhǔn)確度也越高。

綜合第4．1節(jié)和第4．2節(jié)的結(jié)果可以看出，對于所研究的缺陷類型，材料內(nèi)部缺陷深度和材料的熱導(dǎo)率是影響表面溫差的重要因素，缺陷距離表面的深度越小，材料表面的溫差越大;材料表面溫差隨熱導(dǎo)率增加呈現(xiàn)為先增大后減小的趨勢。最終都體現(xiàn)為材料檢測表面的溫差越大，材料缺陷的檢測越準(zhǔn)確。共軛梯度法的結(jié)論和理論符合比較好，驗證了算法的可行性。

5 結(jié)論

本文通過對含有缺陷的材料試樣建立了二維導(dǎo)熱模型。為了獲取材料內(nèi)部的缺陷信息，我們把它轉(zhuǎn)換成求解二維導(dǎo)熱模型熱傳導(dǎo)的反問題，首先，根據(jù)有限元求解方程(1)獲取材料表面溫度Tj，ANSYS仿真獲取材料表面的溫度 Yj(j=1，2，…，m)，然后確定材料表面溫度的計算值和仿真結(jié)果的誤差平方和的函數(shù)，最后借助于共軛梯度法來優(yōu)化誤差平方和的函數(shù)，從而達(dá)到求解熱傳導(dǎo)反問題的目的，上述結(jié)果驗證了算法的可行性。從第4節(jié)的討論中得出:對于所研究的缺陷類型，材料檢測表面的溫差越大，材料缺陷的檢測越準(zhǔn)確。對于溫差較小情況，相較之前的結(jié)果，材料中缺陷檢測的準(zhǔn)確度已經(jīng)有了顯著提高。

文中通過ANSYS仿真可以判定該算法可達(dá)到無損檢測的目的。另外，此算法不僅可獲得缺陷的深度，還能檢測缺陷大小，為紅外成像無損檢測提供了很好的方法。目前文中的算法尚未對缺陷導(dǎo)熱未知的情況進(jìn)行討論。

［1］ N P Avdelidis，D P Almond．Through skin sensing assessment of aircraftstructures using pulsed thermography［J］．NDT ＆ E International，2004，37(5):353 －359．

［2］ C Ibarra-Castanedo，N P Avdelidis，X Maldague．Qualitative and quantitative assessment of steel plates using pulsed phase thermography［J］．Materials Evaluation，2005，63(11):1128 －1133．

［3］ Ravibabu Mulaveesala，Sanjay Awasthi，Suneet Tuli．Infrared non-destructive characterization of boiler tube［J］．Sensor Letters，2008，6(2):312 －318．

［4］ G Busse，P Eyerer．Thermal wave remote and nondestructive inspection of polymers［J］．Applied Physics Letters，1983，43(4):355 －357．

［5］ Liu Guangting，Qiu Delong．Indirect boundary element method for three-dimensional heat-conductive engineering problems［J］．Journalof Tsinghua University:Sci＆ Tech，1996，36(1):8 －12．(in Chinese)劉光廷，邱德隆．三維熱傳導(dǎo)問題的間接邊界元法［J］．清華大學(xué)學(xué)報:自然科學(xué)版，1996，36(1):8－12．

［6］ Edward Hensel，Richard Hills．Steady-state two-dimensional Inverse heat conduction［J］．Numerical Heat Transfer，1989，15(2):227 －240．

［7］ A Pourshaghaghy，F(xiàn) Kowsary，A Behbahaninia．Comparison of four different versions of the variablemetricmethod for solving inverse heat conduction problems［J］．Heat Mass Transfer，2007，43(3):285 －294．

［8］ Chenghung Huang，Chengchia Chiang，Hsimei Chen．Shape identification problem in estimating geometry of multiple cavities［J］．Journal of Thermophysics and Heat Transfer，1998，12(2):270 －277．

［9］ Li Lichao，Yang Lu，Zhang Yanhua．Defect testofmaterial using infrared image processingmethods［J］．Infrared and Laser Engineering，2010，39(2):372 －376．(in Chinese)李立超，楊錄，張艷花．利用紅外圖像處理方法檢測材料缺陷［J］．紅外與激光工程，2010，39(2):372 －376．