基于细观模拟的软土导热系数数值预测模型*

周殷康 阎长虹 谢胜华 郑 军 方 圆 刘 鹏

( ①安徽工业大学建筑工程学院 马鞍山 243002)

( ②南京大学地球科学与工程学院 南京 210023)

( ③南京地铁集团有限公司 南京 210008)

0 引 言

随着整体经济水平的稳步发展，我国的城市化进程也日益加快，然而也出现了一些负面的问题:如建设用地紧缺、交通拥堵、绿化面积不足等“城市综合症”，严重制约着城市的稳健发展( 周殷康，2019) 。城市地下空间的“绿色开发”是解决“城市综合症”的有效途径，也是保障城市可持续发展的重要举措。

围岩的导热系数一直是地下空间“绿色开发”( 安全、舒适、经济、节能) 的重要影响因素。尤其是在地下空间热环境控制方面，围岩的导热系数是确定地下空间环控系统热负荷的基础( Kajtar et al.，2015) 。以地铁工程环控系统热负荷为例，研究表明通过围岩传递的热量占系统总产热的25%～40%左右，而在寒冷地区这一比例还会增加。若地下围岩的传热比例增加10%，则地铁环控设备的投资将减少108～126×104元·站－1，同时每站用电消耗将缩减1380 kW·月－1。对土体导热系数的科学评价是合理估计地下空间环控系统热负荷的依据。若估值保守，会造成很大的资源浪费; 若估计偏高，很可能导致设备超负荷运转，甚至达不到设计要求。可见，围岩导热系数是地下空间热环境控制的基础，对围岩导热系数的科学预测和合理评价在地下空间的工程投资、运行能耗、环控效果等方面( 李晓昭等，2012) 都具有重要的意义。

软土是地下空间开发中一类常见的围岩类型，尤其是在软土广泛分布的东南沿海、长江中下游等高度城市化的地区: 软土分布广泛，且土层较厚，埋深也较浅。在该区地下空间热环境控制中，对软土导热系数科学可靠的评价预测是一个研究重点。

软土属于细粒土，细粒土导热系数的影响因素众多，如矿物组成、含水率、干密度、孔隙比、温度、粒度、微结构等。其中含水率由于极强的环境敏感性，能导致土体导热系数发生剧烈波动，因而普遍认为含水率是土体导热系数的首要影响因素( Likos，2014) ，并且提出了大量相关的预测模型( Tarnawski et al.，2000; Tang et al.，2008) 。然而，这些模型大多都是建立在对大量实测数据统计分析之上的经验拟合模型。而现有细粒土导热测试大多偏向于制样测试都比较方便的中低含水率( 肖琳等，2009) 范围内的土体:如黄土( 王铁行等，2007; 刘鹏等，2016) 和膨润土( 叶为民等，2010; 陈航等，2014) ，比较缺乏高孔隙比、高含水率、低强度软土的导热系数测试数据。而且，基于中低含水率土体实测导热系数拟合得到的经验模型通常是不能直接外推到高含水率条件下的，因而对高含水率软土导热系数的预测比较困难。除了经验模型外，还有从基本热流传输过程出发建立的等效导热系数解析预测模型。例如基于热电比拟的Wiener 模型、H－S 模型以及EMT 模型以及基于以上3 类基础模型的组合模型。这类等效解析模型概念清晰、物理意义明确。然而需要指出的是，这类等效解析模型仅能对应一些比较简单的结构，相比真实土体结构有很大差异。因此，一般不用来直接预测土体导热系数，而作为土体导热系数的上下确界使用。进一步地，Maxwell 模型最早给出了均匀流体介质中圆球形夹杂无相互作用下等效导热系数的计算公式。Gruescu et al. ( 2007) 将孔隙内流体当作椭球形夹杂，基于细观均匀化理论以及单个夹杂的基本解推导了3 类不同情况下( 如Dilute、Mori-Tanaka 和PCW) 土体等效导热系数的均匀化格式，可有效反映夹杂形状和分布以及相互作用对介质导热系数的差异影响。此外，Chen et al.( 2015) 也基于夹杂理论推导了压实膨润土的等效导热系数的预测模型，尤其是考虑了压实膨润土双重孔隙结构对水分分布的影响。需要指出的是，基于“夹杂理论”的解析模型仍然只适用于一些规则形状( 圆、椭圆以及圆柱等) 的介质，与复杂结构土体还有一段距离。可见，不论是经验模型还是解析模型，都不适合直接预测高含水率软土导热系数的预测。

细观数值模拟一直是联系微观结构和宏观特性之间的有效桥梁。通过细观数值模拟计算土体导热系数一般称为数值预测模型。数值模型不仅能反映宏观因素如含水率、孔隙比的作用，而且还可以灵活考虑细观结构( 颗粒大小、形状、定向性以及空间分布) 的影响。Wang et al. ( 2008) 联合四参数随机生长法( QSGS) 和格子波尔兹曼( LBM) 导热模拟研究了多孔介质导热系数的变化规律，该方法在预测多孔介质导热系数方面具有很大优势。对于土体，Zhou et al. ( 2019) 基于QSGS 方法和嵌套策略提出了双峰孔隙结构非饱和土的等效导热系数数值预测模型，并通过大量实测数据对比验证了该数值模型的可靠性。然而，双重孔隙结构在中低含水率下比较普遍，例如在最优含水率干侧压实时的细粒土呈现典型的双峰孔径分布。对于高含水率软土，多是单峰孔隙结构。由于结构上的差异，Zhou et al.( 2019) 的模型无法预测高含水率软土的导热系数。

综上所述，对软土导热系数的科学预测是高度城市化地区地下空间热环境控制的一大关键，而目前适用于软土导热系数预测的可靠模型比较缺乏。基于此，本文尝试提出一种适用于软土导热系数的细观预测模型。

1 软土三相数字建模

本文采用QSGS 方法作为构建软土结构和三相分布的基本算法。该方法内含4 个参数( Cd，Di，Pn，) 以控制生成对象的结构。其中，Cd控制生成对象的几何大小，在本文中以颗粒作为生长相，则Cd主要用于调控颗粒粒径使其尽可能与实际软土粒径特征一致。为了简化建模，仅将实际软土的中值粒径指标( d50) 作为建模的匹配标准，两者之间的关系如式( 1) 所示:

其中，V 为代表性体积单元REV 的面积( 二维条件) ，而N 为REV 所剖分的网格数。V/N 即为分辨率。Di代表8 个方向的生长概率( 图1) 。一般地，通过Di可控制生成对象分布的定向性。若对象定向性不明显时，通常4 个主方向( Di=1，2，3，4) 和4 个角方向上( Di=5，6，7，8) 的生长概率都不变，但Di=1，2，3，4=4Di=5，6，7，8。Pn为建模对象的体积比例，当以孔隙作为建模目标时，则Pn直接对应孔隙率。

图1 生长核的八方向生长示意图Fig. 1 Eight growing directions of each core to its neighbor ones

对于软土，采用QSGS 方法进行结构建模可分为如下5 步，基本流程如下:( 1) 首先以固相为初始生长相，在给定的孔隙网络区域根据前述生长核分布概率Cd确定固相生长核的数量以及分布。具体如下: 每个网格都分配一个［0，1］之间的随机数，当该网格内随机数小于Cd将自动转为生长核网格。一般地，固相生长核的生长概率Cd必须小于固相的体积分数。( 2) 根据各方向上的生长概率Di，生长核逐次向四周8 个方向上的相邻网格扩张。同样地，给每个位于生长方向上的相邻网格赋予一个［0，1］之间的随机数。当该网格随机数小于生长核在该方向上的Di时，该网格即转化为生长相( 固相) 。

( 3) 重复第二步( 2) 直到固相体积比例Pn与实测值吻合，此时其余网格即为颗粒间的孔隙。固相所占的区域定义为颗粒骨架( Ωs) ，而孔隙区域则以Ωp表示( 图2) 。

图2 软土两相结构图:固( 黑) 和气( 白)Fig. 2 Two-phase model of soft soils:solid( black) and gas( white)

( 4) 在孔隙区域( Ωp) ，接着以水相作为生长相( 第二生长相) 完成孔隙水的生成与分布。根据和在Ωp内任何一个网格内都赋予一个［0，1］之间的随机数。如果该孔隙网格刚好临近固相边界，且格子内随机数刚好小于，则该网格变为水相网格。若该孔隙网格处于孔隙内部且已有水相网格相邻，则当该网格内随机数小于时，该网格自动转变为水相，否则维持不变。( 5) 重复第四步( 4) 直至Ωp内孔隙水的饱和度达到软土给定饱和度为止。由于软土含水率通常较高，所以孔隙内水分较多，饱和度大多超过75%( 图3) 。至此，关于软土孔隙结构和三相分布的建模即告完成。

图3 软土三相模型图:液( 灰色)Fig. 3 Multiphase model of soft soils: liquid( gray)

2 导热模拟与数值预测模型

为了计算软土导热系数，必须对前节生成的软土数字模型进行稳态导热模拟。本文采用热流平衡法求解网格内部节点温度的离散方程，并据此确定有效导热系数。模拟中边界条件的设置如图4 所示，上下为绝热边界，左右为高低温边界。

图4 导热模拟边界条件示意图Fig. 4 Schematic view about boundary conditions

对于任意两相邻网格i 和j，从相邻格子i 流入格子j 的热流密度为:

其中，Lij和qij分别为相邻格子中心点间的距离和热流密度，格子i 中心的温度为Ti，格子j 中心的温度为Tj。λij网格i、j 间的等效导热系数，取决于两个网格各自的相态及其导热系数值。液相导热系数为λw= 0.5974 W·( m·K)－1，气相导热系数为λg=0.0243 W·( m·K)－1。固相的导热系数( λs) 取决于固相中石英矿物的体积含量ζ，其导热系数由下式决定。

式中，i，j 间的等效导热系数λij为各自导热系数( λi和λj) 的调和平均值( Huai et al.，2007) 。

在稳态导热条件下，每个网格的温度维持恒定。同时流入该网格以及从网格流出的热流密度相等，也即热流达到平衡:

根据式( 5) 可建立计算区域内各网格中心温度的离散方程组。求解采用高斯－塞德尔迭代法，前后两次迭代的截断误差eps 设置为10－9i℃。根据计算出的稳态温度场，等效导热系数由下式给出:

其中，q 为从左侧每个格子流出的热流密度; Tc为右侧低温边界温度; Th为左侧高温边界温度; L 为计算区域的宽度。

本文采用方形网格系统进行孔隙建模和二维稳态导热模拟。纵横向网格的数量相等。模拟前，对网格剖分进行了优化设计( 图5) 。横轴网格数代表剖分精度，数值越大，则精度越高。当横轴网格数低于200 时，导热系数结果波动较大。当网格数超过200 后，导热系数波动较小并趋于稳定。同时，图5还展示了计算耗时与网格数之间的关系。计算时间随着网格数的增加而增加，尤其是当网格数超过200 之后，在Inter( R) CoreTMi3－4710 运算平台下，单次计算时间已超过1 h。因此，本文中采用的建模网格数为200×200。

图5 导热系数随横向网格数关系Fig. 5 Thermal conductivity versus the numbers of horizontal grids

3 模型验证

3.1 解析解验证

本节首先采用解析解来验证数值模拟的可靠性。平行模型和垂直模型( 图6) 是两种最简单的解析模型，一直作为材料传热的上下界限( Wiener 界限) 。平行模型和垂直模型下的模拟值与解析解如表1 所示，其中，λ1=1.0 W·( m·K)－1。

由表1 可知，本文导热模拟结果与平行、垂直模型的解析解非常接近，表明本文的导热系数预测模型是准确的。但对于结构复杂的软土，还需要跟实测导热系数数据进行对比作进一步的验证。

表1 平行模型、垂直模型模拟值与解析解Table 1 Simulated and analyzed data for series and parallel models

图6 平行模型、垂直模型示意图Fig. 6 Schematic diagram of series and parallel models

3.2 实测值验证

苏通GIL 综合管廊工程穿越软土区，本文选取该工程中20 个高含水土样来验证本文的模型。由于建模中需要软土的一些结构和组成等信息: 如通过粒度确定中值粒径以计算建模参数Cd。同时，固体颗粒的导热系数λs的计算也需要土体的矿物组成等信息。因此，笔者联合激光粒度分析和X 射线衍射分析对20 个土样的粒度和矿物组成等信息进行了一系列测试和分析。粒径分布结果见图7，矿物组成等如图8 所示。

图7 土体粒径分布图Fig. 7 Particle size distribution of the soils

图8 土体矿物组成分析Fig. 8 Minerals composition analysis of soils

表2 软土孔隙比、含水率、石英含量、中值粒径以及导热系数实测表Table 2 Void ratio，water content，quartz fraction，median particle diameter and measured thermal conductivity of soils

实际上，粒度分布获取容易，激光粒度分析仪可直接给出各土样的中值粒径d50。而对于矿物组成信息则要复杂得多。由于λs的确定不仅取决于矿物组成，而且与各矿物含量有关，尤其受石英的相对含量控制。因此，需要对X 射线衍射的矿物组成结果开展进一步的分析。本文采用Topas 软件对主要矿物的相对含量进行了定量拟合分析，然后根据式( 5) 即可得到颗粒骨架λs的导热系数。最后采用稳态热流法( 刘鹏等，2016) 测试了所有软土样的导热系数值。各土样的孔隙比、含水率、石英含量以及中值粒径和导热测试等所有测试结果见表2。

通过中值粒径确定Cd，固相和水相等生长相的体积比例可由含水率和孔隙比简单换算求得。同时假设软土内土颗粒定向性较差，因此Di=1，2，3，4=4Di=5，6，7，8。由于软土属于细粒土，内部水化膜较厚，所以假定孔隙内水以“液膜”形式存在，这里取=( Wang et al.，2008) 。4 个建模参数都确定后，根据前述建模流程构建每个软土样的数字模型。然后，再根据石英含量计算固相导热系数λs，其他两相导热系数已知。据此可直接进行稳态热传导模拟，进而预测出等效导热系数。为检验数值模拟的可靠性，图9 将导热系数实测值与模拟值进行了对比。

图9 软土导热系数预测值与实测值对比图Fig. 9 Modeled data versus experimental values for the thermal conductivity of soft soils

由图9 可知，20 个原状软土导热系数测试值和模拟计算结果吻合较好，大多数模拟值都处在实测值±20%范围内波动。首先验证了本文提出的数值模型的可靠性。同时也表明该模型在高含水率软土导热系数预测方面的潜力，有一定的工程应用价值。

4 细观热流路径分析

除了能有效预测软土的等效导热系数之外，本文所提出的数值模型还有另外一个重要作用，即通过细观模拟可以直观地展示土内各处局部热流路径的分布特征。如图10 所示，其中上部为土体的三相分布图，其下为对应的细观热流分布图。

图10 土体三相分布( 上) 及其对应的热流路径( 下)Fig. 10 Multiphase distribution of soft soils( above)and its corresponding heat flow paths( below)

由图10 可知，在细观尺度上，土内各处的局部热流分布是不均匀的。同时，在代表气相的低导热空白区域，代表高热流路径的红色极小。且在一些连通性较差的高导热区域，也缺乏高通量红色热流经过。红色区域基本上都是在固－水连通性较好的区域分布。可见，土内细观热流具有优势流特征:总是优先沿着连续分布的高导热区域传导，且在低导热的孔隙气区域选择绕行。

5 讨 论

考虑细观结构特征的土体导热系数预测模型是一大难点。在众多土体导热系数预测模型中，一直是以基于实测数据和统计分析的经验预测模型处于主导地位。例如Johansen 模型( Johansen，1975) ，Cote-Konrad 模型( Cote et al.，2005) 以及Lu 模型( Lu et al.，2007) 等。对于经验模型，其参数取值通常比较困难，且大多无明确物理意义。还有一种是解析预测模型，例如Wiener 模型、H－S 模型以及EMT 模型等。这类模型适用于结构简单的介质，而土体结构非常复杂，解析模型通常仅作为实际土体导热系数的界限值，用来圈定真实土体导热系数取值的范围。土体导热系数影响因素众多，既取决于含水率、孔隙比因素，同时也受到粒度、矿物组成以及微结构等因素的影响。明显地，经验模型属于宏观统计模型，无法反映结构等因素的影响。同时，解析模型等仅能反映一些结构简单的介质，无法应用到结构复杂的土体中。对于本文提出的基于细观模拟的导热系数数值预测模型，不仅能反映传统的含水率、孔隙比等影响因素，而且还引入了粒度、矿物等因素的影响。同时，导热系数实测结果与模拟值两者比较吻合，有效验证了本文模型的可靠性。可见，本文模型相对于传统的经验模型能更加接近真实土体的结构和组成等特征。

另外，对于结构灵敏、含水率高且大多处于流塑状态下的软土导热系数的试验测试是比较困难的，而且测试耗时良久、费用昂贵。本文提出的细观预测模型，无需进行导热测试，仅需要一些必要土性信息以确定建模、模拟参数。虽然有些参数是需要试验提供，而且这些信息通常对取样要求不高。例如，即使土样破坏也不会影响矿物组成以及粒度分析的精度。

需要指出的是，本文模型目前尚不成熟，仍旧需要更多的实测数据的进一步验证。同时，该模型参数较多，而且参数获取比较复杂，会提高模型使用的门槛。虽然采用激光粒度仪测试非常方便，但X 射线衍射分析以及矿物的定量拟合比较耗时。这在一定程度上会影响模型的应用。此外，本文模型仅针对常温条件下土体导热系数的预测，对于低温冻土以及中高温下土体( 徐云山等，2017) 的导热系数的预测是不适用的。

土体由大量颗粒连结而成，实际上土体导热主要受土内颗粒间热接触控制。当颗粒间接触更加紧密，热接触条件更好，则土体导热性能越好。同时，当颗粒间接触处湿润时，也会大大提高该处热流通过的效率，导热系数迅速上升。可见，土体整体导热性能的提升实际上与内部局部的热流强化密切相关。而导热系数只是综合各处热流的宏观平均化指标，不能反映土内热流路径等细节信息。然而，目前对土内局部热流的运动及分布特征的参考数据极为匮乏。本文所提出的数值模型恰好可直观地反映真实土内热流的实际分布特征。这就为从热流分布角度深入分析土体导热系数变化奠定了基础: 通过分析和比较土内细观热流在空间分布上的细节差异，可以从源头上分析土体导热系数变化的控制机理。

6 结 论

软土导热系数的科学预测是高度城市化地区地下空间热环境控制的重要基础。已有土体导热系数预测模型在高含水率软土范围内适用性不大，基于此，本文提出了一种可考虑粒度以及矿物组成等因素的土体导热系数预测模型。并通过室内导热测试、激光粒度分析仪以及X 射线衍射分析对高含水率原状软土导热系数、粒度分布以及矿物组成的分析，对实际软土导热系数模型预测值与试验实测值进行了对比，结果比较吻合，这有效验证了本文模型的可靠性。本文主要结论如下:

( 1) 系统提出了基于细观模拟、可考虑粒度和矿物组成的土体导热系数数值预测模型。并借助激光粒度分析仪和X 射线衍射分析仪测试了20 个原状软土的粒度和矿物组成信息。代入到模型中得到其导热系数模拟值，并与实测值进行对比。结果显示模拟值大多在实测值±20%范围内，有效验证了本文模型在土体导热系数预测中的可靠性。

(2) 从热流分布角度展示了土内各处局部的热流分布情况。其结果表明，只有连续的高导热区域才是构成优势流的主要通道; 而低导热区域以及连通性较差的高导热区域对整体传热贡献都不大。