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

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

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

3、cess.mclear;InFileName=input('請輸入基礎(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);%計算剛度Node_Disp=Node_Disp(Joint,StructStiff,Joint.Load,fidout);%計算位移InFL,F=Struss_EndInf(Elem,Len,Orient,Node_Disp,fidout);%計算內(nèi)力*function varargout = LCS( iDef,Coord,varargin )%該LCS函數(shù)計算單元的長度和位向（方向）%輸入：iDef=單元桿端結(jié)點號% Coord=結(jié)點坐標(biāo)x,y% TypeNo=varargin1表示輸出類型的編號：輸出單元長度1、輸出單元位向2、輸出單

5、元長度及位向3%輸出：Len1=單元長度% Orient2=單元位向%輸入?yún)?shù)處理if length(varargin)=0 TypeNo=3;elseif length(varargin)=1 TypeNo=varargin1;else error('調(diào)用函數(shù)LCS時，輸入?yún)?shù)數(shù)目有誤！')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=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('沒有這個基礎(chǔ)數(shù)據(jù)文件');endwhile feof(fidin); %測試文件的結(jié)尾 line=fgetl(fidin);%按行讀取字符串 line=deblank(line(end:-1:1);line(end:-1:1)=line;%去掉每行字符前的空格 ifisempty(line)&strncmp(line,'%',1)%排除空行和注釋行 %-讀取結(jié)點和單元數(shù)據(jù) KeyWord=strtok(l

9、ine,',');%取第一個逗號“，”前的子串，即關(guān)鍵字 dotsuffix=find(line=',');%提取逗號在字符串中的下標(biāo) NumVec=str2num(line(dotsuffix(1)+1:end);%提取第一個逗號后的子串并數(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('沒有這種數(shù)據(jù)類型標(biāo)識！') end endend fclose(fidin);%關(guān)閉文件 *function Joint,Elem,JointDef = PackData( Joint,Elem,JointDef )%整理輸入數(shù)據(jù)% 整理結(jié)點坐標(biāo)數(shù)據(jù)Joint.Coord=JointDef(:,2:3);%結(jié)點坐標(biāo)Joint.DOF=JointDef(:,4:5);%

11、結(jié)點自由度%整理單元數(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é)點數(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é)點坐標(biāo)及自由度');fprintf(fidout,'%sn','%結(jié)點號 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','%單元編號、截面幾何和材料常數(shù)');fprintf(fidout,'%sn','%單元號 始端號 末端號 軸向剛度 線膨脹系數(shù)');for i=1:Elem.NElem fprintf(fidout,'%5d%10d%12g%15g%15gn',Elem.Def(i,:);

15、end%輸出結(jié)點荷載fprintf(fidout,'n');fprintf(fidout,'%sn','%結(jié)點荷載');if isempty(Joint.Load) fprintf(fidout,'%sn','%結(jié)點號 方向 荷載值'); 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','%-無結(jié)點荷載');endend3、單元分析求出單元在整體坐標(biāo)系中的剛度矩陣（這里直接求整體坐標(biāo)下的剛度矩陣還是比較方便的，不需要先求單元坐標(biāo)下的矩陣，再轉(zhuǎn)換成整體坐標(biāo)下的）。 4、整體分析采用單元集成法，分別考慮每個單元對F的單獨貢獻(xiàn)，然后進行疊加。單元集成法求整體剛度矩陣的步驟可表示為： 。這里，在單元剛度矩陣與整體剛度矩陣之間增添了一個中間環(huán)節(jié)單元貢獻(xiàn)矩陣。function StructStiff = strus_Stiff( Joint,Elem,Len,Orient )%計算單元剛度矩陣，集成結(jié)構(gòu)剛度矩陣% Joint=結(jié)點信

17、息，Elem=單元信息% StructStiff=結(jié)構(gòu)剛度矩陣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_E

18、lem_StiffGlobal函數(shù)是計算一個單元在結(jié)構(gòu)坐標(biāo)系（不是局部坐標(biāo)系）中的單元剛度% iDef(i,1:3)=單元號i，始、末端結(jié)點號、軸向剛度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

19、)=-K(2,2);K(4,3)=K(2,1);K=K+triu(K',1);%疊加單元剛度矩陣上三角元素end*function B = AssemElem( iDef,A,B )%AssemElem集成結(jié)構(gòu)剛度矩陣和等效結(jié)點荷載列陣% 輸入：iDef=i單元定義% %A單元剛度矩陣%B結(jié)構(gòu)剛度矩陣LocVec=iDef(2)*2-1,iDef(2)*2,iDef(3)*2-1,iDef(3)*2;ij=find(LocVec=0);%這個這里沒有用，因為我都假設(shè)自由度都有IJ=LocVec(LocVec=0);B(IJ,IJ)=B(IJ,IJ)+A(ij,ij);end5、結(jié)構(gòu)坐標(biāo)

20、系中結(jié)點位移計算function Node_Disp = Node_Disp( Joint,StrucStiff,TotalLoad,fidout )%該函數(shù)求解結(jié)點位移% v=find(transpose(Joint.DOF)=0);v2=find(transpose(Joint.DOF)=1);StrucStiff(v,:)=;StrucStiff(:,v)=;%輸出總剛fprintf(fidout,'n');fprintf(fidout,'%sn','%-計算結(jié)果:若計算結(jié)果與經(jīng)驗、其他專業(yè)軟件相差較多，請仔細(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é)點位移：水平位移u、豎向位移v');fprintf(fidout,'%sn','%結(jié)點號 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、桿件軸力計算function InFL,F = Struss_EndInf( Elem,Len,Orient,Node_Disp,fidout)%該函數(shù)計算桿端內(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','%單元號 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('計算完畢，要查看結(jié)果請打開【基礎(chǔ)數(shù)據(jù)文件名+_Out.txt】文件');end三、算例如上圖，將每個結(jié)點、單元編號（帶括號的數(shù)字為單元編號）。Matlab中該程序執(zhí)行過程：>> process請輸入基礎(chǔ)數(shù)據(jù)文件名(默認(rèn)為：model data.txt）：ex1.txt計算完畢，要查看結(jié)果請打開【基礎(chǔ)數(shù)據(jù)文

27、件名+_Out.txt】文件ex1_Out.txt文件中的部分內(nèi)容：%結(jié)點號 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