节点应力和单元应力的关系

一：用三结点三角形平面单元计算平面结构的应力和位移。1，说明计算简图，网格划分，单元及结点的编号如下图所示。由于结构对称，去四分之一结构分析。其中E=2e10pa,mu=0.167,h=1m.变量注释：Node-------节点定义gElement----单元定义gMaterial---材料定义，包括弹性模量，泊松比和厚度gBC1--------约束条件gNF---------集中力gk------------总刚gDelta-------结点位移子程序注释：PlaneStructualModel———定义有限元模型SolveModel———————求解有限元模型DisplayResults——————显示计算结果k=StiffnessMatrix(ie)———计算单元刚度AssembleStiffnessMatrix(ie,k)—形成总刚es=ElementStress(ie)————计算单元应力functionexam1%输入参数：无%输出结果：节点位移和单元应力PlaneStructualModel;%定义有限元模型SolveModel;%求解有限元模型DisplayResults;%显示计算结果return;functionPlaneStructualModel%定义平面结构的有限元模型%输入参数：无%说明：%该函数定义平面结构的有限元模型数据：%gNode-------节点定义%gElement----单元定义%gMaterial---材料定义，包括弹性模量，泊松比和厚度%gBC1--------约束条件%gNF---------集中力globalgNodegElementgMaterialgBC1gNF%节点坐标%xygNode=[0.0,2.0%节点10.0,1.0%节点21.0,1.0%节点30.0,0.0%节点41.0,0.0%节点52.0,0.0];%节点6%单元定义%节点1节点2节点3材料号gElement=[3,1,2,1%单元15,2,4,1%单元22,5,3,1%单元36,3,5,1];%单元4%材料性质%弹性模量泊松比厚度gMaterial=[1e0,0,1];%材料1%第一类约束条件%节点号自由度号约束值gBC1=[1,1,0.02,1,0.04,1,0.04,2,0.05,2,0.06,2,0.0];%集中力%节点号自由度号集中力值gNF=[1,2,-1];returnfunctionSolveModel%求解有限元模型%输入参数：无%说明：%该函数求解有限元模型，过程如下%1.计算单元刚度矩阵，集成整体刚度矩阵%2.计算单元的等效节点力，集成整体节点力向量%3.处理约束条件，修改整体刚度矩阵和节点力向量%4.求解方程组，得到整体节点位移向量globalgNodegElementgMaterialgBC1gNFgKgDelta%step1.定义整体刚度矩阵和节点力向量[node_number,dummy]=size(gNode);gK=sparse(node_number*2,node_number*2);f=sparse(node_number*2,1);%step2.计算单元刚度矩阵，并集成到整体刚度矩阵中[element_number,dummy]=size(gElement);forie=1:1:element_numberk=StiffnessMatrix(ie);AssembleStiffnessMatrix(ie,k);end%step3.把集中力直接集成到整体节点力向量中[nf_number,dummy]=size(gNF);forinf=1:1:nf_numbern=gNF(inf,1);d=gNF(inf,2);f((n-1)*2+d)=gNF(inf,3);end%step4.处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法[bc_number,dummy]=size(gBC1);foribc=1:1:bc_numbern=gBC1(ibc,1);d=gBC1(ibc,2);m=(n-1)*2+d;f(m)=gBC1(ibc,3)*gK(m,m)*1e15;gK(m,m)=gK(m,m)*1e15;end%step5.求解方程组，得到节点位移向量gDelta=gK\f;returnfunctionDisplayResults%显示计算结果%输入参数：无globalgNodegElementgMaterialgBC1gNFgKgDeltafprintf('节点位移\n');fprintf('节点号x方向位移y方向位移\n');[node_number,dummy]=size(gNode);fori=1:node_numberfprintf('%6d%16.8e%16.8e\n',...i,gDelta((i-1)*2+1),gDelta((i-1)*2+2));endfprintf('\n\n单元应力\n');fprintf('X-STRY-STRXY-STR\n');[element_number,dummy]=size(gElement);forie=1:element_numberes=ElementStress(ie);fprintf('单元号%6d%16.8e%16.8e%16.8e\n',...ie,es(1),es(2),es(3));endreturnfunctionk=StiffnessMatrix(ie)%计算单元刚度矩阵%输入参数:%ie-------单元号globalgNodegElementgMaterialk=zeros(6,6);E=gMaterial(gElement(ie,4),1);mu=gMaterial(gElement(ie,4),2);h=gMaterial(gElement(ie,4),3);xi=gNode(gElement(ie,1),1);yi=gNode(gElement(ie,1),2);xj=gNode(gElement(ie,2),1);yj=gNode(gElement(ie,2),2);xm=gNode(gElement(ie,3),1);ym=gNode(gElement(ie,3),2);ai=xj*ym-xm*yj;aj=xm*yi-xi*ym;am=xi*yj-xj*yi;bi=yj-ym;bj=ym-yi;bm=yi-yj;ci=-(xj-xm);cj=-(xm-xi);cm=-(xi-xj);A=(ai+aj+am)/2;B=[bi0bj0bm00ci0cj0cmcibicjbjcmbm];B=B/2/A;D=[1mu0mu1000(1-mu)/2];D=D*E/(1-mu^2);k=transpose(B)*D*B*h*abs(A);returnfunctionAssembleStiffnessMatrix(ie,k)%把单元刚度矩阵集成到整体刚度矩阵%输入参数:%ie---单元号%k---单元刚度矩阵globalgElementgKfori=1:1:3forj=1:1:3forp=1:1:2forq=1:1:2m=(i-1)*2+p;n=(j-1)*2+q;M=(gElement(ie,i)-1)*2+p;N=(gElement(ie,j)-1)*2+q;gK(M,N)=gK(M,N)+k(m,n);endendendendreturnfunctiones=ElementStress(ie)%计算单元的应力%输入参数%ie-----节点号%es-----单元应力globalgElementgNodegDeltagMateriales=zeros(1,6);de=zeros(6,1);forj=1:1:3de(2*j-1)=gDelta(2*gElement(ie,j)-1);de(2*j)=gDelta(2*gElement(ie,j));endE=gMaterial(gElement(ie,4),1);mu=gMaterial(gElement(ie,4),2);h=gMaterial(gElement(ie,4),3);xi=gNode(gElement(ie,1),1);yi=gNode(gElement(ie,1),2);xj=gNode(gElement(ie,2),1);yj=gNode(gElement(ie,2),2);xm=gNode(gElement(ie,3),1);ym=gNode(gElement(ie,3),2);ai=xj*ym-xm*yj;aj=xm*yi-xi*ym;am=xi*yj-xj*yi;bi=yj-ym;bj=ym-yi;bm=yi-yj;ci=-(xj-xm);cj=-(xm-xi);cm=-(xi-xj);A=(ai+aj+am)/2;B=[bi0bj0bm00ci0cj0cmcibicjbjcmbm];B=B/2/A;D=[1mu0mu1000(1-mu)/2];D=D*E/(1-mu^2);S=D*B;es(1:3)=S*de;es(6)=0.5*sqrt((es(1)-es(2))^2+4*es(3)^2);es(4)=0.5*(es(1)+es(2))+es(6);es(5)=0.5*(es(1)+es(2))-es(6);return3，数据文件：输入数据：gNode=[0.0,2.00.0,1.01.0,1.00.0,0.01.0,0.02.0,0.0];gElement=[3,1,2,15,2,4,12,5,3,16,3,5,1];gMaterial=[1e0,0,1];gBC1=[1,1,0.02,1,0.04,1,0.04,2,0.05,2,0.06,2,0.0];gNF=[1,2,-1];输出数据：节点位移节点号x方向位移y方向位移1-8.79120879e-016-3.25274725e+00028.79120879e-017-1.25274725e+0003-8.79120879e-002-3.73626374e-00141.17216117e-016-8.35164835e-01651.75824176e-001-2.93040293e-01661.75824176e-0012.63736264e-016单元应力X-STRY-STRXY-STR单元号1-8.79120879e-002-2.00000000e+0004.39560440e-001单元号21.75824176e-001-1.25274725e+0002.56410256e-016单元号3-8.79120879e-002-3.73626374e-0013.07692308e-001单元号40.00000000e+000-3.73626374e-001-1.31868132e-001（注：可编辑下载，若有不当之处，请指正，谢谢!）开始输入节点坐标进行单元定义给出材料性质定义有限元模型引入约束条件给定荷载信息计算单刚形成总刚形成结点力求解有限元模型引进位移约束求解平衡方程显示结点位移显示计算结果显示单元应力结束