基于Bathe隱式算法的結構動力學顯式算法

孟凡濤， 趙建鋒

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

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

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

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

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

1 本文顯式算法的提出

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

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

(1)

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

(2)

vi+1=vi+α1Δtai

(3)

由式(1)～式(3)按單自由度體系進行離散系統的z變換，得到的傳遞函數如式(4)所示。對于Bathe隱式算法的放大矩陣為式(5)所示。放大矩陣所對應的特征值的行列式如式(6)所示。將式(6)化簡后可得到式(7)。聯立式(4)中的分母項中的一次項系數和常數項與式(7)中的A1和A2可以得到方程組如式(11)所示。利用Matlab可求得參數α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 線性與非線性穩定性分析

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

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

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

圖3 剛度變化時閉環系統圖示Fig.3 Closed-loop block diagram for system with nonlinear stiffness

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

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

式中：Δri=ktΔdi根據離散控制理論，閉環系統的傳遞函數可寫為

(14)

其中，

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

(15)

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

(16)

圖4 非線性系統的根軌跡Fig.4 Root locus of nonlinear system

圖5 隨剛度比變化不同阻尼比時的穩定界限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 超調性分析

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

(25)

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

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

圖8 非線性系統下的PEFig.8 Variation of relative period error in nonlinear systems

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

3 應用于MDOF時[α2]和[α1]的確定

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

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

(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 數值算例

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

(54)

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

(55)

(56)

圖11 質點1的位移響應時程曲線Fig.11 Time histories of displacement response of particle 1

圖12 質點2的位移響應時程曲線Fig.12 Time histories of displacement response of particle 2

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

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

(57)

(58)

圖13 結構的位移反應Fig.13 Displacement response of the structure

圖14 結構的恢復力曲線Fig.14 Restoring force curve of the structure

5 結 論

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