武 亮，王 菁，糜凯华

（昆明理工大学 电力工程学院，云南 昆明 650500）

混凝土在细观层面上是由砂浆基质、粗骨料和两者间的黏结界面组成的三相复合材料.混凝土材料的各组分对其宏观力学性能均有直接影响.全级配混凝土因骨料粒径较大，致使其裂缝扩展路径与普通混凝土有较大差异.由于以宏观混凝土材料进行力学性能试验，需受试验条件、环境条件及人为因素等影响，所得到的试验结果离散性较高.为此，诸多研究者提出了多种细观模型来研究混凝土材料的力学特性，如基于有限元理论提出的格构模型［1］、MH 细观模型［2］、随机粒子模型［3］、刘光廷模型［4］以及MFPA 细观模型［5］等；基于离散元理论提出的颗粒－界面元模 型［6］和梁－颗粒模型［7］等.这些模型弥补了混凝土宏观试验分析上的很多缺陷和不足.因此，在细观层面上建立全级配混凝土仿真模型对其裂缝扩展规律全过程进行模拟计算，特别是对全级配混凝土大尺寸试件有重要意义.

随着大容量高速计算机的出现及计算机技术的发展，许多研究者利用各种编程语言和一些商业软件完成了混凝土骨料的随机投放和数值模拟.如基于传统的V－B、C语言和Fortran语言编写的随机骨料程序［8－10］；基于CAD 软件下开发的Autolisp语言编写的混凝土骨料模型［11］和用MATLAB 编写的混凝土数值模拟模型［12－14］.然而现有基于MATLAB建立的随机骨料模型通用性较差，且随机骨料的投放量不能满足全级配混凝土的要求，也不能直接导入其他有限元软件进行分析计算.为此，本研究将有限元软件COMSOL 与MATLAB 软件相结合，利用COMSOL 软件自带的功能简化了全级配混凝土细观几何模型的生成程序，高效完成了有限元网格模型的自动生成.

1 随机骨料模型的自动生成

1.1 骨料粒径和投放量的确定

混凝土粗骨料按表面形状，可分为卵石、碎石和人工破碎卵石.

为了使混凝土细观分析能在二维平面上进行，Walraven等［15］在球形骨料与等概率空间分布假定的基础上，建立了混凝土试件空间内骨料含量与其内截面的骨料面积之间的关系，即二维截面试件内任一点骨料直径D＜D0（筛孔直径）的概率Pc为：

式中：Pk为骨料（包括粗骨料和细骨料）体积占混凝土总体积的百分比；Dmax为最大骨料颗粒粒径；D 的上标均表示指数.

1.2 骨料库的生成

关于凸多边形骨料的形成，文献［10］先随机投放三角形或四边形骨料，再进行随机延凸来限定骨料的凸性；文献［16－17］在圆形骨料圆周上随机生成3个或4个点以形成多边形骨料基框架，然后在以多边形最长边为直径的圆内生成多边形的新顶点进行延凸；文献［18］用消除骨料凹面的方法来保持凸性.本文在随机生成的圆形或椭圆形骨料基础上，提出一种简单的凸多边形骨料生成方法，能在圆或椭圆形边界上随机生成凸多边形的顶点以保证其凸性，凸多边形顶点的位置由极角θ（j）确定，如图1所示.

图1 多边形骨料极坐标表示Fig.1 Polygonal aggregate in polar coordinates

凸多边形任意一边的对角φ（j）＝θ（j＋1）－θ（j）.因它有h条边，故φ（j）的平均值为2π／h.假定φ（j）根据平均值进行变化，则可表示为：

为满足所获得的所有φ（j）之和等于2π这一条件，需对φ（j）作如下修正：

修正后，各顶点的极角为：

式中：rand为0～1 间按均匀分布产生的随机数；δ表示φ（j）的波动程度，取值为0～1间的任意数，这里取0.5；α为骨料方位的相角，取值为0～2π间均匀分布的随机数；j为凸多边形的顶点编号.

顺时针连接各个顶点，即生成圆基或椭圆基凸多边形，调用MATLAB 中polyarea函数可计算该凸多边形的面积.骨料库的生成步骤概述如下：

（1）计算各级配段骨料生成的控制参数.根据实际需要选取全级配混凝土二维试件尺寸作为骨料投放区域的大小，将混凝土三维试件骨料体积分数Pk和骨料粒径最大值Dmax代入式（1），分别得到粒径小于某级配段粒径上限（Du）和下限（Dd）的骨料在二维试件截面上出现的概率Pc（D＜Du）和Pc（D＜Dd），两者差值即为该级配段骨料投放面积占混凝土二维试件截面的面积率，将其作为各级配段骨料生成的控制参数.

（2）分级配段生成骨料库.在某级配段内随机生成骨料并计算其面积率.若累积面积率小于该级配段骨料生成控制参数，则继续该级配段骨料的生成；若累积面积率大于该级配段控制参数，则放弃最后这颗骨料的生成，重新在该级配段内生成一颗更小的骨料，直到当前级配段生成的骨料累积面积率达到该级配段骨料生成控制参数的0.99～1.01倍，即可终止骨料生成.

1.3 基于数值图像处理的随机骨料投放

就全级配混凝土而言，由于骨料颗粒较多，未定位的骨料易与前面已定位的骨料发生冲突.为了提高骨料的投放效率，用数值图像处理技术将骨料在区域内服从均匀分布转化为服从条件均匀分布，即在投放区域中若某些区域被前面投放的骨料所占据，则给这些被占区域一个0的选择概率，而给那些自由区域一个较大的选择概率.

对一幅图像采样时，若每行采样数为M，每列采样数为N，则图像大小为M×N 个像素.如果每个像素点的颜色用红、绿、蓝3 种颜色的强度值表示，取值为0～255，则MATLAB 读取一幅M×N个像素点的RGB 图像时可用三维数组X 表示.若图像中某像素点P 的颜色强度值为（255，0，0），则该像素点为红色；若某像素点的颜色强度值为（255，255，255），则该像素点为白色.

基于上述成像原理，以一幅与骨料投放区大小一致的白色图像作为骨料投放背景区，以图形窗口左上角为坐标原点，投放的骨料用红色表示.初始投放时，图像中各像素点的颜色强度值均为（255，255，255）.开始投放时，从数组X 的第3面X（∶，∶，3）中找出数据为255的元素所在的行号和列号（即取样点的坐标），将它们分别储存在矩阵H 和L 中.从矩阵H 中随机选取1个元素s作为待投放骨料形心的横坐标，在矩阵L 中选取1个与矩阵H 中元素s相对应的元素t 作为待投放骨料形心的纵坐标，这样便随机选择到一个未被骨料占据的点作为骨料投放的形心位置.为加强骨料位置的随机性，按下式修正骨料形心位置：

投放骨料成功后，调用MATLAB 中的patch函数将其所占据的区域填充为红色，像素点颜色强度值变为（255，0，0），如图2所示.后继骨料的投放只需寻找颜色强度值为（255，255，255）的像素点作为随机样本进行遴选，重复上面的工作即可.

图2 数值图像处理骨料投放示意图Fig.2 Placing process represented by numerical image processing technology

骨料投放时，骨料影响范围系数的大小直接影响骨料投放的成功率：骨料影响范围系数越大，骨料间距就越大，骨料投放率会变得很低，反之，骨料投放间距越小，骨料投放就越容易.Schlangen等［19］给出的骨料影响范围系数为1.1；Wittmann等［20］认为骨料影响范围系数与骨料含量有关，骨料含量越高，骨料影响范围系数越小；Wang等［21］通过观察混凝土中大量骨料交界面，发现界面厚度与骨料粒径有关.本文将骨料影响范围系数取为1.05.

圆形骨料相互嵌入的判断条件相对简单，判断两圆心之间的距离是否大于它们的半径之和乘以骨料影响范围系数即可.

关于椭圆形骨料相互嵌入的判断条件，文献［22］利用椭圆的外接多边形逼近椭圆，把椭圆之间的嵌入判断转化为多边形之间的嵌入判断，此方法的缺点是：若多边形边数太少，将影响逼近精度而出现误判；若边数太多，又会影响判断效率.文献［23］利用两椭圆的代数条件来判断是否相互嵌入，该方法虽然在效率上得到了优化，但是所建立的广义特征多项式为一元三次方程，求解方程时如果精度设定不当，也会出现误判.本文提出一种新的判断两椭圆间位置关系的算法，其基本原理是根据椭圆－椭圆外切的代数条件，利用仿射变换和逆变换方法，推导出椭圆－椭圆的不适合边界（NFB）轨迹方程［24］来判断二者的关系.两椭圆的平面位置关系如图3所示，椭圆E（n）和E（i）的不适合边界轨迹方程可表示为：

图3 两椭圆的平面位置关系Fig.3 Planar position relationship of the two ellipses

函数F＞0表示椭圆E（i）的中心点在NFB外，两椭圆不相互嵌入；函数F＝0表示椭圆E（i）的中心点在NFB上，两椭圆不相互嵌入；F 为虚数表示椭圆E（i）的中心点在NFB内，两椭圆相互嵌入.

凸多边形骨料相互嵌入情况有两种，如图4所示.目前的判断方法主要有3种，分别是夹角之和测试法［16，18，25］、面积判别准则［10，17］和基于背景网格的侵入判断［26］.本文直接调用MATLAB中的inpolygon函数来判断凸多边形骨料是否相互嵌入.

图4 凸多边形的嵌入形式Fig.4 Intersection forms of the convex polygons

2 有限元网格的自动剖分

2.1 黏结界面和砂浆基质的生成

为了在细观层面上进行全级配混凝土的数值模拟，需要在生成的细观随机骨料模型基础上生成黏结界面和砂浆基质.假定黏结界面等厚度包裹在骨料外面，对圆形、圆基凸多边形、椭圆形和椭圆基凸多边形骨料以形心为缩放中心，以骨料影响范围系数为缩放因子进行缩放，可生成随机骨料的黏结界面.砂浆基质是除骨料及黏结界面以外的投放区域.

2.2 细观有限元网格的剖分

对细观几何模型进行有限元网格剖分常用方法有：映射网格法［17］、改进的推进波前法（AFT）［23］、Delaunay三角剖分法［27］和改进的Delaunay三角剖分法（渐变网格剖分法）［28］等.本文采用的COMSOL软件可以自动进行结构化或非结构化网格剖分，并支持移动网格、虚拟几何剖分网格、稳态及瞬态自适应网格.COMSOL软件中非结构化三角形网格剖分内置了2种算法，即Delaunay算法和前沿算法.Delaunay三角剖分使得所有三角形互不重叠，完整地覆盖了整个问题域，所有的节点都成为Delaunay三角形的顶点，且能够尽可能地避免病态三角形的出现，生成比较规则的三角形.应用COMSOL软件及MATLAB软件，可以方便地完成混凝土二维细观模型的自动生成.

3 数值模型实例

3.1 骨料投放实例

为了验证随机骨料模型算法的有效性，以截面尺寸为450mm×450mm 的混凝土试件为例，混凝土中骨料级配设定为：m（ELA（特大石））∶m（LA（大石））∶m（MA（中石））∶m（SA（小石））＝3∶3∶2∶2，将每个级配段粒径等分为两级，根据式（1）确定试件内截面各级配段骨料应投放的面积占试件截面面积的百分比，体积分数为72%的目标骨料二维随机骨料投放统计结果见表1，所建立的随机骨料模型如图5所示.用配置为AMD Athlon（tm）ⅡX2 215Processor 2.70GHz，内存为2.0GB的计算机进行计算，完成圆形、椭圆形、圆基凸多边形和椭圆基凸多边形模型分别耗时19.890，20.670，20.417，22.026s.

表1 二维随机骨料投放统计结果Table 1 Statistic results of placing aggregates in two－dimensional models

图5 4种随机骨料几何模型Fig.5 Geometric models of four kind of random aggregates

二维试件骨料实际投放量在70%左右，完全可以满足全级配混凝土试件数值仿真模拟的要求.这里的模型是按各级配骨料粒径比例和相应的骨料含量，在每一粒径范围内等机率选取粒径大小而生成连续粒径的随机骨料模型.然而本文提出的随机骨料投放方法还可以生成任意代表粒径的随机骨料模型.如取各级配粒径范围的平均值作为代表粒径，按照上述各级配粒径骨料比例和目标投放量，可分别得到特大石4颗、大石15颗、中石44颗、小石127颗.所建立的代表粒径随机骨料模型见图6.

图6 4种代表粒径骨料几何模型Fig.6 Geometric models of four kind of representative aggregate particles

3.2 网格划分实例

选取试验编号为QS1的某大坝全级配混凝土数据建立劈拉试验二维细观模型，试件尺寸为450mm×450mm×450mm，混凝土骨料用量2 212 kg／m3，密度2 717.5kg／m3.通过式（1）将三维骨料级配转换为二维级配，得到试件截面内骨料粒径115～150 mm，80～115 mm，60～80 mm，40～60mm，30～40mm，20～30mm，15～20mm，10～15mm 的面积率分别为0.070 5，0.112 5，0.082 2，0.099 7，0.059 8，0.071 1，0.066 3，0.026 4，骨 料总体积分数为58.85%.

进行细观有限元网格剖分时黏结界面最大单元尺寸取界面厚度的1／3.骨料区域单元尺寸相对较大，砂浆区域单元尺寸较小，划分网格后的二维有限元模型见图7.由图7 可以看出，这些网格质量较好，能满足骨料和黏结界面、黏结界面和砂浆基质交界处网格一致性的要求.

图7 二维有限元模型中的4种随机骨料网格Fig.7 Meshes of four kind of random aggregates in two－dimensional finite element models

4 结论

全级配混凝土二维细观模型自动生成算法，在细观层面上为全级配混凝土数值模拟提供了一种有效的建模手段.在圆或椭圆上随机取点生成凸多边形保证了骨料的凸性，省去了繁琐的延凸编程.投放骨料时采用数值图像处理技术随机选取形心，分级配逐个随机投放，保证了各级配骨料体积分数能满足投放要求，提高了骨料的投放效率.将细观几何模型导入COMSOL 软件进行细观有限元网格剖分，可以对不同区域剖分不同尺寸的网格.在进行材料属性赋值时，可以按材料性质成集后直接赋值.数值实例表明，所建立的全级配混凝土细观模型，对4种形状骨料的投放率都能达到70%左右，生成的有限元网格能满足进一步计算的需要.

［1］SCHLANGEN E，GARBOCZI E J.Fracture simulations of concrete using lattice models：Computational aspects［J］.Engineering Fracture Mechanics，1997，57（2／3）：319－332.

［2］MOHAMED A R，HANSEN W.Micromechanical modeling of concrete response under static loading－part 1：Model development and validation［J］.ACI Materials Journal，1999，96（2）：196－203.

［3］BAZANT Z P，TABBARA M R.Random particle model for fracture of aggregate or fiber composites［J］.Journal of Engineering Mechanics，1990，116（8）：1686－1705.

［4］刘光廷，王宗敏.用随机骨料模型数值模拟混凝土材料的断裂［J］.清华大学学报：自然科学版，1996，36（1）：84－89.LIU Guangting，WANG Zongmin.Numerical simulation study of fracture of concrete materials using random aggregate model［J］.Journal of Tsinghua University：Science and Technology，1996，36（1）：84－89.（in Chinese）

［5］唐春安，朱万成.混凝土损伤与断裂——数值试验［M］.北京：科学出版社，2003：24－43.TANG Chun'an，ZHU Wancheng.Damage and fracture for concrete—Numerical tests［M］.Beijing：Science Press，2003：24－43.（in Chinese）

［6］ZUBELEWICZ A，BAZANT Z P.Interface element modeling of fracture in aggregate composites［J］.Journal of Engineering Mechanics，1987，113（11），1619－1630.

［7］邢纪波.梁－颗粒模型导论［M］.北京：地震出版社，1999：15－38.XING Jibo.Introduction to the beam－particle model［M］.Beijing：Seismological Press，1999：15－38.（in Chinese）

［8］展辰辉.混凝土及其组成材料动力特性研究［D］.南京：河海大学，2005.ZHAN Chenhui.The study of dynamic properties on concrete and its component materials［D］.Nanjing：Hohai University，2005.（in Chinese）

［9］徐波.基于材料细观结构的混凝土数值模拟与性能分析［D］.杭州：浙江大学，2008.XU Bo.Numerical simulation and performance studies of concrete based on material meso－structures［D］.Hangzhou：Zhejiang University，2008.（in Chinese）

［10］高政国，刘光廷.二维混凝土随机骨料模型研究［J］.清华大学学报：自然科学版，2003，43（5）：710－714.GAO Zhengguo，LIU Guangting.Two－dimensional random aggregate structure for concrete［J］.Journal of Tsinghua University：Science and Technology，2003，43（5）：710－714.（in Chinese）

［11］梁建，娄宗科，韩建宏.基于AUTOCAD的混凝土骨料建模分析［J］.水利学报，2011，42（11）：1379－1383.LIANG Jian，LOU Zongke，HAN Jianhong.Modelling analysis of the concrete aggregate based on AUTOCAD［J］.Journal of Hydraulic Engineering，2011，42（11）：1379－1383.（in Chinese）

［12］程伟峰.混凝土架构模型的数值模拟研究［D］.大连：大连理工大学，2008.CHENG Weifeng.Research on concrete frame and structure by numerical simulation［D］.Dalian：Dalian University of Technology，2008（in Chinese）

［13］高利甲，刘锡军，王玉梅.基于matlab－混凝土二维细观结构数值模拟骨料随机投放［J］.湖南工程学院学报，2011，21（1）：81－84.GAO Lijia，LIU Xijun，WANG Yumei.Numerical simulation for random embarking of aggregate with two－dimensional meso－structure of concrete based on Matlab［J］.Journal of Hunan Institute of Engineering，2011，21（1）：81－84.（in Chinese）

［14］伍君勇.混凝土细观结构的自动生成［D］.大连：大连理工大学，2006.WU Junyong.Automatic generation of mesostructure of concrete［D］.Dalian：Dalian University of Technology，2006.（in Chinese）

［15］WALRAVEN J C，REINHARDT H W.Theory and experiments on the mechanical behavior of cracks in plain and reinforced concrete subjected to shear loading［J］.Heron，1991，26（1A）：26－35.

［16］孙立国，杜成斌，戴春霞.大体积混凝土随机骨料数值模拟［J］.河海大学学报：自然科学版，2005，33（3）：291－295.SUN Liguo，DU Chengbin，DAI Chunxia.Numerical simulation of random aggregate model for mass concrete［J］.Journal of Hohai University：Natural Science，2005，33（3）：291－295.（in Chinese）

［17］马怀发，芈书贞，陈厚群.一种混凝土随机凸多边形骨料模型生成方法［J］.中国水利水电科学研究院学报，2006，4（3）：196－201.MA Huaifa，MI Shuzhen，CHEN Houqun.A generating approach of random convex polygon aggregate model［J］.Journal of China Institute of Water Resources and Hydropower Research，2006，4（3）：196－201.（in Chinese）

［18］张剑，金南国，金贤玉，等.混凝土多边形骨料分布的数值模拟方法［J］.浙江大学学报：工学版，2004，38（5）：581－585.ZHANG Jian，JIN Nanguo，JIN Xianyu，et al.Numerical simulation method for polygonal aggregate distribution in concrete［J］.Journal of Zhejiang University：Engineering Science，2004，38（5）：581－585.（in Chinese）

［19］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（9）：534－542.

［20］WITTMANN F H，ROELFSTRA P E，SADOUKI H.Simulation and analysis of composite structures［J］.Materials Science and Engineering，1984，68（2）：239－248.

［21］WANG Z M，KWAN A K H，CHAN H C.Mesoscopic study of concrete I：Generation of random aggregate structure and finite element mesh［J］.Computers and Structures，1999，70（5）：533－544.

［22］李友云，焦俊芝.具有大量椭圆颗粒／孔洞随机分布区域的计算机模拟及其改进三角形自动网格生成算法［J］.计算力学学报，2004，21（5）：540－545.LI Youyun，CUI Junzhi.Computer simulation method for the domain with large numbers of random ellipse grains／cavities and the improving automatic triangle mesh generation algorithm［J］.Chinese Journal of Computational Mechanics，2004，21（5）：540－545.（in Chinese）

［23］高巧红，关振群，顾元宪，等.混凝土骨料有限元模型自动生成方法［J］.大连理工大学学报，2006，46（5）：641－646.GAO Qiaohong，GUAN Zhenqun，GU Yuanxian，et al.Automatic generation of finite element model for concrete aggregate［J］.Journal of Dalian University of Technology，2006，46（5）：641－646.（in Chinese）

［24］曾龙，陈羽，腾宏飞，等.椭圆－椭圆不适合边界的解析方程及其干涉判断［J］.中国图像图形学报，2007，12（5）：917－921.ZENG Long，CHEN Yu，TENG Hongfei，et al.No fit boundary analytic equation and interference detection for ellipse－ellipse［J］.Journal of Image and Graphics，2007，12（5）：917－921.（in Chinese）

［25］杜成斌，孙立国.任意形状混凝土骨料的数值模拟及其应用［J］.水利学报，2006，37（6）：662－667.DU Chengbin，SUN Liguo.Numerical simulation of concrete aggregates with arbitrary shapes and its application［J］.Journal of Hydraulic Engineering，2006，37（6）：662－667.（in Chinese）

［26］秦川，郭长青，张楚汉.基于背景网格的混凝土细观力学预处理方法［J］.水利学报，2011，42（8）：941－948.QIN Chuan，GUO Changqing，ZHANG Chuhan.A pre－processing scheme based on background grid approach for mesoconcrete mechanics［J］.Journal of Hydraulic Engineering，2011，42（8）：941－948.（in Chinese）

［27］王宝庭，宋玉普，赵国藩.混凝土随机颗粒模型的网格自动剖分方法［J］.大连理工大学学报，1999，39（3）：445－450.WANG Baoting，SONG Yupu，ZHAO Guofan.Automatically dividing way of mesh for random particle model of concrete［J］.Journal of Dalian University of Technology，1999，39（3）：445－450.（in Chinese）

［28］夏晓舟.混凝土细观数值仿真及宏细观力学研究［D］.南京：河海大学，2007.XIA Xiaozhou.The meso numerical simulation and the macromeso mechanics research for concrete material［D］.Nanjing：Hohai University，2007.（in Chinese）