在建筑结构节能降耗分析，高温或火灾下混凝土结构的破坏分析以及大体积混凝土、核反应防护设施等的温度应力计算中，混凝土的导热系数是表征其导热能力、确定温度场的关键参数[1―4]。多年来，国内外研究者对其进行了大量研究。然而，与混凝土的力学性能[5]、渗透[6]及扩散性能[7]等的研究类似，现有研究多集中于无应力状态下混凝土结构[8]。实际上，几乎所有混凝土在其服役期内会受到力学荷载、化学侵蚀和温度变化等多种因素的共同作用[7,9―10]。其中，力学荷载会改变混凝土内部孔隙结构，产生微/宏观裂纹，引起混凝土损伤劣化，从而导致包括导热系数在内的一系列材料参数发生变化。换言之，忽略荷载作用得到的混凝土导热性能研究结果具有一定的片面性。因此，研究荷载作用下混凝土中的温度传导行为具有重要的科学意义。
近年来，国内外研究者从试验研究、理论分析和数值模拟角度陆续开展了荷载作用下混凝土导热行为的研究。试验研究方面，C?erny等[11]研究了压缩应力对砂浆导热系数的影响，发现压缩应力为混凝土抗压强度 90%的试件的导热系数高于无应力试件。这一现象可归因于温度和压力梯度导致的对流效应。然而，Perkowski[12]、Vejmelková等[13]和Zhang等[1]的试验结果均表明，力学荷载作用会降低混凝土的导热能力。
理论解析方面，Hasselman[14]考虑到温度膨胀引起的微裂纹多为“硬币形”(penny-shaped)，以裂纹密度、尺寸和方向为参数，提出了开裂固体材料导热系数的公式。Bary[15]对材料微观结构的简化描述，基于显式有效介质理论预测了各向同性开裂非饱和水泥浆的导热性能。Zhang等[1]在试验研究基础上，假定荷载作用引起的损伤各向同性，提出了预测损伤混凝土导热系数的细观模型。
数值模拟方面，Tang等[2]以Weibull分布描述混凝土的非均质性，建立二维模型，研究了混凝土细/宏观结构对砂浆和混凝土有效导热系数的影响，指出随着损伤累积，混凝土导热系数减小。Shen等[16]采用等效塑性应变表征裂纹，采用串联模型获得开裂区域的有效导热系数，研究了拉伸开裂混凝土的宏观导热性能。
上述研究促进了对于力学荷载作用下混凝土中温度传导行为的认识，然而仍存在一定不足，如：由于试验设备及条件等的限制，试验研究还较少，并且试验结果离散性较大，尚未取得比较一致的定量结果；理论分析则对混凝土的微/细观结构作了较多简化，不能反映其内部真实结构间的相互作用；数值方面，力学-热学耦合研究目前仍处于起步阶段。因而，需要对力学荷载作用下混凝土的导热行为进行更加深入的研究，以期揭示热-力耦合作用的微/细观物理机制。
鉴于此，本文从细观角度出发，结合混凝土非均质性的影响，将其视为由骨料、砂浆和二者之间的界面过渡区(ITZ)组成的三相复合材料，建立混凝土热-力耦合相互作用研究的细观数值分析模型。力学荷载作用后混凝土热传导行为研究的思想为：首先对混凝土试件进行力学作用分析，获得其内部的混凝土开裂或损伤分布；将力学模拟结果作为热传导行为模拟的初始条件，结合复合材料均匀化理论得到开裂或损伤单元的有效导热系数，进而数值获得开裂或损伤后混凝土的宏观导热特性。
1 基本理论
1.1 损伤单元的导热系数
在服役过程中，由于力学荷载等的作用，混凝土结构内部不可避免地会产生一定程度的损伤。对于各向同性材料中的一个细观单元，采用损伤因子d来表征其损伤程度，变量d取值为0和1分别表示材料完好和完全损伤。相应地，材料完好(d=0)时，其导热性能与原材料相同；而材料完全损伤(d=1)时，其导热性能与空气相同(假定损伤开裂的区域由空气填充)。当损伤因子d的值介于0～1时，可将该单元视为由完好材料基质相和裂纹(损伤)相(包含微/宏观裂纹)组成的两相复合材料。此时，可认为该单元的表观导热性能与损伤因子d直接相关。问题的关键转变为确定单元的表观导热系数ka与损伤因子d的定量关系。实际上，可以将损伤因子d等效为单元中损伤区域(裂纹)的体积分数，从而可以运用复合材料均匀化方法获得该单元的表观物理性能[17]。如图1所示，忽略损伤区域的尺寸及分布影响，将损伤区域与完好材料分别集中连续均匀的裂纹相和基质相。其中，单元的总体面积为A，裂纹相的体积分数为d，相应地基质相体积分数为(1-d)。显然，当裂纹相与基质相采用不同的排列方式时，将获得不同的表观物理性能。
图1 复合材料均匀化模型示意图
Fig.1 Schematic of homogenization models of composite materials
如图1(a)所示，当裂纹相与损伤相与温度传导的方向平行时，根据复合材料均匀化理论，损伤单元当前状态表观导热系数ka可由并联模型获得[18]：
其中，kc和km分别表示裂纹相与基质相的导热系数。
当裂纹相与损伤相与温度传导的方向平行时(如图 1(b))，损伤单元的表观系数ka由串联模型得到[18]：
1.2 分析步骤
在外荷载作用下，由于孔隙变形和微/宏细观裂纹的产生，混凝土的细观结构将发生变化，在本文中表现为各单元的损伤因子d将随着荷载水平的变化而不断改变，进而导致各单元的导热系数随之不断发生变化，最终影响混凝土试件的宏观整体有效导热性能。因而，本文首先对混凝土在外荷载作用下的力学行为进行模拟分析；然后参照?avija等[19]和杜修力等[9]的处理方法，将荷载作用后的“结果输出”作为混凝土导热行为分析的“初始输入”条件。这样处理实际上是一种单向耦合行为，即假定混凝土力学行为的退化影响温度传导行为，但是温度传导行为不会影响其力学行为。
具体分析步骤如下，计算流程见图2：
1) 在荷载作用(F)下，计算细观混凝土内任意单元的损伤状态，得到其损伤因子d；
2) 根据各单元的初始导热系数(如砂浆、界面和裂纹三者导热系数分别表示为kMor、kITZ、kc)和损伤因子d，采用式(1)或式(2)计算各单元当前表观导热系数ka；
3) 基于各单元当前导热性能，已知试件的初始及边界条件，模拟分析温度在混凝土中的传导行为。
2 力学及热学耦合计算模型
随机骨料模型是一种描述混凝土细观结构的典型的唯象模型，可表征骨料颗粒在空间的随机分布情况，概念比较清晰，常被用来研究混凝土的细观断裂破坏行为、宏观力学性能及扩散性能等[20―23]。本文首先应用该模型获得混凝土在荷载作用下的响应，进而模拟分析温度在混凝土中的传导行为。
图2 压缩荷载作用下混凝土中温度传导行为计算流程图
Fig.2 Procedure for the calculation of thermal conduction in concrete subjected to compressive loadings
2.1 混凝土细观结构及组分本构关系
如图3所示，方便起见，假定骨料为圆形[9,20]，采用Monte Carlo方法[20,23]确定每颗骨料的位置，建立二维混凝土细观随机骨料结构。图3中，混凝土试件尺寸为150 mm×150 mm，骨料体积分数约为47%，其中小石(直径12 mm)数目为56颗，中石(直径6 mm)数目为6颗。为了节省计算量，界面区的厚度取为1 mm[9,19]。在图中，圆形颗粒表示骨料，其周围环形区域表示界面区，其他区域表示砂浆，不同材料拥有不同的力学性质。当试件承受压缩荷载后，由于损伤分布不同，试件导热性能变得更加不均匀。
图3 混凝土细观尺度分析模型及网格划分
Fig.3 Mesoscopic analysis model of concrete and meshing
对于普通混凝土，由于骨料的抗拉/抗压强度远大于砂浆和界面区，因而假定骨料力学行为为弹性。对于砂浆和界面区，采用由Lubliner等[24]提出并由Lee等[25]改进的混凝土塑性损伤本构关系模型来描述其力学行为。各细观组分的力学参数见表1。基于这些力学参数对图3所示的混凝土试件进行单轴压缩破坏行为数值研究，获得混凝土的单轴抗压强度为35.0 MPa，这与Zhang等[1]的试验数据十分接近(Zhang等的试验结果将被用来验证本文细观方法的有效性)。
表1 混凝土细观组分力学参数
Table 1 Mechanical parameters for the meso-constituents of concrete
混凝土各细观组分的初始热学参数如表 2。由于界面区的物理性能参数很难直接测得，这里根据试验测得的混凝土宏观导热系数对其进行反演试算。基于以上参数，可以获得完好混凝土试件的有效导热系数为 1.74 W/(m·K)，与 Zhang等[1]的实测值一致。对于裂纹相，由于裂纹一般会充满空气，因而其热学参数取值同空气[16]。
表2 完好混凝土细观组分热学参数
Table 2 Thermo-physical parameters for the sound meso-constituents of sound concrete
注：标有“*”的数据取自文献[1]，标有“**”取自文献[26]。
2.2 力学及热学加载及边界条件
分析单轴压缩力学响应时的边界条件为：试件顶面为荷载输入边界，采用位移加载控制；底边采用法向固定约束，底边中点采用水平向及竖向同时约束；两侧为自由边界。
计算混凝土的宏观有效导热系数和温度场分布时，试件内部的初始温度均设为20℃(即室温)。前者边界条件：顶面和底面分别为100℃和60℃恒温边界，两侧为绝热边界；后者边界条件：顶面施加ISO834升温曲线[27]，其余边界为绝热边界。
2.3 有限元网格及计算
采用有限元法对荷载作用下非均质混凝土中的温度传导行为进行研究，采用四节点单元对混凝土试件进行网格划分，平均网格尺寸为 1 mm。单元形式在力学行为模拟时，为线性等参数单元；在温度传导计算时，为线性热传导单元。
3 数值计算结果及分析
3.1 模拟结果与试验结果对比
为了验证本文力学-热学单向耦合细观分析方法的合理性，这里以单轴压缩加载为例，获得了荷载作用后混凝土的导热性能。当荷载水平(试件当前应变与混凝土达到峰值应力时对应应变之比)σr=0.7时，模拟所得的试件的宏观有效导热系数与Zhang等[1]的试验数据的对比如图4所示。为了考察该方法是否具有网格敏感性，同时给出了平均网格尺寸为0.5 mm、1 mm、2 mm和3 mm时的计算结果。在Zhang等[1]的试验中，粗骨料类型为石灰石，水灰比为0.5。养护28 d后，混凝土的平均抗压强度与弹性模量分别为34.77 MPa和36.2 GPa。由图4可知，在四种网格尺寸下，混凝土在平行与垂直荷载方向上的有效导热系数均变化很小，可以认为本文细观方法对网格不敏感。对比模拟结果与试验数据可知，试验测得的混凝土在两个方向上的有效导热系数几乎全部介于基于串联和并联方法所获得的模拟结果之间，且更接近于串联模型的结果。然而，由于缺乏相关试验数据，对于其他情形仍需深入研究。尽管如此，仍可认为本文细观数值方法能够有效地模拟荷载后混凝土的导热行为。
图4 细观模拟结果与Zhang等[1]试验结果对比
Fig.4 Comparison between the simulation results with Zhang et al’s[1]experimental results
3.2 力学行为分析
图5、图6分别给出了不同荷载水平σr下，混凝土中损伤分布云图及宏观应力、损伤水平随应变的变化关系曲线。由图5可知，在低荷载水平(如σr=0.25或σr=0.40)下，混凝土的反应十分接近于弹性，其内部产生的微裂纹很少，因而宏观损伤水平非常低(见图6)。随着荷载水平增大，混凝土薄弱区(即界面区)微裂纹增多，其宏观损伤水平迅速增大。当混凝土宏观应力达到其抗压强度(σr=1.0)后，损伤在界面区继续发展并且进入砂浆，并逐渐贯通形成宏观裂纹，其宏观损伤水平继续增大。最终，混凝土中裂纹发展减缓，相应地宏观损伤水平逐渐趋近于1。
图5 不同荷载水平下混凝土中损伤分布云图
Fig.5 Damage distribution within concrete under different loading levels
图6 混凝土压缩应力、宏观损伤水平随应变的变化曲线
Fig.6 Numerical simulated relationships among the compressive stressσ, macroscopic damage level and the strainεof concrete
3.3 热流密度与有效导热系数
基于已验证的细观数值方法，获得不同压缩荷载水平σr下混凝土中的稳定热流密度分布如图 7，另外，取混凝土试件中心位置，作出不同荷载(损伤)水平下热流密度随深度的变化曲线如图8。
由热流密度与导热系数的关系可知，对于稳态下的给定尺寸的截面，热流密度仅随导热系数变化，因而热流密度的分布实际了反映了混凝土内部不同位置的导热性能。考察图7不难发现，在低荷载水平(如σr= 0或σr=0.40)下，混凝土宏观损伤很小，力学荷载对混凝土导热性能的影响很小，此时混凝土中热流密度分布与无荷载作用时几乎相同。这一点在图 8中体现尤为明显。由于骨料导热性能优于砂浆及界面，热流优先通过骨料，从而使得骨料热流密度明显大于砂浆，骨料之间形成“热桥”[16]。因而，图8中热流密度随深度的变化曲线十分曲折。随着应力进一步增大，损伤亦进一步发展。当荷载水平σr＞0.7时，混凝土中损伤迅速发展，与完好混凝土相比，采用串联方法获得的热流密度分布变得更加均匀，说明“热桥效应”明显减弱。当混凝土宏观应力达到其抗压强度(σr= 1.0)后，宏观裂纹形成，由于裂纹被空气充满，导热能力很低，从而形成了明显的热量阻隔带。这在并联模型结果中体现十分明显，对于串联模型结果，由于其对损伤较为敏感，由于裂纹的阻隔，此时混凝土内部热流密度均较低，整体导热性能明显下降。最终，混凝土中损伤发展变缓，相应热流密度分布趋于稳定。
图7 不同荷载水平下混凝土中热流密度分布
Fig.7 Evolution of heat flux distribution within concrete under different compressive loading levels
由图 8明显可知，随着荷载(损伤)水平增大，混凝土中不同深度处热流密度减小。当荷载水平σr> 1.0时，由于混凝土中损伤发展明显，相比于并联模型，串联模型得到的热流密度明显较低。
图8 不同损伤水平下试件内部热流密度随深度的变化曲线
Fig.8 Heat flux-depth curve within the concrete specimen under different damage levels
图9为由并联模型和串联模型获得的混凝土试件的有效导热系数随宏观损伤水平的变化。由图可知，当混凝土反应接近于弹性时，由于其宏观损伤水平很低，串联模型与并联模型结果相差很小。随着荷载水平增大，混凝土宏观损伤水平增大，两种模型获得的有效导热系数均降低。其中，并联模型有效导热系数随宏观损伤水平近似线性降低，而串串联模型有效导热系数在损伤水平小于0.5时下降剧烈，之后逐渐趋于稳定。
3.4 混凝土内部温度场
将混凝土的顶面边界替换为 ISO834升温曲线[27]，其余各边设为绝热边界，得到1 h后不同荷载水平下混凝土试件的温度场如图10，并绘制试件中部处温度-深度变化曲线如图11。由图10及图11可知，随着荷载水平升高，混凝土内的温度场变得越来越不均匀。由于力学损伤减缓了热量传导，对于试件内某一位置，其温度随荷载水平升高而降低。另外，由于串联模型对于损伤更加敏感，峰值应力后，基于其获得的温度随荷载水平下降程度更大。
图9 混凝土试件宏观损伤水平对其有效导热系数的影响
Fig.9 Effect of the macroscopic damage level on its ETCkeff
图10 不同荷载水平下混凝土内部温度分布
Fig.10 Temperature field within concrete under different loading levels
图11 不同荷载水平下混凝土内部温度-深度曲线
Fig.11 Temperature-depth curves within concrete under different loading levels
4 结论
本文考虑混凝土细观结构非均质性的影响，将其视为由砂浆、骨料和界面组成的三相复合材料，提出并建立了力学-热学耦合作用的细观分析方法与模型。首先对试件进行单轴压缩力学行为模拟分析，获得混凝土试件中的损伤分布情况；进而结合复合材料均匀化方法获得损伤单元的导热特性，基于此对荷载作用后混凝土的导热行为进行了细观数值研究，得出了如下结论：
(1) 细观模拟结果与试验结果吻合良好，验证了本文方法的合理性与有效性；
(2) 对于本文研究工况，串联模型效果优于并联模型；然而，对于其他情况，需要更多物理试验验证；
(3) 随着荷载水平升高，混凝土有效导热性能降低，混凝土中温度场分布更加不均匀。
需要说明的是，本文研究工作仅针对压缩荷载情形，因此需要开展更多的物理试验，明确荷载作用对损伤后混凝土导热性能的影响，以期指导复合材料力学模型的选取，建立更加准确、适用范围更广的数值方法。另外，本文方法作为一种单向耦合方法，仅考虑了混凝土力学性能退化对其导热性能的影响，实际上，温度传导产生温度梯度，会进一步引起混凝土开裂，从而影响其力学性能。相关研究在后续工作中分析与讨论。
参考文献：
[1]Zhang W, Min H, Gu X, et al. Mesoscale model for thermal conductivity of concrete [J]. Construction and Building Materials, 2015, 98: 8―16.
[2]Tang S, Tang CA, Liang Z, et al. Numerical study of the influence of material structure on effective thermal conductivity of concrete [J]. Heat Transfer Engineering,2012, 33(8): 732―747.
[3]英明鉴, 李易, 陆新征, 等. 某超高层混凝土结构火灾作用下的力学响应分析[J]. 工程力学, 2015, 32(5):19―27.Ying Mingjian, Li Yi, Lu Xinzheng, et al. The mechanical response analysis of a super tall concrete structure under fire [J]. Engineering Mechanics, 2015,32(5): 19―27. (in Chinese)
[4]李志卫, 肖建庄, 谢青海. 高温后高强混凝土受压动态损伤[J]. 工程力学, 2017, 34(2): 78―84.Li Zhiwei, Xiao Jianzhuang, Xie Qinghai. Compressive dynamic damage of high-strength concrete after elevated temperatures [J]. Engineering Mechanics, 2017, 34(2):78―84. (in Chinese)
[5]Erdem S, Dawson A R, Thom N H. Impact load-induced micro-structural damage and micro-structure associated mechanical response of concrete made with different surface roughness and porosity aggregates [J]. Cement and Concrete Research, 2012, 42(2): 291―305.
[6]Wyrzykowski M, Lura P. The effect of external load on internal relative humidity in concrete [J]. Cement and Concrete Research, 2014, 65: 58―63.
[7]金浏, 张仁波, 杜修力. 低应力水平下混凝土中氯离子扩散行为多尺度分析方法[J]. 工程力学, 2017,34(3): 84―92.Jin Liu, Zhang Renbo, Du Xiuli. Multi-scale analysis for the chloride diffusivity in concrete subjected to low-level stress [J]. Engineering Mechanics, 2017, 34(3): 84―92.(in Chinese)
[8]Wyrwa? J, Marynowicz A, ?wirska J. Effective thermal conductivity of porous building materials - analysis and verification [J]. Bauphysik, 2008, 30(6): 431―433.
[9]杜修力, 金浏, 张仁波. 压缩荷载作用下混凝土中氯离子扩散行为细观模拟[J]. 建筑材料学报, 2016,19(1): 65―71.Du Xiuli, Jin Liu, Zhang Renbo. Meso-scale simulation of chloride diffusivity in concrete subjected to compressive stress [J]. Journal of Building Materials,2016, 19(1): 65―71. (in Chinese)
[10]杜修力, 金浏, 张仁波. 力学荷载对混凝土中氯离子渗透扩散行为影响评述[J]. 建筑结构学报, 2016,37(1): 107―125.Du Xiuli, Jin Liu, Zhang Renbo. Review on effect of external mechanical loadings on chloride penetration and diffusion into concrete [J]. Journal of Building Structures, 2016, 37(1): 107―125. (in Chinese)
[11]C?erny R, Mad Ra J, Pod Bradská J, et al. The effect of compressive stress on thermal and hygric properties of portland cement mortar in wide temperature and moisture ranges [J]. Cement and Concrete Research, 2000, 30(8):1267―1276.
[12]Perkowski Z. Change of thermal conductivity of concrete caused by brittle damage evolution [J]. Bauphysik, 2008,30(6): 434―437.
[13]Vejmelková E, Padevět P, ?erny R. Effect of cracks on hygric and thermal characteristics of concrete [J].Bauphysik, 2008, 30(6): 438―444.
[14]Hasselman D P H. Effect of cracks on thermal conductivity [J]. Journal of Composite Materials, 1978,12(4): 403―407.
[15]Bary B. Estimation of poromechanical and thermal conductivity properties of unsaturated isotropically microcracked cement pastes [J]. International Journal for Numerical and Analytical Methods in Geomechanics,2011, 35(14): 1560―1586.
[16]Shen L, Ren Q, Xia N, et al. Mesoscopic numerical simulation of effective thermal conductivity of tensile cracked concrete [J]. Construction and Building Materials, 2015, 95: 467―475.
[17]Xi Y, Nakhi A. Composite damage models for diffusivity of distressed materials [J]. ASCE Journal of Materials in Civil Engineering, 2005, 17(3): 286―295.
[18]Kaddouri W, El Moumen A, Kanit T, et al. On the effect of inclusion shape on effective thermal conductivity of heterogeneous materials [J]. Mechanics of Materials,2016, 92: 28―41.
[19]?avija B, Pacheco J, Schlangen E. Lattice modeling of chloride diffusion in sound and cracked concrete [J].Cement and Concrete Composites, 2013, 42: 30―40.
[20]Jin L, Zhang R, Du X, et al. Investigation on the cracking behavior of concrete cover induced by corner located rebar corrosion [J]. Engineering Failure Analysis, 2015,52: 129―143.
[21]杜修力, 韩亚强, 金浏, 等. 骨料空间分布对混凝土压缩强度及软化曲线影响统计分析[J]. 水利学报, 2015,46(6): 631―639.Du Xiuli, Han Yaqiang, Jin Liu, et al. Statistical investigation on effects of aggregate distribution on concrete compression strength and the descending part of stress-strain curve [J]. Journal of Hydraulic Engineering,2015, 46(6): 631―639. (in Chinese)
[22]Du X, Jin L, Ma G. A meso-scale numerical method for the simulation of chloride diffusivity in concrete [J].Finite Elements in Analysis and Design, 2014, 85: 87―100.
[23]Wriggers P, Moftah S O. Mesoscale models for concrete:homogenisation and damage behaviour [J]. Finite Elements in Analysis and Design, 2006, 42(7): 623―636.
[24]Lubliner J, Ollivier J, Oller S, et al. A plastic-damage model for concrete [J]. International Journal of Solids and Structures, 1989, 25(3): 299―326.
[25]Lee J, Fenves G. Plastic-damage model for cyclic loading of concrete structures [J]. ASCE Journal of Engineering Mechanics, 1998, 124(8): 892―900.
[26]Zhao J, Zheng J, Peng G, et al. A meso-level investigation into the explosive spalling mechanism of high-performance concrete under fire exposure [J].Cement and Concrete Research, 2014, 65: 64―75.
[27]International Organization for Standardization. ISO 834-1 Fire resistance test on elements of building construction [S]. 1999.