一种冻土温度场方程的解析解
2024-01-01王东源张科陈熹武小鹏
摘要： 冻土在我国分布十分广泛，冻胀现象是冻土区经常出现的问题，其导致的工程病害屡见不鲜。冻胀主要是由土体内部温度变化及水分迁移造成的，是一类极其复杂的温度、渗流及应力多场耦合问题。通过对冻土中的温度场进行分析，求解含冰量及孔隙率随温度变化条件下的冻土一般性瞬态热传导方程，获得其近似解析解并应用于冻土热传导过程的预测。通过与有限元计算所得精确解对比，验证其有效性。研究成果对冻土冻胀问题的深入研究具有一定的理论和实践意义。
关键词： 冻土; 温度场方程; 含冰量; 孔隙率; 温度变化; 解析解
中图分类号： TU44 文献标志码：A 文章编号： 1000-0844（2024）06-1373-07
DOI：10.20000/j.1000-0844.20230106001
Analytical solution for the temperature
field equation of frozen soil
WANG Dongyuan1， ZHANG Ke1， CHEN Xi2， WU Xiaopeng2
（1. China Petroleum Pipeline Engineering Corporation， Langfang 065000， Hebei， China;
2. College of Civil Engineering and Architecture， Jiaxing University， Jiaxing 314001， Jiangsu， China）
Abstract： Frozen soil is widely distributed in China. As a common problem in frozen soil areas， frost heave frequently causes engineering diseases. Frost heave， which is mainly formed by the temperature change and water migration in the soil during freezing and thawing， is an extremely complex multifield coupling problem of temperature， seepage， and stress. In this paper， the general transient heat conduction equation of frozen soil under the changes of ice content and porosity with temperature was solved by analyzing the temperature field in frozen soil， and its approximate analytical solution was obtained and applied to predict the heat conduction process of frozen soil. Compared with the exact solution obtained from finite element calculation， the effectiveness of the analytical solution was verified. The research results show certain theoretical and practical significance for the intensive study of the problem of frost heave in frozen soil.
Keywords： frozen soil; temperature field equation; ice content; porosity; temperature change; analytical solution
0 引言
我国是一个冻土大国，多年冻土主要分布在高纬度的东北大小兴安岭、长白山，以及西部的青藏高原、祁连山、天山和阿尔泰山等高山、高原地区，多年冻土区约占中国国土面积的21.5%；季节冻土主要分布在黑龙江、吉林、辽宁、内蒙古、甘肃、宁夏、青海、新疆北部等地，面积约占中国国土面积的53.5%。多年冻土和季节性冻土面积总和超过国土总面积的2/3［1］。在冻土区，土体会发生冻胀现象［2］，导致桥梁、道路、管线、水渠等基础设施频受冻害困扰。
土体冻胀通常会使得桥梁桩基、路基边坡，油气管道及水渠发生较大的冻胀变形，导致各种开裂和破坏。此外，冻胀还会引起冻胀丘、冰锥［3］，并造成冻土区铁轨大变形、基坑支护结构破坏、隧道衬砌挂冰、排水系统瘫痪等病害［2，4-5］。大量事实表明［2］，冻土区发生的土体冻胀现象会造成各类公用及民用设施破坏，影响其正常运营，严重时甚至会影响到人民的公共财产及生命安全。针对冻土冻胀问题，国内外学者开展了大量的研究工作，相关资料显示［2，6-7］，土体的冻胀是指土壤中温度降低，孔隙水凝结成冰引起的土体体积膨胀现象。受控于温度状况，土体中的水分会向冻结锋面迁移、聚集并结晶形成冰透镜，是土体冻胀量的主要来源。可见，土壤中的热状况、水分状况及其变化规律是引起冻害严重与否的主要因素。土体的冻胀问题归根结底是热质迁移问题。因此，本文主要对冻土中的温度场进行分析，探讨含冰量及孔隙率随温度变化条件下的冻土一般性瞬变热传导方程的求解问题，研究成果有利于推动冻土冻胀问题的深入研究。
在冻胀过程中，土体内温度分布通常是不均匀的，会形成温度梯度，如油气管道温度往往高于附近土体的温度，在管壁附近的土体中水分会朝着温度降低的方向迁移，引起对流放热。此外，水分在冻结锋面处由于温度降低会发生相变，凝结成冰放出热量，影响着冻结锋的移动速度。这些迁移的水分会聚集在冻结锋面附近，形成冰透镜体［2，8］，使管道附近的冻结区土体发生膨胀，而未冻结区水分的减少又使土体发生固结，改变了土体原有的孔隙结构，进而对土体的渗透及变形过程产生影响。可见，土体的冻胀过程跟温度变化及水分迁移密切相关，而水分的迁移又受到温度变化的影响。土体冻胀主要是由温度的变化引起的。根据Li等［9］、Nassar等［10］及Hansson等［11］的研究成果，土体冻胀的温度场方程可表示为：
式中：ρ为土体密度;C为比热容;T为土体温度;λ为热传导系数;L为土体潜热;ρi为土体中冻结冰的密度;θi为土体孔隙中冰的含量;
式（1）是一个带有源项的抛物型偏微分方程，可以定量地描述冻土中热传导过程以及土体孔隙中含冰量的变化规律。近年来，大量学者对其进行研究，获得了许多近似解：如Lunardini［12］给出了同时考虑土体冻融表层热传导和相变的解析式;Lunardini等［13］利用热平衡积分法提出寒带土壤冻融的Neumann问题的精确解;Klement'Ev等［14］对固结冻土融化过程进行求解，获得一维土体相变问题的自相似解析表达式;Shao等［15］利用傅里叶变换，得到一维土体热传导-对流方程的解析解;李方政等［16］利用指数积分函数，获得冻土壁扩展深度和速度与冻结时间的理论函数表达式;Kurylyk等［17］基于Lunardini解，对一维扩散-对流冻土溶解温度方程进行求解。最近，Hu等［18］开始推导冻结土壤双排及三排管道附近温度场分布公式。
上述求解大多数是基于线性化热传导方程，将含冰量θi或者冻结锋面视为时间变量的函数，采用积分、分离变量、函数展开或傅里叶变换等方法进行的。然而在实际工程中，含冰量通常为温度的函数［9］：
式中：n为土体孔隙率;Si为孔隙中含冰量，且有：
式中：a为经验系数;Tf为冰冻温度。
由式（3）～（4）可见，θi不仅与时间相关，还与空间位置相关。式（1）是一个分段非线性偏微分方程，这给理论解的寻求造成较大的困难。本文基于玻尔兹曼变换，对一类工程中常见的Dirichlet边界冻土温度场公式（1）～（4）进行求解，通过修正Chen等［19］的显式级数解，构造了在定常温Dirichlet定解条件下式（1）～（4）的高精度近似解析解，并与精确解及有限元高精度数值解对比，验证其正确性。
1 冻土初-边界条件
本文主要研究的冻土温度场初始-边界条件如下：
式中：T1gt;T0，均为已知定值。式（5）～（6）为一类定常温Dirichlet定解条件，该初-边界在冻土温度场问题中较为常见，如季节性冻土土层热传导问题、冻结凿井问题、冻土区隧道衬砌隔热层边界土体的热传导问题等［20-22］。
2 冻土孔隙率的简化
当只考虑温度对冻土孔隙率的影响时，根据Lewis等［23］的研究成果，有：
式中：KT为土体骨架的体积模量;Ks为土颗粒的体积模量;βs为温度耦合系数。对式（7）关于时间t进行积分，有：
当初始孔隙率n0为定值且土体颗粒不可压缩时，有：
将式（2）～（4）及式（9）代入式（1），有：
其中，A为：
上述即为冻土孔隙率的简化情况。
3 冻土温度场方程的求解
在式（10）中，令：
其中：
将式（12）代入式（10），则有：
其中：
则初-边界条件式（5）～（6）简化为：
其中：
上述式（14）～（17）为关于U的非线性偏微分方程。近年来，学者们对其进行了大量研究，例如，Parlance等［24］给出其积分形式为：
式中：ϕ=x/t1/2为玻尔兹曼变量。Parlance等［24］基于玻尔兹曼变换，提出了扩散系数D（U）为幂律和指数律形式的半无限边界的二阶隐式解。Witelski［25］基于摄动法，对Newman边界给出了摄动解。Omidvar等［26］采用微分变换及同伦摄动法，给出无限边界的解析解。最近，Chen等［19］通过玻尔兹曼和级数展开技术，给出初-边界条件式（16）～（17）的显式级数解：
式中：ui（i=1，2，3，…，n）为计算参数。
将式（20）代入式（19），在U=U1处逐项求导，有：
式中：j=0，1，2，…，n-1。
联立式（21）即可求解得ui。该近似解对高非线性问题具有较好的收敛性［27］，然而对于线性或近似线性问题，其收敛情况并不理想，例如式（22）中无量纲化之后的热传导问题：
此时，U=T，利用级数展开技术［27］进行计算，并与精确解对比（图1）。
当i为偶数阶时，ui和f值均为复数，上述计算失效。由图1可见，伴随近似阶数n的增加，其相对误差逐步增加，该级数解是发散的。
本文对显式解式（20）进行修正，构造玻尔兹曼变量的显式解为：
式中：ui（i=1，2，3，…，n）及c为计算参数。
将式（19）代入式（18），在求解域内进行配点得：
式中：j=0，1，2，3，…，n-1。
此外，令Lf为：
此为积分方程式（19）在全域内带权ϕ的积分，当ϕ为精确解时，有：
当ϕ为近似解时，总存在一个c=c*值，使得Lf取最小，即：
式中：min（X）为对X的最小值。
联立式（20）及式（22），即可求解ui及c。具体求解步骤为：
（1） 预先设定若干c值，分别代入式（24），求解对应的ui;
（2） 将步骤（1）中的c及ui代入式（25），计算对应的L;
（3） 比较步骤（2）中的Lf值，选择最小的Lf，其对应的c及ui即可满足式（27）中的Lf|c=c*。
4 冻土温度场数值模拟
对一维冻土温度场Dirichlet定解条件进行分析，采用第3节提出的计算方法分别对其进行计算，藉此验证本文方法的正确性。
4.1 冻土温度场线性热传导问题模拟
当Tgt;Tf时，冻土温度场方程退化为普通线性热传导方程，将其无量纲化之后即为式（22），采用式（23）～（27）计算1～3阶近似解ui和c值，如表1所列。
将表1中计算参数代入式（23），可得玻尔兹曼变量1～3阶近似解，并与精确解对比（图2）。
由图2可见，随着玻尔兹曼变量的增加，温度在降低，当ϕ=x/t1/2在4.65附近时，温度趋近于0。本节计算了1～3阶近似解，伴随着阶数增加，其计算结果是逐步收敛于精确解的。表2给出了3阶近似解与精确解计算结果及其相对误差。
在表2中，3阶近似解与精确解的相对误差从T=1 ℃开始逐步增大，在T=0.7 ℃处达到最大，为4.072 9%，而后随着温度的降低逐渐减小，在靠近0处回落到3.159 13%。由表2计算结果可知，目前的结果与精确解较为吻合。对比图1可知，引入参数c之后，对于线性问题可以改善玻尔兹曼变量显式解式（20）的收敛性，加快其收敛速度。
4.2 冻土分段非线性热传导问题模拟
对于冻土热传导问题，根据Li等［2，9］的研究成果，选择土体参数（表3）。
选择边界温度T1为20 ℃，初始温度T0为-10 ℃，初始孔隙率为0.33，则式（10）可变为：
其中，A为：
则式（28）的近似解析解为：
式中：U及U1详见式（12）～（13）和式（18）。采用式（23）～（27）计算1～3阶近似解，ui和c值如表4所列。
将表4中计算参数代入式（23），可得玻尔兹曼变量1～3阶近似解，并与有限元（FEM）解（划分100 000个单元，可视为精确解）对比，如图3所示。
由图3可见，对于分段非线性的冻土热传导问题，本文提出的玻尔兹曼变量ϕ=x/t1/2解析解是随着阶数n的增加逐步收敛于FEM解的，当n=3时，已经与FEM精确解较为吻合。表5给出了3阶近似解的相对误差。
由表5中可知，3阶近似解析解的计算结果与FEM解较为接近，其最大相对误差为4.776 9%，位于T=12 ℃位置附近，在冰冻温度分界点-0.35 ℃处的误差为1.993 6%。此外，计算10 d，20 d，40 d，60 d，80 d及100 d的温度分布如图4所示。
由图4可见，加入了计算参数c的3阶近似解析解可以用于预测不同时间的温度分布。以表3给出的计算参数可知，在100 d时，在x=0.22 m附近温度开始升高，相对于初始时刻，冰冻临界点T=-0.35 ℃向右推进了约0.055 m。目前的计算结果与FEM高精度数值解吻合较好。
5 结语
本文对含冰量及孔隙率同时变化时冻土热传导方程进行求解，主要结论如下：
（1） 基于玻尔兹曼变换，对Chen等［19］的显式解析解进行修正，构造出收敛性较好的冻土热传导方程级数解析解;
（2） 将上述解应用于线性无量纲化的热传导方程及含冰量为分段非线性函数的冻土热传导方程的求解，获得1～3阶近似解析解，并与精确解及高精度FEM解进行对比，对比结果证实，本文的求解方法可以同时求解线性及非线性热传导方程。
参考文献（References）
［1］王澄海，靳双龙，施红霞.未来50 a中国地区冻土面积分布变化［J］.冰川冻土，2014，36（1）：1-8.
WANG Chenghai，JIN Shuanglong，SHI Hongxia.Area change of the frozen ground in China in the next 50 years［J］.Journal of Glaciology and Geocryology，2014，36（1）：1-8.
［2］李智明.基于复合混合物理论的冻土多场耦合研究［D］.哈尔滨：哈尔滨工业大学，2021.
LI Zhiming.Research on coupling of multifield in frozen soil based on hybrid mixture theory［D］.Harbin：Harbin Institute of Technology，2021.
［3］王小军，曾辉辉，贾海锋，等.青藏铁路多年冻土区次生不良冻土现象的调查分析［J］.铁道工程学报，2003，20（增刊1）：1-4.
WANG Xiaojun，ZENG Huihui，JIA Haifeng，et al.Investigation and analysis of seocndary unfavourable phenomenon of permafrost in permafrost area on Qinghai-Tibet railway［J］.Journal of Railway Engineering Society，2003，20（Suppl01）：1-4.
［4］沈宇鹏，王笃礼，林园榕，等.越冬基坑水平冻胀的防治措施效果分析［J］.岩土力学，2021，42（5）：1434-1442.
SHEN Yupeng，WANG Duli，LIN Yuanrong，et al.On the effect of prevention measures against horizontal frost heave of foundation pits over winter［J］.Rock and Soil Mechanics，2021，42（5）：1434-1442.
［5］欧尔峰，黄志军.寒冷地区水渠的热力耦合数值分析［J］.兰州交通大学学报，2010，29（4）：66-71.
OU Erfeng，HUANG Zhijun.Numerical analysis of the temperature and stress fields of ditch in cold regions［J］.Journal of Lanzhou Jiaotong University，2010，29（4）：66-71.
［6］何平，程国栋，朱元林.土体冻结过程中的热质迁移研究进展［J］.冰川冻土，2001，23（1）：92-98.
HE Ping，CHENG Guodong，ZHU Yuanlin.The progress of study on heat and mass transfer in freezing soils［J］.Journal of Glaciology and Geocryology，2001，23（1）：92-98.
［7］程国栋，何平.多年冻土地区线性工程建设［J］.冰川冻土，2001，23（3）：213-217.
CHENG Guodong，HE Ping.Linearity engineering in permafrost areas［J］.Journal of Glaciolgy and Geocryology，2001，23（3）：213-217.
［8］BESKOW G.Soil freezing and frost heaving with special application to roads and railroads［J］.Cold Regions Research and Engineering，1947，26（3）：410-416.
［9］LI Z M，CHEN J，SUN K，et al.Numerical simulation and experimental validation of moisture-heat coupling for saturated frozen soils［J］.Sciences in Cold and Arid Regions，2017，9（3）：250-257.
［10］NASSAR I N，HORTON R.Simultaneous transfer of heat，water，and solute in porous media：I.theoretical development［J］.Soil Science Society of America Journal，1992，56（5）：1350.
［11］HANSSON K，SIMUNEK J，MIZOGUCHI M，et al.Water flow and heat transport in frozen soil：numerical solution and freeze-thaw applications［J］.Vadose Zone Journal，2004，3（2）：693-704.
［12］LUNARDINI V J.Heat transfer with freezing and thawing［M］.New York：Elsevier Science Publishing Company，1991.
［13］LUNARDINI V J，VAROTTA R.Approximate solution to Neumann problem for soil systems［J］.Journal of Energy Resources Technology，1981，103（1）：76-81.
［14］KLEMENT'EV A F，KLEMENT'EVA E A.Self-similar solution of the problem of consolidation and thawing of frozen soil［J］.Journal of Engineering Physics，1988，54（4）：463-466.
［15］SHAO M G，HORTON R，JAYNES D B.Analytical solution for one-dimensional heat conduction-convection equation［J］.Soil Science Society of America Journal，1998，62（1）：123-128.
［16］李方政，夏明萍.基于指数积分函数的人工冻土温度场解析研究［J］.东南大学学报（自然科学版），2004，34（4）：469-473.
LI Fangzheng，XIA Mingping.Study on analytical solution of temperature field of artificial frozen soil by exponent-integral function［J］.Journal of Southeast University （Natural Science Edition），2004，34（4）：469-473.
［17］KURYLYK B L，MCKENZIE J M，MACQUARRIE K T B，et al.Analytical solutions for benchmarking cold regions subsurface water flow and energy transport models：one-dimensional soil thaw with conduction and advection［J］.Advances in Water Resources，2014，70：172-184.
［18］HU X D，HAN L，HAN Y G.Analytical solution to temperature distribution of frozen soil wall by multi-row-piped freezing with the boundary separation method［J］.Applied Thermal Engineering，2019，149：702-711.
［19］CHEN X，HE D D，YANG G Y，et al.Approximate analytical solution for Richards' equation with finite constant water head Dirichlet boundary conditions［J］.Computational and Applied Mathematics，2021，40（7）：236.
［20］罗凯乔.季节性冻土热传导性质研究［D］.哈尔滨：东北林业大学，2016.
LUO Kaiqiao.Study on thermal conductivity property of seasonal frozen soil［D］.Harbin：Northeast Forestry University，2016.
［21］吴嘉林，辛德林，张建平.井工开采技术的创新与发展［J］.煤炭工程，2014，46（1）：4-8.
WU Jialin，XIN Delin，ZHANG Jianping.Innovation and development of underground coal mining technology［J］.Coal Engineering，2014，46（1）：4-8.
［22］夏才初，范东方，李志厚，等.隧道多年冻土段隔热层厚度解析计算结果的探讨［J］.土木工程学报，2015，48（2）：118-124.
XIA Caichu，FAN Dongfang，LI Zhihou，et al.Discussion on analytical calculation for thermal-insulation layer thickness of tunnel in permafrost area［J］.China Civil Engineering Journal，2015，48（2）：118-124.
［23］LEWIS R W，SHREFLER B A.The finite element method in the static and dynamic deformation and consolidation of porous media［M］.2nd ed.New York：John Wiley，1998.
［24］PARLANCE M B，PRASAD S N，PARLANGE J Y，et al.Extension of the Heaslet-Alksne technique to arbitrary soil water diffusivities［J］.Water Resources Research，1992，28（10）：2793-2797.
［25］WITELSKI T P.Motion of wetting fronts moving into partially pre-wet soil［J］.Advances in Water Resources，2005，28（10）：1133-1141.
［26］OMIDVAR M，BARARI A，MOMENI M，et al.New class of solutions for water infiltration problems in unsaturated soils［J］.Geomechanics and Geoengineering，2010，5（2）：127-135.
［27］CHEN X，DAI Y.An approximate analytical solution of Richards equation with finite boundary［J］.Boundary Value Problems，2017，2017（1）：167.
（本文编辑：张向红）
