规则外形的平面应力问题有限元分析规则外形的平面应力问题有限元分析 程序说明本程序适用于规则外形受集中载荷作用的平面应力问题有限元分析，采用三节点三角形单元，不包括网格划分，因此前处理需要自行处理数据；根据最小位能原理进行求解，求解结果为节点位移，不包括后置处理。程序应用举例如下： 程序举例如图为矩形简支板，板厚为 t=1m，板长为 18m，板高为 3m，受图示作用力，弹性模量，泊松比，容重。求位移？2/20mGNE 167. 00 程序求解过程 划分单元划分单元，标出单元号码及节点，选取坐标。 输入数据弹性模量：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:dii=ndz(i,3);if ndz(i,1)=0K(:,2*ii-1)=0;K(2*ii-1,:)=0;K(2*ii-1,2*ii-1)=1;P(2*ii-1,1)=0;elseK(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);endif ndz(i,2)=0K(:,2*ii)=0;K(2*ii,:)=0;K(2*ii,2*ii)=1;P(2*ii,1)=0;elseK(2*ii,2*ii)=a* K(2*ii,2*ii);P(2*ii,1)=a* K(2*ii,2*ii)*ndz(i,2);endKZ=K;KP=P;UV=zeros(2*nnd,1);endfunction 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:nneG=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;endfunction 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:wii=nwz(i,3);P(2*ii-1,1)=nwz(i,1);P(2*ii,1)=nwz(i,2);end 程序求解结果ans =00-0.0003-0.0021-0.0003-0.0036-0.0001-0.00420.0001-0.00360.0001-0.0022000.0013-0.00040.0011-0.00210.0006-0.0036-0.0001-0.0042-0.0007-0.0037-0.0012-0.0023-0.0015-0.0004