1 Introduction

Lithium-ion batteries, since first commercialized by Sony in 1990s, have been widely used as power sources for mobile and stationary consumer electronics, and are perceived to be promising alternatives to power sources of electric vehicles and other medium- to large-sized power or energy storage devices, owing to relatively high energy and power densities, almost no memory effect, and low self-discharge, etc. However, some fundamental mechanisms of lithium-ion battery processes still remain unknown or have not yet reached a widely-accepted cognition [1, 2]. Safety issues, poor cycle life, and high cost etc. continually remind us that the related technologies are still far away from a mature state [1–4].

While experiments are necessary to study the underlying physical chemistry of materials used in lithium-ion batteries and evaluate the battery performance, computational models are useful in providing a fundamental understanding of internal transport mechanisms, and are playing an increasingly important role in research, design and optimization of lithium-ion batteries. In the past two decades, various numerical models have been developed for lithium-ion batteries. Early isothermal electrochemical models were reported by Atlung et al. [5, 6] and West et al. [7]. Doyle et al. [8], Fuller et al. [9], Darling et al. [10], and Doyle et al. [11] developed an isothermal electrochemical model for a dual-insertion lithium-ion cell and found that the model predictions were in good agreement with the experimental data. By this model or adapted variations, battery behaviors and cell/electrode design were extensively explored by Doyle and Fuentes [12], and Nagarajan et al. [13]. Arora et al. [14, 15] even used this class of models to study diffusion limit, side reactions, and capacity fade mechanisms in lithium-ion batteries. Thermal behaviors of lithium-ion batteries affects battery cycle life and safety operation. Pals and Newman [16, 17] developed a one-dimensional thermal model for a single-cell lithium/polymer battery and later extended for a cell stack. Using this thermal model, Botte et al. [18] studied the influence of some design parameters on the thermal behavior of lithium-ion batteries. Chen and Evans [19] developed a two-dimensional model, later extended to three-dimensional [20], to study the effects of various cell components, stack size and cooling conditions on the performance of Li-polymer electrolyte batteries under different discharge rates. Gu and Wang [21] and Srinivasan and Wang [22] developed a thermal-electrochemical coupled model for lithium-ion cells. Fang et al. [23] used this model to predict the performance of an automotive Graphite/NMC lithium-ion cell, including overall cell performance and individual electrode performance, at various operating temperatures.

The above-summarized models can be generally classified as macroscopic models. The macroscopic electrochemical models solve simultaneously the charge (Li+ in electrolyte phase and electron in solid phases) and species (Li+ in electrolyte phase and Li in solid active materials) transport equations. The macroscopic thermal- electrochemical models are developed based on the electrochemical model coupled with an energy conservation equation, and has been validated to be useful in the research and development of lithium-ion batteries. The macroscopic models provide some insight into the limiting mechanisms and predict behaviors of batteries under a wide range of operating conditions, thus a powerful tool for design and optimization of lithium-ion battery energy systems. However, the macroscopic models generally assume the electrodes and separator are homogenous porous media with all the detailed morphology about microscopic pore-structures neglected. That is to say, everywhere inside the electrodes and separator the electrolyte and solid matrix coexist and the only characteristic parameter used to delineate the porous configuration is the component fraction (or porosity). Physical properties are approximated in terms of the component fraction based on a concept of equivalence. This macroscopic treatment to the electrodes and separator of lithium-ion batteries may lose validity for real lithium-ion batteries as the solid particles of various shape and size are not uniformly distributed inside the electrodes and separator. For better describing the cell behaviors and engineering the pore-structures of electrodes and separator, it is necessary to resort to pore-scale models [24].

In recent years, a few researchers attempted to do pore-scale modeling to the multi-disciplinary processes coupled with electrochemical reactions in lithium-ion batteries. Zhang et al. [25, 26] considered a single particle and simulated intercalation-induced stresses and heat generation induced stresses inside spherical or ellipsoidal LiMn2O4 particles together with a surrogate-based analysis. Wang and Sastry [27] developed a mesoscopic pore-scale model based on a microstructural system of regularly or randomly packing spherical or ellipsoidal particles, which represent the solid active materials in electrodes. Yi et al. [28] presented a pore-scale model to study the effect of compression of particulate systems in graphitic lithium-ion anodes on electrode properties and performance. Garcia et al. [29], Garcia and Chiang [30], and Smith et al. [31] proposed a pore-scale electrochemical model and used it to investigate the electrode microstructural effects and to design electrodes. To gain insight into the collector-electrode interactions, including mechanical deformation and contact resistance, Awarke et al. [32] established a three-dimensional mesoscopic particle-scale model, in which a representative volume element was generated using either a random packing or a dynamic collision algorithm considering different phases inside the electrodes. In summary, mesoscopic pore-scale modeling has the potential to unravel the underlying mesoscopic pore-scale mechanisms of multi-disciplinary transport coupled with electrochemical reactions and thus can provide insightful suggestions or even ultimate solutions to the issues that retard further development and application of lithium-ion batteries.

In the aforementioned mesoscopic pore-scale models, the finite element method, which is a traditional Eulerian-based numerical method, has been generally used to solve the model equations. A relatively novel numerical technique—smoothed particle hydrodynamics (SPH) method will be adapted for pore-scale modeling of lithium-ion battery processes in the present work. SPH method is a meshless particle-based Lagrangian fluid dynamics simulation technique and well-suited to pore-scale simulations [33]. In SPH, the simulation domain is represented by a collection of discrete elements or pseudo-particles. These elements or particles are initially distributed with a specified density distribution and evolve in time according to the conservation equations (mass, momentum, energy, etc.). Flow and energy transport properties at a certain position are determined by an interpolation or smoothing of the nearby particle distribution using a weighting function (the smoothing kernel). The SPH method and its extensions have been successfully implemented in simulating diverse problems of fluid dynamics and heat transfer [34–40]. It has been demonstrated that computational techniques dealing with multi- disciplinary transport processes involving complex physico-chemistry can be incorporated in the SPH as well. For example, for simulations dealing with intricate interfaces, diffuse interface and reaction-diffusion technique can be embedded in the SPH method [41, 42].

The SPH method has two defining characteristics, meshfree and Lagrangian-based, which make it more advantageous than the traditional Eulerian-based numerical techniques in what concerns the following aspects: (1) particular suitability to tackle multi-disciplinary transport processes, (2) ease of handling complex free surface and material interface, (3) relatively simple computer codes and ease of machine parallelization, and (4) intrinsic transient essence. Application of the SPH method to the pore-scale modeling of lithium-ion battery can fully benefit from these advantages.

Mesoscopic pore-scale modeling to multi-disciplinary processes coupled with electrochemical reactions in lithium-ion battery consists of two steps. The first step is the microstructural construction for electrodes and separator comprising of various phases. The second step is the establishment and solution of mesoscopic pore-scale mathematical model. The present work will adapt the SPH method for mesoscopic pore-scale modeling for lithium-ion batteries. The main goal is to demonstrate the reasonability and suitability of SPH.

2 Model development

2.1 Physicochemical description

We take the LiCoO2|LiPF6|LiC6 battery as an example. It consists of an anode current collector (Cu), an anode electrode, a separator, a cathode electrode, and a cathode current collector (Al), schematically shown in Fig. 1. The electrodes and separator are all composite. From focused ion beam-scanning electron microscope (FIB/SEM) images, we can distinguish the components that make up the electrodes and separator. Three components are distinguished in the anode electrode: Li x C6 as active materials, pores full of electrolyte, and additives (CMC + SBR etc.). The three distinguishable phases in the cathode electrode are the active materials Li y−x CoO2, the electrolyte filling the pores, and additives (SP + PVDF etc.). The composite separator consists of polymer matrix and electrolyte. The electrolyte, a good ionic conductor but an electronic insulator, provides a medium for Li+ to travel between the two electrodes and keeps electrons flowing in the solid phases (except polymer matrix) and along the external circuit. Electrochemical reactions occurring at the solid active materials/ electrolyte interfaces are as follows:

Fig. 1
figure 1

Schematic of a lithium-ion battery. On the second row are presented typical FIB/SEM images of electrodes and separator

Full size image
$$ {\text{For the graphite electrode}}{:}\,{\text{Li}}_{x} {\text{C}}_{ 6} \underset{\text{charge}}{\overset{\text{discharge}}{\longleftrightarrow}}{\text{Li}}_{ 0} {\text{C}}_{ 6} { + }x{\text{Li}}^{ + } { + }x{\text{e}}^{ - } ; $$
(1)
$$ {\text{For the Li}}_{y - x} {\text{CoO}}_{ 2} {\text{electrode}}{:}\,{\text{Li}}_{{y{ - }x}} {\text{CoO}}_{ 2} { + }x{\text{Li}}^{ + } { + }x{\text{e}}^{ - } \underset{\text{charge}}{\overset{\text{discharge}}{\longleftrightarrow}}{\text{Li}}_{y} {\text{CoO}}_{ 2} . $$
(2)

During discharge, Li extracts from Li x C6 particles in the anode electrode (Eq. (1)), travels in the electrolyte, passes the separator and inserts into Li y−x CoO2 particles in the cathode electrode (Eq. (2)). Simultaneously, electrons released travel in the solid phases (including additives and active materials), flow through an external circuit into the cathode, and are then consumed by the cathode reaction. During charge the reverse processes occur.

2.2 Mathematical description

2.2.1 Model assumptions

As a first step toward establishing the model equations, we made the following assumptions.

  1. (i)

    Volume change of material during cell operation is not taken into account, i.e., constant porosity in electrodes and separator;

  2. (ii)

    Charge transfer kinetics follows the Butler-Volmer equation;

  3. (iii)

    Electrochemical reactions occur at the solid active materials/electrolyte interfaces;

  4. (iv)

    Interfacial electrical equilibrium exists in both the electrolyte and solid active materials and interfacial chemical equilibrium exists in the electrolyte phase;

  5. (v)

    Side reactions do not occur;

  6. (vi)

    Ionic species transport in the electrolyte is by diffusion and migration, and lithium transport inside the solid active material particles by diffusion;

  7. (vii)

    The solid and electrolyte phases are immobile during battery operations.

2.2.2 Governing equations

The multi-disciplinary transport processes coupled with electrochemical reactions in lithium-ion battery can be described by a series of conservation equations.

(i) Charge conservation equations The electron charge conservation in solid phases (including additives and active materials) is described by

$$ \nabla \cdot \left( {\sigma \nabla \phi_{\text{s}} } \right) - j^{\text{Li}} = 0, $$
(3)

where σ is the electronic conductivity. Boundary conditions for constant current (I) discharge/charge are:

$$ - \sigma \left. {\frac{{\partial \phi_{\text{s}} }}{\partial x}} \right|_{x = 0} = - \sigma \left. {\frac{{\partial \phi_{\text{s}} }}{\partial x}} \right|_{x = X} = \frac{I}{A}, $$
(4)

where X is the thickness of lithium-ion battery, A is the surface area of current collector. The Li+ charge conservation in electrolyte phase is described by

$$ \nabla \cdot \left( {\kappa \nabla \phi_{\text{e}} } \right) + \nabla \cdot \left( {\kappa_{\text{D}} \nabla \,\ln c_{\text{e}} } \right) + j^{\text{Li}} = 0. $$
(5)

The ionic conductivity κ is given by [43]

$$ \kappa = 4.1253 \times 10^{ - 4} + 5.007c_{\text{e}} - 4.7212 \times 10^{3} c_{\text{e}}^{ 2} + 1.5094 \times 10^{6} c_{\text{e}}^{ 3} - 1.6018 \times 10^{8} c_{\text{e}}^{ 4} . $$
(6)

The ionic diffusional conductivity, κ D, is given by

$$ \kappa_{\text{D}} = \frac{2RT{\upkappa} }{F}\left( {t_{ + }^{0} - 1} \right)\left( {1 + \frac{{{\text{d}}\,\ln f_{ \pm } }}{{{\text{d}}\,\ln c_{\text{e}} }}} \right), $$
(7)

where F is the Faraday’s constant, R the universal gas constant, and T the absolute temperature in Kelvin. \( t_{ + }^{ 0} \) denotes the transference number of Li+. Depending on the combination of electrolyte and solvent, the transference number \( t_{ + }^{ 0} \) can be a function of Li+ concentration in electrolyte. The mean molar activity coefficient of the electrolyte f ± is assumed constant in the present work. In the two current collectors, there exist only electronic current (ionic current equals zero). Accordingly, we obtain boundary conditions for electrolyte potential as

$$ \left. {\frac{{\partial \phi_{\text{e}} }}{\partial x}} \right|_{x = 0} = \left. {\frac{{\partial \phi_{\text{e}} }}{\partial x}} \right|_{x = X} = 0. $$
(8)

The transfer current density j Li occurring only at the solid active material/electrolyte interface is determined from the Butler-Volmer equation, which describes the reaction flux on the particle surface as a function of the exchange current density i 0 and the surface overpotential η int, as

$$ j^{\text{Li}} = \frac{{A_{lm} }}{V}i_{0} \left[ {\exp \left( {\frac{{\alpha_{\text{an}} F\eta_{\text{int}} }}{RT}} \right) - \exp \left( { - \frac{{\alpha_{\text{ca}} F\eta_{\text{int}} }}{RT}} \right)} \right], $$
(9)

where l, m denotes solid active material and electrolyte, respectively. \( \alpha_{\text{an}} \) and \( \alpha_{\text{ca}} \) are the anodic and cathodic transfer coefficients of electrode reaction j Li. As A lm is the interfacial area of phase l and m, and V is the total volume of the control unit, the ratio of A lm to V is actually the local specific surface area of solid active materials. Exchange current density, i 0, exhibits modest dependency on the Li+ concentration in electrolyte and Li concentration on the solid active material surface, respectively.

$$ i_{ 0} = kFc_{{{\text{s,}}\hbox{max} }} \left( {c_{\text{e}} } \right)^{{\alpha_{\text{an}} }} \left( {\frac{{c_{\text{s}} }}{{c_{\text{s,max}} }}} \right)^{{\alpha_{\text{ca}} }} \left( {1 - \frac{{c_{\text{s}} }}{{c_{\text{s,max}} }}} \right)^{{\alpha_{\text{an}} }} , $$
(10)

where k is the rate constant for the electrochemical reaction. The overpotential η int is associated with the solid active material/electrolyte interface, and is defined as the solid phase potential in the solid active material side (ϕ s,l ) deducted by the electrolyte phase potential in the electrolyte side (ϕ e,m ) and the open-circuit potential (U), i.e.

$$ \eta_{\text{int}} = \phi_{{{\text{s,}}l}} - \phi_{{{\text{e,}}m}} - U. $$
(11)

The open-circuit potential U is a function of the local state of charge, θ [43].

$$ \begin{gathered} U_{\text{an}} = 0.722 + 0.1387\theta + 0.029\sqrt \theta - \frac{0.0172}{\theta } + \frac{0.0019}{{\theta^{1.5} }} \hfill \\ \, + 0.2808\exp \left( {0.90 - 15.0\theta } \right) - 0.7984\exp \left( {0.4465\theta - 0.4108} \right), \hfill \\ \end{gathered} $$
(12)
$$ U_{\text{ca}} = \frac{{ - 4.65 + 88.669\theta^{2} - 401.119\theta^{4} + 342.909\theta^{6} - 462.471\theta^{8} + 433.434\theta^{10} }}{{ - 1.0 + 18.933\theta^{2} - 79.532\theta^{4} + 37.311\theta^{6} - 73.083\theta^{8} + 95.96\theta^{10} }}. $$
(13)

The local state of charge, θ, is the Li concentration on the solid active material/ electrolyte interface evaluated by the maximum Li concentration that the solid active materials can accommodate, i.e.

$$ \theta = \frac{{c_{\text{s}} }}{{c_{\text{s,max}} }}. $$
(14)

Therefore, compared with the macroscopic models, the mesoscopic pore-scale model more accurately resolves the complex electrochemical reactions, which is prescribed to occur only at the interfaces of solid active materials and electrolyte.

(ii) Species conservation equations Diffusion is the major mechanism for transport of Li inside solid active materials. Therefore, the Li species conservation equation in solid active materials (not including the additives) is formulated as the following:

$$ \frac{{\partial c_{\text{s}} }}{\partial t} = \nabla \cdot \left( {D_{\text{s}} \nabla c_{\text{s}} } \right), $$
(15)

with boundary condition

$$ - D_{\text{s}} \frac{{A_{lm} }}{V}\left. {\frac{{\partial c_{\text{s}} }}{{\partial \varvec{n}}}} \right|_{\varOmega } = \frac{{j^{\text{Li}} (t)}}{F}, $$
(16)

where n denotes the external normal direction of solid active materials, Ω represents the interface of solid active material particles and electrolyte. Inside the two current collectors, there is no flux Li species, the boundary conditions for Li species transport thus yield as

$$ \left. {\frac{{\partial c_{\text{s}} }}{\partial x}} \right|_{x = 0} = \left. {\frac{{\partial c_{\text{s}} }}{\partial x}} \right|_{x = X} = 0. $$
(17)

The Li+ species conservation equation in electrolyte phase is formulated as the following:

$$ \frac{{\partial c_{\text{e}} }}{\partial t} = \nabla \cdot \left( {D_{\text{e}} \nabla c_{\text{e}} } \right) + \frac{{1 - t_{ + }^{ 0} }}{F}j^{\text{Li}} - \frac{{\varvec{i}_{\text{e}} \cdot \nabla t_{ + }^{ 0} }}{F}, $$
(18)

where i e is the current density vector and D e the mass diffusion coefficient of Li+ in the electrolyte. In this work, a constant value of the transference number of Li+ is assumed [21], and thus the last term on the right hand side of Eq. (18) is omitted. Inside the two current collectors, the flux density of Li+ species must be zero. Therefore, we obtain the boundary conditions for Li+ species transport as

$$ \left. {\frac{{\partial c_{\text{e}} }}{\partial x}} \right|_{x = 0} = \left. {\frac{{\partial c_{\text{e}} }}{\partial x}} \right|_{x = X} = 0. $$
(19)

(iii) Energy conservation equation The energy conservation equation is formulated as

$$ \rho C_{\text{P}} \frac{\partial T}{\partial t} = \nabla \cdot \left( {\lambda \nabla T} \right) + Q. $$
(20)

The heat generation rate Q is expressed as [21]

$$ Q = j^{\text{Li}} \left( {\phi_{\text{s}} - \phi_{\text{e}} - U + T\frac{\partial U}{\partial T}} \right) + \left( {\kappa_{\text{D}} \nabla \,\ln c_{\text{e}} \nabla \phi_{\text{e}} + \kappa \nabla \phi_{\text{e}} \nabla \phi_{\text{e}} } \right) + \sigma \nabla \phi_{\text{s}} \nabla \phi_{\text{s}} + Q_{\text{c}} . $$
(21)

The first term on the right hand side of Eq. (21) represents the electrochemical reactions heating, the second and third term come from the joule heating in the electrolyte and solid phases, respectively, and the last term represents the heat generation due to contact resistance and the resistance of SEI film. We consider only convective heat release to the environment and write the thermal boundary conditions as

$$ - \lambda \left. {\frac{\partial T}{\partial x}} \right|_{x = 0} = - \lambda \left. {\frac{\partial T}{\partial x}} \right|_{x = X} = h_{\text{E}} (T - T_{\text{E}} ), $$
(22)

where h E denotes the convective heat transfer coefficient and T E the temperature of environment.

In the above equations, all the physicochemical properties are location-dependent and component-wise. For example, the species transport of Li+ only occurs in the electrolyte, thus the component-wise Li+ diffusivity field is defined as

$$ D_{\text{e}} = \left\{ {\begin{array}{*{20}c} {D_{\text{e,0}} \left( {c_{\text{e}} ,T} \right),} & {\text{in electrolyte phase,}} \\ {0,} & {\text{in solid phase and current collector}}, \\ \end{array} } \right. $$
(23)

where D e,0 is the inherent diffusivity of the location component.

(iv) Thermal and electrochemical coupling To accommodate the thermal and electrochemical coupling effects, temperature-dependent physicochemical properties are considered and the temperature-dependence is described by the Arrhenius’ equation [44]

$$ \varPhi = \varPhi_{\text{ref}} \exp \left[ {\frac{{E_{{{\text{act,}}\varPhi }} }}{R}\left( {\frac{1}{{T_{\text{ref}} }} - \frac{1}{T}} \right)} \right], $$
(24)

where Φ is a general variable representing σ, κ, D s, D e, λ, etc., Φ ref is the value at a reference temperature T ref, and E act,Φ the activation energy of the evolution process of Φ, the magnitude of which indicates the sensitivity of Φ to temperature.

2.2.3 SPH formulation and methodology

In SPH, the quantities at time t is represented by a collection of N particles located at position r l (t) (l = 1, 2, …, N), as shown in Fig. 2. The “smoothed” value of any field quantity q(r, t) for the particle at spatial position r and at time t is a weighted sum of all contributions from its neighboring particles

$$ \left\langle {q(\varvec{r},t)} \right\rangle = \sum\limits_{m = 1}^{N} {\frac{{M_{m} }}{{\rho (\varvec{r}_{m} )}}} q\left( {\varvec{r}_{m} ,t} \right)w\left( {\left| {\varvec{r} - \varvec{r}_{m} } \right|,h} \right), $$
(25)

where M m and ρ(r m ) denote the mass and density of particle m, respectively. w(|r|, h) is the weight or smoothing function with h being the smoothing length. The smoothing kernel has an influencing distance of ah. The SPH formalism also provides a means to determine the gradient of q(r, t) as

$$ \left\langle {\nabla q(\varvec{r},t)} \right\rangle = \sum\limits_{m = 1}^{N} {\frac{{M_{m} }}{{\rho (\varvec{r}_{m} )}}} q\left( {\varvec{r}_{m} ,t} \right)\nabla w\left( {\left| {\varvec{r} - \varvec{r}_{m} } \right|,h} \right), $$
(26)

where \( \nabla w \) is the gradient of the interpolating kernel function centered at position r with respect to the location of neighboring particles, r m

$$ \nabla w\left( {\left| {\varvec{r} - \varvec{r}_{m} } \right|,h} \right) = \frac{{\varvec{r} - \varvec{r}_{m} }}{{\left| {\varvec{r} - \varvec{r}_{m} } \right|}}\frac{\partial w}{\partial r}. $$
(27)
Fig. 2
figure 2

A schematic of SPH interpolation

Full size image

Cleary et al. [45] gave the SPH form of the second order diffusion term. Taking the charge conservation equation for electrons in the solid phases, eq. (3), as an example, the electron diffusion term is expressed as

$$ \nabla \cdot (\sigma \nabla \phi_{\text{s}} ) = \sum\limits_{m} {V_{m} } \left( {\sigma_{l} + \sigma_{m} } \right)\phi_{{{\text{s,}}lm}} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }}, $$
(28)

where V m is the volume of SPH particle, V m  = M m /ρ(r m ), and the subscripts l and m are identifications of the local particle l and its neighboring particles m, respectively. r lm  = r l  − r m , w lm  = w(|r l  − r m |, h), and ϕ s,lm  = ϕ s(r l , t) − ϕ s(r m , t). In order to enhance SPH for the solution of transport problems with discontinuous conductivity, eq. (28) is upgraded by replacing the electronic conductivity of individual SPH particle [45] with their harmonic average, as

$$ \nabla \cdot (\sigma \nabla \phi_{\text{s}} ) = \sum\limits_{m} {4V_{m} \frac{{\sigma_{l} \sigma_{m} }}{{\sigma_{l} + \sigma_{m} }}\phi_{{{\text{s,}}lm}} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }}} . $$
(29)

By substituting Eq. (29) back into Eq. (3), the SPH form of the charge conservation equation in the solid active materials yields, as

$$ \sum\limits_{m} {4V_{m} } \frac{{\sigma_{l} \sigma_{m} }}{{\sigma_{l} + \sigma_{m} }}\phi_{{{\text{s,}}lm}} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }} - j^{\text{Li}} = 0. $$
(30)

Similarly, the SPH form of the charge conservation equation for Li+ in the electrolyte phase is formulated as

$$ \sum\limits_{m} {4V_{m} } \frac{{\kappa_{l} \kappa_{m} }}{{\kappa_{l} + \kappa_{m} }}\phi_{{{\text{e,}}lm}} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }} + \sum\limits_{m} {4V_{m} } \frac{{\kappa_{{{\text{D,}}l}} \kappa_{{{\text{D,}}m}} }}{{\kappa_{{{\text{D,}}l}} + \kappa_{{{\text{D,}}m}} }}\frac{2}{{c_{{{\text{e,}}l}} + c_{{{\text{e,}}m}} }}c_{{{\text{e,}}lm}} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }} + j^{\text{Li}} = 0, $$
(31)

where \( c_{{{\text{e,}}lm}} = c_{{{\text{e,}}l}} - c_{{{\text{e,}}m}} \), and \( \phi_{{{\text{e,}}lm}} = \phi_{{{\text{e,}}l}} - \phi_{{{\text{e,}}m}} \).

The relaxation iterative method is employed to solve the charge conservation equations in solid phases and electrolyte phase (Eqs. (30) and (31)). The iterative formulas are

$$ \left\{ {\begin{array}{*{20}c} {\phi_{\text{s}}^{n} = \omega_{\text{s}} \phi_{\text{s}}^{n} + \left( {1 - \omega_{\text{s}} } \right)\phi_{\text{s}}^{{n{ - 1}}} ,} \\ {\phi_{\text{e}}^{n} = \omega_{\text{e}} \phi_{\text{e}}^{n} + \left( {1 - \omega_{\text{e}} } \right)\phi_{\text{e}}^{{n{ - 1}}} ,} \\ \end{array} } \right. $$
(32)

where n is the number of iteration steps, ω s and ω e are the relaxation factors. If relaxation factors ω s and ω e take values greater than 1.0 and less than 2.0, over-relaxation is implemented, if they take values between 0 and 1.0, under-relaxation is implemented.

The species conservation equations in solid active materials and in electrolyte phase can evolve into the following SPH formulations, respectively.

$$ \frac{{\partial \left( {c_{{{\text{s,}}l}} } \right)}}{\partial t} = \sum\limits_{m} {4V_{m} } \frac{{D_{{{\text{s,}}l}} D_{{{\text{s,}}m}} }}{{D_{{{\text{s,}}l}} + D_{{{\text{s,}}m}} }}c_{{{\text{s,}}lm}} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }} + \frac{1}{F}j^{\text{Li}} , $$
(33)
$$ \frac{{\partial \left( {c_{{{\text{e,}}l}} } \right)}}{\partial t} = \sum\limits_{m} {4V_{m} } \frac{{D_{{{\text{e,}}l}} D_{{{\text{e,}}m}} }}{{D_{{{\text{e,}}l}} + D_{{{\text{e,}}m}} }}c_{{{\text{e,}}lm}} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }} + \frac{{1 - t_{ + }^{0} }}{F}j^{\text{Li}} , $$
(34)

where \( c_{{{\text{s,}}lm}} = c_{{{\text{s,}}l}} - c_{{{\text{s,}}m}} \). We take \( t_{ + }^{ 0} = 0. 3 6 3 \) in the present work.

The SPH form of the energy conservation equation is

$$ \left( {\rho C_{\text{P}} } \right)_{l} \frac{{\partial \left( {T_{l} } \right)}}{\partial t} = \sum\limits_{m} {4V_{m} } \frac{{\lambda_{l} \lambda_{m} }}{{\lambda_{l} + \lambda_{m} }}T_{lm} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }} + Q. $$
(35)

Here \( T_{lm} = T_{l} - T_{m} \). Contact resistance and the resistance of SEI film are neglected in the present work, then the heat generation rate Q is calculated by

$$ \begin{gathered} Q = j^{\text{Li}} \left( {\phi_{\text{s}} - \phi_{\text{e}} - U + T\frac{\partial U}{\partial T}} \right) + \left\{ {\left( {\sum\limits_{m} {4V_{m} \frac{{\kappa_{{{\text{D,}}l}} \kappa_{{{\text{D,}}m}} }}{{\kappa_{{{\text{D,}}l}} + \kappa_{{{\text{D,}}m}} }}\frac{{c_{{{\text{e,}}ml}} }}{{c_{{{\text{e,}}l}} + c_{{{\text{e,}}m}} }}\frac{{\partial w_{lm} }}{{\partial r_{l} }}} } \right)\left( {\sum\limits_{m} {V_{m} \phi_{{{\text{e,}}ml}} \frac{{\partial w_{lm} }}{{\partial r_{l} }}} } \right) + } \right. \hfill \\ \left. { \, \left( {\sum\limits_{m} {2V_{m} \frac{{\kappa_{l} \kappa_{m} }}{{\kappa_{l} + \kappa_{m} }}\phi_{{{\text{e,}}ml}} \frac{{\partial w_{lm} }}{{\partial r_{l} }}} } \right)\left( {\sum\limits_{m} {V_{m} \phi_{{{\text{e,}}ml}} \frac{{\partial w_{lm} }}{{\partial r_{l} }}} } \right)} \right\}{ + }\left( {\sum\limits_{m} {2V_{m} \frac{{\sigma_{l} \sigma_{m} }}{{\sigma_{l} + \sigma_{m} }}\phi_{{{\text{s,}}ml}} \frac{{\partial w_{lm} }}{{\partial r_{l} }}} } \right)\left( {\sum\limits_{m} {V_{m} \phi_{{{\text{s,}}ml}} \frac{{\partial w_{lm} }}{{\partial r_{l} }}} } \right), \hfill \\ \end{gathered} $$
(36)

where \( c_{{{\text{e,}}ml}} = c_{{{\text{e,}}m}} - c_{{{\text{e,}}l}} \), \( \phi_{{{\text{e,}}ml}} = \phi_{{{\text{e,}}m}} - \phi_{{{\text{e,}}l}} \) and \( \phi_{{{\text{s,}}ml}} = \phi_{{{\text{s,}}m}} - \phi_{{{\text{s,}}l}} \).

The arrangement of SPH particles, including boundary particles is illustrated in Fig. 3. The SPH particles (white) are uniformly distributed within the computational domain. The separation distance between neighboring particles is Δx in the x direction and Δy in the y direction; and Δx = Δy = ΔS. The inner node approach is used to position the SPH particles, meaning that the SPH particles adjoining the boundaries are not set on the boundary surfaces, but have half ΔS distance away from the respective adjoining boundary.

Fig. 3
figure 3

Illustration of the initial particle arrangement and boundary treatment

Full size image

The Leapfrog algorithm [46], due to its relatively low requirement on storage memory and high computational efficiency, is adopted to approximate the transient terms in Eqs. (33), (34) and (35). The time step Δt is determined by [47]

$$ \Delta t = \mathop {\hbox{min} }\limits_{l} \left( {\frac{{\zeta \left( {\rho C_{\text{P}} } \right)h^{2} }}{{\lambda_{l} }},\frac{{\zeta h^{2} }}{{D_{{{\text{e,}}l}} }},\frac{{\zeta h^{2} }}{{D_{{{\text{s,}}l}} }}} \right). $$
(37)

The calculation stability requires ζ to take a sufficiently small value. In the present work ζ = 0.1 is employed. In the present SPH implementation, the B-splines smoothing kernel is employed [48]. For two-dimensional problems, the kernel is defined as

$$ w_{lm} = \frac{15}{{7\pi h^{2} }} \times \left\{ {\begin{array}{*{20}c} {\frac{2}{3} - s^{2} + \frac{{s^{3} }}{2},} & {0 \le s < 1,} \\ {\frac{{\left( {2 - s} \right)^{3} }}{6},} & {1 \le s < 2} \\ {0,} & {s \ge 2,} \\ \end{array} } \right., $$
(38)

where s = |r lm |/h. The smoothing length h is specified to be ΔS in the present SPH implementation. Generally, a larger radius of support region ah can enhance the accuracy of calculation, while making the calculation far more costly. We take a = 2 and the interpolation of the field quantity with respect to one SPH particle extends to a circular region where 12 neighboring particle are included.

2.2.4 Boundary treatment

In the present SPH implementation, boundary conditions, particularly for the periodic boundary and flux boundary, require to be properly handled. The strategy of “ghost” particles is employed, as schematically shown in Fig. 3. Taking into account the influence region (a circle of radius 2 h) of the employed B-splines smoothing kernel, we arrange two layers of boundary “ghost” particles (black) outside of the computational domain and parallel to the corresponding boundary with the first layer half ΔS away from the boundary plane. These “ghost” particles have the same spatial separation ΔS as real SPH particles.

The periodic boundary condition is imposed in the y direction, meaning that the particles close to the top and bottom boundaries interact with each other in terms of Eqs. (30), (31), (33), (34) and (35) due to the periodical replication of the simulated domain. The real SPH particle and its image particle (“ghost” particle), such as particle A and particle B share the same physical properties. The major variables, i.e. ϕ s, ϕ e, c s, c e, and T for particle B are not calculated from the SPH smoothing interpolation, but artificially adjusted to have the same values of particle A throughout the calculation.

For the treatment of the flux boundary in the horizontal direction, we employ a novel approach [49], which implements a flux boundary condition as a source term added to the governing equation. For example, the charge conservation equation in solid phases is a steady-state second order diffusion equation. When solving it, the anode electrode is defined as a reference electrode (the solid phases potential takes value of −3.02 V). We evolve the constant current flux imposed at the right boundary plane into a source term, and thus the following governing equation yields.

$$ \sum\limits_{m} {4V_{m} } \frac{{\sigma_{l} \sigma_{m} }}{{\sigma_{l} + \sigma_{m} }}\phi_{{{\text{s,}}lm}} \frac{1}{{\left| {\varvec{r}_{lm} } \right|}}\frac{{\partial w_{lm} }}{{\partial r_{l} }} - j^{\text{Li}} + Q_{{{\text{f - }}\phi_{\text{s}} }} = 0, $$
(39)

where the last source term is calculated by

$$ Q_{{{\text{f - }}\phi_{\text{s}} }} = \left\{ {\begin{array}{*{20}c} {\frac{I}{A\Delta x},} & {\text{for the real SPH particles adjoining the right boundary,}} \\ {0,} & {{\text{for else SPH particles}} .} \\ \end{array} } \right. $$
(40)

The convective heat transfer boundary condition (Eq. 22) can be treated in the same way. The additional source term added to the energy conservation equation is expressed by

$$ Q_{{{\text{f - }}T}} \approx \left\{ {\begin{array}{*{20}c} {\frac{{2\lambda Nu(T - T_{\text{E}} )}}{Y\Delta x},} & {\text{for the real SPH particles adjoining boundaries,}} \\ {0,} & {{\text{for else SPH particles}} .} \\ \end{array} } \right. $$
(41)

Using this approach, all the non-zero flux and convective boundary conditions are transformed into zero-flux boundary conditions. It is easy to treat zero-flux boundary in SPH. We assign “ghost” particles (like particle D in Fig. 3) with the same properties as the corresponding real SPH particles (i.e. particle C in Fig. 3). To ensure zero flux boundary condition, major variables for the “ghost” particles are just artificially specified to have the same values as for the corresponding real SPH particles.

2.2.5 Solution procedure

The solution procedure of the mesoscopic pore-scale SPH model of multi-disciplinary transport processes coupled with electrochemical reactions in lithium-ion battery is described as follows.

  1. Step 1

    Position real and “ghost” particles, as shown in Fig. 3.

  2. Step 2

    Initialize variables and physicochemical properties.

  3. Step 3

    Implement boundary conditions for species and energy conservation equations.

  4. Step 4

    Search interacting particle pairs, calculate the kernel function and its derivative.

  5. Step 5

    Calculate physicochemical properties.

  6. Step 6

    Solve charge conservation equations, Eqs. (30) and (31) by internal iterations.

    1. Step 6(i)

      Implement boundary conditions for charge conservation equations.

    2. Step 6(ii)

      Update the values of electrolyte phase potential and solid phases potential.

    3. Step 6(iii)

      Repeat step 6(i) and step 6(ii) until certain convergence criteria is satisfied.

  7. Step 7

    Solve species and energy conservation equations by external iterations.

  8. Step 8

    Update time, repeat step 1 to step 7 until the end of discharge/charge process.

3 Simulation to discharge processes

The main goal of the present work is to demonstrate the viability and suitability of SPH. In this particular work, we focus solely on isothermal discharge processes, i.e. the energy conservation equation, eq. (35), is actually not solved. Later work will deal with the SPH simulations of discharge/charge processes and explore the coupled thermal-electrochemical effects.

We performed two-dimensional simulations with respect to a battery of idealized micro-pore structures in electrodes and separator, shown in Fig. 4. On the anode Cu current collector is pasted the graphite anode electrode, and on the cathode Al current collector the Li y−x CoO2 cathode electrode, in-between a separator is sandwiched. The anode electrode consists of four graphite particles connected by four additive particles, and the cathode electrode contains four Li y−x CoO2 particles connected by four additive particles. In Fig. 4, the x-direction represents the battery thickness direction (perpendicular to the electrode) and in the y-direction displays a periodically repeating unit of the lithium-ion battery.

Fig. 4
figure 4

Component distribution in an idealized single-cell battery. 0-Anode current collector (Cu); 1-Anode active materials (Li x C6); 2-Anode additives (CMC + SBR etc.); 3-Electrolyte (Lithium salts); 4-Polymer matrix of separator; 5-Cathode additives (SP + PVDF etc.); 6-Cathode active materials (Li y−x CoO2); 7-Cathode current collector (Al)

Full size image

The geometry of the simulated battery is a 60 μm × 21 μm rectangular domain, i.e. X = 60 μm (X = L an + L sep + L ca + 2L c), Y = 21 μm (see Fig. 4). The SPH particle separation is set as, ΔS = Δx = Δy = 1 μm. Therefore, including the “ghost” particles (see Fig. 3), totally 64 × 25 (= 1600) SPH particles are used in the simulations. Initially, T = T 0 = 25 °C, c e = c e,0, c s = c s,0. The values of major parameters and important constitutive relationships are summarized in Table 1. The physicochemical properties and constitutive relationships are either directly taken from literatures or best estimates based on available literature data.

Table 1 Parameters used in the simulations
Full size table

We simulated 1 and 2 C discharge processes. The predicted cell voltages are compared in Fig. 5. As evident by Fig. 5, the cell voltages continuously decrease with the increase of discharge time. The slower discharge process gives out higher discharge voltage. If the cut-off voltage is set at 3.2 V, the discharge durations are 3616 and 1815 s for 1 and 2 C processes, respectively. That is to say, the two discharge processes both can fully discharge the total energy stored. No evident transport limiting presents due to the very thin electrodes and the idealized microstructures assumed for the electrodes and separator (see Fig. 4).

Fig. 5
figure 5

Evolution of battery voltage during 1 and 2 C discharge processes

Full size image

In order to test the effect of the number of SPH particles on the simulation results, the dependence of the simulation results on the SPH particle separation was studied. For the 1 C discharge process, we performed two more simulations with total SPH particles (including the ghost particles) being 792 (i.e., ΔS = 1.5 μm) and 5704 (i.e., ΔS = 0.5 μm), respectively. We examined and compared the predicted cell voltage values by the three simulations. A Richardson extrapolation, which assumed the predicted voltage was linearly dependent on 1/(total SPH particles), was computed. In the regime of Richardson extrapolation, once the linear relationship is determined, it is easy to extrapolate the accurate value of a variable just by simply letting the total SPH particles be infinite, i.e. 1/(total SPH particles) = 0. As a result, the predicted voltage based on 1600 SPH particles only deviates by 6 % at most to the Richardson extrapolated accurate value.

Benefiting from the mesoscopic treatment, we can prescribe electrochemical reactions to occur exactly at the solid active materials/electrolyte interfaces. In numerical modeling, only the numerical element (SPH particle) adjoining the solid active materials/electrolyte interfaces are imposed with a transfer current density. Figure 6 displays the calculated transfer current density distribution in the interfacial solid active material particles for both the two discharge processes at different discharge stages. Under a same level of the depth of discharge, the 1 C discharge process generally shows approximately 1/2 magnitude of transfer current density compared with the 2 C discharge process, and the two cases have similar pattern of transfer current density distribution, indicating again, no evident transport limiting phenomena are present.

Fig. 6
figure 6figure 6

The transfer current density (A m−3) distribution during different stages of the 1 and 2 C discharge processes. a 2 C, t = 10 s and b 1 C, t = 20 s are the early period of the discharge process; c 2 C, t = 830 s and d 1 C, t = 1660 s are the mid-term of the discharge process; e 2 C, t = 1650 s and f 2 C, t = 3300 s are the late period of the discharge process

Full size image

The idealized battery has four solid active material particles in each electrode. Seen from Fig. 6, for both the anode and cathode electrode, the three solid active material particles that are closer to the separator have similar transfer current distribution while at the surface of the rest one solid active material particle that is interfacing the current collector the transfer current density is generally larger. This can also be ascribed to the idealized electrode microstructures assumed. As the solid active material particle interfacing the current collector has smaller contact interface to the electrolyte, intuitively, to achieve an approximately synchronous discharge process with the other solid particles the transfer current density imposed at its surface must be of larger magnitude.

The pore-scale SPH model is able to give detailed insights into the multi-disciplinary transport process inside the lithium-ion battery as well. As the calculated species concentration and the charge potential distributions of the 1 and 2 C discharge processes are similar, we only present the results from the 1 C discharge process hereinafter.

In Fig. 7 is displayed the calculated Li concentration distribution inside solid active materials at discharge time instants of 20, 1660 and 3300 s, respectively. As evident by Fig. 7, the Li concentration continuously decreases in the anode electrode and increases in the cathode electrode with the discharge process. Moreover, for a same solid active material particle, an evident Li concentration difference is seen if one compares the Li concentration at its surface with that in the interior. These are expected results as the electrochemical reactions only occur at the solid active materials/electrolyte interfaces. The anode electrochemical reaction consumes Li atoms at the surfaces of graphite particles and produce Li+, which dissolve and transport in the electrolyte, penetrate the separator, and finally reach the cathode. The cathode electrochemical reaction changes the Li+ in the electrolyte at the Li y−x CoO2 particle/electrolyte interfaces into Li atoms, which insert into the Li y−x CoO2 particles. We can see additionally from Fig. 7 that solid active material particles locating at different positions along the battery thickness direction, particularly the solid active material particle interfacing the current collector compared with the other three solid active material particles, are having different Li concentration distributions, indicating the electrochemical reactions are unevenly occurring (Refer to Fig. 6).

Fig. 7
figure 7

Li concentration (mol m−3) distribution in solid active materials during 1 C discharge process. a t = 20 s; b t = 1660 s; c t = 3300 s

Full size image

Further, we averaged the Li concentration in a same solid active material particle at discharge time instants of 20, 1660 and 3300 s, respectively, and the results were plotted in Fig. 8. For both the anode and cathode electrode, the averaged Li concentrations in the three solid active material particles that closer to the separator are of approximately same magnitude and evolution tendency; the averaged Li concentration in the solid active material particle that interfacing the current collector has some difference to that in the other three solid active material particles.

Fig. 8
figure 8

The averaged Li concentration (mol m−3) in solid active material particles during 1 C discharge process

Full size image

In Fig. 9 we depict the calculated Li+ concentration distribution in electrolyte at discharge time instants of 20, 1660 and 3300 s, respectively. Upon loading of discharge current, the anode electrochemical reaction produces Li+, which are dissolved in the electrolyte and transported to cathode; the cathode consumes Li+ in the electrolyte and turns them into Li atoms to insert into the Li y−x CoO2 particles. We see from Fig. 9 that upon loading a Li+ concentration gradient is established in the electrolyte, which drives the Li+ to transport from anode to cathode. We see also the Li+ concentration distribution show little change throughout the discharge process. This is due to no transport limiting presents.

Fig. 9
figure 9

Li+ concentration distribution (mol m−3) in electrolyte during 1 C discharge process. a t = 20 s; b t = 1660 s; c t = 3300 s

Full size image

During discharge process, electrons produced in the anode electrode are transported to the cathode electrode by passing the external circuit. Inside the battery, a solid phase potential distribution will be established to support the flow of electrons. Figure 10 presents the calculated solid phase potential at three discharge time instants: 20, 1660 and 3300 s. In anode, the solid phase potential decreases towards the separator side, which drives the electrons to flow from the anode electrode to the current collector. In cathode, the solid phase potential increases towards the separator side, meaning the electrons are flowing from the cathode current collector to the electrode. In anode, the solid phase potential drop is ~0.06 mV, much smaller than that in cathode, which is ~1.2 mV. This is due to the much smaller electron transport resistance in anode. As the reference potential is set in the anode current collector, the magnitude of solid phase potential in anode is almost unchanged throughout the discharge process. However, in cathode the solid phase potential is decreasing with the discharge process, meaning a diminishing output voltage (see Fig. 5).

Fig. 10
figure 10

Solid phase potential (V, vs. Li/Li+) during 1 C discharge process. a t = 20 s; b t = 1660 s; c t = 3300 s

Full size image

Likewise, during discharge process charge transport of Li+ needs a distribution of phase potential to be established in electrolyte. Figure 11 shows the calculated electrolyte phase potential distribution at discharge time instants of 20, 1660 and 3300 s. The electrolyte phase potential decreases along the battery thickness direction, indicating Li+ charge transport from the anode electrode to the cathode electrode. The electrolyte phase potential is diminishing with the discharge time, while the total electrolyte phase potential drop along the battery thickness direction remains within 1.8–2.2 mV, no significant changes throughout the whole discharge process.

Fig. 11
figure 11

Electrolyte phase potential (V, vs. Li/Li+) during 1 C discharge process. a t = 20 s; b t = 1660 s; c t = 3300 s

Full size image

4 Concluding remarks

We have implemented the SPH numerical technique for modeling the inter-coupled multi-disciplinary processes in lithium-ion batteries. The developed SPH model is based on microstructures of electrodes and separator, enabling mesoscopic pore-scale modeling and analyses of lithium-ion battery processes.

As a first step to test the developed model, we simulated two isothermal discharge processes of 1 and 2 C discharge rates, respectively, with respect to a battery of two-dimensional idealized micro-pore structure in the electrodes and separator. Results indicate that this model is not only able to predict macroscopic quantities such as the output voltage, but also capable of giving detailed distribution information of all the primary and participating quantities, including solid/electrolyte phase potential, Li+ concentration in electrolyte, Li concentration in solid active materials, and transfer current density at the solid active materials/electrolyte interfaces, etc. The reasonability and capability of the developed SPH model are thus preliminarily demonstrated as well.

Next step we will conduct model validation/calibration work, further the model development and apply the model to simulate more realistic battery discharge/charge processes. Future simulations will fully account for the three-dimensional micro- structural configurations in electrodes and separator [50, 51], and the non-isothermal effects etc. However, there exist still grand challenges, such as, the very short time step (calculated by Eq. (37)) due to the constraint imposed by the energy equation calculation, and the extremely intensive computational cost due to considering three-dimensional geometry. Possible solutions to these challenges include the development of implicit SPH iteration algorithm and the parallelization of three-dimensional computational code etc.