熱電偶分?jǐn)?shù)階傳遞函數(shù)階次和參數(shù)估計(jì)

李文軍， 王瀟楠， 鄭永軍

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

0 引言

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

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

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

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

1 熱電偶時(shí)間分?jǐn)?shù)階熱傳導(dǎo)模型

1.1 分?jǐn)?shù)階微積分定義和分?jǐn)?shù)階拉普拉斯變換

(1)

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

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

(2)

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

(3)

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

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

熱電偶由一對有共同接點(diǎn)的材料互異的熱電偶偶絲組成，如圖1所示。圖1中，熱電偶偶絲的直徑為d，熱電偶偶絲1與熱電偶偶絲2之間存在一個(gè)夾角[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)式做分?jǐn)?shù)階拉普拉斯反變換，得到時(shí)間分?jǐn)?shù)階微分形式的關(guān)系式為

(9)

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

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

(10)

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

(11)

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

(12)

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

(13)

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

(14)

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

1.3 露端式熱電偶測量接點(diǎn)的傳熱方程

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

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

(15)

(16)

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

(17)

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

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

(18)

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

(19)

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

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

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

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

(20)

(21)

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

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

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

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

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

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

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

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

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

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

3.1 實(shí)驗(yàn)系統(tǒng)

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

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

3.2 實(shí)驗(yàn)數(shù)據(jù)

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

表6 實(shí)驗(yàn)配置

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

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

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

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

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

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

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

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

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

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

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

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

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

(22)

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

4.2 改進(jìn)Luss-Jaakola算法的驗(yàn)證

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

(23)

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

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

β=[a1a2a3μ1μ2],

向量真實(shí)值為

β=[1.20 1.30 1.00 2.10 0.50],

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

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

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

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

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

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

(24)

(25)

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

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

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

(26)

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

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

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

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

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

5 結(jié)論

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

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

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

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