富友 牛雪天 敏政

(蘭州理工大學 能源與動力工程學院，甘肅 蘭州 730050)

在輸油、輸水管道系統中，由于油泵、水泵機組及閥門眾多，存在發(fā)生事故斷電、閥門間啟閉不協調的可能，這將導致管路系統形成復雜的瞬變流動現象。這種流動現象發(fā)生、發(fā)展的過程會導致整個系統瞬態(tài)過程中水力特性的改變，從而可能導致管道的爆破或塌陷、水力元件的嚴重破壞、輸油/水管線大范圍泄漏，甚至發(fā)生爆炸與火災。鑒于這些嚴重的后果，需從機理上了解這一瞬變現象的本質即彌合性水擊。當水擊壓力波呈現與流動相反方向傳播時，會出現水頭壓力的降低，當下降至飽和蒸汽壓附近時，介質會發(fā)生汽化，形成汽、液兩相流。一旦形成的氣(汽)相發(fā)生壓縮潰滅，會產生明顯的局部射流高壓，這種伴隨氣(汽)相的發(fā)展與潰滅的水擊現象被稱為彌合性水擊。

Baltzer[1]、Adamkowski等[2]分別對彌合性水擊現象進行了理論分析及實驗拍攝，分析中發(fā)現這種現象從物理過程上歸屬于空化范疇，但其發(fā)生、發(fā)展的過程更為復雜：由于液體始終處于被動拉伸狀態(tài)，其空穴生長、壓縮潰滅過程的時間更為短暫，同時在壓力作用下空泡被壓縮潰滅后形成射流所造成的破壞性更大。Tanaka等[3]發(fā)現：當彌合性水擊出現時，壓力上升階段的次級壓力波峰值將比首級峰值更高。Adamkowski等[4]發(fā)現，在彌合性水擊與非彌合性水擊的過渡區(qū)域，會產生一個時間極短的壓力脈沖，這種現象出現時的邊界條件較為特殊，時間較短，極易被忽視。隨著研究的進一步深入，發(fā)現彌合性水擊空穴壓縮潰滅時所產生的壓力激升，并非在第二個壓力周期初期就呈現，而是在高壓平穩(wěn)末端形成更高的射流高壓[5]，該壓力峰值會超過Joukowsky理論計算值[6]，這也正是射流水擊現象。對于這種射流流動現象的詳細研究，是全面揭示彌合性水擊內部機理及流動規(guī)律中急需解決的問題。

為此，本文采用雙流體模型，通過對弛豫過程的相間作用機理進行分析，提出了一種改進的、可用于描述泡粒變化的弛豫過程計算方法，通過氣液兩相激波管算例驗證該計算方法的穩(wěn)定性及準確性，并針對彌合性水擊的內部流動過程進行了分析，揭示了彌合性水擊中的射流-緩沖機理。

1 數學計算模型

本文的數學計算模型采用氣-液兩相雙流體模型進行構建。雙流體模型由一個空穴的輸運方程、兩個連續(xù)性方程、兩個動量方程、兩個能量方程組成[7- 9]。模型中認為雙相均為可壓縮介質，并認為氣相、液相分別具有不同的壓力、密度?？昭ㄝ斶\方程(1)的存在保障了方程組呈現雙曲型特性，模型方程如下：

(1)

(2)

(3)

(4)

式中：k=1代表氣相，k=2代表液相；I代表液相與氣相的過渡區(qū)域；α為空穴率；p為壓力；ρ為密度；單位體積下的總能量E=e+u2/2，e為內能，u為速度；μ為壓力弛豫系數。

1.1 相間參數的構造

由于兩相之間的密度差異較大，因此需要一個相間區(qū)域(相間邊界)來進行過渡。如圖1所示，這是一個使氣體和液體處于平衡狀態(tài)的區(qū)域。該區(qū)域的相間參數可以用如下方程來表述[8,10- 11]：

(5)

(6)

1.2 狀態(tài)方程

狀態(tài)方程建立了壓力、密度和能量之間的關系。該方程使用了加強氣體方程的修正形式[10,12]：

(7)

(8)

Tref為環(huán)境溫度，psat為飽和蒸汽壓力。

圖1 氣液相間結構的示意圖

2 數值計算方法

2.1 改進的Godunov- HLL數值計算方法

把方程(1)-(4)整理成如下形式，忽略空穴輸運方程中由于弛豫現象產生的源項：

(9)

(10)

式中，U=(α1ρ1,α2ρ2,ρu,α1ρ1E1,α2ρ2E2)，F(U)=(α1ρ1u,α2ρ2u,ρu2+α1p1+α2p2,(α1ρ1E1+α1p1)u,(α2ρ2E2+α2p2)u)，Q(U)=(0,0,0,pIU,-pIU)。

在方程求解中，用到了HLL近似黎曼求解器。HLL近似黎曼求解器的數值表達式如下：

(11)

式中：FL、FR、UL、UR分別代表網格節(jié)點i(i=±1/2)左、右兩端的邊界值；S+和S-分別代表不同的波速取值，其表達式為[7,13]

(12)

數值計算中的差分格式的構建方法如下：

(13)

(14)

(15)

由于空穴輸運方程不能寫成守恒形式，因此本文采用改進的Godunov方案進行構建：

(16)

2.2 邊界構建方法

(17)

對全微分方程進行積分并取一階近似，可以得到如下表達式：

T-1(UP-UA)Δt

(18)

(19)

式中，P、A、B的取值如圖2所示。

圖2 特征值及特征線在空間坐標的分布趨勢

進一步將方程(18)、(19)改寫為如下形式：

T-1(UP-UA)Δt=T-1DAΔt

(20)

3條特征線上數值采用以下形式進行插值計算：

(21)

式中：UB代表當特征波速為負值時，向上游傳遞的信息源節(jié)點；UA代表當特征波速為正值時，向下游傳遞的信息源節(jié)點；UC代表A、B間網格上的數值節(jié)點。

2 弛豫計算模型

本文通過結合有限弛豫法和無限弛豫法的計算過程[14- 15]，推導出一種能夠計算泡粒變化的雙流體相間弛豫構建方法[16]。根據空泡靜力學原理，當壓力降低至飽和蒸汽壓力附近時，流動區(qū)域并不一定發(fā)生空化，僅當液體壓力降至psat-2S/R(S為表面張力，R為泡粒半徑)以下時，空泡才可能開始繼續(xù)增長。通常在室溫下，直徑為10 μm的空泡可以承受14 kPa的張力[17]。因此，結合泡動力學相關特性，將-2S/R引入弛豫構建方法中。由于弛豫過程與空間變量無關，它僅是時間的函數[7]，因此，可以將原始方程式

(22)

(23)

(24)

變換為如下形式：

(25)

(26)

(27)

式中，上標*表示弛豫變量，上標0表示計算變量。

根據空泡靜力學原理[7,17]，氣相和液相壓力之差可以表述為

(28)

將方程進一步整理，可以得到：

(29)

(30)

(31)

3 驗證算例

3.1 水氣兩相激波管

為了研究弛豫模型及數值算法的穩(wěn)定性及準確性，選用水氣兩相激波管算例[7]進行驗證。算例中，水氣兩相激波管被分成兩部分，管內相對長度0～0.7范圍內充滿水，管內相對長度0.7～1.0范圍內充滿空氣。由于雙流體模型的數學特性，需要默認兩側分別存在體積分數在10-8以下量級的第二相。水氣兩相激波管的具體參數如表1所示。數值計算解與解析解在t=2.29 ms時的對比如圖3所示，數值計算中CFL值定為0.8。從圖中可以看出，在t=2.29 ms時，激波管內的混合密度、壓力、速度分布與解析解具有良好的一致性。根據混合密度在0.8～0.9相對位置處臺階式分布匹配關系可以看出，該弛豫模型及數值算法具有良好的穩(wěn)定性及準確性。

表1 水氣兩相激波管參數

3.2 彌合性水擊實驗分析-辛普森水擊實驗

為驗證計算模型對帶有相間質量變化物理過程計算的有效性，選用典型的彌合性水擊實驗-辛普森水擊實驗[5]進行對比分析。辛普森水擊實驗主要由4部分組成：上游水箱、管道、球閥和下游水箱。管道長度為36 m，內徑為0.019 05 m，壁厚為0.001 588 m。管道在12.5 m位置有一個90°的彎頭。

在下游水箱附近的36 m位置處存在一處球閥，球閥可以通過手動操作實現快速關閉。整個管道的坡度為1 m，沿管道有3個壓力變送器，分別位于

圖3 水氣兩相激波管數值解與解析解對比

9、27和36 m處，系統中波速為1 280 m/s。實驗管道系統如圖4所示，實驗基本參數如下：流量0.332 m/s，上游水頭34.54 m，水頭損失0.33 m，摩擦系數0.031 5，理論水頭變化43.332 m，關閥時間0.036 08 s。狀態(tài)方程中的修正熱力學參數見表2，Q為單位質量的能量。

通過表2的修正熱力學參數來表達方程組中的狀態(tài)方程。整個瞬變過程通過尾部球閥處的信息輸入來計算整個管路中的瞬變特征，其方法為：首先在尾閥處測量壓力變化，然后根據式(32)計算速度邊界條件。當計算的速度變?yōu)?時，速度邊界將保持為0[5]。

圖4 實驗系統示意圖

表2 狀態(tài)方程中的修正熱力學參數[12]

(32)

從表2可知，實驗中閥門關閉過程所消耗的時間小于壓力波傳播回閥門所需的時間(2L/c)，因此可以將其歸類為直接水擊實驗。由于水擊造成的壓力變化值大于上游邊界壓力，因此它是具有水柱分離功能的直接水擊。

弛豫模型與辛普森水擊實驗之間的水頭變化比較如圖5所示，從圖中可知，計算值與實驗值具有很好的一致性。低壓(0.15 s)后的二次和三次壓力波的響應時間和幅度都具有良好的精度。數值計算經網格無關性驗證后，考慮計算成本，使用了720個網格單元，計算時間步長為10-7s。

3.3 彌合性水擊射流-緩沖機理分析

由茹科夫斯基直接水擊計算公式可知，當發(fā)生直接水擊時，壓力的升高值是一定的，即等同于第一個水擊壓力激升的數值60.32 m。從圖5中可以看出，低壓后的二次和三次壓力波的數值遠高于60.32 m。為闡明這一流動現象的內部機理，取不同時刻下的管路速度變化與壓力變化特征進行對比分析。為了減少管路高程及管路摩擦損失對流動特性的影響，將辛普森水擊實驗中的高程、管路摩擦損失忽略，并將管路簡化為無量綱的相對長度，0代表上游起點，1代表下游終止位置。不同時間節(jié)點下的管路速度與壓力變化趨勢如圖6、7所示。

圖5 辛普森水擊實驗水頭變化信息

圖6 管路的內流速度變化趨勢

圖7 管路的壓力變化趨勢

經分析發(fā)現，對于產生的二次壓力高于理論計算值的彌合性水擊而言，其流動規(guī)律可以分為直接水擊階段(0.01～0.08 s)、緩沖水擊階段(0.09～0.16 s)、射流水擊階段(0.17～0.24 s)。

1)在直接水擊階段下，速度的變化趨勢為：先下降至0，再呈現反向的流動(速度為負值)，如圖 6(a)所示。從圖7(a)可知：在0.01～0.05 s期間，由于邊界處的流體突然停止運動，管路流體速度逐漸降低，動能轉化為壓力能，壓力則隨時間的延長而升高，使液體被壓縮，形成一股壓縮波向上游進行傳播；在0.06～0.08 s時，由于壓差作用，流動開始自下游至上游反向流動，此時形成膨脹波向下游傳遞，使壓力恢復。

2)在緩沖水擊階段的0.09～0.11 s時段，速度開始反向減小，逐漸接近0，但速度變化并未繼續(xù)遵從直接水擊階段中速度恢復至0的流動特征，而是仍具有由下游向上游的方向(負值)，如圖6(b)所示。這一流動變化的產生，是由于此時空穴泡粒在管路末端處不斷析出，形成反向緩沖作用。從圖7(b)可以看出，相對長度在0.4～1.0內，管路內流體絕對壓力降低至飽和蒸汽壓力以下，驗證了空穴泡粒存在并起到緩沖反向作用的說法。

在0.12～0.14 s時段，流動趨勢開始變化為由上游向下游進行，速度變化為正值。在相對長度為0.9～1.0的范圍內，0.14 s時的速度出現較大的梯度變化，此區(qū)域伴隨發(fā)生了空穴泡粒的壓縮潰滅，但此時速度并未恢復至初始邊界速度(0.332 m/s)，而同時期的壓力變化信息顯示，此區(qū)域壓力由低壓逐漸恢復至初始邊界壓力，驗證了存在空穴泡粒壓縮潰滅及其對流動起到緩沖作用的說法。

在0.15～0.16 s時段，速度相繼出現較大的梯度變化，且出現兩種速度變化疊加傳播：由下游向上游的增速波和由下游到上游的減速波。其中增速波的傳播要早于減速波。這種現象是由于管路末端的空穴泡粒發(fā)生壓縮潰滅，造成整個管道中的流體產生瞬態(tài)加速，形成增速波由下游向上游逐漸傳遞。而當邊界處的流體受到阻礙突然停止運動，壓強升高使液體被壓縮，形成升壓波、減速波向上游傳播。由圖7(b)可以看出，同時期存在明顯的壓力升高，但其升高值約為60.32 m，遠小于直接水擊所產生的壓力值77.87 m，因此將這一區(qū)間定義為水擊緩沖區(qū)。

3)由圖6(c)可以看出，在射流水擊階段的0.17～0.20 s時段，速度出現兩種變化的疊加傳播：由上游向下游的增速波和由上游到下游的減速波。其中增速波的傳播要早于減速波。上一階段的增速波在上游邊界處反饋后形成反向傳播(即上游向下游傳播)。與此同時，壓強在管內逐漸升高，由于上下游存在較大壓差(邊界壓力為34.54 m)，導致流體反向流動，形成減速波由上游傳播至下游。而由于增速波對邊界的沖擊形成了射流作用，導致壓力的進一步上升，在0.20 s時形成了90 m左右的二次升高值(如圖7(c)所示)，這一數值的出現，影響了采用直接水擊方法計算壓力升高值的準確度，因此這也是彌合性水擊中最應該注意的部分。

在0.21～0.24 s時段，開始重復0.09～0.12 s 時段的內部流動變化規(guī)律，但由于二次射流壓力波的存在，同時忽略了流動中的能量損耗，形成了與直接水擊壓力波不同的壓力表現規(guī)律。

4 結論

本文采用雙流體模型，通過對弛豫過程的相間作用機理進行分析，提出了一種改進的可用于描述泡粒體積變化的弛豫模型，采用氣液兩相激波管算例驗證了該計算方法的準確性，并使用該方法對彌合性水擊中的射流-緩沖現象進行了內部流動分析，得到以下結論：

(1)對于低壓過后所產生的二次壓力高于理論計算值的彌合性水擊而言，其流動按規(guī)律可以分為直接水擊階段、緩沖水擊階段和射流水擊階段。

(2)對于緩沖水擊階段，由于空穴泡粒的生成，降低了初次水擊對邊界處的壓力沖擊，使得水擊呈現被緩沖狀態(tài)，壓力升高值小于理論計算值，此時射流作用消失。

(3)對于射流水擊階段，空穴潰滅時產生的射流加速并未直接作用于初次水擊對邊界處的壓力沖擊，而是經過反向傳遞后，使得本應恢復至初始速度值的速度獲得了一個加速作用，這一加速過程產生明顯的波動向水擊波傳播的反方向傳播，在初次水擊基礎上造成二次壓力升高，此升高值要高于理論計算值。