Characteristic Prediction and Temperature-Control Strategy under Constant Power Conditions for Lithium-Ion Batteries

Abstract

Accurate characteristic prediction under constant power conditions can accurately evaluate the capacity of lithium-ion battery output. It can also ensure safe use for new-energy vehicles and electrochemical energy storage. As the battery voltage continues to drop under constant power conditions, the battery current output will accordingly increase, which brings a risk of thermal runaway in instances of weak heat dissipation. Therefore, knowing how to control the battery temperature is very critical for safe use. At present, the model-based method for characteristic prediction and temperature control has been used by most scholars, and that is also the key to this method. This work firstly extends a cell model to a pack-based electrochemical two-dimensional thermal coupling model, considering the heterogeneity of different cells inside the pack, and obtains the model parameters for a prismatic lithium-ion battery with a rated capacity of 42 Ah. Characteristic prediction under constant power conditions is then conducted based on an iterative solution method. Validations of characteristic prediction indicate the convenience of the developed models, with average absolute errors of voltage and temperature less than 36 mV and 0.4 K, respectively, and power error less than 0.005%. Finally, two model-based temperature feed-forward control strategies with lower cooling costs and shorter prediction times were developed based on the battery characteristic predictions, which leaves room for further controller development.

1. Introduction

With the increasingly serious problems of environmental pollution and energy shortages, lithium-ion batteries, as one of the clean new energy sources, have been widely used in many industrial fields such as new-energy vehicles, energy storage, etc. It is necessary to predict battery characteristics under constant power conditions for their safe use, due to the potential of overheating.

Currently, the methods for battery power estimation mainly include model-based methods and machine-learning-based methods. Liu et al. [1] used a fractional equivalent circle model (ECM) to estimate the state of power (SOP) and used small-order calculus with state of charge (SOC), battery voltage, and current as constraints. It could be seen from the experiment that the maximum relative error of SOP estimation results was 1.34%. Feng et al. [2] used a novel ECM which added moving average noise to a resistor–capacitor circuit model to accurately capture battery dynamics. A recursive extended least squares algorithm was used to identify the ECM parameters online and showed high accuracy in the experiment. Waag et al. [3] proposed a nonlinear model-based method for predicting the available power in li-ion battery packs using the current dependence of the battery’s resistance. Accurate power prediction was possible at room lower temperatures and where the current dependence was larger. Pan et al. [4] proposed a multi-constraint power capacity forecasting method, which could overcome the shortcomings of three single constraints, including voltage, current, and SOC, and an extended Kalman filter algorithm was subsequently used for power state prediction. Most of the model-based methods are predicted by the ECM of the power battery, and the prediction effects of these methods depend on the model used. Guo et al. [5] experimentally studied the polarization characteristics of the battery in a series of decreasing pulse tests; constructed a feedforward neural network with inputs of SOC, discharge rate, and pulse running time, and characterized the polarization voltage by modeling the polarization resistance excited by the current; and then proposed a data-model fusion method to make accurate online SOP estimation in a prediction window of 30 s to 120 s. This machine-learning-based method had high accuracy but required a large amount of data to train the algorithm.

As the battery voltage gradually decreases with the descending of the SOC and the voltage drops sharply at the end of discharge, the battery current will gradually increase to meet the constant power output requirement. Due to the tight arrangement of the battery pack, there is a risk of thermal runaway under poor heat dissipation conditions. It is thus necessary to predict the power characteristics of the battery in advance and control the temperature of the battery pack. Current methods for heat dissipation mainly include air cooling [6], liquid cooling [7], phase-change material heat dissipation [8], heat pipe heat dissipation [9], and so on.

As there exists an unpredictable and uncertain time delay before a temperature-control system functions, a reasonable temperature-control strategy is thus essential to realizing rapid regulation of battery temperature. General battery system temperature-control strategies include: PID-based control, fuzzy-algorithm-based control, model-based predictive control, and coupling control in several ways. Cen et al. [10] used a PID algorithm to design an air-conditioning system for an electric vehicle to accomplish air circulation in the vehicle and the battery pack. Extensive experimental studies indicated that the optimized values could keep the maximum temperature difference below 2 °C at discharge rates of 0.5 C and 1 C. The control strategy was easy to implement, but its response speed was slow. Afzal et al. [11] developed a new technique for the multi-objective optimization of battery system thermal management using a hybrid genetic algorithm. Jiaqiang E. et al. [12] used the fuzzy grey relational theory and the orthogonal experimental design method to discuss the influence of the optimal forced air-cooling model parameters on the temperature-control ability of the battery module. Fuzzy-based control had strong robustness; however the required expert experience is difficult to obtain and the fuzzification of the collected information would probably lead to a decrease in accuracy. Zhu et al. [13] proposed a model-based predictive control strategy based on finite sets that could simultaneously minimize the energy consumption of cooling or heating. Model-based predictive control is suitable for the control of nonlinear and hysteretic systems, and has the capabilities of advanced prediction, time-varying online iterative optimization, and accurate output feedback correction, and its applicability is related to the accuracy of the battery model.

Currently, the commonly used models include ECM and the electrochemical model. Electrochemical models are favored by many scholars because of their high accuracy in battery state estimation at various conditions [14,15,16]. J. Newman et al. first established a pseudo-two-dimensional (P2D) model [17]. Its accuracy and versatility are good, but the calculation process is complicated and difficult to implement. Some scholars have added descriptions of temperature distribution onto the P2D model and others have proposed thermal coupling models [18,19]. The developed model could be used to optimize the design of thermal management systems and perform aging analysis [20,21,22,23,24]. It could also be used to analyze the battery’s internal heat generation and simulate the surface temperature change of different materials [25,26,27,28]. Romero-Becerril et al. [29] developed a single-particle (SP) model, which used the behavior of a single active particle to substitute for the chemical reaction behavior of the entire electrode. On the basis of the SP model, Luo Weilin et al. [30] proposed an improved SP model, which comprehensively balanced the accuracy and operation time. How to ensure simulation accuracy at a full-charge state operating range and realize rapid calculation becomes the key point in terms of characteristic prediction and temperature control.

This paper’s contribution is as follows: (1) an electrochemical-two-dimensional thermal coupling model is developed for a typical battery pack configuration, considering the differences among single cells, (2) a model-driven battery characteristic prediction algorithm is proposed to predict the voltage, current, and temperature of the battery pack under constant power conditions, and (3) two effective temperature-control strategies are proposed to keep the temperature within a certain range. The rest of this work is as follows: Section 2 introduces the battery pack modeling method and characteristic prediction method under constant power conditions for a battery pack; Section 3 introduces the experiments and verifies the accuracy of the prediction algorithm; Section 4 introduces the model-based thermal control strategy, followed by conclusions in Section 5.

2. Battery Model and Characteristic Prediction Method

2.1. Modeling for a Single Cell

The active positive and negative electrode materials of a lithium-ion battery are regarded as composed of spherical active particles. During the use of the battery, there are several physical and chemical processes: solid-phase diffusion, electrochemical reactions on particle surfaces, liquid-phase diffusion and internal heat generation, heat transfer from the interior of the battery to the outer shell, and heat dissipation on the battery surface. To simplify the calculation of the electrochemical characteristics of the battery, this work assumes that the reaction distribution at each position inside the battery is uniform. Based on Fick’s second law, the change law of solid-phase lithium-ion concentration in the direction of the radius of active battery material is established; the concentration polarization is described by parameters of P_{coni, i = n,p} and parameters of τ_e that reflect the degree of polarization and the time length of the establishment of the concentration difference, which participate in the calculation of the polarization overpotential represented by η_con; the difficulty of the electrochemical reaction inside the battery is described by parameters of the Pact polarization coefficient and it participates in the calculation of the reaction polarization overpotential, represented by η_act; the ohmic polarization overpotential, represented by η_ohm, follows Ohm’s law, and then the battery terminal voltage U_app can be calculated from η_con, η_act, η_ohm, and the open circuit voltage, which is represented by E_ocv.

To quickly predict the change of battery temperature characteristics, this work assumes that the temperature distribution inside the battery is uniform. The heat-generation rate equation can be obtained based on the Bernardi heat generation model. In this paper, the thermal resistance model is used to describe the heat exchange process between the electrode winding body and the battery shell, as well as the battery shell and the external environment. Due to the thermal anisotropy of the prismatic battery, the temperature gradient distribution exists in two directions: x and y, which are shown in Figure 1.

Change in the internal temperature affects the electrochemical reaction processes, and it will affect the model parameters for electrochemical behaviors. It is well-known that the influence of the thermal change on the parameters can be described by the Arrhenius equation. Open-circuit voltage E_ocv and other overpotentials η_act and η_con are dependent on temperature. The Nernst equation correction equation for E_ocv is shown in Equation (1), where T_ref is the reference temperature and is the open-circuit potential in electrochemistry:

$$E_{\mathit{ocv}}{=E}_{\mathit{ocv}}^{\mathit{ref}}{+(T - T}_{\mathit{ref}})\frac{\mathrm{d}{E}_{\mathit{ocv}}}{\mathrm{d}T}$$

(1)

Detailed descriptions of simplified electrochemical-thermal coupling models have been introduced in our previous work [29]. The main numerical equations of the model are shown in

Table 1. Numerical equations of the model.

Correlation Mechanism	Equations
Terminal voltage	${\mathit{U}}_{\mathit{app}}{\left(\mathit{k}\right)=\mathit{E}}_{\mathit{ocv}}{\left(\mathit{k}\right) - \eta }_{\mathit{con}}{\left(\mathit{k}\right) - \eta }_{\mathit{act}}{\left(\mathit{k}\right) - \eta }_{\mathit{ohm}}\left(\mathit{k}\right)$
Open circuit potential correction	${\mathit{E}}_{\mathit{ocv}}{\left(\mathit{k}\right)=\mathit{E}}_{\mathit{ocv}}^{\mathit{ref}}{\left(\mathit{k}\right)+\left(\mathit{T}\right(\mathit{k}) - \mathit{T}}_{\mathit{ref}})\frac{\mathrm{d}{E}_{\mathit{ocv}}}{\mathrm{d}T}$
Basic working principle	${\mathit{x}}_{\mathit{avg}}{\left(\mathit{k}\right)=\mathit{x}}_0 - {\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_{\mathrm{k}}}{I\mathrm{d}t/\mathit{Q}}_{\mathit{n}}$, ${\mathit{x}}_{\mathit{surf}}{\left(\mathit{k}\right)=\mathit{x}}_{\mathit{avg}}\left(\mathit{k}\right) - \Delta \mathit{x}\left(\mathit{k}\right)$, $\Delta \mathit{x}\left(\mathit{k}\right)=\Delta {\mathit{x}}_1\left(\mathit{k}\right)+\frac{2}{7}{\tau }_{\mathit{n}}{\mathit{I}\left(\mathit{k}\right)/\mathit{Q}}_{\mathit{n}}$, ${\mathit{y}}_{\mathit{avg}}{\left(\mathit{k}\right)=\mathit{y}}_0+{\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_{\mathrm{k}}}{I\mathrm{d}t/\mathit{Q}}_{\mathit{p}}$, ${\mathit{y}}_{\mathit{surf}}{\left(\mathit{k}\right)=\mathit{y}}_{\mathit{avg}}\left(\mathit{k}\right)+\Delta \mathit{y}\left(\mathit{k}\right)$, $\Delta \mathit{y}\left(\mathit{k}\right)=\Delta {\mathit{y}}_1\left(\mathit{k}\right)+\frac{2}{7}{\tau }_{\mathit{p}}{\mathit{I}\left(\mathit{k}\right)/\mathit{Q}}_{\mathit{p}}$, ${\mathit{E}}_{\mathit{ocv}}{\left(\mathit{k}\right)=\mathit{U}}_{\mathit{p}}\left[{\mathit{y}}_{\mathit{surf}}\left(\mathit{k}\right)\right]{ - \mathit{U}}_{\mathit{n}}\left[{\mathit{x}}_{\mathit{surf}}\left(\mathit{k}\right)\right]$
Solid diffusion	$\Delta {\mathit{x}}_1(\mathit{k}+1)=\Delta {\mathit{x}}_1\left(\mathit{k}\right)+\frac{1}{{\tau }_{\mathit{n}}}[\frac{{12\tau }_{\mathit{n}}\mathit{I}\left(\mathit{k}\right)}{{7\mathit{Q}}_{\mathit{n}}} - \Delta {\mathit{x}}_1\left(\mathit{k}\right)]\Delta \mathit{t}, \Delta {\mathit{y}}_1(\mathit{k}+1)=\Delta {\mathit{y}}_1\left(\mathit{k}\right)+\frac{1}{{\tau }_{\mathit{p}}}[\frac{{12\tau }_{\mathit{p}}\mathit{I}\left(\mathit{k}\right)}{{7\mathit{Q}}_{\mathit{p}}}-\Delta {\mathit{y}}_1\left(\mathit{k}\right)]\Delta \mathit{t}$
Liquid phase diffusion	$\Delta c_{\mathit{n}}(\mathit{k}+1)=\Delta c_{\mathit{n}}\left(\mathit{k}\right)+\frac{1}{{\tau }_{\mathit{e}}}\left[{\mathit{P}}_{\mathit{conn}}\mathit{I}\left(\mathit{k}\right) - \Delta c_{\mathit{n}}\left(\mathit{k}\right)\right]\Delta \mathit{t}, \Delta c_{\mathit{p}}(\mathit{k}+1)=\Delta c_{\mathit{p}}\left(\mathit{k}\right)+\frac{1}{{\tau }_{\mathit{e}}}\left[{\mathit{P}}_{\mathit{conp}}- \Delta c_{\mathit{p}}\left(\mathit{k}\right)\right]\Delta \mathit{t},$ ${\eta }_{\mathit{con}}\left(\mathit{k}\right)=\frac{2\mathit{RT}\left(\mathit{k}\right)}{\mathit{F}}\left({1 - \mathit{t}}_{+}\right)\ln \left[\frac{{\mathit{c}}_0+\Delta c_p\left(\mathit{k}\right)}{{\mathit{c}}_0 - \Delta c_{\mathit{n}}\left(\mathit{k}\right)}\right]$
Reactive polarization	${\eta }_{\mathit{act}}\left(\mathit{k}\right)=\frac{2\mathit{RT}\left(\mathit{k}\right)}{\mathit{F}}\left[\ln \right(\sqrt{{\mathit{m}}_{\mathit{n}}^2\left(\mathit{k}\right)+1}{+\mathit{m}}_{\mathit{n}}\left(\mathit{k}\right)) - \ln (\sqrt{{\mathit{m}}_{\mathit{p}}^2\left(\mathit{k}\right)+1}{+\mathit{m}}_{\mathit{p}}\left(\mathit{k}\right)\left)\right]$ ${\mathit{m}}_{\mathit{p}}\left(\mathit{k}\right)=\frac{1}{{6\mathit{Q}}_{\mathit{p}}{\mathit{c}}_0^{0.5}}\frac{1}{{{(1 - \mathit{y}}_{\mathit{surf}}\left(\mathit{k}\right))}^{0.5}{{(\mathit{y}}_{\mathit{surf}}\left(\mathit{k}\right))}^{0.5}}{\mathit{P}}_{\mathit{act}}\mathit{I}\left(\mathit{k}\right)$,${ \mathit{m}}_{\mathit{n}}=\frac{1}{{6\mathit{Q}}_{\mathit{n}}{\mathit{c}}_0^{0.5}}\frac{1}{{{(1- \mathit{x}}_{\mathit{surf}}\left(\mathit{k}\right))}^{0.5}{{(\mathit{x}}_{\mathit{surf}}\left(\mathit{k}\right))}^{0.5}}{\mathit{P}}_{\mathrm{act}}\mathit{I}\left(\mathit{k}\right)$
Ohmic polarization	${\eta }_{\mathit{ohm}}{=\mathit{R}}_{\mathit{ohm}}\mathit{I} \left(\mathit{k}\right)$
Heat production	$\dot{\mathit{Q}}\left(\mathit{k}\right)=-\mathit{I}\left(\mathit{k}\right)\mathit{T}\left(\mathit{k}\right)\frac{\mathrm{d}{E}_{\mathit{ocv}}}{\mathrm{d}T}{+\mathit{I}\left(\mathit{k}\right)(\mathit{E}}_{\mathit{ocv}}{\left(\mathit{k}\right) - \mathit{U}}_{\mathit{app}}\left(\mathit{k}\right))$
Internal temperature	${\mathit{T}}^{\mathit{k}+1}{=\mathit{T}}^{\mathit{k}}+\frac{\mathsf{\Delta }\mathit{t}}{{\mathit{m}}_{\mathit{roll}}{\mathit{C}}_{\mathit{p}}}\left(\dot{\mathit{Q}}\right(\mathit{k}) - \left(\frac{1}{{\mathit{R}}_{\mathit{condx}}}+\frac{{8\mathit{V}\lambda }_{\mathit{x}}}{{\mathit{l}}_{\mathit{x}}^2}\right){(\mathit{T}}^{\mathit{k}}{ - \mathit{T}}_{\mathit{surfx}}^{\mathit{k}}) - (\frac{1}{{\mathit{R}}_{\mathit{condy}}}+\frac{{8\mathit{V}\lambda }_{\mathit{y}}}{{\mathit{l}}_{\mathit{y}}^2}{\left)\right(\mathit{T}}^{\mathit{k}}{ - \mathit{T}}_{\mathit{surfy}}^{\mathit{k}}\left)\right)$
Shell temperature	${\mathit{T}}_{\mathit{surfx}}^{\mathit{k}+1}{=\mathit{T}}_{\mathit{surfx} }^{\mathit{k}}+\frac{\mathsf{\Delta }\mathit{t}}{{\mathit{m}}_{\mathit{can}}{\mathit{C}}_{\mathit{can}}}(\frac{1}{{\mathit{R}}_{\mathit{condx}}}{(\mathit{T}}^{\mathit{k}+1}{- \mathit{T}}_{\mathit{surfx}}^{\mathit{k}}) - \frac{{\mathit{Vh}}_{\mathit{x}}}{{\mathit{l}}_{\mathit{y}}}{(\mathit{T}}_{\mathit{surfx}}^{\mathit{k}}{ - \mathit{T}}_{\mathit{a}}\left)\right)$ ${\mathit{T}}_{\mathit{surfy}}^{\mathit{k}+1}{=\mathit{T}}_{\mathit{surfy}}^{\mathit{k}}+\frac{\mathsf{\Delta }\mathit{t}}{{\mathit{m}}_{\mathit{can}}{\mathit{C}}_{\mathit{can}}}(\frac{1}{{\mathit{R}}_{\mathit{condy}}}{(\mathit{T}}^{\mathit{k}+1}{ - \mathit{T}}_{\mathit{surfy}}^{\mathit{k}}) - \frac{{\mathit{Vh}}_{\mathit{x}}}{{\mathit{l}}_{\mathit{y}}}{(\mathit{T}}_{\mathit{surfy}}^{\mathit{k}}{ - \mathit{T}}_{\mathit{a}}\left)\right)$

.

2.2. Parameter Identification

To realize the model simulation, the model parameters needed to be obtained in advance. Except for inherent characteristic parameters such as weight and length, model parameters which appear in Table 1 needed to be obtained individually via the identification method.

There was an obvious time difference among the electrochemical processes of the battery, and the response degrees of different processes caused by the same current excitation were also different, so the decoupling of each part could be completed [30]. The current used for identification is shown in Figure 2.

The thermal parameters that needed to be determined mainly included C_p, λ_x, λ_y, R_condx, and R_condy. C_p, λ_x, and λ_y can be calculated according to the ratio of different components of the battery. The equivalent specific heat capacity C_p,battery can be calculated by Equation (2), where m_i and c_p,i are the mass and specific heat capacity of the corresponding component material. The calculation of λ_x and λ_y was performed as shown in Equation (3), where L_a and_i are the thickness of the component material in the corresponding axial direction, and λ_a and _iare the thermal conductivity of the corresponding component material.

$${\mathit{C}}_{\mathit{p},\mathit{battery}}=\frac{{\sum }_{\mathit{i}=1}^{\mathit{n}}{\mathit{C}}_{\mathit{p},\mathit{i}}{\mathit{m}}_{\mathit{i}}}{{\sum }_{\mathit{i}=1}^{\mathit{n}}{\mathit{m}}_{\mathit{i}}}$$

(2)

$${\lambda }_{\mathit{a}}=\frac{{\sum }_{\mathit{i}=1}^{\mathit{n}}{\mathit{L}}_{\mathit{a},\mathit{i}}}{{\sum }_{\mathit{i}=1}^{\mathit{n}}{(\mathit{L}}_{\mathit{a},\mathit{i}}{/\lambda }_{\mathit{a},\mathit{i}})}, \mathit{a}=\mathit{x}, \mathit{y}$$

(3)

R_condx and R_condy were used to characterize the heat transferred from the inside to the outside in two directions, respectively. They can be obtained iteratively through Equation (4). The initial values of R_condx and R_condy are 0.01.

$${\mathit{R}}_{\mathit{condx}}=\frac{{\mathit{T}}^{\mathit{k}+1}{-\mathit{T}}_{\mathit{surfx}}^{\mathit{k}}}{\left({\mathit{T}}_{\mathit{surfx}}^{\mathit{k}+1}{-\mathit{T}}_{\mathit{surfx}}^{\mathit{k}}\right)\frac{{\mathit{m}}_{\mathit{can}}{\mathit{C}}_{\mathit{can}}}{\mathsf{\Delta }\mathit{t}}+\frac{\mathit{V}{h}_x}{{\mathit{l}}_{\mathit{y}}}{(\mathit{T}}_{\mathit{surfx}}^{\mathit{k}}{-\mathit{T}}_{\mathit{a}})}, {\mathit{R}}_{\mathit{condy} }=\frac{{\mathit{T}}^{\mathit{k}+1}{-\mathit{T}}_{\mathit{surfy}}^{\mathit{k}}}{\left({\mathit{T}}_{\mathit{surfy}}^{\mathit{k}+1}{-\mathit{T}}_{\mathit{surfy}}^{\mathit{k}}\right)\frac{{\mathit{m}}_{\mathit{can}}{\mathit{C}}_{\mathit{can}}}{\mathsf{\Delta }\mathit{t}}+\frac{\mathit{V}{h}_y}{{\mathit{l}}_{\mathit{x}}}{(\mathit{T}}_{\mathit{surfy}}^{\mathit{k}}{-\mathit{T}}_{\mathit{a}})}$$

(4)

2.3. Modeling a Battery Pack

Taking a three-series battery-pack structure as the modeling object, the electrical connection and the thermal connection of each single lithium-ion battery in the battery pack were studied. The series of functional modules is shown in Figure 3.

The electrochemical reactions of the cells are independent of each other; the external characteristic behaviors of the battery pack can be regarded as the linear superposition of the cells connected in series, as shown in Equation (5), where U_{app, i}, and _{i = total, 1,2,3} are the total voltage of the battery pack and each cell, and I_i and _{i = total, 1,2,3} are the total current and the current of each cell. The behaviors of the external characteristics of the cell were determined by the model parameters, and the corresponding model parameter sets of different cells are also different. The parameter acquisition of each cell was the same as described above.

$${\mathit{U}}_{\mathit{app},\mathit{total}}{=\mathit{U}}_{\mathit{app},1}{+\mathit{U}}_{\mathit{app},2}{+\mathit{U}}_{\mathit{app},3}{, \mathit{I}}_{\mathit{total}}{=\mathit{I}}_1{=\mathit{I}}_2{=\mathit{I}}_3$$

(5)

The multi-physics model of the battery pack uses a thermal resistance network method to establish the thermal connection between the cells, assuming that the temperature fields inside the cells are independent of each other and the thermal connection between the cells is established through the thermal gap between the cells. The thermal resistance network diagram of the battery pack is shown in Figure 4. In the process of battery grouping, the batteries are assembled by bonding their larger surfaces with the purpose of improving space utilization. Gap in Figure 4 is the air interstice inside the battery pack.

The symmetrical end faces have the same heat transfer coefficient for a single cell. For a pack as in Figure 4, the heat transfer coefficients h of the inner and outer faces for cell 1 and cell 3 are different, while the h of the two inner faces for cell 2 remains identical. Calculation of the internal temperatures of cell 1 and cell 2 is shown in Equation (6), where $({\mathit{T}}_{\mathit{surfx},\mathit{L}}^{\mathit{k}}+{\mathit{T}}_{\mathit{surfx},\mathit{R}}^{\mathit{k}})/2$ represents the ambient temperature of the gap.

$$\begin{array}{c}{\mathit{T}}_1^{\mathit{k}+1}{=\mathit{T}}_1^{\mathit{k}}+\frac{\mathsf{\Delta }\mathit{t}}{{\mathit{m}}_{\mathit{cell},1}{\mathit{C}}_{\mathit{p},1}}{(\dot{\mathit{Q}}}_1-\left(\frac{1}{{\mathit{R}}_{\mathit{condx},1}}+\frac{{8\mathit{V}\lambda }_{\mathit{x},1}}{{\mathit{l}}_{\mathit{x}}^2}\right){(\mathit{T}}_1^{\mathit{k}}{-(\mathit{T}}_{\mathit{surfx},\mathit{L}1}^{\mathit{k}}{+\mathit{T}}_{\mathit{surfx},\mathit{R}1}^{\mathit{k}})/2)-(\frac{1}{{\mathit{R}}_{\mathit{condy},1}}+\frac{{8\mathit{V}\lambda }_{\mathit{y},1}}{{\mathit{l}}_{\mathit{y}}^2}{\left)\right(\mathit{T}}_1^{\mathit{k}}{-\mathit{T}}_{\mathit{surfy},1}^{\mathit{k}}\left)\right)\\ {\mathit{T}}_3^{\mathit{k}+1}{=\mathit{T}}_3^{\mathit{k}}+\frac{\mathsf{\Delta }\mathit{t}}{{\mathit{m}}_{\mathit{cell},3}{\mathit{C}}_{\mathit{p},3}}{(\dot{\mathit{Q}}}_3-\left(\frac{1}{{\mathit{R}}_{\mathit{condx},3}}+\frac{{8\mathit{V}\lambda }_{\mathit{x},3}}{{\mathit{l}}_{\mathit{x}}^2}\right){(\mathit{T}}_3^{\mathit{k}}{-(\mathit{T}}_{\mathit{surfx},\mathit{L}3}^{\mathit{k}}{+\mathit{T}}_{\mathit{surfx},\mathit{R}3}^{\mathit{k}})/2)-(\frac{1}{{\mathit{R}}_{\mathit{condy},3}}+\frac{{8\mathit{V}\lambda }_{\mathit{y},3}}{{\mathit{l}}_{\mathit{y}}^2}{\left)\right(\mathit{T}}_3^{\mathit{k}}{-\mathit{T}}_{\mathit{surfy},3}^{\mathit{k}}\left)\right)\end{array}$$

(6)

Obviously, cell 1 and cell 3, as the two outermost cells of the battery pack, have different heat transfer coefficients on the outer and inner end faces. The heat transfer coefficient of the inner end gap decreases with the thinness of the gap space, and h_d represents the heat transfer coefficient of the gap. Since the internal temperature fields of the cells are assumed to be independent of each other, the establishment of the thermal connection is based on the heat exchange between the shell temperature and the gap of each cell. Therefore, the difference in the group arrangement position of the different cells can be indirectly characterized by the difference in the heat transfer coefficient.

The assumptions of the thermal model are shown in

Table 2. Assumptions of the thermal model.

Assumptions	Explanations
The temperature conduction in z direction which is shown in Figure 1 is ignored.	The laminated structure inside the battery only exists in the x and y directions, but not in the z direction. The temperature difference in the z direction is ignored.
The temperature distribution inside the battery is uniform.	In order to reduce the complexity of battery modeling and simulation time, this work ignores the temperature difference at different positions inside the battery, referring to a lumped-parameter thermal model.
The resistance of wires in the battery pack is ignored.	Compared with the internal resistance of the battery, the wire resistance is too small and is ignored.

.

2.4. Characteristic Prediction under Constant Power Conditions

Equation (7) can be used to calculate the power of the battery pack, where P_total is the power of the battery pack:

$${\mathit{P}}_{\mathit{total}}{=\mathit{U}}_{\mathit{total}}{\mathit{I}}_{\mathit{total}}$$

(7)

A current iterative solution method was used to simulate battery characteristic performance when a cell is discharging at constant power. The processes of the prediction algorithm are shown in Figure 5.

As shown in Figure 5, the parameters were first initialized, including constant power, initial voltage, current variations, and allowable error value of power output; secondly, initial approximate current was determined by Equation (8), where I is the initial approximate current, P is the set constant power value, and U_total,0 is the initial terminal voltage; thirdly, the battery power and current were calculated repeatedly until the error was within the allowable range; finally, if the voltage of any cell exceeded the upper and lower cut-off voltage, the above processes ended.

$${\mathit{I}=\mathit{P}/\mathit{U}}_{\mathit{total},0}$$

(8)

3. Experiment

In this paper, MATLAB and Python were used for calculation and simulation. The data collection of terminal voltage, current, and temperature was achieved using equipment from Neware Co. Ltd., and the environmental temperature was controlled by a thermal chamber provided by Dongguan Bell Experiment Equipment Co., Ltd. During the testing processes, the above data were acquired synchronously. Note that the measurement errors of the equipment were 0.1 mV, 0.1 mA, and 0.1 K, respectively. A prismatic battery with a rated capacity of 42 Ah was selected as the research object. The specifications of this battery are shown in Figure 6 and

Table 3. Battery specification.

Specification	Parameter	Specification	Parameter
Height	10 ± 2 mm	Weight	about 0.721 kg
Width	25 ± 2 mm	Positive electrode material	NCM
Length	150 ± 2 mm	Negative electrode material	Graphite

.

Table 4. Battery parameters.

Parameters	Cell 1	Cell 2	Cell 3
y0, x0	0.4023, 0.7144	0.4073, 0.7549	0.403, 0.7196
Qp, Qnl	348200, 222740	320500, 212700	335900, 221000
τp, τn	10, 10	40,40	10,10
Pact	453910	334060	361500
Pconp, Pconn	219.1852, 115.7789	721.0612, 20.4525	219.7746115.3661
Rohm	0.0007	0.00015229	0.0007
τe	50	410.6422	100
hx, hy, hd	5, 30, 2	5, 30, 2	5, 30, 2
Cp	1160.17	1160.17	1160.17
λx, λy	2.35, 0.1526	2.35, 0.1526	2.35, 0.1526

and Figure 7 show the battery parameters. As shown in Figure 7, there is a gap between the two directions of R_cond, which is determined by the different properties of the two directions of the battery.

To verify simulation accuracy, a 20 A constant-current discharge experiment was conducted at a constant temperature of 293.15 K, and the connection mode and temperature acquisition points are shown in Figure 8.

The comparison of terminal voltage and temperature is shown in Figure 9. The mean absolute errors are shown in

Table 5. Mean absolute errors of temperature.

State Estimator	Mean Errors	State Estimator	Mean Errors
Total terminal voltage (V)	0.036	Point 4 temperature (K)	0.141
Point 1 temperature (K)	0.361	Point 5 temperature (K)	0.167
Point 2 temperature (K)	0.164	Point 6 temperature (K)	0.069
Point 3 temperature (K)	0.222	Point 7 temperature (K)	0.276

.

It can be concluded from Figure 9 that the temperature rise fluctuated noticeably in the middle period of discharge (about 3000 s–6000 s), and errors between the simulation and experiment are larger. During the whole discharge process, the temperature curves of cell 1 and cell 3 were slightly different, and cell 2 is obviously different from the others, which is because cell 1 and cell 3 were fresh batteries, and cell 2 was an aged battery used for many tests. Because these batteries had different historical usage times, different health conditions, and different capacities, for the same charging and discharging conditions the heat production of the battery was naturally different. In addition, due to their arrangement positions, the heat dissipation condition of cell 2 was not as good as the others—the temperature rise was larger. In the x direction, which is shown in Figure 4, cell 1 and cell 3 were symmetrically distributed in the arrangement structure, and the temperature rise was nearly the same. The temperature difference between shells at the middle and the end of the discharge (after 3000 s) is in contact with the gap which is shown in Figure 4. The results of comparing the accuracy of this model with other simulations are shown in

Table 6. Comparison results with other authors regarding temperature-simulation accuracy.

Others					This Work		Type
Reference	Type of Battery Pack	Working Condition	Ambient Temperature	Values	Working Condition	Values
[31]	7S4P	0.17 C (6.03 A) discharge	25 °C	0.99	0.5 C (20 A) discharge	0.9938	Maximum accuracy (Accuracy is 1-mean relative errors)
		0.4 C (14.444 A) charge	25 °C	0.98
		0.3 C (10.18 A) discharge	25 °C	0.96
		0.4 C (14.444 A) charge	25 °C	0.97
		0.4 C (14.44 A) discharge	25 °C	0.96
		0.4 C (14.44 A) charge	25 °C	0.95
[23]	5S1P	0.9 C discharge	23.3 °C	0.95
[32]	5S1P	WLTP class 3 drive cycle	5 °C	6.7%	0.5 C (20 A) discharge	1.35%	Maximum relative errors
			25 °C	1.5%
			45 °C	1.5%
[33]	1S1P	40 C (44 A) discharge	25 °C	6.22%	0.5 C (20 A) discharge
[34]	1S1P	1 C (2.2 A) discharge	About 12.5 °C	3.5%
		2 C (4.4 A) discharge
		3 C (6.6 A) discharge
		0.9 C discharge	23.3 °C

. It can be seen from Table 6 that the accuracy of the model simulations has been further developed.

Results of the comparison of measurements and predictions under constant power conditions of 120 W and 135 W are shown in Figure 10. Due to the intrinsic error of the experimental equipment, the power of the experimental equipment is greater than the expected value. The comparison of voltage and current at constant power discharge between the prediction and the experiment is shown in Figure 11. The predicted and simulated discharge times are shown in

Table 7. Time of constant power discharge.

Power (W)	Predicted Discharge Time (s)	Experimental Discharge Time (s)
135	12594	12621
120	14188	14168

.

As shown in Figure 11, when discharging at constant power, the voltage decreased, and then the current needed to increase to keep the power constant. Therefore, there was a short-term output of considerable current for the discharge terminal under constant power conditions. It can be concluded from Figure 11 that the current should increase by about 7 A. In addition, when the required power increased, the voltage decreased faster and the current increased faster compared to lower power conditions. With the increasing current, the temperature increased. The maximum temperature changes under several constant power conditions can be found in Figure 12 and specific maximum temperature rise values are shown in

Table 8. Specific temperature-rise values.

Constant Power (W)	Temperature Rise (K)	Constant Power (W)	Temperature Rise (K)
120	2.96	165	4.10
135	3.33	180	4.50
150	3.72	195	4.90

. It can be seen that the temperature rose rapidly at the end of discharge and the maximum temperature rise values increased with the increase in power. So, in order to ensure the safety of the batteries and prevent the temperature from getting out of control, the temperature rise needs to be controlled within a certain threshold.

4. Temperature-Control Strategies

The basic idea of a cooling method is to change the surface h and further reduce the battery temperature. Without discussing the specific cooling methods, this work developed a temperature-control strategy to keep battery temperature within a certain threshold on the basis of model prediction. According to the specific scale of the battery pack, the maximum h was considered to be 100 (W m⁻² K⁻¹) with reference to the air-cooling mode. Based on the above model, with the purpose of reducing the cooling cost or shortening the predicted consumption time, two temperature-control strategies were proposed to change the heat transfer coefficient (h) by iterative solution according to the predicted temperature. Two brief flow charts of these methods are shown in Figure 13, with the variable search noted by dotted lines.

A larger required h means more work needs to be performed by the cooling system and takes heat away from the battery, which will increase the cost. Thus, it is necessary to properly adjust the cooling system to obtain a suitable h. As for the strategy of using a constant search range for h, at each iterative calculation, h is searched in a constant large search range to obtain the most suitable h. As shown in Figure 13a, battery temperature is first predicted after a certain interval time of t_p, which is defined as the period of model-based temperature prediction. If the predicted temperature reaches the threshold, it is necessary to adjust the temperature-control means at this time, such as increasing the coolant flow rate for an air-cooled system or the wind speed for a liquid-cooled system. If the obtained h exceeds the upper limit assumed in advance, it is assumed that the system cannot meet the cooling requirements, and the batteries will be forced to stop discharging.

To reduce the calculation cost of the developed strategy and increase the system’s rapid-response capability, an optimal h is supposed to be selected at a variable search range. As shown in Figure 13b, the initial lower limit h in each iteration process is replaced by an h obtained in the previous process, which can reduce the search range, and the t_p of this strategy is thus smaller, which can reduce the time to obtain the maximum mean temperature. The maximum temperature of the battery pack is the internal temperature of cell 2, which can be seen for the simulation and the experiment in Figure 9. To sum up, Strategy 1 can obtain a more suitable heat transfer coefficient and cost less, and Strategy 2 needs shorter calculational time and has a steadier heat transfer coefficient.

The control temperature threshold in this work was selected as 296.15 K. The maximum temperature variation curves at different constant powers of 165 W, 180 W, and 195 W are shown in Figure 14a, Figure 15a and Figure 16a, and figures that can be used to prove that the battery was discharged at constant power during temperature-control are shown in Figure 14b, Figure 15b and Figure 16b, as well as the simulation time for the two strategies shown in

Table 9. Simulation time per cycle for two strategies.

Power	165 W	180 W	195 W
Simulation time for Strategy 1 (s)	0.505	0.597	0.718
Simulation time for Strategy 2 (s)	0.327	0.382	0.399

. It should be noted that the simulation time shown in Table 9 is the consumption time of each iterative computation for temperature control.

Taking the control strategy at constant 165 W power for example, variation results of the required h for temperature-control with the two control strategies are shown in Figure 17. The battery temperature was controlled with an increasing h when it reached the threshold at about 4500 s, which could be seen from the variations in battery temperature in Figure 14a. Different control effectiveness can be observed with these two control strategies: (1) the amplitude of h in Strategy 1 was smaller than that of Strategy 2, because Strategy 1 had a longer control period of t_p; (2) when the battery temperature was effectively controlled, the heat transfer coefficient of Strategy 1 decreased compared with an unchanged value in Strategy 2. In addition, when the temperature was controlled in a certain control period near the end of discharge in Strategy 1, the temperature of the battery soon reached the threshold again because of the obvious temperature rise. In order to achieve the temperature-control effect, Strategy 1 had to increase the heat transfer coefficient, thus resulting in a fluctuated h within a large range.

5. Conclusions

In this work, the heterogeneity of a battery pack of three cells connected in series was characterized using an electrochemical-thermal coupling model. The simulation accuracy of terminal voltage and temperature was verified at discharge as a 20 A constant-current load. Validations showed that the average absolute errors of battery terminal voltage and temperature were less than 36 mV and 0.4 K, respectively. Based on the model, an iterative solution method was used to predict the characteristics under constant-power conditions. Experimental verification showed that the power error was below 0.005%, and the discharge time deviation was within 30 s. Two temperature-control strategies were developed with the battery characteristic prediction method, and the effectiveness of the two methods was proven by simulation. A temperature-control strategy based on variable search range of heat transfer coefficient needs a shorter calculational time, while a temperature-control strategy based on constant search range of heat transfer coefficient costs less.