On the Existence of a Weak Solution of a Half-Cell Model for PEM Fuel Cells

Abstract

A nonlinear boundary value problem (BVP) from the modelling of the transport phenomena in the cathode catalyst layer of a one-dimensional half-cell single-phase model for proton exchange membrane (PEM) fuel cells, derived from the 3D model of Zhou and Liu (2000, 2001), is studied. It is a BVP for a system of three coupled ordinary differential equations of second order. Schauder's fixed point theorem is applied to show the existence of a solution in the Sobolev space 𝐻¹.

1. Introduction

The modelling of fuel cells has been an attractive topic in the field of electrochemical theory. In the last decade, models for proton exchange membrane (PEM) fuel cells have been formulated by many scientists (see, e.g., [1]). Among these models, some complicated systems of partial differential equations (PDEs) were constructed from principles of fluid mechanics, electrostatics, and heat transfers; however, most of them were solved by numerical simulations only. We are interested in the mathematical analysis of the system of differential equations and the discussion is restricted on the transport phenomenon of a single-phase model given by [2]. The more complicated two-phase models, like those mentioned in [1, 3], are not in the scope of this paper.

In [4], by reducing space variables to one dimension and making several assumptions, a system of PDEs in [5] was simplified to a boundary value problem (BVP) for a linear system of decoupled ordinary differential equations (ODEs), and an exact solution was constructed. In [6], a 1D half-cell model reduced from [5] is considered; that model is a BVP for a nonlinear system of three ODEs of second order which are no longer decoupled and it seems to be hard to find an exact solution. By Schaefer's fixed point theorem, the study in [6] is able to show the existence of a solution in the space of continuously twice differentiable functions. In this paper, motivated by [4, 6], we will derive a 1D half-cell model from the 3D model of [2]; it is still a BVP for a nonlinear system of three ODEs of second order; however, the nonlinearity is different from that of [6] and an alternative strategy will be applied; namely, a weak formulation of the BVP will be considered. In this weak formulation, the function space is replaced by the Sobolev space 𝐻¹ and an iteration process associated with Schauder's fixed point theorem will be adopted. The result of this paper indicates a direction of attacking the complicated system of PDEs for the modelling of PEM fuel cells.

Now, we briefly describe the contents of this paper. In Section 2, the governed equations and boundary conditions in the cathode catalyst layer for the 1D half-cell model of PEM fuel cells are derived. In Section 3, the weak form of a linear generalized Neumann problem is described. Existence and uniqueness of the generalized Neumann problem is guaranteed by the Lax-Milgram theorem and it will be shown that the solution for the linear problem has an a priori bound. In Section 4, Schauder's fixed point theorem is applied to prove the existence of an 𝐻¹ solution for the nonlinear system of ODEs.

2. The Model

In this section, we will reduce a 3D model of Zhou and Liu [2, 7] to a 1D half-cell model. This 3D model was a modification of the 2D model given by Gurau et al. [5], so the derivation of the 1D model is quite the same with what we did in [6], we describe the derivation here for the reader's convenience.

Recall (e.g., see [8]) the species equations are0=𝜀𝜌𝐷^eff_𝑘∇²𝑌_𝑘+0,channelanddiffusionlayers,𝜀𝜌𝑆_𝑘,catalystlayer, (2.1) where 𝑌_𝑘 is the concentration of 𝑘th component gas mixture, and 𝐷^eff_𝑘 is the effective diffusivity of the 𝑘th component in the gas mixture, which is given by𝐷^eff_𝑘=𝐷_𝑘𝐷,channel,_𝑘𝜀^1.5,porousmedia. (2.2) At the cathode, the mass generation source terms𝑆_𝑘 for oxygen, water, and protons are 𝑗_𝑐/(2𝐹𝑐),−𝑗_𝑐/(2𝐹𝑐), and 𝑗_𝑐/(𝐹𝑐), respectively. At the anode, the source terms for hydrogen molecules and protons are −𝑗_𝑎/(2𝐹𝑐) and 𝑗_𝑎/𝐹𝑐, where 𝑗_𝑎 and 𝑗_𝑐 are the transfer current density at anode and cathode, which represent the reaction rates. Note that the value of 𝑗_𝑐 is negative and 𝑗_𝑎 is positive. The relationships between 𝑗_𝑎, 𝑗_𝑐 and the species concentration (𝑌_𝐻2 and 𝑌_𝑂2) are given by the Butler-Volmer equations𝑗_𝑎=(𝑎𝑖^ref₀)_𝑎(𝑌_𝐻2𝑌^ref_𝐻2)^1/2𝛼[exp(^𝑎𝐹𝜂𝑅𝑇^𝑎)−exp(−(1−𝛼^𝑎)𝐹𝜂𝑅𝑇^𝑎𝑗)],_𝑐=(𝑎𝑖^ref₀)_𝑐(𝑌_𝑂2𝑌^ref_𝑂2𝛼)[exp(^𝑐𝐹𝜂𝑅𝑇^𝑐)−exp(−(1−𝛼^𝑐)𝐹𝜂𝑅𝑇^𝑐)], (2.3) where 𝑎 is the active catalyst surface area per unit volume of the catalyst layer, 𝑖^ref₀ is the exchange current density under the reference conditions, 𝑇 is the absolute temperature, 𝑅 is the universal gas constant, 𝛼^𝑎 and 𝛼^𝑐 are symmetric factors, and 𝜂_𝑎 and 𝜂_𝑐 are the corresponding overpotentials.

The energy equation is0=𝑘_eff∇²𝑇+0,channelsanddiffusionlayers,𝑄,catalystlayerandmembrane, (2.4) with 𝑘_eff=𝑘_𝑔,channels,−2𝑘_𝑠+1𝜀/(2𝑘_𝑠+𝑘_𝑔)+(1−𝜀)/3𝑘_𝑠,porousmedia, (2.5) where 𝑘_𝑔 is the thermal conductivity of the gas while 𝑘_𝑠 is the thermal conductivity of the solid matrix of the porous media. The heat generation rates 𝑄 in different regions are given by𝑗𝑄=^𝑎⋅𝜂_𝑎+𝑖²𝜎_cl𝑖,anodecatalystlayer,₂𝜎_𝑚𝑗,membrane_𝑐⋅𝜂_𝑐+𝑖²𝜎_cl,cathodecatalystlayer. (2.6) Note that this is a main difference between the 3D model of [2] and the 2D model of [5].

The phase potential satisfies𝑗∇⋅(𝜎∇𝜙)=^𝑐𝑗,cathodecatalystlayer,0,membrane,_𝑎,anodecatalystlayer, (2.7) where 𝜙 is the phase potential, and 𝜎 is the ionic/electric conductivity which depends on 𝑇:𝜂𝜎(𝑇)=^𝑐𝑘_eff𝑙_cl1(0.005139𝜆−0.00326)exp[1268(−1303𝑇)]. (2.8) The current density is given by𝑖=−𝜎∇𝜙. (2.9)

Next, we assume that 𝑇,Φ,𝑌_𝑘 depend on one space variable and restrict to the cathode side of the catalyst layer, following [4], only one species (the oxygen) (i.e., 𝑘=1, and let 𝑌₁=𝑌.) For simplicity, the derivative with respect to 𝑥 is denoted by ^′=𝑑/𝑑𝑥.

From (2.4), the equation for energy becomes𝑇^′′−𝑘(𝑇)𝑌+𝜆𝑓(𝑇)(Φ^′)²=0,𝑥∈(𝑎,𝑏), (2.10) where Φ is the catalyst layer phase potential, 𝑇 is the energy, and 𝑌 is the oxygen mass fraction; 𝑓(𝑇)∈𝐶¹_𝑏(ℝ) is a regularization of 1/𝜆⋅(𝜎(𝑇)/𝑘_eff) away from 𝑇=0 so that 𝑓≥𝛿₁>0 is required, and 𝑘(𝑇)∈𝐶¹_𝑏(ℝ) is a regularization of−𝜂_𝑐𝑘_eff(𝑎𝑖^ref₀)_𝑐(1𝑌^ref_𝑂2𝛼)[exp(^𝑐𝐹𝜂𝑅𝑇^𝑐)−exp(−(1−𝛼^𝑐)𝐹𝜂𝑅𝑇^𝑐)]. (2.11)

By (2.7), in the cathode catalyst layer, we have the following equation for the phase potential(𝑓(𝑇)Φ^′)^′+𝑔(𝑇)𝑌=0,𝑥∈(𝑎,𝑏), (2.12) where 𝑔(𝑇)∈𝐶¹_𝑏(ℝ) is a regularization of−𝜂_𝑐𝑗_𝑐𝑘_eff𝑙_cl𝑌𝜂=−^𝑐𝑘_eff𝑙_cl(𝑎𝑖^ref₀)_𝑐(1𝑌^ref_𝑂2𝛼)[exp(^𝑐𝐹𝜂𝑅𝑇^𝑐)−exp(−(1−𝛼^𝑐)𝐹𝜂𝑅𝑇^𝑐)]. (2.13) And for the oxygen mass fraction, via (2.2), we obtain the equation in the cathode catalyst layer:𝑌^′′−ℎ(𝑇)𝑌=0,𝑥∈(𝑎,𝑏), (2.14) where ℎ(𝑇)∈𝐶¹_𝑏(ℝ) is a regularization of−𝑗_𝑐2𝐹𝑐𝐷^eff_𝑘𝑌=−1𝐷^eff_𝑘⋅12𝐹𝑐⋅(𝑎𝑖^ref₀)_𝑐(1𝑌^ref_𝑂2𝛼)[exp(^𝑐𝐹𝜂𝑅𝑇^𝑐)−exp(−(1−𝛼^𝑐)𝐹𝜂𝑅𝑇^𝑐)]. (2.15) We can assume that ℎ≥𝛿₂>0. The boundary conditions for this 1D model are𝜇₁𝑇(𝑎)−𝜇₂𝑇^′(𝑎)=1,𝑇^′(𝑏)=0,Φ(𝑏)+𝛽𝑓(𝑇(𝑏))⋅Φ^′(𝑏)=0,Φ^′𝛼(𝑎)=0,₁𝑌(𝑎)−𝛼₂𝑌^′(𝑎)=1,𝑌^′(𝑏)=0, (2.16) where 𝜇₁,𝜇₂,𝛼_1,𝛼₂,𝛽>0. It is convenient to let 𝑙=𝑏−𝑎 in the following discussions.

Note that the derivation of these boundary conditions can be found in [6]; therefore we do not repeat here.

Now, we formulate a weak form of the boundary value problem (2.10)∼(2.16).

Let Ω=(𝑎,𝑏) and consider (𝑇,Φ,𝑌)∈(𝐻¹(Ω))³; thus it is a weak solution of (2.10)∼(2.16) if the following equations hold:𝜇−(¹𝑇(𝑎)−1𝜇₂)⋅𝜑(𝑎)−^𝑏_𝑎𝑇^′⋅𝜑^′𝑑𝑥=^𝑏_𝑎𝑘(𝑇)𝑌𝜑𝑑𝑥−𝜆^𝑏_𝑎𝑓(𝑇)(Φ^′)²𝜑𝑑𝑥,∀𝜑∈𝐻¹−1(Ω),𝛽Φ(𝑏)⋅𝜑(𝑏)−^𝑏_𝑎𝑓(𝑇)(Φ^′)⋅𝜑^′𝑑𝑥+^𝑏_𝑎𝑔(𝑇)𝑌𝜑𝑑𝑥=0,∀𝜑∈𝐻¹𝛼(Ω),−(¹𝑌(𝑎)−1𝛼₂)⋅𝜑(𝑎)−^𝑏_𝑎𝑌^′⋅𝜑^′𝑑𝑥=^𝑏_𝑎ℎ(𝑇)𝑌𝜑𝑑𝑥,∀𝜑∈𝐻¹(Ω). (2.17) For (2.17), we have the following existence theorem.

Theorem 2.1. There exists at least one solution (𝑇,Φ,𝑌) of (2.17) in (𝐻¹(Ω))³.

3. Linear Results

Before we prove Theorem 2.1, some linear results should be proved and we still use the notation (𝑇,Φ,𝑌) for the solution of the following (weak) linear generalized Neumann problem:𝛼−(¹𝑌(𝑎)−1𝛼₂)⋅𝜑(𝑎)−^𝑏_𝑎𝑌^′⋅𝜑^′𝑑𝑥=^𝑏_𝑎ℎ(𝑇_∗)𝑌𝜑𝑑𝑥,∀𝜑∈𝐻¹−1(Ω),𝛽Φ(𝑏)⋅𝜑(𝑏)−^𝑏_𝑎𝑓(𝑇_∗)(Φ^′)⋅𝜑^′𝑑𝑥+^𝑏_𝑎𝑔(𝑇_∗)𝑌_∗𝜑𝑑𝑥=0,∀𝜑∈𝐻¹𝜇(Ω),−(¹𝑇(𝑎)−1𝜇₂)⋅𝜑(𝑎)−^𝑏_𝑎𝑇^′⋅𝜑^′𝑑𝑥=^𝑏_𝑎𝑘(𝑇_∗)𝑌_∗𝜑𝑑𝑥−𝜆^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝜑𝑑𝑥,∀𝜑∈𝐻¹(Ω), (3.1)(3.2)(3.3) where 𝑇_∗, 𝑌_∗, Φ_∗∈𝐻¹(Ω). Since the equations for (𝑇,Φ,𝑌) are decoupled, they can be treated separately.

The existence and uniqueness of the solution for the linear generalized Neumann problem (3.1)∼(3.3) is guaranteed by the following Lax-Milgram Theorem (see [9]).

Theorem 3.1 (A theorem on linear monotone operators). Let 𝒜∶𝑋→𝑋^∗ be a linear continuous operator on the real Hilbert space 𝑋. Suppose that 𝒜 is strongly monotone, that is, there is a 𝑐>0 such that ⟨𝒜𝑢,𝑢⟩≥𝑐‖𝑢‖²∀𝑢∈𝑋, (3.4) then for each given 𝑏̃∈𝑋^∗, the operator equation 𝒜𝑢=𝑏̃,𝑢∈𝑋, (3.5) has a unique solution.

Next, we show that the solution for the linear problem has an a priori bound which can be shown to be independent of (𝑇_∗,Φ_∗,𝑌_∗) so that a domain for the iteration process exists.

Theorem 3.2. Suppose that (𝑌,Φ,𝑇) is a weak solution for (3.1)∼(3.3), then one has ‖𝑌‖_𝐻¹_(Ω)≤𝑁₁,‖Φ‖_𝐻¹_(Ω)≤𝑁₂‖𝑌_∗‖_∞,‖𝑇‖_𝐻¹_(Ω)≤𝑁₃‖Φ^′_∗‖²₂+𝑁₄‖𝑌_∗‖_∞+𝑁₅, (3.6)(3.7)(3.8) where 𝑁₁,𝑁₂,𝑁₃,𝑁₄,𝑁₅ are positive constants, and they depend on ‖𝑓‖_∞,‖𝑔‖_∞,‖ℎ‖_∞.

Proof. (1^∘) Equation (3.6) holds. Since 𝛼−(¹𝑌(𝑎)−1𝛼₂)⋅𝜑(𝑎)−^𝑏_𝑎𝑌^′⋅𝜑^′𝑑𝑥=^𝑏_𝑎ℎ(𝑇_∗)𝑌𝜑𝑑𝑥,∀𝜑∈𝐻¹(Ω), (3.9) let 𝜑=𝑌∈𝐻¹(Ω); therefore we have −𝛼₁𝛼₂𝑌²1(𝑎)+𝛼₂𝑌(𝑎)=^𝑏_𝑎|𝑌^′|²𝑑𝑥+^𝑏_𝑎ℎ(𝑇_∗)𝑌²𝑑𝑥≥min(𝛿₂,1)‖𝑌‖²_𝐻1. (3.10)
By (3.10), we can get that‖𝑌‖²_𝐻1≤1min(𝛿₂𝛼,1)(−¹𝛼₂𝑌²1(𝑎)+𝛼₂1𝑌(𝑎))≤min(𝛿₂⋅1,1)4𝛼₁𝛼₂. (3.11) Set 𝑁₁=((1/min(𝛿₂,1))⋅(1/4𝛼₁𝛼₂))^1/2, then (3.6) is proved.
(2^∘) Equation (3.7) holds. From   (3.2),−1𝛽Φ(𝑏)⋅𝜑(𝑏)−^𝑏_𝑎𝑓(𝑇_∗)(Φ^′)⋅𝜑^′𝑑𝑥+^𝑏_𝑎𝑔(𝑇_∗)𝑌_∗𝜑𝑑𝑥=0,∀𝜑∈𝐻¹(Ω). (3.12) Thus let 𝜑=Φ∈𝐻¹(Ω) so that −1𝛽Φ²(𝑏)−^𝑏_𝑎𝑓(𝑇_∗)(Φ^′)²𝑑𝑥+^𝑏_𝑎𝑔(𝑇_∗)𝑌_∗Φ𝑑𝑥=0. (3.13) It follows that ^𝑏_𝑎𝑓(𝑇_∗)(Φ^′)²1𝑑𝑥=−𝛽Φ²(𝑏)+^𝑏_𝑎𝑔(𝑇_∗)𝑌_∗Φ𝑑𝑥≤‖𝑔‖_∞‖𝑌_∗‖_∞‖Φ‖₂⋅𝑙^1/2. (3.14) Since 0<𝛿₁≤𝑓, we have 𝛿₁^𝑏_𝑎(Φ^′)²𝑑𝑥≤‖𝑔‖_∞‖𝑌_∗‖_∞‖Φ‖₂⋅𝑙^1/2. (3.15)
To prove (3.7) we need a lemma (see [10]).
Lemma 3.3. Let 𝑋 denote a real Banach space, and let 𝑢∈𝑊^1,𝑝([0,𝜏];𝑋) for some 1≤𝑝≤∞. Then (i)𝑢∈𝐶([0,𝜏];𝑋) (after possibly being redefined on a set of measure zero), and(ii)𝑢(𝑡)=𝑢(𝑠)+∫^𝑡_𝑠𝑢^′(𝑥)𝑑𝑥forall0≤𝑠≤𝑡≤𝜏.
By Lemma 3.3, we get‖Φ‖²_𝐻1=^𝑏_𝑎|Φ(𝑥)|²𝑑𝑥+^𝑏_𝑎|Φ^′(𝑥)|²𝑑𝑥≤^𝑏_𝑎[^𝑏_𝑥|Φ^′(𝑡)|𝑑𝑡+|Φ(𝑏)|]²𝑑𝑥+^𝑏_𝑎|Φ^′(𝑥)|²𝑑𝑥. (3.16) By (3.15), for (3.16), we arrive at ‖Φ‖²_𝐻1≤2^𝑏_𝑎[^𝑏_𝑥|Φ^′(𝑡)|𝑑𝑡]²𝑑𝑥+2^𝑏_𝑎|Φ(𝑏)|²1𝑑𝑥+𝛿₁‖𝑔‖_∞‖𝑌_∗‖_∞‖Φ‖₂⋅𝑙^1/2≤(^𝑏_𝑎|Φ^′(𝑡)|𝑑𝑡)²⋅2𝑙+2(Φ(𝑏))²1𝑙+𝛿₁‖𝑔‖_∞‖𝑌_∗‖_∞‖Φ‖₂⋅𝑙^1/2. (3.17) On the other hand, for (3.2), we take 𝜑=1∈𝐻¹(Ω), so we know that Φ(𝑏)=𝛽^𝑏_𝑎𝑔(𝑇_∗)𝑌_∗𝑑𝑥. (3.18) Substituting (3.18) into (3.17), we arrive at ‖Φ‖²_𝐻1≤2^𝑏_𝑎|Φ^′(𝑡)|²𝑑𝑡⋅𝑙²+2𝛽²(^𝑏_𝑎𝑔(𝑇_∗)𝑌_∗𝑑𝑥)²1⋅𝑙+𝛿₁𝑙^1/2‖𝑔‖_∞‖𝑌_∗‖_∞‖Φ‖₂≤2𝛿₁𝑙^5/2‖𝑔‖_∞‖𝑌_∗‖_∞‖Φ‖_𝐻¹+1𝛿₁𝑙^1/2‖𝑔‖_∞‖𝑌_∗‖_∞‖Φ‖_𝐻¹+2𝛽²𝑙²‖𝑔‖²_∞‖𝑌_∗‖²_∞. (3.19) Set 𝑥=‖Φ‖_𝐻¹, 𝐵=(2/𝛿₁)𝑙^5/2‖𝑔‖_∞‖𝑌_∗‖_∞+(1/𝛿₁)𝑙^1/2‖𝑔‖_∞‖𝑌_∗‖_∞,𝐶=2𝛽²𝑙²‖𝑔‖²_∞‖𝑌_∗‖²_∞ so that we have 𝑥²−𝐵𝑥−𝐶≤0. (3.20) Hence, we get that 𝑥≤𝐵+𝐵²+4𝐶2=12(2𝛿₁𝑙^5/2+1𝛿₁𝑙^1/2+(2𝛿₁𝑙^5/2+1𝛿₁𝑙^1/2)²+8𝛽²𝑙²)‖𝑔‖_∞‖𝑌_∗‖_∞. (3.21) Take 𝑁₂=(1/2)((2/𝛿₁)𝑙^5/2+(1/𝛿₁)𝑙^1/2+((2/𝛿₁)𝑙^5/2+(1/𝛿₁)𝑙^1/2)²+8𝛽²𝑙²)‖𝑔‖_∞, and from (3.21), we get (3.7).
(3^∘) Equation (3.8) holds. For (3.3), take 𝜑=𝑇∈𝐻¹(Ω); thus we have−𝜇₁𝜇₂𝑇²1(𝑎)+𝜇₂𝑇(𝑎)−^𝑏_𝑎(𝑇^′)²𝑑𝑥=^𝑏_𝑎𝑘(𝑇_∗)𝑌_∗𝑇𝑑𝑥−𝜆^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑇𝑑𝑥. (3.22)
Now we introduce the following lemma (see [9]).
Lemma 3.4. Let G be a bounded region in ℝ^𝑁 with 𝑁≥1, and 𝑢∈𝐻¹(𝐺), set ‖𝑢‖_𝐻¹=(_𝐺(𝑢²+Σ^𝑁_𝑗=1(𝐷_𝑗𝑢)²)𝑑𝑥)^1/2,‖𝑢‖^∗_𝐻1=(_𝐺Σ^𝑁_𝑖=1(𝐷_𝑖𝑢)²𝑑𝑥+_𝜕𝐺𝑢²𝑑𝜎)^1/2, (3.23) then the two norms in (3.23) are equivalent on 𝐻¹(𝐺).
Hence we consider the 𝐻¹(𝐺) norm with the type‖𝑢‖^∗_𝐻1=(_𝐺Σ^𝑁_𝑖=1(𝐷_𝑖𝑢)²𝑑𝑥+_𝜕𝐺𝑢²𝑑𝜎)^1/2, (3.24) and estimate by ^𝑏_𝑎(𝑇^′)²𝑑𝑥+𝜀^′𝑇²(𝑏)+𝜀^′𝑇²(𝑎)(𝜀^′𝜇isasmallpositivenumbertobedetermined)=−¹𝜇₂𝑇²1(𝑎)+𝜇₂𝑇(𝑎)+𝜆^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑇𝑑𝑥−^𝑏_𝑎𝑘(𝑇_∗)𝑌_∗𝑇𝑑𝑥+𝜀^′𝑇²(𝑏)+𝜀^′𝑇²𝜇(𝑎)≤−¹𝜇₂𝑇²1(𝑎)+𝜇₂𝑇(𝑎)+𝜆^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑇𝑑𝑥+(𝜀^𝑏_𝑎𝑇²(𝑥)𝑑𝑥+𝐶_𝜀^𝑏_𝑎|𝑘|²|𝑌_∗|²𝑑𝑥)+𝜀^′𝑇²(𝑏)+𝜀^′𝑇²𝜇(𝑎)≤−¹𝜇₂𝑇²1(𝑎)+𝜇₂𝜀𝑇(𝑎)+𝜆(2‖𝑇‖²_∞+1(2𝜀^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑑𝑥)²)+𝜀‖𝑇‖²_𝐻1+𝐶_𝜀‖𝑘‖²_∞‖𝑌_∗‖²_∞⋅𝑙+𝜀^′𝑇²(𝑏)+𝜀^′𝑇²(𝑎). (3.25) By Lemma 3.3, we have that |𝑇(𝑥)|≤^𝑥_𝑎|𝑇^′(𝑡)|𝑑𝑡+|𝑇(𝑎)|,∀𝑥∈[𝑎,𝑏]. (3.26) Substituting (3.26) into the bound (3.25), we arrive at ^𝑏_𝑎(𝑇^′)²𝑑𝑥+𝜀^′𝑇²(𝑏)+𝜀^′𝑇²𝜇(𝑎)≤−¹𝜇₂𝑇²1(𝑎)+𝜇₂𝜀𝑇(𝑎)+𝜆(2‖𝑇‖²_∞+12𝜀(^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑑𝑥)²)+𝜀‖𝑇‖²_𝐻1+𝐶_𝜀‖𝑘‖²_∞‖𝑌_∗‖²_∞⋅𝑙+𝜀^′(2(^𝑏_𝑎|𝑇^′|𝑑𝑡)²+2|𝑇(𝑎)|²)+𝜀^′𝑇²(𝑎)(𝜀≪𝜀^′)𝜇≤−¹𝜇₂𝑇²1(𝑎)+𝜇₂𝑇(𝑎)+𝜆𝜀2𝑐²₁‖𝑇‖²_𝐻1+𝜆(2𝜀^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑑𝑥)²+𝜀‖𝑇‖²_𝐻1+𝐶_𝜀‖𝑘‖²_∞‖𝑌_∗‖²_∞⋅𝑙+2𝜀^′𝑏_𝑎(𝑇^′(𝑡))²𝑑𝑡⋅𝑙+2𝜀^′𝑇²(𝑎)+𝜀^′𝑇²(𝑎). (3.27) By (3.27), we have that (1−2𝜀^′𝑙)^𝑏_𝑎(𝑇^′(𝑡))²𝑑𝑡+𝜀^′𝑇²(𝑏)+𝜀^′𝑇²𝜇(𝑎)≤(−¹𝜇₂+3𝜀^′)𝑇²1(𝑎)+𝜇₂𝑇(𝑎)+𝜆𝜀2𝑐²₁‖𝑇‖²_𝐻1+𝜆2𝜀(^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑑𝑥)²+𝜀‖𝑇‖²_𝐻1+𝐶_𝜀‖𝑘‖²_∞‖𝑌_∗‖²_∞⋅𝑙, (3.28) where 𝑐₁ is the constant that appeared in the Sobolev inequality (see [11] and note that Ω⊂ℝ) ‖𝑢‖_∞≤𝑐₁‖𝑢‖_𝐻¹_(Ω) (3.29) for all 𝑢∈𝐻¹(Ω).
Hence,𝜀^′‖𝑇‖²_𝐻1≤𝐶₁(𝜀^′𝑏_𝑎(𝑇^′(𝑡))²𝑑𝑡+𝜀^′𝑇²(𝑏)+𝜀^′𝑇²(𝑎))≤𝐶₁((1−2𝜀^′𝑙)^𝑏_𝑎(𝑇^′(𝑡))²𝑑𝑡+𝜀^′𝑇²(𝑏)+𝜀^′𝑇²(𝑎))≤𝐶₁𝜇[(−¹𝜇₂+3𝜀^′)𝑇²1(𝑎)+𝜇₂𝑇(𝑎)+𝜆𝜀2𝑐²₁‖𝑇‖²_𝐻1+𝜆(2𝜀^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑑𝑥)²+𝜀‖𝑇‖²_𝐻1+𝐶_𝜀𝑙‖𝑘‖²_∞‖𝑌_∗‖²_∞], (3.30) where 𝐶₁>0, and 𝜀^′ is chosen to satisfy 𝜀^′<min{(1/1+2𝑙),(𝜇₁/3𝜇₂)}.
From (3.30), we arrive at(𝜀^′−𝜆𝜀2𝑐²₁𝐶₁−𝜀𝐶₁)‖𝑇‖²_𝐻1≤𝐶₁𝜇(−¹𝜇₂+3𝜀^′)𝑇²1(𝑎)+𝜇₂𝐶₁𝜆𝑇(𝑎)+𝐶2𝜀¹(^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑑𝑥)²+𝐶_𝜀𝐶₁‖𝑘‖²_∞‖𝑌_∗‖²_∞⋅𝑙. (3.31)
Now we choose 𝜀>0 such that 𝜀^′−(𝜆𝜀/2)𝑐²₁𝐶₁−𝜀𝐶₁>0. Note that, in (3.3), if we choose test function 𝜑=1∈𝐻¹(Ω), then𝜇−(¹𝑇(𝑎)−1𝜇₂)=^𝑏_𝑎𝑘(𝑇_∗)𝑌_∗𝑑𝑥−𝜆^𝑏_𝑎𝑓(𝑇_∗)(Φ^′_∗)²𝑑𝑥. (3.32) Hence, we obtain that 𝑇(𝑎) is bounded by a number depending on ‖𝑌_∗‖_∞, 𝑙, ‖Φ_∗‖_𝐻1(Ω), ‖𝑘‖_∞ and ‖𝑓‖_∞, and is independent of 𝑇_∗. So ‖𝑇‖²_𝐻1≤𝐶₁(𝜀^′−(𝜆𝜀/2)𝑐²₁𝐶₁−𝜀𝐶₁)[𝜆2𝜀‖𝑓‖²_∞‖Φ^′_∗‖⁴₂+𝐶_𝜀‖𝑘‖²_∞‖𝑌_∗‖²_∞1⋅𝑙−−4𝜇₁𝜇₂+12𝜀^′𝜇²₂], (3.33) where −(1/𝜇²₂)/4(−𝜇₁/𝜇₂+3𝜀^′) is the maximum of (−𝜇₁/𝜇₂+3𝜀^′)𝑇²(𝑎)+(1/𝜇₂)𝑇(𝑎). By (3.33), we have that ‖𝑇‖_𝐻¹≤𝐶₂[𝜆2𝜀‖𝑓‖_∞‖Φ^′_∗‖²₂+𝐶_𝜀𝑙‖𝑘‖_∞‖𝑌_∗‖_∞1+(4𝜇₁𝜇₂−12𝜀^′𝜇²₂)^1/2], (3.34) where 𝐶₂=(𝐶₁/(𝜀^′−(𝜆𝜀/2)𝑐²₁𝐶₁−𝜀𝐶₁))^1/2.
Let𝑁₃=𝐶₂𝜆2𝜀‖𝑓‖_∞,𝑁₄=𝐶₂𝐶_𝜀𝑙‖𝑘‖_∞,𝑁₅=𝐶₂(14𝜇₁𝜇₂−12𝜀^′𝜇²₂)^1/2, (3.35) so (3.8) is proved.

4. Proof of Theorem 2.1

Now, we show the proof of Theorem 2.1.

Step 1. Under the assumptions made in Section 3, for each (𝑇_𝑛,Φ_𝑛,𝑌_𝑛)∈(𝐻¹(Ω))³,𝑛∈ℕ∪{0}, we first consider the linear generalized Neumann problem, 𝜇−(¹𝑇_𝑛+1(𝑎)−1𝜇₂)⋅𝜑(𝑎)−^𝑏_𝑎𝑇^′_𝑛+1⋅𝜑^′𝑑𝑥=^𝑏_𝑎𝑘(𝑇_𝑛)𝑌_𝑛𝜑𝑑𝑥−𝜆^𝑏_𝑎𝑓(𝑇_𝑛)(Φ^′_𝑛)²−1𝜑𝑑𝑥,𝛽Φ_𝑛+1(𝑏)⋅𝜑(𝑏)−^𝑏_𝑎𝑓(𝑇_𝑛)(Φ^′_𝑛+1)⋅𝜑^′𝑑𝑥+^𝑏_𝑎𝑔(𝑇_𝑛)𝑌_𝑛𝛼𝜑𝑑𝑥=0,−(¹𝑌_𝑛+1(𝑎)−1𝛼₂)⋅𝜑(𝑎)−^𝑏_𝑎𝑌^′_𝑛+1⋅𝜑^′𝑑𝑥=^𝑏_𝑎ℎ(𝑇_𝑛)𝑌_𝑛+1𝜑𝑑𝑥, (4.1) for all 𝜑∈𝐻¹(Ω).
Denote by 𝑅𝑣≡(𝑌_𝑛+1,Φ_𝑛+1,𝑇_𝑛+1)=(𝑅𝑌_𝑛,𝑅Φ_𝑛,𝑅𝑇_𝑛) the unique solution of problem (4.1), where 𝑣=(𝑌_𝑛,Φ_𝑛,𝑇_𝑛)∈(𝐻¹(Ω))³.

From Theorem 3.2, we know that ‖𝑌_𝑛+1‖_𝐻1(Ω)≤𝑁₁≡𝑀₁,‖Φ_𝑛+1‖_𝐻1(Ω)≤𝑁₂‖𝑌_𝑛‖_∞≤𝑁₂𝑐₁‖𝑌_𝑛‖_𝐻1(Ω)≤𝑁₂𝑐₁𝑀₁≡𝑀₂,‖𝑇_𝑛+1‖_𝐻1(Ω)≤𝑁₃‖Φ^′_n‖²₂+𝑁₄‖𝑌_𝑛‖_∞+𝑁₅≤𝑁₃‖Φ_𝑛‖²_𝐻1+𝑁₄𝑐₁‖𝑌_𝑛‖_𝐻1(Ω)+𝑁₅≤𝑁₃𝑀²₂+𝑁₄𝑐₁𝑀₁+𝑁₅≡𝑀₃. (4.2)

Now we consider the convex set𝑆={(𝑣_1,𝑣₂,𝑣₃)∈(𝐻¹(Ω))³∶‖𝑣₁‖_𝐻1≤𝑀₁,‖𝑣₂‖_𝐻1(Ω)≤𝑀₂,‖𝑣₃‖_𝐻1(Ω)≤𝑀₃}. (4.3) By the estimate (4.2), we know that 𝑅 maps 𝑆 into 𝑆.

Step 2. We show that 𝑅 is continuous on 𝑆, that is, lim_𝑗→∞(‖𝑅𝑇_𝑗−𝑅𝑇‖_𝐻1+‖𝑅Φ_𝑗−𝑅Φ‖_𝐻1+‖𝑅𝑌_𝑗−𝑅𝑌‖_𝐻1)=0 (4.4) if (𝑇_𝑗,Φ_𝑗,𝑌_𝑗)→(𝑇,Φ,𝑌) in 𝑆, as 𝑗→∞.
Consider the equations𝛼−(¹𝑅𝑌_𝑗(𝑎)−1𝛼₂)𝜑(𝑎)−^𝑏_𝑎(𝑅𝑌_𝑗)^′𝜑^′𝑑𝑥−^𝑏_𝑎ℎ(𝑇_𝑗)(𝑅𝑌_𝑗𝛼)𝜑𝑑𝑥=0,−(¹𝑅𝑌(𝑎)−1𝛼₂)𝜑(𝑎)−^𝑏_𝑎(𝑅𝑌)^′𝜑^′𝑑𝑥−^𝑏_𝑎ℎ(𝑇)(𝑅𝑌)𝜑𝑑𝑥=0, (4.5) for all 𝜑∈𝐻¹(Ω). The difference (4.5) gives −𝛼₁𝛼₂(𝑅𝑌_𝑗−𝑅𝑌)(𝑎)⋅𝜑(𝑎)+^𝑏_𝑎(𝑅𝑌−𝑅𝑌_𝑗)^′𝜑^′𝑑𝑥+^𝑏_𝑎(ℎ(𝑇)(𝑅𝑌)−ℎ(𝑇_𝑗)(𝑅𝑌_𝑗))𝜑𝑑𝑥=0. (4.6) By (4.6), we have −𝛼₁𝛼₂(𝑅𝑌_𝑗−𝑅𝑌)(𝑎)⋅𝜑(𝑎)+^𝑏_𝑎(𝑅𝑌−𝑅𝑌_𝑗)^′𝜑^′+𝑑𝑥_𝑏𝑎[ℎ(𝑇_𝑗)(𝑅𝑌−𝑅𝑌_𝑗)+𝑅𝑌(ℎ(𝑇)−ℎ(𝑇_𝑗))]𝜑𝑑𝑥=0. (4.7) Set 𝐷_𝑗=𝑅𝑌−𝑅𝑌_𝑗, and let 𝜑=𝐷_𝑗 in (4.7), then ^𝑏_𝑎(𝐷^′_𝑗)²𝑑𝑥+𝛿₂^𝑏_𝑎𝐷²_𝑗𝑑𝑥≤‖𝑅𝑌‖_∞‖ℎ^′‖_∞‖𝐷_𝑗‖_∞‖𝑇−𝑇_𝑗‖₂𝑙^1/2. (4.8) Since 0<𝛿₂≤ℎ, we know that min(𝛿_2,1)‖𝐷_𝑗‖²_𝐻1≤𝑐₁‖𝑅𝑌‖_∞‖ℎ^′‖_∞‖𝐷_𝑗‖_𝐻1‖𝑇−𝑇_𝑗‖_𝐻1𝑙^1/2. (4.9) Since lim_𝑗→∞‖𝑇_𝑗−𝑇‖_𝐻1=0, hence we have lim_𝑗→∞‖𝑅𝑌_𝑗−𝑅𝑌‖_𝐻1=0. (4.10) The proof of lim_𝑗→∞‖𝑅Φ_𝑗−𝑅Φ‖_𝐻1=0,lim_𝑗→∞‖𝑅𝑇_𝑗−𝑅𝑇‖_𝐻1=0 (4.11) is similar and is omitted. Hence, 𝑅 is continuous on 𝑆.

Before the next step, we first state a regularity theorem, (see [12]). Consider the operator 𝕃𝑢=𝐷_𝑖(𝑎^𝑖𝑗(𝑥)𝐷_𝑗𝑢+𝑏^𝑖(𝑥)𝑢)+𝑐^𝑖(𝑥)𝐷_𝑖𝑢+𝑑(𝑥)𝑢, (4.12) whose coefficients 𝑎^𝑖𝑗,𝑏^𝑖,𝑐^𝑖,𝑑(𝑖,𝑗=1,…,𝑛) are continuous on a domain Ω⊂ℝ^𝑛.

Theorem 4.1. Let 𝑢∈𝐻¹(Ω) be a weak solution of the equation 𝕃𝑢=𝑝 in Ω where 𝕃 is strictly elliptic on Ω, the coefficients 𝑎^𝑖𝑗,𝑏^𝑖(𝑖,𝑗=1,….,𝑛) are uniformly Lipschitz continuous on Ω, the coefficients 𝑐^𝑖(𝑖=1,…,𝑛),𝑑 are essentially bounded on Ω, and the function 𝑝 is in 𝐿²(Ω). Also, assume that 𝜕Ω is of class 𝐶² and that there exists a function 𝜑∈𝐻²(Ω) for which 𝑢−𝜑∈𝐻¹₀(Ω). Then one has also 𝑢∈𝐻²(Ω) and ‖𝑢‖_𝐻²_(Ω)≤𝐶(‖𝑢‖_𝐿²_(Ω)+‖𝑝‖_𝐿2(Ω)+‖𝜑‖_𝐻2(Ω))for𝐶=𝐶(𝑛,𝜆,𝐾,𝜕Ω), (4.13) where 𝐾=max{‖𝑎^𝑖𝑗,𝑏^𝑖‖_∞,‖𝑐^𝑖,𝑑‖_∞}.

Remark 4.2. Theorem 4.1 continues to hold for sufficiently smooth 𝜕Ω with 𝜑∈𝐻²(Ω) if we assume only that the principal coefficients 𝑎^𝑖𝑗 are in 𝐶⁰(Ω).

Step 3. 𝑅(𝑆) is precompact.
By Theorem 4.1, one has the following:‖𝑌_𝑛+1‖_𝐻2(Ω)≤𝐶(‖𝑌_𝑛+1‖_𝐿2(Ω)+‖𝑌^𝜑_𝑛+1‖_𝐻2(Ω)), (4.14) where 𝑌^𝜑_𝑛+1∈𝐻²(Ω) is a function for which 𝑌_𝑛+1−𝑌^𝜑_𝑛+1∈𝐻¹₀(Ω), ‖Φ_𝑛+1‖_𝐻2(Ω)≤𝐶(‖Φ_𝑛+1‖_𝐿2(Ω)+‖𝑔(𝑇_𝑛)𝑌_𝑛‖_𝐿2(Ω)+‖Φ^𝜑_𝑛+1‖_𝐻2(Ω)), (4.15) where Φ^𝜑_𝑛+1∈𝐻²(Ω) is a function for which Φ_𝑛+1−Φ^𝜑_𝑛+1∈𝐻¹₀(Ω), ‖𝑇_𝑛+1‖_𝐻2(Ω)≤𝐶(‖𝑇_𝑛+1‖_𝐿2(Ω)+‖𝑘(𝑇_𝑛)𝑌_𝑛+𝜆𝑓(𝑇_𝑛)(Φ^′_𝑛)²‖_𝐿2(Ω)+‖𝑇^𝜑_𝑛+1‖_𝐻2(Ω)), (4.16) where 𝑇^𝜑_𝑛+1∈𝐻²(Ω) is a function for which 𝑇_𝑛+1−𝑇^𝜑_𝑛+1∈𝐻¹₀(Ω).
For (4.14)∼(4.16), set𝑌^𝜑_𝑛+1𝑌(𝑥)=^𝑛+1(𝑏)−𝑌_𝑛+1(𝑎)𝑏−𝑎(𝑥−𝑎)+𝑌_𝑛+1Φ(𝑎),_𝜑𝑛+1Φ(𝑥)=^𝑛+1(𝑏)−Φ_𝑛+1(𝑎)𝑏−𝑎(𝑥−𝑎)+Φ_𝑛+1𝑇(𝑎),_𝜑𝑛+1𝑇(𝑥)=^𝑛+1(𝑏)−𝑇_𝑛+1(𝑎)𝑏−𝑎(𝑥−𝑎)+𝑇_𝑛+1(𝑎), (4.17) for all 𝑥∈[𝑎,𝑏].
By (4.17), we obtain that|𝑌^𝜑_𝑛+1𝑌(𝑥)|=|^𝑛+1(𝑏)−𝑌_𝑛+1(𝑎)𝑏−𝑎(𝑥−𝑎)+𝑌_𝑛+1≤(𝑎)|_𝑏𝑎|𝑌^′_𝑛+1|𝑑𝑥+|𝑌_𝑛+1(𝑎)|≤(𝑙^1/2+𝑐₁)𝑀₁, (4.18) for all 𝑥∈[𝑎,𝑏].
From (4.14), (4.18) and by a simple calculation, we obtain‖𝑌_𝑛+1‖_𝐻2(Ω)≤𝐶(‖𝑌_𝑛+1‖_𝐿2(Ω)+‖𝑌^𝜑_𝑛+1‖_𝐻2(Ω))≤𝐶(𝑀₁+𝑀₁𝑙^1/2+𝑐²₁𝑙+4𝑐²₁𝑙)≡𝛾₁. (4.19)
Using the same formula, we also obtain that‖Φ_𝑛+1‖_𝐻2(Ω)≤𝐶(‖Φ_𝑛+1‖_𝐿2(Ω)+‖𝑔(𝑇_𝑛)𝑌_𝑛‖_𝐿2(Ω)+‖Φ^𝜑_𝑛+1‖_𝐻2(Ω))≤𝐶(𝑀₂+‖𝑔‖_∞𝑀₁+𝑀₂𝑙^1/2+𝑐²₁𝑙+4𝑐²₁𝑙)≡𝛾₂,‖𝑇_𝑛+1‖_𝐻2(Ω)≤𝐶(‖𝑇_𝑛+1‖_𝐿2(Ω)+‖𝑘(𝑇_𝑛)𝑌_𝑛+𝜆𝑓(𝑇_𝑛)(Φ^′_𝑛)²‖_𝐿2(Ω)+‖𝑇^𝜑_𝑛+1‖_𝐻2(Ω))≤𝐶(𝑀₃+‖𝑘‖_∞𝑀₁+𝜆‖𝑓‖_∞(𝑐⁴₁𝛾⁴₂𝑙)^1/2+𝑀₃(𝑙^1/2+𝑐₁)²𝑙+4𝑐²₁𝑙)≡𝛾₃. (4.20) Thus, we have ‖𝑅𝑌_𝑛‖_𝐻2(Ω)≤𝛾₁,‖𝑅Φ_𝑛‖_𝐻2(Ω)≤𝛾₂,‖𝑅𝑇_𝑛‖_𝐻2(Ω)≤𝛾₃,∀𝑛∈ℕ. (4.21)
So, 𝑅 maps 𝑆 into a bounded set in (𝐻²(Ω))³, since 𝐻²(Ω) is compactly imbedded in 𝐻¹(Ω) (see, e.g., [9]); hence, 𝑅(𝑆) is precompact in (𝐻¹(Ω))³.

Hence by Schauder's fixed point theorem, 𝑅 has a fixed point and there exists a (𝑇,Φ,𝑌)∈(𝐻¹(Ω))³ satisfyng (2.17). Thus, we complete the proof of Theorem 2.1.