truss.m

% Example: truss
% ~~~~~~~~~~~~~~
% This example analyzes the effects of a load
% applied to the point of interconnect for
% two rods forming a truss.  The method
% is based on iterating until the solution
% reaches a prescribed convergence criteria.
%
% Data is defined in the declaration statements
% below, where:
%
% P        - downward vertical load
% L_ab_o   - original horizontal span for AB
% L_ac_o   - original horizontal span for AC
% A_ab     - area of rod AB
% A_bc     - area of rod AB
% E_ab     - modulus of elasticity for AB
% E_bc     - modulus of elasticity for AC
% eps      - solution tolerance
% max_iter - max # iterations
%
% User m functions required: 
%    genprint
%----------------------------------------------

clear;
%...Input definitions
Problem=1;
if Problem == 1 
  P=1000; L_ab_o=50; L_ac_o=50;
  A_ab=2.0; A_bc=2.0; E_ab=3000; E_bc=3000; 
  eps=1.0e-9; max_iter=20;
end

%...Initialize
raddeg=180/pi; pi2=pi/2;

%...Calculate initial values
L_bc_o=sqrt(L_ab_o^2+L_ac_o^2); 
Theta_a_o=pi2; Theta_c_o=acos(L_ac_o/L_bc_o);
Theta_b_o=pi2-Theta_c_o;

%...Initialize loop parameters
Alpha(1)=0; Beta(1)=pi2-Theta_c_o;
error=10^10; i=0;

fprintf('\n\nIterative Analysis of Truss');
fprintf(  '\n---------------------------');
fprintf('\n\nOriginal Properties');
fprintf(  '\n  Load (P):    %g',P);
fprintf(  '\n  Area AB:     %g',A_ab);
fprintf(  '\n  Area BC:     %g',A_bc);
fprintf(  '\n  E for AB:    %g',E_ab);
fprintf(  '\n  E for BC:    %g',E_bc);
fprintf(  '\n  Length AB:   %g',L_ab_o);
fprintf(  '\n  Length AC:   %g',L_ac_o);
fprintf(  '\n  Length BC:   %g',L_bc_o);
fprintf('  \n  Angle at A:  %g (degrees)', ...
  Theta_a_o*raddeg);
fprintf('  \n  Angle at B:  %g (degrees)', ...
  Theta_b_o*raddeg);
fprintf('  \n  Angle at C:  %g (degrees)', ...
  Theta_c_o*raddeg);
fprintf('  \n  Angle Alpha: %g (degrees)', ...
  Alpha(1)*raddeg);
fprintf('  \n  Angle Beta:  %g (degrees)', ...
  Beta(1)*raddeg);
fprintf(  '\n  Epsilon:     %g',eps);
fprintf(  '\n  max iter:    %g',max_iter);
fprintf('\n');

%..................................
%...Loop until convergence achieved
%..................................
while error>eps
  i=i+1;
  if i>max_iter 
    fprintf('\n\n Max iterations exceeded\n');
    error('\nProgram aborted');
  end
  %...Set up simultaneous equations and solve
  a=[-cos(Alpha(i)) cos(Beta(i)); ...
      sin(Alpha(i)) sin(Beta(i))];
  b=[0;P]; F=a\b;
  F_ab(i)=F(1); F_bc(i)=F(2);
  %...Calculate new lengths
  delta_ab(i)=F_ab(i)/(A_ab*E_ab);
  delta_bc(i)=-F_bc(i)/(A_bc*E_bc);
  L_ab(i)=L_ab_o*(1+delta_ab(i));
  L_bc(i)=L_bc_o*(1+delta_bc(i));
  %...Use triangle definitions
  S=0.5*(L_ab(i)+L_bc(i)+L_ac_o);
  K=sqrt(S*(S-L_ab(i))*(S-L_bc(i))*(S-L_ac_o));
  Hc=2*K/L_ac_o;
  %...New angles
  Alpha(i+1)=acos(Hc/L_ab(i));
  Beta(i+1)=acos(Hc/L_bc(i));
  Theta_a(i)=pi2-Alpha(i+1);
  Theta_b(i)=Alpha(i+1)+Beta(i+1);
  Theta_c(i)=pi-Theta_a(i)-Theta_b(i);
  %...New location of pin B
  delta_x(i)=L_ab(i)*cos(Alpha(i+1))-L_ab_o;
  delta_y(i)=-(L_ac_o-L_bc(i)*sin(Beta(i+1)));
  error=abs(Alpha(i)-Alpha(i+1));
  err(i)=error;

  fprintf('\n\nIteration # %g',i);
  fprintf('\n  Force AB : %g',F_ab(i));
  fprintf('\n  Force BC : %g',-F_bc(i));
  fprintf('\n  Length AB: %g',L_ab(i));
  fprintf('\n  Length BC: %g',L_bc(i));
  fprintf('\n  Delta AB:  %g',delta_ab(i));
  fprintf('\n  Delta BC:  %g',delta_bc(i));
  fprintf('\n  Delta x:   %g',delta_x(i));
  fprintf('\n  Delta y:   %g',delta_y(i));
  fprintf('\n  Theta A:   %g (degrees)', ...
          Theta_a(i)*raddeg);
  fprintf('\n  Theta B:   %g (degrees)', ...
          Theta_b(i)*raddeg);
  fprintf('\n  Theta C:   %g (degrees)', ...
          Theta_c(i)*raddeg);
  fprintf('\n  Alpha:     %g (degrees)', ...
          Alpha(i+1)*raddeg);
  fprintf('\n  Beta:      %g (degrees)', ...
          Beta(i+1)*raddeg);
end
fprintf('\n\n');

%...Plot results
clf; iter=1:i; F_bc=-F_bc;
subplot(2,1,1);
  plot(iter,F_ab(1:i),'o',iter,F_ab(1:i),'-')
  title('Force in AB');
  xlabel('Iteration'); ylabel('Force in AB');
subplot(2,1,2);
  plot(iter,F_bc(1:i),'o',iter,F_bc(1:i),'-')
  title('Force in BC');
  xlabel('Iteration'); ylabel('Force in BC');
  drawnow;
% genprint('force');
  disp('Press key to continue'); pause;

clf; Alpha=Alpha*raddeg;
subplot(2,1,1);
  plot(iter,Alpha(1:i),'o', ...
       iter,Alpha(1:i),'-')
  title('Angle Alpha'); xlabel('Iteration'); 
  ylabel('Alpha (degrees)');
subplot(2,1,2);
  semilogy(iter,err(1:i),'o',iter,err(1:i),'-')
  title('Error Measure');
  xlabel('Iteration'); ylabel('Error');
  drawnow;
% genprint('error');
  disp('Press key to continue'); pause;

clf; 
plot(delta_x(1:i),delta_y(1:i),'o', ...
     delta_x(1:i),delta_y(1:i),'-') 
  axis('equal');
  title('Displacement of Point B');
  xlabel('Delta in x'); ylabel('Delta in y');
  drawnow;
% genprint('deltab');