計(jì)算力學(xué)平面桁架程序說(shuō)明

1、題目：編制一個(gè)能夠計(jì)算任意形狀的平面桁架（線彈性）的有限元程序一、基本思想1、先建立一個(gè)平面桁架的力學(xué)模型，再把結(jié)構(gòu)整體拆開，分解成若干個(gè)單元（在桿件結(jié)構(gòu)中就把每個(gè)桿件取作一個(gè)單元，基礎(chǔ)數(shù)據(jù)有：結(jié)點(diǎn)坐標(biāo)、結(jié)點(diǎn)編號(hào)、單元編號(hào)、單元的節(jié)點(diǎn)號(hào)、單元的屬性、節(jié)點(diǎn)荷載、節(jié)點(diǎn)自由度）。2、單元分析：采用剛度影響系數(shù)建立單元節(jié)點(diǎn)的平衡方程（即單元?jiǎng)偠确匠蹋?、整體分析：把所有單元的剛度方程組合成整體的剛度方程，這是一組以節(jié)點(diǎn)物理量為未知量的線性方程組，引入邊界條件求解該方程組。4、計(jì)算位移和內(nèi)力二、程序?qū)崿F(xiàn)整個(gè)程序的流程：基本數(shù)據(jù)（結(jié)點(diǎn)坐標(biāo)、結(jié)點(diǎn)荷載、結(jié)點(diǎn)自由度、單元定義）寫入文件InfileName.t

2、xt供程序讀取，并將這些數(shù)據(jù)寫入_Out.txt文件中，然后調(diào)用函數(shù)LCS ()計(jì)算每個(gè)單元的長(zhǎng)度和方向調(diào)用函數(shù)Struss_ReadData()求結(jié)構(gòu)的整體剛度矩陣，該函數(shù)中還有兩個(gè)子函數(shù)Struss_Elem_StiffGlobal()計(jì)算每個(gè)單元的剛度（單剛），AssemElem()將單剛集成到總剛度矩陣調(diào)用函數(shù)Node_Disp()計(jì)算結(jié)點(diǎn)位移（需刪除自由度為0時(shí)，總剛中的對(duì)應(yīng)行和列），并將位移輸出到文件InfileName_Out.txt中調(diào)用函數(shù)Struss_EndInf()計(jì)算桿件軸力，并將軸力輸出到文件InfileName_Out.txt中1、主程序%整個(gè)執(zhí)行過(guò)程（總控）pro

3、cess.mclear;InFileName=input('請(qǐng)輸入基礎(chǔ)數(shù)據(jù)文件名(默認(rèn)為：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_Stif

4、f( Joint,Elem,Len,Orient);%計(jì)算剛度Node_Disp=Node_Disp(Joint,StructStiff,Joint.Load,fidout);%計(jì)算位移InFL,F=Struss_EndInf(Elem,Len,Orient,Node_Disp,fidout);%計(jì)算內(nèi)力*function varargout = LCS( iDef,Coord,varargin )%該LCS函數(shù)計(jì)算單元的長(zhǎng)度和位向（方向）%輸入：iDef=單元桿端結(jié)點(diǎn)號(hào)% Coord=結(jié)點(diǎn)坐標(biāo)x,y% TypeNo=varargin1表示輸出類型的編號(hào)：輸出單元長(zhǎng)度1、輸出單元位向2、輸出單

5、元長(zhǎng)度及位向3%輸出：Len1=單元長(zhǎng)度% Orient2=單元位向%輸入?yún)?shù)處理if length(varargin)=0 TypeNo=3;elseif length(varargin)=1 TypeNo=varargin1;else error('調(diào)用函數(shù)LCS時(shí)，輸入?yún)?shù)數(shù)目有誤！')end%單元長(zhǎng)度和位向計(jì)算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=d

6、xy(1)/Len,dxy(2)/Len;elseif TypeNo=3 varargout1=Len; varargout2=dxy(1)/Len,dxy(2)/Len;end2、從基礎(chǔ)數(shù)據(jù)文件讀取數(shù)據(jù)賦值給數(shù)組function Joint,Elem,fidout = Struss_ReadData( InfileName )%從基礎(chǔ)數(shù)據(jù)文件讀取數(shù)據(jù)賦值給數(shù)組% Joint=struct('NJoint','NDOF','Coord','DOF','Load',);Elem=struct('NElem&#

7、39;,'Def',);fidout=0;%從基礎(chǔ)數(shù)據(jù)文件讀取數(shù)據(jù)Joint,Elem,JointDef=ReadData(Joint,Elem,InfileName);%整理輸入數(shù)據(jù)Joint,Elem,JointDef=PackData(Joint,Elem,JointDef);%把基礎(chǔ)數(shù)據(jù)寫入輸出文件OutFileName,fidout=WriteData(Joint,Elem,InfileName);end*function Joint,Elem,JointDef = ReadData( Joint,Elem,InFileName )%從基礎(chǔ)數(shù)據(jù)文件讀取數(shù)據(jù)fidin=

8、fopen(InFileName,'r');%以只讀方式打開格式文件JointDef=;%設(shè)置初值if fidin=-1 error('沒(méi)有這個(gè)基礎(chǔ)數(shù)據(jù)文件');endwhile feof(fidin); %測(cè)試文件的結(jié)尾 line=fgetl(fidin);%按行讀取字符串 line=deblank(line(end:-1:1);line(end:-1:1)=line;%去掉每行字符前的空格 ifisempty(line)&strncmp(line,'%',1)%排除空行和注釋行 %-讀取結(jié)點(diǎn)和單元數(shù)據(jù) KeyWord=strtok(l

9、ine,',');%取第一個(gè)逗號(hào)“，”前的子串，即關(guān)鍵字 dotsuffix=find(line=',');%提取逗號(hào)在字符串中的下標(biāo) NumVec=str2num(line(dotsuffix(1)+1:end);%提取第一個(gè)逗號(hào)后的子串并數(shù)值化 switch KeyWord case 'Contl' Joint.NJoint=NumVec(1);Joint.NDOF=NumVec(2);Elem.NElem=NumVec(3); case 'Joint' JointDef=JointDef;NumVec; case '

10、Elem' Elem.Def=Elem.Def;NumVec; case'JointLoad' Joint.Load=Joint.Load;NumVec; otherwise error('沒(méi)有這種數(shù)據(jù)類型標(biāo)識(shí)！') end endend fclose(fidin);%關(guān)閉文件 *function Joint,Elem,JointDef = PackData( Joint,Elem,JointDef )%整理輸入數(shù)據(jù)% 整理結(jié)點(diǎn)坐標(biāo)數(shù)據(jù)Joint.Coord=JointDef(:,2:3);%結(jié)點(diǎn)坐標(biāo)Joint.DOF=JointDef(:,4:5);%

11、結(jié)點(diǎn)自由度%整理單元數(shù)據(jù)if length(Elem.Def(1,:)=4 Elem.Def(:,5)=1e-5;%設(shè)置線膨脹系數(shù)默認(rèn)值end*function OutFileName,fidout = WriteData( Joint,Elem,InfileName )%WriteData把基礎(chǔ)數(shù)據(jù)寫入輸出文件%設(shè)置輸出文件名，把.m'替換為_Out.txt'OutFileName=strrep(InfileName,'.txt','_Out.txt');fidout=fopen(OutFileName,'wt');%以只寫方式

12、打開格式文件%基礎(chǔ)數(shù)據(jù)輸出到文件fprintf(fidout,'%sn','%平面桁架靜力分析數(shù)據(jù)');fprintf(fidout,'%sn','%輸入數(shù)據(jù)文件：',InfileName);fprintf(fidout,'n');fprintf(fidout,'%sn','%-輸入數(shù)據(jù)-');fprintf(fidout,'%sn','%控制數(shù)據(jù)');fprintf(fidout,'%sn','%單元數(shù) 結(jié)點(diǎn)數(shù) 自由度數(shù)

13、9;);fprintf(fidout,'%5d%10d%10dn',Elem.NElem,Joint.NJoint,length(find(Joint.DOF);fprintf(fidout,'n');fprintf(fidout,'%sn','%基礎(chǔ)數(shù)據(jù)');fprintf(fidout,'%sn','%結(jié)點(diǎn)坐標(biāo)及自由度');fprintf(fidout,'%sn','%結(jié)點(diǎn)號(hào) X Y DOF1 DOF2 ');for i=1:Joint.NJoint fprint

14、f(fidout,'%5d%9g%9g%9d%9dn',i,Joint.Coord(i,:),Joint.DOF(i,:);endfprintf(fidout,'n');fprintf(fidout,'%sn','%單元編號(hào)、截面幾何和材料常數(shù)');fprintf(fidout,'%sn','%單元號(hào) 始端號(hào) 末端號(hào) 軸向剛度 線膨脹系數(shù)');for i=1:Elem.NElem fprintf(fidout,'%5d%10d%12g%15g%15gn',Elem.Def(i,:);

15、end%輸出結(jié)點(diǎn)荷載fprintf(fidout,'n');fprintf(fidout,'%sn','%結(jié)點(diǎn)荷載');if isempty(Joint.Load) fprintf(fidout,'%sn','%結(jié)點(diǎn)號(hào) 方向 荷載值'); 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 fprintf(f

16、idout,'%sn','%-無(wú)結(jié)點(diǎn)荷載');endend3、單元分析求出單元在整體坐標(biāo)系中的剛度矩陣（這里直接求整體坐標(biāo)下的剛度矩陣還是比較方便的，不需要先求單元坐標(biāo)下的矩陣，再轉(zhuǎn)換成整體坐標(biāo)下的）。 4、整體分析采用單元集成法，分別考慮每個(gè)單元對(duì)F的單獨(dú)貢獻(xiàn)，然后進(jìn)行疊加。單元集成法求整體剛度矩陣的步驟可表示為： 。這里，在單元?jiǎng)偠染仃嚺c整體剛度矩陣之間增添了一個(gè)中間環(huán)節(jié)單元貢獻(xiàn)矩陣。function StructStiff = strus_Stiff( Joint,Elem,Len,Orient )%計(jì)算單元?jiǎng)偠染仃?，集成結(jié)構(gòu)剛度矩陣% Joint=結(jié)點(diǎn)信

17、息，Elem=單元信息% StructStiff=結(jié)構(gòu)剛度矩陣StructStiff=zeros(2*Joint.NJoint);for i=1:Elem.NElem %計(jì)算單元?jiǎng)偠染仃?ElemStiff=Struss_Elem_StiffGlobal(Elem.Def(i,:),Len(i),Orient(i,:); %集成單元?jiǎng)偠染仃?StructStiff=AssemElem(Elem.Def(i,:),ElemStiff,StructStiff);endend*function K = Struss_Elem_StiffGlobal( iDef,L,Orient )%Struss_E

18、lem_StiffGlobal函數(shù)是計(jì)算一個(gè)單元在結(jié)構(gòu)坐標(biāo)系（不是局部坐標(biāo)系）中的單元?jiǎng)偠? iDef(i,1:3)=單元號(hào)i，始、末端結(jié)點(diǎn)號(hào)、軸向剛度EA%L=單元長(zhǎng)度,Orient=位向EA=iDef(4);K=zeros(4);%初始化單元?jiǎng)偠染仃嘖(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

19、)=-K(2,2);K(4,3)=K(2,1);K=K+triu(K',1);%疊加單元?jiǎng)偠染仃嚿先窃豦nd*function B = AssemElem( iDef,A,B )%AssemElem集成結(jié)構(gòu)剛度矩陣和等效結(jié)點(diǎn)荷載列陣% 輸入：iDef=i單元定義% %A單元?jiǎng)偠染仃?B結(jié)構(gòu)剛度矩陣LocVec=iDef(2)*2-1,iDef(2)*2,iDef(3)*2-1,iDef(3)*2;ij=find(LocVec=0);%這個(gè)這里沒(méi)有用，因?yàn)槲叶技僭O(shè)自由度都有IJ=LocVec(LocVec=0);B(IJ,IJ)=B(IJ,IJ)+A(ij,ij);end5、結(jié)構(gòu)坐標(biāo)

20、系中結(jié)點(diǎn)位移計(jì)算function Node_Disp = Node_Disp( Joint,StrucStiff,TotalLoad,fidout )%該函數(shù)求解結(jié)點(diǎn)位移% v=find(transpose(Joint.DOF)=0);v2=find(transpose(Joint.DOF)=1);StrucStiff(v,:)=;StrucStiff(:,v)=;%輸出總剛fprintf(fidout,'n');fprintf(fidout,'%sn','%-計(jì)算結(jié)果:若計(jì)算結(jié)果與經(jīng)驗(yàn)、其他專業(yè)軟件相差較多，請(qǐng)仔細(xì)檢查基本數(shù)據(jù)文件-');r,v

21、3=size(StrucStiff);fprintf(fidout,'%s%d%s%dn','%結(jié)構(gòu)剛度',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

22、(i,2)=TotalLoad(i,3); if length(TotalLoad(i,:)>3 JointLoad(TotalLoad(i,1)*2-2+TotalLoad(i,4)=TotalLoad(i,5); endendJointLoad(v)=;Disp=StrucStiffJointLoad;Node_Disp(v2)=Disp;Node_Disp(v)=0;Node_Disp=transpose(Node_Disp);fprintf(fidout,'n');fprintf(fidout,'%sn','%-結(jié)構(gòu)反應(yīng)-');fp

23、rintf(fidout,'%sn','%結(jié)構(gòu)坐標(biāo)系中的結(jié)點(diǎn)位移：水平位移u、豎向位移v');fprintf(fidout,'%sn','%結(jié)點(diǎn)號(hào) u v ');for i=1:Joint.NJoint fprintf(fidout,'%5d%19.5e%21.5en',i,Node_Disp(2*i-1),Node_Disp(2*i);endend5、桿件軸力計(jì)算function InFL,F = Struss_EndInf( Elem,Len,Orient,Node_Disp,fidout)%該函數(shù)計(jì)算桿端內(nèi)力

24、% InFL為桿端力（分x、y方向），F(xiàn)為單元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(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)=s

25、qrt(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','%結(jié)構(gòu)坐標(biāo)系中的單元桿端內(nèi)力、單元軸力');fprintf(fidout,'%sn','%單元號(hào) F1x F1y F2x F2y F ');for i=1:Elem.NElem fprintf(

26、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('計(jì)算完畢，要查看結(jié)果請(qǐng)打開【基礎(chǔ)數(shù)據(jù)文件名+_Out.txt】文件');end三、算例如上圖，將每個(gè)結(jié)點(diǎn)、單元編號(hào)（帶括號(hào)的數(shù)字為單元編號(hào)）。Matlab中該程序執(zhí)行過(guò)程：>> process請(qǐng)輸入基礎(chǔ)數(shù)據(jù)文件名(默認(rèn)為：model data.txt）：ex1.txt計(jì)算完畢，要查看結(jié)果請(qǐng)打開【基礎(chǔ)數(shù)據(jù)文

27、件名+_Out.txt】文件ex1_Out.txt文件中的部分內(nèi)容：%結(jié)點(diǎn)號(hào) 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.83621e-004 -7.50616e-004 10 -6.02462e-004 -1.31308e-003 11 -6.15483e-00