混凝土三维细观模型的建模方法与力学特性分析*

张煜航，陈青青，张 杰，王志勇，李志强，王志华

（1.太原理工大学机械与运载工程学院应用力学研究所，山西 太原 030024；2.太原理工大学材料强度与结构冲击山西省重点实验室，山西 太原 030024）

混凝土材料从细观层次上由骨料、砂浆以及两者之间的黏接界面层（interfacial transitional zone，ITZ）组成，各组分的细观力学性质直接或间接地决定了混凝土材料的宏观力学性能和破坏模式[1-2]。根据混凝土细观结构特点及其力学性能，深入研究细观结构和宏观力学性能之间的关系，是混凝土力学性质的重要研究内容之一。

Wittmann 等[3]首先将多尺度的研究方法应用到混凝土材料的研究中，认为混凝土是一种非均质复合材料，提出了用随机分布的多边形来模拟混凝土骨料的细观力学模型。该方法一直被后来的研究者所采用。Zhou 等[4-5]建立了圆形随机骨料模型，利用材料损伤本构关系，系统研究了混凝土细观模型在不同受力条件下的拉伸、压缩和爆炸问题，指出骨料和砂浆之间的黏接界面层在压缩时先出现破坏裂纹，在拉伸时对混凝土失效机理和拉伸强度影响最大。Wang 等[6]在椭圆形和多边形骨料的基础上引入孔隙，并研究了骨料形状、骨料体积分数和孔隙率对抗拉强度的影响。吴成等[7]对刚性弹丸侵彻细观混凝土靶进行了数值模拟，分析了砂浆种类、粗骨料种类和粗骨料体积分数等对靶板抗侵彻能力的影响，并通过扩展Forrestal 阻力方程，建立了细观混凝土侵彻深度模型。相比于二维模型，三维细观模型更接近于实际混凝土的真实形态，因此能较好地反映混凝土材料的实际变形与损伤破坏情况。Wang 等[8]提出了“投放算法”用于产生随机分布的三维骨料模型，并在单元之间嵌入零厚度黏接单元，以此来模拟混凝土在拉应力作用下的开裂行为。Zhang 等[9]通过对生成的骨料进行平移和旋转来控制骨料的移动，提出了“随机爬行算法”，该算法可以实现更高骨料含量的混凝土建模。邓勇军等[10]基于骨料随机投放的思想建立了混凝土三维细观几何模型，分析了刚性弹正侵彻过程中发生弹道偏转的原因及可能的影响因素，定量研究了混凝土细观因素对弹道偏转的影响。然而，骨料投放往往需要复杂的算法对骨料相互之间是否发生侵入进行判断，同时需要花费较多的时间，建模效率较低。三维Voronoi 图形中胞元形状与混凝土骨料具有良好的相似性，其胞元可实现随机的空间形状和空间分布，可作为建立骨料的一种方法。

本文中基于三维Voronoi 图形，提出一种简单、高效的混凝土细观模型生成方法。在传统Voronoi 图形的基础上，通过控制多面体的随机度和引入缩放因子等参数，得到具有级配的随机骨料模型。然后，对骨料几何体的外表面进行延伸，得到一定厚度的黏接界面层，从而得到完整的混凝土细观有限元模型。采用连续介质损伤模型分析混凝土材料的静态和动态力学特性，拟为该模型的进一步应用提供基础。

1 混凝土三维细观模型

1.1 细观几何模型

三维Voronoi 图形是一组由连接两顶点直线的垂直平分面形成的连续多面体。在一个特定的空间内，随机分布N 个种子点，Si代表第i 个胞元的核心坐标，相应的胞体由顶点P 围成，顶点P 到所属胞元的核心Si的距离小于或等于到其他胞元核心Sj的距离：

为了避免产生过于畸形的胞元，通过下式控制两核心之间的最小距离：

式中： q( 0＜q＜1)为缩放因子，通过改变q，可以实现对多面体缩放大小的控制，从而得到满足一定级配的随机骨料模型（见图3）。

在生成的骨料模型基础上，利用下式对多面体顶点沿着向量 v进 行延伸，并保留多面体顶点，得到顶点（见图4）：

式中： η( 0＜η＜1)为控制黏接界面层厚度的延伸因子。图5 为利用此方法生成的具有不同界面层厚度的骨料示意图。

图 1 K=0.2 时三维Voronoi 多面体Fig.1 3D Voronoi polyhedron when K=0.2

图 2 骨料缩放示意图Fig.2 Schematic diagram for single aggregate

图 3 不同大小骨料模型Fig.3 Aggregate cells with grading sizes

图 4 ITZ 产生示意图Fig.4 Schematic diagram of generating ITZ

图 5 具有不同ITZ 厚度的骨料外轮廓图Fig.5 Geometry outlines of different ITZ layer thicknesses

1.2 细观有限元模型

上述生成细观模型的方法具有过程简单、随机性可控制以及生成速度快[11]等特点，但多面体缩放会降低初始生成骨料的体积含量，达不到实际混凝土对体积分数的要求。通过重力下落过程[12]不仅可以提高模型中骨料的体积分数，同时可以改善骨料分布，增强随机性，使模型更接近真实混凝土。一定区域内不同骨料体积分数的模型如图6 所示。由于实际混凝土材料中黏接界面层厚度仅为10～50 μm，远小于骨料尺寸，采用四面体划分网格数量巨大，严重耗费计算时间。因此，对于ITZ 的网格划分采用楔形单元，在保留模型计算精度的同时，可以最大程度地降低网格数量，提高计算效率。图7 为试样尺寸为25 mm×25 mm×25 mm、骨料体积分数为37.38%的有限元模型，其中绿色区域为骨料（四面体单元72 114 个），红色区域为黏接界面层（楔形单元25 408 个），界面层厚度0.1 mm，蓝色区域为砂浆（四面体单元36 727 个）。

图 6 含有不同骨料体积分数的立方体试样Fig.6 Cubic specimens with different volume fractions of aggregate

ABAQUS 有限元软件中塑性损伤模型主要用于模拟砂浆、岩石和混凝土等材料在低静水压力下由损伤引起的不可恢复的材料性能，该模型考虑了拉压性能的差异，可以较好地描述砂浆和黏接界面层的力学特性[13-16]。将图7 模型导入ABAQUS 有限元软件，采用塑性损伤模型模拟混凝土细观组分中砂浆和黏接界面层的力学性能，砂浆塑形损伤模型曲线如图8 所示，界面层相应曲线由砂浆曲线按强度进行等比例折减得到。通常骨料强度远高于砂浆强度，因此骨料采用线弹性模型。具体材料参数见表1。

图 7 三维混凝土细观模型有限元模型Fig.7 Meshing results for concrete specimen

图 8 砂浆单轴准静态压缩应力应变曲线Fig.8 Strain-stress curves of mortar phase under uniaxial load

表 1 细观组分材料参数Table 1 Material parameters of three-phase materials

2 三维细观模型静力学分析

2.1 单轴静力特性分析

图 10 单轴压缩与拉伸应力应变曲线Fig.10 Comparison of strain-stress curves between experiment and simulation

图9 为试样单轴加载示意图，试样上下两端为刚性平面，上刚性板采用位移加载，下刚性板完全固定。对于准静态问题，隐式算法不受试样惯性效应以及加载动能的影响，更能反映试样真实受力状态，本文中采用ABAQUS/standard 隐式计算模块模拟准静态加载下的压缩和拉伸。单轴压缩和拉伸应力应变曲线以及相应的实验对照如图10 所示。由图10可知，数值模拟与实验结果[17]曲线趋势相同，峰值误差在10%以内，吻合较好。对于单轴压缩，刚性板与试样之间的摩擦条件的不同，对混凝土试样的破坏模式有较大影响。不同摩擦条件下混凝土单轴压缩破坏模式的对比如图11 所示。在低摩擦条件下，混凝土裂缝主要为平行于加载方向的竖向裂缝；在高摩擦条件下，裂缝转变为八字形裂缝，与通常实验结果[17]相似。图12 为单轴拉伸载荷作用下试样的破坏模式。裂缝主要为垂直于加载方向的环向裂缝，符合实际实验结果。

图 11 不同摩擦条件下试样破坏模式Fig.11 Failure patterns of different frictional conditions

图 12 单轴拉伸试样破坏模式Fig.12 Failure patterns under uniaxial tension

2.2 双轴静力特性分析

尚世明[18]采用定测压加载方式，得到了标准立方体试样（100 mm×100 mm×100 mm 的混凝土试件）在不同侧压下的宏观应力应变曲线和相应的破坏模式。图13 为数值模拟试样双轴加载示意图，试样上下两端为刚性平面，下端固定，上端位移加载。侧向方向施加均匀的压应力，分别为0、8、16、24 MPa。图14为数值模拟与实验在位移加载方向的应力应变对比图。由于试样尺寸、骨料级配以及应力加载等原因的影响，导致数值模拟与实验结果存在误差，但两者趋势相同，峰值应变也较接近，可以认为误差在合理范围之内。不同侧压应力条件下试样的破坏模式如图15 所示。随着侧向压应力不断增大，对试样的约束作用不断增强，受压面裂缝逐渐沿垂直于自由面方向产生，破坏形态转变为层状破坏。由图15 可知，数值模拟与实验破坏模式基本吻合，说明该模型能有效反映混凝土试样在双轴受压载荷下的力学特性。

图 13 定测压加载示意图Fig.13 3D model under biaxial simulation

图 14 双轴压缩应力应变曲线Fig.14 Comparison of strain-stress curves between experiment and simulation

2.3 三轴静力特性分析

为了进一步研究复杂应力状态下混凝土材料的静力特性，基于本文的三维细观模型，对常规三轴条件下混凝土的应力应变关系进行了研究。图16 为不同围压下试样位移加载方向应力应变曲线。如图16 所示，随着围压不断增大，峰值应力增大明显，软化段逐渐模糊，当围压达到30 MPa 时，会出现应力随应变增大而增大的强化段，与文献[19]的实验曲线趋势相似。

图 15 不同侧向压应力下试样破坏模式Fig.15 Failure patterns of different frictional conditions

图 16 三轴压应力作用下应力应变曲线Fig.16 Comparison of stress-strain curves between experiment and simulation

3 三维细观模型动力学分析

混凝土作为一种典型的率相关材料，在动态载荷作用下表现出不同于静态载荷的力学特性[20-22]。运用ABAQUS 显式算法，模拟低应变率下混凝土SHPB 冲击破坏实验，其中入射杆和透射杆为弹性金属材料，采用弹性模型，弹性模量为210 GPa，泊松比为0.28。细观混凝土模型仍采用与准静态模拟相同的材料模型与材料参数。图17 为计算得到的4 种应变率下抗压强度动态增强因子。准静态模拟得到的抗压强度为24 MPa（无摩擦边界），即 fc=24 MPa。从图17 可以看出，抗压强度随着应变率的增大而增大，这与常规混凝土SHPB 实验所得结论一致，说明模型在动态载荷作用下同样具有较好的有效性。细观模型在不同应变率下的破坏模式如图18所示，应变率较低时，宏观裂纹发展不明显，结构无明显破坏，随着应变率的增大，试样表面逐渐出现宏观裂纹，且裂纹逐渐贯穿，导致试样发生破坏。

图 17 不同应变率下的动态增强因子Fig.17 Dynamic increasing factor for concrete at different strain rates

图 18 不同应变率下混凝土试样破坏模式Fig.18 Failure patterns at differentstrain rates

4 结 论

在传统Voronoi 图形的基础上，提出了一种简单高效的混凝土三维细观模型建立方法：即通过引入缩放因子产生具有一定级配的随机凸多面体作为骨料，随后对多面体的外表面进行延伸，得到具有一定厚度的黏接界面层，最后通过布尔运算得到包裹骨料和界面层的砂浆部分。利用塑性损伤模型对该细观模型进行准静态和动态加载下的有限元分析，得到如下结论。

（1）利用本文的三维细观模型模拟混凝土单轴、双轴以及三轴应力状态下的静力学性能，得到的应力应变曲线和破坏模式与实验结果基本一致，说明该模型能较好地反映实际混凝土在准静态载荷下的力学特性，这为进一步从细观角度研究混凝土损伤演化规律和破坏机理，提供了模型基础。

（2）混凝土在单轴压缩下的破坏模式与接触面摩擦条件有较大关系，在摩擦较低的条件下，主要产生与受载方向平行的竖向裂纹；在高摩擦条件下，破坏形态为八字型裂缝。拉应力作用下，主要产生垂直于载荷平面的环状裂纹。

（3）动态SHPB 数值模拟结果表明，混凝土的抗压强度随着应变率的增大而提高。在较低应变率（50 s-1以下）下，试样外表面未发生明显损伤，随着应变率的提高，试样损伤愈发明显，外表面产生明显裂纹，导致混凝土失效破坏。