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

周宇，汪利，劉祚秋

中山大學(xué)航空航天學(xué)院應(yīng)用力學(xué)與工程系，廣東 深圳 518107

時間有限元法是一種在時間域內(nèi)用有限單元離散的Galerkin 算法。時間有限元法的發(fā)展可追溯到20 世紀(jì)60 年代末，時間有限元法一詞最早于20世紀(jì)70年代被Clough提出［1］，F(xiàn)ried及Argyris等最早作了這方面的工作［2］。時間有限元法由于其精度高、計算量適中等特點，在結(jié)構(gòu)動力響應(yīng)的數(shù)值分析領(lǐng)域得到了廣泛應(yīng)用。

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

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

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

1 時間有限元算法原理

在結(jié)構(gòu)動力分析中，需要求解二階常微分方程組

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

其中L2(IT)表示在區(qū)間上平方可積的函數(shù)，H2(IT)表示二階導(dǎo)數(shù)平方可積函數(shù)。然后，建立結(jié)構(gòu)動力問題的變分公式

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

有了位移響應(yīng)函數(shù)，根據(jù)變分公式建立時間有限元方程

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

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

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

2 算法的穩(wěn)定性與周期延長率分析

2.1 穩(wěn)定性分析

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

其中A為傳遞矩陣，要保證算法穩(wěn)定，需滿足

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

2.2 周期延長率

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

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

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

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

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

3 算例驗證

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

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

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

圖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，且固定系統(tǒng)周期為Tn= 4.考慮系統(tǒng)的初始位移為1、初始速度為0，此時位移響應(yīng)的解析解為u(t) = cos(πt/2)。時間有限元法得到的計算結(jié)果見圖3。從圖中可以看出，改變時間步長基本不影響體系的計算振動周期，且計算振動周期與理論自振周期基本一致。為了對周期延長率進行定量分析，用隨機子空間法得到計算數(shù)值解的工程頻率fn，畫出穩(wěn)定圖，如圖4 所示?？梢钥闯?，當(dāng)模型階次為2時，隨機子空間法便能得到滿意的頻率結(jié)果（接近0.25）。按照此思路，改變時間步長，取模型階次為2，可得到周期延長率PE 隨時間步長的變化曲線，如圖5 所示。從圖中可以看出，不同時間步長下，周期延長率幾乎為0 （全小于0.01），表明時間有限元不會改變解的周期，在長時間計算中，不會引起周期錯配。

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

圖4 隨機子空間法的穩(wěn)定圖Fig.4 Stabilization diagram of stochastic subspace identification

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

4 應(yīng)用于梁的振動分析

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

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

式中l(wèi)是單元長度。類似的，單元質(zhì)量矩陣為

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

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

圖7 τ=1/8時間有限元求解的跨中節(jié)點橫向位移和速度響應(yīng)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法的位移和速度響應(yīng)中的最大誤差分別列于表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 結(jié) 論

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

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

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

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