SpectrochzmicaAcm, Vol 399, No. 8, pp 1005-1010.19W Printed m Great Bntain.
0584-8547184 SO3.00+ .OO
0 1984. Pergamon Press Ltd.
A method for the solution of the mass transport equation in a free burning d.c. arc M. PERIL and J. RADIC-PERIL University of Belgrade,Faculty of Sciences, Institute of Physical Chemistry,P.O. Box 550, 11001, Belgrade, Yugoslavia (Received 21 February 1984; in revised form 4 May 1984)
Abstract-A method is developed for the solution of a generalform of the mass transportequation for a free burningd.c.arc.The spatialparticledensitydistributionfunctionis representedin form of a linear combination of appropriate basis functions in which the expansion coefficients are calculated using the variation principle. The method allows calculation of particle distribution at various levels of approximation, and also includes cases in which both radial and axial dependences of the diffusion coefficient as well as the directed transport velocity are taken into account. The results of some simple test calculations are presented.
1. INTRODUCTION THE STUDYof the mass transport in excitation sources is motivated by the fact that the intensity of spectral lines and bands is dependent on the mean time the particles of the investigated substance spend in various regions of the discharge zone. There are a number of papers, e.g. Refs [l-5], devoted to the combined theoretical and experimental investigation of transport phenomena in the plasma of a free burning vertical d.c. arc. BOUMANS[l, 23 developed a mass transport model based on experimental findings of DE GALAN [3]. The diffusion and directed axial velocity caused by the convection and the presence of an electric field were considered, the arc being treated as an infinitely long cylindar with a homogeneous core. The assumption was made that the outward diffusion of the substance from the core was precluded and the diffusion coefficient as well as the axial velocity were taken to have constant values in the core. The spatial distribution of the particle density, obtained by solving the corresponding differential equation, could describe the main features of the experimental facts. The model of BOUMANSwas improved by KRINBERGand SMIRNOVA[4]. They took into account the change in the velocity along the arc axis and investigated the effect of the upper electrode as a mechanical obstacle, as well as the effects of carriers on spatial distribution. The influence of the cathode as an electric obstacle on the axial concentration distribution was investigated by HOLCLAJTNER-ANTUNOVIC et al. [S]: All the mentioned methods, though successful in interpreting the main experimental facts, have two serious drawbacks: Firstly, due to some inherent assumption, e.g. D = const., they are restricted to the description of the phenomena in the central regions of the arc, and thus the calculation of the spatial particle density distribution in the outer regions-in which many of the spectral lines and most of the molecular bands are excited-is less accurate. Secondly replacement of any of the assumptions or of the boundary conditions by others requires a new approach for the solution of the mass transport equation. Therefore we want to develop an approach which makes a reliable treatment of the transport phenomena in the colder arc zones possible and that is general enough to enable calculations at various levels of approximation. In the present paper this method, based on an application of the variation principle, is described and illustrated by some model calculations. In the next paper [6] the experimentally obtained [l] [2] [3] [4]
P. W. J. M. BOIJMANS,Theory of Spectrochemical Excitation. Hilger and Watts, London (1966). P. W. J. M. BOUMANSand L. DE GALAN, Anal. Chem. 38,674 (1966). L. DE GALAN, J. Quant. Spectrosc. Radiative Transfer 5, 735 (1965). I. A. KRINBERGand E. V. SMIRNOVA,Zh. Prikl. Spektrosk. 10,400, (1969);13, 5 (1970);16, 17 (1972);18,970
(1973). [5] 1. HOLCLAJTNER-ANTUNOVI~,V. JADE?, V. VUKANOVICand V. GEORGIJEVIC,Spectrochim. Acta 378, 889 (1982). [6] J. RADI&PERIC, to be published. 1005
1006
M.
F%RI(Sand
J.
R,uW-but
spatial distributions of diffusion coefficients and axial transport velocities are employed to calculate particle density distributions which are then used to interpret measured spectral line and band intensities. 2. METHOD FOR SOLUTION OF THE MASSTRANSPORT EQUATION 2.1. Variational method The variational (Ritz) principle is widely applied in quantum mechanics where it represents a basis for the most powerful method for an approximate solution of the Schriidinger equation. However, it seems that its application to problems of classical physics has not been so common in spite of the fact that the mathematical structure of the equations to be solved is often very similar to that of the Schriidinger equation. In the present paper a method is developed for the solution of the mass transport equation of a system-a vertical free burning d.c. arc-characterized by an axial flow of the matter due to the presence of external forces (convection and electric field) and a simultaneous lateral spreading of the substance due to diffusion. The method is based on Galerkin’s variance of the variational calculus [7,8] in which the solution, i.e. the spatial distribution of the concentration of the substance evaporating from the lower electrode, is developed into a series of appropriate basis functions. A variational problem can be defined in the following way: let it be desired to find the function y(x) which makes the value of the integral b
UYI =
s0
F 6, Y, Y’) dx
(1)
to be stationary, in which some boundary conditions must be satisfied. It can be proved [7,8] that the function y is the solution of the Euler equation, having in this case the form
F&F;<=0.
(2)
According to Ritz’s method an approximate solution of the problem (1) or/and (2) can be found by developing y in a series (3) where g&) represents a set of basis functions and the coefficients ck are chosen to minimize the value of the functional Z(y). Because of the boundary conditions which have to be satisfied, there are generally some constraints on the values of these coefficients-they are not all linearly independent of one another. The Ritz’s method can be used to solve the differential equations since the differential equation to be solved can be treated as an Euler equation of a functional. However, in order to apply this method for the solution of differential equations it is not necessary to seek the corresponding functional. It is assumed that the variational problem can be written in the form (Galerkin’s method): 6Z {Y, 6Y} = (L Y96Y)
(4)
where the parenthesis on the right side denotes the scalar product, and L represents some operator, in our case Ly = F;--&F;s
(5)
then, with y given by (3), the Euler equation becomes equivalent to Lz j =
1, 2, . . . n.
{cckgk}
=
(Lcckgks
gj)
=
0
(6)
I
[7] H. MARGENAUand G. M. MURPHY, The Mathematics ofPhysics and Chemistry. D. Van Nostrand, New York (1956). [8] G. STRANG and G. J. FIX, An Analysis ofthe Finite Elements Method. Prentice-Hall, Englewood Cliffs, N.J. (1973).
Solution of the mass transport equation
1007
Equation (6) shows that the required differential equation can be solved by projecting it onto each of the basis functions and solving the system of algebraic equations so obtained. 2.2. Mass transport equation The general form of the diffusion equation is an - = V(DVn) = DAn + VDVn
at
(7)
with n and D representing the concentration and diffusion coefficient of the substance respectively. For a description of the mass transport in a free burning d.c. arc, it is useful to transform (7) into cylindrical coordinates. Assuming axial symmetry, one obtains: a2n aDan ;$+DG+asaz=O.
Usually the assumption is made that D has a constant value, with the consequence that the last term in Eqn (7) and all terms including derivatives of D in (8) are equal zero. However, this assumption is not reliable in an arc plasma because at least there is a strong radial dependence of temperature and consequently of D. Besides diffusion the particles in a free burning d.c. arc are subjected to the influence of the axial electric field and convection forcing them to move upwards with a velocity u which depends, in principle, on both the r and z coordinates. Some other effects, e.g. thermal diffusion, will not be considered in the present paper. Therefore, the steady mass transport equation we want to be solved is:* (9) We assume that the boundary conditions this equation should satisfy are n=Oif
z=oO r= co
(10.1)
and -Dg+vn
=f(r)
for 2 = 0.
(10.2)
The first two conditions are trivial. The second is equivalent to that used by BOUMANS [ 1,2] to regulate the flux of the substance from the lower electrode-the only difference is that we introduce the functionf(r) being a generalization of the Heaviside step function employed by BOUMANS.The third boundary condition introduced in [l, 23, implying the existence of homogeneous arc core with impenetrable walls or, equivalently, assuming a rectangular profile of the diffusion coefficient is rejected in the present paper: instead the radial variation of D is taken explicitly into account. 2.3. Method for the calculation of the spatial particle distribution If the usual assumption D = const. was made it would be possible to separate the variables in Eqn (9) assuming the solution n(r, z) to be of the form R(r)Z(z). Such a separation into two differential equations, one depending on r and the other on z only, can not be made in our case. We seek an approximate function n in the form
(11) withf, and c#+,, denoting the basis functions depending on the r and z coordinates respectively.
*If the study of the mass transport is not restricted to the central zones of the arc, the radii compqnent of the directed velocity u can become significant and the terms - (l/r) [a(w)/&] - u(an/&) should be added to the left side of Eqn (9).
1008
hf. PERICS and J. RADI&PERIL
The functionf(r)
in (10.2) can be represented in the form f(r) = f
(12)
d._L(r)
n=1 so
that the boundary condition (10.2) reduces to
i.e. due to the linear independence of allf, functions:
IZ= 1, 2, . . . N with z=o
(13)
-e#4n)z=o.
In (13) it is assumed that the diffusion coefficient has a constant value over the whole top of the crater in the lower electrode (z = 0) through which the substance comes into the plasma. Equation (13) represent N constraints which couple the coefficients c,,. They can be taken into account applying the method of Lagrange multipliers-instead of 61(y) one seeks for variation of the expression (14) Equating the expression (14), i.e. all partial derivatives a/&i, { } with zero, a system of N x M equations is obtained F
5
C,,Lij,nm+Aihj=O
i= 1,. . . N,
j =
1,. . . M
n=lm=l
(15)
in which L,.,
E (_MjJLl_L&) E
a fS0
mX(r) ~j(z)Lf,(r)~,(z)rdrdz. 0
(16)
Equations (13) + (15) represents a system of N(M + 1) linear equations with the N x M unknown coefficients c,, and N Lagrange multipliers rZi. 2.4. Basis functions As basis functions j,,(r) we choose the normalised products of an exponential function exp ( - 1/21r2), with L being a parameter which should be optimized for a given problem and the confluent series ,F, (a, 1, lr2) [9], (a = 0, - 1, -2, . . . ). f,(r) = N, 1F, (a, 1, Ir2) exp (- 1/2Jr2) lF,(a,
’ a(a + 1) A2r4 1, Ir’) = 1 +y$+r~+
a(a + 1) (a + 2) L3r6
1.2.3
N,=fi.
where
(17) 3I + .**
(18) (19)
These basis functions automatically ensure a fulfilment of the condition (lO.l), n(r = 00) = 0. The confluent series have the advantage over simple polynomials of r in that they are orthogonal to one another (d V = r dr) thereby resulting in a faster convergence of the trial function n(r, z) to the exact solution. Basis functions 4,(z) are taken in the form $J,&) = NJ&
(,&)
exp (- l/2/1r2)
[9] S. FLUGGE,Practical Quantum Mechanics. Springer, Berlin (1971).
(20)
Solution of the mass transport equation with I$,,(&)
1009
representing Hermite polynomials [7,9] and
(21) the sub-basis {&,,(z)} also being orthogonal (0 < z < cc) and ensuring the proper asymptotic behaviour of n (n, = m= 0). The diffusion coefficient depends in principle on both r and z coordinates but because of the fact that in a vertical free burning d.c. arc the axial temperature gradient is much less pronounced than the radial gradient and consequently aD/ar is usually appreciably larger than aD/az we neglect in actual calculations the last term entirely. The diffusion coefficient is represented in form of a polynomial in rz with an exponential function exp( - 1/2yrZ). The exponent and the coefficients in that polynomial expansion are chosen to ensure a good fit of the curve D(r) which in principle is experimentally obtained. The radial distribution of the flux of the evaporating substances at z = 0, f(r), has to be represented as a linear combination of basis functions f,(r). Axial dependence of the axial velocity u is approximated by a polynomial in z. Such an expansion should be flexible enough to take into account the observed variation of u along the arc axis [lo, 111. For a better description of the region in the vicinity of the upper electrode, where the value of the axial velocity can considerably decrease [S], u can be expressed in the form of a polynomial with an appropriate exponential function. The radial dependence of u can be represented in the same way as for D, the most general form of u being the product of both zand r-dependent functions. With the basis functions and the representation of D and u chosen, all the integrals (16) are reduced to simple sums of I functions. 3. RESULTSOF MODELCALCULATIONS AND DISCUSSION In order to illustrate how the proposed method works we present in this chapter the results of some simple model calculations. Since there is no need to insist on very accurate numerical calculations, only three basis functions for r and three for z-coordinates have been used. The complete basis consists in this way of nine (3 x 3) functions due to the three constraints (Eqn (13)) only six of them are linearly independent. The spatial particle density distributions are calculated for various forms of the particle’s source, the diffusion coefficient and the axial velocity. The results are presented in Fig. 1. Though only of a qualitative character, the results obtained give a picture of the influence of the mentioned parameters on the particle distribution. First of all, as already stated in previous experimental and theoretical studies [l, 31, there exists generally a decrease of the particle concentration along the z-axis and a flattening of the radial distribution with increasing distance from the lower electrode. The diameter of the crater of the lower electrode, acting as a particle’s source, influences strongly the shape of the distribution curves as is evident by comparing the upper and lower parts of Fig. 1. Secondly, using a smaller value for the diffusion coefficient e.g. 5 cm2 s- ’ compared to the 30 cm2 s- ’ used generally in these test calculations, results in a more rapid radial decrease in the particle density. The same effect, but even more pronounced, is caused by the use of a radially dependent diffusion coefficient represented by a strongly falling curve (D = 30 exp (- 6r’)). An increase in the value of the axial velocity causes a proportional decrease of the concentration along the arc axis. A change in the value of the diffusion coefficient does not have such a direct influence on the axial distribution but, as already stated by BOUMANS [l, 21, if by increasing the D value, the ratio u/2D is kept constant, the absolute concentration diminishes by the same factor. Finally, if the axial velocity is taken to be z-dependent, the axial particles distribution is perturbed compared to the case of a constant u-value-an increase of the axial velocity with increasing distance from the lower electrode causes a relative decrease in the particle concentration, while a decrease of u with increasing z can, on the other hand, result in an enrichment of the investigated substance in corresponding regions. [lo] W. HAGENACH, Z. Physik128, 278 (1950). [11] P. S. TOD~ROVI~,V. M. VUKANOVI~,M. M. SIMIC? and M. N. PERIC,Spectrochim. Acta 31B, 103 (1976).
1010
M. PEW and J. RADICPERIC
nrel 2
10
I D (cm2s-1)
1
30
---_-_.
5 30 60
~#DODDO 30er$6r2) _---. 30
I
v(cms-1) 200 200 600 400 200+300 z 200
I
200-100.1
hlvn)
rkm)
I (up)
Fig. 1. Spatial distributions of the particles’ density calculated for various forms of D, u and f(r). Top: a moderately wide disc source (j(r) = exp( -r*)). Bottom: a very narrow particles’ source (f(r) = exp (- 6r*)).
The results of these test calculations indicate that the method proposed in the present paper can be used for the description not only of the coarse structure of the spatial particle distribution in an arc plasma, but also of the finer details caused by the existence of inhomogenities in the system studied. The results presented in this chapter were obtained by the solution of a system of twelve linear algebraic equations with twelve unknowns (nine coefficients in the expansion of particles density distribution function and three Lagrange multipliers). A real calculation with more complicated forms of D, u andf(r) would necessarily require a larger functional basis and consequently the order of the system of equations to be solved would be greater. However, the increase in the number of elementary algebraic operations is the only price that must be paid when considering and describing finer effects.