diff --git a/src/models/britter_mcquaid_puff.jl b/src/models/britter_mcquaid_puff.jl index ebe991e..8e9d15a 100644 --- a/src/models/britter_mcquaid_puff.jl +++ b/src/models/britter_mcquaid_puff.jl @@ -21,7 +21,7 @@ function puff(s::Scenario, ::Type{BritterMcQuaidPuff}, eqs=DefaultSet) end """ - puff(scenario::Scenario, ::BritterMcQuaidPuff[, equationset::EquationSet]) + puff(scenario::Scenario, ::BritterMcQuaidPuff[, equationset::EquationSet]; kwargs...) Returns the solution to the Britter-McQuaid instantaneous ground level release model for the given scenario. @@ -30,16 +30,19 @@ The `equationset` is used to calculate the windspeed at 10m, all other correlations are as per the Britter-McQuaid model. Unless otherwise specified a default power-law wind profile is used. +# Keyword Arguments +- `temp_correction=true`: Adds a correction for non-isothermal releases, per workbook section 5.5 + # References + Britter, Rex E. and J. McQuaid. 1988. *Workbook on the Dispersion of Dense Gases. HSE Contract Research Report No. 17/1988* + AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers """ -function puff(scenario::Scenario, ::BritterMcQuaidPuff, eqs=DefaultSet) +function puff(scenario::Scenario, ::BritterMcQuaidPuff, eqs=DefaultSet; temp_correction=true) Q = _release_flowrate(scenario) - ṁ = _mass_rate(scenario) + m = _release_mass(scenario) t = _duration(scenario) - ρⱼ = _release_density(scenario) + ρᵣ = _release_density(scenario) Tᵣ = _release_temperature(scenario) u₁₀ = windspeed(scenario, 10.0, eqs) @@ -50,15 +53,22 @@ function puff(scenario::Scenario, ::BritterMcQuaidPuff, eqs=DefaultSet) V₀ = Q*t # initial concentration - Qi = ṁ/ρⱼ - c₀ = min(Qi/Q,1.0) + Vi = m/ρᵣ + c₀ = min(Vi/V₀,1.0) # temperature correction - T′ = Tᵣ/Tₐ + if temp_correction + # non-isothermal correction, per workbook section 5.5 + T′ = Tᵣ/Tₐ + else + # no temperature correction + # T′ = 1.0 + T′ = 1.0 + end # relative density g = 9.80616 # gravity, m/s^2 - gₒ = g * ((ρⱼ - ρₐ)/ ρₐ) + gₒ = g * ((ρᵣ - ρₐ)/ ρₐ) # correlation parameter α = 0.5*log10( gₒ * cbrt(V₀) / u₁₀^2 ) @@ -95,7 +105,7 @@ function puff(scenario::Scenario, ::BritterMcQuaidPuff, eqs=DefaultSet) end -function (pf::BritterMcQuaidPuffSolution{F,I})(x,y,z,t) where{F,I} +function (pf::BritterMcQuaidPuffSolution{F,I})(x,y,z,t)::F where{F<:Number,I} R₀ = cbrt(3*pf.V₀/4π) xc = 0.4*pf.u₁₀*t @@ -104,13 +114,13 @@ function (pf::BritterMcQuaidPuffSolution{F,I})(x,y,z,t) where{F,I} #domain check if r² > R² || z < 0 - return 0.0 + return zero(F) end - x′ = x/cbrt(pf.V₀) + x′ = (xc+√(R²))/cbrt(pf.V₀) if x′ ≤ pf.xnf # use near-field correlation - c′ = x′ > 0 ? 3.24 / (3.24 + x′^2) : 1.0 + c′ = x′ > 0 ? 3.24 / (3.24 + x′^2) : one(F) # don't drop below the first correlation concentration c′ = max(c′, maximum(pf.itp.u)) @@ -133,6 +143,6 @@ function (pf::BritterMcQuaidPuffSolution{F,I})(x,y,z,t) where{F,I} if z ≤ H return c else - return 0.0 + return zero(F) end end diff --git a/test/models/britter_mcquaid_puff_tests.jl b/test/models/britter_mcquaid_puff_tests.jl index 074bdfd..f3eb518 100644 --- a/test/models/britter_mcquaid_puff_tests.jl +++ b/test/models/britter_mcquaid_puff_tests.jl @@ -18,11 +18,10 @@ end @testset "Britter-McQuaid puff tests" begin # Britter-McQuaid example, *Workbook on the Dispersion of Dense Gases* # Potchefstroom Accident - s = Substance(name = :ammonia, - molar_weight=0, + s = Substance(name = :pseudo_ammonia_air_mix, + molar_weight=0.017031, vapor_pressure=0, - gas_density = 1.434, # initial density of cloud - # gas_density = 0.88997, # Ammonia, NIST Webbook + gas_density = 1.434, liquid_density = 681.63, # Ammonia, NIST Webbook reference_temp=(273.15-33.316), # boiling point of Ammonia, NIST Webbook reference_pressure=101325, @@ -31,7 +30,7 @@ end latent_heat = 8.17/0.0160425, # Ammonia, NIST Webbook gas_heat_capacity = 1.6151, # Ammonia, NIST Webbook liquid_heat_capacity = 2.0564) # Ammonia, NIST Webbook - r = HorizontalJet( mass_rate = 38e3, # 38 tonnes of Ammonia + r = HorizontalJet( mass_rate = 7.25e4*1.434, duration = 1, diameter = 1, velocity = 7.25e4/(0.25*π), @@ -43,29 +42,35 @@ end scn = Scenario(s,r,a) # known answers # initial concentration - c₀ = 38e3/7.25e4 + c₀ = 1/16.4 # in the near-field - x₁, t₁, c₁ = 160, 200, 0.11109705839813348/1.434 - # example, in the interpolation region - x₂, t₂, c₂ = 355, 200, 0.04368259217437896 - # far field - x₃, t₃, c₃ = 3133, 3000, 0.0004464242323595325 + x₁, t₁, c₁ = 160, 160/(0.4*2), 0.0030830245167018195 + # examples, in the interpolation region + x₂, t₂, c₂ = 333, 124, 0.006097560975609757 + x₃, t₃, c₃ = 333, 1396, 0.00018497619507887978 + # example in the far-field + x₄, t₄, c₄ = 2500, 2500/(0.4*2), 0.0007259996788424369/16.4 # test overall solution - pf = puff(scn, BritterMcQuaidPuff()) + pf = puff(scn, BritterMcQuaidPuff(); temp_correction=false) @test isa(pf, GasDispersion.BritterMcQuaidPuffSolution) @test isa(pf, Puff) - @test pf(x₁,0,0,t₁) ≈ c₁ - @test pf(x₂,0,0,t₂) ≈ c₂ - @test pf(x₃,0,0,t₃) ≈ c₃ + @test c₀*pf(x₁,0,0,t₁) ≈ c₁ + @test c₀*pf(x₂,0,0,t₂) ≈ c₂ + @test c₀*pf(x₃,0,0,t₃) ≈ c₃ + @test c₀*pf(x₄,0,0,t₄) ≈ c₄ # test puff extent R² = cbrt(3*pf.V₀/4π)^2 + 1.2*√(pf.gₒ*pf.V₀)*t₁ - @test pf(x₁, √(R²), 0, t₁) ≈ c₁ + @test c₀*pf(x₁, √(R²), 0, t₁) ≈ c₁ @test pf(x₁, √(R²) + 1, 0, t₁) == 0.0 - H = (pf.c₀*pf.V₀)/(c₁*π*R²) - @test pf(x₁, 0, H, t₁) ≈ c₁ + H = (c₀*pf.V₀)/(c₁*π*R²) + @test c₀*pf(x₁, 0, H, t₁) ≈ c₁ @test pf(x₁, 0, H + 1, t₁) == 0.0 + # test temperature correction + pf = puff(scn, BritterMcQuaidPuff(); temp_correction=true) + @test pf.T′ ≈ (273.15-33.316) / 293.15 + end