SIMULATION OF SEMICONDUCTOR DEVICES AND PROCESSES Vol. 4 Edited by W. Fichtner, D. Aemmer - Zurich (Switzerland) September 12-14,1991 131 - Hartung-Gorre Transformation Methods for Nonplanar Process Simulation K a r l W i m m e r , R o b e r t B a u e r , Stefan H a l a m a , G e r h a r d Hobler*, a n d Siegfried S e l b e r h e r r *Institut Institut fur Mikroelektionik, far Allgemeine Elektrotechnik und Elektronik, Technical University of Vienna, A u s t r i a Abstract We present a transformation method for 2D process simulation in arbitrary structures. Several grid generation techniques have been investigated systematically on their influence on convergence properties. A mapping function strategy applicable to process simulation problems is introduced. New formulations of the discretized equations based on box integration are introduced in order to satisfy global conservation laws automatically. 1 Introduction Transformation methods are appropriate for the numerical simulation in nonplanar 2D and 3D domains [1], [2], a n d are widely used in computational fluid dynamics [3]. The transformation is accomplished by specifying a curvilinear coordinate system which m a p s the nonplanar physical domain (a;, y) onto a rectangular computational domain (u, v) (see Fig. 1). W i t h this approach the simulation domain is simplified at the expense of complicating the equations and boundary conditions. North r{u,vN) l* c ^f" North East r(u„,v) v=vN 3 3 y» CO G O \» Sc utl 1 r(u,' l ) Wes West r(upv) South u Figure 1: Transformation from physical space (a;, y) to computational space (u, v) by a curvilinear coordinate system For application in our 2D process simulator P R O MIS [4], the methods must be suited for moving boundary problems. Additionally, an unknown variety of geometries has to be treated a n d adaptive transient grid strategies with heavy changes in grid spacing must be supported. The methods presented in this paper fulfill the above requirements. 132 2 Grid Generation and A d a p t i o n The generation of the grid used for the solution of the physical equations is separated into two tasks. The first task is the generation of a mapping function. An approximately equidistant grid in the computational domain which resolves all geometric details (Fig. 2a) is mapped onto the physical domain (Fig. 2b). This grid ("geometric grid") is designed such that the points at the boundaries are approximately equidistributed along the arc length. It is only used to establish the mapping function (u,v) +-+ (x,y), and has to be calculated just once for each geometry. The second task is the generation and adaption of the grid for solving the physical equations ("physical grid"). Two criteria control the adaption of the grid according to the evolving dopant profiles. A gradient criterion guarantees the resolution of steep gradients in the dopant profiles. A dose conservation criterion minimizes the local dose error (1), which is expressed consistently with the discretization of the diffusion equation. After each timestep during a transient simulation the adaption is achieved by inserting and deleting grid lines in the computational domain (Fig. 2c) based on the above criteria. The corresponding grid lines in the physical domain are obtained by interpolation in the geometric grid (Fig. 2d). geometric grid physical grid Figure 2: Grid generation and adaption: The "geometric grid" (left) is used to establish the mapping function from computational domain (a) to physical domain (b). The "physical grid" (right) is used for solving the physical equations. For the generation of the mapping function, algebraic, elliptic or variational methods can be used. Algebraic methods are based on an interpolation between the boundaries [5]. Using the boundary points and, if necessary, positions at internal interfaces, e.g. Si/SiO^, we get the positions of the inner grid points from a transfinite interpolation scheme (2) with Lagrange polynomials (3). 133 Boundaries: Interfaces: South: f[u,vi), North: f ( u , u ^ ) , West: r(ui,v), f(«/,v), 1<I<M, r{u,vj), 1<J<N M N r M v East: f{uM,v) N f f{u, v) = Y, *.• (« ) • "T«i. ) + Z) *,• (« ) • («, «i) - ] £ S *<(«) • *;(«) •tf(«i>v i) «=i j=i (2) «=i j = i «s" Mi Elliptic methods are based on the solution of the boundary value problem [2] Au = P, Av = Q. In the computational domain, the corresponding transformed equations (4)-(5) have to be solved. Subscripts denote partial derivatives, and g = xuyv~xvyu is the Jacobian determinant. a • xuu - b • xuv + c • xvv = -g • (P • xu + Q • xv) (4) a • Vuu ~ b • yuv + c • Vn, = -g • (P • yu + Q . yv) (5) a = *l + vli b = *««» + yuVv, c = xl + yl (6) This system of nonlinear elliptic PDEs is discretized on the "geometric grid" using 9-point finite differences. The resulting system of coupled nonlinear algebraic equations is solved with a nonlinear SOR Algorithm [6]. A superimposed iteration scheme (7)-(9) determines the source functions P and Q at the boundaries in order to avoid clustering of gridlines at concave corners, and to allow orthogonality control at the boundaries. P° = Q° = 0 (7) pk+i = pk± aj-ctan g fc + 1 = Qk ± € , • arctan Here a denotes the angle of intersection of gridlines (10) (which should be equal to the required value ar = TT/2), and 6 resembles the gridspacing (11). a = arccos S = \f[u,vi+1) - r[u,Vi)\ (11) The upper signs in (8) and (9) relate to the South boundary, the lower signs to the North boundary. For the East and West boundaries, the correction terms for Pk and Qk are exchanged. In the interior, P and Q are obtained by transfinite interpolation. Variational methods offer a direct influence on certain grid properties [7]. We minimize a linear combination of three integrals (12) to control smoothness 1$, cell area IA and orthogonality Io of the numerically generated grid. / = A s • h + Ao • Io + *A • I A As > 0, \0> 0, A,i > 0, As + A 0 + A,i = 1 (12) (13) 134 IS = J((Vx)2 + (Vy)2)/gdudv IA = - 1 / 2 J Io = - 1 / 2 fw(xuxv g2dudv + yuyv)2dudv (14) (15) (16) In the orthogonality functional Io a weighting function w(u, v) enables additional orthogonality control at the boundary. The minimization problem for the overall criterion J is solved by calculating the EulerLagrange (EL) equations for the variational problem. The EL equations are given in the appendix. They are discretized using 9-point finite differences, the nonlinear algebraic system is solved with a nonlinear SOR algorithm [6]. Fig. 3 shows grids for a test geometry, generated with algebraic, elliptic and variational methods. Figure 3: Grids for a test geometry generated with different methods: left: transfinite interpolation middle: elliptic grid with source function P = Q = 0 right: variational method, \s = 0.3, \A = 0.7, ^o = 0; w = 10 at the boundaries, w = 1 else. 3 Transformation of Physical Equations The general structure of the diffusion equations (17)-(18) and the boundary conditions (19) in PROMIS allows the treatment of a wide variety of diffusion phenomena, including point defect assisted diffusion under oxidizing conditions [8], clustering effects [9], point defect and pair formation kinetics [10], [11], and grain boundary diffusion in polysilicon [12]. N (17) j=l 0 t N Ji = E (°*i • S"" 1 Cj + bij • Cj • grad * + c-,- • C}) + d; (18) J=I N (19) 135 A simple mathematical transformation from cartesian coordinates («, y) to curvilinear coordinates (u, v) using the Jacobian matrix of the transformation leads to discretized equations which violate global conservation laws considerably. By means of box-integration we have derived a conservative form of the transformed equations (20)-(23). They satisfy global conservation properties automatically without any additional computational expense. The "conservative" transformation of the equations from physical to computational domain is performed automatically in the program. In the most general case of a time variant physical domain (moving boundaries) we have to insert (20)-(23) into (17) to get the transformed equations in the fixed computational domain. dCj dCj (*.») ,v) - ^ • ((y»C-,-)« - (VuCj)v) - 2f • UxuCj)v 9 \ / 9 V div Ji = I • ((yvJ? - xvJ?)u + {-yuJ? + - (*„<?;)„) / (20) (21) »MJf)v\ N •? = E [ y ((*<?A. - (2/uCi)„) + ^ i ( ( y « * i ) u - (yu*,),,) + c%C^ + d? + <% (22) (23) The boundary conditions for the North and South boundaries (v = const) and West and East boundaries (u = const) are transformed according to the normal vectors. Jj • n = ± Jj • n = ± 1 VW+vl 4 (xujy - yuJJ) North, South (24) • (yvJJ - zvJ?) West, East (25) Evaluation of Grid Generation Strategies A recessed local oxidation serves as a test example. A structure with 0.3 fim etched recess, 20 nm native oxide and a 0.1 fim nitride mask is oxidized for 90 min at 1000°C in wet ambient. Fig. 4 shows the initial structure and the final oxide shape with an algebraically generated grid. The same structure with a grid generated by a variational method (\A = 0.5, \o = 0.2, Xs = 0.3) is shown in Fig. 2d. Figure 4: Recessed LOCOS: initial structure (left), and final oxide shape after 90 min oxidation at 1000°C in wet ambient (right). 136 The quality of the generated grid is judged upon the convergence properties of the solution of the discretized physical equations. The time for the grid generation is negligible in all cases. In Fig. 5 the relative computation time for the solution with respect to the algebraic case is shown. The elliptic grid saves about 50% of cpu time, if at least 5 iterations are used. With suitable parameters A .4, A<p, A5 a variational grid saves up to 70%. 100% 90% 01 1 ) 1 ! 1 i i i i _ E '3 80% d o 70% - "3 80% - 1UUA _ Xs=0.1 - X A +X 0 =09 - variational method 90% a> § B0% d o 70% -£ 80% (c) -*? a- 50% - ( W _ _ _ ^ E - O- 50% 6 § 40% 01 > 30% as *5 20% . («) w=l=const _ 10% - (b) w=10 at the boundaries - g v > s . l b ) ^ ^ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 a r e a weight XA 40% 30% cd % 20% elliptic method with tupiiimpotoil itertUae ichema for th* M>urce function* P, Q 10% 1 2 3 4 5 s o u r c e function 6 7 8 9 iterations Figure 5: Relative computation time required for solving the physical equations on the different grid types (reference: algebraic grid S100%): (a) variational grids with different area weighting parameters \A (b) same as (a) except for an additional orthogonality at the boundaries. (c) elliptic grids obtained after different numbers of P, Q iterations. From Fig. 5 and other examples it can be seen that orthogonality at the boundary and area control are most important, followed by some smoothness control, and then orthogonality in general. Note that, historically, orthogonality has received the most attention. 5 Conclusion Transformation methods have been applied successfully to process simulation problems. The computational expense for the generation of a "high quality" grid is typically only a few seconds on commonly used workstations and results in a speed up as high as a factor 3 for the solution of the transformed physical equations. An integral attempt for the transformation of the physical equations lead to conservative formulations of the discretized equations, and therefore global conservation properties are satisfied automatically without any additional computational effort. 137 Acknowledgement This project is supported by the research laboratories of: AUSTRIAN INDUSTRIES - AMS International at Unterpremstatten, Austria; DIGITAL EQUIPMENT Corporation at Hudson, USA; SIEMENS Corporation at Munich, Germany; and SONY Corporation at Atsugi, Japan. References [1] H. Umimoto, and S. Odanaka, Three-Dimensional Numerical Simulation of Local Oxidation of Silicon, IEEE Trans. Electron Devices, ED-38, pp. 505-511, 1991 [2] A. Hilgenstock, A Fast Method for the Elliptic Generation of Three-Dimensional Grids with Full Boundary Control, On Numerical Grid Generation in Computational Fluid Dynamics, S. Senguta, J. Hauser, P.R. Eiseman, and J.F. Thompson (Eds.), Pineridge Press, Swansea, U.K., pp. 137-146, 1988 [3] J.F. Thompson, A.U.A. Warsi, and C.W. Mastin, Numerical Grid Generation, Foundations and Applications, North Holland, New York, 1985. [4] P. Pichler, W. Jiingling, S. Selberherr, E. Guerrero, and H. Potzl, Simulation of Critical IC-Fabrication Processes, IEEE Trans. Electron Devices, ED-32, pp. 1940-1953, 1985 [5] P.R. Eiseman, and G. Erlebacher, Grid Generation for the Solution of Partial Differential Equations, ICASE Report No. 87-57, NASA Langley Research Center, Hampton, Virginia, 1987 [6] J.M. Ortega, and W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, San Diego, 1970. [7] J.U. Brackbill, and J.S. Saltzman, Adaptive Zoning for Singular Problems in Two Dimensions, J. Comp. Phys., 46, pp. 342-368, 1982 [8] M.D. Giles, Defect-Coupled Diffusion at High Concentration, IEEE Trans. Computer-Aided Design, CAD-8, pp. 460-467, 1989 [9] M.Y. Tsai, F.F. Morehead, J.E.E. Balgin, and A.E. Michel, Shallow Junctions by HighDose As Implants in Si: Experiments and Modeling, J. Appl, Phys., 51, pp. 3230-3235, 1980. [10] T.L. Crandle, W.B. Richardson, and B.J. Mulvaney, A Kinetic Model for Anomalous Diffusion during Post-Implant Annealing, IEEE IEDM 1988 Technical Digest, pp. 636-639, 1988 [11] G. Hobler, S. Halama, K. Wimmer, S. Selberherr, and H. Potzl, RTA-Simulations 2-D Process Simulator PROMIS, Nupad III Technical Digest, pp. 13-14, 1990 with the [12] F. Lau, Modeling of Polysilicon Diffusion Sources, IEEE IEDM 1990 Technical Digest, pp. 737-740, 1990 138 Appendix: Euler-Lagrange Equations The Euler-Lagrange Equations for the smoothness functional Is (14) are (26)-(27) using the abbreviations (28)-(29). For A2 - BC ^ 0 they reduce to (30)-(31). B{dxuu - 2exuv + fxvv) - A{dyuu - 2eyuv + fyvv) = 0 (26) A{dxuu - 2exuv + fxvv) - C{dyuu - 2eyuv + fyvv) =0 (27) A = xuyu + xvyv 2 d = (x v + y2v)/g3 C = x2u + x2v B = y2u + y2 2 3 f = (x u + e = (xuxv + y»yv)/g y2u) /g 3 (28) (29) dxuu - 2exuv + fxvv = 0 (30) dyUu - ley™ + fyw = 0 (31) As Euler-Lagrange equations for the area functional IA (15) we obtain (32)-(33). ( z u - xv)g + gvxu- guxv = 0 (32) (yv - yu)g + guyv - gvyu = o (33) The Euler-Lagrange Equations for the orthogonality functional with included boundary control Io (16) are (34)-(35) using the substitutions (36)-(40). w • (hxuu + b2xuv + b3xvv + a i ^ u u + a2yUv + a33fo») = ~n w • (aiz u u + a2xuv + a3xvv + cii/ uu + c2yuv + c3yvv) = - r 2 <*i = xvyv a2 = xuyv + xvyu 03 = xuyu h = x2v b2 = 2(2xuxv + yuyv) h = x2 63 = xu u (34) (35) cx = y2 c2 = 2(xuxv + 2yuyv) (36) (37) c _= yl2 (38) c 33 - yu ri = (wuSo - wvyu) • (xuxv + yuyv) /2g r2 = (wvxu - wuxv) • (xuxv + yuyv)2/2g (39) (40)