界面、滞留气团三大部分，并建立 相应的控制方程，根据工程实例确定初始条件以及边界条件;根据有限体积法划分计算网 格，并建立离散方程;采用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格式和