基于時(shí)間有限元的彈簧擺動(dòng)力學(xué)新解法

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

(1. 哈爾濱工業(yè)大學(xué)機(jī)電工程學(xué)院，哈爾濱 150001；2. 軍事科學(xué)院國(guó)防工程研究院，洛陽(yáng) 471023；3. 北京起重運(yùn)輸機(jī)械設(shè)計(jì)研究院有限公司，北京 100007；4. 西安建筑科技大學(xué)機(jī)電工程學(xué)院，西安 710055)

彈簧擺問(wèn)題是一種剛?cè)狁詈系姆蔷€性動(dòng)力學(xué)問(wèn)題[1]，具有強(qiáng)耦合、強(qiáng)非線性等諸多特征。近年來(lái)隨著國(guó)家電力建設(shè)的進(jìn)一步推進(jìn)，輸電塔線系統(tǒng)的在惡劣環(huán)境下[2]，特別是在地震[3]，下?lián)舯┝鱗4]情況下的減振要求成為重點(diǎn)。為了研制較懸掛擺更好的減震器，催生了針對(duì)彈簧擺問(wèn)題的深入研究。張鵬等[5]發(fā)現(xiàn)彈簧擺的減振能力優(yōu)于懸掛擺，并對(duì)利用彈簧擺內(nèi)共振特性設(shè)計(jì)減震器進(jìn)行了探討。霍林生等[6]、王奇等[7]、黃正等[8]在彈簧擺內(nèi)共振特性上也開(kāi)展了一定程度的研究。此外Wu等[9]利用彈簧擺的能量耦合特性研制了低頻能量收集器。隨著彈簧擺應(yīng)用的深入，其模型自身的特性也受到了研究者多方面的關(guān)注。

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

本文主要運(yùn)用時(shí)間有限元以及矩形法作用量積分對(duì)小擺角彈簧擺進(jìn)行處理，采用保辛遞推，通過(guò)部分降低求解精度提高了計(jì)算速度，并能夠保證積分過(guò)程的長(zhǎng)時(shí)間穩(wěn)定，保辛的特性能夠避免系統(tǒng)誤差隨時(shí)間積分產(chǎn)生積累。

1 彈簧擺模型

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

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

取彈簧擺的端點(diǎn)為原點(diǎn)，在此情況下建立系統(tǒng)動(dòng)能和勢(shì)能的表達(dá)式為：

作用量積分，即拉格朗日函數(shù)對(duì)于時(shí)間的積分，表示為：

這時(shí)可以采用時(shí)間有限元對(duì)于該問(wèn)題進(jìn)行描述，對(duì)于兩個(gè)自由度x,y分別在時(shí)間區(qū)段ta,tb上進(jìn)行線性插值可以構(gòu)造一階的時(shí)間有限元函數(shù)：

進(jìn)行有限元變換的目的是，時(shí)間有限元變換是保辛的，能夠保證系統(tǒng)的能量守恒。離散化后的作用量積分可以寫成：

由于存在平方積分項(xiàng)，求解非線性方程組會(huì)耗費(fèi)大量的時(shí)間。為了提高求解速度，本文在構(gòu)造遞推式時(shí)對(duì)平方積分進(jìn)行線性化處理，采用矩形積分近似以求得近似解。單元作用量函數(shù)為：

其中，主要是對(duì)于平方積分項(xiàng)中的部分進(jìn)行處理，用于保證后文的求解過(guò)程是前項(xiàng)遞推的，因而只保留ta時(shí)刻的相關(guān)項(xiàng)，采用矩形積分方法對(duì)于該方程進(jìn)行部分離散，有：

其余部分進(jìn)行常規(guī)的有限元積分離散，有離散后的作用量格式：

其中：

采用保辛迭代格式求解，首先構(gòu)造對(duì)偶變量：

引入狀態(tài)向量并寫成遞推形式：

其中：

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

一階時(shí)間有限元與差分格式類似，不能顯著提高有限元方法的精確度。可以通過(guò)改進(jìn)時(shí)間有限元格式，增加時(shí)間有限元節(jié)點(diǎn)個(gè)數(shù)來(lái)提高計(jì)算精度，故而產(chǎn)生了二階或者階數(shù)更高的時(shí)間有限元。在維度更高的時(shí)間有限元推導(dǎo)過(guò)程中，采用矩陣形式和張量形式更為便捷，本方法中采用矩陣方法對(duì)于二階時(shí)間有限元進(jìn)行推導(dǎo)。

采用基函數(shù)N=[t2t1]進(jìn)行時(shí)間有限元插值，即：

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

由于引入了中間節(jié)點(diǎn)qi對(duì)應(yīng)，采用二階時(shí)間有限元的推導(dǎo)需要進(jìn)行節(jié)點(diǎn)凝聚以滿足遞推格式，而且可以推廣到更高階數(shù)的時(shí)間有限元上。需要注意，對(duì)于不采用有限元離散的非線性項(xiàng)部分，這部分在計(jì)算過(guò)程中不能引入中間節(jié)點(diǎn)進(jìn)行節(jié)點(diǎn)凝聚的方法進(jìn)行計(jì)算。因?yàn)橄挛慕o出的凝聚格式不能針對(duì)于非線性項(xiàng)進(jìn)行處理，故而應(yīng)采用兩端時(shí)間節(jié)點(diǎn)位置進(jìn)行離散。

節(jié)點(diǎn)凝聚格式由變分原理推導(dǎo)，由于Ak2、Ak3兩作用量?jī)H利用端點(diǎn)離散，與qi無(wú)關(guān)，故Ak的節(jié)點(diǎn)凝聚格式同Ak1的相一致，可以采用齊次方程下的節(jié)點(diǎn)凝聚法進(jìn)行推導(dǎo)。通過(guò)可知對(duì)于形式為：

的作用量，可以轉(zhuǎn)換為：

其中：

節(jié)點(diǎn)凝聚可以保證本方法構(gòu)造成保辛遞推格式，完成遞推求解。在節(jié)點(diǎn)凝聚后仍然需要構(gòu)造對(duì)偶變量并確立保辛遞推格式以完成求解。

由于該計(jì)算方法的時(shí)間有限元部分Ak1保辛，誤差的產(chǎn)生主要源于梯形積分在Ak2、Ak3處產(chǎn)生的誤差。這里引入代數(shù)精度的概念，根據(jù)文獻(xiàn)[28]，梯形積分僅具有一階精度，是導(dǎo)致整個(gè)積分算法精度較低的一個(gè)關(guān)鍵因素。時(shí)間有限元又因?yàn)椴捎玫惴ú荒茉谄涞仃囍幸敕蔷€性的vk，僅可選用vk-1相關(guān)值進(jìn)行積分，故對(duì)積分精度有一定影響。本方法中構(gòu)造了一種梯形算法以提高數(shù)值積分的積分精度。該積分算法采用ta時(shí)刻q(t)及其導(dǎo)數(shù)構(gòu)造積分格式，即：

對(duì)于Ak2、Ak3采用該方程構(gòu)造積分格式，由代數(shù)精度的定義可知其具有二階代數(shù)精度，較矩形公式有進(jìn)一步的精度提高。

2 模擬分析

為了檢驗(yàn)方法的準(zhǔn)確性，本文采用數(shù)值模擬同實(shí)驗(yàn)數(shù)據(jù)[15]進(jìn)行對(duì)比。通過(guò)計(jì)算可知λ=2，處于內(nèi)共振區(qū)間。同時(shí)與廣泛使用的諧波平衡法[22]進(jìn)行了對(duì)比，共進(jìn)行2個(gè)狀態(tài)下的數(shù)值模擬。

2.1 無(wú)阻尼情況

通過(guò)已有文獻(xiàn)[10 ? 11, 21 ? 22, 29 ? 31]，對(duì)于內(nèi)共振現(xiàn)象的數(shù)值解法通常采用不考慮阻尼的方程進(jìn)行求解。采用表1的實(shí)驗(yàn)參數(shù)進(jìn)行了模擬，采用傳遞辛矩陣法，單步時(shí)間為0.002 s時(shí)的模擬結(jié)果圖2、圖3所示。

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

圖3 擺球位移隨時(shí)間變化關(guān)系Fig.3 Relationship between swinging ball's displacement and time

表1 模擬參數(shù)1Table 1 Simulation parameters No.1

由圖2～圖4可以觀察到以下幾個(gè)現(xiàn)象：

圖4 系統(tǒng)能量偏差隨時(shí)間變化關(guān)系Fig.4 Relationship between system energy error and time

1) 對(duì)于振幅的分析：在初始狀態(tài)1下出現(xiàn)了明顯的共振現(xiàn)象，而初始狀態(tài)2下內(nèi)共振現(xiàn)象不明顯，這與實(shí)驗(yàn)現(xiàn)象一致，與文獻(xiàn)[22]中提到的只有在小角度擺動(dòng)才出現(xiàn)強(qiáng)烈耦合的現(xiàn)象是一致的，而且該共振顯現(xiàn)的周期與幅值同文獻(xiàn)[22]中采用高階Runge-Kutta法的結(jié)果也一致。這說(shuō)明，本算法能夠?qū)?nèi)共振現(xiàn)象進(jìn)行初步的預(yù)測(cè)。但實(shí)驗(yàn)中振動(dòng)幅度存在明顯的衰減現(xiàn)象，且內(nèi)共振現(xiàn)象發(fā)生的周期與實(shí)驗(yàn)不一致。認(rèn)為運(yùn)動(dòng)過(guò)程中存在阻尼對(duì)于內(nèi)共振現(xiàn)象的發(fā)生是有影響的，而實(shí)驗(yàn)過(guò)程中的阻尼主要由空氣阻尼或?qū)嶒?yàn)器材本身阻尼影響產(chǎn)生。需要進(jìn)行進(jìn)一步分析。

2) 對(duì)于頻率的分析：40 s內(nèi)的平均垂直振動(dòng)頻率為ωs=7.41Hz ， 水平振動(dòng)頻率ωp=3.65Hz，與實(shí)驗(yàn)中測(cè)定的ωs=7.27Hz ，ωp=3.63Hz分別存在1.93%與0.54%的誤差。相較于其進(jìn)行簡(jiǎn)諧運(yùn)動(dòng)近似的相對(duì)偏差4.9%與5.5%有較大提高。

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

2.2 有阻尼情況

由于考慮引入阻尼力后，考慮阻尼系數(shù)γ=0.0013的情況對(duì)于表1中模擬次號(hào)1的初始狀態(tài)進(jìn)行分析，分析結(jié)果如下：

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

圖5 擺動(dòng)曲線與方向位移同時(shí)間關(guān)系Fig.5 The relationship between direction displacement swing curve and direction displacement with time

2.3 時(shí)間有限元方法之間的對(duì)比

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

從表2可看出四種方法精度相近，采用增加階數(shù)與增加代數(shù)精度的方法并沒(méi)有顯著增加計(jì)算精度，反而顯著增加了計(jì)算時(shí)間。綜上，采用一階矩形法求解該問(wèn)題即可獲得理想的結(jié)果。

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

關(guān)于本方法的計(jì)算精度，整個(gè)推導(dǎo)過(guò)程存在2處近似，但僅采用梯形法積分是存在能量誤差的。這部分會(huì)導(dǎo)致能量計(jì)算不夠精確，是誤差的主要來(lái)源。因此對(duì)總能量偏差與步長(zhǎng)關(guān)系進(jìn)行研究。

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

圖6 總能量偏差同步長(zhǎng)關(guān)系(模擬次號(hào)1)Fig.6 Relationship between total energy deviation and step size (simulation No.1)

3 同其他非線性方程解法的對(duì)比

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

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

4 對(duì)強(qiáng)非線性問(wèn)題的適應(yīng)性

表4 大角度擺動(dòng)初始條件Table 4 Initial condition of wide-angle swing

在文獻(xiàn)[11]中提到，在兩振動(dòng)方向頻率比為1∶2的情況下，方程的解是初值敏感的，在不同的初值條件下能夠表現(xiàn)出準(zhǔn)周期性或混沌的特性。這里通過(guò)模擬次號(hào)2與模擬次號(hào)3相對(duì)照以檢驗(yàn)該算法對(duì)于混沌問(wèn)題的適應(yīng)性。采用傳遞辛矩陣法，步進(jìn)時(shí)間為0.002 s時(shí)的模擬結(jié)果如下：

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

圖7 大角度擺動(dòng)模擬狀態(tài)Fig.7 Wide-angle swing simulation

5 結(jié)論

本文通過(guò)模擬證明了時(shí)間有限元能夠處理彈簧擺的內(nèi)共振問(wèn)題。具有收斂較快、精度較高，并在長(zhǎng)尺度積分下保持穩(wěn)定的特性。

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

由于該方法體現(xiàn)的剛?cè)狁詈闲圆恍枰ㄟ^(guò)相對(duì)坐標(biāo)進(jìn)行參與，在絕對(duì)坐標(biāo)下即可實(shí)現(xiàn)。雖相較于保辛攝動(dòng)迭代法[26]精度有所下降，但更適于絕對(duì)節(jié)點(diǎn)坐標(biāo)法的求解。