晁奮進， 王少紅， 周福強， 徐小力

(北京信息科技大學 現(xiàn)代測控技術(shù)教育部重點實驗室，北京 100192)

0 引 言

爐膛溫度場是鍋爐燃燒過程中需要監(jiān)測的重要參數(shù)，直接關(guān)系到鍋爐的燃燒安全與效率，影響污染 物的生成和排放量[1]。鍋爐爐膛由燃燒室、煙囪、圈管、燃燒機等組成，如圖1 所示，是非常復雜的工業(yè)設備，廣泛應用于各種領域。大型鍋爐由于燃燒不均勻會產(chǎn)生很多安全事故，因此監(jiān)測爐內(nèi)溫度場分布對于更好預防安全事故發(fā)生具有重要意義。超聲波測溫技術(shù)是一種高效的非接觸測溫技術(shù)，利用超聲波測溫技術(shù)對爐膛進行溫度監(jiān)測，將大大減少由接觸式測溫帶來的不便。

圖1 鍋爐結(jié)構(gòu)與原理圖

超聲波測溫是近年來國內(nèi)外研究的熱點，國內(nèi)外學者對超聲測溫技術(shù)的研究主要集中在高溫爐內(nèi)氣體溫度場的測量、醫(yī)學超聲熱療法、流體溫度和超聲測溫裝置等方面的研究。開發(fā)新型燃燒氣體超聲測溫技術(shù)，OUK H[2]等通過瑞利指數(shù)或火焰?zhèn)鬟f函數(shù)等信息協(xié)同提高燃燒不穩(wěn)定性診斷的準確性。安連鎖、沈國清等[3-4]運用不同算法重建爐膛內(nèi)三維溫度場，驗證了超聲波測溫的可行性，但并沒有對飛行時間進行準確的分析。耿帥、甄麗等[5-6]基于有限元法模擬了聲場能量分布情況，但不能直觀的反映溫度場的分布情況。

此前很多研究根據(jù)測量的聲波傳播時間重建爐內(nèi)溫度場，這是典型“反演”問題[7]，不能直觀的觀察聲波在爐膛中的傳播情況，本研究通過有限元仿真[8]，研究聲波在爐內(nèi)溫度場中的傳播時間，并實現(xiàn)爐膛內(nèi)部溫度場可視化。由于不具備在現(xiàn)場布置傳感器的條件，因此采用仿真研究進行模型溫度場的重建。COMSOL 具有很好處理多物理場耦合的能力[9]，采用聲、電、結(jié)構(gòu)模塊進行超聲波在爐膛溫度場模型中的傳播仿真，根據(jù)互相關(guān)算法計算出超聲波換能器之間的飛行時間，再由溫度場重建算法重建出爐膛內(nèi)溫度場的分布情況。

1 超聲波測溫原理

氣體介質(zhì)中聲波的傳播速度是該氣體介質(zhì)溫度的函數(shù), 同時與該氣體的組分有關(guān)。在工程應用上，通?？梢哉J為氣體溫度是聲速的第一函數(shù)[10]，如下式所示：

式中：T——氣體溫度， K；

c——聲波在介質(zhì)中傳播速度；

M——氣體摩爾質(zhì)量，kg/mol；

κ——氣體絕熱指數(shù)，本文介質(zhì)為空氣，且不考慮溫度變化的影響， κ取1.4；

R——理想氣體常數(shù)，取8.314 5 J/(mol·K)，與氣體狀態(tài)和種類無關(guān)。

由(1)式可知，對于特定的介質(zhì)，聲速只與溫度有關(guān)。對于單路徑測溫，在待測溫度場周圍對稱安裝兩個超聲波換能器，一個做為聲波發(fā)射器，一個做為聲波接收器，由于兩個超聲換能器之間的距離l已知，再通過互相關(guān)算法計算出換能器之間聲波的飛行時間t，即可計算出聲波在此路徑上的傳播速度c=l/t，再由(1)知，進而計算出聲波傳播路徑上的平均溫度。針對本文研究的多路徑測溫，在爐膛截面四周均勻布置超聲波換能器，每個聲波既是發(fā)射器又是接收器，根據(jù)測溫原理，計算出每條聲波傳播路徑的平均溫度之后，再由溫度場重建算法，即可重建出爐膛溫度場分布情況。多路徑測溫示意圖如圖2 所示。

圖2 聲學法測溫示意圖

2 有限元模型建立

2.1 模型描述

根據(jù)爐膛幾何特點，選取爐膛截面作為幾何模型，截面半徑為3 m。假定爐膛內(nèi)溫度場為單峰偏斜溫度場，爐膛四周對稱安裝八個超聲波換能器1、2、3、4、5、6、7、8，如圖3 所示，采用對稱布置，有利于聲波信號的接收，并且有利于提高超聲換能器之間距離l計算的精度。

圖3 爐膛二維平面模型

爐膛模型溫度分布為單峰偏斜溫度場T(x,y)：

聲波在爐膛氣體介質(zhì)溫度場中的傳播由下式得：

式中： ±——時間， s；

v——環(huán)境流速，m/s；

cmf——聲速，m/s；

ρ——爐膛介質(zhì)密度，kg/m3；

φ——聲波的速度勢，m2/s。

2.1.1 邊界選擇

爐膛外壁采用金屬材料，由于發(fā)射波和反射波交叉會影響原始聲波信號，因此選用能夠充分吸收壁面聲波的邊界條件——平面波輻射[11]方程為：

式中：

2.1.2 網(wǎng)格劃分和求解

有限元模型的建立過程中，網(wǎng)格形狀和網(wǎng)格尺寸對計算準確度有直接的影響，為了保證網(wǎng)格幾何形狀的合理性，選取自由三角形網(wǎng)格作為爐膛截面的網(wǎng)格劃分形式，四個換能器選用自由四邊形網(wǎng)格，針對網(wǎng)格劃分大小，可以選取COMSOL 軟件提供的細化、粗化、極細化等多種網(wǎng)格劃分尺寸的方式。研究聲波在爐膛中的傳播問題，模型網(wǎng)格劃分的大小選用極細化的方式。

在爐膛空氣域的數(shù)值模擬過程中，時間步長CFL是計算收斂的條件，見式(5)，直接影響計算超聲波飛行時間的準確度，在COMSOL 平臺上進行聲學問題研究時，為了保證飛行時間計算精度的準確性，為了能夠加速收斂，節(jié)省仿真計算時間，需要保證CFL小于1[12]。

式中： ?t——時間步長；

c——聲速， m/s；

h——最小網(wǎng)格尺寸。

2.2 飛行時間數(shù)值模擬

為了驗證使用COMSOL 求解飛行時間的準確性，超聲波換能器采用一發(fā)多收的方式，換能器1發(fā)射聲波，換能器3、4、5、6、7 接收聲波，換能器2發(fā)射聲波，換能器4、5、6、7、8 接收聲波，換能器3發(fā)射聲波，換能器4、1、2 接收聲波。圖4 為換能器1 和2 發(fā)射聲波時的波形圖。

圖4 脈沖波發(fā)射圖

2.3 飛行時間計算

2.3.1 互相關(guān)算法

飛行時間的不同反映了兩換能器之間的距離和溫度的差異，由于模型中換能器之間距離l已知，在計算出飛行時間之后，利用超聲波測溫原理就能得到單一路徑上的平均溫度，互相關(guān)算法的測量精度高[13]，利用互相關(guān)算法即可求得換能器之間的飛行時間。

互相關(guān)函數(shù)表示的是兩個時間序列之間的相關(guān)程度，描述信號x(t) ，y(t)在任意兩個不同時刻t1，t2的取值之間的相關(guān)程度。x(t)，y(t)兩個信號之間互相關(guān)函數(shù)表達式為：

當互相關(guān)函數(shù)R12取得極值時，此時對應的時間為兩個不同超聲波換能器之間的聲波飛行時間。

2.3.2 互相關(guān)計算結(jié)果

在發(fā)射和接收期間，發(fā)射和反射聲波信號是周期性的，也有反射波信號，為了保證飛行時間計算的準確性，選取換能器第一次接收到的信號來計算飛行時間。當換能器1 發(fā)射信號，換能器3、4、5、6、7 來接收信號，飛行時間結(jié)果如圖5 所示。

圖5 各換能器飛行時間圖

根據(jù)表1 中計算的結(jié)果，COMSOL 仿真結(jié)果和理論計算結(jié)果得出的飛行時間相對偏差最大值為0.81%，結(jié)果表明，通過互相關(guān)算法計算得出的飛行時間和理論計算的飛行時間基本一致。

表1 不同傳播路徑飛行時間仿真計算和理論計算值及其偏差

圖5 為當換能器1 作為發(fā)射端，換能器3, 4, 5,6, 7 作為接收端時的飛行時間傳播圖，由（6）式可知，當電壓信號第一次取得極值時，兩信號之間的時間間隔就為飛行時間。

3 最小二乘法的超聲波測溫原理

飛行時間是針對于某一傳播路徑，從超聲波傳感器發(fā)射聲波信號開始，到該路徑接收信號的超聲波傳感器第一次接收聲波信號的時間間隔。記為t：

式中：c——聲速；

t——飛行時間；

l——傳播路徑距離；

dl——聲波傳播路徑線性微分元；

a——聲波傳播路徑上聲速的倒數(shù)。

根據(jù)爐膛二維模型的區(qū)域劃分和聲波傳播路徑情況，由于爐膛內(nèi)單峰偏斜溫度場每個子區(qū)域i的溫度是均勻梯度分布的，所以某一特定路徑的超聲波飛行時間可以表示為：

式中：?Ski——路徑lk在子區(qū)域i的區(qū)間長度；

Ai—— 子區(qū)域i超聲波傳播平均速度的倒數(shù)，i=1,2,3,···,M，k=1,2,3,4,···,N。

在第k條路徑上的理論飛行時間tk1與實際超聲波飛行時間tk2之間的誤為 εk：

記S為飛行路徑與各測溫區(qū)間關(guān)系矩陣：

記A為超聲波在待測區(qū)域中的空間特性矩陣：

記t為超聲波飛行時間矩陣：

根據(jù)最小二乘法原理，令

取最小值，可解得

再由聲速和溫度之間的關(guān)系，根據(jù)上式得到每一子區(qū)域的平均溫度為：

在被測區(qū)域建立直角坐標系，計算得到每個小區(qū)間的平均溫度，再運用插值算法，得到被測爐膛區(qū)域的溫度場可視化二維平面圖。

4 仿真驗證

為了更直觀的觀察換能器安裝和區(qū)域劃分情況，在半徑為3 m 的二維爐膛區(qū)域，按如圖6 所示方式在爐膛四周均勻布置8 個換能器，待測區(qū)域被劃分為13 個待測區(qū)間，形成20 條有效聲波傳播路徑。

圖6 二維模型區(qū)域劃分圖

在鍋爐爐膛模型中，由于模型內(nèi)部介質(zhì)為氣體，由式(12)可知，可以求出13 個區(qū)域的平均溫度，再由插值算法，即可得到爐膛模型的二維等溫線圖。單峰偏斜模型溫度場和重建溫度場如圖7 所示。

圖7 單峰偏斜模型溫度場與重建溫度場結(jié)果圖

5 誤差分析

由于實際爐膛環(huán)境比較復雜，有空氣、煤粉以及煤粉燃燒產(chǎn)物，介質(zhì)多樣且不均勻，本文為模擬聲波在爐膛中準確傳播情況做了簡化處理，只研究聲波在靜態(tài)溫度場中的傳播特性，選取COMSOL材料庫中的空氣作為爐膛內(nèi)部空間的材料。

由表1 可知，利用互相關(guān)算法計算的仿真飛行時間比理論計算值均偏大，且其誤差范圍均在1%以內(nèi)，是因為聲源信號存在偽前移，由于聲源處網(wǎng)格密度大，所以偽前移現(xiàn)象明顯，從而導致誤差的產(chǎn)生。

另外邊界條件、聲源選擇、網(wǎng)格劃分、時間步長都是直接影響溫度場重建準確性的因素。為了盡可能準確的重建溫度場，換能器發(fā)出的聲波需具有較強的能量，因此聲源形式采用高斯脈沖波，頻率為4 kHz。邊界條件采用平面波輻射，網(wǎng)格數(shù)量、網(wǎng)格密度、單元階次、單元形狀和單元協(xié)調(diào)性都是評價網(wǎng)格劃分質(zhì)量好壞的重要指標，兼顧這些因素，網(wǎng)格劃分采用自由三角形網(wǎng)格，時間步長如式（5）所示，需小于1。

針對溫度場重建質(zhì)量問題，本文通過使用最大絕對誤差Emax、最小絕對誤差Emin、平均絕對誤差Emean三個誤差數(shù)據(jù)來評估溫度場重建效果，誤差函數(shù)定義如下：

單峰偏斜溫度場仿真重建誤差見表2。

表2 溫度場重建誤差分析

6 結(jié)束語

本文通過理論分析與仿真模擬相結(jié)合的方式，基于COMSOL 平臺，建立鍋爐爐膛二維物理模型，模擬聲波在單峰偏斜溫度場的傳播情況，直觀的觀察到了聲波在模型中的發(fā)射情況，利用互相關(guān)函數(shù)法求得了聲波傳播時間，運用最小二乘法和插值算法重建出溫度場。

研究結(jié)果表明：

1）運用互相關(guān)算法計算聲波的飛行時間，理論結(jié)果和仿真結(jié)果的飛行時間最大相對誤差為0.81%，具有很高的吻合度。

2）依據(jù)超聲波測溫原理，基于互相關(guān)算法計算出的飛行時間，結(jié)合最小二乘法和插值算法，對爐膛單峰偏斜溫度場模型進行了還原，由誤差結(jié)果可知，最大絕對誤差為50 K，最小絕對誤差為17 K，平均絕對誤差為31 K，有較好的重建效果。

3）不同于溫度場“反演”研究，本文采用正問題求解方法，能更清晰理解溫度場重建過程，并能得到較好結(jié)果，為后來者提供了另外的可行思路，拓寬了求解渠道，期待超聲波測溫方法更多的應用于各個領域。