一族基于模型的積分算法在負剛度時的數(shù)值穩(wěn)定性

傅 博,馬 巖,陳 瑾

(長安大學(xué) 建筑工程學(xué)院,陜西 西安 710061)

0 引言

在進行結(jié)構(gòu)地震反應(yīng)的時程分析時,積分算法通常被用來求解結(jié)構(gòu)的運動方程。積分算法根據(jù)位移的表達形式可分為顯式和隱式。在進行非線性結(jié)構(gòu)動力學(xué)分析時,使用隱式算法需要進行迭代求解,而顯式算法則不需要。然而,顯式算法通常為條件穩(wěn)定。為了滿足穩(wěn)定性要求,條件穩(wěn)定算法要求積分時間步長與結(jié)構(gòu)體系最高自振頻率成反比,而無條件穩(wěn)定算法不需要對時間步長進行限制。為了實現(xiàn)較高的計算效率,學(xué)者們提出了一類新型積分算法—基于模型的積分算法[1],該類算法可以同時滿足無條件穩(wěn)定和顯式2個特性。例如,CHEN等[2]使用離散控制理論提出一種無條件穩(wěn)定且顯式的CR(Chen-Ricles, CR)算法。CHANG[3]通過改變CR算法的積分參數(shù),提出了一種改進的MCR(modified Chen-Ricles, MCR)算法。GUI等[4]在CR算法的積分參數(shù)中引入了一個新的參數(shù)λ,形成GUI算法。FU等[5]將2個新的參數(shù)κ1、κ2引入CR算法的積分參數(shù)中,形成一族新型積分算法(簡稱GCR算法),CR、MCR以及GUI算法均是GCR族算法的特例。

需要注意的是,基于模型的積分算法的無條件穩(wěn)定性是針對正剛度條件而非負剛度條件。而在結(jié)構(gòu)動力分析過程中,剛度可能出現(xiàn)負值的情況,這是因為許多結(jié)構(gòu)或材料在加載至屈服后都存在著變形增加荷載反而減小的現(xiàn)象。當(dāng)采用增量剛度的形式求解動力方程時,其恢復(fù)力骨架曲線常常出現(xiàn)負剛度的現(xiàn)象。程民憲[6-8]研究了在動力非線性分析中,負剛度條件下數(shù)值積分算法的收斂性和穩(wěn)定性問題,給出了相應(yīng)的判斷準則,并研究了平均常加速度法、二步Adams顯式方法、中心差分法、Wilson-θ法和Z-變換法等幾種數(shù)值積分算法在負剛度時的收斂性和穩(wěn)定性。吳云芳[9]給出了在負剛度條件下顯式積分法的收斂性的和無條件穩(wěn)定性。吳云芳[10]、XIE等[11]、李爽等[12]研究了Newmark族算法在負剛度條件下的收斂性和穩(wěn)定性,雖然研究結(jié)論存在一定的矛盾,但均證明平均常加速度法在負剛度條件下并非無條件穩(wěn)定。

分析積分算法的穩(wěn)定性時,通常采用單自由度體系,依據(jù)朱鏡清[13]和杜修力等[14]提出的結(jié)構(gòu)動力分析中數(shù)值穩(wěn)定性準則,算法的穩(wěn)定性可以從單自由度體系推廣到多自由度體系中。但在負剛度條件下,即便是單自由度體系為無條件穩(wěn)定的積分算法,用于多自由度體系時也可能會出現(xiàn)不穩(wěn)定。這是因為負剛度條件下,剛度矩陣為非正定矩陣,而非正定特性會引發(fā)一系列數(shù)值計算方面的問題。對于線性方程組的求解,杜修力等[15]指出:結(jié)構(gòu)進入負剛度后,結(jié)構(gòu)反應(yīng)將表現(xiàn)出失穩(wěn)狀態(tài)特征,剛度矩陣變成病態(tài)矩陣,即系數(shù)矩陣的微小誤差將會引起計算結(jié)果巨大的偏差,導(dǎo)致計算不收斂。文獻[16-20]給出了解決病態(tài)矩陣的方法,但這些方法每步都需要進行大量的迭代,這無疑嚴重影響了計算效率。因此,提出一種多自由度體系負剛度條件下無條件穩(wěn)定且無需求解病態(tài)方程的積分算法是很有意義的。

基于上述研究背景,本文分析了GCR族算法在負剛度條件下的數(shù)值穩(wěn)定性。同時,研究了GCR族算法的積分參數(shù)、阻尼比和時間步長等參數(shù)對算法數(shù)值穩(wěn)定性的影響規(guī)律。然后,分析了多自由度體系負剛度條件下隱式算法和GCR族算法的數(shù)值穩(wěn)定性,并提出適用于正負剛度條件下GCR族算法的分析流程和策略。最后,通過典型算例驗證了本文所提策略的有效性。

1 GCR族積分算法

單自由度結(jié)構(gòu)體系運動方程的一般形式為:

(1)

對于線彈性單自由度結(jié)構(gòu)體系,式(1)進行時間離散化后,可得到i+1時刻的運動方程:

(2)

式中,k為剛度。

為求解上述運動方程,可采用不同的積分算法。例如,經(jīng)典的Newmark族算法的表達式為:

(3a)

(3b)

式中:Δt為時間步長;γ和β為積分參數(shù)。Newmark族算法有幾個典型算法:[γ,β]=[1/2,1/4]時為平均常加速度法;[γ,β]=[1/2,1/6]為線性加速度法;[γ,β]=[1/2,0]為Newmark顯式算法。

本文研究的GCR族算法沿用CR算法的位移和速度表達式:

(4a)

(4b)

由式(4)發(fā)現(xiàn),GCR族算法的位移和速度表達式均為顯式。式(4)中α1和α2為積分參數(shù),其表達式為:

α1=(m+κ1Δtc+κ2Δt2k)-1m

(5a)

α2=(1/2+κ1)α1

(5b)

當(dāng)參數(shù)κ1和κ2滿足2κ2≥κ1≥1/2時,GCR族算法在正剛度條件下為無條件穩(wěn)定[5]。GCR族算法有幾個典型算法:當(dāng)[κ1,κ2]=[1/2,1/4]時為CR算法;[κ1,κ2]=[1/2,1/2]為MCR算法;[κ1,κ2]=[1/2,1/λ]為GUI算法。值得注意的是,GCR族積分算法與Newmark族算法除了表達式不同外,具有相同的數(shù)值特性,對應(yīng)的規(guī)則為:[γ,β]=[κ1,κ2]。

2 負剛度時GCR族算法的穩(wěn)定性分析

根據(jù)文獻[6],穩(wěn)定性問題包括系統(tǒng)的穩(wěn)定性和積分算法的數(shù)值穩(wěn)定性。系統(tǒng)受到初始擾動后,如果擾動隨時間增長而衰減,則為穩(wěn)定系統(tǒng);反之,當(dāng)擾動隨時間增長而增長,則為不穩(wěn)定系統(tǒng)。在負剛度條件下,擾動會隨時間增長而增長,系統(tǒng)是不穩(wěn)定的。

對于穩(wěn)定系統(tǒng),積分算法需要絕對穩(wěn)定,而對于不穩(wěn)定系統(tǒng),積分算法需要相對穩(wěn)定。積分算法的穩(wěn)定性、絕對穩(wěn)定和相對穩(wěn)定的定義如下:

對于穩(wěn)定系統(tǒng),當(dāng)j→∞,xj有界,且νj≤νi也有界,則該積分算法絕對穩(wěn)定。

對于不穩(wěn)定系統(tǒng),當(dāng)j→∞,xj無界,νj也無界,但|νj/xj|有界,則該積分算法相對穩(wěn)定。本文研究的負剛度條件下GCR族算法的穩(wěn)定性問題屬于相對穩(wěn)定。

對于非線性單自由度結(jié)構(gòu)體系,將恢復(fù)力分段線性化后,運動方程如式(6)所示:

(6)

式中:ξ為阻尼比;ω為自振頻率;δ=k1/k0為剛度比,k1和k0分別為屈服后剛度和初始剛度,線彈性階段δ=1,負剛度時δ<0;a為外激勵引起的加速度。

式(4)和式(6)可以用式(7)的矩陣遞推關(guān)系表示:

Yi+1=AYi+[0 0 -Δt2]Tai+1

(7)

(8)

式中,Ω=ωΔt。

矩陣A的特征值可以通過求其特征方程得到:

|A-λI|=λ3-2A1λ2+A2λ-A3=0

(9)

式中:I為3×3的單位矩陣;A1=[2-(δα2Ω2+2α1Ω)]/2=1/2矩陣A的跡;A2=(α1-α2)δΩ2-2α1ξΩ+1=矩陣A的余子式之和;A3=0=矩陣A的行列式值。A1和A2中的積分參數(shù)α1和α2分別為:α1=1/(1+2κ1ξΩ+δκ2Ω2),α2=(1/2+κ1)α1。

式(9)可以簡化為:

|A-λI|=λ2-2A1λ+A2=0

(10)

(11)

因為λ3=0,因此GCR族積分算法在負剛度穩(wěn)定性準則:ρ(A)=max(|λ1|,|λ2|)≤eμΩ。

因此,為了保證GCR族積分算法的負剛度下穩(wěn)定,應(yīng)滿足以下要求:

(12)

由式(12)可以給出GCR族算法在負剛度時無條件穩(wěn)定的積分參數(shù)取值范圍:

0≤κ1≤1/2;κ2≤0

(13)

因為GCR族算法與Newmark族算法具有相同的數(shù)值特性,所以以上結(jié)論可以推廣至Newmark族算法,即:當(dāng)0≤γ≤1/2、β≤0時, Newmark族算法在負剛度時無條件穩(wěn)定。這與前期研究的結(jié)論一致:當(dāng)γ=1/2、β=0時, Newmark顯式算法在負剛度時為無條件穩(wěn)定[8];當(dāng)γ=1/2、β=1/4時,平均常加速度法在負剛度時并非無條件穩(wěn)定[10-12]。

為了驗證GCR族算法在負剛度條件下滿足積分參數(shù)0≤κ1≤1/2、κ2≤0時為無條件穩(wěn)定,第3節(jié)將進行穩(wěn)定性影響參數(shù)分析,通過選取不同的積分參數(shù)、阻尼比和時間步長,并繪制譜半徑圖,從而驗證算法穩(wěn)定性。

3 穩(wěn)定性影響參數(shù)分析

3.1 積分參數(shù)的影響

在正剛度情況下,穩(wěn)定性分析時時間步長一般選用Δt/T=1/20～1/10,T為結(jié)構(gòu)周期,所以本文也采用該時間步長范圍的上下限。由圖1可知,當(dāng)κ2=0、ξ=0.05時,ρ(A)/eμΩ隨κ1(0≤κ1≤1/2)和δ的變化規(guī)律,可以發(fā)現(xiàn),當(dāng)Δt/T=1/20和Δt/T=1/10時,ρ(A)/eμΩ均小于等于1,說明GCR族算法無條件穩(wěn)定。在κ2、ξ和Δt相同的條件下,κ1越趨近于1/2,ρ(A)/eμΩ越趨近于1。

圖1 κ2=0、ξ=0.05時ρ(A)/eμΩ隨κ1(0≤κ1≤1/2)和δ的變化

由圖2可知:κ2=0、ξ=0.05時,ρ(A)/eμΩ隨κ1(κ1>1/2)和δ的變化規(guī)律;κ1>1/2時,存在ρ(A)/eμΩ大于1的情況,所以當(dāng)κ1>1/2時,GCR族算法在負剛度的條件下不是無條件穩(wěn)定。

圖2 κ2=0、ξ=0.05時ρ(A)/eμΩ隨κ1(κ1>1/2)和δ的變化

由圖3可知:κ1=1/2、ξ=0.05時,ρ(A)/eμΩ隨κ2(κ2≤0)和δ的變化規(guī)律;當(dāng)Δt/T=1/20和Δt/T=1/10時,ρ(A)/eμΩ均小于等于1,說明GCR族算法無條件穩(wěn)定。在κ1、ξ和Δt相同的條件下,κ2越趨近于0,ρ(A)/eμΩ越趨近于1。

圖3 κ1=1/2、ξ=0.05時ρ(A)/eμΩ隨κ2(κ2≤0)和δ的變化

由圖4可知:κ1=1/2、ξ=0.05時,ρ(A)/eμΩ隨κ2(κ2>0)和δ的變化規(guī)律;κ2>0時,存在ρ(A)/eμΩ大于1的情況,所以,當(dāng)κ2>0時,GCR族算法在負剛度的條件下不是無條件穩(wěn)定。

圖4 κ1=1/2、ξ=0.05時ρ(A)/eμΩ隨κ2(κ2>0)和δ的變化

3.2 阻尼比影響

由圖5可知,κ1=1/2、κ2=0時ρ(A)/eμΩ隨ξ和δ的變化規(guī)律?？梢钥闯霾煌枘岜葧r,ρ(A)/eμΩ均小于等于1,說明GCR族算法無條件穩(wěn)定。隨著ξ的增大,ρ(A)/eμΩ有略微增大的趨勢。

圖5 κ1=1/2、κ2=0時ρ(A)/eμΩ隨ξ和δ的變化

3.3 時間步長影響

由圖6可知,κ1=1/2、κ2=0和ξ=0.05時,ρ(A)/eμΩ隨Δt/T和δ的變化規(guī)律?？梢钥闯霾煌/T時,ρ(A)/eμΩ均小于等于1,說明GCR族算法無條件穩(wěn)定。隨著Δt/T的增大,ρ(A)/eμΩ越來越小。

圖6 κ1=1/2、κ2=0、ξ=0.05時ρ(A)/eμΩ隨Δt/T和δ的變化

4 負剛度時多自由度體系算法穩(wěn)定性

4.1 隱式算法穩(wěn)定性分析

隱式算法用于非線性動力分析時,通常有2種分析方法:一是采用迭代分析,這種方式耗費額外的計算時間,影響計算效率;二是將隱式改為顯式,這樣轉(zhuǎn)換的目的是采用其他反應(yīng)量來表示加速度,該方法在多自由度體系正剛度條件下應(yīng)用是沒有問題的,但是負剛度條件下可能會發(fā)生不穩(wěn)定現(xiàn)象。

以經(jīng)典Newmark族隱式算法(β≠0)為例,解釋上述不穩(wěn)定的原因。Newmark族算法用于多自由度體系時的運動方程及速度、位移增量表達式為:

(14a)

(14b)

(14c)

(15a)

(15b)

再將式(15)代入增量平衡方程式(14a),可得:

(16)

(17a)

(17b)

4.2 GCR族算法穩(wěn)定性分析

GCR族算法用于多自由度體系時的速度、位移及加速度表達式為:

(18a)

(18b)

(18c)

式中,R為恢復(fù)力向量。

GCR族算法的位移表達式為顯式,在動力分析過程中根據(jù)式(18)依次可以求出第i+1步的速度、位移和加速度,不會因為剛度矩陣K在負剛度條件下為非正定矩陣而產(chǎn)生病態(tài)矩陣求解過程,因此不會出現(xiàn)不穩(wěn)定現(xiàn)象。此外,前文已經(jīng)證明GCR族算法在積分參數(shù)為0≤κ1≤1/2;κ2≤0時,負剛度條件下為無條件穩(wěn)定。

當(dāng)2κ2≥κ1≥1/2時,GCR族算法在正剛度條件下無條件穩(wěn)定;而當(dāng)0≤κ1≤1/2;κ2≤0時,GCR族算法在負剛度條件下為無條件穩(wěn)定。不難發(fā)現(xiàn),并不存在一個參數(shù)取值范圍可以使得GCR族算法同時滿足正、負剛度無條件穩(wěn)定。對于多自由度體系剛度軟化過程,可能會出現(xiàn)正負剛度共存的情況,此時需要優(yōu)先滿足負剛度條件下的無條件穩(wěn)定,而正剛度條件下只能滿足條件穩(wěn)定。根據(jù)文獻[5],當(dāng)κ1≥1/2;2κ2≤κ1時,GCR族算法在正剛度時為條件穩(wěn)定,穩(wěn)定限值Ωcrit的表達式為:

(19)

圖7 時程分析流程圖

5 數(shù)值算例

5.1 無阻尼單自由度結(jié)構(gòu)體系

圖8 結(jié)構(gòu)位移響應(yīng)

5.2 有阻尼多自由度結(jié)構(gòu)體系

采用GCR族算法和綜合法[21]進行求解,綜合法視為精確解。GCR族算法的時間步長采用Δt=0.002 s,綜合法時間步長采用Δt=0.001 s。GCR族算法在正剛度條件下均采用κ1=1/2,κ2=1/4;負剛度條件下κ1、κ2的取值如圖9所示。結(jié)構(gòu)進入負剛度后,因負剛度的解是發(fā)散的,位移會快速增大,從而可能導(dǎo)致模型破壞,所以當(dāng)兩層框架結(jié)構(gòu)其中一層的剛度比δ小于-10時就終止分析。圖9為結(jié)構(gòu)每一層的位移響應(yīng)及實時剛度比?？梢钥闯?通過本文的策略可以使得GCR族算法在正負剛度均保持穩(wěn)定。

圖9 位移響應(yīng)和實時剛度比(θ=-3,κ1=1/2)

6 結(jié)論

1)當(dāng)積分參數(shù)滿足0≤κ1≤1/2,κ2≤0時,GCR族算法在負剛度條件下為無條件穩(wěn)定。

2)GCR族算法和Newmark族算法是對應(yīng)的關(guān)系,所以GCR族算法在負剛度無條件穩(wěn)定的參數(shù)條件亦可以用于Newmark族算法中,即:當(dāng)0≤γ≤1/2、β≤0時,Newmark族算法在負剛度條件下為無條件穩(wěn)定。

3)在多自由度體系負剛度條件下,隱式算法轉(zhuǎn)換為顯式再分析的方式會出現(xiàn)病態(tài)方程,算法可能會不穩(wěn)定。

4)GCR族算法的位移表達式是顯式,因此在多自由度體系負剛度條件下不會出現(xiàn)病態(tài)方程。

5)在多自由度體系同時存在正負剛度的情況下,通過判斷實時剛度比的正負,來選取合適的κ1和κ2的值,必要時需要犧牲正剛度時的無條件穩(wěn)定,優(yōu)先滿足負剛度時的無條件穩(wěn)定。