有限元基础教程——2D杆单元（Matlab算例）

---

前言

有限元基础教程（曾攀）Matlab3.2.5算例

函数文件

---

将该文件保存为Bar2D2NodeFunc.m文件

function Bar2D2Node=Bar2D2NodeFunc
Bar2D2Node.Stiffness=@Bar2D2Node_Stiffness;
Bar2D2Node.Assembly=@Bar2D2Node_Assembly;
Bar2D2Node.Stress=@Bar2D2Node_Stress;
Bar2D2Node.Forces=@Bar2D2Node_Forces;
end
%%%%% Bar2D2Node %%begin%%%%
function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)
%该函数计算单元的刚度矩阵
%输入弹性模量E、横截面积A
%输入第一个节点坐标(x1,y1)、第二个节点坐标(x2,y2)和角度alpha(单位是度)
%输出单元刚度矩阵k(4*4)
%----------------
L=sqrt((x2-x1)^2+(y2-y1)^2);
x=alpha*pi/180;
C=cos(x);
S=sin(x);
k=E*A/L*[C*C C*S -C*C -C*S;C*S S*S -C*S -S*S;-C*C -C*S C*C C*S;-C*S -S*S C*S S*S];
end
%%%%%%%%%%%%%%%%%%%
function z=Bar2D2Node_Assembly(KK,k,i,j)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k，单元的节点编号i，j
%输出整体刚度矩阵KK
%-----------------
DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
for n1=1:4
    for n2=1:4
        KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
    end
end
z=KK;
end
%---------------
function stress=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)
%该函数计算单元的应力
%输入弹性模量E、第一个节点坐标(x1,y1)和第二个节点坐标(x2,y2)
%输入角度alpha(单位是度)和单位节点位移向量u
%返回单元应力值stress
%------------------------
L=sqrt((x2-x1)^2+(y2-y1)^2);
x=alpha*pi/180;
C=cos(x);
S=sin(x);
stress=E/L*[-C -S C S]*u;
end
%--------------------
%%%%%%%%%%%%%%%%%%%%%%%%%
function forces=Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u)
%该函数计算单元节点力向量
%输入弹性模量E和横截面积A
%输入第一个节点坐标(x1,y1)、第二个节点坐标(x2,y2)和角度alpha(单位是度)
%输入单元节点位移向量u
%返回单元节点力forces
%----------------------------
L=sqrt((x2-x1)^2+(y2-y1)^2);
x=alpha*pi/180;
C=cos(x);
S=sin(x);
forces=E*A/L*[-C -S C S]*u;
end
%%%%%%%%% Bar2D2Node %%end%%%%%%%%%

      

函数调用

E=2.95e11;
A=0.0001;
x1=0;
y1=0;
x2=0.4;
y2=0;
x3=0.4;
y3=0.3;
x4=0;
y4=0.3;
alpha1=0;
alpha2=90;
alpha3=atan(0.75)*180/pi;
Bar2D2Node=Bar2D2NodeFunc;
k1=Bar2D2Node.Stiffness(E,A,x1,y1,x2,y2,alpha1)
k2=Bar2D2Node.Stiffness(E,A,x2,y2,x3,y3,alpha2)
k3=Bar2D2Node.Stiffness(E,A,x1,y1,x3,y3,alpha3)
k4=Bar2D2Node.Stiffness(E,A,x4,y4,x3,y3,alpha1)
KK=zeros(8);
KK=Bar2D2Node.Assembly(KK,k1,1,2);
KK=Bar2D2Node.Assembly(KK,k2,2,3);
KK=Bar2D2Node.Assembly(KK,k3,1,3);
KK=Bar2D2Node.Assembly(KK,k4,4,3)

k=KK([3 5 6],[3 5 6])
p=[20000;0;-25000];
u=k\p

q=[0 0 u(1) 0 u(2) u(3) 0 0]';
P=KK*q

u1=[q(1);q(2);q(3);q(4)]
stress1=Bar2D2Node.Stress(E,x1,y1,x2,y2,alpha1,u1)
u2=[q(3);q(4);q(5);q(6)]
stress2=Bar2D2Node.Stress(E,x2,y2,x3,y3,alpha2,u2)
u3=[q(1);q(2);q(5);q(6)]
stress3=Bar2D2Node.Stress(E,x1,y1,x3,y3,alpha3,u3)
u4=[q(7);q(8);q(5);q(6)]
stress4=Bar2D2Node.Stress(E,x4,y4,x3,y3,alpha1,u4)

运行结果

除此之外，若运行过程可以循环，可以通过eval函数简化（该算例不适用）

总结

将多个函数文件保存在同一个m文件中，同时若计算过于繁琐且有规律则可以使用eval函数