Timoshenko組合梁動力特性與瞬態(tài)響應(yīng)的求積元分析

李曉偉， 何光輝

(1. 上海開放大學(xué) 松江分校, 上海 201600; 2.浙江海洋大學(xué) 港航與交通運輸工程學(xué)院, 浙江 舟山 316022)

組合梁結(jié)構(gòu)具有較高的承載力-自重比，如鋼-混凝土組合橋梁，木-混凝土樓板，可以較合理地發(fā)揮其各組成材料的力學(xué)性能，在土木工程中具有廣泛的應(yīng)用。

雙層組合梁的力學(xué)模型與計算問題備受關(guān)注。其中，較為經(jīng)典的是Newmark[1]模型。該模型將細(xì)長組合梁的各層子梁簡化為Euler-Bernoulli梁，考慮子梁界面間的彈性相互作用(下文中簡稱“部分作用”)，至今仍具有廣泛的應(yīng)用。此外，基于Newmark模型的改進(jìn)仍然活躍，如Timoshenko梁模型開始被用于描述各子梁，對應(yīng)的計算方法也同樣得到了發(fā)展。Xu等[2]不僅得到了部分作用Timoshenko組合梁靜力彎曲和動力特性分析的變分原理及近似解析解，而且Xu等[3]利用梁的Euler-Bernoulli和Timoshenko模型解析解之間的關(guān)系，建立Newmark模型與Timoshenko組合梁模型位移解之間的關(guān)系，得到了簡便的Timoshenko組合梁解析解法。Ecsedi等[4]和Schnabl等[5]均得到了Timoshenko組合梁靜力彎曲的閉合解析解。Girhammar等[6]利用Hamilton原理得到了彈性邊界Newmark模型動力模態(tài)分析的精確解析解法。然而，當(dāng)組合梁結(jié)構(gòu)或邊界復(fù)雜時，基于微分方程或變分原理的解析解法顯得心有余而力不足，為此，Ranzi等[7]利用Newmark模型靜力問題的通解構(gòu)造廣義形函數(shù)，得到了Newmark模型靜力問題的直接剛度分析法。由于其形函數(shù)精確，它的分析結(jié)果與閉合解析解等價。此外，Nguyen等[8]建立了部分作用Timoshenko組合梁的直接剛度法，且進(jìn)行了組合梁的線性靜力分析，而后考慮組合梁界面的離散剪力鍵，得到了精確的Timoshenko組合梁直接剛度離散方程。

對于組合梁的非線性力學(xué)問題，解析解法與直接剛度法都存在較大的應(yīng)用難度，而有限元法在組合梁結(jié)構(gòu)分析方法中仍處于主導(dǎo)地位。有限元法在組合梁的極限承載能力分析[9]、大撓度分析[10]等方面均具有大量的應(yīng)用。雖然有限元法元法取得了空前成功，但由于有限單元本身的節(jié)點和自由度為常數(shù)，目前有限元商業(yè)軟件提供的單元庫一般僅包含物理量從線性到三次的插值，因為每增加一類插值形式就意味著重新推導(dǎo)和開發(fā)新的有限單元。對于變化迅速物理場的插值，一般需要提高網(wǎng)格密度或提高單元的插值階次，然而，對于最高僅三階多項式插值的有限單元庫，往往傾向于網(wǎng)格加密。

微分求積元法[11-12]不僅繼承有限元法的優(yōu)點，而且由于采用了Shu等[13]建立的快速遞推公式逼近物理場，而非有限元法中以顯式形函數(shù)近似插值基本未知量的方法，具備可變的單元內(nèi)節(jié)點分布與數(shù)目[14]，可等價為任意高階插值有限單元。目前微分求積(元)法在結(jié)構(gòu)力學(xué)領(lǐng)域內(nèi)已得到廣泛關(guān)注，如肖乃佳[15]和王宇[16]分別應(yīng)用微分求積元法于框架和板的結(jié)構(gòu)力學(xué)分析；楊驍?shù)萚17]利用微分求積法進(jìn)行了Euler-Bernoulli組合梁的動力模態(tài)分析；申志強等[18]應(yīng)用微分求積元法進(jìn)行了Newmark模型的大撓度靜力分析。

本文以雙層Timoshenko組合梁自由振動和瞬態(tài)響應(yīng)分析為例，建立其求解的微分求積單元法與常規(guī)有限單元法，并開發(fā)相應(yīng)的數(shù)值程序，通過對多階固有頻率及地震時程響應(yīng)的分析所需CPU時間測試，檢驗求積元法對常規(guī)有限元法的計算效率提高程度。

1 Timoshenko組合梁的運動學(xué)方程

如圖1(a)所示，具有對稱橫斷面的雙層組合梁由上層的c子梁和下層的s子梁通過剪力鍵連結(jié)協(xié)同工作，且在各層子梁的橫截面形心分別設(shè)置笛卡爾坐標(biāo)系xyizi,i=c,s；上、下層子梁的形心分別距離界面h1和h2。當(dāng)組合梁承受關(guān)于xzi面對稱的荷載之后，子梁i的橫斷面逆時針轉(zhuǎn)動θi，形心產(chǎn)生x軸正向的位移ui0。圖1(b)為組合梁軸向和橫向位移示意。假定組合梁的任意一層子梁i的軸向位移滿足平截面假定，即變形前物質(zhì)點(x,zi)具有x軸正向位移

Ui(x,zi)=ui0(x)-ziθi(x),i=c,s

(1)

組合梁的任意一層子梁i橫斷面內(nèi)的物質(zhì)點(x,zi)具有相同的zi向位移

Wi(x,zi)=W(x),i=c,s

(2)

(a) 組合梁立面和橫斷面示意圖

(b) 軸向和橫向位移示意圖圖1 Timoshenko組合梁的結(jié)構(gòu)示意圖和變形假定Fig.1 Diagram of Timoshenko composite beam structure and displacements assumption

2 弱式微分求積元與有限元

2.1 虛功原理

組合梁受到x-zc面內(nèi)zc坐標(biāo)軸負(fù)向，大小為q(x,t)分布荷載,如圖2所示。假定：①剪力鍵的分布密度足夠大，其剪切內(nèi)力近似為沿組合梁界面連續(xù)分布的剪力流Fsh及法向約束力流Fw；②剪力鍵在zc方向剛度足夠大，忽略子梁間的橫向相對運動；③剪力鍵具有大小為kcs的剪切剛度，即Fsh=kcsucs，其中，橫坐標(biāo)x處，兩層子梁的界面滑移為

ucs(x)=Uc(x,-h1)-Us(x,h2)

(3)

圖2 面內(nèi)荷載作用下雙層組合梁的剪力鍵內(nèi)力示意圖Fig.2 Two-layer composite beams subjected to in-plane load and its shear connectors’ internal forces

組合梁動力平衡時，存在虛功原理

(4)

式中：σi和εi分別為i層子梁的應(yīng)力和應(yīng)變向量

σi=[σiτi]T,εi=[εiγi]T

(5)

ui和FIi分別為i層子梁的位移和慣性力向量

(6)

2.2 弱式微分求積元離散方程

Shu等基于Bellman等[19-20]提出的微分求積離散原理，建立了任意連續(xù)函數(shù)F(x,t)在離散節(jié)點xi上的數(shù)值微分公式

(7)

(8)

其中，

(9)

δuk0=[uk0(x1,t)uk0(x2,t) …uk0(xN,t)]T

(10)

圖3 自然坐標(biāo)與物理坐標(biāo)之間的線性映射Fig.3 A mapping relating the natural and physical coordinates for grid points of the composite beam

考慮圖3所示的線性映射

(11)

式(8)可等價為

(12)

類似地,θc，θs和W在xi處關(guān)于x的m階導(dǎo)數(shù)可表示為

(13)

其中，k=c,s，及

(14)

δθk=[θk(x1,t)θk(x2,t)θk(x3,t)…θk(xN,t)]T

(15)

δW=[W(x1,t)W(x2,t)W(x3,t)…W(xN,t)]T

(16)

利用式(12)和式(13)，各子梁位移可近似為

(17)

組合梁在小位移、小變形階段，具有幾何方程

(18)

式(17)代入式(18)，利用式(12)和式(13)，得

(19)

組合梁的界面滑移

(20)

為了使所有的基本未知量均由統(tǒng)一的單元自由度向量δe表示，可做如下線性變換

(21)

式中：變換矩陣Ck,k= 1, 2, …, 5,為5N×5N的矩陣，且Ck的第i行j列元素Ck(i,j)=1，若j=5(i-1)+k，否者Ck(i,j)=0。單元自由度向量δe的成分如下

δe=[Φ1Φ2Φ3…ΦN]T

(22)

Φi=[uc0(xi)θc(xi)us0(xi)θs(xi)W(xi)],
i=1,…,N

(23)

利用線性變換式(21)，式(17)，式(19)和式(20)可轉(zhuǎn)化為

(24)

其中，

(25)

假定組合梁的各層子梁均為各向同性線彈性體，則k(k=c,s)層子梁有

σk=Ekεk

(26)

式中：彈性矩陣Ek由楊氏模量Ek、剪切模量Gk和切應(yīng)力修正系數(shù)κk(對于矩形截面一般取κk=5/6)構(gòu)成

(27)

將式(24)和式(26)代入式(4)，得弱式微分求積元離散方程

(28)

式中：質(zhì)量矩陣Me，剛度矩陣Ke，荷載向量Re定義如下

(29)

(30)

(31)

式中:Bm為m層梁的寬度；xi和zmk分別為x軸方向和zm軸方向的積分點坐標(biāo)；wi和wmk為積分點xi和zmk分別對應(yīng)的積分權(quán)系數(shù)；Nx和Nmz分別為x軸方向和zm軸方向的積分點數(shù)目。便于表達(dá)，式中引入了中間變量

(32)

當(dāng)組合梁受到地表豎直地震加速度激勵時，將組合梁的位移視為相對地表的位移，則單元荷載向量可表示為

(33)

式中：ag(t)為t時刻地表的豎向地震加速度；Ac和As分別為c層與s層子梁的橫截面積。

為了盡可能地提高M(jìn)e，Ke和Re的數(shù)值積分精度，且捕捉到單元的x坐標(biāo)端點，如圖4所示。這里x坐標(biāo)積分采用高斯-羅伯特積分公式，且微分求積離散過程中將高斯-羅伯特積分節(jié)點與求積法節(jié)點設(shè)置重合；y和z坐標(biāo)方向采用高斯-勒讓德數(shù)值積分公式。

圖4 積分點分布示意圖Fig.4 Schematic diagram that illustrates the integral points

當(dāng)結(jié)構(gòu)以頻率ω?zé)o阻尼自由振動時，有

(34)

從而得單元自由振動廣義特征值問題，

(35)

與有限元方法一致，弱式微分求積方程式(28)和式(33)可被組集為整體方程

(36)

(37)

若考慮結(jié)構(gòu)阻尼效應(yīng)，式(36)可修正為

(38)

式中：Cg為結(jié)構(gòu)阻尼矩陣，此處認(rèn)為其符合瑞利阻尼假定，即

Cg=αMg+βKg

(39)

在已知結(jié)構(gòu)第r和s階固有頻率ωr和ωs及對應(yīng)阻尼比ξr和ξs的條件下，瑞利比例系數(shù)α和β可由下式確定

(40)

2.3 有限元離散方程

為了便于弱式微分求積元法與有限元法的對比研究，這里簡要地給出Timoshenko組合梁的有限元離散方程。

圖5所示長度為Le的Timoshenko組合梁有限元在x坐標(biāo)軸方向均勻設(shè)置3個節(jié)點，具有單元自由度向量

(41)

節(jié)點自由度向量

(42)

圖5 Timoshenko組合梁的節(jié)點位移Fig.5 Nodal displacements of the Timoshenko composite beams

根據(jù)圖5中節(jié)點坐標(biāo)，x1=0，x2=Le/2和x3=Le，不難構(gòu)造基本未知量uc0,θc,us0,θs和W的拉格朗日插值

(43)

(44)

式中:I5為5階級單位矩陣；3個節(jié)點的拉格朗日插值基分別為

(45)

利用式(3)和幾何方程式(18)，可分別得界面相對滑動位移與組合梁應(yīng)變的插值

(46)

(47)

式(43)和式(46)代入虛功方程式(4)，得有限單元方程

(48)

式中：

(49)

(50)

(51)

當(dāng)組合梁受到地表豎向地震加速度激勵時，可得與式(33)相應(yīng)的地震荷載向量

(52)

同樣地，式(48)經(jīng)過單元組集，且可以考慮瑞利阻尼，修正為

(53)

對比式(28)和式(48)，不難發(fā)現(xiàn)弱式微分求積元方程和有限元法方程的形式一致。這里以3節(jié)點有限單元為例進(jìn)行了推導(dǎo)，而上節(jié)中的微分求積元的單元節(jié)點數(shù)N可以是大于1的任意整數(shù)，當(dāng)N=3時，不計積分誤差，兩種方程等價。當(dāng)N>3時，有限元方法需要從式(44)開始重新建立，而弱式微分求積元法自動完成這類工作。

3 數(shù)值校驗與收斂效率

3.1 程序校驗

對于圖6所示四組邊界條件的組合梁，均具有圖7所示的橫斷面與材料特性。對于該模型，Huang[22]基于Newmark模型假定解析地得到了其固有頻率，而文獻(xiàn)[23]同時得到了基于Timoshenko梁假定和平面應(yīng)力假定的雙層組合梁固有頻率有限元解；還給出了圖6所示組合梁受圖8所示豎直地震加速度激勵下的瞬態(tài)響應(yīng)。通過分析結(jié)果的對比，本節(jié)將校驗本文建立的求積元法及有限元法數(shù)值程序的正確性。

(a)固支-自由CS>(b)簡支-簡支SS

(c)固支-簡支CS>(d)固支-固支CC圖6 實際工程中的四組邊界條件Fig.6 Four combinations of boundary conditions in practical engineering

圖7 組合梁橫斷面與材料特性Fig.7 Dimensions and material properties of the composite beams

圖8 地表豎向地震加速度Fig.8 Vertical seismic acceleration of the ground

3.1.1 固有基頻

本節(jié)采用100個等長的常規(guī)有限單元和50個等長的5節(jié)點求積單元，即兩種方法采用的總自由度數(shù)均為1 005，分別離散四種邊界條件的組合梁，求解對應(yīng)的基頻。為了考慮長細(xì)比L/H(此處H表示組合梁的總梁高，為0.2 m，跨徑為L)對計算結(jié)果的影響，表1～表3分別給出了長細(xì)比為20、5和2情況下組合梁基頻分析值。從表1～表3可知，不管組合梁細(xì)長與否，本文采用的兩種算法計算結(jié)果與何光輝的分析基本一致，表明了本文算法與程序的正確性。Newmark模型較平面應(yīng)力模型存在較大的差異，且該差異隨著組合梁長細(xì)比的減小而放大，當(dāng)L/H=2時，Newmark模型解已扭曲至不可信，而Timoshenko模型與平面應(yīng)力模型分析值仍較為接近。

表1 基于不同方法的細(xì)長組合梁基頻,L/H=20

表2 基于不同方法的深組合梁基頻,L/H=5

表3 基于不同方法的極深組合梁基頻, L/H=2

3.1.2 地震響應(yīng)

假定圖6所示的簡支組合梁具有圖9所示的橫斷面尺寸及各層子梁和剪力鍵的物理參數(shù)，此外，考慮組合梁各階陣型具有相同的阻尼比0.05，該簡支梁受到圖8所示豎向地表加速度的激勵?；赥imoshenko梁假定，文獻(xiàn)[23]采用等時間步長Δt=0.001 s的Newmark-β直接積分法分析了這一簡支組合梁的瞬態(tài)撓響應(yīng)問題。為了對比校驗本文求積元和有限元程序，本節(jié)繼續(xù)采用“3.1.1”節(jié)中所采用的單元和節(jié)點數(shù)量，分別采用常規(guī)有限元法和求積元法均勻地離散組合梁。對于式(28)與式(53)中的時間導(dǎo)數(shù)，本節(jié)仍以時間步長為0.001 s的Newmark-β法對其離散。

圖9 地震分析中的組合梁橫斷面與材料特性Fig.9 Dimensions and material properties of the composite beams for seismic analysis

圖10(a)和圖10(b)分別給出了細(xì)長和短粗簡支組合梁地震發(fā)生6 s和12 s內(nèi)跨中相對地面的撓度記錄。圖10(a)曲線表明，本文基于Timoshenko假定的常規(guī)有限元法和求積元法所得分析值具有高度一致性，且與平面應(yīng)力模型分析值差異較小。此外，圖10(b)顯示由于短粗的組合梁較細(xì)長組合梁具有更大的橫向剛度，從而表現(xiàn)出更劇烈的振動響應(yīng)，這對時間導(dǎo)數(shù)的逼近存在一定難度，但本文算法與何光輝的分析仍然具有理想的一致性。

3.2 計算效率

提高物理場的逼近速率是微分求積法建立的初衷之一，而本文所進(jìn)行的求積元和有限元分析對比目的之一就是檢驗該速率提升程度。為此，本節(jié)將進(jìn)行圖6(c)所示固支-簡支組合梁的前10階級固有頻率分析及其地震時程響應(yīng)的計算，對比常規(guī)有限元法與求積元法的計算速度。本算例繼續(xù)采用圖9所示的橫斷面幾何尺寸與物理參數(shù)，且假定各陣型阻尼比均為0.05及梁長L= 4 m。所有計算均在具有以下主要配置的PC中完成，CPU：Intel i7 2600，內(nèi)存：DDR3 16 GB，外存：固態(tài)硬盤128 GB，所有數(shù)值程序均由Mathematica軟件編寫。

(a)L/H=20(b)L/H=3圖10 地震激勵下，不同長細(xì)比簡支組合梁的跨中撓度響應(yīng)Fig.10 Mid-span deflection responses of the SS composite beams with different slenderness ratio to the seismic excitation

3.2.1 固有頻率

本算例采用有限單元及具有5節(jié)點、6節(jié)點、7節(jié)點和8節(jié)點的求積單元分別均勻離散組合梁，同時分別進(jìn)行對應(yīng)的固有頻率分析值關(guān)于離散自由度數(shù)的收斂曲線分析，從而可得分析精度下限。

鑒于該精度估計方法，表4列出了有限元法與求積元法分析組合梁前10階具有0.01 rad/s精度的固有頻率所需的結(jié)構(gòu)最低總自由度數(shù)。表中可見，隨著固有頻率階次的提高，分析所需的總自由度數(shù)逐漸增加；在同一階頻率分析中，求積元法需要的自由度數(shù)顯著少于有限元法，此外，隨著求積單元節(jié)點數(shù)的增加，即單元插值階次的提高，所需的結(jié)構(gòu)總自由度數(shù)呈下降趨勢。作為表4的補充，圖11給出了對應(yīng)總自由度的離散后，固有頻率分析所需的CPU運行時間對比。顯而易見，總自由度數(shù)的增加必將引起顯著的運行時間需求。

表4 固有頻率分析所需自由度

為了更好地評價求積元法相對有限元法的計算效率提升程度，表5列出了有限元法與求積元法所需CPU時間比。表5顯示隨著頻率階次的提高或求積單元階次的提高，求積元法較有限元法計算效率的提升越顯著，當(dāng)計算第10階固有頻率時，8節(jié)點求積單元的分析效率達(dá)到了有限元的480倍。

圖11 基于有限元法與求積元法的固有頻率分析CPU時間Fig.11 CPU time required for the natural frequency analysis using the FEM and QEM

頻率階次NGPx=5NGPx=6NGPx=7NGPx=8118.615.816.619.8228.339.157.963.4340.661.484.184.1467.691.1113.2116.4567.691.1113.2116.4673.9119.8161.8166.37108.4157.3237.9237.98125.4215.6336.9394.49125.5294.6356.3480.310125.5294.6356.3480.3

3.2.2 時程響應(yīng)

結(jié)構(gòu)地震時程響應(yīng)的直接積分法有限元分析一般較為耗時，本節(jié)嘗試對比考察有限元法與求積元法在這一方面的效率表現(xiàn)，試圖提出一種更好的時程分析解決方案。

本計算過程采用四種等長單元離散方案，分別為：離散方案①228個有限單元，離散方案②25個5節(jié)點求積單元，離散方案③11個7節(jié)點求積單元，離散方案④20個9節(jié)點求積單元分別離散組合梁，這樣可使離散方案①～離散方案③分別具有總自由度數(shù)2 285，505和335，從而均能保證前10階固有頻率至少精確至0.01 rad/s，而離散方案④中9節(jié)點單元具有較離散方案③高的逼近階次且總自由度遠(yuǎn)超離散方案③，因此，近似視其為基準(zhǔn)離散方案。本算例均采用等步長為0.005 s的Newmark-β法進(jìn)行6 s內(nèi)組合梁跨中相對地面的撓度時程記錄。

圖12(a)～圖12(c)給出了離散方案④所得6 s內(nèi)跨中撓度分析值與離散方案①、離散方案②、離散方案③分析結(jié)果之差，圖中分別以圖例“QEM9-FEM”、“QEM9-QEM5”和“QEM9-FEM7”表示。曲線表明雖然離散方案③采用的總自由度最少，但是其分析精度不僅超過總自由度略多于離散方案②，而且分析精度遠(yuǎn)超過總自由度遠(yuǎn)多于離散方案①。

圖12 離散方案(1)～離散方案(3)相對離散方案(4)關(guān)于跨中撓度的誤差Fig.12 Mid-span deflection difference between results yielded by scheme (4) and those yielded by scheme (1)-(3)

為了更好地評價離散方案①～離散方案③的計算性能，表6總結(jié)了各項指標(biāo)。表中可見，在時程分析中，求積元法顯著節(jié)約了時程分析所需的CPU分析時間，而且隨著求積元階次的提高時間經(jīng)濟性更加顯著，如離散方案②的分析效率達(dá)到了有限元法離散方案①的18.7倍，而離散方案③將有限元法提速至43.2倍。此外，在平均誤差上，離散求積元法也顯著小于有限元法。

表6 四種離散方案的計算情況

4 結(jié) 論

基于Shu等的微分近似離散原理及組合梁動力問題的虛功原理，本文推導(dǎo)了部分作用雙層Timoshenko組合梁自由振動與瞬態(tài)響應(yīng)分析的微分求積單元法與常規(guī)有限單元法離散方程，開發(fā)了對應(yīng)的數(shù)值程序。通過數(shù)值校驗和計算效率分析得到以下主要結(jié)論：

(1)不論固有頻率特征值分析，還是直接積分法時程分析，高階求積元法的求解效率遠(yuǎn)超過拋物線常規(guī)有限單元。

(2)隨著求解問題所需總自由度數(shù)的提高，求積元的求解效率優(yōu)勢更加顯著。

(3)隨著求積元插值階數(shù)的提高，如從5節(jié)點求積元提高至8節(jié)點求積單元，其固有頻率與瞬態(tài)響應(yīng)的求解效率依次提高。