平面桁架结构matlab

桁架结构计算第四章P56
******************************************************************************* function y=plane_truss_element_stiffness(E,A,L,theta) %平面桁架单元刚度
x=theta*pi/180;
C=cos(x);
S=sin(x);
y=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 y=plane_truss_assemble(K,k,i,j) %平面桁架组装
K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);
K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2);
K(2*i-1,2*j-1)=K(2*i-1,2*j-1)+k(1,3);
K(2*i-1,2*j)=K(2*i-1,2*j)+k(1,4);
K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);
K(2*i,2*i)=K(2*i,2*i)+k(2,2);
K(2*i,2*j-1)=K(2*i,2*j-1)+k(2,3);
K(2*i,2*j)=K(2*i,2*j)+k(2,4);
K(2*j-1,2*i-1)=K(2*j-1,2*i-1)+k(3,1);
K(2*j-1,2*i)=K(2*j-1,2*i)+k(3,2);
K(2*j-1,2*j-1)=K(2*j-1,2*j-1)+k(3,3);
K(2*j-1,2*j)=K(2*j-1,2*j)+k(3,4);
K(2*j,2*i-1)=K(2*j,2*i-1)+k(4,1);
K(2*j,2*i)=K(2*j,2*i)+k(4,2);
K(2*j,2*j-1)=K(2*j,2*j-1)+k(4,3);
K(2*j,2*j)=K(2*j,2*j)+k(4,4);
y=K;
*******************************************************************************
function y=plane_truss_element_force(E,A,L,theta,u)%力的表达式
x=theta*pi/180;
C=cos(x);
S=sin(x);
y=E*A/L*[-C -S C S]*u;
*******************************************************************************
function y=plane_truss_element_stress(E,L,theta,u) %应力表达式
x=theta*pi/180;
C=cos(x);
S=sin(x);
y=E/L*[-C -S C S]*u;
***************************************************************************************************** *****************************************************************************************************
clear;
clc;
E= ;
A= ;
L1= ;%杆长度
L2= ;
L3= ;
L4= ;%杆长度
L5= ;
L6= ;
L7= ;
L8= ;
t1= ;%杆角度
t2= ;
t3= ;
t4= ;%杆角度
t5= ;
t6= ;
t7= ;
t8= ;
k1=plane_truss_element_stiffness(E,A,L1,t1) % 单位刚度矩阵
k2=plane_truss_element_stiffness(E,A,L2,t2)
k3=plane_truss_element_stiffness(E,A,L3,t3)
k4=plane_truss_element_stiffness(E,A,L4,t4);% 单位刚度矩阵
k5=plane_truss_element_stiffness(E,A,L5,t5); % 单位刚度矩阵
k6=plane_truss_element_stiffness(E,A,L6,t6)
k7=plane_truss_element_stiffness(E,A,L7,t7)
k8=plane_truss_element_stiffness(E,A,L8,t8) % 单位刚度矩阵
K=zeros( );%组装矩阵大小,节点位移个数Q，节点数乘以2
K=plane_truss_assemble(K,k1, , ); %刚度矩阵的组装
K=plane_truss_assemble(K,k2, , ); %后边为单元连接关系图中杆的指向顺序K=plane_truss_assemble(K,k3, , );
K=plane_truss_assemble(K,k4, , );
K=plane_truss_assemble(K,k5, , ); %刚度矩阵的组装
K=plane_truss_assemble(K,k6, , ); %后边为单元连接关系图中杆的指向顺序
K=plane_truss_assemble(K,k7, , );
K=plane_truss_assemble(K,k8, , ) %注意( )中的数字参数
B=K([ , , , , , , , ,],:) ; %取特定的行5 6 7 8 9 10 11 12
k=B(:,[ , , , , , , , ]) %取特定的列5 6 7 8 9 10 11 12
%得到k矩阵
f=[ ]';
u=k\f
***************************************************************************** ***************************************************************************** 求得位移列矩阵为：
u = [0 0 0 0 0.2133 0.4083 -0.1600 0.4617 0.4267 1.5008 -0.0533 1.6608]’;
以下为求各个杆的应力大小：
u1=[u( );u( );u( );u( )];%后边为单元连接关系图中杆的指向顺序
sigma1=plane_truss_element_stress(E,L1,t1,u1) %element stress，u1是Q列向量
u2=[u( );u( );u( );u( )];
sigma2=plane_truss_element_stress(E,L2,t2,u2)
u3=[u( );u( );u( );u( )];
sigma3=plane_truss_element_stress(E,L3,t3,u3) %element stress，u1是Q列向量
u4=[u( );u( );u( );u( )];
sigma4=plane_truss_element_stress(E,L4,t4,u4)
u5=[u( );u( );u( );u( )];
sigma5=plane_truss_element_stress(E,L5,t5,u5) %element stress，u1是Q列向量
u6=[u( );u( );u( );u( )];
sigma6=plane_truss_element_stress(E,L6,t6,u6)
u7=[u( );u( );u( );u( )];
sigma7=plane_truss_element_stress(E,L7,t7,u7) %element stress，u1是Q列向量
u8=[u( );u( );u( );u( )];
sigma8=plane_truss_element_stress(E,L8,t8,u8)
*******************************************************************************
R=K*u; %求得支反力，从中选取支反力处的值
R=[R( );R( );R( );R( )] %填入位置坐标。