% Axial1, stepped bar, 2 elements, 3 nodes % clear % variables format short e % E format % format long % double precision p1 = 1000 % load at node 1 % define elastic modulus, length, and diameter % the following 3 vectors need one entry per element ex = [ 29 11 ]*1e6 % elastic modulus len = [10 12] % element lengths dia = [1.1 1.6] % element diameters a = pi/4*dia .^2 % area % define the stiffness matrix for each element k12 = a(1)*ex(1)/len(1) % stiff const k23 = a(2)*ex(2)/len(2) k22 = k12 + k23 beta = k12*1e7 % Zienkiewicz factor stiff = [ k12 -k12 0 % stiffness matrix -k12 k22 -k23 0 -k23 k23 ] df = [ p1; 0; 0*beta ] % combined d, F zienk = stiff; zienk(3,3)=zienk(3,3)+beta % Zienkiewicz method displ = inv(zienk) * df; % displacements force = stiff * displ; % force vector for i = 1:length(df)-1 % element displacements eldis(i) = displ(i+1) - displ(i); end stress = ex .* eldis ./ len; % element stresses % print the results fprintf('\n Axial1 ') fprintf('\n Your name, MENG 421, ') disp(date) fprintf('\n Node Displacement Force') fprintf('\n Inches Lb \n') for i = 1:length(df) fprintf(' %2.0f %11.3e %7.0f \n',i,displ(i),force(i)) end fprintf('\n Elem Stress') fprintf('\n Psi \n') for i = 1:length(df)-1 fprintf(' %2.0f %7.0f \n',i,stress(i)) end % End of Matlab code % The following is the output from this program % Axial1 % Your name, MENG 421, 05-Feb-2004 % Node Displacement Force % Inches Lb % 1 9.054e-004 1000 % 2 5.426e-004 0 % 3 3.628e-011 -1000 % Elem Stress % Psi % 1 -1052 % 2 -497