基于图像八叉树的三维比例边界有限元多面体网格生成算法

章 鹏,杜成斌,赵文虎

(1.南京工程学院土木工程与智慧管理研究所,江苏 南京 211167;2.河海大学力学与材料学院,江苏 南京 211100; 3.南昌大学工程建设学院,江西 南昌 330031)

有限元作为主流的固体结构数值分析方法已广泛应用于各种工程问题中[1-3]。在复杂的三维几何结构中,标准有限元网格一般由几种简单的几何单元组成,包括四面体单元、五面体单元及六面体单元[4-5]。四面体网格相较于六面体网格容易生成,然而四面体作为常应变单元,分析计算的精度相对较差。由于有限元法中单元网格边界的协调性要求,对于复杂几何形状(边界)的结构,建立高质量的六面体网格极其困难,在边界处会出现与实际边界形状不符的“台阶”或“锯齿”形状,有时需要大量的人工干预。此外,针对移动边界和裂纹扩展等问题,网格剖分甚至无能为力[6]。

近些年来学者们提出了很多数值方法以解决网格剖分及重剖分的问题,如扩展有限元法[7](extended finite element method,XFEM)、等几何分析法[8](iso-geometric analysis,IGA)、无网格方法[9](mesh free method,MMs)、多面体有限元法[6](polyhedral finite element method,PFEM)等。扩展有限元法通过在有限元位移场中引入表征位移不连续的富集函数,将不连续问题嵌入到有限元中,在裂纹扩展时,不需要进行网格重剖分工作,但在复杂三维问题中,空间裂纹面描述困难、裂尖单元构造复杂给扩展有限元法的三维结构开裂模拟带来了困难[10]。由Hughes等[11]提出的等几何分析法,基于有限元分析方法的等参单元思想,将计算机辅助几何设计(computer aided geometric design,CAGD)中用于表达几何模型的非均匀有理B样条(non uniform rational B-spline,NURBS)的基函数作为形函数,实现了计算机辅助设计(CAD)和计算机辅助工程(CAE)的无缝结合。然而传统CAD曲面中的NURBS曲面缺乏水密性,大量的裂缝、重叠等几何瑕疵使得NURBS表示很难用于等几何分析[8]。无网格法根据一些任意分布的坐标点构造插值函数离散控制方程,在数值计算中不需要生成网格,具有很强的灵活性,然而由于近似函数一般均很复杂,计算量较大,计算效率缺少优势,且边界条件的施加较为烦琐[12]。多面体有限元法对于复杂几何结构的网格剖分更加灵活,同时具有更高的数值精度,该方法的障碍是难以构造满足单元协调性的多项式位移插值形函数[13]。

Song等[14]在研究无限域动力相互作用时提出了比例边界有限元法(scaled boundary finite element method,SBFEM),该方法结合了有限元法(finite element method,FEM)和边界元法(boundary element method,BEM)的优点。与边界元法相比,SBFEM不需要基本解,只需对边界进行离散,降低了数值模拟的维度;它在径向是解析的,而在环向的离散采用与有限元同样的插值方法,具有半解析的特性,相比常规有限元法,具有更高的计算精度[15-16]。基于上述优点,SBFEM在二维问题中得到广泛应用[17-26]。相对来讲,三维比例边界有限元的研究发展较为缓慢,主要原因之一为三维复杂结构的离散难以实现,直到近年来才取得一些进展。Saputra等[27]根据X-射线扫描图像,采用平衡八叉树自动网格划分三维混凝土数字图像,然而平衡八叉树方法在离散边界为曲线或者倾斜的结构时,会出现“台阶”或“锯齿”形结构边界,导致网格难以准确反映真实结构边界。Liu等[28]针对三维STL(standard tessellation language)格式的几何模型提出了一种多面体网格划分方法,通过裁剪结构边界形成多面体SBFEM网格,然而该方法仅针对三维STL格式,对于三维图片等格式不具有通用性。孔宪京等[29]提出了用Voronoi多面体方法离散三维结构,用多边形插值函数代替四边形插值函数,求解环向多边形边界面,但对复杂边界的裁剪方法未作深入介绍。

本文在八叉树-多面体网格自动剖分算法基础上,基于Matlab软件,自主研制了基于图像八叉树三维SBFEM多面体网格的算法程序,并通过数值算例,对本文算法生成网格的SBFEM计算结果的精度和边界适应性进行验证。

1 三维比例边界有限元法简介

在比例边界坐标系(ξ,η,ζ)中,ξ为径向坐标,η、ζ为环向坐标。在比例边界有限元法中需要设置比例中心,比例中心对子域内任意一点可见,在边界上进行离散化,在多面体上的每个面都要离散成三角形和四边形,如图1所示,子域内任意一点b,该点笛卡尔坐标系的坐标(x,y,z),可由ξ、η和ζ表示[30]:

x(ξ,η,ζ)=ξN(η,ζ)xb

(1)

y(ξ,η,ζ)=ξN(η,ζ)yb

(2)

z(ξ,η,ζ)=ξN(η,ζ)zb

(3)

式中:N(η,ζ)为形函数;xb、yb、zb为边界上结点坐标集。

设位移解形式为

u(ξ,η,ζ)=N(η,ζ)u(ξ)

(4)

式中u(ξ)为径向上各点位移。

单元中任意一点的应变[30]、应力分别为

ε(ξ,η,ζ)=B1(η,ζ)u(ξ),ξ+B2(η,ζ)u(ξ)/ξ

(5)

σ(ξ,η,ζ)=D(B1(η,ζ)u(ξ),ξ+B2(η,ζ)u(ξ)/ξ)

(6)

式中:B1(η,ζ)、B2(η,ζ)为单元应变矩阵;D为材料的弹性矩阵。

由域Ω中的平衡方程可以推导得到位移解[31]:

u(ξ)=φ11ξλ-0.5Ic

(7)

式中:φ11为位移模态向量;λ为位移模态对应的特征值向量;c为位移模态常数。

将式(7)代入到式(4)、式(6)即可得出单元中任意一点的位移及应力。三维比例边界有限元法可以方便地处理边界上的悬挂节点,即仅需在单元表面中心添加节点并连接悬挂节点,使单元表面由三角形及四边形组成[32],所以可以很好地与八叉树网格剖分相结合。

2 基于图像八叉树SBFEM多面体裁剪技术和网格生成算法

本文的网格生成主要包括八叉树网格自动生成和多面体裁剪两部分,后者主要针对不同材料的过渡区和复杂的几何边界,为适应边界几何形状,对生成的规则单元进行多面体裁剪。同时,为适应模拟尺度的需要,结合最新数字图像技术的发展,本文针对图像技术的八叉树网格生成和多面体裁剪技术进行研究。

2.1 图像八叉树数据结构及生成技术

2.1.1 结构几何图像像素的确定

首先根据求解精度要求,将结构图像的分辨率设置为h×h×h像素,形成涵盖整个结构图像信息的三维矩阵[33]。矩阵中每个元素均代表大小为1×1×1的立方体像素,其颜色取值范围为0～255,不同的取值代表不同的颜色。在图像八叉树中,h需满足h=2n(n为正整数),还需要设置最大八叉树单元尺寸Smax,该最大单元尺寸为边界框长度的1/(2m)(m正整数),此外需要设置最小单元尺寸Smin,Smin=Smax/Sratio,其中Sratio为最大、最小单元尺寸比值。

本文采用Matlab开源inpolyhedron函数进行像素单元划分,如果像素单元所有的节点皆在结构体内部则为结构内像素,如果像素单元所有的节点皆在结构体外部则为结构外像素,如果像素单元有部分节点在结构体内部,有部分节点在结构体外部则为结构边界像素。将3种不同单元像素分别进行赋值,本文将结构内像素颜色设置为1,将结构边界像素颜色设置为100,将结构外单元像素颜色设置为255。

2.1.2 图像平衡八叉树剖分

图像八叉树的剖分原理为:首先根据Smax将根单元均匀分割为(h/Smax)×(h/Smax)×(h/Smax)个小立方体,生成初始八叉树单元网格,如果一个单元格中最大和最小像素颜色间的差值大于预设的颜色阈值,则该单元格进一步分为8个相等大小的单元格。以图2所示的二维图片为例,对包含结构单元信息的minx≤xp≤maxx像素图片(图2(a))进行八叉树二维对应的四叉树分解,其中颜色阈值取为0,Smin=1、Smax=4,分解结果如图2(b)所示。

图2 二维图像剖分示意图Fig.2 A schematic diagram of 2D image decomposition

图像八叉树网格剖分完成之后,需要采用平衡八叉树(2∶1规则)剖分,以保证相邻单元边长比例不超过2∶1,最后将结构外单元(像素为255的单元)进行删除。图3为一个1/8圆球八叉树分解的过程,其中初始图片大小为128×128×128像素,Smax=32,Smin=4。

图3 1/8圆球八叉树分解过程示意图Fig.3 A schematic diagram of octree decomposition process of 1/8 sphere

2.2 图像结构八叉树网格的多面体裁剪

平衡八叉树网格剖分完成之后,在边界单元中由于八叉树网格只有水平及竖直2种边界,在边界上会存在着“锯齿”或“台阶”现象,因此对八叉树网格的边界单元进行裁剪非常必要。

2.2.1 图像八叉树网格与结构面的交点计算

边界单元与结构面相交必然经过边界单元的棱边,因此计算边界单元棱的交点可得到图像八叉树网格与结构面的交点。

如图4所示,由点La(xa,ya,za)及点Lb(xb,yb,zb)形成的八叉树边界单元棱边与由点F0(0,0,0),F1(x1,y1,z1)及点F2(x2,y2,z2)形成的结构面相交于点P(xp,yp,zp)。

图4 某一棱边与三角形相交示意图Fig.4 Intersection diagram of an edge and a triangle

因为点P在单元棱边上,因此可表示为

(8)

同时点P在结构面上,因此也满足:

(9)

式中t、g、v为待求解参数。

因为每一个八叉树的棱边都平行于一个坐标轴,因此(xb-xa,yb-ya,zb-za)中只有一个非零值,本文以平行于x轴的棱边为例:

(10)

由此,可分别求解得到t、g、v的值,将t代入到式(8),即可得到点P在棱边的位置。采用旋转方法,即可得到平行于y轴及z轴的棱边的交点位置。

2.2.2 图像八叉树网格单元与可能相交的结构面的筛选

一个复杂的几何结构,通常有上千个几何结构表面个数,如果直接采用式(10)对每个结构面与每个边界单元棱边进行交点计算,计算量很大。因此在进行交点计算前,提出采用面-立方体相交判断方法进行边界单元与结构边界相交面的筛选。

设定边界单元x、y和z轴坐标区间分别为[minx,maxx]、[miny,maxy]、[minz,maxz],结构面坐标区间分别为[x1′,x2′]、[y1′,y2′]、[z1′,z2′],则结构面若与边界单元相交,假定相交坐标点为(xp,yp,zp)。

由于该点在边界单元和结构面上,则在x轴必满足:

minx≤xp≤maxxx1′≤xp≤x2′

(11)

由上式可知,x1′、x2′分别满足:

x1′≤maxxx2′≥minx

(12)

同理,y轴、z轴坐标同时满足

y1′≤maxyy2′≥miny

(13)

z1′≤maxzz2′≥minz

(14)

因此,结构面若与边界单元相交,其坐标区间则必须同时满足式(12)～(14)。结构面坐标区间与边界单元坐标关系为面-单元体可能相交的筛选判别条件。

2.2.3 图像八叉树网格切割处理

2.2.3.1 图像八叉树边界单元表面裁剪

裁剪边界单元表面的关键是裁剪表面的棱边,首先确定单元表面节点与结构体的关系,如果在结构外部(用“-”表示),节点需要删除,如果在结构内部(用“+”表示),保留节点,然后连接边界单元与结构体的相交点,完成单元表面的裁剪,如图5所示,分别为三角形、四边形、五边形表面的裁剪过程。

图5 边界单元表面裁剪过程Fig.5 Surface trimming process of boundary element

2.2.3.2 图像八叉树边界单元切割面裁剪

首先确定边界单元中与结构体相交的所有点,如图6所示,a、b、c、d为边界单元的4个交点,将交点以外的单元节点裁剪掉,然后将交点进行连接,为了防止交点顺序混乱,需要按照一定顺序进行连接,如图6(b)所示a、b、c、d以逆时针顺序进行连接,形成表面单元切割面。由相交点位置和个数决定切割面的形态,如3个相交点切割面为三角形,4个相交点为四边形,5个相交点为五边形。

图6 边界单元切割面裁剪及表面三角化过程Fig.6 Cutting surface trimming and triangulation of boundary element

边界单元多面体修剪完成后,多面体网格表面需采用前文所述的三角形及四边形2种表面单元形函数,将多于4个节点的单元表面进行三角化,即在该多边形的形心处插入一个节点,并通过每个节点与形心相连创建三角形。此外由四点组成的切割面可能不在同一个平面内,因此本文统一将切割面进行三角化。如图6(c)所示,为三角化后的多面体网格表面。

2.3 图像八叉树-多面体网格生成流程伪代码实现

为了更好地描述本文基于图像八叉树-多面体网格裁剪及生成的过程,给出流程伪代码如下:

步骤1输入结构几何信息,根据结构尺寸,输入满足精度的像素大小。

步骤2for每一个像素do:

根据像素立方体8个端点是否在结构体内部进行像素单元性质划分,并赋值像素信息。

end for

步骤3根据像素信息,建立图像平衡八叉树单元。

步骤4for 每一个边界单元中do:

①筛选边界单元与结构边界相交面; ②进行边界单元棱边与结构面交点计算,并确定边界单元节点与结构面关系;③删除边界单元的结构外节点,进行图像八叉树边界单元表面裁剪;④通过相交点的逆时针排序连接,形成图像八叉树边界单元切割面;⑤八叉树边界单元表面统一变换为三角形及四边形两种单元形式。

end for

步骤5由结构内单元以及修剪的边界单元共同组成结构多面体单元网格。

其中步骤1、步骤2对应2.1.1节,步骤3对应2.1.2节,步骤4对应2.2节。

通过该算法生成的八叉树网格和修剪后的多面体网格可直接与三维比例有限元计算软件无缝对接,作为输入的基本数据,进而快速地计算出三维结构的位移和应力。

3 数值算例验证

3.1 受压空心球体

一空心球体如图7所示,内径r1=10mm,外径r2=50mm,球体内侧,受到一均布外向压力P=10GPa。球体的材料参数:弹性模量E=200GPa,泊松比ν=0.3。该问题关于球体坐标(r,θ,φ)的位移具有解析解[34]。考虑到几何结构及荷载的对称性,本文采用1/8球体进行模拟,根据位移解析解,在球体的内表面以及外表面施加位移边界条件。

图7 受压空心球体Fig.7 Compressed hollow sphere

采用八叉树-多面体自动网格剖分算法对该算例进行网格剖分,其中图像像素设定为128×128×128,Smin=r2/16,Sratio=8,图8给出了多面体网格。

图8 受压空心球体网格生成Fig.8 Mesh generation of compressed hollow sphere

为了对不同网格的算例进行分析,分别取Smin=r2/16、r2/32、r2/64、r2/128共4种网格进行计算模拟,为便于比较,取Sratio为定值8,4种网格对应的节点数分别为4094、17150、68987和278154。

采用位移相对误差进行精度比较,4种网格相对误差分别为0.0176、0.0051、0.0032和0.0019,可知随着网格密度的增加,相对误差逐渐减小。

3.2 Cook’s membrane 问题

为了检验三维网格剖分的普适性,本文选择一个经典非规则的三维Cook’s membrane平板问题,该算例被多位学者用于数值验证[35]。几何尺寸如图9所示,其中厚度为5mm,该结构的材料参数:弹性模量E=1000MPa,泊松比ν=0.33,结构的右侧边缘作用τ=10MPa的竖向均布荷载。

图9 Cook’s membrane问题几何尺寸(单位:mm)Fig.9 Geometry size of Cook’s membrane problem (unit: mm)

采用八叉树-多面体自动网格剖分算法进行网格剖分,其中图像设定为像素128×128×128,分别对Smin=b/16、b/32、b/64、b/128(b为最大边长,b=60mm)的4种网格进行计算模拟,其中Sratio=8,4种网格节点个数分别为642、2637、10366和41163。图10为该板其中3种不同密度的图像八叉树网格。

图10 Cook’s membrane问题不同密度的图像八叉树网格Fig.10 Image octree mesh of Cook’s membrane problem with different densities

图11为不同网格下沿着AC方向各点的y方向位移,可以看出,不同网格下的位移结果皆较为相近。图12为不同网格节点时节点C的y向位移,图中同时给出了Ooi等[35]提出的双重比例边界有限元结果。由图12可知,随着网格密度的增加,位移表现出快速收敛的趋势,当单元节点数为10336,位移为3.9822mm,随着网格进一步加密,结果基本保持稳定,与Ooi等[35]结果较为吻合,且收敛速度略好于后者。

图11 沿直线AC方向的各点y向位移Fig.11 y-direction displacement of points along the AC direction

图12 不同网格节点时节点C的 y 向位移Fig.12 y-direction displacement of node C with different element nodes

4 结 语

基于图像八叉树方法,提出了针对三维结构的平衡八叉树和多面体网格修剪算法,基于Matlab软件,编制了相应的网格算法程序。本文算法具有以下优势:①基于图像像素方法进行平衡八叉树-多面体的网格剖分与裁剪,仅需提前设置图像像素大小及单元尺寸等参数,整个网格剖分过程不需要任何人工干预,提高了网格的剖分效率。②生成的网格单元可以与比例边界有限元计算软件无缝对接,由于SBFEM可方便处理悬挂节点问题,因而可方便实现粗细网格的快速过渡,在保证精度情况下,相较于有限元法,自由度数可大大减少。③对于复杂的结构边界问题,该算法可以较准确得到与实际结构边界接近的数值模型。因此该算法在三维结构数值模拟问题中,具有较大的优势及潜力,为工程结构的全自动化分析提供了一个新的途径。