康競瀾， 劉 昉， 侯慶志， 何軍齡， 林 磊

(1. 天津大學 建筑工程學院， 天津 300350;2. 西北農(nóng)林科技大學 水利與建筑工程學院， 陜西 楊凌 712100;3. 蘇州熱工研究院， 江蘇 蘇州 215004)

在諸如火電站、核電站等電力或化工工業(yè)系統(tǒng)中，會存在大量用于輸送高溫高壓蒸汽的管道。當管路系統(tǒng)由于維護、檢修或其他原因而停止工作時，管道中的高溫蒸汽因溫度降低而凝結成的小水團會在管路系統(tǒng)低洼處逐漸積聚存留，形成靜態(tài)的液體段塞。一旦系統(tǒng)重新啟動，管道上游的瞬時高壓蒸汽會將這些段塞瞬間加速，像槍膛內(nèi)的子彈一樣快速射出[1]。當高速運動的段塞撞擊到彎管、三通或閥門等管端結構時，將產(chǎn)生巨大的沖擊力，可能會對管路造成嚴重的破壞，從而影響整個系統(tǒng)的正常工作。因此，對管道關鍵部分(如彎管處)的高速運動段塞的沖擊研究尤為重要。目前關于在水平管道或者垂直立管中的穩(wěn)態(tài)段塞流研究已經(jīng)有了比較成熟的體系和方法[2-3]，而對于彎管中單個段塞流的研究尚未完全展開，尤其是段塞流在彎管處的高維沖擊現(xiàn)象值得深入研究。

由于管道中存在的段塞的高風險性，因此需要對其運動狀態(tài)進行分析。之前的研究表明，作用在彎管上的力主要是由流體運動方向改變所造成的動量轉換產(chǎn)生的。段塞流在管道中的流動實際是三維的，但以往的研究主要是基于一維模型，具有較大的局限性。因此為了精準地描述段塞流的運動狀態(tài)，本文建立了單個段塞在三維空管中的運動及其在彎管處的沖擊模型，并將數(shù)值模擬結果與試驗結果進行了比較，真實準確地分析了段塞在空管中的運動過程以及高速運動段塞對彎管的沖擊。

1 段塞流沖擊彎管研究現(xiàn)狀

最近，針對段塞流沖擊彎管的問題進行了綜述研究，本節(jié)將從試驗研究與理論研究兩個方面簡要闡述對本次研究比較重要的成果。

1.1 試驗研究

Fenton 等[4]通過試驗最早研究了直徑為25 mm的傾斜管道中，彎管因單個液體段塞沖擊而產(chǎn)生的作用力。當段塞運動距離大于其自身長度5倍時，氣體將擊穿段塞從而導致沖擊力大幅降低。Neumann等[5]發(fā)現(xiàn)當管道氣液比小于20%時，便不會發(fā)生分層流向段塞流轉變的過程，而由此對彎管處的沖擊力也可忽略不計。Bozkus等[6]在改進的Fenton 試驗裝置基礎上，發(fā)現(xiàn)相對較短的段塞沖擊壓力曲線呈現(xiàn)單峰，而相對較長的段塞呈現(xiàn)雙峰，并把此現(xiàn)象的產(chǎn)生歸結于閥門打開時的水錘效應。Owen等[7]進行了段塞流沖擊孔板的物理試驗，記錄了孔板處的壓力變化，測得了不同長度段塞流到達孔板時的速度范圍在30 m/s到60 m/s不等。

1.2 理論研究

對于段塞流沖擊速度的計算，目前已有多種理論模型。Fenton 等通過穩(wěn)定狀態(tài)下段塞流的名義速度予以描述，其表達式為

(1)

式中：ρ表示下游的氣體密度；A表示管道的橫截面積；m為質量流量。Neumann等提出一個相對復雜模型，考慮了驅動氣體的動力學性質，卻忽略了段塞流的質量損失。Kayhan等[8]運用單步特征線法，分別考慮了恒定滯留系數(shù)以及通過試驗數(shù)據(jù)擬合得到的滯留系數(shù)，并由此計算了質量損失。最近，Tijsseling等[9]通過考慮段塞前后氣液界面的速度差，將段塞的速度分成三個部分，改進了Bozkus的數(shù)值模型，取得了較好的效果。

穩(wěn)態(tài)段塞流對固壁的沖擊與穩(wěn)態(tài)射流的沖擊[10]類似，其表達式為

(2)

其中ρs為段塞密度，Vs為段塞沖擊速度。但該模型計算的沖擊力遠低于試驗測得的實際沖擊力。Fenton[11]在其1989年論文中提到，在工程應用中彎管處的沖擊力可由以下公式計算所得

(3)

上式中第一項為以聲速c傳播的壓力波在彎管處沖擊后產(chǎn)生的力，其持續(xù)時間極短，試驗中難以被觀測，計算時一般被忽略。第二項為流體在彎管處的摩擦力，數(shù)值較小也常被忽略不計。第三項為段塞上下游的壓差，即段塞的驅動力。第四項為由瞬時沖擊速度對壁面產(chǎn)生的作用力。

最近，Hou等[12]基于控制體理論，基于彎管內(nèi)側發(fā)生的流動分離現(xiàn)象，提出了一種新模型，沖擊力的幅值與變化趨勢與試驗結果基本一致，其在x方向的沖擊力表達式為

(4)

式中：Pe為上下游的壓差；Cc為彎管處的流體收縮系數(shù)；Ke為水頭損失系數(shù)。

綜上所述，目前普遍認可的段塞流對彎管沖擊力的簡化計算公式為

(5)

等式右邊兩項分別為段塞受到的驅動力以及動量轉換產(chǎn)生的沖擊力。

2 試驗裝置及數(shù)值模擬

為了盡可能模擬實際工況，Bozkus等[13]借鑒了以往研究者[4-6]的試驗裝置，進行了更大尺度管徑(10 cm)的試驗。閥門下游主管道被設置成傾斜的，注水后可在彎肘處形成截留段塞，準確模擬了真實蒸汽管道中的截留段塞分布。相對于段塞的運動距離，段塞長度被設計得更長(超過其運動距離的五分之一)，以避免其在還未到達彎頭時就已經(jīng)破碎。與以前研究結果相比，壓力測量系統(tǒng)的時間響應表現(xiàn)良好，動態(tài)壓力測量顯示了相似趨勢。試驗得到了彎頭處的沖擊壓力時程曲線。

2.1 試驗裝置

文獻[13]中的試驗裝置如圖1所示，該裝置包括位于上游容積為0.5 m3的圓柱形儲氣罐、一個DN100的球閥，一段內(nèi)徑為10 cm，長85 cm的垂直管段用于連接球閥，以及一根相對水平面呈4.6度，內(nèi)徑為10 cm、末端帶有彎頭的12 m長鋼管。通過壓縮機使儲氣罐充滿壓縮空氣，以模擬發(fā)電廠管道中的壓縮蒸汽。在試驗過程中，僅控制初始段塞長度以及儲氣罐壓力，當儲氣罐中空氣溫度及管道中水溫達到環(huán)境溫度時，對所有工況進行試驗。當球閥快速打開時(開啟時間約為16 ms)，上游的壓縮氣體將驅動段塞在管道中高速運動直至撞擊到彎管處。位于彎管處的壓力傳感器記錄了沖擊壓力的時程曲線，數(shù)據(jù)采樣頻率為1 000 Hz。

(a) 試驗裝置原理圖

2.2 數(shù)學模型

針對段塞流問題，假設液體不可壓縮，其連續(xù)性方程為

?·V=0

(6)

動量方程為

(7)

假定氣體與液體互不相容且彼此占據(jù)連續(xù)空間，因此變化的氣液界面可以采用VOF模型予以描述。為真實還原段塞流在管道中的運動形態(tài)，提高模擬精度，本文采用RNGk-ε湍流模型。RNGk-ε模型的湍流動能方程和湍流動能耗散方程分別表述如下

(8)

(9)

式中：k為湍流動能，表示速度波動的變化量；ε為湍流動能的耗散，表示速度波動對時間的耗散率；μ為運動黏度；μt為湍流黏度；αk=αε=1.39，Cε1=1.42，Cε2=1.68。

2.3 三維CFD模型

針對2.1中的試驗工況，本節(jié)構建了一個三維CFD模型，運用有限體積法離散求解控制方程，壓力與速度的耦合使用PISO算法，求解器采用Fluent 17.0版本。

前處理采用ICEM軟件劃分出O-grid型的結構化網(wǎng)格，這種網(wǎng)格容易實現(xiàn)流動區(qū)域的邊界擬合，比非結構化網(wǎng)格更適用于管道中段塞流的模擬計算。同時進行了邊界層網(wǎng)格劃分，所得y+值主要分布在32～184。圖2顯示了管道整體及局部的網(wǎng)格劃分圖。

(a) 管道整體模型

由于管道傾斜將在段塞頭部產(chǎn)生較長的自由液面，因此為了準確計算段塞初始長度，采用了如圖3所示的方法，利用三角函數(shù)關系進行了簡化。段塞長度L0為圖中x、y的長度之和。

壓力監(jiān)測點設置在彎管中心軸線方向上。通過對驅動壓力為5bar，段塞長度為3 m的算例進行試算，比較了80萬～230萬等7種不同網(wǎng)格數(shù)目的結果(見表1)，確定150萬為合適的計算網(wǎng)格數(shù)。

圖3 段塞標記方法(不按比例)

表1 網(wǎng)格數(shù)量無關性驗證

將帶邊界層與無邊界層的網(wǎng)格算例進行試算，比較了兩組具有代表性的算例結果，如表2所示，結果發(fā)現(xiàn)無論是壓力峰值還是峰值時刻，相對誤差均在2%以內(nèi)。這個誤差相較于模擬誤差是可以接受的，而計算耗時大約是3倍的關系，這將顯著占用計算資源以及花費更長的無必要時間。因此綜合考慮后，決定采用無邊界層的網(wǎng)格進行了數(shù)值模擬。

表2 邊界層網(wǎng)格無關性驗證

針對之前部分學者提到的試驗過程中氣罐內(nèi)壓力的衰減情況[15]，本文考慮了恒定驅動壓力以及隨時間衰減的驅動壓力兩種情況。由于無法準確估計試驗中氣罐壓力的衰減情況，根據(jù)以往經(jīng)驗[16]，假定壓力為線性衰減且衰減率為30%，該步驟通過UDF完成。壓力衰減的表達式為

Ps=P0(1-0.3ts/t0)

(10)

式中：Ps為衰減后的實時驅動壓力；P0為初始驅動壓力；ts為段塞的運動時間；t0為未考慮壓力衰減情況下，段塞即將撞擊到彎管處的時間，通過試算迭代得到。

由試驗得到的時程曲線可知，達到壓力峰值的沖擊時間極短，為毫秒量級。因此，為了精確捕捉壓力峰值，綜合考慮試驗數(shù)據(jù)以及試算結果后，確定時間步長為10-4s。試驗中數(shù)據(jù)采樣頻率為1 000 Hz，為避免采樣頻率不同產(chǎn)生的影響，將上文提到的基于兩種驅動壓力情況計算得到的數(shù)據(jù)分別每隔10組取平均值，組成新的壓力時程數(shù)據(jù)。

3 數(shù)值結果與分析

3.1 模擬結果驗證分析

圖4分別顯示了L0=4 m，P0=2 bar工況下的段塞沖擊彎管前后時刻的二維及三維云圖，可以清晰地觀察到段塞尾部氣體對段塞內(nèi)部的侵蝕以及段塞在彎管處的流動分離情況。同時發(fā)現(xiàn)段塞在沖擊彎管前后，即使是在主體核心區(qū)，液柱與管道接觸面附近仍然摻雜著大量氣體。

表3對比了數(shù)值模擬結果與物理試驗結果，其中L0表示段塞初始長度，P0表示段塞驅動壓力，計算值表示每隔10組取平均值后的峰值壓力；UDF計算值表示引入壓力衰減后的峰值壓力；UDF計算值取均值表示引入壓力衰減后，再每隔10組取平均值后的峰值壓力。表3中同時給出了各計算值與試驗值的相對誤差。對比分析表明，大部分數(shù)值模擬結果與試驗結果具有良好的一致性，并且氣罐壓力衰減速率大小對模擬結果有較大影響。對于段塞沖擊彎管的12個工況，采用恒定驅動壓力模擬得到的壓力峰值相較于試驗值偏大。引入UDF模塊進行氣罐壓力衰減設置后，結果更接近于試驗值。另外，對計算數(shù)據(jù)取均值后，得到的預測值變小，這是因為壓力峰值持續(xù)時間僅為兩三個時間步長且衰減極快。

部分數(shù)值模擬結果與試驗結果差異的原因主要有以下兩個方面：

(1) 對于L0=6 m的工況，不考慮氣罐壓力減小時，計算所得壓力值與試驗值吻合較好，而設置壓力衰減后模擬值變小。這是因為對于初始長段塞，氣罐用于驅動其撞擊到彎頭處所需氣體體積小，氣罐壓力衰減值遠低于采用的30%壓力衰減率。對于初始長度較短，驅動壓力較大的情況，模擬得到的壓力峰值相較于試驗值偏大。例如，L0=4 m，P0=5bar工況下，不論經(jīng)過何種處理，計算值依然偏大較多，主要是因為假定的氣罐壓力衰減不夠大且不夠快。當L0=5 m時，模擬值與試驗值比較接近。

(2) 由于壓力峰值持續(xù)時間極短，在壓力傳感器的數(shù)據(jù)采樣頻率為1 000 Hz的條件下，會存在某些持續(xù)時間極短的壓力峰值捕捉不到的情況。而當數(shù)值模擬的時間步長為10-4s，即采樣頻率為10 000 Hz時，則避免了上述缺點，同時計算更易收斂。尤其是當段塞初始長度較短，而驅動壓力較大時，此時段塞速度極快，峰值壓力的衰減也更快，捕捉到的試驗值均小于模擬值。而當段塞初始長度較長，驅動壓力較小時，此時峰值壓力的衰減相對慢一些，這也是在此情況下，試驗值與數(shù)值模擬值基本吻合的原因。

(a) 沖擊前

表3 CFD與試驗壓力峰值對比

對于3種不同初始長度液柱，數(shù)值模擬得到的壓力時程曲線與文獻結果對比如圖5～圖7所示。為與Bozkvs等提出壓力時程曲線的單位保持一致，這里采用psi作為壓強單位。結果表明，模擬所得整個沖擊過程(幅值和時程變化趨勢)與試驗結果大致相同。同時，以L0=4 m為例，進行了頻響分析，如圖8所示，結果表明，監(jiān)測點壓力主頻均位于0.5～2 Hz，壓力幅值隨著頻率的增加而快速減小。P0=5 bar時，壓力脈動的主頻幅值最大，表明較大的驅動壓力將會導致更劇烈的管道振動。

(a) 試驗

(a) 試驗

(b) 數(shù)值模擬

3.2 段塞沿程質量脫落

段塞在管道中運動時，由于摩擦力及重力的存在，會發(fā)生質量脫落，導致段塞在沖擊彎管前所具有的質量遠小于其初始質量。為了進一步研究其質量脫落規(guī)律，定義段塞質量脫落率α=1-Ls(t)/L0，其中Ls(t)為運動段塞長度。

利用后處理軟件CFD-POST，對Ls(t)進行了分析，并畫出了質量脫落率隨段塞運動距離的變化曲線(見圖9)，其中橫坐標代表段塞尾部距彎肘中心的長度。圖示結果表明：對于本文所研究的工況，給定初始液柱長度時，段塞質量脫落率與其運動距離成正比，其曲線擬合的函數(shù)關系為：y=0.033x+0.113,R2=0.996。脫落率與驅動壓力P0基本無關，即同樣初始長度的段塞，不論其后方的驅動壓力是多少，最終到達彎頭時的液柱質量都是定值。

圖8 L0=4 m工況下段塞沖擊彎管壓力頻響曲線

圖9 L0 =6 m的沿程脫落率曲線

段塞到達彎管處的脫落率與初始長度的關系如圖10所示。顯然，對于初始長度越短的段塞，因運動距離長，其到達末端彎管處的質量脫落率越大，但該脫落率與段塞初始長度并非線性關系，而是二次曲線關系，擬合曲線的函數(shù)關系為：y=0.012 9x2-0.276x+1.584,R2=0.997。對于短的段塞，質量脫落占整體段塞的比例接近90%，導致段塞在撞擊時被加速到極高的速度，這是其峰值力較大的根本原因。

圖10 不同初始段塞長度的脫落率曲線

3.3 彎管處段塞沖擊模式

如前文所述，彎管處的沖擊壓力值與段塞沖擊速度呈正相關。為了準確估計段塞對末端彎管的沖擊作用，本文通過進一步的數(shù)值模擬，分析了段塞對彎管的沖擊模式。對于2 bar驅動壓力下，不同初始長度的段塞，結果顯示沖擊壓力的變化趨勢存在三種情況：

(1) 壓力遞減(L0=4 m)，如圖11所示。段塞到達彎管時具有極高的速度，導致壁面對段塞的反作用力遠遠大于驅動力，加速度方向改變。速度峰值和壓力峰值都出現(xiàn)在初始時刻。本文前述模擬的工況，都屬于該沖擊模式。

圖11 L0= 4 m彎管處的壓力時程圖

(2) 壓力遞增(L0= 11 m)，如圖12所示。段塞到達彎管的速度較低，反作用力小于驅動力，且隨著液體由出口處不斷流出，其速度持續(xù)增加。雖然彎頭處的壁面反力增加，但段塞的加速度方向不發(fā)生改變。因此，速度和壓力峰值都出現(xiàn)在沖擊過程的末尾時刻。此外，在沖擊過程中，因部分氣體截留在彎頭處，隨著氣團被逐漸沖刷釋放，產(chǎn)生了小幅壓力振蕩(見圖12)。

圖12 L0=11 m彎管處的壓力時程圖

(3) 壓力振蕩(L0= 7.5 m)，如圖13所示。段塞對彎管沖擊后，壁面反力與驅動力基本相當。此外，與第二類模式類似，因彎管處截留氣團的沖刷釋放，沖擊壓力產(chǎn)生了大幅振蕩。

圖13 L0=7.5 m彎管處的壓力時程圖

對于2～20 bar之間的其他五種驅動壓力，所得沖擊模式跟2 bar條件下的模式基本相同。這三種沖擊模式可通過公式(5)予以進一步闡釋。當段塞運動到彎管時，因方向改變，在彎管處產(chǎn)生沖擊壓力，其大小取決于段塞的沖擊速度。當段塞運動到彎頭時的速度較大時，因速度引起的壁面反力遠大于段塞驅動力，段塞具有反向加速度，所以速度逐漸降低，沖擊壓力降低。當段塞運動到彎頭時的速度較小時，因速度引起的壁面反力小于段塞驅動力，段塞加速度方向不發(fā)生改變，速度繼續(xù)增加，沖擊壓力增加。當段塞運動到彎頭時的速度適中時，因速度引起的壁面反力與段塞驅動力相當，段塞加速度為零，所以速度和沖擊壓力基本不變。

對于實際工程應用而言，我們主要關注的是第一種沖擊模式，即壓力峰值出現(xiàn)在沖擊的初始時刻，此時段塞對管端結構的沖擊載荷大，破壞力強，這是應該避免的。

3.4 壓力峰值計算方法驗證

為了探究高速運動段塞對彎管沖擊壓力的計算方法，提取了12個工況下的段塞沖擊速度，并依據(jù)式(5)進行了數(shù)值計算，計算結果如表4所示，表中的Pe值為壓力衰減30%后的驅動壓力。表4中同時給出了理論計算值與試驗值及數(shù)值模擬值之間的相對誤差。整體來看，理論計算結果與數(shù)值模擬結果相吻合。由于段塞的瞬時速度值是通過等距離構造管道截面，再取截面的平均速度所得到的，考慮到其誤差，因此可以認為通過該方法能夠相對準確地計算瞬時沖擊壓力。

表4 數(shù)值計算的壓力峰值

4 結 論

本文針對單個段塞在空管里的運動及其對管端彎頭的高速沖擊問題，在總結分析以往理論與試驗研究的基礎上，結合已有試驗數(shù)據(jù)，建立了段塞流瞬態(tài)沖擊三維數(shù)值模型，并利用該模型模擬了文獻中不同工況下段塞的運動及其對彎管的沖擊過程，得到了相應的沖擊壓力時程曲線。對于不同初始長度段塞，采用合適的壓力衰減率，數(shù)值模擬結果具有較高的準確度，沖擊壓力時程曲線與試驗結果基本吻合。在此基礎上，通過數(shù)值試驗，分析了段塞運動過程中的質量脫落率，提出了段塞對彎管的3種沖擊模式，以及段塞對彎管沖擊壓力的簡化計算公式。所得主要結論包括：

(1) 對于初始長度給定的段塞，在驅動壓力主導時，運動單位距離時的質量脫落率相同，最終到達彎頭時的液體質量為定值。

(2) 段塞沖擊過程中彎管處的壓力變化有三種模式：(i) 壓力峰值出現(xiàn)在沖擊初始時刻，然后快速減小；(ii) 壓力突增后，隨速度的增大，壓力增大，峰值出現(xiàn)在沖擊末尾時刻；(iii) 壓力突增后，在該值附近上下波動。

(3) 對于單個運動段塞對彎管的瞬態(tài)沖擊壓力，簡化計算公式(5)能夠準確地計算沖擊壓力峰值，理論計算值與試驗及數(shù)值模擬結果基本吻合，該計算方法在實際工程應用上具有指導意義。