singspan.m

% Example: singspan
% ~~~~~~~~~~~~~~~~~
% This example calculates the shear, moment,
% rotation, and deflection for a single span
% beams using discontinuity functions.
%
% Date is defined in the declaration statements
% below, where:
%
% Blen      - beam length
% Emod      - modulus of elasticity
% Binert    - beam moment of inertia
% No_segs   - number of locations to evaluate.
%             Calculations will be for every
%             (Blen/No_segs) location
% S_type    - type of beam supports
%             =1, pin-pin
%             =2, fixed-fixed
% Load_type - vector containing type of load
%             for each loading
%             =1, concentrated moment
%             =2, concentrated force
%             =3, step function, continuous
%             =4, step function, finite
%             =5, ramp function, continuous
% Load      - vector containing load value
%             for each loading
% Start_x   - vector containing starting
%             location for each loading
% End_x     - vector containing ending
%             location for each loading
%
% User m functions required: 
%    seval, constant, genprint
%----------------------------------------------

clear; 
Problem=3;
if Problem == 1
  %..pin-pin example
  Blen=192; Emod=29e6; Binert=394; 
  No_segs=16; S_type=1;
  Load_type=[2       4         ]; 
  Load=     [-12000  -416.66667];
  Start_x=  [132     60        ]; 
  End_x=    [132     132       ];
elseif Problem == 2
  %...fixed-fixed example
  Blen=300; Emod=29e6; Binert=60; 
  No_segs=15; S_type=2;
  Load_type=[3]; Load=[-20];
  Start_x=[150]; End_x=[300];
elseif Problem == 3
  %...Multiple loading example, pin-pin
  Blen=192; Emod=30e6; Binert=100; 
  No_segs=48; S_type=1;
  Load_type=[4    2     1      5  ]; 
  Load=     [-50  2400  24000  -75];
  Start_x=  [0    60    96     120]; 
  End_x=    [36   60    96     192];
elseif Problem == 4
  %...Boundary moments only, pin-pin
  Blen=10; Emod=30e6; Binert=100; 
  No_segs=10; S_type=1;
  Load_type=[1  1 ]; 
  Load=     [10 10];
  Start_x=  [0  10]; 
  End_x=    [0  10];
end
No_loads=length(Load_type);

%...Check valid ranges
if S_type < 1 | S_type > 2
  fprintf('\n\nInvalid support type error: ');
  fprintf('%g \n\n',S_type);
  error('Program aborted');
end
for i=1:No_loads
  if Load_type(i) < 1 | Load_type(i) > 5
    fprintf('\n\nInvalid load type error: ');
    fprintf('%g \n\n',Load_type(i));
    error('Program aborted');
  end
end

%...Define the coefficients of the 
%...discontinuity functions and handle 
%...the loads which may be applied
%...on the boundaries
M1=0; M2=0; 
Wl=0; Vl=0; Ml=0; Thetal=0; Deltal=0;
for j=1:No_loads
  %...Handle each load type
  if Load_type(j) == 1
    %...Concentrated moment
    if Start_x(j) == 0 
      %...Left boundary load
      M1=Load(j); Nterms(j)=0;
    elseif Start_x(j) == Blen
      %...Right boundary load
      M2=Load(j); Nterms(j)=0;
    else
      W(j,1)=0; V(j,1)=0; M(j,1)=Load(j);
      Theta(j,1)=Load(j); 
      Delta(j,1)=Load(j)/2;
      Nterms(j)=1; Power(j,1)=-2;
      A(j,1)=Start_x(j);
    end
  elseif Load_type(j) == 2
    W(j,1)=0; V(j,1)=Load(j); 
    M(j,1)=Load(j);
    Theta(j,1)=Load(j)/2; 
    Delta(j,1)=Load(j)/6;
    Nterms(j)=1; Power(j,1)=-1;
    A(j,1)=Start_x(j);
  elseif Load_type(j) == 3
    %...Step function, continuous
    W(j,1)=Load(j); V(j,1)=Load(j); 
    M(j,1)=Load(j)/2;
    Theta(j,1)=Load(j)/6; 
    Delta(j,1)=Load(j)/24;
    Nterms(j)=1; Power(j,1)=0;
    A(j,1)=Start_x(j);
  elseif Load_type(j) == 4
    %...Step function, finite
    W(j,1)=Load(j); V(j,1)=Load(j); 
    M(j,1)=Load(j)/2;
    Theta(j,1)=Load(j)/6; 
    Delta(j,1)=Load(j)/24;
    W(j,2)=-Load(j); V(j,2)=-Load(j); 
    M(j,2)=-Load(j)/2;
    Theta(j,2)=-Load(j)/6; 
    Delta(j,2)=-Load(j)/24;
    Nterms(j)=2; Power(j,1)=0; Power(j,2)=0;
    A(j,1)=Start_x(j); A(j,2)=End_x(j);
  elseif Load_type(j) == 5
    %...Ramp function 
    %...  (continuous and increasing)
    dx=End_x(j)-Start_x(j);
    W(j,1)=Load(j)/dx; V(j,1)=Load(j)/(2*dx); 
    M(j,1)=Load(j)/(6*dx);
    Theta(j,1)=Load(j)/(24*dx); 
    Delta(j,1)=Load(j)/(120*dx);
    Nterms(j)=1; Power(j,1)=1; 
    A(j,1)=Start_x(j);
  end
end

%...Evaluate the discontinuity functions at
%...x=Beam_length
[Wl,Vl,Ml,Thetal,Deltal]= ...
  seval(Blen,Nterms,No_loads,Power,A, ...
        W,V,M,Theta,Delta);

%...Determine constants of integration
[Cv,Cm,Ctheta,Cdelta]= ...
  constant(M1,M2,Blen,S_type, ...
           Vl,Ml,Thetal,Deltal);

%...Evaluate the discontinuity functions 
%...along the beam
N1=No_segs+1; x=linspace(0,Blen,N1);
WX=0; VX=0; MX=0; ThetaX=0; DeltaX=0;
for j=1:N1
  [WX,VX,MX,ThetaX,DeltaX]= ...
    seval(x(j),Nterms,No_loads,Power,A, ...
    W,V,M,Theta,Delta);
  %...Combine with boundary contributions
  Wx(j)=WX;
  Vx(j)=VX+Cv;
  Mx(j)=MX+Cv*x(j)+Cm;
  Thetax(j)=(ThetaX+Cv*x(j)^2/2+ ...
            Cm*x(j)+Ctheta)/(Emod*Binert);
  Deltax(j)=(DeltaX+Cv*x(j)^3/6+ ...
            Cm*x(j)^2/2+Ctheta*x(j)+ ...
            Cdelta)/(Emod*Binert);
end

%...Output results
fprintf('\n\n  Single Span Beam Analysis');
fprintf('\nUsing Discontinuity Functions');
fprintf('\n-----------------------------');
fprintf('\n\nProblem Specifications:');
fprintf('\n  Beam length:        %g',Blen);
fprintf('\n  E:                  %g',Emod);
fprintf('\n  I:                  %g',Binert);
fprintf('\n  # of beam segments: %g',No_segs);
fprintf('\n  # of loads:         %g',No_loads);
fprintf('\n  Type of support:    %g',S_type);
if S_type == 1
  fprintf(' (pin-pin)');
elseif S_type == 2
  fprintf(' (fixed fixed)');
end
fprintf('\n\nSummary for each load:');
for i=1:No_loads
  fprintf('\n  Load %g',i);
  fprintf('\n    Magnitude:  %g',Load(i));
  fprintf('\n    Starting x: %g',Start_x(i));
  fprintf('\n    Ending x:   %g',End_x(i));
  fprintf('\n    Load type:  %g',Load_type(i));
  if Load_type(i) == 1
    fprintf(' (concentrated moment)');
  elseif Load_type(i) == 2
    fprintf(' (concentrated force)');
  elseif Load_type(i) == 3
    fprintf(' (step function, continuous)');
  elseif Load_type(i) == 4
    fprintf(' (step function, finite)');
  elseif Load_type(i) == 5
    fprintf(' (ramp function, continuous)');
  end
end
fprintf('\n\nBoundary results:');
fprintf('\n  Left end of beam');
fprintf('\n    Shear:      %g',Vx(1));
fprintf('\n    Moment:     %g',Mx(1));
fprintf('\n    Rotation:   %g',Thetax(1));
fprintf('\n    Deflection: %g',Deltax(1));
fprintf('\n  Right end of beam');
fprintf('\n    Shear:      %g',Vx(N1));
fprintf('\n    Moment:     %g',Mx(N1));
fprintf('\n    Rotation:   %g',Thetax(N1));
fprintf('\n    Deflection: %g',Deltax(N1));
fprintf('\n\n');

%...Plot the results
clf;
plot(x,Wx,'-',x,Wx,'o');
  xlabel('x'); ylabel('Load');
  title('Load Diagram'); drawnow;
% genprint('bload');
  disp('Press a key to continue'); pause;

clf;
subplot(2,1,1);
  plot(x,Vx,'-',x,Vx,'o');
  xlabel('x'); ylabel('Shear Force');
  title('Shear Diagram'); 
subplot(2,1,2);
  plot(x,Mx,'-',x,Mx,'o');
  xlabel('x'); ylabel('Moment');
  title('Moment Diagram'); drawnow;
% genprint('bshrmom');
  disp('Press a key to continue'); pause;

clf;
subplot(2,1,1);
  plot(x,Thetax,'-',x,Thetax,'o')
  xlabel('x'); ylabel('Rotation (radians)');
  title('Rotation Diagram'); 
subplot(2,1,2);
  plot(x,Deltax,'-',x,Deltax,'o')
  xlabel('x'); ylabel('Deflection');
  title('Deflection Diagram'); drawnow;
% genprint('brotdef');