積分點的案例
關于有限單元法中節（結）點與積分點的幾點釋疑
節點和積分點是有限單元法（FEM）的兩個基本概念，初涉有限元計算的同志往往在這點上產生混淆，假設導師面試的時候，問單元應力是什么，若回答不慎，將貽笑大方，得不償失。本文試圖以簡略易懂的說法來闡述節點和積分點的區別。
1.節點位移是有限元法的基本未知量。節點構筑了問題域的幾何離散化形狀，節點是形函數的零點，通常形函數是以節點為依據進行假設的。形函數決定了單元內部各點運動的位移模式（常用帕斯卡三角形來選擇單元位移模式），這樣就形成了數學上所說的插值。
有限元法的原理就是將問題域分割成N多小單元，在每個單元內采用簡單的函數來近似表達單元的真實位移，將各單元再連接起來，就可以近似描述整個問題域的運動。因此，有限元法從根本上就是精確的，而不是準確的。
2.積分點是單元進行數值積分的已知量。有限元法中一般采用高斯積分，但是積分方法不限于高斯積分，如果有人用了Irons積分或者Hammer積分，請不要驚訝。在形成單元剛度矩陣和進行節點應力磨平的時候，需要高斯積分。
以等參單元為例，其剛度矩陣
，這個就需要數值積分來快速計算，高斯點坐標及權系數如表4.2[王勖成]所示。
老師授課時一般對常應力單元進行推導，而常應力單元只有一個積分點，被積函數是常數，因此體現不出高斯積分來。很多老師對高斯積分在單元剛度矩陣的應用不予細述，導致部分同學對單元積分點認識不足。
3.單元應力指的是高斯積分點的應力，而非節點上的應力。有了位移模式，再通過虛功原理得到單元剛度矩陣，然后聚合總剛，求解平衡方程，就會把基本未知量——節點位移求出來了。通過節點位移得到單元應變結果，利用物理方程求得單元應力結果。
在等參元中，單元中n+1階(n=p-m)高斯積分點上的應變或應力近似解比其它部位具有較高的精度，因此我們稱(n+1)階高斯積分點是等參元中的最佳應力點。
展開 Abaqus中獲取積分點坐標的三種方法
經常有小伙伴問獲取積分點坐標的方法，今天給大家介紹三種獲取積分點坐標的方式，希望能給你們帶來幫助。
1 通過abaqus子程序獲取積分點坐標
Abaqus一些子程序中可以直接獲取積分點坐標，例如我們熟知的UMAT子程序中包含COORD參數，即為積分點坐標。順帶一提的是，當打開了幾何非線性時，該積分點是當前構形下的坐標，如果未打開幾何非線性則為初始坐標。
2通過history output輸出積分點坐標
Abaqus可以直接在歷程變量history output中輸出積分點坐標。直接在history output中勾選COORD選項，但是這里需要注意的是，Domain中的Set集合如果是node set，這里輸出來的是節點坐標，當這里是element set的時候，輸出來的才是積分點坐標。
3通過等參單元映射函數計算
等參元中，為了方便計算，把整體坐標映射到自然坐標，然后在自然坐標下進行高斯積分。如果知道了自然坐標下的高斯積分點，通過映射函數反算，便能得到整體坐標下的高斯積分點坐標。以四邊形等參單元為例，其以自然呢坐標為變量的插值形函數如下
坐標變換采取同樣的插值函數(叫做等參的原因)，整體坐標和自然坐標的關系式如下，如果知道自然坐標下的高斯積分點，直接通過此公式計算其在整體坐標下的坐標。
展開 有限元中單元積分點與節點應力相互轉換(CPE4為例)（ABAQUS）
（注意：變量是a,b,c,d,而不是x,y.所以方程組是線性的）
第一個積分點的應力和坐標：S1，(X1，Y1);
第二個積分點的應力和坐標：S2，(X2，Y2);
第三個積分點的應力和坐標：S3，(X3，Y3);
第四個積分點的應力和坐標：S4，(X4，Y4);
現在的問題是：應力分量S1，S2，S3，S4是已知的，我們需要知道真實的積分點的坐標信息嗎？
答案：不需要，只需要知道積分點在整個單元相對位置即可。即等參元中的坐標。（教材中有）
等參元的長和寬都為2.
而有限元中的積分是高斯積分，積分點的位置是固定的。由查表可知：
上表是一維的高斯積分點的坐標，后面的加權系數不用管（我們不求積分）。由一維可以猜出二維（兩個一維）。二維有4個積分點，所以我們對應一維選第二行的數據。
展開 有限元計算過程中積分點應力如何外插至節點處？【公式推導篇】
注：由于技術鄰排版風格有限，故部分內容顯示不全，感興趣的小伙伴可點擊原文進行閱覽：
有限元計算過程中積分點應力如何外插至節點處？【公式推導篇】
https://mp.weixin.qq.com/s/47byQ3b3e5UpbUp7Krs2mQ
本次分享的是：有限元計算過程中，單元積分點應力如何外推至節點？
有關積分點與節點的概念可點擊跳轉閱讀歷史推文：有限元基本概念-【節點和積分點】，現科普一下Q4單元、Q8單元、Q9單元的形函數和高斯積分方案。
Q4單元
Q8/9單元
應力外插
核心理念：坐標系的轉換。
假設是母單元的自然坐標系，是由高斯積分點控制的坐標系（術語可能不專業），假設高斯積分方案為。坐標系轉換關系：
單元內任一點的應力，由4個高斯積分點應力進行插值時，可表示為
其中，是基于高斯積分點的形函數，第一個積分點的坐標在母單元坐標系下為（-1，-1），根據上述的坐標系轉換的方式，在高斯積分點的坐標系下，第一個單元節點在高斯積分點坐標系下坐標為，將此坐標值代入第一個形函數，得，相同的道理，可推導至四個節點在4個形函數下的外插矩陣：
對于Q8、Q9單元，依然可采用高斯積分方案（減縮積分）。
展開 
ABAQUS輸出單元積分點坐標
方法
在ABAQUS CAE的場輸出中選擇的坐標點是節點的坐標，而節點是從積分點插值出來的，單元積分點的信息相對真實。所以最好是獲取積分點的信息，其中積分點的坐標無法在CAE中獲取，需要在關鍵字中添加。具體在每個分析步的單元輸出下面添加COORD，如果需要輸出節點的坐標也可以在節點場輸出下面添加COORD（這和CAE中場輸出選擇節點坐標的效果是一致的）。具體如下圖：
2.注意
在ODB結果中創建場輸出時會附帶著一份XYZ坐標，這個應該也可以當做單元的坐標，，但是我比較過這個附帶的坐標和單元的COORD輸出的坐標，有時候有點差別，可能是數據精度的問題。
展開 基于無網格（mesh-free）策略實現單積分點幾何必須為錯（GND）的計算
參考文獻：《Physically based crystal plasticity FEM including geometrically necessary dislocations: Numerical implementation and applications in micro-forming》
GND 演化方程依賴依賴于剪切應變率的梯度或者塑性變形梯度的旋度，而標準FEM/VUMAT 只告訴你每個積分點本身的 γ˙a、Fp，不會直接給梯度。以往廣泛應用的數值方案通常是：先把 積分點的數據外推到節點，再用線性形函數求梯度，然而這類方案只能用特定單元（如 C3D8），對自適應網格、復雜接觸不友好。
該文章提出的一個mesh-free的方案，該方案的主要優勢是不改單元、不加 DOF，只在材料子程序內部，用鄰近積分點的數據做一次局部重構，就算出梯度，該策略對某個積分點 x，附近有一團“鄰居積分點” xI，作者把它們當成 mesh-free 的“節點”，對每個場變量 u(x)（可以是 γ˙a，Fp 的分量）做 MLS 擬合，如下圖所示：
權函數使用立方樣條，有緊支撐，距離越近權越大：
在實現上作者提到，立方支撐三個方向尺寸約為5個單元尺寸，最多取最鄰近60個（3D）或者30個（2D）積分點，作者指出：當鄰域尺寸比網格尺寸還小的時候，這個非局域模型就自然退化為傳統的局域模型。也就是說，鄰域尺寸本身就扮演了“材料內在長度標度”的角色。
為了提高計算效率秒作者使用了一個“時間滯后 + 公共塊”的策略對GND進行更新。
展開 單元積分點應力如何外插至節點上 | 數值實現篇
繼上次的推文：有限元計算過程中積分點應力如何外插至節點處？【公式推導篇】，本次分享單元積分點應力外插至節點處的數值實現過程。
數值實現
借助以上理論，我們可以基于matlab平臺編制以下代碼段：
% 將積分點應力外插至單元節點上，這里只列舉了Q4的情況
for i = 1:3
StressElem(e,:,i) = [1+0.5*sqrt(3) -0.5 1-0.5*sqrt(3) -0.5;
-0.5 1+0.5*sqrt(3) -0.5 1-0.5*sqrt(3);
1-0.5*sqrt(3) -0.5 1+0.5*sqrt(3) -0.5;
-0.5 1-0.5*sqrt(3) -0.5 1+0.5*sqrt(3)]*...
[stress(e,1,i);stress(e,2,i);stress(e,3,i);stress(e,4,i)];
end
對標Abaqus
模型材料參數為普通的線彈性材料，單元類型選擇CPS4，網格劃分及邊界條件設置如下：
在結果對標過程中，可以先對比自研程序與Abaqus的節點位移場：
Abaqus位移場結果
自研程序位移場結果
在位移場一致的前提下，我們再來對標應力結果。以常見的mises應力為例：
Abaqus位移應力場結果
自研程序應力場結果
結果是一致的，說明了程序的正確性。
展開 abaqus C3D8 單元 計算中采用了多少個積分點？
按照正常的理解，毫無.疑問，abaqus 全積分一定是采用了2x2x2=8個積分點。
從后處理結果來看，似乎也是如此，每個單元存在8個積分點。
然而，如果自己動手跑一遍程序，就會發現事實遠非如此，采用全積分計算得到的結果與abaqus 存在差異，原因何在？
事實賞，abaqus C3D8 采用的選擇積分方式（selective intergation schema），即對于偏應變，采用8個積分，對于球應變，采用中心點積分。這樣計算得到的結果才能與abaqus 完全對標，亦可從abaqus 幫助文檔得到答案。
展開 ABAQUS 輸出節點坐標和積分點坐標
總結inp中添加關鍵字
輸出單元的積分點坐標：*EL FILE
COORD
輸出節點坐標：*NODE FILE
COORD
原貼出處：https://www.researchgate.net/post/How-to-find-integration-point-coordinates-in-Abaqus-CAE
這是帖子討論的，但是我的嘗試是兩個COORD生成的結果文件是一樣的，都是節點坐標
Abaqus中平面應力單元高斯積分點的順序
可以輸出umat接口中的變量coords進行查看
write(*,"(A,I4)") "npt = ", npt
write(*,"(A,3ES16.8)") "coords = ", coords
結果為：
npt = 1
coords = -5.77350269E-01 -5.77350269E-01 1.00000000E-02
npt = 2
coords = 5.77350269E-01 -5.77350269E-01 1.00000000E-02
npt = 3
coords = -5.77350269E-01 5.77350269E-01 1.00000000E-02
npt = 4
coords = 5.77350269E-01 5.77350269E-01 1.00000000E-02
因此Abaqus中平面應力單元高斯積分點的順序為：
展開 07-通過多重積分法求解點接觸彈性變形Fortran和MATLAB程序 ￥19.89
07-通過多重積分法求解點接觸彈性變形Fortran和MATLAB程序，程序請見下文附件及百度網盤鏈接
Abaqus與Dyna電池包沖擊分析結果對比 ￥20
模型設置對比
電池包模型
計算結果對比
對比一：Dyna中不同積分方式及積分點數量對比
Dyna中殼單元積分方式有2種，一種為默認的高斯積分方式，代數精度為2n-1；另一種改進的高斯積分方式Gauss-Lobatto，代數精度為2n-3。所以對于單元形函數階次在3以下的單元，高斯積分在殼單元厚度方向上2個積分點已經能滿足精度要求，而Gauss-Lobatto則需3個以上積分點才能滿足精度要求。需要注意的是，這里的精度僅是積分點位置上的數值精度，而對于殼單元，我們更關注的是殼單元上下表面的數值精度。
不同積分方式積分點位置
對于彈性問題，采用默認積分點數基本能滿足精度要求。
對于彈塑性問題，由于殼單元厚度方向上不再是線性變化，需要更多積分點數來描述殼單元厚度方向上的彈塑性變化狀態，一般需至少5個積分點才能保證計算精度。
在同樣采用2號殼單元類型的情況下，分別采用殼單元厚度方向默認的積分點2+高斯積分，5積分點+高斯積分，5積分點+改進的高斯積分Gauss-Lobatto進行對比。
計算結果
2積分點與5積分點對比（虛線為2積分點）
Gauss 與Gauss-Lobatto積分點對比
由計算結果可知：
對于本例而言，采用默認的2積分點最大值雖然與5積分點差別不大，但是局部應力曲線差別很大。對于5積分點而言，Gauss積分與Gauss-Lobatto積分差別不大，但由于Gauss-Lobatto積分點最外層在殼表面，所以結果更精確。
展開 梁單元的彈塑性-彈塑性梁單元在長度上任意位置都會考慮塑性嗎
在大部分的有限元軟件中，在梁長度方向上會設置若干個積分點，計算時僅僅會捕捉積分點的應力判斷是否進入塑性。
例如，對于abaqus的B33單元，在長度方向上有3個高斯積分點。其具體位置為：(0.1127016L，0.5L，0.887298L)；對于B31，在長度方向上僅一個高斯積分點，位置為中點處。
以下圖的B33為例：
長度為1m，截面為0.1m*0.1m的梁采用1個B33單元，左端約束，右端施加豎向荷載Fz=1N.
計算完成后查詢積分點的S11應力值：
按照前文提到的長度方向積分點的位置為：(0.1127016L，0.5L，0.887298L)，則三個積分點處的應力（截面頂或者底）計算為:
同理可計算M2和M3，結果均與abaqus查詢的結果一致。
在彈塑性材料本構計算中，軟件就是通過捕捉積分點的應力，從而依據積分點的應力值判定是否進入塑性而進行應力更新的，而不是“在整個長度上各個位置都考慮了塑性”。當然，在長度方向上的每個積分點的位置，軟件不僅僅是計算橫截面的頂或者底的應力，實際上橫截面上又會采用多個積分點進行辛普森積分，在abaqus中，這種橫截面的積分點叫section points。一個矩形截面共有25個section points。軟件對長度方向上每個積分點的25個section points均計算應力，并進行彈塑性判斷與應力更新。因此，對于一個B33單元，積分點總數實際上是25*3=75個。
以上，即是本文的全部內容，感謝您的閱讀。歡迎關注公眾號 有限元術
展開 ANSYS中的節點解與單元解是怎么回事？下次別說你還不懂
無疑，這樣沒法在軟件中顯示結果，因此單元解需要確定一些積分點（高斯點），通過積分得到這些積分點的解，這些積分點的解代表單元解。
積分點通常和單元的節點位置不重合，因此想要得到單元節點的解，需要將積分點的解根據某種規則外推，以一種近似的方法得到單元節點的解。由于每個單元外推得到的單元節點解并不完全一致，因此，最初外推得到的單元的節點解不連續，為了讓其連續，將不同單元之間的節點外推得到的節點解進行算術平均，這樣在連續節點處的節點解僅有一個數值，這樣便得到實際在軟件中顯示的節點解。
簡短一點來說：單元解是積分點的解，節點解是外推后平均的解。很明顯，從數值精度上來講，單元解是高于節點解的。
采用ANSYS計算了一個簡單的模型，分別采用solid185單元和solid186單元，185單元是8節點單元，186單元是20節點單元，分別計算后查詢；
最終，單元總數185為256個，186為256個，單元劃分一樣，但是節點數不一樣，185單元劃分的模型節點數為459個，186單元劃分的為1605個。
查看ANSYS計算輸出的單元解，當單元為185時查詢兩個挨著的單元應力解如圖1所示：
圖1
當單元為186時查詢兩個挨著的單元應力解如圖2所示：
圖2
經過以上計算可以看出：
（1）無論是185單元還是186單元，計算后的單元解只輸出8個節點的值，這個非常奇怪，因為185單元和186單元的積分點數目不一樣，185為8個積分點，186為27個積分點；
（2）相鄰單元的共同節點的應力值不一樣，這個是合理的，因為每一個單元的節點解是根據各自的形函數計算并且外推的，有差別。
這里就留下一個問題，為什么186單元也只輸出八個節點的值？
展開 ANSYS中的節點解與單元解是怎么回事？附solid186與solid185單元結果對比文檔下載
無疑，這樣沒法在軟件中顯示結果，因此單元解需要確定一些積分點（高斯點），通過積分得到這些積分點的解，這些積分點的解代表單元解。
積分點通常和單元的節點位置不重合，因此想要得到單元節點的解，需要將積分點的解根據某種規則外推，以一種近似的方法得到單元節點的解。由于每個單元外推得到的單元節點解并不完全一致，因此，最初外推得到的單元的節點解不連續，為了讓其連續，將不同單元之間的節點外推得到的節點解進行算術平均，這樣在連續節點處的節點解僅有一個數值，這樣便得到實際在軟件中顯示的節點解。
簡短一點來說：單元解是積分點的解，節點解是外推后平均的解。很明顯，從數值精度上來講，單元解是高于節點解的。
采用ANSYS計算了一個簡單的模型，分別采用solid185單元和solid186單元，185單元是8節點單元，186單元是20節點單元，分別計算后查詢；
最終，單元總數185為256個，186為256個，單元劃分一樣，但是節點數不一樣，185單元劃分的模型節點數為459個，186單元劃分的為1605個。
查看ANSYS計算輸出的單元解，當單元為185時查詢兩個挨著的單元應力解如圖1所示：
圖1
當單元為186時查詢兩個挨著的單元應力解如圖2所示：
圖2
經過以上計算可以看出：
（1）無論是185單元還是186單元，計算后的單元解只輸出8個節點的值，這個非常奇怪，因為185單元和186單元的積分點數目不一樣，185為8個積分點，186為27個積分點；
（2）相鄰單元的共同節點的應力值不一樣，這個是合理的，因為每一個單元的節點解是根據各自的形函數計算并且外推的，有差別。
這里就留下一個問題，為什么186單元也只輸出八個節點的值？
展開 