四節點平面應變單元的案例
2階8節點減縮積分平面應變單元子程序UELMAT ￥1
2階8節點減縮積分平面應變單元子程序UELMAT源代碼及計算算例
基于ABAQUS的UEL子程序定義4節點平面應變等參單元的剛度問題
摘要：
采用基于ABAQUS的UEL子程序開發4節點平面應變等參單元，采用雙線性形函數，4點高斯積分，本構關系為線彈性各向同性材料，得到的單元剛度矩陣和ABABUS自帶的CPE4單元的單元剛度矩陣（剛度矩陣輸出方式為*element matrix output, elset= ALLE, stiffness=yes, OUTPUT FILE=USER DEFINED）不同；對比ANSYS的單元剛度矩陣，結果顯示兩者也不相同。問題出在哪里呢？本文檔將對此問題進行回答。
本文可以作為ABAQUS高級子程序UEL的入門級教程，做UEL的應該關注下！
基于ABAQUS的UEL子程序定義4節點平面應變等參單元的剛度問題（技術鄰 藍牙）.pdf
展開 平面四節點等參單元程序
用C語言編寫的程序，計算結果寫成兩個文件，一個輸出節點位移和單元應力，另一個可以用tecplot打開，看位移云圖和應力云圖
平面四節點等參單元.part3.rar
平面四節點等參單元.part1.rar
平面四節點等參單元.part2.rar
有限元2D單元妙用 平面應力與平面應變 廣義平面應變 硬干涉 ￥10
平面應力單元還可以跟軸對稱單元結合，模擬出變厚度模型。比如對葉盤的分析。需要注意的是，在ANSYS里面，當我們將平面應力和軸對稱單元結合的時候，平面應力單元的厚度應該設置為所有圓周分布葉片厚度的總和。如下圖。
平面應變單元：
平面四邊形四節點單元計算程序與ANSYS結果對比
為什么要導出單元剛度矩陣
在學習有限元方法時，我們會需要編寫程序計算結構的單元剛度矩陣。此外，當我們需要做有限元軟件二次開發時，我們也需要驗證所做的開發是否正確。為了驗證程序正確性，我們可以從商業有限元軟件中導出單元剛度矩陣來驗證程序的計算結果。下面簡單介紹從ansys軟件中導出平面四邊形四節點單元的單元剛度矩陣。
平面四邊形四節點單元示例
如圖所示，計算這兩個單元組成單元剛度矩陣，并組裝成整體剛度矩陣，求解各個節點的位移。
展開 有限元模擬條形基礎持力層，abaqus 地基初始應力場設置（二）
有限元模擬條形基礎持力層
模型概況
基礎形式：條形基礎
基底摩擦條件：完全粗糙
荷載情況：基礎承受豎向荷載
模擬的目標
1、地基初始應力狀態
2、條形基礎持力層在極限狀態的位移場
3、地基極限承載力
模型的注意事項
1、 基礎簡化為剛性基礎
2、 該問題簡化為平面應變問題，采用 CPE4 四節點平面應變單元
3、 基底“完全粗糙”在模型中的體現：約束基底范圍的水平位移
4、 彈性階段采用線彈性本構模型，塑性階段采用莫爾-庫倫本構模型
5、 基礎埋深范圍內的土層的重力以等效荷載替代
6、 獲取基礎持力層“荷載沉降曲線”的方法：指定基底范圍的沉降，沉降值要足夠大，確保持力層進入極限破壞狀態。通過給定的沉降求解基底范圍節點的豎向約束力。
有限元模型
在加載分析步中，指定基底范圍的沉降為 y 方向 -0.3m（見下圖），以此確保地基達到破壞狀態。
基礎持力層極限狀態下的位移場
通過 YZ 平面鏡像，得到左部分的位移場。
基底壓力與沉降關系曲線
從關系曲線拐點處可以得到基礎持力層的極限承載力：320.7 kPa。
地基初始應力場設置
本例的地基初始應力場是由自重、基礎埋深范圍內土體等效荷載產生的。
在 Initial 初始步中定義 Geostatic stress ,分別設置地基模型頂面、底面的豎向坐標和對應的應力（如下圖），ABAQUS 會根據兩端的應力進行線性插值構建應力場。
此外，要保證初始應力場的平衡，需要在 geostatic 分析步正確施加地基的重度（Body force）、外荷載（基礎埋深范圍土體的等效荷載）。
展開 擴展黃永剛原始晶體塑性程序加入AF背應力模擬金屬疲勞問題
參考文獻：《Low-cycle fatigue life prediction of a polycrystalline nickel-base superalloy using crystal plasticity modelling approach》
在原始程序中修改流動方程，加入背應力項，引入運動硬化項，從而可以描述多晶金屬循環加載中的包辛格效應
背應力的演化遵循AF模型
并使用原始的PAN模型描述滑移系統的硬化行為
為了表征多晶的疲勞壽命，引入兩類疲勞指示因子分別為
一：累計塑性滑移
二：累計能量耗散
以文獻的例，驗證修改模型的準確性，其中文獻作者的幾何模型和材料參數如下
依據該模型，作者模擬得到單調拉伸以及循環加載下材料的宏觀應力應變響應為
微觀響應結果為
基于兩類疲勞指示因子，作者通過線性外推得到了基于模擬的壽命預測結果:
基于作者提供的思路和參數，對黃永剛原始程序進行修改，考慮背應力效應，并進行簡單的數值驗證
1，建立包含200晶粒的二維多晶模型（0.1*0.03mm），并使用四節點平面應變單元進行網格劃分，如下圖
2，施加正弦形式的循環拉壓的位移載荷（1%），引力比為-1
3，模擬結果如下：
第一個滑移系統的背應力：
累計塑性剪切：
累計能量耗散：
宏觀應力應變響應：
展開 ABAQUS UEL - 損傷材料本構簡單應用于4節點平面單元 ￥300
利用ABAQUS自定義單元子程序，既可以開發新的單元，同時也可以定義新的材料本構模型。本文以損傷模型簡單應用于4節點平面單元為案例，介紹ABAQUS UEL的開發和使用。
如上圖所示，該單元包含4個節點，每個節點有兩個自由度，分別在水平（X）和垂直（Y）方向運動。節點1的兩個自由度被固定，節點4的水平自由度被固定，節點2的垂直自由度被固定。節點3和節點4在垂直方向上向上運動，位移為0.1mm。該正方形單元的邊長為100mm。在input文件里，坐標表示為，
定義節點組合與邊界條件為，
為了讓模型收斂性更好，采用quasi-newton 求解器。時間步設置為，
在文件夾中通過Powershell提交job和子程序，
單個單元的變形為，
采用不同的 ??
，在后處理中得到損傷因子的變化，
相對應的力-時間關系為，
對于多個單元的情況，比如9單元組成的模型，
具體介紹見知乎：ABAQUS UEL - 損傷材料本構簡單應用于4節點平面單元 - 知乎 (zhihu.com)
相對應的UEL代碼和input文件在付費內容中，
展開 平面四節點靜力分析程序源碼
平面四節點靜力分析程序源碼.rar
空間8節點的單元剛度的計算源程序.rar
平面應變單元CPE4R齒輪傳動接觸應力計算 ￥49.9
厚齒輪的應力符合平面應變狀態，可以采用平面應變單元CPE4R來進行快速接觸應力計算。
在sketch模塊建立非對稱結構齒輪的草圖，然后建立part，并在assembly模塊進行裝配。
非對稱齒輪草圖
齒輪裝配體
通過適當的結構設計，非對稱齒輪可以在定速轉動的情況下獲得按某規律的變化轉速，在工程上經常會用到。
非對稱齒輪傳動分析結果
非對稱齒輪應力云圖
非對稱齒輪齒合區域局部應力云圖
三角形常應變單元解平面問題實施步驟與注意事項
利用上面討論的三角形常應變單元解平面問題，其具體步驟可歸納如下:
1)將要計算的彈性體劃分成三角形單元。對結點進行編號，列出結點坐標作為輸入信息。
(2)對單元進行編號，列出單元三個結點的號碼作為輸入信息。
(3)計算載荷的等效結點力，把等效結點力作為輸入信息。
(4)按照(6)式計算各單元的常數bi、ci、bj、cj、bm、cm，再按照(4)計算2A。
(5)按照(35)式計算各單元的剛度矩陣。
(6)形成整體剛度矩陣。
(7)處理約束及消除剛體位移。
(8)解線性方程組(32)式，求結點位移。
(9)按照(20)式計算應力矩陣，再按(18)式計算單元應力。根據需要計算主應力和主方向。
通常步驟(4)至(9)均由計算機來完成，而步驟(1)至(3)可以用手工完成，也可由計算機來完成。在實現以上各步驟時，為了達到一定的計算精度，節約計算機存儲量，縮短計算機運行時間等目的，還需要注意下列事項。
1、利用對稱性
在劃分單元前要研究一下，計算對象是否有對稱變形或反對稱變形存在，從而確定是否需要取整個物體，還是取部分物體作為計算模型。例如圖8a所示受純彎曲的梁，它對于x,y軸都對稱，而載荷對于y軸對稱，對于x軸反對稱。可見，應力和應變亦將具有同樣的對稱和反對稱特性，所以我們只需計算1/4梁就行了。分離體如圖8b所示。對于刪去部分結構的影響可以這樣考慮：對于處于y軸對稱面內各結點的x方向位移和y方向分布力都應等于零，而對于處在x軸反對稱面上的各結點的x方向位移和y方向分布力亦都應等于零。這些條件相當于安置如圖8b中的約束。圖中o點上安置y方向的約束是為了消除剛體位移而設置的。又例如在分析圖9中所承受均勻壓力的厚壁圓筒時，根據結構和載荷軸對稱的性質，我們可以取出一個小扇形(圖中陰影部分)進行計算。
展開 
有限元模擬臨坡地基，abaqus 從外部導入初始應力場（三）
其位移方向不確定，需要在模型中包含基礎
2、 該問題簡化為平面應變問題，采用 CPE4 四節點平面應變單元
3、 基底“完全粗糙”，在模型中采用簡化的方式實現：基底與地基共節點
4、 地基土體在彈性階段采用線彈性本構模型，塑性階段采用莫爾-庫倫本構模型
5、 為確保基礎持力層達到極限狀態，在基礎頂面中點施加足夠大的豎向集中荷載
有限元模型
在加載分析步中，在基礎頂面施加集中荷載 2000kN（等效為每延米2000kN）（見下圖），以此確保地基達到破壞狀態。
臨坡基礎持力層極限狀態下的位移場、PEMAG云圖
分析步時間與基底豎向位移的關系曲線
曲線在分析步時間為 0.74 時，發生明顯轉折，可以判斷基礎持力層達到極限破壞的狀態。
對應的地基極限承載力為 0.74 × 2000 = 1480 kN/m。
臨坡地基初始應力場設置
本例中，由于地基邊坡的地表不是水平的，所以初始應力場不能采用“有限元模擬條形基礎持力層，abaqus 地基初始應力場設置（二）”中的方法設置。
需要采用從外部導入的方法設置初始應力場（設置過程復雜，后面將會制作視頻進行介紹）。以下只做簡單介紹（如急需知道詳細操作請留下郵箱，我會發送詳細教程）。
1、在建模的最初階段，把所有材料都按地基土設置，先不添加塑性部分的本構模型。計算臨坡地基在自重作用下的應力分布，如下圖所示。并將計算結果數據庫文件做好備份，后面用于導入。
2、在模型中設置各種材料的參數，添加地基土的塑性本構模型。添加 Geostatic 分析步，然后分別修改“輸出請求控制”、“荷載生效時間”，如下圖。
3、從第 1 步備份的數據庫文件導入初始應力場，如下圖。
最后，檢查應力分布云圖、位移分布云圖，確認初始應力場是否設置成功。
展開 ABAQUS UEL-梯度損傷模型應用于4節點平面單元 ￥600
傳統損傷模型對于單元的尺寸十分敏感，不同單元尺寸會導致有限元模型精度出現明顯偏差。針對該問題，梯度損傷（Gradient-damage）模型的概念被提了出來。
本文詳細介紹了如何將梯度損傷模型應用于4節點平面單元，并在有限元模型中進行模擬。
ABAQUS提供了UEL（user defined element）給使用者進行開發。筆者利用UEL開發4節點平面單元，其邊界條件如下圖所示。其中，節點1的X、Y方向被限制住，節點2的Y方向被限制，節點4的X方向被限制，節點3、4的Y方向有豎向位移0.1mm。單元為100*100mm的二維正方形。
每個節點除了X和Y方向的位移，還帶有非局部應變（nonlocal strain）。
單個單元模型，
多個單元模型，
具體內容可參見知乎文章：
ABAQUS UEL-梯度損傷模型應用于4節點平面單元 - 知乎 (zhihu.com)
相應的input文件和uel代碼付費可見，
展開 晶體塑性每日文章推薦（七）
然而，晶粒取向對宏觀應力-應變曲線的影響很小
基于作者提供的思路和參數，對黃永剛原始程序進行修改，考慮背應力效應，并進行簡單的數值驗證
1，建立包含200晶粒的二維多晶模型（0.1*0.03mm），并使用四節點平面應變單元進行網格劃分，如下圖
2，施加正弦形式的循環拉壓的位移載荷（1%），引力比為-1
3，模擬結果如下：
第一個滑移系統的背應力：
累計塑性剪切：
累計能量耗散：
宏觀應力應變響應：
整體上于作者的模擬情況比較接近，做疲勞分析的同學可以進行類似嘗試
平面四節點靜力分析程序源碼
平面四節點靜力分析程序源碼.rar
空間8節點的單元剛度的計算源程序.rar
