汪怡然，俞曉東，石 林，張 健

(河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098)

1 研究背景

近年來，由于水力元件眾多、沿線地形地質(zhì)情況復(fù)雜等原因，長距離供水工程的運(yùn)行事故頻發(fā)，管路水錘是產(chǎn)生爆管事故的主要原因之一[1]。目前，長距離供水管路水錘問題的研究主要針對單相流或氣液兩相流(液柱彌合)[2-4]，關(guān)于固液兩相流的研究較少。而我國主要河流的含沙量較高(例如黃河的平均含沙量高達(dá)26.5 kg/m3[5])，含沙水的流動(dòng)特性又與清水相差甚遠(yuǎn)。如何求解含沙水錘壓力、評(píng)估泥沙顆粒參數(shù)對高含沙流域長距離供水管路運(yùn)行安全的影響，成為一個(gè)不可忽視的問題。

由于泥沙在液體中受慣性力、重力、阻力的共同作用，其滾動(dòng)、跳躍、懸移的運(yùn)動(dòng)過程難以描述，學(xué)者們通常傾向于在宏觀上劃分固液混合物(本文為含沙水)類型，但基于固相濃度、顆粒直徑、沉降速度等參數(shù)制定的判定標(biāo)準(zhǔn)尚未統(tǒng)一。本文借鑒費(fèi)祥俊[6]關(guān)于管道斷面垂向濃度的Wasp[7]定量指標(biāo)C/CA(C為管頂下0.08倍管徑處的固相體積分?jǐn)?shù)；CA為管路中心處的固相體積分?jǐn)?shù))為判定標(biāo)準(zhǔn)：定義C/CA>0.8即為均質(zhì)流(偽一相流)，反之則為非均質(zhì)流。盡管錢寧和張瑞瑾均認(rèn)為高含沙水懸移質(zhì)屬于受湍流支持的非牛頓流體，一旦流速下降便會(huì)沉降；費(fèi)祥俊認(rèn)為以含沙量高低來區(qū)分含沙水毫無必要，即便是含沙量900 kg/m3的高含沙水懸移質(zhì)(C/CA=0.50～0.73)仍屬于非均質(zhì)流[6]。

目前針對含沙水錘的研究較少，但在水力輸運(yùn)領(lǐng)域涉及固液混合物水錘(以下簡稱漿體水錘)的研究成果較為豐富，這對構(gòu)建含沙水錘數(shù)值計(jì)算方法起到了重要的參考作用。韓文亮等[8-11]結(jié)合連續(xù)性方程、動(dòng)量方程、管路和固相顆粒的彈性模量與變形量方程，推導(dǎo)了均質(zhì)流與非均質(zhì)流的水錘壓力波速公式和壓力增量公式，采用特征線法計(jì)算得到了漿體水錘的混合物波速和瞬態(tài)壓力，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合程度較高。這一方法得到了國內(nèi)外專家學(xué)者的廣泛運(yùn)用[12-14]，成為最主流的漿體水錘數(shù)值計(jì)算方法。此外，周云龍等[15]基于CFD(計(jì)算流體力學(xué))技術(shù)動(dòng)網(wǎng)格方法探究閥門關(guān)閉時(shí)間對漿體水錘壓力增量的影響規(guī)律，結(jié)果顯示最大壓力增量隨著關(guān)閥時(shí)間的延長而減小，這與清水間接水錘的規(guī)律一致。李艷等[16]基于AMESim軟件探究粗顆粒均質(zhì)流垂直管道漿體水錘特性，得出了隨著管道長度的增加、管道直徑的增加、顆粒濃度的減小，垂直管路漿體水錘最大壓力增量減少的結(jié)論。盡管漿體水錘與含沙水錘的原理一致，但含沙水(多為低濃度非均質(zhì)流)與漿體(多為高濃度均質(zhì)流)的流動(dòng)特性存在明顯差異；且已有研究表明，非均質(zhì)流的壓力波速高于均質(zhì)流，且隨著固相濃度的增加差值逐漸增大[8]。因此，含沙水錘的數(shù)學(xué)描述不能直接套用漿體水錘的計(jì)算模型，應(yīng)當(dāng)結(jié)合含沙水的管道流動(dòng)特性與非均質(zhì)流水錘壓力波速公式進(jìn)行研究。

目前，針對管路的清水水錘計(jì)算和漿體水錘計(jì)算大多采用特征線法。本文基于費(fèi)祥俊阻力模型和漿體水錘波速公式，推導(dǎo)了非均質(zhì)流含沙水錘特征線方程；同樣采用特征線法求解了RPV系統(tǒng)的關(guān)閥水錘壓力，探究了含沙量、顆粒直徑與初始流速對含沙水錘壓力的影響，總結(jié)了一般性規(guī)律。研究結(jié)論可為含沙流域長距離供水管路水錘計(jì)算與防護(hù)設(shè)計(jì)提供參考，指導(dǎo)工程實(shí)踐。

2 含沙水流動(dòng)特性與含沙水錘求解

2.1 含沙水流動(dòng)特性相較于含沙水流型、流態(tài)及速度分布等流動(dòng)特性，管路水頭損失對于水錘過程的影響較大。在計(jì)算含沙水管路的水頭損失之前，仍需對本文界定的含沙水類型做出一定的假設(shè)：

(1)假設(shè)不同固相濃度的含沙水均為非均質(zhì)流，忽略顆粒沉積對管道沿程阻力系數(shù)的影響。

(2)假設(shè)不同固相濃度的含沙水均為牛頓流體，忽略高含沙水懸移質(zhì)對流動(dòng)特性的影響。

(3)假設(shè)含沙量S、顆粒平均密度ρs、顆粒平均直徑ds可以表征含沙水的固相濃度、密度和顆粒直徑，忽略顆粒形狀(假設(shè)為球形顆粒)、顆粒級(jí)配、固相濃度及密度分布對含沙水流動(dòng)特性的影響。

目前，基于經(jīng)驗(yàn)或半經(jīng)驗(yàn)理論，專家學(xué)者們提出了諸多的管路阻力模型和水頭損失計(jì)算公式[17]，在此不做贅述。本文選用兼顧懸移質(zhì)和推移質(zhì)作用的費(fèi)祥俊阻力模型[18]，其中非均質(zhì)流含沙水管路水頭損失hm的計(jì)算式為：

(1)

式中：h0為清水的水頭損失；hs為顆粒運(yùn)動(dòng)引起(以沉降為主)的附加水頭損失；在清水項(xiàng)h0中，f0為清水沿程阻力系數(shù)(在工程中常用管道糙率n表示)；g為重力加速度；D為管徑；v為含沙水流速；ρL為液相(水)的密度；ρm為含沙水平均密度；α為考慮懸移質(zhì)對紊流抑制作用的阻力影響系數(shù)，可以基于下述表達(dá)式計(jì)算[19]：

α=1-0.4lgμr+0.2(lgμr)2

(2)

式中μr為混合物黏度與清水黏度的比值，可根據(jù)Einstein稀懸液黏度公式μr=1+2.5S/ρs計(jì)算得出。

(3)

由于阻力系數(shù)本身不易測得，前人基于大量的試驗(yàn)數(shù)據(jù)，構(gòu)建了阻力系數(shù)CD和顆粒雷諾數(shù)Red的關(guān)系圖，其中顆粒雷諾數(shù)Red關(guān)于液相黏度μL的計(jì)算式為：

(4)

(5)

(6)

這樣，聯(lián)立式(1)(2)(6)即可計(jì)算得到長距離供水管路非均質(zhì)含沙水流的水頭損失。

2.2 含沙水錘求解方法關(guān)于清水水錘的特征線法與管路水力元件及邊界條件的數(shù)學(xué)描述在文獻(xiàn)[4]中有詳細(xì)描述，本文不再贅述。若想采用特征線法求解含沙水錘方程，需做出以下假設(shè)：

(1)假設(shè)管路截面面積保持不變，忽略泥沙淤積對于管路截面形狀的影響；

(2)假設(shè)水錘過程中壓力波速保持不變，管路截面壓力不受固相顆粒分布的影響；

(3)忽略泥沙對泵、閥等水力元件性能的影響，假設(shè)管路水力元件的數(shù)學(xué)描述與文獻(xiàn)[4]一致；

(4)忽略瞬態(tài)摩阻變化對含沙水錘過程的影響，管路糙率按穩(wěn)態(tài)流動(dòng)的水力損失計(jì)算。

基于上述假設(shè)，借鑒韓文亮的非均質(zhì)流漿體水錘壓力波速公式[8]計(jì)算含沙水錘壓力波速am的數(shù)值，其表達(dá)式如下：

(7)

式中：Ew為清水的彈性系數(shù)；Es為固相顆粒的彈性模量；e為管路壁厚；E為管路材料的彈性模量。結(jié)合式(1)、式(7)與文獻(xiàn)[4]中簡化后的用于清水水錘計(jì)算的運(yùn)動(dòng)方程和連續(xù)性方程，可以得出含沙水錘的運(yùn)動(dòng)方程L1和連續(xù)方程L2表達(dá)式如下：

(8)

(9)

由于非恒定摩阻在清水水錘過程中對水錘壓力第一個(gè)周期的影響較小，本文忽略瞬態(tài)摩阻變化對含沙水錘過程的影響，定義當(dāng)量摩阻fm表征非均質(zhì)流含沙水管路穩(wěn)態(tài)流動(dòng)的水頭損失：

(10)

式中當(dāng)量摩阻fm可由穩(wěn)態(tài)流動(dòng)初始條件的f0、v0計(jì)算得出。將當(dāng)量摩阻fm代入方程(8)中得到改進(jìn)后的含沙水錘的運(yùn)動(dòng)方程L1如下：

(11)

通過引入未知因子λ關(guān)聯(lián)方程(11)與方程(9)構(gòu)建線性組合公式L=L1+λL2=0，求解得到λ的兩組特征根代入線性組合公式中得到兩組方程分別用C+和C-表示，其具體表達(dá)式如下：

(12)

(13)

引入管路截面面積A與含沙水流量Q表征混合物流速v，對特征線方程進(jìn)行積分，令v=Q/A，B=am/gA，R=fmΔx/2gDA2得到差分形式的代數(shù)方程如下：

C+：HP=HA-B(QP-QA)-RQA|QA|=0

(14)

C-：HP=HB+B(QP-QB)+RQB|QB|=0

(15)

在方程(14)和方程(15)的基礎(chǔ)上引入邊界條件與水力元件代數(shù)方程，即可求解得到特定時(shí)間位置參數(shù)P點(diǎn)的壓力HPi與流量QPi的數(shù)值解。

3 理論分析

3.1 清水模型驗(yàn)證現(xiàn)以Bergant等[22]RPV系統(tǒng)(水箱-管路-閥門)的清水關(guān)閥水錘實(shí)驗(yàn)數(shù)據(jù)為參考，采用前節(jié)所述含沙水錘數(shù)值計(jì)算特征線方法進(jìn)行關(guān)閥瞬變流分析。該實(shí)驗(yàn)包含高為32 m的上游水箱，通過長度為37.2 m、直徑為0.022 m的管道連接至末段閥門，管道清水條件下的糙率n為0.037。在初始流速v0為0.3 m/s、關(guān)閥時(shí)間為0.009 s的實(shí)驗(yàn)條件下，測量得到了管道末端與中點(diǎn)的清水水錘瞬態(tài)壓力數(shù)據(jù)。為了更好地對比清水水錘與含沙水錘的數(shù)值計(jì)算結(jié)果，針對ρm=2650 kg/m3的泥沙顆粒，在顆粒直徑ds=0.01～0.5 mm、含沙量S=1～100 kg/m3的范圍內(nèi)，選取三組不同參數(shù)介質(zhì)(其中一組為清水，即含沙水顆粒參數(shù)設(shè)定為S=0 kg/m3，ds=0 mm)進(jìn)行對比分析如圖1所示。此外，由于含沙水與清水的密度不同，本文將含沙水壓力水頭Hm統(tǒng)一換算至清水壓力水頭Hw，其表達(dá)式如下：

圖1 含沙水錘數(shù)值計(jì)算瞬態(tài)壓力與清水水錘實(shí)驗(yàn)壓力對比

(16)

從圖中可以看出，清水水錘壓力的瞬態(tài)數(shù)值計(jì)算結(jié)果與Bergant實(shí)驗(yàn)數(shù)據(jù)在水錘壓力波的第一個(gè)周期內(nèi)吻合程度較好，其最大壓力值為72.32 m與實(shí)驗(yàn)值偏差小于0.5%；而第一個(gè)周期后仿真數(shù)據(jù)與實(shí)驗(yàn)結(jié)果的偏差較大，這種壓力波衰減速率和波速變化的不同主要是受到實(shí)驗(yàn)誤差與瞬態(tài)摩阻效應(yīng)影響。對于另外兩組泥沙而言，S=20 kg/m3，ds=0.05 mm時(shí)計(jì)算得到的最大瞬態(tài)壓力為73.45 m，較清水水錘高出1.55%；而S=100 kg/m3，ds=0.5 mm時(shí)最大壓力值為77.91 m，較清水高出7.71%。對于最小壓力值而言，含沙水錘與清水水錘數(shù)值計(jì)算結(jié)果的偏差量小于1.62 m，較最大壓力而言偏差較小。無論是管道末端還是中點(diǎn)的含沙水錘壓力波動(dòng)規(guī)律與清水仿真一致，其波速存在一定的偏差但差別較小?？梢钥闯?，低含沙量、低顆粒直徑的含沙水錘壓力與清水水錘更為接近，而與高含沙量、高顆粒直徑的偏差較大。這一方面論證了含沙水錘數(shù)值計(jì)算結(jié)果隨顆粒濃度和直徑變化趨勢的合理性，另一方面說明了一定規(guī)模的顆粒存在對于含沙管路水錘壓力產(chǎn)生一定的影響，需要進(jìn)一步具體地分析不同的顆粒參數(shù)對含沙水錘過程壓力波動(dòng)的影響。

圖2 不同顆粒參數(shù)與波速、糙率、最大壓力關(guān)系對比

3.3 流動(dòng)參數(shù)分析計(jì)算得到不同初始流速條件下相同顆粒參數(shù)(S=20 kg/m3，ds=0.05 mm)的瞬態(tài)壓力變化如圖3所示。結(jié)果表明不同流速條件下的壓力增量不同而波速相同，這與文獻(xiàn)[4]中清水水錘的規(guī)律一致；此外，不同流速條件下的壓力波衰減過程不同，需要進(jìn)一步分析相關(guān)因素糙率的變化規(guī)律。圖4的計(jì)算結(jié)果表明，含沙水水錘過程的壓力增量與初始流速線性相關(guān)，這與清水水錘Joukovsky方程[4]描述的規(guī)律是一致的。然而，本文給定管路條件的清水波速為1319 m/s，給定顆粒參數(shù)下的含沙水錘波速為1325.5 m/s，其最大壓力計(jì)算結(jié)果與壓力增量變化率存在明顯偏差。進(jìn)一步分析不同初始流速條件下的糙率變化可以看出，當(dāng)流速小于0.05 m/s時(shí)糙率與流速呈現(xiàn)出明顯負(fù)相關(guān)；而當(dāng)流速大于0.05 m/s時(shí)逐漸趨于0.0396不變，這與給定的管道清水條件糙率0.037增大約7%。

圖3 不同初始流速瞬態(tài)壓力波動(dòng)對比

圖4 不同管道流速最大壓力與糙率對比

綜上所述可以得出含沙量、顆粒直徑與初始流速共同影響非均質(zhì)流含沙水穩(wěn)態(tài)流動(dòng)的管路糙率，而初始流速與含沙量同時(shí)也會(huì)對水錘過程的壓力增量產(chǎn)生影響，且不同顆粒參數(shù)與流動(dòng)參數(shù)對水錘過程壓力波動(dòng)影響的程度和相關(guān)性不一。

4 結(jié)論

本文基于含沙水費(fèi)祥俊阻力模型、韓文亮非均質(zhì)流漿體水錘方程和特征線方法，定義了恒定當(dāng)量摩阻，構(gòu)建了含沙水錘的特征線方程，建立了RPV系統(tǒng)關(guān)閥的數(shù)值計(jì)算模型，并結(jié)合清水實(shí)驗(yàn)驗(yàn)證。結(jié)果表明：不同顆粒參數(shù)含沙水錘計(jì)算結(jié)果最大壓力差值明顯，而最小壓力的差別相對較??；含沙量與顆粒直徑越小，含沙水錘瞬態(tài)壓力計(jì)算結(jié)果與清水實(shí)驗(yàn)數(shù)據(jù)越接近；隨著含沙量的增加，含沙水錘的壓力增量增大，含沙量100 kg/m3時(shí)的最大壓力較清水水錘高出7.71%，但顆粒直徑對壓力增量的影響較?。划?dāng)量摩阻換算得到的管路糙率與含沙量和顆粒直徑呈正相關(guān)，與初始流速呈負(fù)相關(guān)，給定顆粒參數(shù)范圍內(nèi)糙率變化超過4.7%。

含沙水錘過程中顆粒對摩阻影響的作用機(jī)理十分復(fù)雜，本文采用恒定流當(dāng)量摩阻代換，存在一定的誤差。此外，對于長期運(yùn)行的供水管路而言，顆粒沉積甚至堵塞現(xiàn)象十分嚴(yán)重，非均質(zhì)流固液混合物的水錘壓力波速公式往往不能適用，還需要研究顆粒沉積對波速的影響。上述兩方面問題將結(jié)合模型實(shí)驗(yàn)開展進(jìn)一步研究。