%主程序

function Planetruss44

%计算平面桁架的单元应力和节点位移

Definemodel     %定义模型

Solvemodel      %求解模型

Displayresult   %显示结果

return

%子程序1

function Definemodel

global E m n YSZS jdls theta YS Element F L A

open('tbb.txt');

V=load('tbb.txt');

E=V(1);m=V(2);n=V(3);YSZS=V(4);jdls=V(5);

theta =zeros(1,m);YS= zeros(YSZS,3);Element=zeros(m,2);F=zeros(jdls,3);

L=zeros(1,m);A=zeros(1,m);

theta=[V(6):V(5+m)];YS=[V(6+m):V(6+m+3*YSZS-1)];

Element=[V(6+m+3*YSZS):V(6+m+3*YSZS+2*m-1)];

F=[V(6+m+3*YSZS+2*m):V(6+m+3*YSZS+2*m+jdls*3-1)];

L=[V(6+m+3*YSZS+2*m+jdls*3):V(6+m+3*YSZS+2*m+jdls*3+m-1)];

A=[V(6+m+3*YSZS+2*m+jdls*3+m):V(6+m+3*YSZS+2*m+jdls*3+m+m-1)];

return

%子程序2

function Solvemodel

global YS F K u E m n YSZS jdls theta  Element  L A

%step 1.定义整体刚度矩阵和节点力列阵

K=zeros(2*n);

f=zeros(2*n,1);

%step2.计算单刚并分块投放形成总刚

K=PlaneTrussAssemble(K,m)

function  K=PlaneTrussAssemble(K,m)

global  K  YS F u E m n YSZS jdls theta  Element  L A

for i=1:m

a=theta(i)*pi/180;

c=cos(a);

s=sin(a);

k=E*A(i)/L(i)*[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];

for j=1:2

for k=1:2

for h=1:2

for w=1:2

p=(j-1)*2+h;

q=(k-1)*2+w;

M=(Element(i,j)*2+h);

N=(Element(i,k)*2+w);

K(M,N)=K(M,N)+k(p,q);

end

end

end

end

end

%step3.计算节点力列阵

f=zeros(2*n);

f=PlaneTrussElementForce(n,F)

function f=PlaneTrussElementForce(n,F)

%计算节点力列阵

global f  K  YS F u E m n YSZS jdls theta  Element  L A

f=zeros(2*n,1)

for i=1:n

h=F(i,1);

k=F(i,2);

f((h-1)*2+k)=F(i,3);

end

%step4.

%用置大数法对总刚及节点力列阵进行处理

[f,K]=resolve(YS)

function [f,K]=resolve(YS)

global K f

for  a=1:2*n

b=YS(a,1);

c=YS(a,2);

d=(b-1)*2+c;

f(d)=YS(a,3)*K(m,m)*1e15;

K(m,m)=K(m,m)*1e15;

end

%step5,计算节点位移

u=K\f;

%step6,计算单元应力

es=PlaneTrussElementStress(E,L,theta,u)

global u

for i=1:m

X=theta(i)*Pi/180;

C=cos(x);

S=sin(x);

es(i)=E*A(i)/L(i)*[-C -S C S]*u;

end

return

%子程序3

%输出计算结果

function Displayresult

global Element E A L K YS F u E m n YSZS jdls theta

fprintf('节点位移：\n');

fprintf('节点号  X方向位移  Y方向位移\n');

[n,dummy]=size(2*n);

for i=1:n

fprintf('%6d  %16.8e  %16.8e\n',i,u((i-1)*2+1),u((i-1)*2+2));

end

fprintf('\n\n单元应力\n');

[m,dummy]=size(m);

for ie=1:m

es=PlaneStrussElementStruss(E,L,theta,u);

fprintf('单元号%6d  %16.8e  %16.8e  %16.8e\n',ie,es(1));

end

return