基于離散元法的桿系結(jié)構(gòu)幾何非線(xiàn)性大變形分析

齊 念 葉繼紅

(東南大學(xué)混凝土及預(yù)應(yīng)力混凝土結(jié)構(gòu)教育部重點(diǎn)實(shí)驗(yàn)室，南京 210096)

在工程實(shí)際分析時(shí)，經(jīng)常會(huì)遇到由幾何非線(xiàn)性所引起的結(jié)構(gòu)大變形問(wèn)題，這類(lèi)問(wèn)題的特點(diǎn)是伴隨大位移和大轉(zhuǎn)動(dòng)，因此必須考慮變形對(duì)平衡的影響.對(duì)于桿系結(jié)構(gòu)幾何非線(xiàn)性大變形問(wèn)題，有限元法是目前最為常用的分析方法，學(xué)者們也提出了多種計(jì)算格式，如完全拉格朗日(TL)法、更新拉格朗日(UL)法和協(xié)同轉(zhuǎn)動(dòng)法等［1-3］.但是，這些傳統(tǒng)方法在處理幾何非線(xiàn)性行為時(shí)，需要不斷集成和修正切線(xiàn)剛度矩陣，經(jīng)常會(huì)遇到剛度矩陣奇異或非線(xiàn)性方程迭代求解不易收斂等問(wèn)題，需要對(duì)方法本身進(jìn)行特殊處理，導(dǎo)致計(jì)算效率偏低.

離散元法(discrete element method，DEM)最早由美國(guó)學(xué)者Cundall等［4］提出，現(xiàn)已發(fā)展成為計(jì)算散體力學(xué)領(lǐng)域中一種新的數(shù)值方法.該方法直接應(yīng)用牛頓第二定律，單元之間采用彈簧系統(tǒng)連接，接觸力與接觸位移之間的關(guān)系構(gòu)成了DEM的接觸本構(gòu)模型.由于不要求滿(mǎn)足位移連續(xù)和變形協(xié)調(diào)條件，因此它特別適用于各種非連續(xù)、非均勻以及結(jié)構(gòu)大變形和失效破壞等復(fù)雜過(guò)程及其機(jī)理的研究［5］.最初DEM 法的研究對(duì)象主要是散土體材料，如今已在許多工程領(lǐng)域如巖體邊坡滑動(dòng)［6］、鋼筋混凝土梁的斷裂模擬［7］等方面取得了較好的應(yīng)用.

DEM法在桿系結(jié)構(gòu)中的應(yīng)用還較少.本文提出將DEM法用于求解桿系結(jié)構(gòu)(二維、三維)的幾何非線(xiàn)性大變形問(wèn)題.對(duì)3個(gè)經(jīng)典算例的靜力大變形問(wèn)題及非線(xiàn)性動(dòng)力行為進(jìn)行了模擬，結(jié)果表明DEM法適宜于處理結(jié)構(gòu)大變形問(wèn)題.

1 DEM法

DEM法的基本思想是把研究對(duì)象離散成剛性單元的集合.任取一個(gè)單元α，設(shè)有n個(gè)單元與其相鄰，作用在單元 α上的外力為 Fext，外力矩為Mext.根據(jù)牛頓第二定律，其運(yùn)動(dòng)控制方程為

式中，m，J分別為單元α的質(zhì)量和轉(zhuǎn)動(dòng)慣量;r，ω分別為單元α的位置矢量和角速度矢量;分別為相鄰單元j對(duì)單元α產(chǎn)生的接觸力和接觸力矩;t為時(shí)間.

1.1 單元間的接觸本構(gòu)模型

圖1 相鄰單元之間的接觸模型

在接觸平面內(nèi)將接觸力和接觸力矩按矢量法則進(jìn)行分解，得

將F和M分別向3個(gè)坐標(biāo)軸進(jìn)行投影，得

利用離散元法計(jì)算接觸力和接觸力矩時(shí)采用增量形式.在時(shí)間步長(zhǎng)Δt內(nèi)，由位移增量所引起的接觸力增量為

式中，Kn，Ks分別為彈簧法向剛度和切向剛度系數(shù);ni(i=x，y，z)為 n 在3個(gè)坐標(biāo)軸上的分量;分別為法向和切向相對(duì)位移增量，可通過(guò)接觸點(diǎn)C處的速度計(jì)算得到.

接觸力矩的計(jì)算過(guò)程與接觸力的計(jì)算過(guò)程相似.在時(shí)間步長(zhǎng)Δt內(nèi)，球A與球B之間的相對(duì)轉(zhuǎn)角增量所引起的接觸力矩增量為

式中，kn，ks分別為彈簧轉(zhuǎn)動(dòng)法向剛度和切向剛度系數(shù);分別為法向和切向相對(duì)轉(zhuǎn)角增量.

1.2 時(shí)間步長(zhǎng)的確定

在求解運(yùn)動(dòng)方程時(shí)，DEM通常采用顯示中心差分算法.中心差分法是條件穩(wěn)定算法，其計(jì)算時(shí)間步長(zhǎng)Δt必須小于由該問(wèn)題所決定的某個(gè)臨界值Δtcr，否則算法將不穩(wěn)定.

臨界值Δtcr的確定應(yīng)同時(shí)滿(mǎn)足結(jié)構(gòu)物理步長(zhǎng)和數(shù)學(xué)計(jì)算步長(zhǎng)的要求.物理步長(zhǎng)是指結(jié)構(gòu)中彈性波來(lái)回一次所需要的時(shí)間，框架結(jié)構(gòu)的物理步長(zhǎng)約為1 ms［9］.數(shù)學(xué)計(jì)算步長(zhǎng)是指積分計(jì)算所需的時(shí)間步長(zhǎng)，對(duì)中心差分法而言，解的穩(wěn)定性條件為

式中，ωn為結(jié)構(gòu)的最高階固有振動(dòng)頻率;m為結(jié)構(gòu)質(zhì)量;k為結(jié)構(gòu)剛度.

本文在確定數(shù)學(xué)計(jì)算步長(zhǎng)時(shí)，按如下方法粗略估計(jì)臨界值Δtcr:將式(7)中的k值取為1.1節(jié)中提到的4個(gè)彈簧剛度系數(shù)中的最大值，即

最后，通過(guò)比較物理步長(zhǎng)與數(shù)學(xué)計(jì)算步長(zhǎng)的臨界值，取二者較小值，以確定時(shí)間步長(zhǎng)Δt的大小.

2 動(dòng)力與靜力解及大變形問(wèn)題

與傳統(tǒng)的結(jié)構(gòu)靜力分析不同，離散元法本質(zhì)上是求解結(jié)構(gòu)的動(dòng)力問(wèn)題，因?yàn)槭?1)的解本身就是結(jié)構(gòu)的動(dòng)力反應(yīng).

對(duì)于需要計(jì)算靜力解的問(wèn)題，可采用阻尼耗能或緩慢施加外荷載的方式來(lái)逼近結(jié)構(gòu)的靜止?fàn)顟B(tài).本文是在運(yùn)動(dòng)方程式中添加一項(xiàng)阻尼力，利用阻尼產(chǎn)生的逆向運(yùn)動(dòng)來(lái)達(dá)到削減結(jié)構(gòu)動(dòng)力響應(yīng)的目的.這里采用局部非線(xiàn)性黏滯阻尼形式［8］，其本質(zhì)是在單元的合外力基礎(chǔ)上，總體施加一個(gè)與其方向相反的作用力.單元α的阻尼力可表示為

式中，為阻尼系數(shù)，取值范圍為0～1;F*為單元α的廣義力;V*為單元α的廣義速度.需要說(shuō)明的是，此處的阻尼并非結(jié)構(gòu)的真實(shí)阻尼，它實(shí)際上是一個(gè)虛擬的行為，其目的是為了獲得結(jié)構(gòu)的靜力解而選擇的一種耗能機(jī)制.

分析傳統(tǒng)結(jié)構(gòu)時(shí)，小變形和大變形問(wèn)題是需要進(jìn)行嚴(yán)格區(qū)分的.前者的幾何方程為線(xiàn)性的，而后者的幾何方程則為非線(xiàn)性的.利用有限元法處理大變形問(wèn)題時(shí)，一般是對(duì)單元切線(xiàn)剛度矩陣進(jìn)行修正，采用增量迭代的方法求解結(jié)構(gòu)的響應(yīng).而利用DEM法求解此類(lèi)問(wèn)題時(shí)，并不需要刻意區(qū)分是小變形還是大變形，因?yàn)樵诮⑦\(yùn)動(dòng)控制方程時(shí)并不涉及到幾何方程，也不要求位移連續(xù)，無(wú)需組集剛度矩陣和迭代求解，這與有限元法有著本質(zhì)的區(qū)別.因此，將DEM法用于結(jié)構(gòu)分析時(shí)，可以采用統(tǒng)一的步驟分析結(jié)構(gòu)的小變形和大變形問(wèn)題，其計(jì)算流程圖見(jiàn)圖2.

3 彈簧接觸剛度系數(shù)的確定

圖2 DEM法的計(jì)算流程圖

根據(jù)DEM法的基本理論，利用式(4)和(5)計(jì)算接觸力和接觸力矩增量時(shí)，需要事先給定彈簧接觸剛度系數(shù)的值.本文的研究對(duì)象主要是連續(xù)體——框架結(jié)構(gòu)，而傳統(tǒng)的DEM法彈簧接觸剛度系數(shù)計(jì)算公式是基于散粒體的，不適用于框架結(jié)構(gòu)［10］.基于簡(jiǎn)單梁理論，將通過(guò)彈簧連接的2個(gè)圓球比擬成1根梁(見(jiàn)圖3).圖3中，RA，RB分別為球A和球B的半徑;L=RA+RB為梁的長(zhǎng)度.設(shè)結(jié)構(gòu)構(gòu)件的截面積為S，截面慣性矩為I，扭轉(zhuǎn)慣性矩為J，彈性模量和剪切模量分別為E和G.

圖3 彈簧接觸剛度的等價(jià)梁模型

根據(jù)結(jié)構(gòu)力學(xué)中桿端力與桿端位移之間的關(guān)系，可以得到彈簧接觸剛度計(jì)算公式.梁的軸向剛度系數(shù)Kn和切向剛度系數(shù)Ks分別為

2個(gè)單元之間發(fā)生相對(duì)轉(zhuǎn)動(dòng)產(chǎn)生轉(zhuǎn)角Δθ，Δθ可分解成平行于接觸面內(nèi)的切向分量Δθs和垂直于接觸面內(nèi)的法向分量Δθn.切向分量Δθs所引起的接觸力矩，與一端固定一端滑動(dòng)的梁所產(chǎn)生的桿端彎矩等效.因此，轉(zhuǎn)動(dòng)法向剛度系數(shù)為

垂直于接觸面內(nèi)的轉(zhuǎn)角法向分量Δθn引起的則是扭矩.若是平面問(wèn)題，扭矩為0;對(duì)于一般空間問(wèn)題，扭矩可根據(jù)截面扭轉(zhuǎn)剛度與扭轉(zhuǎn)角相乘得到.因此，轉(zhuǎn)動(dòng)切向剛度系數(shù)為

利用1個(gè)空間懸臂梁對(duì)該方法進(jìn)行驗(yàn)證.梁長(zhǎng)L=200 mm，其截面呈圓環(huán)形，外徑φ=8 mm，壁厚 δ=0.5 mm，彈性模量 E=200 kN/mm2，剪切模量G=80 kN/mm2.梁端部承受集中力F=10 N，彎矩M=100 N·mm，扭矩T=2 kN·mm，材料為線(xiàn)彈性.懸臂梁DEM計(jì)算模型如圖4所示.在模擬時(shí)將最左端圓球固定，因?yàn)槭乔蠼忪o力問(wèn)題，取阻尼系數(shù) =0.7，時(shí)間步長(zhǎng) Δt=1 ms.

圖4 懸臂梁DEM計(jì)算模型

懸臂梁的變形和內(nèi)力計(jì)算結(jié)果見(jiàn)表1.由表可知，根據(jù)DEM法計(jì)算的懸臂梁變形及結(jié)構(gòu)內(nèi)力與理論解相比，最大誤差僅為0.3%，說(shuō)明本文確定彈簧接觸剛度系數(shù)及Δt的方法是正確合理的.

表1 集中力和力矩作用下的懸臂梁計(jì)算結(jié)果

4 數(shù)值算例

根據(jù)DEM法的基本原理及計(jì)算流程，編制了空間框架結(jié)構(gòu)幾何非線(xiàn)性大變形分析程序.通過(guò)對(duì)文獻(xiàn)中常被引用的3個(gè)典型算例進(jìn)行模擬和分析，驗(yàn)證本文方法的正確性和適用性.

4.1 平面正方形剛架的大變形分析

如圖5(a)所示，平面正方形剛架各構(gòu)件的材料和尺寸均相同，相關(guān)幾何及物理參數(shù)為:l=10 cm，S=0.5 cm2，I=0.0104 cm4，E=16 MN/cm2，材料為線(xiàn)彈性，P為外力.DEM計(jì)算模型如圖5(b)所示，球半徑取為1 cm，時(shí)間步長(zhǎng)Δt=0.5 ms，阻尼系數(shù) =0.7，材料密度取為單位1.

圖5 平面正方形剛架及離散元模型

圖6為利用DEM法計(jì)算得到的荷載-變形曲線(xiàn).與文獻(xiàn)［11］中采用橢圓積分計(jì)算的結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)位移及轉(zhuǎn)角的計(jì)算值均與文獻(xiàn)解相吻合，最大誤差不超過(guò)1%.

圖6 平面正方形剛架對(duì)邊受拉下的荷載-變形曲線(xiàn)

4.2 懸臂梁的大彎曲問(wèn)題

如圖7所示，在懸臂梁端部作用一集中彎矩M(t).彎矩隨時(shí)間逐漸增加，梁將會(huì)由原來(lái)的靜止?fàn)顟B(tài)彎成圓形或螺旋形.該典型算例常被用于結(jié)構(gòu)大變形問(wèn)題的分析中［12-13］.梁的相關(guān)幾何參數(shù)及物理參數(shù)采用文獻(xiàn)［12］中的數(shù)據(jù)，彎矩作用方式采用如圖8所示的線(xiàn)性漸增荷載，DEM模擬時(shí)將梁離散成9個(gè)半徑為1.25 mm的圓球.

圖7 受端彎矩作用的懸臂梁

下面分2種工況對(duì)該梁的大彎曲行為進(jìn)行數(shù)值模擬.

圖8 端部彎矩作用方式

工況1 對(duì)梁進(jìn)行緩慢加載，如圖8(a)所示.取Δt=1 ms，=0.7.在該彎矩作用下梁最終將彎成一個(gè)圓.由于加載速率很慢，同時(shí)又有阻尼的耗能作用，因此實(shí)際上是模擬結(jié)構(gòu)的近似靜力變形行為.表2為利用DEM法得到的梁自由端計(jì)算結(jié)果，與解析解［13］相比，兩者最大誤差為 0.6%，說(shuō)明本文方法具有較高的精度.

表2 彎矩作用下懸臂梁自由端計(jì)算結(jié)果

工況2 梁的幾何及材料參數(shù)與工況1相同，但梁端彎矩加載方式如圖8(b)所示，阻尼系數(shù) =0.1.利用DEM 法對(duì)其進(jìn)行模擬，圖9是懸臂梁分別于時(shí)間 t=9，10，11，13 s時(shí)的形狀，發(fā)現(xiàn)梁最后卷曲成螺旋狀，而非工況1中的圓形.這是因?yàn)榇藭r(shí)的計(jì)算結(jié)果為動(dòng)態(tài)解，與所采用的阻尼系數(shù)和加載速率有關(guān)，若保持最終加載力不變，隨著時(shí)間的延續(xù)，由于阻尼耗能，動(dòng)態(tài)解最終也會(huì)趨于靜力解，螺旋狀也將趨于圓狀.

圖9 不同時(shí)刻懸臂梁的形狀

4.3 六角形空間剛架的大變形分析

對(duì)于如圖10所示的六角形空間剛架，結(jié)構(gòu)的6個(gè)支座均為鉸約束，各構(gòu)件截面均為正方形，其幾何物理參數(shù)的取值與文獻(xiàn)［14］保持一致.利用DEM法進(jìn)行模擬時(shí)，將構(gòu)件離散成24個(gè)球單元，取 Δt=0.1 ms， =0.7，材料密度取為單位 1.

圖10 六角形空間剛架(單位:mm)

對(duì)該結(jié)構(gòu)進(jìn)行幾何非線(xiàn)性靜力大變形分析.表3為利用DEM法得到的剛架頂點(diǎn)荷載-位移計(jì)算結(jié)果.由表可知，與文獻(xiàn)［14］中利用有限元法分析的結(jié)果進(jìn)行對(duì)比，兩者最大誤差不超過(guò)0.9%.

表3 六角形空間剛架頂點(diǎn)的荷載-位移結(jié)果

5 結(jié)論

1)本文將DEM法推廣應(yīng)用于桿系結(jié)構(gòu)(二維、三維)的幾何非線(xiàn)性大變形問(wèn)題.基于簡(jiǎn)單梁理論，推導(dǎo)了適用于桿系結(jié)構(gòu)分析的彈簧接觸剛度系數(shù)計(jì)算公式，并利用算例對(duì)其正確性進(jìn)行了檢驗(yàn).

2)基于本文推導(dǎo)的剛度系數(shù)計(jì)算公式，對(duì)桿系結(jié)構(gòu)的幾何非線(xiàn)性大變形問(wèn)題進(jìn)行分析.列出了3個(gè)典型數(shù)值算例，即2個(gè)平面框架和1個(gè)空間網(wǎng)格結(jié)構(gòu)，分別對(duì)其靜力和動(dòng)力大變形行為進(jìn)行了模擬，并將結(jié)果與其他計(jì)算方法的結(jié)果進(jìn)行比較，兩者吻合良好.此外，離散元法在處理幾何非線(xiàn)性時(shí)無(wú)需組集剛度矩陣，也不用迭代求解非線(xiàn)性方程，故該方法適宜于處理?xiàng)U系結(jié)構(gòu)的大變形問(wèn)題.

3)DEM法本質(zhì)上是求解結(jié)構(gòu)的動(dòng)力行為，對(duì)于需要計(jì)算靜力解的問(wèn)題則是通過(guò)增加阻尼項(xiàng)來(lái)進(jìn)行模擬的.阻尼系數(shù) 的選取應(yīng)根據(jù)分析問(wèn)題的性質(zhì)合理確定，對(duì)于靜力解問(wèn)題，綜合考慮數(shù)值精度和計(jì)算效率，建議可取 =0.7.由于DEM 法是采用中心差分法求解運(yùn)動(dòng)方程，計(jì)算時(shí)間步長(zhǎng)Δt的值一般取0.1～1 ms.所提的時(shí)間步長(zhǎng)臨界值估算方法，可供結(jié)構(gòu)分析時(shí)借鑒.

References)

［1］Teh L H，Clarke M J.Co-rotational and lagrangian formulations for elastic three-dimensional beam finite element［J］.Journal of Constructional Steel Research，1998，48(3):123-144.

［2］Yang Y B，Lin S P，Leu L J.Solution strategy and rigid element for nonlinear analysis of elastically structures based on updated Lagrangian formulation［J］.Engineering Structures，2007，29(6):1189-1200.

［3］Thanh N L，Jean M B，Mohammed H.Efficient formulation for dynamics of co-rotational 2D beams［J］.Computational Mechanics，2011，48(2):153-161.

［4］Cundall P A，Strack O D L.A discrete numerical model for granular assemblies［J］.Geotechnique，1979，29(1):47-65.

［5］劉凱欣，高凌天.離散元法研究的評(píng)述［J］.力學(xué)進(jìn)展，2003，33(4):483-490.Liu Kaixin，Gao Lingtian.A review of the discrete element method ［J］.Advances in Mechanics，2003，33(4):483-490.(in Chinese)

［6］James F H，David S C，William S P.Simulation of unstable fault slip in Granite using a bonded-particle model［J］.Pure and Applied Geophysics，2002，159(1/2/3):221-245.

［7］Azevedo N M，Lemos J V，Almeida J R.A discrete particle model for reinforced concrete fracture analysis［J］.Structural Engineering and Mechanics，2010，36(3):1-19.

［8］Cundall P A.PFC user's manual(version 3.1)［M］.Minneapolis，MN，USA:Itasca Consulting Group，Inc.，2004.

［9］喻瑩.基于有限質(zhì)點(diǎn)法的空間鋼結(jié)構(gòu)連續(xù)倒塌破壞研究［D］.杭州:浙江大學(xué)建筑工程學(xué)院，2010.

［10］徐小敏，凌道盛，陳云敏，等.基于線(xiàn)性接觸模型的顆粒材料細(xì)宏觀彈性常數(shù)相關(guān)關(guān)系研究［J］.巖土工程學(xué)報(bào)，2010，32(7):991-998.Xu Xiaomin，Ling Daosheng，Chen Yunmin，et al.Correlation of microscopic and macroscopic elastic constants of granular materials based on linear contact model［J］.Chinese Journal of Geotechnical Engineering，2010，32(7):991-998.(in Chinese)

［11］Kjell M.Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals［J］.International Journal for Numerical Methods in Engineering，1981，17(1):145-153.

［12］Wu Tongyue，Wang Renzuo，Wang Chungyue.Large deflection analysis of flexible planar frames［J］.Journal of the Chinese Institute of Engineers，2006，29(4):593-606.

［13］丁承先，段元鋒，吳東岳.向量式結(jié)構(gòu)力學(xué)［M］.北京:科學(xué)出版社，2012:273-278.

［14］Shugyo M.Elastoplastic large deflection analysis of three-dimensional steel frames［J］.Journal of Structural Engineering，ASCE，2003，129(9):1259-1267.