材料二维微结构仿真随机概率圆优化填充算法

陈　震， 周金宇

(江苏理工学院常州市装备再制造工程重点实验室，江苏 常州 213001)

材料二维微结构仿真随机概率圆优化填充算法

陈震， 周金宇

(江苏理工学院常州市装备再制造工程重点实验室，江苏 常州 213001)

几何建模是影响有限元法计算材料性能精确性的主要因素之一。利用圆等简单的几何实体近似模拟离散化实体效果良好，不仅能大量简化几何计算，还能构造出误差在可接受范围内的近似模型。材料特性不仅取决于平均颗粒尺寸，受颗粒尺寸分布的影响也十分显著。本文设计的随机圆填充算法保证了填充圆半径服从某一概率分布，重点分析了常规圆产生方法和依据干涉圆的干涉类型生成新圆。该算法不仅完全杜绝了圆的干涉还极大的减小了圆的相离，相对于现有的填充算法，填充度有了显著提升。

随机半径；材料仿真；颗粒填充；有限元法

使用离散化单元法分析问题首先要采用离散化实体划分的形式近似代替各种实际情况[1-2]。对于许多复杂的物理现象，利用圆等简单的几何实体近似模拟离散化实体效果良好。圆填充问题也由其学术价值和工业应用的重要性成为几个世纪以来长期研究的课题。17世纪英国数学家牛顿就曾研究过关于等尺寸圆彼此不干涉时在某几何形状中的填充以获得最大填充度的最优几何填充问题。该类算法不仅能简化几何计算，还能构造出误差区间在可接受范围内的模型。模型使得复杂的几何结构变成简单的几何填充问题。随着20世纪50年代计算机技术的迅速发展，计算机技术在材料仿真建模方面得广泛应用[3]。

有限元技术作为解决材料力学问题最有力的数字化方法之一，其精确性主要与被分析材料的几何模型有关[4]。在实际分析应用中，尽量提高离散化实体填充度的同时，还需控制其概率分布对材料性能的影响[5]。

材料微结构模型包括二维和三维两类。二维模型采用圆填充、直线划分圆心边界的形式；三维模型则采用球填充、平面划分球心边界的形式，当前由空间点集生成三维模型的方法主要有增量扩展算法[6]和基于 Delaunary三角剖分算法[7]。二维微结构模型是三维模型的基础，通过对二维模型的有效扩展可以构建出三维模型。基于二维微结构模型的细观有限元分析，对评估材料表面微裂纹的萌生和扩展，推演二维微结构复合材料损伤演化行为具有重要意义[8-9]。

Voronoi算法是常见的随机晶粒生成算法[10]，其定义为：设P={p1, p2, ···, pn}是二维欧氏几何空间平面上的n个互不相同的点集，d(p, pi)表示点p与pi之间的欧几里德距离，由划分的平面称作以为点集的 Voronoi图，V(pi)为Voronoi晶胞，pi为各晶胞的生成子。Voronoi多采用在某平面根据生成晶胞的平均面积确定生成子数量，然后随机撒点并划分点的边界生成晶胞，该方式比较容易控制平均晶胞面积尺寸，但无法控制尺寸分布使其满足多样化的统计规律[5]，具有单一性和不可调节性，且主要适应于晶粒均匀形核并等速生长的退火状态纯金属，应用范围较窄。

故有学者采用 Voronoi图的推广 Laguerre图[5,11]进行建模。Laguerre图多采用增量式算法生成点集，然后通过对这些点集赋予权值形成权圆，但由于点位置的随机性就决定了权圆会发生干涉或相离的情况。干涉会减小权圆面积，相离则会产生过多的未填充区域，导致后继的多边形边界划分失真，产生较大的材料建模误差，如图 1所示。

图1　二维Laguerre图[11]

本文提出的圆填充算法能根据给定的概率密度函数生成随机圆填充矩形区域。文中通过生成随机半径、排除干涉逐层生成的方式，可以实现半径依据要求概率密度函数生成概率圆且没有干涉。同时，由于每个生成的新圆至少与已存在的两个圆或边界相切，故填充度明显高于传统Laguerre权圆填充度，只要根据相关材料的边界性质进行边界划分，即可构建较高精度的二维材料微结构仿真模型。

1　常规颗粒圆生成及活动链更新

该算法采用的是自底向上逐层生成的方法。假设填充区域的左、右、上、下边界分别为Bl，Br，Bt，Bb。首先根据指定的概率密度函数取随机半径r1产生圆C11，使C11同时与边界Bl和Bb相切，然后依次生成C12, C13, ···, C1n，使新生成的圆分别与前一个圆和边界Bb相切直至新圆Cn与边界Br干涉，再以r1n为半径生成新圆使其同时与C1(n-1)和边界 Br相切取代圆 C1n，至此第一层圆生成完毕，同时形成初始链 C11→C12→C13→…→C1(n–1)→C，如图2(a)所示。该链是由具有确定方向的线段 C11→C12, C12→C13,…, C1(n–1)→C1n’组成并可以随着新圆的产生而更新。

初始链形成之后，根据指定的概率密度函数产生随机半径r21作新圆C21，使其同时与C11和边界Bl相切，然后添加新的有向线段C21C11，更新成新的活动链C21→C11→C12→C13→…→C1(n–1)→C1n’，如图2(b)所示。然后按序首先选取C11C12作为当前线段，产生随机半径 r22作新圆 C22，要求该圆与C11和C12相切并居于有向线段C11C12的左侧，然后从初始活动链上删除线段C11C12，并添加新的有向线段 C11C22和 C22C12，更新成新的活动链C21→C11→C22→C12→C13→…→C1(n–1)→C1n’，如图2(c)所示。下一步选取线段C12C13作为当前线段，以相同的方式产生C23，并从初始活动链上删除线段C12C13，再添加新的有向线段C12C23和C23C13，形成新的活动链C21→C11→C22→C12→C23→C13→…→C1(n-1)→C1n’，如图 2(d)所示。依此逐步生成新圆并更新活动链直至本层圆生成完毕。

图2　新圆生成及活动链更新

2　颗粒圆干涉的处理方法

2.1新开始层第1个圆的干涉及活动链更新

新开始层的第 1个圆可能会与已生成的圆存在干涉。如图 3(a)所示，活动链为 C11→C12→C13→…，以r*为半径产生的圆Ci同时与C11和边界 Bl相切且与 C12干涉，则再次以随机半径 r*生成新圆Cz，要求该圆同时与C12和边界Bl相切。再次判断该圆是否与活动链上的圆干涉。如果不干涉，则当前活动链更新为Cz→C11→C12→C13→…；如果依然存在干涉，则活动链不更新。然后按序首先选取 C11C12作为当前线段，进一步生成新圆。

2.2与边界Br干涉及活动链更新

随着沿活动链正方向按序选取当前线段形成新圆，该圆可能会与边界Br干涉。如图3(b)所示，活动链为…→Cp(k–2)→Cp(k–1)→Cpk→C2m→C1n，以C2m→C1n为当前线段，r*为半径产生的圆Ci与边界Bl干涉，则以随机半径 r*生成新圆 Cz，首先要求该圆与C2m和边界Bl相切，判断该圆是否与活动链上的圆干涉。如果不干涉，则当前活动链更新为…→Cp(k–2)→Cp(k–1)→Cpk→C2m→Cz，开始新层；如果存在干涉，再次以随机半径 r*生成新圆 Cz，并要求该圆与 C1n和边界 Bl相切，判断该圆是否与活动链上的圆干涉。如果不干涉，则当前活动链更新为…→Cp(k–2)→Cp(k–1)→Cpk→C2m→Cz，开始新层；如果依然存在干涉，则活动链不更新，按序选取下一线段作当前线段。由于图 3(b)中C2m→C1n为活动链的最后一个线段，所以本层结束，开始新层。

因此，中国版南海争端国际意象的塑造和推广还需通过具体的外交政策和精准的国际舆论引导才能真正获得相关国家的认同。一方面，使建立在客观、尊重历史源流、符合国际法的公平公正基础上的“中国版”南海问题国际意象为各方所周知乃至认同；另一方面，对于在南海问题上的误导性意象必须予以坚决回应乃至驳斥以正视听。

图3　新层首圆及新圆与边界Br的干涉处理

2.3根据当前线段生成的新圆只与已存在圆的干涉及处理

当沿活动链正方向前进，根据当前线段生成的新圆可能会只与已生成的圆存在干涉。此时沿当前线段正方向前进，定义当前线段前方的圆为前置圆，其后方的圆为后置圆。处理此干涉类型需要选取新的当前线段替代原当前线段，可分为以下2种干涉情况：

(1) 与前置圆干涉。如图4(b)所示，以CkCk+1为当前线段随机生成的圆 Ci(半径为 r*)与前置圆Ck+n干涉，则依次生成当前线段CkCk+n，Ck+1Ck+n。根据新的当前线段 CkCk+n以随机半径 r*生成新圆Cz，判断该圆是否与活动链上的圆或边界干涉。如果不干涉，则在当前活动链中删除线段 CkCk+1, Ck+1Ck+2,…, Ck+n–1Ck+n，并添加新的有向线段CkCz和CzCk+n，形成新的活动链；如果干涉，则再选取当前线段 CkCk+n，同理生成新圆并检查干涉，如果不干涉，则在当前活动链中删除线段Ck+1Ck+2,…, Ck+n–1Ck+n，并添加新的有向线段Ck+1Cz′和Cz′Ck+n，形成新的活动链；如果依然干涉，则跳过原当前线段 CkCk+1，选取下一当前线段Ck+1Ck+2。

图4　干涉情况

(2) 与后置圆干涉。如图4(c)所示，以CkCk+1为当前线段随机生成的圆 Ci(半径为 r*)与前置圆Ck–n干涉，则将当前线段更新为 Ck–nCk+1，根据新的当前线段 Ck–nCk+1再次以随机半径 r*生成新圆Cz，判断该圆是否与活动链上的圆或边界干涉。如果不干涉，则在当前活动链中删除线段Ck–nCk–n+1,…, Ck–1Ck, CkCk+1，并添加新的有向线段Ck-nCz和CzCk+1，形成新的活动链；如果干涉，则再选取当前线段 Ck–nCk，同理生成新圆并检查干涉，如果不干涉，则在当前活动链中删除线段Ck–nCk–n+1,…, Ck–1Ck，并添加新的有向线段Ck–n和′Ck，形成新的活动链；如果依然干涉，则跳过原当前线段CkCk+1，选取下一当前线段Ck+1Ck+2。

3　扩展链算法

构造随机权圆扩展链算法的主流程如图 5所示，具体步骤如下：步骤1. 以initial_chain子程序生成初始圆。步骤2. 判断是否满足区域的填充要求。如果满足则进入步骤17。

步骤3. 判断是否是新层。如果不是则进入步骤15。

步骤4. 以random_r子程序生成随机半径r*，生成新层首圆。

步骤6. 更新活动链。

步骤7. 接受新圆并用circle子程序生成该圆。步骤8. 按序选取活动线段，回到步骤2。

步骤9. 更新干涉矩阵，即把与新圆干涉的圆数据依次存入干涉矩阵，数据包括各圆心坐标及其半径，如图4(b)所示，以CkCk+1为当前线段生成的新圆Ci与Ck+n等s个圆干涉，则把这s个圆Ck+n, Ck+n–1, … , Ck+n–s+1的数据依次存入干涉矩阵。

步骤10. 按序读取干涉矩阵中各圆心坐标元素。

步骤11. 根据第1节中干涉类型生成新圆并判断其是否与活动链上已存在的圆或边界干涉，如果不干涉则进入步骤14。

步骤12. 判断标示量t是否为1，如果是，则进入步骤8。

步骤13. 判断干涉矩阵是否读取完毕，如果没有，回到步骤10；如果读完，回到步骤8。

步骤14. 根据第1节中各干涉类型和情况更新活动链，回到步骤7。

步骤15. 生成新圆并判断其是否与边界Br干涉，如果干涉则另标示量t=1，进入步骤11。

步骤 16. 判断新圆是否与活动链上圆干涉，如果干涉进入步骤9；否则进入步骤6。

步骤17. 程序结束。

图5　随机颗粒圆生成算法流程图

算法采用自下向上逐层生成的方法产生概率圆进行区域填充。算法执行产生 n层随机圆；每层产生随机圆的计算次数为mi(i=1, 2, 3,…, n)，mi为第i层概率圆数量；其中第i层有ki个圆需要进行干涉处理，则运行上限为n*mmax*kmax次；常规概率圆的生成运行上限为n*mmax次；此外，循环外附加的初始链生成算法、原始条件赋予和存储矩阵建立等虽然会带来额外的计算量，但相对于循环内考虑概率圆干涉情况的计算量几乎可以忽略，不影响实际复杂度。综上所述，扩展链算法的时间复杂度为O(n3)。

4　数值算例

基于 4个指定概率密度函数进行概率圆填充实验，概率密度函数分别为：

根据指定概率密度函数生成的概率圆分别如图6(a)～(d)所示。各圆之间完全杜绝了干涉现象，且每个圆都至少与两个圆相切，填充度高。表 1中的数据1～4依次对应4种概率密度函数。期望值为概率密度函数确定的圆半径统计参数，实际值为根据相应概率密度函数生成的圆半径的统计参数，填充度为区域内圆的面积之和与矩形区域面积的百分比值。

图6　根据概率分布生成的概率圆

表1　概率圆统计参数表

5　结 束 语

本文详细叙述了圆填充算法的实现过程，该算法可以根据给定的概率密度函数生成随机圆填充区域。根据给定的概率密度函数生成随机半径、排除干涉逐层生成新圆。算法避免了生成过程中圆的干涉情况，并且生成圆的半径完全符合要求概率密度函数。同时，由于新圆在生成时至少与已生成的两个圆或边界相切，故填充度较现有算法有明显提升，提高了仿真建模应用的可靠性和真实性。

[1] 高井望, 徐佩华, 黄润秋, 等. DEM 拟合级配碎石材料大三轴试验的颗粒粗糙度效应研究[J]. 工程地质学报, 2013, 21(4): 560-569.

[2] Gong G B, Zha X X, Wei J. Comparison of granular material behaviour under drained triaxial and plane strain conditions using 3D DEM simulation [J]. Engineering Computation, 2012, 25(2): 186-196.

[3] 秦湘阁, 刘国权. 颗粒增强复合材料显微组织的计算机仿真[J]. 佳木斯大学学报, 2001, 19(1): 1-3.

[4] 崔文凯, 冯仰德, 纪国良, 等. 两种复合材料几何建模算法[J]. 计算机辅助设计与图形学学报, 2015, 27(1): 46-50.

[5] 李俊琛. 复合材料微观组织结构的仿真与计算[D].兰州: 兰州理工大学, 2009.

[6] 孙殿柱, 刘健, 李延瑞，等. 三维散乱点集Voronoi图快速生成算法研究[J]. 武汉科技大学学报, 2010, 35(8): 909-912.

[7] 范志刚, 吴裕功, 赵选贺. Laguerre-Voronoi图软件包的设计和实现[J]. 天津大学学报, 2003, 36(6): 747-752.

[8] 任淮辉, 李旭东. 二维材料微结构设计与数值模拟软件系统[J]. 计算机辅助设计与图形学学报, 2009, 32(1): 45-50.

[9] 任淮辉, 李旭东, 李俊琛, 等. 二维多晶体材料微结构的力学响应计算[J]. 武汉科技大学学报, 2008, 31(4): 427-431.

[10] 司良英, 邓关宇, 吕程, 等. 基于Voronoi图的晶体塑性有限元多晶几何建模[J]. 材料与冶金学报, 2009, 8(3):193-216.

[11] 张赋, 李旭东. 带权点集Laguerre图的增量算法与软件设计研究[J]. 计算机工程与应用, 2012, 48(30): 10-18.

Random Circles Optimization Filling Algorithm of Material Two Dimensional Microstructures Simulation

Chen Zhen,Zhou Jinyu
(Changzhou Key Laboratory of Equipment Remanufacturing Engineering, Jiangsu University of Technology, Changzhou Jiangsu 213001, China)

Geometric modeling is one of the main factors affecting the accuracy of material properties calculation by finite element method. Using simple geometric entities, such as circles, to simulate discrete objects is demonstrated well. Beyond the great simplifications in the geometrical calculation, it can provide an approximate model within acceptable error. Material properties not only depend on the average grain size, but also are influenced by the grain size distribution significantly. Random circles packing algorithm designed in this paper ensured these circles′ radius following a certain probability distribution. Normal circles generation method and generating new circle based on the interference types were analyzed. This algorithm not only avoided the interference completely, but reduced the separation greatly. Compared with existing filling algorithms, the filling density of this algorithm was improved obviously.

random radius; material simulation; sphere packing; finite element method

TP 391.9

A

2095-302X(2015)06-0944-06

2015-05-19；定稿日期：2015-08-03

国家自然科学基金资助项目(51275221)；江苏省产学研联合创新资金资助项目(BY2014038-04)

陈震(1990–)，男，江苏徐州人，硕士研究生。主要研究方向为先进(再)制造技术与装备。E-mail：ericchen31@163.com

周金宇(1973–)，男，江苏扬州人，教授，博士。主要研究方向为结构可靠性、机械疲劳与延寿、现代设计方法。E-mail：yuhangyuan888@sina.com