Added periodic mixing layer

Added mixing layer model based on an infinite sum of cosines. This converges very slowly (in my testing) but is used in some places.

Author: aefarrell

docs/src/plume.md

@@ -222,23 +222,25 @@ The gaussian mixing layer model allows for the plume to be contained entirely wi
 c(x, y, z) = \frac{m_i}{u} \frac{1}{\sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z
 ```
 
-Where $F_z$ is the vertical dispersion term. For the method of images this takes the form of the infinite sum:
+Where $F_z$ is the vertical dispersion term. For the method of images this takes the form of the infinite sum ([Beychok 1994](references.md) p 123)
 
 ```math
 F_z = \frac{1}{\sqrt{2\pi} \sigma_z} \left( \exp \left( -\frac{1}{2} \left( { z -h } \over \sigma_{z} \right)^2 \right)
 + \exp \left( -\frac{1}{2} \left( { z + h } \over \sigma_{z} \right)^2 \right) \right) \\
-\sum_{i=1}^{n} \left[ \exp\left(-\frac{1}{2}\left(\frac{z - h + 2i h_m}{\sigma_z}\right)^2\right) 
-+ \exp\left(-\frac{1}{2}\left(\frac{z - h - 2i h_m}{\sigma_z}\right)^2\right) \\
-+ \exp\left(-\frac{1}{2}\left(\frac{z + h + 2i h_m}{\sigma_z}\right)^2\right) \\
-+ \exp\left(-\frac{1}{2}\left(\frac{z + h - 2i h_m}{\sigma_z}\right)^2\right) \right]
+\sum_{n=1}^{\infty} \left[ \exp\left(-\frac{1}{2}\left(\frac{z - h + 2n h_m}{\sigma_z}\right)^2\right) 
++ \exp\left(-\frac{1}{2}\left(\frac{z - h - 2n h_m}{\sigma_z}\right)^2\right) \\
++ \exp\left(-\frac{1}{2}\left(\frac{z + h + 2n h_m}{\sigma_z}\right)^2\right) \\
++ \exp\left(-\frac{1}{2}\left(\frac{z + h - 2n h_m}{\sigma_z}\right)^2\right) \right]
 ```
 
-For the periodic boundary it takes the form of another infinite sum:
+For the periodic boundary it takes the form of another infinite sum ([Seinfeld and Pandis 2006](references.md) p 858)
 
 ```math
-F_z = ...
+F_z = \frac{2}{h_m} \left( \frac{1}{2} + \sum_{n=1}^{\infty} cos\left( {n\pi z} \over h_m \right) cos\left( {n\pi h} \over h_m \right) \exp\left(-\frac{1}{2}\left(\frac{n\pi \sigma_z}{h_m}\right)^2\right) \right)
 ```
 
+The periodic boundary layer approach can be very slow to converge, in which case the number of terms given by the keyword argument `n_terms` should be increased. By default the sums terminate early if the solution converges, e.g. if `n_terms=100_000_000` or some other huge number but the sum converges after 5 terms, only 5 terms will be calculated. The simple mixing layer approach converges much more quickly, and typically 10 terms are more than enough.
+
 As `SimpleAtmosphere`s do not define mixing height, one is calculated based on the following:
 - for stable atmospheres (class E and F) the mixing height is assumed to be infinite, and the model defaults back to a [Simple Gaussian Plume](@ref)
 - for unstable and neutral atmospheres (classes A through D) the mixing height is calculated from the friction velocity $u^{*}$ and the coriolis parameter $f$: $h_m = 0.3 \frac{u^{*}}{f}$ where $f$ is calculated at a default latitude of 40°N (consistent with [US EPA 1995](references.md)).

docs/src/references.md

@@ -2,6 +2,8 @@
 
 AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers
 
+Beychok, Milton R. 1994. *Fundamentals of Stack Gas Dispersion* 3rd Ed. Irvine, CA: Milton R. Beychok.
+
 Bakkum, E.A. and N.J. Duijm. 2005. "Chapter 4 - Vapour Cloud Dispersion" in *Methods for the Calculation of Physical Effects, CPR 14E* (TNO Yellow Book) Edited by C.J.H. van den Bosch and R.A.P.M. Weterings. The Netherlands.
 
 Businger, J. A., J. C. Wyngaard, Y. Izumi, and E. F. Bradley. 1971. "Flux-Profile Relationships in the Atmospheric Surfaace Layer." *Journal of the Atmospheric Sciences*. 28, 181-189. [doi:10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2](https://doi.org/10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2)
@@ -30,7 +32,7 @@ Paulson, C. A. 1970. "The Mathematical Representation of Wind Speed and Temperat
 
 Spicer, Thomas O. and Jerry A. Havens. 1988. *Development of Vapor Dispersion Models for Non-Neutrally Buoyant Gas Mixtures--Analysis of TFI/NH3 Test Data*. Tyndall Air Force Base, FL: USAF Engineering & Services Laboratory
 
-Seinfeld, John H. 1986. *Atmospheric Chemistry and Physics of Air Pollution*. New York: John Wiley and Sons
+Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons.
 
 Turner, D. Bruce. 1970. *Workbook of Atmospheric Dispersion Estimates*. United States Environmental Protection Agency.
 

src/models/gaussian_mixing_layer.jl

@@ -1,7 +1,7 @@
 # for dispatch
 struct GaussianMixingLayer <: PlumeModel end
 
-# new vertical term
+# mixing layer vertical term -- method of images
 struct SimpleMixingLayer{N<:Integer,F<:Number} <: GaussianVerticalTerm
     nterms::N
     mixing_height::F
@@ -23,27 +23,50 @@ function vertical_term(z, h, σz, ml::SimpleMixingLayer)
     return Fz/(√(2π)*σz)
 end
 
+# mixing layer vertical term -- periodic
+struct PeriodicMixingLayer{N<:Integer,F<:Number} <: GaussianVerticalTerm
+    nterms::N
+    mixing_height::F
+end
+
+function vertical_term(z, h, σz, ml::PeriodicMixingLayer) 
+    Fz = 1/2
+    for n=1:ml.nterms
+        next_term = cos(n*π*z/ml.mixing_height)*cos(n*π*h/ml.mixing_height)*exp(-0.5*(n*π*σz/ml.mixing_height)^2)
+        Fz += next_term
+        if next_term ≈ 0
+            break
+        end
+    end
+    return 2*Fz/ml.mixing_height
+end
+
 @doc doc"""
     plume(::Scenario, GaussianMixingLayer[, ::EquationSet]; **kwargs...)
 
 Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer.
 
 ```math
-c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z
+c(x, y, z) = \frac{m_i}{u} \frac{1}{ \sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z
 ```
 
-where $ Fz $ is the vertical dispersion term, a function of the mixing height $ h_m $, *n* is the number of image terms, and other symbols are as defined for a Gaussian
+where $ F\_z $ is the vertical dispersion term, a function of the mixing height $ h\_m $, and other symbols are as defined for a Gaussian
 plume model.
 
+There are two methods for the mixing layer:
+- `:simplemixinglayer` uses a simple method of images, a series of reflections off of the mixing height
+- `:periodicmixinglayer` uses an infinite series of cosine terms to calculate the vertical dispersion, which is more accurate for large mixing heights
+
 # Keyword Arguments
 - `h_min=1.0`:  Minimum height for windspeed calculations.
 - `method=:simplemixinglayer`: The method used for the mixing layer.
-- `n_terms=10`: Number of image terms for the mixing layer reflection.
+- `n_terms=10`: Number terms in the series calculation of $ F_z $.
 - `mixing_limit=10_000.0`: Limit for the mixing height, in meters. Mixing heights above this are treated as infinite.
 
 # References
 + AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers
 + US EPA. 1995. *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models EPA-454/B-95-003b, vol 2*. Research Triangle Park, NC: Office of Air Quality Planning and Standards
++ Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons.
 
 """
 function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet(); 
@@ -63,6 +86,9 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet()
     hₘ = _mixing_height(scenario, eqs)
     ρₐ = _gas_density(scenario.substance,Tₐ,Pₐ)
     Qᵒ = Q*ρⱼ/ρₐ
+    if hᵣ>hₘ
+        error("Release height $hᵣ exceeds mixing height $hₘ. This is not supported by the Gaussian mixing layer model.")
+    end
 
     # max concentration
     c_max = m/Qᵒ
@@ -71,9 +97,12 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet()
     # mixing layer
     if hₘ > mixing_limit
         # infinite mixing height
+        @warn "Mixing height $hₘ exceeds limit $mixing_limit, using infinite mixing height."
         ml = SimpleVerticalTerm()
     elseif method == :simplemixinglayer
         ml = SimpleMixingLayer(n_terms, hₘ)
+    elseif method == :periodicmixinglayer
+        ml = PeriodicMixingLayer(n_terms, hₘ)
     else
         error("Unknown mixing layer method: $method")
     end
@@ -100,24 +129,28 @@ end
 Returns the solution to a Gaussian plume dispersion model with a simple reflective mixing layer.
 
 ```math
-c(x, y, z) = \frac{m_i}{u} { \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) \over { \sqrt{2\pi} \sigma_y} } F_z
+c(x, y, z) = \frac{m_i}{u} \frac{1}{ \sqrt{2\pi} \sigma_y} \exp\left(-\frac{1}{2}\left(\frac{y}{\sigma_y}\right)^2\right) F_z
 ```
-
-where $ Fz $ is the vertical dispersion term, a function of the mixing height $ h_m $, *n* is the number of image terms, and other symbols are as defined for a Gaussian
+where $ F\_z $ is the vertical dispersion term, a function of the mixing height $ h\_m $, and other symbols are as defined for a Gaussian
 plume model.
 
+There are two methods for the mixing layer:
+- `:simplemixinglayer` uses a simple method of images, a series of reflections off of the mixing height
+- `:periodicmixinglayer` uses an infinite series of cosine terms to calculate the vertical dispersion, which is more accurate for large mixing heights
+
 # Keyword Arguments
 - `downwash::Bool=false`: Include stack-downwash effects if true.
 - `plumerise::Bool=true`: Include plume-rise effects using Briggs' model if true.
 - `h_min=1.0`: Minimum height, in meters, for windspeed calculations.
 - `method=:simplemixinglayer`: The method used for the mixing layer.
-- `n_terms=10`: Number of image terms for the mixing layer reflection.
+- `n_terms=10`: Number terms in the series calculation of $ F_z $.
 - `mixing_limit=10_000.0`: Limit for the mixing height, in meters. Mixing heights above this are treated as infinite.
 
 # References
 - AIChE/CCPS. 1999. *Guidelines for Consequence Analysis of Chemical Releases*. New York: American Institute of Chemical Engineers
 - Briggs, Gary A. 1969. *Plume Rise* Oak Ridge: U.S. Atomic Energy Commission
 - US EPA. 1995. *User's Guide for the Industrial Source Complex (ISC3) Dispersion Models EPA-454/B-95-003b, vol 2*. Research Triangle Park, NC: Office of Air Quality Planning and Standards
+- Seinfeld, John H. and Spyros N. Pandis. 2006. *Atmospheric Chemistry and Physics* 2nd Ed. New York: John Wiley and Sons.
 
 """
 function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere}, 
@@ -140,6 +173,9 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
     stab = _stability(scenario)
     Γ = _lapse_rate(scenario)
     hₘ = _mixing_height(scenario, eqs)
+    if hᵣ>hₘ
+        error("Release height $hᵣ exceeds mixing height $hₘ. This is not supported by the Gaussian mixing layer model.")
+    end
 
     # jet at ambient conditions
     Tₐ = _atmosphere_temperature(scenario)
@@ -171,9 +207,12 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
     # mixing layer
     if hₘ > mixing_limit
         # infinite mixing height
+        @warn "Mixing height $hₘ exceeds limit $mixing_limit, using infinite mixing height."
         ml = SimpleVerticalTerm()
     elseif method == :simplemixinglayer
         ml = SimpleMixingLayer(n_terms, hₘ)
+    elseif method == :periodicmixinglayer
+        ml = PeriodicMixingLayer(n_terms, hₘ)
     else
         error("Unknown mixing layer method: $method")
     end

src/models/palazzi_puff.jl

@@ -37,14 +37,15 @@ downwind dispersion, σx, is calculated:
 - `: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
+- `plume_model::Type{Plume} = GaussianPlume` : the plume model $\chi$
+- additional keyword arguments are passed to the plume model
 
 # 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)
+function puff(scenario::Scenario, ::Type{Palazzi}, eqs=PalazziDefaultSet(); plume_model=GaussianPlume, disp_method=:default, kwargs...)
 
     stab = _stability(scenario)
     Δt = _duration(scenario)
@@ -55,7 +56,7 @@ function puff(scenario::Scenario, ::Type{Palazzi}, eqs=PalazziDefaultSet(); plum
         scenario,  #scenario::Scenario
         :palazzi, #model::Symbol
         disp_method, #dispersion model::Symbol
-        plume(scenario, plume_model, eqs), # plume model
+        plume(scenario, plume_model, eqs; kwargs...), # plume model
         Δt,   # duration
         u,    # windspeed
         stab, # stability class

test/models/gaussian_ml_tests.jl

@@ -13,57 +13,77 @@
                     latent_heat=1,
                     gas_heat_capacity=1,
                     liquid_heat_capacity=1)
-    rel = HorizontalJet(mass_rate=0.1,
+    hj = HorizontalJet(mass_rate=0.1,
                   duration=Inf,
                   diameter=1,
                   velocity=1,
                   height=0,
                   temperature=298.,
                   pressure=101325.,
                   fraction_liquid=0)
-    atm1 = SimpleAtmosphere(temperature=298,
-                 pressure=101325,
-                 windspeed=2,
-                 windspeed_height=1,
-                 stability=ClassF)
 
-    atm2 = SimpleAtmosphere(temperature=298,
-                 pressure=101325,
-                 windspeed=2,
-                 windspeed_height=1,
-                 stability=ClassD)
-    scn = Scenario(sub,rel,atm1)
-    
-    # test default behaviour and type inheritance
-    pl1 = plume(scn, GaussianMixingLayer)
-    @test pl1.verticalterm isa GasDispersion.SimpleVerticalTerm
+    hj_toohigh = HorizontalJet(mass_rate=0.1,
+                  duration=Inf,
+                  diameter=1,
+                  velocity=1,
+                  height=10_000,  # too high for the mixing layer
+                  temperature=298.,
+                  pressure=101325.,
+                  fraction_liquid=0)
 
-    scn2 = Scenario(sub,rel,atm2)
-    pl2 = plume(scn2, GaussianMixingLayer)
-    @test pl2.verticalterm isa GasDispersion.SimpleMixingLayer
-
-    # test vertical jet with mixing layer
-    rel = VerticalJet(mass_rate=0.1,
+    vj = VerticalJet(mass_rate=0.1,
                   duration=Inf,
                   diameter=1,
                   velocity=1,
                   height=10,
                   temperature=298.,
                   pressure=101325.,
                   fraction_liquid=0)
-    scn3 = Scenario(sub,rel,atm1)
-    pl3 = plume(scn3, GaussianMixingLayer)
-    @test pl3.verticalterm isa GasDispersion.SimpleVerticalTerm
+
+    vj_toohigh = VerticalJet(mass_rate=0.1,
+                  duration=Inf,
+                  diameter=1,
+                  velocity=1,
+                  height=10_000,  # too high for the mixing layer
+                  temperature=298.,
+                  pressure=101325.,
+                  fraction_liquid=0)
+
+    stbl = SimpleAtmosphere(temperature=298,
+                 pressure=101325,
+                 windspeed=2,
+                 windspeed_height=10,
+                 stability=ClassF)
+
+    neut = SimpleAtmosphere(temperature=298,
+                 pressure=101325,
+                 windspeed=2,
+                 windspeed_height=10,
+                 stability=ClassD)
+
+    
+    # horizontal jet test default behaviour and type inheritance
+    @test plume(Scenario(sub,hj,stbl), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleVerticalTerm
+    @test plume(Scenario(sub,hj,neut), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleMixingLayer
+    @test plume(Scenario(sub,hj,neut), GaussianMixingLayer; method=:periodicmixinglayer).verticalterm isa GasDispersion.PeriodicMixingLayer
+    @test_throws ErrorException plume(Scenario(sub,hj,neut), GaussianMixingLayer; method=:someothermethod)
+    @test_throws ErrorException plume(Scenario(sub,hj_toohigh,neut), GaussianMixingLayer)
 
-    atm3 = SimpleAtmosphere(temperature=298,
-                  pressure=101325,
-                  windspeed=2,
-                  windspeed_height=10,
-                  stability=ClassD)
-    scn4 = Scenario(sub,rel,atm3)
-    pl4 = plume(scn4, GaussianMixingLayer; downwash=true, plumerise=true)
-    @test pl4.verticalterm isa GasDispersion.SimpleMixingLayer
-    @test pl4.verticalterm.mixing_height ≈ 640.0311954829524
-    @test pl4(500, 0, 0) ≈ 1.7194314353343084e-5
+    # vertical jet test default behaviour and type inheritance
+    @test plume(Scenario(sub,vj,stbl), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleVerticalTerm
+    @test plume(Scenario(sub,vj,neut), GaussianMixingLayer).verticalterm isa GasDispersion.SimpleMixingLayer
+    @test plume(Scenario(sub,vj,neut), GaussianMixingLayer; method=:periodicmixinglayer).verticalterm isa GasDispersion.PeriodicMixingLayer
+    @test_throws ErrorException plume(Scenario(sub,vj,neut), GaussianMixingLayer; method=:someothermethod)
+    @test_throws ErrorException plume(Scenario(sub,vj_toohigh,neut), GaussianMixingLayer)
+    
+    # integration test with a scenario -- simple mixing layer
+    pl_smpl = plume(Scenario(sub,vj,neut), GaussianMixingLayer; downwash=true, plumerise=true)
+    @test pl_smpl.verticalterm.mixing_height ≈ 640.0311954829524
+    @test pl_smpl(500, 0, 0) ≈ 1.7194314353343084e-5
+
+    # integration test with a scenario -- periodic mixing layer
+    pl_per = plume(Scenario(sub,vj,neut), GaussianMixingLayer; downwash=true, plumerise=true, method=:periodicmixinglayer, n_terms=100_000_000)
+    @test pl_per.verticalterm.mixing_height ≈ 640.0311954829524
+    @test pl_per(500, 0, 0) ≈ 1.7194314353343084e-5
 
 end