[发明专利]耦合能量方程的输水管道气-液两相均质流的模拟方法

权利要求书

1.一种基于有限体积法耦合能量方程的输水管道气-液两相均质流的模拟方法，其特征在于，包括如下步骤：

S1：通过引入能量方程构建含自由气体均质两相流结合能量方程的双曲偏微分方程系统模型，根据模拟工况确定计算域、初始条件以及边界条件；

S2：构建出双曲偏微分方程系统模型下各基础波的存在条件以及判定模型；

S3：构建双曲偏微分方程系统模型下断点处三种近似黎曼求解器，分别为一种迭代和两种非迭代近似黎曼求解器；

S4：通过重建流动变量，使得系统模型在空间上达到二阶精度；

S5：根据已有的边界条件，通过牛顿迭代法利用黎曼不变量迭代得到二阶的边界值；

S6：引入源项，通过基于二阶Longge kuta离散格式的时间分裂法，求得时间梯度上二阶有限体积法Godunov格式离散方程，得到二阶精度；

S7：建立稳定性约束条件，更新初始变量进行下一时步计算，获取到最大允许时间步长；

所述步骤S1中双曲偏微分方程系统模型的构建是通过在双曲偏微分方程系统体系下，构建含自由气体的两相均质流基本控制方程；

所述步骤S1中含自由气体的两相均质流基本控制方程具体如下：

包含能量方程的水锤基本方程：

其中，矩阵代表气液混合流体的特征变量；矩阵代表通量；沿管线距离x与时间t是自变量；ρ_m(x,t)是平均截面气液混合流体的密度；V(x,t)是气液混合体平均截面速度；P(x,t)是气液混合流体的绝对压力；E(x,t)是气液混合流体单位体积总能量；

气-液两相流水锤波速与压力的关系方程：

其中，c_m为气液混合流体的波速；C为纯液体下水锤波速；P_ref为参考压力；为参考压力下气液两相混合流体的密度；ψ_ref为参考压力下混合流体内气体的初始体积分数；

气-液两相流压力与密度的关系方程：

其中，ρ_m为气液混合流体的密度，P为气液混合流体的压力；C为纯液体下水锤波速；P_ref为参考压力；为参考压力下气液两相混合流体的密度；ψ_ref为参考压力下混合流体内气体的初始体积分数；

所述步骤S2中基础波包括接触波、激波、稀疏波，构建方法具体如下：

接触波：

对于两相流系统中的接触波关系，利用欧拉方程中的特征结构分析断点处的初值，利用一般黎曼不变量，穿过波λ＝V，得到，

这里Pr代表参考压力；u₁、u₂、u₃见公式(2)；ρ_mr为参考压力下气液混合流体的密度；c_m为气液混合流体的波速；e_r为参考压力下气液混合体的内能；进一步简化后，黎曼不变量写做：

这里星号*代表中间状态，P、V、ρ_m分别代表混合流体的压力、速度、密度，L和R分别代表黎曼问题中接触波左右区间，在等熵条件下，穿过接触波，流体特征不发生变化；

稀疏波：

穿过波λ＝V-c_m，利用一般黎曼不变量，得到

进一步简化后，该微分方程写做

dV+c_md(lnp_m)＝0   (9)

对公式(9)积分后，得到

其中，混合流体平均波速同样地，穿过右侧的稀疏波λ＝V+c_m，得到

激波：

对于左侧激波，假设常波速为S_L，以此为参考系，流体速度V_*、V_L转化为相对速度和

结合兰金-雨果关系，进一步写做

将公式(13)带入(14)得到

引入热力学参数比焓h＝e+Pν和比容ν，得到方程

已知单位质量混合流体内能的表达式

这里ρ_mr,e_r分别代表参考压力P_r下气液混合体的密度和内能，结合公式(17)、(18)，得到

通过公式(16)消去压力项P_*、P_L，得到左侧激波波速

同样地，穿过右激波得到波速

这里L和R分别对应波的左右侧；

所述步骤S2中判定模型的构建方法为：

波型判定：

穿过左右侧非线性波f_L(ρ_m*，U_L)＝V_L-V_*和f_R(ρ_m*，U_R)＝V_*-V_R，得到有关混合流体中间状态处密度ρ_m*代数方程

f(ρ_m*，U_L，U_R)≡f_L(ρ_m*，U_L)+f_R(ρ_m*，U_R)+ΔV＝0       (22)

接下来给出不同波型下f_L(ρ_m*，U_L)和f_R(ρ_m*，U_R)的表达式；

对于左侧波，当出现激波时，基于兰金-雨果关系，公式(13)写做

同时公式(14)转化为

这里Q_mL为混合物质量通量，将公式(23)带入(24)得到

联合公式(24)(25)，消除质量通量Q_mL，得到f_L表达式

进一步地，得到

将公式(27)转化为更一般形式

这里定义参数对于左激波，跨过激波，满足f_L＞0,P_L＜P_*和ρ_mL＜ρ_m*；

当左侧出现稀疏波时，由公式(10)，得到公式f_L满足

同样地，将其简化为

这里定义参数穿过稀疏波，满足f_L0,P_LP_*andρ_mlρ_m*；

同理，得到右侧波f_R满足公式，汇总左右侧代数方程，得到有关黎曼问题中中间状态处气液混合流体速度和密度的表达式

或者

其中K＝L或者R，

接下来，定义

结合方程(22)，得到以下波型判定公式

对于双稀疏波和双激波类型，满足f_min·f_max0，既

结合方程(32)，得到

这里这里引入变量K＝1/|ρ_R-ρ_L|，公式(36)转化为

这里当ρ_mRρ_mL时，当ρ_mRρ_mL时，展开公式(37)，得到

最终，给出以下判定公式

所述步骤S3中三种黎曼求解器的构建方法为：

迭代黎曼求解器：

基于密度方程(22)，利用牛顿迭代法得到下式

这里给出初值ρ_m*⁰的计算公式

两种非迭代黎曼求解器：

对公式(31)进行泰勒展开得到第一个一阶近似的非迭代黎曼求解器

第一个非迭代求解器为两稀疏波求解器，在激波出现时会出现计算误差，构建第二个非迭代求解器

这里K＝L或R，利用第一个非迭代求解器得到T_K的初值T_K^*，并将其带入公式(62)，得到

所述步骤S4中通过含有Total Variation Diminishing形式的MUSCL-Hancoke格式，得到重建的流动变量，具体如下：

S4-1：引入斜率限制器

其中，为斜率限制器参数；

S4-2：流体变量重构

S4-3：进一步变量重构

所述步骤S5具体为：

当为上游边界时，当已知气液混合体压力时，未知量气液混合体密度和波速和由物理方程(3)和(4)得到，结合公式(31)，混合流体的速度写做

这里上标n+1代表计算时步；

当已知混合流体的速度时，假设压力通过公式(4)利用牛顿迭代法得到，通过公式(3)进一步得到波速结合通过以下迭代步骤得到一个新的密度

以及

这里j代表迭代次数；

所述步骤S6中引入源项格式离散方程的求解过程为：

S6-1：对流项

S6-2：Δt/2时步更新

S6-3：Δt时步更新

源项S中梯度(S₀)_i定义为

其中，均为时间分裂法中不同时步计算的中间变量；

所述步骤S7具体为：

S7-1：由CFL得到：

S7-2：显式二阶龙格-库塔离散化约束：

S7-3：给出包括对流部分和源项的最大允许时间步长：