==概述==

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

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

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

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

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

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

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

http://www.yqgqt.org.cn/college/video/c12884 章節(jié)3

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

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

（1）單元的節(jié)點上的坐標(biāo)只要網(wǎng)格劃分完后就已經(jīng)確定了，那么單元內(nèi)的任意一點的坐標(biāo)可以表示為節(jié)點上的坐標(biāo)值。

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

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

（2）有了坐標(biāo)的表達(dá)式，由等參的定義，位移的表示和坐標(biāo)完全一致，那么就可以得到單元內(nèi)任意一點位移的表示：

（3）位移對坐標(biāo)的微分就得到應(yīng)變：

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

（4）由輸入的材料得到應(yīng)力與應(yīng)變的關(guān)系，即本構(gòu)關(guān)系

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

1.2 Abaqus的修正方式

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

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

1.2.1 沙漏產(chǎn)生的原因

單元受到彎曲作用時，正常情況下形狀應(yīng)該出現(xiàn)如下圖所示：

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

對線性減縮單元，將出現(xiàn)如下情況：

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

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

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

從剛度的數(shù)學(xué)表達(dá)式也可以解釋這個沙漏現(xiàn)象。

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

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

    

2X2的矩陣滿秩的話秩是2，但上面這樣的矩陣秩只有1，在加了載荷后是無解的，只不過在計算機中由于精度誤差，計算中的剛度矩陣將會變?yōu)轭愃葡路降木仃嚕?/p>

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

1.2.2 控制沙漏的算法

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

                                                                        [0.001, 0

                                                                        0, 0.001]

增加這個沙漏剛度單元剛度矩陣將變?yōu)橄路降木仃嚕?/p>

                                                                        [ 1.001 ，1

                                                                        2， 2.001]

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

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

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

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

對應(yīng)的沙漏位移函數(shù)為：

后面就是怎么取這個γ的算法了，與Γ和應(yīng)變應(yīng)力矩陣Bi有關(guān)，實際值為：

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

V是體積，Dhg為沙漏的本構(gòu)關(guān)系，F(xiàn)actor來控制沙漏剛度的大小。

1.2.3 Abaqus的沙漏控制方法

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

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

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

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

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

（2）              沙漏剛度會導(dǎo)致產(chǎn)生人工的沙漏應(yīng)變能，該能量在實際中并不存在，這點要求這個因子越小越好。

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

 

1.2.4 控制沙漏的實際意義

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

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

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

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

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

第二種情況：當(dāng)加入沙漏因子后。

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

1.3 算例驗證

Abaqus的修正方式僅僅是猜測，下面我們在自研求解器iSolver中對八節(jié)點六面體單元按照上面猜測的修正方法，和理論和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個不規(guī)則的斜六面體組成，內(nèi)部的四個節(jié)點任意，結(jié)果都是一樣的。在Abaqus中建模和網(wǎng)格劃分，采用C3D8R單元：

輸入材料和載荷如下：

1.3.1.2 分析結(jié)果

1.3.1.2.1 理論結(jié)果

1.3.1.2.2 iSolver結(jié)果

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

     

1.3.2 彎曲梁的算例

1.3.2.1 模型描述

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

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

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

幾何：內(nèi)徑4.12，外徑4.32，厚度h=0.1。

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

網(wǎng)格：整體劃分0.05大小的結(jié)構(gòu)化網(wǎng)格，同時單獨設(shè)置厚度方向為4個網(wǎng)格，C3D8R單元。

    

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

1.3.2.2 分析結(jié)果

1.3.2.2.1 理論結(jié)果

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

1.3.2.2.2 Abaqus結(jié)果

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

1.3.2.2.3 iSolver結(jié)果

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

1.3.2.2.4 結(jié)果分析

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

1.3.3 任意六面體組成的結(jié)構(gòu)

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

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

==總結(jié)==

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

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

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

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

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

snowwave02 From www.yqgqt.org.cn

email: snowwave02@qq.com

==以往的系列文章==

第一篇：S4殼單元剛度矩陣研究。介紹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單元如何來消除剪切自鎖以及S4R如何來抑制沙漏的。

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

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

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

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

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

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

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