Merge pull request #25 from aefarrell/britter-mcquaid-plume

Britter mcquaid plume

src/GasDispersion.jl

@@ -4,7 +4,7 @@ module GasDispersion
 
 # imports
 using Markdown
-using Interpolations: Extrapolation, Line, LinearInterpolation
+using Interpolations
 using SpecialFunctions: erf
 using RecipesBase
 

src/models/britter_mcquaid_plume.jl

@@ -1,111 +1,151 @@
+include("../utils/britter_mcquaid_correls.jl")
+
 struct BritterMcQuaidPlume <: PlumeModel end
 
 struct BritterMcQuaidPlumeSolution <: Plume
     scenario::Scenario
     model::Symbol
-    jet_density::Number
-    temperature_correction::Number
-    critical_distance::Number
-    interpolation::Extrapolation
+    c₀::Number # initial concentration, kg/m³
+    T′::Number # temperature correction
+    D::Number  # critical length
+    lb::Number # plume dimension parameter, m
+    itp::Interpolations.GriddedInterpolation
+    xnf::Number # near field distance
+    xff::Number # far field distance
+    A::Number  # far field constant
 end
 
-Britter_McQuaid_correlations = Dict(
-    0.010 => ( αs=[-1.0, -0.7, -0.29, -0.2, 1.0],
-               βs=[2.25, 2.25, 2.45, 2.45, 1.83]),
-    0.005 => ( αs=[-1.0, -0.67, -0.28, -0.15, 1.0],
-               βs=[2.4, 2.4, 2.63, 2.63, 2.07]),
-    0.020 => ( αs=[-1.0, -0.69, -0.31, -0.16, 1.0],
-               βs=[2.08, 2.08, 2.25, 2.25, 1.62]),
-    0.002 => ( αs=[-1.0, -0.69, -0.25, -0.13, 1.0],
-               βs=[2.6, 2.6, 2.77, 2.77, 2.21]),
-    0.100 => ( αs=[-1.0, -0.55, -0.14, 1.0],
-               βs=[1.75, 1.75, 1.85, 1.28]),
-    0.050 => ( αs=[-1.0, -0.68, -0.29, -0.18, 1.0],
-               βs=[1.92, 1.92, 2.06, 2.06, 1.4]),
-)
-
-
 """
-    plume(scenario::Scenario, BritterMcQuaidPlume)
+    plume(::Scenario, BritterMcQuaidPlume)
 
 Returns the solution to a Britter-McQuaid dispersion model for the given
 scenario.
 
-Currently only implements the max concentration at a downwind distance x, the
-other coordinates are ignored.
+# References
++ Britter, R.E. and J. McQuaid, *Workbook on the Dispersion of Dense Gases* HSE Contract Research Report No. 17/1988, 1988
 
 """
 function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume})
 
-    Q = _mass_rate(scenario)
-    h = _release_height(scenario)
+    Q = _release_flowrate(scenario)
+    ṁ = _mass_rate(scenario)
     ρⱼ = _release_density(scenario)
     Tᵣ = _release_temperature(scenario)
 
     u₁₀ = _windspeed(scenario, 10.0)
     ρₐ = _atmosphere_density(scenario)
     Tₐ = _atmosphere_temperature(scenario)
 
-    # Setting up the Britter-McQuaid curves
-    britter_interps = [ ]
-    concs = sort(collect(keys(Britter_McQuaid_correlations)), rev=true)
-
-    for conc in concs
-        αs, βs = Britter_McQuaid_correlations[conc]
-        f = LinearInterpolation(αs, βs, extrapolation_bc=Line())
-        push!(britter_interps, (c=conc, it=f))
-    end
+    # initial concentration
+    c₀ = ṁ/Q
 
     # relative density
     g = 9.80616  # gravity, m/s^2
     gₒ = g * ((ρⱼ - ρₐ)/ ρₐ)
 
     # critical distance
-    Vr = Q/ρⱼ # volumetric release rate
-    D = √(Vr/u₁₀)
+    D = √(Q/u₁₀)
 
     # temperature correction
-    T′ = Tₐ/Tᵣ
+    T′ = Tᵣ/Tₐ
 
     # correlation parameter
-    α = 0.2*log10( gₒ^2 * Vr * u₁₀^-5 )
+    α = 0.2*log10( gₒ^2 * Q * u₁₀^-5 )
 
     # setting up correlation for constant α
     # calculates the points for the linear interpolation
-    concs = [ elem.c for elem in britter_interps ]
-    βs = [ elem.it(α) for elem in britter_interps ]
+    concs, βs = _bm_pl_c(α)
+
+    # near-field location
+    xnf = 30
+    xmin = 10^(minimum(βs))
+    if xnf < xmin
+        # linear interpolation between the near-field correlation
+        # and the main correlation
+        βnf = log10(xnf)
+        cnf = 306*xnf^-2 / (1+ 306*xnf^-2)
+        βs = [ βnf; βs]
+        concs = [ cnf; concs ]
+    else
+        # for α > 0.605 the near-field overlaps with the correlation
+        # stop the near-field correlation early
+        # (this can lead to discontinuities)
+        xnf = xmin
+    end
+
+    # linear interpolation
+    itp = interpolate((βs,), concs, Gridded(Linear()))
 
-    # linear interpolation between the short distance correlation
-    # and the main correlation
-    βs = [ log10(30); βs]
-    concs = [ 306*30^-2 / (1+ 306*30^-2); concs ]
+    # far field correlation
+    # starts at last interpolation point and decays like x′^-2
+    xff = 10^(maximum(βs))
+    A = minimum(concs)*xff^2
 
-    # linear interpolation, extrapolates past the end with a straight line
-    interpolation = LinearInterpolation(βs, concs, extrapolation_bc=Line())
+    # plume dimension parameters
+    lb  = Q*gₒ/(u₁₀^3)
 
     return BritterMcQuaidPlumeSolution(
         scenario, #scenario::Scenario
         :brittermcquaid, #model::Symbol
-        ρⱼ,    #jet_density::Number
-        T′,    #temperature_correction::Number
-        D,     #D::Number
-        interpolation #interpolation::Extrapolation
+        c₀,    # initial concentration, kg/m³
+        T′,    # temperature_correction
+        D,     # critical length, m
+        lb,    # length parameter, m
+        itp,   # interpolation::Extrapolation
+        xnf,   # near-field distance
+        xff,   # far-field distance
+        A      # far-field constant
     )
 end
 
-function (b::BritterMcQuaidPlumeSolution)(x, y=0, z=0)
-    ρⱼ = b.jet_density
-    T′ = b.temperature_correction
-    D  = b.critical_distance
-    x′ = x/D
-    if x′ < 30
-        # use short distance correlation
-        c′ = x′ > 0 ? 306*x′^-2 / (1+ 306*x′^-2) : 1.0
+function (b::BritterMcQuaidPlumeSolution)(x, y, z)
+
+    # plume dimension parameters
+    Lu  = 0.5*b.D + 2*b.lb
+    Lh0 = b.D + 8*b.lb
+
+    # domain check
+    if (x < -Lu) || (z < 0)
+        return 0.0
+    end
+
+    # plume dimensions
+    if x < 0
+        # Britter-McQuaid model does not give a profile for x<0
+        # assuming crosswind length is Lh0
+        Lh = Lh0
     else
+        Lh = Lh0 + 2.5*cbrt(b.lb*x^2)
+    end
+
+    # within the crosswind extent of the plume
+    if (y < -Lh) || (y > Lh)
+        return 0.0
+    end
+
+    x′ = x/b.D
+    if x′ ≤ b.xnf
+        # use near-field correlation
+        c′ = x′ > 0 ? 306*x′^-2 / (1+ 306*x′^-2) : 1.0
+    elseif b.xnf < x′ < b.xff
         # use linear interpolation
         β = log10(x′)
-        c′ = b.interpolation(β)
+        c′ = b.itp(β)
+    else
+        # use far-field correlation
+        # where A is a function of α
+        c′ = b.A/(x′^2)
+    end
+
+    # non-isothermal correction
+    c′ = c′ / (c′ + (1 - c′)*b.T′)
+    c = b.c₀*c′
+
+    # vertical extent, from continuity
+    Lv = (b.D^2)/(2*Lh*c′)
+    if z ≤ Lv
+        return c
+    else
+        return 0.0
     end
-    c = ( ρⱼ*c′*T′) / (1 - c′ + c′*T′)
-    return c
 end

src/models/gaussian_plume.jl

@@ -18,9 +18,6 @@ end
 
 Returns the solution to a Gaussian plume dispersion model for the given scenario.
 
-The Gaussian plume model is per *Guidelines for Consequence Analysis of Chemical
-Release*, CCPS, New York (1999)
-
 ```math
 c\left(x,y,z\right) = { \dot{m} \over { 2 \pi \sigma_{y} \sigma_{z} u } }
 \exp \left[ -\frac{1}{2} \left( y \over \sigma_{y} \right)^2 \right] \\
@@ -30,9 +27,13 @@ c\left(x,y,z\right) = { \dot{m} \over { 2 \pi \sigma_{y} \sigma_{z} u } }
 
 where the σs are dispersion parameters correlated with the distance x
 
+# References
++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999)
+
 # Arguments
 - `downwash::Bool=false`: when true, includes stack-downwash effects
 - `plumerise::Bool=false`: when true, includes plume-rise effects using Briggs' model
+
 """
 function plume(scenario::Scenario, ::Type{GaussianPlume}; downwash::Bool=false, plumerise::Bool=false)
     # parameters of the jet

src/models/gaussian_puff.jl

@@ -12,13 +12,10 @@ struct GaussianPuffSolution{S<:StabilityClass} <: Puff
 end
 
 @doc doc"""
-    puff(scenario::Scenario, GaussianPuff)
+    puff(::Scenario, GaussianPuff)
 
 Returns the solution to a Gaussian puff dispersion model for the given scenario.
 
-The Gaussian puff model is per *Guidelines for Consequence Analysis of Chemical
-Release*, CCPS, New York (1999)
-
 ```math
 c\left(x,y,z,t\right) = \dot{m} \Delta t
 { { \exp \left( -\frac{1}{2} \left( {x - u t } \over \sigma_x \right)^2 \right) } \over { \sqrt{2\pi} \sigma_x } }
@@ -27,6 +24,9 @@ c\left(x,y,z,t\right) = \dot{m} \Delta t
 + \exp \left( -\frac{1}{2} \left( {z + h} \over \sigma_z \right)^2 \right) } \over { \sqrt{2\pi} \sigma_z } }
 ```
 
+# References
++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999)
+
 """
 function puff(scenario::Scenario, ::Type{GaussianPuff})
 

src/models/intpuff.jl

@@ -13,7 +13,7 @@ struct IntPuffSolution{T<:Number,S<:StabilityClass} <: Puff
 end
 
 @doc doc"""
-    puff(scenario::Scenario, IntPuff; kwargs...)
+    puff(::Scenario, IntPuff; kwargs...)
 
 Returns the solution to an integrated Gaussian dispersion model, where the
 release is modeled as a sequence of Gaussian puffs, for the given scenario.

src/models/simple_jet.jl

@@ -14,14 +14,11 @@ struct SimpleJetSolution <: Plume
 end
 
 @doc doc"""
-    plume(scenario::Scenario, SimpleJet; kwargs...)
+    plume(::Scenario, SimpleJet; kwargs...)
 
 Returns the solution to a simple turbulent jet dispersion model for the given
 scenario.
 
-The turbulent jet model is per Long, V.D., "Estimation of the Extent of Hazard
-Areas Around a Vent", *Chem. Process Hazard*, II, 6, 1963
-
 ```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)
@@ -31,10 +28,14 @@ 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., "Estimation of the Extent of Hazard Areas Around a Vent" *Chem. Process Hazard*, II:6 (1963)
+
 # Arguments
 - `release_angle::Number=0`: 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, ::Type{SimpleJet}; release_angle::Number=0.0, k2::Number=6.0, k3::Number=5.0)
     # Density correction

src/source_models/jet_source.jl

@@ -7,8 +7,8 @@ Returns returns a scenario for a simple jet from a circular hole. The
 jet can either be a liquid or a gas (in which case it is assumed to be an ideal
 gas and the jet is isentropic).
 
-Liquid and gas discharge models are per *Guidelines for Consequence Analysis of
-Chemical Release*, CCPS, New York (1999)
+# References
++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999)
 
 # Arguments
 - `phase=:liquid`: the phase, either :liquid or :gas