Thermomechanics of hydrogen storage in metallic hydrides: modeling and analysis

A thermodynamically consistent mathematical model for hydrogen adsorption in metal hydrides is proposed. Beside hydrogen diffusion, the model accounts for phase transformation accompanied by hysteresis, swelling, temperature and heat transfer, strain, and stress. We prove existence of solutions of the ensuing system of partial differential equations by a carefully-designed, semi-implicit approximation scheme. A generalization for a drift-diffusion of multi-component ionized"gas"is outlined, too.


Introduction
Hydrogen can be produced from a variety of renewable resources or in modern 4th-generation nuclear reactors operating at high temperatures where hydrogen production by water hydrolysis advantageously serves also to their cooling during periods when electricity cannot be produced. Then it is utilized in high-efficiency power generation systems with no emission of pollutants based on thermo-chemistry (burning directly hydrogen) or electro-chemistry (using fuel cells, cf. Section 5 for little more details). Hydrogen contains more energy per unit mass than any other available fuel. However, being the lightest element of the Periodic Table, it is highly volatile. Thus, in order to be compactly stored, standardly it is compressed in heavy high-pressure tanks or liquefied with recourse to expensive cryogenic systems. The lack of an efficient and economical way to store hydrogen is the major barrier to the massive commercial implementation of hydrogen-based technologies, especially in the automotive sector [11]. A promising alternative to cryogenic and high-pressure hydrogen storage option is provided by solid-state storage, a technology which exploits the property of certain metals and alloys to accommodate hydrogen atoms in their interstitial sites [18]. We propose a mathematical model for hydrogen adsorption in metals. Beside diffusion, the model accounts for phase transformation, temperature, strain, and hysteresis, cf. e.g. [10]. Thus, our model, based on a conventional rational mechanics (cf. Remark 2 below), extends those proposed and analyzed in [2,3,5], Since the modeling is entirely new, a detailed derivation is presented in the next Section 2 where also the model is a bit reformulated to facilitate mathematical analysis; of course, various simplifications had to been adopted, cf. Remark 1 below. Mathematical results as far as existence of weak solutions are summarized in Section 3 while their proof by a carefully designed semi-implicit discretisation in time is done in Section 4. Eventually, in Section 5, we briefly sketch the augmentation of the model for a multicomponent, charged (i.e. ionized) chemically reacting medium instead of mere single-component electro-neutral hydrogen, having in mind e.g. application to the mentioned fuel cells or to elastic semiconductors.

Model derivation
We consider a solid body, which we identify with a domain Ω of the threedimensional space. We regard Ω as a platform for several mutually interacting processes and phenomena affecting the kinetics of hydrogen adsorption/desorption [17,18]: • Phase transformation: at a low concentration, hydrogen atoms form a dilute interstitial solid-solution (α-phase). Increasing the hydrogen concentration causes parts of the solid solution to precipitate into a β-phase of larger interstitial concentration and lower density. Further addition of hydrogen takes place at constant pressure, until the metal is entirely converted into hydride.
• Temperature variation: hydrogenation is exothermic and reversible: when the metal is exposed to hydrogen, heat is generated; conversely, heating the hydride drives the reaction in the reverse direction.
• Strain and stress: hydrogenation is accompanied by large expansion of the unit cell volume of the crystal. Within this "swelling", volume changes between the two phases can vary from 8% to 30% and it may cause large stresses.
• Spatial distribution and transport : in addition, an important feature is distributed-parameter character of such storage devices. In particular, the motion of H atoms after dissociation of the H 2 molecule on the surface is diffusion driven by gradient of chemical potential, and heat transfer and force equilibrium must be properly counted.
• ϑ, the temperature field; The microstructural field is a collection of scalar variables which contains information concerning phase transformation and damage. We now derive a system of partial differential equations ruling the evolution of the primary fields. We do this in two steps.
Step 1: Balance laws. We invoke certain well-accepted thermomechanical principles, whose statement requires the introduction of some auxiliary fields: σ stress, f bulk force, f s surface force, s internal microforce, S microstress, g bulk microforce, g s external surface microforce, e internal energy, ε(u) = 1 2 ∇u + ∇u ⊤ small-strain tensor, ψ free energy, µ chemical potential, h hydrogen flux, h bulk hydrogen supply, h s surface hydrogen supply, q heat flux, q bulk heat supply.
Each particular specification of space-time evolution of primary and auxiliary fields constitutes a dynamical process. We require that every dynamical process comply with the following balance equations: . .
where the dot denotes the time derivative. The statements contained in (1) are, in the order: the standard-force balance, the microforce balance, the balance of mass for hydrogen, and the balance of internal energy. The corresponding natural conditions on ∂Ω are: Although the number of balance equations equals that of primary fields, the system (1) and (2) is under-determined. Such indeterminacy reflects the fact that these laws are common to a wide spectrum of thermomechanical systems. Thus, they cannot single out the particular mathematical model that best fits the system under investigation. One needs indeed additional conditions which can distinguish one particular material from another. These are called constitutive prescriptions.
Step 2: Second law. A constitutive prescription is typically expressed as a relation between the instantaneous value of a secondary field at a given point and that of a so-called constitutive list, a list of quantities obtained by taking the values of primary and secondary fields, or their space/time derivatives. A basic principle that guides the formulation of constitutive prescriptions is the requirement that every conceivable dynamical process be consistent with the entropy inequality: irrespectively of the practical difficulties involved in realizing such a process. Thus, unlike the balance laws (1), the imbalance (3) is not explicitly stated in the mathematical model, but it is implicitly enforced through a suitable choice of constitutive prescriptions. The entropy inequality is best exploited by replacing, in the list of fields to be specified constitutively, the internal energy with the free energy: Rewriting (1d) in terms of ψ and s, and substituting it into (3), one arrives at: where we have used the shorthand notation ε = ε(u). A standard argument due to Coleman and Noll [6] allows us to conclude that the free energy may depend at most on (ε, χ, ∇χ, c, ϑ): Moreover, if one assumes that entropy and chemical potential depend on the same list, one obtains The dissipation inequality can further be written in a more compact form by introducing the splitting: With that splitting, one indeed obtains: Step 3: Constitutive equations. To facilitate mathematical analysis but still capturing desired features, we restrict our attention to the following special constitutive ansatz: where λ > 0 is a length-scale parameter. This ansatz ensures, e.g., the heat capacity independent of the variables whose gradient is not directly controlled, i.e. ε and c, and also the chemical potential independent of ε and ϑ. The constitutive equations for entropy and chemical potential are On defining ω := ϕ − ϑ∂ ϑ ϕ, in view of (9) we have the constitutive equation for internal energy e = ψ + sϑ, cf. (4), is As constitutive prescriptions for the dissipative parts of the auxiliary fields we choose Here D is a 2nd-order positive-definite viscosity-moduli tensor, α > 0 counts for rate effects in evolution of χ, ζ is a convex function homogeneous of degree one; note that ζ is typically nonsmooth at 0, which counts for activation of evolution of χ. Moreover, K and M are respectively 2nd-order positive-definite heat-conductivity and hydrogen-diffusivity tensors. We also eventually set h = 0 and g = 0. We therefore arrive at the following system: . .
where ω is from (11). We make the following natural assumption which, in fact, says positivity of the heat capacity: Then, the inverse to ω(χ, ·) does exist and we can express ϑ as which allows us to eliminate temperature ϑ from the system (13). The symbol ∋ appearing in formula (13b) (and in formulas (18b), (28b), and (38b) below) means that the right-hand side is included in the left-hand side, which is a set, since ζ and ϕ 2 are non-smooth (see also (16b) below). Moreover, we will be a bit more specific in (9). A typical contribution to the free energy is where C is a 4th-order elastic-moduli tensor, E is a 2nd-order tensor which incorporates the effect of dilation due to the microstructural parameter χ within metal/hydride phase transformation, α is a tensor accounting for thermal dilation, φ 1 and φ 3 are the simplest variant of the contribution to the chemical potential and the heat capacity, respectively, and δ K is an indicator function of a convex set K ⊂ R N from which the phase-field χ is assumed to range. Assuming χ is scalarvalued, E the unit matrix, and k large, we get essentially an isotropic swelling controlled nonlinearly by the hydrogen concentration by χ ∼ a(c) while allowing still ϕ 1 (χ, ·) to be uniformly convex, as needed later in (45). Obviously, in view of (9), the specific choice (15) leads to Note that, in (16c), ∂ 3 ϑϑχ ϕ 3 ≡ 0, which makes the heat capacity independent of χ, but we can consider more generally the contribution φ 3 dependent also on χ to reflect different heat capacity of metal and of hydride, and therefore we do not restrict ourselves to a particular form of ϕ 3 in (9) but make only certain technical assumptions below, cf. (23k). Similarly, we keep the treatment of ϕ 1 in a nonspecified way. In fact, the specific form (16b) is a simplified linearization of the so-called St.Venant-Kirchhoff potential but, when derived from the St.Venant-Kirchhoff form by linearizing the stress response, it results still in some other terms, cf. [9, Sect. 5.4], which here was neglected rather for notational simplicity.
where N K = ∂δ K denotes standardly the normal cone to the convex set K. The boundary conditions (2) now take the form: Using the convention like u(·, t) =: u(t), we complete the system by the initial conditions: where we have set w 0 = ω(χ 0 , ϑ 0 ). We henceforth shall use the abbreviation for the so-called stored energy, i.e. the temperature independent part of the free energy: By testing (18a,b,c,d) respectively with . u, .
χ, µ, and by a constant ν, integrating by parts in time, using Green's formula with the boundary conditions (19), and taking into account (18e) so that we obtain the following identity: For ν = 0, the identity (22) represents the mechanical energy balance. For 0 < ν < 1, both the internal energy and dissipative terms are seen; henceforth, a discrete version of this estimate will be used in the proof of Lemma 4.2 for ν = 1/2. Eventually, for ν = 1, we recover the standard total-energy balance; note that the dissipative terms (and also adiabatic terms) then vanish.
Remark 1 Of course, the above model adopted a lot of simplifications of the actual situations in hydride storage. In particular the concept of small strains may be questionable at some situations, possible damage is essentially neglected, although formally it can be involved in a general form of ϕ 2 in (9) below but a lot of analytical considerations seem to be difficult to be straightforwardly adapted. Also temperature dependence of the chemical potential of hydrogen is neglected, i.e. φ 1 in (9) does not depend of θ. Further, the chemical reaction in the multi-component system metal/hydrid/hydrogen is basically neglected and hydride is modeled as a mixture of metal and hydrogen with essentially the possibility to obtain the same thermomechanical response of the phase transformation as the corresponding chemical reaction (and, in addition, we can easily model the activated hysteretic response related with this phase transformation).

Remark 2
The thermodynamics of our model follows essentially a classical approach based on rational mechanics and Clausius-Duhem inequality, cf. e.g. [4].
There are some variants of this general scenario [7,14,15,20] which are to some extent equivalent under some simplifications like those in Remark 1 above.

Weak solutions and their existence
Let us summarize the qualification on the data on which we will rely during the analysis of the initial-boundary-value problem (18). We confine ourselves to the (physically relevant) three-dimensional case. We consider a fixed time horizon T and abbreviate Q := Ω×(0, T ) and Σ := Γ×(0, T ), with Γ a boundary of the domain Ω ⊂ R 3 assumed Lipschitz. In the following, we shall use some classical notation for function spaces, in particular the Lebesgue spaces L p , the Sobolev spaces W k,p and, in particular, H k = W k,2 , and vector-valued functions. We suppose that Moreover, we need the qualification of the right-hand sides and the initial data: We note that (23l) will be used in the derivation of the estimate on ∇w, see (40f) below. Also note that (23a,d,e) is not in conflict with the example (16a) where Definition 3.1 (Weak solutions) We say that the six-tuple (u, χ, c, w, µ, ξ) with is a weak solution to the initial-boundary-value problem (18) The above definition obviously arises from (18)-(19) by using standard concept. The inequality (25b) arises by using additionally the identity Q ∆χ· .
and this solution is consistent with the energy conservation equation (22) with ν = 1, as well as with c ≥ 0 and w ≥ 0 on Q.
We will prove this theorem in Section 4 by a semi-implicit time discretisation in a more or less constructive manner, except the fixed-point argument behind the boundary-value sub-problems (28c)-(19c) and (18d)-(29d) and selection of converging subsequences. The additional properties (27) follow from (40f) and (46). The H 2 -regularity of χ is a standard consequence of (27d). For more detailed modes of convergence of the mentioned approximate solutions we refer to (49) and (50) below.

Analysis of (18)-(19) by semidiscretisation in time
We will prove existence of a weak solution to the initial-boundary-value problem (18) by a carefully constructed semi-implicit discretisation in time which, at the same time, will decouple the system to a sequence of convex minimization problems combined with a diffusion equation, and provide thus a rather efficient conceptual numerical strategy.
In comparison with the fully implicit time discretisation (i.e. the so-called Rothe method), our discretisation will allow for a simpler (and constructive) proof of existence of the discrete solutions, weaker assumptions about convexity mode of the stored energy, but we need to impose a bit stronger growth qualification of the data than required by the nature of the continuous system (18).
We use an equidistant partition of the time interval I = [0, T ] with a time step τ > 0, assuming T /τ ∈ N, and denote {u k τ } T /τ k=0 an approximation of the desired values u(kτ ), and similarly χ k τ is to approximate χ(kτ ), etc. Further, let us abbreviate by d k τ the backward difference operator, i.e. e.g. d k , etc. Then, using also notation (20), we devise the following semi-implicit discretisation: for k = 1, ..., T /τ , together with the boundary conditions starting from k = 1 by using An important feature of the scheme (28) is that it decouples to three boundaryvalue problems, which (after a further spatial discretisation) can advantageously be used in a numerical treatment and which is advantageously used even to show existence of approximate solutions: Proof. The first boundary-value problem arising by the decoupling is (28a,b) with (29a,b), and it leads to the minimization of the functional: Due to (16b) and (23d), ϕ 12 (·, ·, c k−1 τ ) and thus the whole functional (31) are strictly convex, for τ > 0 sufficiently small, namely for Therefore, there exists a unique minimizer (u k τ , χ k τ ) ∈ H 1 (Ω; R 3 )×H 1 (Ω; R N ). By standard arguments, cf. e.g. [24,Chap. 2 and 4], this minimizer when completed by ξ k τ ∈ N K (χ k τ ) is a weak solution to (28a,b)-(29a,b). Moreover, ξ k τ ∈ L 2 (Ω; R N ) can be shown by the same arguments as in the Lemma 4.3 below.
Then one can solve (28c)-(29c), which represents a semi-linear diffusion equation. We observe that we can eliminate µ k τ because obviously which leads us to abbreviate Then (28c)-(29c) transforms to the semi-linear boundary-value problem: on Ω together with the boundary condition on Γ: Due to the fully implicit discretisation of ∂ c ϕ 1 (χ, c) which is needed for the a-priori estimates in Lemma 4.2 below, via (34) we inevitably obtain the dependence of ) on c k τ so that the problem (35) unfortunately does not have any potential. Anyhow, thanks to (23a), the 2nd-order tensor M 1 is uniformly positive definite. Also, by (23e) and (23f), M 2 is bounded. As a consequence, the nonlinear operator A : is coercive. Thanks to Assumption (23c), we have that A is also weakly continuous, existence of a solution c k τ ∈ H 1 (Ω) can be obtained using the Galerkin method and Brouwer fixed-point theorem; of course, the obtained solution needs not be unique. Testing (35) by [c k τ ] − and using (23e) implies c k τ ≥ 0 provided c k−1 τ ≥ 0, which yields non-negativity of the hydrogen concentration recursively for any k = 1, ..., T /τ by using (23o).
Let us also note that from (33) we obtain also ∇µ k τ ∈ L 2 (Ω; R 3 ). In particular, we have simply both Dε(d k τ u):ε(d k τ u)/(1+τ |ε(d k τ u)| 2 ) ∈ L ∞ (Ω) and M(ε(u k τ ), χ k τ , c k τ , w k−1 τ )∇µ k τ ·∇µ k τ /(1+τ |∇µ k τ | 2 ) ∈ L ∞ (Ω), and thus the right-hand side of (28d) is in L 2 (Ω). Therefore, eventually, we are to solve (28d)-(29d), which represents a semilinear heat-transfer equation with the right-hand side in H 1 (Ω) * . The only nonlinearity is due to the w-dependence of L(ε(u k τ ), χ k τ , c k τ , w), σ a (χ k−1 τ , w), and s a (χ k−1 τ , w). The later two are needed to guarantee w k τ ≥ 0. Anyhow, since this nonlinearity is of lower order, we can pass through it by compactness and strong convergence. Thus, it suffices for us to check coercivity of the underlying operator. To this aim, we test (28d) by w k τ . The terms on the right-hand side containing w k τ are estimated standardly by using Hölder's and Young's inequalities, and using the qualifications (23l,m). Having coercivity, we see that there exists at least one solution. Moreover, this solution satisfies w k τ ≥ 0, which can be seen by testing (28d) by the negative part of w k τ and using that σ a (χ, w) = 0 and s a (χ, w) = 0 for w ≤ 0. Note that to have such non-negativity it is important to have w k τ in the K-term on the left-hand side, as well as in the nonlinear terms σ a and s a on the right-hand side.
In terms of these interpolants, we can write the approximate system (28) in a more "condensed" form closer to the desired continuous system (18), namely: . .

Lemma 4.2 (First estimates) Let again the assumptions of Lemma 4.1 hold.
Then, for some C and C r independent of τ > 0, Proof. The strategy is to test the particular equations in (28) respectively by d k τ u, d k τ χ, µ k τ , and 1 2 . For (28a,b), we note that a standard convexity argument yields: Owing the our equi-semiconvexity assumption (23d) on ϕ 1 (·, c), we also have: provided 0 < τ ≤ τ 1 with τ 1 from (32). Now, we add the tested equations together. We then use (41)-(42), and still the convexity (23g) of c → ϕ 1 (χ, c) to deduce the estimate where we used the abbreviation ϕ 12 from (20) with ϕ 2 from (16b). We remark that our semi-implicit scheme has benefited from the cancellation of the terms ± 1 τ ϕ 1 (χ k τ , c k−1 τ ) under this test by time-differences, which is a more general phenomenon to be understood as the fractional-step method, here combined with the semiconvexity, cf. [24,]. Now, by (23j), using Hölder inequality, and recalling that w k τ ≥ 0, we obtain the estimate: where ǫ is an arbitrarily small number, and C ǫ depends on ǫ. A similar estimate holds also for the terms multiplying d j τ χ on the right-hand side of (43). We now can adsorb the discrete time derivatives into the left-hand side, and use a discrete Gronwall inequality. Thanks to (23a) and (23g), we have ϕ 1 (χ, c) ≥ ǫc 2 and the a-priori bound χ k τ ∈ K. This gives the estimates (40a-c,e). The estimate (40d) follows from the relation (cf. also (33)): Eventually, we observe that the right-hand sides of (38d) and (39c) are bounded in L 1 (Q) and L 1 (Σ), respectively. For this, we need, in particular, assumptions (23j) and (23k). Then one can use the L 1 -theory for the heat equation to obtain the remaining estimate (40f), see [1]. To this goal, like in [12], one is to perform the test of (28d) by ̟ ′ (w k τ ), where ̟(w k τ ) := ((1+w k τ ) 1−ε − 1)/(ε−1) with ε > 0, and sum it for k = 1, ..., T /τ . Comparing to standard technique, see for instance [21,26,25], the only non-standard estimates is due to the L-term, as in [27], which requires assumption (23l), cf. also [24,Sect.12.9]. . .
Proposition 1 (Convergence for τ → 0) Let again the assumption of Lemma 4.1 hold. Then there is a subsequence such that and any (u, χ, c, w, µ, ξ) obtained in this way is a weak solution to the system (18)- (19) in accord with Definition 3.1 which also preserves the total energy as well as c ≥ 0 and w ≥ 0 as claimed in Theorem 3.2. Moreover, if Ω is smooth, then alsō c τ → c strongly in L 2 (I; H 1 (Ω)). (50b) Proof. For lucidity, we divide the proof into nine steps.
Step 3: Strong convergence of ε(u τ ). As the nonlinearities ∂ χ ϕ 2 , M, K, and L may depend on ε, we need first to prove strong convergence of ε(u τ ). Using the equation for the discrete approximation, and the limit equation proved in Step 2, we can write: where we used the monotonicity of u → Dε( . u) if the initial condition is fixed; cf. [24] for details about the time discretisation. The goal is to pass the right-hand side of (52) to 0. We use Aubin-Lions' Theorem to obtain strong convergence of . u τ in L 2 (Q; R 3 ), and by the Rellich compactness theorem also strong convergence of u τ (T ) in L 2 (Ω; R 3 ), which allows us to pass to the limit: the former equation is by discrete by-part summation, cf. e.g. [24,Remark 11.38], while the later equality simply follows by reversing integration by parts; here (51) is used. By (40c) and (46c), we have c τ → c weakly* in L ∞ (I; H 1 (Ω)), and .
Step 4: Limit passage in the micromechanical inequality. From (38b), we have ∀v ∈ L 2 (I; H 1 (Ω; R N )): The limit passage in (53b) is easy becauseξ τ → ξ weakly in L 2 (Q; R N ) andχ τ → χ strongly in L 2 (Q; R N ) has already been proved in Step 2; thus ξ ∈ N K (χ) is shown. Now we can make a limit passage in (53a). Here on the left-hand side we have collected all terms that need to be handled through a continuity or a weak upper semicontinuity arguments, while the right-hand side is to be treated by weak lower semicontinuity. We benefit from the strong convergence of χ τ (and similarly also of χ τ ) shown in Step 2 and of c τ andε τ proved in Step 3. The only nontrivial limit passage which, however, leads us directly to (25b) is: Eventually, the limit in Qξ τ · . χ τ dxdt is simple because, for anyξ τ ∈ N K ( . χ dxdt since we already know ξ ∈ N K (χ); note that (23n) has been used here. Thus (25b) is proved.
Step 5: Limit passage in the diffusion equation. Again, all strong convergences we already have proved in Steps 2 and 3 are to be used. Then the limit passage in the semi-linear equation (28c) is simple.
Step 6: Mechanical/chemical energy preservation. We balance the kinetic and stored energy integrated over the domain: as we actually did in (22) with ν = 0, i.e.
This is standardly achieved by testing the mechanical-chemical equations (18a,c), and the inclusion (18b), respectively by . u, µ, and .
χ, and by using the chain rule to integrate with respect to t. In order for these tests to be legal, we need ̺ .. u to be in duality with .
Step 8: Limit passage in the heat equation (38d). Having proved the strong convergence in Steps 3 and 7, the right-hand side of (38d) converges strongly in L 1 (Q) and this limit passage towards the weak solution to (18d)-(19d) is then easy.
w ∈ L 1 (I; H 3 (Ω) * ), cf. (46b) and realize the already proved identity (18d), which is in duality with the constant 1, we can perform rigorously this test and sum it with mechanical/chemical energy balance obtained already in Step 6.

Generalization of the model for electrolytes and fuel cell modeling
In a very basic scenario, the above presented model allows for a relatively simple generalization for a multicomponent, charged (i.e. ionized), chemically reacting medium undergoing electro-diffusion in elastic medium. When having in mind hydrogen fuel cells, the specific situation involves elastic polymeric porous layer with negatively-charged dopands and undergoing mechanical deformation/stresses e.g. due to swelling through which water H 2 O, hydrogen ions H + (i.e. protons), and hydronium ions H 3 O + (or, in general, ions of the type H 2n+1 O + n also with n ≥ 0) move by drift and diffusion; we speak about a polymeric electrolyte membrane (=PEM). This membrane is surrounded by two thin layers where catalyzed chemical reactions take place and another electron-conductive layers called an anode and a cathode. The mentioned reactions are H 2 →2H + + 2e − (on the layer between the anode and membrane) and O 2 +4H + +4e − →2H 2 O (on the layer between the cathode and membrane). There is vast amount of literature about such fuel cells and their modeling, cf. e.g. [13,16,22] for a survey. Similar scenario applies for methanol or ethanol fuel cells, except a different chemistry on the anode. All the models seem however to be focused on electro-chemistry without taking (thermo)mechanical interactions properly into account, sometimes being rather one-dimensional and mostly not accompanied by any mathematical analysis. Of course, the following outlined model can serve only as an ansatz to which a lot of concrete data is to be supplied.
The generalization of the above presented model to m diffusive constituents consists in taking the concentration c and the (now electro-)chemical potential µ vector valued. Moreover, we consider a vector of electric charges z of the m constituents; some components of z can be zero. Further, we consider the vector of chemicalreactions rate r = r(x, c), electrostatic potential φ of the self-induced electric field, an electric permittivity ǫ = ǫ(x), and d = d(x) concentration of dopands. The x-dependence allows for distinction of particular layers composing the mentioned fuel cells. The above outlined chemistry was only an example -some alternative similar application might be e.g. methanol fuel cells, cf. [8] for a simplified model.
The mass balance (18c) together with (13e) augments to . c − div M(ε(u), χ, c, ϑ)∇µ = r(c), (57a) where φ is to solve the (rest of Maxwell system for) electrostatics balancing the electrical induction ǫ∇φ as − div ǫ∇φ = z·c + d (57c) with some (here unspecified) boundary conditions. In (57a), M is now a 4-th order symmetric tensor and thus [div(M∇µ)] i := 3 k=1 ∂ ∂x k 3 l=1 m j=1 M ijkl ∂µ j ∂x l . The essence of this electro-statical augmentation is that the test of (57a) by µ now relies on the modification of (21) if integrated over Ω as follows: together with (for simplicity unspecified) term arising from boundary conditions for (57c). The energy balance (22) now involves also the energy of the electrostatic field 1 2 Ω ǫ|∇φ| 2 dx. The term r(c)·µ = r(c)·∂ c ϕ 1 (χ, c) + r(c)·zφ arising by the mentioned test of (57a) is to be treated by Gronwall's inequality under some growth qualification on the chemical-reaction rates r. The convergence analysis imitates [23] as far as the electrostatic part concerns. Standard modeling of electro-chemical devices with sizes substantially larger than the so-called Debye length like fuel cells however simplifies the model by considering local electroneutrality, arising as an asymptotics for ǫ → 0 in (57c), cf. e.g. [13].
Final note is that, considering m = 2 and forgetting ε, χ, and ϑ, the system (57) itself represents the classical Roosbroeck's drift-diffusion model for semiconductors [28], the components of c being then interpreted as concentrations of electrons and holes; cf. e.g. [24,Sect. 12.4]. The generalization presented in this section can thus be also interpreted in its very special case m = 2 as a model for thermodynamics of elastic semiconductors, especially if the mobility tensor M would be allowed to depend also on the intensity of the electric field ∇φ.