SciFEA软件模拟非饱和土壤渗流

SciFEA 饱和—非饱和土壤水分运动计算

计算机是现代科学技术的重要支柱，它不仅能带来巨大的经济效益，而且能带来深远的社会效益，对科学技术的发展起着巨大的推动作用。其在农业领域求解难于得到严格解析解的非线性偏微分方程的广泛应用成为现代化中计算问题的重要标志。国内外计算机已在农业科学实验、农业政策制定、资源普查与监测、植物和环境保护、土地规划及合理利用、农业信息处理与预测、自动化温室控制、农业气候与作物种植、土壤墒情监测与农田灌溉、农业机械化检测以及农药化肥生产控制等领域广迂使用。以节水灌溉为目标的现代灌溉管理需要科学地进行用水管理和定量描述并预测灌溉—蒸发条件下的土壤水分在土壤耕作层内的运动规律。土壤水分运动同时也是制约土壤溶质迁移的主要因素，对于干旱、半干旱的次生盐碱化威胁较严重的地区，这是至关重要的。因此，对于土壤水分的模拟研究不仅有助于现代节水灌溉管理，而且是土壤溶质迁移规律研究的基础。 随着科学的进步与发展，土壤水分运动的研究已由定性描述的形态学观点逐步发展成为定量研究阶段。土壤水分运动基本方程是一个非线性偏微分方程，在许多定解条件下，很难求得解析解。近年来．国内外许多学者都对求解非饱和土壤水分运动方程的数值方法进行了研究。

1问题描述

土壤水分运动基本方程

以土壤水势为因变量的一维饱和—非饱和等温、均质、各向同性土壤水分运动基本方程(Richards，1931)如下：

()()()C K Z t z z ∂Ψ∂∂⎡⎤Ψ=ΨΨ−⎢⎥∂∂∂⎣⎦ (1) 其中：Ψ是土壤水势(cm)；()C Ψ是容水度；()K Ψ是非饱和导水率(cm/min)；

Z是垂向坐标(cm)，向下为正；t 是时间(min)。并假设Ψ和K 均为土壤含水率θ的

单值函数。

初始条件

0(,0)(),0Z Z Z H Ψ=Ψ≤≤ (2) 其中,H 是下边界深度，H →∞。

边界条件

假设地表处于湿润状态下的入渗，且地表处基质势维持不变时；或在蒸发条件下，地表处于风干状态时，均可将地表处土壤水势作为已知水势处理。

01(,)(,),Z t Z t Z Ψ=Ψ∈Γ (3)

计算参数

考虑从地表到地下120cm处的土壤，初始条件为一分段分布的土壤水浓度函数，具体表达式如下：

0.15/1200

(,0)0.2z z θ+⎧=⎨⎩ (4)

在地表进行灌溉时，地表含水率保持为0.38, 地表一下120cm处仍然保持0.2

的含水率。导水率535.8() 3.2410K e θθ−=×(单位cm/min)；扩散率特征曲线为如下

的分布形式：

125.30.84.4210 0.35()1885 0.350.38e D e θθθθθ−−⎧×≤=⎨<≤⎩ (单位cm2/min) (5)

计算t=120min和t=540min时土壤水浓度。

2求解步骤

2.1模块挂载

(1) 启动SciFEA ，选择“工具”—>“模块挂载”菜单弹出如图1-1所

示的对话框，在对话框当中选择“非饱和土壤渗流模拟”。

图1-1 选择项目类型对话框

(2)点击“SciFEA模块挂载”栏中的“非饱和土壤渗流模拟”选项。如图1-1所示。

(3)点击“确定”按钮完成项目类型的选择。

(4)单击工具条中的按钮弹出如图1-2所示，并且选择打开udm.pro文件。

图1-2

(5)点击“打开”按钮之后，弹出如图1-3：

图1-3

(6)最后点击“OK”按钮。

2.2设置材料参数和边界条件

(1)单击工具条中的按钮弹出如图1-4和图1-5所示材料参数数据输入表

格。

图1-4 渗流参数输入

图1-5初始参数

（2）单击工具条中的按钮弹数据输入表格。按照问题描述中的参数依次填入边界条件数据表格。填写完成后如图1-6所示。

图1-6 填写完成渗流边界

3建模、设置材料属性和施加边界条件

3.1启动GID以创建模型。单击工具条中的按钮弹出如图1-7所示前处理初始化窗口。

图1-7前处理初始化窗口

3.2建模

3.2.1点击菜单【Geometry】-【Create】-【Line】，在命令栏依次输入点坐标：

图1-8建成的模型图

3.2.2定义材料特性、施加载荷约束。点击菜单【Data】-【Conditions】弹出“Conditions”对话框。

①定义材料。点击表示面单元的按钮，选中下拉菜单中的“line－warka，最后得到如图1-9所示：

图1-9

②施加载荷。点击线单元按钮，选中下拉菜单中的“surface-aeq4”，最后得到如图1-10所示。

图1-10 施加载荷

需要说明的是，这几步定义材料，施加载荷的操作没有先后顺序要求，但必须保证不能遗漏某个条件的定义。

3.3划分网格和导出数据

(1).划分网格。点击菜单【Meshing】-【Generate】，要求定义单元尺寸大小，点击“OK”按钮，此时弹出“Progress in meshing”进度条，网格划分完毕后弹出消息框显示总的单元数和节点数，如图1-11所示。

图1-11 划分网格得到的单元节点信息

点击“OK”按钮后即可得到划分完的网格模型，如图1-12

图1-11 网格图

(2)把前处理数据转化为SciFEA所需计算格式。点击【Files】-【Save】菜单，

保存为GID文件。点击菜单【Calculate】-【Calculate】，弹出“process info”

消息框，如图1-12所示，点击“ok”按钮，然后退出Gid，至此前处理工作结束。

图1-12 转化数据消息框

4有限元计算

(1) 单击工具条中的按钮弹出如图1-12所示计算模拟窗口。

图1-12 计算模拟窗口

5结果分析

对计算结果进行分析属于后处理，GID提供了丰富的后处理操作，可以从不同角度，不同方式来输出计算结果。