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

==概述==

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

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

iSolver介紹視頻：

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

==第14篇：殼的應力方向 ==

有限元中，物理量用的最多的是標量、矢量和二階張量。其中位移、坐標等都是矢量，而應變、應力等都是二階張量。矢量很容易理解，體的應力等二階張量直接就采用了全局坐標系的也不會有方向理解問題，但梁殼的應力結果很容易搞錯，后處理結果中的S11、S12等的方向有時會覺得和預想的不一致但又不明所以。同時，這個方向也是單元材料的方向，在自編程序時，如果一開始坐標系的定義就弄錯了，那么將直接導致和材料相關的剛度矩陣的錯誤，所以弄清應力的方向定義對自編程序和理解有限元結果都相當重要。

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

具體的驗證詳見下方視頻（帶配音）：

http://www.yqgqt.org.cn/college/video/c12884 20.14 理論系列文章14：殼的應力方向

1.1 數學上的張量方向

矢量的方向是一定的，但它的分量都是基于某個坐標系定義的，坐標系不同，那么分量結果也會不同。

矢量可以表示為：                     

顯然，分量和坐標系的選取有關。譬如我們一般的直角全局坐標系如下，那么分量就是普通的x、y、z三個分量值。

和矢量類似，二階張量可以表示如下，當然也可以用一個更簡單的3X3的矩陣表示，顯然，二階張量的分量等也與坐標系的取值有關。

1.2 Abaqus殼的應力方向

Abaqus后處理中殼的應力會輸出S11，S22，S12等分量，分別對應上面二階張量a的a11、a22、a12等分量，其它分量不輸出，這三個量與殼的坐標系的選取密切相關。Abaqus中，在《anlysis user manual 12.1》的29.6.7 Three-dimensional conventional shell element library的Element output說明了殼采用的是局部坐標系，具體定義如下：

S11：Local 11 direct stress.

S22：Local 22 direct stress.

S12：Local 12 shear stress.

（1）   local 1（以下稱為T1方向）默認情況下為x方向在表面上的投影。

如果x方向正好垂直于表面(Abaqus取cos夾角<0.1度）,那么取z方向投影，所以下方的圓柱面上當x軸和表面垂直時，S11的方向明顯和上下不一致了。

（2）n方向（即T3方向）垂直于殼的表面，至于是往上還是往下，由節點順序決定。

（3）第三個方向local 2（以下稱為T2方向）和T1，T3滿足右手定則。顯然在殼的表面繞T3轉動90度。

由上面的定義，可以看出應力的方向的S11，S22都是沿著表面的，而不是全局坐標系下的。

1.3 應力方向選取原因

殼的應力方向為何要取成上面的T1，T2和T3？有兩個原因：

1.3.1 插值函數

T1，T2在面內才能使得坐標和位移可以使用二維平面內的插值公式。

位移的插值是有限元的基礎，插值方法如下：

  殼的插值方式有兩種，一種是平面殼的理論，先按平面殼來計算K，然后再將K做坐標變換，此時單元內部任意點的坐標值只有x，y，沒有z，插值為四個節點的坐標。這種方法就要求在殼的坐標系下，四個節點的z方向的坐標zi=0，只剩xi和yi，這種方法計算精度不是很準，因為一個四邊形的四個點不能總是保證在一個平面上的，只要是曲面劃分成四邊形很有肯定就會是這種情況。另一種更精確的算法是當做曲面殼，abaqus就是這么做的。此時坐標的插值函數不變，但換為了三個坐標x、y、z的分別插值，四個節點的z方向的坐標zi不需要強制要求是0了。

不管是哪種方法，插值函數都是二維的，如下：

只有T1、T2在殼平面內實際坐標才能映射到等參上。

1.3.2 本構關系

殼的應變方向為何要取成T3和面垂直的另一個原因是因為只有這樣用殼來簡化分析對象時材料的本構關系才是最簡單的。因為殼是對體的簡化，當體的厚度遠小于面內尺寸時，那么可以用殼的理論來近似，此時可以把殼的剛度分為薄膜效應剛度、面外彎曲剛度、面外橫向剪切剛度等部分，每一部分本構關系由相應簡單的材料相關的D矩陣決定，譬如薄膜效應和面外彎曲的D矩陣都是下方矩陣：

上面的本構關系是各向同性材料的，面內的T1和T2不影響D矩陣，但各向異性就不一樣，此時D矩陣將與面內的T1，T2相關，這也是應力的坐標系也叫做材料坐標系的原因。

1.4 程序內部實現方法

其實，原理非常簡單，但程序實現并不像看起來的那么簡單。無論是Abaqus還是自編程序內部想要實現殼的應力的方向，和前面討論的應力方向選取原因是一致的，主要是兩點：

（1）坐標和位移用局部坐標系表示，也就是下面等式：

S1，S2，S3分別對應T1，T2，T3三個方向的坐標分量。

（2）本構關系用局部坐標系做旋轉變化。

1.5 iSolver驗證

為了證明殼的S11、S22是沿表面方向而不是全局坐標系方向，我們用一個球殼加均勻壓力，可發現S11和S22在abaqus中是球對稱的，只有沿著表面才會是應力相等，如果沿著全局的xyz，那么顯然不應該球對稱。

為了進一步驗證，我們在iSolver中按照Abaqus的T1、T2、T3定義殼單元方向，并通過一個簡單的例子來驗證Abaqus殼的應力方向。

1.5.1 算例

只取一個簡單的長方形單元。

參數如下：

尺寸：5X1，厚度0.1。

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

右側兩個節點固支。

左側兩個節點每個加集中力1e5，在殼的平面內往外拉伸。

一開始在XY平面內，長度5方向在X軸，寬度方向在Y軸上。此時將得到一個應力的二階張量S。

1.5.2 Y軸轉動45度

沿Y軸轉動45度，按照Abaqus的定義方式，顯然，二階張量S完全不變。

1.5.3 Z軸轉動45度

沿Z軸轉動45時，顯然，殼的應力將不同，需要坐標變換。

設R為從原始模型的到新模型的變換矩陣，那么從數學的角度如果是一個列矢量a，那么變換后的矢量為：

A=R*a

而如果是二階張量S，在原始坐標系下二階張量可以表示為

S=a*b’

b’為列矢量b的轉置。

顯然變換后的張量為：

SNew=(R*a)*(R*b)’=R*a*b’*R’=R*S*R’

其中R’表示轉置。

1.5.4 iSolver驗證演示

iSolver中詳細的驗證證明請參考演示視頻（帶配音）：

http://www.yqgqt.org.cn/college/video/c12884 20理論系列文章14：殼的應力方向

==總結==

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

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

snowwave02 From www.yqgqt.org.cn

email: snowwave02@qq.com

以往的系列文章：

第一篇：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

第十一篇：自主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