基于計(jì)算擾動(dòng)修正的Wilson-θ法穩(wěn)定性和精度分析

胡林烽，雷家艷,2*，吳嘉偉，高 婧

(1.廈門大學(xué)建筑與土木工程學(xué)院，福建 廈門 361005；2.廈門市交通基礎(chǔ)設(shè)施智能管養(yǎng)工程技術(shù)研究中心，福建 廈門 361005)

直接時(shí)間積分法是結(jié)構(gòu)動(dòng)力分析中最常用的積分算法之一.該算法將結(jié)構(gòu)的激勵(lì)和響應(yīng)離散為一系列小的時(shí)間間隔，并假設(shè)加速度、速度和位移在每個(gè)時(shí)間步長(zhǎng)內(nèi)的變化滿足一定的關(guān)系，然后在離散時(shí)間域逐步求解運(yùn)動(dòng)方程.直接時(shí)間積分法既可應(yīng)用于線性系統(tǒng)，又可應(yīng)用于非線性系統(tǒng)，并且具有良好的計(jì)算精度，因此，它在近20年來(lái)得到了廣泛的應(yīng)用[1-4].

傳統(tǒng)的Wilson-θ法是在線性加速度假設(shè)的基礎(chǔ)上發(fā)展起來(lái)的一種隱式直接時(shí)間積分方法.當(dāng)θ≥1.37 時(shí)，Wilson-θ法具有無(wú)條件穩(wěn)定性，因此該方法的時(shí)間步長(zhǎng)不受結(jié)構(gòu)自振周期的限制.但是隨著對(duì)Wilson-θ法研究的逐漸深入，發(fā)現(xiàn)這種方法存在不可避免的缺點(diǎn)：一個(gè)時(shí)間步長(zhǎng)內(nèi)Wilson-θ法求得的下個(gè)時(shí)刻的位移、速度、加速度并不滿足該時(shí)刻的動(dòng)力平衡方程，從而會(huì)導(dǎo)致較為嚴(yán)重的振幅衰減和周期延長(zhǎng).因此，Wilson-θ法的應(yīng)用受到了一定的限制[5-7].

為了解決Wilson-θ法的缺陷，許多學(xué)者進(jìn)行過(guò)深入研究，對(duì)原有的Wilson-θ法進(jìn)行了一定改進(jìn)[8-10].Clough等[7]將時(shí)間步長(zhǎng)終點(diǎn)的位移和速度代入動(dòng)力平衡方程重新求解加速度進(jìn)行修正，稱為加速度修正的Wilson-θ法.但是這種方法雖然能夠滿足動(dòng)力平衡方程，卻也會(huì)使求得的位移、速度和加速度無(wú)法滿足運(yùn)動(dòng)約束條件[11].黃慶豐等[12]在加速度修正方法的基礎(chǔ)上定義不平衡加速度為原加速度減去修正加速度，并定義計(jì)算擾動(dòng)為不平衡加速度產(chǎn)生的附加作用，利用計(jì)算擾動(dòng)構(gòu)造一個(gè)新的動(dòng)力增量平衡方程并繼續(xù)使用Wilson-θ法進(jìn)行求解，再將結(jié)果與原方程相加進(jìn)行修正.方德平[13]提出了一種修正荷載的Wilson-θ法，在t時(shí)刻增加一個(gè)不平衡荷載以滿足t時(shí)刻的動(dòng)力平衡方程，稱之為修正荷載，然后在t+Δt時(shí)刻的計(jì)算中加入該荷載并乘以修正系數(shù)θ′以抵消修正荷載的效應(yīng).黃慶豐[14]還通過(guò)定義不平衡加速度，并通過(guò)補(bǔ)充能量方程導(dǎo)出θΔt時(shí)間段內(nèi)系統(tǒng)計(jì)算位移修正量與不平衡計(jì)算加速度之間的關(guān)系，以此確定加速度和速度修正量.除此以外，還有其他一些改進(jìn)的Wilson-θ法通過(guò)應(yīng)用高階動(dòng)力學(xué)方程來(lái)消除傳統(tǒng)方法的缺點(diǎn)[15].但是，現(xiàn)有的改進(jìn)方法仍存在一些不足.

本文基于計(jì)算擾動(dòng)[12]的概念，將計(jì)算擾動(dòng)與Newmark-β法相結(jié)合而非繼續(xù)采用Wilson-θ法求解，提出了一種新的Wilson-θ改進(jìn)方法，將其稱為Wilson-Newmark法，使得這種方法能夠同時(shí)滿足動(dòng)力平衡方程和運(yùn)動(dòng)約束條件，在理論上消除了原始Wilson-θ法的計(jì)算誤差，并結(jié)合算例針對(duì)修正方法的穩(wěn)定性、精度以及計(jì)算效率進(jìn)行了分析，綜合評(píng)價(jià)了改進(jìn)方法的優(yōu)缺點(diǎn).

1 Wilson-Newmark法

考慮多自由度體系的動(dòng)力系統(tǒng)，其動(dòng)力平衡方程如下所示：

(1)

傳統(tǒng)的Wilson-θ法假定加速度在[t,t+θΔt]時(shí)間段內(nèi)為線性變化，并結(jié)合位移增量方程，可以得到加速度、速度和位移增量：

(2)

(3)

(4)

其中τ為時(shí)間周期θΔt.為了克服傳統(tǒng)的Wilson-θ法在t+Δt時(shí)刻不滿足動(dòng)力平衡方程的缺點(diǎn)，引入計(jì)算擾動(dòng)[12]的概念，將Wilson-θ法求得的速度和位移增量重新代入動(dòng)力平衡方程中，可以得到修正的加速度：

(5)

定義計(jì)算擾動(dòng)為Wilson-θ法計(jì)算的加速度減去修正加速度并乘以質(zhì)量矩陣：

(6)

將式(5)代入式(6)，可以重新得到一個(gè)完整的動(dòng)力方程：

K[Δz(t+Δt)+z(t)]=F(t+Δt)+ΔFv(t).

(7)

因此，為了消除計(jì)算擾動(dòng)的影響，可以構(gòu)建一個(gè)新的t+Δt時(shí)間段內(nèi)的修正方程：

-ΔFv(t),

(8)

在此基礎(chǔ)上，應(yīng)用Newmark-β法進(jìn)行動(dòng)力增量平衡方程的求解，與Wilson-θ方法不同的是, Newmark-β方法假設(shè)加速度在時(shí)間[t,t+θΔt]是一個(gè)常數(shù).為了分配計(jì)算擾動(dòng)，假定在當(dāng)前時(shí)間段內(nèi)計(jì)算擾動(dòng)的修正方程的位移、速度及加速度的初值為零，則該修正方程的等效荷載即為-ΔFv(t).計(jì)算Newmark-β法的等效剛度為：

(9)

其中γ和β為Newmark-β法的計(jì)算參數(shù).利用動(dòng)力增量平衡方程，將等效荷載除以等效剛度可以得到修正位移增量為：

(10)

根據(jù)平均常加速度假設(shè)，修正速度增量和修正加速度增量分別為：

(11)

(12)

利用求得的各增量可以重新構(gòu)建動(dòng)力增量平衡方程，將修正方程代入式(11)中可得：

K[Δzv(t+Δt)+Δz(t+Δt)+z(t)]=

F(t+Δt).

(13)

因此，時(shí)間步長(zhǎng)終點(diǎn)t+Δt的位移、速度和加速度可以寫成如下形式：

(14)

根據(jù)式(8)～(14)可以對(duì)Wilson-θ法動(dòng)態(tài)平衡方程中存在的誤差進(jìn)行修正，并對(duì)其穩(wěn)定性和精度進(jìn)行分析.

2 穩(wěn)定性分析

直接時(shí)間積分算法的穩(wěn)定性是指隨著積分時(shí)間步長(zhǎng)的增大，計(jì)算中求得的數(shù)值解不會(huì)產(chǎn)生過(guò)大誤差的性質(zhì).當(dāng)直接時(shí)間積分算法的穩(wěn)定性條件不滿足時(shí)，求得的計(jì)算結(jié)果將產(chǎn)生發(fā)散，使得計(jì)算誤差呈指數(shù)增長(zhǎng).因此，穩(wěn)定性是積分計(jì)算過(guò)程中需要考慮的一個(gè)重要因素[16].以圖1所示的典型單自由度(SDOF)系統(tǒng)為例，其動(dòng)力平衡方程為：

圖1 SDOF結(jié)構(gòu)系統(tǒng)Fig.1 The SDOF structural system

(15)

式中，zt、zt、zt分別為結(jié)構(gòu)的位移、速度和加速度，ω為固有頻率，ζ為阻尼比，ft為激振力.

針對(duì)該系統(tǒng)，可將遞推方程形式表示為[13]：

(16)

ρ(A)=max|λi|, (i=1,2,3),

(17)

其中λi表示傳遞矩陣的第i個(gè)特征值.結(jié)合式(15)及式(16)可知，逐步積分法的穩(wěn)定性取決于固有頻率ω、阻尼比ζ和時(shí)間步長(zhǎng)Δt.此外，由于頻率ω可以改寫為2π /T，因此穩(wěn)定性可以認(rèn)為和Δt/T有關(guān).傳遞矩陣的推導(dǎo)過(guò)程如下.

首先，根據(jù)Wilson-θ法的積分公式，式(15)可以改寫為

(18)

(19)

(20)

上式中b0、b1、b2均為Wilson-θ法的計(jì)算系數(shù)，b0=6/(θΔt)2，b1=3/(θΔt)，b2=6/(θΔt)，fn為等效荷載中的荷載分量，可以看出該分量在計(jì)算過(guò)程中并不會(huì)影響到傳遞矩陣，因此在后續(xù)推導(dǎo)中將其略去.將等效荷載除以等效剛度，可以得到t+θΔt時(shí)刻的位移：

(21)

根據(jù)Wilson-θ法的基本原理，由t+τ時(shí)刻的加速度值即可求得t+Δt時(shí)刻的位移、速度和加速度：

(22)

其中b3=θΔt/2，b4=6/(θ3Δt2)，b5=-6/(θ2Δt)，b6=1-3/θ，b7=Δt/2，b8=Δt2/6.其他各參數(shù)由下式給出：

D31=b4(E1-1)，D32=b4E2+b5，D31=b4E3+b6，

D21=b7D31，D22=b7D32+1，D23=b7(D33+1)，

D11=b8D31+1，D12=b8D32+Δt，

D13=b8(D33+2).

因此，傳統(tǒng)Wilson-θ法的狀態(tài)傳遞矩陣可以寫成如下形式：

(23)

修正方法需要在Wilson-θ法的基礎(chǔ)上增設(shè)修正方程.修正方程的形式同式(15)，其外荷載可以根據(jù)式(8)得到：

Δt)]=-[(ω2·D11+2ζω·D21+D31)z(t)+

(24)

由于修正方程的初始位移、速度、加速度均為零，其外荷載即為等效荷載.接下來(lái)引入Newmark-β的計(jì)算參數(shù)a0和a1計(jì)算其等效剛度為：

(25)

其中，a0=1/(βΔt2)，a1=γ/(βΔt).由此可以求得Δt時(shí)間段內(nèi)的修正位移、修正速度及修正加速度增量為：

(26)

(27)

(28)

因此，根據(jù)增量可以得到Wilson-Newmark法中的狀態(tài)傳遞矩陣：

(29)

由于該狀態(tài)傳遞矩陣的譜半徑ρ(A)與參數(shù)θ有關(guān)，為了估計(jì)θ的取值對(duì)于該方法無(wú)條件穩(wěn)定性的影響，有必要計(jì)算譜半徑的變化與θ取值之間的函數(shù)關(guān)系.在計(jì)算過(guò)程中，取阻尼比為0，且時(shí)間步長(zhǎng)趨近于無(wú)窮大，以消除其他參數(shù)的影響.譜半徑隨θ的變化如圖2所示.可以看出，在1≤θ≤104的范圍內(nèi)，本文提出的Wilson-Newmark方法ρ(A)始終小于1，可以認(rèn)為該改進(jìn)方法在該范圍內(nèi)具有無(wú)條件穩(wěn)定性，相對(duì)于傳統(tǒng)Wilson-θ法拓寬了其θ的取值范圍.

圖2 Wilson-Newmark法譜半徑隨θ變化曲線Fig.2 The variation of the spectral radius as a function of θ of Wilson-Newmark method

3 精度分析

3.1 振幅衰減和周期延長(zhǎng)

由于直接時(shí)間積分法只是一種近似方法，在求解過(guò)程中必然會(huì)出現(xiàn)計(jì)算誤差，而在隱式積分中這種誤差通常表現(xiàn)為振幅衰減和周期延長(zhǎng).振幅衰減又稱數(shù)值耗散，是指在迭代過(guò)程中波動(dòng)幅度逐漸減小的一種現(xiàn)象.除此以外，在迭代過(guò)程的每個(gè)時(shí)間子步都會(huì)產(chǎn)生一定的相移，這被稱為周期延長(zhǎng)[15].無(wú)阻尼的單自由度系統(tǒng)，其運(yùn)動(dòng)方程為：

(30)

假定該運(yùn)動(dòng)方程的初始位移為1，初始速度為0，則振幅衰減率(RAD)和周期延長(zhǎng)率(RPE)可以定義為[17-18]：

(31)

圖3為不同時(shí)間步長(zhǎng)下Wilson-θ法(θ=1.4)、Newmark-β法(β=1/4,γ=1/2)和Wilson-Newmark法在SDOF系統(tǒng)中的曲線變化趨勢(shì).可以看出隨著dt/T逐漸增加，除了Newmark-β法在理論上沒(méi)有振幅衰減以外，相對(duì)于精確解，3種方法都存在一定的振幅衰減和周期延長(zhǎng).但在相同的dt/T條件下，Wilson-Newmark法的周期延長(zhǎng)率明顯低于其他兩種方法.此外，與Wilson-θ法相比，Wilson-Newmark法在振幅衰減方面也具有明顯優(yōu)勢(shì)，完全達(dá)到了提高精度的效果.精確解采用模態(tài)疊加法和Duhamel積分進(jìn)行計(jì)算，各算法振幅衰減率和周期延長(zhǎng)率的數(shù)值比較如表1所示.

圖3 不同積分方法的精度比較Fig.3 Precision comparison of different direct time integration methods

表1 不同直接時(shí)間積分方法的RAD和RPE比較

3.2 超調(diào)分析

超調(diào)是指無(wú)條件穩(wěn)定的直接時(shí)間積分算法的計(jì)算結(jié)果在步長(zhǎng)較大的情況下，在最初的幾個(gè)子步會(huì)被放大的現(xiàn)象[19-20].研究表明，現(xiàn)有的直接時(shí)間積分方法中，除Newmark-β法外，幾乎所有方法都存在一定的超調(diào)現(xiàn)象[21].算法的超調(diào)通常是通過(guò)檢查第一步計(jì)算的結(jié)果來(lái)確定的.采用3.1節(jié)中的無(wú)阻尼SDOF系統(tǒng)(初始位移為1，初始速度為0)分析超調(diào)問(wèn)題，該系統(tǒng)的自振圓頻率取為1.是否存在超調(diào)現(xiàn)象可以通過(guò)觀察算法在不同Δt條件下的計(jì)算結(jié)果，當(dāng)步長(zhǎng)增大到一定程度時(shí)，第一步計(jì)算的結(jié)果將會(huì)被放大.圖4給出了不同積分方法第一步計(jì)算位移隨步長(zhǎng)的變化曲線.可以看出，Wilson-θ法的第一步計(jì)算位移隨著步長(zhǎng)增大存在明顯的超調(diào)現(xiàn)象，而Wilson-Newmark法和Newmark-β法則無(wú)此缺陷.

圖4 不同積分方法第一步計(jì)算位移隨步長(zhǎng)的變化曲線Fig.4 The relationship between the first step displacement and step size calculated by different direct time integration methods

4 算例分析

為了對(duì)計(jì)算方法的精度以及效率進(jìn)行分析，本節(jié)采取實(shí)例計(jì)算的方式來(lái)對(duì)比分析不同情況下Wilson-Newmark法、Wilson-θ法(θ=1.4)以及Newmark-β法(β=1/4,γ=1/2)的計(jì)算結(jié)果，并定義了新的評(píng)價(jià)參數(shù)，即時(shí)間增加率和誤差縮減率.

4.1 計(jì)算效率

計(jì)算效率是評(píng)判動(dòng)力分析方法的標(biāo)準(zhǔn)之一，在本文中，通過(guò)記錄程序運(yùn)行的迭代過(guò)程所花費(fèi)的時(shí)間，并通過(guò)時(shí)間增加率計(jì)算其與對(duì)比方法的相對(duì)值來(lái)表示該方法的計(jì)算效率，時(shí)間增加率的表達(dá)式為：

tr=(tk-t0)/t0×100%,

(32)

其中，t0是傳統(tǒng)Wilson-θ法以及Newmark-β法運(yùn)行算例過(guò)程所花費(fèi)的時(shí)間，tk是修正的Wilson-Newmark法運(yùn)行算例過(guò)程所花費(fèi)的時(shí)間.

4.2 計(jì)算精度

計(jì)算精度是評(píng)判動(dòng)力分析方法的另一個(gè)標(biāo)準(zhǔn).在修正方法的計(jì)算過(guò)程中，記錄迭代過(guò)程中得到的每一步的位移計(jì)算值以及位移差值，并進(jìn)行累加計(jì)算.計(jì)算誤差的表達(dá)式定義為[11]：

(33)

其中，zc(i)和zt(i)分別為計(jì)算值和理論解，n為計(jì)算步數(shù).將各種情況下計(jì)算得到的計(jì)算誤差進(jìn)行統(tǒng)計(jì)，并引入誤差減小率來(lái)說(shuō)明與傳統(tǒng)Wilson-θ法的相對(duì)值，誤差縮減率的表達(dá)式為：

er=((e0-ek)/e0)×100%,

(34)

其中，e0是傳統(tǒng)Wilson-θ法以及Newmark-β法運(yùn)行算例遞推過(guò)程產(chǎn)生的計(jì)算誤差，ek是經(jīng)過(guò)修正的Wilson-Newmark法運(yùn)行算例遞推過(guò)程產(chǎn)生的計(jì)算誤差.

4.3 數(shù)值算例

算例采用多自由度簡(jiǎn)諧荷載體系進(jìn)行計(jì)算分析，如圖5所示.其計(jì)算參數(shù)如下：結(jié)構(gòu)質(zhì)量為m1=0.2 kg，m2=0.15 kg，m3=0.1 kg，剛度為k1=1 800 N/m，k2=1 200 N/m，k3=600 N/m，簡(jiǎn)諧荷載為F1=F2=F3=5sin(20t) N，時(shí)間步長(zhǎng)為Δt=0.01 s，總時(shí)長(zhǎng)為Tn=1.5 s.運(yùn)動(dòng)方程的初始位移、速度和加速度均為0.算例中分別考慮了有阻尼(ζ=0.006 4)和無(wú)阻尼兩種情況下的系統(tǒng)振動(dòng).

圖5 三自由度結(jié)構(gòu)系統(tǒng)Fig.5 The 3-DOF structural system

基于模態(tài)疊加法和Duhamel積分，給出了該問(wèn)題的解析解.以第一個(gè)自由度為例，計(jì)算結(jié)果對(duì)比如圖6所示.圖6(a)畫出了無(wú)阻尼情況下各方法的計(jì)算結(jié)果，可以清楚地看到，在無(wú)阻尼條件下，采用Wilson-Newmark法的振幅衰減比Wilson-θ法和Newmark-β法都要小，計(jì)算結(jié)果能夠很好地貼合精確解.此外，對(duì)于無(wú)阻尼作用的系統(tǒng)，Wilson-θ法很難反映曲線的微小變化趨勢(shì).圖6(b)為有阻尼情況下3種方法的計(jì)算結(jié)果，可以看出無(wú)論是否存在阻尼，Wilson-Newmark法得到的響應(yīng)始終比其他兩種傳統(tǒng)方法更接近精確解.表2列出了各個(gè)方法分別在有阻尼和無(wú)阻尼情況下計(jì)算精度和計(jì)算效率的求解結(jié)果.其中，時(shí)間增加率和誤差縮減率是將Wilson-Newmark法與原始的Wilson-θ法和Newmark-β法進(jìn)行對(duì)比得到的，因此表格中Wilson-Newmark法這兩項(xiàng)數(shù)據(jù)為0.通過(guò)對(duì)比各項(xiàng)數(shù)據(jù)可以發(fā)現(xiàn)，3種方法在迭代過(guò)程中都會(huì)存在一定的誤差，但由表2對(duì)3種方法的精度和效率分析可以看出Wilson-Newmark法明顯具有最好的收斂性和精度.

圖6 各方法多自由度簡(jiǎn)諧振動(dòng)系統(tǒng)計(jì)算結(jié)果比較Fig.6 Results comparison of multiple degrees of simple harmonic vibration calculated by different methods

表2 多自由度系統(tǒng)計(jì)算結(jié)果的精度和效率分析

續(xù)表2

為了橫向?qū)Ρ缺疚母倪M(jìn)方法與其他文獻(xiàn)中的改進(jìn)方法的計(jì)算精度，同樣以該算例的第一個(gè)自由度為例，計(jì)算結(jié)果對(duì)比如圖7所示.綜合對(duì)比有阻尼和無(wú)阻尼條件下各方法的計(jì)算結(jié)果發(fā)現(xiàn)，本文的Wilson-Newmark法與文獻(xiàn)[14]的改進(jìn)方法計(jì)算結(jié)果相近，都能夠較好的貼合精確解曲線；而文獻(xiàn)[12]及文獻(xiàn)[13]的改進(jìn)方法計(jì)算精度要稍遜于前兩種方法，且難以反映曲線細(xì)微的變化.此外，根據(jù)表3中不同方法計(jì)算精度和效率的對(duì)比可以看出，Wilson-Newmark法的計(jì)算誤差與文獻(xiàn)[14]改進(jìn)方法的計(jì)算誤差基本一致，但相對(duì)的計(jì)算花費(fèi)時(shí)間要少于文獻(xiàn)[14].文獻(xiàn)[15]的計(jì)算精度同樣要稍低于本文所述方法，不過(guò)該方法的計(jì)算時(shí)間較短.因此，可以認(rèn)為本文的改進(jìn)方法具有一定的優(yōu)越性，能夠有效減少計(jì)算誤差，提高計(jì)算精度.

圖7 不同文獻(xiàn)方法多自由度簡(jiǎn)諧振動(dòng)系統(tǒng)計(jì)算結(jié)果比較Fig.7 Calculation results comparison of multiple degrees of simple harmonic vibration calculated by different literature methods

表3 各文獻(xiàn)計(jì)算結(jié)果的精度和效率分析

續(xù)表3

5 結(jié) 論

本文對(duì)傳統(tǒng)Wilson-θ法進(jìn)行改進(jìn)，以擴(kuò)大其在結(jié)構(gòu)動(dòng)力問(wèn)題求解中的適用范圍，提高其精度和穩(wěn)定性.所提出的改進(jìn)法基于計(jì)算擾動(dòng)的假設(shè)，并應(yīng)用Newmark-β法，修正了傳統(tǒng)Wilson-θ法不能滿足動(dòng)平衡方程的缺陷.穩(wěn)定性分析表明，改進(jìn)后的方法擴(kuò)展了原始方法無(wú)條件穩(wěn)定性的取值范圍.精度分析表明，與傳統(tǒng)的Wilson-θ法和Newmark-β法相比，改進(jìn)的Wilson-Newmark法具有振幅衰減較小、周期延長(zhǎng)也較小的優(yōu)點(diǎn)，計(jì)算精度大大提高；與其他參考文獻(xiàn)中的修正方法相比，Wilson-Newmark法也能保持一個(gè)較高的精度.算例結(jié)果驗(yàn)證了Wilson-Newmark法的計(jì)算時(shí)間雖然比原方法略有增加，但其計(jì)算精度能夠得到較大的改善.