考慮剪切和扭轉(zhuǎn)變形的有限質(zhì)點(diǎn)法纖維梁單元研究

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

(浙江大學(xué)空間結(jié)構(gòu)研究中心，浙江，杭州 310058)

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

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

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

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

1 有限質(zhì)點(diǎn)法基本方法

1.1 質(zhì)點(diǎn)運(yùn)動方程

1.2 質(zhì)點(diǎn)內(nèi)力

有限質(zhì)點(diǎn)法通過利用虛擬逆向運(yùn)動排除結(jié)構(gòu)剛體位移，得到單元的純變形，將單元內(nèi)力的求解限制在小變形小轉(zhuǎn)動的范圍內(nèi)，從而可以方便地求得質(zhì)點(diǎn)內(nèi)力。以空間梁單元AB為例，建立以A為原點(diǎn)的單元局部變形坐標(biāo)系，x軸的方向?yàn)閱卧拈L度方向，y軸和z軸在垂直于x軸的平面內(nèi)，如圖1 所示?？紤]該單元在時間段ta～tb內(nèi)從AB運(yùn)動到A′B′，質(zhì)點(diǎn)A′和B′經(jīng)過虛擬逆向運(yùn)動，分別運(yùn)動到A′′和B′′，從而可以求得梁單元AB的純變形分量ud，具體的推導(dǎo)過程可參見文獻(xiàn)[19]。單元純變形分量ud可表達(dá)為：

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

由單元的純變形可求得單元在局部變形坐標(biāo)系下的相應(yīng)的單元質(zhì)點(diǎn)內(nèi)力，經(jīng)過坐標(biāo)變換即可得到全局坐標(biāo)系下的質(zhì)點(diǎn)內(nèi)力。有限質(zhì)點(diǎn)法具有統(tǒng)一的計(jì)算求解框架，不同單元類型的質(zhì)點(diǎn)內(nèi)力求解方法的實(shí)現(xiàn)只需在相應(yīng)的求解模塊中進(jìn)行修改，易于在程序中實(shí)現(xiàn)。下面將介紹TM 纖維梁單元的質(zhì)點(diǎn)內(nèi)力求解方法。

2 基于TM 纖維梁單元模型的質(zhì)點(diǎn)內(nèi)力求解

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

在纖維梁單元模型中，梁截面被離散為一定數(shù)量的纖維，幾種典型截面的纖維離散方式如圖2所示。纖維梁單元模型假定截面在單元變形過程中始終保持平面，假設(shè)每個纖維截面上的應(yīng)變均勻分布，可取每根纖維中心點(diǎn)的應(yīng)變作為其代表值。在得到各個纖維的應(yīng)變后，結(jié)合纖維材料的本構(gòu)關(guān)系可以計(jì)算得到截面各個位置的應(yīng)力狀態(tài)，對截面內(nèi)力進(jìn)行積分計(jì)算得到單元內(nèi)力。纖維梁單元模型通過獨(dú)立求解截面上各纖維的應(yīng)力狀態(tài)實(shí)現(xiàn)了對梁單元截面的塑性發(fā)展過程的模擬。

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

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

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

2.2 單元位移描述

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

2.3 纖維狀態(tài)的確定

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

材料達(dá)到屈服后應(yīng)變增量可分為彈性部分和塑性部分，即：

其他材料的彈塑性本構(gòu)關(guān)系可根據(jù)具體情況作相應(yīng)的推導(dǎo)與修改。

2.4 單元質(zhì)點(diǎn)內(nèi)力的求解

單元純變形分量和與之對應(yīng)的單元質(zhì)點(diǎn)的基本內(nèi)力之間的關(guān)系由虛功原理確定，它們滿足以下方程：

3 算例分析

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

3.1 空間直角剛架

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

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

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

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

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

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

3.2 空間框架

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

圖7 空間框架Fig. 7 Space frame

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

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

3.3 六角空間剛架

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

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

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

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

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

3.4 單層球面網(wǎng)殼結(jié)構(gòu)

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

表1 計(jì)算機(jī)硬件參數(shù)Table 1 Computer hardware parameters

圖11 單層球面網(wǎng)殼結(jié)構(gòu)Fig. 11 Single-layer spherical reticulated shell structure

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

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

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

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

4 結(jié)論

本文將纖維梁模型基本原理與有限質(zhì)點(diǎn)法相結(jié)合，并考慮纖維材料的三維本構(gòu)關(guān)系，提出了考慮剪切和扭轉(zhuǎn)變形的有限質(zhì)點(diǎn)法纖維梁單元模型，基于Timoshenko 梁理論推導(dǎo)了相應(yīng)的單元質(zhì)點(diǎn)內(nèi)力求解公式，為處理考慮大位移、大轉(zhuǎn)動的彈塑性問題提供了一種有效的解決方案。將纖維梁模型與有限質(zhì)點(diǎn)法相結(jié)合，可以充分發(fā)揮兩者的優(yōu)勢，在處理幾何和材料雙非線性問題中具有求解簡單、易于實(shí)現(xiàn)的特點(diǎn)。數(shù)值算例分析表明，本文提出的纖維梁單元能較好地模擬結(jié)構(gòu)彈塑性變形及結(jié)構(gòu)大變形失穩(wěn)等問題，且能夠考慮桿件截面的塑性發(fā)展過程，相較于FPM 塑性鉸法具有更高的精度，相較于傳統(tǒng)有限元法具有更高的計(jì)算效率。