e=2e11;

area = 1e-4;

l1=0.4;

l2=0.3;

l3 = 0.5;

%==========基本刚度矩阵=============

k1 = e*area/l1* [1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0];

k2 = e*area/l2*[ 0 0 0 0;0 1 0 -1;0 0 0 0;0 -1 0 1];

k3 = e*area/l3*[0.64 -0.48 -0.64 0.48;-0.48 0.36 0.48 -0.36;-0.64 0.48 0.64 -0.48;0.48 -0.36 -0.48 0.36];

%==========组装刚度矩阵=============

k1_temps = e*area/l1* [1 0 -1 0 0 0;

    0 0 0 0 0 0;

    -1 0 1 0 0 0;

    0 0 0 0 0 0;

    0 0 0 0 0 0;

    0 0 0 0 0 0]; % 扩展矩阵

k2_temps = e*area/l2*[0 0 0 0 0 0;

    0 0 0 0 0 0;

    0 0 0 0 0 0;

    0 0 0 1 0 -1;

    0 0 0 0 0 0;

    0 0 0 -1 0 1];

k3_temps = e*area/l3*[0.64 -0.48 0 0 -0.64 0.48;

    -0.48 0.36 0 0 0.48 -0.36;

    0 0 0 0 0 0;

    0 0 0 0 0 0;

    -0.64 0.48 0 0 0.64 -0.48;

    0.48 -0.36 0 0 -0.48 0.36];

k = k1_temps+k2_temps+k3_temps;

a= [k(1,1) k(1,2) k(1,6);k(2,1) k(2,2) k(2,6);k(6,1) k(6,2) k(6,6)];

b=[0;-10;0];

x=a\b*1000;   % u1 v1 v3

% yingli1 = e/l1*[-1 1]*zhuanzhi(-pi)*[x(1);x(2);0;0]/10^6; %杆1的应力 单位MPA

yingli1 = e/l1*[-1 1]*zhuanzhi(0)*[0;0;x(1);x(2)]/10^6; %杆1的应力 单位MPA

yingli2 = e/l2*[-1 1]*zhuanzhi(pi/2)*[0;0;0;x(3)]/10^6;

yingli3 = e/l3*[-1 1]*zhuanzhi(-asin(0.6))*[0;x(3);x(1);x(2)]/10^6;

node_weiyi = strcat('节点1 x方向位移为',num2str(x(1)),'mm,y方向位移为',num2str(x(2)),'mm,节点2 x方向位移为0mm,y方向位移为0mm,','节点2 x方向位移为0mm,','y方向位移为',num2str(x(3)),'mm');

disp(node_weiyi)

node_weiyi = strcat('杆1的应力为',num2str(yingli1),'MPa,杆2的应力为',num2str(yingli2),'MPa,杆3的应力为',num2str(yingli3),'MPa,');

disp(node_weiyi)

disp('应力拉为正,压为负')


matlab平面桁架受力分析_有限元