Kinetic-Fluid Coupling in the Field of the Atomic Vapor Laser Isotopic Separation: Numerical Results in the Case of a Monospecies Perfect Gas DELLACHERIE Stéphane Commissariat à l’Énergie Atomique 91191 Gif sur Yvette, France E-mail: [email protected] Abstract. To describe the uranium gas expansion in the field of the Atomic Vapor Laser Isotopic Separation (AVLIS; SILVA in french) with a reasonable CPU time, we have to couple the resolution of the Boltzmann equation with the resolution of the Euler system. The resolution of the Euler system uses a kinetic scheme and the boundary condition at the kinetic - fluid interface – which defines the boundary between the Boltzmann area and the Euler area – is defined with the positive and negative half fluxes of the kinetic scheme. Moreover, in order to take into account the effect of the Knudsen layer through the resolution of the Euler system, we propose to use a Marshak condition to asymptoticaly match the Euler area with the uranium source. Numerical results show an excellent agreement between the results obtained with and without kinetic - fluid coupling.

July 2002

INTRODUCTION The aim of the Atomic Vapor Laser Isotopic Separation (AVLIS) is to separate uranium235 from uranium-238 to obtain the fuel for nuclear plants (cf. [1]). From this point of view, the AVLIS process vaporizes uranium by using an intense electronic beam which heats an uranium liquid source up to 3000 Kelvin (the uranium output is of some kilogrammes per hour). Then, the uranium vapor is irradiated by a laser beam which ionizes the uranium-235 (and, ideally, not the uranium-238) further collected as a liquid on collectors which are negative electrods. To simulate the uranium expansion, we have to evaluate the distribution function f x v of the uranium expansion which is the stationary solution of the classical mono-species Boltzmann equation




1 Q f f (1) ε because the gas is almost rarefied (in (1), ε defines the order of the mean free path). Let us note that some papers have already focuss on the simulation of this Boltzmann equation for AVLIS applications: see [2], [3], [4] and [5]. Nevertheless, near the source of uranium, it exists a tiny area where the vapor is very dense and near the thermodynamic

∂t f



v ∇x f 








CP663, Rarefied Gas Dynamics: 23rd International Symposium, edited by A. D. Ketsdever and E. P. Muntz © 2003 American Institute of Physics 0-7354-0124-1/03/$20.00

947

equilibrium (then ε  1) which makes the CPU time of the Monte-Carlo technics used to discretize the Boltzmann equation dramaticaly increase. To diminish this CPU time, we discretize the fluid limit of the Boltzmann equation in the dense area where the gas is at the thermodynamic equilibrium – i.e. in the fluid area –, limit which is the Euler system ∂t   ∇x      0 (2) closed with the equation of perfect gas (with γ  5  3), knowing that   ρ
ρ u
ρ E  and that     ρ u
ρ u u  P
ρ E  P  u  by using the classical notations. Of course, in the remaining part of the physical domain – i.e. in the rarefied or kinetic area –, we solve the Boltzmann equation. In other words, we have to solve a domain decomposition problem which is named here kinetic - fluid coupling. Moreover, between the uranium source and the fluid area, the gas is not at the thermodynamic equilibrium although it is very dense: this very tiny area is called Knudsen layer (see the figure 1). To optimize the gain in CPU time, we would like to asymptoticaly match the fluid area where the Euler system is solved with the uranium source. The technics chosen to obtain a boundary condition between the kinetic and fluid domains was previously used in [6] for classical aerodynamics problems (see also [7] and the references herein). This technics uses a kinetic scheme (cf. [8]) in the fluid domain to discretize the Euler system since it allows to define a natural boundary condition at the kinetic - fluid interface with no overlaping between the kinetic and the fluid domains. To define the boundary condition for the asymptotic matching of the fluid domain with the uranium source, we use a Marshak condition which, when it is coupled with the previous kinetic scheme, defines a condition similar to the previous one applied to the kinetic - fluid interface. This Marshak condition is well adapted since, for example, it does not suppose the value of the Mach number at the exit of the Knudsen layer as it should be done if we applied a Dirichlet condition deduced from the study of the half space problem (cf. [12], [13] and [14]). The plan of this paper is the following: in the following section, we recall the basic properties of the kinetic schemes. Then, we describe the kinetic - fluid coupling algorithm, the boundary condition at the kinetic - fluid interface and the Marshak condition designed for the asymptotic matching of the fluid area with the uranium source. In the last section, we present numerical results. Let us note that to simplify the notations, the algorithms are written in monodimensional cartesian geometry although the numerical results are obtained for a bidimensional axisymmetrical geometry (see also [11]).

KINETIC SCHEME Let us define the monodimensional conservative scheme  i

n 1



ρin 

∆t ∆x 

n i 1 2 



n i  1 2 

(3)

for the Euler system (2) (i and n are respectively the space and time subscripts). The kinetic schemes are constructed from the following lemma due to B. Perthame (cf. [8]):

948

Lemma (B. Perthame) Let us define the initial conditions ρ 0
x  , u 0
x  and E 0
x  of the Euler system (2) which are supposed to be regular and let us define the function

χ vx  0 such that



1
v2x  χ vx  dvx

IR



1
1  and χ

 vx  χ vx  . Let h t
x
vx  be

solution of the monodimensional transport equation



∂t h  vx ∂x h  0
h t  0
x
vx   M ρ
u
P  vx 

(4)

u2 3 P    χ   vx  u  0  x  with P such that E    .       P  0  x  ρ  0  x  2 2 ρ Then ρ t
x  , u t
x  and E t
x  defined by !!! ! ρ t
x $% h t
x
vx  dvx
IR " !!! ρ t
x  u t
x &' h t
x
vx  vx dvx
!# IR ρ t
x  E t
x (  h t
x
vx *) v2x  2  P t
x + ρ t
x -, dvx

where M ρ
u
P  vx 



ρ 0x P 0x ρ 0x

IR

is an approximation in ∆t 2 of the solution of (2) when t

. ∆t (in 1D cartesian geometry).

Then, by using an upwind scheme to solve (4) and by taking χ vx / 0 we obtain a first order numerical scheme defined by the numerical flux

1 2π

exp

 i 1 2 ' i  1 2 1 i  1 2

2  vx  2  ,

(5)

where the positive and negative half fluxes are defined by

 i  1 2 ρi
ui
Ei    vx 43 vx 2 0

1 vx v2x  2  Pi  ρi

56

M ρi
ui
Pi  vx  dvx

and by vx 43  i  1 2 ρi  1
ui 1
Ei 1 $7 vx 8 0

1 vx v2x  2  Pi 

56 1

 ρi 

M ρi 

1

9

ρ exp 2π P  ρ

:


ui  1
Pi 1  vx  dvx


1

(7)

M ρ
u
P  vx  being the classical monodimensional maxwellian M ρ
u
P  vx  

(6)

vx  u  P ρ

2

;=<

(8)

An important property of this scheme is that it is possible to prove that it is positive and entropic under a classical CFL criteria (cf. [9]).

949

>@?ACBD?+EF? GHBJIK>

hji_kClnmpo

L Q PR

O M N

SUTWV*X

Y[Z

\[]_^*`

‡

a[b_cCdfe[g

‡

‡

}ƒ~CtFwH„…|-~ rDx{†|-z

qsrutFvDwyx{zW| x

}Uvu~|{€JvD‚x{zW|+x

Figure 1

KINETIC - FLUID COUPLING When ε ˆ 0 in (1), we can formally prove that f t
x
v  converges to the maxwellian distribution function (v ‰ vx
vy
vz [Š IR3 ) ‹

ρ
u
P  v Œ

ρ 2π P  ρ 

3

exp 2

 

vx  u 

2

v2y  v2z

P ρ Ž

(9)

where  ‰ ρ
ρ u
ρ E  is solution of (2), except in some areas as in the Knudsen layer. We now use the following notations (see the figure 1): Sq  1  2 is the uranium source, Xq  1 is the mesh having a frontier on the interface Sq 1 2 , Im  1 2 is the interface which separates the kinetic mesh Xm and the fluid mesh Xm  1 , and n in the normal entering in the fluid domain. The algorithm of the kinetic - fluid coupling is the following: - We initialize the kinetic and fluid domains by solving everywhere the transport equation ∂t f  v  ∇x f  0; - First stage: knowing the boundary condition Γ f luid  kinetic on the interface Im  1 2 , we solve the Boltzmann equation (1) in the kinetic area; - Second stage: knowing the boundary condition Γkinetic  f luid on the interface Im  1 2 , we solve the Euler system (2) with the previous kinetic scheme in the fluid area; - If the gloval level of convergence is not enough, we come back to the first stage. We now define the boundary conditions Γ f luid 

kinetic

and Γkinetic 

f luid .

Boundary conditions at the kinetic - fluid interface We define these conditions on the interface Im  nx  0 (see the figure 1).

950

1 2

where the normal n is such that

Boundary condition Γ f luid  kinetic : We now suppose that ρ
u
P  x ‘ fluid domain is known. Since the gas is supposed to be at the thermodynamic equilibrium in the fluid mesh Xm  1 , we impose that

where

‹

f x Š Im 

‹


v $‰

1 2

x  Xm 

1


v

if

is the maxwellian defined by (9) with ρ
u
P &

vx

. 0

ρ
u
P  x’

(10) Xm

“ 1.

Boundary condition Γkinetic  f luid : We now suppose that f x Š kinetic domain
v  is known. We use the kinetic decomposition of the macroscopic flux  mn  1  2 used in the numerical scheme (3) and deduced from the previous kinetic scheme built with the particular choice χ vx ” 0 12π exp  v2x  2  . Then, this condition is defined with  mn  1 2 % m  1 2 • & n where

 

m 1 2

vx 34  m  1 2 %–7 vx 2 0 and where 

1 vx v2  2

56

f x  Xm
v  dv

(11)

& is defined by the formula (7).  

n m 1 2

It is easy to prove that the boundary conditions (10) and (11) make conservative the kinetic - fluid coupling algorithm.

Asymptotic matching: the Marshak condition The boundary condition at the uranium source Sq f x Š Sq 


v @‰ Φ v 

1 2

 is defined by

1 2

if vx

 0

(12)

where Φ v  is a physical data coming from the modelization of the interaction of the electronic beam with the uranium source. This boundary condition induces that the gas is not at the thermodynamic equilibrium between the evaporation source Sq  1  2 and the fluid domain: this defines the Knudsen layer (cf. figure 1). Then, to asymptoticaly match the fluid domain on the source Sq 1 2 , we have to define a good boundary condition for the macroscopic flux  qn 1 2 . A priori, this can be done by studying the half space problem and by supposing that the exit of the Knudsen layer is sonic, and, then, by using the results of [12], [13] and [14]: this boundary condition defines a (Dirichlet) sonic condition. But, in AVLIS applications, there are situations where the exit of the Knudsen layer is not sonic: these situations appear when the temperature of the source is not uniform. Nevertheless, an other reason is that the results of [12], [13] and [14] can only be applied for the classical mono-species Boltzmann equation (1) and we would like to extend the results to the semi-classical multispecies Boltzmann equations

951

Ê

—Œ˜š™_›œŸžF ¢¡¢ž

ÓUÔÖÕØ×

£¥¤š¦¨§© ¤ª«H¬*§ «-­®°¯²±´³µ®·¶¹¸»º½¼°¾®·¿®°¯³n±ÁÀ´ºnÂj¶Á±

ÅÃÄÅÆ

ÅÃÇÅÆ

È É-ÉjÉ Ù{ښÛJÛ²ÜÞÝÏßWà…á…âÏÚDãHäšÝÏàWáWå æpå Ý´à

­ Ë

ÌWÍÏ΅ÐJÑJÒ

Figure 2: Physical domain and mesh

(cf. [10]). An other idea is to apply a boundary condition deduced from the Marshak condition proposed for the Navier-Stokes system (see [15]: in the case of the NavierStokes system, the Marshak condition defines a Robin condition). Let us note that the Marshak condition was firstly used for radiative transfer and neutron transport problems around the year 1940 in Los Alamos. It is possible to show (cf. [10]) that this Marshak condition when it is coupled with the previous kinetic scheme is simply defined by the macroscopic flux 

n q  1  2 %

 q  1  2 1

& n q 1 2

where  q  1 2 is defined by (11) by replacing f with Φ and where  the formula (7).

( n q 1 2

is defined by

NUMERICAL RESULTS FOR A SUBSONIC KNUDSEN LAYER The physical domain, the mesh and the kinetic - fluid interface are defined on the figure 2 whose geometry is bidimensional and axisymmetrical. The evaporation condition defined by Φ v  on the source S is a centered maxwellian being at the non uniform temperature defined on the figure 2 (3400, 3200 and 3000 Kelvin: the first temperature is due to the impact of the electronic beam). The density of the maxwellian is given by the saturation density of the uranium. The figures 3 and 4 show the uranium density with and without kinetic - fluid coupling, and with the sonic condition (figure 3) or with the Marshak condition (figure 4) for the asymptotic matching of the fluid domain on the source S. We can see that the kinetic - fluid coupling algorithm is correct only with the Marshak condition: indeed, the exit of the Knudsen layer is not sonic everywhere (and not monodimensional) because of the discontinuity 3400 / 3200 Kelvin at R  R1  8 5 < 10  3 m (cf. figure 5). Nevertheless, near the X axis, the sonic condition gives good results since the exit of the Knudsen layer is almost sonic at this place (cf. figure 5).

952


  
 
   



pì ð°ñpìpì ê é ípð°ñpìpì ê é ìpð°ñpìpì ïpð°ñpìuê î°ð°ñpìuê ò_ó ô¢õ

÷²ø½ùÏú û ü ý²ó ô=þ ÿ õ

@BA C D"EFCHGA I"JC A KMLHN O FA P2KEFQ"O A I R


  
 
   



@BA CSD2GA I JC A KTLHN O FA P2K EUF"Q"O A I"R +, )  "!$# % &' (*

:$>8?8:$: 7$6 ;$> ?$:$: 7$6 :$> ?$:$: =$>8?8:7 <8>8?8:7

è

5

ç è·é èpèpö ë ç è·é èuê ó ô»õ

è·é èpèpì è·é èpèpí è ·è é èuê ·è é èuê î è·é èuê ï è·é èuê ë è·é èuê è·é èpèpë

3 . 021

56 5$5$:

Figure 3: Sonic condition for the asymptotic matching

56 5$5$;

56 57

56 57 <

56 57 =

56 57 9

56 5$5$9

56 57

5

4 56 585$9

Figure 4: Marshak condition for the asymptotic matching

]_^a`bdcTeTfdgihj \ WY \ WX \ VW [ VW Z VW Y VW X V

VW VV X

VW VVY

VW VVZ

VW VV[

VW V\

VW V\8X

VW V\$Y

VW V\ Z

V

kml fdn

Figure 5: Mach number at the exit of the Knudsen layer (x o



9 10 

3

m)

CONCLUSION To diminish the CPU time used to simulate the uranium gas expansion in the AVLIS process, we have proposed to couple the resolution of the Boltzmann equation with the resolution of the Euler system by using the technics initialy used in [6] for classical aerodynamics problems. Moreover, we have shown that the Marshak condition coupled with a kinetic scheme is well adapted to asymptoticaly match the Euler area with the uranium source in order to take into account the effect of the Knudsen layer in the resolution of the Euler system. Indeed, it gives a boundary condition without supposing the value of the Mach number at the exit of the Knudsen layer and without solving any non linear system. Numerical results show that the results obtained with and without kinetic - fluid coupling are quasi similar and that the Marshak condition takes into account without difficulties

953

4 56 57

-/. 021

a Knudsen layer having a subsonic exit which is not the case with a sonic boundary condition deduced from the study of the half space problem (cf. [12], [13] and [14]). Finally, let us note that the technics proposed in that paper are extended in [10] to the semi-classical case and to the multispecies case where we take into account quantified energy transfers between the electronic metastable energy levels of atoms of different species, and that these technics can be applied to other vaporization problems as, for example, for the description of the gas expansion in the tail of a comet produced by sun radiations (cf. [16]).

REFERENCES 1. 2. 3. 4.

5.

6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.

Clerc M. and Plurien P. - Advanced uranium enrichment processes - Science, Research and Development, edited by the Commission of the European Communities, report EUR 10743 EN, 1986. Anderson J.B., Foch J.D., Shaw M.J., Stern R.C. and Wu B.J. - Statistical theory of electronic energy relaxation - In Rarefied Gas Dynamics, edited by V. Boffi and C. Cercignani, volume 1, Teubner, Stuttgart, p. 413-421, 1986. Roblin P., Rosengard A. and Nguyen T.T. - Model of electronic energy relaxation in the test-particle Monte-Carlo method - In Rarefied Gas Dynamics, edited by J. Harvey and G. Lord, volume 2, Oxford University Press, p. 522-528, 1995. Roblin P., Gonella C. and Chatain S. - Experimental measurements in gadolinium and copper atomic gas mixture flow during electron beam evaporation. Comparison with test-particle Monte-Carlo simulations - In Rarefied Gas Dynamics, edited by C. Shen, Bejing University Press, p. 940-945, 1997. Nishimura A. - Application of direct simulation Monte-Carlo method for analysis of AVLIS evaporation process - Proceedings of the 6th International Symposium on Advanced Nuclear Energy Research, Mito (Japan), edited by Japan Atomic Energy Research Institute, JAERI-CONF 95-005, p. 375-383, 1995. Bourgat J.F., Le Tallec P., Perthame B. and Qiu Y. - Coupling Boltzmann and Euler equations without overlapping - Sixth conference IUTAM on domain decomposition methods for partial differential equations, AMS Providence, Como, Italy, june 1992. Klar A., Neunzert H. and Struckmeir J. - Particle methods and domain decomposition - In Rarefied Gas Dynamics, edited by C. Shen, Bejing University Press, p. 263-272, 1997. Perthame B. - Second order Boltzmann schemes for compressible Euler equations in one and two space dimensions - SIAM J. Numer. Anal., 29, no. 1, p. 1-19, 1992. Villedieu P. and Mazet P.A. - Schémas cinétiques pour les équations d’Euler hors équilibre thermochimique - La Recherche Aérospatiale, edition Gauthier-Villars, 2, p. 85-102, 1995. Dellacherie S. - Coupling of the Wang Chang-Uhlenbeck equations with the multispecies Euler system - Submitted to Journal of Computational Physics, 2002. Dellacherie S. - About kinetic schemes built in axisymmetrical and spherical geometry - In Godunov Methods: theory and Applications, edited by E.F. Toro, Kluwer Academic / Plenum Publishers, p. 225-232, 2001. Sone Y. - Kinetic Theoretical Studies of the Half-Space - Problem of Evaporation and Condensation - Transport Theory and Statistical Physics, 29, p. 227-260, 2000. Ytrehus T. - Theory and experiments on gas kinetics in evaporation - In Rarefied Gas Dynamics, edited by J.L. Potter and M. Summerfield, AIAA, p. 1197-1212, 1977. Knight C.J. - Theoretical Modeling of Rapid Surface Vaporization with Back Pressure - AIAA Journal, 17, no. 5, p. 519-523, 1979. Golse F. - Applications of the Boltzmann equation within the context of upper atmosphere vehicle aerodynamics - Computer Methods in Applied Mechanics and Engineering, 75, p. 299-316, 1989. Pyarnpuu A.A., Shematovitch V.I., Svirschevsky S.B. and Titov E.V. - Nonequilibrium jet flows in the coma of a comet - In Rarefied Gas Dynamics, edited by C. Shen, Bejing University Press, 1997.

954