地下水数值模拟能解决什么问题

数值模拟（数值法）是对数学模型的一种近似解法，它仅能求出计算域内有限点某个时刻水头的近似值，这个值在实际应用中可以满足精度要求。数值法可以解决许多复杂水文地质条件下的渗流计算问题，应用十分广泛。如用于大中型水源地、地下水的补径排条件复杂、渗流区形状不规则、含水介质为非均质各向异性等条件下，确定水头分布和流量计算。（一）渗流区域离散化（以二维流为例）采用数值模拟技术研究地下水的运动，首先将要研究的水文地质模型内的含水层离散化。所谓离散化，就是将要研究的渗流区非均质各向异性含水层，按照一定的方式剖分（分割）成许多相互联系的小均衡区，在每个小均衡区内是均质各向同性的。在每个小的均衡区内，其含水层参数视为常数；其中心水头值或有条件下的平均水头值视为小均衡区内水头代表值。剖分通常采用两种形式（矩形、多边形）进行。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）均衡时段为Δtn+1：Δtn+1=tn+1-tn0 表示点（i，j）上tn时刻的水头。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.水文地质概念模型①含水层底板为隔水粘土层；②含水层主要为非均质各向同性的潜水含水层；③计算区的边界三面为已知水头的一类边界，另一面为不同程度的弱透水层，计算区面积近600km2；④区内有开采井；⑤地下水流为非稳定平面流，水流符合达西流。2.数值模型1）微分方程：地下水动力学2）初始条件：H（x，y，0）=H0（x，y）3）一类边界条件：H（x，y，t）｜Γ1=H1（x，y，t）4）二类边界条件：地下水动力学式中：W为汇源项，由降水入渗量和井的开采量代数和求出；n为内法线；其他符号同前。3.剖分计算区并整理计算资料将计算区剖分为506个小区、230个节点，其中第一类边界点40个，二类14个，取旱季为模型校正时段，给出10个分区参数并经过试验给出参数初值。4.校正数值模型校正结果表明，微分方程和边界条件吻合。5.验证数值模型取雨季水位资料，分7个时段进行水位验证。根据验证资料绘制高低水位拟合图以及其他所需拟合图件，证明拟合程度良好，符合规范要求。6.模型使用利用验证过的符合实际的模拟模型，根据设计水位预计开采量，或根据设计的开采量预计不同时段的水位降低，尤其是漏斗中心的水位降低。

数值模拟是一种数学模型，通过对某个系统或过程进行建模和计算机模拟，从而得出该系统或过程的一些特定结果。在做好数值模拟时，可以按照以下步骤进行：1.确定模拟目标：首先需要确定模拟的目标，即想要得到哪些结果或解决哪些问题。确定模拟目标后，可以选择适合的数值模型进行建模。2.收集数据和确定参数：在建立数值模型前，需要收集与该系统或过程相关的数据，并且需要确定参数。数据可以从实验中获取，也可以从已有的文献中收集，参数的选取需要经过合理的分析和推断。3.建立数值模型：根据实际情况，选择适合的数值模型进行建立，比如基于微分方程、偏微分方程、最小二乘法等。在建立数值模型时需要保证模型的准确性和可靠性。4.数值计算：通过计算机进行数值计算，将建立好的数值模型进行计算。计算需要注意参数选取的合理性和计算精度的控制，同时需要进行计算结果的验证和比对，以保证计算结果的可信度。5.分析和解释结果：对计算结果进行分析和解释，找出计算结果中的规律和特点，以及对结果进行合理的解释和解读。6.结果应用和优化：根据结果进行应用和优化，比如优化系统设计、改进工艺流程、提高产品性能等。总之，做好数值模拟需要在建立数值模型、计算结果、结果分析和解释等方面进行科学合理的设计和操作，从而得出有价值的数值模拟结果。

一、数值模拟的特点：1、借助于计算机实现。在计算机上实现一个特定的计算，非常类似于履行一个物理实验。这时分析人员已跳出了数学方程的圈子来对待物理现象的发生，就像做一次物理实验。2、运算高效且精准。因为数值模拟利用的是计算机系统中稳定的计算。二、数值模拟与理论研究、实验研究的关系：1、数值模拟是以理论研究为基础进行的模拟研究，没有科学的理论支持，无法做到高效的数值模拟。2、数值模拟类似于实验研究，是利用计算机的实验研究。即利用计算机来进行实验。扩展资料数值模拟的步骤为：1、首先要建立反映问题（工程问题、物理问题等）本质的数学模型。具体说就是要建立反映问题各量之间的微分方程及相应的定解条件。这是数值模拟的出发点。2、数学模型建立之后，需要解决的问题是寻求高效率、高准确度的计算方法。计算方法不仅包括微分方程的离散化方法及求解方法，还包括贴体坐标的建立，边界条件的处理等。3、在确定了计算方法和坐标系后，可以开始编制程序和进行计算。实践表明这一部分工作是整个工作的主体，占绝大部分时间。4、在计算工作完成后，大量数据只能通过图像形象地显示出来。因此数值的图像显示也是一项十分重要的工作。利用录像机或电影放映机可以显示动态过程，模拟的水平越来越高，越来越逼真。参考资料来源：百度百科-数值模拟

1、首先要建立反映问题（工程问题、物理问题等）本质的数学模型。具体说就是要建立反映问题各量之间的微分方程及相应的定解条件。这是数值模拟的出发点。没有正确完善的数学模型，数值模拟就无从谈起。牛顿型流体流动的数学模型就是著名的纳维—斯托克斯方程（简称方程）及其相应的定解条件。2、寻求高效率、高准确度的计算方法由于人们的努力，目前已发展了许多数值计算方法。计算方法不仅包括微分方程的离散化方法及求解方法，还包括贴体坐标的建立，边界条件的处理等。这些过去被人们忽略或回避的问题，现在受到越来越多的重视和研究。3、开始编制程序和进行计算实践表明这一部分工作是整个工作的主体，占绝大部分时间。由于求解的问题比较复杂，比如方程就是一个非线性的十分复杂的方程，它的数值求解方法在理论上不够完善，所以需要通过实验来加以验证。正是在这个意义上讲，数值模拟又叫数值试验。应该指出这部分工作决不是轻而易举的。扩展资料：数值模拟的发展史：1955年Peaceman与Rachford研发的交替隐式解法(ADI)是数值模拟技术的重大突破。该解法非常稳定,而且速度快，所以迅速在包括石油,核物理，热传导等领域得到广泛应用。1958年Douglas,Jim和Blair,P.M第一次进行了考虑毛管压力效果的水驱模拟。60年代数值模拟技术的发展主要在数值解法，第一个有效的数值模拟解法器是1968年Stone推出的SIP(StrongImplicitProcere)。该解法可以很好地用来模拟非均质油藏和形状不规则油藏。Stone在70年代发表了三相相对渗透率模型，由油水和油气两相相对渗透率计算油、气、水三相流动时的相对渗透率，该技术现在还广为应用。70年代另一项主要成就是Peaceman提出的从网格压力来确定井底流压的校正方法。参考资料来源：百度百科—数值模拟

不同的软件进行数值模拟时所需的参数、计算方法、剖分格式等不尽相同，数值模拟的过程也不同，但大致相同，本文以TOUGHREACT为例介绍CO2地质储存数值模拟的流程。（一）研究范围的确定一般情况下，独立的天然水文地质系统是计算区最好的选择，它具有自然边界，便于较准确地利用其真实的边界条件，避免人为边界在资料提供上的困难和误差。但是在实际工作中，常常不能完全利用自然边界，这就需要充分利用勘察和长期观测资料等建立人为边界。在确定计算区域时，除了保证范围足够大以外，还应使假定的边界条件尽可能接近真实状态。计算范围的划定应充分考虑研究目的、区域地质构造、储层岩性、储层岩石矿物组成及地下水化学成分等多方面因素。数值模拟时间根据研究目的不同具有不同的时间尺度。就CO2地质储存数值模拟而言，如果不考虑地球化学作用，封存系统在1000年数量级的模拟时间内基本上已达到平衡或稳定。在划定边界时还应考虑CO2在储层中的扩散距离，与研究区地质模型的孔隙度，渗透率等参数关系密切。为了保证所选模型范围边界在模拟期内不影响模拟结果，尽量通过具有相同地质条件的天然CO2气田（藏）进行类比，确定大体的计算范围的边界。如果考虑地球化学反应，由于CO2注入引发的水－岩－气反应对围岩岩性改变较显著，制约着CO2注人的速度和径向运移的距离等。（二）明确研究目的在进行数值模拟以前首先要明确利用数值模拟技术要解决什么样的问题。对于CO2地质储存工程而言，进行数值模拟的目的主要是在CO2地质储存工程实施前，通过数值模拟技术对工程的选址、方案设计进行优化，工程实施期技术指导、运行期监测及后期CO2泄漏的风险评估等进行预测，以指导项目科学、合理地实施，将CO2泄漏风险降至最低。研究目的决定着前期资料的收集类型、地质建模的侧重点、地质模型离散的精密度以及初始、边界条件的处理方式等过程。（三）资料的收集整理1）通过遥感、综合地质调查、物探、钻探和各类样品测试分析等手段获取场地深部地层岩性、地质构造、水文地质、水文地球化学、岩石矿物资料和数据；2）搜集和分析CO2地质储存场地地质岩性、区域构造格架、活动断层与地震活动情况等；3）采用钻井岩心、测井和地震反射方法，调查CO2地质储存场地目标储层和盖层的空间分布形态，埋深、厚度和规模等；4）使用X射线衍射、扫描电镜等方法研究分析封存场地岩石矿物组成、孔隙结构特征及其物理化学性质；5）通过采取浅部、深部含水层水样进行水质全分析，获得储盖层地层水及浅部含水层初始水化学成分。不同的数值模拟软件其数学模型的数值解法不同，空间离散方式也不尽相同，所需的模型参数也有一定的差异，表9-1即为TOUGHREACT数值模拟所需要的主要参数。表9-1 CO2地质储存模拟过程中需要的主要参数（以TOUGHREACT为例）（四）概念模型的建立利用所掌握研究区的地形地貌、气象水文、地层岩性、地质构造、水文地质、水文地球化学等资料，根据研究目的初步建立地质模型雏形。为了解决实际问题，往往要对所建模型及初始条件、边界条件等进行必须要的概化，概化内容包括：1）含水层系统结构的概化：应依据含水层的类型、岩性、厚度、渗透系数等，将内部结构概化为均质、非均质各向同性或者异性的含水层；2）侧向边界和顶，底部边界的概化；根据研究区边界的划定，将侧向边界和顶、底部边界概化为一、二和三类边界条件；3）根据所要解决问题的难易程度，将地质模型概化为一维、二维或者三维模型，并进行合理精密级的网格剖分；4）源汇项的处理：储存库一般位于地表800m以下，上覆有致密的盖层岩石，CO2难以通过越流方式得到补充和排泄，在理想的条件下整个地质储库系统的CO2通过径向对流或抽取井的方式得到补充和排泄。（五）模拟方案的设计针对不同的研究问题，可以根据设计者技术水平、场地经验等设计出不同的方案。利用数值模拟技术对多种方案进行模拟分析，评估方案实施的可行性，并对方案进行优化，最终得出经济、合理的方案。如进行CO2灌注场地的遴选时，要对几个目标场区进行灌注能力，储存潜力、CO2在储层中的扩散速度及距离以及后期CO2泄漏风险等环节进行模拟，然后对模拟输出的数据进行后期处理。通过分析、统计这些数据确定最佳的灌注场地，同时合理划分储盖层的层位、厚度等。对场区灌注能力和储存潜力评估是数值模拟技术所能解决的CO2地质储层领域的又一难题。不同的灌注方式导致CO2进入储层的速度、流量等不同。根据设计者所设计的所有可行的注入方式进行模拟，分析比较不同方案下注入能力和注入量的大小，确定最优的灌注方式。模拟方案的设计要根据研究问题不同视情况而定，而且方案设计的合理与否依赖于设计者自身理论和实际场地经验，不同的研究着对于同一个问题可能设计出不同的方案。我们可以针对这些方案建立不同的模型，通过模拟技术判断其合理性与可行性，最后确定最佳的方案。（六）数值模型及模拟软件的选择数值模拟的关键是地质模型概化、计算精度和计算速度。由于计算的精度取决于离散的程度，而离散的程度又决定了计算的速度，这是一对矛盾，要根据解决问题的需要来选择离散化的程度和计算速度。CO2在储层中的运移、溶解以及与围岩的化学反应形成了一个多相、多组分的反应体系，涉及的主要数学方程有超临界CO2－水的两相流体运动控制方程、溶质运移控制方程和化学反应方程等。在建立数值模型时，通常采用的方法有有限差分法、有限元法和积分有限差法等。由于实际应用时多采用已有的数值模拟软件对CO2地质储存的全过程进行模拟.不涉及软件的开发及程序代码的编写，只需根据研究的需要选择合适的软件进行模拟预测，而软件一旦选定，数学模型和数值模型基本上已经确定。以TOUGHREACT为例，基于上述所建的概念模型建立数学模型。气相和液相统一的偏微分方程为（9-1），咸水的偏微分方程为（9-2），超临界方程为（9-3）。方程中所涉及的相关字符和角标见表9-2。二氧化碳地质储存技术方法概论二氧化碳地质储存技术方法概论表9-2 数学模型中所涉及的符号的意义（七）数值模型的建立1.网格剖分建立了地质模型以后要对研究区离散化，即进行网格的剖分。首先确定离散点，即把所研究的区域按照某种几何形状（如矩形、任意多边形等）分割成网格系统。研究区的边界可以用最接近它的格线近似表示。当网格划分得足够小时，曲折的格线也能够很好地刻画出边界的形状。此过程又称研究区域的离散化（剖分）。在离散化时要遵循以下两条基本原则。1）几何相似。要求物理模拟模型从几何形状方面接近真实被模拟体。2）物理相似。要求离散单元的特性从物理性质方面（含水层结构、水流状态）近似于真实结构在这个区域的物理性质。网格类型大体上可分为规则和不规则两种。规则网格剖分包括矩形、三角形和其他规则图形的剖分（图9-3），而不规则网格则包括不规则多边形等。网格的形状主要根据研究区形状而定。图9-3 网格剖分网格剖分对计算的精度及计算的效率有很重要的影响。精度越高对模拟结果刻画的越精细，但是数据的计算量越大，对计算机的要求也越高。建议在进行地质模型剖分时先采用较粗的网格剖分，如果这种剖分方式下模拟结果合理然后再进行精细化剖分，用于对模拟结果更加详细的刻画。2.参数和初始条件初始条件是指在初始时刻（t=0）时研究区内求解数学模型主要状态变量的初始值。选择的应用软件不同所需的状态变量数量、种类不同。如TOUGHREACT所需的初始主要状态变量包括压力、温度和组分浓度的空间分布。地质参数包括孔隙度、渗透率、密度、压力、温度、毛细压力等参数值。这些数值一部分采用室内实验测得，另一部分采用参考文献的经验值；地层水的化学成分的初始值采用实际地层水的化学分析，主要是8大离子的浓度、盐度和pH 等。如果研究区深部地层中的水样难以获得，如盖层，则采用静态平衡的方法，利用具有与储层相同盐度的咸水与含有原生矿物的地层岩石在原地层环境下进行化学反应，获取平衡状态下的地层水化学成分的初始值；通过岩矿分析、电子扫描、Ⅹ衍射等手段，获得组成CO2地质储层盖层的原生矿物成分体积含量初始值，并根据原生矿物的组成合理判断次生矿物。从原则上讲，初始时刻是可以任意取定的，只要该时刻所需的参数和状态变量值已知即可。因此我们不应该把初始条件理解为研究系统的初始状态。具体如何取，应该视问题的需要、资料来源、计算方便与否等因素而定。3.边界条件边界条件是某一实际问题数学模型具有定解的必要条件之一。地下水流问题和溶质运移问题边界条件的定义不尽相同，但一般概化为以下三种。（1）一类边界条件（Dirichlet条件）解决水流问题时，此类边界条件为在边界上所有点的水头是给定的；对于溶质运移问题，一类边界条件是指研究区边界上的溶质浓度分布已知。解决CO2—水两相流动问题时，此类边界条件为在边界上所有点的压力是给定的。（2）二类边界条件（Neumann条件）当已知某一边界的单位面积流入或流出的流量时，可视作解决流动问题的二类边界；相对溶质运移来讲，此类边界又称给定弥散通量边界，即边界上的弥散通量随时间变化规律已知。（3）三类边界（Cauchy条件）当研究区一部分满足一类Dirichlet条件，而另一部分满足二类Neumann条件时，这类问题称为混合边界问题，称为三类边界。对溶质运移而言，此类边界为边界上溶质通量随时间变化规律已知。在CO2地质储存数值模拟过程中，由于储层地层多在800m以下，地质模型的顶部和底部根据实际需要可以处理为不透边界；为了避免边界对模拟结果的影响，研究区的范围一般比实际CO2所能运移到的范围大得多，因此，在处理四周边界时一般设置为无穷一类边界或不透边界。在确定边界条件时，应根据水文地质条件以及现有的资料来综合考虑。4.源汇项处理在多孔介质中流动和溶质运移的问题中，对流、水动力弥散和溶质源或／和汇，是决定含水层中任一内点上溶质质量时变率的两大因素。源汇项问题在水质与水量计算中以及正确处理对流－弥散方程和渗流基本微分方程中占有重要地位。作为源汇项的方式很多，如越流补给、含水层弹性释放补给以及抽（注）井的补给等。对于深部咸水层CO2地质储存系统而言，系统的顶部一般为具有低渗、低孔的泥岩、页岩等致密性岩层，越流补给较难发生。整个CO2地质储库系统的源汇项主要指对流（如侧向边界）和抽（注）井。（八）模型的校正与验证模型识别是建立地下流体数值模型最重要的环节之一，正确理解和进行拟合对于提高数值模型的仿真性是至关重要的。在有实测结果的情况下如示范工程，可将模拟结果与实测结果进行比较，对相关参数进行适当合理的调整，使模拟结果在给定的误差范围内与实测结果吻合。若误差较大，应该重新检验概念模型的可靠性，甚至重新建立概念模型。在识别校正以后，应采用校正好的模型继续计算，并与未用来识别校正的实际数据比较，验证模型的准确性和可靠性。若存在较大误差，需重复前面的过程。在没有实测结果的情况下，数值模型的可靠性可通过类比相关资料或根据个人经验和理论判断。（九）模拟预测模型预测是实施数值模拟技术的主要目的。对于CO2地质储存工程而言，由于CO2地质储存技术的提出为时尚短，针对CO2在深部咸水层中的运移、扩散、与地层水和围岩产生的化学反应，以及由于CO2灌注引起的储盖层物理、化学性质变化研究均处于研究和发展阶段。因此，在工程实施过程中急需具有技术指导性的工具产生，避免造成投资的浪费和CO2泄漏等风险的出现。利用经过识别校正与验证过的数值模型对CO2地质储存过程进行模拟预测，有针对性地对模拟数据进行后期处理，如统计分析、比较等手段对结果进行解译，以此达到场地的优选，目标储层灌注能力、储存潜力的评估，CO：扩散运移途径和速度、不同捕集方式封存量及它们之间的时空转化等过程的详细刻画与模拟仿真等目的。同时可以预测CO2在已有、重新激活或新生成的裂隙中逃逸的可能性及时间、CO2泄漏风险评估以及评价CO2泄漏对浅层地下水的水质、水量及对地表环境的影响等。上述结果的分析只是数值模拟技术所能解决问题的冰山一角。对于数值模拟结果的处理要根据所研究的目的进行有针对性的提取和解译。通过对处理后的数据进行总结分析，发现问题从而解决问题，并掌握内在规律，为CO2地质储存工程的前期设计、工程实施、中期监测管理提供理论支持和科学的技术指导，并可以提前开展风险预测，尽早制定预案防范CO2地质储存工程实施及运行过程中可能出现的隐患。

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

数值模拟过程中,面临的最大的问题可能是收敛性问题。模型不收敛的原因很多，网格参数，属性参数，流体PVT参数，岩石相渗曲线，毛管压力曲线，相渗曲线端点标定，初始化，井轨迹，垂直管流表都会造成模型不收敛。模拟计算的时间取决于时间步的大小，如果模型没有发生时间步的截断而且能保持长的时间步，那表明该模型没有收敛性问题，反之如果经常发生时间步截断，那模型计算将很慢，收敛性差。时间步的大小主要取决于非线形迭代次数。如果模型只用一次非线形迭代计算就可以收敛，那表明模型很容易收敛，如果需要2到3次，模型较易收敛，如果需要4到9次，那模型不易收敛，大于10次的化模型可能有问题，如果大于12次，时间步将截断。

一朝去京国，十载客梁园。猛犬吠九关，杀人愤精魂。

（不是专业人士,说说想法……） 这种问题提出,大约是在工程中遇到了求解数学问题的情况（偏微,常微等） 解析分析——用数学分析的方法,比如微分、积分、特殊方程等,对实际情况进行模拟,列出方程,用解析的方式,求出“比较正规”的函数解（这样可以解决的问题是很有限的,因为在实际问题中,有很多是不能用函数简单模拟的；而且这样的方式,解决方法灵活多变,不利于机器模拟——虽然有软件可以做到） 数值模拟——以不用具体的函数表达式,而是用多个点的数值表示函数的方法,来解实际问题的解法,恩,一般来说,差分法是最常用的（这样做的在于应用范围很广,但是计算误差必须估计） 盼有所帮助. 补充一点：另外一种可能,解析分析就是用解析的方法求解,数值模拟用于检验.

如前所述，基于镜质组化学反应动力学的煤热演化史数值模拟方法经历了简单函数关系模式、受热时间－经验法模式、反应活化能－温度函数模式和平行反应化学动力学模式四个发展阶段（秦勇等，1995）。其中，目前在化石能源地质界最为常用的是TTI法、LOM法和Tissot法，近年建立起来的EASY%Ro方法也已成为数值模拟的重要发展方向。（一）时间－温度指数（TTI）法这种方法的基本思想由Lopatin（1971）首先提出，后经Waples（1980）根据31个盆地402个样品的实测资料加以完善，建立起了其数值模拟的方法体系，故又称为Lopatin-Waples法。TTI法遵循化学反应动力学的基本法则来衡量温度和时间两个因素对煤化作用的贡献，即煤化作用速度随受热温度的增高而增大，温度每增加10℃煤化反应速度增加一倍。由此，可将TTI数学模式用连续函数加以表示：山西南部煤化作用及其古地热系统：兼论煤化作用的控气地质机理T=Ts＋G·D式中：T——煤层受热的古地温温度，℃；Ts—－古地表温度，℃；G——古地温梯度，℃/m;D—煤层埋藏深度，m。上述连续函数模式不便于常规运算，故在实际工作中通常采用分段积分方式，根据如下方程进行计算：TTI＝∑γn·△t式中：γ——温度效应因子，考虑到温度每增高10℃煤级量值增加一倍的关系，取r=2;n——温度指数，其大小取决于温度（表5-1）；△t——煤层在某一温度段的受热时间，Ma。表5-1　TTI数学模式中温度指数（n）与煤层受热温度段的对应关系注：温度每增（减）10℃，温度指数增（减）1。TTI值本身可以作为煤级或有机质成熟度指标，Waples（1980）也提供了TTI值与镜质组最大油浸反射率的关系（表5-2）。此外，Wood（1988）通过回归拟合，将二者之间关系采用下式表示：lg%Ro＝—0.006（lgTTI）3+0.042（lgTTI）2+0.162（lgTTI）—0.397表5-2　时间温度指数（TTI）与镜质组反射率（Ro）的相关关系应予指出，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）Ⅰ—褐煤；Ⅱ1—亚烟煤；Ⅱ2—高挥发性烟煤；Ⅱ3—中挥发性烟煤；Ⅱ4—低挥发性烟煤；Ⅲ1—半无烟煤；Ⅲ2—无烟煤1Btu=1055.06J在实际工作中，LOM法采用图解方法求取有机质成熟度或镜质组反射率（图5-2）。国内外应用效果表明，利用这种方法推算出的盆地古地温温度是比较可靠的，与由综合研究所得出的结果以及盆地实际情况较为吻合（Vote,1981；周中毅等，1983,1984,1985）。

经济学的数值模拟一般分为计量和非计量模拟。计量模拟和一般的统计学模拟非常类似。一般都要假设一个data generating process (DGP)，也就是数据生成过程。一般常见的DGP都是用随机生成的正太分布随机变量来模拟。然后根据希望检验的计量模型，或者计量理论，用随机生成的变量进行计量分析。非计量的数值模拟种类很多。一般都和经济学理论相类似。其具体过程就是，首先要制定一个量化经济学理论的规则。比如，检验人的风险偏好对投资结果的影响。假设人的偏好为0到1.0为风险完全规避，1为风险喜好，0.5为风险中性。之后可以从uniform distribution中，随机生成大量变量，并以此变量为根据设定人的风险偏好，通过所要检验的公式，求出最终的投资结果。并通过对投资结果和相对应风险偏好的分析，得出结论。要知道，所有数值模拟都是主观设定的。所以，学术上，最多只能对已经被从数学理论上证明了的经济学原理，进行再次验证。数值模拟，是不能作为理论论证的依据的。最多，对理论论证进行更加形象的表述。

数值模拟也叫计算机模拟。依靠电子计算机，结合有限元或有限容积的概念，通过数值计算和图像显示的方法，达到对工程问题和物理问题乃至自然界各类问题研究的目的。在工程应用中，很多情况下无法进行实验，如采矿问题等，数值模拟内部程序有相应的计算方法，能模拟较复杂过程。数值模拟成本低廉，只需在计算机上进行模拟和数据处理。数值模拟可以任意施加各种方向的载荷，可以施加实验方法达不到的条件。数值模拟促进了试验的发展，对试验方案的科学制定、试验过程中测点的最佳位置、仪表量程等的确定提供更可靠的理论指导。因此数值模拟方法在监测、设备开发、优化、效果预测方面体现了重要价值，数值模拟计算已经与理论分析、试验研究成为科学技术探索研究的三个相互依存、不可缺少的手段。

你好，很高兴回答你的问题，三维数值模拟软件有哪些？1、MODO12集3D绘画，动画与渲染于一体。可以感受高级多边形细分曲面、建模、雕刻、3D绘画、动画与渲染于一体的便利性。2、Blender是开源的多平台全能三维动画软件。喜欢3D绘图的玩家们可以使用Blender制作出自己喜爱的3D模型。3、Maya功能完善，渲染真实感极强。一般影视动画是用maya来做的，制作效率高。4、RhinoRhinocero，简称Rhino，又叫犀牛，是一款三维建模工具，拥有AutoCAD基础的初学者更易于掌握犀牛。5、Silo：一款专注且闪电般快速的独立3d建模器，着重于3D设计、动画、录影游戏制作和传达想法等领域的模型建造和塑形。

数值模拟也叫计算机模拟。依靠电子计算机，结合有限元或有限容积的概念，通过数值计算和图像显示的方法，达到对工程问题和物理问题乃至自然界各类问题研究的目的。

6.2.2.1 地下水系统的模拟模型潜水含水层的模拟模型：含有协变量的地下水动态规划管理模型研究式中：D——地下水系统的模拟渗流区域；h——地下水位（m）；M——含水层厚度（m）；K——渗透系数（m/d）；T——导水系数（m2/d）；t——时间（d）；μ——给水度；μ*——承压含水层贮水系数；Γ1，Γ2——一类及二类边界；n→——边界的外法线方向；h1（x，y，t）——第一类边界的地下水位（m）；h0（x，y）——初始地下水位（m）；q（x，y，t）——第二类边界的单宽流量（m2/d）；P——地下水开采量（m3/d）；ε——天然条件的源汇项（m3/d），如降水入渗补给、蒸发排泄等；Q——与地下水位有关的源汇项流量（m3/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.1653×108m3 ，2008年地下水开采量为5.2201×108m3 ，平均年增长率定为3.83%。按该增长率对2012年开采量进行预测，即6.067×108m3。6.2.2.5.3 预报结果及分析根据地下水开采量平均年增长率3.83%对开采量进行预测，即计算模拟区2012年地下水开采量为6.067×108m3。从预报结果的水位等值线图（图6.14和图6.15）可以看出，潜水水位和承压水位都有一定程度的下降，部分地区水位下降明显。因此，对水位持续下降地段，应采取科学管理和一定的控制措施。此外，由于地下水开采量的增加以及当地混合开采地下水，致使潜水层和承压水层直接产生了水力联系，承压含水层的承压性质正在减弱，向潜水转化。图6.14 2012年9月1日潜水等水位线图针对上述情况，在地下水开发利用中应采取一定的对策：①要科学地布井，有计划地进行浅、中、深层分层开采地下水；②开展地下水长期监测，以便掌握区域水情的变化趋势，发现问题及时解决；③加强水资源管理工作，合理开发利用地下水资源，确定开采量的大小和开采的时段长度，控制超量开采局面。图6.15 2012年9月1日承压水等水位线图

在工程应用上,通常数值模拟结果和实验（试验）结果之间的误差在±20%以内认为是可以接受的.误差分析的方法很多,曾见过一篇采用函数分析法来分析误差的方法,好像可以分析试验曲线与数值模拟曲线之间线型的差别.

常用的数值计算方法有有限差分法和有限单元法。由于有限单元法中的集中储量有限元方法较通常的有限元法具有更多的优点，而且在边界条件的处理上，集中储量有限元法比有限差分法更符合实际，它考虑了边界节点的均衡单元的储水量变化（吴金全，1989）。图1.4.3 均衡区域示意图（一）集中贮量有限元公式推导取单位水平面积、高度为计算层厚度的土柱进行研究（图1.4.3），将土柱（计算区域）垂直向上剖分为n个单元，空间步长为Δz，节点编号为0，1，2，…，n-1，n，在Δt时段内（Δt=tj+1>-tj），对任一内节点i所代表的均衡区zi-1/2到zi+1/2（图1.4.3）之间的土体列水量均衡方程（暂先不考虑极系吸水项）。1.内节点（i=1，2，…，n-1）由达西定律 （h 方程）得：（1）通过zi-1/2断面的水流通量（流入量）为：土壤水盐运移数值模拟（2）通过zi+1/2断面的水流通量（流出量）为：土壤水盐运移数值模拟（3）均衡区域Δz内储水量的变化量（增量）为：土壤水盐运移数值模拟根据质量守恒原理（流入量-流出量=储存量的变化量）得：Δqi=qi-1/2-qi+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）变为：土壤水盐运移数值模拟令： ， ，r3=Δt得：土壤水盐运移数值模拟式中：i=1，2，…，n-1。令土壤水盐运移数值模拟将式（1.4.24）代入式（1.4.23）得：土壤水盐运移数值模拟2.边界节点的处理（1）上边界节点i=0处的方程为：土壤水盐运移数值模拟式中： 为 ，Rj+1/2，为时段平均入渗强度，Ej+1/2为时段平均表土蒸发强度。令（1.4.26）式中：土壤水盐运移数值模拟则（1.4.27）式变为：土壤水盐运移数值模拟（2）下边界节点i=n为第一类边界节点，hn已知，故不需列方程计算，这样第n-1个方程可简化为：土壤水盐运移数值模拟式中：土壤水盐运移数值模拟3.方程组综合内节点和边界节点方程，从而得如下代数方程组：土壤水盐运移数值模拟方程组式（1.4.31）中：b0、c0、f0按式（1.4.27）式计算，fn-1按式（1.4.30）式计算，其余αi、bi、ci按式（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.模型验证进行数值模拟，首先进行模型验证。模型验证时，上边界条件表达式中的θ10由实际观测资料给出。根据有作物生长条件下土壤水分运动的基本方程和差分方程，在已知初始条件和边界条件时，模型验证可以通过以下步骤进行：①根据实测初始负压剖面的分布，用三次样条插值给出各节点上的初始值；②计算蒸发量E；③计算根系吸水层厚度Lr及吸水率S；④根据差分方程计算时段末负压值h。模型验证时，以计算起始时刻的实测负压剖面（或含水率剖面）作为初始剖面，空间步长选用1cm，根据最底部负压计测点和中子仪测点，剖面深度为120cm（含水率剖面为130cm），时间步长选用1h，迭代相对误差e≤0.01。计算中输入的大量信息，如各节点的初始负压值、降雨量、降雨日期、水面蒸发强度、根层土壤含水率等均以数据文件的形式提供。由于大田盖200kg/亩和盖600kg/亩只进行了含水率观测，计算时先将含水率θ转化为负压h，计算结束后再将负压h转化为θ。大田盖400kg/亩有负压资料，可直接用负压h计算。上边界条件由E/E0-θ关系给出。数值模拟按覆盖量（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/E0-θ经验公式来反映覆盖对水分运移的影响。因此，本文提出的模型和数值方法可以用来模拟秸秆覆盖条件下田间土壤水分的运动，可对田间土壤水分动态作中短期预报。图1.4.4 实测值与模拟值对比图3.模型的应用——预报数值模拟的目的之一就是进行预报，根据气象部门提供的降雨量及水面蒸发强度等气象资料，使用验证过的模型进行田间土壤水分动态预报。本文使用实际发生过的降雨量及水面蒸发强度系列资料进行预报，用实测负压资料检验预报结果。程序的运行见图1.4.5。检验的实测资料选用大田覆盖400kg/亩的资料，分生育阶段（苗期、拔节、灌浆-成熟）进行。从图1.4.6可以看出，预报值和实测值吻合较好。图1.4.5 双层结构有根系吸水项垂向一维土壤水数值模拟框图图1.4.6 预报值与实测值对比图（大田盖400kg/亩）

Plaxis、ABAQUS软件。PLAXIS2D/3D程序是由荷兰PLAXISB.V.公司推出的一系列功能强大的通用岩土有限元计算软件，已广泛应用于各种复杂岩土工程项目的有限元分析中，如:大型基坑与周边环境相互影响、盾构隧道施工与周边既有建筑物相互作用、大型桩筏基础(桥桩基础)与邻近基坑的相互影响、板桩码头应力变形分析、库水位骤升骤降对坝体稳定性的影响、软土地基固结排水分析、基坑降水渗流分析及完全流固耦合分析、建筑物自由振动及地震荷载作用下的动力分析、边坡开挖及加固后稳定性分析等等。PLAXIS系列程序以其专业、高效、强大、稳定等特点得到世界各地岩土工程专业人员的广泛认可，日渐成为其日常工作中不可或缺的数值分析工具。尤其在欧洲、新加坡、马来西亚、香港等地应用广泛，PALXIS2D甚至用于常规的二维设计计算中。截至2012年初，世界范围内PLAXIS用户多达16000多家;其中中国用户已有百余家，涵盖了铁路、电力、石化、建筑、航务、冶金等行业设计院、高校、科研院所及少量施工单位。ABAQUS是一套功能强大的工程模拟的有限元软件，其解决问题的范围从相对简单的线性分析到许多复杂的非线性问题。ABAQUS包括一个丰富的、可模拟任意几何形状的单元库。并拥有各种类型的材料模型库，可以模拟典型工程材料的性能，其中包括金属、橡胶、高分子材料、复合材料、钢筋混凝土、可压缩超弹性泡沫材料以及土壤和岩石等地质材料，作为通用的模拟工具，ABAQUS除了能解决大量结构(应力/位移)问题，还可以模拟其他工程领域的许多问题，例如热传导、质量扩散、热电耦合分析、声学分析、岩土力学分析(流体渗透/应力耦合分析)及压电介质分析。

历史拟合是一个相当复杂的工作,历史拟合的过程需要消耗很多人力和机时,它常常占用了模拟研究的大部分费用。油藏模型被建立起来以后，它是否完全反映油气藏实际，并未经过检验。只有利用将生产和注入的历史数据输入模型并运行模拟器，再将计算的结果与油气藏的实际动态相比，才能确定模型中采用的油气藏描述是否是有效的。合肥辰工就是做石油工程软件这一块的，跟中科大渗流实验室做合作的。他们好像也做油藏数值模拟的历史拟合这方面的，你可以咨询他们。

肯定是试验好发论文啊数值模拟不真实，而且有许多的假设，不如试验客观和现实。尤其是混凝土结构，混凝土是各项异性材料，在模拟中要等价于各项同性材料。而且假如试验发论文的时候，你介绍一下试验原理和试验准备工作和试验步骤就能用去大部分的篇幅。模拟也是站在试验的基础上才能发论文的，没有试验单靠模拟站不住脚。

没有准确的范围，有的相差很大，关键看假设的条件是否合适。

数值模拟方法有很多种，最常用的是有限差分法和有限元方法。

这个问题。。。这么说吧，数值这个词，是数值计算方法的简称。那么，仿真模拟就是采用仿真的手段模拟真实场景。比如虚拟练车？数值模拟，指用数值计算方法来首先对物理场景的模拟，以研究一些实验仪器不易实现的物理现象。数值仿真，指用数值计算方法实现对物理场景的仿真。比如，采用数值模拟的方法模拟核电站调试，或者说演练吧。其实后面这两个词差不多，一个是模拟，即场景再现，一个是仿真，即模仿真实，一个从过程角度讲，一个是从意义角度讲。胡乱打字，楼下高手补充。

　　工程数学·计算机

(1)数学模型根据上述非均质、水平方向各向同性、垂直方向存在变异的空间三维结构、非稳定热流体系统,可用如下微分方程的定解问题来描述:图5-16 Ng组三维实体模型建模演示图沉积盆地型地热田勘查开发与利用式中:Ω为渗流区域;h为热流体的水位标高(m);Kx,Ky,Kz分别为热流体储积层中x,y,z方向的渗透系数(m/d);Kn为边界面法向方向热流体储积层渗透系数(m/d);Kz→为边界内侧热流体储积层垂向渗透系数(m/d);S为储水率(1/m);σ为通用水头边界的阻力系数,σ=L/K1,L为模型边界到通用水头边界的水平距离(m),K1为模型边界到通用水头边界之间热流体储积层平均渗透系数(m/d);ε为源汇项(1/d);h0为初始水位分布(m),h0=h0(x,y,z);hb为通用水头边界上的水位分布(m); 为边界面的法线方向; 为热储层顶、底边界面的法线方向;Γ1为通用水头边界;Γ2为渗流区域的侧向隔水边界;Γ3为渗流区域顶、底边界;(2)数值模型绝大部分数学模型是无法用解析法求解的。数值化就是将数学模型转化为可解的数值模型。常用数值化有限单元法和有限差分法。本次工作根据热储层地热地质概念模型特征、实际资料丰富程度以及精度要求,应用有限差分法对馆陶组热储进行评价。模型软件选用模拟地下水流动和污染物运移等特性的计算机程序——Modflow (The modular finite-difference groundwater flow model)。该系统是三维地下水流和溶质运移模拟评价软件。(3)计算区单元剖分在进行热流体动态模拟前,首先必须对模拟区进行单元剖分,剖分时除了遵循一般的剖分原则外,还应充分考虑如下实际情况:①计算区的边界、岩性分区边界、行政分区边界等;②观测孔尽量放在剖分单元的结点上;③集中开采区域,由于水力坡度较大,流场变化趋势较大,所以剖分时适当加密。因此在平面上采用不等距矩形剖分法,在开采井集中区域适当加密剖分间距。沿边界将外部和缺失区域定义为非活动单元格。剖分为60行、55列,其中有1580个非活动单元。最大剖分单元约1.8×1.8km2正方形单元,最小剖分单元为0.28×0.38km2长方形单元。而在垂向上热储层按NgⅠ、中间弱透水层NgⅡ和NgⅢ进行剖分。则Ng有(60×55-1580)×3三层一共是5160个剖分单元,具体情况见图5-17和图5-18。(4)边界条件Ng热储层可以概化出3种类型的边界条件:即①自然流量边界;②零流量隔水边界;③沿沧东断裂形成的内部补给边界(图5-19)。现对不同类型边界条件的处理方法分述如下:A.自然流量水头边界自然流量边界是利用MODFLOW软件中的变水头边界程序包(GHB)中的通用水头边界来处理。使用这个边界条件的目的是避免不必要的向外扩展模型区域,以满足要素影响模型中的水头。因此,变水头边界条件通常被指定在模型域的外面。从外部源进入或者流出计算单元的热流体量与该计算单元水头和外部源水头之差成正比。在通用水头边界单元上,从单元外部的热流体进入或流出计算单元的热流体量与该计算单元水头h和通用水头边界热流体水头hb之差成正比。据此,可确立计算单元热流体量与计算单元水头间的线性关系,即:图5-17 Ng平面单元剖分图图5-18 Ng3D单元剖分图图5-19 Ng计算边界条件图沉积盆地型地热田勘查开发与利用式中:Qb为从通用水头边界进入计算单元的流量(m3/d);Cb为通用水头边界与计算单元间的水力传导系数(m2/d);(L×W)为边界与网格单元的截面积(m2);K为从模型网格中分离源/汇的含水层物质的平均渗透系数(m/d);D为边界源/汇到模型网格间的距离(m);hb为通用水头边界上的热流体水头(m);h为计算单元热流体水头(m)。当Cb很大时,GHB子程序包的作用与定水头边界相同。因此通用水头边界上的水位和水力传导系数Cb决定了边界的流入流出量Qb的大小。B.零流量隔水边界除了新近系馆陶组西侧边界和部分缺失区域定为零流量隔水边界,各热储层顶底板也定为隔水边界。毫无疑问,对于隔水边界,其计算单元与外部域的流入流出量Qf值为零。C.沿沧东断裂形成的内部补给边界由沧东断裂特性,沿沧东断裂形成了连通上下各热储层地热流体的一个通道。滨海塘沽区、大港区对Ng热流体的过量开采可导致该层热流体压力水头的过量下降,从而激发其下部各热储层对该层的补给,因此,沿沧东断裂一带应考虑其他含水层的垂向补给。本次工作通过数值模型软件 MODFLOW提供的河流子程序包来模拟沧东断裂对计算区热流体的补给量。其计算公式如下:沉积盆地型地热田勘查开发与利用式中:Qfau为从断裂进入计算单元的流量(m3/d),进入取正值,流出取负值;Cfau为断裂带外部源与断裂带计算单元间的水力传导系数(m2/d);L、W为通过计算单元网格的断层长度和宽度(m);K为断裂破碎带间计算单元与下层热储层间平均渗透系数(m/d);M为断裂破碎带间计算单元与下层热储层间厚度(m);hfau为断裂破碎带间计算单元下层热储层热流体水头(m);h为断裂破碎带间计算单元水头(m)。沧东断裂带下层基岩热流体储积层水头一般都高于Ng热流体水头,断裂带的渗流补给量一般是正值。即下层基岩热流体通过断裂带顶托补给上层Ng热储层。但由于受原始资料限制,所求的流入流出量Qfau精度较粗,不能准确给出。只能在流场模拟中,根据流场模拟结果来进行校准。(5)模拟期的选择和时间离散数值模拟自2002年1月始,至2006年10月结束,共5年10个月时间作为该层模型识别整个模拟时间长度。根据地热流体的实际开采利用特点,将一个开采周期年分为两个应力期(时段):即取暖期(11月10日至次年3月25日,共135天)和非取暖期(3月26日至11月9日,共230天),每个应力期又分成基本以15天为单位的时间步长。(6)补给、径流、排泄条件对于地热流体,由于埋藏较深,相应补排项比较单一:补给来源主要包括侧向径流补给、随着热储层热流体压力变化引起的热储层岩石骨架释水、部分同组层间越流补给量以及人工回灌量。新近系馆陶组过量开采可能通过沧东断裂激发部分基岩补给量;排泄方式主要为下游单元的侧向径流排泄、人工开采以及越流补给同组其他超采层。对于侧向径流以及内部补给边界的处理前文已有介绍,在此不再赘述。以下主要分析其他补排项的处理。A.弹性释水引起岩石骨架释水的主要原因是热储层热流体压力变化,而岩石骨架释水量的大小,不但受压力变化的影响,热储层弹性释水系数的大小也是一个至关重要的参数。其值由以下公式计算而得Qi=L×W×S*×Δh 5-25式中:Qi为单元弹性释水量(m3);L为单元长度(m);W为单元宽度(m);S*为单元所在地热流体富集段弹性释水系数(无量纲);Δh为单元热流体压力水头变化量(m)B.同组层间越流补给、排泄量由于采用三维水流模型,层间越流量由MODFLOW模型软件根据所赋相关参数自动进行处理。C.泥岩释水泥岩释水是指土体中无所不在的结合水和毛细水由于土体压缩而释放的水量,泥岩释水可引起地面沉降。很明显,由于地热开采引起的地面沉降有限,泥岩释水量也不会太大,但也不应对其忽视。据1998年12月《天津市塘沽区地热流体回灌研究》(李明朗,1998)指出,热流体开采量中黏性土释水量约占总释水量(弹性释水量+粘土释水)的1/5;而总释水量高时约占到热流体开采量40%左右。但随着资源开采的延续,黏性土层释水后压密度增加,其释水量应呈逐步衰减之势。而随着热储层压力下降,总释水量也会逐步减少。D.人工开采、回灌量的处理对于人工开采、回灌量,在模型中都是以开采井的形式给入,开采量为正值而回灌量为负值。(7)水文地质参数水文地质参数主要是指渗透系数K和弹性释水系数S*,前者K体现了热储层流体的渗透能力;而后者S*则体现了热储层储水(释水)能力。水文地质参数选用是依据三维模型的计算需要,主要采取分区赋值的方法。根据Ng热储层的埋藏分布和控热构造分布特征,利用所收集的降压试验资料进行综合统计分析,按所求水文地质参数进行初始分区(见图5-20)。有条件的区域项目前期补充地热井降压试验工作,空白区域根据邻区资料结合地层构造变化特征给出估计值。

直接数值模拟不需要对湍流建立模型，对于流动的控制方程直接采用数值计算求解。由于湍流是多尺度的不规则流动，要获得所有尺度的流动信息，对于空间和时间分辨率需求很高，因而计算量大、耗时多、对于计算机内存依赖性强。直接数值模拟只能计算雷诺数较低的简单湍流运动，例如槽道或圆管湍流，现如今它还难以预测复杂湍流运动。一个直接数值模拟（DNS）是一个模拟的计算流体力学中的Navier - Stokes方程进行数值求解无任何湍流模型。 这意味着整个范围的空间和时间尺度湍流时，必须同时解决。 所有的空间尺度湍流必须在计算网格中解决，从最小的耗散尺度（ Kolmogorov microscales），到最大的积分尺度L，L与大部分包含动能的运动有关。 Kolmogorov尺度，η定义为由其中ν为运动粘度和ε是速度动能耗散。 另一方面，整体规模取决于通常在空间尺度上的边界条件。为了满足这些解的要求，节点数N和网格步长h必须满足使整体规模是包含在计算域内，并且这样可以得到Kolmogorov尺度。由于其中u"是均方根（RMS）的的速度，以前的关系，意味着一个立体的DNS需要满足数量网点北区

1.数学模型描述含水系统地下水渗流的数学模型为：华北煤田排水供水环保结合优化管理式中：H为含水层水位（m）；H0为含水层初始水位（m）；T为含水层的导水系数（m2/d）；μ为潜水含水层给水度；承压含水层为贮水系数；q为二类边界的单宽流量（m3/d·m）；ω为降雨入渗补给强度（m/d）；Q为水井开采量、突水点突水量（m3/d）；n为边界外法线；G为计算区域；Γ2为第二类边界。2.计算区域剖分与时间步长选用三角形网格有限差分法求解数学模型。将模拟计算区域剖分成358个单元，共210个节点，其中内节点150个，边界节点60个。剖分时使各种参数分区界线及地下水动态观测孔和开采井全部落在节点上；对枣沟和观孟前水源地及六斜坡突水点等开采较大的地区加密网格剖分，见图6-8。考虑到大多数水文要素数据均按月统计，所以选取的计算时间单位为月，将所有的源、汇项及边界的数据均逐月整理成月平均值。3.数值模拟模型含水层第i节点均衡域的渗流差分方程为：华北煤田排水供水环保结合优化管理式中：fβ为属于第i均衡域的第β三角形单元的小四边形面积；Δβ是以i节点为公共顶点的第β三角形单元的面积；华北煤田排水供水环保结合优化管理图6-8 研究区网格剖分图以下为几何量：华北煤田排水供水环保结合优化管理式中：下标i、j、m为三角形单元以逆时针排序的三个节点编号，k为时阶。按同样形式写出剖分网格上所有节点均衡域的差分方程，则构成一个庞大的代数方程组。这就是计算区含水层渗流系统的数值模拟模型。该数值模型的矩阵形式为：华北煤田排水供水环保结合优化管理式中：［A］为导水矩阵；［D］为贮水矩阵；Ht为已知水头向量；Ht+1为待求水头向量；F为已知常数向量。该模型系统采用迭代法求解［150～151］。其计算框图见图6-9。图6-9 水资源模拟系统结构图

1、 复杂系统的设计，例如飞机和物流系统。2、 电力系统仿真。3、 工业过程的设计，例如化学加工厂。4、 过程工程仿真工具。5、 用于设计机器人和机器人控制算法的机器人模拟器。6、 模拟汽车碰撞以测试新车型的安全机制。人们对计算机仿真的可靠性和信任取决于仿真模型的有效性，因此验证和确认对于计算机仿真的发展至关重要。计算机模拟的另一个重要方面是结果的可重复性，这意味着模拟模型不应为每次执行提供不同的答案。尽管这看起来很明显，但这是随机模拟中的一个特别注意点，其中随机数实际上应该是半随机数。汽车制造商利用计算机模拟来测试新设计的安全特性。通过在物理模拟环境中构建汽车副本，他们可以节省数十万美元，否则这些成本将用于构建和测试独特的原型。工程师可以一次完成几毫秒的模拟，以确定施加在原型每个部分上的确切应力。计算机图形可用于显示计算机模拟的结果。动画可用于实时体验模拟，例如在训练模拟中。在某些情况下，动画也可能在比实时更快甚至比实时模式更慢的情况下有用。例如，在模拟人类撤离建筑物时，比实时动画更快的动画可以用于可视化队列的建立。此外，模拟结果通常使用各种科学可视化方式聚合成静态图像。可视化以前，计算机模拟的输出数据有时会以表格或矩阵的形式呈现，显示数据如何受到模拟参数的大量变化的影响。矩阵格式的使用与数学模型中矩阵概念的传统使用有关。然而，心理学家和其他人指出，人类可以通过查看图表，甚至是由数据生成的运动图像或电影来快速感知趋势，如计算机生成图像(CGI) 动画所显示的那样。尽管观察者不一定能读出数字或引用数学公式，但通过观察移动的天气图，他们可能能够比扫描雨云表更快地预测事件（并“看到雨正朝着他们的方向前进”）坐标。这种超越数字和公式世界的强烈图形显示有时也会导致输出缺少坐标网格或省略时间戳，就好像偏离数字数据显示太远了。今天，天气预报模型倾向于平衡移动雨/雪云的视图与使用数字坐标和事件数字时间戳的地图。同样，CAT 扫描的CGI 计算机模拟可以模拟肿瘤在长时间的医疗过程中如何缩小或变化，随着肿瘤的变化，将时间的流逝呈现为可见人头的旋转视图。

在三层模型计算中，共划分节点1659 个，划分单元2316 个。图7-5a、b 为三层模型结构和约束状态图。其中选取模型下底部为南，顶部为正北方向，图中的实体用深浅不同颜色来区分不同的材料性质。底部的箭头方向为施加的压应力方向。图7-5 第一个计算模型结构计算结果表明，岩石底部的塑性层向北挤压变形将对上层物质产生明显影响，从图7-6a、b模型网格位移结果可看出，模型的三层材料均发生了变形。模型在 x 方向（EW向）的位移表现出上下均匀分布的特点（图7-6c）。由于模型在 EW两侧未加约束，因此位移向两侧移动，在数值上关于 y 轴对称，但方向向反。y 方向的位移在受力底部为最大，向上向北逐渐递减（图7-6d）。图7-6 第一个模型位移计算结果在应力分布特征上，模型在空间各方向上表现出不同的性质，模型 x 方向（EW 向）应力整体近于均匀分布（图7-7a、b），应力最大部位在受力的南部边界附近，在其上部为自由边界，出现 x 方向极小值。在 y 方向的应力分布与此特征不同，其应力分布在受力的底部塑性层中由南向北逐渐渐少，反映出塑性物质应力分布的特点，上层应力在最南部自由边界处最小为零，向北逐渐加大，反映出了塑性层向上传递应力的特征（图7-7c）。z方向应力由于模型底部得到约束，为近匀分布特征，只在加力边界和约束点处出现局部扰动（图7-7d）。这一结果说明，如果大陆地壳挤压力源来自其底部的塑性层界面摩擦拖拽，那么其上部弹性层的应力分布将受塑性层位应力分布特点的制约，应力分布具有塑性物质特征，从而在变形的总体形态上产生类似塑性物质受挤压的特点。但其本身变形仍遵循弹性定律。图7-7 第一个模型应力分布计算结果模型的应变计算结果与应力分布类似，这里只给出 x 方向和y 方向的应变结果（图7-8a、b、c），从 y 方向应变结果可以明显看出塑性底层应变往北渐渐减弱的过程，与此相应，在表层的弹性层也随之出现应变呈 EW带状分布的特点。若此时表层应变达到其极限强度，其发育位置应在一条横向分布的带中，其破裂无论如何不会出现锐角指向南北的情况。从模型整体位移矢量分布特征上看（图7-8d），模型呈现整体向北运动的趋势，向北位移减小并有向两侧运动的特征。本区的一个显著特点是在表层SN向断裂极其发育，呈现出密集的SN向断裂束，而物探资料反映出本区深部有呈EW向带状展布的特点，并有EW向隐伏断层存在。随着深度增加，这一特征反映越显著，并同青藏高原南缘喜马拉雅山脉东南缘龙门山脉底部EW向构造带相连，反映出其形成机制与青藏高原的隆升有关。具有明显的区域活动特征，而不是一种局部表现。本区的这一构造特征用上述模型显然无法作出完整的解释。说明本区除了有规模巨大的底部塑性挤压驱动外还有自己独特的结构特点，地震测深和电法资料反映在本区中、下地壳之间有低速高导体发育，发育深度在26～30 km，其下部对应天然地震多发带（图4-3）。图7-8 第一个模型应变分布及位移矢量计算结果地壳中的软弱夹层、高压流体对水平拆离滑脱面及逆冲推覆构造产生的作用已经成为共识，同时对由构造挤压导致软弱岩层出现高压流体的研究正逐渐得到重视。构造强烈发育并伴随有低速高导体的广泛发育是本区构造活动的一个显著特点。因此第二个模型就是在前述模型的基础上考虑高导体的存在而模拟了各圈层相互作用的情况。这里将物探资料得到的低速高导体作软弱带处理，在两弹性层之间加入一软弱层，其物性见表7-2。在保持其他条件不变的情况下讨论各圈层的作用关系，模型结构和约束条件见图7-9，其中约束条件与第一个模型相同。图7-9 第二个计算模型结构计算结果反映出软弱层的存在对底部塑性层应力和位移的向上层传递有很大的屏蔽作用。模型在 x、y 方向上的位移分布与模型1大致相同（图7-10a、b），但在 y 方向上的位移分布向上减弱分带特征明显，经过软弱夹层出现明显位移间断。出现与地壳表层薄皮构造运动状态相反的底部层位位移大于上部层位的滑脱现象。从模型 x 的方向（EW方向）应力分布图中（图 7-10c）可以看到，应变的总体形态与模型1相似，为近似均匀的应变分布，但在软弱夹层处出现明显应力变化，在软弱带上应力为一低值区。在 y 方向应力分布图中（图7-10d）可以看出，经过软弱层的上下层位应变分布明显出现不连续，上部层位应变明显减弱，在软弱带中可见低的应变带分布。在前面两个模型计算结果的基础上，考虑到本区地壳表层受印度板块的强烈挤压作用，并考虑到本区主要的断裂控制作用，在建立第三个模型时，在表层选取 SN向怒江断裂、NW向德钦-维西-红河断裂以及NEE向木里-丽江断裂的作用，断裂切割深度考虑在莫霍面下1 000 m处逐渐消失。除保留原有模型的边界条件外，在模型的西部，大致在本区构造发育的“蜂腰”部位，表层施加代表印度板块自西向东的挤压应力，相应在东部边界的表层加入 EW向约束，如图7-11b所示。计算得到的构造特征与实际资料符合很好。在模型的网格位移图（图7-12a）中可以看出，模型整体除 y 方向有位移外还存在明显的x 方向位移，即 EW向的挤压位移，x 方向位移呈以加力部位为核心向外发散状分布（图7-12b），至东部边界应力逐渐过渡为均匀分布状态。图7-10 第二个模型位移及应力分布计算结果图7-11 第三个计算模型结构图7-12 第三个模型位移计算结果y 方向位移分布特征为整体向北移动，但在深部和浅部相差很大，同时在 SN向断裂的蜂腰部位出现扰动（图7-12c）。z 方向的位移分布较为均匀（图7-12d），只在加力和约束部位出现较大变化。在模型的 x 方向应力图中可以看出（图7-13a），应力在 EW向加载处附近集中，向外发散，断层的存在对发散形态的扰动不大，反映了本区长期经受构造应力作用，其断裂分布形态受到改造，使应力作用趋于均匀。在 y 方向应力分布图中（图7-13b）可以看到 SN向的挤压作用，但在上部趋于减弱。从 x、y 方向的应变图中可以看出与应力分布相似的分布特征（图7-14a、b）。图7-13 第三个模型应力分布计算结果图7-14 第三个模型应变分布计算结果图7-15 第三个模型不同加载阶段应力分布计算结果图7-16 第三个模型位移矢量分布计算结果由于塑性变形与加载时间有关，从不同阶段的 y 方向应力分布图中（图7-15）可以看出，在加载过程中应力分布在中部出现低值区，反映东西两侧应力、应变大于中部盆地发育区的特点，中部地块相对两侧向南挤出。从模型位移矢量分布图中也可以明显看出这一特征，在模型顶部位移矢量图中（图7-16a）可以看到模型整体有向北运动的趋势，但在 EW向受到扰动，在 SN向断裂“蜂腰”处最为强烈。在模型侧面位移矢量图中（图 7-16b）可以看到位移随深度存在明显变化，越趋于上部，位移越弱。

下面以某油田为例说明时移地震数据基于油藏数值模拟解释技术的应用和优点。此油田在1988年、1994年分别进行三维地震数据的采集及处理，而且油藏是1988年底开始生产。两次地震数据的采集方向 非常相近，非常适合进行时移地震研究。该油藏属于上新世的中等砂岩油层，以细砂岩为主，砂体向上尖灭于构造高点的盐丘上，其平均孔 隙度为31％，平均渗透率为500mD。油藏为天然水驱开采，用Gassman方法得知油藏含油饱和度从80％ 降低到20％时，将引起声阻抗11％的变化。在本应用实例中，用数值模拟方法对时移地震数据解释的流程框图见8.22。图8.22 时移地震基于数值模拟解释实际应用流程图8.2.3.1 油藏描述数据的收集由于油藏没有数值模拟模型，所以首先收集了油藏建模所需的各类油藏描述数据。三维储层建模方 面采用随机建模的方法。在建模过程中，要应用地质、物探、测井、测试、开发等多学科数据，通过精 细的地层对比和沉积微相分析，应用地质统计学随机建模的方法以及先进的储层建模理论和技术（如序 贯高斯模拟等），建立工区三维地质模型，包括三维构造模型、三维储层格架模型、三维储层参数（孔 隙度、渗透率、含油饱和度）模型，并进行三维体积计算和模型粗化。因此，收集的数据包括：1）用井标定过的两次采集的地震体；2）深度、厚度解释；3）在测井信息约束下的声阻抗反演；4）测井数据及岩心数据；5）生产数据，包括相渗、PVT信息。8.2.3.2 建立油藏数值模拟模型油藏模型主要包括四个部分，地层模型、构造模型、储层模型和流体模型。其中储层模型是建模的 核心部分。储层模型包括储层格架模型和（渗流）参数模型。储层格架模型即储层的空间展布，由储层埋藏深 度和储层厚度两个变量控制，它可以清楚地展示出目的层段储层的空间立体展布状况；储层（渗流）参 数模型反映了储层储集能力、渗流能力以及孔隙中流体的种类与含量，包括孔隙度、渗透率、含油（气，水）饱和度等。（1）建立储层格架模型（图8.23）图8.23 建立的储层格架模型（2）建立油藏流体（渗流）参数模型油藏流体参数模型包括：(1)从油藏特征中得到的深度、厚度、孔隙度的最优值；(2)从渗透率与孔隙 度或与其他参数的相互关系中得到的油藏绝对渗透率的数值和分布；(3)从试井、测井和岩心分析数据中 得到相对渗透率曲线、压力和饱和度参数值；(4)从油藏工程中得到油藏的压缩系数参数。8.2.3.3 前期模型及在动态数据约束下的数值模拟在收集了以上所有数据的基础上，建立油藏数值模拟的初始模型，模型面元为600ft×600ft，网格划 分为40×65×1（Huang，1997）。对模型进行模拟运算，在运算的过程中，手工对模型进行初步修改，修改集中在对油藏的压缩性和 相对渗透率曲线的调整上，并用从单井测井等方面获取的数据对井点所在处的模型参数进行校准。但这种校准只能集中在近井地带，由于油藏非均质性的存在，井间的油藏参数还是存在不可靠之处。为此，对油藏模拟进行以动态生产数据为目标函数的模拟运算，通过反复修改模型和运算，得到模拟生 产60个月的压力和含水率，拟合曲线见图8.24。图8.24 动态数据约束下的压力和含水率拟合曲线图8.25 在动态数据约束下模拟的含水变化拟合成功后油藏含水饱和度的变化如图8.25所示，含水饱和度的变化显示，经过60个月的模拟生产 后，模型中边水明显向油藏中推进。8.2.3.4生产、地震数据共同约束下的历史匹配在以上工作的基础上，进行了在地震和生产数据共同约束下的历史匹配。双重约束下的模拟循环流 程见图8.26。图8.26 地震和生产数据共同约束下的模拟循环流程图油藏在1988年和1994年分别进行了一次地震测量。因此在合成时移地震属性时，也分别用模型模 拟计算出的1988年和1994年的油藏动态模型，以及相应的油藏静态模型，分别合成1988年和1994年 的地震响应，最后将合成的两次地震响应求差，得到由模型合成的地震差异性属性体。采用模拟计算出的饱和度、压力和Gassmann模型来合成弹性模量，其计算公式为海上时移地震油藏监测技术式中：Ks为岩石颗粒弹性模量；Kd为干岩石弹性模量；Φ为孔隙度；Kf为流体系数，可由式： 计算；Kw，Ko，Kg分别是水、油、气的弹性模量；Sw，So，Sg分别是水、油、气的饱和度。在每一次循环运算中，利用以上给出的公式，合成1988年和1994年的地震差异性属性（图8.27）。图8.27 由初始模型合成的1988年和1994年地震差异性属性的平面分布8.3.3.5 时移地震数据处理首先对1988年和1994年的地震数据进行面元重置、归一化处理、属性提取和差异性分析，得到测 量的地震差异。实际上，流体接触面的运移、各种外界干扰、测量区域的误差都有可能造成测量的地震 数据体的差异。因此，要对不同时间采集到的数据进行归一化处理，使得那些与油气藏无关的地震响应 具有可重复性，保留与油气藏有关的地震响应的差异。将处理得到的实测地震差异性属性与合成的地震差异性属性进行对比。图8.28是由初始模型合成的 时移地震属性与实测数据的对比。根据对比结果进行模型修正，进入下一次模拟循环，直到动态历史拟 合和时移地震数据匹配都满足目标函数的要求，即得到油藏的最优化模型。图8.28 合成时移地震数据与实测数据的对比（据Huang，1998）8.2.3.6 时移地震数据基于数值模拟的解释结果图8.29是由最优化模型得到的1988年和1994年的油藏动态模型合成时移地震差异性属性分布。地 震差异的对比图显示，实测和合成的地震差异基本一致。因此，这两个时间点的油藏模型中参数的变化，就是对时移地震数据的最优解释，并且是结合了动态信息的结果。图8.29 测量地震差异与最优模型合成地震差异对比（据Huang，1998）图8.30反映了该模型含水饱和度参数的变化，包括生产前原始含水饱和度、由最初模型和最优模型 模拟得到的1994年的含水饱和度。由于最优模型即满足动态开发历史又能和时移地震数据进行很好的匹 配，因此，由最优模型得来的1994年的含水饱和度分布的置信度是比较高的。图8.30 在时移地震数据约束下得到的油藏参数动态变化（据Huang，1998）8.2.3.7 不同目标函数约束下最优模型对比由于约束历史拟合的目标函数不同，最终得到的最优模型会存在较大差异（图8.31）。比较证明，仅用动态数据或仅用地震数据做目标函数得到的模型拟合含水的误差都要大于用动、静态数据共同约束 下得到的含水拟合误差。图8.31 不同约束下得到的最优含水拟合曲线比较（据Huang，1998）8.2.3.8 结论及认识1）比较证明，用数值模拟的方法对时移地震数据进行解释不但是可行的，而且由于有了时移地震数 据作历史拟合的目标函数，极大地提高了历史拟合的精度和模型调整的可信度。2）本次时移地震技术的应用还证明：用随机最优化的方法缩小测量数据和合成数据差异的方法是有 效的。3）时移地震方法可以用来进行油藏动态描述，用时移地震方法监测油藏流体运动的准确度要超过许 多其他方法。

数值模拟技术诞生于1953年Bruce G.H和PeacemanD.W模拟了一维气相不稳定径向和线形流。受当时计算机能力及解法限制，数值模拟技术只是初步应用于解一维一相问题。两相流动模拟诞生于1954年，West W J和Garvin W.W模拟了油藏不稳定两相流。1955年Peaceman与Rachford研发的交替隐式解法(ADI)是数值模拟技术的重大突破。该解法非常稳定,而且速度快，所以迅速在包括石油,核物理，热传导等领域得到广泛应用。1958年Douglas,Jim和Blair,P.M第一次进行了考虑毛管压力效果的水驱模拟。1959年，Douglas Jim和Peaceman D.W第一次进行了两维两相模拟，这标志着现代数值模拟技术的开始。在他们的模拟器里全面考虑了相对渗透率，粘度，密度，重力及毛管压力的影响。60年代数值模拟技术的发展主要在数值解法，第一个有效的数值模拟解法器是1968年Stone推出的SIP(Strong Implicit Procedure)。该解法可以很好地用来模拟非均质油藏和形状不规则油藏。另一个突破是时间隐式法，该方法可以用来有效的解高流速问题，比如锥进问题。60年代其他方面的发展还有1967年Coats K.H和Nielsen R.L首次进行了三维两相模拟，而且提出了垂直平衡和拟相对渗透率及毛管压力方法。1968年Breitenbach E.A发表了三维三相模拟解法。Stone在70年代发表了三相相对渗透率模型，由油水和油气两相相对渗透率计算油、气、水三相流动时的相对渗透率，该技术现在还广为应用。70年代另一项主要成就是Peaceman提出的从网格压力来确定井底流压的校正方法，及现在通用的Peaceman方程。在解法方面的发展是采用了正交加速的近似分解法。70年代在组分和热采模拟方面也取得很大进展，1973年Nolen J.S描述了考虑油气中间组分分布的组分模拟，Cook提出变黑油模拟来进行组分模拟。Shutler在1970年发表了对两维三相模型的蒸气注入模拟。70年代在EOR方面也取得了极大进展。80年代最大的成就是Appleyyard J R和Cheshire I.M发表了嵌套因式分解法，该解法非常稳定而且速度快，是目前最为广泛应用的解法。正是基于该解法，Cheshire I.M于1981年同John Appleyard和Jon Holmes成立ECL公司，开始研发后来主导数值模拟软件市场的ECLIPSE软件。80年代见证的另一个主要发展是组分模型，虽然组分模型在60年代就已经推出，但很不稳定。80年代提出的体积平衡和Yong-Stephenson方程解决了组分模型稳定问题，使组分模型可以广为应用。Ponting D.K提出了角点网格来模拟模型，这样可以真实地描述油藏。90年代数值模拟的进展主要在粗化技术，并行计算，PEBI网格等方面。Zoltan E.Heinemann提出了PEBI网格，PEBI网格结合了正交网格和角点网格的优点，现在正逐渐成为主流数值模拟网格体系。VIP于1994年推出并行算法，ECLIPSE于1996年推出并行算法。CMG于2001年推出并行算法。粗化技术的难点在于渗透率的粗化，基于流动计算进行的渗透率粗化可以较真实的符合地质模型，现在新的粗化技术还在发展。　 21世纪数值模拟技术发展体现在两方面，一方面是一体化模拟技术，数值模拟将不只是对油藏的模拟，数值模拟将对油藏，井筒，地面设备，管网以及油气处理厂进行一体化模拟，从而最优化管理油田。另一方面是定量进行属性不确定性分析，定量分析属性不确定性对计算结果的影响。

通过对滑坡变形破坏现场的地质分析和相似材料模型模拟试验研究，对滑坡变形破坏及滑坡对管道的作用机理有了一定的认识。为了更深入地研究滑坡与管道相互作用的机理，验证前述认识，采用有限元计算方法对管道沿线典型滑坡进行了数值模拟研究。5.4.1　研究方法对于典型滑坡的滑体、滑动面、滑床、滑坡周界等可以用有限元方法进行离散，着重考虑滑坡潜在滑动面的摩擦接触情况，采用摩擦接触力学方法，对接触单元进行改进，建立接触力学模型，对滑坡进行研究分析。在考虑可能滑坡面的法向抗张强度和切向抗剪强度下，通过多步迭代计算的方法，获得非线性接触问题的收敛解。通过接触力学计算分析，模拟滑坡的孕育；从滑动面的接触状态的改变（粘着、滑动、未接触状态），获得滑坡的滑动情况，滑坡的动态过程（启动、滑动、停止），同时还进行坡体和管道的应力、变形，以及边坡滑面接触关系分析。在这一过程中重点研究管道在滑坡发生过程中受力破坏的力学机理。5.4.2　数值模拟5.4.2.1　改进接触单元数值模拟研究中，接触单元的法向与切向模量不一样的，法向模量EN可由两接触物间不可能相互挤入这一事实来确定。在三维情况下，如取N表示接触单元的法向、T1表示切向方向1、T2表示切向方向2，则可推导出：（1）法向模量EN可由接触物间不可相互挤入的条件确定；（2）而模量 、剪切模量 ；（3）在切向平面内的剪切模量 可由接触状态（粘着、滑动等）确定。于是可确定接触单元的弹性矩阵[D]，进而可获得三维接触单元的刚度矩阵为：山区油气管道地质灾害防治研究式中：[B]是接触单元的几何矩阵。对包含潜在滑面的整个边坡结构进行有限元离散，可得到边坡结构的刚度矩阵：山区油气管道地质灾害防治研究5.4.2.2　潜在滑面的强度对于未出现任何破坏的边坡，其潜在滑面是具有一定强度的（实际上滑面还未形成。一般情况下是具有法向强度 和切向强度τ，的，一般情况下 。及 。法向强度的数值可取为边坡潜在滑面处的材料抗张强度，而对于切向强度通常在岩土力学中是取Mohr-Coulomb模型，认为切向强度τs与法向压应力之间σN是一个线性关系：山区油气管道地质灾害防治研究式中：c是边坡潜在滑面处的内聚力；ψ是内摩擦角。由上式可确定边坡潜在滑面处的切向强度。而在滑坡发生后，潜在滑面已经变为真实滑面，此时真实滑面位置退化为纯摩擦接触问题。5.4.2.3　斜坡潜在滑面的接触状态在滑坡发生前，潜在滑面处于完好的结合状态，在斜坡自重及外部条件的作用下，接触状态会改变。非接触状态；当接触单元的法向应力是张应力，且该张应力大于潜在滑面处材料的抗张强度，即 ，此时，滑面将被拉开。接触单元的法向间隙会大于零，即Zg>0。接触单元的模量需取为一个很小的值EN=10-6～10-4Ee,Ee是接触单元的初始模量山区油气管道地质灾害防治研究当接触单元的法向应力是张应力，且当接触单元的切向应力τ大于潜在滑面材料的剪切强度τs时，潜在滑面将会破坏，此时Zg>0且EN=10-6～10-4Ee。山区油气管道地质灾害防治研究粘着状态：当接触单元的法向应力为压应力，且接触单元的剪应力小于潜在滑面处的抗剪强度时，此时Zg≤0。这意味着接触单元的法向隙Zg是一个很小的近似于零的值山区油气管道地质灾害防治研究当接触单元的法向应力为张应力，且接触单元的剪应力小于潜在滑面处的抗剪强度τs，此时山区油气管道地质灾害防治研究滑动状态：当接触单元的法向应力为压应力，且接触单元的剪应力大于潜在滑面处的抗剪强度时，此时潜在滑面将破坏，变为真实滑面，而Zg≤0。山区油气管道地质灾害防治研究由于剪应力的作用，潜在滑面变为真实滑面，滑坡发生了。接触单元的剪切模量 由下式决定：山区油气管道地质灾害防治研究5.4.2.4　非线性迭代过程的收敛计算接触状态的最终确定，需满足以下条件：非接触状态应满足Zg>0且EN=10-6～10-4Ee，而接触状态（粘着、滑动）应满足Zg=0（在指定精度ε下）。每两次迭代循环计算时，接触单元模量的改变应满足山区油气管道地质灾害防治研究其中ΔE、ΔG是相对的模量改变量，ΔEp、ΔGp是许可的相对模量改变量，山区油气管道地质灾害防治研究对于滑动状态的接触单元，当下式满足时，滑动才处于稳定状态：山区油气管道地质灾害防治研究在以上3个条件同时满足时，该接触问题的计算才是收敛的。5.4.2.5　滑坡过程的逐步计算在边坡滑动过程的数值模拟中，使用了逐步计算的方法。用一系列时间步t1,t2，…，tK，…，tn来分解边坡的滑动过程。而对其中的某一个时间步tK，根据边坡下滑的加速度，对边坡的滑动部分施加相应的惯性力，所以，边坡将保持平衡状态。在某一个时间步tK时，滑坡系统的势能改变ΔΠ可写为：山区油气管道地质灾害防治研究其中：山区油气管道地质灾害防治研究在tk时间步，上式中ΔWK是体积力bi和边界力Fs（不包含接触边界—滑面）所做的功，重力势能的改变量为ΔwE,ui为边坡系统的位移量，ΔUK是应变能改变量，σij和εij是应力和应变，Ω为边坡系统的总体积，ΔAK为滑动面上接触应力所做的功，p。和pt为接触面上的法向接触压力和切向力，r为接触面 和 分别为接触面的上、下边界， 和 分别为接触面上的上、下边界的切向位移，ΔKK为动能改变量。由于问题的复杂性，研究中未考虑材料破坏所消耗的能量。在时间步从t1=0到tn计算完成后，可获得整个滑坡动态过程的应力、应变、接触状态改变情况、滑动距离等重要参数。

开放边界。求采纳

从上面的例子可以看到，数值模拟包含以下几个步骤：首先要建立反映问题（工程问题、物理问题等）本质的数学模型。具体说就是要建立反映问题各量之间的微分方程及相应的定解条件。这是数值模拟的出发点。没有正确完善的数学模型，数值模拟就无从谈起。牛顿型流体流动的数学模型就是著名的纳维—斯托克斯方程（简称方程）及其相应的定解条件。数学模型建立之后，需要解决的问题是寻求高效率、高准确度的计算方法。由于人们的努力，目前已发展了许多数值计算方法。计算方法不仅包括微分方程的离散化方法及求解方法，还包括贴体坐标的建立，边界条件的处理等。这些过去被人们忽略或回避的问题，现在受到越来越多的重视和研究。在确定了计算方法和坐标系后，就可以开始编制程序和进行计算。实践表明这一部分工作是整个工作的主体，占绝大部分时间。由于求解的问题比较复杂，比如方程就是一个非线性的十分复杂的方程，它的数值求解方法在理论上不够完善，所以需要通过实验来加以验证。正是在这个意义上讲，数值模拟又叫数值试验。应该指出这部分工作决不是轻而易举的。在计算工作完成后，大量数据只能通过图像形象地显示出来。因此数值的图像显示也是一项十分重要的工作。目前人们已能把图作得像相片一样逼真。利用录像机或电影放映机可以显示动态过程，模拟的水平越来越高，越来越逼真。

物理模型->数学模型（控制方程和边界条件）->数值离散方法->网格剖分->程序编制->设定参数->结果分析参考：http://blog.sciencenet.cn/blog-577059-538449.html

你说的是哪方面的数值模拟呢？问题不是很明确啊

数值解是在一定条件下通过某种近似计算得出来的一个数值，能在给定的精度条件下满足方程解析解为方程的解析式（比如求根公式之类的），是方程的精确解，能在任意精度下满足方程

必须验证。数值经工作人员模拟出来以后，还要检查它的正确性。验证方法就是把数值进行倒着运算，看一看结果是不是最初值，简单又快捷。数值模拟也叫计算机模拟。依靠电子计算机，结合有限元或有限容积的概念，通过数值计算和图像显示的方法，达到对工程问题和物理问题乃至自然界各类问题研究的目的。

（不是专业人士,说说想法……） 这种问题提出,大约是在工程中遇到了求解数学问题的情况（偏微,常微等） 解析分析——用数学分析的方法,比如微分、积分、特殊方程等,对实际情况进行模拟,列出方程,用解析的方式,求出“比较正规”的函数解（这样可以解决的问题是很有限的,因为在实际问题中,有很多是不能用函数简单模拟的；而且这样的方式,解决方法灵活多变,不利于机器模拟——虽然有软件可以做到） 数值模拟——以不用具体的函数表达式,而是用多个点的数值表示函数的方法,来解实际问题的解法,恩,一般来说,差分法是最常用的（这样做的在于应用范围很广,但是计算误差必须估计） 盼有所帮助. 补充一点：另外一种可能,解析分析就是用解析的方法求解,数值模拟用于检验.

数值模拟也叫计算机模拟。它以电子计算机为手段，通过数值计算和图像显示的方法，达到对工程问题和物理问题乃至自然界各类问题研究的目的。　　在计算机上实现一个特定的计算，非常类似于履行一个物理实验。这时分析人员已跳出了数学方程的圈子来对待物理现象的发生，就像做一次物理实验。　　数值模拟实际上应该理解为用计算机来做实验。比如某一特定机翼的绕流，通过计算并将其计算结果在荧光屏上显示，可以看到流场的各种细节：如激波是否存在，它的位置、强度、流动的分离、表面的压力分布、受力大小及其随时间的变化等。通过上述方法，人们可以清楚地看到激波的运动、涡的生成与传播。总之数值模拟可以形象地再现流动情景，与做实验没有什么区别。　　从上面的例子可以看到，数值模拟包含以下几个步骤：　　首先要建立反映问题（工程问题、物理问题等）本质的数学模型。具体说就是要建立反映问题各量之间的微分方程及相应的定解条件。这是数值模拟的出发点。没有正确完善的数学模型，数值模拟就无从谈起。牛顿型流体流动的数学模型就是著名的纳维—斯托克斯方程（简称 方程）及其相应的定解条件。　　数学模型建立之后，需要解决的问题是寻求高效率、高准确度的计算方法。由于人们的努力，目前巳发展了许多数值计算方法。计算方法不仅包括微分方程的离散化方法及求解方法，还包括贴体坐标的建立，边界条件的处理等。这些过去被人们忽略或回避的问题，现在受到越来越多的重视和研究。　　在确定了计算方法和坐标系后，就可以开始编制程序和进行计算。实践表明这一部分工作是整个工作的主体，占绝大部分时间。由于求解的问题比较复杂，比如 方程就是一个非线性的十分复杂的方程，它的数值求解方法在理论上不够完善，所以需要通过实验来加以验证。正是在这个意义卜讲，数值模拟又叫数值试验。应该指出这部分工作决不是轻而易举的。　　在计算工作完成后，大量数据只能通过图像形象地显示出来。因此数值的图像显示也是一项十分重要的工作。目前人们已能把图作得像相片一样逼真。利用录像机或电影放映机可以显示动态过程，模拟的水平越来越高，越来越逼真。

1、首先模拟地震波在不同土层中的传播过程。2、其次对不同的地震波在不同土层中的传播过程的数值进行统计。3、最后根据统计数据进行对比得出即可。

海洋环流模拟，海洋波浪模拟，海洋污染模拟，海洋气候模拟。1、海洋环流模拟：数值模拟可以对海洋中的环流进行模拟和预测，以便更好地了解海洋环流的形成机制、演变规律和影响因素。2、海洋波浪模拟：数值模拟可以对海洋中的波浪进行模拟和预测，以便更好地了解海洋波浪的形成机制、传播规律和影响因素。3、海洋污染模拟：数值模拟可以对海洋中的污染物输运进行模拟和预测，以便更好地了解污染物的扩散规律和影响因素，为海洋环境保护和治理提供科学依据。4、海洋气候模拟：数值模拟可以对海洋中的气候变化进行模拟和预测，以便更好地了解气候变化的机制、趋势和影响因素，为应对气候变化提供科学依据。

近20年来计算速度和功能在不断改进，新的模拟技术也在不断发展，但无论如何，数值模型是以地质模型为基础的，它通过网格化将一个地质体雕刻出来，网格分得越细，网格节点越多，刻画得就越逼真，计算的工作量也越大。就像数码照相机一样，相素越高，照出来的照片越逼真，计算机处理的容量也越大。数值模拟受到两方面限制。常规模拟与随机模拟两种方法的比较一是受到对油藏认识的限制。任何一个地质体，其在空间上任何一个方向都有着非均质性。有些参数（如含油饱和度、毛细管压力）至今人们还不能充分认识其内在规律，而目前所有的探测手段的精度还远远不能满足我们的要求，因此，我们使用数值模型方法对油藏进行研究所得到的结果还只能说是相对的，或者说是半定量的，是某种尺度上平均的结果。针对这种情况，更需要油藏地质师和油藏工程师大力合作才能取得较好的效果。过去常规的模拟方法往往是油藏地质师把地质模型交给油藏工程师后就基本上不管了，剩下的问题是油藏工程师去建立数值模拟的模型并进行历史拟合，拟合不好往往也是由油藏工程师去修改数值模拟的模型，要么修改相对渗透率，要么修改毛细管压力等参数，致使一些实验数据面目全非。实际上，除了井点上的一些确定性的数据不能修改以外，其他网格节点上的插值数据都是随机的，拟合不好不单是油藏工程的数据问题，更重要的是地质模型的非均质性如何处理，所以，油藏地质师和油藏工程师大力合作十分重要。按照随机建模的理念，油藏地质师和油藏工程师都应该共同研究，分工合作，反复修改地质模型和数值模拟的模型，以更逼近油藏静态的和动态的实际情况。二是受到计算机的限制，计算机处理的容量越大，计算速度越慢、成本越高。过去超百万节点模型运算很难实现，直到近十年来才有了飞快的发展，由于采用了并行处理技术，加快了速度，使得过去串行处理技术难以实现的问题得到了解决。这些新技术包括区域分解计算理论；大型稀疏线性方程组的预条件求解技术；迭代加速计算方法；并行数据流结构设计及处理器之间的数据通讯技术。尽管如此，大容量超百万节点的模拟无论在硬件或软件上还不成熟，我们不可能不论油田大小、不分研究问题的复杂程度，把数值模拟看作神乎其神，甚至顶礼膜拜。要用好数值模拟，必须从油气田实际出发，针对研究的问题，多做些机理性的研究和预测，模型的选择能用二维的模型解决问题不必非要用三维的模型，能用两相的模型解决问题不必非要用三维三相的模型。模型不在于如何复杂和多大容量，关键在于能否代表油气藏的特性，妙算在于人，在《汉书》颜师古注：“务得事实，每求真是也”。使用数值模拟时，如果要取得更符合实际的结果，开发家和油田地质家对油藏及其各种参数有更深入的了解这才是最重要的。对油藏了解得越深刻，提供的参数越接近实际，生产的历史数据拟合越好，用油藏数值模拟得到的结论就越接近实际情况。百万网格节点测试结果比较（Power Challenge并行机6CPU）

有限单元法是应用于构造应力场模拟的最广泛的数学模拟方法。其基本思想是将所研究的地质体以一定的方式（单元形状和节点个数）简化为有限个单元组成的离散化模型，再用相应的计算程序求出数值解答。利用有限元法数值模拟，可以利用地质调查和构造解析获得的较少地质应力状态的资料来反演区域内各点的应力状态，从而获得区域的构造应力场特征，加深认识区域内的构造演化。目前有限单元法的应用已由弹性力学的平面问题扩展到空间问题、板壳问题，分析对象从弹性材料扩展到塑性、粘弹性、粘塑性和复合材料。有限元法数值模拟随着计算机技术的发展在科学计算领域得到广泛应用，20世纪80年代以来，国际上已有较大型的有限元计算程序达几百种，其中较著名的有：ANSYS、NASTRAN、ASK、ADINA、SAP等。以ANSYS为代表的数值模拟软件将有限元分析、计算机图形学和优化技术相结合，已成为科学计算领域不可缺少的有力工具。基于本区岩石圈的三维结构特点，我们首先对本区的三层结构相互作用关系进行了模拟。对本区的物理模拟研究，前人已经做过很多工作，其中在对印度板块挤压下亚洲中东部的构造模拟中，有的反映出大型走滑断裂、裂谷和张性盆地以及压性逆冲断裂等构造现象，有的反映出多层构造中网络状流动现象，认为板内变形受塑性流动网络控制（Tapponnier et al.，1982；李建国等，1997）。这些工作往往只反映了本区的某一方面的特性，而无法对本区的构造形态做出动力学的完善解释。因此在前人的工作基础上，我们首先建立了本区的一个三层结构模型，其中中上地壳深度根据天然地震资料定为30 km，下地壳以莫霍面为其底界，根据地震测深资料取50 km。因为本模型建立的主要目的是确定岩石圈各圈层之间的作用关系，因此模型底部只考虑到100 km的深度。

2137西下工作面上方为1237西下面老塘,该工作面于2002年回采结束。在回采过程中曾有过最大1.6m3/min的涌水,目前上水平12煤层西正洞仍有少量渗水现象,经分析1237西下面老塘内仍有积水。防止上工作面老塘水突入,防水煤柱留设应兼考虑上工作面底板导水裂隙带深度与下工作面顶板导水裂隙带高度。依据数值模拟计算结果,在开切眼处,煤层抽冒最大高度为35.1m；1237工作面煤层及围岩条件与2137类似,用2137工作面底板最大导深近似上工作面底板导深,为15.2m；合计为50.3m。在远离开切眼处,煤层抽冒不发育,顶板最大导高为17.9m,与上工作面底板最大导深合计为33.1m。因此,无论在近开切眼或远离开切眼的位置,60m煤柱足以形成安全防水煤柱。开滦矿区依据其煤层顶板岩石物理力学性质特点,煤层顶板皆属于中硬强度,依据煤矿《矿井水文地质规程》附录七中冒落带与导水裂隙带最大高度的经验公式,裂隙带发育高度为：当煤层倾角为55°～85°时,Hf=100hM煤层开采顶板导水裂隙带高度预测理论与方法式中：h为小阶段垂高；M为煤层厚度；Hf为导水裂隙带高度。依照该公式进行计算,2137工作面最大导水裂隙带高度为75.68m。依据经验公式计算,在2137西下面回风道距1237老半道60m垂高作为防水煤柱,不足以作为安全的防水煤柱。2137西下工作面自2006年1月至2007年12月实现安全回采,60m防水煤经实践检验是足够安全的。这说明与经验公式计算结果相比,采用FLAC3D数值模型确定的顶板导水裂隙带高度更符合实际。开滦矿区在急倾斜煤层开采实践中,形成挡水煤柱留设的经验为：采一层(约30m)留两层(约60m)。数值模拟的结果,理论上支持了这一经验认识。综合以上分析,得出以下几点结论：(1)采动影响下,冒落带、裂隙带形成过程可视为：似连续介质条件下覆岩变形破坏、非连续介质条件下自重作用下的冒落 裂隙扩展、冒落带压密 裂隙扩展三阶段。采FLAC3D软件模拟,将岩体视为大变形连续介质,可以实现冒落带、裂隙带最终形态的近似模拟。走向开挖15m开采后,冒落带碎胀系数1.3；走向开挖30m后,冒落带碎胀系数1.15。工作面顶板冒落带、导水裂隙带在走向开挖15m后,达到其最大值,位于切眼处。(2)采动影响下,急倾斜煤层工作面围岩由内向外产生压缩变形区—扩容区—压缩区。内部压缩区可视为冒落带；依据体积应变1×10-2等值线为导水裂隙带界线；扩容区内局部压缩区视为煤层潜在抽冒区。走向开挖30m后,冒落带达到其最大高度10.4m,煤层潜在抽冒区产生抽冒,当其高度超过体积应变1×10-2等值线时,作为导水裂隙带最大高度。则导水裂隙带最大高度为35.1m。(3)开切眼处2137工作面顶板煤层抽冒带与上部1237工作面底板导水裂隙带构成2137工作面上方透水区,合计垂高50.3m；远离开切眼处,2137工作面顶板导水裂隙带与上部1237工作面底板导水裂隙带共同构成2137工作面上方透水区,合计垂高33.1m。2137工作面设计60m垂高防水煤柱,实现安全回采,说明该数值模拟计算结果的合理性。

算实验，是一种计算机数学试验！

不一样吧……似乎后者包括的范围更大一些……

通过Mechanical APDL Product Launcher 19.0，对轴承座铸件进行浇注结束后，凝固冷却这段时间的温度场、应力场数值模拟，可以发现轴承座在凝固冷却成型阶段温度最高、应力最大、最容易失效的部位还是轴承座的底座。所以综上所述，应当降低底座的厚度，提高其散热性能，避免因温度过高而引起的缩孔缩松，同时增大底座的长度宽度，目的是为了使得螺栓孔的间距变大，从而减少应力集中。必要情况下，可以采取自然时效、热时效、振动时效等方式消除应力，但是自然时效周期长不建议使用，热时效成本太大并且能耗不低，所以可以采用振动时效设备，对铸造过程中的等效应力进行一定量的清除。本次研究旨在通过结果分析，来探明轴承座铸件在凝固冷却成型过程中温度场和应力场的一个变化的情况，并根据这种变化，对铸件本身的结构加以一定的改进，以此来提升轴承座铸件的表面质量、结构精度、工作效率。因此，可以依据本次研究的结果，在以后轴承座的铸造生产过程中，加以运用，以此来优化轴承座凝固冷却的效果，具有一定的意义。

1、首先，需要掌握地球科学、物理学、数学等多个学科的基础知识。2、其次，通过实践操作和练习，增强对于数值模拟方法和工具的熟练掌握，为解决实际问题打下基础。3、然后，参与科研项目或团队可以帮助我们更深入地了解地球动力学数值模拟的应用和发展趋势，并与其他领域的学者进行交流和合作。4、最后，持续关注最新技术和进展，不断更新自己的知识和技能。

一、水文地质概念模型（一）模拟范围黑河流域中游地下水流场模拟范围，东起民乐总寨-山丹祁家店水库，西至嘉峪关市的河口-吕家庄（嘉峪关大断裂），夹峙于祁连山、龙首山、合黎山、北部夹山和南山之间的走廊平原，主要包括张掖盆地和酒泉东盆地两部分，面积约11300 km2。（二）地下水系统结构概化及边界条件研究区是一个只有侧向流入而没有侧向流出的闭型山间断陷水文地质盆地，其间填充的巨厚第四纪松散沉积物，构成连续的、统一的、相互叠置的、横向为盆地边界所限的地下水系统。周边地质构造或山体为天然地质边界。在研究区东部张掖盆地，地下水自南东向北西运动，最终排泄于黑河干流，流出区外。在研究区西部酒泉东盆地，地下水自南西向北东运动，部分排泄于蒸发，部分排泄于黑河干流，流出区外。区内潜水系统接受大气降水的入渗补给。在南部洪积扇地带接受河水和引灌渠水垂向淋滤渗漏补给。在中游冲积细土平原接受沟渠、水库和灌溉水入渗补给。在洪积扇扇缘和与之毗邻冲积细土平原带，地下水以泉的形式大量溢出而排泄。在冲积细土平原黑河河床地带，地下水以泉的形式集中排泄于黑河。在中游盆地北半部地下水埋藏较浅区，地下水以蒸发方式排泄。在人口稠密、农业集中开发区，人工开采是一种重要的地下水排泄方式。在与洪积扇毗邻的地带，承压水系统得到了扇群带地下水侧向径流补给。在人口集中居住的地区和农业集中开发区，人工开采是一种重要排泄方式。在黑河流域中游山前冲洪积扇及黑河河床地带皆为单一的潜水层，而在冲积细土平原地带则逐渐过渡为双层结构或多层结构的潜水-承压水统一体系。在张掖与平原堡之间的黑河河床在140 m之内都是透水性较强河流冲积物，切断了冲积细土平原北半部承压水之间的联系，从而使张掖与临泽形成两个各自独立的承压水地段。在具有多层结构含水层地带，由于下部的多层承压含水层之间隔水层较薄，且不连续，二者之间具有一定的水力联系，因此，可将下部承压含水层组视为统一含水系统。在冲积细土平原区的潜水与承压水系统之间存在0.9～10.5 m的粘性土弱透水层，承压水水位一般高于潜水，但是在张掖地区的古城、碱滩附近潜水水位高于承压水水位，二者之间存在一定水位差，含水层间存在越流水量交换。对于洪积扇扇群带、黑河河床地带等单一潜水层含水层地带，可在潜水系统中人为划定一薄层“弱透水层”，假设“弱透水层”的渗透系数很大，这样就可以把整个流域含水层结构概化为潜水-承压水双层结构。在黑河流域中游南部的祁连山区、东部的民乐总寨-山丹祁家店水库一带，为二类透水流量边界。在西部的嘉峪关大断裂一线、南部的龙首山、合黎山和北部的北山一带，为二类弱透水流量边界。在张掖盆地与酒泉东盆地之间分界处，虽然两侧地下水运动方向不同，但是可近似为连续含水系统，只是导水性差。经过以上分析，黑河流域中游区地下水系统被概化为具有大气降水、河渠水、灌溉水入渗补给，通过泉、河流、蒸发和人工开采排泄，为第二类边界条件，含水层间通过越流联系的双层结构拟三维流含水层系统。地下水流的数学模型可用式（9-3）表示。（三）交换量处理1.河流入渗补给量主要分布于山前洪积扇地带，补给潜水系统。由于河水水位远高于地下水水位，河水以铅直下渗（淋滤渗漏）的方式大量补给地下水，入渗量的大小取决于河流流量、河床的岩性和地形地貌条件。计算入渗量时，参考了甘肃省第二水文地质队1967、1985年实测资料，并在模拟过程中给予识别和确定。2.灌溉水入渗补给量在黑河流域中游农业集中开发区，存在大量灌溉水的入渗，补给量大小取决于灌溉田地的岩性、潜水水位埋深和灌溉水量的大小。补给量按面状补给处理，在模拟过程中识别和确定。3.渠系水入渗补给量主要指黑河流域中游区内的引灌干渠渗漏补给，入渗量的大小由渠系衬砌类型和输水量的大小决定。较小渠系的入渗量，在灌溉水入渗中一起考虑。4.大气降水、凝结水入渗补给量入渗量的大小取决于降水量、潜水水位和地表浅层岩性等因素。据甘肃省第二水文地质队地渗仪实测资料表明，在荒漠区，只有地下水水位埋深小于2 m的地段才接受到降水入渗有效补给，而在灌区，地下水水位埋深小于5 m的地段，也观测到降水入渗有效补给。对于凝结水渗入补给，主要发生在地下水水位埋深小于2 m的地段。该项补给量按面状补给处理，由模型根据分区识别，确定各区入渗系数。5.潜水蒸发量潜水蒸发量大小主要与包气带岩性、地下水水位埋深、植被发育和气候干湿冷暖变化等因素有关，该量可按面状排泄处理。根据甘肃省第二水文地质队地渗仪观测结果，在黑河流域中游区潜水蒸发极限深度为6.5m，潜水年蒸发强度（qe）与地下水水位埋深（g）之间关系曲线如图9-1所示，关系式为西北内陆黑河流域水循环与地下水形成演化模式6.越流量根据弱透水层的岩性、厚度和上下层间地下水水位差条件，由模型识别，确定分区越流系数。7.开采量主要分布在人口集中地区和农业集中开发区，可按点井和面积井处理。8.泉水溢出量图9-1 黑河流域中游区潜水年蒸发强度与地下水水位埋深关系泉水溢出量与泉水的出露高度、地下水水位差呈正比关系，主要参考1984、1986、1999和2001年泉水调查资料。二、数值模拟与结果（一）计算区域剖分选用三角形网格进行剖分，剖分时充分考虑含水层岩性及富水性、河流的分布、水位观测孔位置、区域研究程度和土地开发利用状况、开采井位置等因素。在单层地下水系统中，共剖分成431个节点，747个单元（图9-2）；在双层地下水系统中，共剖分成862个节点，1494个单元，每层各有111个二类流量边界。图9-2 黑河流域中游区地下水系统单元剖分图在进行三角形网格剖分时，采用垂向切开的方法，即每个三角形对垂向两层都适用，包括三角形单元的编码，但是各节点的编码在两层中不同。网格剖分只对第一层潜水系统的三角形单元和节点编码。根据这个信息，计算机通过程序控制对下层承压水系统自动地进行新的一套编码。这样，对下部承压水系统数值模拟计算时，采用的是与上部潜水系统剖分完全不同的新编码。（二）边界条件黑河流域中游区地下水系统边界都是二类流量边界。根据调查资料、裂隙发育程度状况和现今山区降水情况，分段进行计算祁连山区基岩裂隙水的流入量。对于龙首山和合黎山山区基岩裂隙水的侧向径流流入量，按水利部门提供的径流深度，参照该区基岩裂隙垂向发育程度资料分段计算。对于夹山和南山基岩裂隙径流的流入量，参照合黎山资料进行处理。在嘉峪关大断裂，侧向流入量是参考酒泉盆地讨赖河水源地的勘察资料确定。（三）模拟时期选择模拟时期为1987年9月至1988年8月，共分为12个时间段，每个时段的步长均为1个月，潜水以计算区内24个观测孔的地下水水位作为模型识别的主要依据，承压水以计算区内21个观测孔的地下水水位作为模型识别依据。（四）参数分区与模型识别黑河流域中游模拟区的水文地质参数估值，采用分区法确定。根据区内水文地质勘探和现场抽水试验资料，通过数值模拟，将计算区的潜水系统和承压水系统参数初步分区，并假定每个分区具有相同的参数。水文地质参数主要有：①降水入渗系数以及地表水单长渗漏率；②潜水系统给水度与渗透系数；③承压水系统渗透系数与储水系数；④弱透水层越流系数。模型识别要求符合下列标准：①地下水模拟流场特征与实际地下水流场基本一致；②地下水模拟动态变化过程与实际地下水动态变化过程基本相似（图9-3和图9-4）；③地下水模拟均衡变化与实际地下水均衡变化相符；④识别的水文地质参数符合实际水文地质条件。图9-3 模型拟合阶段黑河流域中游区潜水系统初始流场平面图图9-4 模型拟合阶段黑河流域中游区承压水系统初始流场平面图模型的识别采用解逆法反求参数。即将参数的初值带入模型，应用上节所述的针对多层结构地下水系统数学模型，采用有限元法和波前法计算出潜水和承压水的地下水水位，将观测孔的计算水位与实测水位进行比较，其拟合程度采用评价函数来度量。评价函数的表达式为西北内陆黑河流域水循环与地下水形成演化模式式中：M——时段总数；N——观测孔总数；Wj——各观测孔的权函数；———i时段末j号观测孔的计算水头；———i时段末j号观测孔的实测水位。用试估校正法调整参数，并且根据不同时段对不同参数进行识别。当目标值E达到尽量小时，则认为达到了要求。识别后的分区参数如表9-1～表9-6和图9-5～图9-9所示。表9-1 黑河流域中游区潜水系统渗透系数分区参数表9-2 黑河流域中游区潜水系统给水度分区参数表9-3 黑河流域中游区降水入渗系数分区参数表9-4 黑河流域中游区承压水系统渗透系数分区参数表9-5 黑河流域中游区承压水系统储水系数分区参数表9-6 黑河流域中游区弱透水层越流系数分区参数图9-5 黑河流域中游区潜水系统渗透系数分区图图9-6 黑河流域中游区潜水系统给水度分区图图9-7 黑河流域中游区承压水系统渗透系数分区图表9-7 拟合阶段黑河流域中游区潜水系统水量均衡表（108m3/a）表9-7和表9-8表明，模型拟合期的潜水和承压水都是负均衡，地下水水位总体呈下降态势。由于基础资料序列连续性较差，特别是在酒泉东盆地的东部，可供利用的资料少，以至数值模拟结果存在一定误差。但是从拟合结果来看，控制区的观测孔地下水水位的计算值与观测值比较接近（图9-10），无论从整个拟合校正期的过程来看，还是从整个流场来看，拟合效果都较为理想。因此，可以认定所建的模型能够较好地反映黑河流域中游区地下水系统变化规律，服务于不同开采方案条件下地下水动态变化的预测研究。表9-8 拟合阶段黑河流域中游区承压水系统水量均衡表（108m3/a）图9-8 黑河流域中游区承压水系统储水系数分区图图9-9 黑河流域中游区弱透水层越流系数分区图图9-10 黑河流域中游区部分观测孔地下水水位拟合曲线三、地下水水位动态变化模拟在以上研究的基础上，利用求得的水文地质参数，对流域中游未来不同开采条件下的地下水水位动态进行预测研究。预测方案有两种：（A）现时开采方案：不改变现有开采条件，地下水开采量保持在1987年9月至1988年8月期间水平；（B）规划开采方案：不改变现有开采井布局和其他水文地质条件，地下水（潜水和承压水）开采量增加1倍。预测时，以1999年9月初潜水和承压水系统等水位线作为初始流场（图9-11和图9-12），时间步长为1年，预测至2009年9月初潜水和承压水系统地下水水位。方案A预测结果如图9-13和图9-14所示。从图9-13和图9-14可见，预测期内各时段都为负均衡。在张掖盆地东南部的民乐、六坝、毛城子和安阳一带，潜水水位基本保持不变。在张掖市和平原堡一带，潜水水位略有上升，附近河段地下水溢出量略有增加。其他地段，潜水水位都有所下降，平原堡到高台河段地下水溢出量均有所减少，其中临泽县城、临泽农场一带水位下降很大，年均降幅达0.6m。在酒泉东盆地，潜水水位呈上升趋势，年均升幅0.2m，高台-正义峡河段地下水溢出量有所增加。在张掖盆地的张掖市附近，承压水水位基本保持平稳。在临泽县城、临泽农场一带，承压水水位下降较大，年均降幅达0.9m。在酒泉东盆地酒泉市附近，承压水水位略有上升。其他地段承压水水位略有下降。图9-11 黑河流域中游区潜水系统预测初始等水位线图9-12 黑河流域中游区承压水系统预测初始等水位线图9-13 方案A 黑河流域中游区潜水系统等水位线图9-14 方案A 黑河流域中游区承压水系统等水位线方案B预测结果，如图9-15和图9-16所示。从图9-15和图9-16可见，预测期内各时段也是负均衡，地下水水位呈总体下降趋势。除张掖盆地东南部民乐、六坝、毛城子和安阳一带潜水水位基本保持不变之外，张掖市西北部和平原堡一带的潜水水位略有上升，在附近河段河水溢出量略有增加。在其他地段，潜水水位普遍下降，平原堡到高台河段地下水溢出量均有所减少。其中在临泽县城、临泽农场一带，潜水水位下降很大，年均降幅达1.1 m。在酒泉东盆地，潜水水位持续上升，年均升幅0.2 m，在高台-正义峡河段地下水溢出量有所增加。在张掖盆地，承压水水位普遍下降。其中在张掖市附近，承压水水位年均降幅达0.8 m，在临泽县城、临泽农场一带承压水水位年降幅达1.2 m。在酒泉东盆地，除酒泉市附近承压水水位略有上升之外，其他地段的承压水水位呈下降趋势，年降幅达0.3m。图9-15 方案B 黑河流域中游区潜水系统等水位线图9-16 方案B 黑河流域中游区承压水系统等水位线