Journal Menu >>

 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 2222=, ,txz yuuuvExyv DxyuxuDy  (1) 00,,=, ,uxytuxy (2)   ,=,,, 0,xzHxy E,xyt Tty (3)   ,=,,, 0,yzHxy E ,1=,,,0, ,yzxsourceHExy HJtxyxyt T (5) where u is the concentration of the gaseous species, zE is the electric field and ,xyHH is the corresponding magnetic field in two dimensions. Further txyv=,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 ,xyt Ttx (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,Ert (6) 0=E,Hrt (7) where is the electric field, =,, ,,ExyzEEE xyt=H ,, ,,xyzHHH xyt is the magnetic field, =,rxy=1r 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 =0Er=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 =0t0=yxzEHEtyz (8) 0=yxzHEEtzx (9) 0=yxzrHHEtxy (10) Let x, y are spatial discretizations, and t is a time step. We use the following notation  ,=, ,n.Fij 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 yE of E and H of H at the following spatial locations and temporal instants: ::=:spatial coordinateinteother spatial coordinategtime integerhalfte in=:ntegergererE (12) ::spatial coordinateiHotherspatial coordinatalfertimehalf 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: ,z121201,211=, ,1,2xzHijt nxnnnHijE ijijy    E(17) ,z121201,211=,1 ,2yzHi jt nyn,nnHijEij ijx    E(18) 1112,,2011220,=,1112211,,22nnyyrnnxxrEijEijt nznz1.HijHijxtHij Hijy (19)    Copyright © 2011 SciRes. AJCM J. GEISER 165In 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  211 1tc2xy  (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,Auftt (21) =nnuut, (22) where A is the operator that we implicit solve in the equation and  =ftBut 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: 121=1 =I,Htx, =yxzxyy xHEJyHHJBH BHJ (23) 11==II=xz1zzHEECtyE (24) 21== III=y1zzzHECEEtx (25) For the boundary conditions we assume periodic boundary conditions. That means we use the identifica-tion ,= 1,zzENi Ei (26) ,=,1 =1,,zzEiN 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. 20022211=0=0==rrct c □ (28) 22 222221=0xyct   (29) 22221=0xy tDD Dc  (30) 2222221c11 1=0yxttyxttDcDD DDcDDcD  (31) Copyright © 2011 SciRes. AJCM 166 J. GEISER 22()1111xt xtDDVDDVcc  =0 (32) =0xx□□ (33) Waves that satisfy only propagate in =0x□x-direction and those that satisfy only propagate in =0x□x-direction. An analogous formulation can be given for the and direction. yyTo handle  it is comfortable to do a Taylor ex-pansion around 0. 232223452211=1 00121 01VVVVVVVOVV 2V (34) 21 =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=xx) For the left boundary we have do discretize the fol-lowing equation: 1=.xct  (38) This can be done with a FDM-scheme as follows.  11221=02112=1=,2133=,,22nnxx xnnnttjjxxjjtt     ,1;(39) with     1121121,2 =,2,2;21,1 =,1,12nnnnnnjjjjjj (40) and   11131,=,2 ,12231,=,2,122nnnnnnjjjjj;j (41) this leads to   ,1n11,1 =,2,2nn nct xjjjct 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 tbecause it onas to be applied only to the zE field because xH and yH 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 ofe on the next simulations with less syis discretised and is calcu-lammetry 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  =iijjtt tc. t Let The ste from to pit1it ng way. in (23) - (written i25) can now be n the followi 1II112 233=iizzEEtbkbkbk  I (44) II1II II112233=iixxHHtbkbk )  bk (45III III III2 23 3111=iiyyHHtb kbkbk    (46) where iiiiM=,,,jjxyzjjjk MHHEJI, II, III M for and 1, 2,3j. With  3I=1==iiizjz jllzjlEEtE tak   (47)  II=1==iiix lxjl3xj jlHHt Ht  ak)  (48 3II=1==iiiyj yjllyjlIHHtHt ak  (49) Copyright © 2011 SciRes. AJCM J. GEISER 167For a better legibility and because the focused point of time does not change, we write =which is known(in our caiijjJJtse) (50) ,,,jjjzxyEHH - (50) give 9 equatijJ in-stead of Comb 7) ,,,iiiijzxyjjjEHHJ . ining (23) - (25) and (4ons 12=1=jlllzzjy xlEEtal BHBHJ   (51) 3131=1=1jlxxjzlHHtalCE  (52 ) 32=11=1jlyyj zlHH talCE j ,2,3 Remark 3. (53) 1 and 1 are trices, such that also realized as ma-1=1 ,xy and 1=1,xy. Remark 4. The scheme above is only correct for iso-tropic media because in the not isotropic case iessary to consider t is nec-=,xy. of eqith equatstTaking 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. 1QE23zzzQEQER11 12 132122233132 33123,,,=,,,,,, AA AA AAAA AA AAAA AA AAzzzCSIRCSRC SRCSRCSI RCSRC SRC SRC SIEEE (54) 123122211 1213222221 2223322231 3233=zzzzzzQEaSIaS aSEQEaSaSIaSEQEaSaSaSIE     (55) 21122300 00=000000 003zzzzzzQESI EQESI EQESIE       (56) 121:= ,1,jA AjyxjQtRBHBHRJ  (57) A (59) := -AjRjthrowof (58) := -AjCjthcolumnofA1:= 1,1,1 (60) 123:=, ,JJJJ (61) 221BC (62) 2111:=St BC22:==,AAi jRC (63) ij ijaA (65) 2211 1213222221 222322 231 3233:=aaaa aaa a 2a (64) 2220=0ijijijnmtimesanmtimesaa :=IIdentity nm (66) wher and e ,ija ibjc selare the Runge-KutTopreciy: 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 toRemark 5. For an optimiz, n of DIRK fastert reducing th n sco 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 , withoue accuracy. For higher time-discretizations we have to taken into account also higher spatial discretizatioheme. 2.2.3. Stability Analysis of the Implicit Discretizations We deal with the following discretized equation systems:  11 1,121=nnnnnnziz y xiEIAECHCHJuJu  (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 wheretes the transpose of z. T>0zAznz,non-zero vectTz 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 ijiijaath n for =1, ,in and nmber of discretization points. We have the following assum is the nuptions: 1) We assumAssumption 1. e A is positive definite, and therefore we have Copyright © 2011 SciRes. AJCM 168 J. GEISER 1IAse1 (68) e [19]. 2) We assume  112 1nn nnnnzyxizECHCHJu 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 bouninverse matrix, also the previous solution is bouThfollowing Theorem: we have thThe 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,.nnzi zEE (70) So we have an upper bound of the iterative results, gi Convection-Diffusion Equation Fo schee in time. ven by the previous solution at time nt. 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222222 = ,xyzvvvDxyz00,= ,uuu u=,uvu DutxuuDDyzuxtu x We apply dimensional splitting to our problem  =xyzuAuAuAut where 22=.xxuuAv Dxx We use a 1st order upwind scheme for x and a 2nd order central difference scheme for 22x . By cing the artificial diffusion consintrodutant =xDD 2xvx we achieve a 2nd order fini scheme te difference  2=2 xxxLu x vxuxxux uxD.ux uxxxxbecause tor (i.e. the numerical viscosity) of the Taylor expansion of the upwind scheme.    he new diffusion constant eliminates the first order erryLu and zLu are derived in the same way. For the discretization in time we use several elicit Runge-Kutta and Adam-Bashforth methods, this leads to xprestrictions 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=,snn jnjnjjyyhbfty (71)  10=0,1=jsjbd,=0,,.!! uiujsjs j (72) s = 1 (first order) iijWe consider here  1131=,,22nnnnnnyyhfty fty 1 (73) and ) s = 2 (second order 112223 161,,12 125 ,12nn nnnnfty ftyft y 2.3.2. Explic it Runge-Kutta Methods In general a s-stage Runge-Kutta method can be written in the following way: jk=nnyyh (74) 1=1=snn jjyyhb (75) where =1=,sjnjn jllkfthcyhak (76) lWe will take into accounHeun’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 splittiWe 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 =nncct. The opera-s A and B are spatially discretized operattocodiffusion 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,nntt consecutively for =1,3, ,2im,21]. 1 [20 1=, with=,innii ictActBctctct (80)   1=,ict11=,nnii iAc tBctwith ctc(81) depends twhere 00c and nc is the known split approximation at time level =ntt. The split approximation at time level =tefined as 1122=nnmcct. (Clearly, the function 1ice interval 1t1nt is dt on th,nnt, we too, bof siomit 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 ,:ABXrators and their sum areX, where we assume that these ope ors ofarce is enach spac. Let uchy probgenerat the 0C 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