考慮應變率效應的混凝土隨機損傷本構模型研究1)

虢成功 李杰

(同濟大學土木工程學院,上海 200092)

引言

混凝土材料是土木工程應用最廣泛的建筑材料.但由于混凝土材料的高度復雜性,一些關鍵的科學問題尚未得到完善解決,混凝土材料的受力本構關系即是其中之一[1].

在受力過程中,混凝土材料表現出復雜的非線性力學行為,如剛度退化、強度軟化、單邊效應、卸載后存在不可恢復變形等[2-3].同時,由于多相介質材料的隨機分布及各相材料自身力學性質的隨機性,混凝土的受力力學行為具有不可避免的隨機性[4-5].試驗表明:隨著加載速率的增加,混凝土材料強度會有不同程度的提升,相較靜力作用,動力作用下的裂紋形態往往出現分岔、彌散分布的特點[6-8].非線性、隨機性和應變率效應,構成了混凝土材料受力力學行為的三大基本特征[1].

損傷力學的出現與發展,為科學反映混凝土受力力學性質提供了基礎[9].從20世紀80年代開始,經過如文獻[10-14]等一批代表性學者的創造性工作,已經形成了確定性的損傷力學基本理論[15-18],為混凝土結構的非線性受力力學行為分析提供了科學基礎.21 世紀初以來,李杰等[19-24]引入隨機介質的基本概念,逐步建立了混凝土隨機損傷力學的基本理論.在這一研究中,將混凝土材料中的骨料、砂漿、界面過渡區等視為隨機介質,用隨機介質假定替代經典連續介質力學中的均勻性假定,從而打破了混凝土多相復合介質材料中的“相”概念,實現了單一的綜合介質材料的隨機描述.基于隨機介質和隨機損傷的觀點,從混凝土的微-細觀缺陷入手,發展了微-細觀隨機斷裂模型[19-22],為科學反映混凝土受力力學行為的非線性與隨機性提供了基礎.然而,無論是確定性損傷力學[25-28],還是隨機損傷力學[29-30],對混凝土材料的應變率效應問題,均沒有提供完整的解決方案.

鑒于上述背景,本文試圖通過新的努力,基于速率過程理論分析混凝土代表性體積元RVE 受力過程中的能量耗散關系,通過假定裂紋間的相互作用與應變率相關,在微-細觀隨機斷裂模型的基礎上建立適用于中低應變率范圍的隨機損傷本構模型,以期為混凝土結構的動力非線性分析提供基礎.

1 連續彈塑性損傷力學框架

彈塑性損傷本構模型可以表述為[1]

式中,σ,ε和εp分別為應力張量、應變張量和塑性應變張量;I,D和E0分別為4 階單位張量、損傷張量和彈性剛度張量.

借助應變等效假定[31],可以將損傷和塑性解耦

針對混凝土受拉和受壓的差異,將有效應力進行譜分解

式中,σi和n(i) 分別是有效應力張量的第i個特征值和對應的特征向量,H為Heaviside 函數.

同時,將損傷分為受拉損傷d+和受剪損傷d-兩個部分,可將損傷張量表示為損傷標量與對應投影算子的乘積之和

為了描述損傷演化這一不可逆的能量耗散過程,引入材料的Helmholtz 自由能Ψ

損傷能釋放率Y±為Helmholtz 自由能關于損傷的對偶量

Wu等[14]基于Drucker-Prager 型勢函數建立了的解析表達形式,并由此得到了損傷能釋放率的顯式表達

式中,E0為彈性模量,C0為4 階柔度張量,和分別為有效應力張量的第一和第二不變量.

利用損傷能釋放率建立的損傷準則具有熱力學基礎,損傷演化過程可以表述為損傷能釋放率的函數

2 微-細觀隨機斷裂模型

在上述確定性損傷力學的理論中,對于式(13)中的損傷演化函數G往往采用理性猜測或經驗推廣方式確定.這形成了確定性損傷力學理論中最大的缺陷.與之不同,隨機損傷力學引入隨機介質的概念反映非均質介質,開辟了綜合反映混凝土受力力學行為非線性與隨機性的科學道路.

從微觀斷裂的抽象物理模型出發,文獻[19-22]逐步建立了基于抽象物理模型的隨機損傷演化法則,形成了微-細觀隨機斷裂模型.這一模型為綜合反映準脆性材料的非線性與隨機性打開了方便之門[1].

如圖1 所示,將代表性體積單元抽象成微-細觀并聯彈簧模型.在該模型中,裂紋的產生、擴展和損傷的發展用微彈簧的隨機斷裂來表示.且每個微彈簧的應力-應變曲線假定服從彈-脆性關系.依據隨機介質假定[1,23-24],隨機損傷演化可表示為

圖1 受拉微-細觀隨機斷裂模型Fig.1 Tensile micro-meso stochastic fracture model

式中,H為Heaviside 函數; εe±為彈性應變,Δ±(y)為微彈簧斷裂應變,可以假定為平穩隨機場,其一維概率密度函數服從對數正態分布;y表示微彈簧的空間坐標.

令Z±(y)=lnΔ±(y),其均值為λ±,標準差為ζ±,則 λ±,ζ±與Δ±(y) 的均值和標準差滿足如下換算關系

隨機場Δ±(y) 的相關結構可以用指數型相關函數描述

式中,ω±為相關尺度參數,?=|y1-y2|.

利用大量試驗數據,文獻[32]識別給出了不同等級混凝土對應的具體分布參數.結合概率密度演化理論[33],可以分析計算在給定應變時的應力概率密度演化過程[34].

3 納-微-細觀隨機斷裂模型

在上述微-細觀隨機斷裂模型中,假定了微彈簧的應力-應變關系服從線彈性-斷裂關系,即在微彈簧斷裂前,不存在能量耗散,這并不符合真實的物理過程.事實上,關于微彈簧斷裂物理的研究表明[35-36],在微觀單元斷裂前,存在不同尺度的能量耗散過程.

3.1 損傷發展的多尺度能量耗散

以單軸受拉為例,考慮微-細觀隨機斷裂模型中的一個微彈簧單元.顯然,微觀單元內部存在更為細小的次級耗能單元.為此,首先考慮納觀尺度下一各向同性擴展的理想圓盤狀裂紋(圖2).由線彈性斷裂力學可知,在加載過程中,納觀裂紋擴展單位距離δa的耗能ΔQ為[37-38]

圖2 理想圓盤狀裂紋Fig.2 Ideal planar crack

式中,ηa=2πra,Ga為納觀尺度的能量釋放率.

假定存在與 δa相對應的損傷增量 δd,在這個過程中的能量耗散為[35]

式中,Vd為損傷體積.

由于式(17)和式(18)描述的是同一個過程,顯然有

裂紋的擴展速率a˙ 記為單位擴展距離 δa和凈能量勢壘跨越頻率f的乘積

如圖3 所示,f可由速率過程理論[39]給出

圖3 能量勢壘跨越過程Fig.3 Energy barrier crossing process

式中,kB為Boltzmann 常數,h為Planck 常數,T為絕對溫度,Q0為勢壘高度.

速率過程理論適用的條件之一為ΔQ?kT,所以可以對式(21)的雙曲函數線性化.

設微觀單元中各個尺度的裂紋分布可以用層級模型(圖4)描述,且各尺度裂紋存在自相似特性,則納觀裂紋總數N可以表示為[35,38]

圖4 裂紋層級模型Fig.4 Crack hierarchy model

式中,s表示從納觀到微觀的裂紋層級數,ni表示第i個層級中的裂紋數量.

在不同尺度上,裂紋數量應該是裂紋驅動力的函數.在損傷力學框架內,裂紋驅動力即損傷能釋放率Y

據此,文獻[35]從函數 F 具有尺度不變性出發,推出其具體形式為冪函數

式中,qi為第i個尺度的分形參數.由此,裂紋總數N可以表示為

事實上,由于材料力學性能、缺陷和孔隙的隨機分布,在外部作用下,準脆性材料的斷裂在高應力區首先出現.隨著載荷的增加,在已有的裂紋之間形成新的裂紋.而由于應力重分布,已有的裂紋也可能會發生閉合.當裂紋區域開始匯聚并形成局部化效應時,不同尺度的裂紋貫通、形成上一層級的裂紋.為考慮不同尺度裂紋之間的相互作用,可引入有效裂紋數Neff

式中,κ為表示應力屏蔽效應強弱的參數.

忽略納米層級裂紋個體差異性,微觀單元總能量耗散率可以用統計平均方法獲得

結合式(19)、式(20)和式(25),可求得微觀單元總能量耗散率表達為

注意到損傷能釋放率與應變之間的關系,結合式(28)可以發現微觀單元的耗能與加載應變之間具有高度非線性關系.

3.2 隨機損傷表達

將上述分析結果與微-細觀隨機斷裂模型相結合,可以建立混凝土代表性體積單元的隨機彈塑性損傷本構關系為

其中,隨機損傷演化法則為

由此,便給出了可以反映混凝土多尺度損傷演化的納-微-細觀隨機斷裂模型.

借助應變等效假定,可以實現損傷子空間和有效應力子空間的解耦.由此,可以在有效應力空間引入塑性勢函數,選用合適的流動法則計算塑性內變量的演化.為了避免在塑性子空間和損傷子空間迭代導致計算量過大的情況,也可以引入經驗塑性模型[29,40]以適當簡化計算,經驗塑性模型一般表述為

通常可選用文獻[40]提出的經驗塑性函數形式

3.3 考慮率效應的動力擴展

事實上,速率過程理論可以提供混凝土材料在不同加載速率下力學行為變化的解釋.該理論認為材料的破壞問題可以視為粒子從亞穩態的逃逸問題:少數位于勢阱底部的粒子通過熱激活后跨越一個與外力和溫度有關的能量勢壘,從而對應粒子間鍵的破壞.在動力荷載作用下,強度的提高源自加載速率和原子鍵斷裂速率的競爭機制[41-43].

以文獻[41]經粗粒化處理后提出的倒N 型勢為例(圖5),求解不同加載速率下鍵的斷裂速率.圖中U為勢能,δ表示原子鍵拉伸長度,k為鍵的剛度,A為初始無外力作用下的勢壘高度,Eb為有效勢壘高度,x表示反應坐標,當x=0時表示原子鍵完好,x=1 時完全斷裂.

圖5 倒N 型勢壘Fig.5 Inverse N potential

勢能U記為

鍵之間的力為

有效勢壘高度為

當有效勢壘為Eb=0時,可以求得臨界伸長量δc=

假定x隨時間的演化由Langevin 方程控制[41]

式中,η為黏性系數,U′(x,δ) 表示勢能對反應坐標的導數,ξ(t)為隨機力,通常可以假定為白噪聲,其前二階統計特征滿足

式中,δD為Dirac 函數.

由于隨機力的存在,需要求解x的概率密度函數p(x,t)隨時間的演化.這一演化服從如下的Fokker-Planck-Kolmogorov(FPK)方程

由反應坐標x的物理意義可知:x=0處為反射邊界,x=1處為吸收邊界.

定義F(t)為被邊界吸收的概率密度函數

顯然,F(t) 的物理意義對應t時刻鍵斷裂的概率.

利用特征值展開近似求解F(t).為此,對式(39)左右兩邊關于時間求導可得

將式(40)中概率密度函數關于時間的偏導數用一階特征值近似[41,44]

式中,rx為該勢壘形式對應的跨越頻率.

將式(41)代入式(40),化簡得到

跨越頻率可以通過計算平均首次穿越時間(mean first passage time)獲得[39].對于圖5 的勢壘形式

由此,跨越頻率

對式(33)、式(35)和式(44)進行無量綱化,得到

從圖6 可以看出:隨著加載速率的增大,F的演化變慢.由此表明,在同等拉伸長度下,加載速率高的鍵斷裂概率小.顯然,由于加載速率和鍵斷裂速率存在競爭機制,在不同加載速率條件下,不同尺度裂紋的發展速度和閉合速度也是不一樣的.依據這一分析,加載速率越高、裂紋擴展越慢,由此導致能量耗散相對減少,從而使得材料強度得到提高.進一步,依據相關試驗[45-46],可假設式(26)中反映應力屏蔽效應強弱的參數 κ 與相對應變率的對數線性相關

圖6 不同加載速率下F 隨 δ/δc的變化曲線Fig.6 Evolution of F with δ/δcunder different loading rates

通過引入了應力屏蔽參數 κ 與應變率的相關關系,前述納-微-細觀隨機斷裂模型擴展到可以適用于不同應變率條件下的分析.

值得指出,上述分析的基礎,是將式(40)中概率密度函數關于時間的偏導數取一階特征值近似.當加載速率較大時,高階特征值的影響不可忽略[41],所以本文結果僅適用于中、小速率加載情況,不適用于高速率加載情況.

4 本構關系的數值算法

多維彈塑性損傷本構關系的數值計算方法采用Simo等[47]發展的算子分離法,在有限元分析框架中予以實現.這類算法一般分為3 個步驟:彈性預測、塑性修正和損傷修正.算法的基本任務是已知上一時刻tn的基本狀態變量的條件下,給定應變增量Δε,求解tn+1時刻的基本狀態變量σn+1,

4.1 彈性預測

在彈性預測步中,假定沒有塑性演化和損傷演化,按彈性計算試算有效應力

利用損傷能釋放率的單調遞增標量函數g(Y) 表示損傷面,分別建立受拉和受剪損傷發生準則

式中,r±決定了當前損傷面大小,表示從0時刻到當前時刻最大損傷能釋放率.式(50)表明,只有當損傷能釋放率Y±超過歷史最大損傷能釋放率時,才會導致材料的進一步損傷.

4.2 塑性修正

如果n+1 步的試算損傷能釋放率大于歷史最大損傷能釋放率r±,則進行塑性修正.由于塑性應變和塑性應力間具有等價性,可以通過求解n+1步的塑性應力進行塑性修正,塑性應力定義為

n+1步的有效應力為

將式(55)代入式(54),可解得

由于當前步塑性應力的求解需要當前步的塑性因子,故式(56)為隱式方程,需迭代求解.考慮在每個時間增量步上損傷變量變化很小,可采用前進歐拉算法減少計算量.基于此,塑性應力求解可以表達為

4.3 損傷修正

事實上,在[tn,tn+1] 時間段內Ef的計算中損傷變量和損傷能釋放率都是時變的,本質屬于隱式方程,可以用Newton-Raphson 方法進行求解.但由于時間增量很小,損傷變量在這一時間增量的變化很小,故可以采用前一步的損傷變量代入.即分析中僅考慮損傷能釋放率的時變性.由于這一處理具有顯式表達,不需要進行迭代計算,具有明顯的算法優勢.

5 模型驗證

5.1 單軸受拉

取混凝土代表性體積單元進行單軸受拉數值模擬.實驗數據取自文獻[48].混凝土強度為C50,彈性模量E0=3.5×104N/mm2,試驗中加載應變率分別為=1×10-5s-1,=1×10-4s-1和=0.01 s-1.微-細觀隨機斷裂應變隨機場分布參數取值來自文獻[32],具體為λ+=4.869 6,ζ+=0.582 8,ω+=62.本文模型參數取值為=1×103,=15,=1,p+=18,ξ+=0.3,np+=3.分析中,采用概率空間剖分方法選取100個樣本,數值模擬結果與試驗數據對比見圖7.圖例中 m odel-mean和m odel-std 分別表示模型計算應力-應變曲線的均值和標準差.exp-mean和e xp-std 分別是試驗數據的均值和標準差.

圖7 不同應變率下單軸受拉應力應變曲線均值和標準差對比Fig.7 Comparison of mean and standard deviation of uniaxial tension stress-strain curves under different strain rates

在不同應變率條件下一個典型微彈簧單元載加載過程中的耗能與應變的關系曲線如圖8 所示.可見隨應變率的增加,微彈簧單元耗能速率降低,從而使得微彈簧斷裂發生滯后,導致材料強度提高.

圖8 不同應變率下一典型樣本Ef隨應變的演化Fig.8 Evolution of Efwith strain of a tipical sample under different loading rates

為驗證數值算法,選取加載速率為 1×10-5s-1的一個RVE 樣本進行了隱式和顯式兩種數值格式求解結果的對比(圖9).可見兩種解法給出的結果基本一致,說明顯式求解可以兼顧計算精度和求解效率.

圖9 隱式求解與顯式求解對比Fig.9 Comparison between implicit solution and explicit solution

5.2 單軸受壓

對單軸受壓代表性體積單元進行與上述類似的數值模擬.試驗數據取自文獻[32].混凝土強度為C50,彈性模量E0=3.5×104N/mm2,試驗中加載應變率分別為=1×10-5s-1,=1×10-4s-1和=3.5×10-2s-1.文獻[32] 識別的C50微-細觀隨機斷裂模型單軸受壓參數為λ-=7.566 8,ζ-=0.254 6,ω-=84.本文模型參數取值微:=1×10-29,=11,α0-=1,p-=26,ξ-=0.3,np-=2.分析中,采用概率空間剖分方法選取100個樣本.數值分析計算的統計特征與試驗數據對比見圖10,可見理論結果與試驗數據符合良好.

圖10 不同應變率下單軸受壓應力應變曲線均值和標準差對比Fig.10 Comparison of mean and standard deviation of uniaxial compression stress-strain curves under different strain rates

上述分析表明:本文建議模型可以良好地反映混凝土受力力學特征地非線性、隨機性與應變率效應.

5.3 動力強度提高因子

為了明確模型的適用范圍,針對單軸受拉和單軸受壓計算了動力強度提高因子DIF 并與試驗數據點進行對比.DIF 定義為動力強度與靜力強度的比值,圖11 中試驗數據點分別取自文獻[49]和文獻[50],計算中擬靜力強度取為 1×10-5s-1.從圖11 可以看出,模型不僅能反映強度隨應變率增大而提高的趨勢,而且計算的均值加減兩倍標準差范圍能覆蓋大部分的離散數據點.從DIF 的對結果來看,模型的適用應變率范圍大致在 1×10-7～10s-1,足以覆蓋地震中經常出現的應變率范圍.

圖11 DIF 對比Fig.11 Comparison of the DIF

6 結論

基于速率過程理論分析混凝土材料納-微觀裂紋擴展中的能量耗散過程,經由層級模型與微-細觀隨機斷裂模型相銜接,并引入應變加載速率對裂紋擴展速度的影響,建立了能同時反映混凝土非線性、隨機性和率敏感性的混凝土隨機損傷本構關系.通過數值模擬與實驗結果的對比分析,驗證了本文模型的有效性.本文提出模型的適用應變率范圍約為 1×10-7～10s-1,可為混凝土結構動力非線性分析、尤其是在地震作用下的分析提供基礎.

應當指出:當加載的時間尺度與熱激活的時間尺度相當的時候,速率過程理論所計算的逃逸速率不再準確.因此,速率過程理論在高應變速率條件下不再適用.在高應變速率條件下的混凝土動力本構關系,尚需進一步研究.