有限元的matlab编程

1、有限元编程示例题目描述题目描述: :如下图所示的平面桁架，杆件长度、弹性模量、截面积以如下图所示的平面桁架，杆件长度、弹性模量、截面积以及所受节点力及所受节点力P的大小可以自行定义。求节点位移及杆件轴的大小可以自行定义。求节点位移及杆件轴力。力。例一：桁架例一：桁架解题思路：解题思路： 建立模型建立模型 集成总刚集成总刚 求解位移求解位移 求解杆件轴力求解杆件轴力 输出结果输出结果建立模型：建立模型：定义节点坐标定义节点坐标Node = zeros(10,2) ;x=-1*L; %L为横杆长度为横杆长度for i=1:2:10 x=x+L; Node(i,:)=x 0;endx=-1*L;fo

2、r i=2:2:10 x=x+L; Node(i,:)=x H;%H为竖杆长度为竖杆长度end模型相关参数输入模型相关参数输入H=input(竖杆长度(m):); L=input(水平杆长度(m):);E=input(杆件弹性模量(Gpa):); A=input(杆件截面积(m2):);a=input(节点力P(kN)：);节点编号方式节点编号方式定义单元，即储存单元两端的节点号定义单元，即储存单元两端的节点号Element=zeros(21,2);for i=1:2:7 Element(5/2*i-3/2,:)=i,i+1; Element(5/2*i-1/2,:)=i,i+2; Eleme

3、nt(5/2*i+1/2,:)=i,i+3; endfor i=2:2:8 Element(5*i/2-1,:)=i,i+1; Element(5*i/2,:)=i,i+2;endElement(21,:)=9,10;加下划线的为单元编号加下划线的为单元编号集成总刚：集成总刚：xi = Node( Element( ie, 1 ), 1 ) ;%ie为单元号，以下相同为单元号，以下相同yi = Node( Element( ie, 1 ), 2 ) ;xj = Node( Element( ie, 2 ), 1 ) ;yj = Node( Element( ie, 2 ), 2 ) ;获取单元

4、两端节点坐标获取单元两端节点坐标L = ( (xj-xi)2 + (yj-yi)2 )(1/2) ;计算杆件长度计算杆件长度形成等效荷载列阵形成等效荷载列阵f=0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;%每个节点两个自由度，每个节点两个自由度，a为之前输入的节点力为之前输入的节点力计算从局部坐标到整体坐标的坐标转换矩阵计算从局部坐标到整体坐标的坐标转换矩阵Tfunction T = TransformMatrix( ie )%ie为单元号为单元号c = (xj-xi)/L ;s = (yj-yi)/L ;T= c -s 0 0 s c 0 0 0 0 c

5、-s 0 0 s c ;计算单元刚度矩阵计算单元刚度矩阵kk = E*A/L 0 -E*A/L 0 0 0 0 0 -E*A/L 0 E*A/L 0 0 0 0 0 ; T = TransformMatrix( ie ) ; k = T*k*transpose(T) ;% transpose(T) 为为T的转置矩阵的转置矩阵2集成整体刚度矩阵集成整体刚度矩阵Kfor ie=1:1:21 %按单元顺序进行循环按单元顺序进行循环 k=PlaneTrussElementStiffness(ie); %计算第计算第ie个单元的单刚个单元的单刚 m=Element(ie,1); %ie单元的首节点号单元

6、的首节点号 n=Element(ie,2); %ie单元的末节点号单元的末节点号 K(2*m-1,2*n-1)=k(1,3); K(2*m-1,2*n)=k(1,4); K(2*m,2*n-1)=k(2,3); K(2*m,2*n)=k(2,4);K=zeros(20,20);%用来存储整体刚度矩阵用来存储整体刚度矩阵集成总刚的非对角线元素集成总刚的非对角线元素(这里的元素指这里的元素指2*2的小矩阵）的小矩阵）在下面的集成中，将在下面的集成中，将总刚看成总刚看成10*10的矩阵，每个元素为的矩阵，每个元素为2*2的小矩阵的小矩阵 m=Element(ie,2); %ie单元的末节点号单元的末

7、节点号 n=Element(ie,1); %ie单元的首节点号单元的首节点号 K(2*m-1,2*n-1)=k(3,1); K(2*m-1,2*n)=k(3,2); K(2*m,2*n-1)=k(4,1); K(2*m,2*n)=k(4,2);end 集成总刚的对角线元素（这里的元素指集成总刚的对角线元素（这里的元素指2*2的小矩阵）的小矩阵）for i=1:1:10 %按节点的顺序循环按节点的顺序循环 for j=1:1:21 %对于每个节点，再按单元的顺序循环对于每个节点，再按单元的顺序循环 k=PlaneTrussElementStiffness(j); if Element(j,1)=

8、I %如果如果i节点为节点为j单元的首节点单元的首节点 K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1); K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2); K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1); K(2*i,2*i)=K(2*i,2*i)+k(2,2); end if Element(j,2)=i %如果如果i节点为节点为j单元的末节点单元的末节点 K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(3,3); K(2*i-1,2*i)=K(2*i-1,2*i)+k(3,4); K(2*i,2*i-1)=K(2

9、*i,2*i-1)+k(4,3); K(2*i,2*i)=K(2*i,2*i)+k(4,4); end endend求解位移：求解位移：u=zeros(20);根据约束情况修改总刚，采用对角元素置根据约束情况修改总刚，采用对角元素置1法法for i=1:1:20 K(1,i)=0; K(2,i)=0; K(18,i)=0; K(i,1)=0;K(i,2)=0; K(i,18)=0; end %自由度自由度1、2、18被约束了被约束了,所在的行和列的其他元素都改为所在的行和列的其他元素都改为0K=K*1e15;%乘以一个大数，减小计算误差乘以一个大数，减小计算误差f=f*1e15;u=Kf;求解

10、求解K(1,1)=1;%对角线元素置对角线元素置1K(2,2)=1;K(18,18)=1;求解轴力：求解轴力：获取单元两端的节点号获取单元两端的节点号i = Element( ie, 1 ) ;%ie为单元号为单元号j = Element( ie, 2 ) ;获取单元两端的节点位移获取单元两端的节点位移uElement = zeros( 4, 1 ) ;uElement( 1:2 ) = u( (i-1)*2+1:(i-1)*2+2 ) ;uElement( 3:4 ) = u( (j-1)*2+1:(j-1)*2+2 ) ;k = PlaneTrussElementStiffness( ie

11、 ) ;nodef = k *uElement ;%整体坐标下的节点力整体坐标下的节点力T = TransformMatrix( ie ) ;%计算坐标转换矩阵计算坐标转换矩阵nodef = transpose( T ) * nodef ;%从整体坐标转换到局部坐标从整体坐标转换到局部坐标计算单元的节点力计算单元的节点力输出求解结果：输出求解结果：输出位移输出位移fprintf( 节点位移n ) ;for i=1:1:10 disp(节点号,num2str(i), x方向位移：,num2str(u(2*i-1,1), y方向位移：,num2str(u(2*i,1);end输出节点力输出节点力f

12、printf( nn节点力n ) ; for ie=1:1:21 nodef=NodeForce( ie ); disp( 单元号: ,num2str(ie), 节点号：,num2str(Element(ie,1), 节点号：,num2str(Element(ie,2), 轴力:,num2str(nodef(1) ) ; end例二：网例二：网架架思路分析思路分析网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成网架的基本单元有三角锥，三棱体，正方体，截头四角锥等。鉴于网架的形鉴于网架的形式较多，本程序提供一种通用的网架输入方法，但录入较为繁琐，同时提供式较多，本程序提供一种通

13、用的网架输入方法，但录入较为繁琐，同时提供一种正放四角锥网架的简易输入方法作为典型。一种正放四角锥网架的简易输入方法作为典型。考虑几何非线性考虑几何非线性。本程序采用荷载分级加载的方式考虑网架的几何非线性。将总荷载分成1000份分步施加，求解各荷载步下的节点位移，修改网架相应节点坐标以及刚度矩阵，依次迭代求出网架的总位移。本程序的网架位移求解函数附在主程序后面，主程序运行时调用该函数。几点说明几点说明用户自定义输入用户自定义输入几何建模几何建模正放四角锥网正放四角锥网架简易输入架简易输入定义荷载定义荷载定义边界定义边界条件条件网架分析网架分析位移位移应变应变应力应力位移求解函数位移求解函数单刚

14、矩阵单刚矩阵荷载矩阵荷载矩阵约束矩阵约束矩阵总刚矩阵总刚矩阵求解位移求解位移&分级加载，分级加载，通过修改通过修改节点坐标，节点坐标，迭代求解迭代求解几何非线性几何非线性e=input(选择网架类型，0代表自由定义网架，1代表四角锥网架) %网架类型的选择网架类型的选择网架类型的选择网架类型的选择用户自定义网架（网架信息的录入，包括节点、单元、截面、弹性模量等）用户自定义网架（网架信息的录入，包括节点、单元、截面、弹性模量等）if e=0%选择自定义网架选择自定义网架Node=input(定义节点编号及对应坐标，按1 x1 y1 z1;2 x2 y2 z2;.输入);%形成节点储存矩阵

15、形成节点储存矩阵 Men=input(定义单元与节点的关系，按1 node1 node2;2 node3 node4;.输入,node1node2,依次类推);%形成单元储存矩阵形成单元储存矩阵 Msum=length(Men);%查找网架录入的单元数查找网架录入的单元数 Cont1=input(定义单元实常数，若所有杆件截面面积和弹性模量不变，则输入0，否则输入1);定义单元属性的输入方式定义单元属性的输入方式 if Cont1=0 AE1=input(请输入统一的截面面积与弹性模量，按A E输入); AE=zeros(Msum,3); AE(:,1)=1:Msum;AE(:,2)=AE1(

16、1,1);AE(:,3)=AE1(1,2); else AE=input(请输入相应单元的截面面积与弹性模量，按1,A1 E1;2,A2 E2;.输入); end P=input(定义节点荷载，按node1 P1;node2 P2;.输入); %网架荷载输入网架荷载输入 BC=input(定义边界约束，按node1 Conx Cony Conz;node2 Conx Cony Conz);.输入,Con代表x、y、z方向约束，取0为约束，取1无约束); %网架边界条件网架边界条件end单元属性相同单元属性相同单元属性不同单元属性不同荷载及边界条件荷载及边界条件正放四角锥网架定义正放四角锥网架定

17、义if e=1hu=input(输入网架上层节点行数); %定义网架上层节点的行数定义网架上层节点的行数lu=input(输入网架上层节点列数); %定义网架上层节点的列数定义网架上层节点的列数dis_xu=input(输入网架上层节点列间距); %定义网架上层的行间距定义网架上层的行间距dis_yu=input(输入网架上层节点行间距); %定义网架上层的列间距定义网架上层的列间距hd=hu-1; %网架下层节点的行数网架下层节点的行数ld=lu-1; %网架下层节点的列数网架下层节点的列数dis_xd=dis_xu; %网架下层的行间距网架下层的行间距dis_yd=dis_yu; %网架下

18、层的行间距网架下层的行间距dis_z=input(输入网架上下层间距); %网架上下层间距网架上下层间距定义网架上层节点定义网架上层节点定义网架下层节点定义网架下层节点定义网架高度定义网架高度for i=1:hu for j=1:lu Node(i-1)*lu+j,2)=(j-1)*dis_xu; Node(i-1)*lu+j,3)=(i-1)*dis_yu; Node(i-1)*lu+j,4)=dis_z; endendfor i=1:hd for j=1:ld Node(i-1)*ld+j+hu*lu,2)=(j-1+0.5)*dis_xd; Node(i-1)*ld+j+hu*lu,3)

19、=(i-1+0.5)*dis_yd; Node(i-1)*ld+j+hu*lu,4)=0; endend网架上层节点编号与对应坐标网架上层节点编号与对应坐标网架下层节点编号与对应坐标网架下层节点编号与对应坐标Nsum=length(Node); %查询网架的节点数查询网架的节点数for i=1:Nsum %将节点编号录入节点矩阵将节点编号录入节点矩阵 Node(i,1)=i;endfor i=1:hu for j=1:lu-1 Men(i-1)*(lu-1)+j,2)=(i-1)*lu+j; Men(i-1)*(lu-1)+j,3)=(i-1)*lu+j+1; endendfor i=1:lu

20、 for j=1:hu-1 Men(i-1)*(hu-1)+j+(lu-1)*hu,2)=(j-1)*lu+i; Men(i-1)*(hu-1)+j+(lu-1)*hu,3)=j*lu+i; endend节点编号的录入节点编号的录入网架上层横向单元的拓扑网架上层横向单元的拓扑网架上层纵向单元的拓扑网架上层纵向单元的拓扑 for i=1:hd for j=1:ld-1 Men(i-1)*(ld-1)+(lu-1)*hu+(hu-1)*lu+j,2)=(i-1)*ld+j+hu*lu; Men(i-1)*(ld-1)+(lu-1)*hu+(hu-1)*lu+j,3)=(i-1)*ld+j+hu*l

21、u+1; endendfor i=1:ld for j=1:hd-1 Men(i-1)*(hd-1)+(ld-1)*hd+(lu-1)*hu+(hu-1)*lu+j,2)=(j-1)*ld+i+hu*lu; Men(i-1)*(hd-1)+(ld-1)*hd+(lu-1)*hu+(hu-1)*lu+j,3)=j*ld+i+hu*lu; endend网架下层纵向单元的拓扑网架下层纵向单元的拓扑网架下层横向单元的拓扑网架下层横向单元的拓扑网架腹杆单元的拓扑网架腹杆单元的拓扑for i=1:hd for j=1:ld Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(h

22、d-1)*ld+(ld-1)*hd+1,2)=(i-1)*lu+j; Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+1,3)=(i-1)*ld+hu*lu+j;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+2,2)=(i-1)*lu+j+1;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+2,3)=(i-1)*ld+j+hu*lu;Men(i-1)*ld+j-1)*4+(hu-1)*

23、lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+3,2)=i*lu+j;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+3,3)=(i-1)*ld+j+hu*lu;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+4,2)=i*lu+j+1;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+4,3)=(i-1)*ld+j+hu*lu; endend腹杆腹杆N腹杆腹杆N+1腹

24、杆腹杆N+2腹杆腹杆N+3单元编号录入单元储存矩阵单元编号录入单元储存矩阵 Msum=length(Men); %查询网架单元数 for i=1:Msum %将单元编号录入单元将单元编号录入单元储存储存矩阵矩阵 Men(i,1)=i;end定义截面属性定义截面属性E=2.1e11;%默认材料为steelA1=input(请输入网架上层单元的截面面积); %默认网架上层单元截面尺寸相同默认网架上层单元截面尺寸相同A2=input(请输入网架下层单元的截面面积); %默认网架下层单元截面尺寸相同默认网架下层单元截面尺寸相同A3=input(请输入网架腹杆单元的截面面积); %默认网架腹杆单元截面尺

25、寸相同默认网架腹杆单元截面尺寸相同AE=zeros(Msum,3); %定义单元属性矩阵定义单元属性矩阵m1=(hu-1)*lu+(lu-1)*hu; %上层单元截止编号上层单元截止编号m2=m1+(hd-1)*ld+(ld-1)*hd; %下层单元截止编号下层单元截止编号AE(1:m1,2)=A1; %将上层单元尺寸录入将上层单元尺寸录入AE矩阵矩阵AE(m1+1):m2,2)=A2; %将下层单元尺寸录入将下层单元尺寸录入AE矩阵矩阵AE(m2+1):Msum,2)=A3; %将腹杆单元尺寸录入将腹杆单元尺寸录入AE矩阵矩阵AE(:,1)=1:Msum; %将单元编号录入将单元编号录入AE

26、矩阵矩阵AE(:,3)=E; %将材料弹性模量录入将材料弹性模量录入AE矩阵矩阵定义荷载定义荷载cont2=input(定义节点荷载，若网架上层节点力与下层节点力均布，则输入0，否则输入1);if cont2=0 P1=input(请输入网架上层节点荷载); P2=input(请输入网架下层节点荷载); m3=hu*lu; P(1:Nsum,1)=1:Nsum; P(1:m3,2)=P1;P(m3+1):Nsum,2)=P2;else P=input(定义节点荷载，按node1 P1;node2 P2;.输入); end定义边界条件定义边界条件cont3=input(定义边界约束，若网架上层周

27、边节点全约束，则输入0，若下层周边节点全约束，输入1，否则输入2);if cont3=0 n1=2*(hu+lu-2); BC=zeros(n1,4); BC(1:lu-2,1)=2:lu-1; BC(lu-1):(2*lu-4),1)=lu*(hu-1)+2:lu*hu-1; BC(2*lu-3):(2*lu-4+hu),1)=1:lu:lu*(hu-1)+1; BC(2*lu-3+hu):n1,1)=lu:lu:hu*lu;elseif cont3=1 n1=2*(hd+ld-2); BC=zeros(n1,4); BC(1:ld-2,1)=2:ld-1; BC(ld-1):(2*ld-4

28、),1)=ld*(hd-1)+2:ld*hd-1; BC(2*ld-3):(2*ld-4+hd),1)=1:ld:ld*(hd-1)+1; BC(2*ld-3+hd):n1,1)=ld:ld:hd*ld; for i=1:n1 BC(i,1)=BC(i,1)+hu*lu; endelse BC=input(定义边界约束，按node1 Conx Cony Conz;node2 Conx Cony Conz);.输入,Con代表x、y、z方向约束，取0为约束，取1无约束);endendNsum=length(Node);Msum=length(Men);Psum=length(P);BCsum=l

29、ength(BC); %提取各矩阵的行数提取各矩阵的行数考虑几何非线性分析网架考虑几何非线性分析网架for i=1:Psum %将力分为将力分为1000份份 P(i,2)=P(i,2)/1000;endU=zeros(3*Nsum,1); %总位移矩阵总位移矩阵for i=1:1000 u,L1,Kz = grid(Node,Men,AE,P,BC,Nsum,Msum,Psum,BCsum); for j=1:Nsum %根据节点位移修改网架的节点坐标根据节点位移修改网架的节点坐标 Node(j,2)=Node(j,2)+u(3*j-2,1); Node(j,3)=Node(j,3)+u(3*

30、j-1,1); Node(j,4)=Node(j,4)+u(3*j,1); end U=U+u; %每次迭代位移的叠加每次迭代位移的叠加end迭代法修正刚度矩阵和网架位移迭代法修正刚度矩阵和网架位移求解网架杆件的应力应变求解网架杆件的应力应变L0=zeros(Msum,1); %所有根杆的最初长度所有根杆的最初长度for i=1:Msum %单元两端的节点编号单元两端的节点编号 p=Men(i,2); q=Men(i,3); X1=Node(p,2); %单元端节点的坐标单元端节点的坐标 Y1=Node(p,3); Z1=Node(p,4); X2=Node(q,2); Y2=Node(q,3

31、); Z2=Node(q,4); L0(i,1)=sqrt(X2-X1)2+(Y2-Y1)2+(Z2-Z1)2); %网架杆件的初始长度网架杆件的初始长度endLt=L1-L0; %所有杆件长度的增量所有杆件长度的增量strain=zeros(Msum,1); %定义应变矩阵定义应变矩阵stress=zeros(Msum,1); %定义应力矩阵定义应力矩阵for i=1:Msum E=AE(i,3); strain(i,1)=Lt(i,1)/L0(i,1); %第第i根杆件应变根杆件应变 stress(i,1)=E*strain(i,1); %第第i根杆件应力根杆件应力enddisp(U);

32、%输出网架节点位移输出网架节点位移disp(stress); %输出网架杆件应力输出网架杆件应力网架节点位移求解函数网架节点位移求解函数function u,L1,Kz = grid(Node,Men,AE,P,BC,Nsum,Msum,Psum,BCsum)Kz=zeros(3*Nsum,3*Nsum); %定义刚度矩阵定义刚度矩阵L=zeros(Msum,1);for i=1:Msum %单元两端的节点编号单元两端的节点编号 p=Men(i,2); q=Men(i,3); A=AE(i,2); E=AE(i,3);% 单元两端节点坐标单元两端节点坐标 X1=Node(p,2);Y1=Node(p,3);Z1=Node(p,4);X2=Node(q,2);Y2=Node(q,3);Z2=Node(q,4); %单元长度单元长度 L(i,1)=sqrt(X2-X1)2+(Y2- Y1)2+(Z2-Z1)2); %单元的方向余弦单元的方向余弦 l=(X2-X1)/L(i,1); m=(Y2-Y1)/L(i,1); n=(Z2-Z1)/L(i,1); %定义定义转化矩阵转化矩阵 t=sqrt(l2+n2); if t=0 r=0 m 0;-m 0 0;0 0 1; else r=l m n;-l*m/t t -m*n/t;-n/t 0 l