The computation of steady state arcs in nozzle flow

Computer Physics Communications 13 (1978) 363—370 © North-Holland Publishing Company

THE COMPUTATION OF STEADY STATE ARCS IN NOZZLE FLOW M.T.C. FANG, S.K.

and R.D. WRIGHT

Department of Electrical Engineering and Electronics, University of Liverpool, UK Received 8 July 1977

PROGRAM SUMMARY

Title of program: DCANF

Nature of physical problem DCANF finds the supercritical.solution of a high pressure arc burning in a converging—diverging nozzle where the primary standard shape factors [1,2] are assumed constant and the external flow one-dimensional.

Catalogue number: ABUS Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Computer: Installation: CDC7600/ICLL9O6A

Method of solution The three first-order ordinary differential equations, [Q]dx/d~ b [1], describing the d.c. arcs in nozzle flow are solved for a given solution parameter, K assuming constant standard primary shape factors [2]. To obtain the supercritical solution, it is necessary to ensure that Q I = 0 and the compatibility conditions are satisfied simultaneously at the critical point, ~. This is ensured by guessing the arcing conditions at the critical point (i.e. the velocity magnitude factor gc and p~,the ratio of the arc displacement area to the nozzle area). The governing equations are integrated forward from the upstream electrode and backward from the critical point to the socalled patching point. The sum of the squares of the differences of the dependent variables at the patching point is then minimised with respect g~and P~using Powell’s method [3].

University of Manchester, Regional Computer Centre

Operating system: Scope 2.1. Programming language used: FORTRAN High speed storage used: 12000 words Number of bits in a word: 60 bits Overlay structure: none Numberof magnetic tapes required: none

Typical running time With reasonable guesses for g~and ~ a solution can usually be found within 1.5 s on CDC7600.

Other peripherals used: card reader; line printer. No. of cards in combined program test deck: 931

References [1] S.K. Chan, M.T.C. Fang and M.D. Cowley, University of Liverpool, Dept. of Elec. Eng., Arc Research Project Report No. ULAP-T.44. [2] S.K. Chan, M.D. Cowley and M.T.C. Fang, J. Phys. D. 9 (1976) 1085. [3] M.J.D. Powell, Comp. J. 7 (1964) 155.

Gird punchingcode: ICL1900 Keywords: arc, plasma physics, electric, gas flow, thermal plasma, nozzle

*

Engineering Department, Cambridge University. Now at: Mechanical Engineering Department, University of Hong Kong, Hong Kong. 363

M. T.C. Fang et al. I Computation of steady state arcs in nozzle flow

364

LONG WRITE-UP 1. Introduction

area is equal to displacement area. The other quantities in the equations are given by:

A mathematical model for d.c. arcs in nozzle flow (fig. 1) has recently been established using the boun dary layer integral method for axisymmetric arc analysis of Cowley [1]. Three general equations describing the axial variations of the external flow Mach number, M~,the velocity magnitude factor, g, and the arc displacement area, °d’have been derived in terms of the nozzle area variation and power input due to Ohmic heating [2]. With the assumption of constant standard primary shape factors [3], these equations can then be reduced to: F 2+(’y—l)MI

dMI_ MI

da —(C a

—

—

—

j3

C=(5Ah +AkeMI)/(5 +MI), D = M~,1(5 + MI),

F(Ai1~~, g,

~l)

213(Ah

(1

—

MIXC

Ake)D

—

+ MI),

~Hg)

—

—

1)/g],

--

G = (1 —MIXC

—

—

a (1)

1(5

H = (5 ~-~!i + MI ~ie)

2D(Ah

—

~Hg)

Ake) —

~H1(g2 1)/g], —

(4)

and 12zt

F

[2D(A~ Ake) + C(g2 —1 )/g2]

=

—

K

—

~

K[l —MI +

—

l)/g2],

(2)

a F d(113)

=~

a 1 —f3

(F— G)

—

a

K

(3

o —MIs

(3) where ~ = Z/Z~ is the with axial zdistance normalised to that of the nozzle throat = 0 at the nozzle inlet (fig.

=______

thhoAta~A~ where I is the current, ~h the mass flow rate passing the nozzle, h 0 the stagnation enthalpy, and A~the conductance shape factor [1]. The shape factors, Ah and Ake, are related to the primary standard shape factors (denoted by the subscript “1”) and the velocity magnitude factor through the following relations [1]: 2Ak~l+ 3g2(1 Ake g Ah = gA~ 1+ 1 g. —

—

1), a = A/As the nozzle area normalised to that of the throat A~,(3 = Od/A with °d being the arc displacement area [1], and. ‘y = CP/CV = 1.4 for air. In deriving eqs. (1)—(3), it has been assumed that the arc thermal

NOZZLe

.

z.’~”~’

I

I~entr~ion ExternaL

For a given solution parameter K and the nozzle shape, solutions to eqs. (1—3) are uniquely determined once the initial conditions,Mj, g1 and (3~are specified. Of the three initial conditions,g1 and (3~can be determined by considering the physical conditions at the upstream electrode [3]. The initial Mach number is in general unknown and is in turn determined by simultaneous satisfaction of the necessary and sufficient conditions for the existence of the critical flow [3]. These conditions state that:

flow

I—

and the compatibility condition Throat At

Mmimum effective flaw area

R(M,,,,,g,(3,~)O,

(5) (6)

where F is given by eq. (4) and R is one of the right Fig. 1. Model of the nozzle arc.

hand sides of eqs. (1—3).

M.T.C. Fang et al. I Computation of steady state arcsin nozzle flow

2. Method of solution A particular nozzle shape, a= + ~ ~

solution of eqs. (1—3) arises from the numerical matching of 0/0 (i.e. R/F) at the critical point. A slight mismatching leads to the breakdown of numerical solutions of eqs. (1—3) at a point where F = 0. This usually occurs because of inaccurate knowledge of M1 (i.e. coefficient r1) for a given value ofK and

(7\

~ ‘

‘

365

~‘

has been chosen for the present studies although the method of solution is general except that the solution behaviour near the nozzle inlet (~= 0) can be different [2]. = 0 is the stagnation point of the flow for the nozzle defined by eq. (7) since the nozzle area tends to infinity there. The series solution of eqs. (1—3) near the2 stagnation point is given by [2]: v.3

the accumulated computational error.This difficulty is overcome by using a ‘patching’ method, which ensures the exact matching of eqs. (5) and (6) at the critical point. The numerical procedures are: i) Guess the arcing conditions at the critical point, i.e. g, j3~,where subscript c indicates ii) and Calculate the critical Mach number the M critical from point; eq.

a—

(5) and the location of critical point ~ from eq. (6); iii) Calculate the Mach number at ~ (a very small

g

=

M

=

~ ~ + P2 sv. + + ..., +q 2+ ..., 1 ~ + q2~ +r 2 ~ + r3 ~ + ...,

(8)

number) through the continuity equation: 3 M..0ca~(l (3~)[l3k) + ~Ii + ~(‘y l)M11] l)M!~~j3, = M.~,jaj(1 ( using Newton—Raphson’s method, where (3, and g

wherer the coefficients are tabulated in table 1. Note that 1 is an arbitrary constant be and detetmined by simultaneous satisfaction of eqs.to(5) (6). The computational difficulty in finding the critical

—

—

—

—

1 are found from the senes solution of table 1; iv) Calculate the slopes of g, (3 and M0. at the critical point using L’Hospital’s rule [2]; v) A fourth-order Runge—Kutta method is used to Table 1. Series solutions for constant primary standard shape factors upstream electrode at the stagnation point.

: Pt

=

+p 2 3+ P1l~~P2l~ 3~ r~+r2~2+r3;3+...,

Ad~J~K/Ah; Ad and Ah are calculated at ~ = 0,

integrate backward from the critical point to a conveniently chosen ‘patching’ point (fig. 2) and then integrate forward to the same point from upstream electrode; vi) Formulate the square sum of the difference:

—

.

e(S~,j3~)= (MF MB)2 + ($F (3B)2 + (SF and minimise c using Powell’s method [4Jto improve —

—

—

the initial guesses ofg~and (3~.When e ~ 10~,a numerical solution is considered to be found.

3(3a+4bq 0) P2=—

=

q1

6(2a =

1—aq0

piqi,

(1 —aq0)(l+cq0) 3bq0)(1 + cq0) — c(1

+

—

aq0)

arbitrary constant

r2 =piri r3(p2+pj—1)ri+~(y—1)r~. __________________

wherea=1 —Adi, l~2Ad1—1, c=Ahl —1.

,—.—--.-

I

Ferwerd

mtagratà

::;;~—, :m~mtion I i,~[09.~I

0 Fig. 2. Thepatchingprogressandintegratjonscheme.

M. T.C. Fang et aL I Computation of steady state arcs in nozzle flow

366

3. Program description

4. List of subroutines and common variables

The detailed working of the program is described on comment cards within the deck. The main program includes the facility for performing a series of calculations of increasing arc currents (i.e. increasing K’), where the initial guesses of g~and (3~are the estimates obtained by examination of previous trends. A block diagram (fig. 3) is given for the main steps in the evaluation of the residual sum of the squares. The values of the primary standard shape factors Ahi and Akel and the conductance shape factor A~ can easily be altered by changing DATA statements in the program.

4.1. Subroutines

Start

(B

and g

are given) C

IPRIN

-

Calculate from

Calculate

.~‘

4

CR11 (BC, GC, ZC, RMC): evaluates the location of the critical position (ZC) and3~) Mach and No. GC (RMC) (i.e. ge).at theMONIT critical point from BC (i.e. i (N, X, F, T): routine required by NAG

IPRIN-1 and ~ and ~c

d

at Z c

YES

write

library routine for optional print-out of progress toward minimum: not used. Z~,lt c’ 8c’ d8

~

4.2. Common variables These are described in comment cards in the deck.

NO Calculate

EO4CAF (N, XYZ, RESID, E, ESCALE, ICON, W, 1W, CALC2, MONIT, IP, MI, IFAIL): NAG library routine that minimises a function of several variables by Powell’s method. CALC2 (N, XYZ, RESID): calculates the value of the function to be minimised i.e. sum of squares of differences of interactions to the patching point. INTEG (ZI, ZF, SOL, Z, DIVS): controls the integration from ZI to ZF. RUNGE (H, X, Y, E, DY): integrates the equations given by DY from X to X + H by a fourth-order Runge—Kutta method: The local truncation error is estimated in E. DYO (X, Y, DERIV): returns in DERIV the values of the derivatives of Y at X. VALUES (EMSQ, 0MB, G, Z, ALPHA, DALF, F, GEE, EZ, FMGB GRI): evaluates variables convenient in the calculations in DYO.

at Z~

Forward integration of the differential equations from Z~to the patching point

CARD 1 Backward integration ing point

: FORj4AT (13) contains value of

from Z~ to patch-

NMAX CARD 2 RNC

Form the sum of squares of differences Cr) of the two integrations at the patching point YES

: : :

Write N. (Z~), results

II’

CARDS (3—6) : XYZ (1), :

N RETURN

IFRIN

-

10

Fig. 3. The block diagram for the evaluation of residual sum of the squares called by EO4CAF.

XYZ (2) ZMID, ZMAX

NTs4AX number of cases to be attempted. FORIvIAT (F20.10). previous value of current parameter RNC = Ah/Kwhere Ah is defined at the nozzle inlet. FORMAT (2F20.10). previous values of(3 andg at the

critical point (corresponding to RNC) or guesses. : the patching point and the length of the nozzle.

M. T.C. Fang et aL I Computation ofsteady state arcs in nozzle flow

ORNC, DRNC : last but one value of current parameter and required change in current parameter. OX (1), OX (2) : values of(3 andg at the critical point corresponding to ORNC. CARD 7 : FORMAT (Li). LL : a logic variable indicating whether the calculation is one of a series. If it has the value TRUE then XYZ is treated as a guess for (RNC-DRNC). OX and ORNC are not used. However, if it is set to FALSE the calculation is assumed to be the resumption of a series with reducing RNC (i.e. increasing current) and the initial guesses are formed from XYZ, ox, RNC and ORNC.

367

(e) the residual sum of the squares at the patching point. When the minimisation has been satisfactorily cornpleted, CALC2 is called from the main program. In addition to (a)—(e) above, the values of Z, Mach number, (3 and g are printed at each step.

7. Description of test run The test run is for the case Ahi = Akel = 10 and A~= 0.4. The results shown are part of a series, where calculation has been restarted using the previously calculated g~and ~ to estimate the critical conditions for a new value of solution parameter K. The value of K for which the results are calculated is 8.1822 X 10—2.

Acknowledgement 6. Output

Only one output channel, a line-printer is used. For each current level, the current parameter, two values derived from it and the initial estimates of critical values of (3 and g are output by the main program. Then on the first and every tenth subsequent call of the routine CALC2, various values are output as a check on the progress and to facilitate easy restart in case of failure through lack of time. These are: (a) current estimates of critical position, Mach number,M,,,~,velocity magnitude factor,g~and (3~ (b) the corresponding derivatives at the critical point; (c) the Mach number at the tip of upstream electrode see (i.e. mitial Mach number); (d) the results of the forward and backward iritegrations at the patching point;

The authors thank Numerical Algorithms Group Ltd. (NAG) for kindly allowing us to include their library subroutine (EO4CAF) in our program listing. NAG wishes it to be known that this subroutine has now been withdrawn and has been replaced by a new version. References [1] M.D. Cowley, J. Phys. D7 (1974) 2218. [2] S.K. Chan, M.T.C. Fang and M.D. Cowley, Liverpool University, Arc Research Project No. ULAP-T44 (also Proc. of tEE mt. Couf. on Gas Discharges (pubhcation No. 143) p. 5—8). [3] S.K. Chan, M.D. Cowley and M.T.C. Fang, J. Phys. D9 (1976) 1085. [41M.J.D. Powell, Comp. J. 7 (1964) 155.

M. T.C. Fang et aL I Computation ofsteady state arcs in nozzle flow

368b

TEST RUN OUTPUT NO, fl~ CASES PATChING ~‘0TNT RFV 1 RNC

1

LC’GTCAL

.IUUODF,.O1

NOZZ~,j: LE~~CIrN

.1UUOOCaO3

OFLTA ~‘IC

‘

~

53)t)OE4O1

G

~ icr*~c,

r~c

,11~ME+IJ0

Zr

Mr

GC •141 ~‘4~.s;):j

,t7~Sri~+irt DI;

~ TN1~T.~L r’~fH

43~4U(

7

•c’i

UIJOUL+fl1 • I UUI)’~ •~I

• !

*

*

OF .4

*

*

*

*

*

Zr I ~ 1 t~ ~ + ‘~

*

*

.4

Tt1I~r?AL. r.~CH

*4*

*

4

4

*

=

?‘~~

.00

*

4

,

~/ ~

*

4

*

*

+

U0

37 nS “~— UI

— ,

~ II~I)F*i’I

0f~ ~I I

~qu4r’FS *

~ ‘~111 F+0(I

*

*

*

*

*

4

4

4

•

*

*

*

*

*

,

*

c.c

U 3~U ~ + )1

•1 ouon~eui

*

I

,

1 41 ~ 7 ~ +

, ~

.~~

•

Ni’,

)FtI~1) I’ ~rl~(’A~L’ (‘~

*

Mr. 1

•

,i ~?‘P4~~’) • I .U3~Ee1

(~,~4.’).I

‘,I~4.1

,1Ufv4?E—~)~

e~ii

cun

•Q1

S~1A~f~ t)~ ~ *

•

•4.(~r_U?

‘~.

Ir~Tl:C.P4TIL1I4:

~1Ir;

3~3~—~1

—,

*

*

*

4*

*

,1 ,~?Mi

*

4

*

4

*

*

~

I ~

~4.~7iF—O9 *4*4*4*

*

44

*

*

GC •

1 I~9Øf+~

851

TNITTAL PIACH NO.

,12211 E+~’

,01 (.,~4~IC) ,

zr .121 Z3f411~

•‘fl ~e~F#~T’)

V’ E+PO

,4u~t6r.OZ

•c7SSSE+0~)

—, 37)41 ~—0I

,141 ~7E.flfl

• 8386~E—01

*

*

*

*

*

*4

*

*

M. T.C. Fang et aL I Computation of steady state arcs in nozzle flow 1NTE(PATt0~

Z 1 0000C+01 ,lu000r.01

FORIJARO PACKIJARD

Stir “r *

*

*

*

U ,91 (‘~6E.ufl ,91(’~E+OC)

•

SQUAR~S OF DIFhEP~FJCES *

*

4*4

*

ZC ,‘2123E+O1

*

*

*

*

*

*

META

*

,b31 29E.00

~ORWA~D INTFGWATIO”J

,1~2~2E.00

*

*

*

*

*

*

—

3/)41 E—01

,t)”~47E’ Z 081 ~

•bSdllEOl ••13136F.0O Cl ~ 33!. 00 ~94~L.00 ,3f~4nr.00 ,4~~2O~i00 • set s7r.0O • 60 v 2 ~F~* 00

,‘~~‘Eø0.? ,,1~S1~—Q1 301 ~4 F ~ 01 ,‘.1191~—Q1 •‘)1~4’.~~ØI ,a~12.~.OI ,7~?01OE—01 • 81 50 ; E —01

,d4562

~ 3 i’

U

1,

•

f3b~1 r.OQ

•dUOIS~lr.00 • 81 #, 5 ~ * 00 •91 ~M~e00

i1049M

M

•

*

*

141~?E+flO

• ~

rIb

•

7371 IF .~‘(

74’TstU ••fl7 PJ? ~ .5 E + ~ 0

•70#14~.CO •7045AE+N)

•7026I~i’Cfl • 700 33~.(.’i’ ,69? ~‘ E + (0

•90(.03c.OI ,9’~31~i~~OI ,1076o~.0O • 11 ~ ~‘.E +00 • I ~3~3~*C)fl

RACKWARD taJTEGRAT1r)I~ FPm’

Z

*

1

,1e~33PE~0’

•8899”i.00

*

rf’y.

• ‘,Q(’2?4~—0~

• ‘P65~4

*

To !)%TCIIl~f prI~’~’

.1108u1_01

065 ~

*

r.c

• ~u.~$~6F—0~

~‘

*

•o7~SSE*0fl

‘l~4fl(’ ,O10~7

,36754 ,4472~

6

, 1 22~1F+()Q

,36?45E—10

MC ,lu)9flEsOI

M

369

CPIT,rSL

![Tj

1,1998’ 1,1Q??n

,1u)9”r,Ol ,1V~11E,01 ,1u4.93F.O1

•14167E,00 ,160b3E+aD •14045+00

1,1934~

,1u~7,1!,0I

,14D09~,00

1,1~4QA 1,1~79?

,1U~16!,0I ,10307r.O1

,I393’P.(~O ,I3?O1EeO0

,IV1Q~,!*0I ,9V(1B~.00

,13t’4S~,0Q ,1334~E.00

1,08305

•9f~2?c#00

,13fl4AE.oQ

1,0490? ~ ,01 Sin c 81 i 3

,9)091!.O0 ,9efO?.O0 • 9 1) ~75 ~* 00

,1~P~’?E.O0 • 1 ~4~4Ee00 • I ~1 06 ~ * C) U

, 69??SE.(Q • e ~ ~ 3 .~ + 00 • 6.1331 F.

PCI~4T

C

,67609F.flo ,676?4L.rvo ,6P6’i6E.0Q ,6??IQE,OC)

68C’37E+r~u 6M1f,SE.00 • , 68441 E + 0 u

*

4

*

*

*

*

*

M. T. C. Fang et al. I Computation of steady state arcsin nozzle flow

370

jRWARDjPJTE’CJPATION

£

M

1,21~3? I,~48? I ,24416 1,34?i~

•~U~90E*Q1 •1066”C.01 • I tI(8~!*0i ,i14o.~reo1

1,45017

•iVv94~*O1

1,5531h

,1~)4?E.OI

I , ~‘,56I9 1,75921

•

1,96526 ~,O68Z5 ~,17127 £,37?29 ~,4803’ INITIAL

MACH NO,

FI~0’1 CRITICIL

POIrJT

C

rET..~ ,14I6PE~OrI ,14272Ee00 • I 4432E+O0 •lS2e’OE,OO ,1604SE,00 ,I6791r 4Qfl • I ?S0I ~.O()

•675S~E,D~~

,I8i7~E,o0

,656SU*0u

•14U57F,01 •14485F.0I

•i8~424Ee00 ,19441E.0Q

•653~~F+i’0 •6501!E+(1o

•14’~i2E*QI

,~U~33EeoO

,i!’.~20E*01 •1~f1oE*01 •IAU~3E*0I •18441F*0I

~ ,~1144E+Q) ,I~1667~è0) ,22i71!4i~~

,64414E’nO ,64120E+00 •6585~E.C)0 ,635R?E+OQ

M 9I’~26E.~0 ,~?I A26E.Q0 , ,36236E.ilO

RET~ ~ • I 2~81E+00

I 3V70E,0I ,13~66~.O1

•675~gE+oo • 87436E’i~’O

•67060E,U0

,

,4u5,36r~o2

INTEr~RATION~ FORWARD , I 0000E+01 RACKI’A~D ,10000E.01 SUN OF squAprS O~ DtF~EPENCES

6 -

A*~5nFOC)