数值模拟显式计算方法

‘壹’ 数值模拟模型

6.2.2.1 地下水系统的模拟模型

潜水含水层的模拟模型：

含有协变量的地下水动态规划管理模型研究

式中：D——地下水系统的模拟渗流区域；

h——地下水位（m）；

M——含水层厚度（m）；

K——渗透系数（m/d）；

T——导水系数（m²/d）；

t——时间（d）；

μ——给水度；

μ^*——承压含水层贮水系数；

Γ₁，Γ₂——一类及二类边界；

n→——边界的外法线方向；

h₁（x，y，t）——第一类边界的地下水位（m）；

h₀（x，y）——初始地下水位（m）；

q（x，y，t）——第二类边界的单宽流量（m²/d）；

P——地下水开采量（m³/d）；

ε——天然条件的源汇项（m³/d），如降水入渗补给、蒸发排泄等；

Q——与地下水位有关的源汇项流量（m³/d），即协变量。

6.2.2.2 模拟软件选取

本次工作采用Visual MODFLOW软件对水流进行模拟，该软件具有处理含有协变量的地下水模拟模型的功能。

6.2.2.2.1 河水与地下水的交换量

Visual MODFLOW河流程序包输入文件要求河流边界条件的每个网格单元具有以下资料：

（1）河流水位：地表水体的自由水面标高，这个标高可以随时间变化。

（2）河床底板：地表水体的渗流层（夹层）的层底板标高。

（3）水力传导系数：一个表示由渗流层（河床）引起的地表水体和地下水间水流阻力的数值参数。

水力传导系数（C）可以通过单元的长度（L）、单元中河流的宽度（W）、河床厚度（M）以及河床物质的渗透系数（K），利用下面公式计算得出：

含有协变量的地下水动态规划管理模型研究

6.2.2.2.2 排水沟的排泄量

Visual MODFLOW支持MODFLOW包括的标准排水沟边界程序包。排水量正比于含水层的水头与某一固定水头或高程之差。如果含水层水头低于排水沟的固定水头，排水沟程序包假设排水就没有效果。

排水沟程序包要求包含这个边界条件的每个单元具备下列输入资料：

（1）水头标高：排水沟标高或者排水沟内的自由水面的排水沟水头。假设排水沟并没有充满水，这样排水沟管内的水头接近等于排水沟管的中间标高。

（2）水力传导系数：排水沟水力传导系数是一个综合系数，描述了排水沟和地下水系统之间的水头损失。计算排水沟的水力传导系数没有通用的公式，大多数情况下，计算排水沟水力传导系数所需的详细资料对于地下水模型来说无法得到。这些细节包括排水沟周围水头分布的详尽资料、排水沟附近含水层的渗透系数、填充材料的分配、排水管的数目和尺寸大小、堵塞物质的总量和堵塞物质的渗透系数。通常利用水流率测量值和水头差来计算排水沟传导系数。排水沟水力传导系数在模型校正时通常会做调整。

6.2.2.2.3 蒸发量

蒸发蒸腾程序包模拟植物蒸腾、直接蒸发以及由地下水饱水状态移走水而引起的地表面渗漏作用。

蒸发蒸腾程序包需要以下信息：

（1）蒸发蒸腾速度：当水面高程与网格单元高程一致时的蒸发蒸腾速度。

（2）极限深度：网格单元顶板以下的深度，在此深度蒸发蒸腾速度可以忽略不计。

蒸发蒸腾程序包方法基于以下假设：

（1）当水面位于或者高于地表面（第一层顶板）时，水面蒸发蒸腾损失以用户定义的最大量进行。

（2）当水面高程低于极限深度时，或者位于第一层以下时，水面蒸发蒸腾可以忽略不计。

（3）两个极限之间，水面蒸发蒸腾随水位高程呈线性变化。

6.2.2.2.4 泉流量

Visual MODFLOW中没有单独计算泉流量程序包，泉流量在识别、检验和预报阶段均采用排水沟计算模块计算。

应用排水沟程序包计算泉流量时要求输入下列资料：

（1）水头标高：输入泉出露点地表高程。

（2）水力传导系数：输入经识别、检验后与排水沟水力传导系数等效的泉流量的比例系数。

在本次研究过程中，对松原地区地下水系统中的互馈协变关系问题作了一些合理的概化，重点考虑了泉流量、蒸发量、地下水与河水的交换量这三个协变量，并采用Visual MODFLOW软件相应的程序包计算。

6.2.2.3 研究区网格剖分

对松原地区水文地质图扫描，并应用ArcGIS软件矢量化，作为计算区的剖分底图。计算模拟区应用Visual MODFLOW软件共剖分出80行、100列，4071×3个矩形单元体（图6.6）。

图6.6 研究区网格剖分图

6.2.2.4 模型的校正

建立一个描述实际地下水系统的数值模拟模型，实际上就是找一个描述它的偏微分方程并确定其定解条件。一个正确可靠的数学模拟模型应当是实际地下水系统的复制品。也就是说，当施加天然或人为影响时，模拟模型的反应与实际地下水系统的反应应当非常接近，但在实际工作中很难直接达到这一点。首先，选用的方程的类型不一定合适。其次，代入方程中的各种参数，不论从现场还是在实验室都难以准确的获得，边界条件也往往缺乏可靠的资料。此外，在建立水文地质概念模型的过程中，对地下水系统的实际条件又做了许多的简化和假设。所以上述这些因素都可以导致初步建立的模拟模型与实际问题有很大的差距。因此，我们必须根据抽水试验和已有的地下水系统长期动态观测资料对初步建立的模型进行校正。

模型的校正包括两个过程，即识别和检验，通常运用已有的不同时段的两套地下水数据，一套用于识别，另一套用于检验。然后，把数值模拟模型计算出的结果与已有的实际数据进行对比，看两者是否一致。若不一致，就要对模型继续校正，直到满意为止。

6.2.2.4.1 模型的识别

（1）模型识别时段的选择：本次模拟的识别时段选取1999年12月1日到2000年2月29日，共计91d。本次模拟选取枯水期，因为期间源汇项少，地下水位比较稳定，较容易反映含水层的水文地质特性。

（2）参数分区与初值的确定：根据区域水文地质情况和抽水试验成果，潜水含水层水文地质参数分为9个区，弱透水层水文地质参数分为4个区，承压含水层水文地质参数分为8个区。参数分区见图6.7～图6.9。

图6.7 潜水含水层参数分区图

图6.8 弱透水层参数分区图

图6.9 承压含水层参数分区图

（3）源汇项的计算：模型的识别时段正处于冬季，所以模拟区内无大气降水入渗补给、灌溉回渗补给、农业开采和蒸发排泄等，排泄方式主要有开采排泄和泉的排泄，源汇项的计算较简单。

（4）模型识别时段拟合结果：从水位拟合图（图6.10～图6.11）可以看出，潜水观测井水位误差的平均绝对值为0.295m，拟合误差小于0.5m的观测井占总观测井的88.73%；承压水观测井水位误差的平均绝对值为0.302m，拟合误差小于0.5m的观测井占总观测井的86.48%，可见模型识别取得了较好的结果。

表6.3—表6.5为潜水含水层、弱透水层和承压含水层水文地质参数识别结果。

图6.10 潜水含水层识别时段水位拟合图

图6.11 承压含水层识别时段水位拟合图

表6.3 潜水含水层水文地质参数识别结果表

表6.4 弱透水层水文地质参数识别结果表

表6.5 承压含水层水文地质参数识别结果表

6.2.2.4.2 模型的检验

为了进一步验证所建立的数值模拟模型和确定的水文地质参数的真实性，再次利用已有的地下水位动态观测资料对模型进行检验。

（1）检验时段的选择：选择2005年6月1日至2005年8月31日，共计92d作为模型的检验阶段。潜水含水层和承压水含水层分别选取有长期观测资料的观测井进行模型的拟合检验。

（2）源汇项的计算：模型的检验时段正处夏季丰水期，源汇项比较复杂。主要补给来源有大气降水入渗补给、灌溉入渗补给以及侧向径流补给等。主要排泄方式有人工开采排泄、潜水蒸发排泄以及泉的排泄。

（3）模型的检验时段拟合结果：从模型的检验时段拟合结果可以看出（图6.12，图6.13），潜水含水层水位拟合误差小于0.5m的观测井占总观测井的89.53%；承压含水层水位拟合误差小于0.5m的观测井占总观测井的85.28%。模型的检验结果表明，建立的地下水数值模拟模型和确定的水文地质参数是真实可靠的。

6.2.2.5 模型的预报

经过识别和检验的模型，就相当于实际地下水系统的复制品，就可以以模型的行为代表实际地下水系统的行为了。对于各种自然因素、人为因素的影响（比如各种规划开采方案），都可以通过数学模拟模型的运转来预报其后果，这也是建立模拟模型的主要目的之所在。

6.2.2.5.1 预测期、初始条件的确定以及边界条件的处理

本次预测期限为4年，即2008～2012年。根据已有的地下水位动态观测资料，选择地下水位预报的初始时刻为2008年9月1日。

图6.12 潜水含水层检验时段水位拟合图

图6.13 承压含水层检验时段水位拟合图

由于开采水源地距离计算模拟区边界较远，并且考虑到开采时对边界的影响较小，因此本次采用初始时刻的边界条件进行计算。

6.2.2.5.2 源汇项的预报

（1）降水量和水面蒸发量的预报：本次应用蒙特卡罗（Monte-Carlo）方法对研究区内的降水量和水面蒸发量进行随机模拟和预报，这是一种应用随机数来进行模拟试验的方法。该方法把降水量和蒸发量作为一个随机变量来处理，根据过去的资料提取研究区降水量和蒸发量所具有的统计规律性，这些统计规律性代表了该地区降水量和蒸发量的本质特征。

（2）灌溉入渗补给量的预报：灌溉入渗补给量根据研究区内的农业用水额度来进行预报。

（3）开采量的预报：由于经济发展和人口增长等原因，松原地区用水量在逐年增加。2002年地下水开采量为4.×10⁸m³ ，2008年地下水开采量为5.2201×10⁸m³ ，平均年增长率定为3.83%。按该增长率对2012年开采量进行预测，即6.067×10⁸m³。6.2.2.5.3 预报结果及分析

根据地下水开采量平均年增长率3.83%对开采量进行预测，即计算模拟区2012年地下水开采量为6.067×10⁸m³。从预报结果的水位等值线图（图6.14和图6.15）可以看出，潜水水位和承压水位都有一定程度的下降，部分地区水位下降明显。因此，对水位持续下降地段，应采取科学管理和一定的控制措施。此外，由于地下水开采量的增加以及当地混合开采地下水，致使潜水层和承压水层直接产生了水力联系，承压含水层的承压性质正在减弱，向潜水转化。

图6.14 2012年9月1日潜水等水位线图

针对上述情况，在地下水开发利用中应采取一定的对策：①要科学地布井，有计划地进行浅、中、深层分层开采地下水；②开展地下水长期监测，以便掌握区域水情的变化趋势，发现问题及时解决；③加强水资源管理工作，合理开发利用地下水资源，确定开采量的大小和开采的时段长度，控制超量开采局面。

图6.15 2012年9月1日承压水等水位线图

‘贰’ 数值模拟方法

如前所述，基于镜质组化学反应动力学的煤热演化史数值模拟方法经历了简单函数关系模式、受热时间－经验法模式、反应活化能－温度函数模式和平行反应化学动力学模式四个发展阶段（秦勇等，1995）。其中，目前在化石能源地质界最为常用的是TTI法、LOM法和Tissot法，近年建立起来的EASY%R_o方法也已成为数值模拟的重要发展方向。

（一）时间－温度指数（TTI）法

这种方法的基本思想由Lopatin（1971）首先提出，后经Waples（1980）根据31个盆地402个样品的实测资料加以完善，建立起了其数值模拟的方法体系，故又称为Lopatin-Waples法。TTI法遵循化学反应动力学的基本法则来衡量温度和时间两个因素对煤化作用的贡献，即煤化作用速度随受热温度的增高而增大，温度每增加10℃煤化反应速度增加一倍。

由此，可将TTI数学模式用连续函数加以表示：

山西南部煤化作用及其古地热系统：兼论煤化作用的控气地质机理

T=T_s＋G·D

式中：T——煤层受热的古地温温度，℃；

T_s—－古地表温度，℃；

G——古地温梯度，℃/m;

D—煤层埋藏深度，m。

上述连续函数模式不便于常规运算，故在实际工作中通常采用分段积分方式，根据如下方程进行计算：

TTI＝∑γⁿ·△t

式中：γ——温度效应因子，考虑到温度每增高10℃煤级量值增加一倍的关系，取r=2;

n——温度指数冲仿，其大小取决于温度（表5-1）；

△t——煤层在某一温度段的受热时间，Ma。

表5-1TTI数学模式中温度指数（n）与煤层受热温度段的对应关系

注：温度每增（减）10℃，温度指数增（减）1。

TTI值本身可以作为煤级或有机质成熟度指标，Waples（1980）也提供了TTI值与镜质组最大油浸反射率的关系（表5-2）。此外，Wood（1988）通过回归拟合，将二者之间关系采用下式表示：

lg%R_o＝—0.006（lgTTI）³+0.042（lgTTI）²+0.162（lgTTI）—0.397

表5-2时间温度指数（TTI）与镜质组反射率（R_o）的相关关系

应予指出，TTI法尽管几十年来被化石能源地质界所广泛采用，但仍存在某些明显的不足有待于改进，国内外学者为此做过一定努力。例如，温度每碧猜增加10℃煤化速率加快一倍的假设有其局限性，只能在活化能10～25kcal/mol（相当于20～160℃的温度范围，1kcal=4186.8J，下同）区间是可靠的（Magoon,1983）。再如，沉积有机质热降解的反应活化能是受热温度的函数，随成熟度的增加而逐步加大，不同类型有机质达到相同成熟度所需的活化能也有差异，而TTI模式将整个煤化过程中的反应活化能作为常量，显然与实际情况有出入。

（二）悔判型有机成熟度水平（LOM）法

该方法由Hood（1975）建立，后经Bostick等（1979）完善而成为一种国内外广泛采用的数值模拟方法，故又被称为Hood-Bostick方法。

Hood等考虑到总体反应活化能随受热温度而增高的变化趋势（18～33kcal/mo1），采用镜质组油浸反射率标定有机质成熟度水平（LOM）（图5-1）。在模式中采用了有效受热时间的概念，即温度不低于最高受热温度15℃范围内的受热时间，建立起有机质成熟度、温度及有效受热时间之间的相互关系（图5-2）。

图5-1有机质成熟度水平（LOM）与其它有机质成熟度指标的关系（引自周中毅，1990）

Ⅰ—褐煤；Ⅱ₁—亚烟煤；Ⅱ₂—高挥发性烟煤；Ⅱ₃—中挥发性烟煤；Ⅱ₄—低挥发性烟煤；Ⅲ₁—半无烟煤；Ⅲ₂—无烟煤

1Btu=1055.06J

在实际工作中，LOM法采用图解方法求取有机质成熟度或镜质组反射率（图5-2）。国内外应用效果表明，利用这种方法推算出的盆地古地温温度是比较可靠的，与由综合研究所得出的结果以及盆地实际情况较为吻合（Vote,1981；周中毅等，1983,1984,1985）。

‘叁’ 数值模拟技术简介

（一）研究现状

地下多相、多组分流体运移数值模拟是在质量和能量守恒的基础上，建立的多相流体运动以及反映地球化学运移扩散的数学模型，通过离散建立大量的线形或非线形方程组，然后利用计算机计算求解，再通过图像显示模拟结果，达到对工程问题和物理问题乃至相关其他问题研究的目的。CO₂地质封存数值模拟就是利用计算机模拟的方法，来解决CO₂进入地质封存系统后运移、转化、水-岩-气之间的相互反应、CO₂泄漏对浅部含水层影响及诱发的储盖层物性变化等一系列问题，从而指导CO₂地质封存工程的实施。

目前，国内外已开展的关于CO₂地质封存数值模拟的研究工作包括以下几个方面：

1.超临界CO₂-水多相流体运动模拟

Pruess等（2003）模拟了均质各向同性咸水含水层中以恒定流量灌注CO₂条件下，灌注井井周非等温径向流情况。当忽略重力和惯性力效应时，模拟结果中存在相似变量ζ＝R²/t（其中，R为径向流动距离，t为时间），CO₂饱和度、溶解CO₂质量分数、沉淀盐的体积分数及流体压力都是相似变量的函数。这与O' Sullivan（1981）及Doughty和Pruess（1992）的结果一致。两相流的模拟考虑了CO₂和水的相对渗透率及毛细管力作用问题（Van Genuchten，1980），考虑了流体密度、黏度和CO₂溶解性随压力、温度和含盐度的变化，以及盐的沉淀导致含水层渗透率的减小等因素。

Doughty和Pruess（2004）利用Fro咸水含水层封存CO₂监测资料，反推了CO₂灌注后发生的物理和化学过程。他们采用TOUGH2数值模拟软件对两相（液、气）三组分（CO₂、水和溶解NaCl）系统进行模拟。考虑15MPa和65℃条件下，超临界CO₂在咸水中为非混溶流体，并能部分溶解于咸水的情况，分析了多相流系统边界设定的影响及相对渗透率取值问题，即模拟中对侧向边界的设置为均开（或均闭），结果导致压力的模拟结果与实际相比过低（或过高）。研究表明，由于上覆断层对CO₂的封堵作用，侧向边界对CO₂弥散羽的影响不大。模拟结果还显示，相对渗透率函数对CO₂弥散羽的演化有很强的影响。如何确定一个合适的相对渗透率以表征CO₂灌注咸水含水层的变化，仍是亟待解决的问题。Doughty和Pruess模拟了两种CO₂残余饱和度条件下的CO₂羽扩展随时间的变化，发现存在较大差异。残余饱和度较大的情况下，CO₂羽表现出紧缩状，在浮力作用下运移较慢；相反，在残余饱和度较小的情况下，CO₂羽流弥散很快，溶解性显着提高。

2.多组分反应地球化学运移模拟

水-砂岩-CO₂相互作用往往形成一系列次生矿物，或次生矿物组合。Worden et al.（2006）通过岩石学和CO₂灌注长石砂岩的地球化学模拟表明，北海Magnus油田上侏罗统浊积亚长石砂岩中的铁白云石、高岭石和石英可能具有成因联系。其中，铁白云石中的碳来自有机成因的CO₂。Watson et al.（2004）通过CO₂气与CH₄气储集砂岩的比较岩石学研究，证实澳大利亚Otway盆地Ladbroke Grove CO₂气储集砂岩中与CO₂气灌注有关的次生矿物组合为铁白云石-高岭石-次生石英。

Xu et al.（2005）利用一维砂岩-页岩系统模型模拟了储层中灌注的CO₂与矿物发生的化学反应过程，以及对储层环境的影响。模拟显示，在砂岩环境下，CO₂主要被方解石所固定，而方解石的沉淀导致孔隙度减小，进而导致渗透率相应减小。10万年间，砂岩封存能力达到90kg/m³的封存能力，这些被矿物固定下来的CO₂可以永久封存。Zwingmann等运用地球化学模拟软件EQ3/6进行的水-矿物-CO₂相互作用模拟也表明，若将CO₂灌注到日本本州岛中北部新潟盆地更新世灰爪组砂岩，CO₂以溶于水和形成碳酸盐矿物两种形式封存，其中后者封存量最大为21.3mol/kgH₂O，可达总封存量的90％，形成的碳酸盐矿物中也出现了片钠铝石。

3.耦合岩石力学模拟

从目前发表的论文及各国研究计划的综合报告上看，在CO₂咸水含水层封存研究方面，对于CO₂运移机制的分析和模拟很少考虑应力场的耦合作用。事实上，CO₂灌注压力和超临界CO₂的浮力作用将改变地层应力状态，即CO₂在上浮运移和侧向扩散过程中，孔隙压力可能会对原始裂隙和断裂产生影响；CO₂在咸水含水层中的长时期（千年级尺度以上）的封存，将改变含水层的地球化学状态，CO₂-咸水-含水层矿物的化学作用将可能导致岩体力学和水力学性质发生变化。

日本因位于4大板块交界处与环太平洋构造带中，活断层密集发育，地震频繁，地应力分布复杂，在CO₂地质封存评价方面，非常重视CO₂地质封存的力学稳定性研究（李琦等，2002；李小春等，2003）。李琦等（2002；2004；2006）提出了一个考虑初始地应力场、灌注压力、CO₂浮力及含热传导作用的热-水-力（THM）耦合模拟框架，考虑盖层底部附近存在不同倾角断层的二维平面应变地质封存问题。采用有限元算法，对灌注CO₂流体对断层稳定性的影响进行模拟分析。计算结果表明，为了避免断层位移需要特别注意对灌注压力的控制，因为CO₂灌注压力对断层滑动的影响远大于CO₂羽流浮力带来的影响。停止灌注CO₂后，CO₂羽流的上升则成为应力场扰动的主要因素。

（二）主要软件介绍

近年来，计算机模拟技术在许多研究领域得到了广泛的应用，开发出了许多优秀的模拟软件和程序。同样，可用于研究CO₂地质封存的数值模拟软件也很多，主要有PHREEQC、GEM、ECLIPSE、TOUGHREACT、PetroMod、MUFTE-UG和NUFT等，它们都有各自的特点和适用性。在进行数值模拟之前，需对这些数值模拟软件进行评价分析，选择适用于所研究问题的模拟软件。现对国际上常用的几款软件简介如下。

1.PHREEQC

PHREEQC是一款用于计算多种低温水文地球化学反应的计算机软件。以离子缔合水模型为基础，PHREEQC可完成以下任务：(1)计算物质形成种类与矿物的溶解饱和指数；(2)模拟地球化学反演过程；(3)计算批反应与一维运移反应。另外，与多组分溶质-运移模型耦合的PHREEQC可生成PHAST，一个用于模拟地下水流系统的三维反应-运移模拟器。但由于PHREEQC是在单相水流的基础上建立的模型，因而不能模拟超临界CO₂-水的两相流运动。

PHREEQC最简单的应用就是计算溶液中各种化学物质的分布，以及溶液中矿物与气体的饱和状态。反演模拟功能可推导和量化在流动过程中，能够反应化学物质变化的化学反应方程。PHREEQC可处理的反应方程包括建立矿物、表面配合物、阳离子交换剂、土壤溶液、气体组分单位分压、给定压力或给定体积气相间平衡的物质运移反应。在模拟这些均衡反应的同时，PHREEQC还可以模拟动力化学与生物反应，以及模拟从简单的线性衰变（代谢物降解或放射性衰变）到复杂的依赖于溶液化学组成和微生物数量确定的反应速度。这些反应处理功能可在批反应模拟或一维对流、弥散、反应型运移模拟中使用。

2.GEM

GEM v.2009.13（Nghiem et al.，2004）是一款用来模拟利用CO₂和酸性气体提高石油采收率的模拟器，该模拟器完全耦合了地球化学组成状态方程。GEM采用一步求解法进行状态方程的求解。GEM可以用来模拟：对流和弥散流体、油（或超临界CO₂）、气和咸水间的平衡、水相物种间的化学平衡，以及矿物的动态溶解和沉淀。该模拟器采用自适应的隐式离散技术利用一维、二维或者三维模型来模拟孔隙介质中溶质的运移。油相和气相用一个状态方程来模拟，气体在水相的溶解度采用亨利定律模型来计算。水通过蒸发进入到气相、盖层的穿透、热效应和裂隙的封闭作用也可以利用GEM来模拟。

3.ECLIPSE

ECLIPSE是一个并行化的可以模拟黑油、组分、热采等问题的成熟软件。1994年，胜利石油管理局引进了ECLIPSE油藏数值模拟串行软件，广泛开展了从油藏到气藏，从常用油田到特殊油气田、从常规模拟研究到特殊模拟研究等多方面的应用。主要模块有主模型、黑油、组分、热采、流线法、运行平台和ECLIPSE Office等。

ECLIPSE是一个商业软件，在使用中其内核部分是封闭的，使用者只能将其作为一个“黑箱”来操作。其不足之处有：不可能免费的获得和随意地使用和修改；无法耦合最前沿的地质流体热力学模型；无法加入更多影响因素来研究具体问题。因此，ECLIPSE不适宜用于前沿科学研究。

4.TOUGH2/TOUGHREACT

TOUGH2是Transport of Unsaturated Groundwater and Heat（非饱和地下水流及热流传输）的英文缩写，是一个模拟一维、二维和三维孔隙或裂隙介质中，多相流、多组分及非等温的水流及热量运移的数值模拟程序。TOUGH2使用积分有限差分（Integral Finite Differences，IFD）（图3-8）的方法来解决多相流动和多组分化学运移模拟中的空间离散化问题（Pruess et al.，1999s；Xu et al.，2004）。为了满足大规模计算需要，Zhang et al.（2008）开发了TOUGH2的平行计算版本，即TOUGH2-MP。

该方法在对地质介质的离散化上是比较灵活的，允许使用不规则的网格，十分适合对多区域非均质系统和裂隙岩石系统中流体流动、运移和水岩相互作用的模拟。而对于规则的网格剖分，积分有限差分方法相当于传统的有限差分法。其中，对于任意区域Vn，它的质量（对于水、气体和其他化学组分）和能量（对于热）守恒方程可以用积分的方式（式3-5）表达：

图3-8 积分有限差分法中的空间离散化和几何参数数据构成图

中国二氧化碳地质封存选址指南研究

式中：下角标n为表示一个单元格；下角标m为表示和单元格n互相连接的网格m；Δt为时间步长；M_n为单元格n的平均质量或能量密度；A_nm为单元网格n和m交界的面段；F_nm为通过面积为A_nm的质量或能量通量；q_n为单元格n内单位体积的平均源汇率。

许天福等（1998）在TOUGH2的框架基础之上，增加了多组分溶质运移和地球化学反应的模拟功能，形成了一套较为完善的可变饱和地质介质中非等温多相流体反应地球化学运移模拟软件——TOUGHREACT。该软件不仅包括了TOUGH2的全部功能，而且适用于不同温度、压力、水饱和度、离子强度、pH值和氧化还原电位（Eh）等水文地质和地球化学条件下的热-物理-化学过程。还可以应用于一维、二维或三维非均质（物理和化学的）孔隙或裂隙介质中的相关数值模拟研究。在理论上可以容纳任意数量的以固相、液相或气相存在的化学组分（但是在实际模拟中会受到计算能力和计算时间等硬件条件的限制），并且考虑了一系列化学平衡反应，如溶液中的配合反应、气体的溶解或脱溶、离子吸附作用、阳离子交换及受平衡控制或反应动力学控制的矿物溶解或沉淀反应等。可以说TOUGHREACT、是TOUGH2的升级版，近年来在世界范围内CO₂地质封存研究和工程实践中得到了广泛的应用。

除包含TOUGH2所有的功能外，TOUGHREACT还可以应用于一系列的反应性流体和地球化学迁移问题。比如：(1)伴随Kd线性吸附和放射性衰变的污染物迁移问题；(2)在周围环境条件下，自然界中地下水的化学演变；(3)核废料处置地点评估；(4)深部岩层的沉积成岩作用；(5)CO₂地质处置。多相流体运动，多组分反应地球化学，各种封存形式封存量以及随时间空间变化；(6)矿物沉积（如表生铜矿富集）；(7)自然和补给环境下热水系统中的矿物变化。

通过最近几年相关研究者的不懈努力，TOUGHREACT在实际应用中得到了进一步完善和提高，增加了部分新功能，如水相内部反应动力学和生物降解作用，改进了矿物-水反应表面积计算方法，以及气-水反应中气的活度系数的修正等。

5.PetroMod

由德国IES（Integrated Exploration System）公司研究开发的PetroMod多组分、多相态的多维含油气系统模拟软件综合平台已被世界石油业所公认。该软件融入了断层活动性、盐丘上涌和刺穿、火山岩的侵入、气体扩散效应、油气水三相运移和油气吸附模型等相关技术。

该模拟软件平台推出和采用的油气运移组合模拟算法（Hybird）是当今最先进的油气运移模拟算法，既可以保证模拟的精度，又可以极大地提高模拟的运算速度。其中的PetroFlow3D用于油气运移、聚集、圈闭和散失等情况的模拟，同时PetroCharge Express为我们提供了基于图件的油气运移和圈闭模拟的快速分析工具。

6.MUFTE-UG

MUFTE-UG是MUFTE和UG.MUFTE的结合。MUFTE即多相流（Muliphase Flow）、运移（Transport）和能量（Energy）模型。该软件包主要包括物理模型概念和孔隙裂隙介质中等温和非等温多相多组分流动和运移过程的离散方法（Helmig，1997；Helmig et al.，1998）。它能对裂隙孔隙介质进行离散性描述（Dietrich et al.，2005）。UG是非结构性网格（Unstructured Grid）的缩写，它提供的数据结构能快速解算以平行、自适应多网格法为基础的离散型偏微分方程。具有模块化结构的MUFTE-UG很容易解决各种有特殊要求的问题。

模块化结构的MUFTE-UG具有许多不同的环境与技术应用。例如，在环境应用领域，MUFTE-UG能够模拟如下两个问题。

（1）NAPL（非液相流体）向饱和与非饱和土壤的渗流。优化改进的修复技术在MUFTE中具有广泛的研究和发展空间。

（2）地下CO₂的消散。CO₂以高温高压灌注地表以下几百米的地层中，MUFTE-UG可用于非均质含水层（对流和弥散运移）中羽状体演化评价，伴随温度效应（由于膨胀和压缩）和组分间相互溶解（卤水和CO₂）。

7.NUFT

NUFT（Nonisothermal Unsaturated-Saturated Flowand Transport model）是一套用来解决在多孔介质中多相、多组分非等温流动和溶质运移过程中地下污染物运移的数值解法器。此软件利用简单的代码来利用通用的实用程序和输入文件的格式。最近，此代码在Unix和DOS系统下运行成功。

该程序利用一套完整的有限差分空间离散法求解平衡方程组。每一个时间步长内利用Newton-Raphson方法求解非线性方程组，而在每一步迭代过程中利用直接解法和预共轭梯度法求解线性方程组。该模型可以解决一、二和三维水流及溶质运移问题。将来该模型会耦合进毛细滞后、非正交网格离散、有限单元剖分和固体非线性等温吸附等功能。

（三）研究方法

通常情况下，CO₂地质封存数值模拟包括以下主要过程。

（1）建立概念模型：根据各种方法获取的实际资料来概化和建立CO₂地质封存概念模型，包括边界范围、地层或储盖层高程、储盖层确定、参数及分区、源汇项、主要物理化学过程以及模型维度（一维、二维和三维）。

（2）建立数学模型：建立一套描述深部咸水层中多相流动和多组分反应性溶质运移的偏微分方程组，包括初始条件和边界条件问题。

（3）模型离散化：把概念模型中的各种信息通过网格剖分进行离散，形成大量的网格单元，然后通过有限差分、有限单元和积分有限差分等方法转化成单元的质量和能量守恒方程组，再用多种方法将非线性方程组线性化，形成线性代数方程组，然后求解方程组。

（4）模型识别和校正：根据模型计算结果和实际监测数据进行对比拟合，适度合理调整参数，使模型能够综合反映实际情况。在历史拟合过程中出现较大误差，应重新检查概念模型，修正概念模型。对所建模型进行参数敏感性分析，对于较敏感的参数应该慎重选取，甚至需要做大量的试验来确定。

（5）模型预测：建立了可靠的模型后，便可以进行模拟预测。

数值模拟的关键是地质模型概化、计算精度和计算速度。由于计算的精度取决于离散的程度，而离散的程度又决定了计算的速度，这是一对矛盾，要根据解决问题的需要来选择离散化的程度和计算速度。

CO₂在储层中的运移、溶解以及与围岩的化学反应形成了一个多相、多组分的反应体系，涉及的主要数学方程有超临界CO₂-水的两相流体运动控制方程、溶质运移控制方程和化学反应方程等。建立数值模型时，通常采用有限差分法、有限元法和积分有限差分法等。

由于实际应用时多采用已有的数值模拟软件对CO₂地质封存的某一过程进行模拟，不涉及软件的开发及程序代码的编写，只需根据研究的需要选择合适的软件进行模拟预测，而软件一旦选定，数学模型和数值模型基本上已经确定。

‘肆’ 数值模拟的计算机方法

有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将 求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级 数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而 建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数 问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分 的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可 以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式 的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步 长一般根据实际地形的情况和柯朗稳定条件来决定。

构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达 式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等， 其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几 种不同差分格式的组合，可以组合成不同的差分计算格式。

有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分 方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式 ，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形 网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合 同样构成不同的有限元计算格式。对于权函数，伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数 ；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域 内选取N个配置点 。令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。

对于有限元方法，其基本思路和解题步骤可归纳为
(1)建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。
(2)区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。
(3)确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插值条 件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元 具有规则的几何形状，在选取基函数时可遵循一定的法则。
(4)单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将 近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数(即单元中各节点 的参数值)的代数方程组，称为单元有限元方程。
(5)总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进 行累加，形成总体有限元方程。
(6)边界条件的处理：一般边界条件有三种形式，分为本质边界条件(狄里克雷边界条件 )、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件， 一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法 则对总体有限元方程进行修正满足。
(7)解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭 方程组，采用适当的数值计算方法求解，可求得各节点的函数值。

有限体积法（Finite Volume Method）又称为控制体积法。其基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分，必须假定值在网格点之间的变化规律，即假设值的分段的分布的分布剖面。从积分区域的选取方法看来，有限体积法属于加权剩余法中的子区域法；从未知解的近似方法看来，有限体积法属于采用局部近似的离散方法。简言之，子区域法属于有限体积发的基本方法。

有限体积法的基本思路易于理解，并能得出直接的物理解释。离散方程的物理意义，就 是因变量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小的控 制体积中的守恒原理一样。 限体积法得出的离散方程，要求因变量的积分守恒对任意一组控制体积都得到满足，对整个计算区域，自然也得到满足。这是有限体积法吸引人的优点。有一些离散方法，例如有限差分法，仅当网格极其细密时，离散方程才满足积分守恒；而有限体积法即使在粗网格情况下，也显示出准确的积分守恒。就离散方法而言，有限体积法可视作有限单元法和有限差分法的中间物。有限单元法必须假定值在网格点之间的变化规律（既插值函数），并将其作为近似解。有限差分法只考虑网格点上的数值而不考虑值在网格点之间如何变化。有限体积法只寻求的结点值，这与有限差分法相类似；但有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布，这又与有限单元法相类似。在有限体积法中，插值函数只用于计算控制 体积的积分，得出离散方程之后，便可忘掉插值函数；如果需要的话，可以对微分方程 中不同的项采取不同的插值函数。

‘伍’ 数值模拟流程

不同的软件进行数值模拟时所需的参数、计算方法、剖分格式等不尽相同，数值模拟的过程也不同，但大致相同，本文以TOUGHREACT为例介绍CO₂地质储存数值模拟的流程。

（一）研究范围的确定

一般情况下，独立的天然水文地质系统是计算区最好的选择，它具有自然边界，便于较准确地利用其真实的边界条件，避免人为边界在资料提供上的困难和误差。但是在实际工作中，常常不能完全利用自然边界，这就需要充分利用勘察和长期观测资料等建立人为边界。在确定计算区域时，除了保证范围足够大以外，还应使假定的边界条件尽可能接近真实状态。

计算范围的划定应充分考虑研究目的、区域地质构造、储层岩性、储层岩石矿物组成及地下水化学成分等多方面因素。数值模拟时间根据研究目的不同稿谨具有不同的时间尺度。就CO₂地质储存数值模拟而言，如果不考虑地球化学作用，封存系统在1000年数量级的模拟时间内基本上已达到平衡或稳定。在划定边界时还应考虑CO₂在储层中的扩散距离，与研究区地质模型的孔隙度，渗透率等参数关系密切。为了保证所选模型范围边界在模拟期内不影响模拟结果，尽量通过具有相同地质条件的天然CO₂气田（藏）进行类比，确定大体的计算范围的边界。如果考虑地球化学反应，由于CO₂注入引发的水－岩－气反应对围岩岩性改变较显着，制约着CO₂注人的速度和径向运移的距离等。

（二）明确研究目的

在进行数值模拟以前首先要明确利用数值模拟技术要解决什么样的问题。对于CO₂地质储存工程而言，进行数值模拟的目的主要是在CO₂地质储存工程实施前，通过数值模拟技术对工程的选址、方案设计进行优化，工程实施期技术指导、运冲森行期监测及后期CO₂泄漏的风险评估等进行预测，以指导项目科学、合理地实施，将CO₂泄漏风险降至最低。

研究目的决定着前期资料的收集类型、地质建模的侧重点、地质模型离散的精密度以及初始、边界条件的处理方式等过程。

（三）资料的收集整理

1）通过遥感、综合地质调查、物探、钻探和各类样品测试分析等手段获取场地深部地层岩性、地质构造、水文地质、水文地球化学、岩石矿物资料和数据；

2）搜集和分析CO₂地质储存场地地质岩性、区域构造格架、活动断层与地震活动情况等；

3）采用钻井岩心、测井和地震反射方法，调查CO₂地质储存场地目标储层和盖层的空间分布形态，埋深、厚度和规模等；

4）使用X射线衍射、扫描电镜等方法研究分析封存场地岩石矿物组成、孔隙结构特征及其物理化学性质；

5）通过采取浅部、深部含水层水样进行水质全分析，获得储盖层地层水及浅部含水层初始水化学成分。

不同的数值模拟软件其数学模型的数值解法不同，空间离散方式也不尽相同，所需的模型参数也有一定的差异，表9-1即为TOUGHREACT数值模拟所需要的主要参数。

表9-1 CO₂地质储存模拟过程中需要的主要参数（以TOUGHREACT为例）

图9-3 网格剖分

网格剖分对计算的精度及计算的效率有很重要的影响。精度越高对模拟结果刻画的越精细，但是数据的计算量越大，对计算机的要求也越高。建议在进行地质模型剖分时先采用较粗的网格剖分，如果这种剖分方式下模拟结果合理然后再进行精细化剖分，用于对模拟结果更加详细的刻画。

2.参数和初始条件

初始条件是指在初始时刻（t=0）时研究区内求解数学模型主要状态变量的初始值。选择的应用软件不同所需的状态变量数量、种类不同。如TOUGHREACT所需的初始主要状态变量包括压力、温度和组分浓度的空间分布。地质参数包括孔隙度、渗透率、密度、压力、温度、毛细压力等参数值。这些数值一部分采用室内实验测得，另一部分采用参考文献的经验值；地层水的化学成分的初始值采用实际地层水的化学分析，主要是8大离子的浓度、盐度和pH 等。如果研究区深部地层中的水样难以获得，如盖层，则采用静态平衡的方法，利用具有与储层相同盐度的咸水与含有原生矿物的地层岩石在原地层环境下进行化学反应，获取平衡状态下的地层水化学成分的初始值；通过岩矿分析、电子扫描、Ⅹ衍射等手段，获得组成CO₂地质储层盖层的原生矿物成分体积含量初始值，并根据原生矿物的组成合理判断次生矿物。

从原则上讲，初始时刻是可以任意取定的，只要该时刻所需的参数和状态变量值已知即可。因此我们不应该把初始条件理解为研究系统的初始状态。具体如何取，应该视问题的需要、资料来源、计算方便与否等因素而定。

3.边界条件

边界条件是某一实际问题数学模型具有定解的必要条件之一。地下水流问题和溶质运移问题边界条件的定义不尽相同，但一般概化为以下三种。

（1）一类边界条件（Dirichlet条件）

解决水流问题时，此类边界条件为在边界上所有点的水头是给定的；对于溶质运移问题，一类边界条件是指研究区边界上的溶质浓度分布已知。解决CO₂—水两相流动问题时，此类边界条件为在边界上所有点的压力是给定的。

（2）二类边界条件（Neumann条件）

当已知某一边界的单位面积流入或流出的流量时，可视作解决流动问题的二类边界；相对溶质运移来讲，此类边界又称给定弥散通量边界，即边界上的弥散通量随时间变化规律已知。

（3）三类边界（Cauchy条件）

当研究区一部分满足一类Dirichlet条件，而另一部分满足二类Neumann条件时，这类问题称为混合边界问题，称为三类边界。对溶质运移而言，此类边界为边界上溶质通量随时间变化规律已知。

在CO₂地质储存数值模拟过程中，由于储层地层多在800m以下，地质模型的顶部和底部根据实际需要可以处理为不透边界；为了避免边界对模拟结果的影响，研究区的范围一般比实际CO₂所能运移到的范围大得多，因此，在处理四周边界时一般设置为无穷一类边界或不透边界。在确定边界条件时，应根据水文地质条件以及现有的资料来综合考虑。

4.源汇项处理

在多孔介质中流动和溶质运移的问题中，对流、水动力弥散和溶质源或／和汇，是决定含水层中任一内点上溶质质量时变率的两大因素。源汇项问题在水质与水量计算中以及正确处理对流－弥散方程和渗流基本微分方程中占有重要地位。作为源汇项的方式很多，如越流补给、含水层弹性释放补给以及抽（注）井的补给等。

对于深部咸水层CO₂地质储存系统而言，系统的顶部一般为具有低渗、低孔的泥岩、页岩等致密性岩层，越流补给较难发生。整个CO₂地质储库系统的源汇项主要指对流（如侧向边界）和抽（注）井。

（八）模型的校正与验证

模型识别是建立地下流体数值模型最重要的环节之一，正确理解和进行拟合对于提高数值模型的仿真性是至关重要的。在有实测结果的情况下如示范工程，可将模拟结果与实测结果进行比较，对相关参数进行适当合理的调整，使模拟结果在给定的误差范围内与实测结果吻合。若误差较大，应该重新检验概念模型的可靠性，甚至重新建立概念模型。在识别校正以后，应采用校正好的模型继续计算，并与未用来识别校正的实际数据比较，验证模型的准确性和可靠性。若存在较大误差，需重复前面的过程。在没有实测结果的情况下，数值模型的可靠性可通过类比相关资料或根据个人经验和理论判断。

（九）模拟预测

模型预测是实施数值模拟技术的主要目的。对于CO₂地质储存工程而言，由于CO₂地质储存技术的提出为时尚短，针对CO₂在深部咸水层中的运移、扩散、与地层水和围岩产生的化学反应，以及由于CO₂灌注引起的储盖层物理、化学性质变化研究均处于研究和发展阶段。因此，在工程实施过程中急需具有技术指导性的工具产生，避免造成投资的浪费和CO₂泄漏等风险的出现。

利用经过识别校正与验证过的数值模型对CO₂地质储存过程进行模拟预测，有针对性地对模拟数据进行后期处理，如统计分析、比较等手段对结果进行解译，以此达到场地的优选，目标储层灌注能力、储存潜力的评估，CO：扩散运移途径和速度、不同捕集方式封存量及它们之间的时空转化等过程的详细刻画与模拟仿真等目的。同时可以预测CO₂在已有、重新激活或新生成的裂隙中逃逸的可能性及时间、CO₂泄漏风险评估以及评价CO₂泄漏对浅层地下水的水质、水量及对地表环境的影响等。

上述结果的分析只是数值模拟技术所能解决问题的冰山一角。对于数值模拟结果的处理要根据所研究的目的进行有针对性的提取和解译。通过对处理后的数据进行总结分析，发现问题从而解决问题，并掌握内在规律，为CO₂地质储存工程的前期设计、工程实施、中期监测管理提供理论支持和科学的技术指导，并可以提前开展风险预测，尽早制定预案防范CO₂地质储存工程实施及运行过程中可能出现的隐患。

‘陆’ 模拟方法简介

4.4.1.1 FLAC^3D的基本原理

数值计算一直是工程地质学分析方法中的一个重要手段，特别是近年来，随着计算机技术和数值计算方法的迅速发展，从二维到三维，从静态到动态，计算精度、可靠度不断提高，加之大量工程实例的验证，使数值模拟计算方法已日益成为岩体稳定评价中不可或缺的重要方法之一。本研究利用国际工程地质界公认的最新通用软件FLAC^3D2.0作为工具，对典型岩溶塌陷类型的工程环境效应进行评价。

FLAC^3D分析在求解中使用如下3种计算方法：

（1）离散模型方法，连续介质被离散为若干六面体单元，作用力均被集中在节点上。

（2）有限差分方法，变量关于空间和时间的一阶导数均用有限差分来近似表达。

（3）动态松弛方法，由质点运动方程求解，通过阻尼使系统运动方程衰减至平衡状态。

4.4.1.1.1 空间导数的有限差分

快速拉格朗日分析采用混合离散方法，将区域离散为常应变六面体单元的集合体，又将每个六面体看作以六面体角点为角点的常应变四面体的集合体，应力、应变、节点不平衡力等变量均在四面体上进行计算，六面体单元的应力应变取值为其内四面体的体积加权平均，这种方法既避免了常应变六面体单元常会遇到的位移剪切锁死现象，又使得四面体单元的位移模式可以充分适应一些本构的要求，如一四面体，节点编号为1到4，第n面表示与节点n相对的面，设其内一点的速率分量为υ_i，由高斯公式得

岩溶塌陷机理及其预测与评价研究

式中：V为四面体的体积；S为四面体的外表面；n_j为外表面的单位法向向量分量。对于常应变单元，υ_i为线性分布，n_j在每个面上为常量，由式（4-1）可得：

岩溶塌陷机理及其预测与评价研究

式中：上标l表示节点l的变量；上标（l）表示面l的变量。

4.4.1.1.2 运动方程

快速拉格朗日分析以节点为计算对象，在时域内求解，节点运动方程可表示如下：

岩溶塌陷机理及其预测与评价研究

式中

（t）为在t时刻l节点在i方向的不平衡力分量，可由虚功原理导出；m^l为l节点的集中质量。对于静态问题，采用虚拟质量以保证数值稳定，而对于动态问题则采用实际的集中质量。

将上式左端用中心差分来近似，则可得：

岩溶塌陷机理及其预测与评价研究

4.4.1.1.3 应变、应力及节点不平衡力

快速拉格朗日分析由速率来求某一时步的单元应变增量，即

岩溶塌陷机理及其预测与评价研究

有了应变增量，即可由本构方程求出应力增量，进而得到总应力。

4.4.1.1.4 阻尼力

对于静态问题，在式（4-3）的不平衡力中加入了非粘性阻尼，以使系统的振动逐渐衰减直至达到平衡状态（即不平衡力接近零），此时式（4-3）变为：

岩溶塌陷机理及其预测与评价研究

阻尼为：

岩溶塌陷机理及其预测与评价研究

式中：α为阻尼系数，其默认值为0.8；而符号函数

岩溶塌陷机理及其预测与评价研究

4.4.1.1.5 计算循环

由以上可以看出快速拉格朗日分析的计算循环如图4-2所示。

图4-2 FLAC^3D软件计算循环图

4.4.1.2 FLAC^3D的功能和应用范围

尽管FLAC^3D的计算公式源于有限差分方法，但其计算结果与有限元方法的计算结果（对于常应变四面体）相同，而且它与现行的数值方法相比有着明显的优点（黄润秋等，1991、1995、1999a、1999b），其特点包括：

（1）“混合离散化”（mixed discretization）技术的应用，更能精确和有效地模拟计算材料的塑性破坏和塑性流动，在力学上比常规有限元的数值积分更为合理。

（2）全部使用动力运动方程，即使在模拟静态问题时也如此。因此，它可以较好地模拟系统的力学不平衡到平衡的全过程，实现动态的模拟过程。

（3）求解中采用“显式”差分方法大大节约了计算时间，特别对求解任意的非线性应力-应变问题尤为重要。同时它不需要存储较大的刚度矩阵，因此，它与一般的差分分析方法相比，既节约了计算机的内存空间又减少了运算时间，大大提高了解决问题的速度。

（4）物体由多面体单元所表示，可以通过调整三维网格的方法，以适应研究体真实的形状，每个单元力学行为是对应力-应变法则和边界力、约束条件的响应。材料能产生屈服和流动，而且网格也能变形（in large-strain mode），并随材料移动。

（5）强大的后处理功能。能根据需要输出设定阶段的应力应变成果，提取各工程部位的应力、应变的量值，同时还可以根据追踪功能（history）获取所需变量值的历时变化曲线。

研究中使用的FLAC^3D所使用的模型主要包括近10种模型，常用的有以下几种：

（1）空模型（null）。用空材料来表示开挖后的材料，在空单元中的应力被自动地赋零。

即：

岩溶塌陷机理及其预测与评价研究

用此种模型模拟实际中的洞穴。

（2）弹性均质模型（elastic，isotropic）。在FLAC^3D中，对于基岩主要采取了这种模型，而对于非均质的介质再根据其不均匀情况进行特殊的定义。在这种模型中，据虎克（Hooke）定律，有：

Δσ_ij=2GΔε_ij+α₂Δε_kkδ_ij

在此Einstein求和中，δ_ij为位移，α₂为与材料的体积模量及剪切模量相关的材料系数：

岩溶塌陷机理及其预测与评价研究

新的应力值通过下式获得：

岩溶塌陷机理及其预测与评价研究

（3）弹性非均质模型（elastic，orthotropic），主要用于非均质情况的模拟。

（4）Drucker-Prager塑性模型。在FLAC^3D中此模型使用的破坏准则主要为Drucker-Prager拉张破坏准则，用以评价计算中的塑性破坏现象。

（5）摩尔-库仑模型（Mohr-Coulomnb）。在整个数值模拟过程中，对于土层、砂土、粉土等均采用了摩尔-库仑模型作为计算模型，即将其视为弹塑性材料来解决。在此模型中，以Mohr-Coulomnb准则作为材料的破坏准则。