From a2bcfba091b77d367bf3488b30f633d70fc00182 Mon Sep 17 00:00:00 2001 From: aefarrell Date: Fri, 3 Feb 2023 15:50:15 -0700 Subject: [PATCH 1/3] Britter-McQuaid puff model Added correlations for the Britter-McQuaid instantaneous ground level release puff model. --- src/models/britter_mcquaid_plume.jl | 13 ++- src/models/britter_mcquaid_puff.jl | 123 ++++++++++++++++++++++++- src/utils/britter_mcquaid_correls.jl | 133 ++++++++++++++++++++++++++- src/utils/utils.jl | 3 + 4 files changed, 257 insertions(+), 15 deletions(-) diff --git a/src/models/britter_mcquaid_plume.jl b/src/models/britter_mcquaid_plume.jl index 1d5621b..4cbaad1 100644 --- a/src/models/britter_mcquaid_plume.jl +++ b/src/models/britter_mcquaid_plume.jl @@ -1,5 +1,3 @@ -include("../utils/britter_mcquaid_correls.jl") - struct BritterMcQuaidPlume <: PlumeModel end struct BritterMcQuaidPlumeSolution <: Plume @@ -18,11 +16,12 @@ end """ plume(::Scenario, BritterMcQuaidPlume) -Returns the solution to a Britter-McQuaid dispersion model for the given -scenario. +Returns the solution to the Britter-McQuaid continuous ground level release +model for the given scenario. # References + Britter, R.E. and J. McQuaid, *Workbook on the Dispersion of Dense Gases* HSE Contract Research Report No. 17/1988, 1988 ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) """ function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume}) @@ -63,7 +62,7 @@ function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume}) # linear interpolation between the near-field correlation # and the main correlation βnf = log10(xnf) - cnf = 306*xnf^-2 / (1+ 306*xnf^-2) + cnf = 306 / (xnf^2 + 306) βs = [ βnf; βs] concs = [ cnf; concs ] else @@ -91,7 +90,7 @@ function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume}) T′, # temperature_correction D, # critical length, m lb, # length parameter, m - itp, # interpolation::Extrapolation + itp, # interpolation xnf, # near-field distance xff, # far-field distance A # far-field constant @@ -126,7 +125,7 @@ function (b::BritterMcQuaidPlumeSolution)(x, y, z) x′ = x/b.D if x′ ≤ b.xnf # use near-field correlation - c′ = x′ > 0 ? 306*x′^-2 / (1+ 306*x′^-2) : 1.0 + c′ = x′ > 0 ? 306 / ( 306 + x′^2 ) : 1.0 elseif b.xnf < x′ < b.xff # use linear interpolation β = log10(x′) diff --git a/src/models/britter_mcquaid_puff.jl b/src/models/britter_mcquaid_puff.jl index 8ddeeaf..4ea7866 100644 --- a/src/models/britter_mcquaid_puff.jl +++ b/src/models/britter_mcquaid_puff.jl @@ -1,14 +1,127 @@ struct BritterMcQuaidPuff <: PuffModel end +struct BritterMcQuaidPuffSolution <: Puff + scenario::Scenario + model::Symbol + c₀::Number # initial concentration, kg/m³ + T′::Number # temperature_correction + V₀::Number # initial volume, m³ + gₒ::Number # reduced gravity, m/s² + u₁₀::Number # referebce windspeed, m/s + itp::Interpolations.GriddedInterpolation + xnf::Number # near-field distance + xff::Number # far-field distance + A::Number # far-field constant +end + """ puff(scenario::Scenario, BritterMcQuaidPuff) -Generates a Britter-McQuaid dispersion model on the given scenario and returns a -callable struct giving the concentration of the form -c(x, y, z, t) +Returns the solution to the Britter-McQuaid instantaneous ground level release +model for the given scenario. + +# References ++ Britter, R.E. and J. McQuaid, *Workbook on the Dispersion of Dense Gases* HSE Contract Research Report No. 17/1988, 1988 ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) """ function puff(scenario::Scenario, ::Type{BritterMcQuaidPuff}) - # TODO - error("Britter-McQuaid puff dispersion model is not currently implemented") + + Q = _release_flowrate(scenario) + ṁ = _mass_rate(scenario) + t = _duration(scenario) + ρⱼ = _release_density(scenario) + Tᵣ = _release_temperature(scenario) + + u₁₀ = _windspeed(scenario, 10.0) + ρₐ = _atmosphere_density(scenario) + Tₐ = _atmosphere_temperature(scenario) + + # initial cloud dimensions + V₀ = Q*t + + # initial concentration + c₀ = ṁ/Q + + # temperature correction + T′ = Tᵣ/Tₐ + + # relative density + g = 9.80616 # gravity, m/s^2 + gₒ = g * ((ρⱼ - ρₐ)/ ρₐ) + + # correlation parameter + α = 0.5*log10( gₒ * cbrt(V₀) / u₁₀^2 ) + + # setting up correlation for constant α + # calculates the points for the linear interpolation + concs, βs = _bm_pf_c(α) + + # near-field location + xnf = 10^(minimum(βs)) + + # linear interpolation + itp = interpolate((βs,), concs, Gridded(Linear())) + + # far field correlation + # starts at last interpolation point and decays like x′^-2 + xff = 10^(maximum(βs)) + A = minimum(concs)*xff^2 + + return BritterMcQuaidPuffSolution( + scenario, #scenario::Scenario + :brittermcquaid, #model::Symbol + c₀, # initial concentration, kg/m³ + T′, # temperature_correction + V₀, # initial volume, m³ + gₒ, # reduced gravity, m/s² + u₁₀, # referebce windspeed, m/s + itp, # interpolation::Extrapolation + xnf, # near-field distance + xff, # far-field distance + A # far-field constant + ) + +end + +function (pf::BritterMcQuaidPuffSolution)(x,y,z,t) + + R₀ = cbrt(3*pf.V₀/4π) + xc = 0.4*pf.u₁₀*t + R² = R₀^2 + 1.2*√(pf.gₒ*pf.V₀)*t + r² = (x-xc)^2 + y^2 + + #domain check + if r² > R² + return 0.0 + end + + x′ = x/cbrt(pf.V₀) + if x′ ≤ pf.xnf + # use near-field correlation + c′ = x′ > 0 ? 3.24 / (3.24 + x′^2) : 1.0 + + # don't drop below the first correlation concentration + c′ = max(c′, maximum(pf.itp.coefs)) + elseif pf.xnf < x′ < pf.xff + # use linear interpolation + β = log10(x′) + c′ = pf.itp(β) + else + # use far-field correlation + # where A is a function of α + c′ = pf.A/(x′^2) + end + + # non-isothermal correction + c′ = c′ / (c′ + (1 - c′)*pf.T′) + c = pf.c₀*c′ + + # vertical extent, from continuity + H = pf.V₀/(c′*π*R²) + if z ≤ H + return c + else + return 0.0 + end end diff --git a/src/utils/britter_mcquaid_correls.jl b/src/utils/britter_mcquaid_correls.jl index 5fd4aa7..c288b76 100644 --- a/src/utils/britter_mcquaid_correls.jl +++ b/src/utils/britter_mcquaid_correls.jl @@ -1,4 +1,4 @@ - +# plume correlations function _bm_pl_c_1(α) # corresponds to (Cm/C0) = 0.10 if α ≤ -0.55 @@ -96,10 +96,10 @@ end """ _bm_pl_c(α) -Britter-McQuaid plume correlations as digtized in TSCREEN +Britter-McQuaid plume correlations # References -+ EPA, *User's Guide to TSCREEN* U.S. Environmental Protection Agency EPA-454/B-94-023 (1994) ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) """ function _bm_pl_c(α) @@ -114,3 +114,130 @@ function _bm_pl_c(α) ] return concs, βs end + + +############################################################################### +# puff correlations + +function _bm_pf_c_1(α) + # corresponds to (Cm/C0) = 0.10 + if α ≤ -0.44 + return 0.70 + elseif -0.44 < α ≤ 0.43 + return 0.26*α + 0.81 + elseif 0.43 < α ≤ 1.0 + return 0.93 + else + @warn "α= $α is out of range" + return 0.93 + end +end + +function _bm_pf_c_05(α) + # corresponds to (Cm/C0) = 0.05 + if α ≤ -0.56 + return 0.85 + elseif -0.56 < α ≤ 0.31 + return 0.26*α + 1.00 + elseif 0.31 < α ≤ 1.0 + return -0.12*α + 1.12 + else + @warn "α= $α is out of range" + α = min(α,10.0) + return -0.12*α + 1.12 + end +end + +function _bm_pf_c_02(α) + # corresponds to (Cm/C0) = 0.02 + if α ≤ -0.66 + return 0.95 + elseif -0.66 < α ≤ 0.32 + return 0.36*α + 1.19 + elseif 0.32 < α ≤ 1.0 + return -0.26*α + 1.38 + else + @warn "α= $α is out of range" + α = min(α,10.0) + return -0.26*α + 1.38 + end +end + +function _bm_pf_c_01(α) + # corresponds to (Cm/C0) = 0.01 + if α ≤ -0.71 + return 1.15 + elseif -0.71 < α ≤ 0.37 + return 0.34*α + 1.39 + elseif 0.37 < α ≤ 1.0 + return -0.38*α + 1.66 + else + @warn "α= $α is out of range" + α = min(α,10.0) + return -0.38*α + 1.66 + end +end + +function _bm_pf_c_005(α) + # corresponds to (Cm/C0) = 0.005 + if α ≤ -0.52 + return 1.48 + elseif -0.52 < α ≤ 0.24 + return 0.26*α + 1.62 + elseif 0.24 < α ≤ 1.0 + return -0.30*α + 1.75 + else + @warn "α= $α is out of range" + α = min(α,10.0) + return -0.30*α + 1.75 + end +end + +function _bm_pf_c_002(α) + # corresponds to (Cm/C0) = 0.002 + if α ≤ 0.27 + return 1.83 + elseif 0.27 < α ≤ 1.0 + return -0.32*α + 1.92 + else + @warn "α= $α is out of range" + α = min(α,10.0) + return -0.32*α + 1.92 + end +end + +function _bm_pf_c_001(α) + # corresponds to (Cm/C0) = 0.001 + if α ≤ -0.10 + return 2.075 + elseif -0.10 < α ≤ 1.0 + return -0.27*α + 2.05 + else + @warn "α= $α is out of range" + α = min(α,10.0) + return -0.27*α + 2.05 + end +end + +""" + _bm_pf_c(α) + +Britter-McQuaid puff correlations + +# References ++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999) + +""" +function _bm_pf_c(α) + concs = [0.10, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001] + βs = [ + _bm_pf_c_1(α), # (C_m/C_0) = 0.1 + _bm_pf_c_05(α), # (C_m/C_0) = 0.05 + _bm_pf_c_02(α), # (C_m/C_0) = 0.02 + _bm_pf_c_01(α), # (C_m/C_0) = 0.01 + _bm_pf_c_005(α), # (C_m/C_0) = 0.005 + _bm_pf_c_002(α), # (C_m/C_0) = 0.002 + _bm_pf_c_001(α) # (C_m/C_0) = 0.001 + ] + return concs, βs +end diff --git a/src/utils/utils.jl b/src/utils/utils.jl index 15147e4..14b9052 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -79,3 +79,6 @@ include("pasquill_gifford.jl") # Briggs' model for plume rise include("plume_rise.jl") + +# model correlations +include("britter_mcquaid_correls.jl") From a2e3c71ad63fa3eaaf24a6f577d807b342131774 Mon Sep 17 00:00:00 2001 From: aefarrell Date: Wed, 22 Feb 2023 15:24:07 -0700 Subject: [PATCH 2/3] Adding unit tests to Britter-McQuaid Puff model --- src/models/britter_mcquaid_puff.jl | 2 +- test/models/britter_mcquaid_puff_tests.jl | 89 +++++++++++++++++------ 2 files changed, 67 insertions(+), 24 deletions(-) diff --git a/src/models/britter_mcquaid_puff.jl b/src/models/britter_mcquaid_puff.jl index 4ea7866..6c6136a 100644 --- a/src/models/britter_mcquaid_puff.jl +++ b/src/models/britter_mcquaid_puff.jl @@ -92,7 +92,7 @@ function (pf::BritterMcQuaidPuffSolution)(x,y,z,t) r² = (x-xc)^2 + y^2 #domain check - if r² > R² + if r² > R² || z < 0 return 0.0 end diff --git a/test/models/britter_mcquaid_puff_tests.jl b/test/models/britter_mcquaid_puff_tests.jl index 05d89c9..06ccd85 100644 --- a/test/models/britter_mcquaid_puff_tests.jl +++ b/test/models/britter_mcquaid_puff_tests.jl @@ -1,26 +1,69 @@ +@testset "Britter-McQuaid puff correlations" begin + include("../../src/utils/britter_mcquaid_correls.jl") + + @test _bm_pf_c(-0.75)[2] ≈ [0.7, 0.85, 0.95, 1.15, 1.48, 1.83, 2.075] + @test _bm_pf_c(0.0)[2] ≈ [0.81, 1.0, 1.19, 1.39, 1.62, 1.83, 2.05] + @test _bm_pf_c(0.75)[2] ≈ [0.93, 1.03, 1.1849999999999998, 1.375, 1.525, 1.68, 1.8474999999999997] + + # test domain warning + @test_logs (:warn, "α= 1.1 is out of range") _bm_pf_c_1(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pf_c_05(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pf_c_02(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pf_c_01(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pf_c_005(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pf_c_002(1.1) + @test_logs (:warn, "α= 1.1 is out of range") _bm_pf_c_001(1.1) + +end + @testset "Britter-McQuaid puff tests" begin - sub = Substance(name="test gas", - gas_density=1.2, - liquid_density=1000., - reference_temp=298., - reference_pressure=101325.0, - boiling_temp=100., - latent_heat=1., - gas_heat_capacity=1., - liquid_heat_capacity=1.) - rel = Release(mass_rate=1.0, - duration=Inf, - diameter=1.0, - velocity=1.0, - height=1.0, - temperature=298., - pressure=101325., - fraction_liquid=0.0) - atm = DryAir(temperature=298.0, - pressure=101325.0, - windspeed=2.0, - stability=ClassF) - scn = Scenario(sub,rel,atm) - @test_throws ErrorException puff(scn, BritterMcQuaidPuff) + # Britter-McQuaid example, *Workbook on the Dispersion of Dense Gases* + # Potchefstroom Accident + s = Substance(name = :ammonia, + gas_density = 1.434, # initial density of cloud + # gas_density = 0.88997, # Ammonia, NIST Webbook + liquid_density = 681.63, # Ammonia, NIST Webbook + reference_temp=(273.15-33.316), # boiling point of Ammonia, NIST Webbook + reference_pressure=101325.0, + boiling_temp = (273.15-33.316), # Ammonia, NIST Webbook + 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 = Release( mass_rate = 38e3, # 38 tonnes of Ammonia + duration = 1.0, + diameter = 1.0, + velocity = 7.25e4/(0.25*π), + height = 0.0, + pressure = 101325.0, + temperature = (273.15-33.316), # Ammonia, NIST Webbook + fraction_liquid = 0.0) + a = DryAir(windspeed=2.0, temperature=293.15, stability=ClassF) + scn = Scenario(s,r,a) + # known answers + # initial concentration + c₀ = 38e3/7.25e4 + # in the near-field + x₁, t₁, c₁ = 160, 200.0, 0.11109705839813348 + # example, in the interpolation region + x₂, t₂, c₂ = 355, 200.0, 0.06264084534270946 + # far field + x₃, t₃, c₃ = 3133, 3000.0, 0.0006403658290584752 + + # test overall solution + pf = puff(scn, BritterMcQuaidPuff) + @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 puff extent + R² = cbrt(3*pf.V₀/4π)^2 + 1.2*√(pf.gₒ*pf.V₀)*t₁ + @test 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₁ + @test pf(x₁, 0, H + eps(Float64), t₁) == 0.0 end From e3c7e44d3834e27faed1a2ddf9b4a21c15f721ff Mon Sep 17 00:00:00 2001 From: aefarrell Date: Thu, 23 Feb 2023 06:31:08 -0700 Subject: [PATCH 3/3] Added Britter McQuaid puff model to documentation --- docs/src/index.md | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/src/index.md b/docs/src/index.md index 049ad58..3087a18 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -167,3 +167,14 @@ default behaviour is to take the limit $n \to \infty$. ```@docs puff(::Scenario, ::Type{IntPuff}) ``` + +### Britter-McQuaid Model + +The Britter-McQuaid model is an empirical correlation for dense cloud +dispersion. The model generates an interpolation function for the average cloud +concentration and the cloud is rendered as a cylinder. + + +```@docs +puff(::Scenario, ::Type{BritterMcQuaidPuff}) +```