Merge pull request #26 from aefarrell/britter-mcquaid-puff

Britter mcquaid puff

docs/src/index.md

@@ -167,3 +167,14 @@ default behaviour is to take the limit $n \to \infty$.
 ```@docs
 puff(::Scenario, ::Type{IntPuff})
 ```
+
+### Britter-McQuaid Model
+
+The Britter-McQuaid model is an empirical correlation for dense cloud
+dispersion. The model generates an interpolation function for the average cloud
+concentration and the cloud is rendered as a cylinder.
+
+
+```@docs
+puff(::Scenario, ::Type{BritterMcQuaidPuff})
+```

src/models/britter_mcquaid_plume.jl

@@ -1,5 +1,3 @@
-include("../utils/britter_mcquaid_correls.jl")
-
 struct BritterMcQuaidPlume <: PlumeModel end
 
 struct BritterMcQuaidPlumeSolution <: Plume
@@ -18,11 +16,12 @@ end
 """
     plume(::Scenario, BritterMcQuaidPlume)
 
-Returns the solution to a Britter-McQuaid dispersion model for the given
-scenario.
+Returns the solution to the Britter-McQuaid continuous ground level release
+model for the given scenario.
 
 # References
 + Britter, R.E. and J. McQuaid, *Workbook on the Dispersion of Dense Gases* HSE Contract Research Report No. 17/1988, 1988
++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999)
 
 """
 function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume})
@@ -63,7 +62,7 @@ function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume})
         # linear interpolation between the near-field correlation
         # and the main correlation
         βnf = log10(xnf)
-        cnf = 306*xnf^-2 / (1+ 306*xnf^-2)
+        cnf = 306 / (xnf^2 + 306)
         βs = [ βnf; βs]
         concs = [ cnf; concs ]
     else
@@ -91,7 +90,7 @@ function plume(scenario::Scenario, ::Type{BritterMcQuaidPlume})
         T′,    # temperature_correction
         D,     # critical length, m
         lb,    # length parameter, m
-        itp,   # interpolation::Extrapolation
+        itp,   # interpolation
         xnf,   # near-field distance
         xff,   # far-field distance
         A      # far-field constant
@@ -126,7 +125,7 @@ function (b::BritterMcQuaidPlumeSolution)(x, y, z)
     x′ = x/b.D
     if x′ ≤ b.xnf
         # use near-field correlation
-        c′ = x′ > 0 ? 306*x′^-2 / (1+ 306*x′^-2) : 1.0
+        c′ = x′ > 0 ? 306 / ( 306 + x′^2 ) : 1.0
     elseif b.xnf < x′ < b.xff
         # use linear interpolation
         β = log10(x′)

src/models/britter_mcquaid_puff.jl

@@ -1,14 +1,127 @@
 struct BritterMcQuaidPuff <: PuffModel end
 
+struct BritterMcQuaidPuffSolution <: Puff
+    scenario::Scenario
+    model::Symbol
+    c₀::Number    # initial concentration, kg/m³
+    T′::Number    # temperature_correction
+    V₀::Number    # initial volume, m³
+    gₒ::Number    # reduced gravity, m/s²
+    u₁₀::Number   # referebce windspeed, m/s
+    itp::Interpolations.GriddedInterpolation
+    xnf::Number   # near-field distance
+    xff::Number   # far-field distance
+    A::Number     # far-field constant
+end
+
 """
     puff(scenario::Scenario, BritterMcQuaidPuff)
 
-Generates a Britter-McQuaid dispersion model on the given scenario and returns a
-callable struct giving the concentration of the form
-c(x, y, z, t)
+Returns the solution to the Britter-McQuaid instantaneous ground level release
+model for the given scenario.
+
+# References
++ Britter, R.E. and J. McQuaid, *Workbook on the Dispersion of Dense Gases* HSE Contract Research Report No. 17/1988, 1988
++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999)
 
 """
 function puff(scenario::Scenario, ::Type{BritterMcQuaidPuff})
-    # TODO
-    error("Britter-McQuaid puff dispersion model is not currently implemented")
+
+    Q = _release_flowrate(scenario)
+    ṁ = _mass_rate(scenario)
+    t = _duration(scenario)
+    ρⱼ = _release_density(scenario)
+    Tᵣ = _release_temperature(scenario)
+
+    u₁₀ = _windspeed(scenario, 10.0)
+    ρₐ = _atmosphere_density(scenario)
+    Tₐ = _atmosphere_temperature(scenario)
+
+    # initial cloud dimensions
+    V₀ = Q*t
+
+    # initial concentration
+    c₀ = ṁ/Q
+
+    # temperature correction
+    T′ = Tᵣ/Tₐ
+
+    # relative density
+    g = 9.80616  # gravity, m/s^2
+    gₒ = g * ((ρⱼ - ρₐ)/ ρₐ)
+
+    # correlation parameter
+    α = 0.5*log10( gₒ * cbrt(V₀) / u₁₀^2 )
+
+    # setting up correlation for constant α
+    # calculates the points for the linear interpolation
+    concs, βs = _bm_pf_c(α)
+
+    # near-field location
+    xnf = 10^(minimum(βs))
+
+    # linear interpolation
+    itp = interpolate((βs,), concs, Gridded(Linear()))
+
+    # far field correlation
+    # starts at last interpolation point and decays like x′^-2
+    xff = 10^(maximum(βs))
+    A = minimum(concs)*xff^2
+
+    return BritterMcQuaidPuffSolution(
+        scenario, #scenario::Scenario
+        :brittermcquaid, #model::Symbol
+        c₀,    # initial concentration, kg/m³
+        T′,    # temperature_correction
+        V₀,    # initial volume, m³
+        gₒ,    # reduced gravity, m/s²
+        u₁₀,   # referebce windspeed, m/s
+        itp,   # interpolation::Extrapolation
+        xnf,   # near-field distance
+        xff,   # far-field distance
+        A      # far-field constant
+    )
+
+end
+
+function (pf::BritterMcQuaidPuffSolution)(x,y,z,t)
+
+    R₀ = cbrt(3*pf.V₀/4π)
+    xc = 0.4*pf.u₁₀*t
+    R² = R₀^2 + 1.2*√(pf.gₒ*pf.V₀)*t
+    r² = (x-xc)^2 + y^2
+
+    #domain check
+    if r² > R² || z < 0
+        return 0.0
+    end
+
+    x′ = x/cbrt(pf.V₀)
+    if x′ ≤ pf.xnf
+        # use near-field correlation
+        c′ = x′ > 0 ? 3.24 / (3.24 + x′^2) : 1.0
+
+        # don't drop below the first correlation concentration
+        c′ = max(c′, maximum(pf.itp.coefs))
+    elseif pf.xnf < x′ < pf.xff
+        # use linear interpolation
+        β = log10(x′)
+        c′ = pf.itp(β)
+    else
+        # use far-field correlation
+        # where A is a function of α
+        c′ = pf.A/(x′^2)
+    end
+
+    # non-isothermal correction
+    c′ = c′ / (c′ + (1 - c′)*pf.T′)
+    c = pf.c₀*c′
+
+    # vertical extent, from continuity
+    H = pf.V₀/(c′*π*R²)
+    if z ≤ H
+        return c
+    else
+        return 0.0
+    end
 end

src/utils/britter_mcquaid_correls.jl

@@ -1,4 +1,4 @@
-
+# plume correlations
 function _bm_pl_c_1(α)
     # corresponds to (Cm/C0) = 0.10
     if α ≤ -0.55
@@ -96,10 +96,10 @@ end
 """
     _bm_pl_c(α)
 
-Britter-McQuaid plume correlations as digtized in TSCREEN
+Britter-McQuaid plume correlations
 
 # References
-+ EPA, *User's Guide to TSCREEN* U.S. Environmental Protection Agency EPA-454/B-94-023 (1994)
++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999)
 
 """
 function _bm_pl_c(α)
@@ -114,3 +114,130 @@ function _bm_pl_c(α)
     ]
     return concs, βs
 end
+
+
+###############################################################################
+# puff correlations
+
+function _bm_pf_c_1(α)
+    # corresponds to (Cm/C0) = 0.10
+    if α ≤ -0.44
+        return 0.70
+    elseif -0.44 < α ≤ 0.43
+        return 0.26*α + 0.81
+    elseif 0.43 < α ≤ 1.0
+        return 0.93
+    else
+        @warn "α= $α is out of range"
+        return 0.93
+    end
+end
+
+function _bm_pf_c_05(α)
+    # corresponds to (Cm/C0) = 0.05
+    if α ≤ -0.56
+        return 0.85
+    elseif -0.56 < α ≤ 0.31
+        return 0.26*α + 1.00
+    elseif 0.31 < α ≤ 1.0
+        return -0.12*α + 1.12
+    else
+        @warn "α= $α is out of range"
+        α = min(α,10.0)
+        return -0.12*α + 1.12
+    end
+end
+
+function _bm_pf_c_02(α)
+    # corresponds to (Cm/C0) = 0.02
+    if α ≤ -0.66
+        return 0.95
+    elseif -0.66 < α ≤ 0.32
+        return 0.36*α + 1.19
+    elseif 0.32 < α ≤ 1.0
+        return -0.26*α + 1.38
+    else
+        @warn "α= $α is out of range"
+        α = min(α,10.0)
+        return -0.26*α + 1.38
+    end
+end
+
+function _bm_pf_c_01(α)
+    # corresponds to (Cm/C0) = 0.01
+    if α ≤ -0.71
+        return 1.15
+    elseif -0.71 < α ≤ 0.37
+        return 0.34*α + 1.39
+    elseif 0.37 < α ≤ 1.0
+        return -0.38*α + 1.66
+    else
+        @warn "α= $α is out of range"
+        α = min(α,10.0)
+        return -0.38*α + 1.66
+    end
+end
+
+function _bm_pf_c_005(α)
+    # corresponds to (Cm/C0) = 0.005
+    if α ≤ -0.52
+        return 1.48
+    elseif -0.52 < α ≤ 0.24
+        return 0.26*α + 1.62
+    elseif 0.24 < α ≤ 1.0
+        return -0.30*α + 1.75
+    else
+        @warn "α= $α is out of range"
+        α = min(α,10.0)
+        return -0.30*α + 1.75
+    end
+end
+
+function _bm_pf_c_002(α)
+    # corresponds to (Cm/C0) = 0.002
+    if α ≤ 0.27
+        return 1.83
+    elseif 0.27 < α ≤ 1.0
+        return -0.32*α + 1.92
+    else
+        @warn "α= $α is out of range"
+        α = min(α,10.0)
+        return -0.32*α + 1.92
+    end
+end
+
+function _bm_pf_c_001(α)
+    # corresponds to (Cm/C0) = 0.001
+    if α ≤ -0.10
+        return 2.075
+    elseif -0.10 < α ≤ 1.0
+        return -0.27*α + 2.05
+    else
+        @warn "α= $α is out of range"
+        α = min(α,10.0)
+        return -0.27*α + 2.05
+    end
+end
+
+"""
+    _bm_pf_c(α)
+
+Britter-McQuaid puff correlations
+
+# References
++ CCPS, *Guidelines for Consequence Analysis of Chemical Releases*, American Institute of Chemical Engineers, New York (1999)
+
+"""
+function _bm_pf_c(α)
+    concs = [0.10, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001]
+    βs = [
+            _bm_pf_c_1(α),   # (C_m/C_0) = 0.1
+            _bm_pf_c_05(α),  # (C_m/C_0) = 0.05
+            _bm_pf_c_02(α),  # (C_m/C_0) = 0.02
+            _bm_pf_c_01(α),  # (C_m/C_0) = 0.01
+            _bm_pf_c_005(α), # (C_m/C_0) = 0.005
+            _bm_pf_c_002(α), # (C_m/C_0) = 0.002
+            _bm_pf_c_001(α)  # (C_m/C_0) = 0.001
+    ]
+    return concs, βs
+end

src/utils/utils.jl

@@ -79,3 +79,6 @@ include("pasquill_gifford.jl")
 
 # Briggs' model for plume rise
 include("plume_rise.jl")
+
+# model correlations
+include("britter_mcquaid_correls.jl")