（原創(chuàng)，轉(zhuǎn)載請注明出處）

1 概述

本系列文章研究成熟的有限元理論基礎(chǔ)及在商用有限元軟件的實(shí)現(xiàn)方式，通過

（1）   基礎(chǔ)理論

（2）   商軟操作

（3）   自編程序

三者結(jié)合的方式將復(fù)雜繁瑣的結(jié)構(gòu)有限元理論通過簡單直觀的方式展現(xiàn)出來，同時(shí)深層次的學(xué)習(xí)有限元理論和商業(yè)軟件的內(nèi)部實(shí)現(xiàn)原理。

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

                                            

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

通用結(jié)構(gòu)有限元軟件iSolver介紹視頻：

        

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

==第34篇：非線性瞬態(tài)動(dòng)力學(xué)==

在系列文章第三十三篇時(shí)，我們講了線性瞬態(tài)動(dòng)力學(xué)的原理，瞬態(tài)分析有線性和非線性之分，如果系統(tǒng)有材料、大變形、邊界等非線性效應(yīng)，那么就是非線性瞬態(tài)分析，而瞬態(tài)分析往往都有物體的大轉(zhuǎn)動(dòng)大變形問題，也會(huì)涉及材料的損傷破壞、物體的撞擊接觸等，譬如沖擊爆炸就是典型的強(qiáng)非線性瞬態(tài)動(dòng)力學(xué)的過程，所以在實(shí)際工程中的瞬態(tài)分析都是以非線性為主。本章我們介紹一下非線性瞬態(tài)動(dòng)力學(xué)的求解公式，并以上次線性瞬態(tài)動(dòng)力學(xué)中的單擺例子來說明非線性瞬態(tài)動(dòng)力學(xué)在有限元軟件中的內(nèi)部實(shí)現(xiàn)原理。

                                             

1.1 非線性瞬態(tài)動(dòng)力學(xué)理論

1.1.1 非線性瞬態(tài)動(dòng)力學(xué)方程

瞬態(tài)動(dòng)力學(xué)的運(yùn)動(dòng)方程包括了質(zhì)量陣M相關(guān)的慣性力項(xiàng)和阻尼陣C相關(guān)的阻尼力項(xiàng)，如下：

一般質(zhì)量矩陣與時(shí)間和位移無關(guān)，而K和C如果和位移和時(shí)間相關(guān)，那么就是非線性瞬態(tài)問題。為簡單起見，我們不考慮阻尼項(xiàng)。得到如下表達(dá)式：

和線性瞬態(tài)動(dòng)力學(xué)類似，在理論上當(dāng)前時(shí)刻的加速度就應(yīng)該只和當(dāng)前時(shí)刻的位移相關(guān)了，但計(jì)算的數(shù)值解把時(shí)刻分成了多個(gè)，方程如下：

其中

可以取當(dāng)前時(shí)間步的初始時(shí)刻的剛度，那么就是修正的牛頓迭代方法，如果取結(jié)束時(shí)刻的剛度，就是完全牛頓迭代方法。

而加速度一定是要通過起碼兩個(gè)時(shí)刻才能計(jì)算出來的，取前時(shí)刻還是后時(shí)刻的位移來決定加速度就有隱式和顯式之分，任何一種分析理論上都可以得到正確解，只不過顯式不用求剛度，隱式需要求解剛度陣，而通過剛度的解釋可以更容易理解Abaqus內(nèi)部的實(shí)現(xiàn)原理，在本章中，我們選用隱式非線性瞬態(tài)分析。在隱式方法中，將時(shí)間離散，加速度表示為位移的表達(dá)式后，得到下方的公式

這個(gè)就是標(biāo)準(zhǔn)的newton增量迭代法的公式，其中K和F在靜力學(xué)的增量迭代法中表示切線剛度陣和內(nèi)力，只不過這邊需要加上M項(xiàng)：

其中A依然是只與時(shí)間增量步相關(guān)的量，對不同的隱式算法A的值不同，譬如對最簡單的Newmark方法中的梯形隱式算法，加速度在增量步內(nèi)線性變化，譬如下方紅線是在輸入的載荷轉(zhuǎn)換為加速度的值，但計(jì)算中沒法處理連續(xù)的曲線，所以Abaqus或者iSolver中實(shí)際上只會(huì)認(rèn)為增量步內(nèi)加速度時(shí)線性變化的：

此時(shí)得到的修正剛度為：

而F與上一迭代的位移、速度、加速度相關(guān)，可以表示為：

顯然當(dāng)時(shí)間增量步固定時(shí)，此時(shí)的K也是隨增量步和迭代步變化的值，顯然該方程在每個(gè)時(shí)間步長內(nèi)都只能通過迭代。因?yàn)橛械訤也無法約去右側(cè)純的應(yīng)力導(dǎo)致的內(nèi)力項(xiàng)。

1.1.2 穩(wěn)定時(shí)間步

最后我們再討論一下穩(wěn)定時(shí)間步的問題，在系列文章13：顯式和隱式的區(qū)別和25：顯式分析的穩(wěn)定時(shí)間增量中，我們都討論過，CAE求解方法一般有兩種，分別為顯式（Explicit）和隱式（Implicit）。顯式分析時(shí)間步長不穩(wěn)定的，就是增量步長也不能過大，一般不超過模型最小自由振蕩周期的1/10，否則容易導(dǎo)致計(jì)算結(jié)果發(fā)散。

那么隱式就一定是時(shí)間步長穩(wěn)定的嗎？如果是穩(wěn)定的，那為何大家用Abaqus的Standard隱式求解器的時(shí)候處理后屈曲、材料軟化或者沖擊爆炸這種問題經(jīng)常會(huì)碰到很難收斂的問題呢？

大家使用軟件的經(jīng)驗(yàn)只能說明隱式算法也不一定穩(wěn)定的。Abaqus的Standard求解器都是把有限元的問題最終轉(zhuǎn)化為增量迭代法求解數(shù)值問題的，時(shí)間步長是否穩(wěn)定就看增量迭代法求解的過程。

（1）對于線性瞬態(tài)動(dòng)力學(xué)，無需迭代，一次迭代就能求出問題解，且是正確解，自然是穩(wěn)定的。

（2）對于非線性瞬態(tài)動(dòng)力學(xué)，需要多次迭代完成，只要迭代求解曲線問題的算法，就可能存在無法收斂的問題，譬如坍塌問題，也就是力一開始隨著位移上升，在某個(gè)極值點(diǎn)時(shí)力隨著位移下降，曲線有下面類似的馬鞍點(diǎn)。譬如下面的A點(diǎn)，在接近接近極值點(diǎn)A時(shí)，此時(shí)斜率非常小設(shè)為，如果第一個(gè)t=0時(shí)刻的K是K0的話，此時(shí)斜率可以認(rèn)為a*K0，a非常小。

 

1.2 非線性瞬態(tài)動(dòng)力學(xué)的單擺模型

1.2.1 模型介紹

本章我們依然采用線性瞬態(tài)動(dòng)力學(xué)中的單擺例子

這次我們在Abaqus中將幾何非線性打開，此時(shí)軟件將按非線性瞬態(tài)動(dòng)力學(xué)進(jìn)行分析：

1.2.2 計(jì)算結(jié)果

0.5s時(shí)刻的位移（左側(cè)是Abaqus，右側(cè)是iSolver）

1s時(shí)刻的位移（左側(cè)是Abaqus，右側(cè)是iSolver）

1.5s時(shí)刻的位移（左側(cè)是Abaqus，右側(cè)是iSolver）

4.12s時(shí)刻的位移（左側(cè)是Abaqus，右側(cè)是iSolver）

整個(gè)運(yùn)動(dòng)過程動(dòng)畫如下（左側(cè)是Abaqus，右側(cè)是iSolver）

 

查看右端節(jié)點(diǎn)位移曲線如下：

 

顯然，單擺隨著時(shí)間在紙面內(nèi)左右擺動(dòng)，和真實(shí)情況一致，且可以發(fā)現(xiàn)一個(gè)周期和理論值4.12s吻合。

 

1.2.3 結(jié)果解析

1.2.4 迭代步次數(shù)

在系列文章32章線性瞬態(tài)動(dòng)力學(xué)中，我們可以發(fā)現(xiàn)如果幾何非線性沒打開，單擺的迭代iter在Monitor中都是1，也就是沒有迭代：

本章中，將幾何非線性打開后，可以發(fā)現(xiàn)Monitor中，iter已經(jīng)不再是1了，說明了非線性瞬態(tài)動(dòng)力學(xué)需要迭代才能執(zhí)行。

1.2.4.1 第一個(gè)增量步第一次迭代

K在第一個(gè)增量步第一個(gè)迭代步時(shí)，是水平方向初始長度的剛度：

在第一個(gè)增量步，方程就變?yōu)椋?/p>

初始條件下：

（1）R只有Y方向恒值

（2）初始位置時(shí)位移、速度都是0，但Y方向存在加速度

（3）

顯然，上述三個(gè)方程變?yōu)椋?/p>

U1、U2、U3沒有耦合項(xiàng)，可以直接解耦，分別得到三個(gè)方程。

由第一個(gè)和第三個(gè)方程可得位移U1和U3在第一個(gè)增量步結(jié)束時(shí)=0。在程序中強(qiáng)制設(shè)置迭代步=1，在0.1s時(shí)刻可得U1=0，和上述理論一致：

第二個(gè)方程可得在第一個(gè)增量步結(jié)束時(shí)的位移是R/(A*M)，將Newmark的梯形算法的A帶入后得到：

這個(gè)公式顯然和我們熟知的位移和加速度的關(guān)系式完全一致。

在0.1s，得到位移U2=-4.9，在程序中強(qiáng)制設(shè)置迭代步=1，在0.1s時(shí)刻可得U1=0，和上述理論一致：

1.2.4.2 第一次迭代步收斂判斷

在第一次迭代后，程序?qū)⑴袛嗳齻€(gè)方向的不平衡力情況

此時(shí)的F要用此刻的位移、速度、加速度等構(gòu)型計(jì)算，在這個(gè)時(shí)刻，由于單擺已經(jīng)擺動(dòng)了一個(gè)角度，因此，F(xiàn)計(jì)算出的內(nèi)力必然和KU已經(jīng)不一致了，導(dǎo)致也和外力R不一致，只能進(jìn)入下一次迭代。

1.2.4.3 穩(wěn)定時(shí)間檢驗(yàn)

在前面理論中，我們也提到了，只要是迭代，就有不收斂的可能性，在該例中，我們將步長設(shè)為0.2s，此時(shí)得到的位移曲線如下：

位移看起來很正常，但我們要是查看桿中的應(yīng)變曲線時(shí)，可以發(fā)現(xiàn)應(yīng)變在半個(gè)周期內(nèi)還是正常的，但在半個(gè)周期后面明顯發(fā)散了。

1.3 視頻講解和操作驗(yàn)證演示

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

http://www.yqgqt.org.cn/college/video/c12884 20理論系列文章34-非線性瞬態(tài)動(dòng)力學(xué)理論.mp4

以往的系列文章：

1.6.1 ========第一階段========

第一篇：S4殼單元?jiǎng)偠染仃囇芯?/strong>。

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

第二篇：S4殼單元質(zhì)量矩陣研究。

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

第五篇：單元正確性驗(yàn)證。

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用戶子程序開發(fā)步驟。

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

1.6.2 ========第二階段========

第十一篇：自主CAE開發(fā)實(shí)戰(zhàn)經(jīng)驗(yàn)第一階段總結(jié)。

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

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

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

第十三篇：顯式和隱式的區(qū)別。

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

第十四篇：殼的應(yīng)力方向。

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

第十五篇：殼的剪切應(yīng)力。

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

第十八篇：幾何非線性的應(yīng)變。

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

第十九篇：Abaqus幾何非線性的設(shè)置和后臺(tái)。

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

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

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

1.6.3 ========第三階段========

第二十一篇：自主CAE開發(fā)實(shí)戰(zhàn)經(jīng)驗(yàn)第二階段總結(jié)。

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

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

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

第二十三篇：編寫簡單面內(nèi)拉伸問題UEL Step By  Step

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

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

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

第二十五篇：顯式分析的穩(wěn)定時(shí)間增量。

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

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

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

第二十七篇：Abaqus內(nèi)部計(jì)算和顯示的應(yīng)變。

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

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

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

第二十九篇：幾何非線性的T.L.和U.L.轉(zhuǎn)換關(guān)系

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

第三十篇：諧響應(yīng)分析原理

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

第三十一篇：自主CAE開發(fā)實(shí)戰(zhàn)經(jīng)驗(yàn)第三階段總結(jié)

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

第三十二篇：諧響應(yīng)分析算法

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

第三十三篇：線性瞬態(tài)動(dòng)力學(xué)

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