基于離散元的桿件斷裂全過程數值模擬

齊 念， 鄭開啟

(1.南京交通職業技術學院，南京 211188；2.南京林業大學 土木工程學院，南京 210037)

1 引 言

結構在服役期內會遭遇多種荷載和作用，尤其是在強震、強風以及爆炸等極端荷載作用下，組成結構的部分桿件會進入塑性，甚至會出現損傷和裂紋；隨著荷載作用的不斷加強或持續，微細裂紋不斷擴展，當裂紋貫通桿件全截面時，構件發生斷裂，最終引起結構部分或完全倒塌[1]。整個過程中，桿件經歷了彈性小變形、塑性開展、損傷累積直至最后斷裂。采用數值方法對破壞全過程進行模擬，需實時模擬結構的幾何大變形、材料非線性以及構件斷裂失效不連續等問題，這對現有的數值計算方法是一個巨大的挑戰[2,3]。如Lynn等[4]通過對有限元方法進行修正，將自適應變換積分高斯技術引入到有限元方法中，處理單元的斷裂和接觸問題，對飛機沖撞世貿大樓倒塌破壞問題進行了模擬分析；Zhang等[5]將粘結力模型和拓撲數據結構優化引入有限元，用以模擬混凝土結構中的動力破壞和裂紋擴展問題。但有限元法模擬裂紋擴展問題時需不斷進行網格重劃分，使得計算非常復雜并且效率較低。

與有限元法的理論基礎不同，離散元法DEM屬于非連續介質數值方法，現已成為研究巖石巖土、顆粒散粒體以及混凝土等材料非線性力學行為的常用方法[6]。由于不要求滿足位移連續和變形協調條件，因此DEM能方便地應用于各類材料非連續、非均勻以及結構大變形和斷裂破壞等復雜過程及其失效機理的研究。如Potyondy[7]基于應力腐蝕理論提出了平行黏結應力腐蝕模型，以模擬應力或水溶液等對巖石膠結的應力腐蝕作用；楊升等[8]基于PFC3D顆粒流程序，對砂土在剪切過程中的體積變化及力學行為進行了宏細觀分析；劉志林等[9]利用顆粒流離散元方法建立了一種含骨料、砂漿和過渡層的混凝土細觀力學模型，對彈丸侵徹混凝土進行了研究。為了拓展離散元法的應用范圍，在前期研究工作中，本文以顆粒離散元為基礎，提出了一種適于桿系結構分析的DEM計算模型，并成功地應用于平面框架結構、大跨空間網格結構的彈性、動靜力響應分析和幾何大變形分析[10-12]。與有限元法不同，構建框架結構的DEM模型時，是將梁柱等桿件離散為剛性單元集合體，單元與單元之間采用特定的接觸彈簧連接，每個單元的運動遵循牛頓第二定律。從方法的本身屬性而言，DEM處理斷裂等不連續問題在行為模式上更加直觀，因此較有限元法更具優勢。

本文基于已有研究，提出了一種適于結構彈塑性分析的DEM纖維梁模型，然后構建了桿件斷裂模擬算法，并自編了構件斷裂計算程序；最后將其應用于懸臂梁結構和單層網殼振動臺倒塌試驗兩個算例，對桿件的彈性和彈塑性直至斷裂全過程進行模擬，通過數值分析驗證該方法的正確性和適用性。

2 DEM纖維梁模型

在桿系結構的彈性DEM方法中，顆粒球元之間的接觸粘結看作是一根彈性梁，兩球心的距離即為彈性梁的長度L，且采用構件實際截面的幾何物理參數計算各彈簧接觸剛度系數[12]。彈塑性分析時，借鑒纖維模型的研究思路，將兩相鄰顆粒球元之間的接觸粘結用分布式彈簧進行等效，分布彈簧看作是梁截面的若干根纖維，每一根彈簧對映著梁單元上的一根纖維，如圖1所示。

圖1 DEM纖維梁模型

通過定義彈簧(纖維)的本構關系描述截面的變形與受力，進而可確定顆粒球元之間的粘結彈塑性接觸本構方程。假設每根纖維處于單向受力狀態，以平面問題為例，則梁單元截面上任一點處的應變為

ε=εN+εMz

(1)

式中εN為軸力產生的應變，εMz為繞z軸彎曲產生的應變。

εMz=κzyi

(2)

式中κz為繞z軸彎曲引起的曲率，yi為第i個纖維到截面中和軸的距離。

DEM方法采用增量形式對運動方程進行求解。在一個時步Δt內，首先可求得相鄰球元A和球元B的法向相對位移增量ΔUn以及相對轉角增量Δθ，則時步Δt內纖維的應變增量以及曲率增量可表示為

(3)

式中ΔεN為軸向應變增量，Δκz為繞z軸彎曲引起的曲率增量。

將式(3)代入式(1，2)，可得

Δε=ΔεN+Δκzyi

(4)

由式(4)可知，已知截面各纖維的應變增量，根據材料的彈塑性本構方程、加卸載準則和增量理論，即可進一步求得截面纖維的應力增量 ，其具體實現過程與通常的彈塑性有限元分析完全一致[13]，不再贅述。

為了滿足計算精度并節省計算量，將構件截面離散為若干個纖維塊，再將彈簧置于纖維塊的高斯積分點處。通過對截面所有纖維進行應力分析，可得到截面的接觸內力增量為

(5)

式中Ai為第i個纖維的面積，ΔFn為截面法向接觸力增量，ΔMz為繞z軸引起的接觸力矩增量。

上述DEM纖維梁模型，能充分考慮截面塑性的發展過程。對于破壞前塑性階段不可忽略的材料或結構而言，模擬分析結果更加精細化和更接近于工程實際。

3 構件斷裂模擬算法及實現

(6)

如前文所述，構件截面由多根彈簧組成，那么當某一時刻截面上所有彈簧均已發生斷裂，定義為全截面開裂，斷裂模式則表現為兩相鄰球元脫離分開；初始時刻，球元之間的粘結處于完好狀態。以圖2(a)所示的懸臂梁結構DEM模型(離散為7個球元)為例,在力P作用下，靠近固定端處截面所受的內力最大，故球元1和球元3之間的粘結梁截面應變(應力)最先滿足斷裂準則，在某一時刻該部位首先發生斷裂，球元3與球元1發生脫開，接觸粘結消失，如圖2(b)所示。

圖2 懸臂梁DEM斷裂過程模擬

當球元之間的粘結發生斷裂和消失后，需要將該截面的接觸力Fn和接觸力矩Mz設為零，從而使球元在下一個時步的計算中獲得內力和外力的更新。此外，斷裂過程具體模擬時，本文不考慮球元自身的破裂，同時假定構件在斷裂前后球元的總數量保持不變。

根據上述理論，基于Fortran語言開發了DEM纖維梁模型彈性-彈塑性以及構件斷裂分析計算程序。

4 懸臂梁算例

圖3 自由端受集中力作用的懸臂梁

用DEM方法模擬時，將懸臂梁離散為11個球元，球元粘結數量為10，計算時步Δt取為1.0×10-4s。當力P取值較大時，該懸臂梁會因為幾何非線性和材料非線性產生極端大變形，采用考慮截面塑性開展的纖維梁模型更能真實地模擬這種力學行為。

采用纖維梁模型時，如圖4所示，沿截面高度方向劃分為10個纖維塊，每個纖維塊上有兩根彈簧并置于高斯積分點上。

圖4 截面纖維劃分及高斯點

構件斷裂模擬時，采用式(6)作為斷裂準則，取纖維的極限應變εu=0.05。計算過程中每個加載時步均需對截面各纖維進行斷裂判別，當全截面發生斷裂后，該截面處的粘結消失，將球元之間的接觸力和接觸力矩設為零并移除外力，從而實現桿件斷裂過程的數值模擬。

圖5為不同時刻下懸臂梁的變形圖。當外力P=130.6 N時，認為懸臂梁固定端全截面發生斷裂(僅剩1根纖維未失效)，隨后撤去自由端外力且設定斷裂不再發生，懸臂梁脫離固定支座發生自由運動(圖5(a))。不考慮斷裂時，許多學者對該算例進行了研究，如Yang等[14]用有限元法對該算例進行了彈塑性大變形分析。圖5(b)為不考慮斷裂時懸臂梁在不同外力值下的變形圖，隨著荷載的增加梁自由端位移不斷增大，同時梁的變形程度也愈加顯著，可以看出，本文方法與文獻結果吻合較好，驗證了該方法的正確性和有效性。

圖5 不同時刻下懸臂梁的變形

圖6為梁固定端處截面纖維在不同荷載值下的逐漸開裂過程?？梢钥闯?，當力P=99.4 N時，截面最外層纖維(即1號和10號)首先發生斷裂；隨著荷載值的增加，截面塑性區域逐漸擴大，9號、2號、8號和3號等纖維由于達到極限應變相繼失效，截面破壞程度逐漸加??；截面由外到內逐漸開裂，最終當外力P達到130.6 N時，截面除5號纖維外其他纖維已全部失效，在此基礎上荷載稍有增加即發生了全截面斷裂。

圖6 梁固定端處截面纖維開裂過程

圖7為梁固定端處截面纖維在斷裂過程中截面抵抗彎矩與曲率的變化關系曲線。當外力P從0逐漸增加到99.4 N之前，截面各纖維均處于完好狀態，截面抵抗彎矩隨曲率增加而不斷增大，其中當截面部分纖維發生屈服后，由于剛度下降導致曲線斜率減緩；當P=99.4 N時，由于截面最外層纖維發生開裂失效，截面抵抗彎矩由最高點(452.3 Nm)急劇下降到257.4 Nm，然后又線性遞增直至下一個纖維開裂；當截面纖維逐個發生開裂并退出工作后，纖維應力設置為0，截面抵抗彎矩不斷減小并最終趨于0。

圖7 梁固定端處截面抵抗彎矩與曲率關系曲線

圖8是梁固定端相接觸的球元速度時程曲線。在截面未開裂之前，因緩慢加載，球元的速度幾乎為零，可認為是靜態的；而一旦全截面發生了開裂，球元速度瞬間急劇增大，然后隨著脫離體的自由運動不斷發生動態變化。上述計算分析過程驗證了本文所提桿件斷裂模擬算法的有效性和合理性。

圖8 梁固定端處球元速度時程曲線

5 單層網殼結構算例驗證與分析

為驗證本文方法在大跨空間網格結構中的適用性，以文獻[15]開展的單層球面網殼模型結構振動臺倒塌試驗為例進行對比分析。文獻中設計了3個K6型單層球殼結構縮尺模型，跨度均為7.5 m，矢跨比均為1/2。選取其中的模型3作為研究對象，試驗模型的幾何物理參數及地震波加載工況等參見文獻。本文用DEM模擬時，將網殼的各根桿件根據其長度分別離散成5～8個球元，離散元模型如圖9所示，共有球元4315個，球元粘結4698個。試驗模型所用的桿件均為不銹鋼鋼管，截面均為Φ14×0.6。采用纖維梁模型時，其截面纖維劃分如圖10所示。模擬焊縫處斷裂時各鋼管纖維對應的極限應力取為560 MPa[16]。

圖9 單層球殼離散元模型

圖10 圓管截面纖維劃分

圖11為單層球殼模型振動臺試驗破壞現象，當輸入峰值加速度PGA達到2108 gal時，模型從下至上第1圈和第2圈的斜桿發生彎曲并且數量持續增多，部分桿件與球節點處發生斷裂(圖11(a))；繼續增大PGA至2268 gal時，第1圈和第2圈部分桿件和球節點完全與結構脫離，網殼模型發生整體倒塌，如圖11(b)所示。

圖11 試驗模型破壞現象

圖12為本文方法模擬得到的網殼試驗模型破壞過程數值仿真結果，與振動臺試驗現象[15]比較一致?？梢钥闯?，當PGA小于1200 gal時，結構整體基本未發生變形；當PGA提高至1600 gal時，模型底部第1圈少量斜桿首先發生彎曲(圖12(a))；當PGA達到1800 gal時，由于部分焊縫先滿足斷裂準則，該部位斜桿與焊接球節點斷開(圖12(b))，但數量較少；隨著PGA的繼續增大，模型從下至上第1圈和第2圈的斜桿相繼發生彎曲，當PGA增大至2100 gal時，支座附近第1圈有近1/2斜桿與球節點脫離，節點位移顯著增大，但結構并未發生倒塌(圖12(c))；當PGA加至2300 gal時，模型底部第1圈的桿件幾乎全部與球節點發生斷裂，同時第2圈也有部分桿件與球節點脫開，結構因失去了底部桿件支撐而無法繼續承載，最終整個結構向下坍塌(圖12(d))。

圖12 試驗模型的破壞過程數值仿真結果

6 結 論

在彈性DEM方法基礎上，將粘結分布彈簧用梁纖維進行等效，推導了可考慮截面塑性開展的DEM纖維梁模型，實現了桿系結構彈塑性分析；在此基礎上，提出了桿件斷裂模擬算法，并編制了計算程序；將這一模型用于懸臂梁算例，模擬了懸臂梁從彈性小變形到彈塑性階段，再到斷裂和破壞的全過程，結果較為合理；最后將該方法應用于單層網殼倒塌破環行為模擬，分析結果表明數值仿真得到的網殼桿件斷裂過程及結構倒塌模式與網殼振動臺倒塌試驗現象吻合，進一步驗證了該方法的正確性。結論如下。

(1) 采用DEM纖維梁模型模擬桿件的斷裂過程非常直觀，只需定義彈簧的破壞準則即可實現球元之間的粘結消失與否；較FEM方法而言，不存在網格畸變和網格重劃分及迭代不收斂等問題，更具計算優勢。

(2) 本文方法雖可以描述桿件截面的漸進開裂過程，但只能描述結構的宏觀破壞模式，還無法模擬材料內部裂紋的產生、開裂及擴展。后續還需進一步將DEM方法拓展用于精細化建模及仿真模擬，研究結構的細觀或微觀破壞模式，從而可以更深入地揭示結構的倒塌破壞機理。