% iTFM code for calculations of gas temperature in improved formulation of travelling fires  
% by Egle Rackauskaite and Guillermo Rein, Imperial College London, UK.
% Based on the journal paper (doi:10.1016/j.istruc.2015.06.001): E Rackauskaite, C Hamel, A Law, G Rein, Improved formulation of travelling fires and application to concrete and steel structures, Structures, 2015.
% Contact authors at g.rein@imperial.ac.uk and reingu@gmail.com
% Work funded by Engineering and Physical Sciences Research Council and Arup
% Code published under a Creative Commons license CC BY 4.0 

clear all
close all
syms R

% % =======================================================================
% % =============================             =============================
% % =============================    INPUT    =============================
% % =============================             =============================
% % =======================================================================

% === ROOM PARAMETERS ===

WIDTH = 18;       % WIDTH OF THE COMPARTMENT (m)
HEIGHT = 3.5;     % HEIGHT OF THE COMPARTMENT (m)
L = 24;           % LENGTH OF THE COMPARTMENT (m) - FIRE PATH

dx = 0.25;        % LENGTH - CALCULATE TEMPERATURES EVERY dx (m)
dt = 10;          % TIME STEP (s)

% === FIRE PARAMETERS ===

FSIZE = 0.1;      % DIMENSIONLESS FIRE SIZE ( = LENGTH OF FIRE / LENGTH OF FIRE PATH)
QII = 500;        % HEAT RELEASE RATE PER UNIT AREA (kW/m2)
QF = 570000;      % FUEL LOAD DENSITY (kJ/m2)
T_ROOM = 20;      % ROOM TEMPERATURE (°C)
T_NF = 1200;      % PEAK FIRE TEMPERATURE (°C)

% === FLAPING ANGLE ===

FANGLE = 6.5;     % FLAPPING ANGLE (°)

% === FIGURES ===

ti = [30 60 120]; % TIMES OF INTEREST (min)
xi = [0 13 20];   % LOCATIONS OF INTEREST (m)

% % =======================================================================
% % =====================                             =====================
% % =====================    PARAMETER CALCULATIONS   =====================
% % =====================                             =====================
% % =======================================================================

% === FIRE PARAMETERS ===

TB = QF / QII;                   % LOCAL BURNING TIME (s)
TTOTAL = TB * ( 1 / FSIZE + 1 ); % TOTAL FIRE DURATION (s)
LF = L * FSIZE;                  % LENGTH OF THE FIRE (m)
S = LF / TB;                     % SPREAD RATE OF THE FIRE (m/s)
TLIMIT = 1.3 * TTOTAL;           % DURATION FOR WHICH TEMPERATURES ARE CALCULATED (s)

% === NUMBER OF LOCATIONS OF INTEREST AND TIME STEPS ===

tt = 0:dt:TLIMIT;  % ARRAY OF TIME STEPS (s)
xx = 0:dx:L;       % ARRAY OF LOCATIONS (m)

% % =======================================================================
% % ============================               ============================
% % ============================    FLAPPING   ============================
% % ============================               ============================
% % =======================================================================

F = LF + 2 * HEIGHT * tan( FANGLE * pi / 180 ); % FLAPPING LENGTH (m)
Q = WIDTH * LF * QII;                           % TOTAL HEAT RELEASE RATE(kW)
T = 5.38 / HEIGHT * ( Q / R )^(2/3) + T_ROOM;   % ALPERT'S EQUATION (1972) doi: 10.1007/BF02590543

R0 = Q * ( 5.38 / ( HEIGHT * ( T_NF - T_ROOM ) ) )^(3/2); % INTERSECTION POINT OF Tmax (ALPERT'S EQUATION) AND T_NF (e.g. 1200°C)
if R0 > F / 2
    T_NEAR_FIELD = T_NF;
else
    RX1 = max( 0 , R0 - LF / 2 ); 
    RX2 = max( LF / 2 , R0 );
    R2 = F / 2; % DISTANCE FROM THE CENTRE OF THE FIRE TO FLAPPING RANGE EDGE
    T2 = int(T,R,RX2,R2); % SUM OF TEMPERATURE EXCLUDING NEAR-FIELD
    Ttotal = 2 * T2 + 2 * T_NF * RX1 + T_NF * LF; % SUM OF TEMPERATURES OVER FLAPPING LENGTH
    T_NEAR_FIELD = double( Ttotal / F ); % MAX NEAR-FIELD TEMPERATURE INCLUDING FLAPPING
end

% % =======================================================================
% % =========================                      ========================
% % =========================   GAS TEMPERATURES   ========================
% % =========================                      ========================
% % =======================================================================

T_MAX = zeros(numel(xx),numel(tt)); % GAS TEMPERATURE T_MAX(location,time) (°C)
for x = 1:numel(xx)
    for t = 1:numel(tt)
        L_TRAILING = S * tt(t);
        if L_TRAILING < L
            Lt = L_TRAILING;                      % TIME DEPENDENT FIRE FRONT (m)
            Ft = min([ FSIZE , L_TRAILING / L ]); % TIME DEPENDENT FIRE SIZE (dimensionless)
        elseif abs( L_TRAILING - L ) < 1e-14      % [L_TRAILING(t) = L] TO ACCOUNT FOR FLOATING POINT CALCULATION i.e. 1-1 = 0. 1e-14 = tolerance
            Lt = L_TRAILING;                      % TIME DEPENDENT FIRE FRONT (m)
            Ft = min([ FSIZE , L_TRAILING / L ]); % TIME DEPENDENT FIRE SIZE (dimensionless)
        else 
            Lt = L;                               % TIME DEPENDENT FIRE FRONT (m) when fire reaches the end of the fire path
            Ft = 1 + ( LF - L_TRAILING ) / L;     % TIME DEPENDENT FIRE SIZE (dimensionless) when fire reaches the end of the fire path
            if  abs( 1 + ( LF - L_TRAILING ) / L ) < 1e-14 % TO ACCOUNT FOR FLOATING POINT CALCULATION i.e. 1-1 = 0. 1e-14 = tolerance
                     Ft = 0;
            end
        end
        r = abs( xx(x) - Lt + 0.5 * L * Ft ); % DISTANCE FROM FIRE CENTRE
        
        %=========== CALCULATION OF TMAX ================================
        if tt(t) > TTOTAL             % ROOM TEMPERATURE AT THE END OF FIRE
            T_MAX(x,t) = T_ROOM;
        elseif r < 0.5 * L * Ft       % NEAR-FIELD
            T_MAX(x,t) = T_NEAR_FIELD;
        elseif ( r - 0.5 * L * Ft ) < 1e-14 % [r = 0.5*L*Ft] TO ACCOUNT FOR FLOATING POINT CALCULATION i.e. 1-1 = 0. 1e-14 = tolerance
                    T_MAX(x,t) = T_NEAR_FIELD;     
        else                          % FAR-FIELD
            T_MAX(x,t) = T_ROOM + ( 5.38 / HEIGHT ) * ( L * Ft * WIDTH * QII / r )^(2/3);
                if T_MAX(x,t) > T_NEAR_FIELD 
                    T_MAX(x,t) = T_NEAR_FIELD;
                end
        end
        %= END ==========================================================
    end
end

% % ======================================================================= 
% % ============================================================== FIGURES
% % =======================================================================

t_plot = dsearchn(tt',(ti*60)');
x_plot = dsearchn(xx',xi');

LTU = max([numel(t_plot),numel(x_plot),7]);
colors = colormap(hsv(LTU));

% === VARIATION OF GAS TEMPERATURE WITH TIME ===

figure(1)
hold on
for i = 1:numel(x_plot)
    h1(i) = plot( tt/60,T_MAX( x_plot(i) , : ),...
                 'Color', colors(i, :),...
                 'DisplayName',sprintf('%d m',xi(i)) );
end
legend(h1)
axis([0 TLIMIT/60 0 1400])
grid off;
xlabel('Time (min)','Fontsize',10)
ylabel('Temperature (deg.C)','Fontsize',10)

% === VARIATION OF GAS TEMPERATURE ALONG THE FIRE PATH ===

figure(2)
hold on
for i = 1:numel(t_plot)
    h2(i) = plot( xx,T_MAX( : , t_plot(i) ),...
                 'Color', colors(i, :),...
                 'DisplayName',sprintf('%d min',ti(i)) );
end
legend(h2)
axis([0 L 0 1400])
grid off;
xlabel('Location along the fire path (m)','Fontsize',10)
ylabel('Temperature (deg.C)','Fontsize',10)