% Moving force across a single span beam, solved analytically
clear all
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% School of Mechanical and Manufacturing Engineering, University of New
% South Wales
% Author:   Gareth Forbes
% Date created: 1/2/08
% Date last modified: 10/11/08
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
% Revision number: 1
% Description of changes to latest revision:
% Updated to allow animated plotting for any input values. Change of plot
% y-axis labels to show relative displacement. Solution for all input
% values with 101 steps in direction and time plane.
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
% Description of Script: 
% Creates the response of a single span beam under the influence of a
% moving force according to the analytical equations derived in [1]. Note,
% the formulation used here is only valid for light damping at various
% speeds, and with no damping at non-critical speeds. Two movies, at 
% different moving load speeds have been created and uploaded at [2] and
% [3]
%
%
% References:
% [1]   Fryba, L., Vibration of solids and structures under moving loads. 
% 3rd ed. 1999, London: Thomas Telford. xxvii,494 p.
% [2]   http://www.youtube.com/watch?v=XhU_IGx2CdU
% [3]   http://www.youtube.com/watch?v=jXIiQlnzChY
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input varibales for analysis
k = 5; %number of terms in expansion
P = 100; %load (N)
beta = 0.1; % non-dimensional damping parameter
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% material and geometric constants
L = 30480; %length of beam (mm)
E = 200000; %youngs modulus (MPa)
b = 2438; % width of cross section
h = 20; % height of cross section
A = b*h;% cross sectional area of beam (mm^2)
rho = 7.8E-9;% density of beam material (kg/mm^3)
I = (b*h^3)/12;% second moment of area (mm^4)
% critical speed (resonance of first mode)
cr = (pi/L)*(sqrt(E*I/rho/A));
% speed of load
c = 0.5*cr;% (mm/s)
alpha = c/cr; % non-dimensional speed parameter
omega = pi*c/L;
omegaj = ([1:k]*pi/L).^2*(sqrt(E*I/rho/A));
omegab = beta*omegaj(1);
% length vector;
x1 = 0:L/100:L;
x1 = x1';
x = repmat(x1,1,length(x1));
% time vector
tt = L/c;
step = tt/(length(x1)-1);
t1 = 0:step:tt;
t = repmat(t1,length(t1),1);
% static deflection of beam at mid span
u0 = (P*L^3)/(48*E*I); % (mm)
n = round(alpha);
if beta == 0;
    u2 = 0;
else
    if n > 0;
        if abs(n-alpha)<0.01;
            u2 = (u0/2/n^4)*(exp(-omegab*t).*sin(n*omega*t)-(n^2/beta)*cos(n*omega*t).*(1-exp(-omegab*t))).*sin(n*pi*x/L);
        else
            u2 = 0;
        end
    else
        u2 = 0;
    end
end
u = 0;
for j = 1:k;
    if abs(j-alpha)<0.01;
        u1 = 0;
    else
        u1 = u0*(1/(j^2*(j^2-alpha^2)))*(sin(j*omega*t)-(alpha*exp(-omegab*t)/j).*sin(omegaj(j)*t)).*sin(j*pi*x/L);
    end
    u = u +u1;
end
u3 = u +u2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creation of movie
for i = 0:1:length(x1)-1;
    ax1 = (0:100)*L/100;
    subplot(2,1,1), plot(ax1,-u3(:,i+1)/u0,'linewidth',1.5,'color','r');
    xlabel('beam length (mm)')
    ylabel('deflection (u/u0)')
    title('Beam deflection under load')
    axis([-0.03*L 1.03*L -1.7 0.6])
    line([i*L/100 i*L/100],[0-(u3(i+1,i+1)/u0) 0.5-(u3(i+1,i+1)/u0)],'Color','k','linewidth',2)
    line(i*L/100,(-u3(i+1,i+1))/u0,'Color','k','marker','V','MarkerSize',3,'linewidth',3)
    % location of first supports
    a = [1 0];
    b = [L+1 0];
    % height of support
    h = 0.5;
    % support 1
    line([a(1) a(1)+h*L*0.03],[a(2) -h+a(2)],'linewidth',1);
    line([a(1)-h*L*0.03 a(1)+h*L*0.03],[-h+a(2) -h+a(2)],'linewidth',1);
    line([a(1)-h*L*0.03 a(1)],[-h+a(2) a(2)],'linewidth',1);
    % support 2
    line([b(1) b(1)+h*L*0.03],[b(2) -h+b(2)],'linewidth',1);
    line([b(1)-h*L*0.03 b(1)+h*L*0.03],[-h+b(2) -h+b(2)],'linewidth',1);
    line([b(1)-h*L*0.03 b(1)],[-h+b(2) b(2)],'linewidth',1);
    %
    ax2 = (0:100)*tt/100;
    subplot(2,1,2), plot(ax2(1:i+1),-u3(51,1:i+1)/u0,'linewidth',1.5,'color','r');
    xlabel('time (s)')
    ylabel('deflection (u/u0)')
    title('Deflection at the centre of the beam')
    axis([-tt/20 21*tt/20 -1.7 0.6])
    M(i+1) = getframe(gcf);
end
%movie2avi(M,'a','compression','none','quality',100)