張 宇， 盧廣達(dá)， 夏曉舟

（1.河海大學(xué)力學(xué)與材料學(xué)院工程力學(xué)系，南京 211100； 2.同濟(jì)大學(xué)土木建筑工程學(xué)院，上海 200092）

邊坡滑移破壞的數(shù)值模擬一直是工程界和學(xué)術(shù)界關(guān)注的焦點(diǎn). 除有限元之外，發(fā)展了很多專門模擬破壞的數(shù)值方法，如界面元［1］、塊體元［2］、流形元［2］、光滑粒子流體動(dòng)力學(xué)（SPH）［3］、粒子有限元［4］、光滑粒子有限元［5］、物質(zhì)點(diǎn)法［6］、無網(wǎng)格法［7］、擴(kuò)展有限元法［8］以及近場動(dòng)力學(xué)［9-10］. 離散元和塊體元需提前定義滑移面，如果滑移面事先未知，則需要不斷搜索，通過比較最速能量釋放來確定最危險(xiǎn)滑移面，在搜索過程中，網(wǎng)格需不斷重構(gòu)，非常繁瑣. 流形元?jiǎng)t是通過定義在覆蓋上的函數(shù)進(jìn)行加權(quán)平均來定義區(qū)域的位移模式，因而形式上比較靈活，能夠構(gòu)造出非連續(xù)的位移模式. 光滑粒子流體動(dòng)力學(xué)（SPH）方法是針對(duì)每個(gè)粒子進(jìn)行牛頓力學(xué)分析，以追蹤粒子的運(yùn)動(dòng)軌跡. 粒子無網(wǎng)格法則是在局部范圍引入徑向基的核函數(shù)，因而能夠反映應(yīng)變局部化的破壞現(xiàn)象，但計(jì)算量較大，且本征邊界條件不能直接施加. 擴(kuò)展有限元是在傳統(tǒng)的位移模式中，利用單位分解的思想，疊加一個(gè)反映裂紋處位移間斷（或界面處應(yīng)變間斷）的附加不連續(xù)函數(shù)［11］以及裂尖處應(yīng)力奇異的附加三角函數(shù)［12］，所以無須重構(gòu)網(wǎng)格就能實(shí)現(xiàn)裂紋的追蹤模擬，但裂紋的擴(kuò)展機(jī)制受限于斷裂力學(xué)理論，這對(duì)于裂紋分叉和動(dòng)態(tài)裂紋發(fā)展仍然是一個(gè)挑戰(zhàn)，且裂紋在發(fā)展過程中，盡管網(wǎng)格不需要重構(gòu)，但附加自由度則需要隨時(shí)進(jìn)行更新以匹配當(dāng)下的拓?fù)錁?gòu)型. 近場動(dòng)力學(xué)某種意義上是一種超越了尺度概念的分子動(dòng)力學(xué)模型，數(shù)值方法上類似于光滑粒子動(dòng)力學(xué)，但物理內(nèi)涵更為豐富，具體體現(xiàn)在結(jié)構(gòu)上能夠從連續(xù)狀態(tài)自然地過渡到非連續(xù)狀態(tài). 在近場動(dòng)力學(xué)模型中，每個(gè)點(diǎn)的力狀態(tài)由其鄰近點(diǎn)與該點(diǎn)的相對(duì)變形所決定，當(dāng)該點(diǎn)與鄰近所有的點(diǎn)發(fā)生關(guān)聯(lián)，則表示連續(xù)，當(dāng)該點(diǎn)與鄰近所有點(diǎn)切斷了聯(lián)系，則表示在該點(diǎn)處不連續(xù)，處在這兩種狀態(tài)之間的為過渡態(tài)，即損傷介于0到1之間，因而能夠突破連續(xù)介質(zhì)力學(xué)中發(fā)生奇異或間斷的微分運(yùn)算，如涉及的非局部變形梯度. 因而近場動(dòng)力學(xué)模型能夠追蹤裂紋的發(fā)展而無須更改初始構(gòu)型的拓?fù)浣Y(jié)構(gòu)；盡管近場動(dòng)力學(xué)模型是基于積分型的表達(dá)，某種程度上抹掉了某些奇異和間斷的劇烈差異，但仍然保持著物理上的變形演變勢特征，以驅(qū)動(dòng)裂紋的發(fā)生. 實(shí)際中，理論上的應(yīng)力奇異是不存在的，裂紋處的應(yīng)變也是不存在的，因而非局部變形梯度的引入大大拓展了變形場的適用范圍，使得拓?fù)浣Y(jié)構(gòu)不斷變化的問題（如裂紋發(fā)展）能夠在統(tǒng)一的框架下（不增加自由度）進(jìn)行. 當(dāng)然，近場動(dòng)力學(xué)模型的主要問題是計(jì)算量大，存在邊界效應(yīng)，積分項(xiàng)中需要引入高階微分算子才能達(dá)到精度上的要求，否則存在一定的數(shù)值振蕩［13］，而與有限元混合建模是一個(gè)很不錯(cuò)的出路.

本構(gòu)方面，主要分為局部型本構(gòu)模型和非局部型的本構(gòu)模型兩類. 傳統(tǒng)的塑性模型、損傷模型和彈塑性損傷模型都屬于局部型本構(gòu)模型，材料參數(shù)都是基于一定尺寸的試樣測試的，因而參數(shù)和賦予的應(yīng)力-應(yīng)變關(guān)系曲線都是基于試樣平均化獲得的. 因此，在處理材料屈服后的力學(xué)行為時(shí)，局部化模型難以準(zhǔn)確地反映應(yīng)力釋放和轉(zhuǎn)移機(jī)制，也就很難捕捉應(yīng)變局部化的效果. 應(yīng)變梯度理論［14］（包括偶極理論，如Cosserat理論［15］）中物質(zhì)點(diǎn)具有一定的尺寸大?。ㄍǔＨ〔牧系膬?nèi)稟尺度），能夠反映物質(zhì)點(diǎn)的微彎曲和微扭轉(zhuǎn)機(jī)制，屬于非局部化模型. 應(yīng)變梯度理論模型對(duì)于裂尖場的變形能夠精細(xì)地反映，但該理論模型最大的問題是物質(zhì)點(diǎn)的內(nèi)稟尺度不好確定，且在廣義應(yīng)力空間（包括偶應(yīng)力）中構(gòu)建屈服面時(shí)很難通過實(shí)驗(yàn)進(jìn)行標(biāo)定. 另一種非局部模型，就是通過引入非局部積分算子，對(duì)要關(guān)注的物理量（如塑性應(yīng)變、損傷）進(jìn)行加權(quán)平均化處理［16-17］，這種積分型的非局部化模型不僅能夠克服網(wǎng)格依賴性問題，且某種程度上能夠反映破壞過程中的應(yīng)變局部化現(xiàn)象.

本文借鑒近場動(dòng)力學(xué)中點(diǎn)對(duì)的思想，把微損傷定義在點(diǎn)對(duì)上，并在考察點(diǎn)的近場范圍內(nèi)通過加權(quán)平均化處理對(duì)該點(diǎn)的整體損傷進(jìn)行評(píng)估，最后引入損傷導(dǎo)致的剛度退化模式，構(gòu)建一個(gè)非局部化的損傷本構(gòu)模型，并嵌入到有限元的計(jì)算框架下，實(shí)現(xiàn)邊坡的滑移破壞模擬.

文章分為5個(gè)部分，在接下來的第2部分詳細(xì)介紹非局部損傷本構(gòu)模型的構(gòu)建過程，包括鍵損傷模型、點(diǎn)的整體損傷評(píng)估和損傷導(dǎo)致的剛度退化模式；第3部分將展示多尺度有限元建模理念及通過四點(diǎn)剪切梁進(jìn)行剪切荷載導(dǎo)致的破壞現(xiàn)象進(jìn)行驗(yàn)證，以便應(yīng)用到邊坡的滑移破壞模擬，即第4部分的內(nèi)容；第5部分則針對(duì)本文構(gòu)建的模型進(jìn)行深入的討論和總結(jié).

1 非局部損傷本構(gòu)建模

模型任意一點(diǎn)的受力狀態(tài)應(yīng)該由其周邊點(diǎn)與該點(diǎn)的相對(duì)位移所決定. 另外，實(shí)驗(yàn)給定的應(yīng)力-應(yīng)變關(guān)系曲線也是基于試樣均勻化獲得的，因此當(dāng)試樣內(nèi)部出現(xiàn)微裂紋時(shí)，所發(fā)生應(yīng)力釋放和轉(zhuǎn)移在局部模型中是反映不出來的，而這正是應(yīng)變局部化的根源，所以在傳統(tǒng)的損傷本構(gòu)模型中呈現(xiàn)的損傷分布不夠集中. 為此，本文受近場動(dòng)力學(xué)非局部作用思想的啟發(fā)，構(gòu)建一個(gè)基于鄰域點(diǎn)對(duì)（鍵）變形的非局部損傷本構(gòu)模型.

1.1 鍵損傷模型

對(duì)于任一點(diǎn)x，考察其周邊點(diǎn)x′與點(diǎn)x 形成的鍵的變形損傷情況. 定義點(diǎn)對(duì)x′x 的正伸長率ρ+為

圖1 物質(zhì)點(diǎn)對(duì)的相對(duì)運(yùn)動(dòng)Fig.1 The relative motion of a pair of points

在鍵的損傷演化方程中，損傷的發(fā)展速率受歷史最大超越伸長率及控制參數(shù)γ影響，γ值取得越大，損傷發(fā)展就越快，圖2 給出了不同γ 所對(duì)應(yīng)的損傷隨歷史最大超越伸長率的演化曲線.

1.2 點(diǎn)的損傷狀態(tài)評(píng)估

由點(diǎn)影響域內(nèi)鍵損傷情況，對(duì)該點(diǎn)的整體損傷Ω 進(jìn)行加權(quán)評(píng)估如下：

圖2 鍵損傷演化曲線Fig.2 Damage evolution curves of bond

式中：V′為x點(diǎn)所占據(jù)的空間大小；φ( x′,x )為鍵損傷的權(quán)函數(shù)，且要求∫H(x)φ( x′,x )dV′=1.0. 由正伸長率的定義式（1）可知，對(duì)于均質(zhì)材料，影響域內(nèi)距離點(diǎn)x 最遠(yuǎn)的點(diǎn)x′形成的點(diǎn)對(duì)x′x 鍵更容易斷裂，即影響域內(nèi)遠(yuǎn)處點(diǎn)較近處點(diǎn)優(yōu)先脫離影響范圍，這與近場動(dòng)力學(xué)的鍵伸長率定義不一樣.

1.3 剛度退化模式

由胡克定律得：

式中：De為初始彈性矩陣；Ded為彈性損傷矩陣；εE為彈損傷應(yīng)變；εe為彈性應(yīng)變；εed為微裂紋閉合的那部分應(yīng)變. 如果把介質(zhì)區(qū)分為損傷部分和未損傷部分，則對(duì)于未損傷部分，有

由式（5）、（7）～（9），導(dǎo)出彈損傷矩陣與彈性矩陣之間的關(guān)系為

這就是二次剛度退化的模式，但實(shí)際情況是，隨著損傷的發(fā)展，材料點(diǎn)的影響域范圍在縮小，即影響半徑l 隨著影響域內(nèi)部鍵的斷裂而減小，可根據(jù)有效面積的減小定義影響半徑隨損失的演化關(guān)系為

式中：l0為點(diǎn)x 的初始影響域半徑. 由公式（1）可知，影響半徑的減小導(dǎo)致伸長率的增大，從而由鍵損傷判據(jù)（3）式和宏觀損傷（4）式可獲得這樣的結(jié)論：隨著損傷的發(fā)展，鍵更加容易斷裂，即材料的抗力隨著損傷的發(fā)展而逐步衰減，從而加速了剛度的退化，即宏觀損傷Ω 與鍵損傷ω 是相互促進(jìn)的.

2 模型驗(yàn)證

2.1 多尺度有限元建模

在大部分模型中，出現(xiàn)破壞的地方總是集中在很狹小的局部帶中，因此只需在容易出現(xiàn)破壞的范圍內(nèi)劃分足夠密的網(wǎng)格，且采用提出的非局部損傷本構(gòu)模型；其他部分的網(wǎng)格剖分只要保持工程精度即可，且采用彈性本構(gòu)模型. 這樣既能保持足夠的精度，又能提高計(jì)算效率.

對(duì)于非局部損傷本構(gòu)模型，為確保點(diǎn)的整體損傷評(píng)估精度，在給定近場范圍，點(diǎn)對(duì)數(shù)量越多越好，至少保證每個(gè)點(diǎn)的近場范圍內(nèi)有40 個(gè)點(diǎn)對(duì)以上.另外，由于單元?jiǎng)偠染仃囀腔诜e分點(diǎn)進(jìn)行加權(quán)求和的，因此，本文的點(diǎn)對(duì)直接通過單元積分點(diǎn)進(jìn)行構(gòu)建，為方便示意，以三角形單元為例，且每個(gè)單元布置一個(gè)積分點(diǎn)，如圖3所示.

2.2 四點(diǎn)剪切梁

下面以四點(diǎn)剪切梁為例說明本文提出的非局部損傷本構(gòu)模型是否適用剪切型的破壞模式. 梁模型、約束和加載如圖4所示，梁的長度L=0.8 m，梁高H=0.2 m，開口裂縫長a=0.04 m，支座間距b=0.08 m.估算損傷破裂范圍，并進(jìn)行多尺度網(wǎng)格剖分，如圖5所示. 對(duì)細(xì)網(wǎng)格單元，應(yīng)用構(gòu)建的非局部損傷本構(gòu)模型，其中材料彈性模量E=40 GPa，泊松比ν=0.2，臨界伸長率ρc=1.5e-4，控制參數(shù)γ=20.采用弧長法迭代，直到梁結(jié)構(gòu)不能繼續(xù)承載時(shí)結(jié)束計(jì)算，獲得的損傷發(fā)展如圖6所示，從圖中可以看出，損傷從開口裂縫處發(fā)展，發(fā)展路徑和試驗(yàn)破壞形態(tài)（圖7）是一致的，且損傷分布比較集中，說明本文提出的非局部損傷本構(gòu)模型是合理可靠的，能夠模擬出剪切荷載的破壞形態(tài).

圖3 有限元空間離散以及物質(zhì)點(diǎn)對(duì)示意Fig.3 Discretization of finite element space and point pairs of matter

圖4 四點(diǎn)剪切梁Fig.4 Four-point shear beam

圖5 梁的多尺度網(wǎng)格剖分Fig.5 Multi-scale meshing of beams

圖6 四點(diǎn)剪切梁的損傷發(fā)展Fig.6 Damage development of four-point shear beam

圖7 試驗(yàn)破壞后的形態(tài)Fig.7 Shape after test failure

3 邊坡失穩(wěn)模擬

某土質(zhì)邊坡，幾何尺寸如圖8所示，多尺度網(wǎng)格剖分如圖9所示，彈性模量E=1 000.0 MPa，泊松比ν=0.45，臨界伸長率ρc=1e-4，控制參數(shù)γ=10 . 文獻(xiàn)［18］采用通過輔助剛塊施加壓力荷載P，本文采用位移加載模式，即對(duì)剛塊接觸面的節(jié)點(diǎn)施加相同豎向位移荷載，計(jì)算后獲得的損傷分布如圖10 所示，從圖中可以看出一條明顯的損傷滑移帶，呈弧形狀，與文獻(xiàn)結(jié)果相同. 說明本文提出的方法完全適用于剪切滑移破壞問題，而且實(shí)現(xiàn)起來簡單、高效.

圖8 某土質(zhì)邊坡約束及受載Fig.8 Restraint and load of a soil slope

圖9 土質(zhì)邊坡的多尺度網(wǎng)格剖分Fig.9 Multi-scale meshing of soil slope

圖10 邊坡的滑坡模擬Fig.10 Slope landslide simulation

4 討論與結(jié)語

本文受近場動(dòng)力學(xué)的思想啟發(fā)，構(gòu)建了一個(gè)非局部的損傷本構(gòu)模型，該模型突破了連續(xù)介質(zhì)力學(xué)中應(yīng)力-應(yīng)變關(guān)系均勻化標(biāo)定的限制，這種均勻化標(biāo)定對(duì)于材料出現(xiàn)破裂后是不太適合的，因?yàn)樯婕皯?yīng)力釋放和轉(zhuǎn)移，這就是為什么傳統(tǒng)損傷本構(gòu)模型存在網(wǎng)格依賴性問題的原因. Bazant曾通過非局部化算子對(duì)損傷進(jìn)行加權(quán)平均處理［17］，使得任意點(diǎn)的力學(xué)狀態(tài)同其周邊點(diǎn)聯(lián)系起來，這樣試樣的均勻化應(yīng)力-應(yīng)變曲線賦予了該點(diǎn)的非局部范圍，而不是該點(diǎn)自身，但在Bazant模型中容易出現(xiàn)應(yīng)力鎖死，即殘余應(yīng)力由于局部平均化的處理，即便損傷達(dá)到1.0時(shí)也消除不了. 本文提出的非局部損傷模型，其損傷本質(zhì)上是不連續(xù)性程度的表征，因而損傷分布代表的就是缺陷的發(fā)展，而在傳統(tǒng)唯象損傷模型［19］中，所謂損傷的分布實(shí)際上是退化剛度分布. 因而基于非局部損傷的嚴(yán)密定義，結(jié)合宏觀損傷（非局部）與鍵損傷相互促進(jìn)的影響機(jī)制，建立了隨損傷發(fā)展的剛度退化模型和基于影響域范圍逐步縮小的抗力衰減模式，很好地反映了結(jié)構(gòu)在破壞過程中的應(yīng)力釋放和轉(zhuǎn)移機(jī)制，使得一旦某個(gè)單元的損傷發(fā)生，其對(duì)節(jié)點(diǎn)力的貢獻(xiàn)瞬間減弱很多，導(dǎo)致節(jié)點(diǎn)失衡速率加快，向損傷單元反向運(yùn)動(dòng)，形成對(duì)損傷單元的進(jìn)一步拉扯或錯(cuò)動(dòng)，達(dá)到應(yīng)變局部化的效果.

相比于近場動(dòng)力學(xué)模型，本文保留了連續(xù)介質(zhì)力學(xué)框架的一些特性，比如彈性矩陣（包括彈損傷矩陣）天然地反映點(diǎn)對(duì)之間的相互作用，這是鍵基近場動(dòng)力學(xué)模型所忽視的，盡管后面的近場動(dòng)力學(xué)通過各種辦法突破了泊松比的某些限制，但對(duì)于高泊松比材料如本文算例中的土質(zhì)邊坡，仍然不適用. 當(dāng)然，如果采用態(tài)基的近場動(dòng)力學(xué)模型，把連續(xù)介質(zhì)力學(xué)的材料非線性本構(gòu)模型轉(zhuǎn)化到近場動(dòng)力學(xué)中，也能達(dá)到邊坡的滑移破壞效果［20］，但模型復(fù)雜，且不能如本文方法在統(tǒng)一的數(shù)值方法中（如有限元）實(shí)現(xiàn)多尺度模擬，達(dá)到高效計(jì)算的目的.

本文方法簡單高效，容易拓展到多物理場的耦合計(jì)算；對(duì)于黏彈性損傷和蠕變問題，也容易拓展，只需對(duì)鍵損傷模型進(jìn)行改造即可，因而本文方法具有重要的理論意義.