用細觀單元等效化方法模擬混凝土細觀破壞

杜修力，金 瀏

（北京工業(yè)大學 城市與工程安全減災(zāi)教育部重點實驗室，北京100124）

混凝土是典型的非均勻材料，其復(fù)雜的破壞機理及宏觀力學特性是研究者一直關(guān)注的問題。混凝土類非均質(zhì)材料的破壞過程實際上就是微裂紋萌生、擴展、貫通，直到最后宏觀裂紋產(chǎn)生導致混凝土失穩(wěn)破裂的過程。

對于非均勻混凝土材料問題，學術(shù)界進行了長期的探索，在近30年來已取得了很大的進展，提出了很多研究混凝土損傷斷裂過程的細觀力學模型［1］。Wang等［2－3］提出了隨機骨料模型，采用 非線性有限元技術(shù)，研究了混凝土試件的損傷斷裂全過程；Schlangen等［4］采用簡單格構(gòu)模型研究了混凝土等脆性材料的典型破壞機理；Mohamed等［5］提出了混凝土破壞機制研究的微觀結(jié)構(gòu)模型，考慮混凝土細觀組分力學性質(zhì)的隨機性，引入斷裂能概念，采用彌散裂紋模型描述單元受拉行為，對混凝土試件的破壞機理進行了研究；唐春安等［6］借助統(tǒng)計學及數(shù)值模擬方法，建立了混凝土損傷斷裂研究的隨機力學特性模型；Lopez等［7］采用“零厚度”界面單元法對混凝土試件在單軸拉伸下的斷裂破壞過程進行模擬，獲得的宏觀力學特性與試驗結(jié)果吻合良好；杜修力等［8］提出了高計算效率的細觀單元等效化模型，對混凝土材料的力學性能進行研究，數(shù)值結(jié)果與試驗結(jié)果吻合良好。這些研究工作，大多采用連續(xù)介質(zhì)力學與不可逆熱力學理論，導出相應(yīng)的連續(xù)損傷力學本構(gòu)方程和演化方程，進而研究混凝土斷裂破壞機理及宏觀力學性能。

筆者擬將細觀單元等效化模型在連續(xù)介質(zhì)力學領(lǐng)域的應(yīng)用延伸或擴展到非連續(xù)介質(zhì)力學范疇，從而更便于宏觀地認識混凝土材料內(nèi)部裂紋的擴展演化過程?；谠摾砟?，筆者首先簡要闡述細觀單元等效化模型及擴展有限元法的基本原理；其次，對細觀單元的等效力學性能作理論推導并進行數(shù)值驗證；最后采用擴展有限元法來研究混凝土材料的細觀斷裂破壞過程及其宏觀力學性能。

1 基本思想及原理

1.1 細觀單元等效化模型簡介

細觀單元等效化模型，從描述混凝土材料的細觀尺度入手，認為混凝土材料的宏觀非線性及尺度效應(yīng)根源于其內(nèi)部細觀組成的非均質(zhì)性，細觀組成的非均勻性導致宏觀力學行為的非線性。為反映混凝土材料非均勻性的本質(zhì)特征，文獻［9］基于概率統(tǒng)計學理論，提出以細觀單元材料力學特性（彈性模量）的變異系數(shù)作為材料非均勻性量度的指標，并給出了混凝土材料的特征單元尺度。采用特征單元尺度對混凝土細觀隨機骨料模型劃分網(wǎng)格后，得到的等效力學模型能準確地表征混凝土材料的非均勻特性。簡單的講，采用“材料特征單元尺度”來剖分有限元試件網(wǎng)格后獲得的等效力學系統(tǒng)，足以能夠準確獲得混凝土材料的宏觀力學性能。

如圖1所示的混凝土材料細觀單元等效化模型，首先在混凝土細觀隨機骨料模型的基礎(chǔ)上，采用“材料特征單元尺度”來剖分有限元網(wǎng)格并投影到已建立的隨機骨料模型上，得到的各細觀單元等效力學特性采用復(fù)合材料力學等效化方法來確定。非均質(zhì)混凝土材料力學體系最終等效成單元內(nèi)均勻且各向同性、單元間性質(zhì)各異的非線性力學模型。細觀單元等效化模型，本質(zhì)上是以材料分布的不均勻性來體現(xiàn)整個試件的宏觀非線性。

圖1 細觀單元等效化模型

1.2 擴展有限元（X－FEM）法

Belytschko等［10－12］基于單位分解（PUM）［13－14］概念，利用有限元法的形函數(shù)作為一組單位分解函數(shù)，并將非連續(xù)函數(shù)引入定義在覆蓋（節(jié)點影響域）上的位移模式中，即

式中：I為離散域中所有節(jié)點的集合，J為影響域被非連續(xù)界面穿過的節(jié)點的集合，Ni（x）和Nj（x）分別為節(jié)點i、j影響域上的形函數(shù)，函數(shù) H（x）為Heaviside函數(shù)，ui和uj均為常規(guī)的連續(xù)位移模式的系數(shù)，dj為非連續(xù)位移模式的系數(shù)。引入了非連續(xù)位移模式，使得擴展有限元在對不連續(xù)位移場描述時不再依賴單元邊界，也不需要對網(wǎng)格進行重剖分，可以用來求解函數(shù)不連續(xù)（如裂紋）［15］、導數(shù)不連續(xù)（如夾雜、移動相邊界）［16］、函數(shù)切向分量不連續(xù)（如剪切帶、流固耦合）和概貌加強等問題。文獻［15］證明了擴展有限元在非均勻混凝土材料細觀破壞過程及宏觀力學性能研究中的可行性，但其計算效率低下問題也將限制其在三維力學計算模型中的更多應(yīng)用。

2 細觀單元等效力學性質(zhì)

2.1 細觀單元等效力學性能

采用細觀單元等效化思想來分析混凝土材料細觀斷裂破壞及研究宏觀力學特性時，需要對細觀單元的力學性質(zhì)進行等效化分析。在細觀尺度上，可以認為混凝土是由骨料、砂漿及界面組成的三相復(fù)合材料。采用特征單元尺度剖分混凝土隨機骨料模型后獲得的細觀單元內(nèi)骨料、砂漿及界面相組成的某空間分布，如圖2（a）所示，將其進行力學等效化后形成的均勻各向同性的單元如圖2（b）所示。與圖2（a）宏觀力學特性相等效的圖2（b）的力學性質(zhì)如何？下面，理論推導分析細觀單元的等效體力學性質(zhì)。

圖2（a）中混凝土細觀各相材料（骨料、砂漿及界面）采用如圖3（a）、（b）、（c）所示的本構(gòu)關(guān)系，即：細觀組分在達到峰值應(yīng)力前，應(yīng)力 應(yīng)變關(guān)系是線彈性的；達到強度后，采用線性應(yīng)力 裂縫寬度關(guān)系來表征混凝土材料的軟化行為。ft為極限抗拉強度；Gf為斷裂能；W為張開位移；WS為極限張開位移。

圖2 混凝土細觀單元及其等效體

圖3 混凝土細觀組分本構(gòu)關(guān)系

假定等效體圖2（b）的本構(gòu)行為與各組分本構(gòu)關(guān)系形式類似，如圖3（d）所示。那么在對混凝土細觀單元力學特性進行等效化時，需要獲得的力學參數(shù)有等效彈性模量Eeq，等效泊松比υeq，等效斷裂強度feqt及單元等效斷裂能Geqf。對這些參數(shù)的等效化分析，采用如下處理。

1）等效體彈性模量Eeq及υeq

采用Vogit并聯(lián)模型［17］對彈模及泊松比進行等效化處理，得：

式中：上標“eq”表示細觀等效體，上標“a”、“m”及“i”分別表示骨料、砂漿及界面，C為體積分數(shù)。

2）等效斷裂強度feqt

當細觀等效體達到其斷裂強度feqt時，認為在等效體內(nèi)產(chǎn)生的內(nèi)能（應(yīng)變能）是其細觀組成（骨料、砂漿及界面）的總和，等效體內(nèi)儲存的應(yīng)變能Ωeq為：

即：

對于細觀單元等效體，知：

式中Ωa、Ωm及Ωi分別為骨料、砂漿及界面分別達到其抗拉強度時所存儲的應(yīng)變能。

將式（6）代入式（5），得：

3）斷裂能Geqf

這里假定細觀等效體的斷裂能Geqf與其細觀組分的斷裂能關(guān)系如下：

2.2 細觀單元等效力學特性驗證

為驗證2.1節(jié)所提出的等效化方法的準確性，對細觀單元的等效力學特性進行驗證。

選擇5組細觀單元網(wǎng)格中骨料及砂漿的隨機分布情況，如圖4（a）～（e）所示，將其作為不同的隨機樣本來進行數(shù)值驗證，5組樣本中骨料所占據(jù)的體積分數(shù)均取為30%。采用擴展有限元方法對其進行數(shù)值分析（以單軸拉伸行為為例），5組樣本得到了5組不同的宏觀應(yīng)力－應(yīng)變關(guān)系曲線，如圖5所示。獲得的宏觀應(yīng)力－應(yīng)變關(guān)系曲線與筆者力學等效化理論方法（圖5中“粗線”所示）吻合非常好，從而驗證了筆者等效化力學方法的準確性。筆者進行了大量的隨機樣本數(shù)值試驗，均發(fā)現(xiàn)得到的結(jié)果與2.1節(jié)提出的等效力學特性相吻合，由于篇幅所限，這里不做多述。

3 算例分析

3.1 二維混凝土單軸拉伸破壞過程

圖4 典型樣本及等效體

圖5 等效化方法與數(shù)值結(jié)果對比

圖6（a）為混凝土試件的細觀隨機骨料模型，試件尺寸為150mm×150mm，其細觀組分采用圖3所示的本構(gòu)模型，其力學參數(shù)如表1所示。圖6（b）為采用混凝土材料特征單元尺度對圖6（a）進行網(wǎng)格劃分后，采用前面所述方法對細觀單元進行等效化而獲得的細觀單元等效化模型，不同的單元擁有不用的顏色，表示不同的力學特性（如彈模、泊松比、強度及斷裂能等）。混凝土試件底面豎向固定約束，側(cè)面為自由面，頂面為載荷輸入面，這里采用位移加載控制，如圖6（b）所示。將最大主應(yīng)力準則作為混凝土細觀單元的斷裂準則。

基于細觀單元等效化模型思想理念，采用擴展有限元數(shù)值方法，對單軸拉伸行為下混凝土隨機骨料模型的兩組樣本試件的斷裂破壞過程進行數(shù)值研究，得到2組樣本及其對應(yīng)細觀單元等效化分析方法的裂紋擴展路徑，并獲得其對應(yīng)的混凝土宏觀應(yīng)力 應(yīng)變關(guān)系曲線，如圖7所示。

表1 混凝土各組分力學性質(zhì)

圖6 混凝土細觀有限元模型

圖7 混凝土試件2組樣本裂紋路徑及宏觀p ε關(guān)系

從圖7中可以看出：1）采用筆者提出的細觀單元等效化模型，能很好地模擬混凝土試件的裂紋擴展過程及裂紋模式，與采用細觀隨機骨料模型所得到的最終裂紋路徑基本相同；2）相對于細觀隨機骨料模型（2個樣本的單元數(shù)分別為3 826和3 625）而言，采用細觀單元等效化模型，在使用很少的網(wǎng)格單元下（該算例模型為225個單元）便可以較為準確的獲得混凝土材料的宏觀力學性能，獲得的混凝土材料的宏觀彈模及強度等均與隨機骨料模型結(jié)果吻合良好，僅下降段曲線略有區(qū)別。

3.2 三維混凝土單軸拉伸破壞過程

細觀單元等效化模型的高效性使得三維混凝土細觀裂紋擴展過程及宏觀力學性能的研究成為可能。采用Monte Carlo法建立如圖8（a）所示的二級配混凝土三維隨機骨料模型，并得到與之力學性能相對應(yīng)的三維細觀單元等效化模型，圖8（b）所示，模型尺寸為150mm×150mm×150mm。細觀單元等效化模型中，各細觀單元的尺寸為10mm×10mm×10mm，共3 375個單元，各細觀單元采用圖3（d）所示的本構(gòu)關(guān)系來表征其力學行為。

對混凝土2組樣本試件進行分析研究，獲得其宏觀應(yīng)力 應(yīng)變關(guān)系曲線如圖9所示，分別得到混凝土單軸抗拉強度為3.45、3.38MPa。

圖8 細觀混凝土有限元計算模型

圖9 混凝土材料單軸拉伸宏觀p ε關(guān)系曲線

圖10 為混凝土一組隨機樣本試件在單軸拉伸條件下的最大主應(yīng)力變化云圖及試件的裂紋擴展演化過程。材料的不均勻分布導致混凝土內(nèi)部應(yīng)力分布的不均勻，最終導致混凝土材料的非線性。裂紋擴展狀態(tài)A、B、C、D、E、F下，混凝土宏觀力學性能對應(yīng)于圖9中應(yīng)力 應(yīng)變關(guān)系曲線的幾個狀態(tài)點。

從圖10中可以看出，在外荷載作用下，混凝土試件首先在內(nèi)部薄弱的區(qū)域產(chǎn)生應(yīng)力集中，進而達到細觀單元的抗拉強度而產(chǎn)生斷裂；隨著外載增大，裂紋面不斷擴展演化，較大的區(qū)域產(chǎn)生宏觀裂紋面，當其達到C狀態(tài)時，混凝土宏觀平均應(yīng)力達到最大值，繼而試件整體剛度突降，裂紋面不斷擴展演化，經(jīng)歷D、E狀態(tài)而直至最終的完全斷裂，即F狀態(tài)，混凝土試件應(yīng)力全部釋放?；炷猎嚰?nèi)應(yīng)變能、內(nèi)能及斷裂耗散能與宏觀應(yīng)變的關(guān)系曲線如圖11所示?；炷翍?yīng)變能隨外載荷的增大呈現(xiàn)先增大后減小直至為零的趨勢；體系的斷裂耗散能在試件達到一定的宏觀應(yīng)變時突然產(chǎn)生，并逐漸增大直至達到一恒定值。

圖10 裂紋帶擴展演化過程

圖11 混凝土細觀結(jié)構(gòu)體系能量與應(yīng)變ε的關(guān)系

圖12 為2組不同隨機樣本（骨料空間分布形式不同）的細觀單元等效化模型，在單軸拉伸條件下的最終裂紋模式圖。2個裂紋帶的不同，說明了骨料空間分布的不同，影響混凝土材料的裂紋擴展路徑，這從圖9的宏觀應(yīng)力 應(yīng)變關(guān)系曲線下降段的區(qū)別也可以得到體現(xiàn)。

圖12 兩不同樣本最終裂紋帶

這些研究結(jié)果表明，采用筆者提出的細觀單元等效化方法來研究混凝土細觀斷裂破壞過程，在較低計算量的情況下便可以獲得較為準確的計算結(jié)果，這解決了其它細觀力學模型（如隨機骨料模型，隨機力學特性模型等）計算效率低下的“瓶頸”問題，使得從細觀層次上研究更為復(fù)雜的大體積混凝土材料及混凝土構(gòu)件的破壞機理及力學性能成為可能。

4 結(jié)語

從非連續(xù)介質(zhì)力學角度出發(fā)，在細觀單元等效化力學模型的基礎(chǔ)上，采用擴展有限元方法來研究混凝土材料的細觀裂紋擴展過程。推導了細觀單元的等效力學特性，并進行了數(shù)值驗證；最后對單軸拉伸條件下濕篩混凝土的二維及三維試件進行了數(shù)值分析，研究其裂紋擴展過程及宏觀力學性能，獲得與細觀隨機骨料模型相吻合的結(jié)果，驗證了細觀單元等效化力學模型的準確性和高效性。實現(xiàn)了細觀單元等效化力學方法在非連續(xù)介質(zhì)力學領(lǐng)域的應(yīng)用，從而更清晰地模擬了混凝土在外荷載作用下的斷裂破壞過程。

［1］杜修力，金瀏.混凝土靜態(tài)力學性能的細觀力學方法述評［J］.力學進展，2011，41（4）：411－426.DU Xiuli，JIN Liu.A review on meso－mechanical method for studying the static－mechanical properties of concrete［J］.Advances in Mechanics，2011，41（4）：411－426.

［2］Wang Z M，Kwan A K H，Chan H C.Mesoscopic study of concrete I：generation of random aggregate structure and finite element mesh［J］.Computers and Structures，1999，70（5）：533－544.

［3］Kwan A K H，Wang Z M，Chan H C.Mesoscopic study of concrete II：nonlinear finite element analysis［J］.Computers and Structures，1999，70（5）：545－556.

［4］Schlangen E，Garboczi E J.Fracture simulations of concrete using lattice model computational aspects［J］.Engineering Fracture Mechanics，1997，57（2／3）：319－332.

［5］Mohamed A R，Hansen W.Micromechanical modeling of crack－aggregate interaction in concrete materials［J］.Cement and Concrete Composites，1999，21（5／6）：349－359.

［6］唐春安，朱萬成.混凝土損傷與斷裂 數(shù)值試驗［M］.北京：科學出版社，2003.

［7］Lopez C M，Carlo I，Aguado A.Meso－structural study of concrete fracture using interface elements I：numerical model and tensile behavior［J］.Materials and Structures，2008，41（3）：583－599.

［8］杜修力，金瀏.混凝土材料宏觀力學特性分析的細觀單元等效化模型［J］.計算力學學報，2012，29（5）：654－661.DU Xiuli，JIN Liu.Meso－element equivalent model for macro－scopic mechanical properties analysis of concrete materials ［J］.Chinese Journal of Computational Mechanics，2012，29（5）：654－661.

［9］杜修力，金瀏.混凝土材料細觀單元彈模非均勻統(tǒng)計特性研究［J］.工程力學，2012，29（10）：106－115.DU Xiuli，JIN Liu.Research on the heterogeneous clouds［J］.Computer Methods in Applied Mechanics and Engineering，1996，139：237－262.

［15］黃景琦，金瀏，杜修力.界面特性及骨料分布對混凝土破壞模式的影響［J］.土木建筑與環(huán)境工程，2011，33（S2）：38－41.HUANG Jingqi，JIN Liu，DU Xiuli.The impact of interface characteristic and aggregate distribution on failure modes of concrete ［J］.Journal of Civil，Architecture and Environmental Engineering，2011，33（S2）：38－41.

［16］應(yīng)宗權(quán)，杜成斌，王友元.顆粒增強復(fù)合材料的擴展有限元模擬方法［J］.水利學報，2011，42（2）：198－203.YING Zongquan，DU Chengbin，WANG Youyuan.Numerical simulation of particle reinforced composite using extended finite element method［J］.Journal of Hydraulic Engineering，2011，42（2）：198－203.

［17］Voigt W.Uber die Beziehung zwischen den beiden Elastizit?tskonstanten isotroper K?rper［J］.Wied Ann，1889，274（12）：573－587.statistical properties of elastic modulus of concrete meso－scale element［J］.Engineering Mechanics，2012，29（10）：106－115.

［10］Belytschko T，Black T.Elastic crack growth in finite element with minimal remeshing ［J］.International Journal for Numerical Methods in Engineering，1999，45（5）：610－620.

［11］Moes N，Dolbow J，Belytschko T.A finite element method for crack growth without remeshing ［J］.International Journal for Numerical Methods in Engineering，1999，46（1）：131－150.

［12］Daux C，Dolbow J，Sukumar N，et al.Arbitrary branched and intersecting cracks with the extended finite element method ［J］.International Journal for Numerical Methods in Engineering，2000，48（12）：1741－1760.

［13］Melenk J M，Babuka I.The partition of unity finite element method：basic theory and applications ［J］.Computer Methods in Applied Mechanics and Engineering，1996，139：289－314.

［14］Duarte C A，Oden J T.An H－P adaptive method using