單步無條件穩(wěn)定結(jié)構(gòu)動(dòng)力求解隱式算法

李常青 ，鐘志強(qiáng)，蔣麗忠 ，聶磊鑫

(1.中南大學(xué) 土木工程學(xué)院，湖南 長沙 410075；2.高速鐵路建造技術(shù)國家工程研究中心，湖南 長沙 410075)

在結(jié)構(gòu)動(dòng)力問題中，直接積分法是一個(gè)被廣泛研究的課題，并在結(jié)構(gòu)動(dòng)力反應(yīng)計(jì)算中得到廣泛應(yīng)用[1]。根據(jù)計(jì)算步的多少，直接積分法分為單步法、兩步法和多步法。單步法只需要知道上一步的運(yùn)動(dòng)狀態(tài)，即可求得本步的運(yùn)動(dòng)狀態(tài)，而兩步法和多步法的求解分別需要2個(gè)時(shí)刻和多個(gè)時(shí)刻的運(yùn)動(dòng)情況。根據(jù)運(yùn)動(dòng)方程的耦合情況，直接積分法分為顯式法和隱式法。顯式法不需要聯(lián)立求解方程組，根據(jù)已知的信息就能求得該時(shí)刻的解，具有更高的計(jì)算效率；而隱式法需要求解耦聯(lián)的方程組，計(jì)算工作量大。但是顯式法通常是條件穩(wěn)定的，當(dāng)時(shí)間步長超過某一臨界值時(shí)，數(shù)值解的計(jì)算誤差會(huì)瞬間增大并趨于無窮值；而隱式法大多是無條件穩(wěn)定的，對時(shí)間步長的選擇沒有限制，適用范圍更廣。中心差分法是一種經(jīng)典的顯式法，在很多情況下仍有較好的計(jì)算效果，并催生出許多優(yōu)秀的顯式算法[2-5]。本文關(guān)注隱式算法的發(fā)展、研究和應(yīng)用。比較有代表性的隱式算法有Newmark-β法[6]，Houbolt 法[7]，Wilson-θ法[8]，HHT-α法[9]，WBZ-α法[10]，Generalized-α法[11]。這些算法在穩(wěn)定性、精度、數(shù)值阻尼和超調(diào)特性上的表現(xiàn)各不相同。Newmark 平均加速度法具有較好的計(jì)算精度，但是沒有數(shù)值阻尼，在計(jì)算時(shí)由于高頻響應(yīng)的影響可能會(huì)產(chǎn)生數(shù)值振蕩。Houbolt法的高頻響應(yīng)總是漸進(jìn)衰減為0，不能調(diào)控?cái)?shù)值阻尼，同時(shí)需要借助其他方法來起步。Wilson-θ法的譜半徑在中頻段會(huì)發(fā)生突變，使得中頻段的數(shù)值阻尼大于高頻段的數(shù)值阻尼，同時(shí)在位移和速度上會(huì)產(chǎn)生明顯的超調(diào)現(xiàn)象。HHT-α法對數(shù)值耗散的調(diào)控比較有限，不能最大程度地濾除高頻分量。WBZ-α法在滿足2 階精度時(shí)是條件穩(wěn)定的，并且具有不利的超調(diào)特性。Generalized-α法能有效地調(diào)控高頻耗散并對低頻響應(yīng)的耗散較小，但是加速度是1階收斂的。上述算法有些雖是自啟動(dòng)的，但在給定初始位移和初始速度的基礎(chǔ)上，需要額外進(jìn)行初始加速度的計(jì)算，而初始加速度的計(jì)算對加速度精度有較大的影響[12-13]。LI 等[14]認(rèn)為真正自啟動(dòng)的算法是不需要初始加速度計(jì)算的，目前此類算法的研究取得了一定的進(jìn)展。SOARES 等[15-16]提出的算法具有無條件穩(wěn)定性和可控的數(shù)值阻尼，但存在不利的超調(diào)特性，無法輸出加速度響應(yīng)。KIM[17]的算法即使對于無阻尼結(jié)構(gòu)求得的加速度都只有1階精度。這些算法都有一些不理想的數(shù)值性質(zhì)，因此有必要繼續(xù)開展數(shù)值算法的研究，新算法應(yīng)在穩(wěn)定性、精度、數(shù)值阻尼、超調(diào)和自啟動(dòng)上有良好的表現(xiàn)。本文基于量綱分析的思想，基于相關(guān)文獻(xiàn)[18-20]中的思路，提出一種單步無條件穩(wěn)定隱式算法，通過理論分析和數(shù)值算例揭示其數(shù)值特性，為結(jié)構(gòu)動(dòng)力問題的求解提供新方法。

1 隱式算法的格式

1.1 隱式算法的初步格式

直接積分法研究的是離散時(shí)間點(diǎn)上的值，單自由度線性結(jié)構(gòu)的動(dòng)力方程為：

其中：m，c和k分別為結(jié)構(gòu)的質(zhì)量、阻尼和剛度；at，vt，ut分別表示結(jié)構(gòu)在t時(shí)刻的加速度、速度和位移；ft表示結(jié)構(gòu)在t時(shí)刻所受的外力。

構(gòu)造位移和速度的函數(shù)關(guān)系式，假設(shè)：

其中：Δt表示積分時(shí)間步長；ut+Δt，vt+Δt和ft+Δt分別為結(jié)構(gòu)在t+Δt時(shí)刻的位移、速度和所受的外力。根據(jù)國際單位制中的7 個(gè)基本單位，選取t時(shí)刻位移ut，時(shí)間步長Δt和質(zhì)量m作為基本量，并根據(jù)Π定理[21]，對方程(2)和(3)進(jìn)行無量綱化處理，表示為：

其中：cΔt/m和kΔt2/m是表示結(jié)構(gòu)材料屬性的參數(shù)?；谖墨I(xiàn)[18-20]的思路，引入系數(shù)ri(i=1，2，…，11) 和si(i=1，2…，11)，將方程(4)和(5)表示為：

對方程(6)和(7)化簡可得：

其中：常系數(shù)bi=ri+1/r1(i=1，2，…，10)，hi=si+1/s1(i=1，2，…，10)，方程(8)和(9)即為新算法的基本格式。

1.2 待定系數(shù)的初步確定

位移和速度在t+Δt時(shí)刻的局部截?cái)嗾`差可表示為：

其中：u(t+Δt)和v(t+Δt)分別表示結(jié)構(gòu)位移和速度在t+Δt時(shí)刻的精確解；ut+Δt和vt+Δt表示位移和速度在t+Δt時(shí)刻的數(shù)值解?；诰_解得到結(jié)構(gòu)的動(dòng)力平衡方程為：

將t+Δt時(shí)刻的位移，速度和加速度在t時(shí)刻泰勒展開，可得：

對方程(8)和(9)等式右邊的項(xiàng)進(jìn)行處理，令：

將方程(8)和(9)以及方程(12)～(17)代入方程(10)和(11)，并化簡為：

其中：di(i=0，1，…)和qi(i=0，1，…)是關(guān)于m，c，k等量的函數(shù)關(guān)系式，有：

當(dāng)位移和速度有n階精度時(shí)，則有：

根據(jù)方程(22)，當(dāng)位移和速度具有2 階精度時(shí)，有：

方程(23)只有滿足以下條件時(shí)才成立，即：

解方程(24)可得：

則方程(8)和(9)可化為：

根據(jù)方程(26)和(27)可知，位移和速度遞推格式中不含加速度項(xiàng)，其在t+Δt時(shí)刻的響應(yīng)由t時(shí)刻的位移、速度和外力以及t+Δt時(shí)刻的外力所決定，因此可以說本算法是真正自啟動(dòng)的算法。此外，如果需要得到加速度響應(yīng)，可以根據(jù)各個(gè)時(shí)刻的動(dòng)力平衡條件求得。t+Δt時(shí)刻的平衡條件為：

當(dāng)加速度有n階精度時(shí)，有：

加速度在t+Δt時(shí)刻的局部截?cái)嗾`差可表示為：

其中：ji(i=1，2，3，4)是關(guān)于m，c，k的函數(shù)關(guān)系式，則：

根據(jù)方程(29)，因此加速度具有2 階精度，方程(26)～(28)即為本算法的3 個(gè)基本公式，方程中的參數(shù)取決于算法的數(shù)值特性。

2 數(shù)值特性

2.1 收斂性

2.1.1 穩(wěn)定性

以單自由度系統(tǒng)為基本對象，分析算法的穩(wěn)定性，其t時(shí)刻的運(yùn)動(dòng)方程表示為：

其中：ξ為結(jié)構(gòu)的阻尼比；w為結(jié)構(gòu)的固有頻率。將方程(26)和(27)表示為遞推形式：

其中：A和L分別表示算法的放大矩陣和荷載矩陣。矩陣A和L中的系數(shù)表示如下：

其中：Ω=wΔt。放大矩陣A的特征方程為：

其中：E為2 階單位矩陣；λ為放大矩陣A的特征值；A1=(A11+A22)/2；A2=A11A22-A12A21。A1和A2具體如下：

放大矩陣A的譜半徑為：

算法穩(wěn)定性條件為ρ(A)≤1。HILBER 等[22]給出了算法實(shí)現(xiàn)無條件穩(wěn)定的條件，即：

求解方程(39)，可得：

2.1.2 相容性

直接積分法的相容性可以通過其局部截?cái)嗾`差[23]來分析，定義為：

將u(t+Δt)和u(t-Δt)在t時(shí)刻泰勒展開，并根據(jù)t時(shí)刻的平衡方程ma(t)+cv(t)+ku(t)=0，方程(41)可化為：

其中：gi(i=1，2，3，4)為w和ξ的函數(shù)關(guān)系式。

由方程(42)可知，e(t)為Δt的2 階小量，說明本算法是2 階相容的。根據(jù)Lax 定理，當(dāng)算法既具有相容性又具有穩(wěn)定性，則算法收斂。因此，本算法收斂，算法的精度階為相容的階數(shù)，即具有2階收斂精度。

2.2 計(jì)算精度

在結(jié)構(gòu)動(dòng)力問題中，直接積分法對結(jié)構(gòu)體系的運(yùn)動(dòng)進(jìn)行離散化，使結(jié)構(gòu)的運(yùn)動(dòng)方程在離散的時(shí)間點(diǎn)上滿足即可。然而，結(jié)構(gòu)在離散化過程中所產(chǎn)生的高階模態(tài)往往是不準(zhǔn)確的，不能代表結(jié)構(gòu)本身的運(yùn)動(dòng)狀況，因此有必要引入數(shù)值阻尼濾除高頻響應(yīng)的貢獻(xiàn)，同時(shí)也要盡可能地降低對低頻分量的影響。文獻(xiàn)[11]中指出，當(dāng)放大矩陣的特征值為一對共軛復(fù)根時(shí)，算法對高階模態(tài)的耗散最大，同時(shí)隨著Ω 趨于無窮大，特征值應(yīng)為實(shí)數(shù)。為了使本算法具有耗散特性，放大矩陣A的特征值可表示為：

求解方程(44)，可得：

放大矩陣A的譜半徑可表示為：

為了調(diào)節(jié)高頻耗散，引入?yún)?shù)：

其中：λ∞為非正數(shù)。ρ∞的取值影響著算法對高頻耗散的能力，ρ∞越小，算法對高頻耗散越大，當(dāng)ρ∞=0，算法對高頻分量的衰減最大。

在2.1 節(jié)通過對收斂性分析可知，本算法具有2 階理論精度。劉祥慶等[24]指出在實(shí)際計(jì)算中算法所表現(xiàn)出來的精度稱為計(jì)算精度，可以通過振幅衰減AD和周期延長PE來衡量。本文以單自由度有阻尼自由振動(dòng)為例，給出了其計(jì)算公式：

為了分析參數(shù)b5和h8對計(jì)算精度的影響，考慮無阻尼條件下，ρ∞分別取0.25，0.5 以及0.75，比較Δt/T=0.3 時(shí)振幅衰減AD和周期延長PE的值，將振幅衰減和周期延長為最小值時(shí)所對應(yīng)的b5和h8的取值記錄在表1 中。從表1 數(shù)據(jù)可知，當(dāng)ρ∞取不同值時(shí)，振幅衰減和周期延長總是同時(shí)達(dá)到最小值，此時(shí)b5和h8的值相等，因此為了使算法具有最高的計(jì)算精度，取h8=b5。

表1 振幅衰減和周期延長取最小值時(shí)b5和h8的取值Table 1 Values of b5 and h8 when the amplitude decay and period elongation take the minimum value

將h8=b5代入方程(40)和(45)中，可得：

對有阻尼情況下的周期延長進(jìn)行分析可知，當(dāng)b4=1/(4+8b5)，b10=1/4 時(shí)，周期誤差最小。因此，將各參數(shù)的取值用ρ∞表示為：

2.3 頻譜分析

對算法在無阻尼和有阻尼情況下的頻譜特性進(jìn)行分析，如圖1 和圖2 所示。從圖1 和圖2 中可以看到無論是否存在阻尼，本算法都是無條件穩(wěn)定的，并能夠調(diào)節(jié)數(shù)值阻尼。在圖1中，對無阻尼情況下本算法與Newmak 平均加速度法、Houbolt法和Wilson-θ法的頻譜特性進(jìn)行比較。其中，Newmark 平均加速度法與本算法在ρ∞=1 時(shí)的數(shù)值特性完全相同；Houbolt 法雖然具有耗散能力，其譜半徑曲線與本算法ρ∞=0 時(shí)接近，其高頻譜半徑趨于0，但對低頻分量的耗散較大，同時(shí)振幅衰減和周期延長也較大，計(jì)算精度最差；Wilson-θ法在中頻段的數(shù)值阻尼較大，對低頻分量的影響較大，對高頻分量的耗散能力較差，但其計(jì)算精度優(yōu)于Houbolt 法。本算法通過參數(shù)ρ∞調(diào)節(jié)高頻耗散，在引入數(shù)值阻尼的同時(shí)又能較大程度地降低數(shù)值阻尼對低頻響應(yīng)的影響，在不需要數(shù)值阻尼時(shí)與Newmark 平均加速度法等價(jià)，實(shí)現(xiàn)零振幅衰減和最小的周期誤差。當(dāng)ρ∞=1 時(shí)，本算法的位移、速度和加速度遞推格式與Newmark 平均加速度法完全相同。

圖1 不同算法的頻譜特性(ξ=0)Fig.1 Spectrum characteristics of different algorithms (ξ=0)

圖2 新算法與Generalized-α法的頻譜特性(ξ=0.05)Fig.2 Spectrum characteristics of the new algorithm and the Generalized-α method (ξ=0.05)

在圖2 中，對有阻尼情況下Generalized-α法與本算法的頻譜特性進(jìn)行比較。在ρ∞=1 時(shí)，本算法與Generalized-α法的譜半徑相接近，計(jì)算精度完全相同；在ρ∞=0.5 時(shí)，Generalized-α法的頻譜特性優(yōu)于本算法；在ρ∞=0 時(shí)，本算法對低頻分量的影響更小，同時(shí)具有更好的計(jì)算精度。因此，對多自由度體系而言，如果需要最大程度地衰減高頻分量，獲得更精確的結(jié)果，本算法優(yōu)于Generalized-α法。

2.4 超調(diào)特性

當(dāng)選取較大的時(shí)間步長，使得初始幾個(gè)時(shí)間步的計(jì)算結(jié)果被明顯放大的現(xiàn)象稱為超調(diào)[23]。譜半徑穩(wěn)定性條件決定算法在長時(shí)間內(nèi)的穩(wěn)定性，超調(diào)特性會(huì)影響算法在短時(shí)間的穩(wěn)定，必須防止這種現(xiàn)象的發(fā)生。本文通過對第1個(gè)時(shí)間步長內(nèi)的結(jié)果進(jìn)行分析，說明本算法的超調(diào)特性。

考慮單自由度有阻尼體系做自由振動(dòng)，對于給定的位移和速度初始條件，當(dāng)Ω →∞，本算法在第1個(gè)時(shí)間步的位移、速度和加速度可表示為：

由方程(54)可知，第1 個(gè)時(shí)間步的計(jì)算結(jié)果與時(shí)間步長無關(guān)，說明本算法的位移、速度和加速度都不會(huì)發(fā)生超調(diào)。

為了比較算法的數(shù)值特性，定義絕對誤差：

式中：x(i)表示i時(shí)刻的精確解；xi表示i時(shí)刻的數(shù)值解；x可取為u，v，a。

對該單自由度體系，阻尼比取ξ=0.5，自振頻率ω=2π，初始條件u0=1 m，v0=1 m/s，第1 個(gè)時(shí)間步長的計(jì)算結(jié)果如圖3 所示。從圖3 可知，對于給定的ρ∞，隨著Δt的增大，位移、速度和加速度的絕對誤差都趨于有限值，說明本算法不會(huì)出現(xiàn)超調(diào)行為，與理論分析結(jié)果相一致。

圖3 第1個(gè)時(shí)間步相對于積分步長Δt的絕對誤差Fig.3 Absolute error in the first time step versus the integration stepΔt

2.5 計(jì)算流程

對于多自由體系，本算法的具體計(jì)算過程如下。

1) 選取時(shí)間步長Δt和譜半徑趨于無窮大的值ρ∞，計(jì)算參數(shù)：

2) 確定運(yùn)動(dòng)初始值{u}0和{v}0。

3) 形成剛度矩陣[K]，質(zhì)量矩陣[M]，阻尼矩陣[C]。

4) 形成等效矩陣，即：

5) 計(jì)算t+Δt時(shí)刻的等效荷載，為：

6) 計(jì)算t+Δt時(shí)刻的位移和速度，即：

7) 根據(jù)需要計(jì)算t+Δt時(shí)刻的加速度：

循環(huán)計(jì)算上述步驟5～7，即可得到線彈性體系在任一時(shí)刻的位移、速度和加速度。

3 算例分析

在本節(jié)中，通過3個(gè)數(shù)值算例對本算法的數(shù)值特性進(jìn)行分析，揭示本算法在精度和可控?cái)?shù)值阻尼的優(yōu)越性。

3.1 算例1

以單自由度無阻尼結(jié)構(gòu)為例，分析算法在無阻尼條件下的數(shù)值行為。取結(jié)構(gòu)質(zhì)量m=1 kg，結(jié)構(gòu)剛度k=16 N/m，結(jié)構(gòu)自振周期T=π/2 s，初始條件u0=1 m，v0=1 m/s，時(shí)間步長Δt=1/20T，分析不同算法的絕對誤差，如圖4 所示。由于Houbolt法無法實(shí)現(xiàn)自啟動(dòng)，本文均采用中心差分法計(jì)算其前2 個(gè)時(shí)間步的解。從圖4 可知，Newmark 平均加速度法與本算法在ρ∞=0 時(shí)的計(jì)算結(jié)果完全相同；Houbolt 法的位移、速度和加速度絕對誤差都是最大的，其計(jì)算精度最差；Wilson-θ法得到的結(jié)果介于本算法ρ∞=0.25 和ρ∞=0.5 之間，其誤差較大。圖4的結(jié)果與頻譜分析的理論結(jié)果相一致，表明本算法通過調(diào)節(jié)參數(shù)ρ∞，不僅能使算法轉(zhuǎn)化為Newmark 平均加速度法，同時(shí)優(yōu)于Houbolt 法和Wilson-θ法。

圖4 不同算法的絕對誤差Fig.4 Absolute error of different algorithms

3.2 算例2

以單自由度有阻尼結(jié)構(gòu)為例，比較本算法與Generalized-α 法的數(shù)值特性。取結(jié)構(gòu)質(zhì)量m=1 kg，結(jié)構(gòu)剛度k=5 N/m，結(jié)構(gòu)自振周期阻尼比，結(jié)構(gòu)外力ft=sin(2t)，初始條件u0=57/65 m，v0=2/65 m/s。結(jié)構(gòu)的位移精確解可表示為：

為了比較算法的數(shù)值特性，定義全局誤差：

式中：Nt表示總計(jì)算時(shí)間；x(i)表示i時(shí)刻的精確解；xi表示i時(shí)刻的數(shù)值解。x可取為u，v，a。

總計(jì)算時(shí)間取Nt=10T，得到全局誤差隨時(shí)間步長的變化關(guān)系，如圖5 所示。從圖5 可知，對于有阻尼結(jié)構(gòu)，無論ρ∞取何值，本算法的位移、速度和加速度都具有2 階收斂精度，同時(shí)Generalized-α法的位移和速度也具有2 階收斂精度，而Generalized-α法只有當(dāng)ρ∞=1 時(shí)加速度才具有2 階收斂精度，其他情況下只有1 階收斂精度。此外，當(dāng)ρ∞=0，本算法計(jì)算得到的全局誤差明顯小于Generalized-α法；當(dāng)ρ∞=0.5時(shí)，2種算法的位移和速度全局誤差相差不大，而Generalized-α法的加速度全局誤差明顯偏大；當(dāng)ρ∞=1 時(shí)，2 種算法的位移、速度和加速度全局誤差十分接近。整體而言，本算法計(jì)算得到的加速度精度更高，誤差更小，同時(shí)當(dāng)ρ∞的取值接近于0，本算法的數(shù)值性能明顯優(yōu)于Generalized-α法。

圖5 新算法與Generalized-α法的全局誤差Fig.5 Global error of the new algorithm and the Generalized-α method

3.3 算例3

以懸臂桿[25]為對象，在t=0 時(shí)施加端部荷載F(t)，如圖6所示。其位移u(x，t)運(yùn)動(dòng)方程為：

圖6 受端部荷載的懸臂桿Fig.6 A clamped-free bar excited by end load

式中：E為彈性模量；A為截面面積；ρ為質(zhì)量密度。約束條件為：

式中：L為桿件長度。將桿件離散為N個(gè)等長的2節(jié)點(diǎn)單元，單元長度le=L/N，采用一致質(zhì)量矩陣，其單元?jiǎng)偠染仃嚭蛦卧|(zhì)量矩陣為：

此算例中L=200 m，F(xiàn)(t)=1×104N，E=5×107Pa，A=1 m2，ρ=8×10-4kg/m3，N=1 000，時(shí)間步長Δt=4.72×10-7s。

采用Newmark 平均加速度法求解，得到桿件中間節(jié)點(diǎn)位移和速度隨時(shí)間的變化情況，如圖7所示，圖中Reference 表示理論解。圖7 中位移和速度均產(chǎn)生數(shù)值振蕩，有必要引入數(shù)值阻尼來抑制雜散模態(tài)[26]。采用Generalized-α法和本算法在ρ∞=0.75時(shí)求解，結(jié)果如圖8所示。與Newmark平均加速度法相比，2 種具有數(shù)值阻尼的算法都能提供更精確的位移解，都能在較大程度上抑制數(shù)值振蕩，明顯提高速度的求解精度，但本算法在位移和速度的求解中表現(xiàn)出更好的計(jì)算精度，計(jì)算結(jié)果更接近理論解。說明本算法通過引入數(shù)值阻尼，在濾除高頻響應(yīng)的同時(shí)，更大程度地保留低階振型的貢獻(xiàn)，具有更優(yōu)的數(shù)值特性。

圖7 Newmark平均加速度法的響應(yīng)Fig.7 Response of the Newmark average acceleration method

圖8 新算法與Generalized-α法的響應(yīng)Fig.8 Response of the new algorithm and the Generalized-α method

4 結(jié)論

1) 本文算法具有無條件穩(wěn)定、2 階精度和無超調(diào)特性，并通過調(diào)整參數(shù)值能轉(zhuǎn)變?yōu)镹ewmark 平均加速度法，同時(shí)優(yōu)于Houbolt法和Wilson-θ法。

2) 本文算法不需要加速度參與，能實(shí)現(xiàn)真正的自啟動(dòng)，并能提供2階精度加速度輸出。

3) 本文算法具有可控?cái)?shù)值阻尼特性，與Generalized-α法相比，該算法在有阻尼條件下，位移、速度和加速度均有2階精度，并能提供更好的計(jì)算精度。