Two-dimensional MHD boundary layers in argon-cesium plasmas

(1)

Two-dimensional MHD boundary layers in argon-cesium

plasmas

Citation for published version (APA):

Arts, J. G. A., & Merck, W. F. H. (1983). Two-dimensional MHD boundary layers in argon-cesium plasmas. (EUT

report. E, Fac. of Electrical Engineering; Vol. 83-E-139). Technische Hogeschool Eindhoven.

Document status and date:

Published: 01/01/1983

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be

important differences between the submitted version and the official published version of record. People

interested in the research are advised to contact the author for the final version of the publication, or visit the

DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page

numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

(2)

Department of

Electrical Engineering

Two-dimensional MHD Boundary Layers

in Argon-Cesium Plasmas

.

By

J.GA Arts and

W.FH

Merck

EUT Report 83-E-139

ISBN 90-6144-139-0

ISSN 0167-9708

July 1983

(3)

EINDHOVEN UNIVERSITY OF TECHNOLOGY

Department of Electrical Engineering

Eindhoven The Netherlands

TWO-DIMENSIONAL MHD BOUNDARY LAYERS

IN ARGON-CESIUM

PLAS~~S

By

J.G.A. Arts

and

W.F.H. Merck

EUT Report 83-E-139

ISBN

90-6144-139~0

ISSN 0167-9708

Eindhoven

July

1983

(4)

CIP-gegevens

Arts, J.G.A.

Two-dimensional MHD boundary layers in argon-cesium plasmas /

by

J. G.

A. Arts and

~I.

F • H.

Merck. - Eindhoven: Uni versi ty of

Technology. - Fig. - (Eindhoven University of Technology

research report / Department of Electrical Engineering,

ISSN

0167-9708, 83-E-139)

Met lit. opg., reg.

ISBN

90-6144-139-0

SISO

661 UDC

621.313.52 UGI

650

(5)

This work was performed as a part of the research program

of the division Direct Energy Conversion of the Eindhoven

University of Technology, the Netherlands.

The research project is supported by the Ministry of Economic

Affairs of the Netherlands, The Hague.

Further informations can be obtained from:

Dr.ir. W.F.H. Merck,

Division Direct Energy Conversion,

Department of Electrical Engineering,

Eindhoven UniVersity of Technology,

P.O. Box 513,

5600 MB

EINDHOVEN,

The Netherlands

(6)

-2-1. SUMMARY

This report deals with the description of a calculation method for

tWQ-dimensional MHD electrode boundary layers in an Argon-Cesium plasma, based

upon the Patankar scheme. The original Patankar scheme for gasdynamic

cal-culations is extended with appropriate electro-magnetic terms in the

gas-dynamic equations. Special attention is given to an improvement of Patankar's

Couette layer approach and to the calculation of the pressure gradient for

duct flows.

The streamer like discharge structure in the Argon-Cesium plasma is not

accounted for. Two approaches have been chosen. A rather ·coarse model

intro-duced by Sherman and Reshotko, imposing a chosen

j

and

j

distribution

x

y

within the boundary layer, and a two dimensional current potential

distri-bution as a stationary solution of Maxwell's equations for each MHD

gene-rator segment, giving a neat current distribution within the boundary layer.

Currents, measured in the experiments of the Eindhoven Blow Down facility,

are used in the calculations in order to check for the gasdynamic generator

performance in both aforementioned approaches.

At high electro-magnetic interaction the coarse model indicated boundary

layer separation in agreement with the sharp pressure rises measured in the

downstream half of the generator duct. The calculations with neat

current-potential distributions did not predict boundary layer separation.

This discrepancy between the two current distribution models proves that an

accurate and detailed description of the current distribution is

indispen-sable to predict the gasdynamic behaviour accurately.

Programs that give a sufficient approximation of the streamer like discharge

structures in Argon-Cesium plasmas are not yet available nowadays.

Arts, J.G.A. and W.F.H. Merck

TWO-DIMENSIONAL MHD BOUNDARY LAYERS IN ARGON-CESIUM PLASMAS.

Department of Electrical Engineering, Eindhoven University of

Technology, The Netherlands, 1983.

(7)

2. CONTENTS

1. SUMMARY

2. CONTENTS

3. NOMENCLATURE

4. LIST OF FIGURES

5. INTRODUCTION

6. GASDYNAMIC EQUATIONS

6.1. Basic equations

6.2. Patankar's solution method

6.3. Couette-flow analysis of turbulent MHD flows

6.4. Pressure gradient

7. CURRENT DISTRIBUTION

7.1. Sherman-Reshotko-model

7.2. Houben-model

8. COMPUTING PROGRAMS

8.1. Patankar program

8.2. Houben-program

9. RESULTS AND EXPERIMENTS

9.1. Introduction to the Blow Down Experiments

9.2. Confrontation of calculations and experiments.

10. CONCLUSIONS AND RECOMMENDATIONS

11. REFERENCES

12. APPENDICES

Al. Couette flow analysis for a turbulent layer

A2. The pressure gradient.

A3.

r'~inite

difference equations of the core flow.

A4. Parameters of the computer programs

A5. Use of the computer programs.

A6. Listings of the computer programs.

Pag.

2

3

4

5

10

14

15

17

20

21

23

26

28

31

35

61

63

66

70

73

75

79

81

(8)

-4-3. NOMENCLATURE

A

0

geometrical cross-section of channel inlet

A

effective cross-section of channel

e

B

z

magnetic induction

H

total enthalpy

J

heatflux

J

h, s

heatflux through the wall

m

E

mass flux through external surface of B.L.

m

0

r

mass flux through wall surface

Pr

molecular Prandtl-number

Pr

eff

effective Prandtl-number

pr

t

turbulent Prandtl-number

R

universal gasconstant

T

gas temperature

u

X-component of gas velocity

u

0

velocity of core-flow at the channel inlet

U

00

velocity of core flow

v

y-component of gas velocity

W

molecular weight

'I'

streamfunction

'I'

E

streamfunction at external surface B.L.

'I'r

streamfunction at the wall

(= 0)

\l

molecular viscosity

\It

turbulent viscosity

\leff

effective viscosity

P

density

Po

density of core-flow at channel inlet

P

_oo

density of core flow

T

shear stress

T

shear stress at the wall

S

(9)

4. LIST OF FIGURES

Fig. 6.1-1 Coordinate system for MHD boundary layer problem.

Fig. 6.2-1

Stream function coordinates for boundary layers:

'¥ - '¥

I

Fig. 6.2-2

Finite difference grid system.

Fig. 6.2-3

Control volume near the wall; Couette layer

approxi-Fig. 6.4-1

Fig. 7.1-1

mation.

Channel geometry and flow area.

Current density distribution

(j )

and electric field

y

distribution (E ) within the boundary layer.

x

Fig. 7.1-2

Current density distribution in a generator segment.

Fig. 8.1-1

Flow scheme of the computer program MHDMAI.

Fig. 8.1-2

Flow scheme of the computer program GASELMACOM.

Fig. 8.2-1

Lines of integration in Houben-program.

Fig. 8.2-2

Flow scheme of the Houben-program.

Fig. 9.2-1

Velocity-profiles along electrode wall at different

x-locations, with Sherman-Reshotko model and

ofh

jxdy

magnetic induction B

=

5.3 T (Group I, table 9.2-1).

Fig. 9.2-2 Temperature-profiles along electrode wall at different

x-locations. Conditions see fig.· 9.2-1.

0,

page.

15.

17.

18.

19.

22.

24.

25.

29.

32.

33.

34.

40.

(10)

-6-Fig. 9.2-3

Stagnation-enthalpy profiles along electrode wall, at

different x-locations. Conditions see fig. 9.2-1.

Fig. 9.2-4

Free stream values of velocity u, pressure P, stagnation

pressure p , temperature

T

and stagnation temperature T .

s

Conditions see fig. 9.2-1.

Fig. 9.2-5

Boundary-layer thickness do' displacement-thickness d

1 ,

momentum-thickness d

2 , energy-thickness d

3 and

enthalpy-thickness d

4 as a function of x. Conditions see fig. 9.2-1.

Fig. 9.2-6

Friction coefficient C

f

' wall heatflux qw' compressible

shape factor d

l

/d

2 and y-component of the current density

in the BL

j

as a function of x. Conditions see fig. 9.2-1.

y

Fig. 9.2-7

Free stream Mach number M and stagnation enthalpy H.

Conditions see fig. 9.2-1.

page.

41.

42.

43. Fig. 9.2-8

Velocity-profiles along electrode wall at different x-locations

43. with Sherman-Reshotko model and

j

=

- j

over the insulator,

x

oo

magnetic induction B

=

5.3 T (Group 2, table 9.2-1) .

Fig. 9.2-9 Temperature-profiles along electrode wall at different

x-lo-cations. Conditions see fig. 9.2-8.

44. Fig. 9.2-10 Stagnation-enthalpy-profiles along electrode wall at different

44. x-locations. Conditions see fig. 9.2-8.

Fig. 9.2-11

Free stream values of velocity u, pressure p, stagnation

pressure

function

p , temperature

T

and stagnation temperature T

s

of x. Conditions see fig. 9.2-8.

as a

Fig. 9.2-12

Boundary-layer thickness do' displacement-thickness d

1 ,

momentum-thickness d

2 , energy-thickness d

3 and

enthalpy-thickness d

4 as a function of x. Conditions see fig. 9.2-8.

45.

(11)

Fig. 9.2-13

Friction coefficient C

f

' wall heatflux

~,

compressible

shape factor d

1 /d

2 and y-component of the current density

in the BL

j

as a function of x. Conditions see fig. 9.2-8.

y

Fig. 9.2-14

Free stream Mach number M and stagnation enthalpy H as a

function of x. Conditions see fig. 9.2-8.

Fig. 9.2-15 Velocity-profiles along electrode wall at different

x-lo-cations with Sherman-Reshotko model and

j

_x

= -j

_xm

over the

insulator, Run 205 at t

=

31 s, magnetic induction

=

4.46 T

(Group 3, table 9.2-1).

Fig. 9.2-16 Temperature-profiles along electrode wall at different

x-locations. Conditions see fig. 9.2-15.

Fig. 9.2-17

Stagnation-enthalpy-profiles along electrode wall at

different x-locations. Conditions see fig. 9.2-15.

Fig. 9.2-18 Free stream values of velocity u, pressure p, stagnation

pressure p , temperature

T

and stagnation temperature T

s

as a function of x. Conditions see fig. 9.2-15.

Fig. 9.2-19 Boundary-layer thickness do, displacement-thickness d

1 ,

momentum-thickness d

2 , energy-thickness d

3 and

enthalpy-thickness d

4 as a function of x. Conditions see fig. 9.2-15.

page.

46.

47.

48.

49. Fig. 9.2-20

Friction coefficient C

f ,

wall heatflux qw' compressible shape

49. factor d

1 /d

2 and y-component of the current density in the

BL

j

as a function of x. Conditions see fig. 9.2-15.

y

Fig. 9.2-21

Free stream Machnumber M and stagnation enthalpy

H. Conditions

50. see fig. 9.2-15.

Fig. 9.2-22

Velocity-profiles along electrode wall at different

x-locations with Sherman-Reshotko model and

j

=

-j

over

x

xm

the insulator, Run 205 at t

=

36's, magnetic induction B

=

5.3 T (Group 4, table 9.2-1).

(12)

-8-Fig. 9.2-23 Temperature-profiles along electrode wall at different

x-locations. Conditions see fig. 9.2-22.

Fig. 9.2-24 Stagnation-enthalpY-profiles along electrode wall at

different x-locations. Conditions see fig. 9.2-22.

Fig. 9.2-25

Free stream values of velocity u, pressure p, stagnation

pressure p , temperature

T

and stagnation temperature T

s

as a function of x. Conditions see

fig.

9.2-22.

page.

51.

52. Fig. 9.2-26

Boundary-layer thickness do' displacement d

1 , momentum-thick-

52. ness d

2 , energy-thickness d

3 and enthalpy-thickness d

4 as a

function of x. Conditions see fig. 9.2-22.

Fig. 9.2-27 Friction coefficient C

f

' wall heatflux

~,

compressible

shape factor d

1 /d

2 and y-component of the current density

in the BL

j

as a function of x. Conditions see fig. 9.2-22.

Y

Fig. 9.2-28 Free stream Mach number M and stagnation enthalpy H as a

function of x. Conditions see fig. 9.2-22.

53.

53. Fig. 9.2-29 Velocity-profiles along electrode wall at different x-locations

54. with Sherman-Reshotko model and

j = -j

Over the insulator

x

oo

Run 205 on t

=

41 s, magnetic induction B

=

4.98 T (Group 5,

table 9.2-1).

Fig. 9.2-30 Temperature-profiles along electrode wall at different

x-locations. Conditions see fig. 9.2-29.

54. Fig. 9.2-31

Stagnation-enthalpy-profiles along electrode wall at different

55. x-locations. Conditions see

fig.

9.2-29.

Fig. 9.2-32

Free stream values of velocity u, pressure p, stagnation

pressure p , temperature

T .

Conditions see fig. 9.2-29.

s · s

Fig. 9.2-33 Boundary-layer thickness d , displacement-thickness d

1 ,

o

.

momentum-thickness d

2 , energy-thickness d

,

3 and

enthalpy-thickness d

4 as a function of x. Conditions see fig. 9.2-29.

55.

(13)

Fig. 9.2-34 Friction coefficient C

f

' wall heatflux

~,

compressible shape

factor d

1 /d

2 and y-component of the current density in the

BL jy as a function of x. Conditions see fig. 9.2-29.

Fig. 9.2-35

Free stream Mach number M and stagnation enthalpy H as a

function of x. r.onditions see fig. 9.2-29.

page.

56.

57. Fig. 9.2-36

Current distribution in a generator segment with dimensions

57.

2 0.150 x 0.025 m , current

I

=

70 A/m, induction

B

=

5.3 T,

gridsize n(x)

=

50 and m(y)

=

60 steps.

Fig. 9.2-37

Velocity-profiles along electrode wall at different

x-locations with Houben model, Run 205 at t

=

36 s,

magnetic induction

B

=

5.3 T (Group 6, table 9.2-1).

Fig. 9.2-38 Temperature-profiles along electrode wall at different

x-locations. Conditions see fig. 9.2-37.

Fig. 9.2-39

Stagnation-enthalpy-profiles along electrode wall at

different x-locations. Conditions see fig. 9.2-37.

Fig. 9.2-40

Free stream values of velocity

ll,

pressure p, stagnation

pressure p , temperature

T

and stagnation temperature T

s

as a function of

x.

Conditions see fig. 9.2-37.

Fig. 9.2-41

Friction coefficient C

f

' wall heatflux

~,

compressible

shape factor d

1 /d

2 and y-component of the current density

in the BL jy as a function of x. Conditions see fig. 9.2-37.

Fig. 9.2-42

Free stream Mach number M and stagnation enthalpy H as a

function of x. Conditions see

fig.

9.2-37.

Fig. A1-1

Method of integration over the Couette layer.

Fig. A3-1

Schematic diagram of core flow region (shaded area) •

58.

59.

60.

68.

73.

(14)

-10-5. INTRODUcrION

During the last 15 years many papers and reports On the subject of MHD

boundary layers have been published, both upon theoretical and

experimen-tal aspects.

Many papers in this field, especially concerning computational methods,

have been published by the STD-Company (refs. 1 through 4).

The experimental approach has been emphasized by several papers of

Stan-ford University (refs. 5 through 11). Both these groups deal with

com-bustion MHD generators.

The STD-Company has proven its capability to perform 3-dimensional

calcu-lations and obtains good results in explaining experimental results on

very specific problems (ref.

1,

3, 4). They claim the capability to

predict the performance of commercial scale MHD generators (ref. 2) and

show that for large scale generators 2-dimensional calculations cannot

give account for very important effects that originate from 3-dimensional

MHD inhomogeneities (ref. 1).

The Stanford group has performed a lot of experiments on insulator and

electrode wall boundary layers. Velocity, gas temperature, electron density

and electron temperature profiles have been measured by different kinds

of sophysticated measuring techniques (refs. 5 through 11). They also obtain

goOd agreement between experimental results and theoretical predictions.

A specific field of interest in electrode boundary layers is the occurence

of inter-electrode breakdown, or otherwise mentioned, Hallfield-shorting.

This phenomena deteriorates the performance of any kind

of MHD generator.

Oliver (refs. 12, 13) and Russo (ref. 14) have designed time dependent

2-dimensional calculation methods for this specific problem, enabling us to

study the development of an inter electrode breakdown. The Stanford group

has performed some experimental studies in this field, distinguishing

between plasma breakdown and insulator breakdown (ref. 6).

The above mentioned short and incomplete list of publications shows anyhow

that much effort has been put into the study of combustion gas MHD boundary

layers. Much less attention has been given to the field of noble gas MHD

boundary layers. Merck (ref. 15) and High (ref. 16) dealt with the insulator

wall boundary layers, whereas Doss (ref. 17), Koester (ref. 18), Pian

(ref. 19) and Lindhout (ref. 20) dealt with the electrode wall boundary

layer. Hall-field shortings in electrode wall boundary layers were observed

(15)

in the experiments performed by Kerrebrock at M.I.T. (ref. 21). This effects

hampered the production of sufficient electric power and the experiments were abandoned.

The smaller progress in the field of noble gas MHD generators is mainly due

to the following aspects:

- The USA and USSR governements have laid emphasis upon the development of

coal fired combustion gas MHD

- The behaviour of noble gas-cesium plasmas is very complex due to ionization

instabilities, the phenomenon of streamers creation and complicated

plasma-electrode interaction (Koester effect).

Lengyel (ref. 22) was the first to predict the creation of streamer-like

structures by ionization instabilities. Hara (ref. 23) presents a more

sophysticated model for a flowing plasma, showing the time dependent growth

of streamers in an A-Cs discharge. These facts are confirmed by early

expe-riments of Brederlow (ref. 24) and later by Sens (ref. 25) and Hellebrekers

(refs. 26, 27) in shock tube experiments performed at the Eindhoven

Universi-ty of Technology (E.U.T.). With a high speed came"a complicated streamer

structures were photographed, connecting several electrode pairs for short

time periods, moving with plasma velocity and jumping to the downstream

electrode pairs. The time intervals between passing streamers show a

some-how stochastic character.

It is evident that an exact analysis of MHD boundary layers (BL) with

stochastic discharges asks for a high investment of manhours and computing

time. Some time-averaged model of the streamer behaviour should be used as

an approach. On this very moment no sufficient amount of data upon

time-dependent streamer behaviour is available to achieve a realistic mean time picture. So the authors have chosen a steady state solution of the Maxwell equations, delivering a current density distribution according to Houben

(ref. 28). Like most authors, Doss (ref. 17), Blom (ref. 29), he has used

3

small current densities (j < 0.5 A/cm ) to prevent numerical instabilities

due to the exponential character of the Saha coefficients.

It is supposed that in this way a more realistic estimate of the current

and potential distribution is obtained than the ones used by Doss (ref. 17)

and Pian (ref. 19). The calculations show that the BL thickness cannot be

regarded small referred to the axial electrode length, so a neat current

(16)

-12-~ith Houben's program. For higher current densities a linear extrapolation

is used.

Further the calculation of the plasma-electrode interaction causes much

trouble. Koester (ref. 18) showed that a specific interval of electron

temperature and Cs-density exists where high diffuse current densities

can

be

realized. This was confirmed by measurements of Blom (ref. 30) at

E.U.T. To take these effects into account, complete solving of the electron gas equations within the BL is necessary and very exact local electrode

temperatures should be defined. Even at constant electrode temperatures

and slightly simplified boundary conditions, as defined by Sanders (ref.

31), this would lead to extreme computing times (whenever possible) ,

without having the assurance of obtaining better results.

So within the BL the electron energy and Saha equations are solved to find

the local conductivity in order to calculate the local ohmic dissipation.

In Chapter 6 the gasdynamic equations are treated and extended with electro

magnetic terms. They are solved with the Patankar method (ref. 32) analogous

to the STAN 5 program (ref. 33). The pressure gradient term and the Couette

flow region get special attention.

In Chapter 7 the Maxwell equations are solved with the Houben program.

Further the simple Sherman-Reshotko current density model is mentioned and

used as a reference.

The gasdynamic Patankar program and the electro magnetic Houben program

are treated in Chapter 8. To combine these programs special problems had tobe

overcome because they were written in different languages (Algol and Fortran

5). An interpolation method to change the coarse current-potential grid to

the much finer BL grid is treated.

Computation results were confronted with experiments of the Eindhoven Blow

Down Facility (EBDF) in Chapter 9. As explicit calculation of streamers is

not possible nowadays, no complete theoretical prediction of the generator

performance can be made. So the program is used to check the gasdynamical

behaviour of the EBD?-experiments by introducing the measured electrode

currents into the program and check for the pressures and eventual BL

(17)

In Chapter 10 some conclusions and recommendations are given for future

use of the presented two-dimensional calculation method.

(18)

-14-6. GASDYNAMICAL EQUATIONS

In general, the boundary layer problem at the electrode wall of an MHD

generator can only be solved taking into account three sets of equations.

- Maxwell's equations

-

Gasdynamical equations

- Electrongas equations.

Sanders (ref. 31) has derived a set of gasdynamical equations and electron

gas equations taking turbulence, ambipolar diffusion and turbulent

cross-correlation terms into account. Boundary conditions for the electron gas

equations, related to the method given by Koester (ref. 18) were considered.

It

was

realized that the electrongas

BL - equations do not have a

parabo-lic character, due to high axial derivatives at the edges of electrodes.

This fact, mixed up with the unknown cross-correlation terms, the complicated

solution method necessary to find just the boundary conditions at the

elec-trode surface and the fact that the streamer structure introduced still other

complications made the authors decide not to tempt to solve the electrongas

BL-equations completely but just take the electron thermal-non-equilibrium

into account within the gasdynamic BL

In the coming sections the gasdynamical equations and the method of solution

will be considered. Special attention will be given to the region near the

wall, with the Couette flow approximation, and to the pressure gradient. To

calculate the flow field the BL-equations are solved with the free-stream

values of velocity and enthalpy as part of the boundary conditions.

The BL-equations include the continuity, momentum and stagnation enthalpy

equations for the heavy particle gas (ref. 32). They describe a turbulent

MBD flow with the following assumptions:

- magnetic Reynolds number is small:

a

- steady, turbulent flow; at - 0

- two-dimensional flow;

a!

=

0 - radiation losses neglected.

B

=

B

=

constant

z

-

frozen electron temperature and non-equilibrium.

The coordinate system is shown in figure 6.1.

(see next page)

(19)

V

---

'\---~

.BL

x

~

el elrod ewall

Fig. 6.1-1 Coordinate system for MHO boundary layer problem.

The averaging method for turbulent flows as described by Pian (ref. 19)

is used, yielding the following set of equations:

-

continuity:

a

[pu]

+..l.-

[pu]

=

0 ax

ay

- momentum equation:

pu

~

+

pv au

ax

ay

=

+

j Y B

z

(6.1-1)

(6.1-2)

where

~ =

molecular viscosity and

~t =

turbulent part of effective

viscosity

~eff

=

~

+

~t

- stagnation enthalpy equation:

-a

+ -

_ay

with stagnation enthalpy H defined by (Pian, ref. 19)

-H

(6.1-3)

(20)

-16-The following notation has been used:

\leff

-1!..+

\lt

Pr

eff

Pr

t

(6.1-5)

( 1-

1

\l (1-

1

\l ; _{\l(1- - )}

+

e ff

Pr

_eff

Pr

t

Pr

_t

(6.1-6)

with Pr

;

molecular

Prandtl-number

("

0.7)

Pr

eff

=

effective

Prandtl-number

("

0.9)

Pr

t

turbulent Prandtl-number ("

0.9)

The effective viscosity \leff is evaluated by the Prandtl mixing-length

hypothesis.

).Ieff

(6.1-7)

where 1 is the mixing length expressed by: {ref.

32]

1 ; Ky

1 ;

_(6.1-8)

K

and

A

are constants

(K; 0.435, A ; 0.09),

y is the distance from the

wall and Yl is a characteristic thickness of the layer (ref. 32,

p.

20).

Near the wall special expressions for \leff and Pr

t

are used which are

presented in chapter

6.2. The boundary conditions of the boundary-layer equations are:

- non-slip condition u(x,o) ; v(x,o) ;

0 - fixed temperature at the wall T(x,o) ; T (x)

w

(6.1-9)

(6.1-10)

- the free-stream conditions given by the core-flow equations:

du

~

dp

mom. eq.:

p

u - - ; - -

+

j •

B

(6.1-11)

00 00

dx

yen

dH

enth. eg.: Poouoo dX

oo

=

j

-+ E

cont. eq.:

p

u A ; constant

~ ~

e

pRT

perfect gas eq.: p

=-W-where R

=

universal gas constant

(6.1-12)

(6.1-13)

(21)

The pressure-gradient is assumed to be constant over the cross-section

of the channel and is evaluated by expression 6.4.1. The mass-flow is

constant and depends on the inlet conditions.

effective cross-section in which displacement

(App.3 ).

The surface A is the

e

values are included

The initial profiles of velocity and enthalpy are given by the 1/7 power

law and the Croco-Busemann relation expressed in

streamfunction-coordinates

(ref. 19) . We note that the downstream field of a turbulent flow is not

sensitive to details of the initial-profiles.

6.2. Patankar's solution method

To compute the flow-field the method of Patankar-Spalding is applied

(ref. 32). This method is a marching-integration-procedure for which the

parabolic boundary-layer equations are transformed to the

streamfunction-coordinate-system4

The streamfunction

¥

is defined by:

a¥

₌

_m

pu =

ay

a¥

₌

_-

_pv

ax

and the non-dimensional streamfunction w:

w

=

¥ -¥

E

r

(0 ~

w

:<

1)

with ¥r

streamfunction at the wall

(6.2-1)

(6.2-2)

¥E

=

streamfunction at the outer part of the boundary-layer.

free. stream

boundary layer

W·O~---L---~I

x

Fig. 6.2-1

Stream function coordinates for boundary layers:

¥ -

¥

1

w -

(22)

-18-By use of

these expressions

the conservation-equations of momentum and

enthalpy are transformed from the (x,y) to the (x,w)-coordinate system

(Von Mises transformation ref. 19) and written in the general form:

where:

a

=

b =

c

o~

0

~

(a+bw) -

-

(c - ) + d

ow - ow

ow

'I' -'I'

E

I

m -m

E

I

'I' -'I'

E

I

]Jeff Pu

'(6.2-4)

(6.2-5)

(6.2-6)

here 0eff

=

pr

_eff

in the enthalpy equation and 0eff

=

1 in the

momentum equation;

1

d =

pu

-(_ dp

+ .

B )

dx

J y

2 o

1 d

=

ow []Jeff (1-

Pr )

eff

enthalpy equation

in the

~omentum

equation

in the

(6.2-7)

This equation is used to develop a finite difference scheme. Values of

~

at a downstream station are obtained from values of

~

at an upstream

station by use of a finite difference formulation. This procedure can be

repeated until the whole field of interest is covered. To obtain a finite

difference formulation the basic equation of

~

is considered over a

con-trol volume defined by a chosen grid (fig. 6.2.2.).

E-surface

1 i+l,.---t-i.1

-

f-.

L~

l

. -1 \

...

_W&.6

~vo

ume

- ;

o

I - surface

;_I..L.. _ _ _ _

+_

x

(23)

The grid is equidistant in the w-direction. At the I and E boundaries

a different control volume is considered (fig.

6.2.3.).

At the wall a

special approximation is used to formulate the boundary-conditions in

agreement with wall functions calculated from a Couette-flow

approxi-mation. The relation between the different values of

$i

can be written

in tridiagonal-matrix-form:

The values of coefficients Ai' B

i ,

station. Using the

~-values

at the

be solved successively.

w =1

w3

w,o

r- - - . ,

,

'-_ 1 r -

_,

-

-~

w)

w

_Z

.5

~

w,=O

(i

=

1, N)

(6.2-8)

C. are evaluated at the upstream

~

boundaries this set of equations can

~3

$2.5

}

approximation

couette

-P,

$2

Fig. 6.2-3 Control volume near the wall: Couette layer

approxi-mation.

(24)

-20-6.3. Couette-flow analysis of turbulent MHD flows

---The Couette flow analysis is based upon the assumption that the flow can

be defined as a function of the y-coordinate only. This reduces the

governing equations to the form:

momentum -

pv

o~

=

- + -

op

0

_[Ileff

o~]

+

_J. B

(6.3-1)

y z

oy

ox

oy

Il

_oR

-2

enthalpy -

pv -

oH

= 0

[ eff

_{oy + Ileff}

(1

-

1

)

1 ~]

+

oy

pr

eff

pr

eff

2 oy

j E

+

. E

(6.3-2)

x x

J

_y

Y

flow -

0

(pv) 0

m

(6.3-3)

mass

=

pv

₌

oy

s

These equations can be integrated with respect to y and

non-dimensiona-lized, according to ref.

19 and demonstrated in App.

1. In the case of

MHO

duct flows the massflow across the Couette layer is

zero. Hence we obtain the following expression for shearstress

T

and

heatflux

J:

(6.3-4)

= l + m $

-Wu+-oy

+ +

+

(6. 3-5)

The index

"+" refers to the nondimensional quantities. These expressions

can be equated with the well known transport laws for momentum and heat

flux:

=

11+

(6.3-6)

and

Il+

₁

"

_"+

du+

2

"2

W

=p-r-e-f -f

dy

+ '

(6.3-7)

Pr

eff

Yielding expressions for the velocity and enthalpy derivatives within

the Couette layer:

1 [1 + p + y

+ m u ]

Il+

+

+ +

(6.3-8)

(25)

and

pr

eff

= ---''--

[1

+

m <!> \l

+

+ +

- 0

+ +

y

1 +

1-Pr ff

( e )

2 (6.3-9)

The

~rnD

effects are incorporated in the terms p+ containing the Lorentz

force and

0+

containing the ohmic heating. The term

~+

contains the

tur-bulent friction coefficient and Van Driest damping term (ref. 34).

These derivatives can be integrated with aspect to

y across the Couette

layer, provided the constitutive relations are honour (App. 1

). In this

way the u and

~

at the flow side edge of the Couette layer are found

to-gether

with

friction coefficient, shear stress, wall heat flux and

Stan-ton number.

Patankar (ref. 32) has developed an elegant method to find the pressure

gradient in confined flows. The pressure gradient is a part of the

comple-te solution of boundary layers and bulk flow.

At any forward step the pressure gradient is estimated first. Then the

forward step calculations are performed, yielding solutions for the

gas-dynamic quantities. The criterion is the area occupied

by

the flow A

f ,

found through

ff

pu dA

=

m,

and the actual area provided by the duct Ad.

If they are not equal then the pressure gradient has to be corrected in

a sense that the difference Ad-A

f

is counteracted.

Combination of (6.1.-1), (6.1.-2), (6.1.-3), (6.1.-14) yields (see App. 2):

where

dp

=

dx

u W

s

T C

s p

IBl

IB2

<!>

=

u

W

s

T C

s P

I (I -F u )

-2

Bl

w"

(IE+<l»

_{+ IB2}

- F }

W

,

1 _fl

j

Budy

1 fl

=

₁

₀

IE

=

_e

₀ J

E dy

1 _fl

j

B dy

1

<f,

1

0

FW

A

T

s

ds

1

<f, J

ds

A

h,s

(6.4-1)

(6.4-2)

In finite difference approximation the change in cross-section is defined

(26)

-22-(6.4.-3)

The last term a(Ad,u - Af,U) is the correction term for the pressure

gradient expression (6.4.-1).

o

~

a

~

1 is introduced to prevent overshooting.

Figure 6.4.-1 gives an exaggerated impression of the correction effect.

-/ I 1

__

~

_

. ___

~

_~

1id.JL

,<\1 ___ _

I

1

1 I

(27)

7. CURRENT DISTRIBUTION

To solve the boundary-layer equations we need to know the components of

the local electric field and current-density to evaluate the

saurce-terms of the basic equations. Two models are considered in our

calcula-tions. The first model avoids to solve the Maxwell-equacalcula-tions. This

simple model is used by several authors like Sherman-Reshotko, Doss

and others (ref. 19). The second model solves the Maxwell-equations for

one segment (Houben-program) and finds a current distribution which is

used for all segments in the gasdynamical program.

This procedure does not solve the electro-magnetic equations and relies

on the condition that the boundary-layer thickness is small compared to

the length of a segment.

The assumptions are:

Over

conductor:

_jy

=

_{jy (x)}

E

0 x

Over

_{insulator: jy}

=

0 E

E (x)

x

(7.1-1)

With Ohm's law: jx

a

[E -S (E -uB)]

1+S2

x

y

\

=

a

[E -uB + SEx]

2

Y

l+S

(7.1-2)

...

The term j.E can be approximated by the following expressions:

j2 (x)

.

(1+S2)

Over

conductor:

_{E xjx + E yjy}

= Y

_{+ uB jy (x)}

(7.h3)

a

Over

_{insulator: E xjx + E yjy}

_a

E (x)

2

(7.1-4)

x

The term j B is a function of x. OVer the conductor the function

y

j (x) is expressed by:

y

j (x)

(28)

-24-.

- Ex max

-J

_ymax

I~

I

a

f4-

b

I I

Xins_o

Fig. 7.1-1

current density distribution (j ) and electric field

y

distribution (E ) within the boundary layer.

x

and E (x) is lineair in x over the insulator. To avoid discontinuities

x

the insulator is divided in two areas, (see fig.

7.1-1):

I: E

(x)

=

E

max.

x

for 0

~

x.

~ns ~ cS II: E (x)

x

b-x

E

_x

_max.

(

_b-6

ins)

o

~ x. { b ~ns

where 6 is taken 0.05 b.

(7.1-6)

The value of E

max depends on the E

in the bulk which is evaluated by

x

the assumed current distriubtion, the effective conductivity and

Ohm's-law.

For the bulk the current distriubtion is assumed to be super-imposed by

the segmentation ratio. The effective and apparent conductivities are

calculated fram a simple network where external load, voltage drop and

(29)

jMWA_

_FaA

1 ___ I

r----:---r

I JxeUlK I

:

j

[S]

:

I

you",

1 BULK

1 1

I

1

1 1

do

- -=-=.-

=-=-- -

~

_

r

I boundary - layer

Fig. 7.1-2

Current density distribution in a generator segment.

->->

To calculate the j.E-term in the boundary layer over an insulator several

1

current distributions are tried. A first one assumes

f

j dy

=

O.

o

x

So

where d

o

(h-2d )

Q

2d

o

boundary-layer thickness (see fig. 7.1-2)

(7.1-7)

Other calculations are performed with smaller values of jx in the boundary

layer to lower the ohmic hea,ting over the insulator:

J'

=

j

~

x

x bulk

-> -> J

To evaluate the term j.E

= -

+

j uB over an electrode the electrical

con-a

ductivity is calculated from 5aha-equilibrium and a simplified form of the

electron-energy-equation (ref. 35).

It is assumed that the ohmic heating is in equilibrium with the elastic

losses and radiation:

,2

l k

2

v

2

v

_J_

_{m n}

_{(T -T) (~+}

ec)

₊

_Rad,

a

2

e e

e

m

(7.1-8)

a

c

(30)

-26-if

S

<

S

_crl.

' t '

a

_e

ff

a

(7.1-9)

The electron-energy equation is solved itteratively. For a more detailed analysis of the influence of the electron-gas see section 7.2.

7.2. Houben-model

7.2.1. Basic equations

Assumptions

a

- Two-dimensional model,

az

=

0 (ref. 28).

Distribution function of each species is Maxwellian

- Magnetic Reynolds number is small

- Plasma is electrically neutral

-

Stationary, diffuse discharge

Contribution of inelastic collisions is neglected in comparison with

contribution of elastic collisions due to momentum and heat transport.

-

AgT

=

0,

thermal conduction of electron-gas is neglected.

e

Velocity of heavy particles is supposed to be in x-direction.

-+ -+

v

=

(u(y),

0, 0)

-+

g.v

=

0 ..-Heat transport by the heatflux-vector qe

glected

gp

len

is neglected.

e

(~kT

+

E

k

)

j/e

is

ne-2

e

Neglection of the heat-flux-vector

ties

«

5 x 10

3 A/m

2 ,

ref. 28, 29,

and Vb

_{. e}

len

_e

assumes low

current-densi-31 Houben p. 55, Blom p. 57, Sanders

p. 103). This assumption is not applicable to the experimental data from

the Blow-Down experiment, of which the current-densities can be much higher

(: 5

x

10

4

A/m

2 ).

..-Neglection of the electron heat conduction qee

=

AeVT

e is discutable. The influence depends on the boundary-conditions. The Houben program shows

relatively low gradients in T

(q

« q

,«

1 %)). The NLR-program

e

ec

tot

which uses other boundary conditions, shows sharp gradients in T along

e

the electrode-wall. In that case the electron-conductivity is locally

con-sider able (ref.

2~.

With this assumptions the expressions of the basic equations are:

conservation equation electron-gas: v

on

e

(31)

energy-equation electron-gas:

3 kvn

2 e

aT

e

ax

,2

= L -

3nmk

a

e e

(T -T)

e

(El.'

+

l

kT )

2 e

Maxwell and momentum-equation:

a'

_J

a'

x

J=

0 +

ax

ay

a'

_J

a'

x

_J

+

p

(x,y)

ax

:ly

where: P{x,y)

Q{x,y)

=

7.2.2. Method of solution

jx - Q{x,y)jy

a

(~)

ay

a

v,

J J

L

, 1

m,

J=

J

=

0 - Rad -

(I-R).

(7.2-2)

The method to solve the set of equations comprehends two parts:

1)

For given current-distribution a solution of nand T

is obtained

e

from integration of the continuity and energy equations by a

Runga-Kutta method in

x~direction

(direction of heavy particles velocity) .

Boundary conditions of nand T

must be given at x

=

o.

e

2) When n

and T

are known the current distribution is found from the

e

elliptical set (7.2.-3) using a streamfunction

~.

This is done by

a method of successive overrelaxation (ref.

36).

The boundary

condi-tions for this equation are:

- periodicity

- j y

:l'l'

ax

o

on the insulator

- E

x

i !

ay

-s

a'l'

ax

o

on the electrode

The procedure is repeated until a declared accuracy of the

streamfunc-tion is obtained.

(32)

-28-8. COMPUTING PROGRAMS

8.1. ~~~~:2~

In this chapter the gasdynamical programs are described, represented

by the flow schemes in Figs. 8.1.-1 and 8.1.-2. The flow schemes present

the gasdynamical programs in agreement with the two electromagnetic

models presented in chapter 7.2.1 and 7.2.2 respectively. The differences

are concentrated in chapter 9A of the programs where the electromagnetic terms are calculated.

8.1.1. Flow scheme of MHDMAl

The initial part of the main-program contains file-declarations:

file

6: output file (printer or remote)

file 10: contains data of radial profiles after running

file 11: contains data of axial profiles after running

file

9: input file, input of currents and apparent conductivities.

Chapter 0: Declarations of most variables ·of the main-program and

decla-ration of external procedures.

Chapter 1-3: Contains the most important parameters that have to be

adapted

if

a new situation will be calculated, e.g. magnetic field,

ini-tial values of pressure, velocity and temperature, seedfraction, wall-temperature at the inlet, end-value of x, geometrical conditions.

Chapter 4: Specific constants, input of currents and apparent

conductivi-ties.

Chapter 5: Calculation of initial profiles.

Chapter 6: Start of main-loop.

At the x-station x

=

xu

pressure, geometry,

temperature-profiles and BL thickness are calculated. After calculation of averaged flow values used in the expression of the

9ressure

qra9ient the subroutine STRIDE (1) is called. In this first part of stride just some general grid parameters are evaluated.

Chapter 7: Forward step.

In this chapter the length of the forward step is evaluated. The length

depends on the place at the electrode or insulator. At the end of an

(33)

elec-START

J

0

1

2 STRlO[

3 CD

4

5

6 AUX

7 - - - -

WF

FUN

(])

8

9

I

r-~

~

I

-

--I

I

TTER

_Ie-o

J

L-FUNT

1

10 -CD

"l(>

l

STOP

J

(34)

-30-trode or insulator the stepsize is adapted in agreement with the

geometry of the segment.

Chapter 9:

The laminar viscosity is evaluated. Subroutine AUX is called to apply

to the mixing-length hypothesis followed by an evaluation of the

en-trainment which depends on the velocity gradient in the outer part of

the boundary

~ayer.

The stepsize is adapted

if

the entrainment exceeds

a specified value.

In chapter 9A the electromagnetic terms and the pressure gradient are

calculated. In this part of the program we find the main differences

between the flow schemes. In the Sherman-Reshotko model the core flow

values of the electric field components EXMEAN and EYMEAN are expressed

by Ohm's-law. Along the electrode the procedure TTER is called (inside

procedure C05AAF) to find, local values of conductivity by

Saha-equili-brium in agreement with the assumed current -distribution. To save

process-time the procedure TTER is not called in every gridpoint but for most

points a lineair interpolation is applied. Along the insulator a constant

current-density JX is assumed. After calculation of the pressure-gradient

the second part of STRIDE is called. In this part of STRIDE boundary

coefficients for velocity and enthalpy are calculated. Inside this

proce-dure the subroutine WF is called two times to calculate velocity and

ent-halpy wall-functions respectively. The Van Driest relation is included in

the external function FUN. Integration of du/dy in nondimensional form is

done by adding. The wall-function of the enthalpy is produced by

inte-gration of the external function FUNT by use of a standard inteinte-gration

procedure D01ABF or D01BDF.

Chapter 10-10D:

Here the output is produced.

File

6:

printer output

File 10, radial profiles

File 11 : axial profiles

The files 10 and 11 can be read by the plotprogram MHDPLOT in order to

make graphs of several gasdynarnical quantities.

Chapter 11:

(35)

of the coefficients of the tridiagonal matrix (6.2.-8) The core flow

value of the velocity is also calculated in this part of STRIDE. After

STRIDE the program returns to the begin of the loop in chapter 6 until

XU

=

XULAST.

8.1.2. Flow scheme of GASELMACOM

The main-differences with regard to MHDMAl are concentrated in chapter

9A. At every segment the subroutine CURDIS is called to calculate the

current-distribution. This procedure is an ALGOL-program stored in a

library and declared at the begin of GASELMACOM. The procedure

calcu-lates a current-distribution from a streamfunction calculated by Houbens

solution procedure (section 7.2 and 8.2).

In this way the current distribution is the same for every segment, just a specified constant factor adapts the streamfunction to the

experimen-tal value of toexperimen-tal segment current. After calling CURDIS the components

of current densities inside the boundary-layer over one segment are stored in the one-dimensional arrays JXF and JYF. In the x-direction a

lineair interpolation is applied to find current densities from the grid

of the Houben-program to the grid of GASELMACOM. Within CURDIS the

inter-polation in the y-direction from the grid of Houben-program to the grid

6f GASELMACOM in performed by means of polynomials.

To evaluate the JE-term the local conductivity is calculated with the

assumption of Saha-equilibrium (Subroutine TTER).

8.2. ~~:22

The program (filename Houben/Arts) starts with reading the

input-para-meters like geometrical measures of the grid over one segment, accuracy,

magnetic field, ratio of the densities of seed and Argon, pressure, ve-locity, temperature, minimal stepsize, current density, critical Hall-parameter, initial streamfunction-matrix and initial electron-density and electron-temperature at the upstream boundary of the segment.

In the procedure SIGMAS the electron-energy and continuity-equations are integrated simultaneously in x-direction over the grid lines by procedure STIFF resulting in a distribution of electron-density and electron-tempe-rature. It is a well known fact that the current distribution is rather homogeneous in the bulk and only irregular close to the electrode walls

so a fine grid near the walls (MR+1) and a coarse grid in the bulk (DMB)

has been chosen, as shown in Figure 8.2-1, to perform the integration of

(36)

-32-I

START

I

0

1

2 STRIDE

3 CD

4

5

6 AUX

7 - - - -

WF

FUN

8 CD

9

I

r-I

~

- - -

-(

RDI~ I I

J

_L

-'TTl

RI

FUNT

~

10 --CD--"0

I

STOP

I

(37)

u

-

P, T

Te

Ne

X=

~

---0

f0%:I

-J

MR·' LINES

) OMB

LINES

Fig. 8.2-1 Lines of integration in Houben-program.

For the rest of the gridlines, the equidistant grid, a lineair

inter-polation is applied to calculate local values of the quantities electron

density, electron temperature, conductivity, collission frequency and

Hall-parameter. Much processtime is saved in this way. After this

proce-dure conductivity and Hall-parameter are corrected for instabilities.

The procedures MAAKPi and MAAKQi calculate coefficients of the matrices

Pi and Qi needed to solve the elliptical part of the equations. The

elliptical part is solved in procedure FUNCTIEF11 where a new stream

function matrix F1I is calcualted by a method of successive

over-relaxa-tion (ref. 36C). This process is repeated 50 times and will be repeated