Hydrodynamic and electrochemical modeling of vanadium redox ﬂow battery

– Two and three dimensional modeling of a single cell of vanadium redox ﬂow battery has been done thoroughly according to electrochemical and ﬂuid mechanic equations in this study. The modeling has been done in stationary state and its results have been presented in three chemical, electrical and mechanical sub models. The parametric analysis on some of important factors in cell operation demonstrated that increase in electrode and membrane conductivity and electrode porosity contributes to electric potential increase in cells. Also operational temperature increase leads to decrease in cells‘ voltage. Better ﬂuid distribution on the electrode surface area results in better cell operation, therefore the electrolyte ﬂow distribution form in cell has been studied by designing diﬀerent ﬂow frames. Modiﬁed Navier-Stokes equations have been used in these calculations for porous media. The most coverage on electrode surface and low pressure loss had been the best case criteria.


Introduction
Energy storage technologies are going to have a more extensive operation in efficient energy providing systems in near future.Vanadium redox flow battery is one of the most exclusive technologies to be used in renewable energy generators as peak shaving and electric load leveling in energy storage field.
Consequently VRB integration with wind and photovoltaic power plants of which production rate is variable due to environmental, seasonal and climatic conditions could increase grid reliability in additional to electric output stability.Utilizing electric photovoltaic electricity and wind production methods along with VRB energy storage recourses deducts unnecessary development in transmission and distribution systems.Since electrical current is utilized near the location, it decreases electricity transmission from the grid to remote areas and finally deducts excessive expenses and investment in transmission and distribution [1][2][3].
a Corresponding author: a.ozgoli@srbiau.ac.irSome of the advantages of VRB can be mentioned as follow: high energy efficiency, low and controllable losses, not terminating of battery life due to discharge, cyclic operation without decline in energy storage ability, not having chemical declination because of corrosion; possibility of electrolyte provision without oxidation, possibility of probable oxidized electrolyte recovery using chemical/electrochemical processes, unlimited electrolyte lifetime without elimination, possibility of electrolyte recovery for other applications [4][5][6].
Electrolytes do not exist inside the cells and are stored outside the cells, the tanks and their electrolyte solution can be easily replaced so the characteristic of this kind of battery is too flexible.It is found simple to increase the electrolyte amount or exchange it.In addition it is possible to optimize the entered electrolyte to the system considering the resource size if required [7,8].
One of the most important technologies in VRB is that the power of vanadium battery has to do with the whole electrolyte flow on stack electrodes' surface, while the accessible stored energy has to do with charged electrolyte Article published by EDP Sciences volume.Therefore it is possible to increase the stored energy by merely increasing the volume of the charged electrolyte and without changing the size of the stack.
In addition the stacks do not have to be discharged by the same voltage they had been charged.The possibility of discharge of a stack still exists while another stack is being charged [9][10][11][12][13].
The research done by Weber et al. [14] according to Frias-Ferrer and Codina [15,16] studies which demonstrate the relation between the fluid flow and VRB efficiency, can be mentioned in hydrodynamic field.
More precise studies have been about VRB by Ma et al. [17] and Secanell et al. [18].These studies have focused on details such as influence on electrode effective area, fluid flow velocity of electrolyte on battery performance, power and efficiency, battery's lifetime and electrochemical polarization.
The effect of electrolyte flow on performance in the active layer has been noticed in the study done by Jeorisson et al. [3].
On the other hand the importance of uniformity of electrolyte distribution in cell to energy efficiency increase has been presented by Moyabayashi et al. [19].
Although even when the flow is homogenous, some of the effective changes on velocity are noticed on electrodes' surface.This can cause extreme changes in PH along with other factors.
Bengoa et al. [20] and Wragg and Leontaritis [21] demonstrated the importance of fluid flow velocity.The results of these studies indicated that reversed flow zones are created around the cell's inlets.
In addition, Bengoa et al. [20] indicated that the flow distribution is asymmetric due to inlet design in the active zone.Therefore this zone has the least homogeneity rate in the flow.
The mass value coefficient has been measured in cells bearing baffles in variaus geometries by Wragg and Leontaritis [21].These baffles have been utilized in order to uniform the flow distribution.
Measurement and simulation of the fluid flow distribution are the most difficult parameters in VRBs.These parameters are affected by various factors such as compact configuration of battery, materials used in electrode and acidic characteristic of electrolyte.Also the review studies done in this field have provided the best results in homogeneous distribution in the battery and its effect on efficiency [22].
Gonzalez et al. [23] have studied the adjustment of hydrodynamic redox cells modeling by the CFD model and measured results of the experimental case.The way of fluid flow distribution in order to determine its effective factors has been one of the most important purposes in this study.
Modeling of a VRB cell and the effect of fluid flow rate in inlets on performance were noticed in the study performed by Shah et al. [24].
The determined mathematical framework was used in formulating purposed modeling for the steady state model by You et al. [25,26].Determination of the current density and state of charge in influence on battery efficiency have been carried out in this study.
The two dimensional model by Vynnycky [27] has presented the VRB fluid mechanics and electrochemistry in a simple model.Reducing the modeling complication and presenting approaches of VRB modeling by finite element method have been its most crucial purpose in larger scales.
Kazacos [28] demonstrated that the temperature of electrolyte solution can effectively change the VRB efficiency in another study.The results on this study showed that increase in flow rate can raise the stack and electrolyte temperature.
Mechanical and electrochemical modeling of a VRB using a CFD simulation is the main purpose of this study.So that considering the governing equation the chemical reactions within the cell, the potential rate and obtained electricity current and their changes during the process have been presented.On the other hand carrying out fluid mechanic analysis the fluid distribution way in cell and its being compared in various inlet and outlet cell designs have been studied.

Model approach
As shown in Figure 1, two reactions occur on two sides of the membrane simultaneously during the discharge process.While discharging, the electrons detach the anolyte and enter the catholyte part passing the external circuit, while charging the movement of electrons changes direction.At this time reduction in anolyte part and oxidation in catholyte part occur.VRBs are achieved by vanadium ability in case of two different oxidized states.
V 4+ and V 5+ ions are indeed the oxidized vanadium VO 2+ and VO + 2 ions.So the equations of reactions are written as follow [29]: (1)

Modeling assumptions
1.The fluid flow has been considered incompressible.It has been assumed that the fluid density remains unchanged during the operation.
The electrode, electrolyte and membrane physical characteristics have been considered isentropic and homogenous and their changes in different directions have been eliminated.3. Isothermal condition exists in all directions.4. Side reactions such as oxygen and hydrogen changes have not been taken in to account. 5. Effects caused by gravity have been eliminated.6.Effects caused by water permeation to membrane have been eliminated.7. The model has been considered in stationary state.8.The electrode materials have been made of carbon.9.The membrane material has been made of Nafion [25,29].

Governing equations of porous electrode
The conservation equations have been presented here as follow [24]: Species conservation: v denotes the intrinsic fluid velocity vector; P , indicates the fluid pressure; c i , denotes the concentration of species i, i ∈ { V 2+ , V 3+ , VO 2+ , VO + 2 , H + }; S i , denotes the source term for species i; φ s and φ l denote the potential of the solid and liquid phase, respectively; S φ , denotes the source term for charge conservation.
The source terms for species conservation and charge conservation in the positive and negative electrodes have been shown in Table 1.The source terms in the membrane are all set to zero.The conservation equations cannot be solved without coupling the Butler-Volmer law.The Butler-Volmer law explains the electrochemical reactions occurring on the surface of the porous carbon electrode [25,30,31]: j neg and j pos are considered as the transfer current densities for the positive and negative electrodes respectively.The exchange current densities for the positive and negative electrode reactions are [24,32]: The over-potentials for the positive and negative electrode reactions are defined as: where E pos and E neg can be estimated from the relevant Nernst equations [33]: The values of the equilibrium potentials E 0 pos and E 0 neg have been categorized in Table 5.
The local mass transfer coefficient k m has been introduced to describe the effect of species transport between the bulk solution and the liquid-solid interface, thus the local flux at the surface of the positive electrode during charge will be defined as [27,32]: Combining equations ( 14) and ( 15), the vanadium concentrations at the surface of the positive electrode will be [24,25]: where A pos and B pos have the following expressions [26][27][28]: The mass transfer coefficient is approximately calculated by [33]: Similar expressions are applied to the species at the negative electrode as follows: where The local mass transfer coefficient k m,neg has the same expression as equation (20).
The permeability of the porous media is often described by the Carman-Kozeny equation [34]: where k CK is the Carman-Kozeny constant, which depends on the type of media and is used as a fitting parameter and d f is the fiber diameter, ε is the porosity, the value of the specific surface area relates to ε according to the calculation of the specific area in reference [35].Equation below shows the correction of specific surface area "a" in Table 1: Values of porosity and specific surface area are shown in Table 2 where ε 0 and A e are the reference.The effective diffusion coefficient D eff i and the effective conductivity of the porous carbon electrode are calculated by the Bruggemann correction [24] This equation shows the effective conductivity of the electrolyte [25] The following equation shows the effective conductivity in the membrane which has a special expression.Reference [24] indicates the detail of the deducing process where z f is the fixed site charge in the membrane and c f is the fixed charge site concentration.

Boundary conditions
Boundary conditions are required at all boundaries of the computational domains, as well as at internal interfaces require.

Boundary conditions for charge conservation
Since the battery is operated in galvanostatic mode, the flux conditions for potential distribution of the porous electrode are defined as follows (during charge): The charge leaving the solid phase is balanced by the charge entering the electrolyte where I is the applied current density.For discharge the signs are reversed.Therefore, the boundary conditions on the potential distribution for electrolyte during charge are specified as:

Boundary conditions for momentum balance
Where the velocity boundary conditions are used at the inlets, pressure boundary conditions are used at the outlets and on all walls, the no-slip boundary condition is applied for the momentum equations.The detailed expressions are as follows: The inlet velocity can be calculated by the equation below where the outlet pressure P out is usually set to zero: The diffusive fluxes of all species are set to zero where Q is the volumetric flow rate and A is the cross-sectional area.
At the outlets (y = H) [25]: when all the other boundaries are set to walls' meaming that the fluxes are zero:

Porous media fluid mechanic equations
Navier-Stokes equations are required for modeling of flows through a porous medium.The incompressible Navier-Stokes equations can be written as [36,37]: u denotes the velocity vector, P the pressure field, ρ the constant density, μ the dynamic viscosity coefficient and f represents the external body forces acting on the fluid (i.e.gravity).
Although the general form of the Navier-Stokes equation is valid for the flow inside pores of the porous medium, its solution cannot be generalized to describe the flow in porous region.Therefore, the general modification form of Navier-Stokes equation should take place to describe the flow through porous media.To describe the linear relation between the velocity u and gradient of pressure P in the porous medium Darcy's law is used to this aim.It defines the permeability resistance of the flow in a porous media [37]: r is considered as the Darcy's law resistance matrix and u the velocity vector.r is a diagonal matrix with coefficients 1/β, where β is the permeability coefficient in the case of considering a homogeneous porous medium.It is possible to define the Reynolds number associated to the pores: δ is the characteristic pore size.Whereas Darcy law is reliable for values of Re p < 1, otherwise it is required to consider a more general model which accounts also for the inertial effects, such as: C is considered as the inertial resistance matrix.Considering a modified Navier-Stokes equation in the whole domain the momentum equations become as equation (48) including the two source terms associated to the resistance induced by the porous medium (linear Darcy and inertial loss term): A modified Darcy's resistance matrix should be used in simulation software, as follows where C is a diagonal matrix: It should be noted that in laminar flows through porous media, the pressure P is proportional to velocity u and C can be considered zero (r * = r) where X is the identity matrix.Therefore, the Navier-Stokes momentum equations can be rewritten as:

Numerical approach
The VRB calculation model is a mathematical demonstration of mechanical phenomena and also electrochemical processes concerning battery operation.This model consists of a set of partial differential equations along with boundary conditions which explain transferring processes.It also includes relations indicating connections among transferring equations.In addition the equations which explain physical and electrochemical materials used in VRB, are required for modeling.
About explained equation set have been solved by COMSOL Multiphysics r software.The acceptable relative error tolerance has been considered 2.5E-6 here.
The dimensions of the base 3 dimensional cell have been contemplated 35 × 35 × 4 mm for positive and negative carbon felt and 35 × 35 × 0.203 mm for polymer membrane.Also 10 mm height with 3 mm diameter tubes have been considered.
Defined electrode, membrane and electrolyte properties in VRB modeling have been demonstrated in Table 2, Tables 3 and 4 respectively.
In Tables 5 and 6 have been presented defined kinetic parameters and operational parameters respectively in modeling VRB and their amounts based on experimental references.

Results and discussion
Results achieved from the model have been evaluated as follow.

Chemical sub-model
Changes in concentration according to the two dimensional model consisting V 2+ and V 3+ in negative part and V 4+ , V 5+ in positive part have been presented in As noticed in Figures 2 and 3, V 2+ 's concentration is reduced from the inlet to outlet of the cell.The concentration of V 3+ increases.This is evident in battery process.As mentioned before, according to electrochemical battery cell reactions, V 2+ converts in to V 3+ while being discharged in the negative part.The maximum and minimum of V 2+ concentrations are equal to  On the opposite side, the positive side of the battery, V 5+ converts in to V 4+ so the V 4+ concentration increases as shown in Figures 4 and 5

Electrical sub-model
Changes relating electric potential have been shown in Figures 6 and 7.The electric potential reduces from the negative side to the positive side as it is concluded from these figures.This is caused by the process of electricity production in battery cells.This fact is also obvious about electric potential of the membrane because more electric potential is noticed on the side touching the negative electrolyte.
As notified, the electrolyte potential value reduces from the negative to the positive side.This occurs because of electric potential production while the cell is being discharged on the negative side.Maximum amount of electrolyte potential has been calculated as 0.19 V on the negative side.
The modeling results demonstrate that the maximum voltage produced in cell electrode is 1.143 V. Also, the current density has been calculated 105.55 mA.cm −2 at max.Increasing the electric potential in the cell the current density decreases.This fact has to do with cell energy balance of which power production must be in a certain amount.The logical selected level of which the battery operation must not be lowered has been defined as 50 mA.cm−2 .The reason of this fact is deduction in fuel charge concentration and the battery's uselessness.

Parametric study
In order to perform the parametric analysis of the model, various models have been created based on the base model in each of which sensitive parameters have been selected.Their change effect has been studied on results from the model.These parametric analyses have been taken in logical and reasonable ranges.
Electrode conductivity rate in the base model was equal to 66.7 S.m −1 .This parameter is one of the most effective factors because of materials used in battery components.The results achieved in this study have been shown in Figure 8.The effect of increase in electrode conductivity, showed increase in electrode electric potential.Sigma e defines the membrane conductivity (S.m −1 ) in this figure .Change has been performed in membrane conductivity in range between 1 and 19 S.m −1 .The achieved results demonstrate that increase in membrane conductivity more than 10 S.m −1 , doesn't have a significant effect on the whole cell potential, whereas its reduction causes undesirable effects on the whole cell voltage rate.This has been mentioned in Figure 9. Sigma m defines the membrane conductivity (S.m −1 ) in this figure .In this parametric study, the porosity range has been considered between 0.7 and 0.95.The achieved results indicate that increase in porosity, has a desirable effect on the cell electric potential.This has been shown in Figure 9.It is noted that increase in porosity rate, is obtained by selecting a better electrode for the battery cells.This verifies the importance of applying better materials.Epsilon signifies porosity in Figure 10.The temperature in which the system is applicable can be mentioned as one of the most sensitive model parameters.This means that the system works in non-isothermal conditions.The result of the effect of change in temperature on electric potential has been shown in Figure 11.Obviously the undesirable increase in cell temperature is trivial.Configuration constraints must be taken into account whilst determination of temperature changes.
The range of temperature changes for parametric study has been considered between 283.15 and 303.15K in the present model.Also the cell's temperature has been equal to 293.15 K in the base model.The sign T signifies the cell temperature (K) in Figure 11.

Mechanical sub-model
In this sub-model, by considering the equations of fluid dynamics in porous media along with electrochemical equations, displaying streamlines in the battery electrodes and also obtaining the pressure drop due to fluid flow in the cell were performed.The fluid velocity of inlet tubes has been considered 0.006 m.s −1 and pressure loss changes have been achieved in two sides of the cell.In Figure 12, the profile of fluid flow distribution has been demonstrated in the positive side of the cell.The way of electrolyte stream line's distribution in the cell is the significant point in this figure .As it is noticeable, although the flow distribution seems almost complete, it is considered zero in some parts especially in the corners of the cells.

Design of flow frame models
The process of entrance of the fluid to the cell is designated by using the flow frame.Flow frame design deals with kind and number of inlet channels of fluid flow to the cell.This is highly important because better design of using all surfaces of the cell can be achieved as increase of designating the fluid distribution way to cell by various flow frame models.The two factors of fluid velocity and pressure loss must be considered in mentioned channel's design.This mechanical modeling has shown another aspect of itself in this case.The number of flow channels in the frame has been considered 3, as it is noticeable in Figure 13   in the diameter of inlet and outlet tubes in the flow frame has shown increase in its pressure loss.
Fluid flows enter the cell as two perpendicular streams as it is observed in Figure 15  Sensitivity analysis on presented models has been performed according to fluid flow velocity to the battery in this field.The reason is the difference in inlet and outlet diameter of cell in different models.So the value of the fluid velocity changes considering the constant flow rate.Results considering with pressure loss on two sides of the cell and produced electric potential with the same velocity and also modified velocities based on the inlet flow rate have been presented in Table 7.
Best cases can be selected according to acquired results from different flow frame modelings and designed channels in its side.The criteria of this choice are its practical possibility from configurative point of view, better flow distribution, maximum electrolyte coverage of cell surface and low pressure loss between the inlet and outlet of the cell.So model No. 5 has the first selection priority according to presented results in Figures 12 to 16 and also Table 7.On the other hand the produced voltage in this model has been obtained equal to 1.25 V in state of constant flow rate, which is higher than other models.
It should be noted that implementing any channels in the flow frame requires manufacturing cost and its fluid mechanic calculation.Change implementations without this calculation and mere based on assumptions and imitating similar system's pattern aren't reasonable.

Conclusion
Modeling of a VRB has been performed in order to achieve a comprehensive view point of a VRB single cell operation.Thus a couple of 3 and 2 dimensional models have been designed and executed with the same geometries.According to available equations in this field, studies have initially been presented on two chemical and electric sub models.Results indicate that while discharging the cell and because of reduction reaction V 2+ concentration has had the maximum value at the cell inlet and it has decreased to 108.53 mol.m −3 at the cell outlet.
On the other hand V 5+ concentration has been reduced from 156 mol.m −3 to 113.13 mol.m −3 due to oxidation reaction on the negative side of the cell.
Maximum electrode voltage of 1.143 V and current density of 105.55 mA.cm −2 have been achieved in the studied cell.
Parametric studies on some of the crucial factors on VRB cell operation, demonstrated that increase in electrode conductivity from 60 to 85 S.m −1 raises electrode electric potential from 1.1422 to 1.1452 V.This analysis about membrane conductivity indicated that its increase from 1 to 10 S.m −1 is highly effective on electric potential optimization.Its effectiveness afterwards drops significantly.Also the results of this study showed that the influence of increase in porosity from 0.7 to 0.95 range has been trivial on the produced voltage.On the other hand the temperature changes of cell operation can't be selected widely.However the undesirable influence of temperature increase on electric potential has been noticed in this study.
Then the results of mechanical sub model were analyzed focusing on state of electrolyte flow distribution.Since the more electrolyte coverage has on cell surface the better it functions, various models were developed with this purpose.Different designings of flow frame for VRBs contribute to different ways of fluid flow in cells.Therefore 5 models of flow frames were designed and the results of their flow distribution were studied for this purpose.
Model No. 5 which has electric potential of 1.25 V and the pressure loss of 12.64 Pa on the positive side of the battery and in constant flow rate, has been selected as the best case because of suitable electrolyte surface coverage and low pressure loss.

c s 5 c 5
pos (c 4 − c s 4 ) = k pos (c 4 ) αpos,c (c 5 ) pos,a ) F η pos RT − exp − α pos,c F η pos RT (14) Figures 2 to 5. It should be noted that the chemical model includes the following results.This is caused by reactions happened in cell and changes acquired in battery electrolyte's content.

Fig. 6 .
Fig. 6.Electrolyte potential change in the central cross section of the cell.

Fig. 7 .
Fig. 7. Membrane electrolyte potential change in the central cross section of the cell.

Fig. 11 .
Fig. 11.Electrode potential change due to cell temperature changes.
considering model No. 2. The state of fluid flow distribution was studied. 1 mm diameter tubes have been selected in this model.Results indicated that the fluid parallelly flows in the cell and exits in this model.The pressure loss is almost about 12.25 Pa in the positive side of the cell which indicates a high pressure loss in this particular design.Model No. 3 has three inlet and outlet channels similar to model No. 2 as it is seen in Figure 14.The difference is borne in one of the inlet and outlet diameter channels with the other two channels.So that the diameter of the first channel has been considered equal to 1.5 mm and the other two are equal to 0.75 mm diameter.Pressure loss shows some increase in comparison with model No. 2. The pressure loss rate on the positive side of the cell has been calculated about 21.20 Pa.Therefore the effect of change
for model No. 4. The tubes considered for this model have 1.5 mm diameter.Results achieved from this model show that fluid distribution occurs in a way that excepts some of the parts of the corners of the cell, the other parts can be under the electrolyte flow coverage.The pressure loss is low in this model which can be considered as an advantage for model No. 4. The pressure loss amount on the positive side on the cell has been calculated about 9.50 Pa.Perpendicular flows enter and exit the cell from two opposite directions in model No. 5 as shown in Figure 16.Consequently it has shown differences in acquired results in comparison with model No. 4. The results of fluid flow distribution in cell have illustrated two sides of fluid flow which have covered more surface in comparison with model No. 4. The pressure loss has been obtained 8.12 Pa in this case.

Table 1 .
Source terms for species and charge conservation.

Table 2 .
Considered electrode properties in VRB modeling.

Table 3 .
Considered membrane properties in VRB modeling.

Table 4 .
Considered electrolyte properties in VRB modeling.

Table 5 .
Considered kinetic parameters in VRB modeling.

Table 6 .
Considered operational parameters in VRB modeling.

Table 7 .
Comprise the main results of flow frame models.Base model (No. 1) Model No. 2 Model No. 3 Model No. 4 Model No. 5 Negative side pressure loss-same velocity (Pa)