有限元分析基础教程-曾攀 ex1

 

 

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;


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)
k=[k1;k2;k3;k4]
csvwrite('D:\Matlab\FEM1.csv',k,1,1)

KK=zeros(8,8)
KK=Bar2D2Node_Assembly (KK,k1,1,2)
csvwrite('D:\Matlab\FEM01.csv',KK,1,1)
KK=Bar2D2Node_Assembly (KK,k2,2,3)
csvwrite('D:\Matlab\FEM02.csv',KK,1,1)
KK=Bar2D2Node_Assembly (KK,k3,1,3)
csvwrite('D:\Matlab\FEM03.csv',KK,1,1)
KK=Bar2D2Node_Assembly (KK,k4,4,3)
csvwrite('D:\Matlab\FEM04.csv',KK,1,1)

k=KK([3,5,6],[3,5,6])

csvwrite('D:\Matlab\k.csv',k,1,1)

p=[20000,0,25000]';
u=k\p

q=[0 0 0.0002712 0 0.0000565 -0.0002225 0 0]';
P=KK*q ;

u1=[q(1);q(2);q(3);q(4);];
u2=[q(3);q(4);q(5);q(6);];
u3=[q(1);q(2);q(5);q(6);];
u4=[q(7);q(8);q(5);q(6);];

stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha1,u1)
stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,alpha2,u2)
stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,alpha3,u3)
stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,alpha1,u4)

 

function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)
    
L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
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];

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;

 

function stress= Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)
%该函数计算单元的应力
%输入弹性模量 E，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2）
%输入角度 alpha（单位是度）和单位节点位移矢量 u
%返回单元应力标量 stress
L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
x=alpha*pi/180;
C=cos(x);
S=sin(x);
stress=E/L*[-C -S C S]*u;