二维有限元分析的MATLAB仿真研究

登录后复制

%Xiamen University, FEM Class - Fall, 2006 %Patch test of 2D codeclear all; close all; % Geometry properties for a rectangular shape length = 10; height = 10; % number of elements in each direction ndivl = 4; ndivw = ndivl;% x,y:the coordinate of every node; node:the relationship of every node; numele;number of elements; total numnod:number of nodes [x,y,x1,y1,node,numele,numnod] = fem2d_mesh(length,height,ndivl,ndivw); y1=y1*1.9; for i=1:ndivw+1 A=y1(i,i); B=(120-(3+2)*A)/10; x1(i,:)=B*x1(i,:)+3*A; end for i = 1:(ndivl+1) for j=1:(ndivw+1) x((ndivw+1)*(i-1)+j) = x1(j,i); y((ndivw+1)*(i-1)+j) = y1(j,i); end end % Material properties% Force and Displacement BC'S [ifix,disp] = fem2d_ebcs(x,y,numnod,ndivl,ndivw);% Construct Stifffness ndof = 2; %degrees of freedom per node % Guass integration points and weights ksi(1)=-1/sqrt(3); ksi(2)=1/sqrt(3); weight(1)=1; weight(2)=1; % numequns:total number of equations; bigk:global stiffness; force:global force numeqns = numnod*ndof; bigk = zeros(numeqns); force = zeros(numeqns,1);% Loop over elements % nen is number of nodes per element nen = 4; for e = 1:numele % ke:element stiffness [ke] = fem2d_stiffness(node,x,y,ksi,weight,e); % assemble ke into bigk n1 = ndof-1; for i=1:nen; for j=1:nen; rbk = ndof*(node(i,e)-1) + 1; % row number of bigk cbk = ndof*(node(j,e)-1) + 1; % colunm number of bigk rbk1 = ndof*node(i,e); % row number of bigk cbk1 = ndof*node(j,e); % colunm number of bigk re = ndof*(i-1)+1; % row number of ke ce = ndof*(j-1)+1; % colunm number of ke re1 = ndof*i; % row number of ke ce1 = ndof*j; % colunm number of ke bigk(rbk:rbk+n1, cbk:cbk+n1) = bigk(rbk:rbk+n1, cbk:cbk+n1) + ke(re:re+n1, ce:ce+n1); end end end% Apply zero essential boundary conditions for n=1:numnod if (ifix(n) == 1) force(:) = force(:) - disp(2*n-1)*bigk(:,2*n-1); bigk(2*n-1,:) = 0; bigk(:,2*n-1) = 0; bigk(2*n-1,2*n-1) = 1.0; force(:) = force(:) - disp(2*n)*bigk(:,2*n); bigk(2*n,:) = 0; bigk(:,2*n) = 0; bigk(2*n,2*n) = 1.0; end end for n=1:numnod if (ifix(n) == 1) force(2*n) = disp(2*n); force(2*n-1) = disp(2*n-1); end end% Solve stiffness equations disp = bigk\force; % Put the x,y & disp into matrix, calculate the exact solution for i=1:(ndivl+1) for j=1:(ndivw+1) xh(j,i) = x((i-1)*(ndivw+1)+j); yh(j,i) = y((i-1)*(ndivw+1)+j); ux(j,i) = disp(2*((i-1)*(ndivw+1)+j)-1); uy(j,i) = disp(2*((i-1)*(ndivw+1)+j)); end end% plot mesh figure for i=1:numele plot(x(node([1:4,1],i)),y(node([1:4,1],i)),'k-o','linewidth',2,'markersize',12) hold on end for i=1:numnod text(x(i)+0.2,y(i)-0.3,[num2str(i)],'fontsize',12) endfigure surf(xh,yh,ux) shading interptitle('Patch Test','fontsize',12) xlabel('x(m)','fontsize',12);ylabel('y(m)','fontsize',12);zlabel('uh(m)','fontsize',12) set(gca,'fontsize',12)figure surf(xh,yh,uy) shading interptitle('Patch Test','fontsize',12) xlabel('x(m)','fontsize',12);ylabel('y(m)','fontsize',12);zlabel('uh(m)','fontsize',12) set(gca,'fontsize',12)1.2.3.4.5.6.7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.23.24.25.26.27.28.29.30.31.32.33.34.35.36.37.38.39.40.41.42.43.44.45.46.47.48.49.50.51.52.53.54.55.56.57.58.59.60.61.62.63.64.65.66.67.68.69.70.71.72.73.74.75.76.77.78.79.80.81.82.83.84.85.86.87.88.89.90.91.92.93.94.95.96.97.