李文軍， 王瀟楠， 鄭永軍

(中國計量大學(xué) 計量測試工程學(xué)院, 浙江 杭州 310018)

0 引言

熱電偶被廣泛應(yīng)用于工業(yè)測量領(lǐng)域。隨著工業(yè)測量技術(shù)水平的提高，測量瞬態(tài)溫度的需求也日益迫切和廣泛[1-2]。在已有的測量瞬態(tài)溫度方法中，常用的方法是利用多偶組合實現(xiàn)補償式瞬態(tài)測溫，以及利用實驗建模技術(shù)近似得到系統(tǒng)模型實現(xiàn)瞬態(tài)測溫。

Enrico等[3]從理論和實驗兩個角度，對3個溫度傳感器組合測量瞬態(tài)溫度和進行動態(tài)校準(zhǔn)做了對比，討論了產(chǎn)生精確溫度梯度的問題以及保持均勻穩(wěn)定溫度場的問題。吳朋等[4]對于鎧裝熱電偶傳感器測量瞬態(tài)溫度問題，采用支持向量機方法辨識建模，結(jié)合量子遺傳算法對核函數(shù)參數(shù)進行優(yōu)化以減小建模誤差，對鎧裝熱電偶傳感器測量瞬態(tài)溫度做精度補償。Villafae等[5]借助計算流體動力學(xué)，用數(shù)值實驗方法描述了熱電偶內(nèi)部熱擴散的演化過程，分析了外部特定測量環(huán)境下動態(tài)測量誤差的產(chǎn)生原因。Doghmane等[6]采用激光器對K型熱電偶測量接點產(chǎn)生激勵并獲得熱電偶的響應(yīng)，分別利用Yag激光器加熱獲得熱電偶在單位脈沖激勵下的響應(yīng)，利用Argon激光器加熱獲得熱電偶在周期性激勵下的響應(yīng)。郝曉劍等[7]針對鎧裝K型熱電偶動態(tài)響應(yīng)誤差補償，采用CO2激光器加熱產(chǎn)生激勵，并引入紅外探測器做輔助測量，獲得熱電偶的頻率特性對熱電偶進行補償。楊慶濤等[8]針對高超聲速飛行器地面試驗中測量壁面溫度的需求，建立了熱電偶傳感器有限元數(shù)值模型，通過傳感器內(nèi)部傳熱計算，考慮溫差項和儲能項，分析了熱電偶響應(yīng)特性。上述工作在分析熱電偶動態(tài)特性時，都采用了1階或2階整數(shù)階傳熱模型，而整數(shù)階導(dǎo)熱模型中的時間導(dǎo)數(shù)由局部極限定義，表示溫度在某時刻的變化，無法表征整個溫度的變化歷史。

從熱電偶與被測介質(zhì)之間的傳熱過程看，由于熱電偶測量接點存在熱慣性，且接點與被測介質(zhì)之間存在接觸熱阻和復(fù)雜換熱[9]，測量接點與熱電偶偶絲之間也存在熱量傳遞[10]，故熱電偶測量瞬態(tài)溫度的過程構(gòu)成了一個物理上的延遲過程或者叫做遺傳過程[11]。根據(jù)分數(shù)階微積分理論，分數(shù)階微積分具有記憶特點[12]，采用分數(shù)階微積分建立模型能更準(zhǔn)確地描述過程的動態(tài)特性[13-14]。

本文利用分數(shù)階微積分建立了熱電偶偶絲的時間分數(shù)階導(dǎo)熱模型，描述了露端式熱電偶測量接點與被測介質(zhì)之間的換熱過程。給出了測量固體介質(zhì)表面的瞬態(tài)溫度時模型參數(shù)的理論值，利用實驗系統(tǒng)采集了熱電偶輸入和輸出。最后通過實驗數(shù)據(jù)建模，用改進隨機數(shù)直接搜索法對傳遞函數(shù)的參數(shù)以及階次進行了估計。

1 熱電偶時間分數(shù)階熱傳導(dǎo)模型

1.1 分數(shù)階微積分定義和分數(shù)階拉普拉斯變換

(1)

式中：f(t)為t的函數(shù)；Γ(·)為Γ函數(shù)；n為大于μ的最小整數(shù)；τ為被積變量。

設(shè)算子下限a=0，對于零初始狀態(tài)，Riemann-Liouville分數(shù)階微積分對應(yīng)的拉普拉斯變換為

(2)

式中：L表示拉普拉斯變換；s為拉普拉斯算子；F(·)為s的函數(shù)。

(3)

(1)式、(2)式和(3)式給出了分數(shù)階微積分算子和對應(yīng)的拉普拉斯變換。

1.2 熱電偶半無限固體傳熱模型

熱電偶由一對有共同接點的材料互異的熱電偶偶絲組成，如圖1所示。圖1中，熱電偶偶絲的直徑為d，熱電偶偶絲1與熱電偶偶絲2之間存在一個夾角[16]。

由于熱電偶偶絲長度相對于其直徑而言很長，且兩根熱電偶偶絲之間相互絕熱，單根熱電偶偶絲導(dǎo)熱可以視為半無限固體導(dǎo)熱問題[17]。熱電偶偶絲的半無限固體導(dǎo)熱模型如圖2所示，圖2中：q(t)為熱流密度(W/m2)；x為發(fā)生熱傳導(dǎo)方向上的長度(m)；Ts(t)為半無限固體溫度分布函數(shù)T(t,x)在x=0處的特例。

熱電偶偶絲的導(dǎo)熱視為半無限固體導(dǎo)熱，其熱傳導(dǎo)方程為

(4)

式中：c為熱電偶偶絲材料比熱(J/(kg·K))；ρ為熱電偶偶絲材料密度(kg/m3)；k為導(dǎo)熱系數(shù)(W/(m·K))；T0為初始狀態(tài)溫度(K)。

引入過余溫度θ(t,x)=T(t,x)-T0，并代入(4)式，得到：

(5)

(6)

當(dāng)x→-∞， (6)式的邊值解為

(7)

(7)式對x微分并取x=0，有

(8)

根據(jù)(2)式和(3)式，對(8)式做分數(shù)階拉普拉斯反變換，得到時間分數(shù)階微分形式的關(guān)系式為

(9)

根據(jù)所定義的過余溫度，有

θ(t,0)=T(t,0)-T0=Ts(t)-T0.

(10)

將(10)式代入(9)式并整理得到：

(11)

(11)式等號左側(cè)即熱流密度q(t)：

(12)

再定義熱擴散系數(shù)α(m2·s)：

(13)

并記熱電偶測量接點溫度與初始狀態(tài)溫度的差Tb(t)=Ts(t)-T0，則(11)式化為

(14)

(14)式給出了半無限固體傳熱模型下熱電偶偶絲熱流密度與溫度的時間分數(shù)階微分關(guān)系式。

1.3 露端式熱電偶測量接點的傳熱方程

當(dāng)利用露端式熱電偶測量介質(zhì)溫度時，傳熱過程如圖3所示，圖3中：Tg(t)為被測介質(zhì)溫度與初始狀態(tài)溫度的差；Qi(t)為被測介質(zhì)到露端式熱電偶測量接點的熱流量；Q1(t)和Q2(t)分別為熱電偶偶絲1、熱電偶偶絲2的熱流量；A1、k1、α1分別為熱電偶偶絲1的截面積、導(dǎo)熱系數(shù)和熱擴散系數(shù)；A2、k2、α2分別為熱電偶偶絲2的截面積、導(dǎo)熱系數(shù)和熱擴散系數(shù)。

根據(jù)(14)式，分別有

(15)

(16)

露端式熱電偶測量接點的換熱方程表示[18]為

(17)

式中：h為熱電偶測量接點與被測流體之間的換熱系數(shù)(被測介質(zhì)為流體)或者接觸系數(shù)(被測介質(zhì)為固體)(W/(m2·K))；A為熱電偶測量接點與被測物體之間的接觸面積(m2)；m為熱電偶測量接點質(zhì)量(kg)。

(17)式做拉普拉斯變換并整理，得到：

(18)

(18)式給出了熱電偶時間分數(shù)階模型下的傳遞函數(shù)，其中s的階次μ1、μ2分別為1.0和0.5，它是一種固定階次的分數(shù)階傳遞函數(shù)。由于熱電偶在到達穩(wěn)態(tài)溫度過程中，并不是一個嚴格的絕熱過程，s的階次μ1、μ2之間并不一定存在嚴格比例關(guān)系[19]，傳遞函數(shù)可表示為更一般的形式：

(19)

式中：參數(shù)a1為瞬態(tài)傳熱過程中測量接點熱慣性引起的延遲；a2為瞬態(tài)傳熱過程中溫度激勵引起的熱擴散中單相延遲[20]，也就是熱流梯度相比溫度梯度的延遲；a3為常數(shù)項。

2 傳遞函數(shù)參數(shù)的理論區(qū)間

2.1 熱電偶的熱物性參數(shù)和幾何參數(shù)

對于(19)式，在μ1=1.0和μ2=0.5時，參數(shù)a1、a2的理論值可表示為

(20)

(21)

在被測介質(zhì)溫度變化幅值較小時，測量接點的導(dǎo)熱系數(shù)、密度和比熱k、ρ、c都可近似視為常數(shù)[21]。以一種常用的鎳鉻- 鎳硅熱電偶為例，表1給出了熱電偶偶絲的熱物性參數(shù)。

表1 熱電偶偶絲的熱物性參數(shù)

表2 熱電偶測量接點的等效熱物性參數(shù)

表3 熱電偶偶絲的幾何參數(shù)

表4 熱電偶測量接點的等效幾何參數(shù)

2.2 固定階次分數(shù)階模型參數(shù)的理論區(qū)間

熱電偶測量接點與被測固體表面的接觸可以簡化為常壓下間隙介質(zhì)為空氣的金屬表面間接觸[24]，其接觸系數(shù)h理論值介于1 000 W/(m2·K)與5 000 W/(m2·K)之間[25]。根據(jù)(20)式，并代入表2中熱電偶測量接點熱物性等效參數(shù)和表4中熱電偶測量接點幾何等效參數(shù)以及接觸系數(shù)理論值，得到a1的理論值區(qū)間。根據(jù)(21)式并代入表1中熱物性參數(shù)和表3中幾何參數(shù)以及接觸系數(shù)理論值，得到a2的理論值區(qū)間。計算結(jié)果見表5.

表5 參數(shù)的理論值

雖然理論值存在誤差，但是當(dāng)利用改進隨機數(shù)直接搜索法等方法估計參數(shù)時，上述理論值仍然提供了計算所必須的初始值取值區(qū)間。

3 實驗系統(tǒng)和實驗數(shù)據(jù)

3.1 實驗系統(tǒng)

為了獲得熱電偶在溫度激勵下的動態(tài)響應(yīng)，建立了一種實驗系統(tǒng)，結(jié)構(gòu)如圖4所示，主要包括計算機、美國Measurement Computing公司生產(chǎn)的MCC-USB-2408數(shù)據(jù)采集卡、日本歐姆龍公司生產(chǎn)的CPM1A可編程控制器、機械往復(fù)運動機構(gòu)以及寧波元創(chuàng)機電公司生產(chǎn)的M12 NPN激光對射開關(guān)和日本安立公司生產(chǎn)的ACSⅡ-2000溫度標(biāo)準(zhǔn)源恒溫面。

實驗原理是將熱電偶固定在運動機構(gòu)上，由可編程控制器控制運動機構(gòu)，驅(qū)動熱電偶進入定位點，使得熱電偶測量接點與恒溫面接觸并保持，熱電偶與恒溫面的接觸時間用激光對射開關(guān)采集，熱電偶輸出信號由數(shù)據(jù)采集卡采集和記錄。實驗中先對熱電偶做靜態(tài)響應(yīng)標(biāo)定，再參考溫度傳感器動態(tài)響應(yīng)校準(zhǔn)規(guī)程，最后用實驗系統(tǒng)獲得熱電偶在激勵下的輸出。

3.2 實驗數(shù)據(jù)

實驗設(shè)定的激勵溫度區(qū)間設(shè)定為300～400 K之間，表6給出了實驗時的3組實驗配置數(shù)據(jù)。

表6 實驗配置

表6的實驗配置對應(yīng)實驗數(shù)據(jù)集散點圖如圖5所示。數(shù)據(jù)集時間步長為0.25 s.

對3組數(shù)據(jù)做歸一化處理并平均，取激勵發(fā)生時間點為起始時間，截取后的數(shù)據(jù)集長度為101，測量值平均值的散點圖如圖6所示。

上述測量值的平均值可作為待辨識數(shù)據(jù)集。為了簡便，以下討論中測量值的平均值仍用測量值指代。

4 熱電偶分數(shù)階傳遞函數(shù)辨識

4.1 分數(shù)階傳遞函數(shù)的辨識方法

分數(shù)階傳遞函數(shù)的辨識方法有多種，其中改進Luus-Jaakola算法是Luus-Jaakola算法的優(yōu)化，它是一種基于隨機搜索的優(yōu)化方法[26]。采用這種方法，首先在參數(shù)搜索區(qū)間內(nèi)采用隨機搜索的方式獲取區(qū)間次優(yōu)值；然后在次優(yōu)值的基礎(chǔ)上自適應(yīng)改變搜索區(qū)間，做下一步隨機搜索，最終達到最優(yōu)結(jié)果。這種方法的關(guān)鍵步驟是采用一個可變收縮系數(shù)，在每次迭代后，用收縮系數(shù)縮小搜索空間以減少搜索時間并達到搜索目標(biāo)。改進Luus-Jaakola算法的步驟如下：

1)輸入要搜索的向量β(模型參數(shù)和階次)；

2)確定評估函數(shù)J；

3)產(chǎn)生隨機搜索向量初始值β0；

4)迭代計算，直到誤差滿足要求或者迭代次數(shù)達到限定值；

5)選擇收縮系數(shù)φ，將搜索區(qū)間乘以收縮系數(shù)；

6)反復(fù)迭代到滿足評估函數(shù)J所確定的性能指標(biāo)；

在本文算例中，評估函數(shù)J采用輸出誤差平方和形式[27]為

(22)

式中：M為數(shù)據(jù)集長度；y為預(yù)報值；y0為測量值；i為離散時間變量，i=0，1，2，…，M.

4.2 改進Luss-Jaakola算法的驗證

為了驗證改進Luss-Jaakola算法，首先對于已知的模型結(jié)構(gòu)進行辨識。假設(shè)已知模型為

(23)

根據(jù)文獻[28]提出的Luss-Jaakola算法驗證方法，首先由(23)式生成真實值數(shù)據(jù)集，并增加高斯分布的隨機噪聲，得到一組模擬測量值數(shù)據(jù)集；然后根據(jù)模擬測量值數(shù)據(jù)集，用改進Luss-Jaakola算法進行辨識，得到參數(shù)和階次的辨識結(jié)果。

上述驗證算例中，生成的模擬測量值數(shù)據(jù)集離散時間步長為0.125 s，數(shù)據(jù)集長度為201. 定義待搜索參數(shù)a1、a2、a3和階次μ1、μ2構(gòu)成的向量為

β=[a1a2a3μ1μ2],

向量真實值為

β=[1.20 1.30 1.00 2.10 0.50],

辨識結(jié)果如表7所示。表7中真實值為已知模型給出的參數(shù)和階次，估計值為通過模擬測量值辨識得到的參數(shù)和階次。

表7 參數(shù)真實值和估計值

輸出的模擬測量值、分數(shù)階模型估計值和真實值的比較如圖7所示，評估函數(shù)J為0.041.

從表7和圖7可以看到，Luaas-Jaakola算法的辨識結(jié)果較好地逼近了真實值，表明方法有效。

4.3 固定階次分數(shù)階傳遞函數(shù)參數(shù)的辨識

圖7所示的測量值數(shù)據(jù)集，其時間步長為0.25 s，數(shù)據(jù)集長度為101. 根據(jù)(19)式，模型階次已經(jīng)確定，是一種呈比例階次的分數(shù)階結(jié)構(gòu)，其形式為

(24)

(25)

式中：參數(shù)a1為2.24，位于理論值區(qū)間[0.48, 2.38]內(nèi)，表明在固定階次模型下，引入的測量接點熱慣性項與理論估計較為一致，能夠刻畫熱慣性引起的延遲；參數(shù)a2為0.10，偏離了理論值區(qū)間[0.14, 0.72]，小于理論值區(qū)間的最小值0.14，表明在固定階次模型下，引入的熱電偶絲熱擴散項與理論估計不一致，雖然表征了熱電偶絲熱擴散引起的延遲，但是只是延遲的一種弱化表示。輸出測量值和固定階次模型估計值的比較如圖8所示。

從圖8可以看出，固定階次分數(shù)階模型在時間軸0～4 s區(qū)間時與實驗數(shù)據(jù)不能很好地吻合。從物理過程分析，0～4 s區(qū)間內(nèi)，激勵溫度與環(huán)境溫度差較大，熱電偶測量接點熱物性參數(shù)和幾何參數(shù)本身也在隨溫度變化而變化，因此參數(shù)a1和a2并不是定值，這使得分數(shù)階模型在這個時間區(qū)間內(nèi)與測量值之間仍存在較大誤差。

4.4 變階次分數(shù)階傳遞函數(shù)參數(shù)和階次的辨識

(26)

式中：參數(shù)a1和μ1組成項即0.86s2.13項，以及參數(shù)a2和μ2組成項即2.28s0.92項，組合表示了測量接點熱慣性和熱電偶偶絲熱擴散所引起的延遲。輸出測量值和變階次模型估計值的比較如圖9所示。

從圖9可以看出，變階次分數(shù)階模型在時間軸0～4 s區(qū)間與實驗數(shù)據(jù)能較好地吻合。

4.5 固定階次模型和變階次模型擬合結(jié)果比較

圖10是輸出測量值、固定階次模型估計值與變階次模型估計值的比較，其中固定階次模型的評估函數(shù)J為0.043，變階次模型的評估函數(shù)J為0.020.

從圖10可以看到，變階次分數(shù)階模型克服了固定階次分數(shù)階模型在時間軸0～4 s區(qū)間與實驗數(shù)據(jù)不能很好吻合的缺點，又能反映熱電偶整體溫度響應(yīng)的變化特征。由于分數(shù)階算子可視為傳統(tǒng)算子微積分階數(shù)在相鄰整數(shù)值之間的插值，變階次分數(shù)階模型表達了分數(shù)階微分算子的階數(shù)插值性質(zhì)。分數(shù)階微積分運算可以被視為冪函數(shù)和算子作用函數(shù)的卷積，因而更好地描述了熱電偶溫度分布特征的全局相關(guān)性，以及溫度分布的記憶和延遲性質(zhì)。

5 結(jié)論

本文采用分數(shù)階微積分研究了熱電偶的動態(tài)特性。對于熱電偶測量固體表面瞬態(tài)溫度的問題，利用分數(shù)階建模和實驗建模相結(jié)合的方法，建立了一種熱電偶分數(shù)階傳遞函數(shù)的結(jié)構(gòu)，給出了參數(shù)和階次的估計方法，在300～400 K溫度區(qū)間給出了一個鎳鉻- 鎳硅熱電偶的算例。得到如下結(jié)論：

1)分數(shù)階傳遞函數(shù)能夠表征熱電偶測量接點熱慣性延遲和熱電偶偶絲熱擴散單相延遲。

2)固定階次分數(shù)階模型中s0.5表征了熱電偶偶絲熱擴散引起的延遲只是延遲的一種弱化表示。

3)變階次分數(shù)階模型能夠通過階次的插值，用兩種分數(shù)階次組合表示測量接點熱慣性和熱電偶偶絲熱擴散所引起的延遲。