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

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

• 單元的集成（裝配）；

• 處理位移邊界條件；

• 計算支反力；

• 計算單元的其他物理量（應(yīng)力應(yīng)變）。

• 節(jié)點描述

• 場描述

• 單元剛度方程。

圖2-1 平面四節(jié)點矩形單元的映射關(guān)系

? ? ?



(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

同理平面八節(jié)點矩形單元如圖2-3所示，單元共有8個節(jié)點，16個自由度。單元的形函數(shù)定義如下：

? (8)

? ? ? ?(9)

? (10)

? ? ?(11)

? ? ? ?(12)

?? (13)

function N=ShapeFun(s,t)            %% 四邊形八結(jié)點等參單元形函數(shù)矩陣 % 角點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,:);                 %獲取單元結(jié)點 xe = Nodes(ENodes(:),:);                 %獲取節(jié)點坐標 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,:);                 %獲取單元結(jié)點 xe = Nodes(ENodes(:),:);                 %獲取結(jié)點坐標 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

為了求出上述平面四節(jié)點和八節(jié)點單元的單元剛度矩陣，需要借助能量原理（虛功原理、最小勢能原理）進行推導(dǎo)，能量原理的推導(dǎo)過程大家可以參考任意一本有限元理論書籍，都會有詳細的推導(dǎo)過程，這里就不做進一步推導(dǎo)講解，直接給出物理坐標和幾何坐標系下的剛度矩陣的公式

? ?(19)

? ?(20)

其中B矩陣為應(yīng)變矩陣，如下式；D矩陣為材料剛度矩陣，如公式（1）所示，是物理方程中表征應(yīng)力應(yīng)變關(guān)系的矩陣。從上述剛度矩陣的表達式可以看出，自然坐標和物理坐標間要完成坐標映射、偏導(dǎo)映射、面積隱射三個部分，具體映射公式已在上一節(jié)的等參單元講解中詳細給出。

? ? ? ?(21)

4、高斯積分

公式（20）中的單元剛度矩陣通過數(shù)值積分求得，本案例中的四節(jié)點和八節(jié)點四邊形等參單元均采用2*2個積分點的高斯積分即可求得精確結(jié)果。高斯積分點的坐標具體如圖所示。

4-1 Gauss積分點示意圖

公式(20)寫成數(shù)值積分的形式為

? ? ? ? (22)

對于8節(jié)點單元實現(xiàn)上述數(shù)值積分的代碼如下所示：

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));        % 形函數(shù)關(guān)于自然坐標s,t求導(dǎo)         Nxy = zeros(8,2);         for j=1:8           Nxy(j,:) = (JNst(j,:)')';             % 形函數(shù)關(guān)于 x,y 求導(dǎo)=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;  %數(shù)值積分  end

5、均布荷載的施加

在有限元中分布力要轉(zhuǎn)為等效節(jié)點荷載，而確定等效節(jié)點荷載的方法也是通過能量原理推導(dǎo)得到

? ? ? ?(22)

上式中，第一項代表體積力的等效荷載，第二項代表面積力的等效荷載，這個案例我們只考慮面力荷載。實現(xiàn)公式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,:);                    %獲取單元結(jié)點號Pe=zeros(16,1);                          %生成臨時單元節(jié)點力零列向量x1=Nodes(ENodes(1),1);x6=Nodes(ENodes(4),1);L16=abs(x6-x1);                          %單元長度for i=1:2                                 %用于高斯積分的求和循環(huán)    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有限元編程精品課

網(wǎng)格劃分及變形結(jié)果如圖3-1所示。本案例的詳細視頻教程和對應(yīng)的matlab源碼，請關(guān)注我的仿真秀官網(wǎng)和APP精品課程《Matlab有限元編程從入門到精通10講》。

圖3-1 梁變形結(jié)果