pb_dia

Experimental data - Ours - data at 10 mm, 13mm and 20 mm at 30 C n-heptane

% Code to obtain the pan diameter effect of the mass burn rate of pool
% fires

clc; close all;
%opengl('save','software');
% This program calculates the pan fire bur rate given all the geometric and
% operational parameters
% Fuel = n-Heptane
rhofu = 680; % kg/m3, liquid fuel density
Tbfu = 369; % K, boiling point
cpfu = 2.1; % kJ/kg K, liquid fuel specific heat
Lfu = 322;  % kJ/kg, Latent heat of vaporization of the fuel
kfu = 0.00015; % kW/mK, fuel thermal conductivity
alphafu = kfu/rhofu/cpfu; % m2/s,  thermal  diffuisivity of the fuel
S = 15.2;
Tfm = 1200;

% Pan related data
% material MS
%dpan = 2.0; % m, pan diameter
tw = 0.003; % m, pan wall thickness
kw = 0.032; % kW/m K, wall thermal conductivity
rhow = 7800; % kg/m3, wall material density
cpw = 0.46;  % kJ/kg K, wall material specific heat
alphaw = kw/rhow/cpw; % m2/s, thermal diffusivity of wall material
rhocptw = rhow*cpw*tw;  % kJ/m2 K, wall thermal capacity
cpg = 1.0; % kJ/kg K,  gas phase specific heat
hpan = 0.06; % m, pan depth
g = 9.81; % m/s2, acceleration due to gravity
futh(1) = 0.01; futh(2) =0.02; futh(3) = 0.013;
condflx = [];convflx=[];radflx=[];
pan(1)=0.2; pan(2)=0.3; pan(3) =0.4; pan(4) =0.5;pan(5) =2; mfuel=[];

for j = 1:5

% Initial conditions and others
dpan = pan(j);
hfu = 0.020;
T0 = 300; % K, initial temperature
uwind = 0;  % wind speed, m/s
m0 = rhofu*0.7854*dpan^2*hfu; % kg, initial mass of fuel
deltal = 0.002; % m, top fuel thickness of uniform temperature
evf0 = 0.018; % emissivity * view factor
q_ffb = 0;  % flux from flame to free board
deltat=1;    % s, time step
tmax = 800;   % maximum time
% Various heat transfer coefficients
hwr = 0.0;  % water thickness
hfb = hpan - hfu - hwr;% m, free board
hgcv0 = 0.0045;% kW/m2 K, Initial gas phase convective heat transfer coefficient
hgwfu0 = 0.0045; % kW/m2 K, Wall to the fuel heat transfer coefficient
P1=((kw/(hpan*hgcv0))*(hfu/hfb))^(1/4);
P2=1 - exp(-0.25*(dpan/0.21)^(1.5)/P1);
P3 = (((Tbfu-T0)/(Tbfu-300))*(300/Tbfu))^(-0.5);
Mpc = P1*P3*(1.5+8.5*P2);
Mpc1=(kw/hpan/hgcv0)*(4*tw/dpan)*(hfu/hpan*Lfu/cpfu/(Tbfu-T0))^0.25*(1 - 0.3*(hpan/dpan)^0.125);
W = (Tbfu/T0)^0.5*(dpan*0.0045/kw)^0.25*(hfu*hfb/hpan/hpan)^0.1*dpan^0.25; % dpan in m only
W1 =  (1/222)*(kw/hfb/hgcv0)^0.5*(hpan/dpan)^0.5*(Tbfu/T0-1)^(-0.25);
Tpm = 320 +410*W;
if (dpan < 0.199) Tpm = 190 +1200 *W; end
cTp = (Tfm - Tpm)/(Tpm - T0);

C2 = 0.7;
if(Mpc1 >= 6.0 && Mpc1 < 10)C2 = 0.6 + 0.0575*(Mpc1 - 6); end
if(Mpc1 > 10) C2 = 0.83 - 0.3*(Mpc1^0.25 - 1.78)/Mpc1^0.25; end
if(Mpc1 > 6) Slope =0.1*Mpc1 - 0.26 - 0.0017*Mpc1^2; end
if(Mpc1 > 12) Slope = 0.85; end
fprintf('Mpc1=  %5.1f \n', Mpc1);

dTpdt0 = 2.8;
C3 = 2200*(W1-0.026);
if(dpan < 0.2 && W1 > 0.036) C3 = 1180*W1 - 37.3;  end
if (W1 < 0.036) C3 = 465*W1 - 3.75; end
mmflx1 = ((hgcv0*1000*(1200-369))/(Lfu))*(Mpc/4);

fprintf('dpan \t hfu  \t hfb  \t tw  \t kw \t T0 \t  Mpc \t cTp \t  \t W \t W1 \t\t C3  burn time  \n');
fprintf(' m \t     m \t \t   m  \t m \t \t kW/mK \t K  \t\t - \t \t - \t \t - \t -  \t\t  -  \t  s  \n');

t = 0:deltat:tmax; % t in seconds

Mpc1=    4.9 
dpan 	 hfu  	 hfb  	 tw  	 kw 	 T0 	  Mpc 	 cTp 	  	 W 	 W1 		 C3  burn time  
 m 	     m 	 	   m  	 m 	 	 kW/mK 	 K  		 - 	 	 - 	 	 - 	 -  		  -  	  s

Mpc1=    3.3 
dpan 	 hfu  	 hfb  	 tw  	 kw 	 T0 	  Mpc 	 cTp 	  	 W 	 W1 		 C3  burn time  
 m 	     m 	 	   m  	 m 	 	 kW/mK 	 K  		 - 	 	 - 	 	 - 	 -  		  -  	  s

Mpc1=    2.5 
dpan 	 hfu  	 hfb  	 tw  	 kw 	 T0 	  Mpc 	 cTp 	  	 W 	 W1 		 C3  burn time  
 m 	     m 	 	   m  	 m 	 	 kW/mK 	 K  		 - 	 	 - 	 	 - 	 -  		  -  	  s

Mpc1=    2.0 
dpan 	 hfu  	 hfb  	 tw  	 kw 	 T0 	  Mpc 	 cTp 	  	 W 	 W1 		 C3  burn time  
 m 	     m 	 	   m  	 m 	 	 kW/mK 	 K  		 - 	 	 - 	 	 - 	 -  		  -  	  s

Mpc1=    0.5 
dpan 	 hfu  	 hfb  	 tw  	 kw 	 T0 	  Mpc 	 cTp 	  	 W 	 W1 		 C3  burn time  
 m 	     m 	 	   m  	 m 	 	 kW/mK 	 K  		 - 	 	 - 	 	 - 	 -  		  -  	  s

Experimental data - Ours - data at 10 mm, 13mm and 20 mm at 30 C n-heptane

hsme = 'Mass loss data of pool fires'; sheet = 2;
tM1 = xlsread(hsme,'A4:A60'); mM1 = xlsread(hsme, 'Z4:Z60');
tM2 = xlsread(hsme, 'A4:A56'); mM2 = xlsread(hsme, 'AA4:AA56');
tM3 = xlsread(hsme, 'A4:A49'); mM3 = xlsread(hsme, 'AB4:AB49');
tM4 = xlsread(hsme, 'A4:A41'); mM4 = xlsread(hsme, 'AC4:AC41');
tM5 = xlsread(hsme, 'A4:A28'); mM5 = xlsread(hsme, 'AD4:AD28');


%Variable declaration
Ts = T0*ones(numel(t),1); Tp = Ts; Twfu = Ts; Tbot = Ts; Tf = Ts; Tp(1) = T0; tb=T0;
dTpdt = zeros(numel(t),1); rhofurdot = dTpdt; reg = zeros(numel(t),1); mfu =zeros(numel(t),1);
Term1 = zeros(numel(t),1); Tw38mm = Ts; evf = evf0;
qcond = dTpdt; qconv = dTpdt; qrad = dTpdt; qtot = qcond + qconv + qrad; ql = dTpdt; dTsdt = dTpdt;
etafu = dTpdt; Twb=T0*ones(numel(t),1); hgcv = hgcv0*ones(numel(t),1); y1 = zeros(numel(i), 1);
y2 = zeros(numel(i), 1); y3 = zeros(numel(i), 1); Tw1 = T0*ones(numel(t), 1); BTO = Tw1;
hgwfu = zeros(numel(t),1);  % initializing hgfu
mfu = zeros(numel(t),1);
mfu (1) = m0;        % initializing mfu
condscle = sqrt(alphafu*(t+0.01));  %  m, conduction thickness
etal = deltal./condscle;  % dimensionless top fuel thicknes5
etafu(1) = (hfu - reg(1))./condscle(1); ff = 3.*exp(-etafu.^0.66); m = m0;

for i = 2:numel(t)
    if (m/m0) > 0.001
        Tf(i) = T0 + (Tfm-300)*(1 - exp(-t(i)/5)) +15*sin(3.1415*t(i)/30+ t(i).*rand); % K  Flame temperature,
        dTpdt(i) = dTpdt0*((Tf(i)-Tp(i-1))/(Tf(i) - T0)- cTp*hpan/(hfb+reg(i-1))*(Tp(i-1)-T0)/((Tf(i) - T0)));
        Tp(i) = Tp(i-1) + dTpdt(i)*deltat; % K Variation of pan tip temperature with time
        BTO(i) = BTO(i-1)+deltat*16/dpan/rhow/cpw*kw/(hpan+dpan/4)*(Tp(i) - BTO(i-1));
        Tw1(i) = Tp(i) - (Tp(i) -BTO(i))*(hfb+reg(i-1))/(hpan+dpan/4);
        Twb(i) = Tp(i) - (Tp(i) -BTO(i))*(hpan)/(hpan+dpan/4);

        hgwfu(i) = hgwfu0*(1 + C3*reg(i-1)/hfu);
        if((Mpc1 > 1) && (Mpc <= 6.5)), hgwfu(i) = hgwfu0*(1 + C3*(reg(i-1)/hfu)); end % SS only
        if(Mpc1 < 1), hgwfu(i) = hgwfu0; end %  GL only
        if(Mpc1 > 6.5 && mfu(i-1)/m0 < C2), hgwfu(i)= Slope; end  % Al and MS

        qcond(i) = hgwfu(i)*(Twb(i-1)-Tbot(i-1)+ 4*(hfu - reg(i-1))/dpan*((Tw1(i)+Twb(i))/2 - (Ts(i-1)+Tbot(i-1))/2));
        hgcv(i) = hgcv0; % Effects of flame growth due to increase in burn rate and wind accounted
        evf(i) = 0.20*(1 - exp(-0.00045*(1 - exp(-0.4*(dpan/0.2)^3))*(rhofurdot(i-1)*dpan/0.000018)));
        qconv(i) = hgcv(i)*(Tf(i) - Ts(i-1));  %  kW/m2   Convective flux
        qconv(i) = hgcv(i)*(Tf(i) - Ts(i-1));  %  kW/m2   Convective flux
        qrad(i) = evf(i)*(56.78).*(Tf(i)/1000).^4; % kW/m2 Radiative flux
        qtot(i) = qcond(i) + qconv(i)+ qrad(i);  % kW/m2,  Total flux to the liquid
        Tbot(i) = T0 + (Ts(i-1) - T0)*ff(i-1);  % Bottom fuel temperature calculation
        if (Tbot(i)>=Tbfu), Tbot(i)=Tbfu; end
        ql(i) = 3*0.66*kfu.*(Ts(i-1)-Tbot(i)).*exp(-(etal(i)).^(0.66))/deltal^(0.34)/(condscle(i).^0.66); % heat flux into the liquid
        dTsdt(i) = (qtot(i) - ql(i))./(rhofu*cpfu*deltal);
        Ts(i) = Ts(i-1) + dTsdt(i)*deltat;
        if (Ts(i) >=Tbfu), Ts(i) = Tbfu; end
        rhofurdot(i) = qtot(i)./(Lfu + cpfu*(Ts(i)-Tbot(i)));
        reg(i) = reg(i-1) + rhofurdot(i)./rhofu*deltat;
        etafu(i) = (hfu - reg(i))./condscle(i); ff(i) = 3.*exp(-sqrt(etafu(i)*etafu(i)).^0.66);
        m = m - rhofurdot(i)*0.7854*dpan^2*deltat;
        mfu (i) = m;
        mfrac = mfu(i)/m0;
        if ((m > 0) && (mfrac < 0.012)) tb = t(i); end
        if(m >0) && (mfrac < 0.012) kk(j) =i;  end
        mmflx2 = hfu*rhofu*1000/tb;
        mmflxm = (mmflx1 + mmflx2)/2;
        pdvn = 100*sqrt((mmflx2 -mmflx1)^2)/mmflxm;
    end
end
fprintf('%5.3f  %6.3f  %6.3f  %6.3f   %6.3f  %6.1f   %6.1f  %6.3f    %6.3f %6.3f  %6.1f  %6.1f \n', dpan, hfu, hfb, tw, kw, T0, Mpc, cTp, W, W1, C3, tb)
fprintf('mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = %5.2f, %5.2f, %5.2f %5.2f\n', mmflx1, mmflx2, mmflxm, pdvn);
mfuel = [mfuel mfu];condflx = [condflx qcond];convflx = [convflx qconv];radflx = [radflx qrad];

0.200   0.020   0.040   0.003    0.032   300.0      6.7   6.079     0.261  0.048    47.3   566.0 
mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = 19.50, 24.03, 21.76 20.80

0.300   0.020   0.040   0.003    0.032   300.0      8.3   4.952     0.320  0.039    28.1   504.0 
mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = 24.23, 26.98, 25.61 10.76

0.400   0.020   0.040   0.003    0.032   300.0     10.1   4.247     0.370  0.034    11.9   458.0 
mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = 29.42, 29.69, 29.56  0.94

0.500   0.020   0.040   0.003    0.032   300.0     12.0   3.752     0.413  0.030    10.2   376.0 
mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = 34.80, 36.17, 35.49  3.86

2.000   0.020   0.040   0.003    0.032   300.0     28.9   1.508     0.826  0.015     3.2   206.0 
mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = 83.96, 66.02, 74.99 23.93

end


    figure(1); clf
    set(gca,'FontSize', 12)
    plot(t, mfuel(:,1), 'g',t, mfuel(:,2), 'c',t, mfuel(:,3), 'r',t, mfuel(:,4), 'k',tM1,mM1,'+',tM2,mM2,'*',tM3,mM3,'>',tM4,mM4,'<','markersiz',3,'linewidth',2);
    set(gca,'ylim',[0 3],'xlim',[0 600],'FontSize', 12)
    legend('0.2m-pred','0.3m-pred','0.4m-pred','0.5m-pred','0.2m-Expt','0.3m-Expt','0.4m-Expt','0.5m-Expt', 'location', 'northeast');
    title('Expt vs predicted mass loss in 0.2, 0.3, 0.4 & 0.5 m pans');
    xlabel('t, s'); ylabel('mass kg (pred, expt)');

    figure(2); clf
    plot(t, mfuel(:,5), 'g',tM5,mM5,'+','markersiz',5,'linewidth',4);
    set(gca,'ylim',[0 45],'xlim',[0 400],'FontSize', 12)
    legend('2m-pred','2m-Expt', 'location', 'northeast');
     title('Expt vs predicted mass loss in 2 m pan');
    xlabel('t, s'); ylabel('mass kg (pred, expt)');

    figure(3); clf % Conduction flux for 0.2, 0.3, 0.4, 0.5 and 2 m pans.
    plot(t, condflx(:,1),'--g',t, condflx(:,2),':b',t, condflx(:,3),'--k', t,condflx(:,4),'--r',t,condflx(:,5),'--c','LineWidth', 2, 'MarkerSize', 2);
    set(gca,'FontSize', 12)
    xlabel('t, s'); ylabel('Fluxes, kW/m^2');
    title('Conduction heat feed back in 0.2, 0.3, 0.4, 0.5 & 2 m pans');
    legend('condflx-0.2m','condflx-0.3m','condflx-0.4m','condflx-0.5m','condflx-2m', 'location', 'northwest');

    figure(4); clf % Convection flux for 0.2, 0.3, 0.4, 0.5 and 2 m pans.
    plot(t, convflx(:,1),'--g',t, convflx(:,2),':b',t, convflx(:,3),'--c', t,convflx(:,4),'--r',t,convflx(:,5),'--k','LineWidth', 2, 'MarkerSize', 2);
    set(gca,'FontSize', 12)
    xlabel('t, s'); ylabel('Fluxes, kW/m^2');
    title('Convection heat feed back in 0.2, 0.3, 0.4, 0.5 & 2 m pans');
    legend('convflx-0.2m','convflx-0.3m','convflx-0.4m','convflx-0.5m','convflx-2m', 'location', 'northeast');


    figure(5); clf % Radiation flux for 0.2, 0.3, 0.4, 0.5 and 2 m pans.
    plot(t, radflx(:,1),'--g',t, radflx(:,2),':b',t, radflx(:,3),'--c', t,radflx(:,4),'--r',t,radflx(:,5),'--k','LineWidth', 2, 'MarkerSize', 2);
    set(gca,'FontSize', 12)
    xlabel('t, s'); ylabel('Fluxes, kW/m^2');
    title('Radiation heat feed back in 0.2, 0.3, 0.4, 0.5 & 2 m pans');
    legend('radflx-0.2m','radflx-0.3m','radflx-0.4m','radflx-0.5m','radflx-2m', 'location', 'northeast');

Contents

Experimental data - Ours - data at 10 mm, 13mm and 20 mm at 30 C n-heptane