有限元的matlab编程

有限元编程示例题目描述:如下图所示的平面桁架，杆件长度、弹性模量、截面积以及所受节点力P的大小可以自行定义。求节点位移及杆件轴力。例一：桁架解题思路：•建立模型•集成总刚•求解位移•求解杆件轴力•输出结果建立模型：定义节点坐标Node=zeros(10,2);x=-1*L;%L为横杆长度fori=1:2:10x=x+L;Node(i,:)=[x0];endx=-1*L;fori=2:2:10x=x+L;Node(i,:)=[xH];%H为竖杆长度end模型相关参数输入H=input('竖杆长度(m):');L=input('水平杆长度(m):');E=input('杆件弹性模量(Gpa):');A=input('杆件截面积(m^2):');a=input('节点力P(kN)：');节点编号方式定义单元，即储存单元两端的节点号Element=zeros(21,2);fori=1:2:7Element(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];endfori=2:2:8Element(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);获取单元两端节点坐标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为之前输入的节点力计算从局部坐标到整体坐标的坐标转换矩阵TfunctionT=TransformMatrix(ie)%ie为单元号c=(xj-xi)/L;s=(yj-yi)/L;T=[c-s00sc0000c-s00sc];计算单元刚度矩阵kk=[E*A/L0-E*A/L00000-E*A/L0E*A/L00000];T=TransformMatrix(ie);k=T*k*transpose(T);%transpose(T)为T的转置矩阵2集成整体刚度矩阵Kforie=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集成总刚的对角线元素（这里的元素指2*2的小矩阵）fori=1:1:10%按节点的顺序循环forj=1:1:21%对于每个节点，再按单元的顺序循环k=PlaneTrussElementStiffness(j);ifElement(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);endifElement(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);endendend求解位移：u=zeros(20);根据约束情况修改总刚，采用对角元素置1法fori=1:1:20K(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=K\f;求解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=TransformMatrix(ie);%计算坐标转换矩阵nodef=transpose(T)*nodef;%从整体坐标转换到局部坐标计算单元的节点力输出求解结果：输出位移fprintf('节点位移\n');fori=1:1:10disp(['节点号',num2str(i),'x方向位移：',num2str(u(2*i-1,1)),'y方向位移：',num2str(u(2*i,1))]);end输出节点力fprintf('\n\n节点力\n');forie=1:1:21nodef=NodeForce(ie);disp(['单元号:',num2str(ie),'节点号：',num2str(Element(ie,1)),'节点号：',num2str(Element(ie,2)),'轴力:',num2str(nodef(1))]);end例二：网架思路分析网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成网架的基本单元有三角锥，三棱体，正方体，截头四角锥等。鉴于网架的形式较多，本程序提供一种通用的网架输入方法，但录入较为繁琐，同时提供一种正放四角锥网架的简易输入方法作为典型。考虑几何非线性。本程序采用荷载分级加载的方式考虑网架的几何非线性。将总荷载分成1000份分步施加，求解各荷载步下的节点位移，修改网架相应节点坐标以及刚度矩阵，依次迭代求出网架的总位移。本程序的网架位移求解函数附在主程序后面，主程序运行时调用该函数。几点说明用户自定义输入几何建模正放四角锥网架简易输入定义荷载定义边界条件网架分析位移应变应力位移求解函数单刚矩阵荷载矩阵约束矩阵总刚矩阵求解位移&&&分级加载，通过修改节点坐标，迭代求解几何非线性e=input('选择网架类型，0代表自由定义网架，1代表四角锥网架')%网架类型的选择网架类型的选择用户自定义网架（网架信息的录入，包括节点、单元、截面、弹性模量等）ife==0%选择自定义网架Node=input(‘定义节点编号及对应坐标，按[1x1y1z1;2x2y2z2;...]输入’);%形成节点储存矩阵Men=input(‘定义单元与节点的关系，按[1node1node2;2node3node4;...]输入,node1node2,依次类推’);%形成单元储存矩阵Msum=length(Men);%查找网架录入的单元数Cont1=input('定义单元实常数，若所有杆件截面面积和弹性模量不变，则输入0，否则输入1');定义单元属性的输入方式ifCont1==0AE1=input('请输入统一的截面面积与弹性模量，按[AE]输入');AE=zeros(Msum,3);AE(:,1)=1:Msum;AE(:,2)=AE1(1,1);AE(:,3)=AE1(1,2);elseAE=input('请输入相应单元的截面面积与弹性模量，按[1,A1E1;2,A2E2;...]输入');endP=input(‘定义节点荷载，按[node1P1;node2P2;...]输入’);%网架荷载输入BC=input(‘定义边界约束，按[node1ConxConyConz;node2ConxConyConz);...]输入,Con代表x、y、z方向约束，取0为约束，取1无约束’);%网架边界条件end单元属性相同单元属性不同荷载及边界条件正放四角锥网架定义ife==1hu=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('输入网架上下层间距');%网架上下层间距定义网架上层节点定义网架下层节点定义网架高度fori=1:huforj=1:luNode((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;endendfori=1:hdforj=1:ldNode((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;endend网架上层节点编号与对应坐标网架下层节点编号与对应坐标Nsum=length(Node);%查询网架的节点数fori=1:Nsum%将节点编号录入节点矩阵Node(i,1)=i;endfori=1:huforj=1:lu-1Men((i-1)*(lu-1)+j,2)=(i-1)*lu+j;Men((i-1)*(lu-1)+j,3)=(i-1)*lu+j+1;endendfori=1:luforj=1:hu-1Men((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节点编号的录入网架上层横向单元的拓扑网架上层纵向单元的拓扑fori=1:hdforj=1:ld-1Men((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;endendfori=1:ldforj=1:hd-1Men((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网架下层纵向单元的拓扑网架下层横向单元的拓扑网架腹杆单元的拓扑fori=1:hdforj=1:ldMen(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*h