本發(fā)明涉及一種基于有限體積法的輸水管道中空穴流的模擬方法，屬于水電站(泵站)水力學(xué)數(shù)值模擬計(jì)算技術(shù)領(lǐng)域。

背景技術(shù)：

當(dāng)輸水管道系統(tǒng)中壓力低于液體的汽化壓力時，液體可能汽化形成空化泡。這些空化泡可能均勻分布在管道中，形成均質(zhì)空穴流；也可能聚集合并成一個或多個大空穴，發(fā)生水柱分離現(xiàn)象；甚至兩種現(xiàn)象同時存在于系統(tǒng)中?？栈莸漠a(chǎn)生使得系統(tǒng)的瞬變響應(yīng)變得更加復(fù)雜且危害更大，很多學(xué)者均指出分離液柱再彌合后可能產(chǎn)生高于Joukowsky壓力的短期高壓脈沖，導(dǎo)致水力系統(tǒng)破壞。

目前，針對輸水管道中的空化現(xiàn)象，其模擬方法主要為基于特征線法(MOC，Method of Characteristics)提出的離散蒸汽-空穴模型(DVCM,discrete vapor-cavity model)，以及考慮少量自由氣體的離散氣體-空穴模型(DGCM，discrete gas cavity model)。MOC因其計(jì)算簡單，且可以較好地預(yù)測空穴第一次潰滅后的壓力，而在工程計(jì)算中得到了廣泛的應(yīng)用。然而，經(jīng)典的MOC-DVCM模型，在采用細(xì)密網(wǎng)格時，多空穴的潰滅會導(dǎo)致計(jì)算結(jié)果中出現(xiàn)不真實(shí)的壓力脈沖。盡管自由氣體的引入(MOC-DGCM)，在一定程度上抑制了這一脈沖，但計(jì)算精度并未提升很多。此外，由于MOC本身的局限，使得該類方法的解局限于一階精度。而現(xiàn)有的可獲得高階精度解的方法，一方面算法過于復(fù)雜(如Chaiko基于FVM方法提出的均質(zhì)連續(xù)模型)，另一方面計(jì)算耗時太久(如王歡等提出的二維CFD方法)。周領(lǐng)等基于一維FVM方法提出的DVCM模型，盡管相比MOC-DVCM模型明顯抑制了非真實(shí)脈沖并取得了二階精度，但其結(jié)果在最高壓力附近仍存在一定程度的波動。因此，進(jìn)一步地尋求抑制非真實(shí)脈沖的方法仍然非常必要。

技術(shù)實(shí)現(xiàn)要素：

發(fā)明目的：為彌補(bǔ)現(xiàn)有技術(shù)在模擬輸水管道中空穴流時存在的虛假震蕩以及計(jì)算精度低等不足，本發(fā)明基于FVM，通過引入自由氣體，提供了一種算法簡單，易于實(shí)現(xiàn)的模擬方法，以在采用細(xì)密網(wǎng)格時不出現(xiàn)非真實(shí)脈沖并取得二階精度解。

技術(shù)方案：一種基于有限體積法的輸水管道中空穴流的模擬方法(簡稱FVM-DGCM)，具體步驟如下：

步驟1：在FVM體系下，構(gòu)建含自由氣體的瞬變流基本控制方程，根據(jù)工程實(shí)例確定計(jì)算域、初始條件以及邊界條件；

步驟2：根據(jù)FVM劃分計(jì)算網(wǎng)格，并建立離散方程；

步驟3：采用時空均為二階精度的Godunov格式計(jì)算純對流時控制單元界面處的通量；

步驟4：通過基于二階Runge-Kutta離散格式的時間算子分裂法，引入源項(xiàng)，求得離散方程解的二階顯式FVM-Godunov格式；

步驟5：根據(jù)計(jì)算得到的壓力修正壓力并計(jì)算氣穴體積。

進(jìn)一步地，步驟1中，在FVM體系下，構(gòu)建含自由氣體的瞬變流基本控制方程，需在經(jīng)典水錘問題的基礎(chǔ)上假定：(a)引入的自由氣體含量確定，且遵循理想氣體等溫變化規(guī)律；(b)自由氣體集中于控制單元中心，而其余部分仍為純水體；(c)為計(jì)算氣穴體積，每一個控制單元等分為兩個小控制體；(d)每個等分體內(nèi)的壓力及流量為FVM下純水錘問題的解，自由氣體以及空穴對壓力的影響通過引入壓力修正系數(shù)考慮；(e)氣穴體積的變化量可由等分體內(nèi)流量或壓力的變化量計(jì)算。

更進(jìn)一步地，步驟1中，經(jīng)典水錘的基本微分方程為：

其中，沿管線距離x與時間t是自變量；H(x,t)是測壓管水頭；V(x,t)是平均截面速率；g是重力加速度；a是波速；f為達(dá)西-威斯巴哈摩阻系數(shù)；D為管徑。

更進(jìn)一步地，步驟1中，自由氣體遵循理想氣體等溫變化規(guī)律：

其中，T為溫度，M_g為自由氣體質(zhì)量，R_g為氣體常數(shù)，為氣水混合物體積，為標(biāo)準(zhǔn)條件下的氣體空隙率；為氣體壓強(qiáng)絕對值。

更進(jìn)一步地，步驟1中，氣穴體積的連續(xù)性方程為：

其中，為氣穴體積，Q_out和Q_in分別為氣穴體積上下游側(cè)的流量。

進(jìn)一步地，步驟2中，將空間域x離散為N_s個長度為2Δx、時間域t離散為間隔為Δt的控制單元，自由氣體集中于每個控制單元的中心。對第j個控制單元，定義其上、下游界面編號分別為j-1/2、j+1/2。

更進(jìn)一步地，步驟2中，為考慮自由氣體對瞬變流動的影響，并計(jì)算氣穴體積，構(gòu)建其等效網(wǎng)格系統(tǒng)：將每一個控制單元j在空間域上等分為2個小控制體i、i+1，時間域與原始網(wǎng)格相同。每個小控制體長度為Δx，等效網(wǎng)格系統(tǒng)的網(wǎng)格數(shù)為N＝N_s。對第i個小控制體，定義其上、下游界面編號分別為i-1/2、i+1/2。

更進(jìn)一步地，步驟2中，在含自由氣體的FVM網(wǎng)格系統(tǒng)基礎(chǔ)上，建立離散方程的步驟為：

(a)微分方程(1)(2)的黎曼問題可近似為含常系數(shù)的線性雙曲系統(tǒng)的黎曼問題：

其中是V的平均值，為一常數(shù)。

(b)在小控制體i(從界面i-1/2到i+1/2)及時間段Δt(從t到t+Δt)內(nèi)對方程(5)積分，得到流動變量u的離散方程：

其中，上標(biāo)n和n+1分別代表t和t+Δt時步。為u在整個控制體的平均值；f為界面處的通量。

進(jìn)一步地，步驟3中，計(jì)算純對流時界面處的通量的步驟為：

(a)計(jì)算內(nèi)部控制單元界面處通量

根據(jù)Godunov方法，純對流時，對任一控制體i(1<i<N)，界面i+1/2處的通量為：

其中，為在n時步時，u分別到界面i+1/2左、右側(cè)兩側(cè)的平均值。和的計(jì)算方法決定了計(jì)算精度。本發(fā)明中，通過引入MUSCL(Monotone Upstream-centred Scheme for Conservation Laws)Hancock格式來得到空間和時間上的二階精度，同時引入MINMOD限制器以保證解中不出現(xiàn)虛假震蕩。

(b)計(jì)算邊界控制單元界面處通量

為取得二階精度，分別在起始控制體1上游側(cè)、終點(diǎn)控制體N下游側(cè)構(gòu)建兩個虛擬控制體I_-1、I₀，以及I_N+1、I_N+2，并假定在虛擬體處的流動信息與邊界處是一致的。從而可求解邊界黎曼問題，且相應(yīng)的Godunov通量f_1/2和f_N+1/2也可像內(nèi)部單元那樣進(jìn)行計(jì)算。

更進(jìn)一步地，步驟3中，純對流時對流項(xiàng)需滿足CFL(Courant-Friedrichs-Lewy criterion)條件。

進(jìn)一步地，步驟4中，引入源項(xiàng)后，含自由氣體的瞬變流基本微分方程解的二階FVM-Godunov格式為：

其中，為n+1時步，控制單元i在純對流時，流動變量u的通量；為采用時間分裂法第一次更新后的通量。

更進(jìn)一步地，步驟4中，本發(fā)明求解的黎曼問題忽略對流項(xiàng)，故柯朗數(shù)滿足CFL條件。

進(jìn)一步地，步驟5中，根據(jù)計(jì)算得到的壓力判定采用哪種方法修正壓力并計(jì)算氣穴體積：若步驟4中求得的壓力低于液體汽化壓力，則采用方法I修正壓力并計(jì)算氣穴體積；反之，則采用方法II計(jì)算氣穴體積并修正壓力；兩種修正壓力并計(jì)算氣穴體積的方法分別為：

(a)方法I，壓力高于汽化壓力：

對任一控制單元j(1<j<N_s)，其壓力根據(jù)等效網(wǎng)格系統(tǒng)中對應(yīng)的兩個等分控制體i、i+1計(jì)算得到的二階FVM-Godunov格式解計(jì)算：

其中，表示控制單元j的壓力，分別表示控制體i、i+1的壓力；

將方程(9)代入方程(3)，可計(jì)算氣穴體積：

其中，ρ_l表示液體密度；z(j)表示控制單元j處位置高程；H_v表示流體汽化壓力。

為考慮自由氣體對等分控制體的影響，需對等分控制體內(nèi)的壓力進(jìn)行修正。為此，引入壓力修正系數(shù)C_-ap，從而得到修正后的等分控制體內(nèi)的壓力：

其中，0≤C_-ap≤1：C_-ap＝0意味著兩個等分體內(nèi)的壓力均根據(jù)氣穴壓力進(jìn)行修正；C_-ap＝1意味著兩個等分體內(nèi)的壓力不進(jìn)行修正，即忽略氣穴的影響。

(2)方法II，壓力等于或小于汽化壓力：

一旦任一等分體內(nèi)壓力等于或小于液體汽化壓力，則認(rèn)為液體汽化形成蒸汽空穴。此時，氣穴體積根據(jù)氣穴體積連續(xù)性方程(方程(4))計(jì)算。FVM離散格式下，其表達(dá)式為：

其中，分別為控制單元j內(nèi)空穴上、下游側(cè)的流量；表示控制單元j內(nèi)在第n時步氣穴體積。

考慮空穴對等分控制體內(nèi)壓力的影響后，等分控制體內(nèi)壓力修正為：

有益效果：與現(xiàn)有技術(shù)相比，本發(fā)明具有如下優(yōu)點(diǎn)：

(1)本發(fā)明提供的FVM-DGCM模型成功地克服了一維FVM方法在考慮自由氣體并追蹤氣穴、預(yù)測其發(fā)育-潰滅方面的難題，并取得二階精度，同時如MOC類方法一樣簡單且易于實(shí)現(xiàn)；(2)壓力修正系數(shù)C_-ap的引入，成功考慮了自由氣體、空穴對壓力的影響，是采用FVM實(shí)現(xiàn)DGCM模型的關(guān)鍵；(3)壓力修正系數(shù)C_-ap的合理取值使得不管是稀疏網(wǎng)格還是細(xì)密網(wǎng)格，結(jié)果中均不會出現(xiàn)人為的脈沖峰值；(4)在極限情形C_-ap＝1下，該模型忽略了自由氣體的影響，可認(rèn)為是FVM-DVCM，模型適應(yīng)性更強(qiáng)；(5)該模型還可模擬不產(chǎn)生空穴的含氣瞬變流，模型應(yīng)用范圍更廣。

附圖說明

圖1為本發(fā)明的基本流程圖；

圖2為實(shí)施例FVM下的初始網(wǎng)格離散系統(tǒng)圖；

圖3為實(shí)施例FVM下的等效網(wǎng)格離散系統(tǒng)圖；

圖4為實(shí)施例下，壓力修正系數(shù)C_-ap＝1.0時，閥門末端的壓力曲線圖；

圖5為實(shí)施例下，壓力修正系數(shù)C_-ap＝0.9時，閥門末端的壓力曲線圖；

圖6為實(shí)施例下，壓力修正系數(shù)C_-ap＝0.5時，閥門末端的壓力曲線圖；

圖2中：1-自由氣體。

具體實(shí)施方式

下面結(jié)合具體實(shí)施例，進(jìn)一步闡明本發(fā)明，應(yīng)理解這些實(shí)施例僅用于說明本發(fā)明而不用于限制本發(fā)明的范圍，在閱讀了本發(fā)明之后，本領(lǐng)域技術(shù)人員對本發(fā)明的各種等價形式的修改均落于本申請所附權(quán)利要求所限定的范圍。

基于有限體積法的輸水管道中空穴流的模擬方法，按以下步驟進(jìn)行：1.在FVM體系下，構(gòu)建含自由氣體的瞬變流基本控制方程，根據(jù)工程實(shí)例確定計(jì)算域、初始條件以及邊界條件；2.根據(jù)FVM劃分計(jì)算網(wǎng)格，并建立離散方程；3.采用時空均為二階精度的Godunov格式計(jì)算純對流時控制單元界面處的通量；4.通過基于二階Runge-Kutta離散格式的時間算子分裂法，引入源項(xiàng)，求得離散方程解的二階顯式FVM-Godunov格式；5.根據(jù)計(jì)算得到的壓力修正壓力并計(jì)算氣穴體積。圖1給出了本發(fā)明的基本流程圖。

下面將結(jié)合附圖和實(shí)施例對本發(fā)明技術(shù)方案做進(jìn)一步的詳細(xì)描述。

實(shí)施例:為了驗(yàn)證并分析本發(fā)明提供的FVM-DGCM模擬空穴流的模擬效果，選取Simpson于1986年設(shè)計(jì)搭建的水柱分離實(shí)驗(yàn)裝置系統(tǒng)用于驗(yàn)證本發(fā)明方法的有效性。整個系統(tǒng)由上游壓力水箱，銅管，下游球閥，下游壓力水箱組成。銅管總長36m，內(nèi)徑19.05mm，管道斜率1/36。水錘波速為1280m/s。空穴流瞬變由突然關(guān)閉下游球閥引起。實(shí)驗(yàn)工況參數(shù)為：初始流速0.332m/s，上游水庫靜壓水頭23.41m，沿程水頭損失0.33m，水溫24.4℃，閥門關(guān)閉時間0.022s。

本發(fā)明的具體步驟為：

步驟1：在FVM體系下，構(gòu)建含自由氣體的瞬變流基本控制方程，根據(jù)工程實(shí)例確定計(jì)算域、初始條件以及邊界條件。

(a)建立基本控制方程

瞬變流基本微分方程：

其中，沿管線距離x與時間t是自變量；H(x,t)是測壓管水頭；V(x,t)是平均截面速率；g是重力加速度；a是波速；f為達(dá)西-威斯巴哈摩阻系數(shù)；D為管徑。

自由氣體等溫變化狀態(tài)方程：

其中，T為溫度，M_g為自由氣體質(zhì)量，R_g為氣體常數(shù)，為氣水混合物體積，為標(biāo)準(zhǔn)條件下的氣體空隙率，為氣體壓強(qiáng)絕對值。本實(shí)施例下取MOC-DGCM建議值

氣穴體積的連續(xù)性方程：

其中，為氣穴體積，Q_out和Q_in分別為氣穴體積上下游側(cè)的流量。

(b)確定計(jì)算域、初始條件及邊界條件

計(jì)算域：從上游水庫出水口至閥門之間的管道；

初始條件：球閥全開時，初始流速為0.332m/s，各節(jié)點(diǎn)初始壓力根據(jù)上游水庫靜壓水頭減去相應(yīng)的穩(wěn)態(tài)摩阻損失得到；

邊界條件為：管道入口處，水庫提供恒壓邊界，且等于上游水庫靜壓水頭；下游球閥快速關(guān)閉引發(fā)空穴流瞬變。

步驟2：根據(jù)FVM劃分計(jì)算網(wǎng)格，并建立離散方程。

將空間域x離散為N_s個長度為2Δx，時間域t離散為間隔為Δt的控制單元，自由氣體集中于每個控制單元的中心。對第j個控制單元，定義其上、下游界面編號分別為j-1/2、j+1/2。圖2為本實(shí)施例下基于FVM的初始網(wǎng)格離散系統(tǒng)圖。

為考慮自由氣體對瞬變流動的影響，并計(jì)算氣穴體積，構(gòu)建其等效網(wǎng)格系統(tǒng)(見圖3)：將每一個控制單元j在空間域上等分為2個小控制體i、i+1，時間域與原始網(wǎng)格相同。每個小控制體長度為Δx，等效網(wǎng)格系統(tǒng)的網(wǎng)格數(shù)為N＝N_s。對第i個小控制體，定義其上、下游界面編號分別為i-1/2、i+1/2。

在含自由氣體的FVM網(wǎng)格系統(tǒng)基礎(chǔ)上，建立離散方程的具體步驟為：

(a)微分方程(1)(2)的黎曼問題可近似為含常系數(shù)的線性雙曲系統(tǒng)的黎曼問題：

其中是V的平均值，為一常數(shù)。

(b)在小控制體i(從界面i-1/2到i+1/2)及時間段Δt(從t到t+Δt)內(nèi)對方程(5)積分，得到流動變量u的離散方程：

其中，上標(biāo)n和n+1分別代表t和t+Δt時步。為u在整個控制體的平均值；f為單元界面處的通量。

步驟3：采用時空均為二階精度的Godunov格式計(jì)算純對流時控制單元界面處的通量。

進(jìn)一步地，步驟3包含以下子步驟：

步驟3.1：計(jì)算純對流時，內(nèi)部控制單元界面處通量

根據(jù)Godunov方法，純對流時，對任一控制體i(1<i<N)，界面i+1/2處的通量為：

其中，為在n時步時，u分別到界面i+1/2左、右側(cè)兩側(cè)的平均值。和的計(jì)算方法決定了計(jì)算精度。本發(fā)明中，通過引入MUSCL(Monotone Upstream-centred Scheme for Conservation Laws)Hancock格式來得到空間和時間上的二階精度，同時引入MINMOD限制器以保證解中不出現(xiàn)虛假震蕩。

步驟3.2：計(jì)算純對流時，邊界控制單元界面處通量

為取得二階精度，分別在起始控制體1上游側(cè)、終點(diǎn)控制體N下游側(cè)構(gòu)建兩個虛擬控制體I_-1、I₀，以及I_N+1、I_N+2，并假定在虛擬體處的流動信息與邊界處是一致的。從而可求解邊界黎曼問題，且相應(yīng)的Godunov通量f_1/2和f_N+1/2也可像內(nèi)部單元那樣進(jìn)行計(jì)算。

本實(shí)施例中，包含兩種邊界條件：

(a)上游水庫恒水位邊界：

在上游邊界處，與負(fù)特征線相關(guān)的黎曼不變量為：

其中，H_1/2＝H_res，H_res為上游水庫靜壓水頭。

從而可推得變量u_1/2(t)＝(H_1/2,V_1/2)。根據(jù)虛擬單元處的假定，臨近管道進(jìn)口的虛擬單元I_-1和I₀的相應(yīng)值為：

(b)下游閥門邊界：

在下游邊界處，與正特征線相關(guān)的黎曼不變量為：

下游邊界條件為：閥門在閥門關(guān)閉時間T_c內(nèi)關(guān)閉。水頭、流速之間的關(guān)系為：

V_N+1/2＝0,t<T_c (12)

其中，H_D為出口壓力水頭，ξ_v為閥門處損失系數(shù)，變量u_N+1/2(t)＝(H _N+1/2,V _N+1/2)。根據(jù)虛擬單元處的假定，臨近I_N的虛擬單元I_N+1和I_N+2處的相應(yīng)值為：

需注意，步驟3中，對流項(xiàng)應(yīng)滿足CFL條件(Courant-Friedrichs-Lewy criterion)：

C_r≤1 (14)

其中，C_r為柯朗數(shù)。

步驟4：通過基于二階Runge-Kutta離散格式的時間算子分裂法，引入源項(xiàng)，求得離散方程解的二階顯式FVM-Godunov格式。

具體實(shí)現(xiàn)過程為：

(a)純對流時：

(b)利用源項(xiàng)乘以Δt/2進(jìn)行更新：

(c)利用源項(xiàng)乘以Δt再次更新：

此外，本發(fā)明涉及的黎曼問題忽略對流項(xiàng)，故柯朗數(shù)滿足CFL條件。

步驟5：根據(jù)計(jì)算得到的壓力判定采用哪種方法修正壓力并計(jì)算氣穴體積：若步驟4中求得的壓力低于液體汽化壓力，則采用方法I修正壓力并計(jì)算氣穴體積；反之，則采用方法II計(jì)算氣穴體積并修正壓力。

兩種方法分別為：

(a)方法I，壓力高于汽化壓力：

對任一控制單元j(1<j<N_s)，其壓力根據(jù)等效網(wǎng)格系統(tǒng)中對應(yīng)的兩個等分控制體i、i+1計(jì)算得到的二階FVM-Godunov格式解計(jì)算：

其中，表示控制單元j的壓力，分別表示控制體i、i+1的壓力；

將方程(9)代入方程(3)，可計(jì)算氣穴體積：

其中，ρ_l表示液體密度；z(j)表示控制單元j處位置高程；H_v表示流體汽化壓力。

為考慮自由氣體對等分控制體的影響，需對等分控制體內(nèi)的壓力進(jìn)行修正。為此，引入壓力修正系數(shù)C_-ap，從而得到修正后的等分控制體內(nèi)的壓力：

其中，0≤C_-ap≤1：C_-ap＝0意味著兩個等分體內(nèi)的壓力均根據(jù)氣穴壓力進(jìn)行修正；C_-ap＝1意味著兩個等分體內(nèi)的壓力不進(jìn)行修正，即忽略氣穴的影響。

(b)方法II，壓力等于或小于汽化壓力：

一旦任一等分體內(nèi)壓力等于或小于液體汽化壓力，則認(rèn)為液體汽化形成蒸汽空穴。此時，氣穴體積根據(jù)氣穴體積連續(xù)性方程(方程(4))計(jì)算。FVM離散格式下，其表達(dá)式為：

其中，分別為控制單元j內(nèi)空穴上、下游側(cè)的流量；表示控制單元j內(nèi)第n時步的氣穴體積。。

考慮空穴對等分控制體內(nèi)壓力的影響后，等分控制體內(nèi)壓力修正為：

通過以上方法編程計(jì)算后，將FVM-DGCM計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)做對比。圖4-6分別給出了本實(shí)施例下網(wǎng)格數(shù)Ns＝{32,256}，壓力修正系數(shù)C_-ap＝{1.0,0.9,0.5}時，閥門處的壓力曲線圖?？梢钥闯?，壓力修正系數(shù)C_-ap的取值對FVM-DGCM結(jié)果中非真實(shí)壓力脈沖的影響很大。合理選擇C_-ap的取值，可使得不論網(wǎng)格數(shù)為多少，均可有效抑制非真實(shí)壓力脈沖，但過低的取值會導(dǎo)致壓力衰減過大。因此，必須均衡考慮非真實(shí)壓力脈沖、壓力幅值和時間響應(yīng)之間的比重，選取與實(shí)驗(yàn)結(jié)果最接近，且虛假震蕩最少的C_-ap值。綜合考慮，該實(shí)施例下，C_-ap＝0.9且在細(xì)密網(wǎng)格下，結(jié)果最優(yōu)。針對輸水管道中的空穴流現(xiàn)象，相比MOC類方法以及FVM-DVCM，本發(fā)明提出的FVM-DGCM方法，在細(xì)密網(wǎng)格下可以獲得更準(zhǔn)確、穩(wěn)定的結(jié)果。