2.1. Governing Equations
The model presented in this work consists of three domains: positive electrode, ion exchange membrane, and negative electrode. The current work utilizes a Nafion type cation exchange membrane. This choice aligns with the conventional approach in the literature and is consistent with the experimental validation. However, newly developed anion exchange membranes show promising results in terms of lower vanadium crossover and cost and enhanced H
+ permeability [
37]. The following assumptions were made in the numerical model:
Stationary conditions;
Incompressible electrolytes;
The fluids were assumed to be completely diluted;
Side reactions were neglected;
Both electrodes and the membrane were considered isothermal;
The properties of the electrodes, electrolyte, and membrane were isotropic;
Changes in the z-direction of the cell were ignored (depth in
Figure 1).
The model was based on works published by Shah et al. [
38] and Knehr et al. [
39]. The conservation of mass of the charged species can be expressed by Equation (4).
Variable
ε refers to the electrode porosity,
ci is the concentration of species
i, and
Si is the source term of the species (listed in
Table 1).
is the charged species flux described by the Nernst–Planck equation, as shown in Equation (5).
The first term accounts for diffusion, where is effective diffusivity. In the second term (migration term), zi represents the charge of species i, ui is the ionic mobility, F is the Faraday constant, and φl is the liquid potential. In the third term (convection), represents electrolyte velocity.
The effective diffusivity,
, was obtained from the Bruggemann correlation as shown in Equation (6).
The ionic mobility
ui was calculated using the Nernst–Einstein equation, as shown in Equation (7), where
R is the universal gas constant and
T is the temperature.
The electrolyte velocity, represented by
in Equation (5) in the convection term, is calculated by means of Darcy’s Law, as shown in Equation (8), where
p is the pressure and
µ is the dynamic viscosity of the electrolyte, as indicated in
Table 2, among the other electrolyte properties.
K represents the porous electrode permeability calculated using the Kozeny–Carman equation, as shown in Equation (9), where
df is the fibre diameter and
kck is the Kozeny–Carman constant.
To fulfil the condition of electroneutrality, Equation (10) was solved for all charged species, except for SO
42−:
To integrate the species balance with electrochemical reactions and the current flow within the electrode during the charge/discharge processes, the charge conservation equation was solved, as shown in Equation (11).
Equation (11) indicates that the electrochemical reaction rate (
iR) directly corresponds to the charges leaving the electrolyte
which in turn equate to the charges entering the electrode
Both the liquid and solid current densities are expressed by Equations (12) and (13).
The term
, which corresponds to the effective conductivity of the porous electrode, was calculated using Equation (14), where
σs is the electrode bulk conductivity, listed in
Table 3, with other parameters related to the electrodes.
An integration with the Butler–Volmer law, which characterizes the electrochemical reactions occurring at the surface of the porous carbon electrode, is utilised to converge the conservation equations. Following this, the electrochemical reaction rate (
iR) was calculated for both electrodes, positive (“+”) and negative (“−”), as shown in Equation (15) and Equation (16), respectively.
The specific surface area of the electrode was represented by
a, α is the charge transfer coefficient, and
η denotes the overpotential.
and
, the exchange current densities, are expressed as shown in Equations (17) and (18), where
k+ and
k− are the reaction rate constants for the positive and negative side, respectively.
Table 4 lists the kinetic parameters used in Equations (17) and (18).
The overpotential was determined for the positive and negative reactions using Equations (19) and (20).
The standard equilibrium potentials
E+ and
E− were obtained using the Nernst equation, as shown in Equations (21) and (22).
The effective conductivity of the membrane
can be modelled as shown in Equation (23), where
denotes the proton effective diffusion coefficient.
2.2. Boundary Conditions
Figure 1 x and
y coordinates were taken as references for the boundary condition description. At
x =
x0, the anode external boundary is set as an electrical ground, that is, the solid potential is equal to zero:
The species fluxes at the top and bottom of the membrane and the external boundaries of the electrodes are zero, aside from the inlets and outlets:
At
y = 0, a boundary was set for the flux entering the cell through the electrodes.
where
Q is the volumetric flow rate,
wcell is the cell width, and
Le is the electrode thickness. Analogously, the electrodes have a pressure outlet at
y =
hcell, and the flux of the species caused by diffusion is neglected.
A constant current density was applied to the external boundary of the cathode using Equation (28).
where
iavg denotes the user-defined current applied to the boundary. The sign of this parameter determines whether the cell is in charge or discharge. This leads to the application of an electrical insulation to the upper and lower boundaries of the membranes and electrodes.
The geometrical dimensions and operating conditions of the cell are listed in
Table 5.
2.4. Neural Network
In the present study, a multilayer perceptron with backpropagation (MLP-BP) is used, which is a multilayer model with hidden layers. In this model, the output
is estimated using Equation (30). The output of each hidden neuron is calculated with the sigmoid function defined in Equation (31), which receives as input the postsynaptical
of each
neuron from the previous layer, calculated with the linear combination defined in Equation (32), where
represents the inputs to the layers,
the weights of the layers, and
the biases.
The commercial software MATLAB 2022a [
45], commercial code with its Deep Learning toolbox [
46], was used to design and train the ANN.
The number of hidden layers and neurones in each network depended on the complexity of the magnitude considered. Hence, two different networks were trained for the two targeted ANN models: one for voltage prediction, as shown in
Figure 2, and the other for overpotential prediction, as shown in
Figure 3. The training data were split into 70% training, 20% validation, and 10% testing for both networks.
The first network (ANN1) predicts the voltage of the cell under charge and discharge regimes for different charge and current densities. ANN1 is a relatively simple ANN with a single input layer comprising three parameters: SoC, Current density, Charging or Discharging, a single hidden layer with eight nodes, and an output layer with a single node.
Figure 2 shows the architecture of ANN1.
The second network (ANN2) aims to predict the spatial overpotential of the cell for a constant State of Charge of 50% which is expected to be substantially less continuous compared to ANN1. Hence, added complexity is required in ANN2 which contains three hidden layers (with 8, 16, and 8 nodes) and an output layer with a single node for the overpotential.
Figure 3 shows the architecture of the ANN2.
The ANNs are benchmarked using the standard Pearson product-moment correlation coefficient R.