计算力学平面桁架程序说明.docx
《计算力学平面桁架程序说明.docx》由会员分享,可在线阅读,更多相关《计算力学平面桁架程序说明.docx(16页珍藏版)》请在冰豆网上搜索。
计算力学平面桁架程序说明
题目:
编制一个能够计算任意形状的平面桁架(线弹性)的有限元程序
一、基本思想
1、先建立一个平面桁架的力学模型,再把结构整体拆开,分解成若干个单元(在杆件结构中就把每个杆件取作一个单元,基础数据有:
结点坐标、结点编号、单元编号、单元的节点号、单元的属性、节点荷载、节点自由度)。
2、单元分析:
采用刚度影响系数建立单元节点的平衡方程(即单元刚度方程)。
3、整体分析:
把所有单元的刚度方程组合成整体的刚度方程,这是一组以节点物理量为未知量的线性方程组,引入边界条件求解该方程组。
4、计算位移和内力
二、程序实现
整个程序的流程:
1、主程序
%整个执行过程(总控)process.m
clear;
InFileName=input('请输入基础数据文件名(默认为:
modeldata.txt):
','s');
ifisempty(InFileName)
InFileName='modeldata.txt';
end
[Joint,Elem,fidout]=Struss_ReadData(InFileName);
fori=1:
length(Elem.Def(:
1))
[Len(i),Orient(i,:
)]=LCS(Elem.Def(i,:
),Joint.Coord);
end
StructStiff=strus_Stiff(Joint,Elem,Len,Orient);%计算刚度
Node_Disp=Node_Disp(Joint,StructStiff,Joint.Load,fidout);%计算位移
[InFL,F]=Struss_EndInf(Elem,Len,Orient,Node_Disp,fidout);%计算内力
******************************************************************************
function[varargout]=LCS(iDef,Coord,varargin)
%该LCS函数计算单元的长度和位向(方向)
%输入:
iDef=单元杆端结点号
%Coord=结点坐标x,y
%TypeNo=varargin{1}表示输出类型的编号:
输出单元长度1、输出单元位向2、输出单元长度及位向3
%输出:
Len{1}=单元长度
%Orient{2}=单元位向
%输入参数处理
iflength(varargin)==0
TypeNo=3;
elseiflength(varargin)==1
TypeNo=varargin{1};
else
error('调用函数LCS时,输入参数数目有误!
')
end
%单元长度和位向计算
xy1=Coord(iDef
(2),:
);
xy2=Coord(iDef(3),:
);
dxy=xy2-xy1;
Len=sqrt(sum(dxy.*dxy));
ifTypeNo==1
varargout{1}=Len;
elseifTypeNo==2
varargout{1}=[dxy
(1)/Len,dxy
(2)/Len];
elseifTypeNo==3
varargout{1}=Len;
varargout{2}=[dxy
(1)/Len,dxy
(2)/Len];
end
2、从基础数据文件读取数据赋值给数组
function[Joint,Elem,fidout]=Struss_ReadData(InfileName)
%从基础数据文件读取数据赋值给数组
%
Joint=struct('NJoint',[],'NDOF',[],'Coord',[],'DOF',[],'Load',[]);
Elem=struct('NElem',[],'Def',[]);
fidout=0;
%从基础数据文件读取数据
[Joint,Elem,JointDef]=ReadData(Joint,Elem,InfileName);
%整理输入数据
[Joint,Elem,JointDef]=PackData(Joint,Elem,JointDef);
%把基础数据写入输出文件
[OutFileName,fidout]=WriteData(Joint,Elem,InfileName);
end
******************************************************************************
function[Joint,Elem,JointDef]=ReadData(Joint,Elem,InFileName)
%从基础数据文件读取数据
fidin=fopen(InFileName,'r');%以只读方式打开格式文件
JointDef=[];%设置初值
iffidin==-1
error('没有这个基础数据文件');
end
while~feof(fidin);%测试文件的结尾
line=fgetl(fidin);%按行读取字符串
line=deblank(line(end:
-1:
1));line(end:
-1:
1)=line;%去掉每行字符前的空格
if~isempty(line)&~strncmp(line,'%',1)%排除空行和注释行
%-----读取结点和单元数据
KeyWord=strtok(line,',');%取第一个逗号“,”前的子串,即关键字
dotsuffix=find(line==',');%提取逗号在字符串中的下标
NumVec=str2num(line(dotsuffix
(1)+1:
end));%提取第一个逗号后的子串并数值化
switchKeyWord
case'Contl'
Joint.NJoint=NumVec
(1);Joint.NDOF=NumVec
(2);Elem.NElem=NumVec(3);
case'Joint'
JointDef=[JointDef;NumVec];
case'Elem'
Elem.Def=[Elem.Def;NumVec];
case'JointLoad'
Joint.Load=[Joint.Load;NumVec];
otherwise
error('没有这种数据类型标识!
')
end
end
end
fclose(fidin);%关闭文件
******************************************************************************
function[Joint,Elem,JointDef]=PackData(Joint,Elem,JointDef)
%整理输入数据
%整理结点坐标数据
Joint.Coord=JointDef(:
2:
3);%结点坐标
Joint.DOF=JointDef(:
4:
5);%结点自由度
%整理单元数据
iflength(Elem.Def(1,:
))==4
Elem.Def(:
5)=1e-5;%设置线膨胀系数默认值
end
******************************************************************************
function[OutFileName,fidout]=WriteData(Joint,Elem,InfileName)
%WriteData把基础数据写入输出文件
%设置输出文件名,把‘.m'替换为’_Out.txt'
OutFileName=strrep(InfileName,'.txt','_Out.txt');
fidout=fopen(OutFileName,'wt');%以只写方式打开格式文件
%基础数据输出到文件
fprintf(fidout,'%s\n','%平面桁架静力分析数据');
fprintf(fidout,'%s\n',['%输入数据文件:
',InfileName]);
fprintf(fidout,'\n');
fprintf(fidout,'%s\n','%-----输入数据---------');
fprintf(fidout,'%s\n','%控制数据');
fprintf(fidout,'%s\n','%单元数结点数自由度数');
fprintf(fidout,'%5d%10d%10d\n',Elem.NElem,Joint.NJoint,length(find(Joint.DOF)));
fprintf(fidout,'\n');
fprintf(fidout,'%s\n','%基础数据');
fprintf(fidout,'%s\n','%结点坐标及自由度');
fprintf(fidout,'%s\n','%结点号XYDOF1DOF2');
fori=1:
Joint.NJoint
fprintf(fidout,'%5d%9g%9g%9d%9d\n',i,Joint.Coord(i,:
),Joint.DOF(i,:
));
end
fprintf(fidout,'\n');
fprintf(fidout,'%s\n','%单元编号、截面几何和材料常数');
fprintf(fidout,'%s\n','%单元号始端号末端号轴向刚度线膨胀系数');
fori=1:
Elem.NElem
fprintf(fidout,'%5d%10d%12g%15g%15g\n',Elem.Def(i,:
));
end
%输出结点荷载
fprintf(fidout,'\n');
fprintf(fidout,'%s\n','%结点荷载');
if~isempty(Joint.Load)
fprintf(fidout,'%s\n','%结点号方向荷载值');
[r,v]=size(Joint.Load);
fori=1:
r
forj=1:
v
fprintf(fidout,'%d\t',Joint.Load(i,j));
end
fprintf(fidout,'\n');
end
else
fprintf(fidout,'%s\n','%--------无结点荷载');
end
end
3、单元分析
求出单元在整体坐标系中的刚度矩阵(这里直接求整体坐标下的刚度矩阵还是比较方便的,不需要先求单元坐标下的矩阵,再转换成整体坐标下的)。
4、整体分析
采用单元集成法,分别考虑每个单元对F的单独贡献,然后进行叠加。
单元集成法求整体刚度矩阵的步骤可表示为:
。
这里,在单元刚度矩阵
与整体刚度矩阵
之间增添了一个中间环节——单元贡献矩阵
。
function[StructStiff]=strus_Stiff(Joint,Elem,Len,Orient)
%计算单元刚度矩阵,集成结构刚度矩阵
%Joint=结点信息,Elem=单元信息
%StructStiff=结构刚度矩阵
StructStiff=zeros(2*Joint.NJoint);
fori=1:
Elem.NElem
%计算单元刚度矩阵
ElemStiff=Struss_Elem_StiffGlobal(Elem.Def(i,:
),Len(i),Orient(i,:
));
%集成单元刚度矩阵
StructStiff=AssemElem(Elem.Def(i,:
),ElemStiff,StructStiff);
end
end
******************************************************************************
function[K]=Struss_Elem_StiffGlobal(iDef,L,Orient)
%Struss_Elem_StiffGlobal函数是计算一个单元在结构坐标系(不是局部坐标系)中的单元刚度
%iDef(i,[1:
3])=单元号i,始、末端结点号、轴向刚度EA
%L=单元长度,Orient=位向
EA=iDef(4);
K=zeros(4);%初始化单元刚度矩阵
K(1,1)=EA/L*Orient
(1)*Orient
(1);
K(2,2)=EA/L*Orient
(2)*Orient
(2);
K(3,3)=K(1,1);
K(4,4)=K(2,2);
K(2,1)=EA/L*Orient
(1)*Orient
(2);
K(3,1)=-K(1,1);
K(4,1)=-K(2,1);
K(3,2)=K(4,1);
K(4,2)=-K(2,2);
K(4,3)=K(2,1);
K=K+triu(K',1);%叠加单元刚度矩阵上三角元素
end
******************************************************************************
function[B]=AssemElem(iDef,A,B)
%AssemElem集成结构刚度矩阵和等效结点荷载列阵
%输入:
iDef=i单元定义
%
%A单元刚度矩阵
%B结构刚度矩阵
LocVec=[iDef
(2)*2-1,iDef
(2)*2,iDef(3)*2-1,iDef(3)*2];
ij=find(LocVec~=0);%这个这里没有用,因为我都假设自由度都有
IJ=LocVec(LocVec~=0);
B(IJ,IJ)=B(IJ,IJ)+A(ij,ij);
end
5、结构坐标系中结点位移计算
function[Node_Disp]=Node_Disp(Joint,StrucStiff,TotalLoad,fidout)
%该函数求解结点位移
%
v=find(transpose(Joint.DOF)==0);
v2=find(transpose(Joint.DOF)==1);
StrucStiff(v,:
)=[];
StrucStiff(:
v)=[];
%输出总刚
fprintf(fidout,'\n');
fprintf(fidout,'%s\n',['%-------------计算结果:
若计算结果与经验、其他专业软件相差较多,请仔细检查基本数据文件-------------------']);
[r,v3]=size(StrucStiff);
fprintf(fidout,'%s%d%s%d\n','%结构刚度',r,'×',v3);
fori=1:
r
forj=1:
v3
fprintf(fidout,'%19.5e',StrucStiff(i,j));
end
fprintf(fidout,'\n');
end
JointLoad=zeros(Joint.NJoint*2,1);
fori=1:
length(TotalLoad(:
1))
JointLoad(TotalLoad(i,1)*2-2+TotalLoad(i,2))=TotalLoad(i,3);
iflength(TotalLoad(i,:
))>3
JointLoad(TotalLoad(i,1)*2-2+TotalLoad(i,4))=TotalLoad(i,5);
end
end
JointLoad(v)=[];
Disp=StrucStiff\JointLoad;
Node_Disp(v2)=Disp;
Node_Disp(v)=0;
Node_Disp=transpose(Node_Disp);
fprintf(fidout,'\n');
fprintf(fidout,'%s\n',['%-------------结构反应-------------------']);
fprintf(fidout,'%s\n','%结构坐标系中的结点位移:
水平位移u、竖向位移v');
fprintf(fidout,'%s\n','%结点号uv');
fori=1:
Joint.NJoint
fprintf(fidout,'%5d%19.5e%21.5e\n',i,Node_Disp(2*i-1),Node_Disp(2*i));
end
end
5、杆件轴力计算
function[InFL,F]=Struss_EndInf(Elem,Len,Orient,Node_Disp,fidout)
%该函数计算杆端内力
%InFL为杆端力(分x、y方向),F为单元i的轴力
InFL=zeros(Elem.NElem,4);
F=zeros(Elem.NElem,1);
fori=1:
Elem.NElem
ElemStiff=Struss_Elem_StiffGlobal(Elem.Def(i,:
),Len(i),Orient(i,:
));
InFL(i,:
)=ElemStiff*[Node_Disp(Elem.Def(i,2)*2-1);Node_Disp(Elem.Def(i,2)*2);Node_Disp(Elem.Def(i,3)*2-1);Node_Disp(Elem.Def(i,3)*2)];
F(i)=sqrt(InFL(i,1)*InFL(i,1)+InFL(i,2)*InFL(i,2));
ifOrient(i,1)*InFL(i,1)>0||Orient(i,2)*InFL(i,2)>0
F(i)=-F(i);
end
end
fprintf(fidout,'\n');
fprintf(fidout,'%s\n','%结构坐标系中的单元杆端内力、单元轴力');
fprintf(fidout,'%s\n','%单元号F1xF1yF2xF2yF');
fori=1:
Elem.NElem
fprintf(fidout,'%5d%19.5e%19.5e%19.5e%19.5e%19.5e\n',i,InFL(i,1),InFL(i,2),InFL(i,3),InFL(i,4),F(i));
end
fclose(fidout);
disp('disp('计算完毕,要查看结果请打开【基础数据文件名+_Out.txt】文件');
end
三、算例
如上图,将每个结点、单元编号(带括号的数字为单元编号)。
Matlab中该程序执行过程:
>>process
请输入基础数据文件名(默认为:
modeldata.txt):
ex1.txt
计算完毕,要查看结果请打开【基础数据文件名+_Out.txt】文件
ex1_Out.txt文件中的部分内容:
%结点号uv
1-2.73438e-0050.00000e+000
20.00000e+0000.00000e+000
3-1.26168e-0042.14844e-005
4-4.28117e-004-1.17187e-005
5-6.63660e-004-4.49219e-005
6-9.33779e-0043.20623e-004
7-1.21309e-003-7.42187e-005
8-1.27168e-003-2.85156e-004
9-9.83621e-004-7.50616e-004
10-6.02462e-004-1.31308e-003
11-6.15483e-004-7.62335e-004
12-6.28504e-004-2.50000e-004
13-4.28117e-004-2.03125e-004
14-1.26168e-004-1.01563e-004
%结构坐标系中的单元杆端内力、单元轴力
%单元号F1xF1yF2xF2yF
1-2.33333e+0030.00000e+0002.33333e+0030.00000e+0002.33333e+003
20.00000e+000-1.83333e+0030.00000e+0001.83333e+0031.83333e+003
32.33333e+0032.33333e+003-2.33333e+003-2.33333e+003-3.29983e+003
40.00000e+0008.66667e+0030.00000e+000-8.66667e+003-8.66667e+003
50.00000e+0000.00000e+0000.00000e+0000.00000e+0000.00000e+000
60.00000e+0008.66667e+0030.00000e+000-8.66667e+003-8.66667e+003
72.33333e+003-2.33333e+003-2.33333e+0032.33333e+0033.29983e+003
80.00000e+0002.83333e+0030.00000e+000-2.83333e+003-2.83333e+003
90.00000e+0000.00000e+0000.00000e+0000.00000e+0000.00000e+000
100.00000e+000-4.00000e+0030.00000e+0004.00000e+003-4.00000e+003
112.33333e+0032.33333e+003-2.33333e+003-2.33333e+003-3.29983e+003