三角形单元有限元计算程序(matlab)

%此程序计算杆的总刚度矩阵及节点位移，分五步完成：输入各单元数据、计算单元刚度矩阵
%组集总刚度矩阵、计算输出总刚度矩阵 、计算输出节点位移
%2011.4

%输入个单元数据
%输入单元节点编号矩阵，每一行代表该单元的节点编号，即用分号将各单元分开
%用逗号将每个单元内的节点分开，注意输入节点顺序要相同，即各单元要么同为顺时针
%，要么同为逆时针，按单元编号顺序排列
%如[1,2,3;4,3,1]表示两个三角形中的节点编号分别为（1,2,3）、（4,3,1）
cod=input('please input the node of each element in order:');
%计算单元个数，nm为单元个数
[nm,nmn]=size(cod);
%输入节点坐标，每一行代表该节点的坐标，按节点编号顺序排列，即用分号将节点分开
%用逗号将每个节点的坐标分开，按单元编号顺序排列
%如[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 (GPa) of each element in order:');
%输入泊松比miu
miu=input('please input the poison ratio miu of the material:');
%输入厚度t
t=input('please input thickness t (m) of the material:');
%输入节点载荷及相应的节点坐标
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) of the node(if the displacement is unknown please input nan):');
%输入解决问题的性质，平面应力问题性质输入1，平面应变问题输入2
disp('please input the character of the problem:')
ID=input('For Plane stress problem input 1,for Plane strain problem input 2:');

%计算单元刚度矩阵
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);
x3=con(cod(i,3),1);
y3=con(cod(i,3),2);
kele(:,:,i)=Triangle_Stiffness(E,miu,t,x1,y1,x2,y2,x3,y3,ID);%ID输入解决问题的性质
end

%组集总刚度矩阵
%定义空的总刚矩阵
kg=sparse(2*nn,2*nn);%定定义稀疏矩阵，只储存非零元素
%组集各单元矩阵
for n=1:nm
num1=cod(n,1);%提取节点坐标
num2=cod(n,2);
num3=cod(n,3);
kg=Triangle_Assembly(kg,kele(:,:,n),num1,num2,num3);%此函数中实现变半带宽储存
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(Displacement_zero(i))=0;%対应载荷置零
end
Displacement=kg\P';

%输出节点位移
disp('Displacement=')
disp(Displacement)



主函数用到的函数：


计算单元刚度矩阵：
%该函数计算单元的刚度矩阵，输入弹性模量E，泊松比NU，厚度t，三个节点i，j，m
%的坐标xi,yi,xj,yj,xm,ym,平面问题性质参数ID=1时为平面应力问题，ID=2时
%为平面应变问题，最后输出单元刚度矩阵k(6x6)
function k=Triangle_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)
A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
gammam=xj-xi;
B=[betai 0 betaj 0 betam 0;
0 gammai 0 gammaj 0 gammam;
gammai betai gammaj betaj gammam betam]/(2* A);
if ID==1
D=(E/(1-NU*NU))*[1 NU 0 ;NU 1 0 ;0 0 (1-NU)/2];
elseif ID==2
D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0 ;NU 1-NU 0 ;0 0 (1-2*NU)/2];
end
k=t*A*B'*D*B;
k=tril(k);
k=sparse(k);


组集总刚度矩阵：
%该函数进行单元共度矩阵的组装，输入单元刚度矩阵k，单元节点编号i,j,m,
%输出整体刚度矩阵KK
function z= Triangle_Assembly(KK,k,i,j,m)
DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
DOF(5)=2*m-1;
DOF(6)=2*m;
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);