虢成功 李杰

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

引言

混凝土材料是土木工程應(yīng)用最廣泛的建筑材料.但由于混凝土材料的高度復(fù)雜性,一些關(guān)鍵的科學(xué)問題尚未得到完善解決,混凝土材料的受力本構(gòu)關(guān)系即是其中之一[1].

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

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

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

1 連續(xù)彈塑性損傷力學(xué)框架

彈塑性損傷本構(gòu)模型可以表述為[1]

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

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

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

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

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

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

損傷能釋放率Y±為Helmholtz 自由能關(guān)于損傷的對偶量

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

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

利用損傷能釋放率建立的損傷準(zhǔn)則具有熱力學(xué)基礎(chǔ),損傷演化過程可以表述為損傷能釋放率的函數(shù)

2 微-細觀隨機斷裂模型

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

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

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

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

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

令Z±(y)=lnΔ±(y),其均值為λ±,標(biāo)準(zhǔn)差為ζ±,則 λ±,ζ±與Δ±(y) 的均值和標(biāo)準(zhǔn)差滿足如下?lián)Q算關(guān)系

隨機場Δ±(y) 的相關(guān)結(jié)構(gòu)可以用指數(shù)型相關(guān)函數(shù)描述

式中,ω±為相關(guān)尺度參數(shù),?=|y1-y2|.

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

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

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

3.1 損傷發(fā)展的多尺度能量耗散

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

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

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

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

式中,Vd為損傷體積.

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

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

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

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

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

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

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

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

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

在不同尺度上,裂紋數(shù)量應(yīng)該是裂紋驅(qū)動力的函數(shù).在損傷力學(xué)框架內(nèi),裂紋驅(qū)動力即損傷能釋放率Y

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

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

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

式中,κ為表示應(yīng)力屏蔽效應(yīng)強弱的參數(shù).

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

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

注意到損傷能釋放率與應(yīng)變之間的關(guān)系,結(jié)合式(28)可以發(fā)現(xiàn)微觀單元的耗能與加載應(yīng)變之間具有高度非線性關(guān)系.

3.2 隨機損傷表達

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

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

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

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

通?？蛇x用文獻[40]提出的經(jīng)驗塑性函數(shù)形式

3.3 考慮率效應(yīng)的動力擴展

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

以文獻[41]經(jīng)粗?；幚砗筇岢龅牡筃 型勢為例(圖5),求解不同加載速率下鍵的斷裂速率.圖中U為勢能,δ表示原子鍵拉伸長度,k為鍵的剛度,A為初始無外力作用下的勢壘高度,Eb為有效勢壘高度,x表示反應(yīng)坐標(biāo),當(dāng)x=0時表示原子鍵完好,x=1 時完全斷裂.

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

勢能U記為

鍵之間的力為

有效勢壘高度為

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

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

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

式中,δD為Dirac 函數(shù).

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

由反應(yīng)坐標(biāo)x的物理意義可知:x=0處為反射邊界,x=1處為吸收邊界.

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

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

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

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

式中,rx為該勢壘形式對應(yīng)的跨越頻率.

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

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

由此,跨越頻率

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

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

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

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

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

4 本構(gòu)關(guān)系的數(shù)值算法

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

4.1 彈性預(yù)測

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

利用損傷能釋放率的單調(diào)遞增標(biāo)量函數(shù)g(Y) 表示損傷面,分別建立受拉和受剪損傷發(fā)生準(zhǔn)則

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

4.2 塑性修正

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

n+1步的有效應(yīng)力為

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

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

4.3 損傷修正

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

5 模型驗證

5.1 單軸受拉

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

圖7 不同應(yīng)變率下單軸受拉應(yīng)力應(yīng)變曲線均值和標(biāo)準(zhǔn)差對比Fig.7 Comparison of mean and standard deviation of uniaxial tension stress-strain curves under different strain rates

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

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

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

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

5.2 單軸受壓

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

圖10 不同應(yīng)變率下單軸受壓應(yīng)力應(yīng)變曲線均值和標(biāo)準(zhǔn)差對比Fig.10 Comparison of mean and standard deviation of uniaxial compression stress-strain curves under different strain rates

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

5.3 動力強度提高因子

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

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

6 結(jié)論

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

應(yīng)當(dāng)指出:當(dāng)加載的時間尺度與熱激活的時間尺度相當(dāng)?shù)臅r候,速率過程理論所計算的逃逸速率不再準(zhǔn)確.因此,速率過程理論在高應(yīng)變速率條件下不再適用.在高應(yīng)變速率條件下的混凝土動力本構(gòu)關(guān)系,尚需進一步研究.