COMPUTER MEHODS NORTHHOLLAND
IN APPLIED
MECHANICS
AND ENGINEERING
75 (1989) 493514
NUMERICAL SIMULATION OF SEMICONDUCTOR DEVICES Franc0 BREZZI”,
Luisa Donatella
MARINI,
Istituto di Analki Numerica de1 C. N.R.,
Paola PIETRA
Pavia, Italy
The theoretical analysis of the mixed and hybrid exponential fitting methods, introduced in [I], is extended to the case of a nonzero source term. The methods are then studied from the numerical point of view and their upwinding features are analysed.
1. Introduction The equations commonly used for semiconductor device simulation are the Van Roosbroeck equations [241. They are considered to be a reasonably satisfactory model for materials like silicon or germanium and within a certain range of applications. For the use of other possible models better suited to describe the behaviour of other semiconductor materials such as gallium arsenide, we refer for example to [571. The Van Roosbroeck equations involve a number of physical quantities and parameters with enormously different orders of magnitude (from the 10l’ A of the electron charge to the 1017 cmm3 of a realistic doping profile). It is clear that in this situation a scaling is convenient (and recommended). We refer to [3] and [8] for two different types of possibe scalings and we present here the Van Roosbroeck scaled equations obtained following [3] in the special case of a stationary problem with constant coefficients: IAQ!l=np+L), div (V‘piz  nVa,!t)= R , div(Vp+pV$)=R, + boundary conditions
in f2CR!“, in f2 , inR,
(1.1)
.
The unknown JI, IE and p represent the (scaled) electric potential and the (scaled) concentration of negative and positive charges, respectively. D describes the doping profile and R is the socalled generationrecombination term which takes into account the possible ‘birth’ and ‘death’ of charges. The domain 0 is scaled to have diameter equal to 1. The first equation of (1.1) is the classical Poisson’s equation for the potential, and the second and third equations are the continuity equations for negative and positive current densities J”, Jp: J”=VnnVJi,
J”=Vp+pV$.
The choice of boundary conditions depends on the applications. As an example, we can assume the boundary to be split into two disjoint parts: afl = r, U r,. On r,, +, n and p are * Also at Dipartimento 00457825/&9/$3.50
0
di Meccanica Strutturaie
1989, Elsevier
dell’Universit5,
Science Publishers
Pavia, Italy.
B.V. (NorthHolland)
494
F. Bred
et al., Numerical simulation of semiconductor
devices
prescribed, while on r, vanishing outward normal components for J”, J’, and for the electric field E = VI& are assumed. The parameter 1 is the dimensionless scaling parameter, whose value lies typically between low3 and 10? For results and additional references on the limit behaviour of the problem for different scalings we refer to [3] and [9]. Existence results for the stationary problems hold under fairly general assumptions [lo, 3,4], while uniqueness is not true in general. It holds only in some special cases [3,4, 10,9]. When dealing with system (1.1) from a numerical point of view, one key point is the coupling of the equations. The unknowns are strongly coupled and the choice of the iterative method to decouple them is a difficult task. We recall here that the most used iterative algorithms are variants of the Newton method, plus some other tricks, such as continuation in some parameter. We refer to [11,12,3] for detailed descriptions of possible algorithms. A second problem to be faced is the choice of a proper discretization for the single equations. The Poisson equation is not a source of problems; every reasonable approximation scheme will work. Difficulties arise when discretizing the current equations, which are advection dominated (due to the large size of Vy5in a consistent portion of the domain). The most popular and effective method for onedimensional problems is the socalled ScharfetterGummel method [13], better known in other fields as the exponential fitting method [14,15]. The main feature of this scheme is to be currentpreserving. This property is particularly important in the applications, as the current is possibly the most relevant unknown of the problem. In two dimensions many methods have been proposed: Mock’s method [16], the Boxscheme (which is the most widely used in the codes) [17] and, more recently, MarkowichZlamal [18] and the modified method of characteristics (the latter being essentially introduced for transient problems) [19]. All except the latter, are attempts to extend the ScharfetterGummel method to 2D problems. In the present paper we concentrate our analysis on the socalled 2D exponential fitting methods recently introduced in [20, l]. These new schemes, essentially based on mixed or hybrid finite element methods, are both current preserving, and appear to be a very natural and satisfactory extension of the ScharfetterGummel method. In order to describe the methods, we shall focus our attention on one equation only, the current equation for positive charges. We shall denote by u the unknown and assume + to be a given function (piecewise linear in the computations). Find u E Hi(a), div(Vu+uV$)=f, u=g,
such that in Rc:iR2, on ToCall,
(l2)
The structure behind the methods is the following. First, we introduce variables from u to the Slotboom variable p u = p em* .
In terms of p, problem
the classical change of
(1.3) (1.2) can be written in the symmetric
form:
F. Brezzi et al., Numerical simulation of semiconductor
Find p E H’(R), div(
[email protected]) = f ,
such that in R ,
p=,y:=eJg,
on G ,
devices
495
(1.4)
ap
0, anThen, we discretize problem (1.4) with mixed or hybrid schemes. Finally, we go back to the variable u by means of (1.3) and we solve in U. If a proper decomposition is introduced, the discretization of (1.4) provides a symmetric positivedefinite Mmatrix and the Mmatrix property is preserved in the final problem for the variable U. An outline of the paper is as follows. In Sections 2 and 3 we introduce the mixed and hybrid approximations of problem (1.4)) respectively and we extend the analysis of [l] to the case f # 0. Such an extension is straightforward for the mixed method, while it needs some care when dealing with the hybrid approach. In Section 4 we quit the theoretical framework of [l] and we go into more numerical details. In particular, we analyse the qualitative behaviour of the schemes for large ]VI,!J]an d underline their upwinding features. We also point out some numerical difficulties that can arise. Section 5 is devoted to the presentation of numerical results obtained using both schemes on a number of examples.
2. Mixed methods Let us consider (see [21]) a regular sequence {T,}, of decompositions of 0 into triangles T (for the sake of simplicity we shall assume that J2 is a polygon). According to [22], we define the following set of polynomial vectors: RT(T) = {T = ( 71772) and, for every triangulation V, = {T E
171 = Q + Px,
a, P, Y E R> ,
&
Then, the mixed discretization
E P,,(T) VT E Th} .
In
ph
(2.2) (2.3)
of (1.4) is defined as follows:
Find J,, E V,, and ph E IV’, , dx +
(2.1)
Th, we set:
[L’(L’)]’ 1d’lv7EL2(R),7.n=OonT,,71.ERT(T)VTET,},
I+‘/,= {cp E L’(fl)]
e'Jh a 7
72 = Y + PY,
such that
div T dx =
It is clear that p,, will be an approximation approximation of the current
xr.ndT
VTEV,,
of the solution
,
p
(2.4)
of (1.4) while Jh will be an
496
F. Brezzi et al., Numerical simulation of semiconductor
devices
J=e*Vp=Vu+uV+.
(2.5)
In particular, the first equation of (2.4) is a discretized version of (2.5), and the second equation of (2.4) is a discretized version of div J = f. In the computations one usually assumes that $ is piecewise linear in each T, so that the integrals in (2.4) can be computed exactly. Uniqueness results for problem (2.4) follow easily from the general theory of [23]. We remark that the condition div T E L2(6!) in the definition (2.2) of V, implies that every r in V, has a continuous normal component when passing from one element to another. This means, in particular, that the current J,, is preserved. The algebraic treatment of the system (2.4) needs some care. The matrix associated with (2.4) has the form A
B
[ B*
0
1
P6)
and it is not positivedefinite. A way to avoid this inconvenience was introduced by Fraeijs de Veubeke [24] for a similar problem. The idea is to relax the continuity requirement in the space definition (2.2) and to enforce it back by using interelement Lagrange multipliers. We give now a more detailed description of the procedure. For this, we first set ch = {rE[L’(R)]‘]
&ERTVTE
Th}.
(2.7)
Then, calling E, the set of edges of Th, we define, for every function 5 E L2(&): (2.8) We consider now the mixedequilibrium
discretization
of problem
(1.4):
(2.9)
It is easy to see that the linear system (2.9) has a unique solution and that & = J,, ,
6, =
Ph
*
Moreover, h, is a good approximation of p at the interelements The matrix associated with (2.9) now has the structure
(2.10) (see [25] for detailed proofs).
(2.11)
F. Brezzi et al., Numerical simulation of semiconductor
devices
497
where A is blockdiagonal (each block corresponding to a single element). In other words, Jh can be eliminated, element by element, by static condensation. This leaves us with the matrix B*A‘B C*AlB
B*A‘C C*A‘C
(2.12)
’
where B* A‘B is blockdiagonal (actually, in the present case, diagonal) and nonsingular. Hence, ph can also be eliminated, at the element level, by static_condensation. This leaves us, finally, with a matrix (acting only on the interelement multiplier Ah) which can be proved to be a symmetric positivedefinite matrix, which is. an Mmatrix if the triangulation is of weakly acute type (each angle of every T E T,, is SIT/~). This general strategy for dealing with problems of the type (2.4) can be simplified in the present case. Let us see it in more detail. Denoting by & the mean value of f in T and using definition (2. l), the second equation in (2.9) implies that, for every T in T,, we have divj~l~~2P=f,:=~~~fdx.
(2.13)
Hence, & can be split in two parts J,=J”+Jf,
(2.14)
Joi, = (a, r>,
(2.15)
where
(x being the coordinate vector x = (xi, x2)). Taking now in the first equation of (2.9) T = (1,0) and 7 = (0,l) in T, (T = 0 elsewhere), we obtain, with the notation of Fig. 1 e’ dx = JOIT
i
le,liknCk)  IT e’Jf dx .
(2.16)
k=l
We make for the basis functions p in _4h,othe obvious choice F = 1 on one edge and p = 0 on the others. Then, inserting (2.16) in the third equation of (2.9), one can easily compute the contributions of the triangle T to the stiffness matrix and to the righthand side: v(i)
. .(i)
a; =
7
IT
@” =

(2.17)
e’ dx
e* dx
where we set vCk)= lek(nCk) ,
k = 1,3.
(2.18)
As already pointed out, the integrals appearing in (2.17) can be computed exactly, since 4 has been assumed linear in each triangle. It is now easy to check that the final matrix, obtained by
F. Brezzi et al., Numerical simulation of semiconductor
498
devices
Fig. 1.
assembling the elementary cont~butions (2.17), is a symmetric positivedefinite matrix and it is an Mmatrix if the triangulation is of weakly acute type. Note that the matrix (2.17) corresponds to a linear nonconforming finite element approximation of problem (1.4) where the harmonic average of the coefficient e’ is taken instead of the usual average. We are now ready to introduce an approximation of our original unknown u. For this we define, for every 5 E L2(E,), 5’ to be the L2projection of 4 onto Ah,r, that is: (2.19) Using now the change of variable & = te~)~~~, we can rewirte (2.9) in the form: FindJ,EFh,
P~EI+‘~,u~EA~,~,
such that
(2.20)
where we went back to the original names for the variables Jh and ph in virtue of (2.10). If one eliminates Jh and ph in (2.20) by static condensation, the resulting matrix, acting on the unknown u,, only, is an Mmatrix, if the decomposition is of weakly acute type. Actually, it can be obtained from the matrix (2.17), after assembling, by multiplying each column by the value of (e’)’ on the corresponding edge. ~~~A~K
2.1. Different choices of mixed finite element methods can be made. We refer to [l] for the use of quadrilateral decompositions of ft in this context and to 1261 for general references on mixed methods.
We conclude this section with an error estimate result on the current. The proof is a slight modification of that given in [I, Theorem 3.11.
F. Brezzi et al., Numerical simulation of semiconductor
devices
499
THEOREM 2.1. Let u be the solution of (1.2), let J : =
[email protected](e’u) =Vu + uV$ and (J,, p,,) be the solution of (2.4). Then
IIJ J,k Q IIJ 
4JII.
(2.21)
7
where I(Jlli = Jn e*J  J d x, and II,,J E V, is defined locally as
I
e
(l&JJ)*nds=O
VeEE,.
(2.22)
3. Hybrid methods We consider now another class of nonstandard methods for dealing with problem (1.4), namely, the class of hybrid methods. Let again {T,,},, be a regular sequence of decompositions of 0 into triangles T [21]. For each T,, in the sequence and for every function 5 E Co(&) we introduce the sets V~={~E[L’(~)]*~~~T=(CY,Y)
VTET,,
[email protected]&
(3.1)
W h,S={~ECo(~)I~ITEP1(T)
VTET,,cp=tatnodesonT,}.
(3.2)
Then the hybrid discretization
of (1.4) can be written as [26]:
Find (Jh, ph) E V, x W,,x, I
I
OJ,,V~dx=
such that
nr.Vphdx=O cpJf*nds
VrEVh,
(3.3)
VP E W/z,0*
In (3.3) Jf should be a solution of the equation div Jf = fin T for all T. Without substantial loss of accuracy, one can choose instead Jf satisfying div Jf = fT, where fr is the mean value of f in T, so that Jf I T can be chosen, for instance, as in (2.15). The integrals appearing in (3.3) can be computed exactly, since @ is piecewise linear. However, for reasons that will be made clear in the next section, it is convenient to use a suitable numerical integration. For the moment, let us denote by I, the quadrature formula
(3.4) Note that I, is still to be chosen. However, we can already assume (3.4) to be exact for a constant cp. In order to simplify the notation, we can introduce the piecewise constant function _ $ defined in each T by IT
ei = Z,(e’)
,
and then use 4 instead of I,!I in the stiffness matrix of (3.3).
(3.5)
F. Bretzi et al., Numerical simulation of semiconductor
500
devices
It is easy to check that problem (3.3) has a unique solution [23]. It is be an approximation of the solution p of (1.4), while J,, f Jf will be an current J. In particular, the first equation of (3.3) is the discrete version equation of (3.3) ensures a weak conservation property of the current. parts, it implies Jh
+Jf)*
ncp
ds = 0
Vq E W,,,
also clear that ph will approximation of the of (2.5). The second After integration by
(3.6)
.
Then, with the obvious choice for the basisfunctions in W,,, : cp = 1 at a node (=vertex) P and cp = 0 at the other nodes, one obtains from (3.6) that the jumps of the flux of Jh + J* across the edges having P in common sum up to zero. For hybrid schemes, the solution of the algebraic system associated with (3.3) is easier than for mixed methods. As before, the matrix form of system (3.3) is the following: (3.7) where here A is a blockdiagonal matrix, since no continuity requirement is made in V,. One can then eliminate J,, by static condensation. The resulting matrix H = B*A‘B, acting on the nodal values of p,,, is symmetric and positivedefinite. Also, it will be an Mmatrix if the triangulation is of weakly acute type. For the convenience of the reader, we write explicitly the contributions of an element T to the stiffness matrix and to the righthand side
a*;=
v(i)
. ,(i)
,
(3.8)
J
Te*Jf*V(i)dX
br=

2/T
+ e* 
dx
J
T Jf.~“‘dx 2lTI
;
tTifT>
(3.9)
where
[email protected]) is defined as in (2.18). One can then go back to the unknown u,, via the discrete analogue of (1.3); this corresponds to multiply the jth column of the matrix H by eJs’and does not change the Mproperty of the matrix. REMARK
3.1. Note that the matrix (3.8) corresponds to a linear conforming finite element approximation of problem (1.4), where a sort of harmonic average of the coefficient e’ is taken instead of the usual average. More precisely, since (3.5) implies: e&l, =
a linear conforming follows:
$$j , T approximation
(3.10) of problem (1.4) with ‘harmonic’ average can be written as
F. Brezzi et al., Numerical simulation of semiconductor
501
devices
(3.11)
It is clear that the matrix associated with (3.11) coincides with the matrix H given by (3.8), as we said before, while the contribution of each triangle to the righthand side of (3.11) is given bY (3.12) which is much simpler than (3.9). On the other hand, (3.11) will produce an approximation the current
of
(3.13) which will no longer have continuity properties. However, since the difference between the two problems only consists of a different handling of the righthand side, the two solutions will coincide if we can choose Jf in (3.9) such that
I
e+Jf.
y(i)
IT
,ci)
Jf.
dx
T +
2 Te”dx
dx
=0
WI
VT E T,,
VY(')
.
(3.14)
I
This is equivalent
 I
to choosing Jf such that
T e’Jf dx
2 T e’dx I
I
+
TJfdx
2lTI
=0
VTET,.
(3.15)
Note that in the hybrid formulation (3.3) Jf is any solution of the equation div Jf = fr in T for all T. Hence, we can look for Jf = (.I{, Ji) of the form
(3.16)
wherex”=xx,,r,= coordinates of the barycenter of T. In (3.16) 0, and f3,are to be chosen in order to verify (3.15). By inserting (3.16) in (3.15) we obtain a system of two equations in the unknowns til and 19,which is always solvable. Also, we point out that
e: + ei = (fT)2/4.
(3.17)
F. Brezzi et al., Numerical simulation of semiconductor
502
devices
With this choice for Jf in (3.3) one clearly obtains Ph =
6,
7
(3.18)
I
J,, =jh 
e'Jf
dx
T
IT
. e’ dx 
In view of these considerations, we therefore suggest the following approach: to use in the computations the simplified form (3.12) for the righthand side, instead of the original one (3.9). Then, if one is interested in having an approximation of J with reasonable continuity properties (of the type (3.6)), use (3.18) as a postprocessing. We explicitly point out that the computational cost of (3.18) matches closely the difference in computational cost between (3.9) and (3.12). The advantage here is that, in many applications, one would use (3.18) only one might choose to in a few triangles. Besides, when using an iterative procedure, postprocess only the last iteration and not the intermediate ones. REMARK
3.2.
For hybrid methods
over quadrilateral
decompositions
of 0 in this context,
we refer to [l]. We conclude this section with some error estimates on the current density J. First, let us note that from (3.16), (3.17) and (3.18) one easily obtains ]]Jh j&(T)
S chfT
V’T E Th .
(3.19)
Hence, we shall prove estimates for J  & only. For that, we assume f to be piecewise constant and, in order to simplify the notation, we introduce the bilinear form: a(J, 7) =
IR
e’J.Tdx,
(3.20)
and the associated norm (3.21) With this notation,
(3.11) can be written as
Find (jh, ph) E V, X wh,x,
such that
dh,
VTEV~,
T, 
jhVpdx=
R fqdx
I
Note that the first equation theorem.
(3.22)
kpSl&.
of (3.22) is nothing
but (3.13). Then, we have the following
F. Brezzi et al., Numerical simulation of semiconductor
503
devices
THEOREM 3.1. Let u be the solution of (J.2) and let J :=
[email protected](e”u) =Vu + uV$. Let (J”,, p,,) be the solution of (3.22) and let P,J and J be the approximations of J defined by and
P,JEV,
and
J” E V,
(J
phJ,~)=o
STEVE
(3.23)
7
J”= e~V(~~~)~(=~~V~~) ,
(3.24)
where z’ E Wh,X denotes the piecewise linear interpolant of z. Then (3.25)
PROOF.
We have
llJjJ~ = a(JJ,, J&h) = a(J&,
J
P,J) + a(J&,
P,J&) (3.26)
=1+11
[email protected]&,
&JJ,)
=a(JJ”,P,Jj*)+a(jJ;,,P,JJo) (3.27)
= a(J  .k PhJ  &) . In fact a(j
&, p,J  j,) = R e4(epaVp’  e%‘p,)(P,J I =
f
a
(‘pb’v~~)(~J
 jh) dx
 .&) dx
= I a (vi+Vp,)(Jj,z)dx = !a fb’p,)dxjfl(Vp’Vp/,)J,dx=O.
(3.28)
In (3.28) we used the de~nition (3.23) of P,J, the fact that div J = f and pf  p,, E W~,~. From inequality (3.26) and (3.27) we have then, via the CauchySchwartz
llJ~~ll~~llJ~~ll.llJ p,Jll, + ~lJ~ll~ll~~j~ll~ s IIJ ~ll,(llJ  ~~11,+ IIJ  jll,> + IIJ  3”ll,llJfrom which (3.25) follows easily.
Cl
p,~ll,
(3.29)
F. Brezzi et al., Numerical simulation of semiconductor
504
devices
4. Numerical aspects As already pointed out, in the applications to semiconductor device simulation, VI) can be quite large in a part of the domain. This is clearly a potential source of numerical difficulties. More generally, the presence of a big VJ/ can change the qualitative behaviour of the scheme. In the present section we shall deal with some numerical aspects related to this problem. To this aim, it will be more convenient to set
and assume that Vet is smooth everywhere
and 1 is smafl. Accordingly,
(1.2) becomes
=f_
(4.2)
In what follows, we shall first justify, in (a), the need of a numerical integration already pointed out for the hybrid scheme of Section 3. Then, in (b) and (c) we shall assume V&, to be constant and show the qualitative behaviour of the mixed and hybrid scheme, respectively, for 1 << h. Finally, in (d) and (e) we shall briefly discuss the numerical difficulties that can arise in the two schemes when V&, changes abruptly from one triangle to another. (a) Let us first consider the problem of the choice of the integration formula (3.4). We recall that the cont~butions of each triangle T to the stiffness matrix, acting on the final unknown uh, and to the righthand side are respectively a$ =
y(i) * Y(j) e3 (4.3)
4Me’>
and (4.4) A quick glance at (4.2) suggests that, as I+ 0, the behaviour of a; should be of the order of magnitude of I‘. Assume now that I& is linear in 7’ and assume that I CC [V&,jh,. If you neglect the very special cases in which $ reaches its maximum +M on one edge, you are left with the generic case in which @ reaches its maximum only at one vertex. A simple computation shows that, in this case, (4.5) As a consequence, +
the use of exact integration
2~~‘)* ~“‘141~1Tf , 0,
which is not very satisfactory.
in (4.3) would give
if J/i = J/M , otherwise ,
Hence we need an integration
(4.6) formula such that
F. Brezzi et al., Numerical simulation of semiconductor
devices
505
(4.7) For that, we made the following choice: let V, and V, be (any) two vertices where $ reaches its minimum and (respectively) its maximum (you can have a nonunique choice for them, but this does not matter). Let, accordingly, e be the edge of T joining VM to V,. Then, we set
z,(
[email protected]) =#
1e’ dx ,
(4.8)
It is easy to check that choice (4.8) fulfills requirement (4.7). On the other hand, in the triangles T where \VI,!I/ h is not too big, formula (4.8) gives a reasonably good approximation of the exact integral over T, so that we suggest to use it for every triangle T in T,, , as we did in our computations. In order to convince the reader that the use of exact integration in (4.3) leads (for f # 0) to the wrong result, we refer to the numerical experiments of the next section. (b) We now want to analyse the qualitative behaviour of the mixed scheme (2.20) (in the unknown u,,) as I* 0. For this, we consider as a model case the situation of Fig. 2, where P is the midpoint of the common edge of the two equilateral triangles T, and T,. It is clear that the row of the matrix associated with the node P will have only five coefficients that (a priori) are different from zero. We denote them by ap,p and ap,p, (k = 1,2,3,4) with the notation of Fig. 3. For every 8 in [0,2n[ we consider a constant value of VI& = (cos 8, sin 0) on T, U T2, and we compute the coefficients of the row associated with P for 1 << h. In Fig. 4 we can see the behaviour (as functions of 0) of the coefficients ap,p and ap.p, for h = 1 and? = 10m4.The other three coefficients can be obtained from ap,p, by the formulae
+4T). ( 5>=af,f,(e + 4 = af,Je
a,,,,@)= af,f2 0 +
(4.9)
Varying h, the corresponding coefficients vary linearly with h, as far as 1 <
Fig. 2.
506
F. Brezzi et al., Numerical simulation of semiconductor
devices
Fig. 3.
Now, if VI,!+, points into a region R;, only the nodes Pj_l, Pi, Pj+l will be active (i.e. the corresponding coefficient will be different from 0). If instead, VI),,points into a dashed region which is between R, and Ri+l, then only the nodes Pi and Pi+l will be active. Hence the scheme produces a sort of upwinding effect, It can be seen that the matrix produced by the scheme corresponds roughly to the use of an artificial viscosity of the order of h/4. (c) An analogous upwind behaviour is exhibited by the hybrid scheme of Section 3. To see that more clearly, let us consider the simplified situation of Fig. 5, where six equilateral triangles join at a node P. The row of the matrix associated with P will have seven coefficients (a priori) different from zero: u~,~, u,,,~ (k = 1, . . . , 6). As before for every 8 E [0,21r[ we consider a constant value of VI,!+, = (cos 8, sin 0) and we compute the coefficients of the row corresponding to P for I CCh. Figure 7 shows the behaviour of the coefficients up +,and u~,~, as functions of 0 for h = 1 and I = 10e4. The other coefficients are obtainable from u~,~, by the formula aP,P). 2
1
w\,
k
188 I
368 Tbm\ <\\
2
aP,P I 1 Fig. 4.
',
0
F. Brezti et al., Numerical simulation of semiconductor
devices
507
Fig. 5.
%,P,W= UP,P,(e+(kI);), k=2
,...,6.
We can see that the scheme has always at least three coefficients equal to zero, corresponding to the three downwind nodes. Figure 6 helps in understanding how the upwind nodes change with Vt,bo.If VI,!+,points into a region Ri,then the upwind nodes are P,_i, Pi,Pi+l. If instead V& lies along RiflRi+l,then the upwind nodes are Piand Pi+l. In this respect, the scheme exhibits the same upwind nodes of other upwind schemes [27], although with a different conception and different numbers. As such, it can be considered to have an artificial viscosity roughly of the order of h/3. (d) We assumed so far that V$ does not change from one triangle to another. As we have seen, with this assumption both schemes produce, even for very small I, a quite reasonable matrix that corresponds roughly to the use of some kind of artificial viscosity. We shall see now that the situation is less comfortable if we admit the possibility of a VI) which changes abruptly from one triangle to the neighbouring one. In order to see what can happen it is
Fig. 6.
F. Brezzi et al., Numerical simulation of semiconductor
508
0
I 0
188
360
devices
ec
1
Fig. 7.
convenient to have one more look at the contributions of each triangle to the final matrix in the unknown u,. Let us consider the mixed scheme first. Remember that each of the final degrees of freedom is associated with one edge of Th (or, if you prefer, with the midpoint of the edge). The contributions of T to the coefficients aij of the final matrix are (4.11)
where we recall that here (4.12) A simple computation shows that, if I,+~’is the maximum value of + = ++,/I over ei (and if + is not constant on e,), then (e’);,
= I eSM’ .
It is now clear from (4.11)(4.13)
(4.13) and (4.5) that (4.14)
F. Brezzi et al., Numerical simulation of semiconductor devices
509
where eM denotes the maximum value of + in T as in (4.5). Hence, the scheme annihilates the contributions a; corresponding to the (nodes on) the edges ei where @ does not reach its maximum. In particular, when Vt,b changes from one triangle to another, the annihilation effect can lead to a column (corresponding to an edge ei of Th) of very small coefficients. This happens if +“j << $? k for ail edges ek (k # j) sharing a same triangle with ei. On the other hand, a row of very smalf coefficients can never occur. Actually, it is easy to see that equality eMj = eM occurs, in each triangle, for at least two different js. This implies, using (4.14), that a,: will be nonnegligible for at least one value of j. (Clearly, Y(~)+ v(j) cannot vanish for two different values of j). (e) Analogous considerations hold for the hybrid scheme of Section 3, if Vt,l~changes from one triangle to another. Looking at the cont~butions (4.3) of a triangle T to the stiffness matrix in the unknown uh, and using the integration formula (4.8), it is easy to see that
a;=
Y(~)* v(j)/41 1TI ,
if J/i= t,/.tM ,
0,
otherwise .
(4.15)
It is clear from (4.15) that the scheme annihilates the contributions a; corresponding to the nodes j where (tri<< e”. This can lead to a final matrix having either a row or a column of very small coefficients. In particular, assume that a node P is such that 4(P) << y$; t,b(x)for every triangle T having a vertex in P (i.e. $ has a very deep minimum in P, compared with the local mesh size). Then, the scheme will produce a column of very small coefficients. In order to see what can happen to a row, we first remark that, in general, for a given triangle T equality I,?J~ = eM will hold for one value of j only. If the decomposition is of weakly acute type, the corresponding contribution af can be zero if v(~) * v(j) = 0. Suppose now that you have a node Pi in the following unlucky situation: for every triangle T having Pi as a vertex jVt,!~lis big and the maximum of $I over T is reached at a point Pi such that ,@I * ycif = 0. Then, the ith row will be made of very small coefficients.
A
Gl Fig. 8.
F. Brezzi et al., Numerical simulation of semiconductor
510
devices
Fig. 9.
5. Numerical results We present here some numerical results obtained using the methods described in the previous sections over uniform triangular decompositions of the unit square (see Fig. 8). The function +(x) is assumed to be piecewise linear and given by the expression 4(x) = I&(X)/I, with I = 10m4 in all the examples. In the first example we show the resolution of a boundary layer (see Figs. 9,10). The function +,,(x) has the constant gradient VI&= (2, l), and f = 0. The boundary conditions for u are: u = 1 on the edges AB and AD, u = 0 elsewhere. In Fig. 9 the solution has been computed with the hybrid method (3.3) on a 20 X 20 mesh, and in Fig. 10 with the mixed method (2.20) on a 10 x 10 mesh. The boundary layer is catched in one element and no smoothing effect is introduced. Figure 11 shows how an internal layer is transported: I& is chosen as in the previous case, and f = 0. The boundary conditions for u are: u = 0 on AB, AF, and BF,, u = 1 elsewhere. The solution is computed with the hybrid method on a 40 x 40 mesh. (For the mixed method the situation is the same). In this case, a smoothing effect in the crosswind direction appears, as pointed out in Section 4. The example given in Fig. 12 refers
(1.0)
(O,l)
Fig. 10.
F. Brezzi et al., Numerical simulation of semiconductor
devices
511
Fig. 11.
Fig. 12.
to a favourable situation where the mesh is well oriented with respect to the transport. The function +,, has the constant gradient V& = ( 1, 0), f = 0 and the boundary conditions for u are: u = 0 on Al?, AE, and BE,, u = 1 elsewhere. Figure 12 shows the mixed solution on a 10 x 10 mesh. In this case the internal layer is as sharp as it can be. An analogous behaviour is exhibited by the hybrid scheme, if exact integration is used in (3.4). On the contrary, the use of a numerical integration as in (4.8) produces, even in this favourable case, a smoothing effect as that of Fig. 11. Figures 13 and 14 show the results obtained using the hybrid scheme when the source termf is different from zero. The function I,!Qhas the constant gradient VI),, = (2,0) and f = 2. The
Fig. 13.
F. Brezzi et al., Numerical simulation of semiconductor
512
devices
Fig. 14.
(0.0)
.
Fig. 15.
boundary conditions for u are: u = 1 on AD and BC, and homogeneous Neumann boundary conditions elsewhere. In the first case, the quadrature rule (4.8) is used, while in the second one the exact integrals are computed. As already pointed out in Section 4, the second choice leads to a discrete solution which is not an approximation of the solution of the continuous problem. In the last example, I& is chosen with the expression if 0 < p G 0.8 , if0.8sps0.9, if 0.9 Q p ,
Fig. 16.
F. Brezzi et al., Numerical simulation of semiconductor
devices
513
Fig. 17.
Fig. 18.
conditions for u are: u = 1 on AF, and on Neumann boundary conditions are assumed elsewhere. In Figs. 1618 we show u,, as computed, respectively, with hybrid methods (3.3) on a 20 x 20 mesh, mixed methods (2.20) on a 10 x 10 and 20 X 20 mesh. and p = (XT+ x~)“~ (see Fig. 15). The boundary
AG,, u = 0 on F2C and G,C. Homogeneous
REMARK 5.1.For graphical purposes, the solutions computed with the mixed scheme on the N x N mesh have been interpolated at the vertices and presented on a 2N x 2N grid. References exponential fitting and applications to driftdiffusion PI F. Brezzi, L.D. Marini and P. Pietra, Twodimensional models, SIAM J. Numer. Anal. 26 (1989). PI WV. Van Roosbroeck, Theory of flow of electrons and holes in germanium and other semiconductors, Bell Syst. Tech. J. 29 (1950) 560607. 131P.A. Markowich, The Stationary Semiconductor Device Equations (Spinger, Berlin, 1986). [41 M.S. Mock, Analysis of Mathematical Models of Semiconductor Devices (Boole Press, Dublin, 1983). PI H.U. Baranger, Ballistic electrons in a submicron semiconducting structure: A Boltzmann equation approach, Ph.D. Thesis, Cornell Univeristy, 1986. of electron transport properties in if51W. Fawcett, A.D. Boardman and S. Swain, Monte Carlo determination GaAs, J. Phys. Chem. Solids 31 (1970). [71 B. Niclot, P. Degond and F. Poupaud, Deterministic particle simulation of the Boltzmann transport equations of semiconductors, J. Comp. Phys. 78 (1988) 313349.
F. Brezzi et al., Numerical simulation of semiconductor
514 [8] F. Brezzi,
[9] [lo] [ll] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
[23] [24] [25] [26] [27]
devices
A.C. Capelo and L.D. Marini, Singular perturbation problems in semiconductor devices, in: Numerical Analysis, Lecture Notes in Mathematics, Vol. 1230 (Springer, Berlin, 1986) 191198. F. Brezzi, A.C. Capelo and L. Gastaldi, A singular perturbation analysis for reversebiased diodes, to appear in SIAM J. Math. Anal. 19 (5) (1989). J.W. Jerome, Consistency of semiconductor modeling: an existence/stability analysis for the stationary Van Roosbroeck System, SIAM J. Appl. Math. 45 (1985) 565590. H.K. Gummel, A selfconsistent iterative scheme for one dimensional steady state transistor calculations, IEEE Trans. Electron. Devices 11 (1964) 455465. S.J. Polak, C. Den Heijer, W.H.A. Schilders and P.A. Markowich, Semiconductor device modelling from the numerical point of view, Internat. J. Numer. Methods Engrg. 24 (1987) 763838. D. Scharfetter and H. Gummel, Largesignal analysis of a silicon Read diode oscillator, IEEE Trans. Electron. Devices 16 (1969) 6477. K.W. Morton, Generalised Galerkin methods for steady and unsteady problems, in: K.W. Morton and M.J. Baines, eds., Numerical Methods for Fluid Dynamics (Academic Press, New York, 1982) l32. D. Allen and R. Southwell, Relaxation methods applied to determining the motion, in two dimensions, of a viscous fluid past a fixed cylinder, Q. J. Mech. Appl. Math. 8 (1955) 129145. M.S. Mock, Analysis of a discretisation algorithm for stationary continuity equations in semiconductor device models II, COMPEL 3 (1984) 137149. R.E. Bank, D.J. Rose and W. Fichtner, Numerical Methods for Semiconductor Device Simulation, IEEE Trans. Electron. Devices 30 (9) (1983) 10311041. P.A. Markowich and M. Zlamal, Inverseaveragetype finite element discretisations of selfadjoint second order elliptic problems, Math. Comp. 51 (1988) 431449. J. Douglas, jr., I. Martinez Gamba and M.C.J. Squeff, Simulation of the transient behavior of a onedimensional semiconductor device, Mat. Apl. Comput. 5 (1986) 103122. F. Brezzi, L.D. Marini and P. Pietra, Methodes d’elements finis mixtes et schema de ScharfetterGummel, C.R. Acad. Sci. Paris Serie I 305 (1987) 599604. P.G. Ciarlet, The Finite Element Method for Elliptic Problems (NorthHolland, Amsterdam, 1978). P.A. Raviart and J.M. Thomas, A mixed finite element method for second order elliptic problems, in: Mathematical Aspects of the Finite Element Method, Lecture Notes in Mathematics, Vol. 606 (Springer, Berlin, 1977) 292315. F. Brezzi, On the existence uniqueness and approximation of saddlepoint problems arising from Lagrangian multipliers, RAIRO 8R2 (1974) 129151. B.X. Fraeijs de Veubeke, Displacement and equilibrium models in the finite element method, in: O.C. Zienkiewicz and G. Hollister, eds., Stress Analysis (Wiley, New York, 1965). D.N. Arnold and F. Brezzi, Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, M2AN 19 (1985) 732. F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, submitted to Springer. K. Baba and M. Tabata, On a conservative upwind finite element scheme for convective diffusion equations, RAIRO Anal. Numer. 15 (1981) 325.