基于三維固定裂縫模型的混凝土開裂數(shù)值模擬

王碗琴，強(qiáng) 晟，郭曉娜，周蘭庭

(1.河海大學(xué)水利水電學(xué)院，南京 210098；2.南水北調(diào)中線干線工程建設(shè)管理局，北京 100038)

作為一種建筑材料[1-3]，混凝土因其取材廣泛、價(jià)格低廉、抗壓強(qiáng)度高等特性成為使用量最大、運(yùn)用最廣的材料。但是混凝土材料本身的弊端也逐漸暴露出來，其組織結(jié)構(gòu)是一個(gè)多尺度多水平體系，內(nèi)部含有大量形狀、大小、位置都不同的微裂縫，在外荷載的作用下，微裂縫易擴(kuò)展至宏觀裂縫導(dǎo)致結(jié)構(gòu)的破壞失穩(wěn)。自1960年Kaplan[4]首次將斷裂力學(xué)應(yīng)用于混凝土開裂問題以來，混凝土斷裂力學(xué)取得了很大進(jìn)展。目前混凝土斷裂模型主要可分為2大類[5]：①以雙參數(shù)模型和雙K模型為典型的等效彈性裂縫模型，此類模型忽略裂縫面上分布的黏聚力，將其等效為線彈性材料，用線彈性力學(xué)來判定裂縫發(fā)展，用P～CMOD,P～CTOD,P～δ曲線來表征混凝土的非線性特性；②以虛擬裂縫模型和裂縫帶模型為典型的黏聚裂縫模型，此類模型便于數(shù)值分析，摒棄了傳統(tǒng)力學(xué)認(rèn)為裂縫尖端應(yīng)力是無窮大的說法，考慮了裂縫面分布的黏聚力，主要是研究裂縫的形態(tài)和裂縫面非線性區(qū)域的軟化特性。

虛擬裂縫模型由Hillerborg提出，考慮微裂面上的拉應(yīng)力σ與裂縫張開度 之間的軟化關(guān)系，但是隨著裂縫面的擴(kuò)展，需要重新劃分網(wǎng)格。而Bazant和Oh[6,7]提出的裂縫帶模型將斷裂能彌散于裂縫帶中，用應(yīng)力～應(yīng)變軟化曲線σ=σ(ε)來描述裂縫面的軟化關(guān)系，通過改變剛度矩陣來模擬開裂，不需要重新劃分網(wǎng)格,通過調(diào)整應(yīng)力～應(yīng)變曲線來適應(yīng)不同網(wǎng)格大小，以保證斷裂能的唯一。但是斷裂帶模型的缺點(diǎn)[8]在于模型是用單元內(nèi)平行密集的裂縫帶來模擬裂縫，無法確定裂縫的寬度和位置。

根據(jù)對(duì)裂縫方向處理方法的不同，裂縫帶模型又可以分為固定裂縫模型，多向固定裂縫模型和旋轉(zhuǎn)裂縫模型。固定裂縫模型是假定裂縫方向固定，不隨著裂縫擴(kuò)展而變化，允許同一點(diǎn)產(chǎn)生多條裂縫，但同時(shí)也存在一定的缺陷[9,10]：一方面它沒有考慮到主應(yīng)力軸的旋轉(zhuǎn)問題，會(huì)造成結(jié)果的突變；另一方面，沒有考慮到與裂縫平行的混凝土?xí)霈F(xiàn)一定的彈性回彈，而認(rèn)為裂縫附近的單元均超過抗拉強(qiáng)度，造成開裂單元數(shù)的誤判。旋轉(zhuǎn)裂縫模型的開裂方向則隨著主應(yīng)力軸或者主應(yīng)變軸方向的變化而變化，但是在任意一點(diǎn)只允許產(chǎn)生一條裂縫。多向固定裂縫模型[11]則允許在一點(diǎn)形成任意條裂縫，且將裂縫角度限制在一定范圍內(nèi)，應(yīng)用方面沒有固定裂縫模型和旋轉(zhuǎn)裂縫模型廣泛。本文基于固定裂縫模型理論，編制Fortran程序來模擬混凝土開裂。

1 三維固定裂縫模型理論

1.1 本構(gòu)關(guān)系

固定裂縫模型基于應(yīng)變分解[12]思想，將材料應(yīng)變分解成彈性和裂縫應(yīng)變，即：

ε=εe+εcr

(1)

在開裂面建立局部坐標(biāo)系n,s,t，局部坐標(biāo)系與整體坐標(biāo)系之間應(yīng)力關(guān)系如下：

εcr=Necr

(2)

scr=NTσ

(3)

式中：σ、εcr分別為整體坐標(biāo)系下單元應(yīng)力和裂縫應(yīng)變；scr、ecr分別為局部坐標(biāo)系下裂縫應(yīng)力和裂縫應(yīng)變；N為坐標(biāo)變換矩陣：

(4)

忽略裂縫面法向和切向應(yīng)力與應(yīng)變之間的相互影響[13]，定義局部坐標(biāo)系下的應(yīng)力和應(yīng)變關(guān)系如下：

scr=Dcrsecantecr

(5)

(6)

式中：DⅠsecant、DⅡsecant分別表示Ⅰ型、Ⅱ型裂縫割線模量。

混凝土彈性本構(gòu)關(guān)系為:

σ=Dεe

(7)

式中：D為彈性矩陣。

則綜合式(1)、式(2)、式(3)、式(5)、式(7)可得三維混凝土開裂本構(gòu)關(guān)系[14]：

σ=[D-DN(Dcrsecant+NTDN)-1NTD]ε

(8)

考慮到混凝土并不是完全的脆性斷裂，而是存在一個(gè)軟化過程，引入折減系數(shù)μ來體現(xiàn)混凝土從微裂到完全開裂這一過程的軟化特性，引入剪切保留系數(shù)β來表示裂縫面的抗剪能力，一般取0.2。將混凝土和裂縫面的剛度視為串聯(lián)關(guān)系[15],如圖1所示。

圖1 混凝土與裂縫面剛度的串聯(lián)關(guān)系Fig.1 The series relationship of concrete and fracture stiffness

可得：

(9)

(10)

假定裂縫局部坐標(biāo)系與整體坐標(biāo)系重合，則n,s,t軸方向矢量分別為：

(11)

(12)

(13)

代入式(4)可得：

(14)

將式(9)、式(10)、式(14)代入式(8)計(jì)算可得裂縫局部坐標(biāo)系與整體坐標(biāo)系重合時(shí)混凝土開裂本構(gòu)關(guān)系：

1.2 應(yīng)力～應(yīng)變曲線

應(yīng)力～應(yīng)變曲線是固定裂縫模型的一個(gè)重要參數(shù)，有線性、雙線性、指數(shù)型等，其中線性應(yīng)用最多，本文選取線性軟化關(guān)系，如圖2所示。

圖2 應(yīng)力～應(yīng)變曲線Fig.2 The stress and strain curve

經(jīng)推導(dǎo)可得折減系數(shù)的公式[16]如下：

(其中Et為負(fù)值)

(16)

式中：εult為極限拉應(yīng)變；ε為整體坐標(biāo)系下的全應(yīng)變；如果ε>εult，則μ=0；如果ε<εn，則μ=1；否則μ則采用公式(16)計(jì)算。

1.3 斷裂能Gf

斷裂能是試件發(fā)生斷裂外荷載所做的功，即試件斷裂所釋放的應(yīng)變能。斷裂能越大，則材料的延性越高，反之材料越脆。固定裂縫模型中，斷裂能和應(yīng)力～應(yīng)變曲線有如下關(guān)系[15]：

Gf=Wareah

(17)

式中：Warea為應(yīng)力～應(yīng)變曲線所包圍的面積；h為裂縫帶寬度，即裂縫方向網(wǎng)格大小。

則線性軟化關(guān)系中有：

(18)

傳統(tǒng)彌散裂縫模型[17-20]，應(yīng)變軟化的結(jié)果對(duì)網(wǎng)格大小具有敏感性，如果帶有應(yīng)變軟化特性的材料僅僅按照固定的應(yīng)力～應(yīng)變關(guān)系定義，或者說由固定的應(yīng)力～應(yīng)變關(guān)系確定，那么分散在斷裂區(qū)的斷裂能大小將隨著單元加密而減小。而固定裂縫模型中，通過調(diào)整應(yīng)力～應(yīng)變曲線的下降段來保證斷裂能的唯一。

2 算例驗(yàn)證

本文基于三維固定裂縫模型，編譯Fortran程序，對(duì)3點(diǎn)彎曲梁斷裂進(jìn)行數(shù)值模擬，程序總體思路見圖3。試件長2 m,寬0.2 m，高0.45 m。采用分步加載的形式，200個(gè)荷載步，每個(gè)荷載步大小為0.5 kN。混凝土彈性模量E=28.00 GPa,泊松比υ=0.2,剪切保留系數(shù)β=0.2，抗拉強(qiáng)度ft=2.5 MPa。選用3種不同網(wǎng)格尺寸，裂縫方向網(wǎng)格尺寸h分別為0.10，0.04，0.02 m。

圖3 程序計(jì)算流程圖Fig.3 The flow chart of calculation

2.1 算例1

本算例驗(yàn)證斷裂能的守恒準(zhǔn)則。選用Gf=187.5 N/m，由式(18)可計(jì)算不同網(wǎng)格大小所對(duì)應(yīng)的極限抗拉應(yīng)變，3種網(wǎng)格參數(shù)具體見表1。

表1 不同網(wǎng)格參數(shù)Tab.1 The parameters of different grid models

選取特征點(diǎn)(1,0.1,0)進(jìn)行計(jì)算，繪制特征點(diǎn)應(yīng)力～應(yīng)變曲線如圖4～圖6所示。

圖4 h=0.10 m的特征點(diǎn)應(yīng)力～應(yīng)變曲線Fig.4 h=0.10 m the stress and strain curve of keypoint

圖5 h=0.04 m的特征點(diǎn)應(yīng)力～應(yīng)變曲線Fig.5 h=0.04 m the stress and strain curve of keypoint

圖6 h=0.02 m的特征點(diǎn)應(yīng)力～應(yīng)變曲線Fig.6 h=0.02 m the stress and strain curve of keypoint

從圖4～圖6可以看出：計(jì)算所得應(yīng)力～應(yīng)變曲線所圍成的面積和假定應(yīng)力～應(yīng)變曲線所圍成的面積誤差很小，即計(jì)算所得斷裂能基本符合預(yù)先設(shè)定的斷裂能參數(shù)。所以固定裂縫模型中，根據(jù)網(wǎng)格大小來調(diào)整應(yīng)力～應(yīng)變曲線，可以保證斷裂能的一致性，克服網(wǎng)格敏感性問題。由于該模型將斷裂區(qū)之外的混凝土視為彈性材料，未考慮其非線性變形,且計(jì)算求解存在一定誤差，所以計(jì)算所得下降段與假定的應(yīng)力～應(yīng)變曲線并不是完全吻合。

2.2 算例2

本算例選用不同斷裂能對(duì)同一個(gè)網(wǎng)格尺寸模型進(jìn)行計(jì)算，具體參數(shù)及結(jié)果見表2，其中極限荷載是指特征點(diǎn)達(dá)到極限拉應(yīng)變?chǔ)舥lt時(shí)所對(duì)應(yīng)的荷載大小。

表2 計(jì)算參數(shù)Tab.2 The parameters of calculation

不同斷裂能工況下特征點(diǎn)(1,0.1,0)的荷載-位移曲線和應(yīng)變分布云圖見圖7～圖11。

圖7 不同斷裂能下特征點(diǎn)的荷載～位移曲線Fig.7 The load displacement curves of keypoint with different fracture parameters

圖8 不考慮開裂的應(yīng)變包絡(luò)圖Fig.8 The strain envelope diagram regardless of crack

圖9 Gf=130.0 N/m應(yīng)變包絡(luò)圖Fig.9 Gf=130.0 N/m strain envelope diagram

圖10 Gf=150.0 N/m應(yīng)變包絡(luò)圖Fig.10 Gf=150.0 N/m strain envelope diagram

圖11 Gf=187.5 N/m應(yīng)變包絡(luò)圖Fig.11 Gf=187.5 N/m strain envelope diagram

由荷載～位移曲線圖7和應(yīng)變分布圖8～圖11可以看出，在相同荷載條件下，斷裂能越小，節(jié)點(diǎn)位移、應(yīng)變值越大，斷裂能越大，節(jié)點(diǎn)位移、應(yīng)變值越小。因此斷裂能越小，試件斷裂時(shí)外荷載所做的功越小，即材料越脆，在相同荷載條件下，越容易產(chǎn)生裂縫，而斷裂能越大，即材料延性越好，試件承載力越大，符合實(shí)際情況。

3 結(jié) 語

本文基于三維固定裂縫理論，編制了相應(yīng)的計(jì)算程序，通過2個(gè)算例驗(yàn)證了程序的正確性，并得到如下結(jié)論。

(1)固定裂縫模型根據(jù)網(wǎng)格大小調(diào)整應(yīng)力～應(yīng)變曲線，可以使計(jì)算結(jié)果不會(huì)因?yàn)榫W(wǎng)格的疏密而不同，克服網(wǎng)格敏感性問題。

(2)斷裂能作為一種材料特性參數(shù)對(duì)混凝土結(jié)構(gòu)性能的影響較大，斷裂能越小，即混凝土材料脆性越大，可承載的荷載值就越小。

(3)混凝土彈性理論應(yīng)力隨著荷載線性增加，而固定裂縫模型作為一種非線性斷裂力學(xué)模型，能夠考慮斷裂區(qū)的三維應(yīng)力狀態(tài)，考慮混凝土材料所具有的軟化特性，即應(yīng)力達(dá)到抗拉強(qiáng)度之后是會(huì)隨著應(yīng)變?cè)黾佣饾u減小，直至完全開裂，應(yīng)力降為0，比較符合實(shí)際情況。

□

[1] Taliercio A L F, Gobbi E.Fatigue life and change in mechanical properties of plain concrete under triaxial deviatoric cyclic stress[J].Magazine of Concrete Research,1998,50(3):247-255.

[2] Yin W, Hsu Thomas T C.Fatigue behavior of steel fiber reinforced concrete in uniaxial and biaxial compression[J].ACI Materials Journal,

1995,92(1):71-81.

[3] 劉智光.混凝土破壞過程細(xì)觀數(shù)值模擬與動(dòng)態(tài)力學(xué)特性機(jī)理研究[D]. 大連：大連理工大學(xué),2012.

[4] 蔣梅玲.基于斷裂力學(xué)和損傷理論的混凝土斷裂模型試驗(yàn)和數(shù)值研究[D]. 杭州：浙江大學(xué),2010.

[5] 郭少華.混凝土破壞理論研究進(jìn)展[J].力學(xué)進(jìn)展,1993,(4):520-529.

[6] Bazant Z P. Crack band theory for fracture of concrete[J].Materials and Structures,1983,16(93):155-177.

[7] 李 翠.混凝土直接剪切試驗(yàn)的非線性有限元分析[D].大連：大連理工大學(xué),2010.

[8] 張楚漢.論巖石、混凝土離散-接觸-斷裂分析[J].巖石力學(xué)與工程學(xué)報(bào),2008,(2):217-235.

[9] 姜慶遠(yuǎn),葉燕春,劉宗仁.彌散裂縫模型的應(yīng)用探討[J].土木工程學(xué)報(bào),2008,(2):81-85.

[10] 韓 濤,安雪暉.鋼筋混凝土三維多向固定裂縫本構(gòu)模型[J].清華大學(xué)學(xué)報(bào)(自然科學(xué)版),2008,(6):947-950.

[11] 曲 蒙.早齡期混凝土損傷斷裂的數(shù)值分析[D]. 杭州：浙江大學(xué),2006.

[12] De Borst R,Nauta P. Non-orthogonal cracks in smeared finite element model[J].Engineering Computation,1985,2(3):35-46.

[13] 唐國斌,王統(tǒng)寧,項(xiàng)貽強(qiáng).基于應(yīng)變分解方法的鋼筋混凝土開裂模型[J].工程力學(xué),2013,30(3):174-180.

[14] 周元德,張楚漢,金 峰.混凝土開裂的三維非線性數(shù)值模擬[J].清華大學(xué)學(xué)報(bào)(自然科學(xué)版),2003,43(11):1 542-1 545.

[15] Rots J G, de Borst R. Analysis of mixed-mode fracture in concrete [J]. J of Eng Mech,1987,113(11):1 739-1 758.

[16] 徐世烺.混凝土斷裂力學(xué)[M].北京:科學(xué)出版社,2011.

[17] Cervenka Jan，Bazant Zdenek P. Equivalent localization element for crack band approach to mesh-sensitivity in microplane model[J].International Journal for Numerical Methods in Engineering,2005,62(5):700-726.

[18] Rafi Muhammad Masood，Nadjai Ali.Comparison of numerical behaviors of FRP reinforced concrete beams using three smeared crack models[J].Materials and Structures/Materiaux et Constructions,2012,45(1-2):93-106.

[19] Mostofinejad D,Talaeitaba S B. Nonlinear modeling of RC beams subjected to torsion using the smeared crack model[J].Procedia Engineering,2011,14:1 447-1 454.

[20] Ozbolt Josko, Bazant Zdenek.Numerical smeared fracture analysis: nonlocal microcrack interaction approach[J].International Journal for Numerical Methods in Engineering,1996,39(4):635-661.