基于修正四次樣條插值法的結構動力學顯式算法

趙 旭，趙建鋒

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

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

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

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

1 算法推導

1.1 運動微分方程

在第i時刻，線性單自由度(SDOF)系統的運動微分方程如式(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為積分時間步長，ω為系統固有頻率。

1.2 離散控制理論

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

(10)

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

(11)

式中：n，d為系數。

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

表1 開環傳遞函數的系數

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

(12)

式中：ξ為系統的阻尼比。

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

2 算法特性分析

2.1 穩定性分析

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

(13)

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

(14)

表2 閉環傳遞函數的系數

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

(15)

圖1 本文算法根軌跡

(16)

2.2 精度分析

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

(17)

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

(18)

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

(19)

(20)

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

2.3 超調性分析

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

(21)

3 算例

算例1 單自由度系統自由振動

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

算例2 雙自由度系統受恒載激勵

(22)

式中：x的下標1，2分別代表點1、點2；系統初始位移x1=x2=0。

由以上條件可得精確解為

(23)

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

算例3 非線性系統

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

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

4 結論

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