吳宇航， 肖映雄， 徐亞飛

(湘潭大學(xué) 土木工程與力學(xué)學(xué)院，湘潭 411105)

1 引 言

混凝土在細(xì)觀層次上是由粗骨料、砂漿及兩者間過(guò)渡區(qū)(界面層)組成的三相復(fù)合材料，各相材料具有不同的力學(xué)特性，其中骨料的體積分?jǐn)?shù)、形狀、空間分布和取向及界面層對(duì)其宏觀性能影響很大，以往主要基于宏觀層次的混凝土力學(xué)研究已不能充分解釋混凝土材料的損傷和破壞機(jī)理。由于骨料形態(tài)的復(fù)雜性和空間分布的隨機(jī)性，需要采用有限元方法對(duì)這種復(fù)合材料進(jìn)行數(shù)值模擬。

建立能反映混凝土實(shí)際骨料級(jí)配、含量及其形態(tài)的隨機(jī)骨料模型是進(jìn)行混凝土細(xì)觀力學(xué)數(shù)值模擬的前提和基礎(chǔ)。自Roelfstra等[1]最早提出混凝土細(xì)觀模型的概念以來(lái)，廣大學(xué)者對(duì)混凝土細(xì)觀結(jié)構(gòu)三維模型的建立展開(kāi)了進(jìn)一步的研究[2-13,16,17]。劉光延等[2]率先提出隨機(jī)骨料模型，并將其運(yùn)用于混凝土材料的斷裂模擬；覃偉平等[3]進(jìn)行了三維凸型骨料隨機(jī)投放算法研究，并提出投放區(qū)域均勻化思想；關(guān)振群等[4]運(yùn)用橢球形狀骨料顆粒模擬多相復(fù)合材料模型，并提出基于映射法的骨料表面有限元網(wǎng)格生成算法；秦武等[5,6]運(yùn)用CT切片技術(shù)對(duì)三維模型進(jìn)行重構(gòu)及損傷過(guò)程分析；方秦等[7,8]基于MATLAB自編程序生成混凝土細(xì)觀模型，并嵌入COMSOL和 LS -DYNA 等軟件進(jìn)行前處理和數(shù)值模擬；Caballero等[9，10]基于Voronoi 圖形中的胞元形狀模擬混凝土細(xì)觀三維結(jié)構(gòu)；郭瑞奇等[11]利用Fortran程序結(jié)合ANSYS的APDL命令流建立混凝土細(xì)觀模型。除此之外，還有基于Visual C++編寫(xiě)的骨料模型[12]和基于AutoCAD內(nèi)置VBA語(yǔ)言編寫(xiě)的骨料模型[13]等。綜上，混凝土三維細(xì)觀建模的相關(guān)研究很多，但缺點(diǎn)也很明顯，研究人員或困于程序編寫(xiě)的復(fù)雜性，或困于軟件功能的局限性，往往需要熟練運(yùn)用多個(gè)軟件并配合使用，這樣的處理方式帶來(lái)了諸多不便。

本文運(yùn)用Python腳本對(duì)Abaqus軟件進(jìn)行二次開(kāi)發(fā)，實(shí)現(xiàn)了隨機(jī)骨料模型建模、劃分網(wǎng)格、加載和提交計(jì)算等功能的一體化、自動(dòng)化處理，極大地簡(jiǎn)化了混凝土細(xì)觀力學(xué)研究的工作量。二次開(kāi)發(fā)結(jié)合了參數(shù)化建模的思想，提出球形-橢球一體化生成算法，生成了球形和橢球骨料模型。提出頂點(diǎn)控制算法，生成了隨機(jī)凹凸型多面體骨料模型。同時(shí)，提出了布爾切割入侵判別算法，分別對(duì)橢球骨料和多面體骨料進(jìn)行了投放試驗(yàn)。該算法避免了傳統(tǒng)骨料相交判別時(shí)復(fù)雜的數(shù)學(xué)運(yùn)算，且提高了骨料體積含量。本文建模算法可以做到骨料、界面層和基體之間相互獨(dú)立，在進(jìn)行網(wǎng)格剖分時(shí)可避免復(fù)雜的單元屬性判別，得到的網(wǎng)格剖分滿(mǎn)足骨料、砂漿及界面層網(wǎng)格協(xié)調(diào)性要求。最后，以多面體骨料模型為例進(jìn)行有限元數(shù)值模擬，進(jìn)一步驗(yàn)證本文提出的建模方法的有效性。

2 三維球形和橢球形骨料模型的生成

參數(shù)化建模的優(yōu)勢(shì)在于參數(shù)的可控性，首先設(shè)計(jì)了一種球形-橢球一體化生成算法，即通過(guò)控制橢球的參數(shù)-長(zhǎng)短軸之比，來(lái)調(diào)控生成的顆粒種類(lèi)是球體或是橢球，避免了重新編寫(xiě)算法和建模。其次，提出布爾切割入侵判別算法，提高了判別效率及骨料體積含量，并對(duì)橢球型骨料分別進(jìn)行了區(qū)域?yàn)?50 mm×450 mm×450 mm的四級(jí)配、300 mm×300 mm×300 mm的三級(jí)配和150 mm×150 mm×150 mm的二級(jí)配投放試驗(yàn)。參見(jiàn)文獻(xiàn)[7]，各骨料粒徑范圍和占比列入表1。

表1 不同級(jí)配混凝土粒徑范圍Tab.1 Particle size range of concrete with different gradations

2.1 骨料基本參數(shù)的確定

在混凝土試件中，骨料各種參數(shù)都是隨機(jī)的，Python語(yǔ)言的random模塊擁有多種生成隨機(jī)數(shù)的函數(shù)以供使用。調(diào)用隨機(jī)數(shù)，生成外接球的半徑和空間坐標(biāo)分別為

Rn=random.uniform(Rmin,Rmax)

(1)

(2)

(3)

式中a為橢球沿x方向長(zhǎng)半軸,b為橢球沿y方向短半軸；由于Abaqus建模是從二維草圖得到三維部件，所以無(wú)論球體或者橢球都是從半圓或圓弧繞構(gòu)造線旋轉(zhuǎn)得到的，因此橢球沿z方向極半徑c=b；ratio為短軸與長(zhǎng)軸之比，通過(guò)控制此參數(shù)，即當(dāng)ratio=1時(shí)，生成的骨料顆粒為球體；當(dāng)ratio≠1，生成的骨料顆粒為短長(zhǎng)軸之比在一定范圍內(nèi)隨機(jī)的不同形態(tài)橢球。

2.2 骨料的生成

球形和橢球形顆粒的體積可通過(guò)Abaqus腳本內(nèi)置方法getVolume()直接獲得，每生成一個(gè)骨料顆粒，將其體積計(jì)入該粒徑范圍內(nèi)的體積Vi(i=1～4)中，總體積為VT=∑Vi，參考表1的數(shù)據(jù)可知，不同粒徑范圍生成的骨料占比一定，骨料顆粒中半徑(長(zhǎng)半軸、短半軸和極半徑)數(shù)據(jù)按照骨料種類(lèi)從大到小順序生成，詳細(xì)流程如圖1所示。

圖1 骨料生成流程Fig.1 Generating process of random aggregate particles

2.3 骨料的投放與界面層處理

根據(jù)體分比和級(jí)配生成骨料顆粒后，進(jìn)行投放試驗(yàn)，骨料投放的關(guān)鍵在于入侵判別，對(duì)于球形骨料而言，可直接利用式(4)進(jìn)行判別，即

(4)

利用循環(huán)語(yǔ)句(foriin range (n):)對(duì)前n-1個(gè)骨料顆粒進(jìn)行遍歷，當(dāng)?shù)趎個(gè)骨料與前n-1個(gè)骨料都不相交時(shí)，此骨料投放成功。球形骨料投放結(jié)果總結(jié)如圖2所示，結(jié)果表明，在三級(jí)配下可投放球形骨料的體分比超過(guò)55%。

對(duì)于橢球骨料也可以運(yùn)用式(4)，通過(guò)外接球來(lái)進(jìn)行入侵判別，但此方法對(duì)投放空間使用不充分，體積損失較大，導(dǎo)致投放體分比不高。本文提出了一種可提高骨料體積含量的方法，為敘述方便，將其簡(jiǎn)稱(chēng)為布爾切割入侵判別算法。該方法的具體思路如下。

圖2 球形骨料投放結(jié)果三視圖Fig.2 3D views of spherical aggregate after packing process

(1) 按從大到小順序依次選取生成好的骨料。

(2) 隨機(jī)位置和旋轉(zhuǎn)角度，將選取的第n個(gè)骨 料投放到指定區(qū)域中。

(3) 對(duì)第n個(gè)骨料與前n-1個(gè)進(jìn)行布爾切割操作，利用Python語(yǔ)言中的try…except語(yǔ)句測(cè)試切割是否異常，若沒(méi)有異常，則該骨料與之前的骨料存在相交(投放失敗)，若出現(xiàn)異常，則該骨料與之前的骨料不存在相交情況(投放成功)。

(4) 投放失敗的骨料回到第二步繼續(xù)投放，投放成功的骨料與前n-1個(gè)骨料布爾合并。

(5) 將該骨料投放成功時(shí)的坐標(biāo)和角度數(shù)據(jù)存放待用。

另外，為減小高體分比投放時(shí)的計(jì)算量，還可以在進(jìn)行布爾切割入侵判別之前進(jìn)行一次簡(jiǎn)單的初判，判別方法是將式(4)稍作改動(dòng)為

(5)

式中bn為第n個(gè)橢球骨料的短半軸，即先用短半軸進(jìn)行初步判斷，若第n個(gè)橢球骨料的短半軸與前n-1個(gè)骨料的短半軸存在相交情況，則無(wú)需進(jìn)行布爾切割入侵判別，直接重新投放。

界面過(guò)渡區(qū)ITZ(Interfacial Transition Zone)是混凝土材料三相之一，其模量小于其他兩相。Yang[14]研究發(fā)現(xiàn)，不同厚度范圍內(nèi)界面層的彈性模量不同，可將其視為混凝土材料的強(qiáng)度極限相，是建立混凝土細(xì)觀模型的關(guān)鍵。受制于計(jì)算機(jī)計(jì)算能力和時(shí)間代價(jià)，可以取偏大一些的界面層厚度值，也能滿(mǎn)足計(jì)算精度的要求[15]，一般為0.05 mm～0.5 mm。本文取特大石為1.0 mm、大石為0.5 mm、中石為0.3 mm及小石為0.1 mm。

采用布爾運(yùn)算建立獨(dú)立的界面層部件，首先將生成好的球形(橢球)骨料由實(shí)體轉(zhuǎn)化為殼體，再將殼體向外整體偏移，偏移距離為界面層的厚度，然后將新的殼體轉(zhuǎn)化為實(shí)體，再與之前的球形(橢球)骨料做布爾切割運(yùn)算，即可生成其界面層。對(duì)于投放成功的骨料，會(huì)將其位置坐標(biāo)和隨機(jī)旋轉(zhuǎn)角度存入數(shù)組中，該骨料界面層也會(huì)按此坐標(biāo)和旋轉(zhuǎn)角度投放到指定區(qū)域中與之匹配。骨料投放流程如 圖3 所示，可以生成不同級(jí)配對(duì)應(yīng)的橢球形骨料和界面層，圖4給出了三級(jí)配對(duì)應(yīng)的橢球形骨料和界面層及整體分布。由此可見(jiàn)，本文方法可有效生成各種級(jí)配下含橢球骨料及界面層的細(xì)觀隨機(jī)模型，驗(yàn)證了所編Python腳本程序的正確性和適應(yīng)性。

圖3 骨料投放流程Fig.3 Packing process of random aggregate particles

圖4 三級(jí)配橢球形骨料、界面層及整體分布Fig.4 Three -graded ellipsoidal aggregates,ITZs and overall distribution

3 隨機(jī)凹凸型多面體骨料模型的生成

本文提出一種基于凹凸型24面體的多面體(碎石)骨料模型的建模算法，同時(shí)采用布爾切割入侵判別算法進(jìn)行投放試驗(yàn)，得到的骨料模型相較于單純的凸骨料更接近于真實(shí)形態(tài)。

3.1 單個(gè)多面體骨料的生成

利用Abaqus腳本接口提供的豐富的方法與函數(shù)，可以完成單個(gè)多面體骨料從頂點(diǎn)到邊到面到體的生成過(guò)程。生成思路如下。

(6)

(7)

(2) 將頂點(diǎn)連接成邊。通過(guò)調(diào)用WirePolyLine()來(lái)實(shí)現(xiàn)頂點(diǎn)連接，該方法的必需參數(shù)為兩端點(diǎn)的坐標(biāo)，其可通過(guò)索引數(shù)據(jù)庫(kù)中的頂點(diǎn)數(shù)據(jù)獲取。為避免出現(xiàn)多余的線段，頂點(diǎn)之間相連順序遵循以下原則，先連接豎向邊，如P1到P5、P5到P9等；再連接橫向邊，如P2到P3、P3到P4等；最后連接斜向邊，如P2到P7、P3到P8以此類(lèi)推，如圖5(b)所示。

(3) 將各邊覆蓋為面。當(dāng)14個(gè)頂點(diǎn)連接完成后，其表面由各邊分割為24個(gè)區(qū)域。如圖5(b)所示，由于頂點(diǎn)的隨機(jī)波動(dòng)，如P2，P5，P6和P9四個(gè)頂點(diǎn)并不一定在一個(gè)平面上，所以相鄰兩個(gè)面S5和S6之間存在隨機(jī)凹凸性。然后調(diào)用方法CoverEdges()，將各邊封閉為三角平面，如圖5(c)所示。

(4) 將各面填充為體。當(dāng)外表面全部生成完畢后，調(diào)用方法AddCells()，將由外表面構(gòu)成的殼體轉(zhuǎn)化為實(shí)體，即得到單個(gè)隨機(jī)凹凸型多面體骨料。圖6給出了不同角度波動(dòng)范圍下生成的形態(tài)各異的幾個(gè)典型的多面體骨料。

圖5 單個(gè)多面體骨料生成流程Fig.5 Generation process of single polyhedral aggregate

圖6 不同角度波動(dòng)范圍下形態(tài)各異的多面體骨料Fig.6 Polyhedral aggregates with different shapes in different angle ranges

3.2 多面體骨料的投放與界面層處理

根據(jù)圖1流程生成多面體骨料后，進(jìn)行投放試驗(yàn)，投放流程如圖3所示。同樣，在高體分比投放時(shí)，可以在布爾切割入侵判別法之前先進(jìn)行簡(jiǎn)單的初步判別，以減小計(jì)算量，提高投放效率。具體判別方法是將式(4)改動(dòng)為

(8)

式中μ=0.5～0.8，即兩多面體骨料中心點(diǎn)之間的距離小于多面體粒徑之和乘以一個(gè)系數(shù)(系數(shù)大小可根據(jù)需要選取)，即可進(jìn)行下一步布爾切割入侵判別，否則重新選擇位置和旋轉(zhuǎn)角度繼續(xù)投放，如此，可以很大程度上減少布爾切割運(yùn)算次數(shù)，進(jìn)一步控制計(jì)算量。

多面體骨料的界面層生成步驟具體是在已有的頂點(diǎn)數(shù)據(jù)庫(kù)基礎(chǔ)上，將多面體骨料的頂點(diǎn)向外延伸為

(9)

圖7 含界面層的骨料截面圖Fig.7 Section of some aggregate with ITZ

圖8 不同厚度的ITZFig.8ITZ with different thickness

圖9 三級(jí)配多面體骨料、界面層及整體分布Fig.9 Three -graded Polyhedral aggregates,ITZs and overall distribution

4 三維細(xì)觀模型靜力學(xué)數(shù)值模擬

受限于計(jì)算機(jī)性能，本文取體分比為35%的二級(jí)配多面體骨料模型，進(jìn)行單軸壓縮靜力學(xué)數(shù)值模擬。ABAQUS的混凝土塑性損傷模型是基于塑性的連續(xù)介質(zhì)損傷模型，該模型可用于單向加載、多向加載、循環(huán)加載及動(dòng)態(tài)加載等情況,具有較好的收斂性。故而選用塑性損傷模型(CONCRETE_DAMAGE_PLASTICITY)模擬砂漿基體和界面層的力學(xué)性能，用彈性模型模擬粗細(xì)骨料的力學(xué)性能，更符合材料非線性研究要求。各項(xiàng)材料的基本參數(shù)列入表2。

表2 混凝土細(xì)觀三相材料參數(shù)Tab.2 Three -phase material parameters of concrete

采用C3D10四面體單元對(duì)細(xì)觀模型進(jìn)行網(wǎng)格劃分的結(jié)果如圖10所示。由于程序已將骨料、砂漿和界面層自動(dòng)分離，在進(jìn)行網(wǎng)格剖分時(shí)避免了復(fù)雜的單元屬性判別，提高了計(jì)算效率。可以看出，所得網(wǎng)格剖分完全能滿(mǎn)足有限元分析的需要。

為提高細(xì)觀模擬的效率和保證收斂，采用準(zhǔn)靜態(tài)分析(Quasi-static Analysis)方法，即用動(dòng)態(tài)顯式算法(Dynamic-Explicit)來(lái)求解復(fù)雜的混凝土細(xì)觀靜力學(xué)問(wèn)題。單軸壓縮損傷云圖如圖11(a～c)所示，應(yīng)力應(yīng)變曲線及與文獻(xiàn)[10]中圖10(a)的結(jié)果比較如圖11(d)所示。從計(jì)算結(jié)果得到的混凝土損傷形態(tài)和應(yīng)力-應(yīng)變曲線可以看出，基于本文算法生成的細(xì)觀模型得到的應(yīng)力應(yīng)變曲線與混凝土標(biāo)準(zhǔn)曲線基本吻合，能夠較好地將初始線彈性階段、中間應(yīng)力強(qiáng)化階段和末期應(yīng)力軟化階段三個(gè)部分體現(xiàn)出來(lái)。破壞形態(tài)與實(shí)驗(yàn)室標(biāo)準(zhǔn)試件的試驗(yàn)結(jié)果相似，也表現(xiàn)出八字型裂縫，說(shuō)明本文提出的混凝土細(xì)觀模型在模擬混凝土力學(xué)特性方面，具有一定的可靠性。

圖10 多面體骨料網(wǎng)格剖分Fig.10 Mesh generation of polyhedron random aggregates

圖11 混凝土單軸壓縮細(xì)觀模擬結(jié)果Fig.11 Mesoscopic simulation results of concrete uniaxial compression

5 結(jié) 論

在進(jìn)行混凝土細(xì)觀力學(xué)數(shù)值模擬分析時(shí)，首先需要建立能反映混凝土實(shí)際骨料級(jí)配、含量及其形態(tài)的隨機(jī)骨料模型。本文利用Python-Abaqus二次開(kāi)發(fā)平臺(tái)，對(duì)三維混凝土隨機(jī)骨料模型的參數(shù)化建模與數(shù)值模擬進(jìn)行了研究，并得到以下結(jié)論。

(1) 利用提出的球形-橢球一體化生成算法和頂點(diǎn)控制算法,可有效生成各種級(jí)配下含球形、橢球形及任意多面體骨料及界面層的細(xì)觀隨機(jī)模型。

(2) 利用提出的能提高骨料體積含量的布爾切割入侵判別算法成功地對(duì)橢球骨料和多面體骨料進(jìn)行了投放試驗(yàn)，獲得的混凝土細(xì)觀模型其骨料、界面層及基體已自動(dòng)分離開(kāi)來(lái)，便于后續(xù)材料屬性的賦予及網(wǎng)格剖分與加載。

(3) 基于Abaqus顯式求解器對(duì)混凝土試件進(jìn)行單軸壓縮靜力學(xué)數(shù)值模擬，驗(yàn)證了本文建模算法得到的混凝土細(xì)觀模型的可靠性。

利用Python對(duì)Abaqus進(jìn)行二次開(kāi)發(fā)優(yōu)勢(shì)明顯，其可以方便多次重復(fù)或局部修改試算等，為進(jìn)行混凝土細(xì)觀力學(xué)幾何建模、劃分網(wǎng)格、指定材料屬性、提交作業(yè)及后處理分析等提供自動(dòng)化處理程序設(shè)計(jì)，顯著提高了工作效率。另外，本文僅做了簡(jiǎn)單的單軸拉壓仿真試驗(yàn)，還需要對(duì)更多問(wèn)題進(jìn)行更廣泛的數(shù)值試驗(yàn)，以進(jìn)一步測(cè)試Python編程在顯著提高Abaqus仿真效率上的有效性和適應(yīng)性，這些將是今后進(jìn)一步研究的內(nèi)容和方向。