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

關(guān)鍵詞： Euler梁；Timoshenko梁；精細梁；向量式有限元；高跨比

中圖分類號：TU311.4 文獻標(biāo)志碼：A 文章編號：1674-4764（2015）02-0001-07

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

陳 沖，等：基于精細梁模型的向量式有限元分析

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

1 向量式有限元及梁單元

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

1.1 向量式有限元基本概念

1）點值描述 將構(gòu)件離散成有限的空間點，并通過一組內(nèi)插函數(shù)來描述構(gòu)件上其他點的空間位置。經(jīng)過簡化以后，構(gòu)件的質(zhì)量可以按照一定的方法分配到各點并且可以用牛頓運動定律分析結(jié)構(gòu)。

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

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

1.2 Euler梁單元

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

2 精細梁模型

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

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

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

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

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

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

2.2 修正位移表達式

根據(jù)精細梁的理論可得各項修正位移為：

1）橫向彎曲導(dǎo)致截面轉(zhuǎn)動產(chǎn)生的軸向位移ub[15] 。彎曲導(dǎo)致截面轉(zhuǎn)動從而產(chǎn)生的沿x軸方向的位移為

uvb=-yθzb，uwb=-zθyb（7）

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

則由于梁橫向彎曲時截面轉(zhuǎn)動而產(chǎn)生的軸向位移為

ub=uvb+uwb（8）

2）橫向剪切變形影響截面抗彎剛度后產(chǎn)生的軸向位移ubs和橫向位移vbs、wbs[15]?？紤]梁的抗彎剛度受橫向剪切變形影響后，產(chǎn)生的附加撓度而使得梁截面轉(zhuǎn)動產(chǎn)生的軸向位移ubs為

ubs=uvs+uws=βysuvb+βzsuwb（9）

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

vbs=βysvb，wbs=βzswb（10）

式中：βys=-φy1+2φy，βzs=-φz1+2φz。

φy、φz是由于剪切影響對y軸和z軸方向抗彎剛度的影響系數(shù)，φy=12KsEIzGAsyL2，φz=12KsEIyGAszL2、Asy，Asz為桿截面沿y軸和z軸的有效抗剪面積；Iy，Iz是y軸和z軸主慣性矩，L是相當(dāng)長度，Ks為剪切效應(yīng)的截面影響系數(shù)。

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

us，y=-y343b2vsxus，z=-z343h2wsx （11）

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

us=us，y+us，z（12）

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

3 應(yīng)用于向量式有限元的精細梁模型

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

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

u=uT+ub+ubs+us

v=vb+vbs+vs（13）

假設(shè)位移形函數(shù)： uT=a1+a2x，vb=-a3x36+a4x22+a5x+a6，vs=Ωa3xl2[12]。其中：Ω=EIaGAaKsl2，G為剪切彈性模量。根據(jù)式（8）、（9）、（10）可得ub=φy，ubs=βysub，us=-4y33b2Ωa3l2，vbs=βysvb。

則有：

u=a1+a2x+1+βysφy-4y33b2Ωa3l2

v=-1+βysa3x36+a4x22+a5x+a6+Ωa3xl2（14）

當(dāng)取時間較小的途徑單元時轉(zhuǎn)角φ近似等于-vbx，

φ=a3x22+a4x+a5

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

x=0，y=0，u=0，v=0，φ=θ1x=l，y=0，u=Δe，v=0，φ=θ2 （15）

代入式（17）可得6個常數(shù)：

a1=0，a2=Δel，a3=61+βysl2（1+βys+12Ω）θ1+θ2

a4=-3（1+βys）l（1+βys+12Ω）θ1+θ2+θ2-θ1l，a5=θ1，a6=0（16）

將所得的常數(shù)代入原方程式，可得

uvφ=1+βys1+βys+12ΩR11R12R130R22R230R32R33Δeθ1θ2（17）

R11=x1+βys+12Ωl1+βys，

R12=y31+βysx2l2-[4（1+βys）+12Ω]xl+

1+βys+6Ω-8y2Ωl2b2

R13=y31+βysx2l2-21+βys-12Ωxl+

6Ω-8y2Ωl2b2，

R22=-1+βysx3l2+21+βys+3Ωx2l-1+βys+6Ωx

R23=-1+βysx3l2+1+βys-6Ωx2l+6Ωx，

R32=3x2l2-41+βys+3Ωxl1+βys+1+βys+12Ω1+βys

R33=3x2l2-21+βys-6Ωxl1+βys進一步可得到單元的軸向應(yīng)變和剪切應(yīng)變?yōu)?/p>

εx=dΔη2dds=duds=

1l1y1+βys1+βys+12Ω6s1+βys-41+βys+12Ωy1+βys1+βys+12Ω6s1+βys-21+βys-12Ω

Δeθ1θ2=B1de

式中s=xl（18）

γxy=uy+vx=6Ω1+βys1+βys+12Ω01-4y2b21-4y2b2Δeθ1θ2=B2e（19）

3.2 單元軸向應(yīng)變和剪應(yīng)變對應(yīng)的等效內(nèi)力

積分可得單元軸向應(yīng)變對應(yīng)的等效內(nèi)力

f2xMm1zMm2zM=Ea（1+βys）2l（1+βys+12Ω）2

Aa（1+βys+12Ω）2（1+βys）20004（1+βys）2+24Ω（1+βys）+144Ω2Ia2（1+βys）2-24Ω（1+βys）-144Ω2Ia02（1+βys）2-24Ω（1+βys）-144Ω2Ia4（1+βys）2+24Ω（1+βys）+144Ω2Ia

Δeθ1θ2（20）

如梁截面為矩形（寬b，高h），則梁對應(yīng)的內(nèi)力增量為

f2xVm1zVm2zV=36GKslΩ21+βys2（1+βys+12Ω）2

0 0 00 6Aa5-8Iab2 6Aa5-8Iab20 6Aa5-8Iab2 6Aa5-8Iab2Δ eθ1θ2（21）

將兩部分內(nèi)力疊加可得總的內(nèi)力增量：

f2xm1zm2z=f2xMm1zMm2zM+f2xVm1zVm2zV=∫VaBT1EB1dedVa+∫VaBT2GKsB2dedVa（22）

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

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

算例1

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

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

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

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

算例2

考慮兩端為固定端的深梁如圖8所示[16]。彈性模量E=2×1011N/m2，單位長度質(zhì)量m=2.0 kg，剪切效應(yīng)的截面影響系數(shù)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é)點橫坐標(biāo)，v為單元節(jié)點y向位移。

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

算例3

門式框架如圖10所示。彈性模量E=2.1×106，桿件截面寬b=0.3，高h=0.3，剪切效應(yīng)的截面影響系數(shù)Ks=5/6，泊松比υ=0.3，剪切模量G=8.077×105，橫梁中間受一豎向力作用，底端兩節(jié)點固定。分別采用Euler梁（VE）、Timoshenko梁（VT）和精細梁（VP）模型對不同高跨比（h/L=1/30，1/6）的門式框架進行分析，得到圖11、圖12所示荷載橫梁中點豎向位移曲線。

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

5 結(jié)論

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

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

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

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

參考文獻：

[1]丁承先，王仲宇，吳東岳，等. 運動解析與向量式有限元[R]. 中國，臺灣：中央大學(xué)工學(xué)院橋梁工程研究中心，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]喻瑩.基于有限質(zhì)點法的空間鋼結(jié)構(gòu)連續(xù)倒塌破壞研究[D]. 浙江大學(xué)，2012.

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

[7]梁育銘，高樓施工之向量式有限元模擬[D]. 中原大學(xué)，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]李賢華.向量有限元方法應(yīng)用于Timoshenko梁分析之研究[D].臺灣高雄：臺灣國立中山大學(xué)，2009.

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

[14]丁承先，段元鋒，吳東岳. 向量式結(jié)構(gòu)力學(xué)[M]. 北京：科學(xué)出版社，2008.

[15]馬春來.考慮剪切位移和多因素耦合影響下的空間梁非線性分析模型研究[D]. 浙江大學(xué)，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.

（編輯 胡 玲）