單步無條件穩定結構動力求解隱式算法

李常青 ，鐘志強，蔣麗忠 ，聶磊鑫

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

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

1 隱式算法的格式

1.1 隱式算法的初步格式

直接積分法研究的是離散時間點上的值，單自由度線性結構的動力方程為：

其中：m，c和k分別為結構的質量、阻尼和剛度；at，vt，ut分別表示結構在t時刻的加速度、速度和位移；ft表示結構在t時刻所受的外力。

構造位移和速度的函數關系式，假設：

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

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

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

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

1.2 待定系數的初步確定

位移和速度在t+Δt時刻的局部截斷誤差可表示為：

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

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

對方程(8)和(9)等式右邊的項進行處理，令：

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

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

當位移和速度有n階精度時，則有：

根據方程(22)，當位移和速度具有2 階精度時，有：

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

解方程(24)可得：

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

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

當加速度有n階精度時，有：

加速度在t+Δt時刻的局部截斷誤差可表示為：

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

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

2 數值特性

2.1 收斂性

2.1.1 穩定性

以單自由度系統為基本對象，分析算法的穩定性，其t時刻的運動方程表示為：

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

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

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

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

放大矩陣A的譜半徑為：

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

求解方程(39)，可得：

2.1.2 相容性

直接積分法的相容性可以通過其局部截斷誤差[23]來分析，定義為：

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

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

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

2.2 計算精度

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

求解方程(44)，可得：

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

為了調節高頻耗散，引入參數：

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

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

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

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

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

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

2.3 頻譜分析

對算法在無阻尼和有阻尼情況下的頻譜特性進行分析，如圖1 和圖2 所示。從圖1 和圖2 中可以看到無論是否存在阻尼，本算法都是無條件穩定的，并能夠調節數值阻尼。在圖1中，對無阻尼情況下本算法與Newmak 平均加速度法、Houbolt法和Wilson-θ法的頻譜特性進行比較。其中，Newmark 平均加速度法與本算法在ρ∞=1 時的數值特性完全相同；Houbolt 法雖然具有耗散能力，其譜半徑曲線與本算法ρ∞=0 時接近，其高頻譜半徑趨于0，但對低頻分量的耗散較大，同時振幅衰減和周期延長也較大，計算精度最差；Wilson-θ法在中頻段的數值阻尼較大，對低頻分量的影響較大，對高頻分量的耗散能力較差，但其計算精度優于Houbolt 法。本算法通過參數ρ∞調節高頻耗散，在引入數值阻尼的同時又能較大程度地降低數值阻尼對低頻響應的影響，在不需要數值阻尼時與Newmark 平均加速度法等價，實現零振幅衰減和最小的周期誤差。當ρ∞=1 時，本算法的位移、速度和加速度遞推格式與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-α法與本算法的頻譜特性進行比較。在ρ∞=1 時，本算法與Generalized-α法的譜半徑相接近，計算精度完全相同；在ρ∞=0.5 時，Generalized-α法的頻譜特性優于本算法；在ρ∞=0 時，本算法對低頻分量的影響更小，同時具有更好的計算精度。因此，對多自由度體系而言，如果需要最大程度地衰減高頻分量，獲得更精確的結果，本算法優于Generalized-α法。

2.4 超調特性

當選取較大的時間步長，使得初始幾個時間步的計算結果被明顯放大的現象稱為超調[23]。譜半徑穩定性條件決定算法在長時間內的穩定性，超調特性會影響算法在短時間的穩定，必須防止這種現象的發生。本文通過對第1個時間步長內的結果進行分析，說明本算法的超調特性。

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

由方程(54)可知，第1 個時間步的計算結果與時間步長無關，說明本算法的位移、速度和加速度都不會發生超調。

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

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

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

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

2.5 計算流程

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

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

2) 確定運動初始值{u}0和{v}0。

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

4) 形成等效矩陣，即：

5) 計算t+Δt時刻的等效荷載，為：

6) 計算t+Δt時刻的位移和速度，即：

7) 根據需要計算t+Δt時刻的加速度：

循環計算上述步驟5～7，即可得到線彈性體系在任一時刻的位移、速度和加速度。

3 算例分析

在本節中，通過3個數值算例對本算法的數值特性進行分析，揭示本算法在精度和可控數值阻尼的優越性。

3.1 算例1

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

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

3.2 算例2

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

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

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

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

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

3.3 算例3

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

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

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

式中：L為桿件長度。將桿件離散為N個等長的2節點單元，單元長度le=L/N，采用一致質量矩陣，其單元剛度矩陣和單元質量矩陣為：

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

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

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

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

4 結論

1) 本文算法具有無條件穩定、2 階精度和無超調特性，并通過調整參數值能轉變為Newmark 平均加速度法，同時優于Houbolt法和Wilson-θ法。

2) 本文算法不需要加速度參與，能實現真正的自啟動，并能提供2階精度加速度輸出。

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