高斯積分的案例
有限元中的高斯-勒讓德積分
當函數表達式比較復雜時，f(x)的原函數可能難以求出，而采用高斯積分，其省去了求f(x)原函數，只需要將數值代入f(x)的表達式即可求解。</p><p>到目前為止，高斯積分的公式已經介紹完成，那么有兩個最直接最現實的問題出現了：（1）f(x)的表達式是什么形式時適合采用高斯積分，精度怎么樣；（2）xi和wi的取值是多少。</p><p>關于（1），實踐表明，當f（x）的表達式為多項式時，高斯積分是合適的，并且，n點高斯積分可以準確積分2n-1次多項式。</p><p>關于（2），xi和wi的取值一般較多的有限元教科書中會給出數值，如果沒有給出數值，也可以用多項式手動算出具體值，另外，scipy庫，PETSc庫也直接給出了高斯積分的值和權重。</p><p><img src="https://img.jishulink.com/msimage/202308/f410d6e3bb89c8459660277304592181.png"></p><p>以下是高斯積分點積分多項式的一個例子：</p><div contenteditable="false" width="100%">
<img src="https://img.jishulink.com/upload/202308/865d98a129374c668e89080010b652c9.jpg" title="圖片4.jpg" alt="圖片4.jpg" style="max-width:760px;" data-mobile-src="https://img.jishulink.com/upload/202308/865d98a129374c668e89080010b652c9.jpg?
展開 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中平面應力單元高斯積分點的順序為：
展開 有限元中的高斯點與積分詳解(上)_《數值計算與程序設計》系列課程之五 ￥599
本課從實際問題出發，帶著問題去講解有限元中的高斯點與數值積分。一開始拋出了以下3個關鍵問題：
1.對于一個任意函數怎么去得到它的積分？
2.數值積分的本質是什么？為什么簡單地取幾個點就可得到積分值？此種方法的立足點在哪？
3.很多資料上都說“有限元求解精度嚴重依賴于網格質量，過度扭曲的單元會導致結果不收斂或者精度極度惡化”，這只是為什么呢？扭曲單元到底影響的是有限元方法中的哪一步？
圍繞這3個問題，本課分別講了一下三個內容：
1. 數值積分基本方法。
2. 有限元單元積分。
3. 誤差分析。
本次課程分為上下兩課，第一課講了第一和第二個內容。關鍵詞是：數值積分的本質，有限元高斯積分（附件中包含1個小時的詳細課程視頻以及PPT）。
在第二課中，再繼續展開第三部分內容，誤差分析，解決問題“扭曲單元到底影響的是有限元方法中的哪一步”。
希望有興趣的同學多多支持下，你們的支持是我更新的動力
展開 有限元計算過程中積分點應力如何外插至節點處？【公式推導篇】
注：由于技術鄰排版風格有限，故部分內容顯示不全，感興趣的小伙伴可點擊原文進行閱覽：
有限元計算過程中積分點應力如何外插至節點處？【公式推導篇】
https://mp.weixin.qq.com/s/47byQ3b3e5UpbUp7Krs2mQ
本次分享的是：有限元計算過程中，單元積分點應力如何外推至節點？
有關積分點與節點的概念可點擊跳轉閱讀歷史推文：有限元基本概念-【節點和積分點】，現科普一下Q4單元、Q8單元、Q9單元的形函數和高斯積分方案。
Q4單元
Q8/9單元
應力外插
核心理念：坐標系的轉換。
假設是母單元的自然坐標系，是由高斯積分點控制的坐標系（術語可能不專業），假設高斯積分方案為。坐標系轉換關系：
單元內任一點的應力，由4個高斯積分點應力進行插值時，可表示為
其中，是基于高斯積分點的形函數，第一個積分點的坐標在母單元坐標系下為（-1，-1），根據上述的坐標系轉換的方式，在高斯積分點的坐標系下，第一個單元節點在高斯積分點坐標系下坐標為，將此坐標值代入第一個形函數，得，相同的道理，可推導至四個節點在4個形函數下的外插矩陣：
對于Q8、Q9單元，依然可采用高斯積分方案（減縮積分）。
展開 
你不知道的CAE小常識（二十一）（二十二）
在等參元中，單元中n+1階(n=p-m)高斯積分點上的應變或應力近似解比其它部位具有較高的精度，因此我們稱(n+1)階高斯積分點是等參元中的最佳應力點。因此有限元法中是以通過節點位移根據位移模式推出高斯點位移，然后求得應變和應力結果（苦逼的常應力單元，單元內各點應變、應力是一樣的，so，你懂的）。
4.高斯點是等參單元內最佳應力點。高斯積分點雖然是單元內最佳應力點，但在節點上并非最佳，節點應力有一種Newton-cotes積分方法。因此，有的FEM程序會采用其他積分方法以獲得單元邊界或者節點上的精確應力。
5.節點力是高斯點應力在單元內積分，再累加其他單元的相同節點結果得到節點力;而節點應力是通過高斯點應力外推通過平均各單元外推應力而得（應力磨平）。
6.常見單元的高斯積分點位置示意。
歡迎關注微信公眾號：DR有限元仿真
展開 關于有限單元法中節（結）點與積分點的幾點釋疑
節點和積分點是有限單元法（FEM）的兩個基本概念，初涉有限元計算的同志往往在這點上產生混淆，假設導師面試的時候，問單元應力是什么，若回答不慎，將貽笑大方，得不償失。本文試圖以簡略易懂的說法來闡述節點和積分點的區別。
1.節點位移是有限元法的基本未知量。節點構筑了問題域的幾何離散化形狀，節點是形函數的零點，通常形函數是以節點為依據進行假設的。形函數決定了單元內部各點運動的位移模式（常用帕斯卡三角形來選擇單元位移模式），這樣就形成了數學上所說的插值。
有限元法的原理就是將問題域分割成N多小單元，在每個單元內采用簡單的函數來近似表達單元的真實位移，將各單元再連接起來，就可以近似描述整個問題域的運動。因此，有限元法從根本上就是精確的，而不是準確的。
2.積分點是單元進行數值積分的已知量。有限元法中一般采用高斯積分，但是積分方法不限于高斯積分，如果有人用了Irons積分或者Hammer積分，請不要驚訝。在形成單元剛度矩陣和進行節點應力磨平的時候，需要高斯積分。
以等參單元為例，其剛度矩陣
，這個就需要數值積分來快速計算，高斯點坐標及權系數如表4.2[王勖成]所示。
老師授課時一般對常應力單元進行推導，而常應力單元只有一個積分點，被積函數是常數，因此體現不出高斯積分來。很多老師對高斯積分在單元剛度矩陣的應用不予細述，導致部分同學對單元積分點認識不足。
3.單元應力指的是高斯積分點的應力，而非節點上的應力。有了位移模式，再通過虛功原理得到單元剛度矩陣，然后聚合總剛，求解平衡方程，就會把基本未知量——節點位移求出來了。通過節點位移得到單元應變結果，利用物理方程求得單元應力結果。
在等參元中，單元中n+1階(n=p-m)高斯積分點上的應變或應力近似解比其它部位具有較高的精度，因此我們稱(n+1)階高斯積分點是等參元中的最佳應力點。
展開 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中獲取積分點坐標的三種方法
經常有小伙伴問獲取積分點坐標的方法，今天給大家介紹三種獲取積分點坐標的方式，希望能給你們帶來幫助。
1 通過abaqus子程序獲取積分點坐標
Abaqus一些子程序中可以直接獲取積分點坐標，例如我們熟知的UMAT子程序中包含COORD參數，即為積分點坐標。順帶一提的是，當打開了幾何非線性時，該積分點是當前構形下的坐標，如果未打開幾何非線性則為初始坐標。
2通過history output輸出積分點坐標
Abaqus可以直接在歷程變量history output中輸出積分點坐標。直接在history output中勾選COORD選項，但是這里需要注意的是，Domain中的Set集合如果是node set，這里輸出來的是節點坐標，當這里是element set的時候，輸出來的才是積分點坐標。
3通過等參單元映射函數計算
等參元中，為了方便計算，把整體坐標映射到自然坐標，然后在自然坐標下進行高斯積分。如果知道了自然坐標下的高斯積分點，通過映射函數反算，便能得到整體坐標下的高斯積分點坐標。以四邊形等參單元為例，其以自然呢坐標為變量的插值形函數如下
坐標變換采取同樣的插值函數(叫做等參的原因)，整體坐標和自然坐標的關系式如下，如果知道自然坐標下的高斯積分點，直接通過此公式計算其在整體坐標下的坐標。
展開 圓形銅柱Taylor沖擊測試仿真的EFG算法實現
(1)根據分析對象的形狀布置節點和邊界，并將分析對象所在空間劃分背景積分網格,在背景網格內布置高斯積分點，記錄下高斯積分點的位置和積分參數。
(2)在每步計算中,完成以下步驟:
1)掃描高斯積分點，如果高斯積分點在分析域內，則執行步驟2) ~4)。
2)根據高斯積分點的坐標以及該高斯積分點的影響半徑，找到位于該高斯點影響域內的所有節點。
3)根據高斯點和影響域內節點的坐標進行最小二乘運算，得到形函數矩陣以及形函數矩陣的導數矩陣B。
4)計算并存儲該高斯積分點的剛度矩陣和載荷向量。
5)集成并求解整體剛度方程，得到位移向量u'
6)根據1'，由最小二乘法擬合出節點的真正位移u。
7)計算積分點的應變、應力等。
(3)更新節點位置和邊界位置，重復執行步驟(2)，直至計算結束。
圖1 EFG算法實現流程圖
3圓形銅柱Taylor沖擊計算模型
3.1基本建模流程
在WB LSDYNA中建立圓柱的幾何模型(WB LSDYNA只在ANSYS19.0版本及以上才有，本文選用的版本為ANSYS19.0)，圓柱底面半徑設為5cm，圓柱高度設為50cm。圓柱材質為銅，銅柱材料模型為*MAT PLASTIC KINEMATIC，銅柱的材料參數見3.2關鍵字所示。約束及初速度設置：銅柱除底部節點設置為全約束外（關鍵字*SPC_SET），所有節點具有負Z方向初始速度200m/s。計算時間為0.9ms，整個沖擊測試模型采用EFG算法（*Solid_EFG）進行計算，本文建模不考慮熱力學因素，故建模單位采用國際單位制kg-m-s，計算模型如圖2所示。
展開 ANSYS Workbench 接觸高級選項詳解（1）
3.Detection Method
這個是用來選擇接觸點之間的探測方法，ANSYS workbench 提供了四種探測方法：
1) On Gauss Point（高斯積分點）
我們知道，有限元計算時，首先求解總體方程矩陣得到節點位移，然后通過幾何方程和物理方程求解積分點上的應變和應力，最后通過形函數計算得到節點上的應力。
因此積分點上的應變和應力相對于節點上是更加準確的。
高斯積分點探測法
在面對面的接觸計算時，相比節點探測法來說，使用高斯積分點探測方法可以得到更加準確的計算結果。
在面對面的接觸中，如使用罰函數接觸算法或增廣拉格朗日算法時，Workbench默認的探測方法是高斯積分點法。
如上圖所示由于這種探測方法會造成接觸面之間的穿透，因此不適用于MPC 算法和一般拉格朗日算法。
2）Nodal-Normal From Contact
探測與接觸面垂直的接觸法線上的節點
3）Nodal-Normal to Target
探測與目標面垂直的接觸法線上的節點，這種方法是在使用MPC 和一般拉格朗日算法時采取的探測方法
4）Nodal - Projected Normal From Contact（投影法）
探測穿透接觸面和目標面的節點，這個一般不常用，不做過多介紹。
在經典版ANSYS 中提供了更多的探測方法，不過對于工程應用來說，Workbench提供的這幾種方法也足夠了。對于新手來說，使用Workbench計算最佳的方式就是先使用默認的接觸設置先計算，然后在通過計算結果和計算時間，根據自己的需求，調整設置，每次最好只改動一種設置嘗試。
展開 有限元中的高斯點與積分詳解(下)_《數值計算與程序設計》系列課程之九 ￥599
本課從實際問題出發，帶著問題去講解有限元中的高斯點與數值積分。一開始拋出了以下3個關鍵問題：
1.對于一個任意函數怎么去得到它的積分？
2.數值積分的本質是什么？為什么簡單地取幾個點就可得到積分值？此種方法的立足點在哪？
3.很多資料上都說“有限元求解精度嚴重依賴于網格質量，過度扭曲的單元會導致結果不收斂或者精度極度惡化”，這只是為什么呢？扭曲單元到底影響的是有限元方法中的哪一步？
圍繞這3個問題，本課分別講了一下三個內容：
1. 數值積分基本方法。
2. 有限元單元積分。
3. 誤差分析。
希望有興趣的同學多多支持下，你們的支持是我更新的動力
展開 
基于abaqus的有限元理論詳解
這種積分方案對于提高有限元位移解的精度是有益的。因為減縮積分較精確積分所得的積分值偏小。而位移有限單元法中位移函數用有限自由度來逼近無限自由度，使單元的剛度值擴大了。從而上述兩種因素引起的誤差被部分的抵消了。
精確積分通常由形函數中非完全多項式的最高階次所要求，而決定有限元精度的通常是完全多項式的階次。非完全的高次項往往不能提高精度，反而帶來不利的影響，也就是說，積分點的選擇只要能保證形函數中完全多項式部分的精確積分就可以了。不會因積分誤差帶來對有限元計算精度的影響。非完全的高次項在積分時得不到保證相當于對原位移函數作了調整，改善了單元分析精度。這種采用減縮積分保證完全多項式的積分精度來選擇積分點的積分方案稱為優化積分方案。
為使求解方程組K*d=R成為可能，引入邊界條件后K必須是非奇異的。K非奇異的條件是K的行列式不等于0。
若采用精確積分方案計算K，則其非奇異性要求總能得到滿足，因為任何非剛體位移模式對應的精確應變能總是大于0的。其剛度矩陣K必然正定。然而，當采用減縮積分時，K的非奇異性并不是必然的。例如在模型情況下，對應于某種非剛體位移模式，減縮積分時高斯積分點上的應變正好等于0，此時的應變能為0，這種非剛體位移模式稱為零能量模式。因此采用減縮積分時，要檢查K的非奇異性。一般一階減縮積分在某些特殊情況下會出現剛度矩陣的奇異。而二階減縮積分不會出現上述情況。也即是二階減縮積分抵抗沙漏的能力要強于一階減縮積分。
總之，通常有限元網格中的單元均是等參單元，要和標準單元之間進行坐標變換。這使得計算單元剛度矩陣是要使用高斯積分，不同階次的單元，位移函數的階次不同，精確積分所需的積分點個數也不同。剛度矩陣是通過高斯積分（單元上積分點上的函數值乘以權系數）得到的。求解得到的位移是節點的。
展開 四節點/八節點四邊形單元懸臂梁的Matlab有限元編程——《Matlab有限元編程從入門到精通》系列
(21)
4、高斯積分
公式（20）中的單元剛度矩陣通過數值積分求得，本案例中的四節點和八節點四邊形等參單元均采用2*2個積分點的高斯積分即可求得精確結果。高斯積分點的坐標具體如圖所示。
ABAQUS中的單元選擇-理解剪切自鎖和沙漏
通常ABAQUS會默認采用縮減積分單元（Reduced Integration），本文主要介紹相關概念。以下內容參考了《塑性非線性分析原理》等書籍。
圖1
一、完全積分Full Integration和縮減積分Reduced Integration的區別
根據有限元理論，對平面4節點單元，如果在一個方向上采用2個高斯積分點，其所計算的剛度矩陣是精確的。因此，4節點4邊形的2*2=4個高斯積分點方案為完全積分方案，而縮減積分方案（4節點4邊形采用1*1個高斯積分點）不能對剛度矩陣進行精確積分。（圖2）
圖2
二、完全積分方案中的剪切自鎖shear locking現象
當模擬具有彎曲特征的問題時，如果采用完全積分方案單元（如平面4節點應變單元CPE4，平面4節點應力單元CPS4，三維9節點單元C3D8等）進行網格劃分，會出現剪切自鎖的現象，此時加密網格也不能解決問題。為了說明問題，考察如圖3所示的純彎梁。由于平面4節點單元采用線性的位移插值函數，因此變形后的網格必為一梯形。變形后單元的內角不再保持為直角，也即產生了剪應變，這與實際情況不同。換句話說，這類單元在剪切時過于剛硬，與實際情況不同。由于采用了完全積分方案，剛度矩陣被精確積分，剪切自鎖現象將對計算結果產生影響。即使加密網格，計算結果也不會有明顯的改善。讀者可以自行對一懸臂梁算例進行驗證。
需要指出如果采用平面8節點單元，如CPE8，位移在邊上2次分布，可以反映彎曲特征，不會產生剪切自鎖現象。
圖3
三、縮減積分方案中的沙漏現象Hourglass
如果采用縮減積分方案，平面4節點單元只在單元中心點設置一個積分點，該點變形后的剪切應變為0，不會出現剪切自鎖現象。
展開 有限元計算的節點解與單元解續
數值計算時的積分方法有Newton-Cotes積分和高斯積分，兩種積分方法都有其積分點位置，對一個固定的單元有其固定的積分點，并且高斯積分相對來說達到的精度更高。
在書中經過驗證，高斯積分點即是單元的最佳應力點。
采用位移場得到的位移解在全局都是連續的，應變和應力在單元內部是連續的，而在單元之間則一般不連續，也即在單元邊界上發生跳突，因此同一個節點，由圍繞它的不同單元計算得到的應變值和應力值通常不同。
圖1 常見的平面單元最佳應力點位置
為了得到連續的應力場，需要做一些處理，處理的方法有好幾種：
（1）單元平均
這種方法比較簡單，取節點的相鄰單元應力做平均值：
平均應力=（單元1應力+單元2應力）/2
也可以采用精確一點的面積加權平均，也即是將每個單元的應力乘以面積后相加，最后除以所有相鄰單元的面積之和。
（2）單元平均
這種處理方法將圍繞某個節點的所有單元在這個節點處的應力值加和求平均。
（3）總體應力磨平
總體應力磨平方法是構造一個改進的應力解，此改進的解在全局是連續的，改進后的解與有限元法求得的應力解應滿足加權最小二乘的原則，對于由4個單元組成的網格，總體應力磨平后示意圖如圖2.
圖2 總體應力磨平示意圖
這種方法的主要缺點是計算量很大，基本上相當于做了兩次有限元計算。
（4）單元應力磨平
這是為減少改進應力結果的工作量而采用的方法。單元磨平方法的主要缺點是仍然未能充分利用單元的最佳應力點具有高一階精度的特性。
（5）分片應力磨平
該方法是能夠有效解決總體應力磨平和單元應力磨平的缺點。
展開 