MATLAB计算平面桁架体系的程序

%用于计算平面桁架结构体系，2013年3月16日10:49
function mainprogram
%输入结构信息，包括GEOArray,ID-Array,IEN-Array,ndofs,MAT,FORCE。
%GEOArray为节点几何参数
%IDArray为节点自由度表
%IENArray为单元的节点参数及单元特性（MAT）号
%ndofs为自由度个度
%MAT为材料矩阵
%FORCE为力矩阵,列向量
%LMArray为单元对应的自由度，行：[单元号 ix iy jx jy]...
%ix,iy,jx,jy对应单元i，j节点四个自由度，数值表示相应的结构自由度。
%KGlobal为整体刚度矩阵
%DISPLACEMET为位移矩阵
%ELEMENTDIS为单元节点力矩阵，行：[单元号 Uix Uiy Ujx Uiy]
%ELEMENTFORCE为单元节点力矩阵，行：[单元号 Fix Fiy Fjx Fiy]
clc;
global GEOArray IENArray ndofs MAT FORCE LMArray KGlobal DISPLACEMET ELEMENTDIS ELEMENTFORCE
%---------------------------------------------------------------------
%输入格式说明
%GEOArray行：[节点号 x y],节点号为1，2，3...
%IDArray行：[节点号 DOFx DOFy]，节点号为1，2，3...
%IENArra行：[单元号 jointi jointj 材料号],单元号为1，2，3...
%MAT行：[材料号 E A],材料号为1，2，3...
%FORCE为列向量
%---------------------------------------------------------------------
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
GEOArray=[1 2 3 4;
0 10 20 10;
0 10 0 0]';
IDArray=[1 0 0;
2 1 2;
3 0 0;
4 3 4];
IENArray=[1 2 3 4 5;
1 1 4 3 4;
2 4 2 2 3;
1 1 1 1 1]' ;
ndofs=4;
MAT=[1 210e9 0.025];
FORCE=[40e6 -30e6 0 0]';
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%---------------------------------------------------------------------
GENLMA(IDArray);%生成单元-自由度矩阵
KGlobal=zeros(ndofs,ndofs);
DISPLACEMET=zeros(ndofs);
[m,n]=size(IENArray);%获取单元个数
%---------------------------------------------------------------------
%绘制结点图形
jointx=GEOArray(:,2);%取出节点坐标
jointy=GEOArray(:,3);
[mn,mm]=max(jointx);%确定节点坐标范围
xmax=jointx(mm);
[mn,mm]=min(jointx);
xmin=jointx(mm);
[mn,mm]=max(jointy);
ymax=jointy(mm);
[mn,mm]=min(jointy);
ymin=jointy(mm);
plot(jointx,jointy,'ro');%绘制节点
hold on;
grid on;
axis equal;%锁定纵横比
axis([xmin-0.1*xmin xmax+0.1*xmax ymin-0.1*ymin ymax+0.1*ymax]);%设定显示坐标范围
%---------------------------------------------------------------------
%生成单元刚度矩阵并集成
for idm=1:m
ke=ESTIFFNESS(idm);%生成第idm个单元的单元刚度矩阵
ASSEMBLAGEke(idm,ke);%将生成的单元刚度矩阵集成到整体刚度矩上三角元素中
GetKGlobal;%生成完整的整体刚度矩阵
end
%求解位移
FindDis;
%计算单元节点位移
findELEMENTDIS;
%计算单元节点力
findELEMENTFORCE;
%打印结果
printfRESULT;
return;

function GENLMA(IDArray)
%生成单元-自由度矩阵
global LMArray IENArray
LMArray(:,1)=IENArray(:,1);
L=length(LMArray)

;
for i=1:L
LMArray(i,2)=IDArray(IENArray(i,2),2);
LMArray(i,3)=IDArray(IENArray(i,2),3);
LMArray(i,4)=IDArray(IENArray(i,3),2);
LMArray(i,5)=IDArray(IENArray(i,3),3);
end
return;

function ke=ESTIFFNESS(idm)
%提取单元节点及节点坐标
global IENArray MAT GEOArray
jointi=IENArray(idm,2);
jointj=IENArray(idm,3);
x1=GEOArray(jointi,2);
y1=GEOArray(jointi,3);
x2=GEOArray(jointj,2);
y2=GEOArray(jointj,3);
%------------------------------------------------------------------
%绘制单元
jointx=[x1 x2];
jointy=[y1 y2];
plot(jointx,jointy,'k-','LineWidth',2);
hold on;
%------------------------------------------------------------------
%得到单元倾角和整体坐标系下单元刚度矩阵
xy=[x2-x1 y2-y1];%生成单元向量
length=sqrt((xy(1))^2+(xy(2))^2);%单元长度
sinfi=xy(2)/length;
cosfi=xy(1)/length;
%不考虑材料特性的单元刚度举证
ke=[cosfi^2 cosfi*sinfi -cosfi^2 -cosfi*sinfi
cosfi*sinfi sinfi^2 -cosfi*sinfi -sinfi^2
-cosfi^2 -cosfi*sinfi cosfi^2 cosfi*sinfi
-cosfi*sinfi -sinfi^2 cosfi*sinfi sinfi^2];
%提取材料特性
E=MAT(IENArray(idm,4),2);
A=MAT(IENArray(idm,4),3);
%生成整体坐标系下单元刚度矩阵
ke=E*A/length*ke;
return;

function ASSEMBLAGEke(idm,ke)
global KGlobal LMArray
%将单元对非0的自由度贡献的刚度集成到整体刚度矩阵中
for i=2:5
ixy=LMArray(idm,i);%按顺序搜索非0自由度
if ixy~=0
for j=i:5
jxy=LMArray(idm,j);%后续非0自由度
if jxy~=0
if ixy>jxy
KGlobal(jxy,ixy)=ke(i-1,j-1)+KGlobal(jxy,ixy);
else
KGlobal(ixy,jxy)=ke(i-1,j-1)+KGlobal(ixy,jxy);
end
end
end
end
end
return;

function GetKGlobal
global KGlobal
[m,n]=size(KGlobal);
for i=2:m
for j=1:i-1
KGlobal(i,j)=KGlobal(j,i);%将上三角元素复制到下三角
end
end
return;

function FindDis
global KGlobal DISPLACEMET FORCE
DISPLACEMET=KGlobal\FORCE;%KD=F;D=K\F。求解位移
return;

function findELEMENTDIS
%计算单元节点位移
global LMArray DISPLACEMET ELEMENTDIS
ELEMENTDIS(:,1)=LMArray(:,1);
L=length(ELEMENTDIS);
for i=1:L
for j=2:5
degj=LMArray(i,j);%提取单元节点的自由度
if LMArray(i,j)~=0
ELEMENTDIS(i,j)=DISPLACEMET(degj);%提取对应自由度的位移
else
ELEMENTDIS(i,j)=0;
end
end
end

return;

function findELEMENTFORCE
%计算单元节点力
global LMArray ELEMENTDIS ELEMENTFORCE
ELEMENTFORCE(:,1)=LMArray(:,1);
L=length(ELEMENTFORCE);
for idm=1:L
ke=ESTIFFNESS(idm);%生成第idm个单元的单元刚度矩阵
ELEMENTFORCE1=ke*(ELEMENTDI

S(idm,2:5))';%生成第idm个单元的节点力
ELEMENTFORCE(idm,2:5)=ELEMENTFORCE1';%集成
end
return;

function printfRESULT
global KGlobal DISPLACEMET ELEMENTDIS ELEMENTFORCE

% 检查数据文件是否已存在，确保无重复文件名，避免数据丢失
% 创建模型数据库文件
% 检查文件是否创建成功
flag=0;%可输入状态标志，0为不可输入，1为可以输入
blag=fopen('work.dat','r');
if blag==-1
fid_in=fopen('work.dat','w');
if fid_in==-1;
sprintf( '错误：文件创建失败！！！' )
end
fclose(fid_in);
flag=1;
else
blag2=2;
blag2=input( sprintf( '错误：work.dat文件已存在，是否覆盖？ ( Yes=1 / No=2 ): ') ) ;
if blag2==1
fid_in=fopen('work.dat','w');
if fid_in==-1;
sprintf( '错误：文件创建失败！！！' )
end
else
disp(sprintf( '\n请移出work.dat文件，然后再次运行本程序 \n\n' ))
end
fclose(fid_in);
flag=1;
end
if flag==1
fid_in=fopen('work.dat','w');
if fid_in==-1;
sprintf( '错误：文件打开失败！！！' )
else
fprintf(fid_in,'%s\n','用于计算平面桁架结构体系(MATLAB2007b)，2013年3月16日10:49');
fprintf(fid_in,'%s\n\n\n\n','---------------------------------------------------------');
fprintf(fid_in,'%s\n','整体刚度矩阵');
fprintf(fid_in,'%s\n','KGlobal=');
%打印整体刚度矩阵
[imax,jmax]=size(KGlobal);
for i=1:imax
for j=1:jmax
fprintf(fid_in, '%16.5e ', KGlobal(i,j));
end
fprintf(fid_in,'\n');
end
%打印整体位移矩阵
fprintf(fid_in, '\n\n\n%s\n', '整体位移矩阵');
fprintf(fid_in, '%s\n', 'DISPLACEMET= [自由度号 U]');
imax=length(DISPLACEMET);
for i=1:imax
fprintf(fid_in, '%6d %16.5e\n', i, DISPLACEMET(i));
end
%打印单元节点力
fprintf(fid_in, '\n\n\n%s\n', '单元节点位移矩阵');
fprintf(fid_in, '%s\n', 'ELEMENTDIS= [单元号 Uix Uiy Ujx Ujy]');
[imax,jmax]=size(ELEMENTDIS);
for i=1:imax
fprintf(fid_in, '%6d ', ELEMENTDIS(i,1));
for j=2:jmax
fprintf(fid_in, '%16.5e ', ELEMENTDIS(i,j));
end
fprintf(fid_in,'\n');
end
%打印单元节点力
fprintf(fid_in, '\n\n\n%s\n', '单元节点力矩阵');
fprintf(fid_in, '%s\n', 'ELEMENTFORCE= [单元号 Fix Fiy Fjx Fjy]');
[imax,jmax]=size(ELEMENTFORCE);
for i=1:imax
fprintf(fid_in, '%6d ', ELEMENTFORCE(i,1));
for j=2:jmax
fprintf(fid_in, '%16.5e ', ELEMENTFORCE(i,j));
end
fprintf(fid_in,'\n');
end
end
fclose(fid_in);

end
return;