基于預(yù)插內(nèi)聚力單元的混凝土模型生成方法與單軸拉伸細(xì)觀模擬研究

馬 峰，田 羽，張小飛

(1.甘肅省水利水電勘測設(shè)計研究院有限責(zé)任公司，甘肅 蘭州 730000；2.中國電建集團(tuán)北京勘測設(shè)計研究院有限公司，北京 100024；3.廣西大學(xué)土木建筑工程學(xué)院，廣西 南寧 530004)

0 引 言

混凝土由于具有良好的可塑性、耐久性及經(jīng)濟(jì)性等優(yōu)點，被廣泛應(yīng)用于水利、土木、道路等工程領(lǐng)域中，其性能的優(yōu)劣很大程度會影響到工程整體的安全性，因此有關(guān)混凝土力學(xué)性能的研究一直被人們所重視。傳統(tǒng)的混凝土宏觀模型一般在研究過程中會忽略混凝土內(nèi)部的復(fù)雜結(jié)構(gòu)，將其視為均質(zhì)且各向同性的連續(xù)性材料，這種模型并不能準(zhǔn)確描述混凝土在受到外荷載時其內(nèi)部從微裂紋的萌生、擴(kuò)展直至最終產(chǎn)生宏觀裂紋的全過程。

隨著計算機(jī)技術(shù)的不斷發(fā)展，基于細(xì)觀層次的數(shù)值模擬成為一種研究混凝土力學(xué)性能的新途徑[1]基于細(xì)觀尺度上對不同材料組分提出的不同假設(shè)，學(xué)者們提出了不同的細(xì)觀力學(xué)模型及研究方法，目前較為常用的細(xì)觀模型有格構(gòu)模型[2]、M-H模型[3]、隨機(jī)骨料模型[4]、隨機(jī)力學(xué)模型[5]、內(nèi)聚力模型[6]等。而近年來在有關(guān)混凝土裂紋擴(kuò)展方面的研究中，基于內(nèi)聚力模型的研究方法更加被人們所重點關(guān)注，該模型將混凝土視為由骨料、砂漿以及骨料-砂漿間的粘結(jié)面組成的多相型復(fù)合材料，主要通過在相鄰實體單元間插入內(nèi)聚力單元來模擬混凝土內(nèi)部的損傷開裂，采用牽引力分離準(zhǔn)則來描述內(nèi)聚力單元的損傷和演化過程，該模型的優(yōu)點在于不僅可以避免裂紋尖端的應(yīng)力奇異性，同時支持多條裂紋的擴(kuò)展，更加符合實際情況。

本文基于MATLAB軟件對前人隨機(jī)骨料投放算法中單個骨料面積計算及骨料間交疊判斷方法加以改進(jìn)，編寫了更加高效簡便的骨料生成與投放程序，同時通過自編C++程序，將零厚度內(nèi)聚力單元批量插入到實體單元邊界處，并選取合理的本構(gòu)關(guān)系，建立了二級配混凝土內(nèi)聚力模型，之后以物理實驗結(jié)果中抗拉強(qiáng)度和抗拉彈性模量為指標(biāo)，對模型在單軸拉伸情況下開展數(shù)值模擬研究，從細(xì)觀角度分析混凝土在單軸受拉時其內(nèi)部從起始損傷直至破壞的全過程。從而對本文提出的混凝土模型生成方法的可行性與有效性加以驗證。

1 內(nèi)聚力模型的生成方法

1.1 骨料的生成與投放

在實際工程中所使用的粗骨料一般是通過人工破碎生成的，大多數(shù)呈現(xiàn)為凸多面體的形狀，通過X射線對混凝土進(jìn)行掃描也可以發(fā)現(xiàn)[7-8]，混凝土在任意截面內(nèi)的骨料形狀大致呈凸多邊形。因此對于二維的數(shù)值模擬來講，采用凸多邊形的骨料相較圓形骨料能更加真實地反應(yīng)混凝土的細(xì)觀結(jié)構(gòu)。

有關(guān)骨料的投放算法目前已有很多研究，馬懷發(fā)等[9]提出將骨料形狀按圓形考慮并全部投放，之后分別對每個骨料在圓周上取若干點并按逆時針依次連接形成“基骨料”，再以每個“基骨料”的較長邊為直徑向外延拓生成多邊形骨料，該方法雖然可以很好地控制骨料的投放位置及顆粒數(shù)，但所延拓的多邊形骨料的粒徑可能會超出該級配粒徑的范圍；任志剛等[10]對該方法提出了改進(jìn)，對“基骨料”延拓產(chǎn)生的新節(jié)點到圓心的距離不能超出該級配的粒徑范圍且滿足“凸”形的條件；汪奔等[11]通過將生成的骨料進(jìn)行適當(dāng)?shù)目s放也達(dá)到了類似的效果。以上學(xué)者的算法雖能滿足骨料的“凸”性及各級配的占比要求，但缺點是算法均較為復(fù)雜，對于編程初學(xué)者不太友好。基于此，本文對以上學(xué)者的算法加以借鑒和改進(jìn)，編制了一套用于二維隨機(jī)凸多邊形骨料的生成與投放算法，該方法不僅可以滿足粗骨料的占比及級配分布的合理性，同時也更加簡單易行、便于編程實現(xiàn)。

1.1.1 隨機(jī)凸多邊形骨料的生成

本文將混凝土視為由骨料、砂漿及骨料-砂漿界面過渡區(qū)組成的三相復(fù)合性材料，其中直徑D<5 mm的骨料視作砂漿考慮。首先利用瓦拉文公式[12]建立空間內(nèi)骨料占比與二維某一截面骨料面積之間的關(guān)系，使模型可應(yīng)用于二維平面狀態(tài)的分析中，平面內(nèi)任意一點處含有骨料顆粒直徑D

(1)

式中，Pc為不同級配骨料在二維某截面內(nèi)出現(xiàn)的概率；Pk為骨料占混凝土截面總面積的百分比；Dmax為最大骨料的顆粒直徑；D0為某一級配骨料的顆粒直徑。

對于任意一個隨機(jī)凸多邊形骨料來講，其大小和形狀都是隨機(jī)的，本文通過在極坐標(biāo)下生成隨機(jī)的極半徑R以及極角alpha的方法確定凸多邊形骨料的每個頂點坐標(biāo)，之后將其轉(zhuǎn)換為直角坐標(biāo)，其表達(dá)式為

(2)

式中，R為骨料第i個頂點的極半徑；r為某級配下骨料的代表粒徑；alpha為骨料第i個頂點的極角；n為骨料的邊數(shù)，本文取5～12之間的任意整數(shù)；i為骨料的第i個頂點。

1.1.2 隨機(jī)凸多邊形骨料的投放

在完成單個骨料算法的編制后，即可開始投放，投放之前需要確定二維混凝土模型的長和寬、骨料在各個級配的代表粒徑及相應(yīng)占比，由瓦拉文公式即可計算出各個級配骨料的占比，并將其按模型的尺寸轉(zhuǎn)化為面積。具體的投放步驟為

(1)根據(jù)極坐標(biāo)的方法生成各個級配的骨料庫，并記錄下每個骨料的頂點及相應(yīng)外圍節(jié)點的坐標(biāo)和骨料的面積。

(2)按照骨料粒徑的大小在指定區(qū)域內(nèi)逐個投放骨料，并判斷新生成的骨料與之前的骨料是否產(chǎn)生重疊，若重疊則重新投放新的骨料，反之記錄下相應(yīng)的骨料頂點坐標(biāo)，并將該骨料的面積累加至已完成投放的骨料面積中。

(3)判斷各個級配已投放的骨料面積是否達(dá)到或接近所計算的值，由于每次投放很難達(dá)到完全精確，故設(shè)置了3%的容錯率，若未達(dá)到預(yù)期值，則繼續(xù)投放，反之則結(jié)束該級配骨料的投放，進(jìn)行下一級配，直至完成所有級配骨料的投放。

本文隨機(jī)凸多邊形骨料的生成與投放的流程見圖1。

圖1 骨料生效與投放流程

本文借鑒馬懷發(fā)[9]、任志剛[10]、汪奔[11]等學(xué)者有關(guān)隨機(jī)骨料投放算法的研究，并對算法加以改進(jìn)。改進(jìn)有：①對單個骨料面積計算的方法進(jìn)行了改進(jìn)，采用MATLAB中的polyarea(x，y)函數(shù)，將單個骨料的頂點坐標(biāo)按逆時針方向依次存儲在所創(chuàng)建的x和y向量中，然后利用該函數(shù)即可返回該骨料的面積，該方法相較馬懷發(fā)、任志剛等將多邊形按頂點分割成三角形累加計算的方法更加簡便且更加高效，也節(jié)省了編制相應(yīng)代碼的時間。②對骨料的交疊判斷采取了不同的算法，如圖2所示，首先對每個骨料頂點坐標(biāo)擴(kuò)大1.1倍并按逆時針方向依次將相鄰兩點之間均分成8等分，形成一圈骨料外圍節(jié)點，然后利用inpolygon(x，y，xv，yv)函數(shù)判斷每個新生成的骨料外圍節(jié)點坐標(biāo)集合(x，y)是否在每個已生成的多邊形骨料外圍節(jié)點坐標(biāo)集合(xv，yv)的邊界內(nèi)，如果函數(shù)返回值為1，則骨料頂點產(chǎn)生交疊，投放失敗，如果不在邊界內(nèi)，則表明該骨料投放成功，并記錄下相應(yīng)的頂點坐標(biāo)。該方法運行效率高且代碼易于編制，僅需記錄下各骨料的頂點及外圍節(jié)點坐標(biāo)即可完成交疊判斷。③取消將投放骨料顆粒個數(shù)作為最終某一級配下所有骨料投放成功的判定條件，對于隨機(jī)多邊形骨料來講，每個骨料的面積都是不等的，如果將骨料的個數(shù)列為判定條件，不僅會增加判定次數(shù)，降低運行效率，同時也變相加入了人為的干擾因素。

圖2 骨料交疊判斷示意

1.2 內(nèi)聚力單元的批量插入

雖然通常情況下混凝土內(nèi)部的開裂主要產(chǎn)生于骨料-砂漿交界面、水泥砂漿內(nèi)部，但仍有部分骨料可能會產(chǎn)生斷裂，為了更加精確的描述混凝土內(nèi)部的損傷開裂，將骨料內(nèi)部也視為潛在開裂區(qū)。本文通過自編的C++算法對inp文件進(jìn)行批處理，實現(xiàn)了零厚度內(nèi)聚力單元精確且高效的批量插入。具體的算法步驟為

(1)通過上述算法建立隨機(jī)骨料模型，并進(jìn)行網(wǎng)格劃分，本文采用自由劃分網(wǎng)格法，對骨料和砂漿區(qū)域分別創(chuàng)建集合，然后導(dǎo)出該模型的inp文件。

(2)讀取inp文件中所有的節(jié)點及單元信息，并按照骨料與砂漿建立相應(yīng)集合進(jìn)行存儲，之后分別對各單元集合內(nèi)的節(jié)點進(jìn)行遍歷，記錄下每個節(jié)點出現(xiàn)的次數(shù)，每個節(jié)點重復(fù)出現(xiàn)了幾次就表明需要分裂幾次，并根據(jù)原節(jié)點信息在相同的位置生成新的節(jié)點，如圖3所示，節(jié)點3由于被6個單元所公用，因此需要將其分裂5次。

圖3 算法示意

(3)將新產(chǎn)生的節(jié)點重新分配給相應(yīng)的實體單元，并在各實體單元間創(chuàng)建零厚度的內(nèi)聚力單元，并分別創(chuàng)建相應(yīng)的集合，如圖3所示，綠色為骨料-砂漿界面單元，橙色為骨料內(nèi)部單元，紫色為砂漿內(nèi)部單元。

(4)將各集合內(nèi)新的節(jié)點、單元信息重新覆蓋之前的inp文件，并導(dǎo)入ABAQUS，完成零厚度內(nèi)聚力單元的插入。

以建立的二級配混凝土模型(尺寸為150 mm×150 mm，骨料占比約為46.9%)為例，根據(jù)本文的骨料生成與投放算法及內(nèi)聚力單元的批量插入算法所建模型如圖4所示，分別按集合顯示了各單元分布情況。圖中模型共37 614個節(jié)點，由18 657個三角形平面應(yīng)力單元(CPS3)以及12 538個內(nèi)聚力單元組成，其中采用本文C++程序插設(shè)內(nèi)聚力單元僅需8 min左右，說明本文編制算法具有高效性和可行性。

圖4 二級配混凝土模型

2 混凝土單軸受拉破壞細(xì)觀模擬

2.1 模型的生成

為驗證本研究提出的內(nèi)聚力模型生成算法的可行性和有效性，利用王仲華等[13]以灰?guī)r作為粗、細(xì)骨料時所開展的二級配混凝土(28 d齡期)單軸拉伸試驗結(jié)果作為數(shù)值模擬的對象，試驗所得混凝土的抗拉強(qiáng)度平均值為2.36 MPa，抗拉彈性模量為37.7 GPa。

本試驗采用開發(fā)的骨料生成與投放算法及內(nèi)聚力單元的批量插入算法，生成150 mm×150 mm的二級配混凝土內(nèi)聚力模型，選用全局尺寸為2 mm的三角形網(wǎng)格進(jìn)行單元劃分，荷載施加方式為位移勻速加載，對試件的底邊施加豎向的位移約束，同時最中間節(jié)點施加水平位移約束，圖5為模型荷載施加示意。

圖5 二級配混凝土細(xì)觀模型

2.2 本構(gòu)關(guān)系的選取

本文模型中，骨料、砂漿單元采用三角形平面應(yīng)力實體單元(CPS3)劃分網(wǎng)格，為提高計算效率，采用線彈性的本構(gòu)關(guān)系。骨料-砂漿界面單元、骨料內(nèi)部單元、砂漿內(nèi)部單元采用內(nèi)聚力單元(COH2D4)模擬，內(nèi)聚力單元描述裂縫尖端的斷裂方法為：假設(shè)在真實裂紋尖端之前存在一個“斷裂過程區(qū)”，該區(qū)域裂尖的張開位移與應(yīng)力的關(guān)系采用牽引力分離準(zhǔn)則來表示。目前常用的牽引力分離準(zhǔn)則有雙線型、梯型、指數(shù)型[14]等，雙線型牽引力分離準(zhǔn)則適用于描述脆性材料的斷裂，因此本文采用該準(zhǔn)則來描述混凝土的斷裂。

目前ABAQUS中共有4種損傷起裂準(zhǔn)則：最大名義應(yīng)力準(zhǔn)則、最大名義應(yīng)變準(zhǔn)則、平方名義應(yīng)力準(zhǔn)則和平方名義應(yīng)變準(zhǔn)則?；炷猎趩屋S受壓條件下的破壞可能受到法向和切向應(yīng)力的共同作用，同時存在I型(張拉)和II型(剪切)兩種斷裂模式，因此本文考慮了兩種裂紋模式的組合，采用了平方名義應(yīng)力準(zhǔn)則來判定單元的損傷起裂，即

(3)

采取了基于能量演化準(zhǔn)則中的Benzeggagh-Kenane(B-K)混合斷裂準(zhǔn)則來描述單元的損傷演化，該準(zhǔn)則適用于包含I型、II型裂紋的混合型斷裂，表達(dá)式如下：

(4)

2.3 參數(shù)校調(diào)

一般混凝土內(nèi)部各相組分的力學(xué)性能關(guān)系為：骨料>水泥砂漿>骨料-砂漿交界面[15]，因此本著這一原則對砂漿內(nèi)部及骨料-砂漿粘結(jié)面單元的參數(shù)進(jìn)行校調(diào)，通過不斷修正單元材料參數(shù)直至與物理試驗結(jié)果相吻合。為使結(jié)果更加精確，共建立3個150 mm×150 mm的二級配混凝土試件，求解后取其結(jié)果的平均值，最終確定的材料參數(shù)見表1。

表1 材料參數(shù)

對3個150 mm×150 mm的二級配混凝土模擬試件分別進(jìn)行求解，得到圖6的單軸拉伸應(yīng)力應(yīng)變曲線，具體計算結(jié)果為：二級配混凝土中3個試件的抗拉強(qiáng)度最大值分別為2.26、2.31、2.35 MPa，抗拉強(qiáng)度平均峰值應(yīng)力為2.29 MPa；抗拉彈性模量為35.38、35.12、35.11 GPa，平均值為35.20 GPa，與王仲華[13]的試驗數(shù)據(jù)對比，抗拉強(qiáng)度最大誤差不超過2.8%，抗拉彈性模量最大誤差不超過6.6%。

圖6 單軸拉伸應(yīng)力應(yīng)變曲線

2.4 單軸受拉破壞全過程模擬

以試件1為例，對混凝土在單軸受拉時的全過程進(jìn)行分析，圖7給出了試件的單軸受壓應(yīng)力-應(yīng)變曲線，A、B、C、D等4個點分別對應(yīng)圖8中4個具有代表性的損傷演化階段。為了清楚的顯示不同類型單元的損傷破壞，用紅色代表骨料-砂漿粘結(jié)面單元，綠色代表砂漿內(nèi)部單元，黃色代表骨料內(nèi)部單元，c表示試件在受拉方向所放大的倍數(shù)。

圖7 試件1單軸拉伸應(yīng)力應(yīng)變曲線

圖8 試件1損傷演化示意

(1)微裂紋成核階段。此過程發(fā)生于應(yīng)力上升階段，該階段混凝土試件已經(jīng)出現(xiàn)了塑性變形，由于骨料與砂漿單元在力學(xué)性能上差異較大，在其交界面處容易產(chǎn)生應(yīng)力集中，導(dǎo)致骨料-砂漿粘結(jié)面上已經(jīng)產(chǎn)生損傷，在骨料表面生成了一些與受拉方向接近垂直的微裂紋。

(2)微裂紋擴(kuò)展階段。隨著混凝土在受拉方向的應(yīng)變不斷增大，這些微裂紋向相鄰骨料的方向逐漸擴(kuò)展，即擴(kuò)展至了砂漿內(nèi)部單元區(qū)域，由圖8b可以發(fā)現(xiàn)，在微裂紋的擴(kuò)展過程中，相鄰粒徑較大的骨料之間更容易產(chǎn)生連接貫通，這可能是由于隨著骨料粒徑的增大，其表面潛在的微裂紋會越多，因此相鄰兩顆大粒徑骨料之間的砂漿內(nèi)部單元更傾向于連接貫通，有逐漸形成裂紋貫通帶的趨勢。

(3)微裂紋連接貫通階段。此時微裂紋逐漸從骨料表面及其附近區(qū)域擴(kuò)展貫通，在混凝土內(nèi)部大致形成了一條主裂紋帶，由圖8c可知，此時的骨料-砂漿粘結(jié)面單元及砂漿內(nèi)部單元仍然存在，即代表宏觀裂紋還沒有真正產(chǎn)生，在主裂紋帶逐漸形成的過程中，部分區(qū)域的微裂紋由于主裂紋帶的擴(kuò)展而受到限制，甚至?xí)S著主裂紋帶的繼續(xù)增加而逐漸呈現(xiàn)閉合的趨勢。

(4)宏觀裂紋形成階段：該階段大部分骨料-砂漿粘結(jié)面單元及砂漿內(nèi)部單元已經(jīng)完全損傷而被刪除，最終形成了一條肉眼可見的宏觀裂紋帶，裂紋帶近似垂直于受拉方向，但這條裂紋帶并非是一條連續(xù)的曲線，這可能是由于粗骨料的形狀影響了裂紋的擴(kuò)展路徑所導(dǎo)致的。在整個單軸拉伸模擬過程中，骨料內(nèi)部并未發(fā)現(xiàn)黃色單元，表明骨料內(nèi)部單元沒有產(chǎn)生損傷。

3 結(jié) 論

(1)基于MATLAB對前人隨機(jī)骨料投放算法進(jìn)行改進(jìn)，優(yōu)化了單個骨料面積計算及骨料間交疊判斷的方法，編寫了更加高效簡便的骨料生成與投放算法，之后采用C++編程，實現(xiàn)在相鄰實體單元間快速、成批次地插入內(nèi)聚力單元，并根據(jù)單元所處的部位進(jìn)行集合劃分，完成內(nèi)聚力模型的建立。

(2)以王仲華所做的物理試驗結(jié)果中抗拉強(qiáng)度和抗拉彈性模量為指標(biāo)，對本文建立的內(nèi)聚力模型進(jìn)行校調(diào)，與試驗數(shù)據(jù)對比，抗拉強(qiáng)度最大誤差不超過2.8%，抗拉彈性模量最大誤差不超過6.6%。

(3)從單軸拉伸模擬結(jié)果來看，混凝土在均勻受拉時大致經(jīng)歷了4個階段：微裂紋成核階段、微裂紋擴(kuò)展階段、微裂紋連接貫通階段和宏觀裂紋形成階段，這與實際中混凝土單軸受拉全過程大致相似，從而驗證了本研究提出的內(nèi)聚力模型生成算法的可行性和有效性。