基于Voronoi結構的多邊形單元網格生成方法

石怡婧，邵國建，丁勝勇

(河海大學 力學與材料學院，江蘇 南京 210098)

?

基于Voronoi結構的多邊形單元網格生成方法

石怡婧，邵國建，丁勝勇

(河海大學 力學與材料學院，江蘇 南京 210098)

摘要：針對多邊形單元網格難以生成的問題，建立了基于Voronoi結構的多邊形單元網格生成方法。該方法通過區域內的一組初始點來構造Voronoi結構。對初始種子點的分布進行優化，以達到控制多邊形單元網格密度的目的。利用Voronoi結構的特點，通過添加種子點對應目標區域邊界的映射點對邊界進行擬合，實現Voronoi結構對復雜邊界的逼近。對Voronoi結構進行質心化，改善生成多邊形單元形態。給出了該方法的程序實施步驟，并結合網格生成實例，驗證了其合理性和可行性。

關鍵詞：有限元法；多邊形單元；網格生成；Voronoi結構

0引言

在有限元的分析過程中，由于研究對象幾何結構和材質的復雜性，勢必會增加建立有限元網格模型的難度。尤其對于形狀比較復雜的幾何形狀，其網格剖分耗時較多，并且容易出錯。因此，網格剖分技術是網格模型建立的核心內容，一個好的有限元網格對提高計算效率和計算精度起著關鍵作用。目前，對于任意形狀的分析區域，三角形及四邊形單元的網格生成技術已經較為成熟。而多邊形單元作為一種新型的單元形式，由于其邊數及形狀的任意性，與常規單元相比有獨特的優勢。針對復雜幾何體進行網格劃分更加靈活、方便[1-2]，近年來受到越來越廣泛的關注[3]。文獻[4]提出了Delaunay多邊形化生成多邊形網格的技術。文獻[5]將一種類似于金字塔型的方法用于劃分復雜三維模型多邊形網格中。文獻[6]利用傳統三角形有限元網格生成形狀合理的多邊形單元網格。文獻[7]利用多邊形比例邊界有限元方法生成多邊形網格。但是，目前多邊形單元的網格生成算法還比較粗糙，生成的多邊形單元網格非常畸形，處理復雜邊界處的網格時，構造能較好滿足多邊形有限元分析的單元形狀也有一定困難。

本文以Voronoi結構為基礎，通過對Voronoi結構的質心化來彌補現有多邊形單元網格生成方法中單元形態較差的缺陷。同時，探討非均勻網格生成的解決方案及復雜邊界的擬合問題。給出多邊形單元網格自動生成程序的實施步驟，并結合具體實例來驗證本文方法的合理性和可行性。

1Voronoi結構及其質心化

Voronoi圖是針對平面離散點集而言，將平面分成若干個區域，各區域僅有一個點，該點所在的區域是距離該點最近點的集合[8]。Voronoi圖的數學定義：記R2空間上任意兩點xi和xj的歐氏距離為d(xi,xj)，P={P1,P2,…,Pn}是R2空間內互異離散點集合，設

(1)

則稱V(Pi)為任意點Pi(也可稱為種子)的Voronoi區域，又稱為Voronoi單胞。Voronoi單胞的邊是集合P中某兩個點連線的垂直平分線，相鄰Voronoi單胞的種子點連線構成Voronoi圖的對偶圖，被稱為Delaunay三角化。圖1為Voronoi圖及其Delaunay三角化。如圖1所示，“+”號表示種子點，實線表示Voronoi圖，虛線表示Delaunay三角化。在多邊形有限元法中，可以將Voronoi單胞視為多邊形單元。

圖1　Voronoi圖及其Delaunay三角化

一般在Voronoi圖中，Voronoi單胞的種子點與其質心并不重合，此時Voronoi單胞的形狀較為畸形。而質心化Voronoi圖(centroidalVoronoitessellation，CVT)是Voronoi圖的一種特殊形式[9-11]，其種子點與對應的Voronoi單胞的質心重合，即：

(2)

(3)

若該聚類能量函數ε(P)的值逐漸減小并收斂到一個固定值，即：

(4)

此時，生成的Voronoi圖為CVT。

Lloyd算法[10]是目前廣泛使用的計算CVT的一種算法，是一種必然收斂的迭代算法。由于CVT形式不唯一，所以Lloyd算法是局部收斂的。

圖2為基于Lloyd算法的Voronoi結構優化過程。圖2中，“+”號表示種子點，“⊙”號表示Voronoi單胞質心，在單位正方型區域內隨機分布20個種子點，然后利用Lloyd算法生成CVT，優化步數k=20。圖2a為初始的Voronoi結構，此時ε(P)=4.45×10-3，種子點和質心點位移差異較大，Voronoi單胞形狀很差。圖2b為k=2的Voronoi結構，此時ε(P)=2.74×10-3，種子點和質心點位移差異明顯減小，但Voronoi單胞形狀仍然比較差。圖2c為k=20的Voronoi結構，此時ε(P)=2.31×10-3，種子點和質心點基本重合，Voronoi單胞形狀良好。圖2的優化結果顯示：優化過程中ε(P)逐漸減小并趨于穩定，Voronoi單胞形狀由畸形逐漸優化成比較規則的多邊形。

圖2　基于Lloyd算法的Voronoi結構優化過程

2初始種子點分布優化

在有限元分析中，通常需要對一些重點關心區域采用更精細的網格，對一些非重點區域采用較為粗略的網格。一般初始種子點的生成都是采用隨機布點的方法，這種方法會使多邊形單元密度分布呈現隨機化，難以滿足有限元分析的要求。初始種子點的密度會在很大程度上影響最后CVT的密度，所以本文對初始種子點的分布進行優化。

本文初始種子點分布優化采用文獻[12]提出的三角形網格生成算法，其基本思想是將三角形網格類比為平面彈簧結構，三角形的三條邊分別對應連接種子點的各彈簧，每根彈簧的現有長度l和期望長度l0遵循力-位移關系f(l,l0)，并通過改變種子點的位置使整個結構達到平衡。圖3為帶圓孔的正方形區域內初始種子點優化前后的分布情況。通過圖3可以看出：隨著優化過程的進行，種子點分布越來越規律，步驟1到步驟2中種子點分布呈現出一定的規則性，步驟2到步驟3過程中則刪除了分布在邊界上的種子點。

圖3　帶圓孔的正方形區域內初始種子點的優化分布

3Voronoi結構邊界處理

因此，可通過下式得到種子點對應的映射點：

(5)

圖5為Voronoi結構的邊界處理圖。圖5中，“+”為未映射的種子點，“°”為映射點。對大多數種子點而言，對應映射點的存在不一定有助于最終Voronoi結構的生成，因此，需要對種子點進行篩選。將種子點Delaunay三角化，選取Delaunay三角形網格邊界上的種子點即最外層種子點作為待映射的種子點，并生成其對應的最近邊界處的映射點。由于上文種子點對應的Delaunay三角形網格中三角形都為銳角三角形，很容易證明只有最外層種子點與其對應的最近邊界處的映射點生成的Voronoi單胞存在公共邊界。

圖4 種子點及其對應的映射點圖5 Voronoi結構的邊界處理

Voronoi結構邊界生成具體步驟為：(Ⅰ)在目標區域生成種子點集Pin；(Ⅱ)篩選種子點，并進行種子點對應于目標區域邊界?Ω的映射點集Pout的生成，這些映射點位于分析區域的外部；(Ⅲ)進行點集P=Pin∪Pout對應的Voronoi結構的生成，種子點集Pin對應的Voronoi單胞集為目標區域的網格。

在步驟(Ⅲ)生成的Voronoi結構中，如果種子點和其對應的映射點有共同的Voronoi單胞邊界，則所有共同邊界的集合構成了目標區域邊界。以帶圓孔的正方形區域為例，步驟(Ⅲ)中的點集P如圖6所示。圖6中，“+”號表示未映射的種子點，“◇”號表示需映射的種子點，“°”號表示映射點。

圖6　篩選后的種子點和映射點分布情況

4程序實施及實例

基于Voronoi結構的多邊形單元網格程序實施方案的步驟如下：(Ⅰ)在目標區域內隨機生成n個初始種子點，利用初始種子點優化方法對種子點進行分布優化；(Ⅱ)利用邊界處理方法對步驟(Ⅰ)獲得的種子點Pin進行篩選，生成對應的映射點Pout；(Ⅲ)生成點集P=Pin∪Pout對應的Voronoi結構，利用Lloyd算法對Voronoi結構進行進一步優化，最后輸出種子點Pin對應的Voronoi單胞即為最終的多邊形單元網格。

實例在單位圓形區域內部挖去一個較小的圓形區域，其邊界方程表達式為：

利用上述方案生成的多邊形單元網格如圖7所示，其中，μ表示單元內角角度的平均值，σ表示角度的標準差。

圖7　基于Voronoi結構的多邊形單元網格實例

圖7a和圖7c表明：通過本文多邊形網格生成方法得到的均勻網格單元及非均勻網格單元形狀都較為良好。圖7b和圖7d是均勻網格單元和非均勻網格單元情況下的內角角度分布規律圖，不難發現單元內角角度平均值約為120°，標準差為5°～8°。通過對比均勻和非均勻網格中多邊形單元內角角度的統計值可以發現：均勻網格的角度分布更集中于120°，而非均勻網格則相對略為分散，但都近似正多邊形并基本符合正態分布的規律，從而驗證了本文多邊形網格生成方法的合理性和可行性。

5結束語

本文利用三角形網格生成算法對初始種子點分布進行合理優化，以控制網格中單元密度的分布，實現非均勻網格的生成。通過添加種子點對應目標區域邊界的映射點，基于種子點與其對應的映射點具有的共同Voronoi單胞邊界，構成目標區域的邊界，簡化了Voronoi結構對復雜邊界的逼近過程。并對Voronoi結構進行質心化處理，顯著改善了多邊形單元形態。最后通過實例驗證了本文方法可有效控制多邊形網格中單元的密度分布，且獲得的多邊形單元形態良好并近似正多邊形。

參考文獻：

[1]王兆清.有理單元法研究[D].上海:上海大學,2004.

[2]LV J,ZHANG H W,YANG D S.Multiscale method for mechanical analysis of heterogeneous materials with polygonal microstructures[J].Mechanics of materials,2013,56(1):38-52.

[3]盛國雨.基于徑向積分法的多邊形及多面體有限元方法[D].大連:大連理工大學,2015.

[4]王兆清,李淑萍.多邊形單元網格自動生成技術[J].中國圖象圖形學報,2007,12(7):1307-1311.

[5]SOHN D,CHO Y S,IM S.A novel scheme to generate meshes with hexahedral elements and poly-pyramid elements:the carving technique[J].Computer methods in applied mechanics & engineering,2012,201/204:208-227.

[6]蔡永昌,郭盛勇,楊健,等.多邊形單元的構造新方法及實現[J].同濟大學學報(自然科學版),2009,37(7):883-887.

[7]暴艷利,鐘紅,林皋.基于多邊形比例邊界有限元的重力壩地震斷裂模擬[J].水電能源科學,2015,33(4):72-75,42.

[8]王兆清,馮偉.自然單元法研究進展[J].力學進展,2004,34(4):437-445.

[9]DU Q,FABER V,GUNZBURGER M.Centroidal Voronoi tessellations:applications and algorithms[J].SIAM review,1999,41(4):637-676.

[10]DU Q,EMELIANENKO M,JU L.Convergence of the Lloyd algorithm for computing centroidal Voronoi tessellations[J].SIAM journal on numerical analysis,2006,44(1):102-119.

[11]WANG J,WANG X.Edge-weighted centroidal Voronoi tessellations[J].Numerical math (theory,methods and applications),2010,3:223-244.

[12]PERSSON P,STRANG G.A simple mesh generator in MATLAB[J].SIAM review,2004,46(2):329-345.

[13]TALISCHI C,PAULINO G H,PEREIRA A,et al.Polygonal finite elements for topology optimization:a unifying paradigm[J].International journal for numerical methods in engineering,2010,82(6):671-698.

基金項目：國家自然科學基金項目(51278169)

作者簡介：石怡婧(1992-),女,江蘇南通人,碩士生;邵國建(1962-),男,浙江臺州人,教授,博士,博士生導師,主要研究方向為計算力學與工程仿真.

收稿日期：2016-01-21

文章編號：1672-6871(2016)05-0051-05

DOI:10.15926/j.cnki.issn1672-6871.2016.05.012

中圖分類號：TP391

文獻標志碼：A