数值方程与数值模拟

Question

如题所述

Answer 1

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