袁衛(wèi)國 李亞春 李 賁 曾祥茂 朱克東 羅祝君

（廣州市第一市政工程有限公司，廣州510060）

0 引 言

ABAQUS CDP（Concrete Damage Plasticity）模型被廣泛地應(yīng)用于混凝土及鋼-混凝土組合結(jié)構(gòu)的非線性有限元數(shù)值模擬分析中［1-3］，該模型將混凝土宏觀特征的非線性本構(gòu)關(guān)系作為有限元單元的性質(zhì)（細(xì)觀尺度特征）而參與迭代分析計(jì)算，這種處理方法沒有考慮組成模型的單元的個(gè)體性差異（非均質(zhì)性），所以有限元數(shù)值模擬計(jì)算結(jié)果不存在宏細(xì)觀尺度上的差異［4-5］，無法得到混凝土材料非均質(zhì)性對(duì)某一變量的定量影響結(jié)果?；炷敛牧贤嬖诟鞣N自身缺陷，例如孔洞和微裂紋等，而且各組分材料的力學(xué)性質(zhì)參數(shù)具有一定的隨機(jī)性和離散性，而不是通常計(jì)算時(shí)所采用的定值［6-7］。大多數(shù)文獻(xiàn)研究將混凝土視為三相復(fù)合材料，對(duì)每一相材料賦予單一的材料屬性，且主要利用二維隨機(jī)骨料模型模擬混凝土構(gòu)件的斷裂破壞［8-12］，很多學(xué)者對(duì)三維隨機(jī)骨料模型的建立提出了新的不同的方法［13-15］。到目前為止，尚未有學(xué)者對(duì)基于塑性損傷和三維隨機(jī)骨料模型的混凝土彈性模量概率密度分布進(jìn)行研究?；诖?，本文將混凝土視為一相非均質(zhì)材料，通過ABAQUS 自帶的CDP 模型模擬混凝土的塑性損傷開裂破壞，通過三維隨機(jī)骨料模型模擬混凝土的非均質(zhì)性，對(duì)兩種概率密度分布下的彈性模量進(jìn)行分析研究，通過與文獻(xiàn)試驗(yàn)對(duì)比，給出了兩種分布下描述混凝土非均質(zhì)性的參數(shù)取值范圍，為研究混凝土的非均質(zhì)性提供了一定的參考。

1 理 論

1.1 正態(tài)分布隨機(jī)變量X的概率密度

正態(tài)分布隨機(jī)變量X的概率密度為

式中，均值μ及標(biāo)準(zhǔn)差σ(σ ＞0)都是常數(shù)。

均值μ=100 時(shí)不同標(biāo)準(zhǔn)差對(duì)應(yīng)的概率密度函數(shù)如圖1所示。

1.2 Weibull分布隨機(jī)變量X的概率密度

Weibull分布隨機(jī)變量X的概率密度為

式中，m 和μ0都是常數(shù)，m 表示W(wǎng)eibull 概率密度分布的均質(zhì)度，反映隨機(jī)變量x 取值的離散程度；μ0是一個(gè)與自變量平均值有關(guān)的參數(shù)，表示W(wǎng)eibull概率密度分布的數(shù)值尺度，其數(shù)值并不等于自變量平均值。μ0=100 時(shí)不同均質(zhì)度密度函數(shù)如圖2所示。

圖1 μ=100時(shí)不同標(biāo)準(zhǔn)差對(duì)應(yīng)的概率密度函數(shù)Fig.1 Probability density function of different standard deviations with μ=100

圖2 μ0=100時(shí)不同均質(zhì)度概率密度函數(shù)Fig.2 Probability density function of different homogeneities with μ0=100

1.3 混凝土應(yīng)力應(yīng)變關(guān)系

ABAQUS中的損傷模型是一個(gè)基于塑性力學(xué)的連續(xù)介質(zhì)損傷模型，其中混凝土應(yīng)力和應(yīng)變關(guān)系如下所示：

1.4 塑性流動(dòng)法則

ABAQUS 塑性勢(shì)函數(shù)采用了雙曲線數(shù)學(xué)模型，該模型基于非關(guān)聯(lián)的塑性流動(dòng)準(zhǔn)則，其數(shù)學(xué)關(guān)系式如下所示：

式中，參數(shù)φ(θ，fi)是一個(gè)表示偏心的參數(shù)，ABAQUS中其默認(rèn)值為0.1；σt0(θ，fi)=σt|-εplt=0為計(jì)算之前設(shè)定的表示初始損傷時(shí)的單軸拉伸應(yīng)力值；ψ(θ，fi)為在三維應(yīng)力空間中p-q所確定的平面上在高側(cè)壓力作用下的以角度表征的擴(kuò)容角=表示第四強(qiáng)度理論的Mises 等效有效應(yīng)力，-S= -p·I+ -σ 表示應(yīng)力空間中有效應(yīng)力張量的偏分量；=表示三維應(yīng)力空間有效靜水壓力，其中I為應(yīng)力不變量。

1.5 三維塑性屈服空間函數(shù)

屈服函數(shù)F(-σ，ε?pl)表征一個(gè)三維空間函數(shù)，它可以計(jì)入材料在拉伸和壓縮狀態(tài)下材料強(qiáng)度值的不同變化規(guī)律，有效應(yīng)力三維空間屈服函數(shù)的屈服面演化規(guī)律由等效塑性應(yīng)變-εtpl和-εcpl來表征控制，其關(guān)系式如下式所示：

式中，α 由Kupfer 曲線折算而求出，其值等于雙軸抗壓與單軸抗壓極限強(qiáng)度的比值第四強(qiáng)度理論Mises 等效應(yīng)力，其中為的應(yīng)力偏張量=-：I為有效應(yīng)力，I 為應(yīng)力不變量；其中和分別為塑性拉應(yīng)變和塑性壓應(yīng)變和分別為有效抗拉強(qiáng)度和有效抗壓強(qiáng)度；符號(hào)表示=(|x+x)，表示-σ代數(shù)特征值的最大值；γ表示第二應(yīng)力不變量在三維空間拉伸子午屈服面和壓縮子午屈服面取值的比值。

2 MATLAB模型和ABAQUS模型

為了保證本文數(shù)值模擬數(shù)據(jù)的可靠性和有效性，本文按照文獻(xiàn)［16］建模，如圖3所示為試驗(yàn)三點(diǎn)彎曲鋼筋混凝土梁，梁尺寸為3 m×0.2 m×0.45 m，梁底端向上偏移40 mm 處放置有3 根直徑為20 mm的鋼筋。按文獻(xiàn)［16］設(shè)置材料屬性，鋼筋彈性模量為2.1×105MPa，泊松比為0.3，受拉屈服強(qiáng)度為500 MPa，不考慮受壓和受拉損傷；混凝土宏觀尺度上的彈性模量為3.35×1010Pa，骨料彈性模量由概率密度函數(shù)隨機(jī)生成，骨料泊松比取0.2，混凝土每一顆骨料都采用塑性損傷本構(gòu)模型，粘滯系數(shù)設(shè)置為0.000 5，斷裂能為76.5 N/m，抗壓強(qiáng)度和抗拉強(qiáng)度分別為為38 MPa 和3.5 MPa。利用MATLAB 中normrnd 和wblrnd 命令隨機(jī)生成服從正太分布和Weibull分布的數(shù)據(jù)作為彈性模量，利用PYTHON 語言將彈性模量數(shù)據(jù)導(dǎo)入ABAQUS中，將每一顆骨料賦予不同的彈性模量，混凝土塑性損傷參數(shù)如圖4 所示；為了減小跨中加載點(diǎn)處的應(yīng)力集中，在ABAQUS 有限元模型加載點(diǎn)處設(shè)置一剛性墊塊，尺寸為0.2 m×0.2 m×0.2 m，彈性模量為2.1×106MPa，泊松比為0.3。鋼筋采用T3D2單元，混凝土隨機(jī)骨料和墊塊采用C3D8R 單元。在墊塊正上方15 mm 處建立參考點(diǎn)，對(duì)墊塊上表面和參考點(diǎn)建立耦合約束，由于塑性損傷本構(gòu)迭代運(yùn)算極不容易收斂，因此本文對(duì)參考點(diǎn)進(jìn)行位移加載。鋼筋與混凝土的接觸設(shè)置為embedded，墊塊與混凝土的接觸設(shè)置為tie，三點(diǎn)彎曲梁兩端支座約束為固定約束，采用自行編制的MATLAB程序建立模型，MATLAB 中三維隨機(jī)骨料模型如圖5所示（骨料最小粒徑10 mm，最大粒徑20 mm，共計(jì)8 000 顆），與MATLAB 模型對(duì)應(yīng)的ABAQUS有限元模型如圖6所示。

圖3 文獻(xiàn)［16］三點(diǎn)彎曲梁（單位：mm）Fig.3 Three-point bending beam in Ref.［16］

利用MATLAB 語言編制三維隨機(jī)骨料投放程序，生成骨料的空間位置和粒徑信息并輸出到TEXT 文檔，利用PYTHON 語言參數(shù)化建模流程如下：

（1）利用PYTHON 語言將骨料信息導(dǎo)入到ABAQUS有限元軟件中；

（2）利用PYTHON 語言中的For 循環(huán)和Range（）遍歷函數(shù)建立每一顆骨料的part模型；

（3）利用PYTHON 語言中的For 循環(huán)和Range（）遍歷函數(shù)對(duì)每一顆骨料賦予材料屬性；

圖4 ABAQUS塑性損傷參數(shù)Fig.4 Plasticity damage parameters in ABAQUS

（4）利用PYTHON語言中的For循環(huán)和Range（）遍歷函數(shù)將每一顆骨料裝配到assembly模型中；

（5）由于骨料與梁構(gòu)件的空間重合，利用PYTHON 語言中的For 循環(huán)和Range（）遍歷函數(shù)依次用每一顆骨料進(jìn)行布爾cut操作，即在有骨料的位置處切割形成孔洞；

（6）利用PYTHON語言中的For循環(huán)和Range（）遍歷函數(shù)依次將每一顆骨料重新裝配到assembly模型相應(yīng)的位置；

（7）利用PYTHON語言中的For循環(huán)和Range（）遍歷函數(shù)將骨料和切割后的梁構(gòu)件進(jìn)行布爾merge操作。

圖5 MATLAB隨機(jī)骨料模型Fig.5 Random aggregate model in MATLAB

圖6 ABAQUS隨機(jī)骨料模型Fig.6 Random aggregate model in ABAQUS

3 數(shù)值模擬與試驗(yàn)對(duì)比分析

如圖7 所示為不同工況下?lián)p傷云圖，括號(hào)內(nèi)的數(shù)據(jù)為本文數(shù)值模擬得到的三維隨機(jī)骨料鋼筋混凝土梁的極限承載力。

本文數(shù)值模擬梁跨中右側(cè)斜裂縫基本呈45°擴(kuò)展，通過圖7 和圖8 對(duì)比可以看出，本文數(shù)值模擬和文獻(xiàn)［16］裂縫擴(kuò)展形態(tài)基本吻合。從圖9 可以看出荷載位移曲線有明顯的突變點(diǎn)，各種工況荷載-撓度曲線在突變點(diǎn)之前基本完全一致，突變前為直線，突變后曲線斜率越來越小，圖9 給出了試驗(yàn)六根梁1/2 結(jié)構(gòu)的荷載-撓度曲線，其平均值為369.4 kN，圖7 所示工況最小極限荷載為335.792 kN，最大極限荷載為337 kN，最小誤差為8.8%，最大誤差為9.1%，說明本文建立的基于非均質(zhì)性和塑性損傷模型的有限元數(shù)值計(jì)算模型是可行和有效的。

雖然正態(tài)分布能表征混凝土的非均質(zhì)性，但是由于其標(biāo)準(zhǔn)差取值往往在107數(shù)量級(jí)左右才能表征混凝土的非均質(zhì)性，因此筆者建議考慮混凝土非均質(zhì)性時(shí)，采用Weibull 概率密度函數(shù)，下面對(duì)Weibull均質(zhì)度取值進(jìn)行分析。

圖7 不同工況損傷云圖Fig.7 Damage contour of different cases

圖8 文獻(xiàn)［16］試驗(yàn)裂縫擴(kuò)展形態(tài)Fig.8 Crack propagation mode of experiment in Ref.［16］

4 Weibull分布均質(zhì)度取值分析

從圖10 和圖11 可以看出，隨著Weibull 分布均質(zhì)度m 取值的不斷增大，混凝土宏觀尺度上的彈性模量值將接近于細(xì)觀尺度上的彈性模量值的平均值。將MATLAB wblrnd 命令生成的混凝土彈性模量值用相應(yīng)Weibull 分布的均值進(jìn)行均一化處理，不斷調(diào)整均質(zhì)度m的取值范圍后發(fā)現(xiàn)，當(dāng)均質(zhì)度在1到200之間時(shí)，用單調(diào)增函數(shù)指數(shù)函數(shù)擬合均一化后彈性模量與均質(zhì)度m 的關(guān)系，可以得到如圖10 所示的擬合關(guān)系式（6）（E0為Weibull分布的均值），但擬合效果較差。當(dāng)均質(zhì)度取值范圍在1～50 之間時(shí)，可以得到如圖11 所示的擬合關(guān)系式（7）（E0為Weibull 分布的均值），此時(shí)擬合效果較好，這兩個(gè)擬合方程分別如下所示。

在進(jìn)行混凝土非均質(zhì)性數(shù)值模擬時(shí)，混凝土彈性模量均值和均質(zhì)度可以根據(jù)上述擬合公式進(jìn)行選取。

圖9 文獻(xiàn)［16］與本文荷載-撓度曲線對(duì)比Fig.9 The comparison of the load-deflection curve

圖10 均質(zhì)度在1～200之間的擬合曲線Fig.10 Fitting curve with homogeneity between 1 and 200

圖11 均質(zhì)度在1～50之間的擬合曲線Fig.11 Fitting curve with homogeneity between 1 and 50

5 結(jié) 論

（1）通過ABAQUS 自帶的CDP 模型模擬混凝土的塑性損傷，通過三維隨機(jī)骨料模型模擬混凝土的非均質(zhì)性可以精確地模擬混凝土的斷裂破壞過程。

（2）傳統(tǒng)的有限元模型一個(gè)單元只有一種材料屬性，本文通過PYTHON 語言進(jìn)行切割合并布爾操作，給不同的骨料賦予不同的材料屬性，劃分網(wǎng)格后一個(gè)單元可能被賦予多種材料屬性，可以更大程度地模擬非均質(zhì)性，比傳統(tǒng)方法有較大的改進(jìn)。本文提出的這種建模方法可以作為模擬非均質(zhì)性三維建模的一種有效的新方法。

（3）由于混凝土骨料數(shù)目巨大且涉及高度非線性運(yùn)算，采用程序語言循環(huán)函數(shù)建模可以提高建模效率。對(duì)于混凝土非均質(zhì)性的模擬，將混凝土視為一相非均質(zhì)材料比視為三相非均質(zhì)材料節(jié)約建模時(shí)間，而且通過合理取值可以得到較高的數(shù)值精度。

（4）混凝土彈性模量非均質(zhì)特征包括兩方面內(nèi)容：第一，不同的Weibull 參數(shù)必定代表不同的非均質(zhì)材料；第二，相同的Weibull分布參數(shù)，由于數(shù)據(jù)是隨機(jī)產(chǎn)生的，材料參數(shù)在數(shù)值上往往不同。

（5）彈性模量非均質(zhì)性用正態(tài)分布表征時(shí)，正態(tài)分布的標(biāo)準(zhǔn)差表征非均質(zhì)度，標(biāo)準(zhǔn)差最合理取值為107～108；彈性模量非均質(zhì)性用Weibull 分布表征時(shí)，均質(zhì)度最合理取值為50，Weibull 分布比正態(tài)分布能更好地表征混凝土彈性模量的非均質(zhì)性，標(biāo)準(zhǔn)差和均質(zhì)度取值不合理會(huì)導(dǎo)致ABAQUS CDP本構(gòu)自動(dòng)檢測(cè)通不過而無法計(jì)算。

（6）混凝土其他材料參數(shù)（泊松比、抗壓強(qiáng)度、抗拉強(qiáng)度）的非均質(zhì)性概率密度分布類型將是下一個(gè)研究的課題。