《有限元的matlab编程.ppt》由会员分享，可在线阅读，更多相关《有限元的matlab编程.ppt（36页珍藏版）》请在三一办公上搜索。
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;for i=2:2:10 x=x+L;Node(i,:)=x H;%H为竖杆长度end,模型相关参数输入,H=input(竖杆长度(m):);L=input(水平杆长度(m):);E=input(杆件弹性模量(Gpa):
2、);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;Element(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
3、,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为之前输入的节点力,计算从局部坐标到整体坐标的坐标转换矩阵T,function T=TransformMatrix(ie)%ie为单元号,c=(xj-xi)/L;s=(yj-yi)/L;T=c-
4、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,集成整体刚度矩阵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-
5、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,集成总刚
6、的对角线元素（这里的元素指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单
7、元的末节点 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 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被约束了,所在的行和列的其他元素都改为0,K=K*1e15
8、;%乘以一个大数，减小计算误差f=f*1e15;u=Kf;,求解,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);nodef=k*uElement;%整体坐标下的节点力T=Trans
9、formMatrix(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,输出节点力,fprintf(nn节点力n);for ie=1:1:21 nodef=NodeForce(ie);disp(单元号:,num2str(ie),节点号：,num2str(Element(ie,1),节点号：,num
10、2str(Element(ie,2),轴力:,num2str(nodef(1);end,例二：网架,思路分析,网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成网架的基本单元有三角锥，三棱体，正方体，截头四角锥等。鉴于网架的形式较多，本程序提供一种通用的网架输入方法，但录入较为繁琐，同时提供一种正放四角锥网架的简易输入方法作为典型。,考虑几何非线性。本程序采用荷载分级加载的方式考虑网架的几何非线性。将总荷载分成1000份分步施加，求解各荷载步下的节点位移，修改网架相应节点坐标以及刚度矩阵，依次迭代求出网架的总位移。,本程序的网架位移求解函数附在主程序后面，主程序运行时调用该函数
11、。,几点说明,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(定义单元实常数，若所
12、有杆件截面面积和弹性模量不变，则输入0，否则输入1);,定义单元属性的输入方式,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
13、 Conx Cony Conz);.输入,Con代表x、y、z方向约束，取0为约束，取1无约束);%网架边界条件,end,单元属性相同,单元属性不同,荷载及边界条件,正放四角锥网架定义,if e=1,hu=input(输入网架上层节点行数);%定义网架上层节点的行数lu=input(输入网架上层节点列数);%定义网架上层节点的列数,dis_xu=input(输入网架上层节点列间距);%定义网架上层的行间距dis_yu=input(输入网架上层节点行间距);%定义网架上层的列间距,hd=hu-1;%网架下层节点的行数ld=lu-1;%网架下层节点的列数,dis_xd=dis_xu;%网架下层的行
14、间距dis_yd=dis_yu;%网架下层的行间距,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;endend,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_y
15、d;Node(i-1)*ld+j+hu*lu,4)=0;endend,网架上层节点编号与对应坐标,网架下层节点编号与对应坐标,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;endend,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
16、-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*lu+1;endend,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
17、)*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+(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
18、(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+
19、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,腹杆N+2,腹杆N+3,单元编号录入单元储存矩阵,Msum=length(Men);%查询网架单元数,for i=1:Msum%将单元编号录入单元储存矩阵 Men(i,1)=i;end,定义截面属性,E=2.1e11;%默认材料为stee
20、l,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;%下层单元截止编号,AE(1:m1,2)=A1;%将上层单元尺寸录入AE矩阵AE(m1+1):m2,2)=A2;%将下层单元尺寸录入AE矩阵AE(m2+1):Msum,2)=A3;%将
21、腹杆单元尺寸录入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,定义边界条件,cont3=inp
22、ut(定义边界约束，若网架上层周边节点全约束，则输入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-
23、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无约束);endend,Nsum=length(Node);Msum=length(Men);Psum=length(P)
24、;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
25、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);Z2=Node(q,4);,L0(i,1)=sqrt(X2-X1)2+(Y2-Y1)2+(Z2-Z1)2);%网架杆件的初始长度end,Lt=L1-L0;%所有杆件长度的增量,strain=zeros(Msu
26、m,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);%输出网架杆件应力,网架节点位移求解函数,function u,L1,Kz=grid(Node,Men,AE,P,BC,Nsum,Msum,Psum,BCsum),Kz=zeros(3*Nsum,3*Nsum);%定义刚度矩阵L=zeros(Msum,1);,f
27、or 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
28、 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;,%单刚矩阵的扩充，使之行数、列数与总刚对应 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
29、;%形成总刚矩阵,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;endend,u=KzRz;%节点位移L1=L;%单元长度矩阵,end,求解网架位移,具体算例,选择网架类型，0代表自由定义网架，1代表四角锥网架：1输入网架上层节点行数：4输入网架上层节点列数：6输入网架上层节点列间距：2输入网架上层节点行间距：2输入网架上下
30、层间距：2请输入网架上层单元的截面面积：0.0003请输入网架下层单元的截面面积：0.0003请输入网架腹杆单元的截面面积：0.0004定义节点荷载，若网架上层节点力与下层节点力均布，则输入0，否则输入1：0请输入网架上层节点荷载：10000请输入网架下层节点荷载：0定义边界约束，若网架上层周边节点全约束，则输入0，若下层周边节点全约束，输入1，否则输入2：0,截面输入部分,作业：,编制框架结构在外荷载作用下结构的应力、位移等响应。Level 1（60%）：平面框架+静力荷载Level 2（70%）：平面框架+动力荷载Level 3（80%）：空间框架+静力荷载Level 4（90%）：空间框架+动力荷载Level 5（100%）：空间框架+静力荷载+动力荷载,提交作业要求：,程序源文件报告（含计算算例，计算结果的截图等）请于2015年6月24日前将所有作业发至邮箱(这个邮箱只用于收作业）有其他问题请发邮件至：caijg_,
