一种页岩油藏超临界CO2压裂精细模拟方法

一种页岩油藏超临界co2压裂精细模拟方法
技术领域
1.本发明涉及油藏数值模拟技术领域，更具体地，涉及一种页岩油藏超临界co2压裂精细模拟方法。


背景技术：

2.页岩油藏储层致密，必须采用增产改造技术才能实现商业开发，其中水力压裂技术较成熟，但存在水资源消耗大和环境污染等缺点。超临界co2压裂技术引起了广泛的关注，目前研究主要集中在室内和矿场试验分析，如何精细模拟超临界co2压裂缝网扩展过程尚未解决。超临界co2压裂与水力压裂存在本质区别，超临界co2具有独特的物理化学性质，表面张力接近于零，很容易进入页岩孔隙，使岩石更容易破裂，表明co2浸泡后岩石断裂韧性减小；超临界co2流体物性参数对压力温度变化敏感，滤失进入储层与原油产生多相流动传热现象，从而导致相态变化复杂的难题；常规水力压裂模型尚未考虑超临界co2与原油混相流动传热的影响，不再适用于模拟超临界co2压裂裂缝扩展过程。
3.因此，亟需发展一种页岩油藏超临界co2压裂精细模拟流程与方法，综合考虑co2浸泡对岩石断裂韧度的影响、超临界co2与原油流动传热对其相态变化的影响以及页岩储层油气特殊流动机理，此方法填补了超临界co2压裂数值模拟技术的空缺，具有计算精度高和可靠性高的优势。


技术实现要素：

4.提供了本发明以解决现有技术中存在的上述问题。本发明提供一种页岩油藏超临界co2压裂精细模拟方法，填补当前技术存在的空缺，获得一种计算精度高和可靠性高的模拟流程和方法。
5.本发明具体采用如下技术方案：
6.一种页岩油藏超临界co2压裂精细模拟方法，所述方法包括：
7.获取岩石断裂韧度参数；
8.基于天然裂缝分布情况，建立页岩油藏几何模型，并对含天然裂缝的页岩油藏几何模型进行网格划分，以保证天然裂缝位于网格边界单元；
9.插入全局内聚区单元，其中天然裂缝所在网格边界单元标记为天然裂缝内聚区单元，其余网格边界单元标记为基质内聚区单元，并根据岩石断裂韧度参数赋予基质内聚区单元属性；
10.基于离散裂缝模型，建立超临界co2多相流动组分模型，采用有限体积法对控制方程进行数值离散化，应用牛顿迭代法对离散得到的流体压力、相饱和度和组分摩尔浓度非线性方程组进行求解；
11.基于离散裂缝模型，结合热对流和热传导两种方式，建立超临界co2‑ꢀ
岩石热量传输模型，采用有限体积法对控制方程进行数值离散化，采用牛顿迭代法对离散得到的温度场方程进行求解；
12.建立岩体变形和裂缝扩展模型，结合流体压力和热应力的影响，采用有限元方法对控制方程进行数值离散化，采用牛顿迭代法对离散得到的应力场方程进行求解；
13.基于固定应力分裂方法，构建温度场-渗流场-应力场的迭代耦合模型，采用mpi并行计算求解各物理场，循环迭代至迭代耦合模型收敛，根据内聚区损伤演化准则判断裂缝是否扩展，以更新裂缝网格信息，得到裂缝分布形态，实现页岩油藏超临界co2压裂精细数值模拟。
14.进一步地，通过开展超临界co2浸泡后页岩断裂韧度测试实验获取岩石断裂韧度参数。
15.进一步地，所述开展超临界co2浸泡后页岩断裂韧度测试实验的具体步骤包括：
16.将页岩制作成多个半圆形试件，沿圆心将岩石切割得到短槽；
17.将页岩试件浸泡在超临界co2；
18.试件浸泡结束后，开展三点弯曲试验，获得不同浸泡时间下页岩断裂韧度参数。
19.进一步地，所述不同浸泡时间包括2h，4h，8h，12h和24h。
20.进一步地，通过如下方法获取天然裂缝分布情况：
21.根据天然裂缝统计信息，采用随机方法生成天然裂缝分布；
22.所述对含天然裂缝的页岩油藏几何模型进行网格划分，包括：
23.二维模型采用三角形单元剖分，三维模型采用四面体单元剖分，且以天然裂缝为内边界约束，获得匹配型网格。
24.进一步地，所述插入全局内聚区单元的方法为，对每个单元的所有边界进行循环插入，每条边界拆分成两边组成一个内聚区单元，初始网格节点自由度包括压力和温度，更新后网格节点自由度为位移值；所述根据岩石断裂韧度参数测试结果赋予内聚区单元属性方法为，利用岩石断裂韧度与断裂能之间的关系，求解得到岩石断裂能，将其赋给基质内聚区单元断裂能属性；
25.进一步地，所述超临界co2多相流动组分模型包括烃类组分和co2组分质量守恒方程、水组分质量守恒方程、油水三相运动方程、气相运动方程及辅助方程，其中，
26.所述烃类组分和co2组分质量守恒方程表示为：
[0027][0028]
所述水组分质量守恒方程表示为：
[0029][0030]
所述油水三相运动方程表示为：
[0031][0032][0033]
所述气相运动方程表示为：
[0034]
[0035][0036]
所述辅助方程表示为：
[0037][0038]
zi＝noxi+ngyi，no+ng＝1
ꢀꢀꢀ
(8)
[0039]
so+sg+sw＝1
ꢀꢀꢀ
(9)
[0040]
pw＝p
g-p
cgw
(sw)，po＝p
g-p
cgo
(sw,so)
ꢀꢀꢀ
(10)
[0041][0042]
式中，ρo、ρg和ρw分别表示油气水三相密度；μo、μg和μw分别表示油气水三相黏度；vo、vg和vw分别表示油气水三相流动速度；po、pg和pw分别表示油气水三相压力；so、sg和sw分别表示油气水三相饱和度；xi、yi和zi分别表示烃类组分i或者co2组分在油相、气相以及油气两相中摩尔分数；no和ng表示油气相的摩尔分数；nc表示组分个数；φ
l
表示考虑热流固耦合效应的拉格朗日孔隙度；qi和qw表示烃类组分、co2组分和水组分源汇项；d
eff,i
表示气相中组分i的分子扩散系数；q
ads,i
表示气相中组分i的摩尔吸附量； k
ro
、k
rg
和k
rw
分别表示油气水三相相对渗透率；k为绝对渗透率；ka为表观渗透率；b为滑脱因子；αk和kn表示稀薄系数和knudsen数；p
cgw
和p
cgo
分别表示气水和油气的毛管力；表示逸度系数；l和v表示液相和气相；
·
表示点乘符号，
▽
表示梯度算子，t为时间。
[0043]
进一步地，所述超临界co
2-岩石热量传输模型表示为：
[0044][0045]
式中，ρf、cf、vf和pf分别表示当前位置流体相的密度、比热容、渗流速度和压力；t为温度；p、c和λ为饱和流体多孔介质等效密度、比热容和热传导系数；αs表示岩石骨架线膨胀数；αm表示热流体质量交换系数；kd表示岩石体积模量；εv表示体积应变；t0表示初始温度；q表示热源汇项；
·
表示点乘符号，
▽
表示梯度算子，t为时间。
[0046]
进一步地，所述建立岩体变形和裂缝扩展模型，包括：
[0047]
确定岩体变形平衡方程：
[0048]
▽
·
σ+f＝0
ꢀꢀꢀ
(13)
[0049]
考虑流体压力和热应力影响，确定应力增量方程：
[0050]
dσ＝c:dε-αdpi-3αskddti
ꢀꢀꢀ
(14)
[0051]
确定裂缝扩展准则方程：
[0052][0053]
式中，σ为总应力张量；ε为应变张量；c为弹性刚度矩阵；α为比奥系数；i为单位张量；f为体积力张量；tn、t
t
和ts分别表示内聚区单元垂向和两个剪切方向上的牵引力；和分别表示内聚区单元垂向和两个剪切方向上临界破坏牵引力；《
·
》表示运算符，即《*》＝(*+|*|)/2；
·
表示点乘符号，
▽
表示梯度算子。
[0054]
进一步地，所述基于固定应力分裂方法，构建温度场-渗流场-应力场的迭代耦合模型，包括：
[0055]
同一时步内热流固迭代过程中假设总体积应力不发生改变，即
[0056][0057]
确定体积应变更新方程：
[0058][0059]
式中，σv表示总体积应力；l表示热流固迭代次数；表示热流固迭代 l+1次时的总体积应力；表示热流固迭代l次时的总体积应力；表示热流固迭代l+1次的体积应变；表示热流固迭代l次的体积应变；p
l+1
表示热流固迭代l+1次的流体压力；p
l
表示热流固迭代l次的流体压力；t
l+1
表示热流固迭代l次的温度；t
l
表示热流固迭代l次的温度；α表示比奥系数；αs表示岩石骨架线膨胀数；kd表示岩石体积模量。
[0060]
本发明至少具有以下技术效果：
[0061]
本发明公开了一种页岩油藏超临界co2压裂精细模拟方法，综合考虑 co2浸泡对岩石断裂韧度的影响、超临界co2与原油流动传热对其相态变化的影响以及页岩储层油气特殊流动机理，此方法填补了超临界co2压裂数值模拟技术的空缺，具有计算精度高和可靠性高的优势。
附图说明
[0062]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0063]
图1示出了根据本发明实施例的一种页岩油藏超临界co2压裂精细模拟方法的流程示意图；
[0064]
图2示出了根据本发明实施例的一种页岩油藏超临界co2压裂精细模拟方法的岩石断裂韧度测试试样和三点弯曲试验示意图；
[0065]
图3示出了根据本发明实施例的一种页岩油藏超临界co2压裂精细模拟方法的页岩油藏几何模型与天然裂缝分布示意图；
[0066]
图4(a)示出了根据本发明实施例的一种页岩油藏超临界co2压裂精细模拟方法的几何模型初始网格示意图；
[0067]
图4(b)示出了根据本发明实施例的一种页岩油藏超临界co2压裂精细模拟方法的插入全局内聚区单元示意图；
[0068]
图4(c)示出了根据本发明实施例的一种页岩油藏超临界co2压裂精细模拟方法的内聚区单元网格节点自由度示意图；
[0069]
图5(a)示出了根据本发明实施例的一种页岩油藏超临界co2压裂精细模拟方法的水力压裂后的页岩油藏压裂缝网分布结果图；
[0070]
图5(b)示出了根据本发明实施例的一种页岩油藏超临界co2压裂精细模拟方法的
超临界co2压裂后的页岩油藏压裂缝网分布结果图。
具体实施方式
[0071]
为使本领域技术人员更好的理解本发明的技术方案，下面结合附图和具体实施方式对本发明作详细说明。本文中所描述的各个步骤，如果彼此之间没有前后关系的必要性，则本文中作为示例对其进行描述的次序不应视为限制，本领域技术人员应知道可以对其进行顺序调整，只要不破坏其彼此之间的逻辑性导致整个流程无法实现即可。
[0072]
本发明实施例提供一种页岩油藏超临界co2压裂精细模拟方法，请参阅图1所示，是该方法的流程图。
[0073]
该方法包括步骤(1)-(7)，具体如下：
[0074]
(1)开展超临界co2浸泡后页岩断裂韧度测试实验，获取岩石断裂韧度参数。
[0075]
在一些实施例中，所述开展超临界co2浸泡后页岩断裂韧度测试实验具体步骤包括：
[0076]
①
将页岩制作成多个半圆形试件，沿圆心将岩石切割得到短槽；
[0077]
②
将页岩试件浸泡在超临界co2，浸泡时间可取为2h，4h，8h，12h 和24h；
[0078]
③
试件浸泡结束后，开展三点弯曲试验，获得不同浸泡时间下页岩断裂韧度结果，岩石试样和三点弯曲试验示意图如图2所示。
[0079]
(2)获取天然裂缝分布情况，建立页岩油藏几何模型，采用非结构性剖分方法对含天然裂缝的几何模型进行网格划分，以保证天然裂缝位于网格边界。
[0080]
在一些实施例中，所述获取天然裂缝分布情况的方法为，根据天然裂缝统计信息，采用随机方法生成天然裂缝分布，页岩油藏几何模型与天然裂缝分布示意图如图3所示；
[0081]
所述对含天然裂缝的几何模型进行采用非结构性网格剖分方法为，二维模型采用三角形单元剖分，三维模型采用四面体单元剖分，且以天然裂缝为内边界约束，获得匹配型网格。
[0082]
(3)插入全局内聚区单元，其中天然裂缝所在网格边界单元标记为天然裂缝内聚区单元，其余网格边界单元标记为基质内聚区单元，并根据岩石断裂韧度参数测试结果赋予内聚区单元属性。
[0083]
在一些实施例中，所述插入全局内聚区单元的方法为，在初始网格的基础上，如图4(a)所示，对每个单元的所有边界进行循环插入，每条边界拆分成两边组成一个内聚区单元，如图4(b)所示，初始网格节点自由度包括压力和温度，更新后网格节点自由度为位移值，如图4(c)所示。
[0084]
所述根据岩石断裂韧度参数测试结果赋予内聚区单元属性方法为，利用岩石断裂韧度与断裂能之间的关系，求解得到岩石断裂能，将其赋给内聚区单元断裂能属性。
[0085]
(4)基于离散裂缝模型，建立超临界co2多相流动组分模型，采用有限体积法对控制方程进行数值离散化，应用牛顿迭代法对离散得到的流体压力、相饱和度和组分摩尔浓度非线性方程组进行求解；
[0086]
在一些实施例中，所述建立超临界co2多相流动组分模型包括烃类组分、 co2组分和水组分的质量守恒方程、油气水三相的运动方程及辅助方程，具体方程为：
[0087]
烃类组分和co2组分质量守恒方程，
[0088][0089]
水组分质量守恒方程，
[0090][0091]
油水三相运动方程，
[0092][0093][0094]
气相运动方程，
[0095][0096][0097]
辅助方程，
[0098][0099]
zi＝noxi+ngyi，no+ng＝1
ꢀꢀꢀ
(8)
[0100]
so+sg+sw＝1
ꢀꢀꢀ
(9)
[0101]
pw＝p
g-p
cgw
(sw)，po＝p
g-p
cgo
(sw,so)
ꢀꢀꢀ
(10)
[0102][0103]
式中，ρo、ρg和ρw分别表示油气水三相密度；μo、μg和μw分别表示油气水三相黏度；vo、vg和vw分别表示油气水三相流动速度；po、pg和pw分别表示油气水三相压力；so、sg和sw分别表示油气水三相饱和度；xi、yi和zi分别表示烃类组分i或者co2组分在油相、气相以及油气两相中摩尔分数；no和ng表示油气相的摩尔分数；nc表示组分个数；φ
l
表示考虑热流固耦合效应的拉格朗日孔隙度；qi和qw表示烃类组分、co2组分和水组分源汇项；d
eff,i
表示气相中组分i的分子扩散系数；q
ads,i
表示气相中组分i的摩尔吸附量； k
ro
、k
rg
和k
rw
分别表示油气水三相相对渗透率；k为绝对渗透率；ka为表观渗透率；b为滑脱因子；αk和kn表示稀薄系数和knudsen数；p
cgw
和p
cgo
分别表示气水和油气的毛管力；表示逸度系数；l和v表示液相和气相；
·
表示点乘符号，
▽
表示梯度算子，t为时间。
[0104]
(5)基于离散裂缝模型，考虑热对流和热传导两种方式，建立超临界 co2-岩石热量传输模型，采用有限体积法对控制方程进行数值离散化，采用牛顿迭代法对离散得到的温度场方程进行求解。
[0105]
在一些实施例中，所述建立的超临界co2-岩石热量传输模型为：
[0106]
考虑热对流和热传导能量守恒方程，
[0107]
[0108]
式中，ρf、cf、vf和pf分别表示当前位置流体相的密度、比热容、渗流速度和压力；t为温度；p、c和λ为饱和流体多孔介质等效密度、比热容和热传导系数；αs表示岩石骨架线膨胀数；αm表示热流体质量交换系数；kd表示岩石体积模量；εv表示体积应变；t0表示初始温度；q表示热源汇项；
·
表示点乘符号，
▽
表示梯度算子，t为时间。
[0109]
(6)建立岩体变形和裂缝扩展模型，考虑流体压力和热应力的影响，采用有限元方法对控制方程进行数值离散化，采用牛顿迭代法对离散得到的应力场方程进行求解。
[0110]
在一些实施例中，建立岩体变形和裂缝扩展模型为：
[0111]
岩体变形平衡方程，
[0112]
▽
·
σ+f＝0
ꢀꢀꢀ
(13)
[0113]
考虑流体压力和热应力影响，应力增量方程：
[0114]
dσ＝c:dε-αdpi-3αskddti
ꢀꢀꢀ
(14)
[0115]
裂缝扩展准则方程：
[0116][0117]
式中，σ为总应力张量；ε为应变张量；c为弹性刚度矩阵；α为比奥系数；i为单位张量；f为体积力张量；tn、t
t
和ts分别表示内聚区单元垂向和两个剪切方向上的牵引力；和分别表示内聚区单元垂向和两个剪切方向上临界破坏牵引力；《
·
》表示运算符，即《*》＝(*+|*|)/2；
·
表示点乘符号，
▽
表示梯度算子。
[0118]
(7)基于固定应力分裂算法，构建温度场-渗流场-应力场的迭代耦合算法，采用mpi(message passing interface)并行计算求解各物理场，循环迭代至三场耦合模型收敛，根据内聚区损伤演化准则判断裂缝是否扩展，从而更新裂缝网格信息，最终求解获得裂缝分布形态，实现页岩油藏超临界co2压裂精细数值模拟，模拟得到水力压裂和超临界co2压裂缝网分布结果分别如图5(a)和图5(b)所示，由图可知超临界co2压裂改造缝网体积大于水力压裂改造缝网体积。
[0119]
在一些实施例中，所述基于固定应力分裂算法，构建温度场-渗流场-应力场的迭代耦合算法，具体实现方法为，同一时步内热流固迭代过程中假设总体积应力不发生改变，即
[0120][0121]
体积应变更新方程：
[0122][0123]
式中，σv表示总体积应力；l表示热流固迭代次数；表示热流固迭代 l+1次时的总体积应力；表示热流固迭代l次时的总体积应力；表示热流固迭代l+1次的体积应变；表示热流固迭代l次的体积应变；p
l+1
表示热流固迭代l+1次的流体压力；p
l
表示热流固迭代l次的流体压力；t
l+1
表示热流固迭代l+1次的温度；t
l
表示热流固迭代l次的温度；α表示比奥系数；αs表示岩石骨架线膨胀数；kd表示岩石体积模量。
[0124]
以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实
施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。