雙軸受壓狀態下的高延性纖維增強水泥基復合材料本構模型

陳善富，陳靜芬，楊鳳祥，劉志明

(暨南大學力學與建筑工程學院“重大工程災害與控制”教育部重點實驗室，廣州510632)

高延性纖維增強水泥基復合材料(Engineered Cementitious Composites，簡稱ECC)是一種高性能水泥基復合材料。ECC纖維摻量低且在受拉和受壓時具有高延性的特點，同時因其良好的準應變硬化特性和多縫開裂特性，被廣泛應用于橋面板、橋面連接板、建筑防震抗震構件、混凝土保護層等實際工程中[1]。

1992年，密歇根大學的Li等[2]根據細觀力學和斷裂力學的理論，在水泥基體中添加短纖維研制出具有高延性的水泥基復合材料ECC。在我國，眾多研究ECC的學者將其命名為高延性纖維增強水泥基復合材料[3?7]。

自ECC面世以來，國內外學者對其力學性能進行了大量的試驗研究[8?17]，主要集中于單軸拉伸、單軸壓縮、彎曲和剪切等方面。試驗研究[8?17]表明，ECC因纖維橋聯作用而具有良好的拉伸，壓縮和彎曲韌性，同時具有良好的裂縫控制能力和耐久性。

有限元分析是對工程結構進行設計和研究的重要手段，要對一種材料的工程結構進行有限元分析需要該材料的本構模型。但是，目前的大型通用有限元程序(如ANSYS,ABAQUS等)的材料庫中，尚未擁有ECC的本構模型。在對ECC 結構構件進行有限元分析時，一些學者仍然采用有限元程序材料庫中現有的混凝土本構模型[5,18?19]。因此，開發一種能夠較準確地體現ECC力學行為的本構模型對推動ECC在工程結構中的廣泛應用具有重要意義。

為了對ECC結構構件進行有限元分析，許多學者提出了能夠用于預測ECC 力學行為的本構模型。Han 等[20]提出了一種基于總應變的正交各向異性二維旋轉裂縫模型模擬ECC在循環荷載下的力學響應，該模型分別定義了受拉和受壓時的滯回行為，且被吸收進OpenSees地震工程模擬軟件的材料庫中，便于對ECC構件在承受地震作用時的力學行為進行模擬。但是，二維模型存在受形狀、加載方式和邊界條件限制的缺點。為了克服上述缺點，Hung等[21]在Han 模型[20]的基礎上發展了一種ECC的三維本構模型，并編寫了本構模型的用戶自定義材料子程序UMAT嵌入LS-DYNA中描述ECC的力學行為，同時建立ECC沖切板與兩跨連續梁的三維有限元模型進行數值計算，數值計算結果與試驗結果較吻合。張曉悅等[22]基于Rots[23]旋轉裂縫理論，提出了一種ECC的三維本構模型，該模型通過斷裂能來判斷裂縫平面的旋轉，通過模擬ECC單軸受拉和四點彎曲試驗，得到了與試驗較吻合的結果。

上述三種模型都是基于彌散裂縫理論，假定ECC 由于在加載方向開裂后的材料性質與未加載方向的材料性質不同而轉化為正交各向異性的材料，較好地模擬出ECC在多軸受力狀態下的力學響應。但由于目前對ECC的壓縮試驗研究還是主要集中于單軸壓縮試驗，上述模型判定材料受壓破壞的強度參數也僅來源于從ECC單軸壓縮試驗得到的單軸應力-應變曲線中：Han 等[20]根據Kesner等[24]的ECC單軸壓縮試驗結果采用雙線性的應力-應變關系描述ECC的單軸受壓行為；Hung 等[21]根據ECC單軸壓縮試驗結果[25]采用三線性的應力-應變關系描述ECC的單軸受壓行為；張曉悅等[22]用一條非線性的Hognestad 拋物曲線[26]來描述ECC的單軸壓縮應力-應變關系的上升段，其強度參數取自Yoo等[27]的ECC單軸壓縮試驗。

目前，對于ECC受壓行為的試驗研究主要集中于單軸壓縮試驗。除了上述三種模型所引用的ECC 單軸壓縮試驗外，一些學者對不同配合比的ECC 進行了單軸壓縮試驗，得到了不同配合比的ECC 單軸壓縮應力應變曲線，并采用不同的曲線方程進行擬合：徐世烺等[28]采用二次拋物線形式描述其上升段，采用雙折線形式描述其下降段；姜海軍等[29]采用過鎮海[30]提出的混凝土本構方程描述其上升段，采用Sargin[31]提出的混凝土本構方程描述其下降段；李艷等[11]用兩個非線性的方程分別擬合了其上升段和下降段；Zhou 等[15]用線性和非線性組合的分段形式描述其上升段，用雙折線形式描述其下降段。

上述三種模型中判定ECC 受壓破壞的強度參數皆取自ECC單軸壓縮試驗，認為ECC在某一方向達到單軸抗壓強度則該方向發生受壓破壞。然而，現有的ECC雙軸試驗結果都表明[32?34]，ECC在承受雙軸壓應力時，一個方向上的抗壓強度是受其垂直方向上壓應力的大小影響的。這是由于ECC 在雙軸受壓時，一個方向上的壓應力會對其垂直方向上的壓應力所引起的側向變形產生一定程度的約束，限制了內部微裂縫的發展，從而提高了其垂直方向的抗壓強度，導致其垂直方向的應力應變關系與單軸受壓時不同。另外，文獻[33]研究了ECC 在單軸受壓和雙軸受壓時的應力-應變曲線的非線性特性，并擬合了ECC在不同雙軸壓應力比下的應力-應變全曲線方程。因此，若在本構模型中不考慮ECC在雙軸受壓時的非線性應力-應變關系和抗壓強度變化，就不能準確地預測ECC在承受不同雙軸壓應力時的力學行為。

要在本構模型中同時考慮ECC在雙軸受壓狀態下因應力比的不同而產生的非線性應力-應變關系和抗壓強度的變化，顯然各向同性本構模型無法展現上述行為。Darwin 和Pecknold[35?36]為了展現混凝土在雙軸受壓時的非線性應力-應變關系和抗壓強度變化，提出了一種混凝土的二維非線性正交各向異性本構模型。該模型假定混凝土在相互垂直的兩個加載方向上為正交各向異性材料，并采用Kupfer 雙軸強度包絡線[37]來預測不同雙軸應力比下混凝土兩個正交各向異性方向上的抗壓強度的變化，通過引入等效單軸應變的概念將混凝土單軸的應力-應變關系擴展到正交各向異性兩個方向上，能夠同時考慮混凝土在雙軸應力狀態下的強度變化和非線性應力-應變關系，從而能準確地描述混凝土在雙軸應力狀態下的力學行為。文獻[23]的試驗研究表明，在雙軸受壓狀態下ECC 兩個加載方向上的非線性應力-應變關系和抗壓強度是隨著雙軸壓應力比的變化而變化的，且變化規律與混凝土相似。所以，為了在本構模型中展現ECC上述力學行為，可以建立ECC在兩個相互垂直的加載方向上的正交各向異性本構關系，然后通過引入雙軸強度準則的方式考慮ECC兩個加載方向上的強度變化。

因此，本文的主要目的是開發一種能夠同時體現ECC 在雙軸受壓狀態下的非線性應力-應變關系和抗壓強度變化的本構模型，嵌入有限元分析程序ABAQUSv6.14中，用于預測ECC在雙軸受壓狀態下的非線性應力-應變關系和抗壓強度變化。本文建立了本構模型的控制方程，推導了本構模型的顯式數值算法，采用Fortran 語言編寫該算法的用戶自定義材料子程序UMAT，并嵌入大型通用有限元分析軟件ABAQUSv6.14中，通過對ECC試件的雙軸受壓試驗進行有限元分析驗證了所開發的本構模型的有效性。

1 雙軸受壓狀態下的ECC本構模型

本文提出的雙軸受壓狀態下的ECC本構模型能夠描述ECC在雙軸受壓時的非線性力學行為和抗壓強度變化的特點。模型包含：1)ECC在兩個相互垂直的加載方向的二維增量型正交各向異性本構關系，適合于在有限元程序中采用增量-迭代法求解非線性問題[38]；2)一條非線性的受壓應力-應變曲線用于描述ECC雙軸受壓時兩個加載方向上展現的非線性力學行為；3)一條雙軸強度包絡線用于描述ECC 雙軸受壓時兩個加載方向上的抗壓強度變化。

1.1 ECC二維增量型本構關系

針對ECC雙軸受壓時展現的非線性特性，在每一增量步內建立ECC 在兩個垂直加載方向上的二維增量型正交各向異性本構關系為：

為了將單軸應力-應變曲線擴展為雙軸應力狀態下的應力-應變曲線，這里采用文獻[35]的方法，引入等效單軸應變的概念。等效單軸應變是由實際應變消除泊松比效應得到的假定值，沒有實際的物理意義。但有必要指出的是，本研究通過引入等效單軸應變有利于采用應力-等效單軸應變的形式來分別描述ECC 在1、2方向上的受壓應力-應變關系，并通過該關系式求解任意應力狀態下1、2方向的切線模量。

式中，k為增量步總數。

式(9)表明在模型的數值算法中，1、2方向上等效單軸總應變 εui是加載過程中每一增量步等效單軸應變增量dεui的總和。

1.2 ECC受壓應力-應變曲線

圖1 ECC 壓縮應力-應變曲線示意圖Fig.1 Illustration of compressive stress-strain curve of ECC

本文將文獻[33]擬合的ECC受壓應力-應變多項式曲線方程擴展到引入等效單軸應變后各方向上的應力-等效單軸應變關系，如圖2所示，表達式為：

通過將式(11)的 σi對 εui求導可求得式(2)中1、2方向的切線模量Ei(i=1,2)：

圖2 ECC應力-等效單軸應變曲線示意圖Fig.2 Illustration of stress-equivalent uniaxial strain curve of ECC

1.3 ECC雙軸強度包絡線

圖3 ECC 雙軸強度包絡線示意圖Fig.3 Illustration of biaxial strength envelope of ECC

為了展示在本文提出的雙軸受壓狀態下的ECC本構模型(下文標記為“模型1”)中考慮雙軸受壓狀態下ECC抗壓強度變化的必要性，本文還建立了不考慮雙軸受壓狀態下ECC抗壓強度變化的本構模型(下文標記為“模型2”)，即模型2沒有引入第1.3節中所述的理論，然后將模型1和模型2應用于不同配合比的ECC試件在不同應力比下的雙軸受壓有限元分析中，并對采用模型1和模型2進行數值計算得到的結果進行對比分析。

2 數值實現方法

本文分別推導了模型1和模型2的數值算法，并采用Fortran 語言編寫了兩個模型的用戶自定義材料子程序UMAT以更新應力、1和2方向的彈性模量、剪切模量以及相關狀態變量，將子程序嵌入ABAQUSv6.14中對ECC 雙軸受壓試驗進行有限元分析。

模型1的數值算法流程如下所示(上標t表示當前增量步，t?1表示上一增量步)：

a)初始條件：

3 模型驗證

本文通過對文獻[33]中兩種配合比的ECC試件在不同應力比下的雙軸加載試驗進行有限元分析驗證模型1的有效性。

兩種配合比的ECC 試件均由普通硅酸鹽水泥、粉煤灰、精制石英砂、水、高效減水劑和PVA 纖維制作而成。ECC-I 試件的配合比為水泥∶粉煤灰∶石英砂∶水∶高效減水劑=1∶2∶0.6∶0.66∶0.024，纖維體積率為2.0%；ECC-II試件的配合比為水泥∶粉煤灰∶石英砂∶水∶高效減水劑=1∶3∶0.8∶0.92∶0.024，纖維體積率為2.0%。

通過單軸和雙軸受壓試驗獲得的材料屬性見表1。本文假定ECC-I和ECC-II試件的等效泊松比為ν=0.2；模型參數l、m和n取自文獻[35]，分別為?1.6、2.25和0.35。試件尺寸為100 mm×100 mm×100 mm，分別進行了應力比為0、0.25、0.5、0.75 和1時的雙軸受壓試驗。

本文利用有限元程序ABAQUSv6.14進行數值分析，建立了圖4所示的有限元模型分別進行ECC-I和ECC-II試件在上述五種應力比下的雙軸受壓分析。模型尺寸為100 mm×100 mm。在模型的左側和下側分別施加1、2方向的約束，在模型的右側和上側分別施加均布應力 σ1和 σ2，并令σ2/σ1=0、0.25、0.5、0.75和1以模擬ECC試件在五種應力比下的雙軸受壓試驗。為了保證在各應力比下均能達到抗壓強度，本文在ECC-I試件和ECC-II試件的1方向施加的初始荷載分別為75 MPa 和50 MPa，根據應力比可得2 方向上施加的初始荷載。如圖4所示，為了考察不同網格劃分方案對模型計算結果的影響，以ECC-I試件在單軸受壓(β=0)和等壓雙軸受壓(β=1)時為例，分別采用兩種網格劃分方案進行計算：圖4(a)1×1網格和圖4(b)5×5網格。有限單元類型采用線性減縮積分平面應力單元(CPS4R)。

表1 材料屬性[33]Table 1 Material properties

圖4 ECC試件雙軸受壓有限元模型Fig.4 Finite element model of ECCspecimens under biaxial compression

如圖5所示為ECC-I試件在應力比β=0和β=1時分別采用1×1和5×5網格劃分方案的主壓應力方向的應力-應變曲線計算結果，可以看出在應力比β=0和β=1下，兩種網格劃分方案計算得到的主壓應力方向的應力-應變曲線重疊在一起，并且在β=0時兩種網格劃分方案計算得到的極限抗壓強度和極限抗壓應變均一致，分別為?57.657 MPa和?0.298%；在β=1時也一致，分別為?66.527 MPa和?0.327%。由此可知，采用圖4所示的不同網格劃分方案并不影響數值計算結果。這是由于圖4中的模型在兩個方向上施加均布應力(σ1,σ2)的情況下，當采用1×1網格劃分方案，即只有1個單元時，顯然該單元的外力為(σ1,σ2)；當采用5×5網格劃分方案時，根據靜力平衡條件，每個單元的外力同樣為(σ1,σ2)，因此可以得到相同的計算結果。為了節省數值計算時間，后文算例統一采用1×1的網格劃分方案。

圖5 應力比β=0和β=1的ECC-I 試件采用1×1和5×5網格劃分方案的主壓應力方向的預測應力-應變曲線Fig.5 Predicted stress-strain curvesof ECC-I specimens with stress ratios β=0 and β=1 in the major compressive stress direction with 1×1 and 5×5 meshing scheme

如圖6所示為ECC-I試件在應力比β=0和β=1時分別采用三種最大時間增量步長計算得到的主壓應力方向的預測應力應變曲線，表明所采用的三種最大時間增量步長均獲得高度一致的計算結果。表2所示為ECC-I試件在應力比β=0和β=1時分別采用三種最大時間增量步長時的極限抗壓強度和極限抗壓應變計算結果，可知隨著增量步長的減小，計算得到的極限抗壓強度和極限抗壓應變越接近于精確解。在應力比β=0和β=1時，采用2×10?4s 進行計算時得到的極限抗壓強度和極限抗壓應變與精確解的誤差分別小于0.1%和1%。為了保證其余算例與精確解的誤差也能在可接受的范圍內，本文將最大時間增量步長統一設為1×10?4s。

圖6 應力比β=0和β=1的ECC-I 試件在不同最大時間增量步長下的主壓應力方向的預測應力-應變曲線Fig.6 Predicted stress-strain curvesof ECC-I specimens with stress ratios β=0 and β=1 in the major compressive stress direction under variousmaximum time increments

表2 應力比β=0和β=1的ECC-I 試件在不同最大時間增量步長下的主壓應力方向預測極限抗壓強度和極限抗壓應變Table 2 Predicted ultimate compressive strengthsand ultimate compressive strains of ECC-I specimens with stress ratios β=0 and β=1 in the major compressive stress direction under variousmaximum time increments

圖7和圖8分別為ECC-I和ECC-II 試件在雙軸受壓時不同應力比β(β=σ2/σ1, σ1<σ2<0)下主壓應力方向，即1方向的σ1?ε1曲線的試驗和數值計算結果。

由圖7和圖8可知，在單軸受壓(β=0)時，采用模型1和模型2分別進行有限元分析得到的兩種ECC試件的σ1?ε1曲線與試驗曲線吻合良好，且模型1和模型2的計算曲線重合。這是由于當試件處于單軸受壓狀態，模型1用雙軸強度準則可確定σcp1=fc′p和εcp1=ε′cp，而模型2中令σcp1=fc′p和εcp1=ε′cp，因此每一增量步中模型1和模型2的E1相等，從而計算得到了重合的σ1?ε1曲線。

對于雙軸受壓(β≠0)的情況，采用模型1計算得到的σ1?ε1曲線與試驗結果吻合良好，而采用模型2時則吻合程度較差。試驗曲線和采用模型1進行有限元分析得到的σ1?ε1曲線都表現為在加載開始階段，曲線的斜率變化幅度較小，近似為直線，隨后越接近破壞點曲線的斜率下降越明顯，到達破壞時曲線斜率為零。表明模型1 能夠較準確地預測各雙軸壓應力比下ECC試件在主壓應力方向加載過程中的剛度變化特征。而模型2在各雙軸壓應力比下的σ1?ε1曲線的初始斜率均相同，說明模型2無法準確預測ECC試件在不同雙軸壓應力狀態下的剛度變化趨勢。

圖7 各雙軸壓應力比下ECC-I 試件主壓應力方向的應力-應變曲線Fig.7 Stress-strain curvesof ECC-I specimensin the major compressive stress direction under different biaxial stress ratios

表3為ECC-I和ECC-II試件在各應力比下主壓應力方向的試驗和預測極限抗壓強度和極限抗壓應變值。由表3可知模型1預測得到的ECC-I試件的極限抗壓強度值與試驗值相比，在β=1時最小，僅為0.03%，當β=0.25時的誤差亦小于5%，模型1預測得到的ECC-II試件的極限抗壓強度值與試驗值相比，在β=1時誤差最小為?0.00%，在β=0.5時最大，僅為?6.17%。表明模型1能夠準確地預測雙軸受壓時ECC主壓應力方向的抗壓強度。

模型1預測得到的ECC-I和ECC-II試件在各應力比下的極限抗壓應變也同樣與試驗結果較吻合，但預測到得到的ECC-I試件在β=0.25和0.5時和ECC-II試件在β=0.25時的誤差分別達到了16.07%，15.34%和?14.89%，這是由于試驗數據離散性較大而式(17)假定抗壓強度σcpi與等效單軸抗壓應變εcpi為線性關系，在上述三種情況下對應的抗壓強度與等效單軸應變的坐標點剛好位于線性關系直線的兩側，故僅此三組有限元預測值與試驗結果相差較大。

圖8 各雙軸壓應力比下ECC-II 試件主壓應力方向的應力-應變曲線Fig.8 Stress-strain curves of ECC-II specimens in the major compressive stress direction under different biaxial stress ratios

由表3可知模型2預測得到的ECC-I和ECCII試件在各應力比下主壓應力方向極限抗壓強度和極限抗壓應變除了在單軸受壓(β=0)時較準確，在雙軸受壓時皆與試驗值相比相差巨大。表明模型2無法預測ECC抗壓強度在雙軸壓應力狀態下的變化。

從數值分析算例可知，模型1同時考慮了ECC在雙軸受壓時的非線性應力-應變關系和應力比的變化對抗壓強度的影響，能更準確地描述不同配合比的ECC在不同應力比雙軸受壓狀態下的力學行為。

表3 各雙軸壓應力比下ECC-I 和ECC-II 試件主壓應力方向的試驗和預測極限抗壓強度和極限抗壓應變值Table 3 Experimental and predicted ultimatecompressive strengthsand ultimatecompressive strains of ECC-I and ECC-II specimensin the major compressivestress direction under various biaxial stress ratios

4 結論

(1)建立了一個能夠同時考慮ECC在雙軸受壓時的非線性力學行為和抗壓強度變化的二維本構模型，開發了本構模型的顯式積分數值算法，采用Fortran 語言將模型算法編寫成用戶自定義材料子程序UMAT，并嵌入有限元分析軟件ABAQUS v6.14中，用于預測ECC在雙軸受壓時的非線性力學行為及抗壓強度變化。

(2)由數值計算結果可知，本文模型數值計算得到的兩種不同配合比的ECC 試件在五組雙軸壓應力比下的應力-應變曲線與試驗曲線吻合良好，并且能夠較準確地預測各應力比下兩種不同配合比的ECC試件在主壓應力方向上的抗壓強度和抗壓應變，表明本文開發的ECC 本構模型能夠較準確地預測不同配合比ECC試件在雙軸受壓狀態下的力學行為。

(3)作為對比，建立了不考慮ECC在雙軸受壓狀態下抗壓強度變化的本構模型，并對上述兩種不同配合比的ECC試件在五組應力比下進行數值分析，結果表明為了準確描述ECC雙軸受壓狀態下的力學行為，有必要考慮其雙軸受壓狀態下的抗壓強度變化。本文提出的本構模型為準確預測ECC 在雙軸受壓狀態下的非線性應力-應變關系和破壞強度提供了一種有效的分析方法。