Exponential decay of a finite volume scheme to the thermal equilibrium for drift--diffusion systems

0.45 MB·

by  Marianne Bessemoulin-Chatard

Exponential decay of a finite volume scheme to the thermal
equilibrium for drift–diffusion systems
M. Bessemoulin-Chatard
Univ. Nantes, CNRS, UMR 6629 - Laboratoire Jean Leray, F-44000 Nantes
C. Chainais-Hillairet
6
Univ. Lille, CNRS, UMR 8524 - Laboratoire Paul Painlev´e, F-59000 Lille
1
0
2 January 6, 2016

n
a
J Abstract

5
Inthispaper,westudythelarge–timebehaviorofanumericalschemediscretizingdrift–

diffusion systems for semiconductors. The numerical method is finite volume in space,
]
A implicit in time, and the numerical fluxes are a generalization of the classical Scharfetter–
N Gummel scheme which allows to consider both linear or nonlinear pressure laws.
We study the convergence of approximate solutions towards an approximation of the ther-
.
h mal equilibrium state as time tends to infinity, and obtain a decay rate by controlling the
t discrete relative entropy with theentropy production. This result is proved underassump-
a
m tionsofexistenceanduniform-in-timeL∞ estimatesfornumericalsolutions,whicharethen
discussed. Weconclude bypresenting some numerical illustrations of the stated results.
[

1
1 Introduction
v
3
1 1.1 Aim of the article
8
0 The Van Roosbroeck’s drift–diffusion system is a fundamental model for the mathematical
0 description and numerical simulation of semiconductor devices. It consists of two parabolic
.
1 convection–diffusion–reaction equations for the carrier densities (electrons and holes), and a
0 Poisson’s equation for the electrostatic potential. Global existence and uniqueness results have
6
been obtained for this model under natural assumptions [12, 14, 35]. Moreover, this system
1
is shown to be dissipative. Indeed, it admits a Lyapunov functional, which may be physically
:
v interpreted as an energy. Furthermore, it has been proved using this energy functional that the
i
X solution of the Van Roosbroeck system converges at an exponential rate to the thermal equilib-
rium state if the boundary conditions are in thermal equilibrium [13, 15].
r
a
Theclassicaldrift–diffusionmodelisbasedonBoltzmannstatistics. Moreprecisely,itmeans
that the statistical distribution function describing the dependence of the carrier densities on
the chemical potentials is the exponential function. However, this choice may fail to describe
relevantly the physical reality in some cases (for example in case of high carrier densities).Then
other statistics have to be considered, like Fermi-Dirac statistics for instance [15]. This leads
to a modification of the diffusive terms, which become nonlinear. Existence and uniqueness
of weak solutions to a nonlinear drift-diffusion model have been proved in [24]. Dissipativity
and convergence to the thermal equilibrium for large time have also been established for this
1
generalized model [26]. Let us underline that more recently, driven by applications like organic
semiconductors, there is an increasedinterest in drift–diffusion models with arbitrary statistical
distribution functions [11, 40].
From a numerical point of view it is essential to consider numerical schemes which preserve
the main qualitative properties of the continuous system, such as positivity of the densities,
dissipativity and consistency with thermal equilibrium. In the case of Boltzmann statistics, the
Scharfetter–Gummel scheme [23, 38] is widely used. It exploits the exponential dependence on
the chemical potential and allows to recover the correctlarge time behavior. Various extensions
of the Scharfetter–Gummel scheme have been suggested to account the diffusion enhancement
inducedbynonBoltzmannstatistics[25,37,39]. Unfortunately,theyarenotthermodynamically
consistent. More recently, a consistent generalization in the spirit of the original Scharfetter–
Gummel scheme was proposed [9] and applied to simulate organic semiconductors’ behavior
[28, 29]. This method leads to solve a nonlinear boundary value problem at each interface.
Inthispaper,wepreferentiallyfocusonanotherextensionoftheScharfetter–Gummelscheme
usingaproperaverageofthenonlineardiffusion[3]whichguaranteesthermodynamicconsistency.
An alternative interpretation of this scheme is given in [30] and applied to a very general class
of statistical distribution functions arising in organic semiconductors modeling.
Here our aim is to study the large–time behavior of an implicit in time and finite volume in
space discretization of the drift–diffusion system, with a Scharfetter–Gummel approximation of
the convection–diffusion fluxes. Our proof is based on the entropy–dissipationmethod [1]. This
point of view was already adopted in several articles [7, 13, 18, 19, 20]. The crucial point to
obtain an exponential decay rate of the approximate solutions towards the equilibrium is the
control of the relative entropy by the entropy production. The choice of Scharfetter–Gummel
type fluxes for the discretization of the convection–diffusionfluxes is essential at this step.
1.2 The drift–diffusion system and the thermal equilibrium
LetΩbeanopenboundedsubsetofRd(d 1)correspondingtothegeometryofasemiconductor
≥
device and T > 0. This device can be described by the so-called drift–diffusion system. This
system consists of two continuity equations for the electron density N and the hole density P,
and a Poissonequation for the electrostatic potential Ψ. It writes for all (x,t) Ω [0,T]:
∈ ×
∂ N +div(µ ( r(N)+N Ψ))= R(N,P), (1)
t N
−∇ ∇ −
∂ P +div(µ ( r(P) P Ψ))= R(N,P), (2)
t P
−∇ − ∇ −
λ2∆Ψ=P N +C. (3)
− −
The given function C(x) is the doping profile describing fixed background charges. The dimen-
sionlessphysicalparametersµ , µ andλarethe rescaledmobilitiesofelectronsandholes,and
N P
the rescaled Debye length respectively. The definition of r depends on the statistics chosen to
describetherelationbetweenthedensitiesandthechemicalpotentials. Theusualconsiderations
on which the isentropic hydrodynamic model are based suggest a pressure of the form:
r(s)=sα, α 1.
≥
The linear case, where α = 1, is the isothermal model, corresponding to Boltzmann statistics.
The system (1)–(3) is supplemented with initial conditions:
N(x,0)=N (x), P(x,0)=P (x), x Ω, (4)
0 0
∈
2
and with mixed boundary conditions: Dirichlet boundary conditions on the ohmic contacts
and homogeneous Neumann boundary conditions on the insulated boundary segments. More
precisely, the boundary ∂Ω is split into ∂Ω = ΓD ΓN with ΓD ΓN = , and the boundary
∪ ∩ ∅
conditions write:
N(γ,t)=ND(γ), P(γ,t)=PD(γ), Ψ(γ,t)=ΨD(γ), (γ,t) ΓD [0,T], (5)
∈ ×
( r(N) ν)(γ,t)=( r(P) ν)(γ,t)=( Ψ ν)(γ,t)=0, (γ,t) ΓN [0,T], (6)
∇ · ∇ · ∇ · ∈ ×
whereν isthe unit normalto∂Ω outwardto Ω. Inthis paper,weassumethat the mobilities are
constant and equal: µ =µ =1, and we need the following general assumptions:
N P
Hypotheses1. Thedomain Ωisanopenboundedpolygonal (orpolyhedral) subsetofRd (d 1)
≥
and ∂Ω=ΓD ΓN, with ΓD ΓN = and m(ΓD)>0. The doping profile C belongs to L∞(Ω).
∪ ∩ ∅
The boundary conditions ND, PD and ΨD are traces of some functions defined on the whole
domain Ω, still denoted by ND, PD and ΨD. Furthermore, we assume that
N , P L∞(Ω), (7)
0 0
∈
ND, PD L∞ H1(Ω), ΨD H1(Ω), (8)
∈ ∩ ∈
M >0, m 0 such that m N ,P ,ND,PD Ma.e. on Ω. (9)
0 0
∃ ≥ ≤ ≤
For the drift–diffusion model with linear pressure r=Id, the recombination–generationrate
can usually be written under the following form [32]:
R(N,P)=R (N,P)(NP 1). (10)
0
−
This general form includes in particular the Shockley–Read–Hallterm:
NP 1
R (N,P)= − , τ , τ , τ >0,
SRH P N C
τ N +τ P +τ
P N C
or the Auger recombination:
R =(C N +C P)(NP 1).
AU N P
−
Inthecaseofanonlinearpressure,theShockley–Read–Halltermcannotbetakenanymore. Some
recent works about organic electronic devices (see [22] for instance) include several propositions
for the modeling of generationandrecombinationprocessesin the nonlinear case. However,this
leads to highly nonlinear and intricate source terms, and we choose to consider either the linear
case with R defined by (10), or the nonlinear case with R = 0. More precisely, we will need to
assume that either
Hypotheses 2.
r =Id, (11)
R(N,P)=R (N,P)(NP 1), with R continuous and nonnegative, (12)
0 0
−
NDPD =1, (13)
or
Hypotheses 3.
r 1(R), r(0)=r′(0)=0, r′(s) c sα−1, α>1, (14)
0
∈C ≥
R=0. (15)
3
Existenceanduniquenessofweaksolutionstothedrift–diffusionsystemhavebeenstudiedin
[12, 14, 35] for the isothermal model, whereas the nonlinear case is consideredin [24]. The large
time behavior of the isothermal drift–diffusion system (1)–(6) has been studied in [13]. It has
beenproventhatthe solutionto the transientsystemconvergesto the thermalequilibriumstate
as t + if the boundary conditions (5) are in thermal equilibrium. This result was extended
→ ∞
to the degenerate case with nonlinear diffusivities in [26].
More precisely, the thermal equilibrium is a particular steady–state for which electron and hole
currents vanish, namely
r(N)+N Ψ= r(P) P Ψ=0.
−∇ ∇ −∇ − ∇
Theexistenceofathermalequilibriumhasbeenstudiedinthecaseofalinearpressurein[32,34]
and in the nonlinear case in [33]. Let us introduce the enthalpy function h defined by:
s r′(τ)
dτ,
τ
Z1
and the generalized inverse g of h, defined by:
h−1(s) if h(0+)<s< ,
g(s)= ∞
0 if s h(0+),
(cid:26) ≤
where we have implicitly assumed that h(+ ) = + . In the isothermal case, we simply have
∞ ∞
h=log and g =exp.
If the Dirichlet boundary conditions satisfy ND, PD >0 and
h(ND) ΨD =α and h(PD)+ΨD =α on ΓD, (16)
N P
−
the thermal equilibrium is defined for x Ω by:
∈
λ2∆Ψeq =g(α Ψeq) g(α +Ψeq)+C, (17)
P N
− − −
Neq =g(α +Ψeq), (18)
N
Peq =g(α Ψeq), (19)
P
−
with the boundary conditions (5)–(6).
As mentioned in Hypotheses 2, we furthermore assume in the linear case that the Dirichlet
boundary conditions satisfy the mass action law (13): NDPD =1. At the thermal equilibrium,
R(Neq,Peq) = 0 must hold, which implies in view of the form of R (12) that NeqPeq = 1, and
finally that α +α =0.
N P
The proof of convergenceto the thermal equilibrium is based on the entropy method, described
forinstanceinthereviewpaper[1]. ThismethodconsistsoflookingforanonnegativeLyapunov
functional,calledentropy,anditsnonnegativeproduction,connectedwithin anentropy–entropy
production estimate. It provides the convergence in relative entropy of the evolutive solution
towards the equilibrium state. Moreover, if the relative entropy is controlled with the entropy
production,one cancompute a convergencerate. This method hasbeen widely appliedto many
different systems ; see for instance [2] for Fokker–Planck type equations, or [6] for degenerate
parabolic problems.
Here the relative entropy functional is the deviation of the total energy (sum of the internal
energies for the electron and hole densities and the energy due to the electrostatic potential)
4
from the thermal equilibrium:
E(t)= H(N(t)) H(Neq) h(Neq)(N(t) Neq)
− − −
ZΩ(cid:18)
+H(P(t)) H(Peq) h(Peq)(P(t) Peq)
− − −
λ2
+ (Ψ(t) Ψeq)2 dx, (20)
2 |∇ − |
(cid:19)
x
with H(x)= h(s)ds, and the entropy production functional is given by
1
R
I(t)= N (h(N) Ψ)2+P (h(P)+Ψ)2 dx
|∇ − | |∇ |
ZΩ(cid:18) (cid:19)
+ R(N,P)(h(N)+h(P) h(Neq) h(Peq))dx. (21)
− −
ZΩ
In the nonlinear case, we assume that R=0, andthen the lastterm of I(t) vanishes, whereasin
the linear case with recombination–generationrate of the form (12), we obtain that
I(t)= N (log(N) Ψ)2+P (log(P)+Ψ)2 dx+ R (N,P)(NP 1)log(NP),
0
|∇ − | |∇ | −
ZΩ(cid:18) (cid:19) ZΩ
where the last term is clearly nonnegative.
The entropy–entropy production inequality writes:
t
0 E(t)+ I(s)ds E(0). (22)
≤ ≤
Z0
Exponentialdecaytowardsthethermalequilibriumforthedrift–diffusionmodelhasbeenproved
in[14]fortheBoltzmannstatistics,andextendedtotheFermi-Diracstatisticsin[15]. In[13],the
long–timebehaviorofthemodelwithmagneticfieldisstudied. Large–timebehaviorofreaction–
diffusion systems for a finite number of charged species has been investigated in [16, 21, 20].
Finally, the convergence towards the thermal equilibrium state for the drift–diffusion system in
the nonlinear degenerate case is proved in [26], but without any rate.
1.3 Outline of the paper
The outline of the paper is as follows. In Section 2, we introduce the discrete framework. It in-
cludesthedescriptionoftheconsiderednumericalschemesaswellasthedefinitionofthediscrete
relative entropy and the corresponding entropy production. We conclude this section by estab-
lishing a technical property of the generalized Scharfetter–Gummel fluxes in Lemma 1. Section
3is devotedto the detailedproofofourmainresult,whichisstatedinTheorem1. Thedecayof
numerical solutions towards the equilibrium is studied either under Hypotheses 2 corresponding
totheisothermalcasewithnonzerorecombination–generationrate,orunderHypotheses3corre-
spondingtoarathergeneralnonlinearpressure,withoutrecombination–generationrate. Wealso
assume existence and uniform L∞ estimates for numerical solutions to establish the large–time
behaviorofthescheme. TheseassumptionsarethendiscussedinSection4,wherewedistinguish
the nonlinear case from the isothermal case. It appears that the uniform-in-time L∞ estimates
needed to prove Theorem 1 are only obtained in the case of a zero doping profile. However in
the last section, we present some numerical results and observe an exponential convergence to
steady-state even when this condition is not satisfied.
5
2 Presentation of the discrete setting
2.1 Definition of the numerical schemes
In this subsection, we present the finite volume schemes for the time evolution drift–diffusion
system (1)–(6) and for the thermal equilibrium (17)–(19). The mesh = ( , , ) of the
M T E P
domain Ω is given by a family of open polygonal (or polyhedral in 3-D) control volumes, a
T
family of edges (or faces), and a family =(x ) of points. As it is classical in the finite
K K∈T
E P
volumediscretizationofdiffusivetermswithtwo-pointsfluxapproximations,weassumethatthe
mesh is admissible in the sense of [10, Definition 9.1]. It implies that the straight line between
two neighboring centers of cells (x ,x ) is orthogonalto the edge σ =K L.
K L
|
In the set of edges , we distinguish the interior edges σ = K L and the boundary edges
int
E | ∈E
σ . Within the exterioredges,we distinguish the Dirichletboundary edges included inΓD
ext
∈E
from the Neumann boundary edges included in ΓN: = D N . For a control volume
Eext Eext ∪Eext
K , we define the set of its edges, which is also split into = D N .
∈ T EK EK EK,int∪EK,ext∪EK,ext
Foreachedgeσ ,thereexistsatleastonecellK suchthatσ ,whichwillbedenoted
K
∈E ∈T ∈E
K . In the case where σ =K L , K can be either equal to K or L.
σ int σ
| ∈E
For all σ , we define d = d(x ,x ) if σ = K L and d = d(x ,σ) if σ , with
σ K L int σ K ext
∈ E | ∈ E ∈ E
σ . Then the transmissibility coefficient is defined by τ =m(σ)/d , for all σ .
K σ σ
∈E ∈E
We assume that the mesh satisfies the following regularity constraint:
ξ >0 such that d(x ,σ) ξd , K , σ . (23)
K σ K
∃ ≥ ∀ ∈T ∀ ∈E
Let ∆t>0 be the time step. We assume that there exists ∆t >0 such that
max
0 ∆t ∆t . (24)
max
≤ ≤
We set N = E(T/∆t) and tn = n∆t for all 0 n N . The size of the mesh is defined by
T T
≤ ≤
size( )=max diam(K), and we denote by δ =max(∆t,size( )) the size of the space–time
K∈T
T T
discretization.
Afinitevolumeschemeforaconservationlawwithunknownuprovidesavectoru =(u )
T K K∈T
Rθ (withθ =Card( ))ofapproximatevaluesandtheassociatepiecewiseconstantfunction,sti∈ll
T
denoted u :
T
u = u 1 ,
T K K
K∈T
X
where 1 denotes the characteristic function of the cell K. However, since there are Dirichlet
K
conditions on a part of the boundary, we also need to define approximate values for u at the
corresponding boundary edges: u = (u ) RθD (with θD = Card( D )). Therefore,
ED σ σ∈EeDxt ∈ Eext
the vector containing the approximate values both in the control volumes and at the Dirichlet
boundary edges is denoted by u = (u ,u ). We denote by X( ) the set of the discrete
M T ED M
functions u and by X ( ) the subset of X( ) of the functions vanishing at the boundary
M 0
M M
ΓD: u =(u ,0 ).
M T ED
For any vector u =(u ,u ), we define for all K and all σ
M T ED ∈T ∈EK
u if σ =K L ,
L K,int
| ∈E
u = u if σ D ,
K,σ  σ ∈EK,ext
u if σ N ,
 K ∈EK,ext
and 
Du =u u , D u= Du .
K,σ K,σ K σ K,σ
− | |
6
We can now define the discrete L2 norm and the discrete H1-seminorm on X( ) by
1,M
|·| M
u 2 = τ (D u)2, u X( ),
| M|1,M σ σ ∀ M ∈ M
σ∈E
X
u 2 = u 2 = m(K)(u )2, u X( ),
k Mk0,M k Tk0 K ∀ M ∈ M
K∈T
X
where is the classical L2 norm for piecewise constant functions.
0
k·k
LetusrecallinProposition1thediscretecounterpartofthePoincar´einequalityforpiecewise
constant functions. We refer to [4, Theorem 4.3] for a proof of this result in the case of mixed
boundary conditions.
Proposition 1. Let Ω be an open bounded polyhedral domain of Rd (d 1) and ∂Ω=ΓD ΓN
≥ ∪
with ΓD ΓN = and m(ΓD) > 0. Let = ( , , ) an admissible finite volume mesh of Ω
∩ ∅ M T E P
which satisfies (23). There exists a constant C depending only on Ω such that
P
C
P
u u u X ( ). (25)
k Mk0,M ≤ ξ1/2| M|1,M ∀ M ∈ 0 M
The scheme for the transient model. We have to define at each time step 0 n N
T
≤ ≤
the approximate solution un = (un) for u = N, P, Ψ and the approximate values at the
T K K∈T
boundaryun =(un) (whichinfactdoesnotdepend onnsince the boundarydatado not
ED σ σ∈EeDxt
depend on time). First of all, we discretize the initial and boundary conditions:
1
N0,P0 = (N (x),P (x)) dx, K , (26)
K K m(K) 0 0 ∀ ∈T
ZK
(cid:0) (cid:1) 1
ND,PD,ΨD = ND(γ),PD(γ),ΨD(γ) dγ, σ D , (27)
σ σ σ m(σ) ∀ ∈Eext
Zσ
(cid:0) (cid:1) (cid:0) (cid:1)
and we define
Nn =ND, Pn =PD, Ψn =ΨD, σ D , n 0. (28)
σ σ σ σ σ σ ∀ ∈Eext ∀ ≥
Then, we consider a backward Euler in time and finite volume in space discretization of the
drift–diffusion system (1)–(3). The schemes writes:
Nn+1 Nn
m(K) K − K + n+1 = m(K)R(Nn+1,Pn+1), K , n 0, (29)
∆t FK,σ − K K ∀ ∈T ∀ ≥
σX∈EK
Pn+1 Pn
m(K) K − K + n+1 = m(K)R(Nn+1,Pn+1), K , n 0, (30)
∆t GK,σ − K K ∀ ∈T ∀ ≥
σX∈EK
λ2 τ DΨn =m(K)(Pn Nn +C ), K , n 0. (31)
− σ K,σ K − K K ∀ ∈T ∀ ≥
σX∈EK
It remains to define the numerical fluxes n+1 and n+1 which are approximations of
FK,σ GK,σ
( r(N)+N Ψ) ν and ( r(P) P Ψ) ν onthe interval[tn,tn+1). We choose
K,σ K,σ
−∇ ∇ · −∇ − ∇ ·
Zσ Zσ
to discretizesimultaneouslythe diffusive partandthe convectivepartofthe fluxes. Inthe linear
case r = Id, we use the classical Scharfetter–Gummel fluxes. For all K , for all σ , we
K
∈T ∈ E
7
set:
n+1 =τ B DΨn+1 Nn+1 B DΨn+1 Nn+1 , (32)
FK,σ σ − K,σ K − K,σ K,σ
h (cid:16) (cid:17) (cid:16) (cid:17) i
n+1 =τ B DΨn+1 Pn+1 B DΨn+1 Pn+1 , (33)
GK,σ σ K,σ K − − K,σ K,σ
h (cid:16) (cid:17) (cid:16) (cid:17) i
where B is the Bernoulli function defined by
x
B(0)=1 and B(x)= x=0. (34)
exp(x) 1 ∀ 6
−
These fluxes, introduced in [23, 38], are widely used to approximate the drift–diffusion system
in the linearcase. They aresecondorderaccuratein space [31] andthey preservesteady–states.
In [13], the dissipativity of the implicit Scharfetter–Gummel scheme is proved. A proof of the
exponential decay of the free energy to its equilibrium value is given in [20] for an implicit time
discretizationofelectro–reaction–diffusionproblems,andthisresultisextendedtoafulldiscrete
problem in [18, 19]. In [17], some bounds for discrete steady–states solutions obtained with the
Scharfetter–Gummelschemeareestablished. Moreover,adiscreteanalogoftheentropy–entropy
production inequality (22) is proved in [8], yielding the long–time behavior of the Scharfetter–
Gummel scheme for the linear drift–diffusion system.
In the case of a nonlinear pressure function r satisfying (14), we use a generalization of the
Scharfetter–Gummel fluxes defined in [3]. For all K , for all σ , we set:
K
∈T ∈E
DΨn+1 DΨn+1
n+1 =τ dr(Nn+1,Nn+1) B − K,σ Nn+1 B K,σ Nn+1 ,
FK,σ σ K K,σ " dr(NKn+1,NKn+,σ1)! K − dr(NKn+1,NKn+,σ1)! K,σ#
(35)
DΨn+1 DΨn+1
n+1 =τ dr(Pn+1,Pn+1) B K,σ Pn+1 B − K,σ Pn+1 , (36)
GK,σ σ K K,σ " dr(PKn+1,PKn+,σ1)! K − dr(PKn+1,PKn+,σ1)! K,σ #
where
h(b) h(a)
− if a, b>0, a=b,
log(b) log(a) 6
dr(a,b)= − (37)
a+b
 r′ elsewhere.
2
(cid:18) (cid:19)
This particular choice of drensures the preservation of the thermal equilibrium at the discrete
level, which is crucial to have a good long time behavior. We notice that in the isothermal
case r = Id, we recover exactly the classical Scharfetter–Gummel fluxes (32) and (33). These
fluxeswerealsogeneralizedtoalargerclassofstatisticaldistributionfunctionsarisinginorganic
semiconductors modeling [30].
Other extensions of the Scharfetter–Gummel scheme have been proposed. A scheme valid in
the case where both convective and diffusive terms are nonlinear is studied in [9] and applied
to organic semiconductors models in [28, 29], but this method leads to solve a nonlinear elliptic
problem at each interface. We also mention [25, 27] where fluxes (35)–(36) are considered, but
with another definition of dr which does not allow to preserve the thermal equilibrium at the
discrete level. A finite volume scheme preserving the long–time behavior of the solutions of the
nonlinear drift–diffusion model is introduced in [7], and the convergence of this scheme towards
the equilibrium state is proved, based on the control of the discrete energy production.
8
The scheme for the thermal equilibrium. We compute an approximation(Neq,Peq,Ψeq)
T T T
of the thermal equilibrium defined by (17)–(19) with the following finite volume scheme:
λ2 τ DΨeq =m(K)(g(α Ψeq) g(α +Ψeq)+C ), K , (38)
− σ K,σ P − K − N K K ∀ ∈T
σX∈EK
Neq =g(α +Ψeq), K , (39)
K N K ∀ ∈T
Peq =g(α Ψeq), K . (40)
K P − K ∀ ∈T
Existence and uniqueness of solution to this nonlinear scheme is studied in [7].
The discrete entropy and entropy production For n N, the discrete relative entropy
∈
functional is defined by:
En = m(K)(H(Nn) H(Neq) h(Neq)(Nn Neq)
K − K − K K − K
K∈T
X
+H(Pn) H(Peq) h(Peq)(Pn Peq))
K − K − K K − K
λ2
+ Ψn Ψeq 2 (41)
2 | M− M|1,M
and the discrete entropy production functional is defined by:
In = τ min(Nn,Nn )(D (h(Nn) Ψn))2
σ K K,σ σ −
σX∈E h
(K=Kσ)
+ min(Pn,Pn )(D (h(Pn)+Ψn))2
K K,σ σ
+ m(K)R(Nn,Pn)[h(Nn)+h(Pni) h(Neq) h(Peq)], (42)
K K K K − K − K
K∈T
X
where means a sum over all the edges σ and K inside the sum is replaced by K
σ
∈ E
σ∈E
X
(K=Kσ)
(therefore σ is an edge of the cell K =K ).
σ
As already mentioned in the continuous framework, the last term in the definition of In is
reduced to zero in the nonlinear case where R = 0 and is nonnegative in the linear case with a
recombination–generationrate under the form (12).
2.2 Properties of the numerical fluxes
In Section 3, we study the large time behavior of the approximate solution by adapting the
entropy–dissipation method to the discrete level. We start by establishing a lemma which will
be useful to provethe discretecounterpartof (22). It is a generalizationto the nonlinearcase of
Corollary A.2 in [5].
Lemma 1. For all K and all σ , the fluxes n+1 and n+1 defined by (35) and (36)
∈T ∈EK FK,σ GK,σ
satisfy:
n+1D h(Nn+1) Ψn+1 τ min(Nn+1,Nn+1) D h(Nn+1) Ψn+1 2, (43)
FK,σ − K,σ ≤ σ K K,σ σ −
n+1D (cid:0)h(Pn+1)+Ψn+1 (cid:1) τ min(Pn+1,Pn+1) D(cid:0) h(cid:0)(Pn+1)+Ψn+1 (cid:1)2(cid:1). (44)
GK,σ K,σ ≤ σ K K,σ σ
(cid:0) (cid:1) (cid:0) (cid:0) (cid:1)(cid:1)
9
Proof. We focus on the proof of (43) because (44) will be proved in the same way, replacing N
by P and Ψ by Ψ. Let us define
−
n+1 = n+1D h(Nn+1) Ψn+1 τ min(Nn+1,Nn+1) D h(Nn+1) Ψn+1 2,
RK,σ FK,σ − K,σ− σ K K,σ σ −
=D h(Nn(cid:0)+1) Ψn+1 (cid:1) n+1 τ min(Nn+1,Nn+1(cid:0))D h(cid:0)(Nn+1) Ψn+1 (cid:1)(cid:1) .
− K,σ FK,σ − σ K K,σ − K,σ
(cid:0) (cid:1) (cid:16) (cid:0) (cid:1) (cid:17)
Our aim is to prove that n+1 0. We may write:
RK,σ ≤
Dh(Nn+1) DΨn+1
D h(Nn+1) Ψn+1 =dr(Nn+1,Nn+1) K,σ K,σ .
− K,σ K K,σ dr(NKn+1,NKn+,σ1) − dr(NKn+1,NKn+,σ1)!
(cid:0) (cid:1)
Using the fact that x=B( x) B(x) for all x R, it rewrites:
− − ∈
Dh(Nn+1)
D h(Nn+1) Ψn+1 =dr(Nn+1,Nn+1) B − K,σ
− K,σ K K,σ dr(NKn+1,NKn+,σ1)!
(cid:0) (cid:1)
Dh(Nn+1) DΨn+1 DΨn+1
B K,σ B − K,σ +B K,σ .
− dr(NKn+1,NKn+,σ1)!− dr(NKn+1,NKn+,σ1)! dr(NKn+1,NKn+,σ1)!!
But, the definition (37) of dr implies:
Dh(Nn+1)
K,σ =Dlog(Nn+1)
dr(Nn+1,Nn+1) K,σ
K K,σ
and the definition (34) of B ensures:
B Dlog(Nn+1) Nn+1 B Dlog(Nn+1) Nn+1 =0.
− K,σ K − K,σ K,σ
Therefore, we obtai(cid:0)n: (cid:1) (cid:0) (cid:1)
n+1 =τ D h(Nn+1) Ψn+1 dr(Nn+1,Nn+1)
RK,σ σ − K,σ K K,σ
(cid:0) DΨn+1 (cid:1) Dh(Nn+1)
B − K,σ B − K,σ Nn+1 min(Nn+1,Nn+1)
×" dr(NKn+1,NKn+,σ1)!− dr(NKn+1,NKn+,σ1)!! K − K K,σ
(cid:16) (cid:17)
DΨn+1 Dh(Nn+1)
B K,σ B K,σ Nn+1 min(Nn+1,Nn+1) .
− dr(NKn+1,NKn+,σ1)!− dr(NKn+1,NKn+,σ1)!! K,σ − K K,σ #
(cid:16) (cid:17)
Since B is nonincreasing on R and dr(a,b) 0 for all a, b 0, we conclude that n+1 0.
≥ ≥ RK,σ ≤
3 Exponential decay to the discrete thermal equilibrium
Inthis section,weestablishthe mainresultofthis article,namelythedecayrateofapproximate
solutions given by the Scharfetter–Gummel scheme towards an approximation of the thermal
equilibrium. Assumptions concerning existence and L∞ estimates will be discussed in the next
section. The main theorem is the following:
Theorem1(Exponentialdecay). LetHypotheses1befulfilledwithm>0. Let =( , , )be
M T E P
an admissible mesh of Ω satisfying (23) and ∆t be the time step verifying (24). We also assume
10