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 u Dy (1) 00 ,,=, ,uxytuxy (2) ,=,,, 0, xz Hxy E , yt T ty (3) ,=,,, 0, yz Hxy E ,1 =, ,,0, , y zx ource H Exy HJ txy xyt T (5) where u is the concentration of the gaseous species, E is the electric field and , y 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, rt (6) 0= E, rt (7) where is the electric field, =,, ,,Exyz EEE xyt = ,, ,, xyz HH xyt is the magnetic field, =, r 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= xzE HE ty z (8) 0= yx HEE tz x (9) 0=y z r H E tx y (10) Let , are spatial discretizations, and t is a time step. We use the following notation ,=, , n. ij Fixjynt (11) Let represents a spatial coordinate such as , . The goal of Yee’s scheme is to compute the approxima- tions for the various components y E of and of at the following spatial locations and temporal instants: : : =: patial coordinateinte other spatial coordinateg time integer half te in =: nteger ger erE (12) : : spatial coordinatei 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 ijE ijij y E (17) , z 1 2 1 2 0 1, 2 11 =,1 , 2 yz Hi j t n y n, nn 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 . 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 Yee’s algorithm: ● The CFL stability condition for the Yee’s FDTD method is 2 11 1 tc2 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 Yee’s 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, uft t (21) = n n uut, (22) where is the operator that we implicit solve in the equation and = 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 HJBH BHJ (23) 1 1 ==II= xz1 z HEEC ty E (24) 2 1 == III= y1 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 □ -direction and those that satisfy only propagate in =0 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 =12VOV4 (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) For the left boundary we have do discretize the fol- lowing equation: 1 =. 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 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 E field because and 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 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 Htbkbk ) bk (45 III III III 2 23 3 1 11 = ii yy 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 jz jll zjl EEtE tak (47) II =1 == iii x l xjl 3 xj jl Ht Ht ak ) (48 3II =1 == iii yj yjll yjl I 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) ,,, 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 l xjz l 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 , y and 1=1, 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- =, 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 QEaSIaS aSE QEaSaSIaSE QEaSaSaSIE (55) 2 11 22 3 00 00 =0000 00 003 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 :=, , 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 := Identity nm (66) wher and e , ij a i b 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 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 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, 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 uu DD yz uxtu x We apply dimensional splitting to our problem = yz u uAuAu t where 2 2 =. xx uu Av D x We use a 1st order upwind scheme for and a 2nd order central difference scheme for 22 . 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 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 Lu and 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 njn jll kfthcyhak (76) l We will take into accoun Heun’s third-order t the following two: (77) and Kutta’s 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 ctBctctc t (80) 1=, i ct 11 =, nn ii i 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 ,: 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 ,, : 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 ct ctBct ctctc t(83) where =,,, xyz cHHEu , with , is the mag- netic field, E is the of the species. electric fihe concen- ation eld and u is t tr The main idea is to bound the operators 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 ut utBut utwithutc t(84) where ,: uBuXX are linear and den fined in the real Banach space sely de- , 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 ututButut dt utc (85) 1 111 1 d=, with= , iii ii nn i ut 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 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 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 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 yt T ty (87) ,=,,, 0,, yz HxyE yt T tx (88) ,1y zx ,, 0, , =, ource J tx y H Exy H (89) e xyt T wher ,=sin source 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 : 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 , y 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 t here , 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 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 Yee’s 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 xy for tstart. 3) Solve Convection Diffusion equation with and use start tt , z xy for t for the unknown. 4) Dotand go to 2.) till start = start start tt = 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 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 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 = y u uAu where t 2 2 =. xx uu Av D x We use a 1st order upwind scheme for and a 2nd order central difference scheme for 2 2 . 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 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 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= ==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 y uu DD xy uxytu xy ,Hxy E ,, 0, , =,,, 0,, ,=,,, xz yz0,, ,1 =, y zxsource 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, tart tt 2 1 ,=exp 4 a vt uxt tDt which we will use as a convenient initial function: 00 ,= , a uxtu xt Further the function: (, )=1for0, (, )=(, )for xz start 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 ec . Further we have =0.07 , and r =0.5 start t= y v0 fo 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.
|