顆粒離散元法在框架結(jié)構(gòu)中的應(yīng)用

齊 念

(中冶華天南京工程技術(shù)有限公司，江蘇 南京 210019)

顆粒離散元法(Discrete Element Method，DEM)是20世紀(jì)70年代由美國(guó)學(xué)者Cundall教授[1]首先提出的，該方法最初應(yīng)用于土體、巖石等材料的力學(xué)性能分析。其基本思想是把研究對(duì)象離散為剛性單元的集合，單元與單元之間通過(guò)彈簧連接，接觸力與接觸位移之間的關(guān)系構(gòu)成了DEM的接觸本構(gòu)模型。各個(gè)剛性單元滿足牛頓運(yùn)動(dòng)方程，DEM法分析時(shí)允許單元發(fā)生相對(duì)運(yùn)動(dòng)，不要求必須滿足位移連續(xù)和變形協(xié)調(diào)條件，尤其適合非線性問(wèn)題的研究。目前DEM法已廣泛用于計(jì)算散體力學(xué)領(lǐng)域[2,3]，如粉末的加工、研磨，藥品的運(yùn)輸?shù)?。本文通過(guò)理論推導(dǎo)和數(shù)值分析，將DEM法推廣用于框架結(jié)構(gòu)幾何非線性大變形問(wèn)題。然后通過(guò)一個(gè)經(jīng)典算例的大變形行為進(jìn)行模擬和分析并與其他方法的計(jì)算結(jié)果進(jìn)行比較。結(jié)果表明，DEM方法適用于框架結(jié)構(gòu)大變形問(wèn)題的分析。

1 顆粒DEM法基本原理

1.1 單元之間的接觸力分析

顆粒離散元，以顆粒為基本研究對(duì)象，將物體離散作為具有代表性的數(shù)個(gè)單元，利用顆粒流模型構(gòu)建物體的力學(xué)性質(zhì)。以平面問(wèn)題為例，假定單元形狀為圓盤(pán)，單元之間的接觸采用柔性接觸，即允許接觸處產(chǎn)生重疊，圖1為任意兩個(gè)相鄰圓A和圓B間的接觸模型，C點(diǎn)為接觸中心。

當(dāng)圓A與圓B之間產(chǎn)生接觸時(shí)，接觸重疊量Un可表示為：

Un=R[A]+R[B]-D

(1)

其中，R[A]，R[B]，D分別為圓A、圓B的半徑及兩圓之間的中心距。

在離散元方法中，相鄰圓之間的相互接觸作用力是通過(guò)接觸彈簧來(lái)實(shí)現(xiàn)[4]。如圖1所示，在接觸平面內(nèi)，接觸力F可分解成垂直于接觸面的法向力Fn和平行于接觸面的切向力Fs，即：

F=Fn+Fs=Fn×n+Fs×s

(2)

其中，F(xiàn)n，F(xiàn)s分別為法向力和切向力的大小；n，s分別為接觸平面的法向和切向單位矢量。

此外，接觸點(diǎn)C處因兩圓盤(pán)發(fā)生相對(duì)轉(zhuǎn)動(dòng)還會(huì)引起大小為Mz的接觸力矩，這與散粒體離散元中的接觸模型不同。

將接觸力分別向坐標(biāo)軸進(jìn)行投影，得：

Fi=Fn×ni+Fs×si(i=1,2)

(3)

根據(jù)力系平移定理，將Fi移至各單元的中心處，有：

(4)

(5)

1.2 接觸力和接觸力矩的計(jì)算

接觸力和接觸力矩的計(jì)算一般是通過(guò)力與位移之間的關(guān)系來(lái)表達(dá)。對(duì)于法向接觸力Fn：

ΔFs=-KnΔUn

(6)

切向接觸力Fs，用增量的形式進(jìn)行計(jì)算：

ΔFs=-KsΔUs

(7)

ΔUs=VsΔt

(8)

其中，Kn,Ks分別為彈簧法向剛度和切向剛度系數(shù);ΔUs為相對(duì)切向位移增量;Δt為計(jì)算時(shí)步;Vs為接觸點(diǎn)C處的切向速度。

對(duì)于接觸力矩M3的計(jì)算，也采用增量形式來(lái)表達(dá)：

(9)

其中，Kθ為彈簧轉(zhuǎn)動(dòng)剛度;Δθ為相對(duì)轉(zhuǎn)角增量;ω[B],ω[A]分別為圓A和圓B的轉(zhuǎn)動(dòng)角速度。

然后，按以下方式進(jìn)行接觸作用力的更迭計(jì)算：接觸形成時(shí)，對(duì)Fn,Fs,M3初始化為零，然后進(jìn)行下一時(shí)步的計(jì)算，將由相對(duì)位移引起的接觸力增量、相對(duì)轉(zhuǎn)角引起的接觸力矩增量累加到現(xiàn)值上，即：

(10)

1.3 運(yùn)動(dòng)控制方程

(11)

對(duì)于運(yùn)動(dòng)方程的求解，DEM通常采用顯示積分算法，本文采取中心差分格式。

2 算例驗(yàn)證

平面正方形剛架大變形分析。

如圖2所示，一平面正方形剛架在對(duì)邊受拉力作用下發(fā)生大位移大轉(zhuǎn)動(dòng)行為，采用DEM法模擬并與文獻(xiàn)[5]中采用橢圓積分計(jì)算的結(jié)果進(jìn)行對(duì)比。剛架各構(gòu)件的材料和尺寸均相同，相關(guān)幾何及物理參數(shù)為：L=10 cm，A=0.5 cm2，I=0.010 4 cm4，E=1.6×107N/cm2，材料為線彈性。DEM建模時(shí)，將剛架各邊構(gòu)件均離散成11個(gè)半徑為1.0 cm的圓盤(pán)。計(jì)算時(shí)步Δt=1.0×10-3s，阻尼系數(shù)ζ取0.6。

表1為用DEM法計(jì)算得到的剛架在對(duì)邊受拉情形下的荷載—變形結(jié)果，定義荷載無(wú)量綱因子：λ=PL2/EI。與文獻(xiàn)中KJELL采用橢圓積分分析進(jìn)行對(duì)比發(fā)現(xiàn)，兩種方法計(jì)算的結(jié)果相當(dāng)接近，最大誤差不超過(guò)0.9%，表明本文方法具有較高的數(shù)值精度。

表1 平面正方形剛架對(duì)邊受拉下的荷載—變形結(jié)果

圖3為不同荷載大小作用下平面正方形剛架的最終變形圖。從圖3中看出，荷載值越大，結(jié)構(gòu)變形越顯著，用DEM分析時(shí)很容易獲得結(jié)構(gòu)的最終變形，因?yàn)樗?jì)算時(shí)會(huì)自動(dòng)考慮幾何非線性的影響，而不去區(qū)分是大變形還是小變形，因此求得的解就是它真實(shí)的狀態(tài)。

3 結(jié)語(yǔ)

本文將散體力學(xué)中應(yīng)用廣泛的顆粒DEM法推廣用于求解平面框架結(jié)構(gòu)幾何非線性大變形問(wèn)題。算例分析結(jié)果表明：本文方法計(jì)算精度較高，對(duì)框架結(jié)構(gòu)幾何非線性大變形(大位移、大轉(zhuǎn)動(dòng))問(wèn)題能很好的處理。

DEM法是以牛頓第二定律為力學(xué)基礎(chǔ)，在求解結(jié)構(gòu)大變形問(wèn)題時(shí)，不用組集單元的剛度矩陣，不需要迭代，流程相對(duì)簡(jiǎn)單。另外，該方法不要求滿足位移連續(xù)和變形協(xié)調(diào)條件，突破了傳統(tǒng)非線性有限元在結(jié)構(gòu)非線性問(wèn)題中遇到的困難，為今后求解結(jié)構(gòu)大變形這類(lèi)問(wèn)題提供了一種新的研究思路和方法。

[1] Cundall P A,Strack ODL.A discrete numerical model for granular assemblies[J].Geotechnique,1979,29(1):47-65.

[2] 徐 泳，孫其誠(chéng)，張 凌.顆粒離散元法研究進(jìn)展[J].力學(xué)進(jìn)展，2003，33(2)：251-260.

[3] 唐志平.激光輻照下充壓柱殼失效的三維離散元模擬[J].爆炸與沖擊，2001，21(1)：1-7.

[4] Itasca Consulting Group Inc.,Particle Flow Code in 2-Dimensions(PFC2D),Version 3.10,Minneapolis,Minnesota,2004.

[5] KJELL MATTIASSON.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.