邱志平 姜 南

(北京航空航天大學(xué)航空科學(xué)與工程學(xué)院,北京 100191)

引言

隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展,研究人員越來越追求更高效、更穩(wěn)定和長(zhǎng)時(shí)間模擬能力更強(qiáng)的數(shù)值算法[1].同一力學(xué)問題有牛頓力學(xué)、拉格朗日力學(xué)、哈密頓力學(xué)體系三種表示形式,其中一切可忽略耗散的物理過程都可以表示為某種哈密頓形式,但由于求解途徑不同,產(chǎn)生的計(jì)算結(jié)果可能是不等效的[2].傳統(tǒng)算法除少數(shù)例外,都不是辛算法,不可避免地帶有人為耗散性等歪曲體系特征的缺陷,用于短期模擬尚可,用于長(zhǎng)期跟蹤則會(huì)導(dǎo)致結(jié)果嚴(yán)重失真[3-4].而哈密頓系統(tǒng)辛算法卻有保持體系結(jié)構(gòu)的優(yōu)點(diǎn),在結(jié)構(gòu)對(duì)稱性和守恒性方面優(yōu)于傳統(tǒng)算法,特別在穩(wěn)定性和長(zhǎng)期跟蹤能力上具有獨(dú)特優(yōu)越性[2-5].

哈密頓系統(tǒng)辛算法的思想始于20 世紀(jì)80 年代Ruth[6]和馮康[7]的工作.1983 年,Ruth[6]構(gòu)造了可分哈密頓系統(tǒng)的前三階辛差分格式.1984 年,馮康[7]首次系統(tǒng)地提出了哈密頓系統(tǒng)的辛算法,開創(chuàng)了哈密頓力學(xué)數(shù)值計(jì)算的新領(lǐng)域.1988 年,Sanz-Serna[8]、Lasagni[9]和Suris[10]分別從不同角度給出了辛Runge-Kutta 方法的判定條件.1993 年,Sun[11]對(duì)一般的哈密頓系統(tǒng)構(gòu)造了辛分塊Runge-Kutta 方法.2000 年前后,Bridges[12]、Reich[13]和Marsden等[14]構(gòu)建了多辛算法的理論框架.之后,變分積分子、離散梯度法、投影算法、分裂組合算法、對(duì)稱算法等[15]一系列方法不斷提出并發(fā)展.近年來,無網(wǎng)格格式[16]、連續(xù)級(jí)Runge-Kutta 方法[17]等辛算法也相繼被提出.深入的理論分析和大量的數(shù)值實(shí)驗(yàn)令人信服地表明,辛算法在數(shù)值計(jì)算中具有顯著優(yōu)越性.

然而,動(dòng)力系統(tǒng)中不可避免地存在大量的、不同程度的不確定性[18].例如,工程結(jié)構(gòu)中的典型不確定性有:建模過程中的簡(jiǎn)化操作導(dǎo)致所建模型存在誤差,制造環(huán)境、材料多相特征等因素使彈性模量、泊松比等材料參數(shù)具有分散性,制造及安裝誤差使結(jié)構(gòu)幾何尺寸具有不確定性,測(cè)量條件、外部環(huán)境等因素使外載荷具有不確定性等[19-20].這些不確定性將引起動(dòng)力響應(yīng)變化,響應(yīng)不確定量有時(shí)甚至能達(dá)到參數(shù)不確定量的數(shù)倍[21].此外,不確定性的存在還會(huì)使系統(tǒng)的性質(zhì)發(fā)生改變,保守系統(tǒng)若存在不確定性則不再保守.根據(jù)產(chǎn)生機(jī)理不同,不確定性可分為隨機(jī)不確定性和認(rèn)知不確定性[22-23]:隨機(jī)不確定性是由自然變異和隨機(jī)性而導(dǎo)致的不確定性,以隨機(jī)模型定量化;認(rèn)知不確定性是指受知識(shí)水平和社會(huì)環(huán)境等因素制約而產(chǎn)生的認(rèn)知上的不確定性,通常以區(qū)間模型定量化.因此,需要在哈密頓系統(tǒng)中考慮隨機(jī)和區(qū)間不確定性的影響,以確保動(dòng)力學(xué)分析計(jì)算的合理有效性.

1984 年,Zambrini 首先基于變分原理提出了Nelson 隨機(jī)力學(xué)在哈密頓力學(xué)下的運(yùn)動(dòng)方程,逐步建立了哈密頓力學(xué)下的隨機(jī)力學(xué)體系.2002 年,Milstein 等[24]對(duì)一般的隨機(jī)哈密頓系統(tǒng)構(gòu)造了幾類辛Runge-Kutta 方法,開創(chuàng)了隨機(jī)哈密頓系統(tǒng)辛算法這一全新的研究領(lǐng)域.2007 年,王麗瑾[25]提出了隨機(jī)變分積分子理論,使利用隨機(jī)生成函數(shù)構(gòu)造隨機(jī)辛算法成為可能.2009 年,Bou-Rabee 和Owhadi[26]對(duì)隨機(jī)哈密頓系統(tǒng)構(gòu)造了隨機(jī)變分積分子.近幾年,丁效華課題組[27-28]也對(duì)隨機(jī)哈密頓系統(tǒng)進(jìn)行了研究,并構(gòu)造了隨機(jī)辛Runge-Kutta、辛分塊Runge-Kutta 方法等.此外,朱位秋[29-30]研究了多自由度非線性隨機(jī)動(dòng)力學(xué)系統(tǒng),在國際上首次提出了隨機(jī)激勵(lì)的耗散的哈密頓系統(tǒng)理論,構(gòu)建了非線性隨機(jī)動(dòng)力學(xué)與控制的哈密頓理論框架.上述隨機(jī)哈密頓系統(tǒng)相關(guān)研究主要聚焦隨機(jī)白噪聲激勵(lì),但沒有考慮系統(tǒng)本身的參數(shù)隨機(jī)性,并且考慮區(qū)間不確定性的哈密頓系統(tǒng)也未見有人研究.

哈密頓系統(tǒng)辛算法的工程結(jié)構(gòu)應(yīng)用方面,羅恩等[31]建立了多自由度系統(tǒng)彈性動(dòng)力學(xué)的相空間非傳統(tǒng)哈密頓變分原理,提出了稱之為辛?xí)r間子域法的辛算法,精度和效率都具有明顯優(yōu)越性.邢譽(yù)峰課題組[32-33]針對(duì)結(jié)構(gòu)動(dòng)力學(xué)方程構(gòu)造了多種辛差分格式,得到了令人滿意的結(jié)果.高強(qiáng)等[34]用辛算法求解拉壓剛度不同桁架的非線性動(dòng)力問題,Li 等[35]構(gòu)造了一類動(dòng)力學(xué)初值問題的辛算法并應(yīng)用于簡(jiǎn)諧振子和簡(jiǎn)支梁,Yang 等[36]用辛算法對(duì)超長(zhǎng)細(xì)桿進(jìn)行動(dòng)力學(xué)數(shù)值模擬.也有學(xué)者對(duì)多辛算法的工程結(jié)構(gòu)應(yīng)用開展了相關(guān)研究,如梁[37]、板[38]、桿[39]的動(dòng)力響應(yīng)等.然而,用哈密頓系統(tǒng)求解含不確定性的結(jié)構(gòu)動(dòng)力響應(yīng)還少有人研究.

本文針對(duì)隨機(jī)和區(qū)間不確定性,對(duì)含參數(shù)不確定性的非齊次線性哈密頓系統(tǒng)的動(dòng)力響應(yīng)分析進(jìn)行首次嘗試,提出兩種不確定性哈密頓系統(tǒng)的參數(shù)攝動(dòng)法,突破傳統(tǒng)哈密頓系統(tǒng)只適用于保守系統(tǒng)的局限性,并開展兩種哈密頓系統(tǒng)不確定性響應(yīng)的相容性研究,以期為結(jié)構(gòu)動(dòng)力響應(yīng)評(píng)估提供更加有效穩(wěn)定的數(shù)值計(jì)算方法.

1 確定性非齊次線性哈密頓系統(tǒng)

考慮哈密頓正則方程

其中,In是n階單位矩陣;J是標(biāo)準(zhǔn)單位辛矩陣,滿足J?1=JT=?J;H稱為該系統(tǒng)的哈密頓函數(shù).

哈密頓系統(tǒng)(1)稱為線性的,如果哈密頓函數(shù)H是z的二次型

其中C為對(duì)稱矩陣,即CT=C,于是正則方程(1)能夠表示為

其中B=J?1C是無窮小辛陣,即滿足JB+BTJ=0.

非齊次線性哈密頓系統(tǒng)是在線性哈密頓系統(tǒng)基礎(chǔ)上增加了非齊次項(xiàng)

其中F(t)是與時(shí)間t有關(guān)的向量.

此時(shí),哈密頓函數(shù)H′可以表示為

從而,非齊次線性哈密頓系統(tǒng)(5)可以表示為哈密頓正則方程(1)的形式.

對(duì)于哈密頓系統(tǒng),其相空間配備著一個(gè)標(biāo)準(zhǔn)辛結(jié)構(gòu),也就是一個(gè)閉的微分2 形式

哈密頓系統(tǒng)相流保持相空間的辛結(jié)構(gòu)不變,即

2 含擾動(dòng)非齊次線性哈密頓系統(tǒng)的參數(shù)攝動(dòng)法

假設(shè)非齊次線性哈密頓系統(tǒng)(5)中的矩陣B和向量F(t)中元素與系統(tǒng)參數(shù)b=(b1,b2,···,bm)T有關(guān),此時(shí),哈密頓系統(tǒng)(5)可以寫為

當(dāng)系統(tǒng)參數(shù)b=(b1,b2,···,bm)T存在擾動(dòng)時(shí),矩陣B(b)和向量F(b)由標(biāo)稱值變化到擾動(dòng)系統(tǒng),分析擾動(dòng)對(duì)哈密頓系統(tǒng)響應(yīng)z的影響.

根據(jù)攝動(dòng)理論,引入小參數(shù)ε、參數(shù)b、矩陣B、向量F和哈密頓系統(tǒng)的解z可以分別寫為攝動(dòng)級(jí)數(shù)展開形式

其中,ε 是小于單位1 的小量,b0,B0,F0和z0分別是b,B,F和z的標(biāo)稱值,bi,Bi,Fi和zi(i=1,2,···,p)分別是其第i階攝動(dòng)量.含ε 的項(xiàng)表示其與標(biāo)稱值相比是一個(gè)很小的量,在這種情況下,哈密頓系統(tǒng)的解z只有小變化.通過選擇使|ε| 足夠小的參數(shù)可以使級(jí)數(shù)收斂.

將式(10)代入式(9),得

將式(11)展開,比較ε 的同次冪系數(shù),可得

運(yùn)用辛算法求解式(12)中第1 式可以得到z的標(biāo)稱部分,即z0.本文中采用的辛算法均為歐拉中點(diǎn)格式.對(duì)式(12)中其他式子的Bi和Fi(i=1,2,···,p),可以通過Taylor 展開獲得,即將B和F在b=b0處分別進(jìn)行Taylor 展開,得到

其中bij(i=1,2,···,p;j=1,2,···,m)是bi的分量.

由式(10)和式(13),可知

將式(12)中第1 式求得的z0和式(14)求得的B1,F1代入式(12)中第2 式并利用辛算法求解,可以求得z的第1 階攝動(dòng)量z1;進(jìn)而可以通過依次求解式(12)中其他式子得到zi(i=2,3,···,p)的值.從而,z就可以按式(10)求得.在每一步都采用辛算法求解保證計(jì)算結(jié)果能夠保持體系結(jié)構(gòu)特征,避免傳統(tǒng)算法帶有的人為耗散性等歪曲體系特征的缺陷.在實(shí)際計(jì)算中,為了方便求解,常常展開到第1 階攝動(dòng).

當(dāng)參數(shù)b存在的擾動(dòng)為隨機(jī)或區(qū)間不確定性時(shí),上述參數(shù)攝動(dòng)法可以推廣至求解隨機(jī)或區(qū)間非齊次線性哈密頓系統(tǒng),詳細(xì)求解過程如下面第3、4 節(jié)所述.

3 隨機(jī)非齊次線性哈密頓系統(tǒng)的參數(shù)攝動(dòng)法

當(dāng)系統(tǒng)參數(shù)b=(b1,b2,···,bm)T是隨機(jī)變量時(shí),矩陣B、向量F和哈密頓系統(tǒng)的解z也是隨機(jī)的,它們可以分別看作圍繞確定性部分(均值)有一個(gè)隨機(jī)小擾動(dòng).因此,基于前述攝動(dòng)理論,同樣引入小參數(shù)ε,將b,B,F和z分別表示為

其中,bd,Bd,Fd和zd分別是b,B,F和z的確定性部分;bs,Bs,Fs和zs分別是其隨機(jī)部分,且它們的均值均為0.

將式(15)代入式(9),得

將式(16)展開,忽略O(shè)(ε2)高階項(xiàng),并比較ε 的同次冪系數(shù),可得

利用辛算法求解式(17)中第1 式可以求得z的確定性部分zd,即為z的均值.但無法由式(17)直接確定z的隨機(jī)部分zs,需要進(jìn)行變換后加以求解.

對(duì)式(15)求取數(shù)學(xué)期望,有

由于zd是一個(gè)確定性的向量,所以zd與z,zd與zs均相互獨(dú)立,從而可得

因此,z的協(xié)方差矩陣可以寫為

式(20)中的E[zzT]可以寫為

將式(21)代入式(20),可得協(xié)方差矩陣為

該協(xié)方差矩陣的對(duì)角元素表示各點(diǎn)的方差,其他非對(duì)角元素表示各點(diǎn)間的協(xié)方差.因此,z的各分量zi(i=1,2,···,2n)的方差為

同理,b的各分量bj(j=1,2,···,m)的方差為

將z在z=zd處進(jìn)行一階Taylor 展開得到

其中,bsj(j=1,2,···,m)是bs的分量.

從而,z的隨機(jī)部分zs可以表示為

同理,Bs和Fs可以表示為

將式(26)～式(27)代入式(17)中第2 式,得

其中Cov(bk,bl)表示bk和bl的協(xié)方差.

若參數(shù)bj是相互獨(dú)立的,Cov(bk,bl)為0,則式(29)可以簡(jiǎn)化為

將式(24)代入式(30),就可以得到z的各分量zi(i=1,2,···,2n)的方差

因此,zi(i=1,2,···,2n)的標(biāo)準(zhǔn)差為

其中σ[bj]是bj(j=1,2,···,m)的標(biāo)準(zhǔn)差.

在計(jì)算求解過程中,求解式(17)中第1 式求得z的均值z(mì)d和求解式(28)求得z的隨機(jī)部分zs都采用了辛算法,確保計(jì)算結(jié)果能夠保持體系結(jié)構(gòu)特征.

令k為正整數(shù),則對(duì)zi(i=1,2,···,2n)而言,距其均值的距離為k倍標(biāo)準(zhǔn)差的范圍為

其中范圍的上界和下界分別為

4 區(qū)間非齊次線性哈密頓系統(tǒng)的參數(shù)攝動(dòng)法

當(dāng)系統(tǒng)參數(shù)b=(b1,b2,···,bm)T是區(qū)間變量時(shí),即參數(shù)b在一個(gè)區(qū)間向量?jī)?nèi)取值,矩陣B、向量F和哈密頓系統(tǒng)的解z也分別在一個(gè)區(qū)間范圍內(nèi)取值

其中,bI,FI和zI是區(qū)間向量,BI是區(qū)間矩陣,分別是分別是其上界.

此時(shí),哈密頓方程(9)可寫為

式(36)是區(qū)間非齊次線性哈密頓方程.

bI,BI,FI和zI的中值和半徑分別為

利用區(qū)間中心表示法,bI,BI,FI和zI可以分別表示為

其中

如果將?bI,?BI,?FI和?zI分別看作圍繞bc,Bc,Fc和zc的擾動(dòng),則可以采用第2 節(jié)所述參數(shù)攝動(dòng)法求解區(qū)間哈密頓方程(40).

按照區(qū)間的含義,引入小參數(shù)ε,式(40)可以表示為:由于參數(shù)b存在小擾動(dòng)δb,導(dǎo)致B,F和z產(chǎn)生小擾動(dòng)δB,δF和δz,且滿足

條件下的擾動(dòng)方程的形式

式(41)和式(42)所表示的問題可以理解為:在參數(shù)中值bc已知,從而能夠確定中值Bc和Fc,而小擾動(dòng)δb的具體取值未知但其取值范圍(41)已知,小擾動(dòng)δB和δF的具體取值也未知但其取值范圍(41)可以確定的情況下,確定哈密頓系統(tǒng)的解z的界限.展開式(42)并比較ε 的同次冪系數(shù),可得

運(yùn)用辛算法求解式(43)中第1 式可以求得z的中值z(mì)c.由區(qū)間擴(kuò)張,式(43)中第2 式可寫為

將z在z=zc處進(jìn)行一階Taylor 展開得到

其中,δbj(j=1,2,···,m)是δb的分量.

從而可得δz的表達(dá)式

式(46)和式(47)的區(qū)間擴(kuò)張形式為

其中,?bj(j=1,2,···,m)是?b的分量.

將式(48)代入式(44),得

在計(jì)算求解過程中,求解式(43)中第1 式求得z的中值z(mì)c和求解式(49)求得z的半徑?z都采用了辛算法,同樣確保計(jì)算結(jié)果能夠保持體系結(jié)構(gòu)特征.

由區(qū)間相等的定義可得z的上界和下界分別為

z的分量形式zi(i=1,2,···,2n)的上界和下界分別為

5 隨機(jī)和區(qū)間非齊次線性哈密頓系統(tǒng)比較

不確定性是客觀存在的,無論是隨機(jī)方法還是區(qū)間方法,只是描述不確定性的不同形式,不能從本質(zhì)上改變不確定性對(duì)響應(yīng)的影響規(guī)律,利用兩種方法得到的不確定性響應(yīng)理應(yīng)具有相容性.因此,本節(jié)對(duì)隨機(jī)和區(qū)間非齊次線性哈密頓系統(tǒng)的分析結(jié)果進(jìn)行比較,探究?jī)烧唔憫?yīng)界限的包含關(guān)系.

假定參數(shù)b的區(qū)間范圍由概率統(tǒng)計(jì)信息獲取,即可以表示為距其均值距離為k倍標(biāo)準(zhǔn)差的形式

其中,k為正整數(shù),σ[b]是b的標(biāo)準(zhǔn)差,b的上界和下界分別為

由式(54)～式(55),可以得到參數(shù)b的區(qū)間中值和半徑與隨機(jī)確定性部分和標(biāo)準(zhǔn)差之間存在關(guān)系

哈密頓系統(tǒng)的解z的區(qū)間中值與隨機(jī)確定性部分同樣存在關(guān)系

將式(57)的分量形式代入式(53),得到zi(i=1,2,···,2n)的上界和下界

對(duì)于參數(shù)bj,切比雪夫不等式成立

由不等式(61),可以得到由隨機(jī)方法確定的上下界(34)和由區(qū)間方法確定的上下界(59)存在關(guān)系

式(62)表示,在由概率統(tǒng)計(jì)信息確定不確定性參數(shù)的區(qū)間范圍的情況下,對(duì)于不確定性非齊次線性哈密頓系統(tǒng),由區(qū)間方法獲得的哈密頓系統(tǒng)的解的范圍比由隨機(jī)方法獲得的范圍大,即區(qū)間方法得到的上界比隨機(jī)方法得到的上界大,而區(qū)間方法得到的下界比隨機(jī)方法得到的下界小.

6 數(shù)值算例

為了驗(yàn)證所提方法在結(jié)構(gòu)動(dòng)力響應(yīng)中的可行性和有效性,本節(jié)給出兩個(gè)數(shù)值算例,包括懸臂梁和復(fù)合材料層合板,并將本文所提隨機(jī)、區(qū)間方法(分別簡(jiǎn)記為SHPM、IHPM)計(jì)算結(jié)果與傳統(tǒng)隨機(jī)、區(qū)間方法(分別簡(jiǎn)記為TSM、TIM)計(jì)算結(jié)果相比較.

6.1 懸臂梁

考慮如圖1 所示11 節(jié)點(diǎn)、10 單元懸臂梁在正弦激勵(lì)作用下的動(dòng)力響應(yīng).梁長(zhǎng)L=1 m,橫截面積A=2 cm2,橫截面的慣性矩Iz=2 cm4,材料泊松比ν=0.3.正弦激勵(lì)P(t)=?psin(1600πt)N 作用在節(jié)點(diǎn)3 的豎直方向,初始條件為˙x(0)=0,x(0)=0.

圖1 11 節(jié)點(diǎn)、10 單元懸臂梁Fig.1 A cantilever beam with 11 nodes,10 elements

由于材料中不可避免的分散性及測(cè)量誤差,材料的彈性模量、密度和正弦激勵(lì)幅值均具有不確定性,假設(shè)它們所在的區(qū)間范圍I如表1 所示;同時(shí)假設(shè)它們?cè)趨^(qū)間范圍內(nèi)服從正態(tài)分布,均值μ 和標(biāo)準(zhǔn)差σ 也如表1 所示.

表1 材料彈性模量、密度和正弦激勵(lì)幅值的不確定性Table 1 Uncertainties of elastic modulus,density of the material and the amplitude of the harmonic excitation

懸臂梁整體運(yùn)動(dòng)微分方程為

其中,M是懸臂梁整體質(zhì)量矩陣,與材料密度ρ 有關(guān);K是整體剛度矩陣,與材料彈性模量E有關(guān);F(t)是整體載荷向量,與正弦激勵(lì)P(t)有關(guān).

令y=則整體運(yùn)動(dòng)方程(63)可以表示為

首先考察在時(shí)間步長(zhǎng)?t=40μs 下利用辛算法求解哈密頓系統(tǒng)(64)和直接求解方程(63)方法計(jì)算節(jié)點(diǎn)11 的響應(yīng)標(biāo)稱值.利用辛算法求解哈密頓系統(tǒng)(64)計(jì)算得到的t=0～0.1 s 內(nèi)的響應(yīng)標(biāo)稱值曲線如圖2(a)所示,整體呈周期性變化,但直接求解方程(63)得到的響應(yīng)標(biāo)稱值很短時(shí)間內(nèi)即發(fā)散,如圖2(b)所示,只有當(dāng)時(shí)間步長(zhǎng)足夠小,如?t=2μs 時(shí),直接求解方程(63)才能得到和利用辛算法求解哈密頓系統(tǒng)(64)相同的結(jié)果.這一穩(wěn)定性差異反映了哈密頓系統(tǒng)辛算法能夠保持體系結(jié)構(gòu)特征,體現(xiàn)出利用哈密頓系統(tǒng)辛算法求解微分方程的優(yōu)越性.

圖2 時(shí)間步長(zhǎng)?t=40μs 下利用不同算法計(jì)算得到的節(jié)點(diǎn)11 的響應(yīng)標(biāo)稱值Fig.2 The nominal value of the response of node 11 obtained by different algorithms with time step ?t=40μs

圖2 時(shí)間步長(zhǎng)?t=40μs 下利用不同算法計(jì)算得到的節(jié)點(diǎn)11 的響應(yīng)標(biāo)稱值(續(xù))Fig.2 The nominal value of the response of node 11 obtained by different algorithms with time step ?t=40μs(continued)

圖3 本文所提方法和傳統(tǒng)方法計(jì)算得到的節(jié)點(diǎn)11 的響應(yīng)曲線Fig.3 The response curve of node 11 obtained by SHPM,IHPM,TSM and TIM

在同一時(shí)間步長(zhǎng)?t=2 μs 下本文所提隨機(jī)、區(qū)間非齊次線性哈密頓系統(tǒng)的參數(shù)攝動(dòng)法和傳統(tǒng)隨機(jī)、區(qū)間方法計(jì)算得到的節(jié)點(diǎn)11 在t=0～5.0 ms 內(nèi)的響應(yīng)曲線如圖3 所示,其中圖3(a)～圖3(b)分別為兩種隨機(jī)、區(qū)間方法得到的響應(yīng)曲線,圖3(c)為本文所提隨機(jī)、區(qū)間方法計(jì)算得到的響應(yīng)曲線,響應(yīng)標(biāo)稱值也繪制于圖3 中.本文所提隨機(jī)、區(qū)間方法和傳統(tǒng)隨機(jī)、區(qū)間方法計(jì)算得到的節(jié)點(diǎn)11 在t=3.5 ms 時(shí)刻的位移上下界如表2 所示.

表2 本文所提隨機(jī)、區(qū)間方法和傳統(tǒng)隨機(jī)、區(qū)間方法計(jì)算得到的節(jié)點(diǎn)11 在t=3.5 ms 時(shí)刻的位移上下界Table 2 The upper and lower bounds on the displacement of node 11 at t=3.5 ms obtained by SHPM,IHPM TSM and TIM

由圖3 和表2 可知,本文所提隨機(jī)、區(qū)間方法得到的響應(yīng)上下界曲線分別與傳統(tǒng)隨機(jī)、區(qū)間方法所得響應(yīng)上下界曲線均幾乎完全重合,結(jié)果非常相近,驗(yàn)證了所提隨機(jī)、區(qū)間方法的準(zhǔn)確性和有效性.此外,本文所提區(qū)間方法得到的響應(yīng)區(qū)間范圍包含本文所提隨機(jī)方法得到的響應(yīng)區(qū)間范圍,即區(qū)間方法所得響應(yīng)上界大于隨機(jī)方法所得響應(yīng)上界,而區(qū)間方法下界小于隨機(jī)方法下界,這一現(xiàn)象與前述理論推導(dǎo)相符.

前述內(nèi)容驗(yàn)證了在較大時(shí)間步長(zhǎng)?t=40μs 下哈密頓系統(tǒng)辛算法相較于直接求解運(yùn)動(dòng)微分方程方法計(jì)算響應(yīng)標(biāo)稱值的有效性和優(yōu)越性,也驗(yàn)證了在較小時(shí)間步長(zhǎng)?t=2μs 下本文所提隨機(jī)、區(qū)間方法所得響應(yīng)區(qū)間范圍的準(zhǔn)確性.為了進(jìn)一步檢驗(yàn)本文所提隨機(jī)、區(qū)間方法在較大時(shí)間步長(zhǎng)下所得響應(yīng)區(qū)間范圍也能保持較高精度,下面將本文所提隨機(jī)、區(qū)間方法在時(shí)間步長(zhǎng)?t=40μs 下得到的響應(yīng)區(qū)間范圍與蒙特卡洛模擬方法(簡(jiǎn)記為MCM)得到的響應(yīng)區(qū)間范圍進(jìn)行比較.其中,蒙特卡洛模擬在較小時(shí)間步長(zhǎng)?t=2μs 下進(jìn)行,樣本點(diǎn)數(shù)目設(shè)置為1.0×105,對(duì)于每一樣本點(diǎn)都采用辛算法求解,從而可將蒙特卡洛模擬結(jié)果視為精確值.本文所提隨機(jī)、區(qū)間方法和蒙特卡洛模擬方法計(jì)算得到的節(jié)點(diǎn)11 在t=0～5.0 ms 內(nèi)的響應(yīng)曲線如圖4 所示.

圖4 本文所提隨機(jī)、區(qū)間方法和蒙特卡洛模擬方法計(jì)算得到的節(jié)點(diǎn)11 的響應(yīng)曲線Fig.4 The response curve of node 11 obtained by SHPM,IHPM and MCM

由圖4 所示,本文所提隨機(jī)方法與蒙特卡洛模擬方法得到的響應(yīng)上下界曲線差別微小,盡管由于較大時(shí)間步長(zhǎng)原因?qū)е戮扔幸欢ㄏ陆?但仍得到了較為滿意的結(jié)果.這一現(xiàn)象體現(xiàn)出本文所提方法在較大時(shí)間步長(zhǎng)下也能得到具有較高精度的結(jié)果,再次驗(yàn)證了本文所提方法的優(yōu)勢(shì).而對(duì)于本文所提區(qū)間方法,由于區(qū)間擴(kuò)張效應(yīng),得到的響應(yīng)區(qū)間范圍包含本文所提隨機(jī)方法和蒙特卡洛模擬方法得到的響應(yīng)區(qū)間范圍,與前述理論推導(dǎo)相一致.

6.2 復(fù)合材料層合板

考慮如圖5 所示的邊長(zhǎng)為L(zhǎng)=1.0 的四邊固支的復(fù)合材料層合板結(jié)構(gòu),其由5 層正交各向異性材料鋪設(shè)而成,鋪設(shè)角度為(0?,90?,0?,90?,0?),每層厚度為t=0.012,質(zhì)量密度ρ=1.0.板中心處受一垂直于板的正弦激勵(lì)作用P(t)=?psin(200πτ),初始條件為(0)=0,x(0)=0.材料屬性和正弦激勵(lì)幅值同樣均具有不確定性,假設(shè)其區(qū)間范圍I如表3 所示;同時(shí)假設(shè)它們?cè)趨^(qū)間上服從正態(tài)分布,均值μ和標(biāo)準(zhǔn)差σ也如表3 所示.采用4 結(jié)點(diǎn)矩形薄板單元,將板劃分為相等的16 個(gè)單元.為方便,所有數(shù)據(jù)均采用無量綱量.

圖5 四邊固支復(fù)合材料層合板Fig.5 A composite laminate with fully clamped boundary conditions

表3 材料屬性和正弦激勵(lì)幅值的不確定性Table 3 Uncertainties of material properties and the amplitude of the harmonic excitation

考慮板中心處z方向的動(dòng)力響應(yīng).在同一時(shí)間步長(zhǎng)?τ=1.0×10?5下,本文所提隨機(jī)、區(qū)間方法和傳統(tǒng)隨機(jī)、區(qū)間方法在τ=0～0.05 內(nèi)計(jì)算得到的響應(yīng)范圍以及響應(yīng)標(biāo)稱值曲線如圖6 所示,其中圖6(a)～圖6(b)分別為兩種隨機(jī)、區(qū)間方法得到的響應(yīng)曲線,圖6(c)為本文所提隨機(jī)、區(qū)間方法計(jì)算得到的響應(yīng)曲線.本文所提隨機(jī)、區(qū)間方法和傳統(tǒng)隨機(jī)、區(qū)間方法計(jì)算得到的τ=0.025 時(shí)刻位移上下界如表4 所示.

由圖6 和表4 可知,本文所提隨機(jī)、區(qū)間方法得到的響應(yīng)上下界曲線分別與傳統(tǒng)隨機(jī)、區(qū)間方法所得響應(yīng)上下界曲線同樣近乎完全重合,再次驗(yàn)證了所提方法的有效性.同時(shí),本文所提區(qū)間方法得到的響應(yīng)上界比隨機(jī)方法得到的響應(yīng)上界大,區(qū)間方法得到的響應(yīng)下界比隨機(jī)方法得到的響應(yīng)下界小,也與理論推導(dǎo)和算例6.1 現(xiàn)象一致.

圖6 本文所提方法和傳統(tǒng)方法計(jì)算得到的板中心處z 方向的響應(yīng)曲線Fig.6 The response curve of the center of the plate in z-direction obtained by SHPM,IHPM,TSM and TIM

下面進(jìn)一步比較在時(shí)間步長(zhǎng)?τ=2.5×10?4下利用哈密頓系統(tǒng)辛算法和直接求解運(yùn)動(dòng)微分方程計(jì)算得到的響應(yīng)標(biāo)稱值,曲線分別如圖7(a)～圖7(b)所示.由圖7 可知,利用哈密頓系統(tǒng)辛算法求解能夠得到滿意的響應(yīng)標(biāo)稱值,然而在該時(shí)間步長(zhǎng)下直接求解運(yùn)動(dòng)微分方程得到的響應(yīng)標(biāo)稱值同樣在很短時(shí)間內(nèi)發(fā)散.這一現(xiàn)象再次體現(xiàn)了哈密頓系統(tǒng)辛算法由于能夠保持體系結(jié)構(gòu)特征而具有很好的穩(wěn)定性,凸顯了哈密頓系統(tǒng)辛算法在數(shù)值模擬中的顯著優(yōu)越性.

表4 本文所提隨機(jī)、區(qū)間方法和傳統(tǒng)隨機(jī)、區(qū)間方法計(jì)算得到的板中心處z 方向在τ=0.025 時(shí)刻的位移上下界Table 4 The upper and lower bounds on the displacement of the center of the plate in z-direction at τ=0.025 obtained by SHPM,IHPM TSM and TIM

圖7 時(shí)間步長(zhǎng)?τ=2.5×10?4 下利用不同算法計(jì)算得到的板中心處z 方向的響應(yīng)標(biāo)稱值Fig.7 The nominal value of the response of the center of the plate in z-direction obtained by different algorithms with time step ?τ=2.5×10?4

7 結(jié)論

本文考慮非齊次線性哈密頓系統(tǒng),基于攝動(dòng)理論發(fā)展了含小參數(shù)擾動(dòng)的哈密頓系統(tǒng)的參數(shù)攝動(dòng)法,以此為基礎(chǔ)分別將小擾動(dòng)看作隨機(jī)、區(qū)間擾動(dòng),提出了分別針對(duì)隨機(jī)、區(qū)間不確定性哈密頓系統(tǒng)的參數(shù)攝動(dòng)法,突破了傳統(tǒng)哈密頓系統(tǒng)的限制,并由切比雪夫不等式證明了兩者響應(yīng)分析結(jié)果的相容性關(guān)系,通過兩個(gè)綜合考慮外載荷和不確定性的結(jié)構(gòu)動(dòng)力響應(yīng)數(shù)值算例驗(yàn)證得到了如下結(jié)論:

(1)在同一較大時(shí)間步長(zhǎng)下,由直接求解運(yùn)動(dòng)微分方程的方法得到的響應(yīng)標(biāo)稱值很可能會(huì)發(fā)散,而將運(yùn)動(dòng)微分方程引入到哈密頓系統(tǒng)后再利用辛算法計(jì)算仍然可以得到令人滿意的結(jié)果,體現(xiàn)了哈密頓系統(tǒng)辛算法具有很好的穩(wěn)定性,保持體系結(jié)構(gòu)特征的優(yōu)勢(shì)明顯,在數(shù)值仿真模擬方面的顯著優(yōu)越性;

(2)在同一較小時(shí)間步長(zhǎng)下,本文所提隨機(jī)不確定性哈密頓系統(tǒng)的參數(shù)攝動(dòng)法得到的響應(yīng)上下界分別與傳統(tǒng)隨機(jī)方法所得響應(yīng)上下界結(jié)果近乎完全一致,差別微小,本文所提區(qū)間方法得到的響應(yīng)上下界也分別與傳統(tǒng)區(qū)間方法所得響應(yīng)上下界極為相近,驗(yàn)證了所提方法的有效性;

(3)本文所提隨機(jī)方法在較大時(shí)間步長(zhǎng)下得到的響應(yīng)上下界與蒙特卡洛模擬方法在較小時(shí)間步長(zhǎng)下得到的響應(yīng)上下界基本吻合,說明本文所提方法在較大時(shí)間步長(zhǎng)下也能得到具有較高精度的結(jié)果,再次驗(yàn)證了本文所提方法的優(yōu)越性;

(4)當(dāng)參數(shù)的區(qū)間范圍由概率統(tǒng)計(jì)信息確定,即區(qū)間半徑表示為距其均值距離為標(biāo)準(zhǔn)差的整數(shù)倍時(shí),本文所提區(qū)間方法計(jì)算得到的響應(yīng)區(qū)間范圍包含本文所提隨機(jī)方法計(jì)算得到的響應(yīng)區(qū)間范圍,即區(qū)間響應(yīng)上界大于隨機(jī)響應(yīng)上界、區(qū)間響應(yīng)下界小于隨機(jī)響應(yīng)下界.