• 幾何域離散，獲得標準化的單元；

• 通過能量原理（虛功原理或最小勢能原理，獲得單元剛度方程；

• 單元的集成（裝配）；

• 處理位移邊界條件；

• 計算支反力；

• 計算單元的其他物理量（應力應變）。

• 節點描述

• 場描述

• 單元剛度方程。

圖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];

? ???(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,:) = (JNst(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講》。

圖3-1 梁變形結果