基于預插內聚力單元的混凝土模型生成方法與單軸拉伸細觀模擬研究

馬 峰，田 羽，張小飛

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

0 引 言

混凝土由于具有良好的可塑性、耐久性及經濟性等優點，被廣泛應用于水利、土木、道路等工程領域中，其性能的優劣很大程度會影響到工程整體的安全性，因此有關混凝土力學性能的研究一直被人們所重視。傳統的混凝土宏觀模型一般在研究過程中會忽略混凝土內部的復雜結構，將其視為均質且各向同性的連續性材料，這種模型并不能準確描述混凝土在受到外荷載時其內部從微裂紋的萌生、擴展直至最終產生宏觀裂紋的全過程。

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

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

1 內聚力模型的生成方法

1.1 骨料的生成與投放

在實際工程中所使用的粗骨料一般是通過人工破碎生成的，大多數呈現為凸多面體的形狀，通過X射線對混凝土進行掃描也可以發現[7-8]，混凝土在任意截面內的骨料形狀大致呈凸多邊形。因此對于二維的數值模擬來講，采用凸多邊形的骨料相較圓形骨料能更加真實地反應混凝土的細觀結構。

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

1.1.1 隨機凸多邊形骨料的生成

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

(1)

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

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

(2)

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

1.1.2 隨機凸多邊形骨料的投放

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

(1)根據極坐標的方法生成各個級配的骨料庫，并記錄下每個骨料的頂點及相應外圍節點的坐標和骨料的面積。

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

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

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

圖1 骨料生效與投放流程

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

圖2 骨料交疊判斷示意

1.2 內聚力單元的批量插入

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

(1)通過上述算法建立隨機骨料模型，并進行網格劃分，本文采用自由劃分網格法，對骨料和砂漿區域分別創建集合，然后導出該模型的inp文件。

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

圖3 算法示意

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

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

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

圖4 二級配混凝土模型

2 混凝土單軸受拉破壞細觀模擬

2.1 模型的生成

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

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

圖5 二級配混凝土細觀模型

2.2 本構關系的選取

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

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

(3)

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

(4)

2.3 參數校調

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

表1 材料參數

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

圖6 單軸拉伸應力應變曲線

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

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

圖7 試件1單軸拉伸應力應變曲線

圖8 試件1損傷演化示意

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

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

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

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

3 結 論

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

(2)以王仲華所做的物理試驗結果中抗拉強度和抗拉彈性模量為指標，對本文建立的內聚力模型進行校調，與試驗數據對比，抗拉強度最大誤差不超過2.8%，抗拉彈性模量最大誤差不超過6.6%。

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