材料二維微結構仿真隨機概率圓優(yōu)化填充算法

陳　震， 周金宇

(江蘇理工學院常州市裝備再制造工程重點實驗室，江蘇 常州 213001)

材料二維微結構仿真隨機概率圓優(yōu)化填充算法

陳震， 周金宇

(江蘇理工學院常州市裝備再制造工程重點實驗室，江蘇 常州 213001)

幾何建模是影響有限元法計算材料性能精確性的主要因素之一。利用圓等簡單的幾何實體近似模擬離散化實體效果良好，不僅能大量簡化幾何計算，還能構造出誤差在可接受范圍內的近似模型。材料特性不僅取決于平均顆粒尺寸，受顆粒尺寸分布的影響也十分顯著。本文設計的隨機圓填充算法保證了填充圓半徑服從某一概率分布，重點分析了常規(guī)圓產生方法和依據(jù)干涉圓的干涉類型生成新圓。該算法不僅完全杜絕了圓的干涉還極大的減小了圓的相離，相對于現(xiàn)有的填充算法，填充度有了顯著提升。

隨機半徑；材料仿真；顆粒填充；有限元法

使用離散化單元法分析問題首先要采用離散化實體劃分的形式近似代替各種實際情況[1-2]。對于許多復雜的物理現(xiàn)象，利用圓等簡單的幾何實體近似模擬離散化實體效果良好。圓填充問題也由其學術價值和工業(yè)應用的重要性成為幾個世紀以來長期研究的課題。17世紀英國數(shù)學家牛頓就曾研究過關于等尺寸圓彼此不干涉時在某幾何形狀中的填充以獲得最大填充度的最優(yōu)幾何填充問題。該類算法不僅能簡化幾何計算，還能構造出誤差區(qū)間在可接受范圍內的模型。模型使得復雜的幾何結構變成簡單的幾何填充問題。隨著20世紀50年代計算機技術的迅速發(fā)展，計算機技術在材料仿真建模方面得廣泛應用[3]。

有限元技術作為解決材料力學問題最有力的數(shù)字化方法之一，其精確性主要與被分析材料的幾何模型有關[4]。在實際分析應用中，盡量提高離散化實體填充度的同時，還需控制其概率分布對材料性能的影響[5]。

材料微結構模型包括二維和三維兩類。二維模型采用圓填充、直線劃分圓心邊界的形式；三維模型則采用球填充、平面劃分球心邊界的形式，當前由空間點集生成三維模型的方法主要有增量擴展算法[6]和基于 Delaunary三角剖分算法[7]。二維微結構模型是三維模型的基礎，通過對二維模型的有效擴展可以構建出三維模型。基于二維微結構模型的細觀有限元分析，對評估材料表面微裂紋的萌生和擴展，推演二維微結構復合材料損傷演化行為具有重要意義[8-9]。

Voronoi算法是常見的隨機晶粒生成算法[10]，其定義為：設P={p1, p2, ···, pn}是二維歐氏幾何空間平面上的n個互不相同的點集，d(p, pi)表示點p與pi之間的歐幾里德距離，由劃分的平面稱作以為點集的 Voronoi圖，V(pi)為Voronoi晶胞，pi為各晶胞的生成子。Voronoi多采用在某平面根據(jù)生成晶胞的平均面積確定生成子數(shù)量，然后隨機撒點并劃分點的邊界生成晶胞，該方式比較容易控制平均晶胞面積尺寸，但無法控制尺寸分布使其滿足多樣化的統(tǒng)計規(guī)律[5]，具有單一性和不可調節(jié)性，且主要適應于晶粒均勻形核并等速生長的退火狀態(tài)純金屬，應用范圍較窄。

故有學者采用 Voronoi圖的推廣 Laguerre圖[5,11]進行建模。Laguerre圖多采用增量式算法生成點集，然后通過對這些點集賦予權值形成權圓，但由于點位置的隨機性就決定了權圓會發(fā)生干涉或相離的情況。干涉會減小權圓面積，相離則會產生過多的未填充區(qū)域，導致后繼的多邊形邊界劃分失真，產生較大的材料建模誤差，如圖 1所示。

圖1　二維Laguerre圖[11]

本文提出的圓填充算法能根據(jù)給定的概率密度函數(shù)生成隨機圓填充矩形區(qū)域。文中通過生成隨機半徑、排除干涉逐層生成的方式，可以實現(xiàn)半徑依據(jù)要求概率密度函數(shù)生成概率圓且沒有干涉。同時，由于每個生成的新圓至少與已存在的兩個圓或邊界相切，故填充度明顯高于傳統(tǒng)Laguerre權圓填充度，只要根據(jù)相關材料的邊界性質進行邊界劃分，即可構建較高精度的二維材料微結構仿真模型。

1　常規(guī)顆粒圓生成及活動鏈更新

該算法采用的是自底向上逐層生成的方法。假設填充區(qū)域的左、右、上、下邊界分別為Bl，Br，Bt，Bb。首先根據(jù)指定的概率密度函數(shù)取隨機半徑r1產生圓C11，使C11同時與邊界Bl和Bb相切，然后依次生成C12, C13, ···, C1n，使新生成的圓分別與前一個圓和邊界Bb相切直至新圓Cn與邊界Br干涉，再以r1n為半徑生成新圓使其同時與C1(n-1)和邊界 Br相切取代圓 C1n，至此第一層圓生成完畢，同時形成初始鏈 C11→C12→C13→…→C1(n–1)→C，如圖2(a)所示。該鏈是由具有確定方向的線段 C11→C12, C12→C13,…, C1(n–1)→C1n’組成并可以隨著新圓的產生而更新。

初始鏈形成之后，根據(jù)指定的概率密度函數(shù)產生隨機半徑r21作新圓C21，使其同時與C11和邊界Bl相切，然后添加新的有向線段C21C11，更新成新的活動鏈C21→C11→C12→C13→…→C1(n–1)→C1n’，如圖2(b)所示。然后按序首先選取C11C12作為當前線段，產生隨機半徑 r22作新圓 C22，要求該圓與C11和C12相切并居于有向線段C11C12的左側，然后從初始活動鏈上刪除線段C11C12，并添加新的有向線段 C11C22和 C22C12，更新成新的活動鏈C21→C11→C22→C12→C13→…→C1(n–1)→C1n’，如圖2(c)所示。下一步選取線段C12C13作為當前線段，以相同的方式產生C23，并從初始活動鏈上刪除線段C12C13，再添加新的有向線段C12C23和C23C13，形成新的活動鏈C21→C11→C22→C12→C23→C13→…→C1(n-1)→C1n’，如圖 2(d)所示。依此逐步生成新圓并更新活動鏈直至本層圓生成完畢。

圖2　新圓生成及活動鏈更新

2　顆粒圓干涉的處理方法

2.1新開始層第1個圓的干涉及活動鏈更新

新開始層的第 1個圓可能會與已生成的圓存在干涉。如圖 3(a)所示，活動鏈為 C11→C12→C13→…，以r*為半徑產生的圓Ci同時與C11和邊界 Bl相切且與 C12干涉，則再次以隨機半徑 r*生成新圓Cz，要求該圓同時與C12和邊界Bl相切。再次判斷該圓是否與活動鏈上的圓干涉。如果不干涉，則當前活動鏈更新為Cz→C11→C12→C13→…；如果依然存在干涉，則活動鏈不更新。然后按序首先選取 C11C12作為當前線段，進一步生成新圓。

2.2與邊界Br干涉及活動鏈更新

隨著沿活動鏈正方向按序選取當前線段形成新圓，該圓可能會與邊界Br干涉。如圖3(b)所示，活動鏈為…→Cp(k–2)→Cp(k–1)→Cpk→C2m→C1n，以C2m→C1n為當前線段，r*為半徑產生的圓Ci與邊界Bl干涉，則以隨機半徑 r*生成新圓 Cz，首先要求該圓與C2m和邊界Bl相切，判斷該圓是否與活動鏈上的圓干涉。如果不干涉，則當前活動鏈更新為…→Cp(k–2)→Cp(k–1)→Cpk→C2m→Cz，開始新層；如果存在干涉，再次以隨機半徑 r*生成新圓 Cz，并要求該圓與 C1n和邊界 Bl相切，判斷該圓是否與活動鏈上的圓干涉。如果不干涉，則當前活動鏈更新為…→Cp(k–2)→Cp(k–1)→Cpk→C2m→Cz，開始新層；如果依然存在干涉，則活動鏈不更新，按序選取下一線段作當前線段。由于圖 3(b)中C2m→C1n為活動鏈的最后一個線段，所以本層結束，開始新層。

因此，中國版南海爭端國際意象的塑造和推廣還需通過具體的外交政策和精準的國際輿論引導才能真正獲得相關國家的認同。一方面，使建立在客觀、尊重歷史源流、符合國際法的公平公正基礎上的“中國版”南海問題國際意象為各方所周知乃至認同；另一方面，對于在南海問題上的誤導性意象必須予以堅決回應乃至駁斥以正視聽。

圖3　新層首圓及新圓與邊界Br的干涉處理

2.3根據(jù)當前線段生成的新圓只與已存在圓的干涉及處理

當沿活動鏈正方向前進，根據(jù)當前線段生成的新圓可能會只與已生成的圓存在干涉。此時沿當前線段正方向前進，定義當前線段前方的圓為前置圓，其后方的圓為后置圓。處理此干涉類型需要選取新的當前線段替代原當前線段，可分為以下2種干涉情況：

(1) 與前置圓干涉。如圖4(b)所示，以CkCk+1為當前線段隨機生成的圓 Ci(半徑為 r*)與前置圓Ck+n干涉，則依次生成當前線段CkCk+n，Ck+1Ck+n。根據(jù)新的當前線段 CkCk+n以隨機半徑 r*生成新圓Cz，判斷該圓是否與活動鏈上的圓或邊界干涉。如果不干涉，則在當前活動鏈中刪除線段 CkCk+1, Ck+1Ck+2,…, Ck+n–1Ck+n，并添加新的有向線段CkCz和CzCk+n，形成新的活動鏈；如果干涉，則再選取當前線段 CkCk+n，同理生成新圓并檢查干涉，如果不干涉，則在當前活動鏈中刪除線段Ck+1Ck+2,…, Ck+n–1Ck+n，并添加新的有向線段Ck+1Cz′和Cz′Ck+n，形成新的活動鏈；如果依然干涉，則跳過原當前線段 CkCk+1，選取下一當前線段Ck+1Ck+2。

圖4　干涉情況

(2) 與后置圓干涉。如圖4(c)所示，以CkCk+1為當前線段隨機生成的圓 Ci(半徑為 r*)與前置圓Ck–n干涉，則將當前線段更新為 Ck–nCk+1，根據(jù)新的當前線段 Ck–nCk+1再次以隨機半徑 r*生成新圓Cz，判斷該圓是否與活動鏈上的圓或邊界干涉。如果不干涉，則在當前活動鏈中刪除線段Ck–nCk–n+1,…, Ck–1Ck, CkCk+1，并添加新的有向線段Ck-nCz和CzCk+1，形成新的活動鏈；如果干涉，則再選取當前線段 Ck–nCk，同理生成新圓并檢查干涉，如果不干涉，則在當前活動鏈中刪除線段Ck–nCk–n+1,…, Ck–1Ck，并添加新的有向線段Ck–n和′Ck，形成新的活動鏈；如果依然干涉，則跳過原當前線段CkCk+1，選取下一當前線段Ck+1Ck+2。

3　擴展鏈算法

構造隨機權圓擴展鏈算法的主流程如圖 5所示，具體步驟如下：步驟1. 以initial_chain子程序生成初始圓。步驟2. 判斷是否滿足區(qū)域的填充要求。如果滿足則進入步驟17。

步驟3. 判斷是否是新層。如果不是則進入步驟15。

步驟4. 以random_r子程序生成隨機半徑r*，生成新層首圓。

步驟6. 更新活動鏈。

步驟7. 接受新圓并用circle子程序生成該圓。步驟8. 按序選取活動線段，回到步驟2。

步驟9. 更新干涉矩陣，即把與新圓干涉的圓數(shù)據(jù)依次存入干涉矩陣，數(shù)據(jù)包括各圓心坐標及其半徑，如圖4(b)所示，以CkCk+1為當前線段生成的新圓Ci與Ck+n等s個圓干涉，則把這s個圓Ck+n, Ck+n–1, … , Ck+n–s+1的數(shù)據(jù)依次存入干涉矩陣。

步驟10. 按序讀取干涉矩陣中各圓心坐標元素。

步驟11. 根據(jù)第1節(jié)中干涉類型生成新圓并判斷其是否與活動鏈上已存在的圓或邊界干涉，如果不干涉則進入步驟14。

步驟12. 判斷標示量t是否為1，如果是，則進入步驟8。

步驟13. 判斷干涉矩陣是否讀取完畢，如果沒有，回到步驟10；如果讀完，回到步驟8。

步驟14. 根據(jù)第1節(jié)中各干涉類型和情況更新活動鏈，回到步驟7。

步驟15. 生成新圓并判斷其是否與邊界Br干涉，如果干涉則另標示量t=1，進入步驟11。

步驟 16. 判斷新圓是否與活動鏈上圓干涉，如果干涉進入步驟9；否則進入步驟6。

步驟17. 程序結束。

圖5　隨機顆粒圓生成算法流程圖

算法采用自下向上逐層生成的方法產生概率圓進行區(qū)域填充。算法執(zhí)行產生 n層隨機圓；每層產生隨機圓的計算次數(shù)為mi(i=1, 2, 3,…, n)，mi為第i層概率圓數(shù)量；其中第i層有ki個圓需要進行干涉處理，則運行上限為n*mmax*kmax次；常規(guī)概率圓的生成運行上限為n*mmax次；此外，循環(huán)外附加的初始鏈生成算法、原始條件賦予和存儲矩陣建立等雖然會帶來額外的計算量，但相對于循環(huán)內考慮概率圓干涉情況的計算量幾乎可以忽略，不影響實際復雜度。綜上所述，擴展鏈算法的時間復雜度為O(n3)。

4　數(shù)值算例

基于 4個指定概率密度函數(shù)進行概率圓填充實驗，概率密度函數(shù)分別為：

根據(jù)指定概率密度函數(shù)生成的概率圓分別如圖6(a)～(d)所示。各圓之間完全杜絕了干涉現(xiàn)象，且每個圓都至少與兩個圓相切，填充度高。表 1中的數(shù)據(jù)1～4依次對應4種概率密度函數(shù)。期望值為概率密度函數(shù)確定的圓半徑統(tǒng)計參數(shù)，實際值為根據(jù)相應概率密度函數(shù)生成的圓半徑的統(tǒng)計參數(shù)，填充度為區(qū)域內圓的面積之和與矩形區(qū)域面積的百分比值。

圖6　根據(jù)概率分布生成的概率圓

表1　概率圓統(tǒng)計參數(shù)表

5　結 束 語

本文詳細敘述了圓填充算法的實現(xiàn)過程，該算法可以根據(jù)給定的概率密度函數(shù)生成隨機圓填充區(qū)域。根據(jù)給定的概率密度函數(shù)生成隨機半徑、排除干涉逐層生成新圓。算法避免了生成過程中圓的干涉情況，并且生成圓的半徑完全符合要求概率密度函數(shù)。同時，由于新圓在生成時至少與已生成的兩個圓或邊界相切，故填充度較現(xiàn)有算法有明顯提升，提高了仿真建模應用的可靠性和真實性。

[1] 高井望, 徐佩華, 黃潤秋, 等. DEM 擬合級配碎石材料大三軸試驗的顆粒粗糙度效應研究[J]. 工程地質學報, 2013, 21(4): 560-569.

[2] Gong G B, Zha X X, Wei J. Comparison of granular material behaviour under drained triaxial and plane strain conditions using 3D DEM simulation [J]. Engineering Computation, 2012, 25(2): 186-196.

[3] 秦湘閣, 劉國權. 顆粒增強復合材料顯微組織的計算機仿真[J]. 佳木斯大學學報, 2001, 19(1): 1-3.

[4] 崔文凱, 馮仰德, 紀國良, 等. 兩種復合材料幾何建模算法[J]. 計算機輔助設計與圖形學學報, 2015, 27(1): 46-50.

[5] 李俊琛. 復合材料微觀組織結構的仿真與計算[D].蘭州: 蘭州理工大學, 2009.

[6] 孫殿柱, 劉健, 李延瑞，等. 三維散亂點集Voronoi圖快速生成算法研究[J]. 武漢科技大學學報, 2010, 35(8): 909-912.

[7] 范志剛, 吳裕功, 趙選賀. Laguerre-Voronoi圖軟件包的設計和實現(xiàn)[J]. 天津大學學報, 2003, 36(6): 747-752.

[8] 任淮輝, 李旭東. 二維材料微結構設計與數(shù)值模擬軟件系統(tǒng)[J]. 計算機輔助設計與圖形學學報, 2009, 32(1): 45-50.

[9] 任淮輝, 李旭東, 李俊琛, 等. 二維多晶體材料微結構的力學響應計算[J]. 武漢科技大學學報, 2008, 31(4): 427-431.

[10] 司良英, 鄧關宇, 呂程, 等. 基于Voronoi圖的晶體塑性有限元多晶幾何建模[J]. 材料與冶金學報, 2009, 8(3):193-216.

[11] 張賦, 李旭東. 帶權點集Laguerre圖的增量算法與軟件設計研究[J]. 計算機工程與應用, 2012, 48(30): 10-18.

Random Circles Optimization Filling Algorithm of Material Two Dimensional Microstructures Simulation

Chen Zhen,Zhou Jinyu
(Changzhou Key Laboratory of Equipment Remanufacturing Engineering, Jiangsu University of Technology, Changzhou Jiangsu 213001, China)

Geometric modeling is one of the main factors affecting the accuracy of material properties calculation by finite element method. Using simple geometric entities, such as circles, to simulate discrete objects is demonstrated well. Beyond the great simplifications in the geometrical calculation, it can provide an approximate model within acceptable error. Material properties not only depend on the average grain size, but also are influenced by the grain size distribution significantly. Random circles packing algorithm designed in this paper ensured these circles′ radius following a certain probability distribution. Normal circles generation method and generating new circle based on the interference types were analyzed. This algorithm not only avoided the interference completely, but reduced the separation greatly. Compared with existing filling algorithms, the filling density of this algorithm was improved obviously.

random radius; material simulation; sphere packing; finite element method

TP 391.9

A

2095-302X(2015)06-0944-06

2015-05-19；定稿日期：2015-08-03

國家自然科學基金資助項目(51275221)；江蘇省產學研聯(lián)合創(chuàng)新資金資助項目(BY2014038-04)

陳震(1990–)，男，江蘇徐州人，碩士研究生。主要研究方向為先進(再)制造技術與裝備。E-mail：ericchen31@163.com

周金宇(1973–)，男，江蘇揚州人，教授，博士。主要研究方向為結構可靠性、機械疲勞與延壽、現(xiàn)代設計方法。E-mail：yuhangyuan888@sina.com