一種基于Voronoi圖的多晶體有限元建模方法

鄭戰(zhàn)光，汪兆亮，馮　強(qiáng)，袁　帥，王佳祥

(1.廣西大學(xué)機(jī)械工程學(xué)院, 廣西南寧530004；2.廣西制造系統(tǒng)與先進(jìn)制造技術(shù)重點(diǎn)實(shí)驗(yàn)室， 廣西南寧530004)

?

一種基于Voronoi圖的多晶體有限元建模方法

鄭戰(zhàn)光1,2，汪兆亮1，馮強(qiáng)1，袁帥1，王佳祥1

(1.廣西大學(xué)機(jī)械工程學(xué)院, 廣西南寧530004；2.廣西制造系統(tǒng)與先進(jìn)制造技術(shù)重點(diǎn)實(shí)驗(yàn)室， 廣西南寧530004)

摘要:建立多晶體的細(xì)觀有限元模型是研究多晶體材料局部塑性變形不均勻性的前提與基礎(chǔ)，為了靈活地構(gòu)建高可靠度的材料微結(jié)構(gòu)模型，在前人研究成果的基礎(chǔ)上提出一種基于Voronoi圖并結(jié)合單元編號(hào)區(qū)域排布特點(diǎn)，能直接根據(jù)模型中得到的單元編號(hào)順序依次求取單元形心坐標(biāo)的構(gòu)圖方法就顯得極為關(guān)鍵。該方法首先是生成特定平面或空間域里的隨機(jī)點(diǎn)與Voronoi圖基本信息，再結(jié)合單元編號(hào)區(qū)域排布特點(diǎn)依次直接求取中心點(diǎn)坐標(biāo)，接著判斷單元?dú)w屬于距離最近的晶核所在的晶粒內(nèi)，并將所得晶粒編號(hào)及單元編號(hào)以set集合形式添加到INP文件的Part部分，進(jìn)而得到Voronoi多晶體有限元模型，最后以構(gòu)建含10個(gè)晶粒的二維和三維多晶體模型為例和文獻(xiàn)對(duì)比分析來進(jìn)行實(shí)現(xiàn)與驗(yàn)證。結(jié)果顯示：該方法可以依據(jù)單元編號(hào)區(qū)域排布特點(diǎn)直接得到單元編號(hào)且更加容易實(shí)現(xiàn)依次求取單元形心坐標(biāo)，并在一定程度上降低了單元形心坐標(biāo)處理的數(shù)據(jù)量和單元?dú)w屬判斷的難度，通過對(duì)比分析該方法建立的模型精度更接近于文獻(xiàn)中的精確模型，它們之間的最大偏差僅為25.47 MPa，較對(duì)比文獻(xiàn)的簡(jiǎn)化模型最小偏差還要低0.07 MPa。表明該方法可為研究人員快速構(gòu)建多晶材料的Voronoi細(xì)觀有限元模型提供一定的技術(shù)參考。

關(guān)鍵詞:Vonronoi圖；多晶體建模；晶體塑性；有限元；

0引言

隨著計(jì)算幾何學(xué)和計(jì)算機(jī)圖形技術(shù)的快速發(fā)展以及國(guó)內(nèi)外在細(xì)觀尺度下研究金屬不均勻變形與內(nèi)部晶粒演變之間關(guān)系等熱點(diǎn)問題的興起，學(xué)者們紛紛提出了一些方法來模擬多晶體位相分布的方法。目前應(yīng)用較多的方法有Voronoi圖方法、蒙特卡羅方法和元胞自動(dòng)機(jī)方法等，其中，Voronoi圖方法可以與有限元方法相耦合來建立多晶有限元模型，所建立的模型能靈活地分析不同材料微結(jié)構(gòu)演變的物理機(jī)制，而且存儲(chǔ)數(shù)據(jù)少，計(jì)算速度快，因此Voronoi圖方法在國(guó)內(nèi)外得到了廣泛的應(yīng)用[1-5]?；赩oronoi圖的多晶體有限元建模常采用如下兩類方法：①先生成Voronoi晶胞的基本信息，將其按一定順序保存為數(shù)據(jù)文件，然后通過Python腳本語(yǔ)言對(duì)ABAQUS進(jìn)行二次開發(fā)就可以生成Voronoi多晶材料的細(xì)觀有限元模型。如司良英等[6]通過先在MATLAB生成Voronoi晶胞，將其按一定的順序保存為數(shù)據(jù)文件，然后在ABAQUS/CAE中用其與Python接口將其讀入，調(diào)用ABAQUS/CAE中Python程序接口，最終生成所需的有限元多晶幾何模型；樊黎霞等[7]通過MATLAB里的Voronoi函數(shù)把各晶粒的頂點(diǎn)和晶粒順序號(hào)存儲(chǔ)到文件中，通過Python語(yǔ)言讀入到ABAQUS有限元軟件建立的幾何模型，按照順序依次連接晶粒的各頂點(diǎn)，然后采用分割面的功能生成Voronoi多晶有限元幾何模型；張豐果等[8]先在區(qū)域R內(nèi)生成Voronoi晶胞的隨機(jī)種子點(diǎn)，接著對(duì)區(qū)域R預(yù)先劃分的每一個(gè)網(wǎng)格單元采用最短距離的方法判斷其歸屬于哪個(gè)晶粒，將得到的基本信息添加到Python腳本文件，最終建立了基于Voronoi圖的粗糙多晶模型；汪凱[9]采用Python腳本語(yǔ)言getNodes命令，先獲得每個(gè)單元所包含的節(jié)點(diǎn)坐標(biāo)，然后計(jì)算每個(gè)單元形心與所有種子點(diǎn)的距離，通過排序函數(shù)getMin得到最小數(shù)的編號(hào)，就將該單元放到該種子所對(duì)應(yīng)的集合(晶粒)中，即可得到三維Voronoi多晶簡(jiǎn)化模型。②先生成Voronoi圖的基本信息，再將獲得的Voronoi幾何信息添加到ABAQUS預(yù)先生成的單晶模型INP文件中，最終得到多晶體塑性變形有限元模型。如伊興華[10]利用Voronoi圖表技術(shù)，先通過MATLAB編程得到多晶的拓?fù)湫畔?，然后在ABAQUS軟件的前處理模塊CAE中建立模型的幾何信息部分，并結(jié)合MATLAB程序得到的Voronoi信息，最終得到晶體塑性模型的ABAQUS輸入文件。

上述兩類方法在獲取Voronoi圖的基本信息時(shí)均涉及到對(duì)單元形心坐標(biāo)的處理，這樣因選擇的參考坐標(biāo)系不同或因直接采用先求單元節(jié)點(diǎn)坐標(biāo)、再求單元形心坐標(biāo)的方法均會(huì)加大數(shù)據(jù)的處理量，同時(shí)也會(huì)對(duì)后面單元?dú)w屬的判斷增加一定的難度。據(jù)此，本文提出一種基于Voronoi圖并結(jié)合單元編號(hào)區(qū)域排布特點(diǎn)，直接根據(jù)模型中得到的單元編號(hào)順序，并依次求取單元形心坐標(biāo)的構(gòu)圖方法。該方法是先生成特定平面或空間域里的隨機(jī)點(diǎn)與Voronoi圖基本信息，接著結(jié)合單元編號(hào)區(qū)域排布特點(diǎn)依次直接求取中心點(diǎn)坐標(biāo)，判斷單元?dú)w屬，將所得晶粒編號(hào)及單元編號(hào)以set集合形式添加到ABAQUS的INP文件Part部分，進(jìn)而得到Voronoi多晶體有限元模型。

1Voronoi多晶體模型的構(gòu)建

金屬材料晶粒實(shí)際上大小不一、形狀各不相同，在傳統(tǒng)的晶體塑性有限元模擬中常用一個(gè)規(guī)則的多面實(shí)體單元代表一個(gè)晶粒(單胞)，這種模型假設(shè)無法真實(shí)描述材料的微觀織構(gòu)，且不能夠反應(yīng)晶體內(nèi)部的不均性變形[11]。鑒于Voronoi圖技術(shù)算法已經(jīng)非常成熟，在MATLAB的Multi-Parametric Toolbox(MPT)工具箱有專門構(gòu)建Voronoi圖二維與多維的函數(shù)。但是如何通過MATLAB-MPT工具箱并結(jié)合MATLAB編程得到帶有邊界限制的完整的二維與三維Voronoi圖的幾何拓?fù)湫畔⑹墙oronoi多晶有限元模型最為關(guān)鍵的技術(shù)難題；同時(shí)晶粒數(shù)目較多時(shí)，如何將處理好的晶粒編號(hào)及所包含的單元信息以建立集合的方式寫進(jìn)INP文件中來獲得Voronoi多晶有限元模型也需要技術(shù)攻關(guān)。為此，本文提出一種基于Voronoi圖并結(jié)合單元編號(hào)區(qū)域排布特點(diǎn)，直接根據(jù)模型中得到的單元編號(hào)順序，并依次求取單元形心坐標(biāo)的構(gòu)圖方法，為研究人員在ABAQUS平臺(tái)下進(jìn)行Voronoi多晶材料的細(xì)觀有限元模型的建模提供一些技術(shù)參考。

1.1Voronoi圖的簡(jiǎn)介

圖1　Voronoi方法的原理圖Fig.1　The principle diagram of Voronoi method

Voronoi方法廣泛應(yīng)用于氣象學(xué)、地理學(xué)、圖像處理、結(jié)晶學(xué)、微觀結(jié)構(gòu)模擬和城市規(guī)劃等領(lǐng)域[1-6]。它是由俄羅斯數(shù)學(xué)家Voronoi于1908年提出了n維Voronoi圖的定義,其主要思路是通過臨近原則將n維空間體進(jìn)行剖分，使其成為無數(shù)多面體的集合體(如二維情況即將空間平面剖分成無數(shù)多邊形的集合體)，每個(gè)多面體事實(shí)上是依靠其內(nèi)部的一個(gè)核心點(diǎn)控制生成的，每個(gè)核心點(diǎn)的影響域是由至該點(diǎn)的距離最近的點(diǎn)組成的集合(臨近原則)，該核心點(diǎn)的影響域即構(gòu)成相應(yīng)的多面體[1-6]。生成Voronoi圖的初始點(diǎn)稱作Voronoi圖的發(fā)生元。從定義中可以看出，Voronoi圖就是基于距離最近原則對(duì)空間劃分而成的多邊形(二維)或多面體(三維)圖形集合，每個(gè)多邊形內(nèi)部任意一點(diǎn)到該多邊形對(duì)應(yīng)種子的距離都比到其他種子點(diǎn)的距離近，其主要描述了空間點(diǎn)的鄰近區(qū)域或者影響區(qū)域的邊界。Voronoi方法的原理圖如圖1所示。

1.2二維Voronoi模型的構(gòu)圖方法

按照如圖1所示的Voronoi方法的原理圖，在MATLAB中先生成種子點(diǎn)的隨機(jī)坐標(biāo)，依次求解每個(gè)單元的中心點(diǎn)坐標(biāo)到每個(gè)種子點(diǎn)(晶核)的距離，并將該單元?dú)w屬到距離最短的晶核所在的晶粒內(nèi)，從而得到每個(gè)晶粒的完整拓?fù)湫畔?，接著將得到的每一個(gè)晶粒編號(hào)及所對(duì)應(yīng)的單元編號(hào)，以set集合的形式寫進(jìn)ABAQUS的INP文件中Part部分最終實(shí)現(xiàn)整個(gè)二維Voronoi多晶模型的建立。

二維Voronoi多晶模型完整建模的程序流程圖如圖2所示，具體實(shí)施步驟如下：

圖2　二維Voronoi多晶模型建立的流程圖

①先在特定空間剖分得到晶核坐標(biāo)的隨機(jī)數(shù)，利用MATLAB中MPT工具箱提供的mpt_voronoi命令，生成所需大小和晶粒數(shù)的Voronoi圖基本信息并畫出二維模型圖。

②結(jié)合單元編號(hào)區(qū)域排布特點(diǎn)，尤其是單元區(qū)域分布及編號(hào)順序，編程直接求取每一個(gè)單元的中心點(diǎn)坐標(biāo)，單元是通過在ABAQUS/CAE中建立一個(gè)與二維模型圖相同尺寸的2D/Deformable/Shell模型劃分網(wǎng)格之后所得，單元具體信息可以在生成的INP文件中查看。

③根據(jù)ABAQUS/CAE中劃分網(wǎng)格后實(shí)體模型的單元空間排布特點(diǎn)，依次求解每個(gè)單元的中心點(diǎn)坐標(biāo)到每個(gè)晶核的距離，并將該單元?dú)w屬到距離最短的晶核所在的晶胞內(nèi)。

④將得到的每一個(gè)晶粒編號(hào)及其對(duì)應(yīng)的所有單元編號(hào)，最終以多級(jí)列表形式保存到txt文件中。

通過以上MATLAB編程，本文僅得到的是二維Voronoi圖完整的拓?fù)湫畔ⅰＯ旅孢€需要在ABAQUS的INP文件中的Part部分編寫多個(gè)set集合，其中每個(gè)set集合為一個(gè)晶粒，將txt文件中對(duì)應(yīng)的單元編號(hào)依次寫入相應(yīng)的set集合實(shí)現(xiàn)建模。利用INP文件進(jìn)行二維模型建立的過程如下：

①在ABAQUS/CAE中建立一個(gè)2D/Deformable/Shell，裝配，劃分網(wǎng)格單元(四邊形網(wǎng)格單元)，并在Job模塊中生成INP文件。

②將上一步寫好的set集合添加到建立的INP文件中的Part部分。在ABAQUS中打開修改后的INP文件，并在Color Code Dialog中設(shè)置為Color code by：Sets。

1.3三維Voronoi模型的構(gòu)圖方法

以上建立的二維Voronoi模型中每個(gè)晶粒為一個(gè)多邊形，所以二維Voronoi多晶模型的幾何拓?fù)錁?gòu)圖方法相對(duì)簡(jiǎn)單，而三維模型是一個(gè)多面體，并且金屬多晶材料是由多個(gè)不規(guī)則晶粒的空間構(gòu)成，因而其拓?fù)浣Y(jié)構(gòu)比二維情況要復(fù)雜得多，更為糟糕的是在三維模型下單元按晶粒分組和拓?fù)鋽?shù)據(jù)處理在MATLAB中均沒有函數(shù)直接實(shí)現(xiàn)。為此，本文緊緊圍繞Voronoi圖是由一組由連接兩鄰點(diǎn)直線的垂直平分線組成的連續(xù)多邊形，并根據(jù)ABAQUS/CAE中劃分網(wǎng)格后實(shí)體模型的單元編號(hào)空間排布特點(diǎn)，在MATLAB中編程依次求解每個(gè)單元的中心點(diǎn)坐標(biāo)到每個(gè)晶核的距離，并將該單元?dú)w屬到距離最短的晶核所在的晶胞內(nèi)，從而得到每個(gè)晶粒的完整拓?fù)湫畔?，接著將得到的每一個(gè)晶粒編號(hào)及對(duì)應(yīng)的所有單元編號(hào)，以set集合的形式寫入ABAQUS的INP文件中的Part部分，最終實(shí)現(xiàn)整個(gè)三維Voronoi模型的建立。

三維Voronoi多晶模型完整建模的程序流程圖如圖3所示，具體實(shí)施步驟如下：

①先在特定空間剖分得到晶核坐標(biāo)的隨機(jī)數(shù)，利用MATLAB中MPT工具箱提供的 mpt_voronoi命令，生成所需大小和晶粒數(shù)的Voronoi圖基本信息并畫出三維模型圖。

②結(jié)合單元編號(hào)區(qū)域排布特點(diǎn)，尤其是單元區(qū)域分布及單元編號(hào)順序，編程直接求取每一個(gè)單元的中心點(diǎn)坐標(biāo)，單元是通過在ABAQUS/CAE中建立的一個(gè)與三維模型圖相同尺寸的3D/Deformable/Shell模型劃分網(wǎng)格后所得，單元具體信息可以在生成的INP文件中查看。

③根據(jù)ABAQUS/CAE中劃分網(wǎng)格后實(shí)體模型的單元空間排布特點(diǎn)，依次求取每個(gè)單元的中心點(diǎn)坐標(biāo)到每個(gè)晶核的距離，并將該單元?dú)w屬到距離最短的晶核所在的晶胞內(nèi)。

④將上一步得到的每一個(gè)晶粒編號(hào)及其對(duì)應(yīng)的所有單元編號(hào)，最終以多級(jí)列表形式存儲(chǔ)到txt文件中，以供后續(xù)修改INP文件的調(diào)用。

通過以上MATLAB編程，本文僅得到了三維Voronoi圖完整的拓?fù)湫畔ⅰＯ旅孢€需要在ABAQUS的INP文件中Part部分編寫多個(gè)set集合，每個(gè)set集合為一個(gè)晶粒，并將txt文件中對(duì)應(yīng)的單元編號(hào)寫入相應(yīng)的set集合實(shí)現(xiàn)建模。利用INP文件進(jìn)行三維模型建立的過程如下：

①在ABAQUS/CAE中建立一個(gè)3D/Deformable/Shell，裝配，劃分網(wǎng)格單元(四邊形網(wǎng)格單元)，并在Job模塊中生成INP文件。

②將上一步寫好的set集合添加到建立的INP文件中的Part部分。在ABAQUS中打開修改后的INP文件，并在Color Code Dialog中設(shè)置為Color code by：Sets。

圖3　三維Voronoi多晶模型建立的流程圖

2案例分析

鑒于工程中使用最多的是金屬材料，而金屬材料又大多是晶粒的集合，即多晶體；并且多晶體材料的微觀組織結(jié)構(gòu)復(fù)雜，因此在建立Voronoi多晶有限元模型時(shí)，很難完全真實(shí)地反映出晶粒的各個(gè)方面的特征，必須對(duì)其進(jìn)行合理假設(shè)[12-13]。本文考慮到多晶材料在晶粒大小、形狀上均具有很大隨機(jī)性的特點(diǎn)，基于Voronoi圖原理通過MATLAB-MPT工具箱結(jié)合MATLAB編程，分別得到帶有邊界限制的完整的二維、三維Voronoi圖幾何拓?fù)湫畔?，并以建立集合的方式將處理好的晶粒編?hào)及其包含的所有單元信息寫入ABAQUS的INP文件中實(shí)現(xiàn)模型的建立。該模型不僅能夠十分形象真實(shí)地反映實(shí)際多晶體的構(gòu)成形態(tài)，而且還為后續(xù)的材料參數(shù)、取向的賦予、實(shí)施不同類型的加載等[14]以編寫INP文件的方式實(shí)現(xiàn)及其細(xì)觀尺度下金屬變形與內(nèi)部晶粒演變的關(guān)系分析提供前期的數(shù)據(jù)準(zhǔn)備。下面分別以二維和三維多晶體有限元模型構(gòu)建為例對(duì)本文所提出的方法進(jìn)行實(shí)現(xiàn)與驗(yàn)證。

2.1建立二維多晶體有限元模型

為了實(shí)現(xiàn)多晶體晶粒呈不同形態(tài)與尺寸大小的二維Voronoi圖，根據(jù)圖2所示的二維Voronoi多晶模型建立的流程圖，本文首先在MATLAB中隨機(jī)生成任意分布類型的種子點(diǎn)，再根據(jù)隨機(jī)種子點(diǎn)生成二維Voronoi圖。MATLAB中的二維Voronoi圖如圖4所示。

通過MATLAB編程雖然得到了二維Voronoi圖完整的拓?fù)湫畔?，并?shí)現(xiàn)了二維Voronoi圖，但它僅僅是帶有二維Voronoi圖完整拓?fù)湫畔⒌亩嗑w示意圖，并不是一個(gè)帶有完整邊界限制的二維Voronoi多晶體有限元模型。同樣，根據(jù)圖2所示流程圖的INP文件部分的分析流程，本文以建立的一個(gè)包含10個(gè)晶粒的二維多晶體有限元模型為例，其多晶體晶粒的相關(guān)參數(shù)與圖4 所示在MATLAB中生成的二維Voronoi圖尺寸及晶粒位置、個(gè)數(shù)均相同，并通過編寫ABAQUS的INP文件的形式進(jìn)行實(shí)現(xiàn)，最終結(jié)果如圖5所示。

圖4MATLAB中生成的二維Voronoi多晶體示意圖

Fig.4Two dimensional Voronoi polycrystalline

model generated in MATLAB

圖5ABAQUS中生成的二維Voronoi多晶體模型

Fig.5Two dimensional Voronoi multi crystal

model generated in ABAQUS

2.2建立三維多晶體有限元模型

根據(jù)圖3所示的三維Voronoi多晶模型建立的流程圖，本文首先在MATLAB中隨機(jī)生成三維空間任意分布類型的種子點(diǎn)，再根據(jù)隨機(jī)種子點(diǎn)生成的三維Voronoi圖。MATLAB中的三維Voronoi圖如圖6所示。

同樣，通過MATLAB編程得到的三維Voronoi圖也僅僅是帶有三維Voronoi圖完整拓?fù)湫畔⒌亩嗑w示意圖，并不是一個(gè)帶有完整邊界限制的三維Voronoi多晶體有限元模型，因此本文還需要再根據(jù)圖3所示流程圖中的INP文件部分的分析流程進(jìn)行三維Voronoi多晶體有限元模型構(gòu)建。本文以建立的一個(gè)包含10個(gè)晶粒的三維多晶體有限元模型為例，其多晶體晶粒的相關(guān)參數(shù)與圖6所示在MATLAB中生成的三維Voronoi圖尺寸及晶粒位置、個(gè)數(shù)均相同，并通過編寫ABAQUS的INP文件的形式進(jìn)行實(shí)現(xiàn)，最終結(jié)果如圖7所示。

圖6MATLAB中生成的三維Voronoi多晶體示意圖

Fig.6The three dimensional Voronoi

generated in MATLAB

圖7ABAQUS中生成的三維Voronoi多晶體模型

Fig.7Three dimensional Voronoi multi

crystal model generated in ABAQUS

3驗(yàn)證與討論

為了驗(yàn)證本文建模方法的正確性，本文建立了與文獻(xiàn)[9]相同晶粒數(shù)100個(gè)、相同網(wǎng)格類型C3D8I，相同基本尺寸1000 mm×400 mm×400 mm的Voronoi多晶有限元模型，并采用相同的本構(gòu)關(guān)系、相同的約束與相同的加載方式如圖8所示，還采用了相同的文獻(xiàn)[15]中所提供的材料參數(shù)，進(jìn)行了10%拉伸量的拉伸有限元模擬，結(jié)果如圖9與圖10所示。

圖9為進(jìn)行了拉伸有限元模擬之后的變形圖，其表面出現(xiàn)了較為明顯凹凸現(xiàn)象，這是由于晶體塑性本構(gòu)關(guān)系中考慮了晶體的隨機(jī)取向而導(dǎo)致的材料塑性變形不均勻性。

圖8　ABAQUS中生成的100個(gè)晶胞的Voronoi多晶模型

圖9　ABAQUS中有限元模擬后的Voronoi多晶模型變形圖

圖10　相同拉伸量下不同模型單元應(yīng)力—應(yīng)變曲線對(duì)比

圖10為不同模型單元應(yīng)力—應(yīng)變曲線對(duì)比圖。其中簡(jiǎn)化模型一的應(yīng)力—應(yīng)變曲線與精確模型的應(yīng)力—應(yīng)變曲線均出自于文獻(xiàn)[9]，而簡(jiǎn)化模型二的應(yīng)力—應(yīng)變曲線即為本文建模方法所建立的如圖8所示模型，并在相同拉伸量下進(jìn)行有限元模擬所得應(yīng)力—應(yīng)變曲線。從圖10中可以看出，本文所建立的Voronoi多晶簡(jiǎn)化模型整體應(yīng)力水平居于文獻(xiàn)[9]中簡(jiǎn)化模型一與精確模型之間，在彈性區(qū)域與精確模型大部分重合但屈服點(diǎn)較高，并且拉伸曲線整體偏向精確模型。為了進(jìn)一步深入了解本文拉伸曲線的偏向程度，本文分別在橫坐標(biāo)上取0.002、0.002 84、0.008、0.016、0.032、0.044、0.05與0.06等最大對(duì)數(shù)應(yīng)變所對(duì)應(yīng)的Mises應(yīng)力值來做偏差分析。本文假設(shè)以精確模型為參考標(biāo)準(zhǔn)，以簡(jiǎn)化模型與精確模型的差值進(jìn)行偏差分析來評(píng)估本文建模方法的優(yōu)劣。

表1　相同應(yīng)變下不同模型單元Mises應(yīng)力數(shù)據(jù)對(duì)比

從表1中可以得到本文所建立的Voronoi多晶簡(jiǎn)化模型的應(yīng)力水平與精確模型的偏差基本分布在20 MPa左右，最小偏差達(dá)為12.91 MPa，最大偏差僅為25.47 MPa，平均偏差值為18.3 MPa；而文獻(xiàn)[9]的簡(jiǎn)化模型一與精確模型的偏差基本分布在50 MPa左右，最小偏差達(dá)到27.54 MPa(比本文所建立的簡(jiǎn)化模型最大偏差25.47 MPa還大0.07 MPa)，同時(shí)其最大偏差達(dá)到60.13 MPa，平均偏差值為50.27 MPa。因此，說明本文建立的Voronoi多晶有限元模型更接近于文獻(xiàn)[9]的精確模型，可以用來作為金屬微觀結(jié)構(gòu)有限元分析的模型。

4結(jié)語(yǔ)

①本文提出了一種基于Voronoi圖并結(jié)合單元編號(hào)區(qū)域排布特點(diǎn)的方法，該方法可以依據(jù)單元編號(hào)區(qū)域排布特點(diǎn)直接得到單元編號(hào)，且更加容易實(shí)現(xiàn)依次求取單元形心坐標(biāo)。

②二維與三維Voronoi多晶體有限元模型在判斷單元?dú)w屬時(shí)，均采用單元距離最近的晶核所在晶粒內(nèi)的原則，可以通過改變隨機(jī)數(shù)所在的區(qū)域大小，靈活地改變Voronoi多晶體有限元模型的大小。

③一旦種子點(diǎn)的隨機(jī)坐標(biāo)生成無論其順序怎樣排列，最終得到的Voronoi多晶體有限元模型均完全相同；并可以通過改變隨機(jī)數(shù)的個(gè)數(shù)及分布類型，靈活地改變Voronoi多晶體有限元模型晶粒的數(shù)目及分布形態(tài)，一般晶粒數(shù)目大于50個(gè)，能夠較為真實(shí)描述材料的微觀組織結(jié)構(gòu)。

④在相同條件下通過有限元數(shù)值分析與文獻(xiàn)[9]的單調(diào)拉伸曲線對(duì)比，結(jié)果顯示本文所建立的Voronoi多晶簡(jiǎn)化模型整體應(yīng)力水平居于文獻(xiàn)[9]中簡(jiǎn)化模型一與精確模型之間，在彈性區(qū)域與精確模型大部分重合但屈服點(diǎn)比較高，并且拉伸曲線整體偏向精確模型，其平均偏差值僅為18.3 MPa，而文獻(xiàn)[9]中簡(jiǎn)化模型一的平均偏差值卻高達(dá)50.27 MPa，甚至連其最小偏差也比本文所建立的簡(jiǎn)化模型最大偏差還要高0.07 MPa。從而說明本文所提出的建模方法有較高的精度與可靠度。

參考文獻(xiàn)：

[1]趙昊,丁淑蘭.ABAQUS二次開發(fā)在金屬微觀組織模擬中的應(yīng)用[J]. 熱加工工藝,2014,43(24):87-90，97.

[2]王克廷,陳忠家,郭煜澤,等.溫度對(duì)合金塑性變形行為影響的蒙特卡羅模擬[J]. 有色金屬加工, 2014,43(1):9-11，40.

[3]陳飛,崔振山,董定乾.微觀組織演變?cè)詣?dòng)機(jī)模擬研究進(jìn)展[J]. 機(jī)械工程學(xué)報(bào),2015,51(4):30-39.

[4]聞瑤,薛克敏,李萍,等.TA15鈦合金高溫變形多晶體塑性有限元模擬[J]. 塑性工程學(xué)報(bào), 2014,21(6):86-90，101.

[5]萬朔,何力軍,張偉,等.多晶鈹微觀彈性失配行為的數(shù)值分析[J]. 稀有金屬, 2015,39(8):764-768.

[6]司良英, 鄧關(guān)宇 ,呂程,等.基于Voronoi圖的晶體塑性有限元多晶幾何建模[J]. 材料與冶金學(xué)報(bào), 2009,8(3):193-197，216.

[7]樊黎霞, 趙軻, 董雪花.身管徑向精密鍛造的塑性應(yīng)變分析與鍛造比研究[J]. 精密成形工程, 2014, 1(1):1-8.

[8]張豐果, 董湘懷, 章海明,等.微鐓粗過程晶體塑性模型研究[J]. 塑性工程學(xué)報(bào), 2011, 3(3):5-8.

[9]汪凱.多晶體材料加工的細(xì)觀塑性有限元模擬[D]. 昆明：昆明理工大學(xué), 2013.

[10]伊興華.面心立方多晶體塑性變形有限元模擬[D]. 哈爾濱：哈爾濱工業(yè)大學(xué), 2007.

[11]謝韶，唐斌，韓逢博，等.材料微結(jié)構(gòu)三維建模及其晶體塑性有限元模擬[J]. 塑性工程學(xué)報(bào)， 2014,21(1):65-70.

[12]QUEY R, DAWSON P R, BARBE F.Large-scale 3D random polycrystals for the finite element method: Generation, meshing andremeshing[J]. Computer Methods in Applied Mechanics & Engineering, 2011, 200(17-20):1729-1745.

[13]FRITZEN F, B?HLKE T, SCHNACK E.Periodic three-dimensional mesh generation for crystalline aggregates based on Voronoi Tessellations[J]. Computational Mechanics, 2009, 43(5):701-713.

[14]胡桂娟,張克實(shí),莫智莉.45號(hào)鋼后繼屈服與塑性流動(dòng)試驗(yàn)的Chaboche模型分析[J]. 廣西大學(xué)學(xué)報(bào)(自然科學(xué)版),2014,39(1):171-179.

[15]HUANG Y G.A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program[D]. Cambridge :Harvard University,1991.

(責(zé)任編輯梁健)

A method of polycrystal finite element modeling based on Voronoi diagram

ZHENG Zhan-guang1,2, WANG Zhao-liang1, FENG Qiang1, YUAN Shuai1, WANG Jia-xiang1

(1.College of Mechanical Engineering ,Guangxi University, Nanning 530004,China；2. Guangxi Key Laboratory of Manufacturing System & Advanced Manufacturing Technology,Nanning 530004,China)

Abstract:The mesoscopic polycrystal finite element (PFE) model is the premise and foundation to study the local plastic deformation of polycrystalline materials. In order to construct a reliable microstructure model of materials flexibly, a modeling method had been proposed on the basis of Voronoi diagram and regional distribution charact-eristics of units. The areal coordinates of the units could be calculated in sequence directly through the units’ num-bers by the method. Firstly, the random points and the Voronoi information in a specific plane or domain were generated. Then, the coordinates of the centers were directly calculated according to the distribution characteristics of the units’ numbers. Secondly, the nearest element was attributed to the nearest grain, the grains’ numbers and the units’ numbers were added to the part of the ABAQUS INP file by sets. Then, the PFE model was built. Lastly, taking 2D and 3D PFE models of 10 grains as examples, the stress analyses of the PFE models were compared with the results of the literature to verify the modeling method. The results reveal that the modeling method can not only directly get the units numbers and the areal coordinates of the units, but also reduces the data quantity on processing the areal coordinates of the units and the difficulties in determining the ownership of units. With the comparison, it is found that the accuracy of the model is closer to the exact model in the literature. The maximum deviation is only 25.47 MPa between the proposed model and the exact model, and the minimum deviation of the proposed model is lower than the simplified model in the literature by 0.07 MPa. It is indicated that the method can provide some references for building the Voronoi finite element model of polycrystalline materials quickly.

Key words:vonronoi diagram; polycrystal modeling; crystal plasticity; finite element method

中圖分類號(hào)：TP391.41

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1001-7445(2016)02-0460-10

doi:10.13624/j.cnki.issn.1001-7445.2016.0460

通訊作者：鄭戰(zhàn)光(1975—), 男,江西臨川人，廣西大學(xué)副教授,博士；E-mail：zhenglight@126.com。

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(5146002)；廣西自然科學(xué)基金項(xiàng)目(2012GXNSFBA053145)；廣西制造系統(tǒng)與先進(jìn)制造技術(shù)重點(diǎn)實(shí)驗(yàn)室項(xiàng)目(14-045-15S04)；玉林市科技攻關(guān)項(xiàng)目(玉市科攻1421001)

收稿日期：2015-11-05；

修訂日期：2015-12-17

引文格式：鄭戰(zhàn)光，汪兆亮，馮強(qiáng),等.一種基于Voronoi圖的多晶體有限元建模方法[J].廣西大學(xué)學(xué)報(bào)(自然科學(xué)版)，2016，41(2)：460-469.