如何求解差分方程？

隐式差分格式建立后，就可在每个节点上生成一个差分方程，构建与求解域中内节点数目相同的差分方程，组成一个线性方程组。接下来就是求解方程组，获得所需结果。求解方法一般可以分为直接求解法和迭代求解法。直接求解法通过消元直接求得线性方程组的解，常用的直接求解法有高斯消元法、逆矩阵法、主元素法等。该法需要大量内存，因而工程计算中较少采用；迭代求解法就是通过一系列迭代运算，使每次迭代得到的近似解逐渐逼近于真实解，当解的变化量（有时为残差变化量）小于某个设定的指标时，则称迭代已经收敛，此时得到的结果即为原线性方程组的解。迭代法是数值求解中最常用的方法，许多成熟的计算软件，如地下水流三维有限差分模拟软件Visual Modflow便采用了迭代法求解线性方程组。Visual Modflow中的“Modflow”是Visual Modular Three-Dimensional Finite-difference Ground-water Flow Model的简称，是目前广泛应用的地下水流和溶质运移计算软件。永兴岛淡水透镜体的数值模拟要处理大量的网格单元，加之边界情况远比二维复杂，若用人工编程处理，必然耗费大量的时间和精力。因此，永兴岛淡水透镜体的数值模拟也采用Visual Modflow计算软件。

（1）Visual Modflow软件简介

Visual Modflow软件是加拿大Waterloo水文地质公司在美国地质调查局发布的Modflow软件基础上，运用可视化技术开发研制的一款功能强大的三维地下水流和溶质运移计算软件，于1994年公开发行，是目前国际上最为流行且被各国同行一致认可的标准专业软件，已在全世界几十个国家使用。根据操作界面与版本更新程度，Visual Modflow可分为4.0、4.1以及4.2等多个版本，最新一个版本为Visual Modflow Flex 2014.1。永兴岛淡水透镜体的数值模拟采用Visual Modflow 4.1。

Visual Modflow由一系列的软件包组成，包括：

①MODFLOW（三维地下水流动模拟）；

②SEAWAT 2000（三维变密度地下水流动模拟）；

③MT3D（多组分溶质运移模拟）；

④MODPATH（溶质质点示踪分析）；

⑤Zone Budget（区域水均衡分析）；

⑥Win PEST（自动校正和预测分析）；

⑦MGO（抽水井优化）；

⑧3D-ExPlorer（三维显示和动画）。

其中，MODFLOW和MT3D是Visual Modflow软件中求解数学模型的两个重要解算器。Modflow解算器包括：预处理共轭梯度法程序包（PCG，Preconditioned Conjugate Gradient）、强隐式处理程序包（SIP，Strongly Implicit Procedure）和WHS解算器（WHS）等。MT3D解算器有隐式共轭梯度（GCG）解算器。

PCG也称为“预处理共轭斜量法”，它是一种对大型线性和非线性方程组迭代求解的方法。PCG通过建立预优矩阵，达到简化计算提高收敛速度的目的。根据建立预优矩阵方法的不同及在不同类型计算机上的适应情况，PCG法包含两种不同的修正PCG的算法：MICCG法和POLCG法，这两种方法各有其优缺点，MICCG法较适用于标量型计算机，而POLCG较适用于向量型计算机。

SIP亦称强隐式迭代法，是一种通过迭代求解一组大型联立线性方程组的方法。其原理是通过引入多余的两个未知量，如水头（变量），将原有的五对角系数矩阵改造成七对角系数矩阵，从而实现原系数矩阵的LU分解，再通过逐步迭代来求解有限差分方程组。SIP的优点是迭代的收敛速度对结点数目及待解问题的性质不敏感，只需要较少的内存就可计算出最终结果，但是求解速度慢，不如PCG方法快。

WHS解算器运用双共轭梯度稳定（Bi-CGSTAB）加速度程序。这种解算器通过一个近似解，反复迭代使其接近一组大型偏微分方程的解。WHS解算器在同一个时间步长内使用内外部双重迭代求解，外部迭代用于改变因式分解的水文地质参数（渗透系数、弥散度、贮水率）矩阵，内部迭代用于迭代求解外部迭代中建立的矩阵。

GCG又称“广义共轭梯度法”，是软件中溶质运移通用的迭代法。GCG法有两个迭代循环：一个内循环和一个外循环。在内循环过程中，当达到用户指定的收敛精度或达到用户指定的最大循环次数时，停止运算；当新的外部迭代开始时，利用最近计算的浓度将这些系数更新，再次循环计算。

Visual Modflow软件可进行三维水流模拟、溶质运移模拟和反应运移模拟。三维水流模拟包括饱和常密度、饱和变密度、变饱和度和水汽流模拟。流动类型有稳定流和非稳定流。溶质运移模拟包括非吸附模拟、线性等温吸附、Freundlich等温吸附、Langmuir等温吸附、一阶动力吸附模拟等；反应类型包括无反应项和一阶不可逆衰变。Visual Modflow特别适用于多孔介质中地下水流的三维有限差分数值计算。程序设计包括一个主程序和若干相对独立的子程序包，每个子程序包中有数个模块，每个模块完成数值模拟的一部分。Visual Modflow引入应力期（Stress Period）的概念。应力期是指用户根据模拟目标的不同自己定义的模拟周期。在每个应力期内，所有外部源汇项强度（如补给量、抽水量、蒸发量等）视为不变。整个模拟时间按计算要求可以分成若干个应力期，每个应力期又分为若干个时间段（Time Step），同一应力期，各时间段既可以等步长，也可以按一个规定的几何序列逐渐增长，每个应力期长度、时段数目以及时段增加因子可以根据计算需要自行设定。如由于季节降雨量不同，而将一个季节作为一个应力期，将天的若干倍数作为计算时段。整个模拟过程就是在一系列应力期的若干个时间段内，通过对有限差分方程组的迭代求解，得到每个时间段结束时的节点值，每次模拟包括三大循环：应力期循环、时间段循环以及迭代求解循环，程序计算流程图如图6.10所示。

Visual Modflow的界面设计包括三大相对独立又相互联系的模块，即前处理模块、运行模块和后处理模块。前处理模块允许用户直接在计算机上赋值所有必要的参数，以便构建一个新的三维模型。用户可以直接在计算机上定义和剖分模拟区域，任意增减剖分网格和模拟层数，确定边界几何形态和边界性质，定义抽水井的空间位置和出水层位，以及非稳定抽排水量。参数菜单允许用户直接圈定各个水文地质参数的分区范围并赋值相应参数，同时，上下层所有参数可相互拷贝。用户还可预先定义水位校正观测孔的空间位置和观测层位，并输入其观测数据，以便在后续的模型识别工作中模拟使用。

(www.xing528.com)

图6.10　Visual Modflow程序计算流程图

运行模块允许用户修改程序包中的各类参数与数值，包括初始估计值、运行时间、各种计算方法的控制参数、激活干湿交替软件包和设计输出控制参数等。用户可以单独或共同执行水流模型（MODFLOW）、流线示踪模型（MODPATH）和溶质运移模型（MT3D）进行计算。各部分均设计了模型识别和校正的菜单。模型校正既可用手工进行，也可用WinPEST自动进行。

后处理模块可自动地阅读每次模拟结果，输出等值线图、流速矢量图、水流路径图，并可借助软件组合的绘图功能模块Visual Groundwater进行三维显示和输出（如三维等值面图的输出等）。后处理模块允许用户以三种方式展示模拟结果。第一种方式，是在屏幕上直接彩色立体显示所有模拟结果；第二种方式，是在打印机上输出模拟评价的成果表格和成果图件；第三种方式，是以图形或文本文件格式输出模拟结果。输出和显示的图形包括可以标记显示水头、降深、浓度、含水层顶底板标高、含水层厚度、渗流速度矢量等的平面、剖面等值线图和平面、剖面示踪流线图，以及局部区域水均衡图等。

Visual Modflow软件的主要特点是：

①具有合理新颖的菜单结构、友好的界面和功能强大的图形可视化特征，用户可方便地建立新模型和修改已有的模型。

②能以平面和剖面两种方式彩色立体显示模型的剖分网格、输入参数、输出结果。

③将数值模拟过程中各个步骤有机地结合在一起，从开始建模、各类水文地质参数的输入与修改，到运行模型、参数的识别，再到输出结果，整个过程非常系统化、规范化。

④Visual Modflow采用可视化数据处理手段并支持TXT、DAT、EXCEL、MAPINFO及CAD等数据格式，能够克服各种数值计算产生的许多弊端，确保数据的安全性、通用性和标准化。

（2）SEAWAT软件包的数学模型

永兴岛淡水透镜体的数学模型为变密度三维地下水流模型，所以，采用Visual MODFLOW中美国地质调查局的SEAWAT 2000程序。SEAWAT合并了MODFLOW水流模型和MT3DMS溶质运移模型，形成了一个能解决耦合水流和溶质运移方程组的程序。

SEAWAT程序中采用的三维变密度水流主管方程为：

式（6.92）中，α，β，γ为直角坐标系，分别对应x，y，z，与渗透系数的主轴同向；hf为淡水当量水头，即测压管水头；Kf为渗透系数（以淡水为参考）；Sf为贮水率（以淡水为参考）；θ为有效孔隙度；ρ是混合流体的密度；ρf是淡水的密度；qs是单位体积多孔介质源（或汇）的流量；为源（或汇）的密度。

比较式（6.92）与式（6.50）：将式（6.92）中的α，β，γ换为x，y，z，并将式（6.92）中的Kf、Sf转化为以混合流体ρ参考条件下的K、Ss，则式（6.92）与式（6.50）相同。

SEAWAT程序中采用的溶质运移主管方程：

式中，D为水动力弥散系数张量，v为地下水平均流速，Cs为源（或汇）流体的浓度，Rk为化学反应项。

将式（6.93）与式（6.59）进行比较：若式（6.93）不考虑化学反应项，则式（6.93）与式（6.59）相同。

经过以上分析，前面导出的三维变密度地下水流和溶质运移的主管方程与SEAWAT程序包采用的主管方程形式一致，可以采用Visual Modflow软件对淡水透镜体进行数值模拟。