潘天成，呂中榮，汪利

(中山大學(xué)航空航天學(xué)院，廣東 廣州 510006)

在實(shí)際工程中，存在著許多熱現(xiàn)象，及時監(jiān)測結(jié)構(gòu)系統(tǒng)的溫度，識別熱源的位置和強(qiáng)度，對整個熱傳導(dǎo)系統(tǒng)的安全監(jiān)測和溫度控制具有重要意義。例如，通過測量核反應(yīng)堆外部的溫度，從而得到其內(nèi)壁溫度，實(shí)時監(jiān)測和控制整個反應(yīng)的溫度；在電路中容易產(chǎn)生局部的熱源，識別熱源的位置可以用于對電路進(jìn)行檢修；航天飛機(jī)在進(jìn)入大氣層時，經(jīng)常會產(chǎn)生氣動加熱，航天器的局部結(jié)構(gòu)產(chǎn)生熱源，因此根據(jù)局部溫度場來識別熱源的位置和強(qiáng)度就顯得尤為重要。熱源識別屬于熱傳導(dǎo)反問題的研究領(lǐng)域[1]，目的在于識別熱源的位置和強(qiáng)度，及時掌握系統(tǒng)中熱源的屬性。研究表明，熱源識別問題一般是非適定的[2]，導(dǎo)致識別結(jié)果對溫度測量誤差十分敏感。

實(shí)際上，熱源識別的研究主要分為兩類：一是給定熱源的位置或者空間分布，識別熱源的強(qiáng)度；二是同時識別熱源的空間位置和強(qiáng)度。過去幾十年，很多學(xué)者已經(jīng)對熱源識別問題做了大量的研究。為求解此類非適定問題，BECK提出在最小二乘方程中添加正則化項(xiàng)求解，發(fā)展并使用了Tikhonov正則化方法。在一維的熱傳導(dǎo)問題中，BLACK WELL提出了序列估計(jì)方法[3]。HUANG等通過共軛梯度法來研究了一維和二維模型中的一個或兩個點(diǎn)熱源的強(qiáng)度識別問題[4]。但是由于算法迭代的特性，計(jì)算時間比較長。GENG等采用變分迭代算法求得拋物型熱方程無離散化問題的數(shù)值解并確定熱源的強(qiáng)度[5]。NETO等完成了矩形域內(nèi)特定位置的線熱源隨時間變化的強(qiáng)度識別[6]。LI等提出了一種并行蟻群的優(yōu)化算法來識別二維域中的一個點(diǎn)熱源[7]。第一類熱源識別的研究工作比較多，但對于第二類熱源的識別研究卻十分有限。主要難點(diǎn)在于測量得到的溫度與熱源強(qiáng)度和位置之間呈復(fù)雜、難以顯式表達(dá)的非線性關(guān)系，而且在有限元的框架下，測量數(shù)據(jù)量將小于未知參數(shù)數(shù)目，導(dǎo)致問題的嚴(yán)重非適定性。幸運(yùn)的是，稀疏正則化可以很好地解決這一問題[8-9]。

本文將提出一種新的方法來研究熱源識別第二類問題。考慮到點(diǎn)熱源在空間上的分布是稀疏的，引入稀疏正則化[10-12]來對目標(biāo)函數(shù)進(jìn)行約束, 并且采用交替優(yōu)化方法[13]分別對溫度和熱源兩個分離的變量進(jìn)行迭代求解。在尋求正則化參數(shù)的過程中，提出了閾值法。本文所提算法能快速準(zhǔn)確地識別熱源位置和強(qiáng)度，這給實(shí)際工程的應(yīng)用提供了一定的參考。

1 稀疏正則化目標(biāo)函數(shù)建立

1.1 熱傳導(dǎo)問題有限元格式

考慮一個包含邊界?Ω的二維域Ω的熱傳導(dǎo)問題，一般的控制方程為

+Q(x,y,t),

x,y∈Ω,t≥0

(1)

其中，x,y是笛卡爾坐標(biāo)，T為溫度，α為不隨時間變化的熱傳導(dǎo)率，Q是熱源項(xiàng)，本文主要考慮為點(diǎn)熱源的情形。

邊界條件為

=h(Te(x,y,t)-T(x,y,t))

在?qΩ邊界上T(x,y,t)=g(x,y,t),在?TΩ邊界上

?qΩ∪?TΩ=?Ω,

?qΩ∩?TΩ=φ,?TΩ≠φ

(2)

其中，nx、ny分別是二維域邊界的法向方向余弦，h是對流換熱系數(shù)，Te(x,y,t)是周圍的環(huán)境溫度，?TΩ和?qΩ分別為本質(zhì)邊界和對流邊界，g(x,y,t)是在本質(zhì)邊界?TΩ給定的溫度分布。

初始條件為

T(x,y,0)=T0(x,y)

(3)

其中，T0(x,y)為二維域內(nèi)的初始溫度分布。

+Q(x,y,t)=0,

x,y∈Ω,t≥0

(4)

有限元方法是求解熱傳導(dǎo)問題(1)-(4)的有效數(shù)值方法，將二維域Ω劃分為有限個微小的單元dΩ，通過分段線性插值得到的有限元形函數(shù)，可以將域內(nèi)的溫度用各節(jié)點(diǎn)的溫度表示。那么，穩(wěn)態(tài)情況下的有限元格式如下：

(5)

1.2 建立目標(biāo)函數(shù)

(6)

(7)

1.3 稀疏正則化

(8)

其中，λ≥0是一個正則化參數(shù)。至此，基于稀疏正則化的有限元格式穩(wěn)態(tài)熱源識別目標(biāo)函數(shù)已經(jīng)建立。

2 交替優(yōu)化方法

2.1 交替優(yōu)化方法求解目標(biāo)函數(shù)

2.1.1 溫度恢復(fù)步 完整的溫度數(shù)據(jù)從測量得到的溫度中得到恢復(fù)。

(9)

+(F(k-1))TF(k-1)

+(F(k-1))TF(k-1)

其中，公式(10)中的上標(biāo)S和U分別對應(yīng)測量的節(jié)點(diǎn)和未測量的節(jié)點(diǎn)，上標(biāo)US對應(yīng)耦合的部分，I是一個單位矩陣。因此可以得到公式(9)的解，如下

2.1.2 熱源識別步

從恢復(fù)得到的完整溫度數(shù)據(jù)中識別熱源位置和強(qiáng)度，如下

(12)

(13)

其中，(ν)j表示向量ν的第j個分量。顯然，問題(13)是具有顯式解。因此，熱源識別步，即公式(12)的解為

(14)

可以發(fā)現(xiàn)，公式(14)容易得到稀疏正則化目標(biāo)函數(shù)的最優(yōu)解，而且目標(biāo)函數(shù)是解耦的，無需迭代。進(jìn)一步的研究表明，稀疏正則化參數(shù)λ可能會對熱源的強(qiáng)度識別造成一定的偏差，因此應(yīng)避免，本文對公式(14)進(jìn)行了修正，消除λ對識別的熱源強(qiáng)度的影響。

=fm(λ;aj)

(15)

關(guān)于λ的選取將在接下來進(jìn)行討論。

2.2 閾值法尋求正則化參數(shù)

在執(zhí)行每一熱源識別步之前，都必須選取一個合適的正則化參數(shù)。根據(jù)公式(15)，可選取一個閾值集τ。

τ={λjcr:=2|aj|,j=1,2,…n}

(16)

當(dāng)λ≥λjcr時，公式(15)的解為0，即表示在第j個節(jié)點(diǎn)內(nèi)部沒有點(diǎn)熱源存在；λ<λjcr時，結(jié)果表明點(diǎn)熱源在相應(yīng)的節(jié)點(diǎn)內(nèi)部。將閾值集進(jìn)行降序排序

(17)

那么就可以通過閾值法來選取一個正則化參數(shù)，基本算法如圖1所示。其中，lmax和γ是兩個閾值參數(shù)，lmax表示點(diǎn)熱源存在的最大個數(shù)；γ是一個判別比，通常在lmax比較小的情況下，取值較大。

圖1 閾值法確定正則化參數(shù)流程圖Fig.1 The flow chart for determining regularization parameters by threshold setting method

3 數(shù)值算例

3.1 模型建立

本文考慮一個二維穩(wěn)態(tài)的數(shù)值算例。通常，測量的溫度是從有限元模擬并添加一定的噪聲得到。在穩(wěn)態(tài)情況下，噪聲施加下

(18)

其中，Tj是從公式(5)求得的在第j個節(jié)點(diǎn)的溫度數(shù)據(jù)。ε是噪聲水平，一般為0%、2%、5%等等。r是一個標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)。

如圖1所示，薄板為鋁合金材料，長和寬均為0.5 m，密度ρ=2 787 kg/m3，比熱容c=883 J/(kg·K)，熱傳導(dǎo)率α=194 W/(m·℃)， 對流換熱系數(shù)h=10 W/(m2·K)，環(huán)境溫度為0 ℃，一個或多個熱源在薄板的內(nèi)部。每個熱源的面積為0.004 m×0.004 m，相對于板的面積是非常小的，因此可以看作是點(diǎn)熱源。

薄板的有限元單元和節(jié)點(diǎn)信息，如圖2所示。

圖2 二維穩(wěn)態(tài)模型圖Fig.2 Two-dimensional steady-state model

考慮以下幾種工況，來驗(yàn)證所提方法的有效性和準(zhǔn)確性，熱源的強(qiáng)度Q15=1.2×106W/m2，Q20=1.0×106W/m2。

3.2 熱源識別結(jié)果

在識別表1四種工況的過程中，考慮噪聲的隨機(jī)性，進(jìn)行蒙特卡洛(monte carlo)試驗(yàn)，每種工況取100組數(shù)據(jù)進(jìn)行研究，最后取其均值和標(biāo)準(zhǔn)差作為識別結(jié)果。在算法實(shí)現(xiàn)過程中，閾值法選取的參數(shù)分別為lmax=4,γ=5,β=1×106，fb的值則在每一次迭代的過程中選取矩陣ASS的2范數(shù)，即fb=norm(ASS)。

四種工況稀疏正則化與無正則化的識別結(jié)果如圖3-6所示，每種工況第20節(jié)點(diǎn)內(nèi)部的熱源強(qiáng)度識別結(jié)果及相對誤差如表2所示。在無噪聲的情況下，目標(biāo)函數(shù)中有稀疏正則化與無正則化均能準(zhǔn)確識別熱源的位置和強(qiáng)度，無正則化識別的強(qiáng)度誤差也為0。當(dāng)有噪聲水平的時候，稀疏正則化的識別結(jié)果明顯比無正則化的結(jié)果要好，而且能保證解的稀疏性，而無正則化的時候不能保證。同時發(fā)現(xiàn)，噪聲水平的提高，蒙特卡洛實(shí)驗(yàn)得到的熱源識別標(biāo)準(zhǔn)差也會增大，但不影響熱源的定位。工況4表明，僅通過邊界點(diǎn)的測量，本文所提方法也能準(zhǔn)確定位結(jié)構(gòu)內(nèi)部的熱源，且識別的強(qiáng)度誤差不大，誤差為3.62%，而無正則化的算法甚至無法準(zhǔn)確定位。在測量數(shù)據(jù)量比較少的情況下，稀疏正則化的優(yōu)勢明顯地體現(xiàn)出來了。

圖3 工況1有無稀疏正則化識別結(jié)果對比Fig.3 Identified results without and with regularization for scenario 1

表2 第20節(jié)點(diǎn)熱源強(qiáng)度識別結(jié)果Table 2 The identified results of heat source strength at 20th node

圖4 工況2有無稀疏正則化識別結(jié)果對比Fig.4 Identified results without and with regularization for scenario 2

圖5 工況3有無稀疏正則化識別結(jié)果對比Fig.5 Identified results without and with regularization for scenario 3

圖6 工況4有無稀疏正則化識別結(jié)果對比Fig.6 Identified results without and with regularization for scenario 4

4 結(jié) 論

本文考慮到點(diǎn)熱源在空間分布上的稀疏性，將稀疏正則化與交替優(yōu)化方法結(jié)合，提出了一種新的熱源識別方法，很好地克服了熱源識別非適定性的問題。數(shù)值算例的結(jié)果表明：該方法能同時識別熱源的位置和強(qiáng)度，稀疏正則化能加強(qiáng)解的稀疏性，具有較好的抗噪性，并且交替優(yōu)化方法的求解和閾值法尋求正則化參數(shù)使得識別快速準(zhǔn)確，甚至僅通過邊界的測量也可識別內(nèi)部的熱源，這給實(shí)際工程提供了一定的參考意義。