曹 丹，屈惠明

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

1 引言

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

文獻［1］～［3］通過脈沖加熱的方法來獲得物體內(nèi)部缺陷信息，該技術(shù)需要高峰值的能量熱源，同時它對于被檢試樣表面的發(fā)射率和非均勻加熱反饋敏感，而且使用該方法探測到的缺陷深度及表面溫差都較小。利用鎖相加熱法［4］進行紅外熱成像無損檢測的方法雖然可以克服脈沖加熱法的某些缺點，但是鎖相加熱法局限性在于:特定頻率的熱源只能探測到特定深度的缺陷，因此，為了探測材料內(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，進而得到內(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)達到最小值時的缺陷信息的參數(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ù)已知缺陷的材料試樣進行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，誤差也較大，此時的缺陷深度可認為已經(jīng)分辨不出了。

在程序中，缺陷的信息使用參量x1，x2，x3表示，x1和x2分別表示缺陷方向兩個邊界的坐標(biāo)，x3表示缺陷y方向上的坐標(biāo)，使用共軛梯度算法對于不同深度缺陷進行反演，從程序運行的結(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é)果已取得一定的進步。目前對于缺陷的檢測也多局限于距離表面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ǎo)反問題的目的，上述結(jié)果驗證了算法的可行性。從第4節(jié)的討論中得出:對于所研究的缺陷類型，材料檢測表面的溫差越大，材料缺陷的檢測越準(zhǔn)確。對于溫差較小情況，相較之前的結(jié)果，材料中缺陷檢測的準(zhǔn)確度已經(jīng)有了顯著提高。

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

［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)劉光廷，邱德?。S熱傳導(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．