波譜單元法在空間桁架地震響應(yīng)分析中的應(yīng)用

張俊兵，朱宏平，王丹生，閤東東

(1.華中科技大學 土木工程與力學學院，武漢 430074;2.華中科技大學 控制結(jié)構(gòu)湖北省重點實驗室，武漢 430074)

波譜單元法在空間桁架地震響應(yīng)分析中的應(yīng)用

張俊兵1，2，朱宏平1，2，王丹生1，2，閤東東1，2

(1.華中科技大學 土木工程與力學學院，武漢 430074;2.華中科技大學 控制結(jié)構(gòu)湖北省重點實驗室，武漢 430074)

針對傳統(tǒng)波譜單元法(SEM)只能用于求解節(jié)點集中荷載作用下結(jié)構(gòu)動力響應(yīng)問題的不足，提出了一種通過計算地震等效波譜節(jié)點荷載求解桁架結(jié)構(gòu)地震響應(yīng)的方法。基于虛功原理，利用波譜形函數(shù)積分得到地震等效波譜節(jié)點荷載的顯式表達式，通過修改波譜單元法中單元剛度矩陣的波數(shù)，考慮了阻尼對結(jié)構(gòu)動力特征的影響，采用數(shù)值拉普拉斯(Laplace)變換替換快速傅里葉(FFT)變換，回避了傳統(tǒng)波譜單元法中FFT的周期性問題。利用地震荷載等效后的波譜節(jié)點荷載對三維空間桁架結(jié)構(gòu)進行地震響應(yīng)分析，結(jié)果表明，采用本文的方法能方便的計算桁架結(jié)構(gòu)的地震等效波譜節(jié)點荷載，精確求解結(jié)構(gòu)的地震響應(yīng)，與傳統(tǒng)有限元法(FEM)相比，大大減少計算單元數(shù)量，提高計算精度，且便于編程計算。

波譜單元法;阻尼;波數(shù);桁架;有限元;地震響應(yīng)

計算機技術(shù)的發(fā)展使得有限元法在各類工程問題中得到廣泛應(yīng)用。為滿足工程計算精度，采用有限元法對實際結(jié)構(gòu)進行分析時需要大量劃分單元，從而導(dǎo)致計算量大大增加［1］。相比于傳統(tǒng)的有限元法，基于連續(xù)質(zhì)量體系動力分析的波譜單元法，對連續(xù)均勻的構(gòu)件進行分析，不論構(gòu)件長度如何，在構(gòu)件的連續(xù)均勻部分只需要一個單元［2］。該方法利用連續(xù)質(zhì)量體系平衡微分方程的精確波動解作為形函數(shù)，形成與頻率相關(guān)的動態(tài)剛度矩陣并進行求解，因而能夠精確描述結(jié)構(gòu)的動力特征［3］，在不增加單元數(shù)量的情況下大大提高計算精度。Doyle［2，4，5］和 Lee［1，6］等采用基于快速傅里葉變換(FFT)的波譜單元法對集中動荷載作用下的桿件、梁和板進行分析，成功求得結(jié)構(gòu)的動力響應(yīng)。由于FFT具有周期性，利用基于FFT變換的波譜單元法進行結(jié)構(gòu)動力響應(yīng)分析時，必須在波譜分析中考慮無限長單元［4］(throw-off element)。然而，無限長單元在計算時會將系統(tǒng)能量導(dǎo)出［4］，從而導(dǎo)致計算誤差。日本學者Igawa和Komatsu等［7］利用數(shù)值Laplace變換代替FFT變換，求解有限長系統(tǒng)，取得了很好的效果。傳統(tǒng)的波譜單元法只能對集中荷載作用下的結(jié)構(gòu)進行分析，Lee等［8］利用線性疊加原理得到了分布荷載作用下梁的動力響應(yīng)，但該文的方法計算分布荷載時采用數(shù)值積分，不但計算誤差大，而且針對不同的結(jié)構(gòu)需要采用不同的積分表達式，因而不適合大型結(jié)構(gòu)的編程計算。結(jié)構(gòu)阻尼是結(jié)構(gòu)的固有特性，直接影響到結(jié)構(gòu)的動力特征，但到目前為止，大多數(shù)學者［1－9］利用波譜單元法進行動力分析時仍未考慮阻尼的影響。

本文采用基于Laplace變換的波譜單元法對桁架結(jié)構(gòu)進行分析，回避了傳統(tǒng)波譜單元法中FFT變換的周期性問題。基于虛功原理，通過對形函數(shù)進行積分計算，推導(dǎo)了地震荷載作用下桁架結(jié)構(gòu)的地震等效波譜節(jié)點荷載的顯式表達式，通過修改波譜單元剛度矩陣的波數(shù)，考慮了結(jié)構(gòu)阻尼的影響。利用等效后的地震波譜節(jié)點荷載對鋼桁架橋具有200根桿的空間桁架進行地震響應(yīng)分析，數(shù)值算例結(jié)果表明，采用本文的方法的能十分方便的進行編程計算，修正后的波譜單元法能有效求解桁架結(jié)構(gòu)的地震響應(yīng)，減少計算單元數(shù)，提高計算精度。

1 桿的平衡微分方程

考慮如圖1所示均勻線彈性桿受到軸向力作用及地震產(chǎn)生的支座激勵，根據(jù)微段受力平衡可得［10］

式中N=N(x，t)表示桿件軸向力，可以寫為:

式中:u=u(x，t)為桿件相對于單元坐標的軸向位移;E為彈性模量;A為桿件截面面積;α1為桿件內(nèi)部粘彈性阻尼系數(shù);q(x，t)為沿桿件長度方向作用的外荷載，若桿件僅受地震作用，則外荷載q(x，t)=0。

式中fI和fc分別為單位長度的慣性力和粘滯阻尼力。分別寫為:

其中:ρA為單位長度質(zhì)量;ut=ut(x，t)為桿件相對于固定參考軸的軸向位移;ug=ug(x，t)為地震作用產(chǎn)生的支座位移;c為與速度相關(guān)的外部粘滯阻尼系數(shù)。

將式(2)－式(4)代入式(1)中并化簡，得到地震作用下桿的平衡微分方程為［10］:

式中peff(t)為地震等效支座激勵荷載，寫為:

式中ag為地面加速度。

圖1 支座激勵下桿受力平衡計算簡圖Fig.1 Model of bar subjected to support excitation

2 桿的波譜單元剛度矩陣

將式(5)兩邊進行數(shù)值Laplace變換，并考慮桿件端部的邊界條件，則頻域內(nèi)桿的平衡微分方程可以簡寫為:

式中:U=U(x，w)為桿件在頻域內(nèi)的軸向位移表達式;s為Laplace變換參數(shù);Peff=Peff(x，w)為頻域內(nèi)的支座激勵荷載;x為節(jié)點坐標位置函數(shù)，如圖2所示;符號(″)表示對x求二階偏導(dǎo)。

圖2 桿的波譜節(jié)點荷載和位移Fig.2 Spectral nodal forces and displacements

為得到桿件自由振動下的波譜關(guān)系，假定式(7)右側(cè)的支座激勵荷載Peff=0。此時式(7)寫為:

式中kr為波數(shù)，寫成:

式中的波數(shù)包含了阻尼項，在既有波譜單元法的基礎(chǔ)上考慮了結(jié)構(gòu)阻尼。

式(8)的解可以表示為:

考慮相應(yīng)的位移邊界條件:

將式(11)中的位移邊界條件代入式(10)，則式(8)的解可以用單元的節(jié)點位移表示為:

式中d為節(jié)點位移向量，表示為:

式中N(x，w)為與頻率相關(guān)的形函數(shù)，表示為:

用虛功原理表示式(7)，可以寫為:

式中符號(')表示對x求偏導(dǎo)，δ為變分符號。

將式(12)代入式(15)，則桿的平衡微分方程可以表示為向量形式［10］

式中:Ke(x，w)為桿的波譜單元剛度矩陣，寫為:

式中f為波譜等效節(jié)點荷載向量，寫為:

將式(14)中的形函數(shù)N(x，w)代入式(17)，把桿的波譜單元剛度矩陣Ke(x，w)寫成矩陣形式［1，6，7］:

3 地震等效波譜節(jié)點荷載

3.1 承受均布荷載桿的等效波譜節(jié)點荷載

假設(shè)均勻線彈性桿受與沿桿軸向成α角的均布荷載 Peff(x，t)作用，通過 Laplace 變換，Peff(x，t)在頻域內(nèi)表示為Peff(x，w)，考慮其均勻分布，將 Peff(x，w)簡寫為Peff。如圖3所示，將Peff分解為沿桿軸向的分量P∥=Peffsinα和垂直于桿軸向的分量P⊥=－Peffcosα?？紤]到桿只能承受軸向力，聯(lián)立式(14)和式(18)，即可得到桿的等效波譜節(jié)點荷載為:

3.2 空間桁架桿的地震等效波譜節(jié)點荷載

考慮圖4所示的空間桁架桿受到水平地震作用。圖中x－y－z坐標表示整體坐標系，x'－y'－z'坐標表示桿的單元坐標系，桿單元的兩個節(jié)點分別記為i和j，空間位置坐標分別為(xi，yi，zi)和(xj，yj，zj)。L 表示桿長，CX，CY，CZ分別表示桿的軸向方向與 x，y，z軸夾角的余弦，可以用節(jié)點坐標表示為:

圖4中α為桿與其在yz平面內(nèi)投影的夾角，可以寫為:

對比圖3與圖4，可以看出，如果空間桁架桿僅受水平地震作用，則桿的地震等效波譜節(jié)點荷載就可以通過式(20)計算得到。

考慮實際計算時，桿的節(jié)點荷載通常在整體坐標下給出，對式(20)中的地震等效波譜節(jié)點荷載進行進一步分解(如圖5所示)，整體坐標系下地震等效波譜節(jié)點荷載各分量可以表示為:

利用式(23)得到的整體坐標系下的地震等效波譜節(jié)點荷載作為集中荷載，代入傳統(tǒng)的波譜單元法中，即可求得空間桁架結(jié)構(gòu)的地震響應(yīng)。

4 波譜單元法的計算框架

4.1 荷載向量Laplace變換及位移向量Laplace逆變換

由于FFT變換的周期性，采用基于FFT的波譜單元法主要用來求解半無限長和無限長構(gòu)件［7］。本文采用文獻［7］中的方法，利用數(shù)值Laplace變換代替FFT變換，回避了FFT變換周期性問題。

利用FFT實現(xiàn)數(shù)值Laplace變換荷載，頻域內(nèi)的荷載向量可以表示為［7］

式中P(t)為時域內(nèi)的節(jié)點荷載向量;P(s)為頻域內(nèi)的節(jié)點荷載向量;L為Laplace變換算子;F為FFT變換算子;s=σ+iω為Laplace變換參數(shù);其中σ為Laplace變換的實常數(shù);ω為Laplace變換后對應(yīng)的角頻率。

由式(24)可知，時域內(nèi)節(jié)點荷載向量P(t)的Laplace變換可以利用調(diào)制后的荷載向量P(t)e－σt通過FFT變換實現(xiàn)。

頻域內(nèi)的位移向量通過Laplace逆變換即可得到時域內(nèi)的位移向量。Laplace逆變換過程可以表示為［7］:

式中U(s)為頻域內(nèi)的節(jié)點位移向量;U(t)為時域內(nèi)的節(jié)點位移向量;L－1為 Laplace逆變換算子;F－1為快速傅里葉逆變換(IFFT)。

由式(25)可知，時域內(nèi)的節(jié)點位移向量U(t)的Laplace逆變換可以先將頻域內(nèi)的節(jié)點位移向量U(s)進行IFFT變換，然后再乘以eσt得到。

4.2 整體坐標系下空間桿的波譜剛度矩陣

式(19)給出了承受軸向節(jié)點荷載的平面桿的單元剛度矩陣，考慮到空間桁架桿的節(jié)點位移有三個方向(如圖5所示)，計算時需要對式(19)中桿的波譜單元剛度矩陣進行擴展，擴展后的單元剛度矩陣為:

對式(26)中的波譜單元剛度矩陣進行坐標變換，即可得到整體坐標系下的波譜單元剛度矩陣:

式中RT為文獻［11］中給出的空間桁架的轉(zhuǎn)軸變換矩陣，Ke(x，w)為擴展后的波譜單元剛度矩陣，K(x，w)為整體坐標系下的波譜單元剛度矩陣。

4.3 整體波譜矩陣方程

把式(23)中的單元節(jié)點荷載寫成向量形式:

同理，單元節(jié)點位移也可寫成向量形式:

則桿的波譜單元矩陣方程可以表示為:

得到桿的波譜單元矩陣方程后，采用與有限元法相同的分析思路，在頻域內(nèi)將單元剛度矩陣組裝成整體剛度矩陣，將單元節(jié)點位移向量和單元節(jié)點荷載向量分別組裝成整體節(jié)點位移向量和整體節(jié)點荷載向量，組裝后的整體波譜矩陣方程可以表示為:

4.4 波譜單元法計算流程

波譜單元法中單元剛度矩陣的組裝及矩陣方程的求解都與有限元法的計算過程類似，波譜單元法的計算思路可以用以下的流程圖概括:

圖6 波譜單元法計算流程圖Fig.6 A typical flowchart for spectral element method

5 數(shù)值算例

5.1 數(shù)值算例1

采用波譜單元法和劃分不同單元數(shù)的有限元法計算圖7所示的鋼桁架橋中點P點的豎向位移響應(yīng)，有限元法采用的是Newmark積分［12］。計算時選取的地震波為EI Centro波。Laplace變換的采樣時間取為地震波的時間步長ΔT=0.02 s，采樣點數(shù)N=2 048，計算時長為 T =NΔT =40.96 s，變換實常數(shù) σ=2πNΔT［7］。計算時假定桿的外部粘滯阻尼系數(shù)為 α1=0.01，內(nèi)部粘彈性阻尼系數(shù)為 c=0.000 1。桿的尺寸如圖7所示，截面尺寸及參數(shù)見表1。

表1 鋼桁架橋桿截面尺寸及材料參數(shù)Tab.1 Size and material properties of steel truss bridge

圖7 承受豎向地震荷載的鋼桁架橋Fig.7 Steel truss bridge subjected to vertical seismic load

圖8為波譜單元法和有限元法計算得到的15 s～20 s內(nèi)豎向位移響應(yīng)的對比圖。波譜單元法每根桿劃分1個計算單元，有限元法中每根桿件分別劃分為1個、5個、20個、50個單元。

從圖8(a)可以看出，有限元法每根桿劃為1個單元時的計算結(jié)果與波譜單元法的計算結(jié)果相差很大。但隨著有限元法中劃分單元數(shù)的增加，有限元法計算結(jié)果越來越接近波譜單元法的結(jié)果，如圖8(a)－圖8(d)所示。從圖8(d)可以看出，桁架結(jié)構(gòu)每根桿件劃分為50個單元時有限元法的計算結(jié)果與波譜單元法的結(jié)果已經(jīng)比較接近了。

為便于比較波譜單元法計算結(jié)果與有限元法計算結(jié)果的差值，引進均方根指標(RMSD)進行評價，均方根指標可以定義為［13］:

式中uSEMi為波譜單元法計算的位移;uFEMi為有限元法計算的位移，i表示第i個離散時間點。

為了比較波譜單元法與有限元法的計算速度，在同一計算機，相同的運行環(huán)境下，對兩種方法的計算時間進行測試。表2給出了波譜單元法與劃分不同單元數(shù)的有限元法計算時間及相應(yīng)的均方根指標?？梢钥吹剑S著有限元法劃分單元數(shù)的增加，均方根指標逐漸變小，說明有限元分析的結(jié)果隨著劃分單元數(shù)的增加在逐漸逼近波譜單元法的計算結(jié)果，與傳統(tǒng)有限元法相比波譜單元法具有很高的計算精度。

圖8 波譜單元法與不同單元數(shù)的有限元法計算豎向位移響應(yīng)圖Fig.8 Comparison of vertical displacements calculated by SEM and FEM meshed with different number of elements

由于基于連續(xù)質(zhì)量體系的波譜單元法能大大減少計算單元數(shù)量［2］，因而與傳統(tǒng)有限元法相比，波譜單元法能有效的縮短計算時間。桁架橋在劃分6 800個單元的情況下，采用有限元法計算時間為波譜單元法計算時間的59.36倍，但仍未達到波譜單元法的計算精度。即使與傳統(tǒng)有限元法劃分相同單元數(shù)的情況下，波譜單元法由于計算時無需按時間步長進行積分計算，仍然能減少計算時間。對于本算例中的桁架橋，在劃分單元數(shù)相同的情況下，波譜單元法計算時間約為有限元法計算時間的1/11。

表2 譜單元法與有限元法計算豎向位移結(jié)果對比Tab.2 Comparison of vertical displacements calculated by SEM and FEM

圖9 承受橫向地震荷載的空間桁架Fig.9 Space truss subjected to horizontal seismic load

圖10 波譜單元法與不同單元數(shù)的有限元法計算橫向位移響應(yīng)圖Fig.10 Comparison of horizontal displacements calculated by SEM and FEM meshed with different number of elements

5.2 數(shù)值算例2

分析圖9所示的空間桁架結(jié)構(gòu)在水平地震荷載作用下的響應(yīng)。此處地震荷載仍輸入EI Centro波。波譜單元法采用的Laplace變換參數(shù)及假定的桁架桿的阻尼系數(shù)都與算例1相同。桁架桿截面相關(guān)參數(shù)見表3。

表3 空間桁架桿截面尺寸及材料參數(shù)Tab.3 Size and material properties of Space truss

利用波譜單元法對圖9所示的空間桁架結(jié)構(gòu)進行分析，每根桿件作為1個計算單元。采用有限元法分析時，每根桁架桿分別劃分為1個、3個、5個、20個計算單元。圖10對比了通過波譜單元法和有限元法計算桁架下弦中心O點水平位移前5s響應(yīng)?？梢钥闯?，隨著有限元法中劃分單元數(shù)的增加，有限元法計算結(jié)果越來越接近波譜單元法的結(jié)果。桁架結(jié)構(gòu)每根桿件劃分為20個單元時計算得到的O點水平位移曲線與波譜單元法的計算曲線大致重合，如圖10(d)所示。

表4給出了波譜單元法與劃分不同單元數(shù)的有限元法計算時間及相應(yīng)的均方根指標。與算例1一樣，隨著有限元法劃分單元數(shù)的增加，均方根指標逐漸變小，表示有限元法的計算結(jié)果將隨計算單元數(shù)量增加逐漸逼近譜單元法結(jié)果。對于算例2，劃分相同單元數(shù)時，譜單元法計算時間約為有限元法的1/6，而每根桿劃分150個單元時，有限元法計算的時間約為譜單元法計算時間的101倍。算例2再次有效的說明了采用本文改進后的波譜單元法來計算桁架的地震響應(yīng)，不僅能大大改善計算精度，而且能有效提高計算效率，減少計算時間。

表4 譜單元法與有限元法計算水平位移結(jié)果對比Tab.4 Comparison of horizontal displacements calculated by SEM and FEM

6 結(jié)論

本文基于虛功原理，通過波譜形函數(shù)積分推導(dǎo)了桁架結(jié)構(gòu)在地震荷載作用下等效波譜節(jié)點荷載的顯式表達式。在波譜單元法中利用等效后的節(jié)點荷載對EI-Centro波作用下的空間桁架進行地震響應(yīng)分析。分析結(jié)果表明，通過荷載等效的方式，能十分方便進行編程計算，精確求解結(jié)構(gòu)的地震響應(yīng)。采用基于Laplace變換的波譜單元法能有效回避FFT的周期性問題，無需加上無限長單元即可以精確求解有限長結(jié)構(gòu)。通過修改波譜單元法中的波數(shù)，可以簡便的在波譜單元法中考慮阻尼的影響。與傳統(tǒng)有限元法相比，波譜單元法能大大減少單元數(shù)量并提高計算精度。

［1］Lee U.Vibration analysis of one-dimensional structures using the spectraltransfermatrix method ［J］. Engineering Structures，2000，22:681－690.

［2］Dolyle J F，F(xiàn)arris T N.A spectrally formulated finite element for flexural wave propagation in beams［J］.The international journal of Analytical and Experimental Modal Analysis，1990，5(2):99－107.

［3］Lee U，Kim J H，Leung A Y T.The spectral element Method in Structural dynamics［J］.The Shock and Vibration Digest，2000，32(6):451－465.

［4］Dolyle J F，F(xiàn)arris T N.A spectrally formulated finite element for wave propagation in 3 －D frame structures［J］.The international journal of Analytical and Experimental Modal Analysis，1990，5(4):223 －237.

［5］Dolyle J F.A spectrally formulated finite element for longitudinal wave propagation［J］.The international journal of Analytical and Experimental Modal Analysis.1988，3(1):1－5.

［6］Cho J Y，Go H S，Lee U.Dynamic response of the spectral element model by using the FFT［J］.2007，345－346:845－848.

［7］Igawa H，Komatsu K，Yamaguchi I et al.Wave propagation analysis of frame structures using the spectral element method［J］. Journalofsound and vibration， 2004， 277:1071－1081.

［8］Lee U，Lee J K.Spectral element analysis of the structure under dynamic distributed loads［J］.Journal of Mechanical Science and Technology，1998，12(4):565－571.

［9］Lee U，Cho J Y.FFT-based spectral element analysis for the linear continuum dynamic systems subjected to arbitrary initial conditions by using the pseudo-force method ［J］.International Journal for Numerical Methods in Engineering，2007，74(1):159－174.

［10］R.克拉夫，J.彭津.結(jié)構(gòu)動力學［M］.王光遠 等譯.北京:高等教育出版社，2006.

［11］劉樹棠.桿系結(jié)構(gòu)有限元分析與matlab應(yīng)用［M］.北京:中國水利水電出版社，2007.

［12］王元漢，李麗娟，李銀平.有限元法基礎(chǔ)與程序設(shè)計［M］.廣州:華南理工大學出版社，2001.

［13］ Soh C K，Tseng K K H，Bhalla S，et al.Performance of smart piezoceramic patches in health monitoring of a RC bridge［J］.Smart Materials and Structures，2000，9(4)，533－542.

Application of spectral element method in dynamic analysis of a space truss subjected to seismic load

ZHANG Jun-bing1，2，ZHU Hong-ping1，2，WANG Dan-sheng1，2，GE Dong-dong1，2

(1.School of Civil Engineering＆ Mechanics，Huazhong University of Science＆ Technology，Wuhan 430074，China;2.Hubei Key Laboratory of Control Structure，Huazhong University of Science＆ Technology，Wuhan 430074，China)

An extended spectral element method(SEM)was established to get dynamic responses of a space truss subjected to seismic load.The seismic load was equivalent to concentrated node forces by integrating the shape function in SEM based on the principle of virtual work.Both internal viscoelastic damping and external viscous damping of the truss bar were considered by just simply modifying the wave number.Laplace transformation instead of fast Fourier transformation(FFT)was utilized in SEM to avoid the periodicity of FFT.To evaluate the accuracy of Laplace-based SEM，the dynamic responses of a space truss under seismic load was analyzed as a numerical example.The numerical results obtained using FEM were compared with those using SEM.It was found that SEM provides good dynamic results under seismic load;the equivalent seismic node forces are convenient to get by programming calculation;SEM is proved to be an efficient method to analyze dynamic responses of structures accurately while the number of element greatly decreases.

spectral element method(SEM);damping;wave number;space truss;FEM;seismic responses

TU311.3;TU323.4

A

國家杰出青年基金(50925828);國家自然科學基金資助項目(50778077);高校博士學科點專項科研基金(20070487099)

2009－12－18 修改稿收到日期:2010－03－05

張俊兵 男，博士，1985年2月生