准脆性材料（如混凝土、岩石）的热传导问题在工程应用中常见且重要 ^[1]. 在实际工程中，由于温差普遍存在，热传导几乎时刻在发生^[2]. 传热引起温度的变化可能会改变材料的性能与力学行为^[3]，甚至导致材料失效破坏. 在混凝土结构中，温度由于与早龄期的水化过程息息相关而影响混凝土性能演化^[4]当混凝土结构遇到温度在短时间变化剧烈^[5] 的情况时，过大的温度梯度会导致结构热开裂. 在地热能源开发中，将冷水循环注入储热层引发岩石裂缝^[6]，增加了储层的渗透性，对生产有着积极影响. 在核废料处理中，由于衰变释放的热量导致岩石产生裂缝可能会引起放射性物质泄漏^[7].

学者们采用多种试验技术研究准脆性材料的热力学性能^[8]，如声发射技术（acoustic emission technology，AE）^[9]、扫描电镜（scanning electron microscope，SEM）^[10]、热应力装置（thermal stress device，TSD）^[11]等. 物理试验耗费时间且难以观察热开裂过程中裂纹的萌生与扩展. 为了解决这些问题，数值模拟方法被运用于探索准脆性材料的热开裂机理中. 热开裂数值模拟方法大致分为连续介质模型、非连续介质模型和连续−离散耦合模型. 连续介质模型扩展了有限元法（extended finite element method，XFEM）^[12-13]与颗粒有限元法（particle finite element method，PFEM）^[14]，在解决模拟热开裂问题中被广泛采用. 虽然包括相场方法^[15]和近场动力学（peridynamics，PD）^[16]在内的先进连续介质方法被运用到准脆性材料热开裂模拟中，但是多数连续介质模型难以处理多裂纹问题. 如离散元法（discrete element method，DEM）^[17]、非连续变形分析（discontinuous deformation analysis，DDA）^[18-19]的非连续介质模型和以连续−离散耦合方法（finite and discrete element method，FDEM）为代表的连续−离散耦合模型^[20-22]能够处理多裂纹问题，缺点是计算量较大，效率低^[23]. 格构离散单元法（lattice discrete element method，LDEM）运用格构模型改进标准DEM，大大提高了计算效率^[24]. 基于Hrennikoff^[25]提出的格构框架思想，Nayfeh等^[26]将杆单元构成的基本模块离散连续体结构为LDEM原型. LDEM经过不断改进^[27]，已经被成功运用于混凝土^[28]、岩石^[29]、复合材料^[30]等结构开裂的计算中. Grassl^[31]基于格构思想提出开裂混凝土中的物质流动模型，Wang等^[32-33]运用经典格构模型研究氯离子扩散行为，吴迪^[34]采用二维格构模型分析混凝土中的热传导问题，均说明格构模型具有计算物质传输与传热的潜力.

本研究提出LDEM热力耦合模型，通过杆单元基本模块与连续体的热传导等效换算，实现热传导模拟；并在LDEM求解方程中集成热应力，实现准脆性材料热开裂过程的模拟.

LDEM的基本原理是将连续体离散为质量集中于两端的、在空间周期性排列的杆单元，本研究采用Nayfeh等^[26]提出的LDEM模式. 如图1所示，离散化采用的基本立方体模块由20根杆单元（分为2类：纵向杆、对角杆）与9个节点组成，相邻2个基本模块的体心节点由1根纵向杆单元连接. 杆单元仅在轴向承受荷载，每个节点有3个自由度，表示为位移矢量的3个分量. 假定立方体模块对应的连续体边长为L，纵向杆单元长度L_l=L，对角杆单元长度 $ {L_{\text{d}}} = \sqrt 3 L/2 $. 基本模块整体质量m_all与连续体质量相等，即m_all=ρL³，其中ρ为连续体密度由于体心节点在模块内部，角节点为8个模块共有，得到体心节点质量m_cen=ρL³/2，角节点质量m_cor=ρL³/16. 杆单元的横截面积由整体刚度等效计算，在各项同性材料中，纵向杆单元与对角杆单元的等效轴向刚度^[26]分别为

(1) $ \left.\begin{array}{l}{E}_{\text{l}}^{\text{A}}={A}_{\text{l}}E=\phi E{L}^{2}\text{，}\\ {E}_{\text{d}}^{\text{A}}={A}_{\text{d}}E=\dfrac{2}{\sqrt{3}}\delta \phi E{L}^{2}\text{；}\\ \phi =(9+8\delta )/(18+24\delta )\text{，}\\ \delta =9\nu /(4-8\nu ).\end{array}\right\} $

式中：A_l、A_d分别为纵向杆单元与横向杆单元的横截面积，E为被离散连续体的弹性模量，ν为泊松比. 由式（1）得到纵向杆单元与对角杆单元的横截面积关系式为

(2) $ {A_{\text{d}}} = \frac{2}{{\sqrt 3 }}\delta {A_{\text{l}}}{\text{ = }}\frac{2}{{\sqrt 3 }}\delta \phi {L^2} . $

LDEM适用于多相混合的含孔材料，其离散连续体的过程基于连续体的各相空间分布进行^[35]. 以混凝土为例，如图2所示，将连续体划分为空间中的立方体元胞，并依据元胞的空间位置判断其材料属性，将整体区分为骨料、砂浆、界面过渡区（ITZ）与孔隙组成的多相结构. 在元胞的角点与体心生成节点，并将节点集按照LDEM基本模块的排布方式通过杆单元进行连接，根据杆单元所属元胞赋予相应的材料参数. 位于孔隙元胞内的节点不参与单元划分，在模型中形成初始缺陷. LDEM模拟裂纹的形成与扩展是通过在一定开裂准则下的杆单元损伤去除进行的. 在每步计算中，达到极限条件的杆单元将从整体结构中删除，相邻节点不再连续，之后再进入下步计算.

LDEM基本模块的传热过程由组成模块的所有杆单元的传热过程共同决定. 根据傅里叶定律，杆单元横截面上的热流量q_i计算式为

(3) $ {q_i} = - {k_i}\frac{{\Delta {\theta _i}}}{{{L_i}}}. $

式中：k_i为杆单元的热传递系数，Δθ_i为杆单元两端的温差，L_i为杆单元的长度. 传热的等效过程如图3所示. 连续体截面在LDEM模块中等效为4个对角杆单元斜截面与2个纵向杆单元横截面. 连续体单一方向上的温差Δθ转化为LDEM模块对应节点的温差. 为了保证基本模块内传热的均匀性，纵向杆单元与对角杆单元中传输的热量相等，即q_lA_l=q_dA_d. 由Δθ_d=0.5Δθ_l=0.5Δθ，2L_d= $\sqrt 3 $L_l，A_d= ${2}\delta {A_{\text{l}}}/{\sqrt 3}$，得到2类杆单元热传导系数的关系为

(4) $ {k_{\text{d}}} = \dfrac{1}{{2\delta }}{k_{\text{l}}}. $

在单一热流方向上连续体与LDEM传热的等效关系为

(5) $ k\frac{{\Delta \theta }}{L}{L^2} = 2 {A_{\text{l}}}{k_{\text{l}}}\frac{{\Delta \theta }}{L}+4 \sqrt 3 {A_{\text{d}}}{k_{\text{d}}}\frac{{\Delta \theta }}{L}. $

式中：k为连续体热传导系数， $\sqrt 3 $为对角杆单元斜截面积与横截面积比值. 杆单元与连续体导热系数的关系为

(6) $ \begin{array}{l} {k_{\text{d}}} = \dfrac{k}{{12\phi \delta }}, \; {k_{\text{l}}} = \dfrac{k}{{6\phi }}. \end{array} $

连续体在给定的质量m条件下的温度变化表达式为

(7) $ \Delta \theta = \frac{{{Q_{_{\text{all}}}}}}{{c m}}. $

式中：Q_all为流入连续体的热量，c为连续体比热. LDEM基本模块质量m_mod=ρL³=m，其中m为连续体质量，流入模块的热量Q根据式（4）有Q=Q_all，且根据式（3）得到纵向杆单元与对角杆单元中的热量均匀分布，为了保证LEDM基本模块与连续体的温度变化相等，杆单元比热满足关系式

(8) $ {c_{\text{l}}} = {c_{\text{d}}} = c. $

式中：c_l为纵向杆单元比热，c_d为对角杆单元比热.

当各向同性的连续体内受到均匀变温Δθ作用时，若不受外部约束，各点产生正应变：

(9) $ {\varepsilon }_{x}={\varepsilon }_{y}={\varepsilon }_{z}=\alpha \Delta \theta \text{，}{\gamma }_{yz}={\gamma }_{zx}={\gamma }_{xy}=0. $

式中：ε_x、ε_y、ε_z分别为x、y、z方向上的正应变，γ_yz、γ_zx、γ_xy分别为x、y、z方向上的切应变，α为连续体线膨胀系数. 当LDEM中受到同样的Δθ作用时，杆单元发生轴向应变ε=αΔθ. 根据几何构型，当LDEM基本模块与连续体变形一致时，线膨胀系数关系式为

(10) $ {\alpha _{\text{l}}} = {\alpha _{\text{d}}} = \alpha . $

式中：α_l为纵向杆单元线膨胀系数，α_d为对角杆单元线膨胀系数. 当受到外部约束，上述变形不能自由发生时，杆件内部产生热应力：

(11) $ {\sigma _i} = - E {\alpha _i} \Delta \theta . $

原型LDEM采用牛顿第二定律对每个节点进行求解^[27]，在此基础上引入杆单元中由于温度产生的内力，在节点上有

(12) $ {\boldsymbol{M\ddot u}}+{\boldsymbol{C\dot u}}+{\boldsymbol{F}}{\text{(}}t{\text{)}}+{\boldsymbol{W}}(t) - {\boldsymbol{P}}{\text{(}}t{\text{)}} = {\mathbf{0}}， $

(13) $ {\boldsymbol{C}} = 2{\text{π }}\xi {f_{\text{p}}}{\boldsymbol{M}}. $

式中：M、C分别为质量矩阵与阻尼矩阵， $ {\boldsymbol{\ddot u}} $、 $ {\boldsymbol{\dot u }}$分别为节点的加速度与速度，向量F(t)为非温变内力，W(t)为热应力积分项，P(t)为外力项，ξ为阻尼系数，f_p为能谱峰值频率.

LDEM使用Kosteski等^[28]引入的双线性本构模型，以体现拉伸时微裂纹形成与扩展过程中杆单元承载能力的不可逆损伤. 在压缩载荷的情况下，材料表现为线性弹性. 压缩载荷导致材料失效视为由间接牵引作用造成，其抗压极限强度通常比抗拉强度大5~10倍. 双线性本构模型的荷载−应变曲线如图4所示. 图中，段OA斜率杆单元等效刚度EA_i；点A对应的ε_p为临界应变，表示微裂纹开始形成前的最大应变；点B对应的ε_r为极限应变，即杆单元失效去除时的应变. 在AB上取任意点P表示杆单元发生损伤但未完全失效，区域OAP的面积E_D为发生损伤时所耗散的能量密度，区域OPC的面积E_s为杆单元中剩余的弹性能量密度. 当总能量完全耗散，即OAP与OAB重合时，杆单元去除，宏观裂纹产生.

在断裂过程中，LDEM所耗散的断裂能须与连续体等效，当发生平行面开裂时，连续体中耗散的能量为

(14) $ {\textit{Γ}} = {G_{\text{f}}}{L^2}. $

式中：G_f为材料的断裂能. LDEM模块中平行面等效为4根边缘纵向杆单元、1根中心纵向杆单元和4根对角杆单元的断裂，其耗散的能量为

(15) $ {{\textit{Γ}} _{{\text{LDEM}}}} = {c_{\text{A}}}\left( {4 \times \frac{1}{4}+1+4 \times \left( {\frac{2}{{\sqrt 3 }}\delta } \right)} \right){G_{\text{f}}}{L^2}. $

式中： $2\delta/\sqrt 3$为对角杆单元与纵向杆单元的横截面积比，c_A为引入的能量等效系数，由式（14）与（15）相等计算得到c_A≈0.139. 杆单元的等效断裂面积为

(16) $ \left. \begin{array}{l} A_{\text{l}}^{\text{f}} = {c_{\text{A}}}{L^2}{\text{ = }}0.139\;0{L^2}, \\ A_{\text{d}}^{\text{f}} = {c_{\text{A}}}\dfrac{2}{{\sqrt 3 }}\delta {L^2}{\text{ = }}0.180\;6{L^2}. \\ \end{array} \right\} $

由线性断裂力学得到双线性本构模型所围区域为

(17) $ \int_0^{{\varepsilon _{\text{r}}}} F (\varepsilon ){\rm{d}}\varepsilon = \frac{{{\varepsilon _{\text{r}}}{\varepsilon _{\text{p}}}E{A_i}}}{2} = \frac{{{G_{\text{f}}}A_i^{\text{f}}}}{{{L_i}}}. $

极限应变ε_r的计算式为

(18) $ {\varepsilon _{\text{r}}} = \frac{{{G_{\text{f}}}}}{{{\varepsilon _{\text{p}}}E}}\left( {\frac{{A_i^{\text{f}}}}{{{A_i}}}} \right)\left( {\frac{2}{{{L_i}}}} \right). $

应力强度因子

(19) $ {K_{\text{c}}} = \sqrt {{G_{\text{f}}}E} = {\sigma _{\text{p}}}Y\sqrt {{\text{π}}d} = E{\varepsilon _{\text{p}}}Y\sqrt {{\text{π}}d} . $

式中：σ_p为临界应力，Y为取决于试样边界条件和临界裂纹尺寸d的无量纲参数.

假定材料特征长度d_eq=dπY²，当裂纹长度大于d_eq时，将发生不稳定扩展. 由式（18）得到

(20) $ {d_{{\text{eq}}}} = \dfrac{{{G_{\text{f}}}}}{{\varepsilon _{\text{p}}^{\text{2}}E}}. $

临界应变与极限应变的表达式分别为

(21) $ {\varepsilon _{\text{p}}} = \sqrt {\frac{{{G_{\text{f}}}}}{{{d_{{\text{eq}}}}E}}} ， $

(22) $ {\varepsilon _{\text{r}}} = {\varepsilon _{\text{p}}}{d_{{\text{eq}}}}\left( {\frac{{A_i^{\text{f}}}}{{{A_i}}}} \right)\left( {\frac{2}{{{L_i}}}} \right). $

在Abaqus环境中，以顺序热力耦合方式实现LDEM热力耦合模拟. 采用DC1D2单元进行热传导计算，T3D2单元进行热应力与开裂计算. 由于上述单元无法将质量集中于节点，须对杆单元密度进行等效换算. 根据LDEM模块体心节点与顶角节点的质量分配以及杆的连接关系，质量等效关系式为

(23) $ \left. \begin{array}{l} \dfrac{1}{2}\rho {L^3} = 6 \times \dfrac{1}{2}{\rho _{\text{l}}}{A_{\text{l}}}L+8 \times \dfrac{1}{2}{\rho _{\text{d}}}{A_{\text{d}}}{L_{\text{d}}}, \\ \dfrac{1}{{16}}\rho {L^3} = 3 \times \dfrac{1}{2}{\rho _{\text{l}}}\dfrac{1}{4}{A_{\text{l}}}L+\dfrac{1}{2}{\rho _{\text{d}}}{A_{\text{d}}}{L_{\text{d}}}. \\ \end{array} \right\} $

考虑到模块内部质量的均匀性，假定ρ_d=ρ_l，计算得到杆单元的密度为

(24) $ {\rho _{\text{l}}} = {\rho _{\text{d}}} = \frac{\rho }{{6\phi {\text+}8\delta \phi }}. $

当LDEM离散连续体时，边缘缺少部分纵向杆单元，须对边缘面的节点进行附加质量补偿，补偿质量m_c的计算式为

(25) $ {m_{\text{c}}} = \frac{n}{2}{\rho _{\text{l}}}{A_{\text{l}}}L. $

式中：n为该节点缺少的纵向杆单元数量. 材料库中提供脆性开裂本构模型（brittle cracking constitutive model），以材料在I型下的断裂能G_f为断裂参数，且允许计算过程中去除失效单元，该模型与LDEM的双线性本构模型类似. 因此，可以将LDEM热力耦合模型引入Abaqus中进行计算.

以平板的冷却为例，对LDEM热力耦合模型在不同边界条件下的热传导与热弹性应力进行验证.

平板的初温为50 ℃，在一侧施加恒定边界温度0 ℃对平板进行降温，上下为绝热边界并完全约束. 该问题的数学表达为

(26) $ \left. \begin{array}{l} \dfrac{{\partial \theta (x,t)}}{{\partial t}} = D\dfrac{{\partial \theta {{(x,t)}^2}}}{{\partial {x^2}}}, \\ \theta (x > 0,t = 0) = {\theta _0}, \\ \theta (x = 0,t > 0) = {\theta _{\text{c}}}. \\ \end{array} \right\} $

式中：D为导温系数，D=k/ρc. 该问题的理论解^[34]为

(27) $ \theta (x,t) = {\theta _0}+\left( {{\theta _{\text{c}}} - {\theta _0}} \right)\left( {1 - {\rm{erf}}\left( {\frac{x}{{2\sqrt {Dt} }}} \right)} \right). $

以y=0为监测线，监测线x方向的应力解析解为

(28) $ {\sigma _x} = - E\alpha \left( {\theta (x,t) - {\theta _0}} \right). $

平板的热力学参数如下. 导热系数k=0.94 W/(m·℃)，密度ρ=2 440 kg/m³，比热c=970 J/(kg·℃)，线膨胀系数α=1.5×10⁻⁵，弹性模量E=30 GPa，泊松比υ=0.167. 离散的LDEM模块L=1.0 mm. 如图5所示为t=10、100、300 、600 s时，监测线上温度和应力的解析解与LDEM数值解. 对比可知，平板在不同时刻和不同位置时，温度与应力的理论曲线与数值曲线基本重合.

为了研究LDEM网格尺寸对温度和应力结果的影响，对平板问题进行不同网格尺寸的划分（LDEM模块L=0.5、1.0、2.0、5.0 mm），对比监测线上的温度与应力结果. 如图6所示为t=600 s时不同网格尺寸下的温度分布情况，如图7示为t=600 s时，监测线上温度与应力解析解与数值解的对比. 可以看出，数值结果随着网格尺寸的减小收敛于解析解.

为了验证LDEM热力耦合方法在第三类边界条件问题中的适用性，调整平板问题的边界条件：平板初温为50 ℃，左右两侧环境温度为0 ℃，表面热交换系数h=14.7 W/(m²·K)，上下边绝热且完全约束，力学参数保持一致. 采用有限元（finite element method，FEM）、LDEM方法分别求解上述问题，网格划分尺寸均为1 mm. 如图8所示为t=600 s时，2种方法得到的温度分布，为了观察LDEM变形的合理性，将2种方法得到的变形同时放大了4 000倍. 可以看出，FEM、LDEM得到的温度分布与温度影响下产生的变形基本一致. 如图9所示，以y=0为监测线，对比2种方法在监测线上的温度与应力. 可以看出，两者的温度与应力曲线基本重合，说明LDEM热力耦合方法在第三类边界条件问题中同样适用.

为了验证LDEM热力耦合模型在热开裂模拟中的适用性，考虑以下2类问题：1）在均匀温度变化下，由于材料不均匀性产生的热失配应力开裂；2）在均匀的材料中，由温度梯度产生的热应力开裂.

数值模拟以Abdalla^[36]进行的钢筋混凝土结构在高温影响下的热开裂物理试验为基础. 混凝土中钢筋通常由线膨胀系数较高的钢材制成，在均匀温升时混凝土保护层产生压力，引起开裂. 试验采用由2种材料构成的同轴圆柱，如图10所示. 图中，a为内芯半径，b为钢筋混凝土整体半径，r为保护层内任意一点到圆心的距离. 模拟采用的热力学参数如表1所示，建立LDEM计算模型，采用0.001 m的LDEM模块进行网格划分. 本研究将模拟以1 ℃/s的速率使试件从0 ℃均匀加热至100 ℃的过程，并将LDEM计算得到的开裂模式对比试验的开裂模式. Abdalla^[36]给出的钢筋混凝土材料不均匀性热开裂问题环向应力解析解为

(29) $ \left. \begin{array}{l} {\sigma _{\text{τ }}} = \dfrac{{{a^2}\left( {{b^2}+{r^2}} \right)}}{{{r^2}\left( {{b^2} - {a^2}} \right)}}p, \\ p = \dfrac{{\left( {{\alpha _{\text{a}}} - {\alpha _{\text{b}}}} \right)\Delta \theta {E_{\text{a}}}}}{{{E_{\text{a}}}/{E_{\text{b}}}\left( {\beta +{\mu_{\text{b}}}} \right)+\left( {1 - {\mu_{\text{a}}}} \right)}}, \\ \beta = \dfrac{{{b^2}+{a^2}}}{{{b^2} - {a^2}}}. \end{array} \right\} $

式中：α_a、α_b分别为内芯、保护层的线膨胀系数，E_a为内芯的弹性模量，μ_a、μ_b分别为内芯、保护层的泊松比.

如图11所示为LDEM计算得到的在不同温升时，混凝土保护层沿半径方向的环向应力σ_τ 与解析解的对比，两者一致性较好. 可以看出，混凝土保护层r=a处应力最大，为试验中预期断裂开始的位置. 热开裂模拟结果如图12所示. 可以看出，开裂模式与试验结果一致，即裂纹在混凝土与钢筋的界面处产生，沿径向往外边界扩展.

验证基于Jansen等^[37]进行的方形花岗岩热开裂试验. 试件为边长15 cm的正方形，中心钻孔直径为1 cm，采用0.05 cm的LDEM模块划分网格. 法向约束试件的上下边界，左右为自由边界. 试件的初始温度以及四周的温度边界条件均为20 ℃，加热中心孔，使其温度在1 h内由20 ℃均匀上升至200 ℃. 热力学参数与试验中花岗岩参数保持一致，如表2所示. 计算结果如图13所示，随着中心孔温度的升高，热由孔向四周边界传递. 孔附近的岩石由于温度较高，发生较大的热膨胀变形，外边界附近岩石温升幅度较小，体积膨胀较小. 由于上下边界的法向约束作用，岩石内部的体积膨胀引起的拉伸变形只能在水平方法发生，在上下边界产生较大的拉应力. 随着加热过程的进行，先是在试件上下边缘产生拉伸破坏裂纹，之后上下裂纹同时向着中心孔延伸. 由于中心孔为加热源，附近温度梯度较大，产生的变形不一致，孔周围产生细小的裂纹. 开裂模式与如图14所示的方形花岗岩热开裂试验结果^[37]和颗粒流模拟^[38]得到的开裂过程一致.

在工程应用中，混凝土在自身细观结构与外界温度变化的共同影响下，极易产生温度裂缝，温度−应力试验及其数值仿真能够为混凝土结构设计与温度裂缝预防提供重要依据，为此，研究LDEM热力耦合模型在细观混凝土温度−应力试验仿真中的应用.

温度应力试验机（temperature-stress testing machine，TSTM）是测试混凝土在单轴约束下抗裂性能的装置，卢春鹏等^[39]采用TSTM分析早龄期混凝土的温度应力，本研究基于此实验进行仿真. 混凝土配合比与文献[39]的温度应力试验保持一致，水灰比为0.43，单位立方米材料用量分别为水125 kg、水泥189 kg、粉煤灰102 kg、砂672 kg、小石637 kg、中石 637 kg. 混凝土骨料体积分数ϕ ≈40%. 为了研究细观结构对混凝土热开裂的影响，设置不同骨料体积分数ϕ =0%、10%、20%. 考虑到计算成本与骨料投放，试件尺寸选取试验试件中部位置，为600 mm×150 mm×150 mm. 模型骨料投放如图15所示. 以3 mm的LDEM模块划分混凝土模型，并将骨料外层模块作为界面过渡区（interface transition zone, ITZ）. 仿真过程：当混凝土7 d养护龄期时，法向约束混凝土试件上下边界，在初温30 ℃的试件四周进行强制降温使其内部温度以1 ℃/h的速率降温直至开裂破坏. 热力学参数参考文献[40]，如表3所示.

如图16所示为各个试件在中心温降Δθ=20 ℃（未开裂）时的温度与变形云图. 可以看出，由于骨料的热力学性能与砂浆基质以及ITZ存在差别，导致混凝土内部骨料周围的温度与变形分布出现明显“弯折”，随着骨料体积分数的上升，分布的不均匀现象愈加明显. 如图17所示为试件上下约束面在温降直到开裂过程中的法向位移U. 可以看出，在同样的温降幅度下，随着骨料体积的增加，试件的变形依次减小，说明骨料能够降低混凝土在温降作用下的收缩变形. 如图18所示为各个试件上下面约束应力对比. 可以看出，采用试验级配即ϕ=40%的应力曲线与试验结果符合较好，能够反映TSTM试验规律. 试验曲线下降段仅表示断裂发生，没有记录开裂后继应力. 对比不同的试件，其中纯砂浆试件并未发生开裂，含骨料试件随着骨料体积分数的增加，开裂时间提前，即开裂温度上升. 这是由于混凝土的抗拉强度主要受砂浆与ITZ的强度控制，不同试件的强度相近，骨料的增加提高了混凝土的弹性模量，加快了温降过程中的应力增长，导致达到混凝土抗拉强度所需的温降幅度减小. 如图19所示，不同试件的断裂模式是相似的，裂纹从骨料周围的ITZ产生，经由砂浆基质联通，形成平行于约束面的裂纹面，与试验相符. 当骨料体积分数ϕ=40%时，ITZ间距较小，裂纹较为集中，形成1个裂纹面；当ϕ=10%时，在多个部位发生断裂，形成3个较为明显的裂纹面.

（1）通过平板算例验证LDEM热力耦合模型在不同热边界条件下，模拟热传导过程与热弹性应力的准确性及其对网格大小的不敏感性.

（2）通过钢筋混凝土同心圆柱与带孔岩板的热开裂试验模拟，验证LDEM热力耦合模型模拟热失配以及温度梯度产生热开裂的适用性.

（3）将LDEM热力耦合模型应用于混凝土TSTM试验的数值模拟，结果表明应力曲线以及开裂模式与试验结果吻合较好. 骨料体积分数对混凝土热开裂的影响研究表明：随着骨料体积分数的增加，混凝土收缩性减小，抗拉强度基本不变，开裂所需温降幅度降低.

（4）LDEM热力耦合模型能够较好处理混凝土、岩石、陶瓷等准脆性材料中的复杂热开裂问题. 格构离散单元法采用的立方体基本模块假设，只能进行六面体网格划分，因此本研究对象的形状较为简单. 为了更好地模拟真实工程结构的热开裂过程，计划开展非规则基本模块的相关研究.