matlab计算桁架,MATLAB计算平面桁架体系的程序

global l;

E = 1e6;

A = 1;

p = fopen('INPUTDATA.txt','r');

NP = fscanf(p,'% d',1);

NF = fscanf(p,'% d',1);

MP = zeros(NP,2);

for i =1;NP

MP(i,:) = fscanf(p,'% d',2);

end

NE = fscanf(p,'% d',1);

ME = zeros(NE,2);

for i =1:NE

ME(i,:) = fscanf(p,'% d',2)

end

NR = fscanf(p,'% d',1);

NRR = zeros(NR,2);

for i =1:NR

NRR(i,:) = fscanf(p,'% d',2);

end

NL = fscanf(p,'% d',1);

LL = zeros(NL,2);

for i =1:NL

LL(i,:) = fscanf(p,'% d',2);

end

NDF = NF * NP;

K = zeros(NDF,NDF);

for i =1:NE

l(i) = sqrt((MP(ME(i,2),1) - MP(ME(i,1),1))^2 + (MP(ME(i,2),2) - MP(ME(i,1),2))^2);

v = [MP(ME(i,2),1) - MP(ME(i,1),1),MP(ME(i,2),2) - MP(ME(i,1),2)]/l(i);

ax = v * [1;0];ay = v * [0;1];

T(:,:,i) = [ax,ay,0,0;0,0,ax,ay];

ke = (E * A/l(i)) * [1,-1;-1,1];

Ke = T(:,:,i)' * ke * T(:,:,i);

K(2 * ME(i,1) - 1:2 * ME(i,1),2 * ME(i,1) - 1:2 * ME(i,1))=...

K(2 * ME(i,1) - 1:2 * ME(i,1),2 * ME(i,1) - 1:2 * ME(i,1))+Ke(1:2,1:2);

K(2 * ME(i,1) - 1:2 * ME(i,1),2 * ME(i,2) - 1:2 * ME(i,2))=...

K(2 * ME(i,1) - 1:2 * ME(i,1),2 * ME(i,2) - 1:2 * ME(i,2))+Ke(1:2,3:4);

K(2 * ME(i,2) - 1:2 * ME(i,2),2 * ME(i,1) - 1:2 * ME(i,1))=...

K(2 * ME(i,2) - 1:2 * ME(i,2),2 * ME(i,1) - 1:2 * ME(i,1))+Ke(3:4,1:2);

K(2 * ME(i,2) - 1:2 * ME(i,2),2 * ME(i,2) - 1:2 * ME(i,2))=...

K(2 * ME(i,2) - 1:2 * ME(i,2),2 * ME(i,2) - 1:2 * ME(i,2))+Ke(3:4,3:4);

end

P = zeros(2 * NP,1);

for i = 1:NL

P(LL(i,1)) = LL(i,2);

end

for i = 1:NR

j = NRR(i,1);

K(j,j) = 10^30;

end

for i =1:NR

j = NRR(i,1);

P(j) = 10^30 * NRR(i,2);

end

u = K^( -1) * P;

fprintf('各自由度位移为:/n');

fprintf('% 1.4f/n',u);

fprintf('各杆单元内力为:/n');

for i = 1:NE

D = (E/(1-0.3^2))*[1,0.3,0;0.3,1,0;0,0,(1 - 0.3)/2];

B = [ -1/l(i),0,1/l(i),0;0,0,0,0;0, -1/l(i),0,1/l(i)];

uii = [ u( 2 * ME(i,1));u(2 * ME(i,1));u(2 * ME(i,2) - 1);u(2 * ME(i,2))];

ui = T(:,:,i) * uii;

F =E * (ui(2)-ui(1))/l(i);

fprintf('% 3.4f/n',F);

end