周領(lǐng)，劉靜，歐傳奇，錢蘇平，盧坤銘

(1. 河海大學(xué)水利水電學(xué)院，江蘇 南京 210098； 2. 國際小水電中心，浙江 杭州 310000； 3. 泰州市水利局，江蘇 泰州 225300； 4. 中國三峽建工(集團(tuán))有限公司，四川 成都 610041)

水錘現(xiàn)象是有壓輸水管道系統(tǒng)運(yùn)行中常見的問題之一.危險(xiǎn)的水錘壓力會影響管道系統(tǒng)正常運(yùn)行，甚至可能導(dǎo)致管道發(fā)生破壞性事故[1-3].因此，針對管道系統(tǒng)內(nèi)可能發(fā)生的水力瞬變問題進(jìn)行精準(zhǔn)的預(yù)測，具有十分重要的意義.

傳統(tǒng)的水錘計(jì)算方法通常采用恒定摩阻或準(zhǔn)恒定摩阻模型，主要可用于計(jì)算輸水管道系統(tǒng)發(fā)生水錘時產(chǎn)生的正、負(fù)壓波的最大峰值，但不能精確地預(yù)測壓力波的衰減過程.為提高計(jì)算精度，對傳統(tǒng)的水錘模型摩阻積分項(xiàng)采用不同階的精度[4]，但是通過試驗(yàn)和理論研究表明，對于高黏度和高速的管道瞬變流，傳統(tǒng)的水錘模型預(yù)測結(jié)果誤差仍較大，主要表現(xiàn)在壓力波的衰減和相位偏差上[5].

傳統(tǒng)的水錘模型通常假設(shè)恒定流摩阻可描述瞬變流摩阻，故模擬結(jié)果與實(shí)際偏差較大.為了能精確地模擬管道內(nèi)水力瞬變過程，在水錘模擬計(jì)算中引入動態(tài)摩阻模型.其中，應(yīng)用最為廣泛的是以BRUNONE等[6]提出的模型為代表的瞬時加速度模型(也稱經(jīng)驗(yàn)修正模型)，及以ZIELKE[7]為代表提出的加權(quán)類模型.BRUNONE類動態(tài)摩阻模型適用于管道紊流瞬變流，模型結(jié)構(gòu)簡單，易于實(shí)現(xiàn)高階精度及其他格式的求解，且模擬精確度較高.加權(quán)類動態(tài)摩阻模型被公認(rèn)為計(jì)算頻率相關(guān)摩阻損失的精確模型[8].ZIELKE模型能夠準(zhǔn)確模擬管道層流和低雷諾數(shù)紊流的瞬變過程.但由于加權(quán)函數(shù)的復(fù)雜性，ZIELKE模型計(jì)算所需存儲空間較大，計(jì)算機(jī)運(yùn)行時間較長，模型計(jì)算效率較低，故該模型只能適用于非常短的瞬變過程模擬.為此，TRIKHA[9]提出用3個指數(shù)函數(shù)之和型式的函數(shù)近似代替ZIELKE模型的加權(quán)函數(shù)，得到的模型計(jì)算結(jié)果在一定范圍內(nèi)與ZIELKE模型結(jié)果相差不大，但是計(jì)算效率得到了較大的提高.隨后，較多學(xué)者也對ZIELKE模型加權(quán)函數(shù)進(jìn)行了相同型式的改進(jìn)[10-14].另外，VARDY等[15-17]通過在剪切層中使用線性變化的凍結(jié)紊流黏度，提出了一種新型的加權(quán)函數(shù)計(jì)算公式.

為進(jìn)一步提高加權(quán)類動態(tài)摩阻模型的計(jì)算效率，擴(kuò)大模型的適用范圍，文中在TRIKHA[9]基礎(chǔ)上嘗試采用不同的擬合方式對ZIELKE加權(quán)函數(shù)進(jìn)行相同函數(shù)型式的擬合，建立一種新的近似加權(quán)函數(shù).為驗(yàn)證新函數(shù)的優(yōu)越性，比較其與ZIELKE加權(quán)函數(shù)及前人提出的改進(jìn)近似加權(quán)函數(shù)的W-τ關(guān)系曲線，利用試驗(yàn)數(shù)據(jù)對比模型計(jì)算結(jié)果，并在相同的運(yùn)行條件下，對比各模型計(jì)算運(yùn)行時長.

1 加權(quán)類動態(tài)摩阻模型

1.1 基本控制方程

一維水錘模型的基本控制性方程

(1)

(2)

式中：H為水頭，m；a為波速，m/s；g為重力加速度，m/s2；v為斷面平均流速，m/s；J為管道內(nèi)單位長度的水頭損失，J=JS+JU，其中JS和JU分別為由恒定摩擦和非恒定摩擦引起的水頭損失.

利用特征線法將標(biāo)準(zhǔn)線性雙曲型微分方程(1)，(2)變換為常微分方程，每個方程都有1個約束.

C+：

(3)

(4)

C-：

(5)

(6)

當(dāng)采用一階有限差分近似時，將式(3)，(5)分別沿著式(4)，(6)特征線進(jìn)行積分，得到控制性方程，即

C+：

(7)

C-：

(8)

式中：下標(biāo)P表示待求變量；下標(biāo)A和B分別表示待求點(diǎn)的前、后計(jì)算節(jié)點(diǎn)處上一時刻的已知量；Δt為時間步長，Δt=Δx/a，Δx為特征線網(wǎng)格長度.

1.2 ZIELKE模型

1968年ZIELKE提出了一種適用于層流和雷諾數(shù)小于108湍流的非恒定摩阻模型，認(rèn)為管道瞬變流中瞬時管壁切應(yīng)力由恒定項(xiàng)和附加項(xiàng)組成，恒定項(xiàng)為恒定狀態(tài)切應(yīng)力值，附加項(xiàng)利用加權(quán)函數(shù)考慮了歷時速度和加速度對當(dāng)前流態(tài)的影響.

非恒定摩阻項(xiàng)計(jì)算式為

(9)

W計(jì)算式為

當(dāng)τ≥0.02時

(10)

當(dāng)τ<0.02時

(11)

式中：mi，ni為系數(shù)，mi={26.374 4; 70.849 3; 135.019 8; 218.921 6; 322.554 4};ni={0.282 095; -1.25; 1.057 855; 0.937 5; 0.396 696; -0.351 563}.

1.3 TRIKHA模型

表1為各模型系數(shù)mi，ni.

表1 各模型系數(shù)mi，ni

1975年，TRIKHA對ZIELKE模型加權(quán)函數(shù)進(jìn)行簡化，提出一種近似加權(quán)函數(shù)Wapp(τ)代替ZIELKE加權(quán)函數(shù)W(τ)，新函數(shù)消除了原函數(shù)在計(jì)算中對大量過去速度變化值的需求，減小了模型計(jì)算所需存儲空間，縮短了計(jì)算機(jī)運(yùn)算時長，故在一定程度上提高了模型的計(jì)算效率.近似加權(quán)函數(shù)Wapp(τ)的函數(shù)型式為

(12)

為了使近似加權(quán)函數(shù)值Wapp(τ)值和ZIELKE實(shí)際加權(quán)函數(shù)值W(τ)值接近，令Wapp(τ)滿足以下3個條件，即

2) 當(dāng)τ=0.000 1,0.001,0.01時，

Wapp(τ)=W(τ).

3) 當(dāng)τ>0.1時,

求得

mi={40.0; 8.1； 1}；ni={8 000; 200； 26.4}.

TRIKHA改進(jìn)模型能夠較準(zhǔn)確地模擬一定頻率范圍內(nèi)管道層流水力瞬變過程中的水錘壓力、壓力波變形以及壓力衰減，并且較大程度地縮短了模型計(jì)算時長，提高了模型的計(jì)算效率.

之后，較多學(xué)者[10-14]對ZIELKE加權(quán)函數(shù)進(jìn)行了相似函數(shù)型式的近似代換，分別采用不同數(shù)量的指數(shù)函數(shù)，對比TRIKHA模型，在不同程度上提高了模型的計(jì)算準(zhǔn)確性，各改進(jìn)加權(quán)函數(shù)系數(shù)mi，ni見表1.

2 改進(jìn)近似加權(quán)函數(shù)

TRIKHA對ZIELKE動態(tài)摩阻模型的加權(quán)函數(shù)進(jìn)行了3個指數(shù)函數(shù)之和型式的近似替換，提高了模型計(jì)算效率.但改進(jìn)的加權(quán)函數(shù)在τ<0.000 05時與原始加權(quán)函數(shù)擬合誤差較大，故取計(jì)算網(wǎng)格較多(即時間步長較小)時，TRIKHA模型計(jì)算偏差較大.之后的學(xué)者[10-14]均是針對減小加權(quán)函數(shù)擬合誤差進(jìn)行研究改進(jìn)，提高了模擬精度的同時降低了模型的計(jì)算效率.文中在TRIKHA模型基礎(chǔ)上，采用不同的擬合方式對ZIELKE加權(quán)函數(shù)進(jìn)行新的改進(jìn)，在保證模型計(jì)算效率的同時提高其適用范圍及計(jì)算準(zhǔn)確性.

2.1 改進(jìn)近似加權(quán)函數(shù)

對ZIELKE加權(quán)函數(shù)在τ=1.0×10-9～0.02的函數(shù)進(jìn)行新函數(shù)的近似替換，τ>0.02部分仍采用原ZIELKE加權(quán)函數(shù).利用MATLAB中的LSQNONLIN函數(shù)模塊對W-τ樣本點(diǎn)曲線進(jìn)行以下指定函數(shù)型式的擬合，該函數(shù)模塊在擬合過程中結(jié)合了遠(yuǎn)離最小值的線性逼近和接近最小值的平方逼近，專門用于解決最小二乘方法的問題，可使得擬合相對誤差最小.指定函數(shù)型式為

(13)

利用上述擬合方式得到指定函數(shù)各系數(shù)值為

mi={8 894.82; 21.51; 4.18};

ni={127 885.9; 1 803.22; 73.07}.

2.2 加權(quán)函數(shù)對比驗(yàn)證

為驗(yàn)證新的改進(jìn)近似加權(quán)函數(shù)的精確性，將新函數(shù)Wnew與ZIELKE加權(quán)函數(shù)及已有近似加權(quán)函數(shù)的W-τ關(guān)系曲線(采用雙對數(shù)坐標(biāo))進(jìn)行對比，如圖1所示.

圖1 近似加權(quán)函數(shù)曲線對比圖

由圖1可見，新的近似加權(quán)函數(shù)曲線與ZIELKE加權(quán)函數(shù)曲線，在τ為1.0×10-9～1.0×10-4內(nèi)，擬合準(zhǔn)確度相比于TRIKHA函數(shù)高，在τ為1.0×10-4～0.02內(nèi)，擬合度相差不大.在τ為1.0×10-9～1.0×10-5內(nèi)，新函數(shù)比SCHOHL，KAGAWA，VITKOVSKY及VARDY & BROWN近似加權(quán)函數(shù)曲線擬合度更高；當(dāng)τ越小時，新的近似加權(quán)函數(shù)相比于相同系數(shù)數(shù)量下的TRIKHA函數(shù)擬合優(yōu)越性越明顯，故新函數(shù)代入模型的計(jì)算結(jié)果應(yīng)更接近ZIELKE模型計(jì)算結(jié)果.

3 模型驗(yàn)證

利用HOLMBOE等[18]及BERGANT等[19]2個試驗(yàn)案例數(shù)據(jù)進(jìn)行模擬計(jì)算，為驗(yàn)證當(dāng)計(jì)算所取量綱一的時間間隔Δτ<1.0×10-4時，文中提出的新模型計(jì)算結(jié)果比TRIKHA更接近ZIELKE計(jì)算結(jié)果，本次的模擬計(jì)算均取Δτ=1.0×10-6.將新模型模擬計(jì)算結(jié)果與ZIELKE模型、TRIKHA模型計(jì)算結(jié)果及試驗(yàn)數(shù)據(jù)進(jìn)行對比，驗(yàn)證新模型計(jì)算的準(zhǔn)確性；并在相同的模擬計(jì)算條件下，比較新模型與各改進(jìn)模型及ZIELKE模型的運(yùn)行時長，驗(yàn)證新模型計(jì)算的高效性.

3.1 HOLMBOE等試驗(yàn)驗(yàn)證

HOLMBOE等[18]在半徑R=0.012 7 m、長度L=36.09 m的銅管上進(jìn)行了測試，銅管上游連接到一個由壓縮空氣保持恒定壓力的儲罐上，上游水頭恒定為H=32 m.試驗(yàn)所用液體為黏度υ=039.67×10-6m2/s的油.實(shí)測壓力波速度a=1 324.36 m/s.在流動過程中，下游閥門在管道中迅速關(guān)閉.在管路端點(diǎn)(閥門附近)測量壓力波動.管道內(nèi)初始流態(tài)為層流，流速v0= 0.128 m/s (Re=82).

利用上述試驗(yàn)數(shù)據(jù)進(jìn)行模型的模擬計(jì)算，圖2為文中模型、ZIELKE模型及TRIKHA模型計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)的對比.

圖2 模型計(jì)算結(jié)果對比(HOLMBOE)

圖2a和圖2b分別為試驗(yàn)系統(tǒng)閥門處和管道中點(diǎn)處的水頭隨時間變化曲線，圖中Ut,UH分別為時間和壓力水頭，且量為綱為一，Ut=ta/L,UH=ΔH/(av0/g),其中ΔH為變化水頭與上游恒定水頭之差.圖2中，文中提出的新模型的計(jì)算結(jié)果曲線與TRIKHA模型計(jì)算結(jié)果曲線近似重合，相差較小.但在局部位置(尤其在水頭突變之后位置)，文中模型的計(jì)算結(jié)果曲線明顯更加接近ZIELKE模型計(jì)算結(jié)果曲線.對比分析文中模型和TRIKHA模型峰谷值處的誤差大小(以最后1個時間周期為例，分別計(jì)算兩模型與ZIELKE模型對應(yīng)峰谷值的相對誤差)：管道閥門處，文中模型波峰波谷水頭計(jì)算相對誤差分別為8.56%和7.02%，TRIKHA模型分別為9.72%和8.33%；管道中點(diǎn)處，新模型波峰波谷水頭計(jì)算相對誤差分別為9.98%和10.64%，TRIKHA模型分別為17.85%和18.37%.可見，新模型的峰谷值計(jì)算誤差小于TRIKHA模型.故分析表明，在該試驗(yàn)案例下，文中模型的計(jì)算結(jié)果誤差更小，模型計(jì)算準(zhǔn)確性更高.

3.2 BERGANT等試驗(yàn)驗(yàn)證

BERGANT等[19]試驗(yàn)參數(shù)如下：上游是恒水位壓力罐，水頭H=32 m，下游閥門在0.009 s內(nèi)線性關(guān)閉.管道長L=37.23 m，管道內(nèi)徑D=0.022 1 m，管壁厚1.63 mm，波的傳播速度a=1 319 m/s，運(yùn)動黏滯系數(shù)υ=1.184×10-6m2/s，管道內(nèi)初始流態(tài)為層流，流速v0= 0.1 m/s.圖3為文中模型、ZIELKE模型及TRIKHA模型計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)的對比圖，其中圖3a和圖3b分別為閥門處和管道中點(diǎn)處的水頭隨時間變化的結(jié)果.

圖3 模型計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)對比

由圖3可見，文中提出的新的加權(quán)函數(shù)模型與ZIELKE模型計(jì)算結(jié)果相差較小，TRIKHA模型較新模型在相同對應(yīng)時間周期內(nèi)水頭峰谷處的模擬誤差較大，且誤差隨時間的增長而增大.對比新模型和TRIKHA模型峰谷值誤差大小(以最后1個時間周期為例，分別計(jì)算兩模型與ZIELKE模型對應(yīng)峰谷值的相對誤差)：管道閥門處，新模型波峰波谷水頭計(jì)算相對誤差分別為0.50%和0.98%，TRIKHA模型分別為1.74%和3.15%；管道中點(diǎn)處，新模型波峰波谷水頭計(jì)算相對誤差分別為0.93%和1.58%，TRIKHA模型分別為1.54%和2.64%.結(jié)果表明，新模型的峰谷值計(jì)算誤差小于TRIKHA模型.因此，相比于TRIKHA模型，文中模型計(jì)算結(jié)果更加接近ZIELKE模型計(jì)算結(jié)果，計(jì)算誤差更小.

3.3 模型計(jì)算效率

為驗(yàn)證文中模型的計(jì)算效率，利用各模型對上述BERGANT試驗(yàn)進(jìn)行數(shù)值模擬，控制模型運(yùn)行條件相同，比較文中模型及其他各模型的計(jì)算運(yùn)行時長.表2為模型計(jì)算采取不同時間網(wǎng)格數(shù)時各模型運(yùn)行所需時長.

表2 各模型運(yùn)行時長

由表2可見，當(dāng)模型計(jì)算所取時間網(wǎng)格數(shù)相同(時間步長相同)時，文中模型計(jì)算時長與TRIKHA模型相近，均小于其他各改進(jìn)模型，ZIELKE模型計(jì)算所需時長最大；各模型計(jì)算時長均隨時間網(wǎng)格數(shù)的增加(時間步長的減小)而增大，文中模型與TRIKHA模型增長速度最慢，ZIELKE模型呈指數(shù)型增長，增長速度最快.綜上，新模型的計(jì)算效率對比已有的改進(jìn)模型在一定程度上有所提高.

4 結(jié) 論

文中建立了加權(quán)類動態(tài)摩阻模型新的近似加權(quán)函數(shù)，并通過已有的試驗(yàn)數(shù)據(jù)對新模型的準(zhǔn)確性及高效性進(jìn)行了驗(yàn)證.

1) 對比相同改進(jìn)函數(shù)型式的TRIKHA模型，文中提出的新模型計(jì)算準(zhǔn)確性更高，計(jì)算結(jié)果誤差更小，且誤差隨模型計(jì)算所取時間步長的減小而減小.

2) 對比其他改進(jìn)模型，在相同的運(yùn)行條件下，該模型縮短了模型模擬所需計(jì)算時長，提高了加權(quán)類動態(tài)摩阻模型的計(jì)算效率，且隨著計(jì)算時間步長減小，計(jì)算時長增長最緩，即新模型高效性的優(yōu)勢隨時間步長的減小而更加明顯.

文中模型可廣泛應(yīng)用于有壓管道系統(tǒng)水力瞬變過程的模擬預(yù)測，能夠準(zhǔn)確、高效地對管道內(nèi)的壓力波動過程進(jìn)行數(shù)值計(jì)算，對有壓管道的設(shè)計(jì)和安全運(yùn)行具有重要意義.