改進雙流體模型計算有液柱分離的管路水錘壓力

富 友，蔣 勁 ，李燕輝，應 銳

（1.蘭州理工大學能源與動力工程學院，蘭州 730050；2.武漢大學水力機械過渡過程教育部重點實驗室，武漢 430072）

0 引 言

當管路中的流動狀態發生變化時，會產生水錘現象。水錘現象的發生、發展程度受管路材料特性和流動介質屬性的影響[1]。就管路特性方面，材料的彈性模量對管路中的波速傳遞有著明顯影響，特別對于薄管壁流動而言，水錘波的傳遞速度將會發生明顯變化[2]，對于介質屬性而言，密度、壓力及熱力學參數均對相應波的傳遞有明顯的影響[3]。近些年在水錘計算研究方面，主要針對水錘波的傳遞過程中受管路幾何結構的影響，管路水錘計算的流固耦合過程及水錘壓力波的衰減過程方面進行了較為詳細的研究[4-7]，但帶有液柱分離現象的水錘計算一直為計算的難點，有較多的學者對該種現象進行過相應的研究分析，并取得了一定的成果[8-12]。當壓力降低至飽和蒸氣壓力以下時，由于相變導致的液柱分離現象將會對壓力響應的準確性造成更大的變化及影響[13]。

傳統的水錘計算方法，通常采用連續性方程和動量方程進行構建計算模型，并采用特征線法進行求解計算。對于計算帶有液柱分離的水錘現象時，通常采用離散空穴模型DVCM或DGCM模型進行計算，計算中當壓力降低至飽和蒸氣壓時，將所計算的區域作為計算邊界來計算相應的空穴率的變化[14-16]。

對于DVCM模型，Simpson 等[17]采用試驗及數值計算進行驗證，結果表明空穴率大于 10%時，壓力峰值及持續時間均出現了明顯的誤差。

Bergant等[13]指出當空穴率大于10%時，由于交界面的問題，將導致泡狀流動下的離散空泡的交界條件不再適用。Wylie曾將一小部分自由氣體引入計算模型中，進而阻止數值震蕩的出現[18]。也有學者對帶有此類不可溶性氣體的計算方法（DGCM 模型）進行研究，但自由氣體體積的含量具有不可確定性[1]，同時體積的量級對于計算的壓力變換有明顯的影響[19]，該方法也受到一定質疑。不可溶性氣體的存在，使得水錘現象可以從兩相流的角度模擬計算[20-21]。在氣液兩相流動計算中，雙流體模型是能夠較為準確描述氣液兩相流動特征的計算方法，該方法認為各相的動力學特性不完全相同，模型由 1個空泡輸運方程、2個連續性方程、2個動量方程及2個能量方程組成，在保證雙相的可壓縮性的前提下，方程具有良好的雙曲性[22]，可為氣液兩相瞬變流動計算提供基礎。

本文選用改進的雙流體模型進行水錘計算模型的構建，并通過構建合理的相間作用形式及高精度的數值解法來描述帶有液柱分離的氣液兩相流動現象。本文首先對于改進的雙流體模型的數學特性進行分析，并選用全空化模型對相間質量變化進行計算，構建出一種考慮相間質量變化的氣液兩相瞬變流的計算方法。其次對方程中的不守恒項進行了單獨處理并選用矢通量分裂(steger-warming flux vector splitting, SW-FVS)與無振蕩，無自由參數，耗散算法 (non-oscillatory containing, no free parameters and dissipative scheme, NND)相結合的時間上為一階，空間上均為二階的數值方法進行數值計算。最后結合2個實驗室的試驗數據對2種弛豫方法對于計算結果的影響進行討論分析。

1 改進的雙流體計算模型

水錘計算模型是基于改進雙流體模型建立，與傳統的雙流體模型的區別在于方程組中添加了管壁截面參數變量，這一改變也使得該計算模型更適用于管路中的氣液兩相流動計算。方程模型對液相和氣相分別采用不同的壓力、空穴率和密度進行描述。模型由 1個空泡輸運方程，2個連續性方程，2個動量方程及2個能量方程組成。改進的模型方程的雙曲型部分與文獻[3，22-23]所示相似，其中一維的數值模型建立如式（1）～式（7）。

式中αk為空穴率；ρk為密度，ρ為混合密度，ρ1為相間密度，kg/m3，k=1代表氣相，k=2代表液相；u為速度，m/s；λ為弛豫系數；Pk為壓力，Pa；E為彈性模量，MPa；e為管壁厚度，m；m ˙g為質量傳遞率，kg/s；ck為介質波速，ckg為管內波速，m/s；D為直徑，m；A為面積，m2；γk為比熱比；τ為剪切應力，N/m2；Ek為總能量，J，Ek=ek+1/2uk2；ek為內能，J；θ 為管路傾斜角，(°)，γk為比熱比。

1.1 閉合條件

由于計算中考慮了能量的變化，因此需要引入狀態方程保證方程的封閉，這里采用修正的氣體狀態方程即SG-EOS[22]。

相間的內能變化如式（10）～式（12）所示。

比熱比γk，參考壓力πk，焓qk，定壓比熱Cvk，定容比熱Cpk參考文獻[23-24]。在計算中π1認為大小為0（等同于理想氣體狀態方程），這一參數的取值帶入波速方程后可得與聲速方程一樣表達式，γk與 π2采用式（13）、（14）進行計算。

式中Tref為溫度，K；psat為飽和蒸汽壓力，Pa。

1.2 管壁變形的計算方法

計算中僅考慮一維流動，不考慮管道軸向的形變，僅考慮管壁徑向上形變。同時管道內的壓力變化僅僅對管壁內部橫截面積造成影響，對于管壁厚度不造成影響。參照文獻[19]，對管壁內的橫截面積采用式（15）～（16）進行計算。

式中 A0為初始管壁面積，m2；P∞為初始壓力，Pa，D0為初始管壁直徑，m；cn1為管壁材質參數，v為泊松比。

2 數值計算方法

為得到更高求解精度下的計算結果，需要選用二階或更高階的求解格式進行數值求解計算，通常對于顯示求解格式而言，高階的求解格式往往會產生數值震蕩，為了保證在處理壓力震蕩時的精確性及穩定性，本文選用NND格式進行求解計算。對于氣液兩相雙流體模型而言，方程呈現非守恒形式，并包含不守恒項，因此采用高階守恒格式的數值求解方法時，需要對相應的不守恒方程及不守恒項進行單獨處理。

2.1 空穴輸運方程的數值計算方法

方程(1)而言，由于無法將其寫成守恒型格式，因此需要對空泡的輸運方程進行單獨處理，參照修正的Godunov scheme方法[22]，對方程（1）進行離散，其離散形式表現為式（17）～（18）。

2.2 雙曲型算子計算方法

對于方程（1）～（7）而言，方程中守恒部分采用SW-FVS與NND相結合的方法進行計算。

NND是張涵信院士提出的一種數值計算方法[25]。其計算公式如式（21）～（24）所示。

2.3 對于不守恒項的數值計算方法

將方程（2）～（3）、（4）～（5）依照式（21）寫成相應的離散格式。

將方程（25）帶入方程（26）中，根據文獻[26]，可得表達式（28）與（29）。此時式（29）即為方程的不守恒項提供了相應的求解方法。

2.4 相間弛豫的計算方法

對于雙流體模型下的氣液兩相瞬變流而言，氣液兩相分別具有各自的狀態方程，并以此來描述各相的流動特性的變化。當產生瞬變波動時，從波速方程中可以看出，不同相間將會產生不同量級的波[27]。

隨著壓力波的傳遞，兩相必將會產生相異的壓力。對于流動而言，介質不可能長時間保持彼此相異的環境壓力，需要一個平衡的狀態來描述該過程，同時這個過程也將導致兩相體積及內能的相應變化。通常對于兩相壓力的平衡過程是一個二維或三維的平衡過程。對于一維管道流動而言，采用更高維度的計算方法將使計算變的更加復雜，因此需要構建適用于一維求解的計算方法。

2.4.1 有限相間弛豫現象計算方法

由方程組（1）～（7）可知，方程右側存在計算壓力弛豫的源項，即λ（P1-P2）。對于采用兩相瞬變流動的計算方法，可采用有限弛豫壓力系數來進行構建[28-30]，其構建形式一般如式（30）形式。

式中AI為單體體積下泡粒表面積，N為泡粒數量。

這種方法需引入泡粒交界面的計算方程，同時需要考慮泡粒的脫落與融合。對于泡粒的直徑，數量密度的假設將會大大限制計算的穩定性。并考慮到泡粒面積的改變，需要引入相間質量的描述方式。對于交界面質量變化的計算方法，一般認為流動中存在 3種物質，即空氣、水蒸氣和水。當壓力降低至飽和蒸氣壓力以下時，控制體內將發生空化現象。選用全空化模型進行計算相間質量變化，全空化模型是由Singhal等[30]提出。模型考慮流動中的不可溶性氣體，同時具有經驗系數少，計算準確等特點。該種模型在CFD商業軟件中有較為廣泛的應用，其蒸騰和凝結過程采用式 （31）～（34）計算。

式中 kt為湍動能，m2/s2；fk為質量分數；Cc、Ce通常取0.01和0.02；σ為張力，N/m。

2.4.2 無限相間弛豫現象計算方法

當認為壓力弛豫現象的計算過程遠遠小于時間步長時，可認為壓力弛豫系數變換為無限大，此時可對每次計算后的結果進行修正，以此來達到對于弛豫現象的描述。在使用無限弛豫現象修正方法的過程中，不需要采用空化模型來描述相間的質量變化，該種方法已對相間的作用做出修正和描述。修正時源項對于空間的偏微分對計算不產生影響，參考相關文獻[22]有式（35）～（39）。

將式（35）帶入式（38）、（39）中，可以得到式（40）。

將式（40）進行積分可以得到式（41）。

式中上標0表示雙曲算子后的計算結果，上標*表示弛豫過程后的計算結果。將方程進行進一步整理可得式（42）。

將狀態方程（9）、相間壓力計算方程進行帶入，同時由于π1通常取為0，因此有式（43）。

同時從式（38）、（39）可知，過大的泡粒直徑將導致弛豫系數增大，進而導致弛豫系數對于相間能量的反向修正[31-32]，也驗證了2.4.1中的結論。

2.5 邊界條件

對于管道氣液兩相流動的上、下游的邊界條件，采用特征線法進行構建，在構建時一般需要知道上游、下游的速度或壓力的變化，以此來計算其他流動參數的變化趨勢。其中相容性方程采用式（45）進行插值計算。

3 數值計算及試驗比對

為驗證該計算模型的適用性及準確性，選用了 2個實驗室的不同試驗數據，分別與數值模擬的結果進行對比，試驗包括辛普森的關閥水錘試驗[33]、武漢大學的水錘試驗。

3.1 辛普森關閥水錘試驗

辛普森關閥試驗由四部分組成，分別為上游穩壓水池、管線、尾閥和下游水池組成。圖 1為辛普森關閥試驗示意圖。

圖1 辛普森試驗示意圖Fig. 1 Sketch of Simpson’s experiment

管線總長為36 m，直徑為 0.019 05 m，厚度為 0.001 588 m，管道在9 m位置處有一90o轉接頭，管道末端連接低位水池，在管線36 m位置安裝有一球閥。管線具有1 m高的坡度，在管線沿線有3處壓力檢測點，分別位于9、27、36 m。計算初始空穴率為10-7，上游穩壓水箱壓力為31.88 m，流速為1.125 m/s，管線阻力損失為2.8 m，飽和蒸氣壓力為3 032.17 Pa，波速為1 280 m/s，流體主要參數均參考文獻[32]。水和空氣的密度分別為998.2、2.52 kg/m3，飽和水蒸氣的密度為0.077 14 kg/m3。表1為所對應的流體固有熱力學參數。模擬計算以尾閥處速度邊界進行讀入，關閥時間為0.043 07 s。圖2為辛普森關閥試驗結果與計算結果比對圖。

表1 辛普森試驗熱力學參數Table 1 EOS parameters for liquid and gas for Simpson’s experiment

圖2 辛普森試驗結果與計算結果比對Fig. 2 Comparison between model and experiments in Simpson’s experiment

由圖 2可知，當采用有限弛豫系數計算方法時，如果不加入空化模型，那么計算結果與試驗數據將出現極大的差異。從采用有限弛豫系數計算方法并引入空化模型的計算中可以看出，當泡粒數量（單位體積下泡粒數）選定為 1013時，計算結果與試驗有較好的相似度，且泡粒數量選擇量級與文獻[34-35]描述一致。

在計算中發現，泡粒數量對于計算結果有明顯的影響，當泡粒數量小于1013時，如圖3中所示會出現次級波峰的響應時間出現偏差。當泡粒數量大于 1013時，會因為單位濃度下的泡粒表面積過大導致相間的動量傳遞超過真實計算值，造成內能和壓力的反向修正，從而導致計算中斷[31]。

圖3 不同泡粒數量對辛普森試驗的影響Fig. 3 Effect of different particle concentration in Simpson’s experiment

從圖 2計算結果中可以看出，雖無限弛豫系數的計算結果較當泡粒數量選定為 1013時的計算結果有一定偏差，但與試驗值的偏差同樣較小。同時通過圖2d空穴率的比對可以看出，2種計算方法計算的空穴率的變化雖然有一定偏差，但空穴率的變化趨勢相似性較高。

圖4 采用不同弛豫計算方法時氣液兩相壓力差Fig. 4 Difference of gas-liquid two-phase pressure difference with different relaxation calculation methods

圖 4為采用無限弛豫系數修正前氣相壓力與液相壓力差，與采用有限弛豫系數并選用空化模型下的壓力差對比圖。由圖 4可知，采用無限弛豫系數修正前，在壓力低于飽和蒸氣壓力的空化初生區，氣相壓力始終大于液相壓力。由式（42）可知，采用無限弛豫系數，正是對相間的能量變化進行修正，能量的變化雖然沒有體現出相變，但從本質上來說，可以看成是伴隨著能量變化的另一種相變表達形式，因此該種修正方法下不需要再次加入相變修正模型。同樣從圖4中可以看出當在0.3 s左右，正壓力返回至末端時，對于采用有限弛豫系數修正的方法依然造成了壓力相修正的偏差。當選用的泡粒數量繼續增大時，將會出現反向修正的情況。

由于無限壓力弛豫計算方法為一種強制性的修正計算方法，即在每個時間步長后，認為氣相、液相壓力均達到平衡狀態，而事實上壓力的弛豫時間通常為1 ms。若計算時間步長小于1 ms時，由于計算方法的強制修正，使得未到達平衡狀態下的兩相流動強制處于平衡狀態，這也導致了空穴體積的修正偏大，使得在低壓區的空化模型作用被大幅度消減。

3.2 管線上游水錘關閥試驗

為模擬離心泵突然失電后關閥導致的管路系統中壓力瞬變，本文在武漢大學水力機械過渡過程實驗室，對此做了相應試驗進行對比驗證。試驗管路系統布置如圖5所示。試驗由水箱、水泵、閥門、壓力傳感器、電磁流量計及水平管線組成，其中試驗段管線總長為 97.6 m，管徑為0.1 m，管壁材料為鍍鋅鋼管，管壁厚為3 mm，折算管壁聲速約為1 260 m/s；末端為高位恒壓水箱，水位高度為 26 m。水和空氣的密度分別為 998.2和 2.52 kg/m3，飽和水蒸氣的密度為 0.077 14 kg/m3。試驗中采取電動蝶閥進行快關試驗。計算中選用無限弛豫計算方法進行計算。試驗數據如表2所示。

圖5 試驗管路系統布置圖Fig. 5 Layout drawing of test piping system

表2 試驗初始數據Table 2 Initial conditions of experiment

對于壓力的測試，采用閥門位置（壓力傳感器P1點所在位置）和0.5倍管線位置（壓力傳感器P2點所在位置）進行數據采集。試驗采用武漢超宇測控技術有限公司生產的 CY 3088型號壓力變送器，測量范圍為 0～5.0 MPa（絕對壓力），測量精度為 0.5%，流量計為KROHNE生產的IFM 4080K電磁流量計，測量范圍為0～80 m3/h，測量精度為0.5%，試驗系統中的儀表精度等級直接決定了試驗臺的測試精度，在本章試驗系統中試驗臺總測試系統的系統誤差為±0.933 8%。試驗中選用的測量儀器輸出均為電流信號，采用數據采集卡將電流設備轉換為數值信號，并通過傳輸存儲于計算機指定位置上。表 3為關閥試驗熱力學參數。將計算結果和試驗結果進行比對，如圖6所示。由圖6可知，計算結果與試驗結果吻合程度較高，驗證了算法的可行性。

表4為試驗對比誤差。由表4可知，試驗1和試驗2壓力波峰值均超過Joukowsky水錘定律（a△u/g）計算值，證明瞬變流動中產生了明顯的液柱分離現象，同時通過壓力峰值與試驗值的計算偏差可以看出計算值與試驗值偏差較小，驗證了數值方法的準確性。

表3 關閥試驗熱力學參數Table 3 EOS parameters for liquid and vapor water for experiment

圖6 關閥試驗結果與計算結果對比Fig. 6 Comparison between experimental data and calculation data in water hammer experiment

表4 試驗對比誤差Table 4 Error comparison between experimental and calculation

綜上 2個實驗室的不同試驗可以看出無限弛豫系數無空化模型的計算方法可以較好的模擬管路系統中，由于關閥所造成的水力過渡的壓力瞬變，并可以較為精準的預測管路系統的過渡過程。同時相對于以往的氣液兩相計算模型[37]而言，本文所選用的計算模型，去掉了泡狀流假設，可以使得計算的穩定性得到更好的保障。

4 結 論

本文基于改進的雙流體氣液兩相流計算模型，對管道中由于水錘造成的液柱分離的氣液兩相瞬變流動進行計算分析比對，計算中考慮了管壁截面隨壓力變化時對波速的影響，并對方程中的不守恒項及不守恒方程進行了單獨處理，以確保模型計算的穩定性。通過對比兩種弛豫計算方法對于計算的影響，并采用辛普森關閥試驗與管線上游水錘關閥試驗進行分析比對論證，結論如下：

1）采用有限弛豫參數計算時，需考慮由于空化引起的相變，同時需選擇合理的泡粒數量以保證計算精確度。

2）采用無限弛豫參數計算時，空化模型對于計算結果影響較小，可直接對水錘現象造成的液柱分離現象進行分析模擬，但會造成空穴率計算結果過大。

3）采用改進的雙流體氣液兩相流計算模型，并構建合理的數值計算方法及弛豫系數，可以對帶有明顯液柱分離現象的管路氣液兩相瞬變現象進行預測計算，該方法計算精度高，對于工業管道中的水錘防護計算具有較好的應用價值。