基于Python的再生混凝土细观性能模拟方法研究

姚泽良，崔婷婷，段东旭，党发宁，闻 硕

(1.西安理工大学 土木建筑工程学院，陕西 西安 710048；2.西北干旱地区生态水利国家重点实验室，陕西 西安 710048)

再生混凝土是指利用再生粗骨料部分或全部代替天然骨料配制而成的混凝土，再生混凝土技术的开发和利用是建筑废弃物资源化利用的关键[1].由于再生骨料细观结构的复杂性和试验条件的差异性，单纯依赖试验的研究方法，使再生混凝土的研究止步不前.合理选取计算模型既可以减少由试验条件引起的系统性误差和由人为因素引起的偶然性误差，又可以将细微观结构和宏观结构性能联系起来.

目前，一些国内外学者在再生混凝土内部应力分布特征、随机骨料模型、不同级配下再生混凝土抗压强度等方面取得了一定的成果.Ma[2]等提出通过运用某种数值分析方法实现骨料的高堆积密度.肖建庄和李文贵[3]等通过对模型再生混凝土的数值模拟，研究了再生混凝土内部的应力分布特征和各项介质对再生混凝土试件力学性能的影响.党娜娜和彭一江[4]等基于蒙特卡罗法建立二维再生混凝土随机骨料模型，研究了再生混凝土的细观损伤机理.汪奔[5]提出了一种基于网格预先生成的二维随机凹凸型骨料的混凝土细观建模方法.然而，现有的再生混凝土细观力学模型还存在诸多不足.固定颗粒数目的计算模型不能体现出骨料的随机性[6]，再生骨料细观结构复杂，在建立细观模型时未能考虑其多界面的特性；大多数的再生混凝土骨料模型是建立在Matlab、Fortran、C++等计算机语言基础的开发，而这些语言并不具备定义材料属性、划分网格等功能，所以还需要建立好骨料模型后再导入有限元软件进行后续步骤，工作量很大，同时实际工程应用操作难度大.因此，再生混凝土细观力学模型有待深入研究.

为了解决构建符合真实骨料细观结构的再生混凝土骨料计算模型的关键技术难题，实现三维再生混凝土试件的仿真分析，本文将再生混凝土看作由原生骨料、老砂浆基体、老界面区、新砂浆基体、新界面区组成的三维细观结构，基于蒙特卡罗法运用具有动态语义特征的Python语言编写生成三维再生混凝土随机骨料模型的算法，通过Python语言对Abaqus件内核进行二次开发，实现模型构建和数值分析合理整合，采用该算法建立不同骨料含量的三维再生混凝土模型，主要分析投放率为45%的计算模型的力学性能，验证了该算法的有效性.

1 再生混凝土三维细观结构构建

1.1 三维细观设计

再生粗骨料来源于原生混凝土，本身性能就有很多差异,粗骨料的含量形状各不相同[7-8]，其细观结构较普通混凝土具有更高的复杂性和随机性.由于再生骨料附着砂浆(即老硬化砂浆)，再生混凝土内部必存在较多过渡界面区.骨料的分配[9-10]，用量和类型[11-12]在细观层次上对混凝土的抗压强度也有一定的影响，这里将再生混凝土看作由新砂浆基体、外界面区、老砂浆基体、内界面区、原生骨料五相材料组成[13].结合相关研究成果[14]，将再生骨料简化为球形.老砂浆基体、内外界面区由骨料进行面平移得来的球体，根据文献[15]，将老砂浆基体、内外界面区依据骨料半径占总骨料直径的比例关系设置合理的厚度值.经试算，选取1 mm厚度作为老砂浆基体的厚度，选取0.5 mm作为界面厚度，图1为再生混凝土剖面的三维细观结构.

图1 再生混凝土细观结构模型

1.2 骨料级配

骨料级配是指骨料中大小不同的颗粒之间的数量比例，骨料级配的选取直接影响到再生混凝土硬化后的各种性能指标，如强度、耐久性等.由于富勒级配曲线可产生最优的结构密度和强度，因此选取此曲线作为最大密实度级配曲线.富勒抛物线方程如下：

p=100(D/Dmax)0.5

(1)

(2)

式(1)中，p为骨料粒径在D以下的质量百分比；Dmax为最大骨料粒径.式(2)中，F(D)为粒径为D的原生骨料出现的概率，D0、Dmax为最小、最大粒径.根据上述理论，计算出随机球型骨料的数量，骨料的取代率为100%，骨料粒径为5～20 mm，定义的模型尺寸x×y×z=100 mm×100 mm×100 mm.

2 再生混凝土三维模型程序研发

2.1 基本原理

利用蒙特卡罗法产生随机数的方法也称为随机抽样法.本文选用乘同余数法产生伪随机数，乘同余数法的一般形式是对任意初值x1由如下递推公式确定：

(3)

其中:xn表示数序列；a表示乘子；M表示模；modM表示除以模后取余数.

通过构造随机模型使得某一随机变量的数学期望等于所求的解，最基本的随机变量在区间[0,1]上服从均匀分布.若产生[0,1]的均匀随机变量x，为了获取[a,b]上的随机变量X，则由式(3)经变换可得到：

X=a+(b-a)x

(4)

其中：X,x均表示随机变量；a,b表示常数.

基于上述的均匀分布函数，随机生成骨料的中心点坐标，以及用均匀分布的随机变量产生一个随机粒径Di.蒙特卡罗法以其特有的优点被广泛应用于计算数学和物理工程领域，成为解决许多复杂问题的重要手段.蒙特卡罗法的优点有：程序结构简单,对计算机性能要求不高，在计算机上容易实现；收敛速度与问题维数无关，受问题的条件限制影响较小.

2.2 程序研发

Abaqus软件具有很强的通用性，拥有前处理器、隐式求解器、显示求解器、后处理器等模块，可以解决许多复杂的非线性问题.Abaqus软件在混凝土和再生混凝土数值模拟方面应用较广，但由于软件本身的局限性，使得自动化建模、参数化分析等研究无法进行，需要通过二次开发或结合其它软件才能实现.结合Python语言，程序设计流程如图2所示.

图2 设计流程

在Python脚本的编制过程中，首先定义模型的尺寸、骨料所占的百分比、骨料界面区的厚度以及老砂浆的厚度，其次根据基于以上骨料级配理论内置一个粒径分布函数，基于蒙特卡罗法内置一个random随机函数，从粒径d1,d2中生成一个随机浮点数，并获取粒径的随机数(d=d2+random()*(d1-d2))，根据设置的投放量的要求将骨料粒径分级并存入到数组中.

投放空间划分为均匀的小格子，调取随机粒径函数依次投放骨料，将骨料按级配随机投放到划分的均匀网格中，记录投放次数以及骨料的投放点的中心位置，设置不同骨料间的净距条件，若新投放的骨料与已投骨料之间无交点，则记录投放的中心位置，存入数组中；若有交点，则删除新投放骨料，重新循环投放直到得到满足要求的骨料模型，再通过面平移得到新旧砂浆的界面区.最后运用Abaqus和Python之间的交互作用[16]，在Abaqus中进行装配以及网格和荷载的布置等，开始计算分析.

2.3 生成模型验证

将骨料粒径分为三级，分别为20～15 mm、15～10 mm、10～5 mm，再生骨料最大粒径为20 mm，最小粒径5 mm.骨料投放的区域为边长为100 mm的立方体，其中投放率设置在45%～60%之间，图3为按投放率阶梯生成的三维细观模型.

图3 不同骨料含量模型

经与文献[3-4]对比分析，本文方法建立的计算模型具有明显优势：本文方法可快速生成不同骨料含量的计算模型，建模效率高；模型骨料符合级配要求，与再生混凝土真实细观结构吻合较好；不同大小的骨料随机排列，分布均匀，满足随机性要求；骨料投放率可达60%，满足数值计算要求.

3 再生混凝土细观性能模拟计算

为了进一步验证再生混凝土三维细观模型有效性，这里对骨料含量为45%的模型进行了单轴压缩和单轴拉伸模拟计算分析.

3.1 计算过程

选用Abaqus软件自带的损伤塑性本构模型(简称CDP模型)，根据CDP模型要求和《混凝土结构设计规范》中相关规定，定义再生混凝土的偏心率(0.1)、屈服面系数(0.666 67)、粘性参数(1×10-5)等相关塑性参数，各参数值存储于mat数组中.根据相关试验资料[14]，设置新砂浆、老砂浆、外界面、内界面、骨料的弹性模量E、抗压强度fc、抗拉强度ft、泊松比γ、密度ρ等参数.相关各相材料参数见表1.

表1 材料参数[14]

模型装配完成并划分网格后，进行创建分析步、修改输出、设置参考点、设置约束等步骤.本文采用刚性位移加载方式，即在再生混凝土计算模型的上表面设置耦合约束，在其顶部参考点分别施加-0.5 mm和1 mm的位移荷载(图4)，对三维再生混凝土模型开展了单轴压缩和单轴拉伸下的数值模拟计算分析.采用静力，通用分析步，总分析步时间为1 s.

图4 模型荷载施加

3.2 结果分析

3.2.1 抗压、抗拉强度和峰值应变

对试件施加位移荷载，加载结束后，为了减少误差，本文分别提取三组随机模型的应力-应变值，图5为计算分析得到的应力-应变曲线.

由图5(a)可知，在单轴受压下，模拟试件的抗压强度为22.56 MPa，对应的峰值应变为0.002 4，抗压强度介于新砂浆强度和老砂浆强度之间，远低于原生骨料的抗压强度.由图5(b)可知，模拟试件的抗拉强度为2.18 MPa，对应的峰值应变为1.13E-5，抗拉强度约为抗压强度(22.56 MPa)的十分之一.

图5 应力-应变曲线

经与文献[14]的试验数据(抗压强度为21.79 MPa，峰值应变为0.001 92；抗拉强度为2.06 Mpa，峰值应变为1.08E-5)对比，计算得到的强度和峰值应变略高于试验值.分析其主要原因为：由于初始真实再生混凝土浇筑完成后是有微裂缝的，而本文的计算模型假设无初始裂缝.所以在力学性能的分析验证上，其抗压、抗拉强度相对于试验均有所提高，计算结果是合理的.

3.2.2 损伤破坏和应力分布

提取x向、y向、z向的应力，分析再生混凝土应力分布特征.图6、图7分别为单轴受压和单轴受拉时计算模型表面的x向水平、y向竖向、z向水平以及损伤云图(图中“+”为拉应力，“-”为压应力).

图6 单轴受压时应力和损伤状态

图7 单轴受拉时应力和损伤状态

分析表明，x向水平应力S11(图6(a))的最大应力值出现在新硬化砂浆区和天然骨料区，其应力值为-0.22～2.43 MPa；y向竖向应力S22(图6(b))的最大应力值出现在老硬化砂浆区以及内、外界面区域，其应力值为-0.245～2.503 MPa；z向水平应力S33(图6(c))的最大应力值分布在内、外界面区中，其应力值在为-0.14～2.266 MPa.由图7可知，S11水平向最大应力值为2.112 MPa；S22竖向最大应力值为22.319 MPa；S33Z向应力为2.101 MPa.将以上结果与文献[14]中的试验数据(单轴压缩荷载作用下的强度为2.44 Mpa)对比，二者吻合较好，说明计算结果合理.加载结束后，试件的各个表面均出现了损伤破坏带(图6(d)和图7(d))，受压最大损伤值为0.857 6，受拉最大损伤值为0.942 1.

以上分析表明，通过Python操纵Abaqus软件内核，将再生混凝土看作由原生骨料、老砂浆基体、老界面区、新砂浆基体、新界面区组成的三维细观结构，实现了模型构建和数值分析合理整合.基于蒙特卡罗法生成三维再生混凝土随机骨料模型的算法较好地解决了构建再生混凝土三维随机骨料模型的关键技术难题，为再生混凝土基本理论研究提供了新方法.同时，该算法只需在脚本中修改相关参数就可进行相应计算，满足工程需要.

4 结论

运用Python语言编写基于蒙特卡罗法的三维再生混凝土骨料模型算法，解决了Abaqus软件在再生混凝土计算方面的技术难题，实现了模型构建和数值分析合理整合，通过模拟计算分析，验证了该算法的准确性.所得的主要结论如下：

(1)通过Python建立的三维再生混凝土随机骨料模型与再生混凝土真实细观结构相吻合，满足生成骨料的随机性要求，模型高效准确.

(2)通过Python语言直接操纵Abaqus内核进行二次开发，利用二者之间的交互作用可快速生成不同骨料含量的模型，建模效率高，骨料投放率可达60%，满足数值计算要求.

(3)三维再生骨料模型的力学性能指标与试验结果吻合较好，表明该方法计算准确.计算所得强度和峰值应变略高于试验值，其主要原因是三维再生骨料模型均在细观结构上进行了简化，忽略了再生混凝土初始细微裂缝，以及浇筑成型后试块表面的微突起.