基于圖像八叉樹的三維比例邊界有限元多面體網格生成算法

章 鵬,杜成斌,趙文虎

(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]?；谏鲜鰞烖c,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可方便處理懸掛節點問題,因而可方便實現粗細網格的快速過渡,在保證精度情況下,相較于有限元法,自由度數可大大減少。③對于復雜的結構邊界問題,該算法可以較準確得到與實際結構邊界接近的數值模型。因此該算法在三維結構數值模擬問題中,具有較大的優勢及潛力,為工程結構的全自動化分析提供了一個新的途徑。