單元剛度矩陣的案例
單元剛度矩陣
由（27a）式定義的單元剛度矩陣，由于應變B對于3結點三角形單元是常量陣，因此有
代入彈性矩陣D和應變矩陣B后，它的任一分塊矩陣可表示成
由(34)式立即可以得到
由此可見單元剛度矩陣是對稱矩陣。
為了進一步理解單元剛度矩陣的物理意義，我們同樣可以利用最小位能原理建立一
個單元的平衡，這就得到
Pe是單元結點載荷，當然應當包括其它相鄰單元對該單元的作用力。現把ae、Pe順序表示為
這是單元結點平衡方程，每個結點在x和y方向上各有一個平衡方程，3個結點共有六個平衡方程。方程左端是通過單元結點位移表示的單元結點內力，方程右端是單元結點外載。
令a1=1(ui=1),a2=a3=…＝a6=0
由（38）式可以得到
（39）式表明，單元剛度矩陣第一列元素的物理意義是：a1＝1。其他結點位移都為零時，需要在單元各結點位移方向上施加結點力的大小。當然，單元在這些結點力作用下處于平衡，因此在x和y 方向上結點力之和為零：
對于單元剛度矩陣中其他列的元素也可用同樣的方法得到它們的物理解釋。因此單元剛度矩陣中任一元素Kij物理意義為：當單元的第j個結點位移為單位位移而其他結點位移為零時，需在單元第j個結點位移方向上施加的結點力的大小。單元剛性大，則使結點產生單位位移所需施加的結點力就大。因此單元剛度矩陣中的每個元素反映了單元剛性的大小，稱為剛度系數。
展開 有限元理論基礎及Abaqus內部實現方式研究系列12： 幾何梁單元的剛度矩陣 ￥1
http://www.yqgqt.org.cn/content/post/343905
第三篇：S4殼單元的剪切自鎖和沙漏控制。介紹Abaqus的S4單元如何來消除剪切自鎖以及S4R如何來抑制沙漏的。
http://www.yqgqt.org.cn/content/post/350865
第四篇：非線性問題的求解。介紹Abaqus在非線性分析中采用的數值計算的求解方法。
http://www.yqgqt.org.cn/content/post/360565
第五篇：單元正確性驗證。介紹有限元單元正確性的驗證方法，通過多個實例比較自研結構求解器程序iSolver與Abaqus的分析結果，從而說明整個正確性驗證的過程和iSolver結果的正確性。
https://www.yqgqt.org.cn/content/post/373743
第六篇：General梁單元的剛度矩陣。介紹梁單元的基礎理論和Abaqus中General梁單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus梁單元完全一致的剛度矩陣。
https://www.yqgqt.org.cn/content/post/403932
第七篇：C3D8六面體單元的剛度矩陣。介紹六面體單元的基礎理論和Abaqus中C3D8R六面體單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus六面體單元完全一致的剛度矩陣。
https://www.yqgqt.org.cn/content/post/430177
第八篇：UMAT用戶子程序開發步驟。介紹基于Fortran和Matlab兩種方式的Abaqus的UMAT的開發步驟，對比發現開發步驟基本相同，同時采用Matlab更加高效和靈活。
展開 有限元理論基礎及Abaqus內部實現方式研究系列6：General梁單元剛度矩陣 ￥1
==總結==
General梁的B31單元的剛度矩陣在Timoshenko梁理論基礎上的修正如下表：
項次
剛度
修正
不修正
說明
1
軸向拉伸剛度
√
2
橫向彎曲剛度
√
采用減縮積分
3
軸向拉伸和橫向彎曲耦合剛度
√
根據形心偏置進行修正
4
軸向扭轉剛度
√
剪切中心的偏置會影響抗扭剛度系數J
5
橫向剪切剛度
√
增加了一個幾何因子，使得細長梁的時候該項趨于0。
6
軸向扭轉和橫向剪切耦合剛度
√
根據剪切中心的偏置進行修正。
7
其它元素
√
都為0
詳細研究方法，見附件：
有限元理論基礎及Abaqus內部實現方式研究系列6：General梁單元剛度矩陣(SnowWave02 20180830).pdf
==以往的系列文章==
第一篇：S4殼單元剛度矩陣研究。介紹Abaqus的S4剛度矩陣在普通厚殼理論上的修正。
http://www.yqgqt.org.cn/content/post/338859
第二篇：S4殼單元質量矩陣研究。介紹Abaqus的S4和Nastran的Quad4單元的質量矩陣。
http://www.yqgqt.org.cn/content/post/343905
第三篇：S4殼單元的剪切自鎖和沙漏控制。介紹Abaqus的S4單元如何來消除剪切自鎖以及S4R如何來抑制沙漏的。
http://www.yqgqt.org.cn/content/post/350865
第四篇：非線性問題的求解。
展開 一個單元也能干大事之單元剛度初探
本公眾號將推出“一個單元也能干大事”系列文章，試圖從最小的單元出發，來為初學者講解一些有限元基本理論。該系列文章力爭篇幅短小，簡單通俗，深入本質，一看就懂，并且希望能跟讀者討論互動。
一個單元能干什么？今天我們來研究一下四邊形等參單元，試圖帶大家探討一下它的單元剛度矩陣。來，先跟我一起玩個小游戲，放松一下。
1個單元，4個點，8個位移，4條邊，哇
2個單元，8個點，16個位移，8條變，哇哇
3個單元，12個點，24個位移，12條變，哇哇哇……
關于四邊形等參單元剛度矩陣的推導，一般有限元教材都有詳細的介紹，這里就不再討論。書看千遍，不如公式推一遍，公式推千遍，不如代碼擼一行，我們直接上代碼，進行計算。
單元剛度矩陣如下，這是一個8*8的矩陣，我們來研究一下，該單元剛度矩陣的特性。
一、單元剛度矩陣是對稱矩陣，即滿足Kij=Kji
二、單元剛陣主對角線元素恒為正值；因為主對角元素表示力的方向和位移方向一致，故總為正值。
三、我們計算一下剛度矩陣的行列式
可以看到，行列式|K|=0，即單元剛陣是奇異陣，從物理意義上來解釋，這是因為計算單元剛陣時沒有對單元的節點加以約束，雖然，單元處于平衡狀態，但容許單元產生剛體位移，故從單元剛度平衡方程不可能得到唯一位移解。
四、我們來計算一下剛度矩陣的特征值
可以看到，矩陣共有八個特征值，其中有三個零特征值。你覺得這是偶然嗎？不，冥冥之中自有天數。剛度陣的秩為5，說明剛度矩陣只有5行是線性無關的，需要約束其中的3個自由度，方程Ku=F才能求解。而約束3個自由度，就是為了消除3個剛體位移。
展開 
有限元理論基礎及Abaqus內部實現方式研究系列1：S4殼單元剛度矩陣研究
具體的剛度矩陣在Mindlin理論基礎上的修正如下表：
項次
剛度
修正情況
說明
修正
不修正
1
薄膜效應剛度
√
Abaqus為了消除殼單元出現的剪切鎖死，計算薄膜效應時把殼作為體來處理
2
面外彎曲剛度
√
3
面外橫向剪切剛度
√
（1）Abaqus為了處理薄板情況，增加了一個幾何因子
（2）為了消除沙漏問題增加相關剛度
4
第6自由度剛度
√
Abaqus取了和對角剛度相關的小量
5
其它非對角元素
√
（1）during analysis選項導致的小量
（2）第6自由度剛度相關小量
詳細研究方法，見附件：
有限元理論基礎及Abaqus內部實現方式研究系列1：S4殼單元剛度矩陣研究(SnowWave02 20170708).pdf
展開 有限元理論基礎及Abaqus內部實現方式研究系列7：C3D8R六面體單元的剛度矩陣
介紹有限元單元正確性的驗證方法，通過多個實例比較自研結構求解器程序iSolver與Abaqus的分析結果，從而說明整個正確性驗證的過程和iSolver結果的正確性。
https://www.yqgqt.org.cn/content/post/373743
第六篇：General梁單元的剛度矩陣。介紹梁單元的基礎理論和Abaqus中General梁單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus梁單元完全一致的剛度矩陣。
https://www.yqgqt.org.cn/content/post/403932
展開 244 基于matlab的實現四節點板單元剛度矩陣求解 ￥24.9
基于matlab的實現四節點板單元剛度矩陣求解，振動模態分析?？勺灾鬏斎氩牧蠀?，板單元長寬厚尺寸。輸出振型結果，輸出多少階可自主設置。最后以可視化的形式展示，程序已調通，可直接運行。
基于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
展開 結構剛度矩陣的特點
由前面的討論可知結構的剛度矩陣K是由單元剛度矩陣集合而成，它與單元剛度矩陣類同也具有明顯的物理意義。有限元的求解方程(32)式是結構離散后每個結點的平衡方程。結構剛度矩陣K的任一元素Kij的物理意義是：結構第j個結點位移為單位值而其它結點位移皆為零時，需在第i個結點位移方向上施加的結點力的大小。與單元不同之處在于結構是單元的集合體，每個單元都對結構起一定的作用。由于單元剛度矩陣是對稱和奇異的，由它們集成的結構剛度矩陣K也是對稱和奇異的，也就是說結構至少需給出能限制剛體位移的約束條件才能消除K的奇異性，以便由(32)式求得結點位移。
連續體離散為有限個單元體，由圖1可見，每個結點的相關單元只是圍繞在該結點周圍為數甚少的幾個，一個結點通過相關單元與之發生關系的相關結點也只是它周圍的少數幾個，因此雖然總體單元數和結點數很多，結構剛度矩陣的階數很高，但剛度系數中非零系數卻很少，這就是剛度矩陣的大型和稀疏性。只要結點編號是合理的，這些稀疏的非零元素將集中在以主對角線為中心的一條帶狀區域內，即具有帶狀分布的特點。如圖7所示。
綜上所述，有限單元法最后建立的方程組的大型系數矩陣K具有以下性質:(1)對稱性（2)奇異性;(3)稀疏性；(4)非零元素呈帶狀分布。由于方程組的大型，在求解方程時，除引入位移邊界條件使奇異性消失外，其他特點都必須在解方程中予以充分的考慮和利用，以提高解題的效率。"
七、實施步驟與注意事項
利用上面討論的三角形常應變單元解平面問題，其具體步驟可歸納如下:
1)將要計算的彈性體劃分成三角形單元。對結點進行編號，列出結點坐標作為輸入信息。
(2)對單元進行編號，列出單元三個結點的號碼作為輸入信息。
(3)計算載荷的等效結點力，把等效結點力作為輸入信息。
展開 ABAQUS 剛度矩陣 ￥1000
是這樣的：
1，TestElement.py 是編寫的測試代碼，可輸出8結點線性單元的剛度矩陣。
2，abaqus文件Job-testing.inp ，運行可以輸出單元剛度矩陣。
問題是：他們的結點坐標/排序，材料參數都一致， 但得到的剛度矩陣就不一樣。
需要：代碼輸出的剛度矩陣與abaqus得到的剛度矩陣一致或者基本一致。
?abaqus提取單元的剛度矩陣和質量矩陣
abaqus提取單元的剛度矩陣和質量矩陣：
*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
有限元理論基礎及Abaqus內部實現方式研究系列22： 幾何非線性的剛度矩陣求解
http://www.yqgqt.org.cn/content/post/360565
第五篇：單元正確性驗證。介紹有限元單元正確性的驗證方法，通過多個實例比較自研結構求解器程序iSolver與Abaqus的分析結果，從而說明整個正確性驗證的過程和iSolver結果的正確性。
https://www.yqgqt.org.cn/content/post/373743
第六篇：General梁單元的剛度矩陣。介紹梁單元的基礎理論和Abaqus中General梁單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus梁單元完全一致的剛度矩陣。
https://www.yqgqt.org.cn/content/post/403932
第七篇：C3D8六面體單元的剛度矩陣。介紹六面體單元的基礎理論和Abaqus中C3D8R六面體單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus六面體單元完全一致的剛度矩陣。
https://www.yqgqt.org.cn/content/post/430177
第八篇：UMAT用戶子程序開發步驟。介紹基于Fortran和Matlab兩種方式的Abaqus的UMAT的開發步驟，對比發現開發步驟基本相同，同時采用Matlab更加高效和靈活。
https://www.yqgqt.org.cn/content/post/432848
第九篇：編寫線性UMAT Step By Step。介紹基于Matlab線性零基礎，從零開始Step by Step的UMAT的編寫和調試方法，幫助初學者UMAT入門。
http://www.yqgqt.org.cn/content/post/440874
第十篇：耦合約束（Coupling constraints）的研究。
展開 平面四邊形單元的剛度矩陣的推導 ￥2
本圖片是取自《有限元方法基礎教程》（國際單位制）（第五版）， 第六章 建立平面應力和平面應變剛度方程。
==》本博客是對于四節點四邊形單元的剛度矩陣的推導，沒有編程的實現。
ANSYS中整體、單元剛度和質量矩陣的提取
指定輸出單元矩陣
/SOLU
SOLVE
finish
/OUTPUT, TERM ! 將輸出信息送到output windows中
! 這時用編輯器打開cp.out文件，可以看到按單元寫出的質量、剛度等矩
陣
ANSYS中整體、單元剛度和質量矩陣的提取.rar
ANSYS模型剛度、質量矩陣快速提取小軟件—km_from_Ansys ￥88
提取方法
Ansys中其實提供了幾種單元剛度矩陣和整體剛度矩陣的提取接口，先做一下簡單介紹：
（1）單元剛度矩陣提取，集成整體剛度矩陣
首先選擇要提取的單元，通過/debug命令輸出選擇集中的單元剛度矩陣，再根據結構的單元連接信息，組合成整體剛度矩陣，命令如下：
/OUTPUT,cp,out,, ! 將輸出信息送到cp.out文件
/debug,-1,,,1 ! 指定輸出單元矩陣
/SOLU
SOLVE
finish
/OUTPUT, TERM ! 將輸出信息送到output windows中
感受：這個方法很直接，但稍顯麻煩，需要提取到結構的單元組成信息、節點信息以及單元剛度矩陣等，對于稍大一些結構來說，可能需要1個小時甚至更長時間才能完成。
（2）超單元方法
/solu
antype,7 !substructuring分析類型
seopt,matname,1 !設置文件名稱和剛度矩陣類型(剛度，質量，阻尼等)
nsel,all !選擇所有節點
m,all,all !定義所有節點自由度為主自由度
solve !求解
selist,matname,3 !列出整體剛度矩陣
這種方法是可以提取到結構整體剛度矩陣的，但是問題在于需要指定主自由度，對于一般結構而言，這個方法還是適用的，問題是如果結構中存在耦合關系，結構自由度存在依賴關系，該方法就不太適用了，可能還有其他解決辦法，但是我后來還是放棄這種方法。
（3）HBMAT命令方法提取整體矩陣
這種方法的適用性比較廣，也是個人比較推崇的一種方法。
展開 