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