混凝土随机颗粒模型的网格自动剖分方法

第39卷第3期1999年5月大连理工大学学报Journa l of Da l i an Un iversity of Technology Vol.39,No .3M ay 1999文章编号:100028608(1999)03206收稿日期:1997209218;　修订日期:1998211210作者简介:王宝庭(1965～),男,博士,现在河海大学进行博士后研究;宋玉普(1944～),男,教授,博士生导师;赵国藩(1925～),男,教授,博士生导师,中国工程院院士.混凝土随机颗粒模型的网格自动剖分方法王宝庭,　宋玉普,　赵国藩(大连理工大学海岸与近海工程国家重点实验室,辽宁大连　116024)摘要:基于计算几何学原理,用一容易被广大工程技术人员接受的方法,实现了对混凝土随机颗粒模型三角形有限元网格的自动剖分;这个网格可以被各种数值方法所采用.同时自动形成V o rono i 图的方法,也可在工程中广泛利用.关键词:网格;三角剖分 随机颗粒模型;V o rono i 图中图分类号:T P391;TU 755文献标识码:A混凝土随机颗粒模型中,在随机地生成颗粒分布之后,对结构进行有限元或离散元分析时,首先要对研究区域进行网格剖分.这将涉及到成千上万个单元和节点,数据准备工作则变得极其冗长而乏味;而且在数据的准备过程中,难免留下一些人为错误,尽管经常核对,但仍无法察觉.特别是对于需要随机改变颗粒组合、反复进行网格剖分的情况下,不去寻求某种自动或半自动的方法,是难以实现的.本文基于计算几何学上V o rono i 图和D elaunay 图的性能,通过D elaunay 三角剖分实现了对随机颗粒模型进行全自动剖分.1　凸壳问题一个多边形是凸多边形,是指多边形内任意两点的连线全部包含在该多边形内.平面上的n 个点P 1,P 2,…,P n 的凸壳是包含这些点的最小凸多边形(为简单起见,假设没有三点共线),如图1所示.凸壳可以用它的顶点的逆时针表列表示,求其表列是最基本的凸壳问题.凸壳的各边是这样一种线段,对于通过P i 、P j 的直线来说,其他的点位于直线的同一边,那么,连接P i 、P j 的线段就是凸壳的一条边.以上判断可通过如下计算实现:x i y i 1x jy j 1x k y k 1>0(1)若上式成立,则表明P i 、P j 、P k 成顺时针排列,即P k 点在线段P i P j 的左边.对于n 个点来说,需要对n (n -1) 2个点对P i ,P j 判断上述条件是否成立来求出凸壳.对于一个点对的判断需要O (n )次判断,因此,用此方法求出凸壳需要的时间为O (n 3). 这个朴素的方法确实简单,但却过于花费时间,因此需寻找某种算法,以期能够在普通微机上实现.以下采用分治法,凸壳问题可用O (n log n )的时间求出.分治法的思考方法是:把问题分成大体相等的两部分,递归求解各部分问题,组合两部分的解得到全体解.解凸壳问题时,问题的分解就是点集的分割.可以把点集任意分成两部分,但为方便,以点的x 坐标值为界将点集左右分割.设有n 点的集合S ,以x 坐标的顺序分割成左右两部分S L 和S R ,分别对于S L 和S R ,递归求出它们的凸壳CH (S L )和CH (S R ).由于把S 关于坐标x 值分割,因而CH (S L )和CH (S R )不重叠,左右分开;如图2所示.以下要从两个凸壳CH (S L )、CH (S R )求出全体点集S 的凸壳CH (S ).为此,由于CH (S L )和CH (S R )不重叠,可以求出图2所示的CH (S L )和CH (S R )的公切线L -、L +,然后除去CH (S )不需要的部分(图中用虚线表示的部分),就可得到CH (S ). 图1　凸壳 图2　CH (S L )和CH (S R )的合并图3　求L -的过程现在要解决L 的求法.凸壳CH (S L )和CH (S R )的端点中y 坐标值最小的点分别为P L 、P R ,假设P L 在P R 的下面,这时L 的求法如图3所示,可由如下步骤求出:(1)判断直线P L P R 是否与CH (S R )相切(即判断CH (S R )的角点是否在直线P L P R 同一边),若不是,则沿CH (S R )的顶点逆时针顺序更新P R ,直到相切为止.(2)判断直线P L P R 是否与CH (S L )相切,若不是,则沿CH (S L )的顶点逆时针顺序更新P L ,直到相切为止.(3)线段P L P R 即为要求的L -.此方法用O (n )时间求出线段L -,线段L +也同样用O (n )的时间求出.求出L -、L +之后,为了求得CH (S ),需将CH (S L )和CH (S R )不需要的部分除去,只要沿凸壳角点搜索即可.总之,从CH (S L )、CH (S R )求出CH (S )需要O (n ).于是,由于分治法大体是2等分,因而,n 点的凸壳需要O (n log n )求出.2　n 个半平面的交现在考虑,给出n 个半平面,求其相交的领域(图4).仍用分治法,把n 个半平面的集合S 分成大体相等的两部分S 1、S 2,对S 1、S 2分别递归求出其相交的领域R 1、R 2.显然R 1和R 2的交是要求的领域,问题是如何求出R 1和R 2的交.一般地,求两个凸多边形的交需要花费与全体点数n 成正比的时间.例如,用通过各点的垂线把平面分成N +1个带状区域,在每个带状区域内求两个梯形(或三角形)的交,将它们644大连理工大学学报第39卷　合并即得,如图5所示.R 1、R 2是半平面的交,因此是凸多边形(有可能是无限的).用上述方法求R 1和R 2的交要用O (n )的时间,因而,半平面的交要用O (n log n )的时间. 图4　n 个半平面的交 图5　求两个凸多边形的交图6　28个点的V o rono i 图3　V o rono i 图V o rono i 图是计算几何学的主要研究领域.给出平面上n 个点P 1,P 2,…,P n ,点P 的V o rono i 区域V (P i )可表示如下:V (P i )=∩j {P d (P ,P i )≤d (P ,P j )}(2)d (P ,P i )表示点P 和P i 之间的欧几里得距离.V (P i )表示到P i 的距离比到其他任何点的距离都近的点的集合,即点P i 的“势力圈”.V (P i )分割平面而形成的图称为V o rono i 图.图6是一示意图.仍用分治法去求V o rono i 图.首先将n 点的集合S 按x 坐标的顺序大体分成左右两部分,即S L 、S R .用求n 个半平面的交的方法,递归求出S L 、S R 的V o rono i 图V (S L )、V (S R ).由于把S 关于x 的坐标左右两分,可以按以下步骤求得全体点集的V o rono i 图V (S ).(1)求满足m in{d (P ,P i ) P i ∈S L }=m in{d (P ,P j ) P j ∈S R }的点P 的轨迹,这个轨迹如图7(a )所示,是关于y 坐标的一条单调连续折线L .(2)把折线L 右边的V (S L )部分除去,折线L 左边的V (S R )部分除去,就得到V (S ),如图7(b ). 图7　V (S L )和V (S R )的合并图8　1000个随机点的V o rono i 图744　第3期 王宝庭等:混凝土随机颗粒模型的网格自动剖分方法用以上方法求得V o rono i 图需要的时间是O (n log n ).图8为1000个点的V o rono i 图.4　D elaunay 三角剖分的计算混凝土区域的剖分可以通过D elaunay 三角剖分的概念得到.D elaunay 三角形的主要性能表述为:通过D elaunay 三角形的三个顶点的外接圆不包含另外的基本点(空圆盘性质).这将产生一个很方便的计算工具,即“圆内测试”(图9).它可以通过以下判断实现:F (A ,B ,C ,D )=x Ay A x 2A +y 2A 1x By B x 2B +y 2B 1x Cy C x 2C +y 2C 1x D y D x 2D +y 2D 1>0(3)仍用分治法,几何剖分n >5的点集为两部分S L 和S R ,递归地计算S L 和S R 的D elaunay 的三角剖分,最后连接S L 和S R 的三角形剖分.如基本点数n ≤5,可以基于D elaunay 三角形剖分的定义,使用“圆内测试”直接计算.以n =5作为将来进一步剖分的最小限的原因是显然的,因为如果5个点被分成两半的话,至少有一半它的点数少于3,比3个点少的D elaunay 三角形剖分在数学上是行得通的,但对于工程技术人员来说是不直观的. 图9　圆内测试图10　S L 和S R 的连接这两个凸壳CH (S L )和CH (S R )的连接可以通过一个冒泡算法和转换线段来实现.由D elaunay 三角形的性质可知,没有D elaunay 三角形的边在基本点的凸壳外面和所有凸壳的边都是D elaunay 三角形的边,因此,由前述方法找到的线段L-即是L R 边(连接CH (S L )和CH (S R )的边),图10显示其底边和顶边.现在,认为L cand 是G 和位于凸壳CH (S L )上但不属于凸壳CH (S L )∪CH (S R )的基本点,换句话说,是凸壳消影算法中应消去的点,R cand 是H 和位于凸壳CH (S R )上但不属于凸壳CH (S L )∪CH (S R )的基本点.现在想象一个胀开圆从底边开始向上冒泡,如图11.这种接连不断的圆将以一个固定的次序遇到L cand 和R cand 的所有点.考虑“圆内测试”和D elaunay 三角形的定义,去发现各个L R 边.如图12(a )所示,设A D 和CD 是被冒泡算法发现的L R 边,显然通过A 、B 、C 的外接圆包含D (即∠A B C +∠A D C >∠B A D +∠B CD ),因此,需要转换边A C 为B D (图12(b )),形成两个新的三角形A B D 和B CD .对于三角形A EB 和D A B 也必须去履行这种转换,直到不再需要新的转换.这个转换绝对不增加新的L 边(在S L 内的边)或R 边(在S R 内的边).在x 或844大连理工大学学报第39卷　 图11　冒泡过程图12　线段转换y 的方向都可根据基本点的坐标划分成L 和R ,但不同的划分需要转换线段的数量将大不相同.可以证明,用冒泡程序与转换线段一起,符合计算几何学原理,能产生一个精确的D elaunay 三角形剖分.5　基本点集的产生及网格剖分本节以随机颗粒模型为例,自动产生基本点集.首先根据级配曲线,计算各级别骨料颗粒在试件内的个数,然后从最大粒径骨料颗粒开始,随机地产生每一颗粒的形心在试件内的位置,已定位的颗粒所占位置及周围一定范围之内,其他颗粒不能侵入,按照这一规则,从大到小随机定位,直至最后一个颗粒进入截面为止.本文将全级配混凝土的骨料分成四种不同粒径,分别为12.0、6.0、3.0、1.5c m ,试件截面尺寸为45c m ×45c m ,按上述方法产生的一个截面,如图13(a )所示. (a )随机颗粒图(b )对应图(a )的1797个节点、3406个单元剖分图图13　随机颗粒图试件截面形成以后,将1.5、3.0、6.0、12.0c m 四种不同粒径的骨料颗粒按一定分布规律分别取不同个节点,硬化水泥浆基体按一定的规则自动生成节点的坐标,然后用前述方法连944　第3期 王宝庭等:混凝土随机颗粒模型的网格自动剖分方法054大连理工大学学报第39卷　接这些点,就可得到如图13(b)所示的一网格剖分图.6　结束语本文基于计算几何学的原理,用一容易被广大工程技术人员接受的方法,实现了对随机颗粒模型三角形有限元网格的全自动剖分.其V o rono i图也有广泛的应用:可用来模拟晶体中的晶粒,以从细观上研究晶体的力学行为;可以用来模拟堆石坝的截面;也可以用来模拟节理岩体.此方法具有以下特点:(1)这个方法基于简单的原理,具体的算法都容易实行.(2)在普通微机上就能得到快速和好的结果.(3)这些网格可以为各种数值方法所采用.参考文献[1]普雷帕拉塔F P,沙莫斯M I.计算几何导论[M].北京:科学出版社,1990.[2]GU I BA S L,STOL F I J.P ri m itives fo r the m ani pu lati on of general subdivisi ons and the compu tati on ofV o rono i diagram s[J].AC M Tran s Graph ics,1985,4(2):742123.[3]W IT TM ANN F H,RO EL FSTRA P E,S ADOU K I H.Si m ulati on and analysis of compo sitestructu res[J].M ater Sc i Eng,198421985,68:2392248.[4]SCHU T T ER G D,TA ERW E L.R andom particle model fo r concrete based on D elaunaytriangu lati on[J].M ater and Struct,1993,26:67273.Automa tica lly d iv id i ng way of m esh forrandom particle m odel of concreteW ANG　B a o2ting,　SONG　Yu2pu,　ZHAO　G uo2fa n(S ta te Ke y La b.of C oa s ta l&O ffsho re Eng.,D a lia n Univ.of Te chno l.,D a lia n116024,C hina)Abstract:A m ethod based on com pu tati onal geom etry,easily accep ted by techn ician s is p resen ted to au tom atically generate the triangle fin ite elem en t m esh fo r the random particle m odel of concrete.T he generated fin ite elem en t m esh is also app licab le to a variety of o ther num erical m ethods.In additi on,the p resen ted m ethod can generate the V o rono i diagram, w h ich has w ide u se in eng ineering.Key words:m esh;triangu lati on random particle m odel;V o rono i diagram。