付相球，潘旦光，2

（1.北京科技大學土木工程系，北京100083；2.北京科技大學城市地下空間工程北京市重點實驗室，北京100083）

引言

非比例阻尼體系廣泛存在于現(xiàn)有的結(jié)構(gòu)體系中，如下部鋼筋混凝土上部鋼結(jié)構(gòu)的豎向混合結(jié)構(gòu)［1-2］、附加阻尼器的振動控制結(jié)構(gòu)［3］，土-結(jié)構(gòu)相互作用體系［4-7］等。雖然直接積分法可用于非比例阻尼體系的動力反應分析，但計算工作量較大。模態(tài)疊加法計算效率高，且可識別對結(jié)構(gòu)動力反應具有顯著貢獻的模態(tài)，因而得到廣泛的應用［8-9］。

對于非比例阻尼體系，F(xiàn)oss［10］提出以狀態(tài)空間法求解非比例阻尼系統(tǒng)的復模態(tài)，該方法可以得到精確的動力反應。但計算維數(shù)擴大了一倍，且將實數(shù)域運算轉(zhuǎn)換到復數(shù)域，計算量約為實模態(tài)分析的8倍［11］。為提高非比例阻尼體系模態(tài)疊加法的計算效率，Cronin［12］通過忽略阻尼矩陣中的非對角元素，強制解耦非比例阻尼體系以應用實模態(tài)疊加法求解體系的近似反應，簡單易行但誤差難以控制［13-14］。Roesset等［15］提出等效阻尼比的概念用于近似解耦非比例阻尼，并常被用于混合結(jié)構(gòu)的地震反應分析［16-17］，這種方法實際上也屬于強制解耦方法，各階模態(tài)的阻尼比在一定程度上考慮了耦合阻尼的影響，提高了計算精度，但不能解決計算誤差無法估計的問題［18］。Ibrahimbegovic等［11］提出了應用實模態(tài)進行非比例阻尼體系動力反應的迭代模態(tài)疊加法，這種方法將耦合阻尼的影響作為非線性荷載，得到高精度的計算結(jié)果。張靜等［19］針對大規(guī)模非對稱實矩陣的標準特征值問題，發(fā)展了適用于大型復特征值問題求解的Lanczos-QR方法，提高了計算效率，同時避免了范數(shù)誤差。事實上，非比例阻尼體系復模態(tài)疊加法的理論是完備的，問題在于狀態(tài)空間復模態(tài)求解及復數(shù)域的疊加計算時間較長。為提高復模態(tài)的計算效率，可將非比例阻尼系統(tǒng)看成由比例阻尼體系攝動后得到的新系統(tǒng)，應用攝動法求解特征 方程［20-21］。Cha［22］提出一階攝動方 法求解弱非比例阻尼系統(tǒng)；Tang和Wang［23］提出一種可以處理重根問題的攝動方法。但一般的攝動法只適用于弱非比例阻尼的體系，且需要用到原系統(tǒng)完整的模態(tài)空間。Lou等［24-25］提出的模態(tài)攝動法利用原系統(tǒng)的非完整模態(tài)即可得到二階攝動精度，并已應用到無阻尼離散和連續(xù)系統(tǒng)特征值問題的分析。Pan等［26］應用模態(tài)疊加法基本原理，構(gòu)建等效比例阻尼系統(tǒng)作為原系統(tǒng)，提出一種可求解強非比例體系特征值問題的計算方法，并以一個兩自由度體系和一個單層框架結(jié)構(gòu)驗證計算方法的精度。

模態(tài)攝動法本質(zhì)上屬于Ritz法，計算精度與附加模態(tài)的數(shù)目有關(guān)，文獻［26］中兩個算例的自由度都比較小，采用完整模態(tài)或接近完整模態(tài)進行計算，結(jié)果的精度都很高。本文在文獻［26］的基礎(chǔ)上，進一步研究自由度較多體系非完整模態(tài)下模態(tài)攝動法的計算精度，在此基礎(chǔ)上，建立一種基于無阻尼體系實模態(tài)的復模態(tài)疊加法。同時，以非比例阻尼體系的地震反應為例，討論復模態(tài)疊加法中的模態(tài)截斷問題，分析了基于累積振型參與質(zhì)量、累積振型貢獻系數(shù)及累積振型加速度貢獻系數(shù)所得模態(tài)數(shù)對計算精度的影響，并對比研究本文方法和狀態(tài)空間復模態(tài)疊加法計算效率的差異。最后，將本文方法應用到一個帶附加阻尼鋼結(jié)構(gòu)框架地震反應分析，驗證了本文方法的適用性。

1 基于實模態(tài)的復模態(tài)疊加法

一個N自由度線性黏滯阻尼體系的強迫振動方程為

式中u，u?和u?分別為N×1階的位移、速度和加速度向量；f(t)為荷載向量；m，c和k分別為N×N階的質(zhì)量、阻尼和剛度矩陣。對于非比例阻尼體系，狀態(tài)空間法把N維2階微分方程轉(zhuǎn)變?yōu)?N維1階微分方程［10］

為應用模態(tài)疊加法求解式（2）的運動方程，須求解特征方程

式中γ和ψ分別為復特征值和相應的特征向量。若已知前r（r≤N）對特征值γj和γj+r=γˉj，以及相應的特征向量ψj和ψj+r=ψˉj（j=1，2，…，r），其中“”表示復數(shù)共軛。則式（2）強迫振動的反應y可表示為

式中zj為廣義坐標，zˉj表示zj的共軛復數(shù)，r為截斷的模態(tài)數(shù)。由模態(tài)疊加法的基本原理，廣義坐標zj（j=1，2，…，r）的運動方程為

從計算過程看，基于狀態(tài)空間的復模態(tài)疊加法與常規(guī)模態(tài)疊加法相同，但是，式（3）特征值求解的維數(shù)增加了1倍且為復數(shù)運算，使計算時間顯著增加。為提高計算效率，關(guān)鍵是減少式（3）和式（4）的計算時間。為此，本文基于模態(tài)攝動法原理，以無阻尼體系的實模態(tài)為基礎(chǔ)，計算非比例阻尼的復模態(tài)，以減少計算時間。

非比例阻尼矩陣可分解為

式中cp為相應的比例阻尼矩陣，Δc為比例阻尼矩陣與原阻尼矩陣的偏差。若已知無阻尼體系特征方程前n階自振頻率ωj和相應的模態(tài)φj（j=1，2，…，n），則比例阻尼矩陣cp滿足正交特性

將m，cp和k組成的系統(tǒng)稱之為等效比例阻尼系統(tǒng)，作為非比例阻尼體系的原系統(tǒng)。則采用狀態(tài)空間法表示的等效比例阻尼體系中的特征值方程為

非比例阻尼體系新系統(tǒng)的特征值γj和模態(tài)ψj可以利用原系統(tǒng)的特征值sj和模態(tài)ηj進行簡單的攝動分析而近似地求得［26］，即設(shè)：

式 中D11，D12，D21，D22，E11，E12，E21和E22都 為n×n階 的 方 陣；R1和R2為n×1階向 量；x為2n×1階向量。各矩陣和向量的元素具體表達見附錄。

式（12）將復特征方程（3）轉(zhuǎn)換成了2n維的非線性代數(shù)方程組，求解代數(shù)方程通常比求解特征方程要簡單得多，本文采用Newton-Raphson迭代方法求解。在得到未知向量x后，非比例阻尼體系的第j個特征值為

采用有阻尼體系的頻率、阻尼比和特征值的表示方法

則：

式中Re代表取元素的實部。為了避免與相應的等效比例阻尼的自振頻率ωj和阻尼比ζj混淆，和表示第j階的擬無阻尼自振頻率和阻尼比。當式（10）中的n=N時，2N個線性無關(guān)的特征向量ηj組成2N維的完備空間，因此，非比例阻尼體系的特征向量Ψj可以在這個2N維的空間中精確地展開，此時攝動解也將收斂到狀態(tài)空間下的精確解。而通常高階模態(tài)相比于低階模態(tài)影響很小，因此方程（12）的維數(shù)n通常比N小得多，這表明本文所采用的攝動法無需求解原系統(tǒng)所有階模態(tài)即可得到高精度的攝動解，從而減少計算時間。

式（12）中的j從1到r依次取值可得非比例阻尼體系前r階復特征值和復模態(tài)。由此，式（5）中pj(t)和aj用實模態(tài)表示為：

方程（5）得到廣義坐標zj的解。令其中quj與qlj分別為qj的前n和后n個元素組成的向量。由式（4）和（10）可得體系的位移為

式中qu和ql分別為r個quj與qlj向量組成的n×r維的矩陣，z={z1z2...zr}T為r×1的廣義坐標。由式（16），式（17）和式（18）可知，基于實模態(tài)的復模態(tài)疊加法將復模態(tài)、廣義單自由度體系的部分系數(shù)和位移表示為實模態(tài)的線性組合，簡化了計算。

注重農(nóng)業(yè)技術(shù)創(chuàng)新，加快完善和填補優(yōu)勢農(nóng)畜產(chǎn)品和特色農(nóng)畜產(chǎn)品標準化生產(chǎn)操作規(guī)程，加強標準推廣和使用指導，大力宣傳培訓農(nóng)牧業(yè)投入品使用規(guī)范，督促生產(chǎn)經(jīng)營主體嚴格落實間隔期休藥期規(guī)定。發(fā)展農(nóng)牧業(yè)適度規(guī)模經(jīng)營，推進農(nóng)畜產(chǎn)品標準化生產(chǎn)基地建設(shè)，提高農(nóng)畜產(chǎn)品生產(chǎn)經(jīng)營主體的專業(yè)化、組織化程度，先行對農(nóng)畜產(chǎn)品新型經(jīng)營主體開展農(nóng)畜產(chǎn)品質(zhì)量安全綜合指數(shù)評價，將是否按標生產(chǎn)作為政策支持的重要條件，從源頭上增強農(nóng)畜產(chǎn)品經(jīng)營主體的安全意識，推進農(nóng)畜產(chǎn)品的安全生產(chǎn)。

2 算例分析

為研究地震作用下基于實模態(tài)的復模態(tài)疊加法適用性，借鑒文獻［11］的算例構(gòu)建圖1所示的框架結(jié)構(gòu)，討論模態(tài)攝動法的精度、模態(tài)截斷指標以及復模態(tài)疊加法的計算精度和效率。該框架有限元模型包括27個梁單元共72個自由度，材料的彈性模量和密度分別為40 N/m2，1 kg/m3，梁截面慣性矩和面積分別為1 m4，1 m2，單元長度為1 m。采用集中質(zhì)量模型，由靜力凝聚消除轉(zhuǎn)動自由度后，體系包含48個有效無阻尼實模態(tài)?？蚣芙Y(jié)構(gòu)在節(jié)點3，8，12，17，20與25的水平方向有附加阻尼器，阻尼器系數(shù)分別為c1=3 N·s/m，c2=5 N·s/m。以耦合系數(shù)α表示非比例阻尼矩陣耦合的程度［27］

圖1 附加阻尼器的三層框架結(jié)構(gòu)Fig.1 Three-layer frame with concentrated dampers

式中Clk表示阻尼矩陣C中的第l行第k列元素。體系的耦合系數(shù)為1，為強非比例阻尼體系。輸入地震時程為1940年5月在加州地震中記錄到的El Centro波的N-S方向的加速度分量，如圖2所示。

圖2 El Centro波加速度時程Fig.2 El Centro ground acceleration history

2.1 模態(tài)攝動法的精度

當模態(tài)攝動法采用完整的無阻尼體系特征向量空間計算時，其結(jié)果收斂到精確解，但大型復雜結(jié)構(gòu)通常僅計算有限階模態(tài)，為保證計算精度，計算非比例阻尼體系第j階的模態(tài)時，式（10）中所需的無阻尼體系的模態(tài)數(shù)n為式中Δn為計算所需的附加模態(tài)數(shù)。以狀態(tài)空間法所得的擬無阻尼頻率與阻尼比為精確解，n階模態(tài)下由模態(tài)攝動法得到的和為近似解，則擬無阻尼頻率與阻尼比的相對誤差為：

圖3所示為模態(tài)攝動法所得前3階模態(tài)相對誤差隨附加模態(tài)數(shù)的變化情況。表1列出前3階采用不同方法所得擬無阻尼頻率與阻尼比的計算結(jié)果，表2為不同方法的計算誤差。其中n=48表示完整模態(tài)空間，強制解耦法為直接去除模態(tài)阻尼矩陣的非對角元素的解耦方法。

表1 框架結(jié)構(gòu)的擬無阻尼自振頻率與阻尼比Tab.1 The quasi-undamped natural frequencies and damping ratios of the frame system

表2 框架結(jié)構(gòu)的擬無阻尼自振頻率與阻尼比的相對誤差/%Tab.2 The errors of quasi-undamped natural frequencies and damping ratios of the frame system/%

圖3 前3階頻率與阻尼比的相對誤差Fig.3 Relative errors of the first three natural frequencies and damping ratios

計算結(jié)果表明：（1）對于模態(tài)攝動法，隨著附加模態(tài)數(shù)的增加，前3階頻率和阻尼比趨近于精確解；采用完整實模態(tài)的模態(tài)攝動法計算結(jié)果與精確解一致。附加模態(tài)數(shù)為8時，前3階頻率的相對誤差小于1%，前3階模態(tài)阻尼比的相對誤差小于2%。這表明對于大型復雜結(jié)構(gòu)，當采用附加模態(tài)數(shù)大于8的非完整模態(tài)進行模態(tài)攝動法計算時，即可得到很高的計算精度。（2）強制解耦方法計算所得的頻率誤差小于5%，但阻尼比誤差可達9.596%。這是因為強制解耦過程改變了阻尼矩陣，直接對模態(tài)阻尼比產(chǎn)生影響，因此阻尼比與精確解存在明顯的差別；而由單自由度體系有阻尼頻率可知，當阻尼比小于0.2時，阻尼引起體系有阻尼頻率的變化很小，因此強制解耦方法計算時頻率與精確解相差不大。

2.2 復模態(tài)疊加法地震反應分析

對于一般結(jié)構(gòu)而言，少數(shù)低階模態(tài)的疊加即可得到滿足工程要求的解，因此存在模態(tài)截斷的問題［28］。然而對于非比例阻尼體系的復模態(tài)疊加法，能否采用實模態(tài)疊加法中的模態(tài)截斷方法需要進一步的研究。針對本文的三層框架結(jié)構(gòu)，三類實模態(tài)的模態(tài)截斷指標：累積振型參與質(zhì)量、首層平動位移的累積振型貢獻系數(shù)以及首層平動的累積振型加速度貢獻系數(shù)如表3所示。以模態(tài)截斷指標超過90%作為模態(tài)選取的依據(jù)，則由累積振型參與質(zhì)量、首層平動位移的累積振型貢獻系數(shù)和首層平動的累積振型加速度貢獻系數(shù)選取的模態(tài)數(shù)分別為5階，3階和20階。

表3 不同模態(tài)截斷指標計算結(jié)果/%Tab.3 Results of different modal truncation indexes/%

式中ei-max表示第i個自由度位移與精確解的峰值誤差，ei-sum表示第i個自由度位移與精確解的累積誤差表示第i個自由度位移的精確解，ui(t)表示第i個自由度位移的近似解，T為積分時長。在El Centro波作用下，前3階、前5階、前20階和全部48階復模態(tài)采用式（18）計算所得頂層、中間層、底層水平位移時程如圖4所示，在采用模態(tài)攝動法計算特征值時，附加實模態(tài)的個數(shù)均取8。同時為了進行比較，將狀態(tài)空間法精確解也一并繪出。各自由度的峰值誤差、累積誤差如表4所示。

由圖4與表4可以看出：1）當采用48階復模態(tài)疊加計算時，計算結(jié)果與狀態(tài)空間法完全重合，說明完整模態(tài)的復模態(tài)疊加法可以得到精確解；2）底層水平位移的模態(tài)截斷誤差最大，頂層水平位移的模態(tài)截斷誤差最小，說明底層水平位移反應對高階模態(tài)更為敏感；3）當采用前3階、前5階復模態(tài)疊加計算時，各層水平位移的誤差較大，其中底層水平位移峰值誤差達到30%，累積誤差超過了50%，說明對于該框架體系，采用累積振型參與質(zhì)量和首層平動位移的累積振型貢獻系數(shù)作為模態(tài)截斷指標難以滿足計算精度要求；4）當采用前20階復模態(tài)疊加計算時，頂層和中間層的水平位移峰值誤差及累積誤差均在10%以內(nèi)，底層水平位移累積誤差偏大，略大于10%，但峰值誤差只有5.8%，結(jié)合圖4（c）的位移時程結(jié)果可以看出前20階復模態(tài)疊加時底層水平位移與精確解結(jié)果較吻合?？傮w而言，累積振型參與質(zhì)量和首層平動位移的累積振型貢獻系數(shù)所取的復模態(tài)數(shù)偏少，基于首層平動的累積振型加速度貢獻系數(shù)的模態(tài)截斷方法可以滿足計算精度要求。

表4 三層框架結(jié)構(gòu)水平位移誤差/%Tab.4 Errors of horizontal displacements of the threelayer frame/%

圖4 地震作用下水平位移計算結(jié)果Fig.4 Horizontal displacement under earthquake excitation

為對比不同方法的計算精度，圖5給出了前20階模態(tài)下本文方法與狀態(tài)空間法和強制解耦非比例阻尼矩陣后的實模態(tài)疊加法（簡稱實模態(tài)疊加法）的水平位移計算結(jié)果，并與精確解相比較。不同方法相應的峰值誤差及累積誤差如表5所示。

由圖5、表5結(jié)果可以看出：1）本文方法與狀態(tài)空間法的前20階計算結(jié)果接近，兩種方法的時程曲線與精確解吻合很好；三層水平位移的峰值誤差均小于10%，除底層位移累積誤差偏大高于10%以外，頂層及中間層位移累積誤差也均小于10%，顯示了良好的精度；2）實模態(tài)疊加法計算結(jié)果與精確解存在明顯誤差，三層水平位移的累積誤差均大于20%，其中底層位移的峰值誤差也高于20%，這說明對于強非比例阻尼的體系，實模態(tài)疊加方法忽略了非對角線阻尼的影響，導致較大的計算誤差。

表5 三層框架結(jié)構(gòu)水平位移誤差/%Tab.5 Errors of horizontal displacements of the threelayer frame/%

圖5 地震作用下三種方法的水平位移計算結(jié)果Fig.5 Horizontal displacement of three methods under earthquake excitation

為對比本文方法與狀態(tài)空間法的計算效率，前20階模態(tài)疊加下計算時間如表6所示，計算過程統(tǒng)一采用2.19 GHz處理器的筆記本，Matlab計算軟件。可以看出，模態(tài)攝動法求解特征方程所耗費的時間小于狀態(tài)空間法直接求解復特征方程的時間，說明模態(tài)攝動法提高了求解特征方程的效率；本文方法由于計算過程中部分采用實模態(tài)線性組合，減少了復數(shù)運算，計算效率也明顯高于狀態(tài)空間方法。綜合計算精度結(jié)果，本文方法兼顧了計算效率與精度，適用于強非比例阻尼體系的地震反應分析。

表6 計算時間/msTab.6 Calculation time/ms

3 五層鋼結(jié)構(gòu)地震反應分析

為驗證本文方法對真實結(jié)構(gòu)的適用性，下面采用不同方法分析北京市通州區(qū)某公司總部辦公樓的地震反應。該辦公樓主體為鋼框架結(jié)構(gòu)，層高17.7 m，建筑長50.4 m，分為7跨，其中一榀橫向框架鋼結(jié)構(gòu)計算簡圖如圖6所示。鋼結(jié)構(gòu)的梁和柱采用Q345鋼，樓板選用C25混凝土。橫向7.2 m，6 m跨度框架梁截面為HN450 mm×200 mm×9 mm×14 mm，橫向3 m跨度框架梁截面為HN300 mm×150 mm×6.5 mm×9 mm（數(shù)字代表工字鋼的截面高度，截面寬度，腹板厚度，翼緣厚度）；底層柱為400 mm×400 mm×16 mm×16 mm箱型截面，其他層柱為400 mm×400 mm×14 mm×14 mm（數(shù)字代表箱型柱（方管）截面高度，截面寬度，上下厚度，左右厚度）箱型截面。辦公樓重力荷載如樓板、墻體等以等效質(zhì)量的方式施加在框架結(jié)構(gòu)上，其中屋面重力荷載取6.38 kN/m2，樓面重力荷載取7 kN/m2。鋼結(jié)構(gòu)阻尼比取0.02，采用瑞利阻尼分析，瑞利阻尼系數(shù)的參考頻率選用基頻與荷載卓越頻率。同時為提高該框架結(jié)構(gòu)的抗震性能，各樓層橫向A-B跨和CD跨之間各設(shè)置一個黏滯阻尼器，采用斜對角布置如圖6所示，附加阻尼器的阻尼系數(shù)為c=1500 kN·s/m。

圖6 辦公樓計算簡圖Fig.6 Calculation diagram of an office building

該框架結(jié)構(gòu)的基本自振周期為1 s，耦合系數(shù)α=0.57。累積振型加速度貢獻系數(shù)超過90%時，所需的最少模態(tài)數(shù)為7階。在El Centro波作用下，該框架結(jié)構(gòu)的底層和頂層位移反應時程如圖7所示，峰值與累積誤差如表7所示。其中模態(tài)攝動法計算特征值時，附加實模態(tài)的個數(shù)仍取8。

圖7 El Centro波作用下三種方法的水平位移計算結(jié)果Fig.7 Horizontal displacement of three methods under El Centro wave

由圖7及表7可以看出，本文方法所得底層與頂層位移的峰值誤差均在1%以內(nèi)，累積誤差不超過5%，具有很高的計算精度。而實模態(tài)疊加法底層位移峰值誤差大于10%，累積誤差達29.40%?？傮w而言，本文方法計算精度高，可用于強非比例阻尼體系的地震反應分析，而強制解耦的實模態(tài)疊加法易造成無法預計的誤差。

表7 結(jié)構(gòu)水平位移誤差/%Tab.7 Errors of horizontal displacements of the frame/%

4 結(jié)論

工程實際結(jié)構(gòu)中的阻尼常具有非比例阻尼的特征，本文基于模態(tài)攝動法提出基于實模態(tài)的復模態(tài)疊加法，用于地震荷載下的強迫振動分析。通過理論分析和數(shù)值計算可以得出以下結(jié)論：

1）當采用無阻尼體系所有模態(tài)構(gòu)成的完備空間進行模態(tài)攝動法計算時，可以得到非比例阻尼體系精確的復特征值和特征向量；對于非完備空間，當附加模態(tài)數(shù)大于8時，模態(tài)攝動法所得的復特征值誤差小于2%；

2）非比例阻尼體系的復模態(tài)疊加法計算時，基于累積振型參與質(zhì)量及累積振型貢獻系數(shù)進行模態(tài)截斷所得的模態(tài)數(shù)偏少，由此導致結(jié)構(gòu)地震反應的誤差較大，建議采用累積振型加速度貢獻系數(shù)作為復模態(tài)疊加法的模態(tài)截斷指標；

3）本文方法將非比例阻尼體系的特征向量和強迫振動的解采用實模態(tài)線性組合進行計算，從而減少了復數(shù)運算，提高了計算效率；

4）對于強非比例阻尼體系的地震反應分析，本文方法兼顧了計算效率與精度，而強制解耦非比例阻尼矩陣的實模態(tài)疊加法忽略了耦合阻尼的影響，計算誤差難以控制。

附錄：