一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法
1.本发明属于干热岩地热能开发的人工热储建造技术领域，具体涉及一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法。
背景技术：
2.随着经济的发展，传统化石能源的供应日益紧缺。为实现节能减排和碳中和目标，可再生能源在能源消费结构中的比重将提高。地热能作为一种储量丰富、分布广泛的可再生清洁能源，在实现碳中和目标中发挥着重要作用。干热岩是地热能最主要的存储介质，其开发主要基于使用增强型地热系统来提取内部热量。利用水力压裂等工程技术，在深部地下低渗透干热岩中形成人工储热结构，从而在长期内经济地产生可观的热能。
3.然而深部岩体完整，水力压裂困难，因此钻井往往选择在天然裂隙发育和断层带区域，这些区域地应力复杂，天然裂隙分布情况未知，导致水力压裂工程设计时裂隙扩展预测困难。如果水力压裂不能准确预测，可能造成压裂液漏失严重。室内试验受限于苛刻的条件无法有效开展，因此干热岩水力裂缝扩展数值模拟是有必要的。
4.高温深层地层进行水力压裂施工时，低温压裂液快速注入高温地层，压裂缝中的液体与地层快速进行热交换，地层岩石的力学性能因水冷冲击作用而劣化，并且因为冷却收缩而产生低温诱导热应力，两者共同作用可以有效抵消部分岩石起裂时阻力。然而，目前的多数水力压裂设计裂缝扩展模型尚未考虑到水冷冲击对岩石力学性能的劣化作用给裂缝扩展所带来的影响，存在一定的局限性。
技术实现要素：
5.本发明的目的是提供一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，从而准确的预测起裂压力及储层水力裂缝扩展形态。
6.本发明所采用的技术方案是，一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，具体按照以下步骤实施：
7.步骤1、收集干热岩地热储层的物理参数、拟注入的压裂液参数和工程参数；
8.步骤2、根据工程参数，利用comsol绘制拟模拟压裂的岩体形状，得到计算域，并利用matlab实现岩石非均匀性的weibull 分布；
9.步骤3、对计算域添加热-流-固多物理场，建立相互的耦合关系，形成热-流-固多物理场耦合模型；
10.步骤4、对热-流-固多物理场耦合模型进行全耦合计算，将结果导入matlab进行损伤变量计算；
11.步骤5、将损伤变量导入comsol中计算损伤对各参数的影响，并重复步骤4，直到该注入压力下损伤变量不再增加，即损伤稳定；
12.步骤6、增加注入压力p
inj
，重复步骤4、5进行循环计算，直至岩体完全破裂即可停止计算，即可确定压裂结束时裂隙扩展形态。
13.本发明的特点还在于，
14.步骤1中，物理参数包括：最大水平主应力、最小水平主应力、岩石弹性模量、岩石泊松比、岩石抗拉强度、岩石孔隙度、岩石密度、地层原始温度、岩石渗透率、岩石热膨胀系数、岩石导热系数、岩石比热容、岩石非均质系数、岩石传热系数；压裂液参数包括：压裂液黏度、压裂液密度、压裂液导热系数、压裂液比热容、压裂液温度；工程参数包括：井筒直径和井间距。
15.步骤3中，热-流-固多物理场耦合模型方程为：
16.固体力学方程：
17.式中，μ和λ为拉梅常数；α为biot系数；α
t
为岩石的热膨胀系数；k为岩石的排水体积模量，k＝(2μ+3λ)/3；ui和fi分别为体力和位移在i方向的分量，其中分别为体力和位移在i方向的分量，其中p
,i
和t
,i
为孔隙压力和温度在i方向的分量；-αp
,i
为有效孔隙压力作用项；-kα
t
t
,i
为热应力项；
18.流体流动方程：
19.式中，ρ
l
为流体的密度；εv为体积应变；φ为岩体孔隙率，其中φ0为初始孔隙度，φr为强应力状态下的残余孔隙度，α
φ
为孔隙度对应力的敏感系数，为平均有效应力；β
l
为流体的压缩系数；k为岩体的渗透率；μ
l
为流体的动力黏滞系数；为梯度算子；为孔隙压力梯度；qm为流体的源汇项；α
l
为流体的热体积膨胀系数；为应力作用项；为温度作用项；
20.传热方程：
21.式中，(ρc)
eff
＝ρ
scs
(1-φ)+ρ
lcl
φ和λ
eff
＝λs(1-φ)+λ
l
φ分别表征干热岩的等效体比热容和导热系数，其中λs和cs分别为岩体的导热系数和比热容，c
l
和λ
l
分别为流体的比热容和导热系数；q
t
为对流热通量，q
t
＝h(t
ext-t)，h为对流换热系数，t
ext
为外部温度；为热弹性阻尼项，为流速影响的热对流项。
22.水冷冲击作用对岩石力学参数的劣化：
23.弹性模量：e(δt)＝e(t0)(1.348-0.347e
0.00215δt
)
24.抗拉强度：f
t
(δt)＝f
t
(t0)(1.31-0.31e
0.002δt
)
25.式中，e(t0)和f
t
(t0)分别为初始弹性模量和抗拉强度；降温幅度δt＝t
0-t，t0为岩石压裂前的初始温度，t为岩石受水冷冲击后的温度。
26.步骤4中，全耦合计算方法为：将干热岩地热储层的物理参数和压裂液参数代入
热-流-固多物理场耦合模型中，上方和右方分别设置最大主应力和最小主应力的均布边界荷载，下方和左方设置为辊支撑，四周为不透水边界，圆孔为压力边界，四周为绝热边界，圆孔为对流热通量边界，外部温度设置为压裂液温度；定义岩石初始渗透压力为零，初始温度为储层温度，输入初始注入压力p
inj
，开始瞬态计算，得到温度、压力、应力和应变分布；
27.损伤变量计算方程为：
28.最大拉应力准则：f＝σ
t-f
t
(δt)
29.损伤变量d：
30.式中，f为拉应力σ
t
和抗拉强度f
t
(δt)的差值，其中抗拉强是关于降温幅度δt的函数；ε
t0
为极限拉应变；ε
t
为拉应变。
31.根据最大拉应力准则，当岩体单元所受拉应力大于抗拉强度时，损伤变量增加，则判定岩石开裂，此时的注入压力p
inj
即为岩体的起裂压力。
32.步骤5中，损伤对各参数的影响的方程包括：
33.弹性模量：e(d,t)＝(1-d)e(δt)；
34.渗透率：
35.导热系数：λs＝λ
s0
exp(d/α
dλ
)；
36.比热容：
37.biot系数：
38.式中，e(d,t)为岩体弹性模量，e(δt)表示水冷冲击作用对岩石弹性模量的影响；k0为无应力状态的初始渗透率，α
dk
为描述损伤诱导的渗透率影响因子；λ
s0
和c
s0
分别为岩体的初始导热系数和比热容；α
dλ
为损伤对导热系数的影响因子；α
dc
为损伤对比热容的影响因子；α0是初始biot系数。
39.本发明的有益效果是，针对干热岩水力压裂过程中低温压裂液注入高温地层引起的水冷冲击作用，建立了热-流-固多物理场耦合模型，提供了一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法；由于该方法充分考虑了低温压裂液注入高温地层引起的水冷冲击影响，可准确预测起裂压力及储层水力裂缝扩展形态，定量地评估水力压裂施工效果，为增强型地热系统的压裂设计提供理论基础和参考。
附图说明
40.图1是本发明一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法的流程示意图；
41.图2为本发明干热岩计算域及边界条件示意图；
42.图3为本发明热-流-固多物理场耦合的方式示意图；
43.图4为本发明干热岩水力压裂的开裂机理示意图；
44.图5为本发明实施例中考虑水冷冲击作用与未考虑水冷冲击作用的水力压裂裂缝形态对比图。
具体实施方式
45.下面结合附图和具体实施方式对本发明进行详细说明。
46.本发明一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，如图1所示，具体按照以下步骤实施：
47.步骤1、收集干热岩地热储层的物理参数、拟注入的压裂液参数和工程参数；
48.其中，物理参数包括：最大水平主应力、最小水平主应力、岩石弹性模量、岩石泊松比、岩石抗拉强度、岩石孔隙度、岩石密度、地层原始温度、岩石渗透率、岩石热膨胀系数、岩石导热系数、岩石比热容、岩石非均质系数、岩石传热系数；
49.压裂液参数包括：压裂液黏度、压裂液密度、压裂液导热系数、压裂液比热容、压裂液温度；
50.工程参数包括：井筒直径和井间距；
51.步骤2、根据工程参数，利用comsol绘制拟模拟压裂的岩体形状，得到计算域，并利用matlab实现岩石非均匀性的weibull 分布；
52.利用comsol绘制拟模拟压裂的岩体形状的具体步骤为：以坐标原点为中心生成井间距大小的正方形以及井筒直径大小的圆形，对两者进行布尔差集运算得到计算域，如图2所示。
53.利用matlab实现岩石非均匀性的weibull分布的具体步骤为：使用三角形有限元的非结构化网格(通过delaunay细分)对域进行离散化，将网格节点坐标导出；利用matlab生成等同于网格节点数的weibull分布随机数并与储层参数相乘，得到不同参数的二维矩阵形式并导出为txt文件；在comsol基于文件定义插值函数来对材料参数赋值，从而实现储层的参数在空间上的不均匀性；
54.步骤3、对计算域添加热-流-固多物理场，建立相互的耦合关系，如图3所示；
55.使用comsol中的多孔介质传热、达西定律和固体力学等物理场方程。分别激活多孔弹性和热膨胀接口以实现固体力学和流体流动、固体力学和传热方程之间的直接耦合。温度改变会改变流体性质进而影响流场的流动，还可以以热应力的形式引起固体变形，而流场的流动引发热对流反过来又会导致温度发生变化，且流体在裂缝里流动时会产生孔隙水压力，从而导致固体发生变形；同样，固体的变形直接改变孔隙度，也会导致流场的流动发生变化，同时由于热扩散传热又会造成温度场发生改变。
56.多物理场耦合方程为：
57.固体力学方程：
58.式中，μ和λ为拉梅常数；α为biot系数；α
t
为岩石的热膨胀系数；k为岩石的排水体积模量，k＝(2μ+3λ)/3；ui和fi分别为体力和位移在i方向的分量，其中
p
,i
和t
,i
为孔隙压力和温度在i方向的分量；-αp
,i
为有效孔隙压力作用项；-kα
t
t
,i
为热应力项。
59.流体流动方程：
60.式中，ρ
l
为流体的密度；εv为体积应变；φ为岩体孔隙率，与岩石所承受的应力状态有关，其中φ0为初始孔隙度，φr为强应力状态下的残余孔隙度，α
φ
为孔隙度对应力的敏感系数，为平均有效应力；β
l
为流体的压缩系数；k为岩体的渗透率；μ
l
为流体的动力黏滞系数；为梯度算子；为孔隙压力梯度；qm为流体的源汇项；α
l
为流体的热体积膨胀系数；为应力作用项；为温度作用项；
61.传热方程：
62.式中，(ρc)
eff
＝ρ
scs
(1-φ)+ρ
lcl
φ和λ
eff
＝λs(1-φ)+λ
l
φ分别表征干热岩的等效体比热容和导热系数，其中λs和cs分别为岩体的导热系数和比热容，c
l
和λ
l
分别为流体的比热容和导热系数；q
t
为对流热通量，q
t
＝h(t
ext-t)，其中h为对流换热系数，t
ext
为外部温度；为热弹性阻尼项，为流速影响的热对流项。
63.干热岩储层的高温和地应力导致水力压裂过程中流体物理性质发生显着变化。为了更好地模拟真实压裂情况，必须充分考虑不同温度和压力下物性参数的变化。从美国国家标准与技术研究所的物性数据库得到不同温度和压力下流体的密度ρ
l
、黏度μ
l
、导热系数λ
l
和比热容c
l
、可压缩系数β
l
等参数，然后导入comsol中通过插值函数定义流体材料参数。
64.高温岩体水力压裂试验表明，干热岩水力压裂过程中，由于流体
ꢀ‑
岩体温度差异引起的水冷冲击作用会降低岩体的力学性能，增强岩石的渗透性，加速热交换过程。对不同高温花岗岩水冷后的力学试验结果进行归一化和平均化，提出了水冷冲击作用对岩石力学参数的劣化程度的数学表达式：
65.弹性模量：e(δt)＝e(t0)(1.348-0.347e
0.00215δt
)
66.抗拉强度：f
t
(δt)＝f
t
(t0)(1.31-0.31e
0.002δt
)
67.式中，e(t0)和f
t
(t0)分别为初始弹性模量和抗拉强度；降温幅度δt＝t
0-t，t0为岩石压裂前的初始温度，t为岩石受水冷冲击后的温度；
68.步骤4、将对热-流-固多物理场耦合模型进行全耦合计算，将结果导入matlab进行损伤变量计算；
69.全耦合计算方法为：步骤1的干热岩地热储层的物理参数和压裂液参数代入热-流-固多物理场耦合模型中，上方和右方分别设置最大主应力和最小主应力的均布边界荷载，下方和左方设置为辊支撑。四周为不透水边界，圆孔为压力边界。四周为绝热边界，圆孔为对流热通量边界，外部温度设置为压裂液温度；定义岩石初始渗透压力为零，初始温度为
储层温度。输入初始注入压力p
inj
，开始瞬态计算，得到温度、压力、应力和应变分布；
70.损伤变量计算方程为：
71.最大拉应力准则：f＝σ
t-f
t
(δt)
72.损伤变量d：
73.式中，f为拉应力σ
t
和抗拉强度f
t
(δt)的差值，其中抗拉强是关于降温幅度δt的函数；ε
t0
为极限拉应变；ε
t
为拉应变。
74.根据最大拉应力准则，当岩体单元所受拉应力大于抗拉强度时，损伤变量增加，则判定岩石开裂，此时的注入压力p
inj
即为岩体的起裂压力；
75.步骤5、将损伤变量导入comsol中计算损伤对各参数的影响，并重复步骤4，直到该注入压力下损伤变量不再增加，即损伤稳定。
76.损伤对各参数的影响的方程包括：
77.弹性模量：e(d,t)＝(1-d)e(δt)；
78.渗透率：
79.导热系数：λs＝λ
s0
exp(d/α
dλ
)；
80.比热容：
81.biot系数：
82.式中，e(d,t)为岩体弹性模量，是损伤和温度的函数，其中 e(δt)表示水冷冲击作用对岩石弹性模量的影响；k0为无应力状态的初始渗透率，α
dk
为描述损伤诱导的渗透率影响因子；λ
s0
和c
s0
分别为岩体的初始导热系数和比热容；α
dλ
为损伤对导热系数的影响因子；α
dc
为损伤对比热容的影响因子；α0是初始biot系数。
83.步骤6、增加注入压力p
inj
，重复步骤4、5进行循环计算，直至岩体完全破裂即可停止计算，即可确定压裂结束时裂隙扩展形态。
84.干热岩水力压裂的起裂原理如图4所示。由于高温环境，干热岩在形成过程中始终保持热膨胀状态。当被水冷却时，岩石的快速收缩会导致显着的热应力。与孔隙水压力结合形成拉应力。由于热交换，围岩温度迅速下降，水冷冲击使岩石的抗拉强度显着下降。干热岩水力压裂过程中高温岩石的抗拉强度、冷收缩引起的热应力和孔隙水压力均发生动态变化。当岩石的拉应力达到其自身的抗拉强度时，井筒围岩破裂。干热岩的开裂是由水冷冲击引起的抗拉强度劣化、冷收缩引起的热应力和水力压裂中的孔隙水压力共同作用引起的。
85.实施例
86.以国内某增强型地热系统为例，利用表1所示的干热岩储层的主要参数模拟压裂结束后的裂缝形态。按如前步骤开展实例计算，并根据计算结果，输出压裂结束时裂缝形
态。同时，对比了考虑水冷冲击作用和未考虑考虑水冷冲击作用下的压裂裂缝形态的差异，如图5所示，可以看出，不考虑冷冲击作用会对起裂压力的估算产生较大误差。考虑水冷冲击作用的裂隙宽度和长度更大，裂隙分布更密集。
87.表1某干热岩储层的勘探参数
[0088][0089]
本发明针对干热岩水力压裂过程中低温压裂液注入高温地层引起的水冷冲击作用，建立了热-流-固多物理场耦合模型，提供了一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法。由于该方法充分考虑了低温压裂液注入高温地层引起的水冷冲击影响，可准确预测起裂压力及储层水力裂缝扩展形态，定量地评估水力压裂施工效果，为增强型地热系统的压裂设计提供理论基础和参考。
技术特征：
1.一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，其特征在于，具体按照以下步骤实施：步骤1、收集干热岩地热储层的物理参数、拟注入的压裂液参数和工程参数；步骤2、根据工程参数，利用comsol绘制拟模拟压裂的岩体形状，得到计算域，并利用matlab实现岩石非均匀性的weibull分布；步骤3、对计算域添加热-流-固多物理场，建立相互的耦合关系，形成热-流-固多物理场耦合模型；步骤4、对热-流-固多物理场耦合模型进行全耦合计算，将结果导入matlab进行损伤变量计算；步骤5、将损伤变量导入comsol中计算损伤对各参数的影响，并重复步骤4，直到该注入压力下损伤变量不再增加，即损伤稳定；步骤6、增加注入压力p
inj
，重复步骤4、5进行循环计算，直至岩体完全破裂即可停止计算，即可确定压裂结束时裂隙扩展形态。2.根据权利要求1所述的一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，其特征在于，所述步骤1中，物理参数包括：最大水平主应力、最小水平主应力、岩石弹性模量、岩石泊松比、岩石抗拉强度、岩石孔隙度、岩石密度、地层原始温度、岩石渗透率、岩石热膨胀系数、岩石导热系数、岩石比热容、岩石非均质系数、岩石传热系数；压裂液参数包括：压裂液黏度、压裂液密度、压裂液导热系数、压裂液比热容、压裂液温度；工程参数包括：井筒直径和井间距。3.根据权利要求1所述的一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，其特征在于，所述步骤3中，热-流-固多物理场耦合模型方程为：固体力学方程：式中，μ和λ为拉梅常数；α为biot系数；α
t
为岩石的热膨胀系数；k为岩石的排水体积模量，k＝(2μ+3λ)/3；u
i
和f
i
分别为体力和位移在i方向的分量，其中分别为体力和位移在i方向的分量，其中p
,i
和t
,i
为孔隙压力和温度在i方向的分量；-αp
,i
为有效孔隙压力作用项；-kα
t
t
,i
为热应力项；流体流动方程：式中，ρ
l
为流体的密度；ε
v
为体积应变；φ为岩体孔隙率，其中φ0为初始孔隙度，φ
r
为强应力状态下的残余孔隙度，α
φ
为孔隙度对应力的敏感系数，为平均有效应力；β
l
为流体的压缩系数；k为岩体的渗透率；μ
l
为流体的动力黏滞系数；为梯度算子；为孔隙压力梯度；q
m
为流体的源汇项；α
l
为流体的热体积膨胀系数；为应力作用项；为温度作用项；
传热方程：式中，(ρc)
eff
＝ρ
s
c
s
(1-φ)+ρ
l
c
l
φ和λ
eff
＝λ
s
(1-φ)+λ
l
φ分别表征干热岩的等效体比热容和导热系数，其中λ
s
和c
s
分别为岩体的导热系数和比热容，c
l
和λ
l
分别为流体的比热容和导热系数；q
t
为对流热通量，q
t
＝h(t
ext-t)，h为对流换热系数，t
ext
为外部温度；为热弹性阻尼项，为流速影响的热对流项；水冷冲击作用对岩石力学参数的劣化：弹性模量：e(δt)＝e(t0)(1.348-0.347e
0.00215δt
)抗拉强度：f
t
(δt)＝f
t
(t0)(1.31-0.31e
0.002
δ
t
)式中，e(t0)和f
t
(t0)分别为初始弹性模量和抗拉强度；降温幅度δt＝t
0-t，t0为岩石压裂前的初始温度，t为岩石受水冷冲击后的温度。4.根据权利要求3所述的一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，其特征在于，所述步骤4中，全耦合计算方法为：将干热岩地热储层的物理参数和压裂液参数代入热-流-固多物理场耦合模型中，上方和右方分别设置最大主应力和最小主应力的均布边界荷载，下方和左方设置为辊支撑，四周为不透水边界，圆孔为压力边界，四周为绝热边界，圆孔为对流热通量边界，外部温度设置为压裂液温度；定义岩石初始渗透压力为零，初始温度为储层温度，输入初始注入压力p
inj
，开始瞬态计算，得到温度、压力、应力和应变分布；损伤变量计算方程为：最大拉应力准则：f＝σ
t-f
t
(δt)损伤变量d：式中，f为拉应力σ
t
和抗拉强度f
t
(δt)的差值，其中抗拉强是关于降温幅度δt的函数；ε
t0
为极限拉应变；ε
t
为拉应变；根据最大拉应力准则，当岩体单元所受拉应力大于抗拉强度时，损伤变量增加，则判定岩石开裂，此时的注入压力p
inj
即为岩体的起裂压力。5.根据权利要求4所述的一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，其特征在于，所述步骤5中，损伤对各参数的影响的方程包括：弹性模量：e(d,t)＝(1-d)e(δt)；渗透率：导热系数：λ
s
＝λ
s0
exp(d/α
dλ
)；比热容：
biot系数：式中，e(d,t)为岩体弹性模量，e(δt)表示水冷冲击作用对岩石弹性模量的影响；k0为无应力状态的初始渗透率，α
dk
为描述损伤诱导的渗透率影响因子；λ
s0
和c
s0
分别为岩体的初始导热系数和比热容；α
dλ
为损伤对导热系数的影响因子；α
dc
为损伤对比热容的影响因子；α0是初始biot系数。
技术总结
本发明公开了一种考虑水冷冲击作用的干热岩水力裂缝扩展计算方法，具体为：收集干热岩地热储层的物理参数、拟注入的压裂液参数和工程参数；根据工程参数，绘制拟模拟压裂的岩体形状，得到计算域；对计算域添加热-流-固多物理场，建立相互的耦合关系，对多物理场耦合模型进行全耦合计算以及损伤变量计算；并计算损伤对各参数的影响，重复上述步骤，直到该注入压力下损伤变量不再增加；之后再增加注入压力，进行循环计算，直至岩体完全破裂即可停止计算。该方法充分考虑了低温压裂液注入高温地层引起的水冷冲击影响，可准确预测起裂压力及储层水力裂缝扩展形态，定量地评估水力压裂施工效果，为增强型地热系统的压裂设计提供理论基础和参考。基础和参考。基础和参考。
