剛度矩陣的案例
如何從Ansys APDL中提取剛度矩陣與質量矩陣？ ￥69
其可以直接導出full文件中的矩陣數據，而通過我下面給大家提供的源代碼便可以直接將導出的質量矩陣與剛度矩陣直接轉化為matlab中可使用變量，實現我們的矩陣提取操作。
Hbmat法：
! 提取剛度矩陣
/AUX2
FILE,'file',full ! 將’file’改為自己路徑下的.full文件名
HBMAT, 'Stiffness_mat', dat, , ASCII, STIFF, YES, YES !剛度矩陣
HBMAT, 'Mass_mat', dat, ,ASCII, MASS, YES, YES ! 質量矩陣
FINISH
成功導出后你會在你的工作路徑中看到儲存在Ansys中的剛度矩陣與質量矩陣。
2.2GUI方法
目前大部分此類教程都忽略了最基本的GUI方法，可能是自動帶入了經驗豐富的工程師角色，但對于初學者而言，GUI方法十分有助于理解與學習有限元軟件，因此我也將介紹提取剛度質量矩陣的GUI方法。
從而我們便能導出所需要的剛度矩陣與質量矩陣源文件了。
3.源代碼的使用與轉換（使用方法）
我們打開所導出剛度矩陣會發現里面的數字比較混亂，一時間分不出哪些是我們需要的數據，也沒有辦法直接拿它用來計算，這時，我們便需要使用我們自己的矩陣轉換m文件，其可以自動幫我們提取出剛度矩陣與質量矩陣，并生成matlab中的矩陣變量形式，我們可以直接運用這段源代碼來進行操作。
展開 MSC Nastran模型剛度矩陣和質量矩陣的輸出方法
1 概述
MSC Nastran模型的剛度矩陣和質量矩陣，可以輸出為文本文件。工程實際中，工程師可以校核、集成矩陣，進行第二次開發，完成商用軟件和自研程序的完美集成。例如：工程師有一個計算線性動力學方程組的瞬態python程序，可以集成MSC Nastran的剛度矩陣和質量矩陣。
2 剛度矩陣和質量矩陣的輸出方法
1) 剛度矩陣和質量矩陣輸出至punch（.pch）文件
如果需要在其他MSC Nastran計算中，重用MSC Nastran模型的矩陣，可以將MSC Nastran矩陣輸出至Punch文件，方法為：
l 在MSC Nastran卡片中，添加參數：PARAM,EXTOUT,DMIGPCH
注：Punch文件中的矩陣，Patran不支持
2) 剛度矩陣和質量矩陣輸出至f06（.f06）文件
如果想直接在f06中查看輸出的矩陣，可以使用如下方法：
l 在執行控制部分（CEND前），添加如下卡片：
COMPILE EXTOUT $
ALTER 'RETURN'(,-1) $
MATPRN KAA,,,,// $
MATPRN MAA,,,,// $
l 添加如下參數（BEGIN BULK），例如：PARAM,EXTOUT,DMIGPCH
實例：
輸入文件：
剛度矩陣
質量矩陣
3 參考信息
適用版本：MSC Nastran 2005及以后版本。
展開 基于 MATLAB 的 ANSYS Harwell-Boeing 格式稀疏矩陣提取工具 —— 剛度矩陣與質量矩陣 ￥30
在有限元分析中，ANSYS 可以導出大規模稀疏矩陣（如剛度矩陣、質量矩陣），通常使用 Harwell-Boeing (HB) CCS 格式。這些矩陣對后續二次開發、動力學分析或自定義求解器非常重要，但由于其稀疏和壓縮存儲形式，直接在 MATLAB 中讀取和使用并不方便。
本文提供了 兩個 MATLAB 函數，可直接從 ANSYS 導出的 HB 矩陣文件中讀取并重構成 MATLAB 稀疏矩陣：
1.剛度矩陣提取函數
輸入：ANSYS 導出的剛度矩陣 HB 文件（stiff.txt）
輸出：MATLAB 稀疏矩陣 K，可直接用于動力學計算或驗證
支持自動對稱化，保證數值正確
2.質量矩陣提取函數
輸入：ANSYS 導出的質量矩陣 HB 文件（mass.txt）
輸出：MATLAB 稀疏矩陣 M
使用與剛度矩陣同樣的解析邏輯，無需額外修改
案例說明：
本文以高速鐵路接觸網結構為例，展示了如何將 ANSYS 中導出的稀疏剛度矩陣和質量矩陣，在 MATLAB 中完整展開，并進行后續動力學分析準備。
通過該方法，可以將大規模有限元矩陣快速轉化為 MATLAB 可操作形式，為自定義振動分析、模態分析及其他科研或工程應用提供基礎。
優勢與應用：
支持大規模稀疏矩陣解析
自動對稱化，保證數值精度
適用于剛度矩陣、質量矩陣、其他 HB 格式矩陣
可作為動力學求解器或后處理工具的基礎模塊
使用方法：
1.使用以下代碼對ansys中生成的質量及剛度矩陣進行提取，file,5,full（5為工作目錄下full文件的文件名，例如：filename.full）。
展開 NASTRAN輸出剛度矩陣的方法
知道含義以后即可借助MATLAB或者其它軟件，讀取pch中的剛度矩陣，并編寫代碼完成剛度矩陣的輸出。圖2是小翼做的一個結構的總剛度矩陣的局部，經過與自己編寫的結構剛度矩陣輸出代碼計算的剛度矩陣對比，發現一致性較好，部分剛度系數誤差在3~4%左右，當然還是以NASTRAN的精度為準。
圖2 某結構總剛度矩陣局部
單元剛度矩陣
由（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個結點位移方向上施加的結點力的大小。單元剛性大，則使結點產生單位位移所需施加的結點力就大。因此單元剛度矩陣中的每個元素反映了單元剛性的大小，稱為剛度系數。
展開 提取ANSYS中的剛度矩陣和質量矩陣到MATLAB中進行二次開發 ￥88
對稱矩陣只存儲下三角元素，如結構剛度矩陣為對稱矩陣，Harwell-Boeing格式則僅記錄下三角元素。
通過上述解釋我們可以知道該剛度矩陣有1027個列指針，有7978個行索引，有1026個右邊項，矩陣的維度為1026×1026。
要將上述文件讀取到MATLAB中并轉化為常規的稀疏矩陣，需要了解Harwell-Boeing壓縮儲存格式及重構方法，這里面采用的是CSC方式，以列指針、行索引的方式去重構稀疏矩陣。
x.y MSC Nastran模型剛度矩陣和質量矩陣的輸出方法
1 概述
MSC Nastran模型的剛度矩陣和質量矩陣，可以輸出為文本文件。工程實際中，工程師可以校核、集成矩陣，進行第二次開發，完成商用軟件和自研程序的完美集成。例如：工程師有一個計算線性動力學方程組的瞬態python程序，可以集成MSC Nastran的剛度矩陣和質量矩陣。
2 剛度矩陣和質量矩陣的輸出方法
1) 剛度矩陣和質量矩陣輸出至punch（.pch）文件
如果需要在其他MSC Nastran計算中，重用MSC Nastran模型的矩陣，可以將MSC Nastran矩陣輸出至Punch文件，方法為：
l 在MSC Nastran卡片中，添加參數：PARAM,EXTOUT,DMIGPCH
注：Punch文件中的矩陣，Patran不支持
2) 剛度矩陣和質量矩陣輸出至f06（.f06）文件
如果想直接在f06中查看輸出的矩陣，可以使用如下方法：
l 在執行控制部分（CEND前），添加如下卡片：
COMPILE EXTOUT $
ALTER 'RETURN'(,-1) $
MATPRN KAA,,,,// $
MATPRN MAA,,,,// $
l 添加如下參數（BEGIN BULK），例如：PARAM,EXTOUT,DMIGPCH
實例：
輸入文件：
剛度矩陣
質量矩陣
3 參考信息
適用版本：MSC Nastran 2005及以后版本。
展開 剛度矩陣非正定
剛度矩陣非正定
結構：鋼桁架簡支梁橋（一端有加勁弦）
有限元模擬：采用桁段有限元法結合對號入座法則，推導出了結構的剛度矩陣
邊界條件：將邊界條件所限制的自由度對應剛度矩陣中的數值修改為一個很大的數
問題：剛度矩陣非正定，出現了復模態。問題可能出現在兩個地方，一個是邊界條件，另一個是矩陣的原素，我都仔細檢查了一邊，還是不行。所以發帖子上來請教論壇里的高手們。感謝。
附件為剛度矩陣圖形和前十階振型圖
剛度矩陣.jpg
振型.jpg
部分源代碼.rar
SIESTA模擬利用聲子色散曲線推導剛度矩陣
材料的彈性特征可以使用剛度矩陣來進行完全描述。剛度矩陣中任意元素都可以根據由第一性原理計算得出的聲子色散曲線來確定。
在本次案例研究中，我們介紹了一個用立方體結構來分析硅、金剛石和銅單晶的剛度矩陣的方法。
如圖1，硅是具有金剛石結構的立方晶體，其剛度矩陣如下所示。
由于立方晶體具有對稱性，使用SIESTA模擬軟件計算可得該模型中有3個獨立元素：C11C11，C12C12，和C44C44。
圖1 硅的金剛石結構
立方晶體的剛度矩陣
硅的聲子色散曲線可以用SEISTA模擬軟件直接計算（如圖2）。
根據硅的晶胞中原子數為2這一事實依據，硅有三種聲振模式和三種光學模式。硅是具有金剛石結構的立方晶體，其色散曲線表現出不同晶體取向的色散特性。
這里我們重點關注從ΓX點指向的(1,0,0)方向，和從ΓL點指向的(1,1,1)方向擴展的聲振模式色散特性。聲振模式的色散特性在近場的長波區呈線性Γ，且該梯度曲線給出了聲速。
圖2 硅的聲子色散曲線
圖3 區域1和2的色散曲線
圖中虛線表示靠近該區域的每個傳播方向上的聲速Γ點，下標L/T代表縱波和橫波。
在（1,0,0）方向傳播的聲波由一個縱波和兩個簡并后的橫波組成，它們各自的聲速和剛度矩陣可以表示為下列關系式，由此結果可得C11C11 和 C44C44。
進一步使用該結果，根據（1,1,1）方向傳播的聲波可以得到c12。在本案例中，梯度（聲速）在Γ點可以通過使用圖中1和2區域的三階最小二乘法構造一條近似曲線獲得，這是SIESTA模塊的聲子分析功能。根據取得的聲速，可以按以下公式計算剛度矩陣。
聲速與(1,0,0)/(1,1,1)方向上剛度矩陣的關系。
展開 提取整體剛度矩陣、質量矩陣及阻尼矩陣的三種方法
指定輸出單元矩陣
/SOLU
SOLVE
finish
/OUTPUT, TERM ! 將輸出信息送到output windows中
! 這時用編輯器打開cp.out文件，可以看到按單元寫出的質量、剛度等矩陣
3.
其原理很簡單，即使用ansys的超單元即可解決問題。定義超單元，然后列出超單元的剛度矩陣即可。
面是一個小例題，自可明白。
/prep7
k,1
k,2,3000
l,1,2
et,1,beam3
mp,ex,1,2e5
mp,prxy,1,0.3
r,1,5000,2e7,200
lesize,all,,,10
lmesh,all
finish
!----以上正常建立模型，不必施加約束和荷載
/solu
antype,7 !substructuring分析類型
seopt,matname,1 !設置文件名稱和剛度矩陣類型(剛度，質量，阻尼等)
nsel,all !選擇所有節點
m,all,all !定義所有節點自由度為主自由度
solve !求解
selist,matname,3 !列出整體剛度矩陣
展開 利用聲子色散曲線推導剛度矩陣
根據取得的聲速，可以按以下公式計算剛度矩陣。
聲速與(1,0,0)/(1,1,1)方向上剛度矩陣的關系。
各個箭頭表示偏振方向，紅色箭頭表示縱波，綠色箭頭表示橫波。ρ表示硅的重量密度。
表1給出分析所得剛度矩陣中每個元素值和文獻中的值之對比。（C. Kittel 《固體物理學導論》，紅色括號）。
表1 剛度矩陣所得值與文獻中的值對比
ABAQUS 剛度矩陣 ￥1000
是這樣的：
1，TestElement.py 是編寫的測試代碼，可輸出8結點線性單元的剛度矩陣。
2，abaqus文件Job-testing.inp ，運行可以輸出單元剛度矩陣。
問題是：他們的結點坐標/排序，材料參數都一致， 但得到的剛度矩陣就不一樣。
需要：代碼輸出的剛度矩陣與abaqus得到的剛度矩陣一致或者基本一致。
有限元理論基礎及Abaqus內部實現方式研究系列22： 幾何非線性的剛度矩陣求解
自主結構有限元求解器iSolver介紹視頻：
http://www.yqgqt.org.cn/college/video/c12884
==第22篇：幾何非線性的剛度矩陣求解==
幾何非線性在界面上是很容易設置的，但商軟內部的處理相當復雜，我們從最基本的剛度矩陣的求解出發，看看在幾何非線性設置后，剛度矩陣具體是怎么實現的。本文首先介紹幾何非線性下的剛度矩陣的理論推導和計算機求解方法，說明理想的求解方式的困難點和猜測Abaqus內部的解決方法。最后利用一個簡單的算例通過對比iSolver和Abaqus的結果，部分驗證我們對Abaqus幾何非線性的剛度矩陣的實現方式的猜測。
1.1 幾何非線性的剛度矩陣推導理論
在前面17章：幾何非線性的物理含義中，我們提到如果是非線性系統，應變能W隨t的變化就是個非線性過程。每個時刻點可以求出一個斜率，這個斜率最終會形成當前時刻點的剛度矩陣。
求導后得到的剛度K：
也就是剛度矩陣將分為兩塊：
（1） 上式的前面一部分稱為材料剛度陣，依然是以前的BDB形式，只不過B換成了當前時刻的應變位移矩陣
（2） 后面新增項一般稱為幾何剛度陣，在Abaqus中稱為初始應力矩陣（initial stress stiffness）。
1.2 幾何非線性的剛度矩陣計算機求解
1.2.1 理想的求解方式
理論上受力曲線是一條光滑曲線，計算機沒法求解曲線上每個時刻點的結果，只能求解部分有限間隔點的結果。非線性問題不是一條直線，所以需要多次迭代才能實現。
展開 結構剛度矩陣的特點
由前面的討論可知結構的剛度矩陣K是由單元剛度矩陣集合而成，它與單元剛度矩陣類同也具有明顯的物理意義。有限元的求解方程(32)式是結構離散后每個結點的平衡方程。結構剛度矩陣K的任一元素Kij的物理意義是：結構第j個結點位移為單位值而其它結點位移皆為零時，需在第i個結點位移方向上施加的結點力的大小。與單元不同之處在于結構是單元的集合體，每個單元都對結構起一定的作用。由于單元剛度矩陣是對稱和奇異的，由它們集成的結構剛度矩陣K也是對稱和奇異的，也就是說結構至少需給出能限制剛體位移的約束條件才能消除K的奇異性，以便由(32)式求得結點位移。
連續體離散為有限個單元體，由圖1可見，每個結點的相關單元只是圍繞在該結點周圍為數甚少的幾個，一個結點通過相關單元與之發生關系的相關結點也只是它周圍的少數幾個，因此雖然總體單元數和結點數很多，結構剛度矩陣的階數很高，但剛度系數中非零系數卻很少，這就是剛度矩陣的大型和稀疏性。只要結點編號是合理的，這些稀疏的非零元素將集中在以主對角線為中心的一條帶狀區域內，即具有帶狀分布的特點。如圖7所示。
綜上所述，有限單元法最后建立的方程組的大型系數矩陣K具有以下性質:(1)對稱性（2)奇異性;(3)稀疏性；(4)非零元素呈帶狀分布。由于方程組的大型，在求解方程時，除引入位移邊界條件使奇異性消失外，其他特點都必須在解方程中予以充分的考慮和利用，以提高解題的效率。"
七、實施步驟與注意事項
利用上面討論的三角形常應變單元解平面問題，其具體步驟可歸納如下:
1)將要計算的彈性體劃分成三角形單元。對結點進行編號，列出結點坐標作為輸入信息。
(2)對單元進行編號，列出單元三個結點的號碼作為輸入信息。
(3)計算載荷的等效結點力，把等效結點力作為輸入信息。
展開 有限元理論基礎及Abaqus內部實現方式研究系列12： 幾何梁單元的剛度矩陣 ￥1
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。介紹了線性UMAT的接口功能和關鍵接口變量的含義，并通過簡單立方體靜力分析的算例詳細說明了基于Matlab線性UMAT的開發步驟。
http://www.yqgqt.org.cn/content/post/440874
第十篇：耦合約束（Coupling constraints）的研究。介紹了耦合約束的定義和用途，具體闡述了Abaqus中運動耦合約束和分布耦合約束的原理。
http://www.yqgqt.org.cn/content/post/531029
第十一篇：自主CAE開發實戰經驗第一階段總結。
展開 