【理论算法】基于Voronoi方法的柱状节理岩生成方法

火成玄武岩在熔岩流冷却过程中通常形成柱状节理结构。柱状节理结构的存在通常使玄武岩表现出较差的工程地质和力学性能。因此，柱状节理玄武岩通常被认为是岩土工程中的一种软弱岩体，需要进一步研究和大量的工程处理来改善其性质。对柱状节理玄武岩力学性能的研究可以追溯到20世纪80年代美国的玄武岩废弃物隔离项目（BWIP）。

生成具有特定拓扑和统计特性的柱状节理结构在建模中具有重要意义。Voronoi图常用于生成节理岩体和结晶岩石结构。有几个参数，如平均边数、两个种子之间的最小距离和能量，来描述Voronoi图的几何性质。然而，这些参数本身并不能代表不同Voronoi多边形的可变性。为了更好地描述，本工作采用了Duyckaerts和Godefroy使用的变异系数[1]。提出了节理密度、变异系数和节理厚度来描述柱状节理结构的拓扑和统计特性。

对于表示柱状节理结构的Voronoi算法，节理密度        定义为

其中        是Voronoi多边形的平均面积，可以计算为

其中        是第        个多边形的面积，        是多边形的数量。

变异系数        被定义为描述不同Voronoi多边形的变异性：

其中        是Voronoi多边形面积的标准变化:

基于柱状节理结构的定量描述，本文提出了一种基于CCVT的改进算法以生成基于Lloyd迭代和平分方法的RVE模型[2]。在传统的Lloyd迭代方法中，会生成质心图，该质心图过于规则，无法描述柱状玄武岩。在这种新方法中，可以通过改变迭代过程来生成具有指定变异系数的Voronoi图。

该过程可分为以下步骤：

（1） 基于现场观测，得到柱状节理结构的参数，即节理密度        、变异系数        和节理厚度        ，并设置阈值，即收敛误差误差errorc和边长阈值误差errore。

玄武岩柱状节理结构（a） 柱状节理的地质图像，（b）柱状节理草图

（2） 基于关节密度        和域面积        生成具有随机位置的N个点。

（3）使用种子        、相应的变异系数        和质心        计算约束Voronoi图。

（4）如果计算的变化系数        和原位值        之间的误差小于errorc，则进入步骤5；否则，转到步骤6。

（5） 如果        大于        ，则将临时位置        更新为        =        ，然后将质心        设置为下一次迭代的新种子。否则，将临时位置        和种子        的中点设置为新种子。然后转到步骤2进行下一次迭代；

（6） 将生成的种子转换到8个相邻域，并使用新的种子和域边界生成Voronoi图；

（7） 使用正则化过程删除比errore短的边，并基于关节厚度偏移Voronoi图的多边形，考虑自相交。

（8） 使用原始边界修剪多边形，获得具有周期性结构的RVE模型。

基上述理论开发了VoroRock程序并在International Journal of Rock Mechanics and Mining Sciences（JCR Q1，IF=6.849）上发表”Numerical homogenization study on the effects of columnar jointed structure on the mechanical properties of rock mass” (doi.org/10.1016/j.ijrmms.2019.104127)。此外，相关代码已在github上开源，地址为https://github.com/GeoGroup/VoroRock.git。

[1] Duyckaerts C, Godefroy G. Voronoi tessellation to study the numerical density and the spatial distribution of neurones. J Chem Neuroanat. 2000;20(1): 83-92.https://doi.org/10.1016/S0891-0618(00)00064-8

[2] Du Q, Faber V, Gunzburger M. Centroidal Voronoi tessellations: Applications and algorithms. Siam Rev. 1999;41(4): 637-676. https://doi.org/10.1137/S0036144599352836