孟凡濤， 趙建鋒

(1. 青島理工大學(xué) 土木工程學(xué)院，山東 青島 266033； 2. 山東華科規(guī)劃建筑設(shè)計有限公司，山東 聊城 252000)

在結(jié)構(gòu)動力計算中，數(shù)值積分方法被廣泛應(yīng)用，其中較經(jīng)典的方法是Newmark-β法[1]、 Wilson-θ法[2]以及Bathe隱式方法[3]等，這些方法精度都很高且均具備較好的穩(wěn)定性。但是由于隱式方法求解大型方程組，尤其對于非線性問題，需反復(fù)迭代計算量較大，因此在實際應(yīng)用中耗時較多，有其局限性。

Zhai[4]在Newmark-β方法基礎(chǔ)上提出了新型快速顯式積分法，該方法在求解大型非線性動力學(xué)問題時，僅要求質(zhì)量矩陣M為對角陣，不限制阻尼矩陣和剛度矩陣的形式，不需要求解高階線性代數(shù)方程組，因而提高了數(shù)值計算的速度，其截斷誤差與Newmark-β法同階，該方法具有較好的穩(wěn)定條件，且當(dāng)參數(shù)ψ和φ相等且均取為1/2的情況下具有無振幅衰減及周期延長率變化小的優(yōu)點。但是以上提到的算法沒有考慮與結(jié)構(gòu)動力特性協(xié)同的問題[5]，所以上述方法不具備結(jié)構(gòu)動力計算的自適應(yīng)性。

Chang等[6-9]提出了考慮結(jié)構(gòu)動力特性協(xié)同的無條件穩(wěn)定的顯式數(shù)值積分方法(簡稱Chang方法)，但其算法是在平均加速度方法的基礎(chǔ)上經(jīng)過一次內(nèi)嵌的牛頓迭代來實現(xiàn)位移顯式表達(dá)，其方法僅僅是位移考慮了與結(jié)構(gòu)動力特性的自適應(yīng)性。Chang方法隨著其位移表達(dá)式中的積分參數(shù)不斷改進(jìn)，其算法性能也再不斷提升，并發(fā)展了一種帶耗散特性的算法，但其速度項仍不具有考慮結(jié)構(gòu)動力特性協(xié)同的自適應(yīng)性。隨著快速擬動力試驗技術(shù)及實時子結(jié)構(gòu)試驗技術(shù)的發(fā)展，顯式數(shù)值積分算法更加突顯出其優(yōu)勢。Chen等[10-11]利用控制理論[12]提出了一種考慮結(jié)構(gòu)動力特性協(xié)同的無條件穩(wěn)定的顯式積分算法。其后Gui等[13-14]在其基礎(chǔ)上引入了控制周期延長率的參數(shù)，使得該種算法更加完善。Rezaiee-Pajand等[15]利用控制理論在翟婉明提出的算法的基礎(chǔ)上發(fā)展了USE(Unconditional Stability Explicit )算法，該種算法可以考慮結(jié)構(gòu)的動力特性協(xié)同問題，但其不能控制周期延長率的變化。遺憾的是Chang方法、Chen等的CRM(Chen and Ricles Method)以及USE算法均存在異常的振幅增長現(xiàn)象，需引入外荷載項的差分來補救和消除這種現(xiàn)象[16]。

Namadchi等[17]利用Bathe隱式算法的放大矩陣提出了半顯式的SEUSBCM算法。本文則結(jié)合控制理論在Bathe隱式算法的基礎(chǔ)上利用CRM的位移和速度表達(dá)式發(fā)展了一種帶耗散特性的雙顯式無條件穩(wěn)定的新算法。因此，可稱本文方法為(Explicit Unconditionally Stable Time Integration Based on the Composite Scheme Method，EUSBCM)。EUSBCM具有與Bathe隱式復(fù)合積分算法相同的數(shù)值特性，但EUSBCM是一種單步積分算法，不同與Bathe隱式方法積分步內(nèi)具有子積分步的復(fù)合積分。本文對EUSBCM的穩(wěn)定性和精度包括數(shù)值耗散和色散均進(jìn)行了分析，以及與其它積分算法的對比工作，并給出了應(yīng)用于MDOF(Multi-Degree-of-Freedom System)時兩個積分參數(shù)α2，α1的推導(dǎo)過程及表達(dá)式。最后線性和非線性數(shù)值算例驗證了本文所給出的方法的有效性和正確性。EUSBCM具有較強(qiáng)的耗散特性，根據(jù)實時子結(jié)構(gòu)試驗的具體試驗研究[18-19]并結(jié)合筆者以前對實時子結(jié)構(gòu)試驗算法的研究[20-21]，本文算法在實時子結(jié)構(gòu)試驗中應(yīng)用具有更強(qiáng)的優(yōu)勢。

1 本文顯式算法的提出

單自由度彈性系統(tǒng)的動力方程如式(1)所示，本文的數(shù)值積分算法位移、速度表達(dá)式分別如式(2)、式(3)所示。對于式(2)和式(3)中的參數(shù)α2，α1，本文將利用Bathe隱式算法的放大矩陣進(jìn)行求解。

mai+1+cvi+1+kdi+1=fi+1

(1)

di+1=di+Δtvi+α2Δt2ai

(2)

vi+1=vi+α1Δtai

(3)

由式(1)～式(3)按單自由度體系進(jìn)行離散系統(tǒng)的z變換，得到的傳遞函數(shù)如式(4)所示。對于Bathe隱式算法的放大矩陣為式(5)所示。放大矩陣所對應(yīng)的特征值的行列式如式(6)所示。將式(6)化簡后可得到式(7)。聯(lián)立式(4)中的分母項中的一次項系數(shù)和常數(shù)項與式(7)中的A1和A2可以得到方程組如式(11)所示。利用Matlab可求得參數(shù)α1，α2如式(12)～式(13)所示。

(4)

(5)

|[A]n+1-λ[I]|=0

(6)

λ3-A1λ2+A2λ-A3=0

(7)

(8)

(9)

A3=0

(10)

(11)

(12)

(13)

2 本文顯式算法的特性分析

2.1 線性與非線性穩(wěn)定性分析

由于EUSBCM的放大矩陣與Bathe隱式方法的相同，因此EUSBCM會有與Bathe隱式方法相同的部分?jǐn)?shù)值特性，由于放大矩陣相同，因此EUSBCM對線彈性系統(tǒng)是無條件穩(wěn)定的，但EUSBCM與Bathe隱式方法最大的不同在于EUSBCM不需要在每個時間步內(nèi)再進(jìn)行子分步的計算；為方便分析問題，定于第n+1時步的剛度kn+1與結(jié)構(gòu)初始時刻的剛度k0具有δn+1=kn+1/k0的關(guān)系，當(dāng)δn+1>1時，即為剛度硬化。對于剛度硬化系統(tǒng)，Bathe隱式方法仍是無條件穩(wěn)定的算法，但EUSBCM是具有條件穩(wěn)定性的。為了便于說明EUSBCM的穩(wěn)定性，將Bathe隱式算法、SEUSBCM，EUSBCM，YU implicit M[22]、MKRM[23]進(jìn)行了比較，如圖1所示，EUSBCM具有與Bathe隱式算法、SEUSBCM方法相同的譜半徑。對于非線性系統(tǒng)中的穩(wěn)定性，按照不同的剛度變化進(jìn)行比較，如圖2所示，對于剛度軟化系統(tǒng)，EUSBCM是無條件穩(wěn)定的算法，對于剛度硬化系統(tǒng)，EUSBCM是條件穩(wěn)定的算法。離散控制理論采用閉環(huán)系統(tǒng)分析非線性問題的穩(wěn)定性，直接積分法閉環(huán)系統(tǒng)的結(jié)構(gòu)圖，如圖3所示。

圖1 不同算法的穩(wěn)定性比較Fig.1 Comparison of spectral curves for various methods

圖2 非線性系統(tǒng)的譜半徑比較Fig.2 Variation of spectral in nonlinear systems

圖3 剛度變化時閉環(huán)系統(tǒng)圖示Fig.3 Closed-loop block diagram for system with nonlinear stiffness

將式(1)按照單自由度增量形式表述為

mΔai+cΔvi=Δfi-Δri=ΔLi

式中：Δri=ktΔdi根據(jù)離散控制理論，閉環(huán)系統(tǒng)的傳遞函數(shù)可寫為

(14)

其中，

式(14)中分母為零對應(yīng)的特征方程為式(15)。

(15)

對于非線性系統(tǒng)，將ξ=0.02，ωn=2π rad/s，Δt=0.5 s，m=1代入式(15)，利用Matlab繪制其根軌跡如圖4所示。兩條根軌跡中有一個分支在z=-1處穿出單位圓并沿著負(fù)實軸方向發(fā)展，而另一個分支則沒有跨越單位圓，將z=-1代入式(15)得到穩(wěn)定界限如式(16)所示。根據(jù)式(16)可以得出隨剛度比δn+1i增大，不同阻尼比時的穩(wěn)定界限如圖5所示，可直觀的得出當(dāng)δn+1i≤1，算法是無條件穩(wěn)定的，當(dāng)δn+1>1時，算法存在穩(wěn)定上限，因此算法變成了有條件穩(wěn)定。

(16)

圖4 非線性系統(tǒng)的根軌跡Fig.4 Root locus of nonlinear system

圖5 隨剛度比變化不同阻尼比時的穩(wěn)定界限Fig.5 Variation of upper stability limit versus δn+1i for different viscous damping ratios

2.2 精度分析

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

2.3 超調(diào)性分析

對于初始條件不為零的單自由度系統(tǒng)(Single-Degree-of-Freedom System,SDOF)，可能會由于時間步長較大而在計算初始時間階段出現(xiàn)系統(tǒng)的響應(yīng)被放大的現(xiàn)象，這種現(xiàn)象稱為超調(diào)。對于一種直接積分算法的超調(diào)特性，一般分析一個不考慮阻尼的自由振動系統(tǒng)在無外荷載、初始位移和初始速度不為零的條件下，位移和速度在第一個時間步長內(nèi)隨Ω的變化趨勢。本文算法在第一個時間步長內(nèi)的位移和速度的表達(dá)式如式(25)所示。當(dāng)Ω→∞時由式(12)、式(13)及式(25)可以得到本文算法位移和速度均無超調(diào)特性。

(25)

圖6 不同算法的PEFig.6 Comparison of relative period error for various methods

圖7 不同算法的算法阻尼比Fig.7 Comparison of algorithmic damping ratio for various methods

圖8 非線性系統(tǒng)下的PEFig.8 Variation of relative period error in nonlinear systems

圖9 非線性系統(tǒng)下的算法阻尼比Fig.9 Variation of algorithmic damping ratio in nonlinear systems

3 應(yīng)用于MDOF時[α2]和[α1]的確定

為了得到多自由度系統(tǒng)計算時的[α2]和[α1]，采用模態(tài)分解技術(shù)將多自由度系統(tǒng)分解為多個單自由度系統(tǒng)，則在模態(tài)空間令由式(26)和式(27)成立。式(26)中分子為式(27)所示，分母為式(28)所示。式(29)中分子為式(30)所示，分母為式(31)所示。在模態(tài)坐標(biāo)系中自由振動方程可寫為式(32)。將式(32)化簡可以得到式(33)。根據(jù)式(33)中的參數(shù)代入式(26)可以得到式(34)和式(35)成立?；喪?34)和式(35)可以得到式(36)和式(37)，進(jìn)一步化簡可以得到式(38)和式(39)成立，經(jīng)進(jìn)一步整理最終得到式(40)和式(41)。將EUSBCM的位移表達(dá)式寫為模態(tài)坐標(biāo)的形式為式(42)，將模態(tài)坐標(biāo)形式轉(zhuǎn)化為普通形式的位移表達(dá)式為式(43),式(43)的參數(shù)可得到式(44)成立，由式(44)可以得到式(45)進(jìn)一步可以得到式(46)，將式(36)和式(37)代入式(46)得到式(47)。同理，將EUSBCM的速度表達(dá)式寫為模態(tài)坐標(biāo)的形式為式(48)，將模態(tài)坐標(biāo)形式轉(zhuǎn)化為普通形式的速度表達(dá)式為式(49)，式(49)的參數(shù)可得到式(50)成立，由式(50)可以得到式(51)，進(jìn)一步可以得到式(52)，將式(40)和式(41)代入式(52)得到式(53)。

由以上推導(dǎo)過程可以得出，在確定多自由度體系的積分參數(shù)[α2]和[α1]時，不需處理特征值的計算問題，只要已知結(jié)構(gòu)的初始剛度矩陣，質(zhì)量矩陣和確定的阻尼比，則可直接利用式(47)和式(53)的結(jié)果，該兩個參數(shù)僅需計算前確定一次即可，即使非線性計算也無需每個時間步內(nèi)每次都要重新確定該兩個參數(shù)，即該兩個參數(shù)僅與結(jié)構(gòu)的彈性初始狀態(tài)的動力特性和所選定的時間積分步長相關(guān)。

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

36Δt([Φ]-1[M]-1[C][Φ])+144[I]

(34)

Δt2[Φ]-1[M]-1[K]0[Φ])×

(16[I]+8Δt[Φ]-1[M]-1[C][Φ]+

Δt2[Φ]-1[M]-1[K]0[Φ])

(35)

36Δt([M]-1[C])+144[I]

(36)

(16[I]+8Δt[M]-1[C]+Δt2[M]-1[K]0)(37)

24Δt([Φ]-1[M]-1[C][Φ])+144[I]

(38)

(16[I]+8Δt[Φ]-1[M]-1[C][Φ]+Δt2[Φ]-1[M]-1[K]0[Φ])

(39)

24Δt([M]-1[C])+144[I]

(40)

(16[I]+8Δt[M]-1[C]+Δt2[M]-1[K]0)

(41)

(42)

(43)

(44)

(45)

(46)

[α2]=[(9[I]+6Δt[M]-1[C]+Δt2[M]-1[K]0)×

(16[I]+8Δt[M]-1[C]+Δt2[M]-1[K]0)]-1×

[2Δt2[M]-1[K]0+36Δt[M]-1[C]+144[I]]

(47)

(48)

(49)

(50)

(51)

(52)

[α1]=[(9[I]+6Δt[M]-1[C]+Δt2[M]-1[K]0)×

(16[I]+8Δt[M]-1[C]+Δt2[M]-1[K]0)]-1×

[Δt2[M]-1[K]0+24Δt[M]-1[C]+144[I]]

(53)

4 數(shù)值算例

算例1 分析一個單自由度無阻尼線性系統(tǒng)的自由振動，其運動微分方程如式(54)所示。自振頻率ω=2π，初始位移為1 m，初始速度為0。位移響應(yīng)的精確解為x(t)=cos(2πt)，其中時間單位為秒。取積分步長Δt=0.01 s，計算結(jié)果與解析解進(jìn)行比較，對比0～20 s位移時程曲線，由圖10能夠很直觀的得到，所得到的位移與解析解吻合較好。

(54)

圖10 Δt=0.01 s時的自由振動位移響應(yīng)時程曲線Fig.10 Time histories of displacement response with Δt=0.01 s

(55)

(56)

圖11 質(zhì)點1的位移響應(yīng)時程曲線Fig.11 Time histories of displacement response of particle 1

圖12 質(zhì)點2的位移響應(yīng)時程曲線Fig.12 Time histories of displacement response of particle 2

算例3 一單自由度體系結(jié)構(gòu)質(zhì)量為39 000 kg 水平剛度為140 000 kN/m，阻尼比為0.05，屈服后與屈服前的剛度比為0.18，恢復(fù)力模型為Bouc-Wen模型，該模型的數(shù)學(xué)表達(dá)式如式(57)、式(58)所示。

fs(u,z)=aku+(1-α)kz

(57)

(58)

圖13 結(jié)構(gòu)的位移反應(yīng)Fig.13 Displacement response of the structure

圖14 結(jié)構(gòu)的恢復(fù)力曲線Fig.14 Restoring force curve of the structure

5 結(jié) 論

本文基于Bathe隱式算法，利用CRM的位移和速度表達(dá)式發(fā)展了一種帶耗散特性的顯式無條件穩(wěn)定的新算法EUSBCM。該種算法對剛度軟化系統(tǒng)是無條件穩(wěn)定的，對于剛度硬化系統(tǒng)是條件穩(wěn)定的。本文基于離散控制理論對于非線性硬化系統(tǒng)的穩(wěn)定界限給出了理論求解公式。本文算法的位移和速度均是顯式的，且僅有兩個與結(jié)構(gòu)初始彈性狀態(tài)及時間積分間隔相關(guān)的參數(shù)，而SEUSBCM具有四個參數(shù)且速度是隱式格式的表達(dá)形式，因此筆者提出的EUSBCM的計算效率要高于SEUSBCM。EUSBCM具有較強(qiáng)的高頻耗散特性，這是CRM算法所不具有的特性，因此EUSBCM更適合應(yīng)用于實時子結(jié)構(gòu)試驗，且EUSBCM不具有異常的超調(diào)現(xiàn)象，無需采用外荷載項的差分進(jìn)行補救。最后數(shù)值算例均表明EUSBCM具有良好的計算穩(wěn)定性且所計算的結(jié)果均與精確解或經(jīng)典隱式算法的計算結(jié)果吻合良好。