齊 念 葉繼紅

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

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

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

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

1 DEM法

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

(1)

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

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

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

(2)

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

(3)

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

(4)

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

(5)

(6)

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

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

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

(7)

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

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

k=max{Kn,Ks,kn,ks}

(8)

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

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

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

(9)

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

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

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

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

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

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

(10)

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ù)為

(11)

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

(12)

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

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

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

4 數(shù)值算例

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

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

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

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

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

4.2 懸臂梁的大彎曲問題

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

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

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

圖8 端部彎矩作用方式

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

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

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

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

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

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

5 結(jié)論

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

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

)

[1] Teh L H, Clarke M J. Co-rotational and lagrangian formulations for elastic three-dimensional beam finite element [J].JournalofConstructionalSteelResearch, 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].EngineeringStructures, 2007,29(6): 1189-1200.

[3] Thanh N L, Jean M B, Mohammed H. Efficient formulation for dynamics of co-rotational 2D beams [J].ComputationalMechanics, 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].AdvancesinMechanics, 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].PureandAppliedGeophysics, 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].StructuralEngineeringandMechanics, 2010,36(3): 1-19.

[8] Cundall P A.PFCuser’smanual(version3.1) [M]. Minneapolis, MN, USA: Itasca Consulting Group, Inc., 2004.

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

[10] 徐小敏,凌道盛,陳云敏，等. 基于線性接觸模型的顆粒材料細(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].ChineseJournalofGeotechnicalEngineering, 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].InternationalJournalforNumericalMethodsinEngineering, 1981,17(1): 145-153.

[12] Wu Tongyue, Wang Renzuo, Wang Chungyue. Large deflection analysis of flexible planar frames [J].JournaloftheChineseInstituteofEngineers, 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].JournalofStructuralEngineering,ASCE, 2003,129(9): 1259-1267.