基于預(yù)報(bào)誤差法的邊界熱流反演方法研究

王首彬,莘嘉慶,周 遠(yuǎn)

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

1 引言

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

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

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

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

2 傳熱學(xué)正問(wèn)題

如圖1所示,二維非穩(wěn)態(tài)傳熱模型Ω為數(shù)值實(shí)驗(yàn)對(duì)象,測(cè)點(diǎn)為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,熱擴(kuò)散系數(shù)為α=(λ/ρc),λ是模型的導(dǎo)熱系數(shù),利用有限差分法[10,11]來(lái)求解傳熱正問(wèn)題。

3 預(yù)報(bào)誤差模型的建立及求解

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

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

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

Δ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)

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

3.2 確立預(yù)報(bào)誤差模型求解反問(wèn)題

根據(jù)傳熱系統(tǒng)的特點(diǎn),應(yīng)用順序函數(shù)思想中由未來(lái)時(shí)刻的輸出值來(lái)預(yù)測(cè)當(dāng)下時(shí)刻的輸出值,具體過(guò)程如下:

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

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

(10)

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

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

(11)

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

(12)

σ是通過(guò)以下公式獲得的溫度測(cè)量誤差的標(biāo)準(zhǔn)差

(13)

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

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

(14)

可以寫(xiě)作

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

=min

(15)

顯然這里的“最佳”溫度預(yù)測(cè)就是模型的最佳輸出。通過(guò)極小化預(yù)報(bào)誤差準(zhǔn)則式子(16),就可以獲得“最佳”的熱通量值:

(16)

其中

(17)

(18)

(19)

(20)

公式中

(21)

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

3.3 反問(wèn)題求解流程

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

4 數(shù)值實(shí)驗(yàn)及討論

4.1 數(shù)值實(shí)驗(yàn)條件

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

(22)

4.2 熱流殘差和相對(duì)平均誤差

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

(23)

(24)

4.3 測(cè)點(diǎn)數(shù)目對(duì)反演熱流的影響

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

表1 不同測(cè)點(diǎn)個(gè)數(shù)下,反演熱流的熱流殘差和相對(duì)平均誤差

圖2 測(cè)點(diǎn)數(shù)目為40時(shí),反演熱流值與真實(shí)熱流值的對(duì)比

圖3 測(cè)點(diǎn)數(shù)目為110時(shí),反演熱流值與真實(shí)熱流值的對(duì)比

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

4.4 測(cè)量誤差對(duì)反演熱流的影響

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

圖4 在σ=0.005時(shí),反演熱流值與真實(shí)熱流值的對(duì)比

圖5 在σ=0.01時(shí),反演熱流值與真實(shí)熱流值的對(duì)比

由圖4、圖5和表3可得,隨著測(cè)量誤差的增加,反演的熱流波動(dòng)逐漸劇烈且波動(dòng)幅度增加,同時(shí)反演熱流的熱流殘差和相對(duì)平均誤差不斷增加,但準(zhǔn)確度依然較高。

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

4.5 未來(lái)時(shí)間步數(shù)對(duì)反演熱流的影響

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

圖6 在r=3時(shí),反演熱流值與真實(shí)熱流值的對(duì)比

圖7 在r=7時(shí),反演熱流值與真實(shí)熱流值的對(duì)比

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

4.6 測(cè)點(diǎn)位置對(duì)反演熱流的影響

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

表4 不同測(cè)點(diǎn)位置和未來(lái)時(shí)間步數(shù)下反演熱流的熱流殘差和相對(duì)平均誤差

圖8 在xp=0.01m和r=10時(shí),反演熱流值與真實(shí)熱流值的對(duì)比

圖9 xp=0.015m、r=19時(shí),反演熱流值與真實(shí)熱流值的對(duì)比

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

5 結(jié)論

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

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