数值模拟计算方法

❶ 数值模拟主要过程和步骤

1、首先要建立反映问题（工程问题、物理问题等）本质的数学模型。

具体说就是要建立反映问题各量之间的微分方程及相应的定解条件。这是数值模拟的出发点。没有正确完善的数学模型，数值模拟就无从谈起。牛顿型流体流动的数学模型就是着名的纳维—斯托克斯方程（简称方程）及其相应的定解条件。

2、寻求高效率、高准确度的计算方法

由于人们的努力，目前已发展了许多数值计算方法。计算方法不仅包括微分方程的离散化方法及求解方法，还包括贴体坐标的建立，边界条件的处理等。这些过去被人们忽略或回避的问题，现在受到越来越多的重视和研究。

3、开始编制程序和进行计算

实践表明这一部分工作是整个工作的主体，占绝大部分时间。由于求解的问题比较复杂，比如方程就是一个非线性的十分复杂的方程，它的数值求解方法在理论上不够完善，所以需要通过实验来加以验证。正是在这个意义上讲，数值模拟又叫数值试验。应该指出这部分工作决不是轻而易举的。

(1)数值模拟计算方法扩展阅读：

数值模拟的发展史：

1955年Peaceman与Rachford研发的交替隐式解法(ADI)是数值模拟技术的重大突破。该解法非常稳定,而且速度快，所以迅速在包括石油,核物理，热传导等领域得到广泛应用。1958年Douglas,Jim和Blair,P.M第一次进行了考虑毛管压力效果的水驱模拟。

60年代数值模拟技术的发展主要在数值解法，第一个有效的数值模拟解法器是1968年Stone推出的SIP(Strong Implicit Procere)。该解法可以很好地用来模拟非均质油藏和形状不规则油藏。

Stone在70年代发表了三相相对渗透率模型，由油水和油气两相相对渗透率计算油、气、水三相流动时的相对渗透率，该技术现在还广为应用。70年代另一项主要成就是Peaceman提出的从网格压力来确定井底流压的校正方法。

参考资料来源：网络—数值模拟

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

一、数值模拟的特点：

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

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

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

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

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

(2)数值模拟计算方法扩展阅读

数值模拟的步骤为：

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

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

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

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

❸ 数值方程与数值模拟

常用的数值计算方法有有限差分法和有限单元法。由于有限单元法中的集中储量有限元方法较通常的有限元法具有更多的优点，而且在边界条件的处理上，集中储量有限元法比有限差分法更符合实际，它考虑了边界节点的均衡单元的储水量变化（吴金全，1989）。

图1.4.3 均衡区域示意图

（一）集中贮量有限元公式推导

取单位水平面积、高度为计算层厚度的土柱进行研究（图1.4.3），将土柱（计算区域）垂直向上剖分为n个单元，空间步长为Δz，节点编号为0，1，2，…，n-1，n，在Δt时段内（Δt=t^j+1>-t^j），对任一内节点i所代表的均衡区z_i-1/2到z_i+1/2（图1.4.3）之间的土体列水量均衡方程（暂先不考虑极系吸水项）。

1.内节点（i=1，2，…，n-1）

由达西定律

（h 方程）得：

（1）通过z_i-1/2断面的水流通量（流入量）为：

土壤水盐运移数值模拟

（2）通过z_i+1/2断面的水流通量（流出量）为：

土壤水盐运移数值模拟

（3）均衡区域Δz内储水量的变化量（增量）为：

土壤水盐运移数值模拟

根据质量守恒原理（流入量-流出量=储存量的变化量）得：

Δq_i=q_i-1/2-q_i+1/2 （1.4.19）

将式（1.4.16）、式（1.4.17）、式（1.4.18）、式（1.4.19）代入得：

土壤水盐运移数值模拟

式（1.4.20）中，负压h及参数C和K在时间上取时段末j+1时刻的值，并整理得：

土壤水盐运移数值模拟

与有限差分方程比较，集中储量有限元推导出的有限元方程式（1.4.21）与隐式差分方程（h方程）是完全一致的。因此，具有无条件稳定和收敛的优良特性，故选用隐式差分格式对数学模型进行数值离散。若在时间上取时段中间j+1/2时刻的负压h及参数C、K，则可得出与Crank-Nicolsen差分格式完全一致的方程。

若考虑源汇项根系吸水项S，则式（1.4.22）变为：

土壤水盐运移数值模拟

令：

，

，r₃=Δt得：

土壤水盐运移数值模拟

式中：i=1，2，…，n-1。

令

土壤水盐运移数值模拟

将式（1.4.24）代入式（1.4.23）得：

土壤水盐运移数值模拟

2.边界节点的处理

（1）上边界节点i=0处的方程为：

土壤水盐运移数值模拟

式中：

为

，R^j+1/2，为时段平均入渗强度，E^j+1/2为时段平均表土蒸发强度。

令（1.4.26）式中：

土壤水盐运移数值模拟

则（1.4.27）式变为：

土壤水盐运移数值模拟

（2）下边界节点i=n为第一类边界节点，h_n已知，故不需列方程计算，这样第n-1个方程可简化为：

土壤水盐运移数值模拟

式中：

土壤水盐运移数值模拟

3.方程组

综合内节点和边界节点方程，从而得如下代数方程组：

土壤水盐运移数值模拟

方程组式（1.4.31）中：b₀、c₀、f₀按式（1.4.27）式计算，f_n-1按式（1.4.30）式计算，其余α_i、b_i、c_i按式（1.4.24）计算。

方程组式（1.4.31）用矩阵表示可简化为：

［A］［H］^j+1=［F］ （1.4.32）

式中：［ A］为系数矩阵；［ F］为常数项列阵；［ H］^j+1为求解未知量的列阵。

这样，通过数值方法将描述土壤水分运动的偏微分方程转化为求解代数方程组的问题。方程组式（1.4.31）系数矩阵元素满足α_ij=0（当｜i-j｜＞1 时），为三对角方程组，所以，可用“追赶法”求解。

（二）方程的线性化与土壤水分运动参数的取值

系数矩阵［A］中的各元素由时段末（j+1）时刻的土壤水分运动参数给出，常数项列阵［F］中的元素除含有已知时段初j时刻的负压h外，还含有时段末（j+1）时刻的土壤水分运动参数。然而土壤水分运动参数本身又是负压h的函数，因而求解方程组原则上说是非线性的。在利用数值方法求解土壤水分运动方程时，必须将方程线性化，使求解方程组成为线性代数方程组。

因迭代法计算的误差可以控制，求得的结果较逼近实际，而且一般可允许选用较大的时间步长（雷志栋等，1988），故选用迭代法进行线性化。

首先取时段初的参数如

作为时段末参数

的预报值，然后解方程组［ A］［H］^j+1=［F］，求得时段末各节点负压h的第一次迭代值

，根据

及K-h曲线可求得土壤水分运动参数

的校正值。以此参数的校正值作为下一次计算的预报值，然后解方程组可得时段末各节点负压h的第二次迭代值

。重复上述步骤，直到前后两次迭代值，第p-1次和第p次迭代值满足下式为止：

土壤水盐运移数值模拟

式中：e迭代误差为任意给定得正的小数，一般取e=0.01。

参数的取值，一般的说，用三点式或几何平均的方法效果较好（雷志栋等，1988），计算也不复杂，这里选用几何平均的方法：

土壤水盐运移数值模拟

同理，根据达西定律（

）（θ方程）可以推导出集中储量有限元公式的θ方程，该方程与隐式差分方程θ方程完全一致，若在时间上取时段中间j+1/2时刻的含水率θ及参数D、K，则可得到与Crank-Nicolsen差分格式（中心差分）完全一致的方程。θ方程与h方程的求解方法完全一样，由于本文所研究的是双层结构问题，而在分界面处θ不连续、h连续，所以选用h方程进行计算，因此对于θ方程这里就不做推导了。

（三）数值模拟

1.模型验证

进行数值模拟，首先进行模型验证。模型验证时，上边界条件表达式中的θ₁₀由实际观测资料给出。根据有作物生长条件下土壤水分运动的基本方程和差分方程，在已知初始条件和边界条件时，模型验证可以通过以下步骤进行：①根据实测初始负压剖面的分布，用三次样条插值给出各节点上的初始值；②计算蒸发量E；③计算根系吸水层厚度L_r及吸水率S；④根据差分方程计算时段末负压值h。

模型验证时，以计算起始时刻的实测负压剖面（或含水率剖面）作为初始剖面，空间步长选用1cm，根据最底部负压计测点和中子仪测点，剖面深度为120cm（含水率剖面为130cm），时间步长选用1h，迭代相对误差e≤0.01。计算中输入的大量信息，如各节点的初始负压值、降雨量、降雨日期、水面蒸发强度、根层土壤含水率等均以数据文件的形式提供。由于大田盖200kg/亩和盖600kg/亩只进行了含水率观测，计算时先将含水率θ转化为负压h，计算结束后再将负压h转化为θ。大田盖400kg/亩有负压资料，可直接用负压h计算。上边界条件由E/E₀-θ关系给出。数值模拟按覆盖量（200、400、600kg/亩），分生育阶段（400kg/亩）进行，模拟计算在PC机上完成。主要模拟的试验处理有；拟合曲线见图1.4.4。

（1）I-2盖200kg/亩，模拟时间为7月31日至8月30日，共31天。

（2）I-3盖400kg/亩，模拟时间为苗期：6月25日至7月30日，共26天；拔节：7月21日至8月10，共21天；灌浆-成熟：8月11日至9月17日，共38天。

（3）I-4盖600kg/亩，模拟时间为7月31日至8月30日，共31天。

2.模型验证结果及讨论

根据描述土壤水分运动的定解问题，通过数值模拟可以得到土壤水分运动的动态过程，并用实测结果对模型进行验证。如果数学模型能够描述实际的物理过程，排除随机因素外，模拟得到的土壤水分动态过程（模拟值）与实际观测得到的土壤水分动态过程（实测值）应该完全吻合。比较图1.4.4，从图中可以看出，模拟值与实测值吻合较好，表明本文提出的考虑秸秆覆盖有作物生长条件下的模型是可靠的，以上讨论的数值方法是可行的。不同覆盖量、不同生育阶段，可以用不同的E/E₀-θ经验公式来反映覆盖对水分运移的影响。因此，本文提出的模型和数值方法可以用来模拟秸秆覆盖条件下田间土壤水分的运动，可对田间土壤水分动态作中短期预报。

图1.4.4 实测值与模拟值对比图

3.模型的应用——预报

数值模拟的目的之一就是进行预报，根据气象部门提供的降雨量及水面蒸发强度等气象资料，使用验证过的模型进行田间土壤水分动态预报。本文使用实际发生过的降雨量及水面蒸发强度系列资料进行预报，用实测负压资料检验预报结果。程序的运行见图1.4.5。检验的实测资料选用大田覆盖400kg/亩的资料，分生育阶段（苗期、拔节、灌浆-成熟）进行。从图1.4.6可以看出，预报值和实测值吻合较好。

图1.4.5 双层结构有根系吸水项垂向一维土壤水数值模拟框图

图1.4.6 预报值与实测值对比图（大田盖400kg/亩）

❹ 数值模拟方法

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

❺ 数值模拟流程

不同的软件进行数值模拟时所需的参数、计算方法、剖分格式等不尽相同，数值模拟的过程也不同，但大致相同，本文以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₂地质储存工程实施及运行过程中可能出现的隐患。

❻ 湍流的数值模拟方法

现有的湍流数值模拟方法有三种：直接数值模拟、雷洛平均模拟和大涡模拟。

计算复杂程度：雷洛平均模拟-大涡模拟-直接数值模拟

湍流的运动可以看做是时均运动与随机运动的叠加，按照雷洛时均法，运动中的变量为 dt

不可压缩时均运动控制方程组之所以出现方程组不封闭，3个速度，物理量，压力，6个雷洛应力（需求解的位置函数较方程数多），在于方程中出现了湍流脉动值的雷洛应力项，要使得方程组封闭，必须对雷洛应力做出某些假设，即建立应力表达式（或者引入新的湍流方程），通过此表达式把湍流的脉动值与时均值等联系起来。基于某些假设所得出的湍流控制方程，称为湍流模型。

基于这些假设有（皮轮世应力：压应力与切应力）

1、雷洛应力模型：构建雷洛应力补充方程；

2、湍流黏度类模型：引入 湍动黏度 或涡旋黏度，然后湍流应力成为湍动黏度的函数，整个计算在于确定湍动黏度；

3、零方程模型：不使用微分方程，采用代数关系，将湍动黏度与时均值联系起来；

4、一方程模型：建立了湍动能的运输方程；

5、两方程模型：补充两个微分方程使湍流时均控制方程组封闭的一类处理方法。

湍流 两方程模型

湍流动能耗散率(turbulent dissipation)是指在分子粘性作用下由湍流动能转化为分子热运动动能的速率。通常以单位质量流体在单位时间内损耗的湍流动能来衡量，以ε表示。湍流速度在空间上存在着随机涨落，从而形成了显着的速度梯度，在分子粘性力作用下通过内摩擦不断地将湍流动能转化为分子运动的动能。大气湍流的动能耗散主要发生在大小为毫米数量级的湍涡。(%)

湍流动能是湍流速度涨落方差与流体质量乘积的1/2。（J）TKE/M=1/2(ui'^2)

K 越大表明湍流脉动长度和时间尺度越大， ε 越大意味着湍流脉动长度和时间尺度越小，它们是两个量制约着湍流脉动。

标准 模型（standard model）： 分别建立了 与 的运输方程， 为湍动能， 为湍动能耗散率，还是用来求解燃肢 。

适应性：（1）模型中的相关系数，主要根据一些特定条件下的试验结果来确定的；（2）给出的 模型是针对湍流发展非常充分的湍流运动来建立的，即针对 高Re湍流模型 。而当Re较低时（例如近壁面区的流动），湍流发展不充分，湍流的脉动影响可能不如分子黏性影响大，在近壁面可能再现层流。常用的解决壁面流动方法有：一种是壁面函数法；一种是采用低Re的 模型。；（3）标准 用于 强旋流、绕弯曲壁面流动或弯曲流线运动时，会产生一定的失真 。在标准 模型中，对于雷洛应力的各个分量，假定湍流黏度是相同的，即是各向同性的标准，但是在弯曲流线的情况下，湍流时各向异性的。

RNG   模型：通过在大尺度运动项和修正后的黏度项中体现小尺度的影响，而使这些小尺度运动系统地从控制方程中去除。建立的 方程与标准 模型中的方程很相似，与标准的 模型的主要变化有：（1）通过修正的湍动黏度，考虑了平均流动中的旋转及旋转流动情况；（2）在 方程中的产生项增加了一项，从而反映了主流时均应变率Eij。这样，RNG  模型产生项不仅仅与流动情况相关，而且在同一问题中还是空间坐标的函数。 从而可以更好的处理高应变率及流线弯曲程度较大的流动 。RNG  模型仍然是针对充分发展的湍流，而对近壁面区（低Re数）的流动与标准 模型有着相同的困扰。

Realizable  模型：标准模型在应对时均应变率特别大的时候，有可能导致负的正应力，为了使流动符合湍流的物理定律，需要对正应力进行某种数学上的约束，为了保证这种约束的实现，湍流黏度计算式中的系数 应与应变率联系起来。提出Realizable  模型。与桐备标准的 模型的主要变化有：（1）湍动黏度计算公式发生变化，引入了旋转和曲率有关的内容；（2） 方程发生大变化；（3） 旋转均匀剪切流、包含有射流和混合流的自由流动、管道内流动、边界层流动，以及带有分离的流动 。

近壁区使用 模型的问题及对策

有固体壁面的充分发展的湍流流动，沿壁面法线的不同距离上，可将流动分为壁面区和核心区两个部分。核心区为完全湍流区。壁面区分为3个子层：（1）黏性底层；（2）过渡层；（3）对数律层。以黏性力与雷洛切应力的相对大小划分的。

壁面函数法：将壁面上的物理量与湍流核心区域内待求解的未知量直接联系。它必须和高Re的 模型配合使用。基本思想为对湍流核心区的流动用 模型求解，而在壁面区不进行求解，直接使用半经验公式将壁面上的物理量与湍流核心区内求解变量联系起来。划分网格时不需要对壁面区加密，只需把第一个内节点布置在对数律成立的区域内，即配置到湍流充分发展的区域。也可以对壁面区网格加密，以得到近壁区物理量分布。

低Re的 模型：为了使得数值计算从低Re区域一直进行到固体壁面上（该处Re=0），主要的方法是在输送方程中考虑：壁面的黏性、流态的不同、湍流动能的耗散不是各向同性。在使用低Re的 模型进行流动计算时，充分发展的湍流核心区及黏性底层均用同一套公式，这一套公式中当Re数较大时，关于部分修正变成1，且由于黏性底层的速度梯度大因而黏性底层的网格密。 当局部湍流的Re数小于150时，就应该使用低Re的 模型 。

雷洛应力模型

两方程模型难以考虑旋转流动及流线曲率变化的影响，直接对Reynolds方程中湍流脉动应力直接建立微分方程，并进行求解，称之为雷洛应力模型。

RSM由连续性方程、动量守恒方程、应力方程、 方程和 方程构成12个方程组构成封闭的三维湍流流动问题的基本控制方程组，可以通过SIMPLE等算法求解。如果需要对能量或组分等进行计算，需要建立针对标量型变量（如温度、组分、浓度）的脉动量的控制方程。

RSM也是针对湍流发展非常充分的湍流运动来建立的，当Re较低时，上述方程不再适用。在近壁区需要采用壁面函数法或低Re的RSM模型。尽管RSM比 模型应用范围广，包含更多的物理机理，但是它有很多缺陷。计算实践表明，RSM虽然能考虑一些各向异性效应，但是不一定比其他模型好，在计算突扩流动分离区和计算湍流运输各向异性较强的流动时，RSM优于两方程模型，但是RSM计算量大于两方程模型。

大涡模拟

大涡模拟的基本假设是：（1）动量、能量、质量及其他标量主要由大涡输运；（2）流动的几何和边界条件决定大涡的特性，而流动特性主要在大涡中体现；（3）小尺度涡旋受几何和边界条件影响较小，并且各向同性。大涡模拟过程直接求解大涡，小尺度涡旋用模型来封闭。

LES的控制方程是对N-S方程在波数空间或者物理空间进行过滤得到的。过滤的过程是去掉比过滤宽度或者给定物理量宽度小的涡旋，从而得到大涡漩的控制方程，在变量的方程中加入过滤函数。FLUENT中，大涡模拟只能针对不可压缩流体的流动。

高雷洛数流动，两者基本差不多。但是在低雷洛区，重整化群理论的压网格模型对流动转捩和近壁流动问题有较好的模拟效果。

转捩，即从层流到湍流的过渡。流体力学名词，表征一种流动现象，英文为transition。转捩点的计算和预估是设计飞行器的关键前提。

❼ 数值模拟的方法有哪些，各自有什么优缺点，谢谢！

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

❽ 有限元数值模拟方法

有限单元法是应用于构造应力场模拟的最广泛的数学模拟方法。其基本思想是将所研究的地质体以一定的方式（单元形状和节点个数）简化为有限个单元组成的离散化模型，再用相应的计算程序求出数值解答。利用有限元法数值模拟，可以利用地质调查和构造解析获得的较少地质应力状态的资料来反演区域内各点的应力状态，从而获得区域的构造应力场特征，加深认识区域内的构造演化。目前有限单元法的应用已由弹性力学的平面问题扩展到空间问题、板壳问题，分析对象从弹性材料扩展到塑性、粘弹性、粘塑性和复合材料。

有限元法数值模拟随着计算机技术的发展在科学计算领域得到广泛应用，20世纪80年代以来，国际上已有较大型的有限元计算程序达几百种，其中较着名的有：ANSYS、NASTRAN、ASK、ADINA、SAP等。以ANSYS为代表的数值模拟软件将有限元分析、计算机图形学和优化技术相结合，已成为科学计算领域不可缺少的有力工具。

基于本区岩石圈的三维结构特点，我们首先对本区的三层结构相互作用关系进行了模拟。对本区的物理模拟研究，前人已经做过很多工作，其中在对印度板块挤压下亚洲中东部的构造模拟中，有的反映出大型走滑断裂、裂谷和张性盆地以及压性逆冲断裂等构造现象，有的反映出多层构造中网络状流动现象，认为板内变形受塑性流动网络控制（Tapponnier et al.，1982；李建国等，1997）。这些工作往往只反映了本区的某一方面的特性，而无法对本区的构造形态做出动力学的完善解释。因此在前人的工作基础上，我们首先建立了本区的一个三层结构模型，其中中上地壳深度根据天然地震资料定为30 km，下地壳以莫霍面为其底界，根据地震测深资料取50 km。因为本模型建立的主要目的是确定岩石圈各圈层之间的作用关系，因此模型底部只考虑到100 km的深度。

❾ 数值模拟

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

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

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

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.模型使用

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

❿ 应力场数值模拟方法

近30年来，人们采用现场测试、实验室试验、理论分析与模型试验等多种方法，使岩土力学研究取得很大进展^{［162～166］}。如今随着计算机技术的快速发展，岩土力学的研究进入了一个新的阶段，其中数值计算方法已成为解决岩土力学问题的重要手段之一。

6.1.1 概述

许多工程分析问题，如固体力学中的位移场和应力场分布分析、电磁学中的电磁场分析、振动特性分析、传热学中的温度场分析以及流体力学中的流场分布等，都可以通过在给定边界条件下对其控制方程进行求解得到，但是利用解析方法只能求出一些方程性质比较简单且几何边界相当规则的极少数问题。对于大多数实际工程技术问题，由于物体的几何形状比较复杂或者问题的某些特性是非线性的，因而一般无解析解。为了解决此类问题，一般采用两种处理方法：一种是进行简化处理，将方程和边界条件简化为能够处理的问题，从而得到在简化情况下的解，但这种方法应用非常有限，且假设过多将会导致错误的解；另一种是在广泛接收现代数学和力学理论的基础上，借助于计算机和计算软件来获得工程上要求的数值解，这就是目前应用非常广泛的数值模拟方法。

目前在工程技术领域内常用的数值分析方法包括：有限单元法、边界元法、离散单元法以及有限差分法。最初常用的是有限差分法，它可以处理一些相当复杂的问题。但对于几何形状复杂的边界条件，其解的精度受到影响。20世纪60年代出现并得到广泛应用的有限单元法，使经典力学解析方法难以解决的工程力学问题都可以用有限元方法求解。它将连续的求解域离散为一组有限个单元的组合体，解析地模拟或逼近求解区域。由于单元能按各种不同的联结方式组合在一起，且单元本身又可有不同的几何形状，所以能适应几何形状复杂的求解域。但有限单元法需要的存贮容量常非常巨大，甚至大得无法计算。由于相邻界面上只能位移协调，对于奇异性问题（应力出现间断）的处理比较麻烦，这是有限单元法的不足之处。70年代末期，出现了另一种重要的数值方法为边界元法。边界元方法是把求解区域的边界剖分为若干个单元，将求解简化为求单元结点上的函数值，通过求解一组线性代数方程实现求解积分方程。上述两种数值方法的主要区别在于，边界元法是“边界”方法，而有限元法是“区域”方法，它们都是针对连续介质，只能获得某一荷载或边界条件下的稳定解。对于具有明显塑性应变软化特性和剪切膨胀特性的岩体，无法对其大变形过程中所表现出来的几何非线性和物理非线性进行模拟，这就使得人们去寻求适合模拟节理岩体运动变形特性的有效数值方法。

1971年Cundall，P.A^［167］提出了一种不连续介质数值分析模型——离散单元法。该方法优点在于适用于模拟节理系统或离散颗粒组合体在准静态或动态条件下的变形过程。离散单元法的基本原理不同于基于最小总势能变分原理的有限单元法，也不同于基于Betti互等定理的边界单元法，而是建立在牛顿第二运动定律基础上。最初的离散元法是基于刚性体的假设，由于没有考虑岩块自身的变形，在模拟高应力状态或软弱、破碎岩体时，不能反映岩块自身变形的特征，使计算结果与实际情况产生较大出入。Maini，T.，Cundall，P.A.^{［168～169］}等人针对刚体单元没有考虑岩块自身变形的缺点，利用差分方法提出了考虑岩石自身变形的改进的离散单元法，编制了通用的离散元程序UDEC（Universal Discrete Element Code），将离散元推广到模拟岩体破碎和变形情况，推动了离散元的进一步发展。我国学者也相继开展这方面的研究，王泳嘉教授^［170］等将离散单元法应用于采矿工程方面的研究。

6.1.2 FLAC数值模拟方法

（1）概述

数值模拟技术通过计算机程序在工程中得到广泛的应用。一直到20世纪80年代初期，国际上较大型的面向工程的通用程序有：ANSYS、NASTRAN、FLAC、UNDEC、ASKS以及ADINA等程序。它们功能越来越完善，不仅包含多种条件下的有限元分析程序，而且带有功能强大的前、后处理程序。

连续介质快速拉格朗日差分法（Fast Lagrangian Analysis of Continua，简写FLAC）是近年来逐步成熟完善起来的一种新型数值分析方法。把拉格朗日法移植到固体力学中，即将所研究的区域划分为网格，节点相当于流体质点，然后按照时步用拉格朗日方法来研究网格节点的运动，这就是固体力学变形研究中的拉格朗日数值研究方法。

FLAC与基本离散元法相似，但它克服了离散元法的缺陷，吸取了有限元法适用于各种材料模型及边界条件的非规则区域连续问题解的优点。FLAC所采用的动态松弛法求解，不需要形成耗机时量较大的整体刚度矩阵，占用计算机内存少，利于在微机的工程问题。同时，FLAC还应用了节点位移连续的条件，可以对连续介质进行大变形分析。

（2）数学模型

显式有限差分法的基本方程主要包括：平衡方程、几何方程、物理方程和边界条件。在FLAC^3D2.0中采用的拉格朗日描述方程，一般规定介质中一点由向量分量x_i，u_i，v_i，dv_i/dt（i=1，2，3）来表征，其分别代表位置、位移、速度和加速度分量。

其基本原理和基本公式简单叙述如下：

空间导数的有限差分近似

三维FLAC方法中采用了混合离散方法，区域被划分为常应变六面体单元的集合体；而在计算过程中，又将每个六面体分为常应变四面体，变量均在四面体上进行计算，六面体单元的应力、应变取值为其四面体的体积加权平均。

如图6.1所示，所研究区域任一四面体，节点编号为1～4，规定与节点n相对的面为第n面，设定其内任一点的速度分量为v_i，则由高斯散度定理得

煤岩动力灾害力电耦合

式中：V——四面体体积，m³；S——四面体外表面，m²；n_j——外表面单位法向向量分量。

图6.1 四面体

对于常应变单元，n_j在每个面上为常量，因此通过上式积分可得

煤岩动力灾害力电耦合

式中上标f表示f面的变量值，对于为线性分布的速率分量，速度分量的平均值为

煤岩动力灾害力电耦合

式中上标l表示节点l的变量值。将（6.3）式代入（6.2）式可得

煤岩动力灾害力电耦合

经过变换可得节点速率计算公式：

煤岩动力灾害力电耦合

1）平衡方程（运动方程）

显式有限差分法采用的平衡方程就是人们熟知的牛顿第二运动定律，即

煤岩动力灾害力电耦合

式中：F_i——节点合力在i方向分力，N；m_i——节点质量，kg；a_i——节点加速度在i方向分量，m/s²。

作用于各个节点的合力：外力（集中力、均布力、重力等）和内力（单元变形引起的应力在单元节点上的分量）。节点质量是根据节点相邻单元的面积（体积）和密度，按照面积（体积）加权求出。

FLAC^3D以节点为计算对象，将力和质量均集中在节点上，然后通过运动方程在时域内进行求解。节点运动方程可以表示为如下形式：

煤岩动力灾害力电耦合

式中：（t）———t时刻l节点在i方向的不平衡力分量，可以由虚功原理导出；m^l———l节点的集中质量，在分析静态问题时，采用虚拟质量；而在分析动态问题时，则采用实际的集中质量。

将（6.7）式左端用中心差分来近似，则可得

煤岩动力灾害力电耦合

2）变形协调方程——几何方程

作为连续介质力学，变形体之间必须满足变形协调方程（几何方程），否则变形体就会出现分离或嵌入。变形协调方程反映了位移与应变间的关系，对于某一时步的单元应变增量可由下式确定：

煤岩动力灾害力电耦合

求出应变增量后，即可由本构方程得到应力增量，各时步的应力增量叠加即可得到总应力，在大变形时，还需根据本时步单元的转角对本时步前的总应力进行旋转修正，然后即可由虚功原理求出下一时步的节点不平衡力，进入下一时步的计算。

3）物理方程——本构关系

物理方程反映应力与应变之间的关系，在程序中通常被称为材料模式或材料模型。在FLAC^3D2.0中提供了10种基本材料模型，它们是：①Null；②Elastic，isotropic；③Elastic，transversely isotropic；④Druck-Prager plasticity；⑤Mohr-Coulomb plasticity；⑥Ubiquitous joint plasticity；⑦Strain-hardening／softening Mohr-Coulomb plasticity；⑧bilinear strain-hardening/softening ubiquitous-joint plasticity；⑨Modified Cam-clay plasticity 和⑩elastic，orthotropic。

本文进行应力场数值模拟时采用的是Mohr-Coulomb应变硬化软化破坏准则，在FLAC^3D2.0中，Mohr-Coulomb 模型的破坏准则以主应力σ₁，σ₂，σ₃来描述，相应的应变为三个主应变ε₁，ε₂，ε₃。根据Hooke定律，应力、应变增量具有如下表达形式：

煤岩动力灾害力电耦合

式中α₁，α₂为材料常数，可以由体积模量K和剪切模量G确定：

煤岩动力灾害力电耦合

不失一般性，令σ₁≥σ₂≥σ₃，摩尔—库仑准则为

其中：

煤岩动力灾害力电耦合

式中C，φ分别为煤岩的粘聚力和内摩擦角。

FLAC^3D2.0的Mohr-Coulomb 破坏准则如图6.2所示。

图6.2 FLAC^3D的Mohr-Coulomb 破坏准则

本着作中就是选用上述的Strain-hardening／softening Mohr-Coulomb plasticity模型，对单轴压缩煤岩以及矿山地下煤岩独巷掘进时围岩的变形破坏过程进行模拟。

4）阻尼力

对于静态问题，FLAC^3D2.0在式（6.7）的不平衡力中加入了非黏性阻尼，以使系统的振动逐渐衰减直至达到平衡状态（即不平衡力接近零），此时节点运动方程变为：

煤岩动力灾害力电耦合

式中阻尼力（t）由下式确定：

煤岩动力灾害力电耦合

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

煤岩动力灾害力电耦合

5）初始条件与边界条件

边界条件包括面积力、集中载荷等应力边界条件和位移边界条件。此外也可加载体力和初始应力。在编写程序代码时，一般所有的应力和节点速度初始化为零，然后指定初始化应力。集中载荷则加载在面节点上，位移边界条件则以运动方程形式施加到相应的边界节点上。

边界条件分为应力边界条件和位移边界条件，应力边界条件为：

煤岩动力灾害力电耦合

式中：F_i———作用于节点i上的力；——作用于边界上的应力；n_j———边界上的法线沿j方向的矢量大小；Δs———边界的长度。

若是位移边界条件，应将边界条件以运动方程的形式施加到相应的边界节点上。

FLAC^3D2.0^［171］与FLAC^2D3.3也是由美国Itasca Consulting Group Inc开发的三维显式有限差分法程序，它可以模拟岩土或其他材料的三维力学行为。FLAC^3D2.0的计算循环过程如图6.3所示。

图6.3 FLAC^3D2.0的计算循环

6.1.3 FLAC数值模拟方法在采矿工程中的应用^{［172～179］}

采矿过程中围岩活动规律及巷道围岩稳定性问题涉及岩体力学特性、围岩压力、支护围岩相互作用关系及巷道与工作面时空关系等一系列复杂力学问题。随着我国经济建设的高速发展，岩土工程稳定性分析问题日益突出，除采矿工程外，在水利、交通（铁道和公路）、高层建筑的地基等行业也都存在着大量的岩土力学数值计算分析问题。能否用计算机数值模拟分析采矿岩层控制问题和岩土工程问题已成为一个大学岩层控制技术和岩土力学学科水平高低的标志之一。

与ANSYS、ADINA相比，FLAC 和UDEC的最大特点是计算分析岩土工程中的物理不稳定问题，因而特别适用于岩土工程中几何和物理高度非线性问题的稳定性分析，如采场的采动影响规律，软岩巷道的大变形问题，采动后的地表沉陷，露天矿的边坡稳定，水坝的稳定性等问题。

从力学计算方法上讲其主要特点

1）可以直接计算非线性本构关系；

2）物理上的不稳定问题不会引起数值计算的不稳定；

3）开放式程序设计（FISH），用户可以根据需要自己设计程序；

4）既可以分析连续体问题（FLAC），也可以分析非连续体问题（UDEC）；

5）可以模拟分析很大的工程问题；

6）高度非线性问题不增加计算时间。

在采矿工程中，许多学者利用FLAC软件对采矿过程中围岩活动规律及巷道围岩稳定性问题涉及到岩体力学特性、围岩压力、支护围岩相互作用关系及巷道与工作面的时空关系等一系列复杂的力学问题进行了一系列的研究，取得了显着的效果。梅松华等以施工期监测结果为基础，在正交设计原理的基础上，选定反演参数与水平，采用二维显式差分法FLAC进行弹塑性位移反分析。朱建明等在分析FLAC有限差分程序的基础上，提出了变弹性模量方法模拟时间因素对巷道围岩稳定性影响的衰减曲线，为揭示巷道围岩变形机理和有效指导围岩支护提供了有效的分析方法。来兴平等探讨了岩石力学非线性计算软件FLAC^2D3.3在地下巷道离层破坏数值计算中的应用。康红普对回采巷道锚杆支护影响因素进行了FLAC分析，认为FLAC^2D3.3在分析几何非线性和大变形问题方面性能优越。

在煤岩动力灾害预测中，这些方法的优点

1）可以提前知道煤与瓦斯突出、冲击矿压等煤岩动力灾害防治的重点区域；

2）可以得到大范围内的空间信息；

3）可以提前预测预报煤岩动力灾害的危险性；

4）可以确定在采掘过程中，应力的分布状况和集中程度。

在煤岩动力灾害预测中，这些方法也具有以下缺点

1）对实际问题均进行了简化处理；

2）对于煤岩体的力学特性，如弹性模量、泊松比等力学参数，也进行了简化，没有考虑其局部非均质性和各向异性；

3）只能作为一种近似方法使用。