導(dǎo)讀:自9月以來(lái),我的《教你Matlab有限元編程對(duì)懸臂梁進(jìn)行受力分析》原創(chuàng)文章發(fā)布以來(lái),我的精品課《Matlab有限元編程從入門到精通10講》受到了越來(lái)越多的工程師朋友的關(guān)注,已超過(guò)50多人已經(jīng)訂閱學(xué)習(xí),歡迎大家加入訂閱用戶交流群一起討論Matlab有限元編程那些事(哈哈哈,為大家答疑和分享資料)。
今天,我接著分享Matlab有限元編程專業(yè)技能。之前的課程我們學(xué)習(xí)了一維梁?jiǎn)卧?,二維平面單元,三維板殼單元的matlab有限元編程,本次案例主要講解如何用matlab實(shí)現(xiàn)針對(duì)四面體單元?jiǎng)澐值娜S結(jié)構(gòu)進(jìn)行有限元編程,具體案例是一個(gè)懸臂梁受集中荷載的問(wèn)題。圖1為本案例Matlab編程計(jì)算得到的結(jié)果。主要內(nèi)容涉及四面體單元的有限元基本理論的推導(dǎo),主要是單元?jiǎng)偠染仃嚨耐茖?dǎo),此外還包括等參單元和Hammer數(shù)值積分以及三維問(wèn)題的后處理計(jì)算。
圖1 懸臂梁受集中荷載的應(yīng)力云圖
一個(gè)完整的有限元程序基本組成部分包括前處理模塊、分析主程序模塊和后處理模塊。在前處理模塊中,實(shí)現(xiàn)節(jié)點(diǎn)坐標(biāo)輸入、單元節(jié)點(diǎn)編號(hào)、網(wǎng)絡(luò)劃分以及邊界條件輸入等工作;在分析主程序模塊中,求解整體剛度方程;在后處理模塊中,實(shí)現(xiàn)結(jié)果顯示、數(shù)據(jù)輸出等工作。
對(duì)應(yīng)的有限元法的基本步驟:
(1)幾何域離散,獲得標(biāo)準(zhǔn)化的單元; (2)通過(guò)能量原理(虛功原理或最小勢(shì)能原理,獲得單元?jiǎng)偠确匠蹋?(3)單元的集成(裝配); (4)處理位移邊界條件; (5)計(jì)算位移場(chǎng);
(6)計(jì)算單元的其他物理量(應(yīng)力應(yīng)變)。
這幾步中,最核心的內(nèi)容是單元研究,具體包括:
(1)節(jié)點(diǎn)描述(不同坐標(biāo)系節(jié)點(diǎn)坐標(biāo)的變化); (2)場(chǎng)描述(位移場(chǎng),應(yīng)變場(chǎng),應(yīng)力場(chǎng),形函數(shù));
(3)單元?jiǎng)偠确匠蹋ɑ谀芰吭硗茖?dǎo))。
需要說(shuō)明的是后文的四面體單元有限元方程的推導(dǎo)過(guò)程是基于等參單元的基本理論從局部坐標(biāo)(自然坐標(biāo)、體積坐標(biāo))出發(fā)來(lái)推導(dǎo)四面體單元的剛度矩陣,因?yàn)檫@樣做比較規(guī)范自然,推導(dǎo)過(guò)程也適用于其他類型單元。但是因?yàn)樗拿骟w單元相對(duì)簡(jiǎn)單也可以直接從直角坐標(biāo)(全局坐標(biāo))進(jìn)行推導(dǎo),具體推導(dǎo)過(guò)程可參考清華大學(xué)曾攀老師的課程,直接從直角坐標(biāo)(全局坐標(biāo))進(jìn)行推導(dǎo)的過(guò)程省去了等參單元雅各比矩陣呀等坐標(biāo)系映射的各種概念,理解起來(lái)相對(duì)容易。
公式(1)-(3)為彈性力學(xué)中三維空間彈性問(wèn)題的完整描述,分別是空間問(wèn)題的平衡方程、幾何方程、物理方程,這是我們推導(dǎo)有限元方程的基礎(chǔ)。這些公式會(huì)在后面的有限元方程推導(dǎo)過(guò)程中用到。
四面體單元的坐標(biāo)描述涉及了等參單元的概念,在有限元方法中,若要離散邊界為曲線或曲面的求解域,需要建立將形狀規(guī)則的單元變換為邊界為曲線或曲面的單元的方法,在有限元法中對(duì)應(yīng)此問(wèn)題所采用的變換方法是等參變換,即單元幾何形狀的變換和單元內(nèi)場(chǎng)函數(shù)采用相同數(shù)目的節(jié)點(diǎn)及相同的插值函數(shù)進(jìn)行變換。四面體單元的參數(shù)坐標(biāo)就是體積坐標(biāo),體積坐標(biāo)的定義如圖2所示。
圖2 四面體單元體積坐標(biāo)的定義
參數(shù)坐標(biāo)系下的形函數(shù)等于對(duì)應(yīng)的體積坐標(biāo),四個(gè)節(jié)點(diǎn)對(duì)應(yīng)四個(gè)形函數(shù),如下式,
這樣就實(shí)現(xiàn)了自然坐標(biāo)系和物理坐標(biāo)系下的坐標(biāo)映射,如下式,
同樣形函數(shù)也可以用于物理場(chǎng)的插值,公式(6)是位移場(chǎng)的插值表達(dá)式。
有了位移場(chǎng)插值的表達(dá)式就可以通過(guò)公式(2)-(3)的幾何方程和物理方程推導(dǎo)應(yīng)變場(chǎng)和應(yīng)力場(chǎng)的有限元表達(dá)式,但是公式(2)的幾何方程中存在對(duì)于xy也就是物理坐標(biāo)系下的偏導(dǎo)數(shù),代入公式(6)后也就是形函數(shù)對(duì)物理坐標(biāo)的偏導(dǎo)數(shù),因?yàn)槲覀兊男魏瘮?shù)是在自然坐標(biāo)系下定義的,是關(guān)于ξ, η, ζ的函數(shù),想要知道他對(duì)x,y的偏導(dǎo),這里利用了鏈?zhǔn)角髮?dǎo)法則,即可建立如下式(7)所示的形函數(shù)對(duì)自然坐標(biāo)的物理坐標(biāo)偏導(dǎo)數(shù)的映射關(guān)系。中間的這個(gè)矩陣就叫Jacob矩陣。
如何求Jacob矩陣呢?將公式(5)代入公式(7)可以得到如下式(8)所示的自然坐標(biāo)偏導(dǎo)矩陣乘以一個(gè)坐標(biāo)矩陣,
得到j(luò)acob矩陣后求逆,再乘以自然坐標(biāo)偏導(dǎo)矩陣就可以得到形函數(shù)對(duì)物理坐標(biāo)的導(dǎo)數(shù)。如下式(9),
上述形函數(shù)對(duì)物理坐標(biāo)的導(dǎo)數(shù)的求解過(guò)程對(duì)應(yīng)的matlab代碼如下:
function [NDxyz,JacobiDET] = ShapeFunction(ElementNodeCoordinate) %計(jì)算形函數(shù)及形函數(shù)對(duì)局部坐標(biāo)ksi eta zeta的導(dǎo)數(shù) NDL = [-1 1 0 0;-1 0 1 0;-1 0 0 1];%3*4 [N1Dksi N2Dksi N3Dksi N4Dksi;N1Deta N2Deta N3Deta N4Deta……] Jacobi = NDL*ElementNodeCoordinate;%計(jì)算雅可比矩陣3*4 4*3 JacobiDET = det(Jacobi);%計(jì)算雅可比行列式3*3 [DxDksi DyDksi DzDksi;DxDeta…… JacobiINV=inv(Jacobi);%對(duì)雅可比行列式求逆3*3 NDxyz=JacobiINV*NDL;%利用雅可比行列式的逆計(jì)算形函數(shù)對(duì)結(jié)構(gòu)坐標(biāo)的導(dǎo)數(shù)[DN1Dx DN2Dx DN3Dx;DN1Dy DN2Dy DN3Dy;……] end這樣求出形函數(shù)對(duì)物理坐標(biāo)的導(dǎo)數(shù)后就可以代入公式(2)幾何方程求出應(yīng)變場(chǎng)矩陣B,經(jīng)過(guò)能量原理的推導(dǎo)可以得到單元?jiǎng)偠染仃嚨谋磉_(dá)式,注意剛度矩陣的表達(dá)式是一個(gè)積分的運(yùn)算,由于被積函數(shù)較為復(fù)雜,如果代數(shù)積分進(jìn)行計(jì)算要消耗大量的計(jì)算資源,因此有限元理論中引入數(shù)值積分,即我們熟悉的高斯積分和hammer積分,對(duì)于直角坐標(biāo)系我們采用高斯積分,對(duì)于面積或者體積坐標(biāo)系我們采用hammer數(shù)值積分,具體這兩類積分的原理大家可以參考相關(guān)數(shù)值積分教材即可,這里我們只利用hammer數(shù)值積分的結(jié)論。四面體單元具體的數(shù)值積分公式如下:
其中,對(duì)應(yīng)的B矩陣如下式所示:
單元?jiǎng)偠染仃嚧_定好之后,將其按照自由度索引組裝到全局剛度矩陣中,組裝的核心思想就是確定好每個(gè)單元中的每個(gè)節(jié)點(diǎn)中的每個(gè)自由度在整體剛度矩陣中的位置即可。
整體剛度矩陣確定之后,就是荷載向量p的定義,這個(gè)很簡(jiǎn)單,只要依次確定每個(gè)節(jié)點(diǎn)每個(gè)自由度的外力荷載即可。
剛度矩陣和荷載向量確定后接下來(lái)就是邊界條件的施加,大家最熟悉的可能就是直接劃行劃列的方法,這個(gè)方法的弊端在于劃去自由度為零的行列之后,整體單元?jiǎng)偠染仃嚨木幪?hào)以及自由度編號(hào)都會(huì)改變,對(duì)于大規(guī)模計(jì)算非常不方便,影響計(jì)算效率。為了解決這個(gè)問(wèn)題出現(xiàn)了乘大數(shù)法和置一法,置一法只能解決零邊界問(wèn)題,乘大數(shù)法可以解決零邊界和非零邊界,應(yīng)用非常廣泛。邊界條件的施加這里我們采用乘大數(shù)法。針對(duì)邊界條件采用乘大數(shù)法施加邊界的方式為:在自由度上乘以一個(gè)很大的數(shù)a,然后在相應(yīng)的外力向量對(duì)應(yīng)的自由度上乘以,具體如下式(22)所示
為了驗(yàn)證乘大數(shù)法的可靠性,將對(duì)應(yīng)行列寫成代數(shù)表達(dá)式的形式,如下式所示
考慮到遠(yuǎn)大于,所以公式(13)可以寫成,
即可得到??梢?jiàn)乘大數(shù)法與實(shí)際邊界條件是等價(jià)的。
施加了邊界條件的剛度矩陣K與荷載矩陣p之后,可直接利用matlab中的高斯消去法 計(jì)算符“”,完成位移的求解,即u=Kp。
求出節(jié)點(diǎn)位移之后通常我們做受力分析時(shí)還會(huì)需要知道應(yīng)力應(yīng)變的分布情況,具體通過(guò)下述公式求解:
對(duì)應(yīng)的提取應(yīng)力和應(yīng)變的后處理代碼如下:
% 計(jì)算形函數(shù)導(dǎo)數(shù) [NDxyz, ~] = ShapeFunction(ElementNodeCoordinate);%[DN1Dx DN2Dx DN3Dx;DN1Dy DN2Dy DN3Dy;……] ElementNodeDisplacement=U(ElementNodeDOF);%12*1 節(jié)點(diǎn)位移列陣 ElementNodeDisplacement=reshape(ElementNodeDisplacement,Dof,4);%3*4 % 計(jì)算單元應(yīng)變 Strain3_3 3*3的應(yīng)變矩陣 Strain3_3=ElementNodeDisplacement*NDxyz';%高斯積分點(diǎn)處應(yīng)變 3*4 4*3 %把單元應(yīng)變矩陣改寫成6*1 Strain=[Strain3_3(1,1) Strain3_3(2,2) Strain3_3(3,3) ... Strain3_3(1,2)+Strain3_3(2,1).... Strain3_3(2,3)+Strain3_3(3,2) Strain3_3(1,3)+Strain3_3(3,1)]'; Stress(1:6,1) = D*Strain;%高斯積分點(diǎn)處應(yīng)變利用繪圖命令Patch我們可以得到如圖1所示的應(yīng)力云圖。另外計(jì)算得到的變形和位移云圖如圖3-4所示。
圖3 變形前后的網(wǎng)格對(duì)比
圖4 位移云圖
審核編輯:劉清
-
matlab
+關(guān)注
關(guān)注
180文章
2957瀏覽量
229937 -
矩陣
+關(guān)注
關(guān)注
0文章
418瀏覽量
34463 -
信號(hào)處理模塊
+關(guān)注
關(guān)注
1文章
6瀏覽量
6990
原文標(biāo)題:案例實(shí)操:四面體單元懸臂梁的Matlab有限元編程過(guò)程講解
文章出處:【微信號(hào):sim_ol,微信公眾號(hào):模擬在線】歡迎添加關(guān)注!文章轉(zhuǎn)載請(qǐng)注明出處。
發(fā)布評(píng)論請(qǐng)先 登錄
相關(guān)推薦
評(píng)論