Handling Complexity of Semi-solid Redox Flow Battery Operation

Oct 2, 2018 - This work attempts to quantify physical complexity through parameter sensitivity analysis and some basic tools of graph theory. The syst...
0 downloads 0 Views 4MB Size
Download PDF
Article Cite This: J. Phys. Chem. C 2018, 122, 23867−23877

pubs.acs.org/JPCC

Handling Complexity of Semisolid Redox Flow Battery Operation Principles through Mechanistic Simulations Garima Shukla†,‡ and Alejandro A. Franco*,†,‡,§,∥

Downloaded via UNIV OF LOUISIANA AT LAFAYETTE on October 26, 2018 at 05:43:35 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

Laboratoire de Réactivité et Chimie des Solides (LRCS), CNRS UMR 7314, Université de Picardie Jules Verne, HUB de l’Energie, 15 Rue Baudelocque, 80039 Amiens Cedex, France ‡ Réseau sur le Stockage Electrochimique de l’Energie (RS2E), FR CNRS 3459, HUB de l’Energie, 15 Rue Baudelocque, 80039 Amiens Cedex, France § ALISTORE-ERI, European Research Institute, FR CNRS 3104, HUB de l’Energie, 15 Rue Baudelocque, 80039 Amiens Cedex, France ∥ Institut Universitaire de France, 103 Boulevard Saint-Michel, 75005 Paris, France S Supporting Information *

ABSTRACT: In this paper, the semisolid redox flow battery (SSRFB) in static mode is investigated through a kinetic Monte Carlo-based mechanistic model. Electrochemical reactions are said to occur in particle suspensions composed of carbon as conductive additive and silicon as active material, where silicon is known to undergo volume expansion on discharge. The coexistence of different physical phenomena in suspension is not trivial and leads to complex behavior. This work attempts to quantify physical complexity through parameter sensitivity analysis and some basic tools of graph theory. The systematic treatment employed herein not only expands the utility of mechanistic models, but also provides a more comprehensive theoretical understanding of these complex systems which can otherwise only be treated as black boxes. It is concluded that the primary source of complexity of the SSRFB is the competition between multiple phenomena and that quantifying the dynamics between parameters is as important as measuring a specific parameter. A systematic method of studying the dynamics is to compartmentalize the complex system by introduction of mesoscopic parameters that emerge as a result of contributing microscale phenomena.

■

INTRODUCTION Ever since the invention of the semisolid redox flow battery (SSRFB) in 2011,1 interest in multiphase particle suspensions for electrochemical energy storage has taken a new dimension after the initial studies on slurry electrodes in 1971.2 Despite the possibility of adopting materials from the vast number of available portable battery chemistries,3,4 commercialization of SSRFBs remains far away. Historically, new electrochemical technologies have taken up to half a century to reach commercialization.5 This is typically attributed to the complex nature of the underlying processes which necessitates trial-anderror-based optimization. The development of engineering science in the 1960s allowed for mathematical frameworks to systematize the scattered knowledge obtained empirically.5 The engineering design of electrochemical systems published by Newman in 1968 was one offshoot of this interdisciplinary union.6 For the seven years that SSRFBs and related devices like semisolid flow capacitors have been studied, seemingly contradictory experimental findings have been persistent. These observations wherein typically conductivity has been studied as a function of singular or combinations of multiple variables like geometry, operation parameters, and electrode © 2018 American Chemical Society

composition, will be reviewed in future publications. To give a few examples, (i) the electrical resistance observed between the current collector and slurry can vary from negligible influence7 to 40% of the total measured resistance;8 (ii) increase in amount of active material in suspension is shown to both increase7,9,10 and decrease11 conductivity; and (iii) it has been found that conductivity through the depth of the electrode can arise from both stress-bearing conducting carbon networks and electron tunneling across conducting carbon clusters;12−15 however, the precise conditions and degree of relevance of each phenomena are not established. While engineering-based approaches are ideal for optimization of systems where the theoretical relationships between significant and important variables are defined, the bottleneck of development occurs when these relationships are unknown, which is the case for SSRFBs. No mathematical model or tool is deemed appropriate or useful unless the right theoretical connections are made. Therefore, the following work explores the width of Received: July 11, 2018 Revised: September 25, 2018 Published: October 2, 2018 23867

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

through an interconnection of empirical expressions representing each phenomenon. The framework allows any parameter to have either a constant value or a dynamic one resulting from feedbacks from any part of the theoretical network. Thus, a large degree of freedom is achieved to test simplifications of empirical expressions and identify their limitations. While it appears counterintuitive to expect generalizations to emerge from a combination of empirical equations, we assume that empiricism is simply an act of reducing the complex reality. If these reduced parts are put together, much like a puzzle, a general picture will emerge. The present framework, whose algorithm is illustrated in Figure 2, was chosen to lay a foundation for the physical phenomena occurring in the system along with a simplified electrochemical model. It is important to establish a better understanding of the physical aspects such that effects of dynamic particle interactions, which eventually impact both physical and electrochemical behavior, can later be appropriately incorporated. This model can also be classified as an agent-based model (ABM), wherein the particles are the agents that are free to move around and possess dynamic properties like radius, volume, density, state of charge, electronic conductivity, lithium insertion capability, and neighboring particles. ABMs offer significant theoretical freedom by serving as a platform for building hybrid multiscale models,23 thus providing a means to develop new hypotheses.24 These models have been criticized for lack of standardization and while some attempts have been made to build mathematical formalism,24,25 there is still no consensus as to how or whether such models should be validated.26 ABMs can also be classified as white-box models which are conventionally used for obtaining new emergent knowledge in a quantifiable yet holistic way, a notable example is the cellular automata. Strictly speaking, within this work, the stochasticity of motion adds a more black box quality to approximate Brownian motion, thus making this a gray box model. Such mechanistic models attempt to represent the system as close as possible to its physical reality to build a tangible theoretical framework to which mathematical abstractions can be applied to later. More effort needs to be placed into developing ways to merge these mechanistic models with the realm of standard mathematical modeling and this work intends to take a step towards bridging this gap. The standard methods of causal discovery have been based on probabilistic methods employed in graphical modeling approaches like Bayesian networks (BN) or structural equation models (SEM). However, the problem arises in laying the foundation necessary to apply these methods; the hardest step is the creation of the first causal flow diagram27 or formally referred to as model specification in SEM.28 It is challenging to build a theoretical framework with a very high degree of flexibility in a mathematically standardized way. Mathematical tools generally require significant simplifications and in a new system such as this, where reproducible experimental data is scarce, it is difficult to know which simplifications are reasonable and which methods should be chosen; hence the persistent lack of standardization of ABMs. Since methods to reduce complexity of systems are not well developed, this is done instinctively based on prior knowledge at the initial stages of the work. Surprisingly, for many complex systems, most of the established modeling frameworks and procedures are ruled out. For instance, BNs have highly advanced tools to generate new knowledge from prior knowledge but necessitate working with directed acyclic graphs (DAGs). In DAGs interdepen-

available tools and their respective shortcomings to create a hybrid methodology that is capable of handling complexity arising from unknown relationships so as to speed up experimental discovery. The lack of concreteness in the state of art suggests the highly sensitive and poorly understood nature of underlying electrode phenomena, mass transfer, and solution phenomena toward standard control variables. On one hand, the difficulty arises from the more profound dependence of electrochemical performance on rheological properties,8,15−18 revealing the poor understanding of such multiphysics systems. Furthermore, complexity lies in the lack of representative and exhaustive empirical theory to describe particle and suspension dynamics. A large number of empirical equations have been reported to calculate the viscosity of solid−liquid suspensions.19 Such mesoscopic effects, although present, are much diminished in standard battery electrodes, allowing them to be approximated as black boxes in continuum models at the bulk scale. While slurries used to tape-cast solid electrodes are equally complex due to the presence of different types of solid materials and polymers, attempts to generate comprehensive understanding are few,20,21 and not surprisingly, these are new strategies in the domain of modeling. For microscale and macroscale, standard modeling methods like analytical and statistical methods, respectively, exist. However, assumptions and empirical expressions that suit other scales do not necessarily hold at the mesoscale and it is difficult to obtain sufficiently generic experimental evidence for such expressions. This implies a breakdown of classical theoretical structures and a need for a radically new approach where connections between parameters and assumptions can be explored. The basis for this need lies expressly for complex systems where human intuition is unable to predict how different phenomena compete and give rise to unexpected observations, thus showing emergent behavior. Theory. This paper discusses the elaboration of a previously published22 in-house, three-dimensional, kinetic Monte Carlobased framework for semisolid redox flow batteries (Figure 1). The model attempts to capture a single galvanostatic discharge at slow C-rates, in a small volume of static slurry in contact with a current collector. The other electrode is assumed to be metallic lithium, capable of providing sufficient lithium ions such that ion transport can be neglected. The model functions

Figure 1. Schematics representing the battery (on the left) with slurry electrode (in black), separator (in yellow), lithium metal (in gray), and current collectors (in white). The blowup of the small region of the slurry close to the current collector is thus simulated in this work (on the right). The electrons are available at the current collector (in orange), and silicon particles (in gray) and carbon particles (in black) are suspended in a lithium ion containing electrolyte (in green). Particles are assumed to experience Brownian motion, whereas electrons can flow only through particle networks. 23868

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

Figure 2. Schematics of the algorithm and physical equations used in the in-house, python-based mechanistic model. Initialization involves introducing carbon and silicon particles into a three-dimensional matrix of dimensions X (20 units), Y (20 units), and Z (variable) to recreate a suspension of particles in electrolyte based on the required volume fractions. The size of the unit cell is specified in the text for different simulations. Within the main algorithm, particle motion (within the kMC framework) and electrochemical discharge (through simplified Butler Volmer) are said to occur simultaneously. The time is incremented by a time step calculated from kinetic rates of motion, after which both motion and discharge modules are implemented repeatedly until the cutoff potential is reached. For details, see previous publication.22

dencies between variables are directed and the causal flow has topological ordering such that it does not loop back to form a cycle. While there have been attempts to reduce directed cyclic graphs (DCGs) to DAGs,29 not all systems can be satisfactorily modeled as DAGs, for instance, those in social sciences, economics, and incidentally SSRFBs. Open questions at this

juncture include whether causality can be discovered in a DCG at all,30 if the DCG in question can be satisfactorily reduced to a DAG, and what kind of information would be lost or misinterpreted through this assumption. In this work, we attempt to quantify the degree of complexity of the problem through the use of such basic concepts of graph theory. 23869

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

Figure 3. Causal flow diagram capturing the model algorithm was prepared intuitively. Root nodes are specified by boxes with dotted lines and the starting input values for particle size and volume fraction of solid material to be introduced in the slurry. Among the many cycles found within this flow diagram, the main one is highlighted in blue suggesting its DCG nature.

The simple causal flow representation of this model is highlighted in Figure 3 to further elaborate the flow of logic in the algorithm. The root nodes (dotted line boxes) are, namely, temperature, simulation box dimensions, C-rates, initial values of particle radius and volume fraction, and other parameters used for estimation of electric potential. All the other nodes (in solid line boxes) appear to be d-connected despite the presence of collider nodes, since directed edges from both directions originate from nodes that share parents. All of these observations indicate that this is a heavily interdependent parameter system. In the presence of the following cycle (in blue in Figure 3): particle radius → viscosity → diffusion coefficient → kinetic rate → total kinetic rate → time step → electronic charge → state of charge → volume expansion → calculation of potential → particle radius, causal flow can be classified as a DCG. Additionally, the last edge is conditional, in that if the potential crosses the cutoff potential, the simulation stops. As mentioned before, standard methods work with DAGs and not DCGs. Cyclic systems are better described by undirected Markov networks than by BNs which handle directed acyclic systems. Within the given kMC framework, the arrow labeled 3 in Figure 2 marks the end of the algorithm processing the state of the system. After this point, all changes in the parameters are saved and the next step is implemented over the new state of the system, thus adopting a Markov property. The final directed edge looping back from volume expansion → particle radius is taken into account indirectly, breaking the DCG into a DAG, thereby allowing for distinct simulation steps and evolution of time. Multiple cycles still exist within the rest of the causal flow but since they are assumed to undergo changes simultaneously, the increments are made in topological order and stored for update in the next simulation step. Additionally, it is interesting to observe that collider nodes are viscosity, particle diffusion rate, mesostructure, and electric charge provided to a single particle. All these aspects of the model are difficult to quantify and generalize theoretically, which is why only empirical observa-

tions are available. These nodes where inputs from different parameters converge are poorly understood because competing ripple effects from different chains make observations nonintuitive. More detailed analysis of this causal flow diagram and modifications is expected in future work.

■

RESULTS AND DISCUSSION At this early stage of model specification, parameter sensitivity analysis (PSA)31,32 has been identified as a useful tool to configure the model to find suitable input parameter ranges, explore capability of the model, and highlight experimentally relevant parameter interdependencies that emerge. In the future, these sensitivities are expected to lay the groundwork for the application of more established graph theory methods where conditional probabilities, correlations, and covariance structures can be employed which are backed by appropriately obtained experimental data sets. The first part of the results section deals with elaboration and optimization of the modeling framework, followed by a section describing the ability of the model to grasp elements of complexity observed in the published experimental literature. The lack of comprehensive experimental studies in any one material or the material used in this study is not a hindrance in validating the relevance of the model since it will be demonstrated that studying complex dynamics between multiple phenomena is as important as studying a single phenomenon in a highly specific battery assembly. Standardizing Model Configuration and Parametrization. Previous studies in the context of PSA for stochastic simulations indicate that, (i) large step sizes in input parameters are useful to counteract the deleterious effect of stochastic fluctuations and that (ii) statistically reliable sensitivity of output from large perturbations is more valuable than unreliable sensitivities arising from small perturbations.33 With this in mind, a primary sensitivity study with 288 cases is simulated and the raw data is provided in Figure S1 (Supporting Information). The depths of the simulation box 23870

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

Figure 4. Results of primary parameter sensitivity analysis represented in terms of sensitivity of carbon content at the threshold for through-going chains (Ctgc) with respect to perturbations in, (a) slurry depth at 1 μm with perturbations of +0.4 and −0.4 μm, (b) C-rate at 0.1C with perturbations of +0.1 and −0.05, (c) silicon content at 10%/v with perturbations of +4 and −4%/v, and (d) silicon content at 14%/v with perturbations for +4 and −4%/v.

studied are 0.6 μm (6 units), 1.0 μm (10 units), and 1.4 μm (14 units), with unit length 100 nm. The C-rates C/5, C/10, and C/20 are chosen to avoid ionic transport limitations. Since the size of the simulation box is fixed, it is difficult to isolate effects of solid−liquid ratio variation from silicon−carbon ratio variation, yet worthwhile mechanistic information can still be obtained. Silicon volume fraction is varied in increments of 4% between 6 and 18%/v and carbon volume fraction in increments of 2% between 3 and 17%/v. As carbon fraction is increased, the capacity increases up to a threshold for through-going chains (referred to henceforth as Ctgc) and plateaus.34 No aggregation mechanisms are assumed to keep particle chains intact; the discharge networks are formed purely by chance and can be broken in the following time step. Drop in capacity is observed for, (i) faster C-rates where discharge durations are consequently shorter and for (ii) increase in silicon volume fraction in the slurry composition. Both these trends appear to arise from poor particle diffusion kinetics, eventually preventing access to the entirety of the active material.

This data was further processed in terms of sensitivity of Ctgc toward slurry depth, C-rate, and silicon volume fraction, as shown in Figure 4a−d. In this data set, only one nominal value was chosen for slurry depth (Z = 1 μm with perturbation of ±0.4 μm) and C-rate (C-rate = 0.1C with perturbation of + 0.1C and −0.05C), while two were chosen for silicon volume fraction (Si = 10%/v with perturbations of ±4%/v, and Si = 14%/v with perturbations for ±4%/v). On average, perturbations in depth induce the biggest variations in capacity, perturbations in C-rate result in the lowest sensitivities, and perturbations in silicon volume fraction show dependence on the nominal value. The nonlinearity indicates interactions between parameters which is expected from the heavily dconnected nature of the causal flow diagram. However, setting this aspect aside, it can be said that: (i) high sensitivity to slurry depth could be a result of finite size effects, (ii) sensitivity toward C-rate could be dependent on the kinetic parameter used in the simplified Butler Volmer expression used, and (iii) the nonlinearity for perturbations in silicon content could point to interference of conductive carbon 23871

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

Figure 5. Model algorithm can be characterized in terms of (a) variation of standard deviation of capacity in terms of normalized lithium inserted into silicon at the end of discharge, as a function of solid content (silicon + carbon), and depth of slurry, and (b) variation of the ratio of time of discharge simulated and the time taken for the simulation to run to indicate speed of the algorithm.

step which is estimated from the total kinetic rate under the variable step size method of the kMC framework. The individual event rates arise from diffusion coefficients which are calculated through the Stokes−Einstein equation and depend on particle radii and viscosity. Finally, viscosity, through the Thomas empirical expression, depends on volume fraction or the total volume of particles and eventually particle radii again.35,36 Thus, it was found that time step and degree of discreteness of capacity depended heavily on the way particle radii were captured, whether from (i) the size observed on discrete three-dimensional grid or from (ii) intermediate values of radii smaller than the resolution of the grid, which is relatively more continuum. The latter case resulted in much smaller time steps and very low fluctuations in capacity. For the sake of simplicity, at the initial stages of gauging the global and local behaviors of the model, the degree of discreteness is minimized for the following results where only three runs for reproducibility are implemented instead of five. The previously observed sensitivity of Ctgc toward slurry depth is expected to arise from two effects: (i) physically relevant phenomena and (ii) finite size effects from the discrete model. With the degree of discreteness thus minimized, the finite size effects can further be isolated in terms of total system size. Another study, highlighted in Figure 6a, shows the impact of increasing the slurry depth. Here, the unit length taken is 313 nm and depth is varied from 6 to 26 units in steps of 4, with only 3 runs for reproducibility, amounting to 120 simulations. Excluding the first case, which may be too shallow for silicon particles, the nature of the capacity-carbon fraction trend does not change significantly on increasing depth. It sharpens between carbon fractions 9 and 13%/v and becomes nearly superimposable at 11%/v carbon. This suggests that the larger simulation boxes would provide a more uniform environment in terms of particle density throughout the depth and allow the behavior to converge. For immediate purposes, a depth of 20 units is suitable, thus rendering the simulation box a cuboid. Exploring Physical Phenomena. Once discretization and finite size effects are minimized, a study on the impact of slurry depth was conducted to see how phenomenological effects are captured (Figure 7). The depth was varied across values 20, 25,

networks in the presence of the poor electron conducting active material, an observation that has been reported in rheoimpedance studies.16 The nature of the model algorithm can also be characterized through this primary sensitivity study. Fluctuations of the capacity obtained for multiple runs at the same conditions are shown in Figure 5a. Data fluctuation is a prominent feature of stochastic processes that are discrete in both time and space. The trends are studied in terms of solid content to normalize the compositions. Each curve represents a silicon fraction and changes along a single curve correspond to changes in carbon content. It was observed that increasing the carbon content resulted in a more dramatic reduction in fluctuation than increasing the silicon content at the same total solid content. Given the particle size disparity between carbon and silicon, the impact of the initial configuration on the final capacity ought to be less when the population of the smaller, highly mobile, and identical carbon particles is higher. The accuracy of this model is thus dynamic and it is worth identifying the conditions under which the results generated are more reliable. Figure 5b shows the algorithm speed in terms of the ratio as followstime of discharge simulated: time taken for the simulation to run. Here, larger is the time ratio, the more efficient is the model. Increase in slurry depth is shown to decrease the time ratio regardless of silicon content. This probably arises from the higher population of particles and the consequent increase in computational processing. The nonlinearity observed at low carbon fractions for Z = 1.4 μm could arise from the poor network that results in very low capacity and very short discharge time. Thus, a smaller time ratio for larger carbon populations indicates that the model algorithm starts becoming inefficient. Not surprisingly, the trend with capacity fluctuations indicates that precisely these conditions correspond to higher reliability in results. The aforementioned observations lead to question whether discreteness can have various degrees and what could be its origin. The overall electrode capacity is obtained from individual capacities of every silicon particle and silicon’s ability to discharge depends on carbon networks that have access to electrons. Furthermore, the finite package of electrons provided to the silicon is determined by the time 23872

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

lower salt concentrations show sharp decline in conductivity while higher salt concentrations do not show a significant correlation.7 While no explicit evidence is given to support the following explanation, the phenomenon is ascribed to the ability of ions to screen interparticle repulsive forces and promote aggregation sufficient enough for formation of carbon chains but not so much that flocculation is observed.7 Thus, it can be said that at appropriate salt concentrations, sensitivity to channel depth can be as subtle as shown by our model. This suggests that more evidence is necessary to prove that the nonlinear sensitivity of LTO-based slurries toward channel depth is not a feature of insufficient salt concentration. There are no studies showing simultaneous impact of carbon, C-rate, and channel depth on capacity. Given the high dependence of sensitivity of capacity on (i) C-rate and channel depth11 and (ii) carbon content12 reported separately, this model predicts that it is expected that their simultaneous impact on capacity will be highly competitive, system specific, and nonintuitive. The aforementioned cases are further analyzed in terms of distribution of state of discharge of active material particles at the end of discharge, as shown in Figure 8a−d, to provide a mesostructural picture of the slurry. Two distinct behaviors are observed: (i) bimodal distribution with peaks at the extremes and (ii) a distribution of smaller peaks around or trailing a main peak. Case (i) is observed for carbon content below Ctgc wherein conducting networks do not reach all the active material particles, so connected particles discharge completely while a majority remain inactive. For carbon content above Ctgc, case (ii) is observed, indicating the presence of more extensive conducting networks that ensure that nearly all active material particles participate in discharge. Higher capacity observed at 13%/v carbon at C/20 than at 1C is indicated in the shift of the main peaks toward higher states of discharge. Dispersion of peaks is higher for slower C-rates due to longer discharge times and more opportunities for active materials to come in contact with conducting chains. A closer examination of peaks for Z = 40 units showing the highest state of discharge in Figure 8c,d reveals that even though overall capacity (Figure 7) decreases at large depths, it is still possible to have some particles possessing a higher state of discharge. If it were simply due to the greater possibilities of network formation of a larger number of particles in larger simulation boxes, Z = 45 units would have shown an even higher state of discharge. Such a counterintuitive observation shows that complex behavior arises even in such theoretically simplified systems, suggesting that concepts like optimal fractal dimension are worth introducing as future perspectives. A combined small angle neutron scattering, oscillatory rheology, and impedance spectroscopy study concluded that carbon cluster size optimization is critical to achieving lower viscosity while maintaining electronic conductivity.14 Other aspects worth exploring are the assumptions of the electrochemical nature of silicon particles. Initially it was assumed that all silicon particles that are connected to a carbon network originating at the current collector can discharge but never participate in propagating conducting particle networks. Based on density functional theory studies, one can assume that fully delithiated silicon, whether crystalline or amorphous, is expected to be a poor electron conductor due to lack of density of states at the Fermi level; however, this property changes in the presence of inserted lithium.37 To systematically add a layer of theoretical complexity to the system, two parallel studies were conducted to study the impact of these shifting

Figure 6. Sensitivity study with slurry depths varying from 6 units to 26 units in steps of 4 units to minimize finite size effects.

Figure 7. Summary of sensitivity study for simultaneous impact of Crate (C/20 and 1C), carbon content (5 and 13%/v), and slurry depth (Z = 20−45 units in steps of 5 units) on capacity.

30, 35, 40, and 45 units, and a fixed silicon fraction of 14%/v was taken for the following cases: 5%/v carbon at 1C, 13%/v carbon at 1C, 5%/v carbon at C/20, and 5%/v carbon at C/20. The 2 values of carbon fraction are chosen above and below the Ctgc and each case is run 3 times, amounting to 72 simulations. Regardless of C-rate, slurries composed of carbon fractions below Ctgc show poor capacity, while above the Ctgc, the average capacity is higher for C/20. Even though the slope of capacity as a function of depth for both 1C and C/20 is the same, it is important to avoid LiB-based bias and treat SSRFBs as an entirely different technology. One experimental study presents electrode capacity as a function of C-rate, channel depth, and active material concentration such as LTO.11 While the electrode capacity is seen to decrease with increasing Crate, the trend is nonlinear and shows dependence on nominal values of channel depth and LTO concentration.11 The same is observed when increasing channel depth; capacity is seen to fall, but its trend is highly dependent on LTO concentration and C-rate.11 This behavior is simply attributed to poor electrical contact, but the nuances of individual trends are not explained. In another study, where slurry conductivity is studied as a function of salt concentration in an electrolyte, 23873

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

Figure 8. Distribution of state of discharge of silicon particles at the end of discharge for variation of depth between 20 and 45 units in steps of 5: (a) 14%/v Si + 5%/v C at C/20, (b) 14%/v Si + 5%/v C at 1C, (c), 14%/v Si + 13%/v C at C/20, and (d) 14%/v Si + 13%/v C at 1C.

Fermi levels on lithiation rate and electronic conductivity. The first study was performed for two simple cases as shown in Figure 9 considering: (i) one-step lithiation rate, which remains constant regardless of state of discharge, and (ii) two-step lithiation rate, which increases 100-fold for nonzero state of discharge. A system consisting of 13%/v carbon and 14%/v silicon at C-rate of C/20 in a simulation box of depth of 20 units is thus studied. The distribution of state of discharge of silicon particles is observed to be wider for case (ii), even though the final capacity is not very different. The difference in lithiation rate provides certain particles a head start in terms of their proximity to an electron conducting carbon chain. This allows them to discharge more while other particles retain very low or zero state of charge. This study, wherein lithiation rate depends on lithiation extent, illustrates the possibility of using PSA on recursive interdependencies or DCGs to obtain quantified insight. Introduction of the emergent concept of lithiation rate for a specific particle is in itself a significant step toward handling complexity. This allows a vast number of microscopic properties to be bundled together and their resultant can be utilized at the mesoscale. Aside from

Figure 9. Impact of different types of lithiation rates on dispersion of state of discharge of particles. One-step lithiation rate is a constant value, while the two-step lithiation rate increases 100 times once lithiation of a pristine silicon particle has started. 23874

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

Figure 10. Impact of an on-switch (indicated by σ in the legend) that permits silicon particles to participate in network formation, studied at varying silicon volume fractions with (a) 5%/v carbon and (b) 13%/v carbon.

suggesting that not all particles have equal opportunity to discharge even though a majority do. Those with off-switch at 5%/v carbon indicate case (iv), suggesting that most particles are unable to undergo any discharge due to poor network and also offer very poor capacity. In the presence of the on-switch at carbon content 5%/v, it is evident that the added advantage of silicon forming networks and poor carbon networks below Ctgc appear to compete. This behavior results in an absence of a dominating state of charge and a scattering described by case (iii) is observed. This probably arises from the ability of the more spatially isolated silicon particles like those in case (iv) to get discharged under some small probability of contact. The nature of the conducting network that is composed of both carbon and silicon particles is expected to be markedly different from the one with carbon alone. Silicon, being a significantly larger particle and of variable size, is capable of providing numerous points of contact for propagating conducting networks and will probably result in significantly different fractal dimensions, chain lengths, and number of individual chains; all of which will be considered in a prospective work. It has only recently been shown that electron tunneling across particle clusters and not stress-bearing particle chains may dominate electronic conductivity;14 this observation questions a very fundamental assumption that the LiB community tends to apply to SSRFBs. Within the spatial resolution of the model, it is difficult to explicitly capture electron tunneling. However, implicitly, the competition between these two electron-conducting phenomena can be introduced by assuming that conductive chain propagation occurs by a certain probability. The results in Figure 10 provide an elaborate picture of degree of heterogeneity within the suspension in terms of: (i) particle size distribution, (ii) state of discharge charge, and (iii) local particle packing density. Slurries are shown to display both shear thinning and shear thickening behavior,15 where the regime of shear thickening38 is known to be dependent on several particle properties including particle size disparity.39 The state of discharge of particles can also suggest extent of modification of particle surface due to formation of solid electrolyte interface (SEI). It has been observed that SEI modifies particle networks by influencing interparticle

conventionally studied parameters like ion transport, lithiation rate is governed by a massive theoretical web of its own wherein simple concepts like electronic and ionic conductivity of a particle arise from a myriad of interlinked physical and chemical aspects like crystal structure dynamics, existence of grain boundaries, and nature of solid electrolyte interface (SEI). Future development of this model along these lines can provide significant insights because the macroscopic behavior depends on this emergent lithiation rate function which can either (i) increase or decrease, based on dominance of constructive or destructive phenomena; (ii) remain constant with mild fluctuations, indicating that contributing phenomena are equally significant and compete; or (iii) be nonlinear, suggesting the dominating microscale phenomenon shifts temporally. The second study concerned with the electron conductivity of silicon was conducted with an on/off switch for participation of lithiated silicon in conducting networks, as shown in Figure 10a,b, with additional graphs showing final capacity, extent of percolation, and amount of electroactive silicon, in Figure S2 in the Supporting Information. As a reference standard to begin with, it is assumed that both carbon and lithiated silicon show identical electronic conductivity. Silicon contents of 6 and 14%/v are considered at both below and above the Ctgc with carbon contents of 5 and 13%/v respectively, at C/20 in a simulation box of depth 20 units under an on-switch and an off-switch condition; all 8 cases were run 3 times for reproducibility. The condition for on-switch permits high capacities and well-developed conductive networks, regardless of carbon content. The distribution of state of charge of silicon shows clear correlation with the nature of conductive networks. There are four signatures of particle distribution: (i) a single sharp peak, (ii) a sharp peak and smaller trailing peaks, (iii) wide dispersion of peaks across the whole range, and (iv) bimodal distribution with one peak at zero state of charge and one at full discharge. The cases with on-switch at composition 13%/v carbon + 6%/v silicon, 13%/v carbon + 14%/v silicon, and 5%/v carbon + 14%/v silicon fall under category (i), indicating that all the particles are well connected and discharge nearly at the same rates. Cases with off-switch at 13%/v carbon can be classified as case (ii), 23875

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C interaction energies and consequently electron conductivity.8 This is expected to further impact viscoelastic behavior of the slurry, which is sensitive to the presence of aggregates, clusters, or jammed clusters.14 Particle packing fraction determines whether diffusion limited or volumetrically dense clusters are formed; this in turn impacts local electrolyte concentration and consequently ionic conductivity. Studies conducted with reference to impacts of particle nature and volume fraction on conductive and rheological properties indicate existence of highly nonintuitive behavior.1,17,40 It has been observed that ionic conductivity and viscosity are positively correlated and their tendency to increase or decrease temporally depends not only on the dominant type of cluster but also on the time scale of the study.17 While all these aspects have not been considered on-the-fly at this moment and are expected to modify the final result, this model provides another emergent mesoscopic concept, i.e., the degree and quality of heterogeneity in the suspension. As upcoming work will show, this can be used to bridge the gap between the understanding of numerous microscopic phenomena and bulk rheo-electrical behavior.

specific parameter. Finally, once a robust theoretical understanding is achieved and there is enough appropriate and comprehensive experimental evidence to quantify the relationships, a major part of scientific discovery is over. From this point onward, electrochemical engineering principles like those used in LiB models today can be applied. Through the basis of our approach, it is also possible to test and find appropriate physical models and mathematical methods to quantify the specific kind of complexity that each system presents. PSA and Bayesian network methods that have succeeded in demystifying complicated biological systems are inspiring and reassuring for material science. Since the development is especially urgent in materials for energy applications, we aspire to help circumvent the slow pace of scientific discovery that results from trial-and-error to accelerate the lab-to-industry journey for new inventions.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b06642.

■

CONCLUSIONS This work explores a new paradigm for handling complex systems, using semisolid redox flow batteries (SSRFBs) and mechanistic modeling as a basis. At the current level of state of art of knowledge, regardless of what trends are reported for other material systems, when working with a new slurry, it is critical to conduct experimental studies at a variety of nominal values of parameters like active material content, carbon content, salt concentration, C-rate, channel depth, and so on until generalizations emerge regardless of type of cell setup. By considering the entire picture of the theoretical interconnection of concepts and reconsidering commonly accepted assumptions, it is easier to navigate through the experimental results that appear to be contradictory but are not, since they simply cannot be explained by state of the art of relevant physical and electrochemical theory. The main caveat with complex systems is that every problem is unique and may require different procedures to manage the complexity. Through this work, it has been shown that complexity of electrochemical systems has common aspects, such as the existence of competing phenomena, which can be systematically handled and generalized. Simple mathematical tools like PSA and graph theory can be used for studying the interconnections between physical parameters. New theoretical concepts that can provide structure by compartmentalizing and condensing impact of multiple microscopic phenomena have been shown to significantly simplify the problem. Two such emergent parameters were introduced through these studies, namely, lithiation rate of active material and heterogeneity of suspension. This mechanistic modeling framework allows systematic and cautious inclusion of new phenomena; in addition, it is possible to consider theoretical feedback from other parameters so as to reduce degree of empiricism. The beauty of this methodology is that even without fully quantifying each microscopic phenomenon with accuracy relative to other phenomena, it is possible to explain observed bulk properties. The semiquantified nature of mesoscopic parameters can show whether some phenomena are dominating or whether there is competition between the phenomena. In complex systems where parameters are highly interdependent, this information can be more useful than quantifying a

■

Primary sensitivity study of 288 cases; considering an on/off switch for silicon particles to participate in conductive networks and its impact at the end of discharge (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 0033 3 22 82 53 36. Fax: 0033 3 22 82 75 90. ORCID

Alejandro A. Franco: 0000-0001-7362-7849 Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors deeply acknowledge the Région Hauts de France and Fonds européen de développement regional (FEDER) for the financial support through the project “WONDERFUL”. A.A.F. gratefully acknowledges Institut Universitaire de France for its support. Programming contributions by Oscar Xavier Guerrero and helpful discussions at the Laboratoire de Réactivité et Chimie des Solides with Theodosios Framprikis, Vigneshwaran Thangavel, Dr. Matias Quiroga, and Dr. Charles Delacourt are deeply appreciated.

■

ABBREVIATIONS BN, Bayesian networks Ctgc, volumetric carbon fraction at threshold for throughgoing-chains DAG, directed acycling graph DCG, directed cycling graph kMC, kinetic Monte Carlo LiB, lithium ion battery PSA, parameter sensitivity analysis SEM, structural equation models

23876

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877

Article

The Journal of Physical Chemistry C

■

(19) Viscosity of Liquids: Theory, Estimation, Experiment, and Data; Viswanath, D. S., Ed.; Springer: Dordrecht, 2007. (20) Ngandjong, A. C.; Rucci, A.; Maiza, M.; Shukla, G.; VazquezArenas, J.; Franco, A. A. Multiscale Simulation Platform Linking Lithium Ion Battery Electrode Fabrication Process with Performance at the Cell Level. J. Phys. Chem. Lett. 2017, 8, 5966−5972. (21) Liu, Z.; Battaglia, V.; Mukherjee, P. P. Mesoscale Elucidation of the Influence of Mixing Sequence in Electrode Processing. Langmuir 2014, 30, 15102−15113. (22) Shukla, G.; del Olmo Diaz, D.; Thangavel, V.; Franco, A. A. Self-Organization of Electroactive Suspensions in Discharging Slurry Batteries: A Mesoscale Modeling Investigation. ACS Appl. Mater. Interfaces 2017, 9, 17882−17889. (23) Gong, C.; Milberg, O.; Wang, B.; Vicini, P.; Narwal, R.; Roskos, L.; Popel, A. S. A Computational Multiscale Agent-Based Model for Simulating Spatio-Temporal Tumour Immune Response to PD1 and PDL1 Inhibition. J. R. Soc., Interface 2017, 14, No. 20170320. (24) Hinkelmann, F.; Murrugarra, D.; Jarrah, A. S.; Laubenbacher, R. A Mathematical Framework for Agent Based Models of Complex Biological Networks. Bull. Math. Biol. 2011, 73, 1583−1602. (25) Laubenbacher, R.; Jarrah, A.; Mortveit, H.; Ravi, S. A Mathematical Formalism for Agent-Based Modeling, arXiv:0801.0249v1 csMA. arXiv.org e-Print archive. https://arxiv.org/ abs/0801.0249v1, 1 23 (submitted Dec 31, 2007). (26) Fagiolo, G.; Moneta, A.; Windrum, P. A Critical Guide to Empirical Validation of Agent-Based Models in Economics: Methodologies, Procedures, and Open Problems. Comput. Econ. 2007, 30, 195−226. (27) Cooley, W. W. Explanatory Observational Studies. Educ. Res. 1978, 7, 9−15. (28) Schumacker, R. E.; Lomax, R. G. A Beginner’s Guide to Structural Equation Modeling, 3rd ed.; Routledge: New York, 2010. (29) Ariffin, W. N. M.; Salleh, S. The Partitioning Technique of Directed Cyclic Graph for Task Assignment Problem. AIP Conf. Proc. 2016, 1750, No. 020010. (30) Richardson, T. A Discovery Algorithm for Directed Cyclic Graphs, Proceedings of 12th Conference on Uncertainty in Artificial Intelligence, 1996; pp 454−461. (31) Morris, M. D. Factorial Sampling Plans for Preliminary Computational Experiments. Technometrics 1991, 33, 161−174. (32) Francos, A.; Elorza, F. J.; Bouraoui, F.; Bidoglio, G.; Galbiati, L. Sensitivity Analysis of Distributed Environmental Simulation Models: Understanding the Model Behaviour in Hydrological Studies at the Catchment Scale. Reliab. Eng. Syst. Saf. 2003, 79, 205−218. (33) Drews, T. O.; Braatz, R. D.; Alkire, R. C. Parameter Sensitivity Analysis of Monte Carlo Simulations of Copper Electrodeposition with Multiple Additives. J. Electrochem. Soc. 2003, 150, C807−C812. (34) Medalia, A. I. Electrical Conduction in Carbon Black Composites. Rubber Chem. Technol. 1986, 59, 432−454. (35) Edward, J. T. Molecular Volumes and the Stokes-Einstein Equation. J. Chem. Educ. 1970, 47, 261−270. (36) Thomas, D. G. Transport Characteristics of Suspension: VIII. A Note on the Viscosity of Newtonian Suspensions of Uniform Spherical Particles. J. Colloid Sci. 1965, 20, 267−277. (37) Kim, H.; Chou, C.-Y.; Ekerdt, J. G.; Hwang, G. S. Structure and Properties of Li−Si Alloys: A First-Principles Study. J. Phys. Chem. C 2011, 115, 2514−2521. (38) Brown, E.; Forman, N. A.; Orellana, C. S.; Zhang, H.; Maynor, B. W.; Betts, D. E.; DeSimone, J. M.; Jaeger, H. M. Generality of Shear Thickening in Dense Suspensions. Nat. Mater. 2010, 9, 220− 224. (39) Barnes, H. A. Shear-Thickening (“Dilatancy”) in Suspensions of Nonaggregating Solid Particles Dispersed in Newtonian Liquids. J. Rheol. 1989, 33, 329−366. (40) Campos, J. W.; Beidaghi, M.; Hatzell, K. B.; Dennison, C. R.; Musci, B.; Presser, V.; Kumbur, E. C.; Gogotsi, Y. Investigation of Carbon Materials for Use as a Flowable Electrode in Electrochemical Flow Capacitors. Electrochim. Acta 2013, 98, 123−130.

REFERENCES

(1) Duduta, M.; Ho, B.; Wood, V. C.; Limthongkul, P.; Brunini, V. E.; Carter, W. C.; Chiang, Y.-M. Semi-Solid Lithium Rechargeable Flow Battery. Adv. Energy Mater. 2011, 1, 511−516. (2) Hurwitz, H. D.; Srinivasan, S.; Litt, M. Kinetics of a Simple Charge Transfer Reaction at a Slurry Electrode. Ber. Bunsenges. Phys. Chem. 1971, 75, 535−542. (3) Jain, A.; Ong, S. P.; Hautier, G.; Chen, W.; Richards, W. D.; Dacek, S.; Cholia, S.; Gunter, D.; Skinner, D.; Ceder, G.; et al. Commentary: The Materials Project: A Materials Genome Approach to Accelerating Materials Innovation. APL Mater. 2013, 1, No. 011002. (4) Qu, X.; Jain, A.; Rajput, N. N.; Cheng, L.; Zhang, Y.; Ong, S. P.; Brafman, M.; Maginn, E.; Curtiss, L. A.; Persson, K. A. The Electrolyte Genome Project: A Big Data Approach in Battery Materials Discovery. Comput. Mater. Sci. 2015, 103, 56−67. (5) Wendt, H.; Kreysa, G. Electrochemical Engineering: Science and Technology in Chemical and Other Industries, Softcover version of original hardcover ed. 1999.; Springer: Berlin, 2010. (6) Newman, J. Engineering Design of Electrochemical Systems. Ind. Eng. Chem. 1968, 60, 12−27. (7) Dennison, C. R.; Beidaghi, M.; Hatzell, K. B.; Campos, J. W.; Gogotsi, Y.; Kumbur, E. C. Effects of Flow Cell Design on Charge Percolation and Storage in the Carbon Slurry Electrodes of Electrochemical Flow Capacitors. J. Power Sources 2014, 247, 489− 496. (8) Narayanan, A.; Wijnperlé, D.; Mugele, F.; Buchholz, D.; Vaalma, C.; Dou, X.; Passerini, S.; Duits, M. H. G. Influence of Electrochemical Cycling on the Rheo-Impedance of Anolytes for Li-Based Semi Solid Flow Batteries. Electrochim. Acta 2017, 251, 388−395. (9) Hamelet, S.; Tzedakis, T.; Leriche, J.-B.; Sailler, S.; Larcher, D.; Taberna, P.-L.; Simon, P.; Tarascon, J.-M. Non-Aqueous Li-Based Redox Flow Batteries. J. Electrochem. Soc. 2012, 159, A1360−A1367. (10) Li, Z.; Smith, K. C.; Dong, Y.; Baram, N.; Fan, F. Y.; Xie, J.; Limthongkul, P.; Carter, W. C.; Chiang, Y.-M. Aqueous Semi-Solid Flow Cell: Demonstration and Analysis. Phys. Chem. Chem. Phys. 2013, 15, 15833−15839. (11) Madec, L.; Youssry, M.; Cerbelaud, M.; Soudan, P.; Guyomard, D.; Lestriez, B. Electronic vs Ionic Limitations to Electrochemical Performance in Li4Ti5O12-Based Organic Suspensions for LithiumRedox Flow Batteries. J. Electrochem. Soc. 2014, 161, A693−A699. (12) Parant, H.; Muller, G.; Le Mercier, T.; Tarascon, J. M.; Poulin, P.; Colin, A. Flowing Suspensions of Carbon Black with High Electronic Conductivity for Flow Applications: Comparison between Carbons Black and Exhibition of Specific Aggregation of Carbon Particles. Carbon 2017, 119, 10−20. (13) Hatzell, K. B.; Eller, J.; Morelly, S. L.; Tang, M. H.; Alvarez, N. J.; Gogotsi, Y. Direct Observation of Active Material Interactions in Flowable Electrodes Using X-Ray Tomography. Faraday Discuss. 2017, 199, 511−524. (14) Richards, J. J.; Hipp, J. B.; Riley, J. K.; Wagner, N. J.; Butler, P. D. Clustering and Percolation in Suspensions of Carbon Black. Langmuir 2017, 33, 12260−12266. (15) Youssry, M.; Madec, L.; Soudan, P.; Cerbelaud, M.; Guyomard, D.; Lestriez, B. Non-Aqueous Carbon Black Suspensions for LithiumBased Redox Flow Batteries: Rheology and Simultaneous RheoElectrical Behavior. Phys. Chem. Chem. Phys. 2013, 15, 14476−14486. (16) Youssry, M.; Madec, L.; Soudan, P.; Cerbelaud, M.; Guyomard, D.; Lestriez, B. Formulation of Flowable Anolyte for Redox Flow Batteries: Rheo-Electrical Study. J. Power Sources 2015, 274, 424−431. (17) Narayanan, A.; Mugele, F.; Duits, M. H. G. Mechanical History Dependence in Carbon Black Suspensions for Flow Batteries: A Rheo-Impedance Study. Langmuir 2017, 33, 1629−1638. (18) Madec, L.; Youssry, M.; Cerbelaud, M.; Soudan, P.; Guyomard, D.; Lestriez, B. Surfactant for Enhanced Rheological, Electrical, and Electrochemical Performance of Suspensions for Semisolid Redox Flow Batteries and Supercapacitors. ChemPlusChem 2015, 80, 396− 401. 23877

DOI: 10.1021/acs.jpcc.8b06642 J. Phys. Chem. C 2018, 122, 23867−23877