四節點/八節點四邊形單元懸臂梁的Matlab有限元編程——《Matlab有限元編程從入門到精通》系列
導讀：大家好，我是SimPC博士，主要從事工程結構抗震及減隔震研究，玻璃成型熱工設備流動及傳熱研究，玻璃材料力學性能研究。精通有限元等數值算法的實現，有限元軟件二次開發，數據處理，偏微分方程求解，優化算法，GUI界面開發等。有多項科研成果，其中SCI論文4篇，EI3篇，專利2篇。
本文的案例主要以受均布荷載和集中荷載的變截面懸臂梁為研究對象，通過matlab編制四節點和八節點四邊形單元有限元程序來對懸臂梁進行受力分析，提供對應有限元基本理論講解的同時展示相應代碼的實現技巧。
幾何域離散，獲得標準化的單元；
通過能量原理（虛功原理或最小勢能原理，獲得單元剛度方程；
單元的集成（裝配）；
處理位移邊界條件；
計算支反力；
計算單元的其他物理量（應力應變）。
節點描述
場描述
-
單元剛度方程。
圖2-1 平面四節點矩形單元的映射關系
(4)
和
為自然坐標系下的節點坐標值。
(7)
function N=ShapeFun(s,t) N1=1/4*(1-s)*(1-t);N2=1/4*(1 s)*(1-t);N3=1/4*(1 s)*(1 t);N4=1/4*(1-s)*(1 t);N=[N1 0 N2 0 N3 0 N4 0;0 N1 0 N2 0 N3 0 N4];end
同理平面八節點矩形單元如圖2-3所示，單元共有8個節點，16個自由度。單元的形函數定義如下：
(8)
(9)
(10)
和
為自然坐標系下的節點坐標值。
(11)
(12)
(13)
function N=ShapeFun(s,t) %% 四邊形八結點等參單元形函數矩陣 % 角點N1=1/4*(1-s)*(1 t)*(-s t-1); N2=1/4*(1-s)*(1-t)*(-s-t-1); N3=1/4*(1 s)*(1-t)*(s-t-1); N4=1/4*(1 s)*(1 t)*(s t-1); % 邊中點 N5=1/2*(1-t^2)*(1-s); N7=1/2*(1-t^2)*(1 s); N6=1/2*(1-s^2)*(1-t); N8=1/2*(1-s^2)*(1 t); N=[N1 0 N2 0 N3 0 N4 0 N5 0 N6 0 N7 0 N8 0;0 N1 0 N2 0 N3 0 N4 0 N5 0 N6 0 N7 0 N8];
(14)
(15)
(16)
(17)
(18)
function J=Jacobi(ie,s,t,Elements,Nodes) ENodes = Elements(ie,:); %獲取單元結點 xe = Nodes(ENodes(:),:); %獲取節點坐標 x1=xe(1,1);y1=xe(1,2); x2=xe(2,1);y2=xe(2,2); x3=xe(3,1);y3=xe(3,2); x4=xe(4,1);y4=xe(4,2); J=1/4*[-(1 t) -(1-t) 1-t 1 t;1-s -(1-s) -(1 s) 1 s]*[x1 y1;x2 y2;x3 y3;x4 y4];end
function J=Jacobi(ie,kesi,yita,Elements,Nodes) ENodes = Elements(ie,:); %獲取單元結點 xe = Nodes(ENodes(:),:); %獲取結點坐標 x1=xe(1,1);y1=xe(1,2); x2=xe(2,1);y2=xe(2,2); x3=xe(3,1);y3=xe(3,2); x4=xe(4,1);y4=xe(4,2); J=1/4*[-(1-yita),(1-yita),(1 yita),-(1 yita);-(1-kesi),-(1 kesi),(1 kesi),(1-kesi)]*[x1 y1;x2 y2;x3 y3;x4 y4];end
為了求出上述平面四節點和八節點單元的單元剛度矩陣，需要借助能量原理（虛功原理、最小勢能原理）進行推導，能量原理的推導過程大家可以參考任意一本有限元理論書籍，都會有詳細的推導過程，這里就不做進一步推導講解，直接給出物理坐標和幾何坐標系下的剛度矩陣的公式
(19)
(20)
其中B矩陣為應變矩陣，如下式；D矩陣為材料剛度矩陣，如公式（1）所示，是物理方程中表征應力應變關系的矩陣。從上述剛度矩陣的表達式可以看出，自然坐標和物理坐標間要完成坐標映射、偏導映射、面積隱射三個部分，具體映射公式已在上一節的等參單元講解中詳細給出。
(21)
4、高斯積分
公式（20）中的單元剛度矩陣通過數值積分求得，本案例中的四節點和八節點四邊形等參單元均采用2*2個積分點的高斯積分即可求得精確結果。高斯積分點的坐標具體如圖所示。
4-1 Gauss積分點示意圖
公式(20)寫成數值積分的形式為
(22)
對于8節點單元實現上述數值積分的代碼如下所示：
r = [-sqrt(1/3) sqrt(1/3)]; % 2*2 高斯積分點 s = [r(1) r(1) r(2) r(2)]; t = [r(2) r(1) r(1) r(2)]; % 高斯積分點坐標for i=1:4 J = Jacobi(E_ID,s(i),t(i),Elements,Nodes); % 雅可比矩陣 Nst = DiffShapeFun(s(i),t(i)); % 形函數關于自然坐標s,t求導 Nxy = zeros(8,2); for j=1:8 Nxy(j,:) = (J\Nst(j,:)')'; % 形函數關于 x,y 求導=inv(J)*Nst end Bm = [Nxy(1,1) 0 Nxy(2,1) 0 Nxy(3,1) 0 Nxy(4,1) 0 Nxy(5,1) 0 Nxy(6,1) 0 Nxy(7,1) 0 Nxy(8,1) 0; 0 Nxy(1,2) 0 Nxy(2,2) 0 Nxy(3,2) 0 Nxy(4,2) 0 Nxy(5,2) 0 Nxy(6,2) 0 Nxy(7,2) 0 Nxy(8,2); Nxy(1,2) Nxy(1,1) Nxy(2,2) Nxy(2,1) Nxy(3,2) Nxy(3,1) Nxy(4,2) Nxy(4,1) Nxy(5,2) Nxy(5,1) Nxy(6,2) Nxy(6,1) Nxy(7,2) Nxy(7,1) Nxy(8,2) Nxy(8,1)]; ke = ke det(J)*Bm'*D*Bm*Width; %數值積分 end
5、均布荷載的施加
在有限元中分布力要轉為等效節點荷載，而確定等效節點荷載的方法也是通過能量原理推導得到
(22)
上式中，第一項代表體積力的等效荷載，第二項代表面積力的等效荷載，這個案例我們只考慮面力荷載。實現公式22的代碼為
function Pe=UniLoad(ie,N_ID_p1,q0,Nodes,Elements) k=-0.625e-3; % 均布荷載值 N/mms = [-sqrt(1/3) sqrt(1/3)]; % 2*2 高斯積分點ENodes = N_ID_p1(ie,:); %獲取單元結點號Pe=zeros(16,1); %生成臨時單元節點力零列向量x1=Nodes(ENodes(1),1);x6=Nodes(ENodes(4),1);L16=abs(x6-x1); %單元長度for i=1:2 %用于高斯積分的求和循環 N_q=ShapeFun(s(i),1); % 4級子程序：ShapeFun(s(i),1) q_x=q0; Pe=Pe N_q'*q_x*[0;L16/2]; endend
三、Matlab有限元編程精品課
網格劃分及變形結果如圖3-1所示。本案例的詳細視頻教程和對應的matlab源碼，請關注我的技術鄰官網和APP精品課程《Matlab有限元編程從入門到精通10講》,點擊試看《Matlab有限元編程從入門到精通》。
圖3-1 梁變形結果
此外，為幫助大家更好的入門學習Matlab有限元編程分析能力，筆者準備了了一些有限元編程書籍資料供大家下載。歡迎大家直接在附件下載。另外本人推出《Matlab有限元編程從入門到精通》視頻教學課程。點擊試看《Matlab有限元編程從入門到精通》。
本課程為matlab有限元編程專題課，課程主要以案例的形式進行講解，中間會穿插案例中所涉及到的有限元基本理論，案例不局限于力學問題的有限元求解，還會涉及傳熱學、電學等問題的有限元求解。
因為固體力學領域我最熟悉，所以我們從固體力學開始，所涉及的單元有桿單元，梁單元，平面三角形單元，薄板單元，厚板單元，四面體實體單元等等，力學問題有靜力學問題，也有動力學問題，后期還會涉及材料非線性、幾何非線性、接觸非線性等非線性問題，內容豐富，不斷更新完善。
此外，筆者為所有訂閱用戶提供知識圈答疑服務和VIP用戶交流群。并附贈課程相關資料等（平臺支持自行開具電子發票）。
-
快速獲得各典型有限元案例的 Matlab代碼； -
學習并掌握有限元基礎理論； -
掌握Matlab編程實現有限元算法的流程； -
掌握多種有限元單元的基本理論Matlab編程實現過程； -
掌握靜力學、動力學、材料非線性、幾何非線性、接觸非線性問題的Matlab編程實現； 為訂閱用戶提供知識圈答疑服務，并建立VIP用戶交流群，后續可根據訂閱用戶需求進行加餐直播。此外還提供課程對應的學習資料模型一份。
理工科院校學生和教師；
學習型仿真設計工程師；
Matlab有限元編程興趣愛好者和應用者。
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
