稳定三维水力裂缝扩展计算并加速的自适应时间步算法

1.本发明涉及水力压裂、酸化压裂过程中三维水力裂缝扩展数值模拟技术领域，具体涉及一种稳定三维水力裂缝扩展计算并加速的自适应时间步算法。


背景技术：

2.随着勘探开发的进步，水力压裂、酸化压裂已经成为油气勘探领域中必不可少的技术手段。由于油气藏条件的复杂性，水力裂缝扩展数值模拟技术是水力压裂、酸化压裂优化设计的重要手段。随着技术的进步，三维水力裂缝扩展模型业已成为水力裂缝扩展数值模拟的主要发展方向。同时，在水力裂缝扩展过程中，注入的流体会由水力裂缝壁面、天然裂缝等多尺度流动介质大量滤失，降低流体造缝效率，准确计算流体滤失对水力裂缝扩展数值模拟至关重要。前期研究表明，为了保证水力裂缝扩展计算的稳定性和精度，在扩展初期需要采用极小的时间步长(一般为0.01～0.1s)，在扩展后期则可适当放宽。而整个水力压裂、酸压压裂施工时间动辄以数小时、甚至以数天计，在模型基质域为三维模型条件下，如此小的时间步会导致裂缝宽度及流体滤失计算过程中存在非常高的计算成本。同时，为了保证裂缝扩展过程中注入流体的质量守恒，还需要进一步对裂缝内流场、基质渗流场、滤失量和裂缝宽度(受到裂缝内流场压力计算影响)开展迭代计算。上述过程使得计算成本过大，普通的计算资源往往难以实现上述过程(或需要花费非常长的时间)。


技术实现要素：

3.为了克服现有技术中的问题，本发明提供稳定三维水力裂缝扩展计算并加速的自适应时间步算法。
4.本发明解决上述技术问题所提供的技术方案是：稳定三维水力裂缝扩展计算并加速的自适应时间步算法，包括以下步骤：
5.步骤s1、获取当前时间步数ti时水力裂缝区域ω
hf
及水力裂缝长度l
hf
、水力裂缝高度h
hf
、水力裂缝宽度分布w(x,z)；
6.步骤s2、计算切割区域大小d
cut
和基质岩体区域各方向保留网格数并确定对应网格保留区间，对三维水力裂缝扩展模型基质岩体区域进行切割；
7.步骤s3、完成基质岩体区域切割后，将水力裂缝所在平面扩充为网格，并计算扩充网格宽度和水力裂缝区域内等效渗透率；
8.步骤s4、根据储层基质平均渗透率选择预选时间步t
ini
；
9.步骤s5、根据预选时间步t
ini
计算水力裂缝在未来注入四个时间步内的流体压力分布p
hf
(x,z)；
10.步骤s6、在水力裂缝尖端网格找到流体压力最大值所对应的网格，定义该网格压力为p
tip,max
，则对应的在未来注入四个时间步裂缝尖端最大流体压力计算结果分别为p
tip,max
(1)、p
tip,max
(2)、p
tip,max
(3)、p
tip,max
(4)；
11.步骤s7、计算水力裂缝扩展临界宽度，并根据水力裂缝扩展临界宽度计算对应的
水力裂缝扩展临界压力p
ic
；
12.步骤s8、根据水力裂缝扩展临界压力p
ic
和p
tip,max
(3)、p
tip,max
(4)判断预选时间步是否符合要求，若不符合要求，则计算新的预选时间步并重复步骤s5-s8，直到预选时间步符合要求，再进行下一步骤；
13.步骤s9、根据当前时间步由于水力裂缝扩展导致的水力裂缝面积增量确定下一时间步是否需要重新计算自适应时间步；若需要重新计算自适应时间步则重复步骤s1-s8获得新的自适应时间步；
14.步骤s10、重复步骤s1-s9，则可实现在保证三维水力裂缝扩展计算稳定性的条件下加速的自适应时间步计算。
15.进一步的技术方案是，所述步骤s2中切割区域大小d
cut
的计算公式为：
[0016][0017]
式中：d
cut
为切割区域大小，m；a
hf
为水力裂缝面积，m2；km为储层基质平均渗透率，m2，p
hf
(x,z)为水力裂缝区域ω
hf
内流体压力，pa；pe为储层流体压力， pa；μ为注入流体粘度，pa
·
s；q为注入排量，m3/s；t
max
为自适应时间步长范围上限，s。
[0018]
进一步的技术方案是，所述s2中基质岩体区域内以注入点为原点，各方向保留网格数的计算公式为：
[0019]rx
＝π(d
cut
/δx)
[0020][0021]rz
＝∏(d
cut
/δz)
[0022]
式中：d
cut
为切割区域大小，m；r
x
为基质岩体区域中x方向保留网格数； ry为基质岩体区域中y方向保留网格数；rz为基质岩体区域中z方向保留网格数；δx为x方向网格长度，m；δy为y方向网格长度，m，此处为保证计算精度，y方向一般采用对数网格；δz为z方向网格长度，m；
[0023]
若以注入点(水力裂缝中心点)为原点，设注入点网格位置编号为(0,0,0)，则切割后以不同方向网格保留区间为：
[0024]
x方向：
[0025]
[-π(l
hf
/2δx)-r
x
，∏(l
hf
/2δx)+r
x
]
[0026]
y方向：
[0027]
[-ry，ry]
[0028]
z方向：
[0029]
[π(h
hf
/2δz)-rz，π(h
hf
/2δz)+rz]
[0030]
式中：r
x
为基质岩体区域中x方向保留网格数；ry为基质岩体区域中y方向保留网格数；rz为基质岩体区域中z方向保留网格数；δx为x方向网格长度， m；δz为z方向网格长度，m；l
hf
为水力裂缝长度，m；h
hf
为水力裂缝高度，m。
[0031]
进一步的技术方案是，所述步骤s3中扩充网格宽度的计算公式为：
[0032][0033]
式中：y
kc
为扩充网格宽度，m；δx为x方向网格长度，m；δz为z方向网格长度，m；a
hf
为水力裂缝面积，m2；w(x,z)为水力裂缝宽度分布，m。
[0034]
进一步的技术方案是，所述步骤s3中水力裂缝区域内等效渗透率的计算公式为：
[0035][0036][0037]
式中：y
kc
为扩充网格宽度，m；k
x,hf
(x,z)为水力裂缝区域ω
hf
内x方向等效渗透率，m2；k
y,hf
(x,z)为水力裂缝区域ω
hf
内y方向等效渗透率，m2；w(x,z)为水力裂缝宽度分布，m。
[0038]
进一步的技术方案是，所述步骤s4中预选时间步t
ini
的选取如下：
[0039]
t
ini
＝0.05skm《1md
[0040]
t
ini
＝0.1s1md≤km《5md
[0041]
t
ini
＝0.5s5md≤km《10md
[0042]
t
ini
＝1s10md≤km[0043]
式中：km为储层基质平均渗透率；t
ini
为预选时间步，s。
[0044]
进一步的技术方案是，所述步骤s5中流体压力分布的计算公式为：
[0045][0046]
式中：μ为注入流体粘度，pa
·
s；k
x
为模型x方向渗透率，m2；ky为模型y方向渗透率，m2；kz为模型z方向渗透率，m2；φ为孔隙度，无量纲；pm为流体压力，pa；t为计算时间，s；qf为由于注入产生的源项，s-1
。
[0047]
进一步的技术方案是，所述步骤s7中水力裂缝扩展临界宽度的计算公式为：
[0048][0049]
式中：e为储层岩石弹性模量，pa；v为储层岩石泊松比，无量纲；k
ic
为储层岩石断裂韧性，pa
·m1/2
；w
tip
为水力裂缝扩展临界宽度，m；δx为x方向网格长度，m；
[0050]
所述水力裂缝扩展临界压力p
ic
的计算公式为：
[0051][0052]
式中：e为储层岩石弹性模量，pa；v为储层岩石泊松比，无量纲；w
tip
为水力裂缝扩展临界宽度，m；p
ic
为水力裂缝扩展临界压力，pa；δz为z方向网格长度，m；σ
oh
为最小水平主应力，pa。
[0053]
进一步的技术方案是，所述步骤s8中p
ic
和p
tip,max
(3)、p
tip,max
(4)的关系若满足下式，则判断当前的预选时间步符合要求，t
ini
则为计算所得的自适应时间步：
[0054]
p
tip,max
(3)≤p
ic
≤p
tip,max
(4)
[0055]
若不满足，则根据p
tip,max
(3)或p
tip,max
(4)与p
ic
的相对关系计算新的预选时间步，
[0056]
若p
ic
≤p
tip,max
(3)，则取：
[0057][0058]
若p
tip,max
(4)＜p
ic
，则取：
[0059][0060]
式中：p
ic
为水力裂缝扩展临界压力，pa；t
ini
为预选时间步，s；ti'
ni
为新的预选时间步，s。
[0061]
进一步的技术方案是，所述步骤s9中自适应时间步重新计算的条件为水力裂缝扩展导致的水力裂缝面积变化符合以下条件：
[0062]ahf
(ti)/a
hf
(tn)≥1.1
[0063]
式中：a
hf
为水力裂缝面积，m2；ti为当前时间步数；tn为获取当前自适应时间步的时间步数。
[0064]
本发明具有以下有益效果：本发明基于水力裂缝扩展准则，以水力裂缝尖端压力作为判断依据，在需要的计算步实时修正计算时间步，实现时间步取值为能精确计算水力裂缝扩展的最大时间步长，以此保证水力裂缝扩展计算的稳定性并实现计算加速。
附图说明
[0065]
图1是三维水力裂缝扩展模型示意图；
[0066]
图2是基质岩体区域切割示意图；
[0067]
图3是裂缝平面扩充为网格示意图；
[0068]
图4是水力裂缝宽度分布图；
[0069]
图5是水力裂缝区域等效渗透率分布图；
[0070]
图6是第四个时间步时水力裂缝内流体压力分布图；
[0071]
图7是自适应时间步计算结果图。
具体实施方式
[0072]
下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0073]
假设在一l
×h×
w尺寸的三维基质岩体区域内计算水力裂缝扩展并耦合流体在水力裂缝内流动、流体在基质岩体内的渗流及流体从水力裂缝向基质岩体的滤失。其主要控制方程包括：
[0074]
①
水力裂缝内流动控制方程：
[0075][0076]
②
基质岩体内流动控制方程：
[0077][0078]
式中：w为水力裂缝宽度，m；p
hf
为水力裂缝内流体压力，pa；v
nf
为水力裂缝内流体滤失速度，m/s；μ为流体在y方向的流速，m/s；t为注酸时间，s； c
t
为油藏综合压缩系数，pa-1
；k
x
为模型x方向渗透率，m2；ky为模型y方向渗透率，m2；kz为模型z方向渗透率，m2；φ为孔隙度，无量纲；pm为流体压力， pa；qf为由于注入产生的源项，s-1
。
[0079]
③
流体由水力裂缝向基质岩体滤失的控制方程：
[0080][0081]
式中：v
hf
为水力裂缝到基质岩体的滤失速率，m/s；ky为模型y方向渗透率， m2；μ为注入流体粘度，pa
·
s。
[0082]
④
水力裂缝宽度计算方程：
[0083]
x方向：
[0084][0085]
y方向：
[0086][0087]
z方向：
[0088][0089]
式中：u
x
、uy、uz分别为基质岩体在x、y、z方向的位移，m/s；c为基质岩体弹性常数，下标代表应力正、切方向，无量纲；α为孔弹性系数，一般取0.6-0.8； ps为孔隙压力，mpa。
[0090]
⑤
水力裂缝扩展判据：
[0091][0092]
式中：e为储层岩石弹性模量，pa；v为储层岩石泊松比，无量纲；k
ic
为储层岩石断裂韧性，pa
·m1/2
；w
tip
为水力裂缝尖端宽度，m；δx为水力裂缝长度方向单元格长度，m；
[0093]
如图1所示，本发明的稳定三维水力裂缝扩展计算并加速的自适应时间步算法，包括以下步骤：
[0094]
步骤1、假设目前为计算时间步数ti，此时裂缝所占的区域为ω
hf
，裂缝长度为l
hf
，裂缝高度为h
hf
，裂缝宽度分布为w＝f(x,z)；
[0095]
如图1所示的三维水力裂缝扩展模型，确定长度l方向的网格数量为numx，宽度w方
向的网格数量为numy，高度h方向的网格数量为numz；
[0096]
步骤2、计算切割区域大小d
cut
和基质岩体区域各方向保留网格数并确定对应网格保留区间，再对三维水力裂缝扩展模型基质岩体区域进行切割；
[0097]
切割区域大小d
cut
的计算公式如下：
[0098][0099]
式中：d
cut
为切割区域大小，m；a
hf
为水力裂缝面积，m2；km为储层基质平均渗透率，m2，p
hf
(x,z)为水力裂缝区域ω
hf
内流体压力，pa；pe为储层流体压力， pa；μ为注入流体粘度，pa
·
s；q为注入排量，m3/s；t
max
为自适应时间步长范围上限，s；
[0100]
基质岩体区域内以注入点为原点，各方向保留网格数的计算公式为为：
[0101]rx
＝π(d
cut
/δx)
ꢀꢀ
(9)
[0102][0103]rz
＝π(d
cut
/δz)
ꢀꢀ
(11)
[0104]
式中：r
x
为基质岩体区域中x方向保留网格数；ry为基质岩体区域中y方向保留网格数；rz为基质岩体区域中z方向保留网格数；δx为x方向网格长度， m；δy为y方向网格长度，m，此处为保证计算精度，y方向一般采用对数网格；δz为z方向网格长度，m；l
hf
为水力裂缝长度，m；h
hf
为水力裂缝高度，m；
[0105]
若以注入点(水力裂缝中心点)为原点，设注入点网格位置编号为(0,0,0)，则切割后以不同方向网格保留区间为：
[0106]
x方向：
[0107]
[-π(l
hf
/2δx)-r
x
，∏(l
hf
/2δx)+r
x
]
[0108]
y方向：
[0109]
[-ry，ry]
[0110]
z方向：
[0111]
[π(h
hf
/2δz)-rz，π(h
hf
/2δz)+rz]
[0112]
式中：r
x
为基质岩体区域中x方向保留网格数；ry为基质岩体区域中y方向保留网格数；rz为基质岩体区域中z方向保留网格数；δx为x方向网格长度， m；δz为z方向网格长度，m；l
hf
为水力裂缝长度，m；h
hf
为水力裂缝高度，m；
[0113]
步骤3、完成基质岩体区域切割后，为了避免裂缝内流场、基质渗流场、滤失量的迭代，将水力裂缝所在平面扩充为一层网格，采用等效渗透率、水力裂缝孔隙度表征水力裂缝流动能力，将裂缝内流动、基质渗流、滤失流动三个过程化为三维等效渗流统一计算，从而计算扩充网格宽度和水力裂缝区域内等效渗透率；
[0114]
扩充网格宽度一般取水力裂缝的平均宽度，即扩充网格宽度的计算公式为：
[0115]
[0116]
式中：y
kc
为扩充网格宽度，m；δx为x方向网格长度，m；δz为z方向网格长度，m；a
hf
为水力裂缝面积，m2；w(x,z)为水力裂缝宽度分布，m；
[0117]
此时忽略短时间(≤t
max
)内流动的裂缝宽度变化，采用基质岩体内三维渗流方程一体化计算水力裂缝、基质岩体及两者之间的流体流动。即：
[0118][0119]
式中：μ为注入流体粘度，pa
·
s；k
x
为模型x方向渗透率，m2；ky为模型y方向渗透率，m2；kz为模型z方向渗透率，m2；φ为孔隙度，无量纲；pm为流体压力，pa；t为计算时间，s；qf为由于注入产生的源项，s-1
；
[0120]
水力裂缝区域内等效渗透率根据立方定律和水力裂缝宽度计算，即：
[0121][0122][0123]
式中：y
kc
为扩充网格宽度，m；k
x,hf
(x,z)为水力裂缝区域ω
hf
内x方向等效渗透率，m2；k
y,hf
(x,z)为水力裂缝区域ω
hf
内y方向等效渗透率，m2；w(x,z)为水力裂缝宽度分布，m；
[0124]
步骤4、根据储层基质平均渗透率设置自适应时间步寻优计算时间，步长范围下限t
min
一般取0.01s，步长范围上限一般取t
max
；具体如下：
[0125]
t
ini
＝0.05skm《1md
[0126]
t
ini
＝0.1s1md≤km《5md
[0127]
t
ini
＝0.5s5md≤km《10md(19)
[0128]
t
ini
＝1s10md≤km[0129]
式中：km为储层基质平均渗透率；t
ini
为预选时间步，s；
[0130]
步骤5、根据建立的切割后区域模型，选取合适的预选时间步，求解式(16)获得水力裂缝在未来注入四个时间步内的流体压力分布p
hf
(x,z)；
[0131]
步骤6、在水力裂缝尖端网格找到流体压力最大值所对应的网格，定义该网格压力为p
tip,max
，则对应的在未来注入四个时间步裂缝尖端最大流体压力计算结果分别为p
tip,max
(1)、p
tip,max
(2)、p
tip,max
(3)、p
tip,max
(4)；
[0132]
步骤7、为了同时保证水力裂缝计算的扩展精度并最小化计算成本，时间步内裂缝尖端网格内的压力变化上限必须小于且接近水力裂缝扩展的临界压力；基于式(7)中计算所得的水力裂缝扩展临界宽度，并根据水力裂缝扩展临界宽度计算对应的水力裂缝扩展临界压力p
ic
；
[0133][0134]
式中：e为储层岩石弹性模量，pa；v为储层岩石泊松比，无量纲；w
tip
为水力裂缝扩展临界宽度，m；p
ic
为水力裂缝扩展临界压力，pa；δz为z方向网格长度，m；σ
oh
为最小水平主应力，pa；
[0135]
步骤8、根据水力裂缝扩展临界压力p
ic
和p
tip,max
(3)、p
tip,max
(4)判断预选时间步是
否符合要求，若不符合要求，则计算新的预选时间步并重复步骤5-8，直到时间步符合要求，再进行下一步骤；
[0136]
若满足下式，则判断当前的预选时间步符合要求，t
ini
则为计算所得的自适应时间步：
[0137]
p
tip,max
(3)≤p
ic
≤p
tip,max
(4)
ꢀꢀ
(21)
[0138]
若不满足，则根据p
tip,max
(3)或p
tip,max
(4)与p
ic
的相对关系计算新的预选时间步；
[0139]
若p
ic
≤p
tip,max
(3)，则取：
[0140][0141]
若p
tip,max
(4)≤p
ic
，则取：
[0142][0143]
式中：p
ic
为临界压力，pa；t
ini
为预选时间步，s；ti'
ni
为新的预选时间步，s；
[0144]
其中采用该时间步长，结合式(1)～(7)则可进行主程序迭代计算；
[0145]
步骤9、根据当前时间步由于水力裂缝扩展导致的水力裂缝面积增量确定下一时间步是否需要重新计算自适应时间步；若需要重新计算自适应时间步则重复步骤1-8获得新的自适应时间步；
[0146]
其中自适应时间步重新计算的条件为水力裂缝扩展导致的水力裂缝面积变化符合以下条件：
[0147]ahf
(ti)/a
hf
(tn)≥1.1
ꢀꢀ
(24)
[0148]
式中：a
hf
为水力裂缝面积，m2；ti为当前时间步数；tn为获取当前自适应时间步的时间步数；
[0149]
步骤10、重复步骤s1-s9，则可实现在保证三维水力裂缝扩展计算稳定性的条件下加速的自适应时间步计算。
[0150]
实施例
[0151]
1、假设在一400
×
200
×
50m尺寸的三维基质岩体区域，长度l方向的网格数量numx＝200，宽度w方向的网格数量为numy＝21(w方向一般取对数网格)，高度h方向的网格数量为numz＝100。水力裂缝长度为l
hf
＝72m，水力裂缝高度为h
hf
＝64m，水力裂缝宽度w(x,z)如图4所示。此时时间步数ti＝24。
[0152]
2、此时水力裂缝面积a
hf
＝4096m2，储层基质平均渗透率km＝10
×
10-15
m2，水力裂缝内最大流体压力p
hf
＝132
×
106pa，储层流体压力pe＝60
×
106pa，注入流体粘度μ＝0.025pa
·
s，注入排量q＝0.1m3/s。
[0153]
则根据式(8)计算得切割区域大小d
cut
为2.36m。
[0154]
则根据式(9)～(11)计算基质岩体区域各方向保留网格数为：r
x
＝2，ry＝ 7,rz＝2。
[0155]
本模型原点网格编号numx0＝101，numy0＝11，numz0＝51，根据式(12) ～(14)则可计算切割后不同方向保留网格的区间为x方向：[81，121]；y方向：[4，18]；z方向：[33，69]。
[0156]
模型切割后网格数量为41
×
15
×
37＝7770个，较原始模型200
×
100
×
21＝ 42000个减少了98％。
[0157]
3、完成基质岩体区域切割后，将水力裂缝所在平面扩充为一层网格，根据式(15)计算扩充网格宽度y
kc
为2.1
×
10-3
m。
[0158]
基于水力裂缝宽度分布w(x，z)及式(17)～(18)，则可计算水力裂缝区域内等效渗透率，计算结果如图5所示。
[0159]
4、储层基质平均渗透率km＝10md，则选取预选时间步t
ini
＝1s，
[0160]
5、根据式(16)求解得水力裂缝在未来注入四个时间步内的流体压力分布 p
hf
(x,z)如图6所示。
[0161]
6、遍寻水力裂缝尖端网格，则对应的在未来注入四个时间步裂缝尖端最大流体压力计算结果分别为p
tip,max
(1)＝1.2221
×
108pa，p
tip,max
(2)＝1.2563
×
108pa， p
tip,max
(3)＝1.2796
×
108pa，p
tip,max
(4)＝1.2967
×
108pa。模型中x/y方向网格长度δx＝δy＝2m，断裂韧性k
ic
＝7
×
105pa
·
m-1
，杨氏模量e＝50000
×
106pa，根据式(7)计算得w
tip
＝0.0011m；进一步根据式(20)计算水力裂缝扩展临界压力 p
ic
＝1.2832
×
108pa。则p
tip,max
(3)、p
tip,max
(4)及p
ic
之间得关系满足式(21)，则本预选时间步t
ini
＝1s。
[0162]
7、结合式(1)～(7)则可进行主程序迭代计算。同时，应当在适当的时间重新根据步骤1～5计算自适应时间步，根据式(24)可计算得变化时间步得n ＝1，则下个时间步需重新计算自适应时间步长。
[0163]
8、重复步骤1～7，则可实现在保证三维水力裂缝扩展计算稳定性的条件下加速的自适应时间步计算。
[0164]
以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，可利用上述揭示的技术内容作出些变动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。