gx(x,t,u,σx) = exp(-0.5*((x-u*t)/σx)^2)/(√(2π)*σx)I was thinking, recently, about how venting from vessel blowdown is modelled for screening purposes and how, more often than not, it does not take into account the blowdown curve. This is something that could easily be incorporated for simple Gaussian dispersion, which is what I examine here.
Background
The standard approach for assessing the consequences of a release from a pressure vessel is to:1
- Identify the source model (gas, liquid, aerosol)
- Calculate the mass release rate
- Model the dispersion of the release
The mass release rate from a vessel blowdown is taken as the max release rate (at the start of the blowdown) and generally assumed to be constant.2 While the standard references do acknowledge that the flow will decrease over time, this is typically not taken into account in the dispersion models. The one exception that I’m aware of is when modelling flaring due to vessel and pipeline blowdowns: sometimes an average flowrate is taken instead of the max, in which case the blowdown curve is used to derive that average. It is still a constant, though, for the purposes of dispersion modelling.
However, if we think back to the development of the Palazzi model3 for short duration releases, a rather obvious path presents itself for the special case of a release of an ideal gas from an isothermal blowdown: integrate the Gaussian puff model over time with an exponentially decaying mass release rate.
Isothermal Blowdown
Recalling the isothermal blowdown of an ideal gas, the mass release rate, \(w\), is given by
\[ w\left(t\right) = w_0 \exp \left( - \frac{t}{\tau} \right) \]
Where4
4 This follows from the definition of \(\tau\): \[ \tau = {m_0 \over w_0} \] \[ w_0 = {m_0 \over \tau} \] \[ w_0 = { {\rho_0 V} \over \tau }\]
\[ w_0 = { {\rho_0 V} \over \tau } \]
For a blowdown through an isentropic nozzle the time constant \(\tau\) is given by
\[ \frac{1}{\tau} = \frac{c_D A}{V} \sqrt{ {k P_0} \over \rho_0 } \left( 2 \over {k+1} \right)^{\frac{k+1}{2 \left( k - 1 \right)} } \]
With:
- \(c_D\) – the discharge coefficient for the blowdown
- \(A\) – the flow area of the orifice through which the blowdown is happening (e.g. a PSV)
- \(V\) – the total volume of the vessel
- \(k\) – the isentropic expansion factor, which for an ideal gas is the ratio of specific heats \(\frac{c_p}{c_v}\)
- \(P_0\) – the initial pressure in the vessel
- \(\rho_0\) – the initial density of the gas in the vessel
The Single Puff Model
For a release centred at the origin with an elevation h, the concentration profile for a single Gaussian puff is given by:5
\[ c \left(x,y,z,t \right) = w \Delta t \cdot g_x(x, t) \cdot g_y(y) \cdot g_z(z) \]
Where the gs are Gaussian functions in the x, y, and z directions
\[ g_x(x,t) = {1 \over \sqrt{2\pi} \sigma_x } \exp \left( -\frac{1}{2} \left( x-u t \over \sigma_x \right)^2 \right) \]
\[ g_y(y) = {1 \over \sqrt{2\pi} \sigma_y } \exp \left( -\frac{1}{2} \left( y \over \sigma_y \right)^2 \right) \]
gy(y,σy) = exp(-0.5*(y/σy)^2)/(√(2π)*σy)\[ g_z(z) = {1 \over \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] \]
gz(z,h,σz) = ( exp(-0.5*((z-h)/σz)^2)
+ exp(-0.5*((z+h)/σz)^2))/(√(2π)*σz)With:
- \(w\) – the constant mass release rate
- \(\Delta t\) – the duration of the release
- \(u\) – the uniform windspeed (acting only in the x direction)
- \(\sigma\)s – the dispersion parameters.
For puff releases, the dispersion parameters are typically given in reference to the centre of the cloud,6 here I have taken some puff dispersion parameters for a class D atmospheric stability.
# Puff dispersion parameters for Class D atmospheres
σx(xc) = 0.06*xc^0.92
σy(xc) = 0.06*xc^0.92
σz(xc) = 0.15*xc^0.70I like to use Unitful to manage units. This can be a little tricky with correlations, so to make that easier I use a simple macro to add a method to each correlation function mapping the correct input units and output units.
import Pkg
Pkg.add(url="https://github.com/aefarrell/UnitfulCorrelations.jl")using Unitful
using UnitfulCorrelations@ucorrel σx u"m" u"m"
@ucorrel σy u"m" u"m"
@ucorrel σz u"m" u"m"A good habit to get into, when developing code in julia, is to collect model parameters into structs. This is what I do here, collecting the parameters for a single Puff into a Puff struct.
struct Puff
m # mass
h # release height
u # velocity
t # release time
endNow I create the concentration function which takes a single puff, and a location in space and time, and returns the concentration. I also check for the special case where the puff hasn’t actually been released yet, and so does not contribute to the concentration.
Since I want this to be unit aware, both return values have to have the same units. I don’t want to hard-code this as I may also want to use this function with simple numeric types, like Float64. By using the unit function I can ensure the zero result has the same dimensions as the correct result, falling back to no units in the case where all inputs are simple numbers.
function c(p::Puff,x,y,z,t)
λ = t - p.t # time since release
xc = p.u*λ # location of cloud center
if λ > 0t
return p.m*gx(x,λ,p.u,σx(xc))*gy(y,σy(xc))*gz(z,p.h,σz(xc))
else # the puff hasn't been released yet
return 0*unit(p.m)/unit(xc)^3
end
endThe Multi-Puff Model
The single puff model assumes all of the mass is released in a single instant. This significantly over-estimates the concentration for longer duration releases, and so an alternative approach is to break up the release into several puffs and sum the result.
\[ c(x,y,z,t) = \sum_{i=0}^{n} w\left( t_i \right) \delta t \cdot g_x(x, t - t_i ) \cdot g_y(y) \cdot g_z(z) \]
Where \(\delta t\) is the duration of each puff and \(t_i\) is the time when puff i was released.
c(ps::Vector{Puff},x,y,z,t) = sum( c.(ps, x, y, z, t) );Taking the limit \(\delta t \to 0\) takes this from a discrete sum to the corresponding integral
\[ c(x,y,z,t) = \int_{0}^{t} w\left( t^{\prime} \right) \cdot g_x(x, t - t^{\prime}) \cdot g_y(y) \cdot g_z(z) dt^{\prime} \]
For the Palazzi7 model \(w(t) = w_0 \left( H\left(t - \Delta t \right) - H\left( t \right) \right)\)8 and, assuming the \(\sigma\)s are independent of time, this can be integrated to give:
8 \(H \left(t \right)\) being the Heaviside function
\[ c(x,y,z,t) = \frac{w_0}{2u} \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) \cdot g_y(y) \cdot g_z(z) \]
using SpecialFunctions: erf, erfcstruct Palazzi
w # mass release rate
h # release height
u # velocity
t_f # end of release
endfunction c(p::Palazzi,x,y,z,t)
Δt = min(t, p.t_f)
w, u = p.w, p.u
xa = u*(t-Δt)
xb = u*t
# n.b. erf(b,a) = erf(a) - erf(b)
return (w/(2u))*erf((x-xb)/(√2*σx(xb)), (x-xa)/(√2*σx(xa))) *
gy(y,σy(x))*gz(z,h,σz(x))
endc (generic function with 3 methods)A Blowdown Dispersion Model
It should be pretty obvious where I am going next: instead of assuming \(w(t)\) is a constant, let it be the exponential decay from an isothermal vessel blowdown. The integration is a little more tedious but it is not really any more difficult than the Palazzi case.
\[ c(x,y,z,t) = \int_{0}^{t} w\left( t^{\prime} \right) \cdot g_x(x, t - t^{\prime}) \cdot g_y(y) \cdot g_z(z) dt^{\prime} \]
\[ c(x,y,z,t) = \int_{0}^{t} \left[ w_0 \exp\left( -{t^{\prime} \over \tau} \right) \right] \cdot \left[ {1 \over \sqrt{2\pi} \sigma_x } \exp \left( -\frac{1}{2} \left( x-u (t - t^{\prime}) \over \sigma_x \right)^2 \right) \right] \cdot g_y(y) \cdot g_z(z) dt^{\prime} \]
Splitting this into elements that depend on time and those that don’t \[ c(x,y,z,t) = w_0 \left[ \int_{0}^{t} {1 \over \sqrt{2\pi} \sigma_x } \exp\left( -{t^{\prime} \over \tau} -\frac{1}{2} \left( x-u (t - t^{\prime}) \over \sigma_x \right)^2 \right) dt^{\prime} \right] \cdot g_y(y) \cdot g_z(z) \]
Letting everything within the integral equal \(I(x,t)\)
\[ c(x,y,z,t) = w_0 \cdot I(x,t) \cdot g_y(y) \cdot g_z(z) \]
It makes the integration a little easier to introduce \(\lambda = t-t^{\prime}\)
\[ I(x,t) = \int_{0}^{t} {1 \over \sqrt{2\pi} \sigma_x } \exp\left( -{{t - \lambda} \over \tau} -\frac{1}{2} \left( x-u \lambda \over \sigma_x \right)^2 \right) d\lambda \]
By expanding everything within the \(\exp(\dots)\), collecting terms and completing the square we arrive at:
\[ I(x,t) = {1 \over \sqrt{2\pi} \sigma_x } \exp \left( {\sigma_x^2 + 2 u \tau (x - u t)} \over {2 u^2 \tau^2} \right) \int_{0}^{t} \exp\left( -\left( {\sigma_x^2 + u \tau ( x - u \lambda ) } \over {\sqrt{2} \sigma_x u \tau} \right)^2 \right) d\lambda \]
\[ I(x,t) = \frac{1}{2u} \exp \left( {\sigma_x^2 + 2 u \tau (x - u t)} \over {2 u^2 \tau^2} \right) \left[ \mathrm{erf}\left( {\sigma_x^2 + u \tau x} \over {\sqrt{2} \sigma_x u \tau} \right) - \mathrm{erf}\left( {\sigma_x^2 + u \tau (x - u t)} \over {\sqrt{2} \sigma_x u \tau} \right)\right] \]
If we evaluate the \(\sigma_x\)s at the end points then, given that \(\sigma_x \to 0\) as \(x_c \to 0\), this simplifies to:
\[ I(x,t) = \frac{1}{2u} \exp \left( {\sigma_x^2 + 2 u \tau (x - u t)} \over {2 u^2 \tau^2} \right) \mathrm{erfc}\left( {\sigma_x^2 + u \tau (x - u t)} \over {\sqrt{2} \sigma_x u \tau} \right) \]
Giving a final concentration of:
\[ c(x,y,z,t) = \frac{w_0}{2u} \exp \left( {\sigma_x^2 + 2 u \tau (x - u t)} \over {2 u^2 \tau^2} \right) \mathrm{erfc}\left( {\sigma_x^2 + u \tau (x - u t)} \over {\sqrt{2} \sigma_x u \tau} \right) \cdot g_y(y) \cdot g_z(z) \]
struct IsothermalBlowdown
w_0 # mass release rate
τ # time constant
h # release height
u # velocity
t_f # end of release
endfunction c(p::IsothermalBlowdown,x,y,z,t)
w₀, u, τ = p.w_0, p.u, p.τ
xb = u*t
xa = t < p.t_f ? 0*xb : u*(t-p.t_f)
return (w₀/(2u))*
exp( (σx(xb)^2 + 2u*τ*(x - xb))/(2*(u*τ)^2) ) *
erf( (σx(xb)^2 + u*τ*(x - xb))/(√(2)*σx(xb)*u*τ),
(σx(xa)^2 + u*τ*(x - xa))/(√(2)*σx(xa)*u*τ) )*
gy(y,σy(x))*gz(z,h,σz(x))
endAn Example Case
Just to have something to look at, suppose an isothermal blowdown from a vessel which starts at an initial release rate of 1kg/s and the vessel contains 1000kg of an ideal gas. The vent stack is 2m above the ground and ambient windspeed is 2m/s.
# The example case
u = 2.0u"m/s"
h = 2.0u"m"
w₀ = 1.0u"kg/s"
m₀ = 1000.0u"kg"
τ = m₀/w₀The mass release rate, per above, is simply the exponential decay.
w(t) = w₀*exp(-t/τ)The total mass released by time t is simply the time-integral:
\[ m(t) = \int_0^t w_0 \exp \left( -\frac{t^{\prime}}{\tau} \right) d t^{\prime} \]
\[ m(t) = w_0 \tau \left( 1 - \exp \left( -\frac{t^{\prime}}{\tau} \right) \right) \]
m(t) = w₀*τ*(1 - exp(-t/τ))Discrete Puffs
Suppose that after \(\tau\) time has elapsed a block-valve shuts and the release abruptly ends. This release can be modelled as a series of discrete puffs by dividing the interval \([ 0, \tau )\) into \(n\) sub-intervals and release a single puff at the start of the interval with a mass \(m_i = w(t_i) \delta t\).
function discrete_puffs(;n=100, t_0=0τ, t_f=τ)
δt = (t_f - t_0)/(n-1)
pfs = Vector{Puff}()
for t_i ∈ range(t_0;stop=t_f,length=n)
m_i = w(t_i)*δt
pf = Puff(m_i,h,u,t_i)
push!(pfs,pf)
end
return pfs
endpfs = discrete_puffs(n=25);
For the purposes of illustration I chose a rather small number of puffs, as shown in Figure 1. However, if we calculate the total mass released we find it isn’t that far off.
m(pfs::Vector{Puff},t) = sum( pf.m for pf in pfs if pf.t < t );After time τ has elapsed, the total released mass is 632 kg, the total mass of the discrete puffs is 645 kg, an excess of only 2.1%.
Comparing Results
It is worth checking how the implementation of the approximate integral compares. Consider a point 500m downwind of the vent stack, at the same release height as the stack. Suppose a sampling station is set up. What would the concentration profile look like?
I think it is reasonable to want to compare this to the Palazzi model – it should already be clear that a single large puff would be a terrible approximation for a release that goes on for several minutes. We expect the concentration to be bounded between the Palazzi case when the mass rate is a constant at \(w_0\) and when the mass rate is \(w(\tau)\). Furthermore it should connect the two curves with something resembling an exponential decay.
bd = IsothermalBlowdown(w₀,τ,h,u,τ)
The results are showin in Figure 2 above, which matches our expectations. The approximate integral model developed here is virtually identical to the discrete puffs model with 100 puffs. For comparison I also included the case where the Palazzi model is used but with a time-averaged constant release rate. This will have the correct total mass in the release, but clearly underestimates the peak concentration.
The ground level concentration also conforms with our expectations, as shown in Figure 3. The region around the vent itself, besides having some artifacts of the discretization and marching squares, is likely quite unreliable. This is the region where the fundamental assumptions – that the releases have zero momentum and buoyancy, are most egregiously violated. I think this model is still reasonable for concentrations far enough from the vent that the windspeed dominates the advection, though an effective release point would need to be used.
A Note on Sources
It is the nature of the universe that the instant I post this I will find where this model was published in the literature. I haven’t found it yet, but I can’t imagine I am first person to come with this. Knowing me, it is probably in one of the references I look at all the time and, somehow, failed to notice.
If this paragraph is still here when you see this, and you know of a published reference for this model, please leave a comment.
References
Center for Chemical Process Safety. Guidelines for Consequence Analysis of Chemical Releases. New York: Center for Chemical Process Safety/AIChE, 1999.
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. https://www.researchgate.net/publication/328744668_DIFFUSION_FROM_A_STEADY_SOURCE_OF_SHORT_DURATION.