李 爽，翟長(zhǎng)海，劉洪波 ，3，謝禮立 ，2

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

Newmark積分方法在負(fù)剛度時(shí)的數(shù)值穩(wěn)定性分析

李 爽1，翟長(zhǎng)海1，劉洪波1，3，謝禮立1，2

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

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

Newmark積分方法;負(fù)剛度;穩(wěn)定性;阻尼;分析步長(zhǎng)

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

本文研究了Newmark積分方法在負(fù)剛度時(shí)的數(shù)值穩(wěn)定性分析方法并對(duì)其穩(wěn)定性進(jìn)行了分析.同時(shí)，詳細(xì)分析了Newmark積分方法中的參數(shù)、阻尼和分析步長(zhǎng)等對(duì)積分方法穩(wěn)定性的影響規(guī)律.

1 負(fù)剛度時(shí)穩(wěn)定性分析方法

結(jié)構(gòu)的荷載－位移曲線可以采用線性化的方法處理為分段線性形式，最不利的情況下也可以假定在每一計(jì)算時(shí)間步長(zhǎng)內(nèi)或迭代步長(zhǎng)內(nèi)結(jié)構(gòu)為一線性體系.因此，對(duì)于方法穩(wěn)定性分析可轉(zhuǎn)化為給定剛度情況下積分方法的穩(wěn)定性分析.使用單自由度線性體系討論積分方法的穩(wěn)定性，動(dòng)力平衡方程為

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

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

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

對(duì)于使用Newmark積分方法求解式(1)的情況，當(dāng)式(1)中r=1時(shí)A矩陣的標(biāo)準(zhǔn)形式已有研究給出［1］．當(dāng) r≠1 時(shí)，可將式(1)寫成如下形式

通過特征多項(xiàng)式|A－λI|= 0，得到A的3個(gè)特征值分別為

積分逼近算子A是3×3型矩陣，且有3個(gè)互不相同的特征值，根據(jù)矩陣?yán)碚?，其所?duì)應(yīng)的初等因子必皆為一次.穩(wěn)定性的判別根據(jù)文獻(xiàn)［10］中給出的負(fù)剛度時(shí)方法穩(wěn)定性的判斷準(zhǔn)則進(jìn)行判斷，滿足下式即認(rèn)為方法穩(wěn)定

2 負(fù)剛度時(shí)Newmark積分方法的穩(wěn)定性分析

2.1 Newmark方法中參數(shù)的影響

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

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

2.2 阻尼的影響

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

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

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

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

2.3 分析時(shí)間步長(zhǎng)的影響

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

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

3 已有研究的說明

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

圖5 參數(shù)的不同取值對(duì)穩(wěn)定性的影響

4 方法的合理使用

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

5 結(jié)語

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

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

［2］ 王勖成.有限單元法［M］.北京:清華大學(xué)出版社，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］ 李小軍.地震工程中動(dòng)力方程求解的逐步積分方法［J］.工程力學(xué)， 1996，13(2):110－118.

［6］ 劉祥慶，劉晶波，丁樺.時(shí)域逐步積分算法穩(wěn)定性與精度的對(duì)比分析［J］.巖石力學(xué)與工程學(xué)報(bào)， 2007，26(S1):3000－3008.

［7］ 程民憲.負(fù)剛度條件下數(shù)值積分的收斂性和穩(wěn)定性［J］.力學(xué)學(xué)報(bào)， 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法在負(fù)剛度條件下的收斂性和穩(wěn)定性［J］.重慶建筑大學(xué)學(xué)報(bào)， 2001，23(1):21－24.

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

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

LI Shuang1，ZHAI Chang－h(huán)ai1，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.

國(guó)家自然科學(xué)基金資助項(xiàng)目( 90705021，51008101).

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

翟長(zhǎng)海(1976—)，男，教授，博士生導(dǎo)師;

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

(編輯 趙麗瑩)