Palazzi short duration releases

Project.toml

@@ -1,7 +1,7 @@
 name = "GasDispersion"
 uuid = "ac6af613-5479-4178-853e-9d66de417193"
 authors = ["Allan Farrell <allanf@protonmail.com> and contributors"]
-version = "0.3.0"
+version = "0.4.0"
 
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"

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})

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)

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")

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

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

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

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

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")