規(guī)則外形的平面應力問題有限元分析.doc

規(guī)則外形的平面應力問題有限元分析 程序說明本程序適用于規(guī)則外形受集中載荷作用的平面應力問題有限元分析，采用三節(jié)點三角形單元，不包括網(wǎng)格劃分，因此前處理需要自行處理數(shù)據(jù)；根據(jù)最小位能原理進行求解，求解結(jié)果為節(jié)點位移，不包括后置處理。程序應用舉例如下： 程序舉例如圖為矩形簡支板，板厚為t=1m，板長為18m，板高為3m，受圖示作用力，彈性模量，泊松比，容重。求位移？ 程序求解過程 劃分單元劃分單元，標出單元號碼及節(jié)點，選取坐標。 輸入數(shù)據(jù) 彈性模量：eo=2000000000；泊松比：co=0.167；梁的厚度：t=1m；節(jié)點總數(shù)：nnd=14；單元總數(shù)：nne=12；節(jié)點坐標：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；單元節(jié)點編碼：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；存在載荷的節(jié)點總數(shù)；w=7；已知位移的節(jié)點總數(shù)：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()%有限元主程序，適用于規(guī)則外形受集中載荷作用的平面應力問題有限元分析。%彈性模量eo%泊松比co%梁的厚度t%節(jié)點總數(shù)nnd%單元總數(shù)nne%節(jié)點坐標xynnd%單元節(jié)點編碼nonne%存在載荷的節(jié)點總數(shù)w%已知位移的節(jié)點總數(shù)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 程序求解結(jié)果ans =00-0