Newmark積分方法在負剛度時的數值穩定性分析

李 爽，翟長海，劉洪波 ，3，謝禮立 ，2

(1.哈爾濱工業大學 土木工程學院，150090 哈爾濱，shuangli@hit.edu.cn;2.中國地震局工程力學研究所，150080 哈爾濱;3.黑龍江大學建筑工程學院，150080 哈爾濱)

Newmark積分方法在負剛度時的數值穩定性分析

李 爽1，翟長海1，劉洪波1，3，謝禮立1，2

(1.哈爾濱工業大學 土木工程學院，150090 哈爾濱，shuangli@hit.edu.cn;2.中國地震局工程力學研究所，150080 哈爾濱;3.黑龍江大學建筑工程學院，150080 哈爾濱)

為了獲得Newmark積分方法在負剛度時的適用性，研究了Newmark積分方法在負剛度時的數值穩定性分析方法并對其穩定性進行了分析，分析了Newmark積分方法中的參數、阻尼和分析步長對積分方法穩定性的影響規律.結果表明其穩定性與正剛度時差異較大，在負剛度時不再具有無條件穩定性.對于正剛度時常用的Newmark積分參數α=1/4、δ=1/ 2，在負剛度時僅當分析步長大于某一值時積分方法才能保證穩定性，因此難以同時滿足穩定性和準確性的要求.推薦了通過選取合適的Newmark積分參數獲得穩定性的方法.為Newmark積分方法用于分析具有負剛度的結構時的適用性提供了依據，也為積分方法的合理選擇提供了參考.

Newmark積分方法;負剛度;穩定性;阻尼;分析步長

在動力反應分析過程中，結構剛度可能出現負定的情況.產生這一現象的原因可源于一些材料在加載至荷載極限點后變形增加荷載反而減少，進而引起結構荷載－位移曲線出現負剛度段;也可源于幾何非線性效應，附加的幾何剛度引起結構荷載－位移反應出現負剛度段.在時域逐步積分方法中，Newmark積分方法是求解結構動力反應的常用方法.關于Newmak積分方法在正剛度時的數值穩定性問題已有大量研究［1－6］，一致的結論是該積分方法在積分參數α≥(0.5+δ)2/4且δ≥1/2時是無條件穩定的.Newmark積分方法的無條件穩定性是在結構具有正剛度時獲得的，負剛度時是否成立需進一步研究.然而，目前對于結構剛度矩陣為負定時積分方法穩定性的探討僅限于少量研究.文獻［7］指出Newmark常平均加速度方法 ( α=1/4、δ=1/2)滿足屈服后負剛度條件下數值積分的穩定性要求.文獻［8］雖沒有直接研究負剛度情況方法的穩定性，但指出在負剛度情況下Newmark方法中參數取值為α=1/4、δ=1/2時精度較低，而取 α=1/12、δ=1/2時得到的結果較好.文獻［9］通過研究得到的結果為負剛度情況下取α≥(0.5+δ)2/4且δ≥1/2時是條件穩定的.從目前的研究水平來看，關于負剛度情況Newmak方法穩定性方面的研究并不全面，而且也沒有形成一致的認識和結論.

本文研究了Newmark積分方法在負剛度時的數值穩定性分析方法并對其穩定性進行了分析.同時，詳細分析了Newmark積分方法中的參數、阻尼和分析步長等對積分方法穩定性的影響規律.

1 負剛度時穩定性分析方法

結構的荷載－位移曲線可以采用線性化的方法處理為分段線性形式，最不利的情況下也可以假定在每一計算時間步長內或迭代步長內結構為一線性體系.因此，對于方法穩定性分析可轉化為給定剛度情況下積分方法的穩定性分析.使用單自由度線性體系討論積分方法的穩定性，動力平衡方程為

式中:t+Δt¨u、t+Δt˙u和t+Δtu分別為t+Δt時刻結構加速度、速度和位移;t+Δt¨ug為t+Δt時刻的地震動加速度;ξ為阻尼比;ω為自振頻率;r為體系的剛度系數，彈性階段r= 1，屈服后r=k1/k0，負剛度時r＜0.k0為初始剛度，k1為第二剛度.

可以將式(1)表示為以下遞推形式

矩陣A和向量L分別為積分逼近算子和荷載算子.

對于使用Newmark積分方法求解式(1)的情況，當式(1)中r=1時A矩陣的標準形式已有研究給出［1］．當 r≠1 時，可將式(1)寫成如下形式

通過特征多項式|A－λI|= 0，得到A的3個特征值分別為

積分逼近算子A是3×3型矩陣，且有3個互不相同的特征值，根據矩陣理論，其所對應的初等因子必皆為一次.穩定性的判別根據文獻［10］中給出的負剛度時方法穩定性的判斷準則進行判斷，滿足下式即認為方法穩定

2 負剛度時Newmark積分方法的穩定性分析

2.1 Newmark方法中參數的影響

在正剛度情況下，一般選用分析步長Δt/T=1/10～1/ 20，T為結構周期，所以本文也給出在此種情況下負剛度時的情況.圖1給出了Newmark參數 δ=1/2、阻尼比 ξ=0.05 時 ρ(A)/eμΔt值隨Newmark參數α和負剛度系數r的變化規律.當Δt/T=1/10時，存在一個峰值帶，在此區域內的ρ(A)/eμΔt值遠大于 1，即處于此區域內 Newmark積分方法不穩定;在其他區域內，ρ(A)/eμΔt值接近1或小于 1，圖1(a)和(c)中非網格區域 ＜1.當剛度系數r接近于0時，α的變化對ρ(A)/eμΔt值的變化影響較小.隨著剛度系數r的減小，α的變化對ρ(A)/eμΔt值變化的影響為先增大然后部分區域會出現峰值帶，而后減小;當剛度系數r較小時，較小或較大的α值對應的方法均是穩定的，而中等大小的α值對應的方法卻是不穩定的.當Δt/T=1/20時，沒有出現像Δt/T=1/10時所出現的峰值帶，隨剛度系數r減小或隨α的增大ρ(A)/eμΔt值呈增大趨勢. 對 Δt/T=1/10 和Δt/T=1/20的兩種情況，α取小值時均是穩定的，并且此時剛度系數r的影響很小.

圖2為Newmark參數α=1/4、阻尼比ξ=0.05時ρ(A)/eμΔt值隨 Newmark參數 δ和剛度系數 r的變化規律.當Δt/T=1/10時且剛度系數r接近于0時，δ的變化對ρ(A)/eμΔt值的變化影響較小;ρ(A)/eμΔt值有隨剛度系數r減小而增大的趨勢，隨δ的增大而增大的趨勢.圖2(a)和(c)中非網格區域＜ 1，因此僅當r較小時才能保證方法的穩定性.Δt/T=1/20時，無論剛度系數r和δ取何值均能保證方法的穩定.ρ(A)/eμΔt值在剛度系數r接近0時接近 1，有隨剛度系數r減小而減小的趨勢，有隨δ增大而增大的趨勢.

2.2 阻尼的影響

圖3為Newmark參數α =1/4、δ=1/2時ρ(A)/eμΔt值隨阻尼比 ξ變化的情況.可以看出如果分析時間步取小一點，阻尼的影響也相應的小.但是，對于 α=1/4、δ=1/2的情況，ρ(A)/eμΔt均是 ＞ 1 的，無論阻尼比取何值方法均是不穩定的.

圖 1 δ=1/2、ξ=0.05 時 ρ(A)/eμΔt隨 α的變化

圖2 α =1/4、ξ=0.05 時 ρ(A)/eμΔt隨 δ的變化

圖 3 α =1/4、δ=1/2 時 ρ(A)/eμΔt隨 ξ的變化

2.3 分析時間步長的影響

圖4為 Newmark參數 α =1/4、δ=1/2、ξ=0.05 時 ρ(A)/eμΔt值隨分析步長 Δt/T 以及剛度系數r的變化情況.圖4(a)中非網格區域＜ 1，因此有趣的是Δt/T較小時方法卻不能保證穩定性，而只有在Δt/T大于某一值時，方法才是穩定的.如果僅從方法穩定性來講，取大一點的分析步是可以的.但是從計算精度來講，Δt/T≤0.5時(即一個周期內有2個分析步)才有可能得到準確的結果.同時，對于多自由度體系的分析，Δt/T大于某一值才穩定的方法顯然不能同時滿足穩定性和準確性的要求，因為對低階模態滿足了穩定性的要求就很難滿足高階模態的準確性要求.此外，穩定區域與不穩定區域之間存在著峰值帶.分析方法的穩定性和步長有關，所以并不像正剛度時的無條件穩定，負剛度時Newmark積分方法取α=1/4、δ=1/2是條件穩定的.

圖4 α =1/4、δ=1/2時ρ(A)/eμΔt隨Δt/T 的變化

3 已有研究的說明

文獻［8］通過算例分析指出在負剛度時Newmark方法中參數取值為α=1/4、δ=1/2時精度較低，而參數取值為α=1/12、δ=1/2時得到的結果較好.圖5為使用本文方法，針對文獻［8］中算例得到的穩定性分析結果，可見對于該算例情況，在 α =1/12、δ=1/2、ξ=0.1、r= －1時積分方法恰好是穩定的，即 ρ(A)/eμΔt值 ＜ 1，而在α =1/4、δ=1/2、ξ=0.1、r=－1時積分方法是不穩定的.圖5中同時還給出了參數變化對積分方法穩定性的影響，可見只有在特定的一些情況積分方法才能獲得穩定.

圖5 參數的不同取值對穩定性的影響

4 方法的合理使用

在正剛度時無條件穩定的Newmark積分方法(α≥(0.5+δ)2/4、δ≥1/2)在負剛度時不再是無條件穩定的.正剛度時被認為精度最好的α=1/4、δ=1/2情況只有在Δt/T大于某一值時才能保證方法穩定.因為難以同時滿足穩定性和準確性的要求，所以不適合做具有負剛度的多自由度體系的非線性動力分析.由圖1可知，當α取值較小時，Newmark積分方法的穩定性可以保證，因此在實際分析中α可取一小值.一個極端的情況是使用Newmark顯式積分方法，即α=0.

5 結語

Newmark積分方法在正剛度時被認為精度最好的α=1/4、δ=1/2情況在負剛度時只有在分析步長大于某一值時才能保證方法穩定性，因此是條件穩定的，同時對于多自由度體系難以同時滿足穩定性和準確性的要求.給出了負剛度時Newmark直接積分方法中的參數、阻尼和分析步長等對積分方法穩定性的影響規律.

［1］ BATHE K J.Finite element procedures［M］.New Jersey:Prentice－Hall， 1996，806 － 810.

［2］ 王勖成.有限單元法［M］.北京:清華大學出版社，2003:492－494.

［3］ ZIENKIEWICZ O C，TAYLOR R L，ZHU J Z.The finite element method，Volumn 1:the basis［M］.6th ed.Oxford:Butterworth－Heinemann，2005:606 －611.

［4］ CHOPRA A K.Dynamics of structures:theory and applications to earthquake engineering［M］.2nd ed.Beijing:Tsinghua University Press，2005:174－178.

［5］ 李小軍.地震工程中動力方程求解的逐步積分方法［J］.工程力學， 1996，13(2):110－118.

［6］ 劉祥慶，劉晶波，丁樺.時域逐步積分算法穩定性與精度的對比分析［J］.巖石力學與工程學報， 2007，26(S1):3000－3008.

［7］ 程民憲.負剛度條件下數值積分的收斂性和穩定性［J］.力學學報， 1988，20(1):31－40.

［8］ XIE Y M.On the accuracy of time－stepping schemes for dynamic problems with negative stiffness［J］.Communications in Numerical Methods in Engineering， 1993，9(2):131－137.

［9］ 吳云芳.Newmark法在負剛度條件下的收斂性和穩定性［J］.重慶建筑大學學報， 2001，23(1):21－24.

［10］ 程民憲.結構動力分析中幾種逐步積分法［J］.工程力學， 1989，6(2):35－47.

On the stability of Newmark integration method for structure with negative stiffness

LI Shuang1，ZHAI Chang－hai1，LIU Hong－bo1，3，XIE Li－li1，2

(1.School of Civil Engineering，Harbin Institute of Technology，150090 Harbin，China，shuangli@hit.edu.cn;2.Institute of Engineering Mechanics，China Earthquake Administration，150080 Harbin，China;3.School of Civil Engineering and Architecture，Heilongjiang University，150080 Harbin，China)

To study the applicability of Newmark integration method for structure with negative stiffness，the method is first investigated，and then the numerical stability of the integration method is analyzed.The influences of Newmark integration parameters，damping and time increment on numerical stability of Newmark integration method are also discussed.The result shows that the integration method，which is very different from the case of positive stiffness，is not unconditional numerical stable in the case of negative stiffness.It is found that the usually used Newmark integration method in the case of positive stiffness with α =1/ 4，δ=1/2 is difficult to achieve a balance between stability and accuracy because it is stable only when the time increment lager than a critical value in the case of negative stiffness.The strategy to obtain numerical stability is proposed by selecting appropriate Newmark integration parameters.This paper presents references for applicability of Newmark integration method and selection of rational integration method when structures with negative stiffness.

newmark integration method;negative stiffness;stability;damping;time increment

O241

A

0367－6234(2011)08－0001－05

2010－01－14.

國家自然科學基金資助項目( 90705021，51008101).

李 爽(1981—)，男，博士，講師;

翟長海(1976—)，男，教授，博士生導師;

謝禮立(1939—)，男，教授，中國工程院院士.

(編輯 趙麗瑩)