陳 沖，袁行飛，段元鋒，錢若軍

(1.浙江大學 空間結構研究中心，杭州 310058；2. 同濟大學 土木工程學院，上海 200092)

?

基于精細梁模型的向量式有限元分析

陳 沖1，袁行飛1，段元鋒1，錢若軍2

(1.浙江大學 空間結構研究中心，杭州 310058；2. 同濟大學 土木工程學院，上海 200092)

精細梁不同于Euler梁和Timoshenko梁，該模型在考慮剪切變形的同時還考慮了橫向彎曲時截面轉動產生的附加軸向位移及橫向剪切變形影響截面抗彎剛度后產生的附加橫向位移。推導了適用于向量式有限元分析的精細梁單元應變和內力表達式，采用FORTRAN自編了向量式有限元程序。對懸臂梁、兩端固支梁和門式框架進行了算例分析，對比了采用不同梁單元模型下結構的豎向位移。結果表明：當高跨比較小時，3種梁單元的豎向位移相差不大；當高跨比較大時，精細梁單元的豎向位移較Euler梁和Timoshenko梁明顯增大，表明剪切變形及剛度折減引起的附加軸向位移、附加橫向位移不能忽略。精細梁單元模型對高跨比較大的梁進行分析可望得到更精確的結果。

Euler梁；Timoshenko梁；精細梁；向量式有限元；高跨比

向量式有限元(Vector Form Intrinsic Finite Element，VFIFE)是由美國普渡大學丁承先教授提出的一種新型有限元方法，已在土木工程領域得到了較好應用[1-4]。目前，向量式有限元中廣泛采用的梁模型是經典的Euler-Bernoulli梁模型，如喻瑩等[5]對Euler梁結構的幾何破壞、材料破壞以及倒塌破壞作了分析，張若京等[6]基于Euler梁分析了懸臂梁的幾何非線性行為，梁育銘等[7]模擬了Euler梁、桿構件組成結構的施工力學行為?？梢钥隙ǖ氖?，當所研究梁的高度遠小于跨度時，采用Euler梁模型進行分析可得到較精確的解。當梁高跨比較大時，由于剪切影響顯著，Euler梁模型不再適用。為考慮剪切變形影響，Timoshenko[8]于1921年提出了剪切變形理論，后人稱其為Timoshenko梁模型。Thomas[9]、Nickell等[10]使用不同階數的多項式描述變形和旋轉的變化，分析深梁得到了較合理結果，但分析細長梁時產生了剪力自鎖現象。Reddy[11]提出IIE(Interdependent Interpolation Element, IIE)法解決剪力自鎖問題。李彥輝等[12]將IIE法應用于向量式有限元，探討了剪切變形對梁結構的影響，但沒有考慮變形間的耦合作用。本課題組錢若軍等[13]基于變形分析建立了考慮變形耦合作用的空間梁任意點位移表達式，并推導了兩節(jié)點十自由度的空間梁單元模型。上述研究主要基于傳統(tǒng)有限元方法展開，在向量式有限元中的有效性和可行性還有待于驗證。

本文以考慮變形耦合作用的空間精細梁理論為基礎，推導了適用于向量式有限元的精細梁單元應變和內力計算公式，采用FORTRAN自編了向量式有限元程序。并通過算例分析驗證本文精細梁模型的正確性。

1 向量式有限元及梁單元

向量式有限元是一種以向量力學為理論基礎，以數值計算為描述方法的分析方法[14]。

1.1 向量式有限元基本概念

1)點值描述 將構件離散成有限的空間點，并通過一組內插函數來描述構件上其他點的空間位置。經過簡化以后，構件的質量可以按照一定的方法分配到各點并且可以用牛頓運動定律分析結構。

2)途徑單元 把構件上空間點的時間軌跡離散成一組時間點，通過一組標準的控制方程，空間點從一個時間點的位置到達另一個時間點的位置，這一時間段稱為途徑單元。簡化后，可以用大變位和小變形理論分析大變形問題。

3)逆向運動 如圖1所示，桿件從12位置平移至1′2′位置，再以點1′為中心逆向旋轉至1d2d，1d2d與1a2a在同一直線的位置，這一過程即逆向運動。向量式分析通過虛擬的逆向運動來處理純變形，進一步由純變形得到單元內力。

圖1 虛擬逆向運動Fig.1 Virtual inverse motion

圖2 平面梁系結構變形Fig.2 Deformation of plane beam

1.2 Euler梁單元

如圖2所示，假設單元由初始位置1a2a運動到下一時間點12，根據Euler梁的撓曲理論，桿件任意截面上一點的變形滿足下列關系[14]：

(1)

式中：um為x軸與截面交點的軸向位移，v為y方向位移，φ為空間轉角。采用傳統(tǒng)有限元的形函數，可得um、v、φ表達式：

(2)

根據虛功原理，可得單元內力與節(jié)點變形關系為

(3)

進一步由節(jié)點力平衡，求得節(jié)點其它內力分量。

(4)

圖3 梁單元節(jié)點力和彎矩Fig.3 Particle force of plane beam

2 精細梁模型

2.1 拉壓、彎、剪、扭、翹曲及其耦合效應作用下空間梁任意一點的位移

在外荷載作用下，空間梁同時發(fā)生軸向拉壓、彎曲、剪切和扭轉。相應會產生軸向、彎曲和剪切變形，以及橫截面的翹曲變形，如圖4[15]。

圖4 梁的受力狀態(tài)Fig.4 the stress state of beam

圖5等截面梁慣性系Fig.5 the coordinate of space beam

在坐標系O-XYZ(如圖5所示)中，等截面梁任意一點的位移為[15]

u=uT+ub+ubs+us+uΔ+uω

v=vb+vbs+vs+vΔ

w=wb+wbs+ws+wΔ

(5)

其中：ξ=ρθx；uT為軸向拉壓產生的軸向位移；ub為在彎矩作用下梁發(fā)生橫向彎曲時截面轉動產生軸向位移；ubs為考慮梁橫向剪切變形影響截面的抗彎剛度后產生的附加撓度而使得梁截面轉動產生的修正軸向位移；us、vs、ws為剪力產生的軸向位移和橫向剪切位移；vb、wb為由于梁在彎矩作用下產生的橫向彎曲位移；vbs、wbs為考慮梁橫向剪切變形影響截面的抗彎剛度后，產生的附加橫向位移；uΔ、vΔ、wΔ為考慮軸力二次效應影響產生的軸向位移及橫向位移；uω、ξ為約束扭轉引起的軸向位移及切向位移；θx為截面扭轉角[15]。

本文暫不考慮二次效應影響以及約束扭轉引起的軸向位移及切向位移，對以上梁單元作適當簡化，等截面空間梁中任意一點的位移為

u=uT+ub+ubs+us

v=vb+vbs+vs

w=wb+wbs+ws

(6)

其中：uT、vb、wb為常規(guī)Euler梁中的軸向位移和橫向彎曲位移，ws、vs為橫向剪切位移，ub、ubs、vbs、wbs、us為考慮變形耦合的修正位移。

2.2 修正位移表達式

根據精細梁的理論可得各項修正位移為:

1)橫向彎曲導致截面轉動產生的軸向位移ub[15]。彎曲導致截面轉動從而產生的沿x軸方向的位移為

uvb=-yθzb，uwb=-zθyb

(7)

其中：θzb與θyb為空間梁彎曲變形后中和軸的轉角。

則由于梁橫向彎曲時截面轉動而產生的軸向位移為

ub=uvb+uwb

(8)

2)橫向剪切變形影響截面抗彎剛度后產生的軸向位移ubs和橫向位移vbs、wbs[15]。考慮梁的抗彎剛度受橫向剪切變形影響后，產生的附加撓度而使得梁截面轉動產生的軸向位移ubs為

ubs=uvs+uws=βysuvb+βzsuwb

(9)

對于y、z軸方向截面抗彎剛度受剪切變形影響后產生的附加撓度而使得梁截面轉動產生的橫向位移分別有

vbs=βysvb,wbs=βzswb

(10)

3)橫向剪切變形引起的軸向位移us[15]。梁的剪切變形引起的沿x軸的位移us,y和us,z分別沿截面高度和寬度按三次曲線變化，即

(11)

剪切變形引起的沿x軸的總位移us:

us=us,y+us,z

(12)

式中：h、b分別為梁橫截面高度和寬度。

3 應用于向量式有限元的精細梁模型

3.1 精細梁單元節(jié)點位移及應變

基于精細梁理論，可得二維梁單元的位移公式為

u=uT+ub+ubs+us

v=vb+vbs+vs

(13)

則有：

(14)

平面梁單元節(jié)點處的變形滿足下列條件

(15)

代入式(17)可得6個常數：

(16)

將所得的常數代入原方程式，可得

(17)

(18)

(19)

3.2 單元軸向應變和剪應變對應的等效內力

積分可得單元軸向應變對應的等效內力

(20)

如梁截面為矩形(寬b，高h)，則梁對應的內力增量為

(21)

將兩部分內力疊加可得總的內力增量：

(22)

進一步由靜力平衡條件可得f1x、f2y、f1y。

4 向量式有限元靜力分析算例

算例1

懸臂梁如圖7所示[12]，長度L=2.0 m，梁寬b=0.2 m，高h=0.2 m，彈性模量E=2×109N/m2，單位長度質量m=2.0 kg，剪切效應的截面影響系數Ks=5/6，泊松比υ=0.25，剪切模量G=8×108N/m2，梁的右端點處承受外力P=2.5×103N，時間步長Δt=10-5s。用VFIFE編寫Euler梁、Timoshenko梁和精細梁程序對懸臂梁進行計算，可得表1所示自由端豎向位移。

表1 不同分析方法所得的懸臂梁自由端豎向位移

由表1可見，采用考慮了剪切變形和剛度折減影響的精細梁模型后得到的自由端節(jié)點豎向位移最大；只考慮剪切影響的Timoshenko梁得到的位移次之；Euler梁得到的位移最小。但由于高跨比(h/L=1/10)較小，屬于工程梁范疇(h/L=1/8～1/12)，剪切變形和剛度折減影響均較小，因此，采用上述3種梁模型所得結果相差不大。此外，VFIFE所得的結果和ANSYS結果較接近，可驗證所編VFIFE程序的正確性。

分別采用ANSYS實體模型、Timoshenko梁模型和向量式精細梁模型對高跨比(h/L=1/10,1/5,1/3)的上述懸臂梁進行分析，可得表2所示自由端豎向位移。

表2 不同分析方法所得的不同高跨比懸臂梁自由端豎向位移

由表2可見，當高跨比h/L=1/10時，3種方法所得的位移值基本一致；隨著高跨比變大，三者之間的差距也變大。當高跨比為1/3時，實體模型所得的結果最大且與精細梁所得的結果相差不大；Timoshenko梁所得的結果最小且與其他兩種方法所得結果相差較大。因此，在進行深梁分析時，精細梁模型所得的結果較Timoshenko梁所得更為準確。

圖6 懸臂梁承受外力作用Fig.6 Force on cantilever beam

圖7 深梁受均布荷載作用Fig.7 Distributed load on deep beam

算例2

考慮兩端為固定端的深梁如圖8所示[16]。彈性模量E=2×1011N/m2，單位長度質量m=2.0 kg，剪切效應的截面影響系數Ks=5/6，泊松比υ=0.3，剪切模量G=7.692×1010N/m2，梁承受均布外力q=6×107N/m2作用，懸臂梁分為10個單元。分別采用Euler梁(VE)、Timoshenko梁(VT)和精細梁(VP)模型對不同高跨比(H/L=1/3,1/2,1)的固端梁進行分析，可得圖9所示梁位移。圖中x為單元節(jié)點橫坐標，v為單元節(jié)點y向位移。

如圖9所示，當高跨比為H/L=1/3,1/2,1時，采用Euler梁(VE)所得結果明顯小于其它兩種模型，表明Euler梁(VE)用于深梁或淺深梁分析存在較大誤差，所得結果偏小。當高跨比為1/3時，Timoshenko梁(VT)結果與精細梁(VP)結果基本接近，但隨著高跨比增大，兩者差距逐漸增大，表明剛度折減影響逐漸增大。當梁較深時，剛度折減影響不容忽略。因此，對于高跨比較大的深梁或淺深梁，采用同時考慮剪切變形和剛度折減影響的精細梁模型是必要的。

算例3

圖8 不同高跨比時梁位移Fig.8 Displacement of the beam with different depth-span ratio

圖9 平面門式框架Fig.9 The portal frame

圖10 高跨比為1/30 的荷載位移曲線Fig.10 Load-displacement curves when h/L=1/30

圖11 高跨比為1/6的荷載位移曲線Fig.11 Load-displacement curves when h/L=1/6

由圖可見，當高跨比為1/30時，采用3種梁模型所得結果基本接近，表明工程梁分析時，剪切變形和剛度折減的影響較小，可以忽略。當高跨比為1/6時，由上述公式可得剪切變形和剛度折減的所產生的位移修正項變得較大。精細梁模型同時考慮了兩者的影響，Timoshenko梁模型只考慮了剪切變形的影響，而Euler梁模型兩者的影響都沒有考慮，所以VE所得的屈曲荷載最大，VP所得的屈曲荷載最小，此時剪切變形和剛度折減的影響不能忽略。

5 結論

以精細梁理論為基礎，綜合考慮了剪切變形及附加軸向位移、附加橫向位移對單元內力和變形的影響，推導了基于精細梁模型的向量式有限元基本公式。通過編寫FORTRAN程序分析算例得到以下結論：

1)對于高跨比較小的工程梁，Euler梁、Timoshenko梁以及精細梁所得的結果十分接近，表明此時剪切變形和剛度折減的影響較小，可以忽略。

2)對于高跨比較大的深梁或淺深梁，精細梁所得的結果比Timoshenko梁和Euler梁所得的結果大得多，此時梁的剪切變形和剛度折減的影響較大，不能忽略。

3)對工程梁進行分析時，由于剪切變形和剛度折減影響較小，選擇較為簡單的Euler梁模型是可行的。但對于超出工程梁范圍的深梁或淺深梁，采用考慮剪切變形和剛度折減影響的精細梁模型可望得到更準確的結果。

[1] 丁承先，王仲宇，吳東岳，等. 運動解析與向量式有限元[R]. 中國，臺灣：中央大學工學院橋梁工程研究中心，2007.

[2] Ting E C，Shih C，Wang Y K. Fundamentals of a vector form intrinsic finite element：Part I. Basic procedure and a plane frame element [J]. Journal of Mechanics，2004，20(2)：113-122.

[3] Ting E C，Shih C，Wang Y K. Fundamentals of a vector form intrinsic finite element：Part II.

[4] Shih C，Wang Y K，Ting E C. Fundamentals of a vector form intrinsic finite element：Part III. Convected material frame and examples [J]. Journal of Mechanics，2004，20(2)：133-143.

[5] 喻瑩.基于有限質點法的空間鋼結構連續(xù)倒塌破壞研究[D]. 浙江大學，2012.

[6] 班自愿,張若京.懸臂梁大變形的向量式有限元分析[J]. 計算機輔助工程,2012,19:25-29.

[7] 梁育銘，高樓施工之向量式有限元模擬[D]. 中原大學，2012.

[8] Timoshenko S P. On the correction for shear of the differential equation for transverse vibrations of prismatic bars [J].Philosophical Magazines，1921，41：744-746.

[9] Thomas J, Abbas B A H. Finite element model for dynamic analysis of Timoshenko beam [J]. Journal of Sound and Vibration，1975，41(3)：291-299.

[10] Nickell R E, Secor G A. Convergence of consistently derived Timoshenko beam finite elements [J]. International Journal for Numerical Method in Engineering，1972，5：243-245.

[11] Reddy J N. On locking-free shear deformable beam finite elements [J].Computer Methods in Applied Mechanics and Engineering，1997，149：113-132.

[12] 李賢華.向量有限元方法應用于Timoshenko梁分析之研究[D].臺灣高雄：臺灣國立中山大學，2009.

[13] 錢若軍，袁行飛，林智斌. 固體和結構分析理論及有限元法[M]. 南京：東南大學出版社，2013.

[14] 丁承先，段元鋒，吳東岳. 向量式結構力學[M]. 北京：科學出版社，2008.

[15] 馬春來.考慮剪切位移和多因素耦合影響下的空間梁非線性分析模型研究[D]. 浙江大學，2008.

[16] Ahmed S R，Idris A B M，Uddin M W. Numerical solution of both ends fixed deep beams [J]. Computer and Structures，1996，61:21-29.

(編輯 胡 玲)

Vector form intrinsic finite element analysis based on fine beam model

ChenChong1，YuanXingfei1，DuanYuanfeng1，QianRuojun2
(1. Space Structures Research Center，Zhejiang University，Hangzhou 310058，P.R.China;2. College of civil Engineering，Tongji University，Shanghai 200092，P.R.China)

The proposed fine beam model is different from that of Euler beam and Timoshenko beam. Some effects were considered in the new model such as shear displacement，the additional axial displacement produced by lateral bending and the additional transverse displacement produced by reduced stiffness due to transverse shear deformation. The formula of strain and internal force of the fine beam model，applied to Vector Form Intrinsic Finite Element (VFIFE) analysis，were derived and corresponding programs were developed by Fortran language. Cantilever beam，both ends clamped beam and portal frame are analyzed and the vertical displacements were compared using different beam element models. Numerical results showed that when the depth-span ratio was relatively small，the vertical displacements obtained by different beam model had slight difference. When the depth-span ratio was larger，the vertical displacement obtained by the fine beam model was obviously larger than that obtained by the Euler beam and Timoshenko beam. Therefore，the shear displacement，the additional axial displacement and the additional transverse displacement caused by stiffness reduction should not be ignored when the deep beam was analyzed. The new beam model proposed in this paper demonstrated more accurate results when the beam had a large depth-span ratio.

Euler beam model；Timoshenko beam model；fine beam model；VFIFE；depth-span ratio

10.11835/j.issn.1674-4764.2015.02.001

2014-10-11 基金項目：國家自然科學基金項目(51278461)；浙江省重點科技創(chuàng)新團隊(2010R50034)

陳 沖(1986-)，男，博士，主要從事大跨空間結構相關研究，(E-mail)lhcc1986@163.com。

Foundation item：National Natural Science Foundation of China(51278461); Zhejiang Province Key Science & Technology Innovation Team(No.2010R50034)

TU311.4

A

1674-4764(2015)02-0001-07

Received：2014-10-11

Author brief：Chen Chong (1986-), doctoral candidate, main research interest: space structures, (E-mail)lhcc1986@163.com。