有限元的matlab编程.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,有限元编程示例,1,题目描述:,如下图所示的平面桁架，杆件长度、弹性模量、截面积以及所受节点力P的大小可以自行定义。求节点位移及杆件轴力。,例一：桁架,2,解题思路：,建立模型,集成总刚,求解位移,求解杆件轴力,输出结果,3,建立模型：,定义节点坐标,Node=zeros(10,2);,x=-1*L;,%L,为横杆长度,for i=1:2:10,x=x+L;,Node(i,:)=x 0;,end,x=-1*L;,for 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)：);,节点编号方式,4,定义单元，即储存单元两端的节点号,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;,Element(5/2*i+1/2,:)=i,i+3;,end,for i=2:2:8,Element(5*i/2-1,:)=i,i+1;,Element(5*i/2,:)=i,i+2;,end,Element(21,:)=9,10;,加下划线的为单元编号,5,集成总刚：,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);,获取单元两端节点坐标,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,为之前输入的节点力,6,计算从局部坐标到整体坐标的坐标转换矩阵,T,function T=TransformMatrix(ie),%ie,为单元号,c=(xj-xi)/L;,s=(yj-yi)/L;,T=c -s 0 0,s c 0 0,0 0 c -s,0 0 s c ;,计算单元刚度矩阵,k,k=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,7,集成整体刚度矩阵,K,for ie=1:1:21,%,按单元顺序进行循环,k=PlaneTrussElementStiffness(ie);,%,计算第,ie,个单元的单刚,m=Element(ie,1);,%ie,单元的首节点号,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,单元的末节点号,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,8,集成总刚的对角线元素（这里的元素指,2,*,2,的小矩阵）,for i=1:1:10,%,按节点的顺序循环,for j=1:1:21,%,对于每个节点，再按单元的顺序循环,k=PlaneTrussElementStiffness(j);,if Element(j,1)=,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*i,2*i-1)+k(4,3);,K(2*i,2*i)=K(2*i,2*i)+k(4,4);,end,end,end,9,求解位移：,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,被约束了,所在的行和列的其他元素都改为,0,K=K*1e15;,%,乘以一个大数，减小计算误差,f=f*1e15;,u=Kf;,求解,K(1,1)=1;,%,对角线元素置,1,K(2,2)=1;,K(18,18)=1;,10,求解轴力：,获取单元两端的节点号,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);,nodef=k*uElement;,%,整体坐标下的节点力,T=TransformMatrix(ie);,%,计算坐标转换矩阵,nodef=transpose(T)*nodef;,%,从整体坐标转换到局部坐标,计算单元的节点力,11,输出求解结果：,输出位移,fprintf(节点位移n);,for i=1:1:10,disp(节点号,num2str(i),x方向位移：,num2str(u(2*i-1,1),y方向位移：,num2str(u(2*i,1);,end,输出节点力,fprintf(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,12,例二：网架,13,思路分析,网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成网架的基本单元有三角锥，三棱体，正方体，截头四角锥等。,鉴于网架的形式较多，本程序提供一种通用的网架输入方法，但录入较为繁琐，同时提供一种正放四角锥网架的简易输入方法作为典型。,考虑几何非线性,。本程序采用荷载分级加载的方式考虑网架的几何非线性。将总荷载分成,1000,份分步施加，求解各荷载步下的节点位移，修改网架相应节点坐标以及刚度矩阵，依次迭代求出网架的总位移。,本程序的网架位移求解函数附在主程序后面，主程序运行时调用该函数。,几点说明,14,用户自定义输入,几何建模,正放四角锥网架简易输入,定义荷载,定义边界条件,网架分析,位移,应变,应力,位移求解函数,单刚矩阵,荷载矩阵,约束矩阵,总刚矩阵,求解位移,&,&,&,分级加载，通过修改节点坐标，迭代求解,几何非线性,15,e=input(,选择网架类型，,0,代表自由定义网架，,1,代表四角锥网架,),%,网架类型的选择,网架类型的选择,用户自定义网架（网架信息的录入，包括节点、单元、截面、弹性模量等）,if e=0,%,选择自定义网架,Node=input(,定义节点编号及对应坐标，按,1 x1 y1 z1;2 x2 y2 z2;.,输入,);%,形成节点储存矩阵,Men=input(,定义单元与节点的关系，按,1 node1 node2;2 node3 node4;.,输入,node1node2,依次类推,);,%,形成单元储存矩阵,Msum=length(Men);,%,查找网架录入的单元数,Cont1=input(,定义单元实常数，若所有杆件截面面积和弹性模量不变，则输入,0,，否则输入,1);,定义单元属性的输入方式,16,if Cont1=0,AE1=input(,请输入统一的截面面积与弹性模量，按,A E,输入,);,AE=zeros(Msum,3);,AE(:,1)=1:Msum;AE(:,2)=AE1(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=1,hu=input(,输入网架上层节点行数,);,%,定义网架上层节点的行数,lu=input(,输入网架上层节点列数,);,%,定义网架上层节点的列数,dis_xu=input(,输入网架上层节点列间距,);,%,定义网架上层的行间距,dis_yu=input(,输入网架上层节点行间距,);,%,定义网架上层的列间距,hd=hu-1;,%,网架下层节点的行数,ld=lu-1,;%,网架下层节点的列数,dis_xd=dis_xu;,%,网架下层的行间距,dis_yd=dis_yu;,%,网架下层的行间距,dis_z=input(,输入网架上下层间距,);,%,网架上下层间距,定义网架上层节点,定义网架下层节点,定义网架高度,18,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;,end,end,for 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)=(i-1+0.5)*dis_yd;,Node(i-1)*ld+j+hu*lu,4)=0;,end,end,网架上层节点编号与对应坐标,网架下层节点编号与对应坐标,19,Nsum=length(Node);,%,查询网架的节点数,for i=1:Nsum,%,将节点编号录入节点矩阵,Node(i,1)=i;,end,for 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;,end,end,for i=1:lu,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;,end,end,节点编号的录入,网架上层横向单元的拓扑,网架上层纵向单元的拓扑,20,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*lu+1;,end,end,for 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;,end,end,网架下层纵向单元的拓扑,网架下层横向单元的拓扑,21,网架腹杆单元的拓扑,for i=1:hd,for j=1:ld,Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-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)*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;,end,end,腹杆,N,腹杆,N+1,腹杆,N+2,腹杆,N+3,22,单元编号录入单元储存矩阵,Msum=length(Men);%,查询网架单元数,for i=1:Msum,%,将单元编号录入单元,储存,矩阵,Men(i,1)=i;,end,定义截面属性,E=2.1e11;%,默认材料为,steel,A1=input(,请输入网架上层单元的截面面积,);,%,默认网架上层单元截面尺寸相同,A2=input(,请输入网架下层单元的截面面积,);,%,默认网架下层单元截面尺寸相同,A3=input(,请输入网架腹杆单元的截面面积,);,%,默认网架腹杆单元截面尺寸相同,AE=zeros(Msum,3);,%,定义单元属性矩阵,m1=(hu-1)*lu+(lu-1)*hu;,%,上层单元截止编号,m2=m1+(hd-1)*ld+(ld-1)*hd;,%,下层单元截止编号,23,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,矩阵,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,24,定义边界条件,cont3=input(,定义边界约束，若网架上层周边节点全约束，则输入,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;,25,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),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;,end,else,BC=input(,定义边界约束，按,node1 Conx Cony Conz;node2 Conx Cony Conz);.,输入,Con,代表,x,、,y,、,z,方向约束，取,0,为约束，取,1,无约束,);,end,end,26,Nsum=length(Node);Msum=length(Men);,Psum=length(P);BCsum=length(BC);,%,提取各矩阵的行数,考虑几何非线性分析网架,for i=1:Psum,%,将力分为,1000,份,P(i,2)=P(i,2)/1000;,end,U=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*j-1,1);,Node(j,4)=Node(j,4)+u(3*j,1);,end,U=U+u;,%,每次迭代位移的叠加,end,迭代法修正刚度矩阵和网架位移,27,求解网架杆件的应力应变,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);,Z2=Node(q,4);,L0(i,1)=sqrt(X2-X1)2+(Y2-Y1)2+(Z2-Z1)2);,%,网架杆件的初始长度,end,28,Lt=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,根杆件应力,end,disp(U);,%,输出网架节点位移,disp(stress);,%,输出网架杆件应力,29,网架节点位移求解函数,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);,30,%,定义,转化矩阵,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/t;,end,O=zeros(3,3);,T=r O;O r;,%,整体坐标下的单刚矩阵,ke=E*A/L(i,1)*1 0 0-1 0 0;0 0 0 0 0 0;0 0 0 0 0 0;-1 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 0 0;,k=T*ke*T;,31,%,单刚矩阵的扩充，使之行数、列数与总刚对应,G=zeros(6,3*Nsum);,I=1 0 0;0 1 0;0 0 1;,G(1:3,3*p-2:3*p)=I;,G(4:6,3*q-2:3*q)=I;,K=G*k*G;,Kz=Kz+K;,%,形成总刚矩阵,end,%,节点力矩阵的扩充,Rz=zeros(3*Nsum,1);,for i=1:Psum,j=P(i,1);,Rz(3*j,1)=P(i,2);,end,%,对角元素置,1,法处理刚度矩阵,与节点力矩阵,for i=1:3*Nsum,if BCz(i)=0,Kz(i,:)=0;,Kz(:,i)=0;,Kz(i,i)=1;,Rz(i)=0;,end,end,u=KzRz;,%,节点位移,L1=L;,%,单元长度矩阵,end,求解网架位移,32,具体算例,33,选择网架类型，,0,代表自由定义网架，,1,代表四角锥网架：,1,输入网架上层节点行数：,4,输入网架上层节点列数：,6,输入网架上层节点列间距：,2,输入网架上层节点行间距：,2,输入网架上下层间距：,2,请输入网架上层单元的截面面积：,0.0003,请输入网架下层单元的截面面积：,0.0003,请输入网架腹杆单元的截面面积：,0.0004,定义节点荷载，若网架上层节点力与下层节点力均布，则输入,0,，否则输入,1,：,0,请,输入网架上层节点荷载：,10000,请输入网架下层节点荷载：,0,定义边界约束，若网架上层周边节点全约束，则输入,0,，若下层周边节点全约束，输入,1,，否则输入,2,：,0,截面输入部分,34,作业：,编制框架结构在外荷载作用下结构的应力、位移等响应。,Level 1,（,60%,）：平面框架,+,静力荷载,Level 2,（,70%,）：平面框架,+,动力荷载,Level 3,（,80%,）：空间框架,+,静力荷载,Level 4,（,90%,）：空间框架,+,动力荷载,Level 5,（,100%,）：空间框架,+,静力荷载,+,动力荷载,35,提交作业要求：,程序源文件,报告（含计算算例,，计算结果的,截图等）,请于,2015,年,6,月,24,日前,将所有作业发至邮箱,caijg_ren,(,这个邮箱只用于收作业）,有其他问题请发邮件至：,caijg_ren,36,