% Moving force across a single span beam, using the finite element method
% in matlab
clear all
close all
tic
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% School of Mechanical and Manufacturing Engineering, University of New
% South Wales
% Author:   Gareth Forbes
% Date created: 1/2/08
% Date last modified: 12/8/08
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
% Description of Script: 
% Solves the problem of a moving load across a single span simply supported 
% beam, with the use of the finite element method. A transient direct 
% central difference explicit integration scheme is used for solution, note 
% this is often not the best solution method due to the amount of 
% computations need for a direct solution method, transformation into 
% the modal domain will often produce greater effeciency in the solution.
%
%
% References:
% [1]	Cook, R.D., et al., Concepts and applications of finite element
% analysis. 4th ed. 2002, New York, NY: Wiley. xvi, 719 p.
% [2]	Wu, J.-J., A.R. Whittaker, and M.P. Cartmell, Use of finite element
% techniques for calculating the dynamic response of structures to moving 
% loads. Computers and Structures, 2000. 78(6): p. 789-799.
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
%speed of moving load mm/s
V = 0.25*9.1845e+004 % change this value to solve for different speeds
% -------------------------------------------------------------------------
% define material properties
rho = 7.8E-9; %density
b = 50; %width (mm)
h = 20; %height (mm)
A = b*h; %cross section 
I = (b*h^3)/12; %second moment of area
E = 200000; %youngs modulus (MPa)
Lt = 1000; %length (mm)
%number of elements
EL = 30; % note the use of more than approx. 50 elements causes the 
%solution time to become very large. The use of vector assignment instead
%of for loops may correct this
% number of nodes
N = EL + 1;
% forceing function in newtons, single pont force
P = -1;
% first critical speed of simply supported beam (mm/s)
SC = pi/Lt*sqrt(E*I/rho/A)
% load speed to critical speed ratio
crat = V/SC
% manual time step
tstepm = 1000; % minimum time step that is allowed
% all elements are defined as having the same length therefore the element
% lenght is;
L = Lt/EL;
% elemental mass matrix
%me = (rho*A*L/420)*[156 22*L 54 -13*L; 22*L 4*L^2 13*L -3*L^2; 54 13*L 156
%-22*L; -13*L -31*L^2 -22*L 4*L^2]; % distributed mass matrix
me = (rho*A*L/2)*[1 0 0 0;0 1/12*L^2 0 0;0 0 1 0; 0 0 0 1/12*L^2]; % lumped mass matrix
% elemental stiffness matrix
ke = (E*I/(L^3))*[12 6*L -12 6*L; 6*L 4*L^2 -6*L 2*L^2; -12 -6*L 12 -6*L; 6*L 2*L^2 -6*L 4*L^2]; % distributed stiffness matrix
% create mass and stiffness matrix, note all propeties are constant
M1 = zeros(2*(EL+1),2*(EL+1));
M2 = zeros(2*(EL+1),2*(EL+1));
for i = 2:2:2*(EL+1);
    M2 = M2 + M1;
    M1 = zeros(2*(EL+1),2*(EL+1));
    M1(i-1:i+2,i-1:i+2) = me;
end
K1 = zeros(2*(EL+1),2*(EL+1));
K2 = zeros(2*(EL+1),2*(EL+1));
for i = 2:2:2*(EL+1);
    K2 = K2 + K1;
    K1 = zeros(2*(EL+1),2*(EL+1));
    K1(i-1:i+2,i-1:i+2) = ke;
end
% enforce boundary conditions
%  delete first and second last row and column
M = zeros(2*EL,2*EL);
K = zeros(2*EL,2*EL);
M(1:2*EL-1,1:2*EL-1) = M2(2:2*EL,2:2*EL);
M(2*EL,1:2*EL-1) = M2(2*(EL+1),2:2*EL);
M(1:2*EL-1,2*EL) = M2(2:2*EL,2*(EL+1));
M(2*EL,2*EL) = M2(2*(EL+1),2*(EL+1));
K(1:2*EL-1,1:2*EL-1) = K2(2:2*EL,2:2*EL);
K(2*EL,1:2*EL-1) = K2(2*(EL+1),2:2*EL);
K(1:2*EL-1,2*EL) = K2(2:2*EL,2*(EL+1));
K(2*EL,2*EL) = K2(2*(EL+1),2*(EL+1));
%create damping matrix, using proportional damping
alpha = 0.0002;
beta = 0.0002;
cm = alpha*M; 
ck = beta*K;
ct = cm + ck;
% calculate the minimum time step
% solve eigen value problem such that deltat < 2/natfmax. see [1]
[V1,D1] = eig(K,M);
ttotal = Lt/V; % length of time record in seconds %0.1
E1 = sqrt(D1);
deltat1 = 2/(max(max(real(E1))));
deltat = (1/(1/deltat1+1));
tstep = max(round(ttotal/deltat),tstepm);

% create force vector for time step. see [2]
for i = 1:tstep;
    for j = 1:EL+1;
        xp = V*i*deltat;
        s = fix(xp/L)+1;
        zeta = (xp - (s-1)*L)/(L);
        N1 = 1 - 3*zeta^2 + 2*zeta^3;
        N2 = (zeta - 2*zeta^2 + zeta^3)*L;
        N3 = 3*zeta^2 - 2*zeta^3;
        N4 = (-zeta^2 + zeta^3)*L;
            if s == EL+1;
                F(j,i) = 0;
                MO(j,i) = 0;
            elseif j == s ;
                F(j,i) = P*N1;
                MO(j,i) = P*N2;
            elseif j == s+1;
                F(j,i) = P*N3;
                MO(j,i) = P*N4;
            else
                F(j,i) = 0;
                MO(j,i) = 0;
            end
    end
end
% set first and last rows of load vectors to zero, so that he load
% gradually enters and exits the structure, ie. does not jump on and off,
% this could be removed if a jump on and off is desired
F(1,:) = zeros;
F(end,:) = zeros;
MO(1,:) = zeros;
MO(end,:) = zeros;
    % create external load vector
    REXT2 = zeros(2*(EL+1),tstep);
    for i = 1:EL+1;
        REXT2(2*i-1,:) = F(i,:);
        REXT2(2*i,:) = MO(i,:);
    end
% enforce boundary conditions
% delete first and second last row and column
REXT = zeros(2*EL,tstep);
REXT(1:2*EL-1,:) = REXT2(2:2*EL,:);
REXT(1:2*EL-1,:) = REXT2(2:2*EL,:);
REXT(2*EL,:) = REXT2(2*(EL+1),:);
% use a central difference direct explicit integration method to find the
% response history
D = zeros(2*EL,tstep);
D2C = ((1/deltat^2)*M + (1/(2*deltat))*ct);
D1C = (2/deltat^2)*M;
D0C = ((1/deltat^2)*M - (1/(2*deltat))*ct);
RINT = K;
for i = 2:tstep-1;
    % enforce boundary conditions
    %D(1,i) = zero;
    %D(2*EL+1,i) = zero;
    D(:,i+1) = D2C\(REXT(:,i) - RINT*D(:,i) + D1C*D(:,i) - D0C*D(:,i-1));
end
% creat matrix of full element displacements
Df = zeros(2*(EL+1),tstep);
Df(2:2*EL,:) = D(1:2*EL-1,:);
Df(2*(EL+1),:) = D(2*EL,:);
% create matrix of just displacements and rotations
for i = 1:EL+1;
        Dx(i,:) = Df(2*i-1,:);
        Dtheta(i,:) = Df(2*i,:);
end
figure
ax = (0:length(Dx)-1)*ttotal/length(Dx);
fe = round(Lt/V/deltat); %time step when load leaves structure
plot(ax(fe),Dx(round(EL/2+1),fe),'ro')
hold
plot(ax,Dx(round(EL/2+1),:))
xlabel('time')
ylabel('displacement (mm)')
title('Displacement of centre element for whole time period')
legend('markes when load leaves structure')
hold
figure
% creation of movie. Note movie script is quite messy. Also the use of less
% than approx. 50 elements causes the movie to be a little stilted, however
% solution time for more than approx. 50 elements is quite long. Vector
% assignment instead of for loops may solve this problem.
for i = 0:1:100;
    ae = fe/100;
    ax1 = (0:EL)*Lt/EL;
    subplot(2,1,1), plot(ax1,Dx(:,round(i*ae)+1),'linewidth',1.5,'color','r');
    xlabel('beam length (mm)')
    ylabel('deflection (mm)')
    title('Beam deflection under load')
    ge(i+1) = Dx(round(round(ae*i)*deltat*V*EL/Lt+1),round(i*ae)+1);
    axis([-30 1030 -0.017/3 0.006/3])
    line([10*i 10*i],[ge(i+1) 0.0015+ge(i+1)],'Color','k','linewidth',2)
    line(i*10,ge(i+1),'Color','k','marker','V','MarkerSize',3,'linewidth',3)
    % location of first supports
    a = [1 0];
    b = [1001 0];
    % height of support
    h = 0.002;
    % support 1
    line([a(1) a(1)+h*6000],[a(2) -h+a(2)],'linewidth',1);
    line([a(1)-h*6000 a(1)+h*6000],[-h+a(2) -h+a(2)],'linewidth',1);
    line([a(1)-h*6000 a(1)],[-h+a(2) a(2)],'linewidth',1);
    % support 2
    line([b(1) b(1)+h*6000],[b(2) -h+b(2)],'linewidth',1);
    line([b(1)-h*6000 b(1)+h*6000],[-h+b(2) -h+b(2)],'linewidth',1);
    line([b(1)-h*6000 b(1)],[-h+b(2) b(2)],'linewidth',1);
    ax2 = (0:100)*Lt/V/100;
    subplot(2,1,2), plot(ax2(1:i+1),Dx(EL/2+1,1:round(ae):i*round(ae)+1),'linewidth',1.5,'color','r');
    xlabel('time (s)')
    ylabel('deflection (mm)')
    title('Deflection at the centre of the beam')
    axis([-Lt/V/20 21*Lt/V/20 -0.017/3 0.006/3])
    Mov(i+1) = getframe(gcf);
end
toc
%movie2avi(Mov,'a','compression','none','quality',100)