冰的导热系数测量方法的研究和计算模拟

Research and Computational Simulation of Measurement Methods for Thermal Conductivity of Ice

摘要: 基于常用的通过“热传导”过程测量导热系数的稳态法，从逆向思维出发，推导和设计了适用于测量受热易融化的低温不良导热体的导热系数的“冷传导”方法，并且，通过实验和计算模拟验证了这种方法的可行性。利用COMSOL多物理场仿真软件计算和分析了待测冰样品的厚度、隔热层和待测冰样品的形状等对冰的导热系数测量的影响。

Abstract: Based on the commonly used steady-state method of measuring thermal conductivity through the “heat conduction” process, a “cold conduction” method suitable for measuring the thermal conductivity of low-temperature poor thermal conductivity materials that are easily melted by heating was derived and designed from reverse thinking. The feasibility of this method was verified through experiments and computational simulations. The influence of the thickness, insulation layer, and shape of the ice sample to be measured on the thermal conductivity measurement of ice was calculated and analyzed using COMSOL multi-physical field simulation software.

1. 引言

导热系数是表征材料导热性能的一个重要参数，它不仅是评价材料热学特征的依据，也是材料在设计应用时的一个重要依据。精确测定材料的导热系数并进一步分析其影响因素对材料的改进具有很大的现实意义 [1] ，并且，不良导热体的导热系数测量是大学物理实验中的基础实验之一 [2] ，测量导热系数的实验一般有两种方法，分别为稳态法和非稳态法。非稳态法是利用边界条件和初始条件求解热传导方程得到的。一般地，对于热的良导体采用非稳态法测量其导热系数，非稳态法包括瞬态热丝法、热线法和探针法等测量方法 [3] 。稳态法是首先借助热源在待测样品内部形成稳定的温度分布，然后对待测样品进行测量 [1] 。对于不良导热体主要采用稳态法测量其导热系数。冰是生活中一种常见的低温不良导热体，但是，因为其很容易融化，因此，测量冰的导热系数在实验中有一定难度，精确测量冰的导热系数是科研工作者十分感兴趣的一个重要项目。大量前人工作验证 [4] [5] [6] ，COMSOL多物理场仿真软件能够在宏观尺度上精确模拟一些物理实验现象。因此，一方面，本文基于稳态法分析和论证测量冰的导热系数的实验方法；另一方面，借助COMSOL多物理场仿真软件来计算模拟冰的导热系数的测量，从而为冰的导热系数的实验研究提供一定的理论参考。

2. 实验原理及实验方法

2.1. 测量低温不良导热体的导热系数的实验原理

在基于稳态法测量导热系数的实验中，首先，令上铜盘成为一个恒温热源，通过“热传导”方式，将热量从上铜盘经过待测样品向下铜盘传导，温度由上至下递减。考虑到冰这种低温不良导热体受热容易融化的特性，我们从逆向思维出发，确定了采用“冷传导”测量冰导热系数的实验方案，即上铜盘提供冷源，热量由下至上传导，温度由上至下递增。常用的“热传导”测量方法的导热系数计算公式如下：

$\lambda =-\frac{mch}{\pi D\left(T_{10}-T_{20}\right)}\times {\frac{\Delta T_2\left(自然\right)}{\Delta t}|}_{T_2=T_{20}}\times \frac{D+4H}{D+2H}$ (1)

其中，m为散热铜盘的质量；c为散热铜盘的比热容；h为待测样品厚度；R为待测样品半径；T₁₀为稳态

时上铜盘(加热铜盘)的开氏温度，T₂₀为稳态时下铜盘(散热铜盘)的开氏温度， ${\frac{\Delta T_2\left(自然\right)}{\Delta t}|}_{T_2=T_{20}}$ 为下铜盘在

空气中自然散热时下铜盘在稳态温度点T₂₀的温度变化梯度(即冷却率)，D为下铜盘的直径，H为下铜盘的厚度。

但是，在“冷传导”的测量方法中，上铜盘由输出热量变为吸收热量，待测样品仍为传热物质，下铜盘由吸收热量变为输出热量。于是，将传统的“热传导”导热系数测量实验中的加热的上盘传递高温的过程转换为制冷的上盘传递低温的过程。在“低温传递过程”中，(T₁₀-T₂₀)将由“热传导”方法中的正

值变成负值；下铜盘的自然散热变为自然吸热，即自然降温变成自然升温，温度变化梯度 ${\frac{\Delta T_2\left(自然\right)}{\Delta t}|}_{T_2=T_{20}}$ 将从负值变成正值。由于温差(T₁₀-T₂₀)和温度变化梯度 ${\frac{\Delta T_2\left(自然\right)}{\Delta t}|}_{T_2=T_{20}}$ 的正负性变化恰好相抵，因此，公

式(1)将仍然适用于“冷传导”的实验方法。

2.2. 测量低温不良导热体的导热系数的实验方法

2.2.1. 实验操作内容

在实验中，在待测样品厚度为40~50 mm、50~60 mm和60~70 mm的每个区间内，我们分别制作了两个待测冰样品，并且测量了它们的导热系数。

实验过程中，设置上铜盘的恒定温度为−10℃，环境温度控制为7℃。实验开始前，先测量冰样品的厚度及直径。实验开始后，连续测量上下铜盘的温度，并实时计算上下铜板的温差即T₁₀-T₂₀，当温差T₁₀-T₂₀的变化值持续十分钟小于等于0.5℃后，可视为系统达到稳态，记录此时上下铜盘的稳态温度T₁₀和T₂₀。

在系统达到稳态后，移除样品，将上下铜盘贴在一起，使下铜盘的温度降至(T₂₀ − 3)℃以上，此时再移除上铜盘，使下铜盘在环境中自然吸热，当下铜盘升温至(T₂₀ + 3)℃以上时，该阶段完成，记录该过程中上下铜盘的温度变化。

2.2.2. 数据处理方法

实验操作中，记录了上盘“冷传导”下盘直至达到稳态的过程与下盘自然吸热两个过程中上下铜盘的实时温度变化数据。在本实验中，前期给定的参数如表1所示。

Table 1. Given physical parameters in the experimental apparatus

表1. 实验装置中的给定物理参数

由此，式(1)可改写为

$\lambda =-7779.054\times \frac{h}{R^2\left(T_{10}-T_{20}\right)}\times {\frac{\Delta T_2\left(自然\right)}{\Delta t}|}_{T_2=T_{20}}$ (2)

通过自然吸热阶段记录的实验数据，利用Excel、Origin、Matlab、Mathematica、SPSS等 [7] [8] 数据

处理软件可以计算 ${\frac{\Delta T_2\left(自然\right)}{\Delta t}|}_{T_2=T_{20}}$ 。下面以Matlab为例，根据下铜盘在环境中自然升温的实验数据，输

入时间矩阵t和温度矩阵T，运行程序后，调用Cftool工具箱，选择Polynomial形式拟合函数，阶数选为4阶。计算得到四阶多项式拟合函数T(t)的解析式。环境温度为7℃时的自然散热拟合函数图像及解析式如图1所示。

Figure 1. Natural heat dissipation fitting function diagram at ambient temperature of 7˚C

图1. 环境温度为7℃时的自然散热拟合函数图

再利用Matlab反解出在下铜盘自然升温曲线上稳态温度T₂₀所对应的稳态时间点t₂₀，并将拟合的4

阶函数对时间t求导后代入稳态时间点t₂₀，可以得到稳态时间点t₂₀对应的斜率 ${\frac{\Delta T_2\left(自然\right)}{\Delta t}|}_{T_2=T_{20}}$ 。

3. 计算模拟方法和条件设置

在采用COMSOL多物理场仿真软件进行的计算模拟中，我们选取固体传热模型对实验过程进行模拟，模型基于传热方程即公式(1)进行设置。

$\rho C_p\frac{\partial T}{\partial t}+\rho C_{p}u\cdot \nabla T+\nabla \cdot \left(-k\nabla T\right)=Q$ (3)

其中， $\rho$ 为材料密度，C_p为材料热容，k为材料导热系数，u为速度矩阵，T为温度，Q为热源项。

计算模拟中的几何模型为一个冰样品由上下两个铜盘夹住。设置铜盘的半径为35 mm，厚度为3 mm，待测冰样品的半径为35 mm。上铜盘被恒定制冷，温度设置为恒定263.15 K (即−10℃)。外部环境温度设置为280.15 K (即7℃)。考虑到环境温度为7℃，并且，下铜盘与冰接触，故下铜盘整体的初始温度设置为277.15 K (即4℃)。分析空间为一有边界热耗散的开放空气域。铜盘边界设置为热通量，由于系统周围由空气包裹，所以对于该热通量的设置，选择对流热通量，传热系数设置为9 W/(m²∙K)。对于冰的侧面环境，采用两种设置。在第一种设置中，待测样品的侧面包裹由保温泡沫Foam构成的1 mm厚的保温隔热层。在第二种设置中，待测冰样品的侧面与空气域直接接触，没有保温隔热层。网格设置采用自由四边形网格，网格大小选择极细化。在上铜盘和下铜盘设置域探针，对于冰样本，分别设置三个域探针并分别选择最大值、最小值和平均值选项。求解时利用瞬态求解器求解，时间步长设置为0.1 s，时间长度设置为3600 s，求解结果显示为模型的三维温度以及五个探针图。同时，在计算模拟中，需要设置的其他物理参数如表2所示。

Table 2. The physical parameters of the material in the computational simulation

表2. 在计算模拟中材料的物理参数

4. 结果与讨论

4.1. 计算模拟结果

4.1.1. 冰样品厚度及稳态过程的计算模拟

本文分别计算模拟了厚度为40 mm、50 mm、60 mm和70 mm的冰样品的上下盘温差及三维温度。其中，上下盘温差变化如图2所示，50 mm的冰样品稳态时的三维温度分布如图3所示。

Figure 2. The temperature difference between the upper and lower copper plates for ice samples with the thickness of 40~70 mm

图2. 厚度为40~70 mm的冰样品的上下铜盘之间的温度差

Figure 3. Three-dimensional temperature distribution of 50 mm thick ice sample wrapped with thermal insulation layer

图3. 包裹隔热层的50 mm厚度的冰样品的三维温度分布

由图2可知，随着冰的厚度增加，下盘达到稳态的起始时间逐渐增加，稳态的起始温度也逐渐增加。由图3可知，稳态时，冰样品的内部温度梯度分布均匀，由上至下温度逐渐升高。

4.1.2. 保温隔层对冰样影响的计算模拟

在计算模拟中，设置了两种冰的侧面环境，其中无隔热层的冰样品的三维温度如图4所示。与图3相比可知，该实验方法下，隔热层的引入不会使冰的整体温度升高，同时会使冰的温度分布更加均匀。但是，在相同冰层厚度下，无隔热层的冰样品边缘温度高于内部温度。由探针数据可知，无隔热层的冰样品在“稳态”时冰的平均温度(264.61 K)高于有隔热层的冰样品在稳态时的平均温度(263.37 K)，并且，无隔热层的冰样品的温度高于0℃的时间比有隔热层的冰样品更长，这说明无隔热层会使冰的温度更高，从而增大冰融化的可能性，进而影响导热系数的测量。因此，在本实验中，我们需要在冰样品侧表面紧密包裹多层毛毡隔热层，从而减少冰样品与环境之间的热传递。

Figure 4. Three-dimensional temperature distribution of an uninsulated ice sample with a thickness of 50 mm

图4. 厚度为50 mm的无隔热层冰样品的三维温度分布

4.1.3. 冰样品形状的计算模拟

在冰侧面有隔热层包裹的前提设置条件下，本文分别对50 mm厚的等截面积的三棱柱、圆柱体和长方体三种形状的待测冰样品进计算模拟。三种冰样品的三维温度分布如图5所示。计算结果表明，三种形状的冰样品在上述实验条件下均可达到稳态。

Figure 5. Three-dimensional temperature distribution of 50 mm thick ice samples of three shapes with equal cross-sectional area

图5. 厚度为50 mm的等截面积的三种形状冰样品的三维温度分布

分别记录三种冰样品的上下表面的温度，并计算上下表面的温差即实验中上下铜盘的温差，仿真数据如图6所示。由图6可知，在截面积与厚度相同的条件下，圆柱体冰样品达到稳态的时间最短。此外，当截面积相同时，圆柱体冰样品在仪器中所占空间最小，并且，圆柱体冰样品更容易放入测量模具中以便于实验操作。

Figure 6. The temperature difference between the upper and lower surfaces of three ice samples varies with time

图6. 三种形状的冰样品的上下表面温差随时间的变化

4.2. 计算模拟与实验对比

在实验中，我们测量了43.6 mm、49.4 mm、52.1 mm、57.6 mm、63.7 mm和67.9 mm厚的冰样品，其导热系数的实验数据结果和相对误差如表3所示。并且，实验所测得的导热系数与纯冰的导热系数的相对误差最大值为6.291%，证明本实验数据在一定测量误差范围内与真实值基本一致，我们采用的实验方法和实验设备能够满足对冰样品导热系数的精确测量。

Table 3. Thermal conductivity and relative error of six groups of ice samples

表3. 六组冰样品的导热系数及相对误差

以43.6 mm厚的冰样品厚度为例，其实验所得上下铜盘温差及计算模拟所得上下铜盘温差如图7所示。

由图7可知，我们的计算模拟仿真数据与实验数据的变化趋势是一致的。在实验中，43.6 mm厚的冰样品的稳态温度为8.486℃，稳态时间为950 s。在计算模拟中，待测冰样品的稳态温度为8.550℃，稳态时间为955 s。对于另外五组冰样品，其实验和计算模拟所得稳态温度和稳态时间如表4所示。

由实验结果可知，随着冰的厚度增加，下盘达到稳态的起始时间逐渐增加，稳态的起始温度也逐渐增加，这与计算模拟所得结论一致，进而验证了本文采用的计算模拟和实验方法的可行性。

通过计算模拟与实验数据的对比可知，在实际实验中，当测量冰等低温不良导热体时，所选的厚度应在40~70 mm范围内，该范围内的样品达到稳态的时间在45分钟以内，避免因为冰样品的厚度过大而导致实验时间过长，从而节省实验时间和提高科研效率。

Figure 7. Comparison of computational simulation and experimental data for an ice sample with a thickness of 43.6 mm

图7. 厚度为43.6 mm的冰样品的计算模拟仿真和实验数据对比

Table 4. Steady-state temperature and steady-state time of ice samples with different thicknesses in experiments and computational simulations

表4. 在实验与计算模拟中不同厚度的冰样品的稳态温度和稳态时间

5. 结果与讨论

本文从逆向思维出发，推导、分析和验证了测量冰等低温不良导热体的导热系数的实验原理和实验方法，并且，利用COMSOL多物理场仿真软件计算和分析了待测冰样品的厚度、隔热层和待测冰样品的形状等对冰的导热系数测量的影响。我们发现随着冰的厚度增加，下盘达到稳态的起始时间逐渐增加，同时，稳态的起始温度也逐渐增加。在实验中待测冰样品厚度范围为40~70 mm时达到稳态时间不超过45分钟，并且提出避免因为冰样品的厚度过大而导致实验时间过长。计算模拟和实验的结论都表明，达到稳态后，冰内温度梯度分布均匀，并且为由上至下温度逐渐升高。在实验中，我们需要在冰样品侧表面紧密包裹多层毛毡隔热层，从而减少冰样品与环境之间的热传递。在截面积与厚度相同的条件下，圆柱体冰样品达到稳态的时间最短，同时，圆柱体冰样品在仪器中所占空间最小并且更容易放入测量模具中以便于实验操作。同时，本文测量得到的所有冰的导热系数均大于液态水的导热系数0.59 W/(m∙K)，这符合相同物质的固态导热系数大于液态导热系数这一普适规律。基于本文的理论推导和计算模拟工作，我们所设计和验证的实验方法可广泛应用于测量受热容易融化的低温不良导热体的导热系数。在未来的工作中，我们将自主设计和制作实验装置，研究温度和杂质等因素对冰的导热系数等物理性质的影响，从而为海冰 [9] [10] 等低温固体的应用提供具有一定价值的参考数据。

致谢

衷心地感谢沈阳航空航天大学理学院物理实验中心给予本文研究工作的支持。

基金项目

2022年度辽宁省普通高等教育本科教学改革研究项目立项优质教学资源建设与共享项目：在混合式实验教学中开展跨校修读的探索与实践。2022年沈阳航空航天大学理学院教改项目：线上和线下混合式实验教学模式的探索和建设。2023年大学生创新创业训练计划项目：基于单片机的多功能导热系数测量装置的设计与研究。

NOTES

^*通讯作者。