考慮剪切和扭轉變形的有限質點法纖維梁單元研究

林賢宏，羅堯治，唐敬哲，汪 偉，鄭延豐，楊 超

(浙江大學空間結構研究中心，浙江，杭州 310058)

纖維梁單元模型為結構的彈塑性分析提供了一種高效通用且精度較高的模擬手段，在有限元中得到廣泛的研究與應用。在纖維梁模型中，梁截面被劃分成若干離散的纖維，且假定截面變形前后保持平面狀態，通過獨立分析截面上各纖維的彈塑性狀態，進而實現對梁截面塑性發展過程的精細化模擬。傳統纖維梁模型[1-2]中假定纖維為單軸受力狀態，僅考慮纖維材料的單軸本構關系，能夠較好地模擬梁單元在軸力和彎矩為主作用下的非線性行為，但在模擬以剪切變形為主的非線性行為時效果不好。Ranzo 等[3]和Marini 等[4]通過引入截面剪力與剪應變的非線性關系來考慮剪切變形的影響，模型簡單且易于實現，但無法完全考慮軸力、彎矩和剪力之間的耦合作用。為了更真實地模擬各種外力作用下構件的彈塑性響應，一些學者通過引入纖維材料的多軸本構關系，推導了能考慮剪切變形[5-10]和扭轉變形[11]影響的纖維梁模型。胡鄭州等[12]推導了考慮剪切效應的三維纖維梁大位移增量非線性有限元UL 列式，將纖維梁單元應用到考慮大位移、大轉動的幾何非線性分析中，但尚未考慮材料非線性的影響。目前，纖維梁單元模型在有限元理論中得到了較好的發展與應用[13]，但是基于有限元理論發展的纖維梁單元在求解考慮大位移、大轉動的幾何非線性問題時需要采用TL 或UL 列式進行迭代求解，在求解結構材料非線性問題時需要進行剛度矩陣的修正與迭代，分析較為繁瑣，求解代價較高。

與有限元方法不同，基于向量式力學發展的有限質點法(finite particle method，FPM)采用離散的質點來描述連續體，通過求解質點的運動方程來獲得結構的響應特征，適合于求解復雜的結構非線性問題[14]。在有限質點法中，各質點運動方程的求解過程相互獨立，不涉及結構剛度矩陣的集成和求逆，無需迭代求解，避免了結構剛度矩陣奇異造成的數值求解困難。該方法在機構運動[15-16]、結構動力非線性分析[17-18]、結構倒塌破壞[19]、薄殼屈曲[20]、膜結構形態分析[21-22]及結構多尺度精細化模擬[23]等問題的分析中均取得了較好的效果。目前，Tang 等[24]針對有限質點法提出了基于GPU 并行加速的計算策略，并應用于有限質點法通用計算平臺的研發中，進一步提高了有限質點法的求解效率。

已有學者推導了有限質點法的梁單元，并圍繞其彈塑性行為展開了研究。有限質點法在求解結構的彈塑性行為時無需進行結構剛度矩陣的修正與迭代，只需要在每個途徑單元內按小變形假定根據本構關系求解單元內力，具有求解簡單與計算效率高的優點。喻瑩[19]推導了梁單元的質點內力求解方法，給出了利用塑性鉸法分析結構的彈塑性行為的計算流程，能簡單高效地實現結構宏觀破壞模擬分析。但塑性鉸法無法考慮塑性在截面上及沿桿長方向上的塑性分布與發展，計算精度不高。Duan 等[25]將傳統纖維梁模型引入向量式有限元中，模擬了斜拉橋的倒塌破壞過程。Liu等[26]將考慮扭轉變形和材料二維本構關系的纖維梁單元引入向量式有限元中，模擬了埋地管道在走滑斷層附近的大變形行為。目前有限質點法中纖維梁的研究和應用還較少，現有的纖維梁單元尚無法考慮剪切變形，具有一定的局限性。

本文推導了可考慮剪切和扭轉變形的有限質點法纖維梁單元。從有限質點法的基本原理出發，將考慮材料三維本構關系的纖維梁模型融入有限質點法的計算求解框架中，基于Timoshenko梁理論推導了能夠考慮剪切變形和扭轉變形的纖維梁單元(TM 纖維梁) 的質點內力求解公式，為處理考慮大位移、大轉動的彈塑性問題提供了一種有效的解決方案，相較于塑性鉸法具有更高的精度，相較于傳統有限元法具有更高的計算效率。

1 有限質點法基本方法

1.1 質點運動方程

1.2 質點內力

有限質點法通過利用虛擬逆向運動排除結構剛體位移，得到單元的純變形，將單元內力的求解限制在小變形小轉動的范圍內，從而可以方便地求得質點內力。以空間梁單元AB為例，建立以A為原點的單元局部變形坐標系，x軸的方向為單元的長度方向，y軸和z軸在垂直于x軸的平面內，如圖1 所示?？紤]該單元在時間段ta～tb內從AB運動到A′B′，質點A′和B′經過虛擬逆向運動，分別運動到A′′和B′′，從而可以求得梁單元AB的純變形分量ud，具體的推導過程可參見文獻[19]。單元純變形分量ud可表達為：

圖1 梁單元的純變形Fig. 1 Pure deformation of beam element

由單元的純變形可求得單元在局部變形坐標系下的相應的單元質點內力，經過坐標變換即可得到全局坐標系下的質點內力。有限質點法具有統一的計算求解框架，不同單元類型的質點內力求解方法的實現只需在相應的求解模塊中進行修改，易于在程序中實現。下面將介紹TM 纖維梁單元的質點內力求解方法。

2 基于TM 纖維梁單元模型的質點內力求解

2.1 纖維梁單元模型基本原理

在纖維梁單元模型中，梁截面被離散為一定數量的纖維，幾種典型截面的纖維離散方式如圖2所示。纖維梁單元模型假定截面在單元變形過程中始終保持平面，假設每個纖維截面上的應變均勻分布，可取每根纖維中心點的應變作為其代表值。在得到各個纖維的應變后，結合纖維材料的本構關系可以計算得到截面各個位置的應力狀態，對截面內力進行積分計算得到單元內力。纖維梁單元模型通過獨立求解截面上各纖維的應力狀態實現了對梁單元截面的塑性發展過程的模擬。

圖2 典型截面的纖維離散Fig. 2 Fiber discretization of typical sections

有限質點法中單元的變形與內力的求解在單元局部變形坐標系中完成，下文將以梁長為L的矩形截面纖維梁單元AB(如圖3 所示)為例，介紹基于Timoshenko 梁理論的有限質點法纖維梁單元模型的描述方式，并給出質點內力求解方法。

圖3 纖維梁單元示意圖Fig. 3 Schematic diagram of fiber beam element

2.2 單元位移描述

空間梁單元AB在局部變形坐標系下的位移場可描述為：

2.3 纖維狀態的確定

在通常的空間梁單元理論中，可以只考慮截面上任一點的3 個相互獨立的應力和應變分量[27]：

材料達到屈服后應變增量可分為彈性部分和塑性部分，即：

其他材料的彈塑性本構關系可根據具體情況作相應的推導與修改。

2.4 單元質點內力的求解

單元純變形分量和與之對應的單元質點的基本內力之間的關系由虛功原理確定，它們滿足以下方程：

3 算例分析

本文算法已在有限質點法計算平臺[24]中編制了相應的程序模塊，下面將通過具體算例來驗證其在解決結構彈塑性問題中的正確性和有效性，并以單層網殼結構的動力響應分析為例比較了有限質點法平臺與通用有限元軟件ANSYS 的計算效率。

3.1 空間直角剛架

圖4 所示的空間直角剛架結構是空間桿系結構彈塑性求解的典型算例。結構的幾何和材料特性見圖中參數，材料為理想彈塑性。豎向荷載F作用于BD桿的中點C，結構承受彎矩、扭矩和剪力的共同作用。分別取桿件截面邊長a=0.25 m和a=4 m 這2 種情況進行計算分析。采用有限質點法求解該結構的彈塑性變形，將結構離散為41個質點(40 個單元)，將桿件截面離散為10×10的纖維截面。計算時采用力控制的加載方式，時間步長取為2×10-5s。

圖4 空間直角剛架Fig. 4 Spatial right-angle frame

當a=0.25 m 時，節點B的豎向位移與荷載F之間的無量綱曲線如圖5 所示，并與Park 等[29]和喻瑩[19]的計算結果進行了對比。從圖中可以看出，由于喻瑩在分析時采用基于塑性鉸法的屈服準則，結果呈現出多段線的特征；而本文所采用的有限質點法纖維梁單元能有效地模擬桿件在彎扭作用下的截面塑性的發展過程，與Park 采用有限元塑性區法的分析結果吻合較好，屈服曲線均較為光滑。

圖5 空間直角剛架的荷載-位移曲線Fig. 5 Load-deflection curve of spatial right-angle frame

當a=4 m 時，節點B的豎向位移與荷載F之間的無量綱曲線如圖6 所示，圖中還給出了通用有限元軟件ANSYS 的Beam188 單元的計算結果，并與a=0.25 m 時的結果進行了對比。結果表明：① 由于截面尺寸的增大，剪切變形的影響更為顯著，該結構的無量綱荷載系數有所下降；② 在彈性階段，有限質點法TM 纖維梁單元與Beam188 單元的模擬結果基本一致，可較好地反映剪切變形的影響；③ 在塑性階段，Beam188 單元的模擬結果要略大于有限質點法TM 纖維梁單元的模擬結果，這是由于ANSYS 中的Beam188單元中假設剪力與剪應變為彈性關系，而有限質點法TM 纖維梁單元考慮了正應力與剪應力的耦合作用，并在計算中采用了材料的三維彈塑性本構關系。

圖6 空間直角剛架的荷載-位移曲線對比Fig. 6 Load-deflection curve comparison of spatial right-angle frame

3.2 空間框架

圖7 所示為柱子底端剛接的空間框架結構，其幾何和材料特性見圖中參數，材料取等向強化模型。水平集中荷載H作用于E點，豎向集中荷載V=1.3H作用于桿件DE的中點G。Harrison[30]對該結構進行了試驗研究，Teh 等[31]采用有限元塑性區法對該結構進行了彈塑性計算分析。采用有限質點法求解該結構的彈塑性變形，將每根桿件離散為11 個質點(10 個單元)，將桿件截面離散為32 根纖維，計算時間步長取為5×10-6s。

圖7 空間框架Fig. 7 Space frame

節點A的水平位移與荷載H之間的關系曲線如圖8 所示，并與Harrison[30]和Teh 等[31]的結果進行了對比。從圖中可以看出，本文方法能較好地模擬該結構的彈塑性變形行為，與Teh 采用有限元塑性區法分析的結果基本一致，且接近于Harrison的試驗結果。

圖8 空間框架的荷載-位移曲線Fig. 8 Load-deflection curve of space frame

3.3 六角空間剛架

圖9 所示的六角空間剛架結構是空間桿系結構彈塑性失穩分析的典型算例，包括了幾何大變形和材料非線性因素的影響。該結構的幾何和材料特性見圖中參數，材料為理想彈塑性。豎向集中荷載P作用于結構頂點。

圖9 六角空間剛架 /mFig. 9 Framed dome

采用有限質點法模擬該結構的彈塑性失穩過程，每根桿件被均勻離散為9 個質點(8 個單元)，將桿件截面離散為10×10的纖維截面。計算時采用位移控制的加載方式跟蹤結構失穩后的大變形行為，時間步長取為6×10-5s，單位位移增量Δu=6×10-6m。同時，本文還采用了通用有限元軟件ANSYS 的Beam188 單元進行了模擬計算。

圖10 給出了本文計算得到結構頂點豎向位移與荷載大小之間的關系曲線，并與Park 等[29]、喻瑩[19]以及通用有限元軟件ANSYS 的模擬結果進行了對比。由圖可見，本文提出的有限質點法TM 纖維梁單元能夠有效地模擬該結構的彈塑性失穩過程，與Park 采用有限元塑性區法分析的結果吻合較好，與ANSYS 軟件的計算結果基本吻合，相較于喻瑩采用塑性鉸法分析的結果具有更高的精度。

圖10 六角空間剛架頂點的荷載-位移曲線Fig. 10 Load-deflection curve of framed dome

3.4 單層球面網殼結構

圖11 所示為周邊鉸接的K6 型單層球面網殼結構，其幾何和材料特性見圖中參數，材料為理想彈塑性。節點A作用有如圖12 所示的沖擊荷載P，分別取荷載幅值P0=10 kN和P0=100 kN進行動力響應分析。本文在有限質點法計算平臺FPM 中采用TM 纖維梁單元對桿件進行模擬，將每根桿件離散為6 個單元，共有5580 個單元。將桿件截面均勻離散為16 根纖維，計算時間步長取2×10-5s，不考慮阻尼的作用。作為對比，選取ANSYS 軟件對該結構進行動力響應分析，采用Beam188 單元對桿件進行模擬，將每根桿件劃分為6 個單元，迭代步的最大時間步長取2×10-3s。兩種計算程序均在同一臺計算機上運行，其基本硬件參數如表1 所示。

表1 計算機硬件參數Table 1 Computer hardware parameters

圖11 單層球面網殼結構Fig. 11 Single-layer spherical reticulated shell structure

圖13 和圖14 分別給出了荷載幅值P0=10 kN和荷載幅值P0=100 kN時，節點A 的豎向位移時程曲線。由圖可知，本文算法與ANSYS 的模擬結果吻合較好。ANSYS 計算時采用的是基于傳統有限元的隱式動力分析，處理非線性問題時需要進行迭代求解。而FPM 采用中心差分法求解運動方程，無需進行迭代，具有較高的計算效率。兩種程序計算所耗費的時間如表2 所示。當P0=10 kN時，FPM 的計算耗時僅為ANSYS 的19%。當P0=10 kN時，結構的非線性特征更為顯著，ANSYS需要耗費更多的時間用于迭代計算以保證結果的收斂性，FPM 的計算耗時僅為ANSYS 的10%，具有較大的優勢。

圖13 節點A 的豎向位移時程曲線(P0=10 kN)Fig. 13 Vertical displacement time history curve of node A (P0=10 kN)

圖14 節點A 的豎向位移時程曲線(P0=100 kN)Fig. 14 Vertical displacement time history curve of node A (P0=100 kN)

表2 FPM 與ANSYS 的計算效率對比Table 2 Calculation efficiency comparison between FPM and ANSYS

4 結論

本文將纖維梁模型基本原理與有限質點法相結合，并考慮纖維材料的三維本構關系，提出了考慮剪切和扭轉變形的有限質點法纖維梁單元模型，基于Timoshenko 梁理論推導了相應的單元質點內力求解公式，為處理考慮大位移、大轉動的彈塑性問題提供了一種有效的解決方案。將纖維梁模型與有限質點法相結合，可以充分發揮兩者的優勢，在處理幾何和材料雙非線性問題中具有求解簡單、易于實現的特點。數值算例分析表明，本文提出的纖維梁單元能較好地模擬結構彈塑性變形及結構大變形失穩等問題，且能夠考慮桿件截面的塑性發展過程，相較于FPM 塑性鉸法具有更高的精度，相較于傳統有限元法具有更高的計算效率。