湿筛混凝土动态受压力学性能及破坏形态细观数值模拟

张清芳，洪鹤轩，沈 璐,2

(1.大连海洋大学海洋与土木工程学院，大连 116023;2.大连海洋大学辽宁省海洋牧场工程技术研究中心，大连 116023)

0 引 言

水工大体积结构物一般采用全级配混凝土，其最大骨料粒径可达150 mm，这就给试件制作及开展物理试验带来很大的困难。目前国内外普遍采用湿筛法来研究水工大骨料混凝土的动、静态力学性能，即将粒径大于40 mm的粗骨料筛除，制作成湿筛混凝土标准试件来进行试验研究[1]。由于最大骨料粒径、骨料占比及试件尺寸的差异，湿筛混凝土的力学性能与大骨料混凝土及普通混凝土相比均有一定的差别，需要专门进行研究。一些学者[2-4]针对湿筛混凝土与全级配混凝土之间力学性能的差异开展试验研究，研究表明骨料含量、骨料级配、特殊骨料对两者力学性能差异有着较为明显的影响。Serra等[5-6]考虑尺寸效应和最大骨料粒径因素，提出了一种基于湿筛混凝土强度预测大坝混凝土强度的方法，并在Baixo Sabor大坝中进行了验证。还有一些学者[7-11]较为系统地对湿筛混凝土进行了单轴、双轴、三轴动静态力学性能试验，得出了多轴动态强度准则及应力-应变曲线。通过物理试验得到的混凝土力学性能结果一般是比较可靠的，但是也存在着一些问题，比如破坏形态及裂缝数量难以进行量化分析，这直接影响到湿筛效应及动态效应的机理研究。目前，针对湿筛混凝土力学性能及破坏形态的数值模拟研究相对较少，因此需要开展相应的机理研究和数值模拟作为补充。

混凝土细观数值模拟是将混凝土看作是由骨料、砂浆和界面组成的三相复合材料，根据混凝土的骨料级配建立随机骨料数值模型，然后将骨料、砂浆和界面分别赋予细观参数后进行动、静态力学加载试验，从而达到获取混凝土宏观力学性能的目的[12-13]。目前，采用细观数值模拟方法研究混凝土力学性能的成果较多，并且普遍认为细观模型的建模方式及细观参数的设置是开展混凝土细观分析的关键[14-15]。与物理试验相比，采用颗粒流离散元软件PFC2D开展细观数值模拟可以定量地分析混凝土开裂情况，也可以通过分析内部接触力来探究混凝土非均质特性对动态效应的影响[16]。

本文利用颗粒流离散元软件PFC2D，对照现有文献中的物理试验，建立湿筛二级配混凝土的随机骨料模型，然后以应变速率为10-5s-1作为拟静态，根据文献中的静态压缩加载物理试验结果标定出砂浆颗粒之间、骨料颗粒之间以及砂浆颗粒与骨料颗粒接触界面之间的细观参数，进而形成湿筛二级配混凝土细观数值试件，然后对数值试件开展地震作用对应的10-4s-1、10-3s-1、10-2s-1应变速率下的单轴动态压缩加载试验，将数值模拟得出的动态力学性能指标与文献中的物理试验结果进行对比分析，最后对破坏形态、裂缝数量及动态强度提高机理开展进一步研究。

1 理论基础

1.1 颗粒离散元基本原理

图1 颗粒离散元中的接触模型Fig.1 Contact model in discrete elements of particles

颗粒离散元是Cundall等[17]于1979年提出的用于描述颗粒之间运动规律的数值计算方法，PFC2D就是基于颗粒离散元理论开发的软件平台。颗粒离散元理论将单个颗粒视为刚体，通过墙体边界对颗粒组成的模型进行加载，模型的变形及其他力学特性通过颗粒之间的接触关系来体现，其中用于表征颗粒接触关系的力学模型非常关键。颗粒离散元中的接触模型如图1所示，基本的接触模型包括刚度模型、滑移模型和粘结模型。

(1)刚度模型，即两个颗粒之间的接触力与其间的相对位移之间的关系，可以表示为：

Fn=knUn

(1)

ΔFs=-ksΔUs

(2)

式中：Fn为法向力；kn为法向刚度；Un为法向位移；ΔFs为切向力增量;ks为切向刚度;ΔUs为剪切位移增量。

(2)滑移模型，即两个颗粒之间的切向接触力与法向接触力之间的关系，f为摩擦系数。

(3)粘结模型，用以描述接触的极限法向力和极限切向力，包括接触粘结模型和平行粘结模型。

1.2 细观参数反演原理

颗粒流离散元软件PFC2D在进行细观数值模拟方面具有一定的优势，其便捷之处在于无需提前获取材料的本构模型，因此细观参数的取值将直接影响数值模拟的结果。目前，国内外普遍采用的确定细观参数的方法是试算法，此法得到的细观参数虽然能够使计算结果与物理试验数据基本一致，但是往往无法定量地对其进行准确度评价。基于响应面法的细观参数反演[18-20]，就是通过构造高维超曲面响应面函数，定量描述及优化细观参数的一种实用方法。响应面函数通常选择不含交叉项的二次多项式，待反演的混凝土细观参数与应力之间的响应面函数可以写成式(3)～(5)的形式。

(3)

(4)

(5)

当k=1时，将式(3)展开得到式(6)～(12)。

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

2 混凝土细观数值模型与动态加载

2.1 湿筛二级配混凝土随机骨料模型的建立

水工混凝土的配合比按照大坝变水位区常用的三级配混凝土进行设计，为了将数值模拟结果与物理试验结果做对比分析，本文配合比采用文献[11]中的配合比，如表1所示。根据配合比可知，1 m3混凝土中含有1 475 kg的粗骨料，按照花岗岩的密度为2 700 kg/m3进行计算，此时粗骨料的体积分数为54.6%，其中粒径范围为5～20 mm、20～40 mm、40～80 mm的三种骨料的体积占比分别为16.4%、16.4%、21.8%。湿筛二级配混凝土就是将拌和好的三级配混凝土拌合物筛除超过40 mm粒径骨料，再将筛下的混凝土拌合物翻拌均匀，最后进行浇筑振捣形成的[1]。从其制作过程可知，总体积为1 m3的三级配混凝土制成湿筛二级配混凝土后，由于筛除了大粒径的骨料，其体积降为0.78 m3，粒径范围为5～20 mm、20～40 mm的两种粗骨料的体积分数均变为21.0%。湿筛混凝土立方体试件边长取150 mm，考虑到计算效率，本文建立平面模型进行分析，并采用基于Monte-Carlo法提出的经典“取-放”方法[13]，建立湿筛二级配混凝土试件的随机骨料模型，如图2所示，其中，粒径范围为20～40 mm的粗骨料用7个直径为30 mm的圆形代替，粒径范围为5～20 mm的粗骨料用42个直径为12 mm的圆形代替，以保证各直径圆形的体积分数与相应粒径范围的粗骨料体积分数相等。

表1 混凝土的配合比[11]Table 1 Mix ratio of concrete[11] /(kg·m-3)

2.2 湿筛二级配混凝土试件数值模型的建立

大量试验现象表明，在物理试验中，混凝土受压破坏时不仅仅只有砂浆会开裂，也会出现骨料断裂的现象，特别是对于动态加载，加载速率较高时，骨料断裂的比例会更大[7-11]。现有的基于PFC2D软件的混凝土细观数值模拟普遍采用单个刚体颗粒来模拟粗骨料，无法实现骨料开裂现象的模拟。本文采用的建模方式为在PFC2D软件的计算区域内建立边长150 mm的正方形墙体作为试件边界，然后在该范围内随机投放不规则排列的颗粒并进行压实，考虑到计算精度和运行速度，本文颗粒半径在1～1.5 mm范围内。根据图2所示的随机骨料位置将颗粒分组，具体做法为将圆形以内和圆形以外两个部分分为两个图集，图形集外部颗粒命名为sand，30 mm直径圆形图形集内部颗粒命名为agg 1，12 mm直径圆形图形集内部颗粒命名为agg 2。颗粒之间的接触采用平行黏结模型。建立的湿筛二级配混凝土试件细观数值模型如图3所示。

图2 湿筛二级配混凝土试件的随机骨料模型Fig.2 Random aggregate model of wet-screened two-graded concrete specimen

图3 湿筛二级配混凝土试件细观数值模型Fig.3 Meso-numerical model of wet-screened two-graded concrete specimen

图4 物理试验与数值模拟的应力-应变曲线Fig.4 Stress-strain curves of physical experiment and numerical simulation

2.3 细观参数的确定

本文将颗粒间的接触参数分为三部分，一部分是砂浆颗粒间的接触参数，一部分是粗骨料颗粒间的接触参数，还有一部分是砂浆颗粒与粗骨料颗粒间的接触参数。采用PFC2D软件开展混凝土细观数值模拟时，细观参数并无实际的物理含义。通过文献[20]中的方法，根据文献[11]中应变速率为10-5s-1时的湿筛二级配混凝土单轴压缩加载试验数据，得到各组颗粒之间的接触参数，如表2所示。该组参数下的单轴静态受压应力-应变曲线与物理试验结果吻合较好，物理试验与数值模拟的应力-应变(σ-ε)曲线如图4所示，可认为该组细观参数基本满足要求。

表2 各分组颗粒之间的接触细观参数Table 2 Meso-parameters of contact between grouped particles

2.4 数值试件动态加载

为了模拟物理试验中减摩垫片的作用，数值试件的上下两面均设置了摩擦系数为0的刚性墙，通过PFC2D伺服机制控制上层墙体移动，使得墙体的应力值缓慢接近目标值，达到压缩混凝土的效果。试验分为四组，采用相同的试件和加载方式，分别以0.02 m/s、0.3 m/s、2.2 m/s、3.8 m/s的速率移动上层墙体，经过换算，得到四种加载速率下的实际应变速率,分别为5.35×10-5s-1、8.03×10-4s-1、5.89×10-3s-1、1.02×10-2s-1，其中5.35×10-5s-1为拟静态加载。加载完成后，导出应力-应变曲线、试件的破坏形态、裂缝分布情况以及内部接触力等试验结果。

3 细观数值模拟结果与讨论

3.1 应力-应变曲线

图5 不同加载速率下数值试件的应力-应变曲线Fig.5 Stress-strain curves of numerical specimens with different loading rates

图5为不同加载速率下的混凝土数值试件单轴动态压缩应力-应变曲线。从图中可以看出，湿筛二级配混凝土在不同加载速率下应力-应变曲线较为相似，均经历线性上升段、非线性上升段、应力达到峰值后的下降阶段，并且下降段近似平行。随着应变速率的增加，动态加载下的峰值应力、峰值应力处应变均较静态加载时有所增大，峰值应力处应变的增大幅度很小。这表明湿筛混凝土的强度具有明显的应变率效应，而在低应变速率下，湿筛混凝土的延性变化不明显。与文献[11]的物理试验结果比较可知，数值模拟与物理试验应力-应变曲线上升段吻合较好，两者的下降段斜率有一定差异，这主要由于物理试验受试验机刚度影响，特别是动态加载对试验机的刚度要求更高，而数值模拟不存在试验机刚度不足的问题。本文的峰值应力、峰值应力处应变与文献[11]的物理试验数据相比较，如表3所示。当加载速率为10-4s-1、10-3s-1、10-2s-1时，其峰值应力分别较静态增加7.3%、23.3%、37.9%，这与文献[11]中的物理试验结果14.7%、22.8%、33.5%比较接近。从表3中可以看出:数值模拟与物理试验在动态峰值应力方面误差较小，基本在±5%左右；在峰值处应变方面，数值模拟的结果总体偏大，但是最大误差不超过8%，并且随应变速率增大的趋势与文献[11]一致。

表3 数值模拟结果与文献[11]物理试验结果对比Table 3 Comparison between numerical simulation results and physical test results in literature [11]

3.2 破坏形态及裂缝分布情况

在加载至峰值应力时刻，湿筛二级配混凝土在四种应变速率下的破坏形态和裂缝分布情况分别如图6、图7所示，文献[11]中物理试验下的试件破坏形态如图8所示。从图6、图7中可以看出，试件破坏时都会出现纵向贯穿试件的主裂缝，随着应变速率的提高，主裂缝和微裂缝的数量也在增多，这与图8物理试验的现象是一致的。数值模拟可以统计出微裂缝的数量，如图7所示，在应变速率为10-4s-1、10-3s-1、10-2s-1时，其微裂缝数量较静态分别增加48%、78%、149%。综合考虑峰值应力的动态特性可知，随着应变速率的增加，裂缝数量的增长与峰值应力的增长基本保持一致，裂缝数量增长率分别为数值模拟中峰值应力增长率的6.6倍、3.3倍、3.9倍，与物理试验数据相比较则分别为3.3倍、3.4倍、4.4倍，除第一个数据离散性较大以外，其余均在3.3～4.4倍范围内，平均值为4.2倍。此外，从图7中也可以看出，静态加载时，裂缝分布不均匀，基本集中在试件局部,随着应变速率的增加，裂缝的分布逐渐趋于均匀。

图6 不同加载速率下数值试件的破坏形态Fig.6 Damage patterns of numerical specimens with different loading rates

注：DFN表示离散裂隙网络；Fractures(1 254)表示裂隙数量为1 254；Facets为图例，表示有颜色的部分为平面裂纹。 图7 加载完成时数值试件的内部裂缝的分布情况Fig.7 Distribution of internal cracks of numerical specimens at the completion of loading

从图6、图7中还可以看出，应变速率分别为10-5s-1、10-4s-1、10-3s-1、10-2s-1时，直径20 mm以上的大粒径骨料开裂数量分别为3个、4个、5个、7个，随着应变速率的增加，大粒径粗骨料开裂的数量也明显增多。许多学者[1,7-11]在物理试验中也观察到粗骨料断裂数量随加载速率的增大而增多的结论，这与本文细观数值模拟的结果非常吻合，与物理试验相比,采用数值模拟能够更加方便地进行定量统计。

图8 物理试验时试件的破坏状态[11]Fig.8 Failure mode of specimens in physical test[11]

3.3 内部接触力

混凝土峰值应力时刻，不同加载速率下内部颗粒接触力分布如图9所示。从图中可以看出，在不同加载速率下，全截面接触力大小总体比较均匀，但是个别部位接触力较大，呈现出局部不均匀的现象，应变速率为10-2s-1时接触力低于3.0×106N(蓝色)箭头的密集程度高于10-5s-1。为了定量统计内部接触力的不均匀程度，本文将所有颗粒之间的接触力由大到小排列，将前10%的接触力平均值与所有接触力平均值的比值作为不均匀性评价指标，四种应变速率下的指标数值分别为19.746、17.183、15.618、15.063，动态加载情况明显小于静态加载，表明加载速率的提高使得混凝土内部应力分布趋于均匀，不容易出现应力集中现象，这也是造成湿筛混凝土动态峰值应力较静态提高的原因之一。

图9 不同加载速率下混凝土峰值应力时刻接触力分布图Fig.9 Contact force distribution at the moment of peak stress of concrete under different loading rates

4 结 论

本文通过对湿筛二级配混凝土在10-5s-1、10-4s-1、10-3s-1、10-2s-1四种应变速率下的受压力学特性及破坏形态进行细观数值模拟，得出以下结论：

(1)根据本文的建模方式和细观参数开展的动态加载数值模拟，其各项力学性能指标与物理试验数据的误差在8%以内，破坏形态与物理试验现象吻合较好，该细观数值模拟方法可用于混凝土材料动态力学性能分析。

(2)随着应变速率的提高，湿筛二级配混凝土的峰值应力不断增大，裂缝数量也逐渐增多；由应变速率10-5s-1增加到10-2s-1，其峰值应力增加7.3%～37.9%，裂缝数量的增长率平均为峰值应力的4.2倍；

(3)从混凝土的非均质特性角度分析，湿筛混凝土内部接触力的非均匀程度对混凝土宏观动态力学性能将产生一定影响，其内部接触力分布的不均匀性随着应变速率的提高而有所降低，这是造成其动态强度增长的原因之一。