蒸汽供熱管網(wǎng)動(dòng)態(tài)仿真模型的研究

蒸汽供熱管網(wǎng)動(dòng)態(tài)仿真模型的研究

目前，城市集中已進(jìn)入自動(dòng)控制和經(jīng)濟(jì)發(fā)展階段。建立供熱管網(wǎng)的在線優(yōu)化運(yùn)行系統(tǒng),故障實(shí)時(shí)診斷系統(tǒng),以及仿真分析和培訓(xùn)系統(tǒng),是保證熱網(wǎng)安全、經(jīng)濟(jì)運(yùn)行的有效手段。統(tǒng)一的、高精度的動(dòng)態(tài)仿真模型是這些系統(tǒng)的重要基礎(chǔ)。模型的精度取決于模型的機(jī)理性、簡(jiǎn)化假設(shè)、特性參數(shù)如管道的阻力的準(zhǔn)確性度。實(shí)際的管網(wǎng)與設(shè)計(jì)一般有很多差別,同時(shí)隨著年代的變化,管網(wǎng)的不斷改造,管網(wǎng)的特性也會(huì)變化。因此,這些管網(wǎng)特性參數(shù)需要不斷的用實(shí)際管網(wǎng)的運(yùn)行數(shù)據(jù)進(jìn)行修正。由于蒸汽的可壓縮性、熱容小、物理性質(zhì)變化大,蒸汽管網(wǎng)的模型比熱水管網(wǎng)的模型要復(fù)雜。本文介紹了一種蒸汽供熱管網(wǎng)的動(dòng)態(tài)仿真模型。該模型將蒸汽處理成為單相可壓縮流體,并采用隱式Euler算法和稀疏矩陣算法求解。證明了所用的算法具有大范圍收斂特性,同時(shí)對(duì)模型的合理性做了分析。1其他簡(jiǎn)化假設(shè)蒸汽供熱一般為過(guò)熱蒸汽,并假設(shè)輸送過(guò)程中沒(méi)有相變,因而可認(rèn)為是單相可壓縮流體。其他簡(jiǎn)化假設(shè)如下:1)網(wǎng)絡(luò)由節(jié)點(diǎn)和支路組成,節(jié)點(diǎn)分為內(nèi)部節(jié)點(diǎn)和外部節(jié)點(diǎn)。內(nèi)部節(jié)點(diǎn)是管網(wǎng)的分叉點(diǎn);外部節(jié)點(diǎn)為熱力站、供熱源。2)節(jié)點(diǎn)內(nèi)的工質(zhì)處于均勻狀態(tài)。3)同一支路的截面積沿管長(zhǎng)不變,工質(zhì)參數(shù)用進(jìn)出口參數(shù)的加權(quán)平均來(lái)表示。基本方程1節(jié)點(diǎn)tVidρidt=Ν∑j=1DijGij.(1)Vidρidt=∑j=1NDijGij.(1)在一般情況下,焓H的變化比壓力p的變化要緩慢,由dρdt=οροpdpdt+οροΗdΗdt≈οροpdpdtdρdt=οροpdpdt+οροHdHdt≈οροpdpdt,得到下式:Cidpidt=Ν∑j=1DijGij?(2)其中:Ci=Viορiοpi,表示工質(zhì)的可壓縮能力。式(1),(2)中的其他符號(hào)說(shuō)明如下,ρ為工質(zhì)密度;t為時(shí)間;N為節(jié)點(diǎn)總數(shù);Gij為節(jié)點(diǎn)i,j之間管路的質(zhì)量流量;Vi為節(jié)點(diǎn)i的容積(包括管道的容積);Dij值為1時(shí),表示節(jié)點(diǎn)i,j之間有管道連接,管道方向從節(jié)點(diǎn)j指向i節(jié)點(diǎn);Dij值為-1時(shí),管道方向從節(jié)點(diǎn)i指向節(jié)點(diǎn)j;Dij值為0時(shí),表示節(jié)點(diǎn)i,j之間無(wú)聯(lián)系。另規(guī)定i=j時(shí),Dij的值為零。2管段蒸汽壓力和阻力ρijLijdUijdt=Dij(pj-pi+DijΗ′ij)-hw?(3)式中Η′=0.5ρ(U2in-U2out)+ρg(Ζin-Ζout)+hp表示宏觀動(dòng)能、勢(shì)能和管路中的泵產(chǎn)生的壓力,其中:Uin和Uout分別是管段入口和出口流速;Zin和Zout分別代表管段入口和出口的相對(duì)高度。由于管路的截面積沿管長(zhǎng)不變,蒸汽密度又小,因此宏觀動(dòng)能和重力勢(shì)能的影響可忽略不計(jì)。蒸汽網(wǎng)中沒(méi)有泵,hp=0,因此H′=0。hw=∑(0.5λLd+ξ)ρU2表示沿程阻力和局部阻力損失,其中:L為管段長(zhǎng)度,d為管段直徑,ξ為管段當(dāng)量阻力系數(shù),U為工質(zhì)流動(dòng)速度。令:管路的慣性系數(shù)Ιij=LijAij;管路的摩阻系數(shù)Rij=hw(AijUijρij)2,將Iij,Rij,Gij代入式(3)得ΙijdGijdt=Dij(pj-pi)-RijG2ij.(4)3節(jié)點(diǎn)間網(wǎng)絡(luò)聯(lián)系Vidρiuidt=Ν∑j=1DinijGijhj-Ν∑j=1DoutijGijhi+Qi,(5)其中:Dinij是內(nèi)部節(jié)點(diǎn)i和j間的只入網(wǎng)絡(luò)聯(lián)系;Doutij是節(jié)點(diǎn)i和j間的只出網(wǎng)絡(luò)聯(lián)系;定義如下:Dinij為1時(shí)表示有流體自節(jié)點(diǎn)j流入節(jié)點(diǎn)i;Dinij為0時(shí)表示兩節(jié)點(diǎn)間無(wú)網(wǎng)絡(luò)聯(lián)系或i=j;Doutij為1時(shí)表示有流體自節(jié)點(diǎn)i流入節(jié)點(diǎn)j;為0時(shí)表示兩節(jié)點(diǎn)間無(wú)網(wǎng)絡(luò)聯(lián)系或i=j;D=Din-Dout.由于ρu=ρ(h-pv)=ρh-p,因此Vdρudt=Vhdρdt+Vρdhdt-Vdpdt。將此式和式(1)代入式(5)得Μidhidt=Ν∑j=1[Dinij(hj-hi)Gij]+Qi+Vidpidt?(6)其中:M=Vρ,Q=γ[Tenvir-0.5(Tin+Tout)];γ是管段上的散熱系數(shù),單位為W/K。4)Gi,Ν=Ci(pi-pa)Τtci表示內(nèi)部節(jié)點(diǎn)與環(huán)境(被定義為第N個(gè)節(jié)點(diǎn))之間的物質(zhì)交換。Ttci為節(jié)i點(diǎn)對(duì)環(huán)境的泄漏常數(shù);pa為環(huán)境壓力。2計(jì)算方法方程(2),(4),(6)構(gòu)成的模型是非線性的代數(shù)/微分方程組。對(duì)于動(dòng)態(tài)仿真,仿真計(jì)算方法應(yīng)既保證計(jì)算的穩(wěn)定性又保證實(shí)時(shí)性。2.1穩(wěn)定的特性分析式(4)中I體現(xiàn)為對(duì)流速的遲滯作用,通?？梢院雎圆挥?jì)。因此式(2)可寫(xiě)作Cidpidt=Ν∑j=1Dij√Dij(pj-pi)Rij+Ci(pi-pa)Τtci.(7)式(6)中的Vdpdt一項(xiàng)表示的是推進(jìn)功對(duì)工質(zhì)溫度的影響,在小擾動(dòng)的情況下,通常比傳熱項(xiàng)的影響要小,因而也可以忽略此項(xiàng)的作用。式(6)也可改寫(xiě)為下面形式:Μidhidt=Ν∑j=1DinijGij(hj-hi)+Qi.(8)為保證計(jì)算的實(shí)時(shí)性和穩(wěn)定性,本文采用具有A穩(wěn)定特性的單步隱式Euler方法。由式(7)和(8)可得:Μi+1iht+1i-htiΔt=Ν∑j=1DinijGtij(ht+1j-ht+1i)+γij[0.5(Τt+1i+Τt+1j)-Τa]?(9)Ct+1ipt+1i-ptiΔt=Ν∑j=1Dij√Dij(pt+1j-pt+1i)Rij+Ct+1i(pt+1i-pa)Τtci.(10)為減少運(yùn)算量,在小擾動(dòng)條件下,Ct+1i用Cti近似,Tt+1i用Tti近似,Mt+1i用Mti近似不會(huì)產(chǎn)生太大的誤差。由此,式(9)和式(10)已經(jīng)完全解耦。在計(jì)算過(guò)程中,先求解由式(10)構(gòu)成的N階非線性代數(shù)方程組,得到各節(jié)點(diǎn)的壓力后,再求出各管路中的流量,最后代入由式(9)構(gòu)成的非線性代數(shù)方程組,得到各節(jié)點(diǎn)的溫度。2.2逆矩陣a-1由式(4)得到,Gt+1ij=Dij(pt+1j-pt+1i)Rij|Gt+1ij|,Ct+1ipt+1i-ptiΔt=Ν∑j=1pt+1j-pt+1iRij|Gt+1ij|+Ct+11pt+1i-paΤtci?可以簡(jiǎn)寫(xiě)為At+1pt+1=Bt+1.(11)其中:pt+1=[pt+11pt+12?pt+1Ν]?At+1=[at+1ij]?Bt+1=[bt+11bt+12?bt+1Ν]?at+1ij=-D2ijRij|Gt+1ij|,i≠j；at+1ii=Ct+1i(1Δt+1Ctci)+Ν∑j=11Rij|Gt+1ij|?i=1?2???Ν；bt+1i=Ν∑j=1pt+1jRij|G(k)ij|+Ct+1i(ptiΔt+paΤtci).建立如下迭代關(guān)系:Ρt+1k+1=(At+1k)-1Bt+1k?(12)Gt+1ij,k=sign[Dij(pt+1j,k-pt+1i,k)]?√Dij(pt+1j,k-pt+1i,k)Rij?其中:at+1ij,k=-D2ijRij|Gt+1ij,k|,i≠j；at+1ii,k=Ct+1i(1Δt+1Ctci)+Ν∑j=11Rij|Gt+1ij,k|,i=1?2???Ν；bt+1i,k=Ν∑j=1pt+1jRij|Gkij,k|+Ct+1i(ptiΔt+paΤtci).下面證明式(12)的迭代公式是大范圍收斂的。引理1存在正定矩陣A和半正定矩陣B,則C=AB仍是半正定的。引理2對(duì)于映射F:D?Rn→Rn于凸開(kāi)集上D0?D上G可導(dǎo),則F于D0上單調(diào)的充分必要條件是對(duì)所有的x∈D0,F′(x)半正定。引理3大范圍收斂的充分條件:若映射F:D?Rn→Rn在D上是保序的,x0≤y0,〈x0,y0〉?D,則序列x(k+1)=F(x(k)),y(k+1)=F(y(k))能滿足x(k)↑x*,y(k)↓y*與x*≤y*;若映射在〈x0,y0〉上連續(xù),x*=F(x*),y*=F(y*),并且映射的任何不動(dòng)點(diǎn)均在〈x0,y0〉內(nèi)。引理4Gerschgorin圓盤定理:n階矩陣A=[aij]的每個(gè)特征值都是復(fù)平面諸圓盤di={λ||λ-aii|≤n∑j=1,j≠i|aij|?i=1???n}中的一個(gè)。引理5實(shí)對(duì)稱矩陣的特征值全部為實(shí)數(shù)。引理6設(shè)映射G(x)=H(x)+C,H(x)是保序的,且H(0)=0。則設(shè)輔助迭代:W(0)=0,W(k+1)=H(W(k))+d其中d≥0,d-c>0,則{W(k)}單增收斂。從而有正整數(shù)N使得0≤W(N+1)-W(N)≤d-c成立。并且,若c≥0則W(1)-W(0)≥d-c。設(shè)有非負(fù)整數(shù)L(≤N)使得下式成立W(L+1)-W(L)≥d-c,則可取x(0)=W(L),y(0)=W(N)作為初始向量。(x(0)=0肯定是可以的。)證明1A是正定矩陣?？疾霢,由[引理4]可知:矩陣A的特征值的實(shí)部均大于Ci(t)(1Δt+1Τtci)。由于Ci(t)恒大于零,因此,矩陣A的特征值的實(shí)部均大于0。另外,很顯然矩陣A為實(shí)對(duì)稱矩陣,由[引理5]得:A的全部特征值均為實(shí)數(shù)。因而A為正定矩陣,且存在正定的逆矩陣A-1。證明2映射F:F(pt+1k)=[A(pt+1k)]-1[B(pt+1k)]是單調(diào)的。令:F′(pt+1k)=A-1Q,考察[qij]:qii=Ν∑l=112RilGt+1il,k(Gt+1il,k+1Gt+1il,k)2?qij=-12RjlGt+1jl,k(Gt+1jl,k+1Gt+1jl,k)2.顯然qij=qji,因此由推論5可知:Q的特征值均為實(shí)數(shù)。易知qii≥Ν∑j=1,j≠i|qij|?i=1?2???Ν,由引理4,矩陣Q的特征值的實(shí)部均不小于0。因此矩陣Q是一個(gè)正半定的矩陣。由引理1可知F′(pt+1k)是正半定的。由[引理2]可得映射是F單調(diào)的。證明3迭代公式(12)取上一時(shí)刻的值作初值時(shí)迭代收斂。由證明2和推論1可知:式(12)定義的映射是保序的。由式(12)的物理意義可知F(0)≥0(否則就會(huì)出現(xiàn)供熱源和用戶站的壓力均為正壓而內(nèi)部節(jié)點(diǎn)卻為負(fù)壓的狀況)。映射可改寫(xiě)為如下形式:F(pt+1k)=H(pt+1k)+C。其中H(pt+1k)滿足H(0)=0的條件;C≥0。容易證得:H(pt+1k)依然保持保序性質(zhì)。對(duì)于初值為上一時(shí)刻的值p*(t),必定存在一個(gè)d滿足引理6中輔助迭代的條件,使得值p*(t)可以成為初始值。由引理3可知,以t時(shí)刻的值作為式(12)的迭代初值,迭代過(guò)程是收斂的。3仿真結(jié)果的計(jì)算本文給出的蒸汽管網(wǎng)的單相可壓縮流體網(wǎng)絡(luò)模型,及實(shí)時(shí)仿真的數(shù)值計(jì)算方法,在實(shí)際應(yīng)用時(shí),需要考慮以下問(wèn)題。因?yàn)楸灸Ｐ褪羌袇?shù)模型,管路的純延遲特性在模型中沒(méi)有體現(xiàn)。通過(guò)增加節(jié)點(diǎn)可以模擬管道的純延遲現(xiàn)象,但這會(huì)增加計(jì)算量,因此應(yīng)根據(jù)管長(zhǎng)適當(dāng)增加節(jié)點(diǎn)數(shù),直至滿足仿真精度。在使用迭代公式(12)時(shí),每次迭代都需要求解線性方程組。普通的消元法求解線性方程組的存儲(chǔ)量為O(N2),計(jì)算量為O(N3),當(dāng)對(duì)象的節(jié)點(diǎn)數(shù)較大時(shí),資源和時(shí)間消耗過(guò)多。實(shí)際上,公式(12)中線性方程組的系數(shù)矩陣是對(duì)角占優(yōu)的實(shí)