蓄水单抗非均质土土壤水分运动的数值模拟

第３３卷第２期　２００２年３月　太原理工太学学报　Ｖｏ【３３　ＮＯ　２　ＪＯＵＲＮＡＬ　ＯＦ　ＴＡＩＹＵＡＮ　ＵＮＩＶＥＲＳＩＴＹ　ＯＦ　ＴＥＣＨＮＯＬＯＧＹ　Ｍａｒ．２００２　文１１ｔｆ￣号：１００７　９４３２（２００２）ＯＺ一０２１４－０５　蓄水单坑非均质土土壤水分运动的数值模拟　郭文华，孙西欢　潘军锋　仇亚琴　（太原理工大学建筑与环境工程学院）　（山西省水利厅）　（天津城建学院）　摘　喜：针砷蓄水坑灌法中非均质单坑土壤水分问题，建立了数学模型，井采用ＡＤＩ（变替方　向臆式）差分击进行求解以模拟非均质土土壤水分运动规律。特理模型试验验证袁明其与教值模拟　蛄果基本相符，说明将此方法用于果树单坑灌溉入渗问题的研究是行之有拉的，时促进果树灌溉技　术的进一步研究具有重要意义。　美■词：非均质土；ＡＤＩ；教值模拟；蓄水坑灌法　中田分类号；￥２７５．９　文馈标识码；Ａ　蓄水坑灌法是适用于山丘区果园灌溉的一种新　的灌水方法。该灌水方式是在环果树树冠投影范围　内的地面上开挖若于个蓄水坑而对果树进行灌溉，　ＡＦＥＧＯ为蓄水坑。根据多年耕作土壤的剖面分析，　土壤自上而下可分为３层，即表土层，犁底层和心土　层。圈中ｄ　为表土层的厚度，ｄ：为犁底层的厚度，＾　是一种立体灌溉方法。且前该方法无论在理论上还　是技术上研究报道较少，尚有待进一步研究“　］。本　文主要侧重以数值模拟方法研究非均质土单坑人渗　同题，并结合试验，分析蓄水坑灌条件下土壤水分运　动特性，为果园制定合理的灌溉技术设计方案提供　理论依据。　为非均质土的总厚度，Ｒ为模型的最大半径，Ｒ　为　蓄水坑的宽度，＾　为蓄水坑的高度。　１．１基本假定　１）各层土壤具有各向同性；　２）人渗水流为连续介质且不可压缩，在土壤水　流运动过程中，土壤骨架不变形。　１．２数学方程　ｌ数学模型的建立　由于单坑土壤水分运动的物理模型具有轴对称　的特点，放可简化为轴对称的二维平面问题。蓄水坑　剖面圜如图ｌ所示。图中ＡＢ为模型的上边界；ＢＣ　本数学模型采用以负压水头＾为因变量的非饱　和土壤水分运动基本方程　，此方程在柱坐标系（规　定以　轴向下为正，ｒ轴向右为正）下描述为　ｃ（＾）　：了１昙，　）　］＋　０　Ａ　Ｂ　Ｊ　。上　ｒ　Ｄ　Ｇ　’昱　ｈ＾　㈤　卜　．　式中：＾一一负压水头　．　㈩　Ｃ　）——比水容量｝　）一一非饱和土壤导水率。　１．３定解条件　Ｉ．３．１初始奈件　ｆｈ。　（０≤ｒ≤Ｒ，Ｏ≤ｚ＜ｄ１）　丘（ｒ，＃，０）＝　＾０２（ｏ≤ｒ≤Ｒ，　ｌ≤　＜ｄ　ｄ　）．　Ｃ　、　【　。ｓ（０≤ｒ≤Ｒ，ｄ　十ｄ　≤。≤Ⅳ）　（２）　垃　为右，．＋　，．　边界；ＣＤ为下边界ｊＯＺ为蓄水坑的中心轴；　．　一。　其中，＾　＾　＾。　分别为第１，第２，第３层土的韧始　。　…　…　。　…　～～…　作者筒卉：郭文华，男．１９７５年７月生．硕士生．研究方向：节球灌溉．太原．￣３９０２４　收稿日期：２００１—０７—１７　维普资讯 http://www.cqvip.com

第２期　郭文华等：蓄水单坑非均质土土壤水分运动的数值模拟　２ｌ　５　负压水头，计算公式　分别为　进　ｆ丁离散后司得三对角方程如下。　ｈ¨ｌ＝一４　Ｏ２％。“．ｈ。２一一８．２２　”，　嘶．　＆－ｔ　，　十６　＾　＋ｃ　＾譬　，　＝ｄ　，　（９）　ｈ　０ａ一一１　５．６　０ｏｓ　９２．　其中，０　０。　，％分别为第１，第２，第３层土的初始含　其中　＝一　（　．　）　；　水率。　６　＝－ｔ＋Ｒ　二　业（Ｋ　．，）　＋　１．３．２边界条件　１）右边界ＢＣ和下边界ＤＣ距蓄水坑较远，可　Ｒ　（Ｋ　Ｋｉ－：ｌ１　）　］；　认为水分无法到达，故按第一类边界考虑，其负压水　头为初始负压水头。即　Ｃｉ，ｊ一一Ｒ　型（Ｋ　Ｋ　）ｗ；　ＢＣ边界：　ｆｈ。　０≥０，０≤　＜ｄ　）　ｄ　＝Ｒ　（Ｋ　．　Ｋ　Ｌ）“　＾　＋　ｈ（Ｒ，　，ｆ）＝　ｈ。２０≥０，ｄ　≤　＜ｄ１＋ｄ　）　１一Ｒ　［（　Ｋ　ｌＪｔ　）　＋（Ｋ　，　一　）“　］）硅　＋　ｌｈｏ３０≥０，ｄ　＋ｄ２≤　≤Ｚ）；　Ｒ　（Ｋ　．　＋　）”　＾皇　＋　一Ｒ２（Ｋ　一Ｋｆ　一１）；　（３）　Ｒ１一　＝　；　ｉ．ａ　ＤＣ边界：　Ａｔｈ（ｒ，ｈ，　）＝ｈ。　０≥０，０≤ｒ≤Ｒ）．　ｆ４）　＝　．　２）蓄水坑底ＧＥ处加了一个不透水板，使得水　在ｒ方向采用显式、　方向采用踌式，对式（１）　流通量为零，故按第二类边界考虑。即　进行离散后可得三对角方程如下。　ｇ一一Ｋ（＾）夏３ｈ十Ｋ（＾）＝０．　（５）　Ⅱ：．　＾２｝．　十　＾　＋ｃｉ．ｊ＾　ｋ＋　２．　一　（１０）　３）措上边界ＡＢ，侧边界ＧＤ，ＡＦ法线方向的　其中，　水流通量为零，故按第二类边界考虑　即　ｄ：　一一Ｒ　（　！　）“。｛　ＡＢ边界：　＝｛Ｉ十Ｒ　（Ｋ等　Ｋ　ｋ＋＋ｚ１）　一　ｑ一一Ｋ（＾）夏３ｈ＋　（＾）＝０（Ｋ　Ｋｋ　＋一ａ　）　］｝；　：　（６）　ｃ：．　一Ｒ　（Ｋ妨　Ｋ　）　；　ＧＤ，ＡＦ边界；　一ｑ＝一　（＾）警一０．　（７）　．　硝　（麟　Ｋｓ　ｉ　）　＾　．　＋　４）蓄水坑边壁ＥＦ按第一类边界考虑，其负压　｛１一　［　－，２（　ｊ　）ｍ＋　水头始终保持各层土的饱和负压水头。即　ｒｒ　ｚ　（　ｊ　譬　，　）”　］Ｊ＾尊　十　ｆ＾ｎ（ｆ≥０，　ｏ≤　＜ｄ】）　ｈ（Ｒ　Ｌ，ｚ，ｆ）一　ｈ　（　≥０，ｄｌ≤　＜ｄ１十ｄ２）　Ｒ　｝　（Ｋ　：Ｋ譬　）　。＾　ｉ　．　一　ｌｈ６３（　≥０，ｄ　Ｌ＋ｄ２≤２＜ｈ）．　Ｒ　（Ｋ　一Ｋ　）；　（８）　Ａ　ｔ一　式中，ｚ。为注水完毕时，坑中水面到上边界的距离；　＝　；　ｈｍｈ　＾　分　Ⅱ为第１，第２，第３层土的饱和负压水　＝　．　头，计算公式　为　ｈ　一一４．０２％　，ｈ　一一８．２２％　”，　２．２计算网格的划分　ｈ　＝一１５．６　Ｏｋ　”．　将土壤水分运动区域划分为正方形的差分网格　其中，　，％，‰分别为第１，第２，第３层土的饱和含　如图２所示。距离步长Ａｒ＝Ａｚ一２　Ｃｍ，时间采用变　水率。　步长。第一个时间步长取０．０１　ｒａｉｎ．以后每个时间　层的步长比前一个时间层的步长增加１０％，当时间　２数学模型的求解　步长增加到１　ｍｉｎ时，以后其值均采用１　ｒａｉｎ．　２．１方程的离散　２．３边界条件的离散　在ｒ方向采用隐式、　方向采用显式，对式ｒ１）　右边界ＢＣ，由式（３）可得　维普资讯 http://www.cqvip.com

２１６　太原理工大学学报　第３３卷　ｉ　ｉ＋１　：。（　≥０　≤Ｊ≤　…）　凸ｒ　日　侧边界ＡＦ，由式（７）对　取向前差分，在边界结点　ｉ—Ｊ　处列差分方程，可得　ｏ（　≥０，０≤Ｊ≤Ｊ２），（１６）　其中，ｌ，　一蓄水坑中水面处的结点数（１，　一＾　／△　，ｈ　为ｔ时刻的坑中水位）。边壁ＥＦ，由式（８）可得　田２计耳啊格ｊ啊丹不Ｉ圈　ｒ＾６１（　≥０，　／ａｚ≤Ｊ＜ｄｌ／△２），　ｆｈ。１（＾≥Ｏ，Ｏ≤　＜ｄ１／　），　＾　＝　．　（　≥０，ｄ　／△２≤Ｊ＜（　】＋ｄ　）／Ａｚ），　硅　＝　＾　（　≥ｏ，ｄ　／△。≤　＜（ｄ　ｒ＋ｄｚ）／　，　【ｈ　（矗≥０，（　＋ｄ２）／△　≤　≤Ｊ　）．　ｌ＾。３（　≥０，（五＋ｄ２）／ｔａｚ≤Ｊ≤ｍ）．　（１　７）　（１１）　２．４计算过程　其中，ｍ为。方向的最大结点数　／△。）。　本数值计算采用交替方向隐式（ＡＤＤ差分方　下边界ＤＯ，由式（４）可得　法。计算时，在每个结点（　．Ｊ）处建立相应于原方程　砖　＝ｈ　（　≥０，０≤ｉ≤　），　（１２）　（１）的两个差分方程。第１个差分方程对ｒ方向是隐　其中，　为ｒ方向的最大结点数（Ｒ／Ａｒ）　式差分，对＃方向是显式差分，见式（９）；第２个差分　蓄水坑底ＧＥ边界，由式（５）对　取向前差分，　方程对＃方向是隐式差分，对ｒ方向是显式差分，见　在边界结点ｊ＝Ｊ　处列差分方程，可得　式（１　ｏ）　具体求解过程如下。　一　一０　从　＝１开始，对ｉ—Ｊ　的边界为二类边界，补列　一　竺　的差分方程见式（１６）；对ｉ一　的边界为一类边界，　（　≥０，０≤ｉ≤Ｉ１），　边界处的负压水头已知，见式（¨）。则由边界条件及　其中，Ｉ　为蓄水坑边壁处的结点数　／△ｒ）；ｌ，　为蓄　式（９）形成求解方程组可简记为　Ａｈ　＝Ｄ．　（１８）　水坑坑底处的结点数（＾－／△。）．　将上式改写为　其中，Ａ——系数矩阵；　——．　．＾　＋Ｅ．　＾　一ｄｉ．　（矗≥ｏ，０≤ｉ≤　Ｉ）ｔ　求解未知量的列阵；　Ｄ——常数项列阵。　（１３）　然后采用迭代法０　特此非线性方程组转化为线　其　一　ｆ　一　．　性代数方程组　这里，取时段初的Ｃ　．　，　作为时段　￣Ｙａ￣－ＡＢ，由式（６）对差取向前差分，在边界　末ｃ譬　，　翁　的预报值。之后，采用“追赶法”解方程　组（１８），求得时段末所求节点负压水头的第一次选　结点ｊ＝ｏ处列差分方程，可得　代值＾　“”．根据ｃ～＾，　～＾的函数关系，可求得　土壤水分运动参数的校正值。以此参数的校正值作　一雕－　＋耽・一。　为下一次计算的预报值，再解方程组可得时段束所　（＾≥０，Ｉ。≤ｉ≤　）．　求结点负压水头的第二次迭代值ｈｆｔ，“　．重复上述　将上式改写为　步骤，直到所求节点前后两次迭代计算所得负压水　６　＾　＋０．。＾　一ｄⅢ（　≥０，Ｉ１≤ｉ≤　），　头值之差小于所规定的允许误差为止。即应满足　（１４）　ｍ　ＪＩ　‘　－ｌ　１０　≤　≤　≤ｎ　其　一　…．。一ｂ　；以。＝一　．　其中，ｐ——迭代计算次数　侧边界ＧＤ，由式（７）对摹取向前差分，在边界　——允许的相对误差，取　一０．０１．　这样，得到了＾＋１时段末所求结点的负压水头　结点ｉ＝０处列差分方程，可得　＾　｛，＾　，…，舷　，…，雕！｝．．．接着，结合边界条　维普资讯 http://www.cqvip.com

第２期　郭文华等：蓄水单坑非均质土土壤水分运动的数值模拟　２ｌ　７　件，用同样的方法解方程组，依次求出Ｊ＝２，Ｊ一３，　时间分别为３０　ｍｉｎ，１　ｈ，３　ｈ时湿润锋的对比图；图　…沿ｒ向所求结点ｋ＋１时段末的负压水头作为过　６是蓄水坑半径为１６　ｃｍ，坑深为６０　ｃｍ，２　ｈ时负压　渡解，直至Ｊ—ｍ一１．ｋ＋１时段计算完毕后，则改用　水头的等值线对比图。　式（１０）进行　＋２时段的计算。　由图３至图６可以看出，湿润土体大体呈椭球　与　＋１时段相同，从ｉ一１开始，由边界（　—ｄ－　体，在坑底处湿润锋的前移距离最大，且在水平方向　和Ｊ—ｍ）条件和式（１０）形成求解方程组，求解得　的前移距离太于垂直方向，在３０　ｍｉｎ以后湿润锋变　方向所求结点　＋２时段末的负压水头。然后，依次　化逐渐缓慢，到５　ｈ以后大体趋于稳定，这与计算结　求解ｉ一２，　３…．沿ｚ方向所求各结点的值，直至　果基本相符。　ｉ一　一１．这样，就得到了　＋２时段末平面上所求结　ｒ／ｃｍ　点负压水头的数值解。交替使用上面两种格式，直至　求出所要求的各个时段的结果为止。　３计算结果与实测结果的比较分析　根据离散方程编写程序可求得不同时段、不同　结点处的负压水头。为了易与实测数据进行比较分　析，本文选取了有代表性的几组结果绘制成图。图　目０、　３，图４，图５是蓄水坑半径为１　６　ｃｍ，坑深为６０　ｃｍ，０　如　　们　∞　帅　卸　ｒ／ｃｍ　Ｏ　０　２Ｏ　３Ｏ　４Ｏ　圈５　３　ｈ的涅满崔　圈６负压木头为４６．１　ｋＰａ　时的对比圈　４　结　语　本文采用ＡＤＩ差分法对蓄水单坑非均质土的　土壤水分运动进行了数值模拟，并用试验结果对此　予以验证，表明两者基本吻合。另外，对多坑之间水　分人渗的扰动问题应进一步研究，并总结出内在规　律，使之更广泛地应用于生产实践。　圈３　３０ｒａｉｎ的沮桶蕾　围４　１　ｈ的湿润聱　参考文献　Ｃ１］王晓红均质土中蓄水单坑水分运琦的数值模拟与试驻分析ＥＤ］ｌ［磋士论文］太原：太原理工大学、２００１　［２］雷志栋，扬诗秀．谢传森土壤水动力学［Ｍ　北京：清华大学出腹社，１９８８、２９９－３０２　［３］唐梁精，扬秀瑗．苹果园土壤管理与节水灌溉技术［Ｍ］北京：金盾出版社．１９９８．１４６－１　５４　维普资讯 http://www.cqvip.com

２１８　太原理工大学学报　第３３卷　Ｎｕｍｅｒｉｃａｌ　Ｓｉｍｕｌａｔｉｏｎ　ｏｎ　Ｓｏｌｌ　Ｍｏｉｓｔｕｒｅ　Ｍｏｖｅｍｅｎｔ　ｏｆ　Ｎｏｎ—ｈｏｍｏｇｅｎｅｏｕｓ　Ｓｏｉｌ　ｉｎ　Ｓｉｎｇｌｅ　Ｗａｔｅｒ　Ｓｔｏｒａｇｅ　Ｐｉｔ　Ｇｕｏ　Ｗｅｎｈｕａ．Ｓｕｎ　Ｘｉｈｕａｎ　（Ｃｏｌｌｅｇｅ　ｏｆ　Ｃｉｖｉｌ　８Ｌ　Ｅｎｖｉｒｏｎｍｅｎｔ　Ｅｎｇｉｎｅｅｒ。ｒ丁　，了１）　Ｐａｎ　Ｊｕｎｆｅｎｇ　Ｑｉｕ　Ｙａｑｉｎ　（Ｓｈａｎｘｉ　Ｄｅｐａｒｔｍｅｎｔ。ｆ　Ｗａｔｅｒ　Ｃｏｎｓｅｒｖａｎｃｙ）（Ｔｉａ，ｑｉｎ　Ｃｏｒｔｅｇｅ。｛Ｃｉｔｙ　Ｃｏｎｓｔｒｕｃｔｕｒｅ）　Ａｂｓｔｒａｃｔ：Ｔｏ　ｓｅｔｔｌｅ　ｔｈｅ　ｐｒｏｂｌｅｍ　ｏｆ　ｎｏｎ—ｈｏｍｏｇｅｎｅｏｕｓ　ｓｏｉｌ　ｍｏｉｓｔｕｒｅ　ｉｎ　ｗａｔｅｒ　ｓｔｏｒａｇｅ　ｐｉｔ，ａ　ｎｕ—　ｍｅｒｉｃａ１　ｍｏｄｅｌ　ｗａｓ　ｆｏｕｎｄｅｄ．ＡＤＩ（Ａｌｔｅｒｎａｔｉｎｇ　Ｄｉｒｅｃｔｉｏｎ　Ｉｍｐｌｉｃｉｔ）ｉｔｅｒａｔｉｏｎ　ｍｅｔｈｏｄ　ｗａｓ　ｕｓｅｄ　ｔｏ　ｏｂｔａｉｎ　ｔｈｅ　ｓｏｉＩ　ｍｏｉｓｔｕｒｅ　ｍｏｖｅｍｅｎｔ　ｒｕｌｅ　ｉｎ　ｎｏｎ　ｈｏｍｏｇｅｎｅｏｕｓ　ｓｏｉｌ．Ｔｈｅ　ｒｅｓｕｌｔｓ　ｏｆ　ｔｈｅ　ｅｘｐｅｒｉｍｅｎｔ　ａｎｄ　ｎｕｍｅｒｉｃａ１　ｓｉｍｕｌａｔｉｏｎ　ａｒｅ　ｂａｓｉｃａｌｌｙ　ｉｎ　ａＣＣＯｒｄａｎｃｅ．ｈ　ｒｅｖｅａｌｓ　ｔｈａｔ　ｔｈｉｓ　ｍｅｔｈｏｄ　ｉｓ　ｅｆｆｅｃｔｉｖｅ　ａｎｄ　ｈａｓ　ａ　ｓｉｇｎｉｆｉｃａｎｔ　ｖａｌｕｅ　ｆｏｒ　ａｃｃｅｌｅｒａｔｉｎｇ　ｔｈｅ　ｆｕｒｔｈｅｒ　ｒｅｓｅａｒｃｈ　ａｂｏｕｔ　ｏｒｃｈａｒｄ　ｉｒｒｉｇａｔｉｏｎ．　Ｋｅｙ　ｗｏｒｄｓ：ｎｏｎ—ｈｏｍｏｇｅｎｅｏｕｓ　ｓｏｉｌ；ＡＤＩ；ｎｕｍｅｒｉｃａｌ　ｓｉｍｕｌａｔｉｏｎ；　ｍｅｔｈｏｄ　ｏ±ｗａｔｅｒ　ｓｔｏｒａｇｅ　ｌｒｒｔｇａｔｌｏｎ　（编辑：任万森）　（上接第２０９页）　Ａｎ　Ａｐｐｒｏａｃｈ　ｔｏ　Ｄｅｖｅｌｏｐ　Ｐａｒａｍｅｔｒｉｃ　Ｇｒａｐｈ　Ｌｉｂｒａｒｙ　ｏｆ　Ａｒｃｈｉｔｅｃｔｕｒｅ　Ｔｉａｎ　Ｍｅｉｌｉｎｇ，Ｈｏｕ　Ｘｉｎｌｕ　（Ｃｏｌｌｅｇｅ　ｏｆ　Ｃｉｖｉｌ＆ＥｎｖｉｒｏｎｍｅｍＭ　Ｅｎｇｉｎｅｅｒｉｎｇ　ｏｆ丁Ｕ丁）　Ａｂｓｔｒａｃｔ：Ｂａｓｅｄ　ｏｎ　ＡｕｔｏＣＡＤ，ａ　ｇｒａｐｈ　ｌｉｂｒａｒｙ　ｏｆ　ａｒｃｈｉｔｅｃｔｕｒｅ　ｗａｓ　ｄｅｖｅｌｏｐｅｄ　ａｐｐｌｙｉｎｇ　ＶＢＡ．　Ｂｙ　ｍｅａｎｓ　ｏｆ　ｔｈｉｓ　ｓｙｓｔｅｍ，ｐａｒａｍｅｔｒｉｃ　ｄｒａｗｉｎｇ　ｃａｎ　ｂｅ　ｃｏｍｐｌｅｔｅｄ．Ｔｈｅ　ｓｙｓｔｅｍ　ｈａｓ　ｆｒｉｅｎｄｌｙ　ｉｎｔｅｒ—　ｆａｃｅ，ｏｐｅｎ　ａｒｃｈｉｔｅｃｔｕｒｅ　ｔｏ　ｍａｋｅ　ｔｈｅ　ｓｙｓｔｅｍ　ｆｕｒｔｈｅｒ　ｄｅｖｅｌｏｐ　ｃｏｎｖｅｎｉｅｎｔｌｙ　ａｎｄ　ｃｏｍｅｓ　ｔｏ　ｐｒａｃｔｉｃａｌ　ＪｅｖｅＪ．　Ｋｅｙ　ｗｏｒｄｓ：ｐａｒａｍｅｔｒｉｃ　ｄｅｓｉｇｎ；ｇｒａｐｈ　ｌｉｂｒａｒｙ　ｏｆ　ａｒｃｈｉｔｅｃｔｕｒｅ；ＶＢＡ；ＡｕｔｏＣＡＤ　（编辑：任万森）