顆粒增強金屬基復合材料參數(shù)化建模研究

信振洋，王悅，苗文成，陳華輝

中國礦業(yè)大學(北京) 機電與信息工程學院，北京 100083

顆粒增強金屬基復合材料一般由顆粒相、基體相和界面相3部分組成，力學性能良好，在航空航天、車輛、礦用設備等領域有著廣泛的應用。隨著計算機技術的發(fā)展，以有限元為代表的數(shù)值模擬方法逐步成為研究復合材料力學性能的有力工具。顆粒增強復合材料具有較為復雜的細觀結構[1]，主要體現(xiàn)在顆粒的隨機分布和顆粒形狀不規(guī)則兩個方面，導致其幾何模型的構建成為難點。目前顆粒增強復合材料細觀模型主要有2大類：單胞模型(Unit Cell，UC)和代表性體積單元模型(Representative Volume Element，RVE)。單胞模型是將增強體視為規(guī)則排布，從而將整個復合材料簡化為只含有一個或者若干個顆粒的最小單元[2-3]；代表性體積單元模型的增強顆粒數(shù)量較多，從幾個到幾百個不等，隨機分布更接近復合材料的實際結構[4]。

增強顆粒常見有球狀和不規(guī)則狀等，相關研究表明，顆粒形貌越接近實際，模擬結果與實驗值的誤差則越小[5-6]。球狀顆粒形貌較為簡單，可直接采用球狀模型建模；不規(guī)則顆粒的建模相對較為復雜，為使顆粒模型接近實際不規(guī)則顆粒，采用基于真實細觀形貌的三維重構方法構建顆粒模型，通過CT掃描或逐層研磨的方式獲得分層照片，并進一步合成為三維幾何模型[8-10]。這類方法獲得的RVE模型最接近實際的細觀結構形貌，但實施限制較多，成本高，工作量大。為構建接近實際的顆粒模型并且避免三維重構方法的不足，采用程序語言構建隨機凸多面體的方法生成接近實際顆粒的RVE模型[12-14]。建模屬于顆粒增強復合材料力學性能有限元仿真的前處理過程，但是目前建模過程與有限元仿真軟件的集成度并不高，一般需要先用建模軟件生成模型，之后導入有限元軟件計算。這無疑降低了整個仿真流程的效率，容易出現(xiàn)幾何模型與有限元軟件的兼容性問題。因此，有必要將建模算法通過程序語言集成到有限元仿真軟件中，構建更加高效用戶界面，實現(xiàn)復合材料建模與有限元計算的無縫連接。

多顆粒隨機分布代表性體積單元模型的構建可采用隨機順序吸附算法和擾動算法等。隨機順序吸附算法是通過向RVE所占據(jù)的幾何空間中逐一投放顆粒，并判斷當前顆粒是否與已存在顆粒之間發(fā)生干涉，如不干涉則保存當前顆粒的相關數(shù)據(jù)，投放終止條件是達到體積分數(shù)要求或循環(huán)次數(shù)達到最大設定值[4，15]；擾動算法是首先賦予所有顆粒規(guī)則的分布位置，然后對每個顆粒施加隨機距離的擾動量，當所有顆粒滿足無干涉隨機分布之后保存位置信息[16]?；陔S機順序吸附算法采用Python語言，在ABAQUS軟件中實現(xiàn)了顆粒增強復合材料幾何模型的參數(shù)化建模，并封裝為插件(CUMTB RVE Modeling Tool)，該套插件擴展了有限元軟件ABAQUS的功能，程序可移植性強，能有效地降低復合材料仿真的建模難度，提高工作效率。

復合材料力學性能有限元仿真的基本作用之一是對其力學性能進行預測，通過模擬的結果指導復合材料的設計。通過粉末燒結方法制備的ZTAp/Fe45復合材料，具有高硬度和高耐磨性，是一類新型的耐磨復合材料，有望應用于煤礦刮板輸送機和破碎機等設備中，可有效提高相關零部件的耐磨性，降低設備損耗。筆者采用有限元仿真方法分析了ZTA增強顆粒的形狀、體積分數(shù)、界面性能等因素對ZTAp/Fe45復合材料彈性模量的影響，并與實驗結果進行對比。仿真結果表明，不規(guī)則狀顆粒對ZTAp/Fe45復合材料彈性模量的增強效果高于球狀顆粒，顆粒體積分數(shù)的增加和界面彈性模量的提升均能提高復合材料的彈性模量。

1 ABAQUS-Python二次開發(fā)原理

作為商用非線性有限元分析軟件的代表，ABAQUS提供了豐富的程序語言接口，如Python、C/C++、FORTRAN等，以方便用戶通過二次開發(fā)實現(xiàn)特有的功能。前處理過程多采用Python語言，而FORTRAN和C/C++可用來開發(fā)用戶子程序，如自定義材料子程序UMAT/VUMAT以及用戶單元子程序UEL/VUEL等，這些功能擴展增加了ABAQUS軟件的靈活性。

ABAQUS軟件與Python語言的用戶接口如圖1所示，可使用用戶圖形界面(Graphical User Interface，GUI)、命令行或者腳本文件3種方式向ABAQUS自帶的Python解釋器輸入命令，進一步生成input文件提交內核運算。對于簡單的代碼或者不需要大量反復使用的代碼，可采用命令行或腳本文件輸入，但對于需要反復使用或者變量較多的程序，這2種方式操作復雜易出錯。而采用GUI界面輸入的方式則可以降低操作錯誤發(fā)生的概率，用戶不需要關注內部的運行代碼，只需要關注操作界面中參數(shù)的輸入或者選取。本文采用自定義GUI插件的方式實現(xiàn)了ABAQUS軟件中的顆粒增強復合材料的參數(shù)化建模。

圖1 Python語言與ABAQUS/CAE接口Fig.1 Interface of Python and ABAQUS/CAE

2 參數(shù)化建模的實現(xiàn)

2.1 單胞模型

顆粒增強復合材料細觀結構的復雜性決定了采用參數(shù)化建模的必要性。參數(shù)化建模是指通過輸入復合材料的特征參數(shù)，如顆粒形狀、尺寸、體積分數(shù)、分布方式以及模型尺寸等數(shù)據(jù)，構建顆粒增強復合材料幾何模型。這一建模過程需要借助程序語言來實現(xiàn)。

為了實現(xiàn)參數(shù)化建模，針對不同類型的計算模型編寫相應的Python語言程序，并將構建不同幾何模型的程序封裝成相應的插件(CUMTB RVE Modeling Tool)，該插件由3個模塊組成，可用來構建軸對稱單胞(Part-1)、三維周期單胞(Part-2)和三維或二維的隨機分布代表性體積單元(Part-3)的模型。在該插件中，通過直接在用戶界面輸入相關參數(shù)，如模型的名稱、體積分數(shù)、顆粒半徑以及點選顆粒形狀參數(shù)，可實現(xiàn)在ABAQUS軟件中生成對應的復合材料幾何模型，從而實現(xiàn)參數(shù)化建模。

軸對稱單胞[17-20]和三維周期單胞[21-23]是基于顆粒規(guī)則分布的假設，并將復合材料簡化為最小結構單元的建模方式，這2類模型的特點是計算規(guī)模小、效率高，但由于對復合材料的結構簡化較多，不能體現(xiàn)隨機分布的特征，有一定局限性。軸對稱單胞建模插件界面如圖2所示。此類單胞中顆粒模型常見的有切邊圓柱狀、圓柱狀和球狀3種。在軸對稱單胞中，球狀顆粒模型用來模擬實際顆粒形狀為球狀的顆粒；多面體狀顆粒由于不能用回轉體來表征多面體結構，而采用柱狀和切邊圓柱表征顆粒的棱邊和尖角。

圖2 軸對稱單胞建模插件界面Fig.2 Plug-in GUI of the axisymmetric unit cell model

根據(jù)不同的單胞選取規(guī)則[23]，三維單胞可簡化為簡單立方(SC)、體心立方(BCC)、面心立方(FCC)、密排六方(HCP)等，為了能夠施加相對簡單的周期性邊界條件，一般選用立方體形狀的單胞[3，22，25]。如圖3所示，僅對常用的簡單立方、體心立方、面心立方3種單胞模型進行編程和插件封裝。

圖3 三維立方單胞建模插件界面Fig.3 Plug-in GUI of the 3D unit cell model

2.2 隨機分布代表性體積單元模型

不規(guī)則多面體顆粒的形狀具有較強的隨機性，為了降低建模難度，在具體的建模算法中將此類顆粒視為隨機凸多面體，其截面形狀為凸多邊形。對于二維模型，采用顆粒平均尺寸所確定的圓作為控制條件，隨機生成位于圓上的若干點，按逆時針方向(或順時針方向)順次連接各個頂點，即可生成內接于圓的不規(guī)則凸多邊形，多邊形的邊數(shù)被設定為3 ～ 6個。圖4為生成的某個隨機四邊形示意圖。

圖4 隨機多邊形示意圖Fig.4 Schematic diagram of a random polygon

隨機凸多面體的生成可通過球體或者規(guī)則凸多面體空間拓展[26-27]以及二維平面圖形拉伸[28]等方法實現(xiàn)。本程序中采用了二維平面圖形拉伸方法得到三維顆粒模型，顆粒構建如圖5所示。

首先在XOY平面內構建形狀隨機的三角形或四邊形，然后在Z軸方向拉伸一定厚度得到凸多面體，再對所得凸多面體進行倒角處理得到所需的隨機形貌的凸多面體。圖5中顆粒的長軸和短軸分別為l0和b0，其比值l0/b0、l1/l0和b1/b0為變形系數(shù)，在長短軸確定的情況下變形系數(shù)為0 ～ 0.5的隨機數(shù)。二維截面構建完成后，在垂直紙面的法向拉伸一個等于短軸長度b0的厚度，從而得到三維顆粒模型，短軸長度b0等于顆粒的平均尺寸。根據(jù)對一定數(shù)量的ZTA顆粒形貌尺寸的統(tǒng)計，其顆粒的長短軸比值為1 ～ 1.5，因此建模過程中l(wèi)0/b0為1 ～ 1.5的隨機數(shù)。

顆粒的隨機分布采用隨機順序吸附算法實現(xiàn)[4，29-30]，該算法的基本流程如圖6所示。

顆粒的無干涉隨機分布是構建隨機分布代表性體積元的關鍵，球狀和圓形的顆粒干涉采用球心距離判斷，而凸多面體和凸多邊形可根據(jù)體積或面積判斷。例如，某多面體一個頂點侵入到另一個多面體內部時，可用該點與被侵入的多面體各個面所組成的多面體體積之和與被侵入多面體體積相比較即可判斷。顆粒隨機分布建模插件操作界面如圖7所示，通過該插件可構建圓形或隨機凸多邊形顆粒的二維RVE和球狀或隨機凸多面體顆粒的三維RVE。

圖7 隨機分布代表性體積單元插件界面Fig.7 Plug-in GUI of RVEs

3 ZTAp/Fe45復合材料力學性能分析

圖8所示為實驗制備的純Fe45和ZTAp/Fe45復合材料拉伸試樣。ZTAp/Fe45試樣的ZTA顆粒體積分數(shù)φ(ZTA)為20%，形狀為不規(guī)則形貌，顆粒平均尺寸為2 mm，顆粒形貌如圖9所示。采用的實驗設備為WDW-200微機控制電子式萬能實驗機，拉伸速度為1 mm/min。

圖8 純Fe45與ZTAp/Fe45拉伸試樣Fig.8 Tensile specimens of pure Fe45 and ZTAp/Fe45

圖9 不規(guī)則ZTA增強顆粒形貌Fig.9 Irregular morphological of ZTA particles

通過單軸拉伸實驗測得該批次試樣中基體材料Fe45平均彈性模量為187.7 GPa，泊松比為0.3；20%體積分數(shù)的ZTAp/Fe45彈性模量為210.7 GPa。

3.1 顆粒形狀對ZTAp/Fe45力學性能的影響

ZTA增強顆粒通常有不規(guī)則狀和球狀2種形狀。ZTA增強顆粒體積分數(shù)φ(ZTA)為20%的2種RVE模型如圖10所示，顆粒平均尺寸為2 mm，代表性體積單元在X、Y、Z的3個方向尺寸均為8 mm。為了降低模型隨機性所引起的誤差，2種模型分別構建10個，共20個代表性體積單元模型，計算每種模型彈性模量的平均值。有限元計算模型中ZTA顆粒彈性模量為338.0 GPa，抗拉強度為344 MPa，泊松比為0.23[31]。

圖10 ZTAp/Fe45代表性體積單元Fig.10 RVEs of ZTAp/Fe45

上述計算模型的有限元網(wǎng)格剖分示意圖和邊界條件如圖11所示，網(wǎng)格控制尺寸為0.3 mm，采用C3D4四面體單元，模型的網(wǎng)格數(shù)量在10×104～16×104個，代表性體積單元網(wǎng)格數(shù)量的差異與顆粒形狀和分布有關。顆粒網(wǎng)格與基體網(wǎng)格之間采用共節(jié)點的方式連接，即理想結合。模型邊界條件為左側面固定，右側面與參考點耦合，對參考點施加應變量為0.3%的X軸正方向微小位移荷載，采用ABAQUS/Standard隱式求解器進行計算。

圖11 代表性體積單元網(wǎng)格與邊界條件Fig.11 RVEs meshes and boundary conditions

不規(guī)則多面體顆粒代表性體積單元和球狀顆粒代表性體積單元計算得到的平均彈性模量分別為213.1 GPa和203.8 GPa，不規(guī)則多面體代表性體積單元所得結果與實驗值210.7 GPa更為接近，但仍然高于實驗值。產(chǎn)生這種誤差的主要原因在于模型采用理想結合界面提高了基體與顆粒之間的載荷傳遞能力，導致計算結果相對實驗值偏大。造成2種RVE模型彈性模量結果差異的主要原因在于顆粒形狀對復合材料內部細觀應力應變場影響不同。應變量為0.09%時ZTAp/Fe45的增強體Mises應力云圖如圖12所示。

圖12 ZTAp/Fe45的增強體Mises應力云圖[φ(ZTA)為20%，應變量為0.09%]Fig.12 Mises stress nephogram of ZTA particles in ZTAp/Fe45 [φ(ZTA)=20%，strain=0.09%]

在應變量相同的情況下，不規(guī)則多面體顆粒高應力區(qū)域大于球狀顆粒模型，承受的荷載更大，在顆粒的尖角和棱邊附近由于應力集中現(xiàn)象出現(xiàn)明顯的高應力區(qū)，圖12(a)應力云圖中高于ZTA抗拉強度(344 MPa)的灰色部分出現(xiàn)在該區(qū)域。隨著形變量的持續(xù)增大，這些灰色區(qū)域將成為顆粒損傷的萌生區(qū)域；在相同應變量的情況下，圖12(b)中的球狀顆粒并未出現(xiàn)高于344 MPa的區(qū)域，最高應力值相對于不規(guī)則顆粒降低約50 MPa，說明在承受相同應變量時，不規(guī)則顆粒更易發(fā)生損傷。

圖13為代表性體積單元與X軸平行截面的Mises應力云圖。由圖13可知，對于2種顆粒代表性體積單元模型，復合材料的高應力區(qū)域均出現(xiàn)在增強顆粒中，基體所承受的應力明顯低于顆粒，說明ZTAp/Fe45復合材料在形變過程中增強顆粒表現(xiàn)出較好的承載能力。

圖13 ZTAp/Fe45截面Mises應力云圖(應變量為0.09%)Fig.13 Mises stress nephogram section of ZTAp/Fe45 (strain=0.09%)

3.2 ZTA顆粒體積分數(shù)對復合材料力學性能的影響

為分析ZTA顆粒體積分數(shù)對復合材料力學性能的影響，構建體積分數(shù)分別為5%、10%和15%的ZTAp/Fe45復合材料RVE模型，顆粒形狀分別采用不規(guī)則多面體和球狀，2種形狀的顆粒均構建10個隨機分布的代表性體積單元，總共60個模型，示例模型如圖14所示。

φa(ZTA)—不規(guī)則顆粒體積分數(shù)；φb(ZTA)—球狀顆粒體積分數(shù)圖14 不同體積分數(shù)代表性體積單元模型Fig.14 RVEs with different volume fractions

在不同體積分數(shù)情況下，不規(guī)則顆粒RVE彈性模量計算結果均高于球狀顆粒RVE計算結果，如圖15所示。但是由于模型中顆粒形狀以及分布方式的隨機性，導致相同體積分數(shù)的模型計算結果亦會出現(xiàn)一定偏差。對于相同類型的顆粒來說，隨著體積分數(shù)的增加會引起復合材料彈性模量的提升。

圖15 顆粒體積分數(shù)與復合材料彈性模量曲線Fig.15 Volume fraction and elastic modulus curve of composites

體積分數(shù)為10%的增強體在應變量0.09%時的應力云圖如圖16所示。對比圖12可知，在相同的應變量下，φ(ZTA)為10%的ZTAp/Fe45中顆粒所承受的應力相對低于φ(ZTA)為20%的ZTAp/Fe45，說明在較高體積分數(shù)模型中顆粒所承受的應力更大。這主要是由于高體積分數(shù)的情況下，顆粒間距更近，更易形成應力集中區(qū)域，同時也意味著在形變過程中顆粒更易發(fā)生破壞。

圖 16 ZTAp/Fe45復合材料增強體Mises應力云圖(φ(ZTA)=10%，應變量為0.09%)Fig.16 Mises stress nephogram of ZTA particles in ZTAp/Fe45 (φ(ZTA)=10%，strain=0.09%)

3.3 界面性質對ZTAp/Fe45力學性能的影響

在復合材料形變過程中，界面的力學性能對基體和增強體之間的荷載傳遞起關鍵作用，并通過這種方式對復合材料宏觀力學性能產(chǎn)生影響。界面力學性能的模擬可采用黏聚力模型(Cohesive Zone Model，CZM)[21]和XFEM[32-33]等方法；另外，在小變形的情況下亦可采用線彈性力學模型近似表征界面的彈性力學性能[34]。本文采用線彈性力學模型分析小變形情況下界面彈性模量對復合材料彈性模量的影響。有限元模型中顆粒和基體采用C3D4單元進行網(wǎng)格剖分，界面單元采用C3D6三棱柱實體單元。界面單元的插入采用基于Python語言開發(fā)的內部程序iCoh實現(xiàn)，模型中界面單元的法向厚度設定為10 μm，其界面單元示意如圖17所示。

圖17 界面單元示意圖Fig.17 Interface element diagram

由于界面性能的測定較為困難，因此通過對界面設定一系列彈性模量值研究界面彈性模量對復合材料彈性模量的影響，界面的泊松比假定與基體材料泊松比一致。圖18為體積分數(shù)20%的ZTAp/Fe45復合材料采用不同彈性模量的界面所得的彈性模量曲線。

圖18 界面彈性模量與ZTAp/Fe45彈性模量曲線Fig.18 Interface Young’s modulus and Young’s modulus curve of ZTAp/Fe45

由圖18可知，在界面彈性模量低于基體彈性模量(187.7 GPa)的情況下，界面模量的增加使復合材料彈性模量提升明顯；當界面模量在187.7～ 676.0 GPa(2倍的增強體彈性模量)時，復合材料彈性模量增長趨緩，此時界面模量的增加并不能顯著提升復合材料的彈性模量。

圖19 ZTAp/Fe45復合材料應力云圖截面(應變量為0.09%)Fig.19 Mises stress nephogram section of ZTAp/Fe45 (strain=0.09%)

圖19所示為界面模量為9.4 GPa(基體彈性模量的1/20)和676.0 GPa情況下復合材料的應力云圖，應變量均為0.09%。圖19(a)中，界面區(qū)域為明顯低應力區(qū)，對荷載的傳遞能力較弱，此時通過提升界面的彈性模量可以改善其傳遞荷載的能力，提高復合材料彈性模量；圖19(b)中，界面區(qū)域為復合材料的高應力區(qū)，表現(xiàn)出較強的承載能力，而作為復合材料主體的基體相和顆粒相則表現(xiàn)出相對低的應力。根據(jù)復合材料細觀力學均勻化理論[35]，雖然界面相為高應力區(qū)，但由于界面相在復合材料中占據(jù)的體積分數(shù)權重較小，此時復合材料彈性模量的主導因素轉化為占更大體積分數(shù)權重的基體相和顆粒相。因此，在圖18曲線中出現(xiàn)當界面彈性模量達到基體彈性模量之后復合材料彈性模量提升不明顯的現(xiàn)象。

4 結 論

(1) 通過ABAQUS軟件的Python語言接口實現(xiàn)了顆粒增強復合材料單胞及隨機分布代表性體積單元等復雜細觀模型的參數(shù)化建模，并封裝為具有操作友好性的插件，解決了顆粒增強復合材料建模難的問題，提升了軟件的易用性和工作效率。

(2) 不規(guī)則多面體顆粒對復合材料彈性模量的增強作用大于球狀顆粒，但由于其不規(guī)則的形貌特點造成了復合材料內部更易出現(xiàn)應力集中現(xiàn)象。ZTA顆粒體積分數(shù)的增加使ZTAp/Fe45復合材料彈性模量呈現(xiàn)上升的趨勢，說明增強顆粒的加入可有效地提升復合材料的彈性模量。

(3) 在低界面模量的情況下，界面結合性能的改善可顯著地提升復合材料的彈性模量，在界面模量達到基體模量之后，繼續(xù)提高界面彈性模量對復合材料彈性模量的提升作用逐漸趨緩。