石怡婧，邵國建，丁勝勇

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

?

基于Voronoi結(jié)構(gòu)的多邊形單元網(wǎng)格生成方法

石怡婧，邵國建，丁勝勇

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

摘要：針對多邊形單元網(wǎng)格難以生成的問題，建立了基于Voronoi結(jié)構(gòu)的多邊形單元網(wǎng)格生成方法。該方法通過區(qū)域內(nèi)的一組初始點來構(gòu)造Voronoi結(jié)構(gòu)。對初始種子點的分布進行優(yōu)化，以達到控制多邊形單元網(wǎng)格密度的目的。利用Voronoi結(jié)構(gòu)的特點，通過添加種子點對應目標區(qū)域邊界的映射點對邊界進行擬合，實現(xiàn)Voronoi結(jié)構(gòu)對復雜邊界的逼近。對Voronoi結(jié)構(gòu)進行質(zhì)心化，改善生成多邊形單元形態(tài)。給出了該方法的程序?qū)嵤┎襟E，并結(jié)合網(wǎng)格生成實例，驗證了其合理性和可行性。

關(guān)鍵詞：有限元法；多邊形單元；網(wǎng)格生成；Voronoi結(jié)構(gòu)

0引言

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

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

1Voronoi結(jié)構(gòu)及其質(zhì)心化

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

(1)

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

圖1　Voronoi圖及其Delaunay三角化

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

(2)

(3)

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

(4)

此時，生成的Voronoi圖為CVT。

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

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

圖2　基于Lloyd算法的Voronoi結(jié)構(gòu)優(yōu)化過程

2初始種子點分布優(yōu)化

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

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

圖3　帶圓孔的正方形區(qū)域內(nèi)初始種子點的優(yōu)化分布

3Voronoi結(jié)構(gòu)邊界處理

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

(5)

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

圖4 種子點及其對應的映射點圖5 Voronoi結(jié)構(gòu)的邊界處理

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

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

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

4程序?qū)嵤┘皩嵗?/p>

基于Voronoi結(jié)構(gòu)的多邊形單元網(wǎng)格程序?qū)嵤┓桨傅牟襟E如下：(Ⅰ)在目標區(qū)域內(nèi)隨機生成n個初始種子點，利用初始種子點優(yōu)化方法對種子點進行分布優(yōu)化；(Ⅱ)利用邊界處理方法對步驟(Ⅰ)獲得的種子點Pin進行篩選，生成對應的映射點Pout；(Ⅲ)生成點集P=Pin∪Pout對應的Voronoi結(jié)構(gòu)，利用Lloyd算法對Voronoi結(jié)構(gòu)進行進一步優(yōu)化，最后輸出種子點Pin對應的Voronoi單胞即為最終的多邊形單元網(wǎng)格。

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

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

圖7　基于Voronoi結(jié)構(gòu)的多邊形單元網(wǎng)格實例

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

5結(jié)束語

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

參考文獻：

[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]王兆清,李淑萍.多邊形單元網(wǎng)格自動生成技術(shù)[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]蔡永昌,郭盛勇,楊健,等.多邊形單元的構(gòu)造新方法及實現(xiàn)[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