数学模型的求解与结果分析

淡水透镜体的数学模型，由于主管方程的非线性和边界的复杂性，通常不能求得解析解，工程上更多地依赖于数值解法。常用的数值解法有有限差分法、有限元法和有限体积法等。以下介绍有限差分法。

有限差分法的基本思想是将求解域格网化，按网格的分布用折线取代不规则的边界线，以差商代替微商，将偏微分方程离散化，从而把求解域上连续的水头分布值离散成全部网格上的有限个水头值，得到离散的只有有限个未知数的差分方程及对应的边界值，并将差分方程的解作为微分方程的数值形式的近似解。这样，用有限差分法求解淡水透镜体数学模型时，就需要将求解域网格化，选择时间、空间步长，构造差分方程，确定解法。

（1）求解域的格网化

根据永兴岛的平面图，建立直角坐标系xOy，Ox轴由西向东，Oy轴由北向南，然后用二簇分别平行于Ox和Oy的线段将永兴岛划分成许多网格，x和y方向上的网格间距Δx和Δy可以不同，但通常为简化计算，均取等间距网格即正方形网格如图5.3所示。

图5.3　计算用网格与编号

差分网格划分后，需要对每个网格进行编号，每个网格的编号以其在x，y方向上网格数来表示，记为（i，j）。其中i为x方向的编号数，自西向东排列，j为y方向编号数，自北向南排列（图5.3）。

（2）时间、空间步长的选取

淡水透镜体的主管方程中包含有对时间的偏微分项。因此，在对空间坐标进行离散的同时，也要对时间进行离散。然而对时间和空间步长的选取，除在显式差分格式中由稳定性的要求可以确定时间和空间步长的相对大小外，目前尚无严格选取时间和空间步长的理论和方法，通常靠试算和经验加以确定。本应用中，时间步长取Δt＝1 d，空间步长取Δx＝Δy＝100 m。

（3）差分方程的构造

由主管方程（5.12）于是主管方程（5.12）转换为：

对方程（5.17）进行离散。为避免离散方程中多孔介质孔隙率n与表示时间层的n混淆，将方程（5.18）中n用μ替换。方程的离散有两种格式：一种是显式，另一种是隐式。

显式：

隐式：

式中，下标i，j表示第i行第j列网格之值，上标n和n＋1分别表示第n和第n＋1时间层之值。

比较显格式和隐格式可以发现，用显格式进行计算时，任一时间层各网格上的值可由前一时间层各网格上的值明显地表示出来。因此，由差分方程逐层求解网格值时，只需将已知网格之值（如初值和边值）代入差分方程，就可逐一地算出各时间层上各网格之值，而不必求解代数方程组，这是一种代数四则运算，求解过程就变得异常简单。但显格式是条件稳定的，时间和空间步长都会受到严格限制，如果取值不当，就有发散的危险。而对于隐式，任一时间层各网格上之值不仅与前一时间层上的网格值有关，而且也与同一时间层上相邻网格之值有关，这样各时间层上的网格值就不能像显式那样由前一时间层上之值明显表示。因而，用隐式差分方程求解，不能仅作代数四则运算来获得结果，必须逐层地求解代数方程组。这种求解方法比较麻烦，但隐格式是无条件稳定的，网格步长不像显式那样会受到严格限制，也有利于减少计算量和提高计算精度。

为了提高差分格式的精度，可以构造一种既有显式又有隐式形式的差分格式。即所谓“加权平均法”，就是差分方程中的空间差商既不全在n时间层上取值，也不全在n＋1时间层上取值，而是在这两个时间层上取加权平均。首先，用正方形网格将方程（5.18）和（5.19）的求解域网格化。

取Δx＝Δy＝a，并令

方程（5.18）为：

则式（5.21）为：

同理，由方程（5.19）有

式中，

设ω为加权系数0≤ω≤1，由ω×式（5.24）＋（1-ω）×式（5.23）有：

整理得：(www.xing528.com)

将式（5.20）代入式（5.27），得：

由式（5.28）解出：

方程（5.29）即为加权平均的差分格式。当ω＝时，该差分格式是无条件稳定的，且有二阶精度。这个格式是Crank-Nicolson提出的，又称“Crank-Nicolson格式”。

（4）差分方程的求解

方程（5.29）是一个线性方程组，有多种解法。本应用中采用Gauss-Seidel迭代法求解。Gauss-Seidel迭代法比简单迭代法求解时收敛更快。

方程（5.29）中的εij是（i，j）微元单位面积上的补给量，主要是降雨，植被蒸腾和渗漏损失则以负补给量包含其中。根据永兴岛的降雨，综合资料值，取垂直补给量为年均降雨量的40%。由该方程可求得Hij。

每一网格上淡水透镜体表面高出海平面的高程hf，i，j，由下式求得：

淡水盐水界面在海平面以下深度hs，i，j，则由下式求得：

淡水透镜体的贮水量V则近似为：

（5）计算结果

用上述淡水透镜体数学模型，通过数值解法对永兴岛淡水透镜体进行模拟，可得到永兴岛上淡水水头分布和淡水贮量。表5.4和表5.5分别是第7行和第9行各列节点处淡水透镜体的厚度h。图5.4—图5.7分别是第7行、第9行、第5列、第9列各节点处海平面以下淡水透镜体的模拟界面与部分节点处透镜体厚度h的值实测值。图5.8淡水透镜体的等深图（以海平面为基准面）。

表5.4　第7行海平面以下淡水透镜体的厚度

表5.5　第9行海平面以下淡水透镜体的厚度

图5.4　第7行模拟界面与部分节点处厚度的值实测值

图5.5　第9行模拟界面与部分节点处厚度的实测值

图5.6　第5列模拟界面与部分节点处厚度的实测值

图5.7　第9列模拟界面与部分节点处厚度的实测值

上述各图中的实测值是指：对部分西沙永兴岛水井中水面距地表距离，以及地表的海拔高程的测量值，求得的淡水水面在海平面以上的高程hf推得的淡水界面在海平面以下的深度hs。各井hs的计算值与实测值见表5.6。计算值与实测值比较接近，大部分测点的误差为3%～10%。图5.8为淡水透镜体的等深线，其外形与岛的平面外形十分相似。这表明淡水透镜体的二维数学模型能较好地反映透镜体的动力学特性。

图5.8　淡水透镜体等深图（单位：m）

表5.6　hs的模拟值与实测值