数值模拟的计算机方法？数值模拟示例

本文目录

• 数值模拟的计算机方法
• 数值模拟示例
• 数值模拟技术简介
• 数值模拟
• 数值模拟流程
• 如何做好数值模拟
• 试说明数值模拟方法的特点，它与理论研究，实验研究有什么关系
• 数值模拟方法
• 什么是数值模拟方法与实验的优缺点呢

数值模拟的计算机方法

有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将 求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级 数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而 建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数 问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 　　对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分 的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可 以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式 的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步 长一般根据实际地形的情况和柯朗稳定条件来决定。　　构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达 式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等， 其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几 种不同差分格式的组合，可以组合成不同的差分计算格式。　　有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分 方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式 ，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形 网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合 同样构成不同的有限元计算格式。对于权函数，伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数 ；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域 内选取N个配置点 。令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。 对于有限元方法，其基本思路和解题步骤可归纳为 (1)建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。 (2)区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。(3)确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插值条 件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元 具有规则的几何形状，在选取基函数时可遵循一定的法则。(4)单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将 近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数(即单元中各节点 的参数值)的代数方程组，称为单元有限元方程。(5)总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进 行累加，形成总体有限元方程。(6)边界条件的处理：一般边界条件有三种形式，分为本质边界条件(狄里克雷边界条件 )、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件， 一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法 则对总体有限元方程进行修正满足。 (7)解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭 方程组，采用适当的数值计算方法求解，可求得各节点的函数值。 有限体积法（Finite Volume Method）又称为控制体积法。其基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分，必须假定值在网格点之间的变化规律，即假设值的分段的分布的分布剖面。从积分区域的选取方法看来，有限体积法属于加权剩余法中的子区域法；从未知解的近似方法看来，有限体积法属于采用局部近似的离散方法。简言之，子区域法属于有限体积发的基本方法。 有限体积法的基本思路易于理解，并能得出直接的物理解释。离散方程的物理意义，就 是因变量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小的控 制体积中的守恒原理一样。 限体积法得出的离散方程，要求因变量的积分守恒对任意一组控制体积都得到满足，对整个计算区域，自然也得到满足。这是有限体积法吸引人的优点。有一些离散方法，例如有限差分法，仅当网格极其细密时，离散方程才满足积分守恒；而有限体积法即使在粗网格情况下，也显示出准确的积分守恒。就离散方法而言，有限体积法可视作有限单元法和有限差分法的中间物。有限单元法必须假定值在网格点之间的变化规律（既插值函数），并将其作为近似解。有限差分法只考虑网格点上的数值而不考虑值在网格点之间如何变化。有限体积法只寻求的结点值，这与有限差分法相类似；但有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布，这又与有限单元法相类似。在有限体积法中，插值函数只用于计算控制 体积的积分，得出离散方程之后，便可忘掉插值函数；如果需要的话，可以对微分方程 中不同的项采取不同的插值函数。

数值模拟示例

（一）研究区概况

1.区域构造

储存工程区地质构造简单，出露地层单一，除了部分地段见侏罗系出露外，全区大部分均为白垩系及其后的沉积地层，地层产状基本水平，为由西向东逐渐缓升。全区未见岩浆活动及变质作用，属于标准的地台相沉积。

2.区域地层

以神华CCS示范项目二叠系上统石千峰组（F3s）为例，研究区内均有分布，总厚度为300m左右，自下至上划分为5段，分别简称为：千5、千4、千3、千2和千1段，与下伏上石盒子组和上覆刘家沟组均呈假整合接触关系。石千峰组的沉积环境为炎热氧化环境，岩性为一套以砖红色、紫红色为特征的陆源碎屑岩夹少量碳酸盐岩及硅质岩组合。上部为紫红、砖红色泥岩夹钙质条带及淡水灰岩薄层或结核层。中下部为灰白色中粗粒长石局部为细砾岩，砖红色泥岩呈不等厚互层。

3.岩石学特征和矿物成分

（1）岩石学特征

通过岩石样品的薄片鉴定分析结果可知，石千峰组储层主要是棕红色、红色和浅红色的砂岩，主要的岩石类型为岩屑长石砂岩，岩石结构为碎屑结构。碎屑颗粒的粒度在0.1～2.2mm 之间，主要在0.7～1.2mm之间，为巨—粗粒碎屑结构。碎屑颗粒的分选较差、磨圆中等，碎屑颗粒之间以点接触为主，其次为线接触。孔隙式胶结，颗粒支撑。

（2）矿物成分

由所收集的石千峰组岩石矿物X－衍射分析数据可知，石千峰组储层的矿物成分为长石（21％～41%）、石英（24％～49%）；碳酸盐矿物主要为方解石（0％～20%），其次白云石（0％～2%）；黏土矿物为绿泥石（5％～10%）、蒙脱石（15％～20%）、伊利石（5％～10%）、高岭石（0％～5%），其他矿物主要有赤铁矿（0％～3％）。

（3）咸水水化学特征

石千峰组咸水层埋深大、咸水的矿化度高，与上下地层的水化学特征差异大，基本呈封闭状态。咸水的水化学类型为Ca-Na-C1型水，总溶解固体为31.2g/L,Ca介质）

（五）模拟结果

以石千峰组作为目标储层，利用数值模拟技术分别模拟了CO₂在地质储存中的储存潜力、扩散运移、不同捕集形式之间的时空转化等过程。同时，对比分析了考虑矿物地球化学反应对CO₂的注入、扩散等的影响。定量分析对CO₂捕集具有贡献的矿物。结论如下：

1）在不同的注入方式下，定流速注入方式注入量较大.CO₂径向扩散距离稍大于定压注入；

2）两种注入方式下，气态CO₂所占比例均从停止注入时的最高值开始下降，停止注入后100年，定速注入方式气态CO₂所占总注入量的比例大；

3）考虑地球化学反应过程时，注入量和定速率方式相同，CO₂在储层中径向扩散距离稍微小于不考虑化学过程的模型；

4）在捕集方式转化过程中，20年之前未发现矿物捕集，20年之后出现，100年时捕集量占总注入量的5.21％;

5）在矿物捕集过程中，起到固定CO₂作用的矿物在20年之前为方解石，20年之后为方解石和片钠铝石。

数值模拟技术简介

（一）研究现状

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

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

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

Pruess等（2003）模拟了均质各向同性咸水含水层中以恒定流量灌注CO₂条件下，灌注井井周非等温径向流情况。当忽略重力和惯性力效应时，模拟结果中存在相似变量ζ＝R的封存能力，这些被矿物固定下来的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₂地质封存的某一过程进行模拟，不涉及软件的开发及程序代码的编写，只需根据研究的需要选择合适的软件进行模拟预测，而软件一旦选定，数学模型和数值模型基本上已经确定。

数值模拟

数值模拟（数值法）是对数学模型的一种近似解法，它仅能求出计算域内有限点某个时刻水头的近似值，这个值在实际应用中可以满足精度要求。数值法可以解决许多复杂水文地质条件下的渗流计算问题，应用十分广泛。如用于大中型水源地、地下水的补径排条件复杂、渗流区形状不规则、含水介质为非均质各向异性等条件下，确定水头分布和流量计算。

（一）渗流区域离散化（以二维流为例）

采用数值模拟技术研究地下水的运动，首先将要研究的水文地质模型内的含水层离散化。所谓离散化，就是将要研究的渗流区非均质各向异性含水层，按照一定的方式剖分（分割）成许多相互联系的小均衡区，在每个小均衡区内是均质各向同性的。在每个小的均衡区内，其含水层参数视为常数；其中心水头值或有条件下的平均水头值视为小均衡区内水头代表值。剖分通常采用两种形式（矩形、多边形）进行。

1.矩形均衡域

它是用两组正交的平行线把均衡区分为许多小的矩形均衡域，如图7-3所示。在剖分时约定：①定水头或已知水头边界（一类边界）应从小均衡域的中心通过；②隔水边界（二类边界）与小均衡域的边界重合。这种剖分方法类似于直角坐标系，用适当的编号标定小区域及节（结）点（小均衡域的中心点）。常用的术语有：

图7-3 渗流区被剖分成矩形小均衡域

（据李俊亭等，1987）

1）点、行、列，点（节点）为小区域的中心点，网格的横向称行，竖向称列。

2）步长，分为空间步长（Δx，Δy，Δz）（图7-3）和时间步长（Δt）。

3）小区域及节点编号统一记为（i，j），表示小区域及节点位于第i行第j列。

2.多边形均衡域

由于多边形均衡域与复杂边界的几何形状比较接近，因此使用较多。它是先按三角形剖分渗流域，再以三角形为基础构成多边形均衡域，见图7-4。常用的术语及注意事项：

1）点元、面元、线元，三角形的边称线元，三角形的顶点称点元（节点或结点），三角形的面积称面元；

2）要求剖分时三角形的单个内角取30°～90°；

3）渗流区剖后的面积与原面积要吻合，既不要重复也不要开裂。

（二）基本均衡离散方程（以规则网格的有限差分方法为例）

将图7-3中的（i，j）的均衡区与相邻均衡域的水量交换关系表示在图7-5上。

图7-4 渗流区域三角形

图7-5（i，j）均衡区的流量关系示意

（据李俊亭等，1987）

1）均衡时段为Δt_n+1：

Δt_n+1=t_n+1-t_n0

表示点（i，j）上t_n时刻的水头。

2）若（i，j）均衡区内不存在垂向水量交替，则依据水均衡原理有：

地下水动力学

在x轴方向上不同均衡时段分别为：

地下水动力学

式中：

地下水动力学

3）考虑到式（7-14）与式（7-15）的不同，会产生不同的计算结果。计算方案（差分格式）将写出如下通式：

地下水动力学

式中：0≤θ≤1。θ常取3种情况：①当θ=0时称有限差分法的显示差分格式；②当θ=1/2时称有限差分法的对称（中心）差分格式；③θ=1称有限差分法的隐式差分格式。

有限差分方程实际上是基本微分方程的近似表达式，其近似程度可用泰勒级数进行分析。通过微分方程的差分表达式，可以看出在利用差分格式代替微分式时，是存在误差的，即用有限差分方程组模拟地下水流系统会产生误差。

（三）对于边界条件和垂向水量交换的处理

不论是已知水头的一类边界或已知流量的二类边界，计算点落在边界上，该点就不需要列入均衡方程。垂向水量交换的处理也是如此，若点与抽水井重合，该点已列入均衡离散方程时，抽水量就直接参与该点所在均衡区的水均衡。

（四）均衡离散方程的解算

显然，在含水层参数和边界条件都给定的条件下，只要知道某时刻流场中所有点的水头值，就可计算出下个时间步长的所有点的水头。即在已知初始条件的基础上，可以计算不同时刻各点水头值、不同时刻的流场。对于这类问题的求解方法，从广泛使用微机处理的角度来看，超松弛迭代为许多研究者所采用。

（五）应用

综上所述，在已知初始条件、边界条件、垂向水量交换以及给定含水层参数的情况下，可计算渗流区内不同时刻、不同节点的水头值。当前，不论是在地下水资源评价的水量计算中，还是在矿山开采地下水的疏干计算或在因大面积地下水位下降引起的地质灾害防治中，数值法都得到了广泛应用。

目前有许多地下水数值法计算软件，适应性强、有较高的仿真性，广为采用，例如，MOP-FLOW（孔隙水三维有限差分法数值计算软件），GWMS-3D（二维或三维地下水流和污染物质运移数值模拟软件）等。

（六）实例

通过实例的学习，使同学们对用数值法求解过程有所了解。这个过程包括：①水文地质条件概化，建立概念模型；②根据水文地质概念模型，建立数值模型；③剖分计算区，整理计算资料；④校正数值模型；⑤验证数值模型；⑥运用模型进行预报。

实例位于太行山东麓冲洪积扇的交界处。含水层为第四纪松散层，上部为细砂和粉砂层，下部为砂卵砾石、粗砂砾石加土层、含粘土砾石层等。上部含水层地下水已被疏干，当前开采层埋深为40～80m，水位埋深多在10m以下，漏斗中心区已达30m。边沿部分地区水位埋深为2～10m。

1.水文地质概念模型

①含水层底板为隔水粘土层；②含水层主要为非均质各向同性的潜水含水层；③计算区的边界三面为已知水头的一类边界，另一面为不同程度的弱透水层，计算区面积近600km；④区内有开采井；⑤地下水流为非稳定平面流，水流符合达西流。

2.数值模型

1）微分方程：

地下水动力学

2）初始条件：H（x，y，0）=H₀（x，y）

3）一类边界条件：H（x，y，t）｜_Γ1=H₁（x，y，t）

4）二类边界条件：

地下水动力学

式中：W为汇源项，由降水入渗量和井的开采量代数和求出；n为内法线；其他符号同前。

3.剖分计算区并整理计算资料

将计算区剖分为506个小区、230个节点，其中第一类边界点40个，二类14个，取旱季为模型校正时段，给出10个分区参数并经过试验给出参数初值。

4.校正数值模型

校正结果表明，微分方程和边界条件吻合。

5.验证数值模型

取雨季水位资料，分7个时段进行水位验证。根据验证资料绘制高低水位拟合图以及其他所需拟合图件，证明拟合程度良好，符合规范要求。

6.模型使用

利用验证过的符合实际的模拟模型，根据设计水位预计开采量，或根据设计的开采量预计不同时段的水位降低，尤其是漏斗中心的水位降低。

数值模拟流程

不同的软件进行数值模拟时所需的参数、计算方法、剖分格式等不尽相同，数值模拟的过程也不同，但大致相同，本文以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为例）

（四）概念模型的建立

利用所掌握研究区的地形地貌、气象水文、地层岩性、地质构造、水文地质、水文地球化学等资料，根据研究目的初步建立地质模型雏形。为了解决实际问题，往往要对所建模型及初始条件、边界条件等进行必须要的概化，概化内容包括：

1）含水层系统结构的概化：应依据含水层的类型、岩性、厚度、渗透系数等，将内部结构概化为均质、非均质各向同性或者异性的含水层；

2）侧向边界和顶，底部边界的概化；根据研究区边界的划定，将侧向边界和顶、底部边界概化为一、二和三类边界条件；

3）根据所要解决问题的难易程度，将地质模型概化为一维、二维或者三维模型，并进行合理精密级的网格剖分；

4）源汇项的处理：储存库一般位于地表800m以下，上覆有致密的盖层岩石，CO₂难以通过越流方式得到补充和排泄，在理想的条件下整个地质储库系统的CO₂通过径向对流或抽取井的方式得到补充和排泄。

（五）模拟方案的设计

针对不同的研究问题，可以根据设计者技术水平、场地经验等设计出不同的方案。利用数值模拟技术对多种方案进行模拟分析，评估方案实施的可行性，并对方案进行优化，最终得出经济、合理的方案。如进行CO₂灌注场地的遴选时，要对几个目标场区进行灌注能力，储存潜力、CO₂在储层中的扩散速度及距离以及后期CO₂泄漏风险等环节进行模拟，然后对模拟输出的数据进行后期处理。通过分析、统计这些数据确定最佳的灌注场地，同时合理划分储盖层的层位、厚度等。对场区灌注能力和储存潜力评估是数值模拟技术所能解决的CO₂地质储层领域的又一难题。不同的灌注方式导致CO₂进入储层的速度、流量等不同。根据设计者所设计的所有可行的注入方式进行模拟，分析比较不同方案下注入能力和注入量的大小，确定最优的灌注方式。

模拟方案的设计要根据研究问题不同视情况而定，而且方案设计的合理与否依赖于设计者自身理论和实际场地经验，不同的研究着对于同一个问题可能设计出不同的方案。我们可以针对这些方案建立不同的模型，通过模拟技术判断其合理性与可行性，最后确定最佳的方案。

（六）数值模型及模拟软件的选择

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

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

由于实际应用时多采用已有的数值模拟软件对CO₂地质储存的全过程进行模拟.不涉及软件的开发及程序代码的编写，只需根据研究的需要选择合适的软件进行模拟预测，而软件一旦选定，数学模型和数值模型基本上已经确定。以TOUGHREACT为例，基于上述所建的概念模型建立数学模型。气相和液相统一的偏微分方程为（9-1），咸水的偏微分方程为（9-2），超临界方程为（9-3）。方程中所涉及的相关字符和角标见表9-2。

二氧化碳地质储存技术方法概论

二氧化碳地质储存技术方法概论

表9-2 数学模型中所涉及的符号的意义

（七）数值模型的建立

1.网格剖分

建立了地质模型以后要对研究区离散化，即进行网格的剖分。首先确定离散点，即把所研究的区域按照某种几何形状（如矩形、任意多边形等）分割成网格系统。研究区的边界可以用最接近它的格线近似表示。当网格划分得足够小时，曲折的格线也能够很好地刻画出边界的形状。此过程又称研究区域的离散化（剖分）。在离散化时要遵循以下两条基本原则。

1）几何相似。要求物理模拟模型从几何形状方面接近真实被模拟体。

2）物理相似。要求离散单元的特性从物理性质方面（含水层结构、水流状态）近似于真实结构在这个区域的物理性质。

网格类型大体上可分为规则和不规则两种。规则网格剖分包括矩形、三角形和其他规则图形的剖分（图9-3），而不规则网格则包括不规则多边形等。网格的形状主要根据研究区形状而定。

图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₂地质储存工程实施及运行过程中可能出现的隐患。

如何做好数值模拟

数值模拟是一种数学模型，通过对某个系统或过程进行建模和计算机模拟，从而得出该系统或过程的一些特定结果。在做好数值模拟时，可以按照以下步骤进行：

1.确定模拟目标：首先需要确定模拟的目标，即想要得到哪些结果或解决哪些问题。确定模拟目标后，可以选择适合的数值模型进行建模。

2.收集数据和确定参数：在建立数值模型前，需要收集与该系统或过程相关的数据，并且需要确定参数。数据可以从实验中获取，也可以从已有的文献中收集，参数的选取需要经过合理的分析和推断。

3.建立数值模型：根据实际情况，选择适合的数值模型进行建立，比如基于微分方程、偏微分方程、最小二乘法等。在建立数值模型时需要保证模型的准确性和可靠性。

4.数值计算：通过计算机进行数值计算，将建立好的数值模型进行计算。计算需要注意参数选取的合理性和计算精度的控制，同时需要进行计算结果的验证和比对，以保证计算结果的可信度。

5.分析和解释结果：对计算结果进行分析和解释，找出计算结果中的规律和特点，以及对结果进行合理的解释和解读。

6.结果应用和优化：根据结果进行应用和优化，比如优化系统设计、改进工艺流程、提高产品性能等。

总之，做好数值模拟需要在建立数值模型、计算结果、结果分析和解释等方面进行科学合理的设计和操作，从而得出有价值的数值模拟结果。

试说明数值模拟方法的特点，它与理论研究，实验研究有什么关系

一、数值模拟的特点：

1、借助于计算机实现。在计算机上实现一个特定的计算，非常类似于履行一个物理实验。这时分析人员已跳出了数学方程的圈子来对待物理现象的发生，就像做一次物理实验。

2、运算高效且精准。因为数值模拟利用的是计算机系统中稳定的计算。

二、数值模拟与理论研究、实验研究的关系：

1、数值模拟是以理论研究为基础进行的模拟研究，没有科学的理论支持，无法做到高效的数值模拟。

2、数值模拟类似于实验研究，是利用计算机的实验研究。即利用计算机来进行实验。

扩展资料

数值模拟的步骤为：

1、首先要建立反映问题（工程问题、物理问题等）本质的数学模型。具体说就是要建立反映问题各量之间的微分方程及相应的定解条件。这是数值模拟的出发点。

2、数学模型建立之后，需要解决的问题是寻求高效率、高准确度的计算方法。计算方法不仅包括微分方程的离散化方法及求解方法，还包括贴体坐标的建立，边界条件的处理等。

3、在确定了计算方法和坐标系后，可以开始编制程序和进行计算。实践表明这一部分工作是整个工作的主体，占绝大部分时间。

4、在计算工作完成后，大量数据只能通过图像形象地显示出来。因此数值的图像显示也是一项十分重要的工作。利用录像机或电影放映机可以显示动态过程，模拟的水平越来越高，越来越逼真。

数值模拟方法

如前所述，基于镜质组化学反应动力学的煤热演化史数值模拟方法经历了简单函数关系模式、受热时间－经验法模式、反应活化能－温度函数模式和平行反应化学动力学模式四个发展阶段（秦勇等，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＝∑γ+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）。

什么是数值模拟方法与实验的优缺点呢

对有条件进行实验的材料，尽量采用实验方法，辅以数值模拟检验。而在工　　程应用中，很多情况下无法进行实验，如采矿问题等，数值模拟内部程序有相应的计算方法，能模拟较复杂过程。　　直观性与求解速度：实验直观性强，数值模拟直观性不如实验方法好，较抽象，但可以　　快速得到结果。实验操作复杂。　　成本：实验成本高，数值模拟成本低廉，只需在计算机上进行模拟和数据处理。　　施加载荷：数值模拟可以任意施加各种方向的载荷，可以施加实验方法达不到的条件。　　因此数值模拟方法在监测、设备开发、优化、效果预测方面体现了重要价值。　　数据采集：实验只能采集到特定点的的应力应变等数据，不能得到整个材料各点的应力　　应变值，而数值模拟方法可以对各个区域、各个测点进行应力分析和位移分析，对实验进行补充。　　数据处理：应将实验方法和数值模拟方法结合起来使用，分别对结果进行分析后，充分　　考虑两种方法各自的优缺点，互相比较印证，结合理论分析，有针对性地进行数据和结果的修正，才能得到一个比较全面、客观的结论。　　结果可靠性：数值模拟方法在模拟分析过程中，往往要对边界条件和材料属性进行简化，　　或多或少对分析结果产生影响，而且结构离散化的形式不同，得到的结果和精度也不同，随机性比较大，可信度降低。而在实验中不可避免的客观、主观因素也会产生误差，但是比数值模拟的误差少得多，可靠性更高。　　两种方法互相检验：合理的数值模拟方法对实验研究和理论分析具有指导作用，可以弥　　补实验工作的不足。实验与数值模拟结果比较，用来判断数值模拟方法的可行性。