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

1 概述

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

（1） 基礎理論

（2） 商軟操作

（3） 自編程序

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

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

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

百度網盤鏈接: https://pan.baidu.com/s/10d6jHdZ01SBY2JxiS6bffw 提取碼: 6fdf

2 屈曲分析

結構失效的方式有兩種：1. 材料的失效，也就是強度的破壞，2.幾何的失效，也就是結構件的失穩，也叫屈曲。在系列文章4：非線性問題的求解中我們提到了弧長法求非線性屈曲，但非線性屈曲只是屈曲分析的一種方法，在工程上更多是一種簡化方法。我們將分兩個章節介紹結構屈曲分析，本章將介紹屈曲的理論，下一章節介紹屈曲的實驗和模擬。

2.1 屈曲的含義

“屈曲”兩個字中：

屈：大丈夫能屈能伸。屈和伸相反，表示結構件受壓的狀態。

曲：九曲十八彎，曲表示結構件彎曲。

屈曲的含義把“屈”和“曲”連起來就是結構件在受壓時彎曲的現象，更具體點就是結構件在壓力或載荷作用下，因穩定性不足而突然發生大變形甚至失效的現象。用大白話解釋“屈曲”就是：結構件在壓力下“突然塌掉”或“彎成奇怪形狀”的現象。

屈曲一般發生在細長壓桿或者薄板等結構件中。生活中有很多這樣的例子，譬如帳篷的支架在大力下或者頂端放個重包突然失去支撐能力，導致帳篷坍塌，又譬如空的易拉罐用手指按壓時，按壓點會癟下去，力比較小時，易拉罐外殼還能恢復，當指力足夠大時，易拉罐外殼就直接現成一個永久的坑了。

工程上也有很多屈曲的現象，譬如火箭發射過程中，由于火箭加速產生的巨大軸向載荷（如發射時的加速度可達數十倍重力加速度），會導致細長結構（如箭體或儲箱）因穩定性不足而彎曲，又譬如船舶在波浪中航行時，船體因波浪起伏產生巨大的縱向彎曲應力，導致中拱或中垂狀態，引發肋骨或甲板的屈曲。

2.2 屈曲分析的目點

從以上生活常識可知，總有一個突然就垮塌的過程，就像上面的易拉罐例子一樣，如果力比較小，放手后易拉罐還能恢復原樣，我們認為這個力對這個結構件是安全的，但如果放手后產生了不可恢復的形變，那么這個力對結構件是不安全的，從安全到不安全的臨界點的力就是臨界載荷。

屈曲這個過程經歷了：線彈性變形過程-》達到屈曲臨界點-》（后屈曲）屈曲后的結構變形過程。

在工程上屈曲分析的主要目點是計算結構在軸向壓力或彎曲荷載作用下發生屈曲失效的臨界載荷值，從而判斷當前設計是否安全。

2.3 屈曲分析的方法

屈曲分析有多種方法：

2.3.1 非線性屈曲分析

非線性屈曲分析是將力隨著位移的關系表達出來，直到能看出哪點是臨界載荷，臨界載荷時位移增加時，力將不再增加，反而下降，也就是臨界載荷就是載荷Vs位移曲線上的馬鞍點位置。

如下面例子：

兩個桿下方固定位移，中間往下加壓力，每次力加的特別緩慢，確保每個時刻力達到平衡狀態。

如果以往下的位移為橫坐標，加的力為縱坐標，那么畫出一條曲線大體如下：

力一開始隨著位移增加而增加，知道頂點A，當過了頂點再往下壓時，生活常識告訴我們不需要那么大的力也能往下壓了，此時力隨位移減小直到在水平位置的力變為0。采用有限元的非線性屈曲分析就是要尋找上述過程中的A點時的壓力，有兩種方式：

（1） 加強制位移約束，然后輸出反力，做出反力和位移的曲線圖，直接在圖上查看馬鞍點位置對應的反力大小。

（2） 上面方法缺點是實際工程我們不知道應該加哪個強制約束點或者如殼的屈曲，強制約束點太多無妨加和實際一樣的強制位移，所以實際工程更多的是加力，做出力和位移的關系，此時由于要越過馬鞍點，一般在有限元中采用弧長法（Risk）替代Newton迭代來計算，具體可看系列文章4：非線性問題的求解。

2.3.2 基于特征值的線性屈曲

非線性屈曲分析可捕捉整個屈曲從發生到后屈曲整個過程，同時如果變形量比較大或者材料達到塑性段時，那必須考慮非線性問題，只能采用非線性屈曲分析。但只要是非線性分析都有共同缺點，需要更大的計算量，而且非線性屈曲存在馬鞍點時收斂更困難。

在實際工程上，很多問題不關心后屈曲的情況，只關心屈曲剛發生時刻點的臨界載荷，而且有相當一大部分屈曲發生時材料還沒達到塑性，變形也沒超過5%，因此此時可以采用線性屈曲分析更快的得到臨界載荷。

線性屈曲最通用的數值計算方法是基于特征值來求解。核心思想是如果我們已知時刻0和1的兩個位移和剛度K0、K1，那么能不能得到t時刻的剛度呢？

結構有限元的剛度陣按照虛功原理得到。

虛功原理可以理解為外力在虛擬位移下做的虛功=內部應變能的一段小時間內對應變能的積分：

S和E分別表示應力和應變。

求導后得到任意時刻的剛度陣K如下：

也就是剛度矩陣將分為兩塊，上式的前面一部分依然是以前的BDB形式，只不過B換成了當前時刻的應變位移矩陣，而后面新增項將將轉換為S*G剛度陣，稱為幾何剛度陣，也稱為初始應力矩陣（initial stress stiffness）。

因為是小變形線性材料，所以B、D、G為常量，BDB我們可以認為是一個常量，也就是K0，則

? 因為壓縮導致剛度減小，我們在上式中假定S和G是正數。

我們任意時刻外力為F，由于線性屈曲假定整個屈曲前過程都是線性問題，那么得到的應力S正比于位移U，也正比于F，設為

帶入上式得到

(1)

上式就是一個一次方程。將K和F的關系是一條直線如下

其中，這條直線與y軸相交于已知點（0，K0），與x軸相交于未知點（Fe，0）。Fe即我們要求的臨界載荷。

下面僅是怎么求出斜率AG，一種簡單的方式就是在直線上找兩個已知點就能求出斜率了。既然已經有一個已知點（0，K0），那么取時刻1作為另一個已知點（F1，K1）

2.3.3 基于歐拉應力理論修正的線性屈曲

非線性屈曲分析和基于特征值的線性屈曲看起來已經把有限元屈曲分析的所有情況覆蓋了，但實際工程上很多行業還是采用基于歐拉應力理論的線性屈曲。

針對一根壓桿的受壓，如下圖，左端簡支，右端約束y、z位移且加壓力F。設此壓桿是完全彈性的，且應力不超過比例極限，若軸向外載荷F小于它的臨界值Fe，此桿將保持直的狀態而只承受軸向壓縮。如果一個擾動(如—橫向力)作用于桿，使其有一小的撓曲，在這一擾動除去后。撓度就消失，桿又恢復到平橫狀態，此時桿的直的形式的彈性平衡是穩定的。若軸向外載荷F大于它的臨界值Fe，柱的直的平衡狀態變為不穩定，即任意擾動產生的撓曲在擾動除去后不僅不消失，而且還將繼續擴大，直至達到遠離直立狀態的新的平衡位置為止，或者彎折。此時，此壓桿失穩或屈曲。

對此簡單問題，我們猜測一下這個Fe大概為多少？

當到達Fe時，壓桿開始便變形，根據生活常識，應該大體變形為如下形狀：

顯然當L足夠小時，一定會超過材料屈服強度也會到時結構件失效。

實際工程材料因此如果將結構件失效應力和長度做一條曲線將會是如下形式

這條曲線在L>Ly時是雙曲線，在L<Ly時是直線，且失效應力恒定為材料屈服強度。基于特征值的線性屈曲分析結合屈服分析就會得到這條曲線。

這個簡單例子用實驗很容易得到，下章我們也將做一個簡單的實驗（估計應該是全網第一個采用簡單的家用工具做出來的屈曲試驗）。做實驗發現是如下修正曲線：

實際沒有一個明顯的轉折點，很多行業規范會人為分為兩段，在L>Lp時還是雙曲線，而在L<Lp時采用拋物線，兩者的相交點為(Lp, )。在船舶行業，取為/2。

這是一根壓桿得到的曲線，模擬的最終目點還是和實驗盡量接近，既然它比基于特征值的線性屈曲分析更接近試驗，那么在實際工程中也更受歡迎。船舶行業的線性屈曲就采用基于歐拉應力理論修正的線性屈曲。長方形殼單元可以看成是壓桿截面的一個維度取為實際平面尺寸的一個應用。同時，為了適用一般的殼形狀，船舶行業的規范規定了三步的模擬：

（1） 先確定板格的位置，周圍由桁材、縱骨或者不在一個平面的面板圍出來的圖形就是板格，如果是有限元模型，板格一般由多個板單元組成。如下便是一個板格的示例，每個不同顏色表示不同的板格：

（2） 板格等效為長方形的規范，上面示例的板格都是長方形，但實際問題很多板格都是不規則形狀，如下分別就是一個帶弧度的殼單元和三角形等效為長方形的規范示例。

（3） 套用基于歐拉應力理論修正的線性屈曲公式分析是否屈曲。

2.4 屈曲的本質

以下面簡支梁為例：

模型尺寸：

長L=240，截面為10X5。E=1.5e6，打開幾何非線性模擬實際情況。

（1）當不加任何力時，打印剛度陣如下，y方向K22=230：

（2）由上面歐拉公式得到歐拉屈曲臨界力Fe=26773，一端施加壓力Fe=26773，

考慮梁的幾何非線性，得到的剛度陣如下，可以看到K11增加，但K22，K33等減少。K22=118.

（3）壓力增大為Te=26773*2，得到剛度陣如下，K22繼續減少為7.3，已經非常接近0了。

由這個例子，可以知道，桿件受加壓載荷時，橫向方向剛度隨壓力反比減小，此時如果有一個小的橫向力，將產生大的位移，導致橫向抗外力能力的下降，導致失穩彎折。這就是屈曲的本質。

? 幾何非線性時，加了線性屈曲歐拉理論的壓縮載荷，Y方向剛度減小了，甚至2倍時也不為0，所以線性屈曲的理論和實際還是有一些差異的，有興趣可以再細致研究一下差異的來源（個人猜測需要研究考慮幾何非線性的B導致的幾何剛度陣的差異）。

2.5 總結

本文總結了屈曲分析的三種方法：

（1） 非線性屈曲分析

（2） 基于特征值的線性屈曲

（3） 基于歐拉應力理論修正的線性屈曲

在實際工程中，線性屈曲更適合于工程問題，同時，由于線性屈曲理論和試驗的差異，部分行業實際上更多采用基于歐拉應力理論修正的線性屈曲。同時說明屈曲的本質還是縱向加壓后橫向剛度變小，導致橫向抗力能力的下降，導致失穩彎折。

2.6 視頻講解和操作驗證演示

如果覺得上面的文字太復雜，也可以看一下視頻的簡要講解和軟件演示，地址如下：

http://www.yqgqt.org.cn/college/video/c12884 20理論系列文章48-屈曲分析(1)-理論

3 以往的系列文章

3.1 ========第一階段========

第一篇：S4殼單元剛度矩陣研究。

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

第二篇：S4殼單元質量矩陣研究。

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

第三篇：S4殼單元的剪切自鎖和沙漏控制。

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

第四篇：非線性問題的求解。

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

第五篇：單元正確性驗證。

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

第六篇：General梁單元的剛度矩陣。

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

第七篇：C3D8六面體單元的剛度矩陣。

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

第八篇：UMAT用戶子程序開發步驟。

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

第九篇：編寫線性UMAT Step By Step。

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

第十篇：耦合約束（Coupling constraints）的研究。

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

3.2 ========第二階段========

第十一篇：自主CAE開發實戰經驗第一階段總結。

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

第十二篇：幾何梁單元的剛度矩陣。

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

第十三篇：顯式和隱式的區別。

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

第十四篇：殼的應力方向。

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

第十五篇：殼的剪切應力。

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

第十六篇：Part、Instance與Assembly。

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

第十七篇：幾何非線性的物理含義。

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

第十八篇：幾何非線性的應變。

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

第十九篇：Abaqus幾何非線性的設置和后臺。

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

第二十篇：UEL用戶子程序開發步驟。

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

3.3 ========第三階段========

第二十一篇：自主CAE開發實戰經驗第二階段總結。

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

第二十二篇：幾何非線性的剛度矩陣求解。

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

第二十三篇：編寫簡單面內拉伸問題UEL Step By Step

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

第二十四篇：顯式求解Step By Step。

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

第二十五篇：顯式分析的穩定時間增量。

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

第二十六篇：編寫線性VUMAT Step By Step。

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

第二十七篇：Abaqus內部計算和顯示的應變。

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

第二十八篇：幾何非線性的T.L.和U.L.描述方法

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

第二十九篇：幾何非線性的T.L.和U.L.轉換關系

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

第三十篇：諧響應分析原理

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

3.4 ========第四階段========

第三十一篇：自主CAE開發實戰經驗第三階段總結

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

第三十二篇：諧響應分析算法

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

第三十三篇：線性瞬態動力學

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

第三十四篇：非線性瞬態分析

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

第三十五篇：接觸求解算法

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

 第三十六篇：DLOAD用戶子程序開發步驟

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

 第三十七篇：梁單元差異（1）-理論基礎

https://jishulink.com/content/post/1872208

 第三十八篇：梁單元差異（2）-梁截面方向

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

 第三十九篇：梁單元差異（3）-剪力和彎矩

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

 第四十篇：梁單元差異（4）-形心、剪心和偏置

http://www.yqgqt.org.cn/post/1888000

3.5 ========第五階段========

 第四十一篇：自主CAE開發實戰經驗第四階段總結

http://www.yqgqt.org.cn/post/1905963

 第四十二篇：聲學分析(1)-有限元

http://www.yqgqt.org.cn/post/1912491

 第四十三篇：聲學分析(2)-邊界元

http://www.yqgqt.org.cn/post/1923936

 第四十四篇：聲學分析(3)-濕模態

http://www.yqgqt.org.cn/post/1928692

 第四十五篇：約束關系（1）-統一形式

http://www.yqgqt.org.cn/post/1933077

 第四十六篇：約束關系（2）-Lagrange因子法求解

http://www.yqgqt.org.cn/post/1945262

 第四十七篇：約束關系（3）-船舶規范約束導致的Max Ratio問題

http://www.yqgqt.org.cn/post/1950396