Quasi-3D modeling of Li-ion batteries based on single 2D image

Electrochemical physics-based simulations of Li-ion batteries using a mesoscale 3D structure of porous electrodes are one of the most effective approaches for evaluating the local Li concentration in active materials and the Li-ion concentration in electrolytes. However, this approach requires considerable computational resources compared with a simple 2D or 1D homogeneous simulation. In this work, we developed an advanced electrochemical physics-based simulation method for Li-ion batteries that enabled a quasi-3D simulation of charge/discharge using only a single 2D slice image. The governing equations were based on typical theories of electrochemical reactions and ion transport. From referencing the 2D plane, the model was able to simulate both the Li concentration in the active material and the Li-ion concentration in the electrolyte for their subsequent consideration in a virtual 3D structure. To confirm the validity of our proposed model, a full 3D discharge simulation with randomly packed active material particles was performed and compared with the results of the quasi-3D model and a simple-2D model. Results indicated that the quasi-3D model properly reproduced the sliced Li and Li-ion concentrations simulated by the full 3D model in the charge/discharge process, whereas the simple-2D simulation partially overestimated or underestimated these concentrations. In addition, the quasi-3D model made it possible to dramatically decrease the computation time compared to the full-3D model. Finally, we applied the model to an actual scanning electron microscopy equipped with a focused ion beam (FIB-SEM) image of a positive electrode.


Introduction
Lithium-ion (Li-ion) batteries are recognized as the most promising technology for energy storage because of their high energy density, lightweight and long cycling life [1,2]. In addition to various studies on novel materials for electrodes/electrolytes and new battery systems, a variety of simulation technologies [2,3] have been proposed to predict charge/discharge performances, stress conditions and cycling life times. Electrochemical physics-based models (physicochemical models)  are useful tools to calculate the various nonlinear resistance components of a battery, including the diffusion of Li, stress in active material particles, electrochemical reactions and ionic transport in electrolytes, while simple equivalent circuit models [26][27][28][29][30] assume that the resistance of the battery is constant or a function of the current and temperature.
Until now, one-dimensional (1D) electrochemical physics-based models, which assume porous electrodes and separators as uniform media, have been widely adopted and developed. These models have been applied to analyze various electrode materials, including Li x Mn 2 O 4 , LiCoO 2 , LiFePO 4 , and Li(NiCoMn)O 2 , in cathodes and graphite, and Li metal in anodes, along with heat generation, cycle degradation and stress analyses [4][5][6][7][8][9][10]26]. In addition, these models are utilized for systematic studies with a large number of parameter datasets because their calculation costs are relatively low. In the 1D models, however, it is difficult to evaluate nonuniform Li and Li-ion concentrations and reactions, especially in the cross section of electrodes.
Recently, many studies of charge/discharge simulations based on the mesoscale three-dimensional (3D) structure of porous electrodes have been reported [11][12][13][14][15][16][17][18][19][20][21][22]. In these cases, porous electrodes were modeled by random packed spheres/hemispheres [11][12][13][14][15] or actual structures based on scanning electron microscopy equipped with a focused ion beam (FIB-SEM) results [17][18][19][20][21][22] to evaluate the 3D distribution of Li in the active material particles, Li-ion concentration in the electrolyte, and the stress distribution and temperature field of electrodes. These models enable charge/discharge simulation with realistic 3D electrode structures and provide insightful information at the mesoscale; however, there are some limitations in these models: (1) they require a large-scale amount of calculation and (2) they require 3D tomography data, including 3D-SEM/X-CT results or particle-packed artificial structures, which also require significant cost. Thus, simple two-dimensional (2D) simulations using the crosssectional structure of a battery electrode have been performed [23][24][25]. Although these approaches decrease the computational time compared to 3D simulations, they still have some validation problems because of ignoring the effect of the Z-direction.
In this study, we propose a new electrochemical physics-based simulation method for Li-ion batteries that enables a quasi-3D calculation of charge/discharge and a dramatic decrease in the amount of calculation by using a single 2D slice image of porous electrodes for consideration in a virtual 3D structure. In the quasi-3D model, the Li concentration in the active materials and the Li-ion concentration in the electrolyte are simulated in the 2D plane, for consideration in the virtual 3D structure using a single particle model. The study involved the following steps. First, we discuss the inference accuracy of the radius of the active material particle from the 2D plane. Next, the validity and advantage of the calculation cost of this model are confirmed by comparing the results of the 3D discharge simulation with random-packed active material particles. Moreover, we applied the model to an actual FIB-SEM image of a positive electrode and evaluated the distributions of Li and Li-ion concentrations in the plane.

Model development
The quasi-3D model is based on electrochemical reaction and ion transport theories, which are widely used in 3D, 2D and 1D models . In this model, the Li concentration in the active materials and the Li-ion concentration in the electrolyte are simulated from a 2D slice image of a mesoscale 3D electrode structure in a half-cell. The key idea of the model is that the Li and Li-ion concentrations are corrected compared to those in the single particle model that is performed in parallel. Here, we introduce the geometric and numerical schemes used in this approach.

Geometry
In this study, we use a half-cell structure of a positive electrode for galvanostatic discharge simulations. The model geometries for full-3D, simple-2D and quasi-3D simulations are constructed as follows. First, the mesoscale structure of the 3D porous electrode is constructed by randomly packing active material particles. Perfect spheres are assumed to be active material particles and overlaps between these particles are permitted. The radii of the spheres R 3D of the spheres are randomly packed sequentially with a uniform distribution in the L x x L y x L z simulation cell (the x-and y-axes are the in-plane directions, and the z-axis is the thickness direction of the electrode) until the volume ratio of the active material reaches AM .
Here, binder and additives are not modeled for the sake of simplicity. Table 1 shows the parameters adopted for generating the 3D porous structure, which are typical values for positive electrodes. This constructed structure is used for a full-3D charge/discharge simulation. Subsequently, slice planes for the simple-2D simulation and quasi-3D simulation are extracted from the middle position of L y . Figure 1 illustrates the 3D porous electrode structure and the extracted plane.

Governing equations of the quasi-3D model
The mathematical model is based on the electrochemistry and Li transport model on the 2D plane for consideration in the virtual 3D structure. The governing equations in the quasi-3D model are listed in Table 2. The Butler-Volmer equation is assumed on the interface between the active particle and electrolyte [4]. The local current density on the interface i loc is described as where F, R, α and i 0 indicate the Faraday constant, gas constant, transfer coefficient and exchange current density, respectively. Additionally, η denotes the overpotential between the active material and electrolyte where φ s , φ l and U are the active material potential, electrolyte potential and open-circuit potential (OCP), respectively. For simplicity, φ s is assumed to be uniform in the active material used in this study [19]. A typical OCP function U for the positive electrode Li(NiCoMn)O 2 is adopted. The polynomial function is used in order to connect the OCP data with a smooth curve, written as where c s,max is the coefficient of the maximum Li concentration of the active material and p i indicates the coefficient of the polynomial OCV function provided in Table 3. In most cases, Fick's law is used in order to model the Li diffusion in the active material particles of an electrode [23,24]. Herein, we propose the modified Fick's law to take into account the effect of the structure in the y-direction. Thus,  where c s,2D and D s are the Li concentration in the 2D plane and the Li diffusion coefficient of the active material, respectively. Moreover, q s,corr is used to correct the 2D Li concentration and is described as the difference between the fluxes of Li in the y-axis direction,

Li-ion concentration in the electrolyte
where q surf s,corr and q btm s,corr are the molar fluxes of the spherical y-axis surface and bottom, respectively. In this study, we assume that the fluxes are proportional to the concentration gradient, which is evaluated by the reference Li concentrations on the spherical surface c ref s,surf and bottom c ref s,btm , and the distances between these positions in the 2D plane (see Fig. 2(a)).
where the d parameters are determined by the geometric structure of the electrode. In this study, the average length (5) q s,corr= q surf s,corr + q btm s,corr are estimated by a single-particle model that is conducted in parallel.
In previous studies [23,24], the Li-ion concentration c l,2D and potential l,2D in the electrolyte are estimated by the mass conservation law and Nernst-Plank equation. Herein, we add the correction term q l to the mass conservation law in consideration of the effect of the flux in the z-direction, where D l , l and t + are the diffusion coefficient, ionic conductivity and transport number in the electrolyte, respectively. At the boundary between the electrode and separator (z = 0 µm), the Dirichlet condition c ref l = c l,0 is applied. The 2 nd term of the right-hand side of Eq. (7) q l,corr is a where d l is the parameter for the distance to the reference concentration. Herein, we assume d l to be 10 µm, which is the representative length in the electrolyte. Figure 2(b) shows a schematic image of the correction of the 2D Li-ion concentration.
As mentioned above, the reference Li and Li-ion concentrations are estimated by the single-particle model [31]  where j ave means the average molar flux associated with the electrochemical reaction between the active material and electrolyte and i app indicates the applied current density. On the other hand, the reference Li-ion concentration in the electrolyte is evaluated by Thus, we ignore the relaxation time of the Li-ion concentration in the electrolyte in this description. The parameters for the typical positive electrode Li(NiCoMn)O 2 and electrolyte 1 M LiPF 6 /EC:DEC = 1:1 by volume are adopted and listed in Table 4.

Estimation of R 3D from the 2D plane
In the quasi-3D model, the actual radius of the 3D active material R 3D is required to calculate the reference Li concentration c ref s . However, R 3D is not equal to R 2D because the "active material disc" on the 2D plane is a result of a randomly sliced "active material sphere." Herein, we use a Bayesian inference to determine R 3D from the 2D plane.
where P( |x) , P( ) and P(x| ) indicate the posterior distribution, prior distribution and likelihood, respectively. The probability of the radius of the active material on the 2D plane (disk shape), R 2D , from the actual radius R 3D in 3D (spherical shape) is described as The derivation of Eq. (14) is shown in Appendix A.  Assuming that we have no prior information about R 3D , the prior can be described as according to the principle of insufficient reason. Therefore, the posterior of R 2D on the 2D plane from the R 3D of the active material sphere can be written as where R 2D,i indicates the ith radius of the active material disk on the 2D plane.

Governing equations of the full-3D and simple-2D models
The governing equations in the full-3D and simple-2D models are based on electrochemistry and Li transport (15) P R 3D = const.
theories without the correction introduced in the quasi-3D mode. The mass conservation equation of the Li concentration is applied to the active material, and the Nernst-Plank equation and mass conservation equation are adopted to evaluate the potential and Li-ion concentration in the electrolyte. The Butler-Volmer equation is used on the interface between the active particle and electrolyte. In addition to the quasi-3D model, the potential in the active material is assumed to be uniform for simplicity. The governing equations in the full-3D and simple-2D models are listed in Table 5.
The modeling and calculation of the three models were performed using COMSOL Multiphysics™ ver. 5.4 with a standard workstation containing a 16-core Intel Xeon™ (2.60 GHz) processor and 128 GB of RAM.

Results and discussion
The proposed quasi-3D model was validated with a full-3D model and compared to a simple-2D model under various applied current densities. First, the 3D active material radius R 3D was estimated from the radius of the active material disk R 2D by a Bayesian inference. Additionally, we carried out galvanostatic discharge simulations at a low current density and evaluated their results, including the Li concentration in the active material, Li-ion concentration in the electrolyte and voltage profile. Finally, simulations at a high current density were conducted, and the model limitations and robustness were discussed.

Estimation of R 3D
The center positions and radius of each active material disk in the 2D slice structure, shown in Fig. 1(b), were detected by the watershed method. The histograms of the detected  Fig. 3(a). Most radii are near 10-11 µm. Using these radii, we estimated the 3D active material radius R 3D by a Bayesian inference as described in Eq. (16). Figure 3(b) shows the posterior distribution of R 3D . Note that the vertical axis is described as an arbitrary unit. The highest probability is located at approximately 10.5 µm, which is close to the actual radius of 11.0 µm. This indicates that the actual 3D radius of the active material sphere can be inferred from the 2D radii of the disks on an arbitrarily sliced image. Although a large number of disks improves the estimation accuracy, the slice image that includes only 16 disks used in this study is sufficient to infer the actual 3D radius. The difference in accuracy of the R 3D inference using the slice position is discussed in Appendix.

Low current density (0.5C)
Galvanostatic discharge simulations at 0.5C (4.81 [A/m 2 ]) were conducted using the full-3D, quasi-3D and simple-2D models. The cell voltage profiles in these models are described in Fig. 4. The discharge voltages within 7000 s in the 2D-simple and quasi-3D models are close to those in the full-3D model. However, the simple-2D model overestimates the capacity (the voltage drops at approximately 9000 s), whereas the quasi-3D model estimates it more accurately than the full-3D model (the voltage drops at approximately 7200 s). Figure 5 shows the Li concentration in the active material at 50% SOC. In the 3D view of the full-3D model ( Fig. 5(a)), one can see the concentration gradient from the surface to the center of the sphere. We compared the slice image of the full-3D concentration with that of the quasi-3D and simple-2D concentrations, as shown in Fig. 5(b). The simple-2D model underestimates the Li concentration inside the large region (point A) and overestimates it in the small region (point B) because it simulates only the 2D plane (x-z plane), ignoring the flux in the y-direction. On the other hand, the Li concentration in the active material of the quasi-3D model agrees well with the full-3D mode because of the correction with the average concentration of the single-particle model. Figure 6(a) illustrates the concentration profiles of these models on the X 1 -X 1 ′ lines. It is apparent that the profile of the quasi-3D model is close to that of the full-3D model, whereas the predicted concentration of the simple-2D model is relatively poor. The time sequence of Li concentrations at points A and B in these models is described in Fig. 6(b) and (c), respectively. Compared with the poor agreement of the simple-2D model at both points, the Li concentration and its increase in the quasi-3D model agree with those in the full-3D model. The reasons that the Li concentration in the simple-2D model differs from the full and quasi-3D models are as follows. First, the specific surface area of the active material in the 2D plane is different from those of 3D structure. This could affect the Li concentration near the surface in the active material because the total Li flux from the surface is under/overestimated. Second, the simple-2D model ignores the Li flux of z-direction. This could increase the difference of the Li concentration between them. In the quasi-3D model, although the specific surface area of the active material has the same value of the simple-2D model, the correction terms in Eq. (6) and (9) prevent the increase of the difference from those of the full-3D model. The Li-ion concentration in the electrolyte at 50% SOC is illustrated in Fig. 7. As seen in Fig. 7(a), the concentration gradient is formed from the boundary of the separator (bottom plane) to the boundary of the current collector (top plane). The slice image of the full-3D concentration is compared with that of the quasi-3D and simple-2D concentrations, as shown in Fig. 7(b). The simple-2D model underestimates the Li-ion concentration over the whole electrolyte region because it ignores the flux in the y-direction. The Li-ion concentrations on the X 2 -X 2 ′ lines are illustrated in Fig. 8(a). Compared to that of the simple-2D model, the concentration profile of the quasi-3D model is relatively close to that of the full-3D model. In Fig. 8(b), the variation of the Li-ion concentrations over time at point C of these models is described. Additionally, compared to those in the simple-2D model at both points, the Li-ion concentration and its increase in the quasi-3D model are relatively close to those in the full-3D model at both points. As in the case of the Li concentration in the active material, the quasi-3D model has the followability of the Li-ion concentration to the full-3D model by taking the correction terms of the Li-ion flux in Z-direction.

High current density (1C and 2C)
Discharge simulations of the full-3D, simple 2D and quasi-3D models were carried out at high current densities of 1 and 2C. A comparison of the discharge curves of the three models is described in Fig. 9. At 1C, the voltages within 2800 s in the simple-2D and quasi-3D models are close to that in the full-3D model. The simple-2D model apparently overestimates the capacity (the voltage drops at approximately 3700 s), whereas the quasi-3D model estimates it more accurately and matches more closely to the full-3D model (the voltage drops at approximately 3100 s). On the other hand, the discharge voltage of the simple-2D model at 2C drops rapidly at the beginning of the simulation due to Li-ion exhaustion near the active material in the electrolyte. Notably, the voltages of the quasi-3D and full-3D models are similar to each other at 0.5 and 1C.
The Li concentrations in the active material at 50% SOC of the models are shown in Fig. 10. Note that the data of 2C in the simple-2D model are not available because of a calculation failure. As evidenced in the 0.5C case, the simple-2D model underestimates the Li concentration inside the large region (point A) and overestimates it in the small region (point B), whereas those in the active material of the quasi-3D model agree well with the full-3D model. To evaluate the prediction accuracy of the models, the deviation between the simulated values of the full-3D model and those of the quasi-3D/simple-2D models is defined as follows: where x i,F3D and x i,m indicate the simulated value of the full-3D model and the quasi-3D/simple-2D models (m = Q3D, S2D). Figure 11(a) shows the deviation of the Li concentration in the active material at 0.5, 1.0 and 2.0C. Not that there are no data of the simple-2D model at the 2C discharge condition due to the rapid voltage drop at the beginning of the simulation. It is apparent that the deviations of both models increase with the C-rate because the large molar flux at the boundary between the active material and electrolyte leads to the nonuniform concentration in the active material having a larger effect. At 1C, the deviation of the quasi-3D model (595 mol/m 3 ) remains lower than the value of the simple-2D model (1518 mol/m 3 ). The same tendency is seen in the Li-ion concentration in the electrolyte, as described in Fig. 11(b). The deviations of both models increase with the C-rate, whereas at 1C the deviation of the quasi-3D model (155 mol/m 3 ) remains lower than the value of the simple-2D model (295 mol/m 3 ).
These results indicate that the quasi-3D model can be a useful tool to predict the Li and Li-ion concentrations in the 2D plane as long as the applied C-rate is not relatively large. One reason for the impact on the accuracy of the quasi-3D model is that the deviation between the full model and the single particle model used by the quasi-3D model could increase with C-rate. In the single particle model, one representative particle is modeled, and it cannot consider the effect of the concentration gradient in the thickness direction, formed in the case of higher C-rate in particular. Although only discharge simulations have been conducted in this study, charge simulations could exhibit

Effect of slice position
In order to investigate the effect of the slice position on the accuracy of the quasi-3D model, the additional studies have been performed using the slice images near both ends. In the same way as Sect. 3.2, the galvanostatic discharge simulations at 0.5C have been conducted by the full-3D, quasi-3D and simple-2D models, respectively. Figure 12 shows the Li concentration in the active material at 50% SOC in each model. Apparently, the simple-2D model overestimates the Li concentration in the small region (shown by arrows) in both slice positions (slice 1 and 2), whereas those of the quasi-3D model agrees well with the full-3D model. This result indicates that the quasi-3D model has higher followability to the full-3D model than the simple-2D model even in other slice positions. Table 6 shows a comparison of the calculation times of the full-3D and quasi-3D models at 0.5C. The calculation time of the quasi-3D model is approximately one hundred times that of the full-3D model. Thus, the quasi-3D model using a single slice image makes it possible to dramatically decrease the calculation cost while still providing accurate cell voltages and Li and Li-ion concentration distributions.

Application to an actual FIB-SEM image
The proposed quasi-3D modeling is expected to be applied to an actual electrode structure. To test the model with an actual electrode structure, a discharge simulation using an FIB-SEM image of a positive electrode Li(NiCoMn) 1/3 O 2 was conducted. In addition to the scheme for the particle packing structure described previously, the 3D active material radius R 3D was estimated from the active material radius R 2D obtained from an SEM image by a Bayesian inference. Subsequently, we carried out galvanostatic discharge simulations at 1C with typical parameters.

FIB-SEM image segmentation and inference of the active material size
Herein, the FIB-SEM image positive electrode Li(NiCoMn) 1/3 O 2 , which was taken in a previous study [31], was used, and image segmentation was conducted to produce two regions (active material and electrolyte). Figure 13 shows the FIB-SEM image and segmentalized image. Subsequently, we detected the radius of each active material on the 2D slice structure. As shown in Fig. 13(a), the active material particles are not spherical but polyhedral, and R 2D is evaluated using their circumscribed circles. A histogram of the radii of 111 active material polygons is shown in Fig. 14(a). Using these radii, we estimate the 3D active material radius R 3D by a Bayesian Fig. 9 Comparison of the discharge curves of the full-3D, quasi-3D, and simple-2D models at a 1.0C and b 2.0C

Galvanostatic simulation based on the FIB-SEM image
Galvanostatic discharge simulations at 2C (19.2 [A/m 2 ]) were conducted using the segmented FIB-SEM image. Although the parameters for the simulation should be determined in each system, we adopted the same parameters used in the spherical packing structure (Table 4) to test the application possibility of the quasi-3D model to

Conclusion
In this study, we proposed a new electrochemical physicsbased simulation method for Li-ion batteries that enables a quasi-3D calculation of charge/discharge and a dramatic decrease in the amount of calculation by using a single 2D slice image of porous electrodes for consideration in a virtual 3D structure. First, the actual 3D radius of the active material sphere was inferred from the 2D radii of disks on an arbitrarily sliced image. Although a larger number of disks improves the estimation accuracy, we concluded that the slice image including only 16 disks used in this study was Fig.12 Li concentration in the active material at 50% SOC in each model using the slice images near both edges

Appendix B
To confirm the difference in the accuracy of the R 3D inference due to the slice position of the actual 3D structure, we compare the estimated R 3D from the slice positions y = 1/3L y , 2/3L y and 3/3L y (see Fig. 16(a)) in addition to the case of y = 1/2L y described in Sect. 3.1. The number of active material disks on the plane is 15 in the case of y = 1/3L y , 9 in the case of y = 2/3L y , and 16 in the case of y = 3/3L y . Figure 16(b) shows the histogram of detected radii in these three cases. Most radii are near 10.5-12.0 µm in all cases. Using these radii, we estimate the 3D active material radius R 3D by a Bayesian inference as described in Eq. (16). Figure 16(c) describes the posterior distribution of R 3D . Although there are some fluctuations between them, the highest probability of R 3D is located at approximately 10-11 µm. Thus, we can conclude that the accuracy of the R 3D inference at various slice positions is almost the same in this structure. Error between the inferred and the actual R 3D is estimated in each slice position, as shown in Fig. 17.
The errors are within 10% in all slice positions, and they are in the acceptance range. However, note that the accuracy may change based on the system size, number of particles and particle size. Namely, the larger the system size is, the larger the number of particles; thus, the inference accuracy increases.