搜索Matlab有限元計算
2017-08-25 by:CAE仿真在線 來源:互聯網
對于開始學習有限元的人來說,簡單了解下有限元計算的基本思路是很有必要的,因此這篇文章僅僅以一個簡單的一維線性彈簧單元為例,簡述一下有限元計算的基本思想,主要參考書目《有限元方法基礎教程》-Daryl L. Logan和《matlab有限元分析與應用》-P.I. Kattan。前者有詳細的剛度矩陣等推導,主要從理論上對有限元方法進行說明,后者主要是結合著名的矩陣實驗室matlab對結構有限元進行分析與計算,兩者可以結合來看。
下面我就《有限元方法基礎教程》一書中的習題2.17為例,使用matlab來說明一下如何對其進行計算。
習題2.17
如圖所示的五個彈簧組裝,確定節點2和節點3的位移、節點1和節點4的反力。假設節點2和節點3處連接彈簧的剛性垂直桿一直保持水平,但是自由滑動、向左或者向右移動。在節點向右施加外力1000N。令k1=500N/mm,k2=k3=300N/mm,k4=k5=400N/mm。
眾所周知,有限元計算的基本流程主要包含六個部分:
①離散化域
②寫出單元剛度矩陣
③集成整體剛度矩陣
④引入邊界條件
⑤解方程
⑥后處理
①離散化域
有限元模型的建立有兩種方法,一種是直接建立,另一種是網格劃分,對于像桿單元,彈簧單元質量單元等直接建立有限元模型比較方便,但是像梁單元,實體單元等人為建立就很復雜,一般就是先建立好幾何模型再進行網格劃分,以此得到有限元模型。
對于此題,使用直接建立法即可。由圖中可以看出可以離散成五個彈簧單元,各彈簧單元以及對應節點如下:
②寫出單元剛度矩陣
對于靜力結構分析來說,最大的難點就在單元剛度矩陣的得出。簡單的單元剛度矩陣像本題中的彈簧單元使用直接平衡法就能得到,復雜的可能就要用到能量法或者加權殘余法等。
本例中直接給出彈簧單元的剛度矩陣形式,這樣,根據該矩陣形式,我們可以構造一個函數SpringElementStiffness專門用于得到彈簧單元的剛度矩陣,下面是《matlab有限元分析與應用》配套給出的彈簧單元剛度矩陣生成函數(以下函數均為書中給出):
這樣,通過調用該函數能夠方便的生成上述5個單元的剛度矩陣:
③集成整體剛度矩陣(直接剛度法)
得到單元剛度矩陣之后,可以直接利用疊加法得到整體剛度矩陣。由于每個節點的力是由分擔該節點的單元節點相加而成,因此能夠直接通過將各單元對應節點疊加,得到總體剛度矩陣,這個也叫直接剛度法,對應的matlab程序如下:
將單元矩陣k的對應元素疊加到總體剛度矩陣對應的節點位置,由于這里總共有五個單元,因此需要調用該函數五次得到總體剛度矩陣:
④引入邊界條件
對于上圖中,邊界條件屬于混合型,既有位移型也有力型,具體有:
u1=0 u4=0 f2=0 f3=1000N
其中
u3 u4 f1 f4是未知的,需要我們求出。引入上述邊界條件得到以下總體方程:
⑤解方程
上述方程組可以直接提取第三行和第二行先進行求解得到u2和u3的值,即求下面方程:
直接使用矩陣反除求出u2與u3的值,再將位移值帶入總體位移矩陣,得到力值。由于后面會多次使用力值計算程序,因此將該程序也編寫為m文件進行調用,如下:
然后計算得到總位移與總力矩陣
⑥后處理
所有求解完畢之后,可以返回來從新得到各個單元內力,具體計算程序如下:
通過上述求解結果可以得到各節點的位移以及單元受力情況。當然這只是對于最簡單的一維線性彈簧單元,如果要擴展到二維,三維,單元更換為桿,梁,面等,單元剛度矩陣顯然會更加復雜,對于具體的邊界可能還要對剛度矩陣進行修正,但是整體的求解思路類似。
相關標簽搜索：Matlab有限元計算 MatLab培訓 MatLab培訓課程 MatLab在線視頻教程 MatLab技術學習教程 MatLab軟件教程 MatLab資料下載 MatLab代做 MatLab基礎知識 Fluent、CFX流體分析 HFSS電磁分析 Ansys培訓 Abaqus培訓
