MODFLOW的特点

5.1.2.1 MODFLOW 的最显著特点是采用了模块化的结构

这里所说的模块化结构是指由一个主程序和一系列具有高度独立性的子程序——即模块所组成的计算机程序结构。它一方面将许多具有类似功能的子程序组合成为子程序包，另一方面用户可以按实际工作需要选用其中某些子程序包对地下水运动进行数值模拟。此外，这种模块化结构使程序易于理解、修改，甚至添加新的子程序包。自从MODFLOW问世以来，已经有许多新的子程序包被开发出来，用来解决MODFLOW本身不能解决的问题。例如，1989年Prudic开发出了模拟河流与含水层之间水力联系的河流子程序包；Leake和Prudic 1988年、Hoffmann和Leake 2003年开发出了模拟由于抽水引起地面沉降的子程序包；Hsieh 和Freckleton 1993年开发出了模拟水平流动障碍（Horizontal flow⁃barri⁃er）的子程序包；Halford和Hanson 2002年开发出了模拟限制降深抽水条件下的混合井流子程序包等。这些新子程序包的加入，大大提高了 MODFLOW 的应用范围。虽然MODFLOW本身伏段仅限于模拟地下水在孔隙介质中的流动，但大量实际工作表明，只要恰当使用，MODFLOW也可以用来解决许多地下水在裂隙介质中的流动问题。不仅如此，经过合理的线性化，Guo 1995年把 MODFLOW 用来解决空气在土壤中的流动问题。Guo 和Bennett 1997年将MODFLOW与其他用于溶质运移模拟的程序结合起来，用来模拟诸如海水入侵等地下水密度为变量的问题。MODFLOW之所以能得到如此广泛的应用，一个很重要的因素就在于其设计上的合理性裂罩。参与MODFLOW设计、编程的工作人员大多都是具有丰富实际经验的著名水文地质学家。因此，MODFLOW的使用说明书本身就是一部关于地下水数值模拟方面的教科书，很值得广大水文地质工作者阅读参考［15，16］。

5.1.2.2 离散方法简单

包括空间上的离散、时间上的离散。

在空间的离散上，MODFLOW对含水层采用三个轴向上不等距正交的长方体剖分网格，这种网格的优点在于用户易于准备数据文件，便于文件的规范化输入。

在时间的离散上，MODFLOW模拟系统引入了应力期的概念，可以把整个模拟期划分为外应力（如抽水量、蒸发量、补给量等）保持不变的若干个应力期，每个应力期可再细划分为若干个时间段。比如以季度长度为应力期，以天的倍数为时间段。在同一个应力期中，各时间段既可以等步长，也可以几何序列逐渐增长；通过对有限差分方程组的求解，可得到每个时间段末刻的水头值［9，17］。

5.1.2.3 求解方法多样化

求解的方法可以分为直接求解方法和迭代求解方法。MODFLOW原来含有两种迭代求解子程序包：SIP方法（或称为强隐式法）、SOR方法（或称为逐次超松弛迭代法）。由于MODFLOW的模块化结构，Mary Hill于1990年设计增加了一种新的迭代子程序包肆厅闹，即PCG子程序包，该子程序包采用PCG方法（或称为预调共轭梯度法）迭代求解。

对于MODFLOW的多个求解子程序包，一方面，用户可以根据问题的实际情况选用比较合适的求解方法；另一方面，对于某一特定的实际问题，由于水文地质条件的复杂性，用户选择不同的求解子程序包可能都会收敛，也可能只收敛于一种（或几种）求解方法而不收敛于另一种（或几种）求解方法。通过我国国内多项实践应用，PCG法和SIP法实用可靠，而运用SOR子程序包求解的结果精度低，不宜采用［9，16］。

粘塑性自洽多晶体塑性模型（Visco-plasitic Self Consistant，VPSC）是由美国Los Alamos国家实验室的C.N. TOME 教授和R.A. Lebensohn教授联合开发的，用于模拟多晶体材料塑性变形过程宏观力学响应和微观结构演化的本构模型。这里，多晶体材料可以是七大晶系中的任一一种或多种复合，比如立方晶系中的铝合金、铜、铁、钢，密排六方晶系中的镁合金、钛合金、锆合金，等等；而塑性变形过程可以是比较简单的变形过程如单向拉伸、单向压缩、剪切，或更加复杂的自定义变形过程。

金属材料的宏观本构模型一般包括3个部分：d性法则，屈服函数和硬化方程，其中屈服函数和硬化方程最为重要。屈服函数又包括各向同性屈服如Mises屈服或各向异性屈服如Hill48屈搜饥服，硬化方程主要包括描述屈服面大小变化的等向强化和描述屈服面位置变化的随动强化，各自都有许多具体的形式。屈服函数和硬化方程的不断发展，形式上越来越复杂，目的都是为了使得本构模型更加“准确”，具体表现就是数值计算结果与实验更接近。宏观本构模型往往嵌入到商业有限元软件如ABAQUS，ANSYS中，模拟搭漏灶实际零件的加工过程，从而对零件设计和工艺参数提供指导意义。宏观本构模型一般仅要求从数学上建立应力应变响应的关系表达式，而不考虑金属材料变形的物理本质，这使得其缺乏科学意义，且对于具有复杂变形机制的材料如镁合金等，宏观本构难以准确描述其变形过程。与宏观本构模型不同，晶体塑性模型从单个晶粒出发，首先建立单个晶粒的本构方程，再通过一定的方式将单晶相互耦合，从而建立多晶体在外加载荷下的力学响应，晶体塑性模型不仅能够描述宏观本构模型中的应力应变曲线和屈服面演化等宏观信息，还可以给出如滑移模式相对开动率和织构演化等微观信息，如这里所提到的VPSC模型。VPSC通过Eshelby夹杂理论建立相互作用方程，从而将单晶相互耦合起来，被成功用于模拟许多金属材料的塑性变形过程。下面通过1个算例来详细说明VPSC的功能。

算例： AZ31镁合金变形过程的VPSC分析

镁合金具有密度低和机械加工性能好的有点，被广泛用于航空航天和汽车交通等领域。从晶体结构，其属于密排六方（hexagonal close-packed structure，HCP）晶体结构，对称性较低；从微观塑性变形机制看，其滑移模式种类多，如图1所示，且变形过程会发生孪晶，如图2所示；从晶粒取向看，由于其初始织构一般存在明显的择优取向，因此力学性能具有很强的各向异性。镁合金的力学行为十分复杂，现象学模型难以对其准确描述，需要从介观尺度考虑每个晶粒的变形及其取向演化，才能准确描述其变形过程，下面介绍VPSC在模拟AZ31镁合金单向拉伸和单向压缩过程中的应用。

在进行模拟计算之前，我们需要进行输入文件的准备，包括一个主输入文件及三个子输入文件，分别是单晶文件、织构文件和加载文件。其中，主输入文件给出了所含相信息及计算参数等方面的设置，单晶文件给出了AZ31镁合金的d性刚度阵、滑移系和孪晶系以及硬化方程所需的参数；织构文件给出了AZ31镁合金的晶粒数量和织构信息，织构信息一般通过实验获得，图3给出了AZ31镁合金的知扮初始织构，可以看出材料初始有很强的{0001}基面织构；加载文件给出了加载过程，例如这里是单向拉伸和单向压缩，不仅要指定加载方向，还要给定加载速率。

图1. 镁合金的滑移变形机制：(a) 基面（BASAL）<a>滑移；(b) 柱面（PRISMATIC）<a>滑移；(c) 锥面（PYRIAMIDAL）<a>滑移；(d) 二阶锥面（PYRAMIDAL-II）<c+a>滑移

图2. 镁合金的孪晶变形机制：(a) 孪晶示意图；(b) 拉伸孪晶；(c) 压缩孪晶

图3. AZ31轧板的初始织构

在输入文件准备好之后，我们就可以通过VPSC程序分别计算AZ31镁合金在单向拉伸和单向压缩的变形过程了，在经过计算之后，我们可以对比两种加载过程的应力应变曲线，如图4所示。可以看出，AZ31镁合金轧板在沿着RD方向进行拉伸和压缩时，应力应变曲线有明显差异：（1）在单向拉伸下，屈服应力约为128 MPa，而在单向压缩下，屈服应力约为107 MPa；（2）在单向拉伸下，流动应力随着变形的增加呈抛物线式增长，而在单向压缩下，流动应力曲线呈S型，这是由于两种变形方式下，微观塑性变形机制的开动情况差异性导致的，如图5所示。

图4. AZ31镁合金轧板沿着RD方向单向拉伸和单向压缩应力应变曲线

(a) (b)

图5. (a) 单向拉伸过程滑移模式开动情况；(b) 单向拉伸过程滑移模式开动情况

由图5可以看出，在单向拉伸和单向压缩变形过程中，AZ31镁合金内不同滑移模式的开动情况存在很大差异：在单向拉伸变形时，基面和柱面滑移模式占主导地位，且随着变形的增加，基面滑移的开动率逐渐减小，柱面滑移开动率逐渐增加；而在单向压缩变形时，基面滑移始终是最主要的塑性变形机制，在小变形时，孪晶也对塑性变形起20%的贡献，当变形量超过10%以后，孪晶开动率只占不到10%，而柱面和二阶锥面滑移模式逐渐开动。正是由于微观塑性变形机制在不同加载过程中的开动情况差异，使得AZ31镁合金在不同加载过程中的应力应变曲线出现巨大差异，而VPSC模型能够很好地模拟并解释这一过程。通过分析AZ31轧制镁合金的初始织构（图3），可知在轧板内部大多数单晶的取向均为<0001>//ND，具有较强的基面织构，如图6所示。那么在两种加载方式下的变形，内部织构会发生如何变化呢？

图6. AZ31轧制镁合金的基面织构示意图

图7分别给出了单向拉伸和单向压缩过程发生20%变形后的织构情况，对比初始织构可知：在单向拉伸过程中，<0001>基面织构沿着RD方向向ND方向收缩，且由于柱面<a>滑移模式的开启，(100)极图出现了6等分强度峰值分布；在单向压缩过程中，由于拉伸孪晶的开启，基面织构完全消失，六方棱柱的c轴逐渐往RD方向偏转，表现出典型的孪晶织构特征。

(a)

(b)

(a)

(b)

图7. AZ31镁合金变形织构：(a) 单向拉伸20%；(b) 单向压缩20%

由此可见，VPSC不仅能够模拟AZ31镁合金在不同加载过程中的宏观应力应变曲线，还能分析出变形过程中的微观变形机制演化和织构演化，这对于深刻理解变形的微观机理具有重要意义。除了这里给出的镁合金变形过程的模拟，VPSC还能够模拟多种晶体结构（如FCC金属、BCC金属等）、多种相组成（如BCC+HCP、FCC+BCC）在多种加载方式（如剪切、双拉、自定义应变路径等）、多种应变率（准静态~高应变率）、多种温度条件下的变形过程，能够有效分析材料变形过程的微观机理。

---