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

石怡婧，邵國建，丁勝勇

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

?

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

石怡婧，邵國建，丁勝勇

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

摘要：針對多邊形單元網(wǎng)格難以生成的問題，建立了基于Voronoi結(jié)構(gòu)的多邊形單元網(wǎng)格生成方法。該方法通過區(qū)域內(nèi)的一組初始點來構(gòu)造Voronoi結(jié)構(gòu)。對初始種子點的分布進行優(yōu)化，以達(dá)到控制多邊形單元網(wǎng)格密度的目的。利用Voronoi結(jié)構(gòu)的特點，通過添加種子點對應(yīng)目標(biāo)區(qū)域邊界的映射點對邊界進行擬合，實現(xiàn)Voronoi結(jié)構(gòu)對復(fù)雜邊界的逼近。對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ì)的復(fù)雜性，勢必會增加建立有限元網(wǎng)格模型的難度。尤其對于形狀比較復(fù)雜的幾何形狀，其網(wǎng)格剖分耗時較多，并且容易出錯。因此，網(wǎng)格剖分技術(shù)是網(wǎng)格模型建立的核心內(nèi)容，一個好的有限元網(wǎng)格對提高計算效率和計算精度起著關(guān)鍵作用。目前，對于任意形狀的分析區(qū)域，三角形及四邊形單元的網(wǎng)格生成技術(shù)已經(jīng)較為成熟。而多邊形單元作為一種新型的單元形式，由于其邊數(shù)及形狀的任意性，與常規(guī)單元相比有獨特的優(yōu)勢。針對復(fù)雜幾何體進行網(wǎng)格劃分更加靈活、方便[1-2]，近年來受到越來越廣泛的關(guān)注[3]。文獻(xiàn)[4]提出了Delaunay多邊形化生成多邊形網(wǎng)格的技術(shù)。文獻(xiàn)[5]將一種類似于金字塔型的方法用于劃分復(fù)雜三維模型多邊形網(wǎng)格中。文獻(xiàn)[6]利用傳統(tǒng)三角形有限元網(wǎng)格生成形狀合理的多邊形單元網(wǎng)格。文獻(xiàn)[7]利用多邊形比例邊界有限元方法生成多邊形網(wǎng)格。但是，目前多邊形單元的網(wǎng)格生成算法還比較粗糙，生成的多邊形單元網(wǎng)格非?；?，處理復(fù)雜邊界處的網(wǎng)格時，構(gòu)造能較好滿足多邊形有限元分析的單元形狀也有一定困難。

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

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

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

(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]，其種子點與對應(yīng)的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ū)域采用更精細(xì)的網(wǎng)格，對一些非重點區(qū)域采用較為粗略的網(wǎng)格。一般初始種子點的生成都是采用隨機布點的方法，這種方法會使多邊形單元密度分布呈現(xiàn)隨機化，難以滿足有限元分析的要求。初始種子點的密度會在很大程度上影響最后CVT的密度，所以本文對初始種子點的分布進行優(yōu)化。

本文初始種子點分布優(yōu)化采用文獻(xiàn)[12]提出的三角形網(wǎng)格生成算法，其基本思想是將三角形網(wǎng)格類比為平面彈簧結(jié)構(gòu)，三角形的三條邊分別對應(yīng)連接種子點的各彈簧，每根彈簧的現(xiàn)有長度l和期望長度l0遵循力-位移關(guān)系f(l,l0)，并通過改變種子點的位置使整個結(jié)構(gòu)達(dá)到平衡。圖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)邊界處理

因此，可通過下式得到種子點對應(yīng)的映射點：

(5)

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

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

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

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

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

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

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

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

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

圖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°，標(biāo)準(zhǔn)差為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)格的生成。通過添加種子點對應(yīng)目標(biāo)區(qū)域邊界的映射點，基于種子點與其對應(yīng)的映射點具有的共同Voronoi單胞邊界，構(gòu)成目標(biāo)區(qū)域的邊界，簡化了Voronoi結(jié)構(gòu)對復(fù)雜邊界的逼近過程。并對Voronoi結(jié)構(gòu)進行質(zhì)心化處理，顯著改善了多邊形單元形態(tài)。最后通過實例驗證了本文方法可有效控制多邊形網(wǎng)格中單元的密度分布，且獲得的多邊形單元形態(tài)良好并近似正多邊形。

參考文獻(xiàn)：

[1]王兆清.有理單元法研究[D].上海:上海大學(xué),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].大連:大連理工大學(xué),2015.

[4]王兆清,李淑萍.多邊形單元網(wǎng)格自動生成技術(shù)[J].中國圖象圖形學(xué)報,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].同濟大學(xué)學(xué)報(自然科學(xué)版),2009,37(7):883-887.

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

[8]王兆清,馮偉.自然單元法研究進展[J].力學(xué)進展,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.

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

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

收稿日期：2016-01-21

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

DOI:10.15926/j.cnki.issn1672-6871.2016.05.012

中圖分類號：TP391

文獻(xiàn)標(biāo)志碼：A