邱志平, 邱 宇

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

1 引 言

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

(4)

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

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

(5)

稱(chēng)方程(3)為線(xiàn)性自治Birkhoff方程,其中G為對(duì)稱(chēng)的常數(shù)矩陣。線(xiàn)性自治Birkhoff方程的形式可以表示為

(6)

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

(7)

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

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

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

(8)

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

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

(9)

(10)

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

(11)

顯然其是線(xiàn)性自治Birkhoff方程的形式。

(12)

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

(13)

注意其中A為反對(duì)稱(chēng)矩陣。

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

B=B0+εB1

(14)

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

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

(15)

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

(16)

對(duì)于0階攝動(dòng)方程,可以轉(zhuǎn)化為一個(gè)線(xiàn)性自治Birkhoff方程

(17)

同樣地,可以將剩余的攝動(dòng)方程通過(guò)變量變換的形式化為一系列線(xiàn)性自治Birkhoff方程,即

(18)

4 線(xiàn)性自治Birkhoff方程的中點(diǎn)格式

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

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

(19)

式中τ為時(shí)間步長(zhǎng),zk和zk+1為第k和k+1個(gè)時(shí)間步z的值。對(duì)式(19)整理可得

(20)

那么步進(jìn)映射zkzk+1的Jacobi矩陣為

(21)

對(duì)于矩陣P=G-1K,顯然有

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

KG-1K-KG-1K=0

(22)

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

(23)

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

(24)

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

(25)

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

(26)

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

5 數(shù)值算例

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

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

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

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

(27)

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

(28)

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

(29)

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

(30)

那么方程形式為

(31)

式中A為反對(duì)稱(chēng)矩陣,其形式為

(32)

B為對(duì)稱(chēng)矩陣,其具體形式為

(33)

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

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

(34)

(35)

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

進(jìn)一步檢驗(yàn)算法在長(zhǎng)時(shí)間計(jì)算下的表現(xiàn),仍取步長(zhǎng)為0.1 s,取仿真時(shí)間為20 s。本文的保辛中點(diǎn)格式的結(jié)果與解析解的對(duì)比如圖4所示。可以看出,在長(zhǎng)時(shí)間計(jì)算下,本文的保辛中點(diǎn)格式仍與解析解保持一致,這說(shuō)明了算法的高精度。

圖2 四階龍格庫(kù)塔方法與解析解對(duì)比

圖4 長(zhǎng)時(shí)間計(jì)算下本文的保辛中點(diǎn)格式與解析解對(duì)比

5.2 兩室建筑振動(dòng)問(wèn)題

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

圖5 兩室建筑

該系統(tǒng)運(yùn)動(dòng)方程的矩陣形式為

(36)

式中

(37)

(38)

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

采用本文保辛中點(diǎn)格式求解,并與四階龍格庫(kù)塔方法比較,取時(shí)間步長(zhǎng)為0.2 s,仿真時(shí)間為2 s。取四階龍格庫(kù)塔方法在2×10-4s步長(zhǎng)下的結(jié)果作為參考解,對(duì)比算法的精度,如圖6所示。可以看出,對(duì)于非保守系統(tǒng),本文的保辛中點(diǎn)格式仍具有高精度。

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

圖6 本文的保辛中點(diǎn)格式與四階龍格庫(kù)塔和參考解的比較

6 結(jié) 論

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

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