基于Matlab PDE 工具箱的土体冻结温度场模拟研究

姚远宏，张泽雨，刘永成

（中国建筑科学研究院天津分院，天津市 300384）

0 引言

随着城市的发展，地铁建设在各地不断兴起，并在城市交通中扮演着日益重要的角色。地铁中的联络通道是连接2 条隧道的逃生通道和排水通道，在施工时主要以冻结法来固结土体后再开挖，因此合理地模拟土体冻结温度场和冻结壁的发展速度，在提高冻结效率、节约成本、降低能耗方面具有一定的现实意义和经济意义。

Matlab PDE tool 是一款求解偏微分方程并且能使解图形化的软件，可以解决复杂几何图形和复杂边界条件的非稳定热传导问题。模型槽土体冻结温度场是三维的热传导问题，而Matlab PDE tool 只能定义在二维空间。徐梓斌等[1]利用Matlab PDE tool解决了一个短圆柱钢锭置于加热炉中的三维瞬态非稳定热传导问题；王志良等[2]通过坐标转换将三维导热问题转化为二维，并得出4 h 时刻燃烧棒的温度分布云图和热流密度云图。本文给出运用Matlab PDE tool 进行模型槽土体冻结温度场模拟的一般步骤并分析验证模拟值的准确性。

1 Matlab PDE tool 求解方法

Matlab PDE tool 利用有限元的方法求解偏微分方程。应用Matlab PDE tool 可以解决下类抛物线方程：

式中：Ω 为平面有界区域；u 为定义在区域Ω 上的实（或复数）函数；d、c、a 和f 为常数或变量；t 为时间变量。

热传导方程的通用偏微分方程为：

式中：u 为预模拟物体内的温度；k 为导热系数；qv为热源；ρ 为密度；cp为比热容。

在求解模拟区域温度u 时，对比式（1）与式（2）中的方程系数可以看出，d=cpρ，c=k，f-au=qv。

Matlab PDE tool 定义了2 类边界条件：

（1）Dirichlet 边界条件：hu=r。

（2）Neumann 边界条件：n·(c▽u)+qu=g。

其中：n 为垂直于边界的单位矢量；h、r、q、c 和g 为常量或与u 有关的变量。热传导问题中的第一类边界条件与Dirichlet 边界条件对应；第二类和第三类边界条件与Neumann 边界条件对应。

运用Matlab PDE tool 求解模型槽土体冻结导热流程一般步骤可按图1 进行。

图1 Ma tla b P DE tool 求解模型槽土体冻结导热流程图

2 冻土模型槽热传导模拟

2.1 模型槽介绍

本次运用Matlab PDE tool 模拟模型槽中-2 ℃的土，经过120 h 后的温度场。模型槽尺寸为2.0 m×1.0 m×1.2 m，放置于10 ℃的室内，其侧面及底部边界装有厚50 mm 的保温板（硬泡沫板），并在四周铺设塑料膜防止跑水。试验土体取自天津某地铁项目中的原状土，填充模型槽后补水饱和。冻结管直径为0.045 m 且水平贯穿宽度方向，在水平方向距离模型箱壁0.8 m，在竖直方向距离模型箱壁0.6 m，两冻结管水平间距0.4 m。冷冻系统由冷冻机、循环水泵、外接盐水箱组成，采用盐水氯化钙水溶液作为低温冷媒剂循环吸收周围土体的热量。冻结时保温板温度为6.6 ℃，冻结管与周围土体冻结时接触面温度为-25 ℃。通过试验测得土的密度ρ=1860 kg/m3；比热容cp=1550 J/（kg·℃）；土体导热系数k=1.26 W/（m·℃）；模型槽土体表面与外界空气的对流换热系数为10 W/（m2·℃）。

测温线采用国标PT1000 铂热电阻测温线，安捷伦采集仪实时采集温度值。6 个测温点在同一平面内且测温点1、2、3、4 所在直线与5、6 所在直线垂直；测温点2、3 距离冻结管外壁10 cm，测温点1、4距离冻结管外壁20 cm，测温点5 距离冻结管外壁20 cm，测温点6 距离冻结管外壁30 cm。测温点布置图见图2。

图2 测温点布置图

2.2 基本假设

土体冻结温度场是一个有相变、移动边界、热源以及复杂边界的导热问题[3-5]，过程比较复杂。为了便于模拟计算，在采用Matlab PDE tool 进行模拟时作如下假设：

（1）土体为均匀连续体，且各向同性。

（2）土体的热物理学性质在所研究区域内均相同。

（3）冻结过程中相变、水分迁移和冻胀变形均不考虑。

（4）冷量仅沿垂直于冻结管的方向传导。

2.3 模拟计算

根据测温点的位置，以冻结管水平贯穿的中点做正截面，以两冻结管水平间距的中点为原点，以测点1、2、3、4 所在直线为X 轴，测点5、6 所在直线为Y 轴，建立平面直角坐标系。

根据式（1）、式（2）中的方程系数对应关系得d=cpρ=2883 000 J/（m3·℃），c=k=1.26 W/（m·℃），a=0，f=0；模型槽侧面及底部有保温板，不考虑热损失，在Neumann 边界条件中设置为q=0，g=0；模型槽土体上部与周围空气存在热交换，因此设为第三类边界条件，即Neumann 边界条件，由对应关系得q=10×3600=36000，g=10×36000=360000；冻结管与周围土体接触面的边界条件为Dirichlet 边界条件，由对应关系得r=-25，h=1。

求得120 h 时模型槽中土体温度分布图，见图3。

图3 120 h 时模型槽中土体温度分布图

2.4 模拟温度与实测温度的对比分析

将测温点2、3 的实测温度与模拟温度进行对比，其瞬态温度变化图见图4。

图4 测温点2、3 的瞬态温度变化图

由图4 可知：位于2 根冻结管中间部位的2 号测温点的温度基本小于与冻结管距离相同的3 号测温点的温度；2 个测温点的温度在刚开始时均呈快速下降趋势，这是因为冻结管周围土体首先吸收冷量，使其自身温度开始降低；随着冻结过程的进行，降温速度逐渐变慢，这是由于土体中水冻结成冰需要释放大量热量，导致土体温度有所升高，所以土体的降温速度开始变慢。在快速降温时，测温点的模拟温度与实测温度相差很大，模拟温度一般都大于实测温度；当降温速度逐渐变慢时，测温点的模拟温度与实测温度在接近，当温度达到一定的低温时，模拟温度与实测温度基本一致。

出现模拟温度大于实测温度的原因一般有以下几方面：（1）测温仪器误差；（2）土的导热系数随冻结时间逐渐增大，而模拟时假设其不变；（3）模型槽保温较差，与环境存在热交换；（4）由于土中水结冰，导致土的比热容随温度降低而减小[6-8]。

将测温点4 的实测温度与模拟温度进行对比，其每隔12 h 的瞬态温度值见表1。

表1 测温点4 每隔12 h 的瞬态温度值

由表1 可知，对于距冻结管较远的测温点4 而言，其模拟温度与实测温度之间的误差较大；随着冻结时间的增加，两者误差有所减小并最终趋于稳定。分析原因可能有：（1）模型槽保温板材料温度会随着冻结时间的增加而降低，但在边界假设时设定的温度为恒温6.6 ℃；（2）模拟假设时忽略水分迁移，而水分迁移对温度变化影响较大；（3）距离冻结管较远之处的制冷效果没有达到预设冻结温度。

3 结语

（1）通过合理的假设和坐标轴建立，可利用Matlab PDE tool 来模拟三维空间的模型槽土体冻结温度场热传导问题，操作简单并且结果可图形化。冻结管附近的模拟温度曲线与实测温度曲线变化趋势基本一致，能有效预测实际工程冻结温度场的趋势和所需时间。

（2）因模拟时不考虑土内相变潜热，使得冻结过程中，实测温度与模拟温度的误差过大。随着冻结时间的增加，土内相变潜热对温度的影响降低，使得模拟温度与实测温度的差别减小。在模拟时，根据水的状态进行分阶段模拟，可以提高模拟的准确性。