基于神經(jīng)網(wǎng)絡的表面熱流辨識三維效應修正

潘學浩陳偉芳彭玉酌楊 華

(1.浙江大學 航空航天學院,杭州 310027；2.航空工業(yè)第一飛機設計研究院,西安 710089)

0 引 言

表面熱流辨識屬于一類熱傳導反問題,它是通過測量導熱材料內(nèi)壁的溫度測點的溫度歷程,反演出外壁受熱面的熱流時間歷程。熱傳導反問題在航空航天、機械制造、車輛工程及生物工程等領(lǐng)域有著廣泛應用[1]。高超聲速飛行器在大氣層內(nèi)飛行時,面臨嚴重的氣動加熱問題,而對飛行器服役過程中的溫度、熱流等關(guān)鍵參數(shù)的測量是評價熱防護材料使用性能、驗證氣動熱模型和算法、指導熱防護設計的必要手段[2]。但對于駐點等熱流密度大的區(qū)域,通常不能直接安置溫度傳感器或熱流傳感器進行測量,一方面是由于結(jié)構(gòu)開孔導致結(jié)構(gòu)強度下降以及縫隙加熱,引起燒蝕不同步等結(jié)構(gòu)匹配問題[3]；另一方面,有些傳感器本體材料不能承受過高熱載荷,而且傳感器的嵌入帶來壁溫的不連續(xù)以及周向的干擾,導致熱流測量結(jié)果并非真實的氣動熱[4],所以需要采用表面熱流辨識手段監(jiān)測熱流。

國內(nèi)外對熱傳導反問題進行了大量的研究,通常做法是選取合適的目標函數(shù),將辨識問題轉(zhuǎn)化為優(yōu)化問題求解。美國在20世紀60年代的Reentry-F飛行試驗項目中,曾在飛行器彈身上布置了21個溫度測點,利用類似順序函數(shù)法的反演算法估算飛行器表面熱流,從而預測邊界層轉(zhuǎn)捩位置[5],與其它方法測得的轉(zhuǎn)捩位置基本一致。錢煒祺[6]分別用順序函數(shù)法和共軛梯度法研究了一維表面熱流辨識方法,并實現(xiàn)了二維和三維非穩(wěn)態(tài)表面熱流辨識[7-8]；薛齊文[9]應用Tikhonov方法研究了一維熱傳導反問題中,內(nèi)熱源強度、導溫系數(shù)及邊界條件的多宗量辨識；DENG[10]和智會強等[11]研究了神經(jīng)網(wǎng)絡和遺傳算法求解熱傳導反問題,但只考慮了簡單的熱流加載情況；CUI Miao[12]采用無量綱化目標方程,對熱流模型進行參數(shù)辨識,但局限于已知的熱流函數(shù)形式；錢煒祺[13]考慮到材料燒蝕后退,利用簡化后的熱解面燒蝕模型,對一維燒蝕表面熱流辨識進行了研究；張聰[14]利用簡化的一維和二維傳熱模型進行了高超聲速燃燒室壁面熱流的辨識,在軸對稱模型下取得了較好的效果。

在實際工程應用中,待辨識的熱流通常是隨時間和空間復雜變化的,已有的辨識方法不僅在三維情況下辨識精度不高,而且計算量龐大,難以滿足實時性要求。針對某些工程問題通常最關(guān)心重點區(qū)域駐點熱流隨時間的變化,本文提出神經(jīng)網(wǎng)絡和順序函數(shù)法結(jié)合的方法。在導熱材料內(nèi)壁布置適量溫度傳感器,對于每個傳感器測得的溫度數(shù)據(jù)采用一維的順序函數(shù)法得到對應的表面熱流辨識數(shù)據(jù),再以每個測點得到的熱流序列作為神經(jīng)網(wǎng)絡的輸入,以外壁受熱面對應時刻的駐點熱流序列作為網(wǎng)絡輸出,并引入粒子群算法優(yōu)化神經(jīng)網(wǎng)絡的初始權(quán)值和閾值,通過數(shù)值仿真或風洞試驗的方法獲得訓練樣本。該方法保留了順序函數(shù)法良好的抗噪性,同時神經(jīng)網(wǎng)絡本身有強大的非線性映射能力,適合用于三維效應的修正。算例驗證表明經(jīng)過訓練的神經(jīng)網(wǎng)絡預測的駐點熱流和實際值十分吻合。

1 非穩(wěn)態(tài)表面熱流辨識方法

一維、三維表面熱流辨識與二維情況類似,以二維為例,現(xiàn)介紹順序函數(shù)法的辨識方法。對于典型的二維矩形域,考慮非線性非穩(wěn)態(tài)熱傳導控制方程為:

邊界條件為:

初始條件:

非線性偏微分方程可以采用擬線性隱式差分或有限體積法求解。

假設溫度測點數(shù)量為M個,由于導熱過程的時間延遲性,測點t時刻溫度受到前r個時刻熱流密度的影響,已知t i前各時刻的熱流值,所以要辨識t i時刻的熱流q(t i,y),即要求如下目標函數(shù)極小:

其中:(x m,y m,t k)是測點m溫度測量值,T(x m,y m,t k)是對應位置溫度計算值。由于計算溫度場需要知道后續(xù)時刻熱流值,所以假設q(t i,y)到q(t i＋r－1,y)為線性變化,即:

順序函數(shù)法引入靈敏度函數(shù)的概念,即:

X(x,y,t)表示溫度場對熱流密度的靈敏度,靈敏度函數(shù)可以通過對熱傳導方程求偏導數(shù)得到:

邊界條件:

初始條件:

式(7)是一個近似方程,忽略了導熱系數(shù)和比熱對熱流的導數(shù)項,由于二者相對較小,出于簡化的目的可以忽略。

然后采用牛頓-拉夫遜算法得到待辨識熱流分量的迭代修正公式:

k為當前迭代次數(shù),γ為防止迭代發(fā)散的收斂因子。每步迭代中,由修正后的熱流可以計算得到新的目標函數(shù)值,當滿足:

可以認為迭代已經(jīng)收斂,ε為很小的正數(shù),此時的q(t i,y j)作為辨識值。沿時間推進,即可得到連續(xù)的熱流辨識序列。

順序函數(shù)法中未來時間步r的選擇直接影響到熱流辨識的精度和魯棒性,目前r的選取具有一定的經(jīng)驗性。為了選取最優(yōu)的未來時間步,文獻[15]提出了一種基于殘差原理的最優(yōu)未來時間步的確定方法,取不同的r,當滿足方程(12)時,此時的r則為最優(yōu)未來時間步。

為了避免r選擇不當帶來的影響,可以采用共軛梯度法(CGM)進行表面熱流的辨識,它將辨識問題分解為熱傳導求解、靈敏度求解和伴隨變量求解這三個問題進行求解,具體方法可參考文獻[6]。

2 順序函數(shù)法非穩(wěn)態(tài)熱流辨識算例

本文采用的算例計算域為:厚度L x＝0.01 m；寬度L y＝0.1 m；材料比熱c p＝500 J/(kg·K)；密度ρ＝8000 kg/m3；導熱系數(shù)k＝80－0.02T,W/(m·K)；初始溫度300 K。測溫點數(shù)量為21個,時間步長0.2 s,時間步r取25,實際加載的熱流函數(shù)為:

總加熱時間為t_total＝20 s,采用數(shù)值仿真的溫度結(jié)果作為實測值,為檢驗算法的抗噪性能,在溫度實測值基礎(chǔ)上分別加入方差D＝0.04和D＝0.25的零均值均勻分布白噪聲,真實熱流,辨識熱流以及加入噪聲之后的辨識熱流結(jié)果分別如圖1-4所示,Q表示熱流密度。無噪聲辨識相對誤差和加入噪聲后的辨識相對誤差分布分別如圖5～圖7所示。

可見,無噪聲時熱流辨識結(jié)果最大相對誤差約3.6%,加入D＝0.04和0.25的噪聲時,相對誤差最大值約為5.7%和8.0%。由于噪聲的增大,放大了反問題的不適定性,所以辨識結(jié)果相對誤差也會增大,但總體精度仍然保持在較高水平,證明了順序函數(shù)法用于非穩(wěn)態(tài)表面熱流辨識具有很好的準確性和抗噪性。

圖1 真實熱流結(jié)果Fig.1 The exact value of heat flux

圖2 無噪聲熱流辨識結(jié)果Fig.2 The identified value of heat flux without noise

圖3 噪聲方差0.04時的熱流辨識結(jié)果Fig.3 The identified value of heat flux with noise of D＝0.04

圖4 噪聲方差0.25時的熱流辨識結(jié)果Fig.4 The identified value of heat flux with noise of D＝0.25

圖5 無噪聲熱流辨識相對誤差分布Fig.5 The relative error of identified value without noise

圖6 噪聲方差0.04時的熱流辨識相對誤差分布Fig.6 The relative error of identified value with noise of D＝0.04

圖7 噪聲方差0.25時的熱流辨識相對誤差分布Fig.7 The relative error of identified value with noise of D＝0.25

3 人工神經(jīng)網(wǎng)絡三維效應修正模型

根據(jù)國際著名神經(jīng)網(wǎng)絡專家Hecht Nielsen的觀點,人工神經(jīng)網(wǎng)絡(Artificial Neural Networks,ANN)是由人工建立的、以有向圖為拓撲結(jié)構(gòu)的動態(tài)系統(tǒng),它通過對連續(xù)或斷續(xù)的輸入作狀態(tài)響應而進行信息處理[16]。為了使神經(jīng)網(wǎng)絡修正模型盡可能簡單穩(wěn)定,本文選擇目前應用廣泛的BP(Back Propagation)神經(jīng)網(wǎng)絡。它是一種多層前饋神經(jīng)網(wǎng)絡,該網(wǎng)絡拓撲結(jié)構(gòu)如圖8所示。

圖8 BP網(wǎng)絡拓撲結(jié)構(gòu)Fig.8 The topological structure of BP network

BP網(wǎng)絡主要特點是信號前向傳遞,誤差反向傳遞。輸入層和隱含層、隱含層和輸出層神經(jīng)元之前通過權(quán)值連接,隱含層和輸出層的每個神經(jīng)元分別有相應的閾值,隱含層和輸出層還要設置相應的激勵函數(shù)類型。一個典型的隱含層神經(jīng)元輸入輸出關(guān)系為:

其中,l為隱含層節(jié)點數(shù),n為輸入層節(jié)點數(shù),f、ω、a分別為激勵函數(shù)、權(quán)值和閾值。輸出層輸入輸出關(guān)系類似。

神經(jīng)網(wǎng)絡初始化后,需要輸入訓練數(shù)據(jù)樣本,若輸出不滿足誤差要求,則誤差轉(zhuǎn)入反向傳遞,調(diào)整相應的權(quán)值和閾值,直到輸出滿足期望輸出要求。經(jīng)過大量的數(shù)據(jù)訓練成熟的神經(jīng)網(wǎng)絡,將具有強大的非線性映射和泛化能力。

BP神經(jīng)網(wǎng)絡的特點,使其能夠用于表面熱流辨識三維效應的修正?；舅枷胧窃陲w行器表面峰值熱流出現(xiàn)的重點區(qū)域內(nèi)壁面,布置M個溫度傳感器,根據(jù)每個傳感器實時測量的溫度數(shù)據(jù),用一維順序函數(shù)法得到M個實時辨識的熱流結(jié)果,數(shù)據(jù)歸一化之后,將其作為神經(jīng)網(wǎng)絡的輸入,神經(jīng)網(wǎng)絡的輸出再反歸一化,就是此刻該區(qū)域?qū)姆逯禑崃鳌?/p>

4 粒子群算法優(yōu)化神經(jīng)網(wǎng)絡

遺傳算法(GA)和粒子群算法(PSO)都是已被廣泛應用的仿生全局優(yōu)化方法。遺傳算法模擬生物種群進化機制篩選最優(yōu)個體,粒子群算法模擬鳥群覓食行為逼近食物位置。鑒于粒子群算法相比于遺傳算法,不需要編碼且沒有交叉操作,使得模型更簡單,大多數(shù)情況下,粒子群的單項信息共享機制通常能更快地收斂于最優(yōu)解,故本文采用PSO算法優(yōu)化神經(jīng)網(wǎng)絡的初始權(quán)值和閾值,得到最優(yōu)的網(wǎng)絡修正模型。

PSO算法首先在問題的可解空間中將一群粒子初始化,用位置、速度和適應度來表示該粒子特征,適應度值是根據(jù)適應度函數(shù)計算得到的,代表著粒子的優(yōu)劣,每個粒子都代表問題的一個潛在解。粒子需要在解空間不停運動,其位置和速度的更新是根據(jù)個體最優(yōu)值和群體最優(yōu)值選擇的,通過不斷的更新群體,使群體最優(yōu)值逼近全局最優(yōu)解。

假設在一個D維的搜索空間中,n個粒子組成的種群為X＝(X1,X2,…,X n),X i＝(X i1,X i2,…,X iD)T代表第i個粒子在搜索空間中的的位置,粒子速度表示為V i＝(V i1,V i2,…,V iD)T,個體極值表示為P i＝(P i1,P i2,…,P iD)T,全局極值表示為P g＝(P g1,P g2,…,P g D)T。則粒子速度和位置更新公式為:

其中:ω為慣性權(quán)重；c1和c2為加速度因子,為非負常數(shù)；r1和r2是介于0和1之間的隨機數(shù)。通常要將粒子的位置和速度限制在一定區(qū)間,避免出現(xiàn)盲目搜索。為防止PSO早熟收斂,借鑒遺傳算法的變異思想,在粒子每次更新之后,以10%的概率重新初始化粒子。

在實際應用中,將神經(jīng)網(wǎng)絡權(quán)值和閾值排列為行向量,對應粒子位置。這樣每個粒子都對應一個神經(jīng)網(wǎng)絡,以訓練樣本的輸出誤差平方和作為粒子的適應度值,顯然適應度值越低,代表該粒子越優(yōu)。PSO算法最后得到的群體最優(yōu)值作為神經(jīng)網(wǎng)絡模型初始權(quán)值和閾值。

5 三維效應修正模型算例驗證

算例模型采用外徑0.2m,厚度0.03m的空心半球殼,如圖9所示。材料比熱c p＝500 J/(kg·K),密度ρ＝8000 kg/m3,導熱系數(shù)k＝80 W/(m·K),初始溫度300 K,時間步長0.2 s,總加熱時間200 s,時間步r取30,測溫點數(shù)量為9個,呈放射狀均布于內(nèi)壁面,布局在底面投影示意圖如圖10所示。測點1位于球殼端部,內(nèi)側(cè)2-5測點距離測點1約0.03 m,外側(cè)6-9測點距離測點1約0.06 m。

圖9 傳熱模型示意圖Fig.9 The heat transfer model

圖10 測溫點布局示意圖Fig.10 The layout diagram of measuring points

在球殼外表面施加隨時間和空間變化的熱流,熱流函數(shù)形式為:

R為任意點與對稱軸之間的距離,當t＝100 s時,沿對稱軸所在平面任意半圓周熱流分布如圖11所示。

圖11 t＝100 s時,熱流沿圓周分布Fig.11 The distributions of heat flux around circumference at t＝100 s

可見,該熱流形式隨空間變化明顯,且與飛行器表面熱流分布相似,可以較好地模擬氣動加熱。采用數(shù)值仿真軟件計算出9個測點的時變溫度數(shù)據(jù),并疊加標準差0.5K的零均值白噪聲作為后續(xù)輸入數(shù)據(jù),峰值熱流為:

在ANSYS APDL中加載(17)形式的熱流邊界條件,端部測點1的溫度隨時間變化如圖12所示。測點溫度先升高后降低,是因為材料內(nèi)部的熱擴散是三維的,表面熱流空間分布極不均勻,且隨時間先增大后減小,測點周圍單元的熱擴散從以徑向擴散升溫為主,逐步轉(zhuǎn)變?yōu)橐灾芟驍U散降溫為主。

圖12 測點1溫度隨時間變化曲線Fig.12 Time history of temperature at measuring point 1

先僅采用端部測點1的溫度數(shù)據(jù)進行一維辨識,結(jié)果如圖13所示,可見三維效應是很顯著的,一維辨識誤差很大,需要對其進行修正。

圖13 三維效應下的一維辨識結(jié)果對比Fig.13 Comparison of exact and estimated value for one-dimension with three-dimensional effect

將9個測點的一維熱流辨識結(jié)果作為神經(jīng)網(wǎng)絡輸入,峰值熱流序列作為網(wǎng)絡輸出。實際得到1000組輸入輸出數(shù)據(jù),隨機選取100個時刻測試數(shù)據(jù),剩下900組數(shù)據(jù)用于網(wǎng)絡訓練。網(wǎng)絡輸入層為9個神經(jīng)元,輸出層為1個神經(jīng)元,對于最佳隱含層節(jié)點數(shù)的選取,通常有如下經(jīng)驗公式:

n,l,m分別為輸入層,隱含層和輸出層節(jié)點數(shù),a為0～10之間的常數(shù),本文分別設置隱含層神經(jīng)元數(shù)目5～17。網(wǎng)絡隱含層和輸出層激發(fā)函數(shù)分別采用logsig和purelin函數(shù),最大迭代訓練次數(shù)為100,粒子群算法中,加速度因子設為1.49,慣性權(quán)重為0.5,粒子群規(guī)模為20個粒子,尋優(yōu)迭代次數(shù)為100。當隱含層節(jié)點數(shù)為16時,粒子群群體最佳適應度值隨尋優(yōu)次數(shù)變化如圖14所示,圖15為訓練后的網(wǎng)絡峰值熱流預測值與真實值比較結(jié)果,圖16為神經(jīng)網(wǎng)絡預測誤差。將測試樣本輸出按時間順序排列后與真實熱流對比結(jié)果如圖17所示。

圖14 粒子群群體最佳適應度值Fig.14 The best fitness value of particle swarm

圖15 網(wǎng)絡峰值熱流預測結(jié)果Fig.15 The estimated value of peak heat flux

圖16 網(wǎng)絡峰值熱流預測相對誤差Fig.16 The relative error of estimated value of peak heat flux

圖17 峰值熱流預測結(jié)果對比Fig.17 Comparison of exact and estimated value of peak heat flux

可見經(jīng)過PSO算法優(yōu)化的BP網(wǎng)絡預測精度極高,隱含層節(jié)點16時,網(wǎng)絡最大預測相對誤差僅有0.20%,表1是不同隱含層節(jié)點數(shù)時,網(wǎng)絡預測最大相對誤差。

表1 不同隱含層節(jié)點數(shù)網(wǎng)絡預測最大相對誤差Table 1 The maximum relative error of estimated value with different hidden layer nodes

可以看出,隱含層節(jié)點數(shù)在較大范圍變化時,神經(jīng)網(wǎng)絡預測精度仍然能保持在較高水平,隱含層節(jié)點數(shù)9時預測精度最優(yōu)。

在工程應用中,通常希望在達到預期辨識精度前提下盡量減少測點數(shù)量。本文只選取1、4、6三個測點溫度數(shù)據(jù),設計神經(jīng)網(wǎng)絡結(jié)構(gòu)為3輸入,單輸出,隱層節(jié)點數(shù)為7,相同數(shù)量樣本訓練后的預測結(jié)果和相對誤差如圖18和圖19所示。

當測點數(shù)目減少到三個時,峰值熱流預測最大相對誤差為2.3%,雖然誤差略有增大,但預測精度仍然較好,所以工程應用中可以適當減少測點數(shù)量以降低成本。

由于熱流的時空分布形式不一,特別是隨時間變化形式多樣,需要驗證基于一種熱流時空分布建立的

圖18 三個測點峰值熱流預測結(jié)果Fig.18 The estimated value of peak heat flux with three measuring points

圖19 三個測點峰值熱流預測相對誤差Fig.19 The relative error of estimated value of peak heat flux with three measuring points

修正模型,在另外一種區(qū)別較大的熱流加載下的表現(xiàn)。因此,本文根據(jù)半球殼在方程(17)的熱流工況,建立并訓練了BP網(wǎng)絡修正模型(12個隱含層節(jié)點),然后設計了具有不同時空分布形式的熱流,其函數(shù)形式如式(22):

其中:R為任意點與對稱軸之間的距離,t為時間。

同樣采用有限元軟件計算出9個測點的時變溫度數(shù)據(jù),并疊加標準差0.5K的零均值白噪聲作為后續(xù)輸入數(shù)據(jù),將一維的熱流辨識結(jié)果輸入前文已經(jīng)建立的BP網(wǎng)絡修正模型,駐點熱流預測結(jié)果與實際值比較如圖20所示。

可見當應用于不同的熱流工況時,修正模型的駐點熱流預測值可以較好地跟隨實際熱流的變化趨勢,但精度會顯著降低,但相對于一維的辨識結(jié)果,仍然有了極大的提高。本文算例采用的材料較厚,測點離表面較遠,同時模型頭部半徑較小,這些導致傳熱的三維效應非常明顯,熱流辨識的不適定性增加,這些也是圖20對比結(jié)果區(qū)別明顯的原因。

圖20 駐點熱流預測結(jié)果對比Fig.20 Comparison of exact and estimated value of peak heat flux

6 結(jié) 論

本文提出神經(jīng)網(wǎng)絡和一維順序函數(shù)法結(jié)合的方法,構(gòu)建表面熱流辨識三維效應修正模型,由內(nèi)壁測點的溫度測量數(shù)據(jù)直接獲得重點區(qū)域表面熱流峰值實時辨識結(jié)果。通過數(shù)值仿真的算例測試結(jié)果可以看出,本文提出的方法對于峰值熱流的辨識結(jié)果準確度高,模型訓練可在線下進行,避免了三維辨識耗時大、精度低、無法實現(xiàn)實時在線辨識的缺點,同時具有良好的抗噪性和穩(wěn)定性,在航天器表面峰值熱流在線實時辨識中有廣闊的應用前景。