有限元理論基礎及Abaqus內部實現方式研究系列33：線性瞬態動力學

SnowWave02

2021年1月26日 18:18

瀏覽：3801 評論：10 收藏：15

（原創，轉載請注明出處）

1 概述

本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過

（1）基礎理論

（2）商軟操作

（3）自編程序

三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。

有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。

有限元理論基礎及Abaqus內部實現方式研究系列33：線性瞬態動力學的圖2

一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元求解器，通過自研求解器和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。

通用結構有限元軟件iSolver介紹視頻：

有限元理論基礎及Abaqus內部實現方式研究系列33：線性瞬態動力學的圖4

http://www.yqgqt.org.cn/college/video/c12884

==第32篇：線性瞬態動力學==

在系列文章第三十篇、三十二篇時，我們介紹了穩態動力學的原理和算法，本章繼續介紹瞬態動力學。穩態和瞬態都是在隨時間變化的激勵作用下產生的運動形式，其中，穩態是時間足夠長后，系統慢慢趨于穩定后的運動，譬如當不考慮阻尼時彈簧受瞬間載荷力作用時間足夠長后最后會做簡諧振動，由于簡諧振動的理論公式比較簡單，一個簡諧振動在時間域上雖然很長，但只要有頻率和振幅參數就能完全表示，所以可以跳過前面的不穩定的過程而直接求解最終的穩態形式。但達到穩態前的瞬態過程中間每個時刻點的運動狀態是不一樣的，無法簡單的用幾個參數表示，當前時刻的運動狀態需要前面運動狀態往前推進得到，此時只能用瞬態分析。

瞬態分析有線性和非線性之分，從線性瞬態動力學出發其實更容易理解瞬態動力學的求解方式。本章我們將簡單介紹一下瞬態動力學的求解公式，并以一個單擺例子來說明線性瞬態動力學在有限元軟件中的內部實現原理。非線性瞬態動力學的原理將在后面章節介紹。

1.1 線性瞬態動力學的理論

1.1.1 線性瞬態動力學的方程

線性瞬態動力學的運動方程和線性靜力平衡方程很類似，對于線性靜力問題載荷F和位移u是直線關系，可以非常簡單的由一個狀態類乘以一個系數推到另一個狀態。因為靜力學可以簡單的表示為下方的線性表達形式：

Ku=R

譬如下面的懸臂梁問題，載荷1N的時候計算得到最大位移時1.177mm。

如果載荷是1000N時，不用計算也知道就是最大位移會變為1000倍=1177mm。

但在線性動力學問題中分析復雜一點就是，如果力擴大了1000倍，顯然最終的位移不是簡單的乘上1000，瞬態動力學和上式原理上很類似，只不過包括了質量陣M相關的慣性力項和阻尼陣C相關的阻尼力項，如下：

一般質量矩陣與時間和位移無關，而K和C如果和位移和時間無關，那么就是線性瞬態問題。為簡單起見，我們不考慮阻尼項。得到如下表達式：

位移和加速度與t相關，上述方程就有兩個未知數位移和加速度，只有一個方程是無法求出的，所以還有一個加速度和位移的約束關系，實際中必然是加速度先和速度相關，速度再和位移相關。

1.1.2 線性瞬態動力學的有限元求解

在理論上當前時刻的加速度就應該只和當前時刻的位移相關了，但在計算機的數值求解時，無法求出無限個連續時刻的結果，所以只能把整個時間人為的分成多個時間段，方程變為如下：

有限元理論基礎及Abaqus內部實現方式研究系列33：線性瞬態動力學的圖16

此時，加速度就必須要通過起碼兩個時刻的速度或者位移才能計算出來的，可以有不同的取法，取前時刻還是后時刻的位移來決定加速度就有隱式和顯式之分（具體可看系列文章13：顯式和隱式的區別），任何一種分析理論上都可以得到正確解，只不過顯式不用求剛度，隱式需要求解剛度陣，而通過剛度的解釋可以更容易理解Abaqus內部的實現原理，在本章中，我們選用隱式瞬態分析。

在隱式方法中，將時間離散，加速度表示為位移的表達式后，得到下方的公式

其中K和F在靜力學的增量迭代法中表示切線剛度陣和內力，只不過這邊需要加上M項：

其中A是只與時間增量步相關的量，對不同的隱式算法A的值不同，譬如對最簡單的Newmark方法中的梯形隱式算法，加速度在增量步內線性變化，下方紅線是在輸入的載荷轉換為加速度的值，但計算中沒法處理連續的曲線，所以Abaqus或者iSolver等有限元程序實際上只會認為增量步內加速度時線性變化的：

此時得到的修正剛度為：

顯然當時間增量步固定時，此時的K也是個定值。而F與上一時刻的包括慣性力在內的內力，與位移、速度、加速度相關，可以表示為：

其中B是只與時間增量步相關的量。雖然固定步長K是常量，但右端量還包括了上一個時刻位移、速度等變量，R和u依然為非線性曲線，如下：

由上面的表達式，K和F只與增量步I,j相關，但實際上只要一次迭代左右就能相等平衡了，這個原理就和線性彈性力學的內力F=KU，到時右端在第一個迭代完成后就和外力R一樣了，具體可看下面的課程推導：

http://www.yqgqt.org.cn/college/video/c14948 深入淺出有限元：基礎理論->Abaqus操作->matlab編程

也就是說線性瞬態動力學方程與迭代步j無關，無須迭代，方程可以簡單的寫成：

而且只要迭代一次，而且內力項在第一次是位移為0，那么右端F就可以簡單的寫成下方表達式：

1.2 線性瞬態動力學的單擺模型

下面我們將以一個簡單的單擺模型例子來說明上述線性瞬態動力學的解法。

1.2.1 實際單擺模型

實際的單擺如下圖，線長L，無質量，末端綁定一個質量塊m，如果沒有空氣阻力，在重力加速度g作用下往復運動。我們取初始時刻繩索在水平位置。假定繩索足夠硬，在小球作用下變形非常小。顯然，這個問題是一個大位移大轉動小應變的幾何非線性問題（具體可看我們的系列文章18：幾何非線性的應變），這個小球位置隨時間變化的運動過程就是非線性的瞬態過程。