A Surrogate Model of the Butler-Volmer Equation for the Prediction of Thermodynamic Losses of Solid Oxide Fuel Cell Electrode

Abstract

Solid oxide fuel cells are becoming increasingly important in various applications, from households to large-scale power plants. However, these electrochemical energy conversion devices have complex behavior that is difficult to understand and optimize. A numerical simulation is a primary tool for analysis and optimization-design. One of the most significant challenges in this field is improving microscale transport phenomena and electrode reaction models. Two main categories of simulation are black-box and white-box models. The former requires large experimental datasets and lacks physical constraints, while the latter inherits the inaccuracy of typical electrochemical reaction models. Here we show a micro-scale artificial neural network-supported numerical simulation that allows for overcoming those issues. In our research, we substituted one equation in the system, an electrochemical model, with an artificial neural network prediction. The data-driven prediction is constrained and must satisfy all reminded balance equations in the system. The results show that the proposed model can simulate an anode-electrode’s thermodynamic losses with improved accuracy compared with the classical approach. The coefficient of determination $R^2$ for the proposed model was equal to 0.8810 for 800 °C, 0.8720 for 900 °C, and 0.8436 for 1000 °C. The findings open a way for improving the accuracy and computational complexity of electrochemical models in solid oxide fuel cell simulations.

1. Introduction

The solid oxide fuel cell (SOFC) has emerged as a highly promising energy conversion technology, boasting energy conversion efficiency. A single solid oxide fuel cell comprises two porous ceramic electrodes separated by a solid ceramic electrolyte. During operation, hydrogen (or a mixture of hydrogen and carbon monoxide from the reforming process of hydrocarbons) is fed into the anode side, and air is fed into the cathode side. The high operating temperature of the solid oxide fuel cell, typically between 600 and 1000 °C, promotes the transfer of oxygen ions from the cathode to the anode, where they oxidize fuel. This electrochemical reaction releases electrons, creating an electric current that flows through an external circuit and provides electrical power. Therefore, the electrochemical reaction domain comprises the contact line between the electron-conducting, ion-conducting, and pore phases, the so-called Triple-Phase Boundary (TPB). The anode must have a high TPB length density value to sustain the electrochemical reaction. At the same time, phases need to be well connected. A reaction domain where at least one phase is isolated, inactive, and does not contribute to power generation. Therefore, these electrodes feature a complex composite structure that determines the electrochemical performance of the individual cell and also the entire stack of cells [1]. The complex pores structure of an electrode limits the gas transport throughout the gas diffusion layer and is critical to determining the limiting current density of a cell [2]. Numerical simulation has demonstrated its worth as an essential tool for designing and optimizing SOFC electrodes, providing deeper insights into the transport processes across a fuel cell. Outside of design and optimization, numerical simulation provides in-depth insights into the transport phenomena occurring across a solid oxide fuel cell. This includes the flow of fuel and oxidant in channels, ion and electrolyte migration across the electrodes and electrolyte, electrochemical reactions at the reaction domain, and heat generation and transfer from the electrodes. Such insights are crucial for understanding operating principles and how to use them to improve performance.

Black-box, white-box, and gray-box models are different approaches to solid oxide fuel cell simulations that vary in their reliance on experimental data and the incorporation of underlying scientific theories. Black-box models are data-driven and do not lean on an understanding of underlying physics or chemistry. A black-box model uses empirical observations to develop relationships between input and output variables. Previously published studies indicate effective capture of complex relationships between variables, often with a precision exceeding that available for complex computational fluid dynamic analyses [3]. Artificial neural networks (ANNs) are one of the most common black box models used for fuel cell simulation, as evidenced by refs. [4,5,6,7]. The shortcomings of using them include requiring a large amount of experimental data and a lack of physical constraints that can limit predictive capabilities. The possibility of interpreting the data is also limited in this case, as the mechanisms that led to the particular output are unknown.

White-box models are physics-based models built on an understanding of the underlying physics, chemistry, and mathematical balance equations governing the system. These models incorporate fundamental equations, such as conservation laws and reaction kinetics, to simulate solid oxide fuel cell behavior. However, those models have numerous assumptions and simplifications, which can introduce inaccuracies. Various constants in white-box models are obtained by fitting them to the experimental data and inheriting experimental inaccuracies, limiting the application of the model to the operating conditions and system parameters used in the experiment. An excellent example is an electrochemical model for solid oxide fuel cell simulations. The commonly used model is the Butler-Volmer equation; the mathematical relation for the volumetric exchange current density $i_{\mathrm{TPB}}$ ($\mathrm{A}{\mathrm{m}}^{-3}$) is presented in the following form:

$$i_{\mathrm{TPB}}=i_0\left[exp\left(\frac{\alpha nF}{RT}{\eta }_{\mathrm{act}}\right)-exp\left(\frac{-\beta nF}{RT}{\eta }_{\mathrm{act}}\right)\right],$$

(1)

where $i_0$ is the equilibrium exchange current density per unit volume ($\mathrm{A}{\mathrm{m}}^{-3}$). $i_0$ is associated with the available reaction domain, which is quantitatively described by TPB length density and estimated from fitting to empirical data [8], $\alpha$ and $\beta$ are the dimensionless charge transfer coefficients (1) that are obtained by fitting to the experiment [9], n is the number of electrons involved in the electrode reaction, $R=8.3144\mathrm{J}{\mathrm{mol}}^{-1}{\mathrm{K}}^{-1}$ is the universal gas constant, T is the temperature ($\mathrm{K}$), ${\eta }_{\mathrm{act}}$ is the activation overpotential ($\mathrm{V}$), and F = 96,485.3415 $\mathrm{s}\mathrm{A}{\mathrm{mol}}^{-1}$ is the Faraday constant.

The model has several empirical constants and simplifications, including implementing the reaction symmetry coefficient. In many research works, authors assume $\alpha =0.5$ and $\beta =0.5$. This assumption is corroborated by Kazempoor and Braun’s study [10], which is grounded in experimental findings from other research [11,12] involving an electrode examined at temperatures of 750 ${}^{\circ }\mathrm{C}$, 800 ${}^{\circ }\mathrm{C}$, and  850 ${}^{\circ }\mathrm{C}$. It is also a common assumption for an analytical solution of the charge equation, as it reduces exponential terms in the Butler-Volmer model to hyperbolic sine [13,14]. This equation is widely employed in many numerical analyses of solid oxide fuel cells, as confirmed by citations [15,16]. Those values have become so common that some research groups adopt them without supplying a source [17,18], or sources refer to an article that only applies the values, such as in the case of refs. [19,20]. An alternative approach dedicated to SOFCs was introduced by Kawada et al. [21] more than three decades ago. They assert that $\alpha =1$ and $\beta =0.5$ most accurately reflect the experimental data gathered for solid oxide fuel cells. The Butler-Volmer equation, as formulated by Kawada et al. [21], continues to be extensively employed in the numerical simulation of SOFCs, as substantiated by citations [22,23,24]. Furthermore, Kawada et al. [21] demonstrate in supplementary experiments that both sintering temperature and nickel source powder (microstructure) can influence the parameters of the Butler-Volmer equation. In recent work, Gnatowski et al. [25] proposed a computational scheme that employs an artificial neural network to update the charge transfer coefficients based on operational conditions and available datasets. The artificial neural network proposed functional relationships between the charge transfer coefficients and temperature as well as withdrawn current. The study showed that using an artificial neural network to update the electrochemical reaction model improved the accuracy of the predictions made by the model in Solid Oxide Fuel Cell modeling. The approach developed by Gnatowski et al. [25] fell into the category of gray-box models, which combine elements of both black-box and white-box modeling, blending data-driven techniques with knowledge of the underlying physics or chemistry. These models might use machine learning algorithms to approximate certain aspects of the fuel cell system while incorporating fundamental equations and constraints from the white-box approach.

This study aims to build on our previous findings [25,26] and propose a gray-box approach in which the entire electrochemical model of the solid oxide fuel cell anode is replaced by an artificial neural network approximation. At the same time, the neural network prediction is constrained by the fulfillment of charge and mass conservation laws in the system, making the simulation results physically feasible.

2. Mathematical Model

A numerical analysis of the transport phenomena red in the depth direction of an anode was conducted using a system of four differential equations, which are grounded in the fundamental conservation laws of mass and electrical charge. Two different equations for electric potential are incorporated: one addressing the electron current in the nickel phase and the other for ion transfer in the yttrium-stabilized zirconia phase. A diffusion equation is also included for each species in the steam-hydrogen gas mixture that is present within the electrode’s open porosity. Given the thinness of the anode to be 50 $\mathrm{\mu }$$\mathrm{m}$ in this study, we assume isothermal conditions. This notion of minor temperature variation within a single electrode is confirmed by our heat transfer simulations, as evidenced in references [8,27,28]. Furthermore, our previous research has shown that the isothermal model can accurately reproduce the experimental results for an actual SOFC system across a broad range of operational conditions [1]. The complete set of equations is provided below:

$$\left\{\begin{array}{cccc}\hfill i_{\mathrm{TPB}}& \hfill =\nabla ·\left({\sigma }_{\mathrm{ele}}^{\mathrm{eff}}\nabla {\varphi }_{\mathrm{ele}}\right),\hfill & \hfill {\varphi }_{\mathrm{ele}}{|}_0& =\eta ,\frac{{\mathrm{d}\varphi }_{\mathrm{ele}}}{\mathrm{d}x}{|}_L=0\hfill \\ \hfill -i_{\mathrm{TPB}}& \hfill =\nabla ·\left({\sigma }_{\mathrm{ion}}^{\mathrm{eff}}\nabla {\varphi }_{\mathrm{ion}}\right),\hfill & \hfill \frac{{\mathrm{d}\varphi }_{\mathrm{ion}}}{\mathrm{d}x}{|}_0& =0,{\varphi }_{\mathrm{ion}}{|}_L=0\hfill \\ \hfill -{\displaystyle \frac{i_{\mathrm{TPB}}}{2F}}& \hfill =\nabla ·\left({\vartheta }_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{eff}}\nabla p_{{\mathrm{H}}_2\mathrm{O}}\right),\hfill & \hfill p_{{\mathrm{H}}_2\mathrm{O}}{|}_0& \hfill =p_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{bulk}},\frac{{\mathrm{d}p}_{{\mathrm{H}}_2\mathrm{O}}}{\mathrm{d}x}{|}_L=0\hfill \\ \hfill {\displaystyle \frac{i_{\mathrm{TPB}}}{2F}}& \hfill =\nabla ·\left({\vartheta }_{{\mathrm{H}}_2}^{\mathrm{eff}}\nabla p_{{\mathrm{H}}_2}\right),\hfill & \hfill p_{{\mathrm{H}}_2}{|}_0& \hfill =p_{{\mathrm{H}}_2}^{\mathrm{bulk}},\frac{{\mathrm{d}p}_{{\mathrm{H}}_2}}{\mathrm{d}x}{|}_L=0,\hfill \end{array}\right.$$

(2)

where ${\sigma }_{\mathrm{ele}}^{\mathrm{eff}}$ is the effective electron conductivity of the nickel phase ($\mathrm{S}{\mathrm{m}}^{-1}$), ${\sigma }_{\mathrm{ion}}^{\mathrm{eff}}$ is the effective ionic conductivity of yttria-stabilized zirconia ($\mathrm{S}{\mathrm{m}}^{-1}$), ${\vartheta }_{\mathrm{H}2}^{\mathrm{eff}}$ is a function of the effective diffusion coefficients in the porous electrode and is estimated using the bulk diffusion coefficients of hydrogen, the pores’ volume fraction, and the pore phase tortuosity factor, ${\vartheta }_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{eff}}$ is a function of the effective diffusion coefficients in the porous electrode and is estimated using the bulk diffusion coefficients of steam, the pore phase volume fraction, and the pore phase tortuosity factor, ${\varphi }_{\mathrm{ele}}$ is the electric potential of nickel phase ($\mathrm{V}$), ${\varphi }_{\mathrm{ion}}$ is the electric potential of yttria-stabilized zirconia phase ($\mathrm{V}$), $p_{{\mathrm{H}}_2}$ is the partial pressure of hydrogen ($\mathrm{Pa}$), $p_{{\mathrm{H}}_2\mathrm{O}}$ is the partial pressure of steam ($\mathrm{Pa}$), $\eta$ is the total anode’s overpotential ($\mathrm{V}$), and $i_{\mathrm{TPB}}$ is volumetric current density and is estimated by the artificial neural network in this study. Cell’s performance is characterized by $j(T,\eta )$, where j ($\mathrm{A}{\mathrm{m}}^{-2}$) is the current density. This value can be directly measured and used for validation. The current density can be calculated as an integral of $i_{\mathrm{TPB}}$:

$$j={\int }_0^L{i}_{\mathrm{TPB}}\mathrm{d}x,$$

(3)

whereas ionic and electronic currents can be evaluated as:

$$\begin{array}{cc}\hfill j_{\mathrm{ele}}& \hfill =-{\sigma }_{\mathrm{ele}}^{\mathrm{eff}}\nabla {\varphi }_{\mathrm{ele}}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill j_{\mathrm{ion}}& \hfill =-{\sigma }_{\mathrm{ion}}^{\mathrm{eff}}\nabla {\varphi }_{\mathrm{ion}}.\hfill \end{array}$$

(5)

As the experimental data does not contain direct information about values of volumetric current density $i_{\mathrm{TPB}}$, which is a scalar field inside an anode. Therefore, the following design is used: the artificial neural network predicts the value of volumetric current density $i_{\mathrm{TPB}}$ for a given temperature, current density, and position required by the numerical model. The numerical model solves the system of differential equations with the help of the data provided by the ANN. The number of ANN predictions required for one simulation depends on the mesh of the numerical model. Potentials obtained with the numerical model are used to evaluate the discrepancy between the experiment and the model prediction $E_{\mathrm{train},i}$. The shape of volumetric current density $i_{\mathrm{TPB}}$ estimated by the ANN is used to evaluate additional error $E_{\mathrm{aug},i}$ related to a set of constraints with a minimum possible value equal to zero. The process is repeated for each training data point i. An arbitrary number of arbitrary points can be used to ensure constraints on the required domain. The errors $E_{\mathrm{train},i}$ and $E_{\mathrm{aug},i}$ are used to find values of weights and biases. A block diagram of the methodology is presented in Figure 1 and a detailed training process is presented in Algorithm 1.

Detailed artificial neural network training process is described in the Algorithm 1.   

Algorithm 1 Neural Network Training with CMA-ES

Require:  Empirical data set $D=\left\{(T_i^{\exp },j_i^{\exp },{\eta }_i^{\exp })|0\le i\le N\right\}$, CMA-ES optimization parameters $p_{\mathrm{CMA}}$, SOFC’s numerical model M, artificial neural network parameters $p_{\mathrm{ANN}}$ 1: $f_{\min }=\infty$ 2: $\left(D_{\mathrm{train}},D_{\mathrm{test}}\right)=\mathrm{divideData}\left(D\right)$// Divide data to training and test sets 3: $D_{\mathrm{aug}}=\mathrm{augmentedDataset}\left(\right)$// Include arbitrary data points within the bounds of experimental data set 4: CMAoptimizer = initializeCMA $\left(p_{\mathrm{CMA}}\right)$ 5: repeat 6:    nextGeneration(CMAoptimizer) 7:    $P=\mathrm{getPopulation}\left(\mathrm{CMAoptimizer}\right)$ // Get current population 8:    for $c_i\in P$ do 9:      $I=\mathrm{network}\left(c_i,p_{\mathrm{ANN}}\right)$ 10:      $IM=\mathrm{integrate}(M,I)$ // integrate the SOFC’s numerical model and the ANN 11:      $E_{\mathrm{train}}=averageSquaredError(IM,D_{\mathrm{train}})$ 12:      $E_{\mathrm{aug}}=averageSquaredError(IM,D_{\mathrm{aug}})$ 13:      $E_{\mathrm{CMA}}=0.5E_{\mathrm{train}}+0.5E_{\mathrm{aug}}$ 14:      setCandidateFitness($c_i$, $E_{\mathrm{CMA}}$) 15:      $E_{\mathrm{test}}=averageSquaredError(IM,D_{\mathrm{test}})$ 16:      $E_{\mathrm{global}}=E_{\mathrm{CMA}}^2+E_{\mathrm{test}}^2$ 17:      if $f_{\min }>E_{\mathrm{global}}$ then 18:         $f_{\min }\leftarrow E_{\mathrm{global}}$ 19:         save($c_i$) 20:      end if 21:    end for 22: until !stopCondition(CMAoptimizer)

The experimental data is related to a power generation experiment conducted on a particular anode-electrode at various temperatures and current densities. Therefore, only those boundaries are supplied to the neural network. It is important to note that, theoretically, microstructural parameters and fuel composition can also be included as artificial neural network inputs. However, in order to do that, experimental data covering the range of desired inputs must be additionally provided. Providing a suitable dataset covering even tens of different microstructures is time-consuming and challenging. The workflow would include manufacturing cells using different anode slurry compositions, sintering, electrochemical testing, and analyzing microstructure using athe combination of focused ion-beam scanning electron microscopes and three-dimensional reconstruction techniques. For the latter, a single microstructure might take as much as one month of operator work time. Alternatively, data points from the three-dimensional heterogeneous numerical model can be used for network training. Having said that, in this research, we focus on the ability of the integrated model to provide higher accuracy than the classical Butler-Volmer model.

The architecture of the ANN was determined by employing a trial-and-error approach with various configurations, ultimately refining our model to produce satisfactory results. This final architecture comprises three input neurons corresponding to temperature T, overpotential $\eta$, and position within the electrode x. The ANN features two hidden layers, each containing three neurons and a single output neuron corresponding to the exchange current density $i_{\mathrm{TPB}}$. Combinations with two and three hidden layers with three or five neurons were tested as alternative architectures. The Bent Identity function was chosen as the activation function for all neurons as its continuous nature helps prevent abrupt jumps or discontinuities in the gradient, resulting in more stable learning and convergence. In comparison to other activation functions (Softplus, hyperbolic tangent, sigmoid function, arctan), the Bent Identity function demonstrated superior performance in our tests, contributing to the overall efficacy of our chosen neural network model in conjunction with the CMA evolution algorithm. Inputs are scaled to the $[0,1]$ range. The final layer is scaled by 1 × 10⁹ due to the scale of exchange current density values. The error function consists of a linear combination of seven error terms:

$$\begin{array}{cc}\hfill e_1(j,T,\eta )& \hfill =j_{\exp }-j(T,\eta )\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill e_2(j,T,\eta )& \hfill =i_{\mathrm{TPB}}(0,T,\eta )\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill e_3(j,T,\eta )& \hfill =\frac{{\partial i}_{\mathrm{TPB}}}{\partial x}(0,T,\eta )\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill e_4(j,T,\eta )& \hfill =\eta -{\eta }^{\prime }(T,\eta )\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill e_5(j,T,\eta )& \hfill =\frac{1}{M}\sum _i^M\frac{{{\partial }^{2}i}_{\mathrm{TPB}}}{\partial x^2}(x_i,T,\eta )\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill e_6(j,T,\eta )& \hfill =\frac{1}{M}\sum _i^M{j}_{\exp }-j(T,\eta )\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill e_7(j,T,\eta )& \hfill =\frac{1}{M}\sum _i^M\mathcal{H}\left(\frac{{\partial i}_{\mathrm{TPB}}}{\partial x}(x_i,T,\eta )\right),\hfill \end{array}$$

(12)

where M is the number of sampling points (50 in this study). Error terms, in the above order, control for current density discrepancy Equation (6), exchange current density at the electrode’s interface Equation (7), exchange current density derivative at the electrode’s interface Equation (8), overpotential discrepancy Equation (9), the convex shape of exchange current density Equation (10), and the sum of electronic and ionic currents across the electrode Equation (11), and the monotonic shape of exchange current density, Equation (12). For the augmented data set, errors $e_1$ and $e_6$ are not evaluated. The function for optimization takes the form:

$$E=\sum _k{\alpha }_k\sum _{i=1}^n\sqrt{e_k^2(j_i,T_i,{\eta }_i)/n},$$

(13)

where $j_i,T_i,{\eta }_i$ with $i=1\dots n$ are data points from corresponding data set and ${\alpha }_k$ are scaling factors for error terms. Results presented in this work were obtained with setting the ${\alpha }_i$ values as presented in

Table 1. Values of scaling factors in Equation (13).

${\mathit{\alpha }}_1$	${\mathit{\alpha }}_2$	${\mathit{\alpha }}_3$	${\mathit{\alpha }}_4$	${\mathit{\alpha }}_5$	${\mathit{\alpha }}_6$	${\mathit{\alpha }}_7$
$1\times {10}^{-4}$	$2\times {10}^{-8}$	$1\times {10}^{-13}$	20	$1\times {10}^{-18}$	$1\times {10}^{-3}$	$1\times {10}^{-12}$

.

The train set contains 14 points for 800 ${}^{\circ }\mathrm{C}$ and 1000 ${}^{\circ }\mathrm{C}$, seven points each, and the test set was made from the remaining seven points for 900 ${}^{\circ }\mathrm{C}$. Test and training sets include points with $j=0\mathrm{A}{\mathrm{m}}^{-2}$. The augmented data set consists of arbitrary points within the bounds of experimental points, with $T_i\in$ (850, 850, 850, 900, 900, 900, 950, 950, 950) $\mathrm{K}$ and ${\eta }_i\in$ (0.05, 0.1, 0.15, 0.04, 0.06, 0.08, 0.02, 0.04, 0.06) V. Augmented points do not include corresponding current densities. The reason to include them is to ensure the required shape of exchange current density in the whole domain.

3. Results

To investigate the response of overpotential to current density, we conducted a comparative analysis between the computational results from the present study, our previous study [25], and experimental data taken from the open literature [29].

Figure 2 illustrates the juxtaposition of our models’ predictions vs. experimental data, representing experiments conducted at 800 ${}^{\circ }\mathrm{C}$, 900 ${}^{\circ }\mathrm{C}$, and 1000 ${}^{\circ }\mathrm{C}$. The abscissa axis represents the anode overpotential measured for a particular electrode at a given current density. The ordinate axis represents the predicted anode overpotential for a given current density from the experiment. Therefore, the points situated on the continuous line present a perfect fit between empirical observation and prediction. As can be seen from Figure 2, most of the points are situated in the vicinity of the line, indicating fair agreement between calculated and measured data. To compare predictions from this study and our previous model [25], Pearson correlation coefficients $r_i$ and coefficients of determination $R_i^2$ were calculated for each temperature and are juxtaposed in

Table 2. Pearson correlation coefficients $r_i$ and coefficients of determination $R_i^2$ for the data presented in Figure 2. Index N corresponds to this study, while index P to our previous model [25].

T (${}^{\circ }\mathbf{C}$)	${\mathit{r}}_{\mathbf{N}}$ (1)	${\mathit{r}}_{\mathbf{P}}$ (1)	${\mathit{R}}_{\mathbf{N}}^2$ (1)	${\mathit{R}}_{\mathbf{P}}^2$ (1)
800	0.9988	0.9856	0.8810	0.9687
900	0.9897	0.9776	0.8720	0.9191
1000	0.9948	0.9960	0.8436	0.9637

.

In comparison to our previous model, the current solution exhibits lower coefficients of determination along with higher correlation coefficients (except for 1000 ${}^{\circ }\mathrm{C}$, where values differ insignificantly). A higher coefficient of determination signifies better overall agreement with the experimental values of our previous study. On the other hand, higher correlation coefficients indicate that the derivative of the overpotential with respect to the current density is more aligned with the observed experimental data in the present study. This essentially means that the present study offers a more accurate characterization of the trend across all investigated temperatures. It is important to note that the experiment presented in Kishimoto et al. [29] was conducted with the use of a reference electrode located at the electrolyte support. In order to exclude from the measured overpotential the effect of the electrolyte’s resistance and leave only the overpotential of the electrode, the ohmic loss was excluded from the measurements as an anode is a good electrode conductor. Despite this simplification, the measured electrode overpotential is considered the ground truth for the artificial neural network training. The conducted comparative analysis indicates that the proposed hybrid model can properly reproduce the experimental measurements for various temperatures and current densities used in validation.

To explore the relationship between the volumetric current density at the triple-phase boundary and the anode’s depth, we conducted a further examination of the results obtained for validation described in the previous paragraph, specifically at temperatures of 800 ${}^{\circ }\mathrm{C}$, 900 ${}^{\circ }\mathrm{C}$, and 1000 ${}^{\circ }\mathrm{C}$. Figure 3 displays continuous lines in blue, red, and green colors, representing the outcomes of simulations carried out at 800 ${}^{\circ }\mathrm{C}$, 900 ${}^{\circ }\mathrm{C}$, and 1000 ${}^{\circ }\mathrm{C}$, respectively. The left end of the figure corresponds to the anode surface, while the right end represents the anode-electrolyte interface. The electrochemical reaction concentrates near the electrolyte due to the limited ion conductivity of the electrolyte in comparison to the electron conductivity and pace of diffusion at the anode electrode. Upon closer inspection, the trends predicted by the artificial neural network accurately capture the underlying physics and chemistry, aligning with previously published data in terms of shape, order of values, and trends [22,30]. Most notably, the electrochemical reactions occur predominantly in the vicinity of the electrolyte, thus avoiding the typical errors associated with applying the Butler-Volmer (BV) equation. These errors include the assumption that the electrochemical reactions take place throughout the entire electrode, that solution of BV leads to a parabolic course with the minimum of current exchange density in the middle of the electrode, and in some instances, resulting in negative values for exchange-current densities at specific points. The current study successfully overcomes these limitations and errors, demonstrating a more reliable approach to understanding the relationship between volumetric current density and anode depth.

Figure 4 illustrates a distribution of electronic (red line) and ionic (blue line) current densities within the anode’s depth. $x=0\mathrm{m}$ corresponds to the anode surface, while $x=5\times {10}^{-5}$$\mathrm{m}$ represents the anode-electrolyte interface. The calculation was performed with the total anode’s overpotential maintained at $0.1$ V, resulting in an average current density of 1469 A ${\mathrm{m}}^{-2}$. The charge transfer revealed by the rate at which current changes its form from ionic to electronic indicates an electrochemically active layer. Evidently, the electrochemical reactions predicted by the artificial neural network are most intensive at the anode-electrolyte interface and rapidly decrease towards the anode surface direction. This observation suggests that the majority of the electrochemical reaction takes place in close proximity to the anode-electrolyte interface, within approximately 10 $\mathrm{\mu }$$\mathrm{m}$. The active reaction region’s thickness in cermet anodes remains an open research question, but 10 $\mathrm{\mu }$$\mathrm{m}$ is a reasonable approximation for an anode with the considered microstructure and is consistent with the previous observation, both experimental [31] and numerical [24,32]. To summarize, the artificial neural network that substituted the electrochemistry model could correctly represent the macroscopic behavior of the anode in the form of a polarization curve and correctly predict the electrochemically active layer.

Figure 5 illustrates the distribution of electric potential for both ion- and electron-conducting phases as they span the depth of the anode. It is essential to note that the relationship between these potentials is governed by the Butler-Volmer equation in the standard approach. When replacing the model with predictions from an artificial neural network, it is important to verify the physicality of the predictions. A comparison of the results shown in the figure with those from existing literature reveals a high degree of consistency in both values and derivatives at the boundaries of the evaluated intervals [33]. This further reinforces the conclusion that the artificial neural network has successfully captured the underlying physical phenomena within its constraints. These constraints are imposed by the need to satisfy all balance equations within the system.

4. Conclusions

In conclusion, this study demonstrated the development and validation of an artificial neural network-supported numerical simulation in which the Butler-Volmer equation was replaced by a data-driven prediction. It was shown that such a hybrid approach correctly reproduces transport phenomena across an anode-electrode. The accuracy and reliability of capturing the relationship between volumetric current density and overpotential were investigated by comparing the proposed model’s predictions with experimental data at three different temperatures ( 800 ${}^{\circ }\mathrm{C}$, 900 ${}^{\circ }\mathrm{C}$, and 1000 ${}^{\circ }\mathrm{C}$). The proposed surrogate model effectively overcomes some common limitations associated with the use of the Butler-Volmer equation while correctly representing the complex relationship between current density and anode depth in solid oxide fuel cells. Future works should include microstructural features as an input to artificial neural networks. Experimental data featuring a range of microstructural inputs must be acquired to achieve this. Alternatively, data points from the three-dimensional heterogeneous numerical model can be employed for network training. To conclude, the primary objective, which was to assess the capacity of the artificial neural network to replace the Butler-Volmer equation while maintaining physically correct prediction, was met. Focusing on this goal aimed to highlight the advantages of using an artificial neural network-supported simulation to predict the behavior of solid oxide fuel cells under various conditions. For cathodes in Solid Oxide Fuel Cells applications, the oxygen transition from gaseous molecules in the pore phase into ions ${\mathrm{O}}^{2-}$ is also described with a Butler-Volmer equation. The form of the Butler-Volmer equation varies between different composites. For example, for two-phase electrodes, such as an LSCF cathode, the BV form is analogous to the anode form presented in this study. We expect that the results will be similar, and we plan to examine this case next. For three-phase microstructures, reactions occurring at the double-phase boundary and triple-phase boundary can often be described with two Butler-Volmer equations, each corresponding to a different reaction place. In such cases, two neural networks or a neural network with two outputs must be used to model the electrochemical reactions. As learning the artificial neural network in this study is a challenging task, we plan to address the application of the presented methodology to other cases involving different anodes and cathodes in our future research.