CN110706352A - 基于多边形随机骨料的混凝土三相细观模型构建及内氯离子侵蚀数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于多边形随机骨料的混凝土三相细观模型构建方法，首先设置混凝土模型的变量；然后根据混凝土模型中骨料粒径范围生成凸多边形随机骨料；并将生成的骨料逐一投放至模型区域范围直至全部结束，并保存所有凸多边形骨料顶点点线数据至元胞数组中；根据界面过渡区厚度将所有骨料边平移、连接生成界面过渡区；最后保存所有骨料及界面过渡区的点线数据。本发明在上述混凝土三相细观数值模型基础上，真实模拟了其内氯离子侵蚀过程，计算过程简单明了，可大幅提升数值模拟计算效率。

Description

基于多边形随机骨料的混凝土三相细观模型构建及内氯离子 侵蚀数值模拟方法

技术领域

本发明涉及混凝土细观数值模拟技术领域，特别是一种基于多边形随机骨料的混凝土三相细观模型构建方法及其内氯离子侵蚀数值模拟方法。

背景技术

钢筋混凝土自19世纪中叶出现，在国家基本建设中有着广泛地应用，但在海洋环境条件下，钢筋混凝土耐久性破坏一直是威胁结构使用寿命的一个重要且亟待解决的问题。世界各国在建筑行业的高速发展中都有因结构耐久性破坏而产生了巨大经济损失，很多建筑结构的后续维护费用远超建造费用，并且随着年份的增加，维护成本还会进一步提升，特别是沿海地区钢筋混凝土结构耐久性破坏现象比较严重，海洋环境下钢筋混凝土结构大部分都属于露天工程，所处环境复杂，所受影响因素较多。

其中，混凝土中的钢筋锈蚀、混凝土碳化、混凝土冻融循环破坏被认为是造成结构耐久性破坏的主要因素，而在三者之中，钢筋锈蚀被认为是造成结构耐久性破坏的最主要的因素。因此，研究在海洋环境下氯离子在混凝土内部的侵蚀传输规律，对开展对应手段延缓钢筋的锈蚀以及提高混凝土结构的耐久性有着重要的意义。

早期的对混凝土中氯离子侵蚀相关研究成果大多是通过物理试验得到的，成果主要集中于描述氯离子侵蚀规律以及修正氯离子扩散模型上，包括不同因素对菲克第二定律中氯离子浸透系数的影响。但进行物理试验往往需要耗费较大的财力、物力、人力，且难以得到长期观测的结果。随着计算机技术以及应用数学理论的发展，数值计算能够精确的得到微分方程解析解。由于采用菲克第二定律描述氯离子侵蚀混凝土的过程较为有效，越来越多的学者开始考虑建立数值模型研究混凝土中氯离子的侵蚀特性。混凝土是一种复杂的多相复合性物质，现阶段数值模拟研究中普遍将混凝土的细观结构考虑为由水泥浆体、骨料以及界面过渡区共同组成的三相非均质材料。其随机骨料的生成及投放技术是混凝土强度、耐久性计算研究的重要部分。

现有的混凝土数值模拟研究成果中，多数学者将混凝土中的骨料简化为圆形、球形或其他规则多边形形状，然而实际工程中大多采用人工碎石作为混凝土骨料，其形状更接近凸多边形或凸多面体。为保证数值模拟的准确性，建立数值模型应尽量符合真实的混凝土中骨料的形状、尺寸以及分布。申请号CN2017101502969，专利名称为一种混凝土骨料单元的二维随机生成方法中未考虑水泥浆体和骨料之间的界面过渡区。

发明内容

有鉴于此，本发明的目的在于提供一种基于多边形随机骨料的混凝土三相细观模型构建方法，该方法将从细观角度出发，根据混凝土中骨料的实际形状开发新型骨料生成及投放算法，建立混凝土随机凸多边形骨料细观模型，用于对混凝土中氯离子的侵蚀规律展开研究。

为达到上述目的，本发明提供如下技术方案：

本发明提供的基于多边形随机骨料的混凝土三相细观模型构建方法，包括以下步骤：

S1：设置混凝土数值模型的变量；

S2：根据混凝土数值模型的变量生成凸多边形随机骨料；

S3：按照以下公式计算已生成骨料的累积体积分数V：V＝V+V_i；其中，V_i表示第i个凸多边形随机骨料的体积分数；

S4：判断已生成的骨料的累积体积分数V是否满足以下关系：V≥r*X*Y；如果否，则返回继续生成凸多边形随机骨料；如果是，则将生成的所有凸多边形骨料粒径d_i及凸多边形骨料顶点的点线数据保存；其中，r表示生成凸多边形骨料的目标体积分数，X表示拟构建混凝土数值模型的长度，Y表示拟构建混凝土数值模型的宽度；

S5：将已生成好的凸多边形骨料按照骨料粒径d_i从大到小的顺序进行降序排列；

S6：随机生成第i个凸多边形骨料外接圆心坐标(x_i,y_i)；

S7：随机生成第i+1个凸多边形骨料外接圆心坐标(x_i+1,y_i+1)；

S8：根据凸多边形骨料外接圆心坐标(x_i,y_i)、(x_i+1,y_i+1)来判断第i个和第i+1个凸多边形骨料是否会发生重叠。若重叠，则重新随机生成新的凸多边形骨料外接圆心坐标；若不重叠，则按顺序保存满足要求的凸多边形骨料顶点坐标并将各顶点相连接，并将这第i个、第i+1个凸多边形骨料按照满足要求的外接圆心坐标(x_i,y_i)、(x_i+1,y_i+1)投放到混凝土数值模型的X*Y区域当中；

S9：将S5中生成好的凸多边形骨料按照S8的重叠判别准则逐一投放直至全部结束后，将所有凸多边形骨料顶点数据保存至元胞数组f_xy中；

S10：将全部已成功投放的凸多边形骨料各边按照界面过渡区的厚度由内向外平移，平移后的各边沿其直线方向相延伸连接，形成紧邻凸多边形骨料各边的等厚界面过渡区相；

S11：将全部生成好的界面过渡区顶点的点、线数据保存至元胞数组f_itz中；

S12：将满足S9和S10步骤的全部凸多边形骨料及其界面过渡区顶点的点、线数据保存至.dxf格式文件中。

进一步，所述步骤S2中采用多边形随机骨料颗粒生成方法，所述多边形随机骨料颗粒生成方法，具体按照以下步骤：

S21：设置混凝土数值模型的初始条件，包括混凝土模型的长宽尺寸X、Y，多边形随机骨料的目标体积分数r，骨料最小粒径D_min，骨料最大粒径D_max；

S22：生成第i个随机粒径为d_i的多边形骨料；

S23：判断多边形骨料粒径d_i是否满足D_min≤d_i≤D_max，若不满足，则重新执行步骤S22重新生成第i个随机粒径为d_i的多边形骨料，若满足，则进行下一步；

S24：生成直径为d_i的多边形骨料的外接圆，在外接圆上任取三点形成锐角三角形；

S25：储存锐角三角形的点线信息，计算锐角三角形的面积s；

S26：随机生成一个外接圆内的点，如果生成的点没有在多边形外部，则重新生成；如果生成的点在多边形外部，则将新生成的点与原多边形顶点联合形成新的多边形，新生成的多边形为凸多边形；

S27：对新生成多边形各顶点坐标进行重新排序、记录各点坐标，并计算多边形骨料的累加面积；

S28：设定多边形骨料边数n，重复循环至步骤S22开始执行，直到完成设置的边数n为止；

S29：将生成的所有多边形骨料的数据保存在对应矩阵中；

S210：将每个多边形的信息矩阵作为单个元素保存至元胞数组中，得到所有骨料信息的元胞数组XY。

进一步，所述步骤S26的凸多边形的判定，具体按照以下步骤进行：

判断新生成的点是否位于原三角形或多边形外部；以及新生成的点与原三角形或多边形组合形成的多边形为凸多边形；

如果否，则将返回继续判断；如果是，则进行下一步；

按照以下公式计算三角形面积：

式中，(x₁,y₁)、(x₂,y₂)以及(x₃,y₃)为三角形的三个顶点的坐标，且三个点是按逆时针排序；

计算新点P与相邻老点按逆时针顺序组合成的两个三角形的面积，如果所述三角形面积的结果都大于0时，则新生成的多边形为凸多边形。

进一步，所述步骤S9中的多边形骨料在进行投放的过程中还包括以下步骤：

S91：计算确定要投放的多边形骨料与已经投放的多边形骨料之间的圆心距；

S92：计算多边形骨料与已经投放的多边形骨料的外接圆半径之和；

S93：如果圆心距均大于外接圆半径之和，则直接投放，否则，再进行多边形重叠判别。

进一步，所述步骤S93中的多边形重叠判别还包括点不完全包含相交判断，具体如下：

新投放的凸多边形骨料的每个顶点位置都不在之前任何一个已经投放的骨料内部。

进一步，所述步骤S93中的多边形重叠判别还包括点不包含相交判断，所述点不包含相交判断是通过判断凸多边形的一条边与其他多边形各个边的线段相交来进行判断，具体过程如下：

假设A_i和A_i+1构成的线段A_i,i+1，B_j和B_j+1构成的线段B_j,j+1，A_i,i+1和B_j,j+1二者的交点坐标(x,y)计算式如下：

其中，线段A_i,i+1的顶点坐标分别为A_i(x_i,y_i)、A_i+1(x_i+1,y_i+1)；线段B_j,j+1的顶点坐标分别为B_j(x_j,y_j)、B_j+1(x_j+1,y_j+1)；

得到交点坐标(x,y)后，判断交点与线段的位置关系，若同时满足以下条件，则可以判定这两条线段相交：

式中：

min表示坐标元素(x_i,x_i+1,x_j,x_j+1,y_i,y_i+1,y_j,y_j+1)最小值；

max表示坐标元素(x_i,x_i+1,x_j,x_j+1,y_i,y_i+1,y_j,y_j+1)最大值；

根据上述两个公式得到交点坐标公式和判定这两条线段相交公式，对投放后外接圆圆心距小于半径之和、满足点不完全包含相交情况的多边形进行逐边判别，若多边形所有边与其他已投放多边形边交点坐标都不能同时满足两条线段相交公式，表示此多边形不与其他任何多边形相交，则按照此坐标位置投放。

进一步，所述步骤S10中多边形随机骨料颗粒还包括界面过渡区的生成过程，具体步骤如下：

将全部凸多边形随机骨料生成、投放完成后，根据各多边形骨料顶点按照逆时针的方向沿着连接的两条边外侧平移一段相等的距离，所述距离为界面过渡区的厚度；

随后将平移后的点连接，得到界面过渡区的顶点分布。

进一步，所述步骤S11是按照以下方式进行保存：

S111：保存所有界面骨料过渡区顶点的点线数据至元胞数组f_itz中；

S112：保存所有骨料、骨料界面过渡区顶点的点线数据至.dxf格式文件中。

本发明提供的基于多边形随机骨料的混凝土三相细观模型中氯离子侵蚀数值模拟方法，包括以下步骤：

S1：确定创建仿真环境；

S2：在仿真环境中建立混凝土试件几何模型，具体如下：导入已经创建的包含全部骨料及其界面过渡区信息的.dxf文件；所述.dxf文件是根据混凝土随机凸多边形骨料模型构建方法中得到的随机凸多边形骨料信息以及相应的界面过渡区信息的.dxf文件；建立混凝土细观数值模型的三相域，并构建成为联合体，完成几何模型的创建过程；

S3：设定数值模拟计算中需要定义的基础参数，所述基础参数包括初始氯离子浓度C₀、水泥浆体氯离子扩散系数D₀、龄期衰减系数m、界面过渡区扩散系数倍数a_itz；

设置解析函数扩散系数时变模型D_t与表面氯离子浓度时变模型C_s，并设定相应的变元，并确定其参数的取值；

S4：设定边界条件：

基础参数设定完成后，对混凝土数值模型的封闭面及扩散面进行设置；所述封闭面为不允许氯离子透过或浸入的表面，所述扩散面为允许氯离子透过或浸入的表面；

S5：设定混凝土细观数值模型中各组成相的传递属性；

S6：划分网格：采用自适应的三角细化网格划分；

S7：数值模拟计算及后处理：

网格划分完成后，选择瞬态计算，设置计算步长，进行计算；

计算完成后，生成整体氯离子浓度分布图。

进一步，所述骨料模型.dxf文件是通过以下步骤来形成的：

基于MATLAB生成骨料模型：录制骨料模型建立过程中的系列操作；

从系列操作中提取骨料模型建立的操作步骤；

对操作步骤进行取舍：保留需要的操作步骤的操作代码，并整理为可重复执行的骨料模型.dxf文件。

本发明的有益效果在于：

1、本发明通过建立随机凸多面体骨料的混凝土三相细观数值模型，真实模拟了混凝土中氯离子侵蚀扩散过程，根据混凝土的多相复合材料特性，详细描述了混凝土内部复杂结构。

2、在建立混凝土三相细观结构模型的基础上，可以详细研究粗骨料和界面过渡区对氯离子侵蚀的三个主要影响效应：稀释、曲折和界面过渡区效应，实现了混凝土骨料体积分数较大的情况下逾渗效应的研究，即：当骨料体积分数较大时，界面过渡区互相连通导致的渗透效应也会对氯离子浸透系数产生较大影响。

3、该模型能对混凝土中氯离子的长期侵蚀扩散特性进行研究，通过模拟真实情况，能对钢筋、荷载、冻融及其他多因素影响下的混凝土中氯离子侵蚀扩散规律进行研究，可挖掘完善氯离子的多重传输机理，建立可用于综合表征混凝土内氯离子侵蚀扩散的数学模型。

4、通过生成随机凸多边形骨料及界面过渡区，以及骨料投放算法，使得骨料形状和分布完全随机，其形状与内部构造都更加接近现实中的混凝土材料。

5、该算法简单明了，计算效率高，生成的效果符合现实混凝土骨料界面形状，在程序的随机数生成过程中，加入了时间作为随机种子参与随机数的生成，达到了真正随机生成的目的，大幅提升计算效率。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究，本发明的优势对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为混凝土随机凸多边形骨料生成流程图。

图2a为点不完全包含相交的判别。

图2b为点不包含相交的判别。

图3为混凝土随机凸多边形骨料模型程序的总体流程。

图4为通过图3程序生成的混凝土三相细观数值模型图。

图5为氯离子扩散边界设定示意图。

图6a为水泥浆体传递属性设定示意图。

图6b为界面过渡区传递属性设定示意图。

图6c为骨料传递属性设定示意图。

图7为模型网格划分示意图。

图8为混凝土细观数值模型中氯离子浓度分布图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好的理解本发明并能予以实施，但所举实施例不作为对本发明的限定。

实施例1

如图1所示，图1为混凝土随机凸多边形随机骨料生成流程图，本实施例提供的多边形随机骨料颗粒生成方法，包括以下步骤：

1、定义初始条件，包括混凝土试件的长宽尺寸X、Y，骨料的目标体积分数r，骨料最小粒径D_min，骨料最大粒径D_max，本实施例的混凝土的试件尺寸多为100mm，因此采用X、Y为100mm的混凝土试件进行数值计算，D_max、D_min、r可根据具体试验条件确定，定义储存骨料点线矩阵信息的初始数组及元胞数组，方便后续的信息存储与调用；

2、开始生成第i个随机凸多边形骨料，结合随机方法按照如下公式生成粒径为d_i的圆；

其中，P_(di)为骨料粒径小于d_i的累积频率；

D_min为骨料的最小粒径；

D_max为骨料的最大粒径；

unifrnd(0,1)为区间范围为(0,1)之间的任意随机数；

3、判断骨料粒径d_i是否满足D_min≤d_i≤D_max，若不满足，则重新执行步骤2，若满足，则进行下一步。

4、以直径为d_i的圆为基础外接圆，采用极坐标方式，利用unifrnd函数在圆上任取三点形成三角形，为了保证骨料实际粒径为d_i，外接圆需要保证为最小外接圆，需要保证此三角形为锐角三角形，如果为钝角三角形，则其实际外接圆直径为最长边，因此需要添加判定，本实施例选用MATLAB内建函数inpolygon进行判定，当圆心位于三角形内部时，三角形为锐角三角形，进行下一步。否则为钝角三角形，将重新生成三角形。

5、储存三角形点线信息，计算三角形面积s。

6、随机生成一个圆内的点，如果生成的点没有在多边形外部，则重新生成，如果在，则将新生成的点与原多边形合成新多边形，新生成的多边形作为凸多边形，即：同时进行以下判定：

a)保证该顶点位于原三角形或多边形外部。

此方法的判断方式同步骤4，采用内建函数inpolygon进行判定。如果判定结果位于原三角形或多边形内部，则重新生成新顶点；如果位于外部，则基于极坐标系找到角度差值分别为正负最小的相邻两旧顶点，将新顶点插入其中，与原多边形组合形成新多边形，并进行下一步判定。

b)保证该顶点与原三角形或多边形组合形成的多边形为凸多边形。

上述条件同时满足，则进行下一步，否则将重新执行本步骤。

为保证多边形为凸多边形，需要对新点与老点进行关系判别，为此，引入线性代数中的三角形面积公式：

式中：(x₁,y₁)、(x₂,y₂)以及(x₃,y₃)为三角形的三个顶点的坐标，并且这三个点是按逆时针排序的，当新点P与相邻老点按逆时针顺序组合成的两个三角形PA_iA_i+1、PA_i-2A_i-1面积按上述的三角形面积公式计算时，所得到的结果都大于0时，则新生成的多边形为凸多边形。

7、满足上述判别条件后，需要对新生成多边形各顶点坐标进行重新排序、记录各点坐标，并计算累加骨料面积，即：将新三角形面积Δs增至骨料面积s当中；

8、设定凸多边形骨料边数n，一般碎石骨料的边数在4～12，为减少骨料棱角，使骨料随机特征显著，取边数为12，实际应用可根据条件增加或减少。因此，当n<12时，i＝i+1，循环至步骤2开始执行。当n＝12时，停止此层循环，执行下一步骤。

9、生成骨料的过程中，需要随时保证点的连接顺序的合理性，这样多边形的信息才能完整地保存下来。采用角度分割、差值排序的方法，在每生成一个新的点后，对所有点线信息重新排序、连接，这样不仅方便进行多边形凸性的判别，同时可以将整个多边形的数据完整的保存在对应矩阵中。

10、完成步骤8后，需要再将每个多边形的信息矩阵作为单个元素保存至元胞数组中，这样最终就能得到保存了所有骨料信息的元胞数XY，通过XY_{i}可调用单个骨料整体矩阵信息，通过XY_{i}{j}可调用指定骨料指定点的坐标数据，方便后续工作。将此多边形的n个顶点坐标矩阵记录至步骤1建立的初始数组中，再将此数组记录至初始元胞数组中。完成此随机凸多边形的创建。

至此，上述9个步骤即为随机凸多边形骨料的生成基本步骤，对上述过程进行循环直至总面积大于设定粗骨料体积分数，

即：S≥r*X*Y，停止循环，进入投放环节，

其中，S表示骨料总面积。

本实施例提供的多边形随机骨料颗粒生成方法，采用随机凸多边形生成方法，相比于常规的延凸算法，具有以下优点：

1.算法简洁明了，计算效率高，生成结果符合现实混凝土骨料及其界面几何形状。

2.在程序的随机数生成过程中，加入了时间作为随机种子参与随机数的生成，达到了真正随机生成的目的。

3.相比常规延凸算法，本发明算法中点、线数据明显减少且可控，减少网格划分复杂度并方便后续转化为有效模型文件(.dxf)载入相关有限元软件中进行计算。

本实施例提供的多边形随机骨料颗粒投放方法，包括以下步骤：

获取通过多边形随机骨料生成方法形成的骨料信息，

判断生成的骨料是否达到预设的目标体积分数，如果否，则返回重新获取骨料；如果是，则对骨料按照以下步骤进行投放：

对骨料按照粒径d_i降序排列；

按照粒径从大到小的顺序投放骨料：

确定要投放的骨料，计算骨料与已经投放的骨料之间的圆心距，计算骨料与已经投放的骨料的外接圆半径之和，如果圆心距均大于外接圆半径之和，则直接投放，否则，再进行多边形重叠判别。

本实施例提供的多边形重叠判别包括三种情况的判断：点完全包含相交、点不完全相交、点不包含相交。

(a)点完全包含相交：该情况已经由降序排序的方法得以避免；

(b)点不完全包含相交：即保证新投放的凸多边形骨料的每个顶点位置都不在之前任何一个已经投放的骨料内部。假设多边形A_1-6后于多边形B_1-6投放，如图2a所示，图2a为点不完全包含相交的判别；采用MATLAB内建函数[in]＝inpolygon(xq,yq,xv,yv)，对A_i(i＝1,2,…,6)是否位于多边形B_1-6内部进行循环逐点判别；

(c)点不包含相交：需要在上述点判别的基础上，加入线相交判别。基本原理是通过判断凸多边形的一条边与其他多边形各个边的线段相交来进行判断，如图2b所示，图2b为点不包含相交判别，具体过程如下：

若线段A_iA_i+1、B_jB_j+1的顶点坐标为A_i(x_i,y_i)，A_i+1(x_i+1,y_i+1)以及B_j(x_j,y_j)，B_j+1(x_j+1,y_j+1)，通过建立线段A_iA_i+1、B_jB_j+1的线性方程，联立求解可得到交点坐标(x，y)为：

其中，x表示交点坐标的横坐标；y表示交点坐标的纵坐标；a_ii+1为第i个点和第i+1个点连线的梯度；b_jj+1为第j个点和第j+1个点连线的梯度。

得到交点坐标(x,y)后，判断交点与线段的位置关系，若同时满足以下条件，则可以判定这两条线段相交：

式中：

min表示坐标元素(x_i,x_i+1,x_j,x_j+1,y_i,y_i+1,y_j,y_j+1)最小值；

max表示坐标元素(x_i,x_i+1,x_j,x_j+1,y_i,y_i+1,y_j,y_j+1)最大值；

根据上述两个公式得到交点坐标公式和判定这两条线段相交公式，对投放后外接圆圆心距小于半径之和、满足点不完全包含相交情况的多边形进行逐边判别，若多边形所有边与其他已投放多边形边交点坐标都不能同时满足两条线段相交公式，表示此多边形不与其他任何多边形相交，可以按照此坐标位置投放。

本方法由于先按照粒径从大到小顺序进行投放，大部分骨料满足圆心距判别故无需再进行点、线相交的判别，这对模型程序整体的计算效率影响并不大。

本实施例提供的多边形随机骨料颗粒投放方法中还包括界面过渡区的生成，具体步骤如下：

将全部凸多边形随机骨料生成、投放完成后，根据各多边形骨料顶点按照逆时针的方向沿着连接的两条边外侧平移一段相等的距离，这个距离就是界面过渡区的厚度；

随后将平移后的点连接，得到界面过渡区的顶点分布。

本实施例提供的界面过渡区用于模拟：在水泥的水化反应，导致混凝土的骨料和水泥基之间形成了一层微观结构与水泥基材料不同的界面区域，有着更大的氢氧化钙含量和孔隙率。在混凝土细观稀物质传递研究过程中，界面过渡区可以被认为是包裹在骨料外部的一层均匀厚度的膜。

如图3所示，本实施例提供的基于多边形随机骨料的混凝土三相细观模型构建方法，

首先输入变量，所述变量包括混凝土试件长度X、混凝土试件宽度Y，骨料的目标体积分数r，骨料最小粒径D_min，骨料最大粒径D_max，界面过渡区厚度I_itz；界面过渡区扩散系数倍数a_itz，a_itz＝D_itz/D₀；其中，itz表示界面过渡区；Ditz为界面过渡区的氯离子扩散系数，D0为水泥浆体的氯离子扩散系数。

设置生成骨料粒径范围，最小粒径D_min，最大粒径D_max；

判断生成骨料粒径d_i的大小是在最大值D_max和最小值D_min之间，如果否，则返回重新生成骨料粒径；如果是，则以此d_i为骨料粒径生成凸多边形随机骨料；

按照以下公式计算骨料的体积分数V：V＝V+V_i；其中，V_i表示第i个骨料的体积分数；

判断骨料的体积分数V是否满足以下关系：V≥r*X*Y；如果否，则返回重新生成骨料粒径；如果是，则将生成的所有凸多边形骨料粒径及其顶点的点、线数据保存；

将d_i按照从大到小的顺序排列；

随机生成第i个骨料外接圆心坐标(x_i,y_i)；

判断生成的骨料是否有重叠发生，如果有，则重新随机生成新的骨料圆心坐标；如果否，则按顺序保存凸多边形骨料顶点坐标并连接；

将生成的骨料逐一投放直至全部结束，保存所有凸多边形骨料顶点点线数据至元胞数组f_xy中；

根据界面过渡区厚度将所有骨料边平移、连接生成界面过渡区；

保存所有界面骨料过渡区顶点点线数据至元胞数组f_itz中；

保存所有骨料、界面过渡区顶点点线数据至.dxf格式文件中。

实施例2

本实施例提供的基于多边形随机骨料的混凝土三相细观模型内氯离子侵蚀数值模拟方法，基于有限元软件COMSOL进行分析模拟，主要用于进行混凝土模型中氯离子扩散的细观数值模拟，因此主要需要用到化学物质传递模块中的“稀物质传递(tds)”模块。“稀物质传递”接口用于计算溶剂中稀溶质的浓度场，可以计算溶解在气体、液体或固体中的物质的传递和反应。传递的驱动力可以是菲克定律描述的扩散，与流体流动耦合的对流，以及与电场耦合的迁移。运用该模块的数值计算功能来研究混凝土中骨料、界面过渡区等影响因素对其内氯离子侵蚀的作用规律。稀物质传递模块的控制方程正是菲克第二定律，能够很好地模拟氯离子在混凝土中的侵蚀过程。混凝土三相细观模型内氯离子侵蚀数值模拟的具体过程如下：

1、创建仿真环境

首先，打开COMSOL Multiphysics软件，在模型向导的空间维度选项中选择二维，在物理场的选项卡中选择化学物质传递模块中的稀物质传递(tds)模块并将其添加到物理场接口中。点击完成，即可完成仿真环境的创建。

2、建立几何模型

在模型开发器中的【组件】/【几何】模块中，添加混凝土试件边界框(100×100mm²)。随后选择【导入】，在【源】选项卡中选择“DXF文件”或者“所有可导入文件”，在随后点击【浏览】，选择在MATLAB中创建的多边形骨料模型的.dxf文件并点击【导入】完成导入过程。重复上述过程，完成界面过渡区的导入过程。所述.dxf文件是根据混凝土随机凸多边形骨料模型程序的总体流程中得到所保存的混凝土随机凸多边形骨料信息的.dxf文件与相应的界面过渡区.dxf文件。

完成几何导入后对各个区域进行显式选择，分别建立混凝土细观模型的三相域，即水泥浆体域、骨料域以及界面过渡区域。最后将上述各部分构建成为联合体，完成几何模型的创建过程。如图4所示。图4为混凝土三相细观数值模型图。

3、设定基础参数

在【全局定义】/【参数】中，建立数值模拟计算中需要定义的基础参数，包括初始氯离子浓度C₀、水泥浆体中的氯离子扩散系数D₀、龄期衰减系数m、界面过渡区氯离子扩散系数倍数r等。

在【组件】/【定义】中，添加解析函数氯离子扩散系数时变模型D_t与表面氯离子浓度时变模型C_s，并设定相应的变元，并确定其参数的取值。

4、设定边界条件

基础参数设定完成后，需要对混凝土的封闭面及扩散面进行设置。如图5所示，图5为氯离子扩散边界设定；在研究一维扩散时，将左边框设定为表面氯离子浓度C_s的扩散面2，另外三个边框设定为无氯离子通量的封闭面1。

5、传递属性

由于骨料本身较水泥浆体和界面过渡区要更为致密，故假定氯离子不可透过骨料侵蚀混凝土，其内氯离子的浓度恒为0；水泥浆体与界面过渡区的氯离子扩散系数分别为D₀与r·D₀。基于COMSOL软件在设定通量条件时，模型本身具有的“层”的概念，因此在对模型的传递属性进行设置时可以通过层覆盖的方式来生成界面过渡区的层膜区域，具体顺序为：水泥浆体(底层)—界面过渡区(中层)—骨料(顶层)。如图6a、图6b、图6c所示，图6a为水泥浆体传递属性设定,图6b位界面过渡区传递属性设定,图6c为骨料传递属性设定。

6、划分网格

由于选择创建类似实际情况的界面过渡区的厚度，因此改变网格的细化程度会极大影响计算时间。本实施例选择自适应的三角细化网格进行划分。如图7所示，图7为混凝土数值模型网格划分图示。

7、计算及后处理

网格划分完成后，选择瞬态计算，设置计算步长，进行计算。计算完成后，COMSOL会根据结果自动生成默认时间下的氯离子浓度分布图，如图8所示，图8为混凝土细观数值模型中氯离子浓度分布图。

图8中展示了计算得到的氯离子浓度分布云图，其中浓度单位为％。可见氯离子浓度等值线在同一深度处不是直线，即同一深度处的氯离子含量并不相同，这与宏观匀质混凝土有显著区别。骨料的存在会使氯离子传输路径复杂化，对氯离子的传输产生明显的阻碍作用。这同样能够解释物理试验中有时测得的结果离散度较大的原因。

本实施例通过提取氯离子沿不同深度的浓度分布进行分析，在物理试验中一般是通过对混凝土对应深度的一层表面进行磨粉的方法来测定不同深度的氯离子平均浓度来实现。为了与试验结果相对比。本实施例采用混凝土的二维细观数值模型，在COMSOL中采取线平均的方法来求解不同深度平行于对应侵蚀面的对应截线的氯离子平均浓度。具体方法为：【结果】/【数据集】/【二维截线】，分别建立不同深度的二维截线，并在【派生值】/【线平均值】命令下进行不同深度二维截线的氯离子浓度线平均值计算。

实施例3

本实施例建立的模型为随机凸多面形骨料的混凝土三相细观数值模型，其点线信息较规则的圆形、正多边形等骨料要更为复杂，基于COMSOL本身支持导入模型文件的功能，通过建立MATLAB生成模型文件与COMSOL软件之间的交互程序。

在实现用MATLAB生成及投放随机凸多边形骨料程序的基础上，对DXFLib函数库中部分主体函数进行了修改和引用。使其更符合二维混凝土随机凸多边形骨料细观模型的保存要求，可实现将随机凸多边形骨料及其界面过渡区模型文件转换为.dxf格式文件。

本实施例中采用COMSOL Multiphysics 5.3版软件来实现有限元模型计算过程，为了提高COMSOL模型方法的数值计算效率，在其内APP方法模块下新增添了模型方法，基础语言环境为Java。通过该方法可以实现对软件中的大部分操作过程进行自定义重复操作，从而大幅加快工作流程，提高计算效率。

建立模型方法是：利用录制方法功能对图形用户界面(GUI)中执行的一系列操作进行记录，随后通过方法重放便可以进行相同的操作步骤，记录保存对模型文件中的所有操作记录，从操作记录的代码中提取相关的操作，对相关操作进行取舍后，通过历史命令压缩后将Java命令保存为记事本文件，插入到新的模型方法中去。

由于COMSOL Multiphysics会对模型文件中的所有操作记录保存为Java代码。，因此可直接从代码中提取相关的操作，进行取舍后插入到新的模型方法中去，具体方法如下：

通过【文件】菜单下的【压缩历史命令】选项来进行实现。对历史命令压缩后，打开保存Java命令的记事本文件，将位于：“Model model＝ModelUtil.create("Model")；”与“return model；”之间的代码进行保存，并复制到COMSOL中新建立的模型方法中。

对其中的代码进行取舍后，留下需要重复执行的操作代码行，对文件名进行变量的赋值，通过循环实现模型文件的自导入、划分网格与计算的自运行、结果的自导出等等一系列操作的重复执行。

通过建立专门的模型方法，仅通过一定的适应性修改，作为一种可重复使用的模板来应用至不同的研究环节中，比如：自定义初始条件、自定义物理场接口、自定义函数表达式等等。又或者是创建自定义的研究目的，自定义几何文件的导入过程、自定义求解设置、自定义表格的创建与计算结果的导出等等，方便重复使用。

本实施例采用MATLAB建立的混凝土随机凸多边形骨料细观模型的算法基础及原理，以及在此基础上结合COMSOL Multiphysics物理仿真计算软件提出了相应的数值模拟方法：

结合开源函数库DXFLib成功将模型文件转换为可被COMSOL软件识别、导入的.dxf格式模型文件，实现使用MATLAB建立数值模型与有限元计算软件的交互；

基于MATLAB生成多边形骨料模型，在COMSOL Multiphysic的稀物质传递模块下进行混凝土三相细观数值模型中氯离子侵蚀数值模拟。在Java语言环境下开发COMSOL模型方法，可根据不同需求对程序进行自定义运行，实现有限元计算自动运行、计算数据自动导出等功能，大幅提升计算及后处理效率。

以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。

Claims (10)

1.基于多边形随机骨料的混凝土三相细观模型构建方法，其特征在于：包括以下步骤：

S1：设置混凝土数值模型的变量；

S2：根据混凝土数值模型的变量生成凸多边形随机骨料；

S3：按照以下公式计算已生成骨料的累积体积分数V：V＝V+V_i；其中，V_i表示第i个凸多边形随机骨料的体积分数；

S4：判断已生成的骨料的累积体积分数V是否满足以下关系：V≥r*X*Y；如果否，则返回继续生成凸多边形随机骨料；如果是，则将生成的所有凸多边形骨料粒径d_i及凸多边形骨料顶点的点线数据保存；其中，r表示生成凸多边形骨料的目标体积分数，X表示拟构建混凝土数值模型的长度，Y表示拟构建混凝土数值模型的宽度；

S5：将已生成好的凸多边形骨料按照骨料粒径d_i从大到小的顺序进行降序排列；

S6：随机生成第i个凸多边形骨料外接圆心坐标(x_i,y_i)；

S7：随机生成第i+1个凸多边形骨料外接圆心坐标(x_i+1,y_i+1)；

S8：根据凸多边形骨料外接圆心坐标(x_i,y_i)、(x_i+1,y_i+1)来判断第i个和第i+1个凸多边形骨料是否会发生重叠。若重叠，则重新随机生成新的凸多边形骨料外接圆心坐标；若不重叠，则按顺序保存满足要求的凸多边形骨料顶点坐标并将各顶点相连接，并将这第i个、第i+1个凸多边形骨料按照满足要求的外接圆心坐标(x_i,y_i)、(x_i+1,y_i+1)投放到混凝土数值模型的X*Y区域当中；

S9：将S5中生成好的凸多边形骨料按照S8的重叠判别准则逐一投放直至全部结束后，将所有凸多边形骨料顶点数据保存至元胞数组f_xy中；

S10：将全部已成功投放的凸多边形骨料各边按照界面过渡区的厚度由内向外平移，平移后的各边沿其直线方向相延伸连接，形成紧邻凸多边形骨料各边的等厚界面过渡区相；

S11：将全部生成好的界面过渡区顶点的点、线数据保存至元胞数组f_itz中；

S12：将满足S9和S10步骤的全部凸多边形骨料及其界面过渡区顶点的点、线数据保存至.dxf格式文件中。

2.如权利要求1所述的方法，其特征在于：所述步骤S2中采用多边形随机骨料颗粒生成方法，所述多边形随机骨料颗粒生成方法，具体按照以下步骤：

S21：设置混凝土数值模型的初始条件，包括混凝土模型的长宽尺寸X、Y，多边形随机骨料的目标体积分数r，骨料最小粒径D_min，骨料最大粒径D_max；

S22：生成第i个随机粒径为d_i的多边形骨料；

S23：判断多边形骨料粒径d_i是否满足D_min≤d_i≤D_max，若不满足，则重新执行步骤S22重新生成第i个随机粒径为d_i的多边形骨料，若满足，则进行下一步；

S24：生成直径为d_i的多边形骨料的外接圆，在外接圆上任取三点形成锐角三角形；

S25：储存锐角三角形的点线信息，计算锐角三角形的面积s；

S26：随机生成一个外接圆内的点，如果生成的点没有在多边形外部，则重新生成；如果生成的点在多边形外部，则将新生成的点与原多边形顶点联合形成新的多边形，新生成的多边形为凸多边形；

S27：对新生成多边形各顶点坐标进行重新排序、记录各点坐标，并计算多边形骨料的累加面积；

S28：设定多边形骨料边数n，重复循环至步骤S22开始执行，直到完成设置的边数n为止；

S29：将生成的所有多边形骨料的数据保存在对应矩阵中；

S210：将每个多边形的信息矩阵作为单个元素保存至元胞数组中，得到所有骨料信息的元胞数组XY。

3.如权利要求2所述的方法，其特征在于：所述步骤S26的凸多边形的判定，具体按照以下步骤进行：

判断新生成的点是否位于原三角形或多边形外部；以及新生成的点与原三角形或多边形组合形成的多边形为凸多边形；

如果否，则将返回继续判断；如果是，则进行下一步；

按照以下公式计算三角形面积：

式中，(x₁,y₁)、(x₂,y₂)以及(x₃,y₃)为三角形的三个顶点的坐标，且三个点是按逆时针排序；

计算新点P与相邻老点按逆时针顺序组合成的两个三角形的面积，如果所述三角形面积的结果都大于0时，则新生成的多边形为凸多边形。

4.如权利要求1所述的方法，其特征在于：所述步骤S9中的多边形骨料在进行投放的过程中还包括以下步骤：

S91：计算确定要投放的多边形骨料与已经投放的多边形骨料之间的圆心距；

S92：计算多边形骨料与已经投放的多边形骨料的外接圆半径之和；

S93：如果圆心距均大于外接圆半径之和，则直接投放，否则，再进行多边形重叠判别。

5.如权利要求4所述的方法，其特征在于：所述步骤S93中的多边形重叠判别还包括点不完全包含相交判断，具体如下：

新投放的凸多边形骨料的每个顶点位置都不在之前任何一个已经投放的骨料内部。

6.如权利要求4所述的方法，其特征在于：所述步骤S93中的多边形重叠判别还包括点不包含相交判断，所述点不包含相交判断是通过判断凸多边形的一条边与其他多边形各个边的线段相交来进行判断，具体过程如下：

假设A_i和A_i+1构成的线段A_i,i+1，B_j和B_j+1构成的线段B_j,j+1，A_i,i+1和B_j,j+1二者的交点坐标(x,y)计算式如下：

其中，线段A_i,i+1的顶点坐标分别为A_i(x_i,y_i)、A_i+1(x_i+1,y_i+1)；线段B_j,j+1的顶点坐标分别为B_j(x_j,y_j)、B_j+1(x_j+1,y_j+1)；a_ii+1为第i个点和第i+1个点连线的梯度；b_jj+1为第j个点和第j+1个点连线的梯度；

得到交点坐标(x,y)后，判断交点与线段的位置关系，若同时满足以下条件，则可以判定这两条线段相交：

式中：

min表示坐标元素(x_i,x_i+1,x_j,x_j+1,y_i,y_i+1,y_j,y_j+1)最小值；

max表示坐标元素(x_i,x_i+1,x_j,x_j+1,y_i,y_i+1,y_j,y_j+1)最大值；

根据上述两个公式得到交点坐标公式和判定这两条线段相交公式，对投放后外接圆圆心距小于半径之和、满足点不完全包含相交情况的多边形进行逐边判别，若多边形所有边与其他已投放多边形边交点坐标都不能同时满足两条线段相交公式，表示此多边形不与其他任何多边形相交，则按照此坐标位置投放。

7.如权利要求1所述的方法，其特征在于：所述步骤S10中多边形随机骨料颗粒还包括界面过渡区的生成过程，具体步骤如下：

将全部凸多边形随机骨料生成、投放完成后，根据各多边形骨料顶点按照逆时针的方向沿着连接的两条边外侧平移一段相等的距离，所述距离为界面过渡区的厚度；

随后将平移后的点连接，得到界面过渡区的顶点分布。

8.如权利要求1所述的方法，其特征在于：所述步骤S11是按照以下方式进行保存：

S111：保存所有界面骨料过渡区顶点的点线数据至元胞数组f_itz中；

S112：保存所有骨料、骨料界面过渡区顶点的点线数据至.dxf格式文件中。

9.基于多边形随机骨料的混凝土三相细观模型中氯离子侵蚀数值模拟方法，其特征在于：包括以下步骤：

S1：确定创建仿真环境；

S2：在仿真环境中建立混凝土试件几何模型，具体如下：导入已经创建的包含全部骨料及其界面过渡区信息的.dxf文件；所述.dxf文件是根据混凝土随机凸多边形骨料模型构建方法中得到的随机凸多边形骨料信息以及相应的界面过渡区信息的.dxf文件；建立混凝土细观数值模型的三相域，包括：水泥浆体、骨料及界面过渡区，并构建成为联合体，完成几何模型的创建过程；

S3：设定数值模拟计算中需要定义的基础参数，所述基础参数包括初始氯离子浓度C₀、水泥浆体氯离子扩散系数D₀、龄期衰减系数m、界面过渡区扩散系数倍数a_itz；

设置解析函数扩散系数时变模型D_t与表面氯离子浓度时变模型C_s，并设定相应的变元，并确定其参数的取值；

S4：设定边界条件：

基础参数设定完成后，对混凝土数值模型的封闭面及扩散面进行设置；所述封闭面为不允许氯离子透过或浸入的表面，所述扩散面为允许氯离子透过或浸入的表面；

S5：设定混凝土细观数值模型中各组成相的传递属性；

S6：划分网格：采用自适应的三角细化网格划分；

S7：数值模拟计算及后处理：

网格划分完成后，选择瞬态计算，设置计算步长，进行计算；

计算完成后，生成整体氯离子浓度分布图。

10.如权利要求9所述的方法，其特征在于：所述骨料模型.dxf文件是通过以下步骤来形成的：

基于MATLAB生成骨料模型：录制骨料模型建立过程中的系列操作；

从系列操作中提取骨料模型建立的操作步骤；

对操作步骤进行取舍：保留需要的操作步骤的操作代码，并整理为可重复执行的骨料模型.dxf文件。

Images