ar X iv :1 50 6. 00 60 5v 1 [ m at h. A P] 2 9 M ay 2 01 5 On the well-posedness of a mathematical model for Lithium-ion batteries∗ Angel Manuel Ramos Instituto de Matemática Interdisciplinar Departamento de Matemática Aplicada Facultad de Ciencias Matemáticas Universidad Complutense de Madrid Plaza de Ciencias 3, Madrid, 28040, Spain angel@mat.ucm.es May 29th, 2015† Abstract In this article we discuss the well-posedness of a mathematical model that is used in the literature for the simulation of Lithium-ion (Li-ion) batteries. First, a mathematical model based on a macro-homogeneous approach is presented, following previous works. Then it is showed, from a physical and a mathematical point of view, that a boundary condition widely used in the literature is not correct. Although these errors could be just sign typos (that can be explained as carelessness over d/dx versus d/dn, with n the outward unit vector) and authors using this model prob- ably use the correct boundary condition when they solve it in order to do simulations, readers should be aware of the right choice. Therefore, the deduction of the correct boundary condition and a mathematical study of the well-posedness of the corresponding problem is carried out here. Keywords: Lithium-ion batteries; mathematical model; partial differential equations, boundary conditions; uniqueness of solution; existence of solution 1 Introduction Lithium-ion (Li-ion) batteries have become very popular in the last years as a source of energy in multiple portable electronic devices. A mathematical model ∗Accepted for publication in Applied Mathematical Modelling †A previous version of this work can be seen here: http://www.mat.ucm.es/deptos/ma/prepublicaciones/2013/2013-08p.pdf 1 http://arxiv.org/abs/1506.00605v1 http://www.mat.ucm.es/deptos/ma/prepublicaciones/2013/2013-08p.pdf showing the key factors of the battery operation can be very helpful for the design and optimization of new models and also for the real time control of its performance. Based on a macro-homogeneous approach developed by Newman (see [8]) several mathematical models have been developed for these purposes (see [2, 4, 5, 3, 11, 12, 10, 1, 6, 7]) which include the main physics present in charge/discharge processes. A fully mathematical model is presented in Sec. 2, including a system of boundary value problems for the conservation of Lithium and conservation of charge in the solid and electrolyte phases, together with an initial value problem for the conservation of energy. In the literature (see [1], [6], [9], [10], [11], [12]) one can find numerical computations of the model (or simplifications of it), with parameters corresponding to actual devices, that help to highlight the structure of this highly coupled model and show the relevance to the applications. Over the last years, several authors have written articles in different journals including a boundary condition that is not correct (see Remark 5). Although the authors probably use the correct boundary condition when solving the model in order to do simulations, the reader should be aware of the right choice. In this article it is shown why that condition is not only physically incorrect (see Sec. 3) but also mathematically incorrect, since it is proved (see Sec. 4, Remark 7) that the corresponding system of boundary value problems modelling the conservation of charge does not have any solution (therefore the system is not well-posed). Some results regarding the uniqueness and existence of solution (with the correct boundary condition) for a simplified version of the model are shown in Sec. 4. 2 Mathematical model 2.1 Generalities A typical Li-ion battery cell has three regions: A porous negative electrode, a porous positive electrode and an electron-blocking separator. Furthermore, the cell contains an electrolyte, which is a concentrated solution containing charge species that can move all along the cell in response to an electrochemical potential gradient. The negative electrode is an intercalated lithium compound usually made from Carbon (typically graphite), with LiyC6 active material. Here y ∈ [0, 1] is the stoichiometry value of the material, which changes during charge and discharge. For instance, during discharge, lithium ions inside of solid LiyC6 particles diffuse to the surface where they react and transfer from the solid phase into the electrolyte phase (see [11, 12]). During charge they follow the opposite way. Theoretically, in a fully charged Li-ion battery this compound is LiC6 (lithiated graphite; Li saturation; y = 1), in a semi charged/discharged battery it is LiyC6, with y ∈ (0, 1) and in a fully discharged battery it is just carbon (Li depletion; y = 0). In practical operating cases y never attains the extreme values y = 0 or y = 1. The corresponding negative reaction equation 2 is the following: yLi+ + ye− + 6C←→ LiyC6 (0 ≤ y ≤ 1). The positive electrode is usually a metal oxide or a blend of multiple metal oxides (see [1]) such as lithium cobalt oxide (Li1−yCoO2), lithium iron phosphate (Li1−yFePO4), or lithium manganese oxide (Li1−yMn2O4), with y ∈ [0, 1]. For instance, during discharge the positively charged ions travel via diffusion and migration through the electrolyte solution to the positive electrode where they react and insert into solid metal oxide particles (see, e.g., [11, 12]). During charge they follow the opposite way. Theoretically, in a fully charged Li-ion battery y = 1 (Li depletion), in a semi charged/discharged battery y ∈ (0, 1) and in a fully discharged battery y = 0 (Li saturation). Again, in practical operating cases y never attains the extreme values y = 0 or y = 1. The corresponding positive reaction equation for the examples showed above are the following: LiCoO2 ←→ Li1−yCoO2 + yLi+ + ye− (0 ≤ y ≤ 1), LiFePO4 ←→ Li1−yFePO4 + yLi+ + ye− (0 ≤ y ≤ 1), LiMn2O4 ←→ Li1−yMn2O4 + yLi+ + ye− (0 ≤ y ≤ 1). Therefore, considering both, negative and positive reaction, the corresponding total reaction equations are the following: LiCoO2 + 6C←→ Li1−yCoO2 + LiyC6 (0 ≤ y ≤ 1), LiFePO4 + 6C←→ Li1−yFePO4 + LiyC6 (0 ≤ y ≤ 1), LiMn2O4 + 6C←→ Li1−yMn2O4 + LiyC6 (0 ≤ y ≤ 1). A 1D electrochemical model is considered for the evolution of the Li con- centration ce(x, t) (mol / m3) and the electric potential φe(x, t) (V) in the elec- trolyte and the electric potential φs(x, t) (V) in the solid electrodes, along the x–direction, with x ∈ (0, L) and L = L1 + δ + L2 being the cell width (m). We assume that (0, L1) corresponds to the negative electrode, (L1, L1 + δ) corre- sponds to the separator and (L1 + δ, L) corresponds to the positive electrode. This is coupled with a 1D microscopic solid diffusion model for the evolution of the Li concentration cs(x; r, t) in a generic solid spherical electrode particle (situated at point x ∈ (0, L1) ∪ (L1 + δ, L)) along the radial r–direction, with r ∈ [0, Rs] and Rs (m) the average radius of a generic solid active material particle. This 1D approximation is valid since the characteristic length scale of a typical Li-ion cell along the x-axis is on the order of 100 µm, whereas the characteristic length scale for the remaining two axes is on the order of 100,000 µm or more (see [1]). Rs can be different at each electrode and therefore we consider Rs = Rs(x) = { Rs,− if x ∈ (0, L1), Rs,+ if x ∈ (L1 + δ, L). 3 The percentage of available local energy at time t (s) and radius r of a generic negative electrode particle situated at point x ∈ (0, L1) in the cell x–direction is the the same as its stoichiometry value y = y(x, r, t), which can be computed as y(x, r, t) = cs(x; r, t) cs,−,max , where cs,−,max (mol m−3) is the maximum possible concentration in the solid negative electrode and cs(x; r, t) (mol m−3) is the the Li concentration at time t, radius r and point x. Therefore, the bulk state of charge (SOC) for the negative electrode (it can be also done for the positive electrode but both are related and, therefore, it suffices to use only one of them) is SOC(t) = 3 L1(Rs)3 ∫ L1 0 ∫ Rs,− 0 r2 cs(x; r, t) cs,−,max drdx. (1) We point out that SOC is an nondimensional quantity that can be used as an indicator of the available energy in the cell. Measuring cs(x; r, t) is not easy, but we can use the mathematical model below (see, e.g., [1, 2, 11, 12]) to compute it. From a theoretical point of view SOC could be a value between 0 and 1 but, as mentioned above, in practical operating cases it never attains the extreme values 0 or 1. Therefore intermediate values y0%, z0%, y100%, z100% are used for each electrode to refer to 0 % or 100 % SOC. 2.2 A full mathematical model Based on the models appearing in the literature (see [2, 4, 5, 3, 11, 12, 10, 1, 6, 7]) and assuming constant diffusion and activity electrolyte coefficients, a full mathematical model for the performance of a battery, also including heat transfer dynamics, is given by system of equations (2)–(6):              εe ∂ce ∂t −De ∂ ∂x ( εpe ∂ce ∂x ) = 1− t0+ F jLi, in (0, L)× (0, tend), ∂ce ∂x (0, t) = ∂ce ∂x (L, t) = 0, t ∈ (0, tend), ce(x, 0) = ce,0(x), x ∈ (0, L), (2)                    For each x ∈ (0, L1) ∪ (L1 + δ, L) : ∂cs ∂t − Ds r2 ∂ ∂r ( r2 ∂cs ∂r ) = 0, in (0, Rs)× (0, tend), ∂cs ∂r (x; 0, t) = 0, −Ds ∂cs ∂r (x;Rs, t) = Rs(x) 3εs(x)F jLi, t ∈ (0, tend), cs(x; r, 0) = cs,0(x; r), (3) 4              For each t ∈ (0, tend) : − ∂ ∂x ( εpeκ ∂φe ∂x ) + (1− 2t0+) RT F ∂ ∂x ( εpeκ ∂ ∂x ln ( ce ) ) = jLi in (0, L), ∂φe ∂x (0, t) = ∂φe ∂x (L, t) = 0, (4)                          For each t ∈ (0, tend) : −εsσ ∂2φs ∂x2 =− jLi in (0, L1) ∪ (L1 + δ, L), εs(0)σ(0) ∂φs ∂x (0, t) = εs(L)σ(L) ∂φs ∂x (L, t) = − I(t) A , ∂φs ∂x (L1, t) = ∂φs ∂x (L1 + δ, t) = 0, (5)    MCp dT dt = −hAs ( T − Tamb ) + qr + qj + qc + qe, t ∈ (0, tend), T (0) = T0, (6) In the above system of equations ce = ce(x, t) (mol / m3) at time t (s), cs = cs(x; r, t) (mol / m3), φe = φe(x, t) (V), φs = φs(x, t) (V), T = T (t) is the temperature of the cell (K), I = I(t) is the applied current (A), εe = εe(x) =    εe,− if x ∈ (0, L1), εe,sep if x ∈ (L1, L1 + δ), εe,+ if x ∈ (L1 + δ, L) is the volume fraction of the electrolyte, p is the Bruggeman porosity exponent (nondimensional constant), De is the electrolyte diffusion coefficient (m2 s−1), t0+ (dimensionless and assumed here to be constant) is the transference number of Li+, Ds = Ds(x) = { Ds,− if x ∈ (0, L1), Ds,+ if x ∈ (L1 + δ, L) is the solid phase Li diffusion coefficient (m2 s−1), εs = εs(x) = { εs,− if x ∈ (0, L1), εs,+ if x ∈ (L1 + δ, L) is the volume fraction of the active materials in the electrodes, κ = κ (ce(x, t), T (t)) is the electrolyte phase ionic conductivity (S m−1), A (m2) is the cross-sectional area (also the electrode plate area), σs = σs(x) = { σs,− if x ∈ (0, L1), σs,+ if x ∈ (L1 + δ, L) 5 is the electrical conductivity of solid active materials in an electrode (S m−1) and jLi = jLi ( x, φs(x, t), φe(x, t), cs(x;Rs(x), t), ce(x, t), T (t) ) is the reaction current resulting in production or consumption of Li (A m−3) at point x and time t. For jLi the Butler-Volmer equation is usually used (see, e.g., [2, 5, 11, 12, 10, 1, 7]): jLi ( x, φs, φe, cs, ce, T ) =            3εs(x) Rs(x) i0 [ exp ( αaF R T η ) − exp ( −αcF R T η )] if x ∈ (0, L1) ∪ (L1 + δ, L), 0 if x ∈ (L1, L1 + δ) (7) (here, for the sake of simplicity, we have considered the solid/electrolyte inter- facial film resistance to be zero and therefore is not included in the above equa- tion), where as(x) = 3εs(x) Rs(x) is the specific interfacial area of electrodes (m−1), i0 = i0(x, cs, ce), i0(x, cs, ce) = { k−(ce) αa(cs,−,max − cs) αa(cs) αc if x ∈ (0, L1), k+(ce) αa(cs,+,max − cs) αa(cs) αc if x ∈ (L1 + δ, L) (8) is the exchange current density of an electrode reaction (A m−2), k−, k+ are kinetic rate constants (A m−2+6αa+3αc mol−2αa−αc), αa, αc (dimensionless con- stants) are anodic and cathodic transfer coefficients for an electrode reaction, η = η ( x, φs(x, t), φe(x, t), cs(x;Rs(x), t), T (t) ) , η ( x, φs, φe, cs, T ) = { φs − φe − U(x, cs, T ), if x ∈ (0, L1) ∪ (L1 + δ, L), 0 if x ∈ (L1, L1 + δ) is the surface overpotential (V) of an electrode reaction and U(x, cs, T ) is the equilibrium potential (V) at the solid/electrolyte interface (i.e. the open circuit voltage - OCV). A way of expressing U is (see [6]): U(x, c, T ) =          U− ( c cs,−,max ) + ∂U− ∂T ( c cs,−,max ) ( T − Tref ) if x ∈ (0, L1), U+ ( c cs,+,max ) + ∂U+ ∂T ( c cs,+,max ) ( T − Tref ) if x ∈ (L1 + δ, L), where U−, U+, ∂U− ∂T ∂U− ∂T are functions tipically obtained from fitting experi- mental data and cs,+,max (mol m−3) is the maximum possible concentration in the solid positive electrode. 6 Remark 1 For low overpotential cases Butler-Volmer equation (7) can be sim- plified to the following linearized version (see [10]): jLi ( x, φs, φe, cs, ce, T ) =      3εs(x) Rs(x) i0 (αa + αc)F R T η, if x ∈ (0, L1) ∪ (L1 + δ, L), 0, if x ∈ (L1, L1 + δ). Regarding heat transfer Eq. (6), M (kg) is the mass of the battery, Cp (J kg−1 K−1) is the specific heat capacity, h (W m−2 K−1) is the heat transfer coefficient for convection, As (m 2) is the cell surface area exposed to the convec- tive cooling medium (typically air), Tamb is the (ambient) temperature of the cooling medium, qr = qr(t) = A ∫ L 0 jLiηdx is the total reaction heat, qj = qj(t), qj(t) = Aεs,−σ− ∫ L1 0 ( ∂φs ∂x )2 dx+Aεs,+σ+ ∫ L L1+δ ( ∂φs ∂x )2 dx +A ∫ L 0 [ εpeκ ( ∂φe ∂x )2 + (2t0+ − 1) R T F εpeκ ( ∂ ln ce ∂x )( ∂φe ∂x ) ] dx is the ohmic heat due to the current carried in each phase and the limited conductivity of that phase, qc = qc(t) = I(t)2 Rf A is the ohmic heat generated in the cell due to contact resistance between current collectors and electrodes, Rf (Ω m2) is the film resistance of the electrodes and qe = qe(t) = TA [ ∫ L1 0 jLi ∂U− ∂T ( cs(x;Rs,−, t) cs,−,max ) dx + ∫ L L1+δ jLi ∂U+ ∂T ( cs(x;Rs,+, t) cs,+,max ) dx ] is the reversible heat caused by the reaction entropy change (see [13] and [6] for a simpler formulation of qe(t)). Heat sources qr, qj and qc are always positive and, as explained in [11], the second term inside the last integral of qj is generally negative. On the other hand, qe can be either positive or negative. Remark 2 The second term on the left hand side of system (4) is often written in the literature using (1− t0+) (see Refs. [2] and [4]) or 2(1− t0+) (see Refs. [5, 7 11, 12, 10, 1, 7]), instead of (1 − 2t0+), which is the correct term (see pag. 250 of Ref. [8]). We remind that we are assuming here that t0+ is constant and it is given (see [8]) by t0+ = D0,+ D0,− +D0,+ , with D0,−, D0,+ transport coefficients related to the diffusion of Li+ cations and the corresponding anions, respectively. Remark 3 Many authors (see, e.g., [11] and [6]) use some physiochemical pa- rameters involved in the model depending on T and use the Arrhenius equation to express such a dependency, which needs the corresponding activation ener- gies. The dependency of some parameters on T can sometimes be estimated. For instance in Ref. [7] the authors give formulae for κ(ce, T ), De(c, T ) and t 0 +(c, T ). Remark 4 This system of equations does not have uniqueness of solution (if φs(x, t) and φe(x, t) are solutions then φs(x, t) + c(t) and φe(x, t) + c(t) are also solutions, for any function c(t)). A way of avoiding that is to set a refer- ence value of φs(x, t) or φe(x, t) at some point x. For instance we can impose φs(0, t) = 0 for any t ∈ [0, tend] (see Remark 12). Some results regarding the existence and uniqueness of solution are showed in Sec. 4. After solving the above model, we can get SOC(t) by using (1) and the cell voltage V (t) (A), at time t, given by V (t) = φs(L, t)− φs(0, t)− Rf A I(t). (9) The battery pack voltage is calculated by multiplying the single cell voltage of Eq. (9) by the number of serially connected cells in the battery. 3 Correct boundary conditions in (5) Remark 5 In [11, 12, 10, 7], authors use the following boundary conditions at x = 0 and x = L, instead of those presented in (5): − εs,−σ− ∂φs ∂x (0, t) = εs,+σ+ ∂φs ∂x (L, t) = I(t) A . (10) In [5] authors use the following boundary conditions at x = 0 and x = L: − εs,−σ− ∂φs ∂x (0, t) = εs,+σ+ ∂φs ∂x (L, t) = − I(t) A . (11) In Sec. 3.1 it is proved that both conditions are not correct from a physical point of view, because the resulting system of equations does not model a battery performance. Then, in Sec. 4.1 it is proved that they are also mathematically incorrect, because, in fact, the corresponding system of boundary value problems modelling the conservation of charge does not have any solution (therefore the system is not well-posed) unless I(t) ≡ 0 (see Remark 7). 8 3.1 Deduction of system (5) Given an external current I(t) (A) applied to the battery (I(t) > 0 when the battery is discharging), according to Kirchoff’s law I = Is + Ie, where Is, Ie are the current (A) in the solid electrode and in the electrolyte, respectively. Now, by Ohm’s law, Is(x, t) = −Aσ(x) ∂φs ∂x (x, t) = I(t)− Ie(x, t). (12) Then, for each t ∈ (0, tend), taking into account the porosity nature of the solid electrodes (only a fraction of its volume contribute to its conductivity) and using that ∂Ie ∂x (x, t) = AjLi, (13) the following equation for the conservation of charge in the electrode solid phase is obtained: −εsσ ∂2φs ∂x2 = −jLi, in (0, L1) ∪ (L1 + δ, L). We point out that Is(x, t) = 0 (and therefore Ie(x, t) = I(t)) for all x ∈ (L1, L1+ δ) and Ie(0, t) = Ie(L, t) = 0. This provides, using (12), the following boundary conditions: εs(0)σ(0) ∂φs ∂x (0, t) = εs(L)σ(L) ∂φs ∂x (L, t) = − I(t) A , ∂φs ∂x (L1, t) = ∂φs ∂x (L1 + δ, t) = 0, which completes the deduction of system (5). Therefore, if for instance the battery is discharging, jLi > 0 (respectively, jLi < 0) if x ∈ (0, L1) (respectively, if x ∈ (L1 + δ, L)) and the graph of φs(·, t) for some time t ∈ (t, tend) could be qualitatively (not necessarily quantitatively) similar to that showed in Figure 1. 3.2 Deduction of system (4) Neglecting deviations of the electrolyte solution from ideal behavior Ie(x, t) = −Aκ ∂φe ∂x (x, t) +A(1 − 2t0+) RT F κ ∂ ∂x ln ( ce(x, t) ) (see [8]). The first term on the right hand side is due to Ohm’s law (as in (12)) and the second term accounts for concentration variations. Then, for each t ∈ (0, tend), taking into account the porosity nature of the solid electrodes and using again (13), the following equation for the conservation of charge in the electrode solid phase is obtained: − ∂ ∂x ( εpeκ ∂φe ∂x ) + (1− 2t0+) RT F ∂ ∂x ( εpeκ ∂ ∂x ln ( ce ) ) = jLi, in (0, L), 9 0 0.2 0.4 0.6 0.8 1 −1 0 1 2 3 4 5 x/L φ s ( x, t) Figure 1: Typical Graph of φs(·, t) for some time t ∈ (t, tend). together with boundary conditions Ie(0, t) = Ie(L, t) = 0. i.e., −κ(ce(0, t)) ∂φe ∂x (0, t) + (1− 2t0+) 2T F κ(ce(0, t)) ∂ ∂x ln ( ce(0, t) ) = 0, −κ(ce(L, t)) ∂φe ∂x (L, t) + (1− 2t0+) 2T F κ(ce(L, t)) ∂ ∂x ln ( ce(L, t) ) = 0, which, taking into account the boundary condtions in (2), are equivalent to ∂φe ∂x (0, t) = ∂φe ∂x (L, t) = 0. 4 About existence and uniqueness of solution of system (4)–(5) Let us study some properties regarding the existence and uniqueness of solution of (4)–(5). We will only focus here on the equations involving conservation of charge. Even if coupling with the concentrations of Lithium is not included in the analytic work, the system under studied is interesting and not trivial at all. It is interesting because suitable numerical schemes of the full model will probably 10 need to solve this kind of systems at each time step of a time discretized version of the full model. It is non trivial because it is a system of three boundary value problems, defined and coupled over 3 different domains. Moreover, as we will see below, the system does not have uniqueness of solution except if the parameters involved in the system satisfy a compatibility condition (which is satisfied when the parameters come from the full model) and we add an extra boundary condition (which is equivalent to set a reference potencial value and, therefore, consistent with the physics of the problem). The existence result showed in Section 4.3 is not standard, since it requires the use of Functional Analysis techniques. 4.1 Preliminaries Let us suppose (ce, cs, φe, φs, T ) is a solution of the full model. Given t ∈ (0, tend) let us denote u(x) = φe(x, t) and v(x) = φs(x, t). Then, if E(x,w) = exp ( αaF RT ( w − f(x) ) ) − exp ( −αcF RT ( w − f(x) ) ) , u and v must satisfy the following system of equations equivalent to (4)–(5):                        − d dx ( r(x) du dx (x) ) + ν d dx ( r(x) d dx ln(c(x)) ) =          k−b−(x)E(x, v(x) − u(x)), if x ∈ (0, L1), k+b+(x)E(x, v(x) − u(x)), if x ∈ (L1 + δ, L), 0, if x ∈ (L1, L1 + δ), du dx (0) = du dx (L) = 0, (14)                    − d2v dx2 (x) =−q−b−(x)E(x, v(x) − u(x)), if x ∈ (0, L1), − d2v dx2 (x) =−q+b+(x)E(x, v(x) − u(x)), if x ∈ (L1 + δ, L), − dv ∂x (0) = d−, dv ∂x (L) = −d+, dv dx (L1) = dv dx (L1 + δ) = 0, (15) r(x) = εe(x) pκ(ce(x, t), T (t)), c(x) = ce(x, t), ν = (1− 2t0+) 2T (t) F , b±(x) = 3εs(x) Rs(x) (ce(x, t)) αa (cs,±,max − cs(x;Rs, t)) αa (cs(x;Rs, t)) αc , q± = k± ε±σ± , d± = 1 εs,±σ± I(t) A , f(x) = U(x, cs(x;Rs, t), T (t)). (16) 11 Remark 6 Applying the divergence theorem it is easy to show that any solution of (14)–(15) must satisfy k− ∫ L1 0 b−(x)E(x, v(x) − u(x))dx + k+ ∫ 1 L1+δ b+(x)E(x, v(x) − u(x))dx = 0, q− ∫ L1 0 b−(x)E(x, v(x) − u(x))dx = d−, q+ ∫ L L1+δ b+(x)E(x, v(x) − u(x))dx = −d+. Therefore, a necessary compatibility condition for the existence of solution of system (14)-(15) is k−d− q− = k+d+ q+ , which is true if (16) holds, because k±d± q± = I(t) A . Remark 7 If we use boundary conditions (10) or (11) the corresponding com- patibility condition would be k−d− q− = −k+d+ q+ , or equivalently I(t) = −I(t), which is only satisfied if I(t) ≡ 0. Since function r(x) may not be smooth (at least in a general situation), system (14)–(15) (or, equivalently, (4)–(5)) may not have a classical solution. Therefore, we look for what is commonly known as a weak solution. In order to find the requirements that this solution shoud satisfy, we use the following functional analysis framework. Let us consider the Hilbert space H = H1(0, L)× ( H1(0, L1) ∩H 1(L1 + δ, L) ) , with the norm given by ‖(u, v)‖2V = ‖u‖2H1(0,L) + ‖v‖ 2 H1(0,L1) + ‖v‖2H1(L1+δ,L) = ∫ L 0 [ u(x)2 + ( du dx (x) )2 ] dx + ∫ L1 0 [ v(x)2 + ( dv dx (x) )2 ] dx+ ∫ 1 L1+δ [ v(x)2 + ( dv dx (x) )2 ] dx. Here H1(a, b) = { φ ∈ L2(a, b) : dφ dx ∈ L2(a, b) } , where Lp(a, b) (with p ≥ 1) is the set of all measurable functions from (a, b) to IR whose absolute value raised to the p-th power has finite integral and dφ dx denotes the derivative of φ in the sense of distributions. 12 Multiplying (14) (duality product < ·, · >(H1(0,L))′×H1(0,L)) by a function ϕ ∈ H1(0, L) and similarly in (15) by a function ψ ∈ H1(0, L1) ∩H 1(L1 + δ, L) we get ∫ L 0 r(x) du dx (x) dϕ dx (x)dx + ∫ L1 0 dv dx (x) dψ dx (x)dx + ∫ L L1+δ dv dx (x) dψ dx (x)dx −k− ∫ L1 0 b−(x)E(x, v(x)−u(x))ϕ(x)dx−k+ ∫ L L1+δ b+(x)E(x, v(x)−u(x))ϕ(x)dx +q− ∫ L1 0 b−(x)E(x, v(x)−u(x))ψ(x)dx+q+ ∫ L L1+δ b+(x)E(x, v(x)−u(x))ψ(x)dx = ν ∫ L 0 r(x) c(x) dc dx (x) dϕ dx (x)dx − d−ψ(0) + d+ψ(1). This equality can be rewritten as a [(u, v) , (ϕ, ψ)] = l [(ϕ, ψ)] , where a [(u, v) , (ϕ, ψ)] and l [(ϕ, ψ)] are the left hand side and right hand side, respectively, of the previous expression. According to this, we have the following definition of a weak solution (that we will call solution in the following). Definition 8 A (weak) solution of system (14)–(15) is a couple of functions (u, v) ∈ H such that a [(u, v) , (ϕ, ψ)] = l [(ϕ, ψ)] , ∀ (ϕ, ψ) ∈ H. (17) 4.2 About the uniqueness of solution of systems (14)–(15) For simplicity let us study system (14)–(15), which is equivalent to system (4)– (5). Proposition 9 If (u1, v1), (u2, v2) are two solutions of (14)–(15), then there exists a constant s ∈ IR such that u2(x) − u1(x) ≡ s, ∀ x ∈ (0, L). and v2(x) − v1(x) ≡ s, ∀ x ∈ (0, L1) ∩ (L1 + δ, L) . Proof: Let us suppose (u1, v1), (u2, v2) are two solutions of (14)–(15), i.e., they satisfy (17). Then, using (ϕ, ψ) = (u2 − u1, v2 − v1) in (17), with vi(x) =        k− q− vi(x) if x ∈ (0, L1) , k+ q+ vi(x) if x ∈ (L1 + δ, L) , (i = 1, 2), 13 we have a [(u2, v2) , (u2 − u1, v2 − v1)]− a [(u1, v1) , (u2 − u1, v2 − v1)] = 0, which is equivalent to ∫ L 0 r(x) ( d(u2 − u1) dx (x) )2 dx + k− q− ∫ L1 0 ( d(v2 − v1) dx (x) )2 dx+ k+ q+ ∫ L L1+δ ( d(v2 − v1) dx (x) )2 dx +k− ∫ L1 0 b−(x) ( E(x, v2(x)−u2(x)) −E(x, v1(x) − u1(x)) )( ( v2(x)− u2(x) ) − ( v1(x)− u1(x) ) ) dx +k+ ∫ L L1+δ b+(x) ( E(x, v2(x)−u2(x)) −E(x, v1(x)− u1(x)) )( ( v2(x)− u2(x) ) − ( v1(x) − u1(x) ) ) dx = 0. Hence, from the monotonicity properties of the exponential function, ∫ L 0 r(x) ( d(u2 − u1) dx (x) )2 dx + k− q− ∫ L1 0 ( d(v2 − v1) dx (x) )2 dx+ k+ q+ ∫ L L1+δ ( d(v2 − v1) dx (x) )2 dx = 0 and ( v2(x) − u2(x) ) − ( v1(x)− u1(x) ) = 0 ∀ x ∈ (0, L1) ∩ (L1 + δ, L) , which implies that there exists a constant s ∈ IR such that u2(x)− u1(x) ≡ s, ∀ x ∈ (0, L) and v2(x) − v1(x) ≡ s, ∀ x ∈ (0, L1) ∩ (L1 + δ, L) . Corollary 10 If (u1, v1), (u2, v2) are two solutions of (14)–(15), then v1(x)− u1(x) = v2(x)− u2(x) for all x ∈ (0, L1) ∩ (L1 + δ, L) and, therefore, E(x, v1(x) − u1(x)) = E(x, v2(x)− u2(x)). Corollary 11 If we add the boundary condition v(0) = 0 (18) (or v(x) = s ∈ IR or u(x) = s ∈ IR, with x and s arbitrarily chosen), then if system (14), (15), (18) has a solution, it is unique. 14 Remark 12 According to the mathematical results showed in Proposition 9, Corollary 10 and Corollary 11, a way of getting the property of uniqueness of solution of system (4)–(5) (see non-uniqueness results in Remark 4) is setting a reference potential value for φs or φe for each t ∈ (0, tend). For instance φs(0, t) = 0 (or φs(x, t) = s(t) ∈ IR, x and s(t) arbitrarily chosen). (19) 4.3 About the existence of solution of system (14), (15), (18)) Let us study the existence of solution for the following linearized approximation (see Remark 1) of system (14), (15), (18):                        − d dx ( r(x) du dx (x) ) + ν d dx ( r(x) d dx ln(c(x)) ) =          λ−b−(x) (v(x)− u(x)− f(x)) , if x ∈ (0, L1), λ+b+(x) (v(x)− u(x)− f(x)) , if x ∈ (L1 + δ, L), 0, if x ∈ (L1, L1 + δ), du dx (0) = du dx (L) = 0, (20)                    − d2v dx2 (x) =−θ−b−(x) (v(x) − u(x)− f(x)) , x ∈ (0, L1), − d2v dx2 (x) =−θ+b+(x) (v(x) − u(x)− f(x)) , x ∈ (L1 + δ, L), − dv dx (0) = d−, dv dx (1) = −d+, dv dx (L1) = dv dx (L1 + δ) = 0, (21) v(0) = 0 (or v(x) = s ∈ IR, with x and s arbitrarily chosen), (22) where λ± = k± (αa + αc)F RT , θ± = q± (αa + αc)F RT . Remark 13 The results in Proposition 9 are also valid for the solutions of (20)–(21). Therefore, if (20)–(22) has a solution, then it is unique. Multiplypling in (20) (duality product < ·, · >(H1(0,L))′×H1(0,L)) by a func- tion ϕ ∈ H1(0, L) and similarly in the first and second equations of (21) by functions λ− θ− ψ and λ+ θ+ ψ, respectively, with ψ ∈ H1(0, L1)∩H 1(L1 + δ, L), the following equality is obtained: ∫ L 0 r(x) du ∂x (x) dϕ dx (x)dx + λ− θ− ∫ L1 0 dv dx (x) dψ dx (x)dx + λ+ θ+ ∫ L L1+δ dv dx (x) dψ dx (x)dx 15 +λ− ∫ L1 0 b−(x) (v(x) − u(x)) (ψ(x)− ϕ(x)) dx +λ+ ∫ L L1+δ b+(x) (v(x) − u(x)) (ψ(x) − ϕ(x)) dx = λ− ∫ L1 0 b−(x)f(x) (ϕ(x) − ψ(x)) dx+ λ+ ∫ L L1+δ b+(x)f(x) (ϕ(x) − ψ(x)) dx +ν ∫ L 0 r(x) c(x) dc dx (x) dϕ dx (x)dx − λ−d− θ− ψ(0) + λ+d+ θ+ ψ(1). This equality can be rewritten as a [(u, v) , (ϕ, ψ)] = l [(ϕ, ψ)] , where a [(u, v) , (ϕ, ψ)] and l [(ϕ, ψ)] are the left hand side and right hand side, respectively, of the previous expression. According to this, we have the following definition of a weak solution (that we will call solution in the following). Definition 14 A (weak) solution of system (20)–(21) is a couple of functions (u, v) ∈ H such that a [(u, v) , (ϕ, ψ)] = l [(ϕ, ψ)] , ∀ (ϕ, ψ) ∈ H. (23) It is easy to show that, in this case, a : H×H → IR is bilinear and continuous and l : H → IR is linear and continuous. In order to be able to apply Lax- Milgram Theorem we need to show that a : H ×H → IR is coercive, but this is not true (otherwise we would get uniqueness of solution of (20)–(21), which, as we pointed out in Remark 13, is not true). Let us try to overcome this problem by changing H by a suitable closed subspace. Let H = {ϕ ∈ H1(0, L) : ∫ L 0 ϕ(x)dx = 0} × ( H1(0, L1) ∩H 1(L1 + δ, 1) ) . It is known that {ϕ ∈ H1(0, L) : ∫ L 0 ϕ(x)dx = 0} is a closed subspace of H1(0, 1) and there exists a constant α > 0 such that, for any function ϕ in that set, ∫ L 0 r(x) ( dϕ(x) dx (x) )2 dx ≥ α ∫ L 0 ϕ(x)2dx. Proposition 15 There exists a unique couple of functions (u, v) ∈ H, satisfy- ing a [(u, v) , (ϕ, ψ)] = l [(ϕ, ψ)] , ∀ (ϕ, ψ) ∈ H. (24) Proof: It is easy to show that a : H × H → IR is bilinear and continuous and l : H → IR is linear and continuous. Hence, in order to be able to apply Lax-Milgram Theorem we need to show that a : H ×H → IR is coercive. Now, a [(ϕ, ψ) , (ϕ, ψ)] = 16 ∫ L 0 r(x) ( dϕ dx (x) )2 dx+ λ− θ− ∫ L1 0 ( dψ dx (x) )2 dx+ λ+ θ+ ∫ L L1+δ ( dψ dx (x) )2 dx +λ− ∫ L1 0 b−(x) (ψ(x) − ϕ(x)) 2 dx+ λ+ ∫ L L1+δ b+(x) (ψ(x) − ϕ(x)) 2 dx ≥ r 2 ∫ L 0 ( dϕ dx (x) )2 dx+ rα 2 ∫ L 0 ϕ(x)2dx+ λ− θ− ∫ L1 0 ( dψ dx (x) )2 dx + λ+ θ+ ∫ L L1+δ ( dψ dx (x) )2 dx+ λ−b− ∫ L1 0 ( ψ(x)2 + ϕ(x)2 − 2ϕ(x)ψ(x) ) dx +λ+b+ ∫ L L1+δ ( ψ(x)2 + ϕ(x)2 − 2ϕ(x)ψ(x) ) dx, with r = inf x∈(0,1) {r(x)} > 0 and b± = inf x∈(0,1) {b±(x)} > 0. Using Young inequal- ity, for any β− > 0 and β+ > 0 we get ( ψ(x)2 + ϕ(x)2 − 2ϕ(x)ψ(x) ) ≥ ( ψ(x)2 − β±ψ(x) 2 + ϕ(x)2 − 1 β± ϕ(x)2 ) . Then, if we take β− and β+ defined by rα 6 + λ±b± = λ±b± 1 β± , we have that β− ∈ (0, 1), β+ ∈ (0, 1) and a [(ϕ, ψ) , (ϕ, ψ)] ≥ r 6 ∫ L 0 ( dϕ dx (x) )2 dx+ rα 2 ∫ L 0 ϕ(x)2dx + λ− θ− ∫ L1 0 ( dψ dx (x) )2 dx+ λ+ θ+ ∫ L L1+δ ( dψ dx (x) )2 dx +λ−b−(1− β−) ∫ L1 0 ψ(x)2dx+ λ+b+(1− β+) ∫ L L1+δ ψ(x)2dx, ≥ min { r 6 , rα 2 , λ− θ− , λ+ θ+ , λ−b−(1− β−), λ+b+(1− β+) } ||(ϕ, ψ)‖2 H , which proves the coercivity property and finishes the proof. Theorem 16 If (u, v) ∈ H is the solution of (24), then (u+ γ, v+ γ) ∈ H is a solution of (20)–(21), for any γ ∈ IR. Proof: We have to prove that (u, v) = (u+γ, v+γ) satisfies (23) for any γ ∈ IR. Now, for any (ϕ, ψ) ∈ H , (ϕ, ψ) = ( ϕ− 1 L ∫ L 0 ϕ(x)dx, ψ − 1 L ∫ L 0 ϕ(x)dx ) ∈ H. 17 Then, using (24) and the fact that λ−d− θ− = λ+d+ θ+ = I(t) A , a[(u+ γ, v + γ), (ϕ, ψ)] = a [ (u, v), ( ϕ− 1 L ∫ L 0 ϕ(x)dx, ψ − 1 L ∫ L 0 ϕ(x)dx )] = l ( ϕ− 1 L ∫ L 0 ϕ(x)dx, ψ − 1 L ∫ L 0 ϕ(x)dx ) = l (ϕ, ψ) + ( λ−d− θ− − λ+d+ θ+ ) 1 L ∫ L 0 ϕ(x)dx = l (ϕ, ψ) , which concludes the proof. Corollary 17 There exists a unique solution of (20)–(22). 5 Conclusions A fully mathematical model for the simulation of a Lithium-ion battery is pre- sented, based on a macro-homogeneous approach developed by Newman (see [8]) It includes a system of boundary value problems for the conservation of Lithium and conservation of charge in the solid and electrolyte phases, together with an initial value problem accounting for the conservation of energy. These model can be very helpful for the design and optimization of new batteries and also for the real time control of its performance. We point out that, over the last years, several authors have used similar models with an incorrect boundary condition. In this article it is showed, from a physical and a mathematical point of view, why such a boundary condition is incorrect. To show that, we prove that the system of equations including that boundary conditions does not have any solution. Then, we deduce the correct boundary condition and we show some results regarding the uniquesness and existence of solution of the corresponding problem. Acknowledgment This work was carried out thanks to the financial support of the Spanish “Min- istry of Economy and Competitiveness” under project MTM2011-22658; the research group MOMAT (Ref. 910480) supported by “Banco Santander” and “Universidad Complutense de Madrid”; and the “Junta de Andalućıa” and the European Regional Development Fund through project P12-TIC301. This pub- lication was also supported by OCIAM (University of Oxford), where the author was collaborating as an academic visitor. 18 References [1] N. A. Chaturvedi, R. Klein, J. Chrisensen, J. Ahmed and A. Kojic, Algo- rithms for Advanced Battery-Management Systems. Modeling, Estimation, and Control Challenges for Lithium-Ion Batteries, IEEE Control Systems Magazine June (2010) 49–68. [2] M. Doyle, T. F. Fuller and J. Newman, Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell, J. Electrochem. Soc. 140 (1993) 1256–1533. [3] T. W. Farrell and C. P. Please, Primary Alkaline Battery Cathodes. A Simplified Model for Porous Manganese Oxide Particle Discharge, J. Elec- trochem. Soc. 152 (2005) A1930-A1941. [4] T. F. Fuller, M. Doyle and J. Newman, Simulation and Optimization of the Dual Lithium Ion Insertion Cell, J. Electrochem. Soc. 141 (1994) 1–10. [5] P. M. Gomadam, J. W. Weidner, R. A. Dougal and R. E. White, Math- ematical modeling of lithium-ion and nickel battery systems, Journal of Power Sources 110 (2002) 267–284. [6] M. Guo, G. Sikha and R. E. White, Single-Particle Model for a Lithium-Ion Cell: Thermal Behaviour, J. Electrochem. Soc. 158 (2011) 122–132. [7] G. H. Kim, K. Smith, J. Ireland and A. Pesaran, Fail-safe design for large capacity lithium-ion battery systems, Journal of Power Sources 210 (2012) 243–253. [8] J. S. Newman, Electrochemical Systems (Prentice Hall, 1973). [9] R. T. Purkayastha amd R. M. McMeeking, An integrated 2-D model of a lithium ion battery: the effect of material parameters and mor- phology on storage particle stress, Comput Mech (2012) 50:209227 DOI 10.1007/s00466-012-0724-8. [10] K. Smith, C. D. Rahn and C. Y. Wang, Control oriented 1D electrochemi- cal model of lithium ion battery, Energy Conversion and Management 48 (2007) 2565–2578. [11] K. Smith and C. Y. Wang, Power and thermal characterization of a lithium- ion battery pack for hybrid-electric vehicles, Journal of Power Sources 160 (2006) 662–673. [12] K. Smith and C. Y. Wang, Solid-state diffusion limitations on pulse op- eration of a lithium ion cell for hybrid electric vehicles, Journal of Power Sources 161 (2006) 628–639. [13] K. E. Thomas, C. Bogatu and J. Newman, Measurement of the Entropy of Reaction as a Function of State of Charge in Doped and Undoped Lithium Manganese Oxide, J. Electrochem. Soc. 148 (2001) A570–A575. 19 1 Introduction 2 Mathematical model 2.1 Generalities 2.2 A full mathematical model 3 Correct boundary conditions in (??) 3.1 Deduction of system (??) 3.2 Deduction of system (??) 4 About existence and uniqueness of solution of system (??)–(??) 4.1 Preliminaries 4.2 About the uniqueness of solution of systems (??)–(??) 4.3 About the existence of solution of system (??), (??), (??)) 5 Conclusions