有限元理論基礎及Abaqus內部實現方式研究系列40：梁單元差異（4）-形心、剪心和偏置

SnowWave02

瀏覽：5176 評論：7 收藏：15

（原創，歡迎轉載，轉載請說明出處）

1 概述

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

（1）基礎理論

（2）商軟操作

（3）自編程序

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

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

有限元理論基礎及Abaqus內部實現方式研究系列40：梁單元差異（4）-形心、剪心和偏置的圖3

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

通用結構有限元軟件iSolver介紹視頻：

有限元理論基礎及Abaqus內部實現方式研究系列40：梁單元差異（4）-形心、剪心和偏置的圖5

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

==第40篇：Abaqus、iSolver與Nastran梁單元差異（4）形心、剪心和偏置==

有限元理論基礎及Abaqus內部實現方式研究系列40：梁單元差異（4）-形心、剪心和偏置的圖7 從NASA在1964年推出Nastran算起，結構CAE軟件已經發展了將近60年，雖然在60年內商用結構軟件層出不窮，但最流行的通用結構CAE軟件現在依然只有Nastran、Ansys和Abaqus三款軟件，iSolver以前的結果精度完全對標Abaqus，線性和材料非線性與Abaqus計算結果對比精度已經沒有誤差，但在推廣過程中，發現在線性問題上，客戶只相信Nastran的結果，因為很多工程的校核規范和經驗方法都以Nastran的結果為基準，如果換成其它軟件，那么那些經驗系數和校核公式都可能會修改，而這些修改還需要做大量的工程驗證，這是客戶承受不了的代價。客戶改不了的習慣只能是自主軟件改，我們的自主軟件只是一個后來者，想要推廣只能按客戶要求修改精度，因此，iSolver開始兼容Abaqus和Nastran的精度，計算出的結果既可以對標Abaqus，修改設置后又可以對標Nastran，因此需要對Abaqus和Nastran的各方面的算法差異進行深入研究。如果是線性問題，那么Nastran和Abaqus的精度誤差主要體現在單元算法、邊界處理、MPC約束關系等，在2017年第二篇：S4殼單元質量矩陣研究文章中我們就曾經分析過Abaqus的S4殼單元和Nastran的Quad4殼單元質量矩陣的內部實現方式和差異，在這里主要研究Abaqus、iSolver與Nastran梁單元差異，由于這三款軟件的梁單元的差異較多，我們分幾篇文章來說明，本篇是Abaqus、iSolver和Nastran梁差異（4）-形心、剪心和偏置。

1.1 形心、剪心和節點坐標

一根梁雖然在有限元內部計算時采用線單元表示，但實際上只是把梁的三維形狀簡化到了參數中，如果顯示一個梁的具體形狀，譬如典型的T型梁：

將會有下面三個特殊點：

（1）形心C

就是的二維截面幾何圖形的重心。

（2）剪心S

剪心為兩個截面方向的集中力都不會導致扭轉的交點。譬如槽鋼的剪心S大體如下，顯然和形心C不是同一點：

如果在y方向加集中力，可發現只有通過剪心才能不發生扭轉，即下方的第二圖：

（3）節點位置N

節點位置就是建模時輸入的位置，譬如下面的三個節點就是0,0,0、100,0,0,和200,0,0

那么在CAE軟件中就是線顯示的位置：

Nastran中節點稱為Grid，對應bdf中的GRID關鍵字

既然梁簡化為線單元的過程中有三個特殊點：形心C、剪心S和節點N，那么問題來了：

（1）輸入的力加載在什么位置？

（2）計算采用的運動的相對位置是什么？

1.2 輸入的集中力和力矩加在哪個點

梁單元截面上的加載就是集中力、yz方向的彎矩，x軸的扭矩，我們分開討論：

1.2.1 集中力

對于集中力，實際生活告訴我們一根手指去壓迫梁截面，那么手指和截面接觸的位置顯然就是加載點。

實際力和力矩可以加載到任意點上，只要和實際的力和位移關系一致就行，但有限元讓用戶輸入集中力和力矩時沒有輸入截面坐標，只可能默認是上面三個特殊點之一。

具體加載到什么位置可以做個簡單的實驗證明：

Example1：Ex1-LoadPosition\

譬如上面的T型材，如果將節點位置N改到剪心S上方（實際不是改節點位置，而是改偏置）

將左端固定，右端在z方向加集中力載荷。

那么，可能如下：

（1）如加載到形心C，右側節點2繞x軸逆時針旋轉，得到扭轉角度UR1<0

（2）如加載到節點位置N，右側節點2繞x軸順時針旋轉，得到扭轉角度UR1>0

（3）如加載到剪心位置N，不旋轉，得到扭轉角度UR1=0

在iSolver中計算如下：

可發現UR1>0，說明是加載到節點上的。

1.2.2 力矩

力矩在有限元中和力一樣都是加載到過節點N的轉動軸上，注意力矩是輸入，與梁實際的轉動軸無關。譬如下圖節點如果在底部，那么力矩M通過力F*到節點的距離應該表示為下圖所示：

但我們找不出像力一樣可以證明力矩加載位置的例子，哪位讀者要能證明一樣也可以聯系我們。

我們查看Nastran的幫助手冊，也可發現Nastran的集中力也是加載在節點Grid上。

所以，集中力和力矩就是加載在節點Node所對應的N（x、y、z）上，只不過N（x、y、z）可能恰好是剪心，也可以是形心，當然，也可能是不在這兩個點上的任意值。

1.3 計算采用的運動的相對位置

力的加載位置和梁的真正后臺計算的相對位置不同，梁的運動分為平動、繞x軸的扭轉、繞y/z軸的彎曲轉動，平動顯然是節點本身的運行，不做討論，但轉動是繞轉動軸運動的，這個轉動軸和輸入的力矩認為轉動軸不是一回事，由上面討論輸入的力矩認為轉動軸總是過節點N的，但這邊運動的轉動軸是軟件內部計算采用的轉動軸，每個軟件還是不同的。

1.3.1 扭轉的扭矩軸

有限元中，如果受到剪心位置的扭矩，那么扭轉角沿過剪心的x方向軸計算。

即J中的距離是相對剪切中心來計算。

當剪心不在節點上時，無論材料力學理論還是Abaqus/iSolver/Nastran軟件會認為扭轉軸還是繞過剪心的軸線，而不是過節點的：

此時，顯然J不變，但需要把節點的扭矩和扭轉角變換到剪心。

當一個集中力矩移動到另一點時，力矩和扭轉角度都不會變，但繞剪心的旋轉對節點導致截面兩個方向的位移v、w的變換，且顯然與剪心和節點的距離NS有關，具體這個v、w的變換如何影響剛度也可看我們以前梁單元剛度陣的文章。

1.3.2 彎曲的轉動軸

彎矩是繞某個軸的旋轉，我們來分析究竟繞過節點還是形心旋轉？

1.3.2.1 一般的材料力學理論的推導過程

我們先從一般的彎矩公式的最原始的推導說起，在受彎矩作用下，梁做純彎曲變形，整個梁的每條平行與軸線方向的纖維都彎曲為圓形。

那么纖維SS’的相對nn’的升長量是

如果nn’是中性軸，也就是nn’的長度就是纖維SS’原來的長度，那么 ${\varepsilon }_{x}$ 就是SS’的應變。由應力應變關系：

內外力矩平衡：

其中：Iz為

得到應力：

材料力學書中一般認為y表示所在截面點到中性軸的y位置，也就是說彎曲是繞中性軸的運動。

但從上面的推導可以看出，要上面的公式成立：需滿足：

1.純彎曲情況。

2.滿足中性軸nn’是纖維SS’的沒有載荷時的原始長度。

3.彎曲后截面依然垂直中性面，否則下面等式不成立，也就是說只適合Euler梁。

1.3.2.2 有限元中的力矩彎曲軸

實際的情況很多時候都不是純彎曲，同時中性軸也會拉伸，所以到底繞哪個軸彎曲每個軟件也不同，Abaqus中取的是過節點的彎曲軸，取過節點的好處是由于上面說到的集中力和位移等都是節點的物理量，這樣所有的物理量都是節點上的，不用變換。相當于截面上力的內矩采用下圖表示：

而Nastran中彎曲軸就不同了，取的是形心的彎曲軸。類似于軟件內部要做兩步操作：

（1）先將計算轉動軸放到中性面：

（2）和上面的扭轉一樣，Nastran的梁單元在后臺需要將形心的位移轉換到節點上，譬如下方在形心處的彎曲角度θ會導致x方向的位移:

1.4 實際工程上的彎曲軸

實際的彎矩外力是加在過節點Node還是形心的軸向的呢？也就是說梁的截面是沿哪個軸彎曲的呢？個人沒做過實際的物理試驗，只是查了一下力矩的加載方法，一種加載方式是在梁的兩個邊緣加力，然后用力乘以離彎曲軸的距離當做力矩。但試驗室截面到底是觀測到過節點還是過形心個人理解很難說的清楚，如果問一個小朋友，看上面的彎曲圖，轉動角度θ是繞哪個軸旋轉，也許很多人就會說是底部，而不是中性軸，也就是說從外形上很難確定哪個軸旋轉，或者你可以認為是繞任意一個軸旋轉。

那么有限元中這兩種方式模擬結果是不是一致呢？我們只以一個簡單的例子來說明一下。

1.4.1 算例2

Example2：Ex2-TBeamNodeInShearCenter\

我們還是以T型材的一個梁單元為例。

先在iSolver中建模，為減少橫向剪切剛度的影響，我們梁長度L取為240，僅取一個單元。

材料楊氏模量=1.5e6，泊松比=0.3。

節點在剪心處，此時T型材的截面尺寸和偏置如下：

1.4.2 工況和實測結果

右端全約束，取三種工況：

工況（1）：左端自由加力，加y-方向集中力F=1：

iSolver導出到Abaqus B31單元和Nastran的CBEAM單元繞z軸轉動角度UR3結果都是5.02e6，如下：

上面說了彎曲的材料力學理論只適合Euler梁，但我們Abaqus還是采用了B31，僅僅因為B31更常用，但就算是采用B33，Abaqus內部計算依然采用過節點的轉動軸。

工況（2）：左端自由加力矩，加z+方向力矩M=1：

iSolver導出到Abaqus B31單元和Nastran的CBEAM單元UR3結果都是2.09，如下：

工況（3）：在工況（2）的基礎上左端簡支

Abaqus B31單元和Nastran的CBEAM單元UR3結果分別是5.229和4.24，如下：

1.4.3 結果分析

明顯工況（1）和（2）采用經過節點的彎曲軸的Abaqus和采用經過形心的彎曲軸的Nastran結果基本一致，但工況（3）明顯差異較大。我們沒有嚴格的證明，我們猜測無論采用哪個彎曲軸都是正確的，只是多一個平動位移的變換，都是嚴格滿足梁的外力矩=內力矩，外力=內力的關系的。所以（1）和（2）的結果一致，但如果像工況（3）把平動位移約束，彎曲軸的不同導致的差異可能就非常大了。

進一步，對于工況（3），Abaqus和Nastran那么哪個更正確呢？或者誰更可信呢？當我們劃分為100個梁單元時，發現Abaqus和Nastran的UR3結果均是4.2左右，我見過這個問題的理論值都是建立在節點在形心位置的，至于現在節點不在形心的情況，沒找到對應理論值，只能猜測4.2就是比較精確的實際工程值，那么回頭看如果僅僅是一個單元，那么Nastran更精確，當單元足夠多時，Nastran和Abaqus結果都一致。

對于Abaqus的B31一個單元的結果不如Nastran的情況我們發現過很多次了，譬如簡單的懸臂梁加彎曲力的例子，猜測Nastran的梁理論是精確解，無論單元個數多少。但Abaqus是近似解，只有單元多了才和理論解一致，這兩種情況工程上都可用。我們沒有找到對應的理論解釋，如果哪位高手知道為何這樣，歡迎聯系我們。