考慮阻尼影響的修正CD-Newmark組合算法數(shù)值特性研究1

蔡新江王曉峰田石柱

1）蘇州科技學(xué)院土木工程學(xué)院，蘇州 215011

2）蘇州科技學(xué)院江蘇省結(jié)構(gòu)工程重點(diǎn)實(shí)驗(yàn)室，蘇州215011

3）哈爾濱工業(yè)大學(xué)土木工程學(xué)院，哈爾濱 150090

考慮阻尼影響的修正CD-Newmark組合算法數(shù)值特性研究1

蔡新江1，2）王曉峰2）田石柱1，2，3）

1）蘇州科技學(xué)院土木工程學(xué)院，蘇州 215011

2）蘇州科技學(xué)院江蘇省結(jié)構(gòu)工程重點(diǎn)實(shí)驗(yàn)室，蘇州215011

3）哈爾濱工業(yè)大學(xué)土木工程學(xué)院，哈爾濱 150090

顯-隱式組合數(shù)值積分算法結(jié)合了顯式算法無(wú)需迭代和隱式算法無(wú)條件穩(wěn)定的各自優(yōu)點(diǎn)，是結(jié)構(gòu)抗震擬動(dòng)力試驗(yàn)順利運(yùn)行的關(guān)鍵。在對(duì)傳統(tǒng)顯式中央差分法和隱式Newmark β組合算法進(jìn)行參數(shù)修正的基礎(chǔ)上，建立了修正CD-Newmark算法，考慮阻尼的影響分析了組合算法的穩(wěn)定性條件、周期失真率和數(shù)值阻尼比，分別得到了試驗(yàn)子結(jié)構(gòu)的穩(wěn)定性條件和計(jì)算子結(jié)構(gòu)無(wú)條件穩(wěn)定的參數(shù)合理取值范圍，并對(duì)計(jì)算精度進(jìn)行了分析。通過算例分析驗(yàn)證了算法的數(shù)值特性，從而初步解決了CD-Newmark算法存在穩(wěn)定性界限過嚴(yán)的問題，為結(jié)構(gòu)抗震擬動(dòng)力混合試驗(yàn)提供了研究參考。

數(shù)值積分算法 組合算法 中央差分法 紐馬克法 穩(wěn)定性 精度

引言

子結(jié)構(gòu)擬動(dòng)力混合試驗(yàn)是研究大型復(fù)雜結(jié)構(gòu)抗震性能的一種先進(jìn)實(shí)驗(yàn)方法，目前日益得到研究者的重視。目前其主要發(fā)展方向?yàn)椤皩?shí)時(shí)擬動(dòng)力試驗(yàn)”和“遠(yuǎn)程協(xié)同擬動(dòng)力試驗(yàn)”。前者對(duì)試驗(yàn)系統(tǒng)的要求很高，當(dāng)前的研究主要集中于小比例試件，而針對(duì)大比例試件試驗(yàn)還有待深入研究（徐國(guó)山，2010；王強(qiáng)等，2010）；后者是在常規(guī)擬動(dòng)力試驗(yàn)的基礎(chǔ)上增加了遠(yuǎn)程協(xié)同的功能，即多個(gè)實(shí)驗(yàn)室的試驗(yàn)設(shè)備基于互聯(lián)網(wǎng)連接起來(lái)協(xié)作完成整體結(jié)構(gòu)試驗(yàn)（范云蕾等，2011；王大鵬等，2010）。此外，還有很多學(xué)者對(duì)常規(guī)擬動(dòng)力試驗(yàn)及振動(dòng)臺(tái)與擬動(dòng)力結(jié)合試驗(yàn)方面進(jìn)行了研究（侯杰等，2006；程紹革等，2008）。數(shù)值積分算法是擬動(dòng)力混合試驗(yàn)的核心部分，分為顯式算法和隱式算法兩類。顯式算法的優(yōu)點(diǎn)是不需要反復(fù)迭代，而缺點(diǎn)是算法為條件穩(wěn)定，其穩(wěn)定性與時(shí)間步長(zhǎng)tΔ以及結(jié)構(gòu)的自振頻率ω有關(guān)。主要包括中央差分法和顯式Newmark-β法等。隱式算法的優(yōu)點(diǎn)是無(wú)條件穩(wěn)定，時(shí)間步長(zhǎng)選擇不受試件特征限制，有利于減小及控制誤差，具有較好的數(shù)值穩(wěn)定性和能量耗散性，而缺點(diǎn)是算法需要反復(fù)迭代求解，因此影響了計(jì)算效率，而且在每個(gè)時(shí)間步長(zhǎng)中導(dǎo)入了卸載滯后效應(yīng)，使得實(shí)現(xiàn)試驗(yàn)較為復(fù)雜并對(duì)硬件要求較高。隱式算法包括隱式Newmark-β法、α方法、Houbolt法和Wilson-θ法等。針對(duì)上述算法的優(yōu)缺點(diǎn)，本文介紹一種考慮阻尼影響的修正CD-Newmark組合算法。

1 考慮阻尼影響的修正CD-Newmark組合算法

1.1 算法概述

日本學(xué)者M(jìn). Nakashima等（1985）將顯式中央差分法和隱式Newmark法進(jìn)行組合，提出了CD-Newmark算法并應(yīng)用到子結(jié)構(gòu)擬動(dòng)力混合試驗(yàn)中；由于該方法對(duì)頻率較大的計(jì)算子結(jié)構(gòu)的穩(wěn)定性要求不高，又將隱式Newmark法中的隱式位移表達(dá)式拆分成顯式和隱式兩部分，分別作為試驗(yàn)子結(jié)構(gòu)和計(jì)算子結(jié)構(gòu)，因此稱為PC-Newmark算法，該方法雖然較優(yōu)于CD-Newmark算法，但對(duì)結(jié)構(gòu)類型要求較高；為解決該問題M. Nakashima等（1985）又進(jìn)一步提出了OS算法，將恢復(fù)力分為線性部分和非線性部分，線性部分采用顯式Newmark法，非線性部分采用預(yù)測(cè)-校正方法；為了減小預(yù)測(cè)-校正方法帶來(lái)的時(shí)滯誤差，R.克拉夫等（2006）引入了α方法，并提出α-OS法，Chung等（1993）將其發(fā)展為Generalized-α-OS算法。除了上述方法，到目前為止針對(duì)子結(jié)構(gòu)擬動(dòng)力混合試驗(yàn)中的其它組合數(shù)值積分研究還較少。

從本質(zhì)上來(lái)說，CD-Newmark算法是一種組合方法，而PC-Newmark算法和α-OS算法是預(yù)測(cè)-校正方法，是將隱式方法拆成兩部分，一部分為顯式，另一部分為隱式，每進(jìn)行一步都要等到顯式部分計(jì)算完畢后才能計(jì)算隱式部分，最后將計(jì)算部分與試驗(yàn)部分進(jìn)行對(duì)比校正，導(dǎo)致時(shí)滯使誤差加大。M.Nakashima等（1985）得出的CD-Newmark算法的穩(wěn)定性條件為ωΔt≤2，但ω對(duì)應(yīng)于所有計(jì)算子結(jié)構(gòu)自由度被嵌固時(shí)對(duì)應(yīng)結(jié)構(gòu)體系的最高階圓頻率，導(dǎo)致對(duì)時(shí)間步長(zhǎng)Δt要求較為嚴(yán)格。

由上述分析可知：CD-Newmark算法時(shí)間步長(zhǎng)受限的原因，是由于計(jì)算子結(jié)構(gòu)的剛度是所有計(jì)算子結(jié)構(gòu)自由度被嵌固后的整體剛度。本文對(duì)CD-Newmark算法的缺點(diǎn)進(jìn)行參數(shù)修正，形成修正CD-Newmark算法，將計(jì)算子結(jié)構(gòu)的所有被嵌固自由度釋放，并考慮阻尼特性研究組合算法的數(shù)值特性。

基于剪切模型的結(jié)構(gòu)整體運(yùn)動(dòng)方程為：

式中，M、C和K分別為結(jié)構(gòu)質(zhì)量陣、阻尼陣和剛度陣；a、v和d分別為結(jié)構(gòu)加速度、速度和位移反應(yīng)向量；f為外荷載向量，上標(biāo)E、I分別代表試驗(yàn)子結(jié)構(gòu)和計(jì)算子結(jié)構(gòu)，下標(biāo)i、i+1分別代表本時(shí)刻和下一時(shí)刻。

試驗(yàn)子結(jié)構(gòu)位移和速度采用顯式中央差分法：

為了尋求穩(wěn)定性更好的計(jì)算子結(jié)構(gòu)隱式數(shù)值積分算法，這里不再采用Newmark-β算法中的平均加速度方法，即先不設(shè)控制參量β和γ的實(shí)數(shù)值，而是通過求解無(wú)條件穩(wěn)定性條件推導(dǎo)穩(wěn)定性更好的參數(shù)取值范圍，這樣計(jì)算子結(jié)構(gòu)為Newmark-β法的表達(dá)式為：

1.2 有阻尼組合算法的穩(wěn)定性分析

本文考慮簡(jiǎn)化后采用瑞利阻尼假設(shè)，由于穩(wěn)定性條件是針對(duì)任意初始條件的，因此考慮無(wú)外荷載的情況，將式（2）—（5）代入式（1）中，經(jīng)整理后得到相鄰兩步的狀態(tài)方程為：

其中，i時(shí)刻系統(tǒng)傳遞矩陣Yi為：

放大矩陣B為：

式中，b、c、d、e、f、g 為放大矩陣的積分系數(shù)，分別表示如下：

上式中ΩE、ΩI、ΩIE分別為試驗(yàn)子結(jié)構(gòu)、計(jì)算子結(jié)構(gòu)和相互作用部分的穩(wěn)定參數(shù)，ΩE=Δt·ωE，ΩI=Δt ·ωI，Ω2IE=Δt2·MI/KE；ξE和ξI分別為試驗(yàn)子結(jié)構(gòu)和計(jì)算子結(jié)構(gòu)的模態(tài)阻尼比，ξE=CE2MEωE，ξI=CI2MIωI。

數(shù)值積分方法的穩(wěn)定性條件為：

式中，ρ(B)是放大矩陣B的譜半徑，ρ(B)=maxλi，λi是放大矩陣B的特征值。

求解方法與無(wú)阻尼算法一致，通過求解B的特征方程得到特征根λ1=λ2=0，λ3=1，均能滿足穩(wěn)定性條件。

當(dāng)試驗(yàn)子結(jié)構(gòu)模態(tài)阻尼比ξE=0時(shí)，與無(wú)阻尼的修正CD-Newmark法得到的穩(wěn)定性條件ΩE≤2相同；當(dāng)ξE≠0時(shí)，為有阻尼的修正CD-Newmark法的試驗(yàn)子結(jié)構(gòu)穩(wěn)定條件。可見在試驗(yàn)子結(jié)構(gòu)中隨著ΩE的增大，阻尼比越小對(duì)穩(wěn)定性越有利。

當(dāng)式（17）中特征根60λ=時(shí)，子結(jié)構(gòu)的無(wú)條件穩(wěn)定性最優(yōu)，經(jīng)計(jì)算可得到控制參量β和γ的關(guān)系：

將00.5γ≤≤代入到式（18）中，可得：

從式（19）的數(shù)學(xué)形式上可以分析到與無(wú)阻尼的修正CD-Newmark法具有相似性，具體體現(xiàn)在控制參量β，γ，ξ1的關(guān)系上，如圖1所示，可以得出00.5γ≤≤，0.51β≤≤的穩(wěn)定性參數(shù)取值范圍。取1.0β=，0.5γ=的試驗(yàn)子結(jié)構(gòu)模態(tài)阻尼比Iξ變化對(duì)修正CD-Newmark法的穩(wěn)定性影響，如圖2所示，從圖中可以看出隨著Iξ的遞減其穩(wěn)定性越好。

圖1 ξI=0.05不同β和γ的譜半徑Fig. 1 Spectral radius of differentβandγatξI=0.05

圖2 不同Iξ的譜半徑Fig. 2 Spectral radius of differentIξ

1.3 有阻尼組合算法的周期失真率分析

修正CD-Newmark法中周期失真率可以表示為：

上式中F為特征根4,5λ的實(shí)部，E為特征根4,5λ的虛部：

將式（21）、（22）、（23）代入式（20）中，整理后可得周期失真率為：

本文研究的是組合算法，式（24）中需要加一個(gè)自定義調(diào)節(jié)參量δ1=ΩΩE，則周期失真率修正為：

從式（25）的數(shù)學(xué)形式可見，與無(wú)阻尼的修正CD-Newmark法相比，周期失真率僅在Ω中多出ξE，因此，δ1的性質(zhì)與無(wú)阻尼的修正CD-Newmark法的自定義調(diào)節(jié)參量δ性質(zhì)相同，均為隨著δ1的遞增周期失真率減小。ξE變化對(duì)有阻尼修正CD-Newmark法的周期失真率的影響如圖3所示。由圖3可知，前期的阻尼比越小周期失真率越低，雖然后期的周期失真率有穩(wěn)定的趨勢(shì)，但穩(wěn)定性范圍比無(wú)阻尼的有所減小，即圖中的橫坐標(biāo)試驗(yàn)子結(jié)構(gòu)穩(wěn)定參數(shù)ΩE≤1.4。為了保證良好的穩(wěn)定性和精度，本文建議在試驗(yàn)子結(jié)構(gòu)模態(tài)阻尼比0＜ξE≤0.05的周期失真率范圍為0%—8%，這樣精度較高。

1.4 考慮阻尼影響的算法數(shù)值阻尼比

數(shù)值模態(tài)阻尼比ξ可造成振幅衰減，和周期失真率一樣都會(huì)導(dǎo)致數(shù)值計(jì)算的不準(zhǔn)確，影響到組合方法的精度，阻尼比ξ的公式表達(dá)為：

圖4為修正CD-Newmark法的（27）式在不同阻尼比下的分析圖，可見阻尼比越小振幅衰減越小。為此本文建議取0＜ξE≤0.05，修正CD-Newmark法的數(shù)值阻尼比在0%—10%，這樣精度較高。

圖3 不同Eξ的11ξ=周期失真率Fig. 3 Period distortion rate of differentEξat11ξ=

圖4 不同Eξ的數(shù)值阻尼比Fig. 4 Numerical damping ratio of differentEξ

2 數(shù)值算例驗(yàn)證

2.1 穩(wěn)定性分析

選用《結(jié)構(gòu)動(dòng)力學(xué)》中的習(xí)題13.2（喬普拉，2009），假定模型為二層框架，一層為試驗(yàn)子結(jié)構(gòu)（質(zhì)量為2m，剛度為2k），二層為計(jì)算子結(jié)構(gòu)（質(zhì)量為m，剛度為k），其中：質(zhì)量m=100千磅力/g，側(cè)向剛度k=3154千磅力/英寸，ξn=5%，輸入Δt=0.02s的El Centro地面加速度記錄，加速度峰值調(diào)為0.4g。經(jīng)計(jì)算和單位轉(zhuǎn)換可得到一層質(zhì)量為44.5t，二層質(zhì)量為22.3t，一層剛度為5.51×105kN/m，二層剛度為2.76×105kN/m，一階自振頻率為78.8 rad/s，二階自振頻率為157.3rad/s。本文的修正CD-Newmark組合方法穩(wěn)定性條件為ΩE=ΔtωE＜2，滿足穩(wěn)定性的范圍；對(duì)于CD-Newmark法的Δtω=0.02×157.3=3.146，超出穩(wěn)定性范圍。圖5和圖6說明了超出有阻尼的算法穩(wěn)定范圍導(dǎo)致發(fā)散。圖7和圖8說明了穩(wěn)定性良好，驗(yàn)證了本文組合方法條件穩(wěn)定的推導(dǎo)結(jié)論。

2.2 精度分析

計(jì)算模型取自郭昕（2005），為兩層剪切型鋼筋混凝土框架結(jié)構(gòu)。質(zhì)量m1=1751kg，m2=3502kg；剛度K1=1400N/mm，K2=700N/mm；假定結(jié)構(gòu)處于彈性狀態(tài)，輸入El Centro地震波，峰值加速度調(diào)整為0.05g，初始位移、速度均為0，阻尼比為0.05。以底層為試驗(yàn)子結(jié)構(gòu)，二層為計(jì)算子結(jié)構(gòu)，經(jīng)計(jì)算可得到一階自振頻率為11.18 rad/s，二階自振頻率為35.65rad/s，時(shí)間步長(zhǎng)為0.01s，通過MATLAB軟件編程可得到該結(jié)構(gòu)的位移時(shí)程。

經(jīng)計(jì)算筆者得到了結(jié)構(gòu)底層位移和二層位移對(duì)比圖，如圖9和圖10所示。在本算例中通過與精確值的對(duì)比，檢驗(yàn)組合方法的精確性，精確解是利用狀態(tài)空間法計(jì)算的數(shù)值，狀態(tài)空間是一種既保證高精度，又完全考慮耦合的計(jì)算方法。從圖中可見，組合方法的精確性較高，驗(yàn)證了本文的組合方法由于考慮了耦合提高了算法的精確性。

圖5 CD-Newmark法的一層位移Fig. 5 Displacement of the first floor with CD-Newmark

圖6 CD-Newmark法的二層位移Fig. 6 Displacement of the second floor with CD-Newmark

圖7 修正CD-Newmark法的一層位移Fig. 7 Displacement of the first floor with MCD-Newmark

圖8 修正 CD-Newmark法的二層位移Fig. 8 Displacement of the second floor with MCD-Newmark

圖9 一層位移精度對(duì)比Fig. 9 Comparison of displacement accuracy for the first floor

圖10 二層位移精度對(duì)比Fig. 10 Comparison of displacement accuracy for the second floor

3 結(jié)論

（1）當(dāng)ξE≠0時(shí)，考慮阻尼的修正CD-Newmark法的試驗(yàn)子結(jié)構(gòu)穩(wěn)定條件為與有阻尼中心差分法的穩(wěn)定性條件相同，其穩(wěn)定性優(yōu)于CD-Newmark法，試驗(yàn)子結(jié)構(gòu)中隨著ΩE的增大，阻尼比越小對(duì)穩(wěn)定性越有利；計(jì)算子結(jié)構(gòu)的無(wú)條件穩(wěn)定性參數(shù)取值范圍和變化規(guī)律與無(wú)阻尼修正CD-Newmark法相同。

（2）考慮阻尼的修正CD-Newmark法在數(shù)學(xué)形式上與無(wú)阻尼的修正CD-Newmark法具有相似性，具體體現(xiàn)在β、γ、ξI的關(guān)系上，筆者建議取0≤γ≤0.5，而0.5≤β≤1是其穩(wěn)定性參數(shù)取值范圍。

（3）隨著ξI的遞減其穩(wěn)定性越好，筆者建議取0＜ξE≤0.05的周期失真率為0—8%，這樣精度較高。

（4）阻尼比越小振幅衰減越小，筆者建議取0＜ξE≤0.05，修正CD-Newmark法的數(shù)值阻尼比在0—10%，這樣精度較高。

（5）由于問題的復(fù)雜性，筆者沒有考慮在實(shí)時(shí)擬動(dòng)力混合試驗(yàn)的推廣，因此有關(guān)這一問題的后續(xù)研究，將在今后的研究工作中進(jìn)一步考慮。

程紹革，馬路，趙鵬飛，2008．基于凝聚技術(shù)的結(jié)構(gòu)抗震混合試驗(yàn)原理．工程抗震與加固改造，30（2）：62—67．

范云蕾，郭玉榮，肖巖等，2011．多層框架結(jié)構(gòu)遠(yuǎn)程協(xié)同擬動(dòng)力試驗(yàn)研究．土木工程學(xué)報(bào)，44（2）：28—35．

郭昕，2005．地震作用下多自由度結(jié)構(gòu)擬動(dòng)力實(shí)驗(yàn)方法研究 [碩士學(xué)位論文]．西安：西安建筑科技大學(xué)，36—39．

侯杰，劉鈞，杜文博等，2006．框架結(jié)構(gòu)擬動(dòng)力試驗(yàn)方法研究．工程抗震與加固改造，28（6）：66—70．

喬普拉，2009．結(jié)構(gòu)動(dòng)力學(xué)：理論及其在地震工程中的應(yīng)用．北京：清華大學(xué)出版社，515—516．

徐國(guó)山，2010．實(shí)時(shí)子結(jié)構(gòu)試驗(yàn)的等效力控制方法 [博士學(xué)位論文]．哈爾濱：哈爾濱工業(yè)大學(xué)，5—13．

王大鵬，田石柱，蔡新江，2010．網(wǎng)絡(luò)化協(xié)同結(jié)構(gòu)擬動(dòng)力試驗(yàn)方法與技術(shù)．應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，18（6）：910—921．

王強(qiáng)，馮帆，Shawn You等，2010．基于有限元軟件的子結(jié)構(gòu)擬動(dòng)力試驗(yàn)技術(shù)．沈陽(yáng)建筑大學(xué)學(xué)報(bào)（自然科學(xué)版），26（6）：1085—1089．

Chung J.，Hulber，G．，1993．A Time Integration Algorithm for Structural Dynamics with Improved Numerical Dissipation：The Generalized-Alpha Method．Journal of Applied Mechanics，60（2）：371—375．

M．Nakashima，H．Takai，1985．Use of Substructure Techniques in Pseudo-dynamic Testing．BRI Research Paper，No．111．

R．克拉夫，J．彭津，2006．結(jié)構(gòu)動(dòng)力學(xué)．北京：高等教育出版社，252—264．

Stability and Accuracy of Modified CD-Newmark Combined Numerical Integral Algorithms with Consideration of Damping Effect

Cai Xinjiang1，2），Wang Xiaofeng2）and Tian Shizhu1，2，3）

1）School of Civil Engineering，Suzhou University of Science and Technology，Suzhou 215011，China
2）Key Laboratory of Structure Engineering of Jiangsu Province，Suzhou 215011，China
3）School of Civil Engineering，Harbin Institute of Technology，Harbin 150090，China

It is the key to the substructure seismic pseudo-dynamic test to combine with advantages of explicit numerical integral algorithms without iteration and inexplicit numerical integral algorithms unconditionally stable．In this paper，a modified CD-Newmark combined numerical integral algorithm was established based on the parameter modification of traditional central difference algorithm and Newmark β algorithm．Considering the effect of damping，the analysis of the stable condition，the cycle distortion rate，the numerical damping ratio were carried out，then the parameters range of stable condition of experimental substructure and the unconditionally stable of calculated substructure are obtained，and the accuracy analysis is also studied．A case study is carried out to verify the algorithm，and then solve the strict stability problem of the CD-Newmark algorithm．Our study provides a reference for structural seismic pseudo dynamic test.

Numerical integration method；Combined algorithm；Central difference method；Newmark method；Stability；Accuracy

蔡新江，王曉峰，田石柱，2014．考慮阻尼影響的修正CD-Newmark組合算法數(shù)值特性研究．震災(zāi)防御技術(shù)，9（4）：898—906．

10.11899/zzfy20140418

國(guó)家自然科學(xué)基金（51308368，51278322）；省高校自然科學(xué)基金重大項(xiàng)目（13KJA560002）；省重點(diǎn)實(shí)驗(yàn)室開放課題（ZD1104）；住建部科技計(jì)劃項(xiàng)目（K2201213）；校科研基金（XKQ201203）；江蘇高校優(yōu)勢(shì)學(xué)科建設(shè)工程資助項(xiàng)目（PAPD）

2014-02-08

蔡新江，男，生于1981年。講師，博士。主要從事結(jié)構(gòu)抗震試驗(yàn)方法研究。E-mail：caixinjiang77@126.com