混凝土MAZARS模型的非局部化及其數值實現與驗證

徐 磊 王紹洲 金永苗

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

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

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

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

1 MAZARS本構模型及其非局部化

在各向同性線彈性本構模型的基礎上,Mazars通過引入損傷變量d在損傷力學框架內建立針對混凝土材料的宏觀尺度本構模型[9],應力-應變關系見式(1)：

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

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

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

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

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

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

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

可以看出,上述MAZARS模型屬于局部本構模型,即一點的應力狀態完全取決于該點的應變狀態,故基于MAZARS模型的有限元計算結果不可避免地存在單元尺寸依賴性,詳見下述.

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

圖1 混凝土單軸受拉桿件及網格剖分示意圖

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

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

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

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

式中：Ne=lT/le.

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

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

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

圖2為在一組特定MAZARS模型參數下(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)的單軸拉伸應力-位移曲線.圖3為拉伸變形量值為0.001 4 cm時,不同單元尺寸條件下混凝土桿拉應變的軸向分布.

圖2 不同單元尺寸下的應力-位移曲線

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

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

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

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

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

2 數值實現方法與程序開發

為便于用戶在ABAQUS平臺上開發未包括在其內置材料庫中的本構模型,ABAQUS提供了用戶材料子程序接口UMAT.UMAT用戶材料子程序需采用FORTRAN語言開發,在從主程序獲取相關數據后,需基于給定的模型參數通過用戶開發的程序更新積分點應力、雅可比矩陣,并可根據需要自定義狀態變量.

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

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

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

為計算該積分點在當前增量步平衡迭代過程中的非局部等效應變,首先基于該積分點當前狀態下的應變全量,按式(5)計算當前狀態下的局部等效應變,在此基礎上,利用式(17)計算出的當前增量非局部化比例因子,計算出與相應的該積分點當前狀態下的非局部等效應變,見公式(18)：

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

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

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

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

(2)形成彈性矩陣.

(3)計算全量應變列陣,并通過調用ABAQUS內置的實用程序SPRINC計算主應變.

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

(5)依據式(18)計算非局部等效應變.

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

(7)進行應力、雅可比矩陣DDSDDE更新.

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

在非局部損傷有限元分析過程中,ABAQUS主程序對上述UMAT子程序的調用是在積分點的層次上進行的,即在每次整體平衡迭代過程中,均需在對單元循環的基礎上對單元中的積分點循環,從而逐一完成每個積分點的應力、雅克比矩陣與狀態變量更新.

3 算例分析與驗證

為了驗證第2節中所提出的非局部MAZARS數值實現方法的可行性與程序開發的正確性,進行了如下算例分析.

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

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

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

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

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

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

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

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

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

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

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

4 結 語

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