| blob | 532d1a78b35dc12105c0bd29afa97def08b50c78 |
1 function ros = fire_ros_balbi(fuel,speed,tanphi,fmc_g)\r
2 % in\r
3 % fuel fuel description structure\r
4 % speed wind speed\r
5 % tanphi slope\r
6 % fmc_g optional, overrides fuelmc_g from the fuel description\r
7 % out\r
8 % ros rate of spread\r
9 \r
10 % given fuel params\r
11 \r
12 windrf=fuel.windrf; % WIND REDUCTION FACTOR\r
13 fgi=fuel.fgi; % INITIAL TOTAL MASS OF SURFACE FUEL (KG/M**2)\r
14 fueldepthm=fuel.fueldepthm; % FUEL DEPTH (M)\r
15 savr=fuel.savr; % FUEL PARTICLE SURFACE-AREA-TO-VOLUME RATIO, 1/FT\r
16 fuelmce=fuel.fuelmce; % MOISTURE CONTENT OF EXTINCTION\r
17 fueldens=fuel.fueldens; % OVENDRY PARTICLE DENSITY, LB/FT^3\r
18 st=fuel.st; % FUEL PARTICLE TOTAL MINERAL CONTENT\r
19 se=fuel.se; % FUEL PARTICLE EFFECTIVE MINERAL CONTENT\r
20 weight=fuel.weight; % WEIGHTING PARAMETER THAT DETERMINES THE SLOPE OF THE MASS LOSS CURVE\r
21 fci_d=fuel.fci_d; % INITIAL DRY MASS OF CANOPY FUEL\r
22 fct=fuel.fct; % BURN OUT TIME FOR CANOPY FUEL, AFTER DRY (S)\r
23 ichap=fuel.ichap; % 1 if chaparral, 0 if not\r
24 fci=fuel.fci; % INITIAL TOTAL MASS OF CANOPY FUEL\r
25 fcbr=fuel.fcbr; % FUEL CANOPY BURN RATE (KG/M**2/S)\r
26 hfgl=fuel.hfgl; % SURFACE FIRE HEAT FLUX THRESHOLD TO IGNITE CANOPY (W/m^2)\r
27 cmbcnst=fuel.cmbcnst; % JOULES PER KG OF DRY FUEL\r
28 fuelheat=fuel.fuelheat; % FUEL PARTICLE LOW HEAT CONTENT, BTU/LB\r
29 fuelmc_g=fuel.fuelmc_g; % FUEL PARTICLE (SURFACE) MOISTURE CONTENT, jm: 1 by weight?\r
30 fuelmc_c=fuel.fuelmc_c; % FUEL PARTICLE (CANOPY) MOISTURE CONTENT, 1\r
31 \r
32 if exist('fmc_g','var') % override moisture content by given\r
33 fuelmc_g = fmc_g;\r
34 end\r
35 \r
36 % Universal constants\r
37 s = 9; % Stoichiometric constant - Balbi 2009\r
38 Chi_0 = 0.3; % Thin Flame Radiant Fraction - ?Balbi 2009?\r
39 a = 0.05; % Constant from Balbi 2008\r
40 A_0 = 2.25; % Constant from Balbi 2008\r
41 eps = 0.2; % Pastor 2002\r
42 B = 5.67e-8; % Stefan-Boltzman \r
43 Deltah_v = 2.257e6; % Water Evap Enthalpy [J/kg]\r
44 C_p = 2e3; % Calorific Capacity [J/kg] - Balbi 2009\r
45 \r
46 % Fuel Constants\r
47 m = fuelmc_g*100; % FUEL PARTICLE MOISTURE CONTENT [%]\r
48 rho_v = fueldens*16.0185; % FUEL Particle Density [Kg/m^3]\r
49 DeltaH = cmbcnst; % Combustion Enthalpy [J/kg]\r
50 sigma = fgi;\r
51 e_delta = (savr/0.3048)*sigma/(4*rho_v);\r
52 rho = 0.25; % Gas Flame Density [kg/m^3]\r
53 \r
54 % Calculate R_0\r
55 ii = 1;\r
56 % USE Roethermel Model for no slope, no wind ros\r
57 if ii == 1\r
58 bmst = fuelmc_g/(1+fuelmc_g); % jm: 1 \r
59 fuelloadm= (1.-bmst) * fgi; % fuelload without moisture\r
60 % jm: 1.-bmst = 1/(1+fuelmc_g) so fgi includes moisture? \r
61 fuelload = fuelloadm * (.3048)^2 * 2.205; % to lb/ft^2\r
62 fueldepth= fueldepthm/0.3048; % to ft\r
63 betafl = fuelload/(fueldepth * fueldens);% packing ratio jm: lb/ft^2/(ft * lb*ft^3) = 1\r
64 betaop = 3.348 * savr^(-0.8189); % optimum packing ratio jm: units?? \r
65 qig = 250. + 1116.*fuelmc_g; % heat of preignition, btu/lb\r
66 epsilon = exp(-138./savr ); % effective heating number\r
67 rhob = fuelload/fueldepth; % ovendry bulk density, lb/ft^3\r
68 rtemp2 = savr^1.5;\r
69 gammax = rtemp2/(495. + 0.0594*rtemp2); % maximum rxn vel, 1/min\r
70 ar = 1./(4.774 * savr^0.1 - 7.27); % coef for optimum rxn vel\r
71 ratio = betafl/betaop; \r
72 gamma = gammax*(ratio^ar)*exp(ar*(1.-ratio)); % optimum rxn vel, 1/min\r
73 wn = fuelload/(1 + st); % net fuel loading, lb/ft^2\r
74 rtemp1 = fuelmc_g/fuelmce;\r
75 etam = 1.-2.59*rtemp1 +5.11*rtemp1^2 -3.52*rtemp1^3; % moist damp coef\r
76 etas = 0.174* se^(-0.19); % mineral damping coef\r
77 ir = gamma * wn * fuelheat * etam * etas; % rxn intensity,btu/ft^2 min\r
78 xifr = exp( (0.792 + 0.681*savr^0.5)...\r
79 * (betafl+0.1)) /(192. + 0.2595*savr); % propagating flux ratio\r
80 roth_r_0 = ir*xifr/(rhob * epsilon *qig); % default spread rate in ft/min\r
81 R_0 = roth_r_0 * .00508; % default spread rate in m/s\r
82 else\r
83 % Calculate R_00 from flame temp\r
84 %R_00 = eps*B*a*T^4/(2*Deltah_v);\r
85 % Use R_00 = 0.05 from Balbi 2009\r
86 R_00 = .05;\r
87 R_0 = fueldepthm*R_00/(sigma*(1+a*m));\r
88 end\r
89 \r
90 % Calculate u_0\r
91 u_00 = 2*(s+1)/rho;\r
92 tau = fct/30; % Burnout time from Rothermel fuel model/3 [s]\r
93 u_0 = u_00*sigma/tau;\r
94 \r
95 % Calculate A\r
96 nu = min(e_delta,1);\r
97 A = nu*A_0/(1+a*m);\r
98 \r
99 % Calculate flame tilt angle (gamma)\r
100 alpha = atan(tanphi); % Slope angle [rad]\r
101 psi = 0; % Angle between wind and flame; assume parallel currently\r
102 phi = 0; % Angle between flame front vector and slope vector\r
103 gamma = atan(tan(alpha)*cos(phi)+speed*cos(psi)/u_0);\r
104 \r
105 if(gamma <= 0)\r
106 R = R_0;\r
107 else\r
108 % Solve for ros quadratically, eq(13b) Balbi 2009\r
109 qa = cos(gamma)/(12*R_0);\r
110 qb = 1 - cos(gamma)/12 - A*(1 + sin(gamma) - cos(gamma));\r
111 qc = -1;\r
112 R1 = (-qb+sqrt(qb*qb-(4*qa*qc)))/(2*qa);\r
113 R2 = (-qb-sqrt(qb*qb-(4*qa*qc)))/(2*qa);\r
114 R = max(R1,R2);\r
115 end\r
116 % Calculate Radiant Fraction\r
117 Chi = Chi_0/(1 + R*cos(gamma)/(12*R_0));\r
118 %Chi = Chi_0;\r
119 \r
120 % Calculate Flame T\r
121 T_a = 300; % Air Temp [K]\r
122 T = T_a + (1-Chi)*DeltaH/((1+s)*C_p);\r
123 \r
124 % Calculate flame height (H)\r
125 gstar = 9.81*(T/T_a - 1);\r
126 H = u_0*u_0/(gstar*cos(alpha)*cos(alpha));\r
127 \r
128 % Calculate radiant flux (Q)\r
129 Q = (eps/2)*B*(T^4);\r
130 \r
131 ros = min(R,6);\r
132 \r
133 % Calculate depth\r
134 L = ros*tau;\r
135 end