蒸汽管道動態水力熱力特性研究

肖木森, 張 沖, 郭磊宏, 齊智猛, 范 巍, 章雨然

(1.天津市津安熱電有限公司,天津 300000;2.天津大學 環境科學與工程學院,天津 300350)

1 概述

蒸汽熱網是集中供熱系統中連接熱源和用戶的紐帶,關系到生產運營安全和系統能源利用效率。蒸汽作為載熱介質,可同時滿足不同壓力和溫度要求的多種用戶的用熱需求,適用性較廣。蒸汽密度小,因此蒸汽系統適用于高層建筑的供暖,不容易使底部設備超壓。蒸汽熱惰性小,供汽時熱得快,停汽時冷得也快,因此在工業生產的各個工藝環節中,蒸汽系統也得到了廣泛的應用[1-2]。

相比于熱水供熱系統,蒸汽供熱系統內蒸汽流動狀態復雜多變,且出現故障危害大,了解蒸汽在管道中流動的水力狀態和熱力狀態不僅能夠促進蒸汽管網調控向智能化、精細化轉型,也能夠提高整個供熱系統的安全可靠性。因此有必要對蒸汽熱網建立準確且高效的水力熱力耦合數學模型,并能在負荷頻繁波動和室外溫度周期性波動運行工況下計算出供熱管網中各管段、各末端蒸汽的水力狀態和熱力狀態[3-5]。

國內外學者針對蒸汽管網的水力模型和熱力模型開展了一系列研究,史琳等人[6]考慮蒸汽為可壓縮理想氣體,建立了蒸汽管網的水力熱力耦合模型,可以計算流量分配以及溫度、壓力分布,并將該模型應用到大型蒸汽管網中。李世武[7]提出了蒸汽管網的模擬計算方法,該方法基于節點模擬法,以節點方程為基礎,提出了蒸汽管網計算的節點殘量修正計算方法。Kiuchi[8]采用全隱式有限差分法計算管網中的非穩態氣體流動,并基于牛頓-拉夫遜法求解管道有限差分方程。Bermudez等人[9]介紹了管道內實際氣體非等溫、非絕熱可壓縮流動數學模型數值求解的有限體積格式,采用了Van Leer的Q-Scheme數值格式對數學模型中的通量進行迎風離散,與標準歐拉方程不同,該模型考慮了壁面摩擦、變高度以及管道與環境之間的傳熱。Marconcini等人[10]利用連續性方程、動量方程和能量守恒方程,編寫了非絕熱蒸汽管網中過熱蒸汽流動狀態參數的計算機程序VAPSTAT1,它通過計算流體在給定節點的流量、壓力和溫度來模擬蒸汽管網。

雖然在蒸汽管道建模及仿真方面已有較多研究,但大部分研究對蒸汽在管道流動過程中的狀態參數變化考慮不足,或將其忽略,或認為蒸汽是理想氣體并進行簡化計算,但蒸汽的狀態參數變化直接影響到蒸汽管道水力熱力模型計算的準確性。因此,本文綜合考慮蒸汽作為非理想氣體的狀態方程和焓方程,以質量守恒、動量守恒和能量守恒定律為基礎,建立蒸汽管道水力熱力耦合特性的非穩態模型,并模擬分析在環境溫度周期性變化、管道入口壓力變化、管道入口溫度變化這3種常見的非穩態工況下,管道內蒸汽參數的變化情況,得到其變化的一般性規律。

2 蒸汽管網動態水力熱力特性研究方法

2.1 蒸汽狀態參數方程

蒸汽為可壓縮流體,在管道內流動過程中,蒸汽壓力、溫度的不斷變化會導致其比體積、比焓的變化,因此蒸汽比體積、比焓的狀態方程與蒸汽在管道中流動的控制方程相互聯系。本文采用蒸汽的狀態方程如下[11]:

(1)

(2)

式中v--蒸汽比體積,m3/kg

p--蒸汽壓力,Pa

Rm--蒸汽氣體常數,J/(kg·K)

T--蒸汽溫度,K

Cp--蒸汽的無量綱壓力,為蒸汽壓力與參考壓力的比值,參考壓力取1 MPa

h--蒸汽比焓,J/kg

CT--蒸汽的無量綱溫度,為蒸汽溫度與參考溫度的比值,參考溫度取540 K

理想及剩余部分的無量綱壓力、溫度關聯式的計算方法依次如下:

(3)

(4)

(5)

(6)

計算系數的取值見表1,式(5)中計算系數的取值見表2。

表1 式(4)與式(6)中計算系數的取值[11]

續表

表2 式(5)中計算系數的取值[11]

蒸汽動力黏度隨蒸汽溫度和壓力變化而變化,本文采用的蒸汽動力黏度計算公式如下[12],此處蒸汽的無量綱溫度的參考溫度取值與前文不同,取647.096 K。

Cμ=Cμ0Cμ1Cμ2

(7)

(8)

(9)

式中Cμ--蒸汽無量綱動力黏度,為蒸汽動力黏度與參考動力黏度的比值,參考動力黏度取1.00×10-6Pa·s

Cμ0--蒸汽在稀釋氣體極限下的無量綱動力黏度

Cμ1--有限的密度對蒸汽無量綱動力黏度的影響

Cμ2--蒸汽無量綱動力黏度的臨界增強,臨界增強僅在臨界點附近的密度和溫度非常小的區域顯著,本文的蒸汽狀態區域與臨界點距離較遠,在此不考慮,取1

Hi、Hij--計算系數

Cρ--蒸汽無量綱密度,為蒸汽密度與參考密度的比值,參考密度取322.0 kg/m3

式(8)中計算系數Hi的取值見表3,式(9)中計算系數Hij的取值見表4。

表3 式(8)中計算系數Hi的取值[12]

表4 式(9)中計算系數Hij的取值[12]

續表

2.2 蒸汽管道動態模型

蒸汽在管道中的實際流動過程十分復雜,為了便于模型的建立與求解,進行了以下兩項假設:在實際管網中,管道長度遠大于管道外直徑,故假設蒸汽管道是一維的,即蒸汽參數沿管道直徑方向保持恒定;蒸汽通過管壁和保溫層向環境傳熱,傳熱系數沿管道變化不大,故將傳熱系數視為常數。對蒸汽管道進行控制體劃分及分析,蒸汽管道外存在保溫層,內部包含蒸汽主流,可壓縮流體(即蒸汽)通過垂直于流向的截面積為A的管道,其受到重力、管壁摩擦、流體壓力等作用,并時刻與外界環境進行熱交換,狀態參數不斷發生變化。

蒸汽在管道中流動滿足質量守恒,其連續性方程如下[13]:

(10)

式中ρ--蒸汽密度,kg/m3

t--時間,s

u--蒸汽速度,m/s

x--某個節點距坐標原點的距離,m

蒸汽在管道中為一維流動,且速度較大,忽略黏性切應力引起的動量交換,則蒸汽在管道中流動的動量守恒方程如下:

(11)

式中f--蒸汽在管道中流動的沿程阻力系數

dpi--蒸汽管道的內直徑,m

g--重力加速度,m/s2,取9.801 1 m/s2

θ--蒸汽管道與水平方向的夾角

沿程阻力系數f是雷諾數Re和管道內表面粗糙度的函數,其關系式是通過實驗或者半理論分析得到的。本文采用阿里特蘇里公式,其計算方便,且適用性較廣。阿里特蘇里公式如下所示:

(12)

式中Δ--管道內表面粗糙度,m

Re--雷諾數

蒸汽在管道中不斷與環境進行熱交換,對于可壓縮流體,蒸汽溫度的變化會導致蒸汽密度的變化,從而影響流場的分布,故需要計算蒸汽沿管道的溫度分布。忽略導熱產生的熱量交換,則蒸汽在管道中流動的能量守恒方程如下:

(13)

式中k--蒸汽與環境的傳熱系數,W/(m2·K)

Ta--環境溫度,K

對于架空管道,蒸汽在管道中流動的傳熱由5部分組成,包括蒸汽與管壁的對流傳熱、管壁的導熱、保溫層的導熱、保溫層外壁與空氣的對流傳熱以及保溫層外壁與空氣的輻射傳熱。在以管道公稱直徑計算傳熱面積的條件下,該過程的傳熱系數計算式如下:

(14)

式中dp--管道公稱直徑,m

hm--蒸汽與管壁的表面傳熱系數,W/(m2·K)

dpo--管道的外直徑,m

λp--管道熱導率,W/(m·K)

dio--保溫層的外直徑,m

λi--保溫層熱導率,W/(m·K)

ha--保溫層外壁與空氣的表面傳熱系數,W/(m2·K)

hr--保溫層外壁的輻射傳熱系數,W/(m2·K)

2.3 動態模型的離散與求解

2.3.1動態模型的離散

首先對動態模型進行離散。本文在上述模型計算中采用交錯網格技術,選取有限容積法中的控制容積積分法,對上述方程進行離散。考慮到蒸汽在管道中流速很快,通過對流傳遞的熱量比通過導熱傳遞的熱量大很多,采用一階迎風格式可以滿足計算精度。離散所采用的交錯網格見圖1。

圖1 交錯網格分布

主控制體以物性節點(空心圓點)為中心,速度節點(實心黑點)位于主控制體的邊界。相鄰兩個速度節點間的距離為一個控制體長度Δx。主控制體內部,儲存了該控制體的蒸汽壓力、蒸汽密度、蒸汽溫度(蒸汽比焓)等信息,而主控制體的邊界儲存速度,同一個主控制體的左右兩個邊界速度不一定相同。

從左邊起,第1個節點為i=1,儲存入口蒸汽速度,第2個節點為j=1,儲存入口蒸汽壓力與入口蒸汽溫度。第3個節點為i=2,第4個節點為j=2,以此類推,倒數第2個節點為j=n-1,最后1個節點為i=n。

采用分段線性假設計算密度及動量,并利用一階迎風格式對上述控制方程進行離散,過程如下。

① 動量方程

對t時刻節點i所在的控制體進行積分,方程如下,式中量的符號上標表示時間,下標表示節點。方程系數除外,其上標表示節點。

(15)

式中 ΔV--控制體體積,m3

Δt--離散時間單位,s

Δx--離散空間單位,m

速度、密度的計算采用一階迎風格式,按下式計算:

(16)

對控制體交界處的密度采用分段線性假設,能夠得到以下關系式:

(17)

(18)

將式(18)整理成下式形式:

(19)

方程系數的計算式為:

(20)

(21)

(22)

(23)

式(19)～(23)構成離散的動量方程。

(24)

(25)

當i=3,…,n-1時,方程系數按式(20)～(23)計算。

當i=n時,動量方程出口的邊界條件采用了局部單向化假設,即計算結果依賴于n-1節點,不依賴于外界環境,控制體長度變為正常的一半。此時,離散的動量方程的系數有兩項產生了改變,式(20)變為式(26),式(23)變為式(27):

(26)

(27)

② 連續性方程

對節點j所在的控制體進行積分,方程如下:

(28)

式(28)可化簡為:

(29)

將式(17)代入式(29)可得:

(30)

此處引入壓力與速度的修正方程。壓力修正方程如下:

p=p*+p′

(31)

式中p*--壓力計算項,Pa

p′--壓力修正項,Pa

速度修正方程如下:

u=u*+u′

(32)

式中u*--速度計算項,m/s

u′--速度修正項,m/s

修正后的速度u應滿足式(30),則將式(32)代入式(30)可得:

(33)

根據式(19)計算出的速度計算項u*和壓力計算項p*滿足動量方程,因此以下關系式成立:

(34)

(35)

速度修正項u′及壓力修正項p′滿足式(19),所以將式(32)代入式(19),可以得到以下關系式:

(36)

(37)

由此,聯立式(34)與式(36),式(35)與式(37),可以得到以下關系式:

(38)

(39)

忽略相鄰節點的影響,規定以下關系式成立:

(40)

將式(39)、(40)代入式(38),得到速度修正項u′與壓力修正項p′的關系:

(41)

將式(41)代入式(33)可得:

(42)

整理得到:

(43)

為了簡化方程的形式,再次將式(17)逆向代入式(43),可以得到:

(44)

將式(44)整理為下式形式:

(45)

方程系數的計算式為:

(46)

(47)

(48)

(49)

式(44)～(48)為離散的連續性方程。

(50)

當j=3,…,n-2時,離散的連續性方程的系數按式(46)～(48)計算。

(51)

(52)

③ 能量方程

對節點j所在的控制體進行積分,方程如下:

(53)

式(54)可化簡為:

(54)

比焓的計算采用一階迎風格式,按下式計算:

(55)

將式(55)代入式(54),可以得到:

(56)

化簡并整理得到:

(57)

將式(57)整理成下式形式:

(58)

方程系數的計算式為:

(59)

(60)

(61)

(62)

(63)

(64)

當j=3,…,n-2時,方程系數按式(59)～(62)計算。

當j=n-1時,能量方程出口的邊界條件采用了局部單向化假設,即計算結果依賴于n-1節點,不依賴于外界環境,此時,離散的能量方程系數有一項產生了改變:

(65)

至此,完成了所有方程的離散。

2.3.2動態模型的求解

計算上述離散方程時,認為管道入口蒸汽狀態(蒸汽溫度、蒸汽速度、蒸汽壓力)已知,且已知管道長度、公稱直徑(管道外直徑、壁厚)、保溫層厚度、管壁粗糙度、環境溫度、蒸汽與環境的傳熱系數、管道坡度,給定離散時間單位和離散空間單位,可以通過以下步驟對上述模型進行求解,得到某一時刻管道出口的蒸汽溫度、蒸汽壓力及蒸汽速度。

① 將計算的初始時刻稱為t時刻。將t時刻各個節點的ρ、u、h設置為與入口邊界相同。以下開始求解t+1時刻各個節點的p、u、T、ρ、h、μ,為簡化符號,省略以下步驟②～⑦中符號的上標t+1。

② 根據入口邊界條件,按照經驗值給出各個節點t+1時刻的速度迭代初始值u0、壓力迭代初始值p*、溫度迭代初始值T。此處經驗值取值不影響計算結果,只影響收斂速度。根據狀態方程和p*、T計算出t+1時刻各個節點的蒸汽密度ρ、比焓h、動力黏度μ。

③ 由各個節點的u0和p*,根據式(20)～(27)計算離散的動量方程式(19)中的系數。方程系數的計算式中ut+1、pt+1、ρt+1代入的參數即為步驟②中給出的迭代初始值u0、p*、ρ。

⑤ 利用各個節點的p′及式(41)求出各個節點的速度修正項u′,并與第③步求解出的各個節點的u*進行疊加,得到修正后各個節點的蒸汽速度u。

⑥ 將以上步驟得到的各個節點的p和u代入式(58)～(65),求解離散的能量方程。由于采用一階迎風格式,能量方程的求解和動量方程類似,可以由入口邊界起逐一計算各個節點的方程系數和比焓h。求出各個節點的h后,根據狀態方程求出各個節點的新的蒸汽溫度T、蒸汽密度ρ、蒸汽動力黏度μ等。

⑦ 此時的壓力p和速度u滿足連續性方程,但不一定滿足動量方程,因此將本次迭代得到的這一組各個節點的p、u、T、ρ、h、μ作為第②步的迭代初始值,重復步驟③～⑥,直至動量方程、連續性方程和能量方程都收斂,收斂條件為本次迭代與上次迭代結果的相對差值小于10-5,此時各個節點的p和u同時滿足連續性方程及動量方程。此時,計算完成,得到t+1時刻各個節點的p、u、T、ρ、h、μ。

⑧ 若邊界條件隨時間變化,則計算t+2時刻時,代入更新后的邊界條件,和t+1時刻的計算結果,重復上述步驟②～⑦。以此類推,直至計算出每一個時間節點的p、u、T、ρ、h、μ。

以上,完成整個方程的計算。

2.4 實驗驗證

蒸汽管網的水力、熱力工況表現為管網中蒸汽溫度和壓力的分布,選擇天津某典型蒸汽管網穩定運行階段內各熱力站和用戶的流量以及熱力站壓力、溫度監測數據對模型進行驗證。典型管網的管道采用架空布置,經過用戶門前或者路口時進行埋地處理。整個管網按照區域劃分為4條管道,管網拓撲結構見圖2。管道總長7.1 km,供熱面積2.1 km2,熱源與最遠熱用戶22距離3 808 m。采用上述模型和管網的水力計算模型結合,計算各位置的蒸汽壓力和蒸汽溫度。

圖2 典型蒸汽管網拓撲結構

由蒸汽管網水力熱力耦合模型計算結果與實測數據對比可知,所有位置的模型計算值與實測值都較為接近,且離熱源越近的位置模擬值與實測值誤差越小,這是由于隨著熱用戶與熱源距離增加,沿途未被考慮的影響蒸汽管網運行的因素逐漸增加。除個別位置外,模擬值和實測值的相對誤差保持在±8%以內,蒸汽溫度相對誤差絕對值的平均值為3.85%,蒸汽壓力相對誤差絕對值的平均值為2.15%,說明本文建立的蒸汽管道動態模型計算結果可靠。

3 蒸汽管道動態水力熱力計算分析

在實際運行過程中,室外條件會發生變化,蒸汽參數會發生變化,還會對閥門進行調整。因此,本文對蒸汽管道非穩態工況的選取考慮了3種實際情況,分別為環境溫度的變化、熱源溫度的變化、閥門調整導致的管道入口蒸汽壓力的變化,相關計算參數見表5。

表5 相關計算參數

3.1 環境溫度的變化

在蒸汽管網實際運行過程中,環境溫度是周期性變化的,而且不同的地理區域和季節,溫度變化幅度不同,有些地點的全天溫度變化幅度可達10 ℃以上。選取天津冬季某日的逐時環境溫度變化作為非穩態變化項,探究其對蒸汽管道內蒸汽參數的影響。在保證表3中其他參數不變,環境溫度周期性變化的情況下,管道出口蒸汽溫度(簡稱出口溫度)變化見圖3。可以看到,出口溫度的變化相對于環境溫度的變化有一定的延遲。對于算例中1.649 km的蒸汽管道,管道出口蒸汽溫度變化延遲時間約為1.5 h。

圖3 環境溫度變化對出口溫度的影響

實際的蒸汽管道管徑各不同,保溫層厚度也有一定的差異。為探究蒸汽管道不同傳熱系數對出口溫度變化的影響,模擬了不同的傳熱系數下出口溫度變化,見圖4。不同傳熱系數下,出口溫度隨時間的變化規律基本一致,但變化幅度產生了較大變化,傳熱系數越大,出口溫度變化幅度越大。通過以上分析,可以得到這樣的結論:在環境溫度發生變化的情況下,傳熱系數的變化對出口溫度的影響體現為出口溫度變化幅度的變化。

圖4 不同傳熱系數下環境溫度變化對出口溫度的影響

3.2 管道入口溫度的變化

在蒸汽管網運行過程中,由于負荷變化或者環境溫度變化,需要在熱源處提高或者降低蒸汽溫度。熱源處的溫度變化會傳遞到管道的不同位置,導致溫度或者壓力的變化。與熱源距離不同,其熱滯后性也不相同。在保證表3中其他條件不變的情況下,使管道入口蒸汽溫度(簡稱入口溫度)在50 s內上升50 K,管道出口蒸汽溫度和管道出口蒸汽壓力(簡稱出口壓力)變化見圖5。可以看出,入口溫度變化會導致出口溫度變化,且出口溫度變化存在一定的延遲。本案例中1 649 m的管道,出口溫度變化延遲時間大約為130 s。出口壓力也會受到入口溫度變化的影響,在提高入口溫度的同時,出口壓力瞬間降低,幾乎不存在延遲。但由于出口溫度存在延遲變化,導致密度的延遲變化,出口壓力也會慢慢升高,達到新的穩態。出口壓力變化延遲時間與出口溫度變化延遲時間基本相同。在這個過程中,出口壓力達到的最低值低于最終達到穩態的壓力。

圖5 入口溫度變化對出口壓力及出口溫度的影響

溫度波在管道內傳遞比壓力波慢得多,沿管道的溫度傳播存在時間延遲和幅度衰減現象。為探究管道長度對出口溫度變化延遲時間的影響,改變了上述算例的蒸汽管道長度,分別變為原長度(記為L)1 649 m的1/4、1/2、3/4,并進行了模擬,結果見圖6。管道長度不同,延遲時間不同,管道長度為0.25L、0.5L、0.75L、L的出口溫度變化延遲時間分別為50、75、100、130 s,即管道長度每增加約400 m,出口溫度變化延遲時間增加25～30 s。

圖6 不同管道長度下改變入口溫度對出口溫度的影響

3.3 管道入口壓力的變化

調節蒸汽管網中的閥門時,閥門后蒸汽壓力會發生變化,在閥門下游其他位置處的蒸汽狀態參數也會動態變化。為了探究調節閥門對下游蒸汽狀態參數的影響,在保證表3中其他參數不發生變化的情況下,使管道入口蒸汽壓力(簡稱入口壓力)在50 s內下降0.5 MPa,探究入口壓力變化對管道出口蒸汽參數的影響。入口壓力變化對出口溫度和出口壓力的影響見圖7。

圖7 蒸汽管道入口壓力降低對出口蒸汽參數的影響

可以看出,隨著入口壓力降低,出口壓力瞬間降低,與入口壓力變化保持同步,幾乎沒有延遲。隨著入口壓力降低,管道出口蒸汽比焓(簡稱出口比焓)受壓力變化影響,也會發生變化,但其變化過程存在一定的延遲。出口溫度變化與出口比焓變化有較大不同,其先是經歷了一個與出口壓力變化同步的快速下降,然后經歷了與出口比焓變化趨勢相同的緩慢上升。在入口壓力發生變化的階段,出口比焓未發生明顯變化,但是出口壓力的變化導致了同等比焓對應的溫度發生變化。在入口壓力不變的階段,出口壓力不再變化,所以出口溫度只受出口比焓的影響。兩個階段不同的影響因素綜合作用導致了出口溫度的變化。因此,當改變入口壓力時,出口壓力會同步隨之變化,而出口溫度會經歷兩個階段的變化,第一個階段與入口壓力變化趨勢相同且時間上保持同步,第二個階段變化與第一個階段變化趨勢相反,并逐漸趨于平緩。在出口溫度的變化過程中,溫度的最小值或者最大值會超過變化的起始值和終止值。因此在實際調控過程中,需要保證設備能夠承受溫度變化的極值。

4 結論

① 通過天津某典型算例對蒸汽管道動態水力熱力耦合模型進行了驗證,證實了模型的可靠性。

② 環境溫度對出口溫度的影響存在延遲現象,且蒸汽管道的傳熱系數越大,環境溫度對出口溫度的影響越顯著。

③ 隨著入口溫度降低,出口溫度和出口壓力隨之降低,同樣存在延遲現象,且每增加400 m的管長,延遲時間增加25～30 s。

④ 隨著入口壓力降低,出口壓力隨之降低,出口溫度先降低后升高,然后趨于穩定,且溫度變化極值遠小于穩定值,因此在實際調控過程中需要保證設備能夠承受溫度變化的極值。