Philippe Bechouche

Laurent Gosse
email: l.gosse@ba.iac.cnr.it

A SEMICLASSICAL COUPLED MODEL FOR THE TRANSIENT SIMULATION OF SEMICONDUCTOR DEVICES

Keywords: Drift-diffusion equation, Schrödinger-Poisson equation, WKB ansatz, open quantum system, Robin boundary condition AMS subject classifications. 76W05, 65J10

We consider the approximation of a microelectronic device corresponding to a n + -nn + diode consisting in a channel flanked on both sides by two highly doped regions. This is modelled through a system of equations: ballistic for the channel and drift-diffusion elsewhere. The overall coupling stems from the Poisson equation for the self-consistent potential. We propose an original numerical method for its processing, being realizable, explicit in time and nonnegativity preserving on the density. In particular, the boundary conditions at the junctions express the continuity of the current and don't destabilize the general scheme. At last, efficiency is shown by presenting results on test-cases of some practical interest.

1. Introduction and modelling.

1.1. Preliminaries. The present paper scopes onto developing a fast, efficient and stable approach to 1D semiconductor components computations without directly processing any kind of kinetic equation; only moment-type approximations are to be used, in "high-field regime" though, in order to cope with "submicronic standards".

It is a well-known fact nowadays that a whole hierarchy of mathematical models ranges from quantum mechanics to transport equations eventually leading to "fluidlike" models, see e.g. [START_REF] Ben Abdallah | On a hierarchy of macroscopic models for semi-conductors[END_REF][START_REF] Poupaud | Mathematical theory of kinetic equations for transport modelling in semiconductors[END_REF]. Two main tools exist for performing the first step, namely WKB and Wigner measures methods, see [START_REF] Ashcroft | Solid-state physics[END_REF][START_REF] Bechouche | Semiclassical limit for the Schrödinger-Poisson equation in a crystal[END_REF][START_REF] Keller | Semiclassical mechanics[END_REF][START_REF] Lions | Sur les mesures de Wigner[END_REF]; the correct way to interpreting the WKB classical system being through a Vlasov equation, see [START_REF] Gosse | Using K-branch entropy solutions for multivalued geometric optics computations[END_REF][START_REF] Gosse | Multiphase semiclassical approximation of an electron in a onedimensional crystalline lattice -I. Homogeneous problems[END_REF]. It is then at this "mesoscopic scale" that one introduces collisions between electrons and impurities (i.e. phonons, quantum representations of the crystal's vibrations) which thermalize the carriers population. This is rendered through a collision operator put on the right-hand side of the kinetic equation which is now meant to relax towards a local equilibrium called Maxwellian by analogy with Boltzmann's theory of rarefied gases, [START_REF] Degond | Mathematical modelling of microelectronics semiconductor devices[END_REF]. The passage from kinetic models to coarser ones always stems from prescribing a peculiar dependence of the kinetic density in its velocity variable; such simplifications can sometimes be justified by means of convenient time/space rescalings, consult [START_REF] Poupaud | Mathematical theory of kinetic equations for transport modelling in semiconductors[END_REF] for a survey on these questions.

At this level, we take advantage of former results, mainly [START_REF] Cercignani | A drift-collision balance for a Boltzmann-Poisson system in boundad domains[END_REF][START_REF] Carrillo | Computational macroscopic approximations to the 1-D relaxation time kinetic system for semiconductors[END_REF][START_REF] Arnold | Low and high field scaling limits for the Vlasov and Wigner-Poisson-Fokker-Planck systems[END_REF]; in these papers, coexistence of several qualitatively different areas has been clearly evidenced inside a standard n + nn + diode, [START_REF] Blotekjaer | Transport equations for electrons in two-valley semiconductors[END_REF], under reasonable voltage biases. To each of these regions can be associated a particular rescaling of the kinetic transport model; roughly speaking, heavily doped zones (the source and the drain) exhibit "low-field" diffusive behaviour whereas inside the channel takes place "high-field" ballistic transport. To some extent, a schematic but acceptable picture goes as follows: carriers move inside the channel according to monokinetic solutions of the Vlasov-Poisson equation (because they accelerate, no shock is likely to appear at least at steady state) whereas they are "nearly Maxwellian" inside the source/drain regions where their position densities can be computed by means of a drift-diffusion approximation. We obviously expect to observe stabilization of transient dynamics leading to convergence onto some steady time-asymptotic regime (see e.g. [START_REF] Degond | A coupled Schrödinger Drift-Diffusion model for quantum semiconductor device simulations[END_REF][START_REF] Ringhofer | Computational methods for semiclassical and quantum transport in semiconductor devices[END_REF] for direct steady-state computations). We stress that our goal here is to devise a time-marching strategy.

Such a program clearly raises up an issue in terms of time-stepping compatibility, and this whatever the boundary conditions on each side of the channel. Indeed, any drift-diffusion equation, being of parabolic type, asks for a CFL restriction like ∆t = O(∆x 2 ) whereas ballistic regime is essentially of hyperbolic nature (signals travel at finite speed), thus requiring ∆t = O(∆x) instead. Of course, imposing the smallest value of ∆t everywhere, despite ensuring stability, would result in an overwhelming smearing of the sharp profiles inside the channel. We thus selected a mixed implicit/explicit numerical strategy as described in §4.1. The channel being flanked on each side by low-field areas, deriving appropriate boundary conditions connecting these regions was necessary, so as oscillations won't develop. It has been observed that Dirichlet conditions aren't conclusive as they didn't allow for the scheme to relax toward a steady-state; however, imposing the current's continuity at the junctions leads to a nonlinear Robin boundary condition, [START_REF] Arrieta | Parabolic problems with nonlinear boundary conditions and critical nonlinearities[END_REF][START_REF] Evans | Regularity properties for the heat equation subject to nonlinear boundary constraints[END_REF][START_REF] Favini | Nonlinear boundary conditions for nonlinear second order operators on C[0,1[END_REF][START_REF] Goldstein | Smoothing for nonlinear parabolic equations with nonlinear boundary conditions[END_REF], for the drift-diffusion equations: it also remains stable and easy to implement, see §4.3. Obviously Dirichlet conditions are imposed on each side of the device and for the Poisson equation.

The remainder of the paper is devoted to numerical experiments as an outcome of the aforementioned general framework. A careful rescaling of a well-known GaAs device has been performed in §3 and practical results are displayed in §5. In particular, we observe stabilization in time toward a steady profile together with a satisfying residues decay. We also show the way current settles inside the device for large times; a reasonably constant current being the main goal in order to achieve reliable intensity/voltage diagrams (see especially Fig. 5.6 in §5.5 where an original derivation leading to a real improvement of ballistic currents is proposed). Anyway, we aren't aware of any time-marching scheme able to stabilize onto steady-states with perfectly flat currents for such models. Presently, we found out Ohm's law, stipulating a linear dependence between currents and biases, being satisfied as soon as potential drops are strong enough. This was to be somehow expected as we developed the present approach for a high-field setting in which low biases don't really fit.

1.2. From quantum models to kinetic equations. We briefly recall the different levels for the modelling of charge carriers; the finest one here is nonrelativistic quantum mechanics. One considers thus a mixed-state, i.e. a statistical repartition of quantum states ψ ℓ ∈ C indexed by the (quantized) energy levels solution of the one-dimensional Schrödinger-Poisson equation in the parabolic band approximation: (m * stands for the effective mass of the electron in the material under consideration)

i ∂ t ψ ℓ + 2 2m * ∂ xx ψ ℓ = eU (t, x)ψ ℓ , ε s ∂ xx U = e(̺ -D). (1.1)
The wave function ψ ℓ (t, x) corresponds to an energy level E ℓ of the homogeneous Hamiltonian, , e and m * stand for Planck's constant, the charge and the effective mass of the electron in the considered material respectively. The Hartree equation renders a mean-field approximation of the Coulombian forces, with ε 0 being the dielectric permittivity of the medium. D is the given doping concentration and ̺(t, x) is electrons' density, which reads,

̺(t, x) = ℓ F (E ℓ )|ψ ℓ (t, x)| 2 , with e.g. F (E) = (1 + exp((E -E F )/k B θ 0 )) -1 in
case of Fermi-Dirac statistics, E F being usually called the Fermi energy1 . θ 0 is the ambient temperature, and k B the Boltzmann's constant. The sum ranges therefore over the whole spectrum of the Hamiltonian consisting in an orthonormal basis of L 2 ; hence (1.1) is indeed a system of equations all coupled by the Poisson term U . With the next level of modelling comes into play the Wigner transform:

w(t, x, p) = 1 2π R z t, x + y 2 , x - y 2 exp(iyp).dy,
where

z(t, x, y) = ℓ F (E ℓ )ψ ℓ (t, x)ψ ℓ (t, y).
One observes at once that:

̺(t, x) = R w(t, x, p).dp = z(t, x, x).
The Wigner transform solves the (pseudo-differential) Wigner equation,

∂ t w + p m * ∂ x w -eΘ[U ]f = 0, Θ[V ] = i U t, x + 2i ∂ x -U t, x -2i ∂ x ,
completed by the Poisson equation as in (1.1). Passing to the (semi)classical limit → 0, w ⇀ f (t, x, ξ) ≥ 0 in the sense of measures [START_REF] Lions | Sur les mesures de Wigner[END_REF] and the pseudo-differential operator Θ becomes the standard Vlasov term ∂ x U.∂ ξ f . It is usually at this intermediate level that electron-phonon collisions are introduced in the models (see however [START_REF] Ashcroft | Solid-state physics[END_REF][START_REF] Gosse | Multiphase semiclassical approximation of the one-dimensional harmonic crystal. I. The periodic case[END_REF]) by means of an integral operator on the right-hand side:

∂ t f + ξ∂ x f + e m * ∂ x U ∂ ξ f = Q(f ), ε s ∂ xx U = e(̺ -D), E = -∂ x U. (1.2)
We shall not enter the details of Q(f ) since different models exist like e.g. BGK or Fokker-Planck; however, they share the interesting property of admitting the same Maxwellian distribution, a fact that will be extensively used in the next section. At this level of description the electron and current density are the moments of f :

̺(t, x) = R f (t, x, ξ)dξ, J(t, x) = R ξf (t, x, ξ)dξ.
2. Low and high field scalings inside a n + nn + device. Our 1D problem combines two distinct zones according to the importance which is given to the phenomena rendered through Q in the equation (1.2); purely ballistic transport means Q ≡ 0. However, we begin with the opposite case and mainly follow [START_REF] Carrillo | Computational macroscopic approximations to the 1-D relaxation time kinetic system for semiconductors[END_REF] in our presentation. Let us first recall some conventions: L > 0 stands for the characteristic device's length whereas ̺ 0 ≥ 0 and θ 0 = k B T /m * ≥ 0, for an average value of the doping profile D(x) and the lattice temperature. U th = k B T /e is the thermal potential. Four (possibly) different velocity magnitudes appear: v th = √ θ 0 is the thermal velocity. If τ 0 stands for an average value of the relaxation times2 τ (x) inside the device, v relax = L/τ 0 is the relaxation velocity. Similarly, if [Φ] is the bias applied at the contacts, one defines the reference magnitude for the drift velocity as follows,

v drif t = τ 0 e[Φ] m * L , (2.1)
whereas the ballistic velocity reads:

v bal = 2e[Φ] m * .
In a 1D framework, it is reasonable to render collisions by mean of a linear relaxation time operator thus we can start from the following kinetic equation to describe the dynamic of the electrons in the device.

∂ t f + ξ∂ x f - e m * E(t, x)∂ ξ f = 1 τ M(θ 0 )̺(x) -f , (2.2)
where the Maxwelian reads

M(θ 0 ) = 1 √ 2πθ 0 exp - ξ 2 2θ 0 θ0→0 → δ,
where δ stands for the Dirac mass in zero; the electric field E(t, x) is coupled to the Poisson equation by

ε s ∂ xx U = e(̺(x) -D(x)), E(t, x) = -∂ x U. (2.3) Remark 1.
Starting from here, we assume that junctions are rendered through sharp interfaces inside the device. Some quantities like e.g. the temperature will obviously be discontinuous across them. This modelling has already been widely used, see for instance [START_REF] Arnold | Low and high field scaling limits for the Vlasov and Wigner-Poisson-Fokker-Planck systems[END_REF][START_REF] Benabdallah | Coupling one-dimensional time-dependent classical and quantum transport models[END_REF][START_REF] Carrillo | Computational macroscopic approximations to the 1-D relaxation time kinetic system for semiconductors[END_REF][START_REF] Cercignani | Device benchmark comparisons via kinetic, hydrodynamic and high-field models[END_REF][START_REF] Degond | A coupled Schrödinger Drift-Diffusion model for quantum semiconductor device simulations[END_REF].

Drift-diffusion approximation.

It is convenient to normalize the doping concentration and the potential bias: D(x) = ̺ 0 D(x) and U (t, x) = [Φ] Û ( t, x) where we make use of rescaled variables defined as follows:

x = Lx, ξ = v th ξ, t = v drif t t L , τ (x) = τ 0 τ (x).
Moreover the kinetic density rewrites:

f (t, x, ξ) = ̺ 0 v th f ( t, x, ξ).
The main assumption in this approximation is v drif t ≪ v th ≪ v relax ; hence two small parameters naturally appear:

ǫ = v th v relax , η = v drif t v th (2.4)
which we suppose to be of the same order of magnitude. Then "low-field" is precisely meant to be the limit η ≃ ǫ → 0 thus (2.2) rewrites:

ǫ∂ t f + ξ∂ x f - e m * E(t, x)∂ ξ f = 1 ǫτ (M(1)̺ -f ). (2.5)
This regime describes the flow in highly doped areas (source and drain); (2.5) admits asymptotic expansions which allow to compute the two first moments of f in closed form. Since collisions are dominant, drift velocity is much smaller than thermal velocity and the mean free path can be considered infinitesimal, i.e. o [START_REF] Arrieta | Parabolic problems with nonlinear boundary conditions and critical nonlinearities[END_REF]. Hilbert expansion leads to:

∂ t ̺ + ∂ x J = 0, J = -d n (x)∂ x ̺ -µ n (x)E̺. (2.6)
There, µ n (x) = eτ (x)/m * is the electron mobility and d n (x) = τ (x)θ 0 = µ n U th ≥ 0 appears to be a diffusion coefficient. In this regime, it is tacitly assumed that the solution of (2.5) is "nearly Maxwellian", i.e.

f (t, x, ξ) = ̺(t, x)M θ0 (ξ) + o(ǫ).
2.2. Shur-Eastman/WKB ansatz for ballistic regime. This is somehow the opposite; we consider now the following variables:

x = Lx, ξ = v bal ξ, t = v bal t L , τ (x) = τ 0 τ (x).
The kinetic density reads accordingly,

f (t, x, ξ) = ̺ 0 v bal f ( t, x, ξ).
The main assumption in this approximation is v th ≪ v bal ≃ v drif t , which leads to the limit β = v th /v bal → 0, while α = v drif t /v bal = O(1). Thus comes out,

∂ t f + ξ∂ x f + e m * ∂ x U.∂ ξ f = 1 2ατ Q(f ), θ 0 → 0, ∂ xx U = γ(̺ -D). (2.7)
Clearly, neither Hilbert nor Chapman-Enskog expansions are likely to lead to a diffusion-like equation. Instead, one obtains formally a Vlasov-Poisson system for which the kinetic density is of the monokinetic form at zero temperature (see [START_REF] Shur | Near ballistic transport in GaAs devices at 77K[END_REF] and [START_REF] Arnold | Low and high field scaling limits for the Vlasov and Wigner-Poisson-Fokker-Planck systems[END_REF], §5). This can be approximated by the so-called "WKB system", as in [START_REF] Gosse | Convergence results for an inhomogeneous system arising in various high frequency approximations[END_REF],

∂ t ̺ + ∂ x ∂ x ϕ m * ̺ = 0, ∂ t ϕ + 1 2m * |∂ x ϕ| 2 + eU = 0, (2.8)
assuming that the wave functions behave like √ ̺ exp(iϕ/ ), or equivalently, that the solution of (2.7) reads 3 and this holds in particular for the accelerating particles ∂ x u(t, x) ≥ 0 inside the channel. In §5 of [START_REF] Arnold | Low and high field scaling limits for the Vlasov and Wigner-Poisson-Fokker-Planck systems[END_REF], the authors suggest that the distribution f should be of the type δ(ξ -Eτ ), which boils down to assuming that u(t, x) ≃ τ (x)E(t, x) at least after some time. This is indeed what we observe in the numerical results (in constrast with [START_REF] Carrillo | Computational macroscopic approximations to the 1-D relaxation time kinetic system for semiconductors[END_REF] where E and u look quite different at steady-state); see §5.4.

f (t, x, ξ) = ̺(t, x)δ(ξ -∂ x ϕ(t, x)), see [31]. Of course, u = 1 m * ∂ x ϕ is meant to remain smooth

Rescaling the device.

First of all let us remark that the device has not the same physical characteristics in the channel (ballistic zone) where there is a low doping concentration relatively to the source and the drain (Drift-diffusion zone). This implies that the electron mobility is higher in the channel as in the source and drain. The same occurs for the relaxation time in the same proportion. More precisely, let µ 1 be the mobility in the source and drain and µ 2 in the channel. We have

µ 1 µ 2 = τ 1 τ 2 = v 2 v 1 (3.1)
if we suppose that τ 1 (resp.τ 1 ) is proportional to L/v 1 (resp. to L/v 2 ). This means that if we want to see phenomena evolve at the same velocity, we have to scale the time in this proportion. An other way to explain this is if we scale the drift-diffusion equation first (with mobility µ 1 ) , then we have to scale the ballistic equation with the same scaling obtained in the drift-diffusion equation but with the mobility µ 2 .

3.1. Rescaling the drift-diffusion equation. According to (2.6), the physical drift-diffusion model reads for x ∈ [0, L], t ≥ 0:

ε s ∂ xx U = e(̺ -D), ∂ t ̺ = ∂ x (µ 1 U th ∂ x ̺ -µ 1 ̺∂ x U ) . (3.2)
The first equation in (3.2) is Poisson's completed by a continuity equation. In the device's variables, the current density is given by: J

(t, x) = -d 1 ∂ x ̺ + µ 1 ̺∂ x U .
The first part of the current is the diffusion term (d 1 = µ 1 U th ) and the second part the drift term. The general way to process (3.2) starts by rescaling space and time variables according to the size of the device and some characteristic velocity which is to be chosen to convenience:

x = x L , t = t L/v 1 , (3.3)
We can deduce an adimensionalized density ρ and a potential Û as follows:

̺(t, x) = ̺ 0 ρ( t, x) U (t, x) = U 0 Û ( t, x). (3.4)
The doping concentration is obviously treated accordingly:

D(x) = ̺ 0 D(x).
Here ̺ 0 stands for a mean density inside the device (obtained by e.g. doping concentration's arithmetic or geometric average) and U 0 is a reference value for the potential. We first rescale the Poisson equation so as to obtain,

ε s U 0 e̺ 0 L 2 ∂ 2 Û ∂ x2 = ρ -D.
If we denote by [Φ 0 ] the applied bias, the scaled Debye length is given by

γ = ε s [Φ 0 ] e̺ 0 L 2 .
At this level, we decide to choose U 0 = e̺0L 2 εs = γ[Φ 0 ] in order to scale the coefficient to 1. Thus the scaled Poisson equation rewrites:

∂ 2 Û ∂ x2 = ρ -D. (3.5)
Let us now rescale the continuity equation in (3.2); the outcome reads:

∂ ρ ∂ t = ∂ ∂ x µ 1 U th v 1 L ∂ ρ ∂ x - µ 1 U 0 v 1 L ρ ∂ Û ∂ x . (3.6)
A relevant quantity is the ratio between the diffusion and drift current: here it is U th /U 0 . We observe that in general, it is not possible to scale to unit simultaneously the scaled Debye length, the drift and diffusion terms. So in the present framework, we chose to scale to unit the scaled Debye length and the scaled drift term. Let us choose v 1 such that the drift velocity coefficient scales to 1 i.e. µ1U0 v1L = 1. In such a way (3.6) transforms into:

∂ ρ ∂ t = ∂ ∂ x U th U 0 ∂ ρ ∂ x - ρ ∂ Û ∂ x . (3.7)
The scaled current density therefore reads

J( t, x) = - U th U 0 ∂ ρ ∂ x + ρ ∂ Û ∂ x , (3.8)
and (3.7)-(3.8) should be completed with the rescaled potential equation (3.5). A short computation with the numerical values given below shows that for high concentrations of doping, the drift term is more important than the diffusion term, for concentration of doping of 10 19 /m 3 , both terms are of the same order of magnitude, and for concentration of doping lower, the diffusion term is dominant.

Rescaling the Ballistic equation.

Let us now scale the equations (2.8) with the same scaling (3.3) in space and time but with v 2 instead of v 1 and (3.4) for the density and potential. Moreover we assume the phase scales in the following way ϕ(t, x) = ϕ 0 φ( t, x). (3.9) With this scaling equations (2.8) become

∂ ρ ∂ t + ∂ ∂ x ϕ 0 m * Lv 2 ∂ φ ∂ x ρ = 0 ∂ ∂ t ∂ ∂ x ϕ 0 φ m * Lv 2 + 1 2 ∂ ∂ x ϕ 2 0 (m * ) 2 L 2 v 2 2 |∂ x φ| 2 + eU 0 m * v 2 2 ∂ Û ∂ x = 0.
Note that we have derived w.r.t. x the second equation. In order to work with the same time scale as in the Drift-Diffusion case, we have to choose v 2 such that (3.1) is verified. Let let us set u = ϕ0 m * Lv2 ∂ x φ. We can choose ϕ 0 such that ϕ0 m * Lv2 = 1. We set for the scaled current Ĵ = ρu. Then the scaled WKB equation (2.8) becomes:

∂ ρ ∂ t + ∂ ∂ x (ρu) = 0, ∂ ∂ t u + 1 2 ∂ ∂ x u 2 + eU 0 m * v 2 2 ∂ Û ∂ x = 0.
From now on, we shall always work with scaled quantities thus it makes sense to remove the "hats". 2 ) = 0.2. 3.3. The kinetic equation. We finally rewrite the kinetic model (2.2); let us rescale the variables as follows:

x = Lx, ξ = v th ξ, t = v t L , τ (x) = τ 0 τ (x)
where v is a velocity that is to be set according to the scaling which has been chosen. Moreover the adimensioned kinetic density rewrites:

f (t, x, ξ) = ̺ 0 v th f ( t, x, ξ).
As usual, we define an adimensioned density ρ and a potential Φ as follows:

̺(t, x) = ̺ 0 ρ( t, x) U (t, x) = Φ 0 Φ( t, x). (3.10)
With this notation the adimensioned form of (2.2) reads:

∂ t f + v th v ξ • ∂ x f - U 0 U th v th v ∂ xΦ • ∂ ξ f = v relax v M(θ 0 )̺(x) -f . (3.11)
Accroding to Remark 2, v can be chosen to be e.g. between 100v th and 1000v th .

4.

A staggered scheme for the coupled model. Hereafter we always consider working with the rescaled equations derived in the preceding section.

4.1.

A semi-implicit discretization. We now aim at describing in full detail the numerical recipe we propose for the solving of the following adimensionalized problems:

(CL) ∂ t ̺ + ∂ x J = 0, J = -d∂ x ̺ -̺∂ x U, ̺(t, x = 0) = 1, x ∈]0, x L [; (B) ∂ t ̺ + ∂ x J = 0, ∂ t u + u∂ x u + ∂ x U = 0, J = ̺u, x ∈]x L , x R [; (CR) ∂ t ̺ + ∂ x J = 0, J = -d∂ x ̺ -̺∂ x U, ̺(t, x = 1) = 1, x ∈]x R , 1[;
where all these PDE's are coupled through the Poisson equation:

(P ) -∂ xx U = ̺ -D(x), U (t, x = 0) = 0, U (t, x = 1) = U 0 ≤ 0. x ∈]0, 1[.
We are given initial data ̺(t = 0, .) = ̺ 0 ≥ 0, D(x) ≥ 0 and u(t = 0, .) we usually take zero. We shall develop on the implementation of boundary conditions in x = x L , x R rendering the junctions in §3.3. First, let us say that (CL), (CR) will be solved implicitly whereas (B) will be processed explicitly in time. Another important feature is that ̺ and u won't be known at the same abscissae (hence the use of "staggered"). The computational domain is the half-stripe R + ×]0, 1[, which is parametrized by ∆x = 1 N > 0; we define,

x j = j∆x = j N ; j = 0, ..., N.
We also define the following domains:

B = {x Q , x Q+1 , ..., 1 -x Q }, x Q = x q + ∆x 2 , 1 -x Q = x N -q -∆x 2 , C L = {x 0 , x 2 , ..., x q }, C R = {x N -q , ..., x N }, Thus ̺ n=0 = (̺ 0 (x j )) j=0,...,N ∈ R N +1 +
. The boundary conditions read at the discrete level: ̺ 0,n = ̺(t n , 0) and ̺ N,n = ̺(t n , x N ) = ̺(t n , 1).

4.2.

Choosing the time-step ∆t. We select a standard centered discretization for (CL), (CR); it reads: (n the time index, j the space one)

̺ j,n+1 -∆t d(̺j+1,n+1-̺j,n+1)-d(̺j,n+1-̺j-1,n+1) ∆x 2 = ̺ j,n + ∆t ∆x F j+ 1 2 ,n , R j+ 1 2 ,n R 2 -F j-1 2 ,n , R j-1 2 ,n R 2 , (4.1)
where the following notations have been used:

F j+ 1 2 ,n := max 0, Uj+1,n-Uj,n ∆x , min 0, Uj+1,n-Uj,n ∆x , R j+ 1 2
,n := {̺ j,n , ̺ j+1,n } . Clearly, <, > R 2 stands for the standard R 2 scalar product. At this stage, we tacitly assume that ̺ n being known, we already have a way to produce the two borderline values ̺ q+1,n+1 and ̺ N -q-1,n+1 . Thus, the scheme boils down to inverting a tridiagonal q × q matrix reading:

M q =      1 + αd∆t ∆x 2 -d∆t ∆x 2 0 • • • -d∆t ∆x 2 1 + 2d∆t ∆x 2 -d∆t ∆x 2 0 0 . . . . . . . . . 0 -d∆t ∆x 2 1 + βd∆t ∆x 2     
, where α, β can be 1 or 2 depending whether (4.1) is meant to approximate (CL) or (CR) because of the boundary conditions. Anyway, it is clearly unconditionally invertible since it is strictly diagonally dominant; thus the discretization (4.1) will remain stable as long as its right-hand side will be a convex combination, i.e. for a time step ∆t = O(∆x). In order to produce ̺ q+1,n+1 and ̺ N -q-1,n+1 , we solve (B) by means of an explicit time-marching scheme already used in [START_REF] Gosse | Numerical approximation of linear one-dimensional conservation equations with discontinuous coefficients[END_REF][START_REF] Gosse | Convergence results for an inhomogeneous system arising in various high frequency approximations[END_REF], except that we set it up now together with a well-balanced strategy, as in [START_REF] Gosse | A well-balanced flux splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF] for the handling of the electric field. More precisely, from [START_REF] Gosse | Numerical approximation of linear one-dimensional conservation equations with discontinuous coefficients[END_REF],

̺ j,n+1 = ̺ j,n - ∆t ∆x A j+ 1 2 ,n , R j+ 1 2 ,n R 2 -A j-1 2 ,n , R j-1 2 ,n R 2 , (4.2)
with j = q, ..., N -q -1 and,

A j+ 1 2 ,n := max 0, u j+ 1 2 ,n , min 0, u j+ 1 2 ,n
.

Now the velocity values u j+ 1 2 ,n are sought on the staggered computational grid by means of a simple upwind well-balanced scheme:

u j+ 1 2 ,n+1 = u j+ 1 2 ,n - ∆t ∆x F(u j+ 1 2 ,n , ũj+ 3 2 ,n ) -F(ũ j-1 2 ,n , u j+ 1 2 ,n ) . (4.3) F(u, v) = 1
2 (max(u, 0) 2 + min(v, 0) 2 ) is the standard upwind scheme and the modified values ũj- 1 2 ,n , ũj+ 3 2 ,n are computed according to the steady-state equation, u∂ x u + ∂ x U = 0 out of resonance points, see Appendix A and [START_REF] Gosse | A well-balanced flux splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF] for details. We know from [START_REF] Gosse | A well-balanced flux splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF][START_REF] Gosse | Numerical approximation of linear one-dimensional conservation equations with discontinuous coefficients[END_REF][START_REF] Gosse | Convergence results for an inhomogeneous system arising in various high frequency approximations[END_REF] that this approach for the solving of (B) will be stable for a time step ∆t = O(∆x).

Hence the algorithm proceeds as follows: knowing all ̺ n and u j+ 1 2 ,n (j = q, ..., Nq -1) at some time t n = n∆t, we first update u j+ 1 2 ,n+1 inside B with (4.3) in order to deduce ̺ j,n+1 by means of (4.2). Then, knowing ̺ q+1,n+1 and ̺ N -q-1,n+1 , we can invert the implicit scheme (4.1) as soon as we have chosen a correct boundary condition to render the junction between the diffusive/ballistic zones.

A (nonlinear) Robin boundary condition.

Let us mention first that junctions have been put at staggered abscissae, x q+ 1 2 , x N -q-1-1 2 , which correspond to the velocity unknowns. At these points, imposing Dirichlet boundary conditions leads to some kind of inconsistency because of a somewhat redundant definition of the current J between e.g. (CL) and (B). Numerically, it produced some spurious oscillations eventually leading to general instability. Moreover, in order to observe a constant current everywhere inside the device at steady state, it sounds reasonable to impose its continuity at the junctions between diffusive and ballistic areas. That is to say, at x q+ 1 2 ,

u q+ 1 2 ,n+1 = -d ̺ q+1,n+1 -̺ q,n+1 ∆x̺ q+ 1 2 ,n+1 - U q+1,n+1 -U q,n+1
∆x , which implies:

-d ̺ q+1,n+1 -̺ q,n+1 ∆x = ̺ q+ 1 2 ,n+1 u q+ 1 2 ,n+1 + U q+1,n+1 -U q,n+1 ∆x ,
that we can approximate by means of

-d ̺ q+1,n+1 -̺ q,n+1 ∆x = θ̺ q+1,n+1 +(1-θ)̺ q,n+1 u q+ 1 2 ,n+1 + U q+1,n+1 -U q,n+1 ∆x , (4.4)
for θ ∈ [0, 1]. Hence one recognizes on the left-hand side of (4.4) the unknown flux which appears in (4.1) and the M q matrix, and the right-hand side contains only values which can be deduced from the outcome of the explicit schemes used inside B. Actually, only the average value ̺ q+ 1 2 ,n+1 has to be arbitrarily chosen as e.g. the arithmetic average 1 2 (̺ q+1,n+1 + ̺ q,n+1 ), ̺ q,n+1 or ̺ q+1,n+1 , corresponding to θ = 1 2 , θ = 0 and θ = 1 respectively.

Remark 3. This choice of boundary conditions corresponds to the solving of the following problem for the drift-diffusion equation (CL):

∂ t ̺ + ∂ x J = 0, J = -d∂ x ̺ -̺∂ x U, x ∈]0, x L [, ̺(t, x = 0) = D(x = 0), -∂ x ln(̺)(t, x = x L ) = u + ∂ x U
where u(t, x = x L ) is deduced from (B). This is a mixed (nonlinear) Robin-type boundary condition for a strictly parabolic equation. The well-posedness theory for problems of this flavour is quite recent; it started with the paper by L. C. Evans [START_REF] Evans | Regularity properties for the heat equation subject to nonlinear boundary constraints[END_REF] and then developed further into [START_REF] Arrieta | Parabolic problems with nonlinear boundary conditions and critical nonlinearities[END_REF][START_REF] Favini | Nonlinear boundary conditions for nonlinear second order operators on C[0,1[END_REF][START_REF] Goldstein | Smoothing for nonlinear parabolic equations with nonlinear boundary conditions[END_REF].

Of course, we presented everything in x = x q+ 1 2 , but the right junction x = x N -q-1-1 2 is completely equivalent (with maybe a different choice of θ however). 4.4. Some stability results. First of all, we wish to state a result from [START_REF] Goldstein | Smoothing for nonlinear parabolic equations with nonlinear boundary conditions[END_REF] which we adapt to our situation involving only a uniformly elliptic linear operator: Theorem 4.1. (Goldstein & Goldstein, [START_REF] Goldstein | Smoothing for nonlinear parabolic equations with nonlinear boundary conditions[END_REF]) Let β j (ξ) be a strictly increasing and continuous function. The problem:

   ∂ t u = ∂ xx u, x, t ∈ (0, 1) × R + , (-1) j ∂ x u(., x) = β j (u), x = j ∈ {0, 1}, u(t = 0, x) = f (x) ∈ L 1 (0, 1)
x ∈ (0, 1), has a unique mild solution u in (0, 1) × [0, T ] for any T > 0 given. Moreover, u(t, .) ∈ L ∞ (0, 1) and ∂ t u(t, .) ∈ L 2 (0, 1).

Clearly in order to be complete, a full study of the coupled problem (even if neglecting the Poisson coupling equation) would be necessary before concluding to the stability of the whole model given by (CL)-(B)-(CR). However, the former theorem is a stepping stone in this direction. An interesting feature of our numerical approach (4.1)-(4.2) lies in a nonnegativity-preserving property for the density.

Theorem 4.2. Let θ = 0 in (4.4), max j,n (|u j+ 1 2 ,n |+|U j+1,n+1 -U j,n+1 |/∆x)∆t ≤ ∆x and ̺ n ∈ (R + ) N , then the whole scheme (4.1)-(4.2) preserves nonnegativity in the sense that

̺ n+1 ∈ (R + ) N .
Actually, we never used (4.1)-(4.2) with θ = 0 because it looks like being the less accurate choice. However we were unable to prove Theorem 4.2 for θ > 0.

Proof. The proof proceeds in 3 steps, namely the explicit part, then the implicit one and what happens at the junctions.

• We first observe that thanks to the prescribed CFL condition, the explicit parts of both schemes (4.1) and (4.2) turn out to be convex combinations of the values showing up on their right-hand sides. This has been precisely shown in previous paper, see e.g. Proposition 4.6 in [START_REF] Gosse | Numerical approximation of linear one-dimensional conservation equations with discontinuous coefficients[END_REF] or §4.3 in [START_REF] Gosse | Convergence results for an inhomogeneous system arising in various high frequency approximations[END_REF]. • The implicit part of (4.1) is given by M q which is a M-matrix in the terminology of [START_REF] Berman | Nonnegative matrices in the mathematical sciences[END_REF], that is to say its entries satisfy:

m i,i > 0, m i,j =i ≤ 0, |m i,i | > j =i |m i,j |.
So let's suppose we solve M q x = b where b has nonnegative components, it is a classical fact that x will have only nonnegative components too. For completeness we recall its (quick) proof: let i 0 corresponds to the smallest component of x, x i0 := min( x). There holds:

m i0,i0 x i0 = b i0 - j =i0 m i0,j x j ≥ b i0 -x i0 j =i0 m i0,j ,
which implies:

x i0 ≥ b i0 m i0,i0 + j =i0 m i0,j ≥ 0.
• It remains to observe that choosing θ = 0 in (4.4) means that the boundary condition shows up only in the implicit part of (4.1). Let us focus on the left junction: inside M q , we have α = 2 and β = 1. The condition (4.4) gives us another term which reads u q+ 1 2 ,n+1 + Uq+1,n+1-Uq,n+1 ∆x , but thanks to the CFL condition, we know that M q remains diagonal-dominant even on this line. Hence the preceding computation garantees the nonnegativitypreserving property and we are done.

These computations are based on convexity arguments hence they somehow open the way to compactness for the whole scheme. However, we do not have any precise theory at the continuous level for (CL)-(B)-(CR) yet. This might constitute a direction for further investigation.

Remark 4. Actually the choice of the parameter θ is quite a delicate task when considering its importance in order to make the numerical scheme stabilize on a steady-state regime with flat currents (as can be seen on the numerical figures). Here is what we implemented in order to produce the numerical results of §5: first, we distinguished between the left junction where the solution is Lipschitz-continuous and the right one where it generally involves a discontinuity. Then we observed that it wasn't possible to code directly ∂ x ln(̺) = u + ∂ x U simply because it would require to apply the change of variable v = ln(̺) while solving the drift-diffusion equation, which isn't easy. Hence it becomes necessary to linearize the ln(.) function and this is where the θ ∈ [0, 1] parameter comes into play through the mid-point rule. In case this linearization occurs in a smooth location, θ = 1 2 gives a second-order accuracy and is therefore sufficient. In case a discontinuity shows up, the upwind value is often the right choice, that is to say, the weight should go to the point from where the signal comes from (generally, it is the left side).

5. Numerical validation with doping 5 × 10 22 m -3 . In this section, we display some numerical results at steady-state for a standard GaAs device as already considered in [START_REF] Carrillo | Computational macroscopic approximations to the 1-D relaxation time kinetic system for semiconductors[END_REF][START_REF] Arnold | Low and high field scaling limits for the Vlasov and Wigner-Poisson-Fokker-Planck systems[END_REF][START_REF] Blotekjaer | Transport equations for electrons in two-valley semiconductors[END_REF][START_REF] Shur | Near ballistic transport in GaAs devices at 77K[END_REF]. We took a discontinuous doping profile equal to 1 inside x ∈ [0, 0.25[∪]0.75, 1] and 0.002 elsewhere. We iterated the marching schemes presented in §3.2 until residues on ̺ n decay lower than 10 -15 . The resulting parameters are d 1 = 0.00118 and eU 0 /(m * v 2

2 ) = 4. We stress that neither "drift-diffusion nor mixed currents" are flat in the whole computational domain at steady-state; so in order to derive the current-voltage relations, we computed arithmetic averages of the results coming from both discretizations; we shall however propose a more interesting derivation in §5.5 relying on Well-Balanced ideas. 255 points have been used to sample the [0, 1] interval and the time-step is chosen adaptively in order to minimize the mixed scheme's numerical diffusion. The main difference between both curves in Fig. 5.1 lies in slightly lower values for the drift-diffusion currents. Apart from that, both approaches appear to agree. However, we are about to see that despite the relations current-voltage may look similar, the densities and velocities are rather different.

Remark 5. It may sound reasonable to think about simulating directly the equation (2.2), at least to validate the outcome of our mixed scheme as in 1-D, the computational cost can be still considered reasonable. Indeed, performing the rescalings as in §3 leads to the following problem: supplemented by (3.5) and v being a free parameter. Actually, it is also possible to evaluate all the parameters involved: ǫ ≃ 0.09, η ≃ 0.3 (tending to zero in §3.1) and β ≃ 0.1, α ≃ 2 (being o(1) and O(1) respectively in §3.2). However, there's at least one obstacle on this way: the problem still involves two different time scales, and two different Maxwellian distributions. Hence a direct upwind discretization of (3.11) is unstable because some boundary conditions are needed at each junction: a way out may consist in implementing numerically the formulas ( 16)-( 17) in [START_REF] Benabdallah | Coupling one-dimensional time-dependent classical and quantum transport models[END_REF] (see also [START_REF] Degond | A smooth transition model between kinetic and hydrodynamic equations[END_REF] for related questions). The schemes should also be asymptotic-preserving in the sense of [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF]. Finally, we notice from [START_REF] Cercignani | Device benchmark comparisons via kinetic, hydrodynamic and high-field models[END_REF] that direct kinetic simulations deliver results quite similar to drift-diffusion ones.

v∂ t f + v th ξ∂ x f + v th U 0 U th ∂ x U ∂ ξ f = v relax (M(θ 0 )̺ -f ),

Weak bias:

[Φ] = -0.2. Essentially, the main difference between the mixed approach and the global drift-diffusion solving is that the last one always produces continuous solutions whereas jumps show up in the other ones. This is clearly demonstrated by the results in Fig. 5.2, even if the corresponding diffusion coefficient turns out to be relatively small. Velocities and densities turn out to be rather different for both approaches, especially when compared to the Poisson potential and the resulting electric field. We can observe that our prescribed boundary condition for the junction works fine for continuous solutions, but is unable to produce a flat current close to a shock. This can come from the fact it involves an average value, ̺ q+ 1 2 ,n+1 which should perhaps be derived more precisely. However, outside the junctions, currents look satisfying.

Moderate bias:

[Φ] = -0.35. In this case, the higher bias accelerates electrons; the outcome of both schemes is presented on Fig. 5.3. The residues decay is satisfying and the steady-state currents are reasonably flat. Densities and velocities are again very different for both approaches, quantitatively and qualitatively.

Higher bias:

[Φ] = -0.5. In this last case, the higher bias tends to flatten more the "mixed current", see Fig. 5.4 despite some sort of boundary layer on the right side. The potentials and the electric fields are quite similar for both approaches even if velocities and densities are still very different. Convergence in time is obtained in half the time it asked for with [Φ] = -0.2. 5.4. Electric field and velocities inside the device. For completeness, we finally plot the fraction u(x)/E(x) at steady-state since this quantity measures the conductivity as a function of the applied bias; see Fig. 5.5. We observe that independently of the bias, this fraction is practically constant in the zones where the drift-diffusion algorithm is used (this reflects the smallness of the diffusion coefficient for high doping). However, it falls down strongly inside the channel even if increasing the bias improves things. Junctions are slightly visible, but the corresponding discontinuities don't seem to have strong effects on the global picture.

5.5.

A different way to compute currents in the device. We finally propose a new way to compute steady-state currents inside the device, owing to the fact that a Well-Balanced discretization is used to solve the ballistic part (B) of the problem; see Appendix A for details. Actually, assuming (B) is stationary, one derives easily that the following equality holds:

J bal (x) = ̺(x)u(x) = ̺(x) -2(U (x) -U (x L )) + u(x L ) 2 . (5.1)
The same way, one actually can compute stationary drift-diffusion like: formulas (5.1)-(5.2). One clearly sees the inside the ballistic channel, the resulting currents are perfectly flat, which has never been obtained before up to the authors' knowledge. Even the drift-diffusion currents are rendered in a slightly better way: for [Φ] = -0.35, there's no jump around x = x R . Nevertheless, this isn't meant to improve the treatment of the junctions where small discontinuities still appear. We believe that only a discretization of the type (6.1) can solve this issue.

J DD = -̺(d∂ x ln(̺) + ∂ x U ). ( 5

Conclusion and outlook.

We presented in this paper a first attempt to develop a robust mixed approach to transient semiconductor's computations relying on a simple time-marching strategy. The main obstruction to circumvent was the choice of a time-step which would ensure the compatibility of both the parabolic (drift-diffusion) and hyperbolic (ballistic) areas, left apart the derivation of a stable boundary condition at the junctions between the source/drain and the channel. For completeness, we display here the current-voltage relation for a doping value of 10 24 m -3 , see Fig. 6.1. With a diffusion coefficient so low, both curves look very similar. There exist several ways of extending such an approach: one could be the inclusion of the energy bands, more accurate than the simple "effective mass" approximation made in (1.1). Indeed, such a simplification makes the whole model very close to the so-called "Drude theory of conduction" (see [START_REF] Ashcroft | Solid-state physics[END_REF]), which is known to be quite crude. In order to improve the overall accuracy, an efficient algorithm to compute the Bloch decomposition in 1-D has already been studied in [START_REF] Gosse | Multiphase semiclassical approximation of an electron in a onedimensional crystalline lattice -I. Homogeneous problems[END_REF]. It may also be interesting to use the sophisticate "multi-branch" routines of [START_REF] Gosse | Using K-branch entropy solutions for multivalued geometric optics computations[END_REF][START_REF] Gosse | Multiphase semiclassical approximation of an electron in a onedimensional crystalline lattice -I. Homogeneous problems[END_REF] in this present context to achieve more precise simulations inside the ballistic channel. A routine able to handle phonons has been recently proposed in [START_REF] Gosse | Multiphase semiclassical approximation of the one-dimensional harmonic crystal. I. The periodic case[END_REF].

From the perspective of the "well-balanced" aspects of things, the big step would be to succeed in working with a good set of variables which would prevent from linearizing at the junctions (recall the delicate choice of the θ parameter, see Remarks 3 and 4). The first idea could be to replace the ̺ variable in (CL) and (CR) by z = ln ̺. However, the price to pay is to make these drift diffusion problems nonlinear; indeed the equations become in the z variable

∂ t z + d ∂ xx z + |∂ x z| 2 + ∂ xx U + ∂ x U ∂ x z = 0, exp(z)(•, x = 0) = 1, ∂ x z(•, x = x L ) = u + ∂ x U, (6.1)
the left junction being rendered through the last boundary condition. However, the nonlinear term is first-order, hence asks only for a mild CFL restriction.

Remark 2 .

 2 (Numerical values) Let us take the same test case as in[START_REF] Baranger | Ballistic structure in the electron distribution function of small semiconducting devices: general features and specific trends[END_REF][START_REF] Carrillo | Computational macroscopic approximations to the 1-D relaxation time kinetic system for semiconductors[END_REF]: that is to say, m * = 0.065 × 9.109 × 10 -31 Kg, e = 1.602 × 10 -19 C, k B = 1.38 × 10 -23 J/K, ε s = 13.2 × 8.85418 × 10 -12 F/m, T 0 = 300K, L = 0.8 × 10 -6 m. Moreover in the classical zone, we take for the doping concentration D = 10 24 /m 3 and for the electron mobility µ 1 = τ 1 e/m * = 0.75. In the ballistic zone, D = 2 × 10 21 /m 3 and µ 2 = τ 2 e/m * = 4. These values imply that ̺ 0 = 5 × 10 23 /m 3 , U 0 = e̺oL 2 εs = 438.6. Concerning the velocity v 1 = µ1U0 L = 4.111 × 10 8 m/s which implies v 2 = v 1 µ 1 /µ 2 = 7.71 × 10 7 m/s. The scaled diffusion coefficient is d 1 = U th U0 = 5.9 × 10 -5 (and the scaled drift coefficient is µ = 1). Finally we have eU 0 /(m * v 2

Fig. 5 . 1 .

 51 Fig. 5.1. Current-Voltage relation with doping 5 × 10 22 m -3 (left, normal scale; right, log scale).

Fig. 5 . 2 .

 52 Fig. 5.2. Carriers density, velocity, potential, electric field, currents and residues for an applied bias of [Φ] = -0.2. (left to right, top to bottom).

Fig. 5 . 3 .

 53 Fig. 5.3. Carriers density, velocity, potential, electric field, current and residues for an applied bias of [Φ] = -0.35. (left to right, top to bottom).

Fig. 5 . 4 .

 54 Fig. 5.4. Carriers density, velocity, potential, electric field, current and residues for an applied bias of [Φ] = -0.5. (left to right, top to bottom).

Fig. 5 . 5 .

 55 Fig. 5.5. Steady state fraction u(x)/E(x) for applied biases [Φ] = -0.2, -0.35, -0.5. (left to right)

Fig. 5 . 6 .

 56 Fig. 5.6. Steady state currents for the mixed scheme with [Φ] = -0.2, -0.35, -0.5. (top to bottom)

Fig. 6 . 1 .

 61 Fig. 6.1. Current-Voltage relation with strong doping (left, normal scale; right, log scale).

it corresponds to an energy level for which the occupation probability is exactly one half.

for instance, in case Q is given by a BGK model.

if not, breakup occurs and multibranch solutions are to be sought instead, like in[START_REF] Gosse | Using K-branch entropy solutions for multivalued geometric optics computations[END_REF][START_REF] Gosse | Multiphase semiclassical approximation of an electron in a onedimensional crystalline lattice -I. Homogeneous problems[END_REF][START_REF] Shi Jin | Booske An Eulerian method for computing multivalued solutions of the Euler-Poisson equations and applications to wave breaking in klystrons[END_REF].

Acknowledgments. The second author warmly greets the Mathematics Department of University of Granada for its hospitality (and Prof. J. Soler in particular). Both authors were partially supported by the EEC project HYKE (Hyperbolic and kinetic equations) #HPRN-CT-2002-00282.

Appendix A. Well-balanced for (B) in a nutshell. In this short section, we follow [START_REF] Gosse | A well-balanced flux splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF] in order to explain quickly how to manage the forced Burgers equation appearing in (B) in such a way to approximate steady-states as best as possible. Let us suppose that u(t, .) > 0 holds (this is completely fair looking at the kind of steady regimes we have in mind to reproduce), we can rewrite the scalar balance law the following way:

This clearly shows that steady states with u = 0 are given by the conservation law

At this point, it is easy to include the Poisson potential U in the upwind fluxes in order to maintain this property: thus we obtain the modified fluxes,

with the upwinded values, ũl = u 2 l -2(U r -U l ) ≥ 0, ũr = -u 2 r + 2(U r -U l ) ≤ 0.


Maxime Montaru
email: maxime.montaru@inrets.fr

Serge Pelissier
email: serge.pelissier@inrets.fr

Frequency and Temporal Identification of a Li-ion Polymer Battery Model Using Fractional Impedance

Keywords: Lithium-ion Polymer, modelling, diffusive impedance, identification, spectroscopy impedance, chronopotentiometry, non linear model

Frequency and Temporal Identification of a Li-Polymer Battery Model Using

Fractional Impedance -The modelling of high power Lithium batteries may work as a powerful tool in the sizing of batteries pack in hybrid vehicles. Moreover, this tool can be used to evaluate accurately the battery characteristics at each level of its lifetime. Until now, the model structure differs according to its end-use due to simulation limitation. Indeed, impedance measurements obtained by electrochemical impedance spectroscopy needs non-integer order impedance in order to correctly model diffusive phenomenon. However, this impedance cannot be simulated easily in software dedicated to hybrid vehicles and needs approximations. Furthermore, impedance measurements are realised with low currents comparing to currents that are forecast in normal uses of high power batteries and can induce errors if the impedance is varying with the current. This paper deals with non linear dynamic models that use band-limited frequency impedance with non-integer order. The parameters are identified partly by impedance measurements for high frequencies and partly by chronopotentiometry measurements for low frequencies. Validation tests are made with staircase of high intensity and current profile simulated on a hybrid vehicle software. The hypothesis of non-linearity is verified.

INTRODUCTION

The battery is one of the essential parts in a Hybrid Electric Vehicle (HEV) and probably the most critical one with respect to reliability and life expectation. It is a primordial element in the concept of hybridation because it enables to reversibly store energy. The energy stored is picked up during deceleration of the vehicle and used for starting the thermal engine and giving a boost during acceleration. The battery can be used in a sustaining mode, i.e. in a small window of State Of Charge (SOC), or in a depleting mode, i.e. in a large window of SOC. The second case enables a succession of high power discharges which permits the vehicle to reduce much more its consumption or even to stop completely the thermal engine. The power supplied and received by the battery is variable and could be very important depending notably on the sizing of the pack of batteries and its management. In order to optimize its sizing and to maintain a correct lifetime, the knowledge of the model of the electric storage components and their dynamic response is necessary.

The modelling of batteries is not an easy task. The reasons are [START_REF] Buller | Impedance-based non-linear dynamic battery modelling for automotive applications[END_REF]: (i) batteries are not stationary due to the SOC variation and diffusion phenomenon, (ii) they are systems highly non-linear regarding to the current, temperature, etc... Two families of identification of the model parameters, temporal and frequency, are currently used but not necessarily for the same purpose.

In one hand, it is possible to use the temporal chronopotentiometry. This technique consists of applying a pulse of current and measuring the voltage of the elements. It's the method used by the Patnership of the New Generation of Vehicle [START_REF]Freedom CAR Battery Test Manual For Power-Assist Hybrid Electric Vehicles[END_REF] who names it High Power Pulse Characterization. The aim is to evaluate the resistance of batteries after a certain time of applied current. It can also be used to identify much more complex model as it was proposed in the study with ADVISOR models [START_REF] Johnson | Battery performance models in ADVISOR[END_REF]. It's also possible to use profile of current to model dynamic behaviour. This profile can be a pseudo random binary sequence [START_REF] Sabatier | Fractional system identification for lead acid battery state of charge estimation[END_REF] or simply represent a profile of current which could be measured on a HEV [START_REF] Plett | Extended Kalman filtering for battery management systems of LiPB-based HEV battery packs, Part 2. Modeling and identification[END_REF]. In the other hand, a well-known technique for impedance characterization is impedance spectroscopy [START_REF] Barsoukov | Impedance Spectroscopy, Second Edition[END_REF]. This technique can be used to identify relationship between material properties and battery performance if just one interface is studied [START_REF] Barsoukov | Parametric analysis using impedance spectroscopy: relationship between material properties and battery performance[END_REF][START_REF] Itagaki | In situ electrochemical impedance spectroscopy to investigate negative electrode of lithium-ion rechargeable batteries[END_REF]. It is risky to do the same with a complete cell because it's not possible to dissociate the effect of the cathode from the one of the anode. However, it is possible to associate the impedance measured by spectroscopy to the state of charge, or even state of health ( [START_REF] Bundy | An electrochemical impedance spectroscopy method for prediction of the state of charge of a nickel-metal hydride battery at open circuit and during discharge[END_REF][START_REF] Blanke | Impedance measurements on leadacid batteries for state-of-charge, state-of health and cranking capability prognosis in electric and hybrid electric vehicles[END_REF]). Moreover, the impedance of a battery cell has been identified with this technique [START_REF] Kuhn | Comparison between two NiMH battery models[END_REF]. Currently, in electrochemical studies, the voltage oscillation is imposed around the Open Circuit Voltage (OCV) but in case of commercial cell, it is possible to impose current oscillation and to control state of charge during the operations. The measurement of impedance can be done with or without current polarization. When the spectrometer can supply sufficient continuous current, it is possible to highlight non linearity of the battery impedance with the current level [START_REF] Buller | Impedance-based non-linear dynamic battery modelling for automotive applications[END_REF][START_REF] Bundy | An electrochemical impedance spectroscopy method for prediction of the state of charge of a nickel-metal hydride battery at open circuit and during discharge[END_REF][START_REF] Blanke | Impedance measurements on leadacid batteries for state-of-charge, state-of health and cranking capability prognosis in electric and hybrid electric vehicles[END_REF][START_REF] Karden | A frequency-domain approach to dynamical modelling of electrochemical power sources[END_REF][START_REF] Thele | Development of a voltage-behaviour model for NiMH batteries using an impedance based modelling concept[END_REF][START_REF] Mauracher | Dynamic modelling of lead/acid batteries using impedance spectroscopy for parameter identification[END_REF]. Another frequency method to identify battery impedance consists of using Laplace transform of temporal solicitations [START_REF] Takano | Simulation study of electrical dynamic characteristics of lithiumion battery[END_REF][START_REF] Tenno | A Method for Battery Impedance Analysis[END_REF].

The advantage of chrono-potentiometry is to allow discharge and charge of high current level as in normal use. However, the acquisition frequency may be limited by the equipment and prevents from identifying correctly the impedance at high frequencies. On the contrary, impedance spectroscopy gives a precise impedance measurement in a wide band of frequencies but only for low current values. In order to conjugate the advantage of each technique, the method used in this study consists of a partial identification of the model by spectroscopy impedance and the rest by chronopotentiometry. It enables a gain of time due to the less time spent in impedance spectroscopy at low frequencies and a more complete model with the nonlinearity of impedance with current. Furthermore the equipment is standard unlike some other works [START_REF] Thele | Development of a voltage-behaviour model for NiMH batteries using an impedance based modelling concept[END_REF][START_REF] Mauracher | Dynamic modelling of lead/acid batteries using impedance spectroscopy for parameter identification[END_REF]. The model used in this study is an equivalent electric circuit with lumped parameters. In such a model, parameters could be tabulated according to SOC, current, temperature, etc… [START_REF] Thele | Development of a voltage-behaviour model for NiMH batteries using an impedance based modelling concept[END_REF][START_REF] Kuhn | Modelling Ni-MH battery using Cauer and Foster structures[END_REF]. A representation of our model depending on SOC and current is given in Figure 1. It consists of an open circuit voltage OCV in series with a complex impedance. The tabulation of the OCV according to SOC could be a very difficult task if hysteresis phenomenon occurs as observed with NiMH batteries [START_REF] Thele | Development of a voltage-behaviour model for NiMH batteries using an impedance based modelling concept[END_REF], but it is not the case in our study since hysteresis is very small in this type of Li-ion battery. The complex impedance is composed of 3 elements in series: a resistance R1 for modelling the ohmic behaviour, a resistance R 2 in parallel with a Constant Phase Element CPE 2 for modelling the charge transfer and double layer processes and another Constant Phase Element CPE 3 associated to the diffusion behaviour. This impedance model is adapted to fit well with impedance spectrum (cf. Figure 1). In many studies, the diffusive impedance is represented by a Warburg impedance which expression, in frequency domain, is [START_REF] Mauracher | Dynamic modelling of lead/acid batteries using impedance spectroscopy for parameter identification[END_REF][START_REF] Kuhn | Modelling Ni-MH battery using Cauer and Foster structures[END_REF]. However, in our model, we used a simplified expression of the Warburg impedance obtained with a constant phase element CPE 3 . The difference between these two expressions occurs at very low frequencies not reached in this study.

) tanh( 1 ) ( τω ω ω j j Q j Z = [

Modelling the Constant Phase Element CPE 3

A CPE consists of non integer order impedance which is expressed in the frequency domain as:

α ω ω ) .( 1 ) ( j Q j CPE = (1)
With ω the pulsation in rad.s -1 , α a real number comprised between 0 and 1 and Q a constant expressed in . In the case of α = 1, CPE is equivalent to a capacitance. In the case of α = 0, CPE is equivalent to a resistance. A CPE is represented as a straight line in Nyquist plan (Fig. 1) and in Bode diagram (Fig. 2). Simulation in temporal mode necessitates finding an equivalent circuit of CPE. Several authors proposed to use a mathematical series of R//C circuits as an equivalent of the diffusive Warburg impedance [START_REF] Mauracher | Dynamic modelling of lead/acid batteries using impedance spectroscopy for parameter identification[END_REF][START_REF] Kuhn | Modelling Ni-MH battery using Cauer and Foster structures[END_REF]. Based on the same principle, a CPE could be approximated by a fractional impedance in a limited frequency band [START_REF] Sabatier | Fractional system identification for lead acid battery state of charge estimation[END_REF][START_REF] Oustaloup | Représentation et identification par modèle non entier[END_REF]. It consists of several R//C circuits which frequency characteristics are judiciously chosen in order to approximate a straight line only in a limited frequency band [ω min , ω max ] as shown in Figure 2. In this work, we used five R//C as a good compromise between accuracy and duration of computation. The expression of the impedance is given by:

1 . - α s F 5 4 2 1 1 1 ' 1 . .... . 1 ' 1 . 1 ) ( ω ω ω ω ω ω ω ω ω ω γ ω j j j j j j Z + + + + + = (2)
The pulsations ω i and ω' i of each R//C circuits are depending on the others by these relations:

i i ω ω β log log log 1 1 - = + (3) i i ω ω β ′ - = + log log log 1 2 (4)
The slope of the approximated impedance in Bode diagram is given by: The impedance of the circuit R 2 //CPE 2 , also named ZARC [START_REF] Barsoukov | Impedance Spectroscopy, Second Edition[END_REF], is expressed by:

2 ) ( 1 ) ( 2 2 2 α ω τ ω j R j Z + = with 2 2 2 Q R = τ (6)
It is possible to extract the value of these parameters from an impedance spectrum (cf. Figure 3). Firstly, the intersection of the impedance spectrum with the abscissa line enables to detect R 1 . Secondly, the spot at top of the semi-circle is characterized by its real part, its imaginary part and its pulsation ( ) ( )

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ℑ 2 . cos 1 2 2 . sin 2 2 2 2 , 2 π α π α ω R Z c m (8)
Then, the parameter 2 τ , and consequently , is deduced from the eq.( 6) of the pulsation of the top of the semi-circle. The phenomenon of non-linearity is taken into account for the resistance R 2 . The formalism of Butler-Volmer is considered for this purpose [START_REF] Thele | Development of a voltage-behaviour model for NiMH batteries using an impedance based modelling concept[END_REF]:

) ( ) 1 ( 2 η η B a aB A e e i - - - = with η = u -OCV. (9)
A hypothesis of symmetry is made concerning the nonlinearity between the charge and the discharge (a = 0.5). So, ) sinh( η

B A i = ( 1 0 )
It is possible to deduce an expression of the equivalent resistance:

) cosh( 1 2 η η B AB d di R = = (11)
Then, an expression of R 2 (i) can be derived:

C A i AB i R + + = 2 2 ) ( 1 1 ) ( (12)
For the need of fitting, it appears that it is necessary to add a constant C to R 2 . The representation of this function is given at Figure 4.

Simplified model

It appears when identifications of these parameters are made that τ 2 is comprised between 20 and 30 ms. Considering that this time constant is very small compared to acquisition period of the solicitation profile (100 ms), the model Z 2 can be simplified. A pure resistance R 2 replaces the R 2 //CPE 2 which results on a simplified model (cf. Figure 5) to be used for simulation.

R1(SOC) i(t) u(t) OCV(SOC) CPE3(SOC) R2(SOC,i)

Figure 5. Simplified model of li-ion battery used in temporal simulation

To conclude, at each SOC, the battery model used in temporal simulation is defined by: -the OCV, -the resistance R 1 , -the parameters of the non linear model of R 2 : A, B,C, -the parameters of the diffusive impedance CPE 3 : ω 1 , ω 5 γ, p.

So, 9 parameters are needed to be identified in order to simulate correctly the model. Moreover, the evaluation of each of them could give great information on battery behaviour according to state of charge, temperature and state of health.

METHODS & EXPERIMENTAL SETUP

The experimental data of this study concern a Li-ion Polymer battery (KOKAM SLPB 702051300P -12Ah) tested at 20°C.

Existing methods for identification of the model of battery

As previously said, the two main possibilities to obtain experimental data for identification of the model of a battery are impedance spectroscopy and chronopotentiometry.

•Impedance spectroscopy is based on the measurement of complex impedance of the battery at various frequencies . The advantage of this method is the full signature obtained for the battery. Difficulties exist to measure impedance at very low frequencies (below 1mHz) because of the duration of the operation, which could last several hours in this case, and the eventual consequently no-stability of the result. Another limitation of this method is the difficulty for measurement of impedance when a DC current is superposed on the AC measurement signal in order to characterize the impedance of the battery during charge or discharge phases. With high values of current, the SOC of the battery could be largely modified during the measurements; in this case the impedance is not constant all along the test and compensation drift procedures must be applied [START_REF] Buller | Impedance-based non-linear dynamic battery modelling for automotive applications[END_REF][START_REF] Jp | EIS study of electrochemical battery discharge on constant load[END_REF].

) ( ) ( ) ( f i f f Z m e ℑ + ℜ =
•Chronopotentiometry is based on the voltage response of the battery when a current profile is imposed. The impedance is obtained with the u(t)/i(t) ratio. By using different levels of signals, it could be easy to determine the impedance of the battery for any value of charge or discharge currents. The main limitation of this method is the need for a high sampling rate if one wants to measure the impedance at high frequencies. In practice, two sorts of equipments can be used: either a battery tester with high power capabilities but usually with a limited synchronous sampling period (for example 100 ms) or an impedance analyser with high frequencies capability but usually with a limited current (for example 20A). A compromise must be found between frequency and power, with a large incidence on the price of the equipment.

Proposed method for identification of power batteries

We propose a two step identification process based on two frequency domains. The identification of R 1 , R 2 and CPE 2 is performed at "high frequencies", in the 1kHz -500mHz band, using an impedance analyser. The identification of the elements of the diffusive impedance is obtained by chronopotentiometry with low frequency pulses supplied by a battery test equipment.

• First step: characterization of R 1 , R 2 and CPE 2 .

The aim is to measure the values of these elements for different values of the SOC of the battery and different values of the charging/discharging current. This current is imposed during a pulse of 20s. An impedance spectroscopy in the 1kHz -500mHz band with 12 s duration, is performed 5 s after the beginning of the pulse. A charging pulse is followed by a discharging pulse to keep SOC constant. This couple of pulse is repeated with several levels of current (0, 6, 12 and 18 A). Measurements are made at different value of SOC (80, 60, 40 and 20 %). The duration of the pulse is shorter enough to consider no variation in the value of SOC (18 Ax20 s = 0.1 Ah is less than 1% of the capacity of the battery). The equipment used for this measurement is a VSP Biologic with a 20 A amplifier kit. An example of the Nyquist graph of impedance for one SOC value and one level of current is shown on the Figure 6. The values of R 1 , R 2 , and α 2 are determined with the methods explained in the previous part. As shown on this figure, one needs to interpolate the point which intercepts the abscissa and the one which has the maximum imaginary part. The Figure 7 gives the spectroscopy results for different levels of to identify the parameters of the theoretical relation R 2 (i) (cf. Eq. 12). As the maximum value of the experimental current is 18 A, it is necessary to obtained at least one data point for higher current in order to have better fitting. This could be achieved by a complementary measurement of chronopotentiometric response to a high intensity pulse. We used a battery tester composed of a DC source (Delta SM30200 30 V -200 A) associated with an electronic load (H&H ZS4206 -60 V -450 A) to generate short pulses with an amplitude between 50 A and 200 A in discharge and between 50 A and 100 A in charge. The charging pulse duration is 2 s. One must calculate the duration of the discharging pulse to balance the SOC. An example of the battery voltage during pulse is given by the Figure 8. We propose to estimate the value of the resistance from the voltage drop. The sampling frequency (10 Hz) is not high enough to have numerous data points during the front voltage and a linear extrapolation must be done. The aim is to measure the values of these parameters for different SOC and with pulse profile capable of giving information in the 1 Hz -1 mHz band. We have chosen to use two 5 min -1C charging and discharging pulses preceded and followed by a 1 hour rest. All these operations are made at different values of SOC (80, 60, 40 and 20 %). The equipment used for this measurement is the same battery tester previously cited. A preliminary test is necessary to determine the evolution of the OCV of the battery versus its SOC. The curve is given by Figure 9. A small hysteresis of 20 mV is observed. We decided to neglect this phenomenon in this work. An average value of this curve is used to compensate the variation of the OCV during a 5 min pulse. An example of the measured and compensated signals is shown on Figure 10. A Matlab-Simulink program optimizes the diffusive impedance parameters to fit the simulated response with the experimental data by minimising a least square criterion. In order to not limit the influence on the optimization of the long time response, the experimental data used for this optimization is selected at instants with a logarithmic distribution after each change in the current. The selected experimental data instants are also shown in Figure 10.

Validation of the model 3 RESULTS

Identification of R 1 and R 2 //CPE 2

In order to validate the model of the battery, the voltage response of the model to a "validation profile" is compared to the real response of the battery. This "validation profile" is composed of a suite of charging/discharging pulses defined from experimental measurement on a real hybrid vehicle adapted to the capabilities of the tested battery.

As explained above, the resistance R 1 and the parameters of R 2 //CPE 2 are determined from impedancemetry data for different SOC. The results are represented in Figure 11 where each parameter is represented for different levels of current. The values obtained at different SOC are superposed on the same graph. We can see that R 1 is relatively constant and will be set to 1,95 mΩ. Concerning the parameters of R 2 //CPE 2 , we can say that R 2 and depend on the current level and SOC. R 2 tends to decrease when the current increases. As tends to have opposite variation, the time constant τ 2 is expected to be rather constant, which is effectively observed. The limited variation of α 2 can be supposed of low impact on the value of the impedance. For these reasons, we will keep α 2 and τ 2 equal to a value which will depend only on SOC. The behaviour of R 2 as a function of the current level is modelling by eq.( 12). The data obtained for different SOC can be fitted with the same parameter value of A and B. There is only C which depends on SOC. The fitting are represented with experimental data for different SOC in Figure 12. The results obtained with impedance spectroscopy at low currents and the data obtained from the high current pulses are in a good agreement (Cf. Figure 12).

Identification of CPE 3

Now, the parameters of the diffusive impedance CPE 3 (ω 1 , ω 5 , γ and p) have to be evaluated from temporal data. In order to assure convergent results from optimization,

we fixed the values of ω The result of the optimization can be seen on the voltage chronograms represented for the different SOC in Figure 14. The good agreement is attested by a relative error about +/-1 %. As the model is identified in two parts, it is essential to verify the coherence between the results. The Figure 15 shows the impedance spectra measured with the impedance analyser and the spectra of the diffusive impedance simulated after identification of CPE 3 .

Validation on current profile

The validation is realised at SOC=40% with ΔSOC=0. The measured and simulated voltages are represented in Figure 16. As in Figure 14, the quality of the result is attested by a relative error smaller than +/-2%.

DISCUSSION

The main advantage of this work is the possibility of rapidly identifying a battery model which varies with the value of current. The use of impedance spectroscopy at "high frequencies" for different values of current is a good alternative to the chronopotentiometry on pulses. Indeed, the determination of resistance from voltage profiles could be disturbed, as it was already underlined, by the rate of data acquisition, but also by the difficulty to separate the "fast" voltage response of the R 2 //CPE 2 from the "slow" one of the CPE 3 . Even with a very fast data acquisition system, the choice of the right value for the voltage used in the resistance calculus needs to have an idea of the time constant of the R 2 //CPE 2 . Otherwise, bad estimations are expected. In this work, we had opportunity to check that the voltage measured on the front of a pulse after a time equals to 5.τ 2 gives values for the resistance in agreement with the values obtained by spectroscopy (Cf. Figure 17). By this way, the voltage variation caused by R 2 //CPE 2 is finished and the one caused by CPE 3 is minimum. In the present method only one or two measurements on high current pulses are required to complete the medium current measurements by spectroscopy. If one has a powerful analyser (with a 100 A amplifier for instance), the spectroscopic measurements could be sufficient to completely identify R 2 (i). Some eventual improvements of the method must be studied. The present work does not take into account the experimental uncertainties in measurements and their influence on the final result. The hysteresis had been neglected. The results (Cf. Figure 9) show a 20 mV gap in OCV curve (for a voltage element of 3.7V) which is far less than hysteresis in NiMH battery for instance (about 60mV for a 1.3 V element). This result is deduced from measurements of the voltage at different values of SOC. Each measurement occurs after one hour of rest which is supposed to be long enough to reach thermodynamic equilibrium. New experiments with longer rest time will enable to verify if this apparent hysteresis is a real one. Furthermore, the influence of the temperature will be a part of a future work.

The good agreement between measured and simulated voltage in the Figure 16 tends to prove the validity of the method, but it would be a good idea to go deeper into the question of coherence between the spectrum measured and the one simulated (Cf. Figure 15). Are the differences due to inappropriate measurement conditions (low data acquisition rate during the pulse for example) or to physical reasons (like influence of the level of current on the diffusive impedance)?

1 Figure 1 .

 11 Figure 1. Simplified electric circuit model and its representation in Nyquist plan

  one can see on a Bode diagram (Fig.2), the gain of the approached line is fixed by γ when the pulsation tends to zero, i.e.

Figure 2 .

 2 Figure 2. Bode diagram of the approximation of the diffusive impedance by 5 circuits R//C. In summary, a CPE which is defined by and α is approximated in the frequency band [ω min , ω max ] by an impedance defined by four parameters (ω 1 , ω 5 , γ, p).

2 )Figure 3 .

 23 Figure 3. Representation of Z2 in the Nyquist plan.

  values can be obtained by measurement. The real part is equal to R 1 + 0.5.R 2 and enables to evaluate R 2 . The characteristic pulsation 2 , c ω is linked to the other parameters [6]:

  obtained from the following expression of the imaginary part of the top of the semi-circle deduced from eq.6 and eq.7:

2 Q 1 . 4

 214 Non-linearity of R 2

Figure 4 .

 4 Figure 4. Representation of the fitting function of R 2 .

2 QFigure 6 .Figure 7 .

 267 Figure 6. Techniques for identification of R 1 and R 2 //CPE 2 .

Figure 8 .

 8 Figure 8. Measured voltage of a lithium battery on the battery tester during a discharge at 150 A for SOC of 80%.

Figure 9 .

 9 Figure 9. Measured open circuit voltage with verification of hysteresis • Second step: characterization of the parameters of the CPE 3 by chronopotentiometry.

Figure 10 .

 10 Figure 10. Voltage drop response measured on battery with and without OCV correction.

Figure 11 .

 11 Figure 11. Values of R 1 , R 2 , Q 2 , CPE 2 , τ 2 and α 2 as function of the current level for different SOC.

2 Figure 12 .

 212 Figure 12. R 2 measured and fitted as function of the current level for different SOC.

2 Q 2 Q

 22 The optimized values for the parameters are A = 35.87 A, B= 10.73 V -1 . Figure13shows the variation of C which precise values are: C(80) = 1.05 mΩ, C(60) = 1.06 mΩ, C(40) = 1.18 mΩ, C(20) = 1.57 mΩ.

1 and ω 5 (

 5 . Moreover, for this first study, we also imposed the value p of the slope at 0.5. Consequently, only the value of γ  is adjusted depending on the value of SOC. The values of γ and C obtained at different SOC are represented in Figure13.

Figure 13 .

 13 Figure 13. Variations of the values of C and γ as function of SOC

Figure 14 .

 14 Figure 14. Simulated and measured voltages after the optimisation on the 5 min pulse at 1C for different SOC.

Figure 15 .

 15 Figure 15. Simulated and measured impedance spectra after the optimisation on the 5 min pulse at 1C.

Figure 16 .Figure 17 .

 1617 Figure 16. Validation of the battery model on a profile of current

ACKNOWLEDGMENTS

The authors want to thank ADEME (Agence pour le Développement et la maîtrise de l'Energie) which financed a PhD grant. Thanks are also due to JP DIARD from INPG (Institut National Polytechnique de Grenoble) and F BADIN from IFP (Institut Français du Pétrole) who gave some useful advice. We also want to thank J SCORDIA from LMS-Imagine for his help. Discussions with other partners of the program SIMSTOCK were helpful too.