0 引言
人工冻结是一项在地基处理和土方挖掘等方面得到广泛应用的技术，冻结形成的人工冻土在地下工程支护中发挥了重要作用。随着人工冻结技术在实际工程中的应用，人工冻土的冻结过程及冻土强度等方面的研究得到了越来越多的关注。
人工冻土的研究旨在揭示冻结发展的规律，分析冻土强度的计算方法。胡向东等［1-4］对此进行了长期的研究，构建了人工冻结温度场模型，分析了巴霍尔金模型的准确性，在巴氏模型的基础上，分析了直线型单排管、双排管冻结壁温度场及平均温度的计算方法;李方政等［5］以指数积分函数的近似级数表达了人工冻土温度场的变化规律;汪仁和等［6］以数值模拟手段分析了双排冻结管作用下冻结壁温度场的形成特征;张晶等［7］和杨更社等［8］分析了深基坑开挖时冻土墙围护厚度的计算方法，并进行了数值模拟计算。在这些研究的基础上，人工冻结方法在地铁工程和地下工程［9-12］等方面得到了不同程度的应用。
然而，现有的研究工作虽然综合应用了数值模拟手段，但是并未指明对于确定的平面区域如何比选并确定最佳的冻结管布置间距，特别是在数值模拟方法的运用过程中，当更改冻结管间距以比较不同的冻结管布置方案时，需要同步更改模型及边界条件，使操作变得复杂;因此，研究一种能够准确有效估算合理冻结管布置间距的前期估算方法具有重要的意义。
本文通过研究认为，运用当量热阻及其相关原理能够有效解决这个问题。本文基于过增元的热势能理论计算了冻结过程中人工冻土的当量热阻，分析了微元体在多个冻结管影响下冻结时间的计算方法，建立了对平面区域求解冻结管最佳布置间距的数学模型。同时，介绍了针对指定工程条件的冻结范围计算方法，分析了不同冻结管布置对冻结过程和冻结需求时间的影响。
1 热势能理论及其在冻结过程分析中的应用
1．1研究假设
截取冻结管影响范围内的一个冻土柱体进行分析。人工冻土的冻结过程中，冻土表面的热传导、冻结管的热传导是2个主要的热交换途径。为简化分析边界条件，将冻结过程拆分为相互独立的2个部分，如图1所示。
图1 冻结模型
1．2冻结管影响下的冻结过程
1．2．1单个冻结管影响下的冻结过程
冻结管的影响分析应采用线热源作用下的柱体模型，冻结管处为恒定的势能场。假设冻土与冻结管之间的热传递效率远大于冻土自身的热传递效率，冻结管表面受到微小持续时间为dt的热流，热流温度为Tf，其热流向量方向为朝向轴心的柱体径向，未冻结区的温度为Tu。现以热势能原理分析线热源作用下冻结深度与冻结时间的关系。假设热量沿冻结管线热源均匀分布，取线热源上单位厚度Δl的微平面进行研究。单个冻结管影响的计算示意图如图2所示。
图2 单个冻结管影响的计算示意图
过增元等［13-14］针对热能传递的问题，提出以热势能及热势能耗散来描述物体传递热量的能力，通过导热与导电之间的状态比拟，可以认为存在描述热能在温度场中所具有的“热势能”状态参量。其定义为
式中:E为热势能，J•K;Q为热力学能，J;T为温度，K。
对比dt作用前后的初始状态与最终状态，结合热势能耗散的定义［13-14］，已冻结区V(Tf，dt)的热力学能及温度均未发生变化，故其热耗散为零，亦即当量热阻为零。
对比dt作用前后的初始状态与最终状态，已冻结区V(Tf，dt)的热势能耗散为
式中:Qu为冻结前V(Tf，dt)所具有的热力学能;Qf为冻结后V(Tf，dt)所具有的热力学能;Tu为冻结前V(Tf，dt)所具有的温度;Tf为冻结后V(Tf，dt)所具有的温度。
考虑解冻过程的能量变化，式(2)可扩展为
式中:ρ为密度，kg/m3;L为相变潜热，J/kg;Wi为含冰量质量分数;Cf为比热容，J/(kg•K);D0为冻结影响半径上的点与轴心的距离，m。
热流向量的大小由导热效率及冻土与冻土管的接触面积约束，故与导热路径的长度、导热系数、热传导时间及冻结管直径相关。
最大热流密度
式中:q•为热流密度，J/(m2•s);λ为导热系数，W/(m•K)。
故热流的热容量
式中:r为冻结管半径，m，可取冻结管直径r=0．2m;Qvh为热容量，J。
根据热势能定义，输入热流所具有的热势能
对于未冻结区和已冻结区组成的冻土系统，热势能耗散与输入热流所具有的热势能平衡。将式(6)代入热势能耗散方程并计算，可得
亦即冻结时间t与冻结影响半径D0的关系式为
同时，在冻结过程中，冻结区V(Tf，dt)的当量热阻
化简可得
1．2．2多个冻结管影响下的冻结过程
实际工程中的人工冻土同时受到多个冻结管的影响。通过反复修改数值模拟中的各项条件以筛选冻结管布置方案，不是一种实用的技术手段。本文通过运用当量热阻及热势能理论建立数学简化模型，能够快速的比较各种冻结管布置方案的合理性。
对于一个给定的单位厚度为Δl的平面区域及其内部任意分布的n个冻结管，在平面上任取一处微元体，目标微元体对第i个冻结管的当量热阻为Ｒi，其与第i个冻结管的距离为Di。
在热势能理论中，热系统参数与电系统参数相互对照。根据热势能、热流与热阻的关系及叠加定理，通过微元体的热流
微元体与冻结管的等效间距
在此基础上，可以得出在n个冻结管的共同作用下指定微元体的冻结时间方程
式(13)亦可展开为F(Di)的形式，亦即冻结时间t与平面距离Di的函数。
1．3冻结管布置及冻结时间的比较方法
对于任意冻结管，其产生向外发散的均匀势场，故可以将冻结管排布间距问题简化为以下数学模型。对于给定边界的平面及其内部任意分布的n个特殊点{Ai}，可以由平面内的任意点Bj，向距离其最近的特殊点做一条直线L，冻结时间方程F(Di)沿路径的积分∫LF(Di)即为在冻结管组合{Ai}的影响下，冻结影响范围到达点Bj时所需要的冻结时间t(Bj，{Ai})的近似值。简化数学模型如图3所示。
图3 冻结时间计算的简化数学模型
对于平面内的点，均求解对应的冻结时间t(Bj，{Ai})，相互比较并得出其中的最大值t0(Bx，{Ai})，该极值即为该种冻结管组合{Ai}作用下，给定平面内的全部点均达到冻结状态时所需要的冻结时间。其中，Bx为该极值所对应的点，也是此种冻结管布置下平面上的最不利冻结点。
在平面为不规则平面及多个冻结管情况时，上述方法能够有效地比较不同布置方案之间冻结时间的差异。
例如，对于图4所示的由Bj1(4，4)、Bj2(-4，2)、Bj3(-2，-2)、Bj4(2，-2)围成的平面，需布置N=4根冻结管。当按{A1}=(±0．1，±0．1)、{A2}=(±0．5，±0．5)、{A3}=(±1．0，±1．0)3种布置方法布置冻结管时，3种布置方法的t0(Bx，{Ai})所对应的Bx点均为B1(4，4)点。
上述3种布置方法对B1(4，4)计算所得的冻结时间分别为:t0(B1，{A1})=1．319tc、t0(B1，{A2})=1.287tc、t0(B1，{A3})=1．171tc。tc为由方程F(Di)产生的共有代数部分。
图4 平面形状及冻结管布置示意图(单位:m)
令A3(±1．0，±1．0)中的其中一根冻结管(1．0，1．0)的布置位置以直线向(4．0，4．0)连续移动。通过编程输出移动过程中坐标值与冻结时间的关系，如图5所示。图5中:横坐标为冻结管的X坐标值，X坐标值等于Y坐标值;纵坐标为冻结时间t0(Bj，{Ai})与tc的比值。
图5 冻结时间与冻结管位置的关系
从图5可以看出，曲线在冻结管移动到X=2.795处出现拐点。分析相关数据可以发现，当X＞2.795时，随着冻结管对平面内部影响的减弱，最不利冻结点Bx由(4．0，4．0)移动至平面内部的其他点。以上数据表明，布置冻结管时，应同时考虑冻结管对平面边界和平面内部的影响，采用整体冻结时间均匀的布置方案。
冻结管最佳布置方案的求解，是寻求一种特殊点分布{Ai—}，使得对于给定的平面，其计算得出的t0(Bx，{Ai—})小于其他任何特殊点分布方式{Ai}所得出的t0(Bx，{Ai})。可以看到的是，上述过程具有循环和穷举的特征。在精度允许范围内，通过编制程序可实现对不规则平面及复杂的多冻结管情况时t0(Bk，{Ai})的自动比较和t0(Bx，{Ai—})的结果筛选。例如，对于前述由B1(4，4)、B2(-4，2)、B3(-2，-2)、B4(2，-2)围成的平面区域，可以通过程序计算得出上述问题的一个近似解
2 冻结壁厚度的计算及冻结时间的数值模拟
现以土的含水量为W0=20%、冻结断面为L=30m、冻结高度H=10m的一处人工冻结砂土为计算对象，分析人工冻土中冻结管布置间距对所需冻结时间的影响。
根据宁建国等［15-16］的研究，冻结砂土的弹性模量以等效夹杂理论计算，可得到
式中:Es为土体弹性模量，MPa;Ei为冰的冲击动态弹性模量，MPa;μs为土体的泊松比;μi为冰的泊松比;cs为土的体积分数;ci为冰的体积分数。
根据文献［15-16］，冻结温度对冻土的弹性模量影响明显。例如:对于含水量W0=20%的冻土，其在冻结温度为-5℃时，E≈1300MPa;其在冻结温度为-8℃时，E≈1400MPa。
根据张晶等［7］和杨更社等［8］的研究，未考虑蠕变时，冻土墙厚度与侧向土压力等效均布值、弹性模量之间的关系，可以依据弹性薄板能量法进行计算。
式中:p为侧向土压力等效均布值，MPa;μ为冻土的泊松比;［u］为最大位移。
在计算冻土墙厚度时，砂土弹性模量保守取E=1200MPa、砂土泊松比取μ=0．3，最大位移由模拟实验等方法确定。张晶等［7］和杨更社等［8］的研究中，悬臂梁法、弹性薄板能量法及模拟实验所得的最大位移与H的比值［u］H≈0．02。将各项参数代入式(15)并计算可得，上述示例中冻土墙厚度e≥5．50m。
设冻结管双排布置，e为人工冻土有效厚度，L为人工冻土横向宽度，2Ｒa为厚度方向的圆心距，2Ｒb为宽度方向的圆心距，如图6所示。
图6 多排冻结管计算简图
对于L=30m、e=5．50m的平面区域，分别取冻结管布置间距(2Ｒa，2Ｒb)=(1．6，1．2)、(1．6，1．4)、(1．6，1．6)，并在L方向上双排连续布置。以编程方法计算各冻结管布置方案所对应的t0(Bk，{Ai})，所得到的冻结时间关系如表1所示。
表1 不同冻结管布置间距的冻结时间比例关系
通过数值模拟方法，可以更为直观地分析上述不同冻结管间距对冻结时间的影响。
对于L=30m、e=5．50m的平面区域，同样取冻结管间距(2Ｒa，2Ｒb)=(1．6，1．2)、(1．6，1．4)、(1．6，1．6)，并在L方向上双排连续布置。对以上3种情况分别进行数值模拟，以5d为分步间隔，取其近似已完全冻结时的分步作为冻结时间，得到如图7—9所示的温度场示意图。
图7 数值模拟温度场示意图(一)(70d，2Ｒa=1．6m、2Ｒb=1．2m)(单位:K)
图8 数值模拟温度场示意图(二)(80d，2Ｒa=1．6m、2Ｒb=1．4m)(单位:K)
图9 数值模拟温度场示意图(三)(95d，2Ｒa=1．6m、2Ｒb=1．6m)(单位:K)
对比图7—9可以看出:当冻结管间距2Ｒb增大时，达到相似温度场分布所需的时间逐渐增加;冻结管间距(2Ｒa，2Ｒb)=(1．6，1．2)所对应的冻结时间为70d，(2Ｒa，2Ｒb)=(1．6，1．4)所对应的冻结时间为80d，(2Ｒa，2Ｒb)=(1．6，1．6)所对应的冻结时间为95d。取tc0=15d以折算比例并进行比较，则折算后的比例关系如表2所示。
表2 不同冻结管布置间距的冻结时间比例关系(折算后)
对比表1和表2可以看出，3种冻结管布置方案的倍数比例相似。数值模拟方法具有分析准确及表现直观的优点，但是针对不同冻结管布置分别更改数值模拟模型及边界条件，增加了较大的工作量;而运用前述简化数学模型，能够快速地分析不同冻结管布置方案之间所对应的冻结时间的差异。
3 结论与讨论
1)本文以热势能及热势能耗散原理计算了人工冻土冻结过程中冻结影响范围与冻结时间的关系，得出了当量热阻的计算方程。计算结果表明，当量热阻与两点距离、温度差、冻结温度等密切相关。
2)本文提出了一种简化数学模型，能够计算多个冻结管共同作用下不规则平面的冻结需求时间，从而分析不同冻结管布置方案的差异。举例分析了平面区域内不同冻结管布置对冻结时间的影响。计算结果表明，冻结时间最短的冻结管布置方案应同时考虑对平面区域内所有点的影响，应是冻结均匀性最好的方案。
3)本文以等效夹杂理论计算弹性模量，以弹性薄板能量法计算冻土墙厚度，并分别采用计算方法和数值模拟方法分析不同冻结管布置对冻结时间的影响。计算结果表明，冻结管间距对冻结时间有明显影响。
4)数值模拟方法在确定温度场分布时，具有精确度高、计算速度快的优点;但是，数值模拟方法在冻结管布置间距发生变化时，需要同步更改模型及边界条件，带来了额外的工作量，因此，寻求一种准确有效的前期方法估算冻结管布置间距，再以数值模拟方法进行确认，就具有重要的工程意义。在此基础上，如何将上述2个步骤以编程等手段整合在同一个可人机交互操作的系统框架内，将是下一步的研究方向。
摘自：隧道建设