diff --git a/Source/CFAST/cfast_data.f90 b/Source/CFAST/cfast_data.f90 index 2d53152da..f3127d877 100644 --- a/Source/CFAST/cfast_data.f90 +++ b/Source/CFAST/cfast_data.f90 @@ -12,6 +12,8 @@ module cenviro real(eb), parameter :: cp = 1012.0_eb real(eb), parameter :: gamma = 1.40_eb real(eb), parameter :: rgas = (gamma-1.0_eb)/gamma*cp + real(eb), parameter :: mw_n2 = 28.01340e-3_eb, mw_o2 = 31.99880e-3_eb, mw_co2 = 44.00950e-3_eb, mw_h2o = 18.01528e-3_eb, & + mw_c = 12.01070e-3_eb, mw_co = 28.01010e-3_eb, mw_hcn = 27.02534e-3_eb, mw_hcl = 36.46094e-3_eb end module cenviro diff --git a/Source/CFAST/fire.f90 b/Source/CFAST/fire.f90 index ae78ef1a4..190e36552 100644 --- a/Source/CFAST/fire.f90 +++ b/Source/CFAST/fire.f90 @@ -8,7 +8,7 @@ module fire_routines use cfast_types, only: detector_type, fire_type, room_type, target_type, vent_type - use cenviro, only: cp + use cenviro, only: cp, mw_n2, mw_o2, mw_co2, mw_h2o, mw_c, mw_co, mw_hcn, mw_hcl use cparams, only: u, l, m, q, ts, fuel, n2, o2, o2f, mxfires, mxrooms, fuel_moles, fuel_q, fuel_n2, fuel_o2, fuel_co2, & fuel_co, fuel_hcn, fuel_hcl, hcl, fuel_h2o, h2o, fuel_soot, co2, co, hcn, soot, soot_flaming, soot_smolder, ct, & mx_hsep, t_max, ns_mass, flaming, smoldering, trigger_by_time, trigger_by_temp, trigger_by_flux, idx_tempf_trg, & @@ -222,7 +222,7 @@ subroutine interpolate_pyrolysis (ifire,tsec,iroom,mdot_t,area_t,height_t,qdot_t if (fireptr%n_hcn>0) then call interp(fireptr%t_hcn,fireptr%y_hcn,fireptr%n_hcn,xxtime,1,y_hcn) else - y_HCN = fireptr%n_N*0.027028_eb/fireptr%molar_mass + y_HCN = fireptr%n_N*mw_hcn/fireptr%molar_mass end if call interp(fireptr%t_trace,fireptr%y_trace,fireptr%n_trace,xxtime,1,y_trace) call interp(fireptr%t_area,fireptr%area,fireptr%n_area,xxtime,1,area_t) @@ -569,11 +569,11 @@ subroutine chemistry (pyrolysis_rate, molar_mass, entrainment_rate, hoc, y_soot, ! now do the chemistry balance with supplied inputs. nu_soot = molar_mass/0.01201_eb*y_soot - nu_soot_flaming = molar_mass/0.01201_eb*y_soot_flaming - nu_soot_smolder = molar_mass/0.01201_eb*y_soot_smolder - nu_hcn = min(n_N,molar_mass/0.027028_eb*y_hcn) + nu_soot_flaming = molar_mass/mw_c*y_soot_flaming + nu_soot_smolder = molar_mass/mw_c*y_soot_smolder + nu_hcn = min(n_N,molar_mass/mw_hcn*y_hcn) nu_hcl = n_Cl - nu_co = molar_mass/0.02801_eb*y_co + nu_co = molar_mass/mw_co*y_co nu_h2o = (n_H - nu_hcl - nu_hcn)/2.0_eb nu_co2 = n_C - nu_co - nu_hcn - nu_soot nu_o2 = nu_co2 + (nu_h2o + nu_co - n_O)/2.0_eb @@ -582,16 +582,16 @@ subroutine chemistry (pyrolysis_rate, molar_mass, entrainment_rate, hoc, y_soot, ! chemistry balance is molar-based so convert back to mass rates. fuel and o2 are consumed, ! so negative. Others are produced, so positive net_fuel = -pyrolysis_rate_constrained - net_o2 = -pyrolysis_rate_constrained*nu_o2*0.032_eb/molar_mass - net_co2 = pyrolysis_rate_constrained*nu_co2*0.04401_eb/molar_mass - net_co = pyrolysis_rate_constrained*nu_co*0.02801_eb/molar_mass - net_h2o = pyrolysis_rate_constrained*nu_h2o*0.018016_eb/molar_mass - net_hcl = pyrolysis_rate_constrained*nu_hcl*0.036458_eb/molar_mass - net_hcn = pyrolysis_rate_constrained*nu_hcn*0.027028_eb/molar_mass - net_soot = pyrolysis_rate_constrained*nu_soot*0.01201_eb/molar_mass - net_soot_flaming = pyrolysis_rate_constrained*nu_soot_flaming*0.01201_eb/molar_mass - net_soot_smolder = pyrolysis_rate_constrained*nu_soot_smolder*0.01201_eb/molar_mass - net_n2 = pyrolysis_rate_constrained*nu_n2*0.02802_eb/molar_mass + net_o2 = -pyrolysis_rate_constrained*nu_o2*mw_o2/molar_mass + net_co2 = pyrolysis_rate_constrained*nu_co2*mw_co2/molar_mass + net_co = pyrolysis_rate_constrained*nu_co*mw_co/molar_mass + net_h2o = pyrolysis_rate_constrained*nu_h2o*mw_h2o/molar_mass + net_hcl = pyrolysis_rate_constrained*nu_hcl*mw_hcl/molar_mass + net_hcn = pyrolysis_rate_constrained*nu_hcn*mw_hcn/molar_mass + net_soot = pyrolysis_rate_constrained*nu_soot*mw_c/molar_mass + net_soot_flaming = pyrolysis_rate_constrained*nu_soot_flaming*mw_c/molar_mass + net_soot_smolder = pyrolysis_rate_constrained*nu_soot_smolder*mw_c/molar_mass + net_n2 = pyrolysis_rate_constrained*nu_n2*mw_n2/molar_mass net_ct = 0.0_eb if (-net_o2 > o2_available) then @@ -800,13 +800,13 @@ subroutine vent_jet_fire (ito, ifrom, netfuel, entrainment_rate, hrr, species_ma y_smolder_soot = 0.0_eb y_co = room2ptr%species_mass(u,fuel_co)/room2ptr%species_mass(u,fuel) y_hcn = room2ptr%species_mass(u,fuel_hcn)/room2ptr%species_mass(u,fuel) - nu_soot = y_soot*xxmolar_mass/0.01201_eb - nu_co = y_co*xxmolar_mass/0.02801_eb - nu_hcn = y_hcn*xxmolar_mass/0.027028_eb - nu_hcl = room2ptr%species_mass(u,fuel_hcl)/room2ptr%species_mass(u,fuel)*xxmolar_mass/0.036458_eb - nu_co2 = room2ptr%species_mass(u,fuel_co2)/room2ptr%species_mass(u,fuel)*xxmolar_mass/0.04401_eb - nu_h2o = room2ptr%species_mass(u,fuel_h2o)/room2ptr%species_mass(u,fuel)*xxmolar_mass/0.018016_eb - nu_o2 = room2ptr%species_mass(u,fuel_o2)/room2ptr%species_mass(u,fuel)*xxmolar_mass/0.032_eb + nu_soot = y_soot*xxmolar_mass/mw_c + nu_co = y_co*xxmolar_mass/mw_co + nu_hcn = y_hcn*xxmolar_mass/mw_hcn + nu_hcl = room2ptr%species_mass(u,fuel_hcl)/room2ptr%species_mass(u,fuel)*xxmolar_mass/mw_hcl + nu_co2 = room2ptr%species_mass(u,fuel_co2)/room2ptr%species_mass(u,fuel)*xxmolar_mass/mw_co2 + nu_h2o = room2ptr%species_mass(u,fuel_h2o)/room2ptr%species_mass(u,fuel)*xxmolar_mass/mw_h2o + nu_o2 = room2ptr%species_mass(u,fuel_o2)/room2ptr%species_mass(u,fuel)*xxmolar_mass/mw_o2 n_n = nu_hcn n_cl = nu_hcl @@ -1136,11 +1136,10 @@ subroutine update_species (deltt) real(eb), intent(in) :: deltt ! molar_masses of the species - real(eb), parameter :: molar_mass(ns_mass) = & - (/0.02802_eb,0.032_eb,0.04401_eb,0.02801_eb,0.027028_eb,0.036458_eb,0.01201_eb,0.018016_eb,0.01201_eb/) + real(eb) :: molar_mass(ns_mass) = (/mw_n2, mw_o2, mw_co2, mw_co, mw_hcn, mw_hcl, mw_c, mw_h2o, mw_c/) ! reciprocal of avagadro's number (so you can't have less than an atom of a species) - real(eb), parameter :: avagad = 1.0_eb/6.022e23_eb + real(eb), parameter :: avagad = 1.0_eb/6.02214076e23_eb real(eb) :: air_moles(2), v(2) integer i, layer, lsp diff --git a/Source/CFAST/initialization.f90 b/Source/CFAST/initialization.f90 index 4162e6b6f..e8a45bf89 100644 --- a/Source/CFAST/initialization.f90 +++ b/Source/CFAST/initialization.f90 @@ -212,7 +212,7 @@ end subroutine initialize_ambient subroutine initialize_leakage integer iroom, counter - real(eb), parameter :: rho_leak = 352.981915_eb/298.15_eb + real(eb), parameter :: rho_leak = default_pressure/rgas/298.15 real(eb) :: area, area2, leak_area type(room_type), pointer :: roomptr type(vent_type), pointer :: ventptr @@ -744,7 +744,7 @@ subroutine initialize_species ! tables in the handbook of physics and chemistry. xt = interior_ambient_temperature xtemp = 23.2_eb - 3.816e3_eb/(xt-46.0_eb) - xh2o = exp(xtemp)/101325.0_eb*(18.016_eb/28.584_eb) + xh2o = exp(xtemp)/default_pressure*(18.016_eb/28.584_eb) initial_mass_fraction(h2o) = relative_humidity*xh2o ! normalize the compartment atmosphere diff --git a/Source/CFAST/output.f90 b/Source/CFAST/output.f90 index eef3e6376..bc61e5733 100644 --- a/Source/CFAST/output.f90 +++ b/Source/CFAST/output.f90 @@ -452,7 +452,7 @@ subroutine results_vent_flows 'Lower Layer Upper Layer Lower Layer',/, & 'Vent From/Bottom To/Top Fraction Inflow Outflow ', & 'Inflow Outflow Inflow Outflow Inflow Outflow',/,153('-')) -5070 format (a4,3x,a12,3x,a12,15x,a) +5070 format (a4,3x,a12,3x,a12,15x,a104) end subroutine results_vent_flows @@ -1161,7 +1161,7 @@ subroutine flwout (outbuf,flow1,flow2,flow3,flow4,flow5,flow6,flow7,flow8) real(eb), intent(in) :: flow1, flow2, flow3, flow4, flow5, flow6, flow7, flow8 character(len=*), intent(out) :: outbuf - real :: flow(8), x1000,x100,x10,x1,x01 + real :: flow(8) integer :: i @@ -1174,34 +1174,14 @@ subroutine flwout (outbuf,flow1,flow2,flow3,flow4,flow5,flow6,flow7,flow8) flow(6) = real(flow6) flow(7) = real(flow7) flow(8) = real(flow8) - x1000 = 1000.0 - x100 = 100.0 - x10 = 10.0 - x1 = 1.0 - x01 = 0.1 do i = 1, 8 - if (flow(i)>=x1000) then - write (outbuf(13*(i-1)+1:13*i),5000) flow(i) - else if (flow(i)>=x100) then - write (outbuf(13*(i-1)+1:13*i),5010) flow(i) - else if (flow(i)>=x10) then - write (outbuf(13*(i-1)+1:13*i),5020) flow(i) - else if (flow(i)>=x1) then - write (outbuf(13*(i-1)+1:13*i),5030) flow(i) - else if (flow(i)>=x01) then - write (outbuf(13*(i-1)+1:13*i),5040) flow(i) + if (flow(i)<=atol) then + outbuf(13*(i-1)+1:13*i) = ' ' else - write (outbuf(13*(i-1)+1:13*i),5000) flow(i) - end if - if (flow(i)<=atol) outbuf(13*(i-1)+1:13*i) = ' ' - if (flow(i).ne.0.0_eb) write (outbuf(13*(i-1)+1:13*i),5050) flow(i) + write (outbuf(13*(i-1)+1:13*i),5050) flow(i) + endif end do -5000 format (2x,1pg11.3) -5010 format (f6.0,7x) -5020 format (f7.1,6x) -5030 format (f8.2,5x) -5040 format (f9.3,4x) 5050 format (2x,e11.4) end subroutine flwout diff --git a/Source/CFAST/radiation.f90 b/Source/CFAST/radiation.f90 index 0eac57dca..e7404cfdc 100644 --- a/Source/CFAST/radiation.f90 +++ b/Source/CFAST/radiation.f90 @@ -767,6 +767,7 @@ end subroutine cross_product ! --------------------------- absorb ------------------------------------------- real(eb) function absorb(cmpt, layer) + use cenviro, only: mw_co2, mw_h2o ! calculates absorbance, due to gases (co2 and h2o) and soot, for the specified compartment and layer. @@ -820,9 +821,9 @@ real(eb) function absorb(cmpt, layer) ! declare module data - ! physical constants [mw in kg/mol; rg in m^3-atm/mol-kelvin] + ! physical constants [rg in m^3-atm/mol-kelvin] - real(eb), parameter :: mwco2 = 44.0088e-3_eb, mwh2o = 18.0153e-3_eb, rg = 82.0562e-6_eb, k = 1195.5_eb, rhos = 1800.0_eb + real(eb), parameter :: rg = 82.0562e-6_eb, k = 1195.5_eb, rhos = 1800.0_eb ! log(t) data for co2 [t in k] @@ -910,7 +911,7 @@ real(eb) function absorb(cmpt, layer) ag = 0.0_eb ! absorbance for co2 - ng = roomptr%species_mass(layer,co2)/mwco2 + ng = roomptr%species_mass(layer,co2)/mw_co2 plg = ng*rtv*l if (plg>0.0_eb) then cplg = log10(plg) @@ -922,7 +923,7 @@ real(eb) function absorb(cmpt, layer) end if ! absorbance for h2o - ng = roomptr%species_mass(layer,h2o)/mwh2o + ng = roomptr%species_mass(layer,h2o)/mw_h2o plg = ng*rtv*l if (plg>0.0_eb) then cplg = log10(plg) diff --git a/Source/CFAST/solve.f90 b/Source/CFAST/solve.f90 index 47e4bad55..756e79e11 100644 --- a/Source/CFAST/solve.f90 +++ b/Source/CFAST/solve.f90 @@ -7,6 +7,7 @@ module solve_routines use conduction_routines, only: conduction use convection_routines, only: convection use debug_routines, only: output_spreadsheet_residuals + use defaults, only: default_pressure use exit_routines, only: cfastexit, delete_output_files use fire_routines, only: fire, vent_jets, integrate_mass, update_species, collect_fire_data_for_smokeview, update_fire_ignition use isosurface, only: output_isodata @@ -25,7 +26,7 @@ module solve_routines use cfast_types, only: fire_type, room_type, target_type, vent_type - use cenviro, only: odevara, odevarb, odevarc, constvar, cp, rgas, gamma + use cenviro, only: odevara, odevarb, odevarc, constvar, cp, rgas, gamma, mw_co2, mw_h2o use cparams, only: u, l, m, q, mxrooms, mxtarg, mxnode, mxbranch, mxdiscon, maxeq, ns, check_state, set_state, update_state, & nwal, ns_mass, vminfrac, n2, o2, co2, co, h2o, w_from_room, w_to_room, w_from_wall, w_to_wall, w_boundary_condition, & radiation_fix @@ -1112,7 +1113,7 @@ subroutine update_data (y_vector,iflag) ! tables in the handbook of physics and chemistry. xt = exterior_ambient_temperature xtemp = 23.2_eb - 3.816e3_eb/(xt-46.0_eb) - xh2o = exp(xtemp)/101325.0_eb*(18.016_eb/28.584_eb) + xh2o = exp(xtemp)/default_pressure*(mw_h2o*1000._eb/28.584_eb) roomptr%species_fraction(u,h2o) = relative_humidity*xh2o roomptr%species_fraction(l,h2o) = relative_humidity*xh2o @@ -1350,7 +1351,7 @@ subroutine update_data (y_vector,iflag) roomptr => roominfo(iroom) isof = isof + 1 if (radi_verification_flag) then - if (partial_pressure_co2+partial_pressure_h2o == 101325._eb) then ! Only H2O and CO2 exist + if (partial_pressure_co2+partial_pressure_h2o == default_pressure) then ! Only H2O and CO2 exist !if (lsp == 1 .or. lsp == 2) then ! N2 and O2 will have to be zero if (lsp == n2 .or. lsp == o2) then ! N2 and O2 will have to be zero ppgas = 0._eb @@ -1358,11 +1359,11 @@ subroutine update_data (y_vector,iflag) else if (lsp == co2) then ! partical pressure = species_mass/its molecular weight*ideal gas constant*gas temperature/volume of the medium roomptr%species_mass(u,lsp) = & - partial_pressure_co2/101325._eb*44.0088e-3_eb/82.0562e-6_eb/roomptr%temp(u)*roomptr%volume(u) + partial_pressure_co2/default_pressure*mw_co2/82.0562e-6_eb/roomptr%temp(u)*roomptr%volume(u) !else if (lsp == 8) then else if (lsp == h2o) then roomptr%species_mass(u,lsp) = & - partial_pressure_h2o/101325._eb*18.0153e-3_eb/82.0562e-6_eb/roomptr%temp(u)*roomptr%volume(u) + partial_pressure_h2o/default_pressure*mw_h2o/82.0562e-6_eb/roomptr%temp(u)*roomptr%volume(u) else ppgas = 0._eb end if @@ -1378,18 +1379,18 @@ subroutine update_data (y_vector,iflag) isof = isof + 1 if (radi_verification_flag) then - if (partial_pressure_co2+partial_pressure_h2o == 101325._eb) then + if (partial_pressure_co2+partial_pressure_h2o == default_pressure) then !if (lsp == 1 .or. lsp == 2) then if (lsp == n2 .or. lsp == o2) then ppgas = 0._eb !else if (lsp == 3) then else if (lsp == co2) then roomptr%species_mass(l,lsp) = & - partial_pressure_co2/101325._eb*44.0088e-3_eb/82.0562e-6_eb/roomptr%temp(l)*roomptr%volume(l) + partial_pressure_co2/default_pressure*mw_co2/82.0562e-6_eb/roomptr%temp(l)*roomptr%volume(l) !else if (lsp == 8) then else if (lsp == h2o) then roomptr%species_mass(l,lsp) = & - partial_pressure_h2o/101325._eb*18.0153e-3_eb/82.0562e-6_eb/roomptr%temp(l)*roomptr%volume(l) + partial_pressure_h2o/default_pressure*mw_h2o/82.0562e-6_eb/roomptr%temp(l)*roomptr%volume(l) else ppgas = 0._eb end if