基于桿系離散元理論的結構屈曲顯式弧長法研究

葉繼紅，許 強

(1.中國礦業大學江蘇省土木工程環境災變與結構可靠性重點實驗室，江蘇，徐州 221116；2.中國礦業大學徐州市工程結構火災安全重點實驗室，江蘇，徐州 221116)

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

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

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

1 桿系離散元法基本方程及內力求解

1.1 顆粒運動方程

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

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

式中：at、vt分別為t時刻顆粒的加速度和速度矢量，包含3 個平動和3 個轉動；M和C分別為顆粒的等效質量矩陣和阻尼矩陣；Pt和Ft分別為顆粒質心所受外力和內力，其中內力Ft由接觸內力通過力系平移定理計算得到。

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

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

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

1.2 接觸單元的內力計算

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

式中：Δf為接觸點處虛擬零長度彈簧系統的內力增量；K為接觸單元的彈性剛度矩陣；Δd為局部坐標系下虛擬零長度彈簧系統兩端結點的位移增量。

2 接觸本構模型

離散元的核心問題是式(4)中接觸剛度系數的求解。文獻[17]基于簡單梁理論，通過應變能等效推導出適用于“梁單元”的接觸剛度計算公式。但桁架結構較為特殊，如圖2 所示，桿件之間為鉸接，在節點處形成的接觸單元(以下簡稱鉸固接觸單元)，其特征為一顆粒固接、一顆粒鉸接，與兩顆粒固接的梁接觸單元運動特性具有顯著區別。

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

為確定鉸固接觸單元剛度系數，本文分別采用離散元法和結構力學方法對鉸固接觸單元在已知端部位移工況下進行受力分析，并通過端部內力相等建立方程，直接建立非連續介質方法與連續介質問題之間的聯系。

因鉸固接觸單元自身運動的特殊性，對其提出以下三點假定：

1) 鉸端顆粒上的接觸點及兩顆粒球心三點共線(運動假定)。

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

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

2) 轉動剛度為0(鉸節點可自由轉動)。

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

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

由于彈簧系統中的彈簧相互獨立，互不耦合，在確定切向剛度系數kτ時，為直觀起見，基于上述3 點假定，對圖3 平面鉸固單元在已知端部位移工況下進行受力分析。

對于圖3(a)中的離散元模型，采用離散元法計算的端部內力為：

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

對于圖3(a)中的有限元模型，采用結構力學方法[32]計算的端部內力為：

令式(5)與式(6)計算的端部內力相等建立方程，即：

求得切向剛度為：

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

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

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

值得注意的是，采用上述方法求得梁接觸單元剛度系數與文獻[17]相同。

3 基于離散元理論的顯式弧長法

3.1 動力阻尼技術

當不考慮粘滯阻尼時，式(3)簡化為：

此時分析過程簡單且變量少，但由于沒有了阻尼的抑制作用，系統會發生簡諧振動，為此，引入動力阻尼技術[15]，即每時步計算后，由下式計算結構的動能：

式中，N為總顆粒數。

對于靜力問題，在分析過程中，結構產生的動能大部分會轉化為勢能，其動能會減小，若此時的動能小于上一時刻計算的動能，則將所有顆粒的速度矢量設置為0 向量，然后進行下時步分析。

利用式(10)可計算出Δt時間步內各顆粒的位移增量：

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

3.2 總位移約束顯式弧長法

在t+Δt時刻，每個顆粒的位移可表示為：

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

式中：

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

式中，ln+1為指定的總弧長值。

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

式中：

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

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

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

3.3 顯式弧長法的具體計算流程

圖4 給出了總位移約束弧長法的基本思想和收斂過程[15]。從圖中可以看出，一個弧長步計算包含2 個計算過程：預測過程和校正過程[8]。

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

預測過程是一個力加載過程：首先確定每時步荷載增量detP；然后進行η 次離散元循環計算，并通過式(17)計算第η 次循環后的弧長ln+1。

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

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

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

圖5 為基于離散元理論的顯式弧長法的計算流程。

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

4 算例分析

4.1 二桿平面桁架結構

圖6 為一桁架結構，其頂端C點受豎向集中荷載P作用。該桁架的屈曲分析是經典的跳躍問題，集中外荷載P與頂點位移x的理論關系為[33]：

圖6 二桿平面桁架結構及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為頂點C的豎向位移，L0=100 mm。

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

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

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

表1 不同分析工況下計算時間Table 1 Calculation time under different cases

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

4.2 24 桿星型網穹結構

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

圖9 24 星型穹頂結構 /mmFig.9 24-bar star dome structure

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

圖10 24 星型穹頂結構離散元模型Fig.10 DEM analysis model of 24-star dome structure

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

圖11 節點1 的荷載-位移曲線Fig.11 Load-displacement curve of node 1

圖12 節點2 荷載-位移曲線Fig.12 Load-displacement curve of node 2

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

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

4.3 K6-4 單層穹頂網殼結構

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

圖14 K6 型單層穹頂網殼結構 /mFig.14 K6 single dome structure

圖15 K6 型單層穹頂網殼結構離散元分析模型Fig.15 DEM analysis of K6 single dome structure

圖16 為節點1～節點3 的荷載-位移曲線，可以發現，結構在第一次屈曲前，本文方法得到的曲線與文獻[27]及ANSYS 分析方法得到曲線基本吻合。在第一次屈曲后，本文方法和文獻[27]基于離散元理論采用位移控制法得到的曲線整體趨勢相近，僅在數值上存在較小的差異；但與ANSYS 分析結果相比，曲線無論在走勢和數值上都存在較大差異，這可能是因為離散元法是一種非連續介質力學方法，有限元法是連續介質力學方法，兩種分析方法的理論基礎和算法模式不同，計算結果會存在一定的差異，當結構發生屈曲后，結構處于不穩定狀態，這種差異性可能更為突顯。

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

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

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

5 結論

本文基于桿系離散元理論，提出了一種顯式弧長法，用于追蹤結構屈曲全過程曲線。對典型算例和單層網殼結構的彈性失穩問題進行了研究，得到主要結論如下：

(1) 建立了新型接觸單元—鉸固接觸單元的本構模型，分別采用離散元法和結構力學方法對鉸固接觸單元在已知端部位移工況下進行受力分析，并通過端部內力相等建立方程，推導了接觸剛度系數，豐富了離散元法的接觸單元類型，擴大了其應用范圍。

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

(3) 在求解顆粒運動方程時，采用動力阻尼技術獲取結構的準靜態解，因無粘滯阻尼項，大大簡化了計算過程，并有效地與總位移約束弧長法銜接。同時，本文方法無需組集單元剛度矩陣和迭代求解方程組，就可以越過結構屈曲的臨界點，與傳統有限元方法相比更具優越性，為結構失穩分析提供了新的算法工具。