張禹，劉如鐵，陳潔，熊翔，李浩，王鑄博

基于Voronoi分形劃分的雙孔徑多孔材料孔道結(jié)構(gòu)參數(shù)化設(shè)計(jì)

張禹，劉如鐵，陳潔，熊翔，李浩，王鑄博

(中南大學(xué) 粉末冶金國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083)

本文提出了一種適用于模擬雙孔徑多孔介質(zhì)的孔道結(jié)構(gòu)模型，該模型可獨(dú)立形成空間均勻分布的大、小尺寸的兩種孔徑，并在此基礎(chǔ)上開(kāi)發(fā)可以生成參數(shù)化孔道模型圖示及構(gòu)成體素坐標(biāo)的面向過(guò)程程序。通過(guò)劃分函數(shù)的迭代，從對(duì)應(yīng)幾何方法的自相似性引出此模型的一般孔隙度、分形維度、間隙度和比表面積方面的討論。最后通過(guò)面向?qū)ο缶幊?，選擇工業(yè)三維建模軟件開(kāi)發(fā)一種可以生成較大規(guī)模幾何模型文件的技術(shù)路線及滲流與電化學(xué)仿真應(yīng)用方法。

Voronoi劃分；分形；二元孔隙度；孔道結(jié)構(gòu)；多孔介質(zhì)；多孔電極

多孔材料中孔隙的空間分布具有隨機(jī)性，孔隙連通方式和孔隙形貌構(gòu)成了材料的孔道結(jié)構(gòu)，所以多孔材料的孔道結(jié)構(gòu)具有隨機(jī)和復(fù)雜的特征。多孔材料比表面積主要由材料內(nèi)部表面貢獻(xiàn)，孔道結(jié)構(gòu)連通了材料部分內(nèi)外表面，外部流體可以進(jìn)入孔隙中充分接觸材料，因此多孔材料在石油催化、吸附分離、電化學(xué)、傳熱與傳質(zhì)方面有著廣泛應(yīng)用[1?4]。設(shè)計(jì)多孔材料的孔道結(jié)構(gòu)有利于研究和生產(chǎn)過(guò)程中材料參數(shù)識(shí)別，材料性能預(yù)測(cè)以及指導(dǎo)新材料的研究[5?6]。

隨著多孔材料的不斷發(fā)展，不同應(yīng)用領(lǐng)域?qū)τ诙嗫撞牧系目椎澜Y(jié)構(gòu)特征提出不同的要求。例如，在環(huán)路熱管的應(yīng)用中，多孔毛細(xì)芯所產(chǎn)生的毛細(xì)壓力是熱管內(nèi)工質(zhì)流動(dòng)循環(huán)的關(guān)鍵驅(qū)動(dòng)力。其中小孔貢獻(xiàn)毛細(xì)壓力，大孔貢獻(xiàn)高孔隙度以保證良好的滲透性[7]；在多孔氧化硅骨架鎳催化裂解甲烷應(yīng)用中，既需要引入大孔徑結(jié)構(gòu)以避免積碳堵塞和促進(jìn)反應(yīng)物流動(dòng)，也需要引入大量小孔結(jié)構(gòu)提升催化劑比表面積[8]；在空氣凈化濾芯應(yīng)用中，家用空氣凈化設(shè)備出于節(jié)能和過(guò)濾效能考慮，要求濾芯同時(shí)具有小孔(典型直徑小于2.5 μm)保證高比表面積，從而有效吸附通過(guò)空氣中的懸浮顆粒，同時(shí)利用大孔結(jié)構(gòu)的高連通性降低風(fēng)阻[9]。綜上，對(duì)多孔材料進(jìn)行微觀孔道結(jié)構(gòu)設(shè)計(jì)，是制備滿足應(yīng)用要求的高性能多孔材料的前提和基礎(chǔ)。

當(dāng)多孔材料中同時(shí)存在分別集中在不同尺寸范圍的兩類孔徑，且兩類孔徑尺寸相差10倍以上時(shí)，稱此類多孔材料的孔道結(jié)構(gòu)為雙孔徑結(jié)構(gòu)。雙孔徑結(jié)構(gòu)的產(chǎn)生歸因于其制備過(guò)程中不同的造孔機(jī)理。例如以聚氨酯發(fā)泡為模板的多孔鎳材料，添加造孔劑燒結(jié)制備的鎳基或銅基粉末冶金材料，添加造孔劑燒結(jié)金屬氈等，其孔徑分布為典型的雙峰結(jié)構(gòu)[10?11]，大孔徑由造孔劑形成，小孔徑由粉末材料在燒結(jié)過(guò)程中的多種遷移機(jī)理共同作用形成。設(shè)計(jì)雙孔徑結(jié)構(gòu)的目的在于利用其不同大小孔隙的相互匹配，實(shí)現(xiàn)高比表面積與連通性、滲透性及其他多種特定應(yīng)用性能的組合，而如何設(shè)計(jì)孔道結(jié)構(gòu)是雙孔徑材料制備的前提。

目前研究者們發(fā)展了多種獲得孔道結(jié)構(gòu)的方法：隨機(jī)幾何方法、模擬沉積方法、統(tǒng)計(jì)學(xué)習(xí)方法，以及micro-CT掃描等[12]。隨機(jī)幾何方法使用二維或三維隨機(jī)過(guò)程算法和輸入?yún)?shù)對(duì)孔道結(jié)構(gòu)進(jìn)行建模生成，輸出幾何模型的統(tǒng)計(jì)不變量。相比其他方法，隨機(jī)幾何方法具有設(shè)備資源要求較低的優(yōu)勢(shì)[13]。Voronoi圖是一種常用的隨機(jī)幾何方法，也是自然界中的常見(jiàn)現(xiàn)象[14]，它的概念于1907年由VORONOI完善，已經(jīng)廣泛用于生物學(xué)、地理學(xué)、材料科學(xué)等方向的研究[15?17]。目前研究應(yīng)用Voronoi圖建立雙孔徑孔道結(jié)構(gòu)方法可總結(jié)為兩類：以PVT(poisson-voronoi tessellation)為代表的調(diào)節(jié)Voronoi生成點(diǎn)空間分布方法和以LVT(laguerre-voronoi tessellation)為代表的度量空間算法。本研究分析以上兩種算法在構(gòu)建雙孔徑結(jié)構(gòu)的孔隙空間分布和數(shù)量分布的不足，從一種二維的Voronoi分形劃分的概念獲得啟發(fā)[18]，提出一種以三維Voronoi劃分迭代地構(gòu)造雙孔徑材料孔道結(jié)構(gòu)的方法，如圖1。

考慮到兩次劃分迭代產(chǎn)生的自相似特性，本研究分析推導(dǎo)了估算迭代PVT的孔隙度、分形維度、比表面積和間隙度的公式。其次，編寫(xiě)輸出孔道結(jié)構(gòu)的迭代算法，并且借助常用的數(shù)學(xué)建模及工業(yè)軟件，開(kāi)發(fā)一種生成模型文件的技術(shù)路線。最后，將程序輸出的模型作為計(jì)算域，展示了本研究提出的模型在滲流仿真和電化學(xué)仿真應(yīng)用的可行性。

1 Voronoi劃分與孔道結(jié)構(gòu)

圖1 雙孔徑模型生成順序圖

(a) Cylindrical region; (b) Poisson points in cylindrical region; (c) Voronoi tessellation based on poisson points; (d) Zoom-in on a thin square for observation; (e) Second order cylindrical region; (f) Voronoi tessellation in each second order cylindrical region

設(shè)置PVT的泊松強(qiáng)度為空間位置的函數(shù)，形成生成點(diǎn)強(qiáng)度的各向異性分布，從而可以調(diào)整生成點(diǎn)的密度，造成PVT單胞的聚集和離散，形成PVT二元孔徑分布。然而這種結(jié)構(gòu)孔隙空間分布繼承了生成點(diǎn)強(qiáng)度的空間各向異性，因此雙孔徑結(jié)構(gòu)模型呈各向異性，不符合雙孔徑孔道結(jié)構(gòu)兩種孔隙分別空間均勻分布的特征。

對(duì)比孔隙尺寸集中單峰分布的氧化鋁泡沫，LVT模型在單胞體積、單胞面的數(shù)目、單胞棱邊數(shù)目上都非常接近實(shí)際材料，但是LVT目前無(wú)法表征雙孔徑多孔材料的孔道結(jié)構(gòu)，原因是LVT的生成點(diǎn)來(lái)自分子動(dòng)力學(xué)模型混合的堆積球體，球體在擾動(dòng)的邊界條件下充分混合，形成致密且隨機(jī)的堆積。如使用LVT方法進(jìn)行雙孔徑模型建模，需要同時(shí)引入尺寸對(duì)應(yīng)雙峰分布的球體，然而尺寸差異較大的顆粒難以在振動(dòng)力場(chǎng)下混合均勻[23]，造成LVT無(wú)法形成兩種尺寸獨(dú)立分布的隨機(jī)球體填充。

借助以往研究者對(duì)PVT模型參數(shù)化形成特定孔隙度的成果[24]和一種二維Voronoi迭代隨機(jī)幾何圖形的算法設(shè)計(jì)[18]，考慮到雙孔徑材料大孔和小孔的形成機(jī)理的獨(dú)立性，首先將總孔隙度total分配為由大孔形成的部分1和小孔形成的部分2，對(duì)應(yīng)實(shí)際測(cè)量雙孔徑材料塊體的孔隙度和壓汞法獲得的孔徑分布，如式(2)：

WEJRZANOWSKI等[21]總結(jié)了多種驗(yàn)證多孔模型的孔隙度公式，LEI等[25]從規(guī)則晶格劃分出發(fā)開(kāi)發(fā)了孔隙度公式，并且實(shí)現(xiàn)了孔道結(jié)構(gòu)的梯度變化。本文利用LEI等開(kāi)發(fā)的式(4)，調(diào)整隨機(jī)點(diǎn)的數(shù)量和棱邊纖維設(shè)置的半徑，獲得多孔模型孔隙度，其中C為高階修正項(xiàng)，為正方體空間的邊長(zhǎng)。

2 算法

2.1 二階Voronoi劃分

為了減少式(4)的變量數(shù)目，考慮到Voronoi劃分適用的空間不限于其形狀，本文認(rèn)為式(4)的邊長(zhǎng)參數(shù)可以推廣定義為與空間體積相等的正方體邊長(zhǎng)，記的體積為V，可得式(5)。

將式(5)和式(6)代入式(4)并忽略高階修正項(xiàng)，得總孔隙度，如式(7)：

記第一階劃分的泊松點(diǎn)集強(qiáng)度參數(shù)為λ1，設(shè)置纖維半徑為1，第二階劃分的泊松點(diǎn)集強(qiáng)度參數(shù)為2，設(shè)置纖維半徑為2。代入式(2)、式(3)、式(8)，得到式(9)，其中=1,2。

研究實(shí)驗(yàn)中通常可以獲得兩種孔徑尺寸d的數(shù)值[8, 10]，本文以每階劃分形成的單胞體積v代表單個(gè)孔隙體積，孔徑尺寸d為等效球形孔隙直徑，如式(10)。而劃分產(chǎn)生在空間Ω中的平均孔隙總數(shù)與生成點(diǎn)數(shù)n相等，結(jié)合式(6)，可得式(11)。

由式(10)和式(11)得式(12)，即每階劃分產(chǎn)生的孔隙度P、孔徑d與輸入泊松強(qiáng)度λ的關(guān)系，結(jié)合式(9)解出建立二階PVT模型需要的所有參數(shù)。

2.2 模型圖示

UHLMANN[26]指出目前Voronoi劃分的學(xué)術(shù)研究通常使用Matlab或“VORO++”的函數(shù)庫(kù)，兩種方法都提供了運(yùn)行時(shí)間與生成點(diǎn)規(guī)模成線性關(guān)系的快速算法。本研究第二次劃分的空間位置由第一次劃分的棱邊決定，需要對(duì)大量坐標(biāo)進(jìn)行平移、投影和旋轉(zhuǎn)，因而選用Matlab作為向量編程的工具。表1所列為多階劃分算法主函數(shù)的偽代碼，圖2所示為分別調(diào)整大孔與小孔泊松點(diǎn)集強(qiáng)度生成的雙孔徑模型，可見(jiàn)本文算法具有獨(dú)立調(diào)整兩種孔徑分布的功能。

2.3 分析

2.3.1 孔隙度

表1 生成N+1階劃分算法的偽代碼

表2 投影法處理圓柱區(qū)域邊界算法的偽代碼

2.3.2 分形維度

由式(16)、(17)可知式(13)中的參數(shù)為平均化圓柱空間的半徑與長(zhǎng)度的比值。

在初始圓柱長(zhǎng)度為單位長(zhǎng)度時(shí)，代入幾何相似式(9)和式(14)，式(18)可化簡(jiǎn)為式(19)和式(20)：

圖2 相對(duì)厚度為0.5不同泊松強(qiáng)度下形成的大小孔模型雙孔隙度圖示

(a) 30, 300; (b) 30, 600; (c) 60, 300; (d) 60, 600

將式(2)，式(3)推廣至多峰分布，可得孔隙度組成的級(jí)數(shù)形式，如式(23)。將式(23)代入式(22)得到(24)：

WEI等[28]研究了分形幾何方法推導(dǎo)的毛細(xì)孔徑分布，如式(25)。式(24)與式(25)存在差別的原因是WEI的分形方法為迭代地劃分孔隙，本研究的方法為迭代地劃分實(shí)體，兩種方法關(guān)于孔隙相-實(shí)體相對(duì)稱。

從式(25)獲得按冪次分布的孔徑密度分布函數(shù)，將此函數(shù)與泊肅葉流動(dòng)壓降進(jìn)行積分，得到預(yù)測(cè)的Kozeny-Carman數(shù)值與乳液滲流實(shí)驗(yàn)結(jié)果相符[28]。式(24)證明Voronoi分形劃分方法可以獲得冪次分布的孔徑，所以式(18)～(24)分析得到的分形維度可用于雙孔徑模型的流體分析。

由于泊肅葉流動(dòng)將復(fù)雜的孔道近似為平行的圓柱流道，所以研究者引入屈曲度，孔道分形維度等幾何不變量改善精度，這些幾何不變量代表相應(yīng)測(cè)量的數(shù)值與分辨率的變化規(guī)律，一般以Hausdorff“數(shù)格子”的方法獲得。在向量編程環(huán)境中將待處理的模型所在空間網(wǎng)格化，然后判斷每一格點(diǎn)與模型的相交狀態(tài)，例如使用Matlab內(nèi)建函數(shù)find()，可以提取待計(jì)數(shù)的格點(diǎn)。

2.3.3 間隙度

間隙度是衡量聚集和不均勻性的物理量[29]，XIA等[30]使用算法生成二維孔隙圖形并以LBM(lattice Boltzmann method，格子玻爾茲曼方法)分析了其滲透系數(shù)，結(jié)果表明對(duì)于固定的孔隙度，滲透系數(shù)以正冪次規(guī)律與間隙度正相關(guān)，與分形維度負(fù)相關(guān)[30]。此結(jié)論也印證了雙孔徑材料引入較大孔徑增加傳質(zhì)能力的機(jī)理，即在介觀尺度下改變材料的間隙度。

計(jì)算間隙度的數(shù)值需首先設(shè)置位于模型內(nèi)邊長(zhǎng)為的立方體窗格空間Box()，滑動(dòng)窗格并獲取窗格空間的孔隙度Box，一般將滑動(dòng)窗格網(wǎng)格化并設(shè)置格點(diǎn)采樣閾值，所以Box為離散型數(shù)值，記滑動(dòng)窗格次統(tǒng)計(jì)Box對(duì)應(yīng)的頻率為(Box,)，則間隙度(Box,)由式(26)計(jì)算[29]。

2.3.4 比表面積

需注意的是，采用式(28)估算實(shí)際材料比表面積需另外考慮材料孔道內(nèi)表面粗糙度的貢獻(xiàn)。

3 模型實(shí)現(xiàn)及應(yīng)用

3.1 幾何模型

所有輸入體素的并集，對(duì)應(yīng)建模軟件的布爾求和運(yùn)算。雖然目前存在可以與數(shù)值分析軟件交互的三維建模工具，如COMSOL Multiphysics，Ansys APDL 等[5, 19?20, 24]，但是其交互的穩(wěn)定性尚待提高。當(dāng)繪制體素?cái)?shù)高于1 000時(shí)，在內(nèi)存有限的桌面級(jí)電腦上，無(wú)論是修改幾何文件容差還是分步存儲(chǔ)，建模進(jìn)程都難以保證平穩(wěn)運(yùn)行。本文選用Solidworks在.NET環(huán)境中與Excel進(jìn)行交互，Excel來(lái)自預(yù)先存儲(chǔ)的Matlab所輸出的坐標(biāo)矩陣，在逐條讀取坐標(biāo)的同時(shí)，.NET環(huán)境調(diào)用Solidworks程序進(jìn)程在相應(yīng)文件中形成體素，循環(huán)至所有體素創(chuàng)建完畢后保存文件，可有效保證交互的穩(wěn)定性。

將面向過(guò)程的矩陣運(yùn)算與面向?qū)ο蟮慕＿^(guò)程獨(dú)立編程，可顯著提升效率和程序穩(wěn)定性。實(shí)時(shí)呈現(xiàn)運(yùn)算圖示可預(yù)先發(fā)現(xiàn)建模的錯(cuò)誤，將坐標(biāo)數(shù)據(jù)預(yù)先存儲(chǔ)在散列表文件中可以節(jié)約運(yùn)行內(nèi)存，集中計(jì)算資源運(yùn)行體素搭建。

體素搭建是一個(gè)循環(huán)輸入體素并與既有體素布爾求和的過(guò)程。本文采用向Solidworks輸入臨時(shí)體(一種按參數(shù)形成的body類型)的方式輸入體素，再轉(zhuǎn)換臨時(shí)體為殼體(surface類型)，最后填充殼體獲得實(shí)心圓柱體，完成體素輸入。執(zhí)行填充時(shí)設(shè)置merge為T(mén)rue，即可在輸入體素的同時(shí)自動(dòng)布爾求和，上述調(diào)用函數(shù)均如圖3所示。

在2016版本后Solidworks不再限制多實(shí)體零件文件的實(shí)體數(shù)上限，此功能強(qiáng)力支持了多孔體的孔道結(jié)構(gòu)設(shè)計(jì)，并且.NET環(huán)境所支持的語(yǔ)言如VB.NET，C#，C++都可以將每個(gè)體素的建立加入Try-Catch語(yǔ)句，增加程序的健壯性，避免某次體素建立失敗導(dǎo)致程序崩潰。

3.2 仿真應(yīng)用

以有限元或有限體積方法為例，多孔介質(zhì)的滲流模擬通常需要耗費(fèi)巨大的計(jì)算資源，且計(jì)算域復(fù)雜的幾何形狀會(huì)導(dǎo)致極其復(fù)雜的網(wǎng)格劃分，又需要耗費(fèi)人力進(jìn)行網(wǎng)格調(diào)整和優(yōu)化。因此，LBM因其易于編程部署，適用并行計(jì)算，前后處理難度低，適合低速流動(dòng)分析的特點(diǎn)被廣泛應(yīng)用于微流動(dòng)仿真分析領(lǐng)域。本文介紹一種通過(guò)圖像處理方法將三維模型再次壓縮為二維圖形，可以達(dá)到繼續(xù)簡(jiǎn)化計(jì)算規(guī)模且能得出定性結(jié)論的目的。如圖4(a)，黑色區(qū)域?yàn)橥ㄟ^(guò)圖像處理采集三維薄片得到的多孔介質(zhì)實(shí)體部分，藍(lán)色區(qū)域?yàn)榱黧w的計(jì)算域，藍(lán)色的深淺代表LBM運(yùn)行4 000周期該處的無(wú)量綱流速，從而根據(jù)實(shí)際的雷諾數(shù)估計(jì)相應(yīng)多孔介質(zhì)的滲流性能?？梢?jiàn)本案例中，高流速集中在大孔區(qū)域的喉部，根據(jù)泊肅葉定律層流壓降正比于流速的估計(jì)，這與研究多孔介質(zhì)壓降的直觀經(jīng)驗(yàn)相符。

同時(shí)，壓縮三維模型得到的二維圖形同樣可以應(yīng)用于電化學(xué)仿真分析，在不考慮電解質(zhì)流動(dòng)和宏觀離子濃度變化時(shí)，使用有限元或有限體積法進(jìn)行電流密度計(jì)算相比流體力學(xué)仿真更加快速。如圖4(b)所示，紅色線代表電解質(zhì)中的電流密度線，可見(jiàn)電流密度分布取決于計(jì)算域的幾何形狀。

本研究開(kāi)發(fā)的多孔模型可應(yīng)用于三種簡(jiǎn)化假設(shè)的電化學(xué)分析模型：

1) 考慮溶液整體不存在活性離子濃度梯度，電解液在恒穩(wěn)電壓下表現(xiàn)為電阻。此應(yīng)用輸出上述電流密度線和電勢(shì)等勢(shì)面，可用于優(yōu)化電解池整體效率。

2) 考慮電極表面存在濃差極化和雙電層電容，結(jié)合本研究模型，可以進(jìn)行循環(huán)伏安或電化學(xué)阻抗譜等分析。例如可通過(guò)提高小孔生成強(qiáng)度，提升模型比表面積，影響雙電層電容數(shù)值。

圖3 一個(gè)3746體素的多實(shí)體SLDPRT型文件(a)及Solidworks內(nèi)建函數(shù)調(diào)用流程(b)

圖4 使用格子玻爾茲曼方法滲流分析(a)及電化學(xué)模擬電解質(zhì)中多孔電極產(chǎn)生的電流密度分布(b)

3)考慮電解液存在流場(chǎng)，且存在非活性離子支持遷移電流，結(jié)合LBM在同一幾何計(jì)算域的流速分布結(jié)論，可實(shí)現(xiàn)流體力學(xué)和活性離子物質(zhì)傳遞的耦合，從而將電化學(xué)反應(yīng)還原為濃度依賴動(dòng)力學(xué)，實(shí)現(xiàn)在流動(dòng)電解質(zhì)中多孔電極的電化學(xué)仿真，即完整描述Nernst-Plank物質(zhì)傳遞方程。

4 結(jié)論

1) 迭代PVT方法可實(shí)現(xiàn)雙孔徑模型的大孔和小孔的尺寸分布獨(dú)立，同時(shí)保持空間分布均勻。通過(guò)控制每次迭代輸入?yún)?shù)可輸出具有期望孔隙度、分型維度、比表面積和間隙度的模型。

2) 使用Solidworks API可構(gòu)建本研究模型的實(shí)體文件，且實(shí)體文件包含的體素?cái)?shù)具有較復(fù)雜的規(guī)模。使用Matlab圖形模塊可以處理本研究模型的二維 投影。

3) 應(yīng)用LBM方法可實(shí)現(xiàn)分析本研究模型的多孔介質(zhì)滲流仿真。在ComsolMultiphysics環(huán)境中利用本研究模型可進(jìn)行多孔電極表面電化學(xué)分析。后續(xù)研究可實(shí)現(xiàn)相同模型計(jì)算域在不同仿真環(huán)境和多種物理場(chǎng)條件下的耦合仿真。

[1] FOROOZESH J, KUMAR S. Nanoparticles behaviors in porous media: application to enhanced oil recovery[J]. Journal of Molecular Liquids, 2020: 113876.

[2] DIAO K K, ZHAO Y Y. Heat transfer performance of sintered Cu microchannelsproduced by a novel method[J]. International Journal of Heat and Mass Transfer, 2019, 139: 537?547.

[3] EGOROV V, O'DWYER C. Architected porous metals in electrochemical energy storage[J]. Current Opinion in Electrochemistry, 2020, 21: 201?208.

[4] WOOD S, HARRIS A T. Porous burners for lean-burn applications[J]. Progress in Energy & Combustion Science, 2008, 34(5): 667?684.

[5] NIE Z W, LIN Y Y, TONG Q B, et al. Numerical investigation of pressure drop and heat transfer through open cell foams with 3D Laguerre-Voronoi model[J]. International Journal of Heat & Mass Transfer, 2017, 42(2): 969?974.

[6] LOTFI N, FARAHANI T S, YAGHOUBINEZHAD Y, et al. Simulation and characterization of hydrogen evolution reaction on porous NiCu electrode using surface response methodology[J]. International Journal of Hydrogen Energy, 2019, 44(26): 13296?13309

[7] LIU R T, CHEN J, XIONG X. Influence of porogen type and copper powder morphology on property of sintering copper porous materials[J]. Journal of Central South University, 2018, 25(9): 2143?2149.

[8] TANGGARNJANAVALUKUL C, DONPHAI W, WITOON T, et al. Deactivation of nickel catalysts in methane cracking reaction: effect of bimodal meso–macropore structure of silica support[J]. Chemical Engineering Journal, 2015, 262:364?371.

[9] FU X W, LIU J J, DING C F, et al. Building bimodal Structures by a wettability difference-driven strategy for high- performance protein air-filters[J]. Journal of Hazardous Materials, 2021, 415(42): 125742.

[10] LI H, LIU R T, CHEN J, et al. Preparation of nickel porous materials by sintering nickel oxalate and sodium chloride after blending and reduction[J]. Journal of Materials Research and Technology, 2020, 9(3): 3149?3157.

[11] WANG Z B, CHEN J, LIU R T, et al. Preparation of an ultrafine nickel powder by solid-phase reduction with a NaCl separator agent[J]. Advanced Powder Technology, 2020, 31(8): 3433?3439.

[12] AL-KHARUSI A S, BLUNT M J. Network extraction from sandstone and carbonate pore space images[J]. Journal of Petroleum Science and Engineering, 2007, 56(4):219?231.

[13] SURI S. Handbooks in Operations Research and Management Science[M]. Amsterdam: IGI Global, 1995: 425?479.

[14] MEBATSION H K, VERBOVEN P, VERLINDEN B E, et al. Microscale modeling of fruit tissue using Voronoi tessellations[J]. Computers and Electronics in Agriculture, 2006, 52(1/2):36?48.

[15] HAN N, GUO R. Two new Voronoi cell finite element models for fracture simulation in porous material under inner pressure[J]. Engineering Fracture Mechanics, 2019,211:478?494.

[16] GATSONIS N A, SPIRKIN A. A three-dimensional electrostatic particle-in-cell methodology on unstructured Delaunay–Voronoi grids[J]. Journal of Computational Physics, 2009, 228(10): 3742?3761.

[17] MISTANI P, GUITTET A, POIGNARD C, et al. A parallel Voronoi-based approach for mesoscale simulations of cell aggregate electropermeabilization[J]. Journal of Computational Physics, 2019, 380: 48?64.

[18] SHIRRIFF K. Generating fractals from Voronoi diagrams[J]. Computers & Graphics, 1993, 17(2):165?167.

[19] XIAO F, YIN X L. Geometry models of porous media based on Voronoi tessellations and their porosity–permeability relations[J]. Computers & Mathematics with Applications, 2016, 72(2): 328? 348.

[20] CAO B W, WANG S L, WEI D, et al. Investigation of the filtration performance for fibrous media: coupling of a semi-analytical model with CFD on Voronoi-based microstructure-Science Direct[J]. Separation and Purification Technology, 2020, 251:117364.

[21] WEJRZANOWSKI T, SKIBINSKI J, SZUMBARSKI J, et al. Structure of foams modeled by Laguerre-Voronoi tessellations[J]. Computational Materials Science, 2013, 67: 216?221.

[22] NIE Z W, LIN Y Y, TONG Q B. Modeling structures of open cell foams[J]. Computational Materials Science, 2017, 131: 160? 169.

[23] KHALA M J, HARE C, WU C Y, et al. Density and size-induced mixing and segregation in the FT4 powder rheometer: an experimental and numerical investigation[J]. Powder Technology, 2021, 390: 126?142.

[24] LIU W, TANG G H, SHI Y. Apparent permeability study of rarefied gas transport properties through ultra-tight Voronoi porous media by discrete velocity method[J]. Journal of Natural Gas Science and Engineering, 2020, 74: 103100.

[25] LEI H Y, LI J R, XU Z J, et al. Parametric design of Voronoi-based lattice porous structures[J]. Materials & Design, 2020, 191:108607.

[26] UHLMANN M. Vorono tessellation analysis of sets of randomly placed finite-size spheres[J]. Physica A: Statistical Mechanics and its Applications, 2020, 555: 124618.

[27] HENDERSON N, BRêTTASJC, SACCO W F. A three- parameter Kozeny–Carman generalized equation for fractal porous media[J]. Chemical Engineering Science, 2010, 65(15): 4432?4442.

[28] WEI W, CAI J C, XIAO J F, et al. Kozeny-Carman constant of porous media: insights from fractal-capillary imbibition theory[J]. Fuel, 2018, 234: 1373?1379.

[29] PINTO E P, PIRES M A, MATOS R S, et al. Lacunarity exponent and Moran index: acomplementary methodology to analyze AFM images and its application to chitosan films[J]. Physica A: Statistical Mechanics and its Applications, 2021, 581: 126192.

[30] XIA Y X, WEI W, LIU Y, et al. A fractal-based approach to evaluate the effect of microstructure on the permeability of two-dimensional porous media[J]. Applied Geochemistry, 2021, 131: 105013.

Parametric design of bimodal porosity structure based on Voronoi polygons and fractals

ZHANG Yu, LIU Rutie, CHEN Jie, XIONG Xiang, LI Hao, WANG Zhubo

(State Key Laboratory of Powder Metallurgy, Central South University, Changsha 410083, China)

This research proposed a bimodal porosity structure type providing the ability to individually form evenly distributed pores of two separated sizes. A method of generating structures parametrically was provided, and a process-oriented program was developed outputting diagrams and coordinates. Iterations of tessellation module, causing geometrical self-resemblance, were further discussed from porosity, fractal dimension, lacunarity and specific surface area aspects. Technical routes were introduced for building 3D model files of large complexity in common industrial software’s object-oriented application programming interface and applications in permeability or electrochemistry simulations.

Voronoi tessellation; fractal; bimodal porosity; porous structure; porous media; porous electrode

10.19976/j.cnki.43-1448/TF.2022007

TB113

A

1673-0224(2022)03-257-10

國(guó)家國(guó)際科技合作專項(xiàng)資助項(xiàng)目(2015DFR50580)

2022?03?23；

2022?04?20

劉如鐵，教授，博士。電話：0731-88876566；E-mail: llrrtt@csu.edu.cn; 陳潔，副教授，博士。電話：0731-88876566；E-mail: chenjiecsu@163.com

(編輯 陳潔)