Notebook Version 2: This notebook has been updated to reflect changes made to the package DifferentialMobilityAnalyzers.jl to work with the Julia v1 series (tested with Julia 1.1.0). To read the original supplement published with the paper please switch to v1.0.0 of the package DifferentialMobilityAnalyzers.jl and/or download the virtual machine on zenodo.org which contains a complete installation that works with Julia 0.6.4
This notebook introduces the Fredholm integral equation and derives the discretized solution via the forward convolution matrix. The notebook is a supplement to the manuscript
Petters, M. D. (2018) A language to simplify interpretation of differential mobility analyzer response functions, Aerosol Science & Technology, 52 (12), 1437-1451.
Figure 1 shows a schematic instrument setup with an array of detectors. The DMA is operated in classifier mode and set to a single voltage. The mobility classified particles are then passed to one or more detectors.
At constant v the DMA selects mobility "z star" ($z^s$) based on DMA geometery and voltage (also see Notebook 1).
The relationship between physical diameter $D[z^s,k]$ and number of charges is (also see Notebook 1)
which leads to the computational identity
Transfer through the DMA for k charged particles is given by
where Z is a vector of transmitted mobilities. For k > 1 $\Omega(Z, z_k^s)$ selects particles of lower mobility and thus larger size. The net transfer through the DMA can thus be written as
In this form, the net transmission vector through the DMA is readily computed. Table 1 summarizes these functions and demonstrates how they are represented in the code.
Quantity | Code | Reference | Notes |
---|---|---|---|
$D[Z,k]$ | ztod(Λ,k,Z) | Nb.S1.B3c | Nb.S1.B3c means compare Notebook S1 Block 3c |
$T_l(D[Z,k])$ | δ.Tl(Λ,D[Z,k]) | Nb.S1.B3d | d is in nm and Λ is a structure that contains flow rates |
$T_c(D[Z,k])$ | δ.Tc(k,D[Z,k]) | Nb.S1.B3e | |
$\Omega(z,z_k^s)$ | δ.Ω(Λ,z,$z_k^s$) | Nb.S1.B3f | with $z^s$ = vtoz(Λ,v) and $z_k^s = \frac{z^s}{k}$ |
A step-by-step illustration of the transmission function $T$ is presented next.
using Plots, Plots.PlotMeasures, DataFrames, DifferentialMobilityAnalyzers
plotlyjs();
┌ Info: Recompiling stale cache file /home/mdpetter/.julia/compiled/v1.1/Plots/ld3vC.ji for Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80] └ @ Base loading.jl:1184 ┌ Info: Recompiling stale cache file /home/mdpetter/.julia/compiled/v1.1/DataFrames/AR9oZ.ji for DataFrames [a93c6f00-e57d-5684-b7b6-d8193f3e46c0] └ @ Base loading.jl:1184 ┌ Info: Precompiling DifferentialMobilityAnalyzers [050a7be8-18f9-52d6-bdba-0a194485d082] └ @ Base loading.jl:1186 ┌ Info: Recompiling stale cache file /home/mdpetter/.julia/compiled/v1.1/PlotlyJS/1r9Ld.ji for PlotlyJS [f0f68f2c-4968-5e81-91da-67840de0976a] └ @ Base loading.jl:1184 WARNING: could not import Base.quit into AtomShell
The Λ structures define the DMA in terms of flow rates and geometry and power supply polarity. Geometry parameters correspond to the TSI 3080 long column. The effective diffusion length describes diffusional loss in the DMA column (Reineking & Porstendörfer, 1986) which is described further below. Diffusional loss is ignored if leff = 0
. The Λ variable is used in various functions to define the state of the DMA.
t,p = 296.15, 1e5 # Temperature [K], Pressure [Pa]
qsa,β = 1.66e-5, 1/5 # Qsample [m3 s-1], Sample-to-sheath ratio,
r₁,r₂,l = 9.37e-3,1.961e-2,0.44369 # DMA geometry [m]
leff = 13.0 # DMA effective diffusion length [m]
m = 6 # Upper number of charges
DMAtype = :cylindrical # Select DMA geometry
Λ = DMAconfig(t,p,qsa,qsa/β,r₁,r₂,l,leff,:-,m,DMAtype) # Specify DMA with negative polarity
bins,z₁,z₂ = 512, dtoz(Λ,1000e-9), dtoz(Λ,10e-9) # bins, upper, lower mobility limit
δ = setupDMA(Λ, z₁, z₂, bins); # Compute matrices
figure("Nimbus Sans L", 2, 5, 2, 8) # Figure setup for all figures below
zˢ = dtoz(Λ,100e-9) # Pick a mobility and evaluate over grid
p1 = plot(δ.Z, δ.Ω(Λ,δ.Z,zˢ), xaxis = :log10, legend = :none, color = :black, left_margin = 10px,
xlabel = "Electrical mobility (m² V⁻¹ s⁻¹)", ylabel = "Transmission probablity (-)")
p1 = plot!(δ.Z, δ.Ω(Λ,δ.Z,zˢ/2), color = RGBA(0.8,0,0,1))
p1 = plot!(δ.Z, δ.Ω(Λ,δ.Z,zˢ/3), color = RGBA(0,0,0.8,1))
p2 = plot(δ.Dp, δ.Ω(Λ,δ.Z,zˢ), xaxis = :log10, legend = :none, color = :black, left_margin = 10px,
xlabel = "Particle Diameter (nm)", ylabel = "Transmission probablity (-)")
p2 = plot!(δ.Dp, δ.Ω(Λ,δ.Z,zˢ/2), color = RGBA(0.8,0,0,1))
p2 = plot!(δ.Dp, δ.Ω(Λ,δ.Z,zˢ/3), color = RGBA(0,0,0.8,1))
plot(p1,p2, grid = (2,1), fmt = :svg)
The charge filter function computes the charge fraction as a function of diameter and number of charges k. Note that the diameters δ.Dp
map directly to the mobility vector Z
assuming k = 1 charge (i.e. Dp = ztod(Λ,1,Z)
).
figure("Nimbus Sans L", 2, 6, 2, 8)
k = 1
p1 = plot(δ.Z, δ.Tc(k,δ.Dp), xaxis = :log10, label = "δ.Tc(k=$k, Dp)", color = :black, left_margin = 25px,
xlabel = "Electrical mobility (m² V⁻¹ s⁻¹)", ylabel = "Charging probablity (-)")
k = 2
p1 = plot!(δ.Z, δ.Tc(k,δ.Dp), xaxis = :log10, label = "δ.Tc(k=$k, Dp)", color = RGBA(0.8,0,0,1))
k = 3
p1 = plot!(δ.Z, δ.Tc(k,δ.Dp), xaxis = :log10, label = "δ.Tc(k=$k, Dp)", color = RGBA(0,0,0.8,1))
k = 1
p2 = plot(δ.Dp, δ.Tc(k,δ.Dp), xaxis = :log10, label = "δ.Tc(k=$k, Dp)", color = :black, left_margin = 25px,
xlabel = "Diameter (nm)", ylabel = "Charging probablity (-)")
k = 2
p2 = plot!(δ.Dp, δ.Tc(k,δ.Dp), xaxis = :log10, label = "δ.Tc(k=$k, Dp)", color = RGBA(0.8,0,0,1))
k = 3
p2 = plot!(δ.Dp, δ.Tc(k,δ.Dp), xaxis = :log10, label = "δ.Tc(k=$k, Dp)", color = RGBA(0,0,0.8,1))
plot(p1,p2, grid = (2,1), fmt = :svg)
The diffusional transmission filter computes the penetration efficiency as a function of diameter and number of charges k. Note that the diameters δ.Dp
map directly to the mobility vector Z
assuming k = 1 charge (i.e. Dp = ztod(Λ,1,Z)
).
T = (zˢ,k,Λ,δ) -> δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tc(k,δ.Dp).*δ.Tl(Λ,δ.Dp); # DMA Transmission (Notebook 2)
p1 = plot(δ.Z, T(zˢ,1,Λ,δ) , xaxis = :log10,label = "T(zˢ,1)", color = :black,
xlabel = "Electrical mobility (m² V⁻¹ s⁻¹)", ylabel = "Net transmission prob. (-)")
p1 = plot!(δ.Z, T(zˢ,2,Λ,δ), xaxis = :log10, label = "T(zˢ,2)", color = :red)
p1 = plot!(δ.Z, T(zˢ,3,Λ,δ), xaxis = :log10, label = "T(zˢ,3)", color = :blue)
p2 = plot(δ.Dp, T(zˢ,1,Λ,δ), xaxis = :log10, label = "T(zˢ,1)", color = :black,
xlabel = "Diameter (nm)", ylabel = "Net transmission prob. (-)")
p2 = plot!(δ.Dp, T(zˢ,2,Λ,δ), xaxis = :log10, label = "T(zˢ,2)", color = :red)
p2 = plot!(δ.Dp, T(zˢ,3,Λ,δ), xaxis = :log10, label = "T(zˢ,3)", color = :blue)
plot(p1,p2, grid = (2,1), fmt = :svg)
Figure 4. Net transmission probablity due to DMA transfer function, charging probability, and transmission probability expressed as the convolution function
T = (zˢ,k,Λ,δ) -> δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tc(k,δ.Dp).*δ.Tl(Λ,δ.Dp)
for k = 1,2,3.
A size distribution is obtained by systematically stepping through voltages (or zˢ). Figure 5 shows an example setup where concentration at the exit of the DMA is measured with a condensation particle counter (CPC), here assumed to operate with unit efficiency. During the experiment the voltage is scanned to cover a range in electrical mobility. The response function is binned, resulting in the vector $R$, where $R_i$ is the response in instrument channel $i$ and is expressed as number concentration per unit volume. The response vector is sometimes plotted against the mobility diameter ($D[Z,1]$) assuming singly charged particle. The response vector can be computed multiplication of the true size distribution vector $N$ with a convolution matrix $\mathbf{A}$ derived below.
The DMA response vector is obtained via the Fredholm integral equation of the first kind (e.g. Kandlikar and Ramachandran, 1999, Stolzenburg and McMurry, 2008, Petters et al. 2007,2009).
The integral is performed over the limits $z_a$ and $z_b$, which corresponds to the upper and lower mobility limit set by the voltage range used to operate the DMA. The function $\frac{dN}{d\ln D}\frac{d\ln D}{dz}dz$ evaluates to the number concentration of particles that lie in the interval $[z,z + dz]$. Note that $D[z,1]$ is used in the charge filter and loss function since the integral is applied over the transform of the selected centroid mobility $Z_{i,k}^s$ (see above). $Z$ is the mobility vector of the grid and the subscript $i$ denotes the channel response channel. Finally, $\epsilon_i$ is the detector error in channel $i$. The convolution integral can be discretized:
which is equivalent to
t,p = 296.15, 1e5 # Temperature [K], Pressure [Pa]
qsa,β = 1.66e-5, 1/5 # Qsample [m3 s-1], Sample-to-sheath ratio,
r₁,r₂,l = 9.37e-3,1.961e-2,0.44369 # DMA geometry [m]
leff = 13.0 # DMA effective diffusion length [m]
m = 6 # Upper number of charges
DMAtype = :cylindrical
Λ = DMAconfig(t,p,qsa,qsa/β,r₁,r₂,l,leff,:-,m,DMAtype) # Specify DMA with negative polarity
bins,z₁,z₂ = 512, dtoz(Λ,1000e-9), dtoz(Λ,10e-9) # bins, upper, lower mobility limit
δ = setupDMA(Λ, z₁, z₂, bins); # Compute matrices
𝕟 = DMALognormalDistribution([[230, 50, 1.3], [280, 140, 1.6]],δ)
figure("Nimbus Sans L", 2, 3.5, 1.75, 8)
p1 = plot(𝕟.Dp, 𝕟.N, xaxis = :log10, xticks = [10, 100, 1000], legend = :none, right_margin = 35px,
ylabel = "Concentration (cm⁻³)", xlim = (8,1000), color = :black, xlabel = "Diameter (nm)", lt = :steppre)
p2 = plot(𝕟.Dp, δ.𝐀*𝕟.N, legend = :none, xaxis = :log10, right_margin = 25px, lt = :steppre,ylim = (0,8.2),
xlim = (8,1000), color = :black, xlabel = "Diameter (nm)")
plot(p1,p2,grid=(1,2), fmt = :svg)
This section demonstrates how the convolution matrices are computed in DifferentialMobilityAnalyzer.jl. For a given z vector, the matrix can be computed essentially with a single command and the user defined function Σ
Σ = (f,i) -> mapreduce(f,+,1:i)
T = (zˢ,k,Λ,δ) -> δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tc(k,δ.Dp).*δ.Tl(Λ,δ.Dp)
𝐀=(hcat(map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)...))'
The segements below explain the commands step-by-step. As demonstrated above
T = (zˢ,k,Λ,δ) -> δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tc(k,δ.Dp).*δ.Tl(Λ,δ.Dp)
selects for a single mobility zˢ and charge k. The function
Σ = (f,i) -> mapreduce(f,+,1:i)```
sums functions that take a single argument and apply it for the array 1:i. For example Σ(x->x^2,4) = 1^2 + 2^2 + 3^2 + 4^2 = 30 or Σ(exp,2) = exp(1) + exp(2) = 10.107. In the context above
```julia
y = Σ(k->T(zˢ,k,Λ,δ),3)
will compute T(zˢ,1,Λ,δ)+T(zˢ,2,Λ,δ)+T(zˢ,3,Λ,δ)
and return an array of transmission efficiencies.
figure("Nimbus Sans L", 2, 5, 2, 8)
y = Σ(k->T(zˢ,k,Λ,δ),3)
p1 = plot(δ.Z, y, xaxis = :log10,
label = "∑T(zˢ,k)", color = :black, xlabel = "Electrical mobility (m² V⁻¹ s⁻¹)",
ylabel = "Net transmission prob. (-)")
p2 = plot(δ.Dp, y, xaxis = :log10, label = "∑T(zˢ,k)", color = :black,
xlabel = "Diameter (nm)", ylabel = "Net transmission prob. (-)")
plot(p1,p2, grid = (2,1), fmt = :svg)
The command
map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)```
computes the above curve for all zˢ/k ∈ Z, i.e. each bin midpoint from the z-array The length of Z is number of bins. For each Z value, the expression produces an array with length equals to the number of bins. The resulting output is an array of arrays.
x = map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)
println(typeof(x)) # What data type
println(size(x)) # Length of x
println(size(x[1])) # Length of x[1], x[2], ...
Array{Array{Float64,1},1} (512,) (512,)
To demonstrate how arrays of arrays work, here is an example.
test = Array{Float64}[]
push!(test,[1,2,3])
push!(test,[5,6,4])
test
2-element Array{Array{Float64,N} where N,1}: [1.0, 2.0, 3.0] [5.0, 6.0, 4.0]
In julia, the splatting operator will pass all elements of an array to a function. For example
fun = (a,b,c)->a+b+c # A function with three inputs
x = [1,2,3] # an array of three numbers
fun(x...) # fun(x...) = fun(x[1],x[2],x[3])
6
Applied to the test array:
hcat(test...) # hcat(test...) = hcat(test[1], test[2]). Note that hcat takes any number of inputs
3×2 Array{Float64,2}: 1.0 5.0 2.0 6.0 3.0 4.0
Above
x = map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)```
x[1], ... x[512] are the rows of the convolution matrix. To arrange the Array{Array} in the form of a matrix, hcat + splatting is used. Finally the ()' transposes the matrix to get it into the correct orientation.
𝐀 = (hcat(map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)...))'
512×512 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}: 0.100079 0.0937158 0.0863586 … 0.0 0.0 0.0934605 0.100466 0.0940936 3.61387e-14 1.96406e-73 0.0861216 0.0938279 0.100865 3.64087e-14 1.97872e-73 0.0788909 0.0864601 0.0942062 3.66802e-14 1.93093e-34 0.0717671 0.079201 0.0868086 0.0 3.6366e-14 0.0647484 0.0720492 0.0795203 … 2.07541e-99 3.66372e-14 0.0578333 0.0650029 0.0723396 3.75078e-14 1.9745e-34 0.0510203 0.0580606 0.065265 2.43823e-19 0.0 0.0443079 0.0512209 0.0582947 0.0 2.06901e-73 0.0376946 0.0444821 0.0514274 3.8354e-14 2.01904e-34 0.031179 0.0378428 0.0446614 … 0.0 3.80252e-14 0.0247595 0.0313015 0.0379954 0.0 2.0493e-34 0.0184348 0.0248568 0.0314277 0.0 0.0 ⋮ ⋱ ⋮ 9.11176e-17 8.51466e-17 5.86685e-18 … 0.00585187 0.0046835 0.0 6.02804e-18 0.0 0.00701823 0.00577335 1.00769e-16 0.0 3.05663e-17 0.00819643 0.00691461 6.23778e-18 6.26498e-17 3.0558e-17 0.00934198 0.00806617 0.0 5.69534e-17 1.42868e-16 0.0104112 0.00918468 6.33125e-18 6.21005e-18 1.76263e-16 … 0.011364 0.0102277 5.63106e-17 1.58167e-17 5.67598e-17 0.0121668 0.0111564 0.0 0.0 0.0 0.0127943 0.0119383 5.71544e-17 3.1693e-17 8.93261e-17 0.0132302 0.0125489 8.94844e-17 8.97369e-17 0.0 0.0134681 0.0129727 0.0 5.82389e-17 1.50237e-17 … 0.0135107 0.0132035 5.84439e-17 3.2408e-17 0.0 0.0133687 0.0132443
The full transmission function through the DMA is given by
δ.Ω(Λ,δ.Z,i/k).*δ.Tc(k,δ.Dp).*δ.Tl(Λ,δ.Dp)
Modifying this expression allows constructing other matrices. For example, consider the case where particles with known mobility are supplied to a second DMA. There is no bipolar charger. In this case the function δ.Tc(k) is not used. Since nothing is known about the particle mass/charge ratio, and the split within the population, all particles are assumed to be +1 charged. The only error in this assumption is that the multiply charged particles hiding in the mobility peak may have a different loss rate Tl, which does not vary strongly with size except for very small particles. To describe this process, the transmission matrix O is defined as
𝐎 = (hcat(map(i->Σ(k->δ.Ω(Λ,δ.Z,i/k).*δ.Tl(Λ,δ.Dp),1),δ.Z)...))'
which is similar to 𝐀 but with Λ.n = 1 and δ.Tc(k) omitted. The matrix 𝐎 is used to describe size selection in tandem DMA systems where no charge neutralization is used between DMA1 and DMA2 (see Notebooks S8,S9).
𝐎 = (hcat(map(i->Σ(k->δ.Ω(Λ,δ.Z,i/k).*δ.Tl(Λ,δ.Dp),1),δ.Z)...))'
512×512 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}: 0.961405 0.896743 0.823017 … 0.0 0.0 0.897824 0.961336 0.896734 8.68042e-13 0.0 0.827322 0.897815 0.961267 8.74522e-13 0.0 0.757862 0.827314 0.897807 8.8105e-13 0.0 0.689427 0.757855 0.827306 0.0 8.80809e-13 0.622002 0.68942 0.757847 … 0.0 8.87384e-13 0.555573 0.621996 0.689414 9.00927e-13 0.0 0.490124 0.555568 0.62199 0.0 0.0 0.425642 0.49012 0.555562 0.0 0.0 0.362112 0.425638 0.490115 9.21253e-13 0.0 0.299519 0.362108 0.425634 … 0.0 9.21001e-13 0.237851 0.299516 0.362105 0.0 0.0 0.177093 0.237849 0.299514 0.0 0.0 ⋮ ⋱ ⋮ 0.0 0.0 0.0 … 0.14056 0.113438 0.0 0.0 0.0 0.168576 0.139835 0.0 0.0 0.0 0.196876 0.167477 0.0 0.0 0.0 0.224392 0.195369 0.0 0.0 1.0659e-15 0.250073 0.22246 0.0 0.0 1.07385e-15 … 0.27296 0.247722 5.40945e-16 0.0 5.40934e-16 0.292244 0.270217 0.0 0.0 0.0 0.307315 0.289155 5.49051e-16 0.0 5.4904e-16 0.317785 0.303944 5.53149e-16 5.53144e-16 0.0 0.3235 0.314208 0.0 5.57273e-16 0.0 … 0.324524 0.319801 5.61438e-16 0.0 0.0 0.321112 0.320787
Kandlikar, M., & Ramachandran, G. (1999) Inverse methods for analysing aerosol spectrometer measurements: A critical review. Journal of Aerosol Science, 30(4), 413-437, DOI: 10.1016/S0021-8502(98)00066-4.
Petters, M. D., A. J. Prenni, S. M. Kreidenweis, P. J. DeMott (2007) On measuring the critical diameter of cloud condensation nuclei using mobility selected aerosol, Aerosol Science and Technology, 41(10), 907-913, doi:10.1080/02786820701557214.
Petters, M. D., C. M. Carrico, S. M. Kreidenweis, A. J. Prenni, P. J. DeMott, J. R. Collett Jr., and H. Moosmüller (2009) Cloud condensation nucleation ability of biomass burning aerosol, Journal of Geophys Research, 114, D22205, DOI:10.1029/2009JD012353.
Reineking, A. & J. Porstendörfer (1986) Measurements of Particle Loss Functions in a Differential Mobility Analyzer (TSI, Model 3071) for Different Flow Rates, Aerosol Science and Technology, 5:4, 483-486, DOI: 10.1080/02786828608959112.
Stolzenburg, M. R. & P. H. McMurry (2008) Equations Governing Single and Tandem DMA Configurations and a New Lognormal Approximation to the Transfer Function, Aerosol Science and Technology, 42:6, 421-432, DOI: 10.1080/02786820802157823