錢煒祺,周 宇,*,邵元培

(1.空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,綿陽 621000;2.中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所,綿陽 621000)

0 引 言

表面熱流辨識(shí)是通過材料邊界或內(nèi)部的溫度測量值來反演表面熱流,是一類典型的熱傳導(dǎo)逆問題,在航空航天、冶金制造等領(lǐng)域有廣泛的應(yīng)用。表面熱流辨識(shí)是數(shù)學(xué)上的不適定問題,不滿足解的穩(wěn)定性,也就是說即使測點(diǎn)溫度存在較小測量白噪聲也可能導(dǎo)致解的“爆破”;具體來說,在頻域,表面熱流的辨識(shí)誤差在白噪聲的影響下隨著頻率的升高逼近無界。這給辨識(shí)精度估計(jì)帶來了難度,而在工程實(shí)際的測熱系統(tǒng)設(shè)計(jì)時(shí),希望能夠在對(duì)熱環(huán)境有大致預(yù)測的情況下對(duì)其辨識(shí)結(jié)果誤差給出較為準(zhǔn)確的估計(jì),進(jìn)而根據(jù)辨識(shí)結(jié)果精度要求來確定傳感器的測量位置和測量精度。文獻(xiàn)[1]針對(duì)這一問題進(jìn)行了探討,結(jié)果顯示,辨識(shí)結(jié)果誤差不僅與防熱材料物性參數(shù)、測量點(diǎn)位置、測量噪聲形式、測量噪聲大小有關(guān),而且還與表面熱流本身的頻域特性有關(guān)。將傅立葉數(shù)作為相似參數(shù),可以對(duì)辨識(shí)結(jié)果誤差進(jìn)行初步的分析。但是這一結(jié)果給出的主要還是定性分析的結(jié)論,在定量評(píng)估上還存在不足。因此,本文工作是文獻(xiàn)[1]的進(jìn)一步深入探索,首先對(duì)給定單一頻率的熱流辨識(shí)誤差進(jìn)行定量分析,然后結(jié)合頻域分析方法[2-3]對(duì)兩個(gè)和多個(gè)給定頻率組合情況下的誤差規(guī)律進(jìn)行分析,初步建立熱流辨識(shí)誤差的定量分析方法。

1 熱流辨識(shí)方法及辨識(shí)結(jié)果誤差定義

對(duì)典型的一維熱傳導(dǎo)問題,采用文獻(xiàn)[1]中的無量綱方法后,傳熱方程可寫為:

初始條件:t＝0:T＝0

其中,xi(i＝1,P)為溫度測點(diǎn)位置;P 為測點(diǎn)數(shù)目;為溫度測量值;v(t)為測量噪聲。通常的熱傳導(dǎo)正問題計(jì)算是給定熱流Q(t),計(jì)算測點(diǎn)的溫度響應(yīng),而熱流辨識(shí)問題則是根據(jù)觀測方程中的測點(diǎn)溫度信息來確定Q(t)。

表面熱流辨識(shí)方法主要有順序函數(shù)法和共軛梯度法[4-5],本文采用共軛梯度法來對(duì)式(1)進(jìn)行表面熱流辨識(shí),共軛梯度法將辨識(shí)問題轉(zhuǎn)化為求合適的Q(t)使如下目標(biāo)函數(shù)達(dá)極小的優(yōu)化問題:

對(duì)上式做變分后可以得到伴隨變量滿足的伴隨方程,同時(shí)可以得到目標(biāo)函數(shù)對(duì)Q 的梯度,然后代入具有“停止準(zhǔn)則”的共軛梯度法進(jìn)行優(yōu)化計(jì)算;“停止準(zhǔn)則”是一種迭代正則化方法,可以有效抑制辨識(shí)結(jié)果的高頻非物理振蕩。具體算法和有效性驗(yàn)證參見文獻(xiàn)[5,6]。圖1給出了典型算例的辨識(shí)結(jié)果,此時(shí)的表面熱流真值取為頻率f＝2 Hz的正弦函數(shù),即圖1中的“Exact”,測點(diǎn)位于絕熱端x＝1,算例中相應(yīng)的無量綱量k 為1?！癊st.(σ＝0)”表示不考慮測點(diǎn)溫度測量誤差情況下的熱流辨識(shí)結(jié)果,“Est.(σ＝0.05)”表示在測點(diǎn)溫度疊加標(biāo)準(zhǔn)差σ＝0.05Tmean白噪聲情況下(圖2中的“Exp.(σ＝0.05)”)的熱流辨識(shí)結(jié)果,Tmean為測點(diǎn)溫度變化平均值。從圖中可以看到,在不考慮測量噪聲和測量噪聲標(biāo)準(zhǔn)差σ＝0.05的情況下,熱流辨識(shí)結(jié)果和真值都符合;在考慮測量噪聲情況下,根據(jù)熱流辨識(shí)結(jié)果計(jì)算出的測點(diǎn)溫度(圖2中“Fitted(σ＝0.05)”)也與測量值符合,表明辨識(shí)算法是有效的。

圖1 熱流真值與辨識(shí)結(jié)果對(duì)比Fig.1 Comparison of exact and estimated value of heat flux

圖2 測點(diǎn)溫度測量值與辨識(shí)擬合值對(duì)比Fig.2 Comparison of measurement temperature

定義辨識(shí)結(jié)果與真值之間的相對(duì)誤差E 為:

從圖1可以看到,由于傳熱的延遲特性,測點(diǎn)溫度信息對(duì)表面熱流尾段值不敏感,導(dǎo)致熱流尾段值與真值存在一定誤差,因而在式(2)計(jì)算時(shí)只取中間一段來進(jìn)行分析,即t1＝0.25,t2＝1.25。此外,由文獻(xiàn)[1]分析知,在白噪聲的影響下辨識(shí)結(jié)果誤差同樣具有一定的隨機(jī)性。為克服這一影響,采用如下方法:首先給定相同均值、相同標(biāo)準(zhǔn)差的N 組白噪聲;對(duì)疊加這些噪聲后的測量結(jié)果分別辨識(shí)得出相對(duì)誤差值Ei(i＝1,N);然后對(duì)這些相對(duì)誤差值取平均,將該平均值作為對(duì)應(yīng)這一噪聲標(biāo)準(zhǔn)差水平的平均相對(duì)誤差值,記為Em,即

2 給定頻率熱流的辨識(shí)結(jié)果誤差建模

采用上節(jié)方法,首先針對(duì)五組不同頻率(f＝1,2,5,8,10 Hz)正弦形式熱流在不同測量噪聲方差σ＝0.01Tmean、0.02Tmean、0.05Tmean、0.1Tmean情況下的辨識(shí)結(jié)果平均相對(duì)誤差進(jìn)行分析,計(jì)算結(jié)果如表1所示,圖3中給出了頻率10 Hz情況下不同噪聲水平對(duì)應(yīng)的熱流辨識(shí)結(jié)果誤差。

圖3 頻率10 Hz不同噪聲水平對(duì)應(yīng)的熱流辨識(shí)結(jié)果誤差Fig.3 Estimation error for 10 Hz heat flux with reapect to noise standard devition

表1 不同頻率下對(duì)應(yīng)不同噪聲水平的熱流辨識(shí)結(jié)果與給定值平均誤差EmTable 1 Mean error of estimated heat flux with different frequencies for different noise standard devition levels

接下來以表1中的f 和σ/Tmean為輸入,以Em為輸出,采用Kriging響應(yīng)面模型[7]進(jìn)行建模,然后利用該模型對(duì)頻率7 Hz熱流在測量噪聲水平σ＝0.02Tmean、0.04Tmean、0.06Tmean、0.1Tmean情況下的辨識(shí)結(jié)果平均相對(duì)誤差值進(jìn)行預(yù)測,并與實(shí)際計(jì)算結(jié)果進(jìn)行對(duì)比。如圖4所示,圖中“Model prediction”為Kriging模型預(yù)測結(jié)果,“Calculated”為實(shí)際計(jì)算結(jié)果。從圖中可以看到,模型預(yù)測結(jié)果與實(shí)際計(jì)算結(jié)果符合較好,數(shù)值略偏低。

圖4 頻率7 Hz熱流的辨識(shí)結(jié)果誤差預(yù)測結(jié)果對(duì)比Fig.4 Comparison of calculated and predicted estimation error for 7 Hz heat flux

3 組合頻率熱流的辨識(shí)結(jié)果誤差估計(jì)

對(duì)于頻率組合熱流,選取如下幾組熱流進(jìn)行分析:

國家提出的“互聯(lián)網(wǎng)+”及“大數(shù)據(jù)”發(fā)展戰(zhàn)略，使得傳統(tǒng)的包裝及印刷行業(yè)主動(dòng)或被動(dòng)地融入其中，為這個(gè)古老的加工產(chǎn)業(yè)帶來了生機(jī)與希望。其優(yōu)勢(shì)主要體現(xiàn)在遙遠(yuǎn)的距離被拉近，且囊括了整個(gè)加工行業(yè)的方方面面。當(dāng)然，實(shí)現(xiàn)“互聯(lián)網(wǎng)+”及“大數(shù)據(jù)”戰(zhàn)略轉(zhuǎn)型的基礎(chǔ)平臺(tái)之一就是網(wǎng)絡(luò)搜索引擎的應(yīng)用。

采用與上一小節(jié)相同的方法來計(jì)算不同噪聲水平σ＝0.01Tmean、0.02Tmean、0.05Tmean情況下的辨識(shí)結(jié)果平均相對(duì)誤差值。圖5、圖6 給出了σ ＝0.05Tmean情況下熱流b和熱流d的理論值、時(shí)域辨識(shí)值結(jié)果對(duì)比和頻域的功率譜對(duì)比。從頻域分析結(jié)果可以看到,在低頻段,辨識(shí)結(jié)果與真值較為一致,差異主要體現(xiàn)在高頻分量上。由此可知,對(duì)于頻率組合的熱流,低頻分量能在辨識(shí)結(jié)果得到較好地復(fù)現(xiàn),高頻分量是導(dǎo)致辨識(shí)結(jié)果與理論值出現(xiàn)差異的主要原因,因此,辨識(shí)結(jié)果精度可以通過最高頻率熱流分量與測量噪聲之間的對(duì)應(yīng)關(guān)系來進(jìn)行大致估計(jì)。選取兩種頻率組合的情況進(jìn)行分析,以熱流b為例,將其重寫為:

圖5 熱流b理論值、時(shí)域辨識(shí)值結(jié)果和頻域功率譜對(duì)比Fig.5 Comparison of exact and estimated value in time and frequency domain for heat flux b

圖6 熱流d理論值、時(shí)域辨識(shí)值結(jié)果和頻域功率譜對(duì)比Fig.6 Comparison of exact and estimated value in time and frequency domain for heat flux d

并記Q1和Q2的辨識(shí)值為和,則由式(1)的線性可疊加性、函數(shù)正交性和式(2)知:

式中E1、E2表示Q1和Q2的辨識(shí)結(jié)果誤差,由前面分析知,E1≈0,所以上式簡化為:

對(duì)于Q1和Q2,在時(shí)域無法分開,但在頻域可分開,因此,由Parseval(帕塞瓦爾)定理[8]及圖5、圖6可知:

式中Δ為離散變換引起的頻率變化量。則式(6)可另寫為:

其中C 的物理含義為頻域中高頻分量與低頻分量的能量值對(duì)比。

對(duì)于N 種頻率組合(N ≥3)的情況,設(shè)f1＜f2＜…＜fN,由圖6知,與兩種頻率組合時(shí)類似,N 種頻率組合情況下的辨識(shí)結(jié)果與真值差異主要集中在最高頻分量上,此時(shí)的誤差值為:

圖7 熱流c理論值、時(shí)域辨識(shí)值結(jié)果和頻域功率譜對(duì)比Fig.7 Comparison of exact and estimated value in time and frequency domain for heat flux c

表2也給出了對(duì)其余幾組頻率組合情況下采用上述方法的估計(jì)結(jié)果(表中“估算值”),并與采用式(3)和疊加64組同分布白噪聲的熱流辨識(shí)統(tǒng)計(jì)平均結(jié)果(表中“計(jì)算值”)進(jìn)行了比較,可以看到,兩組結(jié)果符合較好。

下面考慮更一般的情況,前面的分析中事先已知了熱流中各分量的頻率、幅值、相位等信息,且熱流均值為0。而在工程實(shí)際應(yīng)用中,這些信息,尤其是熱流的幅值和相位信息不是十分明確,此時(shí)的估算方法需更多地借助頻域分析技術(shù)。考慮圖8(a)的熱流,對(duì)其進(jìn)行頻譜分析,得出功率譜如圖8(b),由功率譜知該熱流的最大頻率為10 Hz;由于熱流均值不為0,頻譜中包含低頻分量;同時(shí)熱流中還包含有頻率為4 Hz的分量。由圖中頻譜曲線積分可知,低頻分量對(duì)應(yīng)的C1≈0.125,4 Hz分量對(duì)應(yīng)的C2≈0.445。接下來還需要確定10 Hz熱流分量的幅值,采用的方法是先濾掉熱流中最大頻率f＝10 Hz的成分,再用原始熱流值減去濾波后的熱流,即得對(duì)應(yīng)頻率10 Hz的熱流分量幅值,如圖9中的“Exact”示,可知,其幅值近似等于0.55。因此,當(dāng)測量誤差σ＝0.005Tmean時(shí),由前面分析知,該誤差將全部作用到10 Hz分量上,相當(dāng)于10 Hz分量對(duì)應(yīng)s＝0.005Tmean/(0.55Tm10Hz)情況,Tm10Hz值在之前已計(jì)算出Tm10Hz＝0.01457,測點(diǎn)平均溫升Tmean的值由式(1)計(jì)算出為Tmean＝0.1414。于是,當(dāng)σ＝0.005Tmean時(shí),s＝0.0880,在圖3中插值得出E3＝0.3396;代入式(9)可估算出此時(shí)對(duì)組合熱流的誤差為:

該值與采用式(3)和疊加64組同分布白噪聲的熱流辨識(shí)統(tǒng)計(jì)平均值0.2576符合較好。進(jìn)一步考慮測量誤差σ＝0.01Tmean的情況,可采用類似方法計(jì)算出s＝0.176,E3＝0.8196,E＝0.6541,與采用式(3)和疊加64組同分布白噪聲的熱流辨識(shí)統(tǒng)計(jì)平均值0.6176也符合較好,進(jìn)一步驗(yàn)證了估計(jì)方法的有效性。

綜上所述,對(duì)于一多頻率組合熱流,快速分析其辨識(shí)結(jié)果誤差的具體步驟為:

表2 典型頻率組合熱流的辨識(shí)結(jié)果平均誤差計(jì)算結(jié)果與分析預(yù)測結(jié)果對(duì)比Table 2 Comparison of calculated and predicted mean error for estimated multiple-frequency combined heat flux

圖8 待分析熱流及其功率譜Fig.8 Exemplified heat flux and its power spectrum

圖9 高頻熱流分量的幅值計(jì)算Fig.9 Amplitude calculation of high frequency component of heat flux

1)對(duì)待分析的熱流進(jìn)行頻譜分析,獲知最大頻率以及其余頻率的能量占比;

2)通過濾波計(jì)算獲得對(duì)應(yīng)最大頻率熱流分量的振幅AmaxHz;

3)針對(duì)最大頻率的周期熱流,計(jì)算其測點(diǎn)平均溫升TmaxHz及如圖3所示的“辨識(shí)結(jié)果誤差～測量誤差”對(duì)應(yīng)關(guān)系曲線;

4)由式(1)計(jì)算待分析熱流對(duì)應(yīng)的Tmean;

5)根據(jù)測量誤差及TmaxHz、Tmean、AmaxHz計(jì)算s;

6)由s 值在步驟(3)得出的“辨識(shí)結(jié)果誤差～測量誤差”關(guān)系曲線中插值得EN值;

7)由式(9)估算出待分析熱流的辨識(shí)結(jié)果誤差E。

4 結(jié) 論

本文對(duì)表面熱流的辨識(shí)誤差進(jìn)行了深入分析,首先對(duì)單一頻率的熱流辨識(shí)誤差建立了其與熱流頻率和測量精度之間的響應(yīng)面模型。然后重點(diǎn)對(duì)多個(gè)頻率組合情況下的辨識(shí)結(jié)果誤差進(jìn)行了定量分析。結(jié)果顯示,頻率組合熱流中的低頻分量在辨識(shí)結(jié)果中能較好地復(fù)現(xiàn),辨識(shí)結(jié)果誤差可只考慮高頻分量的辨識(shí)結(jié)果誤差。于是,由辨識(shí)誤差定義出發(fā),通過頻譜分析方法可計(jì)算出組合熱流中高頻分量的振幅及能量占比,進(jìn)而建立了組合熱流的辨識(shí)誤差的估算方法,并通過多個(gè)算例進(jìn)行了驗(yàn)證。這一結(jié)果一方面有較強(qiáng)的理論意義,揭示了頻率組合熱流與單一頻率熱流辨識(shí)誤差之間的內(nèi)在關(guān)聯(lián),深化了對(duì)辨識(shí)結(jié)果時(shí)域、頻域特性的認(rèn)識(shí),并且不僅對(duì)傳熱逆問題適用,對(duì)其余領(lǐng)域的逆問題誤差分析也有一定的參考價(jià)值。另一方面,這一結(jié)果在工程上也有實(shí)用價(jià)值,可減少“暴力”的仿真辨識(shí)計(jì)算,尤其是形成一些數(shù)據(jù)庫后甚至可不進(jìn)行數(shù)值計(jì)算,直接查表即可對(duì)辨識(shí)結(jié)果誤差進(jìn)行大致估算。

下一步,方法還需要在以下三個(gè)方面進(jìn)行完善和推廣,一是目前分析中熱流的組合頻率值較為稀疏,在頻域中表現(xiàn)為離散譜,便于開展分析,而對(duì)于頻譜范圍較寬、頻譜特性為連續(xù)譜的情況,需進(jìn)一步挖掘辨識(shí)計(jì)算的“停止準(zhǔn)則”與熱流功率譜之間的內(nèi)在關(guān)聯(lián),開展深入研究。二是溫度的采樣頻率對(duì)辨識(shí)誤差的影響。采樣頻率提高有利于熱流辨識(shí)精度的提高,但是采樣頻率與熱流辨識(shí)精度的定量關(guān)系還需進(jìn)一步分析。三是目前分析中噪聲形式只考慮了白噪聲形式,對(duì)于有色噪聲的情況還需要分析拓展。