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