American Journal of Computational Mathematics, 2011, 1, 163-175
doi:10.4236/ajcm.2011.13019 Published Online September 2011 (http://www.SciRP.org/journal/ajcm)
Copyright © 2011 SciRes. AJCM
Operator Splitting Method for Coupled Problems:
Transport and Maxwell Equations
Jürgen Geiser
Department of Mat hematics, Humboldt-Universität zu Berlin, Unter den Linden, Berlin, Germany
E-mail: geiser@mathematik.hu-berlin.de
Received March 31, 2011; revised April 13, 2011; accepted May 2, 2011
Abstract
In this article a new approach is considered for implementing operator splitting methods for transport prob-
lems, influenced by electric fields. Our motivation came to model PE-CVD (plasma-enhanced chemical va-
por deposition) processes, means the flow of species to a gas-phase, which are influenced by an electric field.
Such a field we can model by wave equations. The main contributions are to improve the standard discretiza-
tion schemes of each part of the coupling equation. So we discuss an improvement with implicit Runge-
Kutta methods instead of the Yee’s algorithm. Further we balance the solver method between the Maxwell
and Transport equation.
Keywords: Operator Splitting Method, Initial Value Problems, Iterative Solver Method, Stability Analysis,
Beam Propagation Methods, Transport and Maxwell Equations
1. Introduction
We motivate our study by simulating thin film deposition
processes that can be realized by PE-CVD (plasma en-
hanced chemical vapor deposition) processes, see [1,2].
For the deposition process, the influence of the electric
fields to the transported gases in a plasma reactor is very
important, see [3]. Therefore we deal with a simplified
model of a coupled transport and Maxwell equations.
While the transport equations modeled the transport of
gaseous species and the Maxwell equation the influence
of the underlying flow field.
We deal with the following equations


2
2
2
2
=,
,
txz y
uu
uvExyv D
xy
u
x
u
Dy

 

(1)


00
,,=, ,uxytuxy (2)
  
,=,,, 0,
xz
Hxy E
,
x
yt T
ty

 (3)
  
,=,,, 0,
yz
Hxy E
 
,1
=,
,,0, ,
y
zx
s
ource
H
Exy HJ
txy
xyt T






(5)
where u is the concentration of the gaseous species,
z
E
is the electric field and ,
x
y
H
H is the corresponding
magnetic field in two dimensions. Further
t
xy
v=,vv
is the influenced velocity of the transport equation.
We concentrate on the numerical modeling and simu-
lation of electrical fields, which are coupled with trans-
port equations.
Several methods exist to solve electric field and are of
interest.
One method for a stationary case of the electric field is
a propagation method (BPM). This is a powerful tool to
analyze linear and nonlinear light propagation in axially
varying waveguides like directional couplers, tapered
waveguides, S-shaped bent waveguides, and optical fi-
bers [4-7]. The method has its origin in the field of
propagation of electromagnetic beams in atmosphere,
where the multi-physics modeling was done on the as-
sumption that “the continuous gain medium may be ap-
proximated by a series of gain sheets with free propaga-
tion between the sheets” [8,9]. As it will be shown later
on, this method is in fact a Strang-Marchuk operator
splitting method [10,11]. Here we first describe the BPM
[12]. We introduce the iterative splitting idea to couple
,
yt T
tx

 (4)
164 J. GEISER
Maxwell and Transport equations. Further a splitting
analysis is presented. Numerical experiments are pre-
sented with respect to decoupled and coupled differential
equations.
The paper is organized as follows. The discretization
methods are described in Section 2. In Section 3, the
applied operator splitting methods are presented. The
error analysis of the coupled methods is studied in Sec-
tion 4. The experiments of the new discretization meth-
ods and splitting methods are performed in Section 5. At
the end of this paper we introduce future works.
2. Discretization Method of the Maxwell
Equation
In the following we discuss the discretization methods
for the Maxwell equation.
2.1. FDTD Method: Yee’s Scheme
Yee’s scheme is the standard finite difference time-do-
main (FDTD) discretization of the following time de-
pendent Maxwell curl equations
0=
H,
E
rt


(6)
0=
E,
H
rt


(7)
where is the electric field,

=,, ,,Exyz
EEE xyt
=
H


,, ,,
xyz
H
HH xyt is the magnetic field,
=,
r
x
y

=1
r
is the relative permittivity (given data),
(non-
magnetic material) is the magnetic permeability. Here
0
, 0
are constants. It can be shown that if the diver-
gence free conditions and

=0E
r

=0H
are satisfied at , then they are satisfied for all time.
This is the case for our setting. Therefore it is enough to
consider only the above curl equation. Rewriting them
component-wise, we get in our case
=0t
0=
y
xzE
HE
ty
z
(8)
0=
yx
z
HEE
tz

x
(9)
0=y
x
z
r
H
H
E
tx


y
(10)
Let
x
,
y
are spatial discretizations, and t
is a
time step. We use the following notation
 
,=, ,
n.
F
ij Fixjynt
 (11)
Let
represents a spatial coordinate such as
x
, .
The goal of Yee’s scheme is to compute the approxima-
tions for the various components
y
E
of
E
and
H
of
H
at the following spatial locations and temporal
instants:
:
:
=:
s
patial coordinateinte
other spatial coordinateg
time integer
half
te in
=:
nteger
ger
erE
(12)
:
:
spatial coordinatei
H
otherspatial coordinatalfer
timehalf integer
e hinteg
(13)
Thus the distributions/grid of various components are
staggered in space and in time. This is one of the two
unique characteristics of the Yee’s scheme. The second
unique characteristic is that the various spatial deriva-
tives in Equations (8) - (10) are computed across the one
spatial cell, i.e. the difference center for the central dif-
ference approximation of the spatial derivative is the mid
point of one cell length in the corresponding direction of
the derivative. Thus the Yee’s scheme approximates
Equations (8) - (10) at the following points:
Equation (8),12ix jnt,y
 (14)
Equation (9)12,,jixynt
 (15)

Equation (10),,12ixjynt
 (16)
Such a staggered uncollocated arrangement gives the
Yee’s scheme several nice numerical and physical prop-
erties, see [13]. Then we get finite-difference approxima-
tions as:

,
z
1
2
1
2
0
1
,2
11
=, ,1,
2
xz
Hij
t
n
x
nnn
H
ijE ijij
y



 
 
 
 E
(17)

,
z
1
2
1
2
0
1,
2
11
=,1 ,
2
yz
Hi j
t
n
y
n,
nn
H
ijEij ij
x



 
 
 
 E
(18)

1
11
2
,,
2
0
11
22
0
,
=,
11
1
22
11
,,
22
nn
yy
r
nn
xx
r
Eij
Eij
t
n
z
n
z
1
.
H
ijHij
x
tHij Hij
y







(19)
 

 
 
Copyright © 2011 SciRes. AJCM
J. GEISER
165
In Equation (19) the relative permittivity r
is com-
puted at the corresponding difference center as given by
Equation (16). At the interface between two media, r
is approximated by the average value.
Conditions for the Yees algorithm:
The CFL stability condition for the Yee’s FDTD
method is
 
2
11 1
tc2
x
y
 

(20)
where c is the speed of light in vacuum, see [13].
To restrict the unbounded domain to finite domain,
one uses absorbing boundary condition like the perfectly
matched layers, see [14,15].
Remark 1. Often for more accurate problems a Yees
algorithm which is second order in time and second or-
der in space is often to low. For higher order methods in
time and space can be constructed but are often to deli-
cate and expensive to implement, see [16,17]. We pro-
pose to improve with higher order implicit Runge-Kutta
methods with an idea to sparse matrices schemes, which
saves additional memory.
2.2. Improved Time Discretization Methods for
Maxwell Equation
Based on the problem of reconstructing a higher order
Yee’s algorithm, we deal with separate improvement of
the discretization schemes.
While the spatial discretization of the Yee’s algorithm
is a second order difference scheme, the time discretiza-
tion is also only a second order scheme.
Here we see the deficits of only improving the spatial
scheme with higher order schemes and leave the time-
discretization with a second order scheme.
We propose an improved time-discretization scheme
of higher order and apply fine spatial grids, while the time
error is at least larger, see [18].
We deal with higher order time-discretization methods.
Therefore we propose the Runge-Kutta as adapted time-
discretization methods to reach higher order results. For
the time-discretization we use the following higher order
discretization methods.
We deal with the following semi-discretized partial
differential equations, such equations are used in each
iterative splitting step:

=
u,
A
uft
t
(21)

=
n
n
uut, (22)
where
A
is the operator that we implicit solve in the
equation and
 
=
f
tBut
is the explicit operator, with a
previous solution , e.g. last iterative solution.
u
2.2.1. Higher Order Time-Discretization Methods
with Runge-Kutta Methods
We deal with the following Maxwell equation, given as:

12
1
=
1
=I,
H
tx
, =
yx
z
xyy x
H
EJ
y
H
HJBH BHJ







 (23)

1
1
==II=
xz1
z
z
HEEC
ty

E
 (24)

2
1
== III=
y1
z
z
z
HECE

E
tx
 (25)
For the boundary conditions we assume periodic
boundary conditions. That means we use the identifica-
tion

,= 1,
zz
ENi Ei
(26)

,=,1 =1,,
zz
EiN EiiN. (27)
Remark 2. For the stationary field, we apply a peri-
odic boundary condition, which is sufficient. The Mur
absorbing boundary condition, see [5], is used for the in-
stationary field, while respecting the influence of the
changes at the boundaries.
To get a first realization of an open boundary in the
case of the line-source we use symmetry and a combina-
tion of PBC and Mur’s first order ABC. For the bound-
arys orthogonal to the propagation direction of the field
(left-right) it is useful to work with Mur’s ABC.
2.2.2. Mur’s A B C
We can interpret the electromagnetical field as a wave
that has to fulfill the homogeneous wave equation.
2
00
222
11
=0=0==rr
ct c




 




(28)
22 2
2222
1=0
xyct

 


 


(29)
222
2
1=0
xy t
DD D
c

 


(30)
22
2
22
2
1
c


1
1
1=0
y
xt
t
y
xt
t
Dc
DD D
Dc
DD
cD






 


(31)
Copyright © 2011 SciRes. AJCM
166 J. GEISER
22
()
11
11
xt xt
DDVDDV
cc








  =0
(32)
=0
xx

□□ (33)
Waves that satisfy only propagate in
=0
x
x
-direction and those that satisfy only
propagate in
=0
x
x
-direction. An analogous formulation
can be given for the and direction.
yy
To handle
it is comfortable to do a Taylor ex-
pansion around 0.






2
3
222
34
5
22
1
1=1 00
121
0
1
V
VV
VV
VVOV
V
 

2
V
(34)

2
1
=12VOV4
(35)

2
=1OV (36)
Considering (36) equation (32) turns to
11
=0,
xct xct






(37)
which is Mur’s ABC with first order accuracy. As a first
attempt to model an open boundary we will use this.
Left boundary (0
=
x
x)
For the left boundary we have do discretize the fol-
lowing equation:
1
=.
x
ct
 

(38)
This can be done with a FDM-scheme as follows.
 
11
22
1
=02
1
1
2
=
1
=,2
133
=,,
22
nn
xx x
nn
n
tt
jj
xx
jj
tt



 




 
 


 

 

,1;
(39)
with
  

  

1
1
2
1
1
2
1
,2 =,2,2;
2
1
,1 =,1,1
2
nnn
nnn
jj
jjj


j
(40)
and
 

 

111
31
,=,2 ,1
22
31
,=,2,1
22
nnn
nnn
jj
jjj









;j
(41)
this leads to
  
,1
n

11
,1 =,2,2
nn n
ct x
j
jj
ct x

j



(42)
his tool does not satisfy completely
ly has first order accuracy and even more
important it only absorbs the part of the wave that
propagates orthogonal to the boundary.
But there are also a few advantages. Mur’s ABC h
It is easy to see that t
because it on
as to
be applied only to the
z
E field because
x
H
and
y
H
are dealt with automatically through the ordinary update-
step. The second advantage of Mur's ABC is the low
numeric expense.
For the boundaries parallel to the propagation direc-
tion (top and bottom) we use the PBC. Themmetry
our setting garanties that the inflow and the outflow of
the field equalize each other.
But with the ey
sy of
e on the next simulations with less
sy
is discretised and is calcu-
la
mmetry it seams to be necessary to use perfectly
matched layers.
These 3 equations above mark the starting point. The
spatial part of each equation
ted with the help of the matrix-operators 1212
,,,BBCC
(centered differences corresponding to the 2 dimensional
Yee-lattrice).
In the following we are using the general Butcher-ta-
ble for (3-stage) Runge-Kutta-methods to get a clear no-
tation.
(43)
denote the stepping time and
 
=
ii
j
j
tt tc
. t
Let
The ste from to p

i
t

1i
t
ng way.
in (23) - (
written i
25) can now be
n the followi
 
1II
112 233
=
ii
zz
EEtbkbkbk
 
I
(44)
II
1II II
112233
=
ii
xx
H
Htbkbk
 )
 
bk (45
III III III
2 23 3
1
11
=
ii
yy
H
Htb
 kbkbk 
 
(46)
where
iii
i
M=,,,
jj
xy
z
jjj
k MHHEJ
I, II, III M for
and
1, 2,3j. With
 
3I
=1
==

iii
z
jz jll
zjl
EEtE tak
 
(47)
 

II
=1
==
iii
x l
xjl
3
xj jl
H
Ht Ht

ak
)
(48
 


3II
=1
==
iii
yj yjll
yjl
I
H
HtHt ak
 

(49)
Copyright © 2011 SciRes. AJCM
J. GEISER
167
For a better legibility and because the focused point of
time does not change, we write
 

=which is known(in our ca
ii
jj
JJtse) (50)

,,,
j
jj
zxy
EHH
- (50) give 9 equati
j
J in-
stead of
Comb 7)
 
,,,
iii
i
j
zxy
jjj
EHHJ .
ining (23) - (25) and (4ons

12
=1
=
jlll
zzjy x
l
EEtal BHBHJ
  


(51)
31

3
1
=1
=1
j
l
x
xjz
l
H
HtalCE
 

(52
 )
3
2
=1
1
=1
jl
yyj z
l
HH talCE j

 


,2,3
Remark 3.
(53)
1
and 1
are
trices, such that also realized as ma-

1=1 ,
x
y

and
1=1,
x
y

.
Remark 4. The scheme above is only correct for iso-
tropic media because in the not isotropic case i
essary to consider t is nec-

=,
x
y

.
of
eq
ith equatst
Taking the 6uations(52) and (53) and putting
them together w the 3ions of (51) lead to he
following linear equation system which needs to be
solved.
1
QE

2
3
z
z
z
QE
QE
R




11 12 13
212223
3132 33
1
2
3
,,,
=,,,
,,,
AA AA AA
AA AA AA
AA AA AA
z
z
z
CSIRCSRC S
RCSRCSI RCS
RC SRC SRC SI
E
E
E










(54)
1
2
3
1222
11 1213
2222
21 2223
3222
31 3233
=
z
z
z
z
z
z
QEaSIaS aSE
QEaSaSIaSE
QEaSaSaSIE
 

 

 
 

 
(55)
2
11
22
3
00 00
=0000
00 003
z
z
z
z
z
z
QESI E
QESI E
QESIE
 
 

 

 
 


 

 
(56)

12
1
:= ,1,
jA A
jyxj
QtRBHBHRJ

(57)
A
(59)
:= -
A
j
Rjthrowof (58)
:= -
A
j
CjthcolumnofA
1:= 1,1,1 (60)
123
:=, ,
J
JJJ (61)

221
BC
(62)
2
1
11
:=St BC

22
:==,
AA
i j
RC (63)
ij ij
aA
(65)
22
11 1213
2222
21 2223
22 2
31 3233
:=


aa
aa a
aa a
2
a
(64)
2
2
2
0
=
0
ij
ij
ij
nmtimes
a
nmtimes
a





a

:=
I
Identity nm
(66)
wher and e ,
ij
a i
b
j
c
sel
are the Runge-Kut
Topreciy: If we have
region, there are
ta coefficients.
points in our be more nm
3nm
equations te.
With this result we are able to calculate (52) and (53).
So that it is finally possible to
Remark 5. For an optimiz
, n
of DIRK
fastert reducing
th n
sc
o solv
do the step ((44) - (46)).
ation of the time-discretiza-
tion schemewe caneglect some of the outerdiagonals
the RK methods, which leads to S methods. We
have the b enefit in computations , withou
e accuracy. For higher time-discretizations we have to
taken into account also higher spatial discretizatio
heme.
2.2.3. Stability Analysis of the Implicit Discretizations
We deal with the following discretized equation systems:

 
1
1 1
,121
=
nnnnnn
ziz y xi
EIAECHCHJuJu
 


(67)
where i is the iteration index of the coupling scheme.
Definition 1. We have a positive definite matrix M
( nn
real symmetric matrix), if for all
ors z with real entries where
tes the transpose of z.
T>0zAz

n
z,non-zero vect
T
z deno
Example 1. For finite difference discretization, e.g.
1121x
, it is sufficient to show, that the sum of
e outer-diagonals are equ al or less than the diagonal.
,
=1 ijii
jaa
th n
for =1, ,in and nmber of
discretization points.
We have the following assum
is the nu
ptions:
1) We assum
Assumption 1.
e
A
is positive definite, and therefore
we have
Copyright © 2011 SciRes. AJCM
168 J. GEISER

1
IA
se
1
(68)
e [19].
2) We assume
 

1
12 1
nn nnnn
z
yxiz
ECHCHJu JuE
 
 (69)
The stability is given with in the
Theorem 1. Given is the numerical scheme (70) and
e assumptions 1.
table for all iterative steps i.
Proof 1. Based on the assumptions we can boun
inverse matrix, also the previous solution is bou
Th
following Theorem:
we have th
The scheme is sd the
nded.
en scheme is stable.
We have the following proof idea:
Based on the assumption 1,
A
is positive definite
and the estimation of the remaining term, we have :
1
,.
nn
zi z
EE
(70)
So we have an upper bound of the iterative results,
gi
Convection-Diffusion Equation
Fo
sche
e in time.
ven by the previous solution at time n
t.
2.3.
Discretization Methods of the
r the 3 dimensional convection-diffusion equation we
apply a second order finite differenceme in space
nd a higher order discretization schema


2
2
22
22
=
,
xyz
vvvD
xyz
00
,= ,
uuu u
=,
uvu Du
t
x
uu
DD
yz


uxtu x

 
We apply dimensional splitting to our problem



=
x
yz
u
A
uAuAu
t

where
2
2
=.
xx
uu
Av D
x
x


We use a 1st order upwind scheme for
x
 and a
2nd order central difference scheme for 22
x
 . By
cing the artificial diffusion consintrodutant =
x
DD

2
x
vx we achieve a 2nd order fini
scheme
te difference
  

2
=
2
xx
x
Lu x vx
uxxux ux
D.
ux uxx
x
x

because t
or (i.e. the numerical viscosity) of the Taylor
expansion of the upwind scheme.
  
he new diffusion constant eliminates the first
order err
y
Lu and
z
Lu are
derived in the same way.
For the discretization in time we use several elicit
Runge-Kutta and Adam-Bashforth methods, this leads to
xp
restrictions of the step-size in time but on the other hand
the cost of implicit methods is much to high in this
3-dimensional case.
2.3.1. Ad a m-Bashfort h Methods

1
=0
=,
s
nn jnjnj
j
yyhbfty

(71)

 
1
0=0,
1
=
js
j
b
d
,=0,,.
!! uiujs
js j
(72)
s = 1 (first order)
ii
j
We consider here
 
11
31
=,,
22
nnnnnn
yyhfty fty





1
(73)
and
) s = 2 (second order
 

11
22
23 16
1
,
,
12 12
5
,
12
nn nn
nn
fty fty
ft y

2.3.2. Explic it Runge-Kutta Methods
In general a s-stage Runge-Kutta method can be written
in the following way:
j
k
=
nn
yyh

(74)
1
=1
=
s
nn j
j
yyhb
(75)
where
=1
=,
s
j
njn jll
kfthcyhak



(76)
l
We will take into accoun
Heuns third-order
t the following two:
(77)
and
Kuttas classical four th-order
Copyright © 2011 SciRes. AJCM
J. GEISER
169
(78)
3. Splitting Methods to Couple Maxwell and
Convection Diffusion Equation
We concentrate on the splitting methods, which can be
classified as classical and iterative splitti
We propose iterative splitting methods by discussing
the additive iterative splitting methods, see [20,21].
We consider the following the linear problem
(79)
rrs, e.g. they
orrespond in space to the discretized convection and
Thon with
fixed splitting discretization step size
ng methods.

=,
tctActBct
where the initial conditions are

=
nn
cct
. The opera-
s A and B are spatially discretized operatto
c
o
diffusion operators (matrices). Hence, they can be con-
sidered as bounded operators.
Iterative Splitting Methods
e following algorithm is based on the iterati
. On the time
ing subproblems interval we solve the follow
, cf.
1
,
nn
tt


consecutively for =1,3, ,2im,21].

1 [20
 

1
=, with=,
inn
ii i
ct
A
ctBctctc
t
(80)
  
1=,
i
ct

11
=,
nn
ii i
A
c tBct
with ctc

(81)
depends
t
where 00c and n
c is the known split approximation
at time level =n
tt. The split approximation at time
level =tefined as

11
22
=
nn
m
cct

. (Clearly,
the function
1i
c
e interval 1
t
1n
t
is d
t on th,
nn
t
,
we too, bof si
omit the dependence on n).
In the following we analyze the convergence and the
rate of the convergence of the method (80) - (81) for m
tending to infinity for the linear operators
ut for the sake mplicity, in our notation
,:
A
BX
rators and their sum are
X,
where we assume that these ope
ors of
ar
ce is enach spac
. Let uchy prob
generat the 0
C semigroups. We emphasize that
these operators necessarily bounded, thus the con-
vergenxamined in a general Bae setting.
Theorem 2s consider the abstract Cau-
lem in a Banach space X
en't
=
tctActBc

,0<tt T

0
0= ,
cc
,
(82)
where ,, :
A
BA BXX
are given lins
being generators of the 0
C semigroup and 0
cX
ear operator
is a
given element. Thethe iteration process (80) - (81) is
convergent and the rate of the convergence is of higher
order.
n
B are matrices (i.e. (80) - (81)
is a system of ordinary differential equations
growth estimation we can use the concept
rithmic norm, see e.g. [23]. Hence, for many important
of matrices we can p
hat portan
y of the split subproblems-the iterative splitting
mo the exact solution.
Fo
The proof can be found in [22].
Remark 6. When A and
), for the
of the loga-
classes rove the validity.
Remark 7. We note t a huge class of imt dif-
ferential operators generate a contractive semigroup.
This means that for such problems-assuming the exact
solvabilit
ethod converges in higher order t
In the next subsection we present the used time-dis-
cretization methods.
4. Error Analysis: Coupling Methods
r the coupling methods we deal with nonlinear differ-
ential equations of the following type:


 



d=,with=,
d
nn
c
A
ct ctBct ctctc
t(83)
where
=,,,
xyz
cHHEu
, with
x
H
,
y
H
is the mag-
netic field,
z
E is the
of the species.
electric fihe concen-
ation
eld and u is t
tr
The main idea is to bound the operators
A
ct and
Bct in the discretized equation able
hat is discussed in
the following subsection.
Iterative Operator-Splitting Method as a
int Scheme
e
earize the nonlinear operators, see
[2
We re
quations of the form:
to satisfy a st
method.
A first idea is the fix-point scheme, t
Fix-Po
The iterativoperator-splitting method is used as a
fix-point scheme to lin
1,24].
strict our attention to time-dependent partial dif-
ferential e


 



d=,=,
d
nn
u
A
ut utBut utwithutc
t(84)
where
,:
A
uBuXX are linear and den
fined in the real Banach space
sely de-
X
, involving only spatial
derivatives of c, see [25]. In the following we discuss the
standard iterative operator-splitting methods as a fix-
point iteration method to linearize the operators.
Copyright © 2011 SciRes. AJCM
J. GEISER
170
lit our nonlinear differential equation (84) by
ap
We sp
plying:
 

 



111
=,
with= ,
ii ii
nn
i
di
ut
A
ututButut
dt
utc

(85)
 

 


1
111
1
d=,
with= ,
iii ii
nn
i
ut
A
ututButu
dt
ut c

(86)
wm. is the starting solution,
e is near , or
al fix-po
tion at
xel
here the time step is 1
=nn
tt
. The iterations are
=1,3, ,i2 1
where we assum

0=0ut. So we
lem. n
c is th
level =n
tt.
The split appro

=n
c
the solution n
c
o solve the loc
split approxim
tion at time lev
0
ut
have t
e known
ima
1
a
n
c
int
the
1
prob-
time
=n
tt
is de-
fined as . We assumerators

11
22
=
nn
m
cut


11
,:
ii
the ope

A
uBuXX

the real Banach spac
to be linear and de
e X, for m
nse
1,3,, 2ily de-
fined on= 1
.
Here the linearization is done with respect to the itera-
tions, such that

1i
1
,
i
A
u
re at least non-de-
ive equations, and we can
apply the linear theory.
nearization is at least in the first equation
 
ii
Bu
ors in the iterat
a
pendent operat
The li
1
A
uAu
, and in the second equation
1i
Bu
Bu
1i
We have
 

11 11n nn
Au uu
1
n
ii
Autu t
linear
5. Experiments
In the following experiments, first we deal with the de-
coupled equations, means Maxwell and transport equa-
.
Maxwell equations in 2D is given
 

with sufficient iterations

=1,3, ,21im.
Remark 8. The th the fix-point scheme
can be used for smooth or weak nonlinear operators,
otherwise we loose the convergence behavior, while we
di
ization wi
d not converge to the local fix-point, see [21].
tions, to verify our methods
In the third experiment,ple we consider a sim PE-CVD
process and concentrate on the coupled transport and
Maxwell equation.
5.1. Test Experiment 1: Maxwell Equation
he time-dependent T
as:
  
,=,,, 0,,
xz
Hxy E
x
yt T
ty

 (87)
 
,=,,, 0,,
yz
HxyE
x
yt T
tx

 (88)
,1y
zx
 
,, 0, ,
=,
s
ource
J
tx
y




H
Exy H

(89)
e
xyt T
wher
,=sin
source
J
xyt .
have to implement the We outflow condition,
underlying discretization method (we assume fi
ference methods), means how many concentration is
flowing via the time-step
via the
nite dif-
t
to the cell with th
step
e spatial
x
:
The relative spatial step is given as
1=
relativ
tx

The percentage of the outflow is given as:
=
relativ
xrel
r
x

,=,,
z outz
ErelExy
The same is also given fothe ,
x
y
H
H.
Hee apply the FDTDethod of Yee’s algorithm. re w m
it is important to
balance suc
We assume to have finite difference schemes in time
an
evy) condi-
tion is important to balance the schemes:
While we are dealing with wave-equations:
For spatial and time discretization
h schemes.
d space.
Therefore the CFL (Courant Friedrichs L
x
t

here
x
, t
are the spatial anwd time steps.
fo
To control the electric field

,
z
Exy
, we have the
llowing line source:

,=sin

where =0,0,100
source
J
xy t
xy
The control of the particle transport is given by the
electric field in Figure 1. The electric and transport
sitnal model
in Figure 2.
In the following we have the line sources with the re-
FMaxwell
tric field in the reactor. We apply
Yees algorithm to obtain at least a second order scheme
in time and space. Based on the slower time-scales of the
Maxwell equations, which is less stiff than the transport
eq
uation is given with cut of the three dimensio
sults given in igure 3:
Remark 9. We consider the equation, that
models a periodic elec
uations, we have sufficient accuracy in the full coupled
system. A higher order discretization scheme is neces-
sary for the transport part.
Copyright © 2011 SciRes. AJCM
J. GEISER
171
Figure 1. Electric field in the apparatus.
Algorithm 4.4
1) Initialize Convection-Diffusion equation, till tstart.
2) Solve Electric Field equation with, we obtain s
ttart

,
z
E
xy for tstart.
3) Solve Convection Diffusion equation with and use
start
tt
,
z
E
xy
for t for the unknown.
4) Dotand go to 2.) till
start
=
start start
tt =
s
tart end
tt
Figure 2. Electric field in the apparatus.
Figure 3. Line source of the Electric field in the apparatus.
5.2. Test experiment 2: Convection-Diffusion
Equation
We deal with the 2-dimensional advection-diffusion
equation and periodic boundary conditions


22
22
00
=,
=,
,= ,
t
xy
uvuDu
uu u
vvDD
xy
u
x
y
uxtu x
 
 
 

with the parameters
1
0
==
=0.01
= 0.25.
xy
vv
D
t
The given advection-diffusion problem has an ana-
lytical solution
 
2
1
,=exp
x
4
a
vt
uxt 


ction:
tDt


which we will use as a convenient initial fun
00
,= ,
a
uxtuxt
We apply dimensional splitting to our problem
=
x
y
u
A
uAu
where
t
2
2
=.
xx
uu
Av D
x
x

x

We use a 1st order upwind scheme for and a
2nd order central difference scheme for 2
2
x
. By in-
troducing the artificial diffusion constant =
x
DD
2
x
vx we achieve a 2nd order finite
scheme
difference


2
=
xx
Lu x vx
2u
.
x
ux uxx
uxxx uxx
D

x
 
because therst
order error (i.e. the numical viscosity) of the Taylor
expansion of the upwind scheme.
new diffusion constant eliminates the fi
er
y
Lu is derived in the
ay.
We apply a BDF5 method to gai 5th order accuracy
in time. For simplifications, we no that the dependen-
cies of
same w
n
te
,uxt are suppressed as

ut.
 

1137
=55
60
1051
234.
t
Lututtutut
t
ut t uttut t
 
34
5
t

(90)
Copyright © 2011 SciRes. AJCM
J. GEISER
172
To compare the four methods we have the following
general setting. Let
=0,1 0,1 0,1
the initial concentration
, the unit
cube. There we set up
 
2
0=2exp 0.02
t
ux x





(91)

xa


(92)
T
with= 0.5,0.5,0.5a
which is just the analytical solution
 
2
1
,=exp 4
a
x
vt
uxt tDt


(93)
with and at on


=1v=0.01D0
==0.25tt
.
During the following experiments we will se and
consider an equidistant lattice of poi
t =0
nts (
v
3
N=
x

==1 1yz N ).
The result is shong Figures 4
and 5:
wn within the wi
Remark 10. We consider the transport equation, that
models the mass transport of the ion ized species fro
lower-left to the middle of the reactor. We use h
order time and spatial discretization schemes to obtain
higher order solutions. Such methods,
larger time and spatial steps and obta
Based on the fast time-scales of the transport equa-
tions, which is stiffer than the Maxwell equatio
balance the larger time-steps with suffic
solution of the transport regime in the coupled system.
5.
Electric Field Equation
H
ation, see citelieb05.
follo
m the
igher
we can apply with
in sufficient accu-
rate results.
n, we can
ient accurate
3. Test Experiment 3: Coupling
Convection-Diffusion and
s
(Weak Coupling)
ere, we consider a simple PE-CVD process, that an
underlying mass transport of a gaseous species is influ-
enced by an electric field, see [1,3].
For transport in a plasma environment, we assume a
homogeneous medium and that the influence of the elec-
tric field can be simulated by a coupled transport and
Maxwell equ
For simplifications, we deal with the 2-dimensional
advection-diffusion equation and electric field equation:




22
22
00
,
,,=,,
=,
txz y
uu
uvExy v
x
y

 

uu
DD
xy
uxytu xy




 
,Hxy E


 
,, 0, ,
=,,, 0,,
,=,,,
xz
yz0,,
,1
=,
y
zxsource
x
yt T
ty
HxyExyt



T
tx
H
Exy HJ
txy






xyt T
The advection-diffusion problem has an analytical so-
lution at the beginning for

00,
s
tart
tt
 
2
1
,=exp 4
a
x
vt
uxt tDt





which we will use as a convenient initial function:
00
,= ,
a
uxtu xt
Further the function:


(, )=1for0,
(, )=(, )for
xz start
x
zz
vExytt
vExy Exytt
start
where = 0.001
, = 10.0
start
t.
Figure 4. Initial gaseous concentration at t = 05.
.2
Figure 5. Gaseous concentration at t = 0.5.
Copyright © 2011 SciRes. AJCM
J. GEISER
173
Both equations have the same domain
=0,1
0,1
Nu
.
merically we solve the equation, as in th
algorithm 5.3:
The following figures show the developing of the
concentration under the influence of the electric field, we
deal with a normalized time scale in
e following
s
ec . Further we
have =0.07
, and r =0.5
start
t=
y
v0 fo
s
tart
tt.
Figure
Remark 11. Based on the transport of the ionized spe-
cies from the lower-left to the middle of the reactor, we
see and influence of the species. The former circular
concentration is spread out to a diffusive ellipseere,
we ctric
eld
time- and spatial scales of the underlying transport
and Maxwell equation. Via iterative splitting, we could
couple the two equations systems together and reduce
the numerical errors with additional iterative step s.
The results are given in 6 - 11.
. H
an control the species in the reactor with a elec
. Numerically, it is important to deal with the differ-fi
ent
Figure 6. Gaseous concentration after t = 0.833.
Figu
in r first
fluence of the electric field.
e 8. Gaseous concentration after t = 1.162 with a
Figure 9. Electric field after t = 1.162.
Figure 10. Gaseous concentration after t = 1.483 with a first
influence of the electric field.
Figure 7. Electric field after t = 0.833.
Copyright © 2011 SciRes. AJCM
174 J. GEISER
Figure 11. Electric field after t = 1.483.
6. Conclusions
We present a coupled model based on Maxwell and
Transport equations, that can be applied for simplified
transport model for an ionized gaseous species in a
PECVD reator. Based the different scale models, we
have included the optimal discretization methods for
each separate equation. Splitting methods are used to
couple the separate equations together. Further, we d
cussed the splitting analysis. Numerical examples are
presented to discuss the influence of decoupled and cou-
pled systems. In future, we will analyze the validity of
the models with physical experiments.
7. References
ience of Thin Films,” 2nd Edi-
tion, Academic Press, San Diego, New York, Boston,
London, 2002.
[2] J. Geiser and M. Arab, “Simulation of a Chemical Vapor
Deposition: Mobile and Immobile Zones and Homoge-
neous Layers,” Journal of Porous Media, Begell House
Inc., Redding, 2009, Vol. 1, No. 2, pp. 123-143.
[3] L. Rudniak, “Numerical Simulation of Chemical Vapour
Deposition Process in Electric Field,” Computers &
Chemical Engineering, Vol. 22, Supplement 1, 1998, pp.
755-758.
[4] J. Van Roey, J. van der Donk and P. E. Lagasse,
“Beam-Propagation Method: Analysis and Assessment,”
Journal of the Optical Society of America, Vol. 71, No. 7,
1981, pp. 803-810. doi:10.1364/JOSA.71.000803
is-
[1] M. Ohring, “Materials Sc
[5] M. D. Feit and J. A. Fleck Jr., “Analysis of Rib
Waveguides and Couplers by the Propagating Beam
Method,” Journal of the Optical Society of America A,
Vol. 7, No. 1, 1990, pp. 73-79.
doi:10.1364/JOSAA.7.000073
[6] M. D. Feit and J. A. Fl
Graded-Index O
eck Jr., “Light Propagation in
ptical Fibers,” OSA Applied Optics, Vol.
17, No. 24, 1978, pp. 3990- 3998.
doi:10.1364/AO.17.003990
[7] L. Thylen, E. M. Wright, G. I. Stegeman, C. T. Seaton
and J. V. Moloney, “Beam-Propagation Method Analysis
of a Nonlinear Directional Coupler,” OSA Optics Letters,
Vol. 11, No. 11, 1986, pp. 739-741.
doi:10.1364/OL.11.000739
[8] J. A. Fleck, J. R. Morris Jr. and M. D. Feit, “Time-de-
pendent Propagation of High Energy Laser Beams
through the Atmosphere,” Applied Physics, Vol. 10, No.
2, 1976, pp. 129-160. doi:10.1007/BF00896333
[9] M. Lax, J. H. Batteh and G. P. Agrawai, “Channeling of
Intense Electromagnetic Beams,” Journal of Applied
Physics, Vol. 51, No. 1, 1981, pp. 109-125.
328442doi:10.1063/1.
0] G. I. Marchuk, “Some Applications of Splitting-up Meth-
and Comparion of Dif-
[1
ods to the Solution of Mathematical Physics Problems,”
Aplikace Matematiky, Vol. 1, 1968, pp. 103-132.
[11] G. Strang, “On the Constraction
ference Schemes,” SIAM J. Numerical Analysis, Vol. 5,
No. 3, 1968, pp. 506-517. doi:10.1137/0705041
[12] K. Okamoto, “Fundamentals of Optical Waveguides,”
Academic Press, New York, 2005.
[13] A Taflove, “Computational Electrodynamics: The Finite
Difference Time Domain Method,” Arctech House Inc.,
1995.
[14] J. P. Berenger, “A Perfectly Matched Layer for the Ab-
sorption of Electromagnetic Waves,” J. Comp. Phys., Vol.
111, 2005, pp. 185-220.
[15] S. D. Gedney, “An Anisotropic Perfectly Matched
Layer-Absorbing Medium for the Truncation of Fdtd Lat-
tices,” IEEE Tran. Ant. Prop., Vol. 44, No. 12, 1996, pp.
1630-1639.
[16] W. Sha, X. Wu, M. Chen and Z. Huang, “Application of
the High-Order Symplectic Fdtd Scheme to the Curved
ymplectic Integrator
.
Splitting
and C. Kelley, “Convergence of
Three-Dimensional Perfectly Conducting Objects,”
[17] T. Hirono, W. Lui, S. Seki and Y. Yoshikuni, “A Three-
Dimensional Fourth-Order Finite-Difference
Time-Domain Scheme Using a S
Propagator,”
[18] J. Geiser, “Numerical Simulation of a Model for Trans-
port and Reaction of Radionuclides, 2001,” Proceedings
of the Large Scale Scientific Computations of Engineer-
ing and Environme n t a l P r o b l e ms, Sozopol, 2001
[19] W. Hackbusch, “Iterative Losung Groser Schwachbe-
setzter Gleichungssysteme,” Teubner-Verlag, Stuttgart,
1993.
[20] I. Farago and J. Geiser, “Iterative Operator-
Methods for Linear Problems,” International Journal of
Computational Science and Engineering, Vol. 3, No. 4,
2007, pp. 255-263.
[21] J. Kanney, C. Miller
Iterative Split-Operator Approaches for Approximating
Copyright © 2011 SciRes. AJCM
J. GEISER
Copyright © 2011 SciRes. AJCM
175
Nonlinear Reactive Transport Problems,” Advances in
Water Resources, Vol. 26, 2003, pp. 247-261.
doi:10.1016/S0309-1708(02)00162-8
[22] J. Geiser, “Weighted Iterative Operator-Splitting Methods:
Stability-Theory,” Proceedings of the 6th International
Order Time-Integration Methods and Applications
ear Functional Analysis and Its Ap-
for P
Conference, NMA 2006, Lecture Notes in Computer Sci-
ence, Springer, Berlin, Vol. 4310, 2007, pp. 40-47.
[23] W. H. Hundsdorfer and J. G. Verwer, “Numerical Solu-
tion of Time-Dependent Advection-Diffusionreaction
Equations,” Springer, Berlin, 2003.
[24] J. Geiser, “Iterative Operator-Splitting Methods with
Higher
arabolic Partial Differential Equations,” Journal of
Computational and Applied Mathematics, Elsevier, Am-
sterdam, Vol. 217, 2008, pp. 227-242.
[25] E. Zeidler, “Nonlin
plications II/B: Nonlinear Montone Operators,” Springer-
Verlag, Berlin-Heidelberg-New York, 1990.