Multiphysics simulation of added carbon particles within fluidized bed anode zinc-electrode

Batteries will continue to encounter the problem of dendrite formation until a suitable solution is identified to address the problem. Dendrite formation can short circuit batteries cells, reduce their life span, voltages and cause mechanical abrasion to the cells. Batteries electrodes are part of the approaches that can be used to address these problems but depending on the fabrication of these electrodes and dimensions. Before fabricating and incorporating a real anode reactor to a fabricated ZnBr2 cell system, it was necessary to model the behaviour with injected carbon particles in between 254 microns to 354 microns and simulate the geometry in COMSOL to observe their interaction with the electrolyte. This study investigates the performances of a designed anode reactor and to observe within the reactor the effect of having a uniform and non-uniform current density distribution before the fabrication, physically charge and incorporating it to the anode-side of ZnBr2 cell system.


Introduction
Currently, numerical modelling in any research is now important. It is used to assess the technical solution of a design at any choice and at a design stage. Furthermore, through numerical modelling, the final results can be approached using mathematical modelling without presenting the prototypes of the fabricated physical models [1,2]. Charged particles within zinc bromine batteries cells systems anode reactors determine the reduction and the re-oxidizing of these charged metallic zinc on the anode feeder electrode. Current density distribution is a major concern in the design of electrochemical cells in relation to the incorporated anode and cathode electrodes and circulating electrolytes [3]. Choosing a suitable flow rate also determines current density distributions within electrochemical reactors [4]. However, non-uniform current density distribution in electrochemical cells can be detrimental to the state of health of batteries cells because the electrode areas are usually subjected to high current density distribution [5,6].
In many cases, electrodes degrade faster when part of them are exposed to high current density [7,8]. Additionally, having a uniform current density distribution begins in the electrolyte before progressing to solid electrodes surfaces and within reactors (electrodes) having either high or low surface area [9]. So many issues required addressing, such as the optimization and utilization of electrocatalysts and having the knowledge regarding current density distribution. Most electrocatalysts are fabricated with expensive noble metals [10][11][12].
Non-uniform consumption and deposition, and unreasonably high overvoltage, can result to energy loss and possibly a detrimental side-reaction, which may be the other effects that one would like to minimize. To achieve a laminar flow in fluidized bed reactors is also another issue that required addressing and particularly the flow regime with the injected particles. However, to successfully address these challenges required designing a promising anode reactor that can be used to tailor most of these challenges [13][14][15][16][17][18]. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Geometry and boundary conditions
In Solidworks 2019 [19], the CAD model was designed with the consideration that the experimented fluidized bed reactor had to be conductive. In addition, the designed specification in the model used were the length of 100 mm, breadth of 130 mm and thickness of 12 mm, respectively. The added flow paths were made to be 7 mm thick to prevent the reactor from encountering any uncontrolled turbulences since laminar flow was expected within the reactor. The red arrow in figure 1 has identified the created flow paths. Details of the presented geometry is given in the ECS journal, on practical development of a ZnBr 2 flow battery with a fluidized bed anode zinc-electrode [20]. The study has presented the result from the investigated three feeder electrodes, namely ANSYS [21] simulation, COMSOL [22] simulation and laboratory experiment [20] before selecting the most conductive feeder electrodes and the hydrodynamic behaviour of the examined glass beads and carbon particles within the fabricated anode-reactor. However, by further modelling the geometry in COMSOL was to observe the particle trajectories. With the particle trajectories, it was possible to identify the escaping particles and particles that were conductive during the lab experiment and those not properly conductive according to the gradient. Furthermore, and to compare the presented ANSYS numerical results in the ECS journal [20]. The designed fluidized bed reactor geometric model in figure 4 was selected because it has a great influence on the established experimental results. In addition, the designed and modelled high surface area of the fluidized bed anode reactor has demonstrated that it could prevent the issue of dendrites formation during the charge and discharge process.
The anode fluidized bed reactor was capable for fast electron transfer and properly mixed the injected carbon particles before the fabricating the anode reactor and incorporating it to the fabricated ZnBr2 cell system for the laboratory experiment [20,23,24]. The fabrication and investigation of the anode reactor has been presented in literature [20], with details the practical development of a ZnBr2 flow battery using a fluidized bed anode zincelectrode.

Global definitions
In COMSOL Multiphysics, more detailed information regarding the geometry statistics has been provided for the first investigated parameters, as shown in tables 1 and 2. The CAD import module and particle tracing module were the products utilised in this model. The chosen parameters presented in the tables were chosen and modified to the fit the reactor's design. It was shaped to achieve a laminar flow within the designed anode reactor. All units were specified in m (length unit) and in degree (angular unit). Based on the mesh statistics for the geometry, the observed number of boundaries within the reactor was 61, number of edges (174), number of vertices (116), and one as the number of domains.

Laminar flow within the reactor
At the Inlet, a Poiseuille flow was specified with an average velocity of 0.4 cm s −1 . At the outlet, a uniform pressure of 101600 Pa (relative to atmosphere) was also specified. The Laminar Flow interface was used to solve for the fluid velocity and pressure, as shown in equation (1): where·μ is the dynamic viscosity (SI unit: kg/(m·s)),·u is the fluid velocity (SI unit: m s −1 ),·ρ is the fluid density (SI unit: kg m −3 ), and·p if the pressure (SI unit: Pa). The particle positions are computed by solving second-order equations of motion for the particle position vector components, following Newton's second law, where·q is the particle position (SI unit: m),·v is the particle velocity (SI unit: m s −1 ),·mp is the particle mass (SI unit: kg), and·Ft is the total force (SI unit: N). The only force is the drag force Ft (SI unit: N). Because the particles diameter was 3.54E-7 (m) and the particle velocity relative to the fluid is not too large, the Stokes drag law is applicable, where·u is the fluid velocity (SI unit: m),·m is the fluid dynamic viscosity (SI unit: m s −1 ),·d p is the particle diameter (SI unit: kg). In addition to the drag force, the optional virtual mass force F vm and pressure gradient force F p on the particle was also be considered. The virtual mass and pressure gradient forces were neglected since the density of the particle phase is much greater than the density of the fluid phase, as is true for solid particles in a liquid. However, these forces might approach the same order of magnitude as the drag force since the particles are in a liquid. There are 3000 particles released. The density of the particles released is normalized according to the magnitude of the fluid velocity at the inlet.
This means that there are more particles released where the inlet velocity magnitude is highest and fewer particles released where the velocity magnitude is low.
The model was solved using a stationary study step, because of these two stages; fluid velocity and pressure. Then the particle trajectories were computed using a time dependent study step. Exerting drag forces on the injected particles within the reactor required using solution from the stationary study and defining the fluid velocity for this purpose, not considering the presence of particles during modelling the fluid because the modelling was a one-way coupling. Such coupling is valid for sparse flows of particles with small volume fraction in the fluid.
Neglecting the impact of the momentum onto the fluid by the particles was essential during the modelling. Some assumptions regarding the implementation of the particle tracing were considered. Such as, considering that particles will not be displacing in the fluid based on their occupying volume, neglecting the interaction among the modelled particles to see that the particles distances and diameter were less. Furthermore, using the particle coordinates and center when solving for each of these particles equation of motion and making sure that these injected particles are not travelling more than the expected level within the reactor.

Physical model
Regarding the reactors physical model, the inertial term (Stokes flow) was neglected. Concerning the discretization of the fluid (interface settings) and compressibility, incompressible flow was considered since electrolyte flowing from the reactor's inlet to the outlet was not expected to be stationary. Detail of the parameters for the physical model are presented in table 3. The density and dynamic viscosity of the fluid were both selected prior to the modelling from the material.

Wall and inlet and outlet
To achieve a more accurate result, no slip was the selected condition for the wall due to its relation to the fluid viscosity effects during the interaction. On no slip wall, as presented in equation (7), u, v=0 with the boundary layer. However, in slip wall, the normal velocity is zero (v=0, u is nonzero) and no boundary condition. Figure 2(a) has further provided the reactors wall and flow path with the full geometry in figure 2(c).
From equation (7), the fluid velocity is u (SI unit: m s −1 ). From table 4, the average inflow velocity (u_av) was based on the shape and dimension of the reactor. The boundary condition was designed to be a bit complicated, but necessary for the flow profile to be fully developed. The computational fluid dynamics (CFD), has a special laminar inflow boundary condition to ensure a fully developed flow profile at the inlet. However, it is not necessary to enter a complicated expression for the velocity profile, but just for the average velocity or flowrate equation (7) represents the No Slip Wall.
The diameter of the outlet of the reactor was designed as 30 mm to reduce any unwanted pressure condition during the modelling as shown in figure 2(b). Both normal and suppressed back flow were considered to promote encountering laminar flow within the reactor. Regarding the condition of the wall, the freeze option was considered to prevent the injected particles from escaping from the outlet and going beyond the expected fluidization region, as shown in table 5 for the particle properties.

Model validation
The design of the numerical model was validated with experimental model, as presented in literature [20]. As displayed in figure 3, the details of the top mesh, middle left mesh, middle right mesh and bottom middle mesh is  presented. Some mesh analysis was carried out to determine the best mesh type, size and method for the model. The tetrahedron mesh was the explored mesh method on the geometry, before progressing to the simulation. Regarding the mesh statistics, the reactor has a minimum element quality of 0.201, average element quality of 0.6786, mesh generated element of 86571, triangle mesh of 14302, edge element of 1404 and has a vertex element of 116 and a predefined size extremely fine mesh. The first study computational time was a minute and twentyseven seconds (1 min and 27 s) by using an initial damping factor of 0.01, a minimum damping factor of 1.0E-6 and 100 as the maximum number of iteration by making the modelling fully coupled and using a direct linear solver. Four direct solvers were initially considered before choosing the PARDISO option since the PARDISO solver has the tendency to arrive more quickily to solve any finite element problems of well-conditioned and due to a PARDISO solver biggest advantage to address some extremely ill-conditioned difficulties. Apart from the PARDISO solver, all other solver such as SPOOLES, MUMPS also uses LU decomposition. PARDISO solver is the fastest solver compared to SPOOLES and MUMPS. The slowest solver is the SPOOLES. However, all direct solvers consume a lot of RAM to solver simulation problems. Storing a solution out of core can only be possible using a PARDISO and MUMPS solver and offloading onto the hard disk some of the problems. However, MUMPS solver can support computing cluster and allow using typically on a machine more memory [25][26][27].  The Study-2 used 2 min and 57 s with a solving time dependent in the range of (0,0.01,0.3) seconds and a time step of 3.0E-4 (s). The relative tolerance was 1.0E-6. During the study 2, to make the simulation to converge the non-linear method was changed to automatic newton from constant newston to support the Jacobian once per time step.
However, based on the fact that the number of injecting particles would influence on the experimental results, certain steps had to be taken. To determine the appropriate number in the model, as presented in figures 7-9, as obtained from COMSOL, the first consideration is the geometric design of the bed. This was in consideration with what we obtained in the experimental model as depicted in figures 10-12 and detailed in published literature [20].
Secondly, due to the shape of the anode fluidized bed geometry, the allowed maximum number of particles that was be released from the inlet were up 10000, and similar to that obtained in other literature [28].
Lastly, the viscosity of the liquid and the densities of the particles were also put into consideration because they contributed to the behaviour of the experimented injected particles within the anode reactor. As such we had an appropriate guide on the design used to validate the model. The fabrication of the fluidized bed, and its engineering design specification are in line with design practice [20,[28][29][30].

Results and discussion
According to figure 4 at the top, no turbulence was detected at the bottom of the reactor during the numerical simulation, and all injected particles travelled exactly as expected and did not expected the fluidization region. Particles that failed in reaching a preferred level were due to the mass of these particles. None of the injected particles were expected to reach the outlet of the reactor during the numerical modelling and simulation. The idea of using the particle trajectories approach in COMSOL Multiphysics as mentioned earlier was to observe the interaction of these particles with the. The also showed that some of the injected zinc particles were grounded (stacked) in the reactor. A transmission probability defines the ratio number of particles that have made it to the outlet divided by the released particles number. For instances, injecting 3000 particles and seeing 95% of these particles been grounded has shown that a good result was achieved. From figure 5, at the bottom is the slice velocity magnitude used as the cross-sectional surface of the anode-reactor and sometimes on all the geometry to show changes in the specific area of the plot. The radius of the reactor was measured to the 65 mm, while the length of the reactor was 130 mm.
The scalar quantity from the simulation results, as shown in the plot in figure 6(a), were processed using the pressure contour plot. The displayed results in coloured (series) and lines were the contour plots presented in figure 6(b). The contour plot has made it easy to know the encountered stress by the reactor according to the gradient. The mixing of the particles was visualized through plotting a Poincare by placing coloured dots. At where each particle has passed through a cut plane. This approach is used for defining a single multiple or poincare section in parallel planes, as shown in the plot in figure 6(b). The particles final positions were represented by those colours. Particles in red colour are identified by an initial position of x<0 and with blue colour as x>0. The mixing of all injected particles begins as they migrate to the downstream towards the negative direction-y. Within the reactor, particles not properly mixed and completely can also red or blue. Thus, this produces good particle mixing which depends on the applied velocity. The pressure drop at the inlet and out against the velocity are presented in figures 7(a) and (b) and 8(a) and (b) with the fluidization region in figure 9.
According to the particle trajectories, few identified escaping particles were properly conductive during the laboratory experiment according to the following observations after the experiment. During the real laboratory experiment, after the numerical modelling with COMSOL, the deposited charged zinc particles on the added carbon particles were all conductive due to the prepared concentration for the anode and cathode electrolyte solution apart from the charges and discharges rates (amp) [31][32][33][34].  Secondly, the two-electrolyte solution includes 3 moles of ZnBr2, 1 mole of KCl and 1 mole of ZnCl 2 for the anode concentration and the cathode concentration includes 3 moles of KBr and 1 mole of KCl. These chosen and prepared concentrations were considered to prevent having excess salt during charge [35,36].
Thirdly, due to the various masses of the modelled particles, we expected some of these particles to escape during the numerical modelling because of the particle's densities. However, the different masses did not prevent them to not properly charge [37][38][39][40][41][42][43].   Furthermore, during the laboratory experiment as illustrated in figures 10-12, there were considerations for the anodes and cathode locations. Also, during the laboratory experiment, a mesh was tailored to the outlet of the fabricated reactor, which prevented these particles from escaping. Thus, the particles had to fall back to the surface of the reactor to be further charged and fluidized. This method showed to help to control the discharge of the particles in the fluidized bed.

Conclusions
The investigation on the performances of a designed anode reactor has been carried out. The effect of uniform and non-uniform current density distribution has been observed in this model numerically. It was done before the fabrication, then physically charging and incorporating it to the anode-side of ZnBr 2 cell system. Figures 10-12 represent the experimental model used in the study as this as used to validate the numerical model in COMSOL Multiphysics 5.5.
Through the presented results, especially the particle trajectories plotted in figure 4, the positions of all injected particles and their total shear rate have been observed through the colour gradient since these injected particles has mass and not reaching the outlet as expected. Particles crossing the fluid streamline during the modelling were due to particles inertia and their contact with the reactor wall during mixing when they stopped moving because of the applied freeze boundary condition. Particles agglomerates has also made the velocity of the fluid to be extremely slow towards the final time study. The ratio of the injected particles was known through the transmission probability and the released particle total number.
The presented particle trajectories result in this journal has demonstrated the interaction particles within the reactor before it was incorporated to the anode reactor, on the practical development of a ZnBr 2 flow battery with a fluidized bed anode zinc-electrode [38]. Most of the particles were successfully trapped within the reactor by not exceeding the applied minimum fluidization velocity. However, altering the applied velocity could have supported these injected particles to escape and exceed the expected fluidization region. Part of the recommended future work should include modelling and simulating different particles diameter apart from the investigated sizes in this journal and increasing the shape and dimensions of the reactor for result comparison.