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

趙 旭，趙建鋒

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

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

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

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

1 算法推導(dǎo)

1.1 運(yùn)動(dòng)微分方程

在第i時(shí)刻，線性單自由度(SDOF)系統(tǒng)的運(yùn)動(dò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為積分時(shí)間步長(zhǎng)，ω為系統(tǒng)固有頻率。

1.2 離散控制理論

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

(10)

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

(11)

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

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

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

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

(12)

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

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

2 算法特性分析

2.1 穩(wěn)定性分析

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

(13)

對(duì)于非線性系統(tǒng)，系統(tǒng)在t+1時(shí)步的剛度為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 精度分析

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

(17)

可見(jiàn)，算法具有二階精度，且當(dāng)α= 1/2時(shí)，算法最高可滿足三階精度。除截?cái)嗾`差，可通過(guò)振幅衰減率(AD)和周期延長(zhǎng)率(PE)分別反映算法的振幅誤差和周期誤差。將算法極點(diǎn)改寫為如下形式：

(18)

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

(19)

(20)

式中：AD為振幅衰減率；PE為周期延長(zhǎng)率。

2.3 超調(diào)性分析

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

(21)

3 算例

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

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

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

(22)

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

由以上條件可得精確解為

(23)

取X=1，用NDEM法、CR法和中心差分法求解此系統(tǒng)的位移，時(shí)間步長(zhǎng)取Δt=0.1 s，各算法計(jì)算結(jié)果與精確解的對(duì)比如圖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顯式法作為精確解與其余算法對(duì)比，結(jié)果如圖11、圖12所示。

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

4 結(jié)論

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