Yirang YUAN (袁益讓)

Institute of Mathematics, Shandong University, Jinan 250100, China

E-mail : yryuan@sdu.edu.cn

Changfeng LI (李長(zhǎng)峰)?

Shandong Applied Financial Theory and Policy Research Base, Jinan 250100, China School of Economics, Shandong University, Jinan 250100, China

E-mail : cfli@sdu.edu.cn

Huailing SONG (宋懷玲)

College of Mathematics and Econometrics, Hunan University, Changsha 410082, China

E-mail : shling@hnu.edu.cn

1 Introduction

The numerical simulation of a three-dimensional semiconductor device of heat conduction is a fundamental problem in information science. The mathematical model is defined by four nonlinear partial differential equations with initial-boundary conditions: 1) an elliptic equation for electric potential, 2) a convection-diffusion equation for electron concentration, 3) a convection-diffusion equation for hole concentration,4)a heat conduction equation for temperature. The electric potential appears within the concentration equations and heat equation,and plays an important role in the model. The mathematical model is formulated as follows on a three-dimensional domain ? [1–4]:

The electric potential, electron concentration, hole concentration and temperature are the objective functions, denoted by ψ, e, p and T, respectively. All the coefficients of (1.1)–(1.4) are bounded by two positive constants. α =q/ε, where q and ε are positive constants that denote the electronic load and the permittivity, respectively. UTis the thermal voltage. The diffusion Ds(X) depends on the mobility μs(X), i.e., Ds(X) = UTμs(X), (s = e for the electron and s=p for the hole). ND(X)and NA(X)are the donor impurity concentration and acceptor impurity concentration, respectively. N(X), a difference of ND(X) and NA(X), changes rapidly as X approaches the P-N junction. R1(e,p,T) and R2(e,p,T) are the recombination rates of the electron, hole and temperature. ρ(X) is the heat transfer coefficient. It is important to consider a nonuniform partition in numerical simulation [5, 6].

We have the initial conditions

where e0(X), p0(X) and T0(X) are positive given functions.

In this paper, we concentrate on the Neumann boundary conditions

where ?? is the boundary of ?, and γ is the unit outer normal vector of ??.

A compatibility condition is defined by

and the following condition is given to avoid an ambiguous solution:

The numerical simulation of a semiconductor device is an important topic in the manufacturing of modern semiconductors[5–8]. The sequence iteration is adopted to solve the semiconductor problems introduced by Gummel in 1964, and a new problem of numerical simulation in semiconductor devices was proposed [9]. Douglas and Yuan first put forward a simple and useful finite difference method for two preliminary problems (constant coefficients and without temperature effect), discussed the applications and gave rigorous theoretical results in [10, 11].The characteristic finite element method was discussed for a variable coefficient problem in[12].Since the diffusion only includes the electric-field strength ??ψ, the characteristics-mixed finite element was presented and the optimal-order errors in the H1-norm and the L2-norm were concluded in [13, 14]. Considering actual conditions and the effect of heat, Yuan discussed a characteristic finite difference method for a three-dimensional semiconductor device of heat conduction on a uniform partition [4].

The finite volume element is simple and conservative [15, 16], so it becomes an efficient numerical method for solving partial differential equations. A mixed finite element could obtain the potential and strength simultaneously, and the accuracy is improved one order [17–19]. A block-centered numerical method was discussed in [20, 21]by combining the methods of finite volume element and mixed elements; its computational efficiency was tested experimentally in [22, 23]. Convergence analysis was mainly stated for elliptic problems in [24–26], and the general theoretical discussion was established for the block-centered difference. Rui and Pan used this method to discuss numerical computation for Darcy-Forchheimer flow problems in[27,28]. A satisfactory application of this method was shown for the numerical simulation of a semiconductor device while the computational scheme is not conservative[29,30]. The adaptive finite element method on changing meshes is an effective tool for solving partial differential equations with high accuracy, especially for simulating the local property of the P-N joint of a semiconductor device problem. Dawson and Kirby put forward a mixed finite element method to compute the heat equation and convection-dominated diffusion equation on dynamically changing meshes,and obtained optimal order error estimates on a special changing mesh in[31].

Based on the previous studies, in this paper a block-centered upwind difference method on changing meshes (BCUDMCM) is put forward for a three-dimensional semiconductor device of heat conduction. The potential is computed using a conservative block-centered difference scheme, and the accuracy is improved by one order for electric-field strength. BCUDMCM is applied to compute the concentrations and the temperature, where the diffusion is treated by a block-centered difference and the convection is approximated by an upwind scheme. The upwind method eliminates numerical dispersion and nonphysical oscillation, so it solves the convection-dominated diffusion equations well. The concentrations, temperature and their adjoint vector functions are computed simultaneously. The scheme is conservative locally because piece-wise-defined constant functions are used. This conservative nature is important in numerical simulation of semiconductor device problems. By applying the priori estimates theory,the method of duality, and postprocessing special techniques of changing meshes, we obtain optimal-order error estimates. In this paper, numerical experiments are given for a simplified model to show high computational efficiency and accuracy. This method provides an efficient tool for solving challenging benchmark problems [2, 5–8].

Suppose that the problem of (1.1)–(1.8) is regular:

The coefficients are positive definite:

where D?,D?,μ?,μ?,ρ?and ρ?are positive constants. R1(e,p,T)and R2(e,p,T)are Lipschitz continuous on a ε0-neighborhood of X,i.e.,there exists K >0 such that,for|εi|≤ ε0(1 ≤ i ≤ 6),

In the ensuing discussions, the symbols K and ε denote a generic positive constant and a generic small positive number, respectively; they may take different values at different places.

2 Notation and Preliminaries

Two nonuniform partitions are given to define the BCUDMCM. The coarse partition and the refined one dependent on tnare for the potential and electric-field strength and other unknowns, respectively. The coarse partition is considered first.

For simplicity, take ? ={[0,1]}3, and let ?? denote the boundary. Define

? is partitioned by δx× δy× δz. Here Nx, Nyand Nzare positive integers. Let ?ijk={(x,y,z)|xi?1/2< x < xi+1/2,yj?1/2< y < yj+1/2,zk?1/2< z < zk+1/2} denote an element,and define xi= (xi?1/2+ xi+1/2)/2, hxi= xi+1/2? xi?1/2, hx,i+1/2= xi+1? xi, hx=and define yj, zk, hyj, hzk, hy,j+1/2, hz,k+1/2, hyand hzsimilarly. Let hψ=The partition is supposed to be regular. A simplified case of Nx=4,Ny=3 and Nz=3 is illustrated in Figure 1.

Figure 1 The partition

Define an experimental space by= {f ∈ Cl[0,1]: f|Ii∈ pd(Ii),i = 1,2,··· ,Nx},where Ii=[xi?1/2,xi+1/2],and pd(Ii)denotes a space consisting of all the polynomial functions of a degree of at most d constricted on Ii. The function f(x) is possibly discontinuous on [0,1]if l =?1.are defined similarly. Let

Define the inner products and norms by

and define (v,w)y, (v,w)z,in a similar manner. Let

For a vector function w=(wx,wy,wz)T, define its norms by

The inner product and norm in L2(?)are denoted by(·,·)andFor a function v ∈ Sh,it clearly holds that

Define the difference operators

Then several preliminaries follow for numerical analysis.

Lemma 2.1For v ∈ Shand w ∈ Vh,

Lemma 2.2For w ∈Vh,

ProofThis is equivalent to provingFrom the fact that

Lemma 2.3For w ∈Vh,

ProofThe first termis proven, and the others are shown similarly.Note that

and by the Cauchy inequality, we have

Multiplying both sides by hx,i+1/2hyjhzk, and making the summation, we obtain

Thus, the proof is completed.

Another partition is obtained by refining the coarse partition of ?=[0,1]3uniformly,where the step is usually taken as 1/l times the coarse step; that is, hs= hψ/l for a positive integer l. The refined partition changes with respect to tn.

3 The Block-Centered Upwind Difference Scheme on Changing Meshes

Rewrite (1.1) as

The other equations are rewritten in the divergence forms. Letand zT=??T. Then,

The expanded block-centered difference is discussed in[32,33],which approximates the diffusion flux functions ze, zp, zTand the gradient vectorssimultaneously.

The time interval J =[0, ˉT]is partitioned as follows:

(1) For the potential, ?tψ,1is the first time step and ?tψdenotes the time step during the following computations: J is partitioned by 0 = t0< t1< ··· < tM=and the i-th time level is denoted by ti= ?tψ,1+(i ? 1)?tψ,i ≥ 1.

(2) For the concentrations and temperature, the time interval is partitioned by 0 = t0

(3) Suppose that there exists a positive integer n such that tm= tnfor m, andis a positive integer. Let j0= ?tψ,1/?ts, j = ?tψ/?ts.

The discrete variational form of (3.2) is

The block-centered upwind difference scheme is defined for solving (3.2) by

The value of a discrete solution at tn, tm?1

An upwind approximation is used to compute the flux Ge. The mean value of the integral of·γ is supposed to be zero because of ge= μeeu = 0 on ??. Let σ denote the interface of ω1and ω2. Xlis the barycentre and γlis the unit outer normal vector to ω2. Define

Note that the discrete system is nonsymmetric for (3.7). We could obtain the symmetric formulation

Similarly, the variational form of (3.3) is given by

The block-centered upwind difference scheme of (3.3) is stated as follows:

The variation form of (3.4) is

The block-centered scheme is given for solving (3.4):

We have initial approximations

where the Ritz-projection (its definition is given later), the interpolation, or the L2-projection are used.

BCUDMCM runs as follows: {U0,Ψ0}is obtained by using the initial conditions,the conjugate gradient method and (3.5) together. Then (E1,P1,H1),(E2,P2,H2),··· ,(Ej0,Pjo,Hj0)are computed using the procedures of (3.7), (3.12) and (3.14). Let Sj0+(m?1)j= Sm(S =E,P,H) for m ≥ 1. From (3.5), we can obtain {Um,Ψm}. The computations are repeated for n ≥ 2 and m ≥ 1. The values of{Sj0+(m?1)j+1}, {Sj0+(m?1)j+2},···, {Sj0+mj}(S =E,P,H)are obtained. Finally, we can get all the numerical solutions. According to (C), the numerical solutions exist and are unique.

Remark 3.1Since the mesh changes,the term(En?1,vn)ˉmin(3.7a)should be computed,i.e., the L2?projection of En?1should be given. The piecewise-defined constant function onis mapped onto a piecewise defined constant onThe termsin (3.12a) andin (3.14a) are computed similarly.

4 The Principle of Mass Conservation

Suppose that the recombination rates are zero(i.e.,R1(e,p,T)≡0 and R2(e,p,T)≡0)and suppose that the Neumann boundary condition is homogeneous. When the mesh is unchanging at a time level tn, the local conservation of mass for concentrations holds on each element ω = [xi?1/2,xi+1/2]× [yj?1/2,yj+1/2]× [zk?1/2,zk+1/2]. The conservative nature is expressed by

where ω is an element of the partition, ?ω is the boundary, and γ?ωis the unit outer normal vector. The computational scheme of (3.7a) and (3.12a) is as follows:

Theorem 4.1If Ri(e,p,T)≡0,i=1,2, then

ProofWe just prove (4.2a). (4.2b) can be shown in a similar manner. The value of a function v ∈ Shis assigned by 1 on ω = ?ijkand 0 on other elements; i.e.,

Then (3.7a) is rewritten as

Using the notation in Section 2, we have

Substituting (4.4) into (4.3), we complete the proof.

Furthermore, we conclude the conservation on the whole domain.

Theorem 4.2If Ri(e,p,T)≡0,i = 1,2 and the homogeneous Neumann boundary condition holds, then (3.5a) and (3.7a) have the following nature:

Proof(4.5a) is proven first. (4.5b) can be proven similarly. Summing (4.2a) on all the elements, we have

Here ω1intersects ω2at σ. Xlis the barycentre and γlis the unit outer normal vector to ω2.Recalling the definition of convection flux, we have, when EUn·γl(X)≥ 0 on ω1,

where|σl|is the measure of σl. On ω2,the unit normal vector to σlis ?γl,and EUn·(?γl(X))≤0. Then,

The signs of (4.7a) and (4.7b) are contrary, therefore,

Substituting (4.8) and (4.9) into (4.6), we complete the proof of (4.5a).

The conservation of mass is an important physical element in the numerical simulation of semiconductor devices.

5 Convergence Analysis

Elliptic projections are introduced first. Define

where e and p are exact solutions of (1.1)–(1.8).

Lemma 5.1The coefficients and exact solutions of (1.1)–(1.8) are supposed to satisfy(C) and (R), so there exist two positive constantsindependent of h and ?t such that

We estimate π and σ first. By subtracting (5.1a) (t=tm) and (5.1b) (t=tm) from (3.5a)and (3.5b), respectively, we have

Taking v = πmin (5.6a) and w= σmin (5.6b), we have the combination equation

Applying Lemma 2.1–Lemma 5.1, we have

A duality method is introduced to address πm∈ Sh[32, 34]. Consider the following elliptic problem:

It follows from the regularity of (5.9) that

By Lemma 2.3, (5.9), (5.10) and (5.8), we have

Substituting (5.13) into (5.12) and considering (5.8), we have

Now we discuss the concentration equation (1.2). Subtracting (5.2) at t=tnfrom (3.7),

The terms on the left-hand side of (5.16) are estimated as follows:

The right-hand terms are considered. The first one is estimated by

Let σ denote an interface in the partition. γris the outer normal vector to σ, and Xlis the barycenter. Then,

Using the regularity (R) and the mean value theorem for integrals, we have

From the regularity of en, Lemma 5.1 and the discussions in [21, 32, 33],

It follows from (5.18)-(5.21)that

Then we have

Substituting (5.17)–(5.21) into (5.16) gives

Multiplying both sides of (5.24) by ?ts, summing on n, then usingwe have

In a similar manner, (1.3) is estimated by

Finally, the heat conduction equation (1.4) is demonstrated. Subtracting (5.4a) (t = tn)and (5.4b) (t=tn) from (3.14a) and (3.14b), respectively, we obtain

An induction hypothesis is introduced:

Next, the partition parameters are supposed in order to satisfy

Estimating both sides of (5.30) and using Lemma 2.1–Lemma 5.1 and (5.28), we have

It follows from (5.25), (5.26) and (5.31) that

The last right-hand term of (5.32) is considered under two different cases: that of a fixed partition and that of changing meshes. The sum of electron concentration is divided into two parts:

The unchanging partition case is denoted by n=n′, so

Let n=n′′denote another case, so

Assume that the mesh changes at most ˉM times,and ˉM ≤M?,where M?is a positive constant independent of h and ?t. Then, from (5.33) and Lemma 5.1,

The terms of the hole and temperature are treated similarly. Substituting (5.34) into (5.32)and using the discrete Gronwall Lemma and initial approximations, we have

Combining (5.14) with (5.35), we obtain

Finally, the assumption (5.28) is verified. Using (5.36) and (5.29), we can prove (5.28)easily.

The following theorem is concluded from (5.35), (5.36) and Lemma 5.1:

Theorem 5.2Suppose that exact solutions of (1.1)–(1.8) are regular (R), and that the coefficients are positive definite (C). Numerical solutions are computed by using the scheme of(3.5), (3.7), (3.12) and (3.14). The partition is supposed to satisfy (5.29). Then,

6 An Improved Scheme

In the discussions of Section 3 and Section 5,an assumption is introduced that M?is not a large positive integer,i.e., the meshes could not change too many times (not greater than M?).From theorem 5.2, we see that the computational accuracy depends on M?hs. If the meshes change frequently (M?approaches infinity), the computational accuracy is possibly decreased.Thus, we improve the scheme of (3.5), (3.7), (3.12) and (3.14)and ignore the assumption. The computational accuracy of the improved scheme is independent of the changing times. The electric potential is solved by (3.5). The concentrations and temperature are solved by an improved scheme of (3.7), (3.12) and (3.14). A linear approximation is used to replace the L2-projection, and an optimal error estimate is obtained on changing meshes.

The electron concentration equation is considered first. Givena linear approximationon ωn?1is defined by

(3.7b) and (3.7c) are still used, while (3.7a) is improved by

Similarly, we have error (5.15b) and (5.15c). (5.15a) is replaced by

where Π is the mapping operator of (5.2),

The first three terms on the right-hand side of(6.5)are estimated as above. Here we discuss the last two terms. Note that

Recalling the definitions of (6.1), (6.2) and (6.6), at any time level tn, we have

By using (5.29), it obviously holds that hs=O(?ts). Substituting (6.8) into (6.7), we obtain

The last term is considered. Note that

Using the Taylor expansion, for ? X ∈ ωn?1, we have

Substituting (6.11) into (6.10), we get

Collecting (6.9), (6.12) and (6.5), multiplying both sides of the resulting inequality by ?ts,summing on n=1,2,··· ,N, and usingwe obtain

Similarly, we get the estimates of the hole concentration:

Under the assumption (5.28) and the restriction condition (5.29), we obtain error estimates of the heat equation in a similar manner:

Combining error estimates (6.13)-(6.15) and using the discrete Gronwall’s lemma and initial approximations, we have

It is derived by using (5.14) and (6.16) that

The induction hypothesis (5.28) can be tested easily.

From (6.16), (6.17) and Lemma 5.1, a convergence result is concluded.

Theorem 6.1Suppose that the conditions (C) and (R) hold and that the partition satisfies the restriction condition (5.29). The problem of (1.1)–(1.8) is solved by the improved scheme. Then,

where M???is a positive constant only dependent on ψ,u,e,p,T and their derivatives.

7 Numerical Examples

Here we consider experimental tests for the simplified model. Suppose that the electric potential and the strength are known. Consider a convection-dominated diffusion equation for the concentration to show the application of the block-centered upwind difference method:

where a(x) = 1.0 × 10?4(1+x3), b(x) = e1+2x(?0.2x2+2x ? 1), t ∈ (0,0.1], x ∈ [0,1]. The exact solution is defined by c=exp(?0.2t)(sin(πx2? t))8with the proper f.

The concentration function has one sharp front during the interval[0,1], and changes with respect to the variable t. The finite difference method or finite element method has a numerical oscillation. A block-centered upwind difference is used for this problem and numerical results are shown in Figure 2. That the upwind method on changing meshes solves the problem well is seen in Figure 3. Numerical results of finite difference, upwind finite difference and upwind difference on changing meshes are compared in Figure 4, Figure 5 and Figure 6.

Figure 2 BCUFDM solution at t=0.1

Figure 3 BCUFDCM solution

Figure 4 UFD solution

Figure 5 UFDCM solution

Figure 6 FDM solution

During the computations, let h, ?t and hf= h/k denote the space step, time step and refined space step, respectively. Take ?t=2h. The upwind finite difference, upwind finite difference on changing meshes with a positive integer k,block-centered upwind finite difference and block-centered upwind finite difference on changing meshes are denoted by UFD, UFDCM(k),BCUFD and BCUFDCM(k),respectively. AL2and RL2denote the absolute and relative errors at t=0.1 in the L2-norm.

From Tables 1 and 2, we see that UFD and UFDCM(k) have the same computational accuracy, and BCUFD, BCUFDCM(k) have a similar accuracy.

Table 1 Absolute errors, AL2×103

Table 2 Relative errors, RL2×103

Furthermore, the block-centered upwind approximation can compute c and its derivativesimultaneously. The derivatives are approximated by cx,h. The graphical comparisons are illustrated in Figure 7 and Figure 8,and errors are shown in Table 3 and Table 4. Let

Figure 7 BCUFD solution for cx

Figure 8 BCUFDCM solution for cx

Table 3 Absolute errors for cx, AELn2 ×103 N BCUFD BCUFDCM (2) BCUFDCM (3) BCUFDCM (5)100 3.6237 2.0690 1.6836 1.5106 200 1.8391 1.1649 1.1166 1.0397 500 0.7717 0.4459 0.4033 0.3814 1000 0.3991 0.2158 0.1586 0.1173 2000 0.2043 0.1078 0.07548 0.04982 3000 0.1375 0.07235 0.05042 0.03287

L2 ×103 N BCUFD BCUFDCM (2) BCUFDCM (3) BCUFDCM (5)100 90.6742 51.7662 42.0658 37.7801 200 46.2665 28.9277 27.7108 23.7786 500 19.3460 11.0923 10.0529 9.4914 1000 9.9841 5.3880 3.9581 2.9259 2000 5.1037 2.6900 1.8824 1.2420 3000 3.4321 1.8045 1.2571 0.8194 Table 4 Relative errors for cx, REn

From Table 1, numerical approximations at h=1/1000 and at h=1/200, k =5 have the same accuracy. Their running time (denoted by Time) comparisons are seen in Table 5 (unit:s).

Table 5 Running time comparisons

It is easy to see that BCUFDCM has a high computational accuracy,shortens running time and improves the computational efficiency.

Next, a two dimensional example is considered. Suppose that the electric field strength is given:

Let ? =[0,1]× [0,1]. The exact solution is defined by c=exp(?0.2t)(sin(πx ? t)sin(πy ? t))8.Let u=((?0.2x2+2x ? 1)e1+2x,(?0.2y2+2y ? 1)e1+2y), dm=1.0 × 10?4(1+x2+y2), and D(u) = dmI. This convection-dominated problem is solved by the UFD, UFDCM, BCUFD and BCUFDCM, and their errors are compared in Table 6. The partitions are uniform with h=1/N, ?t=2h and hf=h/k. Take=0.1.

Table 6 Errors comparison, AL2×103

From Table 6, the BCUFD and UFD have similar accuracy under the same partitions.Furthermore, BCUFD can compute the concentration and its gradient operator at the same time. The absolute errors and relative errors of ?D?c are given in Tables 7 and 8. The computational accuracy is improved by one order, consistent with theoretical results.

Table 7 Absolute errors for ?D?c, AEL2 × 106

Table 8 Relative errors for ?D?c, REL2 × 103

Based on numerical examples, it is concluded that the present method can solve the convection-dominated diffusion equation well,and so could solve more complicated actual problems.

8 Conclusions and Discussions

Numerical simulation of the three dimensional semiconductor device transient behavior problem of heat conduction is discussed in this paper. A block-centered upwind finite difference scheme is constructed on changing meshes and numerical analysis are presented. In Section 1, the mathematical model is stated, and the physical background and related research are introduced. In Section 2, the notation and preliminary theoretical results are stated, and two different(coarse and fine)partitions are defined. In Sections 3 and 4,we propose the procedure of a block-centered upwind difference method on changing meshes. The electric potential is solved by a conservative block-centered method, and an approximation of electric field strength with one-order accuracy improvement is shown. A block-centered upwind difference scheme is applied to solve the concentration equations and the heat conduction equation on changing meshes. The diffusion and convection are solved by the block-centered scheme and the upwind scheme,respectively. The composite scheme can solve convection-dominated diffusion problems well because it avoids numerical dispersion and nonphysical oscillation. The concentrations and the temperature and their adjoint vector functions are computed simultaneously. The scheme has an element of conservation,which is important in the numerical simulation of semiconductor devices. In Section 5, we give a theoretical analysis using a priori estimates of differential equations and special techniques to derive the optimal-order error estimates. In Section 6, an improved block-centered upwind difference scheme is addressed on changing meshes, and the optimal order convergence result is given. In Section 7, numerical experiments are illustrated to show the efficiency and application of this method, and the challenging problem [2, 5–8]is solved.

Several interesting conclusions are obtained.

(I) The method has the physical element of conservation, which is important in the numerical simulation of semiconductor devices. This improves and generalizes the research on numerical simulation in information science [2, 5–10].

(II)The method combines the block-centered difference,upwind approximation and changing meshes, so it has strong stability and high accuracy and is especially useful in large-scale engineering computations of complicated three-dimensional regions.

(III) Using special postprocessing, this paper provides optimal error estimates. Thus, an improvement is made on the changing meshes for solving parabolic problems in [31].

(IV) Numerical examples illustrate the feasibility and application of the method presented,and give numerical tests for showing how the sharp fronts change with different times. This is important in the simulation of semiconductor devices [2, 5–10].