Numerical Methods of the Maxwell-Stefan Diffusion Equations and Applications in Plasma and Particle Transport

0.77 MB·English

by  Juergen Geiser

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
 