混凝土破坏过程的复合型界面损伤模型与数值模拟

刘智光，陈健云

(大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连，116024)

将混凝土看作宏观均质材料，根据混凝土的变形特点，人们提出了许多宏观断裂模型。在宏观尺度下，忽略混凝土的内部结构，假定为均质材料，反映了一种工程平均，是工程设计所必需的。但如果研究混凝土材料断裂过程中的损伤演化问题和非线性力学行为产生的物理机制，材料细观结构及其力学性质的非均匀性不可忽略。因此，描述混凝土材料断裂过程的一种更为基本的方法是从其细观组成结构进行研究。

自Roelfstra等[1]提出“数值混凝土”的概念以来，基于材料细观层次上的结构特征和简单本构关系，国内外学者相继发展了一系列模拟混凝土断裂的数值模型。在这些模型中，混凝土通常被视为由骨料、砂浆基质及其之间的界面过渡区(ITZ)组成的复合材料。三者的材料力学性能差异很大，如界面过渡区的渗透性明显高于砂浆基质的渗透性，而其弹性模量和强度比砂浆基质低的。相关试验和研究表明，混凝土性能(如强度、弹性模量、断裂性能)很大程度上与界面过渡区的几何和物理性能有关[2−3]。界面过渡区被认为是混凝土的薄弱环节，普通混凝土中裂纹通常在骨料边缘形成的界面过渡区最先出现。然而，由于界面过渡区的厚度仅为0.01～0.10 mm[4]，一般不能直接进行网格剖分，难于在数值模型中直接考虑。因此，如何合理反映界面过渡区的几何形态及其性质对混凝土材料的影响，是细观数值模拟混凝土断裂过程相关研究的重要内容。

Wriggers等[5]在数值模型中不考虑界面过渡区，将混凝土视为由骨料、砂浆基质组成的两相复合材料，忽略了界面材料对混凝土整体性能的影响。Schlangen等[6−8]将界面过渡区的厚度增大以适应网格剖分，然后相应地改变其弹性模量和强度，但是如何改变并没有明确的关系。Cusatis等[9−10]所采用的格构模型以各相材料串联的方式考虑界面过渡区的弹性模量，而界面过渡区的强度并不能很好地反映。一些学者将界面过渡区作为独立的材料组分，采用Goodman单元、零厚度界面单元等独立界面单元技术[11−12]，应用于混凝土破坏过程模拟取得了较好的效果，但是混凝土材料中界面过渡区大量和无序的存在，往往需要进行烦琐的前处理工作，而且有些界面单元的本构关系稍显复杂。于庆磊等[13]基于数字图像的手段表征混凝土的三相细观结构，在数值模型中可以反映界面过渡区真实形状与分布，但目前还难以严格从界面过渡区的尺寸量级上进行分析。

基于此，本文作者建立了复合型界面损伤模型应用于混凝土破坏过程的细观数值模拟。该模型对多边形随机骨料试件进行规则化网格剖分，将内嵌界面过渡区的单元视为一种广义复合材料单元，采用修正的Vogit-Reuss模型[14]描述单元的平均材料性质，单元的损伤通过各组成材料的(弹脆性)损伤表现。本文作者采用材料参数的随机分布考虑混凝土材料力学性能的细观不均匀性[7]。运用该模型，对单轴拉伸和压缩加载形式下的混凝土破坏过程进行了数值模拟。

1 复合型界面损伤模型

图1所示为混凝土随机骨料结构、规则化有限元网格与细观单元。参照从三维Fuller骨料级配曲线到二维平面骨料级配问题的Walraven公式[15]，可以生成二维混凝土多边形随机骨料试件(图1(a))。将规则化有限元网格映射到试件上形成如图1(b)所示的细观单元，包括完全位于骨料内或砂浆基质中的单一材料单元，以及位于骨料与砂浆基质之间、包含界面过渡区的复合材料单元(图1(c)，文中称为内嵌界面材料单元)。细观数值模拟中通常的做法是将界面过渡区的材料性能和本构关系赋予复合材料单元。但是有限元网格尺寸上往往数倍于界面过程区厚度(0.01～0.10 mm)，界面材料在复合材料单元中所占比例很小，如仅以界面过渡区的材料力学性能描述复合材料单元，显然“削弱”了单元性能，界面过渡区发生损伤后更是如此。

图1 混凝土随机骨料结构、规则化有限元网格与细观单元Fig.1 Random aggregate structure of concrete,regular FE mesh and mesoscopic element

为了在细观数值模拟中考虑界面过渡区材料，应用复合材料力学中的匀质化技巧，建立了复合型界面损伤模型，该模型从2方面对内嵌界面材料单元进行了复合：(1)采用修正Vogit-Reuss模型确定内嵌界面材料单元的平均材料性质；(2)单元的损伤通过各组成材料的损伤体现。

1.1 内嵌界面材料单元的平均材料性质

考虑图2所示的标准型内嵌界面材料单元，假设各相材料的应力和应变是均匀的，，和分别为骨料、砂浆基质和界面过渡区材料的应力(应变)张量。

图2 标准型内嵌界面材料单元Fig.2 Standard interface embedded element

和为内嵌界面材料单元的平均应力和平均应变张量，根据修正的Vogit-Reuss模型有：

式中：la，lb和lm为单元中各相材料的组成比例。各相材料均为线弹性各向同性材料，以Cijkl为表示材料的柔度张量，它们的本构关系可表述为：

由式(1)～(9)可以得到

其中：

式中：Ea(Em,Eb)，va(vm,vb)和Ga(Gm,Gb)分别为骨料(砂浆基质、界面过渡区)的弹性模量、泊松比和剪切模量。若界面过渡区材料向单元的边倾斜(图1(c))，则保持材料之间的比例不变，在倾斜界面确定的局部坐标系下求得内嵌界面材料单元的复合弹性矩阵，然后将其转换到整体坐标系。由上述分析可知：(1)内嵌界面材料单元本构关系由各组分的材料力学参数、组成比例以及界面的方向确定；(2)在细观层次上，尽管各组成材料是各向同性，但内嵌界面材料单元表现出各向异性。

此外，若已知内嵌界面材料单元的平均应力、平均应变，利用式(1)～(9)可得到各相材料的局部应力和局部应变。

需要说明的是，骨料形状和位置是随机的，规则化网格剖分会产生少量其他形式的包含多条界面过渡区的单元，这时在程序实现过程中做一些简化处理，只考虑包含一条界面材料的内嵌界面材料单元。

1.2 内嵌界面材料单元的损伤

内嵌界面材料单元的损伤通过其组成材料的依次损伤体现，采用弹脆性损伤本构关系描述材料的细观力学行为，且各组成材料的损伤计算模型是一致的。

采用拉断的Mohr-Coulomb准则作为损伤判据。当材料的应力状态或者应变状态满足损伤阈值时，材料开始损伤，在此之前材料为线弹性。假设材料在单轴受力状态下满足如图3所示的损伤本构关系，图中σ为应力，ε为应变，以拉为正，压为负。

单轴拉伸作用下以理想弹脆性损伤模型描述材料的力学行为，当材料的最大拉应力超过其拉伸强度时发生拉伸破坏而失去承载力，拉伸损伤演化方程如下：

式中：D为损伤变量，εt0为单轴拉伸损伤阈值。单轴拉伸损伤判据为：

图3 单轴载荷作用下弹性损伤本构关系Fig.3 Elastic damage constitutive law of material under uniaxial stress state

式中：1σ为最大主应力，ft为单轴抗拉强度，满足为材料初始弹性模量。

当材料的应力状态满足Mohr-Coulomb准则时发生剪切损伤，即

式中：φ为细观材料的摩擦角；fc为单轴抗压强度(正数)；3σ为最小主应力。单轴受压状态下损伤变量表达式为：

式中：εc0为材料达到其单轴抗压强度时对应的最大压缩主应变，εc0=−fcE0；λ为单元的残余强度系数；fcr为材料的残余抗压强度(正数)，有fcr=λfc。

按照等效应变原理，材料损伤后的应力−应变关系为：

在程序实现中，内嵌界面材料单元的损伤计算可按如下方式进行：得到各组成材料的应力(应变)后计算各组成材料的损伤，而单元的损伤体现在其复合本构关系的弱化中。为避免问题的复杂性，假设内嵌界面材料单元内各组成材料的损伤是各向同性的，但是得到的单元损伤显然是各向异性的。

1.3 单一材料单元的损伤

对于只含骨料和砂浆基质的单一材料单元也采用上述弹脆性损伤本构关系，只是如果单元发生拉伸损伤，则考虑细观单元的正交各向异性损伤，而如果单元发生剪切损伤，其损伤是各向同性的。

当未损伤单元最大拉应力达到拉伸损伤阈值时，在最大拉应力方向发生损伤，材料的主轴之一指向该方向，另一材料主轴方向上弹性模量不变，但剪切模量发生相应的折减；之后单元的应力状态简化为平面应力问题的单轴状态，一旦简化的应力状态满足拉伸损伤准则，则另一轴的损伤就会发生。材料正交异性损伤的本构关系矩阵为

其中：d1=(1−D1)，d2=(1−D2)，D1和D2分别为 2 种材料主轴上的损伤值。损伤后材料主轴与整体坐标系不重合时需将转换到整体坐标系。

2 单轴拉伸和压缩破坏过程模拟

考虑平面应力问题，采用边长为100 mm的正方形混凝土多边形随机骨料结构数值试件(断面骨料填充率为 45%，最大骨料粒径为 10 mm)，界面厚度取0.1 mm，生成100×100个同一尺寸的正方形单元，包括骨料单元、砂浆基质单元和内嵌界面材料单元。这里考虑抗压强度为30.0 MPa的混凝土，假定骨料、砂浆基质和界面材料的强度和弹性模量等力学参数服从 Weibull统计分布。对于材料细观力学参数统计量选定问题，人们已进行了大量的研究[7,17−18]，得到了一系列较为合理、可靠的参数(数值配比方案)，基于此本文取值如表1所示。骨料、砂浆基质和界面材料的泊松比离散性相对较小，其值分别为0.18，0.20和0.25；摩擦角均取为30°，残余强度系数为0.05，单轴抗拉强度为单轴抗压强度的1/15。试件加载过程采用位移控制，单轴拉伸加载位移步长取0.5 μm，单轴压缩为5 μm。试件底面节点为法向约束，其他自由。

表1 混凝土试件的 Weibull 分布力学参数Table 1 Mechanical parameters of Weibull distribution of concrete specimen

2.1 单轴拉伸试验

为验证数值结果的合理性，说明建议的复合型界面损伤模型有效，与文献[19]的物理试验进行了对比。需要指出的是，本文细观数值试验并不是文献[19]试验的模拟，无法与试验结果进行定量比较，但比较两者应力－应变曲线的形状，研究相似特征和变化趋势，验证本文数值模型是合适的。试件在单轴拉伸荷载作用下的应力－应变全曲线如图4所示(单轴抗拉强度为2.59 MPa，峰值点对应的拉伸应变为0.105×10−3，上升段30%峰值应力处的割线模量为34.1 GPa)。曲线包括线性上升段、非线性强化段和应变软化段，较好地反映了试验过程中混凝土材料的宏观非线性力学行为，说明建立的复合型界面损伤模型能较好地模拟混凝土的应力−应变全过程曲线。

单轴拉伸断裂过程如图5所示，图中黑色部分表示发生了损伤的单元，A，B，C，D表示图4中不同加载水平。

图4 混凝土轴拉应力−应变曲线对比Fig.4 Comparison of stress−strain curves between numerical result and laboratory experimental result of concrete specimen under uniaxial tension

图5 混凝土单轴拉伸的损伤演化过程Fig.5 Failure process of concrete under direct tensile test

试件开始加载后，损伤在试件中缓慢、均匀地发展。加载到峰值应力的70%～80%时(如点A)，大量界面过渡区发生了损伤，形成的微裂纹均布在试件中，且都处于独立的发展状态。此后试件损伤发展加快，部分微裂纹出现贯通趋势，在若干薄弱位置形成局部裂纹，应力−应变曲线可看到明显的非线性响应。到峰值应力时(点B)，在试件最薄弱位置形成包含几条局部化裂纹的主裂纹带，但此时并未形成宏观裂纹。随着宏观应变的进一步增加，新增损伤主要集中发生在主裂纹带，其他部分的微裂纹和局部裂纹发展缓慢甚至休止(点C)。在主裂纹带内，多条局部化裂纹绕过骨料、在界面和砂浆基质中扩展、桥接，最终形成一条大致垂直于拉伸荷载方向的宏观裂纹，这是众多裂纹成核过程。

在试样破坏过程中，损伤主要发生在界面过渡区，峰值应力状态后才有少量砂浆基质单元损伤；另外，由于采用材料参数的随机分布考虑了混凝土各相材料的力学性能细观不均匀性，损伤总是发生在较弱的单元材料中。图6所示为峰值应力状态和加载终止时发生损伤的界面过渡区和砂浆基质的抗拉强度直方图，同时还给出了初始未加载时抗拉强度的分布。虽然混凝土被认为是宏观均匀的，但材料细观层次的局部性能是有差异的，通过反映这种细观非均质性，数值模拟能够得到损伤局部化和材料的渐进破坏过程。

图6 不同加载状态发生损伤的界面过渡区和砂浆基质抗拉强度直方图Fig.6 Histogram of tensile strength in specimen of interfacial transition zone and matrix

2.2 单轴压缩试验

图7所示为单轴压缩加载作用下试件的标准化应力−应变全曲线，峰值应力为31.59 MPa，约为单轴拉伸时峰值应力的12倍，峰值应变为1.30×10−3。直到30%峰值应力，曲线上升段均保持线性，此时的弹性模量为34.0 GPa。之后可看到非线性响应，当应力增加到70%～80%峰值应力时，非线性更加明显。从这一应力水平直至峰值应力时，应力－应变曲线显著弯曲，随着宏观应变的增加，进入应变软化段，曲线明显下降。对比van Vliet和van Mier[20]的物理试验结果(峰值应力约为36.5 MPa，对应的应变为1.85×10−3)，除了残余强度段有差别，数值模拟得到的应力−应变曲线形状与试验结果是相似的，而且应力和应变与试验结果都在同一量级。

图7 混凝土单轴压缩应力−应变曲线(轴向)Fig.7 Axial stress−strain curves of concrete specimen under uniaxial compression

图8所示为单轴压缩轴向应力−应变曲线与轴向应力−侧向应变曲线。直至 80%峰值应力时，横向应变(ε2)随纵向应变(ε1)成比例增长，泊松比仍近似保持为常数(0.19～0.21)；之后，表观泊松比开始增加。到峰值应力时，横向应变增加越来越快，试件表观泊松比约为0.58，与试验值接近[21]；在软化段横向应变比轴向应变增加更快，表观泊松比大于 1。文中算例模拟的是加载端光滑情况，故横向变形增加越快。

图8 单轴压缩轴向应力−应变曲线与轴向应力−侧向应变曲线Fig.8 Axial stress−strain curve and axial stress−lateral strain curve of concrete specimen under uniaxial compression

图9 单轴压缩轴向应力−“虚拟”体积应变曲线Fig.9 Axial stress−“fictitious” volumetric strain curves of concrete specimen under uniaxial compression

这里将(ε1+2ε2)定义为“虚拟”体积应变，如图 9所示，从开始加载至应力为80%峰值应力处“虚拟”体积应变几乎为线性减小的。当应力达到90%峰值应力时，“虚拟”体积应变达到最小值，随后，“虚拟”体积应变的变化趋势反向，并在达到峰值应力时发生体积膨胀。因为加载端是光滑的，峰值应力后试件扩容幅度很大。

图10和11所示为试件单轴压缩断裂过程中不同加载点(图7中A，B，C，D和E)的损伤分布和变形示意图。图10中黑色部分为发生损伤的单元，白色多边形表示骨料颗粒；图11中黑色部分代表骨料。

在加载的初始阶段，直到加载到最大压应力的80%(点A)，在试件中只有少量与加载方向平行的界面过渡区材料发生损伤。随着加载的进行，损伤的单元增多，出现新的微裂纹，应力−应变曲线的斜率迅速减小；同时微裂纹不断扩展，形成多处损伤相对集中的平行于加载方向的局部裂纹。在峰值应力点B，试件中部和右部首先出现了贯通的局部裂纹。进入软化段，损伤急剧增加，大量的局部裂纹形成、彼此贯通；随着宏观应变的增加，如图11中加载点C，D和E的变形图所示，在试件中形成3条主要的宏观裂纹，顺着荷载方向(未完全)分割试件，试件最终表现为劈裂破坏。从图10中加载点E的损伤发展可以看出：大部分的单元都发生了破坏，可以解释单轴压缩破坏的试件通常都比较碎的试验现象。图12所示为单轴压缩试件形态试验结果。数值模拟得到的试件破坏形式和裂纹形态与图12的物理试验结果[22]比较，均表现出较好的相似性。

图10 单轴压缩试件的损伤演化过程Fig.10 Failure process of concrete specimen under uniaxial compression

图11 单轴压缩试件变形图Fig.11 Deformed configuration of concrete specimen under uniaxial compression

图12 单轴压缩试件坏形态试验结果Fig.12 Experiments results of specimen failure mode under uniaxial compression

就本文研究内容而言，试件在单轴拉伸和单轴压缩荷载作用下所表现的局部变形的形式和内容各有异同。在单轴拉伸荷载作用下，峰值应力后试件表现出非常明显的局部变形，在与荷载方向垂直的一段临界区域内形成一条宏观主裂纹带。随后，破坏过程中的变形主要集中在主裂纹带，在远离主裂纹带的区域为卸载，但总的变形继续发展。考虑加载端为理想无摩擦的单轴压缩荷载作用，试件发生劈裂破坏，形成 3条明显的宏观劈裂裂纹，与单轴拉伸时一样，也是由材料的拉伸破坏造成的；峰值应力后在3条宏观劈裂裂纹带出现局部变形。因此，这进一步印证了单轴试验中得到的、包含应变软化段的(名义)应力−应变曲线并不能作为纯粹的材料行为[9,19]。同时也表明了通过反映混凝土材料细观结构的随机性和非均匀性，复合型界面损伤模型具有反映材料破坏过程中的局部化变形的内蕴能力，而仅通过改进材料本构关系是难以做到的。

虽然数值模型考虑了受拉损伤和剪切损伤2种基本的机制，但是单轴拉伸和单轴压缩破坏过程表明，在多数情况下，材料发生拉伸损伤而破坏，剪切损伤机制参与工作比例很小。

3 结论

(1)建立复合型界面损伤模型模拟混凝土材料的非线性力学行为和断裂过程。与以往数值模型最大的不同在于：将界面过渡区作为单元的内嵌组分材料，利用复合材料力学中的匀质化技巧，充分考虑了混凝土材料的细观结构特征和材料性能的细观非均质性。

(2)该模型可以较好地模拟混凝土试件的单轴拉伸和压缩破坏特征，包括断裂行为、峰后软化和局部化变形，特别是得到加载端无摩擦轴压破坏的劈裂形式。

(3)基于复合型界面损伤模型的细观数值模拟适应于规则化有限元网格，避免了直接对界面过渡区进行剖分，既简化了前处理工作、减少了单元数量，提高了计算效率，又简单而很合理地将界面过渡区反映在混凝土材料断裂过程的模拟中。因此，该模型适合应用于动态荷载作用下混凝土断裂的数值模拟和处理三维问题。

[1]Roelfstra P E,Sadouki H,Wittmann F H.The numerical concrete[J].Materials and Structures,1985,18(5): 327−335.

[2]应宗权,杜成斌.考虑界面影响的混凝土弹性模量的数值预测[J].工程力学,2008,25(8): 92−96.YING Zong-quan,DU Cheng-bin.A numerical method for effective elastic modulus of concrete with interfacial transition zone[J].Engineering Mechanics,2008,25(8): 92−96.

[3]Rao G A,Prasad B K R.Influence of the roughness of aggregate surface on the interface bond strength[J].Cement and Concrete Research,2002,32(2): 253−257.

[4]Bebtur A.Microstructure,interfacial effects and micromechanics of cementitious composites[J].Ceram Trans,1990,16: 523−550.

[5]Wriggers P,Moftah S O.Mesoscale models for concrete:Homogenisation and damage behavior[J].Finite Elements in Analysis and Design,2006,42(7): 623−636.

[6]Schlangen E,van Mier J G M.Simple lattice model for numerical simulation of fracture of concrete materials and structures[J].Materials and Structures,1992,25: 534−542.

[7]唐春安,朱万成.混凝土损伤与断裂数值试验[M].北京: 科学出版社,2003: 34−44.TANG Chun-an，ZHU Wan-cheng.Numerical test of damage and fracture of concrete[M].Beijing: Science Press,2003:34−44.

[8]刘光廷,王宗敏.用随机骨料模型模拟混凝土材料的断裂[J].清华大学学报,1996,36(1): 84−89.LIU Guang-ting,WANG Zong-min.Simulation of the fracture of concrete with random aggregate model[J].Journal of Tsinghua University,1996,36(1): 84−89.

[9]Cusatis G,Bazant Z P,Cedolin L.Confinement-shear lattice model for concrete damage in tension and compression (Ⅱ):Computation and validation[J].Journal of Engineering Mechanics,2003,129(12): 1449−1458.

[10]Lilliu G,van Mier J G M.3D lattice type fracture model for concrete[J].Engineering Fracture Mechanics,2003,70(7/8):927−941.

[11]Kwan A K H,Wang Z M Chan H C.Mesoscopic study of concrete Ⅱ: Nonlinear finite element analysis[J].Computer and Structures,1999,70(5): 545−556.

[12]Caballero A,Lopez C M,Carol I.3D meso-structural analysis of concrete specimens under uniaxial tension[J].Computer Methods in Applied Mechanics and Engineering,2006,195(52):7182−7195.

[13]于庆磊,唐春安,朱万成,等.基于数字图像的混凝土破坏过程的数值模拟[J].工程力学,2008,25(9): 72−78.YU Qing-lei,TANG Chun-an,ZHU Wan-cheng.et al.Digital image-based numerical simulation on failure process of concrete[J].Engineering Mechanics,2008,25(9): 72−78.

[14]陈惠发,Saleeb A F.混凝土和土的本构方程[M].于天庆,王勋文,刘西拉,等,译.北京: 中国建筑工业出版社,2004:279−286.CHEN Hui-fa,Saleeb A F.Constitutive equations for materials of concrete and soil[M].YU Tian-qing,WANG Xun-wen,LIU Xi-la,et al,trans.Beijing: China Architecture and Building Press,2004: 279−286.

[15]Walraven J C.Aggregate Interlock: A theoretical and experimental Analysis[D].Delft: Delft University of Technology.Faculty of Civil Engineering and Geosciences,1980: 37−62.

[16]江见鲸,陆新征,叶列平.混凝土结构有限元分析[M].北京：清华大学出版社，2005: 232−233.JIANG Jian-jing,LU Xin-zheng,YE Lie-ping.Finite Element Analysis of Concrete Structures[M].Beijing: Tsinghua University Press,2005: 232−233.

[17]方志,杨钻,苏捷.混凝土受压性能的非均质细观数值模拟[J].湖南大学学报: 自然科学版,2010,37(3): 1−6.FANG Zhi,YANG Zuan,SU Jie.Mesoscopic numerical simulation on the compressive behavior of Heterogeneous Concrete[J].Journal of Hunan University: Natural Sciences,2010,37(3): 1−6.

[18]Zhu W C,Tang C A.Numerical simulation on shear fracture process of concrete using mesoscopic mechanical model[J].Construction and Building Materials,2002,16(8): 453−463.

[19]刘西拉,温斌.考虑广义边界条件的混凝土软化[J].岩石力学与工程学报,2008,27(5): 885−892.LIU Xi-la,WEN Bin.Concrete strain softening considering global boundary conditions[J].Chinese Journal of Rock Mechanics and Engineering,2008,27(5): 885−892.

[20]van Vliet M R A,van Mier J G M.Experimental investigation of concrete under uniaxial compression[J].Mechanics of Cohesive-frictional Materials,1996,1(1): 115−127.

[21]van Mier J G M,Shah S P,Arnaud M,et al.Strain-softening of concrete in uniaxial compression[J]Materials and Structures,1997,30(4):195−209.

[22]van Vliet M R A,van Mier J G M.Softening behavior of concrete under uniaxial compression[C]//Wittmann F H.Fracture Mechanics of Concrete Structures.Freiburg: Aedificatio,1995:383−396.