剛度矩陣的推導的案例
平面四邊形單元的剛度矩陣的推導 ￥2
本圖片是取自《有限元方法基礎教程》（國際單位制）（第五版）， 第六章 建立平面應力和平面應變剛度方程。
==》本博客是對于四節點四邊形單元的剛度矩陣的推導，沒有編程的實現。
J-OCTA利用聲子色散曲線推導剛度矩陣
材料的彈性特征可以使用剛度矩陣來進行完全描述。剛度矩陣中任意元素都可以根據由第一性原理計算得出的聲子色散曲線來確定。
在本次案例研究中，我們介紹了一個用立方體結構來分析硅、金剛石和銅單晶的剛度矩陣的方法。
如圖1，硅是具有金剛石結構的立方晶體，其剛度矩陣如下所示。
由于立方晶體具有對稱性，使用SIESTA模擬軟件計算可得該模型中有3個獨立元素：C11C11，C12C12，和C44C44。
J-OCTA利用聲子色散曲線推導剛度矩陣
材料的彈性特征可以使用剛度矩陣來進行完全描述。剛度矩陣中任意元素都可以根據由第一性原理計算得出的聲子色散曲線來確定。
在本次案例研究中，我們介紹了一個用立方體結構來分析硅、金剛石和銅單晶的剛度矩陣的方法。
如圖1，硅是具有金剛石結構的立方晶體，其剛度矩陣如下所示。
由于立方晶體具有對稱性，使用SIESTA模擬軟件計算可得該模型中有3個獨立元素：C11C11，C12C12，和C44C44。
有限元分析的底層邏輯是什么？
由于三角形單元涉及到二維應變，通常通過應變-位移矩陣來計算。 - 剛度矩陣：經過推導后，三角形單元的剛度矩陣為：??=∫??????????,????
其中，??
是應變-位移矩陣，??
是材料的剛度矩陣，??
是單元的面積。
三角形單元剛度矩陣（簡化版）
對于簡單的線性三角形單元，剛度矩陣一般可以通過以下積分公式推導：
??=??4(120?12012?12?12?1212)
其中??
是三角形單元的面積。
2.3 3D 單元（如四面體單元和六面體單元）
在三維問題中，常用的單元有四面體單元和六面體單元。
四面體單元 四面體單元的推導較為復雜，通常需要通過形函數和積分來獲得剛度矩陣。這里給出大致的推導思路。 步驟：
- 位移場：四面體單元的位移場是通過節點的位移進行插值，通常是多項式形式的。
- 應變能：通過應變-位移關系，計算應變能。應變-位移矩陣 ??
是由形函數的導數組成的。
- 剛度矩陣：同樣利用應變-位移矩陣和材料剛度矩陣 ??
，通過積分得到四面體單元的剛度矩陣。
四面體單元的剛度矩陣一般較為復雜，具體計算通常依賴數值積分（例如高斯積分）。
每種單元的剛度矩陣的推導方法都是基于能量原理（如虛功原理、最小勢能原理）或通過變分法進行的。1D單元的剛度矩陣推導較為簡單，2D和3D單元則需要根據單元的具體幾何形狀和物理特性（如材料性質、形函數等）進行推導。
3. 單元剛度矩陣的組裝
在構造出各單元的剛度矩陣之后，需要將它們根據結構中單元之間的連接關系組裝成整體的剛度矩陣。組裝過程通常依據節點的連接情況，將每個單元的剛度矩陣“嵌入”到整體的剛度矩陣中。
4.
展開 
利用聲子色散曲線推導剛度矩陣
根據取得的聲速，可以按以下公式計算剛度矩陣。
聲速與(1,0,0)/(1,1,1)方向上剛度矩陣的關系。
各個箭頭表示偏振方向，紅色箭頭表示縱波，綠色箭頭表示橫波。ρ表示硅的重量密度。
表1給出分析所得剛度矩陣中每個元素值和文獻中的值之對比。（C. Kittel 《固體物理學導論》，紅色括號）。
表1 剛度矩陣所得值與文獻中的值對比
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)方向上剛度矩陣的關系。
展開 有限元理論基礎及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 理想的求解方式
理論上受力曲線是一條光滑曲線，計算機沒法求解曲線上每個時刻點的結果，只能求解部分有限間隔點的結果。非線性問題不是一條直線，所以需要多次迭代才能實現。
展開 鋼筋混凝土結構有限元分析單元類型和分析模型 附混凝土結構有限元分析下載
單元剛度矩陣的推導與一般有限元相同。
1.2 組合式模型
組合式模型是假設鋼筋以一個確定的角度分布在整個單元中，并假設混凝土與鋼筋之間存在著良好的粘結，認為兩者之間無滑移。又分為分層組合方式和帶鋼筋膜的方式等。該單元剛度矩陣推導時分別求出各自的單元剛度，然后組合起來。
1.3 整體式模型
整體式模型是假設鋼筋分布于整個單元中，并把單元視為連續均勻材料（如ANSYS中的四面體等實體單元-solid65單元選擇混凝土材料時），采用混凝土-鋼筋復合的本構關系，把混凝土、鋼筋二者的貢獻組合起來，一次求得綜合的單元剛度矩陣。
后兩種模型共同點是它們的單元剛度矩陣都是反映鋼筋混凝土的綜合剛度。
下載地址：混凝土結構有限元分析
展開 鋼筋混凝土結構有限元分析單元類型和分析模型
單元剛度矩陣的推導與一般有限元相同。
1.2 組合式模型
組合式模型是假設鋼筋以一個確定的角度分布在整個單元中，并假設混凝土與鋼筋之間存在著良好的粘結，認為兩者之間無滑移。又分為分層組合方式和帶鋼筋膜的方式等。該單元剛度矩陣推導時分別求出各自的單元剛度，然后組合起來。
1.3 整體式模型
整體式模型是假設鋼筋分布于整個單元中，并把單元視為連續均勻材料（如ANSYS中的四面體等實體單元-solid65單元選擇混凝土材料時），采用混凝土-鋼筋復合的本構關系，把混凝土、鋼筋二者的貢獻組合起來，一次求得綜合的單元剛度矩陣。
后兩種模型共同點是它們的單元剛度矩陣都是反映鋼筋混凝土的綜合剛度。
展開 提取整體剛度矩陣、質量矩陣及阻尼矩陣的三種方法
指定輸出單元矩陣
/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 !列出整體剛度矩陣
展開 如何從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中的矩陣變量形式，我們可以直接運用這段源代碼來進行操作。
展開 
基于 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 中如何輸出質量矩陣和剛度矩陣
我記得在patran 中輸出的bdf文件中輸入一條命令，運行一下就可以了。但是記不得了。請教高手，先謝謝了。
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及以后版本。
展開 本構模型_Umat子程序中Jacobian矩陣的推導示例
本構模型_Umat子程序中Jacobian矩陣的推導示例
