基于統(tǒng)計線性化的非線性隨機振動虛擬激勵法研究及應(yīng)用

羅曉峰，項貽強

(1.浙江大學(xué) 土木工程系，杭州 310058;2.浙江工業(yè)職業(yè)技術(shù)學(xué)院，紹興 312000)

我國是一個多發(fā)地震國家，確保橋梁結(jié)構(gòu)在地震荷載下的安全運營具有十分重要的意義，以往的經(jīng)驗表明進行科學(xué)合理的抗震設(shè)計是減輕因地震造成損失的有效途徑，目前國內(nèi)外橋梁結(jié)構(gòu)的抗震設(shè)計中，反應(yīng)譜法、時程法以及隨機振動法等傳統(tǒng)的計算方法已經(jīng)得到廣泛應(yīng)用［1-3］。由于地震地面運動實際上是一種隨機激勵，因此結(jié)構(gòu)的抗震設(shè)計以結(jié)構(gòu)的隨機地震響應(yīng)為基礎(chǔ)更為科學(xué)、合理。隨機振動法以功率譜密度作為分析的核心，所以又稱功率譜法，這種方法可以較充分地考慮地震發(fā)生的統(tǒng)計概率特征，并且適合于跨度較大、結(jié)構(gòu)較復(fù)雜的多點支撐的橋梁，被認為是一種較為先進合理的分析工具，1995年，歐洲結(jié)構(gòu)抗震設(shè)計規(guī)范(Eurocode 8)把功率譜方法作為結(jié)構(gòu)抗震設(shè)計計算的工具［4-5］，我國現(xiàn)行的公路橋梁抗震設(shè)計規(guī)范也將功率譜法作為結(jié)構(gòu)抗震設(shè)計計算的方法之一。Penzien等在80年代初曾將核電站管線簡化后，分別從時域和頻域的角度研究了其在不均勻場地激勵下的安全問題，得出忽略參振振型之間的相關(guān)性和場地不均勻性都會導(dǎo)致很大的誤差;Lin等將輸油管線視為有多個支點的連續(xù)梁，用隨機振動法研究了其抗震性能;Perotti、Harichandran等也研究了多點隨機振動方法。上述方法由于計算效率低，且只能處理較簡單的結(jié)構(gòu)，難以在實際的橋梁工程中得到廣泛地應(yīng)用。近年來，林家浩教授［1］從計算力學(xué)的角度提出了簡單易行，計算效率較高，并在理論上屬精確算法的虛擬激勵法，張笈瑋等［6］在傳統(tǒng)虛擬激勵法的基礎(chǔ)上，提出了一種計算非平穩(wěn)地震激勵下結(jié)構(gòu)峰值響應(yīng)均值的簡便算法;江洋等［7］提出了一種只需少量振型即可獲得高精度的虛擬激勵算法—修正的絕對位移法。但是到目前為止，虛擬激勵法只是成功地解決了結(jié)構(gòu)線性隨機地震響應(yīng)計算的工作量問題，對于計入非線性的橋梁結(jié)構(gòu)的抗震分析，國內(nèi)外還還沒有較深入研究，因此探討與發(fā)展非線性體系隨機地震響應(yīng)的虛擬激勵法顯得十分迫切與重要［8］。

本文在統(tǒng)計等效線性化的基礎(chǔ)上，應(yīng)用FAP數(shù)值解法［9-11］，提出了一種簡便、高效、實用的，可以考慮結(jié)構(gòu)非線性的改進的虛擬激勵法，這種方法避開了分解位移和求解擬靜模態(tài)矩陣的復(fù)雜過程，可以在一般的通用有限元軟件中應(yīng)用。作為算例，以跨徑75m+55m的單塔雙索面預(yù)應(yīng)力混凝土斜拉橋為背景進行抗震分析，應(yīng)用不同的抗震分析方法計算其位移、速度、加速度響應(yīng)功率譜密度，及其位移響應(yīng)歷程和軸力響應(yīng)值，說明基于統(tǒng)計線性化的絕對位移直接求解的虛擬激勵法的精確性。

1 統(tǒng)計線性化方法

1954年美國學(xué)者Booten［12］提出了一種可以處理非線性力學(xué)系統(tǒng)隨機響應(yīng)問題的方法—統(tǒng)計等效線性化方法，對于n個自由度的非線性力學(xué)系統(tǒng)可以用如下微分方程式表達，即:

式中［M］是一個n階的結(jié)構(gòu)質(zhì)量矩陣;［G({y}，{})］為非線性恢復(fù)力和阻尼力之和的等效矩陣。

建立與非線性系統(tǒng)等效的線性化方程:

式中［Ce］為等效粘阻系數(shù)矩陣，［Ke］等效剛度矩陣。

將式(1)與式(2)作差，于是得到:

式中{y}={y(cij，kij，t)}，為方程(2)的穩(wěn)態(tài)解，其中cij，kij為系數(shù)矩陣［Ce］，［Ke］的元素。

誤差矩陣［ε］是cij，kij的函數(shù)，統(tǒng)計線性化的實質(zhì)就是要求合理地選擇cij，kij，使

式中E是數(shù)學(xué)期望。

通過數(shù)學(xué)計算，可以求得式(4)的解，進而求得非線性系統(tǒng)的等效線性剛度［Ke］和等效粘滯阻尼［Ce］。

2 絕對位移直接求解的虛擬激勵法

由文獻［1］中傳統(tǒng)虛擬激勵法的求解原理，可知在求解虛擬動態(tài)相對位移和虛擬靜位移時，都需要計算并提取擬靜模態(tài)矩陣A=-·，十分麻煩，而且有些軟件沒有矩陣提取功能，所以一般不能在通用有限元軟件中求解，需編制相應(yīng)的專用求解程序，要解決這一問題，可參考文獻［7］。

為此，本文在文獻［13］研究的基礎(chǔ)上，提出了基于統(tǒng)計線性化的絕對位移直接求解的虛擬激勵法，具體實施如下:

將方程(2)寫成分塊矩陣的形式［1，13］:

式中yb代表N個支座的地面強迫位移，ys代表非支座節(jié)點的位移，F(xiàn)b代表地面作用于N個支座的力，下標(biāo)s，b代表結(jié)構(gòu)內(nèi)部節(jié)點自由度和支座節(jié)點自由度。

構(gòu)造虛擬加速度激勵:

將與Mb的乘積作為Fb代入式(5)，可得:

將上式按第二行展開，得:

將上式兩邊同時乘以M-1b，得:

顯然，當(dāng)Mb為一大質(zhì)量，即Mb→∞時為零，則由式(8)可得支座虛擬加速度激勵為:

由上述可知，只需對支座賦予一大質(zhì)量塊，并對大質(zhì)量塊施加虛擬力激勵，即可實現(xiàn)虛擬加速度激勵的加載，則支座虛擬速度激勵和虛擬位移激勵為:

將式(6)按第一行展開，可得:

將式(10)代入式(11)，可以求得虛擬絕對位移?ys，按照虛擬激勵法，絕對位移的功率譜矩陣為:

關(guān)于傳統(tǒng)的隨機振動求解原理，以及動態(tài)相對位移自功率譜密度、擬靜態(tài)相對位移自功率譜密度及兩者之間的互功率譜密度的求解，很多學(xué)者對其進行了深入的研究，這里不再贅述。

3 算例分析及比較

3.1 工程背景及模型信息

取一典型的橋長130 m、跨徑布置為75+55 m的兩跨單塔雙索面預(yù)應(yīng)力混凝土斜拉橋為例，其塔墩梁為固結(jié)的結(jié)構(gòu)體系。主梁截面為雙實心邊主梁大懸臂截面，主梁中心高1.9 m，頂板寬38 m，懸臂長4.5 m，主梁側(cè)實心梁寬3 m，背跨側(cè)實心梁寬4 m，實心梁間頂板厚0.28 m。背跨部分梁段由于配重的需要而增設(shè)底板形成箱型截面。主梁采用雙向預(yù)應(yīng)力體系，主塔為鋼筋混凝土斜塔，塔中心線與水平夾角為75°，橋面以上垂直高為50.7 m，主塔采用變截面實心矩形，順橋向截面高度從3 m(塔頂部)變化到8 m(橋上塔根部):橫橋向?qū)挾葹?.5 m，用通用有限元軟件ANSYS進行建模，主梁、主塔采用C50混凝土，斜拉索采用高強鋼絲PSEM7-241，典型主梁截面型式用有限元模擬如圖1所示，計算模型對塔用空間梁單元beam4進行離散，主梁用帶剛臂的空間梁單元beam4簡化為魚骨刺形，斜拉索采用link10單元，左橋端給予豎向平移自由度約束，右橋端給予橫向的平移自由度約束，索塔底部完全約束，索單元和梁單元采用鉸約束，如圖2所示。

圖1 典型的主梁截面型式(單位:m)Fig.1 Typical girder section type(unit:m)

圖2 斜拉橋模型Fig.2 Mode of cable-stayed bridge

3.2 計算概要

為了驗證絕對位移直接求解的虛擬激勵法的正確性，本文應(yīng)用統(tǒng)計線性化方法將體系線性化，之后采用有限元軟件中的隨機振動模塊和諧響應(yīng)模塊進行分析，線性化過程及問題的求解可參考文獻［9］，文獻中充分利用計算機的特點并結(jié)合FAP數(shù)值解法較好地解決了多自由度非線性系統(tǒng)在隨機激勵下的響應(yīng)問題，F(xiàn)AP是統(tǒng)計線性化方法的數(shù)值解法，能方便的用各種算法語言編制出相應(yīng)的計算機程序代碼，本文應(yīng)用FORTRAN 90編制了FAP算法的程序，以實現(xiàn)結(jié)構(gòu)的線性化，限于篇幅，不再贅述。

對于本文提出的絕對位移直接求解的虛擬激勵法，可以用簡諧響應(yīng)分析中的Full法，隔振支座采用combin40，剛度為k=6×106N/m，模態(tài)分析時取前150階振型，振型參與質(zhì)量達到97%，以保證算法的足夠精度，在采用絕對絕對位移直接求解的虛擬激勵法時，在模型基底的4個支坐節(jié)點上附加mass21大質(zhì)量單元，大質(zhì)量單元與支座節(jié)點剛接，質(zhì)量取主橋質(zhì)量的106倍，然后釋放x，y方向自由度的約束，按上述理論提出的方法，對每個大質(zhì)量塊施加虛擬力然后進行諧響應(yīng)分析，求得計入非線性后的虛擬激勵響應(yīng)，進而按照式(12)求得絕對位移的功率譜密度矩陣。計算時輸入的功率譜采用由文獻［7］提出的GB50011-2001修正反應(yīng)譜模型迭代得到的實用地震動功率譜［7，14-15］。一致地震動參數(shù):7度設(shè)防，Ⅲ類場地，第一組設(shè)計地震分組，視波速取100m/s，相干函數(shù)采用L-W模型，整個結(jié)構(gòu)采用Rayleigh阻尼，取基底四個節(jié)點進行三向一致激勵(基地的四個節(jié)點分別為:左橋端約束節(jié)點、右橋端約束節(jié)點以及索塔底部橫橋向中心線兩側(cè)的約束節(jié)點)。

3.3 絕對位移法與隨機振動法的結(jié)果對比分析

圖3和圖4分別給出了該斜拉橋主跨和背跨跨中的響應(yīng)功率譜密度圖，包括位移、速度和加速度的響應(yīng)功率譜密度。

圖3 主跨跨中響應(yīng)功率譜密度對比圖Fig.3 Contrast of response power spectral density in the middle of main span

圖4 背跨跨中響應(yīng)功率譜密度對比圖Fig.4 Contrast of response power spectral density in the middle of back span

從計算結(jié)果可以看出，將體系線性化之后，兩種算法的計算結(jié)果基本一致，并且通過研究發(fā)現(xiàn)，其它主梁和主塔節(jié)點的響應(yīng)功率譜密度也有類似規(guī)律，說明絕對位移法的正確性。再者，通過觀察對比，發(fā)現(xiàn)絕對位移法計算得到位移、速度和加速度響應(yīng)功率譜密度值較之于隨機振動法計算得到的值均偏小，從曲線與橫坐標(biāo)圍成的面積來分析，絕對位移法得到的地震隨機振動波的均方值偏小，即施加給斜拉橋的能量偏小，說明隨機振動法是一種偏于保守的抗震分析方法。

由于輸入的加速度功率譜密度的頻率范圍為0.48-15.9 Hz，因此圖中沒有體現(xiàn)頻率在0-0.48 Hz之間的曲線段，實質(zhì)上，在頻率為零的附近，輸入的激勵會導(dǎo)致結(jié)構(gòu)的響應(yīng)功率譜密度急劇增大，一般的功率譜密度模型會對地震動低頻分量進行合理縮減，在f=0時，功率譜密度一般為零［2］。對于位移和速度響應(yīng)功率譜，在f≥3 Hz之后，其值很小，因此在圖中沒有體現(xiàn)，對于加速度反應(yīng)譜，在f≥3 Hz后，其值有一定范圍的波動，變化曲線如圖3、圖4的(c)圖所示。

絕對位移直接求解的虛擬激勵法與傳統(tǒng)的隨機振動法相比，其計算結(jié)果的相對誤差如表1所示，由表中數(shù)據(jù)可以得出如下結(jié)論:

(1)兩種算法的計算結(jié)果基本接近，說明線性化后，絕對位移直接求解的虛擬激勵法在一定程度上可以代替?zhèn)鹘y(tǒng)隨機振動法來研究結(jié)構(gòu)的抗震性能，此方法的優(yōu)越性具有如下幾點:① 計算時不需要將絕對位移分解成擬靜力位移和動態(tài)相對位移兩部分;②構(gòu)造虛擬激勵進行簡諧響應(yīng)分析時不需要計算擬靜模態(tài)矩陣A，簡單易行;③ 其實質(zhì)并沒有改變傳統(tǒng)虛擬激勵法的基本原理，其計算效率也和傳統(tǒng)虛擬激勵法相當(dāng)。

(2)絕對位移法相比于隨機振動法，加速度響應(yīng)功率譜密度的誤差最小，46號節(jié)點最大誤差僅為1.15%，154號節(jié)點最大誤差僅為1.37%，造成這一現(xiàn)象的主要原因是，在模型的前期分析時，輸入的激勵為加速度功率譜密度，所以在經(jīng)過一系列的計算后，其對應(yīng)的加速度響應(yīng)功率譜密度是最接近實際情況的，而位移、速度響應(yīng)功率譜密度經(jīng)過了比加速度響應(yīng)譜密度更多的計算步驟，使其誤差累積，所以46號節(jié)點位移功率譜密度最大誤差達到11.1%，154號節(jié)點位移功率譜密度最大誤差達到9.14%。

表1 兩種算法的誤差對比表Tab.1 Contrast of two algorithms'error

3.4 絕對位移法、時程法及譜分析的對比分析

目前抗振分析中，除上述提出的隨機振動法和諧響應(yīng)分析法之外，譜分析和動力時程法也是很重要的研究方法，動力時程分析法其實質(zhì)為瞬態(tài)分析的Full法，可以考慮各種非線形因素，因此對于本文研究的算例，可以跳過線性化過程直接對其進行動力時程分析，計算時考慮隔陣支座，結(jié)構(gòu)為非正交阻尼體系，輸入的人工波按上述絕對位移法所采用的條件生成，地震作用采用雙向水平一致地震激勵，按照抗震規(guī)范，次方向的地震加速度乘以0.85的系數(shù)。

為了將上述時程分析的非線性結(jié)果與線性結(jié)果做對比分析，此處采用了文獻［7］提出的GB50011-2001修正反應(yīng)譜模型進行譜分析，仍采用Rayleigh阻尼的形式，其它條件同上［16］。

圖5為主跨跨中節(jié)點和背跨跨中節(jié)點的位移時程圖，通過對比，可以發(fā)現(xiàn)如下規(guī)律:

(1)由于譜分析是在線彈性范圍內(nèi)研究斜拉橋的位移響應(yīng)，因此較之于動力時程法及線性化處理后的絕對位移法，其位移響應(yīng)值均偏小，相比于位移響應(yīng)最大的動力時程法，最大偏差達到50%以上，說明在研究結(jié)構(gòu)的抗震性能時，必須考慮結(jié)構(gòu)的非線性，否則會引起較大誤差，使計算結(jié)果偏離真實值。

(2)計入非線性后的絕對位移法與動力時程法的計算結(jié)果非常接近，對于46號節(jié)點X向位移最大誤差僅為7.8%，對于154號節(jié)點Y向最大誤差為7.3%，并且通過研究發(fā)現(xiàn)，其它主梁、主塔及斜拉索節(jié)點的位移響應(yīng)也有類似規(guī)律，說明線性化后的絕對位移法可以充分考慮結(jié)構(gòu)的非線性，可以作為抗震設(shè)計的一種行之有效的方法。

圖5 節(jié)點位移對比圖Fig.5 Contrast of node displacement

圖6 軸力響應(yīng)對比圖Fig.6 Contrast of axial force response

圖6為三種計算方法在一致地震動作用下典型單元的Z向軸力響應(yīng)均值對比圖，從圖中可以看出，線性化后的絕對位移法求得的單元軸力響應(yīng)值最大，并且與動力時程法計算出的單元軸力值基本相同，最大誤差僅為0.23%(主跨跨中附近的單元)，對于譜分析，其Z向軸力響應(yīng)值均偏小，較之于動力時程法，最大誤差達到8%，說明非線性對結(jié)構(gòu)抗震性能有一定得影響，在橋梁結(jié)構(gòu)的抗震分析中必須加以考慮。

4 結(jié)論

(1)相比于傳統(tǒng)的虛擬激勵法，基于統(tǒng)計線性化的絕對位移直接求解的虛擬激勵法，首先充分考慮了結(jié)構(gòu)的非線性，使其抗震分析結(jié)果更加接近于真實情況;其次在計算時，不需要將絕對位移分解成擬靜力位移和動態(tài)相對位移，在構(gòu)造虛擬激勵進行簡諧響應(yīng)分析時，不需要計算擬靜模態(tài)矩陣，可以在一般的通用有限元分析軟件中直接應(yīng)用，方便虛擬激勵法在實際工程中的應(yīng)用。

(2)體系線性化后，絕對位移法和隨機振動法計算得到的節(jié)點響應(yīng)功率譜密度基本相同，說明絕對位移法的正確性，并且計算時輸出的響應(yīng)功率譜密度的精確性依賴于所輸入的激勵的類型。

(3)相比于動力時程法及線性化處理后的絕對位移法，譜分析沒有考慮結(jié)構(gòu)的非線性，其節(jié)點位移響應(yīng)值和單元軸力響應(yīng)值均偏小，說明非線性對結(jié)構(gòu)的抗震性能有較大影響，在橋梁結(jié)構(gòu)的抗震分析中必須加以考慮。

(4)線性化后的絕對位移法與動力時程法計算出的節(jié)點位移響應(yīng)值和單元軸力響應(yīng)值基本一致，說明基于統(tǒng)計線性化的絕對位移直接求解的虛擬激勵法與動力時程法一樣，可以有效地分析和研究非線性體系的抗震性能，唯一要加以區(qū)別的是，采用絕對位移法進行諧響應(yīng)分析時是以頻率為間隔輸入所需的激勵，而動力時程法是以時間為間隔建立地震波數(shù)據(jù)，但其實質(zhì)是相同的。

(5)基于統(tǒng)計線性化的絕對位移直接求解的虛擬激勵法突破了隨機振動計算效率低的瓶頸，并且比常規(guī)的反應(yīng)譜法更精確;其計算時間僅為動力時程法的0.38倍(以本文所提算例為研究對象)，而且動力時程法的計算結(jié)果過分依賴所選取的加速度時間歷程曲線，要使用許多條地震記錄曲線分別計算后再做統(tǒng)計平均，因此本文所提方法比動力時程法更高效，可廣泛應(yīng)用于跨度較大，結(jié)構(gòu)較復(fù)雜的橋梁體系的抗震分析中。

［1］林家浩，張亞輝.隨機振動的虛擬激勵法［M］.北京:科學(xué)出版社，2004.

［2］Ernesto H Z，Vanmarcke E H.Seismic random-vibration analysis of multisupport-structural systems［J］.Journal of Engineering Mechanics，1994，120(5):1107-1128.

［3］Abdullateef M A，Ashraf A.An efficient multi-point supportmotion random vibration analysis technique［J］.Computers ＆Structures，2002，80(8):1689-1697.

［4］Lin J H，Zhang Y H，Li Q S，et al.Seismic spatial effects for long-span bridges，using the pseudo excitation method ［J］.Engineering Structures，2004，26(9):1207-1216.

［5］Leper P，Ide I M，Paultre P.Multiple-support seismic analysis of large structures［J］.Computers＆ Structures，1990，36(6):1153-1158.

［6］張笈瑋.考慮地震空間效應(yīng)的大跨度空間結(jié)構(gòu)抗震分析與設(shè)計方法［D］.天津:天津大學(xué)，2008.

［7］江 洋.大跨建筑結(jié)構(gòu)多點輸入地震響應(yīng)計算與抗震設(shè)計方法研究［D］.北京:清華大學(xué)，2010.

［8］汪夢甫.非比例阻尼線性體系平穩(wěn)隨機地震響應(yīng)計算的虛擬激勵法［J］.計算力學(xué)學(xué)報，2008，25(1):94-99.

WANG Meng-pu.Pseudo-excitation method for stationary random responses of non-proportionally damped MDOF systems［J］.Chinese Journal of Computational Mechanics，2008，25(1):94-99.

［9］唐國梁，姚昌仁.統(tǒng)計線性化方法及其計算機數(shù)值解［J］.北京工業(yè)學(xué)院學(xué)報，1986，1(1):142-152.

TANG Guo-liang，YAO Chang-ren.Statistical linearization technique and its computer numerical solution［J］.Journal of Beijing University of Technology，1986，1(1):142-152.

［10］林建龍，盧翰庭.統(tǒng)計線性化方法在大型貨車垂向振動分析中的應(yīng)用［J］.鐵道學(xué)報，1990，12(3):1-7.

LIN Jian-long，LU Han-ting.Research on application of the statistical linearization technique to analyzing the vertical vibration response of freight cars［J］.Journal of the China Railway Society，1990，12(3):1-7.

［11］潘家強，多自由度非線性振動系統(tǒng)的統(tǒng)計線性化簡化方程［J］.浙江大學(xué)學(xué)報，1983，1(3):90-99.

PAN Jia-qiang，Simplified statistical linearization equations of discrete nonlinear systems［J］.Journal of Zhejiang University，1983，1(3):90-99.

［12］ Iwan W D，Yang I M.Application of statistical linearization technique to nonlinear multidegree-of-freedom systems［J］.Transaction of the ASME-Journal of Applied Mechanics，1972，39(2):545-550.

［13］李永華，李思明.絕對位移直接求解的虛擬激勵法［J］.振動與沖擊，2009，28(10):185-190.

LI Yong-hua，LI Si-ming.Pseudo excitation method based on solving absolute displacement［J］.Journal of Vibration and Shock，2009，28(10):185-190.

［14］GB50011-2001.建筑抗震設(shè)計規(guī)范［S］.北京:中國建筑工業(yè)出版社，2001.

GB50011-2001.Code for Seismic Design of Buildings［S］.Beijing:China Architecture ＆ Building Press，2001.

［15］石永久，江 洋，王元清.直接求解法在結(jié)構(gòu)多點輸入地震響應(yīng)計算中的應(yīng)用與改進［J］.工程力學(xué)，2011，28(1):75-81.

SHI Yong-jiu，WANG Yang，WANG Yuan-qing.Application and improvement of direct solving method in seismic response analysis of structures under multi-support excitations［J］.Journal of Engineering Mechanics，2011，28(1):75-81.

［16］ Kiureghian A D，Neuenhofer A.Response spectrum method for multi-support seismic excitations［J］.Earthquake Engineering＆ Structural Dynamics，1992，21(8):713-740.