王云莉 ，胡 宏

（1.重慶交通大學 西南水運工程科學研究所，重慶 400016；2.內(nèi)河航道整治技術交通行業(yè)重點實驗室，重慶 400016）

科研與管理

有壓管道系統(tǒng)非恒定流模擬研究進展

王云莉1，2，胡 宏1，2

（1.重慶交通大學 西南水運工程科學研究所，重慶 400016；2.內(nèi)河航道整治技術交通行業(yè)重點實驗室，重慶 400016）

介紹了國內(nèi)外有壓管道系統(tǒng)非恒定流數(shù)值模擬研究的發(fā)展現(xiàn)狀，現(xiàn)階段有壓管道系統(tǒng)非恒定流水力特性數(shù)值模擬中存在主要問題，提出了解決這些問題的主要方法及數(shù)值模擬的發(fā)展趨勢，可為類似研究提供參考。

數(shù)值模擬；有壓管道系統(tǒng)；非恒定流；研究進展

1 發(fā)展現(xiàn)狀

20世紀70年代，計算機的計算速度與效率大大提高，數(shù)值模擬方法模擬水擊成為了現(xiàn)實。數(shù)值模擬主要的優(yōu)勢在于對邊界條件的模擬，在水擊與調壓室涌浪的聯(lián)合計算中有明顯優(yōu)勢，解析法和圖解法慢慢被淘汰。一維數(shù)值模擬中運用最廣最多的是電算法，是以特征線法為基礎發(fā)展起來的。但隨著科學技術的不斷發(fā)展，人們越來越重視對流場內(nèi)部的研究，水擊的研究方向也逐漸由一維向三維模擬轉變，三維數(shù)值模擬是目前水擊的發(fā)展方向，三維數(shù)值模擬與成熟的特征線解法比較起來，更能準確地模擬出流場結構與閥門各種特性在關閉過程中的變化。因此，三維數(shù)值模擬也逐漸成為研究有壓管道系統(tǒng)非恒定流水動力特性的主要模擬方法［1-3］。

2 有壓管道系統(tǒng)的模擬

在水利樞紐等有壓管道系統(tǒng)的設計和運用中，水擊是要考慮的重要因素之一，而為了減小水擊的危害，常需要在有壓管道中設置調壓室。當電站處于某一固定負荷時，機組的引用流量固定，調壓室內(nèi)的水位為某一固定水位。隨著電站負荷的改變，機組引用流量也相應發(fā)生改變，此時調壓室內(nèi)發(fā)生質量波動，波動的穩(wěn)定是調壓室研究中的重要內(nèi)容。若調壓室水位過低，空氣從調壓室進入壓力鋼管形成含氣水錘；若調壓室水位過高，水流則從調壓室溢出，兩者都可能導致重大的安全事故，所以工程上對調壓室的最低、最高涌浪水位非常重視［4-6］。

此外，在初期水電站設計中，一般是把調壓室涌浪和壓力管道水擊分開進行計算，忽略其相互影響。伴隨著計算機技術的發(fā)展與成熟，可以將水擊與調壓井進行聯(lián)合計算。數(shù)值模擬是進行有壓管道非恒定流計算的常用方法。目前，有壓管道系統(tǒng)數(shù)值模擬的方法有很多，主要分為一維數(shù)值模擬和三維數(shù)值模擬。

2.1 一維數(shù)值模擬

2.1.1 控制方程

2.1.1.1 水擊基本方程

水擊基本方程主要有：運動方程、連續(xù)方程。以管道進口為原點，指向下游為正，則有：

運動方程：

運動方程是從管道水體中選取隔離體，應用牛頓第二定律建立的。

連續(xù)方程：

考慮到管道中水體的可壓縮性及管壁的彈性，對非恒定流連續(xù)方程做適當變換而得到連續(xù)方程。式中 g為重力加速度（m/s2）；H為壓力水頭（m）；v為流速（m/s）；x為距離（m）；t為時間（s）；f為摩阻系數(shù)；D為管道內(nèi)徑（m）；θ為管道軸線與水平面夾角（°）；a為水擊波的傳播速度（m/s）。

2.1.1.2 調壓井涌浪方程

運動方程：

連續(xù)方程：

式中 z為庫水位與調壓室水位之差；hf為引水隧洞總的水頭損失；L為引水隧洞總的長度；f為引水隧洞斷面面積；F為調壓室斷面面積；v為引水隧洞中的流速；Q為引水隧洞中的流量；Qt為流入或流出調壓室的流量；Q0為一臺水輪機的引用流量；m為機組的數(shù)量。

2.1.2 水擊和調壓室涌浪聯(lián)合運行的一維數(shù)學模型

聯(lián)合運行的數(shù)學模型中，通常以水擊的運動方程和連續(xù)方程作為基本方程，將調壓室的方程組作為邊界條件引入計算。通過方程的離散、網(wǎng)格劃分、邊界條件分析、時間步長的設定等進行計算。將方程組離散、求全倒數(shù)、整理得特征線方程：

2.2 三維數(shù)值模擬

2.2.1 控制方程

2.2.1.1 連續(xù)方程

連續(xù)方程是任何數(shù)值計算必須遵守的方程，即質量守恒方程：

式中 V為控制體；A為控制面，首項為流入控制體內(nèi)質量增量，第二項為流出控制體質量差。其微分表達式為：

對于不可壓縮流體，ρ為常數(shù)。

2.2.1.2 N-S方程

N-S方程是大多流體流動必須遵守的方程，即為Navier-Stokes方程，也稱為動量守恒方程。

式中 fx，fy，fz分別為單位流體質量力在不同方向上的分量；P為流體內(nèi)應力的張量；sx，sy，sz是流體內(nèi)的廣義源項。對于黏性的可壓縮流體，s，f為零。N-S方程則比較系統(tǒng)地描述了流體的實際流動形態(tài)，所有的流動問題的數(shù)值解法都是圍繞該方程組進行求解的。

2.2.2 方程求解

隨著人們對自然界中紊流運動從整體不斷向微觀發(fā)展，數(shù)值分析在求解明渠水流流動問題的精細化水流模擬要求也越來越高。隨著研究者不斷深入的探討研究與發(fā)展，有限差分法、有限單元法、有限體積法、邊界單元法等越來越成熟，在實際過程中發(fā)揮了重要作用。

2.2.2.1 有限差分法

從實際問題引出的微分方程出發(fā)，將連續(xù)的定解區(qū)域離散為合理分布的網(wǎng)格，把微分方程及定解條件中的微商用相應的差商代替，對得到的差分代數(shù)方程組進行求解，國外學者已從數(shù)學上嚴格證明了差分方程收斂性相容性、穩(wěn)定性及其定解條件，其結果就是原方程在定解區(qū)域上的近似解。有限差分法有較成熟的理論基礎，具體算法簡單易行。

2.2.2.2 有限單元法

將流動問題計算域剖分成各個小區(qū)域，對這些小區(qū)域運用不同的插值函數(shù)選擇區(qū)域基函數(shù)，將區(qū)域結果整合至整個計算域。運用有限差分和變分原理，對不同形狀的單元選擇逼近的插值函數(shù)進行積分處理，適用于邊界復雜的問題。

2.2.2.3 有限體積法

使單元中每個節(jié)點周圍都包含一個控制體，假設網(wǎng)格節(jié)點上的因變量在網(wǎng)格節(jié)點之間的分布規(guī)律，將控制方程對每個控制體體積分，組成一系列離散方程，再結合邊界條件進行求解。牛頓守恒定律自然表現(xiàn)形式，適用于任何類型網(wǎng)格單元。有限體積法可以認為是在有限元法基礎上改進的有限差分法，結合了兩者思想，使計算形式更加統(tǒng)一。目前在很多領域，該方法都得到了廣泛應用。

2.2.2.4 邊界元法

它不需要對整個問題域進行劃分，只需要將問題域的邊界劃分為一系列的單元，邊界微分方程劃分為積分方程，再在邊界上轉化為代數(shù)方程而進行求解。該方法可以降低求解問題的維數(shù)，對于明渠流自由表面的處理較簡單，且精度較高。

2.2.3 紊流數(shù)值模擬理論

2.2.3.1 紊流數(shù)值模擬

對于雷諾數(shù)高于2300的流動稱為紊流，自然界中常見的水流問題多數(shù)為紊流，一般的流體力學往往不能有效地解決紊流問題，而常常結合一定的紊流模型才能達到求解問題的精度。

2.2.3.2 紊流數(shù)值模擬方法

紊流運動是非常復雜的非線性流動，數(shù)值模擬常需要很高的精度才能達到流體運動的實際狀態(tài)，經(jīng)過近些年研究學者的努力，對紊流模型的研究取得了非常顯著成就。研究方法分類如圖1。

圖1 紊流模型分類形式

由于本研究主要應用Reynolds平均法（RANS），現(xiàn)對其進行詳細介紹。

（1）Reynolds平均法（RANS）。紊流運動是非常復雜的非線性流動，對于直接求解N-S方程很難達到其精確的解析解，研究者往往研究紊流的瞬時性，所以常把奈流的瞬時值看作是時間平均值和脈動值的疊加。Reynolds平均法在研究三維紊流問題時具有很高的精度，并得到了廣泛應用。在Reynolds平均法中任一瞬時變量的時間平均值定義為：

在Reynolds平均法中用時均值和脈動值代替瞬時值后，可得到Reynolds時均的紊流運動控制方程。

連續(xù)方程：

N-S方程：

由于時均方程中出現(xiàn)了雷諾應力項uiuj導致方程組不封閉，故引入渦黏性假設。

式中 μt為渦黏性系數(shù)；ui為i時刻的時均速度；k為紊動動能。

根據(jù)對Reynolds應力作出的假設或處理方式不同，目前常用的模型有兩大類：Reynolds應力模型和渦粘模型。Reynolds應力模型方法中，直接構建表示Reynolds應力的方程，然后聯(lián)立求解及建立新的Reynolds應力方程。通常情況下Reynolds應力方程是微分形式，稱為Reynolds應力方程模型。若將Reynolds應力方程的微分形式簡化為代數(shù)應力方程，則稱為代數(shù)應力方程模型［12］。

（2）渦粘模型。主要分為：零方程模型、一方程模型和二方程模型。其中：

零方程模型：最出名的是混摻長度模型推導出μt，混合長度lm和時均速度u關系式，則：

一方程模型：為了彌補混合長度的局限性，建議在湍流時均方程和Reynolds方程的基礎上，建立μt，l和k的關系式，則：

二方程模型：隨著研究者對于求解精度的需要，轉而尋求通用性更高、方法更為精細的紊流模型，更有效地確定特征長度分布。

2.2.4 自由水面的模擬

在水利水電工程、船舶、海洋工程中，水流的自由表面流的模擬具有著重大的現(xiàn)實應用價值。目前，常采用歐拉法和拉格朗日法處理自由液面流動問題。歐拉模型主要體現(xiàn)為在固定的坐標系下或坐標系以一定的方式移動來適應自由表面的移動，而拉格朗日模型則是坐標會隨著流體域的變化而變化，該方法能精確地定位自由面的形狀變化，但是算法極不穩(wěn)定。

隨著自由表面模擬的不斷深入發(fā)展，越來越多的技術被應用到實際工程中來，目前應用最多的有對稱性邊界條件、標高函數(shù)法、MAC法、體積分辨率法等。

2.2.4.1 對稱性邊界

對稱性邊界條件俗稱“剛蓋假定”，是最早應用于自由水面模擬的技術，對稱性邊界條件即假設通過此邊界的所有物理梯度量為零，假設自由液面邊界不會隨時間變動而變動，該方法通常應用于自由表面比較規(guī)則且變化幅度不是太大的流動問題。當處理泄洪、水躍、潰壩水流、波浪等水流問題時，顯然該法的精度不能滿足問題需要。

2.2.4.2 標高函數(shù)法

假設水深是關于x，y，z的單值函數(shù)，以此建立自由表面的運動方程：

方程組達到封閉，即可求解表面任一處的水深h。但該方法也有一定局限性，對于復雜的自由表面流問題，即該法很難有效解決水深A不是關于x，y，t的單值函數(shù)的水流問題。

2.2.4.3 MAC法

MAC法的顯著特點就是能有效解決水深的多值坐標函數(shù)問題。主要表現(xiàn)為在計算域空間的網(wǎng)格內(nèi)設置如同質點的“標記點”，這些點只是用來記錄追蹤流體的運動軌跡，不參與計算，根據(jù)流場的計算，可以得到每個標記點的速度，以此確定被標記的網(wǎng)格單元的狀態(tài)，從而確定出下一時刻自由面的位置。

2.2.4.4 體積分率法（VOF）

體積分率法又稱為體積追蹤法，是將流體在一個網(wǎng)格單元中所占的體積比率當作指標，可以求解空氣與水兩相介質的交互界面，計算體積比率的改變來追蹤分離界面的運動，利用體積分率來標記某計算網(wǎng)格為水、空氣、自由液面，再利用幾何觀念求解此流體體積方程式。釆用線性方法，根據(jù)網(wǎng)格內(nèi)水與空氣的比例來代表自由液面的斜率，表現(xiàn)出網(wǎng)格內(nèi)流體間的界面，進而決定自由液面的形狀。如何決定控制體積內(nèi)空氣與水物理性質所占的比例，則需通過體積分率來決定。表達形式為：

式中 δΩcell為控制體積的大??；δΩwater為控制體內(nèi)水的體積。當F=0時，表示控制體內(nèi)充滿空氣；當F=1時，表示控制體內(nèi)充滿水；當0＜F＜1時，表示控制體內(nèi)同時存在水與空氣，也就是自由液面的位置。

2.2.5 三維數(shù)學模型

水擊和調壓室涌浪聯(lián)合運行三維數(shù)學模型中，運用FLUENT軟件進行計算。邊界條件主要包括：固壁邊界條件、進出口邊界條件。

對于本模型而言，給出的具體邊界條件如圖2。

圖2 有壓管道系統(tǒng)三維數(shù)學模型示意圖

進口邊界類型：進口1為速度進口；出口邊界類型：出口1為壓力出口，采用VOF技術處理水/空氣自由面波動問題；出口2為OUTFLOW。當出口2有閘門時，要模擬閘門的動態(tài)啟閉和有壓管道系統(tǒng)的模擬時，就要引入UDF函數(shù)和動網(wǎng)格技術。對各種閘門的啟閉過程進行動態(tài)仿真時，動網(wǎng)格用于調整和更新閘門運動后的網(wǎng)格，以確保計算網(wǎng)格質量；UDF用于定義閥門的開關規(guī)律，使閥門實現(xiàn)連續(xù)的啟閉過程。為了便于計算，本研究中閘門是線性啟閉。

3 結語

（1）有壓管道系統(tǒng)非恒定流一維數(shù)值模擬已相對成熟，其具有代表性的特征線法有計算速度快，效率高，對于復雜性管道的適應性很強的優(yōu)點，流行了一段時間。但是伴隨特征線法的廣泛使用，其多維模擬能力較差的弊端也越來越明顯，無法得到有壓管道系統(tǒng)流場分布和壓力分布，而三維數(shù)值模擬剛好可以滿足這一點，有壓管道系統(tǒng)非恒定流三維數(shù)值模擬將是今后的發(fā)展方向。

（2）有壓管道系統(tǒng)非恒定流三維數(shù)值模擬中涉及到空氣和水流兩種介質不斷交換問題，可采用成熟的VOF模型進行處理。VOF模型是為解決多相流問題而發(fā)展起來的流體體積模型。其原理是根據(jù)網(wǎng)格單元內(nèi)流體體積的變化來追蹤自由面的變化。近年來，VOF模型在計算兩相流中得到應用。技術已趨于成熟。喻宇［14］運用VOF理論成果模擬出有壓管道水流沖擊氣團的局部三維流場。

（3）國內(nèi)外的研究現(xiàn)狀來看，已有調壓室與水擊聯(lián)合進行的一維數(shù)值模擬和有壓管道的局部三維數(shù)值模擬，但尚無一套完整理論進行調壓室與水擊聯(lián)合進行的三維數(shù)值模擬，這是一個巨大的創(chuàng)新點，建議開展此研究，為類似工程提供理論參考，為有壓管道系統(tǒng)數(shù)值模擬提供研究方向。

［1］段昌國.水電站引水管道水擊的解析與計算［J］.科學通報，1980（3）：134-137.

［2］石寶珍，王學芳，葉宏開.長輸管線水擊計算方法的研究［J］.清華大學學報（自然科學版），1986（1）.

［3］張曉宏，金鐘元.碧口電站壓力引水系統(tǒng)非恒定流的分析與計算［J］.陜西水力發(fā)電，1992（4）.

［4］張師華.水擊壓力和調壓室涌浪聯(lián)合計算的程序編制和分析［J］.武漢水利電力學院學報，1985（1）：23-26.

［5］張健，索麗生.阻抗式調壓室甩負荷涌浪計算顯式公式［J］.水力發(fā)電學報，1999（2）：70-74.

［6］彭天玫，楊效亮.阻抗式調壓井穿井壓力計算及最優(yōu)阻抗孔口的研究［J］.水力發(fā)電學報，1986（3）：61-77.

［7］是勛剛.湍流直接數(shù)值模擬的進展與前景［J］.水動力學研究與進展（A輯），2001（4）：21-27.

［8］羅坤，劉小云，樊建人，等.三維混合層的直接數(shù)值模擬和實驗驗證［J］.中國電機工程學報，2005（6）：13-18.

［9］吳騰虎.管道內(nèi)顆粒懸浮流中大顆粒和湍流相互作用的直接數(shù)值模擬研究［D］.杭州：浙江大學，2011.

［10］李斌，頌平.基于湍流尺度的混合RANS/LES模型［J］.北京航空航天大學學報，2008（1）：31-35.

［11］劉小兵，曾永忠.混流式水輪機中旋渦流的LES法預測［J］.大電機技術，2009（2）：22-27.

［12］馬雷.邊界層理論在低比轉速離心泵葉片設計中的應用［J］.蘭州理工大學學報，2009（2）：41-45.

［13］嚴繼松，廖國玲.有壓管道充水過程水力特性三維數(shù)值模擬［J］.水利水電技術，2015， 11（8）： 219-224.

［14］喻宇.有壓管道水流沖擊氣團的局部三維流場數(shù)值模擬研究［D］.哈爾濱：哈爾濱工業(yè)大學，2011.

Research progress in numerical simulation method of unsteady flow in pressure pipeline system

WANG Yun-li1，2，HU-Hong1，2
（1.Key Laboratory of Inland Waterway Regulation Engineering， Ministry of Transport， PRC，Chongqing 400016，China；2.Southwest Hydraulic Institute for Waterways，Chongqing Jiaotong University，Chongqing 400016，China）

This paper introduces the development of numerical simulation method of unsteady flow in pressure pipeline system at home and abroad，puts forward the existing problems in the simulation of pressure pipeline system at the present stage， and analyzes the possible methods to solve these problems， which can provide a reference for similar research.

numerical simulation；pressure pipe system；unsteady flow；research progress

TV131.4 文獻標識碼：B 文章編號：1672-9900（2017）06-0003-05

2017-09-15

內(nèi)河航道整治技術交通行業(yè)重點實驗室開放基金項目（NHHD-201512）

王云莉（1974-），女（漢族），云南永善人，副研究員，主要從事港口及航道專業(yè)研究，（Tel）13271829317。

王艷肖）