规则外形的平面应力问题有限元分析.doc

规则外形的平面应力问题有限元分析 程序说明本程序适用于规则外形受集中载荷作用的平面应力问题有限元分析，采用三节点三角形单元，不包括网格划分，因此前处理需要自行处理数据；根据最小位能原理进行求解，求解结果为节点位移，不包括后置处理。程序应用举例如下： 程序举例如图为矩形简支板，板厚为t=1m，板长为18m，板高为3m，受图示作用力，弹性模量，泊松比，容重。求位移？ 程序求解过程 划分单元划分单元，标出单元号码及节点，选取坐标。 输入数据 弹性模量：eo=2000000000；泊松比：co=0.167；梁的厚度：t=1m；节点总数：nnd=14；单元总数：nne=12；节点坐标：xynnd=0 0;3 0;6 0;9 0;12 0;15 0;18 0;0 3;3 3;6 3;9 3;12 3;15 3;18 3；单元节点编码：nonne=1 2 9;2 3 10;3 4 11;4 5 12;5 6 13;6 7 14;1 9 8;2 10 9;3 11 10;4 12 11;5 13 12;6 14 13；存在载荷的节点总数；w=7；已知位移的节点总数：d=2；载荷值：nwz=0 -200000 8;0 -200000 9;0 -200000 10;0 -200000 11;0 -200000 12;0 -200000 13;0 -200000 14；已知位移值：ndz=0 0 1;0 0 7。 程序框图及程序程序框图源程序代码function UV=yxymainprogram()%有限元主程序，适用于规则外形受集中载荷作用的平面应力问题有限元分析。%弹性模量eo%泊松比co%梁的厚度t%节点总数nnd%单元总数nne%节点坐标xynnd%单元节点编码nonne%存在载荷的节点总数w%已知位移的节点总数d%载荷值nwz%已知位移值ndz%载荷矩阵P%引入位移边界后的刚度矩阵KZ%引入位移边界后的载荷矩阵KP%总纲矩阵K%单元刚度矩阵KE%应变矩阵B%弹性矩阵Dglobal UV;global KZ;global KP;global P;global w;global nnd;global nwz;global K;global B;global D;global A;global KE;global xynnd;global nonne;global eo;global co;global d;global ndz;clear;clc;KZ,KP=getKZP;UV=inv(KZ)*KP; function KZ,KP=getKZP(d)%引入位移边界后的矩阵global UV;global KZ;global KP;global P;global w;global nnd;global nwz;global K;global B;global D;global A;global KE;global xynnd;global nonne;global eo;global co;global d;global ndz;d=input(d=);w=input(w=);nwz=input(nwz=);nne=input(nne=);nnd=input(nnd=);nonne=input(nonne=);xynnd=input(xynnd=);eo=input(eo=);co=input(co=);ndz=input(ndz=);K=getK(nne,nnd);P=getP(w);a=inf;for i=1:d ii=ndz(i,3); if ndz(i,1)=0 K(:,2*ii-1)=0; K(2*ii-1,:)=0; K(2*ii-1,2*ii-1)=1; P(2*ii-1,1)=0; else K(2*ii-1,2*ii-1)=a* K(2*ii-1,2*ii-1); P(2*ii-1,1)=a* K(2*ii-1,2*ii-1)*ndz(i,1); end if ndz(i,2)=0 K(:,2*ii)=0; K(2*ii,:)=0; K(2*ii,2*ii)=1; P(2*ii,1)=0; else K(2*ii,2*ii)=a* K(2*ii,2*ii); P(2*ii,1)=a* K(2*ii,2*ii)*ndz(i,2); end KZ=K; KP=P; UV=zeros(2*nnd,1);end function K=getK(nne,nnd)%建立总纲矩阵global K;global B;global D;global A;global KE;global xynnd;global nonne;global eo;global co;K=zeros(2*nnd,2*nnd);for ne=1:nne G=zeros(6,2*nnd); i=nonne(ne,1); j=nonne(ne,2); k=nonne(ne,3); G(1,2*i-1)=1; G(2,2*i)=1; G(3,2*j-1)=1; G(4,2*j)=1; G(5,2*k-1)=1; G(6,2*k)=1; KE=getKE(ne); k=G*KE*G; K=K+k;end function KE=getKE(ne)%建立单元刚度矩阵global KE;global B;global D;global A;global eo;global co;B,A=getB(ne);D=getD(eo,co);S=D*B;KE=A*B*S; function B,A=getB(ne)% 建立应变矩阵global xynnd;global nonne;global B;global A;B=zeros(3,6);i=nonne(ne,1);j=nonne(ne,2);k=nonne(ne,3);xi=xynnd(i,1);yi=xynnd(i,2);xj=xynnd(j,1);yj=xynnd(j,2);xk=xynnd(k,1);yk=xynnd(k,2);ai=xj*yk-xk*yj;aj=xk*yi-xi*yk;ak=xi*yj-xj*yi;bi=yj-yk;bj=yk-yi;bk=yi-yj;gi=xk-xj;gj=xi-xk;gk=xj-xi;h=(ai+aj+ak);B(1,1)=bi/h;B(1,3)=bj/h;B(1,5)=bk/h;B(2,2)=gi/h;B(2,4)=gj/h;B(2,6)=gk/h;B(3,1)=gi/h;B(3,2)=bi/h;B(3,3)=gj/h;B(3,4)=bj/h;B(3,5)=gk/h;B(3,6)=bk/h;A=0.5*abs(ai+aj+ak); function D=getD(eo,co)%建立弹性矩阵global D;D=zeros(3,3);d1=eo/(1-co*co);d2=co;d3=0.5*d1*(1-d2);D(1,1)=d1;D(1,2)=d1*d2;D(2,1)=D(1,2);D(2,2)=d1;D(3,3)=d3; function P=getP(w)%建立载荷矩阵global P;global w;global nnd;global nwz;P=zeros(2*nnd,1);for i=1:w ii=nwz(i,3); P(2*ii-1,1)=nwz(i,1); P(2*ii,1)=nwz(i,2);end 程序求解结果ans =00-0