NUMERICAL SIMULATIONS OF SEMICONDUCTOR DEVICES BY STREAMLINE-DIFFUSION METHODS By Xunlei Jiang B.Sc. (Mathematics), Fudan University, Shanghai, P.R.China, 1984 M.Sc. (Mathematics), Fudan University, Shanghai, P.R.China, 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF D O C T O R OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES ( INSTITUTE OF APPLIED MATHEMATICS ) DEPARTMENT OF MATHEMATICS We accept this thesis as conforming to the required standard THE U&iyERSITY OF BRITISH COLUMBIA February 1995 © Xunlei Jiang, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of M a t h e m a t i c s The University of British Columbia Vancouver, Canada Date i n r i 1 7? • 1 QQ5 DE-6 (2/88) Abstract Theoretical and practical aspects of the design and implementation of the streamline-diffusion (SD) method for semiconductor device models are explored systematically. Em-phasis is placed on the hydrodynamic (HD) model, which is computationally more chal-lenging than the drift-diffusion (DD) model, but provides some important physical infor-mation missing in the DD model. We devise a non-symmetric SD method for device simulations. This numerical method is uniformly used for the HD model (including a proposed simplification (SHD)) and the DD model. An appropriate SD operator is derived for the general non-symmetric convection-diffusion system. Linear stability analysis shows that our proposed numerical method is stable if the system can be symmetrized. Stability arguments and numerical experiments also suggest that the combination of the method of lines and the semi-discrete SD method may not be appropriate for the transient problem, a fact which often has been ignored in the literature. An efficient method, consistent with the SD method used for conservation laws, is developed for the potential equation. The method produces a more accurate electric field than the conventional Galerkin method. Moreover, it solves for the potential and electric field in a decoupled manner. We apply our numerical method to the diode and MESFET devices. Shocks for the diode in one and two space dimensions and the electron depletion near the gate for the MESFET in two space dimensions are simulated. Model comparisons are implemented. We observe that the difference in solutions between the HD and DD models is significant. The solution discrepancy between the full HD and SHD models is almost negligible in MESFET simulation, as in many other engineering applications. However, an exceptional case is found in our experiments. ii T Table of C o n t e n t s Abs t r ac t i i List of Tables vii List of Figures viii N o m e n c l a t u r e xi Acknowledgement xiii 1 Introduct ion 1 1.1 Preliminary 1 1.1.1 Semiconductor Device Models 1 1.1.2 Numerical Tools 4 1.2 Brief History Review I: Numerical Simulations 6 1.3 Brief History Review II: Streamline-Diffusion Methods 9 1.4 Objectives and Scope of the Thesis 11 2 Semiconductor Dev ice Mode l s 14 2.1 Introduction 14 2.2 Mathematical Equations 15 2.2.1 PDEs for HD Models 15 2.2.2 PDEs for DD Models 17 2.3 Vector Form of HD Models 18 2.3.1 Form of the HD Model 19 2.3.2 Form of the SHD Model 20 iii 2.3.3 Eigenvalue-Eigenvector Problem 21 2.3.4 Properties 23 2.4 Boundary Conditions 24 2.4.1 Ohmic and Schottky Contacts 25 2.4.2 Traditional Boundary Conditions 26 2.4.3 "In-Flow" Boundary Condition 27 3 N o n - S y m m e t r i c Streamline-Diffusion M e t h o d 29 3.1 Introduction 29 3.2 Brief Review 31 3.2.1 Concept of the SD Method 31 3.2.2 SD Methods for Systems 33 3.3 Semi-Discrete SD Formulations 36 3.3.1 One-Dimensional Case 36 3.3.2 Multi-Dimensional Case 39 3.4 Full Discrete SD Formulations 40 3.4.1 Space-Time Finite Element Methods 41 3.4.2 Methods of Lines 45 3.5 Stability Arguments 45 3.5.1 Space-Time Finite Element Methods 45 3.5.2 Methods of Lines 48 3.6 Numerical Tests 49 3.6.1 Test 1 50 3.6.2 Test 2 57 4 Potent ia l Equat ion Solver 62 4.1 Introduction 62 iv 4.2 Previous Work 64 4.2.1 The Div-Rot Method 64 4.2.2 The Augmented Method 65 4.3 A New Method 66 4.4 Numerical Test 69 5 S imulat ion I: the p+-p-p+ D i o d e 72 5.1 Introduction 72 5.2 SD Method for Devices 73 5.2.1 Shock-Capturing SD Formulation 74 5.2.2 Concept of the Shock Capturing SD Method 76 5.3 Predictor/Multi-Corrector Algorithm 78 5.4 Physical Parameters 81 5.5 Numerical Simulations 83 5.5.1 One-Dimensional Simulations 83 5.5.2 Two-dimensional Simulations 90 5.5.3 Comments and Conclusions 92 6 Simulat ion II: the M E S F E T Dev ice 96 6.1 Introduction 96 6.2 Schottky Barrier 98 6.2.1 Formation of Schottky Barrier 98 6.2.2 Mathematical Model 100 6.3 Discretizations for DD Models 103 6.3.1 The SD Method for DD Models 105 6.3.2 A SD Method for the Steady State DDl Model 105 6.3.3 SG Method for the 1-D Steady State DDl Model 106 v 6.4 Numerical Simulations 107 6.4.1 Example 1 108 6.4.2 Example 2 I l l 6.5 Comments and Conclusions 112 7 Conclusion 125 7.1 Summary 125 7.2 Future Research Directions 128 Bibliography 129 Appendices 137 A A Property for the Hydrodynamic Model (2.1) 137 B Stability of the SD-Trapezoidal Scheme 139 C 1-D Test for a Set of Parameters 141 vi List of Tables 4.1 Errors of Eh with different parameter settings 70 4.2 Convergence rates for Eh with different parameter settings 71 5.1 Values of physical parameters 82 vii List of Figures 3.1 A space-time slab 41 3.2 U\ by (a), (b) and (c)-types of operators in test 1 using the constant-in-time space-time SD scheme 52 3.3 U2 by (a), (b) and (c)-types of operators in test 1 using the constant-in-time space-time SD scheme 53 3.4 Solutions by (c)-type operator with Cr = 1 in test 1 using the constant-in-time space-time SD scheme 54 3.5 Solutions by (c)-type operator in test 1 using SD-backward Euler scheme. 55 3.6 Solutions in test 1 using the Galerkin-backward Euler scheme 56 3.7 Solutions by (a) and (c)-types of operators in test 2 using the constant-in-t ime space-time SD scheme 60 3.8 Geometry of a 2-D p+-p-p+ diode in test 2 61 3.9 Electron density (left) and electron temperature (right) in test 2 using the constant-in-time space-time SD scheme 61 5.1 Geometry of p+ — p — p+ diode 74 5.2 HD model at 300K: Electron densities in 101 7cm - 3 ; velocities in I08cm/s; temperatures in Kelvin and electric fields in 104V/cm 85 5.3 HD model at 77K: evolutions of selected solutions 86 5.4 HD model at 77K: Electron density in 1017cra -3; velocity in 107cm/.s and Mach number; electron temperature in Kelvin; electric field in 104V/cm (from left to right and top to bottom) 87 viii 5.5 SHD model at 77K: Electron density in 101 7cm - 3 ; velocity in 108cm/s ; electron temperature in Kelvin; electric field in 104V/cm (from left to right and top to bottom) 88 5.6 Comparisons of J — VD and J — To curves using the HD and SHD models. 89 5.7 Doping profile in 1017 cm~3 91 5.8 HD model at 300K: Electrical field 94 5.9 HD model at 300K: Streamlines of the current 94 5.10 HD model at 300K: Electron density (left) in 1017 cm~3 and electron tem-perature (right) in Kelvin 95 5.11 HD model at 77K: Electron density (left) in 1017 cm~3 (image is inverted) and electron temperature (right) in Kelvin 95 6.1 Schematic of MESFET operation 97 6.2 Formation of Schottky barrier 98 6.3 Electric field near the Schottky contact 100 6.4 p and — ip with cf>B = 0.8V and po = 1017cra3. Schottky contact is at y=0.2.101 6.5 Electron density (left) and potential (right) with <f>g = 0.8V and pu = 1017cra3 in thermal equilibrium 104 6.6 Doping and 2-D geometry of the device for example 1 109 6.7 Doping and 2-D geometry of the device for example 2 112 6.8 Electron densities in 1017cm~3 for example 1 113 6.9 Electron currents for example 1 114 6.10 Mean velocity components for example 1 115 6.11 Electron temperature in Kelvin for example 1 116 6.12 Electric Potential for example 1 116 6.13 Solutions for (f>B - (f>o = 0.8V 117 ix 6.14 Comparison of electron densities at y = 0.175//m for example 1 118 6.15 Electron density and current with the finer mesh for example 1 119 6.16 Electron density and current for example 2: a grid of 96 x 32 points. . . 119 6.17 Electron density and current for example 2: a grid of 192 x 64 points. . . 120 6.18 Evolution of the electron density for example 2 121 6.19 Evolution of the electron current for example 2 122 6.20 Evolution of the electron temperature for example 2 123 6.21 Evolution of the potential for example 2 124 C.l HD model at 300K: using Fatemi et a/.'s parameters fi0 = 0.145, vs = 0.1 and K = ^ ^ 141 x Nomenc la ture Thermal conductivity Dielectric constant Electron mobility Electron density Surface electron density Doping profile Intrinsic electron density Relaxation times Barrier height function Working function Quasi-Fermi level Electron affinity Electro-static potential Built-in potential Device boundary Contact and insulating parts of T Device domain Time step Sound speeds Charge of an electron Spacing parameter Boltzmann constant Effective mass of an electron Time variable xi Mean velocity Saturation velocity Thermionic emission velocity Space variables Intraband collision terms Boundary operator Electron diffusivity Electric field Conduct and valence band edges Fermi level Electron current density Differential operator Effective density of state Electron temperature Lattice or room temperature Space-time Finite element subspaces Solution of PDEs Numerical solution of PDEs Semi-discrete Finite element subspaces Biased voltage at the drain Electron energy density Depletion width xii Acknowledgement I would like to express my deep gratitude to my supervisor Dr. Uri Ascher for his encouragement and guidance throughout the course of this work. I sincerely thank him for the fruitful suggestions and the patient proofreading so that this thesis could be completed. Thanks also go to the members of my supervisory committee: Drs. Robert Miura, Michael Ward and Brian Wetton. Their helpful comments during this research and writing of the thesis are greatly appreciated. xm Chapter 1 Introduction 1.1 Pre l iminary 1.1.1 Semiconductor Dev ice Mode l s Numerical simulation of semiconductor devices provides an effective tool for device design in many cases of practical relevance. This is a thriving research area for scientists and engineers. Reliable numerical simulators must apply robust discretization techniques to a sufficiently accurate device model. This model should reflect important physical features needed, while redundant information is discarded. Semiconductor device modeling started in the early fifties when the fundamental semiconductor device model, the so-called drift-diffusion (DD) model, was formulated by Van Roosbroeck [70]. This remarkable work started what today is an interdisciplinary subject that combines the efforts of physicists, electrical engineers and mathematicians. Much of the progress has been achieved in understanding the characteristics of devices and improving device design using perturbation analysis, numerical computation and solid state physics. In 1970, a sophisticated transport model was established for devices by Blotekaer [6], which is usually referred to as the hydrodynamic (HD) model. This model includes the non-negligible hot carrier effect for submicron devices. As compared with the conventional DD model, much less work has been done on the HD model until recently. The HD model consists of an incompletely parabolic system of conservation laws for 1 Chapter 1. Introduction 2 charge density1 p, momentum J and energy W, and Poisson's equation for the electro-static potential ip. It is a highly nonlinear convection-diffusion system. The conservation laws can be derived from the Boltzmann equation by taking the first three moments with an additional heat conductivity term. The parabolic continuity equation in the DD model, however, is derived by the first moment method (see [60]). The HD model is a cost-affordable extension of the conventional DD model. The recovery of the DD model from the HD model can be attained under assumptions of low carrier densities and small electric fields. Although the DD model constitutes a popular model for the semiconductor device simulation and is mathematically better understood than the HD model, the limitations of this model have been encountered more and more frequently with the miniaturization of devices or with the increment of the field strength. (Reduction of the active device size can be reflected by the increase of the field strength if one scales the equations of the model with suitable scaling parameters.) Alternatively, the HD model represents a reasonable model, which takes into account important physical phenomena such as velocity overshoots that are missing in the DD model [68, 57, 76]. A thorough investigation of carrier dynamics in semiconductor devices may be carried out using, for example, the Boltzmann equation for particle transport. This, however, involves a large amount of computational effort and makes application of numerical meth-ods a formidable task. Moreover, it contains a lot of redundant information. Shur [76] demonstrated that in order to describe electron dynamics, the full solution of the Boltz-mann equation is not necessary. For today's technology, the HD model (and its exten-sions and in many cases its simplifications) represents a reasonable compromise between the contradictory requirements of physical accuracy and computational tractability. It should be possible to develop semiconductor device analysis programs based on the HD 1In the semiconductor context, the notation n, not the p, is commonly used for the electron charge density. However we will use n for another purpose in the thesis. Chapter 1. Introduction 3 model with the same generality as those based on the simpler DD model which are now routinely used in research and design environments. Although advantages of the HD model over the DD model have been reported in the literature (cf. [76, 17]), weaknesses in the HD model also exist. For example, the commonly used model for the thermal conductivity KomfioTo K = e is an empirical one. The lack of a rigorous model for the heat conductivity causes serious problems like the artificial velocity overshoot [22, 51]. Even though some intuitive remedies were proposed (see [30]), the limitations of these models are still unknown. Therefore discrepancies between solutions of the HD model and the real physical behavior are inevitable in some circumstances. The DD model is by no means obsolete. It is valid under certain conditions and is still widely used either to simulate a variety of devices or to develop and test numerical simulation tools [11, 78]. The boundary of validity of the DD model is not entirely clear in applications. The established work on the DD model provides at least a gateway towards understanding the sophisticated HD model. As a compromise, some simplifications of the HD model (see [7, 8, 55]) were also employed and are considered good approximations to the full HD model in the engineering community. These simplifications alleviate some numerical difficulties inherited in the full HD model while still capturing the important physical phenomena missing in the DD model. However, it is questionable whether one can completely replace the full HD model by one of its simplifications and still obtain reasonably accurate physical data of interest, such as the I-V characteristics, for wide parameter ranges. The order of complexity of the addressed models is: the DD model —» the simplifi-cation of the HD model —» the full HD model. Their associated ranges of applicability Chapter 1. Introduction 4 have yet to be determined. Comparisons of the ability of these models to reproduce desirable properties are needed. It is especially interesting to investigate how "close" the simplification of the HD model approximates the full HD model, and in what situations their difference appears. 1.1.2 Numerica l Tools Numerical simulation of semiconductor devices is one of the challenging tasks in the area of scientific computation. Difficulties associated with solving the semiconductor device equations were discussed by many researchers (see e.g., [1, 3, 58, 60, 22, 29, 51]). The primary sources of difficulties are: • Some of the coefficient functions and unknowns appearing in the device equations vary by several orders of magnitude so that their variations cannot be accurately resolved on a mesh of practical size by classical discretizations. • These are singular perturbation problems with possible boundary and interior junc-tion layers. • Extremely large condition numbers are often encountered for the linearized discrete problem. Possible instability of nonlinear iterations may result [3]. • The conservation laws in the HD model possess hyperbolic modes. It is known that they can develop discontinuities or shocks (cf. Gardner et al. [29]). Numerical techniques in handling the transonic electron flow problems are more difficult to devise than for the subsonic case. • The time-dependent problems are stiff and give differential-algebraic systems of index 1 or 2 (see [1]). Chapter 1. Introduction 5 Efforts towards designing reliable numerical methods were made in the last two decades to overcome one or more of the difficulties addressed above. In particular, the successful Scharfetter-Gummel (SG) discretization [71] (in its original form based on the physi-cally motivated approach to the one-dimensional diode, or in the modified and extended form) has been used up to now for the DD model. Comparatively, very slow progress has been made in developing reliable numerical tools for the HD model. Some efficient methods, such as the so-called essential non-oscillatory (ENO) method, were applied to solve the HD model, but they are one-dimensional in nature. Though applications to multi-dimensional problems are possible (see [51]), they lack theoretical backing. At-tempts at designing reliable numerical methods for the multi-dimensional HD models were even less successful than for the one-dimensional ones. Previous work, including at tempts to extend the SG discretization, either was not robust enough to resolve shock discontinuities or poorly smeared the shock fronts [69, 78]. The shock-capturing capability is an important concern in this thesis. Even without shocks, other sharp layers, such as junction layers, still cause difficulty. The classical methods, such as the central difference and Galerkin methods, are not good candidates. They often generate spurious oscillations near sharp layers. Upwinding schemes (which are designed according to the influence of the flow direction and strength) are therefore required, not only for the HD model but for the DD model as well (see [78, 11]). On the other hand, the requirement for higher resolution of sharp profiles demands a special kind of upwinding method. Simple minded upwinding methods, such as the artificial viscosity method, are not suitable. The streamline-diffusion (SD) finite element method, introduced by Hughes and Brooks [34], provides a tool to achieve the goal. It is also called the streamline/upwinding Petrov-Galerkin (SUPG) method. The essence of the SD method is a Petrov modifica-tion of the test functions in the standard finite element subspace along the streamline Chapter 1. Introduction 6 direction. An upwinding effect is then introduced, therefore spurious oscillations can be suppressed in the presence of discontinuities or sharp changes of solutions. The advan-tages of the SD method are: • The algorithm does not need particular adjustments in the case when shocks are present as compared with other shock-capturing methods, so it is a general purpose method suitable for both hyperbolic and parabolic problems. A suitable amount of artificial diffusion (mainly along the streamline direction) is "turned on" when the problem considered is convection-dominated, and is "turned off" when the problem is diffusion-dominated. • Local dissipation (or smearing) does not strongly degrade the sharpness of shocks or interior layers as compared with conventional upwinding schemes. • It is especially suitable for handling multi-dimensional cases. The layers need not align with grid lines. • It has theoretical backing (see e.g., [51, 46, 50] and a series of papers by Hughes et d. [38, 42, 40, 41, 36, 39, 37, 35, 75, 77]). • It is a reasonable extension of the SG discretization (see [78]). Mathematically, both of them are methods of exponential fitting. The mechanism makes it possible for us to solve both the HD and DD models in the same framework, whereas in previous work, researchers more often than not tended to treat them separately. 1.2 Brief History R e v i e w I: Numerica l Simulations The emphasis towards numerical simulation techniques for semiconductor devices began in the early sixties. The breakthrough was made by Scharfetter and Gummel [71] in 1970. Chapter 1. Introduction 7 They proposed a nonstandard discretization based on a physical approach to the one-dimensional silicon diode. In spite of its success, the mechanism was not well understood mathematically until 1980 when Doolan et al. [20] explained it. The mathematical formulation of the SG discretization was termed exponential fitting. It has since been generalized to two- and three-dimensional DD problems. Bank et al. [4] extended the SG method to two-dimensional problems by using the control volume method on rectangular meshes and by using the finite element method on triangular meshes. Polak et al. [65] further generalized the SG method to allow the use of quadrilateral meshes. Similar types of SG versions were analyzed on tetrahedra by Bank el al [5] and by Sever [74]. However, the multi-dimensional SG extensions were not as successful as the original one-dimensional one. Brezzi et al. [11] introduced the exponential fitting idea to the two-dimensional DD model using mixed and hybrid finite element methods. Nice properties of the methods and satisfactory numerical results were presented. However, the mixed or hybrid finite element subspaces should satisfy the so-called Brezzi-Babuska condition (see [10]). Different settings of nodal points (staggered grids) have to be used for different solution components. This usually brings some inconvenience. Moreover, their work was done just for the continuity equation, not for the whole system. For the time-dependent problem, an implicit t ime discretization scheme is recommended. Based on this idea, various time evolutionary schemes for the DD model in two and three dimensions were evaluated by Wu [81]. The linkage of the SD method to the SG method is attributed to the work of Sharma and Carey [78]. They showed that the original form of the SG method coincides with the SD method for the one-dimensional DD model. Therefore, the SD method provides a reasonable extension of the SG method for multi-dimensional problems. Recent literature pertaining to the numerical solution of the HD models includes the studies of Rudan and Odeh [69] and Odeh et al. [62]. The SG discretization was Chapter 1. Introduction 8 extended to the case of HD models. But these efforts, to some extent, were confined to only cover the case of subsonic transport, i.e., solutions were not expected to experience shock discontinuities. The hyperbolic nature of the HD model may cause shocks for the transonic flow. The investigation of this phenomenon has become an interesting topic and has attracted much attention recently. For model studies of the one-dimensional p+ — p — p+ ballistic diode, mathematical analysis by means of viscosity approximations was considered by Gamba [24, 25]. She claimed there was the possibility of boundary layer formation. Ascher et al. [2] presented an analysis of transonic shocks for a current driven problem in the isentropic case. A phase plane analysis gave a full explanation of how and when the shocks occur. A similar analysis was carried out later by Markowich and Pietra [59] in the non-isentropic case. In both cases, numerical computations illustrating the structure of the solutions were also reported. In the simulations of the p+ — p — p+ diode by the full HD model, shocks were found by Gardner [26] and Fatemi et al. [22] in the one-dimensional steady state and non-steady state cases, respectively. A two-dimensional simulation with shocks was done by Jiang [45] using the shock-capturing SD method. Simulations for more realistic devices beyond the simple p+ — p — p+ model comprise the early work by Cook and Frey [17] for the MESFET, and by Kreskovsky et al. [54] for the PBT. The recent work by Jerome and Shu [51], using the operator splitting ENO method, and Chen et al. [13], using the discontinuous Galerkin method for the two-dimensional MESFET, are more method oriented. Although there were no shocks observed, the shock-capturing algorithms proved to be useful when sharp layers were encountered. In many situations, simplifications of the full HD model proved to be easy to use. One of them is to neglect the convective term in the momentum equation of the conservation Chapter 1. Introduction 9 laws. This was used by Lin et al. [55] to express the HD model in self-adjoint form through a set of new Slotboom-like variables like those that are commonly used for the DD model. The so-called energy balance (EB) model is almost exclusively used in the engineering field in place of the full HD model (see the papers by Rahmat et al. [67] and Benvennuti et al. [7, 8]). In the full HD model, the energy W is the sum of kinetic energy and thermal energy, whereas in the EB model, W is assumed to be the thermal energy alone. Although the difference between solutions of the full HD and EB models is small in many cases, exceptional situation do happen in which the kinetic energy is no longer small. 1.3 Brief His tory R e v i e w II: Streamline-Diffusion M e t h o d s The SD method was introduced for the stationary convection-dominated convection-diffusion equation by Hughes and Brooks in 1979 [34]. They proposed a Petrov modifi-cation of the classical Galerkin finite element method, which is also known as the SUPG method. This new method introduces artificial diffusion only along the streamline direc-tion to suppress possible oscillations near discontinuities or sharp layers of the solution without causing overly large crosswind diffusion. Mathematical studies of the method were started by Johnson and Navert [46], and continued by Navert [61] and Johnson et al. [47, 50]. The work includes the extension to the time-dependent problem with the discontinuous Galerkin methodology. Primary stability, convergence and crosswind smear evaluation results were established for the model problem. The SD method gives added stability as compared with the classical Galerkin method. However, the L2—norm error estimation for the former is a half-order lower than for the latter. That the error bound can not be improved was shown by Peterson [64], who presented an example using the discontinuous Galerkin method, a counterpart of the SD method. Chapter 1. Introduction 10 From numerical experiments, a proper amount of crosswind dissipation is needed to entirely control the oscillations near the discontinuities. A further Petrov modification of the SD method, a shock-capturing (SC) modification, was proposed by Hughes et al. [42]. The method is a true nonlinear one even if the problem considered is linear. It was adopted by Johnson and Szepessy [51] for solving the nonlinear Burgers equation. The explanation of such a modification was given by Johnson et al. [49] later. There were many efforts along this line afterwards, such as those of Codina [15] and Carmo and Galeao [12]. Now the shock-capturing SD method has become a standard formulation. Systematic research was done by Hughes and his coworkers in the sequence of ten papers entitled A new finite element formulation for computational fluid dynamics. The cardinal forms of the SD operator (and its variant, the least-squares operator) and the SC operators were well established for the convection-diffusion scalar equation as well as for the case of a symmetric system. Thorough analyses of the stability, convergence and nonlinear iterative technique were carried out. Their work was directly applied to the compressible Euler and Navier-Stokes equations with success. Their contribution therefore laid a foundation for standardizing the SD method. The advantage of the SD method utilizing the symmetric form for a PDE system was addressed in the above mentioned papers. For instance, there is guaranteed entropy stability of the symmetric SD scheme for the compressible Navier-Stokes equations [77]. However, the method depends strongly on symmetrization (that may not be known sometimes) of each individual problem. Even when a symmetrization procedure is known, the use of primitive variables to avoid added complexity is often preferred in applications. In [9] (by Le Beau et al.), SD methods based on both symmetric and non-symmetric forms of Euler equations were investigated and compared intensively for subsonic, transonic and supersonic flow problems. The experiment indicates that the advantage in using the symmetric formulation is invisible computationally. Chapter 1. Introduction 11 It seems that analysis for the non-symmetric case may not be carried out by a straight-forward extension of the symmetric one. Some of the early treatments were inappropriate (see [43]). Awareness of the importance of a choice of the SD operator was reflected in the recent literature comprising the work by Zhou [83] and Le Beau et al [9] for the Euler equations, and by Jiang [45] for the HD model of the semiconductor device equations. In [45], which forms part of this thesis, a proper operator scaling is also presented, which was often ignored and yet is important. The importance of the scaling was stressed by Hughes and Mallet [40]. 1.4 Object ives and Scope of the Thesis The aim of the thesis is to investigate efficient, robust and consistent numerical methods for solving semiconductor device problems as formulated via the HD and DD models, and to explore the physical behaviors of the models through device simulations in applications. In Chapter 2, we present semiconductor device models which we are going to consider in an decreasing order of complexity: the HD model —> a proposed simplified HD (SHD) model —• the mixed formulation DD (DD2) model —> the conventional DD (DD1) model. Some properties of the device models useful for formulation of the non-symmetric SD method are investigated. In Chapter 3, the non-symmetric SD method for convection-diffusion systems is de-signed to avoid the extra complexity of symmetric formulations of the systems involved. We derive a non-trivial SD operator for the non-symmetric system, which extends the previous work in [45]. The design includes two concerns: one is to follow the correct directions of streamlines and another is to provide suitable scaling. Both of them are important as is shown by concrete examples. The linear stability analysis presented shows that the proposed numerical method is stable if the system can be symmetrized. Chapter 1. Introduction 12 This analysis, together with the support of the numerical evidence that we give, suggests that the combination of the method of lines and the SD semi-discretization may not be appropriate for time-dependent problems. In Chapter 4, we propose a consistent Poisson solver which provides a higher-order convergence rate for the electric field E = —Vifi which is the only quantity needed in the conservation laws. The consistency of this Poisson solver is composed of two parts: (1) the order of accuracy of the electric field E is the same as that of solutions for the system of conservation laws when employing finite element subspaces with the same degree of polynomials; (2) it uses the same grid points as used by the SD scheme for the system, and therefore the device simulation algorithm can be more easily coded. Not much attention was previously paid regarding this consistency. This new method also provides a way to solve ip and E in a decoupled nature as compared with its other counterparts. Error estimation and numerical tests are presented. In Chapter 5, the p+ — p — p+ diode is simulated by our numerical method. The effect of the "in-flow" boundary condition is incorporated weakly in the SD formulation for the first t ime in device simulations. This enhances the numerical stability. Shocks are captured at 77K for transonic flows in one and two space dimensions. In these cases, we find from our numerical evidence that the SHD model no longer approximates the full HD model very well as compared with many other cases encountered in applications. In Chapter 6, we consider the simulation of 2-D MESFET (metal semiconductor field effect transistor) devices. Electron depletion near the gate is simulated. The leakage current of electrons is correctly reflected, although it is very small and can be neglected. The method successfully handles the case where the abrupt junctions are not aligned with grid lines. Experiments show that the performance of the method is very satisfactory and the computational results correctly reflect physical expectations. Model comparisons are implemented. The difference in solutions between the HD and Chapter 1. Introduction 13 DD models is significant in simulation of the MESFET device. However, the difference in solutions between the full HD and SHD models is negligible in this simulation. Finally, in Chapter 7, we conclude by summarizing the work that has been done in the thesis and propose future research directions. Chapter 2 Semiconductor Dev ice Mode l s 2.1 Introduct ion Both the hydrodynamic (HD) and drift-diffusion (DD) models have been discussed and used in the literature of semiconductor device simulation. The HD model, however, has received additional attention in recent years with the development of VLSI technology. It is known that the HD model takes into account important physical phenomena missing in the DD model. Though the limitations of the DD model are becoming more and more obvious, studies of this model are by no means useless. Not only is the simpler DD model sufficient in many simulation cases, also many ideas which originated in the study of the DD model are often readily extendable to the study of the HD model (e.g., see the extension of the SG method by Rudan and Odeh [69] and the new Slotboom-like state variable by Lin et al. [55]). In this chapter, mathematical formulations of the HD and DD models are presented. The DD model can be attained as a simplification from the HD model under certain assumptions. We consider, in Section 2.2, a hierarchical deduction from the full HD model to the simplified HD (SHD) model, to the mixed-formulation DD (DD2) model and finally to the conventional DD (DDl) model. Some properties of the formulated systems of PDEs, which are useful in the later chapters, are presented in Section 2.3. The issue of boundary conditions is discussed in Section 2.4. All the models that we consider are for electron carrier transportation only, i.e., 14 Chapter 2. Semiconductor Device Models 15 ignoring the hole carrier transportation. These models can be used to simulate devices such as the p+ — p — p+ diode and the MESFET. For other devices, where the electron and hole carriers are equally important, it is possible to modify the models presented here for the electron by adding counterparts for the hole (see [55] and [54]). From the numerical point of view, we can decouple the effects of the electron and hole iteratively (but with care as discussed in [3]). Therefore, with such a procedure we essentially solve a PDE system for one carrier at each iteration. 2.2 Mathemat ica l Equat ions In this section, we review mathematical formulations for d-dimensional ( d = l , 2) semi-conductor device problems via the HD and DD models. We assume that the formulations have been properly scaled. 2.2.1 P D E s for H D Mode l s Let us consider first the HD model (or full HD model) of semiconductor device equations which can be derived from taking the first three moments of the Boltzmann equation. It expresses the conservation of charge, current and energy (see [6, 60]) H D Model: pt + V • J = 0, (2.1a) Jt + vV • J + J • Vv + V{pT) = -&M. + Cj, (2.1b) Wt + V • {{W + pT)v) = -*&• + Cw + V • {npVT). (2.1c) Here (p, J, W)* is an unknown vector which represents the electron density, "(negative) current density" (which is different from the actual current density by multiplying a negative charge, but we still refer to it as a current density) and total energy density, Chapter 2. Semiconductor Device Models 16 respectively. The conservation laws (2.1) are coupled with Poisson's equation for the electrostatic potential if) by V • (AVt/0 =p-pD. (2.2) In the above equations, E — — V ^ stands for the electric field, v is the mean velocity and T is the scaled temperature (the usual temperature in Kelvin T <— fc&T/m, with kf, the Boltzmann constant). The physical parameters are: e the charge of the electron, m the effective mass of the electron, K the thermal conductivity and A the dielectric constant. The terms Cj and Cw, represent the rate of change of J and W respectively, due to intraband collisions. The forms of Cj and Cw are (see [26]) c, - -L, Tj W - \pT0 Uw = • Tw In the above, TJ and TW are relaxation times which are obtained empirically by mfj,0T0 TJ = eT 3mn0T0T TJ_ TW ~ 2evl{T0 + T) + 1l' where T0 is the lattice temperature, p,o is the electron mobility and vs is the saturation velocity. The thermal conductivity is taken as Kom/jioTo K = . e The so called doping profile is given by po which is a function of only the space variables. Finally, to close the above system, we need the following current-velocity and energy-temperature relations J = pv, W = \p\v\2+\PT, Chapter 2. Semiconductor Device Models 17 which are ensured by the parabolic energy band structure. The full HD model has highly nonlinear hyperbolic modes due to the presence of the convection term vV- J+J- Vi> in (2.1b). This causes difficulties in the numerical treatment of this model. To remove such difficulties, a reduced HD (RHD) model (neglecting the convection term in (2.1b)) is often used [55]. It is considered to be a good approximation of the full HD model in applications. However, we propose a slightly different, still simplified (SHD) model as follows. S H D Model : pt + V • J = 0, (2.3a) Ji + V(PT) = -e-£ + Cj, (2.3b) Wt + V • ({W + pT)v) + A* = -s^- + Cw + V- (KPVT). (2.3c) Here the correction term A* is given by A* = - - . ( t , V - J+J-Vv). P The rationale for adding the correction term A* will be given later in Section 2.3. Note that shocks are supported by the HD model, but not by the SHD model. For the HD model, possible discontinuities caused by nonlinear shock waves (if there are any) are in the velocity v and electron density p [29]. 2.2.2 P D E s for D D M o d e l s The DD models can be derived by taking the zeroth-order moment of the Boltzmann equation. But we can obtain them by simply using the SHD model. If we further assume that T = T0, (2.3) yields the following DD model with a system of two equations (DD2). D D 2 Model : pt + V • J = 0, (2.4a) Chapter 2. Semiconductor Device Models 18 Jt + T0Vp=-e-^ + Cj. (2.4b) This is known to be a suitable approximation for the cases of small electric field or near thermal equilibrium [60]. If we define D0 = mp,0To/e (Einstein relation), the diffusivity, then (2.4b) can be written as ^ Jt + nopE + D0VP = -J. e Ignoring the term Si^s-Jt (which is assumed to be small for ultra-small devices [60]), we are able to express the current J explicitly as J = -(iiopE + DQVp). If we plug the above J into (2.4a), we obtain the conventional DD (DD1) model. D D 1 Model : pt - V • (fiopE + DQVp) = 0. (2.5) The DD2 model can be viewed as a mixed formulation of the DD1 model, but (2.4) is not used as frequently as (2.5). Here we use it to test the capacity of our designed streamline-diffusion scheme for systems, which is generalized from the optimal formulation of the same scheme for scalar equations. This can be done by comparing numerical solutions of the two models in steady state where the analytical solutions of the two models are identical but numerical solutions are no longer the same. 2.3 Vector Form of H D Mode l s To ease our numerical formulation, we change (/?, J, W)1, the vector of primitive variables, to (/?, J, X1)*, a vector of working variables. System (2.1) or (2.3) (assuming d = 2 for the Chapter 2. Semiconductor Device Models moment to simplify our notation) is then transformed to AoUt + A-VU + C(U) = V • (tCVU), where, U = (/?, J, T)f. AQ, A, C and /C are defined below. 2.3.1 Form of the H D Mode l In equation (2.6) for the HD model, the convection matrix is A = with Ax = 1 0 0 T-J± 2 i 0 P2 J i T P P 0 p p T 0 fJx , A2 2 J 1 / 0 J} J?, P2 T2 p* p 0 A p 0 0 1 A p 2J 2 p T The mass matrix A0 = diag(l, 1,1,-p) is symmetric positive definite. The force term is / C = 0 Tj pEi + Z \j? . IP(T-T0)+{; \ prj ' TW and the diffusion matrix K = diag(K, K) Chapter 2. Semiconductor Device Models 20 with the symmetric positive semi-definite matrix K = diag(0,0,0, Kp), (2.12) Since A\ and A2 are not symmetric, we will say that the system (2.6) is non-symmetric. 2.3.2 Form of the S H D Model The matrices A\ and A2 in equation (2.6) corresponding to the SHD model are / Ai = V o T 0 .ML p 1 0 0 0 0 0 T 0 0 \ p 0 3 T { , A2 = V 0 0 T J,.T P 0 0 0 0 1 0 0 T 0 0 p l J 2 J (2.13) We have added the correction term A* in (2.3c). To see why such a term is added, we look into the details of the transformation from (2.1) to (2.6) by changing the variables. It is obvious from the energy-temperature relation that z pz p i We derive (2.8) by replacing pt using (2.3a) and Jt using (2.1b). Note that only the third equation (2.1c) is changed. If we then neglect the convection term in equation (2.1b) of the full HD model in this stage, we obtain (2.13) which corresponds to what we called the SHD model using the variables (p, J, T)*. Changing (/?, J, T) f back to (p, J, W)\ this model turns into (2.3) with the appearance of an extra term A*. The structure of the matrices in (2.13) is simpler than that of the corresponding matrices for the R.HD model. For example, the last row of A\ corresponding to the RHD model is eUx = ( - J i T | J |2 J i _ \J\2 + J? hh 3 •)6 p 2 /r 2 Chapter 2. Semiconductor Device Models 21 while other rows are the same as those of A\ in (2.13). There is another simplification of the full HD model, the so-called energy balance model (EB), which is used widely by engineers [7, 8]. In this model, the kinetic energy is assumed to be negligible. Therefore the total energy is equal to the thermal energy, i.e., It is very simple to express the last row of A\ corresponding to the EB model as e'Ax = ( 0 , T , 0 , ^ J i ) , which is simpler than that of the SHD model. However, the eigenvalue-eigenvector prob-lems for the EB model as well as for the RHD model are very complicated as compared with those for the full HD and SHD models. 2 .3.3 E i g e n v a l u e - E i g e n v e c t o r P r o b l e m Let us define where fi = (^1,^2)* is an arbitrary parameter vector satisfying | JJL | 2 = {H\ + fJ-l)1^2 = 1-Consider the generalized eigenvalue-eigenvector problem (AM - XpAo^i = 0 (2.14) (i = 1 , . . . , 4), then the generalized eigenvalues of AM are: Type HD SHD Eigenvalues v • fX v • fi v • fi — c v • JJL + c 0 v • (J, —c c Chapter 2. Semiconductor Device Models 22 In the above, c = A / 5 T / 3 , the "sound speed" for K = 0. With the order of eigenvalues as above, A^ has the following generalized eigen-matrix Y„. \ 0 "/*2 Ml 0 1 t>l V2 3c2 5p 1 1 V\ — C//i Uj + C[J,\ v2 - CfX2 v2 + C/J,2 2c2. 2j£ 5p 5p (2.15) for the full HD model, and Y„. \ 0 -M2 Ml 0 1 filV • fi H2v • jU M 2 - T 1 -CMi - C ^ 2 2c2 5p 1 C ^ i C/^2 2c2 5p (2.16) / for the SHD model. The matrices Y^ are formed by assembling eigenvectors Y^i (i = 1 , . . . ,4) as columns. The inverse matrices of Y^ for the HD and SHD models are / Y~1 = V X fJ, —\X2 I 0 Mi 0 0 \ _ - £ -2c r 10 .a, 2c H2. _P_ 2c 2c 2 t'-M i _3_ W. 2c T 10 2c where x is the cross product operator, and / Y. - l 0 2T 3 ( C 2 - ( U - M ) 2 ) 5f-M+3C 1 0 ( c + w ^ ) 3c—5t;-^i \ 10(c—v-fj.) M2 0 —2c/«i 2cfii H2. 2c Ml 0 -2cfx2 2cfx2 2c2 ) (2.17) (^•M)2 2c(c+v-fi) (2.18) 2c(c—f-^i) / respectively. Therefore, we can diagonalize the matrix A^A^ to AM with eigenvalues of (2.14) forming its entries. Chapter 2. Semiconductor Device Models 23 2.3.4 Propert i e s We have the following results which enable construction of a suitable non-symmetric SD formulation for the HD models in the next chapter. Propos i t ion 2.1 For the HD model, there exists a diagonal matrix Ap=diag([i\,... ^d\a) with constants /?,• > 0 (i = 1 , . . . ,d + 2) , such that ApY~x AQ1 KY^ is symmetric positive semi-definite. Proof. Here we only prove the typical case where d = 2 without loss of generality. From the definitions of Y^ A0 and K we have ' 0 ^ Y-'A^KY, = 9K To (0,-1, §,f). (2.19) V 2 / which is independent of p. It is easy to see that if we select 4 4 A^ = diag(l,l,-,-), (2.20) then the matrix ApY^ 1A01KY^l is symmetric positive semi-definite with rank 1. • Propos i t ion 2.2 For the SHD model, if \v\ < y/T, then there exists a diagonal matrix Ap = diag(/3i,..., Pd+2) with # > 0 (i = l,...,d + 2), such that ApY~1A01KYli i symmetric positive semi-definite. Proof. Similar to the argument in the proof of Proposition 2.1 (d = 2), / n \ IS Y-'A^KY, 2K 0 1 (y-ii)2— c2 1 2c(c+v-ti) 1 , \ 2c(c-v-ti) } 2c1 2r2 (2.21) Chapter 2. Semiconductor Device Models 24 We select A, , Jiag(l,'" • ">' - C\ «£±llA, SfcHjfl), (3.88) (v • p,y — 1 c c then the matrix K^Y^AQ KY^ is symmetric positive semi-definite. Since |u| < y/T, all fa (i = 1 , . . . , 4) are positive. • For system (2.1) with primitive variables, a similar result to that of Proposition 2.1 can be obtained. This can be found in Appendix A. 2.4 Boundary Condit ions The semiconductor device equations in Section 2.2 are posed in a bounded domain ft representing the device geometry. The boundary T of 0 is piecewise smooth for the two-dimensional problem and trivially consists of two points for the one-dimensional problem. The boundary can generally be split into two disjoint parts: r = r1Ur2. For devices under consideration, such as the p+ — p — p+ and the MESFET, Ti and T2 represent those parts of the boundary which correspond to metal Ohmic contacts and insulating parts, respectively. In the one-dimensional case, it reduces to T = Ti. For the MESFET device, Tj can be further split into two disjoint parts: r1 = rf|jrf, where T° denotes the parts of Ti corresponding to Ohmic contacts and Tf corresponding to Schottky contacts. Chapter 2. Semiconductor Device Models 25 2.4.1 Ohmic and Schottky Contacts The concept of Ohmic and Schottky contacts are very important in modern device tech-nology. The investigation and description of these contacts can be found in many pub-lished works (see, e.g., [72, 58, 66, 79] and the references therein). A Schottky contact is a special metal-semiconductor contact which produces a so-called Schottky barrier with a significantly large barrier height e<f>B (the definition of <f>B will be given in detail in Chapter 6). It influences the device by forming a so-called depletion region, where the current is negligible. An Ohmic contact possesses a negligible resistance. In principle, there is a potential barrier at the contact, with a very low height. The usual formation of the Ohmic contact is to dope the semiconductor heavily, so that it does not significantly perturb the device performance. For the pure voltage driven contacts, the electrostatic potential can be expressed by (see [72]) V>|rf = i>bi-i>so + VG, (2.23) xp\r? = rpu + VD. (2.24) tpbi is a built-in potential, which is given by ^ = ^ l n ^ , (2.25) e pi where pi is the intrinsic electron density, ^so is the potential related to the barrier height and VG, VD are external biased voltages at Tf and Tf, respectively. The typical value of •050 is usually between half a volt and one volt. The thermionic emission and diffusion theory [79] leads to the current expression at the Schottky contact, J = -vth(p ~ Po), (2.26) Chapter 2. Semiconductor Device Models 26 where Vth is the thermionic emission velocity and po is the surface electron density in thermal equilibrium (VQ = VD = 0), which admits the following expression: po = pD e x p ( - ^ 7 ^ 5 o ) . (2.27) When the Schottky contact operates in the reverse biased mode, the leakage current is small, i.e., J « 0, (2.28) at r f . Therefore we can simply use the following boundary condition for p at the Schottky contact: p\rf = po, (2.29) if the contact is in the reverse biased mode [72]. 2.4.2 Traditional Boundary Condit ions We now write the traditional boundary conditions as follows, according to [69], T | r i = To (2.30) and J-n\r2 = 0, (2.31) V T - n | r 2 = 0, (2.32) E-n\r2 = 0. (2.33) In the above, n denotes the unit normal vector on T. Besides these, we specify At the Ohmic contact: p\v? = PD, (2-34) </,|ro = tpu + Vo. (2.35) Chapter 2. Semiconductor Device Models 27 At the Schot tky contact: Hrf = Po, (2.36) V>|ro = *l>bi-il>so + VG. (2.37) We later take VG = 0. ipso will be specified in Chapter 6. Remark 2.1 The boundary condition (2.31) implies (see [69]) that Vp-n\r2 = 0 . (2.38) Conditions (2.31) and (2.38) are equivalent only in the DD models. • The boundary conditions (2.30)-(2.37) are physical boundary conditions. They turn out to be sufficient only for the SHD and DD models [80]. Mathematically, we need an extra "in-flow" boundary condition in case the two-dimensional HD model is used. This issue will be discussed in the next subsection. 2.4.3 "In-Flow" Boundary Condit ion System (2.1) is an incompletely parabolic problem. General well-posed boundary condi-tions for such a system were analyzed in [31] in the fluid dynamics context. In the context of the semiconductor device problem, Thomann and Odeh [80] carried out a discussion on the well-posedness of the boundary conditions for the two-dimensional full HD model. They found that an extra "in-flow" boundary condition according to the eigenvalue v • n of An (taking (i = n at the boundary) is needed in the case that the contact Ti is subsonic, i.e., hllrx^ln, (2-39) where c* = Vf (2.40) Chapter 2. Semiconductor Device Models 28 stands for the "sound speed" for re ^ 0. In order to describe the "in-flow" boundary condition, we define I \ n , a part of T\, as an "in-flow" boundary, i.e., it satisfies J-n\rtn<0. Then the "out-flow" boundary is Tout = Ti\Fin. According to [80] and [73], instead of (2.34) and (2.36), we can take J-r\rin = 0, p\rmt = PD, (2.41) where r is a unit vector tangential to V, and the "in-flow" boundary condition J-n\rin = f(p) (2.42) for a proper function / , derived from the so-called "maximal dissipative" principle. Two observations were made in [80]: 1. The "in-flow" condition (2.42) is satisfied by (2.26) at the Schottky contact with a defined / . 2. When the device size tends to infinity (i.e., the HD model reduces to the DD model), (2.41) and (2.42) lead to the normal Dirichlet boundary condition for p. Unfortunately, since no explicit function / is suggested in general, this approach is im-practical for implementation. Moreover, there is still some doubt about the necessity of this condition in many real case simulations, especially when the convection term in (2.1b) does not play a significant role. So far there is no general agreed-upon "in-flow" boundary condition available in practice. A problem dependent "in-flow" boundary con-dition will be proposed in Chapter 5. Chapter 3 N o n - S y m m e t r i c Streamline-Diffusion M e t h o d 3.1 Introduct ion Let us consider the convection-diffusion system of I partial differential equations in d spatial variables, CU = A0Ut + A-VU-V -}CVU-F = 0 (3.1) where A0 (mass matrix) is an / x / symmetric positive definite matrix, and Ul = (Ui,Ui,..., Ui), (solution vector) A* = (Ai, A 2 , . . . , Ad), (convection matrix) JC = diag(K, • • •, K) (diffusion matrix) with I x I matrices Aj and K. The term F in (3.1) is a source term. System (3.1) is parabolic if K is nonsingular and hyperbolic when K = 0. We are interested here in a more general case where K is allowed to be a singular nonzero matrix. This gives a so-called incompletely parabolic system. Definit ion 3.1 If Aj (j = 1 , . . . ,d) and K are all symmetric, then system (3.1) is said to be symmetric, otherwise it is said to be non-symmetric. • For some physical systems, such as the Euler equations of gas dynamics, the compress-ible Navier-Stokes equations and the hydrodynamic model of semiconductors, a sym-metrization procedure exists such that a symmetric form of (3.1) can be obtained (e.g., 29 Chapter 3. Non-Symmetric Streamline-Diffusion Method 30 [33, 38, 16, 73]). However, even in such cases where a symmetrizing transformation is known, the advantage of using the symmetric form in numerical computation doesn't seem to be obvious. In [9], streamline-diffusion (SD) methods based on both symmetric and non-symmetric forms of Euler equations were investigated and compared intensively for subsonic, transonic and supersonic flow problems. It was shown that solutions ob-tained using the two formulations are very close. In practice, people often prefer not to use the symmetric form due to the added complexity of the resulting formulation. This is the case, for example, in semiconductor device simulations, where the symmetriza-tion is seldom used. For these reasons, it is useful to consider the SD formulation for non-symmetric systems. The SD formulation for the symmetric convection-diffusion system was discussed by Hughes and Mallet [40] for the semi-discrete finite element scheme. An even more gen-eralized design was given by Shakib et al. [77] for the full discrete scheme. The study of one-dimensional model problems is crucial in order to obtain proper formulations since some previous studies failed to treat the problem appropriately even in such a case. It seems that the analysis for the non-symmetric case may not be carried out by a straightforward extension of the argument in [40]. In this chapter, a brief review of the SD method and some previous relevant results are presented in Section 3.2. Semi-discrete SD formulations for non-symmetric convection-diffusion systems of model problems are discussed in Section 3.3. A proper SD operator is derived in such a manner that it gives an optimal behavior if the system can be decoupled. In Section 3.4, the result in Section 3.3 is further generalized to full discretization schemes, including space-time SD formulations and methods of lines. Linear stability of commonly employed transient schemes is discussed in Section 3.5. In Section 3.6, numerical tests are presented, which support our arguments. From now on, we assume throughout this chapter that (3.1) is a linear constant coefficient system unless otherwise specified. Chapter 3. Non-Symmetric Streamline-Diffusion Method 31 3.2 Brief R e v i e w A brief review of SD methods relevant in our context is given in this section. For conve-nience, only problem (3.1) in steady state is considered. 3.2.1 Concept of the SD M e t h o d To introduce the concept, we consider the scalar convection-diffusion equation of (3.1) in the form a-Wu = KAU + / , (3.2) which is defined in 0 with ft C Md. We prescribe mixed boundary condition on the boundary T of the domain 0 «|rx = 9, (3-3) V u - n | r 2 = 0. (3.4) Let Th = {Oe} be a finite element triangulation of Q, where Q,e, (e = l , . . . , n e ; ) , are typically triangles or quadrilaterals. We define the following finite element subspace, with spacing parameter h, Vhg = {vh | vh € C°(n) , vh | f i e e Pk(ne), vh\Tl = g}, (3.5) where Pk is a polynomial space of degree at most k. The essence of the SD method is a Petrov modification of the usual test function vh G VQ in the classical Galerkin method by adding a perturbation term along the streamline direction a, i.e., Vh * Vh + TCt- Vvh, Chapter 3. Non-Symmetric Streamline-Diffusion Method 32 where r > 0 is a scaling parameter. To denote the above perturbation term, we define the SD operator Ph by Let (f,g)ns denote the usual inner product of functions / and g, and | | / | | n s be the corresponding Z2-norm on a space domain O5. We customize (/ , g) = (f,g)n and | | / | | = | | / | | n if tts = O. Then the SD variational form for (3.2) reads: Find uh G Vg such that for all Lp e V# (a • V u \ (p) + {nVuh, Vy>) (Galerkin) + J2(a -Vuh-V- (KVuh) - / , Ph{v)W (SD) = (/,¥>)• (3-6) Note that the second line of (3.6) yields J > • Vuh + • • •, Ph{<p))a. = E ( a • V«fc, ra • V^)n e + • • •. (3.7) e = l e = l streamline-diffusion The first term of (3.7) introduces an upwinding effect so that it enhances the stability of the scheme. The choice of r is important: either too small or too large a value may cause problems of oscillatory or overly diffused behavior of numerical solutions, respectively. We define the element Peclet number ae = -xlahK/K, (3.8) where \a\2 = (Yli=i aX)X^ a n < i he is the element spacing parameter. A suitable choice of T in [42] satisfies: if ae is small (diffusion-dominated), then r | 0 e ~ 0(K)2/K (3.9) Chapter 3. Non-Symmetric Streamline-Diffusion Method 33 and, if ae is large (convection-dominated), then r\Ue ~ 0(he/\a\2). (3.10) For example, when fte is a rectangle element, the optimal choice of r is r = l-he({ae)j\a\2 (3.11) with £(ae) = coth(a e) - ae'1. This leads to the following error estimate for the steady state case. T h e o r e m 3.1 Let amax = max e (a e ) . (1) Ifamax ^> 1 (convection-dominated), then \\uh-u\\ = 0{hk+1l2). (3.12) (2) If amax <C 1 (diffusion-dominated), then \\uh-u\\ = 0(hk+1). O (3.13) The error estimate in (3.12) is a half-order lower than in (3.13). It may not be improved in general (along the line of the argument in [64]), but in many cases, the sharp error estimate 0(hk+1) is often observed. 3.2.2 S D M e t h o d s for Sys t ems We now consider extensions of the scalar SD formulation for the convection-diffusion system (3.1). We assume that the problem is defined in ft = Md to ease our analysis. Remark 3.1 If problem (3.1) is defined in a bounded domain, we need to prescribe boundary conditions. Unfortunately, it is extremely difficult, in general, to find a set of Chapter 3. Non-Symmetric Streamline-Diffusion Method 34 well-posed boundary conditions for such an incompletely parabolic system. Only for some special cases, for instance the Friedrichs' system [47], are we able to specify boundary conditions suitable for a numerical analysis. • For later purposes, let us define [A] the 'pseudo-transpose' of A by [A]* = (i4*,i4*,... ,i4$). (3.14) Introducing AM = jttiAj + fi2A2 + ••• + HdAd, (3.15) where \i = (/zi, (J>2, • • •, fJ-dY with fa > 0 (j = 1 , . . . , d) is an arbitrary parameter vector satisfying |^| = 1, we assume in the following that A^ is diagonalizable, i.e., there exist a non-singular matrix Y^ and a diagonal matrix A^ such that Y;lA^A,Y, = A, (3.16) with diagonal elements forming eigenvalues Aj, (« = 1 , . . . , / ) , of AQ 1AfM. A finite element subspace analogous to (3.5 ) is given by Vh = {Vh | yh e (c°(Q,))\ Vh \Qee (Pk(ne))1}. (3.17) The counterpart of (3.6) for (3.1) in steady state reads: Find Uh € Vh such that for all (A • VU\ * ) + (JCVUh, W ) (Galerkin) nel + p ^ , P \ $ ) K (SD) = {FtV). (3.18) For the coupled system (3.1), we lose the vision of streamline directions. Therefore, an appropriate generalization of the operator Ph is not trivial for system (3.1), even in the Chapter 3. Non-Symmetric Streamline-Diffusion Method 35 symmetric case. Some previous attempts fail to treat each component equation properly even for a decoupled one-dimensional problem [43]. The two common nonzero choices of Ph (when AQ = h, an identity matrix) are Choice (a ) : d Ph(q) = TA-Vy = J2rAjyXj; (3.19) Choice (b): d Ph{q) = T[A]-W = ^TAtjyX], (3.20) 3=1 where r > 0. When (3.1) is symmetric, (a) and (b) are identical. At a first glance, choice (a) of Ph may seem to be more attractive than choice (b). This is because if we take a look at the artificial diffusion term in (3.18), i.e., (A • W \ Phm) = f (VUhyDsV^dx (3.21) Jn using choice (a), we find that the artificial diffusion matrix Ds — ArAi is symmetric positive semi-definite, which is normally expected as in the symmetric system case. In contrast, the artificial diffusion matrix Ds = AT A using choice (b) is indefinite. However, the choice of (a) was criticized, for example, by [42] and [83]. In [83], the choice of (b) was actually used to compute one-dimensional Euler flow. Numerical computations in [43] showed better results with choice (b) than with choice (a). With the detailed investigation that follows, we are able to understand the intrinsic behavior of the above choices. Furthermore, a suitable operator Ph can be generated for the incompletely parabolic system (2.6) that we are concerned with in semiconductor device simulations. For clarity, we carry out the discussion only for the finite element space Vh with uniform rectangular elements Oe. The obtained results can be extended in a straightforward manner to more general cases using element-wise local coordinate transformations (see [77, 40]). Chapter 3. Non-Symmetric Streamline-Diffusion Method 36 3.3 Semi-Discrete SD Formulations Let us consider the convection-diffusion system (3.1). The problem is defined in Q x (0, T] (with O = Md). Then we require that (3.1) be subject to the initial condition U{-,0) = Uo, (3.22) with Uo having compact support in Md. For convenience, we define Z^ = A0' Y^ so that Z^Av1/2A»Ao1/2Z» = A,. (3.23) The following assumption on the matrix K is needed for our discussion. A s s u m p t i o n 3.1 There exists a diagonal matrix Ap = diag(f3i,..., 0i), with /% > 0 (i = I,...,I), such that Aj Z~XAQ KAQ Z^A" 1 ^ is symmetric positive semi-definite. • The validity of this assumption is justified by Propositions 2.1 and 2.2 for the HD and SHD models. 3.3.1 One-Dimens ional Case For the one-dimensional case d = 1, (3.1) can be written as A0Ut + AUX = KUXX + F. (3.24) We first consider a simpler case, i.e., A0 = // and K = nil, where // is the I x I identity matrix. Then (3.24) can be diagonalized by the transformation U = ZV, to read Vt + AVX = KVXX + F, (3.25) where F — Z~^F. It is easy to see that (3.25) is fully decoupled. Each component equation of (3.25) has the form vu + XiVi}X = KvitXX + fi. (3.26) Chapter 3. Non-Symmetric Streamline-Diffusion Method 37 This is a scalar advection-diffusion equation. The semi-discrete SD-type variational form for (3.26) can be obtained using the SD method (3.6) with T,- given by (3.11) (replacing a in (3.11) by A;). Thus the optimal behavior of the scheme is attained. If we assemble (3.26) in vector form for all i (i = 1 , . . . , I), we obtain the SD formulation for (3.25) (AVX\ *) + (KVX\ «.) + £(AK* - KVI -F + Vt\ ATA$,)ne e= l = (F-Vth,$). (3.27) In (3.27), AT = diagfa, r 2 , . . . , n). (3.28) After employing the change of variables Vh = Z~1Uh and $ = Z*\]/, (3.27) converts to the semi-discrete SD formulation for the system (3.24), that is, {AUhx^) + (KU^^x) + ^2(jCUh,(Ar)^x)ne e = l = ( F - « / * , * ) , (3.29) where r = ZATZ~\ (3.30) Thus, the proper choice of Ph in this case should be Ph = {AT)^X. (3.31) It is clear now that choice (a) of Ph is not suitable and choice (b) of Ph, which indeed is better, needs to be further rescaled by a more general r . Next, we will consider the case when AQ ^ Ii and K is more general, possibly even singular. Under Assumption 3.1, our previous argument can be adjusted corresponding Chapter 3. Non-Symmetric Streamline-Diffusion Method 38 to this more general case. Similar to the derivation of (3.25), (3.24) becomes, by the transformation V — AJ Z~XA0 U, Vt + AVX = KVXX + F, (3.32) where F = AlJ2 Z^A^^F and K = A^Z^A^KA^ZA^2 (3.33) is symmetric positive semi-definite. Noting that (3.32) is a special symmetric system, the symmetric SD formulation for (3.32) can be derived according to [40], yielding (AVX\ $) + (KVX\ * , ) + J2(AVxh - KVt -F + Vt\ ArA$x)a. e = l = (F-Vt\$), (3.34) where AT can still be defined by (3.28) and its entries are defined by n = \ht(ai)/\\i\, (3.35) cti = -\\i\h/iti, (3.36) £(<*;) = coth(ai) - afl, (3.37) Ki = elAfz-'A^/2KAo1/2ZA-,1/2et = e\Z-xA^l2KAlxl2Zei. (3.38) In the above, re,- > 0 is ensured by Assumption 3.1. Employing the change of variables Vh = AXJ2Z-1 Al/2Uh and $ = A^1/2Z*Ao/2$, we obtain the semi-discrete SD formulation for system (3.24) {AUhx, tf) + (KUS, 9.) + J2(£U\ (Ar)«*,)n. e = l = (F-AoU?,V). (3.39) Chapter 3. Non-Symmetric Streamline-Diffusion Method 39 The matrix r in (3.39) is given by r = A-ll2ZKTZ-lA-112. (3.40) From the above, we obtain Ph in such a case as Ph(V) = ( A T ) * * . . (3.41) As we have observed, the above definition of r is independent of the particular form of A/?. Therefore, A^ is just an auxiliary matrix which is not involved in actual computations. 3.3.2 Mul t i -Dimens iona l Case In this subsection, we attempt to extend our construction of Ph from the one-dimensional case to the multi-dimensional case. A reasonable extension seems to be Ph(i&) = [A]T • V * . (3.42) The scaling matrix r above is defined by r = A^Z^diagfr, • • •, T O Z " 1 V 2 , (3-43) where for i = 1 , . . . , / , d 1 JT fid = - W l , . -~1)\ (3.44) 1 21 iU&/N2 (3.45) and M a = ( £ ( ^ ) 2 ) 1 / 2 > (3-46) £i = coth(ai) - aj1, (3.47) on = \\k\2hJKi (3.48) Ki = elZ^A^KA^Z^a. (3.49) Chapter 3. Non-Symmetric Streamline-Diffusion Method 40 Here A;j is the ith eigenvalue of AQ 1 AJ. The above K{ > 0 is again ensured by Assumption 3.1. The definition of Ph satisfies three design criteria suggested in [40], i.e., (1) For / > 1 and d = 1, it reduces to the optimal form presented in Subsection 3.3.1. (2) For / = 1 and d > 1, it coincides with that in (3.11). (3) For / > 1, d > 1, and simultaneously diagonalizable coefficient matrices Aj, it satisfies criterion 2, with respect to each uncoupled scalar equation. In addition to the above three criteria, the definition of Ph also satisfies (4) If the system is symmetric, it gives an analogue to those defined in [40] in the sense that AT A* is symmetric positive semi-definite, since the matrix r given by (3.43) reduces to r = Yltddiag(T1,---,Ti)Y^d. The verification of claims (l)-(4) is not difficult and is omitted here. Remark 3.2 For the two-dimensional (d=2) HD model of the semiconductor device equations with the (singular) diffusion matrix K given by (2.12), the corresponding Ki in expression (3.^9) can be obtained explicitly as KX = 0 Ki = I K2 = f (3-50) 3.4 Full Discre te SD Formulations Let us define B(Uh, tf) = (A • VU\ * ) + ()CVUh, Vtf) + £ ( A • VUh- V • (/CVt/'1), P / l(^))n c(3.51) e = l Chapter 3. Non-Symmetric Streamline-Diffusion Method 41 Figure 3.1: A space-time slab, and JF($) = ( F , $ + P /1(^)) (3.52) then the semi-discrete SD formulation (3.18) can be written as (A0Uth, $ + Ph(V)) + B(Uh, tf) = F{V). (3.53) 3.4.1 Space-Time Finite Element Methods Simply consider a uniform partition 0 = to < <i < • • • < % = T of the time interval I(T) = (0, T). Denote by In(T) = (tn, tn+1) the nth time subinterval with At = tn+i — tn Chapter 3. Non-Symmetric Streamline-Diffusion Method 42 being the length of the subinterval. A space-time "slab" is then defined as (see Fig 3.1). Qn = ftxIn(T). (3.54) For the n th space-time slab, we define space-time element domains as Qen = ne x In(T) (3.55) (e = 1 , . . . , nei). We define the following finite element subspace Th = {yh | yh e {C°(Qn))1, Vh | Q e G (Pk(Qen))1}. (3.56) Note that any test function in Th is continuous within each space-time slab Qn, but is discontinuous across the interface of slabs. Considering this, let Vh{t±)=\im±Vh{tn + e), (3.57) and the jump in time of Vh lVh(tn)J = Vh{t+n) - V\rn). (3.58) Given Uh(to), a projection of the initial value UQ on Vh, the space-time SD formulation is as follows. Within each Qn (n = 0 , . . . , N - 1 ) , find Uh G Th such that for all $ G Th, f (A0Uth-y + A-VUh-y + lCVUh-Vi£)dxdt (Galerkin) JQn + I AolUh(tn)} • #(*+) dx (Jump Condition) + J2 £Uh-Ph{®)dxdt (SD) e t l JQen = [ F-^dxdt. (3.59) The operator Ph in (3.59) can be determined using the results in Section 3.3 as we show below. By the change of variables Uh = A^1/2Vh and $ = AQ 1 / 2 $ , (3.59) becomes / (Vth • $ + A • W f c • $ + £ W f c • V $ ) tfocft Chapter 3. Non-Symmetric Streamline-Diffusion Method 43 + / lVh(tn)j • $(*+) dx + J2 [ CVh- Ph{$) dxdt JV e=l JQn = I F • Qdxdt, (3.60) JQn where, Ph = A^ 1,2PhA^1/2, F - Ao1 /2FA^1 /2 and A* = (A,, -.,&) = ( 4 I / J M - 1 / a , • • •, A0-1/2A,A0"1/2) £ = rf*oflF(JKr, • • •, If) = diag(Ao1/2KAo1/2, •••, A^ 1/2K A^1'2). First we consider the case when At = h. Let d+l ^+1 = 7d+iiT^~^)t' (3'61) The operator Ph in Subsection 3.3.2 is only suitable for the steady state case. However, we can incorporate it in the space-time finite element methodology. We take Ph in (3.60) to be Ph($) = [A]? • V$ + f • $ t (3.62) with f = Zmidiag(T1,---,-n)Z£¥1, (3.63) if we considered (3.60) to be a formulation for the "(d+l)-dimensional" steady state convection-diffusion system. In (3.63), r4- (i = 1 , . . . , I) can still be defined by (3.45), but | A,-12 should be changed to |A«|2 = ( l + E ( ^ ) 2 ) 1 / 2 ) (3-64) i= i where A,-j is the zth eigenvalue of Aj. We notice that the matrix (Ii + J2i=i A~i)/Vd + 1 is diagonalizable by its eigen-matrix ZMd+1, as is the matrix (J2i=i A~i)/vd- Hence Z,d+1 = Z.d. (3.65) Chapter 3. Non-Symmetric Streamline-Diffusion Method 44 If At ^ h, we scale the time subinterval by h/At. Therefore, the above formulation remains the same except that |A,-|2 in (3.64) must be revised to l^ -|2 = ((^)2 + E(^)2)1/2-However, we prefer to make the following change heuristically \M2 = mh2+E(h3)2)1/2, (3.66) where Cr = 1 except Cr = 0 for constant-in-time space-time finite elements. Transforming from (3.60) back to (3.59), we obtain the corresponding term Ph in (3.59). To summarize, we give the expression for Ph as follows. ' * ($ ) = r = fid = Ti = |A.-|a = 6 = Cti = ACj = [A]T -V^ + T-^t, A^^Z^diagfa, • • •, T{)Z~]AI1/2, 1 > ' - ' i y " o^&'/Ma (a(^)2 + E(^)2)1/2, coth(aj) — a,"1, x|A8-|2VK*) (3.67) (3.68) (3.69) (3.70) (3.71) (3.72) (3.73) (3.74) Here, A tj is the zth eigenvalue of AQ1AJ. We point out that for general non-symmetric convection-diffusion systems, the so-called Galerkin/least-squares algorithm introduced in [35] is not appropriate, since we have to use [A] in the definition of Ph instead of A in such cases. Chapter 3. Non-Symmetric Streamline-Diffusion Method 45 3.4.2 M e t h o d s of Lines Let JJ be the approximate solution vector of Uh at the time tn = nAt with a uniform time step At. With U being a projection of the initial value UQ on Vh, we are interested in two simple time discretizations below. SD-backward Euler: TTh(n+1) _ Jjhin) (A0- ^ , tf + Ph(V)) + B{Uh{n+1\ tf) = ;F(tf), (3.75) for any $ e V f t . SD-trapezoidal: (A>- ^ , * + P A (*)) + B(-(Uh{n+1) + /7&(B)), $ ) = ;F(tf), (3.76) for any * G Vfe. The above schemes are not strict SD schemes in the sense of the space-time finite element idea, but they are simple to use. 3.5 Stabi l i ty A r g u m e n t s From now on, we consider k = 1 in definitions of the finite element subspaces (see (3.17) and (3.56)). For defmiteness, we set F = 0, i.e., no external force is prescribed. 3.5.1 Space-Time Fini te Element M e t h o d s To simplify our analysis, we carry out the stability discussions only for the constant-in-time space-time SD scheme without loss of generality. It is convenient for our discussion to set Ao = Ii, otherwise we can simply eliminate it by working on the switched formulation (3.60). With AQ = I}, the constant-in-time version of (3.59) is reduced to Chapter 3. Non-Symmetric Streamline-Diffusion Method 46 where Uh(n+l) = Uh(t+) and Uh(n) = Uh{t~). Recall that r is given by (3.68) and ZM is the eigen-matrix of the eigenvalue-eigenvector problem (A M - \I,)ZMj = 0. (3.78) The stability analysis of (3.77) is not easy in general, but at least we can handle it if the following normal symmetrization exists. Suppose there exist non-singular matrices B and C such that A0 = CB (Symm. Pos. Def.) (3.79) Ai = CAiB (Symm.) i = l,...,d (3.80) K = CKB (Symm. Pos. Semi-Def.) (3.81) By the change of variables Uh{n) = BVh(n) and * = (7*$, (3.77) becomes At (A • VVh{n+l\ $ ) + At (JCVVh{n+1\v$) + {AQ(Vh{n+1) - Vh{n)), $ ) + A * ( i - W / l ( " + 1 ) , A f - V $ ) = 0, (3.82) where f = B^rC'1 = B^Z^ diagfa, • • • ,rt) Z^C'1. (3.83) Considering the transformed eigenvalue-eigenvector problem ( i o 1 / 2 A ^ i o 1 / 2 - A J , ) ^ , , = 0, (3.84) we claim that the above f is appropriate, i.e., it can be expressed by f = i 0 - 1/2Z»d diagin, • • •, n) Z;lA-Q112 (3.85) which corresponds to the transformed symmetric system. Thus, f given by (3.83) is a symmetric positive definite matrix. It is easy to see, by comparing (3.78) and (3.84), Chapter 3. Non-Symmetric Streamline-Diffusion Method 47 that Z,d = BAl/2Z,d = B(CB)-^2Z,d. (3.86) Using the above, f given by (3.85) turns into f = GB^I\{CB)^B^Z,d) diagfa, • • •, n) {{CBf^B^Z^Y^CB)-1^ = B-1ZIMddiag{r1,-..,Tl)Z;1dC-\ (3.87) which proves the claim. The stability of scheme (3.82) is now easy to prove. Rewrite (3.82) as (A0(Vhin+1) - Vh(n)), $ ) + AtB(Vh{n+1\ $ ) = 0. (3.88) It is obvious that g(yh(n+l) yh(n+1)\ _ (A.Tjyh(n+1) yh(n+1h i (XVV^" + 1 ^ v y A ^ n + 1 h + (A-VVh(n+1\Af -VVh(n+1))>0. (3.89) If we take $ = Vh{n+1) in (3.88), then l l ^ ( n + 1 ) | U < | |^ ( n ) | | i o , (3.90) where || • \\Ao is a weighted L2-norm defined by \\V\\2AQ = {V, A0V). The stability of (3.82) and therefore the stability of (3.77) are proven. We are interested in two symmetrizing cases: Case 1: C = I{ This symmetrization corresponds, for example, to those for Navier-Stokes equations and the semiconductor HD model (see [72, 77]). C a s e 2: C = B'1 The special situation in this case is that of simultaneously diagonalizable systems. If As-sumption 3.1 holds, then the stability of scheme (3.59) can be obtained for 1-D problems. We summarize the above discussion in the following theorem. Chapter 3. Non-Symmetric Streamline-Diffusion Method 48 Theorem 3.2 If problem (3.1) admits a normal symmetrization, then the space-time SD discretization scheme (3.59), with the operator Ph defined by (3.67)-(3.74), is uncondi-tionally stable. • 3.5.2 Methods of Lines The stability proof for the SD-backward Euler scheme is more difficult than for the space-time SD scheme due to the existence of the term (A0(Uh(n+l) - Uh(n))/At, Ph($)) in (3.75). We can only provide a proof for the one-dimensional case under Assumption 3.1 and r = 8hA$l with the scalar parameter 6 > 0. By transformations Uh{n) = A^1/2ZA^1/2Vh{n) and $ = A„ 1/2Z~*Aj /2$, (3.75) yields yfcOH-l) _ yh(n) yh(n+l) _ yh{n) (- K T — ' *) + Sh (- K T ^ ' A $ * } + {AV* ' $ ) + (ffV?(B+1), **) + ^ (AK" (n+1\ A$,) = 0, (3.91) where, the symmetric positive semi-definite matrix K is defined by (3.33). Taking $ = yh(n+i) i n ^ g i ^ w e h a y e | i y^ ( n + 1 ) | | 2 I I F ^ I I 2 yh("+V _ yhW I ,| III || , sh{V V ^ ^ ( n + D j 2At 2At v At + l l ^ ^ f e + Sh \\AVxh{n+1)\\2 < 0. (3.92) Taking $ = 6h (Vh{n+1) - Vh{n))/At in (3.91), we have _ |^( fc(H+1) _ fc(»)||a (v—_!LA^M, A t 2 " II • v A i , M n,^(^), |2. _ ii^Wig-} + M ! ( | | A ^ ( " + 1 ) | | 2 - ||AV,fc(B)||3) < 0.(3.93) 2Atv" x "K " * "K/ 2At We define 111 • 111, an auxiliary norm, by m i l 2 = ||V||2 + ^ | |AK | | 2 + {6hf\\Vx\\\. (3.94) Chapter 3. Non-Symmetric Streamline-Diffusion Method 49 Then (3.92) + (3.93) yields 2 A i u " ' At +-^-2\\Vh{n+1) - Vh{n)\\2 + Sh | |A l^ ( n + 1 ) | | 2 < 0. (3.95) Applying the inequality a2 + b2 > 2ab to the last three terms at the left hand side of the above, we conclude that | | | ^ f c ( B + 1 ) | | | 2 < | | | y * ( n ) | | | a . (3.96) To summarize, we have the following theorem. T h e o r e m 3.3 If Assumption 3.1 holds and Ph = ShA^1 A^x with a scalar parameter 8 > 0, then the one-dimensional SD-backward Euler scheme (3.75) is unconditionally stable. • We notice that the SD-backward Euler scheme (3.75) is equivalent to the constant-in-time space-time SD scheme (3.59) if we drop the term {A0(Uh{n+1) - Uh{n)), Ph(V)) (3.97) in (3.59). As we can see from the above stability discussion, the term (3.97) is difficult to handle in the analysis. The difficulty caused by this term is also prominent when we try the stability analysis for the SD-trapezoidal scheme (3.76) in the 1-D case. A result similar to Theorem 3.3 can be obtained but under an extra restrictive condition (see Appendix B). Existence of such a term is not only troublesome in the analysis but also exhibits poor behavior in computations, as we will see in the next section. 3.6 Numer ica l Tests For the purpose of numerical comparisons, we use three different choices of P . Choices (a) and (b) are given by (3.19) and (3.20), respectively, with r = ShAg1. We will refer Chapter 3. Non-Symmetric Streamline-Diffusion Method 50 to our generalized Ph (3.67)-(3.74) as the (c)-type. Better performance was observed in [43] for the (b)-type than for the (a)-type. This can be explained now by our analysis in the previous sections: The (b)-type operator does have the upwinding effect in the correct direction. The (c)-type operator is different from the (b)-type in that it also gives a proper scaling. Two examples are tested here. The first is a one-dimensional hyperbolic system with constant coefficients. The next is the nonlinear hydrodynamic semiconductor device model for both one and two space dimensions, in which case K is singular. 3.6.1 Tes t 1 Our test hyperbolic system is given by Ut + AUX where ' l 1 + ^ 0 (3.98) (3.99) / The eigenvalues of A are 1, — e with eigen-matrix Y *1 1 ^ 0 •1 (3.100) The boundary condition is chosen to be periodic. System (3.98) can be diagonalized by the transformation V = Y~1U. The initial-value function UQ is set corresponding to the initial value V0, a pulse function, of the V-system. In this example, we take h = 0.01 and At = 0.005. The SD operator for the (c)-type is given by Ph($) = ( A T ) ' ^ according to (3.67)-(3.74) and Chapter 3. Non-Symmetric Streamline-Diffusion Method 51 Numerical tests are presented with e = 0.01. A typical value 6 = 1 is taken for the (a)- and the (b)-type operators. Comparison of the three types is shown in Fig. 3.2 and Fig. 3.3 using the constant-in-time space-time SD method (3.77). It is evident that the (c)-type gives the best result while the (a)-type gives the worst, with respect to oscillatory behavior. The solutions of the (c)-type are clearly non-oscillatory. We can show that the (c)-type scheme is actually a monotone scheme in such a case. To verify this, it is sufficient to consider a scalar hyperbolic equation ut + aux = 0 ( a > 0 ) (3.102) Mn) u,~ ™i„^ ~f „Mn) since (3.98) can be virtually decoupled. If we define uf , the value of uk at the node Xi, the (c)-type scheme can be easily expressed by the following finite difference scheme employing the trapezoidal quadrature rule h{n+l) . v (1 + vQ)u\ ^(i-^+r-^+«K-, h (n+1) _ uh(n) (3.103) where v — aAt/h and 6 = a/JCr(h/At)2 + a2. If we denote the solution vector U^ = ( M ^ ( n ) , ^ ( r i ) , . . .)*, then, by choosing Cr = 0, (3.103) reads with M ' 1 + v -v 1 + v (3.104) (3.105) The inverse of M is M - i 1 + v — 1 (3.106) Chapter 3. Non-Symmetric Streamline-Diffusion Method 52 5 2.0 -1.5 -1.0 -0.5-0 -Legend Approx. Exact 0 0.5 TYPE-(a) ^ 1.0 X \ 1.5 2.0 *"' 2.0 -1.S -1.0 -0.5 --0.5 -i Legend Approx. Exact TYPE-(b) \f / u \ "' 2.0 -1.5 -1.0 -0.5 -TYPE-(c) Legend Exact : J\ k-; \ Figure 3.2: U\ by (a), (b) and (c)-types of operators in test 1 using the constant-in-time space-time SD scheme. Chapter 3. Non-Symmetric Streamline-Diffusion Method 53 -0.5 -1.0-1,5-Legend Approx. Exact TYPE-(a) F 1 \ eg => 0 --0.5 --1.0-Legend Approx. Exact i ! 0 0.5 TYPE-(b) 1.0 X 1.5 2.0 Legend Approx. Exact ! i Figure 3.3: U2 by (a), (b) and (c)-types of operators in test 1 using the constant-in-time space-time SD scheme. Chapter 3. Non-Symmetric Streamline-Diffusion Method 54 Figure 3.4: Solutions by (c)-type operator with Cr = 1 in test 1 using the constant-in-time space-time SD scheme. Chapter 3. Non-Symmetric Streamline-Diffusion Method 55 2.0 -1.5 -5 1.0-0.5 -( Legend Exact SD-BACKWARD EULER /• i / l l \ ) 0.E 1.0 X 1.5 2.0 SD-BACKWARD EULER -1.0 --1.5 Legend Exact / 1.0 Figure 3.5: Solutions by (c)-type operator in test 1 using SD-backward Euler scheme. Chapter 3. Non-Symmetric Streamline-Diffusion Method 56 Figure 3.6: Solutions in test 1 using the Galerkin-backward Euler scheme. Chapter 3. Non-Symmetric Streamline-Diffusion Method 57 which has non-negative entries. This proves the monotonicity of (3.103). If we use Cr = 1, instead of Cr = 0, the monotonicity is violated and therefore oscillations cannot be avoided. This is demonstrated in Fig. 3.4. Remark 3.3 For multi-dimensional convection-diffusion problems, the monotonicity of the SD method unfortunately does not hold. To further control the oscillatory behavior of the numerical solutions, the so-called nonlinear shock-capturing modification may be needed to introduce an amount of crosswind dissipation (see [15]). • Methods of lines are not suitable in the SD formulations as mentioned in the previous section. In this example, the SD-backward Euler scheme with the optimal Ph is also tested to support our argument. The numerical results are shown in Fig. 3.5, which have oscillations. As a reference, numerical solutions using the Galerkin-backward Euler scheme are provided in Fig. 3.6. 3 .6.2 Test 2 Since the HD model is a system of nonlinear PDEs, a so-called predictor/multi-corrector algorithm with implicit-explicit nonlinear Newton type iterations is used (see [43, 77]). We delay the description of such a procedure to Chapter 5. Physical parameters of the model needed for the following computations are presented in Table 5.1. Extensive simulations of semiconductor devices will be carried out in Chapters 5 and 6. We denote ipass the number of iterations required at each time step to relax to solutions within the design error tolerance, tol = 10 - 4 . We set ipass < 40; if more iterations are required then the scheme is considered to be divergent. In this example, the lattice temperature is set to be 300K. 1-D C a s e : Chapter 3. Non-Symmetric Streamline-Diffusion Method 58 A l.O-fim p+-p-p+ silicon diode is used with a 0.25-//m source, 0.5-//m channel and 0.25-/mi drain regions. The doping profile pD is given by pr> = 101 7cm - 3 , if 0 < x < 0.25 or 0.75 < x < 1.0; pD = 1015cm-3 , if 0.35 < x < 0.65. The above segments of pD are connected smoothly by polynomial fits. The doping in this case is relatively mild. A uniform mesh of 100 points is used for this example. A comparison has been made between results of the (a)-type (S = 1) and the (c)-type. We first apply IV biased voltage to the device. The steady state is reached by solutions of both types around time t = 6. We have observed that the iteration number ipass for the (a)-type is almost as twice as that for the (c)-type. When a 2V biased voltage is applied, the (a)-type diverges around time t = 0.22 and the (c)-type reaches the steady state at t ime t = 6. Numerical solutions for a 2-volt biased voltage are shown in Fig. 3.7. 2-D C a s e : Let us consider a two-dimensional p+-p-p+ diode shown in Fig. 3.8 with the device domain [0, 1.0(/^m)]x[0, 0.5(^m)]. We set a 0.2-fim source and drain and the doping profile is defined by pD = 5 x 1017cm~3 in [0, 0.25(^m)]x[0.25, 0.5(//m)] and in [0.75, 1.0(/«in)]x[0.25, 0.5(/im)] and PD = 2 x 10 1 5cm - 3 elsewhere, with abrupt junctions. A uniform mesh of 80 x 40 points is used. In this case, a IV bias is prescribed. The (a)-type scheme fails to converge (producing an unphysical negative electron density). For the (c)-type, the solution tends to the steady state about t = 10. Typical solutions of this type are shown in Fig. 3.9. The constant-in-time scheme is used in the above. It is first-order in time, but nev-ertheless it can be used to obtain the steady state solution with an overall accuracy of 0(h3/2). Moreover, we observe that in practice errors in spatial discretization often domi-nate those of the time discretization. It may be argued that a linear-in-time scheme would be preferable over the constant-in-time one, despite its additional complexity. However, we use the predictor/multi-corrector algorithm with explicit correction passes, so the Chapter 3. Non-Symmetric Streamline-Diffusion Method 59 scheme is not unconditionally stable any more. The analysis in [75] showed that the linear space-time method might cause instability when the advection dominates. We also have tested the above problems using the SD-backward Euler scheme. The evidence indicates that the SD-backward Euler scheme is even less stable than the clas-sical Galerkin-backward Euler scheme. In contrast, the space-time SD scheme is much better. The method of lines with the semi-discrete SD approach is often used in the lit-erature. For instance, they were employed for scalar equations in [15, 32] and for systems in [43]. However, the discussion and numerical experiments in this chapter suggest that we should use the space-time SD schemes, especially for systems. Chapter 3. Non-Symmetric Streamline-Diffusion Method 60 1.25 -1.00-1 0.75-I eg 0.50-0.25-Legend Velocity Current Temperature 0 F - - * - - - - "*—t-0 0.25 TYPE-(c) \ r~ \ \ \ \ j V'--'z~'--::H ' " ' ' • ! • . . . . • • | - . . . - i . . . - T f t r . . . . . . " * — 1 " " ' 0.50 0.75 X 1. DO Figure 3.7: Solutions by (a) and (c)-types of operators in test 2 using the constant-in-time space-time SD scheme. Chapter 3. Non-Symmetric Streamline-Diffusion Method 61 source P+ P drain P + o Figure 3.8: Geometry of a 2-D p+-p-p+ diode in test 2. Figure 3.9: Electron density (left) and electron temperature (right) in test 2 using the constant-in-time space-time SD scheme. Chapter 4 Potent ial Equation Solver Let us consider the d-dimensional potential equation (2.2) in the form ' V-E = L, (4.1) _ E = - V ^ in i l , assuming for this purpose that p is a known function. The equation (4.1) is subject to the mixed boundary conditions 0 | r , = <?o, (4-2) E • n |r2 = 0. (4.3) Here, E is the electric field and go is given by (2.35) and (2.37). The discussion that follows can be directly extended to the case where E = —Wip with D{x) > 0 known. But this extension is left to the reader because it is not strictly needed for our purposes. 4.1 Introduct ion Upon elimination of E in (4.1) one obtains Poisson's equation for tp. However, in our context, E is a more important physical variable than tp. As far as the typical HD model is concerned for example, E is the only quantity that is needed in the system of conservation laws (2.1). Since the system of conservation laws is sensitive to the force of the electric field, an accurate approximation to E is desired. Finite element methods with high degree piecewise polynomials, of course, will achieve this goal. However, to 62 Chapter 4. Potential Equation Solver 63 avoid the unnecessary complexity of using different settings in the finite element spaces, we prefer to stick with linear finite elements as used for the system of conservation laws. The classical Galerkin method only provides first-order accuracy for E. Thus it will cause a corresponding order reduction for all solution components even if the order of the SD method is higher. Therefore, a more accurate method whose order is compatible with that of the SD method is needed for our special purpose. Normally, in order to obtain an accurate solution of E, we should treat E as an independent variable. It is natural to consider a mixed finite element method for (4.1). Unfortunately, this requires that the finite element spaces for tf> and E be compatible, i.e., that the space pair satisfies the Brezzi-Babuska condition (see [10]). Different settings of nodal points (staggered grids) have to be used for ip and E. To circumvent this constraint, some alternative approaches were proposed, for example, the augmented finite element method in [10] and the div-rot finite element method in [52, 53]. The augmented method solves a coupled system for if) and E simultaneously, so it is costly. The div-rot method solves a system of equations for E only, hence it could be more efficient than the augmented method. However, to guarantee the convergence of the div-rot method, the boundaries Fi and T% must be connected, i.e., the segments of I \ cannot be separated by T2, and vice versa. Unfortunately this does happen in our device simulation problems. We therefore propose a new method that shares the advantages of both the augmented and div-rot methods. The arrangement of the chapter is as follows. In Section 4.2, we provide formulations of the augmented and div-rot methods. We then propose a new formulation in Section 4.3 by combining ideas of these methods and carry out numerical experiments in Section 4.4. Chapter 4. Potential Equation Solver 64 4.2 Prev ious Work 4.2.1 T h e D i v - R o t M e t h o d The div-rot method was initially designed to solve the div-rot system such as the Maxwell equations using the so-called least-squares methodology. If we apply the rotation operator to the second equation of (4.1), we have V x E = 0. Then (4.1) can be converted to the div-rot system V - £ - Z = 0, (4.4) V x E = 0, which eliminates the variable i\). A finite element approximation for (4.4) was discussed by Kfizek and Neittaanmaki (see [52, 53]). Before introducing the method, we define the finite element subspace & = tf | q* € (C°(n))d, qh \Qee (Pi(O e)d , qh X n | r , = 0, qh • n | r , = 0}. (4.5) If we assume that go is constant, then the div-rot formulation reads: Find Eh € Qh such that for all qh G Qh a(Eh,qh) = (L,V-qh), (4.6) where <P,q) = (V -p, V • q) + (V X p, V X q). (4.7) To ensure convergence, the bilinear form a(-,-) must satisfy the Q^-ellipticity or the Lax-Milgram theorem (see [52, 18, 82]). Kfizek and Neittaanmaki showed that for all qh G Qh, ||fffc||<C7(||VVll + l|Vx^||), (4.8) Chapter 4. Potential Equation Solver 65 for some constant C independent of h, if and only if Fi and T2 are connected. 4.2.2 T h e A u g m e n t e d M e t h o d The idea of the augmented method originated from that of Franca and Hughes's Galerkin least-squares method in [23] (which is a counterpart of the SD method). Define a pair of finite element subspaces Vl£ = {uk\wheC0(n),wh\Q.eP1(ne),wl>\r1=g0}, (4.9) Qh = {qh\qhe(C°(n))d,qh\Qee(Pi(Sle))d,qh-n\r2=0}. (4.10) The augmented formulation reads (see [10]): Find (iph,Eh) G Wg0 x Qh such that for all (wh,qh) e w j X Qh a{E\ qh) + b(i/>h, qh) = /3(L, V • qh) - a JTi g0qh • nds, k b(wh,Eh) + c(4>h,wh) = (L,wh), (4.H) where a(p,q) = a ( p , g ) + ) 9 ( V . p , V - g ) , (4.12) b(w,q) = -(l-a)(w,V-q), (4.13) c(v,w) = a(Vv,\7w). (4.14) The parameters a and (3 satisfy 0 < a < 1 and (3 > 0, respectively. Note that we have a coupled system to solve. R e m a r k 4.1 Motivated by the div-rot method, we also can consider the bilinear form a(-, •) as a(p,q) = a(p,q) + P[(V-p,V-q) + (Vxp,Vxq)}, (4.15) which enhances the stability of the method. If a = 0, the div-rot method is recovered. • Chapter 4. Potential Equation Solver 66 4.3 A N e w M e t h o d Based on the above work, we propose a new formulation that can be viewed as a balance of the augmented method and the div-rot method. The new approach reads: Find i>h € Wg0 and Eh G Qh, such that for any wh € WQ and qh G Q&, where a (£* , qh) + 6 t y \ ^ ) = 8(L, V-qh)- fTi g0qh • nds, _ c(iph,wh) = (L,wh), <P,q) = (p,?) + / 9 [ ( V - p , V - 9 ) + ( V x p , V x g ) ] , 6(u;,g) = - ( u ; , V - g ) , c(u,w;) = (VujVio). (4.16) (4.17) (4.18) (4.19) Notice that (4.16) is actually decoupled, therefore it saves computation time (as compared with the augmented method). The parameter 8 determines the order of accuracy of the method as we will see later. R e m a r k 4.2 The formula (4-16) can be interpreted formally as the finite element ap-proximation of the problem where -8(AE - VL) + E + VV> = 0, -At/> = L, AE = V V • E - V x V x E. (4.20) (4.21) If ft —> 0, the classical Galerkin method is obtained. If 8 —> oo, the div-rot method is obtained. In case that I \ and Y2 are connected, (V • p, V • q) + (V x p, V x q) is positive (cf. (4.8)). • Chapter 4. Potential Equation Solver 67 Denoting e,/, = if> — iph and eE = E — Eh, we are ready to state the following theorem. Theorem 4.1 If /3 ~ h and (if>,E) are smooth enough, then \\e4 = 0(h2), (4.22) \\eE\\ = 0{h3'2). (4.23) Proof. Since (4.16) is decoupled, it is a straightforward matter to obtain (4.22) (see, e.g.,[18]). Now we prove (4.23). It is obvious that e^ and eE satisfy a(eE, qh) + b(e^ qh) = 0, Vqh G Qh. (4.24) Let us define TlhE € Qh as the interpolant of E, say at nodal points, and denote nE = E — II/i-B, then the following standard approximation result holds: \\VE\\ + MIIV • VE\\ + ||V x r)E\\) = 0(h2). (4.25) By (4.24) and the Schwarz inequality, a{eE,eE) = a(eE, (E - UhE) + (ILhE - Eh)) = a{eE,riE + {UhE-Eh)) = a(eE,VE)-b(ei),I[hE-Eh) = a(eE,i]E)-b(e^,eE) + b(e^,T]E) < \\eE\\ \\VE\\ + 0 (||V • rjE\\ ||V -eE\\ + ||V x r,E\\ ||V x eE\\) + ||e^|| | |V-eB | | + ||e^|| | | V - J / B | | . (4.26) Using (4.22) and (4.25), we derive estimates for the terms on the right hand sides of the above inequality as follows \\eE\\\\vE\\<Ch4+l-\\eE\\2 (4.27) Chapter 4. Potential Equation Solver 68 P (l|V • r)E\\ IIV • eE\\ + || V x nE\\ IIV x eE\\) < C/3h" + 0 ( i | |V • e £ | | 2 + I | | V x e,;||2), (4.28) II^H HV-e^JI < C ^ + | HV-e^H2 (4.29) and I M | ||V^|| < C7*3. (4.30) The above positive constant C may depend on the solution ifi and E, but it is independent of h. Hence, we obtain h* I a(eE, eE) < C{h4 + /3h2 + — + hz) + -a{eE, eE). (4.31) Therefore a(eE, eE) < 2C{h4 + ph2 + j + h3). (4.32) The theorem is proved if we choose /? ~ h. • Remark 4.3 From the proof of Theorem 4-1, it is easy to see that the convergence result still holds if we consider (4-16) without the rotation term (V x E , V x qh), which is used to stabilize the scheme. • Theorem 4.1 indicates that the convergence rate of the scheme (4.16) is consistent with that of the usual SD method. The scheme (4.16) adopts the idea of the least-squares finite element method, but the methods are different. In a recent paper by Pehlivanov et al. [63], a least-squares method was proposed and analyzed. The preconditioner there, which decouples the calculations for i^h and Eh, is quite similar to our method if we take j3 = 1 and drop the rotation term in (4.16). That method thus gives first-order accuracy in \\eE\\, as was shown in [63] and also can be seen in the proof of Theorem 4.1. Chapter 4. Potential Equation Solver 69 We point out that the proof of Theorem 4.1 can be applied to a more general case in which a higher degree polynomial space Pk is used in the subspace pair (4.9) and (4.10). In such a case, the convergence results will be ||e^|| = 0{hk+1), (4.33) \\eE\\ = 0{hk+x'2). (4.34) 4.4 Numer ica l Test Let O C M2 be a unit square with I \ = T{ UTf, where T\ = {x = 0} and Tj = {x = 1}, and ]?2 = r \ T i . Consider the following model problem: § 7T —Atp = -7r2sin(—x) cos(7ry) (4.35) 4 v2 with mixed boundary conditions ^ In = °- (4-36) if) |ra = cos(7ry), (4.37) V ^ • n |ra = 0. (4.38) The exact solution is provided by T/I = sin(—x) cos(iry). (4.39) We use uniform partitions of fi into rectangular elements. The finite element subspaces Wg0 and Qh consist of piecewise bilinear polynomials. Denoting E = — V ^ , Table 4.1 and 4.2 exhibit the numerical results for Eh using (4.16). They confirm that 0 ~ h gives the best convergence results. We observe that the numerical convergence rates for our uniform mesh are higher than those predicted in Theorem 4.1 for general meshes. This is expected and can be explained along the arguments of Peterson in [64] and Lin Chapter 4. Potential Equation Solver 70 Mesh Size h = 0.2 h = 0.l h = 0.05 h = 0.025 e^ll 0 = 0 9.100D-2 3.125D-2 1.038D-2 3.556D-3 Without Rotation Term 0 = h 2.529D-1 4.280D-2 7.191D-3 1.384D-3 0 = 1 7.523D-1 2.607D-1 7.352D-2 1.934D-2 With Rotation Term 0 = h 1.807D-1 3.439D-2 6.059D-3 1.254D-3 0=1 3.905D-1 1.747D-1 5.577D-2 1.488D-2 Table 4.1: Errors of E with different parameter settings. in [56]. The computed solution can have a higher convergence rate in many situations as compared with that from the theoretical estimation, except for some special mesh constructions. In particular, sharper error estimates are possible for the rectangular meshes by more subtle investigations (see [56]). Moreover, the bilinear mode xy in the finite element subspaces often raises the convergence rate, which is beyond the reach of our error estimations. Chapter 4. Potential Equation Solver 71 Mesh Size h = 0.1 h = 0.05 h = 0.025 Convergence Rate 0 = 0 1.5420 1.5900 1.5454 Without Rotation Term 0 = h 2.5629 2.5733 2.3773 0 = 1 1.5289 1.8261 1.9265 With Rotation Term 0 = h 2.3935 2.5048 2.2725 0 = 1 1.1604 1.6473 1.9061 Table 4.2: Convergence rates for Eh with different parameter settings. Chapter 5 Simulat ion I: the p+-p-p+ D i o d e 5.1 Introduct ion The p+-p-p+ diode, which simulates a channel of a metal oxide semiconductor field ef-fect transistor (MOSFET), is a fundamental device which was used to understand the behavior of the hydrodynamic (HD) model and to optimize the HD modeling by simula-tions (e.g., overcoming the spurious velocity overshoot [30]). Recall that the HD model for semiconductor devices consists of a system of conservation laws for charge density, momentum and energy. Since there are hyperbolic modes (see Gardner et al [29] for a de-tailed discussion), it supports shocks or discontinuities of solutions. As a model problem, the one-dimensional p+-p-p+ diode was considered intensively by many researchers. For instance, uniqueness of subsonic solutions was proved by Degond and Markowich [19]; transonic solutions were analyzed mathematically by means of viscosity approximations by Gamba (see [24], [25]); Ascher et al [2] and Markowich and Pietra [59] studied shock phenomena for a current driven problem using a phase plane analysis, which gave a full explanation of how and when the shock occurs; shocks were found for the full system by Gardner [26] in a steady state case and by Fatemi et al [22] in a non-steady state case with low temperatures. There are many numerical methods designed to handle the occurrence of shocks. For example, the so-called ENO (essentially nonoscillatory) scheme used in [22] is an efficient shock-capturing algorithm. But most of them are one-dimensional in nature although 72 Chapter 5. Simulation I: the p+-p-p+ Diode 73 they can be used for two-dimensional problems by operator splitting techniques. The streamline-diffusion (SD) method is suitable for handling multi-dimensional cases with theoretical backing. In this chapter, the SD method presented in Chapter 3 with a shock-capturing modification is used to simulate the silicon p+-p-p+ diode for both the subsonic and transonic electron flows. Shocks of the transonic flows are captured for one-and two-dimensional diodes. In Section 5.2, a SD method with a shock-capturing modification as well as a nu-merical treatment of the "in-flow" boundary condition is presented. In Section 5.3, we describe the so-called predictor/multi-corrector algorithm applied to our problem. Phys-ical parameters for device simulations are given in Section 5.4 and finally numerical experiments are carried out in Section 5.5. 5.2 S D M e t h o d for Dev ices In this section, we consider the SD formulation for the <f-dimensional (d = 1 or 2) HD model (2.6) in the form CU = AoUt + A-VU + C(U) - V • (KVU) = 0, (5.1) where U = (p, J, T ) ' is the unknown vector of working variables. We define the problem (5.1) in a bounded domain 0 with its boundary T — I \ LHV the union of contact and insulating parts. The domain 0 is divided into three regions by junctions: the source, the channel and the drain regions (traversing Q, from left to right in Fig. 5.1 for example). Chapter 5. Simulation I: the p+-p-p+ Diode 74 IN (fin) r2 P r2 IN (fou.) Figure 5.1: Geometry of p+ — p — p+ diode. 5.2.1 Shock-Capturing SD Formulation The SD formulation for system (5.1) follows that of Chapter 3 with minor adjustments. We use the following finite element subspace {vh | Vh G {C\Qn))\ Vh \Qene (Pi(Qen))1, VVh | r ix/„(D=flr}, (5.2) where I = d + 2. The symbol V is a boundary condition operator and g is an expres-sion standing for essential boundary conditions given by (2.30)-(2.31) and (2.41). It is understood that the subspace TQ is a special case of Th with zero essential boundary conditions. As we have mentioned in Section 2.4, an additional "in-flow" boundary condition is needed mathematically. Although this was noticed and analyzed in [80], no analytic expression was given there because no obvious related physical quantity is available. We therefore seek a reasonable approximation in our cases using the momentum equation (2.1b). We assume the current and velocity are almost constant near I \ n , so the convec-tion term vV • J + J • Vt> is almost zero. This is true in many cases, for example when the doping profile is constant near r , n . If we further ignore J*, we propose the following Chapter 5. Simulation I: the p+-p-p+ Diode 75 approximation of the "in-flow" boundary condition: J • n\rin = -»o(pE + ^ V p + ^ V T ) • n | r m , (5.3) e e in which we assume that the "in-flow" current is a sum of drift, diffusion and thermal currents. This may cause an error initially, but it soon dies out when t increases. For any test function W € Th, if we denote # = (^ 2 , • • •, ^d+i)*, the space-time SD formulation with a shock-capturing modification for the HD model (5.1) is as follows. Within each Qn (n = 0 , . . . , N - 1), find Uh € Tgh such that for all tf G T0h, f [(AQU? + A-VUh + C(Uh))-y + }CVUh-Vy}dxdt (Galerkin) JQn + f A0[Uh(tn)j • *(<+) dx (Jump Condition) " e l f + J2 cuh-phmdxdt (SD) £iJQSn n-el n + J2 / v(Uh)VUh • A 0 V * dxdt (SC) + f Jh • n # • n dsdt (Weak In-Flow) ir inx/„(T) = " / A*o ( P D ^ + — V / + ^^-VTh) • n f • n dad*, (5.4) Jrinxi„(T) e ' e where A0 = diag(Ao, AQ, AQ) and the generalized gradient operator V is given by The second term to the last on the left-hand side of (5.4) is a so-called shock-capturing term which introduces a certain amount of crosswind diffusion near shocks. The coeffi-cient u(Uh) is defined as 1/(17) = 6hmax{ \*%$- , h1'2}, (5.5) Chapter 5. Simulation I: the p+-p-p+ Diode 76 where 8 is a parameter, | • \^ is an Ao-weighted /2-norm and the integer p = 1 or 2. The last terms at the left and right sides of (5.4) constitute a weak formulation of the "in-flow" boundary condition (5.3). This will alleviate the possible discrepancy introduced by (5.3) as compared with (2.1b) at nearby points (though this discrepancy should be small anyway). Moreover, it enhances the stability of the numerical scheme. 5.2.2 Concept of the Shock Capturing S D M e t h o d Although the SD method in its basic form discussed in Chapter 3 gives a dramatic improvement as compared with the standard Galerkin method, oscillations cannot be entirely prevented near the discontinuities or sharp layers of solutions. As a remedy, a shock-capturing (SC) modification is used to achieve the goals: (1) preserving high accuracy in the regions where the solutions are smooth; (2) suppressing oscillations near the discontinuities or sharp gradients of the solutions. The method is a true nonlinear one even if the problem considered is linear. A shock-capturing SD method was first proposed for the steady state convection-diffusion problem by Hughes et al. [42]. The extension to a system was studied in [41]. The method was adapted by Johnson and Szepessy [51] using the space-time finite element discretization for the nonlinear Burgers equation. There were many efforts along this line afterwards (e.g., see [77, 15, 12]). Now the shock-capturing SD method has become a standard formulation. Let us consider again the scalar convection-diffusion problem (3.2) in Section 3.2: Cu = a • Vw — KAU — / = 0. (5.6) Compared to the basic SD Petrov modification, the shock-capturing SD Petrov modifi-cation of the usual test function vh € VQ in the classical Galerkin method reads vh <— vh + ra- Vvh + T\\ <X| • Vvh, Chapter 5. Simulation I: the p+-p-p+ Diode 77 where the parameter T\\ > 0. Similar to the SD operator, we define the SC operator by P,f (vh) = T-i, a„ • V t A Different from the SD operator, the SC operator needs to satisfy only a few design con-ditions [77]: in order to control the oscillations, this operator should act in the direction of the gradient; for consistency it should be proportional to the residual Chuh; and for accuracy it should vanish quickly or remain within the order of accuracy of the basic SD scheme in smooth regions of the solutions. Q u a d r a t i c S C O p e r a t o r : For |Vwfc|2 7^ 0, the projection ay, which projects the residual Cuh onto the direction of the gradient Vu , is given by a |1 ~ \Vu*\l Clearly V u \ (5.7) o r V u * = Cuh, (5.8) arvh = 0, Vvh 6 nuU(V«A). (5.9) The inner product of the SD operator Pk{vh) and Cuh on each element fie gives the following quadratic SC form / J^%Vuh-Vvhdx. (5.10) L i n e a r SC O p e r a t o r : The SC operator satisfying the design conditions is not unique. We can also consider the following linear SC form Lj»wk7uKVvhix- (5-n) Chapter 5. Simulation I: the p+-p-p+ Diode 78 For system (5.1), we can revise the above SC operators according to [77] involving ap-pearances of the matrix AQ and the time derivative. The added h in the denominator of (5.5) is to prevent division by zero in case that IVC/71^ = 0. In (5.5) we require, according to [44], that u{Uh) > 8hz'2 (5.12) which does not degrade the order of accuracy of the scheme. 5.3 Predictor/Multi-Corrector Algorithm For each n, the discretization scheme (5.4) leads to a system of nonlinear algebraic equations. A predictor/multi-corrector algorithm (see [43, 77] for details) is used here to reduce the resulting nonlinear system to a sequence of linear systems. Moreover, we shall incorporate the potential equation solver presented in Chapter 4 within this algorithm. For simplicity, we only consider the constant-in-time finite element subspace Th here. The procedure is similar for the linear-in-time case though it is more complicated. Within each space-time slab Qn, the trial solution Uh can be expressed by uh(n+i) = g vp(x)up-*\ (5.13) p = l where n& is the number of nodal points of the space finite element mesh, Up ' is the vector of nodal unknowns at node P for the nth space-time slab and \?p is the shape function associated with node P. Define u « = ((«i"})*, (t4"})* («£?)*)*• (5.14) Substituting (5.13) into (5.4) and taking $ = $p for P = l,...,nd, (5.4) leads to a nonlinear algebraic system G(/»,u( ! ,+1)) = 0 (5.15) Chapter 5. Simulation I: the p+-p-p+ Diode 79 with uW, associated with the solution in the (n-l)st slab, coming from the jump condition in (5.4). Newton's iteration applied to (5.15) yields G(-, u{i+1)) 3 G(; « ( 0 ) + 5 y A u ( 0 = 0, (5.16) where «(t+l) = ((«l,(i+l))*, («2,(i+l))f, • • • , ("^.(H-l))*)* (5.17) is the (zH-l)st iterative approximation of u("+1) with U(0) = u ^ . Au(,-j is defined by Au(t-) = u ( t + i) - u w . (5.18) Denote by R^ = G(-,u{i)), (5.19) the residual vector. Then R^ = ((R?)\(R^)\...,(R^)\ (5.20) where i?g} = At I [(A • Vt/Jij + C( l^ ) ) ) • $ P + £ W # ) • V * P ] dx (Galerkin) + / A0 [£/$•) - tffc(n)] • # p dx (Jump Condition) + A * £ / r ^ . ^ ^ p ) ^ (SD) e = l " ^ e "e! „ +At £ / v(Ufo)VUfo • AoiUfayWp dx (SC) e = l ^ e + A t / J j ) • n ^ P • n ds (Weak In-Flow) +A* / MO ( ^ 4 ) + — V ^ ) + ^ Vlfo) • n * P • n ds (5.21) Chapter 5. Simulation I: the p+-p-p+ Diode 80 with Uk\, the zth iterative approximation to U , defined by £ p = i tffe = **>(*)«* (0- (5.22) To keep the notation simple, we denote the algebraic system from (4.16) of the Poisson solver at each time step by (tf<n\ £W) = P„,(u< (5.23) where the pair (^("),£(")) here stands for a pair of nodal unknown vectors for the pair as u^ for JJ without causing any confusion. We can now describe the following algorithm. Pred ic tor /Mul t i -Correc tor Algori thm: Given initial data, for n = 0 , 1 , . . . , N — 1, do (Predictor) «(o) = u ( n ) (ip{0),Ei0)) = Psoi(uW) (Multi-corrector loop) For i = 0 , 1 , . . . , ipass — 1, do A P A u w = - i ? ( i ) U(i+1) = U(t) + AtAu(i) (V'(»+l)»-S(i+l)) = Pso;(u(«'+i)) Continue ? .(n+l) _ , , , . . (^(»+i), £(»+*)) = (i>iipass),E(ipass)) End. Chapter 5. Simulation I: the p+-p-p+ Diode 81 In the above box, we use an explicit pass in the multi-corrector loop. The parameter ipass can be either fixed or determined by i?W. The mass lump matrix is M* = (MJK), (5.24) where MJK fQ AoVj-Vjdx HJ = K, 0 if J + K. (5.25) This method is conditionally stable. The algorithm CFL condition is of 0(At/h2) due to the thermal conductivity and the artificial diffusivity [77]. The algorithm is simple in implementation because the matrix M* is easy to invert. 5.4 Physical Parameters For the subsequent device simulations, we adopt the following system of units: Classification Length Time Mass Potential Energy Charge Temperature Capacitance Unit IO-6 m (1 fim ) IO"12 s (1 ps) 10-10 Kg Volt io-18 J io-18 C Kelvin (unsealed) IO"18 F Typical values of physical parameters for silicon devices that we use are given in Table 5.1 under such a unit system. Chapter 5. Simulation I: the p+-p-p+ Diode 82 Physical Parameters me Kb e X K0 Value 0.9109 0.138046 x 10-4 0.1602 11.7 x 8.85418 0.4 Table 5.1: Values of physical parameters. In Table 5.1, m e is the electron mass. The effective mass of the electron m = 0.26me at T0 = 300K and 0.24me at To = 77K; the intrinsic electron density pi = 1.45 x 10 l o cm - 3 at T0 = 300K and 2.84 x 10 2 0 cm' 3 at T0 = 77K. We point out that K0 in Table 5.1 is heuristically the best for the diode we consider. If KQ is away from this value, a spurious velocity overshoot will appear (see [22, 29]). The models for p,o and vs are those used in [26]: A/x where Vo P-min T 1 + (PD/Pref) Pr = 80cm2IV • s x < (300K/To)OA5 To > 200K 1.5°-45(200if/To)0-15 To < 200K, Ap = (300^/To)21430 cm Hr Nref = (To/300i;f)3-21.12 x 1017cm"3, 7 = 0.72(ro/300iO°-°65. (5.26) (5.27) (5.28) (5.29) (5.30) and 2.4 l + 0.8exp(0.5T0/300iO - x 107cm/s. (5.31) Chapter 5. Simulation I: the p+-p-p+ Diode 83 5.5 Numer ica l Simulations In this section, numerical experiments for the one- and two-dimensional p+ — p — p+ silicon diodes are performed for both the subsonic and transonic electron flows. The SD scheme (5.4) with the constant-in-time finite element subspace Th is used throughout the section. This will not degrade the entire order of accuracy of solutions since the explicit predictor/multi-corrector algorithm introduced in Section 5.3 requires that At be proportional to h2. 5.5.1 One-Dimens ional Simulations In this subsection, we use uniform meshes with h = 0.05, At = 0.2 h. The shock-capturing term in (5.4) is not included, i.e., 6 is set to be zero. Case 1: To = 300K A 0.6-pm silicon diode is used with a 0.1-pm source, OA-pm channel and 0.1-pm drain. The doping profile is given by pD(x) = 5 x 1017 cm"3 , if 0 < x < 0.05 or 0.55 < x < 0.6, pD(x) = 2 x 1015 c m - 3 , if 0.15 < x < 0.45. The above segments of pr> are connected smoothly by a scaled version of the polynomial: Q(x) = (x + l ) 4 ( - 5 z 3 + 20x2 - 29a; + 16). (5.32) Steady state solutions are obtained by letting solutions evolve in time. Solutions with biased voltages VD = 1 , 1-5 and IV are shown in Fig. 5.2. Remark 5.1 We set the biased voltage to be zero at the source as a reference. Therefore, whenever we say we apply a biased voltage VD to the device, it means that the bias is prescribed for the drain. • Chapter 5. Simulation I: the p+-p-p+ Diode 84 The Mach number \v\/c* is less than 1 for all x, so the electron flow is a subsonic one. It is worth noticing that no spurious velocity overshoot is observed using the parameters presented in Section 5.4. For comparison, numerical solutions with spurious velocity overshoots are given in Appendix C using a different set of parameters. C a s e 2: T0 = 77 K We take Gardner's example in [26]. A 1.2-pm silicon diode is used with 0.1-pm source, 1.0-pm channel and 0.1-pm drain. The doping profile is given by pD(x) = 1018 cm - 3 , if 0 < x < 0.05 or 1.15 < x < 1.2, pD(x) = 1015 cm"3 , if 0.15 < x < 1.05. The above segments of prj are connected smoothly by a scaled version of the polynomial (5.32). The biased voltage VD = IV is applied. Evolutions of velocity and temperature in time are shown in Fig. 5.3. Solutions reach the steady state at about 11 ps. We notice strong initial layers for the solutions. A steady state shock wave with Mach number 4.1 is shown in Fig. 5.4. A shock wave is developed at x ~ 0.43. It can only be formed at the transition for the electron flow from the supersonic region to the subsonic region. The cooling of electrons, when the flow enters the channel, is stronger in this case as compared with the above To = 300i^ case. The resulting velocity and temperature profiles are typical and are qualitatively the same as those in [26] and [21]. C a s e 3 : H D vs S H D To compare solutions of the full HD and SHD models, we plot the J — VD (at T0 = 300K) and J — T0 (at VD = 2V) curves of both models in Fig. 5.6 using the same 1.2-pm diode as in case 2 except pD(x) = 5 x 1017 c m - 3 , if 0 < x < 0.05 or 1.15 < x < 1.2, pD(x) = 2 x 1015 cm"3 , if 0.15 < x < 1.05. Chapter 5. Simulation I: the p+-p-p+ Diode 85 6-5-4 -3-2-1-0 ' ( 0.1 DENSITY Legend 1.0V0II 1.5 Volt in von 0.2 0.3 j 1 Ii —±=-^ r 0.4 0.5 '. 0.6 1542-1028-514 -Legend 1.0 Volt 1.5Volt 2.0 Volt ) 0.1 TEMPERATURE A 02 0.3 0.4 0 5 0 6 0 --2.5-•5.0--7.5--10.0 -y\ Legend 1.0 Volt 1.5 Volt *** "™ 0.1 ELECTRICAL FIELD '^>\ l 0.2 0.3 0.4 0.5 0 6 Figure 5.2: HD model at 300K: Electron densities in 1017cra 3; velocities in 108cm/s; temperatures in Kelvin and electric fields in 104V/cm. Chapter 5. Simulation I: the p+-p-p+ Diode 86 Evolution in Time: Velocity Profile Evolution in Time: Temperature Profile 12 1.2 Figure 5.3: HD model at 77K: evolutions of selected solutions. Chapter 5. Simulation I: the p+-p-p+ Diode 87 4-3 2-1-VELOCITY AND MACHNO , ) 0.1 ,' ,7 0.2 ,'" ' 0.3 \ ^ | •»-. . 0.4 0.5 0.6 0.7 0.8 Legend Velocity MachNo ~~",l 0.9 1.0 1.1 1.2 1250 -1000-750-500 -250 -) 0.1 0.2 0.3 TEMPERATURE 0.4 0.5 0.6 0.7 0.8 0.9 L_ 1.0 1.1 1 2 0--2--3-) 0.1 0.2 0.3 ELECTRICAL FIELD 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r .1 1 2 Figure 5.4: HD model at 77K: Electron density in 10 i r cm - 3 ; velocity in I07cm/s and Mach number; electron temperature in Kelvin; electric field in 104V/cm (from left to right and top to bottom). Chapter 5. Simulation I: the p+-p-p+ Diode 88 VELOCITY 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1250-1000-750-500-250-0 i 0.1 0.2 0.3 0.4 TEMPERATURE 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1 2 Figure 5.5: SHD model at 77K: Electron density in 1017cm~3; velocity in 108cm/s ; electron temperature in Kelvin; electric field in 104V/cm (from left to right and top to bottom). Chapter 5. Simulation I: the p+-p-p+ Diode 89 J - V CURVES 1.0 • Legend SHD Model ^ ^ ^ 1.5 2.0 Biased Voltage 2.5 3.0 J - T CURVES 100 • v- „ v Legend SHD Model ^ " " ^ ^ ^ ^ . ^ ^ 150 200 Lattice Temperature 250 300 Figure 5.6: Comparisons of J — VD and J — T0 curves using the HD and SHD models. Chapter 5. Simulation I: the p+-p-p+ Diode 90 The difference is not very pronounced in the J — VD curves when the biased voltage Vb is increased. However, the difference is more pronounced in the J — To curves when the ambient temperature To reaches 77K. We have observed that the current density of the full HD model is always smaller than that of the SHD model. The steady state solutions of the SHD model at T0 = 77K, which are noticeably different from those of the HD model, are presented in Fig. 5.5 with exactly the same parameter and geometry settings as in case 2. 5.5.2 Two-dimens ional Simulations Let us consider the HD model of the two-dimensional p+-p-p+ silicon diode shown in Fig. 5.1 with ft = [0, l(fim)] x [0,0.5(fim)]. The two junctions are given by *i(y) = \-\(y-m2, (5-33) *M = l + \y2- (5-34) We define the doping profile (see Fig. 5.7): PD{X) = 5 x 1017 cm - 3 , if 0 < x < x1 or x2 < x < 1.0, pD(x) = 2 x 1015 cm - 3 , if xi + 0.1 < x < x2 - 0.1. The above segments of pp are, as in the one-dimensional case, connected by a scaled version of <5(0 of (5.32), where £ = 200.0 x (xi(y) — x + 0.005) for source-channel junction £ = 200.0 x (x — £2(2/) + 0.005) for channel-drain junction An 80 x 40 uniform mesh is used in the following computations for the steady state solutions for the subsonic flow at To = 300K and the transonic flow at T0 = 77K. The bias is Vb = IV in the following cases. Chapter 5. Simulation I: the p+-p-p+ Diode Figure 5.7: Doping profile in 1017 cm 3 Case 1: To = 300/^ In this case, we set At = 0.002 and 6 = 0. We present the electrical field and streamlines of the current in Fig. 5.8 and Fig. 5.9. Surface plots for the electron density and temperature are shown in Fig. 5.10. The Mach number in this case is 0.7. We have observed oscillation of the current near junctions due to the sharp profile of the doping pr>. This (numerical) phenomenon was also reported in [22]. Since the current is small in comparison with other solutions, it is very sensitive to any small perturbations. Therefore it bears a relatively larger error than other quantities. The results show that the distortion of the current is more prominent at the p-p+ junction than at the p+-p junction. The error in the current is proportional to the order of accuracy of the numerical scheme, so higher order schemes or alternatively finer grids could help. Case 2: T0 = 77 K In this case, we set At = 0.001 and 6 = 1. A shock wave is presented with Mach Chapter 5. Simulation I: the p+-p-p+ Diode 92 number 5.2 at the steady state (see Fig. 5.11). The shock profile develops in the channel and spreads over about 3 grid points. The peak value of Mach number can reach as high as 7.8 in evolution. The cooling effect is milder than that in the 300K case. Numerical computations were applied to different geometries of the diode as well as different doping profiles. For some cases with shocks or abrupt junctions with sharp corners, we need to add the shock-capturing term defined in Section 5.3 with 8 > 0, which stabilizes the numerical scheme. Note, though, that we want to take 8 as small as possible: if 8 is large then it causes an overly diffused behavior of the solutions and, more seriously, it produces a poor current approximation as explained in the previous paragraph. The size of 8 used here is still much smaller than what is required if the classical artificial viscosity method is employed. 5.5.3 C o m m e n t s and Conc lus ions Shocks of HD models are simulated for both one- and two-dimensional p+ — p — p+ diodes by the shock-capturing SD method. Robustness of the employed SD method is evidenced by its performance and will be demonstrated again in the MESFET simulations in Chapter 6. Although in many applications, such as the one in the following chapter, the solutions of the SHD model are close to those of the HD model, we observed an exceptional case. The prominent difference between these two models is that the HD model, but not the SHD model, supports shocks. With the same parameter setting, the behavior of the HD model is substantially different from that of the SHD model when shocks are developed in the HD model. Numerical experiments presented here have shown the influence of the shock on the solutions of the HD model. Failure to capture the shocks might cause serious error in the source-drain current. We also examined different components of the current and found that the strength of Chapter 5. Simulation I: the p+-p-p+ Diode 93 the thermal current is less than that of the drift current or diffusion current. Chapter 5. Simulation I: the p+-p-p+ Diode 94 Figure 5.8: HD model at 300K: Electrical field. Figure 5.9: HD model at 300K: Streamlines of the current. Chapter 5. Simulation I: the p+-p-p+ Diode 95 Figure 5.10: HD model at 300K: Electron density (left) in 1017 cm 3 and electron tem-perature (right) in Kelvin. Figure 5.11: HD model at 77K: Electron density (left) in 1017 cm 3 (image is inverted) and electron temperature (right) in Kelvin. Chapter 6 Simulation II: the M E S F E T Device 6.1 Introduct ion The metal semiconductor field effect transistor (MESFET) is one of the important semi-conductor devices. The physical characteristics of the device are manipulated by a rec-tifying metal-semiconductor contact at the the gate, known as a Schottky contact after Walter Schottky, who first proposed a model for the potential barrier formation in 1930s. Although much progress has been made in the development of metallic contacts in semiconductor technology, a full understanding of the mechanism of barrier formation is far from complete. Therefore, a highly simplified model for simulation is used. The MESFET devices were simulated by Cook and Frey [17], Jerome and Shu [51] and Chen et al. [13] using the hydrodynamic (HD) model. Although there were no shocks observed, the shock-capturing algorithms in [51] and [13] proved to be useful when sharp layers occur. In this chapter, MESFET devices at room temperature are simulated effectively by our streamline-diffusion (SD) finite element methods using the HD models as well as the drift-diffusion (DD) models. We indicate that the SD method (without the shock-capturing term) for the conventional DD model is identical to the well-known Scharfetter-Gummel (SG) discretization [71] in one space dimension. Consequently, SD methods can be viewed as reasonable extensions of the SG formulation for the general convection-diffusion systems in multi-space dimensions. 96 Chapter 6. Simulation II: the MESFET Device 97 Vs = 0 o p+ \ \ \ Depletion Region \ p P+ ' J V # *-Figure 6.1: Schematic of MESFET operation. Numerical experiments are implemented for models in an increasing order of compli-cation: the conventional DD (DDl) model —> the mixed-formulation DD (DD2) model —•»• the simplified HD (SHD) model —• the full HD model. Comparisons are made among solutions of the different models. We find a qualitative discrepancy between the solutions of the DD models and the HD models. A small difference, mainly in the temperature, is also found between the full HD and simplified HD models. The arrangement of the chapter is as follows: In section 2, the physical aspect of the MESFET device (or essentially that of the Schottky diode), i.e., the formation of the Schottky barrier, is modeled. In section 3, the methodology of the SG discretization and close relationship with the SD concept are described. Numerical results and discussions are supplied in section 4. Good performance of the numerical method is demonstrated in this simulation. Chapter 6. Simulation II: the MESFET Device 98 e Gate \ i EFm vletal ^ ^ *— d —» (a) t 1 Semiconductor -E c - E F E v e 1 i <t>B - wd T e,fso - \ TZTi (b) y * \' e \|. 1 .-_eJto.. t Figure 6.2: Formation of Schottky barrier 6.2 Schot tky Barrier A typical MESFET structure and its operation are shown in Fig. 6.1. The source, gate and drain regions are prescribed with external biased voltages Vs, VG and VD, respectively. We assume that the source is grounded, i.e., Vs = 0. The main feature of the device is the formation of the so-called Schottky barrier (of the potential) at the interface of the gate by intimately contacting the metal to the semiconductor. To derive the mathematical model reflecting the formation of the barrier, it is sufficient to consider the thermal equilibrium state of the device. 6.2.1 Formation of Schot tky Barrier The commonly employed energy band diagrams can be used to describe the physical characteristics of the device. Fig. 6.2 illustrates the process of barrier formation near the Chapter 6. Simulation II: the MESFET Device 99 gate (along the y direction) in thermal equilibrium. To depict this process, we introduce a metal working function <f)m with e<f>m being defined as the amount of energy required to raise an electron from the Fermi level Ep to the vacuum level. Vacuum level is the energy level of an electron just outside the metal with zero kinetic energy. Similarly, we can introduce a semiconductor working function <f)s with <f>s < (f>m. The working functions describe material differences. For a semiconductor, since <j>s (actually the Fermi level Ep) varies with the doping, therefore another parameter Xs, the electron affinity, is employed instead. By the electron affinity, exs is defined as the energy difference of an electron between the vacuum level and the lower edge of the conductor band Ec-According to the Schottky model, the barrier results from the difference in working functions of the two materials. When metal and semiconductor get closer than a very short distance d (see energy diagram (a) in Fig. 6.2), an electric field, as shown in Fig. 6.3, exists in the region between the metal and semiconductor due to the higher energy of electrons in the semiconductor side. The electrons flow into the metal until the Fermi levels on the two sides are brought into coincidence. As the electrons move out of the semiconductor into the metal, the free electron density in the semiconductor region near the boundary decreases. Note that in thermal equilibrium the Fermi level remains constant. Since the separation between conductor band edge Ec and the Fermi level Ep increases with decreasing electron density, the conductor band Ec and valence band Ey are bent upward. Typical electron density and potential barrier profiles near the gate are calculated in Fig. 6.4. Energy diagram (b) in Fig. 6.2 shows the case after the full contact is made and thermal equilibrium is reached. To measure how much the energy bands bend up, we define </>#, a quantity called the barrier height function. It is clear that <j>B = <t>m- Xs- (6-1) Chapter 6. Simulation II: the MESFET Device 100 E l e c t r i c F i e l d ft ft / / / I I I I I I I \ \ \ / / / I I I I I I M \ % . i / i i i i i i i i i , Figure 6.3: Electric field near the Schottky contact. The region where the energy bands bend is called the depletion region since the electron density there is less than the doping density. Wd denotes the width of the region. From the thermionic-emission theory, only a very small leakage current exists in the gate region if a reverse bias voltage is applied to the gate (see [79, 66]). In this work, we are interested in the case where <f>s > 0. 6.2.2 Mathemat i ca l Mode l In thermal equilibrium, the external biased voltages are set to be zero (VQ = 0 and VD = 0). In such a case, J = 0 and T = T0. Then all the device models (2.1)-(2.5) in Chapter 2 reduce to the steady state as: T0Vp + epE 0, V ( A - V 0 ) = /o -PD. (6.2a) (6.2b) Solutions of (6.2) are presented in Fig. 6.5 with <f>B = 0.8V and pr> = 10 cm . Chapter 6. Simulation II: the MESFET Device 101 ELECTRON DEPLETION NEAR SCHOTTKY CONTACT 0.20 POTENTIAL BARRIER NEAR SCHOTTKY CONTACT 0.20 Figure 6.4: p and — ip with (f>B = 0.8V and po = 10 cm . Schottky contact is at y=0.2 Chapter 6. Simulation II: the MESFET Device 102 Recall the quasi-Fermi level <f>n, which is a constant in the thermal equilibrium case. Integrating (6.2a) yields ^ ^ - E F / e = V - — I n A (6.3) e pi Equation (6.3) provides a potential-electron density relation ^ = ^ l n A (6.4) e Pi if we set <j>n = 0 as a reference level. For the Schottky contact, the electron density can be expressed by p0 = Ncexp( —<^B) , (6.5) ml0 defined as the surface electron density (see [72, 79]). Nc is termed the effective density of state in the conduction band, which has the typical value of Nc = 4.07 x 101 9cm - 3 for silicon at room temperature. The above admits an alternative expression p0 = PD exp( — ( ^ B - <f>o)), (6.5') mlo where ecf>o measures the energy difference between the conduction band edge and Fermi level, and raTo Nc (R r., (f>0= In . (6.6) e PD Substituting p0 in (6.4) by (6.5) or (6.5'), we obtain the built-in potential at the Schottky contact of the gate tf.^mife (6.7) 6 Pi or e pi Chapter 6. Simulation II: the MESFET Device 103 In the source and drain, we have Ohmic contacts and p = pD. Therefore V»g = ^ l n ^ . (6.8) e Pi Now the boundary conditions defined in Chapter 2 can be completed. ipso in (2.37) is given by ipso = <f>B- fa- (6.9) In order that the semiconductor operates in the non-degenerated parameter domain (i.e., Ec — EF > 2raT0), we require pD < exp(-2)Nc = 5.508 x 1018cnT3 . (6.10) R e m a r k 6.1 The boundary conditions we employed here are slightly different from those given in [51] and [13], where all variables are required to satisfy Neumann conditions on T2. Those conditions are obviously over-determined and may cause formation of artificial boundary layers. • 6.3 Discret izat ions for D D Mode l s It is convenient to express the DD models (2.4) and (2.5) in the uniform form CU = Ut + A-VU + C{U)-V-()CVU) = 0 (6.11) as equation (2.6) for the HD models, where U = (p, J)* for the DD2 model and U — p for the DD1 model. For two-dimensional problems for example, A = (Ai, A2)*, C(U) and )C are given by, for DD2 Model: Chapter 6. Simulation II: the MESFET Device Electron Density Potential Figure 6.5: Electron density (left) and potential (right) with <J>B PD = 1017cm3 in thermal equilibrium. 0.8V Ax ( 0 1 0 X To 0 0 ^ 0 0 0 ) A2 = ' 0 0 1 N 0 0 0 V T0 0 0 J ( c 0 :PE1 + ;PE2 + \ TJ and K = 0; for DD1 model: M = -fJ-opE1, A2 —p0pE2, C = —pLopV • E and K, = D0I2. Chapter 6. Simulation II: the MESFET Device 105 6.3.1 T h e S D M e t h o d for D D Models As in Section 5.2.1 for the HD models, we use the the following finite element subspace for the DD models Th = {Vk | yk e (c°{Qn)y, Vh \Qene (P1(Qen))\ Wh |r lX/„(T)= 9h (6-12) where I — d + 1 for the DD2 model and / = 1 for the DDl model. Here g stands for essential boundary conditions given by (2.34) and (2.36) for the DDl model, and in addition (2.31) for the DD2 model. The DD models can be solved using the same scheme as (5.4) except that we drop the (SC) and (Weak In-Flow) terms, which reads / (Uth + A • VUh + C{Uh)) • * + KVUh • V $ dxdt (Galerkin) + I lUh{tn)\ • V(t+) dx (Jump Condition) J \l + J2 £Uh-Ph(V)dxdt = 0. (SD) (6.13) e = l J Q n 6.3.2 A S D M e t h o d for the Steady State D D l Mode l In the case when ph is taken to be a piecewise linear polynomial approximation to p, the steady state SD scheme for the scalar DDl model reduces to (V • (nophEh), tt) + ( A ) V / , W ) (Galerkin) + ( V • (nophEh), Tn0Eh • V * ) = 0. (SD) (6.14) We see from the above that the streamline is along the direction of the electric field. For-mulation (6.14) is similar to the multi-dimensional SD extension of Scharfetter-Gummel discretization by Sharma and Carey in [78], but the scaling parameter r we use here is different from theirs. In [78] the numerical electric field Eh is assumed to be piecewise Chapter 6. Simulation II: the MESFET Device 106 constant, therefore the (SD) part in (6.14) can be simplified to (V • (fi0phEh), rfi0Eh • V * ) = (fx0Eh • VP\ r/i0Eh • V * ) . (6.15) However, we use a piecewise linear approach to Eh here. 6.3.3 SG M e t h o d for the 1-D Steady State D D 1 Mode l The steady state DD1 model in one space dimension d4- = 0 (6.16) dx yields a constant current J on the whole space domain 0 , where dip dp . J = , 0 p - - D 0 - . (6.17) Employing a uniform partition of ft with grid points Xj = x0 + jh (j = 1 , . . . , iV), we discretize (6.16) and (6.17) by assuming that within each interval (XJ, Xj+i), the potential •0 is linear and the electric field E = —dip/dx is constant. Integrating (6.17) on j t h interval to obtain p explicitly as p(x) = Pj exp(-^(</, - 0,-)) - —f~ [1 - e x p ( ^ ( ^ - ^ ) ) ] , (6.18) where pj and ipj are approximate solutions at Xj, and £ i+ i / 2 = (ipj+i — i/>j)/h. Setting p(xj+i) = pj+i on j t h interval, (6.18) implies J = -D 2 ^ + 1 / 2 \pj+1 e x P ( ^ + 1 / 2 ) ~ Pi e x P(~^ '+ i /2) i / 6 19x h exp(^-+ 1 / 2) - e x p ( - ^ i + 1 / 2 ) where Pj+1/2 = l^j+i^l and a j + 1 / 2 is the cell (element) Peclet number defined by (3.8). With the DD1 model in this case, aj+1/2 = ^-hEj+1/2. (6.20) Chapter 6. Simulation II: the MESFET Device 107 Using (6.19) on j-'th and (j+l)st intervals, the constant J yields the SG discretization Pi+i[&+i/2(coth(/? i+1 /2) + 1)] - pj[/3j+1/2(coth{/3j+1/2) - 1) # _ 1 / a ( c o t h ( # + 1 / 2 ) + 1)] + ^•_1[/3 i_1 / 2[coth(^_1 / 2) - 1)] = 0. (6.21) The above procedure is also known as an exponential fitting. As for the one-dimensional SD method (6.14), the scaling parameter r that follows (3.11) is r = — — [coth(a i + 1 / 2) - aJ*1/2] (6.22) Zfj,0Ej+1/2 y T ' on the jih interval. We note that a i + i / 2 [ c o t h ( a i + 1 / 2 ) - a7+1/2] = /3 j + i / 2 [coth(^+ 1 / 2 ) - ^ f 1 / 2 ] . (6.23) Therefore, the (SD) term in (6.14) can be written as (HoEh • V / , Tfi0Eh • V * ) = (D,Vph, V * ) (6.24) with artificial diffusivity Ds given by D, = D0pj+1/2[coih{PJ+1/2) - / 3 ^ 1 / 3 ] , (6.25) if we use a piecewise linear iph and a piecewise constant Eh. Following the arguments in [78], it is not difficult to show that (6.21) can be recovered by (6.14) with the streamline diffusion (6.24). Therefore the SD methods are reasonable extensions of the SG method in the context of numerical solutions for semiconductor devices. 6.4 Numerica l Simulat ions We simulate two-dimensional silicon MESFET devices at room temperature To = 300K by using the HD, SHD, DD2 and DD1 semiconductor device models given in Chapter 2. Chapter 6. Simulation II: the MESFET Device 108 The HD and SHD models are solved by the shock-capturing SD method (5.4), but the finite element subspace Tgh defined in (5.2) is changed to one with the function g rep-resenting essential boundary conditions (2.30)-(2.31), (2.34) and (2.36) corresponding to the MESFET. The DD2 and DDl models are solved by the SD method (6.13). The finite element subspaces we use are the constant-in-time, linear-in-space types. Numerical simulations are carried out on a device domain fl = [0,0.6/jm] x [0,0.2pm]. The source, gate and drain regions are given by segments of [0,0.1] x {y = 0.2}, [0.2,0.4] x {y = 0.2} and [0.5,0.6] x {y = 0.2}, respectively. The problem that we are interested in is the one with the gate in the state of reversed bias (i.e., <f>B > 0). In such a case, a large depletion region occurs near the gate, such that the carrier movement in this region is negligible due to the formation of the potential barrier. Therefore electrons can flow from the source to drain only through a channel underneath the gate, if an external biased voltage is applied to the drain (see Fig. 6.1 for the depiction of the device operation). Unless otherwise specified, the potential barrier height function is (f>B = 0.8V and a biased voltage VD = 2V is applied to the drain in the following computations. Initial values are solutions of thermal equilibrium cases. 6.4.1 E x a m p l e 1 The doping profile pD is given by (see Fig. 6.6): po = 3 x 10 1 7em - 3 in [0,0.1] x [0.15,0.2] and in [0.5,0.6] x [0.15,0.2], and pD = I x 101 7cra - 3 elsewhere, with abrupt junctions. This was used in [51, 13] Comparisons of numerical solutions of different models are made in this example. Steady state solutions are obtained by letting the solutions evolve in time. A uniform grid of 96 x 32 points is used first. The time step A t= l .D-3 . Results with this grid are shown in Fig. 6.8 to Fig. 6.13. The surface electron density is about po = 1 .567xl0 6 cm - 3 in this case. The surface plots of electron densities of the four models Chapter 6. Simulation II: the MESFET Device 109 Doping profile 2-D MESFET Device Figure 6.6: Doping and 2-D geometry of the device for example 1. are shown in Fig. 6.8. We find that solutions of HD models are different from those of DD models. The HD and SHD models have boundary layers for p near the drain, whereas the DD2 and DD1 models do not exhibit such layers. This is a noticeable phenomenon of the hot electron effect in the hydrodynamic models. Mathematical analysis of the possible formation of a boundary layer near the outlet boundary of the device for the HD model was given by Gamba (see [24, 25]). Also, there are valleys in the solutions of the HD models, which stretch from the gate region through the substrate. The electron current is an important component in device simulation. In Fig. 6.9, we present the current flows of models DD2, SHD and HD. It clearly shows the flow of electrons from the source to the drain. In our simulation, the leakage current of electrons is correctly reflected (from the semiconductor towards the gate, which is shown by the plot of the velocity component in y direction in Fig. 6.10), although it is very small. This is noticeably different from that of [13]. The width of the depletion region is about Wd = 0.1/xm. We have observed that Wd for the HD models is smaller than that for the DD models. We detect an error in the current around the junction near the source, but not near the drain. This might be caused by numerical dissipation, which is mainly Chapter 6. Simulation II: the MESFET Device 110 manipulated by the shock-capturing term. Since the current is very sensitive, it is easily polluted by any perturbation. The parameter 6 = 1 is used in the example. We also experimented with different values of the parameter 6. When 6 decreases, we get a better current approximation. In order to control the possible oscillation of solutions and achieve stability of the numerical computation, we pay the small price of slightly losing accuracy in the sharp change regions. However, the approximated electron current obtained does not exhibit non-physical behavior, i.e., there is no recirculation of the electron flow (cf. [13]). We provide temperature profiles for the HD and SHD models in Fig. 6.11 and poten-tials for the HD and DDl models in Fig. 6.12. There is a bigger cusp in the temperature of the SHD model around the point x = 0.2, y = 0.2 than in the HD model. In Fig. 6.13, we show solutions of the HD model when <f>B — cf)0 = 0 .8^ . Correspond-ingly, the surface electron density is about po = 1.81 x 10 4cm - 3 . The solutions do not seem to have big differences to those in the <f>jg = 0.8V case. We make the following observations: (i) The solutions for the HD model and SHD model are close except for the temperatures. In the device simulation literature, engineers often prefer to use the SHD model, because it approximates the full HD model very well in many realistic situations and takes into account most important physical phenomena missing in the DD model. Such is the case in the simulation reported here, (ii) The SD method (6.13) for the DDl model with r defined in Chapter 3 (1=1), is optimal according to [42]. Comparisons were made by Sharma and Carey [78] for a diode and a MOSFET. Their evidence shows that the SD method causes less smearing than other SG extensions. The SD method for the DD2 model is reasonable although perhaps not optimal unless the system can be decoupled. The numerical experiments show that the results obtained for the DD2 model are no worse than those for the DDl model. In Fig. 6.14, cross-sectional electron densities at y=1.75//ra of DDl vs DD2 models and HD vs SHD models are given Chapter 6. Simulation II: the MESFET Device 111 for the purpose of comparison. A finer mesh of 192 x 64 points is used next. The time step At =5.D-4 and parameter 8 = 1.5 for this grid. The steady state electron density and current flow of the HD model are given in Fig. 6.15. The electron current is slightly improved after the refinement, while the improvement of the electron density is almost invisible. 6.4.2 E x a m p l e 2 The doping p r o f i l e ^ is given by (see Fig. 6.7): po = 3 x l 0 1 7 c m - 3 in (x —0.2)2+y2 < 0.12 and in (x — 0.2)2 + (y — 0.6)2 < 0.12 , and pD = 1 x 101 7cra - 3 elsewhere. Now abrupt junctions are not aligned with mesh lines. Uniform grids of 96 x 32 and of 192 x 64 are used again for this example. Simulated results for the steady state solutions of the HD model are shown in Fig. 6.16 and Fig. 6.17. The numerical solutions behave qualitatively very well near the junctions. In fact, the SD method is so designed that it deals with this situation effectively. We notice in this case that the boundary layer near the drain for the electron density is very mild, in contrast to that in example 1. This suggests that the formation of the boundary layers depends strongly on the geometric setting of the device. A similar phenomenon was simulated and analyzed using a 1-D current driven problem by Ascher et al. in [2]. Evolutions of the solutions are shown in Figs. 6.18- 6.21 using the 96 grid. We increase the biased voltage from Vp = 0 in the thermal equilibrium state, to VD = 2V within O.lps and then stay at this value for the rest of the time. It is found that solutions behave wildly in a short time interval right after t=0.1ps. We use 8 = 1.5 during this time interval in order to enhance the stability. We reduce it to 1 when the solutions are close to the steady state. The boundary layer of the electron density and the strength of the current density at the outlet become milder when the steady state is reached. Chapter 6. Simulation II: the MESFET Device 112 Doping profile 2-D MESFET Device Source V ) \° I Gate f« a2^ J P Substrate Drain + \ P 0.6um r 1 1 Figure 6.7: Doping and 2-D geometry of the device for example 2. 6.5 C o m m e n t s and Conclusions We have applied the SD method to 2-D MESFET devices by using the HD and DD models. The method extends the well-known SG formulation. The performance of the presented numerical schemes is very satisfactory and the computational results meet physical expectations well. Comparisons of different models are made in the numerical simulations. They indicate that the results of the HD and DD models are qualitatively different. However, the results of the SHD and full HD models are relatively close in this simulation, especially for the electron density and current. Chapter 6. Simulation II: the MESFET Device 113 Electron Density of DD1 Model Electron Density of SHD Model Electron Density of DD2 Model Electron Density of HD Model Figure 6.8: Electron densities in 1017cra 3 for example 1. Chapter 6. Simulation II: the MESFET Device E l e c t r o n Current of DD2 M o d e l M \ \N \ NS ' T T T T N S \ \ v l \ \ \ \ \ \ \ \ \ ^ v -I V \ \ \ \ \ \ N N"s >» >• -rr •j . / M l ,'jJJAA. . , , iIII , • / / / / / / / / / / / / •S/S/ / I I •^^y / / / i ---~.--~ ^ / / / E l e c t r o n Current of SHE) M o d e l 7T I l I i \j v s ' ' • • i •» \ \ M \ \» • • • r-r-\—\—\' \ \ \ \ - ' \ \ \ \ \ \ \ \ "> s \ \ \ \ N \ "^ ^ ---/ ,ILJ.L , , I I s , / / / s / / 11 ,s s / I I ^ ^ / / / I E l e c t r o n Current of HD M o d e l 77T -v"v\ \ \ \ \ • • ' \ \ \ \ \ \ \ N N ' \ \ \ \ \ \ \ N ~ -\ \ \ \ N "s X ^ ~-~ • •ULJ.L , , / 1 I , , / / 1 1 s s , / I I ^^ s S / I I • ^^^ s / / I Figure 6.9: Electron currents for example 1. Chapter 6. Simulation II: the MESFET Device 115 Velocity Component Vx 0 . 2 x 0 0 V e l o c i t y Component Vy Figure 6.10: Mean velocity components for example 1. Chapter 6. Simulation II: the MESFET Device Electron Temperature of SHD Model Figure 6.11: Electron Electric Potential of DD2 Model Electron Temperature of HD Model ure in Kelvin for example 1. Electric Potential of HD Model Figure 6.12: Electric Potential for example 1. Chapter 6. Simulation II: the MESFET Device 117 Electron Density E l e c t r o n Ouarirerit I I \ \ \ j \ \ ' ' • ' " T \ - \ ' \ \ V > • • \ \ \ \ \ \ \ - ' \ \ \ \ \ \ s x • . i / l I I ,'u-J.L . , , / / , , / 1 I . , , , / ! I •^s, , t I ^ ^ s / I Electron Temperature Electric Potential Figure 6.13: Solutions for 4>B — <f>o = 0.8V. Chapter 6. Simulation II: the MESFET Device 118 DD1 vs DD2 0.2 0.3 0.4 0.5 0.6 HD vs S H D Figure 6.14: Comparison of electron densities at y — 0.175//ra for example 1. Chapter 6. Simulation II: the MESFET Device 119 Electron Density Electron Current I I I \ \ T T T T N \ \ \ \ • ' i \ \ \ \ \ \ \ \ •> ' l l \ \ \ \ \ \ \ V /u.J.1. ,,/1 • , / / l I - , , , ! I \ Figure 6.15: Electron density and current with the finer mesh for example 1. Electron Density Electron Current l l l l , ! ^ 1 " " i i \ y\ \ \ \ * * • i i \A \\\w- - ' . . , 1 / / . - . \ / / , , NJ / / / / /1 Figure 6.16: Electron density and current for example 2: a grid of 96 x 32 points. Chapter 6. Simulation II: the MESFET Device 120 Electron Density a n t D e n s i t y I I 1 I ,\ n \ \ A N N v' I i \ y'\ \ \ \ * • • I \ \A \ \ \ \ N N ' j..v'\ \ w w w l \ \ \ \ \ \ N W I . . I /1 , , \ / / , , , ru . , , , — ^ x / / Figure 6.17: Electron density and current for example 2: a grid of 192 x 64 points. Chapter 6. Simulation II: the MESFET Device 121 Figure 6.18: Evolution of the electron density for example 2. •g ajdurexs JOJ ^traurto uoipap 9q^ jo uoi'mjoAg :6I*9 aanSt^ 4J'"" \"H"" 11V'' ' / ' V ' » ^ ^ N S N \ \ \ \ \ , W S W W W \ , s W W \ \ \ y i " . . . v \ \ \ \,\ \ II . . . , n \ ) ' U I I s i l l I I s<3 0 - E = 3 he"-\lh<-, ^ - * . W N \ \ \ \ \ \ \ v . x N W \ \ \ \ \ I \ w W W W \ \,<'~\ . > \ \ w \ \ v\ \ I . . . , v \ \ \ \A \ I I * w)' \ n i v i l l i I 3d 0 ' T = 3 I / - - • -/ / / / " • Hu\" , ^ ^ N N \ \ \ \ \ \ A \ W W W \ V . w w \ \ \ \,x'r , s W W \ y'\ I I . . x \ \ \ \/\ \ \ \ ..,, x \ y \ i 11 s<3 9 - 0 = = l \ S N N \ S \ \ V \ W W N \ \ \ \ \ I ^ N W W W W W W \ \ i i v \ \ \ w \ v\ I I i w w / h n . ,x\ y u i i . . w i l t s d S2 -0=3 umTjqTXTnbs xE U Z J 3 x ;Ii sa T-O=3 221 33.rA9(j XHdSHH dm -II uoTivjnuiTS g ^dvifj Chapter 6. Simulation II: the MESFET Device 123 Figure 6.20: Evolution of the electron temperature for example 2. Chapter 6. Simulation II: the MESFET Device 124 o. WSSBK If I P ^«^ A A Figure 6.21: Evolution of the potential for example 2. Chapter 7 Conclusion 7.1 S u m m a r y Theoretical and practical aspects of the design and implementation of the streamline-diffusion (SD) finite element methods to the solutions of the semiconductor device models have been explored systematically in this thesis. The two major device models, i.e., the fundamental drift-diffusion (DD) model and the hydrodynamic (HD) model, which are widely used today, were considered. Emphasis, however, was placed on the HD model, which is computationally more challenging in the context of device simulations than the traditional DD model, but provides some important physical information missing in the DD model [68, 57, 76]. Under certain assumptions, the simplified HD (SHD) model, the mixed-formulation DD (DD2) model and the conventional DD (DD1) model have been deduced hierarchi-cally from the full HD model with gradually reduced orders of complexity. These models were solved by our non-symmetric SD method in a uniform framework. The non-symmetric SD method for the convection-diffusion systems was designed to avoid the extra complexity of symmetric formulations of the systems involved. We derived a non-trivial SD operator for the non-symmetric system. The designing work includes two concerns: one is to follow the correct directions of the streamlines and another is to provide suitable scaling. The proper choice of the operator is important - a careless choice of it leads to serious problems, as was shown in this work. With 125 Chapter 7. Conclusion 126 an exploration of useful properties of the semiconductor device models, we derived the SD operator for these models. Linear stability analysis was carried out for the model problem. We have shown that our proposed numerical scheme is stable if the system can be symmetrized. Stability arguments and numerical experiments presented suggest that the method of lines with SD semi-discretization schemes may not be appropriate for transient problems. Such combinations were used previously by some researchers. The SD method developed was applied to the models of semiconductor device trans-port equations. From the mathematical point of view, an "in-flow" boundary condition is needed for the HD model unlike the other models [80]. The weak treatment of this extra "in-flow" boundary condition was incorporated in the SD formulation for the HD model. This practical consideration alleviates the possible discrepancy introduced by the intuitive artificial boundary condition and enhances the stability of the numerical scheme. Pertaining to Poisson equation for the potential, an efficient method consistent with the SD method is desired. It should provide an accurate electric field, i.e., an accurate derivative of the potential as well. For the HD model, for example, not much attention was previously paid in the method design regarding the consistency of schemes for the conservation laws and the potential equation. We proposed a consistent Poisson solver which provides a higher order convergence rate for the electric field. The order of accuracy of the electric field is the same as that of solutions for the system of conservation laws when employing finite element subspaces with the same degree of polynomials. Different from the commonly used mixed method that uses a staggered grid, this Poisson solver uses the same grid points for the potential and electric field, which are also used by the SD scheme for the system. Therefore the device simulation algorithm can be more easily coded. The new Poisson solver absorbs the idea of the augmented and least-squares methods Chapter 7. Conclusion 127 [10, 52, 53], but it is superior to them. Its decoupled nature saves computation time as compared with the mixed method, the augmented method and the usual least-squares method. Error estimation was presented. The convergence result for such a method is better than the usual least-squares method [63] and the augmented method. The p+ — p — p+ diode and the MESFET device were simulated by using our numerical methods. Shocks were captured for the p+ — p — p+ diode in one and two space dimensions. For the two-dimensional MESFET device, electron depletion near the gate was simulated. The leakage current of electrons was correctly reflected, although it is very small and can be neglected. The method successfully handles the case where the abrupt junctions are not aligned with grid lines. Experiments show that the performance of the method is very satisfactory and the computational results confirm the physical expectations. Model comparisons were implemented. The difference in solutions between the HD and DD models was significant in simulation of the MESFET device. However, the solution difference between the SHD and HD models was negligible in this simulation. In the engineering field, the SHD model and other simplifications of the HD model appearing in the literature (such as the so-called energy balance (EB) model [7, 8]) have often been viewed as good approximations of the full HD model and have been comfortably used in replacement of the full HD model in many applications. However, as was demonstrated in our numerical experiments, if a strong field is applied or shock waves develop in the HD model, the solutions of the SHD model will no longer approximate those of the HD model well. Comparison of steady state numerical solutions between DD1 and DD2 models con-firms that the performance of our non-symmetric extension of the SD operator for the convection-diffusion system is almost as good as that of the optimal one given in [42] for the scalar convection-diffusion equation. In the context of numerical device simulations, the latter is associated with the well-known Scharfetter-Gummel (SG) discretization. Chapter 7. Conclusion 128 Therefore, our SD formulations are reasonable extensions of the SG formulation for gen-eral convection-diffusion systems in multi-space dimensions. 7.2 Future Research Direct ions Although the numerical method that we use in this thesis is robust in simulation for a wide parameter range, it is not entirely trouble-free as is reported in Chapters 5 and 6. From our numerical experience, as well as others' in the literature [22, 13], we found that the solution for the electron current, the important part in the simulation, is not completely satisfactory near sharp layers in comparison with other solution components. Numerical tests in the thesis indicate that the current error is sensitive to perturbations - including perturbations in the strength of the artificial diffusion introduced by the scheme in order to achieve numerical stability and to suppress oscillation at sharp layers. However, even though it is sensitive to perturbations, the approximate electron current obtained does not exhibit non-physical behavior, i.e., there was no recirculation of the electron flow (cf. [13]). It is thus of great interest to design a new approach to obtain a more reliable approximation for the electron current in sensitive regions. As far as the device simulations are concerned, our further research activities will involve introducing quantum corrections in the HD model. Some work has been done recently by Gardner et al (see [27, 28, 14]). It seems also relevant to tunneling devices like the hetero junction bipolar transistor (HBT). The application of the new method in Chapter 4 to the DD model for devices is of interest too. Bibliography [1] U. Ascher, On Numerical Differential Algebraic Problems with Application to Semi-conductor Device Simulation. SIAM Numer. Anal., Vol. 50, N o . 3, pp.517-538, (1989). [2] U. Ascher, P.A. Markowich, P. Pietra and C. Schmeiser, A Phase Plane Analysis of Transonic Solutions for the Hydrodynamic Semiconductor Model. Math. Models and Methods Appl. Sci., Vol. 1, N o . 3, pp.347-376, (1991). [3] U. Ascher, P.A. Markowich, C. Schmeiser, H. Steinriick and R. Weiss, Conditioning of the Steady State Semiconductor Device Problem. SIAM J. Appl. Math., Vol. 49, N o . 1, pp.165-185, (1989). [4] R.E. Bank, D.J. Rose and W. Fichtner, Numerical Methods for Semiconductor Device Simulation. IEEE Trans. El. Dev., Vol. 30, N o . 9, pp.1031-1041, (1983). [5] R.E. Bank, W. Fichtner, D.J. Rose and R.K. Smith, Computational Aspects of Semiconductor Device Simulation. Numer. Anal. Manuscript, 85-3, AT & T Bel Lab., (1985). [6] K. Blotekjaer, Transport Equations for Electron in Two-valley Semiconductors. IEEE Trans. El. Dev., Vol. ED-17. pp.38-47, (1970). [7] A. Benvenuti, W.M. Coughran, Jr., M.R. Pinto and N.L. Schryer, Hierarchical PDE Simulation of Nonequilibrium Transport Effects in Semiconductor Devices. NUPAD IV, pp.155-160, (1992). [8] A. Benvenuti, M.R. Pinto, W.M. Coughran, Jr., N.L. Schryer, C.U. Naldi and G. Ghione, Evaluation of the Influence of Convective Energy in HBTs Using a Fully Hydrodynamic Model. International Electron Device Meeting, pp.499-502, (1991). [9] G.J. Le Beau, S.E. Ray, S.K. Aliabali and T.E. Tezduyar, SUPG Finite Ele-ment Computation of Compressible Flows with the Entropy and Conservation Vari-ables Formulations. Comput. Methods Appl. Mech. Engrg., Vol. 104, pp.397-422, (1993). [10] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods. Springer-Verlag, New York, Berlin, Heidelberg, London, Paris, Tokyo, Hong Kong and Barcelona, (1991). 129 Bibliography- ISO [11] F. Brezzi, L. Marini and P. Pietra, Two-Dimensional Exponential Fitting and Ap-plications to Semiconductor Device Equations. SIAM J. Numer. Anal., Vol. 6, N o . 6, pp.1342-1355, (1989). [12] E. Carmo and A. Galeao, Feedback Petrov-galerkin Methods for Convection-Dominated Problems. Comput. Methods Appl. Mech. Engrg., Vol. 88, pp.1-16, (1991). [13] Z. Chen, B. Cockburn, J.W. Jerome and C.W. Shu, Mixed-RKDG Finite Element Methods for the 2-D Hydrodynamic Model for Semiconductor Device Simulation. J. Comput. Phys., to appear, (1995). [14] Z. Chen, B. Cockburn, C.L. Gardner and J.W. Jerome, Quantum Hydrodynamic Simulation of Hysteresis in the Resonant Tunneling Diode., IMA Preprint Series #1227, University of Minnesota, (1994). [15] R. Codina, A Discontinuity-Capturing Crosswind Dissipation for the Finite El-ement Solution of the Convection-Diffusion Equation. Comput. Methods Appl. Mech. Engrg., Vol. 110, pp.325-342, (1993). [16] F. Chalot, T.J. Hughes and F. Shakib, Symmetrization of Conservation Laws with Entropy for High-temperature Hypersonic Computations. Comput. Sys. En-grg., Vol. 1, N o . 2-4, pp.495-521, (1990). [17] R.K. Cook and J. Frey, Two-Dimensional Numerical Simulation of Energy Trans-port Effects in Si and GaAs MESFET's. IEEE Trans. El. Dev., Vol. ED-29 , N o . 6, pp.970-977, (1982). [18] P.G. Ciarlet, The Finite Element Method for Elliptic Problems. North Holland, Amsterdam, (1978). [19] P. Degond and P.A. Markowich, On a One-Dimensional Steady State Hydrodynamic Model for Semiconductors. Appl. Math. Lett., Vol. 3 , N o . 3 , pp.25-29, (1990). [20] E.P. Doolan, J.J.H. Miller and W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers. Boole Press, Dublin, (1980). [21] E. Fatemi, C.L. Gardner, J.W. Jerome, S. Osher and D.L. Rose, Simulation of a Steady-State Electron Shock Wave in a Submicron Semiconductor Device Using High-Order Upwind Methods. Computational Electronics: Semiconductor Transport and Device Simulation, pp.27-32, Boston, Kluwer Academic Publishers, (1991). Bibliography 131 [22] E. Fatemi, J.W. Jerome and S Osher, Solution of the Hydrodynamic Device Model Using High-Order Non-Oscillatory Shock Capturing Algorithms. IEEE Trans. Comput.-Aided Design, Vol. 10, N o . 2, pp.232-244, (1991). [23] L.P. Franca and T.J. Hughes, Two Classes of Mixed Finite Element Methods. Corn-put. Methods Appl. Mech. Engrg., Vol. 69, pp.89-129, (1988). [24] I.M. Gamba, Stationary Transonic Solutions of a One-Dimensional Hydrodynamic Model for Semiconductors. Technical Report N o . 143, Department of Mathemat-ics, Purdue University, (1990). [25] I.M. Gamba, Boundary Layer Formation for Viscosity Approximations in Tran-sonic Flow. Technical Report N o . 149, Department of Mathematics, Purdue Uni-versity, (1991). [26] C.L. Gardner, Numerical Simulation of a Steady State Electron Shock Wave in a Submicrometer Semiconductor Device. IEEE Trans. El. Dev., Vol. 38 , N o . 2, pp.392-398, (1991). [27] C.L. Gardner, The Quantum Hydrodynamic Model for the Semiconductor Devices. SIAM J. Appl. Math., Vol. 54, pp.409-427, (1994). [28] C.L. Gardner, Resonant Tunneling in the Quantum Hydrodynamic Model. VLSI Design, to appear, (1995). [29] C.L. Gardner, J.W. Jerome and D.J. Rose, Numerical Methods for the Hydrody-namic Device Model: Subsonic Flow. IEEE Trans. Comput.-Aided Design, Vol. 8, pp.501-507, (1989). [30] A. Gnudi, F. Odeh and M. Rudan, An Efficient Discretization Scheme for the Energy Continuity Equation in Semiconductors. Proceeding of SISDP, pp.387-390, (1988). [31] B. Gustafsson and A. Sundstrom, Incompletely Parabolic Problems in Fluid Dy-namics. SIAM J. Appl. Math., Vol. 35 , No .2 , pp.343-357, (1978). [32] J-M. Grygiel and P.A. Tanguy, Finite Element Solution for Advection-Dominated Thermal Flows. Comput. Methods Appl. Mech. Engrg., Vol. 93 , pp.277-289, (1991). [33] A. Harten, On the Symmetric Form of Systems of Conservation Laws with Entropy. J. Comput. Phys., Vol. 49, pp.151-164, (1983). Bibliography 132 [34] T.J. Hughes and A. Brooks, Multidimensional Upwind Scheme with No Crosswind Diffusion. FEM for Conservation Dominated Flows, T.J. Hughes et al eds., AMD, Vol. 34, ASME, New York, pp.19-35, (1979). [35] T.J. Hughes, L.P. Franca and G.M. Hulbert, A New Finite Element Formulation for Computational Fluid Dynamics: VIII. The Galerkin/Least-squares Method for Advective-diffusive Equations. Comput. Methods Appl. Mech. Engrg., Vol. 73 , pp.173-189, (1989). [36] T.J. Hughes, L.P. Franca and M. Balestra, A New Finite Element Formulation for Computational Fluid Dynamics: V. Circumventing the Babuska-Brezzi Condition: A stable Petrov- Galerkin Formulation of the Stokes Problem According Equal-Order Interpolations. Comput. Methods Appl. Mech. Engrg., Vol. 59, pp.85-99, (1987). [37] T.J. Hughes and L.P. Franca, A New Finite Element Formulation for Computa-tional Fluid Dynamics: VII. The Stokes Problem with Various Well-posed Bound-ary Conditions: Symmetric Formulations that Converge for all Velocity Spaces. Comput. Methods Appl. Mech. Engrg., Vol. 65 , pp.85-96, (1987). [38] T.J. Hughes, L.P. Franca and M. Mallet, A New Finite Element Formulation for Computational Fluid Dynamics: I. Symmetric Forms of the Compressible Euler and Navier-Stokes Equations and the Second Law of Thermodynamics. Comput. Methods Appl. Mech. Engrg., Vol. 54, pp.223-234, (1986). [39] T.J. Hughes, L.P. Franca and M. Mallet, A New Finite Element Formulation for Computational Fluid Dynamics: VI. Convergence Analysis of the Generalized SUPG Formulation for Linear-Dependent Multi-Dimensional Advective-Diffusive Systems. Comput. Methods Appl. Mech. Engrg., Vol. 63 , pp.97-112, (1987). [40] T.J. Hughes and M. Mallet, A New Finite Element Formulation for Computa-tional Fluid Dynamics: III. The Generalized Streamline Operator for Multidimen-tional Advective-Diffusive Systems. Comput. Methods Appl. Mech. Engrg., Vol. 58, pp.305-328, (1986). [41] T.J. Hughes and M. Mallet, A New Finite Element Formulation for Computa-tional Fluid Dynamics: IV. A Discontinuity-Capturing Operator for Multidimen-sional Advective-Diffusive Systems. Comput. Methods Appl. Mech. Engrg., Vol. 58 , pp.329-336, (1986). [42] T.J. Hughes, M. Mallet and A. Mizukami, A New Finite Element Formulation for Computational Fluid Dynamics: II. Beyond SUPG. Comput. Methods. Appl. Mech. Engrg. Vol .54, pp.223-234, (1986). Bibliography 133 [43] T.J. Hughes, T.E. Tezduyar, Finite Element Methods for First-Order Hyperbolic system with Particular Emphasis on the Compressible Euler Equations. Comput. Methods Appl. Mech. Engrg., Vol.45, pp.217-284, (1984). [44] J. JafFre, C. Johnson and A. Szepessy, Convergence of the Discontinuous Galerkin Finite Element Method for Hyperbolic Conservation Laws. Technical Report, N o . 1 9 9 3 - 1 1 / I S S N 0347-2809, Department of Mathematics, Chalmers University of Technology, The University of Goteborg, (1993). [45] X. Jiang, A Streamline-upwinding/Petrov-Galerkin Method for the Hydrodynamic Semiconductor Device Model. To appear in Math. Models and Methods in Appl. Sci., (1995). [46] C. Johnson and U. Navert, An Analysis of Some Finite Element Methods for Advection-Diffusion Problems. Analytical and Numerical Approaches to Asymp-totic Problems in Analysis, 0 . Axelsson, L.S. Frank and A. Van der Sluis eds., North-Holland, Amsterdam, pp.99-116, (1981). [47] C. Johnson, U. Navert and J. Pitkaranta, Finite Element Methods for Linear Hyperbolic Problems. Comput.Methods Appl. Mech. Engrg., Vol.45, pp.285-312, (1984). [48] C. Johnson and A. Szepessy, On the Convergence of a Finite Element Method for a Nonlinear Hyperbolic Conservation Law. Math. Comput., Vol. 49, N o . 180, pp.427-444, (1987). [49] C. Johnson, A. Szepessy and P. Hansbo, On the Convergence of Shock Capturing Streamline Diffusion Finite Element Method for Hyperbolic Conservation Laws. Math. Comput., Vol. 54, N o . 189, pp.107-129, (1990). [50] C. Johnson, A.H. Schatz and L.B. Wahlbin, Crosswind Smear and Pointwise Errors in Streamline Diffusion Finite Element Methods. Math. Comp., Vol. 49 , N o . 179, pp.25-38, (1987). [51] J.W. Jerome and C.W. Shu, Energy Models for One-Carrier Transport in Semicon-ductor Devices. ICASE Report, N o . 91-78, Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, Hampton, Virginia, (1991). [52] M. Kfizek and P. Neittaanmaki, Finite Element Approximation for a Div-Rot Sys-tem with Mixed Boundary Conditions in Non-Smooth Plane Domains. Aplikace Matematiky, Vol. 29, pp.272-285, (1984). Bibliography 134 M. Kfizek and P. Neittaanmaki, Finite Element Approximation of Variational Prob-lems and Applications. (1990). J.P. Kreskovsky, M. Meyyappan and H.L. Grubin, The Moments of the Boltzmann Transport Equation as Applied to the Gallium Arsenide Permeable Base Transistor. COMPEL, pp.99-105, (1987). Q. Lin, N. Goldsman and G-C. Tai, Highly Stable and Routinely Convergent 2-Dimensional Hydrodynamic Device Simulation. (1993), to appear in Solid State Electronics. Q. Lin, How to Sharpen the Error of Finite Element Methods. Preprint, Inst, of Systems Science, Academia Sinica, Beijing, (1994). T.J. Maloney and J. Frey, Transient and Steady-State Electron Transistor in GaAs and InP. J. Appl. Phys., Vol. 48, pp.781-787, (1977). P.A. Markowich, The Stationary Semiconductor Device Equations. Springer-Verlag, Wein, New York, (1986). P.A. Markowich and P. Pietra, A Non-Isentropic Euler-Poisson Model for a Colli-sionless Plasma. Preprint, (1991). P.A. Markowich, C.A. Ringhofer and C. Schmeiser, Semiconductor Equations. Springer-Verlag, Wien, New York, (1990). U. Navert, A finite Element Method for Convection-Diffusion Problems. Thesis, Chalmers Univ. of Tech. and Univ. of Gothenburg, (1982). F. Odeh, M. Rudan and J. White, Numerical Solution of the Hydrodynamic Model for a One-Dimensional Semiconductor Device. COMPEL, Vol. 6, pp.151-170, (1987). A.I. Pehlivanov, G.F. Carey and R.D. Lazarov, Least-Squares Mixed Finite Ele-ments for Second Order Elliptic Problems. SIAM J. Numer. Anal., Vol. 3 1 , N o . 5, pp.1368-1377, (1994). T.E. Peterson, A Note on the Convergence of the Discontinuous Galerkin Method for a Scalar Hyperbolic Equation. SIAM J. Numer. Anal., Vol. 28, N o . 1, pp.133-140, (1991). S.J. Polak, C D . Heijer, W.H.A. Schilders and P.A. Markowich, Semiconductor Device Modelling from the Numerical Point of View. Inter. J. Numer. Engrg., Vol. 24, pp.763-838, (1987). Bibliography 135 D.L. Pulfrey and N. G. Tarr, Introduction to Microelectronic Devices. Prentice Hall Series in Solid State Physical Electronics, N. Holonyak, Jr. Editor, Prentice Hall, Englewood Cliffs, New Jersey, (1989). K. Rahmat, J. White and D.A. Antoniadis, Computation of Drain and Substrate Currents in Ultra-Short-Channel NMOSFET's Using the Hydrodynamic Model. IEEE, IEDM, pp.115-119, (1991). J. Ruch, Electron Dynamics in Short Channel Field Effect Transistors. IEEE Trans. El. Dev., Vol. ED-19 , pp.652-654, (1972). M. Rudan and F. Odeh, Multi-dimensional Discretization Scheme for the Hydro-dynamic Model of Semicondctor Devices. COMPEL, Vol. 5, N o . 3, pp.149-183, (1986). W. V. Van Roosbroeck, Theory of Flow of Electrons and Holes in Germanium and Other Semiconductors. Bell Syst. Tech. J., Vol. 29, pp.560-607, (1950). D.L. Scharfetter and H.K. Gummel, Large-Signal Analysis of a Silicon Read Diode Oscillator. IEEE Trans. El. Dev., Vol. ED-16 , pp.64-77, (1969). S. Selberherr, Analysis and Simulation of Semiconductor Devices. Springer-Verlag, Wien, Nwe York, (1984). M. Sever, Symmetric Forms of Energy-Momentum Transport Models, preprint, (1991). M. Sever, Delaunay Partitioning in Three Dimensions and Semiconductor Models. COMPEL, Vol. 5, N o . 2, pp.75-93, (1986). F. Shakib and T.J. Hughes, A New Finite Element Formulation for Computational Fluid Dynamics: IX. Fourier Analysis of Space-Time Galerkin/Least-Squares Al-gorithms. Comput. Methods Appl. Mech. Engrg., Vol. 87, pp.35-58, (1991). M.S. Shur, Influence of Nonuniform field Distribution on Frequency Limits of GaAs Field Effect Transistors. Electron Lett., Vol. 12, pp.615-616, (1976). F. Shakib, T.J. Hughes and Z. Johan, A New Finite Element Formulation for Computational Fluid Dynamics: X. The Compressible Euler and Navier-Stokes Equations. Comput. Methods Appl. Mech. Engrg., Vol. 89, pp.141-219, (1991). M. Sharma and G.F. Carey, Semiconductor Device Modeling Using Flux Upwinding Finite Elements. COMPEL, Vol. 8, N o . 4, pp.219-224, (1989). S.M. Sze, The Physics of Semiconductor Devices. John & Sons, New York, (1981). Bibliography 136 [80] Thomann, E. and Odeh, F., On the Well-Posedness of the Two-Dimensional Hy-drodynamic Model for Semiconductor Devices. COMPEL, Vol. 9, N o . 1, pp.45-57, (1990). [81] C.H. Wu, Numerical Analysis of the Time-Dependent Two and Three Dimensional Semiconductor Device Equations. Ph.D. Thesis, Univ. of Dublin, (1990). [82] K. Yosida, Functional Analysis. Sixth Edition, Springer-Verlag, Berlin, Heidelberg, New York, (1980). [83] G. Zhou, An Adaptive Streamline Diffusion Finite Element Method for Hyperbolic Systems in Gas Dynamics. Thesis, Universitat Heideberg, (1992). Append ix A A Proper ty for the Hydrodynamic Mode l (2.1) For simplicity, only a one-space-dimensional problem is considered here. Multi-space dimensional problems can be treated similarly with a little extra work. We rewrite the conservation laws (2.1) in the following vector form Ut + AUX + C(U) = {KUX (A.l) where u* = (P, J, wy (A.2) and ' 0 0 0 ^ K = 2K 0 0 0 J2 w J \ 7* -J- 1 The matrix A can be diagonalized by a non-singular matrix Y, such that (A.3) Y lAY = A = diag(v - c, v, v + c), (A.4) where Y 1 i ! v — c \ H -cv f H + cv J V V + c „2 (A.5) 137 Appendix A. A Property for the Hydrodynamic Model (2.1) 138 and Y~1 = 1-S 1 u6 - 6 KUS-V \(\-vb) | , (A.6) 3c2 In the above, c2 = f-, # = 3f + f, 5 = |% and 6 = A (see [22]) Proposition A . l There exists Ap = diag(flx,l32,f3z) with fa > 0 (i = 1,2,3) such that ApY_1KY is symmetric positive semi-definite. Proof: Denote the last row of K Note that It is easy to see that J2 W J nu) = (-j- - , - - , i ) . r P P Y-^K = ^b ( i \ 2 v t y ^(C/) ^ F = T ( 1 , - - , 1 ) (A.7) (A.8) (A.9) kp = diag(4,3,4) (A.10) satisfies what we need. • A p p e n d i x B Stabil i ty of the SD-Trapezoidal Scheme We consider the SD-trapezoidal scheme (3.76) in the form TTh(n+1) _ Tjh(n) Jjh^1) _ Tjh(n) (A, 5 i , • ) + ( * Xt ,P"W) + S(i(£/',("+1) + f/''(")),*) = 0, (B.l) for any $ e V l . As we pointed out in section 3.5.2, existence of the second term above restricts us in obtaining a stability result for the general case. A restriction on the diffusion matrix is needed to obtain the following stability result. T h e o r e m B . l If K = KII then the one-dimensional SD-trapezoidal scheme (B.l) is unconditionally stable. Proof: Let Uh{n) = A^1'2 ZA^l,2Vh{n\ * = A^Z^A1/2®, then (B.l) is transformed into T//I("+1) _ yh(n) yh(n+1) _ yh(n) ( a= , •) + ( s i ' A A ^ } +\{(Hvh(n+l) + vh(%^) (K(vh{n+1) + vh(n))x^x) +(A(Vh{n+1) + Vh{n\, AATV>*)} = 0, (B.2) where, AT is given in (3.28). Taking tf = AT(V fc (n+1) - Vh{n))/At in (B.2), we have -^\\A\l\Vh{n+1) -Vh{n))\\2 139 Appendix B. Stability of the SD-Trapezoidal Scheme 140 2At 1 { ( A ( ^ ( n + 1 ) + Vh{n\ , A T ( ^ ( n + 1 ) - V*(B))) +(K(Vh{n+1) + Vh{n\, AT(Vh{n+1} - Vh(n}] +(||ATAyj l (n+1) | |2 - ||ATAV^(?l)||2)} = 0. (B.3) Taking $ = (Vh{n+1) + Vh{n))/2, we have 1 2At \yh(n+1)\\2_Uyh(.n)\\2\ I {yh(n+l) _ yh(n) ^ A T A ( ^ ( - + 1 ) + yh(n))x) 2At 1 + | |A i /2A ( yA(-+i) + yAW) r | | 2 } = 0 - ( B . 4 ) We define ||| • |||, an auxiliary norm, by lll^lll2 = II V||2 + ||ArAK||2 + K(ATVX,VX). (B.5) Then (B.3) + (B.4) yields ^ _ / | | | y / i ( n + 1 ) | | | 2 _ | | | T / f c ( n ) | | | 2 \ 2Ar'" '" '" '" ; l-\\A\l\Vh{n+1) -Vh{n))f At21 +^ | |Ay 2 A(^ ( n + 1 ) + V*(n))«||2 < 0. (B.6) Applying the inequality a2 + 62/4 > ab to the last three terms above, we obtain |||yfc(»+l)|||2 < |||yfcW|||2n ( B J ) Fourier analysis applied to (B.l) for the one-dimensional linear scalar convection-diffusion problem was carried out in [40]. It was shown that (B.l) can be considered as a coun-terpart of the one-pass implicit predictor/multi-corrector algorithm for the linear-linear SD method. Appendix C 1-D Test for a Set of Parameters DENSITY Legend 1.0 Volt 1.5 Volt £.)} VOIl \ I I ^— 1 0 0.1 02 OS 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ELECTRICAL FIELD Legend 1.0 Volt 1.5 Volt ^.u von \ x-;'i 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Figure C.l: HD model at 300K: using Fatemi et a/.'s parameters //0 = 0.145, vs = 0.1 and K = to2k. 2e 141
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical simulations of semiconductor devices by streamline-diffusion...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical simulations of semiconductor devices by streamline-diffusion methods Jiang, Xunlei 1995
pdf
Page Metadata
Item Metadata
Title | Numerical simulations of semiconductor devices by streamline-diffusion methods |
Creator |
Jiang, Xunlei |
Date Issued | 1995 |
Description | Theoretical and practical aspects of the design and implementation of the streamlinediffusion (SD) method for semiconductor device models are explored systematically. Emphasis is placed on the hydrodynamic (HD) model, which is computationally more challenging than the drift-diffusion (DD) model, but provides some important physical information missing in the DD model. We devise a non-symmetric SD method for device simulations. This numerical method is uniformly used for the HD model (including a proposed simplification (SHD)) and the DD model. An appropriate SD operator is derived for the general non-symmetric convection-diffusion system. Linear stability analysis shows that our proposed numerical method is stable if the system can be symmetrized. Stability arguments and numerical experiments also suggest that the combination of the method of lines and the semidiscrete SD method may not be appropriate for the transient problem, a fact which often has been ignored in the literature. An efficient method, consistent with the SD method used for conservation laws, is developed for the potential equation. The method produces a more accurate electric field than the conventional Galerkin method. Moreover, it solves for the potential and electric field in a decoupled manner. We apply our numerical method to the diode and MESFET devices. Shocks for the diode in one and two space dimensions and the electron depletion near the gate for the MESFET in two space dimensions are simulated. Model comparisons are implemented. We observe that the difference in solutions between the HD and DD models is significant. The solution discrepancy between the full HD and SHD models is almost negligible in MESFET simulation, as in many other engineering applications. However, an exceptional case is found in our experiments. |
Extent | 11534583 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-06-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080025 |
URI | http://hdl.handle.net/2429/8795 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1995-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1995-982942.pdf [ 11MB ]
- Metadata
- JSON: 831-1.0080025.json
- JSON-LD: 831-1.0080025-ld.json
- RDF/XML (Pretty): 831-1.0080025-rdf.xml
- RDF/JSON: 831-1.0080025-rdf.json
- Turtle: 831-1.0080025-turtle.txt
- N-Triples: 831-1.0080025-rdf-ntriples.txt
- Original Record: 831-1.0080025-source.json
- Full Text
- 831-1.0080025-fulltext.txt
- Citation
- 831-1.0080025.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080025/manifest