高階單步法的非線性誤差傳遞效應分析

欒強利，陳章位，文 祥

(浙江大學 流體動力與機電系統國家重點實驗室，杭州 310027)

擬動力試驗廣泛應用于各類建筑結構的抗震性能研究，試驗過程誤差的大小會影響到試驗結果的可靠性，因此在試驗過程中我們需要對擬動力試驗的誤差進行有效控制。擬動力試驗過程中，對試驗誤差的有效控制一方面可以通過采用高精度的試驗設備降低試驗誤差，另一方面可以通過選擇合適的數值積分方法減少試驗過程中每一步誤差的累積，因此，一個優秀的數值積分方法應該具有較小的誤差傳遞效應，從而有效控制試驗結果的準確性［1－2］。

高階單步法是由王煥定等［3］于1996年提出，具有高精度、無條件穩定的特性，已經應用于非線性結構動力反應的時程分析;崔雪娜等［4］在高階單步法的研究基礎上，對在結構剛度為時間解析函數時高階單步法的適用性進行了研究。本文在前人研究的工作基礎上，進一步研究了系統結構剛度解析表達時系統剛度及算法參數β的變化對高階單步法穩定性及其數值阻尼的影響，另一方面，通過采用兩步累積誤差迭代的方法，研究高階單步法的誤差傳遞效應，進一步驗證了算法的數值阻尼對系統積分累積誤差的影響。

1 高階單步法

非線性動力方程可以表示為:

將其寫成狀態方程的形式為:

其中:

高階單步β法在t=t0+(n+1)Δt時刻由位移、速度組成的狀態向量Zn+1可由下式計算［3］:

考慮到非線性動力方程中參數的時間相關性，則

其中:

將式(6)、(7)代入式(5)并進行化簡，得到:

其中:

對于離散載荷

將式(13)代入式(8)并整理得:

其中:

根據式(14)，第n+1時間步的位移和速度與第n時間步的位移和速度有關。

2 高階單步法的穩定性及其數值阻尼分析

2.1 算法的穩定性分析

研究高階單步法的穩定性及其數值阻尼在非線性系統中的影響因素，忽略系統的物理阻尼，同時考慮系統剛度的變化，參考文獻［5－7］提出的剛度非線性系數模型:

對于非線性系統其結構剛度k是變化的，為了能夠評價系統中結構剛度的變化，模型中引入非線性系統離散模型的結構剛度變化系數μn+1，同時將系統的結構剛度變化解析方程表示為:

進一步得到離散系統的剛度變化系數

當剛度變化系數μn+1＞1時，結構存在剛度硬化，當剛度變化系數μn+1＜1時，結構發生剛度軟化，當μn+1=1時，結構為剛度常量，本文在研究剛度變化系數對算法穩定性及其數值阻尼的影響過程中，剛度變化系數取值0.8，1.0 和 1.2，分別對應系統結構剛度軟化、剛度常量以及剛度硬化情況。

將系統頻率［8］的概念代入式(17)，整理得到

將式(20)代入式(4)、(7)、(9)和(10)得到

其中:

系統放大矩陣

對應特征方程為

Cij(i，j=1，2)為放大矩陣對應的元素，相應的矩陣譜半徑 ρ(C)=max{|λi|，i=1，2}，是關于 a、δn、δn+1、Δt、β、ω0的函數。通過數值分析驗證，式(25)的判別式

判別式的結果與剛度變化系數μn+1及β等的取值有關，但總是小于0的，因此放大矩陣的譜半徑可表示為:

根據文獻［3］得到的結論，算法無條件穩定的條件是β≥0.5，因此首先考慮在 β =0.5時，剛度變化系數對系統穩定性的影響，三種剛度變化模型對應的各自譜半徑隨結構頻率ω0的變化規律如圖1所示，當μn+1=0.8，δn+1=0.8n+1時，結構剛度軟化，在 n=2 的情況下進行分析，放大矩陣的譜半徑始終大于1，使得計算結果發散，即積分算法不穩定;當 μn+1=1.0，δn+1=1.0時，結構剛度為常量，其譜半徑恒等于1，算法臨界穩定;當 μn+1=1.2，δn+1=1.2n+1，時，結構剛度硬化，對應譜半徑小于1，算法無條件穩定。進一步研究發現在μn+1=0.8，n=2剛度軟化的情況下，調節β值，發現當β取值為0.518 6時，算法臨界穩定，如圖2所示，此時對應的其余兩種剛度變化情況，算法均能實現無條件穩定，因此，對于剛度軟化情況下導致的算法不穩定狀態，可以通過調節β的取值實現算法的穩定。在本文研究的預設條件下，對于系統剛度軟化情況，臨界穩定狀態 β的取值與 n和 μn+1取值的關系如表1，表2所示。

表1列出在μn+1=0.8時，為保證算法的無條件穩定，β的最小取值與對應的n取值的關系，表中n的取值范圍是在0～25之間，臨界穩定狀態下β的變化范圍為0.518 5～0.53，且β的取值隨著n值的增大逐漸增大，當n過大時，結構剛度接近于零，這種情況在實際試驗過程中是不存在的，因此對n值在討論范圍之外的情況未作進一步討論。表2表示在n=2時，為保證算法的無條件穩定，β的最小取值與對應的μ值之間的關系，當μ在0.4～1.0之間變化時，β的取值是隨著μ取值的增大逐漸減小的，當μ的取值為1.0時，系統的結構剛度為常量，β的算法穩定臨界取值為0.5，這也驗證了文獻［3］中關于β在結構剛度為常量情況下，其對于算法穩定性的取值條件，即β≥0.5，在μ的取值小于0.4時，系統結構剛度過小，因此未作進一步討論。

圖1 β=0.5時不同剛度變化系數下ρ－ω0關系Fig.1 Numerical relationship between ρandω0 under different stiffness variation coefficients with β=0.5

圖2 β=0.518 6時不同剛度變化系數下ρ－ω0關系Fig.2 Numerical relationship between ρandω0 under different stiffness variation coefficients with β=0.518 6

圖3 μ=0.8，n=20時不同 β取值下ρ－ω0關系Fig.3 Numerical realationship betweenρandω0 under differentβ with μ =0.8 and n=20

對于表1中在μ=0.8，n=20時，不同的β取值對譜半徑的影響如圖3所示，當ω0的取值大于180時，隨著β取值的增大，譜半徑ρ(C)發生迅速衰減，實驗證明，同樣的變化規律發生在n(n≤25)取不同數值的情況下，即β值越大，譜半徑ρ(C)衰減越快。

表1 μ=0.8時臨界β值與n的對應關系Tab.1 Numerical relationship between the thresholdβand n withμ=0.8

表2 n=2時臨界β值與μ的對應關系Tab.2 Numerical relationship between the thresholdβandμwith n=2

2.2 算法的數值阻尼分析

擬動力試驗過程中的累積誤差是隨著試驗結構頻率的增加而不斷增長的，即高頻響應對試驗的累積誤差影響較大［9－11］，因此帶有數值耗散功能即數值阻尼的積分方法可以用來有效抑制或消除算法中隨結構高頻響應增加的累積誤差即降低算法的誤差傳遞效應。系統放大矩陣的特征值可以表示為［8］:

其中:

自由振動條件下，結構的位移響應可以表示為:

將式(28)代入上式，并令:

則結構的位移響應為:

因此，位移響應的大小主要取決于式(33)中右端的第一項因子(第二項為有限值)，進一步化簡式(33)得到:

其中:

反映了算法對結構位移響應的阻尼作用。根據表2中n與β的關系，考慮n=2，β=0.52時，數值阻尼比隨剛度變化系數μ的變化情況如圖4所示，可見在給定條件下，當結構為剛度常量或硬化情況時，算法中引入了數值阻尼，而在結構剛度軟化情況時，算法中并沒有引入足夠的算法阻尼。因此，為了使高階單步法在結構剛度軟化過程中仍具有一定的算法阻尼，需要進一步考慮關于算法數值阻尼的影響因素，為此，單獨討論結構剛度軟化情況時，β取值的變化對算法數值阻尼的影響，在n=2，μ=0.8時，β取值變化對數值阻尼比的影響如圖5所示，隨著β取值的增大，系統的阻尼比有逐漸增大的趨勢，因此對于結構剛度軟化情況，通過適當增大β值，可以有效調整算法的數值阻尼。

圖4 n=2，β=0.52時剛度變化對數值阻尼比的影響Fig.4 Numerical damping ratio under different stiffness variation coefficients with n=2 and β=0.52

圖5 n=2，μ=0.8時β變化對數值阻尼比的影響Fig.5 Numerical damping ratio under differentβ with n=2 and μ =0.8

3 系統阻尼對誤差傳遞效應的影響

擬動力試驗過程中，當我們使用高階單步法做位移控制時，位移步之間存在誤差的傳遞效應，一方面是由于試驗過程中傳感器的靈敏度為有限值，另一方面是由于計算機的計算精度為有限位數，都將導致試驗的計算位移與反饋位移之間存在誤差，擬動力試驗過程需要對這些本身帶有誤差的數據進行不斷的迭代運算，因此導致試驗誤差的累積，形成累積誤差，對高階單步法誤差傳遞效應的研究，參考文獻［5］中對算法誤差傳遞效應的分析方法，引入下列變量:

dn:第n步精確數值計算位移;

den:第n步精確位移，包括之前步的位移誤差;

dan:第n步實際位移，包括當前步與之前步的位移誤差;

edn:第n步產生的位移誤差。

它們之間存在如下的關系:

根據式(36)，對式(14)進行重寫，則

令

得到算法的累積誤差

其中:

由于算法的累積誤差與其積分步數有關，因此需要以確定的積分步數對算法的誤差傳遞效應進行研究，參考文獻［9－11］，采用兩步累積誤差迭代的方法對誤差傳遞效應進行研究，并進一步化簡式(39)，得到

其中:Dn表示對當前第n+1步之前第n步誤差的誤差傳遞作用，Dn－1表示對當前第n+1步之前第n－1步誤差的誤差傳遞作用，Dn－2表示對當前第n+1步之前第n－2步誤差的誤差傳遞作用。

根據式(40)，當前第n+1步的位移累積誤差為

為進一步研究算法對位移時程的誤差傳遞效應，研究當 n=2，β =0.52 時，Dn－2、Dn－1與 Dn對其各自積分步中位移誤差項的作用，在不同的剛度變化系數下，Dn－2、Dn－1與 Dn隨系統頻率的變化如圖 6 所示。通過對圖6中(a)、(b)和(c)的比較，發現當前積分時間步的累積誤差對過去時間步的誤差有明顯的抑制作用，但抑制作用的強弱有明顯的不同，Dn－2衰減最快，說明的誤差抑制作用顯著，Dn－1次之，Dn最差，說明過去的積分時間步距離當前步的時間越久，其對應的誤差傳遞效應越弱，即當前積分時間步的累積誤差效果主要取決于過去有限時間步的步數。同時，對比在不同剛度系數條件下算法的誤差傳遞效應，發現在剛度軟化過程中，其對應的誤差傳遞效應比在剛度常量及硬化過程中的效應更強，而在剛度硬化情況下，其對應的誤差傳遞效應較小。

前面已經證明，在擬動力試驗過程中，我們可以通過調節β的取值，控制算法的數值阻尼，因此，通過改變β值可以實現對試驗過程誤差傳遞效應的控制，圖7為考慮結構剛度變化系數為0.8的情況下，取不同的β值對 Dn－2、Dn－1和 Dn的影響，對比發現，在確定的結構頻率點處，當β(大于0.5)的取值過小時，系統的誤差傳遞效應較為明顯，尤其是對Dn的影響較大，隨著β取值的逐漸增大，算法的誤差傳遞效應逐漸變小，同時，對于確定的 β 值，Dn－2、Dn－1和 Dn對應于不同結構頻率表現出不同的誤差傳遞效應，在系統頻率較高時，由之前的討論可知，系統的阻尼較大，因此，系統的誤差傳遞效應明顯，Dn－2、Dn－1能夠較快地趨于零，Dn趨于某一較小的確定值(由β值決定)。上述結論是在剛度變化系數為0.8的情況下得到的，經過大量的數值計算，可以證明這個結論同樣適用于剛度常量及硬化過程。

4 誤差傳遞效應仿真實驗

為進一步說明高階單步法的誤差傳遞效應對擬動力試驗過程的影響，對單步積分算法進行數值實驗仿真。數值實驗采用的地震波為1940 EI Centro地震波，峰值加速度為2 g，為便于結果分析，取單自由度非線性系統，質量為40 000 kg，剛度108kN/m，忽略系統的物理阻尼，同時考慮引入誤差對單步積分結果的影響，研究當b取0.9、1.0及1.1時，結構的位移時程響應和累積誤差隨系統剛度的變化情況。三種取值條件下系統對應的結構剛度變化系數分別為0.997 90、1.0和1.001 91，相應的分析結果如圖8～10所示，其中(a)圖為位移時程曲線(理論曲線與實際曲線)，(b)圖為引入的隨機誤差，其幅值為0.12 mm，(c)圖為累積誤差曲線。對比圖8～10中的(a)圖即位移響應時程曲線，發現三種(剛度軟化、剛度常量和剛度硬化)情況下，引入誤差后結構的位移時程曲線(實際曲線)仍與理論曲線有較好的重疊，說明三種情況下對應的系統均有較強抗干擾能力，通過對三種情況下對應的系統累積誤差分析，可以進一步得到上述的結論。同時，分析發現當結構剛度變化系數取值(0.997 90)小于1發生剛度軟化情況時，對比其它兩種情況，在系統的結構剛度較小情況下，系統的累積誤差較大，并有逐漸變大的趨勢，說明系統的誤差傳遞效應較強，而在系統結構剛度常量及硬化過程中，系統的累積誤差均得到有效的控制，誤差被控制在一個穩定的范圍內，并沒有隨積分步數的增多持續增大，誤差傳遞效應較小。三種情況下對應的系統相位角Ω=ωΔt的變化如表3所示，在結構軟化過程中，相位角的變化范圍為1.0～0.418，逐漸變小，由之前討論的關于剛度軟化對系統數值阻尼及累積誤差的影響可知，Ω較小(系統頻率較小)時數值阻尼較小，誤差傳遞效應較強，導致系統累積誤差逐漸增大;另一方面，系統頻率的變小使得系統相較于線性系統(圖9)振動周期變長，有限時間里歷經的循環次數變少(圖8)。在系統結構剛度硬化過程中，系統相位角Ω的變化范圍為1.0～3.14，逐漸增大，數值阻尼增大，誤差傳遞效應變弱;同時，系統頻率的增大導致系統相較于線性系統(圖9)周期的縮短，因而在有限時間里系統歷經的循環次數增多(圖10)，循環次數的增多將導致系統累積誤差的增大，在系統結構硬化過程中，由于大阻尼的存在能夠很好的抑制系統誤差的累積，使系統的誤差累積效應較小(圖10(c))。

圖6 β =0.52時不同剛度變化系數對 Dn－2、Dn－1和 Dn的影響Fig.6 Dn－2，Dn－1 and Dn under different stiffness variation coefficients with β =0.52

圖 7 μ =0.8 時 β 變化對 Dn－2、Dn－1和 Dn的影響Fig.7 Dn－2，Dn－1 and Dn under differentβ with stiffness variation coefficientμ =0.8

圖8 剛度軟化時擬動力試驗仿真Fig.8 Numerical simulation of pseudodynamic testing of consistent softening system

圖9 剛度常量時擬動力試驗仿真Fig.9 Numerical simulation of pseudodynamic testing of linear elastic system

圖10 剛度硬化時擬動力試驗仿真Fig.10 Numerical simulation of pseudodynamic testing of consistent hardening system

表3 i=0和600時對應的Ω1值Tab.3 Values ofΩ1 for i=0 and 600

表4 i=0和600時對應的ωi值Tab.4 Values ofωi for i=0 and 600

圖11 剛度軟化時β對系統累積誤差的影響Fig.11 Influence on the cumulative error under differentβs for the consistent softening system

在系統結構剛度軟化過程中，系統的頻率從50 rad/s減小到20.9 rad/s(表4)，由圖5，圖7的理論分析結果可知，在系統的結構頻率較小時，算法的數值阻尼較小，結果導致系統有較強的誤差累積效應。為進一步研究系統結構剛度軟化過程中β變化對系統累積誤差的作用，討論 β 取值分別為0.4，0.5，0.6 和0.7 情況下，系統的累積誤差，結果如圖11所示，在β取值0.4時，隨著系統剛度的減小，系統累積誤差有逐漸增大的趨勢，導致系統失穩，在β取值0.5情況下，由于算法數值阻尼的增大(圖5)，系統的累積誤差得到有效控制，進一步增大β，系統的累積誤差進一步減小，對比圖11中(b)、(c)、(d)發現隨著β取值的增大，系統的累積誤差有減小的趨勢，繼續增大β(大于0.6)，算法對系統累積誤差的影響已不明顯，系統的累積誤差趨于穩定的變化范圍。

5 結論

本文對高階單步法在非線性系統中的應用進行了研究，討論了剛度變化系數μ和參數β對系統的穩定性及其數值阻尼的影響，證明在系統軟化過程中β取值過小將導致系統不穩定，而通過適當增大β值可以滿足系統在結構剛度軟化過程中的穩定性;同時，在系統剛度軟化過程中系統的數值阻尼較小，這不利于系統對高階模態響應的抑制作用，而通過增大β值能夠使系統即使在剛度軟化的過程中仍具有較好的阻尼特性。系統的誤差傳遞效應對系統結構剛度變化系數μ和參數β的取值比較敏感，在系統結構剛度軟化過程中，系統的誤差傳遞效應比較強烈，導致系統的累積誤差增大，通過適當調節β的取值可以有效減少系統的累積誤差，削弱系統的誤差傳遞效應。

［1］Shing P B，Mahin S A.Cumulative experimental errors in pseudodynamic tests［J］.Earthquake Engineering and Structural Dynamics，1987，15:409 －424.

［2］Shing P B，Mahin SA.Elimination of spurious higher-mode response in pseudodynamic tests ［J］. Earthquake Engineering and Structural Dynamics，1987，15:425－445.

［3］王煥定，張永山，王 偉，等.非線性結構時程分析的高階單步法［J］.地震工程與工程振動，1996，16(3):48－54.WANG Huan-ding，ZHANG Yong-shan，WANG Wei，et al.The high order one step－βmethod for seismic response analysis of non-linear structures［J］.Earthquake Engineering and Engineering Vibration，1996，16(3):48－54.

［4］崔雪娜，王煥定.剛度解析表達時高階單步法的研究［J］.地震工程與工程振動，2006，26(4):63－67.CUI Xue-na，WANG Huan-ding.Research on a high order single step method for analytic expression of stiffness［J］.Earthquake Engineering and Engineering Vibration，2006，26(4):63－67.

［5］Chang S Y.Nonlinear error propagation analysis for explicit pseudodynamic algorithm ［J］. Journal of Earthquake Mechanics，2003，129(8):841 －850.

［6］Chang SY.Performance of HHT－αmethod for the solution of nonlinear systems［J］.International Journal of Structural Stability and Dynamics，2008，8(2):321－337.

［7］Chang S Y.Numerical characteristics of explicit method for nonlinear systems［J］.Journal of Earthquake Engineering，2007，11(5):694－711.

［8］邱法維，錢稼茹，陳志鵬.結構抗震實驗方法［M］.北京:科學出版社，2000.

［9］Chang SY.Nonlinear performance of explicit pseudodynamic algorithm ［J］.Journal of Earthquake Engineering，2010，14(2):211－230.

［10］Chang S Y.Error propagation of HHT － α method for pseudodynamic tests ［J］. Journal of Earthquake Engineering，2005，9(2):223－246.

［11］Chang SY，Sung Y C.Analytical exploration ofγ－method explicit method for pseudodynamic testing of nonlinear system［J］.Earthquake Engineering and Engineering Vibration，2005，4(1):117－127.