pb_mat

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

clc; close all;
opengl('save','software');
% This program calculates the pan fire burn rate given all the geometric and
% operational parameters. The properties of the fuel and the pan used is
% show below.

% 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 = 0.202; % 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.04; % m, pan depth
g = 9.81; % m/s2, acceleration due to gravity
 Kwall(1) =0.016; Kwall(2) = 0.032;  Kwall(3) = 0.06;
mfuel = [];condflx = [];convflx=[];radflx=[]; Tpp=[];TTbot=[];rrhofurdot=[];
%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');


for j = 1:3

% Initial conditions and others
dpan = 0.202; % Pan diameter in m
kw=Kwall(j);
hfu = 0.013; % Fuel depth in m
T0 = 300; % initial temperature in K
m0 = rhofu*0.7854*dpan^2*hfu; % initial mass of fuel in kg
deltal = 0.002; % top fuel thickness of uniform temperature in m
evf0 = 0.018; % emissivity * view factor
q_ffb = 0;  % flux from flame to free board at time t=0.
deltat=1;    %time step used in calculation, sec
tmax = 800;   % maximum time
% Various heat transfer coefficients
hfb = hpan - hfu ;% m, Initial free board is equal to difference b/w the pan depth and fuel depth used in m.
hgcv0 = 0.0045; %  Initial gas phase convective heat transfer coefficient in kW/m2 K.
hgwfu0 = 0.0045; % Pan wall to the fuel heat transfer coefficient in kW/m2 K.

% Correlation based on the Non-dimesional number to calculate mean mass
% bunr rate. This correlation is deduced by considering the geometrical and thermo-chemical properties of pan and the fuel used.

P1=((kw/(hpan*hgcv0))*(hfu/hfb))^(1/4); %P1 accounts for conductive,convective heat feedback to fuel and fuel depth effects.
P2=1 - exp(-0.25*(dpan/0.21)^(1.5)/P1); % P2 accounts for the pan diameter effect.
P3 = (((Tbfu-T0)/(Tbfu-300))*(300/Tbfu))^(-0.5); % P3 accounts for the initial temperature effect of fuel
hm=kw/(hpan);
Mpc1=(hm/hgcv0)*(4*tw/dpan)*(hfu/hpan*Lfu/cpfu/(Tbfu-T0))^0.25*(1 - 0.3*(hpan/dpan)^0.125);
Mpc = P1*P3*(1.5+8.5*P2);
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.45; 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');

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

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

Mpc1=   13.8 
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;


tM2 = xlsread(hsme, 'A4:A51'); mM2 = xlsread(hsme, 'M4:M51');
TpM2 = xlsread(hsme,'J4:J51');
tM3 = xlsread(hsme, 'A4:A43'); mM3 = xlsread(hsme, 'Q4:Q43');
TpM3 = xlsread(hsme,'N4:N43');
tM4 = xlsread(hsme, 'A4:A35'); mM4 = xlsread(hsme, 'AF4:AF35');
TpM4 = xlsread(hsme,'AE4:AE35');

%Variable declaration
Ts = T0*ones(numel(t),1); Tp = Ts; Twfu = Ts; Tbot = Ts; Tf = Ts;
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; BTO = Ts;
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;tb=T0;
hgwfu = zeros(numel(t),1);  % initializing hgfu
mfu = zeros(numel(t),1); mfu(1) = m0;  qcond=mfu; qconv=mfu; qrad=mfu;       % initializing mfu
 m = m0; rhofurdot=mfu;
condscle = sqrt(alphafu*(t+0.01));  %  m, conduction thickness
etal = deltal./condscle;  % dimensionless top fuel thickness
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);


        if((Mpc1 > 1) && (Mpc1 <= 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; % Convective heat transfer coefficient
        evf(i) = 0.20*(1 - exp(-0.00045*(1 - exp(-0.4*(dpan/0.2)^3))*(rhofurdot(i-1)*dpan/0.000018))); % product of constant part of emissivity and view factor
        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]; Tpp =[Tpp Tp];TTbot =[TTbot Tbot]; rrhofurdot=[rrhofurdot rhofurdot];

0.202   0.013   0.027   0.003    0.016   300.0      6.4   5.086     0.312  0.033    15.9   511.0 
mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = 18.52, 17.30, 17.91  6.82

0.202   0.013   0.027   0.003    0.032   300.0      7.2   6.057     0.262  0.047    46.1   369.0 
mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = 20.90, 23.96, 22.43 13.62

0.202   0.013   0.027   0.003    0.060   300.0      8.1   7.043     0.224  0.064    84.3   293.0 
mean mass flux (g/m2s), from corrln, from code, mean, deviation-p = 23.44, 30.17, 26.80 25.12

end %on j



     figure(1); clf % Conduction, convection and radiation fluxes for 10, 20 and 13 mm fuel depths.
    plot(t, condflx(:,1),'--g',t, condflx(:,2),':b',t, condflx(:,3),'--k', 'LineWidth', 2, 'MarkerSize', 2);
    set(gca,'FontSize', 12)
    xlabel('t, s'); ylabel('Fluxes, kW/m^2');
    legend('condflx-SS','condflx-MS','condflx-AL','location', 'northeast');

    figure(2); clf % Experimental vs predicted pan tip temperatures.
   plot(t, Tpp(:,1),'--c',t, Tpp(:,2),'--g',t, Tpp(:,3),'--b',tM2,TpM2,'+',tM3,TpM3,'>',tM4,TpM4,'o','markersiz',2,'linewidth',3);
   set(gca,'ylim',[300 500],'FontSize', 12)
   xlabel('t, s'); ylabel('Pan wall tip temperature, K');
   title('Pan tip Temperature');
   legend('pred-SS','pred-MS','Pred-AL','expt-SS','expt-MS','expt-AL','location', 'northeast');


    figure(4); clf % mass loss expt vs predicted
    set(gca,'FontSize', 12)
    plot(t, mfuel(:,1), 'r',t,mfuel(:,2),'g',t,mfuel(:,3),'k',tM2,mM2,'*',tM3,mM3,'<',tM4,mM4,'o', 'markersiz',2,'linewidth',3);
    set(gca,'ylim',[0 0.3],'FontSize', 12,'Fontweight','bold')
    legend('10mm-pred','13mm-pred','20mm-pred','10mm-Expt','13mm-Expt','20mm-Expt', 'location', 'northeast');
    xlabel('t, s'); ylabel('mass kg (pred, expt)')

Contents

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