热固耦合下裂隙产状对导热系数影响的模拟分析*

渠成堃， 周 辉， 任振群， 程文武

(1.中国科学院岩土力学研究所 岩土力学与工程国家重点试验室， 武汉 430071；2.兖州煤业集团 济宁三号煤矿， 山东 济宁 273500)

热固耦合下裂隙产状对导热系数影响的模拟分析*

渠成堃1， 周 辉1， 任振群2， 程文武2

(1.中国科学院岩土力学研究所 岩土力学与工程国家重点试验室， 武汉 430071；2.兖州煤业集团 济宁三号煤矿， 山东 济宁 273500)

为了研究裂隙产状对岩石导热系数的影响，通过模拟立方定律建立了裂隙产状与导热系数间的关系.利用多场耦合分析软件Comsol Multiphysics建立四种不同的裂隙几何模型，施加温度边界条件后，得到了不同裂隙产状下的模型整体导热系数.结果表明：裂隙对热传导有明显的阻碍效果，增大裂隙开度，模型整体导热系数呈幂函数形式减小；增大裂隙倾角，导热系数呈线性增长规律；减小裂隙数量，增大裂隙表面粗糙度，导热系数呈现增大趋势.

热固耦合；高放废料处置；裂隙；表面粗糙度；立方定律；导热系数；多场耦合分析软件；数值模拟

近年来，关于高放废料处置过程中的热固耦合研究[1-3]成为广大学者研究的重点.考虑多场耦合作用下的岩石长期稳定性，对于高放废料处置有重大意义，然而岩石是一种含有大量裂隙的各向异性体，考虑到裂隙对岩石整体力学特性的影响，研究岩石在热固耦合下的长期力学特性绝非易事.

本文对比渗流场中的立方定律[10-12]，通过多场耦合分析软件Comsol Multiphysics建立不同裂隙产状模型，以此获得等效热力学开度与导热系数间关系.

1 数学模型的构建

1.1 整体导热系数计算

q=-kT

(1)

式中：q为热通量；k为导热系数；T为不同方向的温度梯度.

由于整体热通量q不便计算，而在模型上加入热边界条件时，通过整体的热通量与通过模型某一截面的热通量相等.因此，在计算模型总体导热系数ktot时，依据热通量守恒公式，即

(2)

式中：k1为完整岩石的导热系数；T1为模型一断面的温度梯度；Ttot为模型整体温度梯度.计算过程中模型整体的Ttot由整体温度边界条件得到，而T1可以由某一截面计算数据得到.由于在取截面时仅取岩石基质，因此可以近似认为k1为完整岩石的导热系数.

1.2 导热机理分析

目前研究普遍认为，岩体中热量传递是以热传导形式进行的.因此，在求解多场耦合中温度场时，傅里叶热传导公式被广泛应用，其表达式为

(3)

式中：ρ为岩体密度；cp为岩体比热容；u为热流的对流速度；Q为热量.等号左边前两项为热传导项，第三项为热对流项，由于岩体主要为基质，因此认为其对流速度u为0，即可以忽略对流项影响.

由于岩体内部存在裂隙，热量通过裂隙时，除以热传导形式传播外，还以热辐射形式从裂隙的上表面辐射到下表面，传递形式如图1所示.

图1 裂隙岩体热量传递形式Fig.1 Form of heat transfer in fractured rock

考虑到裂隙热辐射效果，当热量在裂隙中进行传递时，需引入热辐射方程，即

(4)

式中：n为方向向量；ε为表面发射率；σ为相关常数；Tamb为环境温度；T0为沿n方向的温度差.

2 不同裂隙产状模型的建立及计算

在Comsol中建立高10cm，宽5cm的长方形岩石试样.模型上边界施加373K(100 ℃)恒定温度，下边界温度为273K(0 ℃).选取完整花岗岩，其导热系数为3.49W/(m·K)，花岗岩表面发射率为0.85.研究中构建了四种不同形式的裂隙模型，分别为单裂隙平行板模型、单倾斜裂隙模型、多倾斜裂隙模型及粗糙单裂隙模型.

2.1 单裂隙平行板模型

在模型中间构造一条水平裂隙，裂隙开度以0.02mm梯度逐渐由0mm增大至2mm，如图2所示.分析不同开度下温度变化情况，以此研究裂隙开度对导热系数的影响.

图2 2 mm裂隙模型温度分布图Fig.2 Temperature distribution in model with 2 mm fracture

根据图2及式(2)，计算不同开度下的导热系数，得到导热系数与裂隙开度关系曲线如图3所示.

图3 导热系数与裂隙开度关系曲线Fig.3 Relation curve between thermal conductivity and fracture aperture

由图3可以看出，不同裂隙开度下模型整体导热系数k随裂隙开度b的非线性变化规律可表示为

k=1.687 5b-0.359

(5)

由此可见，导热系数随裂隙开度变化呈现出与立方定律类似的幂函数关系.原因在于裂隙开度增大，相当于在岩石中添加了一层导热性差的材料，致使导热系数明显减小.当裂隙开度增大到一定程度后，可认为岩石被裂隙分隔成上下两独立部分，此时热传导对于裂隙上表面至下表面的温度传递贡献较小，热量更多是以热辐射形式传递，因而导热系数变化趋于平缓.

2.2 单倾斜裂隙模型

在长方形模型中绘制一水平裂隙，并以15°梯度逐渐增大裂隙倾角至90°，研究一系列不同倾角下模型导热系数的变化情况，如图4所示.

图4 不同倾角裂隙模型温度分布图Fig.4 Temperature distribution in models with different fracture angles

根据计算得到不同倾角裂隙温度梯度分布，取模型某一截面，计算其温度梯度加权平均值，最后结合式(2)得到导热系数k随裂隙倾角α的变化曲线，如图5所示.

图5 导热系数与裂隙倾角关系曲线Fig.5 Relation curve between thermal conductivity and fracture angle

由图5可知，导热系数k与裂隙倾角α近似呈线性关系，相应关系可表示为

k=0.02α+1.617 3

(6)

随着裂隙倾角的增大，水平截面上裂隙占有面积逐步减小，裂隙对于热传递阻碍效果明显减弱.从模型温度分布可以看出，以45°倾角为界，当裂隙倾角大于45°时，被裂隙分隔的两部分岩石呈现出明显温度变化；而裂隙倾角为90°时，可以认为模型被裂隙分隔成左右两独立部分，相当于两块完整岩块，此时裂隙不再对热传导产生影响.

2.3 多倾斜裂隙模型

在长方形模型中将多条倾斜裂隙组合以研究导热系数随裂隙数量的变化情况，如图6所示.

图6 不同裂隙产状模型温度分布图Fig.6 Temperature distribution in models with different fracture occurrences

由图6可知，模型温度以裂隙为边界被划分为几个部分，裂隙的存在对模型整体温度分布影响较大.取一温度变化不大的水平截面，计算此截面上的温度梯度，并依据式(2)计算不同裂隙产状下模型导热系数变化，如表1所示.

表1 不同裂隙产状模型的导热系数Tab.1 Thermal conductivities in models with different fracture occurrences

结合图6以及计算的导热系数可知，随着裂隙倾角的增大，裂隙对热量传递阻碍作用逐渐减小.

结合图6和表1得到的导热系数可以明显看出，45°裂隙数量增加，岩块被裂隙分割数量增多，每块岩体上温度差异较大.通过计算结果分析，45°单裂隙时，岩块最底部稳定后温度达到296 K，而随着裂隙数量增加，岩块底部最大温度分别为292、289和287 K，表明裂隙对岩块中的温度传递有明显阻碍作用.

2.4 粗糙单裂隙模型

在长方形模型中绘制网格，并定义图7中间部分网格为裂隙，从而得到简单粗糙裂隙模型，以此研究裂隙粗糙度对模型导热系数的影响规律.

图7 不同裂隙粗糙度模型温度分布图Fig.7 Temperature distribution in models with different fracture roughness

由图7可知，粗糙裂隙两侧温度差异明显.同样取一温度变化较小的水平截面，利用式(2)计算不同粗糙度裂隙模型的导热系数变化，结果如表2所示.

分析图7中温度变化情况可知，随着裂隙粗糙度的增加，裂隙处颜色逐渐变浅，裂隙上表面与下表面温差逐渐减小.由表2可知，随着裂隙表面粗糙度的增加，岩体整体导热性明显提升.裂隙表面粗糙度增加，一定程度上增大了裂隙上表面与下表面接触率.另外，增加表面粗糙度可理解为在原平行板裂隙中填充了基质，显然岩石基质导热系数远高于空气，因此，导热性较平行板模型明显提升.

表2 不同裂隙粗糙度模型的导热系数Tab.2 Thermal conductivities in models with different fracture roughness

3 模拟结果分析

根据平行板裂隙模型所得到的计算结果，可以看出导热系数与裂隙开度存在类似于立方定律的幂函数关系.

在2 mm开度的平行板裂隙模型上表面施加80 MPa轴向压力，假设岩石基质不可压缩，受压缩仅为裂隙部分，因此在岩石单轴压缩过程中裂隙开度b可以表示为

b=2-v

(7)

式中，v为模型轴向位移.代入式(5)，则导热系数与裂隙开度关系可表示为

k=1.687 5(2-v)-0.359

(8)

将此导热系数k代入Comsol计算，并将结果与导热系数为常数情况进行对比，结果如图8所示.

图8 单温度场与热固耦合条件下温度对比

Fig.8 Comparison in temperatures between single temperature field and thermal-solid coupling condition

由图8可知，如果在计算温度变化时不考虑外部荷载影响，由于裂隙阻碍了上下表面温度传递，裂隙上表面至下表面温度会出现明显跳跃，两侧温度差异较大.如果考虑热固耦合的情况，由于施加了外部荷载，使得裂隙逐渐闭合，裂隙开度减小，因而裂隙上下表面的温差变小，裂隙岩体整体导热性增加，这更加符合现场情况.

4 结 论

本文利用Comsol Multiphysics建立了含裂隙岩样模型，旨在通过模拟立方定律建立裂隙开度与导热系数的关系，主要得出以下结论：

1) 导热系数随着平行板裂隙开度的增加近似呈幂函数形式递减.裂隙对模型自上而下的热传导有明显阻碍效果，并且裂隙开度越大，阻碍效果越明显.

2) 裂隙倾角增加，模型导热系数呈近似线性递增规律，且当裂隙倾角为90°时，可认为模型被分割为左右两独立部分，其导热系数通过计算与完整岩体一致.

3) 随着模型内裂隙数量的增加，模型整体导热系数呈明显递减趋势.通过对多裂隙模型计算，一方面对倾斜裂隙模拟结果进行了印证，另一方面也说明了裂隙对模型整体传热效果有阻碍作用.

4) 通过对几种具有不同裂隙产状模型计算可知，模型整体导热系数随裂隙粗糙度增加呈递增趋势.

5) 对比了不耦合和耦合情况下裂隙岩体表面温度分布，在仅考虑温度作用时，由于裂隙存在，裂隙上表面与下表面温度呈现较大差异.而施加轴向荷载后，裂隙在荷载作用下逐渐闭合，因而裂隙岩体的导热性明显提升，这比较接近真实情况.

(LIU Xue-wei，LIU Quan-sheng，LU Chao-bo，et al.Rock mechanics coupled thermo-mechanical(TM) fracture propagation numerical manifold method numerical simulation [J].Chinese Journal of Rock Mechanics and Engineering，2014，33(7)：1432-1441.)

(CAI Guo-qing，ZHAO Cheng-gang，TIAN Hui.Nu-merical simulation of coupled thermos-hydro-mechanical behavior for engineered barriers in high-level radio-active waste disposal [J].Chinese Journal of Geotechnical Engineering，2013，35(Sup1)：1-8.)

[3]Li X L.TIMODAZ：a successful international coope-ration project to investigate the thermal impact on the EDZ around a radioactive waste disposal in clay host rocks [J].Journal of Rock Mechanics and Geotechnical Engineering，2013，5(3)：231-242.

(XUE Luan-luan.A composite element model for coupled seepage-heat transfer of fractured rock mass [J].Rock and Soil Mechanics，2016，37(1)：263-268.)

[5]张永利，曹竹，肖晓春，等.温度作用下煤体裂隙演化规律数值模拟及声发射特性研究 [J].力学与实践，2015，37(3)：350-354.

(ZHANG Yong-li，CAO Zhu，XIAO Xiao-chun，et al.Numerical simulation and acoustic emission characte-ristics of the evolution of coal body fissures under diffe-rent temperatures [J].Mechanics in Engineering，2015，37(3)：350-354.)

[6]贾春兰，朱凯.复合温度条件下石灰岩多场耦合裂隙渗透侵蚀试验研究 [J].岩土工程学报，2015，37(7)：1307-1312.

(JIA Chun-lan，ZHU Kai.Evolution of rock fracture permeability in coupled processes with variable temperatures [J].Chinese Journal of Geotechnical Engineering，2015，37(7)：1307-1312.)

[7]许增光，杨荣，柴军瑞.考虑岩体与裂隙水流热量交换作用的温度场有限元数值分析 [J].水资源与水工程学报，2014，25(4)：69-72.

(XU Zeng-guang，YANG Rong，CHAI Jun-rui.Numerical analysis of finite element on temperature field in effect of heat exchange of rock and fissure water [J].Journal of Water Resources and Water Engineering，2014，25(4)：69-72.)

[8]张岩，李宁，于海鸣，等.温度应力对裂隙岩体强度的影响研究 [J].岩石力学与工程学报，2013，32(增刊1)：2660-2668.

(ZHANG Yan，LI Ning，YU Hai-ming，et al.Research on influence of thermal stress on fractured rock mass strength [J].Chinese Journal of Rock Mechanics and Engineering，2013，32(Sup1)：2660-2668.)

[9]于永江，张华，张春会，等.温度及应力对成型煤样渗透性的影响 [J].煤炭学报，2013，38(6)：936-941.

(YU Yong-jiang，ZHANG Hua，ZHANG Chun-hui，et al.Effects of temperature and stress on permeability of standard coal briquette specimen [J].Journal of China Coal Society，2013，38(6)：936-941.)

[10]王道远，刘刚.监测技术在深基坑开挖工程中的应用 [J].地质学刊，2010，34(1)：67-72.

(WANG Dao-yuan，LIU Gang.Application of monitoring technology in excavation engineering of deep foundation [J].Journal of Geology，2010，34(1)：67-72.)

[11]卢占国，姚军，王殿生，等.平行裂缝中立方定律修正及临界速度计算 [J].实验室研究与探索，2010，29(4)：14-16.

(LU Zhan-guo，YAO Jun，WANG Dian-sheng，et al.Correction of cubic law and calculation of critical velocity in parallel fractures [J].Laboratory Research and Exploration，2010，29(4)：14-16.)

[12]徐维生，陈兴周，李娟娟.岩体粗糙裂隙渗流研究 [J].西安理工大学学报，2010，26(2)：218-222.

(XU Wei-sheng，CHEN Xing-zhou，LI Juan-juan.Research on seepage in coarse fractures of rock mass [J].Journal of Xi’an University of Technology，2010，26(2)：218-222.)

(责任编辑：钟 媛 英文审校：尹淑英)

Simulation analysis for influence of fracture occurrence on thermal conductivity under thermal-solid coupling

QU Cheng-kun1,ZHOU Hui1,REN Zhen-qun2,CHENG Wen-wu2

(1.State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan 430071,China;2.Jining III Coal Mine,Yanzhou Coal Company,Jining 273500,China)

In order to study the influence of fracture occurrence on the thermal conductivity,the relationship between the fracture occurrence and thermal conductivity was established through stimulating the cubic law.With the multi-field coupling analysis software Comsol Multiphysics,four different fracture geometry models were established.After the thermal boundary conditions were added,the thermal conductivities of whole model under different fracture occurrences were obtained.The results show that the fracture has obvious blocking effect on the thermal conduction.When the fracture aperture increases,the thermal conductivities of whole model obviously decrease in the form of power function.With increasing the fracture angle,the thermal conductivities will linearly increase.With decreasing the fracture number and increasing the fracture surface roughness,the thermal conductivities show an increasing trend.

thermal-solid coupling;disposal of high-level nuclear waste;fracture;surface roughness;cubic law;thermal conductivity;multi-field coupling analysis software;numerical simulation

2016-03-09.

国家自然科学基金资助项目(51209085，51579093，51479193).

渠成堃(1989-)，男，山东济宁人，博士，主要从事围岩开挖扰动区和岩体多场耦合等方面的研究.

22 17∶40在中国知网优先数字出版.

http：∥www.cnki.net/kcms/detail/21.1189.T.20161222.1740.032.html

10.7688/j.issn.1000-1646.2017.02.18

TU 45

A

1000-1646(2017)02-0219-06