单边简支梁有限元程序_MATLAB_FEM_简支梁

% ------------------------------------------------------------------------| % | % MAE4700-5700, Finite Element Analysis for Mechanical & Aerospace Design | % | % Copyright: Cornell University (this software should not be used without | % written permission) | % | % Authors: N. Zabaras (zabaras@cornell.edu) & Xiang Ma (xm25@cornell.edu) | % | % ------------------------------------------------------------------------| % function postprocessor; include_variables; % In case you need to compare with analytical stiffness calculations, % you can display the global stiffness matrix by converting the sparse % matrix to a full matrix as shown below. % % K = full(K) % % You can also print f, d or r if needed for comparisons % f % d % r % The remaining of this function finishes the postprocessing operations displacement = []; moment = []; shear = []; for e = 1:nel % loop over elements to plot % displacements, moments and shear forces de = zeros(nen*ndof,1); % extract element nodal displacements for i = 1 : nen for a = 1:ndof loci = (i-1)*ndof + a; glbi = (Elems(i,e) - 1) * ndof + a; de(loci) = d ( glbi ); end end IENe = Elems(:,e); % extract element connectivity information xe = Nodes(IENe); % extract element coordinates L = xe(2) - xe(1); % element length J = L / 2; % compute Jacobian for transformation nplot = 10; % number of points to plot % Compute displacements, moments and shear forces x = linspace(xe(1),xe(2),nplot); % equally distributed points within each element xplot = (2*x - xe(1) - xe(2))/L; % transform to local coordinate for i = 1:nplot xi = xplot(i); % current coordinate % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 由于考虑剪切变形的影响,要适当的改动一下形函数. % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 原来的形函数 % N = [1/4*(1-xi)^2*(2+xi), L/8*(1-xi)^2*(1+xi), ... % 1/4*(1+xi)^2*(2-xi), L/8*(1+xi)^2*(xi-1)]; % shape functions % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 考虑剪切影响之后的形函数 N = [1/2*(1/4*(1-xi)^2*(2+xi)+1/2*(1-xi)), L/8*(1-xi)^2*(1+xi), ... 1/2*(1/4*(1+xi)^2*(2-xi)+1/2*(xi+1)), L/8*(1+xi)^2*(xi-1)]; B = [3/2*xi, L*(3/4*xi - 1/4), -3/2*xi,L*(3/4*xi + 1/4)]/J^2; % 2nd derivative of N S = [3/2, 3/4*L, -3/2, 3/4*L]/J^3; % 3nd derivative of N EIe = EI(e); % Young's modulus x I % The following commands append as a new row the vector consisting % of the x location & corresponding displacement/moment/shear force. % This allows easy ploting of the displacements, moments and shear % force as a function of location in the beam. displacement = [displacement;[x(i), N*de]]; % displacement output moment = [moment; [x(i), EIe*B*de]]; % moment output shear = [shear; [x(i), -EIe*S*de]]; % Shear force output end end % plot displacements, moment and shear forces [x_plot,S_ex,M_ex,w_ex]=exact; figure(2) plot(displacement(:,1),displacement(:,2),'-.r'); hold on; plot(x_plot,w_ex,'-k'); legend('FE','Exact Solution','Location','SouthWest'); hold on; ylabel('displacement'); xlabel('x'); title('Displacements: FE versus analytical beam solutions'); axis([0 12 -1.5 0]); figure(3) plot(moment(:,1),moment(:,2),'-.r'); hold on; plot(x_plot,M_ex,'-k'); legend('FE','Exact Solution','Location','NorthWest'); hold on; ylabel('moment'); xlabel('x'); title('Moments: FE versus analytical beam solutions'); axis([0 12 -300 100]); figure(4) plot(shear(:,1),shear(:,2),'-.r'); hold on; plot(x_plot,S_ex,'-k'); legend('FE','Exact Solution'); hold on; ylabel('shear'); xlabel('x'); title('Shear: FE versus analytical beam solutions'); axis([0 12 -35 -10]);