梁单元有限元计算程序(matlab)
%此程序计算杆的总刚度矩阵及节点位移,分五步完成:输入各单元数据、计算单元刚度矩阵
%组集总刚度矩阵、计算输出总刚度矩阵 、计算输出节点位移
%2011.4
%输入个单元数据
%输入单元节点编号矩阵,每一行代表该单元的节点编号,即用分号将各单元分开
%用逗号将每个单元内的节点分开,按单元编号顺序排列
%如[1,2;2,3;1,3]表示三个杆中的节点编号分别为(1,2)、(2,3)、(1,3)
cod=input('please input the node of each element in order:');
%计算单元个数,nm为单元个数
[nm,nmn]=size(cod);
%输入各单元的角度
alpha=input('please input the angle (degree) of each element in order:');
%输入节点坐标,每一行代表该节点的坐标,按节点编号顺序排列,即用分号将节点分开
%用逗号将每个节点的坐标分开,按单元编号顺序排列
%如[1,2;2,3;1,3]表示三个节点的横纵坐标分别为(1,2)、(2,3)、(1,3)
con=input('please input the coordinates (m) of each node in order:');
%计算结点个数,nn为结点个数
[nn,nnn]=size(con);
%输入单元的弹性模量
E=1e9*input('please input E array (GPa) of each element in order:');
%输入单元的截面面积
A=input('please input A array (m^2) of each element in order:');
%输入单元的惯性矩
I=input('please input I array (m^4) of each element in order:');
%输入节点载荷及相应的节点坐标
P=input('please input the force (N) on the node(if the force is unknown please input nan):');
%输入节点位移
Displacement=input('please input the displacement (m or rad) of the node(if the displacement is unknown please input nan):');
%计算单元刚度矩阵
kele=zeros(6,6,nm);
for i=1:nm
x1=con(cod(i,1),1);%提取节点坐标
y1=con(cod(i,1),2);
x2=con(cod(i,2),1);
y2=con(cod(i,2),2);
kele(:,:,i)=Beam_Stiffness(E(i),A(i),I(i),x1,y1,x2,y2,alpha(i));
end
%组集总刚度矩阵
%定义空的总刚矩阵
kg=sparse(3*nn,3*nn);%定定义稀疏矩阵,只储存非零元素
%组集各单元矩阵
for n=1:nm
num1=cod(n,1);
num2=cod(n,2);
kg=Beam_Assembly(kg,kele(:,:,n),num1,num2);%此函数中实现变半带宽储存
end
%输出总刚度矩阵
disp('kg=')
disp(kg)
%计算节点位移
kg=kg+tril(kg,-1)';%把下三角矩阵变为整体矩阵
disp('kg_full=')
disp(full(kg))
%置一赋零法引入边界条件
Displacement_zero=find(Displacement==0);
for i=1:length(Displacement_zero)
kg(Displacement_zero(i),:)=0;%总刚行置零
kg(:,Displacement_zero(i))=0;%总刚列置零
kg(Displacement
_zero(i),Displacement_zero(i))=1;%对角线元素置1
P(Displace
ment_zero(i))=0;%対应载荷置零
end
Displacement=kg\P';
%输出节点位移
disp('Displ