CN111027236B - 岩土体的等效导热系数的细观尺度数值研究方法 - Google Patents

Abstract

本发明涉及一种岩土体的等效导热系数的细观尺度数值研究方法。本发明针对岩土体的非均质特性，从细观尺度角度针对岩土工程中最常见的材料——土体，基于有限元数值方法建立了模拟模型以此估计其等效导热系数。所建立的数值模型通过Monte Carlo算法进行模拟，得到具有空间随机性的土体结构，结合有限元计算方法与基于稳态分析的傅立叶公式，对土体的三种情况下的等效导热系数进行对比分析，并探究土体类型、孔隙率、饱和度、分散相的空间排布等因素对土体等效导热系数的影响。本发明从细观尺度出发，提出了一种可实现模拟岩土体的内部结构非均质性的方法，为分析其导热特性提供了新思路，且模拟方法简单易行，快速方便。

Description

岩土体的等效导热系数的细观尺度数值研究方法

技术领域

本发明属于多相复合材料特性研究领域，特别涉及一种岩土体的等效导热系数的细观尺度数值研究方法。

背景技术

岩土工程中很多材料都被视为是多相复合材料，例如，土通常被视为一个多相***，由土骨架、空气、水以及冰构成。岩土材料在环境和工程应用中发挥着重要作用，其热性能是岩土材料最关键的参数之一。导热系数、对流换热系数和比热容是影响土体热行为的主要因素，其中导热系数是反映传热能力的最关键参数，其在地热能源开发和地下工程建设中有着重要的应用。

随着化石能源燃料的急剧减少，及由此造成的温室效应，人们越来越倾向于开发使用环保绿色、可持续的可再生能源，诸如风能、太阳能、潮汐能、生物能、地热能等。在这些绿色能源中，地热能以其分布广、易开发、技术成熟的优点，近年来以极其迅猛势头高速发展，城市地铁、隧道等地下工程施工方兴未艾，但由于没有统一的国家与行业标准，导致测试方法各行其是，给准确计算、评估、评价导热系数带来了巨大困难。因此，准确测量导热系数对节省投资，合理开发有着重要的意义。

在地下工程建设过程中，其导热系数是计算温度场最主要的热物理参数之一，同时也是热工地面冻结技术施工设计时确定冻结帷幕的重要参数。在地下工程(如地铁隧道和深部矿井等)人工冻结法设计与施工中，准确选取未冻土和冻土的热物理参数也是非常重要的。在实际工程中对施工区域内的所有土体进行测量是不现实的，合理评价热导率对于岩土工程建设具有重要的理论意义与工程实用价值。

导热系数受多种因素共同影响，包括土体的矿物成分、粒度、级配、粒度、干密度、孔隙度、含水量、温度、孔隙大小、孔隙形状、孔隙取向和孔隙空间排列等。确定导热系数的方法主要包括试验测试和经验模型预测两种方法。

试验测量方法一般分为探针法与现场钻孔测量两类。探针法通过设备可以直接测得土体的导热系数，具有测试周期短、精度高、测试环境要求相对较低，设备便携性好等特点，但受限于探针尺寸只能测量浅层的土体的导热系数。现场钻孔试验可以测得综合导热系数，但操作繁琐费用较高，其结果只是当地很小的范围土体值。试验法测得的导热系数精度较高，但与经验公式计算方法相比，相对耗时且比较昂贵。

第二类确定方法是经验公式计算法，可以提供比较粗略的估计值。该方法只需确定公式的输入参数，即可得到等效导热系数，这有利于对各种材料性能的估计，但这种预测方法很大程度上受物理模型的影响，只能作为实际工程前进行预测的一种简单方法。

多相材料具有较强的非均质性，同时受多种因素的影响，因此对多相复合材料的导热系数的确定较为复杂。而以上方法，均不能独立探究某一单一因素对于多相复合材料的等效导热系数的影响。因此如何准确的测定多相复合材料的等效导热系数，同时探究不同因素对导热系数的影响，一直是人们关注的热点问题。

发明内容

本发明的目的在于提供一种岩土体的等效导热系数的细观尺度数值研究方法。本发明所提供的这种细观尺度数值研究方法，是基于有限元方法建模，通过Monte Carlo蒙特卡洛算法实现岩土体内部结构的非均质性，并评估孔隙率、饱和度、土体类型、土体状态等多种因素对等效导热系数的影响。

在本发明中所提出的从细观尺度基于有限元方法研究岩土体的研究方法，具有以下四个显著特点。一是该方法从细观尺度角度出发，通过有限元方法，以岩土工程中的土体为例，建立数值模拟模型。与常规的试验方法、经验公式方法相比，计算更为简单，能够较好地模拟岩土体的非均质性。第二个显著特点是该方法计算原理简单，基于稳态热分析的原理，在材料上下表面施加一定的温度差，通过有限元计算即可得到材料各节点的热流密度与温度场分布情况。将有限元节点计算结果带入傅立叶公式即可得到某一分布情况下的多相材料的等效导热系数。三是该方法借助Monte Carlo蒙特卡洛算法，实现了材料内部的异质性，每次模拟都可得到具有不同空间排布的多相土体结构，并通过Monte Carlo蒙特卡洛结果计算获取等效导热系数的统计规律。四是该方法与试验方法或传统经验公式相比，可改变某一单一因素研究其对导热系数的影响，有利于探究影响导热系数的最主要因素。总而言之，本发明所提出的岩土体的等效导热系数的细观尺度数值研究方法具有以下优点：和传统方法相比，计算原理简单明了，可操作性强，省时省力；通过有限元方法建模，可构建多相材料内部的非均质性；通过Monte Carlo蒙特卡洛模拟，可得到等效导热系数的统计规律；可实现单一变量控制，探究多种因素对整体等效导热系数的影响。

本发明提供的岩土体的等效导热系数的细观尺度数值研究方法，包括如下步骤：

(1)确定研究对象属性。包括所要研究的岩土体模拟试样尺寸、各组分的几何尺寸与材料参数、各组分的百分比等输入参数。各相的体积分数可根据孔隙率(n)和饱和度(Sr)表示，其中孔隙率为分散相占总体积的百分比，饱和度代表水相所占孔隙体积的百分比，其他相可根据以上两个参数转换得到：土骨架体积含量为1-n，孔隙相体积含量为n，水相体积含量为n·Sr，空气相为n·(1-Sr)。

(2)对预模拟的岩土体，构建数值模拟模型。其中分散相(夹杂相)为圆形，随机分布在基质内部，并借助Monte Carlo蒙特卡洛方法形成随机模拟试样。用Monte Carlo蒙特卡洛方法模拟时，首先必须生成满足随机均匀分布的随机变量。

任一分散相颗粒在投放域中的位置都是随机分布的，在笛卡尔坐标系XOY中，任一分散相颗粒的参考点O的坐标可表示为

式中：f₁₁(x，y)、f₁₂(x，y)、f₂₁(x，y)、f₂₂(x，y)分别为分散相投放区域(基质相)的边界曲线函数。

对于矩形截面，任意参考点O的坐标为：

式中：X_max、X_min、Y_max、Y_min分别为分散相投放区域边界上横、纵坐标的最大值和最小值。

考虑到模拟的复合材料的夹杂相为圆形，则有：

式中：Range(1,1)、Range(1,2)、Range(2,1)、Range(2,2)分别为模拟试样的边缘位置横、纵坐标的最小值和最大值。η₁、η₂为[0,1]随机数。

在实际构建细观模型过程中，最复杂、最难以解决的就是重合性判断。由于模拟过程中各分散相之间无相互接触，对此可以在每个分散相颗粒周围设定一个影响区域，在此范围内其他夹杂相不能进入，也即应满足投放条件，这可用相邻夹杂相颗粒中心距来判别，即

式中：x_O、y_O为本次生成的夹杂颗粒的参考点坐标；X_O、Y_O为已经生成并投放的分散颗粒的参考点坐标，并对其进行编号，i∈R；r₁、r₂分别为本次生成的分散相的半径、已经投放的夹杂相的半径。

需要注意的是，以上进行重合性判断的方法，又称“拿和放”方法，先拿一个分散相颗粒，尝试放入模型中。一旦发现该分散相颗粒与模型中已存在的某一分散相重合了，再拿出来，然后再次尝试放入模型，该过程一直持续到找到合适的分散相位置。具体操作过程可参照附图1。

根据以上基本原理的介绍，即可实现具有随机分布的分散相的岩土体模型建立。

(3)对建立的有限元模型进行网格剖分，采用的是***格剖分的方法，通过借助有限元商业软件可以自动实现数值模拟模型的网格剖分。该方法与背景网格剖分方法相比，能够自动与几何体的外形贴合，有较强的适用性，可自动完成网格剖分，不会因为网格映射等问题造成分散相颗粒形状的失真。网格剖分结果详见图3。

(4)施加边界条件。在模型上下表面施加一定的温度差，左右两边为绝热边界。边界条件示意图如图2所示。

(5)有限元计算与后处理。进行稳态热分析，通过有限元后处理模块可以得到最终模型的温度场、热流密度等分布情况，提取模型顶部各节点的热流密度值。具体的温度场分布云图以及热流密度分布图见图4和图5。

(6)计算等效导热系数，根据步骤五得到的模拟试样顶部各节点的热流密度，将其带入傅立叶公式，可以得到等效导热系数。具体表达式为：

其中，ΔQ/Δt为单位时间内通过的热流量，k表示导热系数，A为与热传递方向正交的材料面积，▽T表示温度梯度。考虑施加的边界条件与模型尺寸，上式可表示为：

其中q_i为顶部每个节点的热流密度，n为顶部节点总数，T_top和T_bottom分别为上下表面的温度值，L为沿热流传递方向的试样长度，k_eff为最终非均质材料的等效导热系数。

(7)参数敏感性分析

主要研究孔隙率、饱和度、分散相空间位置、土体类型以及土体状态对等效导热系数的影响。具体实现：

①孔隙率、饱和度，可以通过改变模型的输入参数n和Sr实现；

②分散相空间位置分布，可通过Monte Carlo模拟多次求解探究其对于等效导热系数的影响实现；

③土体类型。由于土体材料由不同的矿物成分构成，其导热性能也不同。土体类型的不同，具体表现为导热系数值的差异。因此通过改变多相材料中基质相的导热系数值，即可探究其对等效导热系数的影响；

④土体状态。一般根据水是否充满土中的孔隙，可将土体分为饱和土与非饱和土；根据是否含有冰相，可分为冻土与非冻土。不同的状态均可通过饱和度以及含冰量来确定。根据含冰量的不同，土体也可归为三种状态：自然状态土(无冰相)、部分冻结土(有水有冰，但实际工程中具体含冰量很难测定)、完全冻结土(无水相)。

总而言之，通过本发明方法与现有技术相比，具有以下显著效果：

1、本方法不同于传统的试验方法与经验公式方法，从细观尺寸角度出发基于有限元方法建立数值模型，并通过基于稳态分析的傅立叶公式可得到等效导热系数。本方法为求解岩土体等效导热系数提供了一个全新思路，方法简单，可操作性强，节省了时间和设备等成本；

2、本方法通过Monte Carlo模拟实现岩土体内部结构的非均质性，各夹杂相随机分布在基质中且互不交叠，能够较为真实模拟复合材料内部各组分的空间随机排布，为探究岩土体材料属性提供了新途径，算法简单明了，易于推广；

3、Monte Carlo蒙特卡洛算法可以对某一孔隙率或饱和度情况下的土体进行模拟，有利于通过多次计算探究其统计规律，从细观尺度分析影响多相材料导热系数的热学机制；

4、本方法同时考虑了若干因素对等效导热系数的影响，涉及孔隙率、饱和度、分散相空间排布、土体状态、土体类型等多种因素。

附图说明

图1是岩土体等效导热系数求解流程图。

从图中可以看出求解过程主要包括有限元建模、网格剖分、有限元求解、导热系数 计算以及影响因素分析等部分。

图2是边界条件示意图。

从图中可以看出，只在模型上下表面施加了一定的温差，左右边界为完全绝热状 态。

图3是有限元网格剖分图。

以部分冻结非饱和土为例，此时土为四相***，分别为土骨架、水、气、冰相。从图 中可以看出，在100mm×100mm的二维有限元模型中，分散相(水相、气相、冰相)随机分布在 基质相(土骨架)中，各分散相之间互不交叠，尺寸大小不一。

图4是温度场分布云图。

从图中可以看出，在一定的温度荷载下，整体的温度场分布比较均匀，成层分布。

图5是热流密度分布云图。

从图中可以看出，导热系数大的组分其通过的热流密度也会更大一些。

图6是针对分散相的空间排布对岩土体导热系数影响进行探究，进行三组600次Monte Carlo模拟计算结果。

其中，(a)为组A累积平均导热系数；(b)为组B累积平均导热系数；(c)为组C累积平 均导热系数；(d)为三组累积平均导热系数；(e)为三组累积标准差。

图7是针对土体类型对岩土体导热系数影响进行探究的模拟结果。

图8为不同状态下的土样有限元网格剖分图。

其中，(a)表示自然状态；(b)表示部分冻土状态；(c)表示完全冻土状态。

图9是针对孔隙率和饱和度对岩土体导热系数影响进行的探究，具体表现为不同状态下导热系数与孔隙率、饱和度的关系图。

其中，(a)表示自然状态；(b)表示部分冻土状态；(c)表示完全冻土状态。

图10为不同饱和度下的导热系数与孔隙率的关系。

其中，(a)Sr＝20％(b)Sr＝40％(c)Sr＝60％(d)Sr＝80％(e)Sr＝100％

具体实施方式

为了使本发明的目的、技术方案及优点更加清晰，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。在进行参数敏感性分析过程中，均会用到本发明提出的基于有限元法对岩土体等效导热系数进行细观尺度的研究方法，因此此处所描述的具体实例为不同影响因素对等效导热系数的影响研究。注意在进行单一变量的敏感性分析时，其他变量应保持一致。

实施例1

本实例根据本发明提出的基于有限元法对岩土体等效导热系数进行细观尺度的数值研究方法，针对分散相的空间排布对岩土体导热系数影响进行探究。具体步骤如下：

(1)确定研究对象属性。

①确定输入参数(孔隙率和饱和度)

对自然状态土体进行模拟，研究表明自然状态下土体的含水率为13％-105％，均值为29％。选取典型值：孔隙率(n)为50％，饱和度(Sr)为78.3％。各相的体积分数根据孔隙率(n)和饱和度(Sr)确定，分散相占总体积的百分比即为孔隙率50％，饱和度代表水相所占孔隙体积的百分比，其他相可根据以上两个参数转换得到：土骨架体积含量为50％，孔隙相体积含量为50％，水相体积含量为39.15％，空气相为10.85％。

②确定几何尺寸

模拟试样尺寸为100mm×100mm。各个分散相视为圆形，半径为1-5mm随机分布在基质中且不相互接触。

③确定各相材料参数

各相的导热系数分别为：土骨架的导热系数为2W/(m·K)，冰的导热系数为2.25W/(m·K)，水相的导热系数为0.5W/(m·K)，空气的导热系数为0.02W/(m·K)。

(2)对预模拟的岩土体，构建数值模拟模型。其中分散相(夹杂相)为圆形，随机分布在基质内部，并借助Monte Carlo方法形成随机模拟试样。用Monte Carlo方法模拟时，首先必须生成满足随机均匀分布的随机变量。

任一分散相颗粒在投放域中的位置都是随机分布的，在笛卡尔坐标系XOY中，对于正方形截面，任意参考点O的坐标为：

式中：X_max、X_min、Y_max、Y_min分别为分散相投放区域边界上横、纵坐标的最大值和最小值，在这里X_max-X_min＝Y_max-Y_min＝100。

考虑到模拟的复合材料的夹杂相为圆形，则有：

式中：Range(1,1)、Range(1,2)、Range(2,1)、Range(2,2)分别为模拟试样的边缘位置横、纵坐标的最小值和最大值。η₁、η₂为[0,1]随机数。

在实际构建细观模型过程中，最复杂、最难以解决的就是重合性判断。该部分可通过“拿和放”算法实现。由于模拟过程中各分散相之间无相互接触，对此可以在每个分散相颗粒周围设定一个影响区域，在此范围内其他夹杂相不能进入，也即应满足投放条件，这可用相邻夹杂相颗粒中心距来判别，即

式中：x_O、y_O为本次生成的夹杂颗粒的参考点坐标；X_O、Y_O为已经生成并投放的分散颗粒的参考点坐标，并对其进行编号，i∈R；r₁、r₂分别为本次生成的分散相的半径、已经投放的夹杂相的半径。

(3)对建立的有限元模型进行网格剖分，采用的是***格剖分的方法，通过借助有限元商业软件可以自动实现数值模拟模型的网格剖分。该方法与背景网格剖分方法相比，能够自动与几何体的外形贴合，有较强的适用性，可自动完成网格剖分，不会因为网格映射等问题造成分散相颗粒形状的失真。

(4)施加边界条件。在模型上下表面施加一定的温度差，左右两边为绝热边界。

(5)有限元计算与后处理。进行稳态热分析，通过有限元后处理模块可以得到最终模型的温度场、热流密度等分布情况，提取模型顶部各节点的热流密度值。

(6)计算等效导热系数，根据步骤五得到的模拟试样顶部各节点的热流密度，将其带入傅立叶公式，可以得到等效导热系数。具体表达式为：

其中，ΔQ/Δt为单位时间内通过的热流量，k表示导热系数，A为与热传递方向正交的材料面积，▽T表示温度梯度。考虑施加的边界条件与模型尺寸，上式可表示为：

其中q_i为顶部每个节点的热流密度，n为顶部节点总数，T_top和T_bottom分别为上下表面的温度值，L为沿热流传递方向的试样长度，k_eff为最终非均质材料的等效导热系数。

(6)重复以上步骤，进行三组600次Monte Carlo模拟，得到统计分布结果。结果表明，导热系数对于分散相的随机位置分布并不敏感。600次Monte Carlo模拟结果的累积平均导热系数介于最大值和最小值之间，如图6所示。偏移量(最大和最小值之间的差值)接近1W/(m·K)，随着模拟次数的增加，计算得到的等效导热系数逐渐趋于稳定。从图6可以看出，90％的模拟结果均位于1.05±0.05W/(m·K)的区域内。因此，分散相的空间随机排布对于岩土体等效导热系数的影响并不显著。

实施例2

本实例根据本发明提出的基于有限元法对岩土体等效导热系数进行细观尺度的数值研究方法，针对土体类型对岩土体导热系数影响进行探究。具体步骤如下：

(1)确定研究对象属性。

①确定输入参数(孔隙率和饱和度)

本实例主要研究土体类型对等效导热系数的影响，因此对孔隙率和饱和度不做具体限制，。各相的体积分数根据孔隙率(n)和饱和度(Sr)确定，其中孔隙率为分散相占总体积的百分比，饱和度代表水相所占孔隙体积的百分比，其他相可根据以上两个参数转换得到：土骨架体积含量为1-n，孔隙相体积含量为n，水相体积含量为n·Sr，空气相为n·(1-Sr)。

②确定几何尺寸

模拟试样尺寸为100mm×100mm。各个分散相视为圆形，半径为1-5mm随机分布在基质中且不相互接触。

③确定各相材料参数

研究表明矿物成分不同，土颗粒的导热系数也略有差异，土颗粒的导热系数为1.2-7.5W/(m·K)。为调研不同土骨架导热系数对整体等效导热系数的影响，各相的导热系数分别设置为：土骨架的导热系数为1、2、4、8W/(m·K)，冰的导热系数为2.25W/(m·K)，水相的导热系数为0.5W/(m·K)，空气的导热系数为0.02W/(m·K)。

(2)对预模拟的岩土体，构建数值模拟模型。其中分散相(夹杂相)为圆形，随机分布在基质内部，并借助Monte Carlo方法形成随机模拟试样。用Monte Carlo方法模拟时，首先必须生成满足随机均匀分布的随机变量。

任一分散相颗粒在投放域中的位置都是随机分布的，在笛卡尔坐标系XOY中，对于正方形截面，任意参考点O的坐标为：

式中：X_max、X_min、Y_max、Y_min分别为分散相投放区域边界上横、纵坐标的最大值和最小值，在这里X_max-X_min＝T_max-Y_min＝100。

考虑到模拟的复合材料的夹杂相为圆形，则有：

式中：Range(1,1)、Range(1,2)、Range(2,1)、Range(2,2)分别为模拟试样的边缘位置横、纵坐标的最小值和最大值。η₁、η₂为[0,1]随机数。

在实际构建细观模型过程中，最复杂、最难以解决的就是重合性判断。该部分可通过“拿和放”算法实现。由于模拟过程中各分散相之间无相互接触，对此可以在每个分散相颗粒周围设定一个影响区域，在此范围内其他夹杂相不能进入，也即应满足投放条件，这可用相邻夹杂相颗粒中心距来判别，即

式中：x_O、y_O为本次生成的夹杂颗粒的参考点坐标；X_O、Y_O为已经生成并投放的分散颗粒的参考点坐标，并对其进行编号，i∈R；r₁、r₂分别为本次生成的分散相的半径、已经投放的夹杂相的半径。

(3)对建立的有限元模型进行网格剖分，采用的是***格剖分的方法，通过借助有限元商业软件可以自动实现数值模拟模型的网格剖分。该方法与背景网格剖分方法相比，能够自动与几何体的外形贴合，有较强的适用性，可自动完成网格剖分，不会因为网格映射等问题造成分散相颗粒形状的失真。

(4)施加边界条件。在模型上下表面施加一定的温度差，左右两边为绝热边界。

(5)有限元计算与后处理。进行稳态热分析，通过有限元后处理模块可以得到最终模型的温度场、热流密度等分布情况，提取模型顶部各节点的热流密度值。

(6)计算等效导热系数，根据步骤五得到的模拟试样顶部各节点的热流密度，将其带入傅立叶公式，可以得到等效导热系数。具体表达式为：

其中，ΔQ/Δt为单位时间内通过的热流量，k表示导热系数，A为与热传递方向正交的材料面积，▽T表示温度梯度。考虑施加的边界条件与模型尺寸，上式可表示为：

其中q_i为顶部每个节点的热流密度，n为顶部节点总数，T_top和T_bottom分别为上下表面的温度值，L为沿热流传递方向的试样长度，k_eff为最终非均质材料的等效导热系数。

(6)重复以上步骤，以此改变土骨架的导热系数，分别进行600次Monte Carlo模拟，得到统计分布结果，如图7所示。结果表明，土骨架导热系数对等效导热系数有显著影响。随着孔隙度的增大，土骨架的体积分数减小，土骨架的导热系数对整体导热系数贡献降低，等效导热系数逐渐降低。同时，土骨架导热系数越高，等效导热系数对孔隙度越敏感。

实施例3

本实例根据本发明提出的基于有限元法对岩土体等效导热系数进行细观尺度的数值研究方法，针对孔隙率和饱和度对岩土体导热系数影响进行探究。具体步骤如下：

(1)确定研究对象属性。

①确定输入参数(孔隙率和饱和度)

本实例主要研究孔隙率和饱和度对等效导热系数的影响，模拟中对孔隙率不做具体限制，随机生成具有一定孔隙率的多相土体结构。饱和度分别设为0％、20％、40％、60％、80％以及100％。各相的体积分数根据孔隙率(n)和饱和度(Sr)确定，其中孔隙率为分散相占总体积的百分比，饱和度代表水相所占孔隙体积的百分比，其他相可根据以上两个参数转换得到：土骨架体积含量为1-n，孔隙相体积含量为n，水相体积含量为n·Sr，空气相为n·(1-Sr)。

②确定几何尺寸

模拟试样尺寸为100mm×100mm。各个分散相视为圆形，半径为1-5mm随机分布在基质中且不相互接触。

③确定各相材料参数

各相的导热系数分别设置为：土骨架的导热系数为2W/(m·K)，冰的导热系数为2.25W/(m·K)，水相的导热系数为0.5W/(m·K)，空气的导热系数为0.02W/(m·K)。

(2)对预模拟的岩土体，构建数值模拟模型。其中分散相(夹杂相)为圆形，随机分布在基质内部，并借助Monte Carlo方法形成随机模拟试样。用Monte Carlo方法模拟时，首先必须生成满足随机均匀分布的随机变量。

任一分散相颗粒在投放域中的位置都是随机分布的，在笛卡尔坐标系XOY中，对于正方形截面，任意参考点O的坐标为：

式中：X_max、X_min、Y_max、Y_min分别为分散相投放区域边界上横、纵坐标的最大值和最小值，在这里X_max-X_min＝Y_max-Y_min＝100。

考虑到模拟的复合材料的夹杂相为圆形，则有：

式中：Range(1,1)、Range(1,2)、Range(2,1)、Range(2,2)分别为模拟试样的边缘位置横、纵坐标的最小值和最大值。η₁、η₂为[0,1]随机数。

在实际构建细观模型过程中，最复杂、最难以解决的就是重合性判断。该部分可通过“拿和放”算法实现。由于模拟过程中各分散相之间无相互接触，对此可以在每个分散相颗粒周围设定一个影响区域，在此范围内其他夹杂相不能进入，也即应满足投放条件，这可用相邻夹杂相颗粒中心距来判别，即

式中：x_O、y_O为本次生成的夹杂颗粒的参考点坐标；X_O、Y_O为已经生成并投放的分散颗粒的参考点坐标，并对其进行编号，i∈R；r₁、r₂分别为本次生成的分散相的半径、已经投放的夹杂相的半径。

(3)对建立的有限元模型进行网格剖分，采用的是***格剖分的方法，通过借助有限元商业软件可以自动实现数值模拟模型的网格剖分。该方法与背景网格剖分方法相比，能够自动与几何体的外形贴合，有较强的适用性，可自动完成网格剖分，不会因为网格映射等问题造成分散相颗粒形状的失真。

(4)施加边界条件。在模型上下表面施加一定的温度差，左右两边为绝热边界。

(5)有限元计算与后处理。进行稳态热分析，通过有限元后处理模块可以得到最终模型的温度场、热流密度等分布情况，提取模型顶部各节点的热流密度值。

(6)计算等效导热系数，根据步骤五得到的模拟试样顶部各节点的热流密度，将其带入傅立叶公式，可以得到等效导热系数。具体表达式为：

其中，ΔQ/Δt为单位时间内通过的热流量，k表示导热系数，A为与热传递方向正交的材料面积，▽T表示温度梯度。考虑施加的边界条件与模型尺寸，上式可表示为：

其中q_i为顶部每个节点的热流密度，n为顶部节点总数，T_top和T_bottom分别为上下表面的温度值，L为沿热流传递方向的试样长度，k_eff为最终非均质材料的等效导热系数。

(6)重复以上步骤，改变不同的土体状态，分别进行600次Monte Carlo模拟。图9(a)、(b)和(c)分别为未冻土、部分冻土和完全冻土的600次蒙特卡罗模拟的计算结果。总的来说，随孔隙率的增加，除完全饱和冻土外，导热系数均呈现降低趋势；相同孔隙度下，不同饱和度的导热系数存在差异。在孔隙率较低的情况下，等效导热系数值对饱和度不敏感，随孔隙率的增大，导热系数不断减小，且相邻两个导热系数值的差值增大。结果显示三种土壤状态的导热系数各不相同。未冻土作为一种三相材料，其孔隙度和饱和度之间存在着较为简单的线性关系。土体的等效导热系数随着水分逐渐填满土壤孔隙线性增加。对于冻结状态的土体在不同饱和状态下，其等效导热系数差异较大。饱和度的变化会导致冻土等效导热系数的变化，孔隙率越大影响愈明显。完全冻结状态土体不同饱和度之间差异最大，部分冻土的等效导热系数变化介于全冻土状态与自然土体状态之间。

实施例4

本实例根据本发明提出的基于有限元法对岩土体等效导热系数进行细观尺度的数值研究方法，针对土体状态对岩土体导热系数影响进行探究。具体步骤如下：

(1)确定研究对象属性。

④确定输入参数(孔隙率和饱和度)

本实例主要研究土体状态对等效导热系数的影响，主要考虑三种土体状态，分别为：自然土体状态、部分冻结状态、完全冻结状态，其中部分冻结状态中水相和冰相的体积比为1:1。对孔隙率不做具体限制，随机生成具有一定孔隙率的多相土体结构。饱和度设为0％、20％、40％、60％、80％以及100％。各相的体积分数根据孔隙率(n)和饱和度(Sr)确定，其中孔隙率为分散相占总体积的百分比，饱和度代表水相所占孔隙体积的百分比，其他相可根据以上两个参数转换得到：土骨架体积含量为1-n，孔隙相体积含量为n，水相体积含量为n·Sr，空气相为n·(1-Sr)。

⑤确定几何尺寸

模拟试样尺寸为100mm×100mm。各个分散相视为圆形，半径为1-5mm随机分布在基质中且不相互接触。

⑥确定各相材料参数

各相的导热系数分别设置为：土骨架的导热系数为2W/(m·K)，冰的导热系数为2.25W/(m·K)，水相的导热系数为0.5W/(m·K)，空气的导热系数为0.02W/(m·K)。

(2)对预模拟的岩土体，构建数值模拟模型。其中分散相(夹杂相)为圆形，随机分布在基质内部，并借助Monte Carlo方法形成随机模拟试样。用Monte Carlo方法模拟时，首先必须生成满足随机均匀分布的随机变量。

任一分散相颗粒在投放域中的位置都是随机分布的，在笛卡尔坐标系XOY中，对于正方形截面，任意参考点O的坐标为：

式中：X_max、X_min、Y_max、Y_min分别为分散相投放区域边界上横、纵坐标的最大值和最小值，在这里X_max-X_min＝Y_max-Y_min＝100。

考虑到模拟的复合材料的夹杂相为圆形，则有：

式中：Range(1,1)、Range(1,2)、Range(2,1)、Range(2,2)分别为模拟试样的边缘位置横、纵坐标的最小值和最大值。η₁、η₂为[0,1]随机数。

在实际构建细观模型过程中，最复杂、最难以解决的就是重合性判断。该部分可通过“拿和放”算法实现。由于模拟过程中各分散相之间无相互接触，对此可以在每个分散相颗粒周围设定一个影响区域，在此范围内其他夹杂相不能进入，也即应满足投放条件，这可用相邻夹杂相颗粒中心距来判别，即

式中：x_O、y_O为本次生成的夹杂颗粒的参考点坐标；X_O、Y_O为已经生成并投放的分散颗粒的参考点坐标，并对其进行编号，i∈R；r₁、r₂分别为本次生成的分散相的半径、已经投放的夹杂相的半径。

(3)对建立的有限元模型进行网格剖分，采用的是***格剖分的方法，通过借助有限元商业软件可以自动实现数值模拟模型的网格剖分。该方法与背景网格剖分方法相比，能够自动与几何体的外形贴合，有较强的适用性，可自动完成网格剖分，不会因为网格映射等问题造成分散相颗粒形状的失真。

(4)施加边界条件。在模型上下表面施加一定的温度差，左右两边为绝热边界。

(5)有限元计算与后处理。进行稳态热分析，通过有限元后处理模块可以得到最终模型的温度场、热流密度等分布情况，提取模型顶部各节点的热流密度值。

(6)计算等效导热系数，根据步骤五得到的模拟试样顶部各节点的热流密度，将其带入傅立叶公式，可以得到等效导热系数。具体表达式为：

其中，ΔQ/Δt为单位时间内通过的热流量，k表示导热系数，A为与热传递方向正交的材料面积，▽T表示温度梯度。考虑施加的边界条件与模型尺寸，上式可表示为：

其中q_i为顶部每个节点的热流密度，n为顶部节点总数，T_top和T_bottom分别为上下表面的温度值，L为沿热流传递方向的试样长度，k_eff为最终非均质材料的等效导热系数。

(6)重复以上步骤，改变不同的土体状态，分别进行600次Monte Carlo模拟。图8为三种不同状态的有限元网格剖分图，图10为数值模拟计算结果，可以看出在某一饱和度下，全冻土的导热系数显著高于其他两种情况，自然状态的导热系数最低。从图10(a)可以看出，土体在自然状态下的导热系数随孔隙率的增加几乎呈线性下降趋势。在饱和度为20％时，含冰相(部分或完全冻土)土体的等效导热系数也有类似的趋势。随着饱和度的增加，特别是当饱和度达到60％以上，含冰状态下的导热系数呈非线性下降趋势。特别地，在饱和80％时，全冻土的导热系数仅降低0.2W/(m·K)。冰的导热系数值接近于土骨架的导热系数，比水的导热系数大一个数量级，因此当冻土中的水相完全变成冰时，导热系数随孔隙率的增加呈线性增加，此时饱和冻土可以看作是由低导热系数的气相与导热系数较大的固相组成的复合材料。

Claims (1)

1.一种岩土体的等效导热系数的细观尺度数值研究方法，其特征在于：包含如下步骤：

1)建立数值模拟模型，通过蒙特卡洛算法实现多相复合材料内部结构的空间分布随机性：

首先确定模拟试样的尺寸，并将土体视为由空气、水、土骨架以及冰构成的多相复合材料，其中土骨架作为基质相，孔隙相作为分散相随机分布在基质中；各相的体积分数根据孔隙率n和饱和度Sr表示，其中孔隙率即为分散相占总体积的百分比，饱和度代表水相所占孔隙体积的百分比，其他相可根据以上两个参数转换得到：土骨架体积含量为1-n，孔隙相体积含量为n，水相体积含量为n·Sr，空气相为n·(1-Sr)；

所述模拟试样尺寸与各分散相具体尺寸如下：模拟试样为100mm×100mm的正方形试样，各个分散相视为圆形，半径为1-5mm，随机分布在模拟试样中且不相互接触，同时冰相尺寸与水相尺寸保持一致，不考虑水变为冰时的体积变化；

所述各相的导热系数分别为：土骨架的导热系数为1.2-7.5W/(m·K)，冰的导热系数为2.25W/(m·K)，水相的导热系数为0.5W/(m·K)，空气的导热系数为0.02W/(m·K)；

2)针对建立的数值模拟模型，进行有限元网格划分：

将建立的数值模拟模型导入到有限元计算软件中，实现***格剖分；有限元软件根据几何体的外形，自动完成网格剖分；

3)基于稳态分析法，对模拟试样施加边界条件：

当通过试样的热通量达到一定条件，即模型在任意时刻的温度与时间一致时，便可测量导热系数；数值模拟过程中，在试样上下表面施加10℃温差，左右表面为完全绝热边界；

4)在给定的边界条件下，进行有限元求解计算，并将有限元计算得到的试样顶部的热流密度，带入傅立叶公式得到等效的导热系数：

经过有限元软件进行稳态分析计算后，可以得到模拟试样顶部节点的热流密度，将其带入傅立叶公式，可以得到等效导热系数；具体表达式为：

其中，ΔQ/Δt为单位时间内通过的热流量，k表示导热系数，A为与热传递方向正交的材料面积，▽T表示温度梯度；考虑施加的边界条件与模型尺寸，上式可表示为：

其中q_i为顶部每个节点的热流密度，n为顶部节点总数，T_top和T_bottom分别为上下表面的温度值，L为沿热流传递方向的试样长度，k_eff为最终非均质材料的等效导热系数；

5)基于以上步骤，考虑了分散相空间排布、土骨架导热系数不同的土体类型、饱和度、孔隙率以及三种不同土体状态对导热系数的影响；

考虑土体类型对导热系数的影响，具体表现为土骨架采用不同的导热系数进行分析，分别为1W/(m·K)、2W/(m·K)、4W/(m·K)和8W/(m·K)；

所述饱和度、孔隙率的范围分别为0-100％和0-65％；

考虑三种不同土体状态对等效导热系数的影响，三种状态分别为：自然土体状态、部分冻结状态、完全冻结状态，其中部分冻结状态中水相和冰相的体积比为1:1。

Images