- Email: [email protected]
Numerical schemes to model condensation and evaporation of aerosols
Pergamon
Atmospheric Environment Vol. 30, No. 6, pp. 919 928, 1996 Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 1352-2310/96 $15.00 + 0.OA)
1352-2310(95) 00288-X
NUMERICAL SCHEMES TO MODEL C O N D E N S A T I O N A N D EVAPORATION OF AEROSOLS S U R E S H D H A N I Y A L A * a n d A N T H O N Y S. W E X L E R Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, U.S.A.
(First received 13 January 1995 and in final form 10 July 1995) Abstract--A positive-definite scheme is developed to solve the aerosol condensation and evaporation problem. Bott's positive-definite advection scheme (Bott, 1989a) is used to model the movement of the distribution in the size coordinate while analytical solutions are derived to model the growth term. The numerical predictions of the positive-definite scheme are compared with that of previous attempts to model this equation, i.e. a finite-element-finite-difference method and accurate space derivatives scheme (ASD). New analytical solutions are derived for the condensation/evaporation equation and these are used to estimate the accuracy of the schemes for modeling single component and multicomponent cases. The scheme is seen to be accurate at small Courant numbers and is numerically very efficient.
Key word index: Aerosols, condensation, evaporation, numerical methods, positive-definite.
INTRODUCTION
In the study of atmospheric aerosols, much effort has been devoted to understanding the physical and chemical mechanisms that evolve the aerosol size distribution. The chemical composition of a particle often affects its growth rate, and thus the computation of the size distribution often may not be performed independently of that for the composition. However, little work has been done on solving the governing equations for the evolution of the multicomponent chemical composition distribution function. Most of the numerical methods that have been proposed are based on dividing the size distribution into elements or sections, and assuming the form of the size distribution within these elements. The assumption of a constant distribution within a section leads to the so-called "sectional method" (Gelbard et al., 1980; Gelbard and Seinfeld, 1980; Bassett et al., 1981). Some researchers have used fixed sections (Warren and Seinfeld, 1985; Tsang and Rao, 1988), while others have used moving sectional methods (Kim and Seinfeld, 1990; Gelbard, 1990) to model condensation, evaporation and coagulation. To avoid numerical difficulties associated with sectional methods, a continuous distribution may be employed; a mass distribution function is defined for each species and the pertinent differential equations solved to obtain the composition distribution (Pilinis,
* Currently at the Department of Mechanical Engineering, University of Minnesota, Minneapolis, MN 55455, U.S.A.
1990). The numerical schemes that have been used to solve condensation and evaporation, accurate space derivatives (ASD) (Wexler et al., 1994) and finite element (FEM) (Pilinis, 1990), have the disadvantage of yielding negative aerosol concentration values under certain circumstances. In this work we have developed a numerical method, based on Bott's advection scheme (Bott, 1989a, b), that yields positive concentration values, is relatively fast, and has limited numerical diffusion and dispersion.
The governing equation The governing equation for condensation and evaporation of a continuously distributed, internallymixed aerosol is
c~pi 1 aHpi ~-7 = n i p - ~ ~-~
lli
where Pi is the mass distribution of the ith species, such that pid/a is the mass in the size range /a to /1 + d/a, # is the log of the particle diameter, Hi is the condensation/evaporation rate of the ith species of the aerosol, with units of s - 1, H = ~ Hi and p = ~ Pi, and s is the number of aerosol compounds simulated (Wexler et al., 1994). An understanding of the role of each term of equation (1) will facilitate the development of numerical schemes that solve the equation. The first term, Hip, describes the effect of condensation or evaporation on aerosol mass. Condensation (or evaporation) of the ith species, results in the growth (or decay) of the Pi distribution. Simultaneously, the mass of the
919
920
S. DHANIYALA and A. S. WEXLER
individual particles increases, such that the Pi distribution moves in the size coordinate, p. This is represented by the second term, -(1/3)OHpJd#. Advection schemes Considering the above differential equation with only the second term on the right-hand side, we have
Opi 10Hpi . . . . . Ot 3 Op
(2)
which resembles the one-dimensional advection equation often written in generic form as Oc &
Ouc Ox
(3)
where c is a quantity being transported and u is a component of the velocity field. Thus the wealth of numerical methods that have been developed to solve equation (3) can be exploited in the solution of equation (2). While equations (2) and (3) are similar in form, the physical processes they model are not identical. The condensation rate, H, in equation (2) is analogous to the velocity, u, in equation (3), but the velocity field can usually be taken as divergence-free while the values of H (or Hi) can vary arbitrarily in sign and magnitude from one size to the next. Thus, schemes derived under a divergence-free assumption for the velocity field cannot be used to accurately solve the condensation/evaporation equation. Several numerical schemes for advection are available in the literature (e.g. MacCormack, 1969; Gazdag, 1973; Roache, 1976; Bott, 1989a) and a number of reviews are available that compare the important schemes (Chock and Dunker, 1983; Chock 1985, 1991; Rood, 1987). An ideal advective scheme transports the quantity c at the correct velocity, maintains the correct amplitudes, and preserves steep gradients at the correct locations. However, numerical modeling of advection is plagued with difficulties. All schemes have numerical errors, and the major errors are caused by discretizing the continuous variables, which introduces phase, amplitude, and Gibbs errors. Efforts to reduce or correct one error may result in the magnification of another. One of the schemes previously considered to model the condensation/evaporation equation is Accurate Space Derivatives (ASD) (Gazdag, 1973; Wexler et al., 1994). This method is based on the Taylor expansion of the advected quantity about the current time. Each of the temporal derivatives is converted to derivatives in the size coordinate using the advection equation itself and the size coordinate derivatives are then calculated with Fourier transforms. This scheme is very accurate, but time consuming and not positive definite, i.e. quantities might become erroneously negative after numerically advecting them. One of the most popular and simple positivedefinite advection schemes is the first-order upwind differencing scheme. The scheme is explicit, numer-
ically efficient and positive-definite. However, it is very diffusive. Smolarkiewicz (1983, 1984) developed a conservative and positive-definite scheme by introducing corrective advection fluxes to reduce the truncation error caused by the upstream method. Based on the integrated flux schemes of Tremback et al. (1987), Bott (1989a, b) presented a conservative, positive-definite algorithm that is computationally efficient. This scheme maintains the positivity of the advected values by normalizing and limiting the advective fluxes, and by employing specific limiters to avoid negative mixing ratios. It has also been shown to be very accurate, so we have considered BoWs advection scheme to model the second term of equation (1).
NUMERICAL METHODS
The positive-definite scheme To employ a Bott-like method, growth and advective-like processes in equation (1) are split and individually represented by the following differential equations (Yanenko, 1971; Oran and Boris, 1987) 6~pi & = Hip
Opi
(4)
1 OHpi 3 Oft
Ot
(5)
This split-step approach has some obvious limitations. The timestep must be small enough that the properties of the sections do not change appreciably. For small enough timesteps the time-splitting approach has been shown to be very accurate and it gives us the advantage of choosing the best schemes to model the two terms individually. Considering Hi and H to be constant for a given time slice, the solution to equation (4) is Pi(#, t + At) = pi(/t, t) + -~ p(~, t) x [exp[(HAt)-l] pi(/g, t + At) =Pi(/~, t) +H~p(/~, OAt
for H 5 0
(6)
for H =0
(7)
where At is the time increment For the case, H 4: 0, we must also solve equation (5). Using the concept of advective fluxes (Bott, 1989a) and combining the solution for the growth term as obtained in equation (6) gives pd/~, t + At) =Pi(P, t) + -~ p(/l, t) [exp(HAt) --1] At ---(FT+~/: - ~L,~)
A~
(8)
Numerical schemes to model condensation which following Bott (1989a) can be discretized to yield • [Hdj
[pi]7+ ~ = [p~]7 + ~ At
[p]7 [exp(HAt) - 1]
,
- ~ (Fj+ 1/z - V~i"-a/2)
(9)
where A/~ is the spacing between the sections, [H~] s and [p~]~ are the values of H~ and Pi at grid point j and timestep n, respectively, and [p]7 = ~ = l [ p l ] ' ] and [ H ] j = Z~=I [Hi]j. FT+ 1/2 and F 7_ 1/z are the pi-fluxes through the right and left boundary of section L respectively. Using the upstream method to solve equation (5) we have that the p;-flux through the right boundary of the section is
921
change of pi is governed only by the advective term, -(1/3)OHp#Ol~. We can define two integrals I].++1/2 and I].~-1/2 to evaluate the fluxes at the right boundary of thejth section. These integrals represent the flux out of thejth section to the (j + 1)th section and the flux from the (j + l)th section to the jth section, respectively (Bott, 1989a), and can be shown to be 1/2 /*
I
IL,2 =
lj-+ 1/2 =
|
(IO)
where cf = max(167 l, 0) and c; = - min(16)' I, 0), and -" = ½[H]jAt/A#. For the the Courant number is cj sake of simplicity we can omit the use of superscript n from now on. For stability, the numerical scheme must satisfy the Courant-Friedrich-Levy (CFL) criterion, which here implies that c f (or cj-) ~< 1 at each timestep (Bott, 1989a). To understand the physical meaning of this flux term, consider thejth and (j + 1)th sections with [H]j and [H]j+ 1 as the condensation or evaporation rates in these sections. In a situation where the particles in the jth section are experiencing net condensation, i.e. [ H ] j is positive, and the particles in the (j + 1)th section are experiencing a net evaporation, i.e. [H]j+ 1 is negative, the particles in section j increase in size, while the particles in sectionj + 1 decrease. Thus, the distribution in section j moves right towards section j + 1, while the distribution in section j + 1 moves to the left. These movements result in a net flux across the boundary between sectionsj a n d j + 1 represented by Fj+a/2. A similar discussion elucidates the net movement across the boundary between sections j - 1 and j giving the flux Fj_ 1:2. The flux F~+ 1/2 is the prflux at the right boundary of the jth section and the flux F7_ 1/z is the prflux at the left boundary of the jth section. To evaluate these fluxes, consider the case where the condensation rate of the ith species, H , is zero. Thus the temporal
(11)
[p,] (p'(j + 1), t) d#'.
(12)
i,/ -
F;+ , / ~ = - ~ ( c+ [ p,37 - c.- ,[p,37+1)
[p,](#'(j),t)d#'
d 1/2 - c +
1/2
Large numerical diffusion in the evaluation of these integrals is generated by the upstream method due to the poor representation of the [p~] values by a step in each grid box. Accurate calculation of the above integrals requires a higher-order fit for the l'p~] distribution. The method of polynomial-fit was developed by Crowley (1968) and Tremback et al. (1987). Using an area preserving polynomial of order l, the distribution is
[p,] (#' (j), t) = [p,]~.l(ll' ( j ) ) l
= ~ aj",klZ'k(j)
(13)
/¢=0
where the coefficients of the lth order polynomial fit, ai~k, are functions of (l + 1) [p~]-values. The functional dependence is tabulated for l = 2 and l = 4 in Table 1 (Bott, 1989a, b). Substituting this assumed distribution into equation (12) reduces the integrals If+ 1/2 and If+ 1/2 to l
It.~+1/2 = k~=O(k +aJ'kl)2k+ 1 (1 -- (1 -- 2 e l ) k+l ) l
aj+l.k
(14) ) k + 1).
ll.j+:/z = ~ (k + 1)2k+a (--1)k(1-- (1 -- 2Cj-+1 k=O
(15) The total area of the discrete [Pi] distribution is [pi] A/~ and that of the continuous polynomial distri-
Table 1. Coefficients a~.kfor the I = 2 and I = 4 versions of Bott's area-preserving flux-form algorithm (after Bott, 1989b) 1=2 1
1=4 -- ll6[pi]~+l + 2134[p~]i- ll6[pi]i-x + 9[pi]i-2)
ai.0
~([Pi],+ 1 --26[p,]/+ [Pi],- 1)
1-~2o (9[pi]i+2
ai.l
1 2(I-Pl]i+ 1 -- [pi],-l)
ai.2
½([Pl]i÷l -- 2[pi]i + [P~]i-I)
( - 5 [.PJi+2 + 34 [p,], + 1 - 34 [p,],_ 1 + 5 [p.,],_ z) ~s(-3[P,],+2 - 36[p,],+~ - 661-p,], + 36[p,],_~ -3[pi],-2)
at. 3 ai. 4
~([P,],+2 -2[pi],+1 + 2[pill-1 - [P,],-2) ~([P,],+z-4[p,],+l + 6 [ p , ] , - 4[p,],_, + [P,],-2)
922
S. DHANIYALA and A. S. WEXLER
Finite element method
bution is given by the integral,
We have also used a finite element method (FEM) for purposes of comparison with the results of Pilinis (1990). Equation (1) can be rewritten as
1/2 ¢*
[,,j = | [p,] (~'(j), t) dr' QI
-- 1/2
api
aj, k = y, (k+1)+2~+,((-1)'+1). l
(16)
k=O
The resulting difference in the areas can be accounted for by considering the weighting factor [pi]jIL~, and multiplying it by the [pi] (#'(J), t) values. A sufficient condition for maintaining positivity in pure advection (Bott, 1989a) is given by +
0 <~IL~+,/2 + It.j-,~ 2 ~ [Pi]j
(17)
which maintains a positive flux out of a section and limits it by the value of [Pi]. However, in the condensation/evaporation equation we have an additional term Hip. When the ith chemical component is evaporating from the aerosol, i.e. Hi is negative, the sum of the advective flux and the flux due to evaporation must not exceed the value of [Pi] in the section. If we are to maintain the positivity of the concentration values we need to take into account this decay term, and the positive-definite condition (17) becomes
+
= o.
(23)
Using a weighted residual formulation and considering the Galerkin method for the weight functions, the F E M formulation gives
(A + ~ OAt(B + C)) P~+1 1 =(A--~At(1-O)(B+C))p, - Atl~=l [pl]~D, (24) where A, B, C and Di are matrices whose elements are #2 71
[Ajd= | [N,]j [Nd~dr #1 /z2
[Bjkl = f [Ni]j[Ni]k~;---~lt d# Itl
0 <<.It~j+1/2 + l~Tj-~/2 <~[Pi]j [Hi]j ~ + ~
H 6~pl Pi aH
0-7 - Hip + 3-
#2 i!
LP-b (exp(HAt) - 1).
(18)
O[Ndj [Cj~] = 3| H al~ [Ni]~ dr l*l
Returning to the advective fluxes, we have
F$+,/2 = ~ L~LJ
[pi]j
]l,j+l [Pilj+l
~2 tt
•
(19)
[Djk]l = I Hi[NIl1 [Nl]kdp lal
The flux Fi+ 1/2 depends on the value and sign of Hj at the right boundary. F~+,/2 is set to zero if it is not of the same sign as Hj as it represents the flux due to Hj at that boundary. Finally the flux Fj+ ,/2 can be written as Fj+ 1/2
A# = "~
+ (~l,j+
1/2~01,j - - ~ l _ j + 1 / 2 ~ 1 , j + 1)
(20)
where 31.J+~/2=max( + 0,IL~+ + ~/2), 3#j+. 1/2 = m a x (0, I1~+,/2) and gdz.j= [pi]j/~Lj. If [Hi] is negative, component i is evaporating and
Hi -~ p(p, t) (exp(nAt) - 1) + e)
(21)
where e is a small number introduced to avoid numerical difficulties when ~,~ is nearly zero. If [Hi] is positive, component i is condensing and 3,,~ = max(I~./, 3,+j+ 1/2
+ ~ / ~ j - 1/2 + e).
for j, k = 1, ..., NEN, where N E N is the number of elemental nodes. In each section, we have considered the shape function to be quadratic. The C r a n k Nicolson scheme (0 = 1/2) is used for all cases, which makes the scheme stable and gives a high order of accuracy.
(22)
The limit for 3~.i in equation (22) considers the concentration present in the previous time slice as the limiting factor. An alternate formulation would relax this limit in consideration of the condensational flux.
NUMERICAL TESTS AND ANALYTICAL SOLUTIONS
The three numerical schemes, the Bott-like positivedefinite method, the ASD method and the FEM, are compared for three cases: pure advection, pure growth and simultaneous growth and advection. The accuracy of the numerical schemes is estimated by evaluating the root mean square (r.m.s) errors for each of the three test cases defined as NN
n
r.m.s, error = x/~/=1 ([P,]j - [pi]~,)2
I~p,(u)du
where NN is the number of nodes or sections, [Pi]7 and [p~]~." are numerical and exact solutions at timestep n and section j, respectively. The accuracy of the numerical schemes is a function of two dimensionless groups--the Courant number and the ratio of the
Numerical schemes to model condensation inter-nodal spacing, A#, to the characteristic width of the distribution. Characteristic width is defined as the width of the distribution where the value of p~ is 1/e of its maximum value. Consider the case where H~ is zero, but H is not; that is component i is not evaporating or condensing but one or more of the other components are, such that there is a net growth or shrinkage of the particle. In this case, the governing equation for the mass distribution function of the ith species is given by equation (5). The performance of the scheme can be evaluated using several standard advection tests. The problems most commonly associated with advective schemes are dispersion and diffusion. Using this scheme to advect sine, square and triangular distributions with varying condensation rates will help elucidate the effect of these errors. Figures 1-3 highlight the dispersive errors in the three schemes during the advection of a square wave at large Courant numbers (Courant number = ~ [H]j At/Al~ = 0.99). ASD over- and under-shoots the discontinuities, which leads to negative concentrations. Bott's scheme
923
overshoots near the discontinuities but retains the positivity of the concentration values. Like ASD, FEM suffers from dispersive errors that lead to negative values. These observations are summarized in Fig. 4 where we plot the r.m.s error against Courant number for each scheme. Reducing the timestep reduces the errors in the Bott and FEM scheme. Typically the largest numerical diffusion is produced during the first few advection steps when the discontinuities of the initial distribution are smoothed. In the special case where the net condensation/ evaporation rate of the aerosol, H, is zero, but the condensation or evaporation rate of any individual chemical component, H~, is non-zero, we have the case of pure growth of the p~ distribution, and the governing differential equation is given by equation (4). The p~ distribution changes such that the area under the curve is increasing or decreasing, representing an increase in mass of the ith species, but there is no accompanied movement of the distribution in the size coordinate, as the net particle size is unchanged. Again, the growth/decay of sine, square and triangular distributions are simulated.
1.2 Initial
analytical solution - ASD solution .....
Final
1
analytical solution - -
'~.
1.2 1.4 i 0.8
0.8
0.6
o.,
0.2
=~
. . . . . . . . . . . Li.= . . . . . . . . . =.~,.
-2
' -1
' 0
1
i
,
r
J
2
3
4
5
t
i
61' :g :: j
:~ e -0,2
o
1--li i
0.2 0
solution
~
0.6
0.4
FEM
~ Initial
-02
i~
-0.4 ~ -2
-1
0
1
2
3
4
5
6
7
u
Fig. 1. Pure advection: a comparison of the ASD's predictions with the analytical solution at a Courant number of 0.99. 1.2
Initial 1
analytical solution - B o l t ' s scheme -~....
Final ........ i
Fig. 3. Pure advection: a comparison of the FEM's results with the analytical solution at a Courant number of 0.99. 1.8 1.6
Bolt's scheme ~ .... FEM scheme --~ .... ASD scheme --
",.
1.4 0.8
1.2 1
0.6 0.8 0.4
0.6 0.4
0.2 0.2
-2
-1
0
2
3
4
5
0 0,1
0.2
0.3
0.4
0.5
0.6
0,7
0.8
0.9
1
Courant number
Fig. 2. Pure advection: a comparison of the results obtained using Bott's scheme with the analytical solution at a Courant number of 0.99.
Fig. 4. Pure advection: a plot of the Courant number vs r.m.s, error-accuracy estimation of the schemes for a square wave.
924
S. DHANIYALA and A. S. WEXLER 200
0.005 BOR'S scheme ~--. FEM scheme .-=.... ASD scheme . . . .
0.0045 0.004
150
0.0035 0.003 ~ 0.0025 E ~: 0.002
1O0
0.0015 S0
0.001 0.0005 0
. . . . .
=
-o.oo;
0
;.o~
. . . .
=,
o.o-~5 -o.oi
=
o.o25
0
o.o3
-2
-I
0 Ix
Courant number
Fig. 5. Pure growth: a plot of the Courant number vs r.m.s. error-accuracy estimation of the schemes for a square wave. 200
,
,
1
2
3
Fig. 7. Condensation/evaporation for single species aerosol: a comparison of the results obtained using BoWs scheme with the first analytical solution [equation (25)]; Courant number = 0.1.
analytical solution - . ASD solution ~ - -
150
200
Final
analytical solution - FEM solution • ~....
/~
Final
150 100
I
!ii 100
50
0
............. i
-2
-1
0 ~t
i
i
1
2
Initial 3
0I -2
Fig. 6. Condensation/evaporation for single species aerosol: a comparison of the ASD's predictions with the first analytical solution [equation (25)]; Courant number = 0.1. For pure growth test (H is zero, but Hi is not), a two component system is used such that one compound is condensing on the particle and the other evaporating such that the condensation/evaporation rate of the particle is zero (i.e., H1 = - H 2 ; H =0). The tests show that the schemes model pure growth very accurately, but ASD and FEM may yield negative values, because they do not limit the mass of the evaporating species to the mass initially present. The flux limiters in the positive-definite scheme maintain positivity. For pure growth of a square wave, the r.m.s, errors are plotted against the Courant numbers in Fig. 5, where Courant number for the pure growth is defined as H~At. In the above two cases, pure growth and pure advection have been considered mainly for model evaluation. Generally in atmospherically relevant processes, the growth term and the advection-like term are both non-zero. Analytical solutions for pi are not available but solutions for p = ~ i pi have been derived: p(#, t) = 3 exp(6g) e x p ( - exp (3#) e x p ( - Ht) - Ht) (25) p(p, t) = exp(3/z)exp(- e x p ( - 3 g ) e x p ( - Ht)).
(26)
:=~==__________=====:
i
i
i
i
-1
0
1
2
3
Fig. 8. Condensation/evaporation for single species aerosol: a comparison of the FEM's results with the first analytical solution [equation (25)]; Courant number = 0.1.
Bolt's scheme
FEM scheme --~.... ASD scheme .....
// 1.5
/)
/-
.................%11 ........................ ..J-.o-_ ......... . ......... . 0.5
0
0
0-1
012
0.3
0.4
0.5
0.6
0.7
0.S
0.g
Courant number
Fig. 9. Condensation/evaporation for single species aerosol: a plot of the Courant number vs r.m.s, error-accuracy estimation of the schemes for the first analytical solution [equation (25)]. Equation (25) was derived by Gelbard and Seinfeld (1980), while equation (26) was derived for use in this work. For a single component aerosol undergoing
Numerical schemes to model condensation 0.12
925
140 .....
Bott ASD ...... FEM ..6....
0.1
Final
120
anal'flical solution - Bott's scheme -+....
100
0.08
80 0.06 ,
iT60
I 0.04
40 0.02
ol
20 :
0.15
...........
0.2
I
~
0.25
0.3
_ i
,
0.35
0.4
0.45
0,5
0.55 -2
-1
0
inter-nodal spacing/characteristic width
1
2
3
F
Fig. 10. Condensation/evaporation for single species aerosol: a plot of the ratio of the characteristic width to section size (A~) vs r.m.s, error-accuracy estimation of the schemes for the first analytical solution [equation (25)].
Fig. 12. Condensation/evaporation for single species aerosol: a comparison of the results obtained using Bott's scheme with the second analytical solution [equation (26)]; Courant number = 0.l. 140
14o [----
,
analytical solution
'
120 ~
-ASD solution . . . .
Final / ~
,o0
Final .
120
analytical solution - FEM solution t~ .
100
/1
80 60 40 20
20
[ni '
_ ol ...... -2
,
-1
,
0
...
-i
1
:~±
i,
0,=-
.... i ......... i ......... l
-20 -2
-1
0 u
2
1
2
3
la
Fig. 11. Condensation/evaporation for single species aerosol: a comparison of the ASD's predictions with the second analytical solution [equation (26)]; Courant number = 0.1.
Fig. 13. Condensation/evaporation for single species aerosol: a comparison of the FEM's results with the second analytical solution [equation (26)]; Courant number = 0.1. 1.2
condensation, results of the numerical predictions for the first analytical solution are illustrated in Figs 6 - 8 and for the second in Figs 11-13. The results shown are for test cases with Courant number of 0.1 (Courant number = ½ [ H ] j At/A/~). Figures 9 and 10 show the r.m.s error dependence on Courant number and internodal spacing, respectively, for the first analytical solution [equation (25)], while Figs 14 and 15 show the same for the second analytical solution [equation (26)]. Since mass distributions given by the two analytical solutions are Coo smooth functions, A S D models the condensation process for this initial distribution very accurately. However, time-splitting and time-stepping errors associated with the positivedefinite scheme and FEM, respectively, reduce their accuracy. Although analytical solutions are unavailable for the distribution of the individual species in a multicomponent aerosol, the total distribution of all the species of the aerosol can be compared with analytical
~o~,,sc.=na- ..... t FEM scheme .-,.--. J ASD scheme -*- ;/t
1
0.8
................ii / ...........
0,6 J
j,,"
0.4
..........i -....... ......... ~
t
.-
m-..0,2
0
0~;"-';.~"-;~;--0.~ o~5-0~.~-'";.V-o~.i ;.8 Courant number
Fig. 14. Condensation/evaporation for single species aerosol: a plot of the Courant number vs r.m.s, error-accuracy estimation of the schemes for the second analytical solution I-equation (26)]. p distributions. The initial distribution of the three species are considered identical and similar to the distribution obtained from the first analytical
926
S. DHANIYALA and A. S. WEXLER 0.05
,
.
'
eott
0.045
! i
ASD - 4 . . . . FEM --~....
0.04 0.035
~ ~
'
h_,,.o~ ~ 8 o t t t o k a ~ .Q--
hal
/
:::SSI......
0.03 0.025
..,.
0.02
.:.-"
,,
;:/ 0
........./," "/
...a ........ ..._.4-
0.015
5
0.01 0.005
~
5 0
0.05
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0.55
inter-nodal spacing/characteristicwidth
t
-5 -2
i
-1
i
o
t
1
2
i.t
Fig. 15. Condensation/evaporation for single species aerosol: a plot of the ratio of the characteristic width to section size (A#) vs. r.m.s, error-accuracy estimation of the schemes for the second analytical solution [equation (26)].
i
'
. -
Fig. 18. Multispecies aerosols: plot of growth of the mass distribution of the third component accompanied by an increase in particle size (H3 = 10; H = 9) as modeled by the three schemes for the initial distribution of the species given by the first analytical solution [equation (25)].
Bott's scheme FEM scheme --=.... ASD scheme ......
1
Fired
0.8
/ #/" /,
0.6
/,,~'-..
0.4
..**'/"" .i
i
-2
-1
1
........ i .........
0.2
2
~nitial
Final
FEM md~Jon -*,-aottm=~n .=..
g
Fig. 17. Multispecies aerosols: plot of the movement of the mass distribution of the second component due to the increase in particle size (H2 = 0; H = 9) as modeled by the three schemes for the initial distribution of the species given by the first analytical solution [equation (25)].
. .e'""'"
-
....,"
"'",0""""
..... ..~r'"
e'"'""
0
Fig. 16. Multispecies aerosols: plot of the decay of the mass distribution of the first component accompanied by particle growth (Ht = - 1; H = 9) as modeled by the three schemes for initial distribution of the species given by the first analytical solution [equation (25)].
//
ml
,./
01 0.2 03 0.4 0.5 0.6 0.7 o.s 0.9 Courant number
Fig. 19. Multispecies aerosols: plot of the numerical errors with Courant number for the three schemes for the first analytical solution [equation (25)].
solqtion. The first species is considered to evaporate with a rate of H t = - 1 , the second species neither evaporates nor condenses, i.e. H2 = 0, and the third species is assumed to grow at a constant condensation rate of H3 = 10. The total condensation/evaporation rate is H = 9. The change of the mass distributions of the three species due to the different condensation and evaporation rates is shown in Figs 16-18. The r.m.s, error dependence on the Courant number for the total mass distribution (p), shown in Figs 19 and 20 indicates that the positive-definite scheme models the multicomponent condensation/evaporation process accurately at small Courant numbers. The runtimes of the schemes for modeling condensation/evaporation of a three-component species whose total mass distribution p is given by the first analytical solution [equation (25)], is shown in Fig. 21. The positive-definite scheme is seen to be the least
Numerical schemes to model condensation 0.7 BoWs scheme
FEM scheme ASD scheme
0.6 0.5
,/"
0.4 0.3 7
a
....
~r
0.2 0.1
927
into Eulerian air pollution models to predict the size and composition distribution of atmospheric aerosols. The Bott's scheme used in this work has been modified by several researchers to improve the scheme's stability and monotonicity. We have, however, considered the original version of Bott's scheme and proven its application to obtain a positive-definite scheme to model condensation/" evaporation in aerosols. For some applications, incorporating these improvements to Bott's scheme will improve performance.
0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Courant number
Fig. 20. Multispecies aerosols: plot of the numerical errors with Courant number for the three schemes for the second analytical solution [equation (26)1.
18 BoWs scheme ASD scheme ---.... FEM scheme .......
16 14 12 k 10 8 6 4
........... e ...............
2 0 0.1
~ 0.2
-,¢----,--._,_;...........,,-.........-_,___;........ 0.3
0.4
0.5
0.6
0.7
0.8
0.9
Courant number
Fig. 21. Runtimes: a comparison of the time taken by the three schemes for advecting a multicomponent species, whose total mass distribution function p is given by the first analytical solution [equation (25)], through one cell width for varying Courant numbers.
cpu intensive. All the tests were conducted on a S U N 4/670MP workstation.
CONCLUSIONS Three numerical methods were investigated for solving the internally mixed condensation/evaporation equation. Tests were performed on the ability of the methods to represent growth and decay of the size distributions, their shift in the size coordinate, and both processes simultaneously. ASD, a spectral method, is shown to be the slowest and generally the most accurate, but it suffers from Gibbs errors and often predicts unphysical negative concentrations. The finite element method suffered from large dispersive errors for all but the smoothest distributions. A new method developed here and based on Bott's advection scheme is relatively fast and accurate especially at low Courant number, and is positive-definite. This method will be incorporated
REFERENCES
Bassett M., Gelbard F. and Seinfeld J. H. (1981) Mathematical model for multi-component aerosol formation and growth in plumes. Atmospheric Environment 15, 2395-2406. Bott A. (1989a) A positive definite advection scheme obtained by nonlinear renormalization of the advective fluxes. Mon. Weath. Rev. 117, 1006-1015. Bott A. (1989b) Reply. Mort. Weath. Rev. 117, 2633-2636. Chock (1991) A comparison of numerical methods for solving the advection equation-II. Atmospheric Environment 25A, 853-871. Chock D. P. (1985) A comparison of numerical methods for solving the advection equation-II. Atmospheric Environment 19, 571-586. Chock and Dunker (1983) A comparison of numerical methods for solving the advection equation-lL Atmospheric Environment 17, 11-24. Crowley W. P. (1968) Numerical advection experiments. Mon. Weath. Rev. 96, 1-11. Gazdag J. (1973) Numerical convective schemes based on accurate computation of space derivatives. J. Comput. Phys. 13, 100-113. Gelbard F. (1990) Modeling multicomponent aerosol particle growth by vapor condensation. Aerosol Sci. Technol. 12, 399-412. Gelbard F. and Seinfeld J. H. (1980) Simulation of multicomponent aerosol dynamics. J. Colloid Interface Sci. 78, 485-501. Gelbard F., Tambour Y. and Seinfeld J. H. (1980) Sectional representations for simulating aerosol dynamics. J. Colloid Interface Sci. 76, 541-556. Kim Y. P. and Seinfeld J. H. (1990) Simulation of multicomponent aerosol condensation by the moving sectional method. J. Colloid Interface Sci. 135, 185-199. MacCormack R. W. (1969) The effect of viscosity in hypervelocity impact cratering. A I A A Pap. 69-354. Oran S. E. and Jay Boris P. (1987) Numerical Simulation of Reactive Flow. Elsevier, Amsterdam. Pilinis C. (1990) Derivation and numerical solution of the species mass distribution equations for multicomponent particulate systems. Atmospheric Environment 24A(7), 1923-1928. Roache P. J. (1976) Computational Fluid Dynamics. Hermosa. Rood, R. B. (1987) Numerical advection algorithms and their role in atmospheric transport and chemistry models. Rev. Geophys. 25(1), 71 100. Smolarkiewicz P. K. (1983) A simple positive-definite advection scheme with small implicit diffusion. Mon. Weath. Rev. 111, 479-486. Smolarkiewicz P. K. (1984) A fully multidimensional positive definite advection transport algorithm with small implicit diffusion. J. Comput. Phys. 54, 325-362. Tremback C. J., Powell J., William R. C. and Pielke A. R. (1987) The forward-in-time upstream advection scheme:
928
S. DHANIYALA and A. S. WEXLER
extension to higher orders. Mon. Weath. Rev. 115, 540-555. Tsang T. H. and Rao A. (1988) Comparison of different numerical schemes for condensational growth of aerosols. Aerosol Sci. Technol. 9, 271-277. Warren D. R. and SeinfeldJ. H. (1985) Simulation of aerosol size distributionevolution in systems with simultaneous
nucleation, condensation and coagulation. Aerosol Sci. Technol. 4, 31-43. Wexler A. S., Lurmann F. W. and Seinfeld J. H. (1994) Modeling urban and regional aerosols-I. Model development. Atmospheric Environment 28(3), 531-546. Yanenko N. N. (1971) The Method of Fractional Steps. Springer, Berlin.
Atmospheric Environment Vol. 30, No. 6, pp. 919 928, 1996 Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 1352-2310/96 $15.00 + 0.OA)
1352-2310(95) 00288-X
NUMERICAL SCHEMES TO MODEL C O N D E N S A T I O N A N D EVAPORATION OF AEROSOLS S U R E S H D H A N I Y A L A * a n d A N T H O N Y S. W E X L E R Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, U.S.A.
(First received 13 January 1995 and in final form 10 July 1995) Abstract--A positive-definite scheme is developed to solve the aerosol condensation and evaporation problem. Bott's positive-definite advection scheme (Bott, 1989a) is used to model the movement of the distribution in the size coordinate while analytical solutions are derived to model the growth term. The numerical predictions of the positive-definite scheme are compared with that of previous attempts to model this equation, i.e. a finite-element-finite-difference method and accurate space derivatives scheme (ASD). New analytical solutions are derived for the condensation/evaporation equation and these are used to estimate the accuracy of the schemes for modeling single component and multicomponent cases. The scheme is seen to be accurate at small Courant numbers and is numerically very efficient.
Key word index: Aerosols, condensation, evaporation, numerical methods, positive-definite.
INTRODUCTION
In the study of atmospheric aerosols, much effort has been devoted to understanding the physical and chemical mechanisms that evolve the aerosol size distribution. The chemical composition of a particle often affects its growth rate, and thus the computation of the size distribution often may not be performed independently of that for the composition. However, little work has been done on solving the governing equations for the evolution of the multicomponent chemical composition distribution function. Most of the numerical methods that have been proposed are based on dividing the size distribution into elements or sections, and assuming the form of the size distribution within these elements. The assumption of a constant distribution within a section leads to the so-called "sectional method" (Gelbard et al., 1980; Gelbard and Seinfeld, 1980; Bassett et al., 1981). Some researchers have used fixed sections (Warren and Seinfeld, 1985; Tsang and Rao, 1988), while others have used moving sectional methods (Kim and Seinfeld, 1990; Gelbard, 1990) to model condensation, evaporation and coagulation. To avoid numerical difficulties associated with sectional methods, a continuous distribution may be employed; a mass distribution function is defined for each species and the pertinent differential equations solved to obtain the composition distribution (Pilinis,
* Currently at the Department of Mechanical Engineering, University of Minnesota, Minneapolis, MN 55455, U.S.A.
1990). The numerical schemes that have been used to solve condensation and evaporation, accurate space derivatives (ASD) (Wexler et al., 1994) and finite element (FEM) (Pilinis, 1990), have the disadvantage of yielding negative aerosol concentration values under certain circumstances. In this work we have developed a numerical method, based on Bott's advection scheme (Bott, 1989a, b), that yields positive concentration values, is relatively fast, and has limited numerical diffusion and dispersion.
The governing equation The governing equation for condensation and evaporation of a continuously distributed, internallymixed aerosol is
c~pi 1 aHpi ~-7 = n i p - ~ ~-~
lli
where Pi is the mass distribution of the ith species, such that pid/a is the mass in the size range /a to /1 + d/a, # is the log of the particle diameter, Hi is the condensation/evaporation rate of the ith species of the aerosol, with units of s - 1, H = ~ Hi and p = ~ Pi, and s is the number of aerosol compounds simulated (Wexler et al., 1994). An understanding of the role of each term of equation (1) will facilitate the development of numerical schemes that solve the equation. The first term, Hip, describes the effect of condensation or evaporation on aerosol mass. Condensation (or evaporation) of the ith species, results in the growth (or decay) of the Pi distribution. Simultaneously, the mass of the
919
920
S. DHANIYALA and A. S. WEXLER
individual particles increases, such that the Pi distribution moves in the size coordinate, p. This is represented by the second term, -(1/3)OHpJd#. Advection schemes Considering the above differential equation with only the second term on the right-hand side, we have
Opi 10Hpi . . . . . Ot 3 Op
(2)
which resembles the one-dimensional advection equation often written in generic form as Oc &
Ouc Ox
(3)
where c is a quantity being transported and u is a component of the velocity field. Thus the wealth of numerical methods that have been developed to solve equation (3) can be exploited in the solution of equation (2). While equations (2) and (3) are similar in form, the physical processes they model are not identical. The condensation rate, H, in equation (2) is analogous to the velocity, u, in equation (3), but the velocity field can usually be taken as divergence-free while the values of H (or Hi) can vary arbitrarily in sign and magnitude from one size to the next. Thus, schemes derived under a divergence-free assumption for the velocity field cannot be used to accurately solve the condensation/evaporation equation. Several numerical schemes for advection are available in the literature (e.g. MacCormack, 1969; Gazdag, 1973; Roache, 1976; Bott, 1989a) and a number of reviews are available that compare the important schemes (Chock and Dunker, 1983; Chock 1985, 1991; Rood, 1987). An ideal advective scheme transports the quantity c at the correct velocity, maintains the correct amplitudes, and preserves steep gradients at the correct locations. However, numerical modeling of advection is plagued with difficulties. All schemes have numerical errors, and the major errors are caused by discretizing the continuous variables, which introduces phase, amplitude, and Gibbs errors. Efforts to reduce or correct one error may result in the magnification of another. One of the schemes previously considered to model the condensation/evaporation equation is Accurate Space Derivatives (ASD) (Gazdag, 1973; Wexler et al., 1994). This method is based on the Taylor expansion of the advected quantity about the current time. Each of the temporal derivatives is converted to derivatives in the size coordinate using the advection equation itself and the size coordinate derivatives are then calculated with Fourier transforms. This scheme is very accurate, but time consuming and not positive definite, i.e. quantities might become erroneously negative after numerically advecting them. One of the most popular and simple positivedefinite advection schemes is the first-order upwind differencing scheme. The scheme is explicit, numer-
ically efficient and positive-definite. However, it is very diffusive. Smolarkiewicz (1983, 1984) developed a conservative and positive-definite scheme by introducing corrective advection fluxes to reduce the truncation error caused by the upstream method. Based on the integrated flux schemes of Tremback et al. (1987), Bott (1989a, b) presented a conservative, positive-definite algorithm that is computationally efficient. This scheme maintains the positivity of the advected values by normalizing and limiting the advective fluxes, and by employing specific limiters to avoid negative mixing ratios. It has also been shown to be very accurate, so we have considered BoWs advection scheme to model the second term of equation (1).
NUMERICAL METHODS
The positive-definite scheme To employ a Bott-like method, growth and advective-like processes in equation (1) are split and individually represented by the following differential equations (Yanenko, 1971; Oran and Boris, 1987) 6~pi & = Hip
Opi
(4)
1 OHpi 3 Oft
Ot
(5)
This split-step approach has some obvious limitations. The timestep must be small enough that the properties of the sections do not change appreciably. For small enough timesteps the time-splitting approach has been shown to be very accurate and it gives us the advantage of choosing the best schemes to model the two terms individually. Considering Hi and H to be constant for a given time slice, the solution to equation (4) is Pi(#, t + At) = pi(/t, t) + -~ p(~, t) x [exp[(HAt)-l] pi(/g, t + At) =Pi(/~, t) +H~p(/~, OAt
for H 5 0
(6)
for H =0
(7)
where At is the time increment For the case, H 4: 0, we must also solve equation (5). Using the concept of advective fluxes (Bott, 1989a) and combining the solution for the growth term as obtained in equation (6) gives pd/~, t + At) =Pi(P, t) + -~ p(/l, t) [exp(HAt) --1] At ---(FT+~/: - ~L,~)
A~
(8)
Numerical schemes to model condensation which following Bott (1989a) can be discretized to yield • [Hdj
[pi]7+ ~ = [p~]7 + ~ At
[p]7 [exp(HAt) - 1]
,
- ~ (Fj+ 1/z - V~i"-a/2)
(9)
where A/~ is the spacing between the sections, [H~] s and [p~]~ are the values of H~ and Pi at grid point j and timestep n, respectively, and [p]7 = ~ = l [ p l ] ' ] and [ H ] j = Z~=I [Hi]j. FT+ 1/2 and F 7_ 1/z are the pi-fluxes through the right and left boundary of section L respectively. Using the upstream method to solve equation (5) we have that the p;-flux through the right boundary of the section is
921
change of pi is governed only by the advective term, -(1/3)OHp#Ol~. We can define two integrals I].++1/2 and I].~-1/2 to evaluate the fluxes at the right boundary of thejth section. These integrals represent the flux out of thejth section to the (j + 1)th section and the flux from the (j + l)th section to the jth section, respectively (Bott, 1989a), and can be shown to be 1/2 /*
I
IL,2 =
lj-+ 1/2 =
|
(IO)
where cf = max(167 l, 0) and c; = - min(16)' I, 0), and -" = ½[H]jAt/A#. For the the Courant number is cj sake of simplicity we can omit the use of superscript n from now on. For stability, the numerical scheme must satisfy the Courant-Friedrich-Levy (CFL) criterion, which here implies that c f (or cj-) ~< 1 at each timestep (Bott, 1989a). To understand the physical meaning of this flux term, consider thejth and (j + 1)th sections with [H]j and [H]j+ 1 as the condensation or evaporation rates in these sections. In a situation where the particles in the jth section are experiencing net condensation, i.e. [ H ] j is positive, and the particles in the (j + 1)th section are experiencing a net evaporation, i.e. [H]j+ 1 is negative, the particles in section j increase in size, while the particles in sectionj + 1 decrease. Thus, the distribution in section j moves right towards section j + 1, while the distribution in section j + 1 moves to the left. These movements result in a net flux across the boundary between sectionsj a n d j + 1 represented by Fj+a/2. A similar discussion elucidates the net movement across the boundary between sections j - 1 and j giving the flux Fj_ 1:2. The flux F~+ 1/2 is the prflux at the right boundary of the jth section and the flux F7_ 1/z is the prflux at the left boundary of the jth section. To evaluate these fluxes, consider the case where the condensation rate of the ith species, H , is zero. Thus the temporal
(11)
[p,] (p'(j + 1), t) d#'.
(12)
i,/ -
F;+ , / ~ = - ~ ( c+ [ p,37 - c.- ,[p,37+1)
[p,](#'(j),t)d#'
d 1/2 - c +
1/2
Large numerical diffusion in the evaluation of these integrals is generated by the upstream method due to the poor representation of the [p~] values by a step in each grid box. Accurate calculation of the above integrals requires a higher-order fit for the l'p~] distribution. The method of polynomial-fit was developed by Crowley (1968) and Tremback et al. (1987). Using an area preserving polynomial of order l, the distribution is
[p,] (#' (j), t) = [p,]~.l(ll' ( j ) ) l
= ~ aj",klZ'k(j)
(13)
/¢=0
where the coefficients of the lth order polynomial fit, ai~k, are functions of (l + 1) [p~]-values. The functional dependence is tabulated for l = 2 and l = 4 in Table 1 (Bott, 1989a, b). Substituting this assumed distribution into equation (12) reduces the integrals If+ 1/2 and If+ 1/2 to l
It.~+1/2 = k~=O(k +aJ'kl)2k+ 1 (1 -- (1 -- 2 e l ) k+l ) l
aj+l.k
(14) ) k + 1).
ll.j+:/z = ~ (k + 1)2k+a (--1)k(1-- (1 -- 2Cj-+1 k=O
(15) The total area of the discrete [Pi] distribution is [pi] A/~ and that of the continuous polynomial distri-
Table 1. Coefficients a~.kfor the I = 2 and I = 4 versions of Bott's area-preserving flux-form algorithm (after Bott, 1989b) 1=2 1
1=4 -- ll6[pi]~+l + 2134[p~]i- ll6[pi]i-x + 9[pi]i-2)
ai.0
~([Pi],+ 1 --26[p,]/+ [Pi],- 1)
1-~2o (9[pi]i+2
ai.l
1 2(I-Pl]i+ 1 -- [pi],-l)
ai.2
½([Pl]i÷l -- 2[pi]i + [P~]i-I)
( - 5 [.PJi+2 + 34 [p,], + 1 - 34 [p,],_ 1 + 5 [p.,],_ z) ~s(-3[P,],+2 - 36[p,],+~ - 661-p,], + 36[p,],_~ -3[pi],-2)
at. 3 ai. 4
~([P,],+2 -2[pi],+1 + 2[pill-1 - [P,],-2) ~([P,],+z-4[p,],+l + 6 [ p , ] , - 4[p,],_, + [P,],-2)
922
S. DHANIYALA and A. S. WEXLER
Finite element method
bution is given by the integral,
We have also used a finite element method (FEM) for purposes of comparison with the results of Pilinis (1990). Equation (1) can be rewritten as
1/2 ¢*
[,,j = | [p,] (~'(j), t) dr' QI
-- 1/2
api
aj, k = y, (k+1)+2~+,((-1)'+1). l
(16)
k=O
The resulting difference in the areas can be accounted for by considering the weighting factor [pi]jIL~, and multiplying it by the [pi] (#'(J), t) values. A sufficient condition for maintaining positivity in pure advection (Bott, 1989a) is given by +
0 <~IL~+,/2 + It.j-,~ 2 ~ [Pi]j
(17)
which maintains a positive flux out of a section and limits it by the value of [Pi]. However, in the condensation/evaporation equation we have an additional term Hip. When the ith chemical component is evaporating from the aerosol, i.e. Hi is negative, the sum of the advective flux and the flux due to evaporation must not exceed the value of [Pi] in the section. If we are to maintain the positivity of the concentration values we need to take into account this decay term, and the positive-definite condition (17) becomes
+
= o.
(23)
Using a weighted residual formulation and considering the Galerkin method for the weight functions, the F E M formulation gives
(A + ~ OAt(B + C)) P~+1 1 =(A--~At(1-O)(B+C))p, - Atl~=l [pl]~D, (24) where A, B, C and Di are matrices whose elements are #2 71
[Ajd= | [N,]j [Nd~dr #1 /z2
[Bjkl = f [Ni]j[Ni]k~;---~lt d# Itl
0 <<.It~j+1/2 + l~Tj-~/2 <~[Pi]j [Hi]j ~ + ~
H 6~pl Pi aH
0-7 - Hip + 3-
#2 i!
LP-b (exp(HAt) - 1).
(18)
O[Ndj [Cj~] = 3| H al~ [Ni]~ dr l*l
Returning to the advective fluxes, we have
F$+,/2 = ~ L~LJ
[pi]j
]l,j+l [Pilj+l
~2 tt
•
(19)
[Djk]l = I Hi[NIl1 [Nl]kdp lal
The flux Fi+ 1/2 depends on the value and sign of Hj at the right boundary. F~+,/2 is set to zero if it is not of the same sign as Hj as it represents the flux due to Hj at that boundary. Finally the flux Fj+ ,/2 can be written as Fj+ 1/2
A# = "~
+ (~l,j+
1/2~01,j - - ~ l _ j + 1 / 2 ~ 1 , j + 1)
(20)
where 31.J+~/2=max( + 0,IL~+ + ~/2), 3#j+. 1/2 = m a x (0, I1~+,/2) and gdz.j= [pi]j/~Lj. If [Hi] is negative, component i is evaporating and
Hi -~ p(p, t) (exp(nAt) - 1) + e)
(21)
where e is a small number introduced to avoid numerical difficulties when ~,~ is nearly zero. If [Hi] is positive, component i is condensing and 3,,~ = max(I~./, 3,+j+ 1/2
+ ~ / ~ j - 1/2 + e).
for j, k = 1, ..., NEN, where N E N is the number of elemental nodes. In each section, we have considered the shape function to be quadratic. The C r a n k Nicolson scheme (0 = 1/2) is used for all cases, which makes the scheme stable and gives a high order of accuracy.
(22)
The limit for 3~.i in equation (22) considers the concentration present in the previous time slice as the limiting factor. An alternate formulation would relax this limit in consideration of the condensational flux.
NUMERICAL TESTS AND ANALYTICAL SOLUTIONS
The three numerical schemes, the Bott-like positivedefinite method, the ASD method and the FEM, are compared for three cases: pure advection, pure growth and simultaneous growth and advection. The accuracy of the numerical schemes is estimated by evaluating the root mean square (r.m.s) errors for each of the three test cases defined as NN
n
r.m.s, error = x/~/=1 ([P,]j - [pi]~,)2
I~p,(u)du
where NN is the number of nodes or sections, [Pi]7 and [p~]~." are numerical and exact solutions at timestep n and section j, respectively. The accuracy of the numerical schemes is a function of two dimensionless groups--the Courant number and the ratio of the
Numerical schemes to model condensation inter-nodal spacing, A#, to the characteristic width of the distribution. Characteristic width is defined as the width of the distribution where the value of p~ is 1/e of its maximum value. Consider the case where H~ is zero, but H is not; that is component i is not evaporating or condensing but one or more of the other components are, such that there is a net growth or shrinkage of the particle. In this case, the governing equation for the mass distribution function of the ith species is given by equation (5). The performance of the scheme can be evaluated using several standard advection tests. The problems most commonly associated with advective schemes are dispersion and diffusion. Using this scheme to advect sine, square and triangular distributions with varying condensation rates will help elucidate the effect of these errors. Figures 1-3 highlight the dispersive errors in the three schemes during the advection of a square wave at large Courant numbers (Courant number = ~ [H]j At/Al~ = 0.99). ASD over- and under-shoots the discontinuities, which leads to negative concentrations. Bott's scheme
923
overshoots near the discontinuities but retains the positivity of the concentration values. Like ASD, FEM suffers from dispersive errors that lead to negative values. These observations are summarized in Fig. 4 where we plot the r.m.s error against Courant number for each scheme. Reducing the timestep reduces the errors in the Bott and FEM scheme. Typically the largest numerical diffusion is produced during the first few advection steps when the discontinuities of the initial distribution are smoothed. In the special case where the net condensation/ evaporation rate of the aerosol, H, is zero, but the condensation or evaporation rate of any individual chemical component, H~, is non-zero, we have the case of pure growth of the p~ distribution, and the governing differential equation is given by equation (4). The p~ distribution changes such that the area under the curve is increasing or decreasing, representing an increase in mass of the ith species, but there is no accompanied movement of the distribution in the size coordinate, as the net particle size is unchanged. Again, the growth/decay of sine, square and triangular distributions are simulated.
1.2 Initial
analytical solution - ASD solution .....
Final
1
analytical solution - -
'~.
1.2 1.4 i 0.8
0.8
0.6
o.,
0.2
=~
. . . . . . . . . . . Li.= . . . . . . . . . =.~,.
-2
' -1
' 0
1
i
,
r
J
2
3
4
5
t
i
61' :g :: j
:~ e -0,2
o
1--li i
0.2 0
solution
~
0.6
0.4
FEM
~ Initial
-02
i~
-0.4 ~ -2
-1
0
1
2
3
4
5
6
7
u
Fig. 1. Pure advection: a comparison of the ASD's predictions with the analytical solution at a Courant number of 0.99. 1.2
Initial 1
analytical solution - B o l t ' s scheme -~....
Final ........ i
Fig. 3. Pure advection: a comparison of the FEM's results with the analytical solution at a Courant number of 0.99. 1.8 1.6
Bolt's scheme ~ .... FEM scheme --~ .... ASD scheme --
",.
1.4 0.8
1.2 1
0.6 0.8 0.4
0.6 0.4
0.2 0.2
-2
-1
0
2
3
4
5
0 0,1
0.2
0.3
0.4
0.5
0.6
0,7
0.8
0.9
1
Courant number
Fig. 2. Pure advection: a comparison of the results obtained using Bott's scheme with the analytical solution at a Courant number of 0.99.
Fig. 4. Pure advection: a plot of the Courant number vs r.m.s, error-accuracy estimation of the schemes for a square wave.
924
S. DHANIYALA and A. S. WEXLER 200
0.005 BOR'S scheme ~--. FEM scheme .-=.... ASD scheme . . . .
0.0045 0.004
150
0.0035 0.003 ~ 0.0025 E ~: 0.002
1O0
0.0015 S0
0.001 0.0005 0
. . . . .
=
-o.oo;
0
;.o~
. . . .
=,
o.o-~5 -o.oi
=
o.o25
0
o.o3
-2
-I
0 Ix
Courant number
Fig. 5. Pure growth: a plot of the Courant number vs r.m.s. error-accuracy estimation of the schemes for a square wave. 200
,
,
1
2
3
Fig. 7. Condensation/evaporation for single species aerosol: a comparison of the results obtained using BoWs scheme with the first analytical solution [equation (25)]; Courant number = 0.1.
analytical solution - . ASD solution ~ - -
150
200
Final
analytical solution - FEM solution • ~....
/~
Final
150 100
I
!ii 100
50
0
............. i
-2
-1
0 ~t
i
i
1
2
Initial 3
0I -2
Fig. 6. Condensation/evaporation for single species aerosol: a comparison of the ASD's predictions with the first analytical solution [equation (25)]; Courant number = 0.1. For pure growth test (H is zero, but Hi is not), a two component system is used such that one compound is condensing on the particle and the other evaporating such that the condensation/evaporation rate of the particle is zero (i.e., H1 = - H 2 ; H =0). The tests show that the schemes model pure growth very accurately, but ASD and FEM may yield negative values, because they do not limit the mass of the evaporating species to the mass initially present. The flux limiters in the positive-definite scheme maintain positivity. For pure growth of a square wave, the r.m.s, errors are plotted against the Courant numbers in Fig. 5, where Courant number for the pure growth is defined as H~At. In the above two cases, pure growth and pure advection have been considered mainly for model evaluation. Generally in atmospherically relevant processes, the growth term and the advection-like term are both non-zero. Analytical solutions for pi are not available but solutions for p = ~ i pi have been derived: p(#, t) = 3 exp(6g) e x p ( - exp (3#) e x p ( - Ht) - Ht) (25) p(p, t) = exp(3/z)exp(- e x p ( - 3 g ) e x p ( - Ht)).
(26)
:=~==__________=====:
i
i
i
i
-1
0
1
2
3
Fig. 8. Condensation/evaporation for single species aerosol: a comparison of the FEM's results with the first analytical solution [equation (25)]; Courant number = 0.1.
Bolt's scheme
FEM scheme --~.... ASD scheme .....
// 1.5
/)
/-
.................%11 ........................ ..J-.o-_ ......... . ......... . 0.5
0
0
0-1
012
0.3
0.4
0.5
0.6
0.7
0.S
0.g
Courant number
Fig. 9. Condensation/evaporation for single species aerosol: a plot of the Courant number vs r.m.s, error-accuracy estimation of the schemes for the first analytical solution [equation (25)]. Equation (25) was derived by Gelbard and Seinfeld (1980), while equation (26) was derived for use in this work. For a single component aerosol undergoing
Numerical schemes to model condensation 0.12
925
140 .....
Bott ASD ...... FEM ..6....
0.1
Final
120
anal'flical solution - Bott's scheme -+....
100
0.08
80 0.06 ,
iT60
I 0.04
40 0.02
ol
20 :
0.15
...........
0.2
I
~
0.25
0.3
_ i
,
0.35
0.4
0.45
0,5
0.55 -2
-1
0
inter-nodal spacing/characteristic width
1
2
3
F
Fig. 10. Condensation/evaporation for single species aerosol: a plot of the ratio of the characteristic width to section size (A~) vs r.m.s, error-accuracy estimation of the schemes for the first analytical solution [equation (25)].
Fig. 12. Condensation/evaporation for single species aerosol: a comparison of the results obtained using Bott's scheme with the second analytical solution [equation (26)]; Courant number = 0.l. 140
14o [----
,
analytical solution
'
120 ~
-ASD solution . . . .
Final / ~
,o0
Final .
120
analytical solution - FEM solution t~ .
100
/1
80 60 40 20
20
[ni '
_ ol ...... -2
,
-1
,
0
...
-i
1
:~±
i,
0,=-
.... i ......... i ......... l
-20 -2
-1
0 u
2
1
2
3
la
Fig. 11. Condensation/evaporation for single species aerosol: a comparison of the ASD's predictions with the second analytical solution [equation (26)]; Courant number = 0.1.
Fig. 13. Condensation/evaporation for single species aerosol: a comparison of the FEM's results with the second analytical solution [equation (26)]; Courant number = 0.1. 1.2
condensation, results of the numerical predictions for the first analytical solution are illustrated in Figs 6 - 8 and for the second in Figs 11-13. The results shown are for test cases with Courant number of 0.1 (Courant number = ½ [ H ] j At/A/~). Figures 9 and 10 show the r.m.s error dependence on Courant number and internodal spacing, respectively, for the first analytical solution [equation (25)], while Figs 14 and 15 show the same for the second analytical solution [equation (26)]. Since mass distributions given by the two analytical solutions are Coo smooth functions, A S D models the condensation process for this initial distribution very accurately. However, time-splitting and time-stepping errors associated with the positivedefinite scheme and FEM, respectively, reduce their accuracy. Although analytical solutions are unavailable for the distribution of the individual species in a multicomponent aerosol, the total distribution of all the species of the aerosol can be compared with analytical
~o~,,sc.=na- ..... t FEM scheme .-,.--. J ASD scheme -*- ;/t
1
0.8
................ii / ...........
0,6 J
j,,"
0.4
..........i -....... ......... ~
t
.-
m-..0,2
0
0~;"-';.~"-;~;--0.~ o~5-0~.~-'";.V-o~.i ;.8 Courant number
Fig. 14. Condensation/evaporation for single species aerosol: a plot of the Courant number vs r.m.s, error-accuracy estimation of the schemes for the second analytical solution I-equation (26)]. p distributions. The initial distribution of the three species are considered identical and similar to the distribution obtained from the first analytical
926
S. DHANIYALA and A. S. WEXLER 0.05
,
.
'
eott
0.045
! i
ASD - 4 . . . . FEM --~....
0.04 0.035
~ ~
'
h_,,.o~ ~ 8 o t t t o k a ~ .Q--
hal
/
:::SSI......
0.03 0.025
..,.
0.02
.:.-"
,,
;:/ 0
........./," "/
...a ........ ..._.4-
0.015
5
0.01 0.005
~
5 0
0.05
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0.55
inter-nodal spacing/characteristicwidth
t
-5 -2
i
-1
i
o
t
1
2
i.t
Fig. 15. Condensation/evaporation for single species aerosol: a plot of the ratio of the characteristic width to section size (A#) vs. r.m.s, error-accuracy estimation of the schemes for the second analytical solution [equation (26)].
i
'
. -
Fig. 18. Multispecies aerosols: plot of growth of the mass distribution of the third component accompanied by an increase in particle size (H3 = 10; H = 9) as modeled by the three schemes for the initial distribution of the species given by the first analytical solution [equation (25)].
Bott's scheme FEM scheme --=.... ASD scheme ......
1
Fired
0.8
/ #/" /,
0.6
/,,~'-..
0.4
..**'/"" .i
i
-2
-1
1
........ i .........
0.2
2
~nitial
Final
FEM md~Jon -*,-aottm=~n .=..
g
Fig. 17. Multispecies aerosols: plot of the movement of the mass distribution of the second component due to the increase in particle size (H2 = 0; H = 9) as modeled by the three schemes for the initial distribution of the species given by the first analytical solution [equation (25)].
. .e'""'"
-
....,"
"'",0""""
..... ..~r'"
e'"'""
0
Fig. 16. Multispecies aerosols: plot of the decay of the mass distribution of the first component accompanied by particle growth (Ht = - 1; H = 9) as modeled by the three schemes for initial distribution of the species given by the first analytical solution [equation (25)].
//
ml
,./
01 0.2 03 0.4 0.5 0.6 0.7 o.s 0.9 Courant number
Fig. 19. Multispecies aerosols: plot of the numerical errors with Courant number for the three schemes for the first analytical solution [equation (25)].
solqtion. The first species is considered to evaporate with a rate of H t = - 1 , the second species neither evaporates nor condenses, i.e. H2 = 0, and the third species is assumed to grow at a constant condensation rate of H3 = 10. The total condensation/evaporation rate is H = 9. The change of the mass distributions of the three species due to the different condensation and evaporation rates is shown in Figs 16-18. The r.m.s, error dependence on the Courant number for the total mass distribution (p), shown in Figs 19 and 20 indicates that the positive-definite scheme models the multicomponent condensation/evaporation process accurately at small Courant numbers. The runtimes of the schemes for modeling condensation/evaporation of a three-component species whose total mass distribution p is given by the first analytical solution [equation (25)], is shown in Fig. 21. The positive-definite scheme is seen to be the least
Numerical schemes to model condensation 0.7 BoWs scheme
FEM scheme ASD scheme
0.6 0.5
,/"
0.4 0.3 7
a
....
~r
0.2 0.1
927
into Eulerian air pollution models to predict the size and composition distribution of atmospheric aerosols. The Bott's scheme used in this work has been modified by several researchers to improve the scheme's stability and monotonicity. We have, however, considered the original version of Bott's scheme and proven its application to obtain a positive-definite scheme to model condensation/" evaporation in aerosols. For some applications, incorporating these improvements to Bott's scheme will improve performance.
0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Courant number
Fig. 20. Multispecies aerosols: plot of the numerical errors with Courant number for the three schemes for the second analytical solution [equation (26)1.
18 BoWs scheme ASD scheme ---.... FEM scheme .......
16 14 12 k 10 8 6 4
........... e ...............
2 0 0.1
~ 0.2
-,¢----,--._,_;...........,,-.........-_,___;........ 0.3
0.4
0.5
0.6
0.7
0.8
0.9
Courant number
Fig. 21. Runtimes: a comparison of the time taken by the three schemes for advecting a multicomponent species, whose total mass distribution function p is given by the first analytical solution [equation (25)], through one cell width for varying Courant numbers.
cpu intensive. All the tests were conducted on a S U N 4/670MP workstation.
CONCLUSIONS Three numerical methods were investigated for solving the internally mixed condensation/evaporation equation. Tests were performed on the ability of the methods to represent growth and decay of the size distributions, their shift in the size coordinate, and both processes simultaneously. ASD, a spectral method, is shown to be the slowest and generally the most accurate, but it suffers from Gibbs errors and often predicts unphysical negative concentrations. The finite element method suffered from large dispersive errors for all but the smoothest distributions. A new method developed here and based on Bott's advection scheme is relatively fast and accurate especially at low Courant number, and is positive-definite. This method will be incorporated
REFERENCES
Bassett M., Gelbard F. and Seinfeld J. H. (1981) Mathematical model for multi-component aerosol formation and growth in plumes. Atmospheric Environment 15, 2395-2406. Bott A. (1989a) A positive definite advection scheme obtained by nonlinear renormalization of the advective fluxes. Mon. Weath. Rev. 117, 1006-1015. Bott A. (1989b) Reply. Mort. Weath. Rev. 117, 2633-2636. Chock (1991) A comparison of numerical methods for solving the advection equation-II. Atmospheric Environment 25A, 853-871. Chock D. P. (1985) A comparison of numerical methods for solving the advection equation-II. Atmospheric Environment 19, 571-586. Chock and Dunker (1983) A comparison of numerical methods for solving the advection equation-lL Atmospheric Environment 17, 11-24. Crowley W. P. (1968) Numerical advection experiments. Mon. Weath. Rev. 96, 1-11. Gazdag J. (1973) Numerical convective schemes based on accurate computation of space derivatives. J. Comput. Phys. 13, 100-113. Gelbard F. (1990) Modeling multicomponent aerosol particle growth by vapor condensation. Aerosol Sci. Technol. 12, 399-412. Gelbard F. and Seinfeld J. H. (1980) Simulation of multicomponent aerosol dynamics. J. Colloid Interface Sci. 78, 485-501. Gelbard F., Tambour Y. and Seinfeld J. H. (1980) Sectional representations for simulating aerosol dynamics. J. Colloid Interface Sci. 76, 541-556. Kim Y. P. and Seinfeld J. H. (1990) Simulation of multicomponent aerosol condensation by the moving sectional method. J. Colloid Interface Sci. 135, 185-199. MacCormack R. W. (1969) The effect of viscosity in hypervelocity impact cratering. A I A A Pap. 69-354. Oran S. E. and Jay Boris P. (1987) Numerical Simulation of Reactive Flow. Elsevier, Amsterdam. Pilinis C. (1990) Derivation and numerical solution of the species mass distribution equations for multicomponent particulate systems. Atmospheric Environment 24A(7), 1923-1928. Roache P. J. (1976) Computational Fluid Dynamics. Hermosa. Rood, R. B. (1987) Numerical advection algorithms and their role in atmospheric transport and chemistry models. Rev. Geophys. 25(1), 71 100. Smolarkiewicz P. K. (1983) A simple positive-definite advection scheme with small implicit diffusion. Mon. Weath. Rev. 111, 479-486. Smolarkiewicz P. K. (1984) A fully multidimensional positive definite advection transport algorithm with small implicit diffusion. J. Comput. Phys. 54, 325-362. Tremback C. J., Powell J., William R. C. and Pielke A. R. (1987) The forward-in-time upstream advection scheme:
928
S. DHANIYALA and A. S. WEXLER
extension to higher orders. Mon. Weath. Rev. 115, 540-555. Tsang T. H. and Rao A. (1988) Comparison of different numerical schemes for condensational growth of aerosols. Aerosol Sci. Technol. 9, 271-277. Warren D. R. and SeinfeldJ. H. (1985) Simulation of aerosol size distributionevolution in systems with simultaneous
nucleation, condensation and coagulation. Aerosol Sci. Technol. 4, 31-43. Wexler A. S., Lurmann F. W. and Seinfeld J. H. (1994) Modeling urban and regional aerosols-I. Model development. Atmospheric Environment 28(3), 531-546. Yanenko N. N. (1971) The Method of Fractional Steps. Springer, Berlin.