基于修正四次樣條插值法的結(jié)構(gòu)動力學(xué)顯式算法

趙 旭，趙建鋒

(青島理工大學(xué) 土木工程學(xué)院，青島 266525)

在結(jié)構(gòu)動力分析中，數(shù)值積分法以其適用性廣和計算機(jī)輔助分析方便而得到廣泛應(yīng)用。Newmark法、Wilson-θ法、中心差分法、HHT-α法和CR法等經(jīng)典算法均屬于數(shù)值積分算法。數(shù)值積分法按照求解方式的不同可分為隱式法和顯式法兩類。隱式算法，如Newmark法和Wilson-θ法，指不能由當(dāng)前已知時間步的狀態(tài)量(位移、速度、加速度)直接計算得到下一時步的狀態(tài)量，需迭代計算；否則為顯式算法，如中心差分法和CR法。隱式算法一般為無條件穩(wěn)定，但增加步長會影響精度，且計算量大。顯式算法不需迭代計算，但一般為條件穩(wěn)定，積分步長受到限制，穩(wěn)定性和精度較差。隨著結(jié)構(gòu)動力分析的復(fù)雜性和實時混合試驗的發(fā)展，數(shù)值積分算法的高穩(wěn)定性、高精度和高計算效率越來越受到重視。因此，國內(nèi)外學(xué)者對無條件穩(wěn)定顯式算法的研究取得較多成果。

CHANG[1-3]提出無條件穩(wěn)定的顯式數(shù)值積分算法，并對積分參數(shù)不斷改進(jìn)，算法性能也不斷提升。文穎等[4]基于加速度泰勒展開式，提出高精度的顯式算法，但該算法為條件穩(wěn)定。楊超等[5]以加速度為基本變量，提出一類具有非耗散特性的顯式積分方法，并且可通過調(diào)整參數(shù)改變算法的穩(wěn)定性。冉田苒等[6]基于隱式HHT-α法提出無條件穩(wěn)定的顯式新算法，但需求解方程組得到不同時刻的狀態(tài)量。KOLAY等[7]利用離散控制理論在廣義α法的基礎(chǔ)上發(fā)展了無條件穩(wěn)定的顯式算法——KR-α法。SHOJAEE等[8]提出了修正的四次樣條插值積分算法(下文稱為MQBS)，該方法為無條件穩(wěn)定，有三階精度，但為隱式算法，對于非線性問題需反復(fù)迭代，計算效率不高，且需要求加速度的高階導(dǎo)數(shù)。

離散控制理論設(shè)計結(jié)構(gòu)動力學(xué)算法與傳遞矩陣分析有相同的效果，且使算法的設(shè)計與推導(dǎo)更加合理簡便[9]。本文基于MQBS法的傳遞矩陣和CR法的位移速度遞推式，利用離散控制理論的Z變換和離散傳遞函數(shù)，設(shè)計了一種無條件穩(wěn)定的顯式算法(稱新算法為NDEM)。NDEM法可通過調(diào)節(jié)系數(shù)改變算法周期延長率進(jìn)而調(diào)節(jié)算法的精度，具有零振幅衰減率，最高可滿足三階精度。從理論上分析新設(shè)計算法的穩(wěn)定性和精度(包括周期延長率和振幅衰減率)，給出了非線性硬化系統(tǒng)的穩(wěn)定性界限，并通過算例證明了所設(shè)計算法的有效性。

1 算法推導(dǎo)

1.1 運動微分方程

在第i時刻，線性單自由度(SDOF)系統(tǒng)的運動微分方程如式(1)所示，通常求解此微分方程需要知道位移和速度遞推式，即式(2)和式(3)。

(1)

(2)

(3)

(4)

det(A-λI)=λ4-A1λ3+A2λ2-A3λ+A4=0

(5)

(6)

式中：A1為傳遞矩陣的跡；A2為傳遞矩陣二階主子式的和；A3為傳遞矩陣三階主子式的和；A4為傳遞矩陣的行列式。

將式(6)代入式(5)可整理得到傳遞矩陣特征值p為

(7)

式中：Ω=Δt·ω，Δt為積分時間步長，ω為系統(tǒng)固有頻率。

1.2 離散控制理論

離散控制理論設(shè)計動力學(xué)算法需要用到Z變換和離散傳遞函數(shù)[10]。式(1)—(3)中各狀態(tài)量的Z變換如式(8)(9)所示，F(xiàn)(z)和X(z)分別為系統(tǒng)不同時刻輸入量和輸出量的Z變換。根據(jù)Z變換的位移定理，可進(jìn)行時步之間的轉(zhuǎn)換，見式(10)。

(10)

離散傳遞函數(shù)G(z)為開環(huán)系統(tǒng)傳遞函數(shù)，由系統(tǒng)輸出量與輸入量的Z變換X(z)和F(z)的比值表示。

(11)

式中：n，d為系數(shù)。

將式(9)和式(10)代入式(1)—(3)可得離散傳遞函數(shù)G(z)的各項系數(shù)，如表1所示。通過離散傳遞函數(shù)可得到基于CR法位移速度遞推式和MQBS法極點的顯式算法參數(shù)，實現(xiàn)隱式算法的顯式化。

表1 開環(huán)傳遞函數(shù)的系數(shù)

式(11)的G(z)分母項等于0時為特征方程，將式(7)所示的極點代入解特征方程即可得到參數(shù)表達(dá)式如下：

(12)

式中：ξ為系統(tǒng)的阻尼比。

新算法計算結(jié)構(gòu)響應(yīng)包括以下步驟：①輸入結(jié)構(gòu)初始特性，及時間步長Δt；②由式(1)求解初始加速度；③選擇系數(shù)X代入式(12)計算參數(shù)；④通過式(2)、式(3)計算第i時步的位移和速度，將結(jié)果代入式(1)得到該時刻的加速度；⑤令i=i+1，重復(fù)第④步的計算，直至結(jié)束。

2 算法特性分析

2.1 穩(wěn)定性分析

線性系統(tǒng)穩(wěn)定性可通過分析傳遞矩陣的特征值和傳遞函數(shù)的極點得到，二者等價[9]。若傳遞矩陣A的譜半徑滿足ρ(A)≤1，則算法是穩(wěn)定的。新算法譜半徑如式(13)所示，取不同的系數(shù)X，作譜半徑圖。由式(13)可得，對于任意系數(shù)X，新算法譜半徑恒等于1，故新算法對線性系統(tǒng)是無條件穩(wěn)定的。

(13)

對于非線性系統(tǒng)，系統(tǒng)在t+1時步的剛度為kt+1，定義kt+1與初始剛度k0的比值δt+1為剛度變化系數(shù)。若δt+1>1則為剛度硬化，反之δt+1<1為剛度軟化。非線性系統(tǒng)穩(wěn)定性由離散控制理論閉環(huán)傳遞函數(shù)進(jìn)行分析[15]。根據(jù)離散控制理論，閉環(huán)傳遞函數(shù)為

(14)

表2 閉環(huán)傳遞函數(shù)的系數(shù)

非線性系統(tǒng)特征方程為式(14)的分母等于0，即

(15)

圖1 本文算法根軌跡

(16)

2.2 精度分析

一種算法的精度可由周期延長率、振幅衰減率和截斷誤差進(jìn)行分析[5, 11-13]。由文獻(xiàn)[11]，本文算法的截斷誤差計算結(jié)果如下：

(17)

可見，算法具有二階精度，且當(dāng)α= 1/2時，算法最高可滿足三階精度。除截斷誤差，可通過振幅衰減率(AD)和周期延長率(PE)分別反映算法的振幅誤差和周期誤差。將算法極點改寫為如下形式：

(18)

AD和PE可表示為式(19)(20)的形式，取不同的系數(shù)值，作新算法與Newmark系列算法的AD和PE對比圖。

(19)

(20)

式中：AD為振幅衰減率；PE為周期延長率。

2.3 超調(diào)性分析

當(dāng)時間步長較大時，非零初始條件下的單自由度系統(tǒng)可能會出現(xiàn)系統(tǒng)位移速度響應(yīng)異常放大現(xiàn)象，即超調(diào)現(xiàn)象[14-15]。研究算法的超調(diào)性通常分析系統(tǒng)在非零初始位移和速度、無外荷載的情況下，算法位移、速度在第一個時間步長內(nèi)隨Ω變化的趨勢。以x0，v0表示初始位移和速度，第一時步的位移和速度可表示為式(21)，由式(21)可知算法在第一時步不會出現(xiàn)位移速度響應(yīng)的放大現(xiàn)象，即NDEM算法無超調(diào)。

(21)

3 算例

算例1 單自由度系統(tǒng)自由振動

圖5和圖6為不同時間步長的響應(yīng)對比，圖7為100個周期內(nèi)的長期響應(yīng)對比。圖中的位移時程曲線表明，CR法和中心差分法均出現(xiàn)周期延長，且隨著計算時步的增長出現(xiàn)誤差累積現(xiàn)象。隨著系數(shù)X取值的增加，NDEM法的精度逐漸下降，與前文分析結(jié)果一致。從圖6可明顯看出，當(dāng)X=1時，NDEM法的精度最高，優(yōu)于CR法和中心差分法，且CR法計算結(jié)果與X=3的NDEM法計算結(jié)果重合。

算例2 雙自由度系統(tǒng)受恒載激勵

(22)

式中：x的下標(biāo)1，2分別代表點1、點2；系統(tǒng)初始位移x1=x2=0。

由以上條件可得精確解為

(23)

取X=1，用NDEM法、CR法和中心差分法求解此系統(tǒng)的位移，時間步長取Δt=0.1 s，各算法計算結(jié)果與精確解的對比如圖8—10所示。從圖中可以看出，三種顯式算法都與精確解基本吻合。由圖9可知，NDEM法與精確解的吻合程度更高。

算例3 非線性系統(tǒng)

三層框架結(jié)構(gòu)各層質(zhì)量為m1= 104kg，m2=m3= 103kg，初始剛度為k01= 105kN/m，k02=100 kN/m，k03=100 kN/m。該結(jié)構(gòu)為非線性剛度硬化系統(tǒng)，k1=k01(1+10·Δu12)，k2=k02(1+0.1·Δu22)，k3=k03(1+0.1·Δu32)，其中Δu為各層的層間位移。在結(jié)構(gòu)底部施加水平正弦加速度100sin(πt)，將Δt=0.001 s的Newmark顯式法作為精確解與其余算法對比，結(jié)果如圖11、圖12所示。

圖11為Δt=0.01 s時各算法得到的頂層位移與精確解的對比，從圖中可以看出，對于非線性系統(tǒng)，新算法計算結(jié)果與精確解吻合良好，具有較好的精度。圖12為Δt=0.02 s時各算法得到的底層位移與精確解的對比，此時中心差分法出現(xiàn)失穩(wěn)，無法得到正確結(jié)果，而新算法與CR法均能得到正確結(jié)果，表明新算法有較好的穩(wěn)定性。

4 結(jié)論

基于離散控制理論和隱式的修正四次樣條插值算法，設(shè)計一種無條件穩(wěn)定的顯式新算法NDEM。新算法具有二階精度，參數(shù)α=1/2可滿足三階精度，無振幅衰減，無超調(diào)，周期延長率隨時間步長Δt的增加而增大，當(dāng)算法調(diào)節(jié)系數(shù)X=1時周期延長率絕對值最小，算法精度最高。新算法對線性系統(tǒng)和剛度軟化系統(tǒng)是無條件穩(wěn)定，對剛度硬化系統(tǒng)為條件穩(wěn)定，其穩(wěn)定界限隨系數(shù)X的增加而增加。算例分析表明，系數(shù)X可調(diào)整NDEM法的周期延長率進(jìn)而控制算法精度，與理論分析結(jié)果一致，驗證了新算法是有效的。