梁单元有限元计算程序(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('Displacement=')

disp(Displacement)

主程序用到的函数：

计算单元刚度矩阵

function k=Beam_Stiffness(E,A,I,x1,y1,x2,y2,alpha)

%计算梁长度

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));

%计算局部单元刚度矩阵

ke=[E*A/L,0,0,-E*A/L,0,0;...

0,12*E*I/(L^3),6*E*I/(L^2),0,-12*E*I/(L^3),6*E*I/(L^2);...

0,6*E*I/(L^2),4*E*I/L,0,-6*E*I/(L^2),2*E*I/L;...

-E*A/L,0,0,E*A/L,0,0;...

0,-12*E*I/(L^3),-6*E*I/(L^2),0,12*E*I/(L^3),-6*E*I/(L^2);...

0,6*E*I/(L^2),2*E*I/L,0,-6*E*I/(L^2),4*E*I/L];

%计算变换矩阵

x=alpha*pi/180;

C=cos(x);

S=sin(x);

Te=[C S 0;-S C 0;0 0 1];

T=[Te,zeros(3,3);zeros(3,3),Te];

k=T'*ke*T;

k=tril(k);

k=sparse(k);

组集总刚度矩阵

function z=Beam_Assembly(KK,k,i,j)

DOF(1)=3*i-2;

DOF(2)=3*i-1;

DOF(3)=3*i;

DOF(4)=3*j-2;

DOF(5)=3*j-1;

DOF(6)=3*j;

for n1=1:6

for n2=1:6

KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);

end

end

z=KK;

%将该稀疏矩阵中元素放入整体刚度矩阵中的相应位置，

%使得整体刚度矩阵仍然为下三角阵形式的稀疏矩阵；

z=z+triu(z,1)'-triu(z,1);

z=sparse(z);