顧 政，鄭建軍，周欣竹

（浙江工業(yè)大學(xué)建筑工程學(xué)院，杭州 310014）

多孔材料普遍存在于自然界中，木材、生物的骨骼以及巖石等都是天然的多孔材料。隨著科學(xué)技術(shù)的發(fā)展，出現(xiàn)了越來(lái)越多用金屬、陶瓷和高聚合物等制成的人造多孔材料［1］。這些多孔材料已廣泛應(yīng)用于各個(gè)工程領(lǐng)域，它們不僅具有多種優(yōu)異的性能，而且制造工藝簡(jiǎn)單。因此，研究其力學(xué)性能具有突出的重要性。

由于多孔材料是由孔隙和固相所組成的復(fù)合體，孔結(jié)構(gòu)（如孔形狀、孔隙率和孔的連通性等）是影響宏觀彈性性能的主要因素，因而細(xì)觀結(jié)構(gòu)與宏觀力學(xué)性能之間的定量關(guān)系成為當(dāng)前國(guó)際工程界的前沿課題之一。Voigt和Reuss分別根據(jù)等應(yīng)變和等應(yīng)力假設(shè)給出了多晶體材料體積模量和剪切模量的近似解［2］，Gibson和Ashby［3］在單孔單元的基礎(chǔ)上建立模型，獲得了蜂窩多孔材料二維彈性參數(shù)，即Gibson－Ashby方程，Roberts［4］研究了開(kāi)孔和閉孔泡沫的彈性性能。隨著數(shù)值方法和計(jì)算機(jī)技術(shù)的日益發(fā)展，有限元法被廣泛用于多孔材料的力學(xué)性能分析［5］，對(duì)“代表性體積元”進(jìn)行數(shù)值求解，獲得宏觀力學(xué)性能。但對(duì)于復(fù)雜的多孔材料，不僅網(wǎng)格劃分極其困難，總剛度矩陣占據(jù)大量?jī)?nèi)存，而且總剛度方程求解花費(fèi)大量時(shí)間，以至于無(wú)法獲得滿意的數(shù)值解。為此，該文在前人工作的基礎(chǔ)上，應(yīng)用快速傅里葉變換法討論了多孔材料彈性模量的計(jì)算方法。

1 數(shù)值方法

1.1 多孔材料模擬

類(lèi)似于“代表性體積元”，可以從細(xì)觀結(jié)構(gòu)圖像的像素點(diǎn)上獲得結(jié)構(gòu)單元的代表性信息［6］。因而為了能較準(zhǔn)確地預(yù)測(cè)多孔材料的彈性模量，首先應(yīng)建立多孔材料模型，作為初步嘗試，文章僅考慮二維模型。多孔材料區(qū)別于普通密實(shí)固體材料的最顯著特點(diǎn)是具有孔隙，多孔材料建模時(shí)應(yīng)著重考慮孔隙分布特點(diǎn)以及孔隙率大小。大量的試驗(yàn)研究表明［7］，多孔材料中的孔隙率、孔徑分布、孔隙位置等均服從一定的統(tǒng)計(jì)分布規(guī)律，如隨機(jī)分布和正態(tài)分布等。因此，建立多孔材料幾何模型的關(guān)鍵是按照一定的概率分布確定孔隙大小和位置，在數(shù)學(xué)上可以通過(guò)各種變換或抽樣來(lái)實(shí)現(xiàn)。

為了便于計(jì)算，該文假設(shè)孔隙為圓形，而且大小相等、互不重疊。在模擬孔隙時(shí)，先選取一邊長(zhǎng)為L(zhǎng)的正方形，在分布第i個(gè)半徑為R的圓孔時(shí)，在正方形區(qū)域內(nèi)生成其圓心坐標(biāo)（xi，yi），如果第i個(gè)圓孔與前面已經(jīng)分布某一個(gè)或幾個(gè)圓孔重疊，則重新生成第i個(gè)圓孔的圓心坐標(biāo)；如果第i個(gè)圓孔不與前面已經(jīng)分布的（i－1）個(gè)圓孔重疊，那么繼續(xù)分布第i＋1個(gè)圓孔，直到達(dá)到給定的孔隙率為止。作為一個(gè)算例，設(shè)L＝100mm，R＝4mm，孔隙率C分別為0.1、0.3和0.5，所獲得的孔分布如圖1所示。

1.2 基本方程

首先在模擬區(qū)域內(nèi)等距離選取像素點(diǎn)，根據(jù)快速傅里葉變換原理，每個(gè)二維細(xì)觀結(jié)構(gòu)分布圖都包含2K×2K個(gè)像素點(diǎn)，這些像素點(diǎn)各自具有力學(xué)性質(zhì)，且相互獨(dú)立，計(jì)算步驟如下：

1）將圖像劃分為2K×2K個(gè)細(xì)胞單元，如圖2所示，取每個(gè)細(xì)胞單元的中心點(diǎn)為像素點(diǎn)，這里稱(chēng)K為像素點(diǎn)個(gè)數(shù)指數(shù)。

2）對(duì)每個(gè)像素點(diǎn)進(jìn)行判斷，如果落在孔隙內(nèi)，彈性模量取為零；如果落在固相內(nèi)，彈性模量取單位值，這樣所獲得的多孔材料的彈性模量為相對(duì)值。

多孔材料各點(diǎn)的彈性張量Cijkl（x）是坐標(biāo)x的函數(shù)，其應(yīng)力－應(yīng)變關(guān)系可表示成

式中，C0ijkl為彈性張量常數(shù)，τij（x）定義為

通過(guò)引入周期性格林張量Γijkl，式（1）的解可表示為［6］

對(duì)式（3）進(jìn)行傅里葉變換有

對(duì)于各向同性材料，（ξ）為［6］

式中，λ和μ為拉梅常數(shù)，εi為傅里葉空間坐標(biāo)。

1.3 迭代求解

方程（1）～方程（4）可以通過(guò)以下迭代方法進(jìn)行求解：

1）給定初始均勻應(yīng)變?chǔ)?ij，由式（1）求得初值應(yīng)力σ0ij；

2）對(duì)于第i＋1次迭代，先由（2）計(jì)算τij（x），對(duì)τij（x）進(jìn)行傅里葉變換求得＾τij（ξ），再檢驗(yàn)收斂性；

3）由式（4）計(jì)算第i＋1次迭代應(yīng)變，再將應(yīng)變進(jìn)行傅里葉反變換；

4）由式（1）計(jì)算應(yīng)力。

一旦迭代收斂，計(jì)算各點(diǎn)的加權(quán)應(yīng)力和加權(quán)應(yīng)變，最后獲得多孔材料的彈性模量。該文以前后兩次迭代值的相對(duì)誤差小于10－3作為收斂準(zhǔn)則。

2 收斂性和有效性驗(yàn)證

2.1 收斂性驗(yàn)證

在下面的計(jì)算中，取正方形邊長(zhǎng)為1 000mm，孔隙率C＝0.2，像素點(diǎn)個(gè)數(shù)為28×28個(gè)，固相材料的泊松比為0.3，計(jì)算所得的相對(duì)彈性模量與孔隙個(gè)數(shù)之間的關(guān)系如圖3所示，其中，虛線表示解析解［8］。從該圖可以看出，當(dāng)孔隙個(gè)數(shù)較少時(shí)，彈性模量上下波動(dòng)，當(dāng)孔隙數(shù)大于100時(shí)，彈性模量基本趨于穩(wěn)定，這與文獻(xiàn)［9］的結(jié)論一致，這是因?yàn)榭紫对蕉?，材料越均勻，像素點(diǎn)所代表的結(jié)構(gòu)單元性質(zhì)越接近于真實(shí)情況。再取孔隙個(gè)數(shù)為100，相對(duì)彈性模量與像素點(diǎn)個(gè)數(shù)指數(shù)之間的關(guān)系如圖4所示。從該圖可以看出，當(dāng)K較小時(shí)，彈性模量上下波動(dòng)幅度較大，隨著K的增大，結(jié)果趨于穩(wěn)定，K＝10與K＝11之間彈性模量的相對(duì)誤差僅為0.3%，表明該方法已經(jīng)收斂。圖4還表示，隨著像素點(diǎn)的增加，數(shù)值解與解析解越接近，這是因?yàn)殡S著像素點(diǎn)的增加，結(jié)構(gòu)單元性質(zhì)得到更細(xì)致的描述，更接近于真實(shí)情況，另一方面，由于傅里葉變換本身具有一定的誤差，當(dāng)像素點(diǎn)超過(guò)一定值后，累積誤差也會(huì)影響計(jì)算結(jié)果。由圖4還可以得出，當(dāng)像素點(diǎn)個(gè)數(shù)為256×256時(shí)，數(shù)值解與解析解最接近。

2.2 有效性性驗(yàn)證

基于前面的討論，取孔隙個(gè)數(shù)為100，像素點(diǎn)個(gè)數(shù)為256×256，相對(duì)彈性模量隨孔隙率變化如圖5所示。從圖5可以看出，數(shù)值解與解析解良好吻合，當(dāng)孔隙率為0.10、0.20、0.30、0.40和0.50時(shí)，兩者之間的誤差分別為6.70%、1.17%、2.12%、3.36%和7.03%，其平均值為4.37%。因此，文中方法的有效性得到初步證實(shí)。

3 結(jié) 論

a.基于Moulinec和Suquet所提出的快速傅里葉變換法，討論了多孔材料彈性模量計(jì)算，通過(guò)與文獻(xiàn)中的解析解比較，初步證實(shí)了該數(shù)值方法的有效性。

b.定量評(píng)價(jià)了孔隙個(gè)數(shù)和像素點(diǎn)個(gè)數(shù)對(duì)計(jì)算結(jié)果的影響，發(fā)現(xiàn)多孔材料越均勻、像素點(diǎn)個(gè)數(shù)越多，數(shù)值解越精確。

［1］ 劉培生.多孔材料引論［M］.北京：清華大學(xué)出版社，2004.

［2］ Torquato S.Random Heterogenerous Materials：Microstructure and Macroscopic Properties［M］.New York：Springerverlag，2001.

［3］ Gibson L J，Ashby M F.Cellular Solids：Structure and Properties［M］.Cambridge：Cambridge University Press，1997.

［4］ Roberts A P，Garboczi E J.Elastic moduli of model random three－dimensional closed－cell cellular solids［J］.Acta materialia，2001，49（2）：189－197.

［5］ Bardenhagena S G，Brydona A D，Guilkey J E.Insight into the physics of foam densification via numerical simulation［J］.Journal of the Mechanics and Physics of Solids，2005，53（3）：597－617.

［6］ Moulinec H，Suquet P.A Numerical Method for Computing the Overall Response of Nonlinear Composites with Complex Microstructure［J］.Computer Methods in Applied Mechanics and Engineering，1998，157（1－2）：69－94.

［7］ 鞠 楊，楊永明，宋振鐸，等.巖石孔隙結(jié)構(gòu)的統(tǒng)計(jì)模型［J］.中國(guó)科學(xué)E輯，2008，38（7）：1026－1041.

［8］ Zheng Q S，Hwang K C.Two－dimensional Elastic Compliances of Materials with Holes and Microcracks［J］.Proceedings of the Royal Society of London，1997，453（1957）：353－364.

［9］ Hu N，Wang B，Tan G W，et al.Effective Elastic Properties of 2－D Solids with Circular Holes：Numerical Simulations［J］.Composites Science and Technology，2000，60（9）：1811－1823.