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

富 友，蔣 勁 ，李燕輝，應(yīng) 銳

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

0 引 言

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

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

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

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

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

1 改進的雙流體計算模型

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

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

1.1 閉合條件

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

相間的內(nèi)能變化如式（10）～式（12）所示。

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

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

1.2 管壁變形的計算方法

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

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

2 數(shù)值計算方法

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

2.1 空穴輸運方程的數(shù)值計算方法

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

2.2 雙曲型算子計算方法

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

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

2.3 對于不守恒項的數(shù)值計算方法

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

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

2.4 相間弛豫的計算方法

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

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

2.4.1 有限相間弛豫現(xiàn)象計算方法

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

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

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

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

2.4.2 無限相間弛豫現(xiàn)象計算方法

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

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

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

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

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

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

2.5 邊界條件

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

3 數(shù)值計算及試驗比對

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

3.1 辛普森關(guān)閥水錘試驗

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

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

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

表1 辛普森試驗熱力學(xué)參數(shù)Table 1 EOS parameters for liquid and gas for Simpson’s experiment

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

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

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

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

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

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

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

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

3.2 管線上游水錘關(guān)閥試驗

為模擬離心泵突然失電后關(guān)閥導(dǎo)致的管路系統(tǒng)中壓力瞬變，本文在武漢大學(xué)水力機械過渡過程實驗室，對此做了相應(yīng)試驗進行對比驗證。試驗管路系統(tǒng)布置如圖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。試驗中采取電動蝶閥進行快關(guān)試驗。計算中選用無限弛豫計算方法進行計算。試驗數(shù)據(jù)如表2所示。

圖5 試驗管路系統(tǒng)布置圖Fig. 5 Layout drawing of test piping system

表2 試驗初始數(shù)據(jù)Table 2 Initial conditions of experiment

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

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

表3 關(guān)閥試驗熱力學(xué)參數(shù)Table 3 EOS parameters for liquid and vapor water for experiment

圖6 關(guān)閥試驗結(jié)果與計算結(jié)果對比Fig. 6 Comparison between experimental data and calculation data in water hammer experiment

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

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

4 結(jié) 論

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

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

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

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