武守信魏吉瑞楊舒蔚

?(西南交通大學(xué)土木工程學(xué)院橋梁系，成都610031)?(交通隧道工程教育部重點實驗室，成都610031)

固體力學(xué)

基于能量等效原理的應(yīng)變局部化分析:II.有限元解法1)

武守信?,?,2)魏吉瑞?,?楊舒蔚?,?

?(西南交通大學(xué)土木工程學(xué)院橋梁系，成都610031)?(交通隧道工程教育部重點實驗室，成都610031)

以非局部塑性理論為基礎(chǔ)，應(yīng)用狀態(tài)空間理論,通過局部和非局部兩個狀態(tài)空間的塑性能量耗散率等效原理，提出了一種求解應(yīng)變局部化問題的新方法，以得到與網(wǎng)格無關(guān)的數(shù)值解.針對二維問題的屈服函數(shù)和流動法則導(dǎo)出了求解非局部內(nèi)變量的一般方程，并提出了在有限元環(huán)境中求解應(yīng)變局部化問題的應(yīng)力更新算法.為了驗證所提出的方法，對1個一維拉桿和3個二維平面應(yīng)變加載試件進行了有限元分析.數(shù)值結(jié)果表明，塑性應(yīng)變的分布和載荷--位移曲線都隨著網(wǎng)格的變小而穩(wěn)定地收斂，應(yīng)變局部化區(qū)域的尺寸只與材料內(nèi)尺度有關(guān)，而對有限元網(wǎng)格的大小不敏感.對于一維問題，當(dāng)有限元網(wǎng)格尺寸減小時，數(shù)值解收斂于解析解.對于二維剪切帶局部化問題，數(shù)值解隨著網(wǎng)格尺寸的減小而穩(wěn)定地向唯一解收斂.當(dāng)網(wǎng)格尺寸減小時，剪切帶的寬度和方向基本上沒有變化.而且得到的塑性應(yīng)變分布和網(wǎng)格變形是平滑的.這說明，所提方法可以克服經(jīng)典連續(xù)介質(zhì)力學(xué)模型導(dǎo)致的網(wǎng)格相關(guān)性問題，從而獲得具有物理意義的客觀解.此模型只需要單元之間的位移插值函數(shù)具有C0連續(xù)性，因而容易在現(xiàn)有的有限元程序中實現(xiàn)而無需對程序作大的修改.

應(yīng)變局部化，非局部塑性，內(nèi)尺度，網(wǎng)格相關(guān)性，有限元

引言

巖土材料在接近破壞或斷裂之前會出現(xiàn)高度局部化的塑性變形.這一現(xiàn)象不僅出現(xiàn)在脆性材料的破壞過程中，也出現(xiàn)在金屬和其他聚合物材料的破壞過程中[111].一般這些材料在破壞之前，局部化的塑性變形使得材料在斷裂之前呈現(xiàn)一條或幾條帶狀區(qū)域，通常叫做剪切帶.在剪切帶形成和擴展的過程中還伴隨有應(yīng)變軟化，即材料在加載方向的承載力隨著變形的增加而減小.對剪切帶局部化的數(shù)值模擬可以預(yù)測材料破壞的位置、最大變形、以及材料在斷裂之前承載力隨變形增大逐漸減小的全過程.這對于理解這類材料的破壞機理、預(yù)測巖土邊坡和路基的失穩(wěn)破壞、提高建筑物基礎(chǔ)和隧道開挖過程中的安全性具有重要的工程意義.

然而，基于經(jīng)典連續(xù)介質(zhì)力學(xué)理論的剪切帶局部化數(shù)值模擬會導(dǎo)致具有網(wǎng)格相關(guān)性的結(jié)果[1216].經(jīng)過30多年的研究，目前廣泛認為經(jīng)典連續(xù)介質(zhì)力學(xué)采用逐點描述材料行為的方式，沒有包含材料的內(nèi)尺度信息，導(dǎo)致描述應(yīng)變軟化過程的微分方程失去強橢圓性[1718]，而相應(yīng)的邊值問題變得不適定.因此經(jīng)典連續(xù)介質(zhì)力學(xué)不適合描述高度不均勻的變形.

為了克服網(wǎng)格相關(guān)性，學(xué)者們提出了很多理論或模型.其中非局部塑性模型通過在屈服函數(shù)中引入材料內(nèi)尺度，克服了連續(xù)介質(zhì)力學(xué)的根本缺陷.然而，一般非局部塑性模型在應(yīng)用中仍然出現(xiàn)了與網(wǎng)格相關(guān)的數(shù)值模擬結(jié)果[19].雖然Vermeer和Brinkgreve[20]、Str¨omberg和Ristinmaa[21]以及L¨u等[2223]先后提出和發(fā)展了混合非局部塑性模型并得到了與網(wǎng)格無關(guān)的解，但到目前為止，對于一般非局部塑性模型導(dǎo)致數(shù)值結(jié)果與網(wǎng)格相關(guān)的內(nèi)在原因仍缺乏進一步的研究.

一般非局部塑性模型由Ba?zant和Lin[24]于1988年首先提出.它源于Eringen[25]的非局部連續(xù)介質(zhì)力學(xué)理論.Nilsson[26],Borino等[27],de Sciarra[28]從熱動力學(xué)和最大塑性耗散功方面論證了非局部塑性模型符合熱力學(xué)第一和第二定律.然而實際應(yīng)用上并沒有達到最初提出這一模型的目的.這說明，基于非局部塑形模型矯正網(wǎng)格相關(guān)性的機理還未能得到深入理解，因而有深入研究的必要.

在文獻[29]中，作者以一般非局部塑性理論為基礎(chǔ)，并應(yīng)用狀態(tài)空間理論，通過局部和非局部兩個狀態(tài)空間的能量等效提出了一種求解應(yīng)變局部化問題的新方法.應(yīng)用這一方法可以得到與網(wǎng)格無關(guān)的結(jié)果.文獻[29]給出了理論框架并得到了一維問題的解析解.本文給出二維問題的一般列式和有限元解法，并給出一維問題的數(shù)值結(jié)果和二維問題的數(shù)值算例.

本文組織如下：第1節(jié)給出二維問題的屈服函數(shù)，軟化規(guī)律和流動法則在兩個狀態(tài)空間的一般列式；第2節(jié)導(dǎo)出非局部內(nèi)變量的求解方法；第3節(jié)給出有限元解法和應(yīng)力積分算法；第4節(jié)給出一維問題的數(shù)值結(jié)果并和解析解比較，并且給出二維剪切帶的數(shù)值算例；第5節(jié)得出結(jié)論.本文的工作基于小變形假設(shè).

1 兩個狀態(tài)空間的屈服函數(shù)和流動法則

在局部狀態(tài)空間中，如果域Ω中一點x的應(yīng)力狀態(tài)滿足屈服準則

則材料在x點開始屈服,即x∈Ωp，其中Ωp代表包括所有屈服材料點的區(qū)域；F(σ)為等效應(yīng)力，σY0表示初始屈服應(yīng)力，σ代表該點的應(yīng)力矢量，例如，對于二維問題，如果初始屈服后緊跟著應(yīng)變軟化，則假定在局部狀態(tài)空間中，塑性變形服從關(guān)聯(lián)的正交流動法則

(1)應(yīng)變軟化是一種滿足下列條件的全局行為

(2)屈服應(yīng)力σY為非局部內(nèi)變量?η在非局部狀態(tài)空間的特征塑性區(qū)域Ωcp上的積分的函數(shù).這里的Ωcp不同于局部狀態(tài)空間中的Ωp.

為了簡便且不失一般性，本文只考慮J2塑性的各向同性軟化，對于隨動軟化和各向異性軟化，只要遵循局部和非局部狀態(tài)空間塑性能量耗散率相等的原則，均可以導(dǎo)出類似的方程.非局部狀態(tài)空間的屈服條件和軟化規(guī)律通過下式引入

其中Vcp為特征塑性區(qū)域Ωcp的體積，并且對于二維區(qū)域，有

其中Acp代表二維特征塑性區(qū)域的面積，t為厚度.

應(yīng)該指出，Ωcp是非局部狀態(tài)空間的特征塑性區(qū)域，它的大小取決于材料內(nèi)尺度和權(quán)函數(shù)的具體形式.

其中，ˉEps是非局部軟化模量，ˉEps＜0.上標“ps”代表塑性軟化(plastic softening).對于非線性軟化，本文暫不涉及.由于本文后面的本構(gòu)關(guān)系方程均為增量形式，因此，文中的方法容易應(yīng)用于非線性軟化.因為非局部狀態(tài)空間的應(yīng)力狀態(tài)與局部狀態(tài)空間相應(yīng)的應(yīng)力狀態(tài)相同，方程(5)還可以寫為

塑性加載和卸載遵循Kuhn-Tucker互補條件[32]

以及一致性要求

將式(5)和式(9)代入到式(12)得到

其中

因為卸載過程中應(yīng)力是減小的，假設(shè)在某一時刻有一彈性試應(yīng)力σtr(參看文獻[32])，則˙σ可以表示為

˙σ=re˙σtr(15)

其中re是一個待確定的系數(shù).將式(15)和式(14)代入式(13)得到

相應(yīng)的彈性應(yīng)變率為

假設(shè)在非局部狀態(tài)空間中，關(guān)聯(lián)正交流動法則仍然適用，則有

將式(2)代入式(19)，對于vonM ises和Drucker-Prager屈服函數(shù)，則有

其中，γ為常數(shù).對于vonM ises屈服函數(shù)，容易證明

對于其他類型的屈服函數(shù)，γ可能并非常數(shù)，但本文的推導(dǎo)過程仍然適用.為方便起見，本文只考慮von M ises屈服準則，在后面的推導(dǎo)中取˙η=˙λ.類似地，在非局部狀態(tài)空間對于von M ises屈服函數(shù)可以證明

根據(jù)文獻[29]中式(7b)可知，在局部狀態(tài)空間的域Ωp內(nèi)的塑性能變化率為

其中正交流動法則式(2)已經(jīng)被引入到上式中.注意到只是流動法則而非軟化律被引入到了局部狀態(tài)空間.由于σ和˙λ(?f/?σ)在區(qū)域Ωp是連續(xù)的，對等式(23)右邊應(yīng)用積分中值定理可得

其中,xm為Ωp中的某點；(xm)代表(x)在Ωp中的平均值；Vp為Ωp的體積；對于二維情況，有Vp=tAp，其中t為材料的厚度，Ap為Ωp的面積.合并式(23)和式(24)可得

其中δ(Ap)為二維空間的Dirac-δ函數(shù)，定義如下

其中R2表示二維歐幾里得空間.因此式(25)可以寫為

其中

根據(jù)式(22b)和式(23)，將式(28)和式(29)代入文獻[29]中的式(19)可得

其中xp代表塑形區(qū)域Ωp中的某點.再將式(22b)和式(30)代入文獻[29]中的式(27)得

其中

從式(29)、式(31)以及式(32)可以看出，當(dāng)局部狀態(tài)空間的塑性乘子為Dirac-δ函數(shù)時，非局部狀態(tài)空間的非局部塑性乘子及其在特征塑性區(qū)域中的平均值具有明確定義.

2 非局部內(nèi)變量的求解

將式(17)和式(18)分別代入文獻[29]中的式(15a)和式(15b),然后再代入文獻[29]中的式(15c)和式(16)，并考慮到本文的式(30)，可以得到非局部狀態(tài)空間的總能量變化率為

令

則有

將式(34)和式(35)代入式(33b),得到

其中,αe和αp分別稱為彈性內(nèi)能比例和塑性內(nèi)能比例.由不等式(3)可知，對于應(yīng)變軟化，有αe＜0和αp＞1.將式(34)～式(36)代入式(29)，可得到〈〉m如下

3 有限元解法

對于一維問題，采用本文的模型和方法可以得出解析解.詳見文獻[29].對于二維及三維應(yīng)變局部化問題，本文所提出的模型和方法只能用有限元方法獲得數(shù)值解.由于有限元的應(yīng)變和位移場列式和平衡方程的建立與常規(guī)過程并無差異，故本文不作詳述，只列出與本構(gòu)方程和能量方程有關(guān)的算法.為了將注意力集中于由應(yīng)變軟化引起的局部化變形過程，假定載荷達到最大值后開始應(yīng)變軟化過程.為了保證整個計算過程中邊界條件的一致性，在彈性階段和塑性軟化階段，加載均由位移控制.彈性階段的響應(yīng)可通過常規(guī)的有限元方法獲得，在彈性分析中并不應(yīng)用本文所提出的方法.

3.1 局部狀態(tài)空間中屈服單元的監(jiān)測

在變形局部化的模擬中，通常假定某些邊界保持固定，而在另外某些邊界施加位移控制的載荷.假定在邊界?u0Ω上位移保持為零(其中u0代表零位移邊界)，在邊界?u1Ω施加位移向量?u?u1Ω(u1代表非零位移邊界，且?uΩ=?u1Ω∪?u0Ω).并且假定在當(dāng)前時間tn，材料體的狀態(tài)在局部空間完全確定.在下一個時間tn+1，即(n+1)級加載步，在邊界?u1Ω上施加位移增量

其中，i代表邊界N?u1Ω上的節(jié)點號，N?u1Ω代表邊界?u1Ω上的節(jié)點總數(shù)，?βn+1是第n+1級加載步的增量載荷參數(shù).在方向?u?u1Ω的合力增量為

當(dāng)滿足條件

則全局應(yīng)變軟化開始，變形局部化開始發(fā)生，具有高度塑性變形的剪切帶將域?Ω劃分為若干部分.令e代表單元編號，如果單元e的高斯積分點在局部狀態(tài)空間的應(yīng)力狀態(tài)滿足條件

則該單元成為屈服單元.其中，下標“GP”代表高斯積分點，q代表單元e中高斯積分點的最大編號.換句話說，只要單元e內(nèi)的任何一個高斯積分點滿足初始屈服條件，它就被視為屈服單元.與局部狀態(tài)空間相關(guān)聯(lián)的塑性功將集中于這些屈服的單元內(nèi)，局部內(nèi)變量將通過非局部權(quán)函數(shù)從這些屈服單元映射到非局部狀態(tài)空間，從而成為非局部內(nèi)變量計算的基礎(chǔ).

3.2 非局部狀態(tài)空間各分量的計算

當(dāng)材料經(jīng)歷應(yīng)變軟化時，假定在時刻tn，材料的狀態(tài)在局部狀態(tài)空間已經(jīng)由

確定，在非局部空間由

確定，并且滿足平衡條件.在時刻tn+1的目標是，當(dāng)在邊界?u1Ω上施加位移增量時，將狀態(tài)更新到tn+1時的狀態(tài)根據(jù)虛功原理，如果將上一時刻tn加載結(jié)束后對應(yīng)于邊界?u1Ω上的位移總量的反力作為已經(jīng)平衡的力(注意，邊界?gΩ上沒有力)，并且把在時間tn+1施加的位移增量作為虛位移，那么局部狀態(tài)空間中外力虛功的變化為

為了得到非局部狀態(tài)空間相應(yīng)的內(nèi)力虛功，首先需要得到在(n+1)加載步的彈性內(nèi)能和塑性內(nèi)能的比例和這可以通過類似于導(dǎo)出式(33)和式(36)的過程來實現(xiàn).

其中Vcp是與非局部狀態(tài)空間相關(guān)聯(lián)的特征塑形區(qū)域的體積，Ncp代表Vcp中單元的總數(shù).Vcp的具體值與權(quán)函數(shù)相關(guān)，將在后面的數(shù)值算例中討論.在tn+1時刻，非局部內(nèi)變量的全量值更新為

將式(45)代入式(18)中得到非局部塑性應(yīng)變增量為

在非局部狀態(tài)空間中，彈性內(nèi)能、塑性內(nèi)能以及總虛功的增量為

其中Nall代表Ω中的單元總數(shù).上式還可以寫為

其中

將式(49)～式(51)代入式(44)～式(48)，可以得到非局部狀態(tài)分量的增量并對全部非局部狀態(tài)和局部狀態(tài)的分量進行更新.系數(shù)re和彈性應(yīng)力增量?n+1可以依據(jù)式(16)和式(15)得到.式(16)和式(15)的離散形式可以分別寫為

其中試應(yīng)力增量?σtr可以通過對邊界?u1Ω施加位移增量??u1Ω，然后通過彈性分析得到.應(yīng)力增量?n+1和全應(yīng)力n需要代入非局部屈服條件式(4)以檢查其是否滿足屈服條件.如果條件

相應(yīng)的等效節(jié)點力為

其中，“A”代表有限元中的組合算子(詳見文獻[32]).單元剛度矩陣Ke的計算和及整體剛度矩陣K的形成可以通過常規(guī)方法得到

與塑性應(yīng)變相應(yīng)的等效位移可由下式得出

在這一過程中，自然邊界?u1Ω上沒有約束，因為邊界?u1Ω上的反力在塑性流動過程中保持不變，而邊界?u0Ω上的位移保持為零.

以上算法總結(jié)于如下：

假定在第n級載荷(應(yīng)變軟化階段)已得到邊界?u1Ω上的狀態(tài)分量

(1)在邊界?u1Ω上施加增量位移?n+1然后計算彈性試應(yīng)力?σtr

(2)計算外力功的變化?Wn+1

(5)計算非局部塑性乘子和塑性應(yīng)變增量

(6)計算系數(shù)re,彈性應(yīng)力和應(yīng)變增量?n+1和

(8)更新局部和非局部狀態(tài)空間的分量

4 算例

4.1 一維應(yīng)變局部化的數(shù)值解

應(yīng)用本文提出的數(shù)值解法對文獻[29]中的一維拉桿算例進行了有限元分析.拉桿幾何形狀和邊界條件參見文獻[29]中圖1.除了內(nèi)尺度以外的其他材料參數(shù)均和文獻[29]中的參數(shù)相同.此處材料內(nèi)尺度取l=6.26mm作為參照內(nèi)尺度，對應(yīng)于文獻[26]中Nilsson的模型中具有不同定義的材料內(nèi)尺度15.7mm.并以l=4.31mm，l=2.35mm作為對比，以考察不同的材料內(nèi)尺度對局部化區(qū)域尺寸的影響.有限元網(wǎng)格采用常應(yīng)變兩節(jié)點單元，每個單元有兩個高斯積分點.共采用的單元數(shù)量分別為5，11，41，101四種網(wǎng)格進行計算，數(shù)值計算結(jié)果示于圖1～圖5中.圖1和圖2表明，隨著網(wǎng)格越來越細小，應(yīng)變局部化過程中的塑性應(yīng)變分布、局部化區(qū)域的尺寸、以及載荷位移曲線都穩(wěn)定地收斂到解析解.圖3表明，隨著拉桿端部位移的增加，分布在局部化區(qū)域內(nèi)的塑性應(yīng)變也增加，局部化區(qū)域雖然有所擴展，但大體上不變.對于l=6.26mm的情況，局部化區(qū)域尺寸為Lcp=6l=6×6.26mm=37.56mm.圖4和圖5給出了塑性應(yīng)變分布和載荷--位移曲線隨著不同材料內(nèi)尺度l的變化情況.從這兩圖可以看出，隨著l的減小，塑性應(yīng)變越來越集中于更小的區(qū)域.當(dāng)l→0，本文的非局部方法得到的數(shù)值結(jié)果接近于局部塑性理論得到的解.這與文獻[29]得到的結(jié)論是一致的.

圖1 一維拉桿在不同有限元網(wǎng)格下的塑性應(yīng)變分布Fig.1 Plastic strain distributionsalong the axisof thebarunder tension for di ff erentmeshes

圖2一維拉桿在不同有限元網(wǎng)格下的載荷--位移曲線Fig.2 Load-displacementcurvesof thebarunder tensionw ith strain softening for di ff erentmeshes

圖3 一維拉桿中塑性應(yīng)變隨著端部位移增大的發(fā)展Fig.3 Developmentof plastic strainsalong theaxisof thebarw ith increased end displacement

圖4 一維拉桿內(nèi)塑性應(yīng)變隨著不同材料內(nèi)尺度的變化Fig.4 Plastic strain distributionsalong theaxisof thebar for di ff erent internal length scales

圖5 一維拉桿在不同的材料內(nèi)尺度下的載荷--位移曲線Fig.5 Load-displacementcurvesof thebarw ith di ff erentinternal length scales

4.2 二維剪切帶局部化的數(shù)值解

4.2.1 平面應(yīng)變加載試件

一個處于平面應(yīng)變加載狀態(tài)的試件的幾何和加載條件見圖6.試件尺寸為60mm×120mm，其底面放置在水平的剛性平面上，兩個側(cè)面受水平方向的剛性約束，頂面保持水平并通過一個剛性壓板緩慢地施加可控的位移.試驗過程中，垂直于XY平面的方向的變形保持為零.材料特性參數(shù)如下：E=11920MPa；泊松比v=0.2；σY0=100MPa；ps=-0.1E=-1192MPa.假定線性應(yīng)變軟化并且采用von M ises屈服準則.為了與一維算例一致，采用二維高斯型權(quán)函數(shù)如下

圖6 平面應(yīng)變加載試驗的試件和加載方式Fig.6 Geometry and loading condition of the specimen under plane strain test

為了觸發(fā)應(yīng)變局部化，在試件左下角指定尺寸7.5mm×7.5mm的方形區(qū)域為弱區(qū)域，它的屈服強度為并采用了3個不同的材料內(nèi)尺度l=6.26mm，l=3.13mm和l=1.565mm以考察剪切帶局部化受內(nèi)尺度影響的程度.有限元采用八節(jié)點等參元，形函數(shù)在單元之間是C0連續(xù).每個單元內(nèi)采用2×2高斯積分點.應(yīng)用3種不同的有限元網(wǎng)格，即128個單元(粗)，512個單元(中等)，和2048個單元(細).這些網(wǎng)格數(shù)量都是8×16網(wǎng)格的整數(shù)倍，每種網(wǎng)格的弱區(qū)域尺寸都相同.

為了確定應(yīng)變軟化過程中特征塑性區(qū)域Vcp的大小,在應(yīng)變局部化開始之后，首先在域Ωp中所有屈服單元的積分點的內(nèi)變量中選出最大值

然后，對每一屈服單元，如果滿足

那么它就會被包含到Vcp中，否則從Vcp除去.這將保證99.7%的塑性變形包含在Vcp中.

計算結(jié)果見圖7～圖10.對于內(nèi)尺度l=6.26mm，3種不同網(wǎng)格的變形情況示于圖7.變形的網(wǎng)格清晰地顯示了應(yīng)變軟化過程中變形集中在一個帶狀區(qū)域(剪切帶).剪切帶的寬度和方向基本上沒有隨著網(wǎng)格的不同而改變，這在圖8的等效塑性應(yīng)變分布中也可以看出.剪切帶的方向與X軸成近似45°.注意到，通過非局部塑性模型得到的剪切帶以外區(qū)域的變形與基于梯度塑性模型得到的結(jié)果略有不同(詳見文獻[18]).這是因為非局部模型的高斯型權(quán)函數(shù)是漸進趨于零的，而且在彈塑性邊界上沒有(不需要)指定邊界條件.圖9表明隨著網(wǎng)格的細化，載荷--位移曲線穩(wěn)定地收斂.采用精細網(wǎng)格得到的峰值后曲線比采用粗網(wǎng)格得到的要陡.

圖7 3種有限元網(wǎng)格的變形(l=6.26mm,?u=3mm)Fig.7 Deformation patterns for three di ff erentmeshes(l=6.26mm,?u=3mm)

圖8 3種有限元網(wǎng)格的等效塑性應(yīng)變分布Fig.8 Contour plotof theequivalentplastic strains for three di ff erentmeshes(l=6.26mm,=3mm)

圖10 給出了在邊界位移?u=1.15mm作用下，3種不同的內(nèi)尺度對應(yīng)的等效塑性應(yīng)變分布.從圖中可看出，剪切帶的寬度隨著內(nèi)尺度的減小而明顯減小，對應(yīng)于3個內(nèi)尺度l=6.26mm，3.13mm和1.565mm的3個剪切帶的寬度分別是大約38mm，19mm，9mm，大體上等于6l.這與文獻[29]中的式(58)預(yù)測的值是一致的.圖12表明隨著內(nèi)尺度的減小，峰值后下降段曲線的向下坡度增加.當(dāng)內(nèi)尺度變得很小時，如l=1.565mm，峰值后曲線就會變得非常陡.這些結(jié)果與一維應(yīng)變局部化的結(jié)果是一致的，并且也和梯度塑性模型預(yù)測的結(jié)果相符[16,18].

圖9 3種不同網(wǎng)格的載荷--位移曲線Fig.9 Load-displacementcurves for three di ff erentmeshes

圖10 不同材料內(nèi)尺度對應(yīng)的等效塑性應(yīng)變分布(?u=1.15mm)Fig.10 Contour plotof theequivalentplastic strains for three di ff erent internal length scale(?u=1.15mm)

圖11 不同材料內(nèi)尺度對應(yīng)的載荷--位移曲線Fig.11 Load-displacementcurves for di ff erentinternal length scales

圖12 邊坡的幾何尺寸和有限元網(wǎng)格Fig.12 Geometry and finit elementmesh of an embankment

4.2.2 邊坡穩(wěn)定性

為了進一步驗證本文提出的模型，考慮圖13所示的邊坡失穩(wěn)問題.邊坡上有一基礎(chǔ)，其剛度比邊坡材料大很多，可視為是剛性的.這一算例曾經(jīng)由Borja和Regueiro[33]用來驗證他們提出的強不連續(xù)局部化模型.本文采用下列材料特性參數(shù)：E=10MPa；v=0.4；σY0=100 kPa；ˉEps=-0.1E=-1.0MPa;容重γ=16kN/m2.仍采用von M ises屈服準則和高斯權(quán)函數(shù)以及線性軟化假定.有限元網(wǎng)格采用416個八節(jié)點等參元.為了使剪切帶的厚度大于網(wǎng)格尺寸，材料內(nèi)尺度取為l=150mm.計算得到邊坡失穩(wěn)后的網(wǎng)格變形和載荷--位移曲線分別如圖13和14所示.由圖13中可見，采用本文模型得到的剪切帶區(qū)域的網(wǎng)格變形是平滑的，而采用一般非局部塑性模型得到的剪切帶是不光滑的，且厚度大約和一個單元的尺寸相近.由圖14也可看出，用本文模型得到的載荷位移曲線在越過峰值點后平滑地逐漸下降，而采用一般非局部模型得到的曲線在越過峰值點后以很陡的斜率突然下降，后者是數(shù)值模擬結(jié)果具有網(wǎng)格相關(guān)性的典型特征.

圖13 邊坡失穩(wěn)后的變形Fig.13 Deformedmesh of the slope

圖14 邊坡失穩(wěn)過程中的載荷--位移曲線Fig.14 Load-displacementcurves for the slope

4.2.3 煤巖的平面應(yīng)變加載試驗

為了與試驗結(jié)果對比以驗證數(shù)值結(jié)果的可靠性，用本文模型對文獻[34]中的煤巖試件在平面應(yīng)變條件下的加載試驗進行了有限元分析.材料參數(shù)為：E=4.0GPa；v=0.19；黏聚系數(shù)c=16.83MPa;內(nèi)摩擦角φ=25°；ps=-400.0MPa.采用Drucker-Prager屈服準則和關(guān)聯(lián)流動法則，并仍用高斯權(quán)函數(shù)以及線性軟化的假定.由于前述算例已經(jīng)表明，采用本文提出的模型可得到與網(wǎng)格無關(guān)的結(jié)果，所以將試件劃分為300個單元以節(jié)約計算時間.材料內(nèi)尺假定為l=0.3mm.圖15是煤巖試件破壞形態(tài)的試驗結(jié)果和數(shù)值模擬結(jié)果的比較；圖16是載荷位移曲線的試驗和計算結(jié)果的比較.由于本試件在加載時沒有側(cè)向的約束載荷，因此試驗時沒有得到載荷--位移曲線的明顯下降段.但是由圖15和圖16可見，計算和試驗得到的剪切帶位置和傾角基本相同，試驗和計算得到的載荷--位移曲線的相當(dāng)接近的.由于試驗結(jié)果沒有給出剪切帶的厚度，所以無法推算出材料內(nèi)尺度，用于數(shù)值模擬所用的材料內(nèi)尺度是根據(jù)單元最小尺寸假定的.因此，計算得到的剪切帶厚度只能是大體的尺寸.值得提到的是，剪切帶厚度是試件斷裂之前塑性變形集中的區(qū)域，它的精確測量較為困難，尤其對于尺寸較小的試件[3435].

圖15 煤巖平面應(yīng)變加載破壞形態(tài)Fig.15 Failure pattern of the coalspecimen under plan strain compression test

圖16 煤巖平面應(yīng)變試驗的載荷--位移曲線Fig.16 Load-displacementcurves for the coalspecimen

5 結(jié)論

本文以非局部塑性理論和狀態(tài)空間理論為基礎(chǔ)，通過局部和非局部兩個狀態(tài)空間能量的等效原理，提出了一種對應(yīng)變局部化進行數(shù)值模擬的新方法.本文所提出的方法可以克服經(jīng)典的局部塑性理論和一般非局部塑形理論導(dǎo)致的網(wǎng)格相關(guān)性問題，從而得到應(yīng)變局部化問題的客觀解.本文導(dǎo)出了二維問題的一般列式，并給出了有限元算法.一維和二維問題的有限元算例表明，應(yīng)用本文的模型和算法得到的塑性應(yīng)變分布和載荷--位移曲線隨著單元網(wǎng)格尺寸的細化而穩(wěn)定地收斂.這說明，用本文的方法可得到與網(wǎng)格無關(guān)的結(jié)果.由于采用了兩個狀態(tài)空間塑性耗散能的等效，因而應(yīng)變局部化過程中能量的耗散總是大于零.算例分析還表明，剪切帶的厚度依賴于材料內(nèi)尺度，當(dāng)材料內(nèi)尺度趨近于零時，本文方法的結(jié)果將退化為經(jīng)典的局部塑性理論的結(jié)果.

與梯度塑性理論相比，在有限元分析中采用本文的模型和方法只需要位移場在單元之間具有C0連續(xù)性，不要求斜率的連續(xù)性，因此本文的算法容易嵌入到現(xiàn)有的有限元程序中而不需要做大的修改.對于大尺寸巖土結(jié)構(gòu)的有限元分析，本文的方法易于在有限元程序中實現(xiàn).

本文僅采用了von M ises和Drucker-Prager屈服準則求解各向同性應(yīng)變軟化導(dǎo)致的應(yīng)變局部化.將來進一步的研究需要考慮應(yīng)用較為復(fù)雜的屈服準則，如Mohr-Coulomb準則以及帽蓋模型(cap model)，并采用非關(guān)聯(lián)流動法則和非線性軟化模型.隨動軟化和各向異性軟化以及大變形情況下如何應(yīng)用本文的方法也有待研究.

致謝本文部分二維計算程序在中國科學(xué)院超級計算環(huán)境(ScGrid)上運行并獲得結(jié)果，特此致謝！

1黃茂松.土體穩(wěn)定與承載特性的分析方法.巖土工程學(xué)報,2016,38(1)：1-34(Huang Maosong.Analysismethods for stability and bearing capacity of soils.Chinese Journal ofGeotechnical Engineering,2016,38(1)：1-34(in Chinese))

2錢建固,呂璽琳,黃茂松.平面應(yīng)變狀態(tài)下土體的軟化特性與本構(gòu)模擬.巖土力學(xué),2009,30(3)：617-622(Qian Jiangu，L¨u Xilin，Huang Maosong.Softening characteristics of soilsand constitutive modeling under plane strain condition.Rock and SoilMechanics，2009，30(3)：617-622(in Chinese))

3徐連民,王天竹,祁德慶等.巖土中的剪切帶局部化問題研究：回顧與展望.力學(xué)季刊,2004,25(4)：484-489(Xu Lianmin，Wang Tianzhu，QiDeqing,etal.Study on geotechnicalshear band localization：retrospectand prospect.Chinese Quarterly ofMechanics,2004,25(4)：484-489(in Chinese))

4張啟輝,趙錫宏.粘性土局部化剪切帶變形的機理研究.巖土力學(xué),2002,23(1)：31-35(ZhangQihui，Zhao Xihong.Study onmechanism of localized shearband formation in clay.Rock and SoilMechanics,2002,23(1)：31-35(in Chinese))

5王杰,李世海,張青波.基于單元破裂的巖石裂紋擴展模擬方法.力學(xué)學(xué)報,2015,47(1)：105-118(Wang Jie,LiSihai,ZhangQingbo.Simulation of crack propagation of rock based on splitting elements.Chinese JournalofTheoreticaland Applied Mechanics,2015,47(1)：105-118(in Chinese))

6李宏儒,胡再強,馮飛等.結(jié)構(gòu)性黃土二元介質(zhì)本構(gòu)模型在局部化剪切帶中的應(yīng)用.巖土力學(xué),2012,33(9)：2803-2844(LiHongru，Hu Zaiqiang，F(xiàn)eng Fei,etal.Application ofstructural loessbinarymedium model to localization shearband.Rockand SoilMechanics,2012,33(9)：2803-2844(in Chinese))

7黃茂松,賈蒼琴,錢建固.巖土材料應(yīng)變局部化的有限元分析方法.計算力學(xué)學(xué)報,2007,24(4)：465-471(Huang Maosong，Jia Cangqin，Qian Jiangu.Strain localization problems in geomaterials using finit elements.Chinese Journal ofComputational Mechanics，2007，24(4)：465-471(in Chinese))

8邵龍?zhí)?劉港,郭曉霞.三軸試樣破壞后應(yīng)變局部化影響的實驗研究.巖土工程學(xué)報,2016,38(3)：385-394(Shao Longtan，Liu Gang，Guo Xiaoxia.E ff ects of strain localization of triaxial samples in post-failure state.Chinese Journal ofRock Mechanics and Engineering,2016,38(3)：385-394(in Chinese))

9李學(xué)豐,黃茂松,錢建固.基于非共軸理論的各向異性砂土應(yīng)變局部化分析.工程力學(xué),2014,31(3)：205-246(LiXuefeng，Huang Maosong，Qian Jiangu.Strain localization analysis of anisotropic sands based on non-coaxial theory.Engineering Mechanics,2014,31(3)：205-246(in Chinese))

10王水林,鄭宏,劉泉聲等.應(yīng)變軟化巖體分析原理及其應(yīng)用,巖土力學(xué),2014,35(3)：609-630(Wang Shuilin，Zheng Hong，Liu Quansheng,etal.Principleofanalysisof strain-softening rockmass and itsapplication.Rock and SoilMechanics,2014,35(3)：609-630(in Chinese))

11 Read HE,Hegemier GA.Strain softening of rock,soil and concrete—a review article.Mechanics and Materials,1984,3(3)：271-294

12 Vardoulakis I.Shear band inclination and shearmodulus of sand in biaxial test.International Journal for Numerical and Analytical Methods in Geomechanics,1980,4(2)：103-119

13 Jir′asek M,Ba?zant ZP.Inelastic analysis of structures.Chichester,England：JohnWiley&Sons,2002

14 Pietruszczak ST,M r′oz Z.Finite element analysis of deformation of strain-softening materials.International Journal for Numerical Methods in Engineering,1981,17(3)：327-334

15 OrtizM,Leroy Y,Needleman A.A finit elementmethod for localized failure analysis.ComputerMethods in Applied Mechanicsand Engineering,1987,61(2)：189-214

16 de Borst R,Sluys LJ,M¨uhlhaus H-B,et al.Fundamental issues in finit elementanalyses of localization of deformation.Engineering Computations,1993,10(2)：99-121

17 Simo JC.Strain softening and dissipation：a unificatio of approaches//Mazars J,Ba?zant JP,eds.Cracking and Damage：strain localization and size e ff ect.London and New York：Elsevier Applied Science,1988：440-461

18 de Borst R,M¨uhlhaus H-B.Gradient-dependent plasticity：formulation and algorithm aspects.International Journal for Numerical Methods in Engineering,1992,35(3)：521-539

19 Luzio GD，Ba?zant ZP.Spectral analysis of localization in nonlocaland over-nonlocalmaterialsw ith softening plasticity ordamage.International JournalofSolidsand Structures,2005,42(23)：6071-6100

20 Vermeer PA,BrinkgreveRBJ.A new e ff ectivenon-localstrainmeasure for softening plasticity//Chambon V,Desrues J,Vardoulakis I,eds.Localisation and Bifurcation Theory for Soilsand Rocks.Rotterdam,Netherland：A.A.Balkema,1994

21 Str¨omberg L,Ristinmaa M.FE-formulation of a nonlocalplasticity theory.ComputerMethods in Applied Mechanicsand Engineering,1996,136(1-2)：127-144

22 L¨u XL,Bardet J-P，Huang MS.Numerical solutions of strain localization w ith nonlocal softening plasticity.ComputerMethods in Applied Mechanicsand Engineering,2009,198(47)：3702-3711

23 L¨u XL,Bardet J-P，Huang MS.Spectralanalysisof nonlocal regularization in two-dimensional finit elementmodels.InternationalJournal of Numerical and Analytical Methods in Geomechanics,2012,36(2)：219-235

24 Ba?zant ZP,Lin F-B.Nonlocalsmeared crackingmodel for concrete fracture.Journal ofstructural Engineering,1988,114(11)：2493-2510

25 Eringen AC.Nonlocal Continuum Field Theories.New York：Springer-Verlag,2002

26 Nilsson C.Nonlocal strain softening bar revisited.International JournalofSolidsand Structures,1997,34(34)：4399-4419

27 Borino G,Fuschi P,Polizzotto C.A thermodynamic approach to nonlocalplasticity and related variationalprinciples.JournalofApplied Mechanics,1999,66(4)：952-963

28 de Sciarra FM.A general theory for nonlocal softening plasticity of integral-type.International Journal of Plasticity,2008,24(8)：1411-1439

29武守信,魏吉瑞,楊舒蔚.基于能量等效原理的應(yīng)變局部化分析I：一維解析解.力學(xué)學(xué)報,2017,49(3)：(Wu Shouxin,WeIJirui,Yang Shuwei.Analysisof strain localization by energy equivalence between local and nonlocal state spaces-part I：one-dimensional analyticalsolution.Chinese JournalofTheoreticaland Applied Mechanics,2017,49(3)：(in Chinese))

30馮春,李世海,劉曉宇.基于顆粒離散元法的連接鍵應(yīng)變軟化模型及其應(yīng)用.力學(xué)學(xué)報,2016,48(1)：76-85(Feng Chun,LiShihai,Liu Xiaoyu.Particle-dem based linked bar strain softening model and itsapplication.Chinese JournalofTheoreticaland Applied Mechanics,2016,48(1)：76-85(in Chinese))

31萬征,姚仰平,孟達.復(fù)雜加載下混凝土的彈塑性本構(gòu)模型.力學(xué)學(xué)報,2016,48(5)：1159-1171(Wan Zheng,Yao Yangping,Meng Da.An elastoplastic constitutivemodel of concrete under complicated load.Chinese JournalofTheoreticaland Applied Mechanics,2016,48(5)：1159-1171(in Chinese))

32 Simo JC,Hughes TJR.Computational Inelasticity.New York：Springer-Verlag,1998

33 Regueiro RA,Borja RI.A finit elementmodel of localized deformation in frictionalmaterials takingastrong discontinuity approach.Finite Elements in Analysisand Design,1999,33：283-315

34 Yum lu M,Ozbay MU.A study of the behavior of brittle rocks under plan strain and triaxial loading condition.International Journal ofRockMechanicsand M ining Science&Geomechanics Abstracts,1995,32(7)：725-733

35韓鐵林,師俊平,陳蘊生等.軸、側(cè)向同卸荷下砂巖力學(xué)特性影響的試驗研究.力學(xué)學(xué)報,2016,48(3)：936-943(Han Tielin，Shi Junping，Chen Yunsheng，etal.Experimental study onmechanics characteristicsofsandstoneunderaxialunloadingand radialunloading path.Chinese Journal of Theoretical and Applied Mechanics,2016,48(3)：936-943(in Chinese))

ANALYSISOFSTRAIN LOCALIZATION BY ENERGY EQUIVALENCE:II.FINITE ELEMENT SOLUTION1)

Wu Shouxin?,?,2)Wei Jirui?,?Yang Shuwei?,??

(DepartmentofBridge Engineering，SchoolofCivil Engineering，Southwest Jiaotong University，Chengdu 610031，China)?(Key Laboratory ofTransportation TunnelEngineering ofMinistry ofEducation，Chengdu 610031，China)

Founded on the nonlocalplasticity and the state space theories,anew approach isproposed to fin themeshindependent solution of the strain localization problems by equating the rates of plastic energy dissipation in the local and nonlocal state spaces.Follow ing the previous paper by the authors,general formulas are developed for the solution of the nonlocal internal variables in the two-andmore than two-dimensional problems.A stress updating algorithm is proposed to integrate the rate form constitutive equations in the finit elementcontext.To verify the proposed approach,a one-dimensionalmodel problem and three two-dimensional plane strain problems are solved numerically by the finit elementmethod.Numerical results show that the plastic strain distributions and the load-displacement curves stablyconvergewithrefinemen of the finit elementmesh.The sizeof the localization zone dependsonly on the internal length scaleand is insensitive to themesh size.For the one-dimensionalproblem,numericalsolutions converge to the analytical ones.For the two-dimensionalproblems,although no analyticalsolutionsareavailable,thenumericalsolutions converge toward the unique ones.Thew idth and the inclination are almost not changed as themesh size is reduced.Also,the distribution of the plastic strainsand the deformation patternsare smooth in the entire domain.A slope stability problem and a plane strain test of a coal specimen are also solved numerically to demonstrate the robustness of the proposed approach.It iswell shown that the proposed approach can overcome the drawbacks of the classical continuum theory and lead to physicallymeaningful,mesh-independentsolution of strain softening problems.Because only C0continuity isneeded between elementboundaries,the proposed approach is easy to be incorporated into the existing finit element codew ithoutsubstantialmodification

strain localization,nonlocalplasticity,internal length scale,mesh-dependence,finit element

O344.3，TU501

A

10.6052/0459-1879-16-330

2016-11-14收稿，2017-03-20錄用,2017-03-21網(wǎng)絡(luò)版發(fā)表.

1)教育部留學(xué)回國人員科研啟動基金(201250300)和西南交通大學(xué)土木工程學(xué)院基礎(chǔ)研究創(chuàng)新計劃基金資助項目.

2)武守信，副教授，博士，主要研究方向：橋梁和巖土結(jié)構(gòu)的有限元分析和本構(gòu)關(guān)系.E-mail：swu@home.sw jtu.edu.cn

武守信，魏吉瑞，楊舒蔚.基于能量等效原理的應(yīng)變局部化分析：II.有限元解法.力學(xué)學(xué)報,2017,49(4)：880-893 Wu Shouxin,WeiJirui,Yang Shuwei.Analysisofstrain localization by energy equivalence：II.Finiteelementsolution.Chinese Journal ofTheoreticaland Applied Mechanics,2017,49(4)：880-893