一種新的混凝土各向異性彈塑性損傷本構(gòu)模型及其數(shù)值實施

焦延濤，程立平

(1. 華北水利水電大學(xué)水利學(xué)院，鄭州 450045；2. 平頂山學(xué)院化學(xué)與環(huán)境工程學(xué)院，平頂山 467000)

目前，國內(nèi)外關(guān)于混凝土塑性-損傷模型的研究成果較多，但對于混凝土這種準(zhǔn)脆性材料，由于模型的復(fù)雜性，對其進行非線性有限元計算分析仍然是相當(dāng)具有挑戰(zhàn)性的任務(wù)。

對于混凝土類材料，其非線性特征主要有彈塑性、粘性、損傷、斷裂等。而損傷對材料的非線性行為起著很大的作用，特別是在受拉情況下，但混凝土在受壓情況下，在宏觀裂縫形成之前，除了破壞外，還表現(xiàn)出明顯的塑性性能。由于損傷和塑性變形均對混凝土的非線性行為均具有貢獻，因此，單純的采用彈性-損傷模型來描述混凝土的非線性行為是不夠的[1]。

目前，現(xiàn)有的彈塑性損傷本構(gòu)模型，大多是以熱力學(xué)為基礎(chǔ)，通過連續(xù)損傷力學(xué)建立相關(guān)的熱力學(xué)勢，按照經(jīng)典塑性力學(xué)的方法，用正交流動法則建立熱力學(xué)基礎(chǔ)的損傷演化方程。文獻[2 - 13]就是在熱力學(xué)框架下通過建立耗散勢，從而建立了混凝土材料的彈塑性損傷模型。另外，近年來由吳建營等[14-17]提出并發(fā)展起來的統(tǒng)一相場損傷理論，也屬于此類模型的范疇。這種方法建立的損傷演化方程符合嚴(yán)密的力學(xué)推理，有堅實的理論基礎(chǔ)，算法的魯棒性也較好，而且模型也能較為完整的描述混凝土材料在復(fù)雜應(yīng)力狀態(tài)下的非線性行為，盡管此類模型的優(yōu)點很多，但損傷耗散勢的建立是一項較為困難的工作，對研究者的力學(xué)功底要求較高，而且損傷耗散勢的試驗驗證也有一定的難度?？紤]到建立損傷耗散勢的難度，一種將基于應(yīng)變的損傷演化方程和基于有效應(yīng)力的塑性方程結(jié)合起來的彈塑性損傷本構(gòu)模型受到了一些學(xué)者的歡迎，文獻[18 - 27]即通過這種方法構(gòu)建了混凝土類材料的彈塑性損傷本構(gòu)模型。通過這種方法建立的彈塑性損傷本構(gòu)模型，不需要通過迭代就可以得到損傷變量的值，這樣就簡化了非線性計算過程，在數(shù)值計算時也能獲得穩(wěn)定的解。但由于考慮各向異性損傷的復(fù)雜性，大多作者均只采用了一個各向同性的損傷變量，這并不能反映混凝土材料在復(fù)雜應(yīng)力狀態(tài)下?lián)p傷各向異性的特性，也不能反映混凝土材料拉壓損傷不等性的特點。

有鑒于此，本文借鑒前人的研究成果，采用兩種不同的損傷演化方程，即由雙曲型函數(shù)表示的拉伸損傷演化方程和由指數(shù)型函數(shù)表示的壓縮損傷演化方程，并由應(yīng)變來驅(qū)動損傷的演化以及控制損傷的各向異性特性。如文獻[19, 28]所述，損傷主要是一種應(yīng)變控制的現(xiàn)象，因此通過應(yīng)變來驅(qū)動損傷的演化是合理可行的。

由于在有效應(yīng)力空間中建立的塑性損傷模型，其塑性部分在數(shù)值計算上更穩(wěn)定(見文獻[2 - 4])，因此本文在有效應(yīng)力空間中建立了該模型的塑性屈服準(zhǔn)則、不同的運動硬化準(zhǔn)則和各向同性流動準(zhǔn)則。此外，由于通過使用應(yīng)變等效假設(shè)，假定有效構(gòu)型中的總應(yīng)變，包括彈性和塑性應(yīng)變，與名義構(gòu)型中的總應(yīng)變相等，因此有效應(yīng)力和損傷演化的計算可以通過解耦的算法來實現(xiàn)。即計算過程中，應(yīng)變的計算無需考慮損傷對其的影響，而有效應(yīng)力的計算在有效構(gòu)型中進行，也無需考慮損傷對其的影響。如此以來，有效應(yīng)力以及應(yīng)變的計算即可按照經(jīng)典塑性理論的方法來實現(xiàn)，從而可簡化本構(gòu)方程的推導(dǎo)，為數(shù)值實現(xiàn)提供優(yōu)勢。另外，采用了解耦的算法后，模型的損傷部分則可進行獨立計算，從而可顯式實現(xiàn)。

另外，如文獻[29]所述，混凝土材料的應(yīng)力-應(yīng)變曲線分為上升段及下降段。在應(yīng)力-應(yīng)變曲線的上升段，剛度矩陣為正定的，進行彈塑性計算分析時，采用常規(guī)的迭代方法可以獲得滿意的解答，但在應(yīng)力-應(yīng)變曲線的下降段，若仍然采用常規(guī)的迭代方法，則計算不能收斂，這是因為此時剛度矩陣不再是正定的，可稱之為“負(fù)剛度”。針對此問題，國內(nèi)外學(xué)者提出了許多相關(guān)的方法，以克服這種迭代不收斂情況。目前，常用的這些方法有逐步搜索法、虛假剛性彈簧法、位移控制法、強制迭代法、硬化剛度法、弧長法等。盡管采用上述幾類方法，可克服混凝土應(yīng)力-應(yīng)變曲線下降段“負(fù)剛度”的影響，從而提升計算效率，但深入研究這幾類方法后，可發(fā)現(xiàn)這幾類方法大多通過控制荷載增量，來實現(xiàn)消除“負(fù)剛度”的影響。而通過控制荷載增量的方法，則需重組整體剛度矩陣或者分解剛度矩陣，從而會增加模型編程實現(xiàn)的難度。這對于一些希望利用通用有限元軟件ABAQUS 的UMAT 用戶子程序，或者ANSYS 軟件的USERMAT 用戶子程序來二次開發(fā)實現(xiàn)混凝土彈塑性損傷模型的一些研究者來說，無疑也增加了其編程難度，因為不管采用UMAT 或USERMAT 用戶子程序來定義材料本構(gòu)，均只需定義出材料的應(yīng)力-應(yīng)變關(guān)系，即Jacobian矩陣，軟件即會自動組合總的剛度矩陣，而無需重組或分解剛度矩陣。顯然采用上述幾種方法會增加程序二次開發(fā)實現(xiàn)的難度，為此本文給出了一種與本文模型相適應(yīng)的迭代算法，以配合ABAQUS 軟件的UMAT 用戶子程序使用。

綜上所述，本文的目的旨在建立一個相對簡單的、能反映混凝土拉壓不等性及各向異性的特點，且能便于編程實現(xiàn)的混凝土塑性損傷本構(gòu)模型。值得一提的是，本文工作只適用于小應(yīng)變條件下的混凝土損傷-斷裂分析。

1 損傷變量定義

1.1 各向同性損傷

1) 各向同性損傷變量的定義為了應(yīng)用連續(xù)介質(zhì)損傷力學(xué)的基本原理，首先要對名義構(gòu)型和有效構(gòu)型的概念進行解釋。名義構(gòu)型也可以稱為損傷構(gòu)型，而有效構(gòu)型是一種虛構(gòu)的存在，也可稱其為未損傷構(gòu)型[20]。如圖1(a)所示，在損傷構(gòu)型中，所有類型的損傷，如空隙和裂紋等都包含在桿件中；而在有效構(gòu)型中，所有類型的損傷都從桿件中移除，如圖1(b)所示。

圖1 單軸拉伸荷載作用下混凝土柱有效構(gòu)型與名義構(gòu)型對比Fig. 1 The comparison of the nominal and effective configurations of cylindrical bar under uniaxial tension

根據(jù)文獻[30]所述，可以得到名義應(yīng)力和有效應(yīng)力之間的以下關(guān)系：

其中：

式中：A、AD和分別為如圖1 所示的桿件截面的整個面積、總損傷面積和有效面積；d為標(biāo)量損傷量，其值在0～1 之間變化；帶上劃線的物理量為有效構(gòu)型中的物理量；不帶上劃線的物理量為損傷構(gòu)型中的物理量。

式(1)所示的關(guān)系只適用于單軸應(yīng)力狀態(tài)，當(dāng)雙軸或三軸應(yīng)力狀態(tài)下，采用各向同性損傷模型時，有效應(yīng)力張量與名義應(yīng)力張量之間的關(guān)系可表示為：

理論上，損傷變量可以由式(2)確定，但由于混凝土材料的不均勻性，通過物理試驗的方法測定材料的損傷面積是一項復(fù)雜的工作[31]。為此文獻[32]提出了應(yīng)變等效假設(shè)，以期通過間接手段確定損傷變量。通過對文獻[3 - 4]的分析，利用應(yīng)變等效假設(shè)，可以得到以下關(guān)系：

利用廣義Hook 定律，考慮塑性變形的應(yīng)力-應(yīng)變關(guān)系可定義為：

式中：E為四階損傷彈性剛度張量；為四階未損傷彈性剛度張量。對于各向同性線性彈性材料，可由以下公式給出：

根據(jù)式(5)，各向同性損傷變量的間接形式可表示為：

對于一維情況，上述方程可改寫為：

式中，E為材料損傷后的彈性模量，而式(8)也即式(2)的等價形式。

2) 應(yīng)變張量的分解

由于混凝土在拉伸和壓縮狀態(tài)下表現(xiàn)出損傷不同的特性，為了充分描述由于拉伸、壓縮和循環(huán)荷載造成的混凝土拉壓損傷不同的這一特性，使用譜分解技術(shù)將應(yīng)變張量分解為正負(fù)兩部分[28]。式(9)中，上標(biāo)“+”和“-”分別表示拉伸和壓縮物理量。將總應(yīng)變張量進行分解，則可得：

式中， ε+和 ε-分別為總應(yīng)變張量的拉伸和壓縮部分。

對稱二階應(yīng)變張量的譜分解形式可表示為：

總應(yīng)變張量的正負(fù)部分可以表示為：

式中，H(x) 為階躍函數(shù)(當(dāng)x＞0 時，H(x)=1；當(dāng)x＜0 時，H(x)=0)。

將式(11)代入式(12)，可得到如下表達式：

式中，I為四階單位張量，而P+可以表示為：

由于彈性應(yīng)變張量的歸一化特征向量與總應(yīng)變張量的歸一化特征向量相同，因此，彈性應(yīng)變張量的正負(fù)部分可用下式表示為：

根據(jù)式(15)可得，彈性應(yīng)變張量可以表示為：

而根據(jù)式(5)和式(15)則可得，兩種構(gòu)型中的正應(yīng)力張量可表示為：

負(fù)應(yīng)力張量則可表示為：

根據(jù)式(16)～式(18)，可得出以下關(guān)系式：

根據(jù)文獻[3]所述，如果假設(shè)式(3)所示的有效應(yīng)力與名義應(yīng)力的關(guān)系，對拉伸和壓縮狀態(tài)均有效，則有：

從而名義應(yīng)力張量可以表示為：

式中，d+和d-分別為拉、壓各向同性損傷變量。值得注意的是，損傷變量不能簡單地分解為d+和d-，即d≠d++d-。

1.2 各向異性損傷

各向同性損傷假設(shè)在損傷演化過程中混凝土材料強度和剛度的退化在不同方向上相同，這與實際是不符的。因此，為了準(zhǔn)確描述混凝土損傷各向異性的特性，本文在假設(shè)損傷主軸與應(yīng)變主軸重合的基礎(chǔ)上，采用了二階對稱張量ω來表示損傷。由于ω是一個二階對稱張量，因此，其也可用譜分解表示為：

式中：vi(i=1,2,3)為損傷張量的歸一化特征向量；為第i個主損傷量。

由于假設(shè)損傷變量的主軸與應(yīng)變張量的主軸重合，因此，式(22)中的向量vi與式(10)中的向量ni可視為相同，則式(22)可重寫為：

由于Ei是個未知量，本文中將由基于應(yīng)變的損傷演化方程確定。關(guān)于損傷演化方程的詳細(xì)描述將在第2.2 節(jié)中給出。

另外，由于定義了二階損傷張量，因此需對式(3)進行擴展，即：

如文獻[20]所述，由式(25)得到的有效應(yīng)力張量并不是對稱的。而非對稱的有效應(yīng)力張量，不管是在本構(gòu)方程的推導(dǎo)或是在有限元的編程計算分析過程中均是不便于應(yīng)用的。為此，許多學(xué)者提出了不同的對稱化方法。本文采用文獻[33]采用的方法，對有效應(yīng)力張量進行對稱化處理。通過使用該方法，式(25)可重寫為：

或：

其中：

基于譜分解的概念，二階損傷張量可以表示為：

結(jié)合式(19)和式(32)，總的名義應(yīng)力張量可表示為：

從式(33)可以看出，損傷影響張量M可以表示為：

且：

2 彈塑性損傷模型

2.1 塑性部分

由于混凝土在拉壓荷載作用下呈現(xiàn)出不同的材料性能，本文采用文獻[34]提出的能同時考慮拉、壓塑性的屈服準(zhǔn)則。另外，對于混凝土類材料，當(dāng)計算其應(yīng)力的下降段時，會出現(xiàn)“負(fù)剛度”，從而影響收斂速度，為此本文將模型的塑性部分建立在有效構(gòu)型中。文獻[2]也表明：將塑性損傷模型的塑性部分建立在有效構(gòu)型中，可與前文所述的應(yīng)變等效假定相互配合，從而使模型的推導(dǎo)過程大大簡化，而數(shù)值實施時，其算法也更高效、更穩(wěn)定。屈服準(zhǔn)則在有效構(gòu)型中可表示如下：

式中，H+和H-定義為：

式中，符號 〈〉表示Macauley 括號，定義為：

和：

式中，Gp為塑性勢函數(shù)，其表達式將在下文給出。

結(jié)合式(40)和式(47)，式(38)可表示為：

參考文獻[3 - 4, 34]，對于混凝土類材料，壓縮狀態(tài)下塑性變形階段會產(chǎn)生明顯的體積改變，而采用非關(guān)聯(lián)的塑性流動準(zhǔn)則可很好的重現(xiàn)這一現(xiàn)象，因此本文選用非關(guān)聯(lián)流動法則來確定塑性應(yīng)變張量，即：

本文中選用Drucker-Prager 屈服函數(shù)作為塑性勢函數(shù)Gp，則：

式中， αp為膨脹系數(shù)，對于混凝土類材料，其取值范圍一般在0.2～0.3，本文取0.2。對式(53)進行求導(dǎo)，則塑性流動方向 ?Gp/可表示為：

2.2 損傷部分

由于損傷是一個不可逆過程，本文通過損傷加載函數(shù)、加載卸載條件和損傷變量的演化規(guī)律來描述模型的損傷部分。

損傷加載函數(shù)和加卸載條件可以表示為：

由于損傷和塑性變形都是導(dǎo)致混凝土非線性響應(yīng)的原因，則在單軸荷載作用下和

i可表示為[35]：

1) 拉伸損傷演化方程

式中，d0為初始損傷。參考文獻[37]的取值，本文取其值為0.063。

2) 壓縮損傷演化方程

參考文獻[3]，混凝土彈性剛度的退化在拉伸和壓縮狀態(tài)下表現(xiàn)出不同的特性(如圖2 所示)，為此需定義出不同的壓縮損傷演化規(guī)律。根據(jù)文獻[38]所述，對于混凝土類材料，在單軸壓縮載荷作用下，其塑性行為比損傷產(chǎn)生的早，即當(dāng)?shù)刃?yīng)力達到壓縮屈服強度時，材料中開始出現(xiàn)塑性行為，而當(dāng)?shù)刃?yīng)力達到峰值應(yīng)力時，則出現(xiàn)新的損傷(即壓縮狀態(tài)下≠?；诖?，單軸壓縮載荷下本文采用指數(shù)軟化規(guī)律，其峰值后的應(yīng)力-應(yīng)變關(guān)系可表示如下：

圖2 單軸荷載作用下混凝土力學(xué)特性Fig. 2 Concrete behavior under uniaxial

3 數(shù)值實施

根據(jù)本文引言部分所述，采用解耦算法，模型的塑性部分和損傷部分可以分開單獨進行求解。另外，如引言部分所述，對于混凝土類材料，在應(yīng)力-應(yīng)變曲線的下降段，由于“負(fù)剛度”的存在，計算很難收斂。本文將模型的塑性部分建立在有效構(gòu)型中，材料屈服后，則根據(jù)式(36)的屈服準(zhǔn)則，經(jīng)過迭代計算，屈服點后拉伸有效應(yīng)力-應(yīng)變曲線如圖2(a)中虛線部分所示，屈服點后壓縮有效應(yīng)力-應(yīng)變曲線則如圖2(b)中虛線部分所示。屈服點前，拉伸/壓縮有效應(yīng)力-應(yīng)變曲線與名義應(yīng)力-應(yīng)變曲線重合。從圖2 可以看出，采用將模型的塑性部分建立在有效構(gòu)型中的方法，有效應(yīng)力不存在下降段，而通過有效應(yīng)力計算Jacobian 矩陣，可有效避免“負(fù)剛度”的出現(xiàn)，從而大大提高收斂速度。但研究混凝土的損傷-斷裂狀態(tài)，實際上研究的是其損傷構(gòu)型中的力學(xué)性態(tài)，需要獲得的應(yīng)力-應(yīng)變曲線實際上是名義應(yīng)力-應(yīng)變曲線，當(dāng)有效應(yīng)力計算出以后，再根據(jù)損傷量的計算結(jié)果，即可根據(jù)式(33)計算出名義應(yīng)力。而由于采用應(yīng)變等效假設(shè)，有效構(gòu)型中的應(yīng)變與損傷構(gòu)型中應(yīng)變相等，從而名義應(yīng)力-應(yīng)變曲線即可獲得。名義應(yīng)力-應(yīng)變曲線的下降段如圖2的實線所示，為避免“負(fù)剛度”的影響，本文將名義應(yīng)力作為狀態(tài)變量，不參與Jacobian 矩陣的計算。

基于ABAQUS 平臺，采用UMAT 子程序?qū)Ρ疚哪Ｐ瓦M行二次開發(fā)，具體的模型數(shù)值求解原理如下。

3.1 塑性部分

本文采用向后歐拉隱式算法進行塑性數(shù)值積分，其包括兩個方面：彈性預(yù)測和塑性修正。增量步n和n-1之間的未知變量可以更新如下：

在 增 量 步n-1 給 出 一 組 變 量(ε(n-1),,) 和 一個應(yīng)變增量變量 Δε=Δtε˙，則式(64)為一個關(guān)于變量 (ε(n+1),,))的非線性方程組。更新變量來自前一個時間歩結(jié)束時的收斂值。非線性方程組的求解通過Newton-Raphson 迭代算法來進行。有效應(yīng)力及相關(guān)變量更新過程如下：

1) 初始化：

2) 彈性預(yù)測：

3) 檢查屈服及收斂條件：

4) 初始迭代：

5) 牛頓迭代：

① 應(yīng)力迭代

② 檢查收斂條件

如果 |l|＜TORL且＜TOLF，進入第6)步；否則進入③。

③ 更新內(nèi)變量：

④k=k+1，返回①。

6) 更新最終迭代結(jié)果：

7) 結(jié)束。

3.2 損傷量及名義應(yīng)力計算

采用解耦算法，本文中損傷量的計算采用顯式算法。損傷量獲得后，結(jié)合有效應(yīng)力的計算結(jié)果根據(jù)式(33)計算名義應(yīng)力。同時為了避免“負(fù)剛度”的影響，將名義應(yīng)力做為狀態(tài)變量進行更新，不參與剛度矩陣的迭代計算，具體算法流程如下：

4 模型驗證

為了驗證本文所提出的塑性損傷模型的適用性和有效性，本節(jié)將給出不同加載條件下由本模型計算的結(jié)果與混凝土試件試驗結(jié)果的對比，具體如下所述。

4.1 單軸壓縮加載-卸載試驗

由于塑性本構(gòu)方程是在有效構(gòu)型中定義的，因此本文中塑性材料參數(shù)的確定采用文獻[4]給出的方法。該方法可概括為：首先，通過連接圖3所示的每個卸載點(A點～E點)和重新加載點(A'點～E'點)，可以確定每個循環(huán)的損傷彈性模量E。進而可由方程=(/E)σ求得有效應(yīng)力。有效應(yīng)力得出后，繪制其與相應(yīng)應(yīng)變的關(guān)系曲線，然后可根據(jù)式(49)和式(50)確定出塑性材料參數(shù)。對于損傷材料參數(shù)，由于本文損傷演化方程是由應(yīng)力-應(yīng)變方程推導(dǎo)而來，因此，其可由式(60)和式(62)與單軸拉、壓應(yīng)力-應(yīng)變曲線擬合確定。

圖3 文獻[39]試驗所得有效構(gòu)型及名義構(gòu)型中應(yīng)力-應(yīng)變曲線Fig. 3 Experimental stress-strain curves in the effective and nominal configurations from the experimental results of literature[39]

此外，如文獻[37]所述，通常由于混凝土試件中初始微裂隙的存在，因此由單軸試驗確定的初始彈性模量E0并不是完全無損的彈性模量，其沒有考慮初始損傷d0的影響。因此，根據(jù)式(8)，可以將初始完全未損傷的彈性模量定義為=E0/(1-d0)。

使用上述方法，根據(jù)文獻[39]提供的單軸壓縮加載-卸載試驗曲線擬合的壓縮塑性和損傷材料參數(shù)數(shù)如表1 所示。數(shù)值模型計算結(jié)果與試驗結(jié)果的對比如圖4 所示。

表1 根據(jù)文獻[39]試驗結(jié)果得到的材料參數(shù)Table 1 Material constants identified from the experimental results[39]

圖4 本文模型計算的單軸壓縮加載-卸載結(jié)果與文獻[39]試驗結(jié)果對比圖Fig. 4 The model responses in uniaxial loading-unloading compression compared to experimental results presented in literature[39]

如圖4 所示，模型計算結(jié)果能很好地描述峰后區(qū)域，但單軸抗壓強度的計算結(jié)果有點偏高。其主要原因是在試驗過程中等效應(yīng)力在達到單軸抗壓強度前，新的損傷已經(jīng)開始。因此，物理試驗時材料的實際損傷是略大于初始損傷d0的，但模型計算時單軸抗壓強度是由公式=(1-d0)來獲得的，而的值是根據(jù)前述方法由試驗曲線擬合得到，其值一定的情況下，(1-d0)的值比實際值大，從而即導(dǎo)致了模型計算結(jié)果比試驗結(jié)果偏大。

4.2 單軸拉伸加載-卸載試驗

按照第4.1 節(jié)中提出的相同方法，由文獻[40]的單軸拉伸加載-卸載試驗擬合得到的材料參數(shù)如表2 所示。將文獻[40]單軸拉伸加載-卸載試驗結(jié)果與模型計算果進行對比，如圖5 所示。從圖5可以看出，模型計算結(jié)果與試驗結(jié)果吻合良好。

表2 根據(jù)文獻[40]試驗結(jié)果得到的材料參數(shù)表Table 2 Material constants identified from the experimental results of literature[40]

圖5 本文模型計算的單軸拉伸加載-卸載結(jié)果與文獻[40]試驗結(jié)果對比圖Fig. 5 The model responses in uniaxial loading-unloading tension compared to experimental results presented in literature[40]

4.3 單軸單調(diào)壓縮試驗

單軸單調(diào)壓縮加載情況，根據(jù)文獻[39]提供的試驗數(shù)據(jù)擬合的材料參數(shù)如表3 所示。模型計算結(jié)果與試驗結(jié)果的對比如圖6 所示。

表3 單軸單調(diào)壓縮試驗所用材料參數(shù)表Table 3 Material properties used for the monotonic uniaxial compressive test

圖6 本文模型計算的單軸單調(diào)壓縮加載結(jié)果與文獻[39]試驗結(jié)果對比圖Fig. 6 The model responses in monotonic uniaxial compression compared to experimental results presented in literature[39]

對于單軸單調(diào)加載情況，由于第4.1 節(jié)采用的方法已不再適用，因此單軸單調(diào)加載情況，采用origin 軟件根據(jù)最佳近似原則對試驗數(shù)據(jù)進行擬合得到材料參數(shù)，下文單軸單調(diào)拉伸加載情況，材料參數(shù)也根據(jù)此方法得到。從圖6 可以看出，模型計算結(jié)果基本能反應(yīng)實際情況。

4.4 單軸單調(diào)拉伸試驗

單軸單調(diào)拉伸加載情況，根據(jù)文獻[41]提供的試驗數(shù)據(jù)擬合得到的材料參數(shù)如表4 所示。模型計算的應(yīng)力-應(yīng)變曲線與試驗的應(yīng)力-應(yīng)變曲線對比如圖7 所示。從圖7 可以看出，模型計算結(jié)果基本能反應(yīng)實際情況。

圖7 本文模型計算的單軸單調(diào)拉伸加載結(jié)果與文獻[41]試驗結(jié)果對比圖Fig. 7 The model responses in monotonic uniaxial tension compared to experimental results presented in literature[41]

表4 單軸單調(diào)拉伸加載試驗所用材料參數(shù)表Table 4 Material parameters used for the monotonic uniaxial tensile test

4.5 雙軸單調(diào)壓縮加載試驗

將文獻[42]的雙軸壓縮試驗結(jié)果與本文模型的計算結(jié)果進行對比，如圖8 所示。材料參數(shù)如表5所示。

表5 雙軸壓縮加載試驗所用材料參數(shù)表Table 5 Material constants used for the biaxial compressive test

從圖8 可以看出，當(dāng)應(yīng)力比為σ1/σ2=-1/0時，采用本文模型得到的計算結(jié)果與試驗結(jié)果吻合很好，但當(dāng)應(yīng)力比為 σ1/σ2=-1/-1 和σ1/σ2=-1/-0.52時，應(yīng)力-應(yīng)變曲線的上升段能與試驗曲線基本吻合，而下降段則與試驗曲線存在一定的差異。究其原因，可能為試驗過程中，當(dāng)應(yīng)力比為σ1/σ2=-1/0時，材料只有一向受壓，另外兩向無約束作用，這樣受壓時，材料的微缺陷與內(nèi)部裂紋很容易就會擴展，從而導(dǎo)致有效承載面積的減小，也即損傷量的增大；當(dāng)保持 σ1與單向受壓時相同，逐漸增加 σ2，此時，與單向受壓時不同，材料在另一向由于受到約束作用，內(nèi)部的微缺陷及微裂紋不容易擴展，甚至隨著 σ2的增大會有閉合的趨勢，這就使得材料的有效承載面積減小很慢甚至?xí)龃?，從而使損傷量保持不變或減小。而多軸加載狀態(tài)下，采用本文模型計算時，雖然考慮了損傷各向異性的特點，但是計算三個主應(yīng)變方向的損傷量，各向?qū)嶋H上均是按照單軸加載狀態(tài)進行計算的，并沒有考慮三個方向之間相互作用對材料損傷的影響。因此，多軸加載狀態(tài)下，各軸之間相互作用對損傷的影響有待進行更深入的研究。

圖8 本文模型計算的單軸及雙軸單調(diào)壓縮加載結(jié)果與文獻[42]試驗結(jié)果對比圖Fig. 8 The model responses in monotonic uniaxial and biaxial compressive loading compared to experimental results reported by literature[42]

另外，需要說明的是本文利用譜分解技術(shù)將應(yīng)變張量分解為正負(fù)兩部分，并將對應(yīng)的應(yīng)力張量分解為正負(fù)兩部分。通過這種處理，不僅可以考慮拉壓損傷不同的特點，也可方便地考慮泊松比的影響。如圖8(c)所示，如果雙軸應(yīng)力比為σ1/σ2=-1/-0.52 ， ε1和 ε2均為壓縮應(yīng)變，考慮泊松比影響，則 ε3為伸長應(yīng)變。而伸長應(yīng)變可能導(dǎo)致拉伸損傷，因此在雙軸壓縮試驗中，考慮拉伸材料常數(shù)是必不可少的。

5 雙邊開口四點彎曲梁加載試驗

5.1 本文模型計算結(jié)果分析

為驗證本文模型的有效性，本節(jié)采用本文提出的混凝土彈塑性各向異性損傷本構(gòu)模型對文獻[43]所做的混凝土雙邊開口四點彎曲梁(DEN)試件的斷裂破壞情況進行數(shù)值模擬。試件厚37.5 mm，截面為400 mm×150 mm 的長方形，槽的尺寸為25 mm×5 mm。試驗采用的加載系統(tǒng)為固定支座加載系統(tǒng)，具體如圖9 所示。

圖9 文獻[43]中雙邊開口混凝土梁試件采用的固定支座加載系統(tǒng)Fig. 9 The fixed loading supports of the DEN specimen in literature[43]

基于ABAQUS 軟件的UMAT 子程序，對本文模型進行二次開發(fā)，計算過程中采用8 節(jié)點六面體線性非協(xié)調(diào)模式單元(C3D8I 單元)。為模擬固定支座對混凝土梁試件的線性約束作用，在有限元模型中，對支撐處節(jié)點的除豎向平動自由度外的其余自由度均約束。本文所用混凝土雙邊開口梁試件三維有限元網(wǎng)格如圖10 所示。

圖10 雙邊開口混凝土梁試件三維有限元網(wǎng)格Fig. 10 3D mesh of the DEN specimen

由于采用力加載方式在數(shù)值模擬過程中容易導(dǎo)致計算不收斂，因此本文計算時，采用位移加載方式在圖9 所示的加載點F1R 及F2 處分別施加-1 mm 和1 mm 的強制位移荷載，而在F2R 和F1 兩個加載點處分別施加-1/15 mm 和1/15 mm 的強制位移荷載。根據(jù)第4.3 節(jié)所述材料參數(shù)識別方法，得到的材料參數(shù)如表6 所示。

表6 雙邊開口混凝土梁試件斷裂數(shù)值模擬材料參數(shù)Table 6 Material parameters used for the DEN specimen fracture simulation

圖11～圖13 為試件完全破壞時，采用本文模型得到的雙邊開口梁試件3 個主應(yīng)變方向上損傷區(qū)分布情況。

圖14 為文獻[43]雙邊開口梁試驗裂紋最終擴展路徑。對比圖11～圖13 的損傷區(qū)分布可知，采用本文彈塑性各向異性損傷本構(gòu)模型得到的損傷區(qū)與圖14 中所示的文獻[43]試驗所得的試件斷裂區(qū)基本吻合。圖11～圖13 中3 個主應(yīng)變方向上損傷值存在大于1 小于0 的情況，分析原因可能為ABAQUS 軟件在計算過程中將積分點處損傷值外插至節(jié)點處，然后進行加權(quán)平均計算，從而造成了損傷值不在0～1 范圍內(nèi)這一現(xiàn)象。

圖11 第1 主應(yīng)變方向?qū)?yīng)的損傷區(qū)分布Fig. 11 Distribution of damage zone corresponding to the first principal strain direction

圖12 第2 主應(yīng)變方向?qū)?yīng)的損傷區(qū)分布Fig. 12 Distribution of damage zone corresponding to the second principal strain direction

圖13 第3 主應(yīng)變方向?qū)?yīng)的損傷區(qū)分布Fig. 13 Distribution of damage zone corresponding to the third principal strain direction

圖14 文獻[43]雙邊開口梁試驗裂紋擴展路徑Fig. 14 Crack modes reported in literature[43]

5.2 ABAQUS 自帶塑性損傷模型(CDP 模型)計算結(jié)果分析

為進一步驗證本文模型的有效性，本節(jié)采用ABAQUS 軟件自帶的CDP 模型對5.1 節(jié)所述的雙邊開口四點彎曲梁進行損傷計算。采用CDP 模型計算時，有限元模型、加載條件等與5.1 節(jié)完全相同。采用CDP 模型計算時所用材料參數(shù)，根據(jù)DEN 試件混凝土的強度等級并參考文獻[44]所給出的方法進行計算確定，具體材料參數(shù)如表7 和表8 所示。

表7 CDP 模型計算時所用混凝土彈塑性參數(shù)Table 7 Elasto-plastic mechanic parameters of concrete used by CDP model

表8 CDP 模型計算時所用混凝土損傷參數(shù)Table 8 Damage parameters of concrete used by CDP model

由CDP 模型計算的DEN 試件斷裂時的損傷云圖，如圖15 和圖16 所示。

從圖15～圖16 可以看出，由CDP 模型計算的拉損傷區(qū)分布，能基本反映試件中裂紋的擴展趨勢，但存在明顯的分叉裂紋，與實際情況存在一定出入，而且壓損傷區(qū)分布占較大一塊區(qū)域，裂紋路徑不明顯，與實際的裂紋分布情況也存在較大差異。

圖15 由CDP 模型計算的拉損傷云圖Fig. 15 Distribution of damage in tension calculated by CDP model

圖16 由CDP 模型計算的壓損傷云圖Fig. 16 Distribution of damage in compression calculated by CDP model

5.3 本文模型與CDP 模型計算結(jié)果的對比分析

對比本文模型及CDP 模型計算的裂紋路徑，可以看出，本文模型計算的裂紋擴展路徑更接近試驗結(jié)果，而且本文模型可以考慮損傷各向異性的特點，而CDP 模型只能考慮混凝土材料拉壓損傷不等的特點。

另外，采用本文提出的材料本構(gòu)模型及數(shù)值計算方法，只需20 min 左右即可完成整個模型的計算；而在有限元模型及加載條件完全相同的情況下，采用ABAQUS 軟件自帶的CDP 本構(gòu)模型則需約70 min 才能完成最終的計算。兩者相比，本文模型可大大提高計算效率。

經(jīng)對比分析可知，本文模型能反映混凝土損傷各向異性的特點，損傷分布情況也與試驗裂紋路徑基本吻合，而且也可大大提高計算效率，節(jié)約計算機時，因此模型可用于混凝土類材料的損傷計算分析。

6 結(jié)論

將損傷力學(xué)與經(jīng)典塑性理論相結(jié)合，提出了一種新的塑性各向異性損傷本構(gòu)模型來模擬混凝土的非線性行為。根據(jù)本文的工作，可以得出以下結(jié)論：

1) 本文建立了確定各向異性損傷變量的兩種不同的損傷演化方程。該方法不僅考慮了混凝土在拉壓荷載作用下的不同損傷機理，而且可考慮混凝土損傷各向異性的特點。

2) 在運用應(yīng)變等效假設(shè)的基礎(chǔ)上，通過解耦算法將塑性與損傷分離開來，可大大簡化本構(gòu)模型的推導(dǎo)，也可提高模型的計算效率。

3) 通過與大量試驗數(shù)據(jù)的比較，結(jié)果表明：該模型能有效地模擬混凝土在不同加載條件下的非線性特性。

4) 對DEN 試件的數(shù)值模擬表明：該各向異性塑性-損傷本構(gòu)模型與ABAQUS 軟件自帶的CDP本構(gòu)模型相比，其不僅能反映混凝土損傷各向異性的特點，而且其計算的損傷路徑更符合實際情況；從計算效率來說，其效率相比CDP 模型也更高。