Numerical Methods of the Maxwell-Stefan Diffusion Equations and Applications in Plasma and Particle Transport 5 1 Ju¨rgen Geiser 0 2 RuhrUniversity of Bochum, n TheInstituteof Theoretical Electrical Engineering, a Universit¨atsstra”se 150, D-44801 Bochum, Germany J [email protected] 3 2 ] Abstract. In this paper, we present a model based on a local ther- A modynamic equilibrium, weakly ionized plasma-mixture model used for N medical and technical applications in etching processes. We consider a . simplified model based on the Maxwell-Stefan model, which describe h multicomponent diffusivefluxesin thegas mixture.Based on additional t a conditions tothefluxes,weobtain an irreducibleand quasi-positive dif- m fusion matrix. Such problems results into nonlinear diffusion equations, [ which are more delicate to solve as standard diffusion equations with Fickian’sapproach.Weproposeexplicittime-discretisationmethodsem- 1 beddedtoiterativesolversforthenonlinearities. Suchacombinational- v lows tosolvethedelicate nonlinear differentialequations moreeffective. 2 We present some first ternary component gaseous mixtures and discuss 9 the numerical methods. 7 5 0 . Keywords:Maxwell-Stefanapproach,Plasmamodel,Multi-componentmix- 1 ture, explicit discretization schemes, iterative schemes. 0 5 1 AMS subject classifications. 35K25,35K20, 74S10, 70G65. : v i X 1 Introduction r a We are motivated to understand the gaseousmixtures of a normal pressure and roomtemperatureplasma.Theunderstandingofnormalpressure,roomtemper- ature plasma applications is important for applications in medical and techni- cal processes. Since many years, the increasing importance of plasma chemistry basedonthe multi-component plasmais a key factorto understandthe gaseous mixture processes, see for low pressure plasma [7] and for atmospheric pressure regimes [8]. We consider a simplified Maxwell-Stefan diffusion equation to model the gaseousmixture of multicomponent Plasma.While most classicaldescription of the diffusion goes back to the Fickian’s approach, see [5], we apply the modern 2 description of the multicomponent diffusion based on the Maxwell-Stefan’s ap- proach, see [6]. The novel approach considers are more detail description of the flux and concentration, which are in deed not only proportional coupled as in the simplified Fickian’s approach. Here, we deal with a inter-species force bal- ance,whichallowstomodelcross-effects,e.g.,so-calledreversediffusion(up-hill diffusion in direction of the gradients). Suchamoredetailedmodelingresultsinirreducibleandquasi-positivediffusion- matrices,whichcanbereducedbytransformingwithreductionsortransforming withPerron-Frobeniustheoremstosolvablepartialdifferentialequations,see[1]. The obtained system of nonlinear partial differential equations are delicate to solve and numerically, we have taken into account linearisation methods, e.g., iterative fix-point schemes. The paper is outlined as follows. In section 2 we present our mathematical model. A possible reduced model for the further approximations is derived in Section 3. In section 4, we discuss the underlying numerical schemes. The first numerical results are presented in Section5.Inthecontents,thataregiveninSection6,wesummarizeourresults. 2 Mathematical Model Forthefullplasmamodel,weassumesthattheneutralparticlescanbedescribed as fluid dynamical model, where the elastic collision define the dynamics and few inelastic collisions are, among other reasons, responsible for the chemical reactions. To describe the individual mass densities, as well as the global momentum and the global energy as the dynamical conservation quantities of the system, corresponding conservation equations are derived from Boltzmann equations. The individual character of each species is considered by mass-conservation equations and the so-called difference equations. The extension of the non-mixtured multicomponent transport model, [7] is done with respectto the collisionintegralsrelatedto the right-hindside sources of the conservation laws. The conservation laws of the neutral elements are given as ∂ ∂ ρ + ·ρ u =m Q(s), ∂t s ∂r s s s n ∂ ∂ ρu+ · P∗+ρuu =−Q(e), ∂t ∂r m (cid:0) (cid:1) ∂ ∂ E∗ + · E∗ u+q∗+P∗·u =−Q(e). ∂t tot ∂r tot E (cid:0) (cid:1) where ρ : density of species i, ρ = N ρ , u : velocity, E∗ : total energy of s i=1 i tot the neutral particles. P 3 Furtherthe variableQ(s) is thecollisiontermofthemassconservationequa- n tion,Q(e) isthecollisiontermofthemomentumconservationequationandQ(e) m E is the collision term of the energy conservation equation. We derive the collision term with respect to the Chapmen-Enskog method, see [3], and achieve for the first derivatives the following results: m Q(s) =−∇·(ρ Vj), (1) s n i i Xj=0 ns Q(e) =− ρ F (2) m i i Xi=1 ns Q(e) =− ρ ρF (u+ V(j)), (3) E i i i Xi=1 Xj=0 where i=1,...,n , F is an external force per unit mass (see Boltzmann equa- s i tion), further the diffusion velocity is given as: V0 =0 (4) i N ∆T V1 =− D (d +k ), (5) i ij j Tj T Xj=1 where N d =0, i=1 i P ∇p ρ i d =∇x +x − F , (6) i i i i p ρ d =d −y d∗, (7) i i i j Xj where xi = nns is the molar fraction of species i. We haveanadditionalconstraintbasedonthe mass fractionof eachspecies: ∂ y +∇y =R (y ,...,y ), (8) i i i 1 N ∂t wherey isthemassfractionofspeciesi,R isthenetproductionrateofspecies i i i due to his reactions. Remark 1. Thefullmodelproblemconsiderafullcoupledsystemofconservation laws and Maxwell-Stefan equations. Each equations are coupled such that the gaseousmixtureinfluencesthe transportequationsandviceverse.Inthefollow- ing, we decouple the equations system and consider only the delicate Maxwell- Stefan equations. 3 Simplified Model with Maxwell-Stefan Diffusion Equations Wediscussinthefollowingamulticomponentgaseousmixturewiththreespecies (ternarymixture).Themodel-problemisdiscussedintheexperimentsofDuncan and Toor, see [4]. 4 Here, they studied an ideal gaseous mixture of the following components: – Hydrogen (H , first species), 2 – Nitrogen (N , second species), 2 – Carbon dioxide (CO , third species). 2 TheMaxwell-Stefanequationsaregivenforthethreespeciesas(seealso[2]): ∂ ξ +∇·N =0, 1≤i≤3, (9) t i i 3 N =0, (10) j Xj=1 ξ N −ξ N ξ N −ξ N 2 1 1 2 3 1 1 3 + =−∇ξ , (11) 1 D D 12 13 ξ N −ξ N ξ N −ξ N 1 2 2 1 3 2 2 3 + =−∇ξ , (12) 2 D D 12 23 where the domain is given as Ω ∈IRd,d∈IN+ with ξ ∈C2. i For such ternary mixture, we can rewrite the three differential equations (9) and (11)-(12) with the help of the zero-condition (10) into two differential equations, given as: ∂ ξ +∇·N =0, 1≤i≤2, (13) t i i 1 N +αN ξ −αN ξ =−∇ξ , (14) 1 1 2 2 1 1 D 13 1 N −βN ξ +βN ξ =−∇ξ , (15) 2 1 2 2 1 2 D 23 where α= 1 − 1 , β = 1 − 1 . (cid:16)D12 D13(cid:17) (cid:16)D12 D23(cid:17) Further we have the relations: – Third mole-fraction: ξ =1−ξ −ξ , 3 1 2 – Third molar flux: N =−N −N . 3 1 2 4 Numerical Methods Inthe following,wediscussthenumericalmethods,whicharebasedoniterative schemes with embedded explicit discretization schemes. We apply the following methods: – Iterative Scheme in time (Global Linearisation, Matrix Method), – Iterative Scheme in Time (Local Linearisation with Richardson’s Method). Forthespatialdiscretization,weapplyfinitevolumeorfinitedifferencemeth- ods. The underlying time-discretization is based on a first order explicit Euler method. 5 4.1 Iterative Scheme in time (Global Linearisation, Matrix Method) We solve the iterative scheme: ξn+1 =ξn−∆tD Nn, (16) 1 1 + 1 ξn+1 =ξn−∆tD Nn, (17) 2 2 + 2 A B Nn+1 −D ξn+1 1 = − 1 (18) (cid:18)C D(cid:19)(cid:18)N2n+1(cid:19) (cid:18)−D−ξ2n+1(cid:19) for j = 0,...,J , where ξn = (ξn ,...,ξn )T, ξn = (ξn ,...,ξn )T and 1 1,0 1,J 2 2,0 2,J I ∈ IRJ+1 × IRJ+1, Nn = (Nn ,...,Nn )T, Nn = (Nn ,...,Nn )T and J 1 1,0 1,J 2 2,0 2,J I ∈ IRJ+1 ×IRJ+1, where n = 0,1,2,...,N and N are the number of J end end time-steps, i.d. N =T/∆t. end The matrices are given as: A,B,C,D ∈IRJ+1×IRJ+1, (19) 1 A = +αξ , j =0...,J, (20) j,j 2,j D 13 B =−αξ , j =0...,J, (21) j,j 1,j C =−βξ , j =0...,J, (22) j,j 2,j 1 D = +βξ , j =0...,J, (23) j,j 1,j D 23 A =B =C =D =0, i,j =0...,J, i6=J, (24) i,j i,j i,j i,j means the diagonal entries given as for the scale case in equation (95) and the outer-diagonalentries are zero. The explicit form with the time-discretization is given as: Algorithm 1 1.) Initialisation n=0: N0 A˜ B˜ −D ξ0 1 = − 1 (25) (cid:18)N20(cid:19) (cid:18)C˜ D˜(cid:19)(cid:18)−D−ξ20(cid:19) where ξ0 = (ξ0 ,...,ξ0 )T, ξ0 = (ξ0 ,...,ξ0 )T and ξ0 = ξin(j∆x), ξ0 = 1 1,0 1,J 2 2,0 2,J 1,j 1 2,j ξin(j∆x), j =0,...,J and given as for the different intialisations, we have: 2 1. Uphill example 0.8 if0≤x<0.25, ξin(x)=1.6(0.75−x)if0.25≤x<0.75, , (26) 1 0.0 if0.75≤x≤1.0, ξin(x)=0.2, for allx∈Ω =[0,1], (27) 2 2. Diffusion example (Asymptotic behavior) 0.8if0≤x∈0.5, ξin(x)= , (28) 1 (cid:26)0.0else, ξin(x)=0.2, for allx∈Ω =[0,1], (29) 2 6 The inverse matrices are given as: A˜,B˜,C˜,D˜ ∈IRJ+1×IRJ+1, (30) 1 A˜ =γ ( +βξ0 ), j =0...,J, (31) j,j j D 1,j 23 B =γ αξ0 , j =0...,J, (32) j,j j 1,j C =γ βξ0 , j =0...,J, (33) j,j j 2,j 1 D =γ ( +αξ0 ), j =0...,J, (34) j,j j D 2,j 13 D D 13 23 γ = , j =0...,J, (35) j 1+αD ξ0 +βD ξ0 13 2,j 23 1,j A˜ =B˜ =C˜ =D˜ =0, i,j =0...,J, i6=J, (36) i,j i,j i,j i,j Further the values of the first and the last grid points of N are zero, means N0 =N0 =N0 =N0 =0 (boundary condition). 1,0 1,J 2,0 2,J 2.) Next time-steps (till n=N ): end 2.1) Computation of ξn+1 and ξn+1 1 2 ξn+1 =ξn−∆tD Nn, (37) 1 1 + 1 ξn+1 =ξn−∆tD Nn, (38) 2 2 + 2 2.2) Computation of Nn+1 and Nn+1 1 2 Nn+1 A˜ B˜ −D ξn+1 1 = − 1 (39) (cid:18)Nn+1(cid:19) (cid:18)C˜ D˜(cid:19)(cid:18)−D ξn+1(cid:19) 2 − 2 where ξn =(ξn ,...,ξn )T, ξn =(ξn ,...,ξn )T and the inverse matrices are 1 1,0 1,J 2 2,0 2,J given as: A˜,B˜,C˜,D˜ ∈IRJ+1×IRJ+1, (40) 1 A˜ =γ ( +βξn+1), j =0...,J, (41) j,j j D 1,j 23 B =γ αξn+1, j =0...,J, (42) j,j j 1,j C =γ βξn+1, j =0...,J, (43) j,j j 2,j 1 D =γ ( +αξn+1), j =0...,J, (44) j,j j D 2,j 13 D D 13 23 γ = , j =0...,J, (45) j 1+αD ξn+1+βD ξn+1 13 2,j 23 1,j A˜ =B˜ =C˜ =D˜ =0, i,j =0...,J, i6=J. (46) i,j i,j i,j i,j Further the values of the first and the last grid points of N are zero, means Nn =Nn =Nn =Nn =0 (boundary condition). 1,0 1,J 2,0 2,J 3.) Do n=n+1 and goto 2.) 7 4.2 Iterative Scheme in Time (Local Linearisation with Richardson’s Method We solve the iterative scheme given in the Richardson iterative scheme: ξn+1,k =ξn−∆tD Nn+1, (47) 1 1 + 1 ξn+1,k =ξn−∆tD Nn+1, (48) 2 2 + 2 An+1,k−1 Bn+1,k−1 Nn+1 −D ξn+1,k−1 1 = − 1 (49) (cid:18)Cn+1,k−1 Dn+1,k−1(cid:19)(cid:18)Nn+1(cid:19) (cid:18)−D ξn+1,k−1(cid:19) 2 − 2 for j = 0,...,J , where ξn = (ξn ,...,ξn )T, ξn = (ξn ,...,ξn )T and 1 1,0 1,J 2 2,0 2,J I ∈ IRJ+1 × IRJ+1, Nn = (Nn ,...,Nn )T, Nn = (Nn ,...,Nn )T and J 1 1,0 1,J 2 2,0 2,J I ∈ IRJ+1 ×IRJ+1, where n = 0,1,2,...,N and N are the number of J end end time-steps, i.d. N =T/∆t. end Further k =1,2,...,K is the iteration index with where ξn+1,0 = (ξn ,...,ξn )T, ξn+1,0 = (ξn ,...,ξn )T and I ∈ IRJ+1 × 1 1,0 1,J 2 2,0 2,J J IRJ+1 is the start solution given with the solution at t=tn. The matrices are given as: An+1,k−1,Bn+1,k−1,Cn+1,k−1,Dn+1,k−1 ∈IRJ+1×IRJ+1, (50) 1 An+1,k−1 = +αξn+1,k−1, j =0...,J, (51) j,j D 2,j 13 Bn+1,k−1 =−αξn+1,k−1, j =0...,J, (52) j,j 1,j Cn+1,k−1 =−βξn+1,k−1, j =0...,J, (53) j,j 2,j 1 Dn+1,k−1 = +βξn+1,k−1, j =0...,J, (54) j,j D 1,j 23 An+1,i−1 =Bn+1,i−1 =Cn+1,i−1 =Dn+1,i−1 =0, i,j =0...,J, i6=J,(55) i,j i,j i,j i,j means the diagonal entries given as for the scale case in equation (95) and the outer-diagonalentries are zero. The explicit form with the time-discretization is given as: Algorithm 2 1.) Initialisation n=0with an explicit time-step(CFL condition is given): N0 A˜ B˜ −D ξ0 1 = − 1 (56) (cid:18)N20(cid:19) (cid:18)C˜ D˜(cid:19)(cid:18)−D−ξ20(cid:19) where ξ0 = (ξ0 ,...,ξ0 )T, ξ0 = (ξ0 ,...,ξ0 )T and ξ0 = ξin(j∆x), ξ0 = 1 1,0 1,J 2 2,0 2,J 1,j 1 2,j ξin(j∆x), j =0,...,J and given as for the different intialisations, we have: 2 1. Uphill example 0.8 if0≤x<0.25, ξin(x)=1.6(0.75−x)if0.25≤x<0.75, , (57) 1 0.0 if0.75≤x≤1.0, ξin(x)=0.2, for allx∈Ω =[0,1], (58) 2 8 2. Diffusion example (Asymptotic behavior) 0.8if0≤x∈0.5, ξin(x)= , (59) 1 (cid:26)0.0else, ξin(x)=0.2, for allx∈Ω =[0,1], (60) 2 The inverse matrices are given as: A˜,B˜,C˜,D˜ ∈IRJ+1×IRJ+1, (61) 1 A˜ =γ ( +βξ0 ), j =0...,J, (62) j,j j D 1,j 23 B =γ αξ0 , j =0...,J, (63) j,j j 1,j C =γ βξ0 , j =0...,J, (64) j,j j 2,j 1 D =γ ( +αξ0 ), j =0...,J, (65) j,j j D 2,j 13 D D 13 23 γ = , j =0...,J, (66) j 1+αD ξ0 +βD ξ0 13 2,j 23 1,j A˜ =B˜ =C˜ =D˜ =0, i,j =0...,J, i6=J, (67) i,j i,j i,j i,j Further the values of the first and the last grid points of N are zero, means N0 =N0 =N0 =N0 =0 (boundary condition). 1,0 1,J 2,0 2,J 2.) Next timesteps (till n = N ) (iterative scheme restricted via the CFL end condition based on the previous iterative solutions in the matrices): 2.1) Computation of ξn+1,I and ξn+1,I 1 2 ξn+1,k =ξn−∆tD Nn+1, (68) 1 1 + 1 ξn+1,k =ξn−∆tD Nn+1, (69) 2 2 + 2 2.2) Computation of Nn+1,k−1 and Nn+1,k−1 1 2 Nn+1 A˜n+1,k−1 B˜n+1,k−1 −D ξn+1,k−1 1 = − 1 (70) (cid:18)Nn+1(cid:19) (cid:18)C˜n+1,k−1 D˜n+1,k−1(cid:19)(cid:18)−D ξn+1,k−1(cid:19) 2 − 2 where ξn =(ξn ,...,ξn )T, ξn =(ξn ,...,ξn )T and the inverse matrices are 1 1,0 1,J 2 2,0 2,J given as: A˜n+1,k−1,B˜n+1,k−1,C˜n+1,k−1,D˜n+1,k−1 ∈IRJ+1×IRJ+1, (71) 1 A˜n+1,k−1 =γ ( +βξn+1,k−1), j =0...,J, (72) j,j j D 1,j 23 Bn+1,k−1 =γ αξn+1,k−1, j =0...,J, (73) j,j j 1,j Cn+1,k−1 =γ βξn+1,k−1, j =0...,J, (74) j,j j 2,j 1 Dn+1,k−1 =γ ( +αξn+1,k−1), j =0...,J, (75) j,j j D 2,j 13 9 D D 13 23 γ = , j =0...,J, (76) j 1+αD ξn+1,k−1+βD ξn+1,k−1 13 2,j 23 1,j A˜n+1,k−1 =B˜n+1,k−1 =C˜n+1,k−1 =D˜n+1,k−1 =0, i,j =0...,J, i6=J(.77) i,j i,j i,j i,j Further the values of the first and the last grid points of N are zero, means Nn+1 =Nn+1 =Nn+1 =Nn+1 =0 (boundary condition). 1,0 1,J 2,0 2,J Further k =1,2,...,K is the iteration index with where ξn+1,0 = (ξn ,...,ξn )T, ξn+1,0 = (ξn ,...,ξn )T and I ∈ IRJ+1 × 1 1,0 1,J 2 2,0 2,J J IRJ+1 is the start solution given with the solution at t=tn. 3.) Do n=n+1 and goto 2.) 5 Numerical Experiments Inthefollowing,weconcentrateonthefollowingthreecomponentsystem,which is given as: ∂ ξ +∂ N =0, 1≤i≤3, (78) t i x i 3 N =0, (79) j Xj=1 ξ N −ξ N ξ N −ξ N 2 1 1 2 3 1 1 3 + =−∂ ξ , (80) x 1 D D 12 13 ξ N −ξ N ξ N −ξ N 1 2 2 1 3 2 2 3 + =−∂ ξ , (81) x 2 D D 12 23 where the domain is given as Ω ∈IRd,d∈IN+ with ξ ∈C2. i The parameters and the initial and boundary conditions are given as: – D = D = 0.833 (means α = 0) and D = 0.168 (Uphill diffusion, 12 13 23 semi-degenerated Duncan and Toor experiment), – D = 0.0833,D = 0.680 and D = 0.168 (asymptotic behavior, Duncan 12 13 23 and Toor experiment, see [4]), – J =140 (spatial grid points), – The time-step-restriction for the explicit method is given as: ∆t≤ (∆x)2 , 2max{D12,D13,D23} – The spatial domain is Ω =[0,1], the time-domain [0,T]=[0,1], – The initial conditions are: 1. Uphill example 0.8 if0≤x<0.25, ξin(x)=1.6(0.75−x)if0.25≤x<0.75, , (82) 1 0.0 if0.75≤x≤1.0, ξin(x)=0.2, for allx∈Ω =[0,1], (83) 2 10 2. Diffusion example (Asymptotic behavior) 0.8if0≤x∈0.5, ξin(x)= , (84) 1 (cid:26)0.0else, ξin(x)=0.2, for allx∈Ω =[0,1]. (85) 2 – The boundary conditions are of no-flux type: N =N =N =0,on∂Ω×[0,1]. (86) 1 2 3 We could reduce to a simpler model problem as: ∂ ξ +∂ ·N =0, 1≤i≤2, (87) t i x i 1 N +αN ξ −αN ξ =−∂ ξ , (88) 1 1 2 2 1 x 1 D 13 1 N −βN ξ +βN ξ =−∂ ξ , (89) 2 1 2 2 1 x 2 D 23 where α= 1 − 1 , β = 1 − 1 . (cid:16)D12 D13(cid:17) (cid:16)D12 D23(cid:17) We rewrite into: ∂ ξ +∂ ·N =0, (90) t 1 x 1 ∂ ξ +∂ ·N =0, (91) t 2 x 2 D113 +αξ2 −αξ1 N1 = −∂xξ1 (92) (cid:18) −βξ2 D123 +βξ1(cid:19)(cid:18)N2(cid:19) (cid:18)−∂xξ2(cid:19) and we have ∂ ξ +∂ ·N =0, (93) t 1 x 1 ∂ ξ +∂ ·N =0, (94) t 2 x 2 N1 = D13D23 D123 +βξ1 αξ1 −∂xξ1 . (95) (cid:18)N2(cid:19) 1+αD13ξ2+βD23ξ1 (cid:18) βξ2 D113 +αξ2(cid:19)(cid:18)−∂xξ2(cid:19) The next step is to apply the semi-discretization of the partial differential operator ∂ . ∂x We apply the first differential operator in equation (93) and (94) as an for- ward upwind scheme given as −1 0... 0  1−1 0 ... 0 ∂ =D+ = 1 · ... ... ... ... ... ∈ IR(J+1)×(J+1) (96) ∂x ∆x    0 1−1 0    0 ... 0 1−1  