计算力学平面桁架程序说明.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
升级会员 