计算力学平面桁架程序说明.docx

1、计算力学平面桁架程序说明题目：编制一个能够计算任意形状的平面桁架（线弹性）的有限元程序一、基本思想1、先建立一个平面桁架的力学模型，再把结构整体拆开，分解成若干个单元（在杆件结构中就把每个杆件取作一个单元，基础数据有：结点坐标、结点编号、单元编号、单元的节点号、单元的属性、节点荷载、节点自由度）。2、单元分析：采用刚度影响系数建立单元节点的平衡方程（即单元刚度方程）。3、整体分析：把所有单元的刚度方程组合成整体的刚度方程，这是一组以节点物理量为未知量的线性方程组，引入边界条件求解该方程组。4、计算位移和内力二、程序实现整个程序的流程：1、主程序%整个执行过程（总控）process.mclear

2、;InFileName=input(请输入基础数据文件名(默认为：model data.txt）：,s);if isempty(InFileName) InFileName=model data.txt;end Joint,Elem,fidout = Struss_ReadData( InFileName );for i=1:length(Elem.Def(:,1) Len(i),Orient(i,:)=LCS(Elem.Def(i,:),Joint.Coord);endStructStiff=strus_Stiff( Joint,Elem,Len,Orient);%计算刚度Node_Disp

3、=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=varargin1表示输出类型的编号：输出单元长度1、输出单元位向2、输出单元长度及位向3%输出：Len1=单元长度% Orient2=单元位向%输入参数处

4、理if length(varargin)=0 TypeNo=3;elseif length(varargin)=1 TypeNo=varargin1;else error(调用函数LCS时，输入参数数目有误！)end%单元长度和位向计算xy1=Coord(iDef(2),:);xy2=Coord(iDef(3),:);dxy=xy2-xy1;Len=sqrt(sum(dxy.*dxy);if TypeNo=1 varargout1=Len;elseif TypeNo=2 varargout1=dxy(1)/Len,dxy(2)/Len;elseif TypeNo=3 varargout1=Le

5、n; varargout2=dxy(1)/Len,dxy(2)/Len;end2、从基础数据文件读取数据赋值给数组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=

6、PackData(Joint,Elem,JointDef);%把基础数据写入输出文件OutFileName,fidout=WriteData(Joint,Elem,InfileName);end*function Joint,Elem,JointDef = ReadData( Joint,Elem,InFileName )%从基础数据文件读取数据fidin=fopen(InFileName,r);%以只读方式打开格式文件JointDef=;%设置初值if fidin=-1 error(没有这个基础数据文件);endwhile feof(fidin); %测试文件的结尾 line=fgetl(f

7、idin);%按行读取字符串 line=deblank(line(end:-1:1);line(end:-1:1)=line;%去掉每行字符前的空格 ifisempty(line)&strncmp(line,%,1)%排除空行和注释行 %-读取结点和单元数据 KeyWord=strtok(line,);%取第一个逗号“，”前的子串，即关键字 dotsuffix=find(line=,);%提取逗号在字符串中的下标 NumVec=str2num(line(dotsuffix(1)+1:end);%提取第一个逗号后的子串并数值化 switch KeyWord case Contl Joint.NJ

8、oint=NumVec(1);Joint.NDOF=NumVec(2);Elem.NElem=NumVec(3); case Joint JointDef=JointDef;NumVec; case Elem Elem.Def=Elem.Def;NumVec; caseJointLoad Joint.Load=Joint.Load;NumVec; otherwise error(没有这种数据类型标识！) end endend fclose(fidin);%关闭文件 *function Joint,Elem,JointDef = PackData( Joint,Elem,JointDef )%整

9、理输入数据% 整理结点坐标数据Joint.Coord=JointDef(:,2:3);%结点坐标Joint.DOF=JointDef(:,4:5);%结点自由度%整理单元数据if length(Elem.Def(1,:)=4 Elem.Def(:,5)=1e-5;%设置线膨胀系数默认值end*function OutFileName,fidout = WriteData( Joint,Elem,InfileName )%WriteData把基础数据写入输出文件%设置输出文件名，把.m替换为_Out.txtOutFileName=strrep(InfileName,.txt,_Out.txt);

10、fidout=fopen(OutFileName,wt);%以只写方式打开格式文件%基础数据输出到文件fprintf(fidout,%sn,%平面桁架静力分析数据);fprintf(fidout,%sn,%输入数据文件：,InfileName);fprintf(fidout,n);fprintf(fidout,%sn,%-输入数据-);fprintf(fidout,%sn,%控制数据);fprintf(fidout,%sn,%单元数 结点数 自由度数);fprintf(fidout,%5d%10d%10dn,Elem.NElem,Joint.NJoint,length(find(Joint.D

11、OF);fprintf(fidout,n);fprintf(fidout,%sn,%基础数据);fprintf(fidout,%sn,%结点坐标及自由度);fprintf(fidout,%sn,%结点号 X Y DOF1 DOF2 );for i=1:Joint.NJoint fprintf(fidout,%5d%9g%9g%9d%9dn,i,Joint.Coord(i,:),Joint.DOF(i,:);endfprintf(fidout,n);fprintf(fidout,%sn,%单元编号、截面几何和材料常数);fprintf(fidout,%sn,%单元号 始端号 末端号 轴向刚度 线

12、膨胀系数);for i=1:Elem.NElem fprintf(fidout,%5d%10d%12g%15g%15gn,Elem.Def(i,:);end%输出结点荷载fprintf(fidout,n);fprintf(fidout,%sn,%结点荷载);if isempty(Joint.Load) fprintf(fidout,%sn,%结点号 方向 荷载值); r,v=size(Joint.Load); for i=1:r for j=1:v fprintf(fidout,%dt,Joint.Load(i,j); end fprintf(fidout,n ); end else fpri

13、ntf(fidout,%sn,%-无结点荷载);endend3、单元分析求出单元在整体坐标系中的刚度矩阵（这里直接求整体坐标下的刚度矩阵还是比较方便的，不需要先求单元坐标下的矩阵，再转换成整体坐标下的）。 4、整体分析采用单元集成法，分别考虑每个单元对F的单独贡献，然后进行叠加。单元集成法求整体刚度矩阵的步骤可表示为： 。这里，在单元刚度矩阵与整体刚度矩阵之间增添了一个中间环节单元贡献矩阵。function StructStiff = strus_Stiff( Joint,Elem,Len,Orient )%计算单元刚度矩阵，集成结构刚度矩阵% Joint=结点信息，Elem=单元信息% St

14、ructStiff=结构刚度矩阵StructStiff=zeros(2*Joint.NJoint);for i=1:Elem.NElem %计算单元刚度矩阵 ElemStiff=Struss_Elem_StiffGlobal(Elem.Def(i,:),Len(i),Orient(i,:); %集成单元刚度矩阵 StructStiff=AssemElem(Elem.Def(i,:),ElemStiff,StructStiff);endend*function K = Struss_Elem_StiffGlobal( iDef,L,Orient )%Struss_Elem_StiffGlobal

15、函数是计算一个单元在结构坐标系（不是局部坐标系）中的单元刚度% 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

16、)=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);end5、结构坐标系中结点位移计算function Nod

17、e_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,%sn,%-计算结果:若计算结果与经验、其他专业软件相差较多，请仔细检查基本数据文件-);r,v3=size(StrucStiff);fprintf(fidout,%s%d%s%dn,%结构刚度,

18、r,v3);for i=1:r for j=1:v3 fprintf(fidout,%19.5e,StrucStiff(i,j); end fprintf(fidout,n);endJointLoad=zeros(Joint.NJoint*2,1);for i=1:length(TotalLoad(:,1) JointLoad(TotalLoad(i,1)*2-2+TotalLoad(i,2)=TotalLoad(i,3); if length(TotalLoad(i,:)3 JointLoad(TotalLoad(i,1)*2-2+TotalLoad(i,4)=TotalLoad(i,5);

19、 endendJointLoad(v)=;Disp=StrucStiffJointLoad;Node_Disp(v2)=Disp;Node_Disp(v)=0;Node_Disp=transpose(Node_Disp);fprintf(fidout,n);fprintf(fidout,%sn,%-结构反应-);fprintf(fidout,%sn,%结构坐标系中的结点位移：水平位移u、竖向位移v);fprintf(fidout,%sn,%结点号 u v );for i=1:Joint.NJoint fprintf(fidout,%5d%19.5e%21.5en,i,Node_Disp(2*i

20、-1),Node_Disp(2*i);endend5、杆件轴力计算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);for i=1:Elem.NElem ElemStiff=Struss_Elem_StiffGlobal(Elem.Def(i,:),Len(i),Orient(i,:); InFL(i,:)=ElemStiff*Node_Disp(

21、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); if Orient(i,1)*InFL(i,1)0|Orient(i,2)*InFL(i,2)0 F(i)=-F(i); end endfprintf(fidout,n);fprintf(fidout,%sn,%结构坐标系中的单元杆端内力、单元轴力);fprintf(fidout,%sn,%

22、单元号 F1x F1y F2x F2y F );for i=1:Elem.NElem fprintf(fidout,%5d%19.5e%19.5e%19.5e%19.5e%19.5en,i,InFL(i,1),InFL(i,2),InFL(i,3),InFL(i,4),F(i);endfclose(fidout);disp(disp(计算完毕，要查看结果请打开【基础数据文件名+_Out.txt】文件);end三、算例如上图，将每个结点、单元编号（带括号的数字为单元编号）。Matlab中该程序执行过程： process请输入基础数据文件名(默认为：model data.txt）：ex1.txt计

23、算完毕，要查看结果请打开【基础数据文件名+_Out.txt】文件ex1_Out.txt文件中的部分内容：%结点号 u v 1 -2.73438e-005 0.00000e+000 2 0.00000e+000 0.00000e+000 3 -1.26168e-004 2.14844e-005 4 -4.28117e-004 -1.17187e-005 5 -6.63660e-004 -4.49219e-005 6 -9.33779e-004 3.20623e-004 7 -1.21309e-003 -7.42187e-005 8 -1.27168e-003 -2.85156e-004 9 -9

24、.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%结构坐标系中的单元杆端内力、单元轴力%单元号 F1x F1y F2x F2y F 1 -2.33333e+003 0.00000e+000 2.33333e+003 0.00000e+000 2.33333e+003 2 0.00000e

25、+000 -1.83333e+003 0.00000e+000 1.83333e+003 1.83333e+003 3 2.33333e+003 2.33333e+003 -2.33333e+003 -2.33333e+003 -3.29983e+003 4 0.00000e+000 8.66667e+003 0.00000e+000 -8.66667e+003 -8.66667e+003 5 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 6 0.00000e+000 8.66667e+003 0.00000e

26、+000 -8.66667e+003 -8.66667e+003 7 2.33333e+003 -2.33333e+003 -2.33333e+003 2.33333e+003 3.29983e+003 8 0.00000e+000 2.83333e+003 0.00000e+000 -2.83333e+003 -2.83333e+003 9 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 0.00000e+000 10 0.00000e+000 -4.00000e+003 0.00000e+000 4.00000e+003 -4.00000e+003 11 2.33333e+003 2.33333e+003 -2.33333e+003 -2.33333e+003 -3.29983e+003