混凝土细观结构的数值模拟研究

宋来忠 廖大乾

(三峡大学 理学院， 湖北 宜昌 443002)

混凝土可视为由大量随机分布的增强颗粒—骨料以及沙子和水泥所组成的混合物—砂浆所组成的多相复合材料．起初，人们先假设增强颗粒形状为球形来进行计算机模拟．由于增强颗粒形状一般都不是球形，人们进一步假设颗粒形状为椭球[1-2]；还有学者将增强颗粒的形状假设为广义椭球体进行数值模拟[3]．为了使所模拟的增强颗粒的形状更加接近自然形状，一些学者[4-6]假设颗粒为多面体，研究了“大量多面体随机分布的空间区域”的计算机模拟．作者[7-8]使用自由曲线曲面变形技术，将椭球变形得到了更加接近自然增强颗粒形状的形体，去掉了一些学者在研究中的“凸性”限定条件，并且给出了这种颗粒形体的参数方程，得到了更加接近实际的混凝土材料的细观结构．但有一个问题：如此模拟的颗粒过于不规则，以至于难以确定它的内部与外部．由于星形不一定是凸的，但却可以判断其内部与外部；所以，本文在文献[7-8]的基础上，对颗粒的形状给出一定的限制，使其成为星形，直接给出空间点与这类星形颗粒位置关系的判别法则、空间点与这类星形颗粒的距离的近似计算与误差估计，使得颗粒的影响区域在试件创建中减小，颗粒在试件中的分布更加合理、密度加大，从而能创建出具有更高颗粒含量的模型试件．

1 基础知识

1.1 星形的定义

定义1如果一个简单闭区域Ω上至少存在一点0∈Ω，使得对任意P∈Ω都有直线段OP∈Ω，则称Ω关于点O是星形．

显然，所有的凸单连通的简单闭区域都是星形．如果用I(Ω)表示区域Ω的内部，E(Ω)表示区域Ω的外部，∂Ω表示区域Ω的边界，假如Ω关于点O是星形的，那么对空间中任意一点，过P的直线与Ω边界的∂Ω交点不超过2．并且易得如下判定定理．

定理1设Ω为关于点O是星形的，若空间任意点P0(x0，y0，z0)与O的连线交边界∂Ω于P1(x1，y1，z1)，P2(x2，y2，z2)两点，那么

(1)

1.2 星形参数化描述及条件

作者在文献[8]引入了一种新的伸缩因子函数，下面将要使用这个伸缩因子将椭球变形．

定义2设0

令

t=(u-u0)2+(v-v0)2(2)

作R2上的连续函数：

(3)

令

(4)

则B(u，v)被称为是R2上的基本伸缩因子．其中Ω1是峰值区域，而Ω2为变形区域．

椭球面的参数方程可表示为：

(5)

其中I=[0，2π]×[0，π]，T=(tij)，i，j=1，2，3，且t11=cosφcosψ-sinφcosθsinψ，t12=cosφsinψ+sinψcosθcosφ，t13=sinφsinθ，t21=-sinφcosψ-cosφcosθsinψ，t22=-sinφsinψ+cosφcosθcosψ，t23=cosφsinθ，a31=sinθsinψ，t32=-sinθcosψ，a33=cosθ.而θ，φ，ψ为欧拉角，且θ∈[0，π)，ψ∈[0，2π)，φ∈[0，2π)．P(u，v)为定义在参数(u，v)平面R2中区域I上的光滑曲面．

设变形中心为O′的取峰值区域Ω1和支区域Ω2(Ω1⊂Ω2⊂Ω)．取伸缩系数α1，α2，α3∈R，伸缩系数矩阵记为D=diag(α1，α2，α3)，丰满指数k取为正整数．那么，变形之后的参数曲面Pd(u，v)与变形之前的椭球面P(u，v)就有如下关系：

Pd(u，v)=(DB(u，v)k+E)(P(u，v)-O′)+O′，

(u，v)∈I(6)

显然，边界曲面可由(7)表示的区域，不一定是星形区域．但实际上，有

定理2如果取O′=O，只要αi≥-1，i=1，2，3．那么，边界为(7)的Ω关于点O为星形区域．

满足定理2条件的方程(7)的空间图形是星形(如图1所示)．

(u，v)∈I(7)

其中，在第i个变形的区域中，峰值区域为Ω1i，支区域为Ω2i，它们分别为：

由于伸缩系数矩阵Di=diag(α1i，α2i，α3i)的所有元素都是大于等于-1的，故对于空间点P′(x′，y′，z′)=P′(x′(u′，v′)，y′(u′，v′)，z′(u′，v′))，可考虑函数：

g(u′，v′)=

那么，有

定理3设P′是空间中任意一点，那么

(10)

1.3 点到星形体的距离计算的搜索算法及误差

定义3空间中任意一点P′(x′，y′，z′)到∂Ω的距离，可以定义为

(11)

(u，v)∈I(12)

其参数区域是闭区域I=[0，2π]×[0，π]，所以可以考虑如下快速的搜索法算法：

第1步：给定自然数n，取ui+1=2iπ/n，vj+1=jπ/n；i=0，…，n-1，j=0，…，n-1

第2步：由式(8)计算星形表面上对应的n2个点：Pd(x(ui，vj)，y(ui，vj)，z(ui，vj))

第3步：将上一步得到的n2个点的坐标，代入目标函数，计算n2个函数值：f(x(ui，vj)，y(ui，vj)，z(ui，vj))

第4步：取min{f(x(ui，vj)，y(ui，vj)，z(ui，vj))}

作为所求距离d(P，∂Ω)的近似值．

定义4空间中任意一点P′到Ω的距离，可以定义为

(13)

如图2所示．误差满足

可得距离的误差范围：

(14)

图1 椭球变为星形 图2 点到骨料的距离误差图

2 基本算法

以试件中需生成的星形总个数的最大值N、或者试件中需生成的星形的总体积数S(星形体积的计算见文献[8])、或者模拟时搜索中心连续被拒绝总次数K，来控制终止．用k记连续选取中心点O而被拒绝的次数；用n记已生成的星形个数；用s记已生成的星形的总体积数．主要步骤如下．

第1步：设定产生随机星形的各种参数及其分布区域、试件中星形的总个数最大值N、星形的体积总数S以及搜索中心连续被拒绝总次数K，给初值n=1，k=0，s=0;

第2步：根据颗粒半径的范围，随机产生一个数a；

第3步：在模拟区域内随机选取一点O(x0，y0，z0)，判断它是否是在所有已生成的星形颗粒的外部；若否，则令k=k+1，转到第4步；若是，由计算距离的搜索算法计算O与所有的已生成星形颗粒之间的距离，并判断最小的距离d是否≥a(注意控制误差)；若否，则令k=k+1，转到第4步；若是，则取椭球的中心为O(x0，y0，z0)，长半轴长为a，生成一个椭球，再将椭球随机变形，使之成为星形颗粒，并计算星形的体积s0，令s=s+s0，n=n+1，到到第4步；

第4步：如果满足n>N，或k>K，或s>S，结束程序；否则，转到第2步．

3 混凝土数值模拟实例

由于界面和骨料的尺度差异太大，受计算工具的限制，难以用有限元方法对其进行力学分析计算．所以，本文采用其他一些学者的做法[4-6，9]，假设混凝土是具有星形增强颗粒以及砂浆所组成的两相复合材料，伸缩系数控制在-0.16～-0.1之间．然后，根据文献[10]所给的级配曲线模拟生成混凝土试件，如图3、图4所示．

图3 骨料含量为55%的 图4 骨料含量为60%的二级配随机模拟 三级配随机模拟

图3是假设增强骨料形状为星形，模拟的骨料含量为55%的二级配混凝土试件，模拟区域为[0，15]×[0，15]×[0，15](单位：cm)，颗粒粒径为0.5∶2，2∶4两种级别．两种骨料的体积分别是：839.165 6 cm3和1 017.202 2 cm3，所有骨料的总体积为1 856.367 8 cm3占到试件体积的55%．两种骨料的配合比大约为4.5∶5.5．每种骨料颗粒数为3 527和130，骨料颗粒总数为3 657．

图4是假设增强骨料形状为星形，模拟的骨料含量为60.04%的三级配混凝土试件，模拟区域[0，30]×[0，30]×[0，30](单位：cm)，颗粒粒径分别为0.5∶2,2∶4和4～8三种．每种骨料的体积分别为：4 050.08 cm3，4 864.87 cm3和7 295.18 cm3，骨料总体积为16 210.15 cm3占试件总体积的60.04%．3种骨料的配合比约为2.5∶3∶4.5．每种骨料个数为5 128和229，33，骨料颗粒总数为5 380．显然，比文献[8]的算法生成试件的骨料含量高．

值得注意的是：误差估计式(14)，必须充分大才能成立，但n若取得太大，会影响计算速度．所以，在作模拟时，一定要注意检查，以便选取适当的n值．

4 有限元网格模型实例

本文选择背景网格剖分，网格全部采用Abaqus软件中的C3D8R类型单元划分．将正方体的8个顶点作为背景网格的基本点来判断，如果8个顶点有5个都属于骨料单元，则这8个点组成的网格单元属于骨料单元，否则，这8个点组成的网格单元属于砂浆单元．

以图3骨料含量为55%的几何模型为例，对15 cm×15 cm×15 cm区域的用背景网格划分，将其边界其以3 mm为步长，分成50等分，共125 000个网格单元．得到的骨料单元51 284个，砂浆单元73 716，其中骨料的网格单元的含量为41%．如图5～7所示．要说明的是：选用背景网格，当几何模型离散成有限元单元之后，一些较小的骨料，将会被砂浆“融掉”所以会降低骨料的含量．

图5 砂浆背景网格图 图6 骨料背景网格图 图7 砂浆和骨料混合的背景网格图

5 力学分析模拟试算

下面采用Abaqus软件中混凝土塑性损伤模型的本构关系[9]，对上述有限元模型，进行力学分析试算．材料参数按文献[10]选取，如表1所示．

表1 混凝土材料力学参数

对上述有限元模型试件底部添加固定约束，如图8所示．为了方便处理，在试件的上端添加了一个参考点RP1，然后对顶部的上表面第一层与该点进行全约束，保证顶部均匀受力，同时通过该点来体现整个试件的应力与应变的关系[11]，如图9所示．

图8 试件底部约束图 图9 试件顶部约束图

然后，在试件的顶部施加的是方向向下的2 mm的位移荷载，得到试件受力过程的应力分布云图结果如图10所示．加载结果的破坏过程如图10所示．图10(a)表示在没有加外力的作用下，混凝土试件的初始状态．图10(b)表示的是混凝土试件局部区域开始出现损伤单元．图10(c)表示混凝土试件的顶面所受的应力表现的不如正面的应力那么明显，其中正面的应力比较集中，逐渐呈现出“X”形状的趋势．从图10(d)可以看出，随着荷载的施加，试件损伤单元沿试件45度方向慢慢扩展，最终形成一条宏观裂缝．图10(e)表示的是在极限应力作用下，两个受力面的三维视图，其中两个应力面均显示了混凝土试件的损伤单元都是沿着45°方向扩展．这与文献[12]报道的真实的混凝土试件在受压状态下的损伤结构基本一致．

图10 三维混凝土试件单轴压缩破坏损伤过程图

对比图11和Abaqus手册[9]中的标准混凝土单轴压缩损伤图像图，可以看出两者变化趋势基本一致．

图11 应力应变关系图

6 结 语

本文在文献[7-8]的基础上，对复合材料的增强颗粒的形状给出一定的限制，使其成为星形，创建出比以往文献中报道的颗粒含量更高的试件．为进一步从细观力学的角度，分析具有增强颗粒的复合材料的性能，提供了创建更加接近真实试件的几何模型方法．从理论上来说，本文的方法也是文献[13]的三维推广．