ABAQUS剛度的案例
(一維彈簧單元的直接剛度方法)Python編程和ABAQUS結(jié)果對比
這個在我去年的時候就已經(jīng)發(fā)過帖子了，但是沒有講解代碼，沒有和ABAQUS有限元軟件對比。這次發(fā)帖子就是想講解代碼以及和ABAQUS結(jié)果對比，發(fā)現(xiàn)我又更深層次理解了有限元方法。
這次寫得是最簡單的模型：一維線彈簧單元。采用的是直接剛度法求解。
直接剛度法的求解思路如下所示，其中粉紅色的是輸入，淡藍(lán)色是輸出。主要是要區(qū)分齊次邊界條件和非齊次邊界條件，非齊次邊界條件的話就要修改【F】。
下面將貼出我用python寫得一維彈簧單元的直接剛度法：
例子計算：
如圖是一個彈簧系統(tǒng)，單元節(jié)點信息如下，5節(jié)點受到一個強制位移20mm，明顯這是一個非齊次邊界條件問題。
Python編程輸入信息如下：
結(jié)果如下：
可以看到，輸出結(jié)果和書上的答案一致。
下面進行ABAQUS模擬：
添加彈簧單元
添加邊界條件
顯示單元編號、節(jié)點編號如下所示，紅色表示單元編號，黃色是節(jié)點編號
ABAQUS計算結(jié)果如下：
首先是變形圖前后對比
反力云圖如下所示，基本和直接剛度法計算的結(jié)果一致
位移云圖如下所示，基本和直接剛度法計算的結(jié)果一致
整體剛度矩陣如下所示，因為ABAQUS彈簧單元是三維的，每一個節(jié)點有3個自由度，15x15，原味的剛度矩陣如下
我們把剛度矩陣轉(zhuǎn)化為一維的，方便和編程的結(jié)果對比
從結(jié)果可以看出，ABAQUS的整體剛度矩陣和直接剛度法計算出來的整體剛度矩陣有些差異，如圖標(biāo)紅的所示。
那么在整體剛度矩陣上為什么ABAQUS會和直接剛度法的整體剛度矩陣有差異呢，到底ABAQUS的整體剛度矩陣對不對呢？答案將在下一期揭曉。歡迎大家積極討論。
展開 Abaqus一鍵輸出剛度矩陣和質(zhì)量矩陣的插件EMM ￥20
=======概述=======
EMM(Export Matlab Matrix)是集成在ABAQUS/CAE中的一個插件，能夠一鍵輸出Abaqus模型的單元及全局剛度、質(zhì)量、載荷矩陣，并自動轉(zhuǎn)換為MATLAB矩陣。
Abaqus可以輸出剛度、質(zhì)量等單元或者全局矩陣。但需要手動添加關(guān)鍵詞，較為麻煩，且輸出的.mtx矩陣格式比較亂，不易閱讀。
本工具能一鍵實現(xiàn)以下功能：
1. 輸出Abaqus模型的剛度、質(zhì)量、載荷矩陣到.mtx，包括單元和全局的；
2. 將上述矩陣轉(zhuǎn)換為MATLAB的.mat文件；
3. 用Matlab打開上述.mat文件。
如果用戶有任何問題或者需要合作，歡迎聯(lián)系我們。email: SnowWave02@qq.com
=======使用過程=======
打開Abaqus模型，切換到Job模塊，在Abaqus菜單欄的Plug-ins里看到iSolver插件的菜單。
點擊iSolver->Export Matlab Matrix…，彈出EMM界面。
勾選需要輸出的矩陣和Set，點擊EMM界面右側(cè)的Submit后將提交Abaqus計算，計算完畢得到如下類似矩陣。
=====演示視頻======
http://www.yqgqt.org.cn/college/video/c11343
=====下載======
用到的模型
Job-Cabin-T.rar
用戶手冊
Export Matlab Matrix工具說明(SnowWave02 20191217).pdf
EMM插件
展開 ABAQUS 剛度矩陣 ￥1000
我有個abaqus的問題，你們幫幫我出出主意。 是這樣的：
1，TestElement.py 是編寫的測試代碼，可輸出8結(jié)點線性單元的剛度矩陣。
2，abaqus文件Job-testing.inp ，運行可以輸出單元剛度矩陣。
問題是：他們的結(jié)點坐標(biāo)/排序，材料參數(shù)都一致， 但得到的剛度矩陣就不一樣。
需要：代碼輸出的剛度矩陣與abaqus得到的剛度矩陣一致或者基本一致。
TransferMatrix：導(dǎo)出Abaqus剛度矩陣
本期給大家推薦一款由木木自研的小工具：TransferMatrix，主要用于導(dǎo)出 Abaqus 中各種剛度矩陣，方便和自研的程序進行實時對比，基于 Pyside6 搭建的軟件界面，用戶僅需導(dǎo)入 inp 文件，就可以自動調(diào)用電腦內(nèi)的 Abaqus 進行計算分析，不受限于 Abaqus 版本，可以導(dǎo)出：
單元剛度矩陣
單元質(zhì)量矩陣
單元分布節(jié)點荷載列陣
整體剛度矩陣
整體質(zhì)量矩陣
可選是否導(dǎo)出后打開 Matlab ，并保存為 mat 文件
默認(rèn)保存的格式為 csv，用戶可以直接在TransferMatrix 中查看
可選是否保存整體剛度/質(zhì)量矩陣為稀疏矩陣/全矩陣形式
軟件界面見下圖：
若勾選導(dǎo)出后打開matlab，可直接查看：
軟件實現(xiàn)原理
使用注意事項
軟件啟動速度較慢，需要有點耐心，鄙人對于 PySide 的打包技術(shù)認(rèn)識有限
若 inp 文件內(nèi)沒有密度參數(shù)，是不能導(dǎo)出質(zhì)量矩陣的，建議在動力類型的分析步中進行導(dǎo)出質(zhì)量矩陣
非協(xié)調(diào)單元和雜交單元只能輸出質(zhì)量矩陣
在輸出載荷列陣時，Abaqus 產(chǎn)生的是單元節(jié)點的分布載荷，如果沒有這中載荷類型，也是無法導(dǎo)出的
建議再導(dǎo)入 inp 文件之前，確保這個文件能在 Abaqus 上跑通，本軟件只是負(fù)責(zé)轉(zhuǎn)換剛度矩陣
不支持在TransferMatrix 中查看整體剛度/質(zhì)量矩陣，因為一般情況下全局矩陣尺寸很大，可選擇導(dǎo)出至 Matlab 中進行查看
如果出現(xiàn)風(fēng)險提示，請選擇信任：
軟件同級目錄放置的config.json文件，里面配置的是本地的 abaqus.bat 和 matlab.exe 目錄，建議用 everything 工具搜索出這兩個的安裝位置，復(fù)制進去即可。
展開 
?abaqus提取單元的剛度矩陣和質(zhì)量矩陣
abaqus提取單元的剛度矩陣和質(zhì)量矩陣：
*Output, history, variable=PRESELECT,
*File Format,ASCII
*Element Matrix Output,Elset=Beam-1.Set-2,
File Name=shuchu,Frequency=1,Output File=User Defined,Stiffness=Yes
*End Step
*Step
*Matrix generate，stiffness
*End Step
有限元理論基礎(chǔ)及Abaqus內(nèi)部實現(xiàn)方式研究系列6：General梁單元剛度矩陣 ￥1
本文首先簡單介紹梁單元的基本理論，分析了每一部分剛度的來源，并研究了Abaqus中General梁的B31單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus完全一致的剛度矩陣，自研軟件的正確性證明是非常困難的，無論你測試多少標(biāo)準(zhǔn)算例或者和試驗結(jié)果對比，都很難讓用戶踏實的相信你的結(jié)果和商業(yè)軟件是一致的，但如果一個模型僅有General梁組成，那么無論這個模型多么的復(fù)雜，我們一般都可以得到和Abaqus完全一致的分析結(jié)果。最后我們使用了一個簡單的算例綜合驗證，針對該算例我們對比采用了同樣修正的自編程序iSolver和Abaqus的結(jié)果，可以發(fā)現(xiàn)結(jié)果如預(yù)期的一樣，沒有任何誤差。
==演示視頻==
https://www.yqgqt.org.cn/college/video/c12884?chapter=1
梁的靜力分析結(jié)果校核視頻演示任意分段梁的靜力分析，證明iSolver的分析結(jié)果和Abaqus沒有誤差，有兩個算例：Generalized梁和L型梁。
==總結(jié)==
General梁的B31單元的剛度矩陣在Timoshenko梁理論基礎(chǔ)上的修正如下表：
項次
剛度
修正
不修正
說明
1
軸向拉伸剛度
√
2
橫向彎曲剛度
√
采用減縮積分
3
軸向拉伸和橫向彎曲耦合剛度
√
根據(jù)形心偏置進行修正
4
軸向扭轉(zhuǎn)剛度
√
剪切中心的偏置會影響抗扭剛度系數(shù)J
5
橫向剪切剛度
√
增加了一個幾何因子，使得細(xì)長梁的時候該項趨于0。
6
軸向扭轉(zhuǎn)和橫向剪切耦合剛度
√
根據(jù)剪切中心的偏置進行修正。
展開 有限元理論基礎(chǔ)及Abaqus內(nèi)部實現(xiàn)方式研究系列12： 幾何梁單元的剛度矩陣 ￥1
如果有任何其它疑問或者項目合作意向，也歡迎聯(lián)系我們：
snowwave02 From www.yqgqt.org.cn
email: snowwave02@qq.com
詳細(xì)研究方法，見附件：
有限元理論基礎(chǔ)及Abaqus內(nèi)部實現(xiàn)方式研究系列12：幾何梁單元的剛度矩陣.pdf
以往的系列文章：
第一篇：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é)果的正確性。
https://www.yqgqt.org.cn/content/post/373743
第六篇：General梁單元的剛度矩陣。介紹梁單元的基礎(chǔ)理論和Abaqus中General梁單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus梁單元完全一致的剛度矩陣。
展開 續(xù)集(一維彈簧單元的直接剛度方法)Python編程和ABAQUS結(jié)果對比
我們編程的結(jié)果如下所示
ABAQUS結(jié)果和編程結(jié)果對比
所以ABAQUS提取的整體剛度矩陣實際是經(jīng)過補償后的剛度矩陣，嚴(yán)格來說是不正確的，但是并不影響力和位移的結(jié)果。
正確的剛度矩陣如下
如果疑問，歡迎交流和指正。
補償法的部分代碼如下
有限元理論基礎(chǔ)及Abaqus內(nèi)部實現(xiàn)方式研究系列1：S4殼單元剛度矩陣研究
具體的剛度矩陣在Mindlin理論基礎(chǔ)上的修正如下表：
項次
剛度
修正情況
說明
修正
不修正
1
薄膜效應(yīng)剛度
√
Abaqus為了消除殼單元出現(xiàn)的剪切鎖死，計算薄膜效應(yīng)時把殼作為體來處理
2
面外彎曲剛度
√
3
面外橫向剪切剛度
√
（1）Abaqus為了處理薄板情況，增加了一個幾何因子
（2）為了消除沙漏問題增加相關(guān)剛度
4
第6自由度剛度
√
Abaqus取了和對角剛度相關(guān)的小量
5
其它非對角元素
√
（1）during analysis選項導(dǎo)致的小量
（2）第6自由度剛度相關(guān)小量
詳細(xì)研究方法，見附件：
有限元理論基礎(chǔ)及Abaqus內(nèi)部實現(xiàn)方式研究系列1：S4殼單元剛度矩陣研究(SnowWave02 20170708).pdf
展開 有限元理論基礎(chǔ)及Abaqus內(nèi)部實現(xiàn)方式研究系列22： 幾何非線性的剛度矩陣求解
我們關(guān)注CAE中的結(jié)構(gòu)有限元，所以主要選擇了商用結(jié)構(gòu)有限元軟件中文檔相對較完備的Abaqus來研究內(nèi)部實現(xiàn)方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數(shù)學(xué)上其實并不嚴(yán)謹(jǐn)，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。
自主結(jié)構(gòu)有限元求解器iSolver介紹視頻：
http://www.yqgqt.org.cn/college/video/c12884
==第22篇：幾何非線性的剛度矩陣求解==
幾何非線性在界面上是很容易設(shè)置的，但商軟內(nèi)部的處理相當(dāng)復(fù)雜，我們從最基本的剛度矩陣的求解出發(fā)，看看在幾何非線性設(shè)置后，剛度矩陣具體是怎么實現(xiàn)的。本文首先介紹幾何非線性下的剛度矩陣的理論推導(dǎo)和計算機求解方法，說明理想的求解方式的困難點和猜測Abaqus內(nèi)部的解決方法。最后利用一個簡單的算例通過對比iSolver和Abaqus的結(jié)果，部分驗證我們對Abaqus幾何非線性的剛度矩陣的實現(xiàn)方式的猜測。
1.1 幾何非線性的剛度矩陣推導(dǎo)理論
在前面17章：幾何非線性的物理含義中，我們提到如果是非線性系統(tǒng)，應(yīng)變能W隨t的變化就是個非線性過程。每個時刻點可以求出一個斜率，這個斜率最終會形成當(dāng)前時刻點的剛度矩陣。
求導(dǎo)后得到的剛度K：
也就是剛度矩陣將分為兩塊：
（1） 上式的前面一部分稱為材料剛度陣，依然是以前的BDB形式，只不過B換成了當(dāng)前時刻的應(yīng)變位移矩陣
（2） 后面新增項一般稱為幾何剛度陣，在Abaqus中稱為初始應(yīng)力矩陣（initial stress stiffness）。
展開 【ABAQUS模態(tài)動力學(xué)】Composite&abaqus 預(yù)應(yīng)力模態(tài)分析&輸出單元剛度矩陣
劃分網(wǎng)格
定義邊界條件
5.2 預(yù)應(yīng)力模態(tài)
預(yù)應(yīng)力模態(tài)，按我的理解就是，假設(shè)t=0時，結(jié)構(gòu)的剛度矩陣、質(zhì)量矩陣為M0,K0；t=t1時，結(jié)構(gòu)（分析對象）收到外部激勵的作用，使得結(jié)構(gòu)的剛度矩陣，質(zhì)量矩陣發(fā)生改變，結(jié)構(gòu)的剛度矩陣、質(zhì)量矩陣為M1,K1；t=t2時，開始進行特征值提取，此時求解的是t1狀態(tài)的結(jié)果。
從上面這個理解出發(fā)，ABAQUS預(yù)應(yīng)力模態(tài)只要在frequency分析步之前進行General,Static分析步,打開NLGeom選項(分析過程中剛度矩陣會不斷變化)。
提取單元剛度矩陣：
【ABAQUS 二次開發(fā)筆記】輸出單元剛度矩陣 - hayden_william - 博客園
以上均為我的一點理解，不一定完全正確，本文僅作為個人學(xué)習(xí)記錄之用，其他概不負(fù)責(zé)。
展開 
基于ABAQUS的鋼筋混凝土梁的剛度分析
ABAQUS 是一套功能強大的工程模擬的有限元軟件，其解決問題的范圍從相對簡單的線性分析到許多復(fù)雜的非線性問題。 ABAQUS 包括一個豐富的、可模擬任意幾何形狀的單元庫。并擁有各種類型的材料模型庫，可以模擬典型工程材料的性能，其中包括金屬、橡膠、高分子材料、復(fù)合材料、鋼筋混凝土、可壓縮超彈性泡沫材料以及土壤和巖石等地質(zhì)材料，作為通用的模擬工具， ABAQUS 除了能解決大量結(jié)構(gòu)（應(yīng)力 / 位移）問題，還可以模擬其他工程領(lǐng)域的許多問題，例如熱傳導(dǎo)、質(zhì)量擴散、熱電耦合分析、聲學(xué)分析、巖土力學(xué)分析（流體滲透 / 應(yīng)力耦合分析）及壓電介質(zhì)分析。
本文采用ABAQUS簡單搭建鋼筋混凝土結(jié)構(gòu)有限元模型，解決工程實際問題。本文對某鋼筋混凝土模型進行模擬，采用實體單元（混凝土）與Beam單元（鋼筋）相結(jié)合的方式進行模擬。首先導(dǎo)入模型，劃分實體網(wǎng)格和一維梁單元。建立材料和截面，分為實體和一維梁單元。然后將材料賦予個模型，進行網(wǎng)格劃分。然后加載，簡支梁模型，約束兩端下部，在上部施加均勻壓力。最后求解計算。
由結(jié)果處理可知，鋼筋混凝土梁的最大位移為2.55cm，可以達到模擬效果。ABQUS可以有效模擬鋼筋混凝土模型，對工程實際有很大的參考指導(dǎo)意義。
設(shè)備基本情況：i5 4核 ；耗時：10分鐘內(nèi)。
圖1 鋼筋混凝土梁CAD模型
圖2 線框模型，中間為鋼筋模型
圖3-4 網(wǎng)格模型
圖5 施加載荷示意圖
圖6 位移云圖
展開 基于ABAQUS的橡膠懸置膠合件剛度仿真計算
本文將以一個懸置膠合件仿真的實例講解一下如何利用ABAQUS來獲取其三個方向的靜態(tài)特性。所用膠合件的數(shù)模圖如圖1所示。其設(shè)計圖紙上標(biāo)注的三向剛度如表1所示，膠料硬度是邵氏50±5度。
圖1 膠合件結(jié)構(gòu)
表1設(shè)計要求
1、 網(wǎng)格劃分
采用HYPERMESH對圖一懸置進行網(wǎng)格劃分到的有限元模型如圖2所示。
2、材料設(shè)置
把劃分好的網(wǎng)格導(dǎo)入ABAQUS中，設(shè)置其材料參數(shù)，由于不同本構(gòu)模型對橡膠懸置膠合件剛度計算結(jié)果有一定的影響。結(jié)合何小靜，上官文斌發(fā)表的《橡膠隔振器靜態(tài)力－ 位移關(guān)系計算方法》一文的研究結(jié)果表明，Mooney-Rivlin 模型的計算精度最高，其相對誤差均小于10%，所以本文采用M-R模型進行計算。50度膠料的M-R材料常數(shù)C10=0.2969，C01=0.0584。
3、剛度求解
3.1求解X方向剛度
按表 1要求，做如下設(shè)置：在Z方向先預(yù)載8mm，再在X向加載500N。取值0~5.6mm,對X向靜剛度進行求解。
求得的力和位移關(guān)系見表2所示，用表中數(shù)據(jù)進行畫圖差值可得到圖3所示的X向靜剛度為38N/mm,與設(shè)計值非常接近，其變形云圖見圖4所示
表2 X向力和位移關(guān)系表
圖3 X向剛度差值結(jié)果
圖4 X向云變形圖
3.2求解Y方向剛度
按表 1要求，做如下設(shè)置：在Z方向先預(yù)載8mm，再在X向加載1000N。取值2~4mm,對Y向靜剛度進行求解。
求得的力和位移關(guān)系見表3所示，用表中數(shù)據(jù)進行畫圖差值可得如圖5到Y(jié)向靜剛度為98N/mm,與設(shè)計值80N/mm有一定差異，見圖4。其變形云圖見圖6.
展開 利用ABAQUS進行屈曲梁負(fù)剛度分析
計算結(jié)果
兩端均壓縮0.75mm后梁的變形為：
中部加載前梁的應(yīng)力分布云圖
加載后梁的應(yīng)力分布云圖：
將數(shù)據(jù)導(dǎo)出，在MATLAB中繪圖，得到力位移曲線：
由圖中數(shù)據(jù)計算可知梁的負(fù)剛度為：
按柔性設(shè)計手冊里的公式：
可見兩者相差極小，可認(rèn)為仿真結(jié)果有效。取上圖局部放大，得到壓力為0時梁的位移情況如下圖所示。
由于仿真計算采樣點的限制且誤差不可避免，通過上圖可知，位移為0時，受力也為接近0。
ABAQUS疑難雜癥之沙漏剛度
源文件可在公眾號:‘易木木響叮當(dāng)’ 內(nèi)回復(fù)“沙漏剛度”自動獲取。
