萬五一，陳瀟逸，張永進

（1.浙江大學 水利工程學系，浙江 杭州 310058；2.浙江省水利水電勘測設計院，浙江 杭州 310002）

1 研究背景

水錘現象可以導致管道的壓力振蕩，造成爆管和壓扁等危害。在有壓管路系統中，閥門的突然啟閉、水泵的事故停機、水輪機導葉的驟然啟閉等均能引起壓力與流量的急劇變化，從而導致水錘危害的發(fā)生。如何模擬、預測和控制水錘是目前水力學研究領域的熱點問題之一[1-3]。水錘的模擬對防止管道的水錘危害與保障管道的安全運行至關重要。隨著計算機技術的發(fā)展，不同的數值方法被用于輸水管路的水錘模擬，主要有特征線法（MOC）[4]、有限元法[5]、有限體積法[6]等。在此基礎上，國內外的學者對水錘數值計算方法不斷進行改進，也將一些新的數值方法引入到水錘模擬中。耿艷芬等[7]采用特征分解構建了一階和二階Godunov 格式的水錘計算模型，并通過通量限制和重構獲取了高精度的格式。趙修龍等[8]建立了Crank-Nicolson 格式的水錘離散方程，該格式具有二階精度且穩(wěn)定性好。劉韓生等[9]將高分辨TVD 格式引入水錘方程的求解中，該格式具有精度高、耗散低、不產生數值振蕩等優(yōu)點。Kim[10]提出采用脈沖響應法模擬配水系統的瞬變流，并提出了一種考慮層流與湍流過渡的修正方法。Budinski[11]采用了格子玻爾茲曼方法求解管道瞬變流問題，該方法具有不受MOC網格限制與計算效率高的優(yōu)點。

盡管如此，MOC 仍為目前最常用的水錘數值計算方法之一，其具有計算簡潔、易于編程、可處理復雜系統等優(yōu)點，廣泛地運用在水力發(fā)電[12]、油氣輸運[13]、核能發(fā)電[14]等各個領域的瞬變流模擬研究。但MOC 通常需要采用規(guī)則的網格，且空間步長與時間步長受到特征線的嚴格限制。在解決特定的復雜管路系統問題時，需要通過調整波速、改變管長或插值等方法來滿足不同管道的相同時步要求。若系統中連接有一段很短的管道，則所有管道需要選取很短的時間步長以滿足同步，可能需要較大的計算時耗[4]。

MacCormack格式是一種求解可壓縮流體流動問題的有限差分格式，在時間和空間上具有二階精度[15]。該方法具有易于理解，便于實現的優(yōu)點[16]，在計算流體力學中廣泛地應用于Navier Stokes 方程的求解[17-20]。但是MacCormack格式目前在瞬變流的數值模擬中運用的相對較少。Chaudhry 等[5]簡要介紹了采用MacCormack格式對水錘的擬線性雙曲型偏微分方程的求解過程。Amara[21]將MacCormack格式引入摩擦管道的水錘數值計算中，但仍基于特征線法建立了邊界條件，從而無法對網格中的時間步長與空間步長進行單獨設定。王勇等[22]基于MacCormack格式建立了長輸液體管道的水力瞬變模型，并分別采用了MOC方法和線性外插法對邊界條件進行處理。在我們之前的研究中，提出了基于MacCormack時間推進格式的水錘計算模型（MTMS）[23]。與常規(guī)的MOC 相比，該模型突破了網格劃分中時間步長與空間步長相互制約的局限，具備網格劃分靈活的優(yōu)點，但其仍然采用了統一的空間步長。

為了分析包含短管情況下的變特性復雜管路系統中的水錘，拓展變空間步長情況下的同步水錘模擬方法，本文通過MacCormack格式建立了變空間步長水錘求解模型（V-MCS）及相應的邊界條件，對變特性管路系統中的水錘現象進行了數值模擬。在V-MCS方法中采用了變空間步長的不規(guī)則離散網格，該網格不再受特征線的約束，其空間步長的選取只需滿足Courant-Friedrichs-Lowry（CFL）準則的要求，在變特性復雜管路系統中，該方法克服了為保證與短管時間步長一致而需采用較密計算網格的問題。通過規(guī)則網格的MOC方法與MTMS方法對V-MCS方法的模擬結果進行了驗證與分析，本文提出的V-MCS方法對變特性復雜管路系統的水錘數值模擬具有較好的參考價值。

2 數學模型

2.1 基本控制方程描述管道有壓流的基本控制方程包括了連續(xù)方程和運動方程[4]，其矩陣形式可寫為：

式中：h為測壓管水頭；t為時間；v為流速；x為節(jié)點位置；a為水錘波速；g為重力加速度；f為摩阻系數；D為管道直徑。

式（1）包括了一對準線性雙曲線偏微分方程，包含了測壓管水頭h和流速v 兩個因變量，為便于利用MacCormack格式進行偏微分方程離散與因變量求解，將基本控制方程改寫為：

在使用MOC模型求解時，通常將摩阻系數f 近似為恒定摩阻系數來考慮。為提高計算精度，本文采用Brunone模型[24]在控制方程中引入非恒定摩阻。Brunone模型中考慮了瞬時當地加速度和對流加速度對非恒定摩阻的影響，改進后的表達式為[25]：

式中： fq為恒定摩阻系數；k為Brunone 摩擦系數。

采用Vardy的剪切衰減系數C*對Brunone 摩擦系數k進行計算[26]：

式中： Re為雷諾數。結合式（2）和式（3）可得控制方程轉化后的形式為：

2.2 基于變空間步長的MacCormack 格式的水錘模型圖1為變空間步長的x-t 計算網格劃分，其中ns為管道的網格劃分節(jié)點數。在MacCormack格式中，時間步長?t與空間步長?x 不受特征線約束條件?x=a?t的限制，可對網格進行更加靈活的劃分或修改。MacCormack格式的約束條件為CFL（Cou?rant-Friedrichs-Lowry）準則，表達式為[16]：

圖1 變空間步長的計算網格劃分

因此，將網格中的空間步長取為滿足CFL 準則條件下的任意值?xi，i=1，2，…，ns-1。如圖1所示，節(jié)點i-1與節(jié)點i之間的空間步長為?xi-1，節(jié)點i與節(jié)點i+1之間的空間步長為?xi。假設t時刻的節(jié)點均為已知，利用MacCormack格式的預測校正二步法進行時間推進，求解出t+?t時刻各個節(jié)點的和。首先在預測步中對空間導數進行向前差分，將t時刻節(jié)點i的時間導數寫為：

作為校正，結合式（7）和式（10）得到空間導數平均值為：

式（12）為基于MacCormack格式的變空間步長水錘計算模型，適用于網格中t+?t時刻所有中間節(jié)點的求解。在該模型中采用了變空間步長的不規(guī)則離散網格，空間步長可在滿足CFL 準則的前提下進行選取。由于預測校正二步法中，在預測步向前差分，在校正步向后差分，因此該格式具有二階精度。由于該模型中包含了對流加速度項，考慮了流速對水錘波速的影響。此外，該模型采用非恒定摩阻項，有利于更準確地模擬水錘波的衰減過程[27]。

2.3 邊界條件

2.3.1 進口水庫邊界條件 在輸水管路系統中，管路進口端通常與水庫相連，在短時間內，可認為上游水庫水位為恒定。因此管路進口端水庫的限制條件為：

整理得，進口水庫邊界的求解模型為：

2.3.2 調壓井邊界條件 調壓井是輸水管路中重要的水錘防護設備，調壓井邊界條件中的連續(xù)微分方程可表示為：

式中：Qu為調壓井上游流量；Qd為調壓井下游流量；Z為調壓井水位；Aw為調壓井斷面面積。

式中，j、k分別為與調壓井上游和下游相鄰的管道節(jié)點。

對式（17）采用向前差分[28]，并結合式（5）進行整理后得到調壓井邊界的求解模型為：

2.3.3 出口閥門邊界條件 輸水系統出口為控制閥門，通過該閥門可對輸水線路的流量大小進行調節(jié)，閥門的關閉可導致輸水線路中水錘現象的發(fā)生。出口端的可通過線性外插方法計算得出，可根據已計算出的及閥門孔口方程得到，因此出口邊界條件可寫為：

式中Zd為閥門下游水位。當閥門完全關閉后，管道出口端的流速限制條件為：

根據上述提出的基于MacCormack格式的變空間步長水錘模型及邊界條件，從已知時刻t 開始時間推進，可求出t+?t，t+2?t，…，t+n?t 等各個時刻所有節(jié)點的水力參數。

3 模型驗證

為了對本文提出的V-MCS 水錘計算模型及邊界條件進行驗證，分別采用V-MCS方法、MOC方法及MTMS方法對常規(guī)的水庫-管路-閥門（RPV）輸水系統的水錘過程進行模擬。在RPV系統中，管道長度L為1000.0 m，直徑D為0.5m，波速a為1000.0 m/s，糙率系數n為0.012。管道進口端水庫水位Zu為30.0 m，管道的初始流量Q0為0.18m3/s。管道出口端閥門關閉時間Tc為4 s。其中MOC方法采用規(guī)則網格，時間步長?t為0.01 s，空間步長?x為10.0 m，網格分段數為100 段，MTMS方法也采用與MOC方法相同的規(guī)則網格。而V-MCS方法采用變步長網格，時間步長?t為0.01 s，在滿足CFL準則的基礎上將空間步長取為10.0 m≤?x ≤30.0 m 范圍內的隨機數，管道分段數為45 段。為了驗證V-MCS方法的正確性，將其與MOC方法和MTMS方法的水錘壓力波動過程模擬結果進行比較。如圖2與圖3展示了采用三種不同方法計算的RPV系統的水錘壓力波的波動過程，結果表明，三種不同方法計算得到的管道中點與末端的壓力波動曲線基本重合，壓力波的振幅大小吻合良好，驗證了V-MCS方法在水錘計算中的準確性。

圖2 管道中點水錘壓力波動過程模擬結果對比

圖3 管道末端水錘壓力波動過程模擬結果對比

4 案例分析

4.1 初始條件在實際輸水工程項目中，輸水線路通常距離較長，管線根據地形變化而常有起伏，管道的材質也常有差異，因此在運行調度過程中的水力條件相對復雜。水錘數值模擬的目的主要是預測系統中可能出現的壓力與流量波動情況，獲得管道的薄弱節(jié)點及出現的水錘極值壓力，為工程設計與運行提供參考。在案例分析將本文提出的V-MCS方法應用于輸水工程管路系統的水錘模擬中，進一步地對水錘模擬結果進行分析與討論。

如圖4所示為變特性輸水管路系統示意圖，其中上游水庫水位Zu為56.5 m，下游水廠水池的水位Zd為22.0 m。管路總長度為24.16 km，由DN2400 輸水隧洞與DN1800 鋼管組成，沿線設置4 座通氣井與1 座調壓井，管線的基本參數如表1所示。管道末端閥門關閉時間Tc為30 s，管線的流量由4.20 m3/s 減小至0 m3/s，最大計算時刻為10 800 s。本文基于V-MCS方法對該管路系統的水錘過程進行模擬與分析，并采用MOC方法和MTMS方法與其模擬結果進行對比。

圖4 變特性輸水管路系統示意圖

表1 輸水管線基本參數

4.2 模擬與分析在采用MOC方法和MTMS方法進行管道的網格劃分時，先選定時間步長再根據?x=a?t 確定空間步長的大小。為滿足網格的分段數為整數的要求，對波速a進行細微調整，則調整后的波速a′為：

通常情況下，時間步長越小，調整后的波速與實際波速越接近，在本例中，當時間步長取為0.01 s時，偏差最大的波速為1040.9 m/s，偏差為2.7%。MOC方法和MTMS方法中管道的總分段數為2343 段。而采用V-MCS方法進行變網格劃分時，空間步長不再固定，可根據計算需求對空間步長進行選擇性劃分。如表1所示，第6、9、12、14 號管段長度小于200 m，對這些短管的空間步長?x 仍取為a?t，而將管路其余部分的空間步長取為a?t ≤?x ≤12.0 m 范圍內的隨機數，管路總分段數為2217 段。

圖5—圖7為采用V-MCS方法對算例中管路系統的水錘過程進行數值模擬的結果。圖5為基于V-MCS方法得到的閥前水錘壓力波動過程。由圖5可知，由V-MCS方法得到的水錘極值壓力大小與MOC方法、MTMS方法的模擬結果相同。閥前的水錘壓力最大值為75.40 m，水錘壓力最小值為42.25 m。隨著時間推進，由V-MCS方法與MTMS方法模擬得到的水錘壓力波衰減過程與常規(guī)MOC方法相比出現了微小偏差，其原因在于V-MCS方法與MTMS方法中均考慮了非恒定摩阻的影響。圖6所示為閥前流量變化過程，基于V-MCS方法得到的流量變化過程與其余兩種方法得到的模擬結果一致。

圖5 閥前壓力波動過程

圖6 閥前流量變化過程

圖7所示為基于V-MCS方法和MOC方法的管道沿線極值壓力分布情況。圖中藍色線為管道初始恒定水頭壓力線，紅色線為最大壓力包絡線，綠色線為最小壓力包絡線。由圖7可知，管道最小壓力包絡線始終位于管道中心線高程以上，管道全線無負壓現象發(fā)生。管道最大內水壓力為68.90 m，位于管道閥前位置?；赩-MCS方法得到的水錘極值壓力與MOC方法得到的極值壓力包絡線吻合良好。

圖7 管道沿線極值壓力分布

圖8 不同分段數下的閥前壓力波動過程對比

表2 V-MCS與MOC 模擬結果對比

為了分析管道分段數對模擬結果的影響，圖8對不同分段數情況下的閥前壓力波動過程進行了比較。如圖8所示，隨著V-MCS方法空間步長取值范圍的增大，管道分段數隨之減小，位于第一相的水錘極值壓力大小基本一致。但選取的空間步長范圍過大時，對水錘波衰減過程的準確模擬有所影響。將V-MCS方法與MOC方法的計算結果整理于表2，其中無量綱時間Ct為V-MCS方法與MOC方法計算耗時的比值。比較兩種方法的Ct可知，在步長相當的情況下，由于V-MCS的單步計算更為復雜，與MOC方法相比計算耗時更長。由于V-MCS中的空間步長可以進行調整，在保證時間步長不變的情況下，通過增大空間步長可以較大程度減小總的計算耗時。由于在水錘模擬中，線路中的壓力極值是計算中尤為關心的問題，因此在采用V-MCS方法時可適當地減少管道分段數來縮短計算耗時。如需獲得精確的水錘壓力波動衰減過程，則可通過增加管道分段數來提高計算的精度。

4.3 討論模擬結果表明，本文提出的V-MCS方法可有效地得到管道水錘的模擬結果。V-MCS方法具有時間空間上的二階精度，計算結構清晰，易于編程實現。V-MCS方法具有以下優(yōu)點。首先，在模型建立中保留了控制方程中的對流加速度項，將管道流速對水錘波傳播的影響考慮了進來，可為高流速情況下的水錘模擬提供參考。在模型中考慮了非恒定摩阻，能夠較好地反映摩阻對水錘波的衰減過程的影響。其次，在網格劃分上，MOC方法采用了固定的離散網格，且時空步長受到特征線路徑的嚴格約束。但V-MCS方法采用了不規(guī)則的離散網格，根據計算需求調整網格的空間步長大小，突破了時空步長一一對應的限制。因此V-MCS方法可以有針對性地僅對管道薄弱區(qū)域或控制節(jié)點進行詳細劃分，相比MOC方法減少了管道全線的總分段數。此外，對于復雜管網系統，V-MCS方法可以方便地對不同的管道選取不等的空間步長，而無需進行波速調整或插值，從而有利于對復雜系統中的各種管道進行同步求解。但在V-MCS方法中若選取的空間步長過大，對水錘壓力波動衰減過程的準確模擬有所影響。而且，由于V-MCS方法將原有的固定步長根據需要改成了特定的步長，在數值計算中可能一定程度上增加了數值編程的難度。

5 結論

建立了基于MacCormack格式的變空間步長水錘計算模型（V-MCS），在模型中采用了變空間步長的離散網格，空間步長可在滿足CFL 準則的基礎上進行選取。此外，結合線性外插法建立了與模型相匹配的邊界條件。

通過與MOC方法及MTMS方法模擬結果的對比，對提出的V-MCS方法及相應邊界條件進行了驗證。驗證結果表明，V-MCS方法可以有效地對管道水錘過程進行模擬，其模擬結果與MOC方法及MTMS方法的模擬結果吻合良好。

采用V-MCS方法對變特性復雜管路系統的水錘現象進行了模擬與分析，結果表明，V-MCS方法可以準確地模擬出控制節(jié)點的水錘極值壓力大小。由于在模型中考慮了非恒定摩阻的影響，從理論上講，得到的水錘壓力波衰減過程模擬結果與實際更為接近?；赩-MCS方法得到的極值壓力與常規(guī)方法得到的極值壓力包絡線吻合良好。

在包含短管等特定的變特性管路中，V-MCS方法可以根據計算需求對不同管道選取不等的空間步長，減小網格分段數以此節(jié)約計算時間。因此，V-MCS方法有利于對復雜管路系統的水錘求解進行同步運算，對輸水工程管路系統的水錘數值模擬具有一定的實際意義。