Axisymmetric accretion flows in astrophysics

2/9

computer rnysics commumcauons ‘i’+ ~ivo F) Z i~—zoo North-Holland, Amsterdam

AXISYMMETRIC ACCRETION FLOw~ir~

iiiix~i~

J.A. ROBERTSON Department of Applied Mathematics, The University, St. Andrews, Fife, Scotland

This paper describes a code that has been written to simulate viscous accretion flows that possess angular momernum. .‘tn explicit Eulerian flux-corrected-transport (FCT) algorithm is used to perform advection of the fluid, which is a perfect gas. An alternating-directions-implicit (ADI) scheme is used to calculate viscous diffusion. The viscous dissipation of energy can be included if desired. The code has been used to model the~viscous boundary layer between a thin Keplerian accretion disc and

I. introaucuon

Accretion processes constitute one of the most exciting areas of study in modern astrophysics, and provide some of the most challenging problems for observers and theoreticians alike. Accretion, arising as it does from gravitational attraction, is likely to be of importance on every astronomical scale, from the formation of planets in the early Solar System to the fuelling of the most luminous known objects, such as quasars and active galactic nuclei. True, direct observational evidence for accretion of material is at present scarce. Most information is to be had from the study of close binary stars, and many detailed observations are available for binaries in which the accreting star is a white dwarf. One star may lose material in the form of a stellar wind or a stream through the inner Lagrangian point. The companion star may then capture this matter directly, or the gas may settle first into a thin equatorial “accretion disc”, through which it is slowly fed to the central star. Any magnetic fields that are present will exercise a strong influence on the manner in which accretion occurs. In the case of quasars and active galactic nuclei, it has been suggested that accretion occurs via thick, pressure-supported tori orbiting the compact central source. Various aspects of accretion are reviewed by Rees [1], Warner [2], Pnngle [3], Begelman, Blandford and Rees [4], Frank, King and Raine [5], Lewm and

Van den Heuvel [6] and Hayakawa [7]. The dynamical stability of thin discs and thick tori is currently being actively investigated. Although it seems probable that they are all linearly unstable to non-axisymmetric modes [8,9], the consequences of such instability may be relatively harmless. It is possible that these linear instabilities may give rise to a saturated regime in which their effect may best be regarded as a sort of viscosity. The high level of viscosity that is known to be present in thin discs around white dwarfs has yet to be satisfactorily accounted for. The code to be described in section 3 below was written to investigate axisymmetric accretion of matter from a disc or torus by a star. One aim was to study the structure and possible time-dependent behaviour of the viscous boundary layer that forms at the inner edge of a thin accretion disc. Of interest in this context are the extent to which the boundary layer spreads in latitude over the stellar surface [10] and the efficiency with which the material of the disc mixes with the outer stellar layers [11]. These points are of importance for studies of the explosions of classical novae. The spectrum of the radiation emitted by the boundary layer, where half of the total accretion energy is dissipated, depends on the layer’s structure, but this in turn depends on the cooling [12]. A complete treatment would require simultaneous solution of the equations of hydrodynamics and radiative transfer. To date, only the hydrody-

~c t.isevier science ruousners is.v. North-Holland Physics Publishing Division)

JU1U-’fO)D/o//bUi.DU

namics have been implemented. Transfer of tion in the diffusion approximation could corporated without difficulty, being math~ cally similar to the diffusion of angular m~ t~im described in section 3.2.

me numerical methods

~

artificial viscosity difference scheme tinuities. Thus, a s~.. curate where the flow is smooth may first-order accurate at a shock. When variable is updated by applying the L fluxes and source terms, a strong diffusion is also applied to that variable, thereby ensuring that the solution is smooth. To the result of this step is --

would be impossible to survey in a short e all the available techniques for numerical odynamics, and the brief discussion here is tly confined to methods that have been used udy accretion flows.

1~-~’’~~J

oped an axisymmetric, general relativistic code to study inviscid accretion flows in the gravitational field of a Kerr black hole. They compare a variety of finite-difference algorithms for ideal hydrodynamics, with particular emphasis on methods that impose a monotomcity condition on advected quantities. Such a condition ensures that quantities that increase or decrease monotomcally before advection also possess this property after advection. Their own technique is based on that of Van Leer [15—18],with the addition of artificial viscosity [19] to enable accurate treatment of shocks. They also discuss advective schemes devised by Wilson [20] and Barton (Centrella and Wilson [21]), and a procedure to improve local conservation of energy and angular momentum [22]. Ró~ycka[23] describes a code written in cylindrical geometry to study the expansion of stellarwind bubbles into the interstellar medium. The code is based on the thesis by Norman [24], and incorporates artificial viscosity and a monotonic transport algorithm derived from the work of Van Leer. Clarke, Karpik and Henriksen [25] present non-relativistic calculations of the evolution of an approximately spherically-symmetric flow (Bondi accretion perturbed by the addition of angular momentum) into a thick disc. They use a flux-corrected-transport (FCT) scheme based on the algorithm SHASTA devised by Boris and Book [26].

-

WIICIV

I.1L~ aILIUUIII

UI

aIILi—UIIIUMUI1

I~ I1ILUL~U LU a

level that ensures monotomc behaviour. The algorithm ETBFCT has been used by Wang and Robertson [31] to study the Rayleigh—Taylor instability. Colella and Woodward [32] have devised a higher-order extension, the Piecewise Parabolic Method (PPM), of Godunov’s method (see also Van Leer [18]), which allows very sharp representation of discontinuities, in particular, contact discontinuities. Extra dissipation, e.g. artificial viscosity, is required by this method. Woodward and Colella [33] compare in detail the performance of a variety of schemes, with emphasis on their ability to model flows with strong shocks in two dimensions. The algorithm ETBFCT is one of those considered. 2.2. Particle and particle—grid schemes Particle methods for hydrodynamics, being Lagrangian, have the advantage that local conservation of quantities such as mass and angular momentum is automatically guaranteed. They have recently been reviewed by Eastwood [40] and Monaghan [41]. A finite-size-particle scheme for 2D and 3D problems was devised by Lucy [34], who used the equations describing a self-gravitating compressible gas to study the fission of proto-stars. His method has since been used by Gingold and

J .11. flUUfl

&JUfl

/

Monaghan [35] to study rotating polytropic by Wood [36] to study collapsing, isotherm clouds, and by Noithemus and Katz [37] to the effect on a star of a close encounter ‘‘~khole. Recently, Anzer, Börner and Mona~ian have used this technique to make a study in )f the accretion by a neutron star of the wind a hot companion star. he obvious disadvantage of pure particle iods is their high cost, involving as they do actions between every pair of particles. The nal Particle-in-Cell (PlC) scheme devised by ow (see, e.g., ref. [39]) combined the adaaec of Pulenan and Laaranaian methods by

garueu as unsui~aoic ior nyurouynaiiucs. 1w.c111 developments, however, have allowed improvements in the basic PlC scheme. Optimized versions of PlC are discussed by Eastwood [40]. He shows that these new Ephemeral Particle-in-Cell (EPIC) schemes, which can be used for ideal and non-ideal hydrodynamics and MHD, are able to out-perform state-of-the-art pure finite-difference methods. They combine the best features of FCT, PPM and particle methods.

.LOI

i~uuloulculurn,

2

(ph) = r2 sin2O tan 0 ap/ao, P whence, using the radial equation,

~

,~_

,(,~3 ~

~

—

r2/(A,f

(~

cornigurauons jor uisvs anu Loll, wiui LUC proviso that the self-gravity of both disc and atmosphere is negligible. The equations of the problem may be written in the following form. Continuity equation:

ap/at+v•(pv)=o.

(3)

Radial momentum equation:

3(pvr)/at+ v•(pv,v) = —ap/ar—GMp/r2÷ ~7i. +p(v+v~~)/r—(T

3. Adopted numerical method One of the most interesting questions is how the material at the inner edge of a thin accretion disc mixes with the outer stellar layers. It is therefore convenient to perform computations in a spherical shell, and tp use polar coordinates (r, 9, 4). The inner boundary may then be chosen at a level within the star to which the disc material cannot penetrate, and which therefore remains essentially undisturbed. It is also convenient to assume that the stellar atmosphere, which has a very small pressure scale-height, fills the computational domain. One thereby avoids the difficulties of tracking the boundary of the disc and of treating regions of vacuum. The initial state for a calculation is one of pure rotation. One first chooses a spherically symmetric distribution of pressure, pa(r), representing the

is,

nents of the invis~ 6-equation, assumii the quantity

09+r~)/r.

(4)

9-momentum equation:

a(pv9)/ar+ =

ç”.(pv9v)

—(ap/ao)/r÷r.~ +p +

( v~cot 6

(~ —

—

VrVO )/r

cot 0)/r.

(5)

Angular momentum equation: a(ph)/at+v’.(phv)=c”.(rsin0T~).

(6)

Internal energy equation:

ap/at+v•(~)=---(y—1)pv’v+(y—l)& (7) In these equations, p, p, v,1,, h = v4,r sin 0, y = and ~‘ are, respectively, the pressure, den-

sity, toroidal velocity, specific angular mome ratio of specific heats and dissipation fur The poloidal velocity, v, is V = vrer + v

9e9

the viscous stress tensor, T, is (9)

~er+ ~e9+ r~e,, ~e,e.g., rrre,.

+

‘r,.9e9

+

‘r~e..

(10)

Sessions for the dissipation function, and omponents of r (in terms of v, v~,v.,. and the ~,

U.ILLIJI

ALL ¼A.JILA}J’JOLhUh.

Evolution of eqs. (3) to (7) from the initial state of rotational equilibrium is carried out using Eulerian finite-difference methods on a non-uniform grid of N~by N9 cells. The complexity of the equations makes it impracticable to treat them as they stand. In particular, disparity between the time-scales for advection and viscous diffusion forces one to consider explicit treatment of the former but implicit treatment of the latter. It is therefore convenient to regard eqs. (3) to (7) as a combination of three distinct processes which may be carried out consecutively, namely advection, diffusion and dissipation, and to split each timestep, in which the solution is advanced from time t to t + ~t, into three parts. This is the same strategy as is adopted by Weber, Boris and Gardner [43] for their resistive MHD code. A minor disadvantage of this procedure is that there is no clear measure of the order of accuracy of the numerical scheme in time. But a numerical method is better judged by its performance in any tests that can be devised than by the ease with which it submits to numerical analysis. What might be regarded as a major disadvantage is the difficulty of posing consistent boundary conditions when the equations are split up. (The advective step involves hyperbolic equations, the diffusive step, parabolic equations.). But consistent treatment is in any case extremely complex even if one neglects practical, numerical considerations, and is

numerical treatmer The great adva adaptability. It is

~

processes, such as thermal diffusion, or duce alternative algorithms for any part. The separate steps are now described. 3.1. Advection

are omiueu irom eqs. (Si) tO ~ /). rossioie scnemes for this step have been mentioned in section 2. The adopted scheme is based on the Eulerian version of the one-dimensional, explicit, conservative, flux-corrected-transport algorithm, ETBFCT, of Boris [301. Because this is a 1D algorithm, the advective step must itself be split into two 1D substeps, the algorithm being first applied along every line of fluid cells in turn in the radial direction, and then in the 0-direction. (The algorithm may be used along any coordinate line of any orthogonal coordinate system.) The splitting of the source terms that appear on the right-hand sides of eqs. (3) to (7) is rather arbitrary, but 0-derivatives must be avoided in the radial substep, and vice versa. Thus, each of eqs. (3) to (7), which have the form

af/at+v•(fv)=s=sr+so,

(ii)

is split into a radial equation, 3f/8t +

—~

-~--(rfVr)

=

S~,

and a 0-equation, 1 af/at+ r 9-~(sin9fv9)=S9.

(12a)

(12b)

To obtain second-order accuracy in time for the radial and 0-substeps separately, a simple predictor—corrector scheme is used:

J.I1.

flUUC~ (hUll

/

25~

fl

Radial substep: (1) calculate sources at time t”. (2) advect every equation of the form (12~ 2= t” + 8t/2. Symbolic from t~ to t”~’~ 8tF1 ~~~(r2f(0~0)) a ~(1)f(O)± + S,(0)].

~1

~alculate sources advect in r from

=

j(O) + ~t

{

t”

at t’~”2. to t’~ = t” + 8t:

Interface of area B

}.

3_1 Fig. 1. Notation for the algorithm for the advei

2f(O)v,(1)) + S~1~ -~

~—

(r

-substep:

(4) apply diffusive flux (unphysical! Suppresses

f~=f~~ ~{r

s~i~ U(sin 0f~2~v~2~)

(5) calculate anti-diffusive flux: A

(ftft

—

j~1/2l5j—1/2~Jj

+s~2)J.

ij—i

(6) apply anti-diffusive flux: (7) calculate sources S 3~. (8) advect in 0 from t”9~to t” ~‘: s~n0 ~(sin

f(4)f(2)+8t[

+

~n+1_~td

if of(2)v~3))

j±1/2

\Jit j—1/21/

j~

step (6) almost cancels step (4). 13

2

(1) calculate advective flux: (fflL

j—1/2~j—1/2~Jj

In

ij—i

(2) apply advective flux and sources (physical!): jt jn + (~-1/2 1~+l/2)/’~f + atS~. —

(3) calculate diffusive flux: ~-1/2(ff

~ t’j_1/2,

~

—

As discussed in section 2, FCT algorithms apply strong diffusion everywhere to inhibit unphysical oscillations of the solution. This diffusion is then removed by applying an anti-diffusion whose magnitude is very similar to that of the diffusion except in the immediate vicinity of flow discontinuities, where the anti-diffusion is reduced to a level that ensures that monotonicity is preserved. The basic algorithm is thus (see fig. 1 for notation):

=

~Jf

If The choice

S~3)j.

L’ _1~ j—l/2 2

A

(A

—f~).

P’f—1/2

—

6

3Ej_i/2,

1

12

6

6~j—1/2’

where =

~tBf_1/2Vf_1/2(T

(is)

+

is particularly favourable in that it minimizes the dispersive error. Between steps (5) and (6), the anti-diffusive flux must be corrected to provide accurate treatment at discontinuities. This is obtained as follows: (5 ‘~

\ I

(Sb)

.....~td.....jtd

A

j—l/2Jj ‘j-1/2

(

ij—i’

1, ~ —1,

if if

~f_1/2<0.

(Sc) corrected anti-diffusive flux: A~_l/ 2=I~_l/2max(0,min( IAJ_l/21, ~_l/2Af~f±l/2,

~_l/2Af_l~J_3/2

}).

.//a.

~ooenson /

/

Step (Sc) has an undesirable side-effect. It sentially a local adjustment of the amot anti-diffusion that may safely be applied adjacent grid-cells experience unrelated ---s. The result is that regions through winch a ontinuity has passed tend to have ragged proof the fluid variables (cf. fig. 4 in Robertson Frank [45] with fig. 12 in Hawley, Smarr and ;on [14] who do not use an FCT scheme). It is that a more sophisticated flux-correction is rable if algorithm smooth results are to in If aiD is applied in be 2D,obtained unphysical [lations will almost inevitably result. See sak [461for a fully 2D FCT algorithm.

ucpal LU1~ 110111 LIII.

the current values cells, so that only the r-scan, one thc~j1

i~a~

of the form + 1/2 ~

1~f—1,k’ vfk,

V

~~j—1.k, all U1 Uf+lk, G,~(uat 9 k’ cells, v at all 9 cells),

=

where u

=

v,., v

=

(20)

v 9, and F’~±1/2 and G,~ are

1’ F’~ mouuy tue velocity field, but the variables p and p remain unaltered. The equations are In

UU5

9

sLep, viscosity acts to

p avr/at= 甕i.— (‘r99+’r,~)/r,

(16)

p 3v9/at = v.r9

(17)

+ (T,~—

cot 0)/r

~U

12f,k—1’ 1~Jk~ ~jk+’

3,j~_p =

Ulk, U~k±1’

G~”2(u at all 9 cells,

v

at all 9 cells).

(21)

Thus, in each scan, one has a 2 X 2 block-tridiagonal system of equations which can readily be solved by a generalization of the standard technique for a tridiagonal system.

and prsin08v,/8t=V(rsin0T~).

(18)

3.3. Dissipation Since ~ depends only on v~,eq. (18) may be solved independently. Eqs. (16) and (17) are coupled, both involving yr and v9, but not v~.If the coefficients of viscosity were very small, an explicit scheme would be satisfactory. In general, however, an implicit scheme is essential to ensure that the time-step is limited only by the Courant condition for the advective step and not by the diffusive step. The ADI scheme of Peaceman and Rachford [47] is satisfactory. For the angular momentum equation (18), one first integrates over the volume, A, of a grid-cell, obtaining

If the kinetic energy that is dissipated by viscosity all goes into heating the gas, the change in gas pressure is

Apr sin 0 av~/at=

4. Tests

sin 0

T~ dS,

(19)

where the surface integral is over the entire cell boundary. From eq. (19) one can easily obtain ADI finite-difference equations that conserve angular momentum. For both the r and 0 scans of the ADI scheme, reduction of a tridiagonal system

=

~

—

1)&

(22)

The new value of the pressure is obtained from ~

(23)

where and ~ 1 are the values of the dissipation function before and after the diffusive step. If so desired, one can alter eq. (23) to allow some or all of the dissipated energy to be radiated away.

There is a lack of known solutions for compressible viscous flow in a spherical shell. It is therefore not possible to check the performance of the complete code. Fortunately, however, there are

285

Fig. 2. Contours of density for an initial equilibrium configuration of disc and atmosphere.

Fig. 3. Evolved state of the configuration in fig. 2, showing a hin Keplerian accretion disc and boundary layer spreading over the stellar surface. Viscous dissipation has been omitted.

~ ‘.

(1% ~ ‘I

I~

c.,c’c’~~

U

~ ~ ~

~

-~--~-~

________

_________ _______________

-

—

—

-~ .

— — __..—...——.——..—-—..—-.-—~s —. —~--_--,--*--,.—-~.

—~r

some solutions for incompressible flow, an provide valuable checks of most of the cod A simple test that can be used to cc algorithms for advection is provided by th sive advection of a scalar function, f, by )ed velocity field: )t + V~vf= 0. (24) --

~

examples, is also Wilson [13,14]. Th cal geometry does but an accurate numerical solution coull structed and used to provide a more check of the advective step. Diffusion of angular momentum is readil) tested by following the evolution to a steady state of toroidal flow between two spheres rotating al ~

ather severe test of this form is described by ~sak[46]. ‘he non-linearity of the equations of hydrodytics means that such a test only provides a er limited idea of the value of an advective rithm. A more complete test is provided by a

/

‘ ~

,

-I

?t?tlttt,,,,

~

I---

/ /

/

ttt?tttt,,t,t,,_~

~ —---

—---

tttt%\\\~tt~tt/~Y

tttttttttt~tt

t

/

—

/

/

—

/

—

/

/

/

~

/

/

/

/

—

—

~ ..-‘

-

—I. — —~ —‘ ~

/

/

---s ~

~

~

~

••~~~ —I.

~

—~

,

/

— —

,

-.

—

-

, —

/ /

—~ —

/

/

/ /

—

•—~ ~

-~

~nmThThEfl9 Fig. 5. Enlargement of fig. 4, showin

— I—

—.

/ / ~ / / / ,~ / ,//

~ ~

—

, 7,

— ~

, /

/

,

/ / // — -.

-~- —

—

,

/

/ /

-~

-~

~

/

~

—-----.-, -----.--— — —.

—~

, /

/

/—

~

/

— —

,

/ /

—,

‘~

-,

/ /

, /

-

-~ —

,

,

1

— ,

/

/

— —

-

, ,

/

I -

,,‘—Il

/ /

I

-.

—

.—

/

/, /, / /

/ —

,

/

/ /

—-—

,

/

/

—

—

‘

,

/

/

—

.—

ftftttHt,,tt tttttittttItt t,tttttt,tttt

I

/

I

/

--

uls at the inner boundary of the disc.

-.~

, -I,.

J .11. iwoerison / /a

287

provided by the following problem. The s equations in cylindrical coordinates for the’ ity (CR’ v5) of an incompressible viscous po flow are —

=

vR/R

=

0

0.

solution of these equations is 2J k 1(kR) cosh kz,

(25)

(26)

(27)

ciiirusion may oe testea simuitaneousiy oy lrnposing the solution (27), (28) on the spherical boundaries of the domain and v =0 in the interior, and evolving to a steady state. If the specificheat ratio, y, is chosen sufficiently large, the resulting flow is almost incompressible. Finally, one may consider incompressible viscous flow between concentric rotating spheres, a problem that has been studied by Pearson [49] and Munson and Joseph [SO]. Pearson provides numerical solutions of the time-dependent behaviour following an impulsive start, while Munson and Joseph give approximate analytical solutions for the steady state.

5. illustrative calculations Fig. 2 shows logarithmically-spaced density contours for an initial equilibrium configuration formed from an isothermal atmosphere and a rotating ring. In fig. 3, density contours are shown at a time when the ring has evolved into a thin disc which has reached the star and spread up to latitude 60°, hugging the stellar surface. In this case all the dissipated energy has been radiated away. Figs. 4 and S show the velocity field at the same time as fig. 3, but for the case in which dissipated energy is entirely absorbed locally. For further details, see ref. [45].

institut rur i~strop with financial sui Gesellschaft. Com~ on the Cray computer at the ~..omputer Garching. I am most grateful to Yi-Mit for introducing me to FCT algorithms.

References

Mitton and J. Whelan (Reidel, Dordrecht, 1976) p. 85. [3] J.E. Pringle, Ann. Rev. Astron. Astrophys. 19 (1981) 137. [4] M.C. Begelman, R.D. Blandford and MJ. Rees, Rev. Mod. Phys. 56 (1984) 255. [5] J. Frank, A.R. King and DJ. Raine, Accretion Power in Astrophysics (Cambridge Univ. Press, Cambridge, 1985). [6] W. Lewin and E. van den Heuvel, eds., Accretion-driven Stellar X-ray Sources (Cambridge Univ. Press, Cambridge, 1983). [7] S. Hayakawa, Phys. Rep. 121 (1985) 317. [8]J. Papaloizou and J.E. Pringle, Mon. Not. Roy. Astron. Soc. 208 (1984) 721. [9] J. Papaloizou and J.E. Pringle, Mon. Not. Roy. Astron. Soc. 213 (1985) 799. [10]R. Kippenhahn and H.-C. Thomas, Astron. Astrophys. 63 (1978) 265. [11] J. MacDonald, Astrophys. J. 273 (1983) 289. [12] A.R. King and G. Shaviv, Nature 308 (1984) 519. [13] J.F. Hawley, L.L. Smarr and J.R. Wilson, Astrophys. J. 277 (1984) 519. [14] J.F. Hawley, L.L. Smarr and JR. Wilson, Astrophys. J. Suppi. Series 55 (1984) 211. [15] B. van Leer, J. Comput. Phys. 14 (1974) 361. [16] B. van Leer, J. Comput. Phys. 23 (1977) 263. [17] B. van Leer, J. Comput. Phys. 23 (1977) 276. [18] B. van Leer, J. Comput. Phys. 32 (1979) 101. [19] J. von Neumann and R.D. Richtmyer, J. AppI. Phys. 21 (1950) 232. [20] J.R. Wilson, in: Sources of Gravitational Radiation, ed. L. Smarr (Cambridge Univ. Press, Cambridge, 1978) p. 423. [21] J. Centrella and J.R. Wilson, Astrophys. J. Suppl. Series 54 (1984) 229. [22] M.L. Norman, J.R. Wilson and R.T. Barton, Astrophys. J. 239 (1980) 968. [23] M. R6±ycka,Astron. Astrophys. 143 (1985) 59. [24] M.L. Norman, Ph.D. Thesis, University of California, Livermore, USA (1980).

~oo

J../i.

twoenson /

[25] D. Clarke, S. Karpik and R.N. Henriksen, Astr Suppl. Series 58 (1985) 81. [26]J.P. Boris and D.L. Book, J. Comput. Phys. 11 ( [27]D.L Book, J.P. Boris and K. Ham, J. Comput. (1975) 248. J.P. Boris and D.L. Book, J. Comput. Phys. 20 (1976) 397. J.P. Boris, Comput. Phys. Commun. 12 (1976) 67. J.P. Boris, Naval Research Lab. Memo. Report 3237 (1976). Y.-M. Wang and J.A. Robertson, Astrophys. J. 299 (1985) 85. P. Colella and P.R. Woodward, J. Comput. Phys. 54 (1984) 174. P.R. Woodward and P. Colella, J. Comput. Phys. 54 (1984) 115. L.B. Lucy, Astron. J. 82 (1977) 1013.

[42]J.-L. Tassoul, Th Press, Princeton,: [43]WJ. Weber, J.P. Conunun. 16 (197 [44]Y.Q. Hu and S.T. Wu, J. Comput. Phys. 55 (: [45]J.A. Robertson and J. Frank, Mon. Not. F Soc. 221 (1986) 279. [46]ST. Zalesak, J. Comput. Phys. 31 (1979) 335. [47]D.W. Peaceman and H.H. Rachford, J. Soc. Ind. App] Math. 3 (1955) 28. [48]G.A. Sod, J. Comput. Phys. 27 (1978) 1.