- Email: [email protected]
An inversion method for the determination of the particle size distribution from diffusion battery measurements
Pergamon
J.Aerosol
SCI , Vol. 26, No. 5, pp. 797X312, 1995 Cowrinht 0 1995 Elsevier Science Ltd Printed’,; &at Bntain. All rights reserved CNl21~8502/95 $9.50 + 0.00
0021-8502(95)00007-0
AN INVERSION METHOD FOR THE DETERMINATION OF THE PARTICLE SIZE DISTRIBUTION FROM DIFFUSION BATTERY MEASUREMENTS D. Lesnic, L. Elliott Department
of Applied
Mathematical
Studies,
and D. B. Ingham University
of Leeds, Leeds LS2 9JT, U.K.
(First received 8 September 1994; and in final form 3 January 1995)
Abstract-The inverse problem of the determination of the aerosol particle size distribution from a set of deposition measurements taken on the screens of the stages of a diffusion battery is analysed. The mathematical formulation reduces the problem to solving a system of integral equations for the unknown particle size distribution. Since the problem is ill-posed a numerical method which is based on minimising the norm of the derivative function subject to simple bounds of the variables and certain linear constraints, and which does not use any a priori information regarding the parametrised form of the solution is employed. For single- and double (normal or lognormal)-mode particle size distributions the numerical method provides a good estimation of the corresponding analytical solutions provided that, in general, the problem is approximated as a continuous Fredholm integral equation of the first kind. For triple-mode distributions a sequential application of the numerical method provides results which give good estimates for the first two modes but the prediction of the third mode is only satisfactory.
1. INTRODUCTION
An important problem in atmospheric aerosol physics relates the size distribution of fine aerosol particles (diameter less than 1 pm) to a set of measurements obtained using a particle sizing instrument. These instruments include diffusion batteries, which consist of a set of rectangular or circular channels (channel-type), or fibrous filters (screen-type) through which an aerosol is passed. At this stage it is worth mentioning that the analysis performed in this paper can in principle be applied to other particle size instruments such as cascade impactors, variable slit impactors, optical counters, electrical mobility analyser, etc., and will be considered in future investigations. Although in all these instruments there is only a limited number of data points which are available to obtain a particle size distribution, the situation appears more difficult for the diffusion battery case since the diffusion battery responses, i.e. the deposition kernels, are, in general, are not sharp and may be non-monotonic (see Wu et al., 1989). It is this unusual behaviour of the diffusion battery kernels compared with those in other spheres of interest that necessitates the more elaborate inversion procedures discussed herein to retrieve the particle size distribution function. If an inversion method developed for a non-monotonic response instrument is shown to perform well, then for a monotonic response instrument the method will be expected to perform even better. In addition, for the diffusion battery, the existence of an analytical form for the deposition kernels enables a more accurate and simple development of the inversion method, whilst for most other instruments the deposition kernels are obtained empirically, hence requiring much interpolation, resulting in the introduction of additional noise into the problem. Some of the previous studies were concerned with the design and calibration of diffusion batteries (see Sinclair and Hopes, 1975; Cheng et al., 1980; Baklanov and Dubtsov, 1993). Also, many other studies, both experimental and theoretical, have dealt with the determination of the deposition functions based on the analysis of particle collection efficiency in channels or filters (see, for example, Emi et al., 1982; Cheng and Yeh, 1983; Cheng et al., 1985). When the particle size distribution is known in advance then one can evaluate directly the deposition data (see Skaptsov and Lipatov, 1994). However, in practice, only the measurements of the fraction of particles collected at each state through the diffusion 797 AS 26-5-G
D. Lesnic et al.
798
battery are available, and this gives rise to a typical inverse problem for the determination of the unknown particle size distribution, and which still remains an unsolved problem (see Yee, 1989). Given a set of fractional measurements, di, for i = 1, . . . , n, and a set of deposition functions, Ki(r), for i = 1, . . . , n, it is desirable to retrieve, as far as possible, the unknown continuous particle size distribution function, f(r), which satisfies the integral equations di =
Im*= Ki(r)f(r)dr, s hi”
i = 1, . . . , n,
(1)
where [rmin, rmax] is the particle size (radius) range of interest. It is clear that the solution of the system of integral equations (1) is not unique, since there are an infinite number of continuous size distribution functions f(r) which will yield the discrete set of deposition measurements, di, for i = 1, . . . , n. In addition, since there are errors associated with each measurement, equations (1) are satisfied only within a tolerance Ei, for i = 1, . . . , n, namely lms. d!“’= d.I + s.I = Ki(r)f@) (r) dr > i=l 3*.. 9n. (2) s hi” The small quantities si, for i = 1, . . . ,n, are introduced in order to simulate the inherent experimental measurement errors that may cause large deviations of the solution f’“‘(r) of the system of equations (2) from its exact solution f(r) of the system of equations (1) (see Wolfenbarger and Seinfeld, 1990), i.e. there is a possibility of a discontinuous dependence of the solution upon the input data. The non-uniqueness and any discontinuous dependence upon the data show that the inverse problem (1) is ill-posed. However, most of the many solutions of problem (1) are physically unacceptable because they are highly oscillatory and include negative values. If a numerical procedure produces negative values for the solution, then simply replacing them with zero values, as suggested by Bashurova et al. (1991), might not be the best option in certain cases. In order to overcome these difficulties, additional physical constraints may be added to the set of equations (1) or (2) in order to restrict the class of possible solutions. In particular, we enforce that the particle size distribution be a positive function, i.e.
rmaxl~ r E Crdn,
f(r) 2 0,
(3)
and also a normalised condition is required which, without loss of generality, may be written as (4) Nevertheless, the natural constraints (3) and (4) are still too weak to enable problem (1) to have a unique solution. Various numerical methods have been proposed for the solution of the problem (l)-(4). These methods include: manual techniques based on trial function analysis (see Soderholm, 1983), regularisation procedures (see Twomey, 1963), numerical filtering (see Twomey, 1965), nonlinear iterative methods (see Twomey, 1975), nonlinear programming using physical constraints (see Cooper and Spielman, 1976), maximal entropy production (see Yee, 1989), extreme value estimation (see Paatero, 1990), etc. Most of these studies were concerned in dealing with the discontinuous dependence upon the data, i.e. stabilising the results by imposing smoothing constraints. However, when dealing with the response of a diffusion battery the highly ill-posed nature of the mathematical problem is mainly due to the non-unicity of the solution rather than to the discontinuous dependence on the data for which a regularisation procedure will offer, in general, a stable solution (see Tikhonov and Arsenin, 1977). A priori parameter estimations of the form of solution have been proposed by assuming that the particle size distribution is a superposition of lognormal kernels (see Heintzenberg, 1975): fcr)
=
i
k=l
hk
exp
r
Determination of the particle size distribution
799
or normal kernels (see Lloyd et al., 1994); f(r) = 2 hk exp -@ - I&)* ( XT-> k=l
(6)
where hkrpk and a&.,for k = 1, . . . ,p, represent the amplitude, the mean and the standard deviation, respectively. In this way the infinite-dimensional space in which the solutions may lie is reduced to a finite-dimensional space but pays the penalty that the solution is a priori assumed analytical, which implies a degree of smoothness which is not satisfied by the real particle size distributions. In addition, other types of possible particle size distributions such as Junge type or Gamma distributions will never be recovered. It is the purpose of this study to introduce a numerical method which does not use any a priori information regarding the form of the solution, i.e. the parameter estimation given by equations (5) or (6), and which is based only on a modified first-order regularisation approach in order to achieve the stability of the results. Also, the degree of non-uniqueness of the problem may be reduced by increasing the number of equations under a piecewise linear approximation of the deposition data and a suitable extension of the deposition functions. This can easily be performed by joining the discrete experimental data by straight lines. Furthermore, a successive application of the numerical method enables good estimations of the available analytical solutions, up to the first two peaks in a triple-mode distribution, to be obtained.
2. DIFFUSION
BATTERY DESIGN
AND DEPOSITION
FUNCTIONS
The diffusion battery design is similar to that analysed by Lloyd et al. (1994), i.e. it consists of eight stages of stacks of 400 mesh woven wire screens, with a ninth stage which is an absolute filter, i.e. there are no particles passing through it. Particles with a radius smaller than r,,,i,, are all deposited on the first screen, whilst particles with a radius greater than rmaxare removed by an impactor at the diffusion battery inlet. In this section the input data which is necessary for the inverse analysis developed in Section 3 is defined. The deposition function, Ki(r), as a function of particle radius, r, on a given stage i, for i=l , . . . ,9, is the difference between the penetration through all stages including itself, Pi, and the penetration up to that stage, Pi-l, namely K,(r) = 1 - PI(r), Ki(r) = Pi_ l(r) - Pi(r), K9(r1
=
(74
i = 2, . . . ,8,
(7b)
(7c)
P8W
The fractional penetration of aerosols at stage i, Pi(r), according to the fan model filtration theory of the wire screen diffusion battery (see Cheng and Yeh, 1980) is expressible as Pi(r) = exp
-
2niah q(r) , 7c(l - a)u >
i=l,...,
8,
where u = 0.3116 is the solid volume fraction of the screen, ni = 2’-’ is the cummulative number of screens from the first stage down to stage i, for i = 1, . . . ,8, h = 76.9 pm is the screen thickness and a = 14.7 pm is the wire radius. The single wire collection efficiency, q(r), can be expressed as the sum of the individual efficiencies due to diffusion, direct interception, inertial impaction and a correction term due to diffusion and interception (see Kirsch and Stechkina, 1978), namely v(r) = v&) + vR(r) + 11(r) + IfDR(r),
(9)
where qD(r) = 2.7 Pe(r)-*13, qR(r) = (2x)-‘[(l
+ R(r))-’ - 1 - R(r) + 2(1 + R(r))ln(l + R(r))],
Wa) (lob)
D. Lesnic et al.
800
q,(r) = (2x)-’ Stk(r) C(29.6 - 28cc0.62)R(r)2 - 27.5R(r)‘.*], qDR(r) = 1.24~~ ‘I2 Pe(r)- l/2 R(r)213,
(1W Wd)
where Pe(r) = 2aU,/D(r) is the Peclet number, U o = 4Q/nd: is the freestream fluid veloflow rate, dr = 3.5 cm is the screen diameter, city, Q = 0.44 lmin- r is the volumetric diffusion coefficient for a spherical particle, D(r) = kTC(r)/6qw is the k = 1.38054x lo-r6ergK-’ is the Boltzmann constant, T = 293 K is the absolute temperature, C(r) = 1 + (1.246 + 0.42 exp (- 0.87 r/A)) A/r is Cunningham slip correction factor, A = 0.0665 pm is the mean free path of the gas molecules, ,U = 0.000184 P is the fluid viscosity, R(r) = r/a is the interception parameter, x = -0.5 In(c) - 0.75 + c - 0.25c2, c = 2c+r, Stk(r) = ppUOr2C(r)/(9pa) is the Stokes number and pp is the particle density. Since only fine aerosol particles are investigated the particle radius range of interest is taken as [r,in, rmax] = [5 nm, 500 nm] and therefore the inertial term, q,(r), in expression (9) can be neglected.
3. NUMERICAL
METHOD
Discretising the interval [rmin, r,,,] into m equidistant intervals (i.e. elements), [rj_ r, rj], for j = 1, . . . , WI,and assuming the function f(r) to be constant over each element and taking the value at its midpoint, r; (i.e. nodes), equations (l), (3) and (4) become di = 2 j=l
where fj =f(r$)
Kijfj,
i = l,...,n,
(11)
fj20,
j=l,...,
m,
(12)
j$l(rj
- rj-l)fj
= 1,
(13)
Ki(r)dr,
i=
and the matrix IJ Kij
=
s
l,...,
n, j=
l,...,
m,
(14)
rj-I
is evaluated using a trapezoidal quadrature rule. Better quadrature rules, e.g. Simpson rule, can be used but it was found not to have any significant effect on the final results. Extending the matrix K by defining Kc,+ l)j = rj - rj_ 1, for j = 1, . . , m, and the vector d by defining d n+1 = 1, equations (11) and (13) result in a linear system of (n + 1) equations with m unknowns f = (fj), for j = 1, . . . , m, namely d=Kf which has to be solved subject
to the positive
(15) constraint
condition
(12), namely
f>O. Similarly, form
when noise is present
in equations
(2) the resulting
(16) system of equations
d(E) = Kf’“’ 9 where d = (dy)) for i = l,...,(n If m = II + 1 and the matrix
+ 1) and dfi, = 1. K is non-singular, then one would f=K-‘d.
takes the
(17) simply have (18)
However, because the problem is ill-posed the system of equations (15) is highly illconditioned (see Cooper and Wu, 1990), and therefore the direct method based on the inversion of equation (15) cannot be applied. If m < n + 1 = 10, the approximation of the integral equations (1) with the discrete version (11) would be nevertheless unrealistic for functions which rapidly change their
Determination
of the particle
size distribution
801
magnitude. In addition, it is desirable to obtain a continuous function rather than at a few discrete points and therefore m should be large, much greater than n + 1. Alternatively, one may still utilise the least-squares method based on minimising the norm J/d - KflJ’, which has the solution f = (K”K)-‘Pd.
(19)
However, it was found that formula (19) produces very inaccurate results when a large number of parameters m are to be estimated. Clearly, when m > n the problem has no unique solution and additional constraints, e.g. smoothing, should be imposed. To obtain a feasible and a well approximate stable solution, the numerical method proposed in this study is based on minimising the functional
s hIa”
Jl(f) = Ilf’ll’ =
Cf’W12dx =
j$2 FjI:-‘--t;
(20)
Phi”
subject to the bounds of the variables (16) and to the linear constraints JKf-dl
QE,
(21)
where E 2 0 is a positive small quantity to be prescribed. When solving the system of equations (15) with exact data, E is taken to be zero, whilst when solving the noisy data approximation system of equations (17), E is taken to be max ;= 1 l&ii.The actual numerical method modifies the first-order regularisation procedure when the amount of noise added to the exact data is known in advance by eliminating the difficulty in choosing the regularisation parameter. The minimisation of the functional Jl(f) given by expression (20), subject to the simple bounds of the variables (16) and to the linear constraints (21), is solved using the NAG routine E04UCF (see Gill et al., 1986). This routine is designed to minimise an arbitrary smooth function subject to constraints, which may include simple bounds on the variables, linear constraints and smooth non-linear constraints. 4.
NUMERICAL
RESULTS
AND
DISCUSSION
The implementation of the numerical method presented in the previous section was performed for a wide variety of superpositions of lognormal and normal Gaussian distributions, typically of the form
/(r)=A[~lh,exp(-(r~~)2) +~lh;exp(-l$~~))],
(22)
where hk and hi were the amplitudes and the constant A was chosen such that the normalised condition (13) is satisfied. Whilst in atmosphere aerosols the distribution of particle size may be approximated by a lognormal distribution it is the intention of the study to develop a method with a wider range of applicability. For instance, that found in an industrial complex where the aerosols arise from a variety of processes and each process is likely to produce a specific particle size distribution. Also, note that in expression (22), a factor of l/r was not included in the definition of the lognormal distribution in order to obtain significant superposed modes. If the factor l/r is included then the lognormal distributions will overshadow the normal distributions and, for example, the case p = p’ = 1 will produce almost the same distribution as the case p = 0, p’ = 1. In other words, instead of a two peak distribution we would be dealing with what is basically a one peak distribution and hence a much simpler case to investigate. Therefore, expression (22) generates more complicated distributions to be retrieved by the inverse analysis. The magnitude of the function f(r), given by equation (22), is according to the condition (4) very large, of 0(104), whilst the definition domain [l,in,l,,,J = [5 nm, 5OOnm] is very small. Consequently, it is to be expected that the data will be highly sensitive to round-off errors and therefore, for high accuracy, double precision was implemented in all the
D. Lesnic et al.
802
numerical calculations performed. A resealing has also been performed on equation (1) but it was found that it produced virtually the same solution to that present in the absence of any scaling. In this study, the discussion is restricted to a maximum of three superimposed modes, i.e. p + p’ d 3, but the extension of the work to a larger number of modes is at present under further investigation. In order to test the numerical method described in Section 3, using expression (22), the set of the discrete deposition data, d, calculated using equation (15) is taken as the input data and the goal is to retrieve the now unknown size distribution function, f(r). The next step in a future work would be to apply the inversion method for real experimental data. The number of values m at which the continuous function is to be retrieved is taken to be m = 100, which was found to be sufficiently large so as to be able to model any distribution function with a maximum of three peaks. Figures la and lb show the numerical results for a single (i.e. p + p’ = 1) normal mode (pl = 240 nm, g1 = 30 nm) and a lognormal mode (pi = 140 nm, a; = 1.5) distributions, respectively, which are in good agreement with the corresponding analytical solutions. It should be noted that other types of functionals to be minimized have been investigated, such as the following: the zero-order regularisation,
Jo(f) = -f f3
(23)
j= 1
the second-order
regularisation, J2V)
=
f
cfi
-
2&l
(241
+.&212
j=3
or the entropy, E(f) = - f sjln(fj), j=
(25)
1
and it was found that the functional Jl(f), given by expression (20) provides the closest agreement with the analytical solutions. Figures 1 also include the numerical results obtained by using the linear constraints (21) when 5% noise is present and this illustrates the stability of the method for a single-mode distribution under small perturbations in the input data. The noise data d(‘) = d + E, where E = (EC)for i = 1, . . . , n, are included by adding to the exact values d, random variables Ei with mean zero and standard deviation taken to be 5% from the maximum value of di, for i = 1, . . . ,n, generated by the NAG routine GOSDDF (see Brent, 1974). This routine returns a random real number taken from a normal distribution with mean p and standard deviation cr.The kernel functions used here are discussed in a more general form below. When considering double-mode distributions, i.e. p + p’ = 2, the method was found not to give the same good results as those obtained for the single-mode distributions, especially in predicting the amplitude of the second peak. This is probably because, when dealing with a more complex distribution which presents more than one minimum or maximum, the number of linear constraints, namely n + 1 = 10, imposed in expression (21), is too small to retrieve completely all the m = 100 unknowns. More information can be supplied if one extends the discrete integral equations (l), for i = 2, . . . ,8, to a continuous Fredholm integral equation of the first kind, namely IllI.. d(s) =
K(s,r)f(r)dr,
s lmin
sE
C&81,
(26)
where the kernel K(s, r) is defined by K(s,r)
With this prolongation
= Pl(r)2s-2 - Pl(r)2s-1,
s E [2,8],
r E [rmin,rmax].
(27)
the kernel K(s, r) is infinitely differentiable in both variables on the
Determination of the particle size distribution
803
WlO’ 14-
12lo-
ilz +
642-
(4
00 0
0
100
100
200
300 r
Iin nm)
200 r
(in nml
300
400
500
400
500
Fig. 1. The variation of f(r) as a function of r obtained using constraints (21) for the exact data -), 5% noisy data (- - -) and the analytical solution (F), for a single-mode distribution (a) a normal mode (pLI= 240 nm, CJ~= 30 nm); (b) a lognormal mode (p’, = 140 nm, CT’=, 1.5).
(-
interval [2,8] x [rmin,r,,J
and it is easy to observe that
K(i, r) = Ki(r)y
i =
2, . . . 787 r E Crmin, ~maxl
(28)
since from equations (7b) and (8) we have K,(r)
= PI(r)“-’
- P,(r)“i
= Pl(r)2’-’
- Pl(r)Zi-‘,
i =
2, . . . ,8.
(29)
Initially, we suppose that the integral equation (26) can be applied at many points s within the real interval [2,8]. Consider for this an equidistant subdivision [sk_ 1,sk], for
D. Lesnic et al.
804
k = 1, . . . , n,, of the interval
[2,8],
given by
sk = 2 + 6k/nl,
k = 0, . . . ,n,
(30)
and taking n, = 6n2,
(31)
we obtain sk = 2 + k/nz,
k = 0 ,..., nl.
(32)
The subdivision (32) was chosen such that the interval subdivision points, the integer values of {2,3, . . . ,7,8> for some values of k, namely for
k E {0,n*,2n2,3n2,4n2,5n2,6nz). Applying
now the integral
equation
(26) at the points
s
Sk, coincide
with
(33)
s = Sk, for k = 0, . . . ,n,, we obtain
rm.i
4SJ = which in discretised
k = 0, . . . ,nl,
K(sk,r)f(r)dr, r”un
form results
(34)
in
d(sk)=
f j=
Kkjfj,
k
=
(35)
O,...,nl,
1
where 5
Kkj
=
K(sk,r)dr,
k = 0, . . . . n,,
j = 1, .
m.
(36)
Equations (11) applied for i = 1 and i = 9, and equations (13) and (35) form a system of (nl + 1) + 1 + 2 = n, + 4 linear equations with m unknowns and the numerical method for determining the unknown vector f applies as described in Section 3. Recall now that from equation (31) n1 should be a multiple of 6 and, in this study a value of n, = 120 was considered sufficiently large such that the number of equations, n, + 4 = 124, exceeded the number of unknowns, m = 100. A further increase in the value of n, may cause incompatibility problems related to the rounding-off errors. Also, we remark that an increase in the number of stages of the diffusion battery, n = 9, is not to be recommended since from equation (29) as PI(r) E [0, l), the values of K;(r) will be virtually zero for i > 9 and will make the system of equations (15) even more singular. In addition, with increasing the number of stages the parasitic particle deposition should be taken into consideration and the analysis becomes more complicated. The matrix K of n + 1 = 10 rows and m = 100 columns defined by equation (14) has the rank 9 and this figure does not increase as the number of stages of the diffusion battery is increased (see Lloyd, 1994), whilst the extended matrix K of n, + 4 = 124 rows and m = 100 columns defined by equation (36) has the rank 19. Therefore, it is to be expected that the extension of the kernel given by equation (27) will improve the estimation of the function to be retrieved. At this stage, it should be noted that equations (34) cannot be used exactly in the absence of the knowledge of function f(r) which is unknown. In fact, only 9 data values are available, passing through the namely di, for i = 1, . . . ,9. However, a piecewise linear interpolation data points di, for i = 2, . . . (8, can be considered. Although some undesirable noise is now added to the data, this penalty was found not to affect significantly the numerical results in comparison with the improvement achieved in the prediction of the analytical solution. Figure 2 shows the piecewise linear approximation passing through 9 exact data values di, for i = 1, . . . ,9, for a double-mode normal distribution (pl = 40 nm, (rl = 10 nm and pL2= 200 nm, g2 = 15 nm) in comparison with the exact values d(s) given by equation (26). Based on this approximation, Fig. 3a shows the variation of f(r) as a function of r obtained by minimising the functional Jl(f) given by expression (20) subject to the bounds of the
Determination
of the particle size distribution
805
0.15-
z n
O.l-
0.05-
i
sb
sb
4b
2b
id0
ii0
S ), the deposition measurements Fig. 2. The deposition measurements values di, for i = 1, .. . ,9, ( ?? function d(s) (--) and its piecewise linear approximation (---) for a double-mode normal distribution (rI = 40 nm, u1 = 10 nm and p2 = 200 nm, o2 = 15 nm).
variables (16) and to the n, + 4 = 124 linear constraints d(SJdi -
f
j=l f
Kkjfj
Kijfj
nl,
(374
i = 1 and i = 9,
WW
G-C,
GE,
k=O ,...,
j=l
1 =
f j=l
&n+l,jfj-
r (in nm) Fig. 3a.
(37c)
D. Lesnic et al.
806
200
300
400
500
400
500
r (in nm)
x10’
6
200
300 r (in nm)
Fig. 3. The variation of f(r) as a function of I obtained using constraints (21) (-- -), using constraints (37) (-+) and the analytical solution (--) for double-mode distributions (a) normal modes (pi = 40 nm, uI = 10 nm and p2 = 200 nm, oz = 15 nm); (b) lognormal modes (pi = 40 nm, u’, = 1.3 and p; = 140 nm, CT;= 1.5); (c) a combined normal mode (pi = 40 nm, u1 = 10 nm) and a lognormal mode (p’, = 140 nm, g’, = 1.5).
This is shown to improve the estimation of the analytical solution in comparison with the numerical results obtained by using only n + 1 = 10 constraints given by expressions (21). Also, the same conclusion may be drawn for a double-mode lognormal OJL;= 40 nm, a; = 1.3 and p; = 140 nm, 05 = 1.5) and both a normal and lognormal mode (pi = 40 nm, e1 = 10 nm and p’r = 140 nm, 0; = 1.5) illustrated in Figs 3b and c, respectively. At this stage it should be noted that the numerical method employed in this study is able to estimate well single- and double-mode distributions which have significantly positive
Determination
of the particle
size distribution
807
values only in the approximate range of [5 nm, 300nm]. The remaining range of particle sizes [I300nm, 500 nm] is associated with the region where, for the actual diffusion battery characteristic design and volumetric flow rate, the fractional penetration functions defined by equation (8) become non-monotonic (see also Wu et al., 1989), and therefore it introduces non-uniqueness features. For other research which correlates the volumetric flow rate to a given particle size range of interest on which the fractional penetration functions become monotonic, see Lesnic et al. (1995). The piecewise linear approximation plotted in Fig. 4 for a triple mode, i.e. p + p’ = 3, normal distribution (,u~ = 40nm, crl = lOnm, pz = 2OOnm, c2 = 15nm and p3 = 360nm, u3 = 20 nm) shows the same characteristics to those presented in Fig. 2 for a double-mode normal distribution. With this approximation the variation of f(r) as a function of r using the new method based on the constraints (37) is shown in Fig. 5a. Although these results are still poor in estimating the analytical solution for the last two peaks, they give more information about the properties of the solution than does the original method, which is based on constraints (21). However, for the triple-mode distribution, better results may be achieved if one uses the numerical method presented previously, but now in a sequential manner. Based on the local recovery ofeach peak of the input distribution it produces better estimations than the global recovery of all the peaks. From Fig. 5a it can be seen that the initial full application of the numerical method, based on constraints (37), provides results which are in very good agreement with the analytical solution for the function f(r) where the first peak is located, and this observation will be further assumed. Thus, if fi(r) defines the numerical solution in an interval [l,in, r*] which extends to zero values in the interval (r*,r,,J, where r* is chosen as the position of a local minimum in the first region, e.g. I* E 100 nm for Fig. 5a, then the function f(r) can be expressed as
r E Crmin, r,,,l,
f(r) = fi (r) + f(r),
(38)
where f(r) is an unknown function to be determined in the interval (r*, rmax]and taken to be zero in the interval [rmin,r*]. On using equation (26) the problem for T(r) reformulates as I* Ima, K(s, r)T(r) dr = d(s) = d(s) - c K(s, r) fi (r) dr, s E C2,Sl. (39) J
Jhl”
I*
0.2
0.15
z ‘0 0.1
I 0
I
I
I
I
I
I
20
40
60
80
100
120
S
Fig. 4. The deposition measurements values di, for i = 1, . ,9, (a), the deposition measurements function d(s) (-) and its piecewise linear approximation (---) for a triple-mode normal distribution (pl = 40 nm, o1 = 10 nm, p12 = 200 nm, CJ*= 15 nm and ps = 360 nm, u3 = 20 nm).
D. Lesnic et al.
808
6.0
260
(4
300
400
500
r (in nm)
6
2
(b)
r (in nm) Fig. Sa and b.
Analogous constraints for f(r) can be deduced from the integral equation (39) as the constraints (37) for the initial function f(r) have been deduced from the integral equation (26). A new application of the method in the interval (r*,r,,,] produced results which estimate well the second region of the function f(r) where the second peak is located. Similarly, as above we can define by f*(r) the numerical solution in the interval @*,I**] which extends to zero values outside of this interval, where r ** is chosen as the position of a local minimum for the second region, e.g. r ** z 250 nm for Fig. 5a. Then the function f(r) can be expressed as f(r) = _Mr) + f’(r),
rE
Crmin, r,,,l~
(40)
Determination of the particle size distribution
200
300 r
O-J
(4
’
809
(in nm)
1
I
I
I
100
200
300
400
1 500
r (in nm) Fig. 5. The variation of f(r) as a function of r obtained using constraints (21) (- --), using constraints (37) (---), the successive algorithm based on the integral equations (26), (39) and (41) (-Cl-) and the analytical solution (---) for (a) three normal modes (pr = 40 nm, o1 = 10 nm, p2 = 200 nm, e2 = 15 nm and p3 = 360 nm, o3 = 20 nm); (b) combined two normal modes (pl = 200 nm, et = 15 nm and JQ = 360 nm, e2 = 20 nm) and one lognormal mode (& = 40 nm, CT;= 1.3); (c) combined two lognormal modes (& = 40 nm, CT;= 1.3 and & = 140 nm, u; = 1.5) and one normal mode (cc1= 200 nm, o1 = 15 nm); (d) three lognormal modes (,u‘, = 40 nm, 6; = 1.3, & = 140 nm, 6; = 1.5 and pi = 240 nm, u; = 1.7).
where the unknown satisfies
function
f(r)
to be determined
in the remaining
r** In?*. K(s, r)f(r) dr = a(s) = a(s) K (s, r) Sz (r) dr, s I. s I** Unfortunately,
a final application
of the method
interval
81.
(41)
equation
(41) in the
s E D,
based on the integral
(r**,r,,,]
810
D. Lesnic et al.
6 z .+ 4
0
100
200
300 r
(in nm)
400
500
Fig. 6. The variation of f(r) as a function of r obtained using constraints (21) (- --), using constraints (37) (-v), the successive algorithm based on the integral equations (26), (39) and (41) (-U-) and the analytical solution (--) for a three normal mode, p1 = 40 nm, o1 = 10 nm, hI = 1, pz = 200 nm, CJ~= 15 nm, hz = 1.5 and bcj = 360 nm, o3 = 20 nm, h, = 2.
interval (r** , r,,,] did not produce a sufficiently well desired estimation-for the third peak. As aforementioned the difficulty in numerically estimating the function f(r) is a result of the fact that the interval (r** ,rmax] = [250 nm, 500 nm] in which the function is sought contains the “problem” range in which the fractional penetration functions become non-monotonic. Therefore, it is to be expected that the solution f=(r) will be a poor approximation since there is insufficient information provided by the integral equation (41) in the interval (r**, I,,, 1. In addition, a large amount of noise might have been accumulated through the successive algorithms procedure. The numerical results obtained by successive application of the initial method based on the integral equations (26), (39) and (41) are included in Fig. 5a, and show a good estimation of the analytical solution at least in the first two regions where the first two peaks are located. The same conclusion may be drawn from Figs 5b and c for the combined triple mode, two normal and one lognormal (pl = 200 nm, c1 = 15 nm, p2 = 360 nm, c2 = 20 nm and p; = 40 nm, a; = 1.3) and the two lognormal and one normal (pl = 200nm, g1 = 15 nm, PL; = 40nm, 0; = 1.3 and p; = 140nm, o$ = 1.5) distributions, respectively. Finally, for completeness, Fig. 5d shows the results for a triple-mode lognormal distribution (pL; = 40 nm, a; = 1.3, p> = 140 nm, a; = 1.5 and p; = 240 nm, o; = 1.7) for which the sequential algorithm applied for other types of triple-mode distributions was unnecessary since in this case the function has actually only two peaks similar to that for a double-mode distribution, In all the calculations undertaken so far, hk, for k = 1, . . . ,p, and hi, for k = 1, . . . , p’, have either been set to zero or unity. However, removing this equal amplitude restriction for any non-zero modes produces, as shown in Fig. 6, where h1 = 1, h2 = 1.5 and h3 = 2, the same success for the sequential method as in Figs 5a-c. From the overall discussion of the results analysed in this section it can be concluded that, for the diffusion battery characteristic design investigated the numerical method introduced in this study estimates well the first two regions, i.e. where the first two peaks are located within the approximate range of [5 nm, 300 nm] for any fine particle size distribution.
Determination of the particle size distribution
811
5. CONCLUSIONS
The retrieval of the unknown aerosol particle size distribution from a set of discrete diffusion battery measurements has been analysed. The numerical method introduced in this study does not make use of any a priori information about the form of the solution and is based on a modified first-order-regularisation procedure by minimising the norm of the first derivative subject to bounds of the variables and linear constraints. In addition, it eliminates the difficulty of choosing the regularisation parameter when the amount of noise included in the data is known in advance. For single-mode distributions (normal or lognormal) having significant positive values only in the range [5 nm, 300 nm] the numerical results agree well with the analytical solutions. When double-mode distributions are present the problem is approximated as a Fredholm integral equation of the first kind for which the numerical solution provided results which are in good agreement with the analytical solutions. Finally, the successive application of the numerical method enabled good results for the first two peaks in a triple-mode distribution to be obtained. A further study, which is at present under investigation, will be concerned with the relations between the diffusion battery characteristic design, the volumetric flow rate, the particle size range of interest and the interval of monotonicity of the fractional penetration functions of aerosols at the stages of the diffusion battery. Acknowledgements-We would like to thank Dr J. J. Lloyd from Manchester Royal Infirmary for all his kind assistance during the preparation of this study and for providing the values of the constants necessary for defining the set of deposition functions in Section 2. Also the very constructive comments of the referees are gratefully acknowledged.
REFERENCES Baklanov, A. M. and Dubtsov, S. N. (1993) An experimental set-up for aerosol spectometers calibration. J. Aerosol Sci. 24, S237-S238. Bashurova, V. S., Koutzenogii, K. P., Pusep, A. Y. and Shokhirev, N. V. (1991) Determination of atmospheric aerosol size distribution functions from screen diffusion battery data: mathematical aspects. J. Aerosol Sci. 22, 373-388. Brent, R. P. (1974), Algorithm 488. Comm. ACM, 704-706. Cheng, Y. S., Keating, J. A. and Kanapilly, G. M. (1980) Theory and calibration of a screen-type diffusion battery, J. Aerosol Sci. 11, 549-556.
Cheng, Y. S. and Yeh, H. C. (1980) Theory of a screen-type diffusion battery. J. Aerosol Sci. 11, 313-320. Cheng, Y. S. and Yeh, H. C. (1983) Performance of a screen-type diffusion battery. In Aerosols in Mining and Industrial Work Enuironments (Edited by Marple, V. A. and Liu, B. Y. H.), pp. 1077-1094. Wiley, New York. Cheng, Y. S., Yeh, H. C. and Brinsko, K. J. (1985) Use of wire screens as a fan model filter. Aerosol Sci. Technol. 4, 165-174.
Cooper, D. W. and Spielman, L. A. (1976) Data inversion using nonlinear programming with physical constraints: aerosol size distribution measurement by impactors. Atmos. Enuir. 10, 723-729. Cooper, D. W. and Wu, J. J. (1990) The inversion matrix and error estimation in data inversion: application to diffusion battery measurements. .I. Aerosol Sci. 21, 217-226. Emi, H., Wang, C. S. and Tien, C. (1982) Transient behaviour of aerosol filtration in model filters. AIChE J. 28, 397-405.
Gill, P. E., Hammarling, S. J., Murray, W., Saunders, M. A. and Wright, M. H. (1986) User’s Guide for LSSOL (Version 1.0). Report SOL 86-1, Stanford University. Heintzenberg, J. (1975) Determination in situ of the size distribution of the atmospheric aerosol. J. Aerosol Sci. 6, 291-303.
Kirsch, A. A. and Stechkina, I. B. (1978) The theory of aerosol filtration with fibrous filters. In Fundamentals of Aerosol Science (Edited by Shaw, D. T.), pp. 165-256. Wiley, New York. Lesnic, D., Elliott, L. and Ingham, D. B. (1995) The inverse determination of the particle size distribution from diffusion battery data. 9th U.K. Ann. Aerosol Sot. Conf., Norwich (submitted). Lloyd, J. J. (1994) Prioate communication. Lloyd, J. J., Taylor, C. J. and Shields, R. A. (1994) Analysis of diffusion battery data using simulated annealing. Pkoc. 8th U.K. Ann. Aerosol Sot. Co& New York, pp. 29-34. Paatero, P. (1990) The extreme value estimation deconvolution method with applications in aerosol research. Univ. Helsinki kept. ‘Ser. Phys., HU-P-250. Sinclair, D. and Hopes, G. S. (1975) A novel form of diffusion battery. Am. Ind. Hyg. Assoc. J. 36, 39-42. Skaptsov, A. S. and Lipatov, G. N. (1994) On the interpretation of rectangular diffusion battery data. J. Aerosol Sci. 25, 341-344.
Soderholm, S. C. (1983) Analyzing screen diffusion battery data. In Aerosols in Mining and Industrial Work Environments (Edited by Marple, V. A. and Liu, B. Y. H.), pp. 1095-1115. Wiley, New York. Tikhonov, A. N. and Arsenin, V. Y. (1977) Solutions of Ill-Posed Problems. W. H. Winston & Sons, Washington.
812
D. Lesnic et al.
Twomey, S. (1963) On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature. J. Assoc. Comput. Mach. 10, 97-101. Twomey, S. (1965) The application of numerical filtering to the solution of integral equations encounted in indirect sensing measurements. .I. Franklin Inst. 279, Y95-109. Twomey, S. (1975) Comparison of constrained linear inversion and an iterative nonlinear algorithm applied to the indirect estimation of particle size distribution. J. Cornput. Plays. 18, 188-200. Wolfenbarger, J. K. and Seinfeld, J. H. (1990) Inversion of aerosol size distribution data. .I. Aerosol Sci. 21, 227-247. Wu, J. J., Cooper, D. W. and Miller, R. J. (1989) Evaluation of aerosol deconvolution algorithms for determining submicron particle size distributions with diffusion battery and condensation nucleus counter. J. Aerosol Sci. Zo, 477-482. Yee, E. (1989) On the interpretation of diffusion battery data. J. Aerosol Sci. 20, 797-811.
J.Aerosol
SCI , Vol. 26, No. 5, pp. 797X312, 1995 Cowrinht 0 1995 Elsevier Science Ltd Printed’,; &at Bntain. All rights reserved CNl21~8502/95 $9.50 + 0.00
0021-8502(95)00007-0
AN INVERSION METHOD FOR THE DETERMINATION OF THE PARTICLE SIZE DISTRIBUTION FROM DIFFUSION BATTERY MEASUREMENTS D. Lesnic, L. Elliott Department
of Applied
Mathematical
Studies,
and D. B. Ingham University
of Leeds, Leeds LS2 9JT, U.K.
(First received 8 September 1994; and in final form 3 January 1995)
Abstract-The inverse problem of the determination of the aerosol particle size distribution from a set of deposition measurements taken on the screens of the stages of a diffusion battery is analysed. The mathematical formulation reduces the problem to solving a system of integral equations for the unknown particle size distribution. Since the problem is ill-posed a numerical method which is based on minimising the norm of the derivative function subject to simple bounds of the variables and certain linear constraints, and which does not use any a priori information regarding the parametrised form of the solution is employed. For single- and double (normal or lognormal)-mode particle size distributions the numerical method provides a good estimation of the corresponding analytical solutions provided that, in general, the problem is approximated as a continuous Fredholm integral equation of the first kind. For triple-mode distributions a sequential application of the numerical method provides results which give good estimates for the first two modes but the prediction of the third mode is only satisfactory.
1. INTRODUCTION
An important problem in atmospheric aerosol physics relates the size distribution of fine aerosol particles (diameter less than 1 pm) to a set of measurements obtained using a particle sizing instrument. These instruments include diffusion batteries, which consist of a set of rectangular or circular channels (channel-type), or fibrous filters (screen-type) through which an aerosol is passed. At this stage it is worth mentioning that the analysis performed in this paper can in principle be applied to other particle size instruments such as cascade impactors, variable slit impactors, optical counters, electrical mobility analyser, etc., and will be considered in future investigations. Although in all these instruments there is only a limited number of data points which are available to obtain a particle size distribution, the situation appears more difficult for the diffusion battery case since the diffusion battery responses, i.e. the deposition kernels, are, in general, are not sharp and may be non-monotonic (see Wu et al., 1989). It is this unusual behaviour of the diffusion battery kernels compared with those in other spheres of interest that necessitates the more elaborate inversion procedures discussed herein to retrieve the particle size distribution function. If an inversion method developed for a non-monotonic response instrument is shown to perform well, then for a monotonic response instrument the method will be expected to perform even better. In addition, for the diffusion battery, the existence of an analytical form for the deposition kernels enables a more accurate and simple development of the inversion method, whilst for most other instruments the deposition kernels are obtained empirically, hence requiring much interpolation, resulting in the introduction of additional noise into the problem. Some of the previous studies were concerned with the design and calibration of diffusion batteries (see Sinclair and Hopes, 1975; Cheng et al., 1980; Baklanov and Dubtsov, 1993). Also, many other studies, both experimental and theoretical, have dealt with the determination of the deposition functions based on the analysis of particle collection efficiency in channels or filters (see, for example, Emi et al., 1982; Cheng and Yeh, 1983; Cheng et al., 1985). When the particle size distribution is known in advance then one can evaluate directly the deposition data (see Skaptsov and Lipatov, 1994). However, in practice, only the measurements of the fraction of particles collected at each state through the diffusion 797 AS 26-5-G
D. Lesnic et al.
798
battery are available, and this gives rise to a typical inverse problem for the determination of the unknown particle size distribution, and which still remains an unsolved problem (see Yee, 1989). Given a set of fractional measurements, di, for i = 1, . . . , n, and a set of deposition functions, Ki(r), for i = 1, . . . , n, it is desirable to retrieve, as far as possible, the unknown continuous particle size distribution function, f(r), which satisfies the integral equations di =
Im*= Ki(r)f(r)dr, s hi”
i = 1, . . . , n,
(1)
where [rmin, rmax] is the particle size (radius) range of interest. It is clear that the solution of the system of integral equations (1) is not unique, since there are an infinite number of continuous size distribution functions f(r) which will yield the discrete set of deposition measurements, di, for i = 1, . . . , n. In addition, since there are errors associated with each measurement, equations (1) are satisfied only within a tolerance Ei, for i = 1, . . . , n, namely lms. d!“’= d.I + s.I = Ki(r)f@) (r) dr > i=l 3*.. 9n. (2) s hi” The small quantities si, for i = 1, . . . ,n, are introduced in order to simulate the inherent experimental measurement errors that may cause large deviations of the solution f’“‘(r) of the system of equations (2) from its exact solution f(r) of the system of equations (1) (see Wolfenbarger and Seinfeld, 1990), i.e. there is a possibility of a discontinuous dependence of the solution upon the input data. The non-uniqueness and any discontinuous dependence upon the data show that the inverse problem (1) is ill-posed. However, most of the many solutions of problem (1) are physically unacceptable because they are highly oscillatory and include negative values. If a numerical procedure produces negative values for the solution, then simply replacing them with zero values, as suggested by Bashurova et al. (1991), might not be the best option in certain cases. In order to overcome these difficulties, additional physical constraints may be added to the set of equations (1) or (2) in order to restrict the class of possible solutions. In particular, we enforce that the particle size distribution be a positive function, i.e.
rmaxl~ r E Crdn,
f(r) 2 0,
(3)
and also a normalised condition is required which, without loss of generality, may be written as (4) Nevertheless, the natural constraints (3) and (4) are still too weak to enable problem (1) to have a unique solution. Various numerical methods have been proposed for the solution of the problem (l)-(4). These methods include: manual techniques based on trial function analysis (see Soderholm, 1983), regularisation procedures (see Twomey, 1963), numerical filtering (see Twomey, 1965), nonlinear iterative methods (see Twomey, 1975), nonlinear programming using physical constraints (see Cooper and Spielman, 1976), maximal entropy production (see Yee, 1989), extreme value estimation (see Paatero, 1990), etc. Most of these studies were concerned in dealing with the discontinuous dependence upon the data, i.e. stabilising the results by imposing smoothing constraints. However, when dealing with the response of a diffusion battery the highly ill-posed nature of the mathematical problem is mainly due to the non-unicity of the solution rather than to the discontinuous dependence on the data for which a regularisation procedure will offer, in general, a stable solution (see Tikhonov and Arsenin, 1977). A priori parameter estimations of the form of solution have been proposed by assuming that the particle size distribution is a superposition of lognormal kernels (see Heintzenberg, 1975): fcr)
=
i
k=l
hk
exp
r
Determination of the particle size distribution
799
or normal kernels (see Lloyd et al., 1994); f(r) = 2 hk exp -@ - I&)* ( XT-> k=l
(6)
where hkrpk and a&.,for k = 1, . . . ,p, represent the amplitude, the mean and the standard deviation, respectively. In this way the infinite-dimensional space in which the solutions may lie is reduced to a finite-dimensional space but pays the penalty that the solution is a priori assumed analytical, which implies a degree of smoothness which is not satisfied by the real particle size distributions. In addition, other types of possible particle size distributions such as Junge type or Gamma distributions will never be recovered. It is the purpose of this study to introduce a numerical method which does not use any a priori information regarding the form of the solution, i.e. the parameter estimation given by equations (5) or (6), and which is based only on a modified first-order regularisation approach in order to achieve the stability of the results. Also, the degree of non-uniqueness of the problem may be reduced by increasing the number of equations under a piecewise linear approximation of the deposition data and a suitable extension of the deposition functions. This can easily be performed by joining the discrete experimental data by straight lines. Furthermore, a successive application of the numerical method enables good estimations of the available analytical solutions, up to the first two peaks in a triple-mode distribution, to be obtained.
2. DIFFUSION
BATTERY DESIGN
AND DEPOSITION
FUNCTIONS
The diffusion battery design is similar to that analysed by Lloyd et al. (1994), i.e. it consists of eight stages of stacks of 400 mesh woven wire screens, with a ninth stage which is an absolute filter, i.e. there are no particles passing through it. Particles with a radius smaller than r,,,i,, are all deposited on the first screen, whilst particles with a radius greater than rmaxare removed by an impactor at the diffusion battery inlet. In this section the input data which is necessary for the inverse analysis developed in Section 3 is defined. The deposition function, Ki(r), as a function of particle radius, r, on a given stage i, for i=l , . . . ,9, is the difference between the penetration through all stages including itself, Pi, and the penetration up to that stage, Pi-l, namely K,(r) = 1 - PI(r), Ki(r) = Pi_ l(r) - Pi(r), K9(r1
=
(74
i = 2, . . . ,8,
(7b)
(7c)
P8W
The fractional penetration of aerosols at stage i, Pi(r), according to the fan model filtration theory of the wire screen diffusion battery (see Cheng and Yeh, 1980) is expressible as Pi(r) = exp
-
2niah q(r) , 7c(l - a)u >
i=l,...,
8,
where u = 0.3116 is the solid volume fraction of the screen, ni = 2’-’ is the cummulative number of screens from the first stage down to stage i, for i = 1, . . . ,8, h = 76.9 pm is the screen thickness and a = 14.7 pm is the wire radius. The single wire collection efficiency, q(r), can be expressed as the sum of the individual efficiencies due to diffusion, direct interception, inertial impaction and a correction term due to diffusion and interception (see Kirsch and Stechkina, 1978), namely v(r) = v&) + vR(r) + 11(r) + IfDR(r),
(9)
where qD(r) = 2.7 Pe(r)-*13, qR(r) = (2x)-‘[(l
+ R(r))-’ - 1 - R(r) + 2(1 + R(r))ln(l + R(r))],
Wa) (lob)
D. Lesnic et al.
800
q,(r) = (2x)-’ Stk(r) C(29.6 - 28cc0.62)R(r)2 - 27.5R(r)‘.*], qDR(r) = 1.24~~ ‘I2 Pe(r)- l/2 R(r)213,
(1W Wd)
where Pe(r) = 2aU,/D(r) is the Peclet number, U o = 4Q/nd: is the freestream fluid veloflow rate, dr = 3.5 cm is the screen diameter, city, Q = 0.44 lmin- r is the volumetric diffusion coefficient for a spherical particle, D(r) = kTC(r)/6qw is the k = 1.38054x lo-r6ergK-’ is the Boltzmann constant, T = 293 K is the absolute temperature, C(r) = 1 + (1.246 + 0.42 exp (- 0.87 r/A)) A/r is Cunningham slip correction factor, A = 0.0665 pm is the mean free path of the gas molecules, ,U = 0.000184 P is the fluid viscosity, R(r) = r/a is the interception parameter, x = -0.5 In(c) - 0.75 + c - 0.25c2, c = 2c+r, Stk(r) = ppUOr2C(r)/(9pa) is the Stokes number and pp is the particle density. Since only fine aerosol particles are investigated the particle radius range of interest is taken as [r,in, rmax] = [5 nm, 500 nm] and therefore the inertial term, q,(r), in expression (9) can be neglected.
3. NUMERICAL
METHOD
Discretising the interval [rmin, r,,,] into m equidistant intervals (i.e. elements), [rj_ r, rj], for j = 1, . . . , WI,and assuming the function f(r) to be constant over each element and taking the value at its midpoint, r; (i.e. nodes), equations (l), (3) and (4) become di = 2 j=l
where fj =f(r$)
Kijfj,
i = l,...,n,
(11)
fj20,
j=l,...,
m,
(12)
j$l(rj
- rj-l)fj
= 1,
(13)
Ki(r)dr,
i=
and the matrix IJ Kij
=
s
l,...,
n, j=
l,...,
m,
(14)
rj-I
is evaluated using a trapezoidal quadrature rule. Better quadrature rules, e.g. Simpson rule, can be used but it was found not to have any significant effect on the final results. Extending the matrix K by defining Kc,+ l)j = rj - rj_ 1, for j = 1, . . , m, and the vector d by defining d n+1 = 1, equations (11) and (13) result in a linear system of (n + 1) equations with m unknowns f = (fj), for j = 1, . . . , m, namely d=Kf which has to be solved subject
to the positive
(15) constraint
condition
(12), namely
f>O. Similarly, form
when noise is present
in equations
(2) the resulting
(16) system of equations
d(E) = Kf’“’ 9 where d = (dy)) for i = l,...,(n If m = II + 1 and the matrix
+ 1) and dfi, = 1. K is non-singular, then one would f=K-‘d.
takes the
(17) simply have (18)
However, because the problem is ill-posed the system of equations (15) is highly illconditioned (see Cooper and Wu, 1990), and therefore the direct method based on the inversion of equation (15) cannot be applied. If m < n + 1 = 10, the approximation of the integral equations (1) with the discrete version (11) would be nevertheless unrealistic for functions which rapidly change their
Determination
of the particle
size distribution
801
magnitude. In addition, it is desirable to obtain a continuous function rather than at a few discrete points and therefore m should be large, much greater than n + 1. Alternatively, one may still utilise the least-squares method based on minimising the norm J/d - KflJ’, which has the solution f = (K”K)-‘Pd.
(19)
However, it was found that formula (19) produces very inaccurate results when a large number of parameters m are to be estimated. Clearly, when m > n the problem has no unique solution and additional constraints, e.g. smoothing, should be imposed. To obtain a feasible and a well approximate stable solution, the numerical method proposed in this study is based on minimising the functional
s hIa”
Jl(f) = Ilf’ll’ =
Cf’W12dx =
j$2 FjI:-‘--t;
(20)
Phi”
subject to the bounds of the variables (16) and to the linear constraints JKf-dl
QE,
(21)
where E 2 0 is a positive small quantity to be prescribed. When solving the system of equations (15) with exact data, E is taken to be zero, whilst when solving the noisy data approximation system of equations (17), E is taken to be max ;= 1 l&ii.The actual numerical method modifies the first-order regularisation procedure when the amount of noise added to the exact data is known in advance by eliminating the difficulty in choosing the regularisation parameter. The minimisation of the functional Jl(f) given by expression (20), subject to the simple bounds of the variables (16) and to the linear constraints (21), is solved using the NAG routine E04UCF (see Gill et al., 1986). This routine is designed to minimise an arbitrary smooth function subject to constraints, which may include simple bounds on the variables, linear constraints and smooth non-linear constraints. 4.
NUMERICAL
RESULTS
AND
DISCUSSION
The implementation of the numerical method presented in the previous section was performed for a wide variety of superpositions of lognormal and normal Gaussian distributions, typically of the form
/(r)=A[~lh,exp(-(r~~)2) +~lh;exp(-l$~~))],
(22)
where hk and hi were the amplitudes and the constant A was chosen such that the normalised condition (13) is satisfied. Whilst in atmosphere aerosols the distribution of particle size may be approximated by a lognormal distribution it is the intention of the study to develop a method with a wider range of applicability. For instance, that found in an industrial complex where the aerosols arise from a variety of processes and each process is likely to produce a specific particle size distribution. Also, note that in expression (22), a factor of l/r was not included in the definition of the lognormal distribution in order to obtain significant superposed modes. If the factor l/r is included then the lognormal distributions will overshadow the normal distributions and, for example, the case p = p’ = 1 will produce almost the same distribution as the case p = 0, p’ = 1. In other words, instead of a two peak distribution we would be dealing with what is basically a one peak distribution and hence a much simpler case to investigate. Therefore, expression (22) generates more complicated distributions to be retrieved by the inverse analysis. The magnitude of the function f(r), given by equation (22), is according to the condition (4) very large, of 0(104), whilst the definition domain [l,in,l,,,J = [5 nm, 5OOnm] is very small. Consequently, it is to be expected that the data will be highly sensitive to round-off errors and therefore, for high accuracy, double precision was implemented in all the
D. Lesnic et al.
802
numerical calculations performed. A resealing has also been performed on equation (1) but it was found that it produced virtually the same solution to that present in the absence of any scaling. In this study, the discussion is restricted to a maximum of three superimposed modes, i.e. p + p’ d 3, but the extension of the work to a larger number of modes is at present under further investigation. In order to test the numerical method described in Section 3, using expression (22), the set of the discrete deposition data, d, calculated using equation (15) is taken as the input data and the goal is to retrieve the now unknown size distribution function, f(r). The next step in a future work would be to apply the inversion method for real experimental data. The number of values m at which the continuous function is to be retrieved is taken to be m = 100, which was found to be sufficiently large so as to be able to model any distribution function with a maximum of three peaks. Figures la and lb show the numerical results for a single (i.e. p + p’ = 1) normal mode (pl = 240 nm, g1 = 30 nm) and a lognormal mode (pi = 140 nm, a; = 1.5) distributions, respectively, which are in good agreement with the corresponding analytical solutions. It should be noted that other types of functionals to be minimized have been investigated, such as the following: the zero-order regularisation,
Jo(f) = -f f3
(23)
j= 1
the second-order
regularisation, J2V)
=
f
cfi
-
2&l
(241
+.&212
j=3
or the entropy, E(f) = - f sjln(fj), j=
(25)
1
and it was found that the functional Jl(f), given by expression (20) provides the closest agreement with the analytical solutions. Figures 1 also include the numerical results obtained by using the linear constraints (21) when 5% noise is present and this illustrates the stability of the method for a single-mode distribution under small perturbations in the input data. The noise data d(‘) = d + E, where E = (EC)for i = 1, . . . , n, are included by adding to the exact values d, random variables Ei with mean zero and standard deviation taken to be 5% from the maximum value of di, for i = 1, . . . ,n, generated by the NAG routine GOSDDF (see Brent, 1974). This routine returns a random real number taken from a normal distribution with mean p and standard deviation cr.The kernel functions used here are discussed in a more general form below. When considering double-mode distributions, i.e. p + p’ = 2, the method was found not to give the same good results as those obtained for the single-mode distributions, especially in predicting the amplitude of the second peak. This is probably because, when dealing with a more complex distribution which presents more than one minimum or maximum, the number of linear constraints, namely n + 1 = 10, imposed in expression (21), is too small to retrieve completely all the m = 100 unknowns. More information can be supplied if one extends the discrete integral equations (l), for i = 2, . . . ,8, to a continuous Fredholm integral equation of the first kind, namely IllI.. d(s) =
K(s,r)f(r)dr,
s lmin
sE
C&81,
(26)
where the kernel K(s, r) is defined by K(s,r)
With this prolongation
= Pl(r)2s-2 - Pl(r)2s-1,
s E [2,8],
r E [rmin,rmax].
(27)
the kernel K(s, r) is infinitely differentiable in both variables on the
Determination of the particle size distribution
803
WlO’ 14-
12lo-
ilz +
642-
(4
00 0
0
100
100
200
300 r
Iin nm)
200 r
(in nml
300
400
500
400
500
Fig. 1. The variation of f(r) as a function of r obtained using constraints (21) for the exact data -), 5% noisy data (- - -) and the analytical solution (F), for a single-mode distribution (a) a normal mode (pLI= 240 nm, CJ~= 30 nm); (b) a lognormal mode (p’, = 140 nm, CT’=, 1.5).
(-
interval [2,8] x [rmin,r,,J
and it is easy to observe that
K(i, r) = Ki(r)y
i =
2, . . . 787 r E Crmin, ~maxl
(28)
since from equations (7b) and (8) we have K,(r)
= PI(r)“-’
- P,(r)“i
= Pl(r)2’-’
- Pl(r)Zi-‘,
i =
2, . . . ,8.
(29)
Initially, we suppose that the integral equation (26) can be applied at many points s within the real interval [2,8]. Consider for this an equidistant subdivision [sk_ 1,sk], for
D. Lesnic et al.
804
k = 1, . . . , n,, of the interval
[2,8],
given by
sk = 2 + 6k/nl,
k = 0, . . . ,n,
(30)
and taking n, = 6n2,
(31)
we obtain sk = 2 + k/nz,
k = 0 ,..., nl.
(32)
The subdivision (32) was chosen such that the interval subdivision points, the integer values of {2,3, . . . ,7,8> for some values of k, namely for
k E {0,n*,2n2,3n2,4n2,5n2,6nz). Applying
now the integral
equation
(26) at the points
s
Sk, coincide
with
(33)
s = Sk, for k = 0, . . . ,n,, we obtain
rm.i
4SJ = which in discretised
k = 0, . . . ,nl,
K(sk,r)f(r)dr, r”un
form results
(34)
in
d(sk)=
f j=
Kkjfj,
k
=
(35)
O,...,nl,
1
where 5
Kkj
=
K(sk,r)dr,
k = 0, . . . . n,,
j = 1, .
m.
(36)
Equations (11) applied for i = 1 and i = 9, and equations (13) and (35) form a system of (nl + 1) + 1 + 2 = n, + 4 linear equations with m unknowns and the numerical method for determining the unknown vector f applies as described in Section 3. Recall now that from equation (31) n1 should be a multiple of 6 and, in this study a value of n, = 120 was considered sufficiently large such that the number of equations, n, + 4 = 124, exceeded the number of unknowns, m = 100. A further increase in the value of n, may cause incompatibility problems related to the rounding-off errors. Also, we remark that an increase in the number of stages of the diffusion battery, n = 9, is not to be recommended since from equation (29) as PI(r) E [0, l), the values of K;(r) will be virtually zero for i > 9 and will make the system of equations (15) even more singular. In addition, with increasing the number of stages the parasitic particle deposition should be taken into consideration and the analysis becomes more complicated. The matrix K of n + 1 = 10 rows and m = 100 columns defined by equation (14) has the rank 9 and this figure does not increase as the number of stages of the diffusion battery is increased (see Lloyd, 1994), whilst the extended matrix K of n, + 4 = 124 rows and m = 100 columns defined by equation (36) has the rank 19. Therefore, it is to be expected that the extension of the kernel given by equation (27) will improve the estimation of the function to be retrieved. At this stage, it should be noted that equations (34) cannot be used exactly in the absence of the knowledge of function f(r) which is unknown. In fact, only 9 data values are available, passing through the namely di, for i = 1, . . . ,9. However, a piecewise linear interpolation data points di, for i = 2, . . . (8, can be considered. Although some undesirable noise is now added to the data, this penalty was found not to affect significantly the numerical results in comparison with the improvement achieved in the prediction of the analytical solution. Figure 2 shows the piecewise linear approximation passing through 9 exact data values di, for i = 1, . . . ,9, for a double-mode normal distribution (pl = 40 nm, (rl = 10 nm and pL2= 200 nm, g2 = 15 nm) in comparison with the exact values d(s) given by equation (26). Based on this approximation, Fig. 3a shows the variation of f(r) as a function of r obtained by minimising the functional Jl(f) given by expression (20) subject to the bounds of the
Determination
of the particle size distribution
805
0.15-
z n
O.l-
0.05-
i
sb
sb
4b
2b
id0
ii0
S ), the deposition measurements Fig. 2. The deposition measurements values di, for i = 1, .. . ,9, ( ?? function d(s) (--) and its piecewise linear approximation (---) for a double-mode normal distribution (rI = 40 nm, u1 = 10 nm and p2 = 200 nm, o2 = 15 nm).
variables (16) and to the n, + 4 = 124 linear constraints d(SJdi -
f
j=l f
Kkjfj
Kijfj
nl,
(374
i = 1 and i = 9,
WW
G-C,
GE,
k=O ,...,
j=l
1 =
f j=l
&n+l,jfj-
r (in nm) Fig. 3a.
(37c)
D. Lesnic et al.
806
200
300
400
500
400
500
r (in nm)
x10’
6
200
300 r (in nm)
Fig. 3. The variation of f(r) as a function of I obtained using constraints (21) (-- -), using constraints (37) (-+) and the analytical solution (--) for double-mode distributions (a) normal modes (pi = 40 nm, uI = 10 nm and p2 = 200 nm, oz = 15 nm); (b) lognormal modes (pi = 40 nm, u’, = 1.3 and p; = 140 nm, CT;= 1.5); (c) a combined normal mode (pi = 40 nm, u1 = 10 nm) and a lognormal mode (p’, = 140 nm, g’, = 1.5).
This is shown to improve the estimation of the analytical solution in comparison with the numerical results obtained by using only n + 1 = 10 constraints given by expressions (21). Also, the same conclusion may be drawn for a double-mode lognormal OJL;= 40 nm, a; = 1.3 and p; = 140 nm, 05 = 1.5) and both a normal and lognormal mode (pi = 40 nm, e1 = 10 nm and p’r = 140 nm, 0; = 1.5) illustrated in Figs 3b and c, respectively. At this stage it should be noted that the numerical method employed in this study is able to estimate well single- and double-mode distributions which have significantly positive
Determination
of the particle
size distribution
807
values only in the approximate range of [5 nm, 300nm]. The remaining range of particle sizes [I300nm, 500 nm] is associated with the region where, for the actual diffusion battery characteristic design and volumetric flow rate, the fractional penetration functions defined by equation (8) become non-monotonic (see also Wu et al., 1989), and therefore it introduces non-uniqueness features. For other research which correlates the volumetric flow rate to a given particle size range of interest on which the fractional penetration functions become monotonic, see Lesnic et al. (1995). The piecewise linear approximation plotted in Fig. 4 for a triple mode, i.e. p + p’ = 3, normal distribution (,u~ = 40nm, crl = lOnm, pz = 2OOnm, c2 = 15nm and p3 = 360nm, u3 = 20 nm) shows the same characteristics to those presented in Fig. 2 for a double-mode normal distribution. With this approximation the variation of f(r) as a function of r using the new method based on the constraints (37) is shown in Fig. 5a. Although these results are still poor in estimating the analytical solution for the last two peaks, they give more information about the properties of the solution than does the original method, which is based on constraints (21). However, for the triple-mode distribution, better results may be achieved if one uses the numerical method presented previously, but now in a sequential manner. Based on the local recovery ofeach peak of the input distribution it produces better estimations than the global recovery of all the peaks. From Fig. 5a it can be seen that the initial full application of the numerical method, based on constraints (37), provides results which are in very good agreement with the analytical solution for the function f(r) where the first peak is located, and this observation will be further assumed. Thus, if fi(r) defines the numerical solution in an interval [l,in, r*] which extends to zero values in the interval (r*,r,,J, where r* is chosen as the position of a local minimum in the first region, e.g. I* E 100 nm for Fig. 5a, then the function f(r) can be expressed as
r E Crmin, r,,,l,
f(r) = fi (r) + f(r),
(38)
where f(r) is an unknown function to be determined in the interval (r*, rmax]and taken to be zero in the interval [rmin,r*]. On using equation (26) the problem for T(r) reformulates as I* Ima, K(s, r)T(r) dr = d(s) = d(s) - c K(s, r) fi (r) dr, s E C2,Sl. (39) J
Jhl”
I*
0.2
0.15
z ‘0 0.1
I 0
I
I
I
I
I
I
20
40
60
80
100
120
S
Fig. 4. The deposition measurements values di, for i = 1, . ,9, (a), the deposition measurements function d(s) (-) and its piecewise linear approximation (---) for a triple-mode normal distribution (pl = 40 nm, o1 = 10 nm, p12 = 200 nm, CJ*= 15 nm and ps = 360 nm, u3 = 20 nm).
D. Lesnic et al.
808
6.0
260
(4
300
400
500
r (in nm)
6
2
(b)
r (in nm) Fig. Sa and b.
Analogous constraints for f(r) can be deduced from the integral equation (39) as the constraints (37) for the initial function f(r) have been deduced from the integral equation (26). A new application of the method in the interval (r*,r,,,] produced results which estimate well the second region of the function f(r) where the second peak is located. Similarly, as above we can define by f*(r) the numerical solution in the interval @*,I**] which extends to zero values outside of this interval, where r ** is chosen as the position of a local minimum for the second region, e.g. r ** z 250 nm for Fig. 5a. Then the function f(r) can be expressed as f(r) = _Mr) + f’(r),
rE
Crmin, r,,,l~
(40)
Determination of the particle size distribution
200
300 r
O-J
(4
’
809
(in nm)
1
I
I
I
100
200
300
400
1 500
r (in nm) Fig. 5. The variation of f(r) as a function of r obtained using constraints (21) (- --), using constraints (37) (---), the successive algorithm based on the integral equations (26), (39) and (41) (-Cl-) and the analytical solution (---) for (a) three normal modes (pr = 40 nm, o1 = 10 nm, p2 = 200 nm, e2 = 15 nm and p3 = 360 nm, o3 = 20 nm); (b) combined two normal modes (pl = 200 nm, et = 15 nm and JQ = 360 nm, e2 = 20 nm) and one lognormal mode (& = 40 nm, CT;= 1.3); (c) combined two lognormal modes (& = 40 nm, CT;= 1.3 and & = 140 nm, u; = 1.5) and one normal mode (cc1= 200 nm, o1 = 15 nm); (d) three lognormal modes (,u‘, = 40 nm, 6; = 1.3, & = 140 nm, 6; = 1.5 and pi = 240 nm, u; = 1.7).
where the unknown satisfies
function
f(r)
to be determined
in the remaining
r** In?*. K(s, r)f(r) dr = a(s) = a(s) K (s, r) Sz (r) dr, s I. s I** Unfortunately,
a final application
of the method
interval
81.
(41)
equation
(41) in the
s E D,
based on the integral
(r**,r,,,]
810
D. Lesnic et al.
6 z .+ 4
0
100
200
300 r
(in nm)
400
500
Fig. 6. The variation of f(r) as a function of r obtained using constraints (21) (- --), using constraints (37) (-v), the successive algorithm based on the integral equations (26), (39) and (41) (-U-) and the analytical solution (--) for a three normal mode, p1 = 40 nm, o1 = 10 nm, hI = 1, pz = 200 nm, CJ~= 15 nm, hz = 1.5 and bcj = 360 nm, o3 = 20 nm, h, = 2.
interval (r** , r,,,] did not produce a sufficiently well desired estimation-for the third peak. As aforementioned the difficulty in numerically estimating the function f(r) is a result of the fact that the interval (r** ,rmax] = [250 nm, 500 nm] in which the function is sought contains the “problem” range in which the fractional penetration functions become non-monotonic. Therefore, it is to be expected that the solution f=(r) will be a poor approximation since there is insufficient information provided by the integral equation (41) in the interval (r**, I,,, 1. In addition, a large amount of noise might have been accumulated through the successive algorithms procedure. The numerical results obtained by successive application of the initial method based on the integral equations (26), (39) and (41) are included in Fig. 5a, and show a good estimation of the analytical solution at least in the first two regions where the first two peaks are located. The same conclusion may be drawn from Figs 5b and c for the combined triple mode, two normal and one lognormal (pl = 200 nm, c1 = 15 nm, p2 = 360 nm, c2 = 20 nm and p; = 40 nm, a; = 1.3) and the two lognormal and one normal (pl = 200nm, g1 = 15 nm, PL; = 40nm, 0; = 1.3 and p; = 140nm, o$ = 1.5) distributions, respectively. Finally, for completeness, Fig. 5d shows the results for a triple-mode lognormal distribution (pL; = 40 nm, a; = 1.3, p> = 140 nm, a; = 1.5 and p; = 240 nm, o; = 1.7) for which the sequential algorithm applied for other types of triple-mode distributions was unnecessary since in this case the function has actually only two peaks similar to that for a double-mode distribution, In all the calculations undertaken so far, hk, for k = 1, . . . ,p, and hi, for k = 1, . . . , p’, have either been set to zero or unity. However, removing this equal amplitude restriction for any non-zero modes produces, as shown in Fig. 6, where h1 = 1, h2 = 1.5 and h3 = 2, the same success for the sequential method as in Figs 5a-c. From the overall discussion of the results analysed in this section it can be concluded that, for the diffusion battery characteristic design investigated the numerical method introduced in this study estimates well the first two regions, i.e. where the first two peaks are located within the approximate range of [5 nm, 300 nm] for any fine particle size distribution.
Determination of the particle size distribution
811
5. CONCLUSIONS
The retrieval of the unknown aerosol particle size distribution from a set of discrete diffusion battery measurements has been analysed. The numerical method introduced in this study does not make use of any a priori information about the form of the solution and is based on a modified first-order-regularisation procedure by minimising the norm of the first derivative subject to bounds of the variables and linear constraints. In addition, it eliminates the difficulty of choosing the regularisation parameter when the amount of noise included in the data is known in advance. For single-mode distributions (normal or lognormal) having significant positive values only in the range [5 nm, 300 nm] the numerical results agree well with the analytical solutions. When double-mode distributions are present the problem is approximated as a Fredholm integral equation of the first kind for which the numerical solution provided results which are in good agreement with the analytical solutions. Finally, the successive application of the numerical method enabled good results for the first two peaks in a triple-mode distribution to be obtained. A further study, which is at present under investigation, will be concerned with the relations between the diffusion battery characteristic design, the volumetric flow rate, the particle size range of interest and the interval of monotonicity of the fractional penetration functions of aerosols at the stages of the diffusion battery. Acknowledgements-We would like to thank Dr J. J. Lloyd from Manchester Royal Infirmary for all his kind assistance during the preparation of this study and for providing the values of the constants necessary for defining the set of deposition functions in Section 2. Also the very constructive comments of the referees are gratefully acknowledged.
REFERENCES Baklanov, A. M. and Dubtsov, S. N. (1993) An experimental set-up for aerosol spectometers calibration. J. Aerosol Sci. 24, S237-S238. Bashurova, V. S., Koutzenogii, K. P., Pusep, A. Y. and Shokhirev, N. V. (1991) Determination of atmospheric aerosol size distribution functions from screen diffusion battery data: mathematical aspects. J. Aerosol Sci. 22, 373-388. Brent, R. P. (1974), Algorithm 488. Comm. ACM, 704-706. Cheng, Y. S., Keating, J. A. and Kanapilly, G. M. (1980) Theory and calibration of a screen-type diffusion battery, J. Aerosol Sci. 11, 549-556.
Cheng, Y. S. and Yeh, H. C. (1980) Theory of a screen-type diffusion battery. J. Aerosol Sci. 11, 313-320. Cheng, Y. S. and Yeh, H. C. (1983) Performance of a screen-type diffusion battery. In Aerosols in Mining and Industrial Work Enuironments (Edited by Marple, V. A. and Liu, B. Y. H.), pp. 1077-1094. Wiley, New York. Cheng, Y. S., Yeh, H. C. and Brinsko, K. J. (1985) Use of wire screens as a fan model filter. Aerosol Sci. Technol. 4, 165-174.
Cooper, D. W. and Spielman, L. A. (1976) Data inversion using nonlinear programming with physical constraints: aerosol size distribution measurement by impactors. Atmos. Enuir. 10, 723-729. Cooper, D. W. and Wu, J. J. (1990) The inversion matrix and error estimation in data inversion: application to diffusion battery measurements. .I. Aerosol Sci. 21, 217-226. Emi, H., Wang, C. S. and Tien, C. (1982) Transient behaviour of aerosol filtration in model filters. AIChE J. 28, 397-405.
Gill, P. E., Hammarling, S. J., Murray, W., Saunders, M. A. and Wright, M. H. (1986) User’s Guide for LSSOL (Version 1.0). Report SOL 86-1, Stanford University. Heintzenberg, J. (1975) Determination in situ of the size distribution of the atmospheric aerosol. J. Aerosol Sci. 6, 291-303.
Kirsch, A. A. and Stechkina, I. B. (1978) The theory of aerosol filtration with fibrous filters. In Fundamentals of Aerosol Science (Edited by Shaw, D. T.), pp. 165-256. Wiley, New York. Lesnic, D., Elliott, L. and Ingham, D. B. (1995) The inverse determination of the particle size distribution from diffusion battery data. 9th U.K. Ann. Aerosol Sot. Conf., Norwich (submitted). Lloyd, J. J. (1994) Prioate communication. Lloyd, J. J., Taylor, C. J. and Shields, R. A. (1994) Analysis of diffusion battery data using simulated annealing. Pkoc. 8th U.K. Ann. Aerosol Sot. Co& New York, pp. 29-34. Paatero, P. (1990) The extreme value estimation deconvolution method with applications in aerosol research. Univ. Helsinki kept. ‘Ser. Phys., HU-P-250. Sinclair, D. and Hopes, G. S. (1975) A novel form of diffusion battery. Am. Ind. Hyg. Assoc. J. 36, 39-42. Skaptsov, A. S. and Lipatov, G. N. (1994) On the interpretation of rectangular diffusion battery data. J. Aerosol Sci. 25, 341-344.
Soderholm, S. C. (1983) Analyzing screen diffusion battery data. In Aerosols in Mining and Industrial Work Environments (Edited by Marple, V. A. and Liu, B. Y. H.), pp. 1095-1115. Wiley, New York. Tikhonov, A. N. and Arsenin, V. Y. (1977) Solutions of Ill-Posed Problems. W. H. Winston & Sons, Washington.
812
D. Lesnic et al.
Twomey, S. (1963) On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature. J. Assoc. Comput. Mach. 10, 97-101. Twomey, S. (1965) The application of numerical filtering to the solution of integral equations encounted in indirect sensing measurements. .I. Franklin Inst. 279, Y95-109. Twomey, S. (1975) Comparison of constrained linear inversion and an iterative nonlinear algorithm applied to the indirect estimation of particle size distribution. J. Cornput. Plays. 18, 188-200. Wolfenbarger, J. K. and Seinfeld, J. H. (1990) Inversion of aerosol size distribution data. .I. Aerosol Sci. 21, 227-247. Wu, J. J., Cooper, D. W. and Miller, R. J. (1989) Evaluation of aerosol deconvolution algorithms for determining submicron particle size distributions with diffusion battery and condensation nucleus counter. J. Aerosol Sci. Zo, 477-482. Yee, E. (1989) On the interpretation of diffusion battery data. J. Aerosol Sci. 20, 797-811.