22

(1.四川大学 水利水电学院，成都 610065； 2.四川大学 水力学与山区河流开发保护国家重点实验室，成都 610065)

1 研究背景

无砂混凝土作为典型的多孔混凝土材料，是由粗骨料、水泥、水及少量外加剂拌合而成的一种轻质混凝土[1]。较高的孔隙率和良好的透气、透水性使其在护岸护坡工程和城市生态透水地基工程方面具有良好的应用[2]。无砂混凝土内部不含细骨料，其自身强度来源于粗骨料在水泥浆的作用下形成大量的点粘结和微小面粘结，与常规混凝土相比，其发达的孔隙结构使得破坏损伤过程更为独特且复杂[3]。Bazant和OH[4]基于细观尺度下的损伤原理提出了微平面模型，该模型认为混凝土的裂纹可以向任意方向扩展。杨延毅和周维垣[5]将断裂力学模型和损伤力学模型进行了结合，提出了混凝土材料的断裂准则，能成功地用于一些实际工程。徐阳[6]通过试验分析与三维细观数值模型的建立，对多孔生态混凝土损伤、裂缝发展规律进行了探讨，结果表明：应力通过骨料纵向传播，上层存在孔隙的骨料区产生的拉应力引发初始裂缝。刘光廷和王宗敏[7]借助蒙特卡罗法利用有限元技术模拟了混凝土单边裂纹的受拉情况，通过观察损伤、断裂破坏的整个过程，可知其模拟结果与实际相接近。

离散元法由Cundall等[8-10]首先提出。PFC颗粒流方法作为离散单元法的一种，多用于研究散粒介质材料的细观力学特性，它弥补了传统连续介质力学模型在计算非连续介质及大变形问题上的不足，在模拟无砂混凝土这种结构离散型较大的材料时具有一定优势。王云飞等[11]利用PFC2D软件建立不同形态粗骨料混凝土数值模型，研究粗骨料含量和形态对强度和损伤演化的影响。宿辉等[12]采用PFC2D创建了含不规则骨料的生态混凝土数值模型，利用控制变量法针对不同细观力学参数对其双轴压缩破坏的影响进行了研究，结果表明利用PFC创建的不规则块体结构能够真实地模拟生态混凝土。张正珺等[13]采用PFC2D对混凝土试件单轴压缩破坏全过程进行模拟，结果表明，基于PFC的混凝土材料破坏模拟具有显著的优点，可以真实地模拟裂纹的生成、扩展以及破坏的全过程。

以往基于PFC的混凝土数值模拟多从二维角度进行建模分析，具有一定的局限性。鉴于此，本文使用三维颗粒流软件PFC3D对无砂混凝土单轴受压试验进行数值模拟，在对典型卵石骨料仿真建模的基础上建立无砂混凝土三维数值模型，从细观角度研究无砂混凝土强度变化及应力-应变曲线特性，分析破坏损伤演化过程及分布特征，以期为无砂混凝土的力学及损伤机理研究提供参考。

2 数值模型建立

2.1 骨料模型生成

无砂混凝土中粗骨料形态差异较大，传统的计算中把骨料全部简化成球型并不能真实地反映实际情况，因此，在骨料形态数目统计的基础上选取椭球状和盘状2种较为典型的骨料形状进行模拟试验，利用clump命令建立了骨料颗粒模型，盘状模型由5个球型颗粒构建而成，椭球状模型由3个球型颗粒构建而成。在加载过程中颗粒模型自身不会发生破坏，破坏只发生在颗粒之间生成的接触粘结。2种形状骨料模型与实际骨料对比如图1所示。

图1 骨料模型与实际骨料对比Fig.1 Comparison of numerical aggregate model with actual aggregate

2.2 细观力学参数选取

PFC3D中颗粒间采用特定的粘结模型进行粘结，如接触刚度模型、线性模型和线性平行粘结模型。其中，采用线性平行粘结模型生成的粘结可以同时传递力和力矩，当法向接触力达到法向粘结强度σn或切向接触力达到切向粘结强度σs时，颗粒粘结发生断裂形成微损伤，进而发展形成破裂面使试样破坏，这与实际混凝土材料的损伤发展相符[14]。因此，模拟中颗粒间接触采用线性平行粘结模型，颗粒与墙体接触采用线性接触。接触模型确定后对计算所需的细观力学参数进行确定，各接触模型的力学参数有法向刚度kn、切向刚度ks、法向粘结强度σn、切向粘结强度σs；颗粒力学参数有颗粒密度、颗粒摩擦系数μ。这些参数除颗粒密度外均无法直接通过材料试验获得，因此，在进行数值计算前，需要对模型细观力学参数进行标定。

由于模拟中所采用的颗粒并不是大小、形状都完全相同的球体，颗粒间的接触刚度kn和ks设置为统一值并不合理，因此，采用设置有效模量emod来对不同接触状态的粘结自动计算其接触刚度。对于有效模量、法向粘结强度、切向粘结强度的标定，参考Fakhimi等[15]提出的量纲分析法及康政等[16]提出的数值仿真反演计算思路，结合无砂混凝土实际试验强度值进行数值模拟试算确定。对于颗粒摩擦系数μ，Potyondy和Cundall[17]认为，摩擦模型是在颗粒粘结发生断裂破坏后才被激活，因此摩擦系数μ对数值仿真试件的峰值强度的影响较小，仅对破坏后卸载段略有影响。目前尚未有明确的宏观量与之对应，本文通过对比实际试验应力-应变曲线卸载段情况，采用试算法来确定合适的μ值，各力学参数最终采用如表1所示的数值。

表1 模型参数设置Table 1 Settings of model parameters

2.3 混凝土模型生成

在标准混凝土模型立方体范围内以指定的孔隙率生成颗粒模型，颗粒体积范围为2～4 cm3，大小服从Gauss分布，2种形状颗粒总体积比为1∶1。生成的无砂混凝土模型及粘结模型如图2所示。

图2 无砂混凝土数值模型Fig.2 Numerical model of no-fines concrete

3 数值模拟分析

3.1 孔隙率对抗压强度及应力-应变曲线的影响分析

在试块模型建立后，通过给墙体赋予指向试块的速度来实现单轴抗压强度试验的模拟。7个不同孔隙率混凝土模型的单轴模拟试验应力-应变曲线如图3所示。

图3 单轴抗压强度试验中试样的应力-应变曲线Fig.3 Stress-strain curves of samples in uniaxial compression test

图3中，各孔隙率试块应力-应变曲线上升段都接近于直线，在达到峰值强度后曲线迅速下降，孔隙率0.32以下的模型在下降段下降速率较快。同时，试块峰值抗压强度及弹性模量随孔隙率的上升逐渐降低，抗压强度从孔隙率0.26时的24.9 MPa降低到孔隙率0.38时的7.97 MPa。将各孔隙率下模型宏观基本参数统计成表2，并绘制各参数相应于孔隙率的变化趋势如图4。

表2 混凝土数值模型宏观参数及抗压强度Table 2 Macroscopic parameters and compressive strength of numerical concrete model

图4 数值模型参数与孔隙率关系Fig.4 Relationship of numerical model parameters against porosity

由表2及图4可知，随着孔隙率的上升，模型峰值强度、组成颗粒数、颗粒间粘结数均呈下降趋势，各参数与孔隙率之间均有较好的线性相关性，峰值强度与孔隙率之间的相关趋势与试验研究结论较为一致[18]。结合模型宏观参数分析认为，随着孔隙率上升，模型本身所包含的颗粒数减少，接触粘结数量也随之下降。颗粒之间的接触粘结对应于无砂混凝土骨料之间由水泥浆胶结形成的微粘结面或粘结点，这些微粘结面和粘结点是维持无砂混凝土强度重要因素。因此，随着孔隙率的上升，颗粒(骨料)之间的接触粘结(微粘结面或粘结点)数量的下降是导致抗压强度下降的主要原因。

3.2 应力-应变全曲线分析

为了分析数值模拟得出的应力-应变曲线的合理性及无砂混凝土应力-应变特性，将试验结果与熊耀清和姚谦峰[19]给出的多孔混凝土应力-应变全曲线模型及过振海教授等[20]给出的常规混凝土应力-应变全曲线模型进行对比，将曲线进行标准化处理得到应力-应变关系曲线簇如图5所示。

图5 混凝土单轴受压应力-应变全曲线簇Fig.5 Full stress-strain curves of concrete under uniaxial compression

图5中采用无量纲坐标，横坐标为ε/εc，纵坐标为σ/σc。其中，σ为混凝土所受压应力，ε为对应的压应变，σc为混凝土受压应力峰值，εc为对应于σc的应变。

由图5可以看出：在曲线上升阶段，数值模拟结果与基于实际试验得出的多孔混凝土应力-应变全曲线模型基本一致。当ε/εc<0.5时，二者接近于线性上升；当0.5<ε/εc<1时，多孔混凝土曲线模型轻微上凸，出现一定的塑形应变。而数值模拟曲线虽在峰值强度前出现小幅波动，但总体仍接近于线性上升，体现出较强的脆性。在下降阶段，模型曲线与数值模拟曲线在峰值强度后下降较快，不同孔隙率试块数值模拟曲线下降速率具有较强的随机性。

另一方面，与多孔混凝土模型曲线及数值模拟曲线相比，常规混凝土模型曲线在上升阶段出现明显的上凸现象，说明常规混凝土随着加载的进行，塑形应变不断增加，刚度下降及软化现象较为明显。在下降阶段，常规混凝土曲线模型下降趋势较为平缓。

混凝土塑性-损伤模型表明：混凝土作为一种准脆性材料，其受压破碎是随机分布的微裂纹演化而成。在宏观现象上，混凝土的软化、永久变形及刚度退化是微裂纹逐渐演化的结果[21]。结合该模型，分析以上规律认为：相比于常规混凝土而言，多孔混凝土具有大量的内部孔隙空间。随着加载的进行，多孔混凝土内部产生的微小断裂面极易在内部高孔隙率的结构下形成贯通断裂面，这使得多孔混凝土在线弹性阶段结束后还未产生明显的软化和刚度退化现象时便达到其峰值强度，混凝土随即发生破坏卸载，体现出低强度高脆性的力学性能。无砂混凝土作为典型的多孔混凝土，其内部更大的孔隙空间使其脆性特点更加显著。

3.3 损伤断裂演化过程分析

为了直观地对模型试块内部的断裂演化情况进行跟踪，通过PFC内置的fish语言编写程序来实现接颗粒间接触断裂的实时监测。该程序的主要功能为：在数值模拟加载过程中，每当颗粒间的接触发生断裂破坏时，该程序自动记录接触断裂位置并生成圆片状断裂面来对断裂点进行标识。通过该程序对整个加载过程断裂发展进行监测，并将断裂位置分布图与实时的应力-应变曲线绘制于同一图中形成参照。断裂演化过程如图6所示。

由图6可以看出：随着加载的进行，在曲线上升阶段，模型边角及加载面处断裂逐渐产生并增多。当应力达到峰值时，模型虽未形成贯通的断裂面，但已不能继续抵抗增加的荷载，模型进入卸载阶段。应力达到峰值后，断裂点数量急剧上升形成宏观的贯通断裂面，模型完全破坏。该演化过程较好地解释了单轴受压试验中宏观贯通断裂面往往出现在峰值强度后的卸载阶段的现象。同时，图6(f)中显示了破坏后模型试块的情况，颗粒颜色异常区域为破坏发生的位置。将该模型破坏情况与实际试验中试块破坏情况(如图7)进行对比发现：数值模型易破坏情况与实际情况较为一致，断裂区域均多发生于边角及加载面处，说明PFC能够较好地模拟无砂混凝土断裂损伤演化过程。

图6 混凝土数值模型断裂损伤演化过程Fig.6 Damage evolution of numerical concrete model

图7 无砂混凝土实际破坏情况Fig.7 Actual failure of no-fines concrete

3.4 峰值强度下模型断裂信息统计及分布特征

无砂混凝土破坏过程中，主贯通面往往出现在峰值强度后的卸载阶段。然而对于混凝土结构承载能力的研究来说，在应力达到峰值强度时内部所产生的断裂才是影响其承载力的关键因素。因此，对应力与断裂数的关系进行监测，结果如图8所示，峰值强度下对应的断裂信息列于表3。

图8 应力与断裂数关系Fig.8 Relationship between stress and fracture number

由表3可知，当断裂数与模型内部初始粘结数之比超过0.032时，部分模型开始达到峰值强度；该比值达到0.086时，各模型均已达到峰值强度，体现出高脆性的特点，即当内部断裂数与模型总颗粒接触数比例不足9%时，模型试块均陆续达到其峰值承载能力而被破坏。

表3 峰值强度下断裂信息统计Table 3 Statistics of fracture information at peak intensity

图9 断裂位置分布统计Fig.9 Statistics of fracture location distribution

为了进一步定量地对混凝土模型易损伤区域进行划分，对峰值强度下断裂位置分布进行统计分析，以模型中心点为原点，将各个断裂位置距中心点的距离进行计算，并绘制断裂点到模型中心的距离分布图，如图9所示。从图9中可以看出，各孔隙率模型峰值强度下断裂数量在距离中心点位置>0.075 m的区域分布较为集中。据此对混凝土模型“易损区”进行定义，即标准混凝土试块除去内接球的边角区域为混凝土试块的“易损区”，该区域所占模型总体积的47.6%。

进一步统计该区域断裂数量与峰值强度下断裂总数量的比值，结果如表4。

表4 峰值强度下易破坏区断裂数比例Table 4 Ratio of fracture number in the vulnerable zone at peak intensity

由表4可知占模型总体积的47.6%的“易损区”包含了峰值强度下69%以上的断裂数。同时，这与实际试验中试块多在边角发生破坏的现象也较为符合，说明所定义“易损区”的合理性。

4 结 论

(1)通过对骨料形态数目进行统计，选取椭球状和盘状2种形态骨料为基础利用clump命令进行仿真建模，2种骨料模型与实际骨料的形态特征基本一致，为后续数值模拟试验的合理性奠定了基础。同时，为进一步提高骨料模拟精度，在后续研究工作中可以考虑引入三维激光扫描技术提取骨料模型，并在clump建模过程中增加组成clump的颗粒数来提高骨料建模精度，使骨料生成的粘结接触特点更符合实际；另外，可选取更多种类形态的骨料作为建模基础，提高骨料模型的代表性，从而保障骨料模拟的精度。

(2)随着孔隙率的上升，无砂混凝土模型内部颗粒数及颗粒粘结数下降，抗压强度逐渐降低，从孔隙率0.26时的24.9 MPa降低到孔隙率0.38时的7.97 MPa。同时抗压强度与孔隙率关系变化趋势与试验研究结论较为一致。

(3)数值模拟得出的应力-应变全曲线与基于试验得出的多孔混凝土应力-应变全曲线模型较为一致，说明PFC能较好地模拟无砂混凝土应力-应变特征。相比于常规混凝土，无砂混凝土在曲线上升阶段变化趋势更接近线性，而在下降段卸载更为迅速，下降趋势也出现了更强的离散型，体现了无砂混凝土高脆性的特点。

(4)通过观察破坏损伤模拟过程发现：无砂混凝土断裂多出现在边角及加载面附近，随着荷载的增加断裂数不断上升，在达到峰值强度后断裂数急剧上升贯通破裂面，数值模拟断裂情况与实际试验试块破坏形态基本吻合。

(5)基于峰值强度下断裂位置分布情况定义了无砂混凝土“易损区”，即除去试块内切球后的剩余部分，该区域体积为试块总体积的47.6%，而包含断裂数的比例却在69%以上。同时，该区域与实际试块受压破坏区域位置较为吻合，说明该定义区域具有一定的合理性。