增强型地热系统热流固耦合模型及数值模拟

孙致学, 徐 轶, 吕抒桓, 徐 杨, 孙 强, 蔡明玉, 姚 军

(1.中国石油大学石油工程学院,山东青岛 266580; 2.武汉大学水资源与水电工程科学国家重点实验室,湖北武汉 430072)

增强型地热系统热流固耦合模型及数值模拟

孙致学1, 徐 轶2, 吕抒桓1, 徐 杨1, 孙 强1, 蔡明玉1, 姚 军1

(1.中国石油大学石油工程学院,山东青岛 266580; 2.武汉大学水资源与水电工程科学国家重点实验室,湖北武汉 430072)

增强型地热系统(EGS)利用水力压裂地下高温岩体形成人工热储，通过载热流体循环提取干热岩(HDR)所存储的地热能，其开采过程包含渗流、热能交换及岩体介质变形，为典型的热流固(THM)耦合问题。将裂隙岩体视作基于离散裂隙网络和基质岩体的双重介质，给出THM耦合的数学模型，基于商业有限元软件COMSOL Mutilphysics进行二次开发，实现裂隙岩体温度场、渗流场和应力场的全耦合求解。利用一个已知解析解的算例验证耦合模型和全耦合求解方法的正确性。最后利用随机生成的二维裂隙网络模型模拟EGS的运行过程，分析干热岩储层内渗流、温度、应力及变形的分布规律。计算结果表明，储层内的贯通裂隙构成主要导水区，水的对流传热作用明显；高压水注入和温度变化导致岩体裂隙发生位移，改变储层的传输特性，进而影响地热能的提取；考虑THM耦合作用对于研究增强型地热系统的的利用效率和运行规律非常有必要。

干热岩； 增强型地热系统； 裂隙岩体； 热流固耦合； 数值模拟

干热岩(hot dry rock,HDR)是指埋藏于地面以下(深度3～10 km,温度150～650 ℃)内部不存在流体或仅有少量地下流体的高温岩体[1-2]。干热岩地热资源由于其清洁可再生性、分布的广泛性和高热储温性而更具开发价值与前景,被广泛认为是21世纪最有潜力的新型能源[3]。由于埋藏较深,干热岩岩体致密,渗透性极低,难以用传统方法开采。增强型地热系统(enhanced geothermal system,EGS)是采用水力压裂人工形成地热储层,通过水在注入井、热储层和生产井之间进行循环,从干热岩中经济地提取热量的工程[4]。EGS是目前干热岩能源开发工程的关键技术,至今已有40多年的研究历史[5]。EGS的开发与运行过程中涉及到裂隙岩体内水动力场、应力场以及温度场的复杂时空演化规律,即热流固(thermo-hydro-mechanical,THM)耦合问题。数值模拟被认为是研究该问题的有效途径之一。通过建立数学模型,模拟裂隙岩体内的传热、渗流与应力变化过程,对于预测EGS系统的运行状态、利用效率、使用寿命以及解决其开发过程中遇到的实际问题都具有重要意义。近年来,裂隙岩体THM耦合数学模型的建立成为研究增强型地热系统数值模拟的热点问题之一。建立的数学模型主要有离散裂隙网络模型和等效连续介质模型两类。离散裂隙网络模型因为更接近于实际储层、能够模拟裂隙中的水热迁移过程,引起了广泛的研究兴趣[6-15]。尽管国内外学者都针对裂隙岩体热流固耦合数值模拟开展了大量有益的探索,但是从整体上看还很不成熟,其中一个关键技术难题是如何模拟大规模复杂裂隙,以及有效处理裂隙介质中传热、渗流与应力过程的多场耦合作用[1,16]。尤其是应力场变化引起裂隙的张开和闭合会改变介质的渗透率,影响传热效率,这些耦合作用对EGS的运行状态将产生较大影响[13],而在以往的研究中被忽略。对于含复杂裂隙网络的高温地热岩体,THM耦合数学模型的建立和求解都因各变量互相作用、影响的复杂性而难以实现。在EGS数值模拟研究中,能较真实地反映其THM耦合实际形态的数学模型和求解方法是亟待解决的问题。笔者在赵阳升等[14-15,17]研究的基础上,将岩体视作由离散裂隙网络和基质岩块组成的双重介质,考虑传热、渗流与应力过程中的相互作用,建立复杂裂隙岩体的THM耦合数学模型，基于商业有限元软件COMSOL Mutilphysics进行二次开发,实现数学模型的全耦合求解;针对饱和土热弹性固结问题,将耦合模型模拟结果与解析解进行对比,以验证模型的可靠性;最后,将该模型应用于简化的增强型地热系统,研究其温度场、渗流场以及应力场的变化规律。

1 热流固耦合模型

1.1 基本假定

EGS裂隙储层可以简化为基质岩块-裂隙双重介质模型。将基质岩块视作孔隙介质,渗透率较低,水力压裂形成的人工裂缝网络构成了水的主要流动通道。开采过程中注入的水同时在基质岩块和裂隙网络中流动。在地下深部岩体中的水由于受高压作用,不可能汽化[17],因此可假定储层为单相水所饱和,水的流动为层流,服从达西定律。裂缝水以对流和传导的方式发生热量交换和热量传递。为了简化计算,认为岩体及裂隙始终处于弹性状态,并且基于小变形假设。

1.2 数学模型

(1)基质岩块应力和位移场方程，

σij,j+Fi=0,

(1)

μui,jj+(λ+μ)uj,ji-αBp,i-βTTs,i+Fi=0.

(2)

其中

μ=E/2(1+ν),λ=Eν/[(1+ν)(1-2ν)],

βT=αTE/(1-2ν).

式中,σij为应力二阶张量分量,Pa;u为位移,m;Fi为体力,Pa;μ与λ为拉梅常数;E为弹性模量,Pa;ν为泊松比;p为水压力,Pa;αB为Biot耦合系数,αB≤1;Ts为岩块温度,K;βT为热膨胀因子;αT为热膨胀系数,K-1;αBp,i和βTTs,i分别为水压力作用项和温度应力作用项。

(2)裂隙变形方程,

(3)

(4)

(3)基质岩块渗流场方程,

(5)

(6)

式中,t为时间,s;u为岩块中的水流速,m/s;κ为岩体渗透率,m2;η为流体的动力黏度,Pa·s;S为岩体的储水系数,Pa-1;e为岩块的体积应变;Q为渗流的源汇项,s-1。

(4)裂隙内渗流场方程,

(7)

其中

式中,Sf为裂隙的储水系数,Pa-1;κf为裂隙内的渗透率,m2;df为裂隙宽度,m;Qf为岩块与裂隙面的流量交换;n为裂隙面法向;ef为裂隙面的体积应变;τ为沿裂隙切向求导。

(5)基质岩块温度场方程,

(8)

式中,ρs为岩块密度,kg/m3;λs为岩块的热传导系数,W/(m·K);cs为岩块比热容,J/(kg·K);W为热源,W/m3。

由于孔隙率低,基质岩块中水的流速较小,可以认为其中的水温等于岩块温度,因此不考虑岩块孔隙中的对流作用。

(6)裂隙水温度场方程,

(9)

式中,ρf为水的密度,kg/m3;cf为水的比热容,J/(kg·K);λf为水的热传导系数,W/(m·K);uf为裂隙内水的流速,m/s;Tf为裂隙内水的温度,K;Wf为裂隙表面水从基岩吸收的热量,W/m2;第二项表示裂隙水对温度场的对流作用。

假定裂隙水与基岩的热交换服从牛顿换热公式,流体和固体在边界处的换热量用对流传热系数计算,即认为裂隙内水的平均温度与基岩表面温度不同,单位面积上从基岩传给裂隙水的热量[8-9]为

Wf=h(Ts-Tf).

(10)

式中,h为换热系数,W/(m2·K)。

1.3 耦合作用

1.3.1 裂隙面渗流与应力的耦合特性

岩体内的孔隙、裂隙介质由于应力作用发生剧烈变形,其渗透系数必然会发生不同程度的变化,并且与岩体的应力水平或应变水平相关;特别是裂隙发育的岩体中,裂隙作为主要的渗流通道,应力变化对其渗透性产生极大影响。研究应力与渗透系数的动态关系是研究裂隙岩体应力场、渗流场耦合特性的重要内容。Louis[18]利用岩体钻孔压水试验对单裂隙面渗流与应力的关系进行了探索,提出了指数型的经验公式,

(11)

式(11)在渗流场与应力场耦合分析中使用简单、方便[19],因此较为常用。本文中假定基质岩块的渗透系数为常数,采用式(11)描述裂隙应力水平与渗透系数的关系。

1.3.2 流体性质变化

流体性质也是影响耦合作用的一个重要因素。在高温高压作用下,水的密度ρf不再是常数,一般可以表示为温度和压力的函数[17],

1/ρf=3.086-0.899 017(4 014.15-T)0.147 166-

0.39(658.15-T)-1.6(p-225.5)+δ.

(12)

式中,T为水的绝对温度;δ为关于水温和压力的函数,一般不超过1/ρf的6%。

此外,水温对水的动力黏度影响较大。水的动力黏度η=υρf,υ为水的运动黏滞系数,m2/s。一般υ可按经验公式[20]计算:

(13)

水的密度和黏度都与温度相关。温度的变化直接影响渗流和传热两个过程,二者的作用是一种强耦合关系。

2 耦合模型的有限元求解

以上述各组方程辅以初始条件和边界条件,就构成了基于离散裂隙网络的裂隙岩体温度场—渗流场—应力场耦合数学模型,一般采用有限单元法进行模型的数值求解。由于耦合模型含有的因变量众多,求解复杂。以往的求解策略多是采用单向耦合或部分耦合(将不同介质的固体变形、热传导方程和渗流看作独立的子系统进行有限元分析,用一个物理场的结果作为另一个物理场的初始值,或在上述的基础上反复交错进行迭代计算)。该类方法计算原理简单,容易实现,应用较为广泛,但是计算工作量大且求解精度不高。

本文中利用商业软件COMSOL Multiphysics进行THM耦合模型的数值求解。该软件基于有限单元法,针对多物理场提供了全耦合求解方法,能够在耦合分析中将耦合的各物理场联立起来形成一个统一的耦合方程组进行求解,同时计算得到各个独立场的因变量。这种方法全面考虑了场与场之间的耦合效应,最接近于真实情况[20]。

对复杂裂隙网络的模拟也一直是裂隙岩体THM数值模拟的难点。一般采用具有一定厚度的薄层单元模拟[13-14],但是在考虑大量裂隙时会造成网格数量巨大而难以计算。借鉴文献[8]的思路,将裂隙视作无厚度的线单元,进一步考虑含无厚度裂隙单元的THM耦合分析,并在COMSOL软件中实现。其中基质岩块利用实体单元离散,软件内置了传热、渗流及应力分析的相应模块;裂隙面则采用软件提供的低维单元——Coefficient Form Boundary PDE Interfaces进行二次开发模拟;通过设置岩块单元与裂隙单元之间的物理量交换实现耦合。该方法能够有效模拟EGS岩体中的大规模复杂裂隙网络。

3 耦合模型的验证

通过算例证明本文中的THM全耦合数学模型以及基于COMSOL的耦合数值求解的正确性和有效性。

白冰[22]利用有限Fourier 变换及其逆变换推导了一维条件下饱和土柱热弹性固结问题的解析解,这是一个典型的热流固耦合问题。利用数值方法对该问题进行求解,模型如图1所示。

图1 饱和土热弹性固结计算模型Fig.1 Model for thermal consolidation problem

土柱高度为7 m,在土柱的顶面向下施加载荷10 kPa,土柱初始温度为10 ℃,初始孔隙水压力为10 kPa,顶面受60 ℃的温度载荷作用。模型的边界条件为:左右两侧为完全不透水和绝热边界,x方向位移约束;下边界为绝热不透水边界,z方向位移约束;上表面为渗流自由边界。为了便于与解析解对比,本例中将固体和流体的热力学参数进行等效,取介质的整体热力学参数进行计算。算例中所用的物理参数如下:土体的弹性模量为60 MPa,泊松比为0.4,孔隙率为0.2,热膨胀系数为3.0×10-7K-1,比热容为167.2 J/(kg·K),热传导系数为0.836 W/(m·K),密度为2 000 kg/m3;水的密度为1 000 kg/m3,黏度为1 mPa·s;Biot水力耦合系数αB为1。

耦合分析结果如图2所示。图2中分别显示了土柱内不同位置处位移、孔隙水压力和温度随时间的变化。由图2可知,本文中的数值解和解析解非常吻合,尽管这一模型中没有涉及岩体裂隙,但能从一定程度上说明本文中的数学模型和求解方法可靠。

图2 数值解与解析解的结果对比Fig.2 Comparison between numerical and analytical solutions

4 算例分析

增强型地热系统采用人工压裂的手段产生大量张开的连通裂隙带,以提高热储的渗透率,经济地开采出地热。由于热储位于地下深部岩体,直接获取这些裂隙带的信息进行模拟目前仍然是非常困难的[8]。本文中利用随机生成的裂隙网络描述EGS储层,模拟地热开采中渗流、应力及温度变化过程。

4.1 计算模型

计算模型如图3所示。假定其为300 m×300 m范围的干热岩热储。模型中主要有2组裂隙,倾角分别为30°和110°,裂隙岩体各向异性明显。裂隙迹长采用正态分布,其长度平均值为30 m,方差为10 m。随机生成的裂隙网络能够从统计意义上代表真实的裂隙岩体,反映裂隙岩体THM耦合过程的某些主要特性。

图3 二维EGS裂隙网络模型Fig.3 2D fracture network model of EGS

4.2 计算参数及计算条件

4.2.1 计算参数

数值模拟的可靠性依赖于计算参数的选取。由于缺乏相应的试验资料,参考已有的数值模拟以及通过假定选取部分计算参数。模型中采用的主要材料参数如下:水的密度为1 000 kg/m3,动力黏度为1 mPa·s,比热容为4 200 J/(kg·K),导热系数为0 W/(m·K);岩块密度为2 700 kg/m3,渗透率为1.0×10-18m2,储水率为1.0×10-8Pa-1,弹性模量为30 GPa,泊松比为0.25,比热容为1 000 J/(kg·K),导热系数为3 W/(m·K),热膨胀系数为5.0×10-6K-1,孔隙率为0.000 1;裂隙宽度为0.000 5 m,渗透率为1.0×10-11m2,储水率为1.0×10-9Pa-1,法向刚度为1 200 GPa/m,切向刚度为400 GPa/m;重力加速度为9.8 m/s2;Biot耦合系数为1.0;换热系数为3 000 W/(m2·K)。

采用式(11)描述裂隙面渗流与应力的耦合特性,为便于分析,将影响系数α简单取为0.2×10-6Pa-1。受到拉应力(拉应力为正)时,裂隙张开,渗透系数增大;受到压应力时,裂隙面闭合,渗透系数减小。

4.2.2 初始条件和边界条件

利用上述模型模拟EGS从开始注水到运行的全过程,须满足的边界条件和初始条件主要有:

(1)渗流场。假定左侧边界为进水井,水压力19 MPa;右侧边界为生产井,水压力6 MPa;上、下边界为不透水边界;注水前储层内初始水压力6 MPa。

(2)温度场。进水井温度边界与井内水温一致,取为20 ℃;上、下边界为绝热边界;储层内岩体和水的初始温度均为200 ℃。

(3)应力场。应力场模拟中,模型四周均设为位移约束边界。由于工程上更关心的是EGS开发和运行过程对岩体应力场的扰动,因此本文中不考虑初始地应力场的影响,仅研究高压水注入和地热开采对裂隙岩体应力场产生的影响,有助于理解渗流场和温度场与岩体应力和位移变化之间的耦合作用关系。

算例中采用瞬态模拟求解,计算总时长为40 a,步长取1 d。计算模型的有限元网格总数为36 668,结点总数为18 491。有限元分析在个人平台Intel(R)Core(TM)i7-4790K CPU@4.00 GHz上完成,计算用时约40 min。

4.3 数值模拟结果

4.3.1 渗流场

渗流场的流速直接影响裂隙岩体中温度场的变化,研究渗流场的变化规律对控制ESG的运行有重要意义。图4为不同时刻储层渗流场水压力的分布。在初始阶段,注水井处的低温水在高压作用下流入裂隙中,沿裂隙网络流向生产井。由于裂隙导压能力更强,裂隙中的水压力升高较快,随后基质岩块中的水压力也逐渐升高,与裂隙中的水压趋于一致。经过大约20 a,整个系统内的水压力分布趋于稳定,认为渗流场基本达到了稳定状态,储层内的水在固定的压差作用下循环流动。

4.3.2 温度场

图5为EGS储层内不同时刻温度分布。初始阶段,随着低温水的注入,基质岩块与裂隙水发生热交换,裂隙水温迅速升高。由于热传导作用,距离裂隙面较远的岩块温度逐渐下降,形成低温区。此时,生产井处岩体的温度基本没有发生变化,系统能保证持续的热量输出。随着开采时间不断增长,储层内的低温区也在随之扩大。可以明显看出,在一些分布密集的裂隙通道附近岩体的温度变化更快。这是由于这些贯通的裂隙构成了主要的导水区,其水流速度更高,裂隙水的热对流作用明显。温度场与渗流场耦合的规律与文献[8]和[10]的结果较为相符。

分析结果表明EGS中裂隙岩体温度场分布的强非均匀性和各向异性。利用水力压裂产生裂隙增强岩体的导流能力、提高采热效率是这一技术的重要特征。本文中的方法能够反映出裂隙岩体中热量主要沿裂隙传输的机制,有助于深入认识EGS中的热流耦合作用。

图4 不同时刻储层内的水压分布Fig.4 Pressure distribution in reservoir at different time

图5 不同时刻储层内的温度分布Fig.5 Temperature distribution in reservoir at different time

4.3.3 应力场

水压和温度的变化影响储层内的应力场和位移场分布。图6和7分别为储层内不同时刻的第一主应力和位移等值线图。可以看出,在系统开始运行阶段,由于高压水的注入,注水井附近的温度很快下降,裂隙岩体的应力和位移发生剧烈变化。同时,裂隙面附近第一主应力以拉应力为主,说明高压水进入裂隙使裂隙面张开,导致裂隙渗透系数增大。渗透系数增大说明裂缝渗透阻力下降,更有利于热能的提取[14]。在系统运行后期(20 a以后),水压变化已经不大,此阶段岩体应力场和位移场仍然发生较大变化,这是由于储层的温度分布一直在发生改变,低温区不断扩大,基质岩块发生收缩变形,导致岩体应力场重新分布。数值模拟结果反映了EGS运行中温度、渗流和应力耦合的基本规律。

图6 不同时刻第一主应力等值线图Fig.6 Contours of the maximal principal stress at different time

图7 不同时刻位移等值线图Fig.7 Contours of displacement at different time

在高压水的作用下,裂隙和基质岩块的变形提高了系统的导流能力,加快了热量的传输速率;温度的剧烈变化也会导致岩体应力分布的改变。在EGS储层中,温度、渗流和应力场之间相互影响和制约。

4.3.4 EGS运行情况预测

生产井出口水温是评估EGS出力和寿命研究的重要指标之一[17]。图8为生产井内的水温分布。在20 a以前生产井内水温始终为200 ℃,说明这一阶段系统可以保持稳定出力。随着储层内温度变化,低温区不断扩大,20 a以后生产井附近的温度开始下降,而且生产井内水温并非均匀下降,在裂隙面附近由于裂隙水的流速高,水温下降更快。

为了进一步评估系统的运行及可持续开发能力,利用文献[8]提供的公式计算出口平均水温,

(14)

式中,加和项为裂隙;积分项为基质岩块。

图9为不同入口水温情况下计算得到的出口平均水温随时间的变化曲线。从图9可以看出,3种情况下EGS的稳定运行时间均约为20 a。在20 a到40 a时间内,储层内低温区开始扩大到生产井处,生产井处的出口平均水温逐渐下降,系统的出力也会不断减小。入口水温越低,则出口平均水温下降越快。到40 a时,出口平均水温为岩体初始状态的75%～85%,将降低系统的运行效率。

图8 生产井内的水温变化情况Fig.8 Temperature distribution in production well

图9 出口平均水温与时间的变化关系Fig.9 Relationship between outflow temperature and time

5 结 论

(1)将岩体视作由离散裂隙网络和基质岩块组成的双重介质,有效考虑裂隙岩体传热、渗流与应力过程中的相互作用,建立了相应的THM耦合数学模型。

(2)基于商业软件COMSOL Mutilphysics二次开发,实现了模型的全耦合求解,求解精度更高。对于饱和土热固结问题,能够获得与解析解较为一致的结果,在一定程度上说明了耦合模型的可靠性。

(3)在增强型地热系统的算例分析中,数值模拟可以得到温度、渗流及应力各物理场的分布特征,且和现有的一般性结论较为吻合。数值模拟反映了三场之间的耦合作用规律,说明干热岩开采过程中存在较强的THM耦合作用。研究增强型地热系统的利用效率和运行控制有必要考虑THM三场耦合作用特点。

[1] 许天福,张延军,曾昭发,等.增强型地热系统(干热岩)开发技术进展[J].科技导报,2012,30(32):42-45. XU Tianfu, ZHANG Yanjun, ZENG Zhaofa, et al. Technology progress in an enhanced geothermal system (hot dry rock) [J]. Science & Technology Review, 2012,30(32):42-45.

[2] 许天福,张延军,于子望,等.干热岩水力压裂实验室模拟研究[J].科技导报,2015,33(19):35-39. XU Tianfu, ZHANG Yanjun, YU Ziwang, et al. Laboratory study of hydraulic fracturing on hot dry rock[J]. Science & Technology Review, 2015,33(19):35-39.

[3] 汪集旸,胡圣标,庞忠和,等.中国大陆干热岩地热资源潜力评估[J].科技导报,2012,30(32):25-31. WANG Jiyang, HU Shengbiao, PANG Zhonghe, et al. Estimate of geothermal resources potential for hot dry rock in the continental area of China [J]. Science & Technology Review, 2012,30(32):25-31.

[4] 王晓星,吴能友,苏正,等.增强型地热系统开发技术研究进展[J].地球物理学进展,2012,27(1):15-22. WANG Xiaoxing, WU Nengyou, SU Zheng, et al. Progress of the enhanced geothermal systems (EGS) development technology [J]. Progress in Geophysics, 2012,27(1):15-22.

[5] 胡剑,苏正,吴能友,等.增强型地热系统热流耦合水岩温度场分析[J].地球物理学进展,2014,29(3):1391-1398. HU Jian, SU Zheng, WU Nengyou, et al. Analysis on temperature fields of thermal-hydraulic coupled fluid and rock in enhanced geothermal system [J]. Progress in Geophysics, 2014,29(3):1391-1398.

[6] 张树光,徐义洪.裂隙岩体流热耦合的三维有限元模型[J].辽宁工程技术大学学报(自然科学版),2011,30(4):505-507. ZHANG Shuguang, XU Yihong. 3D finite element model of fluid-heat coupling heat-transfer in fracture rock [J]. Journal of Liaoning Technical University (Natural Science), 2011,30(4):505-507.

[7] 张树光,李志建,徐义洪,等.裂隙岩体流-热耦合传热的三维数值模拟分析[J].岩土力学,2011,32(8):2507-2511. ZHANG Shuguang, LI Zhijian, XU Yihong, et al. Three-dimensional numerical simulation and analysis of fluid-heat coupling heat-transfer in fractured rock mass [J]. Rock and Soil Mechanics, 2011,32(8):2507-2511.

[8] 陈必光,宋二祥,程晓辉.二维裂隙岩体渗流传热的离散裂隙网络模型数值计算方法[J].岩石力学与工程学报,2014,33(1):43-51. CHEN Biguang, SONG Erxiang, CHENG Xiaohui. A numerical method for discrete fracture network model for flow and heat transfer in two-dimensional fractured rocks [J]. Chinese Journal of Rock Mechanics and Engineering, 2014,33(1):43-51.

[9] CHEN B G, SONG E X, CHENG X H, et al. Plane-symmetrical simulation of flow and heat transport in fractured geological media: a discrete fracture model with comsol [C]//Multiphysical Testing of Soils and Shales. Berlin Heidelberg: Springer, 2013:149-154.

[10] XU C, DOWD P A, TIAN Z F. A simplified coupled hydro-thermal model for enhanced geothermal systems [J]. Applied Energy, 2015,140:135-145.

[11] LEI H, XU T, JIN G. TOUGH2Biot—a simulator for coupled thermal-hydrodynamic-mechanical processes in subsurface flow systems: application to CO2geological storage and geothermal development [J]. Computers & Geosciences, 2015,77:8-19.

[12] RUTQVIST J. Status of the TOUGH-FLAC simulator and recent applications related to coupled fluid flow and crustal deformations [J]. Computers & Geosciences, 2011,37(6):739-750.

[13] BAHRAMI D, DANKO G, FU P, et al. Poroelastic and self-propped single fracture THM models for EGS studies [C]//Proceedings of Fourtieth Workshop on Geothermal Reservoir Engineering. California: Stanford University Press, 2015.

[14] 赵阳升,王瑞凤,胡耀青,等.高温岩体地热开发的块裂介质固流热耦合三维数值模拟[J].岩石力学与工程学报,2002,21(12):1751-1755. ZHAO Yangsheng, WANG Ruifeng, HU Yaoqing, et al. 3D numerical simulation for coupled THM of rock matrix-fractured media in heat extraction in HDR [J]. Chinese Journal of Rock Mechanics and Engineering, 2002,21(12):1751-1755.

[15] ZHAO Y, FENG Z, FENG Z, et al. THM (thermo-hydro-mechanical) coupled mathematical model of fractured media and numerical simulation of a 3D enhanced geothermal system at 573 K and buried depth 6000—7000 M [J]. Energy, 2015,82:193-205.

[16] 丁军锋,王世民.裂隙岩石变形-渗流-传热耦合的三维数值模拟:2014年中国地球科学联合学术年会论文集[C].北京：地质出版社，2014.

[17] 赵阳升.多孔介质多场耦合作用及其工程响应[M].北京:科学出版社,2010.

[18] LOUIS C. Rock hydraulics in rock mechanics [M]. New York:Springer Verlag, 1974:299-387.

[19] 沈振中,徐志英,雒翠.三峡大坝坝基粘弹性应力场与渗流场耦合分析[J].工程力学,2000,17(1):105-113. SHEN Zhenzhong, XU Zhiying, LUO Cui. Coupled analysis of viscoelasticity stress field and seepage field of the three Gorges Dams foundation [J]. Engineering Mechanics, 2000,17(1):105-113.

[20] 吴持恭.水力学(上册)[M].北京:高等教育出版社,1982.

[21] 孙培德,杨东全,陈奕柏.多物理场耦合模型及数值模拟导论[M].北京:中国科学技术出版社,2007.

[22] 白冰.岩土颗粒介质非等温——维热固结特性研究[J].工程力学,2005,22(5):186-191. BAI Bing. One-dimentional thermal consolidation characteristics of geotechnical media under non-isothermal condition [J]. Engineering Mechanics, 2005,22(5):186-191.

(编辑 李志芬)

A thermo-hydro-mechanical coupling model for numerical simulation of enhanced geothermal systems

SUN Zhixue1, XU Yi2, LÜ Shuhuan1, XU Yang1, SUN Qiang1, CAI Mingyu1, YAO Jun1

(1.School of Petroleum Engineering in China University of Petroleum, Qingdao 266580, China;2.StateKeyLaboratoryofWaterResourcesandHydropowerEngineeringScience,WuhanUniversity,Wuhan430072,China)

In the enhanced geothermal system (EGS), an artificial fracture network can be formed in the geo-thermal reservoir by hydraulic fracturing for heat transmission fluids to be circulated in the reservoir to extract heat from hot dry rocks (HDR), which involves complex thermo-hydro-mechanical (THM) coupling processes in the fractured rock matrix. Therefore, THM coupling analysis is of great importance for studying the performance of the EGS. A THM coupling model was presented in this paper, in which the fractured rock was regarded as a dual medium (discrete fracture networks and matrix rock), and a fully-coupling finite element analysis was implemented in a commercial software, COMSOL multiphysics, for analysis of the temperature, pressure and stress distributions in the fractured rock system during geo-thermal recovery. The coupling model and the numerical approach were verified in comparison with a proven analytical solution method. The model was used to simulate the EGS process with a 2D randomly generated fracture network to study the characteristics of flow, heat transfer and mechanical behaviors in a HDR reservoir. The results show that the connected fractures in the reservoir form the main flow pathways for fluid circulation, in which heat convection is essential for heat transfer. High pressure water injection and fracture extension induced by temperature variation can increase the conductivity of the fractured rock matrix, thus affecting the heat extraction. It is important to consider the THM coupling effects for studying the efficiency and performance of the EGS process.

hot dry rock; enhanced geothermal system; fractured rock matrix; thermo-hydro-mechanical coupling; numerical simulation

2015-12-23

国家自然科学基金项目(51404291);山东省青年基金项目(2013ZRE28068); 中央高校基本科研业务费专项(14CX05024A,14CX02045A);国土资源部重点实验室开放基金项目(SHW[2014]-DX-11)

孙致学(1979-),男,副教授,博士,研究方向为地热资源及油气田开发工程。E-mail:szx1979@126.com。

姚军(1964-),男,教授,博士,博士生导师,研究方向为油气田开发工程。E-mail:youcang@upc.edu.cn。

1673-5005(2016)06-0109-09

10.3969/j.issn.1673-5005.2016.06.014

TK 529

A

孙致学,徐轶,吕抒桓,等. 增强型地热系统热流固耦合模型及数值模拟[J].中国石油大学学报(自然科学版),2016,40(6):109-117.

SUN Zhixue, XU Yi, LÜ Shuhuan, et al. A thermo-hydro-mechanical coupling model for numerical simulation of enhanced geothermal systems[J].Journal of China University of Petroleum(Edition of Natural Science),2016,40(6):109-117.