徐 磊 王紹洲 金永苗

（河海大學(xué) 水利水電工程學(xué)院, 南京 210098）

混凝土是典型的準(zhǔn)脆性材料,在破壞階段具有明顯的應(yīng)變局部化特征[1].現(xiàn)階段,以有限元為代表的連續(xù)介質(zhì)數(shù)值分析方法是開展混凝土材料與結(jié)構(gòu)破壞模擬的主要手段[2].研究表明,在基于傳統(tǒng)局部本構(gòu)模型的混凝土材料與結(jié)構(gòu)非線性有限元分析中,計算結(jié)果會依賴于有限元模型的單元尺寸[3],從而喪失客觀性.為解決這一問題,Ba?ant等[4]提出將積分形式非局部理論與傳統(tǒng)局部本構(gòu)模型相結(jié)合,通過引入與細(xì)觀結(jié)構(gòu)相關(guān)的材料特征長度與構(gòu)建體現(xiàn)不同位置材料間力學(xué)狀態(tài)相互影響的非局部變量,實現(xiàn)傳統(tǒng)局部本構(gòu)模型的非局部化,有效克服了基于由局部模型導(dǎo)致的有限元計算結(jié)果單元尺寸依賴性.

另一方面,通用有限元軟件ABAQUS因其強(qiáng)大的非線性分析能力得以在結(jié)構(gòu)分析中被廣泛采用[5].Xu等利用ABAQUS成功實現(xiàn)了混凝土拱壩從局部開裂到整體破壞全過程的有限元模擬[6].然而,ABAQUS中雖內(nèi)置了面向混凝土材料的損傷塑性模型[7]、彌散開裂模型[8]等本構(gòu)模型,但本質(zhì)上均屬局部本構(gòu)模型,故基于ABAQUS的混凝土結(jié)構(gòu)破壞分析難以保證分析結(jié)果的客觀性.理論上雖可通過編制非局部混凝土本構(gòu)模型的UMAT用戶材料子程序?qū)崿F(xiàn)基于ABAQUS平臺的混凝土結(jié)構(gòu)破壞的合理分析,但在具體實施中卻由于受限于ABAQUS中固有的數(shù)據(jù)傳遞方式(在UMAT接口中僅能獲取當(dāng)前積分點(diǎn)的相關(guān)信息,無法調(diào)用該積分點(diǎn)附近其他積分點(diǎn)的信息)而無法完成.

鑒于此,本文首先基于由Mazars提出的宏觀尺度混凝土損傷模型[9](簡稱為MAZARS模型),推導(dǎo)了混凝土單軸受拉條件下以單元尺寸為變量之一的應(yīng)力與變形量解析表達(dá)式,直觀地闡明了單元尺寸依賴性的產(chǎn)生原因,進(jìn)而將MAZARS模型與積分形式非局部理論相結(jié)合,給出了積分形式非局部MAZARS模型理論表達(dá).在此基礎(chǔ)上,提出了與ABAQUS數(shù)據(jù)傳遞方式相適應(yīng)的非局部MAZARS模型數(shù)值實現(xiàn)方法,并完成了相關(guān)程序編制.通過算例分析,驗證了基于ABAQUS平臺的非局部MAZARS模型數(shù)值實現(xiàn)的可行性與程序開發(fā)的正確性.研究成果不僅在一定程度上拓展了ABAQUS在混凝土材料與結(jié)構(gòu)損傷方面的分析功能,亦可為其他非局部本構(gòu)模型在ABAQUS中的數(shù)值實現(xiàn)提供借鑒與參考.

1 MAZARS本構(gòu)模型及其非局部化

在各向同性線彈性本構(gòu)模型的基礎(chǔ)上,Mazars通過引入損傷變量d在損傷力學(xué)框架內(nèi)建立針對混凝土材料的宏觀尺度本構(gòu)模型[9],應(yīng)力-應(yīng)變關(guān)系見式(1)：

式中：E0和v0分別為材料初始彈性模量與泊松比；σij和εij分別為應(yīng)力和應(yīng)變張量；σkk為體積應(yīng)力；δij是Kronecker符號.

在MAZARS模型中,損傷變量d為拉伸損傷變量dt與壓縮損傷變量dc的加權(quán)組合,見式(2)：

式中：αt與αc之和為1,分別表示拉伸損傷權(quán)重系數(shù)與壓縮損傷權(quán)重系數(shù),按式(3)、(4)計算：

式中：〈〉+為Macauley括號；β為模型參數(shù),取值范圍一般為分別為主應(yīng)變中的拉伸、壓縮部分為模型中定義的等效應(yīng)變,可由主應(yīng)變εi(i=1,2,3)值計算：

在MAZARS模型中,以材料變形歷史中產(chǎn)生的最大等效應(yīng)變?yōu)閾p傷演化方程的自變量k且令其初值為k0,拉伸損傷演化方程與壓縮損傷演化方程分別見式(6)、式(7).

式中：k0代表材料進(jìn)入損傷階段的控制閾值(即當(dāng)k=k0時,處于線彈性變形階段),其值可取為混凝土單軸拉伸狀態(tài)下的峰值拉應(yīng)變；At、Bt、Ac、Bc為模型參數(shù),對于一般混凝土材料,0.7≤At≤1.2,1≤Ac≤1.5,104≤Bt≤5×104,103≤Bc≤2×103.

損傷加載函數(shù)f(ε~,k)見式(8)：

可以看出,上述MAZARS模型屬于局部本構(gòu)模型,即一點(diǎn)的應(yīng)力狀態(tài)完全取決于該點(diǎn)的應(yīng)變狀態(tài),故基于MAZARS模型的有限元計算結(jié)果不可避免地存在單元尺寸依賴性,詳見下述.

圖1為一長度為lT的混凝土桿,左端受水平向位移約束,右端受水平向均布拉伸位移作用,處于單軸受拉狀態(tài).為觸發(fā)單軸拉伸狀態(tài)下的應(yīng)變局部化,假定在混凝土桿中部存在初始缺陷,并采用MAZARS模型描述該部位混凝土的力學(xué)特性.

圖1 混凝土單軸受拉桿件及網(wǎng)格剖分示意圖

為分析有限元計算結(jié)果的單元尺寸依賴性,假定混凝土桿沿軸線方向被均勻剖分為Ne個單元(如圖1所示),則任一單元尺寸(長度)le為

由于在單軸拉伸破壞過程中,混凝土桿將在存在初始缺陷的中部發(fā)生損傷開裂(其兩側(cè)將始終處于線彈性變形階段),故包含這一部位的單元拉應(yīng)變單調(diào)增大,而其他單元拉應(yīng)變則為先增大后減小.

依據(jù)式(1)以及中部損傷單元與其他未損傷單元拉應(yīng)力相等條件,可得：

進(jìn)一步,可將混凝土桿拉伸變形量U表示為中部損傷單元的損傷變量d、拉應(yīng)變以及單元尺寸le的函數(shù),見式(11)：

式中：Ne=lT/le.

在單軸應(yīng)力狀態(tài)下,混凝土桿應(yīng)力受控于中部損傷單元拉應(yīng)變及相應(yīng)的損傷變量,考慮式(11),其計算表達(dá)式見式(12)：

式中：k0=ft/E0,ft為混凝土單軸抗拉強(qiáng)度.

依據(jù)式(12)、(13),在給定相關(guān)MAZARS模型參數(shù)的基礎(chǔ)上,即可分別計算出在單元尺寸le特定取值下,對應(yīng)于不同量值下的混凝土桿拉伸變形、應(yīng)力量值.

圖2為在一組特定MAZARS模型參數(shù)下(E0=30 GPa,v0=0.2,At=1.0,Bt=15 000,k0=0.00012),長度為10 cm的混凝土桿在不同單元尺寸下(分別為1.43 cm,2.00 cm,3.33 cm)的單軸拉伸應(yīng)力-位移曲線.圖3為拉伸變形量值為0.001 4 cm時,不同單元尺寸條件下混凝土桿拉應(yīng)變的軸向分布.

圖2 不同單元尺寸下的應(yīng)力-位移曲線

圖3 不同單元尺寸下混凝土桿拉應(yīng)變軸向分布

從圖2可以看出,混凝土桿在損傷開裂階段的應(yīng)力變形響應(yīng)體現(xiàn)出了明顯的單元尺寸依賴性.隨著單元尺寸的減小,軟化段逐漸變陡且未能趨于收斂.從圖3可以看出,軸向單元尺寸越小,中部損傷單元的應(yīng)變集中程度越高.

為消除上述基于MAZARS模型的有限元計算結(jié)果單元尺寸依賴性,可將MAZARS模型與非局部積分理論相結(jié)合,從而實現(xiàn)MAZARS模型的非局部化.在非局部MAZARS模型中,將等效應(yīng)變作為非局部化變量,即將式(8)中的局部等效應(yīng)變替換為非局部等效應(yīng)變.依據(jù)非局部積分理論,材料點(diǎn)x處的非局部等效應(yīng)變(x)可按式(14)、(15)計算：

式中：V是對局部等效應(yīng)變求加權(quán)的空間域,其范圍受控于材料特征長度lc；s為與x相關(guān)的積分域內(nèi)的材料點(diǎn)；ψ(x,s)為非局部積分權(quán)重函數(shù),可按式(16)計算：

從式(14)可以看出,混凝土內(nèi)任一點(diǎn)的應(yīng)力不僅與該點(diǎn)自身的局部等效應(yīng)變相關(guān),還通過損傷變量d與其附近特征長度范圍內(nèi)各點(diǎn)的局部等效應(yīng)變相關(guān),從而可有效避免由局部模型導(dǎo)致的計算結(jié)果網(wǎng)格尺寸依賴性.

2 數(shù)值實現(xiàn)方法與程序開發(fā)

為便于用戶在ABAQUS平臺上開發(fā)未包括在其內(nèi)置材料庫中的本構(gòu)模型,ABAQUS提供了用戶材料子程序接口UMAT.UMAT用戶材料子程序需采用FORTRAN語言開發(fā),在從主程序獲取相關(guān)數(shù)據(jù)后,需基于給定的模型參數(shù)通過用戶開發(fā)的程序更新積分點(diǎn)應(yīng)力、雅可比矩陣,并可根據(jù)需要自定義狀態(tài)變量.

對于局部本構(gòu)模型的開發(fā),一般情況下均可通過UMAT從主程序獲取全部所需的數(shù)據(jù),如當(dāng)前積分點(diǎn)的應(yīng)變、應(yīng)變增量等.但對于積分形式的非局部MAZARS模型而言,任一積分點(diǎn)的非局部等效應(yīng)變計算不僅需要獲取其自身的應(yīng)變,還需要獲取其附近在材料特征長度范圍內(nèi)其他積分點(diǎn)的應(yīng)變數(shù)據(jù).由于ABAQUS在計算中是逐個積分點(diǎn)調(diào)用UMAT,且僅能通過UMAT獲取當(dāng)前積分點(diǎn)的應(yīng)變數(shù)據(jù),故ABAQUS固有的數(shù)據(jù)傳遞方式無法充分滿足開發(fā)非局部MAZARS模型的內(nèi)在要求.

在無法改變ABAQUS主程序與UMAT接口之間數(shù)據(jù)傳遞方式的前提下,為在ABAQUS中實現(xiàn)基于非局部MAZARS模型的混凝土損傷有限元分析,本文借鑒Pereira等[10]在LS-DYNA中開發(fā)非局部本構(gòu)模型時所采用的方法,提出一種用于不需要獲取當(dāng)前迭代步中其他相關(guān)積分點(diǎn)應(yīng)變數(shù)據(jù)的積分點(diǎn)非局部等效應(yīng)變近似方法,詳述如下.

假定需要計算非局部等效應(yīng)變的某積分點(diǎn)在前一增量步(nthIncrement)完成平衡迭代后的非局部等效應(yīng)變與局部等效應(yīng)變分別為與,并定義該積分點(diǎn)在當(dāng)前增量步((n+1)thIncrement)的非局部化比例因子kn+1為

為計算該積分點(diǎn)在當(dāng)前增量步平衡迭代過程中的非局部等效應(yīng)變,首先基于該積分點(diǎn)當(dāng)前狀態(tài)下的應(yīng)變?nèi)?按式(5)計算當(dāng)前狀態(tài)下的局部等效應(yīng)變,在此基礎(chǔ)上,利用式(17)計算出的當(dāng)前增量非局部化比例因子,計算出與相應(yīng)的該積分點(diǎn)當(dāng)前狀態(tài)下的非局部等效應(yīng)變,見公式(18)：

由于式(18)假定該積分點(diǎn)在當(dāng)前增量步中的非局部等效應(yīng)變與局部等效應(yīng)變比值和前一增量步相同,故按式(18)計算出的非局部等效應(yīng)變與按式(14)計算出的非局部等效應(yīng)變一般并不一致.但已有研究表明[11],在當(dāng)前增量步步長較小的條件下,上述二者之間的差異不大,故若分析中設(shè)置較小的增量步長,則由式(18)計算出的非局部等效應(yīng)變可視為實際非局部等效應(yīng)變的近似值.

在上述基礎(chǔ)上,本文提出了可與ABAQUS數(shù)據(jù)傳遞方式相適應(yīng)的的非局部MAZARS模型數(shù)值實現(xiàn)方法,基本思路是：(1)通過UMAT實現(xiàn)在當(dāng)前增量步的任一迭代步中各積分點(diǎn)局部等效應(yīng)變的計算與COMMON BLOCK存儲、非局部等效應(yīng)變計算、應(yīng)力與雅可比矩陣更新等；(2)通過UEXTERNALDB實現(xiàn)在前一增量步完成平衡迭代后各積分點(diǎn)非局部化比例因子的計算與COMMON BLOCK存儲；(3)通過USDFLD實現(xiàn)各積分點(diǎn)位置坐標(biāo)的獲取與COMMON BLOCK存儲；(4)通過共用COMMON BLOCK實現(xiàn)UMAT與UEXTERNALDB、UEXTERNALDB與USDFLD之間的數(shù)據(jù)傳遞.

UMAT用戶材料子程序主要計算步驟如下：

(1)獲取由ABAQUS主程序傳遞至UMAT中的非局部MAZARS模型參數(shù)、應(yīng)變、狀態(tài)變量等數(shù)據(jù),獲取非局部化比例因子(在UEXTERNALDB中計算并存儲于COMMON BLOCK中).

(2)形成彈性矩陣.

(3)計算全量應(yīng)變列陣,并通過調(diào)用ABAQUS內(nèi)置的實用程序SPRINC計算主應(yīng)變.

(4)依據(jù)式(5)計算局部等效應(yīng)變,并更新UMAT與UEXTERNALDB共用的COMMON BLOCK中該積分點(diǎn)的局部等效應(yīng)變.

(5)依據(jù)式(18)計算非局部等效應(yīng)變.

(6)依據(jù)非局部MAZARS模型的損傷加載函數(shù)進(jìn)行加載判斷,若為加載,則更新非局部等效應(yīng)變歷史最大值、損傷變量值,反之,則保持非局部等效應(yīng)變歷史最大值、損傷變量值不變.

(7)進(jìn)行應(yīng)力、雅可比矩陣DDSDDE更新.

對于任一增量步,在其完成平衡迭代后,通過調(diào)用ABAQUS內(nèi)置的實用程序UEXTERNALDB,依據(jù)式(14)與存儲在COMMON BLOCK中各積分點(diǎn)的局部等效應(yīng)變與位置坐標(biāo)(由ABAQUS調(diào)用實用程序USDFLD讀取各積分點(diǎn)位置坐標(biāo)并存儲于COMMON BLOCK中),計算該增量步結(jié)束后各積分點(diǎn)的非局部等效應(yīng)變,進(jìn)而更新各積分點(diǎn)的非局部化比例因子.

在非局部損傷有限元分析過程中,ABAQUS主程序?qū)ι鲜鯱MAT子程序的調(diào)用是在積分點(diǎn)的層次上進(jìn)行的,即在每次整體平衡迭代過程中,均需在對單元循環(huán)的基礎(chǔ)上對單元中的積分點(diǎn)循環(huán),從而逐一完成每個積分點(diǎn)的應(yīng)力、雅克比矩陣與狀態(tài)變量更新.

3 算例分析與驗證

為了驗證第2節(jié)中所提出的非局部MAZARS數(shù)值實現(xiàn)方法的可行性與程序開發(fā)的正確性,進(jìn)行了如下算例分析.

用在ABAQUS中二次開發(fā)的非局部MAZARS模型,對一混凝土單邊切口梁(見圖4)的三點(diǎn)彎曲試驗進(jìn)行數(shù)值模擬,并對切口(5 mm×50 mm)附近一定范圍(40 mm×150 mm)分別采用3種不同單元尺寸,即大尺寸5 mm、中尺寸2.5 mm、小尺寸1.25 mm.

圖4 混凝土單邊切口梁試件簡圖(單位：mm)

為便于對比分析,首先在一組給定的MAZARS模型參數(shù)取值下(E0=30 GPa,v0=0.2,k0=0.000 1,At=1,Bt=15 000,Ac=1.2,Bc=1 500,β=1),對具有不同單元尺寸的數(shù)值試件開展基于MAZARS模型的數(shù)值模擬.模擬中,試件底面左、右側(cè)支撐處分別施加豎向和水平向、豎向位移約束,頂面中部按位移加載方式施加量值為1mm的豎直向下位移荷載.

圖5給出了不同單元尺寸下模擬所得的力與位移關(guān)系曲線,圖6給出了不同單元尺寸下模擬所得的等效應(yīng)變分布.

圖5 不同單元尺寸下荷載-撓度曲線(局部)

從圖5可以看出,在線彈性變形階段,試件的受力變形響應(yīng)基本一致,但在峰后軟化階段,試件的受力變形響應(yīng)呈現(xiàn)出明顯的單元尺寸依賴性,隨著單元尺寸的減小,試件的“脆性”趨于增強(qiáng)；從圖6可以看出,基于MAZARS本構(gòu)模型計算所得的局部等效應(yīng)變分布特征明顯受控于單元尺寸,單元尺寸越小,局部等效應(yīng)變越集中,這實際上也是試件“脆性”增強(qiáng)的原因所在.以上分析表明,基于MAZARS模型的損傷有限元分析無法保證結(jié)果的客觀性.

圖6 不同單元尺寸下局部等效應(yīng)變(變形放大30倍)

針對具有不同單元尺寸的3個數(shù)值試件,分別開展基于非局部MAZARS模型的單軸拉伸模擬(材料特征長度取為15 mm,其他模型參數(shù)與局部分析相同).圖7給出了不同單元尺寸下基于非局部MAZARS模型模擬所得的荷載-撓度曲線,圖8給出了不同單元尺寸下模擬所得的非局部等效應(yīng)變分布.

圖7 不同單元尺寸下荷載-撓度曲線(非局部)

圖8 不同單元尺寸下非局部等效應(yīng)變(變形放大30倍)

從圖7可以看出,不同單元尺寸下,試件的荷載-撓度曲線基本一致,表明采用非局部MAZARS本構(gòu)模型可有效避免結(jié)構(gòu)受力變形響應(yīng)的單元尺寸依賴性；從圖8可以看出,基于本文所開發(fā)非局部MAZARS模型子程序所得的試件非局部等效應(yīng)變分布與量值亦基本一致,這實際上也是試件宏觀力學(xué)響應(yīng)基本一致的原因所在.

4 結(jié) 語

為實現(xiàn)在ABAQUS平臺上基于MAZARS本構(gòu)模型的混凝土損傷有限元客觀分析,本文結(jié)合積分形式的非局部理論,給出了非局部MAZARS本構(gòu)模型的理論表達(dá),并應(yīng)用非局部等效應(yīng)變的近似計算方法,提出了與ABAQUS數(shù)據(jù)傳遞方式相適應(yīng)的非局部MAZARS模型數(shù)值實現(xiàn)方法.在此基礎(chǔ)上,基于ABAQUS提供的UMAT、UEXTERNALDB與USDFLD二次開發(fā)接口,通過編制相關(guān)程序完成了非局部MAZARS本構(gòu)模型在ABAQUS平臺上的數(shù)值實現(xiàn).算例分析表明,基于MAZARS模型的混凝土損傷有限元分析結(jié)果呈現(xiàn)出明顯的單元尺寸依賴性,而利用本文所開發(fā)的非局部MAZARS本構(gòu)模型子程序可保證混凝土損傷破壞有限元分析成果的合理性.本文研究成果不僅拓展了ABAQUS在混凝土材料與結(jié)構(gòu)損傷方面的分析功能,亦可為其他非局部本構(gòu)模型在ABAQUS中的數(shù)值實現(xiàn)提供借鑒與參考.