基于粘結單元的三維隨機細觀混凝土離散斷裂模擬

楊貞軍，黃宇劼，堯 鋒，劉國華

(1. 武漢大學土木建筑工程學院，湖北省巖土與結構安全重點實驗室，武漢，湖北 433000；2. 浙江大學建筑工程學院，杭州，浙江 310058)

傳統的宏觀混凝土模型假設均質材料，不能反映細觀組分空間分布的隨機性和非均質性引起的損傷與斷裂局部化，難以可靠預測材料的力學響應與裂縫擴展過程。因此進行混凝土細觀斷裂模擬研究十分必要[1 ? 5]。在三維細觀混凝土模擬中，假設骨料為球體是一種常用的簡化建模方法[4,6 ?7]，雖然可以高效地生成細觀結構，但無法反映骨料形態的影響。近年來，采用更為復雜的橢球和凸多面體來模擬骨料日趨常見[8 ? 9]。在應用上，復雜骨料模型大多用于幾何特征研究[8,10]、邊界效應分析[11]、滲透性能研究[12 ? 14]以及線彈性應力分析[15 ? 17]；在三維復雜非線性損傷和斷裂的模擬方面仍非常有限，且大多采用連續介質損傷塑性模型[18 ? 21]，難以定量預測裂縫寬度，因此有必要采用離散裂縫模型進行更精細的斷裂模擬。另一方面，在骨料幾何建模方法上，一般通過自編程序進行參數化建模[18,21]，或基于Voronoi 圖生成細觀結構[22]，或基于ANSYS 的ADPL 命令流建模[23]。但是，對三維復雜細觀結構進行網格劃分并不容易，較為直接的方法是將細觀結構映射到規則的立方體背景網格上[4, 18 ? 20, 24]，但要保證精度則會導致巨大的單元數量。ABAQUS 具有強大的前處理模塊，可以靈活有效地對復雜結構進行網格自動劃分，不過尚未見到將其用于三維混凝土細觀結構建模與網格劃分的研究。

本文首先采用Python 編寫程序，利用ABAQUS腳本接口[25]進行前處理二次開發，考慮骨料實際級配，有效地建立多面體隨機骨料幾何模型并進行有限元網格劃分。然后通過自編C++程序在骨料-砂漿界面上以及砂漿中高效插設零厚度的離散粘結裂縫單元[26 ? 28]，模擬受拉試件復雜細觀三維多裂縫的起裂與擴展直至破壞的全過程，并對粘性裂縫單元的主要材料參數(抗拉強度與斷裂能)對荷載-位移曲線、斷裂過程、裂縫面特征的影響進行詳細分析，以加深對混凝土細觀斷裂機理的認識。應該指出，筆者論文[26 ? 28]也采用了插設界面單元的方法，但均非基于本文提出的隨機骨料模型；雖然文獻[26 ? 27]采用隨機場來間接表征混凝土的空間非均質性，但和文獻[28]一樣，僅提出了均質材料中插設二維、三維界面的方法。本文是前述工作的擴展，并將Matlab 程序改進為C++，大大提高了效率，更適合大規模三維有限元網格的處理。

1 基本算法

1.1 隨機多面體骨料建模方法

算法分兩步，第一步建立幾何模型，第二步建立有限元模型。

首先編寫MATLAB 程序進行多面體隨機骨料的生成與投放，建立幾何模型。采用基于球體的方式產生凸多面體，即凸多面體的各頂點位于球面上。該球體稱為基球，其球面上隨機產生8 個～25 個點作為多面體的頂點，頂點i 的笛卡爾坐標用Euler 空間的極角αi和方位角βi表示為：

式中：(x0, y0, z0)是球心的坐標；r0是球半徑；αi在[0, 2π]、βi在[0, π]范圍內按均勻分布隨機產生。然后采用MATLAB 算法“convhulln”基于頂點建立三角形面，這些面構成骨料多面體。為避免頂點距離太近而造成三角形面畸形，需保證球面上的頂點間距不小于ξr0，本文取ξ=0.5。多面體的生成較為耗時；為提高效率，本文不采用邊投放、邊生成骨料的方式，而是通過生成足夠數量(本文為20000)的直徑為單位一的多面體，預先建立骨料形態數據庫，圖1 為一些骨料例子。

圖1 不同形態的多面體骨料Fig.1 Polyhedral aggregates of different shapes

針對給定模擬區域，具體算法如下：

1) 根據所需骨料含量，采用實際骨料級配表1求解各級配的骨料體積，從骨料直徑最大的級配開始分級配生成并投放骨料；

表1 三級配骨料尺寸分布[29]Table1 Three-graded aggregate size distribution[29]

2) 對每個級配的操作步驟：隨機產生一個球心坐標作為投放位置，從骨料數據庫中隨機選取一個骨料，并在此級配區間內按均勻分布隨機賦予其直徑；

3) 若該骨料不與已有骨料相交或重疊，則該骨料投放成功，進行下一個投放，否則重新進行本次投放；

4) 本級配骨料體積達到后，進行下一級配投放；

5) 循環2)步～4)步，直至達到總骨料含量。

在上述算法第3)步中，新生成骨料與已有骨料之間可能存在如圖2 所示的3 種位置關系：不相交也不重疊、相交、重疊(包含)。

圖2 兩個骨料之間的三種位置關系Fig.2 Spatial relations of two aggregates

通過下述算法判斷兩個多面體相交或重疊：如圖2(a)所示，對于新生成的多面體A 和已有多面體B，要使二者不相交且不重疊，則B 的各面需在A 各面的同一側。A 的頂點A1～A3 所構成平面A123 的方程是：

式中：(xAi, yAi, zAi)是頂點Ai(i=1～3)的坐標，則對于平面A123 外的點(x, y, z)，有F(x, y, z) > 0 或F(x, y, z)<0。則上述位置判斷問題可以等效為：遍歷A 的各面，如果存在一個面(例如A123)，使得A 的內部各點(凸多面體可只選取其形心)與B 的任意頂點(xBj, yBj, zBj)位于面A123 的不同側，則可判定A 與B 既不相交也不重疊。相比于遍歷A、B 的頂點以保證各自頂點不在另一個多面體內部的算法，上述位置判斷算法更為直觀和有效。若A 的形心坐標為(xA0, yA0, zA0)，A 與B 既不相交也不重疊時只需滿足下式：

骨料與骨料、骨料與邊界之間預設最小間距(本文設置為較大骨料直徑的0.05 倍)以方便網格生成。

第二步，編寫Python 腳本程序，進行ABAQUS 前處理二次開發，建立多面體骨料并完成模型網格劃分。具體算法如下：

1) 建立模型數據庫；

2) 生成混凝土試件立方體部件(part)；

3) 讀取隨機骨料幾何模型中所有骨料的信息：基球的球心坐標、頂點坐標、各面的頂點構成；

4) 對于每個骨料，生成三維球體部件；

5) 延拓骨料的各面，切割該球體，過程見圖3；

6) 在assembly 模塊中創建部件實例(instance)，采用布爾操作合并生成的多面體與試件立方體；

7) 重復上述步驟，直至創建所有骨料，見圖4；

圖3 在ABAQUS 中通過Python 腳本程序生成多面體Fig.3 Polyhedra generation using Python scripts in ABAQUS

圖4 生成的50 mm 試件多面體骨料有限元模型Fig.4 FEM model of a specimen with polyhedral aggregates

8) 布置全局網格種子，采用“自由網格”技術用四面體單元進行網格劃分。

以邊長為50 mm 立方體試件和骨料含量30%為例，上述算法建立的骨料有限元模型見圖4。該模型含有423 個骨料，網格單元平均尺寸為2 mm。

1.2 插設三維粘結裂縫單元

本文將ABAQUS 零厚度三維六節點粘結裂縫單元COH3D6 插設到砂漿的有限元實體單元之間，以及骨料-砂漿界面上，用于模擬混凝土復雜非線性斷裂過程。采用損傷系數D [0, 1]描述裂縫單元達到強度后的損傷程度，反映剛度退化，該系數是裂縫單元有效相對位移δm的函數：

式中：δn、δt、δs分別是法向、切向的相對位移；< >為Macaulay 括號，表示為：

以線性軟化準則為例，損傷系數可以寫成：

式中：δm,max是加載歷史中的最大有效相對位移；δm0和δmf分別是裂縫起裂和完全破壞時的有效相對位移。

法向剛度kn與切向剛度ks、kt分別表示為：

相應的應力則可以表示為：

采用名義應力平方準則作為起裂準則:

式中，tn、ts和tt分別為法向、切向的應力。

文獻[26 ? 27]報道了在二維、三維實體有限元單元網格中插設粘結單元的算法和MATLAB 程序實施，用于均質材料。本文將該算法加以拓展，通過編寫C++程序，在混凝土細觀結構中插設三維粘結單元，圖5 給出了算法示意圖。假設骨料不會開裂，骨料內無需插設粘結單元，因此只在骨料-砂漿界面以及砂漿中采用粘結單元(圖5(c))，分別標記為CIE_INT 與CIE_CEM。

圖5 細觀混凝土三維粘結單元的插設算法Fig.5 3D inserting algorithm of cohesive elements

具體插設算法如下：

1) 讀取初始網格，圖5(a)顯示了網格中與骨料表面相連的2 個砂漿四面體單元A 與B，以及1 個與砂漿單元B 相連的砂漿單元C。

2) 如圖5(b)所示，搜尋單元的面，骨料的面存入結構體FACE_INT 中，砂漿的面存入FACE_CEM 中。

3) 如圖5(b)所示，將節點按其所在單元的數量拆開(即不改變坐標地進行節點的原位復制)并重新編號，再分配給各單元。

4) 如圖5(c)所示，根據FACE_INT、FACE_CEM 分別建立粘結單元CIE_INT、CIE_CEM，從而完成粘結單元的插設，圖6 給出一個骨料表面插設的零厚度粘結單元CIE_INT 作為示例。

5) 輸出插設粘結單元后的整體網格信息，作為ABAQUS 可讀取的*.inp 輸入文件供計算使用。

對于圖4 中約有3.2 萬個節點、18.1 萬個四面體單元(C3D4)以及36.7 萬個單元面的原始模型，采用原有MATLAB 程序[27 ? 28]進行插設耗時5 h，而本文C++程序不到3 min。插設粘結單元之后，模型大約有53.1 萬個節點數、27.0 萬個粘結單元，其中骨料-砂漿界面上大約有2.6 萬個粘結單元。

圖6 骨料表面的粘結單元(CIE_INT)Fig.6 Cohesive element on the surface of an aggregate

2 單軸受拉模擬

骨料和砂漿假設為線彈性材料，粘結單元采用線性受拉/受剪軟化關系和名義應力平方起裂準則進行模擬。采用文獻[22]中的材料參數(表2)，并假設受剪與受拉斷裂參數相同(偏于保守)。粘結單元的初始彈性剛度應足夠高以模擬未開裂初始狀態，但剛度過高會導致算法不收斂。經過多次嘗試，將其取為6×106MPa/mm。

表2 材料參數[22]Table2 Material parameters[22]

對生成的50 mm 的細觀混凝土立方體試件(圖4)進行單軸受拉模擬。邊界條件如圖7 所示，左端面全約束，右端面受均勻分布水平(x 向)位移荷載控制，最終位移0.06 mm。使用ABAQUS/Explicit 顯式求解器進行準靜態計算分析，設定荷載步時間為0.005 s。

圖7 單軸受拉混凝土試件的邊界條件Fig.7 Boundary conditions under uniaxial tension

圖8 給出了該模型采用兩組不同斷裂能時的平均應力-位移曲線，其中一組采用表2 的參數(即參照組)，另一組Gc0.06Gi0.01 表示界面粘結單元的斷裂能為0.01 N/m，其余參數同表2。圖中也給出了典型單軸受拉試驗數據[30]作為對比，該試驗試件的尺寸是50 mm×50 mm×60 mm，與模型比較接近。可見Gc0.06Gi0.01 得到的模擬結果與試驗吻合較好，表明模型能夠有效預測混凝土力學響應。但由于試驗試件與本文數值模型在細觀組分方面并不相同，該試驗數據僅作參考。

圖8 單軸受拉混凝土試件的宏觀應力-位移曲線Fig.8 Average stress - displacement curves of concrete specimen under uniaxial tension

2.1 開裂過程分析

圖9 顯示了采用參考材料參數(表2)模擬的試件正面和背面視角的開裂過程(對應于圖8 曲線的A 點～F 點)。其中紅色的單元為損傷因子D≥0.9的粘結單元，在本文中假設為已完全損壞的離散裂縫。為了觀察裂縫表面形態，圖10 顯示了去掉這些損壞裂縫單元后的相應變形過程。為了方便觀察，變形放大系數均設為200。

圖9 單軸受拉開裂過程：粘結裂縫單元的受拉變形Fig.9 Fracture process under uniaxial tension: the red cohesive elements with D≥0.9

如圖9 和圖10 所示，在早期加載階段(A 點之前)，由于界面粘結單元的斷裂參數比砂漿低，骨料-砂漿界面上出現大量微裂縫，但其損傷因子D 仍較低，反映在裂縫單元的張開位移較小，以及荷載-位移曲線在加載位移較小時即顯示出非線性特征；隨著加載位移的增加直至試件達到抗拉強度(A 點)，微裂縫寬度緩慢發展，但裂縫并不顯著；當試件進入軟化段時(B 點～D 點)，一些微裂縫的寬度迅速增加進而局部化形成主裂縫，并與砂漿中新生成的裂縫連通，而另外一些微裂縫則逐漸卸載閉合；如圖中階段(E 點～F 點)所示，隨著位移進一步增加，連通的裂縫迅速張開變寬，形成一條貫穿主裂縫，試件拉裂。這和混凝土單向受拉實驗[30]和其他數值模擬[22,28]的斷裂擴展過程一致。

2.2 材料斷裂參數分析

現有研究表明[26 ? 27]，粘結裂縫單元的材料參數，特別是抗拉強度與斷裂能，對模擬結果影響較大。本文對砂漿中與骨料-砂漿界面上的粘結單元選取不同組別的抗拉強度與斷裂能進行研究。保持斷裂能不變(見表2)，粘結單元取不同抗拉強度時，使用TcTi 的編號規則，如Tc4Ti2 表示砂漿和界面粘結單元的抗拉強度分別為4 MPa 和2 MPa；保持表2 中抗拉強度不變，取不同斷裂能，使用GcGi 的編號規則，如Gc0.06Gi0.03 表示砂漿和界面粘結單元的斷裂能分別為0.06 N/mm 和0.03 N/mm。將采用表2 材料參數的模型設置為參考組(編號為Tc4Ti4_REF 與Gc0.06Gi0.03_REF)。

圖10 單軸受拉開裂過程Fig.10 Crack propagation process under uniaxial tension

取8 組粘結單元材料參數進行模擬，對每一組，只變動TcTiGcGi 參數中的一個，其余參數與參考組相同。得到如圖11 和圖12 所示試件宏觀平均應力-位移曲線，以及如圖13 和圖14 所示的試件破壞形式。

圖11 粘結裂縫單元的抗拉強度對宏觀應力-位移曲線的影響Fig.11 Average stress - displacement curves: effects of fracture energy of cohesive elements

圖12 粘結裂縫單元的斷裂能對宏觀應力-位移曲線的影響Fig.12 Average stress - displacement curves: effects of fracture energy of cohesive elements

圖13 粘結裂縫單元的抗拉強度對試件破壞形式的影響Fig.13 Effects of tensile strength of cohesive elements on the fracture surfaces

由圖11 可見，粘結單元的抗拉強度對試件應力-位移曲線峰值以及峰值前非線性段影響顯著：隨著粘結單元強度的增加，試件強度增加；當骨料-砂漿界面裂縫單元強度降低而砂漿裂縫單元強度增加時，試件強度亦增加，表明砂漿粘結單元的抗拉強度對試件承載力起控制作用。圖11 也表明：隨著試件強度的增加，其脆性也逐漸增加。由圖12 可見，粘結單元的斷裂能對試件應力-位移曲線軟化段影響顯著，隨著斷裂能的提高，試件的延性增加，但對試件強度以及曲線峰前段影響不大。

圖14 粘結裂縫單元的斷裂能對試件破壞形式的影響Fig.14 Effects of fracture energy on crack surfaces

由圖13 可見，粘結單元抗拉強度對裂縫面形態有較大影響，該影響主要體現于砂漿、界面粘結單元的強度相對比值γT=Tc/Ti。當γT≠1 時，斷裂面的位置基本相同，γT主要影響斷裂面的形態；當γT=4 時(圖13(c)設計組Tc4Ti1)，斷裂面上出現的大骨料數量最多，其次是γT=3(圖13(b)設計組Tc6Ti2)，然后是γT=2(圖13 參考組Tc4Ti2)。由于本文粘結單元以名義拉應力作為起裂準則，隨著砂漿、界面粘結單元強度差別的增大，界面愈發成為薄弱環節。由于大骨料表面的界面單元面積較大，較容易成為裂縫起裂與擴展通道，使試件傾向于在大骨料附近形成裂縫面。另外，當γT=1 時，如圖13(a)和圖13(d)所示，試件斷裂面較平滑且以小骨料數量居多；其中設計組Tc2Ti2的粘結單元強度較低，離加載端最近一層裂縫單元因應力集中而破壞。

圖14 顯示砂漿、界面粘結單元的斷裂能相對比值γG=Gc/Gi對裂縫形態的影響。當γG≤3 時(圖14(a)、圖14(b)、圖14(d)和參考組)，試件破壞時的位移和斷裂面形態基本一致；隨著斷裂能的提高，裂縫面略趨于曲折；但當γG=6 時(圖14(c)設計組Gc0.06Gi0.01)，砂漿裂縫單元的斷裂能顯著高于界面單元，界面成為薄弱環節，由于大骨料表面的界面單元面積較大，較容易成為裂縫擴展通道因而影響了裂縫擴展路徑，使最終斷裂面的位置與形態與其他組差異較大。

進一步將Gc0.06Gi0.01 與參考組Gc0.06Gi0.03_REF 比較，分析斷裂能、宏觀應力-位移曲線以及開裂過程之間的關聯。如圖15、圖16 所示，在這兩個模型的應力-位移曲線的軟化段中選取4 個應力水平接近的階段(A 點～D 點)進行分析。由圖可見，由于Gc0.06Gi0.01 的裂縫沿著圖14(c)黑色箭頭所指的大骨料(大骨料的界面作為薄弱環節)擴展，導致其裂縫面的形成比Gc0.06Gi0.03_REF 更為曲折復雜，因此其應力-位移曲線的軟化段由于斷裂耗能的需要變得較為平緩，而沒有呈現出Gi減小導致的脆性特征。由此可見，混凝土的力學響應能夠反映其裂縫發展特征，二者既決定于斷裂材料參數，也受到骨料大小、形狀等細觀結構因素的影響，顯示出十分復雜的破壞機理。

圖15 斷裂能對開裂過程的影響Fig.15 Effects of fracture energy on crack propagation processes

圖16 斷裂能對應力-位移曲線的影響Fig.16 Effects of fracture energy on stressdisplacement curves

3 結論

本文采用Python 編程，利用ABAQUS 腳本接口進行前處理二次開發，建立混凝土三維隨機多面體骨料的細觀有限元模型；通過高效插設離散粘結裂縫單元，成功模擬了單向受拉復雜三維多裂縫的起裂與擴展。主要結論如下：

(1) 在早期加載階段，由于界面粘結裂縫單元的斷裂參數比砂漿低，骨料-砂漿界面上出現大量微裂縫，荷載-位移曲線較早顯示出非線性特征；隨著位移的增加直至試件達到抗拉強度，微裂縫寬度緩慢增加；進入軟化段后，一些微裂縫寬度迅速增加，并與砂漿中新生成裂縫連通形成局部化的主裂縫，而另外一些微裂縫則逐漸卸載閉合。

(2) 宏觀應力-位移曲線主要受砂漿、界面粘性裂縫單元的抗拉強度和斷裂能絕對數值的影響：砂漿粘結裂縫單元的抗拉強度對試件承載力起控制作用，試件承載力隨著其強度的提高而增大；而斷裂能對非線性軟化段影響顯著，斷裂能的提高增加了試件的延性。

(3) 裂縫面的位置和形態主要受砂漿、界面粘性裂縫單元的抗拉強度、斷裂能相對比值的影響。比值較大時(γT>2 或者γG>3)，即界面粘結單元的力學性能遠弱于砂漿粘結單元時，由于大骨料表面的界面單元面積較大，較容易成為裂縫起裂與擴展通道，使試件傾向于在大骨料附近形成裂縫面。

(4) 混凝土的力學響應反映其裂縫發展特征，二者既決定于斷裂參數，也受到細觀結構的影響，本文建立的模型能夠有效地描述混凝土復雜三維斷裂過程。