正冻土中水热耦合数学模型及有限元数值模拟

以數值方法探討人工凍土溫度場與滲流場耦合試驗之模型設計黃任中蔡宜書張睦雄吳博凱雲林科技大學營建工程系摘要當地盤存在地下水流情況，因地下水所造成之熱對流將可能使人工凍土壁不易形成，亦或將已成型之凍土壁產生解凍融化。

為探討此一課題，本文將進行二維縮尺模型之設計討論；利用數值分析軟體（TEMP/W、SEEP/W）模擬溫度場與滲流場作用下，人工冰凍模型凍土成長受水流之影響。

設計考量因素包含：冷凍管徑、冷凍管間距、多管組合之配置、及尺寸效應之影響等。

此模型設計與其試驗結果，將可作為後續數理推導之主要參考依據。

A Study on the Design of Physical AGF Models with Coupled Temperature and Seepage Fields through Numerical ApproachIn the case of flowing groundwater, the seepage would carry with heat which could be detrimental to the formation of frozen soils in an artificial ground freezing (AGF) project. In order to study the seepage effect, the paper herein discusses the design of a physical scaled model, under the influence of coupled temperature and seepage fields, through numerical simulations by commercially available software (TEMP/W and SEEP/W). The design parameters considered include: the size and spacing of freezing pipes, pipe arrangement, and boundary effect, etc. The design and associated test results of the AGF model would provide a sound basis for comparison with the theoretical development and numerical simulation of more complicated cases in the future.一、前言人工冰凍工法（artificial ground freezing method, AGF）之原理係於地盤中依設計間隔鑽孔埋設冷凍液循環管（即冷凍管）後，將冷凍液（如氯化鈣溶液、或稱鹵水、或低溫液態氮）注入循環，使得周圍土壤間隙內的水分冷卻降溫，進而將冷凍管周圍之土壤凍結形成柱狀，使冷凍管與冷凍管間之凍土形成凍結土壁或凍結交圈，以抑制地下伏流，強化土壤強度，達到止水效果；待主體結構物完成後，再使其解凍復原。

基于多孔介质理论的冻土水热迁移耦合模型推导李杨【摘要】在多孔介质理论的基础上,基于非线性达西定律并假设水分迁移过程为单向、可逆及水分迁移过程中无溶质迁移,推导得出了无相变以及考虑相变的水分迁移方程；引入土体传热方程和土骨架质量密度变化方程,过程中考虑了冻土中温度变化,扩散和对流,以及水分相变和温度、质量变化之间的相互影响.联立各方程得到了非饱和冻土水热分布控制方程,建立了非饱和冻土水热迁移耦合模型.文中亦对模型中具体参数的确定方法提出了建议.方程属于非线性偏微分方程,无法得出解析解,须采用数值解法.%In this paper, unsaturated frozen soil moisture and heat transfer coupling model were established based on the theory of porous media. Assuming that the water migration process was one - way and reversible and without solute transport, water transport equation of without and with phase transition were derived based on nonlinear Darcy's law. Then the heat transfer equation in soil and the change equation of quality of soil skeleton density were introduced into the model and in this process, the temperature variation during phase change of water in frozen soil and the interaction between temperature and quality changes were considered. At the same time, the effect of diffusion and convection were taken into account too. Eventually the equation of water and thermal distribution of unsaturated frozen soil was got. And the method for determining the model parameters was given. Model e-quations are nonlinear partial differential equations, only numerical method but analytical solution can be obtained.【期刊名称】《河北工程大学学报（自然科学版）》【年(卷),期】2012(029)003【总页数】4页(P11-14)【关键词】多孔介质;季节冻土;水热迁移;耦合模型【作者】李杨【作者单位】福建工程学院土木工程系福建福州350014;吉林大学建设工程学院吉林长春130026【正文语种】中文【中图分类】TU752;TU411.92土体是一种多孔介质，土体中水分的迁移流动属于多孔介质流体流动[1－2]。

冻土区油气管道周围土壤的热水力三场的数学模型作者：暂无来源：《中国储运》 2011年第11期文／薛洪江吕宏庆摘要：针对穿越冻土区埋地管道存在冻害破坏的安全问题，根据冻土区管道周围实际环境的具体情况，分别综述了管道周围土壤温度场、水热耦合场及管道与水热力三场耦合的数学模型，提出在实际工程建立各个数学模型时需考虑的因素，以期为冻土区埋地油气管道的设计、施工及运行提供参考。

关键词：冻土：管道：温度场：热水耦合场：热水力耦合场：数学模型随着冻土区油气资源的开发，使得穿越冻土区的埋地管道越来越多，其最常见的安全问题是冻害破坏。

为此，在管道设计初期，应对不同地质构造下的土体进行计算分析，预测埋地管道周围土壤冻融过程中热水力三场的变化。

目前在道路和桥梁工程方面，已较全面地开展了冻土工程相关技术的研究工作，但对冻土区油气管道工程的研究较少。

国内外学者在此方面的研究主要集中在对埋地管道周围土壤的温度场、冻土融化圈、水分迁移、应力场和变形场等的分析与计算上。

为更好地解决冻土区埋地油气管道可能遇到的问题，预测埋地管道周围土壤冻融过程中热水力三场的变化，有必要对埋地管道与周围环境的耦合关系进行更深的理论研究，以确保管道的安全运行。

1．埋地管道周围温度场的数学模型温度在土壤中各点分布的集合称为土壤的温度场，温度场是研究冻土的冻融对管道影响的前提。

而埋地油气管道周围冻土温度场的变化规律，受到冻土自身的物性、大气环境温度和管内流体温度的影响。

在冻土区，随着地面温度的波动，将会引起温度梯度作用下热量在土壤中的传输，导致地面热量向下传播，同时埋地管道与周围土壤发生热量交换，温度场数学模型如下。

该模型中，在冻结区和融化区均假设土体在空间上是分层均质的，土体中的含水量较小，均匀分布，而且无源补给和排水作用，各土层含水量稳定。

因此，土体的热物理参数分层稳定，从而可求管道周围的温度场。

实际上，冻土和多年冻土区多具有较强的源汇作用，例如高地下水位、大气降水积水等，水分成了一个不可忽略的影响因素。

冻融土壤水热盐运移规律研究及数值模拟土壤季节性冻融作用是影响旱寒区农业发展的一个重要因素,尤其是冻融期土壤水盐向上运移将会导致越冬期农田地表积盐,从而影响春季播种。

由于冻融土壤水热盐运移过程的复杂性,其演化机理及模拟一直是旱寒区土壤水热盐运移过程研究的瓶颈,开展这一问题的研究对于全面深入理解旱寒区水循环规律和土壤盐渍化形成机理,科学地进行水盐调控和土壤盐渍化的防治具有重要的现实意义和学术价值。

本文以内蒙古河套灌区为背景,通过室内和田间试验,测定了不同盐分条件下的冻融土壤水热运动参数、水盐运移过程和蒸发过程,分析了水热盐运移规律。

改进了水热盐运移模型(CoupModel)并引入GLUE(Generalized Likelihood Uncertainty Estimation)不确定性分析方法,对冻融土壤水热盐运移模拟不确定性进行了深入分析,取得了如下成果：1)通过室内及野外原位观测试验,测定了含盐冻融土壤水热特性参数,分析了水热参数随水分、盐分及温度的变化规律,建立了含盐冻融土壤水热参数估算模型。

试验结果表明,冻融土壤的比热容与热传导率与土壤负温之间分别呈幂函数和线性关系,土壤冻结曲线受土壤初始含水量、盐分种类与含盐量共同影响。

2)通过室内和野外土柱及测坑试验,研究冻融土壤水热盐运移过程和规律。

研究发现,地下水位的高低影响土壤冻结锋的向下推进速度,不同地下水位及不同土壤初始盐分条件对冻融期内土壤蒸发有显著影响。

高地下水位及高盐分含量为土壤蒸发提供了便利。

3)建立了冻融土壤水热盐运移模型--CoupModel,并结合不确定性分析方法GLUE对冻融土壤水热盐运移特性进行了模拟分析,深入探讨了模拟结果的不确定性。

冻融土壤水热运移模拟存在较大不确定性,而提高试验观测手段及完善模型结构则是减小模拟结果不确定性的有效途径。

基于试验结果对含盐土壤冻结特性的认识,考虑盐分对土壤冻结特性的影响,对CoupModel中土壤冻结曲线模型进行了改进。

冻土中热水机械蒸汽的多场耦合研究摘要：本文根据各国学者对冻结土多场的研究成果，对冻结土的多场耦合理论和机理以及冻结过程中温度场、水分场和应力场的动态变化过程进行了分析和研究。

冻土多物理场耦合的研究是一个复杂、多物理、多学科的领域，本文主要从水热蒸汽机械（HTVM）场耦合的方面进行了综述。

本综述有助于促进对冷区土壤冻结过程和冻结过程中冻土耦合机理的研究，促进对土壤冻结过程中多场耦合动态过程的深层、多维理解。

1.介绍冻融沉降是寒冷地区冻土最常见的冻土破坏。

这主要是由于冻土的温度、湿度、应力和浓度场的变化，以及多物理场之间的相互作用（Mu，1987）。

冻土的各种霜冻问题本质上源于多孔介质中的多相耦合（包括固体、液体、气体和热）（Li，2001）。

季节性冻土的土壤稳定性主要受冻融循环中传热、水分迁移和相变化的相互作用和相互影响的控制。

例如，路基的温度、湿度和应力场都是动态变化的。

这些油田之间的耦合效应是造成许多冻害问题的直接原因（Lai，1999）。

因此，有必要系统地回顾冻土多物理场耦合背后的机制。

冻土的多物理场耦合是一个同时考虑多个物理场的复杂的多学科研究部门。

冻土多物理场耦合的研究进一步分为以下几个部分：热水（TH）耦合、热水机械（THM）耦合、热水蒸汽机械（THVM）耦合和热水盐力学（THSM）耦合。

此外，一些学者还致力于研究岩石的热-水-化学-机械（THCM）耦合（Su，2010），这是一个复杂而动态的过程。

近年来，有关冻结场的学术兴趣一般集中在宏观强度性质和热-水-力学耦合本构模型上。

考虑了晶体独特的张力行为和压力熔化，以及冰水相变的结晶动力学模型。

然而，由于缺乏特殊设备，目前进行的本构实验很少考虑微观变形的机理。

因此，迫切需要对微小量表进行本性调查。

2.冻土多物理场耦合技术的研究现状2.1冻土高温耦合技术的发展提出了浅层黄土的TH耦合模型，模型结果与实验结果一致。

通过该模型估算的水分和温度的动态变化，验证了参数选择的可行性和预测浅层冻土水热动力行为的准确性[1]。

受高温混凝土的水,热,力耦合的galerkin有限元格式受高温混凝土的水热力耦合的Galerkin有限元格式高温混凝土在各种工程领域中广泛应用，例如核电站和高温炉等。

然而，在高温环境下，混凝土会受到热、水和力的耦合作用，从而引起各种不稳定问题，如开裂和变形等。

为了更好地理解和解决这些问题，研究人员开发了各种数学模型和数值方法，其中Galerkin有限元方法是一种常用的数学工具。

Galerkin有限元方法是一种数值近似方法，旨在通过将连续问题离散化为一组代数方程来解决物理问题。

在高温混凝土中，水、热、力的耦合作用对结构的稳定性和性能有着重要影响。

因此，研究人员提出了一种基于Galerkin有限元方法的水热力耦合模型，来模拟高温混凝土中各个物理过程之间的相互作用。

在建立数学模型时，我们需要考虑高温混凝土在水、热、力三个方面的耦合行为。

首先，水的存在会引起混凝土的膨胀和吸湿，从而导致温度和应力场的变化。

因此，我们需要引入水的吸湿模型和渗流模型来描述水在混凝土中的分布和流动。

其次，高温引起混凝土内部温度场的变化，进而影响混凝土的力学性能。

因此，我们需要考虑温度场的传导和对流传热过程。

最后，由于混凝土的热膨胀和水分迁移，会引起内部应力的分布和变化。

因此，我们需要采用弹性力学模型来描述混凝土的机械行为。

根据以上分析，我们可以得到高温混凝土水热力耦合的Galerkin有限元模型。

首先，我们需要建立水的吸湿和渗流模型，采用质量守恒和达西定律来描述水的分布和流动。

其次，我们需要建立温度场的传导和对流传热模型，采用热传导方程和对流传热方程来描述温度场的变化。

最后，我们需要建立混凝土的弹性力学模型，采用应力平衡和胡克定律来描述混凝土的力学行为。

在数值实现方面，我们需要在有限元网格上离散化上述数学模型，并使用Galerkin方法得到代数方程组。

然后，我们可以使用求解器来解决这个方程组，得到高温混凝土中水、热、力三个物理场的数值解。

第37卷第6期2023年12月水土保持学报J o u r n a l o f S o i l a n d W a t e rC o n s e r v a t i o nV o l .37N o .6D e c .,2023收稿日期:2023-04-28资助项目:国家自然科学基金项目(41977060) 第一作者:蒋瑶钰(1998 ),女,安徽安庆人,在读硕士研究生,主要从事土壤侵蚀与水土保持研究㊂E -m a i l :j i a n g y a o y u 0503@163.c o m 通信作者:王彬(1983 ),男,河南新乡人,教授,博士,博士生导师,主要从事土壤侵蚀与水土保持研究㊂E -m a i l :w a n g b i n 1836@b jf u .e d u .c n 基于C O M S O L 模型的黑土农耕地剖面冻结水热耦合分析蒋瑶钰,王彬,陈祖明,王文刚(北京林业大学水土保持学院,北京100083)摘要:为量化分析典型黑土农耕地在冻结过程中的剖面水热耦合特征,以东北典型黑土区耕层土壤为研究对象,通过室内单土柱模拟试验,并结合C OM S O L M u l t i p h ys i c s 模型P D E 模块,实现黑土冻结过程中温度场与水分场的耦合求解㊂结果表明:(1)C OM S O L 模型可实现黑土剖面冻结过程的水热耦合分析,温度场(R 2=0.83,R M S E =0.91ħ)与水分场拟合效果(R 2=0.88,R M S E=0.02c m /c m3)均较好,可满足拟合精度需要;(2)黑土农耕地剖面温度的变化速率呈先快速冷却后缓慢趋于稳定的变化趋势,且越接近冷端,土壤温度下降至稳定阶段的速率越快;(3)黑土农耕地剖面水分场在冻结过程中呈明显的水分重分布现象,且在冻结峰面呈现暖端聚集速率快[0.00093~0.00100c m /(c m 3㊃h )]而冷端反之[0.00034~0.00042c m /(c m3㊃h )]的冰晶聚集现象㊂研究结果可为黑土冻结水热运移规律分析提供技术支持,并为黑土冻融 水力复合侵蚀机理研究和剖面水热动态模拟提供理论基础㊂关键词:水热耦合;温度场;水分场;冻融侵蚀中图分类号:S 157.1 文献标识码:A 文章编号:1009-2242(2023)06-0187-07D O I :10.13870/j.c n k i .s t b c x b .2023.06.024H y d r o t h e r m a l C o u p l i n g A n a l y s i s o f F r e e z i n g Pr o c e s s o fB l a c kS o i l F a r m l a n dP r o f i l eB a s e do nC O M S O L M o d e lJ I A N G Y a o y u ,WA N GB i n ,C H E NZ u m i n g ,WA N G W e n g a n g(S c h o o l o f S o i l a n d W a t e rC o n s e r v a t i o n ,B e i j i n g F o r e s t r y U n i v e r s i t y ,B e i j i n g 100083)A b s t r a c t :T o q u a n t i t a t i v e l y a n a l y z et h ec o u p l i n g c h a r a c t e r i s t i c so fw a t e ra n dh e a t i nt h e p r o f i l eo f t y pi c a l b l a c ks o i l c u l t i v a t e d l a n dd u r i n g t h e f r e e z i n gp r o c e s s ,t h e t o p s o i l i n t h e t y p i c a l b l a c k s o i l r e g i o n i nN EC h i n a w a s t a k e na s t h e r e s e a r c ho b j e c t .T h r o u g h i n d o o r s i n g l e s o i l c o l u m n s i m u l a t i o n e x pe r i m e n t s a n d c o m b i n a t i o n w i t hC OM S O L M u l t i p h y s i c s m o d e lP D E m o d u l e ,t h ec o u p l i n g s o l u t i o nof t e m pe r a t u r ef i e l da n d m o i s t u r e f i e l dd u r i ng b l a c ks o i l f r e e z i n gpr o c e s sw a s a c h i e v e d .R e s u l t s s h o w e d t h a t :(1)T h eC OM S O L m o d e l c o u l d r e a l i z e t h eh y d r o t h e r m a l c o u p l i n g a n a l y s i s o f t h e f r e e z i n g p r o c e s s o f b l a c k s o i l p r o f i l e ,a n d t h e f i t t i n g e f f e c t o f t e m pe r a t u r ef i e l d (R 2=0.83,R M S E=0.91ħ)a n d w a t e r f i e l d (R 2=0.88,R M S E=0.02c m /c m 3)w a sg o o d ,whi c hc o u l d m e e t t h e f i t t i n g a c c u r a c y r e q u i r e m e n t s .(2)T h et e m p e r a t u r ec h a n ger a t eo fb l a c ks o i l c u l t i v a t e d l a n d p r o f i l ee x h i b i t e dat r e n do f r a p i dc o o l i n g f o l l o w e db y as l o wa p p r o a c ht os t a b i l i t y ,a n dt h e c l o s e r t o t h e c o l d e n d ,t h e f a s t e r t h e r a t e a tw h i c h s o i l t e m p e r a t u r e d e c r e a s e d t o t h e s t a b i l i t y s t a g e .(3)T h e f r o z e nb l a c k s o i l c u l t i v a t e d l a n d p r o f i l e u n d e r w e n t s i g n i f i c a n tw a t e r r e d i s t r i b u t i o n d u r i n g t h e f r e e z i n g p r o c e s s .A n i c e c r y s t a l a g g r e ga t i o n p h e n o m e n o nw a s ob s e r v e d a t t h ew a r me n d ,w h ic h o c c u r ed a t a f a s t r a te [0.00093~0.00100c m /(c m 3㊃h )],w h i l e t h e c o l d e n d s h o w e d t h e o p p o s i t e p h e n o m e n o n [0.00034~0.00042c m /(c m 3㊃h )].T h e r e s e a r c h r e s u l t s c a n p r o v i d e t e c h n i c a l s u p p o r tf o r a n a l y z i ng th ew a t e r h e a t t r a n s po r t l a wo f f r o z e nb l a c k s o i l ,a sw e l l a s a t h e o r e t i c a l b a s i s f o r t h e s t u d y o f t h em e c h a n i s mo f f r o z e n t h a w i n g a n dh y d r a u l i c c o m p o s i t e e r o s i o no f b l a c ks o i l p r o f i l e s a n dd yn a m i c s i m u l a t i o no fw a t e r -h e a t i n p r o f i l e s .K e yw o r d s :h y d r o -t h e r m a l c o u p l i n g ;t e m p e r a t u r e f i e l d ;m o i s t u r e f i e l d ;f r e e z e -t h a we r o s i o n东北黑土区具有明显的季节性冻融特征,冻融 水力复合侵蚀现象明显,导致土壤退化㊁粮食减产,严重威胁国家粮食安全[1]㊂黑土冻结过程的量化表征是冻融 水力复合侵蚀过程描述的关键难点㊂土壤冻结过程中剖面水热状况存在复杂的耦合变化,造成剖面土壤含水量发生重分布[2],继而改变次年地表水文条件,间接影响次年春季农作物的生产力[3],增加水土流失风险[4]㊂研究黑土冻结过程水热耦合特征和机理,可为冻融过程的水热动态模拟提供理论基础,并为构建冻融 水力复合土壤侵蚀机理模型提供科学依据㊂目前,大多学者[5-7]通过野外原位观测与室内模拟的手段,研究不同冻结强度㊁土地利用方式和覆盖条件下(积雪㊁秸秆)黑土坡面的水热动态变化过程㊂部分学者[8]通过控制不同初始含水率㊁降温幅度和地下水条件下的单向冻融试验,定量分析黑土剖面水分迁移的主要影响因素㊂近年来,对于东北黑土区冻结过程中土壤温度㊁水分的变化特征已取得丰富成果;但受限于试验条件等因素影响,土壤剖面水热耦合过程的定量表征仍相对薄弱,一定程度上限制冻融 水力复合侵蚀过程机制的研究㊂数值模型模拟是研究冻土水热耦合规律的重要方式,已被广泛地应用在冻土路基㊁寒区隧道及冻土边坡等工程领域[9-11]㊂随着有限元法的引入,数值模型模拟的精度不断提升㊂相关学者引入水分迁移热熵[12]㊁相变过程水分对流[13],以及冰相水分运移阻抗效应等,提出水热场偏微分方程有限元差分形式[14]㊂P h i l i p等[15]基于热平衡及多孔介质水分黏性流动理论,首次提出土壤水热耦合运动模型㊂部分学者选用徐学祖方程[12]㊁固液比(B I)[16-17]㊁含冰率 温度关系[18]等为联系方程求解水热耦合高度非线性偏微分方程组㊂其中,基于C OM S O L模型的系数型偏微分模块(P D E)被广泛运用于季节性冻土区高度非线性水热耦合偏微分方程的求解㊂东北黑土区是我国重要的粮食生产基地,农耕地土壤冻融 水力复合侵蚀现象明显㊂针对耕作黑土,探究单向冻结模式下黑土温度变化和水分迁移特征,对土壤冻融侵蚀过程机制的深入理解具有重要作用㊂选取东北典型黑土区耕层土壤为研究对象,通过室内单土柱模拟试验,在评价冻土水热耦合模型适用性的基础上,实现C OM S O L M u l t i p h y s i c s模型的自定义数学模块(P D E)对水热场方程进行求解,实现黑土冻结过程中温度场和水分场数值模拟㊂研究结果可为黑土冻结水热运移规律分析提供技术支持,并对进一步研究该区土壤冻融侵蚀过程机制认识等提供理论依据㊂1材料与方法1.1试验土壤与试验设计试验土样采自东北典型黑土区黑龙江省齐齐哈尔市克山县克山农场(125ʎ08'00ᵡ 125ʎ37'00ᵡE,48ʎ12'00ᵡ 48ʎ23'00ᵡN)耕层土壤㊂鉴于冻融 水力复合侵蚀多发生在扰动强烈的耕层,选取研究区15c m耕层土壤为研究对象㊂研究区属季节性冻土区,土壤冻结多发生在11月中旬至翌年3月,平均气温-10~-30ħ,平均地温0~5ħ㊂供试土壤质地为粉质黏土(黏粒29.4%,粉粒61.3%,砂粒9.3%)㊂试验于2021年8 12月在北京林业大学水土保持学院水土保持重点实验室进行,试验共设计3组(预设初始含水率梯度为0.20,0.25,0.30c m/c m3),共设置7次处理(表1)㊂参照自然条件下冻结期实测的黑土农耕地多年平均土壤温度及含水率,设置控制条件试验模拟的冷端温度为(-10ʃ2)ħ,暖端温度为(3ʃ0.5)ħ,试验土柱初始温度为3ħ,冻结方向为垂直向下的单向冻结(表1)㊂2组独立数据分别用于模型参数率定(T1~T6)和验证(T7)㊂表1试验设计试验编号初始重量含水率/(c m3㊃c m-3)暖端温度/ħ冷端温度/ħ试验土柱初始温度/ħT10.2053.0-5.03.0T20.2253.0-9.73.0T30.2732.8-9.43.0T40.2942.7-10.23.0T50.2952.8-8.23.0T60.2792.9-14.53.0T70.2583.0-10.03.0 1.2研究方法1.2.1室内单土柱模拟试验试验参照‘土工试验方法标准“[19](G B/T50123-1999)进行制备㊂土样过2mm筛,按设计含水量分层填装至亚克力透明制样筒(土样高15c m,直径10c m,容重1.28g/c m3),随后密封静置于恒温条件下(3ʃ0.5ħ)12h使土水分均匀分布㊂试验装置由外循环温度控制装置(S H D-6000,中国)㊁高低温恒温液浴循环两用槽(X T5218-B12)㊁边界恒温板(-10ħ冷端㊁3ħ暖端)㊁热电耦式温度探头(J K-24C)和试验土柱组成,整个试验过程采用保温隔热材料包裹试验土柱以保证一维垂直冻结的温度梯度条件(图1)㊂试验土柱内等距设置针式温度探头,每30s测量1次㊂冻结结束稳定阶段,采用烘干法[20]分层测定对应土样重量含水量㊂1.2.2冻土水热耦合微分方程组(1)基本假设㊂根据季节性冻土的主要物理过881水土保持学报第37卷程,设置假设为:①土壤是不可压缩的㊁均匀的㊁各向同性的介质㊂②孔隙水渗流服从达西定律㊂③限制边界条件下季节性冻土发生垂向一维水热运动㊂④土中盐分的影响可忽略㊂⑤冻结过程中主要水分迁移为液态(忽略凝华影响)㊂图1 垂直-维冻结土柱试验示意(2)温度场控制方程㊂基于傅里叶定律[21]并考虑二维水热耦合问题,将相变潜热作为热源进行处理,得到冻土热传导微分方程为:ρC (θ)∂T ∂t =λ(θ)Ñ2T +L ㊃ρI ∂θI∂t(1)式中:ρ为土的密度(k g /m 3);ρI 为冰的密度(k g /m 3);C 为热容量[J /(k g㊃ħ)];λ为导热系数[W /(m ㊃ħ)];L 为相变潜化热,取值334.56(k J /k g );T 为温度(ħ);t 为时间(s );体积含水量为θ=θu +ρI /ρw ㊃θI (k g /m 3),其中,θu 为未冻水的体积含量(%);ρw 为水的密度(k g/m 3)㊂(3)水分场控制方程㊂基于非饱和融土水分运动规律[22],引入水冰相变项,得到冻土中水分迁移R i c h a r d s 方程:∂θu ∂t +ρI ρw ∂θI ∂t=Ñ[D (θu )Ñθu +k (θu )](2)考虑土壤冻结过程中冰晶对水流的阻塞作用[14],冻土中水的扩散率计算公式为:D (θu )=k (θu )c (θu)㊃I (3)I =10-10θI(4)式中:k (θu )为土体渗透率(m /s );c (θu )为比水容量(1/m );I 为描述土中冰对水分迁移阻碍作用的阻抗因子㊂(4)联系方程㊂冻土温度场和水分场中包含温度㊁孔隙冰体积和未冻水体积3个未知量,故须引入联系方程才能实现水热耦合模型求解㊂选用固液比B I 为耦合项,固液比是在徐学祖方程[16]基础上进行的数学形式转换,本质为未冻水含量和温度的关系,在方程求解时具有更易计算的优点,用于表示土中冰与未冻水的体积之比,计算公式为:B I =θIθu =1.1(T T f )B -1 T <T f 0 T >T f ìîíïïï(5)式中:T f 为土体冻结温度(ħ);固液比相关系数B 为随土质与含盐量变化的常数,一般为0.56[23]㊂选用v a nge n u e h t e n (V G )(1980)滞水模型和G a r d n e r 渗透系数模型(1958)[22],并将冻土的相对饱和度(S )定义为:S =θu -θrθu -θr(6)式中:θr 为残余含水率(%);θs 为饱和含水率(%)㊂联立温度场微分方程(1)㊁水分场微分方程(2)以及联系方程B I ,即得到冻土水热耦合求解方程组:ρC (θ)∂T ∂t =λ(θ)Ñ2T +L ㊃ρI ∂θI ∂t∂θu ∂t +ρI ρw ∂θI ∂t =Ñ[D (θu )Ñθu +k (θu )]B I =θIθuìîíïïïïïïïï(7)综上,冻土水热耦合方程中所需的参数见表2㊂981第6期 蒋瑶钰等:基于C OM S O L 模型的黑土农耕地剖面冻结水热耦合分析表2 函数参数符号指标物理意义公式参数文献来源k (θu )渗透率单位水力梯度下流量k s ㊃S l [1-(1-S 1/m )m ]2L ,m G a r d n e r (1958)[22]I 阻抗因子冰对未冻水迁移阻滞作用10-10θIθIT a y l o r (1978)[14]c (θu )比水容量基质势变化引起含水量变化a ㊃m /(1-m )㊃S 1/m (1-S 1/m )ma ,mV G (1980)[22]D (θu )扩散率传播温度的速率k (S )/[(θs -θr )/c (S )]㊃I L uN (2012)[22]1.3 数据处理及模型模拟适用性评价分析采用E x c e l 2003软件对数据进行统计分析,O r i gi n 2019软件作图㊂采用决定系数(R 2)㊁均方根误差(R M S E )评价模型模拟效果,当R 2越接近1,R M S E 越接近0时,证明其误差越小,其计算公式为:R 2=1-ði (S i -M i )ði (S i -M i )2(8)R M S E =1nðn i =1(S i -M i )2(9)式中:S i 为模拟值;M i 为实测值;`M i 平均值;n 为实测样本数(个)㊂2 结果与分析2.1 基于C O M S O L 二次开发的冻土水热耦合模型构建通过C OM S O L 中系数型偏微分模块(P D E )进行二次开发,建立冻土水热耦合数值模型㊂将温度场及水分场高度非线性方程组(式7)化解成系数型标准形式后,完成P D E 模块参数输入㊂系数型偏微分方程的标准形式为:e a∂2u ∂t2+d a ∂u∂t +Ñ(-c Ñu -αu +γ)+βÑu +a u =f (10)式中:u 为所求变量;e a 为质量系数;d a 为阻尼系数;c 为扩散系数;q 为边界上的吸收系数;β为对流系数;ɑ为吸收系数;γ为保守通量源项;f 为源项㊂将冻土水热耦合求解方程组式(7)中温度场及水分场方程分别转化为系数偏微分方程组的形式,并与系数偏微分方程标准形式进行对比[17],分别可得温度场偏微分方程参数为:d a =ρC (θ)c =λ(θ)f =L ㊃ρI ㊃(θs -θr )㊃(∂B (T )∂t ㊃S +B (T )㊃∂S∂t )ìîíïïïï(11)水分场偏微分方程参数为:d a =1+ρI ρw B (T )c =D (S )γ=-K (S )a =ρI ρw ㊃∂B (T )∂tìîíïïïïïïïï(12)将式(11)㊁式(12)参数对应输入P D E 模块,完成温度场方程与水分场方程的联立㊂此外,固液比方程参考白青波等[16]在黏土冻结试验数值模型模拟中算法,采用模拟试验T 1~T 6数据组进行模型参数率定(表1),并通过迭代法[24]得到适用于研究区的固液比方程:B I =1.744T0.56-1.013(13)将温度场以及水分场方程的P D E 模块联立固液比方程(公式13),即可耦合求解水热联合方程组,进而建立冻土水热耦合模型㊂基于前期试验并集成分析现有研究[25]成果,确定温度场和水分场模拟所需参数(表3),其中饱和土渗透系数为9.62ˑ10-5m /s ,土密度为1.28g /m 3,饱和含水率为0.40c m 3/c m 3,固液比相关系数为0.56,相变潜化热为33456J /k g[16,26-27]㊂同时,选用单土柱冻结试验另一组独立数据(T 7)进行模型验证㊂表3 参数汇总参数数值单位土的比热容0.85J /(k g㊃K )水的比热容4200.00J /(k g㊃K )冰的比热容2100.00J /(k g㊃K )土的导热系数0.95W /(m ㊃K )水的导热系数0.62W /(m ㊃K )冰的导热系数2.31W /(m ㊃K )相变潜化热33456.00J /k g 冰密度918.00k g /m 3水密度1000.00k g /m 3土密度1280.00k g /m 3饱和含水率0.40c m 3/c m 3残余含水率0c m 3/c m3土体冻结温度272.46K与固液比相关的系数0.56无量纲饱和土的渗透系数9.62ˑ10-5m /sV G 模型参数l0.50无量纲V G 模型参数m 0.22无量纲V G 模型参数a 2.591/m 2.2 单土柱土壤水热耦合特征2.2.1 温度场模拟 黑土剖面温度的变化速率随时间变化表现为快速冷却阶段(F C )㊁缓慢冷却阶段(S C )和稳定阶段(S S )3个阶段(图2㊁图3)㊂各土层(3,6,9,12c m )实测值稳定温度分别为-0.1,-2.2,-4.4,-7.2ħ㊂土层高度3c m 接近土柱暖端,12c m 接近土柱冷端,温度变化速率随土柱高度的变化091水土保持学报 第37卷而变化,高度越接近冷端,温度下降到稳定阶段的速率越快,即温度到达稳定所需时间越少㊂黑土剖面冻结区和未冻结区温度沿土柱均呈接近线性分布(图4)㊂冻结早期,冻结区存在较大的温度梯度;与未冻区相比,冻结区的导热系数更大,热容更小㊂冻土水热耦合模型的温度场模拟值效果较好(R2=0.83,R M S E=0.91ħ),整体拟合趋势与实测趋势基本一致(图2),平均低估4.9%㊂温度场拟合低估率随土壤剖面高度呈明显差异,靠近控温边界(冷端12c m㊁暖端3c m)处的低估率较小,而远离控温边界处的土壤温度低估幅度明显增加,低估率增加为边界处的1.9~2.1倍㊂综上,基于C OM S O L模型P D E解算的土壤剖面温度的拟合效果达到较好水平,可用于黑土剖面传热过程的模拟㊂图2冻结过程试验土柱各层次土壤温度变化图3温度场随时间变化(模拟值)图4土柱的模拟温度随时间的分布2.2.2水分场模拟冻结作用使黑土剖面土壤水分发生明显的重分布现象(图5)㊂冻结峰面处存在明显的冰晶聚集现象,暖端聚集速率快[0.00093~0.00100c m/(c m3㊃h)],冷端聚集速率慢[0.00034~ 0.00042c m/(c m3㊃h)],最大含水量(高度4.5c m)实测为0.35c m㊃c m3,较初始含水率提高9.1%㊂此外,剖面水分迁移率变化趋势与温度变化相反,暖端水分聚集速率为0.00093~0.00100c m/(c m3㊃h),而冷端较其慢2.2~2.9倍㊂随着冻结时间推移,冻结锋面逐渐下移,水分迁移量逐渐增大(图6)㊂初始含水率为0.2575c m/c m3,冷端温度为-10ħ时,冻结8.2,30,48,60h时,最大含水率分别为0.24,0.28,0.33,0.35c m/c m3,冻结锋面高度分别为9.64,4.96,4.82,4.28c m,冻结锋面高度随时间推移逐渐稳步下移㊂此外,随着冻结时间推移,冻结锋面以下含水率逐渐减小㊂冻结8.2,30,48,60h时,冻结锋面以下即未冻结区,含水率平均为0.17,0.17,0.15,0.15c m/c m3,相对于初始含水率分别减少33.05%,33.98%,40.50%, 41.74%,水分迁移量逐渐增大㊂在温度梯度影响下,冻结与未冻结部分间的土水势差变大,驱使随着水从未冻结区逐渐向冻结区迁移,冻结锋面逐渐下移;同时随着冻结时间推移,冻结锋面迁移量逐渐增大㊂冻土水热耦合模型的水分场模拟值较好,整体拟合趋势与实测趋势整体相符(R2=0.88,R M S E= 0.020c m/c m3)(图5)㊂靠近冷端附近,模拟值与实测值总含水量平均差值为0.034c m/c m3㊂图5试验结束时水分场分布图6不同冻结时刻含水率模拟值沿土柱高度分布3讨论土壤剖面温度梯度对控制条件下单向冻结土柱温度变化速率有显著影响㊂黑土土柱温度变化速率呈现明显的3个变化阶段,且越接近冷端变温速率稳191第6期蒋瑶钰等:基于C OM S O L模型的黑土农耕地剖面冻结水热耦合分析定时间越短㊂Z h o u等[28]㊁H u a n g等[29]通过C OM-S O L和完全耦合的热 水力 力学(T HM)模型模拟饱和冻结土柱温度场的相关研究得到相近的温度变化趋势㊂剖面温度临近0ħ时,由于孔隙水相变成冰释放潜热,造成温度随时间的变化速率相对变缓[29]㊂温度场拟合效果分析表明,远离控温边界处的土壤温度低估幅度大于靠近控温边界处的低估率,低估率增加为边界处的1.9~2.1倍㊂其原因主要是实体模型和模型模拟间的径向传热条件存在差异[30],模型模拟中的径向传热边界为绝热条件,实体模型无法完全避免热交换现象㊂冻结结束后水分场发生明显的重分布现象,该现象与张明礼等[13]通过C OM S O L软件反演的壤土单向冻结试验结束时水分场结果相似㊂冻结锋面处冰晶聚集是由于冻结层向下移动缓慢,使得暖端附近未冻土层有足够的时间迁移更多的水分到冻结锋处[28],而冷端附近的冻土层由于受到冰的阻滞作用致使迁移量减少㊂在土壤冻结过程中,一部分土壤水分由液态变为固态,土壤未冻水量减少,基质势降低,冻结与未冻结部分间的土水势差变大,随着冻结深度的发展,土壤未冻结区的液态水不断向冻结锋面运移,部分在冻结界面附近冻结,部分受到冰的阻滞作用而滞留,两者共同作用导致最大含水量出现在冻结峰面附近[31]㊂同时,在对水分场拟合效果分析中,冷端附近误差较大,一方面由于冷端土样温度梯度大,水分迁移较迅速,冻结速率快,导致液相水分在此聚集时间较短,易引起试验误差[32];另一方面,由于未考虑非饱和土样中的气相凝华现象㊂采用C OM S O L模型的自定义数学模块(P D E)对水热场方程进行求解,实现黑土冻结过程中温度场和水分场数值模拟,温度场与水分场拟合效果均较好㊂与白清波等[16]研究结果相似,但在固液比参数选取上不同,其主要原因为土壤质地不同导致冻结过程中未冻水含量和温度传递特性改变,致使固液比参数存在差异[31,33]㊂此外,模型计算结果整体在精度范围内,但如何缩小误差仍存在一定可提升空间㊂参数准确性将直接影响模拟结果的有效性,土壤冻结过程是土壤性质与外界条件等多种因素综合作用的结果,其水热耦合变化十分复杂[26],在相互作用中相关参数也有所变化,为方便模型求解计算,整个冻结过程中采用相同的参数,可以结合野外实测数据优化模型各土层参数,是未来待继续研究的主要问题㊂4结论(1)C OM S O L模型适用于黑土剖面冻土的水热耦合分析㊂基于模型模拟的温度场(R2=0.83, R M S E=0.91ħ)和水分场拟合效果(R2=0.88, R M S E=0.02c m3/c m3)均可达到满意水平㊂(2)黑土农耕地剖面温度的变化速率随时间变化表现为快速冷却阶段(F C)㊁缓慢冷却阶段(S C)和稳定阶段(S S)3个阶段,剖面高度越接近冷端,温度下降到稳定阶段的速率越快㊂(3)黑土农耕地剖面在冻结过程中呈明显的重分布特征,且在冻结峰面表现出暖端聚集速率快[0.00093~0.00100c m/(c m3㊃h)]而冷端反之[0.00034~0.00042c m/(c m3㊃h)]的冰晶聚集现象㊂参考文献:[1]张光辉,杨扬,刘瑛娜,等.东北黑土区土壤侵蚀研究进展与展望[J].水土保持学报,2022,36(2):1-12. [2]王子龙,付强,姜秋香,等.季节性冻土区不同时期土壤剖面水分空间变异特征研究[J].地理科学,2010,30(5): 772-776.[3] D a i ZG,F e i LJ,H u a n g DL,e t a l.C o u p l i n g e f f e c t s o fi r r i g a t i o na n dn i t r o g e n l e v e l so n y i e l d,w a t e r a n dn i t r o-g e nu s ee f f i c i e n c y o fs u r g e-r o o ti r r i g a t e d j u j u b ei n as e m i a r i dr e g i o n[J].A g r i c u l t u r a l W a t e r M a n a g e m e n t, 2019,213:146-154.[4] B a nYY,L e i T W,L i uZQ,e t a l.C o m p a r a t i v e s t u d y o fe r o s i o n p r o c e s s e s of t h a w e da n dn o n-f r o z e ns o i l b y c o n c e n-t r a t e dm e l t w a t e r f l o w[J].C a t e n a,2017,148:153-159. [5]张科利,彭文英,王龙,等.东北黑土区土壤剖面地温和水分变化规律[J].地理研究,2007,26(2):314-320. [6]张少良,沈庆松,王曜,等.不同冻结强度下容重和含水量对黑土剖面水分变化特征影响[J].东北农业大学学报,2016,47(12):48-55.[7]胡伟,张兴义,严月.不同土地利用方式下冻融期黑土水热过程观测研究[J].土壤与作物,2018,7(3):312-323.[8]刘红希,范昊明,许秀泉.黑土冻融过程中水分垂直迁移模拟研究[J].水土保持学报,2021,35(1):169-173. [9] R u iD H,Z h a i JB,L iG Y,e t a l.F i e l de x p e r i m e n t a ls t u d y o ft h ec h a r a c t e r i s t i c so fh e a ta n d w a t e rt r a n s f e rd u r i n g f r o s t he a v i n g[J].C o l dR e g i o n s S c i e n c e a n dT e c h-n o l o g y,2019,168:e102892.[10] L a iY M,L i uSY,W uZ W,e t a l.N u m e r i c a l s i m u l a-t i o n f o r t h e c o u p l e d p r o b l e mo f t e m p e r a t u r e a n ds e e p-a g e f i e l d s i n c o l dr e g i o nd a m s[J].J o u r n a l o fH y d r a u l i cR e s e a r c h,2002,40(5):631-635.[11]W a n g C,M aZP.M a t h e m a t i c a lm o d e l a n dn u m e r i c a ls i m u l a t i o no fh y d r o t h e r m a lc o u p l i n g f o r u n s a t u r a t e ds o i l s u b g r a d e i n t h e s e a s o n a l f r o z e nz o n e[J].I O PC o n-f e r e n c e S e r i e s:E a r t h a n d E n v i r o n m e n t a l S c i e n c e,291水土保持学报第37卷2021,719(3):e032042.[12] L i uZ,Y uX.C o u p l e d t h e r m o-h y d r o-m e c h a n i c a lm o d e l f o rp o r o u sm a t e r i a l su n d e rf r o s ta c t i o n:T h e o r y a n di m p l e-m e n t a t i o n[J].A c t aG e o t e c h n i c a,2011,6(2):51-65.[13]张明礼,郭宗云,韩晓斌,等.基于C OM S O L M u l t i p h y s-i c s数学模块的冻土水热耦合分析[J].科学技术与工程,2018,18(33):7-12.[14] T a y l o rGS,L u t h i n JN.A m o d e l f o r c o u p l e dh e a t a n dm o i s t u r et r a n s f e r d u r i n g s o i lf r e e z i n g[J].C a n a d i a nG e o t e c h n i c a l J o u r n a l,1978,15(4):548-555.[15] P h i l i p JR,D eV r i e sD A.M o i s t u r em o v e m e n t i n p o r o u sm a t e r i a l su n d e rt e m p e r a t u r e g r a d i e n t s[J].T r a n s a c t i o n s,A m e r i c a nG e o p h y s i c a lU n i o n,1957,38(2):222-232.[16]白青波,李旭,田亚护,等.冻土水热耦合方程及数值模拟研究[J].岩土工程学报,2015,37(增刊2):131-136.[17]王晓刚.冻土区桩土体系冻胀融沉特性研究[D].西安:西安科技大学,2019.[18]周家作,谭龙,韦昌富,等.土的冻结温度与过冷温度试验研究[J].岩土力学,2015,36(3):777-785. [19]赵刚,陶夏新,刘兵.原状土冻融过程中水分迁移试验研究[J].岩土工程学报,2009,31(12):1952-1957. [20]孙蕾,王磊,蔡冰,等.土壤水分测定方法简介[J].中国西部科技,2014,13(11):54-55.[21]陶文铨.传热学[M].西安:西北工业大学出版社,2006.[22] L uN,L i k o s W J.非饱和土力学[M].北京:高等教育出版社,2012:269-287.[23]徐学祖,邓友生.冻土中水分迁移的实验研究[M].北京:科学出版社,1991:21-28.[24]王刚,安琳.C OM S O L M u l t i p h y s i c s工程实践与理论仿真:多物理场数值分析技术[M].电子工业出版社,2012:167-194.[25]顾汪明,周金星,王彬,等.冻融循环作用对黑土水稳性团聚体特征的影响[J].中国水土保持科学,2020,18(4):45-52.[26]李智明.冻土水热力场耦合机理研究与工程应用[D].哈尔滨:哈尔滨工业大学,2017.[27]王冲.施用生物炭对土壤热导率的影响及机制研究[D].沈阳:沈阳农业大学,2018.[28] Z h o uJZ,L iD Q.N u m e r i c a l a n a l y s i so f c o u p l e dw a-t e r,h e a t a n ds t r e s s i ns a t u r a t e df r e e z i n g s o i l[J].C o l dR e g i o n sS c i e n c e a n dT e c h n o l o g y,2012,72:43-49. [29]H u a n g X,R u d o l p h D L.C o u p l e d m o d e l f o r w a t e r,v a p o u r,h e a t,s t r e s sa n ds t r a i n f i e l d s i nv a r i a b l y s a t u-r a t e d f r e e z i n g s o i l s[J].A d v a n c e s i n W a t e rR e s o u r c e s,2021,154:e103945.[30] L a iY M,P e iW S,Z h a n g M Y,e t a l.S t u d y o n t h e o r ym o d e l o fh y d r o-t h e r m a l-m e c h a n i c a l i n t e r a c t i o n p r o c e s si n s a t u r a t e d f r e e z i n g s i l t y s o i l[J].I n t e r n a t i o n a l J o u r n a lo fH e a t a n d M a s sT r a n s f e r,2014,78:805-819. [31]何瑞霞,金会军,赵淑萍,等.冻土导热系数研究现状及进展[J].冰川冻土,2018,40(1):116-126.[32] B a iR Q,L a iY M,Z h a n g M Y,e t a l.S t u d y o nt h ec o u p l e dh e a t-w a t e r-v a p o r-m e c h a n i c s p r o c e s s o f u n s a t u-r a t e d s o i l s[J].J o u r n a l o f H y d r o l o g y,2020,585:e124784.[33] L iK Q,L iD Q,L i uY.M e s o-s c a l e i n v e s t i g a t i o n so nt h e e f f e c t i v e t h e r m a l c o n d u c t i v i t y o fm u l t i-p h a s em a t e-r i a l su s i n g t h ef i n i t ee l e m e n tm e t h o d[J].I n t e r n a t i o n a lJ o u r n a l o f H e a t a n d M a s s T r a n s f e r,2020,151:e119383.(上接第186页)[17]张浩,吕茂奎,江军,等.侵蚀红壤区植被恢复对表层与深层土壤有机碳矿化的影响[J].水土保持学报,2016,30(1):244-249,314.[18]梁振春,李静,吴靖,等.退耕还林还草对黄土高原坡地磷素的影响[J].南方农业学报,2018,49(4):688-694.[19] L i uRS,W a n g D M.S o i l C,N,Pa n dKs t o i c h i o m e t-r y a f f e c t e db y v e g e t a t i o n r e s t o r a t i o n p a t t e r n s i n t h e a l-p i n e r e g i o no f t h eL o e s sP l a t e a u,N o r t h w e s tC h i n a[J].P L o SO n e,2020,15(11):e0241859.[20]何高迅,王越,彭淑娴,等.滇中退化山地不同植被恢复下土壤碳氮磷储量与生态化学计量特征[J].生态学报,2020,40(13):4425-4435.[21]陈晓旋,陈淑云,曾从盛,等.螃蟹对闽江河口湿地土壤碳氮磷含量及其生态化学计量学特征影响[J].环境科学学报,2018,38(3):1179-1188.[22] L i a n g JJ,C r o w t h e rT W,P i c a r d N,e ta l.P o s i t i v eb i o d i v e r s i t y-p r o d uc t i v i t y r e l a t i o n s h i p p r ed o m i n a n ti ng l o b a l f o r e s t s[J].S c i e n c e,2016,354(6309):79-134.[23]吴昊.不同类型群落物种多样性指数的比较研究[J].中南林业科技大学学报,2015,35(5):84-89. [24] S t e v e n sCJ,Q u i n t o n JN,B a i l e y AP,e t a l.T h e e f f e c t s o fm i n i m a l t i l l a g e,c o n t o u rc u l t i v a t i o na n di n-f i e l dv e g e t a t i v eb a r r i e r so ns o i le r o s i o na n d p h o s p h o r u sl o s s[J].S o i la n dT i l l a g eR e s e a r c h,2009,106(1):145-151.[25]杨桦,彭小瑜,杨淑琪,等.滇南喀斯特断陷盆地土地利用方式对土壤有机碳及其活性组分的影响[J].生态学报,2022,42(17):7105-7117.391第6期蒋瑶钰等:基于C OM S O L模型的黑土农耕地剖面冻结水热耦合分析。