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

==概述==

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

（1）   基礎理論

（2）   商軟操作

（3）   自編程序

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

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

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

自主結構有限元求解器iSolver介紹視頻：

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

==第25篇：顯式分析的穩定時間增量==

相對隱式分析，顯式分析在單元函數中的計算要相對簡單，只需計算應力的更新，而少了單元剛度矩陣的求解，但顯式在增量步的自動步長、網格畸變的控制、質量縮放等方面又比隱式要多做許多工作。其中求解過程中的穩定時間增量是每個顯式分析都會遇到的問題，由于沒有迭代，顯式分析都能得到某個結果，不存在收斂問題，而結果的正確性和穩定時間增量密切相關，很多人做顯式分析都一般直接交給商軟去做自動步長，當商軟顯式結果發散或者結果差異很大時也不清楚怎么修改時間步長。同時，對于顯式自定義單元VUEL等子程序，Abaqus要求用戶自己計算穩定時間增量dtimestable，Abaqus會根據用戶計算的穩定時間增量來限制它的增量步長。本文將簡單介紹一下穩定時間增量的概念和理想及工程應用上的兩種計算方式，并用Abaqus中一個簡單的算例來驗證工程上的穩定時間增量的計算公式，便于你對穩定時間增量的理解和自己編程實現。

                                             

1.1.1 穩定性的含義

在本系列13篇：顯式和隱式的區別中提到，顯式分析都是條件穩定的，譬如下面求解一個微分方程：

                                                                        y'(x) = -y+x+1

其中y(0) = 1。顯式分析可得到下面的增量表達式：

當h<2時，y收斂，但h>=2時，y將發散。當h=2.1時，可發現顯式分析的藍色線將和理論值越來越遠。

上面針對的是一個具體的數學函數，對一個任意的有限元動力學系統，只要采用顯式分析，也同樣存在穩定性問題。動力學問題如果把質量導致的慣性力考慮在內，也是在解下面的平衡方程：

如果上式所有量都是t或者所有量都是t+dt，那么上式肯定是平衡的，沒有一點問題。但在有限元求解過程中，只能通過已知時刻點t來求未知時刻點t+dt的值，而某些值在t+dt的值是不知道的，譬如剛度矩陣K，所以，只能退而求其次，上式中部分是時刻t的，譬如剛度矩陣K，而部分是時刻t+dt的值，譬如載荷F和質量陣M。

這么取值后上式左邊就不再=0了，產生了人為導致的誤差。而這種人為的誤差，如果是隱式分析，那么可以通過迭代來解決，迭代多次后最后的時刻的Tol可以達到預定的小量，方程依然是平衡的。但對顯式分析，流程如下所示，不會再次迭代，只會進行下一次的t+2*dt的求解，這樣，就可能導致兩種結果：

（1）   t時刻的誤差將在t+dt時刻累積，此時系統就是不穩定的。譬如上圖中的藍色曲線。

（2）   t時刻的誤差在t+dt時刻不累積，此時系統就是穩定的。譬如上圖中的如果h<2。

1.1.2 穩定時間增量的理想計算方式

總結顯式分析中增量步之間的關系，可以發現增量步和上一個增量步有簡單的關系：

那么當前時刻和第一個時刻的一個簡單關系：

而這個C類似于一個比例系數，與dt的取值有關，譬如上述函數例子中C=1-dt。對這種只有一個自由度的問題，可以發現|C|<1時，系統穩定，否則系統不穩定。在數學上可以證明C與系統的特征值是密切相關的。也就是說dt與系統的特征值有關，對一個有限元系統來說，數學上已經證明dt如果滿足下式就是穩定的：

其中

為有限元系統的最大特征值。我們稱右端為穩定時間增量的理想計算方式。

1.1.3 穩定時間增量簡化的必要條件

在穩定時間增量dt_ideal的理想計算方式中，需要計算有限元系統的最大特征值。而如果每次顯式分析前都要計算一遍最大特征值，且對幾何變化大的問題，在每個增量步都要更新一次，那么計算量將非常大。更何況對非對稱的剛度等系統，還沒有辦法得到特征值。在實際有限元的工程計算中，會找另一種簡單的穩定時間增量dt_engeer的工程計算方式，這個dt_engeer必須滿足兩個條件：

（1）   條件1：這種計算方式并不一定嚴格求出上面的dt_ideal，但求出的穩定時間一定比理想的穩定時間增量更小。即當有限元的時間增量步長dt<dt_engeer，必然有dt<dt_ideal，此時系統依然滿足穩定性要求。

（2）   條件2：dt_engeer比dt_ideal不能小太多，應該是同一量級，否則增量步將大幅增加，導致系統耗時過大。

1.1.4 一般有限元的穩定時間增量的工程計算方法

本著這個原則，一般有限元上做了兩次計算簡化工作：

（1）第一步：是將整個系統的最大特征值取為所有單個單元的特征值,由于系統的約束會壓縮總體頻率，導致dt_element<dt_ideal。也就是：

（2）單個單元依然有特征值求解問題，因此，再進一步簡化，實際取的穩定時間增量不再需要計算特征值，而是估算為：

其中min表示對所有單元遍歷后的最小值，Le為單元特征長度，Cd為材料疏密波速度，對桿件有Cd=sqrt(E/d), E、d分別為桿長度、楊氏模量和密度。按此公式計算將非常簡單，只需遍歷所有單元，將得到的參數代入上式即可。

為了驗證上述簡化滿足前面所述兩點，我們舉一個由桿單元truss組成的簡單系統來驗證一下dt_engeer的具體值。此時，每個桿單元的頻率可以簡單的估算為

其中L、E、d分別為桿長度、楊氏模量和密度。注意上式我們用了集中質量。此時的第一步穩定時間增量

顯然，dt_engeer和dt_element同一量級，且略小。

由于dt_engeer<dt_ideal，可以發現有限元商軟采用的dt_engeer不一定是最優的，如果你能找到一個時間增量步dt，使得dt_engeer<dt<dt_ideal，那么可以比商軟默認的耗時更少。后面的例子中也證明了這點。

1.1.5 Abaqus對穩定時間增量的工程計算方法的修正

Abaqus對穩定時間增量的工程計算還做了兩處不同的修正：

（1）   考慮Bulk viscosity效應修正工程穩定時間增量,主要是利用Linear Bulk Viscosity系數b1減小dt_engeer。即

b1在Step->Other的Linear Bulk Viscosity可修改。

默認為0.06，此時

（2）   為了避免系統對實數的精度或者截斷誤差，加了一個Tolerance，使得

顯然這個Tolerance是個遠小于1的值，Abaqus取為0.01，

• 此值我們沒發現在Abaqus界面上怎么修改，如果誰能找到，也希望能交流一下。

1.2 Abaqus的實現驗證

我們將在Abaqus中采用一個簡單的顯式分析算例，來驗證兩個問題：

（1）   Abaqus采取的穩定時間增量和上述最后的穩定時間增量的計算公式一致。

（2）   對于某些問題，其實Abaqus采取的穩定時間增量是dt_engeer，相對保守，實際上可以取一個更加接近dt_ideal的值，系統依然是穩定的。

1.2.1 模型例子

模型為一根長1m，截面為4平方cm的圓柱，左端5cm由Steel組成，楊氏模量和密度分別是2e11Pa和7800千克每立方米，剩下部分由碎石組成，楊氏模量和密度分別是4.432e6Pa和1560千克每立方米。右端點固定在墻上，左端點面載荷4e6sin(150t)。我們只研究0-0.01s時間內的位移。

在Abaqus中建模如下，我們簡單將模型劃分為20個單元。采用truss單元。

1.2.2 穩定時間增量的理論值

1.2.2.1 穩定時間增量的理想計算的理論值

理想計算方式需要先計算系統最大模態特征，由于是20個單元，采用truss單元，就相當于只有21個自由度，右端約束后，無約束的自由度為20個，得到的K和M矩陣的秩為20，那么無論用哪種模態計算方法，得到的模態最大為20階。在Abaqus中計算，結果如下，可得20階的模態頻率為30864Hz。

 

1.2.2.2 穩定時間增量的工程計算的理論值

最小的工程穩定時間增量顯然是左端的Steel單元，此時為：

1.2.3 自動步長

在Abaqus中選擇顯式分析，dynamic,explicit，同時設置為自動步長。

運行結束后查看.sta文件，Abaqus會在此文件中在第一個增量前記錄前十個最小的單元穩定時間增量。可發現如下所示，最小的單元穩定時間增量為第一個單元，且值和理論完全一致：

此時總共增量步為1087次，得到的左端位移隨時間的變換曲線如圖：

1.2.4 固定步長

Abaqus中改為固定步長：

取固定步長分別為dt2=1e-5和dt3=1.06e-5，即

                                                             dt_engeer<dt2<dt_ideal<dt3

得到的左端位移隨時間變化曲線如下：

可發現dt3已經發散，而dt2和自動步長基本一致，但dt2只計算了1000個增量步，比自動步長少了8.7%。從而驗證了，對某些問題，為了加速，時間增量可以取一個超過Abaqus默認更加接近dt_ideal的值，系統依然是穩定的。

1.3 聯系方式

如果有任何其它疑問或者項目合作意向，也歡迎聯系我們：

snowwave02 From www.yqgqt.org.cn

email: snowwave02@qq.com

以往的系列文章：

1.4.1 ========第一階段========

第一篇：S4殼單元剛度矩陣研究。介紹Abaqus的S4剛度矩陣在普通厚殼理論上的修正。

http://www.yqgqt.org.cn/content/post/338859

第二篇：S4殼單元質量矩陣研究。介紹Abaqus的S4和Nastran的Quad4單元的質量矩陣。

http://www.yqgqt.org.cn/content/post/343905

第三篇：S4殼單元的剪切自鎖和沙漏控制。介紹Abaqus的S4單元如何來消除剪切自鎖以及S4R如何來抑制沙漏的。

http://www.yqgqt.org.cn/content/post/350865

第四篇：非線性問題的求解。介紹Abaqus在非線性分析中采用的數值計算的求解方法。

http://www.yqgqt.org.cn/content/post/360565

第五篇：單元正確性驗證。介紹有限元單元正確性的驗證方法，通過多個實例比較自研結構求解器程序iSolver與Abaqus的分析結果，從而說明整個正確性驗證的過程和iSolver結果的正確性。

http://www.yqgqt.org.cn/content/post/373743

第六篇：General梁單元的剛度矩陣。介紹梁單元的基礎理論和Abaqus中General梁單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus梁單元完全一致的剛度矩陣。

http://www.yqgqt.org.cn/content/post/403932

第七篇：C3D8六面體單元的剛度矩陣。介紹六面體單元的基礎理論和Abaqus中C3D8R六面體單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus六面體單元完全一致的剛度矩陣。

http://www.yqgqt.org.cn/content/post/430177

第八篇：UMAT用戶子程序開發步驟。介紹基于Fortran和Matlab兩種方式的Abaqus的UMAT的開發步驟，對比發現開發步驟基本相同，同時采用Matlab更加高效和靈活。

http://www.yqgqt.org.cn/content/post/432848

第九篇：編寫線性UMAT Step By Step。介紹基于Matlab線性零基礎，從零開始Step by Step的UMAT的編寫和調試方法，幫助初學者UMAT入門。

http://www.yqgqt.org.cn/content/post/440874

第十篇：耦合約束（Coupling constraints）的研究。介紹Abaqus中耦合約束的原理，并使用兩個簡單算例加以驗證。

http://www.yqgqt.org.cn/content/post/531029

1.4.2 ========第二階段========

第十一篇：自主CAE開發實戰經驗第一階段總結。介紹了iSolver開發以來的階段性總結，從整體角度上介紹一下自主CAE的一些實戰經驗，包括開發時間預估、框架設計、編程語言選擇、測試、未來發展方向等。

http://www.yqgqt.org.cn/content/post/532475

第十二篇：幾何梁單元的剛度矩陣。研究了Abaqus中幾何梁的B31單元的剛度矩陣的求解方式，以L梁為例，介紹General梁用到的面積、慣性矩、扭轉常數等參數在幾何梁中是如何通過幾何形狀求得的，根據這些參數，可以得到和Abaqus完全一致的剛度矩陣，從而對只有幾何梁組成的任意模型一般都能得到Abaqus完全一致的分析結果，并用一個簡單的算例驗證了該想法。

http://www.yqgqt.org.cn/content/post/534362

第十三篇：顯式和隱式的區別。介紹了顯式和隱式的特點，并給出一個數學算例，分別利用前向歐拉和后向歐拉求解，以求直觀表現顯式和隱式在求解過程中的差異，以及增量步長對求解結果的影響。

http://www.yqgqt.org.cn/content/post/537154

第十四篇：殼的應力方向。簡單介紹了一下數學上張量和Abaqus中殼的應力方向，并說明Abaqus這么選取的意義，最后通過自編程序iSolver來驗證殼的應力方向的正確性。

http://www.yqgqt.org.cn/content/post/1189260

第十五篇：殼的剪切應力。介紹了殼單元中實際的和板殼近似理論中的剪切應力，也簡單猜測了一下Abaqus的內部實現流程，最后通過一個算例來驗算Abaqus中的真實的剪切應力。

http://www.yqgqt.org.cn/content/post/1191641

第十六篇：Part、Instance與Assembly。介紹了Part、Instance與Assembly三者之間的關系，分析了Instance的網格形成原理，并猜測Abaqus的內部組裝實現流程，隨后針對某手機整機多part算例，通過自編程序iSolver的結果比對驗證我們的猜想。

http://www.yqgqt.org.cn/content/post/1195061

第十七篇：幾何非線性的物理含義。介紹了幾何非線性的簡單的物理含義，并通過幾何非線性的懸臂梁Abaqus和iSolver的小應變情況的結果，從直觀上理解幾何非線性和線性的差異。

http://www.yqgqt.org.cn/content/post/1198459

第十八篇：幾何非線性的應變。首先從位移、變形和應變的區別說起，然后通過一維的簡單例子具體介紹了幾何非線性下的應變的度量方式，并給出了工程應變、 真實應變、Green應變三者一維情況下在數學上的表達方式。

http://www.yqgqt.org.cn/content/post/1201375

第十九篇：Abaqus幾何非線性的設置和后臺。首先介紹了幾何非線性一般的分類，然后詳細說明了Abaqus中幾何非線性的設置方式和常用單元的分類，最后以一個殼單元的簡單算例為對象，可以發現應變理論、Abaqus和iSolver三者在線性、小應變幾何非線性和大應變幾何非線性三種情況下都完全一致，從而驗證Abaqus幾何非線性后臺采用的應變和我們的預想一致。

http://www.yqgqt.org.cn/content/post/1203064

第二十篇：UEL用戶子程序開發步驟。本文首先簡單的討論了UEL的一般含義，并詳細的介紹了基于Fortran和Matlab兩種方式的UEL的開發步驟，對比發現開發步驟基本相同，但Matlab更加高效和靈活。

http://www.yqgqt.org.cn/content/post/1204261

1.4.3 ========第三階段========

第二十一篇：自主CAE開發實戰經驗第二階段總結。從實戰角度介紹自主CAE在推廣和工程化應用的過程中的體會，同時說明一個CAE平臺最重要的兩個特點：可擴展和易維護。

http://www.yqgqt.org.cn/content/post/1204970

第二十二篇：幾何非線性的剛度矩陣求解。介紹幾何非線性下的剛度矩陣的理論推導和計算機求解方法，最后利用一個簡單的算例驗證我們對Abaqus幾何非線性的剛度矩陣的實現方式的猜測。

http://www.yqgqt.org.cn/content/post/1254435

第二十三篇：編寫簡單面內拉伸問題UEL Step By Step。通過簡單面內拉伸問題UEL的編寫，介紹解決有限元問題從理論到算法再到編程實現的一般流程以及面內拉伸問題的基本算法步驟，最后將UEL與Abaqus的S4R單元計算結果進行對比，進行驗證。

http://www.yqgqt.org.cn/content/post/1256835

第二十四篇：顯式求解Step By Step。本文概要性地介紹了Abaqus中心差分法的理論以及算法實現的整體流程，并通過簡單彈簧顯式動力學分析算例與iSolver和Abaqus計算結果進行對比，驗證了算法和整體流程的正確性。

http://www.yqgqt.org.cn/content/post/1261165