邱志平, 邱 宇

(1.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院,北京 100191;2.北京航空航天大學(xué) 沈元學(xué)院,北京 100191)

1 引 言

結(jié)構(gòu)動響應(yīng)預(yù)測是工程中非常重要的問題,同時是結(jié)構(gòu)振動控制和載荷識別的基礎(chǔ)[1]。結(jié)構(gòu)動響應(yīng)預(yù)測主要依靠的是動力學(xué)分析的相關(guān)方法,而傳統(tǒng)的結(jié)構(gòu)動力學(xué)分析主要基于Lagrange體系,以變分原理為基石[2],且基于此已經(jīng)發(fā)展出了成熟且應(yīng)用廣泛的分析方法。Hamilton對經(jīng)典力學(xué)重新表述并基于此建立了Hamilton力學(xué)體系,相較于Lagrange體系其具有明顯的辛結(jié)構(gòu)[3]但是不具有對應(yīng)的變分原理。由Birkhoff[4]提出的Birkhoff系統(tǒng)是Hamilton系統(tǒng)最自然的一般推廣,同時擁有辛結(jié)構(gòu)和相應(yīng)的變分原理[5]。其考慮耗散項,相比Hamilton系統(tǒng)能夠涵蓋更多的實際力學(xué)系統(tǒng)。因此,利用Birkhoff系統(tǒng)的相關(guān)方法和原理求解結(jié)構(gòu)動響應(yīng)問題具有重要的意義。

目前,Birkhoff方程的求解方法大多從Hamilton算法[6-9]推廣而來。張興武等[10]將哈密頓辛差分格式推廣至自治Birkhoff系統(tǒng),利用Cayley變換構(gòu)造了自治Birkhoff系統(tǒng)的歐拉中點等辛差分格式?；诠茴D系統(tǒng)的生成函數(shù)法,蘇紅玲等[11]提出了構(gòu)造Birkhoff系統(tǒng)辛算法的生成函數(shù)法。近年來,文獻[12-14]基于Paff-Birkhoff變分原理,通過離散方法直接獲得了約束Birkhoff系統(tǒng)等的Birkhoff辛算法。

然而,以上方法仍然存在局限性。自治Birkhoff系統(tǒng)的歐拉中點格式要求反對稱系數(shù)矩陣非奇異,因而無法適用于奇數(shù)維廣義Birkhoff系統(tǒng)。生成函數(shù)法存在構(gòu)造的困難,實際應(yīng)用存在障礙。離散的Birkhoff方法存在離散誤差,使得解存在波動性。

本文針對結(jié)構(gòu)動響應(yīng)問題,提出一種Birkhoff形式下的保辛中點格式。首先將保守和非保守情形下的結(jié)構(gòu)動響應(yīng)方程化為線性自治的Birkhoff方程。之后對該線性自治方程進行中心差分,經(jīng)過推導(dǎo)得到其中點格式。該中點格式不需要對Birkhoff方程中的反對稱矩陣求逆,從而不要求該矩陣非奇異,因此適用于奇數(shù)維廣義Birkhoff系統(tǒng)。最后通過兩個數(shù)值算例說明了本文方法的有效性。

2 Birkhoff方程及其辛結(jié)構(gòu)

Birkhoff方程是哈密頓方程的自然推廣,其一般形式為[14]

(1)

式中z=(z1,z2,…,zm)T∈m,m可以是偶數(shù),也可以是奇數(shù)。當(dāng)m為奇數(shù)時,式(1)稱為奇數(shù)維Birkhoff方程[16]。式(1)中B(z,t)和F(z,t)統(tǒng)稱為Birkhoff函數(shù),特別稱B(z,t)為Birkhoff量。K為一個反對稱矩陣,其元素滿足

(2)

當(dāng)Birkhoff函數(shù)B和F都不顯含時間變量t時,稱Birkhoff方程(1)為自治的,其形式為

(3)

當(dāng)Birkhoff函數(shù)F不顯含時間t,而B顯含時間t時,稱Birkhoff方程(1)為半自治的,其形式為

(4)

當(dāng)Birkhoff函數(shù)B和F都顯含時間t時,此時稱為非自治的,即為式(1)的形式。

對于自治Birkhoff方程(3),當(dāng)系數(shù)矩陣K為常數(shù)矩陣,且Birkhoff函數(shù)B為變量z的二次型,即

(5)

稱方程(3)為線性自治Birkhoff方程,其中G為對稱的常數(shù)矩陣。線性自治Birkhoff方程的形式可以表示為

(6)

由Kij可以定義一個2-形式,用局部坐標(biāo)表示為

(7)

閉2-形式(7)包括了自治和半自治Birkhoff方程的幾何結(jié)構(gòu)[14]。

3 結(jié)構(gòu)動響應(yīng)問題的線性自治Birkhoff形式

結(jié)構(gòu)動響應(yīng)問題的控制方程可以表示為

(8)

式中M為質(zhì)量矩陣,D為阻尼矩陣,S為剛度矩陣,其都是對稱矩陣。

首先考慮保守情況,即D=0,且不存在外激勵載荷F(t)。此時動響應(yīng)方程的形式為

(9)

(10)

將其系數(shù)矩陣記為A和C,那么方程(10)化為

(11)

顯然其是線性自治Birkhoff方程的形式。

(12)

將其系數(shù)矩陣分別記為A,B和N,方程(12)化為

(13)

注意其中A為反對稱矩陣。

利用攝動方法[17],可將方程(13)化為一系列線性自治Birkhoff方程[17]。引入小參數(shù)ε,將矩陣B表示為

B=B0+εB1

(14)

式中B0=(BT+B)/2,B1=(BT-B)/2。將自變量z表示為攝動級數(shù)展開形式

z=z0+z1ε+z2ε2+…

(15)

將式(14,15)代入式(13)可得一系列攝動方程

(16)

對于0階攝動方程,可以轉(zhuǎn)化為一個線性自治Birkhoff方程

(17)

同樣地,可以將剩余的攝動方程通過變量變換的形式化為一系列線性自治Birkhoff方程,即

(18)

4 線性自治Birkhoff方程的中點格式

對于形如式(6)的線性自治Birkhoff方程,對其進行中心差分,可得

K(zk+1-zk)/τ=G(zk+1+zk)

(19)

式中τ為時間步長,zk和zk+1為第k和k+1個時間步z的值。對式(19)整理可得

(20)

那么步進映射zkzk+1的Jacobi矩陣為

(21)

對于矩陣P=G-1K,顯然有

KP+PTK=KG-1K+KTG-TK=

KG-1K-KG-1K=0

(22)

因此P為無窮小辛矩陣,同理乘上常數(shù)系數(shù)后2G-1K/τ仍為無窮小辛矩陣,而Φτ為無窮小辛矩陣2G-1K/τ的Cayley變換[10]。于是有

(23)

線性自治Birkhoff方程的辛結(jié)構(gòu)為式(7)的形式,其矩陣形式為

(24)

對于離散的線性自治Birkhoff方程,第k和k+1兩個時間步的辛結(jié)構(gòu)之差為

(25)

將式(21)代入式(25),并利用關(guān)系式(23),可得

(26)

因此,本文提出的中點格式(20)是保辛的。

5 數(shù)值算例

5.1 三自由度彈簧質(zhì)點系統(tǒng)

首先將本文方法應(yīng)用于一個三自由度的彈簧質(zhì)點系統(tǒng),其是一個保守系統(tǒng),如圖1所示。

圖1 三自由度彈簧質(zhì)點系統(tǒng)

每個質(zhì)量塊的質(zhì)量均為m=1 kg,而彈簧的剛度系數(shù)具有較大的分散性,用以驗證算法的穩(wěn)定性,其分別為k1=1000 N/m,k2=100 N/m,k3=10 N/m,k4=1 N/m。該系統(tǒng)運動方程的矩陣形式為

(27)

式中x=(x1,x2,x3)T,質(zhì)量矩陣M和剛度矩陣S分別為

(28)

將初始條件設(shè)置為

(29)

對于式(27)的保守系統(tǒng),可以直接引入新變量化為線性子自治Birkhoff方程。令

(30)

那么方程形式為

(31)

式中A為反對稱矩陣,其形式為

(32)

B為對稱矩陣,其具體形式為

(33)

首先檢驗算法在大步長下的表現(xiàn),取步長為0.1 s,仿真時間為2 s,對比不同算法在x1位移求解的表現(xiàn)。x1位移解析解表達式為

x1(t)=-0.0990cos(ω1t)+0.0979cos(ω2t)+0.0011cos(ω3t)

(34)

(35)

在此步長下,四階龍格庫塔方法很快就發(fā)散了,如圖2所示。而本文所提保辛中點格式與解析解吻合良好,如圖3所示。這說明了本文的保辛中點格式具有良好的穩(wěn)定性,即便在大步長下求解大剛度差異系統(tǒng)響應(yīng),仍能保持高精度。

進一步檢驗算法在長時間計算下的表現(xiàn),仍取步長為0.1 s,取仿真時間為20 s。本文的保辛中點格式的結(jié)果與解析解的對比如圖4所示?？梢钥闯?在長時間計算下,本文的保辛中點格式仍與解析解保持一致,這說明了算法的高精度。

圖2 四階龍格庫塔方法與解析解對比

圖4 長時間計算下本文的保辛中點格式與解析解對比

5.2 兩室建筑振動問題

一些大型建筑,如機場航站樓和大型倉庫等,在研究這些建筑的內(nèi)部壓力時,可以將其簡化為具有開口的兩室建筑。下面考慮這樣一個兩室建筑,如圖5所示,其是一個非保守系統(tǒng)。

圖5 兩室建筑

該系統(tǒng)運動方程的矩陣形式為

(36)

式中

(37)

(38)

各個參數(shù)的取值列入表1[19]。

采用本文保辛中點格式求解,并與四階龍格庫塔方法比較,取時間步長為0.2 s,仿真時間為2 s。取四階龍格庫塔方法在2×10-4s步長下的結(jié)果作為參考解,對比算法的精度,如圖6所示?？梢钥闯?對于非保守系統(tǒng),本文的保辛中點格式仍具有高精度。

表1 兩室建筑參數(shù)取值

圖6 本文的保辛中點格式與四階龍格庫塔和參考解的比較

6 結(jié) 論

針對結(jié)構(gòu)動響應(yīng)問題,本文提出了一種Birkhoff形式下的保辛中點格式。這種方法相較于傳統(tǒng)方法,具有精度高和穩(wěn)定性強的特點。相較于其他隱式的保辛離散方法,本文的顯式方法需要更少的迭代步驟和計算資源,因此在實際應(yīng)用中具有更高的求解效率。此外,由于算法中無需對反對稱系數(shù)矩陣求逆,因此也適用于奇數(shù)維Birkhoff方程。兩個數(shù)值算例表明了算法的有效性,相比于四階龍格庫塔更具有優(yōu)勢。

需要注意的是,本文方法是在線性Birkhoff 系統(tǒng)的框架內(nèi)提出的。對于非線性系統(tǒng),由于其非線性系數(shù)矩陣的存在,本文方法無法直接應(yīng)用。然而借助一些變換方法(如攝動法等),存在將本文方法應(yīng)用于非線性Birkhoff系統(tǒng)的可能性。