Skip to content

Commit

Permalink
Testing correlations
Browse files Browse the repository at this point in the history 
Moved Britter-McQuaid correlations to utils, added tests for correlations.

Changed the definition for the vertical extent, I'm fairly sure there is a mistake in the Britter-McQuaid workbook as the vertical extent formula makes no sense -- the plume shrinks rapidly until it is millimeters tall causing mass to disappear out of the system.
  • Loading branch information
@aefarrell
aefarrell committed Jan 28, 2023
1 parent a641d71 commit f3d905b
Show file tree
Hide file tree
Showing 3 changed files with 144 additions and 117 deletions.
130 changes: 15 additions & 115 deletions src/models/britter_mcquaid_plume.jl
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
include("../utils/britter_mcquaid_correls.jl")

struct BritterMcQuaidPlume <: PlumeModel end

struct BritterMcQuaidPlumeSolution <: Plume
Expand All @@ -13,115 +15,6 @@ struct BritterMcQuaidPlumeSolution <: Plume
A::Number # far field constant
end

# Britter-McQuaid plume correlations
# as digitzed in TSCREEN
function _bm_pl_c_1(α)
# corresponds to (Cm/C0) = 0.10
if α ≤ -0.55
return 1.75
elseif -0.55 < α ≤ -0.14
return 0.24*α + 1.88
elseif -0.14 < α ≤ 1.0
return -0.5*α + 1.78
else
@warn "α= $α is out of range"
return -0.5*α + 1.78
end
end

function _bm_pl_c_05(α)
# corresponds to (Cm/C0) = 0.05
if α ≤ -0.68
return 1.92
elseif -0.68 < α ≤ -0.29
return 0.36*α + 2.16
elseif -0.29 < α ≤ -0.18
return 2.06
elseif -0.18 < α ≤ 1.0
return -0.56*α + 1.96
else
@warn "α= $α is out of range"
return -0.56*α + 1.96
end
end

function _bm_pl_c_02(α)
# corresponds to (Cm/C0) = 0.02
if α ≤ -0.69
return 2.08
elseif -0.69 < α ≤ -0.31
return 0.45*α + 2.39
elseif -0.31 < α ≤ -0.16
return 2.25
elseif -0.16 < α ≤ 1.0
return -0.54*α + 2.16
else
@warn "α= $α is out of range"
return -0.54*α + 2.16
end
end

function _bm_pl_c_01(α)
# corresponds to (Cm/C0) = 0.01
if α ≤ -0.70
return 2.25
elseif -0.70 < α ≤ -0.29
return 0.49*α + 2.59
elseif -0.29 < α ≤ -0.20
return 2.45
elseif -0.20 < α ≤ 1.0
return -0.52*α + 2.35
else
@warn "α= $α is out of range"
return -0.52*α + 2.35
end
end

function _bm_pl_c_005(α)
# corresponds to (Cm/C0) = 0.005
if α ≤ -0.67
return 2.40
elseif -0.67 < α ≤ -0.28
return 0.59*α + 2.80
elseif -0.28 < α ≤ -0.15
return 2.63
elseif -0.15 < α ≤ 1.0
return -0.48*α + 2.56
else
@warn "α= $α is out of range"
return -0.48*α + 2.56
end
end

function _bm_pl_c_002(α)
# corresponds to (Cm/C0) = 0.002
if α ≤ -0.69
return 2.60
elseif -0.69 < α ≤ -0.25
return 0.39*α + 2.87
elseif -0.25 < α ≤ -0.13
return 2.77
elseif -0.13 < α ≤ 1.0
return -0.50*α + 2.71
else
@warn "α= $α is out of range"
return -0.50*α + 2.71
end
end

function _bm_pl_c(α)
concs = [0.10, 0.05, 0.02, 0.01, 0.005, 0.002]
βs = [
_bm_pl_c_1(α), # (C_m/C_0) = 0.1
_bm_pl_c_05(α), # (C_m/C_0) = 0.05
_bm_pl_c_02(α), # (C_m/C_0) = 0.02
_bm_pl_c_01(α), # (C_m/C_0) = 0.01
_bm_pl_c_005(α), # (C_m/C_0) = 0.005
_bm_pl_c_002(α) # (C_m/C_0) = 0.002
]
return concs, βs
end

"""
plume(scenario::Scenario, BritterMcQuaidPlume)
Expand Down Expand Up @@ -212,7 +105,7 @@ function (b::BritterMcQuaidPlumeSolution)(x, y, z)
Lh0 = b.D + 8*b.lb

# domain check
if x < -Lu
if (x < -Lu) || (z < 0)
return 0.0
end

Expand All @@ -225,9 +118,8 @@ function (b::BritterMcQuaidPlumeSolution)(x, y, z)
Lh = Lh0 + 2.5*cbrt(b.lb*x^2)
end

Lv = (b.D^2)/Lh
# within the extent of the plume
if (y < -Lh) || (y > Lh) || (z < 0) || (z > Lv)
# within the crosswind extent of the plume
if (y < -Lh) || (y > Lh)
return 0.0
end

Expand All @@ -246,6 +138,14 @@ function (b::BritterMcQuaidPlumeSolution)(x, y, z)
end

# non-isothermal correction
c = b.c₀*c′ / (c′ + (1 - c′)*b.T′)
return c
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
end
108 changes: 108 additions & 0 deletions src/utils/britter_mcquaid_correls.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,108 @@
# Britter-McQuaid plume correlations
# as digitzed in TSCREEN
function _bm_pl_c_1(α)
# corresponds to (Cm/C0) = 0.10
if α ≤ -0.55
return 1.75
elseif -0.55 < α ≤ -0.14
return 0.24*α + 1.88
elseif -0.14 < α ≤ 1.0
return -0.5*α + 1.78
else
@warn "α= $α is out of range"
return -0.5*α + 1.78
end
end

function _bm_pl_c_05(α)
# corresponds to (Cm/C0) = 0.05
if α ≤ -0.68
return 1.92
elseif -0.68 < α ≤ -0.29
return 0.36*α + 2.16
elseif -0.29 < α ≤ -0.18
return 2.06
elseif -0.18 < α ≤ 1.0
return -0.56*α + 1.96
else
@warn "α= $α is out of range"
return -0.56*α + 1.96
end
end

function _bm_pl_c_02(α)
# corresponds to (Cm/C0) = 0.02
if α ≤ -0.69
return 2.08
elseif -0.69 < α ≤ -0.31
return 0.45*α + 2.39
elseif -0.31 < α ≤ -0.16
return 2.25
elseif -0.16 < α ≤ 1.0
return -0.54*α + 2.16
else
@warn "α= $α is out of range"
return -0.54*α + 2.16
end
end

function _bm_pl_c_01(α)
# corresponds to (Cm/C0) = 0.01
if α ≤ -0.70
return 2.25
elseif -0.70 < α ≤ -0.29
return 0.49*α + 2.59
elseif -0.29 < α ≤ -0.20
return 2.45
elseif -0.20 < α ≤ 1.0
return -0.52*α + 2.35
else
@warn "α= $α is out of range"
return -0.52*α + 2.35
end
end

function _bm_pl_c_005(α)
# corresponds to (Cm/C0) = 0.005
if α ≤ -0.67
return 2.40
elseif -0.67 < α ≤ -0.28
return 0.59*α + 2.80
elseif -0.28 < α ≤ -0.15
return 2.63
elseif -0.15 < α ≤ 1.0
return -0.48*α + 2.56
else
@warn "α= $α is out of range"
return -0.48*α + 2.56
end
end

function _bm_pl_c_002(α)
# corresponds to (Cm/C0) = 0.002
if α ≤ -0.69
return 2.60
elseif -0.69 < α ≤ -0.25
return 0.39*α + 2.87
elseif -0.25 < α ≤ -0.13
return 2.77
elseif -0.13 < α ≤ 1.0
return -0.50*α + 2.71
else
@warn "α= $α is out of range"
return -0.50*α + 2.71
end
end

function _bm_pl_c(α)
concs = [0.10, 0.05, 0.02, 0.01, 0.005, 0.002]
βs = [
_bm_pl_c_1(α), # (C_m/C_0) = 0.1
_bm_pl_c_05(α), # (C_m/C_0) = 0.05
_bm_pl_c_02(α), # (C_m/C_0) = 0.02
_bm_pl_c_01(α), # (C_m/C_0) = 0.01
_bm_pl_c_005(α), # (C_m/C_0) = 0.005
_bm_pl_c_002(α) # (C_m/C_0) = 0.002
]
return concs, βs
end
23 changes: 21 additions & 2 deletions test/models/britter_mcquaid_plume_tests.jl
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,22 @@
@testset "Britter-McQuaid plume correlations" begin
include("../../src/utils/britter_mcquaid_correls.jl")

@test _bm_pl_c(-0.75)[2] ≈ [1.75, 1.92, 2.08, 2.25, 2.4, 2.6]
@test _bm_pl_c(-0.40)[2] ≈ [1.7839999999999998, 2.016, 2.21, 2.3939999999999997, 2.564, 2.714]
@test _bm_pl_c(-0.20)[2] ≈ [1.8319999999999999, 2.06, 2.25, 2.45, 2.63, 2.77]
@test _bm_pl_c(0.0)[2] ≈ [1.78, 1.96, 2.16, 2.35, 2.56, 2.71]

# test domain warning
@test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_1(1.1)
@test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_05(1.1)
@test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_02(1.1)
@test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_01(1.1)
@test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_005(1.1)
@test_logs (:warn, "α= 1.1 is out of range") _bm_pl_c_002(1.1)


end

@testset "Britter-McQuaid plume tests" begin
# Britter-McQuaid example, *Guidelines for Consequence Analysis of Chemical
# Releases* CCPS, 1999, pg 122
Expand Down Expand Up @@ -30,7 +49,7 @@
# far field
x₃, c₃ = 1200.0, 0.008660197082795383

# test type inheritance
# test overall solution
pl = plume(scn, BritterMcQuaidPlume)
@test isa(pl, GasDispersion.BritterMcQuaidPlumeSolution)
@test isa(pl, Plume)
Expand All @@ -49,7 +68,7 @@
@test pl(0 - 2*eps(Float64), Lh0 + 2*eps(Float64), 0) == 0.0
@test pl(0 - 2*eps(Float64), Lh0 - 2*eps(Float64), 0) ≈ c₀

Lv = pl.D^2/Lh0
Lv = pl.D^2/(2*Lh0)
@test pl(0, 0, -1) == 0.0
@test pl(0, 0, Lv + 2*eps(Float64)) == 0.0
@test pl(0, 0, Lv - 2*eps(Float64)) ≈ c₀
Expand Down

0 comments on commit f3d905b

Please sign in to comment.