近年来，随着地球能源开发消耗，环境污染日益加剧，大力开发可再生绿色能源是必然趋势。蕴藏在深部高温裂隙岩体中的地热能将逐渐成为21世纪重要的可再生绿色能源[1－2]。以开采深部岩体地热能为目标的地热对井系统受到广泛重视[3]。地热对井系统运作的基本原理是向储层注入井中灌入低温流体介质，置换出深部岩体中的热能，并由生产井中提取。地热开采过程涉及低温流体与相邻高温裂隙岩体间相互作用，存在复杂的渗流-传热耦合过程[4]：流体运动会影响裂隙及其与相邻岩块间的热传导和热对流行为；同时，温度场变化也会影响流体黏度，甚至改变流体性态，进而反作用于渗流场。由于耦合作用的复杂性及裂隙岩体的高度非均匀性，考虑地热开采中渗流与传热行为的数值模拟是分析热储工程的重要手段。
针对地热开采过程中涉及的裂隙岩体渗流与传热耦合问题，目前采用的数值模型主要包括等效连续介质模型和离散裂隙网络模型。等效连续介质模型视裂隙和岩块为同一均匀介质进行模拟，渗流和传热效应均匀作用于整个研究域[5－7]，忽略岩体中裂隙的复杂分布特征。尽管等效介质模型简单便捷，但其无法描述流体和热能在不同介质(裂隙和岩块)间的传输和交换。离散裂隙网络模型假设岩体由岩块和裂隙网络组成，显式地考虑每一条裂隙的贡献，更接近于实际储层的裂隙岩体结构，能够较为真实地模拟裂隙中的水热迁移过程，引起了广泛的研究兴趣。如张树光等[8]建立了裂隙岩体的流-热耦合数学模型，并采用实体单元模拟单条裂隙，分析了岩体与水流之间的传热过程；陈必光等[9]研究了岩石基质与裂隙水界面热交换，提出了基于离散裂隙网络的二维渗流传热耦合数学模型，其中裂隙采用无厚度单元模拟；Xu等[10]提出了一种简化方法来模拟对井系统地热开采中复杂裂隙岩体的流-热耦合过程；薛娈鸾[11]基于复合单元法建立了裂隙岩体渗流-传热耦合的复合单元模型，采用交叉迭代算法，对含单条裂隙岩体的渗流场和温度场进行耦合分析；黄诗冰等[12]建立了低温冻结条件下裂隙岩体水-热耦合模型，评价了裂隙网络渗流对冻结交圈的影响。目前针对裂隙岩体热流耦合数值模拟已开展了大量有益的探索，但多集中在二维模型或少量裂隙(如单条裂隙)情况。对于地热对井系统，由于相较于储层研究区域，裂隙开度及对井结构孔径尺寸极小，且裂隙数量巨大且分布密集，采用三维离散裂隙网络模型模拟渗流传热耦合过程仍然十分困难。主要难点在于如何有效考虑裂隙介质、岩块及对井中渗流传热交换，以及解决三维大规模复杂裂隙网络模拟中网格离散难度大且计算收敛性差等问题。
本文提出裂隙岩体三维渗流传热耦合的等效模拟方法，将岩体视为由岩块基质和离散裂隙网络组成的双重介质，基于渗流传热耦合理论和离散裂隙网络模型，采用实体单元模拟岩块、无厚度平面单元模拟裂隙网络以及一维线单元模拟对井，推导了相应的有限元格式，并基于商业有限元软件COMSOL Multiphysics进行二次开发，实现热流耦合数学模型的全耦合求解。通过与精细模拟方法相比较，验证方法的有效性。最后将该法应用于含大规模复杂裂隙网络地热对井系统的数值模拟中。
1 热流耦合数学模型
地热对井系统几何模型构成主要包括：注入井、生产井和富含裂隙的地下储层。采用裂隙-孔隙双重介质模型进行模拟，其中包含3个主要的物理过程，裂隙、基岩内的流体运动，裂隙内热量交换与传输，以及周围岩体内的温度变化过程，分别采用相应的控制方程进行描述。由于深部储层中的压力足够高，液体可以避免发生汽化。
1.1 控制方程
参照文献[9—13]中渗流和传热耦合的基本理论，裂隙岩体中热流耦合数学模型的控制方程如下：
1) 岩块基质渗流场方程
式中：p/Pa为水压力；t/s表示时间；u/(m/s)为岩块中水的流速；S/Pa-1为岩块储水系数；Q/s-1为渗流的源汇项；η/(Pa · s )为流体动力粘度；ρf/(kg/m3)为水的密度；g/(m/s2)为重力加速度；κ/m2为岩块渗透率。
2) 裂隙内渗流场方程
式中：Sf/Pa-1为裂隙的储水系数；κf/m2为裂隙的渗透率；df/m为裂隙的开度；Qf/(m/s)表示从岩块流入裂隙的流量，可按下式计算：
其中，n表示裂隙面法向。
真实的岩体裂隙可看作是由两个具有相同或不同粗糙度的壁面叠合而成，其内部流体运动一般认为服从立方定律[9,12]，则渗透率可表示为：
式中：ξ是粗糙度修正系数。
3) 岩块和裂隙温度场方程
式中：T/K为温度；ρ/(kg/m3)为密度；λ/(W/(m·K))为热传导系数；c/(J/(kg·K))为岩块热容。下标m和f分别代表基岩和裂隙；W/(W/m3)为热源；Wf/(W/m2)为单位面积上裂隙表面水从基岩吸收的热量，表达式如下[10]：
式中，h/(W/(m2·K))为换热系数。
1.2 裂隙岩体渗流-传热耦合机制
裂隙岩体中流体及热能的主要传输通道是裂隙网络。结合式(3)和式(7)可知，裂隙岩体中存在渗流时，流体流速的变化会通过热对流作用影响岩体的温度场。同时，温度场的改变又会导致流体性质发生改变，引起裂隙渗透性的变化，以致研究域渗流场的改变。
裂隙中流体的运动黏度与温度相关，即[10,12]：
式中：v/(×10- 6 m2/s)为流体运动粘度，v=η/ρf；T/(℃)为流体温度。
利用热流耦合模型不仅可计算裂隙岩体的渗流-传热过程，而且还可反映裂隙与相邻岩块间的渗流量和热量交换规律。
2 等效模拟方法
本文提出等效模拟方法对地热对井裂隙岩体渗流传热耦合过程进行仿真。基于裂隙岩体随机离散裂隙网络模型，采用无厚度平面单元模拟大规模裂隙、一维线单元模拟对井。与精细模拟相比，该法可以有效降低模型的有限元前处理网格剖分难度，简化求解过程(如图1所示)。
图1 精细模拟与等效模拟示意图
Fig.1 Illustration of refined and equivalent numerical method
2.1 随机离散裂隙网络模拟
由于地热储层位于地下深部岩体，目前直接获取这些岩体裂隙的具体分布信息仍非常困难[9,14]。因此，根据裂隙几何特征的统计规律随机生成离散裂隙网络是精细模拟深部裂隙岩体的有效途径。本文采用文献[15]中的方法，根据裂隙迹长、间距、密度、开度和走向等信息建立相应的裂隙几何特征统计模型，再通过蒙特卡罗方法随机生成三维裂隙网络模型。
离散裂隙网络模型的显著特点是可以精细地考虑裂隙对整个研究域渗流传热的贡献。由于岩体中裂隙数量巨大，且相对于其他两个方向来说，其开度方向尺寸极小，几乎可以忽略不计。但裂隙开度直接影响其渗透性，类似于混凝土中的裂纹网络[16]及界面过渡区结构[17]，裂隙网络作为岩体中的薄弱介质，具备高渗透性及低强度材料属性，在模拟中必须予以充分考虑。一般采用特殊的单元来描述裂隙，例如扩展单元[18]、复合单元[11]、嵌入单元[19－20]以及无厚度单元等[9]。
本文在已有研究[16－17]的基础上，将无厚度单元法进一步推广到裂隙岩体渗流传热耦合的三维数值模拟中。该法将裂隙开度隐含到渗流控制方程中，有限元前处理中无需考虑裂隙开度，可减小有限元网格剖分难度和降低计算量，在模拟大规模裂隙网络时具有优势。
2.2 对井模拟
对于地热对井系统中的注入井和生产井，由于其尺寸相对较小，并且与离散裂隙网络交叉形成复杂的几何拓扑结构，这会给模型的有限元网格离散带来较大困难。为此，本文提出采用线单元模拟对井来解决这一问题。
图2给出了线单元模拟对井的示意图。图中实线所示为线单元，虚线所示为井筒的实际尺寸，井筒半径为r。
图2 线单元模拟对井的示意图
Fig.2 Modeling of inlet/outlet wells using line elements
类似于大体积混凝土中冷却水管的模拟方法[21]，本文提出利用线单元进行等效模拟，可以将对井井筒离散模拟到计算域中，同时井筒半径隐含在线单元内，减小了网格离散的难度。该法将对井井筒作为具有一定渗透和传热特性的介质考虑到有限元模型中，同时将其与周围介质的热流交换以边界源项的方式加入有限元计算中，即：
式中：Ql/(m2/s)表示从岩块流入井筒的流量；q/s-1为流量边界源项；Wl/(W/m)表示从岩块流入井筒的热通量；f/(W/m3)为热通量边界源项。
需要指出的是，由于井筒半径相对于模型计算尺度而言尺寸极小，井筒与周围介质的热流交换速度较快，二者之间的压力及温度迅速趋于一致，可以采用线单元进行有效模拟。同时，在该模型中井壁厚度忽略不计，因此本文的模拟中未考虑井壁的材料特性及其影响。
2.3 有限元离散格式
本文采用有限单元法来求解裂隙岩体的热流耦合控制方程。空间离散采用一般的伽辽金加权余量法，时间离散采用向后欧拉法(Backward Euler)。上述方程的有限元离散格式为：
其中：
式中：Ω是岩块计算域；γ是裂隙计算域；Γ是计算域边界；N是岩块域的形函数；Nf是裂隙域的形函数；
是边界外法向流量；
是边界外法向热通量。
2.4 算法实现
等效模拟方法将裂隙视作无厚度单元，将对井视作线单元，岩块采用实体单元模拟，进一步考虑裂隙岩体热流耦合分析及岩块与裂隙之间的热流交换。本文基于商业有限元软件COMSOL进行二次开发，实现含大规模裂隙岩体热流耦合模型的数值求解。结合渗流传热控制方程、裂隙和对井模拟方法以及有限元求解算法，岩块单元利用软件内置的传热和渗流模块进行分析；裂隙面单元和井单元分别采用软件提供的低维单元—Coefficient Form boundary PDE Interfaces和Coefficient Form edge PDE Interfaces进行二次开发来模拟；通过设置岩块单元与裂隙单元之间的物理量交换来实现耦合。该方法完全考虑了场与场之间的耦合效应，因此最接近于岩体渗流场和温度场真实的耦合作用情况。
3 算例验证
3.1 三维交叉裂隙渗流模拟
本文采用无厚度单元对三维裂隙网络进行了平面等效，为验证该法对三维渗流模拟的适用性，采用如图3所示的交叉裂隙模型进行算例考核。该模型尺寸为0.15 m×0.15 m×0.15 m，其内部含有5条交叉裂隙。假定裂隙开度均为di，其等效水力开度可通过下列公式计算[22]：
式中：deq/m表示等效水力开度；θi为第i条裂隙与水流方向的夹角；Li/m为第i条裂隙的长度；L/m为裂隙网络沿压力梯度方向的总长度。该模型中，前1/3段和后1/3段都为两条斜裂隙并联，中间1/3段为单条水平裂隙，将式(20)和式(21)结合即可得到交叉裂隙的等效水力开度。不考虑裂隙面粗糙度，可根据立方定律得到该模型的渗透率解析解。
图3 三维交叉裂隙岩体几何模型 /m
Fig.3 Geometry of rock mass with crossed fractures
采用数值方法模拟通过该交叉裂隙模型的稳定渗流过程，假定岩块基本不透水，取其渗透率为裂隙的1.0×10-10，模型左(x=0 m)、右(x=0.15 m)边界分别施加1500 Pa和0 Pa的水压力，其余边界为不透水边界。利用数值积分获得通过模型左右边界的流量，即可计算得到含交叉裂隙岩体的等效渗透率。考虑不同裂隙开度，预测该交叉裂隙模型的等效渗透率，并与解析解进行对比，结果如图4所示。可以看出，数值预测的结果与解析解完全吻合，说明了本文模型对复杂三维渗流模拟的适用性。
图4 数值结果与解析解的对比
Fig.4 Comparison of numerical and analytical results.
3.2 单裂隙热流耦合模拟
为了验证所建立方法在热流耦合模拟中的有效性和可靠性，本算例将等效模拟与精细模拟方法进行对比分析。选取15.0 m×7.5 m×5.0 m(长×宽×高)的长方体储层，其中含0.0005 m厚的水平裂隙以及半径为0.1 m的注入井和生产井，井深为4.0 m，图5给出了计算几何模型。材料参数为水的密度1000 kg/m3，粘度0.001 Pa·s，比热容4200 J/(kg·K)，导热系数0 W/(m·K)；岩块的密度2700 kg/m3，渗透率1.0×10-16 m2，比热容1000 J/(kg·K)，导热系数3 W/(m·K)；裂隙开度0.0005 m，粗糙度系数0.00048，根据式(5)得渗透率为1.0×10-11 m2。计算初始条件和边界条件为：① 渗流场：储层内初始水压力5 MPa；假定左侧为注入井，水压力21 MPa；右侧为生产井，水压力5 MPa；模型外部表面为不透水边界。② 温度场：储层内岩体和水的初始温度均为200 ℃；注入井温度边界与井内水温一致，取20 ℃；生产井为自由流出边界；模型外部表面为绝热边界。采用瞬态求解，计算总时长为30 d，时间步长取0.1 d。
两种计算方法有限元网格见图6，精细模拟中裂隙和对井均采用实体单元模拟，单元数为184076，结点数为34946；等效模拟模型单元数为58758，结点数为10612。可以看出，采用等效模拟方法的网格整体要比精细模拟的更加规则，网格剖分难度降低，网格数量、计算量和计算时长大幅减少。
图5 单裂隙岩体几何模型Fig.5 Geometry model of rock mass with a single fracture
图6 精细模拟和等效模拟方法有限元网格
Fig.6 FE meshes for refined and equivalent numerical method
图7给出了不同时刻生产井深度方向及裂隙面水平中心线方向上的温度分布规律。可以看出，等效模拟方法的计算结果与精细模拟结果较为吻合，说明了等效模拟方法的有效性和可靠性。而两种结果存在偏差的原因是：① 图7(a)和图7(b)中等效模拟结果分别为生产井线单元上和裂隙面单元中心线上的温度分布，而精细模拟则为井筒外壁温度场、裂隙上下表面温度场取平均的结果，因而二者存在一定差异，在井筒与裂隙面交叉位置z=2.5 m和x=12.5 m处差异最为明显；② 由于对井尺寸较大且网格并不密集，若将对井尺寸减小或网格进一步加密，两种方法计算偏差会有所减小。
图8给出了10 d时储层中心剖面上温度场分布规律。可以看出，随着低温水的注入，裂隙温度迅速降低。裂隙面上温度变化较快，裂隙水的对流传热作用明显。由于热传导作用，离裂隙面较近的岩块温度逐渐下降，出现低温区。10 d时生产井处岩体的温度开始发生变化。
图7 不同时刻温度分布规律
Fig.7 Temperature distribution at different elapsed times
图8 10 d时储层中心剖面温度场分布规律
Fig.8 Temperature distribution on the central profile of geothermal reservoir at 10 d
总的来说，两种方法的计算结果较为一致，说明了本文提出的等效方法对于地热对井系统的模拟是可行和可靠的。该法基于渗流-传热耦合模型，采用无厚度单元模拟裂隙、线单元模拟对井，可以有效地模拟裂隙中的渗流传热过程以及裂隙与岩块之间的热流交换。
4 地热对井三维热流耦合模拟
采用实际工程进行验证是检验数值模型可靠性的重要手段。但是，目前深部高温地热资源开采难度大，技术仍不成熟，对地热对井系统的研究还处于理论研究阶段，相关工程实践资料缺乏。实际高温地热岩体中裂隙数量巨大，且分布较密集，直接获取裂隙分布信息及岩体物理力学参数进行建模仍是非常困难的。本节通过生成随机裂隙网络模型来模拟含大规模裂隙的地热对井系统，以进一步探讨等效模拟方法在实际工程中应用的可行性及有效性。随机裂隙网络模型可以反映岩体中裂隙分布的高度非均匀性和各向异性，从统计意义上表征裂隙岩体的物理力学特性。通过仿真大规模裂隙岩体中的渗流传热行为，研究离散裂隙网络和岩块中的渗流-传热耦合作用机制及其主要影响因素，可以为预测地热系统的运行状态以及优化深部地热能开发方案提供一定的理论依据。
以某裂隙岩体地热对井系统为例，其几何尺寸为200.0 m×200.0 m×100.0 m(长×宽×高)，包含随机生成的裂隙网络。其中裂隙总数为600条，裂隙平均迹长为10.0 m，服从方差为2.0 m2的对数正态分布。裂隙随机分布在储层中，相互交叉和搭接形成贯通的裂隙网络。图9给出了地热对井几何模型及有限元网格，其中单元总数为2476348，结点总数为426604。岩块渗透率为1.0×10-18 m2；裂隙开度为0.0005 m，粗糙度系数为0.024，渗透率为5.0×10-10 m2，其余材料参数、计算初始条件和边界条件与第3节相同，算例中采用瞬态求解，计算总时长为30 a，时间步长取0.1 a。
1) 储层内温度场分布结果
图9 地热对井裂隙岩体几何模型及有限元网格
Fig.9 Geometry model and FE meshes of the fractured rock mass of doublet system
图10给出了不同时刻地热储层中心剖面上温度场分布规律。从中可以看出，初始阶段(1 a)随着低温水的注入，注入井与其周围裂隙网络和岩块之间快速发生热量交换，裂隙中水温迅速降低，形成低温区；离注入井较远的裂隙网络处温度也开始降低，但此时，生产井处岩体温度基本没有发生变化，系统能保证持续的热量输出。随后(5 a ～20 a)，储层中的低温区不断扩大，尤其是在贯通裂隙网络通道上温度降低较快，说明了贯通裂隙网络是主要的导水通道，水流速度高，进一步加剧裂隙水的热对流作用，导致冷锋面沿裂隙贯通通道发生移动，同时可以明显观察到过早热突破和长尾效应，这与文献[9]中的二维模型计算结果较为一致；生产井处岩体温度也开始发生改变，热量输出逐渐减小。总的看来，由于裂隙网络的非均匀及各向异性分布，导致了裂隙岩体温度场分布的不均匀，裂隙中的热对流作用是影响温度场分布的重要因素，采用本文方法可以有效反映裂隙岩体中的渗流传热耦合机制。
图10 不同时刻地热储层中心剖面温度场分布规律
Fig.10 Temperature distribution on the central profile of geothermal reservoir at different elapsed times
2) 储层内平均温度及采热率
储层内平均温度及采热率是地热系统开采性能评估的重要指标[23]。局部热开采率ζL定义为[24]：
式中：Tr/K和Ts(t)/K分别为初始时刻和t时刻岩块温度；Tinj/K 为注入流体的温度。
整体采热率ζ指在时刻t，热储中已开采热能与可利用总热能的比值：
图11 不同裂隙开度对储层平均温度和整体采热率的影响
Fig.11 Effects of fracture apertures on average temperature and heat extraction ratio
图11给出了不同裂隙开度情况下，地热对井裂隙岩体中平均温度及整体开采率的变化过程。可以看出，平均温度及整体采热率与时间分别呈负和正相关关系。随着裂隙开度的增加，储层平均温度下降，整体采热率不断提高。同时，二者的变化速度随裂隙开度的增大而减小；裂隙开度越大，岩体整体渗透率增加，热能开采的效率大幅提高；当裂隙开度进一步增大，整体采热率增幅减小。这些结果对于进一步优化开采条件及热储建造设计具有一定参考价值。
5 结论
本文提出了地热对井裂隙岩体三维渗流传热耦合的等效模型方法。通过与解析方法和精细模拟方法相比较，验证本文方法的有效性，进一步应用于预测含大规模裂隙岩体中的渗流及传热耦合行为，为深部地热资源的研究和开发利用提供一定参考。
(1) 提出了地热对井裂隙岩体三维渗流传热耦合的等效模拟方法，地热对井视为由离散裂隙网络、岩块和对井组成的系统，采用无厚度平面单元模拟离散裂隙、一维线单元模拟对井，实现热采过程的数值模拟。该算法可以大幅减小前处理网格剖分难度和降低计算量，在含大规模裂隙网络的岩体研究中具有较大优势和应用前景。
(2) 开展算例分析将等效模拟方法与解析方法及精细模拟方法进行了对比，均获得了较为一致的结果，在一定程度上说明了渗流传热耦合模型和等效模拟方法的有效性和可靠性。
(3) 以大规模裂隙岩体地热对井系统为研究背景进行仿真分析，结果表明：由于裂隙网络的非均匀及各向异性分布，导致了裂隙岩体温度场分布出现高度不均匀；冷锋面沿储层内的主要贯通裂隙网络移动，裂隙中对流作用明显；裂隙开度是影响温度场分布的重要因素。采用本文方法可以有效反映裂隙岩体中渗流传热耦合机制。
需要补充的是，本文提出的地热对井裂隙岩体渗流传热耦合等效模拟方法未涉及岩体变形，地热储层开发过程中的流固热三场耦合效应及其引起的裂隙岩体失稳及井壁破坏等问题仍有待进一步研究。
[1]Olsthoorn D, Haghighat F, Mirzaei P A.Integration of storage and renewable energy into district heating systems: A review of modelling and optimization [J].Solar Energy, 2016, 136: 49－64.
[2]Wei G, Meng J, Du X, et al.Performance analysis on a hot dry rock geothermal resource power generation system based on kalina cycle [J].Energy Procedia, 2015,75: 937－945.
[3]Tester J W, Anderson B, Batchelor A, et al.The future of geothermal energy: Impact of enhanced geothermal systems (EGS) on the United States in the 21st century[R].Cambridge, MA, USA: Massachusetts Institute of Technology, 2006.
[4]Samardzioska T, Popov V.Numerical comparison of the equivalent continuum, non-homogeneous and dual porosity models for flow and transport in fractured porous media [J].Advances in Water Resources, 2005,28(3) : 235－255.
[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.(in Chinese)
[6]Jiang F M, Chen J L, Huang W B, et al.A three-dimensional transient model for EGS subsurface thermo-hydraulic process [J].Energy, 2014, 72: 300－310.
[7]Zeng Y C, Su Z, Wu N Y.Numerical simulation of heat production potential from hot dry rock by water circulating through two horizontal wells at Desert Peak geothermal field [J].Energy, 2013, 56(63): 92－107.
[8]张树光, 李志建, 徐义洪, 等.裂隙岩体流-热耦合传热的三维数值模拟分析[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.(in Chinese)
[9]陈必光, 宋二祥, 程晓辉.二维裂隙岩体渗流传热的离散裂隙网络模型数值计算方法[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.(in Chinese)
[10]Xu C S, Dowd P A, Tian Z F.A simplified coupled hydro-thermal model for enhanced geothermal systems[J].Applied Energy, 2015, 140: 135－145.
[11]薛娈鸾.裂隙岩体渗流-传热耦合的复合单元模型[J].岩土力学, 2016, 37(1): 263－268.Xue Luanluan.A composite element model for coupled seepage-heat transfer of fractured rock mass [J].Rock and Soil Mechanics, 2016, 37(1): 263－268.(in Chinese)
[12]黄诗冰, 刘泉声, 程爱平, 等.低温裂隙岩体水-热耦合模型研究及数值分析[J].岩土力学, 2018, 39(2):735－744.Huang Shibing, Liu Quansheng, Cheng Aiping, et al.A coupled hydro-thermal model of fractured rock mass under low temperature and its numerical analysis [J].Rock and Soil Mechanics, 2018, 39(2): 735－744.(in Chinese)
[13]Xu C, Dowd P A, Zhao F T.A simplified coupled hydro-thermal model for enhanced geothermal systems[J].Applied Energy, 2015, 140: 135－145.
[14]李鹏飞, 朱其志, 顾水涛, 等.岩石类材料裂隙形成和扩展的相场方法模拟[J].工程力学, 2018, 35(3): 41－48.Li Pengfei, Zhu Qizhi, Gu Shuitao, et al.A phase field method to simulate crack nucleation and crack propagation in rock-like materials [J].Engineering Mechanics, 2018, 35(3): 41－48.(in Chinese)
[15]Thovert J F, Mourzenko V V, Adler P M.Percolation in three-dimensional fracture networks for arbitrary size and shape distributions [J].Physical Review E, 2017,95(4): 042112.
[16]Li X X, Chen S H, Xu Q, et al.Modeling capillary water absorption in concrete with discrete crack network [J].Journal of Materials in Civil Engineering, 2017, 30(1):04017263.
[17]Li X X, Xu Y, Chen S H.Computational homogenization of effective permeability in three-phase mesoscale concrete [J].Construction and Building Materials, 2016,121: 100－111.
[18]Berrone S, Pieraccini S, Scialo S.On simulations of discrete fracture network flows with an optimization-based extended finite element method [J].SIAM Journal on Scientific Computing, 2013, 35(2):A908－A935.
[19]钱鹏, 徐千军.基于单元嵌入技术和弹性比拟的含裂纹混凝土三维渗流模拟方法[J].工程力学, 2017,34(4): 125－133.Qian Peng, Xu Qianjun.Three-dimensional seepage analysis for cracked concretes based on embedded elements and elastic analogy [J].Engineering Mechanics,2017, 34(4): 125－133.(in Chinese)
[20]钱鹏, 徐千军.不同裂纹分布的孔隙材料渗透系数[J].工程力学, 2017, 34(12): 39－47.Qian Peng, Xu Qianjun.Permeability of porous material with different crack distributions [J].Engineering Mechanics, 2017, 34(12): 39－47.(in Chinese)
[21]张超, 段寅, 刘杏红, 等.基于并层单元的大体积混凝土水管冷却温度场热-流耦合精细计算[J].工程力学,2014, 31(12): 147－154.Zhang Chao, Duan Yin, Liu Xinghong, et al.The precise heat-fluid coupling method of mass concrete with cooling pipes based on layer-merged element [J].Engineering Mechanics, 2014, 31(12): 147-154.(in Chinese)
[22]Sarkar S, Toksöz M N, Burns D R.Fluid flow modeling in fractures [R].Cambridge: Massachusetts Institute of Technology, Earth Resources Laboratory, 2004.
[23]Sanyal S K, Butler S J.An analysis of power generation prospects from enhanced geothermal systems [C].Geothermal Resources Council Transactions 2005,Antalya, Turkey, April 24－29, 2005.
[24]Jiang F M, Chen J L, Huang W B, et al.A three-dimensional transient model for EGS subsurface thermo-hydraulic process [J].Energy, 2014, 72: 300－310.