基于粘结单元的三维随机细观混凝土离散断裂模拟

杨贞军，黄宇劼，尧 锋，刘国华

(1. 武汉大学土木建筑工程学院，湖北省岩土与结构安全重点实验室，武汉，湖北 433000；2. 浙江大学建筑工程学院，杭州，浙江 310058)

传统的宏观混凝土模型假设均质材料，不能反映细观组分空间分布的随机性和非均质性引起的损伤与断裂局部化，难以可靠预测材料的力学响应与裂缝扩展过程。因此进行混凝土细观断裂模拟研究十分必要[1 − 5]。在三维细观混凝土模拟中，假设骨料为球体是一种常用的简化建模方法[4,6 −7]，虽然可以高效地生成细观结构，但无法反映骨料形态的影响。近年来，采用更为复杂的椭球和凸多面体来模拟骨料日趋常见[8 − 9]。在应用上，复杂骨料模型大多用于几何特征研究[8,10]、边界效应分析[11]、渗透性能研究[12 − 14]以及线弹性应力分析[15 − 17]；在三维复杂非线性损伤和断裂的模拟方面仍非常有限，且大多采用连续介质损伤塑性模型[18 − 21]，难以定量预测裂缝宽度，因此有必要采用离散裂缝模型进行更精细的断裂模拟。另一方面，在骨料几何建模方法上，一般通过自编程序进行参数化建模[18,21]，或基于Voronoi 图生成细观结构[22]，或基于ANSYS 的ADPL 命令流建模[23]。但是，对三维复杂细观结构进行网格划分并不容易，较为直接的方法是将细观结构映射到规则的立方体背景网格上[4, 18 − 20, 24]，但要保证精度则会导致巨大的单元数量。ABAQUS 具有强大的前处理模块，可以灵活有效地对复杂结构进行网格自动划分，不过尚未见到将其用于三维混凝土细观结构建模与网格划分的研究。

本文首先采用Python 编写程序，利用ABAQUS脚本接口[25]进行前处理二次开发，考虑骨料实际级配，有效地建立多面体随机骨料几何模型并进行有限元网格划分。然后通过自编C++程序在骨料-砂浆界面上以及砂浆中高效插设零厚度的离散粘结裂缝单元[26 − 28]，模拟受拉试件复杂细观三维多裂缝的起裂与扩展直至破坏的全过程，并对粘性裂缝单元的主要材料参数(抗拉强度与断裂能)对荷载-位移曲线、断裂过程、裂缝面特征的影响进行详细分析，以加深对混凝土细观断裂机理的认识。应该指出，笔者论文[26 − 28]也采用了插设界面单元的方法，但均非基于本文提出的随机骨料模型；虽然文献[26 − 27]采用随机场来间接表征混凝土的空间非均质性，但和文献[28]一样，仅提出了均质材料中插设二维、三维界面的方法。本文是前述工作的扩展，并将Matlab 程序改进为C++，大大提高了效率，更适合大规模三维有限元网格的处理。

1 基本算法

1.1 随机多面体骨料建模方法

算法分两步，第一步建立几何模型，第二步建立有限元模型。

首先编写MATLAB 程序进行多面体随机骨料的生成与投放，建立几何模型。采用基于球体的方式产生凸多面体，即凸多面体的各顶点位于球面上。该球体称为基球，其球面上随机产生8 个～25 个点作为多面体的顶点，顶点i 的笛卡尔坐标用Euler 空间的极角αi和方位角βi表示为：

式中：(x0, y0, z0)是球心的坐标；r0是球半径；αi在[0, 2π]、βi在[0, π]范围内按均匀分布随机产生。然后采用MATLAB 算法“convhulln”基于顶点建立三角形面，这些面构成骨料多面体。为避免顶点距离太近而造成三角形面畸形，需保证球面上的顶点间距不小于ξr0，本文取ξ=0.5。多面体的生成较为耗时；为提高效率，本文不采用边投放、边生成骨料的方式，而是通过生成足够数量(本文为20000)的直径为单位一的多面体，预先建立骨料形态数据库，图1 为一些骨料例子。

图1 不同形态的多面体骨料Fig.1 Polyhedral aggregates of different shapes

针对给定模拟区域，具体算法如下：

1) 根据所需骨料含量，采用实际骨料级配表1求解各级配的骨料体积，从骨料直径最大的级配开始分级配生成并投放骨料；

表1 三级配骨料尺寸分布[29]Table1 Three-graded aggregate size distribution[29]

2) 对每个级配的操作步骤：随机产生一个球心坐标作为投放位置，从骨料数据库中随机选取一个骨料，并在此级配区间内按均匀分布随机赋予其直径；

3) 若该骨料不与已有骨料相交或重叠，则该骨料投放成功，进行下一个投放，否则重新进行本次投放；

4) 本级配骨料体积达到后，进行下一级配投放；

5) 循环2)步～4)步，直至达到总骨料含量。

在上述算法第3)步中，新生成骨料与已有骨料之间可能存在如图2 所示的3 种位置关系：不相交也不重叠、相交、重叠(包含)。

图2 两个骨料之间的三种位置关系Fig.2 Spatial relations of two aggregates

通过下述算法判断两个多面体相交或重叠：如图2(a)所示，对于新生成的多面体A 和已有多面体B，要使二者不相交且不重叠，则B 的各面需在A 各面的同一侧。A 的顶点A1～A3 所构成平面A123 的方程是：

式中：(xAi, yAi, zAi)是顶点Ai(i=1～3)的坐标，则对于平面A123 外的点(x, y, z)，有F(x, y, z) > 0 或F(x, y, z)<0。则上述位置判断问题可以等效为：遍历A 的各面，如果存在一个面(例如A123)，使得A 的内部各点(凸多面体可只选取其形心)与B 的任意顶点(xBj, yBj, zBj)位于面A123 的不同侧，则可判定A 与B 既不相交也不重叠。相比于遍历A、B 的顶点以保证各自顶点不在另一个多面体内部的算法，上述位置判断算法更为直观和有效。若A 的形心坐标为(xA0, yA0, zA0)，A 与B 既不相交也不重叠时只需满足下式：

骨料与骨料、骨料与边界之间预设最小间距(本文设置为较大骨料直径的0.05 倍)以方便网格生成。

第二步，编写Python 脚本程序，进行ABAQUS 前处理二次开发，建立多面体骨料并完成模型网格划分。具体算法如下：

1) 建立模型数据库；

2) 生成混凝土试件立方体部件(part)；

3) 读取随机骨料几何模型中所有骨料的信息：基球的球心坐标、顶点坐标、各面的顶点构成；

4) 对于每个骨料，生成三维球体部件；

5) 延拓骨料的各面，切割该球体，过程见图3；

6) 在assembly 模块中创建部件实例(instance)，采用布尔操作合并生成的多面体与试件立方体；

7) 重复上述步骤，直至创建所有骨料，见图4；

图3 在ABAQUS 中通过Python 脚本程序生成多面体Fig.3 Polyhedra generation using Python scripts in ABAQUS

图4 生成的50 mm 试件多面体骨料有限元模型Fig.4 FEM model of a specimen with polyhedral aggregates

8) 布置全局网格种子，采用“自由网格”技术用四面体单元进行网格划分。

以边长为50 mm 立方体试件和骨料含量30%为例，上述算法建立的骨料有限元模型见图4。该模型含有423 个骨料，网格单元平均尺寸为2 mm。

1.2 插设三维粘结裂缝单元

本文将ABAQUS 零厚度三维六节点粘结裂缝单元COH3D6 插设到砂浆的有限元实体单元之间，以及骨料-砂浆界面上，用于模拟混凝土复杂非线性断裂过程。采用损伤系数D [0, 1]描述裂缝单元达到强度后的损伤程度，反映刚度退化，该系数是裂缝单元有效相对位移δm的函数：

式中：δn、δt、δs分别是法向、切向的相对位移；< >为Macaulay 括号，表示为：

以线性软化准则为例，损伤系数可以写成：

式中：δm,max是加载历史中的最大有效相对位移；δm0和δmf分别是裂缝起裂和完全破坏时的有效相对位移。

法向刚度kn与切向刚度ks、kt分别表示为：

相应的应力则可以表示为：

采用名义应力平方准则作为起裂准则:

式中，tn、ts和tt分别为法向、切向的应力。

文献[26 − 27]报道了在二维、三维实体有限元单元网格中插设粘结单元的算法和MATLAB 程序实施，用于均质材料。本文将该算法加以拓展，通过编写C++程序，在混凝土细观结构中插设三维粘结单元，图5 给出了算法示意图。假设骨料不会开裂，骨料内无需插设粘结单元，因此只在骨料-砂浆界面以及砂浆中采用粘结单元(图5(c))，分别标记为CIE_INT 与CIE_CEM。

图5 细观混凝土三维粘结单元的插设算法Fig.5 3D inserting algorithm of cohesive elements

具体插设算法如下：

1) 读取初始网格，图5(a)显示了网格中与骨料表面相连的2 个砂浆四面体单元A 与B，以及1 个与砂浆单元B 相连的砂浆单元C。

2) 如图5(b)所示，搜寻单元的面，骨料的面存入结构体FACE_INT 中，砂浆的面存入FACE_CEM 中。

3) 如图5(b)所示，将节点按其所在单元的数量拆开(即不改变坐标地进行节点的原位复制)并重新编号，再分配给各单元。

4) 如图5(c)所示，根据FACE_INT、FACE_CEM 分别建立粘结单元CIE_INT、CIE_CEM，从而完成粘结单元的插设，图6 给出一个骨料表面插设的零厚度粘结单元CIE_INT 作为示例。

5) 输出插设粘结单元后的整体网格信息，作为ABAQUS 可读取的*.inp 输入文件供计算使用。

对于图4 中约有3.2 万个节点、18.1 万个四面体单元(C3D4)以及36.7 万个单元面的原始模型，采用原有MATLAB 程序[27 − 28]进行插设耗时5 h，而本文C++程序不到3 min。插设粘结单元之后，模型大约有53.1 万个节点数、27.0 万个粘结单元，其中骨料-砂浆界面上大约有2.6 万个粘结单元。

图6 骨料表面的粘结单元(CIE_INT)Fig.6 Cohesive element on the surface of an aggregate

2 单轴受拉模拟

骨料和砂浆假设为线弹性材料，粘结单元采用线性受拉/受剪软化关系和名义应力平方起裂准则进行模拟。采用文献[22]中的材料参数(表2)，并假设受剪与受拉断裂参数相同(偏于保守)。粘结单元的初始弹性刚度应足够高以模拟未开裂初始状态，但刚度过高会导致算法不收敛。经过多次尝试，将其取为6×106MPa/mm。

表2 材料参数[22]Table2 Material parameters[22]

对生成的50 mm 的细观混凝土立方体试件(图4)进行单轴受拉模拟。边界条件如图7 所示，左端面全约束，右端面受均匀分布水平(x 向)位移荷载控制，最终位移0.06 mm。使用ABAQUS/Explicit 显式求解器进行准静态计算分析，设定荷载步时间为0.005 s。

图7 单轴受拉混凝土试件的边界条件Fig.7 Boundary conditions under uniaxial tension

图8 给出了该模型采用两组不同断裂能时的平均应力-位移曲线，其中一组采用表2 的参数(即参照组)，另一组Gc0.06Gi0.01 表示界面粘结单元的断裂能为0.01 N/m，其余参数同表2。图中也给出了典型单轴受拉试验数据[30]作为对比，该试验试件的尺寸是50 mm×50 mm×60 mm，与模型比较接近。可见Gc0.06Gi0.01 得到的模拟结果与试验吻合较好，表明模型能够有效预测混凝土力学响应。但由于试验试件与本文数值模型在细观组分方面并不相同，该试验数据仅作参考。

图8 单轴受拉混凝土试件的宏观应力-位移曲线Fig.8 Average stress - displacement curves of concrete specimen under uniaxial tension

2.1 开裂过程分析

图9 显示了采用参考材料参数(表2)模拟的试件正面和背面视角的开裂过程(对应于图8 曲线的A 点～F 点)。其中红色的单元为损伤因子D≥0.9的粘结单元，在本文中假设为已完全损坏的离散裂缝。为了观察裂缝表面形态，图10 显示了去掉这些损坏裂缝单元后的相应变形过程。为了方便观察，变形放大系数均设为200。

图9 单轴受拉开裂过程：粘结裂缝单元的受拉变形Fig.9 Fracture process under uniaxial tension: the red cohesive elements with D≥0.9

如图9 和图10 所示，在早期加载阶段(A 点之前)，由于界面粘结单元的断裂参数比砂浆低，骨料-砂浆界面上出现大量微裂缝，但其损伤因子D 仍较低，反映在裂缝单元的张开位移较小，以及荷载-位移曲线在加载位移较小时即显示出非线性特征；随着加载位移的增加直至试件达到抗拉强度(A 点)，微裂缝宽度缓慢发展，但裂缝并不显著；当试件进入软化段时(B 点～D 点)，一些微裂缝的宽度迅速增加进而局部化形成主裂缝，并与砂浆中新生成的裂缝连通，而另外一些微裂缝则逐渐卸载闭合；如图中阶段(E 点～F 点)所示，随着位移进一步增加，连通的裂缝迅速张开变宽，形成一条贯穿主裂缝，试件拉裂。这和混凝土单向受拉实验[30]和其他数值模拟[22,28]的断裂扩展过程一致。

2.2 材料断裂参数分析

现有研究表明[26 − 27]，粘结裂缝单元的材料参数，特别是抗拉强度与断裂能，对模拟结果影响较大。本文对砂浆中与骨料-砂浆界面上的粘结单元选取不同组别的抗拉强度与断裂能进行研究。保持断裂能不变(见表2)，粘结单元取不同抗拉强度时，使用TcTi 的编号规则，如Tc4Ti2 表示砂浆和界面粘结单元的抗拉强度分别为4 MPa 和2 MPa；保持表2 中抗拉强度不变，取不同断裂能，使用GcGi 的编号规则，如Gc0.06Gi0.03 表示砂浆和界面粘结单元的断裂能分别为0.06 N/mm 和0.03 N/mm。将采用表2 材料参数的模型设置为参考组(编号为Tc4Ti4_REF 与Gc0.06Gi0.03_REF)。

图10 单轴受拉开裂过程Fig.10 Crack propagation process under uniaxial tension

取8 组粘结单元材料参数进行模拟，对每一组，只变动TcTiGcGi 参数中的一个，其余参数与参考组相同。得到如图11 和图12 所示试件宏观平均应力-位移曲线，以及如图13 和图14 所示的试件破坏形式。

图11 粘结裂缝单元的抗拉强度对宏观应力-位移曲线的影响Fig.11 Average stress - displacement curves: effects of fracture energy of cohesive elements

图12 粘结裂缝单元的断裂能对宏观应力-位移曲线的影响Fig.12 Average stress - displacement curves: effects of fracture energy of cohesive elements

图13 粘结裂缝单元的抗拉强度对试件破坏形式的影响Fig.13 Effects of tensile strength of cohesive elements on the fracture surfaces

由图11 可见，粘结单元的抗拉强度对试件应力-位移曲线峰值以及峰值前非线性段影响显著：随着粘结单元强度的增加，试件强度增加；当骨料-砂浆界面裂缝单元强度降低而砂浆裂缝单元强度增加时，试件强度亦增加，表明砂浆粘结单元的抗拉强度对试件承载力起控制作用。图11 也表明：随着试件强度的增加，其脆性也逐渐增加。由图12 可见，粘结单元的断裂能对试件应力-位移曲线软化段影响显著，随着断裂能的提高，试件的延性增加，但对试件强度以及曲线峰前段影响不大。

图14 粘结裂缝单元的断裂能对试件破坏形式的影响Fig.14 Effects of fracture energy on crack surfaces

由图13 可见，粘结单元抗拉强度对裂缝面形态有较大影响，该影响主要体现于砂浆、界面粘结单元的强度相对比值γT=Tc/Ti。当γT≠1 时，断裂面的位置基本相同，γT主要影响断裂面的形态；当γT=4 时(图13(c)设计组Tc4Ti1)，断裂面上出现的大骨料数量最多，其次是γT=3(图13(b)设计组Tc6Ti2)，然后是γT=2(图13 参考组Tc4Ti2)。由于本文粘结单元以名义拉应力作为起裂准则，随着砂浆、界面粘结单元强度差别的增大，界面愈发成为薄弱环节。由于大骨料表面的界面单元面积较大，较容易成为裂缝起裂与扩展通道，使试件倾向于在大骨料附近形成裂缝面。另外，当γT=1 时，如图13(a)和图13(d)所示，试件断裂面较平滑且以小骨料数量居多；其中设计组Tc2Ti2的粘结单元强度较低，离加载端最近一层裂缝单元因应力集中而破坏。

图14 显示砂浆、界面粘结单元的断裂能相对比值γG=Gc/Gi对裂缝形态的影响。当γG≤3 时(图14(a)、图14(b)、图14(d)和参考组)，试件破坏时的位移和断裂面形态基本一致；随着断裂能的提高，裂缝面略趋于曲折；但当γG=6 时(图14(c)设计组Gc0.06Gi0.01)，砂浆裂缝单元的断裂能显著高于界面单元，界面成为薄弱环节，由于大骨料表面的界面单元面积较大，较容易成为裂缝扩展通道因而影响了裂缝扩展路径，使最终断裂面的位置与形态与其他组差异较大。

进一步将Gc0.06Gi0.01 与参考组Gc0.06Gi0.03_REF 比较，分析断裂能、宏观应力-位移曲线以及开裂过程之间的关联。如图15、图16 所示，在这两个模型的应力-位移曲线的软化段中选取4 个应力水平接近的阶段(A 点～D 点)进行分析。由图可见，由于Gc0.06Gi0.01 的裂缝沿着图14(c)黑色箭头所指的大骨料(大骨料的界面作为薄弱环节)扩展，导致其裂缝面的形成比Gc0.06Gi0.03_REF 更为曲折复杂，因此其应力-位移曲线的软化段由于断裂耗能的需要变得较为平缓，而没有呈现出Gi减小导致的脆性特征。由此可见，混凝土的力学响应能够反映其裂缝发展特征，二者既决定于断裂材料参数，也受到骨料大小、形状等细观结构因素的影响，显示出十分复杂的破坏机理。

图15 断裂能对开裂过程的影响Fig.15 Effects of fracture energy on crack propagation processes

图16 断裂能对应力-位移曲线的影响Fig.16 Effects of fracture energy on stressdisplacement curves

3 结论

本文采用Python 编程，利用ABAQUS 脚本接口进行前处理二次开发，建立混凝土三维随机多面体骨料的细观有限元模型；通过高效插设离散粘结裂缝单元，成功模拟了单向受拉复杂三维多裂缝的起裂与扩展。主要结论如下：

(1) 在早期加载阶段，由于界面粘结裂缝单元的断裂参数比砂浆低，骨料-砂浆界面上出现大量微裂缝，荷载-位移曲线较早显示出非线性特征；随着位移的增加直至试件达到抗拉强度，微裂缝宽度缓慢增加；进入软化段后，一些微裂缝宽度迅速增加，并与砂浆中新生成裂缝连通形成局部化的主裂缝，而另外一些微裂缝则逐渐卸载闭合。

(2) 宏观应力-位移曲线主要受砂浆、界面粘性裂缝单元的抗拉强度和断裂能绝对数值的影响：砂浆粘结裂缝单元的抗拉强度对试件承载力起控制作用，试件承载力随着其强度的提高而增大；而断裂能对非线性软化段影响显著，断裂能的提高增加了试件的延性。

(3) 裂缝面的位置和形态主要受砂浆、界面粘性裂缝单元的抗拉强度、断裂能相对比值的影响。比值较大时(γT>2 或者γG>3)，即界面粘结单元的力学性能远弱于砂浆粘结单元时，由于大骨料表面的界面单元面积较大，较容易成为裂缝起裂与扩展通道，使试件倾向于在大骨料附近形成裂缝面。

(4) 混凝土的力学响应反映其裂缝发展特征，二者既决定于断裂参数，也受到细观结构的影响，本文建立的模型能够有效地描述混凝土复杂三维断裂过程。