混凝土材料具有制作工艺简单、价格低廉等优点，自问世以来，在桥梁、水利、地下等工程结构中应用广泛，而高寒地区混凝土的使用过程中冻融破坏是主要问题^[1-2]. 混凝土冻害是复杂的物理化学作用过程，其中的内部机理至今没有完全研究透彻^[3].

由于冻融循环中须考虑的因素众多，描述混凝土冻融循环数值模拟方程的建立过程处于探索阶段. Powers等^[4]提出经典静水压假说并定量确定静水压. Powers等^[5-7]发展和完善了渗透压力理论. Bažant等^[8]建立基于孔隙分布的混凝土抗冻耐久性的模型，计算所需参数需要在多种新型试验的基础上进行标定，实际应用效果不佳. Zuber等^[9-10]基于质量守恒定律和达西定律，建立饱和状态下的混凝土水-力耦合方程；Duan等^[11]考虑温度场，建立热-水-力耦合模型，但只考虑混凝土为均质材料时混凝土在受冻情况下的热力学响应情况. Zhou等^[3]在Duan等^[11]的研究基础上，考虑材料物理化学性质在冻融过程中的变化，提出更准确的多物理场耦合模型；模型结晶压力的计算未考虑相变点前、后的变化以及耦合模型所计算的细观模型为二维，与实际实验及工程应用中有一定的差别. Zuber等^[9-11]研究认为混凝土为均质材料，在细观尺度，混凝土砂浆基质与骨料之间的界面存在薄弱区域，称为界面过渡区（interfacial transition zone，ITZ）. 砂浆基质在此界面的结构与远离此处的浆体在成分、密度和其他物理性质方面有较大的不同，因此在细观尺度，应当将混凝土区分为砂浆、骨料、界面过渡区（ITZ）等三相介质^[3].

本文研究混凝土在冻融环境热-水-力多场耦合作用下三维细观尺度的热力学响应机制. 采用蒙特-卡洛方法，结合混凝土骨料富勒级配曲线，应用Matlab与Comsol Multiphysics接口进行二次开发，建立可以区分砂浆-骨料-界面过渡区（ITZ）三相的三维混凝土细观随机骨料模型. 将混凝土冻融损伤的热-水-力耦合模型引入混凝土细观尺度分析，利用有限元方法实现耦合模型，探讨砂浆渗透率、界面过渡区（ITZ）线膨胀系数及骨料体积分数对混凝土抗冻融性能的影响.

混凝土是复杂的多孔材料，具有不均匀性，内部包括水化产物和未水化的水泥颗粒部分、粗骨料和细骨料、添加剂和空气等^[12]. 在细观尺度上，有必要生成能够区分混凝土不同相的模型.

在混凝土的实际生产应用中，采用的骨料大多为卵石和碎石，形状接近球形，因此将骨料简化成球体^[13]. 根据混凝土骨料分布的级配^[14]，可得各个骨料粒径的分布：小石5~20 mm，中石20~40 mm，大石40~80 mm，特大石80~150 mm^[15].

物理试验所用的试样为40 mm×40 mm×160 mm立方体，考虑到计算量和时间成本以及试样的对称性，为了保持试件的长方体构型，避免尖角产生应力集中，数值试样选取物理模型试件尺寸的1/4，即20 mm×20 mm×160 mm立方体. 将粒径小于5 mm的细砂看作基质的一部分. 考虑到试件尺寸的限制，本文的骨料粒径仅考虑小石.

通过应用Matlab与Comsol Multiphysics接口进行二次开发，生成骨料、界面过渡区（ITZ）和砂浆基质三相三维细观模型. 模型的生成过程如下.

1）初始参数输入. 输入试件体积、骨料级配以及不同粒径段的体积分数.

2）骨料生成. 在骨料的体积分数范围内，按照粒径分布，先生成较大级配的骨料，再生成较小级配的骨料. 根据蒙特-卡洛理论，随机生成骨料球心坐标及半径，先判断骨料是否满足试件边界条件. 若满足，则再判断新生成的骨料与原生成骨料两两之间是否接触，引入骨料影响范围系数^[13]以防止骨料之间接触过近而影响计算结果，判断条件为

(1) $ \begin{split} &\sqrt {{{\left( {x\left( n \right) - x\left( i \right)} \right)}^2} + {{\left( {y\left( n \right) - y\left( i \right)} \right)}^2} + {{\left( {z\left( n \right) - z\left( i \right)} \right)}^2}} \geqslant \\ &\;\;\;\;\;\;\;\;\;\;\mu \left( {r\left( n \right) + r\left( i \right)} \right);\;i = 1,2,\cdots,n - 1. \end{split} $

式中： $\left( {x\left( i \right),y\left( i \right),z\left( i \right)} \right)$和 $r\left( i \right)$分别为已生成骨料的形心坐标和半径； $\left( {x\left( n \right),y\left( n \right),z\left( n \right)} \right)$和 $r\left( n \right)$分别为新生成骨料的坐标和半径，用以计算并判断是否与已生成骨料发生接触； $\mu $为骨料影响范围系数，取1.05.

3）界面过渡区生成. 欧阳利军等^[16]的研究表明，ITZ厚度为5~100 μm. 本文在骨料表面取50 μm作为ITZ厚度.

4）整体模型生成. 输入试件砂浆部分坐标的相关信息，砂浆、界面过渡区（ITZ）和骨料三部分装配成联合体，形成细观模型.

考虑计算成本和试件尺寸对骨料投放的限制，选择10%骨料体积分数、20%骨料体积分数等体积分数较低且为一级配的试样进行分析研究. 在骨料体积分数 ${\varphi _{\rm{a}}} $为20%试件的网格模型中，骨料部分网格数为17 604，界面过渡区（ITZ）部分网格数为9 108，砂浆基质部分网格数为63 334，总网格数为90 046. 三相介质网格尺寸为0.05~16.0 mm. 图1中网格质量指数I_q普遍大于0.7，表明所划分网格的质量指数较好，可以用于细观计算分析.

假设混凝土内部孔隙处于完全饱和状态，不考虑水化和溶解物理化学过程，且孔隙中只存在水和冰，即

(2) ${\varphi _{\rm{l}}} + {\varphi _{\rm{i}}} = 1.$

式中： ${\varphi _{\rm{l}}}$和 ${\varphi _{\rm{i}}}$分别为水和冰的体积分数.

Matala^[17]提出在冻结和融化过程中临界孔隙半径的估算方法：

(3) $\begin{array}{*{20}{c}} {{R_{{\rm{pf}}}} = 0.584 + 0.005\;2\theta - {{63.46}}/{\theta }},&{\mathop T\limits^. < 0} . \end{array}$

(4) $\begin{array}{*{20}{c}} {{R_{{\rm{pm}}}} = 0.757 + 0.007\;4\theta - {{33.45}}/{\theta }},&{\mathop T\limits^. > 0} . \end{array}$

式中： ${R_{{\rm{pf}}}}$和 ${R_{{\rm{pm}}}}$分别为冻结和融化过程中临界孔隙半径； $\theta = \left| {T - {T_{\rm{m}}}} \right|$，其中 ${T_{\rm{m}}}$为水的正常融点，273.15 K， $T$为水的温度.

考虑液体超压，孔隙水压力的计算公式^[3]为

(5) $p = {p_{{\rm{atm}}}} + \frac{{{\varphi_{\rm{l}}}\left( {1 - {{{\rho _{\rm{i}}}} / {{\rho _{\rm{l}}}}}} \right)}}{{ {{b / {\left( {n{K_{\rm{0}}}} \right) - {1 / {{K_{\rm{m}}} + {{{\varphi_{\rm{i}}}} / {{K_{\rm{i}}}}}}}}} + {{{\varphi_{\rm{l}}}} / {{K_{\rm{l}}}}}} }}.$

式中： $p$为孔隙水压力， ${p_{{\rm{atm}}}}$为大气压， ${\rho _{\rm{l}}}$和 ${\rho _{\rm{i}}}$分别为水和冰的密度， ${K_{\rm{0}}}$、 ${K_{\rm{m}}}$、 ${K_{\rm{i}}}$和 ${K_{\rm{l}}}$分别为多孔材料、砂浆、冰和水的体积模量.

根据拉普拉斯方程，结晶压力可以表示为^[1]

(6) $\begin{array}{*{20}{c}} {{p_{\rm{i}}} = p + {{f{\gamma _{{\rm{il}}}}}}/{{{R_{\rm{p}}}}}};&{T < {T_{{\rm{tran}}}}} \end{array}.$

式中： ${p_{\rm{i}}}$为结晶压力； $f$为孔隙形状因子，混凝土内部孔隙分布情况可以由压汞法（MIP）进行测量； ${T_{{\rm{tran}}}}$为相变点温度^[1]，273 K.

参考热力学约束平均理论^[18]可知，假设相平均压力等于它们的表面平均压力，可得

(7) ${p^{\rm{*}}} = {p_{\rm{i}}}{\varphi_{\rm{i}}} + p{\varphi_{\rm{l}}}.$

式中： ${p^{\rm{*}}}$为孔隙压力.

混凝土多孔体系平衡方程为

(8) $\nabla {{\sigma }} + {{F}}={{0}}.$

式中：σ为结构所受应力； ${{F}}$为重力，可以忽略.

考虑冻融循环过程中混凝土产生的损伤，引入损伤因子，损伤因子的公式可以由冀晓东等^[19]提出的损伤因子计算公式进行计算：

(9) $d = 1 - {\left( {1 - pN} \right)^q}.$

式中： $d$为损伤因子； $N$为冻融循环次数； $p$和 $q$为宏观材料参数，可由实验确定.

引入热应变 ${\varepsilon _{{\rm{th}}}} = {\alpha _{\rm{l}}}\Delta T$，其中 ${\alpha _{\rm{l}}}$为线膨胀系数，∆T= T − T_ref，T_ref为参考温度，因此有效应力矩阵为

(10) ${{{\sigma}} '} = \left( {1 - d} \right){{D}}\left( {{{\varepsilon}} - {\alpha _{\rm{l}}}\Delta T{{I}}} \right).$

式中： ${{\varepsilon}} $为应变张量，D为刚度矩阵.

最终得到力学平衡方程^[3]为

(11) $\nabla \left[ {\left( {1 - d} \right){{D}}\left( {{{\varepsilon}} - {\alpha _{\rm{l}}}\Delta T{{I}}} \right) - b{p^*}{{I}}} \right] = {\bf{0}}.$

考虑冻融过程中水冰的质量守恒，结合式(5)、(7)，得到混凝土体系冻融循环连续性方程^[3]：

(12) $ \beta \dot{p}=\nabla \cdot \left(\frac{D}{\eta }\nabla p\right)+S-b\frac{\partial \varepsilon _{\rm{v}}}{\partial t}.$

式中：D为渗透率；

(13) $\beta = \frac{{(b - n){\varphi_{\rm{l}}}}}{{{K_{\rm{M}}}}} + \frac{{n{\varphi_{\rm{l}}}}}{{{K_{\rm{l}}}}} + \frac{{n{\varphi_{\rm{i}}}}}{{{K_{\rm{i}}}}};$

$S$为压力源项，

(14) $S = \left( {1 - \frac{{{\rho _{\rm{i}}}}}{{{\rho _{\rm{l}}}}}} \right)\frac{{\partial {V_{{\rm{l}} \to {\rm{i}}}}}}{{\partial T}} + \bar \alpha - \frac{{b - n}}{{{K_{\rm{M}}}}}\frac{{\partial ({\varphi_{\rm{i}}}{p_{\rm{i}}})}}{{\partial T}} - n\frac{{{\varphi_{\rm{i}}}}}{{{K_{\rm{i}}}}}\frac{{\partial \kappa }}{{\partial T}},$

其中

(15) $\overline \alpha = \left( {b - n} \right) \overline {{\alpha _{\rm{0}}}} + n{\varphi_{\rm{l}}}{\overline{\alpha _{\rm{l}}}} + n{\varphi_{\rm{i}}}{\rm{\overline{\alpha _{\rm{i}}}}},$

(16) $\kappa = {p_{\rm{i}}} - {p_{\rm{l}}}.$

第1部分 $\left( {1 - {{{\rho _{\rm{i}}}}}/{{{\rho _{\rm{l}}}}}} \right){{\partial {V_{{\rm{l}} \to {\rm{i}}}}}}/{{\partial T}}$与结冰速率相关，第2部分 $\overline \alpha$与系统中各相的热体积膨胀而形成的压力相关（ $\overline {{\alpha _{\rm{0}}}}{\text{、}} {\overline{\alpha _{\rm{l}}}}{\text{、}}{\overline{\alpha _{\rm{i}}}}$分别为多孔材料、水和冰的体膨胀系数），第3项 $[({{b - n}})/{{{K_{\rm{M}}}}}]{{\partial ({\varphi _{\rm{i}}}{p_{\rm{i}}})}}/{{\partial T}}$与由于孔隙中结冰产生的结晶压力相关，第4项 $(n{{{\varphi _{\rm{i}}}}}/{{{K_{\rm{i}}}}}){{\partial \kappa }}/{{\partial T}}$与水冰界面的表面张力变化相关^[1].

参考Zhou等^[3]将混凝土受压损伤对渗透率的影响^[19]引入冻融循环中，得到渗透率的变化公式：

(17) $D = \varphi_{\rm{l}}^{0.5}{\left[ {1 - {{\left( {1 - \varphi_{\rm{l}}^{\rm{2}}} \right)}^{0.5}}} \right]^2}{D_0}\exp\; \left[ {{{\left( {\omega d} \right)}^\beta }} \right].$

式中： $\omega $和 $\beta $为恒定系数，由Picandet等^[20]通过物理实验得出，分别为11.3和1.64；D₀为初始渗透率.

考虑平流的影响^[3]，热传导方程可以写为

(18) $ \bar \rho \bar c\frac{{\partial T}}{{\partial t}} + {\rho _{\rm{l}}}{c_{\rm{l}}}\left( {{{v}}\nabla T} \right) = \nabla \cdot\left( {\mathop \lambda \limits^ - \nabla T} \right) + L\frac{\partial }{{\partial t}}\left( {n{\rho _{\rm{i}}}{\varphi_{\rm{i}}}} \right).$

式中：v为流体流动速度；L为潜热； $\bar \rho \bar c$和 $\bar \lambda $分别为^{[1, 3, 21]}

(19) $\bar \rho \bar c = \left( {1 - n} \right){\rho _{\rm{m}}}{c_{\rm{m}}} + n{\varphi_{\rm{l}}}{\rho _{\rm{l}}}{c_{\rm{l}}} + n{\varphi_{\rm{i}}}{\rho _{\rm{i}}}{c_{\rm{i}}},$

(20) $\bar \lambda = \frac{{n{\varphi_{\rm{l}}}{\lambda _{\rm{l}}} + n{\varphi_{\rm{i}}}{\lambda _{\rm{i}}} + {\lambda _{\rm{m}}}}}{{n{\varphi_{\rm{l}}} + n{\varphi_{\rm{i}}} + 1}},$

其中 ${\rho _{\rm{m}}}$为砂浆的密度， ${c_{\rm{m}}}$、 ${c_{\rm{l}}}$和 ${c_{\rm{i}}}$分别为砂浆、水和冰的比热容， $n$为孔隙率， ${\lambda _{\rm{m}}}$、 ${\lambda _{\rm{l}}}$和 ${\lambda _{\rm{i}}}$分别为砂浆、水和冰的热导率.

ITZ是骨料与砂浆之间的界面，线膨胀系数可以视为与砂浆相同. 材料均忽略非线性. 具体参数的确定可以参考Zhou等^[12]的计算过程，Zhou等^[12]研究分析得到的结果表明该方法是可行的.

在力学性能方面，Zhou等^[22]通过霍普金森杆层裂试验，参考试验数据^[23]可得弹性模量力学参数. ITZ的弹性模量取为砂浆的85%，渗透率取为砂浆的3倍，泊松比与砂浆相同^[3].

选取水灰质量比为m_w/m_c=0.4的水泥砂浆，相关参数由曾强等^[2] 的物理实验和段安^[1]的数值模拟参数得出. 具体的计算参数如表1、2所示. 表中，ρ为密度，E为弹性模量，K为体积模量，μ为泊松比，n为孔隙率，b为比奥系数，λ为热导率，c为比热容，K₀为多孔骨架体积模量，η为水的动力黏度，κ为水-冰界面张力。

模拟初始为零应变状态，温度为5 °C，孔隙水压力为0 MPa. 在边界条件方面，力学平衡中，将对称面和底面法向位移设置为0 m，其他面为自由面. 在连续性方程对应的渗流场中，将对称面设置为无流动，其他面的孔隙水压力设置为0 MPa. 在热传导方程对应的热学场中，将对称面设置为热绝缘，其他面温度设置为环境温度. 采用第一类边界条件，线性降温，表面温度的变化公式如下。

(21) $\begin{array}{*{20}{c}} {{T_{\rm{b}}} = {T_1} - t/1\;000 + 36 a};&{\mathop T\limits^. } < 0. \end{array}$

(22) $\begin{array}{*{20}{c}} {{T_{\rm{b}}} = {T_2} + 3 t/1\;000 - 108 a};&{\mathop T\limits^. } > 0. \end{array} $

式中： ${T_1}$=278.15 K， ${T_2}$=170.15 K， ${T_{\rm{b}}}$为环境温度， $a$为冻融循环周期，t为时间.

环境温度变化如图2所示.

利用1章生成的有限元网格模型，考虑多物理场求解解析解的难度过大，故采取有限元方法进行离散数值求解. 采用大型有限元软件Comsol Multiphysics进行多物理场耦合求解，选择固体传热、达西定律和固体力学模块，采用瞬态求解器，时间步为60 s，收敛的相对容差为0.01.

将该研究修正模型与Duan等^[11]提出的数值模型进行对比验证分析. 引入冻结应变^[11]，可以表示为

(23) ${\varepsilon _{\rm{f}}} = {\varepsilon _{{z}}} - {\alpha _{\rm{l}}}\Delta T. $

式中： ${\varepsilon _{\rm{f}}}$为冻结应变， ${\varepsilon _{{z}}}$为z方向应变.

模型采用均质试件. 边界温度以4 °C/h的降温速率从2 °C线性降温至−21 °C. 选取试件表面点，选择温度为0 °C至−20 °C，计算所得的冻结应变见图3.

修正模型计算结果与试验所得结果表现出类似的趋势. 通过试验所得的最大冻结应变为6.40×10⁻⁴，Duan等^[11]提出的模型计算所得的最大冻结应变为5.62 ×10⁻⁴，利用该修正模型得到的最大冻结应变为6.12×10⁻⁴. 原因是修正模型考虑了损伤效应，与Duan等^[11]的模型相比，该修正模型的研究结果更准确. 利用修正模型可以更加准确地预测材料的应变演化情况. 然而，模型结果与试验结果存在一定的差别. MIP测量孔隙的试验必须进行干燥处理，该过程将导致孔隙的测量不准确. 利用MIP法只能测量到孔径为15 nm以上的孔隙，这可能导致孔隙的统计不完全. 为了简化计算，假设试件为完全饱和并且孔隙为圆柱形，孔隙中仅存冰和水. 在试验中，试件内部孔隙形状不规则且不可避免地存在非饱和孔隙，这些差别将影响计算结果. 虽然本文考虑了损伤效应，但是基于弹性假设，没有考虑残余变形部分，这一部分是混凝土冻胀变形的重要部分^{[3, 11]}. 总体来说，与Duan等^[11]的模型相比，修正模型拥有更好的预测变形能力.

在数值试验中，研究砂浆和界面过渡区（ITZ）渗透率、界面过渡区（ITZ）线膨胀系数、骨料体积分数等因素的影响. 这些因素可能对混凝土在冻融情况下的热力学响应产生一定影响.

冻融过程中将产生水的迁移，这将对混凝土内部孔隙压力及应力产生重要影响，因此分析渗透率对混凝土在冻融环境下的热力学响应有重要意义. 如图4、5所示分别为砂浆渗透率为3.675×10⁻²¹、3.675×10⁻²⁰和3.675×10⁻¹⁹ m²时的应力分布情况（按照参数确定相关结果，ITZ渗透率取砂浆渗透率的3倍）. 选取第1个冻融循环结束的时刻，图1所示10 h以及第3个循环降温过程的25 h时刻，绘制应力云图. 可以看出，无论是在冻结还是在融化的过程中，随着渗透率的增大，应力逐渐减小. 这是因为根据达西定律可知，压力与渗透率成反比^[1]，随着渗透率的增大，水分的迁移难度减小，因此压力梯度减小，所求应力减小，这一结论与Powers等^[3-4]的理论一致. 此外，与二维模型^[3]相比，随着渗透率的增大，骨料内部和外部第一主应力相差更明显. 可以看出，随着渗透率的增大，试件中出现骨料内部为负值而骨料外部为正值的情况. 这可能是由于骨料的材料属性跟其他组分差别较大，如骨料不透水且线膨胀系数较外部小，随着渗透率的增大，第一主应力减小，骨料外部冻胀效应相较于热膨胀产生的温度应变无绝对优势，因此温度应变在较高渗透率的情况下发挥了较明显的作用. 这导致骨料内、外的形变程度不同，骨料内部受压，第一主应力为负值，且在骨料附近区域应力出现不连续现象.

由图6~8可知，在每个循环的约前1.5 h左右，应力较小，这是因为此段为在相变点之前，孔隙水还未结冰，结晶压力为0，存在较小的静水压力，为负值（以压为正）. 之后大量的大孔隙中的水结冰，孔隙水产生迁移，孔隙水压力上升较快。从达西定律可知，渗透率与水压力成反比，渗透率越大，水分迁移能力越强，液体压力梯度越小，导致水压力越小. 根据拉普拉斯方程可知，结晶压力越小，孔隙压力作用在试件上，使得试件受拉，由此计算所得的应力变小. 结晶压力和应力随着渗透率的改变而产生的变化规律与孔隙水压力变化规律类似. 可见，渗透率对应力有重大影响.

ITZ的孔隙率较高，孔隙孔径较大，未水化的水泥较少. 分析此处的线膨胀系数对整体应力的分布有一定意义.

如图9(a)~(c)、10(a)~(c)所示分别为ITZ原线膨胀系数α_mortar、3倍、5倍时的应力分布情况. 图中，α_itz为线膨胀系数。选取第3个循环降温过程的15 h时刻以及第5个循环结束的50 h时刻，绘制应力云图. 从图9可以看出，在不同的ITZ线膨胀系数下，混凝土应力分布情况相似. 当线膨胀系数增大时，试样骨料与砂浆的界面出现的应力略有增大，这是因为线膨胀系数越大，热膨胀效应越明显，应力增大，说明线膨胀系数增大有可能造成混凝土在该薄弱区域产生损伤.

由图11可知，在不同的ITZ线膨胀系数下，混凝土应力分布情况相似. 在混凝土的冻融循环中，冰晶成核将产生结晶压力，热膨胀引起热应力。在温度变化的过程中，结晶压力产生的冻胀效应远远超过因热应力而引起的收缩和膨胀效应，因此线膨胀系数对应力的影响相对较小.

如图12(a)~(c)、13(a)~(c)所示分别为均质试样，骨料体积分数为10%和20%的试样的温度分布情况. 选取第1个循环降温过程的5 h时刻以及第4个循环结束的40 h时刻，绘制温度云图及等温面图. 由于试件体积较小，在降温阶段和循环结束时刻，不同体积分数的试件温度差别不大. 与二维细观分析结果^[3]相比，骨料的存在使得温度在混凝土内部骨料附近的分布产生更明显的弯折现象，即温度分布相对不均匀，这是因为骨料导热性能与砂浆和ITZ有一定差别，传热中产生了温度分布不均匀的情况. 这对于控制温度影响混凝土内部薄弱面开裂有一定启示.

如图14、15所示，在第一主应力云图中，均质试样的主应力分布相对其他2种情况而言更均匀，骨料对应力分布产生了一定的影响. 无论是冻结还是融化过程，可以看出骨料内部与外部砂浆和ITZ的应力差别较大，这是由于骨料为不透水介质，骨料内部不受结晶压力和水压力. 由参数确定结果可知，在细观参数方面，骨料的细观计算参数与其他组分相比，差别较大. 骨料内部均承受较大应力而外部应力相对降低，骨料的存在导致混凝土内部应力分布产生不均匀现象，除去骨料内部应力部分，混凝土应力整体呈现降低的趋势.

从图16、17可知，随着骨料体积分数的增加，试件所受第一主应力峰值从72 MPa降至66 MPa，最后降至46 MPa. 这是因为随着骨料体积分数的增大，相对减少了砂浆的体积分数，总渗透率降低，第一主应力减小，这说明骨料体积分数越高，混凝土的抗冻融性能越好. 在传热计算中，由等温面结果可知，骨料的存在显著改变了温度的分布情况，但是由于试件尺寸较小，随着骨料体积分数的增加，温度没有很明显的变化，因此温度变化历程曲线基本一致.

（1）渗透率对混凝土抗冻融性能有着重大的影响. 随着渗透率的增大，孔隙水压力、结晶压力和第一主应力减小，且骨料内部和外部第一主应力相差更明显；当渗透率很大时，出现骨料受压而砂浆受拉的情况. 在工程应用中，选择合适强度的水泥砂浆可以降低渗透率，从而降低整体应力.

（2）ITZ线膨胀系数对混凝土抗冻融性能的影响不大. 当ITZ线膨胀系数增大时，骨料与砂浆的界面的应力略有增大，线膨胀系数增大有可能造成混凝土在该薄弱区域产生损伤. ITZ线膨胀系数对试样内部总体应力几乎无影响.

（3）骨料体积分数对混凝土的抗冻融性能有较大的影响. 随着骨料体积分数的增加，第一主应力减小. 骨料可以显著提升混凝土的整体抗冻融性能.