王首彬,莘嘉慶,周 遠

(1. 天津城建大學控制與機械工程學院,天津 300384;2. 哈爾濱工業(yè)大學儀器科學與工程學院,黑龍江 哈爾濱 150001)

1 引言

熱傳導反問題一直以來是一個跨學科難題,其計算推導的病態(tài)性、傳熱系統(tǒng)的阻尼性和延遲性等因素導致熱傳導反問題求解難度較大[1]。在工程實踐中,邊界條件、材料特性和未知熱源等難以獲得,需要通過溫度測量和反問題算法來獲取這些未知條件和參數(shù)的具體信息。

近年來,學者們針對常規(guī)最小二乘類優(yōu)化方法的改進一直在進行,得到了一些較好的結果。文獻[1]提出了帶復雜變量的徑向積分邊界元法和LM算法相結合的方法,來辨識傳熱學反問題中的導熱系數(shù)。Agarwal等[2]開發(fā)了估計模鹽界面熱通量的新方法,對比計算機輔助冷卻曲線分析更加準確。文獻[3]中提出了一種通過瞬態(tài)溫度測量在線估算具有復雜幾何形狀的非線性導熱系統(tǒng)的時變表面熱流的新方法。文獻[4]使用帶有Tikhonov正則化的輸出最小二乘法和改進的共軛梯度法估計了熱通量等參數(shù)。

將最大似然函數(shù)法引入反問題的研究,能夠提高不確定系統(tǒng)反演過程的魯棒性和準確性。華盛頓大學Emery等聯(lián)合莫斯科航空學院空間系統(tǒng)工程系[5]一起通過擴展的最大似然原理來處理已知模型參數(shù)中不確定性問題的逆解,并進行最佳實驗設計,所得結果證明了整個方法的可靠和準確。美國康奈爾大學的王敬波等[6]提出了一種貝葉斯推理方法,該方法可以計算獲得邊界熱通量的后驗概率密度函數(shù)(PPDF),且反演結果準確性較好。本課題組針對熱傳導反問題,通過利用分散模糊推理法[7]、預測模型法[8]和順序函數(shù)法[9]等進行研究,取得了較好的結果。

本文將預報誤差原理應用到傳熱系統(tǒng)中,此方法不僅具有貝葉斯推理以及最大似然法的統(tǒng)計性質(zhì),而且無需確定觀測變量的聯(lián)合概率密度函數(shù),即可辨識整個時間域中的邊界熱流,且所求結果具有良好的魯棒性與準確性。

2 傳熱學正問題

如圖1所示,二維非穩(wěn)態(tài)傳熱模型Ω為數(shù)值實驗對象,測點為S1、S2…Sd,排成一列,模型控制方程為(1)到(6)。圖中加粗右邊界為熱流邊界。

圖1 二維非穩(wěn)態(tài)傳熱模型

二維非穩(wěn)態(tài)傳熱系統(tǒng)控制方程以及定解條件如下

(0≤x≤Lx,0≤y≤Ly,t≥0)

(1)

T(x,y,0)=T0,0≤y≤Ly,0≤x≤Lx

(2)

(3)

(4)

(5)

(6)

邊界熱流為q(t), 初始溫度是T0,熱擴散系數(shù)為α=(λ/ρc),λ是模型的導熱系數(shù),利用有限差分法[10,11]來求解傳熱正問題。

3 預報誤差模型的建立及求解

針對傳熱反問題的不適定性和延遲性,引入順序函數(shù)法的原理;預報誤差法一般用于控制系統(tǒng),因此根據(jù)傳熱系統(tǒng)的性質(zhì)引入階躍響應的概念。由重構的溫度場、單位階躍響應系數(shù)和順序函數(shù)法原理共同構成預報誤差公式,并通過求解相應的準則函數(shù)的最小值得到邊界熱流值。

3.1 傳熱系統(tǒng)的單位階躍響應系數(shù)

單位階躍響應函數(shù)定義為h(x,y,τk)=?T(x,y,τk)/?q(τk)其中τk為時間離散點k處的時間。h(x,y,τk)反映熱流條件?q(τk)作用下?T(x,y,τk)的改變程度。式(7)中利用后一時刻τk的h(x,y,τk)與前一時刻τk-1的h(x,y,τk-1)做差所求得的Δh(x,y,τk-τk-1)就是從τk到τk-1時刻的單位階躍響應系數(shù)。以此可類推得到不同時間段、不同測點處的單位階躍響應系數(shù),如式(7)到式(9)所示,其中,r為未來時間步數(shù),k為離散時間點(k=1,…,H),H為離散時間點數(shù),x,y分別為橫縱坐標位置。

Δh(x,y,τk-τk-1)=h(x,y,τk)-h(x,y,τk-1)

(7)

Δh(x,y,τk+1-τk-1)=h(x,y,τk+1)-h(x,y,τk)

(8)

?

Δh(x,y,τk+r-1-τk-1)=h(x,y,τk+r-1)-h(x,y,τk-1)

(9)

以上設立的傳熱系統(tǒng)模型中,邊界熱流值是未知的,溫度場受邊界熱流的影響而不斷變化,因此可以將傳熱系統(tǒng)視為一個由邊界熱流為輸入,溫度場為輸出的系統(tǒng)。據(jù)此提出單位階躍響應系數(shù):熱流值增加一個單位,溫度場的溫度相應改變一個數(shù)值,即Δh(x,y,τk-τk-1)。簡化公式令Δh(x,y,τi-τk-1)=Δh(d,m)。式中m=i-(k-1),(i=k,k+1,…,k+r-1)。

3.2 確立預報誤差模型求解反問題

根據(jù)傳熱系統(tǒng)的特點,應用順序函數(shù)思想中由未來時刻的輸出值來預測當下時刻的輸出值,具體過程如下:

已知正問題模型中前k-1個時間點的熱流值q1,q2,…,qk-1和之后r個時間點內(nèi)測點上的溫度測量值Tk,Tk+1,…,Tk+r-1,同時假定時間τk,τk+1,…,τk+r-1的熱流變化滿足式(10)

qk=qk+1=…=qk+r-1

(10)

根據(jù)傳熱系統(tǒng)的階躍響應系數(shù)、順序函數(shù)原理和預報誤差法原理建立預報誤差模型如下

T(d,τm,qk)=T(d,τk,0)+Δh(d,m)*qk

(11)

式中T(d,τm,qk)是τm時刻d測點的測量溫度,T(d,τk,0)是τk時刻d測點且假設當qk=0所重構的計算溫度。Δh(d,m)是測點d處對應時刻的階躍響應系數(shù)。根據(jù)順序函數(shù)原理,假設qk為全部相同的辨識參數(shù)。

(12)

σ是通過以下公式獲得的溫度測量誤差的標準差

(13)

根據(jù)在反演熱流值計算求得的重構溫度值T(d,τk,0)和階躍響應系數(shù)值Δh(d,m)的條件下,“最佳”的預報溫度T(d,τm,qk)期望為

(d,τm,qk)=E[T(d,τm,qk)|T(d,τk,0),Δh(d,m)]

(14)

可以寫作

E{‖T(d,τm,qk)-(d,τm,qk)‖2|T(d,τk,0),Δh(d,m)}

=min

(15)

顯然這里的“最佳”溫度預測就是模型的最佳輸出。通過極小化預報誤差準則式子(16),就可以獲得“最佳”的熱通量值:

(16)

其中

(17)

(18)

(19)

(20)

公式中

(21)

為了使得到的二階導數(shù)正定,同時減少計算量,實際中常使用公式(20)作為近似二階導數(shù)。

3.3 反問題求解流程

使用Newton-Raphson 法求邊界熱流的計算步驟如下:

4 數(shù)值實驗及討論

4.1 數(shù)值實驗條件

選擇二維正方形平板為數(shù)值實驗模型,其參數(shù)如下:長度Lx=0.55m;寬度Ly=0.55m;ρc=4×106J/(m3·K);λ=47W/m·℃;傳熱系統(tǒng)的初始溫度是20℃。假定熱流邊界的函數(shù)形式如下式

(22)

4.2 熱流殘差和相對平均誤差

定義熱流殘差如式子(23)所示

(23)

(24)

4.3 測點數(shù)目對反演熱流的影響

圖2和圖3為固定測點位置距離x=Lx邊界xp=0.005m,時間步數(shù)r=4,測量誤差為σ=0.005時,不同測點數(shù)目下,脈沖形式的熱流的反演值與真實值的對比,其中d分別為40和110。表1為不同測點數(shù)目下,熱流殘差和相對平均誤差的對比。

表1 不同測點個數(shù)下,反演熱流的熱流殘差和相對平均誤差

圖2 測點數(shù)目為40時,反演熱流值與真實熱流值的對比

圖3 測點數(shù)目為110時,反演熱流值與真實熱流值的對比

由圖2和圖3以及表1可以得出,測點數(shù)目的增加能夠降低相對平均誤差和熱流殘差。從反演熱流值的振蕩變化可以發(fā)現(xiàn),測點個數(shù)的增長可以抑制振蕩的幅度,這也說明本文所提算法對反演熱流值具有濾波效果。

4.4 測量誤差對反演熱流的影響

圖4和圖5是不同測量誤差下,脈沖形式的熱流的反演值與真實值的對比。參數(shù)如下:測點位置距離x=Lx邊界xp=0.005m,時間步數(shù)r=5,測點個數(shù)d=70,測量誤差分別為σ=0.005和σ=0.01。表2表示不同誤差下,反演熱流的熱流殘差和相對平均誤差的對比。

圖4 在σ=0.005時,反演熱流值與真實熱流值的對比

圖5 在σ=0.01時,反演熱流值與真實熱流值的對比

由圖4、圖5和表3可得,隨著測量誤差的增加,反演的熱流波動逐漸劇烈且波動幅度增加,同時反演熱流的熱流殘差和相對平均誤差不斷增加,但準確度依然較高。

表3 不同未來時間步數(shù)下,反演熱流的熱流殘差和相對平均誤差

4.5 未來時間步數(shù)對反演熱流的影響

在不同未來時間步數(shù)(r=3和r=7)下,固定測點位置距離x=Lx邊界xp=0.005m,固定測點個數(shù)d=60,固定測量誤差為σ=0.01時,給出了脈沖形式的熱流的反演值與真實值的對比(如圖6、圖7所示),以及反演熱流的熱流殘差和相對平均誤差的對比(如表3所示)。

圖6 在r=3時,反演熱流值與真實熱流值的對比

圖7 在r=7時,反演熱流值與真實熱流值的對比

對比圖6和圖7、觀察表3中對應項可知,適當增加未來時間步數(shù)可以有效減少相對平均誤差和熱流殘差。在算法具體實施中控制好未來時間步數(shù),能夠有效降低反演熱流的總體誤差。

4.6 測點位置對反演熱流的影響

測點位置距離x=Lx邊界分別為xp=0.01m和xp=0.015m時, 調(diào)整未來時間步數(shù)可以得到反演熱流值和真實熱流值的對比情況(如圖8、圖9所示),以及不同位置的熱流殘差和相對平均誤差情況(如表4所示)。其中測點個數(shù)為d=50,測量誤差為σ=0.005。

表4 不同測點位置和未來時間步數(shù)下反演熱流的熱流殘差和相對平均誤差

圖8 在xp=0.01m和r=10時,反演熱流值與真實熱流值的對比

圖9 xp=0.015m、r=19時,反演熱流值與真實熱流值的對比

從圖8、圖9以及表4給出的反演結果可以發(fā)現(xiàn),測點位置離熱流表面越遠,所需要的未來時間步數(shù)r越大,即未來時間的溫度信息越多。同時,未來時間步數(shù)的過度增加會導致平均相對誤差的增大。在相同誤差水平下,測點溫度對表面熱流的響應減弱,會導致反演的熱流曲線偏離,使得反演結果誤差增加,但熱流的反演結果依然可觀。

5 結論

本文針對二維非穩(wěn)態(tài)導熱系統(tǒng)邊界熱流反演問題,提出了結合順序函數(shù)法以及階躍響應系數(shù)概念建立預報誤差模型的算法,其準則函數(shù)最大似然原理的應用,提高了反演的精度,增強了反演過程的穩(wěn)定性。數(shù)值算例給出了在不同測點個數(shù)、不同測量誤差、不同未來時間步數(shù)、以及不同測點位置情況下的反演結果,進一步證明了所提算法的魯棒性和準確性。

采用本文方法,根據(jù)溫度測點的位置和測量誤差的水平,選取適當?shù)臏y點個數(shù)和未來時間步數(shù),可以有效地改善邊界熱流的反演效果。當距離熱流邊界較遠或測量誤差較大時,應該增加測點個數(shù)和未來時間步數(shù),增加反演結果的準確性和穩(wěn)定性。