20世纪60年代以来，地质热流体的定量化研究沿着两个方向发展，其一为无反应条件下的传热和流体动力学研究（Donaldson，1962；Elder，1967a，1967b；Combarnous and Bories，1975；Ribando and Torrance，1976；Norton and Knight，1977；Norton and Knapp，1977；Norton and Taylor，1979；Cathles，1981；Norton，1984；Garven，et al.，1984a、1984b；Bethke，1985，1986，1989；Bethke，et al.，1988；Garven，1989；Garven，et al.，1993；Phillips，1990；Lowell，et al.，1993；Raffensperger and Garven，1995a、1995b；陈跃庭，1989；任启江等1994a、1994b；郭国章等，1994，1997）；其二为无流动条件下的多相多组分化学反应研究（Helgeson，1968；Helgeson，et al.，1969；Brimhall，1980；Reed，1982；Brimhall and Ghiorso，1983；Sverjensky，1984；Bowers and Taylor，1985；徐士进，1985；王益锋等，1987；周会群等，1988），这两方面研究都取得了大量成果。80年代，现代地质学者开始定量研究化学反应和物质输运的耦合过程（Rubin，1983；Miller and Benson，1983；Walsh，et al.，1984；Lasaga，1984；Lichtner，1985，1988，1992，1993；Appelo，et al.，1987；Bryant，et al.，1987；Kirkner，et al.，1988；Ague，et al.，1989；Liu，et al.，1989a、1989b；Novak，et al.，1989；Yeh，et al.，1989；Wells，et al.，1991；Phillips，1990；Lichtner，et al.，1992；Steefel，et al.，1992；Sevougian，et al.，1993），这些研究绝大多数建立在平衡近似基础上。
在众多传热和流体动力学研究中，地质学者较早、较多注意的是与岩浆活动有关的热液体系。成矿作用过程中，热量、流量、质量迁移、输运定量模型研究以斑岩型矿床居多（Norton，et al.，1977；Norton，1982；Cathles，1977；Johnson，et al.，1985；Nesbitt，1990；任启江等1994a、1994b；郭国章等，1994，1997）。1977年，美国学者Norton等（1977a、1977b）首先发表了深成环境下热和质量传输的计算模拟成果。同年，Cathles（1977）运用有限差分法计算模拟了火成岩侵入与液体和蒸汽为主导的地热体系的关系，并依次探讨了斑岩铜矿的成因。随后几年间，Norton等学者（Henley，et al.，1983；Knapp，et al.，1981；Norton，1982）又针对亚利桑那东南部一些特定斑岩铜矿区用有限元法模拟了可渗岩体周围热液运移途径和质量传输过程，强调了岩体及其围岩渗透率大小对体系的含矿性有重要影响。90年代开始，国内任启江教授等人（任启江等，1994a、1994b；郭国章等，1994，1997）对陕西金堆城斑岩钼矿、安徽沙溪斑岩铜（金）矿和江西德兴斑岩铜矿等矿床的成矿过程也开展了类似的研究，并发表了相应成果。
与大气降水来源流体及盆地流体大规模运移成矿有关的定量模型研究开展较晚（Garven，1995），而深部和浅部地下水流在矿床形成过程中起着重要作用这一思想已持续约一个世纪。1882年，Chamberlin已识别出地层中金属的侵染状特征并表明Wisconsin西南部碳酸盐岩中的矿床是由下渗的地下水流所形成（Garven，et al.，1993）。Van Hise and Bain（1902），Cox（1911）和Siebenthal（1915）在密西西比河谷和三州（Tri-State）铅锌矿床研究中也提出了相似的概念模型。近年来，区域地下水流体系、盆地构造和层控矿床成矿作用之间的成因联系已经得到较好认识。如，赋存于碳酸盐岩的大型密西西比河谷型方铅矿、闪锌矿、重晶石和萤石矿被认为是经过盆地流体的区域规模运移而形成（Leach，1979；Anderson and Macqueen，1982）。沉积盆地中区域流体运移的水文模拟表明，前陆盆地的构造抬升所引起的地形（或重力）驱动水流体系能够形成大型的层控矿床（Garven，et al.，1984a、1984b）。类似的水文体系在一些世界超大型油田的石油聚集和圈闭过程中也存在（Garven，1989；Bethke，et al.，1990）。深部地下水流作为一种成矿营力已得到地质学家的认可。同时，地球化学、地热、地球物理方面也不断地提出新证据（Hitchon，1984；Oliver，1986；Behke，et al.，1990；Deming，1994；Jessop，et al.，1994）。本文重点阐述区域古流体成矿作用模拟研究中所涉及的传热、流体运移基本控制方程及其数学解法。
一、二维传热、流体运移基本控制方程
1.流体流动方程
在流体—岩石体系中，密度为ρ的单一相稳态水流体系的质量守恒方程为（Garven et al，1984a；Garven，1985，1986；Gvirtzman，et al，1997a，1997b）：
湘中区域古流体及锡矿山锑矿成矿作用模拟
比流量qi可由达西定律求得（Garven，et al，1984a、1984b；Garven，1985，1986，1989；Garven，et al，1993，1994；Forster and Smith，1990）：
湘中区域古流体及锡矿山锑矿成矿作用模拟
式中：kij——内在渗透率张量或称渗透率张量（m2）［L2］，通常采用的单位为cm2或da
（darcy），1da=9.8697×10-9cm2；
u——流体的动力粘滞系数（Pa·s或kgm-1s-1）［ML-1T-1］；
P——流体的压强（Nm-2或kg·m-1·s-2）［ML-1T-2］；
g——重力加速度（ms-2）［LT-2］；
qi——达西渗透速度或称比流量（md-1或ma-1）［LT-1］；
x3——参考点以上的标高（坐标）（m）［L］；
xj——笛卡尔坐标（m）［L］；
i，j——1，2，3 并遵循求和约定。
层控矿床的形成依赖于非等温含水层中变密度流体（盆地卤水）的运移（Garven，et al，1984a，1993）。可定义相对粘滞系数为ur，相对流体密度为ρr，其分别为
湘中区域古流体及锡矿山锑矿成矿作用模拟
湘中区域古流体及锡矿山锑矿成矿作用模拟
上式中的u0和ρ0分别为参考粘滞系数和参考密度。依据等效淡水水头（实际水头）的含义有
湘中区域古流体及锡矿山锑矿成矿作用模拟
变换上式得：
湘中区域古流体及锡矿山锑矿成矿作用模拟
把（5-6）式代入（5-2）式并化简得：
湘中区域古流体及锡矿山锑矿成矿作用模拟
式中：
———多孔介质的渗透系数（水力传导系数）张量（md-1）[LT-1]
———重力方向单位矢量第j个分量，即e1=e2=0，e3=1。
将（5-7）式引入（5-1）式有：
湘中区域古流体及锡矿山锑矿成矿作用模拟
对于剖面二维稳态水流体系（i，j=1，3即x，z）当坐标轴与渗透系数张量的主方向平行时，（5-8）式可化简为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
由（5-9）式结合适当的边界条件，以及模型介质参数可求得渗流区内（计算域）各点的水头h（x，z），并由（5-7）式可进一步求得各点的比流量qi（渗透速度）。
实际平均线速度公式（Freeze and Cherry，1979；Garven，et al.，1984a）为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
式中：φ—多空介质的孔隙度。
据（5-10）式可求出流域内的实际平均流速（线速度）vi。实际平均流速用于预测质量迁移中溶质的运移；渗流速度qi用于计算热迁移中对流项的作用。
2.热流量方程（能量方程）
传导、对流、弥散皆可引起地下热量迁移，从而影响地温场的分布特征。对于沉积盆地中区域热流问题，当水流体系为稳定态以及盆地热流值变化不显著时，稳态分析足以解决温度场的分布特征问题（Garven，et al.，1984a、1984b）。
对于一流体-岩石体系，在稳态条件下，饱水多孔介质的能量守恒方程为（Garven，et al.，1984a；Garven，1989）：
湘中区域古流体及锡矿山锑矿成矿作用模拟
式中：T——温度（℃）［θ］；
cf——流体的比热（J℃-1kg-1）［L2T-2θ-1］；
λij——多孔介质的有效热传导系数（Jm-1d-1℃-1）［MLT-3θ-1］；
Dij——多孔介质的热机械弥散系数（Jm-1d-1℃-1）［MLT-3θ-1］。
（5-11）式中左端第一项代表热传导作用；第二项代表热机械弥散作用；第三项代表热对流作用。对于小规模热迁移问题，相对于热传导而言，热机械弥散较小；而在区域规模上，热机械弥散对温度场的影响变得较为显著（Mercer and Faulst，1980）。数学分析表明，在自然对流热迁移问题中（Tyvand，1981），弥散作用增加了流场的稳定性，并且抑制对流的发生。
定义一个热传导-弥散张量Eij如下：
湘中区域古流体及锡矿山锑矿成矿作用模拟
则（5-11）式化简为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
结合定解条件，由上式就可求得温度场的空间分布特征。
3.状态方程
在沉积盆地中流体密度及动力粘滞系数的变化主要受温度和盐度影响，而压力对其影响甚微。为了模拟盆地中流体的运移过程，必须给定流体密度、动力粘滞系数与各影响因素间的函数关系。
对于液态水而言，当温度低于250℃时，其比热和热传导率可以为常数，而其密度和粘度仅是温度的函数，其表达式为（Forster and Smith，1990）：
湘中区域古流体及锡矿山锑矿成矿作用模拟
在250～400℃，温度、压力均为水特性的影响因素，但温度占主导作用；当温度超过400℃时，压力则成为主要因素。
前面已经提到，对于盆地（地壳）内大规模运移流体而言，其流体密度、粘度的变化均受到流体盐度的影响。Kestin（1978，1981）等给出了20～150℃，0.0～6.0mol（NaCl），0.1～35.0MPa范围内，流体密度动力粘度的变化关系。这种变化关系可用于流体的大规模运移模拟。
二、数值计算方法
地下水流问题有三种解法：即解析法、数值法、模拟法。前面所给出的方程是一组庞大而复杂的非线性代数—偏微分方程组，目前只能借助于数值法求解，并且需要做进一步化简。数值法包含有限差分法和有限单元法（简称有限元法）。
有限差分法是一种古典的数值计算方法。其基本思路是：用渗流区内有限个离散点的集合代替连续的渗流区，在这些离散点上用差商近似地代替微商，将微分方程及其定解条件化为以未知函数在离散点上的近似值为未知量的代数方程（称之为差分方程），然后求解差分方程，从而得到微分方程的解在离散点上的近似值。
鉴于有限差分存在以下缺陷：①难以处理复杂的区域和边界；②只考虑结点的作用而忽略了结点间的单元，因而在某些情况下计算精度远远不如有限元法，故本文采用有限元法进行数值求解。
1.流体运移方程的三角形单元Galerkin有限元法
设所研究的渗流区域D被剖分为M个三角形单元，结点总数为N。其中，设内结点和第二类边界结点共n个，并将它们编号为1，2，…，n，第一类边界结点编号为n+1，n+2，…，N。
用Galerkin法求解二维稳态区域流体运移问题，就是求下列形式的试探函数（陈崇希等，1990；Garven，1989）：
湘中区域古流体及锡矿山锑矿成矿作用模拟
作为方程（5-9）的近似解，并使其满足一定的边界条件。式中：NL（x，z）（L=1，2，…，N）是N个线性无关的函数组，称为基函数，HL是结点L处的水头值。由于h是微分方程的近似解，因此，一般说来，将（5-16）式带入（5-9）式时：
湘中区域古流体及锡矿山锑矿成矿作用模拟
称R（x，z）为误差函数。
在实际计算中，总希望R（x，z）在区域D上的加权积分等于零。即：
湘中区域古流体及锡矿山锑矿成矿作用模拟
式中ωL为权函数。
Galerkin有限元法是一种取基函数为权函数的加权剩余法，即：
湘中区域古流体及锡矿山锑矿成矿作用模拟
应用多边量乘积求导公式，（5-17）式可展开为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
湘中区域古流体及锡矿山锑矿成矿作用模拟
将（5-20）式代入（5-19）式得：
湘中区域古流体及锡矿山锑矿成矿作用模拟
在任何一个三角形单元e上，粘度u、密度ρ和渗透率K均为常数（陈崇希、唐仲华，1990），则（5-21）式可化为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
式中：he——为单元e上的水头近似函数；
———为单元e上结点L的基函数。
下面将会看到，在所研究的区域上，基函数NL（x，z）通常是分段线性地但又连续定义的，而NL的一阶导数可以不连续。由于h（x，z）是NL的线性组合，其一阶导数也是不连续的，从而使得（5-23）式中积分的计算变得困难。这里采用下述方法对被积函数中的二阶导数降阶。
由于
湘中区域古流体及锡矿山锑矿成矿作用模拟
将上述两式相加，在单元e上积分并移项得：
湘中区域古流体及锡矿山锑矿成矿作用模拟
湘中区域古流体及锡矿山锑矿成矿作用模拟
应用Green公式：
湘中区域古流体及锡矿山锑矿成矿作用模拟
（5-25）式代入（5-24）式得：
湘中区域古流体及锡矿山锑矿成矿作用模拟
式中：Be是单元e的边界，曲线积分沿Be的正向积分。（n，x），（n，z）分别表示Be的外法线方向n分别与x轴正向和z轴正向的夹角。
将（5-26）式代入（5-23）式得：
湘中区域古流体及锡矿山锑矿成矿作用模拟
这就是二维稳态区域流体运移的Galerkin方程。
在单元e上，其结点按逆时针顺序编号，依次为i，j，k，规定与三个结点相联系的基函数值为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
其中：
湘中区域古流体及锡矿山锑矿成矿作用模拟
Δ为三角形单元e的面积：
湘中区域古流体及锡矿山锑矿成矿作用模拟
基函数在结点I（I=i，j，k）上的值为1。在其他结点及与结点I相邻的单元以外其值为0，因而在单元e上水头的近似函数可表示为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
在二维稳态区域流体运移的Galerkin方程（5-27）式中，只有n个待求水头值H1，H2，…，Hn（即内结点和第二类边界结点的水头）。因而，求解该方程只需要建立n个方程即可。为此，在Galerkin方程（5-27））式中分别取L=1，2，…，n。其中的均为内结点和第二类边界结点的基函数，且在第一类边界B1上取值为零。根据基函数性质，NL（x，z）只在结点L的子区DL内不为零，而在DL外全取零值。
设子区DL内有mL个单元，则方程（5-27）式为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
其中Be表示单元的边界，B2表示第二类边界，B2∩Be表示B2与Be的公共部分。
依据（5-29）式得：
湘中区域古流体及锡矿山锑矿成矿作用模拟
将（5-31）式和（5-32）式代入（5-30）式：
湘中区域古流体及锡矿山锑矿成矿作用模拟
按Hi、Hj、Hk归并，得：
湘中区域古流体及锡矿山锑矿成矿作用模拟
湘中区域古流体及锡矿山锑矿成矿作用模拟
若记第L个方程中HP的系数为ALp（L=1，2，…，N），则（5-33）式变成
湘中区域古流体及锡矿山锑矿成矿作用模拟
将上式中与第一类边界有关的项（ALPHP，P=n+1，n+2，…，N）移到右边，并与合并，记作QL，用矩阵形式表示为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
2.能量方程的GaLerkin有限元方法
对于能量方程（5-2-13），应用Galerkin法求解，同样是寻找如下形式的试探函数
湘中区域古流体及锡矿山锑矿成矿作用模拟
其分析过程如同上述的流体运移方程Galerkin法求解，此处不再赘述。结果可用矩阵形式表示为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
其中：
湘中区域古流体及锡矿山锑矿成矿作用模拟
单元e内的达西流速为：
湘中区域古流体及锡矿山锑矿成矿作用模拟
式中：j=1，2，3对应单元的三个结点。
三、计算过程及程序编制
上面应用Galerkin有限元法分别对流体运移方程、能量方程（热量迁移方程）进行了离散，最后得到线性代数方程组（5-35），（5-36），即：
湘中区域古流体及锡矿山锑矿成矿作用模拟
对于方程（5-39）式和（5-40）式求解有三种方法：①微分代数法（DAE方法）（Lichtner，1985）；②直接代入法（DSA方法）（Lewis，et al.，1987）；③顺序迭代法（SIA方法）（Garven，et al.，1984a、1984b；1989；Cederberg，1985；Garven，et al.，1993）。目前广泛采用的是SIA方法。
方程（5-39）和（5-40）中的系数矩阵[A]和[R]是对称、正定的稀疏矩阵，采用SIA方法及一维变带宽压缩存储技术可以将计算内存减少到最小。具体的求解上述方程组的步骤如下：
（1）参数输入（包括水动力学参数、热参数）、流场及温度的边界条件。
（2）模拟计算开始，假定不存在盐度、温度梯度，计算各结点水头值压力值。
（3）根据（5-37）和（5-38）式计算网格中每个单元的达西流速。
（4）求解温度方程（5-36）式，得到各结点的温度值。
（5）由压力、温度及特定的盐度新给值，根据状态方程计算流体的密度和粘度。重复（2）～（5）直到迭代收敛于稳定的温度值。当最大温度变化介于网格所有结点的特定公差时，迭代达到收敛，计算完成。
图5-2 计算流程简图
计算流程图如图5-2。
内容全部来源于网络收集，如有侵权，请联系网站删除：QQ：24596024