葉繼紅，許 強(qiáng)

(1.中國(guó)礦業(yè)大學(xué)江蘇省土木工程環(huán)境災(zāi)變與結(jié)構(gòu)可靠性重點(diǎn)實(shí)驗(yàn)室，江蘇，徐州 221116；2.中國(guó)礦業(yè)大學(xué)徐州市工程結(jié)構(gòu)火災(zāi)安全重點(diǎn)實(shí)驗(yàn)室，江蘇，徐州 221116)

有限元分析方法中，結(jié)構(gòu)的屈曲現(xiàn)象一般作為幾何非線性問題處理[1-3]，其荷載-位移曲線具有“位移跳躍”(snap-through)或“荷載跌落”(snap-back)顯著特征，傳統(tǒng)的分級(jí)加載模式無法順利越過臨界點(diǎn)。為此，RIKS[4]和CRISFIELD[5]提出了隱式弧長(zhǎng)法，通過引入弧長(zhǎng)約束方程實(shí)現(xiàn)“自適應(yīng)”加載目的，成功捕捉結(jié)構(gòu)完整的屈曲路徑。HUANG 等[6]采用隱式弧長(zhǎng)法捕捉了雙曲復(fù)合殼體的卡扣現(xiàn)象。PAN 和LIANG[7]基于隱式弧長(zhǎng)法研究了風(fēng)梁等加固裝置對(duì)儲(chǔ)罐薄壁結(jié)構(gòu)屈曲行為的影響。但在隱式算法中需要不斷地集成和修正切線剛度矩陣，在臨界點(diǎn)附近剛度矩陣出現(xiàn)病態(tài)甚至奇異。RAMESH 等[8-10]將柱面弧長(zhǎng)法與動(dòng)力松弛法(DRM)結(jié)合，提出了一種顯式弧長(zhǎng)法，由于不涉及剛度矩陣的集成、求解，更具穩(wěn)定性，得到了廣泛的應(yīng)用與發(fā)展。ZU 等[11-12]基于顯式弧長(zhǎng)法對(duì)索桿張力結(jié)構(gòu)成型過程進(jìn)行了數(shù)值模擬。張鵬飛等[13]將有限質(zhì)點(diǎn)法與顯式弧長(zhǎng)法結(jié)合，對(duì)薄殼屈曲問題進(jìn)行了研究。ZHANG 等[14]采用顯式弧長(zhǎng)法對(duì)殼與桿系結(jié)構(gòu)進(jìn)行了屈曲分析。LEE 等[15-16]在顯式弧長(zhǎng)法中引入動(dòng)力學(xué)阻尼，并采用對(duì)角剛度替換控制方程中的真實(shí)質(zhì)量，簡(jiǎn)化了計(jì)算過程，提高了算法的穩(wěn)定性與收斂速度。

與有限元法不同，傳統(tǒng)離散元法是一種非連續(xù)介質(zhì)力學(xué)方法，最初主要應(yīng)用于研究巖石等非連續(xù)介質(zhì)的力學(xué)行為，因其適合處理大變形、大位移等強(qiáng)非線性問題，所以被逐步引入到結(jié)構(gòu)工程領(lǐng)域?；诖?，葉繼紅教授課題組率先提出桿系離散元法(member discrete element method，MDEM)[17]，并成功將其應(yīng)用于連續(xù)介質(zhì)問題。該法將桿件離散為有限的剛性球顆粒，相鄰顆粒間通過接觸本構(gòu)模型聯(lián)接，然后采用牛頓第二定律確定各顆粒的空間坐標(biāo)，并以顆粒的空間位置描述結(jié)構(gòu)的空間位置和幾何構(gòu)型。桿系離散元不僅擴(kuò)大了傳統(tǒng)離散元的適用范圍，同時(shí)保持了傳統(tǒng)離散元法無需組集整體剛度矩陣和迭代求解運(yùn)動(dòng)方程的優(yōu)勢(shì)。目前已成功用于桿系結(jié)構(gòu)的幾何非線性行為[18-19]、材料非線性行為[20-21]、半剛性連接[22-23]、非線性動(dòng)力響應(yīng)[24-26]、屈曲行為[27]、網(wǎng)殼倒塌破壞模擬中[28-31]。上述研究成果表明，桿系離散元法在模擬結(jié)構(gòu)非線性以及斷裂倒塌等復(fù)雜行為中具有獨(dú)特優(yōu)勢(shì)，是一種簡(jiǎn)單而有效的方法。

本文在桿系離散元理論的基礎(chǔ)上，建立了一種顯式弧長(zhǎng)法以追蹤結(jié)構(gòu)屈曲平衡路徑。該方法無需組集剛度矩陣，求解時(shí)不需要刻意區(qū)分結(jié)構(gòu)是小變形還是大變形行為，與傳統(tǒng)有限元方法相比更具優(yōu)越性。首先介紹了桿系離散元的基本理論；其次，建立了新型接觸單元—鉸固接觸的本構(gòu)模型，分別采用離散元法和結(jié)構(gòu)力學(xué)方法對(duì)鉸固接觸單元在已知端部位移工況下進(jìn)行受力分析，然后通過端部?jī)?nèi)力相等建立方程，確定了接觸彈簧剛度系數(shù)；再者，利用動(dòng)力阻尼技術(shù)獲取結(jié)構(gòu)的準(zhǔn)靜態(tài)阻尼運(yùn)動(dòng)，大大簡(jiǎn)化了顆粒運(yùn)動(dòng)方程的求解過程，然后引入總位移約束弧長(zhǎng)法，并詳細(xì)介紹了與離散元法相結(jié)合的求解策略和實(shí)施過程；最后通過算例驗(yàn)證了該方法對(duì)結(jié)構(gòu)屈曲分析的準(zhǔn)確性和適用性。

1 桿系離散元法基本方程及內(nèi)力求解

1.1 顆粒運(yùn)動(dòng)方程

桿系離散元法將結(jié)構(gòu)離散為若干球顆粒的集合，如圖1(a)所示，相鄰顆粒間(如顆粒A 和B)則通過虛擬的零長(zhǎng)度彈簧系統(tǒng)進(jìn)行聯(lián)接，并形成離散元法的基本分析單元-接觸單元，如圖1(b)所示。在任意時(shí)刻t，桿系離散元模型中所有顆粒的運(yùn)動(dòng)均滿足牛頓第二定律，每個(gè)顆粒的運(yùn)動(dòng)控制方程具體形式為：

圖1 桿系離散元分析模型Fig.1 Member discrete element analysis model

式中：at、vt分別為t時(shí)刻顆粒的加速度和速度矢量，包含3 個(gè)平動(dòng)和3 個(gè)轉(zhuǎn)動(dòng)；M和C分別為顆粒的等效質(zhì)量矩陣和阻尼矩陣；Pt和Ft分別為顆粒質(zhì)心所受外力和內(nèi)力，其中內(nèi)力Ft由接觸內(nèi)力通過力系平移定理計(jì)算得到。

假定在每時(shí)間步Δt內(nèi)速度是線性變化的，由中心差分法，可得t時(shí)刻顆粒的速度和加速度為：

將式(2)代入式(1)，得到顆粒在t+Δt/2 時(shí)刻的速度為：

式中，Rt=Pt-Ft為不平衡力。

1.2 接觸單元的內(nèi)力計(jì)算

要獲取圖1 中各顆粒的內(nèi)力，需先求出相鄰顆粒間接觸點(diǎn)(如接觸點(diǎn)C)的內(nèi)力，然后通過力平移定理，反向疊加到相鄰顆粒(如顆粒A 和B)上。對(duì)于線彈性體系，接觸點(diǎn)C處的內(nèi)力增量可由接觸本構(gòu)關(guān)系計(jì)算[22]，其表達(dá)式為：

式中：Δf為接觸點(diǎn)處虛擬零長(zhǎng)度彈簧系統(tǒng)的內(nèi)力增量；K為接觸單元的彈性剛度矩陣；Δd為局部坐標(biāo)系下虛擬零長(zhǎng)度彈簧系統(tǒng)兩端結(jié)點(diǎn)的位移增量。

2 接觸本構(gòu)模型

離散元的核心問題是式(4)中接觸剛度系數(shù)的求解。文獻(xiàn)[17]基于簡(jiǎn)單梁理論，通過應(yīng)變能等效推導(dǎo)出適用于“梁?jiǎn)卧钡慕佑|剛度計(jì)算公式。但桁架結(jié)構(gòu)較為特殊，如圖2 所示，桿件之間為鉸接，在節(jié)點(diǎn)處形成的接觸單元(以下簡(jiǎn)稱鉸固接觸單元)，其特征為一顆粒固接、一顆粒鉸接，與兩顆粒固接的梁接觸單元運(yùn)動(dòng)特性具有顯著區(qū)別。

圖2 二桿桁架結(jié)構(gòu)及其離散元分析模型Fig.2 Two-bar truss structure and DEM analysis model

為確定鉸固接觸單元?jiǎng)偠认禂?shù)，本文分別采用離散元法和結(jié)構(gòu)力學(xué)方法對(duì)鉸固接觸單元在已知端部位移工況下進(jìn)行受力分析，并通過端部?jī)?nèi)力相等建立方程，直接建立非連續(xù)介質(zhì)方法與連續(xù)介質(zhì)問題之間的聯(lián)系。

因鉸固接觸單元自身運(yùn)動(dòng)的特殊性，對(duì)其提出以下三點(diǎn)假定：

1) 鉸端顆粒上的接觸點(diǎn)及兩顆粒球心三點(diǎn)共線(運(yùn)動(dòng)假定)。

在離散元法中，接觸彈簧位移由接觸點(diǎn)對(duì)描述，其值由兩顆粒相對(duì)運(yùn)動(dòng)關(guān)系確定。如圖3 所示，初始時(shí)，接觸點(diǎn)3、4 重合并分別在兩顆粒邊上，當(dāng)兩顆粒運(yùn)動(dòng)時(shí)，接觸點(diǎn)對(duì)(3、4)發(fā)生相對(duì)位移，考慮鉸節(jié)點(diǎn)處的顆粒1 可自由轉(zhuǎn)動(dòng)，本文在確定接觸點(diǎn)對(duì)相對(duì)位移大小時(shí)，假定鉸接顆粒上的接觸點(diǎn)(3)和兩顆粒球心(1、2)三點(diǎn)共線。

圖3 鉸固單元的受力分析Fig.3 Force analysis of hinged-fixed elements

2) 轉(zhuǎn)動(dòng)剛度為0(鉸節(jié)點(diǎn)可自由轉(zhuǎn)動(dòng))。

3) 利用力系平移定理將接觸力移至兩顆粒質(zhì)心時(shí)，在鉸接顆粒上產(chǎn)生的附加彎矩由固端顆粒承擔(dān)(鉸固單元整體受力平衡)。

在離散元法中，單元的力平衡條件被轉(zhuǎn)化為兩剛體接觸處作用力與反作用力的相互作用，將接觸力移至顆粒質(zhì)心時(shí)，會(huì)在質(zhì)心處產(chǎn)生附加彎矩，由于鉸接處的顆粒不能承受彎矩作用，因此，從單元整體受力平衡的角度考慮，將鉸接顆粒質(zhì)心處產(chǎn)生的附加彎矩移至固端，由固端顆粒承擔(dān)。

由于彈簧系統(tǒng)中的彈簧相互獨(dú)立，互不耦合，在確定切向剛度系數(shù)kτ時(shí)，為直觀起見，基于上述3 點(diǎn)假定，對(duì)圖3 平面鉸固單元在已知端部位移工況下進(jìn)行受力分析。

對(duì)于圖3(a)中的離散元模型，采用離散元法計(jì)算的端部?jī)?nèi)力為：

式中：M11、M22分別為接觸力在顆粒1 和顆粒2 球心處產(chǎn)生的等效彎矩。

對(duì)于圖3(a)中的有限元模型，采用結(jié)構(gòu)力學(xué)方法[32]計(jì)算的端部?jī)?nèi)力為：

令式(5)與式(6)計(jì)算的端部?jī)?nèi)力相等建立方程，即：

求得切向剛度為：

同理，對(duì)于圖3(b)，可求得切向剛度kτ=6EI/L3。

由于鉸固接觸單元拉壓行為與結(jié)構(gòu)力學(xué)中桿件的拉壓行為一致，因此其軸向剛度系數(shù)與結(jié)構(gòu)力學(xué)中桿件軸向剛度相同，即kn=EA/L。

引申至三維空間時(shí)，鉸固接觸單元的剛度系數(shù)矩陣為：

值得注意的是，采用上述方法求得梁接觸單元?jiǎng)偠认禂?shù)與文獻(xiàn)[17]相同。

3 基于離散元理論的顯式弧長(zhǎng)法

3.1 動(dòng)力阻尼技術(shù)

當(dāng)不考慮粘滯阻尼時(shí)，式(3)簡(jiǎn)化為：

此時(shí)分析過程簡(jiǎn)單且變量少，但由于沒有了阻尼的抑制作用，系統(tǒng)會(huì)發(fā)生簡(jiǎn)諧振動(dòng)，為此，引入動(dòng)力阻尼技術(shù)[15]，即每時(shí)步計(jì)算后，由下式計(jì)算結(jié)構(gòu)的動(dòng)能：

式中，N為總顆粒數(shù)。

對(duì)于靜力問題，在分析過程中，結(jié)構(gòu)產(chǎn)生的動(dòng)能大部分會(huì)轉(zhuǎn)化為勢(shì)能，其動(dòng)能會(huì)減小，若此時(shí)的動(dòng)能小于上一時(shí)刻計(jì)算的動(dòng)能，則將所有顆粒的速度矢量設(shè)置為0 向量，然后進(jìn)行下時(shí)步分析。

利用式(10)可計(jì)算出Δt時(shí)間步內(nèi)各顆粒的位移增量：

在式(12)中，結(jié)構(gòu)所受的外力Pt可由一個(gè)未知荷載參數(shù)λt+Δt與施加的靜荷載P0的乘積表示，即：

3.2 總位移約束顯式弧長(zhǎng)法

在t+Δt時(shí)刻，每個(gè)顆粒的位移可表示為：

將式(13)代入式(14)，整理得：

式中：

為確定式(13)中荷載參數(shù)λt+Δt值，引入總位移弧長(zhǎng)約束方程[5]：

式中，ln+1為指定的總弧長(zhǎng)值。

將式(15)代入式(17)中，整理可得：

式中：

求解式(18)得荷載參數(shù)：

當(dāng)式(20)中的根式大于0 時(shí)，方程有兩個(gè)實(shí)根分別為λ1、λ2，與之對(duì)應(yīng)的位移為dt+Δt,1、dt+Δt,2，分別求出此時(shí)兩位移向量與上一時(shí)刻位移向量的“夾角”：λt+Δt的值由以下條件式確定[8]：

式中：Fj和P0,j分別為加載顆粒j在加載方向上所受內(nèi)力和靜荷載。

3.3 顯式弧長(zhǎng)法的具體計(jì)算流程

圖4 給出了總位移約束弧長(zhǎng)法的基本思想和收斂過程[15]。從圖中可以看出，一個(gè)弧長(zhǎng)步計(jì)算包含2 個(gè)計(jì)算過程：預(yù)測(cè)過程和校正過程[8]。

圖4 總位移約束顯式弧長(zhǎng)法示意Fig.4 Explicit arc-length method for total displacement constraint

預(yù)測(cè)過程是一個(gè)力加載過程：首先確定每時(shí)步荷載增量detP；然后進(jìn)行η 次離散元循環(huán)計(jì)算，并通過式(17)計(jì)算第η 次循環(huán)后的弧長(zhǎng)ln+1。

離散元法其本質(zhì)是動(dòng)力分析，對(duì)于靜力問題，可以采用阻尼耗能和緩慢施加外荷載方式逼近結(jié)構(gòu)靜力解，因此每時(shí)步荷載增量detP取值不宜過大，否則可能導(dǎo)致算法失效。若獲得結(jié)構(gòu)準(zhǔn)靜力解所需加載時(shí)間為Ttol(可取1 s～100 s)，則每時(shí)步荷載增量可取式(23)計(jì)算值或更小值。參數(shù)η 決定柱面弧長(zhǎng)的校正頻率，一般可取10～100，當(dāng)detP取值較大時(shí)，結(jié)構(gòu)響應(yīng)可能發(fā)生較小波動(dòng)，此時(shí)η 可取較小值，增加校正次數(shù)，以獲得穩(wěn)定解。

式中，P0為施加的靜荷載。

校正過程：引入弧長(zhǎng)約束方程式(17)，且在整個(gè)校正過程中弧長(zhǎng)值ln+1保持不變。由式(18)～式(22)確定新的荷載系數(shù) λ′，并在該荷載系數(shù)下進(jìn)行離散元分析，當(dāng)不平衡力R的范數(shù)小于容許值χ(其取值比detP小2～3 數(shù)量級(jí))時(shí)，結(jié)束校正過程，然后進(jìn)入下一弧長(zhǎng)步計(jì)算；否則，再由弧長(zhǎng)約束方程確定新的荷載系數(shù)λt+Δt，繼續(xù)循環(huán)迭代，直至滿足收斂要求。

圖5 為基于離散元理論的顯式弧長(zhǎng)法的計(jì)算流程。

圖5 基于離散元理論的顯式弧長(zhǎng)法計(jì)算流程Fig.5 Calculation flow of explicit arc-length method based on discrete element theory

4 算例分析

4.1 二桿平面桁架結(jié)構(gòu)

圖6 為一桁架結(jié)構(gòu)，其頂端C點(diǎn)受豎向集中荷載P作用。該桁架的屈曲分析是經(jīng)典的跳躍問題，集中外荷載P與頂點(diǎn)位移x的理論關(guān)系為[33]：

圖6 二桿平面桁架結(jié)構(gòu)及DEM 模型Fig.6 Two-bar truss structure and DEM analysis model

式中：彈性模量E=200 GPa；桿件橫截面積A=10 mm2；慣性矩I=745.18 mm4；θ=6°；x=V/L0，V為頂點(diǎn)C的豎向位移，L0=100 mm。

ANSYS 分析時(shí)，對(duì)于桁架結(jié)構(gòu)采用Link8 單元，單元?jiǎng)澐植捎靡粭U一單元；對(duì)于剛架結(jié)構(gòu)采用Beam4 梁?jiǎn)卧?，不考慮桿件的剪切變形，單元?jiǎng)澐植捎靡粭U三單元形式，即將網(wǎng)殼結(jié)構(gòu)每根桿件劃分為3 個(gè)單元，并采用弧長(zhǎng)法進(jìn)行結(jié)構(gòu)失穩(wěn)全過程跟蹤。后文算例亦按此方法進(jìn)行ANSYS 有限元分析。建立DEM 分析模型時(shí)，將該結(jié)構(gòu)離散為13 個(gè)相同的顆粒球單元，接觸(單元)數(shù)量為12，計(jì)算時(shí)步取Δt=10-3s。弧長(zhǎng)法分析參數(shù)為，加載點(diǎn)C處的y向靜荷載P0=2500 N，每時(shí)步荷載增量detP=1 N，迭代次數(shù)η=10，循環(huán)最大容許值χ=1×10-2N。圖7 為C點(diǎn)荷載-位移曲線，可以看出，本文方法計(jì)算結(jié)果與理論解及ANSYS 分析結(jié)果基本一致，曲線幾乎完全吻合，說明本文給出的鉸固單元本構(gòu)模型合理性，同時(shí)也驗(yàn)證了本文基于離散元理論建立的顯式弧長(zhǎng)法的正確性。

圖7 C 點(diǎn)荷載-y 方向位移曲線Fig.7 C point load-y direction displacement curve

在離散元顯式弧長(zhǎng)法中，參數(shù)detP的取值對(duì)算法收斂起決定作用，圖8 為不同detP取值下C點(diǎn)的荷載-位移曲線，可以看出，在式(23)確定的估算值范圍內(nèi)，算法均可以得到正確結(jié)果。相比有限元分析軟件中的弧長(zhǎng)法，本文方法參數(shù)少、穩(wěn)定性好，同時(shí)參數(shù)取值容易，無需使用者有足夠經(jīng)驗(yàn)，即可完成對(duì)結(jié)構(gòu)的屈曲分析。但由于離散法中顆粒較多，分析步長(zhǎng)Δt通常取值很小，其計(jì)算成本相對(duì)較高，且detP取值越小，計(jì)算時(shí)間越長(zhǎng)，見表1。

表1 不同分析工況下計(jì)算時(shí)間Table 1 Calculation time under different cases

圖8 不同分析工況下C 點(diǎn)荷載-y 方向位移曲線Fig.8 Load-y direction displacement curve of point C under different analysis cases

4.2 24 桿星型網(wǎng)穹結(jié)構(gòu)

圖9 為空間六角形星型穹頂結(jié)構(gòu)，該結(jié)構(gòu)已被用作桁架模型非線性分析的驗(yàn)證示例。星型穹頂在頂點(diǎn)中心節(jié)點(diǎn)1 處施加z向集中荷載P，構(gòu)件截面均為正方形且尺寸相同，截面積A＝317 mm2，截面慣性Iy=Iz=8370 mm4，扭轉(zhuǎn)慣性矩Ix=14 110 mm4，彈性模量E=3030 MPa，剪切模量G=1262 MPa，材料為線彈性，六個(gè)支座均為鉸接。

圖9 24 星型穹頂結(jié)構(gòu) /mmFig.9 24-bar star dome structure

該網(wǎng)穹結(jié)構(gòu)共有13 個(gè)節(jié)點(diǎn)，24 根桿件，建立DEM 分析模型，節(jié)點(diǎn)球徑取25 mm，對(duì)桿件剩余部分以接近節(jié)點(diǎn)球徑為原則進(jìn)行均等離散，如圖10 所示，結(jié)構(gòu)最終被離散為121 個(gè)球元，共形成132 個(gè)接觸單元，計(jì)算時(shí)步取Δt=10-3s?；￠L(zhǎng)法分析參數(shù)為，節(jié)點(diǎn)1 處的靜荷載P0=900 N，detP=0.5 N，η=10，χ=1×10-3N。

圖10 24 星型穹頂結(jié)構(gòu)離散元模型Fig.10 DEM analysis model of 24-star dome structure

圖11、圖12 分別為記錄點(diǎn)的荷載-位移曲線，可以看出，本文方法計(jì)算結(jié)果與文獻(xiàn)[34]及ANSYS 分析結(jié)果基本相同，曲線幾乎完全吻合，說明本文給出的鉸固單元本構(gòu)模型和基于離散元理論建立的顯式弧長(zhǎng)法的正確性，并能夠有效跟蹤結(jié)構(gòu)彈性屈曲全過程。

圖11 節(jié)點(diǎn)1 的荷載-位移曲線Fig.11 Load-displacement curve of node 1

圖12 節(jié)點(diǎn)2 荷載-位移曲線Fig.12 Load-displacement curve of node 2

圖13 為不同detP取值下節(jié)點(diǎn)1 的荷載-位移曲線，可以看出，曲線幾乎完全重合，說明采用式(23)計(jì)算的detP估算值，算法是穩(wěn)定的。表2為不同分析工況下的計(jì)算時(shí)間，detP取值越小，耗時(shí)越長(zhǎng)，此時(shí)可增大η 取值，以降低柱面弧長(zhǎng)校正頻率，減少分析時(shí)間。

圖13 不同工況下節(jié)點(diǎn)1 荷載-位移曲線Fig.13 Load- displacement curve of node 1 under different analysis cases

4.3 K6-4 單層穹頂網(wǎng)殼結(jié)構(gòu)

如圖14 所示的K6 型單層穹頂網(wǎng)殼結(jié)構(gòu)，跨度L為30 m，矢高為2 m，采用?180×5 鋼管，材料彈性模量E=210 GPa，剪切模量G=85 GPa，網(wǎng)殼僅在節(jié)點(diǎn)1 作用豎向集中荷載P，支座采用周邊鉸接。建立離散元模型時(shí)將各桿件劃分為7 個(gè)球顆粒，結(jié)構(gòu)共劃分為841 個(gè)球顆粒，如圖15 所示，接觸(單元)數(shù)量為936，計(jì)算時(shí)步取Δt=10-4s。弧長(zhǎng)法分析參數(shù)為，節(jié)點(diǎn)1 處靜荷載P0=1000 kN，detP=0.1 kN，η=10，χ=1×10-3kN。

圖14 K6 型單層穹頂網(wǎng)殼結(jié)構(gòu) /mFig.14 K6 single dome structure

圖15 K6 型單層穹頂網(wǎng)殼結(jié)構(gòu)離散元分析模型Fig.15 DEM analysis of K6 single dome structure

圖16 為節(jié)點(diǎn)1～節(jié)點(diǎn)3 的荷載-位移曲線，可以發(fā)現(xiàn)，結(jié)構(gòu)在第一次屈曲前，本文方法得到的曲線與文獻(xiàn)[27]及ANSYS 分析方法得到曲線基本吻合。在第一次屈曲后，本文方法和文獻(xiàn)[27]基于離散元理論采用位移控制法得到的曲線整體趨勢(shì)相近，僅在數(shù)值上存在較小的差異；但與ANSYS 分析結(jié)果相比，曲線無論在走勢(shì)和數(shù)值上都存在較大差異，這可能是因?yàn)殡x散元法是一種非連續(xù)介質(zhì)力學(xué)方法，有限元法是連續(xù)介質(zhì)力學(xué)方法，兩種分析方法的理論基礎(chǔ)和算法模式不同，計(jì)算結(jié)果會(huì)存在一定的差異，當(dāng)結(jié)構(gòu)發(fā)生屈曲后，結(jié)構(gòu)處于不穩(wěn)定狀態(tài)，這種差異性可能更為突顯。

圖16 節(jié)點(diǎn)1～節(jié)點(diǎn)3 的荷載-z 方向位移曲線Fig.16 Load-direction displacement curves of nodes from 1 to 3

圖17 為不同detP取值下節(jié)點(diǎn)1 的荷載-位移曲線，可以看出，三種分析工況得到的曲線基本吻合，說明在式(23)確定的detP估算值下，算法是穩(wěn)定的。detP較小的取值保證了算法的穩(wěn)定性，但需要更多的計(jì)算時(shí)間，見表3。

圖17 不同工況下節(jié)點(diǎn)1 荷載-z 向位移曲線Fig.17 Load-z direction displacement curve of node 1 under different analysis cases

5 結(jié)論

本文基于桿系離散元理論，提出了一種顯式弧長(zhǎng)法，用于追蹤結(jié)構(gòu)屈曲全過程曲線。對(duì)典型算例和單層網(wǎng)殼結(jié)構(gòu)的彈性失穩(wěn)問題進(jìn)行了研究，得到主要結(jié)論如下：

(1) 建立了新型接觸單元—鉸固接觸單元的本構(gòu)模型，分別采用離散元法和結(jié)構(gòu)力學(xué)方法對(duì)鉸固接觸單元在已知端部位移工況下進(jìn)行受力分析，并通過端部?jī)?nèi)力相等建立方程，推導(dǎo)了接觸剛度系數(shù)，豐富了離散元法的接觸單元類型，擴(kuò)大了其應(yīng)用范圍。

(2) 在離散元顯式弧長(zhǎng)法中，每時(shí)步荷載增量detP的取值對(duì)算法收斂至關(guān)重要，采用本文方法確定的detP估算值，可保證算法穩(wěn)定性，其取值越小，算法越穩(wěn)定；參數(shù)η 控制柱面弧長(zhǎng)約束的校正頻率，當(dāng)detP較小時(shí)，其取值可稍大些。相比傳統(tǒng)有限元弧長(zhǎng)法，該算法參數(shù)少、穩(wěn)定性好，且參數(shù)取值容易，對(duì)分析者自身能力與經(jīng)驗(yàn)要求低，但計(jì)算成本相對(duì)較高。

(3) 在求解顆粒運(yùn)動(dòng)方程時(shí)，采用動(dòng)力阻尼技術(shù)獲取結(jié)構(gòu)的準(zhǔn)靜態(tài)解，因無粘滯阻尼項(xiàng)，大大簡(jiǎn)化了計(jì)算過程，并有效地與總位移約束弧長(zhǎng)法銜接。同時(shí)，本文方法無需組集單元?jiǎng)偠染仃嚭偷蠼夥匠探M，就可以越過結(jié)構(gòu)屈曲的臨界點(diǎn)，與傳統(tǒng)有限元方法相比更具優(yōu)越性，為結(jié)構(gòu)失穩(wěn)分析提供了新的算法工具。