[发明专利]一种全质量守恒的多相流数值模拟方法

权利要求书

1.一种全质量守恒的多相流流动数值模拟方法，其特征在于，包括以下步骤：

步骤1，根据所模拟的页岩气藏中气相和水相流动状态，推导出多相流流动控制方程为：

所述式(1)和(2)中，下标g的参数为气相参数，下标w的参数为水相参数；φ为孔隙度；S为饱和度；B为地层体积系数；V_E为单位质量岩石吸附的气体量；ρ_R为岩石密度；u为流体流速；q为气体产速；

所述式(1)和(2)中，单位质量岩石吸附的气体量V_E用Langmuir吸附公式进行计算，流体流速采用Darcy定律进行计算；

步骤2，对所述式(1)和式(2)均进行有限体积离散后分别转化得到式(3)与式(4)的离散格式的饱和度方程，为：

所述式(3)和(4)中，为n+1时间步的气相压力；为n+1时间步的气相饱和度；为n+1时间步的水相压力；为n+1时间步的水相饱和度；A、B、C均为与n或n+1时间步的压力和饱和度相关的系数；

对所述式(3)和(4)进行整合并消去饱和度项得到式(5)的离散格式的压力方程，为：

步骤3，采用以下步骤循环求解所述式(5)离散格式的压力方程和所述式(3)与式(4)的离散格式的饱和度方程，直至得出所有时间步的压力值和饱和度值，包括：

步骤31，在对每一个时间步Δtⁿ迭代计算的第一个迭代步中，将pⁿ值赋予p^v，值赋予其中，n为时间步的步数，v为每一个时间步的迭代计算中的迭代步的步数；pⁿ为计算出的时间步n的气相和水相的平均压力值；p^v为迭代步v计算出的气相和水相的平均压力值；为计算出的时间步n的水相饱和度值；为迭代步v计算出的水相饱和度值；

步骤32，根据所述式(5)求解δp，并根据p^v＝p^v+δp对p^v进行更新；

步骤33，根据所述式(3)或(4)求解δS_w，并根据对进行更新；

步骤34，根据预设的同时包含判断δp和δS_w是否收敛的收敛条件判断当前迭代步是否收敛性，当满足该收敛条件时，进入步骤35进行下一时间步的计算，否则按步骤31至步骤34对当前时间步重新进行迭代计算；

步骤35～37，根据公式Δtⁿ⁺¹＝min(αΔtⁿ,Δt_max)对时间步进行控制，公式中，Δtⁿ⁺¹指下一时间步；α指时间步增长系数；Δt_max指最大时间步长；

步骤38，当迭代步数v达到预定最大迭代步数M仍无法收敛时，计算流程进入步骤39；

步骤39，将当前时间步的时间步长进行减半，并以减半后的时间步按步骤31至步骤38的流程重新进行计算，直至求解得出当前时间步结束后的压力值和饱和度值，进入步骤310；

步骤310，按步骤31至步骤39重复进行后续时间步的求解，直至求解出所有时间步结束后的压力值和饱和度值。

2.根据权利要求1所述的全质量守恒的多相流流动数值模拟方法，其特征在于，所述步骤34中，预设的同时包含判断δp和δS_w收敛的收敛条件为：δppⁿ×10^-6和δS_w1×10^-6。

3.根据权利要求1或2所述的全质量守恒的多相流流动数值模拟方法，其特征在于，所述步骤32中，采用线性求解器求解所述式(5)的离散格式的压力方程。

4.根据权利要求3所述的全质量守恒的多相流流动数值模拟方法，其特征在于，所述线性求解器采用GMRES求解器。

5.根据权利要求1或2所述的全质量守恒的多相流流动数值模拟方法，其特征在于，所述步骤35～37中，

步骤35中，通过公式Δtⁿ⁺¹＝α·Δtⁿ对下一时间步进行计算；

步骤36中，判断Δtⁿ⁺¹是否大于Δt_max，若是则进行步骤37，若否则将n加1，将值赋予p^v值赋予pⁿ⁺¹；

步骤37中，将Δt_max赋值予Δtⁿ⁺¹。