刚架矩阵位移法matlab计算程序

clc;

%输入基本信息

jd=input('请输入节点数量');

free=input('请输入自由度');

gj=input('请输入杆件数量');

P=zeros(free,1);

P=input('请输入外荷载矩阵（[x;y;z]）');

d=zeros(free,1);

member=zeros(1,6,gj);

K=zeros(6,6,gj);

k=zeros(6,6,gj)

for i=1:1:gj

angle(i)=input(sprintf('请输入第%d杆件角度（角度制）',i));

E(i)=input(sprintf('请输入第%d杆件弹性模量（N/mm2）',i));

A(i)=input(sprintf('请输入第%d杆件截面面积（mm2）',i));

L(i)=input(sprintf('请输入第%d杆件长度（mm）',i));

member(:,:,i)=input(sprintf('请输入第%d杆件定位向量([1 2 3 4 5 6])',i)); end

%计算局部坐标系下单元刚度矩阵k

for i=1:1:gj

k=(:,:,i)=E(i)*I(i)/L(i)^3*[A(i)*L(i)^2/I(i) 0 0 -A(i)*L(i)^2/I(i) 0 0;

0 12 6*L(i) 0 -12 6*L(i);

0 6*L(i) 4*L(i)^2 0 -6*L(i) 2*L(i)^2;

-A(i)*L(i)^2/I(i) 0 0 A(i)*L(i)^2/I(i) 0 0;

0 -12 -6*L(i) 0 12 -6*L(i);

0 6*L(i) 2*L(i)^2 0 -6*L(i) 4*L(i)^2 ];

end

%计算各杆件转换矩阵T

T=zeros(6,6,gj);

for i=1:1:gj

T=(:,:,i)=[cosd(angle(i)) sind(angle(i)) 0 0 0 0;

-sind(angle(i)) cosd(angle(i)) 0 0 0 0;

0 0 1 0 0 0;

0 0 0 cosd(angle(i)) sind(angle(i)) 0;

0 0 0 -sind(angle(i)) cosd(angle(i)) 0;

0 0 0 0 0 1];

end

%计算整体坐标系下单元刚度矩阵K

K(:,:,i)=T(:,:,i)'*k(:,:,i)*T(:,:,i);

end

%形成S

Ss=zeros(free);

S=zeros(free);

for i=1:1:gj

for n=1:1:6

for j=1:1:6

if (member(1,n,i)<=free && member(1,j,i)<=free)

Ss(member(1,n,i),member(1,j,i))=K(n,j,i);

end

end

end

S=S+Ss ;

Ss=zeros(free);

end

%计算杆间荷载作用Pf

Ff=[FAbcosd-FSbsind;

FAbsind+FSbcosd;

FMb;

FAecosd-FSbsind;

FAesind+FSecosd;

FMe]

for i=1:1:gj

gjl=input('请输入杆间力数量');

ppp=zeros(1,4,gjl);

for j=1:1:gjl

ppp(:,:,i)=input(sprintf('请输入第%d杆件l1,l2,W，a ([1 2 3 4])',i));

l1=ppp(1,1,i);

l2=ppp(1,2,i);

W=ppp(1,3,i);

a=ppp(1,3,i);

if a==0,l1+l2==L(i)

pd=1

else if a==0,l1+l2

pd=2

else if a>0,l1+l2==L(i)

pd=3

else

pf=4

end

end

end

switch pd

case 1

Qf(:,:,gjl)=[

0;

W*l2^2*(3*l1+l2);

W*l1*l2^2;

0;

W*l1^2*(l1+3*l2);

-W*l1^2*l2/L(i)^2]

case 2

Qf(:,:,gjl)=[

W*cosd(a);

W*sind(a)*l2^2*(3*l1+l2);

W*sind(a)*l1*l2^2;

W*cosd(a);

W*sin(a)*l1^2*(l1+3*l2);

-W*l1^2*l2/L(i)^2]

case 3

Qf(:,:,gjl)=[

0;

W*L(i)/2*(1-l1/L(i)^4*(2*L(i)^3-2*l1^2*L(i)+l1^3)-l2^3/L(i)^4*(2*L(i)-l2));

W*L(i)^2/12*(1-l1^2/L(i)^4*(6*L(i)^2-8*l1*L(i)+3*l1^2)-l2^3/L(i)^4*(4*L(i)-3*l2));

W*L(i)/2*(1-l1^3/L(i)^4*(2*L(i)-l1)-l2/L(i)^4*(2*L(i)^3-2*l2^2+l2^3));

0;

-W*L(i)^2/12*(1-l1^3/L(i)^4*(4*L(i)-3*l1)-l2^2/L(i)^4*(6*L(i)^2-8*l2*L(i)+3*l2^2))

]

case 4

Qf(:,:,gjl)=[

W*cosd(a)*(L(i)-l1-l2);

W*sind(a)*L(i)/2*(1-l1/L(i)^4*(2*L(i)^3-2*l1^2*L(i)+l1^3)-l2^3/L(i)^4*(2*L(i)-l2));

W*sind(a)*L(i)^2/12*(1-l1^2/L(i)^4*(6*L(i)^2-8*l1*L(i)+3*l1^2)-l2^3/L(i)^4*(4*L(i)-3*l2));

W*sind(a)*L(i)/2*(1-l1^3/L(i)^4*(2*L(i)-l1)-l2/L(i)^4*(2*L(i)^3-2*l2^2+l2^3));

W*sind(a)*cosd(a)*(L(i)-l1-l2);