基于離散元法預(yù)測(cè)碳化硅陶瓷的高溫力學(xué)性能

姜?jiǎng)購(gòu)?qiáng)，譚援強(qiáng)，張高峰，彭銳濤

（湘潭大學(xué)機(jī)械工程學(xué)院，湘潭411105）

0 引 言

陶瓷材料由于具有優(yōu)良的力學(xué)性能而被廣泛應(yīng)用于航天、航空、發(fā)動(dòng)機(jī)、機(jī)械電子等領(lǐng)域。以碳化硅（SiC）陶瓷為例，因其具有優(yōu)良的高溫力學(xué)性能以及較低的熱膨脹系數(shù)、較高的熱導(dǎo)率以及密度小等優(yōu)點(diǎn)，而被廣泛用于制造火箭發(fā)動(dòng)機(jī)燃燒室內(nèi)壁、噴嘴及鼻錐等重要零件［1-2］。從文獻(xiàn)報(bào)道可知［3-5］，近年來(lái)諸多研究人員針對(duì)陶瓷等結(jié)構(gòu)材料的高溫力學(xué)性能進(jìn)行了大量試驗(yàn)，雖然獲得了比較可靠的力學(xué)性能參數(shù)，但試驗(yàn)周期長(zhǎng)，需要消耗大量的人力、物力和財(cái)力，且對(duì)于具有復(fù)雜形狀的陶瓷結(jié)構(gòu)件來(lái)說(shuō)，受試驗(yàn)條件限制，難以獲得其力學(xué)性能參數(shù)及高溫破壞機(jī)理。

隨著計(jì)算機(jī)技術(shù)及數(shù)值分析技術(shù)等的飛速發(fā)展，采用數(shù)值模擬的手段來(lái)研究陶瓷材料在不同溫度下力學(xué)性能的方法得到了廣泛應(yīng)用。朱偉等［6］采用分子動(dòng)力學(xué)法模擬了不同溫度下RDX（三亞甲基三硝胺）晶體的力學(xué)性能，王超等［7］采用分子動(dòng)力學(xué)法模擬了防熱陶瓷材料的高溫力學(xué)性能及斷裂模式。

1971年，Cundall［8］首次提出了適用于巖石力學(xué)的離散元法（DEM），目前該方法已在巖石［9］、混凝土［10］和陶瓷［11-13］等材料中得到了廣泛應(yīng)用。離散元法起源于分子動(dòng)力學(xué)，將材料視為離散顆粒的集合體，顆粒之間用鍵連接，用鍵的斷裂來(lái)描述微裂紋的產(chǎn)生，而且模擬過(guò)程中溫度具有可控性，從而使離散元法成為模擬陶瓷等脆性材料高溫力學(xué)性能及高溫?fù)p傷破壞過(guò)程的一種新手段。

作者基于離散元法建立并校準(zhǔn)了碳化硅陶瓷常溫下的力學(xué)性能，并建立了對(duì)塊體材料離散化的熱力學(xué)模型，將其引入單軸壓縮及三點(diǎn)彎曲等試驗(yàn)來(lái)模擬碳化硅陶瓷在不同溫度下的力學(xué)響應(yīng)，并動(dòng)態(tài)地觀察其損傷破壞過(guò)程。

1 離散元模型的建立及校準(zhǔn)

離散元中顆粒與顆粒之間的相互作用力由位移重疊量及相應(yīng)的力-位移準(zhǔn)則來(lái)確定，顆粒在任意時(shí)步內(nèi)的運(yùn)動(dòng)滿足牛頓第二定律。離散元法將裂紋的生成與擴(kuò)展看成是顆粒與顆粒之間連接鍵的斷裂，能直觀地處理裂紋的形成與擴(kuò)展過(guò)程，如圖1所示。

圖1 離散元模型中的裂紋Fig.1 Cracks in DEM model

為了建立與碳化硅陶瓷微觀結(jié)構(gòu)及力學(xué)行為相似的離散元模型，采用PFC2D［14］軟件生成一個(gè)顆粒隨機(jī)分布且緊密排列的顆粒粘結(jié)（BPM）模型，顆粒之間采用平行鍵連接。眾所周知，陶瓷材料屬于典型的非勻質(zhì)體，晶粒尺寸大小不一，且含有隨機(jī)分布的氣孔和雜質(zhì)。為了使離散元模型和實(shí)際材料在微觀角度上相似，采用不同尺寸的顆粒來(lái)描述晶粒（最大半徑與最小半徑之比為1.5），并使整個(gè)模型的孔隙率與材料的孔隙率相符。除此之外，為了使所建立的離散元模型與實(shí)際材料的力學(xué)性能（如彈性模量、抗壓強(qiáng)度、泊松比等）相匹配，需要對(duì)離散元模型中的參數(shù)（如顆粒剛度、顆粒之間的摩擦、顆粒之間的鍵強(qiáng)度等）進(jìn)行反復(fù)調(diào)試和設(shè)置。目前，這些微觀參數(shù)很難通過(guò)試驗(yàn)的方法來(lái)測(cè)試，將陶瓷材料的宏觀力學(xué)性能與微觀參數(shù)相匹配的方法也還在繼續(xù)探索中。因此，采用數(shù)字試驗(yàn)校準(zhǔn)的方法對(duì)離散元模型進(jìn)行不同的試驗(yàn)?zāi)M，以校準(zhǔn)模型中的微觀參數(shù)。采用單軸壓縮、三點(diǎn)彎曲、單邊切口梁等試驗(yàn)來(lái)校準(zhǔn)所建立的碳化硅離散元模型［12］，結(jié)果如表1所示。前述力學(xué)性能試驗(yàn)過(guò)程中的碳化硅陶瓷通過(guò)真空反應(yīng)燒結(jié)工藝制備（碳化硅微粉粒徑8～10μm，反應(yīng)燒結(jié)溫度1 700℃），燒結(jié)成型后的材料中含約10%（質(zhì)量分?jǐn)?shù)）的游離硅，其余為碳化硅。

表1 碳化硅陶瓷離散元模擬與實(shí)測(cè)主要力學(xué)性能Tab.1 Main mechanical properties of SiC ceramic in experiment and in DEM simulation

2 離散體中的熱力學(xué)模型

2.1 熱傳導(dǎo)方程

在PFC2D 中，采用含有熱存儲(chǔ)器（與顆粒相關(guān)）和熱管（與鍵相關(guān)）的微觀結(jié)構(gòu)模型來(lái)模擬溫度場(chǎng)，熱量在通過(guò)熱存儲(chǔ)器之間的激活熱管中傳遞，如圖2所示。熱管與顆粒之間的接觸同時(shí)存在，但只有在兩個(gè)顆粒之間有重疊或者有連接鍵存在時(shí)，熱管才是激活的。如果構(gòu)成材料的微觀模型在加載過(guò)程中有連接鍵產(chǎn)生或者斷裂，都將改變激活熱管的數(shù)量，從而改變材料的導(dǎo)熱能力。

圖2 PFC2D中顆粒、熱存儲(chǔ)器和熱管的示意Fig.2 Abridged general view of particles，thermal reservoirs and heat pipes in PFC2D

對(duì)于連續(xù)體而言，假定物體應(yīng)變的改變不會(huì)產(chǎn)生溫度的變化，則其熱傳導(dǎo)方程為

式中：qi為沿坐標(biāo)軸xi方向的熱流通量；qv為單位體積內(nèi)的熱功率密度；ρ為質(zhì)量密度；Cv為材料的比熱容；T 為溫度；t為時(shí)間，xi為空間坐標(biāo)；下標(biāo)i為方向指標(biāo)（i＝1，2，3）。

根據(jù)傅里葉變換，可以得到連續(xù)體的熱流通量與溫度梯度之間的關(guān)系：

式中：kij為材料在空間方向上的熱傳導(dǎo)率張量，對(duì)各向同性材料，熱傳導(dǎo)率在各個(gè)方向上保持常數(shù)；xj為空間坐標(biāo)，下標(biāo)j為方向指標(biāo)（j＝1，2，3）。

對(duì)于離散的體系而言，需要將式（1）進(jìn)行轉(zhuǎn)換。假定離散體的總體積為V，單位體積內(nèi)的熱流通量以qi的散度來(lái)表示，在體積為V 的熱存儲(chǔ)器上qi的散度平均值為

對(duì)式（3）通過(guò)高斯差分將其轉(zhuǎn)換為對(duì)面積的積分：

式中：ni為表面向外的單位法向矢量；S 為表面積。

假定離散體中的熱量只在N個(gè)熱管中發(fā)生流動(dòng)，從而式（4）可以轉(zhuǎn)換為

式中：上標(biāo)（p）為與熱管p相關(guān)的變量；Q（p）為從熱存儲(chǔ)器中流經(jīng)熱管p的功率；ΔS（p）為熱管p的橫截面積。

將式（6）代入到式（1）中，可以得到離散體的熱傳導(dǎo)方程：

式中：Qv＝V·qv，為離散體的熱源強(qiáng)度；m 為離散體的質(zhì)量；Cv為連續(xù)體材料的比熱容。

如果把單個(gè)熱管看作單位長(zhǎng)度內(nèi)具有熱阻η的一維單元，則每個(gè)熱管的功率為

式中：ΔT 為熱管兩端的溫度差；l為熱管的長(zhǎng)度。

2.2 熱管熱阻與材料熱傳導(dǎo)率的關(guān)系

在離散元中，假定單個(gè)熱管的熱阻為η，體系中包含N個(gè)熱管，總體積為V，則在材料體積V 內(nèi)的平均熱流通量為

同時(shí)，熱量只在熱管內(nèi)流通，從而上式可以轉(zhuǎn)換為：

式中：A（p）為 熱 管p 的 有 效 橫 截 面 積；l（p）為 熱 管p的長(zhǎng)度。

熱管的熱流通量定義為

將式（8）代入式（11）中，可以得到：

式中：nj為熱管的單位法向矢量；ΔT 為熱管兩端的溫度差。

將式（12）、（13）代入到式（10）中，可以得到：

比較式（14）與式（2）可以得到在離散元體系中熱傳導(dǎo)率張量kij：

從式（15）可知，只要給定單個(gè)熱管的長(zhǎng)度、單位法向矢量、熱阻以及總熱管數(shù)量等參數(shù)就可以計(jì)算出離散單元所描述的材料的熱傳導(dǎo)率。

2.3 熱力學(xué)耦合計(jì)算

在離散元中，材料的熱應(yīng)變通過(guò)顆粒以及顆粒之間連接鍵的熱膨脹來(lái)描述，且溫度引起的應(yīng)變通過(guò)改變顆粒的半徑來(lái)實(shí)現(xiàn)。給定一個(gè)溫度的變化ΔT，可以計(jì)算出粒子半徑R 的改變量ΔR，即：

式中：α為與每個(gè)顆粒相關(guān)的線性熱膨脹系數(shù)。

如果顆粒之間存在連接鍵（如平行鍵），且熱管是激活的，那么連接鍵的長(zhǎng)度ˉL 將在溫度變化時(shí)發(fā)生相應(yīng)的改變，從而改變鍵的受力情況。假定在溫度變化時(shí)，ˉL 的改變只影響法向力的改變，即：

通過(guò)以上的推導(dǎo)可知，在離散元體系中材料的熱傳導(dǎo)率通過(guò)給定單個(gè)顆粒的熱阻來(lái)確定，顆粒與顆粒之間的熱傳導(dǎo)以連接鍵的存在為基礎(chǔ)，若產(chǎn)生新的連接鍵或連接鍵發(fā)生斷裂，則將改變材料的熱傳導(dǎo)能力。熱量所產(chǎn)生的應(yīng)變通過(guò)改變顆粒的半徑或連接鍵的長(zhǎng)度來(lái)實(shí)現(xiàn)，進(jìn)而引起材料內(nèi)部接觸力的改變。

3 碳化硅陶瓷高溫力學(xué)性能的模擬與討論

采用前面已校準(zhǔn)的碳化硅陶瓷微觀參數(shù)，通過(guò)施加不同的溫度場(chǎng)以模擬不同溫度下材料力學(xué)性能的變化，通過(guò)模擬單軸壓縮試驗(yàn)和三點(diǎn)彎曲試驗(yàn)來(lái)獲得材料的彈性模量、泊松比、抗壓強(qiáng)度及抗彎強(qiáng)度等力學(xué)性能參數(shù)，模型如圖3所示。單軸壓縮試驗(yàn)中包含的顆粒數(shù)為7 486個(gè)，加載速度為0.2 m·s-1；三點(diǎn)彎曲試驗(yàn)中包含的顆粒數(shù)為10 675個(gè)，加載速度為0.1m·s-1，顆粒的平均半徑為2μm，最大半徑與最小半徑之比為1.5。為了消除模擬的隨機(jī)性所產(chǎn)生的誤差，每個(gè)模型分別選取5個(gè)不同的隨機(jī)數(shù)進(jìn)行模擬，最后取其平均值作為模擬結(jié)果。模擬過(guò)程中選取9個(gè)不同的溫度場(chǎng)（從100 ℃遞增至1 700 ℃，增量為200 ℃）。

圖3 離散元力學(xué)性能測(cè)試模型Fig.3 DEM model for uniaxial compression test（a）and three-point bending test（b）

3.1 溫度對(duì)力學(xué)性能的影響

從圖4中可以看出，當(dāng)溫度為100～500 ℃時(shí)，材料的應(yīng)力-應(yīng)變曲線基本一致，這說(shuō)明在該溫度區(qū)間時(shí)，溫度對(duì)材料的力學(xué)性能基本沒(méi)有影響。隨著溫度的升高，應(yīng)力的峰值逐漸降低，對(duì)應(yīng)的應(yīng)變也相應(yīng)降低，應(yīng)力-應(yīng)變曲線的斜率也逐漸減小，從而表明材料的彈性模量有所降低。對(duì)比高溫和低溫區(qū)間的失效響應(yīng)曲線，可以看到材料由最開(kāi)始的脆性斷裂（500 ℃以下）逐漸向準(zhǔn)脆性斷裂（500～1 100 ℃）甚至塑性失效轉(zhuǎn)變（如1 300 ℃）。應(yīng)力-應(yīng)變曲線的結(jié)果表明，溫度對(duì)材料力學(xué)性能有明顯的影響，隨著溫度的升高，材料開(kāi)始出現(xiàn)一定程度的軟化，其破壞機(jī)理也相應(yīng)發(fā)生改變。

圖4 不同溫度下碳化硅陶瓷單軸壓縮試驗(yàn)過(guò)程中的應(yīng)力-應(yīng)變曲線Fig.4 Stress-strain curves of SiC ceramic during uniaxial compression at different temperatures

由圖5（a）可知，隨著溫度的升高，材料的彈性模量逐漸降低，泊松比逐漸增大，這表明材料的抗變形能力降低，在高溫作用下發(fā)生了一定的軟化現(xiàn)象。由圖5（b）可見(jiàn)，溫度越高，材料的抗壓強(qiáng)度和抗彎強(qiáng)度越低。力學(xué)性能隨溫度的變化曲線同樣表明，在溫度低于500℃時(shí)，材料的力學(xué)性能變化較小，而溫度一旦持續(xù)升高后，其力學(xué)性能將發(fā)生明顯的改變。

圖5 碳化硅陶瓷力學(xué)性能隨溫度變化的模擬結(jié)果Fig.5 Simulated results of mechanical properties as a function of temperature for SiC ceramic：（a）elastic modulus and Poisson′s ratio and（b）compressive strength and bending strength

3.2 不同溫度下材料的破壞形態(tài)

由圖6可見(jiàn)，當(dāng)顆粒之間的連接鍵所承受的外載荷大于鍵本身的強(qiáng)度時(shí)，連接鍵發(fā)生斷裂，此時(shí)便認(rèn)為是新微裂紋產(chǎn)生了。若一系列的微裂紋連接形成曲折的裂紋路徑時(shí)，宏觀裂紋便產(chǎn)生了。圖6中的白色線段表示裂紋，可見(jiàn)當(dāng)溫度較低（低于500 ℃）時(shí)，裂紋形成的材料損傷線與底面成45°～50°的夾角，從而表明此時(shí)碳化硅陶瓷的破壞形態(tài)以脆性斷裂為主；隨著溫度升高（如700，900℃），裂紋的擴(kuò)展趨勢(shì)出現(xiàn)了一定變化，部分裂紋匯合形成的局部宏觀裂紋不再與底面成45°～50°的夾角，但材料的破壞形態(tài)依然以脆性斷裂為主；當(dāng)溫度升至1 100 ℃以上后，裂紋的擴(kuò)展不再有明顯的規(guī)律，且在試樣中間有“墩粗”的跡象，說(shuō)明此時(shí)材料的破壞機(jī)理不再是以脆性斷裂為主，而是呈現(xiàn)了一定的塑性變形，且在高溫作用下材料出現(xiàn)“墩粗”現(xiàn)象。

圖6 不同溫度下碳化硅陶瓷單軸壓縮破壞后的形貌Fig.6 Failure morphology of SiC ceramic after uniaxial compression test at different temperatures

4 結(jié) 論

（1）建立了脆性材料的離散元熱力學(xué)模型，給出了材料微觀熱力學(xué)參數(shù)與宏觀熱力學(xué)參數(shù)之間的關(guān)系。

（2）建立并校準(zhǔn)碳化硅陶瓷材料的離散元模型，模擬了在不同溫度下材料的單軸壓縮試驗(yàn)及三點(diǎn)彎曲試驗(yàn)；當(dāng)溫度小于500℃時(shí)，溫度場(chǎng)對(duì)材料的力學(xué)響應(yīng)及力學(xué)性能參數(shù)基本沒(méi)有影響，但隨著溫度的逐漸升高，材料的彈性模量逐漸降低，泊松比逐漸增大，抗壓強(qiáng)度和抗彎強(qiáng)度也隨之降低。

（3）結(jié)合材料的應(yīng)力-應(yīng)變曲線、力學(xué)性能參數(shù)以及不同溫度下的單軸壓縮破壞形貌可知，隨著溫度的升高，材料的破壞形態(tài)由脆性斷裂逐漸向塑性變形轉(zhuǎn)變，并會(huì)使材料出現(xiàn)一定的軟化現(xiàn)象。

［1］OHRING M.Materials science of thin films：deposition and structure［M］.Singapore：Academic Press，2006.

［2］DING Y，NORTHWOOD D O，ALPAS A T.Fabrication by magnetrnn sputtering of Al-Cu-Fe quasicrystalline films for tribological applications［J］.Surface and Coating Technology，1997，96：140-147.

［3］劉孝福，婁延春，蘇貴橋，等.定向凝固用陶瓷型殼高溫力學(xué)性能研究現(xiàn)狀［J］.特種鑄造及有色合金，2010，30（10）：913-917，983.

［4］黃小蕭，溫廣武，程顯明.A14SiC4陶瓷的高溫抗氧化性能和高溫力學(xué)性能的研究［J］.材料工程，2004（12）：32-35.

［5］白世鴻，喬生儒，舒武炳，等.先進(jìn)結(jié)構(gòu)陶瓷高溫力學(xué)性能測(cè)試與表征［J］.材料工程，2000（10）：45-48.

［6］朱偉，肖繼軍，馬秀芳，等.不同溫度下RDX 晶體力學(xué)性能的MD模擬［J］.火炸藥學(xué)報(bào)，2007，30（4）：17-19，24.

［7］王超，梁軍，欒旭.防熱陶瓷材料高溫力學(xué)性能及斷裂模式的分子動(dòng)力學(xué)模擬［J］.材料工程，2008（增1）：239-242.

［8］CUNDALL P A，STRACK O D L.A discrete numerical model for granular assembles［J］.Géotechnique，1979，299（1）：47-65.

［9］HUANG H.Discrete element modeling of tool-rock interaction［D］.Minnesota：University of Minnesota，1999.

［10］HZENTZ S，DONZE F V，DAUDEVILLE L.Discrete element modeling of concrete submitted to dynamic loading at high strain rates［J］.Computers and Structures，2004，82：2509-2524.

［11］姜?jiǎng)購(gòu)?qiáng)，譚援強(qiáng)，聶時(shí)君，等.碳化硅陶瓷預(yù)應(yīng)力加工的離散元模擬 與 實(shí) 驗(yàn) 研 究［J］.無(wú) 機(jī) 材 料 學(xué) 報(bào)，2010，25（12）：1286-1290.

［12］TAN Yuan-qiang，YANG Dong-min，SHENG Y .Discrete element method（DEM）modeling of fracture and damage in the machining process of polycrystalline SiC［J］.Journal of the European Ceramic Society，2009，29：1029-1037.

［13］TAN Yuan-qiang，YANG Dong-min，SHENG Y .Study of polycrystalline Al2O3machining cracks using discrete element method［J］.International Journal of Machine Tool ＆ Manufacture，2008，48：975-982.

［14］CUNDALL P A.PFC2D user′s manual［M］.Version 3.Lakewood，Colorado，USA：Itasca Consulting Group Inc，2002.