李文軍, 王瀟楠, 鄭永軍
(中國計(jì)量大學(xué) 計(jì)量測試工程學(xué)院, 浙江 杭州 310018)
熱電偶被廣泛應(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)
式中: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)的拉普拉斯變換。
熱電偶由一對有共同接點(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)系式。
當(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)。
對于(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ù)
熱電偶測量接點(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ū)間。
為了獲得熱電偶在溫度激勵(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ì)下的輸出。
實(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ù)集。為了簡便,以下討論中測量值的平均值仍用測量值指代。
分?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.
為了驗(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í)值,表明方法有效。
圖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)與測量值之間仍存在較大誤差。
(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ù)能較好地吻合。
圖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ì)。
本文采用分?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ò)散所引起的延遲。