A cellular automata framework for porous electrode reconstruction and reaction-diffusion simulation

Jindong Dai ,Chi Zhai ,Jiali Ai ,Guangren Yu ,Haichao Lv,Wei Sun,*,Yongzhong Liu

1 College of Chemical Engineering,Beijing University of Chemical Technology,Beijing 100029,China

2 Faculty of Chemical Engineering,Kunming University of Science and Technology,Kunming 650500,China

3 Department of Chemical Engineering,Xi’an Jiaotong University,Xi’an 710049,China

Keywords:Mathematical modeling Lithium-ion battery Cellular automata Diffusion Electrochemistry Dynamic simulation

ABSTRACT Due to the high charge transfer efficiency compared to that of non-porous materials,porous electrodes with larger surface area and thinner solid pore walls have been widely applied in the lithium-ion battery field.Since the capacity and charge-discharge efficiency of batteries are closely related to the microstructure of porous materials,a conceptually simple and computationally efficient cellular automata (CA)framework is proposed to reconstruct the porous electrode structure and simulate the reactiondiffusion process under the irregular solid-liquid boundary in this work.This framework is consisted of an electrode generating model and a reaction-diffusion model.Electrode structures with specific geometric properties,i.e.,porosity,surface area,size distribution,and eccentricity distribution can be constructed by the electrode generating model.The reaction-diffusion model is exemplified by solving the Fick’s diffusion problem and simulating the cyclic voltammetry (CV) process.The discharging process in the lithium-ion battery are simulated through combining the above two CA models,and the simulation results are consistent with the well-known pseudo-two-dimensional (P2D) model.In addition,a set of electrodes with different microstructures are constructed and their reaction efficiencies are evaluated.The results indicate that there is an optimum combination of porosity and particle size for discharge efficiency.This framework is a promising one for studying the effect of electrode microstructure on battery performance due to its fully synchronous computation way,easy handled boundary conditions,and free of convergence concerns.

1.Introduction

Lithium-ion battery (LIB) has become an important power source for portable devices and electrical vehicles,the performance of lithium-ion batteries is related to the microstructure of the porous electrode.Understanding the relationship between electrode morphology and microscopic transport properties would assist in optimizing electrode structure and overcoming limitations on battery capacity and charging rate.Luetal.[1] verified that LiMn2O4cathode particles with large interface areas and small volumes perform higher specific capacity by a set of experiments.Sivakkumaretal.[2] found that the decrease in particle size of graphite electrodes is helpful in improving the charge efficiency.However,exploring a set of electrodes with different microstructures by experiments is expensive in cost and time,moreover,the morphological features of the particles are stochastic and difficult to regulate due to the multistep fabrication process including mixing,coating,evaporation,and pressing.Therefore,model-based methods are preferred to study the influence of electrode microstructure on battery performance.

Most electrode models are described by partial differential equations (PDEs),which are derived from porous electrode theory and concentrated solution theory,the most widely used one is the pseudo-two-dimensional(P2D)model[3].This model is composed of a one-dimensional macro-scale model to describe concentration variation in electrolyte and a micro-scale model to describe concentration distribution in electrode particles.However,due to the assumption of uniformly dispersed spherical particles,the effect of irregular particle shape and heterogeneous particle distribution on battery performance has not been considered.Some researchers introduced heterogeneous distributed electrode particles into the PDE models to overcome this shortage,which involves two main task,e.g.,generate a broad range of electrode microstructures on the computer and simulate the reaction-diffusion process in these virtual electrodes.

To reconstruct irregular electrode morphology,Westhoffetal.[4] adopted the force-biased algorithm to approximate the distribution of particles by spheres,then a Gaussian random field is used to adjust the shape of every particle to make it irregular.Stephensonetal.[5] proposed a stochastic grid (SG) model closely related to the Monte Carlo technique to reconstruct electrode morphology by cooperative and competitive phase behaviors.These two methods have been verified with the actual electrode,however,the intersection and growth of interfaces in both methods are complex in definition.Wangetal.[6]proposed a simple and practical algorithm called quartet structure generation set(QSGS)to reconstruct porous medium morphology,which has been applied for solid electrolyte layer(SEI)reconstruction[7]and LiCoO2cathode reconstruction[8],but only the porosity is considered to correspond the actual material,more geometry properties are needed to be compared between reconstructed and actual electrodes.Smithetal.[9] adopted a more direct way in which a real battery crosssection is utilized as a simulation background,however,one cross-section is inadequate to describe the geometrical property of the whole electrode due to the stereoscopic inhomogeneous structure [10].

Besides the reconstruction of electrode microstructure,reaction-diffusion processes in the disordered porous media are hard to be simulated by PDEs due to the irregular solid-liquid boundary conditions.A conceptually simple and computationally efficient framework is needed and cellular automata (CA) is a proper alternative.As a useful modeling method,CA has been adopted to study the inherent dynamics of systems and provide evolution mechanisms for complex processes,such as reactiondiffusion [11],infection transmission [12],drug release [13] and many others.Compared with PDEs,CA is more convenient to construct the porous structure and simulate the mass transfer under complex boundary conditions.The implementation of CA models is fully synchronous and free of non-convergence problems during the calculation procedure.

CA models follow the regulation claimed by Wolfram [14],‘‘Simple rules give rise to complex pattern and self-organizing behavior”,which makes CA with a lower complexity of modeling,simulation,and higher computation efficiency.Zubeldiaetal.[15] proposed a CA model to reconstruct soil porous medium and verified by microcomputed tomography,but the availability for porous structures with other scales,e.g.,electrode particles in micron-scale has not been proved.Bandman [16] first applied CA to reconstruct electrode morphology,the influence of CA parameters on reconstruction has been studied but the results have not been verified by actual electrodes.The above researches indicate that CA is an appropriate tool to reconstruct electrode microstructure,but a standardized CA model is absent.

Among many studies adopting CA to simulate the diffusion process,the most common example is the naive diffusion [17],wherein a cell interchanges with its neighbor at a certain probability.Cells in this pattern can only be updated asynchronously,simultaneously update leads to vanish or increase of particles,which violates the conservation of mass.Margolusetal.[18]proposed a partitioning CA model to simulate diffusion synchronously,wherein lattice is divided into 2 × 2 blocks and each block is rotated by π/2 either clockwise or counterclockwise with certain probabilities.Bandman[19] has modeled diffusion and convection by assigning different probabilities to different directions of rotationviathe Margolus scheme.Gurikovetal.[20] adopted the Margolus scheme to study transport phenomena in porous media including diffusion and adsorption.However,the relationship between CA simulation results and actual diffusion coefficient has not been clarified,which is focused in this work.

For the electrochemical reaction problem,CA has been used for simulating lithium dendrite growth [21],adsorption of lithiumions to planar electrode [22],and oxidation-reduction reactions on the electrode surface [23].The above researches are all based on the assumption that liquid bulk is much larger than the electrode,i.e.,the effect of porous structure on reaction has not been considered in available CA models.In addition,many CA models are inapplicable in actual situations due to the absence of a direct relation between model parameters and real-world variables.A systematic identification process for CA parameters’ physical meaning is demonstrated in this paper.

To the best of our knowledge,no effort has been made to develop a CA-based framework to reconstruct electrode microstructure and simulate reaction-diffusion processes in porous electrodes.Ionic diffusion in electrolyte always be described by macroscopic mass transfer equations,and reactions at the electrode-electrolyte interface are often described by microscopic electrochemical equations.The aim of this work is to establish a CA framework that can describe this multiscale process with irregular solid-liquid boundary,which is more convenient and efficient in modeling and computation compared with PDE systems.

The rest of this paper is arranged as follows: a CA model is established to reconstruct the porous electrode with the specific geometric property;next,another CA model with the Margolus scheme is established to simulate the reaction-diffusion process;based on the above two CA models,a complete battery model is established and the discharge simulation results is consistent with the P2D model;then,a set of electrodes with different microstructures are constructed and evaluated for the discharge efficiency;finally,the conclusion is given.

2.Method

2.1.CA model for electrode reconstruction

In general,CA models are composed of four parts,i.e.,grid,neighbor,state,and rule,which can be represented as MCA1={G,N,S,R}.The two-dimensional square grid is adopted in this study and the coordinates of cells are correspondingly set as G={(i,j),i=0,1,2,...,I,j=0,1,2,...,J},whereIis the number of rows of the grid andJis the number of columns of the grid.Moore neighbor is adopted,e.g.,the state of center cell att+1 moment is influenced by the states of 8 surrounding cells and itself attmoment,which can be presented as N={(i-1,j-1),(i-1,j),(i-1,j+1),(i,j-1),(i,j+1),(i+1,j-1),(i+1,j),(i+1,j+1)}.Two possible states are defined for each cell,white areas formed by cells in state S(i,j)=1 correspond to the solid electrode phase,and black areas formed by cells with S(i,j)=0 represent the electrolyte,the state function can be represented by S={0,1}.The transformation rules in this part are derived from Wolfram’s pioneering work,CA models with two-dimension and two states are also called the 2D binary CA.For this kind of CA,transform rules can be represented as Eq.(1),

where S(N(i,j,t))are the states of neighbors of the cell (i,j)at timet,there are 10 possible neighbor situations for one center cell at timet,i.e.,the number of cells with value of 1 ranged from 0 to 9 including the center cell,and there are two possible states(0 or 1)for the center cell at timet+1,thus,there are 210possible transform rules.Each term of Eq.(1) represents a transform mode,and a complete rule includes 10 transform modes,which corresponds to 10 possible neighbor situations.The sum result of Eq.(1)is usually used as the serial number of different rules[24],rather than an actual operation of cells.Wolframetal.[14] has divided these rules into different types,and Bandman[16]indicates that CA that may be successfully used to construct porous medium morphology belongs to class 2 of the Wolfram’s classification,i.e.,the class of so-called ‘‘phase separation CA”,in which white cells and black cells both exhibiting aggregated self-organized properties.

The transition rule adopted in this work are presented in Eq.(2).It should be noted that phase separation CA rules are not limited to Eq.(2),different phase separation CA rules present different separation velocities and same self-organized properties.For each time step,transform rule R is applicated to all cells(i,j)∈G,all cells are updated simultaneously.In addition,cell size is set as 0.05 μm×0.05 μm and the region size is set as 50 μm×50 μm,i.e.,I=J=1000.

To link a CA reconstruction electrode with the actual electrode,four parameters are computed,including porosity,specific surface area,size distribution,and eccentricity distribution.In a twodimensional image,porosity is defined as a fraction of the area of voids over the total area;the specific surface area is defined as a fraction of the total particles perimeter over the total particles area;the size distribution is obtained by the diameter of circles with equal area to each particle,and the eccentricity distribution is calculated by the external ellipse of each particle.

Steps of the reconstruction process are shown in Fig.1.The steps on the left are related to the original images,i.e.,(1) acquire the actual electrode images;(2)binary process for the images,pixels representing solid particles are assigned a value of 1,pixels representing electrolytes are assigned a value of 0;(3) compute the four image parameters;(4) set a tolerance for the parameters between actual images and CA simulation results.The steps on the right are related to CA simulation,i.e.,(1) select a CA rule for simulation;(2) set the domain size;(3) set the initial solid ratioPsand the total time stepstsof evolution;(4) after a simulation,judge the similarity between the simulation results and the actual image,if the parameters are unsatisfied for the tolerance,initial solid ratioPsand evolution numbertswill be adjusted to conduct a new simulation,if the tolerance is satisfied,the CA simulation result will be outputted.

Fig.1.Flow chart of CA reconstruction.

2.2.CA model for diffusion and reaction

The same description MCA2={G,N,S,R}is used to represent the reaction-diffusion CA model,two-dimensional square grid and Moore neighbor are also adopted in this part,G={(i,j),i=0,1,2,...I,j=0,1,2,...,J},N={(i-1,j-1),(i-1,j),(i-1,j+1),(i,j-1),(i,j+1),(i+1,j-1),(i+1,j),(i+1,j+1)}.

A set S={0,1,2,3}contains all possible states of each cell,0 for the cells filled with electrolyte,1 for the cells represented the solid electrode phase,2 for the cells occupied with reducer,3 for the cells occupied with the oxidizer,e.g.,the lithium-ion.It should be noted that cell size needs to be small enough to obtain accurate results,and cell size needs to be large enough to improve computational efficiency.In other words,if the space of one cell can only be filled with one ion,each ionic level motion will be considered to accumulate the micro changes to macro-level change,which is too expensive in computation cost,if the space of one cell is filled with too many ions,the aggregated distribution of ions will be against common sense and increase error.As a mesoscopic simulation method,an appropriate cell size makes the model relate microscopic mechanisms to macroscopic results.After determining the cell size,ionic concentration can be calculated by multiplying cell proportion with a fixed coefficient,i.e.the concentration scale factor κcon.

As shown in Eq.(3),concentration scale factor κconis used to relate cellular automata proportions to physical concentrations and its unit is mol.m-3.The value of κconis determined in the initial of the single diffusion simulation to correspond to an appropriate division value,for example,when the cell number of one column is 1000 and κconis 1 mol.m-3,the smallest change in cell proportion is 1/1000 so that the smallest change in concentration is 0.001 mol.m-3.

Diffusion rules are set according to the Margolus scheme,the grid is divided into 2 × 2 blocks and each block is rotated by π/2 either clockwise or counterclockwise with certain probabilities at the same moment.The rotated rule can be represented as Rrotate={θcw,θccw}:

wherepcwrepresents the clockwise rotation probability,andpccwrepresents the counterclockwise rotation probability.In general,these two probabilities are set equal when convection is not considered,and only the blocks without solid phase are allowed to be rotated.

After rotating each block,the division way of blocks is transformed from even division to odd division.The division rule can be represented as Rdivide={θeven,θodd},for even division,diagonal cells with coordinates that satisfied mod((i+j),2)=0 are sorted into one block,for odd division,diagonal cells with coordinates that satisfied mod((i+j),2)=1 are sorted into one block,mod((i+j),2) represents the remainder of (i+j) divided by 2,these two division ways are shown in Fig.2.At each time step,the grid is divided into even mode and each block is rotated according to the possibility,then the grid is divided into odd mode and all blocks are also rotated.

Fig.2.Even division (left) and odd division (right).

For all reducer and oxidizer cells,if their neighbors contain at least one electrode cell,the electroactive cell changes its state from reducer to oxidizer or vice-versa depending on the probability of reaction.This transition rule can be represented as Rreaction={θox,θre}:

the probability functionspoxandpreare set by an electrochemicalfounded approach,which is obtained by the Bulter-Volmer (BV)equation:

whereitotalis the total current,iaandicare the oxidized current and the reduction current respectively,i0,refis the reference exchange current density,which represents the ability of an electrode to gain and lose electrons,cre,refandcox,refare reference concentration corresponding toi0,ref.The exchanged electrons number corresponding to 1 mole reaction is set as 1 in this work.Fis the Faraday constant(9.6458× 104C.mol-1),Ris the gas constant(8.314 J.mol-1.K-1),Tis the temperature,αaand αcare the charge transfer coefficient,η is the over-potential,creandcoxare the reducer concentration and oxidizer concentration,then the probability functions are set as follows:

wherekcarepresents the ability of an electrode to gain or lose electrons in this CA system,only the reactant cells next to the electrode cells need to be determined whether reaction occurs according to Eq.(7).The generated currentitotalin this discrete model can be calculated as follows:

whereNox-cellandNre-cellrepresent the number of cells that be oxidized and be reduced at one iteration,κmolerepresents the number of moles of one cell filled with oxidizer or reducer,Nais the Avogadro constant (6.02 × 1023mol-1),eis the elementary charge(1.6 × 10-19C),κtimerepresents the actual time corresponding to one iteration,in one iteration,each cell state will be updated according to the transform rule R={Rdivide,Rrotate,Rreaction}.In this paper,all the CA models are realized in Matlab 2020a,and all the PDE models are solved by the PDE tool in Matlab 2020a,the PDE models are used for verification purpose to the proposed CA models.

3.Results and Discussion

3.1.Reconstruction of porous electrodes

The reconstruction results obtained by the phase separation CA model are shown in Fig.3.Fig.3(a) presents the actual electrode morphology of a LiCoO2(LCO) anode [25],which is obtained by tomographic technology and binary process,the image size is 633 × 938 pixels,and one pixel is corresponding to 0.0532 μm.Fig.3(b) presents the actual electrode morphology of a LiMn2O4(LMO) anode that taken from the literature [26],the image size is 543 × 715 pixels,and one pixel is corresponding to 0.278 μm.Fig.3(c) presents the actual morphology of a carbon-silicon (C-Si)composite electrode that taken from the literature [27],the image size is 303 × 340 pixels,and one pixel is corresponding to 0.413 μm.

Fig.3.Actual image (a-c) and reconstructed image (d-i).

Fig.3(d),(e),and (f) show the CA reconstruction results corresponding to Fig.3(a),(b),(c),which are obtained by the phase separation rule Eq.(2).The size of each CA image is 1000 × 1000 cell,and the side length of each cell is 0.05 μm.It can be seen intuitively from reconstruction images that LMO electrode particles are much larger than LCO particles and C-Si particles,which is not obvious in actual images due to the different measuring scales.

Quantitative comparisons between actual and reconstructed images are presented in Tables 1-3.The tolerance of the deviation is set as 0.1,and the deviation calculation equations are presented in Eq.(9),

whereSporerepresents the total area of the pores,Sparticlerepresents the total area of the particles,Lparticlerepresents the total perimeter of the particles,εCAand ε0represent the porosity of CA image and actual image respectively,SVCA(μm-1) andSV0(μm-1) represent the specific surface area of CA image and actual image respectively.The distribution of particle size and eccentricity of the actual and reconstructed images are divided into equidistant intervals,the deviation value of distribution proportion between CA images and actual images in each interval are averaged to calculate the final deviation,di,CAanddi,0represents the particle proportion that equivalent diameter is in theith interval.The calculation of eccentricity distribution deviation is same as that of size distribution,eci,CAandeci,0represents the particle proportion that eccentricity is in theith interval.

Table 1 Comparison between actual (a) and reconstructed (d) images of LCO electrode in Fig.3.

Table 2 Comparison between actual (c) and reconstructed (e) images of LMO electrode in Fig.3.

Table 3 Comparison between actual (e) and reconstructed (f) images of C-Si electrode in Fig.3.

The size distribution and the eccentricity distribution of the actual image and the reconstructed image are depicted in Fig.4 by a histogram form.Fig.4(a),(c),(e)show the equivalent diameter distribution of the LCO,LMO and C-Si electrodes respectively,equivalent diameter means the diameter of a circle with the same area of the electrode particle.Fig.4(b),(d),(f)show the eccentricity distribution of the LCO,LMO and C-Si electrodes respectively.The horizontal axis is divided into equidistant intervals,and the ordinate axis indicates the particle proportion in each interval.It can be seen that there is a similar trend between each actual distribution and the corresponding reconstructed distribution.Although there are relatively large deviations in a few intervals,the mean absolute errors of each image are all less than the set value,i.e.,0.1.

Fig.4.Size distribution and eccentricity distribution in the actual and reconstructed images.

To reach the tolerance,initial solid ratioPsand evolution numbertsare adjusted according to Fig.1.FinallyPsis set as 0.27(d),0.256(e),0.263(f),tsis set as 70(d),315(e),110(f) respectively for the phase separation rule.This algorithm reproduces the fabrication process of electrodes to some extent,i.e.,a random packing of initial solid cells is corresponding to the coating process and the parameterPsis related to the solution concentration,the separation of the solid phase cells and the electrolyte phase cells is corresponding to the solvent evaporation process,the parametertsis related to the drying time.

3.2.Diffusion simulation based on Margolus CA

Firstly,a test was performed to verify that the CA model is reliable to correspond the Fick’s Second law.A simple onedimensional diffusion problem is described by Eq.(10),wherecrepresents the concentration,t(s) represents the physical time,x(μm) represents the diffusion distance,D(m2.s-1) represents the diffusion coefficient.The region length is set as 50 μm,initial concentration distribution is uniformly set as 0.1 along thexdirection.The concentration at the left end of the region is stayed at 1,and the concentration at the right end of the region is stayed at 0.

To correspond to this diffusion problem through the CA discrete model described by Eq.(3) and Eq.(4),the coordinates of cells are set as G={(i,j),i=0,1,2,...,I,j=0,1,2,...,J},whereinIandJare both set as 1000,cell side length is set as 0.05 μm,concentration scale factor κconis set as 1 mol.m-3and time scale factor κtimeis set as 2 × 10-5s,the determining process of the cell length is provided in the Supplementary Material.It should be noted that the value ofIis correlative with the stability of the simulation results,through the preliminary exploration of the model,the selection value ofIis 1000,wherein stability and computational efficiency are both acceptable.In these diffusion problems,there are only two possible states for each cell,i.e.,0 for solvent and 2 for solute.The state of cells with coordinate ofj=1 is stayed at 2,i.e.,S(i,1)=2.The state of cells with the coordinate ofj=1000 is stayed at 0,i.e.,S(i,1000)=0.Ten percent of the cells with coordinate 1 ＜j＜1000 are set asS(i,j)=2,the other cells are set asS(i,j)=0.These settings are to match the boundary condition of Eq.(10).The rotated probabilitypcwandpccware both set as 0.5.Through parameterized scanning of the diffusion coefficient with an interval of 1 × 10-12m2.s-1,it can be found that the diffusion coefficient of 1.25 × 10-10m2.s-1in the PDE model is equivalent to the rotated probability of 0.5 in the CA model,Fig.5 presents that simulation results of these two models are in agreement in space-time distribution under this pair of parameters.

Fig.5.Diffusion simulation of CA model and PDE model.

A series of rotated probabilities are set to simulate this diffusion problem and each corresponding diffusion coefficient is identified.The relation curve of rotated probability and diffusion coefficient is depicted in Fig.6,a parabolic relation is founded betweenpcwandD,which is presented in Eq.(11)

Fig.6.Parabolic relation between pcw and DFick.

the correlation between diffusion coefficients and rotated probabilities can be described by Eq.(11)whenLcell=0.05 μm and κtime=2.5×10-5s,the scope of diffusion coefficient is[1.4×10-11m2.s-1,1.25 × 10-10m2.s-1],in this scope,the mean relative error defined by Eq.(B1)(in Supplementary Material)of the CA simulation results can be maintained lower than 0.1.According to this fitting function, more rotated probability corresponding to different actual diffusion can be found to meet the simulation requirement.

Besides rotated probability,cell lengthLcelland time scale factor κtimeare also influential for the diffusion coefficient,a rotation of one cell with a longer side length is corresponding to a longer displacement,i.e.,a larger diffusion coefficient.In contrast,a larger value of κtimemeans a lesser global evolution number in a unit physical time,i.e.,corresponding to a smaller diffusion coefficient.Fitting results of different cell lengths and time scale factors corresponding to the same diffusion coefficient are presented in Table 4.

Table 4 Fitting results of different Lcell and κtime corresponding to the same diffusion coefficient

As shown in Table 4,when the rotation probability is constant,cell lengthLcelland time scale factor κtimeneed to be decreased or increased at varying rates to correspond to the same diffusion coefficient.To investigate the varying rate of these two parameters corresponding to one diffusion coefficient,Fig.7 is plotted,wherein ordinate is the time scale factor κtimeand abscissa is the square of cell lengthLcell.A linear relation is found in Fig.6,the fitting functions are presented in Eq.(12),

Fig.7.Linear relation between κtime and Lcell.

where the constant term is small enough to be ignored,then the relation of slope and diffusion coefficient can be written as:

the cell lengthLcellis determined by balancing the simulation accuracy and the computational efficiency in a specific simulation,an appropriate cell length can be found easier through Eq.(13).

3.3.Electrochemical reaction simulation

In this part,the BV equation is converted to the reaction probability in the discrete model,which is convenient to be integrated with the diffusion model.The converted method is verified by a cyclic voltammetry process.A schematic diagram of the cyclic voltammetry experiment is shown in Fig.8.The electrode reaction typically involves three steps: the reactant moves to the electrode surface,the electron transfer to the reactant and the reaction occurs,finally the product moves away from the electrode.The working electrodes are immersed in electrolyte with sufficient volume,the effect of another counter electrode can be ignored and the concentration at the right end of the solution is constant.The applied voltage on the working electrode is linearly varying over time and will be overturned after a peak.

Fig.8.Cyclic voltammetry schematic diagram.

A typical cyclic voltammetry curve is shown in Fig.9.There are only reduced substances in the initial solution,as the potential scanning from -0.4 V to 0.4 V at a voltage sweep rate of 1 V.s-1,the electrode current first remains constant and then rapidly increases until the reducer at the electrode surface is substantially oxidized,corresponding to a peak of current,then the current decreases and gradually reaches a steady state,in this stage,the arriving reducer at the electrode is limited by the diffusion rate.The potential scan direction is then reversed,as the potential scanning from 0.4 V to -0.4 V,symmetric reaction process occurs for the oxidizer.More detailed description of cyclic voltammetry can be found in the literature [28].

Fig.9.Typical cyclic voltammetry curve.

A cyclic voltammetry process can be described by PDEs as shown in Eq.(14)and Eq.(15).Eq.(14)presents the diffusion part of this system,the length of the solution regionLxis set as 50 μm,the initial reducer concentrationcre0is uniformly set as 0.2 mol.m-3along thex-direction,and there are only reduced substances in the initial solution,i.e.,the initial oxidizer concentrationcox0is 0.The reducer concentration at the right end of the region is stayed at 0.2 mol.m-3.DoxandDreare the diffusion coefficients of the oxidizer and the reducer,coxandcreare the oxidizer concentration and the reducer concentration,Riis the amount of consumed reactant corresponding to the generated current at this moment.

Eq.(15)presents the reaction part of this system,reference concentration of the oxidizercox,refand the reducercre,refare both set as 1 mol.m-3.E0is the equilibrium potential of the electrode,Eis the applied voltage on the working electrode,the scanning speed is 1 V.s-1,the voltage range is [-0.4 V,0.4 V],and the scanning time from -0.4 V to 0.4 V is 0.8 s.

To correspond this cyclic voltammetry model,a reactiondiffusion CA model is established,the model description is shown in Eqs.(5)—(8),the grid settings are the same as the diffusion section,IandJare both set as 1000,cell side length is set as 0.05 μm and κtimeis set as 5 × 10-5s,the rotated probabilitypcwandpccware both set as 0.5,which is corresponding to a diffusion coefficient of 5×10-11m2.s-1.There are four possible states for each cell,i.e.,0 for the solvent,1 for the electrode,2 for the reducer,and 3 for the oxidizer.The state of cells with the coordinate ofj=1 is stayed at 1,i.e.,S(i,1)=1,to represent a planar electrode.20 percent of the cells with coordinatej＞1 are set asS(i,j)=2,the other cells are set asS(i,j)=0,the proportion of cell atS(i,j)=2 in the rightmost column is stayed at 20% to match the boundary condition of Eq.(14).The concentration scale factor κconis set as 1 mol.m-3,which means the cellular ratio is numerically equal to the concentration.The reaction probability of a reactant cell located at the electrode surface is defined by Eq.(7),the over-potential η is calculated byE-E0in this system,kcais a fitting parameter to correspond the ability of an electrode to gain or lose electrons,which is set as 0.54 in presented simulation,a complete list of parameter values is shown in the Supplementary Material.

As shown in Fig.10,oxidation probability rises up exponentially and reduction probability presents a symmetrically decreased trend as the potential scanning from -0.4 V to 0.4 V.It should be noted that only the reducer/oxidizer cells next to the electrode cells are capable to be oxidized/reduced.In the CA model,the generated current is calculated by multiplying the cell number reacting at this moment with the mole number filled in one cell,rather than solving the PDEs coupled with concentration and current.The curve peak position is adjusted bykca,and the peak value is adjusted bykmolein Eq.(8),the fitting result is shown in Fig.11,to achieve this agreement,kmoleis set as 9.34 × 10-11mol,which means there are 9.34 × 10-11mol oxidizer or reducer in corresponding cells.

Fig.10.Oxidation and reduction probability curves.

Fig.11.Cyclic voltammetry curves of CA and PDE model.

As shown in Fig.12,the concentration of species evolves in time and space during the process,the reducer substances near the electrode(x=0 μm)are preferentially oxidized,and a sloping concentration distribution along thex-direction is presented for both reduction and oxidation.The consistent concentration distribution of CA model and PDE model also prove that the diffusion part is appropriate to be integrated with the reaction part in this CA framework,because the parameters related to diffusion including κcon,κtime,pcwetc.are set the same as the previous section.

Fig.12.The concentration of reduction and oxidation of the CA model and the PDE model.

3.4.Battery discharge simulation

In this part,a complete battery model including cathode,separator and anode is established based on the CA reconstruction model and the CA reaction-diffusion model.The P2D model is adopted for comparison and verification.The schematic diagram of the P2D model is presented in Fig.13,spherical electrode particles are evenly distributed in the electrolyte,and a separator isolates the direct electron path between the positive and negative electrodes,but allows lithium-ions to diffuse through.Thexdirection is used to calculate the concentration distribution in the electrolyte phase and the potential distribution in both the electrolyte phase and the electrode phase,ther-direction is used to calculate the concentration distribution in the electrode phase.The governing equations and parameter values used for simulation are listed in the Supplementary Material.

Fig.13.Schematic diagram of the P2D model.

To establish a complete battery modelviathe CA framework,a 5000 × 1250 grid is adopted,cell length is set as 0.1 μm to correspond a 50 μm negative electrode,a 25 μm separator and a 50 μm positive electrode.The negative electrode with porosity of 0.75 and average diameter of 12.5 μm is reconstructed by the CA reconstruction method in section 1.1.For the positive electrode,porosity is 0.75 and average diameter is 8 μm.These parameter values are set the same in the P2D and CA simulations.The concentration scale factor κconis set as 4000 mol.m-3,the initial ratio of lithium-ion cells to electrolyte cells in each column is set as 1:1,which is corresponding to a uniform initial concentration of 2000 mol.m-3.Only the discharge process is considered,thus the negative electrode particles are regarded as a source to release lithium-ions and the positive electrode particles are regarded as a sink to absorb lithium-ions,three states are considered for the cells,i.e.,0 for electrolyte,1 for electrode and 2 for lithium-ion.Fig.14 presents the initial situation of this CA model,lithium-ion cells are represented by yellow grids,electrode cells are represented by light blue grids and electrolyte cells are represented by deep blue grids.

Fig.14.The CA model for battery discharge simulation.

The release rate of lithium ions from the negative electrode is calculated by parabolic simplification [29] of the solid diffusion equation,as shown in Eq.(17).

wherecsrepresents the solid concentration,JLirepresents the flux of lithium ions,i.e.,the number of moles of lithium-ions pass through the unit solid phase area,idisis the constant discharge current,Lis the length of the negative electrode,asis the interfacial surface area[30],Rsis the electrode particle radius and εsis the solid volume fraction,all of these parameters values are from the P2D model.Then the increase of liquid concentration in negative electrode is calculated by the decrease of solid concentration as shown in Eq.(17),

wherecerepresents the liquid concentration,εeis the liquid volume fraction,NLi-cellrepresent the number of lithium-ion cells in liquid phase.At each time step,the increase number of lithium-ion cells will be calculated by Eq.(17) in each column,and the increased lithium-ion cells will be randomly filled in the liquid phase cells.

The absorption of lithium-ions in the positive electrode is simulated through the reaction CA model,because there is only a reduction reaction in the positive electrode during discharge process,the BV equation is used to define the reaction probability at the electrode surface andkcais set as 0.125 in Eq.(18).The spatial distribution of over-potential is assumed as a known variable in CA simulation,which has been calculated by the P2D model.At each time step,the lithium-ion cells next to the electrode cell will be determined whether the absorption occurs according to Eq.(18),and the absorbed lithium-ion cells will be replaced by the electrolyte cells.

The diffusion of lithium-ions in electrolyte is simulated through the diffusion CA model,time scale factor κtimeis set as 2 × 10-4s,the rotated probabilitypcwandpccware both set as 0.5,which is corresponding to a diffusion coefficient of 5 × 10-11m2.s-1.

The constant-current discharge simulation results of CA model and P2D model are shown in Fig.15 and Fig.16,the time-space distributions of solid and liquid concentration simulated by these two models are in good agreement during the discharge stage.The negative solid concentration presents a decreasing trend along time and the negative liquid concentration presents an increasing trend,and there is an opposite trend in positive electrode.The solid concentration spatial distribution presents a platform shape,and a slope shape is formed in liquid concentration spatial distribution.The CA simulation results are less smooth then that in the P2D model,one main reason is that the reconstruction electrode particles in Fig.14 are inhomogeneous and irregular,while the spherical electrode particles are evenly distributed in the P2D model.It has been proved that the integration the two CA model to simulate the reaction-diffusion process in porous electrode is feasible,and the effect of inhomogeneous distribution and irregular shape of electrode particles can be investigated by this CA framework.

Fig.15.The solid concentration of lithium-ion of the CA model and the P2D model.

Fig.16.The liquid concentration of lithium-ion of the CA model and the P2D model.

3.5.Discharge efficiency evaluation of porous electrodes

In this part,a set of positive electrodes with different microstructures are constructed and evaluated for the discharge efficiency under the same negative electrode conditions.Porosity and mean particle diameter are the two parameters forced in this work.A set of positive electrodes with same mean particle diameter and different porosity are constructed,another set of positive electrodes with same porosity and different mean particle diameter are constructed,and the schematic diagram of these two sets of electrodes are presented in Fig.17.

Fig.17.Electrodes with different microstructure.

A reaction-diffusion simulation is performed in these electrodes.The grid settings are the same as in the diffusion section,IandJare both set as 1000,cell side length is set as 0.05 μm and κtimeis set as 5 × 10-5s,the rotated probabilitypcwandpccware both set as 0.5,which is corresponding to a diffusion coefficient of 5 × 10-11m2.s-1.There are three possible states for each cell,i.e.,0 for solvent,1 for electrode,and 2 for lithium-ion.The ratio of lithium ion cells to electrolyte cells in each column is 1:9 and κconis set as 1 mol.m-3,corresponding to an initial concentration of 0.1 mol.m-3.The proportion of lithium-ion cells is stayed at 0.5 at coordinate ofj=1 to correspond a fixed negative electrode condition.The state of cells with the coordinate ofj=1000 is stayed at 0,corresponding to a current collector at the right end of the electrode.The reaction possibility of lithium-ion cells being absorbed on the electrode surface is simplified as a constant,i.e.,0.001,rather than being calculated by the BV equation.

Electrode with porosity ranging from 0.6289 to 0.9795 with the mean particle diameter within the range of [1.5400 μm,1.5450 μm] are applied for this reaction-diffusion simulation,the cumulative reaction cell number in 10 s of these electrodes are presented in Fig.18,and the specific data are presented in Table 5.It can be seen that there is a non-monotonic relationship between porosity and cumulative reaction number,as the porosity decreased from 0.9795 to 0.8058,the cumulative reaction number presents a rapid increase trend,as the porosity decreased from 0.8058 to 0.6289,the growth rate of cumulative reaction number slowed down and then the reaction number presents a decreasing trend.To investigate this phenomenon,the spatial distributions of reaction numbers along thex-direction of the electrodes marked by blue,red and yellow in Fig.18 are depicted in Fig.19.

Table 5 The impact of electrode porosity on cumulative reaction number

Fig.18.Cumulative reaction number of different porosity.

Fig.19.Spatial distribution of reaction number effected by porosity.

As shown in Fig.19,the cumulative reaction number near the left boundary is increased with the decrease of porosity,but the reaction number away from the left boundary is decreased with the decrease of porosity.The left boundary is a source to generate lithium-ions and the right boundary is a sink to take away lithiumions.A lower porosity is corresponding to more electrode particles and higher reaction surface,however,diffusion will be hindered when the solid phase volume is too high.The diffusion-hindered effect will be higher than the surface increase effect on reaction efficiency when porosity is lower than the inflection point in Fig.18.It should be noted that the inflection point only exists in the diffusion-controlled systems,and most of the lithium-ions battery systems are diffusion-controlled [31],diffusion rate of lithium-ions is strongly related to the charging/discharging rate of a lithium-ions battery.

For another set of electrodes,porosity is set within the range of[0.7500,0.7600] and the mean particle diameter is ranged from 0.9615 to 5.3615,the number of cumulative reaction cells in 10 seconds is recorded.A monotone relationship between mean particle diameter and cumulative reaction cell number has been found and depicted in Fig.20,the specific data is presented in Table 6.At the same porosity,cumulative reaction number presents an increasing trend as the decrease of mean diameter,which is caused by the increased reaction surface area.The spatial distributions of reaction numbers of the electrodes marked by blue,red and yellow in Fig.20 are depicted in Fig.21.It can be seen that the diffusionhindered effect is nearly for these three electrodes,as the effective diffusion distance is all around 25 μm.At this porosity,the discharge efficiency is monotonous related to the mean particle diameter.The cumulative reaction cell number can be transformed to reaction mole number or reaction concentration by multiplying this result by scale factors,then this framework can be applied to seek the microstructure with better charging/discharging efficiency of lithium-ions batteries.

Table 6 The impact of electrode particle diameter on cumulative reaction number

Fig.20.Cumulative reaction number of different diameter.

Fig.21.Spatial distribution of reaction number effected by diameter.

4.Conclusions

In this paper,we have introduced a new framework for electrode reconstruction and reaction-diffusion simulation based on CA.In this framework,the computation way is fully synchronous that allows for high-efficient parallel implementations;irregular solid-liquid boundary conditions are easy to handle in this framework;a physical connection between CA parameters and actual variables is established.Through this framework,electrodes with specific microstructure are constructed;quantitative agreement between CA simulations and theoretical solutions are obtained for several reaction-diffusion problems,including Fick diffusion,cyclic voltammetry test and battery discharging process;the discharge efficiencies of electrodes with different morphology are evaluated and analyzed.This CA-based framework provides a simple and effective alternative to investigate the influence of microporous structure on battery performance.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors gratefully acknowledge the National Natural Science Foundation of China (21878012).

Supplementary Material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.cjche.2023.01.022.