[发明专利]基于有限体积法的输水管道中空穴流的模拟方法

权利要求书

1.一种基于有限体积法的输水管道中空穴流的模拟方法，其特征在于，采用二阶有限体积法的Godunov格式来模拟输水管道系统中空穴流现象，具体步骤如下：

步骤1：在FVM体系下，构建含自由气体的瞬变流基本控制方程，根据工程实例确定计算域、初始条件以及边界条件；

步骤2：根据FVM划分计算网格，并建立离散方程；

步骤3：采用时空均为二阶精度的Godunov格式计算纯对流时控制单元界面处的通量；

步骤4：通过基于二阶Runge-Kutta离散格式的时间分裂法，引入源项，求得离散方程解的二阶显式FVM-Godunov格式；

步骤5：根据计算得到的压力修正压力并计算气穴体积；

步骤1中，构建的含自由气体的瞬变流基本控制方程包含：

(1)水锤的基本微分方程：

∂H∂t+V∂H∂x+a2g∂V∂x=0---(1)]]>

∂V∂t+V∂V∂x+g∂H∂x+f|V|V2D=0---(2)]]>

其中，沿管线距离x与时间t是自变量；H(x,t)是测压管水头；V(x,t)是平均截面速率；g是重力加速度；a是波速；f为达西-威斯巴哈摩阻系数；D为管径；

(2)自由气体等温变化状态方程：

其中，T为温度，M_g为自由气体质量，R_g为气体常数，为气水混合物体积，为标准条件下的气体空隙率；为气体压强绝对值；

(3)气穴体积的连续性方程：

d∀gdt=Qout-Qin---(4)]]>

其中，为气穴体积，Q_out和Q_in分别为气穴体积上下游侧的流量；

步骤2中，在FVM方法下建立含自由气体的计算网格的方法为：

(a)建立初始网格：

将空间域x离散为N_s个长度为2Δx、时间域t离散为间隔为Δt的控制单元，自由气体集中于每个控制单元的中心；对第j个控制单元，定义其上、下游界面编号分别为j-1/2、j+1/2；

(b)建立等效网格：

将每一个控制单元j在空间域上等分为2个小控制体i和i+1，时间域与原始网格相同；每个小控制体长度为Δx，等效网格系统的网格数为N＝N_s；对第i个小控制体，定义其上、下游界面编号分别为i-1/2、i+1/2；

步骤2中，在含自由气体的FVM网格系统基础上，建立离散方程的步骤为：

(a)微分方程(1)(2)的黎曼问题可近似为含常系数的线性双曲系统的黎曼问题：

∂u∂t+∂f(u)∂x=s(u),f(u)=A‾u---(5)]]>

其中是V的平均值，为一常数；

(b)在小控制体i(从界面i-1/2到i+1/2)及时间段Δt(从t到t+Δt)内对方程(5)积分，得到流动变量u的离散方程：

Uin+1=Uin+ΔtΔx(fi+1/2-fi-1/2)+ΔtΔx∫i-1/2i+1/2sdx---(6)]]>

其中，上标n和n+1分别代表t和t+Δt时步；为u在整个控制体的平均值；f为界面处的通量；

步骤3中，计算纯对流时界面处的通量的步骤为：

(a)计算内部控制单元界面处通量

根据Godunov方法，纯对流时，对任一控制体i(1<i<N)，界面i+1/2处的通量为：

fi+1/2=A‾i+1/2ui+1/2(t)=12A‾i+1/2{1a/gg/a1ULn--1a/gg/a-1URn}---(7)]]>

其中，为在n时步时，u分别到界面i+1/2左、右侧两侧的平均值；和的计算方法决定了计算精度；通过引入MUSCL Hancock格式来得到空间和时间上的二阶精度，同时引入MINMOD限制器以保证解中不出现虚假震荡；

(b)计算边界控制单元界面处通量

为取得二阶精度，分别在起始控制体1上游侧、终点控制体N下游侧构建两个虚拟控制体I_-1、I₀，以及I_N+1、I_N+2，并假定在虚拟体处的流动信息与边界处是一致的；从而可求解边界黎曼问题，且相应的Godunov通量f_1/2和f_N+1/2也可像内部单元那样进行计算；

步骤5中，根据计算得到的压力判定采用哪种方法修正压力并计算气穴体积：若步骤4中求得的压力低于液体汽化压力，则采用方法I修正压力并计算气穴体积；反之，则采用方法II计算气穴体积并修正压力；两种修正压力并计算气穴体积的方法分别为：

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

对任一控制单元j(1<j<N_s)，其压力根据等效网格系统中对应的两个等分控制体i、i+1计算得到的二阶FVM-Godunov格式解计算：

(Hg)jn+1=(Hin+1+Hi+1n+1)/2---(9)]]>

其中，表示控制单元j的压力，分别表示控制体i、i+1的压力；

将方程(9)代入方程(3)，可计算气穴体积：

(∀g)jn+1=p0*α0∀ρlg[(Hg)jn+1-z(j)-Hv]---(10)]]>

其中，ρ_l表示液体密度；z(j)表示控制单元j处位置高程；H_v表示流体汽化压力；

为考虑自由气体对等分控制体的影响，需对等分控制体内的压力进行修正；为此，引入压力修正系数C_-ap，从而得到修正后的等分控制体内的压力：

Hi(n+1)=C-apHin+1+(1-C-ap)(Hg)jn+1---(11)]]>

Hi+1(n+1)=C-apHi+1n+1+(1-C-ap)(Hg)jn+1---(12)]]>

其中，0≤C_-ap≤1：C_-ap＝0意味着两个等分体内的压力均根据气穴压力进行修正；C_-ap＝1意味着两个等分体内的压力不进行修正，即忽略气穴的影响；

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

一旦任一等分体内压力等于或小于液体汽化压力，则认为液体汽化形成蒸汽空穴；此时，气穴体积根据气穴体积连续性方程(方程(4))计算；FVM离散格式下，其表达式为：

(∀g)jn+1=(∀g)jn+(Qjn+1-Qujn+1)Δt=(∀g)jn+(Qin+1-Qi+1n+1)Δt---(13)]]>

其中，分别为控制单元j内空穴上、下游侧的流量；表示控制单元j内在第n时步的气穴体积；

考虑空穴对等分控制体内压力的影响后，等分控制体内压力修正为：

Hi(n+1)=Hi+1(n+1)=(Hg)jn+1---(14).]]>