Numerical simulation of lonization instability

Numerical simulation of lonization instability

Energy Conversion. Vol. 10, pp. 1-5. Pergamon Press, 1970. Printed in Great Britain Numerical Simulation of Ionization Instability LAJOSL, LENGYELt...

526KB Sizes 0 Downloads 81 Views

Energy Conversion.

Vol. 10, pp. 1-5.

Pergamon Press, 1970. Printed in Great Britain

Numerical Simulation of Ionization Instability LAJOSL, LENGYELt (Received

20 May 1969)

Introduction

Ionization- or electrothermal-instability [1, 2] is a phenomenon observed in low temperature nonequilibrium plasmas used in certain MHD energy conversion experiments. Under the action of a sufficiently strong magnetic field perpendicular to the direction of the applied (or induced) electric field an initially homogeneous plasma structure may undergo complete breakdown leading to the formation of striae consisting of alternating layers of higher and lower current density, electron temperature, electron density, etc. [3, 4]. The effective plasma parameters (electrical conductivity, Hall parameter) may be significantly impaired by the striated structure. The existing perturbation analyses (for a comprehensive review of the literature available on ionization instability the reader is referred to Refs. [5-7]) have rendered valuable information on the stability boundaries, initial growth rates, the effects of various factors on the development of instabilities. However, the applicability of these results is restricted to the region of small disturbances. The information gained from nonlinear analyses is more complete. However, because of the complexity of the problem the introduction of various simplifying assumptions (infinite plasma, plane waves, etc.) becomes unavoidable also in these analyses. This again imposes certain limitations on the applicability of the corresponding results. The purpose of the present study is to follow up numerically the development of electrothermal waves in a domain where the disturbance can no longer be considered as a small one. The analysis is applied to a finite plasma bounded by physical boundaries. Mathematical Model

A brief description of the mathematical model considered here has been given in Ref. [8]. A similar, but somewhat simpler, model has also been used in Ref. [9]. The following case shall be treated here: A plasma is enclosed by a vessel consisting of two infinitely segmented and ideally emitting electrode walls at y = 0, H (Jtan = 0) and by two insulator walls at x = 0, S (Jnorm ~ 0, see Fig. 1). The total current sent through the plasma in the - - y direction is given and kept constant, the current distribution is maintained homogeneous at the electrode walls by connecting each electrode pair through a separate circuit. A homogeneous magnetic field B is applied in the ÷ z direction, all variations in the direction

Y H

S Fig. 1. Assumed geometry and the unperturbed current distribution.

of the magnetic field are neglected, the vessel is assumed to be of unit depth. One will note that the restrictions imposed on the current distribution by the presence of the insulator walls at x = O, S could easily be removed by assuming a sufficiently long duet length in the x direction and a periodic distribution for the section 0 ~ x ~ S considered here. Crucial, however, is the selection of the boundary conditions at y = 0 and y = H: by assuming finite length electrodes [10], for example, the boundary conditions become magnetic field dependent and the current distribution along the boundaries strongly non-uniform. The instabilities developing in such a system are rather difficult to separate from those of a purely numerical nature caused by large geometry-dependent field gradients. In the usual two-temperature plasma approximation, and neglecting radiation, thermal conduction, and like processes, the current distribution is defined for a potassium seeded argon plasma by the following set of equations: J = aE--

V.

×

J = O,

+

Vpe

(1) (2)

eB ,

mere

g

V x E ~- 0

Hee 2 --

t Institut far Plasmaphysik, Garching bei Miinchen, Germany.

B

=

(3)

DleVe

2

L A J O S L. L E N G Y E L

One

at

2 -- ne(neeql

2 v - - h e ) reeomb,

(4)

2

?/eeql ~ fsaha(Te) ?/KO - - ?/eeql

Ot ne

q- E~

= - - -- 3k - - neve(Te -- Tg) a ma q- ~-'e [V(~kTe)--kTeVneJne

4 18kTe\

i

/ - - ~ I ( n A Q e A q- nKQeK q- n,,+Q.),

(5)

(6)

where the usual notation has been used. The collision cross sections QeA, Qe~:, and Qa as well as the recombination rate coefficient Vreeomb are considered as functions of the electron temperature [1 l, 12]. A stream function is introduced for the electric current by means of the relations Jx = aT/3y, Ju = - ~y/ax; the twodimensional current distribution is given by the solution of the nonlinear partial differential equation

its time development. This is precisely the procedure used in the present computations: at time t = 0 a localized perturbation is introduced into the central region of an initially homogeneous field which may be stable or metastable, depending on the applied magnetic field strength. The time history of the disturbance was computed by solving iteratively the above set of equations in a manner similar to the method described in Ref. [13]. Special attention has been paid to the accurate determination of the electron temperature, for this reason Equations (4)-(6) were iterated simultaneously (an approximate integral of Equation (4) is available [13], but the quantities neeql-~-neSaha, and Vreeomb are still functions of the electron temperature). Results and Conclusions

Typical results of numerical computations are given in Figs. 2-6. The input data used in these calculations are as follows: H = S = 7 cm, Jy0 = -- 1 A/cm 2, potassium seeded argon plasma at p ---- 1 atm, Tg = 1800°K, CK = 0.002, subscript 0 denotes initial conditions. At time t = 0 the electron temperature was perturbed in the middle of the field (either in a single point, or in a2~ q- 02y 0~, 0v ~x~2 Oy~ q- P(x, y) ~x q- Q(x, y) ~y = R(x, y), (7) the five central points of a network system consisting of 29 × 29 mesh points) by means of a localized energy input (or extraction) pulse (ATe/Te ~ 0"01). The time with the boundary conditions development of the disturbance was computed with At = 10 t~sec time increment. As was expected, at low ~ Y - - 0 at x = 0, S; y = 0, H; ay magnetic field values the disturbance disappeared from the field within the first 50-80 ~sec. and (8) Machine plots showing the evolution of the current distribution from initially homogeneous metastable -- const at y = O, H. states ale given for B ---- 2 tesla (/3o ~ 11.5) in Figs. 2 ax and 3. Figure 2 corresponds to a positive central Te The coefficients P( x, y), Q(x, y), and R(x, y) in Equation perturbation, Fig. 3 to a negative one. The formation (7) are functions of the spatial derivatives of the electrical of striae consisting of alternating layers with higher and conductivity ~, Hall parameter /3, electron temperature lower current density, electron temperature, and electron density is dearly visible. Instead of diffusing homoTe and electron density ne [13]. In this formulation the magnetic field affects only the geneously, the disturbance develops in well defined current distribution itself [through the t-derivatives in preferential directions, its strength increases with time. P(x, y) and Q(x, y)] but not the boundary conditions. At the beginning, this development remains completely The effects of the boundaries on the electron gas (wall unaffected by the presence of boundaries. The current recombination, etc.) are neglected. The steady-state flows in a characteristic zigzag pattern. For the case solution that satisfies the above equations and boundary shown in Fig. 2, the development of the electron temconditions is 07/Oy = Jx ~ 0 in the entire field. From perature distribution along the vertical axis x ---- S/2 is shown in Fig. 4, while the electron temperature and the condition V • J = 0 it follows that electron density evolutions shown in Fig. 5 correspond to B = 1.5 tesla and a negative temperature perturbaJ = Jy -----const = Iu/S, tion at t-----0. As can be seen, for the initial electron where Iu is the total current sent through the plasma. density considered here the development of the striae is The current stream lines are parallel straight lines (see relatively slow in the first phase of their formation, then Fig. 1). This solution is independent of the magnetic the rate increases rapidly. The time necessary for the field strength. Hence if the assumed ideal conditions formation of the striae seems to be governed here by the could be realized in practice, one could hypothetically electron density relaxation time. The developing waveincrease (infinitely slowly) the magnetic field strength to length (2.0-2.5 cm) is not related to the size of the any value without producing any change in an initially central perturbation introduced into the field. The uniform current distribution. Whether or not the final extrema of the curves shown in Figs. 4 and 5 shift with state is stable (or metastable) can be checked by intro- time in the direction of the current. The propagational ducing a perturbation into the field and following up velocity of the waves calculated from the magnitude

Numerical Simulation of Ionization Instability

(a)

3

(a)

/!l (b)

(b)

(c) Fig. 2. Time development of the current distribution for B = 2 tesla, Juo = 1 A l c m ~, and (/lTe/Te)o ~-0"01. (a) t = 40 /~sec; Co) t = 60 psec; (c) t = 80 psec.

Fig. 3. Time development of the current distribution for B = 2 tesla, Jyo = -- 1 A/cm2, and (ATe/Te)o ~- -- 0"01. (a) t = 60 ~sec; Co) t = 80/~sec; (c) t = 90 ~ c .

4

LAJOS L. LENGYEL

of this shift is about 10-40 m/sec. A comparison of the results corresponding to different magnetic field strengths shows that the wave-length decreases and the rate of growth of the disturbance increases with increasing magnetic field strength. Figure 6 shows the development of the striae for a case where a uniformly distributed perturbation was introduced into the field. An interesting aspect of this pattern is that the current stream-lines become constricted with time also in the cross-flow layers of initially homogeneous low current densities. The resulting cellular pattern is in agreement with certain experimental observations. The wavelength seen here is defined by the periodicity of the applied perturbation. The striae propagate in the negative y direction with a velocity of approximately 10 m/sec. Te

[°K]

(a)

j J S O jJsec

2700

i

.. ~BOpsec

50jusec ~

v-,-".-~0 H sec

2600

/0/~sec

y

/ 2500

2420 ' 0

/

J,, I

2

/ 4

3

5

6 y[cm]

Fig. 4. Time developmentof the electron temperature distribution along t h e x = S / 2 axis for the ease shown in Fig. 2.

/

\,

120 psec

(hi

[o~] 2700

~

~

2600 25OO

2t00

n. x'lO-~,

~J120

/

1.5

0,5

i

i

t

I

2

3

i

/_~

pse¢ ~100 psec

i

t

4

s

(c)

i

y Ccm]

Fig. 5. Time development of the electron temperature and electron density distributions along the x = S / 2 axis for B = 1"5 tesla, Jyo = -- I A/era2, and ( A T e / T e ) o = - - 0"01.

Fig. 6. Time development of the current distribution in the case of a periodic perturbation pattern. B = 2 tesla, Jyo = -- 1 A l c m % ( A T d T e ) o = - - 0"01. (a) t = 75 Fsee; (b) t = 105 t~see; (c) t = 120/~sec.

Numerical Simulation of Ionization Instability

In each case described above, the computations were terminated when the developing large field gradients in combination with the given finite mesh size began to impair the accuracy of the calculations (see, for example, Figs. 2(c) and 3(c). At this point, the electron density fluctuation has reached approximately 100 per cent (ne/neo ~ 2, see Fig. 5). It is noteworthy that some of the final random configurations reported in Ref. [9] correspond to ne/neo = 7 by a mesh system consisting of only 200 (10 × 20) mesh points. Since in our case the computations were terminated intentionally whenever the accuracy of the results became uncertain, and this happened in each case before the development of the striae could be completed, nothing can be said about the final values of the effective plasma parameters. For the case illustrated in Fig. 5, the following effective plasma parameter values were found at the end of the 120/~sec time period: aeff/ao = 0"418,

where treff -~- ( J y ) / ( E v ) = O" 353 mho/cm; fleff/~0 : 0"342, where /3elf ---- ( E x ) / ( E y ) = 2"97. Allowance for dissipative processes (thermal conduction, radiation losses, etc.) would, of course, reduce the maxima of the field gradients, thus making the computations

5

and the results more accurate. This, combined with a mesh size finer than the one used in the present computations could provide valuable information on the basic characteristics of ionization instability, viz. growth time, wave amplitudes, effective plasma parameters, and the effect (if any) of the perturbation pattern on the final wave pattern in particular.

References

[1] E. P. Velikhov, Proc. Symp. MPD Electr. Pow. Gen., Newcastle upon Tyne, 135 (1962). [2] J. L. Kerrebrock, AIAA J. 2, 1072 (1964). [3] V. N. Belousov, V. V. Eliseev and I. Ya Shipuk, Proc. Int. Symp. MHD Electr. Pow. Gen., Salzburg, 2, 323 (1966). [4] W. Riedmueller, Proc. lnt. Syrup. MHD Electr. Pow. Gen., Warsaw, 1, 519 (1968). [5] A. V. Nedospasov, Uspekhi Fiz. Nauk 94, 439 (1968). [6] P. Zettwoog, Rapporteur's Report, Sect. 2c, Proc. lnt. Symp. MHD Electr. Pow. Gen., Salzburg, 2, 303 (1966). [7] G. A. Kasabov, Rapporteur's Report Sect. lf, Proc. Int. Symp. MHD Electr. Pow. Gen., Warsaw, 6, 3555 (1969). [8] L. L. Lengyel, Phys. Lett., 29A, (2), (60) (1969); see also the Late Paper presented at the 10th Symp. Engng Aspects of MHD, M.I.T. (1969). [9] E. P. Velikhov, L. M. Degterev, A. A. Samarskii and A. P. Favorskii, Proc. lOth Symp. Engng Aspects of MHD, M.I.T., 1 (1969). [10] D. A. Oliver, Proc. 9th Syrnp. Engng Aspects of MHD, U.T.S.I., 126 (1968). [11] L. S. Frost and A. V. Phelps, Phys. Rev. 136, 1538 (1964). [12] T. A. Cool and E. E. Zukoski, Phys. Fluids 9, 780 (1966). [13] L. L. Lengyel, Proc. 8th Symp. Engng Aspects of MHD, Stanford University, 105 (1967); see also Energy Conversion 9, 13 (1969).