事故停泵水錘計算模型參數(shù)局部敏感性分析

(山東大學(xué) 土建與水利學(xué)院，山東 濟南 250061)

1 研究背景

隨著我國社會經(jīng)濟的逐步發(fā)展，水資源供需日益緊張。為緩解供需之間的矛盾，大型泵站管道調(diào)水工程修建得越來越多，因此開展保障輸水工程安全運行方面的研究意義重大。水錘的防護(hù)是泵站設(shè)計過程中的重點，國內(nèi)外學(xué)者對此開展了大量的研究[1-3]?；谒N計算模型并結(jié)合數(shù)值模擬的方法，為泵站設(shè)計時如何采取安全經(jīng)濟的水錘防護(hù)措施提供了有效手段[4-6]。然而，由于模型參數(shù)之間的非線性關(guān)系及模型自身的不確定性，給模型的參數(shù)識別帶來了一定的困難，為選取最優(yōu)防護(hù)條件來調(diào)整模型參數(shù)以達(dá)到預(yù)期的防護(hù)效果帶來了阻礙。如何根據(jù)輸水要求來構(gòu)建輸水模型，并高效地實現(xiàn)安全和經(jīng)濟的防護(hù)措施一直是研究的熱點問題。因此，開展對水錘計算模型敏感性的研究，可以加強對模型的認(rèn)識，為水錘防護(hù)情景模擬時高效率地調(diào)整控制參數(shù)以快速達(dá)到預(yù)期防護(hù)效果提供參考。

模型敏感性分析是研究模型輸入?yún)?shù)對輸出結(jié)果影響程度大小的一種方法，它可以對影響模型的關(guān)鍵輸入?yún)?shù)進(jìn)行有效識別,具體可分為局部敏感性分析和全局敏感性分析兩大類[7]?，F(xiàn)有的研究方法主要包括：Sobol法[8]、Morris篩選法[9-10]、傅立葉幅度敏感性分析方法(FAST，及EFAST)[11]等。黃金良[12]、王浩昌[13]、趙磊[10]等分別利用SWMM模型，針對降雨徑流或水質(zhì)的局部或全局的敏感性進(jìn)行了分析。黃清華等[14]采用擾動分析的方法，針對SWAT模型徑流參數(shù)敏感性開展了研究。宋丹明等[15]，將Morris篩選法和EFAST運用于CERES-Wheat模型中進(jìn)行敏感性分析，結(jié)果表明Morris篩選法在該模型中分析效果更好。Xu X等[16]基于LH-OAT的方法，針對耦合模型SWAP-EPIC進(jìn)行了參數(shù)的全局敏感性分析，為提高模型的精度提供了有效的參考途徑。

在水錘計算模型方面，黃凱等[17]針對核電廠用水系統(tǒng)事故停泵工況，分析了管網(wǎng)高程差、水泵轉(zhuǎn)動慣量、閥門關(guān)閉時間變化時水錘峰值(最大水錘壓力)的對應(yīng)變化趨勢，但是并沒有結(jié)合敏感性分析方法開展更進(jìn)一步的研究。Chaudhry M H[18]提出了在變量預(yù)期的范圍內(nèi)，采用一定幅度的擾動進(jìn)行分析，不失為一種分析瞬變流水力參數(shù)敏感性的有效方法。

基于以上研究，本文從水錘運動的基本微分方程——連續(xù)性方程與運動性方程出發(fā)，結(jié)合現(xiàn)有的理論分析，采用修正的Morris篩選法，對泵站事故停泵工況條件開展了水錘計算模型的參數(shù)敏感性分析。

2 水錘計算基本微分方程

2.1 水錘計算基本方程

對于管道水錘計算，以管道水流的基本微分方程——連續(xù)性方程和動量方程來進(jìn)行表達(dá)[6，19-20]。

連續(xù)性方程為

(1)

動量方程為

(2)

式中,H為管道沿線水頭,m；V為平均流速，m/s；g為重力加速度,m/s2；f為摩阻系數(shù);α為管軸中心線與水平線的夾角；D為管道直徑,mm；a為水擊波波速；x為空間步長；t為時間步長。

聯(lián)立式(1)與式(2)，在圖1 所示的網(wǎng)格上采用特征線法，可以求解出管道內(nèi)部任意一點的壓力值，如式(3)所示，具體求解過程參見EB.懷特[19]等人的研究。

圖1 特征線網(wǎng)格Fig.1 The grid of characteristic line

(3)

式中，B=a/gA；CM,CP均為計算過程量；HR,HS分別為R，S兩點的壓力；QR,QS分別為R，S兩點的流量。

2.2 模型關(guān)鍵參數(shù)選取

上述2.1節(jié)中對水錘計算基本微分方程和解進(jìn)行了簡要介紹。根據(jù)管道內(nèi)部水錘壓力計算式求解結(jié)果(即式(3))不難看出，管道內(nèi)部流速(流量)、管徑、摩阻系數(shù)、水擊波速、管道坡度(單位高差內(nèi)管道的長度)等參數(shù)與事故停泵時沿程管路的水錘計算結(jié)果具有一定的關(guān)系。同時，水錘計算結(jié)果與事故停泵時水泵及電機機組的轉(zhuǎn)動慣量[21]、泵后閥關(guān)閉時間相關(guān)[6，19-20]

根據(jù)以上分析，選取了水擊波波速、管道摩阻系數(shù)、機組轉(zhuǎn)動慣量、干管輸水流量、末端閥關(guān)閥時間、輸水管道長度、管道直徑及最大水錘壓力共8個參數(shù)為研究對象，進(jìn)行水錘計算模型事故停泵工況的敏感性分析。其中，事故停泵時產(chǎn)生的最大水錘壓力為研究的重要指標(biāo)，也是本次模型研究的輸出參數(shù)。

3 情景設(shè)定

借助于管道高程、長度、機組、防護(hù)裝置等參數(shù)來構(gòu)建輸水管路的模型，并確定了初始運行時的參數(shù)。初始條件下構(gòu)建的模型管道長為5.4 km，輸水干管管徑為DN1600，泵房內(nèi)設(shè)置6臺機組，單臺機組的額定輸水流量為0.694 m3/s，設(shè)計揚程為60 m。泵站平面布置情況見圖2。穩(wěn)態(tài)運行時進(jìn)出水位分別為85.0，125.0 m，5臺機組穩(wěn)態(tài)時的輸水流量為3.5 m3/s。

圖2 泵站布置示意Fig.2 Layout of pump station

表1中的參數(shù)為運行時初始設(shè)定值。此時，產(chǎn)生的最大水錘壓力值為1.122 MPa。將泵站穩(wěn)態(tài)運行時的值設(shè)為初始值，采用擾動分析的方法，對單個參數(shù)在一定范圍內(nèi)變化時沿程管路最大水錘壓力產(chǎn)生的影響進(jìn)行了分析研究。

表1 水錘計算模型主要參數(shù)基準(zhǔn)值及選取依據(jù)Tab.1 The input parameters of water hammer model

4 水錘計算模型局部參數(shù)敏感性分析

4.1 Morris篩選法

模型參數(shù)敏感性分析包括全局敏感性分析和局部敏感性分析。本文采用使用較為廣泛的Morris篩選法對模型參數(shù)局部敏感性進(jìn)行分析。Morris篩選法的基本思想是評估單個因素產(chǎn)生的微量變化對輸出結(jié)果的響應(yīng)變化。運用Morris篩選法選取模型中的一個變量xi，其余參數(shù)值固定不變，在變量閾值范圍內(nèi)隨機改變xi，通過運行模型即可得到目標(biāo)函數(shù)y(x)=y(x1,x2,x3,…,xn)的值，用影響值ei判斷參數(shù)變化對輸出值的影響程度。其表達(dá)式如下[10，12，22]：

(4)

式中，y*為模型參數(shù)變化后的計算輸出值；y為模型參數(shù)變化前的計算輸出值；Δi為模型參數(shù)的變化值。

修正的Morris篩選法對因子采用初級影響的概念。假設(shè)x1,x2,x3,…,xn是影響模型輸出結(jié)果的n個輸入量，y為模型模擬輸出的結(jié)果，那么認(rèn)為其第i個因子的初級影響值可表示為

(5)

式中，Δ為第i個輸入?yún)?shù)的改變量。

修正的Morris篩選法采用自變量以固定步長變化，靈敏度判別因子取Morris多個平均值，如下式[10,12]：

(6)

式中，S為敏感性判別參數(shù)；Yi為計算模型第i次運行輸出結(jié)果；Yi+1為計算模型第i+1次運行輸出結(jié)果；Y0為模型參數(shù)率定以后計算結(jié)果的初始值(或基準(zhǔn)值)；Pi為第i次模型運算參數(shù)值相對于率定參數(shù)以后參數(shù)值(或基準(zhǔn)值)的變化的百分率；Pi+1為第i+1次模型運算參數(shù)值相對于率定以后初始參數(shù)值(或基準(zhǔn)值)的變化百分率；n為模型運行次數(shù)。

參考趙磊[10]、黃金良[12]等采用的敏感程度劃分的參數(shù)范圍，對敏感性參數(shù)的級別進(jìn)行評價。評價結(jié)果如表2所示。

表2 敏感性層次的劃分Tab.2 The classification criteria for the parameter sensitivity

4.2 模型參數(shù)敏感性分析與討論

在模型輸入?yún)?shù)基準(zhǔn)值既定的基礎(chǔ)上，采用修正的Morris篩選法對水錘計算模型進(jìn)行局部敏感性分析。首先，對輸出結(jié)果進(jìn)行定量的表達(dá)，以固定步長5%一次針對單個變量進(jìn)行擾動分析(擾動范圍在±25%變化幅度內(nèi))。對模型進(jìn)行模擬分析時，單個參數(shù)變化，其他參數(shù)保持不變，基于擾動分析后的各個參數(shù)變化范圍如表3所示。

在各個參數(shù)逐個變化時輸入水錘計算模型，運算并記錄各參數(shù)不同變化幅度時對應(yīng)輸出的最大水錘壓力值(共7×10=70組參數(shù)輸出結(jié)果)，按照式(5)和式(6)先后進(jìn)行計算，從而得到敏感性參數(shù)S的數(shù)值。結(jié)合表2的敏感性劃分方法，可以得到水錘計算模型參數(shù)局部敏感性分析結(jié)果，如表4所示。由圖3可以較為直觀地看出各個參數(shù)對于最大水錘壓力的敏感性程度。

表3 水錘計算模型主要參數(shù)及其取值范圍Tab.3 The main parameters and range of water hammer model

注：參數(shù)6的改變是指在保持各段管道兩端節(jié)點高程不變的情況下，只變化管道的長度，因此其實質(zhì)是改變了管道的坡度。

表4 水錘計算模型參數(shù)局部敏感性分析結(jié)果Tab.4 The sensitivity analysis results of water hammer model parameters

圖3 水錘計算模型參數(shù)局部敏感性分析Fig.3 The sensitivity analysis results of water hammer model parameters

根據(jù)表4和圖3不難看出，在本文的計算中，關(guān)閥時間、管道直徑、管道輸水流量對于水錘計算模型而言是較為敏感的3個參數(shù)，而且敏感性遠(yuǎn)高于剩余的4個參數(shù)。其中，泵后閥關(guān)閥時間、管道直徑為高靈敏度參數(shù)，泵后閥關(guān)閥時間為最敏感的參數(shù)。

泵后閥的關(guān)閥時間快慢，是產(chǎn)生直接水擊和間接水擊的關(guān)鍵。在本研究計算中，水擊波傳遞的一個相長約為9.8 s，而以上的70組數(shù)據(jù)中，部分計算組別在劃分的泵后閥關(guān)閥時間(9.0～15.0 s)內(nèi)，會產(chǎn)生直接的水擊壓強，進(jìn)而使得模型輸出的最大水錘壓力值變化范圍大，增強了該參數(shù)在本模型中的敏感性。

式(1)與式(2)表明，各水力參數(shù)呈現(xiàn)為復(fù)雜的非線性關(guān)系，難以直接對結(jié)果進(jìn)行表征和分析。如果簡單地從水擊壓強計算的基本公式(Δh=-aΔv/g，Δv為流速的變化量)來分析，那么管道的流速變化值增量與水錘壓力呈正相關(guān)關(guān)系，而管道流量、直徑的變化，會直接影響到管道穩(wěn)態(tài)運行時流速的大小，進(jìn)而引起事故停泵時管路內(nèi)部流速變化量的顯著改變以致增大了水錘壓力；況且在工程實踐經(jīng)驗中，認(rèn)為減小管道流速是一種預(yù)防水擊危害的有效方法。單從對這兩點簡要分析的情況來看，管徑和輸水流量作為2個較為敏感的參數(shù)具有一定的理論分析基礎(chǔ)。

5 結(jié) 論

本文基于水錘計算微分方程求解結(jié)果和現(xiàn)有的研究理論，采用修正的Morris篩選法，對事故停泵工況條件下的水錘計算模型進(jìn)行了參數(shù)局部敏感性分析，分析結(jié)果表明：

(1) 模型選取的7個輸入?yún)?shù)靈敏度依次為泵后閥關(guān)閥時間、管道管徑、輸水流量、摩阻系數(shù)、水擊波波速、機組轉(zhuǎn)動慣量、管道長度(坡度)。

(2) Morris篩選法在水錘計算模型中不失為一種較為方便實用的局部敏感性分析方法，可以加強對模型參數(shù)的識別；而且在進(jìn)一步開展水錘防護(hù)情景模擬時，其高效率的控制參數(shù)調(diào)整可以快速達(dá)到預(yù)期的防護(hù)效果。