*Article* **Conversion Kinetics and Ionic Conductivity in Na-**β**"-Alumina + YSZ (Na**β**"AY) Sodium Solid Electrolyte via Vapor Phase Conversion Process**

**Liangzhu Zhu 1,2,\* and Anil V. Virkar 1,\***


**Abstract:** Sodium ion batteries have been receiving increasing attention and may see potential revival in the near future, particularly in large-scale grid energy storage coupling with wind and solar power generation, due to the abundant sodium resources, low cost, and sufficiently high energy density. Among the known sodium ion conductors, the Na-β"-alumina electrolyte remains highly attractive because of its high ionic conductivity. This study focuses on the vapor phase synthesis of a Na-β"-Alumina + YSZ (Naβ"AY) composite sodium electrolyte, which has higher mechanical strength and stability than conventional single phase β"-Alumina. The objectives are the measurement of conversion kinetics through a newly developed weight-gain based model and the determination of sodium ionic conductivity in the composite electrolyte. Starting samples contained ~70 vol% α-Alumina and ~30 vol% YSZ (3 mol% Y2O3 stabilized Zirconia) with and without a thin alumina surface layer made by sintering in air at 1600 ◦C. The sintered samples were placed in a powder of Na-β"-alumina and heat-treated at 1250 ◦C for various periods. Sample dimensions and weight were measured as a function of heat treatment time. The conversion of α-Alumina in the α-Alumina + YSZ composite into Naβ"AY occurred by coupled diffusion of sodium ions through Na-β"-alumina and of oxygen ions through YSZ, effectively diffusing Na2O. From the analysis of the time dependence of sample mass and dimensions, the effective diffusion coefficient of Na2O through the sample, *Deff* , was estimated to be 1.74 <sup>×</sup> <sup>10</sup>−<sup>7</sup> cm<sup>2</sup> <sup>s</sup>−1, and the effective interface transfer parameter, *keff* , was estimated as 2.33 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm s<sup>−</sup>1. By depositing a thin alumina coating layer on top of the bulk composite, the chemical diffusion coefficient of oxygen through single phase Na-β"-alumina was estimated as 4.35 <sup>×</sup> <sup>10</sup>−<sup>10</sup> cm<sup>2</sup> <sup>s</sup><sup>−</sup>1. An AC impedance measurement was performed on a fully converted Naβ"AY composite, and the conductivity of the composite electrolyte was 1.3 <sup>×</sup> <sup>10</sup>−<sup>1</sup> S cm−<sup>1</sup> at 300 ◦C and 1.6 <sup>×</sup> <sup>10</sup>−<sup>3</sup> S cm−<sup>1</sup> at 25 ◦C, indicating promising applications in solid state or molten salt batteries at low to intermediate temperatures.

**Keywords:** sodium β"-alumina; Naβ"AY; sodium electrolyte; sodium solid-state battery; vapor phase process; sodium batteries

### **1. Introduction**

Na-β"-alumina (more commonly referred to as β"-alumina) is a sodium ion conductor, which is used as a solid electrolyte in sodium-sulfur batteries, sodium-nickel chloride batteries, and in alkali metal-based thermoelectric converters [1–3]. While significant attention to Naβ"-alumina dates back as early as in the 1960s when sodium-sulfur batteries were originally developed by Weber and Kummer at the Ford Motor Company [1], the relatively high working temperature in a molten sodium-sulfur battery has significantly limited their applications and development, and soon, they were almost completely replaced by more popular lithium ion batteries in mobile applications. However, with increasing demands on large-scale energy

**Citation:** Zhu, L.; Virkar, A.V. Conversion Kinetics and Ionic Conductivity in Na-β"-Alumina + YSZ (Naβ"AY) Sodium Solid Electrolyte via Vapor Phase Conversion Process. *Membranes* **2022**, *12*, 567. https://doi.org/10.3390/ membranes12060567

Academic Editor: Jijeesh Ravi Nair

Received: 3 May 2022 Accepted: 26 May 2022 Published: 30 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

storage due to climate change and the bloom of a solid-state battery, applications of Na-β" alumina are seemingly rising again because of the abundance of sodium as compared to lithium sources and the high ionic conductivity of about 1 S cm−<sup>1</sup> in single crystal Na-β" alumina and 0.2–0.4 S cm−<sup>1</sup> in polycrystalline β"-alumina at 300 ◦C [4–10]. For example, recently, Fertig et al. have provided a detailed review on the potential revival of Na-β"-alumina for sodium solid-state batteries [4]. Ligon et al. have reported large planar Na-β"-alumina solid electrolytes (150 mm in diameter) for next generation Na-batteries [5]. Lee et al. have demonstrated a full-scale pilot model of a cost-effective sodium-nickel-iron chloride battery over 40 Ah using a Na-β"-alumina electrolyte [9]. Zhu et al. have reported that by ion exchange, a Na-β"-alumina-containing composite electrolyte may be ion exchanged with molten salts to lithium-ion or silver-ion conducting electrolytes [11], which may further expand the application of Na-β"-alumina in other types of batteries.

The conventional process for the fabrication of dense samples of Na-β"-alumina consists of first calcining a mixture of Na2CO3, α-alumina, and LiNO3 (or MgO) at ~1250 ◦C in air, which leads to a powder mixture containing Na-β"-alumina, Na-β-alumina (βalumina), and some NaAlO2. Powder compacts are then sintered in sealed platinum or MgO crucibles at ~1600 ◦C [12]. Sintering in sealed containers suppresses the loss of Na2O through the vapor phase. Densification occurs by a transient liquid phase mechanism. Sintered samples are subsequently heat-treated at a lower temperature (~1400 ◦C) to convert β-alumina into β"-alumina by reacting with NaAlO2. Usually, a small amount of NaAlO2 remains unreacted in β"-alumina as a thin film along grain boundaries, which makes it susceptible to degradation due to moisture in the atmosphere. Thus, Na-β"-alumina made by the conventional process is stored in desiccators.

In concept, if an already-sintered α-alumina can be converted into Na-β"-alumina by reacting it with Na2O, it may be possible to avoid the formation of NaAlO2 along the grain boundaries, thus making it moisture-resistant. A possible vapor process for the conversion of sintered α-alumina into Na-β"-alumina involves exposing it to a vapor containing Na2O. The approximate composition of Na-β"-alumina is Na2O ~6Al2O3. For the conversion of α-alumina into Na-β"-alumina, the reaction is

$$\text{Na}\_2\text{O} + \text{\textdegree-6Al}\_2\text{O}\_3 \rightarrow \text{Na}\_2\text{O} \cdot \text{\textdegree-6Al}\_2\text{O}\_3 \tag{1}$$

If a fully dense sintered Al2O3 is exposed to Na2O vapor, a thin layer of Na-β"-alumina forms on the surface. Further conversion of the interior α-alumina requires the transport of Na2O (as coupled (ambipolar) transport of 2Na<sup>+</sup> and O2−) through the formed Na-β" alumina (Figure A1a in Appendix A). While the Na+ diffusivity through Na-β"-alumina is high, O2<sup>−</sup> diffusivity is very low. Thus, the conversion kinetics are dictated by the diffusion coefficient of oxygen ions in Na-β"-alumina and are very sluggish. Furthermore, the conversion kinetics are diffusion-limited (parabolic). At 1300 ◦C, for example, the thickness of α-alumina converted into Na-β"-alumina in 16 h was only about 25 μm [13].

In the novel vapor phase process, a two-phase composite of α-alumina and YSZ, with both phases being contiguous, is exposed to Na2O vapor [14]. Once α-alumina on the surface is converted to Na-β"-alumina, subsequent conversion of the interior α-alumina to Na-β"-alumina involves the transport of Na2O, such that 2Na<sup>+</sup> transports through the formed Na-β"-alumina, while the O2<sup>−</sup> transports through the YSZ phase (Figure 1b). More specifically, the effective diffusion of Na2O through the two phase mixture occurs by a coupled (ambipolar) transport of 2Na+ through Na-β"-alumina and O2<sup>−</sup> through YSZ. Since the diffusion coefficient of O2<sup>−</sup> through YSZ is much higher than that of O2<sup>−</sup> through Na-β"-alumina, the kinetics of conversion of α-alumina + YSZ into Na-β"-alumina + YSZ is much faster than the conversion of single phase α-alumina into Na-β"-alumina. The resulting two-phase composite of Naβ"AY is also much stronger than conventional Na-β"-alumina and is moisture-resistant.

**Figure 1.** XRD patterns of α-Alumina + YSZ composites (**a**) before and (**b**) after subjecting to vapor phase conversion at 1250 ◦C.

Since the invention of the novel vapor phase process patented by Virkar et al. [14], it has been gaining attention in Na batteries [15–19]. A recent study has shown that a planar sodium nickel chloride battery demonstrated a specific energy density of as high as 350 Wh kg<sup>−</sup>1, operated at 190 ◦C over 1000 cycles. The composite electrolyte adopted in that battery consisted of α-alumina and 8YSZ, also with a volume ratio of 7:3 as starting materials, and then was converted via the vapor phase conversion process. Such a composite electrolyte also offered high mechanical strength. A flexural strength of higher than 300 MPa has been reported in literature [20,21], making them an excellent Na-conducting membrane in Na batteries.

The only available model reported so far on the kinetics of this novel vapor phase conversion process is the one reported by Parathasarathy and Virkar [22], where the kinetics of conversion of α-alumina + YSZ composites of various grain sizes over a range of temperatures between 1250 ◦C and 1400 ◦C were investigated. The experimental procedure involved packing sintered α-alumina + YSZ samples into Na-β"-alumina powder, heat treating at a given temperature for a period of time, cooling down to room temperature, and grinding/polishing the sample to measure (using a microscope) the conversion thickness *x*(*t*) as a function of cumulative heat treatment time *t*. The process thus is tedious, as it requires repeated grinding and polishing. Furthermore, this involves the destruction of the sample. Under some situations, where both conversion fraction and conductivity as a function of conversion time are of interest, the above thickness-based model may be insufficient, since it requires cutting of the sample. Therefore, an alternate model needs to be developed to fulfill the routine examination of this vapor phase process.

The objective of the present work was to measure the sample weight and the external dimensions after various thermal treatments so that the kinetics of conversion could be measured more effectively without destroying the sample and, at the same time, minimizing the effort required in the actual measurements.

An equally important model was also developed and experimentally verified to measure the chemical diffusion coefficient of oxygen through single phase -β"-alumina by depositing a thin layer of α-alumina on the α-alumina + YSZ composite prior to conversion. The α-alumina-coated samples were subjected to the same conversion treatment as the uncoated samples. Both of the newly developed models may be used as effective tools for future studies and applications of the vapor phase process. Note that the concepts and models developed here may well fit to other binary or possibly ternary diffusion system where diffusion kinetics of certain mobile ions, such as oxygen-ion, proton, metal ions, etc., are of interest.

#### **2. Materials and Methods**

#### *2.1. Preparation of α-Alumina + YSZ Samples*

High purity α-alumina (AKP-53, Sumitomo Chemical) and 3YSZ (TZ-3Y, Tosoh Corporation, Tokyo, Japan) were mixed in a volume ratio of alumina:3YSZ of 7:3. A slurry of the powder mixture was made in distilled water with a small amount of ammonium poly methyl methacrylate (DARVAN C-N, R. T. Vanderbilt, Inc., New York, NY, USA), added as a dispersant and as a stabilizer. The slurry was planetary milled for 12 h, dried, and then heated to 500 ◦C for 5 h. The dried powder was sieved through a 70-mesh screen to remove any large agglomerates. A few grams of powder was placed in a circular die and uniaxially pressed under a force of 5 tons, followed by cold isostatic pressing at 30,000 psi. The disc was then pre-sintered at 1200 ◦C for 2 h and then polished on both sides to a finish of ~50 nm using an alumina suspension in water. After polishing, selected samples were either spin-coated or dip-coated with an alumina suspension depending on the desired thickness. After slowly heating the samples to 1000 ◦C to burnout the organic components, all of the samples were sintered in air at 1600 ◦C for 5 h. The weight and the dimensions (diameter and thickness) were measured as the baseline information.

#### *2.2. Formation of Na-β"-Alumina + YSZ Samples via Vapor Phase Conversion*

Sintered discs were then buried in Na-β"-alumina powder (Materials and Systems Research, Inc., Salt Lake City, UT, USA) in an alumina crucible and covered with a lid. The crucible was placed in a furnace in air and heated to 1250 ◦C (5 ◦C min−<sup>1</sup> ramp up rate), maintained for a period of time at temperature, and rapidly cooled to room temperature (>15 ◦C min−<sup>1</sup> ramp down rate). The objective was to ensure that for most of the reaction time, the samples were at the heat treatment temperature. Since the process of vapor phase conversion of α-alumina + YSZ into Na-β"-alumina + YSZ is thermally activated, most of the conversion occurred isothermally at 1250 ◦C, and any conversion during heating to 1250 ◦C and cooling down from the heat treatment temperature could be neglected. This process was repeated several times until a maximum cumulative time of ~84 h. The dimensions of the samples and their weights were carefully recorded after each thermal treatment. For weight measurement, a balance with a resolution of 0.0001 g and repeatability <0.0002 g was used. The digital caliper used for dimension measurements had a resolution of 0.01 mm with negligible deviation in repeat measurements. In order to examine the conversion front at different conversion times, a small part of one of the samples was cut from the edge after each conversion period and fine polished for microstructure characterization. Samples used for weight measurements were not cut.

#### *2.3. Characterization*

Similar samples were prepared for X-ray diffraction (XRD) and scanning electron microscopy (SEM). An Energy Dispersive Spectroscopy (EDS) elemental line scan was performed on a partially converted sample to examine the α-alumina + YSZ/Na-β"-alumina + YSZ interface. Density measurements were performed by the standard fluid immersion method on an as-sintered α-alumina + YSZ sample and a fully converted Naβ"AY sample. Conductivity was measured in air as a function of temperature using an AC Electrochemical Impedance Spectroscopy (EIS) with gold paste as electrodes. Measurements were conducted over a frequency range from 10 Hz to 1 MHz. The high frequency intercept was taken as the measure of conductivity.

#### **3. Results and Discussion**

#### *3.1. X-ray-Diffraction*

Figure 1 shows XRD patterns of an as-sintered α-alumina + YSZ sample and a fully converted Naβ"AY sample. In Figure 1a, all peaks are identified as belonging to either αalumina or YSZ. Peaks belonging to α-alumina have been identified by arrows. In Figure 1b, which is for the converted sample, peaks belonging to Na-β"-alumina have been identified by arrows. No peaks belonging to α-alumina were observed in the converted sample. Some of the peaks belonging to Na-β-alumina overlapped with those of Na-β"-alumina. Thus, the presence of some Na-β-alumina cannot be ruled out.

#### *3.2. SEM and EDS Analysis*

Figure 2a shows an SEM image of a polished section of a sample that had been heattreated at 1250 ◦C for 4 h by packing in Na-β"-alumina powder. The converted region is darker in color. The average conversion thickness measured was ~110 μm. Figure 2b shows an EDS scan from the surface of the sample, across Na-β"-alumina + YSZ/α-alumina + YSZ interface and into the α-alumina + YSZ unconverted region. While the EDS scan was not smooth and showed considerable scatter, the demarcation between the converted region and the pristine region is clear.

Figure 3 shows SEM micrographs of the cross-sections of an uncoated sample, a sample with a ~2.5 μm-thick layer Na-β"-alumina (initially, an alumina layer of 2 μm-thickness) and a sample with a ~15 μm-thick layer of Na-β"-alumina (initially, an alumina layer of 12 μm-thickness) after conversion treatment for 4 h, 4 h, and 84 h, respectively. The micrographs in (b) and (c) thus correspond to samples that were initially coated with αalumina. These layers fully converted to Na-β"-alumina. The interiors of all three samples contained Na-β"-alumina + YSZ. The dark phase is Na-β"-alumina; the light phase is YSZ. The Na-β"-alumina and YSZ phases were both contiguous, as would be required for coupled transport of Na2O to occur through the two phase mixture: Na+ through Na-β"-alumina and O2<sup>−</sup> through YSZ [22].

Figure 4 compares the conversion thicknesses of the alumina-coated (~2.5 μm Na-β" alumina after conversion) portion of the sample, with part of the region that was uncoated, after conversion at 1250 ◦C for 6 h. The conversion thickness corresponding to the noncoated region was about twice that of the coated region. This shows that the presence of a 2.5-μm Na-β"-alumina surface coating substantially lowered the conversion kinetics.

**Figure 2.** (**a**) An SEM image showing the converted and the pristine regions of the initially α-alumina + YSZ sample after 4 h of conversion at 1250 ◦C. (**b**) The corresponding EDS line scan is indicated by the arrow.

**Figure 3.** SEM images showing the enlarged view at the interface of (**a**) The uncoated region, (**b**) A ~2.5 μm Na-β"-alumina-coated sample, and (**c**) A ~15 μm Na-β"-alumina-coated sample.

**Figure 4.** SEM images showing the conversion fronts for the uncoated and 2.5 μm Na-β"-aluminacoated regions after 6 h of conversion.

#### *3.3. Estimation of the Kinetic Parameters on the Uncoated Samples*

Appendix A gives the relevant equations describing the kinetics of conversion and the chemical diffusion coefficients. Equation (A1) gives the conversion thickness *x* corresponding to isothermal treatment time of *t* in terms of the effective diffusion coefficient, *Deff* , and the effective interface transfer parameter, *keff* . Appendix B gives the kinetic equation in terms of the sample dimensions and the increase in weight as a function of time. The following describes an example of the estimation of the kinetic parameters on an uncoated sample.

The density of the as-sintered α-alumina + YSZ samples, *ρAY*, was measured as 4.59 g cm<sup>−</sup>3. The initial thickness of the α-Al2O3 + YSZ disc sample was 3.64 mm, and the initial diameter was 29.54 mm. Figure 5 shows the measured weight, the thickness, and the diameter changes in percent with respect to the initial measurements. The initial weight of the sample, *mo*, was 11.2452 g. After conversion for a cumulative time of 84 h at 1250 ◦C, the sample weight increased to 11.9302 g. That is, the sample weight increased by 685 mg after a cumulative conversion time of 84 h at 1250 ◦C. After 72 h, the weight increased only slightly. This indicated that the sample had fully or very close to fully converted after 84 h.

**Figure 5.** Measured weight and geometry change of an uncoated sample as a function of conversion time.

At 1250 ◦C. The sample cross-section confirmed that it had fully converted to Naβ"-alumina + YSZ. Thus, the final weight of the sample of 11.9302 g was taken as the weight after full conversion, that is, *m*(∞) = 11.9302 g. The density of the fully converted thinner sample prepared using the same procedure, *ρβ <sup>Y</sup>*, was measured as 3.94 g cm−3. From the measured values of *mo*, *<sup>m</sup>*(∞), *<sup>ρ</sup>AY*, and *ρβ <sup>Y</sup>*, the ratio *<sup>x</sup>* (*t*) *<sup>x</sup>*(*t*) corresponding to Equation (A16) was measured as

$$\frac{\mathbf{x}'(t)}{\mathbf{x}(t)} = \frac{l\_0}{l(\infty)} = \frac{\rho\_{AY} m\_o}{\rho\_{\mathbb{R}^{\theta'}Y} m(\infty)} = 0.81\tag{2}$$

The ratio *mo*

$$\frac{m\_0}{m(\infty)} = 0.94\tag{3}$$

is used in Equation (A18).

The thickness measured after conversion (84 h) was 4.45 mm. Thus, *lo <sup>l</sup>*(∞) based on the thickness measurement is

$$\frac{l\_o}{l(\infty)} = 0.82\tag{4}$$

The increase in thickness was ~22%. The corresponding increase in diameter was only ~2.5%. For this reason, the assumption that most of the dimensional change occurs along the thickness direction is reasonable. Based on the weight gain, it is concluded that the stoichiometry of Na-β"-alumina is almost exactly Na2O·6Al2O3.

The kinetic equation in terms of weight change given in Equation (A19) is reproduced here

$$\frac{\left(\Delta m(t)\right)^{2}}{4D\_{eff}\left(A\rho\_{\beta''Y}\left(1-\frac{m\_a}{m(\infty)}\right)\right)^{2}} + \frac{\Delta m(t)}{2k\_{eff}\left(A\rho\_{\beta''Y}\left(1-\frac{m\_a}{m(\infty)}\right)\right)} = t\tag{5}$$

The same equation can be written as

$$\Delta m(t) = 4D\_{eff} \left( A \rho\_{\beta''Y} \left( 1 - \frac{m\_{\mathcal{O}}}{m(\infty)} \right) \right)^2 \frac{t}{\Delta m(t)} - \frac{2D\_{eff}}{k\_{eff}} \left( A \rho\_{\beta''Y} \left( 1 - \frac{m\_{\mathcal{O}}}{m(\infty)} \right) \right) \tag{6}$$

In Equations (5) and (6), the Δ*m*(*t*) is in g. Equation (6) shows that a plot of Δ*m*(*t*) vs. *t* <sup>Δ</sup>*m*(*t*) should be linear, with the slope given by 4*Deff* - *Aρβ <sup>Y</sup>* - <sup>1</sup> <sup>−</sup> *mo m*(∞) <sup>2</sup> and the intercept given by <sup>−</sup>2*Deff keff* - *Aρβ <sup>Y</sup>* - <sup>1</sup> <sup>−</sup> *mo m*(∞) . Thus, from the slope and the intercept, one should be able to estimate *Deff* and *keff* . Figure 6a gives a plot of <sup>Δ</sup>*m*(*t*) *<sup>A</sup>* vs. *<sup>t</sup>* - Δ*m*(*t*) *A* . Thus, the slope is 4*Deff* - *ρβ <sup>Y</sup>* - <sup>1</sup> <sup>−</sup> *mo m*(∞) <sup>2</sup> , and the intercept is <sup>−</sup>2*Deff keff* - *ρβ <sup>Y</sup>* - <sup>1</sup> <sup>−</sup> *mo m*(∞) . As seen in Figure 6a, the plot is linear. The corresponding slope was 4.99 × <sup>10</sup>−<sup>8</sup> g2 cm−<sup>4</sup> <sup>s</sup><sup>−</sup>1, and the intercept was −0.04 g cm<sup>−</sup>2. It is significant that the intercept was negative, as required. The corresponding estimated values of the kinetic parameters are *Deff* ∼= 1.74 × <sup>10</sup>−<sup>7</sup> cm2 <sup>s</sup>−<sup>1</sup> and *keff* ∼= 2.33×10−<sup>6</sup> cm s<sup>−</sup>1. The experimental data were directly fitted to Equation (5) by a polynomial fitting, as shown in Figure 6b. The estimated values were *Deff* ∼= 1.76 × <sup>10</sup>−<sup>7</sup> cm2 <sup>s</sup>−<sup>1</sup> and *keff* ∼= 2.29 × <sup>10</sup>−<sup>6</sup> cm s−1, which were in good agreement with the results of the linear plot in Figure 6a. The estimated value of *Deff* was in good agreement with that measured by Parthasarrathy and Virkar [22] based on conversion thickness, which was *Deff* ∼= 1.5 × <sup>10</sup>−<sup>7</sup> cm2 <sup>s</sup><sup>−</sup>1. The *keff* measured by Parthasarathy and Virkar [22] ranged between ~3.6 × <sup>10</sup>−<sup>7</sup> cm s−<sup>1</sup> for large-grained samples to ~1.0 × <sup>10</sup>−<sup>6</sup> cm s−<sup>1</sup> for fine-grained samples. The *keff* estimated in the present work was thus between 3- and 10-times larger than in the study by Parthasarathy and Virkar. Given the possible differences in microstructures, the agreement is deemed reasonable. The present work thus shows that weight measurements can be used to estimate both kinetic parameters without having to section and polish after each thermal treatment.

**Figure 6.** (**a**) A plot of <sup>Δ</sup>*m*(*t*) *<sup>A</sup>* vs. *<sup>t</sup>* - Δ*m*(*t*) *A* for the uncoated sample. From the slope and the intercept,

both *Deff* and *keff* can be determined. (**b**) A plot of *<sup>t</sup>* vs. <sup>Δ</sup>*m*(*t*) *<sup>A</sup>* for the uncoated sample by a second order polynomial fitting with the intercept set as 0. From the 1st order and 2nd order coefficients, both *keff* and *Deff* can be determined.

#### *3.4. Effects of Na-β"-Alumina Coating on Conversion Kinetics*

Figure 7 compares the measured weight changes for the uncoated, 2.5 μm coated, and 15 μm coated (thicknesses correspond to the formed Na-β"-alumina after conversion) disc samples. For the uncoated and 15 μm coated samples, the figure shows the actual measured weight changes (in percent) as a function of time. For the sample with a 2.5 μm coating, the actual conversion thickness was measured as a function of time. In order to graph the data on the same plot, the expected weight changes in the 2.5 μm coated sample were calculated from the measured conversion thickness using Equation (A15). In Figure 7, these data are shown by a dashed line. The uncoated sample showed a much higher weight percentage change than the coated samples. As seen in the figure, the uncoated sample exhibited the largest weight gain (faster kinetics), and the 15 μm coated sample showed the lowest (slowest kinetics of the three samples tested). Also significant is the observation that over the duration of the tests, the weight gain vs. time plots were linear for both of the coated samples, indicating that the kinetics can be described as interface-controlled.

**Figure 7.** Comparison in weight change for (**a**) uncoated, (**b**) ~2.5 μm Na-β"-alumina-coated, and (**c**) ~15 μm Na-β"-alumina-coated samples as a function of conversion time at 1250 ◦C.

#### *3.5. Estimation of the Kinetic Parameter for the Coated Samples*

Figure 8a,b show plots of <sup>Δ</sup>*mt <sup>A</sup>* vs. *t* for the two coated samples. Note that these plots are linear, indicating that the kinetics is interface-controlled. Equation (A19) then becomes

$$\frac{\Delta m(t)}{A} \stackrel{\sim}{=} 2k\_{eff} \left( \rho\_{\beta''Y} \left( 1 - \frac{m\_o}{m(\infty)} \right) \right) t \tag{7}$$

**Figure 8.** (**a**) A plot of <sup>Δ</sup>*m*(*t*) *<sup>A</sup>* vs. *t* for a sample with a 2.5 μm Na-β"-alumina surface coating. From the slope, the *keff* is obtained. (**b**) A plot of <sup>Δ</sup>*m*(*t*) *<sup>A</sup>* vs. *t* for a sample with a 15 μm Na-β"-alumina surface coating. From the slope, the *keff* is obtained.

Thus, the slope is given by

$$2k\_{eff}\left(\rho\_{\beta''Y}\left(1-\frac{m\_o}{m(\infty)}\right)\right) \tag{8}$$

from which the *keff* can be estimated. The estimated values of *keff* were 6.45 × <sup>10</sup>−<sup>7</sup> cm s−<sup>1</sup> for the sample with a 2.5 <sup>μ</sup>m Na-β"-alumina coating and 2.26 × <sup>10</sup>−<sup>7</sup> cm s−<sup>1</sup> for the sample with a 15 μm Na-β"-alumina coating. This *keff* has a number of series contributions, as discussed below.

Figure 9 shows a schematic variation of the chemical potential of Na2O from the gas phase, *μ<sup>I</sup> Na*2*O*; just inside the single phase Na-β"-alumina close to the gas phase, *<sup>μ</sup>*<sup>1</sup> *Na*2*O*; in the single phase Na-β"-alumina close to the single phase Na-β"-alumina/Na-β"-alumina + YSZ interface, *μ*<sup>2</sup> *Na*2*O*; just inside the Na-β"-alumina + YSZ close to the single phase Naβ"-alumina/Na-β"-alumina + YSZ interface, *μ*<sup>3</sup> *Na*2*O*; inside the Na-β"-alumina + YSZ close to the Na-β"-alumina + YSZ/α-alumina + YSZ interface, *μ*<sup>4</sup> *Na*2*O*; and inside the α-alumina + YSZ close to the Na-β"-alumina + YSZ/α-alumina + YSZ interface, *μI I Na*2*O*. The transfer of Na2O across the three interfaces can be described by three interface kinetic parameters, namely, *kI*1, *k*23, and *k*4*I I*. The transport of Na2O through single-phase Na-β"-alumina is dictated by the chemical diffusion coefficient of Na2O, *<sup>D</sup>β Na*2*O*. However, this transport occurs through a layer of fixed thickness, *δ*. Thus, it reflects as an interface step, with the

corresponding *<sup>k</sup>*<sup>12</sup> given by *<sup>D</sup>β Na*2*O <sup>δ</sup>* . The measured interface parameter, *keff* , is thus related to these parameters by the following equation.

$$\frac{1}{k\_{eff}} = \frac{1}{k\_{I1}} + \frac{1}{k\_{23}} + \frac{1}{k\_{4II}} + \frac{\delta}{\hat{D}\_{NazO}^{\theta''}}\tag{9}$$

**Figure 9.** A schematic showing the variation of the chemical potential of Na2O, *μNa*2*O*, through a sample with surface layer of Na-β"-alumina of thickness, *δ*.

The preceding suggests that a plot of <sup>1</sup> *keff* vs. layer thickness, *<sup>δ</sup>*, should be a straight line with the slope given by <sup>1</sup> *Dβ Na*2*O* and the intercept given by <sup>1</sup> *kI*<sup>1</sup> <sup>+</sup> <sup>1</sup> *<sup>k</sup>*<sup>23</sup> <sup>+</sup> <sup>1</sup> *<sup>k</sup>*4*I I* . Figure 10 is a plot of <sup>1</sup> *keff* vs. *<sup>δ</sup>*. In this study, measurements were conducted on samples with only two different thicknesses. Thus, there are only two data points. The straight line shown is thus through these two points. Furthermore, plotted on the same plot, however, is <sup>1</sup> *keff* for the uncoated sample, for which *δ* = 0. For the uncoated sample, the relevant equation is

$$\frac{1}{k\_{eff}} = \frac{1}{k\_{I3}} + \frac{1}{k\_{4II}}\tag{10}$$

where *kI*<sup>3</sup> corresponds to the interface transfer parameter at the gas phase/Na-β"-alumina + YSZ interface.

**Figure 10.** A plot of 1/*keff* vs. the Naβ"-alumina layer thickness, *δ*. Also shown in the figure is the *keff* for the uncoated sample.

Note that

$$\frac{1}{k\_{I3}} + \frac{1}{k\_{4II}} \neq \frac{1}{k\_{I1}} + \frac{1}{k\_{23}} + \frac{1}{k\_{4II}}\tag{11}$$

From Figure 10, it is observed that the

$$\frac{1}{k\_{I3}} + \frac{1}{k\_{4II}} < \frac{1}{k\_{I1}} + \frac{1}{k\_{23}} + \frac{1}{k\_{4II}}\tag{12}$$

or 
$$\frac{1}{k\_{I3}} < \frac{1}{k\_{I1}} + \frac{1}{k\_{23}}\tag{13}$$

or

$$k\_{I3} > \frac{k\_{I1}k\_{23}}{k\_{I1} + k\_{23}}\tag{14}$$

At this stage, it is not possible to separately determine the different interface transfer parameters.

From the slope, the chemical diffusion coefficient of Na2O through single phase Na<sup>β</sup>"-alumina is given by *<sup>D</sup>β Na*2*<sup>O</sup>* <sup>∼</sup><sup>=</sup> 4.35 <sup>×</sup> <sup>10</sup>−<sup>10</sup> cm<sup>2</sup> <sup>s</sup>−1. By contrast, the *Deff* determined for coupled transport of Na2O as 2Na<sup>+</sup> through Na-β"-alumina and O2<sup>−</sup> through YSZ was ~1.74 <sup>×</sup> <sup>10</sup>−<sup>7</sup> cm2 <sup>s</sup><sup>−</sup>1. That is, *<sup>D</sup>β* /*YSZ Na*2*<sup>O</sup>* was over 400 times larger than *<sup>D</sup>β Na*2*O*. As given in Equation (A6), *Deff* for the two-phase sample is given by

$$D\_{eff} = \breve{D}\_{Na\_2O}^{\otimes \prime \prime / \chi\_{SZ}} \frac{2V\_m}{f} \Delta C\_{Na\_2O}^{\otimes \prime} \tag{15}$$

wherein

$$
\hat{D}\_{Na\_2O}^{\mathcal{S}''/YSZ} = \frac{\mathcal{C}\_{O^{2-}}^{YSZ} D\_{O^{2-}}^{YSZ} V\_{YSZ}}{\mathcal{C}\_{Na\_2O}^{\mathcal{S}''}} \frac{1}{k\_B T} \left(\frac{\partial \mu\_{Na\_2O}}{\partial \ln \mathcal{C}\_{Na\_2O}^{\mathcal{S}''}}\right) \tag{16}
$$

That is, the *<sup>D</sup>β* /*YSZ Na*2*<sup>O</sup>* is proportional to the diffusion coefficient of O2<sup>−</sup> in YSZ, namely, *DYSZ <sup>O</sup>*2<sup>−</sup> . However, in single-phase Na-β"-alumina, the Na2O chemical diffusion coefficient is given by

$$\bar{D}\_{Na\_2O}^{\mathbb{R}^{\otimes^{\prime\prime}}} = \mathbb{C}\_{O^{2-}}^{\mathbb{R}^{\otimes^{\prime}}} D\_{O^{2-}}^{\mathbb{R}^{\prime\prime}} \frac{1}{k\_B T} \left( \frac{\partial \mu\_{Na\_2O}^{\mathbb{R}^{\otimes^{\prime}}}}{\partial \mathbb{C}\_{Na\_2O}^{\mathbb{R}^{\otimes^{\prime}}}} \right) \tag{17}$$

That is, the *<sup>D</sup>β Na*2*<sup>O</sup>* is proportional to the diffusion coefficient of O2<sup>−</sup> in Na-β"-alumina. It is well known that oxygen diffusion is much faster in YSZ, an oxygen ion conductor, than in Na-β"-alumina. The present results are consistent with this expectation. Table 1 lists the measured *keff* and *Deff* for alumina-coated and non-coated samples.

**Table 1.** Summary of measured *keff* and *Deff* for alumina-coated and non-coated samples.


#### *3.6. Measurement of Ionic Conductivity*

Conductivity was measured on a disc sample of 0.9-mm-thickness that had been converted at 1250 ◦C for 36 h and a disc sample of 3.7-mm-thickness that had been converted at 1250 ◦C for 108 h. The heat treatment was long enough to ensure that the entire samples converted into Naβ"AY. Figure 11 shows an Arrhenius plots for the Naβ"AY composite and along with some common Na and Li electrolytes for better comparison [23–29]. Over the range of temperatures measured, the oxygen ion conductivity of YSZ was much lower than the sodium ion conductivity of Na-β"-alumina. Thus, the measured conductivity is attributed exclusively to sodium ion conduction. The measured activation energies of ~23.2 kJ mol−<sup>1</sup> (after 36 h conversion) and ~21.2 k J mol−<sup>1</sup> (after 108 h conversion) were also in accordance with the previously reported activation energy of Na+-ion conduction, especially in fine-grained Na-β"-alumina [10]. At 300 ◦C, the resistivity was measured as ~16 Ω cm (conductivity of 6.3 × <sup>10</sup>−<sup>2</sup> S cm<sup>−</sup>1) and ~8 <sup>Ω</sup> cm (conductivity of 1.3 × <sup>10</sup>−<sup>1</sup> S cm<sup>−</sup>1) for the samples converted for 36 h and 108 h, respectively. The calculated conductivity based on the best linear fitting of the Arrhenius plot at 25 ◦C was 1.4 × <sup>10</sup>−<sup>3</sup> and 1.6 × <sup>10</sup>−<sup>3</sup> S cm−<sup>1</sup> for the samples converted for 36 h and 108 h, respectively.

**Figure 11.** An Arrhenius plot of the measured conductivity of NaβAY and some common Na and Li electrolytes. Conductivity sources: Na3PS4 [23], 60Na2S-40GeS2 [23], β/β"-alumina mixture [15], (1 − x)Li2S − xP2S5 [23], Li3xLa(2/3−x)TiO3 (LLTO, x = 0.17) [27], NaSiCON (Produced by Ceramatec) [23], Li(1+x)Ti(2−x)Alx(PO4)3 (LATP, x = 0.3 [29], LiPF6 (1M in EC/DMC mixture) [26], Liβ"AY (Vapor phase) [25], Li(4−2x)ZnxGeS4 (LiSICON, x = 0.05) [28].

An interesting feature revealed by this combined plot is that the activation energies of the listed electrolytes do not change vastly regardless of the crystal structure or conduction mechanisms. In comparison with a liquid electrolyte, most of the solid electrolytes show better conductivity at a higher temperature regime. The comparison in the solid electrolytes indicates that the Naβ"AY composite electrolyte showed higher conductivity than most of the listed electrolytes, other than the mixed β/β"-alumina electrolyte and the NaSICON-type Na electrolyte. However, when mechanical strength is of a concern, Naβ"AY (~300 MPa flexural strength [20,21]) may be a better choice, due to a higher mechanical strength than that of NaSICON electrolyte, with ~100 MPa flexural strength [30]. Therefore, the Na"βAY composite produced by the vapor phase conversion process seems

a promising solid electrolyte for Na battery applications. It is to be noted that although the Li-β"-alumina + YSZ electrolyte (Liβ"AY) synthesized by the vapor phase conversion process showed lower conductivity as compared with the Naβ"AY electrolyte, its conductivity was still higher than a few of the listed Li solid electrolytes.

#### **4. Conclusions**

The vapor phase conversion of α-alumina + YSZ into Na-β" alumina + YSZ (Naβ"AY) was investigated using the measurement of weight gain and external sample dimensions as a function of time. Vapor phase conversion studies were also conducted on samples coated with thin layers of α-alumina on the one surface, which during the process, first converted into a layer of Na-β"-alumina. In the coated samples, the transport of Na2O occurred through single phase Na-β"-alumina prior to the conversion of the interior α-alumina in the α-alumina + YSZ composite. From the measurement of conversion kinetics of the uncoated sample, the effective diffusion coefficient, *Deff* , and effective interface transfer parameter, *keff* , were determined. At 1250 ◦C, the measured *Deff* was 1.74 × <sup>10</sup>−<sup>7</sup> cm<sup>2</sup> <sup>s</sup><sup>−</sup>1, in agreement with the previously reported results of Parthasarathy and Virkar [22]. The measured *keff* was 2.33 × <sup>10</sup>−<sup>6</sup> cm s<sup>−</sup>1, which is within an order of magnitude of the values reported by Parthasarathy and Virkar [22]. Some difference may be related to the possible differences in microstructural details. The *Deff* is a measure of the chemical diffusion coefficient of Na2O in the two-phase Naβ"AY composite, *<sup>D</sup>β* /*YSZ Na*2*<sup>O</sup>* , which is proportional to the oxygen ion diffusion in YSZ.

The measured kinetics of conversion of the coated samples were found to be linear in time, indicating interface-controlled kinetics. Since the chemical diffusion of Na2O must first occur through the surface Na-β"-alumina layer, it reflects as an interface step which is proportional to the oxygen diffusion coefficient through single-phase Na-β"-alumina. From the dependence of conversion kinetics on the thickness of the surface Na-β"-alumina layer, the chemical diffusion coefficient of Na2O through Na-β"-alumina, *<sup>D</sup>β Na*2*O*, was estimated. The estimated value was 4.35 × <sup>10</sup>−<sup>10</sup> cm<sup>2</sup> <sup>s</sup><sup>−</sup>1, which is over 400-times smaller than *<sup>D</sup>β* /*YSZ Na*2*<sup>O</sup>* . The measured *<sup>D</sup>β Na*2*<sup>O</sup>* was proportional to the oxygen diffusion coefficient through Na-β"-alumina. The present work thus provides a method to determine the chemical diffusion coefficient of Na2O through a Naβ"AY composite, as well as through single-phase Na-β"-alumina. This method may have applicability to other systems.

The measured highest conductivity of the converted Naβ"AY was 1.6 × <sup>10</sup>−<sup>3</sup> S cm−<sup>1</sup> at <sup>25</sup> ◦C and 1.3 × <sup>10</sup>−<sup>1</sup> S cm−<sup>1</sup> at 300 ◦C. These values are higher than many single-phase Li or Na solid electrolytes, indicating the Naβ"AY composite synthesized by the vapor phase conversion process may be used as a promising electrolyte membrane for Na batteries.

**Author Contributions:** Methodology and original draft preparation, L.Z.; Conceptualization, formal analysis, writing-review and editing, L.Z. and A.V.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the National Science Foundation, under Grant No. DMR-1407048, and in part by the US Department of Energy, under Grant Number DE-FG02- 06ER46086.

**Data Availability Statement:** The data that support the plots within this paper are available from the corresponding authors upon reasonable request.

**Acknowledgments:** This work was supported in part by the National Science Foundation, under Grant No. DMR-1407048, and in part by the US Department of Energy, under Grant Number DE-FG02-06ER46086. The authors would like to thank the support from Fuel Cell Research Team and Key Laboratory of Advanced Fuel Cells and Electrolyzers Technology of Zhejiang Province, Ningbo Institute of Materials Technology and Engineering, Chinese Academy of Sciences; Ningbo Municipal Government and the Chinese Academy of Science (Beijing).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Vapor Phase Conversion Kinetics and Transport Parameters**

The vapor phase conversion of α-Alumina + YSZ (or with a thin α-alumina coating layer) to Na-β"-alumina/YSZ consists of three steps in series. (a) Step 1: Na2O from the vapor phase is adsorbed on the Na-β" alumina + YSZ (or pure Na-β"-alumina for an alumina-coated sample) surface. (b) Step 2: coupled diffusion of 2Na<sup>+</sup> and O2<sup>−</sup> through converted Na-β"-alumina + YSZ in the composite or through, first, a thin layer of Na-β" alumina single phase, followed by coupled diffusion through a Na-β" alumina + YSZ composite for the coated sample. (c) Step 3: conversion of α-Alumina + YSZ to Na-β" alumina + YSZ composite at the reaction front; this requires lateral transport of O2<sup>−</sup> ions arriving from the YSZ phase through Na-β"-alumina, α-Alumina, or along the α-Alumina/Na-β"-alumina interface, as discussed by Parthasarathy and Virkar [22]. Since in both coated and non-coated samples, the reaction front is the same, the interfacial kinetic parameter at this interface should be the same. This transport at the reaction front is similar to the conversion of γ-Fe into pearlite (α-Fe + Fe3C) in low carbon steels [31].

The kinetics of the conversion process of α-alumina + YSZ into Na-β"-alumina + YSZ was modeled by Parathasarthy and Virkar [22]. Figure A1b shows a schematic. Isothermal conversion kinetics could be described by an equation of the form

$$\frac{\mathbf{x}^2}{D\_{eff}} + \frac{\mathbf{x}}{k\_{eff}} = \mathbf{t} \tag{A1}$$

where *x* is the thickness converted in time *t*, *Deff* is the effective diffusion coefficient of Na2O for ambipolar transport of Na2O through the composite with 2Na+ transporting through Na-β"-alumina and O2<sup>−</sup> transporting through the YSZ, and *keff* is the interface transfer parameter that includes processes at the vapor phase/Na-β"-alumina + YSZ interface and also at the Na-β"-alumina + YSZ/α-alumina + YSZ reaction front. The *Deff* is given by

$$D\_{eff} = \frac{2V\_{YSZ}C\_{O^{2-}}^{YSZ}D\_{O^{2-}}^{YSZ} \left(\mu\_{NazO}^{I} - \mu\_{NazO}^{II}\right)V\_m}{k\_BTf} \tag{A2}$$

where *Vm* is the molar volume of Na-β"-alumina, *μ<sup>I</sup> Na*2*<sup>O</sup>* is the chemical potential of Na2O in gas phase at the Na2O-gas phase/Na-β"-alumina interface, *μI I Na*2*<sup>O</sup>* is the chemical potential of Na2O in α-alumina at the Na-β"-alumina + YSZ/α-alumina + YSZ interface, and *f* is a factor which includes the dimensional change that occurs when α-alumina in the composite converts into Na-β"-alumina. We will write Equation (A2) as

$$D\_{eff} = \frac{2V\_{YSZ}\mathbb{C}\_{O^{2-}}^{YSZ}D\_{O^{2-}}^{YSZ}\left(\mu\_{Na\_2O}^{I} - \mu\_{Na\_2O}^{II}\right)\left(\mathbb{C}\_{Na\_2O}^{\theta''I} - \mathbb{C}\_{Na\_2O}^{\theta''II}\right)V\_m}{k\_B T \left(\mathbb{C}\_{Na\_2O}^{\theta''I} - \mathbb{C}\_{Na\_2O}^{\theta''II}\right)f} \tag{A3}$$

or

$$D\_{eff} = \frac{2V\_{YSZ}\mathcal{C}\_{O^{2-}}^{YSZ}D\_{O^{2-}}^{YSZ}\left(\mathcal{C}\_{Na\_2O}^{\theta''I} - \mathcal{C}\_{Na\_2O}^{\theta''II}\right)V\_m}{k\_B Tf} \left(\frac{\partial \mu\_{Na\_2O}^{\theta''}}{\partial \mathcal{C}\_{Na\_2O}^{\theta''}}\right) \tag{A4}$$

or

$$D\_{eff} = \frac{2V\_{YSZ}\mathcal{C}\_{O^{2-}}^{YSZ}D\_{O^{2-}}^{YSZ}}{k\_B T} \left(\frac{\partial \mu\_{Na\_2O}^{\beta''}}{\partial \mathcal{C}\_{Na\_2O}^{\beta''}}\right) \frac{V\_m}{f} \Delta \mathcal{C}\_{Na\_2O}^{\beta''} \tag{A5}$$

where Δ*Cβ Na*2*<sup>O</sup>* denotes the difference in Na2O concentration in Na-β"-alumina from the Na2O-gas phase/Na-β"-alumina boundary to the Na-β"-alumina/α-alumina boundary. A comparison of Equations (A2) and (A5) shows that

$$D\_{eff} = \breve{D}\_{Na\_2O}^{\otimes \prime \prime / YSZ} \frac{2V\_m}{f} \Delta C\_{Na\_2O}^{\otimes \prime} \tag{A6}$$

**Figure A1.** Schematics showing vapor phase conversion in (**a**) single phase α-alumina, (**b**) α-alumina + YSZ composite, and (**c**) α-alumina-coated α-alumina + YSZ composite. The schematics do not show grain structure in this one-dimensional model. The widths of YSZ, α-alumina, and Na-β"-alumina are essentially equal to the respective volume fractions.

Thus, the estimated *Deff* measured from conversion kinetics is a measure of the chemical diffusion coefficient of Na2O through the two-phase Na-β"-alumina + YSZ composite, namely, *<sup>D</sup>β* /*YSZ Na*2*<sup>O</sup>* , since <sup>2</sup>*Vm <sup>f</sup>* <sup>Δ</sup>*Cβ Na*2*<sup>O</sup>* is a dimensionless constant. This chemical diffusion coefficient is proportional to the O2<sup>−</sup> diffusion coefficient through YSZ.

The *keff* also depends on the chemical potentials of Na2O. Equation (A1) assumes that the conversion thickness at the beginning of the isothermal process is negligible.

Figure A1a shows a schematic of the conversion of single phase α-alumina into Naβ"-alumina by the vapor phase process. In this case, the effective diffusion coefficient, as determined from conversion kinetics, would given by

$$D\_{eff} = \tilde{D}\_{Na2O}^{\beta''} \frac{2V\_{\text{nr}}}{f'} \Delta C\_{Na2O}^{\beta''} \tag{A7}$$

where *f* is the factor that includes the dimensional change when α-alumina converts into Na-β"-alumina. As given in Equation (A7), the effective diffusion coefficient is proportional to the chemical diffusion coefficient of Na2O through Na-β"-alumina, which is effectively proportional to the O2<sup>−</sup> diffusion coefficient through single phase Na-β"-alumina. This is expected to be very low compared to the O2<sup>−</sup> diffusion coefficient through YSZ. Thus, we expect that

$$
\mathbb{D}\_{Na\_2O}^{\beta^{\nu}/YSZ} > > \mathbb{D}\_{Na\_2O}^{\beta^{\nu}} \tag{A8}
$$

Figure A1c schematically shows the conversion process in thin α-alumina-coated α-alumina + YSZ composites. Initially, the thin α-alumina layer is converted into Na-β" alumina, the kinetics of which is governed by the effective diffusion coefficient, corresponding to transport through single-phase Na-β"-alumina, which is governed by oxygen diffusion through single-phase Na-β"-alumina. Once this layer is fully converted into Na-β"-alumina, the conversion of α-alumina in the interior α-alumina + YSZ composite begins. This conversion of the interior α-alumina into Na-β"-alumina occurs by coupled transport of 2Na+ through Na-β"-alumina and of O2<sup>−</sup> through YSZ. However, this conversion requires a continual supply of Na2O, which must transport through the thin surface layer of the single-phase Na-β"-alumina of a fixed thickness. Because of the fixed thickness of this layer, the overall conversion kinetics contain an additional 'interface' transfer step, corresponding to Na2O diffusion through the single-phase Na-β"-alumina of a fixed thickness, *δ*. That is, the transport of Na2O through the single-phase Na-β"-alumina layer is equivalent to an additional interface step given by

$$k\_{\beta''} \propto \frac{D\_{O^{2-}}^{\beta''}}{\delta} \tag{A9}$$

The larger the layer thickness, the smaller is the *kβ* , and the slower is the kinetics. This suggests that the kinetics of conversion should be slower in samples coated with a Na-β"-alumina layer and that by comparing the kinetics of conversion of samples initially coated with two different thicknesses of alumina layers, it should be possible to estimate the diffusion coefficient of O2- through a single phase Na-β"-alumina.

#### **Appendix B. Study of Kinetics of Vapor Phase Conversion by Weight Measurement**

Figure A2 shows a schematic of the conversion process of an α-alumina + YSZ composite. The following analysis assumes that the samples are in the form of thin discs. When α-alumina is converted into Na-β"-alumina, there is both an increase in weight, since Na2O reacts with Al2O3 to form Na-β"-alumina (Na2O~6Al2O3), and an increase in volume, since the density of Na-β"-alumina is lower than that of Al2O3, and it also includes Na2O, which is an additional reason for the increase in volume. Prior work has shown that if the disks are sufficiently thin, most of the dimensional change occurs in the thin direction, and the disk diameter remains essentially unchanged. If the initial sample thickness is *lo*, the cross-sectional area is *A*; then, the initial sample volume is *Vo* = *Alo*. The density of the α-alumina + YSZ composite is *ρAY*, and the density of the Na-β"-alumina is *ρβ <sup>Y</sup>*. The initial mass of the sample is

$$m\_o = \rho\_{AY} V\_o = \rho\_{AY} A l\_o \tag{A10}$$

**Figure A2.** Schematics showing (**a**) a sintered α-alumina + YSZ sample before conversion and (**b**) after partial conversion. *x* (*t*) is the initial thickness of α-alumina + YSZ which converted into Na-β"-alumina + YSZ of thickness *x*(*t*), in a conversion time of *t*.

After time *t*, the thickness *x* (*t*) of the initial sample on both surfaces converts to Na-β"-alumina + YSZ of thickness *x*(*t*). Thus, the sample thickness at time *t* becomes

$$l(t) = l\_\circ + 2\mathbf{x}(t) - 2\mathbf{x}'(t) \tag{A11}$$

and the volume at time *t* is

$$V(t) = Al(t) = A\left(l\_o + 2x(t) - 2x'(t)\right) \tag{A12}$$

The corresponding mass at time *t* is

$$\rho\_i(\mathbf{t}) = A \left\{ \rho\_{AY} (l\_\vartheta - 2\mathbf{x}'(t)) + \rho\_{\tilde{\beta}''} \mathbf{2} \mathbf{x}(t) \right\} = \rho\_{AY} A l\_\vartheta + 2A \left( \mathbf{x}(t) \rho\_{\tilde{\beta}''} - \mathbf{x}'(t) \rho\_{AY} \right) \tag{A13}$$

or

$$m(t) = m\_o + \Delta m(t)\tag{A14}$$

where

$$
\Delta m(t) = 2A \left( \mathbf{x}(t)\rho\_{\beta''Y} - \mathbf{x}'(t)\rho\_{AY} \right) \tag{A15}
$$

is the increase in sample weight after time *t*. Note that

$$\frac{\mathbf{x}'(t)}{\mathbf{x}(t)} = \frac{l\_o}{l(\infty)} = \frac{m\_o \rho\_{\beta^\nu Y}}{m(\infty)\rho\_{AY}}\tag{A16}$$

where the final thickness after full conversion is *l*(∞), and the final mass is *m*(∞). Thus,

$$
\Delta m(t) = 2Ax(t)\rho\_{\beta''Y}\left(1 - \frac{m\_o}{m(\infty)}\right) = 2Ax(t)\rho\_{\beta''Y}\left(1 - \frac{\rho\_{AY}l\_o}{\rho\_{\beta''Y}l(\infty)}\right) \tag{A17}
$$

The *ρAY*, *lo*, and *mo* are measured on the initial sample. The *ρβ <sup>Y</sup>* is measured on a fully converted sample. Thus, both *l*(∞) and *m*(∞) can be estimated. The preceding implies that it is not necessary to carry out full conversion to achieve the final weight, *m*(∞), and the final thickness, *l*(∞), in order to use the kinetic equation given in (A1). In the present work, however, the sample was heat-treated for a long enough time to achieve the final weight and thickness. Thus, Equation (A17) is used in the analysis.

From Equation (A17), the conversion thickness in terms of weight gain is given by

$$x(t) = \frac{\Delta m(t)}{2A\rho\_{\beta''Y}\left(1 - \frac{m\_a}{m(\infty)}\right)} = \frac{\Delta m(t)}{2A\rho\_{\beta''Y}\left(1 - \frac{\rho\_{\beta'}I\_\sigma}{\rho\_{\beta''}I(\infty)}\right)}\tag{A18}$$

Equations (A17) and (A18) show that the conversion thickness *x*(*t*) at time *t* can be given in terms of change in weight, Δ*m*(*t*), and the densities of the initial sample, (*ρAY*), and after conversion to Na-β"-alumina, (*ρβ <sup>Y</sup>*), the sample cross sectional area, *A*; the initial and final sample thicknesses, *lo* and *l*(∞), respectively; and the initial and final masses, *mo* and *m*(∞), respectively. Therefore, Equation (A1) can be written as

$$\frac{\left(\Delta m(t)\right)^{2}}{4D\_{eff}\left(A\rho\_{\beta''Y}\left(1-\frac{m\_{\sigma}}{m(\infty)}\right)\right)^{2}} + \frac{\Delta m(t)}{2k\_{eff}\left(A\rho\_{\beta''Y}\left(1-\frac{m\_{\sigma}}{m(\infty)}\right)\right)} = t\tag{A19}$$

Thus, by measuring Δ*m*(*t*) as a function of time *t* and fitting experimental data to Equation (A19), it should be possible to obtain both kinetic parameters, *Deff* and *keff* . Alternatively, Equation (A19) may be linearized as follows

$$
\Delta m(t) = 4A^2 (\rho\_{\beta''Y} - \rho\_{AY}\lambda)^2 D\_{eff} \frac{t}{\Delta m(t)} - 2A(\rho\_{\beta''Y} - \rho\_{AY}\lambda) \frac{D\_{eff}}{k\_{eff}} \tag{A20}
$$

A plot of Δ*m*(*t*) versus *<sup>t</sup>* <sup>Δ</sup>*m*(*t*) is thus expected to be linear. Then, from the slope and the intercept, it should be possible to obtain both *Deff* and *keff* .

#### **References**


### *Article* **Fault Structural Analysis Applied to Proton Exchange Membrane Fuel Cell Water Management Issues**

**Etienne Dijoux 1,2,\*, Nadia Yousfi Steiner 2, Michel Benne 1, Marie-Cécile Péra <sup>2</sup> and Brigitte Grondin-Perez <sup>1</sup>**


**Abstract:** Proton exchange membrane fuel cells are relevant systems for power generation. However, they suffer from a lack of reliability, mainly due to their structural complexity. Indeed, their operation involves electrochemical, thermal, and electrical phenomena that imply a strong coupling, making it harder to maintain nominal operation. This complexity causes several issues for the design of appropriate control, diagnosis, or fault-tolerant control strategies. It is therefore mandatory to understand the fuel cell structure for a relevant design of these kinds of strategies. This paper proposes a fuel cell fault structural analysis approach that leads to the proposition of a structural graph. This graph will then be used to highlight the interactions between the control variables and the functionalities of a fuel cell, and therefore to emphasize how changing a parameter to mitigate a fault can influence the fuel cell state and eventually cause another fault. The final aim of this work is to allow an easier implementation of an efficient and fault-tolerant control strategy on the basis of the proposed graphical representation.

**Keywords:** fuel cell fault structural analysis; diagnosis; fault tolerant control; fuel cell cathode water management

### **1. Introduction**

Proton exchange membrane fuel cells (PEMFCs) are efficient and clean energy supply systems. However, they are subject to the occurrence of various faults, which decreases their reliability. Faults are inordinate phenomena that degrade a system's performance more or less rapidly and substantially [1]. Their occurrence can be attributed to several factors (exogenous and endogenous). Both exogenous factors, such as gas purity or demanding load profile, and endogenous factors, such as poor internal design or natural aging, can lead to fault occurrence and, therefore, to fuel cell damage. The operating conditions need to be adjusted to mitigate the faults. Moreover, PEMFCs are nonlinear, multivariate, and strongly coupled systems, which complicates their ability to be maintained under normal operation. Therefore, it is necessary to highlight the coupling of the PEMFCs' parameters to facilitate the understanding of the fault occurrence process. Indeed, an exhaustive analysis of the variables' effects and interactions inside the system is a major issue to be considered to set an efficient fault-tolerant control (FTC) strategy.

The literature provides different approaches to modeling and analyzing a system in order to understand the influences of its variables and their interactions. Noyes [1] highlighted two types of methodologies that allow this analysis: statistical and functional. On the one hand, a statistical analysis consists in the observation of events [2,3] with the aim of making assumptions to predict events in similar situations. For instance, Bayesian statistical analysis [2] aims to process small datasets. Indeed, this approach allows one to obtain relevant information by limiting costly observations. The obtained information is then iteratively refined according to a Bayesian law. On the other hand, functional analysis,

**Citation:** Dijoux, E.; Steiner, N.Y.; Benne, M.; Péra, M.-C.; Grondin-Perez, B. Fault Structural Analysis Applied to Proton Exchange Membrane Fuel Cell Water Management Issues. *Electrochem* **2021**, *2*, 604–630. https://doi.org/10.3390 /electrochem2040038

Academic Editors: Marc Cretin, Sophie Tingry and Zhenghua Tang

Received: 30 September 2021 Accepted: 19 October 2021 Published: 1 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

focuses on all functions ensured by the system and their influence on the occurrence of a fault. These approaches are usually based on a graphical formalism that leads to the specification of good operating conditions.

Several methods for functional analysis can be found in the literature, such as FAST (Function Analysis System Technique), SADT (Structured Analysis and Design Technique), or FTA (Fault Tree Analysis).

FAST is a graphical representation of a system's functions that answers the following questions: why does a function have to be ensured? How is this function ensured? When must the function be ensured? Other functional analysis methods also exist in the literature [4–6], and they consist in analyzing a system with measurements to observe its structural and functional modifications.

The SADT [7,8] method is a graphical tool associated with a top-down analysis method. This method is used to decompose a function with a functionally oriented methodology. It consists of modeling the process by breaking it down into subsets. A data-based or function-based diagram is then created to model the process of each subset.

FTA is also a functionally based analysis that allows one to perform a failure mode analysis. It is a top-down approach: a top event is considered, and all combinations of sub-events that lead to it are determined [9].

The top event is reached with a combination of several sub-events. Sub-events have their origin in the combination of basic events. Fault trees, therefore, bring information about the variables (basic events) that are involved in a specific fault occurrence.

To summarize, on the one hand, FAST uses a diagram to organize the ways of thinking, acting, or talking. It enables the development of technical solutions according to a functional logic. However, this method does not consider the system complexity or the coupling phenomena of internal system variables. On the other hand, the SADT method takes into consideration the complexity and allows the analytical decomposition of the system according to a hierarchical structure. However, this approach does not allow links between transitions of operating conditions, such as from normal condition to a faulty one. On the contrary, FTA is relevant in understanding how a fault can occur with a combination of basic events. Indeed, FTA allows the link between each transition of the operating conditions. With logic gates and sub-events, it describes all paths that lead to the appearance of the fault. FTA highlights the information about fault occurrence and, therefore, the relevant variables for fault mitigation. Fuel cells are strongly coupled systems, and several variables have mutual interactions that are not highlighted by FTA. For this reason, a new analysis methodology, which is called Fault Structural Analysis (FSA), is applied to these systems. Indeed, the FSA allows describing the system's structure and highlighting all variables that influence the fuel cell operating conditions. This approach is also relevant for designing fault tolerant control strategies because it is helpful for fault mitigation and system monitoring process. For instance, as presented in [10], authors present their strategy for fuel cell fault mitigation. Their work consists of gathering information about the PEMFC state of health through the remaining useful lifetime. The approach is based on the analysis of the system nominal and faulty conditions which are provided by a key variable behavior. This strategy is thus highly dependent on the relevant choice of the key variable that should be subject to a study of its field of action in the fuel cell for more efficiency. In [11], Yang et al. try to improve the PEMFC reliability with the implementation of a robust fault observer for air management system fault diagnosis. Once again, the choice of the estimated variable is a key factor for their diagnosis tool. Indeed, the implementation of their strategy depends on sensitivity of the diagnosis variables to the fuel cell functionalities which are subject to faulty conditions. In [12], authors proposed a fuel cell health management system. They used the electrochemical impedance spectroscopy (EIS) in a fault tolerant control strategy in order to diagnose the water management faults. The drawback of EIS diagnosis tools lies in their low computational time, their offline operating mode and the cost of the used equipment. To avoid this problem, a solution for the implementation of a diagnostic tool can be based, for each relevant variable, on identifying the one most influenced by each

faulty condition. Another study proposed by Rubio et al. [13], consist of the implementation of a fuzzy model to determine the water dehydration in a PEMFC. The real-time aspect of the strategy involves the use of fast response time of the control variable. The current, the flow rate and the voltage are thus used in the strategy for the fuel cell hydration characterization. This study only considers fast response time variables for the diagnosis tool, but the studied phenomena have low, medium and high frequency behavior. In the case of the introduction of variables which are influenced on the overall spectrum, authors would improve their strategy efficiency.

The FSA design leads to a graphical representation which is achieved in three steps: (i) identification of the system's functionalities and definition of each control variable, (ii) finding the constraints (restriction of the system functionalities) that are influenced by system's variables and (iii) designing the structural graph.

Water management faults is a recurrent faulty condition for PEMFC systems. Indeed, they can lead to severe performance losses and in some cases to irreversible degradation. Therefore, the FSA approach that leads to a structural graph considers two faulty conditions related to the water management: flooding and membrane drying out. The introduction of these two faults in the structural graph will underline the fuel cell functionalities impacted by their occurrence, and also highlight the available control variables which can be used for their mitigation in the case of an FTC strategy. FSA leads thus to the following contributions:


This paper is organized as follows: The first section is dedicated to the PEMFC water management issues. The second section proposes a structural analysis which allows highlighting all couplings between the system variables. Then, in the third section the structural analysis approach is defined and applied on the PEMFC system. The two last sections discuss and conclude about the structural analysis approach applied to a PEMFC system.

#### **2. Water Management Faults**

PEMFC systems may be subjected to different faulty operating modes. In [14], authors define fault as a decrease in system performance caused by improper major or minor fuel cell operation. These kinds of operations could lead to a permanent loss of the fuel cell performance due to the occurrence of faults or fuel cell ageing. This paper only focuses on faults which lead to fuel cells performance losses. These faults can be classified according to some criteria like: effects, response time, recovery property (the loss of performance can be totally recovered or not) or location.

Fuel cell faults can be detected and isolated, with the use of appropriate fault diagnosis tools. The literature relates several diagnosis techniques and many of them are used for PEMFCs water management issues. For instance, Lu et al. [15] proposed a fault diagnosis based on a fast electrochemical impedance spectroscopy (EIS) measurement. The developed tool allows an on-line flooding and drying out diagnosis but the authors underline that it cannot rules on the fuel cell state of health (SoH) in the case of multiple fault occurrence. Regarding their experiment, the multiple fault occurrence happens when the recovering time between each faulty condition is not enough. Therefore, in order to have a complete fault deletion and to improve the diagnosis tool performance, it is therefore relevant to know what the involved variables during a fault mitigation process are. Another example of a diagnosis tool consists of a model-based observer for fuel cell internal states estimation [16]. Authors aim to resolve the unmeasured internal variables issue thanks to a virtual sensor based on observer. This work highlights the need to understand the fuel cells operation mechanism through the internal state estimation and by the coupled variables involved in the change of the fuel cell SoH. In the follow-up, Alves-Lima et al. [17]

have worked on a quantitative video-rate hydration imaging of Nafion. Their results have shown that membrane water content is correlated to several factors: the membrane thickness, the fuel cell temperature (the room temperature in their study) and the water desorption process are major factors that influence the membrane hydration. Once again, the knowledge of the fuel cell internal variables are major issues to maintain the membrane under nominal hydration. In [18] authors use an EIS measurement to characterize the impact of the membrane water management on the performance of PEMFC commercial stacks. Their goal is to understand what the effects of the inlet gas water content on the fuel cell operations via the analysis of the Nyquist plot for fuel cell stack and single cells are.

These papers underline the need of knowledge about the internal fuel cell variables and how they have mutual influences. Indeed, as shown in [17], a modification of only one internal variable value could lead to the system destabilization. For this reason, the FSA analysis is focused on all the variables involved in the fuel cell to highlight extensively and systematically the multilateral effects of water management faults inside the fuel cell.

Cell flooding and membrane drying out are two possible consequences of improper water management. Flooding is caused by an accumulation of liquid water either in the diffusion layer of the electrodes, the bipolar plate channels or the feeding lines that limits the access of the reactants to the catalyst sites, and then decreasing the electrochemical reaction rates. Membrane drying out is the result of an insufficient hydration of the PEMFC membrane, thus increasing its proton resistivity.

Li et al. [19] identified flooding as the most recurrent PEMFCs' fault and point out that the cathode—being the place of water production—is particularly affected. The antagonist phenomenon is the drying out. It can occur when the gases relative humidity is too low [20,21], the input gas flow rate is too high, the operating temperature is too high or when these improper operating parameters are combined.

To better understand these faults, a functional analysis of PEMFC water management has been investigated through the Fault Tree Analysis in [11], among the possible approaches presented previously. The FTA is reported in Figure 1.

The FTA approach provides information about some coupling phenomena between variables inside the system. For instance, the temperature (T) influences the gas relative humidity. However, it is also linked to other variables which can be used to settle a fault tolerant control law. It is therefore important to understand how the different control variables are coupled to reduce the risk of unexpected phenomena. The proposed FSA approach aims to bring out this information in an explicit and systematic way via the design of diagrams called structural graphs.

**Figure 1.** *Cont*.

**Figure 1.** FTA applied on membrane drying out (**a**) and on a flooding (**b**) reproduced with the permission from Yousfi Steiner et al., 2021 [9].

#### **3. Fault Structural Analysis: Definition and Objectives**

FSA is a low-level representation of a system behavior that allows highlighting the operating conditions that potentially lead to a fault occurrence and the variables that could be used to mitigate this fault. The representation is based on connections between the system variables through a bipartite graph. All system features are described by a set of constraints, viewed as restrictions on the system functionalities, and the violation of one of them that indicates a fault occurrence.

#### *3.1. Dynamical Systems*

Structural analysis considers only the structural information, highlighting the variable involved in the studied phenomenon. A representation is given in Figure 2.

**Figure 2.** Known and unknown variables representation.

Inputs Ur are defined as: Ur with r {1,2, . . . ,n} Outputs Yi are defined as: Yi with i {1,2, . . . ,k} Unknown states, which are not directly measured but could be estimated from the known ones, they are defined by: xj with j {1,2, . . . ,l}.

#### *3.2. Bipartite Graph*

Beauguitte [22] defined the bipartite graph as a structure that displays relationships between two separate sets of vertices, describing system characteristics, variables, and constraints. For this reason, in the next section the bipartite graph will be called structural graph. The author underlines that vertices can be separated according to their contribution to an event such as the occurrence of a particular fault.

The links between vertices are usually not oriented and represent the system's structure which is noted as: G = (C U Z, Γ), where G is the structure of the bipartite graph and U the union operator. Z is the set of characteristic variables: Z {z1, z2, ... , zn}. C is the set of constraints: C = {c1, c2, ... , cm}. Arcs that connect each vertex (constraints/variables) are noted:

$$\Gamma = |(\mathbf{c\_i}, \mathbf{z\_j}) \vdash | \quad \mathbf{z\_j} \text{ exist in } \mathbf{c\_i}, \mathbf{z\_j} \in Z, \mathbf{i} \in [1, \mathbf{m}], \mathbf{j} \in [1, \mathbf{n}]). \tag{1}$$

An illustration of a bipartite graph is given in Figure 3.

**Figure 3.** Illustration of a bipartite graph.

It is possible to create an incidence matrix M that represents the bipartite graph and consists of a Boolean matrix where the rows are the constraints and the columns the characteristic variables:

$$\mathbf{M} = \{ \mathbf{m}\_{i,j} \: \mid \: \mathbf{m}\_{i,j} = 1 \text{ if } (\mathbf{c}\_i, \mathbf{v}\_j) \in \Gamma, 0 \text{ else} \} \tag{2}$$

#### *3.3. Differentiation*

The studied phenomena (faults) are time dependent. Therefore, a dynamic model is mandatory in order to take into account time dependent variables. Three options can be considered:

Option 1: Considering x (variable of the studied system) and . x as the same variable and treating the dynamical equations in the same way as the static ones.

Option 2: Considering x and x = dx/dt as structurally distinct and using the model with the explicit differential equation: dx/dt ([23]).

Option 3: Considering x and x = dx/dt as structurally distinct and proceeding to the structural differentiation of the initial model ([24]).

The choice to use the explicit differential equation . x is preferred in our study to take into consideration the system dynamic behavior with the present and past values (option 2).

#### *3.4. Representation of Faults*

As said before, faults are abnormal phenomena which decrease the system performance and can lead to its degradation. It is thus mandatory to detect their occurrence to proceed to a fault mitigation strategy. In case a diagnosis tool is used, several levels of a priori knowledge must be considered.


For this work, a low level of knowledge is considered. Indeed, accurate fault model is not mandatory. Only the functionalities that are influenced by their occurrence are needed whereas the cause-and-effect relationship is not.

The FSA approach and objectives are now explained. The structural graph design process is therefore described by the following steps: (i) find a model of the system which describes its functionalities, (ii) identify the system constraints in order to get all the variables which have an influence on the system operation, (iii) create an incidence matrix for the design of the structural graph. The next section consists in choosing a PEMFC model to guide the structural graph design.

#### **4. PEMFC Functionalities and Control Variables**

#### *4.1. PEMFC System*

A PEM fuel cell is composed of an anode, supplied with hydrogen, and a cathode, supplied with oxygen (pure or from the air). At the anode side, one molecule of hydrogen is oxidized thanks to a catalyst made of platinum (Pt). It allows the formation of two protons and two electrons:

$$\text{H}\_2 \to 2\text{H}^+ + 2\text{e}^- \tag{3}$$

At the cathode side, oxygen is reduced by the protons thanks to a catalyst made of platinum (Pt). Its reduction allows the formation of water through the electrochemical reaction:

$$\text{H}\_{2}\text{O}\_{2} + 2\text{H}^{+} + 2\text{e}^{-} \rightarrow \text{H}\_{2}\text{O} \tag{4}$$

The electrons move from the anode to the cathode through an external electric circuit. These two electrodes are separated by a proton exchange membrane that is impermeable to reactant gases. The supply of reactants to the PEMFC and the exhausts are carried in and out via channels in bipolar or mono-polar plates as represented in [25].

For the next sections, the studied system is a PEMFC stack. All cells are supposed to have the same mean behavior, and the reactant supply system is represented on the scheme on Figure 4:

**Figure 4.** Example of humidified system synoptic of a PEMFC.

The Figure 4 represents the main elements of the PEMFC system that will be involved in fault generations and then monitored for their mitigation. Indeed, the system is supplied with hydrogen and oxygen and the gas flows are controlled with two controllers which manage to supply the PEMFC with the good ratio of reactants. The reactants are humidified with the humidification system. The PEMFC input gas relative humidity is adjusted with a heated line. The PEMFC pressure is controlled with two back pressure controllers. To simplify the representation, in the following, the considered oxidant is pure oxygen, but the approach would be the same with air.

#### *4.2. PEMFC Modeling with an Energetic Macroscopic Representation (EMR)*

To perform an FSA for PEMFCs water management faults, a zero-D fuel cell model is considered. The energetic macroscopic representation (EMR) proposed in [26] is used to identify the system functionalities. This macroscopic model is chosen because it highlights the system's control variables. Indeed, the structural analysis presented in this paper aims to bring the relevant control variables which could mitigate a fault occurrence. Other multidimensional models would lead the same analysis, but they involve more variables which are not measurable or controllable. Therefore, this kind of model are useless for the aim of fault mitigation.

The EMR representation is a quasi-static model involving differential equations that make possible the study of time-dependent systems. It also specifies the equations that describe the system's functionalities that could be submitted to a faulty condition. The principle of the EMR consists in considering the different power conversions and the causeconsequence relationship. Then, each element is linked to the others through a couple of action/reaction variables which product is a power. The instantaneous conversion and the accumulation of power are distinguished.

The PEMFC is modelled as several sub systems linked together in order to create a complete fuel cell model. This simple model and few additional equations are relevant to perform the ASD because it allows taking into account the variables of interest.

This representation is composed of three parts.

First, the fluidic part is dedicated to the gas channels located between the gas tank and the reaction sites. An electric analogy is used with pressure is assumed to be a potential and volume flow rates are assumed to be currents. Using electrical analogy:


The Figure 5 shows the analogy between the pipe and the RC electrical circuit [27].

**Figure 5.** Analogy between pipe of the hydraulic circuit and an electric RC circuit.

The oxygen circulation qO2 is enabled by fuel cell input/output pressure difference. The consumed (resp. produced) gases' flows qcO2 are calculated from the Faraday law and the number of electrons exchanged nO2 counted negatively (resp. positively). Partial pressures at each electrode are imposed by the gas accumulation represented by qch. Ifc is the electric current of the fuel cell, N is the number of cells.

$$\text{PO2} = \text{PscO2} + \text{RdO21} \text{ qO2} \tag{5}$$

$$\text{qO2out} = (\text{PscO2} - \text{PsO2}) / \text{RdO22} \tag{6}$$

$$(\text{dPsC2})/\text{dt} = 1/\text{Ch} \left( \text{qC2} - \text{qcO2} - \text{qO2out} \right) \tag{7}$$

$$\text{qcO2} = \pm \text{N} \ast \text{Ifc/(nO2} \ast \text{F)} \ast (\text{R} \ast \text{Tfc})/\text{PO2} \tag{8}$$

$$\text{qO2out} = \text{qO2} + \text{qcO2} \tag{9}$$

In Equation (7) the derivative term is assumed to be structurally different from PscO2 (Cf. Section c).

The second part of the EMR is the electrochemical part. The Nernst potential En is based on the computation of a thermodynamic potential E0 and the influence of the partial pressures and the temperature:

$$\mathbf{E}\mathbf{n} = \mathbf{E}\mathbf{O} + \Delta\mathbf{E} \tag{10}$$

with,

$$\text{EO} = \alpha + \beta \, \* \, \text{Tfc} + \gamma \, \* \, \text{Tfc}^\* \text{2} + \delta \, \* \, \text{Tfc}^\* \text{3} + \text{v} \, \* \, \text{Tfc} \* \ln(\text{Tfc}) \tag{11}$$

and,

$$\Delta \mathbf{E} = \text{Accd} \ast \ln\left(\text{PscH2}/\text{P0}\right) + \text{Bcd} \ast \ln(\text{PscO2}/\text{P0}) \tag{12}$$

where α, β, γ, ν, Acd(Tfc) and Bcd(Tfc) are model adjustment variables. Tfc is the fuel cell temperature and the chosen model does consider a cooling system. E0 is the thermodynamic potential of the PEMFC [26]. The voltage drop ΔV calculation (activation, concentration and ohmic losses) is then carried out:

$$
\Delta \mathbf{V} = \mathbb{E} \Delta \mathbf{V} \, \mathbb{Z}\_{-} \mathbf{a} \mathbf{c} \mathbf{t} + \mathbb{E} \Delta \mathbf{V} \, \mathbb{Z}\_{-} \mathbf{c} \mathbf{n} \mathbf{c} + \mathbb{E} \Delta \mathbf{V} \, \mathbb{Z}\_{-} \mathbf{o} \mathbf{m} \, \tag{13}
$$

The PEMFC voltage is thus:

$$\mathbf{V}\_{-}\mathbf{M} = \mathbf{E}\_{-}\mathbf{n} - \Delta\mathbf{V} \tag{14}$$

Each loss is expressed as follows, using In the cross over current, I0 the exchange current, and Il the limit current:

$$
\Delta \text{Vact} = \text{A} \, \ast \, \text{Tfc} \, \ast \, \ln((\text{Ifc} + \text{In})/\text{I0}) \tag{15}
$$

$$
\Delta \text{Vcornc} = \text{B} \ast \text{Tfc} \ast \ln(1 - \text{Ifc}/\text{II}) \tag{16}
$$

$$
\Delta \text{Vohm} = \text{Rm} \ast \text{Ifc} \tag{17}
$$

The third part is the electric impedance of the cell, represented by the block "Charge double layer" but as the fast-electric dynamic is not considered in the analysis, this block is not considered. Considering N the number of stack cells, the stack voltage is expressed as follows:

$$\text{Vfc} = \text{N} \star \text{ (VM} + \text{Vc)} \tag{18}$$

Thanks to this fuel cell model, it is possible to define the control variables for the system functionalities. All variables that influence the PEMFC functionalities are supposed to be known. The next step consists of defining the constraints that can be influenced by a fault occurrence.

#### **5. Constraints of PEMFC**

#### *5.1. Structural Analysis of a PEMFC*

The structural analysis design is focused on the cathode fluidic part because it is the location of the water production. The cathode area therefore represents the major issue regarding fuel cells water management. It is located between the gas tank and the reaction sites. This part is modeled with the electric analogy represented on Figure 5. The oxygen circuit is composed of a fluidic resistance RdO21, an exhaust resistance RdO22, a hydraulic capacity ChO2 and the consumed oxygen flow generator. The water circuit is composed of the same type of elements, with the produced oxygen flow generator. The input gas flow qO2 is imposed by a flow controller, the values of the input water flow qH2Oin is imposed by the humidification system depending on the controlled humidity rate and the temperature of the fuel cell Tfc. The pressure at the exhaust is the atmospheric pressure.

The gas flow inside the PEMFC supply channels is considered as the first event for the normal fuel cell operation. This first event, called constraint for the structural analysis, is the first constraint noted C1. However, this gas circulation can be influenced by the system variables given in the following equations:

$$\rm P\_{O2in} - P\_{O2out} = R\_{dO21} \, \* \, q\_{O2in} + R d\_{O22} \, \* \, q\_{O2in} - q\_{O2in} \tag{19}$$

where PO2in and PO2out are the pressure of the input and output gas respectively. Regarding the water circuit the equation becomes:

$$\rm P\_{\rm H2Con} - P\_{\rm H2Conit} = R\_{\rm dH2O1} \, \ast \, \rm q\_{\rm H2Con} + R\_{\rm dH2O2} \, \ast \, \rm q\_{\rm H2Con} + q\_{\rm H2Oân}$$

The cathode pressure drop is expressed as follows:

$$
\Delta \mathcal{P} = \mathcal{P}\_{\text{tot\\_in}} - \mathcal{P}\_{\text{tot\\_out}}.
$$

Using the humidified input gas flow, which is actually the total input flow (qO2hum\_in), the pressure drop expression becomes:

$$\Delta \mathbf{P} = \mathbf{R}\_{\mathrm{dO2hum}\_1} \mathbf{q}\_{\mathrm{O2hum}\_{\mathrm{in}}} + \mathbf{R}\_{\mathrm{dO2hum}\_2} \left( \mathbf{q}\_{\mathrm{O2hum}\_{\mathrm{in}}} - \mathbf{q}\_{\mathrm{O2}\_{\mathrm{ca}}} + \mathbf{q}\_{\mathrm{H2O}\_{\mathrm{ca}}} \right).$$

Thus:

$$\text{C1}: \,\mathrm{q}\_{\mathrm{O2hum}\_{\mathrm{in}}} = \frac{1}{\mathrm{R}\_{\mathrm{dO2hum}\_{1}} + \mathrm{R}\_{\mathrm{dO2hum}\_{2}}} \left(\mathrm{\Delta P} + \mathrm{R}\_{\mathrm{dO2hum}\_{2}}\mathrm{q}\_{\mathrm{O2\_{ca}}} - \mathrm{R}\_{\mathrm{dO2hum}\_{2}}\mathrm{q}\_{\mathrm{H2O\_{ca}}}\right) \tag{20}$$

The humidified gas inside the channels crosses the gas diffusion layer (GdL) for the purpose of air diffusion to the catalytic site. However, the GdL gas concentration has to be kept at a high value for a good and safe catalytic feeding. The gas amount accumulated in the GdL is thus the third constraint C2. The variables which have an influence on this gas amount inside the GdL appear in the expression of the hydraulic capacity of the GdL Ch:

$$\mathbf{C\_h} = \frac{1}{\mathbf{P\_{O2\_{ca}}} + \mathbf{P\_{H2O\_{ca}}}} \left( \mathbf{q\_{O2}} + \mathbf{q\_{H2O\_{ca}}} - \mathbf{q\_{O2}} + \mathbf{q\_{H2O\_{ca}}} - \mathbf{q\_{O2\_{out}}} - \mathbf{q\_{H2O\_{out}}} \right) \tag{21}$$

where:

$$\mathbf{q}\_{\text{O2hum}\_{\text{out}}} = \frac{\mathbf{P}\_{\text{O2}\_{\text{ca}}} + \mathbf{P}\_{\text{H2Oca}} - \mathbf{P}\_{\text{O2}\_{\text{out}}} - \mathbf{P}\_{\text{H2O}\_{\text{out}}}}{\mathbf{R}\_{\text{dO2hum}\_{2}}} \tag{22}$$

hence:

$$\text{C2: } \mathsf{C}\_{\mathrm{h}} = \frac{1}{\mathsf{P}\_{\mathrm{O}\mathbb{2}\_{\mathrm{a}}} + \mathsf{P}\_{\mathrm{H}\mathrm{2O}\_{\mathrm{a}}}} \left( \mathsf{q}\_{\mathrm{IO}\_{\mathrm{m}}} + \mathsf{q}\_{\mathrm{H}\mathrm{2O}\_{\mathrm{m}}} - \mathsf{q}\_{\mathrm{IO}\_{\mathrm{2}}} + \mathsf{q}\_{\mathrm{H}\mathrm{2O}\_{\mathrm{a}}} - \frac{\mathsf{P}\_{\mathrm{O}\mathbb{2}\_{\mathrm{a}}} + \mathsf{P}\_{\mathrm{H}\mathrm{2O}\_{\mathrm{a}}} - \mathsf{P}\_{\mathrm{O}\mathrm{2}\_{\mathrm{out}}} - \mathsf{P}\_{\mathrm{H}\mathrm{2O}\_{\mathrm{out}}}}{\mathsf{R}\_{\mathrm{dO}\mathrm{2haun}\_{\mathrm{2}}}} - \mathsf{q}\_{\mathrm{H}\mathrm{2O}\_{\mathrm{out}}} \right) \tag{23}$$

It should be noted that the constraint C2 has also the fuel cell temperature as input variable. Indeed, the flowrates are expressed as volume so they depend on Tfc.

At the catalytic sites, and for high current densities, the kinetic of reactions increases and the gas consumption intensifies. During this operating condition, the quantity of O2 species decreases and the steam production increases. Therefore, the fuel cell voltage decreases. The steam partial pressure at the catalytic site is thus the fifth constraint C3 because it has to be under the saturation pressure value to avoid water condensation. The variables which influence this constraint are expressed by the following inequality relationship which represent the steam partial pressure at the catalyst (PH2Oca) and the saturation pressure (Psat):

$$\mathcal{C}3: \text{ } \mathcal{P}\_{\text{H2Oca}} \le \mathcal{P}\_{\text{sat}} \tag{24}$$

Then as depicted in [28], there is a thermodynamic equilibrium between the GdL [29] water content (λGdL) and the membrane water content λ<sup>m</sup> ([30]). The constraint C4 represents this equilibrium. Variables which have an influence on C4 are expressed below ([30]):

$$\mathbf{C4}: \ \lambda\_{\mathrm{GdL}} = \mathbf{a}\_1 + \mathbf{a}\_2 \left(\frac{\mathbf{P}\_{\mathrm{H2O}\_{\mathrm{ca}}}}{\mathbf{P}\_{\mathrm{sat}}}\right) - \mathbf{a}\_3 \left(\frac{\mathbf{P}\_{\mathrm{H2O}\_{\mathrm{ca}}}}{\mathbf{P}\_{\mathrm{sat}}}\right)^2 + \mathbf{a}\_4 \left(\frac{\mathbf{P}\_{\mathrm{H2O}\_{\mathrm{ca}}}}{\mathbf{P}\_{\mathrm{sat}}}\right)^3 \tag{25}$$

Regarding the membrane, it has to be hydrated for an appropriate fuel cell operation. An electro osmotic flow (qosm) allows the membrane water supply via a protonic load flow [10]. This flow must get a sufficient value for a good membrane hydration. This is the fifth constraint C5 which involves the electro osmotic flow. Then, a diffusive flow (qdiff) through the membrane also exists and modifies the water content of the membrane. Like the electro osmotic flow, it has to get a relevant value for a good membrane hydration. This value constitutes the sixth constraint C6. The constraint C5 and C6 are expressed as below:

$$\mathbf{q}\_{\rm H2Oca} = \mathbf{q}\_{\rm osm} + \mathbf{q}\_{\rm diff} \tag{26}$$

with:

$$\text{C5}: \ q\_{\text{osm}} = \lambda\_{\text{m}} \tau\_0 \frac{\mathbf{I}\_{\text{fc}}}{\text{F}} \tag{27}$$

and:

$$\text{C6}: \, \mathbf{q}\_{\text{diff}} = -\mathbf{D}\_{\text{m}} \frac{\rho\_{\text{dry}}}{\mathbf{M}\_{\text{m}}} \frac{\mathbf{d}\lambda\_{\text{m}}}{\mathbf{dx}} \tag{28}$$

An optimal membrane hydration is mandatory for the nominal operation of the PEMFC. Its water content is related to the water concentration inside the membrane. A constraint C7 is thus considered for the membrane water content. This constraint is influenced by the water concentration variation and is expressed as below ([31]):

$$\text{CT}: \,\lambda\_{\text{m}} = \frac{\text{M}\_{\text{m}}}{\rho\_{\text{dry}}} \mathbf{c}\_{\text{H2O}} \tag{29}$$

where cH2O is the membrane water concentration which depends on the fuel cell temperature.

The membrane hydration also depends on the thermodynamic equilibrium; a relationship between the GdL water content (λGdL) and the membrane water content λ<sup>m</sup> exists. The constraints C5 and C6 are therefore linked to the fuel cell temperature.

The last constraint C8 is also about the membrane water content. Indeed, the lower the membrane water content is, the higher the ohmic resistance is. This resistance can be expressed as below ([31]):

$$\text{C8}: \ R\_{\text{m}} = \frac{\text{t}\_{\text{m}}}{\left(\text{b}\_{1} \exp\left(\text{b}\_{2}\left(\frac{1}{303} - \frac{1}{\text{T}\_{\text{fc}}}\right)\right)\right)}\tag{30}$$

where, tm represents the membrane thickness, b1 and b2 are coefficients that depend on fuel cell being tested and b1 depends on the membrane water content ([31]):

$$\mathbf{b}\_1 = \mathbf{b}\_{11}\lambda\_m - \mathbf{b}\_{12} \tag{31}$$

Then, the higher the membrane resistance is, the higher the ohmic losses are. The ohmic losses is written as:

$$\text{C8} \colon \Delta \text{V}\_{\text{ohm}} = \text{R}\_{\text{m}} \text{I}\_{\text{fc}} \tag{32}$$

The incidence matrix can now be set in order to create the PEMFC structural graph.

#### *5.2. Incidence Matrix of the Structural Analysis*

Based on the extraction of the PEMFC constraints, it is possible to create an incidence matrix A (ai,j) that allows to link vertices (variables/constraints) and arcs. It contains n rows and m columns:

	- The incidence matrix is represented in the Table 1.

**Table 1.** Incidence matrix of the structural analysis.




The control variables are separated from the others with a double vertical line. The structural graph is designed on the basis of the incidence matrix of the Table 1. It is represented on the Figure 6 below.

The next step consists of adding the PEMFC faults on the structural graph in order to represent their interaction inside the PEMFC.

#### *5.3. Fuel Cell Flooding Structural Analysis*

A PEMFC flooding can occur in two areas. Both inside the GdL with a water droplet accumulation which reduces the catalytic site reactant feeding, and inside the channels by propagation of water droplet accumulation. This water accumulation can also appear directly inside the supply channels and reduce the fuel cell reactant feeding. To integrate the flooding in the structural analysis, the fault is assumed to be a variable which has

influence on the PEMFC constraints defined above. For this purpose, a new variable (Fflood) which represents a flooding occurrence is added. Fflood can therefore be expressed as below:

$$F\_{\text{flood}} = \frac{V\_{\text{electrode\\_available}}}{V\_{\text{gco\\_electrode}}} \tag{33}$$

with, Velectrode\_available is the global volume available in the compartment (channels + GdL) and Vgeo\_electrode the geometrical volume of the electrode. This variable is set to 1 in case of optimal hydration and to 0 when completely clogged.

The flooding variable has an influence on the supply channels. The new variable is thus added to the constraint C1 which is linked to PEMFC supply channels:

$$\text{C1}: \,\mathrm{q}\_{\mathrm{O2hum}\_{\mathrm{lin}}} = \frac{\mathrm{F\_{\text{flood}}}}{\mathrm{R\_{\mathrm{dO2hum}\_{\mathrm{I}}} + \mathrm{R\_{\mathrm{dO2hum}\_{\mathrm{2}}}}} \left(\mathrm{\Lambda P} + \mathrm{R\_{\mathrm{dO2hum}\_{\mathrm{2}}}} \mathrm{q}\_{\mathrm{O2\_{\mathrm{ca}}}} - \mathrm{R\_{\mathrm{dO2hum}\_{\mathrm{2}}}} \mathrm{q}\_{\mathrm{H2O\_{ca}}}\right) \tag{34}$$

The flooding variable has also an influence on the GdL and thus on the constraints C2 and C4. These two constraints become:

$$\text{C2: C}\_{\text{H}} = \frac{\text{F}\_{\text{Hood}}}{\text{P}\_{\text{O2}\_{\text{ca}}} + \text{P}\_{\text{H2O}\_{\text{ca}}}} \left( \mathbf{q}\_{\text{O2}\_{\text{an}}} + \mathbf{q}\_{\text{H2O}\_{\text{n}}} - \mathbf{q}\_{\text{O2}\_{\text{ca}}} + \mathbf{q}\_{\text{H2O}\_{\text{ca}}} - \frac{\mathbf{P}\_{\text{O2}\_{\text{ca}}} + \mathbf{P}\_{\text{H2O}\_{\text{cat}}} - \mathbf{P}\_{\text{O2}\_{\text{cat}}} - \mathbf{P}\_{\text{H2O}\_{\text{cat}}}}{\text{R}\_{\text{dO2}\text{num}\_{2}}} + \mathbf{q}\_{\text{H2O}\_{\text{cat}}} \right) \tag{35}$$

$$\mathbf{C4}: \ \lambda\_{\mathrm{GdL}} = \mathrm{F\_{flood}}\left(\mathbf{a}\_1 + \mathbf{a}\_2 \left(\frac{\mathrm{P\_{H2O\_{ca}}}}{\mathrm{P\_{sat}}}\right) - \mathbf{a}\_3 \left(\frac{\mathrm{P\_{H2O\_{ca}}}}{\mathrm{P\_{sat}}}\right)^2 + \mathbf{a}\_4 \left(\frac{\mathrm{P\_{H2O\_{ca}}}}{\mathrm{P\_{sat}}}\right)^3\right) \tag{36}$$

The incidence matrix is updated with the flooding variable as represented in the Table 2.

**Table 2.** Incidence matrix updated with the flooding variable.




The structural graph with the flooding variable is represented on the Figure 7:

The next step consists in introducing the membrane drying out variable inside the structural analysis. The goal is to identify the constraints which are influenced with its occurrence.

**Figure 7.** Structural graph with the flooding variable.

#### *5.4. Membrane Drying out Structural Analysis*

The membrane drying out fault results from a decrease of the membrane water content. Therefore, it has an influence on constraints C7 and C8 that involve the membrane water content variable (λ). The fault can thus be represented by a variable which is expressed as below:

$$F\_{\rm dry} = \frac{V\_{\rm abs\\_membrane}}{V\_{\rm abs\\_tot\\_membrane}}\tag{37}$$

where, Vabs\_membrane is the volume absorbed by the membrane and Vabs\_tot\_membrane is the total volume absorbable by the membrane. Fdry has a value between 0 and 1.

This variable can now be introduced in the constraint C9:

$$\mathbf{C}\mathbf{\theta}:\ \lambda\_{\mathbf{m}} = \frac{\mathbf{M}\_{\mathbf{m}} \cdot \mathbf{F}\_{\mathrm{dry}}}{\varrho\_{\mathrm{dry}}} \mathbf{c}\_{\mathrm{H2O}}\tag{38}$$

The incidence matrix is updated with the membrane drying out variable as represented on the Table 3.


**Table 3.** Incidence matrix updated with the membrane drying out variable.

The structural graph with the membrane drying out variable is represented on the Figure 8:

**Figure 8.** Structural graph with the membrane drying out variable.

#### **6. FSA in Experimental Context**

The experimental FSA illustration is made through the mitigation of fuel cell flooding fault. Indeed, an active fault tolerant control (AFTC) strategy which is detailed in another work [32], is applied on a single cell fuel during a flooding.

AFTC strategies, compared to other kinds of FTC strategies, differs mainly by their structure. Indeed, in a previous literature review [1] it has been highlighted that there are two kinds of FTC strategies: Active or Passive ones. Passive strategies (PFTC) consist of a robust controller design for fault mitigation. In this case, the controller is designed by considering the fault as a disturbance. This structure allows to not use diagnosis tools. But as more is important the number of faults that should be mitigated as more is the PFTC design complexity. To reduce this complexity and instead of use a unique robust controller, the AFTC structure decomposed the mitigation process to several steps. For instance, in [33] authors proposed a three-modules fault mitigation process on a powertrain city bus. Authors used as first module a diagnosis tool for fault detection whereas the second module is composed of decision-making part. Finally, the third module is composed by a set of controllers which implement on the powertrain city bus the mitigation strategy computed by the second module. Another work in [23] where authors proposed another application of the AFTC strategy. In this case the purpose is to address the water management issues. They manage to couple a fault detection and isolation (FDI) algorithm with a reconfiguration mechanism and an adjusting controller. The FDI process is a machine learning based whereas a self-tuning PID is implemented as the control part. Authors also highlight the advantages of the self-tuning PID which shows robustness against noise and model uncertainties. Wu et al. proposed in another fault mitigation process which also consists of an implementation of a three-module AFTC to diagnose a PEMFC fault occurrence, to decide on a relevant mitigation action, and to apply the strategy on the fuel cell through a control set. The main purpose of the method is the distribution of the complexity of the AFTC design into the three modules which significantly simplifies its implementation. In [13] authors proposed a model-based AFTC strategy for PEMFC temperature sensor

fault. They used a FDI process with the aim of real time fault diagnosis with a sliding mode controller for the fuel cell thermal management.

The advantage of the use of AFTC structure lies in its modular aspect. Indeed, it allows the decomposition of the mitigation process into several steps. Each step being dedicating to a different task, it reduces the complexity of the strategy.

#### *6.1. Experimental Set Up*

The fault tolerant control strategy is applied on a test bench of the FCT brand [34]. It manages to supply a fuel cell with hydrogen at the anode and oxygen or air at the cathode. The experimentations are carried out only with oxygen. The tested fuel cell is a 50 W unit which is composed of N117 ION POWER single-cell of 50 cm<sup>2</sup> [35]. The synopsis of test bench is given on the Figure 4.

The Figure 4 shows that the test bench is mounted with two flow rate controllers in order to regulate the input reactant flows. Then, two humification systems are placed on O2 and H2 feeding lines in order to control their relative humidity thanks to two gas temperature controllers. At the event of the fuel cell, two back pressures controllers are used at the anode and cathode lines to manage the fuel pressure. An electronic load is finally connected to the FC electric terminals which can absorb a power of 1.8 kW. The whole test bench is monitored with a LabVIEW virtual instrument.

#### *6.2. Flooding Generation Test*

The experiment consists of mitigating a flooding occurrence with an AFTC strategy. It will modify the fuel cell operating conditions by changing some control variables iteratively. The selection of the control variables is based on the previous structural graph and on the available sensors and actuators on the test bench. The selected control variables are the: input cathode gas flow rate (qO2ca); fuel cell temperature (Tfc); inlet gas flow relative humidity (qH2Oca).

The structural graph shows that the flooding affects the fuel gas channels and the GDL on the fuel gas channels and on the GDL. The control variables qO2ca and qH2Oca have an influence on the gas channels whereas Tfc influences the GDL.

In the experiment, the flooding is generated by introducing liquid water in the cell from the canalization which is between the humidifier and the fuel cell. Indeed, on Figure 4 the gas at the outlet of the humidification system is temperature-regulated in order to reach the relative humidity setpoint. With the temperature controller, it is possible to condense the steam. Then the condensed water goes into the fuel cell to cause a flooding inside the electrode.

The used algorithm for the AFTC strategy is based on an iterative modification of the operating condition with the modification of the selected control variables. The Figure 9 is an illustration of the AFTC mechanism.

#### *6.3. AFTC Application to Flooding*

The used AFTC strategy has been selected from a previous literature review [1] which is constituted of three modules: diagnosis; decision; control. The diagnosis module is used for the flooding identification based on some outputs fuel cell measurements. Then, if a fault occurred, a decision process is launched through the decision module. At this step, a mitigation strategy is computed for a fast and sustainable fault mitigation. The decision strategy is finally implemented on the fuel cell through the control module which is composed of a set of controllers.

**Figure 9.** Iterative AFTC strategy.

In this work, the diagnosis consists of identifying a flooding occurrence by monitoring the fuel cell pressure drop and the voltage. The decision process provides a decision regarding the fault occurrence to proceed fast mitigation with minimal change in the operating point. The guidelines which are the output of the decision process are transmitted to all controllers to apply the mitigation strategy on the fuel cell.

The AFTC iterative process leads to the setup of two testing cases for the flooding mitigation. Figure 10 represents these testing cases.

**Figure 10.** Two possible test cases for flooding mitigation.

In both cases, the three selected control variables are used by the AFTC strategy for fault mitigation. If a fault is identified by the diagnosis block, one of the variable values is modified. For instance, in the case 1, the strategy starts with a correction of the gas flow rate qO2ca (1). This is referred to as a Transient action because qO2ca is temporally modified until the flooding is mitigated. The second step of the decision process (2) is an action on Tfc and on qO2ca. Here, Tfc is permanently modified to a new value. In parallel with the action on Tfc, another Transient action on qO2ca is triggered. The third step (3) of the case 1 consists of a permanent modification of the input gas water content (qH2Oca). This action is also supported by a Transient action on qO2ca. The fourth step (4) consists of a permanent change of qO2ca. Here, there is no Transient action.

The test case 2 consists of a different order of the corrective actions. It is composed of the same actions as for the case 1, in a different sequence; except for the first action which is a Transient action on qO2ca.

#### *6.4. Experimental Validation of the AFTC Strategy Based on the Variables Extracted from the Structural Graph*

A flooding mitigation process based on the test case 1 has been tested on a single fuel cell. Results are depicted in Figure 11a–c.

**Figure 11.** *Cont*.

**Figure 11.** Pressure drop (**a**), oxygen flow rate (**b**) and temperature (**c**) superposition with the fuel cell voltage for the test case 1.

The Figure 11a shows the superposition of the fuel cell voltage with the cathode pressure drop (PD). The PD evolution is used for flooding diagnosis. For instance, F1 corresponds to an increase in PD, which means flooding is in progress.

On Figure 11b all the applied Transient actions on O2 flow rate qO2ca are marked with the AFTC. Each increase in qO2ca aims to mitigate the flooding. When the flooding is mitigated (not diagnosed anymore), the qO2ca value is reset to its value before the Transient action triggering.

On Figure 11c, the second mitigation action of the test case 1, that is the fuel cell temperature is triggered (AFTC 2) at the same time as a Transient action

Figure 12 represents a flooding mitigation process based on the test case 2.

The Figure 12 represents the experimentation applied to the test case 2. Here the first mitigation action is always based on a Transient action, on qO2ca. The second action is also based on the same control variable qO2ca. The increase in qO2ca (AFTC 3) with the Transient action manages to mitigate the flooding by causing the liquid water to drain from the fuel cell.

**Figure 12.** *Cont*.

**Figure 12.** Pressure drop (**a**), oxygen flow rate (**b**) and temperature (**c**) superposition with the fuel cell voltage for the test case 2.

The test results are summarized in the Table 4.

**Table 4.** Test results for flooding mitigation for the tests cases 1 and 2.


#### *6.5. Experimental Analysis*

The test of the case 1 and 2 allows the flooding mitigation in about 15 to 18 min and for the same number of Transient actions and permanent actions. The main difference concerns the voltage recovery which is three times fastest in the case of an action on temperature rather than on the cathode flow rate.

Indeed, regarding the structural graph it appears that a flooding acts on the gas feeding channels and on the GDL. Actions on qO2ca allow the water to drain from the fuel cell channels but do not allow the water to drain from the GDL. In the same way, Tfc manages to remove water from the GDL and has no effects on the fuel cell channels. For these reasons, the case 1 which is a mix of actions on the channels and GDL through the qO2ca and Tfc, is a more efficient way for flooding mitigation than case 2 which only acts on the channels.

#### **7. Discussion**

The structural analysis gives a graphical representation of all PEMFC variables that have an influence on the fuel cell operating conditions. It highlights the coupling of variables inside the PEMFC and leads to the design of a structural graph. The structural graph of the PEMFC on the Figure 8 shows that the fuel cell temperature appears as the most coupled variable. It has an influence on the steam partial pressure at the catalyst site (C3), on the thermodynamic equilibrium of the GdL and the membrane water content (C4), on the electro osmotic and diffusive flows (C5 and C6) and on the membrane water content (C7 and C8). Other variables such as the membrane water flow rate and the input steam flow rate appear to be also strongly coupled. This information, brought by the structural graph, is very relevant to understand the PEMFC complexity and why it is so challenging to maintain fuel cell systems under nominal operating conditions.

Then, the flood variable (Fflood), which is considered as a fuel cell variable, is integrated in the structural graph. This variable has an influence on the input gas circulation inside the supply channels (C1) and on the GdL gas volume (C2). In case of a design of a flooding mitigation strategy, the structural graph gives the influenced constraints in order to get the appropriate control variables. The structural graph shows that the input gas pressure (Ptot) and flow rate have an influence on the gas feed channel and can lead to the mitigation of the flooding. Regarding a flooding diagnosis tool inside the gas feed channels, the measurement of the inlet gas pressure at the inlet and at the catalytic sites are relevant variables for the design. Indeed, when there is water accumulation inside the gas feed channels, the gas pressure difference between the inlet and catalytic sites increases. However, the pressure at the catalytic site is not measurable and for this reason the pressure difference is determined with the inlet and outlet pressures.

Regarding the flooding in the gas diffusion layer, the analysis conducted gives no relevant variable to diagnose the GdL flooding. This is due to the scale of the model which does not provides the fuel cell behavior inside the GdL and does not discriminate a clogging in the channels or in the GdL. However, the humidified inlet gas and the fuel cell temperature are relevant control variables for the GdL flooding mitigation with a fault mitigation strategy. Indeed, by decreasing the steam injected inside the fuel cell the water accumulation inside the GdL is restrained. An increase of the fuel cell temperature modifies the relative humidity inside the fuel cell and makes the water droplets evaporate.

Concerning the membrane drying out variable on the Figure 8, the structural graph shows that it is related to the membrane water content. The fuel cell current and temperature appear as relevant control variables on the structural graph for a mitigation strategy. These variables have an influence on the membrane water content and can be used in a fault mitigation strategy to alleviate the drying out occurrence. Indeed, the decrease of the fuel cell temperature leads to rehydrate the membrane by an increase of the relative humidity inside the fuel cell. The change of the current value has also an influence on the electroosmotic flow and thus on the water dragged through the membrane which participates to its rehydration. For this faulty condition, the fuel cell voltage, water content and current are relevant variables to set a membrane drying out diagnosis tool.

Therefore, the PEMFC FSA leads to the highlighting of the internal variable coupling. The integration of the fault variables in the structural graph allows understanding the fault process and their location. The analysis underlines the PEMFC constraints that are directly influenced by their occurrence and shows the relevant variables that can be used in FTC strategies.

#### **8. Conclusions**

The structural analysis approach aims to synthesize the known information about the PEMFC structure and water management issues, leading to flooding and drying out faults. Therefore, the analysis allows describing the system functionalities and their variables. It allows the graphical representation of the fuel cell strong coupling and provides information about relevant variables for the design of a diagnosis tool for flooding and membrane drying out. This analysis also highlights the relevant variables for the design of FTC strategies. Indeed, in case of the development of PEMFC monitoring or mitigation tools, the knowledge of the relevant variables that have to be used and given by the FSA is very helpful. The relevance of the FSA and particularly the graph is shown by the experimental mitigation of a flooding occurrence on a single fuel cell. Indeed, the experiment shows that it is very important to consider which control variable has an influence on which fuel cell functionalities in order to introduce in the FTC strategy enough control variables for fault mitigation.

As discussed in the previous section, the information contained in the FSA depends on the accuracy of the chosen model. But if a higher level was used for the structural graph design, maybe other control variables could appear in the structural graph. For this work and for the used AFTC strategy, only a low level of knowledge was needed because for the FTC design, the fact that not all variables can be measured nor estimated has to be kept in mind to lead to an implementable solution on the used test bench. Hence, there is a trade-off to be found between the different levels of knowledge regarding the selected models.

The literature also reports some other fuel cell faults like CO poisoning, hydrogen and oxygen starvation and low cathode and anode stoichiometry. The next step of the FSA design is therefore to complete the structural graph by adding these faulty conditions.

**Author Contributions:** Conceptualization, E.D., N.Y.S., M.B. and M.-C.P.; methodology, E.D., N.Y.S., M.B. and M.-C.P.; validation, E.D., N.Y.S., M.B. and M.-C.P., B.G.-P.; formal analysis, E.D., N.Y.S. and M.-C.P.; investigation, E.D., N.Y.S., M.B. and M.-C.P.; writing—original draft preparation, E.D.; writing—review and editing, E.D., N.Y.S., M.B. and M.-C.P.; supervision, E.D., N.Y.S., M.B., B.G.-P. and M.-C.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work has been supported by the Région Réunion and the Région Bourgogne Franche-Comté and by the EIPHI Graduate School (contract ANR-17-EURE-0002).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### MDPI

St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

MDPI Books Editorial Office E-mail: books@mdpi.com www.mdpi.com/books

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-6909-3