富 友 解朝蓬 蔣 勁

(1.蘭州理工大學(xué) 能源與動(dòng)力工程學(xué)院,蘭州 730050;2.安徽三聯(lián)泵業(yè)股份有限公司,安徽 馬鞍山 243000;3.武漢大學(xué) 動(dòng)力與機(jī)械學(xué)院,武漢 430072)

在液柱分離水錘現(xiàn)象中,二次升壓現(xiàn)象一直是研究的重點(diǎn)內(nèi)容.這種壓力激升往往伴隨著空泡的生成與潰滅,在閥門處尤為明顯,其壓力激升甚至比根據(jù)Joukowsky公式計(jì)算的最大壓力還要高得多[1].有學(xué)者分析了二次升壓過(guò)程中閥門處空泡的初生—發(fā)育—潰滅的過(guò)程,從而探究空泡潰滅對(duì)流動(dòng)特性的影響[2]；也有試驗(yàn)結(jié)果表明,二次升壓現(xiàn)象在低瞬變強(qiáng)度下會(huì)尤為突出[3],但都對(duì)其內(nèi)部壓力波的變化規(guī)律研究較少.

在伴有液柱分離的瞬變流動(dòng)分析中,空泡的生成使得水錘壓力波的變化程度與單相水錘不同.探索更加精確的瞬變流動(dòng)摩阻模型對(duì)水錘壓力波進(jìn)行捕捉,有助于分析二次升壓壓力波的內(nèi)部變化規(guī)律[4].此外,這種二次升壓的產(chǎn)生會(huì)對(duì)機(jī)組、閥門和儀表產(chǎn)生較大的危害,所以對(duì)其進(jìn)行研究是非常有必要的.本文基于雙流體數(shù)學(xué)模型,采用數(shù)值模擬與試驗(yàn)對(duì)比相結(jié)合的研究方法,構(gòu)建出一種改良的水錘數(shù)值解法,目的是為了探討二次升壓的內(nèi)部流動(dòng)機(jī)理,分析邊界水頭、邊界壓力、水錘波速、管道長(zhǎng)度和管道直徑對(duì)二次升壓的影響.

1 計(jì)算方法

1.1 雙流體數(shù)學(xué)模型

本文采用歐拉-歐拉雙流體模型[5-7],該模型中兩相均為可壓縮.模型包含一個(gè)空泡輸運(yùn)方程、兩個(gè)連續(xù)性方程、兩個(gè)動(dòng)量方程和兩個(gè)能量方程.表達(dá)式如下：

方程中：下標(biāo)1代表氣相；下標(biāo)2代表液相；α為空穴率；u為速度(m/s)；ρ為密度(kg/m3)；p為壓力(Pa)；J u為非恒定摩阻；J s=fu|u|/2D；D為管道直徑(m)；J g為重力項(xiàng),J g=gsinθ；μ為壓力弛豫系數(shù)；λ為速度弛豫系數(shù)；E=e+u2/2；pI為兩相相間壓力,pI=α1p1+α2p2；uI為兩相相間速度,uI=(α1ρ1u1+α2ρ2u2)/(α1ρ1+α2ρ2)；α1+α2=1.

為使上述方程組封閉,引入狀態(tài)方程：

其中：γk為比熱容比,γk=C p,k/C v,k；C p,k為定壓比熱容；C v,k為定容比熱容.π1通常為0,π2表達(dá)式見式(9)：

其中：Tref為試驗(yàn)環(huán)境溫度；psat為飽和蒸汽壓.

波速方程為：

方程(8～10)中：k=1代表氣相；k=2代表液相.

1.1.1 相間弛豫系數(shù)

在使用雙流體模型時(shí),必須考慮瞬態(tài)流動(dòng)中的壓力和速度弛豫問(wèn)題.其決定了模型的計(jì)算穩(wěn)定性[8-9].相間弛豫系數(shù)包括速度弛豫系數(shù)和壓力弛豫系數(shù).

速度弛豫系數(shù)可表達(dá)為：

在壓力弛豫系數(shù)中引入空泡靜態(tài)平衡條件：

其中：S為表面張力；R為泡粒半徑.

進(jìn)而改進(jìn)的壓力弛豫系數(shù)表示為：

1.1.2 非恒定摩阻模型

瞬態(tài)流動(dòng)模型可采用單系數(shù)和雙系數(shù)摩阻模型來(lái)計(jì)算壓力波的衰減.單系數(shù)和雙系數(shù)摩阻模型由于時(shí)間和空間偏微分項(xiàng)的存在,使其在雙流體模型中的應(yīng)用變得復(fù)雜.加權(quán)函數(shù)非恒定摩阻模型[10]是對(duì)Zielke模型和Trikha模型[11]的簡(jiǎn)化.基于該模型的數(shù)值編碼,使用加權(quán)函數(shù),可以節(jié)省大量計(jì)算所需的時(shí)間.

對(duì)Trikha模型簡(jiǎn)化,通過(guò)加權(quán)平均,將瞬態(tài)流動(dòng)阻力模型轉(zhuǎn)化為如下形式：

式中：

其中：υ為運(yùn)動(dòng)粘度(m2/s)；g為重力加速度(m/s2)；y i為速度變量(m/s)；y i中i取1、2、3,對(duì)應(yīng)m i和n i中3個(gè)數(shù).

1.2 求解雙流體模型的數(shù)值方法

1.2.1 Godunov-Rusanov格式

將方程組(1～7)進(jìn)行重新整理,忽略空泡輸運(yùn)方程中的弛豫源項(xiàng).整理后如下：

其中,U=(αρ,αρu,αρE)T,F=(αρ,αρu2+αP,αρE+αP)T,H=(0,pI,pIuI).

礦床內(nèi)發(fā)育了兩組斷裂，一組為NNW向張扭性斷裂，該組斷裂以扭性為主，因扭力是順時(shí)針?lè)较虻?，故多為右斜列式出現(xiàn)，同時(shí)還伴生有小角度的扭裂隙和張裂隙。在這些構(gòu)造裂隙帶中充填有閃鋅礦、輝銻礦、脆硫銻鉛礦和石英、含錳菱鐵礦、含鐵菱錳礦、鐵白云石等金屬礦物和脈石。該組斷裂多成群平行排列出現(xiàn)，地表出露30余條。它是本礦床的主要控礦構(gòu)造，礦脈的產(chǎn)狀及形態(tài)，嚴(yán)格受其控制。第二組為NNE向扭性斷裂，分布在礦床西側(cè)，多充填有乳白色石英脈，局部有輝銻礦化，有切穿NNW向斷裂的現(xiàn)象，為第二期形成的斷裂(圖2)。在礦床主要地段未發(fā)現(xiàn)，對(duì)礦脈無(wú)破壞性。礦區(qū)內(nèi)主要控礦構(gòu)造及性質(zhì)見表1。

求解格式如下：

式中：Δ=?α1/?x,其離散格式為

式中：S=max{|u i±c i|,|u i+1-c i+1|,|u i-1-c i-1|}.

因?yàn)榭张葺斶\(yùn)方程(16)為非守恒型方程,不能采用上述守恒形式格式求解.因此,對(duì)其采用如下求解格式：

1.2.2 邊界條件

如圖1所示,邊界條件的構(gòu)造采用特征線法(圖中橫坐標(biāo)為網(wǎng)格長(zhǎng)度上的網(wǎng)格節(jié)點(diǎn),縱坐標(biāo)為時(shí)間層).

圖1 特征值及特征線在空間坐標(biāo)的分布趨勢(shì)d x/d t=u±c

相容性方程可以表示為：

式中：φR代表當(dāng)特征波速為負(fù)值時(shí),向上游傳遞的信息節(jié)點(diǎn)；φL代表當(dāng)特征波速為正值時(shí),向下游傳遞的信息節(jié)點(diǎn)；φA,φB,φC是不同網(wǎng)格節(jié)點(diǎn)上的值.

2 數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比

管線長(zhǎng)度分為5個(gè)部分,管道長(zhǎng)度在0～0.65 m之間,管道標(biāo)高為0.47 m；管道長(zhǎng)度在0.65～24.83 m 之間,標(biāo)高為0.26 m；管道長(zhǎng)度在24.83～49.17 m之間,標(biāo)高為0.13 m；管道長(zhǎng)度在在73.52～98.56 m之間,標(biāo)高為0 m.計(jì)算時(shí)間步長(zhǎng)為5×10-5s,管道波速為1 240～1 291.8 m/s.在管線沿線有4處壓力檢測(cè)點(diǎn),分別位于管線1/4、2/4、3/4 和4/4 處.文獻(xiàn)[12]給出了相應(yīng)的流體熱力學(xué)參數(shù).使得該模型較好地保證了計(jì)算的穩(wěn)定性和精度.

試驗(yàn)系統(tǒng)原理圖如圖2所示.

圖2 試驗(yàn)系統(tǒng)示意圖

表1 試驗(yàn)具體參數(shù)

數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果對(duì)比如圖3～8所示.

圖3 管線P3處流速為2.74 m/s的水頭隨時(shí)間的變化特征

圖4 管線P1處流速為2.74 m/s的水頭隨時(shí)間的變化特征

圖5 管線P3處流速為1.43 m/s的水頭隨時(shí)間的變化特征

圖6 管線P1處流速為1.43 m/s的水頭隨時(shí)間的變化特征

圖7 管線P3處流速為0.47 m/s的水頭隨時(shí)間的變化特征

無(wú)論是管道末端還是管道中點(diǎn)的模擬結(jié)果都表明：采用非恒定摩阻的新模型,不僅可以捕捉到液柱分離過(guò)程中水擊的最大升壓和降壓值,而且還可以捕捉最大升壓之后的壓力衰減變化.

3 內(nèi)部流動(dòng)特性分析

3.1 水錘二次升壓分析

從圖8可以看出,管道末端的水錘壓力波的二次壓力升高值明顯高于首級(jí)波峰.這種現(xiàn)象是由液柱分離再?gòu)浐袭a(chǎn)生的壓力波反射疊加所導(dǎo)致.

圖8 管線P1處流速為0.47m/s的水頭隨時(shí)間的變化特征

為探究上述二次升壓現(xiàn)象的流動(dòng)機(jī)理,以“閥門-管道-穩(wěn)壓水箱”系統(tǒng)為研究對(duì)象,忽略摩阻影響.選取兩種瞬變流動(dòng)現(xiàn)象,第1種是初始流速為0.47 m/s,初始水頭為70 m 無(wú)二次升壓的水錘現(xiàn)象,其壓力分布如圖9(a)所示；第2種是初始流速為0.47 m/s,初始水頭為50 m 有二次升壓的水錘現(xiàn)象,其壓力分布如圖9(b)所示.

圖9 閥門處兩種瞬態(tài)流動(dòng)壓力分布

圖10和圖11為壓力波在多相瞬變流中的傳播過(guò)程,其中實(shí)線為升壓波,虛線為降壓波.壓力波的傳遞過(guò)程遵循兩個(gè)規(guī)律：(1)壓力波與壁面、閥門等邊界條件相接觸時(shí),其反射波與原波性質(zhì)相同；(2)壓力波與穩(wěn)壓邊界相作用時(shí),其反射波與原波性質(zhì)相異.

圖10 無(wú)二次升壓的水錘現(xiàn)象

圖11 伴有二次升壓的水錘現(xiàn)象

如圖10所示,壓力波1、4、5為升壓波,2、3、6為降壓波.在無(wú)液柱分離情況下,壓力變化按照理論計(jì)算公式(Hmax/min=H0±,H0為邊界水頭)下降與上升,壓力下降與上升的變化范圍如圖12所示(橫坐標(biāo)0代表上游,1代表下游),該條件下并未出現(xiàn)二次升壓現(xiàn)象.

圖12 管線壓力波分布

如圖11所示,當(dāng)末尾閥門快速關(guān)閉,閥門處產(chǎn)生升壓波1向上游傳播,升壓波在上游與穩(wěn)壓罐作用形成降壓波2反射至閥門處,傳播時(shí)間為2L/c.閥門處的低壓波2和3反射導(dǎo)致壓力降低至飽和蒸汽壓,閥門處發(fā)生空化.降壓波3從上游穩(wěn)壓罐處變化為升壓波4反射,升壓波4在管道末端與空泡相遇,此時(shí)波速產(chǎn)生突變,如圖13所示.

圖13 波速隨管線變化

升壓波4在兩相混合區(qū)域產(chǎn)生折射與反射兩種壓力波：第1種壓力波與空泡作用反射形成降壓波5向上游傳遞；第2種壓力波在兩相區(qū)域以升壓波4繼續(xù)向閥門邊界折射傳遞.當(dāng)壓力波4傳遞至閥門處,反射為升壓波6向上游傳遞.在降壓波5和升壓波6的傳遞期間,相同管線位置處降壓波5 早于升壓波6,其變化如圖14(a)所示.

升壓波6傳遞至上游處變化為降壓波8向下游傳遞,降壓波5傳遞至上游變化為升壓波7向下游傳遞.在相同管線位置升壓波7早于降壓波8,其變化如圖14(b)所示.升壓波7在上游處與升壓波6相遇,兩種相同性質(zhì)的正壓波在N 點(diǎn)疊加傳遞至末端閥門處,最終使得末尾閥門處壓力峰值高于理論計(jì)算值.

圖14 管線壓力波分布

3.2 水錘二次升壓隨系統(tǒng)參數(shù)變化

為了判定水錘二次升壓值隨系統(tǒng)參數(shù)的變化,本文選取了邊界水頭、邊界流速、管道直徑、管道長(zhǎng)度、和水錘波速5個(gè)參數(shù)來(lái)進(jìn)行內(nèi)部流動(dòng)特性分析.

3.2.1 邊界水頭對(duì)二次升壓的影響

圖15統(tǒng)計(jì)了二次升壓值(二次峰值)在特定初始流速下隨邊界水頭的變化,圖中橫坐標(biāo)為邊界水頭(H0)與Joukowsky升壓(ΔH=(cv0)/g)的比值,縱坐標(biāo)表示二次升壓值(H2max)與Joukowsky 升壓值(ΔH)的比值.通過(guò)分析發(fā)現(xiàn),在特定的初始流速下,二次升壓值隨著初始水頭的增大呈現(xiàn)出先增大后減小的趨勢(shì),并且在H0≈0.8ΔH處會(huì)取得極值.此結(jié)論應(yīng)用在水錘泵的輸水過(guò)程中,可以使其獲得更高的揚(yáng)程.

圖15 二次升壓值隨邊界水頭的變化

3.2.2 邊界流速對(duì)二次升壓的影響

表2統(tǒng)計(jì)了二次升壓值在特定初始水頭下隨初始流速的變化.可以看出,在初始水頭為24.8 m 條件下,初始流速在0.30～0.82 m/s之間,二次升壓值隨著初始流速的增大而升高,當(dāng)初始流速大于1.02 m/s,二次升壓現(xiàn)象便不再發(fā)生.在初始水頭為41.8 m的條件下,初始流速為0.30 m/s,并沒(méi)有發(fā)生二次升壓現(xiàn)象,因?yàn)樵诖藯l件下的初始邊界壓力大于Joukowsky升壓值,不能發(fā)生液柱分離,從而使得一次峰值與二次峰值相等.而初始流速在0.47～1.14 m/s之間,二次升壓值隨著初始流速的增大而增大,當(dāng)初始流速大于1.35 m/s,二次升壓現(xiàn)象便不再發(fā)生.在初始水頭為71.8 m 條件下,初始流速在0.64～1.35 m/s之間,二次升壓值隨著初始流速的增大而增大,當(dāng)初始流速大于1.70 m/s,二次升壓現(xiàn)象便不再發(fā)生.以上說(shuō)明,在液柱分離過(guò)程中,二次升壓值隨著初始流速的增大而增大,當(dāng)初始流速大于臨界流速,二次升壓便不會(huì)發(fā)生.因此,在工程應(yīng)用中,如果要避免因液柱分離而產(chǎn)生的二次升壓,可以通過(guò)數(shù)值模擬來(lái)選擇最優(yōu)工況.

表2 二次升壓值隨初始流速變化(單位：Pa)

3.2.3 管道直徑、管道長(zhǎng)度和水錘波速對(duì)二次升壓的影響

表3統(tǒng)計(jì)了在初始邊界H0=0.8ΔH條件下,二次升壓相對(duì)值隨管道直徑、管道長(zhǎng)度和管道波速的變化.在理論情況下,壓力波在水中的傳播速度不超過(guò)1 450 m/s,實(shí)際波速值要小于理論值,并且為使得該條件下發(fā)生液柱分離,選取波速最小值為1 100 m/s.

表3 二次升壓相對(duì)值H2max/ΔH 的變化情況

可以看出,隨著管道直徑的增加,二次升壓相對(duì)值保持不變；隨著管道長(zhǎng)度的增加,二次升壓值相應(yīng)增大；隨著水錘波速的增大,二次升壓值隨之增大.

3.2.4 二次升壓試驗(yàn)驗(yàn)證

為了驗(yàn)證二次升壓值隨邊界水頭的變化規(guī)律,本節(jié)對(duì)圖15的模擬結(jié)果進(jìn)行了相應(yīng)的實(shí)驗(yàn)驗(yàn)證.

表4 試驗(yàn)具體數(shù)據(jù)

首先進(jìn)行了初始流速為0.348 3 m/s,邊界水頭為25.428 m 的水錘試驗(yàn).在此基礎(chǔ)上,保持初始流速基本不變,取邊界水頭為H0=0.6ΔH、H0=0.7ΔH、H0=0.8ΔH和H0=0.9ΔH,分別對(duì)相應(yīng)試驗(yàn)2、3、4、5進(jìn)行驗(yàn)證.試驗(yàn)結(jié)果如圖16所示,可以看出,在H0=0.8ΔH處,二次升壓值最大.

圖16 二次升壓值隨邊界水頭變化

圖17 管線中點(diǎn)處水頭隨時(shí)間的變化特征

上述試驗(yàn)結(jié)果對(duì)圖15的模擬結(jié)果進(jìn)行了驗(yàn)證.從水錘壓力波隨管線分布可以看出,由于摩阻的存在,升壓波在閥門處反射向上游傳遞的過(guò)程中持續(xù)衰減,并且在管線中點(diǎn)處二次升壓現(xiàn)象不再發(fā)生.

4 結(jié)論

本文通過(guò)對(duì)瞬態(tài)流動(dòng)中的二次升壓現(xiàn)象進(jìn)行分析,得到了以下結(jié)論：

1)在伴有液柱分離的水錘中,空泡的生成與潰滅會(huì)使得管道內(nèi)部發(fā)生壓力波的疊加.管道中兩種升壓波疊加沿管線傳遞至閥門處,使得閥門處流體進(jìn)一步被壓縮,從而產(chǎn)生二次升壓現(xiàn)象,這種二次升壓值會(huì)高于理論計(jì)算值.

2)二次升壓值在H0≈0.8ΔH處取得極值,隨初始水頭增大呈現(xiàn)先增大后減小的趨勢(shì).二次升壓值隨著初始流速的增大而增大,當(dāng)初始流速大于臨界流速,二次升壓便不會(huì)再發(fā)生.二次升壓值隨水錘波速和管道長(zhǎng)度的增大而增大,管道直徑對(duì)二次升壓值無(wú)影響.

3)通過(guò)5 組試驗(yàn)結(jié)果驗(yàn)證了二次升壓值在H0≈0.8ΔH處取得極值.在升壓波從閥門處向上游傳遞的過(guò)程中,管線阻力對(duì)其衰減影響顯著.