方建银，党发宁

(西安理工大学 岩土工程研究所，陕西 西安 710048)

岩石、混凝土等脆性材料大量应用于土木工程建设，其力学特性受到广泛关注和研究，特别是当前水利大坝建设所需的大体积混凝土，其静、动力学特性更是研究中的重点和难点。当前利用物理试验对其进行的研究，往往由于试验设备的局限性，达不到理想的预期结果。“数值试验”具有灵活性高、可重复操作、可模拟物理试验不能进行的试验以及研究成本低等优点，倍受科研人员的青睐。而如何能够合理地利用数值试验来研究岩石、混凝土等非均匀脆性材料的力学特性是许多研究者亟待解决的难题。为此，以细观层次为基础的数值模拟方法应运而生。目前的研究大多仅从二维和三维模型建立[1-3]、裂纹演化规律[4-5]及尺寸效应[6-7]等方面进行了研究，而三维随机骨料模型的建立、加载方式和试样尺寸效应对脆性材料特性、损伤破裂的影响及脆性材料数值强度定义方面研究较少。

本文基于损伤力学原理，利用均质岩样研究了加载方式和试样尺寸对材料特性和裂纹的影响，确定了适合研究脆性材料强度及裂纹演化的数值方法。将此方法推广到细观混凝土的强度及裂纹演化分析中，从细观层面上研究了混凝土的强度与破裂特性，并利用CT试验验证了本文所确定的数值试验方法的合理性。

1 细观损伤演化模型及试验条件

本文选用弹性损伤本构关系来描述混凝土类脆性材料的力学特性。依据Lemaitre等价应变原理可以将应力应变关系表示为σ=Eε，其中E=E0(1-D)是材料受损后的弹性模量，E0是材料初始弹性模量，D是损伤变量。

对于混凝土类脆性材料，无论是拉破坏还是压破坏，归根到底都是由于材料内部单元所承受的拉应力超出了其抗拉极限强度。因此本文采用最大拉应变强度准则，即当细观材料单元的最大主拉应变值超出了给定的极限拉应变阈值时单元开始损伤，此时单元仍然具有一定的刚度和承载力，损伤演化模型见式(1)。

(1)

式中：λ为强度残余系数；ε0为抗拉强度对应的主拉应变；εr=ηε0为抗拉残余应变；η为残余应变系数；εu=ζε0(ζ>η)为极限拉应变；ζ为极限应变系数；εmax为加载历史上主拉应变的最大值。

对于损伤演化方程(1)中的参数的取值，根据文献[8]，取λ=0.1，η=5，ζ=10。当单元的最大拉应变达到εu时认为单元完全破坏，将该单元杀死。

1.1 数值模型及约束条件

三维混凝土随机骨料模型是从细观层面出发，将混凝土看成是由骨料、砂浆及二者之间的结合面所组成的三相复合材料，在Ansys软件中采用蒙特卡罗(Monte Carlo)方法[9]，利用随机变量得到骨料位置，据此建立起三维“数字混凝土”模型。然后运用骨料投影网格法，通过“多次划分单元”技术来判别出骨料、砂浆及界面单元。并分别对各组分材料赋予弹性模量等参数，具体建模请参考文献[10] 。

均质岩样及数字混凝土试样均采用120 mm×60 mm(高×直径)模型。其中数字混凝土试样如图1所示。约束采用试样底面中心点全约束，底面其它节点采用垂向约束，试样顶面节点采用垂向约束。

图1 混凝土三相材料有限元模型

1.2 数值试验参数

为了保证数值试验与真实试验的一致性，采用了西北水利水电勘测设计研究院测定的混凝土三相组成材料的力学参数(见表1)，均质岩样选用骨料参数。

表1 混凝土各组分材料参数[1]

2 数值试验的加载方法及试样尺寸选择

2.1 加载方式的选择

力学试验和数值试验均可采用以应力控制的柔性加载和以位移控制的刚性加载方式。不同加载方式往往得出不同的试验结果，对于数值试验研究更是如此。由于可以解析计算均值岩石试样单轴拉压时的应力及位移，为此，本文首先针对均值岩石试样进行了3个应力控制和3个位移控制的单轴压缩试验，以研究用数值试验恰当地确定材料强度的方法，试验结果如图2所示。

由两种加载方式下的试验结果可见，应力单调增加法的优点是荷载易于控制，曲线光滑；不足是应力应变曲线单调增加，位移突变点不明显，数值计算中为柔性加载，导致试样顶面各节点位移不相同，与试验室试样顶面位移相同的边界条件不符合。位移控制加载法的优点是应力应变曲线与试验室应力应变全曲线相似，峰值点明显，数值计算中为刚性加载，试样顶面各节点位移相同，与试验室试样顶面位移相同的边界条件相符合；不足是与力学试验习惯不同。应力控制加载的数值试验中，应力只能单调递增而不能衰减；而位移控制加载的试验结果中应力的变化更符合实际。

图2 不同加载方式下的单轴受压试样应力应变曲线

对比试验结果并考虑数值试验与力学试验结果的一致性，以下研究都采用以位移控制的加载方式，并将计算的应力应变曲线的峰值点作为材料的极限强度点。

2.2 试样尺寸的选择

脆性材料力学特性的尺寸效应受到人们广泛关注。为了验证位移控制加载的损伤模型数值试验方法中试件的尺寸效应，利用60 mm×120 mm、120 mm×240 mm和180 mm×360 mm三种尺寸的均质试样，研究了长细比相同而直径不同对试样强度的影响，计算的应力位移曲线和应力应变曲线分别如图3和图4所示。

由图3可以看出，三种试样的应力位移曲线不论位移加载还是应力加载，在弹性阶段和屈服的初期阶段均是完全重合的，只有屈服的后期段由于加载方式不同而导致残余段曲线不重合，原因是应力控制加载的数值试验中，应力只能单调递增而不能衰减。

图3 不同试样应力位移曲线

图4 不同试样应力应变曲线

由图4可看出，数值试验中在相同的长细比下，不论是应力控制加载还是位移控制加载，三种试样的应力应变曲线都几乎重合，即不论采用哪种试样其弹性模量不会发生变化，同样屈服强度也都不会因为试样的大小而不同。即对于均质岩样，长细比一定的情况下数值试验中不考虑其尺寸效应的影响。这与文献[7] 基于最弱链模型和缺陷的Poisson 分布假设建立的准脆性材料统计模型所得的结论是一致的。

本文以下采用长细比为2的60 mm×120 mm试样开展研究。

3 混凝土强度及破坏形态数值试验研究

3.1 损伤区的划分

为了能定量地进行分析研究，就需要对材料进行分区定义。根据式(1)，定义试样各单元的最大拉应变εmax<ε0为弹性区(此时单元未发生损伤);ε0≤εmax≤εr为第一损伤区(此时单元开始损伤)；εr<εmax≤εu为第二损伤区(此时单元发生了比较大的损伤)；εmax>εu为第三损伤区(此时单元主拉应变已经超出了材料的最大拉应变极限值，单元已经破坏)。以下分析直接应用第一损伤区、第二损伤区和第三损伤区进行描述。

3.2 “数字混凝土”强度计算方法研究

从细观力学角度将混凝土看成由骨料、砂浆和界面组成的复合材料，此时的混凝土是一种结构而非一种材料，其整体强度可以由数值计算得到。将利用均质岩样确定的方法应用于“数字混凝土”，通过试样的损伤数值模拟试验可以得出图5所示的单轴拉、压受力状态下的应力应变曲线。

由图5可看出，数值试验所得应力应变曲线不论是拉还是压都与常规力学试验的应力应变全曲线[11]很相似，它们同样有上升段和下降段(既有弹性阶段和损伤屈服阶段，也有残余变形阶段)。

因此，“数字混凝土”试验也采用应力应变曲线的极值点C点作为复合材料的极限强度点。据此可得本次“数字混凝土”试验的抗压强度为7.32 MPa，抗拉强度为1.69 MPa。

图5 应力-应变全曲线

3.3 混凝土损伤特性的细观力学研究

由损伤数值试验可得图6所示的单轴压、拉状态下试样损伤单元柱状图。

图6 单元损伤柱状图

对比图6(a)、(b)可以看出，在单轴压缩试验中，前三步荷载时不论是界面还是砂浆均没有单元损伤，说明试样处于完全弹性阶段，这与图5的OA段相对应；而随着荷载步的增加，界面单元首先开始出现损伤且均为第一损伤区，这一特性与图5的AB段相对应，即试样出现非线性阶段。随着荷载的不断增加，试样界面单元损伤数量急剧增加，当达到第6步荷载时，有大量界面单元进入第二损伤区，且有少量单元损伤破坏，同时此阶段有砂浆单元进入第一损伤区，且随着荷载增大，界面和砂浆进入第一损伤区的单元数也在不断增加。

随着荷载进一步增加，当达到第7步荷载时，开始有试样界面单元进入第二损伤区，而第9步荷载正好是试样的极限荷载步，这与图5的BC段相对应。从第9步荷载开始一直到加载完成，随着荷载的增加，界面第二损伤区的单元在增加的同时开始向损伤破坏区转变，破坏单元急剧增加，而砂浆单元大部分还处于第一损伤区，只有少量单元进入第二、第三损伤区，即只有少量砂浆单元破坏，这与图5的CE段对应。

另外，可以发现，整个加载过程中界面单元的损伤数量及被杀死数量均比砂浆单元数量多，这说明混凝土试件在静压破坏时破裂面追随试件薄弱面(界面)发展。对于单轴静拉试验如图6的(c)、(d)两图所示，其规律和单轴压缩试验一样。说明不论是单轴拉伸还是单轴压缩，试样破坏裂纹都是追随试样中薄弱的界面而发展的，裂纹遵从萌生、发展、贯通的发展过程，其裂纹演化如图8和图9的(a′)、(b′)、(c′)所示。

3.4 “数字混凝土”强度与损伤单元关系的分析

为了更进一步分析强度与损伤面积的关系，可以通过荷载与损伤单元面积的关系曲线来研究。单元总的损伤面积为总的损伤单元个数乘以一个单元的面积，而一个单元的面积仅为0.219 8 mm2，为此可以用荷载与损伤单元总数的关系曲线来代替荷载与损伤单元面积关系曲线进行分析(如图7所示)。

图7 荷载-总损伤单元数关系曲线

由图7可以看出，不论是单轴压缩还是单轴拉伸，随着位移加载的增大，弹性阶段都没有单元发生损伤，随着加载进一步增大，开始有单元发生损伤，应力到峰值强度时损伤单元数量并未达到最大值，而是随着位移加载继续增多，但增加幅度越来越小。当增大到一定值时，出现了损伤单元数随着加载开始减少，直到加载停止。深入分析可以发现，这种现象主要由两方面因素造成。一方面，从开始加载到损伤单元数达到最大这个过程中，随着位移加载的增加试件中大量单元的应变超出了其弹性应变极值而发生损伤，所以损伤单元一直处于增多的趋势。但是，这一过程中，随着位移加载的增大有少量单元的应变超出了其极限应变发生破坏，并释放了少量的应变能，这致使损伤单元数量的增加幅度在这一过程末尾时变缓，这时试样中产生了微裂纹；另一方面，随着位移加载的增加，有越来越多单元应变超出了其极限应变发生损伤破坏，微裂纹相互贯通形成宏观裂纹，试样这时释放出了大量的应变能，致使整体损伤单元数量出现了减少的趋势。

4 数值试验与CT试验对比研究

为了验证位移控制加载的混凝土损伤模型计算结果的合理性，采用单轴拉压状态下的CT试验结果对比研究混凝土破坏裂纹的演化过程。试验采用西门子SOMATOM-plus医用CT扫描仪和西安理工大学开发的CT专用动态三轴仪，CT扫描分辨率为0.356 mm×0.35 mm×1 mm，采用应变控制。试样规格为60 mm×120 mm的一级配C15混凝土圆柱体。

根据CT 机自带的图像重建程序建立的CT试验结果如图8及图9的(a)～(f)所示，由于扫描图片的数量很大，所以本文只选取一个扫描断面有代表性的三次扫描结果和数值试验中对应的三步加载试验结果进行对比分析。

图8 受压试件不同荷载步损伤截面图

图9 受拉试件不同荷载步损伤截面对比图

分析图8及图9可得出以下几点。

1)CT实时扫描作为一种无损探测技术，可以直观地发现混凝土试样内部的裂纹萌生发展情况，探测发现试样在单轴拉压情况下首先界面出现肉眼不可见的损伤微裂纹，随着荷载的增大，裂纹绕过骨料发展贯通，进而形成贯通的裂纹面，反映了静载荷作用下试样界面是其薄弱环节，裂纹追随薄弱环节发展，由于骨料强度高，可以阻碍微裂纹的发展，进而提高混凝土的强度；当荷载达到峰值时，试件中形成的微裂纹不明显，肉眼很难观察到，随后随着应变增大，试件进入残余强度变形阶段，裂纹开始分叉、伸长、大量的微裂纹贯通形成一条主裂纹，肉眼可以直观地分辨出来，表明试件失稳并发生破坏。单轴压缩试验试样破裂面大体与试样成45°，且破裂面积较大，破裂较彻底，而单轴拉伸试样破裂面几乎垂直于试样，只有一道贯通的横向裂纹，由于压试件的裂纹面积远大于拉试件的裂纹面积，根据断裂力学形成单位裂纹面所需的能量是常数这一规律，因此混凝土的单轴压缩试验强度远大于单轴拉伸试验强度，可以从一个侧面解释混凝土压强度高于拉强度的根源所在。

2)“数字混凝土”拉压试样损伤截面如图8和图9的(a′)～(f′)所示，可以看出不论是单轴压缩还是单轴拉伸，界面单元最先开始出现损伤，进而导致整个试件应力重分布，在它们周围容易形成应力集中，促使周围单元发生集中拉应力，产生拉伸损伤，导致单元进一步的损伤破坏；其次是砂浆单元发生损伤进而破坏，骨料单元由于强度大不破坏；另外，随着荷载逐级增大，当超过峰值荷载后，试样裂纹开始发展、贯通，最终形成贯通的破裂面，且都同样追随最薄弱的界面。对于单轴压缩，裂纹大体与试样成45°角，而对于单轴拉伸则形成与试样垂直的单一水平裂纹。

3)CT试验结果与数值模拟的结果在裂纹的萌生、扩展过程上有一定的一致性，从统计意义上讲，其破坏规律是基本相同的。

5 结 论

1)通过数值试验比较，得出本文建立的混凝土三维随机骨料模型、确定的位移控制加载方法及采用的损伤本构模型，可以较好地模拟混凝土的损伤演化特性。

2)可以用数值试验位移控制的应力应变曲线极值点作为数字混凝土材料的强度点。

3)不论是在拉还是压荷载作用下，混凝土裂纹均是先从相对较弱的界面层开始萌生，然后微裂纹绕着骨料扩展、贯通，这与CT实时扫描试验结果具有较好的一致性，说明混凝土的静态裂纹追随结构的弱面发展。

4)通过混凝土单轴拉压试验裂缝形态的比较发现，单轴压缩试验试样破坏彻底，裂纹面积大；单轴拉伸试验试样破坏简单，裂纹面积小，这与物理试验试样破坏形式相同。说明了细观混凝土试样强度取决于其受力状态，应力状态决定了试样的破坏形态，而破坏形态决定了试样的破坏面积，破坏面积的大小正比于试样的强度。初步解释了压强度高于拉强度的实质是破坏时破坏面积不同造成的。

参考文献：

[1] 唐春安，朱万成.混凝土损伤与断裂数值试验[M].北京：科学出版社，2003.

[2] 宋玉普.多种混凝土材料的本构关系和破坏准则[M].北京：中国水利水电出版社，2002.

[3] Mohamed A R, Hansen W. Micromechanical modelling of concrete response under static loading-Part 1: Model development and validation [J] . ACI Materials Journal, 1999, 96(2):196-203.

[4] 梁昕宇.数值混凝土模型的改进及其细观力学特性研究[D].西安：西安理工大学，2008.

Liang Xinyu.The improvement of numerical model on concrete meso-mechanics and the study of its dynamic properties[D].Xi’an: Xi’an University of Technology，2008.

[5] 杨强，张浩，周维垣.基于格构模型的岩石类材料破坏过程的数值模拟[J].水利学报，2002，(4):46-50.

Yang Qiang, Zhang Hao, Zhou Weiyuan. Lattice model for simulating failure process of rock[J].Journal of Hydraulic Engineering, 2002,(4):46-50.

[6] 黎保琨，彭一江.碾压混凝土试件细观损伤断裂的强度与尺寸效应分析[J].华北水利水电学院学报，2001，22(3)：50-53.

Li Baokun, Peng Yijiang. Strength and size effect analyses of meso-level damage fracture for RCC specimens[J].Journal of North China Institute of Water Conservancy and Hydroelectric Power,2001,22(3):50-53.

[7] 张明，卢裕杰，杨强.准脆性材料的破坏概率与强度尺寸效应[J].岩石力学与工程学报.2010,29(9)：1782-1789.

Zhang Ming, Lu Yujie, Yang Qiang. Failure Probability and strength size effect od quasi-brittle materials[J].Chinese Journal of Rock Mechanics and Engineering,2010,29(9):1782-1789.

[8] 马怀发，陈厚群，黎保琨.混凝土试件细观结构的数值模拟[J].水利学报，2004，(10)：27-35.

Ma Huaifa，Chen Houqun, Li Baokun. Meso-structure numerical simulation of concrete specimens[J].Journal of Hydraulic Engineering, 2004,(10):27-35.

[9] Metroplis N, UlamS. Monte Carlo method[J].J Amer States，1949，(44):335-341.

[10] 党发宁,韩文涛,郑娅娜, 等.混凝土破裂过程的三维数值模型[J].计算力学学报，2007，24(6)：829-833.

Dang Faning, Han Wentao, Zheng Yana, et al. 3D numerical simulation of failure process of concrete[J].Chinese Journal of Computational Mechanics,2007,24(6):829-833.

[11] 车宏亚，颜德亘，程文瀼.混凝土结构[M].北京：中国建筑工业出版社，1998.