Matlab有限元計(jì)算

2017-08-25 by:CAE仿真在線(xiàn) 來(lái)源:互聯(lián)網(wǎng)

對(duì)于開(kāi)始學(xué)習(xí)有限元的人來(lái)說(shuō),簡(jiǎn)單了解下有限元計(jì)算的基本思路是很有必要的,因此這篇文章僅僅以一個(gè)簡(jiǎn)單的一維線(xiàn)性彈簧單元為例,簡(jiǎn)述一下有限元計(jì)算的基本思想,主要參考書(shū)目《有限元方法基礎(chǔ)教程》-Daryl L. Logan和《matlab有限元分析與應(yīng)用》-P.I. Kattan。前者有詳細(xì)的剛度矩陣等推導(dǎo),主要從理論上對(duì)有限元方法進(jìn)行說(shuō)明,后者主要是結(jié)合著名的矩陣實(shí)驗(yàn)室matlab對(duì)結(jié)構(gòu)有限元進(jìn)行分析與計(jì)算,兩者可以結(jié)合來(lái)看。

下面我就《有限元方法基礎(chǔ)教程》一書(shū)中的習(xí)題2.17為例,使用matlab來(lái)說(shuō)明一下如何對(duì)其進(jìn)行計(jì)算。

習(xí)題2.17

如圖所示的五個(gè)彈簧組裝,確定節(jié)點(diǎn)2和節(jié)點(diǎn)3的位移、節(jié)點(diǎn)1和節(jié)點(diǎn)4的反力。假設(shè)節(jié)點(diǎn)2和節(jié)點(diǎn)3處連接彈簧的剛性垂直桿一直保持水平,但是自由滑動(dòng)、向左或者向右移動(dòng)。在節(jié)點(diǎn)向右施加外力1000N。令k1=500N/mm,k2=k3=300N/mm,k4=k5=400N/mm。

眾所周知,有限元計(jì)算的基本流程主要包含六個(gè)部分:

①離散化域

②寫(xiě)出單元?jiǎng)偠染仃?

③集成整體剛度矩陣

④引入邊界條件

⑤解方程

⑥后處理

①離散化域

有限元模型的建立有兩種方法,一種是直接建立,另一種是網(wǎng)格劃分,對(duì)于像桿單元,彈簧單元質(zhì)量單元等直接建立有限元模型比較方便,但是像梁?jiǎn)卧?實(shí)體單元等人為建立就很復(fù)雜,一般就是先建立好幾何模型再進(jìn)行網(wǎng)格劃分,以此得到有限元模型。

對(duì)于此題,使用直接建立法即可。由圖中可以看出可以離散成五個(gè)彈簧單元,各彈簧單元以及對(duì)應(yīng)節(jié)點(diǎn)如下:

②寫(xiě)出單元?jiǎng)偠染仃?

對(duì)于靜力結(jié)構(gòu)分析來(lái)說(shuō),最大的難點(diǎn)就在單元?jiǎng)偠染仃嚨牡贸觥：?jiǎn)單的單元?jiǎng)偠染仃囅癖绢}中的彈簧單元使用直接平衡法就能得到,復(fù)雜的可能就要用到能量法或者加權(quán)殘余法等。

本例中直接給出彈簧單元的剛度矩陣形式,這樣,根據(jù)該矩陣形式,我們可以構(gòu)造一個(gè)函數(shù)SpringElementStiffness專(zhuān)門(mén)用于得到彈簧單元的剛度矩陣,下面是《matlab有限元分析與應(yīng)用》配套給出的彈簧單元?jiǎng)偠染仃嚿珊瘮?shù)(以下函數(shù)均為書(shū)中給出):

這樣,通過(guò)調(diào)用該函數(shù)能夠方便的生成上述5個(gè)單元的剛度矩陣:

③集成整體剛度矩陣(直接剛度法)

得到單元?jiǎng)偠染仃囍?可以直接利用疊加法得到整體剛度矩陣。由于每個(gè)節(jié)點(diǎn)的力是由分擔(dān)該節(jié)點(diǎn)的單元節(jié)點(diǎn)相加而成,因此能夠直接通過(guò)將各單元對(duì)應(yīng)節(jié)點(diǎn)疊加,得到總體剛度矩陣,這個(gè)也叫直接剛度法,對(duì)應(yīng)的matlab程序如下:

將單元矩陣k的對(duì)應(yīng)元素疊加到總體剛度矩陣對(duì)應(yīng)的節(jié)點(diǎn)位置,由于這里總共有五個(gè)單元,因此需要調(diào)用該函數(shù)五次得到總體剛度矩陣:

④引入邊界條件

對(duì)于上圖中,邊界條件屬于混合型,既有位移型也有力型,具體有:

u1=0 u4=0 f2=0 f3=1000N

其中

u3 u4 f1 f4是未知的,需要我們求出。引入上述邊界條件得到以下總體方程:

⑤解方程

上述方程組可以直接提取第三行和第二行先進(jìn)行求解得到u2和u3的值,即求下面方程:

直接使用矩陣反除求出u2與u3的值,再將位移值帶入總體位移矩陣,得到力值。由于后面會(huì)多次使用力值計(jì)算程序,因此將該程序也編寫(xiě)為m文件進(jìn)行調(diào)用,如下:

然后計(jì)算得到總位移與總力矩陣

⑥后處理

所有求解完畢之后,可以返回來(lái)從新得到各個(gè)單元內(nèi)力,具體計(jì)算程序如下:

通過(guò)上述求解結(jié)果可以得到各節(jié)點(diǎn)的位移以及單元受力情況。當(dāng)然這只是對(duì)于最簡(jiǎn)單的一維線(xiàn)性彈簧單元,如果要擴(kuò)展到二維,三維,單元更換為桿,梁,面等,單元?jiǎng)偠染仃囷@然會(huì)更加復(fù)雜,對(duì)于具體的邊界可能還要對(duì)剛度矩陣進(jìn)行修正,但是整體的求解思路類(lèi)似。

開(kāi)放分享：優(yōu)質(zhì)有限元技術(shù)文章,助你自學(xué)成才