On the well-posedness of a multiscale mathematical model for Lithium-ion batteries

We consider the mathematical treatment of a system of nonlinear partial differential equations based on a model, proposed in 1972 by J. Newman, in which the coupling between the Lithium concentration, the phase potentials and temperature in the electrodes and the electrolyte of a Lithium battery cell is considered. After introducing some functional spaces well-adapted to our framework we obtain some rigorous results showing the well-posedness of the system, first for some short time and then, by considering some hypothesis on the nonlinearities, globally in time. As far as we know, this is the first result in the literature proving existence in time of the full Newman model, which follows previous results by the third author in 2016 regarding a simplified case.


Introduction
Let us suppose we have a mathematical system of differential equations involving, after a homogenization process (or other averaging methods), two different scales, that we may call macro and micro, for simplicity. A macro-scale domain is given by x ∈ (0, L), where different processes may occur in three relevant subintervals: (0, L 1 ), (L 1 , L 1 + δ) and (L 1 + δ, L). At each point x ∈ (0, L 1 ) ∪ (L 1 + δ, L) we consider that there is a micro-scale domain given by a sphere with radius R(x) > 0. Different processes are modeled in each scale with differential equations, which are coupled at different levels, including boundary conditions. Let us considered, motivated by the modeling of charge transport in Lithium batteries (as we will show below) the following (relatively) abstract framework (which can be generalized to other cases).

) and under which conditions is it global in time?
We study here a particular interesting case arising in the modeling of Lithium batteries. Lithium-ion batteries are currently used extensively for storing electricity in mobile devices, from phones to cars. Nevertheless, they are still many drawbacks for them, as their reduced charge capacity, the long time needed to charge them, thermal runaways, etc. Therefore, an intensive research is being done in order to improve these devices. A good mathematical model of the transport mechanisms within batteries is very important in that research in order to understand the physics involved in the processes and to allow quick numerical experiments. But this models need to be well understood from a mathematical point of view, so that conclusions extracted from them are more rigorous. The well-posedness of the mathematical models being used is, therefore, one of the key points to be considered. This is, indeed, the goal of this paper.
In [31] a full model for Lithium ion batteries was presented, based on the well-known J. Newman model (see [27]), and partial results for the well-posedness were given. In this paper, we intend to complete those well-posedness results and study the regularity of the solutions. Let us write the complete system as presented in [31] but using in the electrolyte the electric potential measured by a reference Lithium electrode (φ e ) instead of its real electric potential. This is done because in electrochemical applications the potential in an electrolyte is typically measured by inserting a reference electrode of a pure compound, typically a Lithium electrode (see [7,33,35]). Local existence means that we are able to prove existence of solutions if the final time t end is "small enough". In a special case, we will show that t end can be taken as large as wanted, making the existence of solutions global in time.
A typical Lithium-ion battery cell has three regions: a porous negative electrode, a porous positive electrode and an electro-blocking separator. Furthermore, the cell contains an electrolyte, which is a con- centrated solution containing charged species that move along the cell in response to an electrochemical potential gradient. Let L be the length of a cell of a battery, let L 1 be the length of the negative electrode and let δ be the length of the separator. We assume the radius of an electrode particle to be R s (x), with A schematic representation of a cell is given in Figure 1. The unknowns in the system of equations that we study are: • the concentration of Lithium ions in the electrolyte c e = c e (x, t) in every macroscopic point x ∈ (0, L), • the concentration of Lithium ions in the electrodes c s = c s (r, t; x), which for every macroscopic point x ∈ (0, L 1 ) ∪ (L 1 + δ, L) is defined in a microscopic ball of radius R s (x), with r indicating the distance to the center. We assume radial symmetry in the diffusion taking place in each particle. Therefore, we do not need to consider angular coordinates, • the electric potential ϕ s = ϕ s (x, t) in the electrodes, • the electric potential measured by a reference Lithium electrode in the electrolyte, φ e = φ e (x, t), • the temperature T(t) in the cell. The system of equations is given by (1.5)-(1.9) below.
The transport of Lithium through the electrolyte can be modeled by the following macro-scale system of equations (conservation of Lithium in the electrolyte): , x ∈ (0, L), is the reaction current resulting from intercalation on Li into solid electrode particles and represents the Lithium flux between the electrodes and the electrolyte (which is a nonlinear function of all the variables, an expression that will be given later), D e > 0 is a diffusion function, α e > 0 is constant and c e,0 is a known initial state. The transport of Lithium through the electrodes can be modeled, at almost every macroscopic point x ∈ (0, L 1 ) ∪ (L 1 + δ, L), by the following micro-scale system of equations (conservation of Lithium in the electrodes), written in radial coordinates due to the radial symmetry: where D s > 0 is a diffusion coefficient and α s > 0, both with a constant value in (0, L 1 ) and another constant in (L 1 + δ, L). Furthermore, c s,0 is a radially symmetric initial state. Notice that problem (1.6) is written in spherical coordinates. The electric potential in the electrodes, ϕ s , and the electric potential measured by a reference Lithium electrode in the electrolyte, φ e , can be modelled, for each t ∈ (0, t end ), by the following equations expressing the conservation of charge in the electrolyte: and the conservation of charge in the electrodes where κ = κ(c e , T), given by a function κ ∈ C 2 ((0, +∞) 2 ), κ > 0, and σ ∈ L ∞ ((0, L 1 ) ∪ (L 1 + δ, L)) are conductivity coefficients (uniformly positive), α φ e ≥ 0 is a constant, f φ e ∈ C 2 (0, +∞), A is the cross-sectional area (also the correct collector area) and I, the input current, is a piecewise constant function defined for t ∈ [0, t I end ]. Naturally, if the input current is defined up to a time t I end , then we can only expect to solve the system up to a time t end ≤ t I end . Finally, we consider the temperature T = T(t) inside the battery, which we consider spatially constant. Its time evolution is given by , where α T is a positive constant, T amb represents the ambient temperature, T 0 ≥ 0 is an initial known temperature and F T is a continuous functional that represents the effect of the other variables over T (an expression that will be given later). F T may depend of the values at time t of functions c e , c s , φ e and ϕ s locally or globally in space (see (1.18)-(1.22)).
In [31], instead of equations (1.5), (1.7) and (1.8) the author considers the equations where ε e > 0 is constant in (0, L 1 ), (L 1 , L 1 + δ), (L 1 + δ, L) and ε s in (0, L 1 ), (L 1 + δ, L). This general case, which is not in divergence form for c e and φ e , introduces an extra level of difficulty we will not consider here. The same techniques we present in this paper apply to that case, but with the introduction of some additional technicalities (see, e.g., [4,16] and the references therein). Here we have considered ε e constant in (0, L), and therefore we can remove it by assimilating ε p−1 e in D e , ε −1 e in α e and ε p e in κ. We simplify the domain notation by introducing the following sets: In fact, we will consider R s constant in both (0, L 1 ) and (L 1 + δ, L) as in (1.4) (see Figure 2). where U is the open circuit potential and η the surface overpotential of the corresponding electrode reaction.
There is no common agreement on the structural assumptions of U. Some papers use a fitting function either polynomial [36] or exponential [6], whereas other authors propose an improved version with logarithmic behavior close to the limit cases [32,35]. We will start working in a general framework, which contains as a particular important example the Butler-Volmer flux (see [27,31,35]), in which the functions are taken as f φ e (c e ) = ln c e , (1.14) α a ∈ (0, 1), (1.15) α c ∈ (0, 1), (1.16) where δ 1 , δ 2 are positive and constant in each electrode, c s,max is a constant that represents the maximum value of c s , and α a and α c (dimensionless constants) are anodic and cathodic coefficients, respectively, for an electrode reaction. The function U was later proposed in [32] as where α, β, γ are positive functions in L ∞ (J δ ) and p is a smooth bounded function.
The following particular choice of F T can be considered (see, e.g., [31]): where ln is the natural logarithm and R f is the film resistance of the electrodes. We point out that T has been assumed to be uniform across the whole cell. Nonetheless, a more complete model (specially one dealing with temperature blow-up) with a heat diffusion equation could be used. The study of blow-up in this type of equations has been largely studied (see, for example, [2,12]). By using this model, the voltage at time t of the cell can be estimated (see [31]) as The model, as written above, presents a number of difficulties: pseudo-two-dimensional (P2D) model, elliptic equations coupled with parabolic equations, discontinuities in space of some of the involved functions, the fact that j Li is not smooth and not monotone, etc. As far as we know, there is no proof in the literature of the global existence of solution before this work.
In order to state the results, we will introduce the following definitions, which we introduce as a mathematical tool: , where t end is some constant t end ≤ t I end . If t end < t I end , then one the following conditions holds as t ↗ t end : T → +∞.
(1. 26) We remark that without Assumption 2.13 existence for potentials φ e , ϕ s can only be established up to a constant (as expected).
In the particular case of j Li being the Butler-Volmer flux, given by (1.10)-(1.17), under some physically reasonable hypothesis on the parameters (see Assumptions 2.19, 2.20, 2.21, 2.24 below), we will prove a general existence result, and characterize the nature of possible blow-up. It states that, under reasonable extra conditions, the non-physical obstructions for well-posedness c s,B → 0, c s,max or c e → 0, +∞ that appear in Theorem 1, are not the cause of the blow-up behavior of local solutions.
In the last part of this paper (see Section 2.6), we give an ad hoc modification of the system for which we can state a global existence theorem by removing, in two steps, the impediments to global existence of solutions given by (1.27). Actually, since the blow-up conditions (1.27) do not correspond to the physical intuition, it is likely that the solutions do not, in fact, develop this kind of behavior.
In a first stage we will find a bound for the flux with respect to ϕ s − φ e,Li , and show:

Regularity assumptions of the nonlinear terms and initial data
Our most general formulation in this paper concerns the case of considering the following regularity conditions on the data: Assumption 2.1. Let us assume the data: (2.1) Notice that this implies that the lateral limits f(t ± i ) exist, but need not coincide.
The key ingredient of our approach is to choose some suitable convex subsets of some appropriate functional spaces. The choice of spaces well adapted to the system is a delicate matter, and errors or loose approaches to this might result in incorrect results (for an explanation of this philosophy see, e.g., [9]). Let us define the following open convex sets where we expect to find the solutions: where H 1 (a, b) is the usual Sobolev space over the interval (a, b) (see, e.g., [10,30]). It is important to point the natural space in which we will look for the pair (φ e , ϕ s ).
Instead of narrowly focusing on (1.12)-(1.16), we shall state an assumption (sufficient to prove Theorem 1) satisfied by a broader family of functions: Finally, on the temperature term F T we will require the following:

Definition of weak solution
We introduce the natural space for radial solutions Remark 2.6. Even though we will always present the variables as (x, t), or (r, t; x), it will sometimes be mathematically advantageous to consider maps t ∈ [0, t 0 ] → u( ⋅ , t) ∈ X, where X will be a space of spatial functions, in which we will use the notations C([0, t 0 ]; X) or L 2 (0, t 0 ; X), depending on the regularity.

Definition of weak-mild solution
Dealing directly with the weak formulation is technically very difficult. On the other hand, solutions in the classical sense (having all the necessary derivatives) may not exist. There is an intermediate type of solutions, known as "mild solutions". As a general rule (and in particular this applies to our problem), any classical solution is a mild solution and any mild solution is a weak solution. Let us introduce this kind of solutions in the simplest case: the heat equation. Consider the problem in a bounded, smooth domain Ω ⊂ ℝ N . Let u 0 ∈ L 2 (Ω) and f ∈ L 2 ((0, t end ) × Ω). One can construct, as an intermediate step, the solution of the problem by considering the decomposition of L 2 (Ω) in terms of eigenfunctions of −∆. Let us write the unique solution of (2.9) as v(t) = S(t)u 0 . The operator S(t) is a semigroup (see [10]), and has some interesting properties we will not discuss. A solution u of problem the non-homogeneous problem (2.8) can be written, for every t ∈ [0, t end ], as This kind of solution is known as a "mild solution". As in [17], one can define the "Green operator" for problem (2.8) as the function In our problem we will need to work with a suitable Green operator associated to each of the equations. Assuming Assumptions 2.2, we will define several Green operators. For any t 0 > 0 we define (see [20]): as the solution of the problem For system (1.6) we will need to do some extra work due to the fact that the equation is only "pseudotwo-dimensional". First we define the solution of problem (1.6) for every x fixed, by solving the corresponding problem The next step is to consider the dependence on x. Therefore, we construct the Green operator associated to problem (1.6) collecting all x ∈ J δ : given by Finally, we consider the Green operator for system (1.9) as the function This operator is well defined and of class C 1 due to Assumption 2.5. It will be useful to introduce the following Nemistkij operators: (c e , c s,B , φ e , ϕ s , T) → j Li ∘ (c e , c s,B , φ e , ϕ s , T), (c e , c s,B , φ e , ϕ s , T) → j Li ∘ (c e , c s,B , φ e , ϕ s , T).
It is well known (see, e.g., [21]) that these operators are locally Lipschitz continuous and C 1 (in the sense of the Fréchet derivative), properties that will be used in the proof of Theorem 1, due to the regularity of the elements of the composition (i.e. (2.2) and (2.3)).  We say that the extension is "proper" if b > a. We say that an extension is "maximal" if it does not admit a proper extension.
Notice that the contribution of c s can be studied, basically, as a one-dimensional behavior on J δ . More precisely, we consider the Green operator on the boundary of the balls B R s (x) , which we will show that contains all the necessary information In this sense we can rewrite (2.10) in terms of the restriction c s,B , instead of c s , as follows: Define c s (x, r, t) to be the solution, for every x ∈ J δ , of the problem Since c s = c s,B on ∂B R s (x) , due to (2.12), we havẽ , it follows that c s and G c s ,t 0j Li are solutions of the same parabolic problem (2.13) and, by the uniqueness of the solutions, Li .
Therefore, (c e , c s , φ e , ϕ s , T) is a solution of (1.5)-(1.9) in the sense of Definition 2.8, up to time t 0 , which can be arbitrarily close to t end .

Assumptions and results regarding Theorem 1
It was first shown in [31] that the uniqueness of solutions (φ e , ϕ s ) of system (1.7)-(1.8) holds up to a constant relating the difference between φ e and ϕ s . To avoid this, we set the following assumption: The idea of the proof of Theorem 1 (which is done in Section 5) is the following. First we will show (see Proposition 2.16 and its proof in Section 3) that we can solve (1.7) and (1.8) to obtain (φ e , ϕ s ) if c s , c e , T are given, extending to the nonlinear case the study for the linearized equation proved in [31]. Then we will apply a fixed point argument to the evolution problems (1.5), (1.6) and (1.9) to obtain the conclusion. Due to this proposition we know that the following functions are well defined: where (φ e , ϕ s ) ∈ H 1 (0, L) × H 1 (J δ ) is the (unique) solution of (2.6)-(2.7) satisfying (2.15). Assuming now that I is a continuous function, i.e. I ∈ C([0, t I end ]), we define the Green operator, for t 0 < t I end ,

), T(t), I(t)).
In Section 3 we will prove that: Remark 2.18. Since we will allow for charge and discharge cycles, we allow for I to be piecewise continuous, and this is why we define the piecewise weak-mild solution (see Definition 2.9).
The proof of the local existence of solutions will be based on finding a unique fixed point, in C([0, t 0 ]; K X ) for t 0 small enough, of the operator problem (c e , c s,B , T) = (G c e ,t 0 (α a N j Li ,t 0 ), G c s,B ,t 0 (α s N j Li ,t 0 ), G T,t 0 ) ∘G ϕ,t 0 (c e , c s,B , T).

Assumptions and remarks regarding Theorem 2
In Theorem 1 one of the reasons of a finite existence time could be that c s,B → 0, c s,max or c e → 0. This conditions do not, a priori, pose a relevant physical problem since the battery may very well be completely full or empty. However, the generality of our setting allows for no better statement. Let us study the case which seems to be the most relevant for the modeling of Lithium-ion batteries by considering new assumptions. This is purely technical, but it seems reasonable since, in empirical cases in the literature, one typically has c s,max ∼ 10 −2 mol cm −3 . In particular, in [36] the authors take c s,max = 1.6 × 10 −2 mol cm −3 . Let us consider the following special nonlinear terms (which, for a broader generality, include some extra constants with respect to (1.12)-(1.16)), which we will take as assumption in Theorem 2.

Assumption 2.21. We assume that
where α a , α s , β a ∈ (0, 1), γ 1 , γ 2 > 0 and δ 1 (x), δ 2 (x) > 0 are constant in each electrode. Notice that γ 1 , γ 2 are constant but not dimensionless. Furthermore, we consider U slightly more general than in (1. On the exponents we will consider the following assumption (see Remark 2.26): Theorem 2 gives a sufficient condition so that, under the structural consideration of a potential U similar to the one considered by Ramos-Please [32], should there be a blow-up of the solution, this is not caused by the non-physical behavior c e → 0, +∞ or c s → 0, c s,max as t → t end . Remark 2.26 gives an intuition of why this happens. On the other hand, if this conditions are not satisfied, then the toy models (2.27)-(2.28) suggest that if the sufficient conditions in Theorem 2 are not satisfied, then it could happen that c e → 0, +∞ or c s → 0, c s,max . The suggestion of the potential in Ramos and Please [32] was motivated by physical consideration, and not by the mathematical theory. Nonetheless, this seems to be the exact right structure for the mathematical theory. This is a notion of robustness to this proposal.
Considering the parameters used by some authors (see, e.g., [36]) (whereα a ,α c are non-dimensional charge transfer coefficients, F is the Faraday's constant and R is the universal gas constant), we have Considering (1.17), as in [32], we have

Remark 2.26.
It is well known that equations of the form with q < 1 can have a solution u ̸ ≡ 0 such that {x ∈ (0, L) : u(x, t) = 0} has positive measure for all time t after an initial time t 0 . These regions are known as "dead cores" and a detailed analysis of this phenomenon can be found, for instance, in [14,15]. On the other hand, if we study the equation we see that, if q > 1, then the solution can blow up u → +∞ for a finite time, whereas there is no blow-up if q ≤ 1. As we will see in the proof of Theorem 2, we can bound c e and c s by solutions of problems of type (2.27), where conditions (2.23)-(2.26) define the exponents q appearing in (2.27) and (2.28). Taking into account these results, it seems that, if these conditions are not satisfied, it is likely that the solutions c e or c s reach 0 in finite time, which is a non-physical behavior.

Assumptions and remarks regarding Theorems 3 and 4
Tackling an example of global existence of solutions for the model seems to be a Herculean mission. However, if we assume some nicer behavior of the nonlinear term U, at least far from the natural "working" conditions of the system, then we can establish a result of global existence of solutions.

Truncated potential behavior
In this subsection we strive to make some small modifications to the problem so that we can show that the charge potentials do not blow up in finite time. First let us assume that for large difference of potentials the flux is constant. We slightly modify (2.21)-(2.22) to the following new hypothesis: c s , T)), (2.30) and H is a bounded smooth cut-off function of the exponential: where s ∞ is an arbitrarily large but fixed cut-off value, and ζ is such that ζ > 0.

Truncated temperature behavior
Due to the delicate interconnectedness of the different terms in F T , it is too difficult to prove a global existence theorem. Nonetheless, since many authors consider that the temperature is constant in the cell (see [13,18,33]), we will allow ourselves a substantial simplification on the structural assumption for F T , in order to obtain the global uniqueness result avoiding the appearance of possible blow-up phenomena, which are related with the potential of Lithium batteries for thermal runaway and explosion under high temperature operation (see [36]).

Green operator for the semilinear elliptic system (1.7)-(1.8). Proof of Propositions 2.16 and 2.17
The idea of the proof of Proposition 2.16 is to follow the idea for "low overpotentials" in [31], but replacing Lax-Milgram's theorem by the Brézis theory of "pseudo-monotone operators" (see [8]). Under some special assumptions, this result is sometimes referred to as Minty-Browder's theorem [11,26] (see, e.g., [34]). Uniqueness up to a constant was already proved in [31]. Let us define, for a given (c e , c s,B , T) ∈ K X , the following function: which corresponds to η| ϕ s −φ e,Li =0 . Since (c e , c s,B , T) is known, we can define, for x ∈ [0, L] and Φ ∈ ℝ, Notice that We also define j Li 0 (x) = j Li (x, 0), which corresponds to j Li | ϕ s −φ e,Li =0 .

Remark 3.2.
Of course, if (φ e , ϕ s ) is a solution of system (1.7)-(1.8) and C is a constant, then (φ e + C, ϕ s + C) is also a solution. However, there exists only one solution in X ϕ .

Remark 3.3.
The main part of the proof above was to apply the monotonicity of j Li with respect to ϕ s − φ e,Li . The idea behind this monotonicity method has to do with the convexity of the associated energy functional.
We have so far proved that the map G ϕ given by (2.16) is well defined. We prove now, applying the Implicit Function Theorem, that this map is C 1 .
Proof of Proposition 2.17. We will apply the Implicit Function Theorem for the Banach space-valued mapping F :X ×Ŷ →Ẑ, to solve for an operatorĜ ϕ : U ⊂X →Ŷ in an expression of the form The choice of functional spaces will bê We will then check that due to the definition of φ e,Li (see (2.14)). Notice thatX is an open set of a Banach space, and therefore we can consider the Implicit Function Theorem (see, e.g., [24]) in this setting. We will use the notation We consider the maps A 1 , A 2 , A 3 , A 4 :X ×Ŷ →Ẑ, given by Our definition of weak solution is precisely F(x,ŷ) = 0Ẑ .
The function A 4 is linear and continuous, therefore C ∞ . It is automatic to see that On the other hand, A 1 , A 2 are multilinear and continuous, and therefore of class C 1 . In particular, for y = (u, v),ȳ = (ū,v ) ∈ X Φ we have Since ∂j Li ∂η is of class C 1 (see (2.2) in Assumption 2.3), so is A 3 and where, for x ∈ J δ , Letx 0 = (κ, c e , c s,B , T) ∈X, and letŷ 0 = (u 0 , v 0 ) = (φ e,Li , ϕ s ) ∈ X ϕ be the solution found in Proposition 2.16.
Then g(x) > 0 is in C(J δ ). Therefore, there exists a constant g 0 such that g(x) ≥ g 0 > 0 in J δ . Then understood as a bilinear form is continuous and coercive in X ϕ × X ϕ . Therefore, by Lax-Milgram's theorem D (u,v) F(x 0 ,ŷ 0 ) is bijective. Then, by the Implicit Function Theorem applied to F, there exists a unique C 1 functionĜ ϕ : U → X ϕ , defined in a neighborhood U ofx 0 inX, such that (3.5) holds.
Since is C 1 . We finally consider the following translation (which is also of class C 1 ): Due to the uniqueness (up to a constant) result we proved in Proposition 2.16, is the map G ϕ | V (as constructed in (2.16) through Proposition 2.16). Thus, G ϕ is of class C 1 in a neighborhood of (c 0 e , c 0 s,B , T 0 , I 0 ). Since this argument holds over any point (c 0 e , c 0 s,B , T 0 , I 0 ) ∈ K X × ℝ, we have shown that G ϕ is of class C 1 over K X × ℝ. This implies thatG ϕ is also C 1 and it concludes the proof.
Therefore, as an immediate consequence of Proposition 2.17, we have the following lemma in terms of the space of piecewise continuous functions C part defined by (2.1):
This operator has a nice representation formula A number of properties can easily be derived from this expression. For instance, the continuous dependence with respect to the data: where ρ is an increasing, continuous function such that ρ(0) = 0. The term G c s,B ,t is a little trickier. First we recall (see [3,28]) that, for any q > 2, and ‖G c s ,R,t 0 (u 0 , g)‖ L ∞ (0,t 0 ;L ∞ (B R )) ≤ C(‖u 0 ‖ L ∞ (B R ) + ‖g‖ L q (0,t 0 ) ).
Therefore, due to the linearity of the equation, for x, y in the same connected component of J δ , we have This solves the problem of continuity of G c s g with respect to x via the continuous dependence of the operator. Since c s,0 is continuous, working in each component we can prove directly that is Lipschitz continuous. Furthermore, it is easy to check that We also have the following time estimate, for t 0 ≥ 0, By defining the vectorial Green operator for the evolutionary part due to the previous results, we have where ρ is a continuous function such that ρ(0) = 0.

Proof of Theorem 1
First, let us assume that I is continuous. Let us define the function f(c e , c s,B , φ e , ϕ s , T) = (α e (x)N j Li (c e , c s,B , φ e , ϕ s , T), α ϕ s N j Li (c e , c s,B , φ e , ϕ s , T), −hA s (T − T amb ) + F T ).
It is clear that f : K X → Y is locally Lipschitz continuous (due to the definition and regularity of N j Li and F T ). Let us define, for any t > 0, This operator is also locally Lipschitz continuous. Let x = (c e , c s,B , T). Then we can rewrite the fixed point problem (2.11) (which was our definition of a weak-mild solution of (1.5)-(1.9)) as is locally Lipschitz continuous (due to Lemma 3.4), we can set a bounded neighborhood U ⊂ K X around x(0) = (c e (0), c s,B (0), T(0)) ∈ K X and a bounded set V ⊂ Y such that the composition is globally Lipschitz continuous. Due to the continuity of G t and the fact that G 0 ≡ x(0) ∈ K X , there exists Therefore, due to (4.1), we obtain that, for t 2 = min{t 0 , t 1 } > 0, is a contracting map. Then we can apply the Banach Fixed Point Theorem to find a unique solution of problem (5.1). If I is C part , then one can simply paste the mild solutions from the different time partitions. In this sense there exists a unique "piecewise weak-mild solution". By applying classical results, there exists a maximal existence time t end ≤ t I end . If t end < t I end , then we have d(x(t), ∂K X ) → 0 or ‖x(t)‖ X → +∞ as t → t end . This is equivalent to (1.26) and the proof is complete.

Proof of Theorem 2
There are a number of papers studying the semilinear equation u t − ∆u = f(u) with Robin-type boundary conditions, and its eventual possible blow-up, depending on the growth of f . Some of the most classical results are due to Amann [1] (for some more recent works, see e.g., [29]).
Let us assume that the solution (c e , c s,B , φ e , ϕ s , T) is defined in [0, t 0 ), where 0 < t 0 < t I end , and assume that (1.27) does not hold as t ↗ t 0 . We will show that (1.26) does not hold, and therefore t end > t 0 . We can think of the problem written in the following way: where C 1 and C 2 are functions which we will show can be estimated. We can define (α φ e + μ(x, T(t)) T(t) ) ≥ 1, and define β 4 as a monotone increasing locally Lipschitz continuous function such that β 4 (0) = 0 and β 4 (s) ≥ max{s α cs +γ 2 λ max , s α cs +γ 2 λ max }, then, by definingc s as the solution of we arrive at the definition of the searched function c s = c s,max −c s . We havec s > 0 (due to Assumption 2.24) andc s is a subsolution of (6.1). Hencec s ≤ c s,max − c s . Therefore, c s ≤ c s < c s,max . Putting these two bounding solutions together, we conclude 0 < c s ≤ c s ≤ c s < c s,max . Hence (1.26) does not hold as t ↗ t 0 . Applying Theorem 1, we obtain that t end > t 0 . Therefore, either t end = t I end or (1.27).

On the truncated potential case. Proof of Theorem 3
Assume that (1.28) does not hold as t → t 0 with 0 < t 0 < t I end . We can substitute the constants in the proof of Proposition 3 and repeat the argument. We get good bounds for c e , c s,B for t ∈ [0, t 0 ], which do not depend on ‖ϕ s −φ e,Li ‖ L ∞ . Hence, applying Lemma 3.4, we can obtain some estimates of ϕ s and φ e,Li in [0, t 0 ]. Therefore, by Theorem 2, we have that t end > t 0 . Hence, by contraposition, if t end < t I end , then (1.28) must hold.

On the truncated temperature case. Proof of Theorem 4
It is immediate to establish global sub-and supersolutions for T, which ensure (1.28) does not happen, and hence we get the global existence of solutions.

Mass conservation and derived properties
We begin by analyzing the compatibility conditions for the elliptic problems with Neumann boundary conditions. The compatibility conditions for the existence of ϕ s are Hence, we deduce that This is also the compatibility condition for system (1.8).
Let us take a look at the mass balance condition. From system (1.6) we observe that Thus, as expected, the mathematical model satisfies that the total concentration in the electrolyte is constant.
Since c e ≥ 0 (we have constructed the convex spaces K X and K Z so that the solution of the mathematical model satisfies this condition, which is consistent with the physics of the problem since it is a concentration), this conclusion also implies that ‖c e ‖ L 1 (0,L) is constant. In other words, as a whole, the electrolyte is saturated of Lithium. Any Lithium contribution from an electrode is immediately compensated somewhere else.
In fact, if one considers the Lithium concentration in the anode and cathode, for the pseudo-two-dimensional (P2D) model considered here, one gets If R 2 s,+ α s,+ = R 2 s,− α s,− , then the P2D model satisfies that any lose of Lithium in its anode is instantaneously received by its cathode and viceversa. Otherwise, there would appear to be a net variation in the total amount of Lithium in the solid phase.
The state of the charge at time t of the cell, SOC(t), can then be estimated as a normalized average of the mass of Lithium in one of the electrodes, typically the negative one (see [31]): c s,max dr dx.

Model validation and identification of parameters
In order to validate the model and identify suitable parameters, it is often necessary to compare the solutions of the model with real-world data and minimize some error functional using suitable parameter identification techniques and minimization algorithms (see, e.g., [19,22,23]). Since internal properties of the battery, such as the spatial distribution of Lithium ions, cannot be measured during operation, it is standard to measure only the temperature T (on the outside of the battery, but it is a safe assumption that the temperature is spatially homogeneous) and the output voltage V (see, for example, [5,37]), which can be estimated by using (1.23).

Conclusions
In this work, we have shown that the system of equations (1.5)-(1.9) has a unique solution (under some mild conditions). In the most general setting, only local in time existence result can be shown, and the nature of the possible blow-up is characterized. The most difficult part of the study of an eventual possible blow-up is the structure of the open circuit potential U. By considering the model of U given by (2.18), which was proposed in [32], we can rule out some non-physical behaviors as the possible causes of blow-up, by making some assumptions on the coefficients. Finally, for the sake of mathematical completeness, we provide some extra conditions that allow for global uniqueness in time and make some comments regarding the conservation of Lithium in the model and the model validation.