Hydrus-2D中文程序说明

HYDRUS_2D中文程序说明李道西HYDRUS_2D是一个可用来模拟地下滴灌土壤水流及溶质二维运动的有限元计算机模型。该模型的水流状态为二维或轴对称三维等温饱和－非饱和达西水流，忽略空气对土壤水流运动的影响，水流控制方程采用修改过的Richards方程，即嵌入汇源项以考虑作物根系吸水。程序可以灵活处理各类水流边界，包括定水头和变水头边界、给定流量边界、渗水边界、自由排水边界、大气边界以及排水沟等。水流区域本身可以是不规则水流边界，甚至还可以由各向异性的非均质土壤组成。通过对水流区域进行不规则三角形网格剖分，控制方程采用伽辽金线状有限元法进行求解。无论饱和或非饱和条件，对时间的离散均采用隐式差分。采用迭代法将离散化后的非线性控制方程组线性化。对于非饱和土壤水力特性，HYDRUS_2D采用VG模型进行描述，嵌入了Scott[1983]、Kool和Parker[1987]经验模型中的假定：吸湿（脱湿）扫描线与主吸湿（脱湿）曲线成比例变化。并运用一个比例程序，将用户定义的水力传导曲线与参考土壤相比较，通过线性比例变换，在给定的土壤剖面近似水力传导变量。第一节程序的基本模块HYDRUS_2D程序模块可以顺序嵌套调用，由以下七个基本模块组成：HYDRUS2D：主程序，定义系统的整个计算机环境。它控制整个程序的运行过程，根据需要调用相应的子程序模块。程序执行前，首先需选定模拟选项，包括水流、溶质运移、热运移或是否考虑根系吸水等；然后给定时空单位、土壤水力参数以及用来模拟的边界条件。程序执行后，可输出一系列土壤水力特性曲线、设定观测点处随时间变化的含水率或负压水头曲线，以及沿边界的实际或累积水通量。输出文件还可提供质量平衡信息和逆向最优结果。.ProjectManager：该模块用来管理已建立的工程数据，包括打开、删除、重命名工程和保存工程的输入输出数据等。每个工程可能是针对不同的具体问题，ProjectManager会自动将每个工程单独建立一个以工程名命名的文件夹保存相应的工程数据。GEOMETRY：该模块是一个可用鼠标或键盘图绘水流区域并输出的CAD程序，也可通过导入二进制文件的方式实现。水流边界可以由直线、圆、弧或多义线等不同曲线组成；内部边界也可由内部曲线给定，如排水沟、井等。另外，还可以对已绘区域进行修改，如删除、复制、移动、旋转等操作。MESHGEN2D：该模块用来将二维的水流区域离散成不规则的三角形网格。第一步：边界离散化，边界结点数和其密度可由用户自行确定。第二步：整个水流区域的基于Delaunay规则的三角形离散化。按照默认的光滑因子，可以将指定的水流区域自动生成最优的三角形有限元网格，例如对于指定的边界结点，它可以生成最小的三角形单元剖分。BOUNDARY：该模块用来让用户给定特定情况的初始和边界条件，以及取定观测点等。HYDRUS2：该模块是一个可用来模拟二维非饱和土壤水运动的FORTRAN程序。模型可求解含根系吸水源汇项的Richards方程，可以灵活处理各类水流边界，包括定水头和变水头边界、给定流量边界、渗水边界、自由排水边界、大气边界以及排水沟等。针对离散化控制方程后的系数矩阵的不同形式，采用了不同的求解方法，例如带状矩阵对应高斯消去法；对称矩阵对应共轭梯度法；非对称矩阵，对应ORTHOMIN法。另外，该程序升级版本还包含了一个参数最优算法，可对各种土壤水力参数从几个观测的数据出发进行逆向估计。对于土壤含水率或负压水头数据，采用了Marquardt_Levenberg非线性最优化技术估算土壤水分特征曲线中的经验参数；对于持水或导水率数据，则将待优参数通过罚函数约束在某个可行区域（贝叶斯估计），然后寻求最优。GRAPHICS：该模块用来将输出结果表示成图形。包括等值线图、光谱图、流速矢量图以及等值线图和光谱图的随时间变化的动画显示等。第二节程序菜单命令说明图5－1是HYDRUS_2D程序的主界面。各菜单的子菜单使用及说明如下：图5－1HYDRUS_2D程序的主菜单2.1Project主菜单ProjectManager：用来建立新的工程，或者管理已建立的工程数据，包括打开、删除、复制、重命名工程等，如图5－2。图5－2ProjectManager对话框SaveData：用来保存工程的输入输出数据。Exit：退出主程序。2.2MainInf.主菜单MainProcesses：通过它，用户可以对工程进行描述，如“InfiltrationfromaSubsurfaceSource”，并选择需模拟的对象（可复选）：水流、溶质运移或根系吸水。当勾选“水流”时，程序会自动考虑瞬态流，否则当稳态流处理，如图5－3。图5－3MainProcesses对话框MainGeometryInformation：包括四个部分，如图5－4：1、单选长度单位（LengthUnits）：毫米、厘米、米。一旦输入数据的长度单位发生改变，程序会自动将所有变量换成对应的新的单位。2、单选水流类型（TypeofFlow）：二维水平流、二维垂直流、三维轴对称流。当为三维轴对称流时，程序自动视z轴为垂向对称轴。3、单选几何类别（GeometryofType）：矩形网格、不规则网格。仅仅是当主程序缺乏外部嵌入式模块GEOMETRY和MESHGEN时才采用简单的矩形网格剖分。4、填写数字框：土壤质地种数（NumberofMaterials）、子区数（NumberofSubregions）。当水流区域有夹层时，土壤质地种数将大于默认值1；子区数决定了用户希望进行质量平衡计算的区域数。图5－4MainGeometryInformation对话框图5－5MainTimeInformation对话框MainTimeInformation：包括三个部分，如图5－5：1、单选时间单位（TimeUnits）：秒、分、小时、天。该项一旦发生改变，程序中的所有相关变量的单位均会自动进行相应变化。2、时间离散选项（TimeDiscretization）：五个数字框分别表示模拟初始时刻、模拟总时间、初始时间步长、允许最小和最大时间步长。需要说明的是，HYDRUS2为了提高计算效率，视求解情况会自动调整时间步长（其调整规则将在菜单Parameters/WaterFlow…/IterationCriteria中进行介绍），所以存在这些选项的给定。3、边界条件选项（BoundaryConditions）：缺省为不勾选“随时间变化的边界条件”，则所有边界条件均视为恒定。否则存在非稳定边界，需填写非稳定边界数目的数字框，如腾发、降雨、灌溉等均可能导致非稳定边界的存在。PrintInformation：如图5－6。通过它，用户可选择决定是否需要输出每个时阶的平均压力水头、平均或累积含水率及溶质通量、迭代信息等（T-LevelInformation）；是否需要输出结点流量等（FluxInformation）；是否需要屏幕输出（ScreenOutput）。若给定输出次数，点击（SelectPrintTimes），则要求用户进一步给定需要输出更详细信息的各个时刻，所有时刻不可大于设定的模拟总时间。5－6PrintInformation对话框2.3Geometry主菜单EditTransportDomainGeometry：用来调用GEOMETRY模块，如图5－7，以便图绘水流区域。其绘制过程与一般的CAD绘图软件差不多。水流边界可以由直线、圆、弧或多义线等不同曲线组成；内部边界也可由内部曲线给定，如排水沟、井等。注意主界面右下角的状态栏，可以粗略估计鼠标所在的坐标位置；通过设置菜单选项Options/GridSettings…还可以精确定位点的坐标。GenerateFiniteElementMesh：用来调用MESHGEN2D模块，如图5－8，以便将已绘制的二维水流区域离散成不规则的三角形网格。第一步：通过菜单BoundaryPoints/NumberofPointsonaCurve…设定边界结点总数，默认值为120。该项决定了离散密度的大小，边界结点总数越多，三角形网格密度越大，计算越精确，但计算机的运算时间会成倍增加。图5－7GEOMETRY模块主界面图5－8MESHGEN2D模块主界面第二步：通过菜单BoundaryPoints/InsertaFixedPoint给定边界离散固定点，以设定局部网格密度。例如，对于在地下滴灌的滴头埋深处至少应插入一个离散固定点，以使该处网格剖分相对加密提高该处计算精度。注：边界角点均默认是离散固定点。第三步：通过菜单BoundaryPoints/MeshDensityataFixedPoint…设定各离散固定点处局部网格密度。离散固定点两侧（按逆时针方向分为左右两侧）均需设置相应的密度值，密度值范围大于0.01小于100。密度值的含义实际上是一个比较量。例如有两个相邻离散固定点A、B，它们的左右密度值分别设定为A：3、2和B：0.2、0.2。若离散后A点左右侧的步长分别为、，B点左右侧的步长分别为、，则1a2a1b2b1231.52aa==，120.210.2bb==，113150.2ab==。而A、B间的网格密度分布按结点数从3向0.2递减。第四步：通过菜单MeshGeneration/AutomaticGeneration或点击左下角的“MakeMesh”按钮，即可自动将整个水流区域剖分成不规则三角形网格。当然，也可以先通过MeshGeneration/EditParameters设定平滑参数，但MESHGEN2D模块按默认值就已经可以使网格剖分达到最佳的平滑效果，所以，没有特殊情况一般不修改默认值，而跳过此步直接生成。DomainandBoundaryParameters：用来调用BOUNDARY模块，以便让用户给定特定情况的初始和边界条件，以及取定观测点。第一步：通过菜单Condition/InitialConditions/WaterFlowIC给定水流初始条件，即整个水流区域的负压水头或含水率初始分布。点击靠近左下角“Select”按钮选定区域，然后点击左下角“SetValue”按钮即可设定初始值，如图5－9。调出的对话框中的单选按钮还可使水头分布随深度线性变化。图5－9BOUNDARY模块初始条件的设定图5－10BOUNDARY模块边界条件的设定第二步：通过菜单Condition/WaterFlowBC给定水流边界条件。点击左侧各类边图5－10BOUNDARY模块中的边界条件的设定界按钮，然后选定相应的结点即可。对于定压力水头边界和定流量边界还可能弹出输入数值对话框，如图5－10。第三步：通过菜单Condition/ObsevationNodes设定观测点，观测点必须是网格结点，可以通过View/Zoom进行放大精确点选，如图5－11。观测点应尽量有代表性，以便进行规律分析。图5－11BOUNDARY模块中的观测点的设定2.4Parameters主菜单WaterFlow…：包含两个子菜单。1、IterationCriteria：包括两个部分：迭代规则（IterationCriteria）和时间步长调整（TimeStepControl），如图5－12。迭代规则数值框依次表示：最大迭代次数、含水率允许误差、负压水头允许误差。关于时间步长调整，可作如下说明：为了提高求解效率，HYDRUS_2D模拟计算过程中，根据指定的初始时间步长（如1E－6h），视求解情况会自动调整时间步长，其规则如下：（1）时间步长不能小于指定的最小时间步长（如1E－6h），也不能大于指定的最大时间步长（如1h）。（2）如果达到收敛时的迭代次数小于等于某指定下限（默认为3），在下一时段内，将时间步长乘以一个大于1的常数（默认为1.3），以加大时间步长；如果达到收敛时的迭代次数大于等于某指定上限（默认为7），在下一时段内，将时间步长乘以一个小于1的常数（默认为0.7），以减小时间步长。（3）如果在某个时间步长内，迭代次数大于指定的最大迭代次数（实际模拟过程中可根据所要求的