From dfc23e8b41b32a4d5d48229b3ae95a15e9f3405b Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Sat, 28 Jun 2025 21:10:13 -0600 Subject: [PATCH 1/8] First pass at mixing layer Mixing layer model from ISC3 user guide, volume 2 equation 1-50. Calculates the full sum, without any checking to see if it is worth continuing. --- src/GasDispersion.jl | 3 +- src/models/gaussian_mixing_layer.jl | 234 ++++++++++++++++++++++++++++ src/utils/equation_sets/briggs.jl | 8 +- src/utils/equation_sets/isc3.jl | 4 +- src/utils/equation_sets/tno.jl | 10 +- src/utils/equation_sets/turner.jl | 4 +- src/utils/mixing_layer.jl | 3 + src/utils/pasquill_gifford.jl | 10 +- src/utils/utils.jl | 3 + 9 files changed, 260 insertions(+), 19 deletions(-) create mode 100644 src/models/gaussian_mixing_layer.jl create mode 100644 src/utils/mixing_layer.jl diff --git a/src/GasDispersion.jl b/src/GasDispersion.jl index 1486b49..b108fab 100644 --- a/src/GasDispersion.jl +++ b/src/GasDispersion.jl @@ -20,7 +20,7 @@ export SourceModel, JetSource # plume models export PlumeModel, Plume, plume -export GaussianPlume, SimpleJet, BritterMcQuaidPlume +export GaussianPlume, GaussianMixingLayer, SimpleJet, BritterMcQuaidPlume # puff models export PuffModel, Puff, puff @@ -90,6 +90,7 @@ plume(s; kwargs...) = plume(s, GaussianPlume; kwargs...) # plume models include("models/gaussian_plume.jl") +include("models/gaussian_mixing_layer.jl") include("models/simple_jet.jl") include("models/britter_mcquaid_plume.jl") diff --git a/src/models/gaussian_mixing_layer.jl b/src/models/gaussian_mixing_layer.jl new file mode 100644 index 0000000..c6f0b37 --- /dev/null +++ b/src/models/gaussian_mixing_layer.jl @@ -0,0 +1,234 @@ +# for dispatch +struct GaussianMixingLayer <: PlumeModel end + +# Solution to the gaussian plume +struct GaussianMixingLayerSolution{F<:Number,N<:Integer,P<:PlumeRise,S<:StabilityClass,E<:EquationSet} <: Plume + scenario::Scenario + model::Symbol + rate::F + max_concentration::F + mass_to_vol::F + windspeed::F + effective_stack_height::F + mixing_height::F + nterms::N + plumerise::P + stability::Type{S} + equationset::E +end + +@doc doc""" + plume(::Scenario, GaussianPlume[, ::EquationSet]; kwargs...) + +Returns the solution to a Gaussian plume dispersion model for the given scenario. + +```math +c\left(x,y,z\right) = {m_{i} \over { 2 \pi \sigma_{y} \sigma_{z} u } } +\exp \left[ -\frac{1}{2} \left( y \over \sigma_{y} \right)^2 \right] \\ +\times \left\{ \exp \left[ -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right] ++ \exp \left[ -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right] \right\} +``` + +where the σs are dispersion parameters correlated with the distance x. The +`EquationSet` defines the set of correlations used to calculate the dispersion +parameters. The concentration returned is in volume fraction, assuming the plume +is a gas at ambient conditions. + +# References ++ + +""" +function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet(); h_min=1.0, n_terms=10) + # parameters of the jet + hᵣ = _release_height(scenario) + m = _mass_rate(scenario) + Q = _release_flowrate(scenario) + ρⱼ = _release_density(scenario) + + # jet at ambient conditions + Tₐ = _atmosphere_temperature(scenario) + Pₐ = _atmosphere_pressure(scenario) + hₘ = _mixing_height(scenario) + ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ) + Qᵒ = Q*ρⱼ/ρₐ + + # max concentration + c_max = m/Qᵒ + c_max = c_max/ρₐ + + return GaussianMixingLayerSolution( + scenario, # scenario::Scenario + :gaussian, # model::Symbol + m, # mass emission rate + c_max, # max_concentration + ρₐ, # mass concentration to vol concentration + windspeed(scenario,max(hᵣ,h_min),eqs), # windspeed + hᵣ, # effective_stack_height::Number + hₘ, # height of mixing layer + n_terms, # number of terms in sum + NoPlumeRise(), # plume rise model + _stability(scenario), # stability class + eqs # equation set + ) +end + +@doc doc""" + plume(::Scenario{AbstractSubstance,VerticalJet,Atmosphere}, GaussianPlume[, ::EquationSet]; kwargs...) + +Returns the solution to a Gaussian plume dispersion model for a vertical jet. By default the Briggs +plume rise model is used. + +```math +c\left(x,y,z\right) = {m_{i} \over { 2 \pi \sigma_{y} \sigma_{z} u } } +\exp \left[ -\frac{1}{2} \left( y \over \sigma_{y} \right)^2 \right] \\ +\times \left\{ \exp \left[ -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right] ++ \exp \left[ -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right] \right\} +``` + +where the σs are dispersion parameters correlated with the distance x. The +`EquationSet` defines the set of correlations used to calculate the dispersion +parameters. The concentration returned is in volume fraction, assuming the plume +is a gas at ambient conditions. + +# Arguments +- `downwash::Bool=false`: when true, includes stack-downwash effects +- `plumerise::Bool=true`: when true, includes plume-rise effects using Briggs' model + +# References ++ AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers ++ Briggs, Gary A. 1969. *Plume Rise* Oak Ridge: U.S. Atomic Energy Commission + +""" +function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, ::Type{GaussianMixingLayer}, eqs=DefaultSet(); downwash::Bool=false, plumerise::Bool=true, h_min=1.0, n_terms=10) + # parameters of the jet + Dⱼ = _release_diameter(scenario) + uⱼ = _release_velocity(scenario) + hᵣ = _release_height(scenario) + m = _mass_rate(scenario) + Q = _release_flowrate(scenario) + ρⱼ = _release_density(scenario) + + # parameters of the environment + u = windspeed(scenario,max(hᵣ,h_min),eqs) + stab = _stability(scenario) + Γ = _lapse_rate(scenario) + hₘ = _mixing_height(scenario) + + # jet at ambient conditions + Tₐ = _atmosphere_temperature(scenario) + Pₐ = _atmosphere_pressure(scenario) + ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ) + Qᵒ = Q*ρⱼ/ρₐ + + # max concentration + c_max = m/Qᵒ + c_max = c_max/ρₐ + + # stack-tip downwash check + if (downwash==true) && (uⱼ < 1.5*u) + Δh_dw = 2*Dⱼ*( (uⱼ/u) - 1.5 ) + else + Δh_dw = 0.0 + end + + hᵣ = hᵣ + Δh_dw + + # plume rise + if plumerise == true + Tᵣ = _release_temperature(scenario) + plume = plume_rise(Dⱼ,uⱼ,Tᵣ,u,Tₐ,Γ,stab) + else + plume = NoPlumeRise() + end + + return GaussianMixingLayerSolution( + scenario, #scenario::Scenario + :gaussian, #model::Symbol + m, #mass emission rate + c_max, #max concentration + ρₐ, #mass to vol, at ambient conditions + u, #windspeed + hᵣ, #effective_stack_height::Number + hₘ, # height of mixing layer + n_terms, # number of terms in sum + plume, #plume rise model + stab, #stability class + eqs #equation set + ) +end + +function (g::GaussianMixingLayerSolution{F,N,NoPlumeRise,S,E})(x, y, z, t=0) where {F,N,S,E} + # domain check + h = g.effective_stack_height + if (x==0)&&(y==0)&&(z==h) + return g.max_concentration + elseif (x≤0)||(z<0) + return zero(F) + else + G = g.rate + u = g.windspeed + hₘ = g.mixing_height + σy = crosswind_dispersion(x,S,g.equationset) + σz = vertical_dispersion(x,S,g.equationset) + + Fy = exp(-0.5*(y/σy)^2)/(√(2π)*σy) + + Fz = exp(-0.5*((z-h)/σz)^2) + exp(-0.5*((z+h)/σz)^2) + for i=1:g.nterms + H₁ = z - (2*i*hₘ - h) + H₂ = z + (2*i*hₘ - h) + H₃ = z - (2*i*hₘ + h) + H₄ = z + (2*i*hₘ + h) + Fz += exp(-0.5*(H₁/σz)^2) + exp(-0.5*(H₂/σz)^2) + exp(-0.5*(H₃/σz)^2) + exp(-0.5*(H₄/σz)^2) + end + Fz /= √(2π)*σz + + c = (G/u)*Fy*Fz + # c is in kg/m^3 + # use density at ambient conditions to convert to vol pct + c_vol = c/g.mass_to_vol + + return min(c_vol,g.max_concentration) + end +end + +function (g::GaussianMixingLayerSolution{F,N,<:BriggsModel,S,E})(x, y, z, t=0) where {F,N,S,E} + # domain check + h = g.effective_stack_height + if (x==0)&&(y==0)&&(z==h) + return g.max_concentration + elseif (x≤0)||(z<0) + return zero(F) + else + G = g.rate + u = g.windspeed + h = g.effective_stack_height + hₘ = g.mixing_height + m = g.plumerise + Δh = plume_rise(x, m) + σy = crosswind_dispersion(x,S,g.equationset) + σz = vertical_dispersion(x,S,g.equationset) + hₑ = h + Δh + σyₑ = √( (Δh/3.5)^2 + σy^2 ) + σzₑ = √( (Δh/3.5)^2 + σz^2 ) + + Fy = exp(-0.5*(y/σyₑ)^2)/(√(2π)*σyₑ) + + Fz = exp(-0.5*((z-hₑ)/σzₑ)^2) + exp(-0.5*((z+hₑ)/σzₑ)^2) + for i=1:g.nterms + H₁ = z - (2*i*hₘ - hₑ) + H₂ = z + (2*i*hₘ - hₑ) + H₃ = z - (2*i*hₘ + hₑ) + H₄ = z + (2*i*hₘ + hₑ) + Fz += exp(-0.5*(H₁/σzₑ)^2) + exp(-0.5*(H₂/σzₑ)^2) + exp(-0.5*(H₃/σzₑ)^2) + exp(-0.5*(H₄/σzₑ)^2) + end + Fz /= √(2π)*σzₑ + + c = (G/u)*Fy*Fz + # c is in kg/m^3 + # use density at ambient conditions to convert to vol pct + c_vol = c/g.mass_to_vol + + return min(c_vol,g.max_concentration) + end +end diff --git a/src/utils/equation_sets/briggs.jl b/src/utils/equation_sets/briggs.jl index d6d3cd5..c6fbf9f 100644 --- a/src/utils/equation_sets/briggs.jl +++ b/src/utils/equation_sets/briggs.jl @@ -6,7 +6,7 @@ struct BriggsRuralσy <: DispersionFunction end struct BriggsRuralσz <: DispersionFunction end """ - crosswind_dispersion(x, Plume, StabilityClass, Briggsσy) + crosswind_dispersion(x, StabilityClass, Briggsσy) Plume crosswind dispersion correlations, for rural terrain @@ -24,7 +24,7 @@ crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{BriggsRuralσy}) = 0.04x/ """ - vertical_dispersion(x, Plume, StabilityClass, CCPSRural) + vertical_dispersion(x, StabilityClass, BriggsRuralσz) Plume vertical dispersion correlations, for rural terrain @@ -46,7 +46,7 @@ struct BriggsUrbanσy <: DispersionFunction end struct BriggsUrbanσz <: DispersionFunction end """ - crosswind_dispersion(x, Plume, StabilityClass, CCPSUrban) + crosswind_dispersion(x, StabilityClass, BriggsUrbanσy) Plume crosswind dispersion correlations, for urban terrain @@ -64,7 +64,7 @@ crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{BriggsUrbanσy}) = 0.11x/ """ - vertical_dispersion(x, Plume, StabilityClass, CCPSUrban) + vertical_dispersion(x, StabilityClass, BriggsUrbanσz) Plume vertical dispersion correlations, for urban terrain diff --git a/src/utils/equation_sets/isc3.jl b/src/utils/equation_sets/isc3.jl index 88d8e7e..e5fbc2b 100644 --- a/src/utils/equation_sets/isc3.jl +++ b/src/utils/equation_sets/isc3.jl @@ -30,7 +30,7 @@ ISC3Rural = BasicEquationSet{IrwinRural,Nothing,ISC3Ruralσy,ISC3Ruralσz} """ - crosswind_dispersion(x, Plume, StabilityClass, ISC3Rural) + crosswind_dispersion(x, StabilityClass, ISC3Ruralσy) Plume crosswind dispersion correlations, for rural terrain @@ -46,7 +46,7 @@ crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{ISC3Ruralσy}) = 465.1162 """ - vertical_dispersion(x, Plume, StabilityClass, ISC3Rural) + vertical_dispersion(x, StabilityClass, ISC3Ruralσz) Plume vertical dispersion correlations, for rural terrain diff --git a/src/utils/equation_sets/tno.jl b/src/utils/equation_sets/tno.jl index 9e9b1b2..28cdc8c 100644 --- a/src/utils/equation_sets/tno.jl +++ b/src/utils/equation_sets/tno.jl @@ -65,7 +65,7 @@ function windspeed(z::Number, u::Number, zR::Number, L::Number, ::Union{Type{Cla end """ - crosswind_dispersion(x, Plume, StabilityClass, TNO) + crosswind_dispersion(x, StabilityClass, TNOPlumeσy) Plume crosswind dispersion correlations @@ -80,7 +80,7 @@ crosswind_dispersion(x::Number, ::Type{ClassE}, ::Type{TNOPlumeσy}) = 0.098x^0. crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{TNOPlumeσy}) = 0.065x^0.902 """ - vertical_dispersion(x, Plume, StabilityClass, TNO) + vertical_dispersion(x, StabilityClass, TNOPlumeσz) Plume vertical dispersion correlations @@ -101,7 +101,7 @@ struct TNOPuffσz <: DispersionFunction end TNOPuff = BasicEquationSet{TNOWind,TNOPuffσx,TNOPuffσy,TNOPuffσz} """ - crosswind_dispersion(x, Puff, StabilityClass, TNO) + crosswind_dispersion(x, StabilityClass, TNOPuffσy) Puff crosswind dispersion correlations @@ -112,7 +112,7 @@ crosswind_dispersion(x::Number, stability::Any, ::Type{TNOPuffσy}) = 0.5*crossw """ - vertical_dispersion(x, Puff, StabilityClass, TNO) + vertical_dispersion(x, StabilityClass, TNOPuffσz) Puff vertical dispersion correlations @@ -127,7 +127,7 @@ vertical_dispersion(x::Number, ::Type{ClassE}, ::Type{TNOPuffσz}) = 0.15x vertical_dispersion(x::Number, ::Type{ClassF}, ::Type{TNOPuffσz}) = 0.12x """ - downwind_dispersion(x, Puff, StabilityClass, TNO) + downwind_dispersion(x, StabilityClass, TNOPuffσx) Puff downwind dispersion correlations diff --git a/src/utils/equation_sets/turner.jl b/src/utils/equation_sets/turner.jl index b19a077..c8cbf8d 100644 --- a/src/utils/equation_sets/turner.jl +++ b/src/utils/equation_sets/turner.jl @@ -4,7 +4,7 @@ struct Turnerσz <: DispersionFunction end Turner = BasicEquationSet{DefaultWind,Nothing,Turnerσy,Turnerσz} """ - crosswind_dispersion(x, Plume, StabilityClass, Turner) + crosswind_dispersion(x, StabilityClass, Turnerσy) Plume crosswind dispersion correlations @@ -20,7 +20,7 @@ crosswind_dispersion(x::Number, ::Type{ClassE}, ::Type{Turnerσy}) = 0.091x^0.91 crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{Turnerσy}) = 0.067x^0.90 """ - vertical_dispersion(x, Plume, StabilityClass, Turner) + vertical_dispersion(x, StabilityClass, Turnerσz) Plume vertical dispersion correlations diff --git a/src/utils/mixing_layer.jl b/src/utils/mixing_layer.jl new file mode 100644 index 0000000..9da8fe1 --- /dev/null +++ b/src/utils/mixing_layer.jl @@ -0,0 +1,3 @@ +_mixing_height(s::Scenario) = _mixing_height(s.atmosphere) + +_mixing_height(a::SimpleAtmosphere) = Inf # just for a test \ No newline at end of file diff --git a/src/utils/pasquill_gifford.jl b/src/utils/pasquill_gifford.jl index 5e1ecd7..ee3a667 100644 --- a/src/utils/pasquill_gifford.jl +++ b/src/utils/pasquill_gifford.jl @@ -6,7 +6,7 @@ crosswind_dispersion(x::Number, s::Any, ::BasicEquationSet{W,SX,SY,SZ}) where {W vertical_dispersion(x::Number, s::Any, ::BasicEquationSet{W,SX,SY,SZ}) where {W,SX,SY,SZ<:DispersionFunction} = vertical_dispersion(x,s,SZ) """ - crosswind_dispersion(x, StabilityClass; avg_time=600.0) + crosswind_dispersion(x, StabilityClass, Defaultσy; avg_time=600.0) Plume crosswind dispersion correlations @@ -52,7 +52,7 @@ end """ - vertical_dispersion(x, Plume, StabilityClass) + vertical_dispersion(x, StabilityClass, Defaultσz) Plume vertical dispersion correlations @@ -106,7 +106,7 @@ struct CCPSPuffσy <: DispersionFunction end struct CCPSPuffσz <: DispersionFunction end """ - crosswind_dispersion(x, Puff, StabilityClass) + crosswind_dispersion(x, StabilityClass, CCPSPuffσy) Puff crosswind dispersion correlations @@ -122,7 +122,7 @@ crosswind_dispersion(x::Number, ::Type{ClassF}, ::Type{CCPSPuffσy}) = 0.02*x^0. """ - downwind_dispersion(x, Puff, StabilityClass) + downwind_dispersion(x, StabilityClass, CCPSPuffσx) Puff downwind dispersion correlations @@ -135,7 +135,7 @@ end """ - vertical_dispersion(x, Puff, StabilityClass) + vertical_dispersion(x, StabilityClass, CCPSPuffσz) Puff vertical dispersion correlations diff --git a/src/utils/utils.jl b/src/utils/utils.jl index 43ead0c..f07f7f8 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -10,6 +10,9 @@ include("pasquill_gifford.jl") # Briggs' model for plume rise include("plume_rise.jl") +# mixing height correlations +include("mixing_layer.jl") + # model correlations include("britter_mcquaid_correls.jl") From f76a2152fbffe85bae6cdaab699b3102b28263c4 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Sun, 29 Jun 2025 09:06:17 -0600 Subject: [PATCH 2/8] Refactor Refactor simple Gaussian plume to be more re-useable, then the mixing layer is just an extension of the existing plume solution --- src/models/gaussian_mixing_layer.jl | 126 +++++++--------------------- src/models/gaussian_plume.jl | 37 +++++--- test/models/gaussian_plume_tests.jl | 2 +- 3 files changed, 56 insertions(+), 109 deletions(-) diff --git a/src/models/gaussian_mixing_layer.jl b/src/models/gaussian_mixing_layer.jl index c6f0b37..101be08 100644 --- a/src/models/gaussian_mixing_layer.jl +++ b/src/models/gaussian_mixing_layer.jl @@ -1,24 +1,30 @@ # for dispatch struct GaussianMixingLayer <: PlumeModel end -# Solution to the gaussian plume -struct GaussianMixingLayerSolution{F<:Number,N<:Integer,P<:PlumeRise,S<:StabilityClass,E<:EquationSet} <: Plume - scenario::Scenario - model::Symbol - rate::F - max_concentration::F - mass_to_vol::F - windspeed::F - effective_stack_height::F - mixing_height::F +# new vertical term +struct SimpleMixingLayer{N<:Integer,F<:Number} <: GaussianVerticalTerm nterms::N - plumerise::P - stability::Type{S} - equationset::E + mixing_height::F +end + +function vertical_term(z, h, σz, ml::SimpleMixingLayer) + Fz = exp(-0.5*((z-h)/σz)^2) + exp(-0.5*((z+h)/σz)^2) + for i=1:ml.nterms + H₁ = z - (2*i*ml.mixing_height - h) + H₂ = z + (2*i*ml.mixing_height - h) + H₃ = z - (2*i*ml.mixing_height + h) + H₄ = z + (2*i*ml.mixing_height + h) + next_term = exp(-0.5*(H₁/σz)^2) + exp(-0.5*(H₂/σz)^2) + exp(-0.5*(H₃/σz)^2) + exp(-0.5*(H₄/σz)^2) + Fz += next_term + if next_term ≈ 0 + break + end + end + return Fz/(√(2π)*σz) end @doc doc""" - plume(::Scenario, GaussianPlume[, ::EquationSet]; kwargs...) + plume(::Scenario, GaussianMixingLayer[, ::EquationSet]; kwargs...) Returns the solution to a Gaussian plume dispersion model for the given scenario. @@ -56,16 +62,16 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() c_max = m/Qᵒ c_max = c_max/ρₐ - return GaussianMixingLayerSolution( + return GaussianPlumeSolution( scenario, # scenario::Scenario - :gaussian, # model::Symbol + :simplemixinglayer, # model::Symbol m, # mass emission rate c_max, # max_concentration ρₐ, # mass concentration to vol concentration windspeed(scenario,max(hᵣ,h_min),eqs), # windspeed hᵣ, # effective_stack_height::Number - hₘ, # height of mixing layer - n_terms, # number of terms in sum + SimpleCrossTerm(), + SimpleMixingLayer(n_terms,hₘ), NoPlumeRise(), # plume rise model _stability(scenario), # stability class eqs # equation set @@ -141,94 +147,18 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere plume = NoPlumeRise() end - return GaussianMixingLayerSolution( + return GaussianPlumeSolution( scenario, #scenario::Scenario - :gaussian, #model::Symbol + :simplemixinglayer, #model::Symbol m, #mass emission rate c_max, #max concentration ρₐ, #mass to vol, at ambient conditions u, #windspeed hᵣ, #effective_stack_height::Number - hₘ, # height of mixing layer - n_terms, # number of terms in sum + SimpleCrossTerm(), + SimpleMixingLayer(n_terms,hₘ), plume, #plume rise model stab, #stability class eqs #equation set ) end - -function (g::GaussianMixingLayerSolution{F,N,NoPlumeRise,S,E})(x, y, z, t=0) where {F,N,S,E} - # domain check - h = g.effective_stack_height - if (x==0)&&(y==0)&&(z==h) - return g.max_concentration - elseif (x≤0)||(z<0) - return zero(F) - else - G = g.rate - u = g.windspeed - hₘ = g.mixing_height - σy = crosswind_dispersion(x,S,g.equationset) - σz = vertical_dispersion(x,S,g.equationset) - - Fy = exp(-0.5*(y/σy)^2)/(√(2π)*σy) - - Fz = exp(-0.5*((z-h)/σz)^2) + exp(-0.5*((z+h)/σz)^2) - for i=1:g.nterms - H₁ = z - (2*i*hₘ - h) - H₂ = z + (2*i*hₘ - h) - H₃ = z - (2*i*hₘ + h) - H₄ = z + (2*i*hₘ + h) - Fz += exp(-0.5*(H₁/σz)^2) + exp(-0.5*(H₂/σz)^2) + exp(-0.5*(H₃/σz)^2) + exp(-0.5*(H₄/σz)^2) - end - Fz /= √(2π)*σz - - c = (G/u)*Fy*Fz - # c is in kg/m^3 - # use density at ambient conditions to convert to vol pct - c_vol = c/g.mass_to_vol - - return min(c_vol,g.max_concentration) - end -end - -function (g::GaussianMixingLayerSolution{F,N,<:BriggsModel,S,E})(x, y, z, t=0) where {F,N,S,E} - # domain check - h = g.effective_stack_height - if (x==0)&&(y==0)&&(z==h) - return g.max_concentration - elseif (x≤0)||(z<0) - return zero(F) - else - G = g.rate - u = g.windspeed - h = g.effective_stack_height - hₘ = g.mixing_height - m = g.plumerise - Δh = plume_rise(x, m) - σy = crosswind_dispersion(x,S,g.equationset) - σz = vertical_dispersion(x,S,g.equationset) - hₑ = h + Δh - σyₑ = √( (Δh/3.5)^2 + σy^2 ) - σzₑ = √( (Δh/3.5)^2 + σz^2 ) - - Fy = exp(-0.5*(y/σyₑ)^2)/(√(2π)*σyₑ) - - Fz = exp(-0.5*((z-hₑ)/σzₑ)^2) + exp(-0.5*((z+hₑ)/σzₑ)^2) - for i=1:g.nterms - H₁ = z - (2*i*hₘ - hₑ) - H₂ = z + (2*i*hₘ - hₑ) - H₃ = z - (2*i*hₘ + hₑ) - H₄ = z + (2*i*hₘ + hₑ) - Fz += exp(-0.5*(H₁/σzₑ)^2) + exp(-0.5*(H₂/σzₑ)^2) + exp(-0.5*(H₃/σzₑ)^2) + exp(-0.5*(H₄/σzₑ)^2) - end - Fz /= √(2π)*σzₑ - - c = (G/u)*Fy*Fz - # c is in kg/m^3 - # use density at ambient conditions to convert to vol pct - c_vol = c/g.mass_to_vol - - return min(c_vol,g.max_concentration) - end -end diff --git a/src/models/gaussian_plume.jl b/src/models/gaussian_plume.jl index bd64e05..0e2f41f 100644 --- a/src/models/gaussian_plume.jl +++ b/src/models/gaussian_plume.jl @@ -1,8 +1,12 @@ +# for modularity and code re-use +abstract type GaussianCrossTerm end +abstract type GaussianVerticalTerm end + # for dispatch struct GaussianPlume <: PlumeModel end # Solution to the gaussian plume -struct GaussianPlumeSolution{F<:Number,P<:PlumeRise,S<:StabilityClass,E<:EquationSet} <: Plume +struct GaussianPlumeSolution{F<:Number,C<:GaussianCrossTerm,V<:GaussianVerticalTerm,P<:PlumeRise,S<:StabilityClass,E<:EquationSet} <: Plume scenario::Scenario model::Symbol rate::F @@ -10,11 +14,20 @@ struct GaussianPlumeSolution{F<:Number,P<:PlumeRise,S<:StabilityClass,E<:Equatio mass_to_vol::F windspeed::F effective_stack_height::F + crossterm::C + verticalterm::V plumerise::P stability::Type{S} equationset::E end -GaussianPlumeSolution(s,m,Q,c,ρ,u,h_eff,pr,stab,es) = GaussianPlumeSolution(s,m,promote(Q,c,ρ,u,h_eff)...,pr,stab,es) +GaussianPlumeSolution(s,m,Q,c,ρ,u,h_eff,cross,vert,pr,stab,es) = GaussianPlumeSolution(s,m,promote(Q,c,ρ,u,h_eff)...,cross,vert,pr,stab,es) + +struct SimpleCrossTerm <: GaussianCrossTerm end +cross_term(y, σy, ::SimpleCrossTerm) = exp(-0.5*(y/σy)^2)/(√(2π)*σy) + +struct SimpleVerticalTerm <: GaussianVerticalTerm end +vertical_term(z, h, σz, ::SimpleVerticalTerm) = ( exp(-0.5*((z-h)/σz)^2) + exp(-0.5*((z+h)/σz)^2) )/(√(2π)*σz) + @doc doc""" plume(::Scenario, GaussianPlume[, ::EquationSet]; kwargs...) @@ -62,6 +75,8 @@ function plume(scenario::Scenario, ::Type{GaussianPlume}, eqs=DefaultSet(); h_mi ρₐ, # mass concentration to vol concentration windspeed(scenario,max(hᵣ,h_min),eqs), # windspeed hᵣ, # effective_stack_height::Number + SimpleCrossTerm(), + SimpleVerticalTerm(), NoPlumeRise(), # plume rise model _stability(scenario), # stability class eqs # equation set @@ -144,13 +159,15 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere ρₐ, #mass to vol, at ambient conditions u, #windspeed hᵣ, #effective_stack_height::Number + SimpleCrossTerm(), + SimpleVerticalTerm(), plume, #plume rise model stab, #stability class eqs #equation set ) end -function (g::GaussianPlumeSolution{F,NoPlumeRise,S,E})(x, y, z, t=0) where {F<:Number,S<:StabilityClass,E<:EquationSet} +function (g::GaussianPlumeSolution{F,C,V,NoPlumeRise,S,E})(x, y, z, t=0) where {F,C,V,S,E} # domain check h = g.effective_stack_height if (x==0)&&(y==0)&&(z==h) @@ -163,9 +180,9 @@ function (g::GaussianPlumeSolution{F,NoPlumeRise,S,E})(x, y, z, t=0) where {F<:N σy = crosswind_dispersion(x,S,g.equationset) σz = vertical_dispersion(x,S,g.equationset) - c = ( G/(2*π*u*σy*σz) - * exp(-0.5*(y/σy)^2) - * ( exp(-0.5*((z-h)/σz)^2) + exp(-0.5*((z+h)/σz)^2) ) ) + Fy = cross_term(y,σy,g.crossterm) + Fz = vertical_term(z,h,σz,g.verticalterm) + c = (G/u)*Fy*Fz # c is in kg/m^3 # use density at ambient conditions to convert to vol pct @@ -175,7 +192,7 @@ function (g::GaussianPlumeSolution{F,NoPlumeRise,S,E})(x, y, z, t=0) where {F<:N end end -function (g::GaussianPlumeSolution{F,<:BriggsModel,S,E})(x, y, z, t=0) where {F<:Number,S<:StabilityClass,E<:EquationSet} +function (g::GaussianPlumeSolution{F,C,V,<:BriggsModel,S,E})(x, y, z, t=0) where {F,C,V,S,E} # domain check h = g.effective_stack_height if (x==0)&&(y==0)&&(z==h) @@ -194,9 +211,9 @@ function (g::GaussianPlumeSolution{F,<:BriggsModel,S,E})(x, y, z, t=0) where {F< σyₑ = √( (Δh/3.5)^2 + σy^2 ) σzₑ = √( (Δh/3.5)^2 + σz^2 ) - c = ( G/(2*π*u*σyₑ*σzₑ) - * exp(-0.5*(y/σyₑ)^2) - * ( exp(-0.5*((z-hₑ)/σzₑ)^2) + exp(-0.5*((z+hₑ)/σzₑ)^2) ) ) + Fy = cross_term(y,σyₑ,g.crossterm) + Fz = vertical_term(z,hₑ,σzₑ,g.verticalterm) + c = (G/u)*Fy*Fz # c is in kg/m^3 # use density at ambient conditions to convert to vol pct diff --git a/test/models/gaussian_plume_tests.jl b/test/models/gaussian_plume_tests.jl index 1aa32f1..e957a7a 100644 --- a/test/models/gaussian_plume_tests.jl +++ b/test/models/gaussian_plume_tests.jl @@ -30,7 +30,7 @@ pl = plume(scn) # test default behaviour and type inheritance - @test GasDispersion.GaussianPlumeSolution(scn,:test,1.0,2,3,4,5,GasDispersion.NoPlumeRise(),ClassA,DefaultSet()) isa GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, Nothing, GasDispersion.Defaultσy, GasDispersion.Defaultσz}} + @test GasDispersion.GaussianPlumeSolution(scn,:test,1.0,2,3,4,5,GasDispersion.SimpleCrossTerm(),GasDispersion.SimpleVerticalTerm(),GasDispersion.NoPlumeRise(),ClassA,DefaultSet()) isa GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, Nothing, GasDispersion.Defaultσy, GasDispersion.Defaultσz}} @test pl isa GasDispersion.GaussianPlumeSolution @test pl isa Plume @test pl(-1,0,0) == 0.0 From 34270c753e681464ea1d26e5e5995d2d1f6b3e88 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Sun, 29 Jun 2025 13:13:51 -0600 Subject: [PATCH 3/8] Unit tests and some refactoring SimpleAtmospheres have a default mixing height calculation, which is unfortunately a function of latitude. I used the default baked into TSCREEN and ISC3 (40\deg N) but this is problematic. --- docs/src/equation_sets.md | 4 +- docs/src/plume.md | 52 +++++++++++----- src/models/gaussian_mixing_layer.jl | 97 +++++++++++++++++++---------- src/utils/mixing_layer.jl | 15 ++++- test/models/gaussian_ml_tests.jl | 69 ++++++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 186 insertions(+), 52 deletions(-) create mode 100644 test/models/gaussian_ml_tests.jl diff --git a/docs/src/equation_sets.md b/docs/src/equation_sets.md index 3669ed5..0d9e6a6 100644 --- a/docs/src/equation_sets.md +++ b/docs/src/equation_sets.md @@ -112,7 +112,7 @@ SimpleAtmosphere atmosphere: The plume using the default equation set is simply this -```jldoctest eqnset_example; output = false, filter = r"(\d*)\.(\d{4})\d+" => s"\1.\2***" + All of these plumes can then be plotted, to better visualize what is going on. These are identical plume models with the only differences being the windspeed correlation and the dispersion correlations. diff --git a/docs/src/plume.md b/docs/src/plume.md index a399f28..4861c08 100644 --- a/docs/src/plume.md +++ b/docs/src/plume.md @@ -7,13 +7,15 @@ independent, this includes, for example, emissions from elevated stacks. plume ``` -## Gaussian Plumes +## Gaussian Plume Models + +### Simple Gaussian Plumes ```@docs plume(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{GaussianPlume}) ``` -The gaussian plume model assumes that concentration profile in the crosswind (y) and vertical (z) directions follow gaussian distributions with dispersions $\sigma_y$ and $\sigma_z$, respectively. This model can be derived from an advection-diffusion equation, or simply taken as a given. +The simple gaussian plume model assumes that concentration profile in the crosswind (y) and vertical (z) directions follow gaussian distributions with dispersions $\sigma_y$ and $\sigma_z$, respectively. This model can be derived from an advection-diffusion equation, or simply taken as a given. The basic gaussian would have the plume expand downward beyond the ground, to correct for this an additional term for *ground reflection* is added. This is equivalent to adding a mirror image source reflected across the x-z plane, and causes mass to accumulate along the ground instead of simply disappearing (as would happen in the naive case). @@ -37,7 +39,7 @@ with Three important parameters are determined from correlations, which are functions of the atmospheric stability: the windspeed at the release point, the crosswind dispersion, and the vertical dispersion. The model converts the final concentration to volume fraction, assuming the plume is a gas at ambient conditions. -### Crosswind dispersion correlations +#### Crosswind dispersion correlations The crosswind dispersion, $\sigma_{y}$ is a function of downwind distance, $x$ as well as stability class @@ -59,7 +61,7 @@ Where $\delta$, $\beta$, and $\gamma$ are tabulated based on stability class ([S | F | 0.0674 | 0.9 | -### Vertical dispersion correlations +#### Vertical dispersion correlations The vertical dispersion, $\sigma_{z}$ is a function of downwind distance, $x$ as well as stability class @@ -80,7 +82,7 @@ Where $\delta$ and $\beta$ are tabulated based on stability class ([Seinfeld 198 | E | 0.02275 | 1.3010 | -0.0450 | | F | 0.01122 | 1.4024 | -0.0540 | -### 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$. This scenario is adapted from CCPS *Guidelines for Consequence Analysis of Chemical Releases*([AIChE/CCPS 1999](references.md), 47) @@ -149,15 +151,6 @@ And then pass it to the `plume` function ```jldoctest gaussplume; output = false, filter = r"(\d*)\.(\d{4})\d+" => s"\1.\2***" g = plume(scn, GaussianPlume) -# output - -GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}()) - -``` - -Where `g` is a callable which returns the concentration (in vol fraction) at any point. For example at 100m downwind and at a height of 2m - -```jldoctest gaussplume; output = true, filter = r"(\d*)\.(\d{4})\d+" => s"\1.\2***" g(100,0,2) # output @@ -166,7 +159,7 @@ g(100,0,2) ``` -Which is ~612ppm (vol). Beyond simply having a number, we may want a plan-view of the plume at a given height, say 2m. +Where `g` is a callable which returns the concentration (in vol fraction) at any point. For example at 100m downwind and at a height of 2m the result is ~612ppm (vol). Beyond simply having a number, we may want a plan-view of the plume at a given height, say 2m. ```@setup gaussplume using ..GasDispersion @@ -217,6 +210,35 @@ plot(g, xlims=(0,50), ylims=(-10,10), height=3.5, clims=(0,LEL), aspect_ratio=:equal) ``` +### Gaussian Plumes with Mixing Layers + +```@docs +plume(::Scenario, ::Type{GaussianMixingLayer}) +``` + +The gaussian mixing layer model allows for the plume to be contained entirely within a mixing layer of a given height. This can be done using either the method of images or a periodic boundary. + +```math +c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z +``` + +Where $F_z$ is the vertical dispersion term. For the method of images this takes the form of the infinite sum: + +```math +F_z = {\left( \exp \left[ -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right] ++ \exp \left[ -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right] \right) \over {\sqrt{2\pi} \sigma_z} } +\sum_{i=1}^{n} \left[ \exp\left(-\frac{1}{2}\left(\frac{z - h + 2i h_m}{\sigma_z}\right)^2\right) ++ \exp\left(-\frac{1}{2}\left(\frac{z - h - 2i h_m}{\sigma_z}\right)^2\right) ++ \exp\left(-\frac{1}{2}\left(\frac{z + h + 2i h_m}{\sigma_z}\right)^2\right) ++ \exp\left(-\frac{1}{2}\left(\frac{z + h - 2i h_m}{\sigma_z}\right)^2\right) \right] +``` + +For the periodic boundary it takes the form of another infinite sum: + +```math +F_z = ... +``` + ## Simple Jet Plumes ```@docs diff --git a/src/models/gaussian_mixing_layer.jl b/src/models/gaussian_mixing_layer.jl index 101be08..0924506 100644 --- a/src/models/gaussian_mixing_layer.jl +++ b/src/models/gaussian_mixing_layer.jl @@ -24,27 +24,33 @@ function vertical_term(z, h, σz, ml::SimpleMixingLayer) end @doc doc""" - plume(::Scenario, GaussianMixingLayer[, ::EquationSet]; kwargs...) + plume(::Scenario, GaussianMixingLayer[, ::EquationSet]; h_min=1.0, n_terms=10) -Returns the solution to a Gaussian plume dispersion model for the given scenario. +Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer. ```math -c\left(x,y,z\right) = {m_{i} \over { 2 \pi \sigma_{y} \sigma_{z} u } } -\exp \left[ -\frac{1}{2} \left( y \over \sigma_{y} \right)^2 \right] \\ -\times \left\{ \exp \left[ -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right] -+ \exp \left[ -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right] \right\} +c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z ``` -where the σs are dispersion parameters correlated with the distance x. The -`EquationSet` defines the set of correlations used to calculate the dispersion -parameters. The concentration returned is in volume fraction, assuming the plume -is a gas at ambient conditions. +where \(F_z\) is the vertical dispersion term, a function of the mixing height \(h_m\), \(n\) is the number of image terms, and other symbols are as defined for a Gaussian +plume model. + +# Keyword Arguments +- `h_min=1.0`: Minimum height for windspeed calculations. +- `method=:simplemixinglayer`: The method used for the mixing layer. +- `n_terms=10`: Number of image terms for the mixing layer reflection. +- `mixing_limit=10_000.0`: Limit for the mixing height, in meters. Mixing heights above this are treated as infinite. # References -+ ++ AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers ++ US EPA. 1995. *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models EPA-454/B-95-003b, vol 2*. Research Triangle Park, NC: Office of Air Quality Planning and Standards """ -function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet(); h_min=1.0, n_terms=10) +function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet(); + method=:simplemixinglayer, + h_min=1.0, + n_terms=10, + mixing_limit=10_000.0) # parameters of the jet hᵣ = _release_height(scenario) m = _mass_rate(scenario) @@ -54,7 +60,7 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() # jet at ambient conditions Tₐ = _atmosphere_temperature(scenario) Pₐ = _atmosphere_pressure(scenario) - hₘ = _mixing_height(scenario) + hₘ = _mixing_height(scenario, eqs) ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ) Qᵒ = Q*ρⱼ/ρₐ @@ -62,6 +68,16 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() c_max = m/Qᵒ c_max = c_max/ρₐ + # mixing layer + if hₘ > mixing_limit + # infinite mixing height + ml = SimpleVerticalTerm() + elseif method == :simplemixinglayer + ml = SimpleMixingLayer(n_terms, hₘ) + else + error("Unknown mixing layer method: $method") + end + return GaussianPlumeSolution( scenario, # scenario::Scenario :simplemixinglayer, # model::Symbol @@ -71,7 +87,7 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() windspeed(scenario,max(hᵣ,h_min),eqs), # windspeed hᵣ, # effective_stack_height::Number SimpleCrossTerm(), - SimpleMixingLayer(n_terms,hₘ), + ml, NoPlumeRise(), # plume rise model _stability(scenario), # stability class eqs # equation set @@ -79,33 +95,38 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() end @doc doc""" - plume(::Scenario{AbstractSubstance,VerticalJet,Atmosphere}, GaussianPlume[, ::EquationSet]; kwargs...) + plume(::Scenario{AbstractSubstance,VerticalJet,Atmosphere}, GaussianPlume[, ::EquationSet]; **kwargs...) -Returns the solution to a Gaussian plume dispersion model for a vertical jet. By default the Briggs -plume rise model is used. +Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer. ```math -c\left(x,y,z\right) = {m_{i} \over { 2 \pi \sigma_{y} \sigma_{z} u } } -\exp \left[ -\frac{1}{2} \left( y \over \sigma_{y} \right)^2 \right] \\ -\times \left\{ \exp \left[ -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right] -+ \exp \left[ -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right] \right\} +c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z ``` -where the σs are dispersion parameters correlated with the distance x. The -`EquationSet` defines the set of correlations used to calculate the dispersion -parameters. The concentration returned is in volume fraction, assuming the plume -is a gas at ambient conditions. +where \(F_z\) is the vertical dispersion term, a function of the mixing height \(h_m\), \(n\) is the number of image terms, and other symbols are as defined for a Gaussian +plume model. -# Arguments -- `downwash::Bool=false`: when true, includes stack-downwash effects -- `plumerise::Bool=true`: when true, includes plume-rise effects using Briggs' model +# Keyword Arguments +- `downwash::Bool=false`: Include stack-downwash effects if true. +- `plumerise::Bool=true`: Include plume-rise effects using Briggs' model if true. +- `h_min=1.0`: Minimum height, in meters, for windspeed calculations. +- `method=:simplemixinglayer`: The method used for the mixing layer. +- `n_terms=10`: Number of image terms for the mixing layer reflection. +- `mixing_limit=10_000.0`: Limit for the mixing height, in meters. Mixing heights above this are treated as infinite. # References -+ AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers -+ Briggs, Gary A. 1969. *Plume Rise* Oak Ridge: U.S. Atomic Energy Commission +- AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers +- Briggs, Gary A. 1969. *Plume Rise* Oak Ridge: U.S. Atomic Energy Commission +- US EPA. 1995. *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models EPA-454/B-95-003b, vol 2*. Research Triangle Park, NC: Office of Air Quality Planning and Standards """ -function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, ::Type{GaussianMixingLayer}, eqs=DefaultSet(); downwash::Bool=false, plumerise::Bool=true, h_min=1.0, n_terms=10) +function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, + ::Type{GaussianMixingLayer}, eqs=DefaultSet(); + downwash::Bool=false, plumerise::Bool=true, + method=:simplemixinglayer, + h_min=1.0, + n_terms=10, + mixing_limit=10_000.0) # parameters of the jet Dⱼ = _release_diameter(scenario) uⱼ = _release_velocity(scenario) @@ -118,7 +139,7 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere u = windspeed(scenario,max(hᵣ,h_min),eqs) stab = _stability(scenario) Γ = _lapse_rate(scenario) - hₘ = _mixing_height(scenario) + hₘ = _mixing_height(scenario, eqs) # jet at ambient conditions Tₐ = _atmosphere_temperature(scenario) @@ -147,6 +168,16 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere plume = NoPlumeRise() end + # mixing layer + if hₘ > mixing_limit + # infinite mixing height + ml = SimpleVerticalTerm() + elseif method == :simplemixinglayer + ml = SimpleMixingLayer(n_terms, hₘ) + else + error("Unknown mixing layer method: $method") + end + return GaussianPlumeSolution( scenario, #scenario::Scenario :simplemixinglayer, #model::Symbol @@ -156,7 +187,7 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere u, #windspeed hᵣ, #effective_stack_height::Number SimpleCrossTerm(), - SimpleMixingLayer(n_terms,hₘ), + ml, plume, #plume rise model stab, #stability class eqs #equation set diff --git a/src/utils/mixing_layer.jl b/src/utils/mixing_layer.jl index 9da8fe1..e99d999 100644 --- a/src/utils/mixing_layer.jl +++ b/src/utils/mixing_layer.jl @@ -1,3 +1,14 @@ -_mixing_height(s::Scenario) = _mixing_height(s.atmosphere) +const ω = 7.2921e-5 # Earth's angular velocity in rad/s -_mixing_height(a::SimpleAtmosphere) = Inf # just for a test \ No newline at end of file +_mixing_height(s::Scenario, eqs::BasicEquationSet) = _mixing_height(s.atmosphere, eqs) + +function _mixing_height(a::SimpleAtmosphere, es::BasicEquationSet; ϕ=40*π/180) + stab = _stability(a) + if stab ∈ (ClassA, ClassB, ClassC, ClassD) + f = 2ω*sin(ϕ) # Coriolis parameter + u_str = friction_velocity(a, es) + return 0.3*u_str/f + else + return Inf # infinite mixing height for stable conditions + end +end \ No newline at end of file diff --git a/test/models/gaussian_ml_tests.jl b/test/models/gaussian_ml_tests.jl new file mode 100644 index 0000000..5be4631 --- /dev/null +++ b/test/models/gaussian_ml_tests.jl @@ -0,0 +1,69 @@ +@testset "Gaussian mixing layer tests" begin + # Gaussian plume example, *Guidelines for Consequence Analysis of Chemical + # Releases* CCPS, 1999, pg 97 + 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=Inf, + diameter=1, + velocity=1, + height=0, + temperature=298., + pressure=101325., + fraction_liquid=0) + atm1 = SimpleAtmosphere(temperature=298, + pressure=101325, + windspeed=2, + windspeed_height=1, + stability=ClassF) + + atm2 = SimpleAtmosphere(temperature=298, + pressure=101325, + windspeed=2, + windspeed_height=1, + stability=ClassD) + scn = Scenario(sub,rel,atm1) + + # test default behaviour and type inheritance + pl1 = plume(scn, GaussianMixingLayer) + @test pl1.verticalterm isa GasDispersion.SimpleVerticalTerm + + scn2 = Scenario(sub,rel,atm2) + pl2 = plume(scn2, GaussianMixingLayer) + @test pl2.verticalterm isa GasDispersion.SimpleMixingLayer + + # test vertical jet with mixing layer + rel = VerticalJet(mass_rate=0.1, + duration=Inf, + diameter=1, + velocity=1, + height=10, + temperature=298., + pressure=101325., + fraction_liquid=0) + scn3 = Scenario(sub,rel,atm1) + pl3 = plume(scn3, GaussianMixingLayer) + @test pl3.verticalterm isa GasDispersion.SimpleVerticalTerm + + atm3 = SimpleAtmosphere(temperature=298, + pressure=101325, + windspeed=2, + windspeed_height=10, + stability=ClassD) + scn4 = Scenario(sub,rel,atm3) + pl4 = plume(scn4, GaussianMixingLayer; downwash=true, plumerise=true) + @test pl4.verticalterm isa GasDispersion.SimpleMixingLayer + @test pl4.verticalterm.mixing_height ≈ 640.0311954829524 + @test pl4(500, 0, 0) ≈ 1.7194314353343084e-5 + +end diff --git a/test/runtests.jl b/test/runtests.jl index 223e87c..402b063 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,6 +30,7 @@ if GROUP == "All" || GROUP == "Model" # testing dispersion models include("models/gaussian_plume_tests.jl") + include("models/gaussian_ml_tests.jl") include("models/gaussian_puff_tests.jl") include("models/intpuff_tests.jl") include("models/palazzi_tests.jl") From c7b532d25b7b6498bf6b65f4633443aaf25b247b Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Sun, 29 Jun 2025 13:57:11 -0600 Subject: [PATCH 4/8] Updated documentation Updated the docs to refer to the Gaussian Mixing Layer model --- docs/src/equation_sets.md | 18 ++++++++-------- docs/src/plume.md | 32 +++++++++++++++++------------ docs/src/puff.md | 5 +++-- src/models/gaussian_mixing_layer.jl | 8 ++++---- src/models/gaussian_plume.jl | 2 +- src/models/palazzi_puff.jl | 2 +- src/models/simple_jet.jl | 23 +++++++++++++++++++++ src/models/slab_puff.jl | 15 ++++++++++++++ 8 files changed, 75 insertions(+), 30 deletions(-) diff --git a/docs/src/equation_sets.md b/docs/src/equation_sets.md index 0d9e6a6..eb6b0e4 100644 --- a/docs/src/equation_sets.md +++ b/docs/src/equation_sets.md @@ -112,12 +112,12 @@ SimpleAtmosphere atmosphere: The plume using the default equation set is simply this - +``` All of these plumes can then be plotted, to better visualize what is going on. These are identical plume models with the only differences being the windspeed correlation and the dispersion correlations. diff --git a/docs/src/plume.md b/docs/src/plume.md index 4861c08..e5919c0 100644 --- a/docs/src/plume.md +++ b/docs/src/plume.md @@ -9,7 +9,7 @@ plume ## Gaussian Plume Models -### Simple Gaussian Plumes +### Simple Gaussian Plume ```@docs plume(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{GaussianPlume}) @@ -210,26 +210,26 @@ plot(g, xlims=(0,50), ylims=(-10,10), height=3.5, clims=(0,LEL), aspect_ratio=:equal) ``` -### Gaussian Plumes with Mixing Layers +### Gaussian Plume within a Mixing Layer ```@docs -plume(::Scenario, ::Type{GaussianMixingLayer}) +plume(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{GaussianMixingLayer}) ``` The gaussian mixing layer model allows for the plume to be contained entirely within a mixing layer of a given height. This can be done using either the method of images or a periodic boundary. ```math -c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z +c(x, y, z) = \frac{m_i}{u} \frac{1}{\sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z ``` Where $F_z$ is the vertical dispersion term. For the method of images this takes the form of the infinite sum: ```math -F_z = {\left( \exp \left[ -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right] -+ \exp \left[ -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right] \right) \over {\sqrt{2\pi} \sigma_z} } +F_z = \frac{1}{\sqrt{2\pi} \sigma_z} \left( \exp \left( -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right) ++ \exp \left( -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right) \right) \\ \sum_{i=1}^{n} \left[ \exp\left(-\frac{1}{2}\left(\frac{z - h + 2i h_m}{\sigma_z}\right)^2\right) -+ \exp\left(-\frac{1}{2}\left(\frac{z - h - 2i h_m}{\sigma_z}\right)^2\right) -+ \exp\left(-\frac{1}{2}\left(\frac{z + h + 2i h_m}{\sigma_z}\right)^2\right) ++ \exp\left(-\frac{1}{2}\left(\frac{z - h - 2i h_m}{\sigma_z}\right)^2\right) \\ ++ \exp\left(-\frac{1}{2}\left(\frac{z + h + 2i h_m}{\sigma_z}\right)^2\right) \\ + \exp\left(-\frac{1}{2}\left(\frac{z + h - 2i h_m}{\sigma_z}\right)^2\right) \right] ``` @@ -239,10 +239,15 @@ For the periodic boundary it takes the form of another infinite sum: F_z = ... ``` -## Simple Jet Plumes +As `SimpleAtmosphere`s do not define mixing height, one is calculated based on the following: +- for stable atmospheres (class E and F) the mixing height is assumed to be infinite, and the model defaults back to a [Simple Gaussian Plume](@ref) +- for unstable and neutral atmospheres (classes A through D) the mixing height is calculated from the friction velocity $u^{*}$ and the coriolis parameter $f$: $h_m = 0.3 \frac{u^{*}}{f}$ where $f$ is calculated at a default latitude of 40°N (consistent with [US EPA 1995](references.md)). + +## Jet Plumes +### Simple Jet Model ```@docs -plume(::Scenario, ::Type{SimpleJet}) +plume(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{SimpleJet}) ``` Simple jet dispersion models are a useful tool for evaluating dispersion near the region where a jet release is occurring. They are based on a simplified model where the air is stationary and all of the momentum needed to mix the release is supplied by the jet. This is in some ways the opposite assumptions than are used in the Gaussian Plume model -- where the release is assumed to have negligible velocity and the momentum is entirely supplied by the wind. @@ -259,7 +264,7 @@ with + $\rho_j$ - initial density of the jet material, kg/m^3 + $\rho_a$ - density of the ambient atmosphere, kg/m^3 -### Model Parameters +#### Model Parameters The model parameters $k_2$ and $k_3$ are per [Long (1963)](references.md) @@ -274,7 +279,7 @@ the initial concentration is calculated from the mass flowrate and volumetric fl c_0 = { Q_i \over Q } = { \dot{m} \over \rho Q } = { \dot{m} \over { \rho \frac{\pi}{4} d^2 u } } ``` -### Example +#### Example Suppose we wish to model the dispersion of gaseous propane using the same scenario, `scn`, worked out above. @@ -294,7 +299,8 @@ plot(j, xlims=(0,100), ylims=(-10,10), height=2) ``` -## Britter-McQuaid Model +## Top-Hat Models +### Britter-McQuaid Model ```@docs plume(::Scenario, ::Type{BritterMcQuaidPlume}) diff --git a/docs/src/puff.md b/docs/src/puff.md index aa20ee1..0835996 100644 --- a/docs/src/puff.md +++ b/docs/src/puff.md @@ -19,7 +19,8 @@ A simple gaussian puff model assumes the release is instantaneous, and all mass ```math c_{puff} = { m_i \over { (2 \pi)^{3/2} \sigma_x \sigma_y \sigma_z } } \exp \left( -\frac{1}{2} \left( {x - ut} \over \sigma_x \right)^2 \right) -\exp \left( -\frac{1}{2} \left( {y} \over \sigma_y \right)^2 \right) \left[ \exp \left( -\frac{1}{2} \left( {z - h} \over \sigma_z \right)^2 \right) +\exp \left( -\frac{1}{2} \left( {y} \over \sigma_y \right)^2 \right) \\ +\left[ \exp \left( -\frac{1}{2} \left( {z - h} \over \sigma_z \right)^2 \right) + \exp \left( -\frac{1}{2} \left( {z + h} \over \sigma_z \right)^2 \right)\right] ``` @@ -320,7 +321,7 @@ in the provided equationset are for windspeed. ### SLAB Jet Model ```@docs -puff(::Scenario, ::Type{SLAB}) +puff(::Scenario{Substance,VerticalJet,Atmosphere}, ::Type{SLAB}) ``` The SLAB jet model is derived from the SLAB software package developed by diff --git a/src/models/gaussian_mixing_layer.jl b/src/models/gaussian_mixing_layer.jl index 0924506..fe86085 100644 --- a/src/models/gaussian_mixing_layer.jl +++ b/src/models/gaussian_mixing_layer.jl @@ -24,7 +24,7 @@ function vertical_term(z, h, σz, ml::SimpleMixingLayer) end @doc doc""" - plume(::Scenario, GaussianMixingLayer[, ::EquationSet]; h_min=1.0, n_terms=10) + plume(::Scenario, GaussianMixingLayer[, ::EquationSet]; **kwargs...) Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer. @@ -32,7 +32,7 @@ Returns the solution to a Gaussian plume dispersion model with a simple reflecti c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z ``` -where \(F_z\) is the vertical dispersion term, a function of the mixing height \(h_m\), \(n\) is the number of image terms, and other symbols are as defined for a Gaussian +where $ Fz $ is the vertical dispersion term, a function of the mixing height $ h_m $, *n* is the number of image terms, and other symbols are as defined for a Gaussian plume model. # Keyword Arguments @@ -95,7 +95,7 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() end @doc doc""" - plume(::Scenario{AbstractSubstance,VerticalJet,Atmosphere}, GaussianPlume[, ::EquationSet]; **kwargs...) + plume(::Scenario{Substance,VerticalJet,Atmosphere}, GaussianMixingLayer[, ::EquationSet]; **kwargs...) Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer. @@ -103,7 +103,7 @@ Returns the solution to a Gaussian plume dispersion model with a simple reflecti c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z ``` -where \(F_z\) is the vertical dispersion term, a function of the mixing height \(h_m\), \(n\) is the number of image terms, and other symbols are as defined for a Gaussian +where $ Fz $ is the vertical dispersion term, a function of the mixing height $ h_m $, *n* is the number of image terms, and other symbols are as defined for a Gaussian plume model. # Keyword Arguments diff --git a/src/models/gaussian_plume.jl b/src/models/gaussian_plume.jl index 0e2f41f..68664c7 100644 --- a/src/models/gaussian_plume.jl +++ b/src/models/gaussian_plume.jl @@ -84,7 +84,7 @@ function plume(scenario::Scenario, ::Type{GaussianPlume}, eqs=DefaultSet(); h_mi end @doc doc""" - plume(::Scenario{AbstractSubstance,VerticalJet,Atmosphere}, GaussianPlume[, ::EquationSet]; kwargs...) + plume(::Scenario{Substance,VerticalJet,Atmosphere}, GaussianPlume[, ::EquationSet]; kwargs...) Returns the solution to a Gaussian plume dispersion model for a vertical jet. By default the Briggs plume rise model is used. diff --git a/src/models/palazzi_puff.jl b/src/models/palazzi_puff.jl index 73a17c4..3e46f65 100644 --- a/src/models/palazzi_puff.jl +++ b/src/models/palazzi_puff.jl @@ -37,7 +37,7 @@ downwind dispersion, σx, is calculated: - `: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$ +- `plume_model::Type{Plume} = GaussianPlume` : the plume model $\chi$ - `disp_method = :default` : the method for calculating the downwind dispersion # References diff --git a/src/models/simple_jet.jl b/src/models/simple_jet.jl index 286186a..2695fbf 100644 --- a/src/models/simple_jet.jl +++ b/src/models/simple_jet.jl @@ -65,6 +65,29 @@ function plume(scenario::Scenario, ::Type{SimpleJet}, eqs::EquationSet=DefaultSe end +@doc doc""" + plume(::Scenario{Substance,VerticalJet,Atmosphere}, SimpleJet; kwargs...) + +Returns the solution to a simple turbulent jet dispersion model for a vertical jet. + +```math +c\left(x,y,z\right) = k_2 c_0 \left( d \over z \right) \sqrt{ \rho_j \over \rho_a } +\exp \left( - \left( k_3 { r \over z } \right)^2 \right) +``` + +where *r* is the radial distance from the jet centerline. Assumes a circular jet +with diameter equal to the jet diameter. Ground-reflection is included by method +of images. + +# References ++ Long, V.D. 1963. "Estimation of the Extent of Hazard Areas Round a Vent." *Chem. Process Hazard*. II:6 + +# Arguments +- `release_angle::Number=π/2`: the angle at which the jet is released, in radians +- `k2::Number=6` parameter of the model, default value is recommended by Long +- `k3::Number=5` parameter of the model, default value is recommended by Long + +""" function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, ::Type{SimpleJet}, eqs::EquationSet=DefaultSet(); release_angle::Number=π/2, k2::Number=6.0, k3::Number=5.0) # Density correction ρj = _release_density(scenario) diff --git a/src/models/slab_puff.jl b/src/models/slab_puff.jl index 0f74ff0..517ab0a 100644 --- a/src/models/slab_puff.jl +++ b/src/models/slab_puff.jl @@ -128,6 +128,21 @@ function puff(scenario::Scenario, ::Type{SLAB}, eqs::EquationSet=DefaultSet(); AkimaInterpolation(out.cc.betax[tperm], out.cc.t[tperm])) end + +@doc doc""" + puff(::Scenario{Substance,VerticalJet,Atmosphere}, SLAB; kwargs...) + +Returns the solution to the SLAB vertical jet dispersion model for the given +scenario. + +# References ++ Ermak, Donald L. 1990. *User's Manual for SLAB: An Atmospheric Dispersion Model For Denser-Than-Air Releases* Lawrence Livermore National Laboratory + +# Arguments +- `t_av::Number=10`: averaging time, seconds +- `x_max::Number=2000`: maximum downwind distance, meters, this defines the problem domain + +""" function puff(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, ::Type{SLAB}, eqs::EquationSet=DefaultSet(); t_av=10, x_max=2000) c_max = 1.0 From 25ac4108cf277a67e42f7b1377655f290f27d571 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Mon, 30 Jun 2025 14:42:26 -0600 Subject: [PATCH 5/8] Added periodic mixing layer Added mixing layer model based on an infinite sum of cosines. This converges very slowly (in my testing) but is used in some places. --- docs/src/plume.md | 16 ++--- docs/src/references.md | 4 +- src/models/gaussian_mixing_layer.jl | 55 +++++++++++++++--- src/models/palazzi_puff.jl | 7 ++- test/models/gaussian_ml_tests.jl | 90 ++++++++++++++++++----------- 5 files changed, 118 insertions(+), 54 deletions(-) diff --git a/docs/src/plume.md b/docs/src/plume.md index e5919c0..975ecd3 100644 --- a/docs/src/plume.md +++ b/docs/src/plume.md @@ -222,23 +222,25 @@ The gaussian mixing layer model allows for the plume to be contained entirely wi c(x, y, z) = \frac{m_i}{u} \frac{1}{\sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z ``` -Where $F_z$ is the vertical dispersion term. For the method of images this takes the form of the infinite sum: +Where $F_z$ is the vertical dispersion term. For the method of images this takes the form of the infinite sum ([Beychok 1994](references.md) p 123) ```math F_z = \frac{1}{\sqrt{2\pi} \sigma_z} \left( \exp \left( -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right) + \exp \left( -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right) \right) \\ -\sum_{i=1}^{n} \left[ \exp\left(-\frac{1}{2}\left(\frac{z - h + 2i h_m}{\sigma_z}\right)^2\right) -+ \exp\left(-\frac{1}{2}\left(\frac{z - h - 2i h_m}{\sigma_z}\right)^2\right) \\ -+ \exp\left(-\frac{1}{2}\left(\frac{z + h + 2i h_m}{\sigma_z}\right)^2\right) \\ -+ \exp\left(-\frac{1}{2}\left(\frac{z + h - 2i h_m}{\sigma_z}\right)^2\right) \right] +\sum_{n=1}^{\infty} \left[ \exp\left(-\frac{1}{2}\left(\frac{z - h + 2n h_m}{\sigma_z}\right)^2\right) ++ \exp\left(-\frac{1}{2}\left(\frac{z - h - 2n h_m}{\sigma_z}\right)^2\right) \\ ++ \exp\left(-\frac{1}{2}\left(\frac{z + h + 2n h_m}{\sigma_z}\right)^2\right) \\ ++ \exp\left(-\frac{1}{2}\left(\frac{z + h - 2n h_m}{\sigma_z}\right)^2\right) \right] ``` -For the periodic boundary it takes the form of another infinite sum: +For the periodic boundary it takes the form of another infinite sum ([Seinfeld and Pandis 2006](references.md) p 858) ```math -F_z = ... +F_z = \frac{2}{h_m} \left( \frac{1}{2} + \sum_{n=1}^{\infty} cos\left( {n\pi z} \over h_m \right) cos\left( {n\pi h} \over h_m \right) \exp\left(-\frac{1}{2}\left(\frac{n\pi \sigma_z}{h_m}\right)^2\right) \right) ``` +The periodic boundary layer approach can be very slow to converge, in which case the number of terms given by the keyword argument `n_terms` should be increased. By default the sums terminate early if the solution converges, e.g. if `n_terms=100_000_000` or some other huge number but the sum converges after 5 terms, only 5 terms will be calculated. The simple mixing layer approach converges much more quickly, and typically 10 terms are more than enough. + As `SimpleAtmosphere`s do not define mixing height, one is calculated based on the following: - for stable atmospheres (class E and F) the mixing height is assumed to be infinite, and the model defaults back to a [Simple Gaussian Plume](@ref) - for unstable and neutral atmospheres (classes A through D) the mixing height is calculated from the friction velocity $u^{*}$ and the coriolis parameter $f$: $h_m = 0.3 \frac{u^{*}}{f}$ where $f$ is calculated at a default latitude of 40°N (consistent with [US EPA 1995](references.md)). diff --git a/docs/src/references.md b/docs/src/references.md index 7ed086d..d5fbe3b 100644 --- a/docs/src/references.md +++ b/docs/src/references.md @@ -2,6 +2,8 @@ AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers +Beychok, Milton R. 1994. *Fundamentals of Stack Gas Dispersion* 3rd Ed. Irvine, CA: Milton R. Beychok. + Bakkum, E.A. and N.J. Duijm. 2005. "Chapter 4 - Vapour Cloud Dispersion" in *Methods for the Calculation of Physical Effects, CPR 14E* (TNO Yellow Book) Edited by C.J.H. van den Bosch and R.A.P.M. Weterings. The Netherlands. Businger, J. A., J. C. Wyngaard, Y. Izumi, and E. F. Bradley. 1971. "Flux-Profile Relationships in the Atmospheric Surfaace Layer." *Journal of the Atmospheric Sciences*. 28, 181-189. [doi:10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2](https://doi.org/10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2) @@ -30,7 +32,7 @@ Paulson, C. A. 1970. "The Mathematical Representation of Wind Speed and Temperat Spicer, Thomas O. and Jerry A. Havens. 1988. *Development of Vapor Dispersion Models for Non-Neutrally Buoyant Gas Mixtures--Analysis of TFI/NH3 Test Data*. Tyndall Air Force Base, FL: USAF Engineering & Services Laboratory -Seinfeld, John H. 1986. *Atmospheric Chemistry and Physics of Air Pollution*. New York: John Wiley and Sons +Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons. Turner, D. Bruce. 1970. *Workbook of Atmospheric Dispersion Estimates*. United States Environmental Protection Agency. diff --git a/src/models/gaussian_mixing_layer.jl b/src/models/gaussian_mixing_layer.jl index fe86085..8f9d18f 100644 --- a/src/models/gaussian_mixing_layer.jl +++ b/src/models/gaussian_mixing_layer.jl @@ -1,7 +1,7 @@ # for dispatch struct GaussianMixingLayer <: PlumeModel end -# new vertical term +# mixing layer vertical term -- method of images struct SimpleMixingLayer{N<:Integer,F<:Number} <: GaussianVerticalTerm nterms::N mixing_height::F @@ -23,27 +23,50 @@ function vertical_term(z, h, σz, ml::SimpleMixingLayer) return Fz/(√(2π)*σz) end +# mixing layer vertical term -- periodic +struct PeriodicMixingLayer{N<:Integer,F<:Number} <: GaussianVerticalTerm + nterms::N + mixing_height::F +end + +function vertical_term(z, h, σz, ml::PeriodicMixingLayer) + Fz = 1/2 + for n=1:ml.nterms + next_term = cos(n*π*z/ml.mixing_height)*cos(n*π*h/ml.mixing_height)*exp(-0.5*(n*π*σz/ml.mixing_height)^2) + Fz += next_term + if next_term ≈ 0 + break + end + end + return 2*Fz/ml.mixing_height +end + @doc doc""" plume(::Scenario, GaussianMixingLayer[, ::EquationSet]; **kwargs...) Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer. ```math -c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z +c(x, y, z) = \frac{m_i}{u} \frac{1}{ \sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z ``` -where $ Fz $ is the vertical dispersion term, a function of the mixing height $ h_m $, *n* is the number of image terms, and other symbols are as defined for a Gaussian +where $ F\_z $ is the vertical dispersion term, a function of the mixing height $ h\_m $, and other symbols are as defined for a Gaussian plume model. +There are two methods for the mixing layer: +- `:simplemixinglayer` uses a simple method of images, a series of reflections off of the mixing height +- `:periodicmixinglayer` uses an infinite series of cosine terms to calculate the vertical dispersion, which is more accurate for large mixing heights + # Keyword Arguments - `h_min=1.0`: Minimum height for windspeed calculations. - `method=:simplemixinglayer`: The method used for the mixing layer. -- `n_terms=10`: Number of image terms for the mixing layer reflection. +- `n_terms=10`: Number terms in the series calculation of $ F_z $. - `mixing_limit=10_000.0`: Limit for the mixing height, in meters. Mixing heights above this are treated as infinite. # References + AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers + US EPA. 1995. *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models EPA-454/B-95-003b, vol 2*. Research Triangle Park, NC: Office of Air Quality Planning and Standards ++ Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons. """ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet(); @@ -63,6 +86,9 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() hₘ = _mixing_height(scenario, eqs) ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ) Qᵒ = Q*ρⱼ/ρₐ + if hᵣ>hₘ + error("Release height $hᵣ exceeds mixing height $hₘ. This is not supported by the Gaussian mixing layer model.") + end # max concentration c_max = m/Qᵒ @@ -71,9 +97,12 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() # mixing layer if hₘ > mixing_limit # infinite mixing height + @warn "Mixing height $hₘ exceeds limit $mixing_limit, using infinite mixing height." ml = SimpleVerticalTerm() elseif method == :simplemixinglayer ml = SimpleMixingLayer(n_terms, hₘ) + elseif method == :periodicmixinglayer + ml = PeriodicMixingLayer(n_terms, hₘ) else error("Unknown mixing layer method: $method") end @@ -100,24 +129,28 @@ end Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer. ```math -c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z +c(x, y, z) = \frac{m_i}{u} \frac{1}{ \sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z ``` - -where $ Fz $ is the vertical dispersion term, a function of the mixing height $ h_m $, *n* is the number of image terms, and other symbols are as defined for a Gaussian +where $ F\_z $ is the vertical dispersion term, a function of the mixing height $ h\_m $, and other symbols are as defined for a Gaussian plume model. +There are two methods for the mixing layer: +- `:simplemixinglayer` uses a simple method of images, a series of reflections off of the mixing height +- `:periodicmixinglayer` uses an infinite series of cosine terms to calculate the vertical dispersion, which is more accurate for large mixing heights + # Keyword Arguments - `downwash::Bool=false`: Include stack-downwash effects if true. - `plumerise::Bool=true`: Include plume-rise effects using Briggs' model if true. - `h_min=1.0`: Minimum height, in meters, for windspeed calculations. - `method=:simplemixinglayer`: The method used for the mixing layer. -- `n_terms=10`: Number of image terms for the mixing layer reflection. +- `n_terms=10`: Number terms in the series calculation of $ F_z $. - `mixing_limit=10_000.0`: Limit for the mixing height, in meters. Mixing heights above this are treated as infinite. # References - AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers - Briggs, Gary A. 1969. *Plume Rise* Oak Ridge: U.S. Atomic Energy Commission - US EPA. 1995. *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models EPA-454/B-95-003b, vol 2*. Research Triangle Park, NC: Office of Air Quality Planning and Standards +- Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons. """ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, @@ -140,6 +173,9 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere stab = _stability(scenario) Γ = _lapse_rate(scenario) hₘ = _mixing_height(scenario, eqs) + if hᵣ>hₘ + error("Release height $hᵣ exceeds mixing height $hₘ. This is not supported by the Gaussian mixing layer model.") + end # jet at ambient conditions Tₐ = _atmosphere_temperature(scenario) @@ -171,9 +207,12 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere # mixing layer if hₘ > mixing_limit # infinite mixing height + @warn "Mixing height $hₘ exceeds limit $mixing_limit, using infinite mixing height." ml = SimpleVerticalTerm() elseif method == :simplemixinglayer ml = SimpleMixingLayer(n_terms, hₘ) + elseif method == :periodicmixinglayer + ml = PeriodicMixingLayer(n_terms, hₘ) else error("Unknown mixing layer method: $method") end diff --git a/src/models/palazzi_puff.jl b/src/models/palazzi_puff.jl index 3e46f65..aea98f7 100644 --- a/src/models/palazzi_puff.jl +++ b/src/models/palazzi_puff.jl @@ -37,14 +37,15 @@ downwind dispersion, σx, is calculated: - `: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 +- `plume_model::Type{Plume} = GaussianPlume` : the plume model $\chi$ +- additional keyword arguments are passed to the plume model # 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=PalazziDefaultSet(); plume_model=GaussianPlume, disp_method=:default) +function puff(scenario::Scenario, ::Type{Palazzi}, eqs=PalazziDefaultSet(); plume_model=GaussianPlume, disp_method=:default, kwargs...) stab = _stability(scenario) Δt = _duration(scenario) @@ -55,7 +56,7 @@ function puff(scenario::Scenario, ::Type{Palazzi}, eqs=PalazziDefaultSet(); plum scenario, #scenario::Scenario :palazzi, #model::Symbol disp_method, #dispersion model::Symbol - plume(scenario, plume_model, eqs), # plume model + plume(scenario, plume_model, eqs; kwargs...), # plume model Δt, # duration u, # windspeed stab, # stability class diff --git a/test/models/gaussian_ml_tests.jl b/test/models/gaussian_ml_tests.jl index 5be4631..f927513 100644 --- a/test/models/gaussian_ml_tests.jl +++ b/test/models/gaussian_ml_tests.jl @@ -13,7 +13,7 @@ latent_heat=1, gas_heat_capacity=1, liquid_heat_capacity=1) - rel = HorizontalJet(mass_rate=0.1, + hj = HorizontalJet(mass_rate=0.1, duration=Inf, diameter=1, velocity=1, @@ -21,29 +21,17 @@ temperature=298., pressure=101325., fraction_liquid=0) - atm1 = SimpleAtmosphere(temperature=298, - pressure=101325, - windspeed=2, - windspeed_height=1, - stability=ClassF) - atm2 = SimpleAtmosphere(temperature=298, - pressure=101325, - windspeed=2, - windspeed_height=1, - stability=ClassD) - scn = Scenario(sub,rel,atm1) - - # test default behaviour and type inheritance - pl1 = plume(scn, GaussianMixingLayer) - @test pl1.verticalterm isa GasDispersion.SimpleVerticalTerm + hj_toohigh = HorizontalJet(mass_rate=0.1, + duration=Inf, + diameter=1, + velocity=1, + height=10_000, # too high for the mixing layer + temperature=298., + pressure=101325., + fraction_liquid=0) - scn2 = Scenario(sub,rel,atm2) - pl2 = plume(scn2, GaussianMixingLayer) - @test pl2.verticalterm isa GasDispersion.SimpleMixingLayer - - # test vertical jet with mixing layer - rel = VerticalJet(mass_rate=0.1, + vj = VerticalJet(mass_rate=0.1, duration=Inf, diameter=1, velocity=1, @@ -51,19 +39,51 @@ temperature=298., pressure=101325., fraction_liquid=0) - scn3 = Scenario(sub,rel,atm1) - pl3 = plume(scn3, GaussianMixingLayer) - @test pl3.verticalterm isa GasDispersion.SimpleVerticalTerm + + vj_toohigh = VerticalJet(mass_rate=0.1, + duration=Inf, + diameter=1, + velocity=1, + height=10_000, # too high for the mixing layer + temperature=298., + pressure=101325., + fraction_liquid=0) + + stbl = SimpleAtmosphere(temperature=298, + pressure=101325, + windspeed=2, + windspeed_height=10, + stability=ClassF) + + neut = SimpleAtmosphere(temperature=298, + pressure=101325, + windspeed=2, + windspeed_height=10, + stability=ClassD) + + + # horizontal jet test default behaviour and type inheritance + @test plume(Scenario(sub,hj,stbl), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleVerticalTerm + @test plume(Scenario(sub,hj,neut), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleMixingLayer + @test plume(Scenario(sub,hj,neut), GaussianMixingLayer; method=:periodicmixinglayer).verticalterm isa GasDispersion.PeriodicMixingLayer + @test_throws ErrorException plume(Scenario(sub,hj,neut), GaussianMixingLayer; method=:someothermethod) + @test_throws ErrorException plume(Scenario(sub,hj_toohigh,neut), GaussianMixingLayer) - atm3 = SimpleAtmosphere(temperature=298, - pressure=101325, - windspeed=2, - windspeed_height=10, - stability=ClassD) - scn4 = Scenario(sub,rel,atm3) - pl4 = plume(scn4, GaussianMixingLayer; downwash=true, plumerise=true) - @test pl4.verticalterm isa GasDispersion.SimpleMixingLayer - @test pl4.verticalterm.mixing_height ≈ 640.0311954829524 - @test pl4(500, 0, 0) ≈ 1.7194314353343084e-5 + # vertical jet test default behaviour and type inheritance + @test plume(Scenario(sub,vj,stbl), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleVerticalTerm + @test plume(Scenario(sub,vj,neut), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleMixingLayer + @test plume(Scenario(sub,vj,neut), GaussianMixingLayer; method=:periodicmixinglayer).verticalterm isa GasDispersion.PeriodicMixingLayer + @test_throws ErrorException plume(Scenario(sub,vj,neut), GaussianMixingLayer; method=:someothermethod) + @test_throws ErrorException plume(Scenario(sub,vj_toohigh,neut), GaussianMixingLayer) + + # integration test with a scenario -- simple mixing layer + pl_smpl = plume(Scenario(sub,vj,neut), GaussianMixingLayer; downwash=true, plumerise=true) + @test pl_smpl.verticalterm.mixing_height ≈ 640.0311954829524 + @test pl_smpl(500, 0, 0) ≈ 1.7194314353343084e-5 + + # integration test with a scenario -- periodic mixing layer + pl_per = plume(Scenario(sub,vj,neut), GaussianMixingLayer; downwash=true, plumerise=true, method=:periodicmixinglayer, n_terms=100_000_000) + @test pl_per.verticalterm.mixing_height ≈ 640.0311954829524 + @test pl_per(500, 0, 0) ≈ 1.7194314353343084e-5 end From 6be5e72cd33604608d21360391a9908a08d001aa Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Mon, 30 Jun 2025 16:24:16 -0600 Subject: [PATCH 6/8] Corrected domain issues Added a `ProblemDomain` type and a domain check function the generalize the GaussianPlumeSolution. Previously it was ignoring the mixing layer when deciding on the domain of the problem. Added a limiter to plume rise, forcing plume rise below the mixing layer. --- src/models/gaussian_mixing_layer.jl | 18 ++++++++++++++---- src/models/gaussian_plume.jl | 27 +++++++++++++++++---------- src/utils/plume_rise.jl | 10 ++++++++++ src/utils/utils.jl | 18 +++++++++++++++++- test/models/gaussian_plume_tests.jl | 2 +- 5 files changed, 59 insertions(+), 16 deletions(-) diff --git a/src/models/gaussian_mixing_layer.jl b/src/models/gaussian_mixing_layer.jl index 8f9d18f..24e24f3 100644 --- a/src/models/gaussian_mixing_layer.jl +++ b/src/models/gaussian_mixing_layer.jl @@ -107,6 +107,9 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() error("Unknown mixing layer method: $method") end + # domain + dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, hₘ) + return GaussianPlumeSolution( scenario, # scenario::Scenario :simplemixinglayer, # model::Symbol @@ -119,8 +122,8 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet() ml, NoPlumeRise(), # plume rise model _stability(scenario), # stability class - eqs # equation set - ) + eqs, # equation set + dom) # problem domain end @doc doc""" @@ -200,6 +203,10 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere if plumerise == true Tᵣ = _release_temperature(scenario) plume = plume_rise(Dⱼ,uⱼ,Tᵣ,u,Tₐ,Γ,stab) + if plume.final_rise > hₘ-hᵣ + @warn "Plume rise $(plume.final_rise) exceeds mixing height $hₘ, plume rise is limited to the mixing height." + plume = _new_plume_rise(plume, hₘ) + end else plume = NoPlumeRise() end @@ -217,6 +224,9 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere error("Unknown mixing layer method: $method") end + # domain + dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, hₘ) + return GaussianPlumeSolution( scenario, #scenario::Scenario :simplemixinglayer, #model::Symbol @@ -229,6 +239,6 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere ml, plume, #plume rise model stab, #stability class - eqs #equation set - ) + eqs, #equation set + dom) #problem domain end diff --git a/src/models/gaussian_plume.jl b/src/models/gaussian_plume.jl index 68664c7..699441e 100644 --- a/src/models/gaussian_plume.jl +++ b/src/models/gaussian_plume.jl @@ -6,7 +6,7 @@ abstract type GaussianVerticalTerm end struct GaussianPlume <: PlumeModel end # Solution to the gaussian plume -struct GaussianPlumeSolution{F<:Number,C<:GaussianCrossTerm,V<:GaussianVerticalTerm,P<:PlumeRise,S<:StabilityClass,E<:EquationSet} <: Plume +struct GaussianPlumeSolution{F<:Number,C<:GaussianCrossTerm,V<:GaussianVerticalTerm,P<:PlumeRise,S<:StabilityClass,E<:EquationSet,D<:ProblemDomain} <: Plume scenario::Scenario model::Symbol rate::F @@ -19,8 +19,9 @@ struct GaussianPlumeSolution{F<:Number,C<:GaussianCrossTerm,V<:GaussianVerticalT plumerise::P stability::Type{S} equationset::E + domain::D end -GaussianPlumeSolution(s,m,Q,c,ρ,u,h_eff,cross,vert,pr,stab,es) = GaussianPlumeSolution(s,m,promote(Q,c,ρ,u,h_eff)...,cross,vert,pr,stab,es) +GaussianPlumeSolution(s,m,Q,c,ρ,u,h_eff,cross,vert,pr,stab,es,dom) = GaussianPlumeSolution(s,m,promote(Q,c,ρ,u,h_eff)...,cross,vert,pr,stab,es,dom) struct SimpleCrossTerm <: GaussianCrossTerm end cross_term(y, σy, ::SimpleCrossTerm) = exp(-0.5*(y/σy)^2)/(√(2π)*σy) @@ -67,20 +68,23 @@ function plume(scenario::Scenario, ::Type{GaussianPlume}, eqs=DefaultSet(); h_mi c_max = m/Qᵒ c_max = c_max/ρₐ + # domain + dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, Inf) + return GaussianPlumeSolution( scenario, # scenario::Scenario :gaussian, # model::Symbol m, # mass emission rate c_max, # max_concentration - ρₐ, # mass concentration to vol concentration - windspeed(scenario,max(hᵣ,h_min),eqs), # windspeed + ρₐ, # mass concentration to vol concentration + windspeed(scenario,max(hᵣ,h_min),eqs), # windspeed hᵣ, # effective_stack_height::Number SimpleCrossTerm(), SimpleVerticalTerm(), NoPlumeRise(), # plume rise model _stability(scenario), # stability class - eqs # equation set - ) + eqs, # equation set + dom) # problem domain end @doc doc""" @@ -151,6 +155,9 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere plume = NoPlumeRise() end + # domain + dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, Inf) + return GaussianPlumeSolution( scenario, #scenario::Scenario :gaussian, #model::Symbol @@ -163,8 +170,8 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere SimpleVerticalTerm(), plume, #plume rise model stab, #stability class - eqs #equation set - ) + eqs, #equation set + dom) #problem domain end function (g::GaussianPlumeSolution{F,C,V,NoPlumeRise,S,E})(x, y, z, t=0) where {F,C,V,S,E} @@ -172,7 +179,7 @@ function (g::GaussianPlumeSolution{F,C,V,NoPlumeRise,S,E})(x, y, z, t=0) where { h = g.effective_stack_height if (x==0)&&(y==0)&&(z==h) return g.max_concentration - elseif (x≤0)||(z<0) + elseif _in_domain(x,y,z,g.domain) == false return zero(F) else G = g.rate @@ -197,7 +204,7 @@ function (g::GaussianPlumeSolution{F,C,V,<:BriggsModel,S,E})(x, y, z, t=0) where h = g.effective_stack_height if (x==0)&&(y==0)&&(z==h) return g.max_concentration - elseif (x≤0)||(z<0) + elseif _in_domain(x,y,z,g.domain) == false return zero(F) else G = g.rate diff --git a/src/utils/plume_rise.jl b/src/utils/plume_rise.jl index a17ba7f..654f44a 100644 --- a/src/utils/plume_rise.jl +++ b/src/utils/plume_rise.jl @@ -31,6 +31,16 @@ Base.isapprox(a::MomentumPlume, b::MomentumPlume) = all([ if typeof(getproperty(a,k))<:Number ]) +# some helper functions to limit the plume rise, e.g. when a mixing layer prevents the plume from rising too high +""" + _new_plume_rise(plume::BriggsModel, max_rise::Number) +Creates a new plume rise model with the same parameters as `plume`, but limits the final +rise to `max_rise`. This is useful when the plume rise exceeds the mixing height or +other constraints in the model. +""" +_new_plume_rise(plume::BuoyantPlume, max_rise::Number) = BuoyantPlume(plume.Fb, plume.xf, plume.u, max_rise) +_new_plume_rise(plume::MomentumPlume, max_rise::Number) = MomentumPlume(plume.Fm, plume.xf, plume.β, plume.u, plume.s, max_rise, plume.stab) + """ plume_rise(Dⱼ,uⱼ,Tᵣ,a::Atmosphere) Implements the Briggs plume rise equations for buoyancy and momentum driven diff --git a/src/utils/utils.jl b/src/utils/utils.jl index f07f7f8..e43fcf4 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -27,4 +27,20 @@ include("equation_sets/isc3.jl") # Default correlations DefaultSet = BasicEquationSet{DefaultWind,Nothing,Defaultσy,Defaultσz} -DefaultPuffSet = BasicEquationSet{DefaultWind,CCPSPuffσx,CCPSPuffσy,CCPSPuffσz} \ No newline at end of file +DefaultPuffSet = BasicEquationSet{DefaultWind,CCPSPuffσx,CCPSPuffσy,CCPSPuffσz} + +# a lazy check for if x,y,z are in the domain +struct ProblemDomain{F<:Number} + xmin::F + xmax::F + ymin::F + ymax::F + zmin::F + zmax::F +end + +function _in_domain(x::Number, y::Number, z::Number, domain::ProblemDomain) + return x ≥ domain.xmin && x ≤ domain.xmax && + y ≥ domain.ymin && y ≤ domain.ymax && + z ≥ domain.zmin && z ≤ domain.zmax +end \ No newline at end of file diff --git a/test/models/gaussian_plume_tests.jl b/test/models/gaussian_plume_tests.jl index e957a7a..9ac09a1 100644 --- a/test/models/gaussian_plume_tests.jl +++ b/test/models/gaussian_plume_tests.jl @@ -30,7 +30,7 @@ pl = plume(scn) # test default behaviour and type inheritance - @test GasDispersion.GaussianPlumeSolution(scn,:test,1.0,2,3,4,5,GasDispersion.SimpleCrossTerm(),GasDispersion.SimpleVerticalTerm(),GasDispersion.NoPlumeRise(),ClassA,DefaultSet()) isa GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, Nothing, GasDispersion.Defaultσy, GasDispersion.Defaultσz}} + @test GasDispersion.GaussianPlumeSolution(scn,:test,1.0,2,3,4,5,GasDispersion.SimpleCrossTerm(),GasDispersion.SimpleVerticalTerm(),GasDispersion.NoPlumeRise(),ClassA,DefaultSet(),GasDispersion.ProblemDomain(0.0,Inf,-Inf,Inf,0.0,Inf)) isa GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, Nothing, GasDispersion.Defaultσy, GasDispersion.Defaultσz}, GasDispersion.ProblemDomain{Float64}} @test pl isa GasDispersion.GaussianPlumeSolution @test pl isa Plume @test pl(-1,0,0) == 0.0 From b49e1665ccc336b067931c7a10509d2d48fe2132 Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Mon, 30 Jun 2025 17:27:01 -0600 Subject: [PATCH 7/8] Updated doc tests and unit tests Updated doc tests and unit tests to ensure better coverage, fix broken things. --- docs/src/equation_sets.md | 14 +++++++------- test/models/gaussian_ml_tests.jl | 14 ++++++++++++++ test/utils/util_tests.jl | 2 ++ 3 files changed, 23 insertions(+), 7 deletions(-) diff --git a/docs/src/equation_sets.md b/docs/src/equation_sets.md index eb6b0e4..ec202cb 100644 --- a/docs/src/equation_sets.md +++ b/docs/src/equation_sets.md @@ -117,7 +117,7 @@ dflt = plume(scn, GaussianPlume, DefaultSet()) # output -GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}()) +GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Defaultσy, Defaultσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf)) ``` @@ -128,7 +128,7 @@ ccps_rurl = plume(scn, GaussianPlume, CCPSRural()) # output -GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinRural, Nothing, BriggsRuralσy, BriggsRuralσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8420321686971456, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinRural, Nothing, BriggsRuralσy, BriggsRuralσz}()) +GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinRural, Nothing, BriggsRuralσy, BriggsRuralσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8420321686971456, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinRural, Nothing, BriggsRuralσy, BriggsRuralσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf)) ``` @@ -137,7 +137,7 @@ ccps_urb = plume(scn, GaussianPlume, CCPSUrban()) # output -GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinUrban, Nothing, BriggsUrbanσy, BriggsUrbanσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.7989729675905327, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinUrban, Nothing, BriggsUrbanσy, BriggsUrbanσz}()) +GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinUrban, Nothing, BriggsUrbanσy, BriggsUrbanσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.7989729675905327, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinUrban, Nothing, BriggsUrbanσy, BriggsUrbanσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf)) ``` @@ -146,7 +146,7 @@ isc3_rurl = plume(scn, GaussianPlume, ISC3Rural()) # output -GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinRural, Nothing, ISC3Ruralσy, ISC3Ruralσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8420321686971456, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinRural, Nothing, ISC3Ruralσy, ISC3Ruralσz}()) +GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{IrwinRural, Nothing, ISC3Ruralσy, ISC3Ruralσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8420321686971456, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{IrwinRural, Nothing, ISC3Ruralσy, ISC3Ruralσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf)) ``` @@ -155,7 +155,7 @@ isc3_urb = plume(scn, GaussianPlume, ISC3Urban()) # output -GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{ISC3UrbanWind, Nothing, BriggsUrbanσy, BriggsUrbanσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.0947417281650496, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{ISC3UrbanWind, Nothing, BriggsUrbanσy, BriggsUrbanσz}()) +GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{ISC3UrbanWind, Nothing, BriggsUrbanσy, BriggsUrbanσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.0947417281650496, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{ISC3UrbanWind, Nothing, BriggsUrbanσy, BriggsUrbanσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf)) ``` @@ -164,7 +164,7 @@ tno = plume(scn, GaussianPlume, TNOPlume()) # output -GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{TNOWind, Nothing, TNOPlumeσy, TNOPlumeσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8753751236458281, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{TNOWind, Nothing, TNOPlumeσy, TNOPlumeσz}()) +GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{TNOWind, Nothing, TNOPlumeσy, TNOPlumeσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 0.8753751236458281, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{TNOWind, Nothing, TNOPlumeσy, TNOPlumeσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf)) ``` @@ -173,7 +173,7 @@ turner = plume(scn, GaussianPlume, Turner()) # output -GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Turnerσy, Turnerσz}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Turnerσy, Turnerσz}()) +GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassF, BasicEquationSet{DefaultWind, Nothing, Turnerσy, Turnerσz}, GasDispersion.ProblemDomain{Float64}}(Scenario{Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}, HorizontalJet{Float64}, SimpleAtmosphere{Float64, ClassF}}(Substance{String, GasDispersion.Antoine{Float64}, Float64, Float64, Float64, Int64, Int64, Int64}("propane", 0.044096, GasDispersion.Antoine{Float64}(9.773719865868816, 2257.9247634130143, 0.0), 1.864931992847327, 526.13, 288.15, 101325.0, 1.142, 231.02, 425740, 1678, 2520), HorizontalJet{Float64}(0.08991798763471508, Inf, 0.01, 208.10961399327573, 3.5, 288765.2212333958, 278.3846872082166, 0.0), SimpleAtmosphere{Float64, ClassF}(101325.0, 298.15, 1.5, 10.0, 0.0, ClassF)), :gaussian, 0.08991798763471508, 0.9999999999999998, 1.8023818673116125, 1.150112899011524, 3.5, GasDispersion.SimpleCrossTerm(), GasDispersion.SimpleVerticalTerm(), GasDispersion.NoPlumeRise(), ClassF, BasicEquationSet{DefaultWind, Nothing, Turnerσy, Turnerσz}(), GasDispersion.ProblemDomain{Float64}(0.0, Inf, -Inf, Inf, 0.0, Inf)) ``` diff --git a/test/models/gaussian_ml_tests.jl b/test/models/gaussian_ml_tests.jl index f927513..f9819b0 100644 --- a/test/models/gaussian_ml_tests.jl +++ b/test/models/gaussian_ml_tests.jl @@ -48,6 +48,15 @@ temperature=298., pressure=101325., fraction_liquid=0) + + vj_rh_high = VerticalJet(mass_rate=0.1, + duration=Inf, + diameter=1, + velocity=100, + height=630, # plume rise will be too high + temperature=500., + pressure=101325., + fraction_liquid=0) stbl = SimpleAtmosphere(temperature=298, pressure=101325, @@ -73,8 +82,13 @@ @test plume(Scenario(sub,vj,stbl), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleVerticalTerm @test plume(Scenario(sub,vj,neut), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleMixingLayer @test plume(Scenario(sub,vj,neut), GaussianMixingLayer; method=:periodicmixinglayer).verticalterm isa GasDispersion.PeriodicMixingLayer + @test plume(Scenario(sub,vj,neut), GaussianMixingLayer; plumerise=false).plumerise isa GasDispersion.NoPlumeRise @test_throws ErrorException plume(Scenario(sub,vj,neut), GaussianMixingLayer; method=:someothermethod) @test_throws ErrorException plume(Scenario(sub,vj_toohigh,neut), GaussianMixingLayer) + + # testing limit on plume rise + pl_rise = plume(Scenario(sub,vj_rh_high,neut), GaussianMixingLayer; downwash=false, plumerise=true) + @test pl_rise.plumerise.final_rise ≈ pl_rise.verticalterm.mixing_height # integration test with a scenario -- simple mixing layer pl_smpl = plume(Scenario(sub,vj,neut), GaussianMixingLayer; downwash=true, plumerise=true) diff --git a/test/utils/util_tests.jl b/test/utils/util_tests.jl index 13766c3..3dd8f22 100644 --- a/test/utils/util_tests.jl +++ b/test/utils/util_tests.jl @@ -81,6 +81,8 @@ end # test null case @test isa(GasDispersion.NoPlumeRise(),GasDispersion.PlumeRise) + @test GasDispersion._new_plume_rise(GasDispersion.BuoyantPlume(50, 100, 10, 20), 200) == GasDispersion.BuoyantPlume(50, 100, 10, 200) + @test GasDispersion._new_plume_rise(GasDispersion.MomentumPlume(50, 100, 1.5, 10, 0.1, 20, ClassA), 200) == GasDispersion.MomentumPlume(50, 100, 1.5, 10, 0.1, 200, ClassA) # test Briggs model of plume rise g = 9.80616 # m/s^2 From 44f1588d400c754ad0d77c809c6e384da8ac3d9d Mon Sep 17 00:00:00 2001 From: allan <75338859+aefarrell@users.noreply.github.com> Date: Mon, 30 Jun 2025 17:39:14 -0600 Subject: [PATCH 8/8] Additional unit tests Added checks to ensure simple mixing layer exits the for loop correctly, either by breaking out early if the sum has converged or by running out the loop --- test/models/gaussian_ml_tests.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/models/gaussian_ml_tests.jl b/test/models/gaussian_ml_tests.jl index f9819b0..db9c10c 100644 --- a/test/models/gaussian_ml_tests.jl +++ b/test/models/gaussian_ml_tests.jl @@ -89,6 +89,8 @@ # testing limit on plume rise pl_rise = plume(Scenario(sub,vj_rh_high,neut), GaussianMixingLayer; downwash=false, plumerise=true) @test pl_rise.plumerise.final_rise ≈ pl_rise.verticalterm.mixing_height + pl_low_terms = plume(Scenario(sub,vj_rh_high,neut), GaussianMixingLayer; downwash=false, plumerise=true, n_terms=1) + @test pl_low_terms(pl_rise.plumerise.xf,0,0.9*pl_rise.verticalterm.mixing_height) < pl_rise(pl_rise.plumerise.xf,0,0.9*pl_rise.verticalterm.mixing_height) # integration test with a scenario -- simple mixing layer pl_smpl = plume(Scenario(sub,vj,neut), GaussianMixingLayer; downwash=true, plumerise=true)