時間有限元的算法穩定性與周期延長率分析*

周宇，汪利，劉祚秋

中山大學航空航天學院應用力學與工程系，廣東 深圳 518107

時間有限元法是一種在時間域內用有限單元離散的Galerkin 算法。時間有限元法的發展可追溯到20 世紀60 年代末，時間有限元法一詞最早于20世紀70年代被Clough提出［1］，Fried及Argyris等最早作了這方面的工作［2］。時間有限元法由于其精度高、計算量適中等特點，在結構動力響應的數值分析領域得到了廣泛應用。

在算法研究領域，Hulbert［3］建立了結構動力學方程的時間不連續Galerkin 有限元法，并證明了算法的穩定性和收斂性。隨后Li 和Wiberg［4］研究了時間不連續Galerkin 有限元方法，提出了一種新的迭代求解算法。French 等［5-6］為求解二階波動方程問題，提出了一種時間連續Galerkin 有限元方法。袁曉彬等［7］采用加權殘值Galerkin法，建立了求解一階線性微分方程組的全域時間有限元法。許偉等［8］基于五次Hermite 插值函數，構造了求解結構動力學二階線性微分方程組的全域時間有限元算法。就在最近，Wang 和Zhong［9］構造出與傳統的強控制方程形式等價的變分公式，以此建立了時間連續Galerkin有限元公式。

在算法實際應用方面，張方等［10］建立了動載荷識別的時間有限元模型，將時間有限元法應用到動載荷識別領域。榮吉利和李瑞英［11］用時間不連續Galerkin 法對水輪發電機軸系橫向振動響應進行了仿真計算。隋永楓等［12］將時間有限元法擴展到陀螺系統，證明了算法的優越性。姜燕等［13］用時間有限元預測了銑削的穩定性，并用實驗證明了時間有限元的有效性。

在作者最近的工作［9］中，已經給出了時間有限元嚴格的先驗誤差分析，并且發現與常規的Newmark法不同的是，時間有限元的計算誤差不會隨時間逐漸增大。這對工程計算十分重要，因為大部分算法中計算時間越長、誤差越大。為了將時間有限元應用于工程分析，本文將進一步研究時間有限元的算法穩定性與周期延長率等，并通過梁的動力分析算例驗證算法的精度。

1 時間有限元算法原理

在結構動力分析中，需要求解二階常微分方程組

其中M為質量矩陣，C為粘滯阻尼矩陣，K為剛度矩陣，f為荷載向量，T為計算的時間長度，u，u?分別為位移、速度，IT為時域區間。

其中L2(IT)表示在區間上平方可積的函數，H2(IT)表示二階導數平方可積函數。然后，建立結構動力問題的變分公式

其中雙線性泛函B：H2(IT) ×H2(IT) →R 和線性泛函l：H2(IT) →R表示為

有了位移響應函數，根據變分公式建立時間有限元方程

變分公式的具體推導過程以及時間有限元法的收斂性分析可參考文獻［9］。然后，計算單元剛度矩陣Di和荷載向量{f}i

實際中，常采用高斯求積法計算Di和{f}i，最后通過組裝得到整體代數方程

由于D是三階對角塊矩陣，在已知初始位移和速度的情況下，可用追趕法快速求解方程（9）。

2 算法的穩定性與周期延長率分析

2.1 穩定性分析

為分析時間有限元法的算法穩定性，針對自由振動問題（荷載f(t) = 0）建立的遞推公式為［14］

其中A為傳遞矩陣，要保證算法穩定，需滿足

其中ρ(A)為譜半徑，λ(A)為傳遞矩陣的特征值。

2.2 周期延長率

當系統存在一定的周期延長時，計算響應與實際響應將在長時間后存在周期錯配，使得長期行為難以預測。周期延長率的計算公式為［15］

式中T為體系的理論自振周期；T′為數值解的振動周期。

實際分析中，得到的數值解并不是完美的簡諧函數，如何獲得數值解的周期是計算周期延長率的關鍵。采用隨機子空間法［16］，根據離散的位移數據直接找出周期。首先，對時間有限元的解插值得到間隔為Δt的位移值-u1，-u2，…，-un。然后，利用離散的位移值，構造一個Hankel矩陣，即

其中m+n為時間節點總個數。對矩陣H進行奇異值分解（SVD），得［17］

其中Si為H的奇異值，r為H的秩。根據隨機子空間識別理論［16］，得到離散系統特征矩陣Q，對Q進行特征值分解，得到特征值λ1，λ2，…，λn。最終，可得到離散數據的頻率或周期

3 算例驗證

考慮如下單自由度體系的自由振動問題

其中ξ為阻尼比，ωn為自振頻率，Tn= 2π/ωn為系統周期，計算時長T=12 s.

一方面，進行算法的穩定性分析。主要考慮時間步長τ，阻尼比ξ，周期Tn（或頻率）等因素的影響：1）固定τ= 0.25 s，阻尼比ξ=0～1%和不同周期Tn下，傳遞矩陣譜半徑隨步長周期比（τ/Tn）的變化見圖1；2） 固定τ= 0.5 s，阻尼比ξ=0～1%和不同周期Tn下，傳遞矩陣譜半徑隨τ/Tn的變化見圖2。可以看出，阻尼比ξ和τ/Tn顯著影響譜半徑的值；而給定ξ和τ/Tn，譜半徑幾乎與時間步長τ沒有關系。當存在阻尼（ξ＞0.05%）時，譜半徑在合理的步長周期比范圍(τ/Tn＜1)內小于1，表明時間有限元計算是穩定的；而對于無阻尼或ξ＜0.05%的系統，只有當τ/Tn＜0.3時，譜半徑小于1，此時屬于條件穩定狀態。

圖1 τ=0.25時譜半徑隨阻尼比ξ和τ/Tn的變化Fig.1 Variation of spectral radius with damping ratio ξ and τ/Tn under τ=0.25

圖2 τ=0.5時譜半徑隨阻尼比ξ和τ/Tn的變化Fig.2 Variation of spectral radius with damping ratio ξ and τ/Tn under τ=0.5

另一方面，為了研究時間有限元法的周期延長率，令阻尼比ξ為0，且固定系統周期為Tn= 4.考慮系統的初始位移為1、初始速度為0，此時位移響應的解析解為u(t) = cos(πt/2)。時間有限元法得到的計算結果見圖3。從圖中可以看出，改變時間步長基本不影響體系的計算振動周期，且計算振動周期與理論自振周期基本一致。為了對周期延長率進行定量分析，用隨機子空間法得到計算數值解的工程頻率fn，畫出穩定圖，如圖4 所示。可以看出，當模型階次為2時，隨機子空間法便能得到滿意的頻率結果（接近0.25）。按照此思路，改變時間步長，取模型階次為2，可得到周期延長率PE 隨時間步長的變化曲線，如圖5 所示。從圖中可以看出，不同時間步長下，周期延長率幾乎為0 （全小于0.01），表明時間有限元不會改變解的周期，在長時間計算中，不會引起周期錯配。

圖3 不同時間步長時單自由度無阻尼自由振動系統的位移響應Fig.3 Displacement response of single degree of freedom undamped vibration system with different time steps

圖4 隨機子空間法的穩定圖Fig.4 Stabilization diagram of stochastic subspace identification

圖5 周期延長率隨時間步長的變化Fig.5 Change trend of period elongation with time step

4 應用于梁的振動分析

考慮承受集中荷載的等截面簡支梁模型（如圖6 所示）。荷載距左端支座的距離為s（s=L/4），大小為P(t) =Fsin(ωt)。研究的響應時間區間為[0，12]，模型的基本參數是：質量分布m= 1 000 kg/m，梁長L= 10 m，彈性模量E= 3 × 109Pa，截面慣性矩I= 6 × 10-7m3，外荷載激振頻率ω=π/2 rad/s以及荷載幅值F= 1 000 N。梁的初始橫向位移和轉角均設為0。為了建立梁的動力方程，首先對梁進行離散化，均勻劃分為8個單元。選取結點的橫向位移和轉角為自由度，單元剛度矩陣為

圖6 受集中荷載的等截面簡支梁模型Fig.6 Model of a beam with equal section under concentrated load

式中l是單元長度。類似的，單元質量矩陣為

將單元剛度矩陣、質量矩陣和外力向量進行組裝，集成結構體系的總體剛度矩陣K、質量矩陣M和外力荷載向量f(t)，再采用阻尼理論假設得到體系的總體阻尼矩陣C，形成如式（1）的梁振動方程。然后，采用時間有限元法進行求解。

取跨中節點的橫向位移為研究對象，可以得到跨中節點的橫向位移的時間有限元解，如圖7所示。從圖中可以看出時間有限元的數值解與方程的解析解十分吻合。當然，為了對比，額外應用Newmark 法（γ=0.5，β=0.25） 對該問題進行求解。時間步長τ=1/8 下，時間有限元和Newmark 法在所有時間節點上位移和速度的逐點誤差如圖8和圖9所示。從圖中可以看出，時間有限元法的逐點誤差明顯比Newmark 法小。Newmark 法的誤差會隨著時間增加，而時間有限元的逐點誤差卻一直保持在較低的水平。

圖7 τ=1/8時間有限元求解的跨中節點橫向位移和速度響應Fig.7 Transverse displacement and velocity response of mid-span nodes solved by time finite element under τ=1/8

圖8 時間有限元法和Newmark法的位移誤差對比Fig.8 Comparison of errors in pointwise displacement by time FEM and Newmark method

圖9 時間有限元法和Newmark法的速度誤差對比Fig.9 Comparison of errors in pointwise velocity by time FEM and Newmark method

為了量化計算誤差，將時間有限元和Newmark法的位移和速度響應中的最大誤差分別列于表1和表2中。可以看出，相同時間步長下，時間有限元法位移的最大誤差僅為Newmark 法的6.67%～9.6%，速度誤差不超過Newmark 法的5.78%。時間有限元法的精度顯然比Newmark法高。

表1 時間有限元和Newmark法跨中位移的最大誤差Table 1 Maximum error in mid-span displacement of time finite element and Newmark method

表2 時間有限元和Newmark法跨中速度的最大誤差Table 2 Maximum error in mid-span velocity of time finite element and Newmark method

5 結 論

本文研究了時間有限元法的算法穩定性和周期延長率問題，得到了以下結論：

1） 對于有阻尼（阻尼比大于0.05%）單自由度系統，當ξ和τ/Tn＜1 取不同值時，譜半徑ρ(A)都始終小于1，表明時間有限元法的穩定性很好。對于無阻尼（或者阻尼比小于0.05%）單自由度系統，τ/Tn＜0.3時系統達到條件穩定。

2） 時間有限元法的周期延長率（幾乎）為0，可以認為計算振動周期等于理論自振周期，不存在周期延長的現象。

3） 時間有限元法在求解多自由度線性系統的動力反應問題時，其計算精度要遠高于一般的Newmark法。