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

==概述==

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

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

iSolver介紹視頻：

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

==第15篇：殼的剪切應(yīng)力 ==

自編有限元應(yīng)力的校核除了Mises等合力外，也應(yīng)該校核各個(gè)應(yīng)力分量。材料力學(xué)中六個(gè)應(yīng)力分量如下：                  

其中Tau11，Tau22，Tau33為正應(yīng)力，Tau12,13,23為三個(gè)剪切應(yīng)力，對(duì)殼來(lái)說(shuō)，Tau33=0，Tau12為面內(nèi)剪應(yīng)力，Tau13,23即為本文所說(shuō)的橫向剪切應(yīng)力。

最近在做iSolver殼的應(yīng)力分量和Abaqus比對(duì)時(shí)，發(fā)現(xiàn)Abaqus的橫向剪切應(yīng)力和預(yù)想的不一致。iSolver按照常用的殼的理論得到的剪切應(yīng)力是個(gè)與厚度無(wú)關(guān)的常量，但Abaqus的橫向剪切應(yīng)力分量TSHR13，TSHR23，在各個(gè)截面方向積分點(diǎn)section point不一樣。

花了點(diǎn)時(shí)間細(xì)致的研究了一下，猜測(cè)Abaqus中剪切應(yīng)力TSHR13、23是真實(shí)應(yīng)力，但有限元理論和iSolver中計(jì)算的是板殼近似理論中平均剪切應(yīng)力。本章將介紹殼單元中實(shí)際的和板殼近似理論中的剪切應(yīng)力，也猜測(cè)了Abaqus的內(nèi)部實(shí)現(xiàn)流程，最后通過(guò)一個(gè)算例來(lái)驗(yàn)算Abaqus中的真實(shí)的剪切應(yīng)力，并通過(guò)iSolver來(lái)計(jì)算板殼理論的平均剪切應(yīng)力。

1.1 殼的真實(shí)的剪切應(yīng)力

剪應(yīng)力是材料由于抗拒面之間的滑動(dòng)而產(chǎn)生的沿表面方向的應(yīng)力。殼的中間層存在剪切應(yīng)力，這個(gè)可以通過(guò)下面簡(jiǎn)單的例子驗(yàn)證。兩塊板疊加在一起，簡(jiǎn)支，中點(diǎn)加力，板間假定無(wú)摩擦，那么將會(huì)得到下面的形狀，中間層表面上梁的伸長(zhǎng)和下梁的縮短完全由x方向應(yīng)力決定，此時(shí)中間層無(wú)抗拒滑動(dòng)的力，也就不存在剪應(yīng)力，。

但如果兩塊板中面部分用膠水粘住，膠水將會(huì)阻礙中間層上下兩個(gè)面的相對(duì)滑動(dòng)，上邊面的纖維長(zhǎng)度會(huì)變少，下邊面的纖維長(zhǎng)度增加，使得中間層上下兩個(gè)纖維長(zhǎng)度相等，也就在中間層將產(chǎn)生剪應(yīng)力。

同時(shí)，很顯然，上圖最上層面任意一點(diǎn)沒(méi)有任何切向的外力，所以不會(huì)有阻礙滑動(dòng)的剪切應(yīng)力了。而從中間層到最上層，可以猜測(cè)剪應(yīng)力將逐步減小。根據(jù)材料力學(xué)的理論，實(shí)際的截面上的剪應(yīng)力分布如下：

其中V為剪力。

顯然與材料點(diǎn)所在截面的y方向坐標(biāo)y1是二次關(guān)系。用圖形化表示為下圖右側(cè)，隨截面厚度方向是拋物線，中面最大，上下表面為0：

在各向同性材料中，剪切應(yīng)力和剪切應(yīng)變也就是剪切角成正比，所以，如果一個(gè)橡皮條做成的懸臂梁，那么原來(lái)在鉛直線上畫的直線受力后將會(huì)變成如下圖所示的m1n1的曲線。

1.2 板殼近似理論的平均剪切應(yīng)力

上面的剪切應(yīng)力表達(dá)式中，需要預(yù)先知道剪力V，但實(shí)際有限元流程中在計(jì)算剛度矩陣K時(shí)并不知道V，K只由兩者決定，一個(gè)是應(yīng)力應(yīng)變關(guān)系，也就是本構(gòu)關(guān)系矩陣C，另一個(gè)是應(yīng)變和位移關(guān)系，也就是常說(shuō)的B矩陣。因此有限元中對(duì)殼做了直線法的近似，認(rèn)為變形前垂直于中面的截面的所有材料點(diǎn)變性后依然位于一個(gè)平面內(nèi)，譬如下圖的紅色箭頭表示原先面上的所有材料點(diǎn)受力后組成新的材料點(diǎn)平面。根據(jù)這個(gè)假設(shè)，那么可以發(fā)現(xiàn)所有的剪切角也就是剪切應(yīng)變是個(gè)恒定值，乘以各向同性的剪切模量G，那么得到的剪切應(yīng)力也是恒定值，相當(dāng)于一個(gè)平均效應(yīng)的應(yīng)力，但后面的例子通過(guò)iSolver的計(jì)算可以看到，這個(gè)平均剪切應(yīng)力并不是簡(jiǎn)單的是剪力V在截面上的平均。

1.3 Abaqus內(nèi)部實(shí)現(xiàn)流程猜測(cè)

Abaqus程序的內(nèi)部流程由增量迭代法實(shí)現(xiàn)，猜測(cè)和普通的有限元理論是一致的，對(duì)靜力分析步驟如下：

1. 根據(jù)本構(gòu)關(guān)系和尺寸得到K。

2. 由外力平衡得到位移d。

3. 由位移d計(jì)算內(nèi)部剪切應(yīng)力S13、S23，從而得到節(jié)點(diǎn)力。

4. 求非平衡力，判斷收斂，如果收斂，那么結(jié)束。

而S13、S23猜測(cè)應(yīng)該是按板殼近似理論得到的與厚度無(wú)關(guān)的平均應(yīng)力，真實(shí)的剪切應(yīng)力并不參與迭代。

在迭代完畢后，如果用戶需要輸出截面真實(shí)的剪切應(yīng)力，那么根據(jù)前面所說(shuō)的剪力V來(lái)計(jì)算TSHR13/23。

1.4 算例

我們只能驗(yàn)算Abaqus真實(shí)剪切應(yīng)力TSHR13和23的結(jié)果，但沒(méi)有找到方法來(lái)驗(yàn)證用于Abaqus流程的剪切應(yīng)力是取的平均應(yīng)力還是真實(shí)應(yīng)力，因?yàn)锳baqus的S13，S23沒(méi)有找到方法輸出。但通過(guò)自編程序iSolver，我們還是計(jì)算了橫向平均剪切應(yīng)力的大小，如果有人也自己編程序，遇到類似問(wèn)題時(shí)可以對(duì)比一下結(jié)果。

1.4.1 算例說(shuō)明

只取一個(gè)簡(jiǎn)單的長(zhǎng)方形。

參數(shù)如下：

尺寸：5X1，厚度0.1。

材料：Young’s Modulus 1e8， Poisson Ratio 0.3。

右側(cè)兩個(gè)節(jié)點(diǎn)固支。

左側(cè)兩個(gè)節(jié)點(diǎn)每個(gè)加集中力1e5，z方向。這樣將產(chǎn)生面外彎曲。

在Step中勾選輸出截面應(yīng)力和輸出截面積分點(diǎn)的變量。

1.4.2 中面真實(shí)剪切應(yīng)力理論結(jié)果

V=100000*2，b=1，h=0.1

TauMax=3/2*V/(b*h)=3e6

1.4.3 Abaqus真實(shí)剪切應(yīng)力結(jié)果

下表面TSHR13=0：

中面為3e6，和理論完全一致。

1.4.4 平均剪切應(yīng)力結(jié)果

Abaqus中無(wú)法輸出S13和S23，也就是平均剪切應(yīng)力，為了便于大家理解S13這個(gè)值大概是什么量級(jí)，我們采用iSolver，在內(nèi)部計(jì)算得到積分點(diǎn)上的剪切應(yīng)力得到S13=8.1e5，而只有采用平均面積得到的V/(b*h)=2e6的2/5左右，就算加上剪切修正因子5/6也不是簡(jiǎn)單的按平均面積得到的結(jié)果，猜測(cè)是因?yàn)樵谶@種情況下，殼只在兩個(gè)節(jié)點(diǎn)處加載荷，而不是在一個(gè)自由邊上均勻加載，而平均面積計(jì)算的剪切應(yīng)力是假定均勻加載得到，所以不一致。當(dāng)然，有限元中采用這個(gè)積分點(diǎn)應(yīng)力來(lái)求節(jié)點(diǎn)力，得到的節(jié)點(diǎn)力依然和外力是平衡的。

==總結(jié)==

本章介紹了殼單元中實(shí)際的和板殼近似理論中的剪切應(yīng)力，也簡(jiǎn)單猜測(cè)了一下Abaqus的內(nèi)部實(shí)現(xiàn)流程，最后通過(guò)一個(gè)算例來(lái)驗(yàn)算Abaqus中的真實(shí)的剪切應(yīng)力。同時(shí)還有兩個(gè)疑問(wèn)也希望與大家討論，得到大家的指點(diǎn)。

（1）我們暫時(shí)沒(méi)有找到方法來(lái)驗(yàn)證用于Abaqus流程的剪切應(yīng)力是取的平均應(yīng)力還是真實(shí)應(yīng)力，因?yàn)锳baqus的S13，S23沒(méi)有找到方法輸出，不知道誰(shuí)了解怎么輸出殼的這兩個(gè)應(yīng)力？

（2）按照下面的公式，Mises力應(yīng)該包括正應(yīng)力和剪切應(yīng)力，但查看殼的Mises應(yīng)力結(jié)果可以看出，殼的Mises力沒(méi)有計(jì)入Tau13，Tau23，不知道為什么這樣取Mises力？這樣用來(lái)校核殼的Mises應(yīng)力達(dá)到屈服是否會(huì)有問(wèn)題？

如果有任何其它疑問(wèn)或者項(xiàng)目合作意向，也歡迎聯(lián)系我們：

snowwave02 From www.yqgqt.org.cn

email: snowwave02@qq.com

 

以往的系列文章：

第一篇：S4殼單元?jiǎng)偠染仃囇芯?/strong>。介紹Abaqus的S4剛度矩陣在普通厚殼理論上的修正。

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

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

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

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

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

第四篇：非線性問(wèn)題的求解。介紹Abaqus在非線性分析中采用的數(shù)值計(jì)算的求解方法。

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

第五篇：單元正確性驗(yàn)證。介紹有限元單元正確性的驗(yàn)證方法，通過(guò)多個(gè)實(shí)例比較自研結(jié)構(gòu)求解器程序iSolver與Abaqus的分析結(jié)果，從而說(shuō)明整個(gè)正確性驗(yàn)證的過(guò)程和iSolver結(jié)果的正確性。

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

第六篇：General梁?jiǎn)卧膭偠染仃?/strong>。介紹梁?jiǎn)卧幕A(chǔ)理論和Abaqus中General梁?jiǎn)卧膭偠染仃嚨男拚绞剑捎眠@些修正方式可以得到和Abaqus梁?jiǎn)卧耆恢碌膭偠染仃嚒?/p>

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

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

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

第八篇：UMAT用戶子程序開發(fā)步驟。介紹基于Fortran和Matlab兩種方式的Abaqus的UMAT的開發(fā)步驟，對(duì)比發(fā)現(xiàn)開發(fā)步驟基本相同，同時(shí)采用Matlab更加高效和靈活。

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

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

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

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

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

第十一篇：自主CAE開發(fā)實(shí)戰(zhàn)經(jīng)驗(yàn)第一階段總結(jié)。介紹了iSolver開發(fā)以來(lái)的階段性總結(jié)，從整體角度上介紹一下自主CAE的一些實(shí)戰(zhàn)經(jīng)驗(yàn)，包括開發(fā)時(shí)間預(yù)估、框架設(shè)計(jì)、編程語(yǔ)言選擇、測(cè)試、未來(lái)發(fā)展方向等。

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

第十二篇：幾何梁?jiǎn)卧膭偠染仃?/strong>。研究了Abaqus中幾何梁的B31單元的剛度矩陣的求解方式，以L梁為例，介紹General梁用到的面積、慣性矩、扭轉(zhuǎn)常數(shù)等參數(shù)在幾何梁中是如何通過(guò)幾何形狀求得的，根據(jù)這些參數(shù)，可以得到和Abaqus完全一致的剛度矩陣，從而對(duì)只有幾何梁組成的任意模型一般都能得到Abaqus完全一致的分析結(jié)果，并用一個(gè)簡(jiǎn)單的算例驗(yàn)證了該想法。

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

第十三篇：顯式和隱式的區(qū)別。介紹了顯式和隱式的特點(diǎn)，并給出一個(gè)數(shù)學(xué)算例，分別利用前向歐拉和后向歐拉求解，以求直觀表現(xiàn)顯式和隱式在求解過(guò)程中的差異，以及增量步長(zhǎng)對(duì)求解結(jié)果的影響。

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

第十四篇：殼的應(yīng)力方向。簡(jiǎn)單介紹了一下數(shù)學(xué)上張量和Abaqus中殼的應(yīng)力方向，并說(shuō)明Abaqus這么選取的意義，最后通過(guò)自編程序iSolver來(lái)驗(yàn)證殼的應(yīng)力方向的正確性。

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