界面、滯留氣團三大部分，并建立 相應的控制方程，根據工程實例確定初始條件以及邊界條件;根據有限體積法劃分計算網 格，并建立離散方程;采用Godunov格式求解水體的離散方程的數值通量，并取得二階精度； 通過基于Runge-Kutta方法對水體的離散方程的數值積分項進行求解，從而得到二階顯式 有限體積法的Godunov格式；給出二階顯式有限體積法的Godunov格式所滿足的CFL條件 (Courant-Friedrichs-Lewy criterion);采用Godunov格式和理想氣體狀態方程聯立實現 水-氣交界面的動態追蹤。
[0050] 下面將結合附圖和實施例對本發明技術方案做進一步的詳細描述。
[0051 ]實施例:如圖2所示的實驗系統來研究水流快速沖擊滯留氣團的瞬變流現象。整個 系統由上游壓力水箱，有機透明管，閥門組成。管道總長8.83m，內徑40mm，管道水平，閥門距 離管道末端3.25m，上游入口壓力為0.16MPa，管道末端封閉。初始時，閥門全關，閥后尾水水 深為20mm。實驗測得水錘波速為400m/s，管道平均摩阻0.075~0.095。水-氣耦合作用過程 由突然開啟下游閥門引起，高速高清攝像機記錄其從全關至全開時間為〇. 〇7s~0.09s。 [0052]本發明的具體步驟為：
[0053]步驟1:將管道系統內瞬變流進行劃分為水體、水-氣交界面、滯留氣團三大部分， 并建立相應的控制方程，根據工程實例確定初始條件以及邊界條件。
[0054]水體的基本微分方程為：
[0055]
[0056]
[0057]其中，沿管線距離X與時間t是自變量;H(x，t)是測壓管水頭;V(x，t)是平均截面速 率;g是重力加速度;a是波速;f為達西-威斯巴哈摩阻系數;D為管徑。
[0058]滯留氣團的控制方程為：
[0059]
[0060] 其中，Ha、Va、匕分別為氣團瞬變壓力、體積和長度;其對應的初始值分別為^〇、Va0、 La0;m為理想氣體狀態方程的多變指數，由于實驗中瞬變過程極為短暫，僅為數秒鐘，氣團的 壓縮膨脹可視為絕熱過程，m為1.4。
[0061 ]水-氣交界面的控制方程為：
[0062]
[0063] Hw=Ha (5)
[0064] 其中，Lf為沖擊水體的長度;1、仏分別水-氣交界面處流速和壓力。
[0065] 本實施例下，計算域為從上游水庫出水口至閥門之間的管道;初始條件為閥門全 關時，初始流速為〇，閥門上游水體初始壓力等于上游水庫靜壓水頭，閥門下游滯留氣團初 始壓力為大氣壓力，尾水壓力考慮重力作用;邊界條件為:管道入口處，水庫提供恒壓邊界， 且等于上游水庫靜壓水頭;下游球閥快速開啟建立水力瞬變，管道末端封閉。
[0066] 步驟2:根據有限體積法劃分計算網格，并建立離散方程。
[0067] 圖3為本實施例下水錘區網格離散系統。對第i個控制體，定義其上、下游界面編號 分別為 i_l/2、i+l/2。
[0068]建立離散方程的具體步驟為：
[0069] (a)微分方程(1)(2)的黎曼問題可近似為含常系數的線性雙曲系統的黎曼問題：
[0070]
(6)
[0071] 其中
F是V的平均值，為一常數。
[0072] (b)在控制單元i(從界面i-1/2到i + 1/2)及時間段At(從t到t+At)內對方程(6) 積分，可以得到流動變量u的離散方程：
[0073]
(7)
[0074] 其中，上標η和n+1分別代表t和t+Δ t時步。
.為u在整個控制體 的平均值;f為單元界面處的通量
為源項；
[0075]步驟3:采用Godunov格式求解水體的離散方程的數值通量，并取得二階精度。
[0076] 進一步地，步驟3包含以下子步驟：
[0077] 步驟3.1:求解水體內部控制單元界面處通量
[0078] 首先，根據Godunov方法，其黎曼問題為以下初值問題：
[0079]
[0080]
[0081] 其中，呢、Uf為在η時步時，u分別到界面i+1/2左、右側兩側的平均值。
[0082] 對矩陣瓦進行簡化計算，可以得到其特征值及特征向量分別為：
[0083]
[0084]
[0085] 因特征向量線性無關，進一步可得到：
[0086]
(12) t=l 1=1
[0087] 求解四個未知系數，β2,可得：
[0088] (Π )
[0089] (14)
[0090] 黎曼問題(方程(8)和(9)) 一般解的原始變量形式為：
[0091] u(x,t)=M(1)+a2K(2) (15)
[0092] 現在，利用方程(15)可得到變量在界面i+1/2處的精確解：
[0093]
(16)
[0094] 從而，對te [tn，tn+1]，任一內部控制單元i(l〈i〈N)，在界面i+1/2處的通量為：
[0095]
[0096] 接著，引入MUSCL-Hancock格式來計算內部單元通量fi+i/2,取得二階精度。具體步 驟為：
[0097] (a)數據重構:采用每個控制單元[Xl-1/2， Xl+1/2]內的分段線性函數代替數據單元 平均值up，，則在極值點處Ilf的值為：
[0098]
(18)
[0099] 其中，4,為選定的適度斜坡向量，用以增加計算格式的精度，并保證解中不出現 虛假振蓀。太另'_由.ΔΗ先檉MTNM⑷限制
[0100]
[0101]
[0102] (b)推算:對每一個單元[Xl-1/2，χ1+1/2]，方程（18)中的邊界外推值Uf，Uf;根據下 式乘以0.5 At推算：
[0103]
[0104] (c)黎曼問題:為計算內部界面通量f1+1/2,需結合下面數據求解傳統的黎曼問題：
[0105]
[0106] 將方程(21)代入方程（17)，即可得到管中充滿水時，通量在時間t=[tn，tn+1]內， 所有內部單元界面i+1/2處的二階格式。
[0107] 步驟3.2:構建虛擬控制單元以求解沖擊水體上下游邊界控制單元界面處通量。為 在邊界面處也取得二階精度，分別在起始控制單元h上游側、終點控制單元N下游側構建兩 個虛擬控制單元1-1、Ιο，以及In+i、In+2,并假定在虛擬單元處的流動信息與邊界處是一致的。 從而可求解邊界黎曼問題，且相應的Godunov通量f 1/2。和fN+i/2也可像內部單元那樣進行計 算。
[0108] 步驟3.3:將負特征線與黎曼向量相結合以求得管道進口邊界控制單元界面處通 量。
[0109] 本實施例中，上游水庫恒水位邊界：
[0110] 在上游邊界處，與負特征線相關的黎曼不變量為：
[0111]
[0112] 其中，HV2 = Hres，Hres為上游水庫靜壓水頭。
[0113] 從而可推彳彳
，:變量111/2(1:) = (!11/2，￥1/2)。根據虛擬單元處 的假定，臨近管道進口的虛擬單元I-dPIo的相應值為：
[0114]
(23)
[0115]步驟4:通過基于Runge-Kutta方法對水體的離散方程的數值積分項進行求解，從 而得到二階顯式有限體積法的Godunov格式。
[0116] 具體實現過程為：
[0117] (a)純對流時：
[0118]
(24)
[0119] (b)利用源項乘以Δ t/2進行更新：
[0120]
(25)
[0121] (c)利用源項乘以At再次更新：
[0122]
(26).
[0123]步驟5 :給出二階顯式有限體積法的Godunov格式所滿足的CFL條件（Courant-Friedrichs-Lewy criterion)。
[0124] CFL條件下的最大時間步長Atmax,CFL:
[0125]
[0126] 其中，Cr為柯朗數。
[0127] 此外，引入的源項滿足以下穩定性約束：
[0128]
[0129]適用于源項的最大時間步長Atmax,s:
[0130]
[0131 ]可推得包含對流項和源項的最大允許時間步長為：
[0132] Δ tmax = m i η ( Δ tmax, cFL , Δ tmax, s ) (30)
[0133] 步驟6:采用Godunov格式和