==概述==

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

一方面我們查閱Abaqus軟件手冊得到修正方法的說明，另一方面我們自己編程實現簡單的結構有限元求解器，通過自研求解器和Abaqus的結果比較結合理論手冊如同管中窺豹一般來研究Abaqus的修正方法，從而猜測商用有限元軟件的內部計算方法。在研究的同時，準備將自己的研究成果記錄下來寫成一個系列文章，希望對那些不僅僅滿足使用軟件，而想了解軟件內部實現方法甚至是做自己的軟件的朋友有些幫助。由于水平有限，里面可能有許多錯誤，歡迎交流討論。

==第七篇：C3D8R六面體單元的剛度矩陣==

前面幾篇文章都是介紹了梁和殼的結構，其實梁和殼只是體的兩種特殊形式，一般情況下體才是一個結構正常的形狀，而且當計算機足夠強大時，僅僅只要用體單元就能求解結構問題，這個時候甚至可以不用做幾何簡化和中面提取等輔助操作了，大大減少了前處理的時間。

前面文章提到梁相對殼來說，商業軟件的修正方式相對較少，如果自己編程序，采用這些修正方式可以得到和商業軟件完全一致的梁單元剛度矩陣，如果剛度矩陣完全一致，那么對任何的梁的算例都可以得到和商業軟件完全一致的結果了。進一步，對體來說，商用軟件對體的修正方式相對梁更少，因為梁對軸向拉伸、軸向扭轉、橫向彎曲、橫向剪切等的剛度都需要不同的修正，但體只需要一個方向的修正就行，只要掌握了體的修正方式，那么自編程序的體單元組成的任意結構分析都可以得到商業軟件完全一致沒有誤差的結果。

在Abaqus的Standard和Explicit求解中默認的體單元都為C3D8R，C表示體單元，3D表示三維單元，8表示八節點，R表示Reduced減縮積分，也就是八節點一次線性六面體單元，該單元在Abaqus中是計算耗時最少的六面體單元，而六面體的計算精度本身就比四面體和鍥形單元高，所以C3D8R在體單元實際應用中最廣泛，我們本章以此為研究對象。

本文先討論一般的六面體單元的基本理論和C3D8R的修正方式，然后在自編有限元程序iSolver實現同樣的修正方式，最后驗證iSolver的結果和Abaqus完全一致，從而證明Abaqus的內部修正和我們設想的一致。具體驗證過程也可以參考我們的演示錄像：

http://www.yqgqt.org.cn/college/video/c12884 章節3

1.1 六面體單元剛度矩陣的基本理論

體單元剛度矩陣的理論相對梁和殼來說要簡單很多。按照有限元的一般的等參理論：

（1）單元的節點上的坐標只要網格劃分完后就已經確定了，那么單元內的任意一點的坐標可以表示為節點上的坐標值。

i=1到節點N，在六面體單元中一次和二次單元節點個數一般為8、20，四面體分別為4、10。20節點六面體的另外12個節點正好位于12條邊的中點，10節點四面體類似。

Ni表示形函數，類似于一個比重，所求的x點約靠近xi點必然這個比重越大，最大顯然應該是x點位于xi點上，此時比重就是1，其它的形函數為0。

（2）有了坐標的表達式，由等參的定義，位移的表示和坐標完全一致，那么就可以得到單元內任意一點位移的表示：

（3）位移對坐標的微分就得到應變：

• 注意Abaqus的剪切應變順序為12，13，23，和這邊的B矩陣順序不同。

（4）由輸入的材料得到應力與應變的關系，即本構關系

（5）體的剛度矩陣就是：

1.2 Abaqus的修正方式

一個[-1,1]內連續的積分此時就轉換為了離散的積分，積分點的數目由gi的個數決定。最簡單的情況就是取一個積分點，就是C3D8R單元，表示減縮積分。Abaqus的C3D8R除了減縮積分外，還做了下面的修正。

只要是減縮積分都會存在沙漏，C3D8R也不例外。

1.2.1 沙漏產生的原因

單元受到彎曲作用時，正常情況下形狀應該出現如下圖所示：

上方材料纖維受拉，下方纖維受壓，導致上方是拉力，下方是壓縮力，這樣使得在左右邊界力矩保持平衡，如下圖所示。

對線性減縮單元，將出現如下情況：

由于單元拉伸應力只計算中心點P的值，而中心點P在X方向的纖維L1的長度沒有變化，使得單元拉伸應力為0，也就是和外力矩M無法平衡，導致在M作用下單元扭曲的非常厲害。

一個單元扭曲的形狀如上圖，當存在多個這樣的單元時，會形成如下的形狀。上下兩個單元就類似一個沙漏。

顯然，沙漏只存在于彎曲現象中，對拉伸沒有此現象。

從剛度的數學表達式也可以解釋這個沙漏現象。

在減縮積分時，只有一個積分點，那么上述值就是

上述應變與位移的關系B0矩陣為中心點P的值，存在多個項是相同的，使得剛度陣K也存在多個相同值的元素，導致剛度矩陣秩的缺失，舉個簡單的例子，就容易導致類似下方的剛度矩陣。

    

2X2的矩陣滿秩的話秩是2，但上面這樣的矩陣秩只有1，在加了載荷后是無解的，只不過在計算機中由于精度誤差，計算中的剛度矩陣將會變為類似下方的矩陣：

這時只要加一個非常小的載荷將導致解非常大。這也是當存在沙漏時，單元變形明顯變大的原因。

1.2.2 控制沙漏的算法

沙漏控制現在最常用的方法是增加一個人工的沙漏剛度，譬如：

                                                                        [0.001, 0

                                                                        0, 0.001]

增加這個沙漏剛度單元剛度矩陣將變為下方的矩陣：

                                                                        [ 1.001 ，1

                                                                        2， 2.001]

這個時候解依然存在，不一定是正確解，但肯定比由于計算機精度造成的誤差小很多了。

那么具體如何增加沙漏剛度和大小應該取多少合適呢？

這個沙漏控制應該只對沙漏模式（對應的位移函數設為Γ）起作用，但對其它模式（對應的位移函數為O）不起作用，也就是在沒有沙漏現象的正常情況這個剛度應該不起作用，這點通過增加一個歸一化的γ矩陣，使得γ*Γ=1，同時對其它模式γ*O=0來實現。

在一次減縮積分的六面體單元中，以x方向為例，存在四個沙漏模式

對應的沙漏位移函數為：

后面就是怎么取這個γ的算法了，與Γ和應變應力矩陣Bi有關，實際值為：

和普通的剛度矩陣一樣，沙漏剛度的值為

V是體積，Dhg為沙漏的本構關系，Factor來控制沙漏剛度的大小。

1.2.3 Abaqus的沙漏控制方法

Abaqus也是通過增加一個沙漏剛度來控制C3D8R的沙漏：

即在Element Type中可以設置Hourglass Control的方法，Standard中default就是stiffness的控制方法，只有Explicit其它幾項才有作用，同時設置了一個縮放因子Scaling factors來控制沙漏剛度的幅度。

這個沙漏剛度在物理上沒有任何的物理意義，是人為加進去的，所以由它計算出來的任何量是沒有意義的，這也是abaqus后處理中ALLAE這個變量中的AE表示Artificial Strain Energy的含義。因為沒有意義，想要結果正確，那么這個人工的應變能相對總的內能ALLIE就需要越小越好。一般認為ALLAE/ALLIE的比值必須控制在5%-10%以下，才認為計算結果可信。

沙漏的因子Factor很多時候需要修改，它的取值受下面兩者影響：

（1）              有效抑制沙漏模式，這點要求這個因子越大越好。

（2）              沙漏剛度會導致產生人工的沙漏應變能，該能量在實際中并不存在，這點要求這個因子越小越好。

上面兩點本身就是矛盾的兩點，所以很多時候這個因子根本不存在，當出現沙漏的時候，修改沙漏因子并不能真的解決沙漏問題。

 

1.2.4 控制沙漏的實際意義

在前面章節S4R的沙漏控制也提到了，通過算法來控制減縮單元譬如C3D8R/S4R的沙漏，這個沙漏因子很多情況下沒法同時滿足不增加太多的人工剛度又控制了沙漏，但Abaqus還是默認在C3D8R和S4R中增加了一個小的沙漏因子，既然這個小的沙漏因子對存在的沙漏無法得到有效控制，那么意義在哪里？

可能是筆者有誤，對Abaqus了解的還不夠深入，也希望跟大家一起探討，只能猜測一種可能是這種沙漏對Explicit分析及其重要，畢竟Abaqus的Explicit對沙漏的控制要比Standard有更多的選擇，另一種可能情況是不在于這個因子的大小，而只在于是否有這個因子，如果沒有沙漏，那么只要涉及到剛度矩陣求逆的操作都直接失敗，而存在這個增加的人工沙漏后可以讓剛度矩陣的秩增加，從而所有矩陣求逆的操作都能順利進行。以殼單元S4R為例，譬如下面的這個例子：

計算一個四周簡支板的模態分析，采用減縮殼單元S4R。

第一種情況：如果沒有沙漏因子，Abaqus就會計算失敗。

采用iSolver自研程序計算，可得到下面的前10階模態，和理論值相比，可以發現多了前三項不存在的偽模態。猜測這可能也是Abaqus計算失敗提示的內部意義。

第二種情況：當加入沙漏因子后。

Abaqus和iSolver都可以計算出前10階模態精確值，說明沙漏的存在使得模態分析求解順利進行。

1.3 算例驗證

Abaqus的修正方式僅僅是猜測，下面我們在自研求解器iSolver中對八節點六面體單元按照上面猜測的修正方法，和理論和Abaqus比較從而驗證我們的猜測。

1.3.1 分片試驗

1.3.1.1 模型描述

模型：該算例取之Abaqus Verification Manual 6.12-1的1.5.2 Patch test for three-dimensional solid elements，是Abaqus用來證明體單元的分片試驗。inp文件為：Job-PatchTest-Solid-7Hexa.inp，我們只取第一個step。

模型由7個不規則的斜六面體組成，內部的四個節點任意，結果都是一樣的。在Abaqus中建模和網格劃分，采用C3D8R單元：

輸入材料和載荷如下：

1.3.1.2 分析結果

1.3.1.2.1 理論結果

1.3.1.2.2 iSolver結果

S11和S12的結果和理論完全一致，證明iSolver的六面體單元通過基本的分片試驗。

     

1.3.2 彎曲梁的算例

1.3.2.1 模型描述

模型：該算例取之MacNeal論文[2]，為當初驗證Nastran體單元的一個例子。

inp文件為Job-CurvedBeam-Z.inp。

模型為一個1/4的半圓形的梁，一端固支，另一端受面外拉力1N，分析受力點在面外的位移。

幾何：內徑4.12，外徑4.32，厚度h=0.1。

材料：楊氏模量1e7，泊松比0.25。

網格：整體劃分0.05大小的結構化網格，同時單獨設置厚度方向為4個網格，C3D8R單元。

    

邊界和載荷：一端固支，另一端受面外拉力1N。

1.3.2.2 分析結果

1.3.2.2.1 理論結果

由論文上可得理論上U3最大位移為0.5022。

1.3.2.2.2 Abaqus結果

Abaqus分析后得到U3最大值=0.5332，相比理論值為1.06：

1.3.2.2.3 iSolver結果

Abaqus分析后得到U3最大值=0.5332，相比理論值為1.06，和Abaqus完全一致：

1.3.2.2.4 結果分析

Abaqus和iSolver的C3D8R的位移結果都比理論值偏大一些，是因為彎曲載荷下存在沙漏作用，沙漏控制只是減弱了位移的幅度，但剛度還是比實際的要偏軟。如果采用C3D8完全積分，可得到下面的位移，明顯比理論值偏小，此時沒有沙漏現象，但存在剪切自鎖使得單元偏剛。

1.3.3 任意六面體組成的結構

對任意六面體單元組成的結構，iSolver的結果都和Abaqus完全一致，具體可看演示錄像：

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

==總結==

本文首先簡單介紹了一下8節點六面體單元的基本理論，分析了剛度矩陣的來源，并研究了Abaqus中應用最廣泛的C3D8R單元的剛度矩陣的修正方式，自編程序如果采用這些修正方式可以得到和Abaqus完全一致的分析結果，并用多個算例驗證了該想法。

Abaqus的C3D8R六面體單元的剛度矩陣在八節點線性單元理論基礎上的修正如下表：

項次	問題	修正情況	說明
修正	不修正
1	積分	√		采用減縮積分
2	沙漏現象	√		增加人工的沙漏剛度對沙漏進行控制

有興趣的可以自行下載iSolver進行驗證，因為看不到Abaqus的源代碼，上述C3D8R的修正方式也僅是猜測，如果你在使用iSolver測試其它的由C3D8R組成的模型結果時發現和Abaqus結果不一致，歡迎聯系我們。

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

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