基于離散元的桿件斷裂全過(guò)程數(shù)值模擬

齊 念， 鄭開(kāi)啟

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

1 引 言

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

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

本文基于已有研究，提出了一種適于結(jié)構(gòu)彈塑性分析的DEM纖維梁模型，然后構(gòu)建了桿件斷裂模擬算法，并自編了構(gòu)件斷裂計(jì)算程序；最后將其應(yīng)用于懸臂梁結(jié)構(gòu)和單層網(wǎng)殼振動(dòng)臺(tái)倒塌試驗(yàn)兩個(gè)算例，對(duì)桿件的彈性和彈塑性直至斷裂全過(guò)程進(jìn)行模擬，通過(guò)數(shù)值分析驗(yàn)證該方法的正確性和適用性。

2 DEM纖維梁模型

在桿系結(jié)構(gòu)的彈性DEM方法中，顆粒球元之間的接觸粘結(jié)看作是一根彈性梁，兩球心的距離即為彈性梁的長(zhǎng)度L，且采用構(gòu)件實(shí)際截面的幾何物理參數(shù)計(jì)算各彈簧接觸剛度系數(shù)[12]。彈塑性分析時(shí)，借鑒纖維模型的研究思路，將兩相鄰顆粒球元之間的接觸粘結(jié)用分布式彈簧進(jìn)行等效，分布彈簧看作是梁截面的若干根纖維，每一根彈簧對(duì)映著梁?jiǎn)卧系囊桓w維，如圖1所示。

圖1 DEM纖維梁模型

通過(guò)定義彈簧(纖維)的本構(gòu)關(guān)系描述截面的變形與受力，進(jìn)而可確定顆粒球元之間的粘結(jié)彈塑性接觸本構(gòu)方程。假設(shè)每根纖維處于單向受力狀態(tài)，以平面問(wèn)題為例，則梁?jiǎn)卧孛嫔先我稽c(diǎn)處的應(yīng)變?yōu)?/p>

ε=εN+εMz

(1)

式中εN為軸力產(chǎn)生的應(yīng)變，εMz為繞z軸彎曲產(chǎn)生的應(yīng)變。

εMz=κzyi

(2)

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

DEM方法采用增量形式對(duì)運(yùn)動(dòng)方程進(jìn)行求解。在一個(gè)時(shí)步Δt內(nèi)，首先可求得相鄰球元A和球元B的法向相對(duì)位移增量ΔUn以及相對(duì)轉(zhuǎn)角增量Δθ，則時(shí)步Δt內(nèi)纖維的應(yīng)變?cè)隽恳约扒试隽靠杀硎緸?/p>

(3)

式中ΔεN為軸向應(yīng)變?cè)隽浚う蕑為繞z軸彎曲引起的曲率增量。

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

Δε=ΔεN+Δκzyi

(4)

由式(4)可知，已知截面各纖維的應(yīng)變?cè)隽?，根?jù)材料的彈塑性本構(gòu)方程、加卸載準(zhǔn)則和增量理論，即可進(jìn)一步求得截面纖維的應(yīng)力增量 ，其具體實(shí)現(xiàn)過(guò)程與通常的彈塑性有限元分析完全一致[13]，不再贅述。

為了滿(mǎn)足計(jì)算精度并節(jié)省計(jì)算量，將構(gòu)件截面離散為若干個(gè)纖維塊，再將彈簧置于纖維塊的高斯積分點(diǎn)處。通過(guò)對(duì)截面所有纖維進(jìn)行應(yīng)力分析，可得到截面的接觸內(nèi)力增量為

(5)

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

上述DEM纖維梁模型，能充分考慮截面塑性的發(fā)展過(guò)程。對(duì)于破壞前塑性階段不可忽略的材料或結(jié)構(gòu)而言，模擬分析結(jié)果更加精細(xì)化和更接近于工程實(shí)際。

3 構(gòu)件斷裂模擬算法及實(shí)現(xiàn)

(6)

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

圖2 懸臂梁DEM斷裂過(guò)程模擬

當(dāng)球元之間的粘結(jié)發(fā)生斷裂和消失后，需要將該截面的接觸力Fn和接觸力矩Mz設(shè)為零，從而使球元在下一個(gè)時(shí)步的計(jì)算中獲得內(nèi)力和外力的更新。此外，斷裂過(guò)程具體模擬時(shí)，本文不考慮球元自身的破裂，同時(shí)假定構(gòu)件在斷裂前后球元的總數(shù)量保持不變。

根據(jù)上述理論，基于Fortran語(yǔ)言開(kāi)發(fā)了DEM纖維梁模型彈性-彈塑性以及構(gòu)件斷裂分析計(jì)算程序。

4 懸臂梁算例

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

用DEM方法模擬時(shí)，將懸臂梁離散為11個(gè)球元，球元粘結(jié)數(shù)量為10，計(jì)算時(shí)步Δt取為1.0×10-4s。當(dāng)力P取值較大時(shí)，該懸臂梁會(huì)因?yàn)閹缀畏蔷€(xiàn)性和材料非線(xiàn)性產(chǎn)生極端大變形，采用考慮截面塑性開(kāi)展的纖維梁模型更能真實(shí)地模擬這種力學(xué)行為。

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

圖4 截面纖維劃分及高斯點(diǎn)

構(gòu)件斷裂模擬時(shí)，采用式(6)作為斷裂準(zhǔn)則，取纖維的極限應(yīng)變?chǔ)舥=0.05。計(jì)算過(guò)程中每個(gè)加載時(shí)步均需對(duì)截面各纖維進(jìn)行斷裂判別，當(dāng)全截面發(fā)生斷裂后，該截面處的粘結(jié)消失，將球元之間的接觸力和接觸力矩設(shè)為零并移除外力，從而實(shí)現(xiàn)桿件斷裂過(guò)程的數(shù)值模擬。

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

圖5 不同時(shí)刻下懸臂梁的變形

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

圖6 梁固定端處截面纖維開(kāi)裂過(guò)程

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

圖7 梁固定端處截面抵抗彎矩與曲率關(guān)系曲線(xiàn)

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

圖8 梁固定端處球元速度時(shí)程曲線(xiàn)

5 單層網(wǎng)殼結(jié)構(gòu)算例驗(yàn)證與分析

為驗(yàn)證本文方法在大跨空間網(wǎng)格結(jié)構(gòu)中的適用性，以文獻(xiàn)[15]開(kāi)展的單層球面網(wǎng)殼模型結(jié)構(gòu)振動(dòng)臺(tái)倒塌試驗(yàn)為例進(jìn)行對(duì)比分析。文獻(xiàn)中設(shè)計(jì)了3個(gè)K6型單層球殼結(jié)構(gòu)縮尺模型，跨度均為7.5 m，矢跨比均為1/2。選取其中的模型3作為研究對(duì)象，試驗(yàn)?zāi)Ｐ偷膸缀挝锢韰?shù)及地震波加載工況等參見(jiàn)文獻(xiàn)。本文用DEM模擬時(shí)，將網(wǎng)殼的各根桿件根據(jù)其長(zhǎng)度分別離散成5～8個(gè)球元，離散元模型如圖9所示，共有球元4315個(gè)，球元粘結(jié)4698個(gè)。試驗(yàn)?zāi)Ｐ退玫臈U件均為不銹鋼鋼管，截面均為Φ14×0.6。采用纖維梁模型時(shí)，其截面纖維劃分如圖10所示。模擬焊縫處斷裂時(shí)各鋼管纖維對(duì)應(yīng)的極限應(yīng)力取為560 MPa[16]。

圖9 單層球殼離散元模型

圖10 圓管截面纖維劃分

圖11為單層球殼模型振動(dòng)臺(tái)試驗(yàn)破壞現(xiàn)象，當(dāng)輸入峰值加速度PGA達(dá)到2108 gal時(shí)，模型從下至上第1圈和第2圈的斜桿發(fā)生彎曲并且數(shù)量持續(xù)增多，部分桿件與球節(jié)點(diǎn)處發(fā)生斷裂(圖11(a))；繼續(xù)增大PGA至2268 gal時(shí)，第1圈和第2圈部分桿件和球節(jié)點(diǎn)完全與結(jié)構(gòu)脫離，網(wǎng)殼模型發(fā)生整體倒塌，如圖11(b)所示。

圖11 試驗(yàn)?zāi)Ｐ推茐默F(xiàn)象

圖12為本文方法模擬得到的網(wǎng)殼試驗(yàn)?zāi)Ｐ推茐倪^(guò)程數(shù)值仿真結(jié)果，與振動(dòng)臺(tái)試驗(yàn)現(xiàn)象[15]比較一致。可以看出，當(dāng)PGA小于1200 gal時(shí)，結(jié)構(gòu)整體基本未發(fā)生變形；當(dāng)PGA提高至1600 gal時(shí)，模型底部第1圈少量斜桿首先發(fā)生彎曲(圖12(a))；當(dāng)PGA達(dá)到1800 gal時(shí)，由于部分焊縫先滿(mǎn)足斷裂準(zhǔn)則，該部位斜桿與焊接球節(jié)點(diǎn)斷開(kāi)(圖12(b))，但數(shù)量較少；隨著PGA的繼續(xù)增大，模型從下至上第1圈和第2圈的斜桿相繼發(fā)生彎曲，當(dāng)PGA增大至2100 gal時(shí)，支座附近第1圈有近1/2斜桿與球節(jié)點(diǎn)脫離，節(jié)點(diǎn)位移顯著增大，但結(jié)構(gòu)并未發(fā)生倒塌(圖12(c))；當(dāng)PGA加至2300 gal時(shí)，模型底部第1圈的桿件幾乎全部與球節(jié)點(diǎn)發(fā)生斷裂，同時(shí)第2圈也有部分桿件與球節(jié)點(diǎn)脫開(kāi)，結(jié)構(gòu)因失去了底部桿件支撐而無(wú)法繼續(xù)承載，最終整個(gè)結(jié)構(gòu)向下坍塌(圖12(d))。

圖12 試驗(yàn)?zāi)Ｐ偷钠茐倪^(guò)程數(shù)值仿真結(jié)果

6 結(jié) 論

在彈性DEM方法基礎(chǔ)上，將粘結(jié)分布彈簧用梁纖維進(jìn)行等效，推導(dǎo)了可考慮截面塑性開(kāi)展的DEM纖維梁模型，實(shí)現(xiàn)了桿系結(jié)構(gòu)彈塑性分析；在此基礎(chǔ)上，提出了桿件斷裂模擬算法，并編制了計(jì)算程序；將這一模型用于懸臂梁算例，模擬了懸臂梁從彈性小變形到彈塑性階段，再到斷裂和破壞的全過(guò)程，結(jié)果較為合理；最后將該方法應(yīng)用于單層網(wǎng)殼倒塌破環(huán)行為模擬，分析結(jié)果表明數(shù)值仿真得到的網(wǎng)殼桿件斷裂過(guò)程及結(jié)構(gòu)倒塌模式與網(wǎng)殼振動(dòng)臺(tái)倒塌試驗(yàn)現(xiàn)象吻合，進(jìn)一步驗(yàn)證了該方法的正確性。結(jié)論如下。

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

(2) 本文方法雖可以描述桿件截面的漸進(jìn)開(kāi)裂過(guò)程，但只能描述結(jié)構(gòu)的宏觀破壞模式，還無(wú)法模擬材料內(nèi)部裂紋的產(chǎn)生、開(kāi)裂及擴(kuò)展。后續(xù)還需進(jìn)一步將DEM方法拓展用于精細(xì)化建模及仿真模擬，研究結(jié)構(gòu)的細(xì)觀或微觀破壞模式，從而可以更深入地揭示結(jié)構(gòu)的倒塌破壞機(jī)理。