靜力彈塑性損傷實(shí)用本構(gòu)模型及其穩(wěn)定性分析

靜力彈塑性損傷實(shí)用本構(gòu)模型及其穩(wěn)定性分析

混凝土結(jié)構(gòu)可承受影響、地震和爆炸等動(dòng)態(tài)載荷?；诨炷领o態(tài)力學(xué)性能的分析結(jié)果與實(shí)際情況之間存在較大差異。因此，考慮到適應(yīng)性率的結(jié)構(gòu)模式在研究中發(fā)揮著重要作用。在動(dòng)態(tài)負(fù)荷的作用下，混凝土材料表現(xiàn)出適應(yīng)性，即材料的損傷和非線性程度隨著適應(yīng)性的增加而降低，抗張力和抗壓強(qiáng)度隨著適應(yīng)性的增加而增加。混凝土材料的適應(yīng)性在影響負(fù)荷（如碰撞和爆炸）和地震負(fù)荷下的作用下是明顯的。研究人員很早就注意到上述現(xiàn)象,提出了多種動(dòng)力作用下的混凝土本構(gòu)模型理論,然而,由于這些理論模型本身以及數(shù)值實(shí)現(xiàn)算法的復(fù)雜性,工程實(shí)際中應(yīng)用最為廣泛的仍然是一些經(jīng)驗(yàn)?zāi)Ｐ?如根據(jù)經(jīng)驗(yàn)提高材料的單軸抗拉強(qiáng)度(或單軸抗壓強(qiáng)度)等簡(jiǎn)化處理方法.考慮應(yīng)變加載速率的混凝土動(dòng)力本構(gòu)模型通常是在率不相關(guān)本構(gòu)模型的基礎(chǔ)上做如下擴(kuò)展得到:①將率不相關(guān)本構(gòu)模型破壞面擴(kuò)大,在實(shí)際工程中應(yīng)用較多,但沒(méi)有物理背景;②將彈塑性應(yīng)力—應(yīng)變關(guān)系粘滯化,使得塑性應(yīng)力—應(yīng)變關(guān)系率相關(guān),在動(dòng)力加載過(guò)程中混凝土微粒塑性應(yīng)變的發(fā)展存在滯后;③考慮損傷演化的率相關(guān)性,考慮了動(dòng)力加載情況下混凝土材料損傷的滯后.文獻(xiàn)通過(guò)將塑性應(yīng)變和損傷變量均粘滯化,將其提出的率不相關(guān)彈塑性損傷本構(gòu)模型進(jìn)行動(dòng)力推廣,建立了相應(yīng)的率相關(guān)本構(gòu)模型,但模型較為復(fù)雜,使用不方便.文獻(xiàn)認(rèn)為:高應(yīng)變率對(duì)微裂紋發(fā)展,即損傷演化的遲滯作用導(dǎo)致了混凝土材料非線性程度降低及動(dòng)力強(qiáng)度提高.文獻(xiàn)通過(guò)對(duì)損傷變量進(jìn)行粘滯化,提出了適合大型鋼筋混凝土結(jié)構(gòu)分析的混凝土率相關(guān)彈塑性損傷本構(gòu)模型,文獻(xiàn)采用和文獻(xiàn)類(lèi)似的方法將提出的靜力彈塑性損傷本構(gòu)模型粘滯化,并通過(guò)算例表明應(yīng)變率效應(yīng)對(duì)實(shí)際結(jié)構(gòu)的動(dòng)力反應(yīng)有一定的影響.損傷力學(xué)能夠從本質(zhì)上描述混凝土材料內(nèi)部微裂紋的開(kāi)裂、發(fā)展所引起的材料宏觀非線性行為,其中損傷變量的發(fā)展演化就是為了模擬材料內(nèi)部微裂紋的開(kāi)裂發(fā)展,從而在損傷力學(xué)中可以將損傷變量粘滯化,通過(guò)延緩損傷變量的演化來(lái)模擬高應(yīng)變率對(duì)微裂紋發(fā)展的遲滯作用.筆者在以前的研究中,建議了一個(gè)混凝土率不相關(guān)彈塑性損傷本構(gòu)模型,該模型能較好地描述混凝土材料在靜力或偽靜力加載下的各種非線性行為.基于以上分析,本文將筆者建議的模型進(jìn)行動(dòng)力推廣,使其能夠描述動(dòng)力加載下混凝土材料的非線性行為.阻尼是結(jié)構(gòu)振動(dòng)過(guò)程中一種特有的能量耗散機(jī)制.目前應(yīng)用最為廣泛的是Rayleigh阻尼,Rayleigh阻尼假定阻尼矩陣C是質(zhì)量矩陣M和剛度矩陣K的線性組合,數(shù)值計(jì)算表明:單純采用質(zhì)量比例阻尼也無(wú)法減小高頻噪音的影響.考慮到質(zhì)量比例阻尼本身并無(wú)物理依據(jù),很多學(xué)者都僅僅考慮剛度比例阻尼的影響.本文將剛度阻尼應(yīng)力引入到建立的彈塑性損傷本構(gòu)模型中,使其能在材料層次考慮阻尼的影響.利用該模型對(duì)Koyna混凝土重力壩進(jìn)行了地震動(dòng)作用下的非線性時(shí)程分析,分析結(jié)果與文獻(xiàn)符合較好,表明了建立模型的有效性.本文首先簡(jiǎn)要介紹筆者建議的靜力彈塑性損傷本構(gòu)模型,然后對(duì)其進(jìn)行擴(kuò)展使之能夠考慮應(yīng)變率效應(yīng)和剛度阻尼的影響,最后通過(guò)算例從結(jié)構(gòu)位移反應(yīng)和數(shù)值穩(wěn)定性兩個(gè)方面探討應(yīng)變率效應(yīng)和剛度阻尼對(duì)結(jié)構(gòu)動(dòng)力反應(yīng)的影響.1本構(gòu)模型的開(kāi)發(fā)損傷力學(xué)能描述材料的各種非線性本構(gòu)行為,如剛度退化應(yīng)變軟化等,目前用結(jié)合塑性力學(xué)和損傷力學(xué)的彈塑性損傷本構(gòu)模型來(lái)模擬混凝土材料得到了較為廣泛的認(rèn)同.由于塑性力學(xué)本身的復(fù)雜性,目前用大多彈塑性損傷本構(gòu)模型數(shù)值處理復(fù)雜計(jì)算,效率低穩(wěn)定性不好,且涉及參數(shù)較多,一般只給出了特定情況下的材料參數(shù)取值,因此較難在實(shí)際工程中應(yīng)用.筆者在以前的研究中提出了一個(gè)率不相關(guān)彈塑性實(shí)用本構(gòu)模型,通過(guò)對(duì)彈性Helmholtz自由能進(jìn)行修正,使得模型能較為準(zhǔn)確地模擬混凝土材料在雙軸、三軸加載下的非線性本構(gòu)行為,其彈塑性部分使用經(jīng)驗(yàn)彈塑性模型降低模型的數(shù)值復(fù)雜性,下面首先對(duì)筆者開(kāi)發(fā)的靜力彈塑性損傷本構(gòu)模型做簡(jiǎn)要介紹.1.1應(yīng)力球量擬合損傷能量釋放率采用如下公式計(jì)算:Y±=-?ψe?d±=12γ|ˉσ±||ˉσ|(ˉσ∶ˉD∶ˉσ-19((ˉσm)2δ∶ˉD∶δ))(1)式中:標(biāo)記“∶”為二階縮并積;Y±表示正、負(fù)損傷能量釋放率;ψe表示修正后的彈性Helmholtz自由能;d±表示受拉、受壓損傷變量;ˉσ表示有效應(yīng)力張量;ˉσ±為ˉσ的正、負(fù)分量;ˉD為有效柔度張量;ˉσm表示有效應(yīng)力球量;‖X‖=X∶X;γ為本文引入的折減系數(shù);δ表示向量(1,1,1,0,0,0).考慮到雙軸受壓加載下(?σ3<?σ2<?σ1=0,其中:?σ1??σ2??σ3分別為第一、第二、第三主應(yīng)力)材料強(qiáng)度的提高和加載應(yīng)力比密切相關(guān),另外不同應(yīng)力狀態(tài)下,應(yīng)力偏量大小相等時(shí)(不考慮方向),應(yīng)力比能在一定程度上反映應(yīng)力球量的大小,本文將γ取為應(yīng)力比的表達(dá)式:γ=1-0.265(1-Η(?σ2))√|?σ2/?σ3|(2)式中,H(x)為Heaviside函數(shù).如果式(2)計(jì)算得γ<0,則取γ=0.式(1)中:ˉσ+=Ρ+∶ˉσ(3a)ˉσ-=ˉσ-ˉσ+=Ρ-∶ˉσ(3b)式中:對(duì)稱(chēng)張量P+、P-分別為ˉσ的正、負(fù)投影張量,表示為ˉσ的特征值^-σi和特征向量Pi的函數(shù):Ρ+=∑iΗ(^-σi)(Ρii?Ρii)(4a)Ρ-=Ι-Ρ+(4b)式中:標(biāo)記“?”為張量積,Pii=Pi?Pi為二階對(duì)稱(chēng)張量;I為四階一致性張量.求出損傷能量釋放率后,采用式(5)可計(jì)算出受拉、受壓損傷變量d+、d-:d±=1-11+(a±[Y±-Y±0])b±(5)式中:a±,b±,Y±0均為材料參數(shù).加權(quán)損傷變量:d=1-(1-sd+)(1-d-)(6)式中:s=|ˉσ+|/|ˉσ|=√ˉσ+∶ˉσ+/√ˉσ∶ˉσ.1.2經(jīng)典塑性理論的解釋文獻(xiàn)給出如下塑性應(yīng)變的經(jīng)驗(yàn)計(jì)算公式:dεp=βpEΗ(dd-)?Ιˉσ∶dε?E-10∶Ιˉσ(7)式中:εp為塑性應(yīng)變,dεp表示對(duì)εp取微分;E表示材料彈性模量;E0表示材料初始彈性剛度張量;Ιˉσ=ˉσ/|ˉσ|為單位有效應(yīng)力張量;參數(shù)βp控制塑性應(yīng)變的大小,βp∈(0,1),0表示線彈性,1表示理想彈塑性.本文給出如下公式計(jì)算βp取值:βp=0.1+0.45(1-Η(?σ2)√|?σ2/?σ3|(8)以上經(jīng)驗(yàn)計(jì)算公式(7)—(8)大大簡(jiǎn)化了塑性應(yīng)變的計(jì)算,而經(jīng)典塑性理論的流動(dòng)法則、塑性硬化法則等都沒(méi)有得到體現(xiàn).對(duì)這種非正統(tǒng)做法的解釋是:①本文模型的建立是為了分析大型實(shí)際工程,這種方法計(jì)算精度雖然不如經(jīng)典塑性力學(xué),但能大大提高計(jì)算效率并增強(qiáng)模型的數(shù)值穩(wěn)定性;②式(7)假設(shè)彈性應(yīng)變方向E0-1∶σˉ為塑性流動(dòng)方向,βpE?Η(dd-)?Ισˉ∶dε?/|σˉ|可視為塑性流動(dòng)因子.式(7)只考慮受壓塑性應(yīng)變,為了考慮受拉塑性應(yīng)變,本文采用如下表達(dá)式:dεp=E(βpΗ(dd-)+βtΗ(dd+))?Ισˉ∶dε?E0-1∶Ισˉ(9)從式(9)中可以看出,參數(shù)βt和βp分別控制受壓、受拉塑性變量發(fā)展,當(dāng)受拉損傷發(fā)展即dd+>0時(shí),dεp>0,塑性應(yīng)變也得到發(fā)展.為了進(jìn)一步提高模型的計(jì)算效率和穩(wěn)定性,在工程實(shí)用中通常將本構(gòu)曲線下降段變得緩和.這里引入受拉塑性應(yīng)變?cè)谝欢ǔ潭壬暇褪菫榱颂岣吣Ｐ头€(wěn)定性,圖1給出了βt對(duì)混凝土材料單軸受拉曲線的影響,可見(jiàn)隨著βt的增大,材料受拉骨架曲線變得更加緩和,且損傷變量的發(fā)展也變得緩慢.在從后文算例可以看出,這樣處理能在一定程度上提高本構(gòu)模型的穩(wěn)定性.2求解算法仿真應(yīng)變率無(wú)關(guān)模型損傷能量釋放率閥值直接取為損傷能量釋放率的歷史最大值,對(duì)于應(yīng)變率相關(guān)模型,則需要對(duì)損傷能量釋放率進(jìn)行Perzyna粘滯規(guī)則化,文獻(xiàn)建議:r˙±=μ±?±(Y±,r±)≥0(10a)?±(Y±,r±)=Y0±(?Y±-r±?r±)a±(10b)式中:r±表示拉、壓損傷能量釋放率閥值;r˙±=dr±/dt表示將r±關(guān)于時(shí)間求導(dǎo);μ±為模型參數(shù).采用無(wú)條件穩(wěn)定,且具有二次收斂精度的梯形算法求解式(10),可得:rn+1=rn+Δtμ?±(Yn+1/2,rn+1/2)(11)式中:Yn+1/2=0.5(Yn+Yn+1);rn+1/2=0.5(rn+rn+1).利用Newton-Raphson算法定義f(r)=-r+rn+Δtμ?±(Yn+1/2,rn+1/2)=0,并給出迭代式:rn+1i+1=rn+1i-f(rn+1i)/f′(rn+1i)(12)圖2給出了不同應(yīng)變率作用下混凝土材料的單軸受拉、單軸受壓應(yīng)力—應(yīng)變?nèi)€數(shù)值模擬結(jié)果,圖3給出了不同ε˙(應(yīng)變率)加載條件下,混凝土材料抗拉極限承載力、抗壓極限承載力數(shù)值模擬結(jié)果與相應(yīng)靜力強(qiáng)度比值(在圖3中用ρ表示).作為對(duì)比,圖3中還給出了實(shí)驗(yàn)結(jié)果,可以看出計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果符合很好.材料參數(shù)為:E=31000MPa,泊松比ν=0.2、單軸抗拉強(qiáng)度f(wàn)t=3.38MPa、單軸抗壓強(qiáng)度f(wàn)c=27.6MPa.對(duì)模型其他參數(shù),采用單軸加載進(jìn)行標(biāo)定,其取值分別為:μ+=2,μ-=10000,a-=a+=4.模型計(jì)算混凝土材料在應(yīng)變速率分別為0.00001、0.0001、0.001、0.01和0.1時(shí)的單軸受壓應(yīng)力—應(yīng)變滯回曲線如圖4所示.3ca多張性材料的粘滯阻尼應(yīng)力為了在材料本構(gòu)模型中直接考慮剛度阻尼的影響,本文在建議的彈塑性損傷本構(gòu)模型中引入剛度阻尼應(yīng)力.在瑞利阻尼中只考慮剛度阻尼,則無(wú)損材料阻尼力表達(dá)式為:σˉvis=βkE0∶ε˙?其中βk為剛度組合系數(shù).Cauchy粘滯阻尼應(yīng)力σvis可表示為:σvis=(1-d)σˉvis=βk(1-d)E0∶ε˙(13)彈塑性損傷本構(gòu)關(guān)系為:σ=(1-d)σˉσˉ=E0∶εe=E0∶(ε-εp)}(14)式中,εe,εp分別為彈性和塑性應(yīng)變.則總應(yīng)力可表示如下:σtot=σ+σvis=(1-d)(σˉ+σˉvis)(15)4應(yīng)變率相關(guān)模型在有限元隱式計(jì)算中,通常采用Newton-Raphson算法求解平衡方程,此時(shí)需要用到一致切線模量,參照文獻(xiàn),給出本文模型一致切線模量的推導(dǎo)過(guò)程.?σtot?ε=(1-d-?d?εe∶E0-1∶σˉ)?σˉ?ε+dσvisdε(16)其中:dσvisdε=βk(1-d)E01dt.將式(7)代入式(14)并對(duì)ε求導(dǎo)可得:?σˉ?ε=E0-E(βpΗ(dd-)+βtΗ(dd+))?Ισˉ?Ισˉ?(17)將式(6)對(duì)εe求導(dǎo)可得:?d?εe=s(1-d-)?d+?εe+(1-sd+)?d-?εe(18)對(duì)于應(yīng)變率相關(guān)模型,在式(5)中將Y±,Y±0替換成r±,r±0則有:{?d±?εe=?d±?r±?r±?εe?d±?r±=a±b±(a±[r±-r0±])(b±)-1(1+(a±[r±-r0±])b±)2(19)根據(jù)式(11)—(15)給出的算法可得dr±=λdY±.其中λ=Δtμar0?Yn+1/2/rn+1/2-1?a-1rn+1/22rn+1/22+Δtμar0?Yn+1/2/rn+1/2-1?a-1Yn+1/2.?r±?εe=λ?Y±?εe(20)?Y±?εe=|σ±||σ|(Eˉ∶εe-19(εmeδ∶Eˉ∶δ))(21)式中:Eˉ表示有效剛張量;εem表示彈性應(yīng)變球量.將式(17)—(21)代入式(16)即可求解切線剛度模量.5接觸模型分析應(yīng)用本文建立的本構(gòu)模型,對(duì)Koyna混凝土重力壩(圖5)進(jìn)行了數(shù)值模擬,混凝土材料參數(shù)取和文獻(xiàn)一致,分別為:密度ρ0=2643kg·m-3,彈性模量E=31027MPa,泊松比ν0=0.2,βp=0.5,單軸抗拉強(qiáng)度f(wàn)t=2.9MPa,單軸抗壓強(qiáng)度f(wàn)c=24.1MPa.材料阻尼取其第一振型臨界阻尼的3%.在分析中,壩體和基礎(chǔ)之間假定為剛性連接,混凝土壩體部分網(wǎng)格有760個(gè)4節(jié)點(diǎn)平面應(yīng)力單元組成,地震動(dòng)引起的水對(duì)壩體的壓力則通過(guò)附加質(zhì)量法模擬.Koyna大壩在自重及水壓力下的自振頻率見(jiàn)表1.從圖6中可以看出本文模型計(jì)算位移反應(yīng)和文獻(xiàn)分析結(jié)果符合較好,本文模型分析結(jié)果略大.本文在分析中發(fā)現(xiàn)如果不考慮剛度阻尼,分析到后期容易出現(xiàn)不收斂.加入剛度阻尼后,本構(gòu)模型的穩(wěn)定性得到顯著提高,可見(jiàn)阻尼對(duì)結(jié)構(gòu)的能量耗散作用能改善隱式動(dòng)力計(jì)算的收斂性.表2給出不同情況下的計(jì)算時(shí)間,從中可知加入受拉塑性應(yīng)變也能在一定程度上縮短計(jì)算時(shí)間,提高模型的穩(wěn)定性,而對(duì)計(jì)算結(jié)果影響不大,本文認(rèn)為這是因?yàn)榭紤]受拉塑性應(yīng)變后材料受拉下降段變得緩和,且受拉損傷演化得到減緩從而提高模型數(shù)值穩(wěn)定性.在實(shí)際分析中,混凝土的受拉性能對(duì)分析結(jié)果影響很小,但對(duì)計(jì)算收斂性的影響卻很大,受拉塑性應(yīng)變的引入能有效提高模型的穩(wěn)定性和收斂性,從而縮短計(jì)算時(shí)間.從圖7可知,考慮應(yīng)變率效應(yīng)后結(jié)構(gòu)頂層位移反應(yīng)峰值比不考慮應(yīng)變率效應(yīng)略大,這點(diǎn)和文獻(xiàn)分析結(jié)果一致.從表2可以看出,考慮應(yīng)變率效應(yīng)后計(jì)算時(shí)間進(jìn)一步縮短,可見(jiàn)考慮率應(yīng)變?cè)诮档筒牧蠐p傷程度的同時(shí)也提高了模型的計(jì)算效率和穩(wěn)定性.圖8為三種情況下(βt均取0.5)模型損傷分布圖,在不考慮剛度阻尼和應(yīng)變率效應(yīng)情況下加載到后期,由于模型損傷過(guò)度發(fā)展導(dǎo)致計(jì)算不收斂,本文在這里給出4.727s時(shí)的損傷分布圖,其他兩種情況均取加載結(jié)束時(shí)(10s)的損傷分布.從圖8中可以看出剛度阻尼、應(yīng)變率效應(yīng)對(duì)材料損傷發(fā)展有顯著影響,都能延緩材料損傷發(fā)展程度從而提高計(jì)算穩(wěn)定性,如果在實(shí)際結(jié)構(gòu)的計(jì)算分析中不考慮其影響會(huì)導(dǎo)致結(jié)構(gòu)損傷過(guò)快發(fā)展導(dǎo)致計(jì)算結(jié)果失真,有時(shí)甚至?xí)霈F(xiàn)計(jì)算不收斂.6基于約束剛度阻尼的接觸本文在提出的靜力彈塑