運(yùn)用MATLAB實現(xiàn)換熱網(wǎng)絡(luò)動態(tài)特性的通用自動求解

姜 濤,崔國民

(上海理工大學(xué)熱工程研究所,上海 200093)

0 引 言

通常,一個換熱網(wǎng)絡(luò)是一個化工系統(tǒng)的子系統(tǒng),所以換熱網(wǎng)絡(luò)物流的供給性質(zhì)例如輸入溫度和流率常常受化工過程中其他因素的影響而偏離它們的標(biāo)準(zhǔn)值。如果標(biāo)準(zhǔn)值要求保持不變或在一個給定的范圍之內(nèi),當(dāng)供給性質(zhì)變化時,就需要一個控制系統(tǒng)[1]。換熱設(shè)備的動態(tài)特性一直以來是一個研究課題,它是分析、設(shè)計和改進(jìn)換熱設(shè)備及其控制系統(tǒng)的重要依據(jù),因此受到眾多學(xué)者的關(guān)注,其中換熱設(shè)備的傳遞特性顯得尤為重要[1～3]。但是到目前為止,研究人員主要把目光集中于單體換熱器的動態(tài)特性研究,這方面的成果層出不窮,推導(dǎo)得到了各種形式的傳遞函數(shù)。而對于換熱網(wǎng)絡(luò)動態(tài)特性的研究,目前主要針對于某些特定的熱力系統(tǒng)或結(jié)構(gòu),而且鮮見關(guān)于換熱網(wǎng)絡(luò)傳遞特性研究方面的資料。李政[4]等在國產(chǎn)200 t/h循環(huán)流化床鍋爐建模與仿真工作中開發(fā)了一種復(fù)雜換熱器系統(tǒng)動態(tài)特性的計算方法;楊德明[5]以能量守恒為基礎(chǔ),建立了換熱網(wǎng)絡(luò)的機(jī)理動態(tài)數(shù)學(xué)模型。目前這些數(shù)學(xué)模型及表達(dá)式的通用性和實用性都較差,尚缺乏一種能普遍廣泛運(yùn)用的通用求解方法,因此換熱網(wǎng)絡(luò)動態(tài)特性的研究具有十分重要的意義。

本文基于信號流圖,提出一種快速求解換熱網(wǎng)絡(luò)動態(tài)特性的通用方法,該方法以單體換熱器為模塊,連接各個換熱器形成換熱網(wǎng)絡(luò)的信號流圖,以換熱器入口信號是否準(zhǔn)備好作為模塊計算依據(jù),采用以模塊為核心的動態(tài)求解策略依次求解各個換熱器,直到所有換熱器計算完畢,最終獲得換熱網(wǎng)絡(luò)流體出口溫度變化對入口參數(shù)變化的傳遞函數(shù)。同時本文運(yùn)用MATLAB編制了自動求解程序,MATLAB強(qiáng)大的符號運(yùn)算功能使得計算方便簡潔,速度快,且通用性強(qiáng)。

1 單體換熱器的處理

在工業(yè)過程中,大部分被控對象為具有純滯后的一階或二階慣性環(huán)節(jié)[6]。根據(jù)換熱器的特點,其動態(tài)數(shù)學(xué)模型可以用具有純滯后的一階慣性環(huán)節(jié)表示,其傳遞函數(shù)表示為

式中:W(s)為出口參數(shù)變化與入口參數(shù)變化的比值;U(s)和Y(s)分別為換熱器輸入和輸出參數(shù)變化量的拉氏變換;T,K,τ為3個決定換熱器動態(tài)響應(yīng)過程的動態(tài)特性參數(shù),其中,T為時間常數(shù),反映流體參數(shù)變化的快慢,K為靜態(tài)放大系數(shù),τ為滯后時間。這3個動態(tài)特性參數(shù)可以通過理論分析計算或曲線擬合得出,對于不同的入口參數(shù)變化,動態(tài)特性參數(shù)也不同,因此有24個不同參數(shù)。當(dāng)入口參數(shù)變化不大時,出口參數(shù)的變化可看成是各入口參數(shù)引起的變化的線性疊加,這樣單體換熱器冷熱流體出口溫度的變化可寫成:

此時換熱器即相當(dāng)于一個4輸入 4輸出系統(tǒng),如圖1(a)所示,將其轉(zhuǎn)換成信號流圖形式如圖1(b)所示。

圖1 單體換熱器傳遞函數(shù)模型及其信號流圖Fig.1 Transfer functions and signal flow graph of heat exchanger

在MAT LAB中創(chuàng)建一個二維數(shù)組DANAMICPA存放這些動態(tài)特性參數(shù)的具體數(shù)值,該數(shù)組分24列,行號代表換熱器號,以供計算所需。

2 換熱網(wǎng)絡(luò)建構(gòu)及識別

如上所述,每個單體換熱器都可轉(zhuǎn)換成信號流圖形式,這樣各個換熱器連接成的換熱網(wǎng)絡(luò)也就可以轉(zhuǎn)換成信號流圖的形式。在換熱網(wǎng)絡(luò)中,同股流體流經(jīng)的各個換熱器,流體的入口參數(shù)變化為上一個換熱器的出口參數(shù)變化,也就是說,上一個換熱器的出口信號變化輸入給了下一個換熱器,依據(jù)這種傳遞關(guān)系,一個2×2一級換熱網(wǎng)絡(luò)可以用信號流圖表示成如圖2所示的結(jié)構(gòu)。

圖2 換熱網(wǎng)絡(luò)超結(jié)構(gòu)模型及其信號流圖Fig.2 Signalflow graph of a 2×2 heat exchanger network

如圖所示,虛框內(nèi)為一個換熱器的信號圖,即一個模塊,圓點表示信號節(jié)點,實線表示節(jié)點間的傳遞函數(shù)。各個模塊間用虛線連接,表示連接管道,此上不考慮熱量損失,只考慮時間滯后。

本文基于單體換熱器的信號圖,采用部件模塊連接的方式構(gòu)建換熱網(wǎng)絡(luò)的信號流圖。將單體換熱器看成一個小模塊,模塊入口信號為流入該換熱器的流體信號變化的s域表達(dá)式,該表達(dá)式乘以傳遞函數(shù)即為出口信號變化的s域表達(dá)式,又作為該股流體流向的下一個換熱器的入口信號。這樣各個換熱器形成一個個小模塊,按照流體間的實際換熱匹配和流體先后順序,將各個單體換熱器的出、入口信號通過連線 (管道)連接起來,形成換熱網(wǎng)絡(luò)的信號流圖。這種方法具有直觀、方便以及通用性強(qiáng)的特點。

網(wǎng)絡(luò)識別采用索引矩陣法,矩陣行號表示換熱器號,即單元模塊號。矩陣分4列,第1列表示該行號的換熱器的1熱流體從哪個單元模塊輸出,第2列表示該行號的換熱器的熱流體輸入到哪個單元模塊去,第3列表示該行號的換熱器的冷流體從哪個單元模塊輸出,第4列表示該行號的換熱器的冷流體輸入到哪個單元模塊去,0表示給流體為換熱網(wǎng)絡(luò)原始輸入流體或最終輸出流體。這樣,換熱網(wǎng)絡(luò)的結(jié)構(gòu)就能自動識別出來。如圖2所示結(jié)構(gòu)的換熱網(wǎng)絡(luò)的索引矩陣的代碼為

3 自動求解策略

換熱網(wǎng)絡(luò)動態(tài)特性的求解建立在獲得換熱網(wǎng)絡(luò)傳遞函數(shù)基礎(chǔ)上,因此問題關(guān)鍵在于如何獲得換熱網(wǎng)絡(luò)的傳遞函數(shù)。

本文采用的求解策略不采用冷流體或熱流體作為模擬路徑,沿流體流動方向循序求解的方法,而是在求解過程中,采用以模塊即換熱器為核心的動態(tài)求解策略,即在求解過程中,以換熱器入口信號是否準(zhǔn)備好作為求解的依據(jù),所謂的準(zhǔn)備好就是指該入口信號是否為已知,隨時計算可以求解的換熱器。這種求解策略首先從入口模塊——熱源和冷源變化信號都已知的換熱器入手,對其進(jìn)行求解,然后將求解得到的出口信號傳遞到下面的換熱器,為下面的換熱器求解做準(zhǔn)備。繼續(xù)對各個換熱器進(jìn)行搜索,如果有入口信號都準(zhǔn)備好而且沒有計算過的換熱器,對其進(jìn)行求解,直到所有的換熱器都求解完畢,完成模擬。

首先創(chuàng)建二維矩陣EACHDECIDE判斷單體換熱器各入口信號是否準(zhǔn)備好,該矩陣行號表示換熱器號,兩列分別表示熱流體和冷流體是否準(zhǔn)備好,矩陣元素為1表示準(zhǔn)備好,矩陣元素為0表示沒有準(zhǔn)備好。若該換熱器的流體為源流體,則元素為1,因為換熱網(wǎng)絡(luò)流體入口參數(shù)變化是已知的。

然后創(chuàng)建二維矩陣DECIDE判斷換熱器是否能計算,同樣行號表示換熱器號,第1列表示該行號換熱器冷、熱流體的入口信號是否已知,只有矩陣EACHDECIDE中該行元素都為 1矩陣DECIDE該行才為1,表示換熱器準(zhǔn)備好可以進(jìn)行計算;第2列表示該行號換熱器之前是否計算過,元素為1表示已經(jīng)計算過,元素為0表示尚未計算過。

再次給出每個換熱器的輸入和輸出信號矩陣EACHINPUT和 EACHOUTPUT,這兩個矩陣為符號矩陣,矩陣元素為有符號 s的表達(dá)式,即為換熱器入口或出口信號的s域表達(dá)式。矩陣有4列,分別表示熱流體溫度、冷流體溫度、熱流體流量和冷流體流量的入口或出口信號表達(dá)式。

圖3給出了程序計算的流程圖,這樣可以得到所有換熱器流體出口溫度對換熱網(wǎng)絡(luò)流體入口參數(shù)變化的傳遞函數(shù)。如果索引矩陣中第3列或者第4列為0,表示該流體為換熱網(wǎng)絡(luò)最終輸出流體,其傳遞函數(shù)即為換熱網(wǎng)絡(luò)該流體的傳遞函數(shù)。得到換熱網(wǎng)絡(luò)的傳遞函數(shù)后,在MATLAB中經(jīng)拉式反變換,得到時域內(nèi)流體溫度關(guān)于時間的表達(dá)式,帶入數(shù)值計算即可觀察換熱網(wǎng)絡(luò)的動態(tài)特性。

圖3 程序計算流程圖Fig.3 Calculation flowchart of the program

4 計算實例

本文以化工廠4熱、5冷的換熱網(wǎng)絡(luò)為例,網(wǎng)絡(luò)結(jié)構(gòu)如圖4所示,物流數(shù)據(jù)如表1所示。

各動態(tài)特性參數(shù)求解采用以下公式。

流體出口溫度對入口溫度變化放大系數(shù)。

圖4 換熱網(wǎng)絡(luò)結(jié)構(gòu)Fig.4 Structure chart of the heat exchanger network

在上面的表達(dá)式中,K的下標(biāo)cth表示冷流體出口溫度變化對熱流體入口溫度變化的放大系數(shù),其它的以此類推。

時間常數(shù):

滯后時間τ:

(1)對流量變化的響應(yīng)

對于不可壓縮流體,當(dāng)換熱器通道內(nèi)流體的流量發(fā)生變化時,試件在整個流程上的換熱狀況都有改變,因此無論是逆流還是順流工況,在出口處的溫度響應(yīng)不存在延滯。

(2)對于溫度變化的響應(yīng)

a.在順流工況下,出口溫度的滯后時間取決于兩股流體中流速最快的流體,它最先將流體進(jìn)口溫度變化的影響帶到出口處。也就是說:τ=L/vmax。

b.在逆流工況下,進(jìn)口溫度發(fā)生變化的流體出口溫度響應(yīng)的滯后時間由 τ=L/va來確定,而未發(fā)生變化的流體響應(yīng)的滯后時間為零。

表1 換熱器及流體參數(shù)Tab.1 Parameter Table of heat exchangers and fluids

當(dāng)換熱網(wǎng)絡(luò)中熱流體1的入口溫度發(fā)生10℃的上升階躍時,利用MATLAB編制的自動求解程序,獲得換熱網(wǎng)絡(luò)流體出口溫度變化對入口溫度變化的傳遞函數(shù),并拉式反變換后得到時域內(nèi)流體變化表達(dá)式如下,同時變化曲線如圖5所示。

從圖5中可以看出,流體C4,H1,C2,C5,H3出口溫度變化響應(yīng)均存在一定滯后,且滯后時間不同,滯后時間長短與流體的流程有關(guān)。由于流體為逆流布置,所以h1入口溫度發(fā)生變化時,從圖4可以看出,流體C3出口溫度立即變化,沒有滯后。上述各股流體均經(jīng)過一段時間后達(dá)到一個新的穩(wěn)態(tài)值,曲線形狀符合實際情況,說明本文提出的方法是可行的。

5 結(jié) 論

本文提出了一種換熱網(wǎng)絡(luò)動態(tài)特性的通用求解方法,該方法以單體換熱器為模塊,連接各個換熱器形成換熱網(wǎng)絡(luò)的信號流圖,以換熱器入口信號是否準(zhǔn)備好作為模塊計算依據(jù),采用以模塊為核心的動態(tài)求解策略依次求解各個換熱器,直到所有換熱器計算完畢,最終獲得換熱網(wǎng)絡(luò)流體出口溫度變化對入口參數(shù)變化的傳遞函數(shù)。同時本文運(yùn)用MAT LAB編制了自動求解程序,MATLAB強(qiáng)大的符號運(yùn)算功能使得計算方便簡潔,速度快,且通用性強(qiáng)。這種求解方法為換熱網(wǎng)絡(luò)的動態(tài)特性研究提供了一種新的思路,同時為換熱網(wǎng)絡(luò)的控制奠定的基礎(chǔ)。

圖5 熱流體、冷流體出口溫度變化Fig.5 Output temperature change of hot and clod fluids

[1]高維平,楊瑩.換熱網(wǎng)絡(luò)的優(yōu)化改進(jìn)和控制[J].化工學(xué)報,2003,54(7):965-971.

[2]崔國民,張勤,陸貞,等.序貫?zāi)K法實現(xiàn)換熱網(wǎng)絡(luò)的模擬[J].工程熱物理學(xué)報,2006,27(4):682-684.

[3]郭佳,崔國民,王方方,等.換熱器動態(tài)設(shè)計方法研究.工程熱物理學(xué)報,2007,28(5):865-867.

[4]李政,倪維斗.復(fù)雜換熱器系統(tǒng)的動態(tài)特性計算.熱能動力工程,1997,12(4):292-296.

[5]楊德明.換熱網(wǎng)絡(luò)動態(tài)數(shù)學(xué)模型及其求解.江蘇工業(yè)學(xué)院學(xué)報,2006,18(1):10-13.

[6]Roetzel W.,Xuan Y.Dynamic Behaviour of Heat Exchangers,Boston:WIT Press,1999.