From 14d0acc528e59bc4ecfecbeca4ae2935038d6800 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Thu, 26 Jun 2025 18:37:30 -0600 Subject: [PATCH 1/6] Lazy first pass at Palazzi Just moving the approximation used in IntPuff to its own file and calling it "Palazzi". This is step 1. Step 2 is adding all the other features Palazzi may require and changing the default behaviour. --- src/GasDispersion.jl | 3 +- src/models/intpuff.jl | 49 ++++++-------------- src/models/palazzi_puff.jl | 86 ++++++++++++++++++++++++++++++++++++ test/models/intpuff_tests.jl | 2 +- 4 files changed, 102 insertions(+), 38 deletions(-) create mode 100644 src/models/palazzi_puff.jl diff --git a/src/GasDispersion.jl b/src/GasDispersion.jl index fc4bf9d..1486b49 100644 --- a/src/GasDispersion.jl +++ b/src/GasDispersion.jl @@ -24,7 +24,7 @@ export GaussianPlume, SimpleJet, BritterMcQuaidPlume # puff models export PuffModel, Puff, puff -export GaussianPuff, IntPuff, BritterMcQuaidPuff, SLAB +export GaussianPuff, Palazzi, IntPuff, BritterMcQuaidPuff, SLAB # equation sets export EquationSet, BasicEquationSet, DefaultSet, DefaultPuffSet @@ -117,6 +117,7 @@ puff(s; kwargs...) = puff(s, GaussianPuff; kwargs...) # puff models include("models/gaussian_puff.jl") +include("models/palazzi_puff.jl") include("models/intpuff.jl") include("models/britter_mcquaid_puff.jl") include("models/slab_puff.jl") diff --git a/src/models/intpuff.jl b/src/models/intpuff.jl index 83f1e2c..a1279de 100644 --- a/src/models/intpuff.jl +++ b/src/models/intpuff.jl @@ -51,7 +51,19 @@ function puff(scenario::Scenario, ::Type{IntPuff}, eqs=DefaultPuffSet(); n::Numb Pₐ = _atmosphere_pressure(scenario) ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ) - if n > 1 + if n == Inf + return PalazziSolution( + scenario, #scenario::Scenario + :intpuff, #model::Symbol + ṁ, # massrate + ρₐ, # mass-to-vol + Δt, # duration + h, # release height + u, # windspeed + stab, # stability class + eqs # equation set + ) + elseif n > 1 return IntPuffSolution( scenario, #scenario::Scenario :intpuff, #model::Symbol @@ -124,38 +136,3 @@ function (ip::IntPuffSolution{F,<:Integer,S,E})(x,y,z,t) where {F<:Number,S<:Sta return min(c,one(F)) end - -function (ip::IntPuffSolution{F,<:AbstractFloat,S,E})(x,y,z,t) where {F<:Number,S<:StabilityClass,E<:EquationSet} - # domain check - if (x<0)||(z<0)||(t<0) - return zero(F) - end - - Qi = ip.rate - Δt = ip.duration - h = ip.height - u = ip.windspeed - - # Only account for puffs that have already been emitted - Δt = min(t,Δt) - - # Gaussian dispersion in the x direction - σx_a = downwind_dispersion(u*(t-Δt), S, ip.equationset) - σx_b = downwind_dispersion(u*t, S, ip.equationset) - a = (x-u*(t-Δt))/(√2*σx_a) - b = (x-u*t)/(√2*σx_b) - ∫gx = erf(b,a)/(2u) - - # Gaussian dispersion in the y direction - σy = crosswind_dispersion(x, S, ip.equationset) - gy = exp((-1/2)*(y/σy)^2)/(√(2π)*σy) - - # Gaussian dispersion in the z direction - σz = vertical_dispersion(x, S, ip.equationset) - gz = ( exp((-1/2)*((z-h)/σz)^2) + exp((-1/2)*((z+h)/σz)^2) )/(√(2π)*σz) - - # concentration - c = Qi*∫gx*gy*gz/ip.mass_to_vol - - return min(c,one(F)) -end diff --git a/src/models/palazzi_puff.jl b/src/models/palazzi_puff.jl new file mode 100644 index 0000000..e8763c2 --- /dev/null +++ b/src/models/palazzi_puff.jl @@ -0,0 +1,86 @@ +# defining type for dispatch +struct Palazzi <: PuffModel end + +struct PalazziSolution{F<:Number,S<:StabilityClass,E<:EquationSet} <: Puff + scenario::Scenario + model::Symbol + rate::F + mass_to_vol::F + duration::F + height::F + windspeed::F + stability::Type{S} + equationset::E +end +PalazziSolution(s,m,r,ρ,d,h,u,stab,es) = PalazziSolution(s,m,promote(r,ρ,d,h,u,)...,stab,es) + +@doc doc""" + puff(::Scenario, Palazzi[, ::EquationSet]; kwargs...) + +... + +# References ++ Palazzi, E, M De Faveri, Giuseppe Fumarola, and G Ferraiolo. “Diffusion from a Steady Source of Short Duration.” *Atmospheric Environment*. 16, no. 12 (1982): 2785–90. + +""" +function puff(scenario::Scenario, ::Type{Palazzi}, eqs=DefaultPuffSet()) + + stab = _stability(scenario) + ṁ = _mass_rate(scenario) + Δt = _duration(scenario) + h = _release_height(scenario) + u = windspeed(scenario,h,eqs) + + # jet at ambient conditions + Tₐ = _atmosphere_temperature(scenario) + Pₐ = _atmosphere_pressure(scenario) + ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ) + + return PalazziSolution( + scenario, #scenario::Scenario + :intpuff, #model::Symbol + ṁ, # massrate + ρₐ, # mass-to-vol + Δt, # duration + h, # release height + u, # windspeed + stab, # stability class + eqs # equation set + ) +end + + +function (ip::PalazziSolution{F,S,E})(x,y,z,t) where {F<:Number,S<:StabilityClass,E<:EquationSet} + # domain check + if (x<0)||(z<0)||(t<0) + return zero(F) + end + + Qi = ip.rate + Δt = ip.duration + h = ip.height + u = ip.windspeed + + # Only account for puffs that have already been emitted + Δt = min(t,Δt) + + # Gaussian dispersion in the x direction + σx_a = downwind_dispersion(u*(t-Δt), S, ip.equationset) + σx_b = downwind_dispersion(u*t, S, ip.equationset) + a = (x-u*(t-Δt))/(√2*σx_a) + b = (x-u*t)/(√2*σx_b) + ∫gx = erf(b,a)/(2u) + + # Gaussian dispersion in the y direction + σy = crosswind_dispersion(x, S, ip.equationset) + gy = exp((-1/2)*(y/σy)^2)/(√(2π)*σy) + + # Gaussian dispersion in the z direction + σz = vertical_dispersion(x, S, ip.equationset) + gz = ( exp((-1/2)*((z-h)/σz)^2) + exp((-1/2)*((z+h)/σz)^2) )/(√(2π)*σz) + + # concentration + c = Qi*∫gx*gy*gz/ip.mass_to_vol + + return min(c,one(F)) +end \ No newline at end of file diff --git a/test/models/intpuff_tests.jl b/test/models/intpuff_tests.jl index 8f3a8ca..427f316 100644 --- a/test/models/intpuff_tests.jl +++ b/test/models/intpuff_tests.jl @@ -32,7 +32,7 @@ @test GasDispersion.IntPuffSolution(scn,:test_promotion,1.0,2,3,4,5,6,ClassA,DefaultPuffSet()) isa GasDispersion.IntPuffSolution{Float64, Int64, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, GasDispersion.CCPSPuffσx, GasDispersion.CCPSPuffσy, GasDispersion.CCPSPuffσz}} @test puff(scn, IntPuff;n=1) isa GasDispersion.GaussianPuffSolution @test puff(scn, IntPuff;n=3) isa GasDispersion.IntPuffSolution{<:Number,<:Integer,<:StabilityClass,<:EquationSet} - @test puff(scn, IntPuff) isa GasDispersion.IntPuffSolution{<:Number,<:Float64,<:StabilityClass,<:EquationSet} + @test puff(scn, IntPuff) isa GasDispersion.PalazziSolution{<:Number,<:StabilityClass,<:EquationSet} @test_throws ErrorException puff(scn, IntPuff; n=0) # testing 3 puffs From 325220576503557f7bed6821b48980a2e4d80570 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Fri, 27 Jun 2025 14:29:25 -0600 Subject: [PATCH 2/6] Generalizing Palazzi is generic in regards to the plume model, it is simply $\Chi$ in the equation --- src/models/intpuff.jl | 4 +--- src/models/palazzi_puff.jl | 40 ++++++++++---------------------------- 2 files changed, 11 insertions(+), 33 deletions(-) diff --git a/src/models/intpuff.jl b/src/models/intpuff.jl index a1279de..9e0400b 100644 --- a/src/models/intpuff.jl +++ b/src/models/intpuff.jl @@ -55,10 +55,8 @@ function puff(scenario::Scenario, ::Type{IntPuff}, eqs=DefaultPuffSet(); n::Numb return PalazziSolution( scenario, #scenario::Scenario :intpuff, #model::Symbol - ṁ, # massrate - ρₐ, # mass-to-vol + plume(scenario, GaussianPlume, eqs), Δt, # duration - h, # release height u, # windspeed stab, # stability class eqs # equation set diff --git a/src/models/palazzi_puff.jl b/src/models/palazzi_puff.jl index e8763c2..1534190 100644 --- a/src/models/palazzi_puff.jl +++ b/src/models/palazzi_puff.jl @@ -4,15 +4,13 @@ struct Palazzi <: PuffModel end struct PalazziSolution{F<:Number,S<:StabilityClass,E<:EquationSet} <: Puff scenario::Scenario model::Symbol - rate::F - mass_to_vol::F + plume::Plume duration::F - height::F windspeed::F stability::Type{S} equationset::E end -PalazziSolution(s,m,r,ρ,d,h,u,stab,es) = PalazziSolution(s,m,promote(r,ρ,d,h,u,)...,stab,es) +PalazziSolution(s,m,p,d,u,stab,es) = PalazziSolution(s,m,p,promote(d,u)...,stab,es) @doc doc""" puff(::Scenario, Palazzi[, ::EquationSet]; kwargs...) @@ -23,26 +21,18 @@ PalazziSolution(s,m,r,ρ,d,h,u,stab,es) = PalazziSolution(s,m,promote(r,ρ,d,h,u + Palazzi, E, M De Faveri, Giuseppe Fumarola, and G Ferraiolo. “Diffusion from a Steady Source of Short Duration.” *Atmospheric Environment*. 16, no. 12 (1982): 2785–90. """ -function puff(scenario::Scenario, ::Type{Palazzi}, eqs=DefaultPuffSet()) +function puff(scenario::Scenario, ::Type{Palazzi}, eqs=DefaultPuffSet(); plume_model=GaussianPlume) stab = _stability(scenario) - ṁ = _mass_rate(scenario) Δt = _duration(scenario) h = _release_height(scenario) u = windspeed(scenario,h,eqs) - # jet at ambient conditions - Tₐ = _atmosphere_temperature(scenario) - Pₐ = _atmosphere_pressure(scenario) - ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ) - return PalazziSolution( scenario, #scenario::Scenario - :intpuff, #model::Symbol - ṁ, # massrate - ρₐ, # mass-to-vol + :palazzi, #model::Symbol + plume(scenario, plume_model, eqs), # plume model Δt, # duration - h, # release height u, # windspeed stab, # stability class eqs # equation set @@ -56,31 +46,21 @@ function (ip::PalazziSolution{F,S,E})(x,y,z,t) where {F<:Number,S<:StabilityClas return zero(F) end - Qi = ip.rate - Δt = ip.duration - h = ip.height - u = ip.windspeed + Χ = ip.plume(x,y,z,t) - # Only account for puffs that have already been emitted + Δt = ip.duration Δt = min(t,Δt) + u = ip.windspeed # Gaussian dispersion in the x direction σx_a = downwind_dispersion(u*(t-Δt), S, ip.equationset) σx_b = downwind_dispersion(u*t, S, ip.equationset) a = (x-u*(t-Δt))/(√2*σx_a) b = (x-u*t)/(√2*σx_b) - ∫gx = erf(b,a)/(2u) - - # Gaussian dispersion in the y direction - σy = crosswind_dispersion(x, S, ip.equationset) - gy = exp((-1/2)*(y/σy)^2)/(√(2π)*σy) - - # Gaussian dispersion in the z direction - σz = vertical_dispersion(x, S, ip.equationset) - gz = ( exp((-1/2)*((z-h)/σz)^2) + exp((-1/2)*((z+h)/σz)^2) )/(√(2π)*σz) + ∫gx = erf(b,a)/2 # concentration - c = Qi*∫gx*gy*gz/ip.mass_to_vol + c = Χ*∫gx return min(c,one(F)) end \ No newline at end of file From 2d6fa9c599c7f219339466c2317f9ec349761776 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Fri, 27 Jun 2025 18:10:40 -0600 Subject: [PATCH 3/6] Palazzi unit tests Quick unit tests to ensure coverage. --- test/models/palazzi_tests.jl | 37 ++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 38 insertions(+) create mode 100644 test/models/palazzi_tests.jl diff --git a/test/models/palazzi_tests.jl b/test/models/palazzi_tests.jl new file mode 100644 index 0000000..9cdf1c2 --- /dev/null +++ b/test/models/palazzi_tests.jl @@ -0,0 +1,37 @@ +@testset "Palazzi short duration puff tests" begin + # Gaussian plume example, *Guidelines for Consequence Analysis of Chemical + # Releases* CCPS, 1999, pg 101 + sub = Substance(name="test gas", + molar_weight=1, + vapor_pressure=nothing, + gas_density=1.2268, + liquid_density=1000, + reference_temp=298, + reference_pressure=101325, + k=0, + boiling_temp=100, + latent_heat=1, + gas_heat_capacity=1, + liquid_heat_capacity=1) + rel = HorizontalJet(mass_rate=0.1, + duration=10, + diameter=1, + velocity=1, + height=10, + temperature=298, + pressure=101325, + fraction_liquid=0) + atm = SimpleAtmosphere(temperature=298, + pressure=101325, + windspeed=2, + stability=ClassF) + scn = Scenario(sub,rel,atm) + x₁, t₁, Δt, h = 500, 250, 10, 10 + + # testing default behaviour + @test puff(scn, Palazzi) isa GasDispersion.PalazziSolution{<:Number,<:StabilityClass,<:EquationSet} + pf = puff(scn, Palazzi) + @test pf(x₁,0,h,t₁) ≈ 0.00035586710616021256/1.2268 + @test pf(x₁,0,h,-t₁) == 0.0 + +end diff --git a/test/runtests.jl b/test/runtests.jl index af95660..223e87c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,6 +32,7 @@ if GROUP == "All" || GROUP == "Model" include("models/gaussian_plume_tests.jl") include("models/gaussian_puff_tests.jl") include("models/intpuff_tests.jl") + include("models/palazzi_tests.jl") include("models/simple_jet_tests.jl") include("models/britter_mcquaid_plume_tests.jl") include("models/britter_mcquaid_puff_tests.jl") From d39267b6cea0ae4cf2b6d33f37e1f812ba758157 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Fri, 27 Jun 2025 21:19:16 -0600 Subject: [PATCH 4/6] Differing behaviour and defaults Palazzi can be calculated differently depending on how $\sigma_x$ is calculated, and references disagree on how this should be done. Added a keyword argument to select the different options. Palazzi uses plume dispersion parameters by default, using $\sigma_x = \sigma_y$ --- src/models/intpuff.jl | 1 + src/models/palazzi_puff.jl | 60 ++++++++++++++++++++++++++++-------- test/models/palazzi_tests.jl | 11 +++++-- 3 files changed, 57 insertions(+), 15 deletions(-) diff --git a/src/models/intpuff.jl b/src/models/intpuff.jl index 9e0400b..d3f7cf5 100644 --- a/src/models/intpuff.jl +++ b/src/models/intpuff.jl @@ -55,6 +55,7 @@ function puff(scenario::Scenario, ::Type{IntPuff}, eqs=DefaultPuffSet(); n::Numb return PalazziSolution( scenario, #scenario::Scenario :intpuff, #model::Symbol + :intpuff, #disp::Symbol plume(scenario, GaussianPlume, eqs), Δt, # duration u, # windspeed diff --git a/src/models/palazzi_puff.jl b/src/models/palazzi_puff.jl index 1534190..945f0d4 100644 --- a/src/models/palazzi_puff.jl +++ b/src/models/palazzi_puff.jl @@ -4,24 +4,47 @@ struct Palazzi <: PuffModel end struct PalazziSolution{F<:Number,S<:StabilityClass,E<:EquationSet} <: Puff scenario::Scenario model::Symbol + disp::Symbol plume::Plume duration::F windspeed::F stability::Type{S} equationset::E end -PalazziSolution(s,m,p,d,u,stab,es) = PalazziSolution(s,m,p,promote(d,u)...,stab,es) + +PalazziDefaultSet = BasicEquationSet{DefaultWind,Defaultσy,Defaultσy,Defaultσz} + +downwind_dispersion(x::Number, stab::Any, ::Type{Defaultσy}) = crosswind_dispersion(x, stab, Defaultσy) @doc doc""" puff(::Scenario, Palazzi[, ::EquationSet]; kwargs...) -... +Returns the solution to a short duration Gaussian puff dispersion model for the given scenario, based on the work of Palazzi *et al*. + +```math +c\left(x,y,z,t\right) = \Chi\left(x,y,z\right) \frac{1}{2} \left( \mathrm{erf} \left( { {x - u (t-\Delta t)} \over \sqrt{2} \sigma_x } \right) - \mathrm{erf} \left( { {x - u t} \over \sqrt{2} \sigma_x } \right) \right) +```` + +where Χ is a Gaussian plume model and the σs are dispersion parameters. The `EquationSet` +defines the set of correlations used to calculate the dispersion parameters and windspeed. +The concentration returned is in volume fraction, assuming the puff is a gas at ambient +conditions. + +There are multiple variations of the Palazzi short duration model, changing how the +downwind dispersion, σx, is calculated: +- `:default` follows Palazzi and calculates σx at the downwind distance *x* +- `:intpuff` calculates σx at the downwind distance to the cloud center *u t* +- `:tno` follows the TNO Yellow Book eqn 4.60b, using the distance *x* while the plume is still attached and the distance to the cloud center, *ut*, afterwards + +# Arguments +- `plume_model::Type{Plume} = GaussianPlume` : the plume model $\Chi$ +- `disp_method = :default` : the method for calculating the downwind dispersion # References + Palazzi, E, M De Faveri, Giuseppe Fumarola, and G Ferraiolo. “Diffusion from a Steady Source of Short Duration.” *Atmospheric Environment*. 16, no. 12 (1982): 2785–90. """ -function puff(scenario::Scenario, ::Type{Palazzi}, eqs=DefaultPuffSet(); plume_model=GaussianPlume) +function puff(scenario::Scenario, ::Type{Palazzi}, eqs=PalazziDefaultSet(); plume_model=GaussianPlume, disp_method=:default) stab = _stability(scenario) Δt = _duration(scenario) @@ -31,6 +54,7 @@ function puff(scenario::Scenario, ::Type{Palazzi}, eqs=DefaultPuffSet(); plume_m return PalazziSolution( scenario, #scenario::Scenario :palazzi, #model::Symbol + disp_method, #dispersion model::Symbol plume(scenario, plume_model, eqs), # plume model Δt, # duration u, # windspeed @@ -40,27 +64,37 @@ function puff(scenario::Scenario, ::Type{Palazzi}, eqs=DefaultPuffSet(); plume_m end -function (ip::PalazziSolution{F,S,E})(x,y,z,t) where {F<:Number,S<:StabilityClass,E<:EquationSet} +function (ps::PalazziSolution{F,S,E})(x,y,z,t) where {F<:Number,S<:StabilityClass,E<:EquationSet} # domain check if (x<0)||(z<0)||(t<0) return zero(F) end - Χ = ip.plume(x,y,z,t) + Χ = ps.plume(x,y,z,t) + Δt = min(t,ps.duration) + u = ps.windspeed - Δt = ip.duration - Δt = min(t,Δt) - u = ip.windspeed + if ps.disp == :default + xa = xb = x + elseif ps.disp == :intpuff + xa = u*(t-Δt) + xb = u*t + elseif ps.disp == :tno + if t < Δt + xa = xb = x + else + xa = xb = u*t + end + end # Gaussian dispersion in the x direction - σx_a = downwind_dispersion(u*(t-Δt), S, ip.equationset) - σx_b = downwind_dispersion(u*t, S, ip.equationset) + σx_a = downwind_dispersion(xa, S, ps.equationset) + σx_b = downwind_dispersion(xb, S, ps.equationset) a = (x-u*(t-Δt))/(√2*σx_a) - b = (x-u*t)/(√2*σx_b) - ∫gx = erf(b,a)/2 + b = (x-u*t)/(√2*σx_b) # concentration - c = Χ*∫gx + c = Χ*erf(b,a)/2 return min(c,one(F)) end \ No newline at end of file diff --git a/test/models/palazzi_tests.jl b/test/models/palazzi_tests.jl index 9cdf1c2..f0686e8 100644 --- a/test/models/palazzi_tests.jl +++ b/test/models/palazzi_tests.jl @@ -30,8 +30,15 @@ # testing default behaviour @test puff(scn, Palazzi) isa GasDispersion.PalazziSolution{<:Number,<:StabilityClass,<:EquationSet} - pf = puff(scn, Palazzi) - @test pf(x₁,0,h,t₁) ≈ 0.00035586710616021256/1.2268 + pf = puff(scn, Palazzi; disp_method=:default) + @test pf(x₁,0,h,t₁) ≈ 1.6373715867403676e-5 @test pf(x₁,0,h,-t₁) == 0.0 + + pf1 = puff(scn, Palazzi, DefaultPuffSet(); disp_method=:intpuff) + @test pf1(x₁,0,h,t₁) ≈ 0.00035586710616021256/1.2268 + + pf2 = puff(scn, Palazzi; disp_method=:tno) + @test pf2(10,0,10,5) ≈ 0.02846251120337852 + @test pf2(x₁,0,h,t₁) ≈ 1.6373715867403676e-5 end From 21477f15b36c00eb6ded24803b0da33131062821 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Fri, 27 Jun 2025 21:33:42 -0600 Subject: [PATCH 5/6] Fixed error in logic --- src/models/palazzi_puff.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/models/palazzi_puff.jl b/src/models/palazzi_puff.jl index 945f0d4..ab45c85 100644 --- a/src/models/palazzi_puff.jl +++ b/src/models/palazzi_puff.jl @@ -80,7 +80,7 @@ function (ps::PalazziSolution{F,S,E})(x,y,z,t) where {F<:Number,S<:StabilityClas xa = u*(t-Δt) xb = u*t elseif ps.disp == :tno - if t < Δt + if t < ps.duration xa = xb = x else xa = xb = u*t From 5551cd2984134e97dcd1c08a976511770219bd25 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Sat, 28 Jun 2025 11:53:34 -0600 Subject: [PATCH 6/6] Palazzi Documentation Documentation of the Palazzi short duration model --- Project.toml | 2 +- docs/src/puff.md | 51 +++++++++++++++++++++++++------------- docs/src/references.md | 2 ++ src/models/palazzi_puff.jl | 12 ++++----- 4 files changed, 43 insertions(+), 24 deletions(-) diff --git a/Project.toml b/Project.toml index f374d7d..ab8d8d2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GasDispersion" uuid = "ac6af613-5479-4178-853e-9d66de417193" authors = ["Allan Farrell and contributors"] -version = "0.3.0" +version = "0.4.0" [deps] DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" diff --git a/docs/src/puff.md b/docs/src/puff.md index 37731e9..aa20ee1 100644 --- a/docs/src/puff.md +++ b/docs/src/puff.md @@ -6,13 +6,15 @@ Puff models are for "instantaneous" releases or other time-dependent releases. puff ``` -## Gaussian Puffs +## Gaussian Puff Models + +### Simple Gaussial Puffs ```@docs puff(::Scenario, ::Type{GaussianPuff}) ``` -A gaussian puff model assumes the release is instantaneous, and all mass is concentrated in a single point. The cloud then disperses as it moves downwind with the concentration profile is given by a series of gaussians with dispersions $\sigma_x$, $\sigma_y$, and $\sigma_z$, which are found from correlations tabulated per stability class. Similarly to the plume model, a ground reflection term is included to correct for the fact that material cannot pass through the ground. +A simple gaussian puff model assumes the release is instantaneous, and all mass is concentrated in a single point. The cloud then disperses as it moves downwind with the concentration profile is given by a series of gaussians with dispersions $\sigma_x$, $\sigma_y$, and $\sigma_z$, which are found from correlations tabulated per stability class. Similarly to the plume model, a ground reflection term is included to correct for the fact that material cannot pass through the ground. ```math c_{puff} = { m_i \over { (2 \pi)^{3/2} \sigma_x \sigma_y \sigma_z } } @@ -32,7 +34,7 @@ with The model assumes the initial release is a single point, with no dimensions. Unlike the plume model, this concentration is a function of time. The model converts the final concentration to volume fraction, assuming the puff is a gas at ambient conditions. -### Downwind dispersion correlations +#### Downwind dispersion correlations The downwind dispersion, $\sigma_{x}$ is a function of downwind distance of the cloud center, $x_c$, as well as stability class @@ -42,7 +44,7 @@ The downwind dispersion, $\sigma_{x}$ is a function of downwind distance of the Where $\delta$ and $\beta$ are identical to those tabulated for the crosswind dispersion. -### Crosswind dispersion correlations +#### Crosswind dispersion correlations The crosswind dispersion, $\sigma_{y}$ is a function of downwind distance of the cloud center, $x_c$, as well as stability class @@ -62,7 +64,7 @@ Where $\delta$, $\beta$, and $\gamma$ are tabulated based on stability class([AI | F | 0.02 | 0.89 | -### Vertical dispersion correlations +#### Vertical dispersion correlations The vertical dispersion, $\sigma_{z}$ is a function of downwind distance of the cloud center, $x_c$, as well as stability class @@ -81,7 +83,7 @@ Where $\delta$ and $\beta$ are tabulated based on stability class([AIChE/CCPS 19 | E | 0.10 | 0.65 | | F | 0.05 | 0.61 | -### Example +#### Example Suppose we wish to model the dispersion of gaseous propane from a leak from a storage tank, where the leak is from a 10mm hole that is 3.5m above the ground and the propane is at 25°C and 4barg. Assume the discharge coefficient $c_{D} = 0.85$. Assume the leak occurs for 10s. This scenario is adapted from CCPS *Guidelines for Consequence Analysis of Chemical Releases*([AIChE/CCPS 1999](references.md), 47) @@ -216,8 +218,27 @@ plot(g, t′; xlims=(85,115), ylims=(-10,10), clims=(0,5e-3), aspect_ratio=:equa end ``` +### Palazzi Short Duration Puff Model +```@docs + puff(::Scenario, ::Type{Palazzi}) +``` +The `Palazzi` model integrates over the Gaussian puff model for a short duration release ([Palazzi 1982](references.md)), where the dispersion parameters $\sigma$ are assumed to be independent time. + +```math +c\left(x,y,z,t\right) = \chi\left(x,y,z\right) \frac{1}{2} \left( \mathrm{erf} \left( { {x - u (t-\Delta t)} \over \sqrt{2} \sigma_x } \right) - \mathrm{erf} \left( { {x - u t} \over \sqrt{2} \sigma_x } \right) \right) +``` + +where $\chi$ is a Gaussian plume model and $\sigma_x$ is the downwind dispersion parameter. The model assumes the initial release is a single point, with no dimensions. Additionally, the model converts the final concentration to volume fraction, assuming the puff is a gas at ambient conditions. -## Integrated Gaussian Puffs +By default the Palazzi model assumes a simple Gaussian plume model, $\chi$, for which this is a "correction", and uses plume dispersion parameters with $\sigma_x \left( x \right) = \sigma_y \left( x \right)$. However the user is free to use *any* plume model which is a subtype of `Plume`, and any equation set that implements downwind dispersion. + +There are multiple variations of the Palazzi short duration model, depending on how the downwind dispersion, $\sigma_x$ , is calculated: +- `:default` follows Palazzi and calculates $\sigma_x$ at the downwind distance *x* +- `:intpuff` calculates $\sigma_x$ at the downwind distance to the cloud center at the start and end of the cloud, $ut$ and $u \left(t-\Delta t\right)$ +- `:tno` follows the TNO Yellow Book eqn 4.60b, using the distance *x* while the plume is still attached to the release point, and the distance to the cloud center, *ut*, afterwards + + +### Integrated Gaussian Puff Model ```@docs puff(::Scenario, ::Type{IntPuff}) @@ -240,14 +261,9 @@ with - $\sigma_y$ - crosswind dispersion, m - $\sigma_z$ - downwind dispersion, m -The model assumes the initial release is a single point, with no dimensions. Additionally, model converts the final concentration to volume fraction, assuming the puff is a gas at ambient conditions. - -### Dispersion Parameters +The model assumes the initial release is a single point, with no dimensions. Additionally, model converts the final concentration to volume fraction, assuming the puff is a gas at ambient conditions. If no number, *n*, is specified the model defaults to an integral approximation similar to the [Palazzi Short Duration Puff Model](@ref). -The dispersion parameters are the same as used for the `GaussianPuff` model. - - -### Example +#### Example Continuing with the propane leak example from above, we now model the release as a sequence of 100 gaussian puffs. Essentially chopping the 10s over which the release happens into 0.1s intervals and releasing one puff per interval at a time for 10s. @@ -287,8 +303,8 @@ ig_inf = puff(scn, IntPuff) plot(ig_inf, 86; xlims=(90,110), ylims=(-10,10), aspect_ratio=:equal) ``` - -## Britter-McQuaid Model +## Box Models +### Britter-McQuaid Model ```@docs puff(::Scenario, ::Type{BritterMcQuaidPuff}) @@ -300,7 +316,8 @@ concentration and the cloud is rendered as a cylinder. The only correlations use in the provided equationset are for windspeed. -## SLAB Jet Model +## Integral Models +### SLAB Jet Model ```@docs puff(::Scenario, ::Type{SLAB}) diff --git a/docs/src/references.md b/docs/src/references.md index 520d70a..7ed086d 100644 --- a/docs/src/references.md +++ b/docs/src/references.md @@ -22,6 +22,8 @@ Lees, Frank P. 1996. *Loss Prevention in the Process Industries, 2nd ed*. Oxford Long, V.D. 1963. "Estimation of the Extent of Hazard Areas Round a Vent." *Chem. Process Hazard*. II:6 +Palazzi, E, M De Faveri, Giuseppe Fumarola, and G Ferraiolo. “Diffusion from a Steady Source of Short Duration.” *Atmospheric Environment*. 16, no. 12 (1982): 2785–90. + Pasquill, Frank. 1974. *Atmospheric Diffusion, 2nd Ed.* New York: Halstead Press, New York Paulson, C. A. 1970. "The Mathematical Representation of Wind Speed and Temperature Profiles in the Unstable Atmospheric Surface Layer." *Journal of Applied Meteorology and Climatology*. 9(6): 857-861. [doi:10.1175/1520-0450(1970)009<0857:TMROWS>2.0.CO;2](https://doi.org/10.1175/1520-0450(1970)009<0857:TMROWS>2.0.CO;2) diff --git a/src/models/palazzi_puff.jl b/src/models/palazzi_puff.jl index ab45c85..73a17c4 100644 --- a/src/models/palazzi_puff.jl +++ b/src/models/palazzi_puff.jl @@ -22,10 +22,10 @@ downwind_dispersion(x::Number, stab::Any, ::Type{Defaultσy}) = crosswind_disper Returns the solution to a short duration Gaussian puff dispersion model for the given scenario, based on the work of Palazzi *et al*. ```math -c\left(x,y,z,t\right) = \Chi\left(x,y,z\right) \frac{1}{2} \left( \mathrm{erf} \left( { {x - u (t-\Delta t)} \over \sqrt{2} \sigma_x } \right) - \mathrm{erf} \left( { {x - u t} \over \sqrt{2} \sigma_x } \right) \right) -```` +c\left(x,y,z,t\right) = \chi\left(x,y,z\right) \frac{1}{2} \left( \mathrm{erf} \left( { {x - u (t-\Delta t)} \over \sqrt{2} \sigma_x } \right) - \mathrm{erf} \left( { {x - u t} \over \sqrt{2} \sigma_x } \right) \right) +``` -where Χ is a Gaussian plume model and the σs are dispersion parameters. The `EquationSet` +where χ is a Gaussian plume model and the σs are dispersion parameters. The `EquationSet` defines the set of correlations used to calculate the dispersion parameters and windspeed. The concentration returned is in volume fraction, assuming the puff is a gas at ambient conditions. @@ -70,7 +70,7 @@ function (ps::PalazziSolution{F,S,E})(x,y,z,t) where {F<:Number,S<:StabilityClas return zero(F) end - Χ = ps.plume(x,y,z,t) + χ = ps.plume(x,y,z,t) Δt = min(t,ps.duration) u = ps.windspeed @@ -91,10 +91,10 @@ function (ps::PalazziSolution{F,S,E})(x,y,z,t) where {F<:Number,S<:StabilityClas σx_a = downwind_dispersion(xa, S, ps.equationset) σx_b = downwind_dispersion(xb, S, ps.equationset) a = (x-u*(t-Δt))/(√2*σx_a) - b = (x-u*t)/(√2*σx_b) + b = (x-u*t)/(√2*σx_b) # concentration - c = Χ*erf(b,a)/2 + c = χ*erf(b,a)/2 return min(c,one(F)) end \ No newline at end of file