Acta Mathematicae Applicatae Sinica, English Series Vol. 33, No. 4 (2017) 1053–1072 DOI: 10.1007/s10255-017-0721-y http://www.ApplMath.com.cn & www.SpringerLink.com

Acta Mathemacae Applicatae Sinica, English Series The Editorial Office of AMAS & Springer-Verlag Berlin Heidelberg 2017

A Mixed-finite Volume Element Coupled with the Method of Characteristic Fractional Step Difference for Simulating Transient Behavior of Semiconductor Device of Heat Conductor And Its Numerical Analysis Yi-rang YUAN1,† , Qing YANG2 , Chang-feng LI1,3 , Tong-jun SUN1 1 Institute

of Mathematics, Shandong University, Jinan 250100, China († E-mail: [email protected])

2 School

of Mathematical Sciences, Shandong Normal University, Jinan 250014, China

3 School

of Economics, Shandong University, Jinan 250100, China

Abstract The mathematical system is formulated by four partial differential equations combined with initialboundary value conditions to describe transient behavior of three-dimensional semiconductor device with heat conduction. The first equation of an elliptic type is defined with respect to the electric potential, the successive two equations of convection dominated diffusion type are given to define the electron concentration and the hole concentration, and the fourth equation of heat conductor is for the temperature. The electric potential appears in the equations of electron concentration, hole concentration and the temperature in the formation of the intensity. A mass conservative numerical approximation of the electric potential is presented by using the mixed finite volume element, and the accuracy of computation of the electric intensity is improved one order. The method of characteristic fractional step difference is applied to discretize the other three equations, where the hyperbolic terms are approximated by a difference quotient in the characteristics and the diffusion terms are discretized by the method of fractional step difference. The computation of three-dimensional problem works efficiently by dividing it into three one-dimensional subproblems and every subproblem is solved by the method of speedup in parallel. Using a pair of different grids (coarse partition and refined partition), piecewise threefold quadratic interpolation, variation theory, multiplicative commutation rule of differential operators, mathematical induction and priori estimates theory and special technique of differential equations, we derive an optimal second order estimate in L2 -norm. This numerical method is valuable in the simulation of semiconductor device theoretically and actually, and gives a powerful tool to solve the international problem presented by J. Douglas, Jr. Keywords

transient behavior of three-dimensional semiconductor device; numerical simulation; mixed finite

volume element; modified characteristic fractional step difference; second order estimate in L2 norm 2000 MR Subject Classification

1

65M015; 65M60; 65N30; 82D37

Introduction

Traditional numerical methods are invalid for precise simulation of the problems from modern development and technical applications of semiconductor device, since a complicated diffusion system of quasi-linear partial differential equations with initial boundary value conditions is involved. Special for a multi-dimensional problem, modern numerical simulation technique is Manuscript received April 27, 2015. Revised January 19, 2016. Project supported by the National Natural Science Foundation of China (Grant Nos. 11101124 and 11271231), the National Tackling Key Problems Program for Science and Technology (Grant No. 20050200069), and the Doctorate Foundation of the Ministry of Education of China (Grant No. 20030422047). † Corresponding author.

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1054

aimed to find new computational schemes of high order accuracy over more irregular geometric regions[8,16,33,34]. The mathematical description of three-dimensional semiconductor device problem is formulated by four quasi-linear partial differential equations[1,6,8,9,16,28,29,33,34], including an elliptic equation for the electric field potential, two convection diffusion equations for the concentrations of electron and hole and a heat conductor equation for the temperature, combined with initial-boundary conditions. On a three-dimensional region, the whole system is stated as follows, X = (x, y, z)T ∈ Ω, t ∈ J = (0, T], (1) − Δψ = ∇ · u = α(p − e + N (X)),   ∂e = ∇ · De (X)∇e + μe (X)eu − Re (e, p, T ), (X, t) ∈ Ω × J, (2) ∂t   ∂p = ∇ · Dp (X)∇p − μp (X)pu − Rp (e, p, T ), (X, t) ∈ Ω × J, (3) ∂t   ∂T − ΔT = (De (X)∇e + μe (X)eu) − (Dp (X)∇p − μp (X)pu) · u, (X, t) ∈ Ω × J. ρ ∂t (4) In the above expressions, the electric potential ψ, the electron concentration e, the hole concentration p, the temperature T , and the electric field intensity u = −∇ψ are objective functions. The electric potential is determined by the intensity and involved in the latter three equations and governs the concentrations and the heat directly. α = q/ε is defined by the quotient of two positive constants, the electronic load q and the dielectric coefficient ε. All coefficients in (2)–(4) have a positive upper bound and a positive lower bound. The diffusion coefficients De (X) and Dp (X), are equal to UT μe (X) and UT μp (X), the products of the mobilities μe (X) or μp (X) and the heat UT . A given function N (X) is defined by ND (X)−NA(X), the difference of impurity concentrations of the donor and the acceptor, whose values change rapidly as X approaches closely the semiconductor P-N junction. The symbols Re (e, p, T ) and Rp (e, p, T ) denote generation-recombination rates of the electron, the hole and the temperature. The coefficient of heat conduction ρ(X) is positive definite. Some necessary conditions are stated to determine the solutions. Initial value conditions: e(X, 0) = e0 (X),

p(X, 0) = p0 (X),

T (X, 0) = T0 (X),

where e0 (X), p0 (X) and T0 (X) are known and positive functions. Neumann boundary value conditions:     ∂ψ  ∂e  ∂p  ∂T  = = = = 0, t ∈ J, ∂γ ∂Ω ∂γ ∂Ω ∂γ ∂Ω ∂γ ∂Ω where γ denotes outer normal vector to the bound ∂Ω of Ω. The compatibility condition:  [p − e + N ]dX = 0, t ∈ J.

X ∈ Ω,

(5)

(6)

(7)

Ω

The uniqueness condition:

 ψdX = 0, Ω

t ∈ J.

(8)

Numerical simulation of semiconductor device originates from a sequence of iterations put forward by Gummel in [7]. Douglas and Yuan presented a kind of applicable difference method

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

1055

for 1-D and 2-D simple cases with constant coefficients and temperature effect being ignored, and gave the theoretical analysis first[6,23] . Then Yuan discussed the methods of characteristics coupled with finite element or mixed finite element for two-dimensional problems with varied coefficients and derived optimal order error estimates in H 1 norm and L2 norm[24,25,28] . To discuss numerical simulation more effectively and more practically, the authors considered the effect of temperature in three dimensional case[12,17,29,33,34] , proposed the fractional step methods combined with characteristic difference and upwind difference, and obtained second order error estimate in L2 norm[30−32] . By consideration of mass conservation in numerical simulation, a combination of mixed finite element method and method of characteristics was discussed by the authors in [4,5,18,20–22], where the potential equation is approximated by a normal lower degree mixed finite element scheme and the convection-dominated diffusion equations in the concentrations are solved by the mixed finite element method and the characteristic method. Successively, another conservative algorithm combining the method of mixed finite volume element and the method of characteristic fractional step difference was studied[11,15,19] . In this paper, the mixed finite volume element method which is conservative in mass and practical computations is used to approximate the potential and makes the potential computation one order accurate. The modified characteristic fractional step difference method is applied to approximate the values of electron concentration, hole concentration and temperature, where a multi-dimensional problem is divided into a sequence of one-dimensional subproblems and the computation is carried out easily by the method of speedup. Therefore, the parallel program developed by the method can shorten the computation time, decrease the complexity and could be used in actual simulation of high order accuracy. In theoretical analysis, by using a coupled partition of coarse grids and refined grids, piecewise threefold quadratic interpolation of product type, calculus of variations, multiplicative operators communication rule of differential type, induction hypothesis and a priori error estimate theory, and other techniques of differential equations, we present an optimal order error estimate in l2 norm. The method discussed in this paper is significant in theoretical analysis and application in model analysis, numerical method and mechanism research as well as software design[1,6−9,16,23,33,34], and also, this method can solve the international outstanding reservoir problem proposed by Jr. J. Douglas. The problem is supposed to be positive definite, 0 < D∗ ≤ Ds (X) ≤ D∗ , s = p, e;

(C)

0 < ρ∗ ≤ ρ(X) ≤ ρ∗ ,

(9)

where D∗ , D∗ , ρ∗ and ρ∗ are positive constants. Exact solutions are assumed to satisfy the regularity, (R)

ψ ∈ L∞ (J; H 3 (Ω)) ∩ H 1 (J; W 4,∞ (Ω)), 2

2

e, p, T ∈ L∞ (W 4,∞ (Ω)),

2

∂ e ∂ p ∂ T , , ∈ L∞ (L∞ (Ω)). ∂t2 ∂t2 ∂t2

(10)

And the recombination rates Re (p, e, T ), Rp (p, e, T ) are Lipschitz continuous with respect to p, e, T in ε0 -domain centered by (p, e, T ), that is to say there exists a positive constant M , for |εi | ≤ ε0 (1 ≤ i ≤ 6), such that   Rs (p(X, t) + ε1 , e(X, t) + ε2 , T (X, t) + ε3 ) − Rs (p(X, t) + ε4 , e(X, t) + ε5 , T (X, t) + ε6 )   ≤M |ε1 − ε4 | + |ε2 − ε5 | + |ε3 − ε6 | , (X, t) ∈ Ω × J. (11) All the above assumptions are reasonable for physical consideration[1,8,9,16]. In this paper, M and ε denote a generic positive constant and a generic positive small constant, respectively, and they have different meanings at different places.

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1056

2

Notations and Preparations

To structure the algorithm by using mixed finite volume element and characteristic fractional step difference, we define two different grids, a coarse partition and a refined partition. The former is nonuniform and the program for computing the electric potential and the intensity run under a large time step. The latter is uniform and the size of the time step for the computation of the concentrations and heat is relatively small. The coarse partition is discussed first. For convenience to illustrate the construction of the algorithm, we consider the threedimensional domain Ω = {[0, 1]}3, and discretize it as: Ixh : 0 < x1/2 < x3/2 < · · · < xNx −1/2 < xNx +1/2 = 1, Iyh : 0 < y1/2 < y3/2 < · · · < yNy −1/2 < yNy +1/2 = 1, Izh : 0 < z1/2 < z3/2 < · · · < zNz −1/2 < zNz +1/2 = 1. The computational domain Ω is divided by Ixh × Iyh × Izh , and for i = 1, 2, · · · , Nx , j = 1, 2, · · · , Ny , k = 1, 2, · · · , Nz , the subdomain corresponding to a fixed point (i, j, k) is defined by Ω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 }. Some other notations are defined as follows, xi = (xi−1/2 + xi+1/2 )/2, yj = (yj−1/2 + yj+1/2 )/2, zk = (zk−1/2 + zk+1/2 )/2, hxi = xi+1/2 − xi−1/2 , hyj = yj+1/2 − yj−1/2 , hzk = zk+1/2 − zk−1/2 , hx,i+1/2 = (hxi +hxi+1 )/2 = xi+1 −xi , hy,j+1/2 = (hyj +hyj+1 )/2 = yj+1 −yj , hz,k+1/2 = (hzk + hzk+1 )/2 = zk+1 − zk , hx = max {hxi }, hy = max {hyj }, hz = max {hzk }, and hψ = 1≤i≤Nx

1≤j≤Ny

1≤k≤Nz

(h2x + h2y + h2z )1/2 . The partition is regular if there exist α1 , α2 > 0, such that α1 hx ,

min {hyj } ≥ α1 hy ,

1≤j≤Ny

min {hxi } ≥

1≤i≤Nx

min {hzk } ≥ α1 hz , min{hx , hy , hz } ≥ α2 max{hx , hy , hz }.

1≤k≤Nz

1 z7/2 z3 z5/2 z2 z3/2 z1

1=y7/2 y3 y5/2

y2 y3/2 y1 0=y1/2

0 z1/2 0 x1/2 x1

x3/2

x2

x5/2

x3

x7/2 x4 x9/2 1

Fig.1. Sketch of Nonuniform Partition

A simplified partition of Nx = 4, Ny = 3, Nz = 3 is shown in Fig. 1. Define a space of test functions by Mld (δx ) = {f ∈ C l [0, 1] : f |Ωi ∈ pd (Ωi ), i = 1, 2, · · · , Nx }, where pd (Ωi ) denotes a space consisting of all the polynomial functions of degree at most d restricted on Ωi = [xi−1/2 , xi+1/2 ]. The function f (x) is possibly discontinuous on [0, 1] as l = −1. The 0 0 0 (δx ) ⊗ M−1 (δy ) ⊗ M−1 (δz ), spaces Mld (δy ) and Mld (δz ) are defined similarly. Let Sh = M−1 x y z x 1 0 0 y 0 1 Vh = {w|w = (w , w , w ), w ∈ M0 (δx ) ⊗ M−1 (δy ) ⊗ M−1 (δz ), w ∈ M−1 (δx ) ⊗ M0 (δy ) ⊗ 0 0 0 M−1 (δz ), wz ∈ M−1 (δx ) ⊗ M−1 (δy ) ⊗ M01 (δz ), w · γ|∂Ω = 0}. For a grid function v(x, y, z), vijk ,

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

1057

vi+1/2,jk , vi,j+1/2,k and vij,k+1/2 denote the values of v(x, y, z) at different points (xi , yj , zk ), (xi+1/2 , yj , zk ), (xi , yj+1/2 , zk ) and (xi , yj , zk+1/2 ), respectively. The inner products and norms are defined by Ny Nz Nx

(v, w)m =

hxi hyj hzk vijk wijk ,

i=1 j=1 k=1

(v, w)x =

Ny Nz Nx

hxi−1/2 hyj hzk vi−1/2,jk wi−1/2,jk ,

i=1 j=1 k=1

(v, w)y =

Ny Nz Nx

hxi hyj−1/2 hzk vi,j−1/2,k wi,j−1/2,k ,

i=1 j=1 k=1

(v, w)z =

Ny Nz Nx

hxi hyj hzk−1/2 vij,k−1/2 wij,k−1/2 ,

i=1 j=1 k=1

||v||2s = (v, v)s , ||v||∞(x) = ||v||∞(z) =

s = m, x, y, z,

||v||∞ =

max

1≤i≤Nx ,1≤j≤Ny ,1≤k≤Nz

max

|vi−1/2,jk |,

max

|vij,k−1/2 |.

1≤i≤Nx ,1≤j≤Ny ,1≤k≤Nz 1≤i≤Nx ,1≤j≤Ny ,1≤k≤Nz

||v||∞(y) =

|vijk |, max

1≤i≤Nx ,1≤j≤Ny ,1≤k≤Nz

|vi,j−1/2,k |,

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

, |||w||| = ||wx ||2x + ||wy ||2y + ||wz ||2z 

x 2 1/2 ||w||m = ||w ||m + ||wy ||2m + ||wz ||2m ,

|||w|||∞ = ||wx ||∞(x) + ||wy ||∞(y) + ||wz ||∞(z) , ||w||∞ = ||wx ||∞ + ||wy ||∞ + ||wz ||∞ .

n

∂ v p Define Wpm (Ω) = {v ∈ Lp (Ω)| ∂xn−l−r ∂y l ∂z r ∈ L (Ω), n − l − r ≥ 0, l = 0, 1, · · · , n; r = 0, 1, · · · , n, n = 0, 1, · · · , m; 0 < p < ∞}, and let W2m (Ω) denoted by H m (Ω). The inner product and norm in L2 (Ω) are denoted by (·, ·) and || · ||, and for every function v ∈ Sh , it holds obviously (12) ||v||m = ||v||.

Introduce the difference operators and other notations as follows, vi+1,jk − vijk vi,j+1,k − vijk , [dy v]i,j+1/2,k = , hx,i+1/2 hy,j+1/2 wi+1/2,jk − wi−1/2,jk vij,k+1 − vijk , [Dx w]ijk = , [dz v]ij,k+1/2 = hz,k+1/2 hx i wi,j+1/2,k − wi,j−1/2,k wij,k+1/2 − wij,k−1/2 , [Dz w]ijk = ; [Dy w]ijk = hy j hz k [dx v]i+1/2,jk =

x w ijk = z = w ijk

x x + wi−1/2,jk wi+1/2,jk

2 z z wij,k+1/2 + wij,k−1/2 2

,

y = w ijk

y y + wi,j−1/2,k wi,j+1/2,k

2

,

,

y x z and w  ijk = (w ijk ,w ijk ,w ijk )T . Noting that the values of electric potential vary slowly with respect to t, we adopt a large step Δtψ but a small step Δts for the concentrations. Let Δts = T /L for a positive integer L, j = Δtψ /Δts , tn = nΔts , tm = mΔtψ , ψ n = ψ(tn ) and

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1058

ψm = ψ(tm ), then we can get the value of v(x, y, z, t) at the values of ψ at the closest two time levels to tn , ⎧ ⎪ ⎨ ψ0 , n Eψ = (1 + r/j)ψm − (r/j)ψm−1 , ⎪ ⎩ n t = tm + rΔts ,

t = tn by an linear extrapolation of tn ≤ t 1 tm < tn < tm+1 , r = 1, 2, · · · , j.

(13)

The notations labeled by superscripts and subscripts are assigned for the potential and the concentrations, respectively. Let dt v n = (v n − v n−1 )/Δts denote the difference quotient of v(x, y, z, t) with respect to t. The following three statements are given for numerical analysis. For every pair of functions v ∈ Sh and w ∈ Vh ,

Lemma 1.

(v, Dx wx )m = −(dx v, wx )x ,

(v, Dy wy )m = −(dy v, wy )y ,

(v, Dz wz )m = −(dz v, wz )z . (14)

For every function w ∈ Vh ,

Lemma 2.

|| w ||m ≤ |||w|||.

(15)

Proof. It is adequate to prove ||w x ||m ≤ ||wx ||x , ||w y ||m ≤ ||wy ||y and ||w z ||m ≤ ||wz ||z . From the fact that Ny Nz Nx

x hxi hyj hzk (w ijk )2

≤

i=1 j=1 k=1

Ny Nz

hy j hz k

j=1 k=1

=

Ny Nz

=

hy j hz k

Nx  h

x,i−1

2

i=2

hy j hz k

j=1 k=1

=

2

i=1

j=1 k=1 Ny Nz

Nx x x (wi+1/2,jk )2 + (wi−1/2,jk )2

Ny Nz Nx

x (wi−1/2,jk )2 +

2

i=2

Nx hx i=1

Nx hx,i−1 + hx

i

hx i

2

i

x (wi−1/2,jk )2



x )2 (wi−1/2,jk

x hx,i−1/2 hyj hzk (wi−1/2,jk )2 .

(16)

i=1 j=1 k=1

we have ||w x ||m ≤ ||wx ||x , and get error estimates of the other two terms. Lemma 3.

For a vector w ∈ Vh , ||wx ||x ≤ ||Dx wx ||m ,

Proof.

2

||wy ||y ≤ ||Dy wy ||m ,

||wz ||z ≤ ||Dz wz ||m .

The proof is finished by an illustration of ||wx ||x ≤ ||Dx wx ||m . By the fact that x = wl+1/2,jk

l

x x (wi+1/2,jk − wi−1/2,jk )=

i=1

l x x wi+1/2,jk − wi−1/2,jk i=1

hx i

and by Cauchy inequality, we have x (wl+1/2,jk )2 ≤ xl

Nx i=1

2 hxi [Dx wx ]ijk .

1/2 h1/2 x i hx i ,

(17)

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

1059

Multiplying both sides by hx,i+1/2 hyj hzk , and summing them up, we have Ny Nz Nx

x (wi−1/2,jk )2 hx,i−1/2 hyj hzk ≤

i=1 j=1 k=1

Ny Nz Nx

2 [Dx wx ]ijk hxi hyj hzk .

(18)

i=1 j=1 k=1

Therefore, the proof of Lemma 3 ends. The refined grids for Ω = {[0, 1]}3 are shown by three different partitions:

2

h

I x : 0 = x0 < x1 < x2 < · · · < xM1 −1 < xM1 = 1, h

I y : 0 = y0 < y1 < y2 < · · · < yM2 −1 < yM2 = 1, h

I z : 0 = z0 < z1 < z2 < · · · < zM3 −1 < zM3 = 1, where M1 , M2 and M3 are positive constants. Three steps in different directions are defined by h1 = 1/M1 , h2 = 1/M2 and h3 = 1/M3 and grids are denoted by xi = ih1 , yj = jh2 and zk = kh3 . Let hs = (h21 + h22 + h23 )1/2 and let Di+1/2,jk = 12 [D(Xijk ) + D(Xi+1,jk )], Di−1/2,jk = 12 [D(Xijk ) + D(Xi−1,jk )], then we can give the definitions of Di,j+1/2,k , Di,j−1/2,k , Dij,k+1/2 and Dij,k−1/2 . Let n n n n δx (Dδx W )nijk = h−2 1 [Di+1/2,jk (Wi+1,jk − Wijk ) − Di−1/2,jk (Wijk − Wi−1,jk )],

δy (Dδy W )nijk

=

n h−2 2 [Di,j+1/2,k (Wi,j+1,k

−

n Wijk )

−

n Di,j−1/2,k (Wijk

(19b)

n n n δz (Dδz W )nijk = h−2 3 [Dij,k+1/2 (Wij,k+1 − Wijk ) − Dij,k−1/2 (Wijk n n n n ∇h (D∇h W )ijk = δx (Dδx W )ijk + δy (Dδy W )ijk + δz (Dδz W )ijk .

3

−

(19a)

n Wi,j−1,k )],

−

n Wij,k−1 )],

(19c) (20)

The Composite Procedures of Mixed Finite Volume Element with Second Order Characteristic Fractional Step Difference

For convenience to apply mixed finite volume element, we give an equivalent formulation of (1), ∇ · u = α(p − e + N ), u = −∇ψ,

X ∈ Ω, t ∈ J,

X ∈ Ω, t ∈ J.

(21a) (21b)

Let Z, U, E, P and H denote numerical approximations of ψ, u, e, p and T , respectively. The electric potential and tensity are approximated by the method of mixed finite volume element as follows[11,15,19] , [Dx U x ]m,ijk + [Dy U y ]m,ijk + [Dz U z ]m,ijk = α[Pm,ijk − Em,ijk + Nijk ], y x = −[dx Z]m,i−1/2,jk , Um,i,j−1/2,k = −[dy Z]m,i,j−1/2,k , Um,i−1/2,jk

(22a)

z Um,ij,k−1/2 = −[dz Z]m,ij,k−1/2 .

(22b)

Noting that the flow moves along the direction of characteristics, we apply the method of characteristics to discretize one-order hyperbolic terms of (2) and (3), and get numerical solutions more exactly [3,4,7−9] . Equations (2) and (3) are reformulated as follows, ∂e = ∇ · (De ∇e) + μe u · ∇e + eu · ∇μe + αμe (X)e(p − e + N (X)) − R(e, p, T ), ∂t ∂p = ∇ · (Dp ∇p) − μp u · ∇p − pu · ∇μp − αμp (X)p(p − e + N (X)) − R(e, p, T ). ∂t

(23) (24)

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1060

Let τe = τe (X, t) represent a unit vector of the characteristics (−μe u1 , −μe u2 , −μe u3 , 1), and 1/2  let τp = τp (X, t) be the unit vector of (μp u1 , μp u2 , μp u3 , 1). Given Φs = 1 + μ2s |u|2 , s = e, p, then the directional derivatives in the characteristics are defined by Φe

∂ ∂ − μe u · ∇, = ∂τe ∂t

Φp

∂ ∂ + μp u · ∇, = ∂τp ∂t

(25)

by which equations (23) and (24) are transformed into ∂e − ∇ · (De ∇e) = αμe (X)e(p − e + N (X)) + eu · ∇μe − Re (e, p, T ), ∂τe ∂p − ∇ · (Dp ∇p) = −αμp (X)p(p − e + N (X)) − pu · ∇μp − Rp (e, p, T ). Φp ∂τp Φe

(26) (27)

Applying the backward difference quotient in the τe -direction to approximate the direction n+1 ∂e derivative ∂e∂τe = ∂τ (X, tn+1 ), e ∂en+1 en+1 (X) − en (X + μe un+1 Δt) (X) ≈ . ∂τe Δts (1 + μ2e |un+1 |2 )1/2 Similarly,

∂pn+1 pn+1 (X) − pn (X − μp un+1 Δt) (X) ≈ . ∂τp Δts (1 + μ2p |un+1 |2 )1/2

Then a product formulation of three differential operators is derived from the electron concentration Equation (26) and differs by a dimensionless  ∂  ∂  ∂  ∂  ∂  ∂  n+1 De 1 − Δts De 1 − Δts De e 1 − Δts ∂x ∂x ∂y ∂y ∂z ∂z 

(28) = en + Δts fe (en+1 , pn+1 , T n+1 ) + O (Δts )2 , where en = e(X + μe un+1 Δts ), fe (en+1 , pn+1 , T n+1 ) = αμe (X)en+1 (pn+1 − en+1 + N (X)) + en+1 un+1 · ∇μe − Re (en+1 , pn+1 , T n+1 ). The values of en+1 are obtained in the following three steps: Step A. to find en+1/3 by  ∂  ∂  n+1/3 1 − Δts De e = en + Δts fe (en+1 , pn+1 , T n+1 ). (29a) ∂x ∂x Step B. to find en+2/3 by  ∂  ∂  n+2/3 1 − Δts De e = en+1/3 . ∂y ∂y Step C. to find en+1 by 

1 − Δts

∂  ∂  n+1 De e = en+2/3 . ∂z ∂z

(29b)

(29c)

It is easy to see that the equations (29a)–(29c) are equivalent to (28) by cancelling the intermediate variables en+1/3 and en+2/3 . In a similar fashion, the hole concentration Equation (27) is divided by  ∂  ∂  ∂  ∂  ∂  ∂  n+1 1 − Δts Dp 1 − Δts Dp 1 − Δts Dp p ∂x ∂x ∂y ∂y ∂z ∂z

 = pn + Δts fp (en+1 , pn+1 , T n+1 ) + O (Δts )2 , (30)

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

1061

where pn = p(X − μp un+1 Δts ), fp (en+1 , pn+1 , T n+1 ) = −αμp (X)pn+1 (pn+1 −en+1 +N )−pn+1 un+1 ·∇μp −Rp (en+1 , pn+1 , T n+1 ), and the values of pn+1 are obtained by three successive computations. The heat conductor Equation (4) is turned into     ∂2 ∂2 ∂2 1 − Δts 2 1 − Δts 2 T n+1 1 − Δts 2 ∂x ∂y ∂z 

n −1 n+1 n+1 n+1 =T + Δts ρ fT (e ,p ,T ) + O (Δts )2 ,

(31)

where fT (en+1 , pn+1 , T n+1 ) = {[De ∇en+1 +μe en+1 un+1 ]−[Dp ∇pn+1 −μp pn+1 un+1 ]}·un+1 , and the divided formulation is equal to the original Equation (4) and different by a dimensionless. Then the values of T n+1 are computed in three steps. Therefore, the method of mixed finite volume element modified with second order characteristic fractional step differences is concluded as follows. Given the approximations {U m , Zm } ∈ Vh × Sh by the method of mixed finite volume element at t = tm and the discrete solun+1 n+1 n+1 tions {Em , Pm , Hm }, we aim to find the difference solutions {Eijk , Pijk , Hijk } at tn+1 = tm + rΔts (r = 1, 2, · · · , j) and to find the approximations {U m+1 , Zm+1 } at the next time level t = tm+1 . Firstly, the values of E are computed by n+1/3

(1 − Δts δx (De δx )) Eijk    n + Δts αμe E n (P n − E n + N ) + E  n EU n+1 · ∇μe − Re (E  n , P n , H n ) =E ijk n+2/3

(1 − Δts δy (De δy )) Eijk

n+1/3

= Eijk

n+2/3

n+1 (1 − Δts δz (De δz )) Eijk = Eijk

, ijk

0 ≤ j ≤ M2 ,

,

(32b)

0 ≤ k ≤ M3 ,

,

0 ≤ i ≤ M1 , (32a)

(32c)

where the second order difference operator in x-direction is defined by  1 De,i+1/2,jk (Ei+1,jk − Eijk ) − De,i−1/2,jk (Eijk − Ei−1,jk ) , h21   1 1 De (Xijk ) + De (Xi+1,jk ) , De (Xijk ) + De (Xi−1,jk ) . = De,i−1/2,jk = 2 2

δx (De δx E)ijk = De,i+1/2,jk

and δy (De δy eh )ijk and δz (De δz eh )ijk are given similarly. EU n+1 is defined similarly to (13). The second order difference operators of P and H in the following formulas are defined similarly. Secondly, the values of P are obtained by n+1/3

(1 − Δts δx (Dp δx )) Pijk  n n + N ) =Pijk + Δts − αμp P n (Pn − E

  n , P n , H n ) , − P n EU n+1 · ∇μp − Rp (E ijk

(1 −

n+2/3 Δts δy (Dp δy )) Pijk

=

n+1/3 Pijk , n+2/3

n+1 (1 − Δts δz (Dp δz )) Pijk = Pijk

,

0 ≤ i ≤ M1 , 0 ≤ j ≤ M2 ,

0 ≤ k ≤ M3 .

Thirdly, the values of H are obtained by n+1/3

(1 − Δts δx (δx )) Hijk ,  n+1 n n =Δts ρ−1 ] − [Dp ∇h P n ijk [De ∇h E + μe E EU

(33a) (33b) (33c)

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1062

 −μp EU n+1 ] ijk · EU n+1 ijk , n+2/3

(1 − Δts δy (δy )) Hijk

0 ≤ i ≤ M1 , n+1/3

= Hijk

n+2/3

n+1 (1 − Δts δz (δz )) Hijk = Hijk

,

,

(34a)

0 ≤ j ≤ M2 ,

(34b)

0 ≤ k ≤ M3 .

(34c)

In the above procedures (32)–(34), the values at outer points are handled by using the reflection principle[5,14,26,27] , which is reasonable because of homogeneous Neumann boundary value condition. The values of E n (X) and P n (X) are assigned by threefold piece-wised n n quadratic interpolations of the values {Eijk } and {Pijk } at twenty-seven points close to X n+1 n n n n n  p,ijk ), X n  = P n (X (see [2]), Eijk = E (Xe,ijk ), Xe,ijk = Xijk + μe,ijk EU ijk Δts , Pijk p,ijk = n+1 Xijk − μp,ijk EU ijk Δts . Initial approximations: 0 = e0 (Xijk ), Eijk

0 Pijk = p0 (Xijk ),

0 Hijk = T0 (Xijk ),

Xijk ∈ Ωh .

(35)

The program runs for a time step as follow. The initial values of electric potential and intensity {U 0 , Z0 } are given by using the method of mixed finite volume element (22) and initial approximations (35). Then by using the characteristic fractional step difference procedures (32), (33) and (34) and using the method of speedup, we can get the values of the other three functions in parallel {E 1 , P 1 , H 1 }, {E 2 , P 2 , H 2 }, · · ·, {E j , P j , H j }. Noting that {E j , P j , H j } = {E1 , P1 , H1 }, and applying the conjugate gradient method, we can get {U 1 , Z1 } from (22). From (32), (33) and (34), we get {E j+1 , P j+1 , H j+1 }, {E j+2 , P j+2 , H j+2 }, · · ·, {E2 , P2 , H2 }. In the above sequence of computations, we can obtain all the numerical solutions. Finally, it is noted that the solutions exist and are sole by the positive definite Condition (C).

4

Convergence Analysis

By the notations and lemmas in Section 2, the procedures (22a) and (22b) are rewritten by the following inner products x y z (Dx Um + Dy Um + Dz Um , v)m = α(Pm − Em + N, v)m ,

−

y z + (Um , wy )y + (Um , wz )z x y (Zm , Dx w + Dy w + Dz wz )m =

∀ v ∈ Sh ,

(36a)

x (Um , wx )x

0,

∀ w ∈ Vh .

(36b)

 ∈ Vh , Z  ∈ Sh , An auxiliary elliptic projection is introduced to determine a pair of function U  x + Dy U  y + Dz U  z , v)m = α(p − e + N, v)m , (Dx U ∀ v ∈ Sh , x x y y z z x y  , w )y + (U  , w )z − (Z,  Dx w + Dy w + Dz wz )m = 0,  , w )x + (U (U

(37a) ∀ w ∈ Vh . (37b)

 ζ = Z  − ψ, ξ = U − U  and β = U  − u are estimated later. Error functions η = Z − Z, Exact solutions of (1)–(8) are supposed to satisfy (9) and (10). By the discussion of Weiser,  , Z)  of (37a) and (37b) exists and is sole, and they are Wheeler[19] , the pair of functions (U bounded as shown in Lemma 4. Lemma 4. If exact solutions and coefficients of (1)–(8) satisfy (9) and (10), then there exist C1 , C2 > 0 independent on h, Δt, t, n such that   ||ζ||m + β  ≤ C1 h2ψ , (38a)      (38b) U  ≤ C2 , ∞

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

1063

where hψ = (h2x + h2y + h2z )1/2 . It is necessary to estimate η and ξ for getting error bound of exact solutions of (21a), (21b) and numerical solutions of (36a), (36b). Subtracting (37a) from (36a), and subtracting (37b) from (36b), we have x y z + Dy ξm + Dz ξm , v)m = α (em − Em − pm + Pm , v)m , (Dx ξm x (ξm , wx )x

+

y (ξm , w y )y

+

z (ξm , wz )z

∀ v ∈ Sh ,

− (ηm , Dx w + Dy w + Dz w )m = 0, x

y

z

(39a) ∀ w ∈ Vh . (39b)

Taking v = ηm and w = ξ m , we have |||ξ m |||2 = α (em − Em − pm + Pm , ηm )m .

(40)

Introducing the duality method to estimate ηm ∈ Sh (see [2,3,10,13]), and defining another elliptic problem as follows, (x, y, z)T ∈ Ω, ∇ · ω = ηm , (x, y, z)T ∈ Ω, ω = ∇ψ,

(41a) (41b)

ω · γ = 0,

(41c)

(x, y, z)T ∈ ∂Ω.

The function ω  ∈ Vh can be found by the following equations  ∂ω  x   ∂ω x  ,v = ,v , ∀ v ∈ Sh , ∂x ∂x  ∂ω  y   ∂ω y  ,v = ,v , ∀ v ∈ Sh , ∂y ∂y  ∂ω  z   ∂ω z  ,v = ,v , ∀v ∈ Sh , ∂z ∂z and its partial derivatives are bounded by  x 2  y 2  z 2  x 2  y 2  z 2       ∂ ω  ∂ω   ∂ω   ∂ω   ∂ ω              +  ∂ ω  ∂y  +  ∂z  ≤  ∂x  +  ∂y  +  ∂z  .  ∂x 

(42a) (42b) (42c)

(43)

From (41), (42) and (43), it holds 2

 ) = (ηm , Dx ω  x + Dy ω  y + Dz ω  z )m ||ηm || = (ηm , ∇ · ω) = (ηm , ∇ · ω     x y z = (ξm ,ω  x )x + (ξm ,ω  y )y + (ξm ,ω  z )z ≤ ξ m  · ||| ω||| . Considering Lemma 3, (44) and the regularity of elliptic problems together, we have  x 2  y 2  z 2  ∂ω   ∂ω   ∂ω  2        ∂x  +  ∂y  +  ∂z  ≤ C ||ηm || .

(44)

(45)

Continue,  x 2  y 2  z 2  ∂ ω  ∂ ω  ∂ ω       ||| ω ||| ≤ ||Dx ω  || + ||Dy ω  || + ||Dz ω  || =  +  +      ∂x ∂y ∂z   x 2  y 2  z 2    ∂ω     +  ∂ω  +  ∂ω  ≤ C ||ηm ||2 . ≤          ∂x ∂y ∂z  2

Therefore,

x 2

y 2

z 2

    ||ηm || ≤ C ξ m  .

(46)

(47)

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1064

It follows from (40) and (47),  2  2      2 2  ξ m  ≤ ε ξ m  + C ||pm − Pm ||m + ||em − Em ||m . Taking ε = 1/2, we have  2    2 2  ξ m  ≤ C ||pm − Pm ||m + ||em − Em ||m ,  2 2 2  ||ηm || ≤ C ||pm − Pm ||m + ||em − Em ||m .

(48)

(49a) (49b)

Let θe = e − E, θp = p − P and χ = T − H denote error functions, where e, p and T are exact solutions, {U, Z} are mixed finite volume element solutions and {E, P, H} are numerical solutions. The partition for computing the concentrations and the temperature is composed of uniform cubes with hs = h = 1/N . Some inner products and norms in discrete spaces l2 (Ω) and h1 (Ω) corresponding to L2 (Ω) and H 1 (Ω) are given as follows[6,29]. v n , wn =

N

n n vijk wijk h3 ,

|v n |0 = v n , v n 1/2 ,

i,j,k=1

[ v n , wn x =

N −1

N

n n vijk wijk h3 ,

[ v n , wn y =

i=0 j,k=1

[ v n , wn z =

N N −1

N N −1 N

n n vijk wijk h3 ,

i=1 j=0 k=1 n n vijk wijk h3 ,

1/2

|δx v n |0 = [ δx v n , δx v n x ,

i,j=1 k=0 1/2

|δy v |0 = [ δy v n , δy v n y ,  2 1/2 |∇h v n |0 = |δs v n |0 . n

1/2

|δz v n |0 = [ δz v n , δz v n z

,

s=x,y,z

Error equations of the electron and the concentrations are argued firstly. The electron concentration Equation (32), cancelled E n+1/3 and E n+2/3 , is turned into the following equivalent expression, n+1 n Eijk −E ijk

Δts

− ∇h (De ∇h E n+1 )ijk

n  n + N )ijk + E  n EU n+1 · ∇μe,ijk =αμe,ijk Eijk (Pn − E ijk ijk  n n n ijk − Re (E , Pijk , Hijk ) − Δts δx (De δx (δy (De δy E n+1 )))ijk

+ δx (De δx (δz (De δz E n+1 )))ijk + δy (De δy (δz (De δz E n+1 )))ijk + (Δts )2 δx (De δx (δy (De δy (δz (De δz E n+1 )))))ijk ,  where ∇h (De ∇h E n+1 )ijk = δs (De δs E n+1 )ijk .



Xijk ∈ Ωh ,

(50)

s=x,y,z

Error equation of electron concentration is derived from (26) (t = tn+1 ) and (50), n+1 n ) θe,ijk − (en (X e,ijk ) − E ijk n

− ∇h (De ∇h θen+1 )ijk Δts 

n+1 

n   n+1 n+1  n+1 n  n  n EU n+1 + eijk uijk − E =αμe,ijk en+1 ijk ijk ijk pijk − eijk + Nijk − Eijk Pijk − Eijk + Nijk  n+1 n+1 n n+1 n n · ∇μe,ijk − Re (en+1 )))ijk ijk , pijk , Tijk ) + Re (Eijk , Pijk , Hijk ) − Δts δx (De δx (δy (De δy θe  n+1 n+1 + δx (De δx (δz (De δz θe )))ijk + δy (De δy (δz (De δz θe )))ijk + (Δts )2 δx (De δx (δy (De δy (δz (De δz θen+1 )))))ijk + εn+1 e,ijk ,

Xijk ∈ Ωh ,

(51)

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

1065

n

n+1 2 n where |εn+1 e,ijk | ≤ M {Δts + h }, X e,ijk = Xijk + μe,ijk uijk Δt, and E (X) denotes the value of the electron concentration E by threefold piecewise quadratic interpolation of the values n {Eijk } at nearby points. The numerator of the first term on the left-hand side is changed into a summation of three brackets decomposed into n n+1  n ) = (ξ n+1 − ξn ) + (en (X n ) − en (X  n )) − (I − I2 )en (X  n ), (52) − (en (X e,ijk ) − E θe,ijk e,ijk ijk ijk e,ijk e,ijk ijk

where I denotes an identity operator and I2 denotes the quadratic interpolation operator. Then error Equation (51) is rewritten by n+1 n − θe,ijk θe,ijk

n  n ) (I − I2 )en (X n ) en (X e,ijk ) − en (X e,ijk e,ijk − ∇h (De ∇h θen+1 )ijk + − Δts Δts Δts 

n+1 

n   n+1 n+1  n+1 n n  n ijk + eijk uijk − E =αμe,ijk en+1 EU n+1 ijk ijk pijk − eijk + Nijk − Eijk Pijk − Eijk + Nijk  n+1 n+1 n n+1 n n )))ijk · ∇μe,ijk − Re (en+1 ijk , pijk , Tijk ) + Re (Eijk , Pijk , Hijk ) − Δts δx (De δx (δy (De δy θe  + δx (De δx (δz (De δz θen+1 )))ijk + δy (De δy (δz (De δz θen+1 )))ijk

+ (Δts )2 δx (De δx (δy (De δy (δz (De δz θen+1 )))))ijk + εn+1 e,ijk ,

Xijk ∈ Ωh .

(53)

n+1 Testing both sides of (53) by θe,ijk , making summation, and using summation by parts, we have

 θn+1 − θn e

e

 , θen+1 + De ∇h θen+1 , ∇h θen+1

Δts   (I − I )en (X   en (X n ) − en (X en ) en ) 2 e + , θen+1 − , θen+1 Δts Δts    n+1 n+1

n+1 n n  n + N , θen+1 =αμe e −e +N −E P −E p    n EU n+1 · ∇μe , ∇θen+1 − R(en+1 , pn+1 , T n+1 ) − R(E  n , Pn , H n ), θen+1 +  en+1 un+1 − E   − Δts δx (De δx (δy (De δy θen+1 ))), θen+1 + · · · + δy (De δy (δz (De δz θen+1 ))), θen+1   + (Δts )2 δx (De δx (δy (De δy (δz (De δz θen+1 ))))), θen+1 + εn+1 , θen+1 . (54) e Partition parameters are supposed to satisfy the constriction, 3/2

Δts = O(h2 ),

h2 = o(hψ ).

Introducing the following induction hypothesis,     sup sup |θen |∞ , θpn ∞ , |χn |∞ → 0, 0≤n≤L

0≤n≤L (m=[n/j])

    ξ m(n) 

(55)

∞

→ 0,

(h, hψ , Δts ) −→ 0.

(56) Firstly error Equation (54) is estimated, and the first term on the left-hand side is decomposed into  θn+1 − θn   θn+1 − θn   θn − θn  e e e e (57) , θen+1 = e , θen+1 + e , θen+1 . Δts Δts Δts The first term on the right-hand side in the above expression is estimated as follows,  θn+1 − θn e

e

Δts

 , θen+1 ≥

By using n θe,ijk

n − θe,ijk =−

 1  n+1 n+1 θe , θe − θen , θen . 2Δts

n  X ijk

Xijk

  ∇θen · EU ijk / EU ijk  dσ,

(58a)

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1066

the induction hypothesis (56), the constriction (55), (49) and negative norm estimates, we can estimate the second term on the right-hand side of (57) by  θn − θn e

e

Δts

 2   2   , θen+1 ≤ ε ∇h θen+1 0 + M θen  .

(58b)

 1  n+1 2 |θe |0 − |θen |20 − M |θen |20 − ε|∇h θen+1 |20 . 2Δts

(59a)

0

It follows from (57) and (58),  θn+1 − θn e

e

Δts

 , θen+1 ≥

For the third term on the left-hand side of (54), using the following result[3,26,27] derived from (38) and (49),   (59b) |||u − U|||2 ≤ M |θe |20 + |θp |20 + h4ψ , we have  en (X n ) − en (X  n) e

e

Δts

   , θen+1 ≤ M (Δts )2 + h4 + h4ψ + |θen+1 |20 + |θe,m |20 + |θe,m−1 |20 .

(59c)

For the fourth term, we have −

 (I − I )en (X   n)   2 e , θen+1 ≤ M (Δts )2 + h4 + |θen+1 |20 . Δts

(59d)

Giving analysis for all the terms on the right-hand side of (54), we have   



    n + N , θn+1  αμe en+1 pn+1 − en+1 + N − E n P n − E e    2 2 2 2 2 2 ≤M (Δts )2 + h4 + h4ψ + θen+1 0 + |θen |0 + |θe,m |0 + |θe,m−1 |0 + |θp,m |0 + |θp,m−1 |0 , (60a)     n+1 n+1  n n+1 n+1  · ∇μe , ∇θe  u − E EU  e     n , P n , H n ), θen+1  + Re (en+1 , pn+1 , T n+1 ) − Re (E  2  2 2 ≤M (Δts )2 + h4 + h4ψ + θen+1 0 + |θen |0 + |θe,m |0  2 2 2 2 (60b) + |θe,m−1 |0 + |θp,m |0 + |θp,m−1 |0 + |χn |0 ,    n+1 n+1    εe , θe  ≤ M (Δts )2 + h4 + θen+1 2 . (60c) 0 Multiplying both sides of (54) by 2Δts , summing on n from 0 to L, and using ξ 0 = ζ 0 = 0, (59) and (60), we have L

|θeL+1 |20 +  ≤M

|∇h θen+1 |20 Δts

n=0 L 

 n+1 2  n+1 2  n+1 2  θ  + θ  + χ  Δts + h4 + h4 + (Δts )2 e p ψ 0 0 0



n=0

− 2Δts

L   δx (De δx (δy (De δy θen+1 ))), θen+1 + · · · + δy (De δy (δz (De δz θen+1 ))), θen+1 Δts n=0

+ 2(Δts )2

L

δx (De δx (δy (De δy (δz (De δz θen+1 ))))), θen+1 Δts .

n=0

(61)

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

1067

Error analysis is going for the terms on the right-hand side of (61). Though the operators −δx (De δx ), −δy (De δy ), · · · are self-conjugate, positive definite and bounded, their products are not generally commutative. Note that δx δy = δy δx , δx δy = δy δx , · · ·, we have − (Δts )2 δx (De δx (δy (De δy θen+1 ))), θen+1

=(Δts )2 De δx (δy (De δy θen+1 )), δx θen+1 =(Δts )2 δy δx (De δy θen+1 ), De δx θen+1

= − (Δts )2 δx (De δy θen+1 ), δy (De δx θen+1 ) = − (Δts )2 De δx δy θen+1 + δx De · δy θen+1 , De δy δx θen+1 + δy De · δx θen+1 

2 = − (Δts )2 De,i,j+1/2,k De,i+1/2,jk δx δy θen+1 ijk h1 h2 h3 Xijk ∈Ωh

 n+1 n+1 n+1 n+1 + De,i,j+1/2,k δy De,i+1/2,jk δx δy θe,ijk · δx θe,ijk + De,i+1/2,jk δx De,i,j+1/2,k δx δy θe,ijk · δy θe,ijk  n+1 n+1 · h1 h2 h3 + δx De,i,j+1/2,k · δy De,i+1/2,jk δx θe,ijk · δy θe,ijk h1 h2 h3

  n+1 2 n+1 De,i,j+1/2,k δy De,i+1/2,jk δx δy θe,ijk δx δy θe,ijk h1 h2 h3 − (Δts )2 ≤ − (Δts )2 D∗2 

Xijk ∈Ωh

·

n+1 δx θe,ijk

+

Xijk ∈Ωh

 n+1 · δy θe,ijk h1 h2 h3  n+1 · δy θe,ijk h1 h2 h3 .

n+1 De,i+1/2,jk δx De,i,j+1/2,k δx δy θe,ijk

n+1 + δx De,i,j+1/2,k · δy De,i+1/2,jk δx θe,ijk

Using the positive definite condition (C) and Cauchy inequality, cancelling high-order difference  2 n+1 operators δx δy θe,ijk , we have −(Δts )2 δx (De δx (δy (De δy θen+1 ))), θen+1 ≤ M Δts ∇h θen+1 0 Δts . Then we have the following estimate for the second term on the right-hand side of (61), − 2Δts

L   δx (De δx (δy (De δy θen+1 ))), θen+1 + · · · + δy (De δy (δz (De δz θen+1 ))), θen+1 Δts n=0

L   ∇h θen+1 2 Δts . ≤M Δts 0

(62)

n=0

The last term of (61) is argued as follows by the collection of (C), (55), inverse estimates and Cauchy inequality, (Δts )3 δx (De δx (δy (De δy (δz (De δz θen+1 ))))), θen+1

= − (Δts )3 δx (δy (De δy (δz (De δz θen+1 )))), De δx θen+1 = − (Δts )3 δy δx (De δy (δz (De δz θen+1 ))), De δx θen+1 =(Δts )3 δx (De δy (δz (De δz θen+1 ))), δy (De δx θen+1 )

=(Δts )3 De δx δy (δz (De δz θen+1 )) + δx De · δy δz (De δz θen+1 ), De δy δx θen+1 + δy De · δx θen+1  =(Δts )3 De δz δx δy (De δz θen+1 ), De δx δy θen+1 + δy De · δx θen+1  + δx De · δz δy (De δz θen+1 ), De δx δy θen+1 + δy De · δx θen+1  

 = − (Δts )3 δx δy (De δz θen+1 ), δz De De δx δy θen+1 + δy De · δx θen+1 

  + δy (De δz θen+1 ), δz δx De De δx δy θen+1 + δy De · δx θen+1  = − (Δts )3 De δx δy δz θen+1 + δx δy De · δz θen+1 , De De δx δy δz θen+1 + δz (De De ) · δx δy θen+1 + δz (De δy De ) · δx θen+1 + De δy De · δx δz θen+1 + De δy δz θen+1 + δy De · δz θen+1 , δx De

 · De δx δy δz θen+1 + δz (δy De · De ) · δx δy θen+1 + δy De δx δz θen+1 + δz (δy De · δy De ) · δx θen+1

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1068



= − (Δts )3

Xijk ∈Ωh





 n+1 2 De,i+1/2,jk De,i,j+1/2,k De,ij,k+1/2 δx δy δz θe,ijk h1 h2 h3

n+1 + De,ij,k+1/2 · δz (De,i+1/2,jk De,i,j+1/2,k )δx δy θe,ijk n+1 + De,ij,k+1/2 De,i,j+1/2,k δy De,i+1/2,jk · δx δy θe,ijk

  n+1 n+1 + De,i+1/2,jk De,ij,k+1/2 De,i,j+1/2,k δy δz θe,ijk + · · · δx δy δz θe,ijk h1 h2 h3 + · · ·   ≤M Δts |∇h θen+1 |20 + |θen+1 |20 Δts .

(63)

Applying (62) and (63), we estimate (61) for Δts suitably small L    L+1 2 ∇h θen+1 2 Δts  θe  + 0

 ≤M

0

n=0 L   n+1 2  n+1 2  n+1 2   Δts + h4 + h4ψ + (Δts )2  θe  +  θp  +  χ 0

0



0

.

(64)

n=0

Cancelling P n+1/3 and P n+2/3 from (3), we have its equivalent formulation: n+1 n Pijk − Pijk

Δts

− ∇h (Dp ∇h P n+1 )ijk

n  n + N )ijk − P n EU n+1 · ∇μp,ijk − Rp (E n , P n , H n ) = − αμp,ijk Pijk (Pn − E ijk ijk ijk ijk ijk

− Δts {δx (Dp δx (δy (Dp δy P n+1 )))ijk + δx (Dp δx (δz (Dp δz P n+1 )))ijk + δy (Dp δy (δz (Dp δz P n+1 )))ijk } + (Δts )2 δx (Dp δx (δy (Dp δy (δz (Dp δz P n+1 )))))ijk , Xijk ∈ Ωh . n+1

Error equation of the hole concentration follows from (27) (t = t n+1 θp,ijk −

n (pn (X p,ijk )

(65)

) and (65),

n − Pijk )

− ∇h (Dp ∇h θn+1 )p,ijk 

n+1 

n  n+1 n  n = − αμp,ijk pn+1 ijk pijk − eijk + Nijk − Pijk Pijk − Eijk + Nijk   n+1 n+1 n+1 n+1 n n n n · ∇μp,ijk − R(en+1 − pn+1 ijk uijk − Pijk EU ijk ijk , pijk , Tijk ) + R(Eijk , Pijk , Hijk )   − Δts δx (Dp δx (δy (Dp δy θpn+1 )))ijk · · · + δy (Dp δy (δz (Dp δz θpn+1 )))ijk Δts

+ (Δts )2 δx (Dp δx (δy (Dp δy (δz (Dp δz θpn+1 )))))ijk + εn+1 p,ijk ,

(66)

n X p,ijk

2 = Xijk − μp,ijk un+1 where |εn+1 p,ijk | ≤ M {Δts + h }, ijk Δts . After a numerical analysis similar to the electron concentration, we have the estimate of the hole concentration L L    n+1 2  |θe |0 +|θpn+1 |20 +|χn+1 |20 Δts +h4 +h4ψ +(Δts )2 . (67) |∇h θpn+1 |20 Δts ≤ M |θpL+1 |20 + n=0

n=0

The heat conductor Equation (34) is equivalent to the following formulation by cancelling H n+1/3 and H n+2/3 , n+1 n Hijk − Hijk

− ∇h (∇h H n+1 )ijk Δts   n+1 n+1 −1 n n n n = [De,ijk ∇h Eijk + μe,ijk Eijk EU n+1 ijk ] − [Dp,ijk ∇h Pijk − μp,ijk Pijk EU ijk ] · EU ijk ρijk   − Δts δx (δx (δy (δy H n+1 )))ijk + · · · + δy (δy (δz (δz H n+1 )))ijk + (Δts )2 δx (δx (δy (δy (δz (δz H n+1 )))))ijk , Xijk ∈ Ωh .

(68)

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

1069

Error equation of the heat conduction is derived from (4) (t = tn+1 ) and (68), n χn+1 ijk − χijk

− ∇h (∇h χn+1 )ijk Δts  n+1 n+1 n+1 n+1 n n n = [De,ijk ∇h θe,ijk + μe,ijk (en+1 ijk uijk − Eijk EU ijk )] − [Dp,ijk ∇h θp,ijk − μp,ijk (pijk uijk   n+1 −1 n+1 n+1 n+1 n − Pijk EU n+1 ijk )] · uijk ρijk + [De,ijk ∇h eijk + μe,ijk eijk uijk ]  n+1 n+1 − [Dp,ijk ∇h pn+1 ijk − μp,ijk pijk uijk ]   n+1 −1 n+1 · (un+1 )))ijk · · · + δy (δy (δz (δz χn+1 )))ijk ijk − EU ijk )ρijk − Δts δx (δx (δy (δy χ + (Δts )2 δx (δx (δy (δy (δz (δz χn+1 )))))ijk + εn+1 T,ijk ,

(69)

2 where |εn+1 T,ijk | ≤ M {Δts + h }. Testing both sides of (69) by χn+1 ijk , making summation and using summation by parts, we have   χn+1 − χn , χn+1 + ∇h χn+1 , ∇h χn+1 Δts =[De ∇h θen + μe (en+1 un+1 − E n EU n+1 )] · un+1 ρ−1 , χn+1

− [Dp ∇h θpn − μp (pn+1 un+1 − P n EU n+1 )] · un+1 ρ−1 , χn+1

+ {[De ∇h en+1 + μe en+1 un+1 ] − [Dp ∇h pn+1 − μp pn+1 un+1 ]} · (un+1 − EU n+1 )ρ−1 , χn+1   − Δts δx (δx (δy (δy χn+1 ))), χn+1 + · · · + δx2 (δy (δz (δz χn+1 ))), χn+1 , χn+1 . + (Δts )2 δx (δx (δy (δy (δz (δz χn+1 ))))), χn+1 + εn+1 T

(70)

Using (55), (56), (38) and (49), we can estimate the terms of (70),  χn+1 − χn

 , χn+1 ≥

 1  n+1 2 |0 − |χn |20 , |χ Δts 2Δts n+1 n+1 ∇h χ , ∇h χ = |∇h χn+1 |20 ,

(71a) (71b)

[De ∇h θen + μe (en+1 un+1 − E n EU n+1 )] · un+1 ρ−1 , χn+1 − [Dp ∇h θpn − μp (pn+1 un+1 − P n EU n+1 )] · un+1 ρ−1 , χn+1

+ {[De ∇h en+1 + μe en+1 un+1 ] − [Dp ∇h pn+1 − μp pn+1 un+1 ]} · (un+1 − EU n+1 )ρ−1 , χn+1      2 2  2 ≤ε |∇h θen |20 + ∇h θpn 0 + M (Δts )2 + h4 + h4ψ + χn+1 0 + |θen |20 + θpn 0  2 2 2 2 (72a) + |θe,m |0 + |θe,m−1 |0 + |θp,m |0 + |θp,m−1 |0 ,    n+1 n+1   2 ε ,χ  ≤ M (Δts )2 + h4 + χn+1 0 . (72b) T Multiplying both sides of (70) by 2Δts , summing on n from 0 to L, and collecting χ0 = 0, (71) and (72) together, we have |χL+1 |20 +

L

|∇h χn+1 |20 Δts

n=0 L   |∇h θen |20 + |∇h θpn |20 Δts ≤ε n=0

+M

L    n+1 2  |0 + |θen+1 |20 + |θpn+1 |20 Δts + h4 + h4ψ + (Δts )2 |χ n=0

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

1070

L   δx (δx (δy (δy χn+1 ))), χn+1 + · · · + δy (δy (δz (δz χn+1 ))), χn+1 − 2(Δts ) 2

n=0

+ 2(Δts )3

L

δx (δx (δy (δy (δz (δz χn+1 ))))), χn+1 .

(73)

n=0

Giving a discussion similar to the procession of (64), we get |χL+1 |20 +

L

|∇h χn+1 |20 Δts

n=0

≤ε

L

  |∇h θen |20 + |∇h θpn |20 Δts

n=0

+M

L    n+1 2  |χ |0 + |θen+1 |20 + |θpn+1 |20 Δts + h4 + h4ψ + (Δts )2 .

(74)

n=0

Considering the resulting estimates (64), (67) and (74) together, we obtain |θeL+1 |20 + |θpL+1 |20 + |χL+1 |20 +

L   |∇h θen+1 |20 + |∇h θpn+1 |20 + |∇h χn+1 |20 Δts n=0

≤M

L  

  |θen+1 |20 + |θpn+1 |20 + |χn+1 |20 Δts + h4 + h4ψ + (Δts )2 .

(75)

n=0

Applying Gronwall inequality, we have |θeL+1 |20 + |θpL+1 |20 + |χL+1 |20 +

L   |∇h θen+1 |20 + |∇h θpn+1 |20 + |∇h χn+1 |20 Δts n=0

4

≤M {h +

h4ψ

2

+ (Δts ) }.

(76)

It remains to testify the induction hypothesis (56). From the that initial errors are zero and the estimate (59c), we see that (56) is true obviously for n = 0. Suppose that the statement is true for 1 ≤ n ≤ L. By (76) and (55), we have      L+1     + θL+1  + χL+1  ≤ M h−3/2 h2 + h2 + Δts ≤ M h1/2 −→ 0, θ e p ψ ∞ ∞ ∞ (77a)      −3/2  2 1/2 h + h2ψ + Δts ≤ M hψ −→ 0, ≤ M hψ (77b) σ [(L+1)/j]  0,∞

then we get that (56) is correct for n = L + 1. From (38), (49) and (76), we conclude the following theorem. Theorem 1. Exact solutions of (1)–(8) are supposed to be regular (R) and the problem is assumed to be positive definite (C). The method of mixed finite volume element modified with second order characteristic fractional step difference (22), (32)–(34) is adopted and the constriction (55) holds, then ||e − E||L∞ ((0,T],l2 ) + ||p − P ||L∞ ((0,T],l2 ) + ||T − H||L∞ ((0,T],l2 ) + ||e − E||L ((0,T],h1 ) 2 + ||p − P ||L ((0,T],h1 ) + ||T − H||L ((0,T],h1 ) + ||ψ − Z||L∞ ((0,T],L2 ) + ||u − U||L∞ ((0,T],V ) 2 2   ∗ 2 2 ≤ M Δts + h + hψ , (78)

Mixed-finite Volume Characteristic Fractional Difference of Semiconductor Device

where ||g||L∞ (J;X) = sup ||g n ||X , nΔt≤T

||g||L2 (J;X) = sup

L 

LΔt≤T

1071

 ||g n ||X Δt ,

n=0

∗

the constant M is dependent on ψ, e, p, T and their derivatives, and Δt =

5

Δts ,

for the concentration

Δtψ ,

for the potential.

Conclusions and Discussions

The main work of this paper is to construct the procedure to approximate the behavior of semiconductor device by combining mixed finite volume element method with characteristic fractional step difference, and to provide numerical analysis on it. In the first section, some mathematical model, research background and international research status are introduced. Preliminary notations, two different partitions of the spatial domain and several lemmas are given in the second section. In the third section, the composite procedure to solve semiconductor device problem is constructed, where the mixed finite volume element is used to approximate the electric intensity, and the characteristic fractional step differences are applied to compute the concentrations and the temperature. By the proposed scheme, the three-dimensional problem is solved by solving three successive one-dimensional problems and computational cost is decreased. In the fourth section error estimate of optimal order in l2 norm for the numerical solution is obtained. In this paper there are serval interesting conclusions stated as follows. (I) The effect of the heat is considered in constructing the mixed finite volume element method of the transient behavior problems. Thus, the numerical approximation can reflect physical natures quite well. (II) This method can be applied in large-scaled numerical computation of three-dimensional problem with complicated geometric domain. (III) A mixed finite volume element method can keep the law of mass conservation, by which we can improve computation accuracy of one order in computing electric field intensity. It is very important to understand numerical simulation of semiconductor device. (IV) The fractional step algorithm modified by characteristic difference is adopted to compute the concentrations and the heat which works on modern parallel computers and gives numerical computation of second-order accuracy and high efficiency in parallel for transient behavior problems. Acknowledgements. The authors express their deep appreciation to Prof. J. Douglas Jr, and Prof. Jiang Lishang for their helpful suggestions in the serial of research on numerical simulation of semiconductor device. References [1] Bank, R.E., Fichtner, W.M., Rose, D.J., et. al. Transient simulation of sillcon devices and circuits. IEEE Computer-Aided Design, 6: 436–451 (1985) [2] Ciariet, P.G. The finite element methods for elliptic problem. North-Holland, Amsterdam, 1978 [3] Douglas, Jr. J. Simulation of miscible displacement in porous media by a modified method of characteristics procedure. In: Numerical Analysis, Dundee, 1981, Lecture Note in Mathematics 912, Springer-Verlag, Berlin, 1982 [4] Douglas, Jr. J., Ewing, R.E., Wheeler, M.F. Approximation of the pressure by a mixed method in the simulation of miscible displacement. RAIRO Anal. Numer., 17(1): 17–33 (1983)

1072

Y.R. YUAN, Q. YANG, C.F. LI, T.J. SUN

[5] Douglas, Jr.J., Yuan, Y.R. Numerical simulation of immiscible flow in porous media based on combining the method of characteristics with mixed finite element procedure. Numerical Simulation in Oil Recovery, 119–132, New York: Springer-Verlag, 1986 [6] Douglas, Jr.J., Yuan, Y.R. Finite difference methods for transient behavior of a semiconductor device. Mat. Apli. Comp., 6(1): 25–38 (1987) [7] Gummel, H.K. A self-consistent iterative scheme for one-dimensional steady-state transistor calculation. IEEE Trans: Electron Device, 11: 455–465 (1964) [8] He Y., Wei T.F. Computer simulation method for semiconductor device. Scicence Press, Beijing, 1989 [9] Jerome, J.W. Mathematical Theory and Approximation of Semiconductor Models. Philadelphia, SIAM, 1994 [10] Jiang L.S., Pang Z.Y. Finite element method and its theory. People’s Education Press, Beijing, 1979 [11] Jones, J.E. A mixed finite volume method for accurate computation of fluid velocities in porous media. Ph. D. Thesis, University of Colorado, Donrer. Co., 1995 [12] Lou, Y. On basic semiconductor equation with heat conduction. J. Partial. Diff. Eqs., 1: 43–54 (1995) [13] Nitsche, J. Linear splint-funktionen and die von Rits for elliptishce randwort probleme. Arch. for Rational Mech. and Anal., 36: 348–355 (1968) [14] Russell, T.F. Time stepping along characteristics with incomplete iteration for a Galerkin approximation of miscible displacement in porous media. SIAM J. Numer. Anal., 22(5): 970–1013 (1985) [15] Russell, T.F. Rigorous black-centered discritization on irregular grids: Improved simulation of complex reservoir systems. Project Report, Research Corporation, Tulsa, 1995 [16] Shi, M. Physics of modern semiconductor device. Science Press, Beijing, 2002 [17] Sun, C.W., Lu, Q.S., Fan, Z.X. Laser irradiation effect. National Defence Industry Press, Beijing, 2002 [18] Sun, T.J., Yuan, Y.R. An approximation of incompressible miscible displacement in porous media by mixed finite element method and characteristics-mixed finite element method. J. Comp. Appl. Math., 228: 391–411 (2009) [19] Weiser, A., Wheeler, M.F. On convergence of block-centered finite difference for elliptic problems. SIAM J. Numer. Anal., 25(2): 351–375 (1988) [20] Yang, Q, Yuan, Y.R. An approximation of semiconductor device by mixed finite element method and characteristics-mixed finite element method. Applied Mathematics and Computation, 225: 407–424 (2013) [21] Yang, Q., Yuan, Y.R. An approximation of semiconductor device of heat conduction by mixed finite element method characteristics-mixed finite element method. Applied Numerical Mathematics, 70: 42–57 (2013) [22] Yang, Q., Yuan, Y.R. An approximation of three-dimensional semiconductor device by mixed finite element method and characteristics-mixed finite element method. Numerical Mathematics: Theory, Methods and Application, 8(3): 356–382 (2015) [23] Yuan, Y.R., Ding, L.Y., Yang, H. A new method and theoretical analysis of numerical analog of semiconductor. Chinese Science Bulletin, 27(7): 790–795 (1982) [24] Yuan, Y.R. Characteristics method with mixed finite element for transient behavior of semiconductor device. Chin. Sci. Bull., 36(17): 1356–1357 (1991) [25] Yuan Y.R. The approximation of the electronic potential by a mixed method in the simulation of semiconductor. J. Systems Sci. Math. Sci., 11(2): 117–120 (1991) [26] Yuan, Y.R. Time stepping along characteristics for the finite element approximation of compressible miscible displacement in porous media. Math. Numer. Sinica, 14(4): 385–400 (1992) [27] Yuan, Y.R. Finite difference methods for a compressible miscible displacement problem in porous media. Math. Numer. Sinica, 15(1): 16–28 (1993) [28] Yuan, Y.R. Finite element method and analysis of numerical simulation of semiconductor device. Acta Math. Sci., 13(3): 241–251 (1993) [29] Yuan, Y.R. Finite difference method and analysis for three-dimensional semiconductor device of heat conduction. Sci. China Math., 39(11): 1140–1151 (1996) [30] Yuan, Y.R. Characteristic finite difference fractional step methods for three-dimensional semiconductor device of heat conduction. Chin. Sci. Bull., 45(2): 123–131 (2000) [31] Yuan, Y.R. Finite difference fractional step method for transient behavior of a semiconductor device. Acta Mathematica Scientia (Series B), 3: 427–438 (2005) [32] Yuan, Y.R. Modification of upwind finite difference fractional step methods by the transient state of the semiconductor device. Numer. Methods Partial Differential, 24: 400–417 (2008) [33] Yuan, Y.R. Recent progress in numerical methods for semiconductor devices. Chinese J. Computational Physics, 26(3): 317–324 (2009) [34] Yuan, Y.R. Theory and application of reservoir numerical simulation, Chapter 7, numercial method of transient behavior of semiconductor device. Science Press, Beijing, 2013