CN113378493A - 一种全质量守恒的多相流数值模拟方法 - Google Patents

Abstract

本发明公开了一种全质量守恒的多相流流动数值模拟方法，包括：步骤1，根据所模拟的页岩气藏中气水两相流动状态，推导出多相流流动控制方程；步骤2，根据多相流流动控制方程得出离散格式的饱和度方程；根据离散格式的饱和度方程得到离散格式的压力方程；步骤3，循环求解离散格式的压力方程和饱和度方程，得出多相流流动数值。通过循环求解压力方程和饱和度方程，并且在计算过程中，与饱和度相关参数的计算不再取决于当前时间步的初始条件，而是依赖于当前或前一迭代步的计算结果，这样与IMPES方法对比，该方法不仅满足全质量守恒限制条件，而且提高了计算稳定性和计算效率。

Description

一种全质量守恒的多相流数值模拟方法

技术领域

本发明涉及地层多相流流动数值模拟领域，尤其涉及一种用于页岩气藏中的气水两相流相流流动计算的全质量守恒的多相流数值模拟方法。

背景技术

目前，有研究表明，页岩气藏中往往共存着大量的初始水，如美国的五大页岩储层：Woodford、Marcellus、Fayetteville、Haynesville和Barnett盆地分别含有0.1,0.12～0.35,0.15～0.35,0.15～0.2and 0.25～0.3的初始水，中国位于四川的页岩储层含有约0.1～0.45的初始水。另外，页岩气的开采，往往需要结合大规模水力压裂施工。水力压裂过程中也需要向地层注入大量的水分。这些水的存在，对页岩气的流动具有较大的影响。

目前，多相流流动数值模拟方法是处理储层中的多相流流动，特别是页岩气气水两相流流动问题的主要方法，主要有全隐式法和隐压显饱法，其中，全隐式法同时处理压力项和饱和度项，具有全质量守恒的优点，但其计算成本和内存需求较高，不适合网格数量较大的情况；隐压显饱方法采用隐式求解压力项及显示求解饱和度项，具有计算效率高、占用内存小的优点，但其计算结果质量不守恒且不稳定。

因此，如何提供一种全质量守恒、计算稳定且效率高的数值模拟计算方法，是处理多相流流动需要解决的问题。

发明内容

基于现有技术所存在的问题，本发明的目的是提供一种全质量守恒的多相流数值模拟方法，能解决现有全隐式法进行多相流流动数值模拟计算，所存在的计算成本和内存需求较高的问题。

本发明的目的是通过以下技术方案实现的：

本发明实施例提供一种全质量守恒的多相流流动数值模拟方法，包括以下步骤：

步骤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^{n +1}指下一时间步；α指时间步增长系数；Δtⁿ指当前时间步；Δt_max指最大时间步长；

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

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

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

由上述本发明提供的技术方案可以看出，本发明实施例提供的全质量守恒的多相流数值模拟方法，其有益效果为：

通过循环求解压力方程和饱和度方程，并且在计算过程中，与饱和度相关参数的计算不再取决于当前时间步的初始条件，而是依赖于当前或前一迭代步的计算结果，这样与IMPES方法对比，该方法不仅满足全质量守恒限制条件，而且提高了计算稳定性和计算效率。该方法能用于处理储层中的多相流流动问题，特别是页岩气的气、水两相流流动问题。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的全质量守恒的多相流流动数值模拟方法流程图；

图2为本发明实施例提供的全质量守恒的多相流流动数值模拟方法的求解式(5)离散格式的压力方程和所述式(3)与式(4)的离散格式的饱和度方程的流程图；

图3为本发明实施例的方法与现有IMPES方法所计算的总气量(累积气产量与地层残存气体的和)随生产时间的关系图；

图4为本发明实施例的方法与现有IMPES方法所计算的累积气产量随生产时间的关系图；

图5为本发明实施例的方法与现有IMPES方法所计算的总水量(累积水产量与地层残存水的和)随生产时间的关系图；

图6为本发明实施例的方法与现有IMPES方法所计算的累积水产量随生产时间的关系图；

图7为本发明实施例中的地层网格划分示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述；显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例，这并不构成对本发明的限制。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

首先对本文中可能使用的术语进行如下说明：

术语“和/或”是表示两者任一或两者同时均可实现，例如，X和/或Y表示既包括“X”或“Y”的情况也包括“X和Y”的三种情况。

术语“包括”、“包含”、“含有”、“具有”或其它类似语义的描述，应被解释为非排它性的包括。例如：包括某技术特征要素(如原料、组分、成分、载体、剂型、材料、尺寸、零件、部件、机构、装置、步骤、工序、方法、反应条件、加工条件、参数、算法、信号、数据、产品或制品等)，应被解释为不仅包括明确列出的某技术特征要素，还可以包括未明确列出的本领域公知的其它技术特征要素。

术语“由……组成”表示排除任何未明确列出的技术特征要素。若将该术语用于权利要求中，则该术语将使权利要求成为封闭式，使其不包含除明确列出的技术特征要素以外的技术特征要素，但与其相关的常规杂质除外。如果该术语只是出现在权利要求的某子句中，那么其仅限定在该子句中明确列出的要素，其他子句中所记载的要素并不被排除在整体权利要求之外。

除另有明确的规定或限定外，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如：可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本文中的具体含义。

当浓度、温度、压力、尺寸或者其它参数以数值范围形式表示时，该数值范围应被理解为具体公开了该数值范围内任何上限值、下限值、优选值的配对所形成的所有范围，而不论该范围是否被明确记载；例如，如果记载了数值范围“2～8”时，那么该数值范围应被解释为包括“2～7”、“2～6”、“5～7”、“3～4和6～7”、“3～5和7”、“2和5～7”等范围。除另有说明外，本文中记载的数值范围既包括其端值也包括在该数值范围内的所有整数和分数。

术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述和简化描述，而不是明示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本文的限制。

下面对本发明所提供的全质量守恒的多相流流动数值模拟方法进行详细描述。本发明实施例中未作详细描述的内容属于本领域专业技术人员公知的现有技术。本发明实施例中未注明具体条件者，按照本领域常规条件或制造商建议的条件进行。本发明实施例中所用试剂或仪器未注明生产厂商者，均为可以通过市售购买获得的常规产品。

如图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时间步的水相饱和度；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^{n +1}指下一时间步；α指时间步增长系数，可根据需要设定；Δtⁿ指当前时间步；Δt_max指最大时间步长；

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

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

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

上述方法步骤32中，采用线性求解器求解所述式(5)的离散格式的压力方程。优选的，线性求解器采用GMRES求解器。

上述步骤34中，预设的同时包含判断δp和δS_w是否收敛的收敛条件为：δp＜pⁿ×10^-6和δS_w＜1×10^-6。

本发明的多相流流动数值模拟方法，用以处理储层中的多相流流动问题，特别是页岩气的气水两相流流动问题，通过循环求解压力方程和饱和度方程，并且在计算过程中，与饱和度相关参数的计算不再取决于当前时间步的初始条件，而是依赖与当前或前一迭代步的计算结果，这样不仅满足全质量守恒条件，同时提高了计算稳定性和计算效率。

下面对本发明实施例具体作进一步地详细描述。

如图1所示，本实施例提供一种全质量守恒的多相流流动数值模拟方法，能在满足全质量守恒条件下，处理页岩气的气、水两相流流动问题，并且计算稳定性和计算效率高，包括如下步骤：

步骤1，本文以页岩气藏中气、水两相流动为例，推导多相流流动控制方程：

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

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

步骤2，对式(1)和(2)进行有限体积离散，将它们分别转化为：

式(3)和(4)中，

为n+1时间步的气相压力；

为n+1时间步的气相饱和度；

为n+1时间步的水相压力；

为n+1时间步的水相饱和度；A、B、C均为与n或(n+1)时间步的压力和饱和度相关的系数；

进一步地，将式(3)和(4)进行整合并消去饱和度项，得到式(5)：

式(5)为离散格式的压力方程，式(3)、(4)为离散格式的饱和度方程；

步骤3，按以下步骤循环求解式(5)的离散格式的压力方程，式(3)、(4)的离散格式的饱和度方程，直至得出所有时间步结合后的压力值和饱和度值，流程参见图2，包括：

步骤31，在对每一个时间步Δtⁿ迭代计算的第一个迭代步中，将pⁿ值赋予p^v，

值赋予

其中，n为时间步的步数，v为每一个时间步的迭代计算中的迭代步的步数；pⁿ为计算出的时间步n的压力值；p^v为迭代步v计算出的压力值；

为计算出的时间步n的水相饱和度值；

为迭代步v计算出的水相饱和度值；

步骤32，依据式(5)求解δp，且依据p^v＝p^v+δp对p^v进行更新；采用线性求解器对式(5)的求解，优选的，本步骤中采用GMRES求解器对式(5)进行求解；

步骤33，依据式3或4求解δS_w，且依据

对

进行更新；

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

步骤35、36、37，通过公式Δtⁿ⁺¹＝min(αΔtⁿ,Δt_max)进行时间步控制，公式中，Δt^{n +1}指下一时间步；α指时间步增长系数；Δtⁿ指当前时间步；Δt_max指最大时间步长；其中，步骤35中，通过公式Δtⁿ⁺¹＝α·Δtⁿ对下一时间步进行计算；步骤36中，判断Δtⁿ⁺¹是否大于Δt_max，若是则进行步骤37，若否则将n加1，将

值赋予

p^v值赋予pⁿ⁺¹；步骤37中，将Δt_max赋值予Δtⁿ⁺¹；

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

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

步骤310，按步骤31至步骤39重复进行后续时间步的求解，直至求解出所有时间步结束后的压力值和饱和度值，即完成多相流流动数值模拟计算。

本实施例的全质量守恒的多相流流动数值模拟方法，能对页岩气藏中经多段水力压裂水平井的气、水两相流动进行模拟，本实施例方法的基本数据如表1所示。

表1为本实施例参数表

本实施例方法与IMPES方法所计算的总气量(累积气产量与与地层残存气体的和)随生产时间的关系如图3所示；

本实施例方法与IMPES方法所计算的累积气产量随生产时间的关系如图4所示。

本实施例方法与IMPES方法所计算的总水量(累积水产量与地层残存水的和)随生产时间的关系如图5所示。

本实施例方法与IMPES方法所计算的累积水产量随生产时间的关系如图6所示。

本实施例中地层网格划分如图7所示。

本实施例方法与IMPES方法相比，计算结果、计算效率上的差异如表2所示，

表2为本实施例方法与IMPES方法有关计算效率的数据对比表

CPU＝CPU时间；NT＝总时间步；NI＝总迭代步。

通过上述可以看出，本发明模拟方法的计算稳定性和计算效率均高于现有的IMPES方法。

本领域普通技术人员可以理解：实现上述实施例方法中的全部或部分流程是可以通过程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (5)

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ⁿ指当前时间步；Δt_max指最大时间步长；

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

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

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

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

Images