Modeling Electro-Chemo-Mechanical Behaviors within the Dense BaZr0.8Y0.2O3−δ Protonic-Ceramic Membrane in a Long Tubular Electrochemical Cell

This paper reports an extended Nernst–Planck computational model that couples charged-defect transport and stress in tubular electrochemical cell with a ceramic proton-conducting membrane. The model is particularly concerned with coupled chemo-mechanical behaviors, including how electrochemical phenomena affect internal stresses and vice versa. The computational model predicts transient and steady-state defect concentrations, fluxes, stresses within a thin BaZr0.8Y0.2O3−δ (BZY20) membrane. Depending on the polarization (i.e., imposed current density), the model predicts performance as a fuel cell or an electrolyzer. A sensitivity analysis reveals the importance of thermodynamic and transport properties, which are often not readily available.


Introduction
Materials such as doped barium zirconates (e.g., BaZr 1−x Y x O 3−δ , BZY) are perovskites that have good proton conductivity at intermediate temperatures (T 500 • C) [1] and chemical stability in environments containing H 2 O and CO 2 [1][2][3][4]. Thus, these materials are suitable as proton-conducting membranes for applications that depend on hydrogen separations (e.g., fuel cells) [5,6]. Protonic ceramic membranes have also been employed to compress hydrogen up to 50 bar. Methane is reformed in the presence of steam on the Ni-BZY based ceramic-metal support, and the produced hydrogen is galvanically driven through the membrane and compressed on the other side of the membrane [7]. Although BZY materials are predominantly proton conductors in moist reducing atmospheres, they are mixed ionic-electronic conductors (MIEC) in oxidizing atmospheres with the presence of O-site polarons. As the temperature increases above 600 • C, the conduction through oxygen vacancies becomes not negligible. The present paper focuses on the chemo-mechanical behavior of BaZr 0.8 Y 0.2 O 3−δ , called BZY20.
The concentration variations associated with charged-defect transport within the membrane cause lattice-scale strain and volume deformation, often referred to as chemical expansion [8,9]. Because ceramic materials are brittle, defect-induced stresses can initiate distortions, cracks, and fracture. Moreover, as discussed in early works by Lee [10] and Larché and Cahn [11,12], local stress gradients also contribute to defect fluxes [13].
Atkinson et al. [14,15] analyzed the effects of defect-induced stress in solid-oxide fuel cells (SOFCs) and oxygen-ion-conducting membranes with alternative geometries and mechanical constraints. Euser et al. developed an extended Nernst-Planck-Poisson (NPP) model to evaluate stress in planar and radial oxygen-separation membranes [16][17][18]. Andersson et al. [19] and Han et al. [20] used X-ray diffraction to measure the chemical expansion of BZY for different dopant levels in moist and dry environments. Marrocchelli et al. [21] and Bishop et al. [22] reported detailed descriptions of chemical expansion in the perovskite materials. Dubois et al. [23] developed a coupled chemo-thermo-mechanical model using an extended Nernst-Planck formulation and evaluated the stresses associated with a protonic-ceramic fuel cell (PCFC) in a button-cell configuration. Ricote et al. [24] used high-temperature X-ray diffraction (HT-XRD) and thermogravimetric analysis (TGA) to study the effects of oxygen partial pressure on lattice parameters and oxygen nonstoichiometry in BaZr 0.9 Dy 0.1 O 3−δ (BZDy10).
The current study develops an electro-chemo-mechanical computational model based on an extended Nernst-Planck (NP) formulation. The one-dimensional radial model considers a tubular configuration with a thin dense BZY20 membrane supported by a relatively thick porous composite electrode ( Figure 1). The model predicts transient profiles of defect concentrations and stress states as functions of applied electrical current. As with all models, results depend on physical parameters, including thermodynamic and transport properties. A sensitivity analysis seeks to identify the most important parameters and properties. Proton-conducting tubular cell consisting of a dense BZY20 electrolyte supported by porous Ni-BZY20 operating as a fuel cell.

Defect Chemistry
The structures and behaviors of proton-conducting ceramics have been documented in numerous reviews [1,2,19,[25][26][27][28][29][30][31]. Oxides with ABO 3 perovskite structures, such as BaZrO 3 , are among the most promising proton conductors. Effective proton conductivity is established by B-site doping with a trivalent cation such as yttrium (Y 3+ ) to replace a fraction of the Zr 4+ , thus introducing oxygen vacancies to preserve local charge neutrality. The present study considers 20% Y doping, BaZr 0.8 Y 0.2 O 3−δ . The needed thermodynamic, transport, and mechanical properties are derived from prior publications [1,2,4,19]. Figure 1 illustrates aspects of the tubular geometry and the transport processes. The profiles and fluxes of three mobile charged defects (protons OH • O , oxygen vacancies V •• O and polarons O • O ) depend on the current density via an external circuit. Depending on the polarization (i.e., the direction of the electrical current), the cell may operate as a fuel cell or as an electrolyzer [32]. As illustrated in Figure 1, the present study considers gas-phase mixtures of 20% H 2 O and 80% O 2 on the exterior electrode (positrode) and 97% H 2 and 3% H 2 O on the interior electrode (negatrode).
Defect-incorporation reactions at the membrane boundaries may be stated in Kröger-Vink notation as 1 2 Corresponding equilibrium constants for the incorporation reactions are where p k are the gas phase partial pressures and the defect concentrations are represented in lattice units. The defect molar concentrations may be evaluated using the molar volume (for BZY20, V m = 4.57 × 10 −5 m 3 mol −1 ) as The present model includes polaron traps [4]. In other words, some fraction of the otherwise highly mobile small polarons can be trapped, thus immobilized in the proximity of the yttrium dopant. Stated as a reaction, The associated equilibrium constant follows as The present model also assumes that the gas phases are equilibrated on both sides of the cell, Because the gas phase is equilibrated, the defect-incorporation equilibrium constants are not all independent; they are constrained as The equilibrium constants for each reaction can be evaluated via thermodynamic properties as where T is temperature, R is the gas constant, ∆G • is the change in Gibbs free energy, ∆S • is the change in defect entropy , and ∆H • represents change in enthalpy. Table 1 lists the thermodynamic properties used in the present study [4]. Table 1. Enthalpy and entropy of defect incorporation reactions in BZY20. Reproduced with permission from Zhu et al. [4]. Copyright The Electrochemical Society, 2018.

Reaction
Assuming the defect reactions at the membrane surfaces are at equilibrium, defect concentrations at surfaces are constants that depend on the gas phase composition. Evaluating the boundary concentrations must consider site and electroneutrality balances. Charge neutrality requires that The ABO 3 pervoskite structure further constrains the oxygen site balance as The yttrium doping is fixed as Y Zr • L = 0.2, which also constrains the trap balance as The defect concentrations at the boundaries are calculated by solving Equations (5), (6), (9), (13)- (15) simultaneously. Algorithmic details about the calculation may be found in Zhu et al. [4].

Extended NP Membrane Model
The defect concentrations are governed by a conservation principle, represented mathematically as where t is time, J k are defect transport fluxes, andω k are the net molar production rates of the polaron traps within the membrane [4]. In the current work, standard Nernst-Planck fluxes are extended to incorporate the contribution from hydrostatic-stress gradients. The defect-transport fluxes include the effects of diffusion J D k , migration, J M k , and hydrostatic stress J S k , The flux may be represented as where D k are defect diffusion coefficients, z k are the defect's charge, Φ e is the electrostatic potential, β k is the coefficient of chemical expansion, and σ h is the hydrostatic stress. The defect charges are Hydrostatic stress within the membrane is defined as The polaron traps are produced and consumed according the reaction represented by Equation (8). The rate of progress of the trap reaction may be evaluated via mass-action kinetics asq where k f and k b are forward and backward rate constants, which are related through the equilibrium constant. The present model assumes k f = 10 6 . The species production rate due to the trap reaction can be expressed asω Y Zr

Chemical-Expansion Coefficient
The chemical-expansion coefficient β k is defined to be the change in lattice volume with respect to the changes in defect concentration. As reported by Andersson et al. [19], the present model uses Ricote et al. [24] showed that the O-site polarons produce negligible chemical expansion. The present model attributes all chemical expansion to the protons, neglecting lattice distortion due to polarons and oxygen vacancies.

Electrostatic Potential Φ e
The electrostatic potential may be calculated by solving the Gauss equation, where the electrostatic-potential field is E = ∇Φ e , ε 0 and ε r are vacuum and relative permittivities, respectively, and i is the external current density. However, the present model implements an approximately equivalent method by enforcing strict local electroneutrality. In this limit, the local charge flux must vanish as The electrostatic potential on one of the boundaries is set to a reference Φ e = 0, while the other boundary is set to impose the net current density i through the membrane. In other words, the charge flux through the membrane is exactly balanced by the current density i through the external circuit, This algebraic equation, together with Equation (16), is sufficient to determine the local electrostatic-potential gradient ∇Φ e .

Defect Diffusion Coefficients
The charged-defect diffusion coefficients are represented in an Arrhenius equation as where D • k and E k are pre-exponential factors and activation energies, respectively. As reported by Zhu et al. [4], Table 2 lists the values of D • k and E k for three mobile charged defect. Table 2. Pre-exponential factors and activation energies for diffusivities of mobile charged defects in BZY20. Reproduced with permission from Zhu et al. [4]. Copyright The Electrochemical Society, 2018.

Stress, Strain, and Displacement
The isotropic stress-strain relationship including chemically induced stress may be written as where σ is the stress tensor, is the strain tensor, I is the unit tensor, and ∆[X k ] represents the change in concentration with respect to a reference zero-strain state. The current study does not incorporate the residual strain related to the fabrication process [5,23]. The Lamé constant λ and shear modulus G are defined as In these expressions, E m is the Young's modulus, and ν m is Poisson's ratio of the membrane. The tube is assumed to be infinitely long, resulting in a plane-strain state. Consequently, the strain components in the axial z direction vanish, The strain-displacement relationships with respect to radial displacement u are expressed as Assuming a quasi-static stress field and the absence of body forces, the divergence of the stress must vanish, In one-dimensional radial coordinates, the stresses can be represented as ∂σ rr ∂r Considering the tube to be long and axisymmetric, variations in axial z and circumferential θ stresses can be neglected. The radial component of the stress equilibrium simplifies to ∂σ rr ∂r In practice, protonic-ceramic membranes should be as thin as possible to promote high proton fluxes (usually on the order of tens of microns). Such very thin membranes must be supported with relatively thick porous support structures that also serve as electrodes (cf., Figure 1). In reducing environments (e.g., anodes of protonic-ceramic fuel cells) the most widely used supports are composites of Ni-BZY. Such structures provide the mechanical strength and stability. The percolating Ni phase serves dual roles as the hydrogen-reduction catalyst and the electrical conductor. Table 3 lists the dimensions of tube used for the present study. The membrane is 10 µm thick and the porous support is 1 mm thick. The outer electrode is specified to be thin, but not explicitly included in the present model.
The current study, which specifically focuses on the electro-chemo-mechanical behavior within the membrane, neglects defect transport and chemical expansion within the porous support. Table 3 also includes mechanical properties of BZY20 and Ni-BZY20. These properties were determined experimentally using the ultrasonic measurements [23].
The mechanical aspects of the present model follow the approach reported by Euser et al. [17]. By substituting the stress-strain relationships (Equation (25)) and strain-displacement relationships (Equation (29)), the equilibrium equation (Equation (32)) can be solved analytically to produce the local displacements u, assuming a perfect bond between support and membrane properly satisfies continuity of displacement and traction at the membrane boundary. The constants of integration can be evaluated assuming inner and outer surfaces of the tube are traction-free. Using these displacements, the local stress components within the membrane and porous support can be evaluated.

Computational Implementation
The computational model is formulated using the method-of-lines and a one-dimensional radial finite-volume mesh network. The model incorporates two-way coupling between the defect transport and stress in the membrane. Transport induced stress is evaluated using the analytical solution. Euser et al. [17] reported the radial displacement and stress for the tubular membrane, respectively, as where C 1 and C 2 are constants of integration, which can be calculated using the mechanical boundary conditions, as discussed earlier. In addition, where r i is the inner radius of membrane. The transport problem is solved numerically. Equation (16) is effectively a parabolic partial differential equation that is solved in Matlab using the ode15i function. The second-order flux derivatives are approximated using conservative central differences. The migration J M k and stress J M k fluxes are discretized using an up-winding scheme, where the coefficient of [X k ] behaves as an artificial velocity [33,34]. The simulations run efficiently and quickly on a typical personal computer.

Proton Flux Profiles
Proton flux is an important measure of membrane performance. Figure 4 shows the predicted proton flux as a function of current density for selected oxygen concentrations at the tube's exterior electrode. Interestingly, under fuel-cell polarization (i > 0), the flux is only weakly affected by the O 2 concentration. The effect is much more significant under electrolysis polarization. As may be anticipated, high steam concentrations (i.e., low O 2 concentration) lead to higher proton fluxes under electrolysis conditions.  At steady state, the defect-flux profiles must be essentially flat (i.e., very little spatial dependence). However, because of the radial coordinates, there will be very slight radial variations in the steady-state fluxes. The net proton flux comprises diffusion, migration, and stress contributions (Equation (17)). These flux contributions do vary spatially, even at steady state. Figure 5 plots separately the individual steady-state flux-contribution profiles under a range of imposed current densities. These results show that the migration contributions generally dominate. Especially under fuel-cell polarization, the diffusion and stress contributions have strong gradients near the outer boundary. The behaviors are consistent with those reported by Euser et al. [16,17] for La 0.6 Sr 0.4 Co 0.8 Fe 0.2 O 3−δ (LSCF) oxygen-transport membranes. there is a small proton flux. Upon polarization, the flux profiles have spatial gradients, with the highest fluxes near the hydrogen fuel boundary. However, within a few seconds, the proton-flux profile becomes spatially uniform again at a significantly higher value, driving protons from the fuel side toward the oxygen side.    Figure 7 shows hydrostatic stress profiles within the membrane for different imposed current densities. The stresses in the fuel cell mode (i > 0) are considerably higher than in electrolysis mode (i < 0). The fuel-cell stresses have local maxima within the membrane, slightly exceeding the high fuel-side boundary stresses. Moreover, the fuel-cell stress gradients are high, especially near the oxygen boundary. The curvature of the hydrostaticstress profiles is positive for electrolysis, while it is negative at open circuit and in fuel-cell operation. Under strong electrolysis polarization, there can be a slight local minimum in the steady-state hydrostatic-stress profile. The present model attributes chemical strain to the protons. Comparing Figures 2a and 7 shows that the correlation between the proton concentration profile and the hydrostatic stress profiles is evident. Figure 8 shows profiles of the normal stress components that contribute to the hydrostatic stress. The radial stresses are compressive while hoop and axial stresses are tensile and nearly identical. The hydrostatic stress distribution is similar to hoop and axial stresses, although with lower magnitude since radial stresses are compressive (Figure 7). The stresses in the porous support are negligible, as chemical expansion in the support is omitted. Figure 9 shows the transient response of the hydrostatic stress as a current density of i = 1 A cm −2 is suddenly imposed from an initial open-circuit steady-state condition. The local stresses decrease at early times, with local upward curvature at a time of around one second. By about five seconds, the steady-state hydrostatic-stress profile is achieved. The transient stress profiles behave similar to the transient proton concentrations profiles (Figure 3a).

Stress Profiles
In a practically operating cell, the gas compositions and temperature can vary owing to reagent depletion, dilution, etc. Cell electrochemical performance generally depends on gas-phase concentrations. Local stresses are also affected by gas-phase composition. Figure 10 shows the maximum hydrostatic stress as a function of oxygen concentration for the imposed current density of i = 1 A cm −2 . As the oxygen partial pressure at the cathode boundary increases, the maximum hydrostatic stress decreases considerably. This curve does not depend significantly on the current density. Consequently, maintaining high oxygen concentration could enhance mechanical integrity.

Sensitivity Analysis
There is considerable uncertainty in some of the physical properties for BZY20. Thus, it is interesting to investigate the sensitivity of predicted results to the model's properties and parameters. Figure 11 is a tornado plot that shows the effects of varying individual properties by ±10% on the predicted proton flux under fuel-cell polarization. The nominal model-predicted proton flux is J OH • O = 0.1036 mol s −1 m −2 . Under these conditions, the highest sensitivities were found to be the activation energy of the proton diffusion coefficient  The present model uses best-available properties [4]. However, results of the sensitivity analysis point to needs for measuring, improving, and validating thermodynamic and transport properties. Such research is ongoing, but is a long-term process.

Summary and Conclusions
The performance of BZY20 as a proton-conducting fuel-cell or electrolysis membrane is simulated using an extended Nernst-Planck model. The model considers the onedimensional radial behavior in a long tubular cell. The Nernst-Planck fluxes depend on local gradients of defect concentration, electrostatic potential, and hydrostatic stress. The model predicts transient and steady state profiles of defect concentrations and stress within a thin (10 µm) BZY20 membrane that is supported on a relatively thick porous Ni/BZY20 composite electrode. The present model concentrates on defect transport within the dense membrane, assuming ideal electrodes. In other words, charge-transfer polarization at the electrodes is neglected. The three charged defects (protons, oxygen vacancies, and small polarons) are incorporated into the membrane via equilibrated defect-incorporation reactions.
An important measure of membrane performance is the proton flux. However, as a practical matter, internal stresses are also important because of potential damage mechanisms with brittle ceramic materials. The present model couples the electrochemical performance and chemo-mechanical performance. A sensitivity analysis reveals an ongoing need to validate material-specific thermodynamic, transport, and mechanical properties, thus improving predictive capabilities of the models.