基于時間有限元的彈簧擺動力學新解法

陳德坤，丁幸波，溫新婕，劉曼蘭，蘭 朋

(1. 哈爾濱工業大學機電工程學院，哈爾濱 150001；2. 軍事科學院國防工程研究院，洛陽 471023；3. 北京起重運輸機械設計研究院有限公司，北京 100007；4. 西安建筑科技大學機電工程學院，西安 710055)

彈簧擺問題是一種剛柔耦合的非線性動力學問題[1]，具有強耦合、強非線性等諸多特征。近年來隨著國家電力建設的進一步推進，輸電塔線系統的在惡劣環境下[2]，特別是在地震[3]，下擊暴流[4]情況下的減振要求成為重點。為了研制較懸掛擺更好的減震器，催生了針對彈簧擺問題的深入研究。張鵬等[5]發現彈簧擺的減振能力優于懸掛擺，并對利用彈簧擺內共振特性設計減震器進行了探討?；袅稚萚6]、王奇等[7]、黃正等[8]在彈簧擺內共振特性上也開展了一定程度的研究。此外Wu等[9]利用彈簧擺的能量耦合特性研制了低頻能量收集器。隨著彈簧擺應用的深入，其模型自身的特性也受到了研究者多方面的關注。

鑒于傳統方法存在一定缺陷，針對于哈密頓系統自身特性求解的時間有限元方法應運而生。由于有限元方法是基于變分原理獲得的，故保辛。時間有限元最早是由Zienkiewicz等[23]提出的對于時間坐標的有限元，但并未要求保辛[24]。由于保守體系可用Hamilton體系描述，而保守體系的差分格式應當保辛[25]。故鐘萬勰[24]提出了保辛時間有限元的計算方法對于線性問題進行處理，體現了該方法的可靠性。吳峰等[26]利用了保辛攝動算法求解了剛-柔性桿的耦合問題，顯示了保辛算法在處理小變形幾何非線性問題的優勢。

本文主要運用時間有限元以及矩形法作用量積分對小擺角彈簧擺進行處理，采用保辛遞推，通過部分降低求解精度提高了計算速度，并能夠保證積分過程的長時間穩定，保辛的特性能夠避免系統誤差隨時間積分產生積累。

1 彈簧擺模型

時間有限元與保辛遞推求解不同于傳統的微分方程求解，其從動力方程的Hamilton變分原理出發，對系統作用量進行時間離散。本文以繞支點無阻尼自由擺動的彈簧擺為對象進行研究，研究模型如圖1所示。采用直角坐標系描述彈簧擺的空間位置，設定二自由度坐標分別為x與y，質點質量m，彈簧在質點自重作用下的長度為l0，彈簧自然長度為r0，彈簧勁度系數為k，重力加速度為g，單位全部采用國際單位制。

圖1 彈簧擺示意圖Fig.1 Schematic diagram of spring pendulum

取彈簧擺的端點為原點，在此情況下建立系統動能和勢能的表達式為：

作用量積分，即拉格朗日函數對于時間的積分，表示為：

這時可以采用時間有限元對于該問題進行描述，對于兩個自由度x,y分別在時間區段ta,tb上進行線性插值可以構造一階的時間有限元函數：

進行有限元變換的目的是，時間有限元變換是保辛的，能夠保證系統的能量守恒。離散化后的作用量積分可以寫成：

由于存在平方積分項，求解非線性方程組會耗費大量的時間。為了提高求解速度，本文在構造遞推式時對平方積分進行線性化處理，采用矩形積分近似以求得近似解。單元作用量函數為：

其中，主要是對于平方積分項中的部分進行處理，用于保證后文的求解過程是前項遞推的，因而只保留ta時刻的相關項，采用矩形積分方法對于該方程進行部分離散，有：

其余部分進行常規的有限元積分離散，有離散后的作用量格式：

其中：

采用保辛迭代格式求解，首先構造對偶變量：

引入狀態向量并寫成遞推形式：

其中：

經過前面的矩形積分近似最終可以處理為S(xa,ya)的形式，可以代入初始狀態向量進行遞推求解。

一階時間有限元與差分格式類似，不能顯著提高有限元方法的精確度?？梢酝ㄟ^改進時間有限元格式，增加時間有限元節點個數來提高計算精度，故而產生了二階或者階數更高的時間有限元。在維度更高的時間有限元推導過程中，采用矩陣形式和張量形式更為便捷，本方法中采用矩陣方法對于二階時間有限元進行推導。

采用基函數N=[t2t1]進行時間有限元插值，即：

采用時間有限元法，即Ak1采用二階時間有限元方法離散，Ak2、Ak3采用其他數值方法離散，這里要求的數值離散方法需要能夠滿足前文提到的遞推要求，因而僅能通過ta端點進行構造。在本解法中采用矩形積分方法進行離散。

由于引入了中間節點qi對應，采用二階時間有限元的推導需要進行節點凝聚以滿足遞推格式，而且可以推廣到更高階數的時間有限元上。需要注意，對于不采用有限元離散的非線性項部分，這部分在計算過程中不能引入中間節點進行節點凝聚的方法進行計算。因為下文給出的凝聚格式不能針對于非線性項進行處理，故而應采用兩端時間節點位置進行離散。

節點凝聚格式由變分原理推導，由于Ak2、Ak3兩作用量僅利用端點離散，與qi無關，故Ak的節點凝聚格式同Ak1的相一致，可以采用齊次方程下的節點凝聚法進行推導。通過可知對于形式為：

的作用量，可以轉換為：

其中：

節點凝聚可以保證本方法構造成保辛遞推格式，完成遞推求解。在節點凝聚后仍然需要構造對偶變量并確立保辛遞推格式以完成求解。

由于該計算方法的時間有限元部分Ak1保辛，誤差的產生主要源于梯形積分在Ak2、Ak3處產生的誤差。這里引入代數精度的概念，根據文獻[28]，梯形積分僅具有一階精度，是導致整個積分算法精度較低的一個關鍵因素。時間有限元又因為采用迭代算法不能在其迭代矩陣中引入非線性的vk，僅可選用vk-1相關值進行積分，故對積分精度有一定影響。本方法中構造了一種梯形算法以提高數值積分的積分精度。該積分算法采用ta時刻q(t)及其導數構造積分格式，即：

對于Ak2、Ak3采用該方程構造積分格式，由代數精度的定義可知其具有二階代數精度，較矩形公式有進一步的精度提高。

2 模擬分析

為了檢驗方法的準確性，本文采用數值模擬同實驗數據[15]進行對比。通過計算可知λ=2，處于內共振區間。同時與廣泛使用的諧波平衡法[22]進行了對比，共進行2個狀態下的數值模擬。

2.1 無阻尼情況

通過已有文獻[10 ? 11, 21 ? 22, 29 ? 31]，對于內共振現象的數值解法通常采用不考慮阻尼的方程進行求解。采用表1的實驗參數進行了模擬，采用傳遞辛矩陣法，單步時間為0.002 s時的模擬結果圖2、圖3所示。

圖2 擺球在兩次模擬下的擺動曲線Fig.2 Swing curve of swinging ball in two simulations

圖3 擺球位移隨時間變化關系Fig.3 Relationship between swinging ball's displacement and time

表1 模擬參數1Table 1 Simulation parameters No.1

由圖2～圖4可以觀察到以下幾個現象：

圖4 系統能量偏差隨時間變化關系Fig.4 Relationship between system energy error and time

1) 對于振幅的分析：在初始狀態1下出現了明顯的共振現象，而初始狀態2下內共振現象不明顯，這與實驗現象一致，與文獻[22]中提到的只有在小角度擺動才出現強烈耦合的現象是一致的，而且該共振顯現的周期與幅值同文獻[22]中采用高階Runge-Kutta法的結果也一致。這說明，本算法能夠對內共振現象進行初步的預測。但實驗中振動幅度存在明顯的衰減現象，且內共振現象發生的周期與實驗不一致。認為運動過程中存在阻尼對于內共振現象的發生是有影響的，而實驗過程中的阻尼主要由空氣阻尼或實驗器材本身阻尼影響產生。需要進行進一步分析。

2) 對于頻率的分析：40 s內的平均垂直振動頻率為ωs=7.41Hz ， 水平振動頻率ωp=3.65Hz，與實驗中測定的ωs=7.27Hz ，ωp=3.63Hz分別存在1.93%與0.54%的誤差。相較于其進行簡諧運動近似的相對偏差4.9%與5.5%有較大提高。

3) 對于機械能守恒效果的分析：由于對保守系統進行分析，觀察總能量的變化規律是評判算法好壞的重要指標。在該算法中能量偏差存在著周期性的變化，在計算步長為0.002 s時可以保證能量偏差與輸入能量初值比維持在5‰以下。在計算過程中可以看出總能量偏差同初始能量之間的變化關系是周期性質的，而不像[10]中提到的對于微分方程求解方法的能量誤差會出現發散的現象。這能夠保證能量積分的相對穩定，對于長歷程積分有很好適應性。

2.2 有阻尼情況

由于考慮引入阻尼力后，考慮阻尼系數γ=0.0013的情況對于表1中模擬次號1的初始狀態進行分析，分析結果如下：

圖5反映了含有阻尼力的擺動曲線與實驗曲線更為類似，也更符合實驗曲線的衰減規律。此外在增加系統阻尼的情況下，模擬軌跡出現了較之前的一些不同。首先，是阻尼系數導致振動頻率的降低，增加阻尼后，兩個方向振動頻率分別降低至3.62 Hz與7.24 Hz，與實驗計算的頻率偏差分別為0.27%與0.41%，相較于不加阻尼分別提高了1.66%與0.13%，效果很好。其次，是阻尼系數的引入能夠導致內共振現象發生頻率產生變化，對于模擬次號1實驗，內共振現象峰值出現頻數由5次降低至4次，說明阻尼對內共振現象確實存在較大影響。最后，通過多次改變初始條件的值可以發現，通過使初值在一定范圍(x方向不變，y方向0.01 m～0.1 m)內浮動可以發現，內共振現象的發生是依賴與初始條件的，且這種依賴現象同阻尼關系不大。在計算中可以發現，在y在區間0.05 m～0.06 m時，內共振現象很難發生，x方向與y方向幾乎相對獨立的進行振動。

圖5 擺動曲線與方向位移同時間關系Fig.5 The relationship between direction displacement swing curve and direction displacement with time

2.3 時間有限元方法之間的對比

除一階矩形法(僅對Ak2使用)外，一階梯形(對Ak2、Ak3使用)法，二階矩形法(僅對Ak2使用)以及二階梯形法(對Ak2、Ak3使用)均給出了類似的結果，并對各組結果進行比較。表2采用計算步長0.2 ms，計算總時長20 s，分別對于第一次實驗和第二次實驗，針對于這四種方法的運算速度以及計算精度通過下表進行比較。

從表2可看出四種方法精度相近，采用增加階數與增加代數精度的方法并沒有顯著增加計算精度，反而顯著增加了計算時間。綜上，采用一階矩形法求解該問題即可獲得理想的結果。

表2 積分方式對于模擬效果的影響Table 2 The influence of integral method on simulation

關于本方法的計算精度，整個推導過程存在2處近似，但僅采用梯形法積分是存在能量誤差的。這部分會導致能量計算不夠精確，是誤差的主要來源。因此對總能量偏差與步長關系進行研究。

由圖6可以看出，計算步長對于總能量偏差影響較大，在采用線性時間有限元的情況下在步長高于0.002 s時總能量偏差即超過5‰，因而在計算過程中需要盡量選取低于0.002 s的步長，以保證積分的準確性。

圖6 總能量偏差同步長關系(模擬次號1)Fig.6 Relationship between total energy deviation and step size (simulation No.1)

3 同其他非線性方程解法的對比

時間有限元方法在計算精度相同的情況下同其他微分方程求解算法進行對比。計算機CPU采用四核i7-4710 MQ，主頻2.5G Hz，內存4 GB。在計算步長為10 μs時計算2 s計算時間為0.64 s，且收斂速度較其他數值方法有明顯提高。同表3中的非線性微分方程求解方法相比較，在λ =1，δ=6.2×10-7時，同4階R-K法的計算速度類似，顯著高于Ode15 s。確實在降低求解步長的情況下提高了計算速度，并保持了積分誤差的穩定性。

表3 不同頻率下數值方法比較[10]Table 3 Comparison of numerical methods at different frequencies[10]

4 對強非線性問題的適應性

表4 大角度擺動初始條件Table 4 Initial condition of wide-angle swing

在文獻[11]中提到，在兩振動方向頻率比為1∶2的情況下，方程的解是初值敏感的，在不同的初值條件下能夠表現出準周期性或混沌的特性。這里通過模擬次號2與模擬次號3相對照以檢驗該算法對于混沌問題的適應性。采用傳遞辛矩陣法，步進時間為0.002 s時的模擬結果如下：

從圖7(b)和圖7(c)的比較中可以發現在不同初值下，相圖x方向出現了一定的混沌特征，這與文獻[11]的結論一致。但從總能能量偏差的角度可以看出能量偏差雖然較小角度擺動的能量偏差更大，但在合理范圍內且并沒有出現發散趨勢，通過前文的模擬也可知通過減小步長可以大幅降低總能量偏差。該模擬表明了時間有限元的保辛遞推算法具有非線性長時間積分下穩定的特點。

圖7 大角度擺動模擬狀態Fig.7 Wide-angle swing simulation

5 結論

本文通過模擬證明了時間有限元能夠處理彈簧擺的內共振問題。具有收斂較快、精度較高，并在長尺度積分下保持穩定的特性。

在對時間有限元步長進行了誤差分析與速度分析后發現，該方法不存在誤差累積的現象，計算精度主要取決于步長，在步長為0.001 s～0.01 s的范圍內最大能量偏差小于2.17%，且能夠保證計算速度。一階矩形時間有限元的計算精度與時間經濟性更好。此外對于更一般的大角度非線性彈簧擺問題也進行了求解，得到了長時間積分下誤差可控的模擬結果。

由于該方法體現的剛柔耦合性不需要通過相對坐標進行參與，在絕對坐標下即可實現。雖相較于保辛攝動迭代法[26]精度有所下降，但更適于絕對節點坐標法的求解。