一种模拟暂堵剂在干热岩粗糙裂隙内输运过程的数值方法

1.本发明涉及数值模拟技术领域，具体涉及一种模拟暂堵剂在干热岩粗糙裂隙内输运过程的数值方法。


背景技术：

2.随着全球经济的迅速发展，人们对能源的需求量越来越大。但是，人们长期所依赖的常规能源，如煤、石油、天然气等都是一次性不可再生能源，而且在利用这些常规能源时不可避免地对人类生存环境产生巨大的污染。因此，学者们都在积极寻找其他清洁的、可再生的新型替代能源。
3.干热岩是一种清洁、可再生，其内部没有或极少量含有水或蒸汽的热岩体能源，主要由各种变质岩或结晶类岩体组成，较常见的岩石有黑云母片麻岩、花岗岩、花岗闪长岩等。其埋深较浅、温度较高、有开发经济价值，保守估计地壳中干热岩(3
‑
10km深处)所蕴含的能量相当于全球所有石油、天然气和煤炭所蕴藏能量的30倍。
4.水力压裂技术是提取干热岩资源地热能的一项重要技术，常规水力压裂形成的缝网结构较为单一，因此干热岩采热效率往往不够理想。暂堵转向压裂是提高干热岩缝网复杂程度的一项新技术手段，通过向人工裂隙内加入暂堵剂，提升缝内净压力，迫使人工裂隙发生转向，以增加人工缝网复杂程度。因此，需要对暂堵剂在裂隙内的输运过程有着精准的认识。
5.目前，针对暂堵剂在干热岩裂隙内输运过程的研究方法主要有实验研究和数值模拟研究两种。
6.现阶段实验研究主要通过影响暂堵剂在裂隙内输运过程的相关参数，如暂堵剂浓度、注入速度等，对暂堵剂对人工裂隙封堵情况进行研究，然而实验过程中无法观察到暂堵剂缝内运移过程。为了克服这一问题，相关学者开发了配有高速摄像机的暂堵实验装置，但该装置无法模拟地层中的高压阻力，难以复现真实工况下暂堵剂在干热岩裂隙内的流动状态。因此，实验研究难以对暂堵剂在裂隙内的运移过程进行精准刻画。
7.数值模拟研究通过将暂堵剂在裂隙内的输运过程抽象为流体
‑
颗粒两相流在光滑平行板裂隙模型内输运过程的问题进行分析，通常采用的数值模拟方法主要为欧拉
‑
欧拉描述方法。欧拉
‑
欧拉描述方法的数学方程中暂堵剂(固体)和流体的运动方程均采用欧拉描述进行计算。该方法中，流体相被处理为连续介质，固体相则被处理为拟流体，使其动力学特性与液相类似，流体相和拟流体相控制方程组表达形式一致，这种方法具有较快的计算速度。但是该方法对暂堵剂的连续性假设削弱了固
‑
液两相流动的真实性，暂堵剂的材料属性及暂堵剂间相互作用等信息无法得到体现，模拟结果无法展示暂堵剂的复杂运动过程，从而无法准确、真实地揭示暂堵剂的运移机理。因此，现阶段对暂堵剂在干热岩缝内输运过程的数值模拟研究主要存在三个问题：首先，现有的暂堵剂运移数学模型未考虑高温作用下暂堵剂力学封堵性能的影响，无法准确刻画干热岩流动换热过程中裂隙内温度变化对封堵效率的影响；其次，现阶段数值模拟研究所用的物理模型通常为光滑平行板模型，未
考虑裂隙粗糙表面对暂堵剂输运过程的影响；最后，基于欧拉
‑
欧拉描述方法的数值模拟研究难以准确揭示暂堵剂缝内运移机理。


技术实现要素：

8.有鉴于此，本发明的目的在于克服现有技术的不足，提供一种模拟暂堵剂在干热岩粗糙裂隙内输运过程的数值方法，以解决现有技术中数值模拟无法真实反映暂堵剂运动状态的问题。
9.为实现以上目的，本发明采用如下技术方案：
10.根据本发明实施例的第一方面，提供一种模拟暂堵剂在干热岩粗糙裂隙内输运过程的数值方法，包括：
11.步骤s1、构建裂隙几何模型，确定所述裂隙几何模型的粗糙度；
12.步骤s2、确定所需数学模型，包括：连续相数学模型、离散相数学模型和交叉耦合方程；将计算区域内的流体介质视为连续相，通过连续相数学模型进行求解；将暂堵剂视为离散相，通过离散相数学模型进行求解；连续相与连续相之间的相互作用，通过simple算法进行求解；连续相与离散相之间相互作用，通过交叉耦合方程进行求解；
13.步骤s3、对所述数学模型进行离散，给出相应边界条件和初始化数值，通过预设的双向耦合计算流程，模拟暂堵剂在干热岩粗糙裂隙内输运过程；
14.步骤s4、输出模拟结果。
15.优选地，所述步骤s1中确定所述裂隙几何模型的粗糙度，通过以下项中的任一种方法确定，包括：
16.节理粗糙度直接对比方法、节理粗糙度经验公式计算方法、分形维数方法和人工设计粗糙形状。
17.优选地，所述步骤s2中，连续相数学模型通过以下方法确定，包括：
18.流体运动采用质量守恒方程、动量守恒方程和能量守恒方程进行求解，其中，质量守恒方程表示为：
[0019][0020]
其中，s
m
表示源项，为离散相添加到连续相中的质量，或其他任何需要考虑到连续相中的相；ρ表示流体密度；t表示时间；表示拉普拉斯算子；表示流体速度；表示偏导数；
[0021]
动量守恒方程表示为：
[0022][0023]
其中，p表示缝内流体压力；表示流体重力加速度；表示外部力，例如，与离散相相互作用产生的力；表示切应力张量，其表达式如下：
[0024]
[0025]
其中，μ表示流体粘度；i表示单位张量；上标t表示流体速度向量的转置；
[0026]
能量守恒方程表示为：
[0027][0028]
k
eff
＝k+k
t
ꢀꢀꢀ
(5)
[0029][0030][0031][0032]
其中，t为流体温度；p为出口压力；e为能量常数；k
eff
为有效导热率；k为流体材料导热率；k
t
为流体湍流导热率；为流体扩散通量；为有效切应力张量；s
h
为体积热源；h为不可压缩流体的焓；y
j
为流体质量分数；h
j
为流体质量分数为j时的焓；t
ref
为298.15k；c
f,j
为流体定压比热容；∫表示积分。
[0033]
优选地，所述步骤s2中，离散相数学模型通过以下方法确定，包括：
[0034]
通过牛顿第二定律对离散相暂堵剂进行力平衡积分，计算暂堵剂的运动轨迹，其表达式为：
[0035][0036][0037][0038]
其中，m
p
为暂堵剂质量；为暂堵剂速度；ρ
p
为暂堵剂密度；d
p
为暂堵剂直径；为流体速度；μ为流体粘度；ρ为流体密度；τ
r
为暂堵剂弛豫时间；re为雷诺数；c
d
为常数，取0.5；d表示导数；为附加力，可根据实际情况选择；
[0039]
暂堵剂旋转通过求解暂堵剂角动量的附加常微分方程：
[0040]
[0041][0042][0043]
其中，为力矩；为流体与暂堵剂之间的相对角速度；为流体与暂堵剂之间相对角速度的绝对值；i
p
为惯性矩；为暂堵剂角速度；c
ω
为旋转阻力系数。
[0044]
优选地，所述步骤s2中，确定交叉耦合方程，包括：
[0045]
考虑连续相对离散相的影响，构建质量、动量和能量交换方程；
[0046]
考虑离散相对连续相的影响，构建修正后的连续性方程、动量方程和能量方程；
[0047]
考虑连续相内部的相互影响，通过simple算法获得速度、压力和温度耦合的连续相场。
[0048]
优选地，所述考虑连续相对离散相的影响，构建质量、动量和能量交换方程，包括：
[0049]
流体和暂堵剂之间的质量交换公式：
[0050][0051]
其中，m表示暂堵剂质量变化值；δm
p
表示控制体内暂堵剂的质量变化；表示跟踪暂堵剂的初始质量流率；m
p,0
表示控制体内暂堵剂初始质量。
[0052]
流体和暂堵剂之间的动量交换公式：
[0053][0054]
其中，f表示暂堵剂动量变化值；表示暂堵剂质量流率；δt表示时间步长；f
other
表示其它相间作用力；
[0055]
流体和暂堵剂之间的能量交换公式：
[0056][0057]
其中，q表示暂堵剂热量变化值；表示控制体内的暂堵剂平均质量；c
p
表示暂堵剂比热容；δt
p
表示控制体内暂堵剂温度变化；h
fg
表示暂堵剂挥发或者析出时的热量；h
pyrol
表示暂堵剂分解时所需热量；t
p
表示离开控制体时暂堵剂的温度；t
ref
表示暂堵剂在发生相变时对应的参考温度；c
p,i
表示暂堵剂挥发或者析出时的比热。
[0058]
优选地，所述考虑离散相对连续相的影响，构建修正后的连续性方程、动量方程和能量方程，包括：
[0059]
通过将连续性方程、动量方程和能量方程进行改写，添加离散相相关信息，完成离散相对连续相影响过程，具体为：
[0060]
连续性方程改写为：
[0061][0062]
其中，α
p
表示离散相体积分数；m
lp
表示连续相到离散相的质量传递；m
pl
表示离散相到连续相的质量传递；
[0063]
动量方程改写为：
[0064][0065]
其中，k
lp
表示动量交换系数；表示连续相相对于离散相的速度；表示离散相相对于连续相的速度,当m
lp
＞0时，当m
lp
＜0时，当m
pl
＞0时，当m
pl
＜0时，0时，表示离散相切应力张量，通过如下公式计算：
[0066][0067]
其中，i表示单位张量；上标t表示离散相速度向量的转置；
[0068]
能量方程改写为：
[0069][0070]
其中，h
p
表示离散相的焓；p
p
表示离散相受到的压力；表示离散相的热通量；s
p
表示离散相的热源项；q
lp
表示连续相与离散相之间的热交换；h
lp
表示连续相的相变焓；h
pl
表示离散相的相变焓；为离散相扩散通量；h
j,p
为离散相质量分数为j时的焓。
[0071]
优选地，所述考虑连续相内部的相互影响，通过simple算法获得速度、压力和温度耦合的连续相场，包括：
[0072]
通过simple算法，计算连续性方程和动量方程获得压力场和速度场；
[0073]
利用速度场的相关参数去计算能量方程，得到能量场，从而得到完整的连续相场；所述连续相场包括速度场、压力场和能量场。
[0074]
优选地，所述预设的双向耦合计算流程，包括：
[0075]
步骤s31、设定输入模块，给出连续相和离散相相关参数，进行预设步数的连续相计算，得到一个连续相场；
[0076]
步骤s32、将当前连续相场相关参数应用到对离散相的计算中，计算离散相质量、动量和能量交换公式，得到连续相场对离散相的影响参数；
[0077]
步骤s33、将当前连续相对离散相的相关影响参数应用到对离散相的计算中，对牛顿第二运动定律相关方程进行求解，得到离散相运动轨迹、运动速度等相关参数；
[0078]
步骤s34、将当前离散相相关参数应用到对连续相的计算中，通过计算修正后的连续性方程、动量方程和能量方程，得到修正后的连续相场，从而完成一次双向耦合计算；
[0079]
步骤s35、判定是否达到步骤s31中设置的计算步数，如不满足条件，返回步骤s32继续进行计算；如果满足条件，结束计算。
[0080]
优选地，所述步骤s31，包括：
[0081]
给定初始假设的速度和压力条件，利用simple算法对流体运动质量守恒方程、动量守恒方程进行求解，得到连续相场的真实压力和速度相关信息，通过能量方程中速度梯度项，将连续相场真实速度相关信息引入能量方程，通过对能量方程进行求解，得到连续相场的真实能量相关信息，进一步得到完整的连续相场，完成连续相部分求解；和/或，
[0082]
所述步骤s32，包括：
[0083]
当暂堵剂穿过控制体时，流体和暂堵剂之间的质量、动量和能量交换过程通过计算暂堵剂的质量交换公式、动量交换公式和热量交换公式来求解连续相传递给离散相的质量、动量和能量值；和/或，
[0084]
所述步骤s33，包括：
[0085]
将所述质量、动量和能量值代入离散相数学模型，通过牛顿第二定律对离散相暂堵剂进行力平衡积分，计算离散相暂堵剂的运动轨迹、运动速度和受力情况等相关信息，描绘出暂堵剂在当前连续相场内的运动状态；和/或，
[0086]
所述步骤s34，包括：
[0087]
对连续性方程、动量方程和能量方程进行修正，得到修正后的连续性方程、动量方程和能量方程，将当前计算所得质量变化值、动量变化值、能量变化值和离散相运动轨迹、运动速度等信息引入，重新计算连续相场，得到修正后的连续相场，完成一次双向耦合；和/或，
[0088]
所述步骤s35，包括：
[0089]
将当前计算步数与步骤s31中设定的计算步数进行比较，如果当前步数小于设定步数，则返回步骤s32，继续进行计算；若果当前步数等于设定步数，则判定计算完成，结束计算。
[0090]
本发明采用以上技术方案，至少具备以下有益效果：
[0091]
本发明提供的技术方案，采用欧拉
‑
拉格朗日描述方法对暂堵剂在干热岩粗糙裂隙内的运移过程进行研究，构建了干热岩粗糙裂隙物理模型，将流体视作连续相，将暂堵剂看作离散相，既考虑流体流动对暂堵剂运动的影响，又考虑暂堵剂运动对流体的影响，还考虑温度变化(通过能量方程实现)和裂隙粗糙表面对暂堵剂力学封堵性能的影响，从而实现连续相和离散相的双向耦合计算，达到同时考虑连续相和离散相相互作用，准确求解流
‑
固双向耦合问题的目的。
[0092]
进一步地，单个暂堵剂的运动遵循牛顿第二运动定律，在拉格朗日坐标系下求解其运动方程，得到不同时刻暂堵剂的位置、速度和受力情况等信息，从而克服了欧拉
‑
欧拉描述方法的缺陷，真实地反映了暂堵剂在流场中的运动过程。
附图说明
[0093]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现
有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0094]
图1为本发明一实施例提供的一种模拟暂堵剂在干热岩粗糙裂隙内输运过程的数值方法的流程图；
[0095]
图2为本发明一实施例提供的节理粗糙度(jrc)直接对比方法确定的标准jrc粗糙程度示意图；
[0096]
图3为本发明一实施例提供的利用人工设计粗糙形状构建的裂隙几何模型的示意图；
[0097]
图4为本发明一实施例提供的双向耦合计算流程的示意图；
[0098]
图5为本发明一实施例提供的模拟暂堵剂在干热岩粗糙裂隙内输运过程的效果图。
具体实施方式
[0099]
为使本发明的目的、技术方案和优点更加清楚，下面将对本发明的技术方案进行详细的描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施方式，都属于本发明所保护的范围。
[0100]
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
[0101]
实施例一
[0102]
参见图1，根据本发明一实施例提出的一种模拟暂堵剂在干热岩粗糙裂隙内输运过程的数值方法，包括：
[0103]
步骤s1、构建裂隙几何模型，确定所述裂隙几何模型的粗糙度；
[0104]
步骤s2、确定所需数学模型，包括：连续相数学模型、离散相数学模型和交叉耦合方程；将计算区域内的流体介质视为连续相，通过连续相数学模型进行求解；将暂堵剂视为离散相，通过离散相数学模型进行求解；连续相与连续相之间的相互作用，通过simple算法进行求解；连续相与离散相之间相互作用，通过交叉耦合方程进行求解；
[0105]
步骤s3、对所述数学模型进行离散，给出相应边界条件和初始化数值，通过预设的双向耦合计算流程，模拟暂堵剂在干热岩粗糙裂隙内输运过程；
[0106]
步骤s4、输出模拟结果。
[0107]
优选地，所述步骤s1中确定所述裂隙几何模型的粗糙度，通过以下项中的任一种方法确定，包括：
[0108]
节理粗糙度(jrc)直接对比方法、节理粗糙度(jrc)经验公式计算方法、分形维数(d
f
)方法和人工设计粗糙形状四种方法。
[0109]
第一种方法：节理粗糙度(jrc)直接对比方法：
[0110]
将岩石裂隙表面节理粗糙度划分为10个标准，jrc的范围从0
‑
2，2
‑
4等，一直到16
‑
18和18
‑
20对应10个标准jrc粗糙程度示意图，如图2所示。图2为10个标准jrc粗糙程度示意图，通过将裂隙粗糙表面轮廓直接与标准等级进行对比的方法，确定干热岩粗糙裂隙表面节理粗糙度，jrc数值越大，则相应裂隙表面粗糙度越大。
[0111]
第二种方法：节理粗糙度(jrc)经验公式计算方法：
[0112]
首先，通过ct扫描技术可以获得干热岩人工裂隙粗糙表面的详细数据信息；接着，利用公式(22)计算裂隙剖面无量纲高度z2，公式如下：
[0113][0114]
式中z2表示裂隙剖面无量纲高度；n表示计算用数据点个数；δx表示两个数据点之间的距离；z
i
表示某一数据点高度；z
i+1
表示与z
i
数据点相邻数据点高度；∑表示求和。
[0115]
将裂隙剖面无量纲高度z2带入公式(23)可以计算得到干热岩粗糙裂隙表面节理粗糙度。
[0116]
jrc＝32.2+32.47logz2ꢀꢀꢀ
(23)
[0117]
式中jrc表示裂隙表面节理粗糙度；log表示以10底的对数。
[0118]
第三种方法：分形维数(d
f
)方法：
[0119]
通过ct扫描技术可以获得干热岩裂隙粗糙表面的详细数据信息，利用公式(24)计算裂隙轮廓相对高度波动σ
s
，公式如下：
[0120][0121]
式中σ
s
表示裂隙轮廓相对高度波动；m表示所有数据点的数量；s表示两个数据点之间间隔数据点的数量；y
i
表示数据点的高度；y
i+s
表示与y
i
相邻数据点的高度。
[0122]
裂隙轮廓相对高度波动σ
s
与亨斯特指数h之间的关系如公式(25)所示。
[0123]
σ
s
∝
s
h
ꢀꢀꢀ
(25)
[0124]
将公式中σ
s
与s对应数据点绘制为双对数曲线图形式，对图中线性较好部分进行拟合，得到一条拟合直线，读取直线上的数据点，求出该直线的斜率，该斜率作为亨斯特指数h。将其带入公式(26)计算可得由分形维数表示的干热岩裂隙表面粗糙度。
[0125]
d
f
＝2
‑
h
ꢀꢀꢀ
(26)
[0126]
式中d
f
表示干热岩裂隙粗糙表面分形维数；h表示亨斯特指数。
[0127]
第四种方法：人工设计粗糙形状生成，如矩形、三角形或者正弦曲线形状，其粗糙度由以下公式计算：
[0128][0129]
式中c表示干热岩裂隙表面粗糙度；h
c
表示粗糙峰值高度；h表示平滑裂隙宽度。
[0130]
优选地，所述步骤s2中，确定所需数学模型包括三个部分，分别为连续相数学模型、离散相数学模型和交叉耦合方程。计算区域内的流体介质视为连续相，通过连续相数学模型进行求解；暂堵剂视为离散相，通过离散相数学模型进行求解；连续相与连续相之间的相互作用，通过simple算法进行求解；连续相与离散相之间相互作用，如质量、动量和能量交换通过耦合方程进行求解。
[0131]
所述数学模型的边界条件和初始条件为：入口处采用速度边界条件，出口处采用压力边界条件，裂隙壁面采用第一类温度边界条件，即壁面温度恒定不变。
[0132]
值得说明的是，温度边界条件可以为第一类边界条件、第二类边界条件和第三类边界条件中的任意一种。
[0133]
值得说明的是，初始条件数值可根据需要设置。
[0134]
优选地，所述步骤s3中预设的双向耦合计算流程，包括：
[0135]
步骤s31、设定输入模块，给出连续相和离散相相关参数，进行预设步数的连续相计算，得到一个连续相场；
[0136]
步骤s32、将当前连续相场相关参数应用到对离散相的计算中，计算离散相质量、动量和能量交换公式，得到连续相场对离散相的影响参数；
[0137]
步骤s33、将当前连续相对离散相的相关影响参数应用到对离散相的计算中，对牛顿第二运动定律相关方程进行求解，得到离散相运动轨迹、运动速度等相关参数；
[0138]
步骤s34、将当前离散相相关参数应用到对连续相的计算中，通过计算修正后的连续性方程、动量方程和能量方程，得到修正后的连续相场，从而完成一次双向耦合计算；
[0139]
步骤s35、判定是否达到步骤s31中设置的计算步数，如不满足条件，返回步骤s32继续进行计算；如果满足条件，结束计算。
[0140]
优选地，所述步骤s4、输出模拟结果，具体为：
[0141]
为帮助分析干热岩粗糙裂隙内暂堵剂输运过程，可以在计算结束后查看所需计算结果，包括裂隙区域内的流体速度云图和裂隙区域内压力云图，当前时刻、当前区域内每个暂堵剂的速度和受力情况，以及暂堵剂输运过程动态图等。
[0142]
可以理解的是，本实施例提供的技术方案，精准刻画了暂堵剂在干热岩粗糙裂隙内的瞬态运移过程，通过交替求解连续相数学模型、离散相数学模型和交叉耦合方程，达到同时考虑连续相和离散相相互作用的目的，进而准确求解了流
‑
固双向耦合的问题。
[0143]
实施例二
[0144]
根据本发明另一实施例提出的一种模拟暂堵剂在干热岩粗糙裂隙内输运过程的数值方法，包括如下五个部分：输入模块初始化、裂隙几何模型构建、数学模型构建、双向耦合计算流程和输出计算结果。
[0145]
1、输入模块初始化：
[0146]
开始计算前，需要进行初始化设置，首先给定流体的初始速度为0.5m/s、流体粘度为3mpa
·
s和流体密度为1000kg/m3；其次给定暂堵剂注入速度为0.5m/s、暂堵剂直径为0.003m、暂堵剂密度为2500kg/m3等信息；接着需要设置边界参数，包括入口速度为0.5m/s、湍流强度为5％和水力直径0.012m；出口边界条件包括出口压力为30mpa、湍流强度为5％和水力直径为0.012m；壁面粗糙，采用第一类边界条件，给定温度500k，当暂堵剂与壁面发生接触时，暂堵剂被反射；最后，给定流体和暂堵剂的时间步长都为δt＝0.001s，计算步数为10000步。流体、暂堵剂和边界条件的有关输入参数见表1
‑
表5。
[0147]
表1流体参数
[0148][0149]
表2暂堵剂参数
[0150][0151]
表3入口边界条件
[0152][0153][0154]
表4出口边界条件
[0155][0156]
表5壁面边界条件
[0157][0158]
2、裂隙几何模型构建
[0159]
通过第四种方法，即通过人工设计粗糙形状生成，设计干热岩裂隙粗糙表面为矩形方波形状，通过粗糙度计算公式(27)计算得到裂隙表面粗糙度，具体如下：
[0160][0161]
式中，h
c
表示粗糙峰值高度，给定初始值为0.9mm；h表示平滑裂隙宽度，给定初始值为6mm；c表示干热岩裂隙表面粗糙度，计算可得干热岩裂隙表面粗糙度为0.15。
[0162]
参见图3，图3表示通过人工设计粗糙形状方法生成的裂隙表面粗糙度为0.15的物理模型，计算区域长400mm，高200mm，宽6mm，裂隙粗糙位置由矩形方波表示，分布在裂隙左右表面。
[0163]
3、数学模型构建
[0164]
计算所需数学模型包括三个部分，分别为连续相数学模型、离散相数学模型和交叉耦合方程。
[0165]
(1)连续相数学模型
[0166]
计算区域内的流体介质视为连续相，通过连续相数学模型公式(1)～(8)进行求解，具体如下：
[0167]
流体运动采用质量守恒方程、动量守恒方程和能量守恒方程进行求解，质量守恒方程表达式为:
[0168][0169]
式中s
m
表示源项，为暂堵剂添加到流体中的质量；ρ表示流体密度；t表示时间；表示拉普拉斯算子；表示流体速度；表示偏导数。
[0170]
动量守恒方程表示为：
[0171][0172]
式中p表示缝内流体压力；表示流体重力加速度；表示暂堵剂间相互作用产生的力；表示切应力张量，其表达式如下：
[0173][0174]
式中μ表示流体粘度；i表示单位张量；上标t表示流体速度向量的转置。
[0175]
能量守恒方程表示为：
[0176][0177]
k
eff
＝k+k
t
ꢀꢀꢀ
(5)
[0178][0179][0180][0181]
式中t为流体温度；p为出口压力；e为能量常数；k
eff
为有效导热率；k为流体材料导热率；k
t
为流体湍流导热率；为流体扩散通量；为有效切应力张量；s
h
为体积热源；h为不可压缩流体的焓；y
j
为流体质量分数；h
j
为流体质量分数为j时的焓；t
ref
为298.15k；c
f,j
为流体定压比热容；∫表示积分。
[0182]
(2)离散相数学模型
[0183]
暂堵剂视为离散相，通过离散相数学模型公式(9)～(15)进行求解，具体如下：
[0184]
通过牛顿第二定律对暂堵剂进行力平衡积分，计算暂堵剂的运动轨迹，其表达式为：
[0185][0186][0187]
[0188]
式中m
p
为暂堵剂质量；为暂堵剂速度；ρ
p
为暂堵剂密度；d
p
为暂堵剂直径；为流体速度；μ为流体粘度；ρ为流体密度；τ
r
为暂堵剂弛豫时间；re为雷诺数；c
d
为常数，取0.5；d表示导数；为暂堵剂与壁面间相互碰撞产生的力。
[0189]
暂堵剂旋转通过求解暂堵剂角动量的附加常微分方程：
[0190][0191][0192][0193]
式中为力矩；为流体与暂堵剂之间的相对角速度；为流体与暂堵剂之间相对角速度的绝对值；i
p
为惯性矩；为暂堵剂角速度；c
ω
为旋转阻力系数。
[0194]
(3)交叉耦合方程
[0195]
流体对暂堵剂的影响通过质量、动量和能量交换公式(15)～(17)进行求解，具体如下：
[0196]
质量交换公式表示为：
[0197][0198]
式中，m表示暂堵剂质量变化值；δm
p
表示控制体内暂堵剂的质量变化；表示跟踪暂堵剂的初始质量流率；m
p,0
表示控制体内暂堵剂初始质量。
[0199]
动量交换公式表示为：
[0200][0201]
式中，f表示暂堵剂动量变化值；表示暂堵剂质量流率；δt表示时间步长；f
other
表示暂堵剂间相互作用力。
[0202]
能量交换公式表示为：
[0203][0204]
式中，q表示暂堵剂热量变化值；表示控制体内的暂堵剂平均质量；c
p
表示暂堵剂比热容；δt
p
表示控制体内暂堵剂温度变化；h
fg
表示暂堵剂析出时的热量；h
pyrol
表示暂堵剂分解时所需热量；t
p
表示离开控制体时暂堵剂的温度；t
ref
表示暂堵剂在发生相变时对应的参考温度；c
p,i
表示暂堵剂析出时的比热。
[0205]
考虑暂堵剂对流体的影响，通过将连续性方程、动量方程和能量方程进行改写，添加暂堵剂相关信息，完成暂堵剂对流体影响过程，暂堵剂对流体的影响通过修正后的连续性方程、动量方程和能量方程(18)～(21)进行求解，具体如下：
[0206]
连续性方程改写为：
[0207][0208]
式中，α
p
表示暂堵剂体积分数；m
lp
表示流体到暂堵剂的质量传递；m
pl
表示暂堵剂到流体的质量传递。
[0209]
动量方程改写为：
[0210][0211]
式中，k
lp
表示动量交换系数；表示流体相对于暂堵剂的速度；表示暂堵剂相对于流体的速度；当m
lp
＞0时，当m
lp
＜0时，当m
pl
＞0时，当m
pl
＜0时，＜0时，表示暂堵剂切应力张量，通过如下公式计算：
[0212][0213]
能量方程改写为：
[0214][0215]
式中，h
p
表示暂堵剂的焓；p
p
表示暂堵剂受到的压力；表示暂堵剂的热通量；s
p
表示暂堵剂的热源项；q
lp
表示流体与暂堵剂之间的热交换；h
lp
表示流体的相变焓；h
pl
表示暂堵剂的相变焓；为暂堵剂扩散通量；h
j,p
为暂堵剂质量分数为j时的焓。
[0216]
4、双向耦合计算流程
[0217]
参见图4，通过交替求解连续相数学模型、离散相数学模型和交叉耦合方程完成双向耦合过程，具体过程包括：
[0218]
给定初始假设的速度和压力条件，利用simple算法(即速度
‑
压力耦合算法)对流体运动质量守恒方程、动量守恒方程进行求解，得到流体场的真实压力和速度相关信息，并通过能量方程中的速度梯度项将速度场相关信息代入能量方程，对能量方程进行离散，速度场相关信息被处理为离散后的能量方程中的系数项，完成动量方程与能量方程的耦合，求解离散后的能量方程，得到流体场的真实能量相关信息，进一步得到完成的流体场，完成流体部分求解。
[0219]
当暂堵剂穿过控制体时，流体和暂堵剂之间的质量、动量和能量交换过程通过计
算暂堵剂的质量交换公式、动量交换公式和热量交换公式求解流体传递给暂堵剂的质量、动量和能量值，这些物理量可作用于随后的暂堵剂计算中去。
[0220]
通过牛顿第二定律对暂堵剂进行力平衡积分，计算暂堵剂的运动轨迹、运动速度和受力情况等相关信息，准确刻画出暂堵剂在当前流体场内的运动状态，至此，完成考虑到流体相关参数对暂堵剂的作用。
[0221]
对连续性方程、动量方程和能量方程进行修正，得到修正后的连续性方程、动量方程和能量方程，将当前计算所得质量变化值、动量变化值、能量变化值和暂堵剂运动轨迹、运动速度等信息引入，重新计算流体场，得到修正后的流体场，完成双向耦合计算。
[0222]
判定是否达到设定计算步数，若没有达到设定计算步数，则继续计算；若达到设定计算步数，则结束计算。
[0223]
5、输出计算结果
[0224]
为帮助分析干热岩粗糙裂隙内暂堵剂输运过程，可以在计算结束后查看所需计算结果，包括裂隙区域内的流体速度云图和裂隙区域内压力云图，当前时刻、当前区域内每个暂堵剂的速度和受力情况，以及暂堵剂输运过程动态图等。
[0225]
流体初始输入参数如表1所示，暂堵剂相关参数如表2所示，边界条件相关参数如表3
‑
表5所示。使用相同的网格尺寸与模拟参数，按上述实施例提供的干热岩粗糙裂隙内暂堵剂流动模拟方法，以暂堵剂直径为自变量，模拟了暂堵剂直径分别为1mm、2mm、3mm、4mm和5mm情况下暂堵剂受力变化趋势，将模拟结果进行比较，如图5所示。
[0226]
由图5分析可知，对于干热岩粗糙裂隙模型而言，随着暂堵剂直径的增大，其所受平均作用力逐渐增大，但受力增量不同，整体相差近40倍；粗糙裂缝内，暂堵剂间相互作用力的影响处于次要地位，颗粒受力主要由颗粒在不同速度下与粗糙壁面的矩形部分碰撞产生。
[0227]
可以理解的是，本实施例提供的技术方案，采用欧拉
‑
拉格朗日描述方法对暂堵剂在干热岩粗糙裂隙内的运移过程进行研究，构建了干热岩粗糙裂隙物理模型，将流体视作连续相，将暂堵剂看作离散相，既考虑流体流动对暂堵剂运动的影响，又考虑暂堵剂运动对流体的影响，还考虑温度变化(通过能量方程实现)和裂隙粗糙表面对暂堵剂力学封堵性能的影响，从而实现连续相和离散相的双向耦合计算，达到同时考虑连续相和离散相相互作用，准确求解流
‑
固双向耦合问题的目的。
[0228]
进一步地，单个暂堵剂的运动遵循牛顿第二运动定律，在拉格朗日坐标系下求解其运动方程，得到不同时刻暂堵剂的位置、速度和受力情况等信息，从而克服了欧拉
‑
欧拉描述方法的缺陷，真实地反映了暂堵剂在流场中的运动过程。
[0229]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。术语“多个”指两个或两个以上，除非另有明确的限定。