在abaqus中生成voronoi多面体的方法

在abaqus中生成voronoi多面体的方法本文转自技术邻张磊这篇小文旨在介绍在abaqus中生成Voronoi凸多面体的方法，在模拟晶体或类似结构时，经常需要生成许多相互连接的Voronoi多面体，如图1所示。通常方法是利用MATLAB或是其他软件生成多面体的空间结构，然后导入abaqus进一步处理。此外还有国外学者开发的neper及其衍生软件可以使用，不过该软件目前仅支持Linux操作系统，而且往往需要用户做进一步的开发才能满足该用户特定的需求，对于使用者的编程能力要求较高。图1Voronoi多边形示意图其实生成Voronoi多面体这一过程不需借助第三方软件，在abaqus中合理利用脚本也可以实现。Abaqus常用的脚本语言是Python，Python是一种简单易学的编程语言，而且abaqus脚本学习起来也很轻松，因为用户在abaqus的图形界面里每做一个操作，rpy文件就会记录下对应的脚本命令（一些特殊的操作除外），再配合上help文档的详细说明，很多新手也可以在短时间内运用脚本来解决一些仿真中比较繁琐的问题。Python库Scipy提供了现成的Voronoi类，和MATLAB里的函数一样，可以提供voronoi各个单元的空间点坐标，有关Voronoi类的调用，所需变量及其属性，可参考如下，其中属性部分我用中文做了介绍：scipy.spatial.Voronoi(points,furthest_site=False,incremental=False,qhull_options=None)Parameters:points:ndarrayoffloats,shape(npoints,ndim)Coordinatesofpointstoconstructaconvexhullfromfurthest_site:bool,optionalWhethertocomputeafurthest-siteVoronoidiagram.Default:Falseincremental:bool,optionalAllowaddingnewpointsincrementally.Thistakesupsomeadditionalresources.Qhull_options:str,optionalAdditionaloptionstopasstoQhull.SeeQhullmanualfordetails.(Default:“QbbQcQzQx”forndim4and“QbbQcQz”otherwise.Incrementalmodeomits“Qz”.)注：对于一般用户，只需要注意points参数即可，代表voronoi的种子点，可以是N维数组，只不过超过3维的数组就没有空间上的意义了。（在我们这个宇宙是这样的，对于高维宇宙可能并不成立^_^）Attributespoints(ndarrayofdouble,shape(npoints,ndim))Coordinatesofinputpoints.核心点坐标，n行2列的数组vertices(ndarrayofdouble,shape(nvertices,ndim))CoordinatesoftheVoronoivertices.多边形顶点的坐标，n行2列的数组，只标出有限的顶点，处于无限远的顶点不会标出ridge_points(ndarrayofints,shape(nridges,2))IndicesofthepointsbetweenwhicheachVoronoiridgelies.输出的是点信息,每条voronoi边所处的两个核心点的序号ridge_vertices(listoflistofints,shape(nridges,*))IndicesoftheVoronoiverticesformingeachVoronoiridge.多边形每条边的端点序号，如果端点是无限远，那么序号就是-1regions(listoflistofints,shape(nregions,*))IndicesoftheVoronoiverticesformingeachVoronoiregion.-1indicatesvertexoutsidetheVoronoidiagram.给出了每个多边形的端点信息，如果端点是无限远，那么序号就是-1，有多个无限远端点的，都只用一个-1point_region(listofints,shape(npoints))IndexoftheVoronoiregionforeachinputpoint.Ifqhulloption“Qc”wasnotspecified,thelistwillcontain-1forpointsthatarenotassociatedwithaVoronoiregion.对于每一个输入点，给出其对应的region序号注：ridge_vertices这个属性的名字有一定误导性，对于2维voronoi来说，这个属性指的是边上的顶点序号，对于3维voronoi来说，这个变量实际上指的是面上的顶点序号，对于N维voronoi而言，这个变量指的是N-1维边界上的各个顶点。参考出处：~pbechler/scipy_doc/generated/scipy.spatial.Voronoi.html#scipy.spatial.Voronoi介绍完voronoi类后，我们就可以愉快地使用它了，在使用之前，首先需要生成voronoi多面体的种子点，通常种子点是空间上随机分布的点，不过如果我们对种子点施加一定的控制，可以生成很多有趣的空间形状，比如蜂窝状的空间结构。生成种子点后，将种子点坐标代入类中，就可以建立一个voronoi类了。随后就可以把这个voronoi类的信息转化为几何信息，利用abaqus的建模功能，生成voronoi多面体实体（也就是软件脚本语言中的cell对象）。Abaqus提供了从点到线到面到体的生成过程，生成voronoi多面体的过程就很好地体现了这一点：1.通过ridge_vertices属性，可以获知一个面上的顶点序号（贴心的是，这个属性提供的点是按照点的相连顺序排布的，而不是乱序，这一点非常重要），而通过顶点序号，就能从vertices属性中获知该点的坐标。2.由此就可以生成voronoi某个多面体中的一个面的各条边：使用脚本命令WirePolyLine将该面的各点顺次连接生成一个空间上的多边形线框。3.然后使用geometryedit的coveredges功能将线框围住的区域生成1个shell，这个shell就是多面体中的一个面了。4.对于该多面体中的其他表面，也可以如法炮制，直到生成这个多面体的全部外表面。然后再使用creatsolid—fromshell功能（对应的脚本方法为AddCells），即可生成该多面体。这就是点—线—面—体生成voronoi多面体的过程。下面我用一个配图例子来说明一下以上过程：首先，在空间内建立一个面上的顶点—边线关系，为了尽可能简单，我们采用三角形线框来构成一个面，生成三角形线框步骤如下：1.使用creatwire功能，输入各个点的坐标，首尾相连（选择chainedwires）形成一个三角形，点击ok即可生成线框：图2由点生成线框2.再将另外三条棱边输入，获得四面体的各条棱边：图3补全四面体的全部棱边3.使用Geometryedit工具中的CoverEdges功能，依次封闭各个线框：图4封闭线框为壳体封闭后的壳体如下图所示，为了表示清楚，这里做了一个剖面：图5壳体示意图4.使用createsolid：fromshell功能，将壳体填充成实体：图6由壳体生成实体填充后的实体如下图所示，为了显示方便，下图依旧使用了一个剖切面：图7实体示意图至此，由点坐标信息生成实体的过程就结束了，相应的Python脚本请见附件，将该脚本复制后粘贴至下图的命令行接口中，即可获得图示的四面体模型。利用这种思路，我们可以在空间中生成任意的凸多面体。无论一个Voronoi多面体有多么复杂，利用脚本就可以方便地生成它。在生成过程中，可以利用脚本，将不同的几何元素添加到set之中，方便后续赋予材料属性，添加约束或者接触，撒布种子点等等功能。通过使用set，从理论上说，使用者可以对每一个点，每一条边，每一个面和体进行控制，这是neper做不到的。值得注意的是，voronoi类的属性提供了三种类型的顶点：1)整个voronoi多面体的全部顶点vertices；2)各个面上的顶点ridge_vertices；3)各个体上的顶点regions。并没有提供哪些面是属于哪一个体的，这一点需要用脚本对坐标序号进行筛选判断，把面和体匹配起来，才能实现上述功能。最后附上成品效果图和脚本，CAE附件：图8Voronoi示意图，图中两种不同的颜色表示两种不同的材料图9对图8中的Voronoi多面体进行网格划分fromabaqusimport*fromabaqusConstantsimport*session.Viewport(name='Viewport:1',origin=(0.0,0.0),width=282.98975533247,height=152.084639698267)session.viewports['Viewport:1'].makeCurrent()session.viewports['Viewport:1'].maximize()fromcaeModulesimport*fromdriverUtilsimportexecuteOnCaeStartupexecuteOnCaeStartup()Mdb()session.viewports['Viewport:1'].setValues(displayedObject=None)p=mdb.models['Model-1'].Part(name='Part-1',dimensionality=THREE_D,type=DEFORMABLE_BODY)p=mdb.models['Model-1'].parts['Part-1']session.viewports['Viewport:1'].setValues(displayedObject=p)p.WirePolyLine(points=(((0.0,0.0,0.0),(5.0,0.0,0.0)),((0.0,5.0,0.0),(0.0,0.0,0.0)),((5.0,0.0,0.0),(0.0,5.0,0.0))),mergeWire=OFF,meshable=ON)e=p.edgesedges=e.getSequenceFromMask(mask=('[#7]',),)p.Set(edges=edges,name='Wire-1-Set-1')p.WirePolyLine(points=(((0.0,0.0,0.0),(0.0,0.0,5.0)),),mergeWire=OFF,meshable=ON)e=p.edgesedges=e.getSequenceFromMask(mask=('[#1]',),)p.Set(edges=edges,name='Wire-2-Set-1')v=p.verticesp.WirePolyLine(points=((v[0],v[2]),(v[3],v[0])),mergeWire=OFF,meshable=ON)p=mdb.models['Model-1'].parts['Part-1']e=p.edgesedges=e.getSequenceFromMask(mask=('[#3]',),)p.Set(edges=edges,name='Wire-3-Set-1')e=p.edgesp.CoverEdges(edgeList=e[0:1]+e[3:5],tryAnalytical=True)e1=p.