曹 丹，屈惠明

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

1 引言

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

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

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

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

2．1 模型的建立

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

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

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

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

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

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

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

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

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

邊界條件為:

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

2．2 模型的求解

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3．2 共軛梯度法(CGM)

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

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

其中，α為搜索步長(zhǎng);d為共軛搜索方向，有:

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

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

搜索步長(zhǎng)αn可通過(guò)缺陷信息的優(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ù)的收斂條件為:

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

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

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

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

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

4．1 缺陷深度對(duì)缺陷檢測(cè)的影響

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

文中通過(guò)ANSYS仿真可以判定該算法可達(dá)到無(wú)損檢測(cè)的目的。另外，此算法不僅可獲得缺陷的深度，還能檢測(cè)缺陷大小，為紅外成像無(wú)損檢測(cè)提供了很好的方法。目前文中的算法尚未對(duì)缺陷導(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)劉光廷，邱德?。S熱傳導(dǎo)問(wèn)題的間接邊界元法［J］．清華大學(xué)學(xué)報(bào):自然科學(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)李立超，楊錄，張艷花．利用紅外圖像處理方法檢測(cè)材料缺陷［J］．紅外與激光工程，2010，39(2):372 －376．