混凝土MAZARS本構(gòu)模型在ABAQUS中的數(shù)值實現(xiàn)及驗證

韓 峰，徐 磊，金永苗，王紹洲，崔姍姍

(1.浙江省水利水電勘測設(shè)計院，浙江 杭州 310002；2.河海大學(xué)水利水電工程學(xué)院，江蘇 南京 210098)

混凝土是由粗骨料、砂漿及界面過渡區(qū)構(gòu)成的具有復(fù)雜力學(xué)特性的準(zhǔn)脆性材料，其宏觀尺度下的開裂破壞直接歸因于細(xì)觀尺度下裂紋的萌生、擴(kuò)展與集聚，并呈現(xiàn)出典型的軟化特征和應(yīng)變局部化[1]。現(xiàn)階段，針對包括混凝土壩在內(nèi)的各類混凝土結(jié)構(gòu)的破壞過程仿真仍主要在單一宏觀尺度下進(jìn)行，而準(zhǔn)確模擬混凝土材料的損傷破壞特性則是保障仿真分析結(jié)果合理性的關(guān)鍵前提[2]。為此，在連續(xù)介質(zhì)力學(xué)框架內(nèi)，相關(guān)學(xué)者提出了包括損傷模型[3]、塑形模型[4]、彌散裂縫模型[5]等諸多宏觀尺度下的混凝土本構(gòu)模型，為開展混凝土結(jié)構(gòu)開裂破壞過程計算分析奠定了理論基礎(chǔ)。其中，由Mazars等提出的混凝土損傷模型[6](以下簡稱“MAZARS模型”)因其具有概念明確、模型參數(shù)易于標(biāo)定等優(yōu)點已得到了較為廣泛的認(rèn)可和應(yīng)用。

另一方面，通用有限元平臺ABAQUS具有強(qiáng)大的非線性分析能力，已在包括水工結(jié)構(gòu)在內(nèi)的各類工程結(jié)構(gòu)受力變形分析中得到了廣泛應(yīng)用[7]。為模擬混凝土在不同變形階段的力學(xué)特性，ABAQUS中內(nèi)置了面向混凝土材料的非線性本構(gòu)模型，主要有彌散開裂模型(Smeared Cracking Model)、塑性損傷模型(Damaged Plasticity Model)等，但包括MAZARS模型在內(nèi)的其他較為經(jīng)典的混凝土本構(gòu)模型尚未被包括在ABAQUS提供的材料庫中，這在很大程度上限制了在ABAQUS平臺上開展基于MAZARS本構(gòu)模型的混凝土結(jié)構(gòu)損傷破壞分析。

為此，本文基于MAZARS本構(gòu)模型基本理論與ABAQUS提供的用戶材料模型接口，給出了基于ABAQUS的MAZARS本構(gòu)模型數(shù)值實現(xiàn)流程，進(jìn)而通過編制UMAT子程序?qū)AZARS本構(gòu)模型進(jìn)行了二次開發(fā)。在此基礎(chǔ)上，通過開展混凝土單軸拉伸破壞過程模擬驗證了程序開發(fā)的正確性，并結(jié)合混凝土重力壩、拱壩損傷破壞分析進(jìn)行了初步應(yīng)用。

1 MAZARS本構(gòu)模型

為描述混凝土在漸進(jìn)破壞過程中材料性能的逐漸劣化，Mazars在損傷力學(xué)框架內(nèi)建立了針對混凝土材料的宏觀尺度本構(gòu)模型[6]。在各向同性線彈性本構(gòu)模型的基礎(chǔ)上，通過引入損傷變量d綜合考慮混凝土材料由各種因素導(dǎo)致的損傷，可得處于不同變形階段的混凝土應(yīng)力-應(yīng)變關(guān)系，即

(1)

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

為區(qū)別混凝土抗拉與抗壓能力的不同，在MAZARS模型中，將損傷變量d表達(dá)為拉伸損傷變量與壓縮損傷變量的加權(quán)組合，即

d=αtdt+αcdc

(2)

式中，dt與dc分別為拉伸損傷變量與壓縮損傷變量；αt與αc之和為1，分別表示拉伸損傷權(quán)重系數(shù)與壓縮損傷權(quán)重系數(shù)，其表達(dá)式如下

(3)

(4)

(5)

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

(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。

(8)

在應(yīng)用MAZARS模型模擬混凝土材料本構(gòu)行為時，需要首先確定上述8個模型參數(shù)，包括E0、v0、At、Bt、Ac、Bc、k0以及β。其中，通過開展混凝土單軸壓縮試驗即可確定E0、v0、Ac與Bc；通過開展混凝土單軸拉伸試驗可確定k0、At與Bt；而β的取值可通過開展混凝土剪切試驗獲取。MAZARS模型在一組確定的參數(shù)取值下(E0=30 GPa，v0=0.2，k0=0.000 1，At=1，Bt=15 000，Ac=1.2，Bc=1 500，β=1)給出的混凝土單軸壓縮與單軸拉伸應(yīng)力應(yīng)變曲線[6]見圖1。從圖1可以看出，MAZARS模型可以較好地模擬出混凝土在破壞階段的非線性軟化特性。

圖1 基于MAZARS模型的混凝土單軸應(yīng)力應(yīng)變曲線

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

近些年來，ABAQUS在混凝土材料與結(jié)構(gòu)破壞分析中的應(yīng)用日益廣泛[8-9]，且用戶可利用其提供的用戶材料子程序接口UMAT對所需的本構(gòu)模型進(jìn)行二次開發(fā)。本文通過編制UMAT子程序，對MAZARS本構(gòu)模型進(jìn)行了數(shù)值實現(xiàn)，子程序主要計算流程如下：

(1)讀取由ABAQUS主程序傳入的應(yīng)變列陣、應(yīng)變增量列陣、狀態(tài)變量(等效應(yīng)變歷史最大值、損傷變量)列陣、模型參數(shù)列陣等數(shù)據(jù)。

(2)依據(jù)給定模型參數(shù)中彈性模量與泊松比計算彈性矩陣。

(3)計算當(dāng)前狀態(tài)下的應(yīng)變列陣(流程(1)中應(yīng)變列陣與應(yīng)變增量列陣之和)。

(4)調(diào)用ABAQUS子程序SPRINC計算當(dāng)前狀態(tài)下的主應(yīng)變。

(5)依據(jù)式(5)計算當(dāng)前狀態(tài)下的等效應(yīng)變。

(6)依據(jù)式(8)判斷材料狀態(tài)，若等效應(yīng)變不大于等效應(yīng)變歷史最大值，則保持等效應(yīng)變歷史最大值、損傷變量值不變；若等效應(yīng)變大于等效應(yīng)變歷史最大值，則更新等效應(yīng)變歷史最大值，并依據(jù)公式(2)、(3)、(4)、(6)、(7)計算并更新?lián)p傷變量值。

(7)依據(jù)流程(6)中給出的損傷變量值，完成應(yīng)力更新，并基于流程(2)中的彈性矩陣計算雅克比矩陣DDSDDE。

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

3 算例分析與驗證

為了驗證MAZARS本構(gòu)模型數(shù)值實現(xiàn)的正確性和有效性，進(jìn)行了如下算例分析。

3.1 算例一

該算例應(yīng)用在ABAQUS中二次開發(fā)的MAZARS模型對一混凝土立方試件(150 mm×150 mm×300 mm)的單軸拉伸全過程進(jìn)行了數(shù)值模擬。采用C3D8空間等參數(shù)單元對混凝土試件進(jìn)行網(wǎng)格剖分，單元與結(jié)點數(shù)量分別為6 750與7 936，有限元網(wǎng)格見圖2。模型底面結(jié)點施加豎向約束，并在模型頂面施加均布拉伸位移，量值為0.012 cm。

圖2 混凝土試件有限元網(wǎng)格

為便于對比分析，模擬中所采用的MAZARS模型參數(shù)與繪制圖1所示曲線所采用的模型參數(shù)相同。圖3給出了模擬所得的在不同加載階段的應(yīng)力應(yīng)變散點圖(圖中亦示出了相應(yīng)的應(yīng)力應(yīng)變理論曲線)。圖4給出了在位移加載過程中損傷變量隨軸向拉伸應(yīng)變增大的變化曲線。

圖3 混凝土試件單軸拉伸應(yīng)力應(yīng)變曲線

圖4 損傷變量變化曲線

從圖3、4可以看出，采用MAZARS本構(gòu)模型子程序可準(zhǔn)確地模擬出混凝土試件單軸拉伸損傷破壞全過程，初步驗證了程序編制的正確性。

3.2 算例二

該算例開展了基于MAZARS本構(gòu)模型的重力壩模型在水壓力超載條件下的破壞過程分析，分析中假定壩基為剛性基礎(chǔ)。重力壩模型高60 m，壩底寬40 m，上游蓄水高度50 m。MAZARS模型參數(shù)同算例一。結(jié)構(gòu)模型荷載主要考慮壩體自重與上游面水壓力，并通過逐步增大水壓力分析重力壩模型的漸進(jìn)破壞過程。

定義實際施加的水壓力與初始水壓力的比值為Kw。圖5給出了Kw在不同取值(2.0及4.0)時的損傷變量分布；圖6給出了Kw在不同取值(2.0及4.0)時的主拉應(yīng)力分布。從圖5、6可以看出，損傷區(qū)首先出現(xiàn)于上游壩踵部位，隨后沿著建基面逐漸向下游擴(kuò)展，與之相應(yīng)的是壩體主拉應(yīng)力量值集中部位逐漸由上游向下游移動，原因在于損傷區(qū)上游部位損傷量值增大導(dǎo)致的應(yīng)力釋放。上述結(jié)果表明，采用本文所開發(fā)的MAZARS模型子程序可有效地模擬出水壓力超載過程中由于壩體底部拉應(yīng)力量值超過混凝土抗拉強(qiáng)度導(dǎo)致的壩體漸進(jìn)破壞。

圖5 損傷變量分布

圖6 主拉應(yīng)力分布

3.3 算例三

為進(jìn)一步驗證基于ABAQUS平臺所開發(fā)的MAZARS本構(gòu)模型子程序在實際工程結(jié)構(gòu)分析中的適用性，以我國西南地區(qū)某高為78.4 m的混凝土雙曲拱壩工程為例，考慮正常蓄水位+溫降工況，進(jìn)行了壩體損傷有限元分析。壩體-壩基系統(tǒng)的有限元網(wǎng)格見圖7。數(shù)值分析中，壩體混凝土采用MAZARS本構(gòu)模型(模型參數(shù)與前述兩個算例相同)，壩基巖體則采用彈塑性本構(gòu)模型(變形模量為6.5 GPa，泊松比為0.25，摩擦系數(shù)為0.82，凝聚力為1.0 MPa)。限于篇幅，僅給出了壩體順河向水平位移、主拉應(yīng)力以及損傷變量分布云圖，見圖8～10。

圖7 壩體-壩基系統(tǒng)有限元網(wǎng)格

圖8 壩體順河向水平位移分布

圖9 壩體主拉應(yīng)力分布

圖10 壩體損傷變量分布

從圖8～10可以看出，壩體順河向位移指向下游，最大順河向位移發(fā)生在拱冠頂部；壩體主拉應(yīng)力分布具有較好的對稱性，在水壓力作用下，壩體拉應(yīng)力區(qū)主要集中在壩體底部及中下部兩岸拱端上游側(cè)；由于壩體底部拉應(yīng)力水平較高，其上游側(cè)出現(xiàn)了受拉損傷區(qū)，但范圍不大，損傷區(qū)最大徑向開展深度約3 m。上述結(jié)果表明，基于所開發(fā)的MAZARS本構(gòu)模型子程序可有效實施針對復(fù)雜混凝土結(jié)構(gòu)的損傷有限元分析。

4 結(jié) 語

MAZARS模型是混凝土經(jīng)典本構(gòu)模型之一，但目前并未被包括在得到廣泛應(yīng)用的非線性有限元分析軟件ABAQUS提供的材料庫中。為此，本文基于MAZARS本構(gòu)模型基本理論與ABAQUS的UMAT用戶材料子程序接口，提出了MAZARS模型的數(shù)值實現(xiàn)方法，并給出了完整的子程序計算流程，進(jìn)而完成了MAZARS模型在ABAQUS中的二次開發(fā)。通過采用所開發(fā)的MAZARS模型程序?qū)炷猎嚰⒔Y(jié)構(gòu)模型以及工程結(jié)構(gòu)3種不同復(fù)雜程度分析對象開展損傷有限元計算模擬，較為系統(tǒng)地驗證了程序開發(fā)的正確性與有效性。研究成果不僅在一定程度上擴(kuò)展了ABAQUS軟件的分析功能，亦可為其他本構(gòu)模型在ABAQUS中的二次開發(fā)提供借鑒與參考。