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

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

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