基于立方定律的断层流—热耦合数值计算方法

doi:10.11872/j.issn.1005-2518.2020.06.122

基于立方定律的断层流—热耦合数值计算方法

陈刚^,¹, 马玲^,², 龚红胜¹

1.昆明理工大学国土资源工程学院，云南昆明 650031

2.昆明理工大学城市学院，云南昆明 650051

Numerical Calculation Method of Fault Flow-Thermal Coupling Based on Cubic Law

CHEN Gang^,¹, MA Ling^,², GONG Hongsheng¹

1.Faculty of Land Resources Engineering，Kunming University of Science and Technology，Kunming 650031，Yunnan，China

2.City College，Kunming University of Science and Technology，Kunming 650051，Yunnan，China

通讯作者: 马玲（1981-），女，云南新平人，工程师，从事水文地质和地下水数值模拟研究工作。maling@kust.edu.cn

收稿日期: 2020-07-09 修回日期: 2020-08-03 网络出版日期: 2021-01-29

基金资助:

国家自然科学基金项目“基于裂隙三维空间分布的矿区地下水流动模拟研究”. 41562017

Received: 2020-07-09 Revised: 2020-08-03 Online: 2021-01-29

作者简介 About authors

陈刚（1981-），男，河南林州人，博士研究生，讲师，从事岩体裂隙渗流和成矿动力学模型研究工作chen_kust@qq.com , E-mail：chen_kust@qq.com

摘要

由于断层宽度远小于其延伸方向的尺寸，造成数值模型建模困难、计算效率低等问题。把断层概化为无几何厚度的空间曲面，并引用裂隙渗流理论中的立方定律对其进行渗流计算，可以有效降低建模难度。本研究旨在验证断层概化方法的可行性和合理性，并解决断层中的裂隙流与基岩中的达西流之间的流—热耦合问题。采用公式推导得出断层与基岩之间的流—热耦合控制方程，并使用数值模型的计算结果进行了验证分析。结果表明：当不考虑断层内部结构影响时，利用立方定律计算断层渗流的方法是可行的；耦合控制方程对基岩与断层之间的流—热耦合计算合理有效；随温度而变化的流体粘滞性对数值模型计算结果影响显著。

关键词： 立方定律 ; 断层 ; 流—热耦合 ; 数值计算 ; 热液型矿床 ; 裂隙渗流

Abstract

For the ore-forming process of hydrothermal deposits，the seepage of fluids in the rock matrix and fissures （faults） produces material and energy transmission，and forms orebodies at specific locations with changes in temperature and pressure.Because the width of the fault is much smaller than the dimension of its extension direction，it causes problems such as difficulty in modeling numerical models and low calculation efficiency.According to the geometric characteristics of the fault，it can be generalized into a space surface to reduce the difficulty of modeling.The generalized fault uses the cubic law in rock mass fracture seepage theory to calculate the fault seepage problem.The seepage of hydrothermal fluid is not limited to faults，but also occurs in bedrock，and this process is calculated using Darcy’s law. The fracture flow in the fault and the Darcy flow in the bedrock interact with each other.In order to ensure the continuity of the pressure，velocity，mass，and energy of the seepage field in the numerical model calculation domain，the flow-heat coupling calculation is required.The purpose of this study is to verify the feasibility and rationality of the generalization method of the fault space surface，and to solve the problem of flow-heat coupling between the fissure flow in the fault and the Darcy flow in the bedrock.The viscosity of fluid has the property of changing with temperature.This article will discuss whether the change of viscosity has an effect on the calculation result of the numerical model initially.Based on the theoretical formula of cubic law，the formula is derived according to the characteristics of small fault thickness，and the flow-heat coupling control equation of fracture flow and Darcy flow is obtained. In order to verify the rationality of the control equations，numerical model experiments are used for verification and analysis.After analysis，it is considered that the method of calculating the seepage of the fault using cubic law is feasible when the internal structure of the fault is not taken into consideration，which can reduce the difficulty of modeling the numerical model.Because the fault uses a spatial surface，the reduction of the dimension compared to the overall model also brings increased computing efficiency. After analyzing the results of the numerical experiments，it is considered that the coupling control equation is reasonable and effective for the calculation of the flow-heat coupling between the bedrock and the fault，which is in accordance with the laws of seepage and heat conduction.Based on the original experimental numerical model，a model in which the viscosity coefficient of the fluid does not change with temperature is established，and the change curves of mass and heat conduction flux are compared.It is found whether the change of the fluid viscosity is considered to have a significant effect on the calculation result of the numerical model.

Keywords： cubic law ; fault ; flow-thermal coupling ; numerical calculation ; hydrothermal deposit ; fissure seepage

PDF (4085KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

陈刚, 马玲, 龚红胜. 基于立方定律的断层流—热耦合数值计算方法[J]. 黄金科学技术, 2020, 28(6): 846-858 doi:10.11872/j.issn.1005-2518.2020.06.122

CHEN Gang, MA Ling, GONG Hongsheng. Numerical Calculation Method of Fault Flow-Thermal Coupling Based on Cubic Law[J]. Gold Science and Technology, 2020, 28(6): 846-858 doi:10.11872/j.issn.1005-2518.2020.06.122

成矿动力学数值模拟方法是了解复杂成矿过程的一种有效方法，近年来诸多学者对该方法进行了研究^［1-5］。利用数值模拟方法，针对成矿动力学过程的一个或几个子过程（如流体流动、热传输、变形和化学反应等）建立相应的数学物理模型时，首先要根据地质和地球化学特征设置相应的初始和边界条件，然后利用计算机求得相应的数学物理方程的数值解，直观地再现成矿系统时空演化细节，从而提高对成矿物理化学过程的认识和找矿预测能力^［4］。成矿动力学模拟涉及矿床学、地球化学、岩石水力学、流体动力学、岩石力学和计算机数值计算等多个学科^［6-7］。金属矿床的成矿过程大多与地质流体有关，可概括为金属的萃取、运移和沉淀3个过程^［8］。这3个成矿过程包括物理和化学过程，但学者们对于化学过程的研究程度高于物理过程^［7］。综合考虑构造、流体和传热等多因素的成矿模拟研究十分必要，有助于更深入地认识热液成矿过程中构造—流体—传热耦合作用^［9］。对于成矿构造而言^［10-11］，断层对成矿模式的控制和影响非常重要。

热液型矿床的物理过程非常复杂，构造往往是其主要的热液渗流通道，成矿动力学数值模拟通常会面临以下问题：（1）岩体中裂隙、断层的渗透性远大于岩石基质，但裂隙、断层数量众多，二者的宽度远小于其延伸方向尺寸，造成模型建立困难；（2）从渗流角度而言，岩石基质中的渗流符合达西定律，而裂隙中的渗流多数情况下符合立方定律，二者相互作用、相互影响形成复杂的流—热耦合问题；（3）流体的动力粘滞性具有随温度和压力变化的特点，粘滞性的变化对模拟成矿结果的影响问题。

岩石节理、岩体裂隙和断层具有尺度联系性和自相似性，可认为是不同尺度裂隙的表现形式^［12］。岩石中的裂隙宽度相对于其延伸方向的尺寸而言要小得多，Juanes等^［13］提出用低维单元来模拟裂隙以减小计算量。在不考虑断层内部结构的情况下，研究断层的渗流和热传导等问题时，断层可以抽象为低维的三维曲面或二维曲线，并借助于裂隙渗流理论进行研究和分析。

岩石裂隙中的渗流不同于孔隙介质的渗流，其具有特定的渗流规律，形成一门独立的学科，即岩石水力学^［14］。裂隙广泛分布于各类地层中，并构成主要渗流通道；而裂隙周围的岩石基质的渗透性很小，但占有地层绝大部分体积和孔隙^［15］。在成矿动力学模拟过程中，对热液的渗流模拟是其中非常重要的一个方面，渗流模拟如果忽略裂隙的渗流影响，其模拟结果可能会存在较大的偏差。裂隙中的流体与基岩之间的流量、热量交换计算十分重要，尤其是在流—热耦合计算过程中。热液渗流对于成矿结果的影响显著，这其中就涉及裂隙与基岩之间的流量、热量交换问题。

对于裂隙中的热流耦合问题，Lauwerier^［16］很早就提出了一定假设条件下单条裂隙中的热传导解析解公式；Pruess等^［17］在Lauwerier研究基础上进一步实现了数值方法求解。Cheng等^［18］在合理忽略裂隙内热存储项的前提下，考虑了基岩中二维热传导得到解析解，说明仅考虑基岩中垂直于裂隙方向的热传导会造成较大的误差。自20世纪 90 年代起，国内外学者初步建立了岩体中热—流—固耦合问题的控制方程^［19-20］。赵阳升等^［20］提出了由基质岩块和裂缝组成的块裂介质模型，通过基质岩块与裂缝之间相互作用的均衡关系建立了固—流—热耦合数学模型，对高温岩体地热资源进行了模拟与评价；张树光等^［21］基于裂隙岩体的流—热耦合数学模型，认为岩体内裂隙水流所引发的热量迁移对裂隙岩体的温度场分布具有重要影响。近年来众多学者^［21-24］利用成熟的数值模拟软件展开了裂隙、断层中流—热耦合和热—流—固耦合等方面的研究工作，并取得了积极成效。

本研究的目的是验证断层概化为空间曲面方法的可行性和合理性，并解决断层中的裂隙流与基岩中的达西流之间的流—热耦合问题。初步探讨流体粘滞性随温度变化这一性质对数值模型计算结果的影响。解决上述问题，可以降低建模难度，在数值模型中增加大量对热液渗流有影响的断层或裂隙成为可能，使热液渗流数值计算更为精细、准确，并扩展数值方法在成矿动力学数值模拟中的应用范围。

1 裂隙渗流理论简介

岩石裂隙形态复杂，必须进行适当的抽象和简化，最简单的方式是简化成由2块光滑平行板构成的缝隙。前苏联学者最早提出了光滑平板流^［14］，随后David^［25］和Louis^［26］对此展开了大量研究，发现通过光滑平板构成的裂隙流量与裂隙宽度的3次方成正比，由此提出了著名的立方定理^［25-26］。之后很多学者展开了更进一步的研究，并取得了显著成果^{［14，27-28］}。

2个光滑平板之间的隙宽为常数 $a$ ，则通过等宽缝隙的流量q与隙宽 $a$ 之间的关系^［14］表示为

q = \frac{g a^{3}}{12 ν} J

（1）

式中：q为流量（m²/s）；g为重力加速度（m/s²）；a为缝隙的隙宽（m）；J为水力梯度；ν为流体的运动粘滞系数（m²/s）。

从式（1）可以看出，通过缝隙的流量大小q与缝隙隙宽a的三次方成正比，即著名的立方定律。为更好地表达裂隙本身的渗透特性，使用渗透率概念κ（m²），其表达式为

κ = \frac{a^{2}}{12}

（2）

渗透率κ表示的是缝隙内在的水力传导性质，只与缝隙隙宽的平方成正比。由于它的概念只关注于缝隙本身，纯粹不受流体因素影响，所以渗透率概念在国内外文献和研究中使用更为广泛。根据前人研究，当缝隙中流体的雷诺数Re<500时，立方定律适用。

裂隙的凸起程度与隙宽的比值对裂隙的渗透性有很大的影响，基于此认识提出了粗糙度的概念，用于修改立方定律^［29］。计算公式如下：

κ = \frac{a^{2}}{12 f}

（3）

式中：f为裂隙表面的粗糙度；κ为裂隙渗透率。

受限于篇幅，本文对于裂隙渗流主要考虑开度和粗糙性2个方面。

2 达西流与裂隙流耦合控制方程

由于裂隙隙宽与岩体尺度之间存在巨大差异，在三维岩体模型中需要将裂隙简化为二维面，在二维岩体模型中将裂隙简化为一维线。此时，需要解决裂隙流与达西流的耦合问题，即不同维度之间相邻网格对象的质量交换、压力传导、速度连续性和能量交换等物理场相互联系问题，否则会出现裂隙与岩体网格内的物理过程各自独立、互不干扰的错误。根据文献^［20-22］整理出三维裂隙中连续性方程和能量方程，并在此基础上进一步分析如下：

三维条件下裂隙内流量方程表示为

q_{F} = - \frac{κ_{F}}{μ} a (\nabla_{T} P_{F} + ρ g \nabla_{T} D)

（4）

裂隙内连续性方程（渗流平衡方程）表示为

a \frac{\partial}{\partial t} (ε_{F} ρ) + \nabla_{T} \cdot (ρ q_{F}) = a Q_{F m}

（5）

式中：q_F表示裂隙中单位长度体积流量（m²/s）；μ为流体动力粘滞系数（Pa·s）；ε_F为裂隙的孔隙度（无量纲）； $\nabla_{T}$ 表示沿裂隙切线方向梯度算子；Q_F_m为质量源项［kg/（m³·s）］； $P$ 为流体压力（Pa）；ρ为流体密度（kg/m³）；κ_F为裂隙的渗透率（m²）， $κ_{F} = \frac{a^{2}}{12 C}$ ；C表示裂隙面粗糙度系数；为了与岩石中的质量源区别，下角标F表示裂隙。

每个裂隙都有2个面，分别称为上、下裂隙面。图1是裂隙流达西流耦合概念模型，可以看出裂隙中流体分别通过上、下裂隙面与基岩产生流量、压力的交换和耦合。由于沿裂隙切向不产生流量，那么从岩石基质中流入裂隙中的流量只能垂直于裂隙面进行，且取决于二者的压差。将裂隙作为岩石基质的边界，且裂隙内部温度和压力连续，可以求出单位长度上从基岩沿裂隙面法向流入裂隙内部的流量为

\begin{array}{l} Q_{F m}^{u p} = - \frac{κ_{F}}{μ} \frac{\partial P^{u p}}{\partial n^{u p}} \\ Q_{F m}^{d o w n} = - \frac{κ_{F}}{μ} \frac{\partial P^{d o w n}}{\partial n^{d o w n}} \end{array}

（6）

式中： $Q_{F m}^{u p}$ 、 $Q_{F m}^{d o w n}$ 表示单位长度裂隙上、下表面从岩石中流入的流量（m²/s）；P^up、P^down分别表示上、下裂隙面压力；n^up、n^down表示裂隙上下表面法向方向；其他符号含义同前文。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 裂隙流—达西流耦合概念模型

Fig.1 Conceptual model of fracture flow-Darcy flow coupling

代入到裂隙的连续性方程式（5）中得出：

a \frac{\partial}{\partial t} (ε_{F} ρ) + \nabla_{T} \cdot (ρ q_{F}) = a Q_{F m}^{u p} + a Q_{F m}^{d o w n}

（7）

由于裂隙隙宽很小，且裂隙内部为连续介质场，所以对一个微小单元格来说，可认为其内部压力一致，表示为

P^up=P^down =P_F （8）

式中： $P_{F}^{u p}$ 、 $P_{F}^{d o w n}$ 分别表示上、下裂隙面压力；P_F表示裂隙内部压力。

上、下裂隙面压力相等，这一特性造成裂隙面两侧的基岩内的压力在靠近同一裂隙时，有压力逐渐相等的趋势。当裂隙流平衡方程为式（7）时，可以实现裂隙流与达西流的耦合。因为是2种物理场的耦合，数值建模时使用不同维度的网格，裂隙网格相较于基质网格而言要低一维度。当裂隙内部压力和岩石基质与裂隙接触面上的压力一致时，可以完成裂隙与基岩、裂隙与裂隙之间的流量交换计算，计算结果可以反映出裂隙对基岩渗流的影响，以及裂隙交叉后的渗流。

裂隙中的能量守恒方程为

a \frac{\partial C_{F} T_{F}}{\partial t} + a ρ_{w} c_{w} {\vec{\vec{u}}}_{w} \nabla_{T} T_{F} - a \nabla_{T} (λ_{F} \nabla_{T} T_{F}) = f^{e}

（9）

式中：a为裂隙的隙宽；C_F表示裂隙内介质的平均体积热容；T_F表示裂隙内的温度；t为时间项；ρ_w表示流体密度；c_w表示裂隙内流体的体积热容；u_w表示裂隙内的流体的渗流速度； $\nabla_{T}$ 表示沿裂隙切线方向梯度算子；λ_F为裂隙介质的热传导系数；f^e表示外源能量的流入。其中：

u_{w} = - \frac{κ_{F}}{μ} (\nabla_{T} P_{F} + ρ g \nabla_{T} D)

（10）

式（9）中等号左边第一项表示能量随时间变化；第二项表示裂隙内流体对流产生的能量增加；第三项表示裂隙内传导方式增加的能量。当裂隙内没有充填物时，C_F和λ_F可取值为裂隙内流体的相应参数。与前文一样，方程等号的右侧能量输入项可以分解为上、下裂隙面2个输入项，表示为

\begin{array}{l} f_{u p}^{e} = ρ_{w} c_{w} (- \frac{κ_{F}}{μ} \frac{\partial P_{u p}}{\partial n_{u p}}) T_{u p} - λ_{e q} \frac{\partial T_{u p}}{\partial n_{u p}} \\ f_{d o w n}^{e} = ρ_{w} c_{w} (- \frac{κ_{F}}{μ} \frac{\partial P_{d o w n}}{\partial n_{d o n w}}) T_{d o w n} - λ_{e q} \frac{\partial T_{d o w n}}{\partial n_{d o w n}} \end{array}

（11）

式中： $f_{u p}^{e}$ 和 $f_{d o w n}^{e}$ 分别表示上、下裂隙面外源能量的流入；λ_eq为等效热传导系数；其他符号含义与前文相同。

由于裂隙的隙宽很小，所以内部温度是相同的，则有：

T_up=T_down=T_F （12）

式（11）代入式（9）后就是裂隙内能量守恒方程。利用式（8）和式（12）作为基质和裂隙的压力、温度边界条件，可以保证裂隙与岩石基质之间的质量和能量连续性；联立式（6）和式（11）方程作用裂隙流量、能量的输入项，可以保证裂隙内的质量、流量守恒，同时保证了达西流和裂隙流2种物理场的耦合。

流体的动力粘滞系数是一个随温度变化的物理量，可表示为μ（T），代入裂隙的渗流平衡方程和能量守恒方程中，可以更真实地反映流体由于温度变化带来渗流影响。成矿热液成分复杂，很难获取各温度下的粘滞系数。为了便于分析，本文使用水代替数值模拟中的热液，并根据文献资料^［30］绘制了水的动力粘滞系数随温度的变化曲线，加入数值模型中进行计算，如图2所示。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 水的粘滞性系数变化曲线^[30]

Fig.2 Variation curve of water viscosity coefficient^[30]

3 数值模型验证及分析

3.1 模型求解方法及验证

流—热耦合是一个复杂的多物理场强耦合过程，本文中的数值模型采用COMSOL软件进行建模和计算，模型中分别使用了达西渗流模块、裂隙流模块和固体传热模块，用于模拟不同的物理过程。为保证模型中裂隙与基质之间的流量、能量的耦合，渗流计算模块中在裂隙上面和下面分别施加式（6）所代表的质量源边界条件，以实现式（7）方程；传热模块中在裂隙上面和下面分别施加式（11）所代表的热能量边界条件。

为了验证本文中热—流耦合公式的合理性和准确性，本文选择2个案例进行分析。第一个案例是一个简单的单裂隙传热模型，边界条件与Lauwerier^［16］一致，对比数值解与解析解。第二个案例为包含断层和3个地层的三维地质模型。

Lauwerier^［16］在研究油层传热问题时，为便于开展研究，提出以下假设条件：

（1）水只在裂隙中进行层流流动，且在裂隙中的水温沿y轴方向保持相同，即裂隙水的温度仅仅与x坐标相关；

（2）裂隙宽度为2b，裂隙内水的流速恒定；

（3）裂隙及上下两层岩石初始温度为T₀，注入水的温度为T_in；

（4）忽略岩石基质中平行于裂隙方向的热传导；

（5）裂隙周围岩石基质为无限厚。

通过分析和推导，给出了特定条件下裂隙水温度随着时间变化的解析解方程；经过时间t后，裂隙内水的温度沿x轴分布表示为

T (x, t) = T_{0} + (T_{i n} - T_{0}) e r f c [\frac{(λ_{x} x) / (ρ_{w} C_{w} d_{f})}{2 \sqrt[]{λ_{s} / (ρ_{s} C_{s}) u_{f} (u_{f} t - x)}}] U (t - \frac{x}{u_{f}})

（13）

式中：erfc为余误差函数；U为单位阶跃函数；u_f为裂隙中流体流速；T₀为初始温度；T_in为入口流体温度；λ_x为裂隙的热传导系数；λ_s为地层的热传导系数；ρ_w为流体密度；C_w为流体热容；C_s为地层热容；t为时间。

第一个案例的数值解与解析解结果对比见图3。从图中可以看出，使用本文的耦合控制方程的数值模型求解结果符合理论解析解，说明本文的耦合控制方程及计算方法是可行的。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 裂隙内部温度变化曲线

Fig.3 Temperature change curve inside the fracture

3.2 案例二数值模型及边界条件

本文建立了一个包含断层的三维地质模型，基岩中为达西渗流，断层采用裂隙流，基岩与断层之间采用前文中的方程进行流—热耦合，根据计算结果观察和分析耦合方程的合理性，并分析断层对流—热耦合的影响。模型包含3个地层和一条陡倾断层（几何模型见图4），每个地层都有不同的物理参数。各地层的属性参数参考已有文献中的数据^{［2， 31］}做适当补充和调整后得到详细属性，如表1所示。模型底部尺寸为500 m×500 m，高度为500 m，地层界面及断层面为三维空间曲面，数值模拟方法采用有限元。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 地质模型结构及测点空间位置

（a）模型结构（b）模型中测点分布位置

Fig.4 Geological model structure and spatial location of measuring points

表1 地层参数列表

Table 1 List of formation parameters

地层	渗透率/m²	孔隙率	导热系数/（W·m^-1·K^-1）	密度/（kg·m^-3）	恒压热容/（J·kg^-1·K^-1）	隙宽/m	裂隙粗糙度
第一层	1.00E-14	0.1	2	2 300	900
第二层	1.00E-10	0.3	3	2 500	850
第三层	1.00E-11	0.3	3.5	2 700	850
断层		1	10	1 000	4 200	0.002	1.6

新窗口打开| 下载CSV

模型初始地下水水头地表为0 m，随深度线性增加。模型四周热力学边界为开放边界，底面为绝热层；顶面为对流热通量边界，以模拟地层与大气之间的热量交换。模型初始环境温度为20 ℃。在底部断层施加随时间变化的热源和压力，模拟高温热液从断层流入，具体见图5温度和压力曲线图。模型中的流体粘滞性设定为随温度变化，采用图2所示曲线。模型采用瞬态求解，总的计算时间为200年。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 模型施加的温度—压力曲线

Fig.5 Temperature and pressure curve applied by the model

3.3 计算结果分析

主要围绕断层通量、传导热通量以及基质中的温度变化对计算结果进行分析。图6（a）是断层的通量曲线，分别展示了断层上方边界通量、下方边界通量、总通量以及上下方边界通量之和随时间变化的情况。需要说明的是，由于上下方边界法向方向相反，所以通量之和实则是上下方边界通量数据之差。从图中可以看出，断层总通量与上下方通量之和保持一致，这表示断层内部的质量守恒，即式（7）是成立的。图6（b）是断层热传导通量曲线图，同样展示了上方边界通量、下方边界通量、总通量以及上下方边界热传导通量之和随时间变化的情况（为了保证图形清晰，仅将80年内的变化反映在图上，后续曲线变化很小）。从图中可以看出，断层总热传导通量与上下方通量之和数量一致，表示断层内部的通量守恒，即式（9）是成立的。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 数值模型断层面通量曲线

Fig.6 Flux curves of fault plane of numerical model

模型只在断层底部施加了温度和压力条件，断层中的裂隙流与基岩中的达西流之间的耦合是否形成，可以观察基岩中的温度是否随断层的温度和压力变化来进行分析。在模型的第二层基岩中，同一水平高程上，取5个与断层面不同距离的空间点［见图4（b），空间点分布在红线上］，观察其温度和压力变化，计算结果如图7所示。从图7（a）中可以看出，距离断层越近，基岩中的温度响应越迅速。由于基岩的热传导性较差，其温度变化具有滞后性，图中温度曲线符合这一情况。观察图7（b）中的压力曲线分布可知，距离断层面越远，同一时间的压力越低；在断层底部施加的压力驱动基岩中的流体向四周压力低的区域流动，并随着距离的增加而降低，符合沿压力梯度方向减小的趋势。此外，压力的传导速率很快，所以各点的压力曲线没有显著的滞后现象。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 第二层基岩中多点温度、压力曲线

Fig.7 Multi-point temperature and pressure curves in the second layer of bedrock

综合图6和图7的曲线分布情况，将上文中的控制方程应用在模型后，可以实现质量和能量的守恒，并保证了数值模型计算域内渗流场压力、速度和质量的连续性，可以说明断层与基岩流—热耦合成立，控制方程是合理的。

为了直接地观察模型中温度场的变化，选择了2个横切断层的剖面，以展示不同时间的温度分布（图8）。从图中可以看出，初期高温区域主要集中在断层及其附近区域。随着断层底部施加的高温高压边界条件的终止，断层内部温度快速下降，高温区域主要分布在基岩中。受断层的热传导作用影响，基岩中输入了热能，但传导速度很慢。基岩中的热能消散主要有3个方面：一是反向向断层进行热传导；二是与地表大气间的热交换；三是向四周开放边界进行热传导。从图中可以看出，基岩向四周开放边界的热传导速度最慢，历经100多年的计算时长，仍不能完全恢复到环境初始温度20 ℃。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 不同时间模型内温度分布切片图（温度单位：℃）

Fig.8 Slices of temperature distribution in different time models（temperature unit：℃）

将模型中的断层抽象为三维曲面裂隙，其隙宽参数为2 mm，考虑裂隙粗糙度因素，其水力隙宽换算后约为1.58 mm。这一宽度数值对自然界中的断层而言也不算大，但其对热能、渗流的传导影响却造成较大范围内基岩流—热产生巨大的变化，从侧面反映出断层流—热的影响巨大，尤其对于热液成矿的矿床，断层影响是不能忽略的重要因素。

为了分析模拟过程中流体粘滞性变化与否对模拟结果的影响，复制一份上述的数值模型，只修改流体粘滞性系数为不随温度变化的固定值0.001 Pa·s（即20 ℃时水的粘滞性系数）并进行计算，将计算结果与前一模型结果进行对比。为便于直观地反映2个模型计算结果的差异，取2个模型中断层的通量、热传导通量的比值的绝对值（ $|粘滞系数随$ 温度变化的模型结果/粘滞系数固定值模型结 $果|$ ）绘制成随时间变化的曲线，如图9所示。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 断层通量、热传导通量的比值曲线

Fig.9 Ratio curves of flux and heat conduction flux

从图9可以看出，当流体粘滞系数随温度变化时，对比粘滞系数为固定值的模型，其通量计算结果最大有6倍的差距，而热传导通量可达到最大60倍的差距，即在质量和能量传输过程中，由于温度的升高造成流体粘滞性的降低，可极大地加强其传输能力。更进一步说明，在热液型矿床成矿动力学模拟时，流体粘滞性变化这一因素是不可忽略的因素之一，否则将会带来很大的偏差。

4 结论

成矿动力学模拟是了解成矿过程的有效方法之一，成矿过程由于其复杂性涉及多个学科领域。本文以裂隙渗流理论为基础，探讨了断层在成矿动力学模拟过程中的影响，综合分析得出如下结论：

（1）裂隙和岩石基岩中的渗流有着不同的渗流理论，本文得出的裂隙流和达西流耦合控制方程，保证了数值模拟过程中需要保证数值模型计算域内渗流场压力、速度、质量和能量的连续性；通过对数值模型的分析，认为耦合控制方程可以合理、准确地实现裂隙流和达西流的耦合。

（2）由于裂隙和断层的空间形态具有宽度远小于延伸尺寸的特性，本文将断层简化为曲面，基于热—流耦合控制方程，使用数值方法计算了涵盖达西渗流、裂隙流和热传导的成矿动力学过程的数值模型；由于裂隙流内渗流速度快、热传导效率高，对模型最终的温度场和能量传导影响十分显著。

（3）经过数值模型的验证认为，在不考虑断层内部结构时，成矿动力学模拟中将断层概化为无几何厚度的空间曲面，使用裂隙渗流理论进行计算是可行的。

http://www.goldsci.ac.cn/article/2020/1005-2518/1005-2518-2020-28-6-846.shtml

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Liu

L M

，Zhao

Y L

，Zhao

C B

Coupled geodynamics in the formation of Cu skarn deposits in the Tongling-Anqing district，China：Computational modeling and implications for exploration

［J］.Journal of Geochemical Exploration，2010，106（1/2/3）：146-155.