Analysis of the unsteady thermal response of a Li-ion battery pack to dynamic loads

It is becoming increasingly apparent that wide application of electric vehicles (EVs) are subject to significant improvements in battery technology. Temperature sensitivity is a major issue adversely affecting battery performance and requiring a robust thermal control. Yet, this is challenged by the large variety of temporal scenarios though which heat is generated in a battery pack, demanding dynamic tools to predict the thermal evolution of batteries. Classical transfer functions provide a low-cost and effective predictive tool. However, they are limited to linear systems, while nonlinear predictive tools can become impractical for EV applications. Therefore, this study provides a methodology to assess the dynamics of battery cooling. This is achieved through conduction of high fidelity modelling of battery cooling exposed to different temporal disturbances on the internal heat generation. The results are then post-processed to evaluate the extent of linearity. A quantitative measure of non-linearity is further applied to clearly determine the degree of nonlinearity in the heat transfer response. It is shown that battery cooling system can be approximated as a linear dynamical system as long as the disturbances are of short duration and relatively low amplitude. Conversely, long and large amplitude temporal disturbances can render strongly nonlinear thermal responses. © 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Electric Vehicles (EV) and Hybrid Electric Vehicles (HEV) have recently gained enormous attention due to rising environmental concerns [1]. Significant efforts are being made worldwide in support of electric vehicles to reduce greenhouse gas emissions [2]. However, a primary challenge facing EVs is to find energy storage solutions capable of high mileage, high performance driving and quick charging [3,4]. Recent advancements in energy storage technologies have further accelerated the battery-driven trend for automobiles. Lithium-ion (Li-ion) batteries have received significant attention for use within EVs because of their high energy density, lighter weight, low self-discharge rate, high recyclability and prolonged life compared to lead-acid or nickel-metal hydride batteries [5,6]. However, the performance e life cycle, discharge capacity e and safety of these batteries are heavily affected by the operational temperature [7e10]. Hence, effective control of the battery temperature has become an essential part of the overall battery management in EVs.
Studies by Bandhauer et al. [7] and Pesaran et al. [11] showed that the optimal temperature range for Li-ion batteries is between 15 C and 35 C. Further, the temperature difference within the battery module should not exceed 5 C [12]. Due to the limited temperature tolerance and their complex and relatively unstable chemistry, exceedingly low or high temperatures will cause irreversible damage to the battery and shorten its lifetime [13]. Importantly, this can become worse when hundreds or even thousands of battery cells are connected in series, parallel or seriesparallel [13,14]. Long-term existence of large temperature difference between battery cells can cause poor consistency within the battery module. This leads to cell by cell variation of the thermal effects, which in turn enlarges the temperature difference and creates a vicious cycle [13]. Further, degradation of a single cell results in poor battery module performance, i.e. the weakest cell may cause inadequate charging and discharging of the battery module. To prevent battery cell damage and maximise their performance, many battery thermal management systems (BTMS) have been investigated from external air and liquid cooling to passive cooling using phase change materials and heat pipes [15e17]. However, designing BTMS requires extensive understanding of their electrical, electrochemical and thermal behaviours. In particular, a deeper insight into heat generation, thermal transport and heat dissipation processes is required. Furthermore, due to the highly dynamic heat generation and heat transfer in Liion batteries, time-dependent thermal management is a crucial issue for them [18,19]. Over the last decade, most theoretical, numerical, and experimental studies on Li-ion batteries have been concerned with the dynamics of heat generation within batteries, see for example [19e22]. A small number of studies have also focused on the heat transfer aspect of battery cells under the dynamic loads set by unsteady charge and discharge of the batteries [19,20,23,24]. However, currently, there is no study on understanding the dynamic response of heat transfer in battery cells with time-varying thermal loads.
Yang et al. [23] used a 2D conjugate heat transfer model coupled with a 1D electrochemical model to numerically study the effect of aligned and staggered battery cell arrangements using forced convection. The study also focused on the effect of longitudinal and transverse spacing of battery cells. A k-ε turbulence model was employed to accurately capture the flow field and the corresponding heat transfer characteristics. The staggered battery alignment lead to cooler overall battery module temperatures, however, the aligned arrangement was more effective for reaching temperature uniformity between battery cells. Furthermore, increasing the transverse gap between battery cells lead to the aligned arrangement dominating in both energy and heat transfer efficiency. Nevertheless, this study revealed that for a densely packed battery module, the transverse and longitudinal gaps between cells should be kept minimum. Wang et al. [19] numerically and experimentally studied the thermal management of Li-ion battery modules under dynamic loads. They provided a detailed computational fluid dynamics (CFD) analysis as well as a reducedorder model (ROM) requiring far less computational power. The experimental model was setup through a fixed fan and four cylindrical cells excited by a dynamic load profile. The temperature profiles of the battery cells were setup using a feedback loop to a fan controller to either reduce or turn off airflow in order to avoid over cooling. This was to keep the battery cell temperatures within the optimal operating range while minimising parasitic power consumption. Wang et al. [19] also found that the use of a control algorithm to actively cool the battery cells reduces the parasitic power consumption by 30% compared with that in continuous cooling.
Li et al. [25] studied the battery characteristics of an LiFePO4/ Graphene (LFP/G) hybrid cathode lithium-ion battery using an electrochemical-thermal coupled model. Moreover, the study also focused on the total heat generation of the battery under varying discharge rates. The analysis of the chemical, thermal and electrical processes was done using an MSMD (Multi-Scale Multi-Dimensional) in Ansys Fluent. The study was then validated using an experimental approach focusing on the current density, positive and negative current potential and thermal distribution within the battery. The validation allowed for an improved Ansys Fluent MSMD electrochemical-thermal model coupled with mass, energy, charge conservation, and electrochemical kinetics. Panchal et al. [26] investigated methods for cooling lithium-ion batteries by utilising cooling plates to cool down a large prismatic LFP/G battery. Using StarCCM þ for numerical modelling and with the addition of laboratory testing, Panchal et al. studied the thermal and flow characteristics of the coolant inside the cooling plate at high discharge rates of 3C and 4C. It was found that as the discharge rate of the battery is increased, the temperature distribution within the cooling plate and the heat flux measured at the anode, cathode and the centre of the surface also increased. These results can be considered for designing more efficient cooling techniques for batteries.
Chu et al. [27] modelled the use of a parallel-series combined liquid cooling system for a 288 V Nickel-metal hydride (Ni-MH) battery pack. Using StarCCM þ for simulating and analysing the heat transfer and fluid dynamics of the cooling system, the study indicated that the cell's temperatures and temperature differences could be kept within an ideal range. The study was validated by a battery pack temperature experiment on a bench and in a vehicle. The results found that by increasing the cross-section ratio of the coolant, the temperature difference across the battery pack decreased by 2.5 C. Moreover, increasing the number of inner supporting walls, can further aid in reducing the maximum temperature by increasing the overall heat transfer; however, this will lead to a greater pressure drop across the pipeline. Increasing the number of supporting walls from 0 to 9 causes the pressure drop to increase by 2.5 kPa while the maximum temperature is decreased by 3 C. Ling et al. [20] investigated hybrid thermal management system using phase change materials (PCM) to ensure each battery cell remains within the operating range and the cell-to-cell temperature difference is minimised. The PCMs were cooled using forced air convection to avoid reaching a critical melting temperature that leads to the failure of the thermal management system. Ling et al. [20] further studied the effects of varying air speed and the thermophysical properties of the PCMs. The use of a passive thermal management system combined with forced air convection provided an efficient means to minimise cell to cell temperature difference while the battery pack was kept at an acceptable temperature range. Mahamud and Park [28] studied the effect of reciprocating airflow on improving the thermal uniformity of Li-ion battery. They used a single fan and flip doors to alternate the flow direction creating a reciprocating flow. These authors concluded that by introducing a reciprocating airflow, the cell temperature difference reduced by 4 C (72% compared to non-reciprocating flow) and the maximum temperature dropped by 1.5 C.
Despite having small heat capacity and low thermal conductivity, air is one of the most common fluid used in BTMS due to its availability and low costs [1]. Normally, air is propelled using fans and driven over the surface of battery cells. This might be acceptable for a small number of battery cells, however, for more densely packed battery cells, the coolant temperature may significantly rise after flow over the first few cells leading to higher temperatures at the outlet of the battery module compared to inlet [19]. As already pointed out, this can severely impact the lifetime and safety of the battery pack. The issue can be addressed by using hybrid thermal management systems, like those discussed by Ling et al. [20], or by increasing the Reynolds number of the coolant flow. Since most secondary systems on an EV are also powered by the battery module, increasing the flow rate would increase parasitic power, reducing the overall range of the EV. For high density battery modules or high cell discharge rates, air is inadequate [15,29]. Liquid coolants (e.g. water and water mixtures) offer many advantages over air, as such liquid-based BTMS can be 3500 times more effective than air and reduce the use of parasitic energy [29,30].
Taking an experimental approach, Zhu et al. [31] studied the effect of immersion cooling on the performance of high concentration photovoltaic (PV) cells. The study utilised de-ionised water as the coolant to ensure the PV cells would remain operational at a low working temperature. Using immersion cooling, the module temperatures peaked at 49 C with the cell to cell temperature distribution being less than 4 C. Although the cooling capacity of liquid immersion cooling is very favourable, the technology is still immature for EV application. Pesaran [15], Chen et al. [30] and Karimi and Li [32] investigated the effects of liquid cooling on batteries. The higher thermal conductivity of liquids leads to higher rates of heat transfer. Yet, due to the higher viscosity of liquids compared to air, liquid cooling often encounters high pressure losses and hence high parasitic power use. As a result, liquid coolants are normally kept at low flow rates. Karimi and Li [32] numerically compared the natural convection and force air and liquid convection to develop an optimised thermal management system for batteries. Forced convection cooling allows the battery pack temperature to remain within the optimal operating range. However, it leads to unbalanced temperature and voltage distributions in the pack which deteriorate the batteries durability. There exist few studies on the effects nanofluids on the thermal behaviour of battery cells and battery modules [33e36]. For example, Wiriyasart et al. [33] studied the temperature distribution and pressure drops of using Titanium Dioxide nanofluids in corrugated mini-channels of the battery module using a computational analysis approach. The study found that in comparison to water, the nanofluid coolant leads to a reduction of 27.6% in the maximum battery surface temperature. Furthermore, the concentration of nanoparticles had a significant effect on the cooling capacity of coolant, resulting in further lowering of the maximum battery temperatures.
The preceding review of literature indicated that the existing BTMSs are primarily concerned with accurately evaluating the temperature field of the cooling fluid and heat generation inside battery cells. Unsurprisingly, this requires the use of very demanding computations to precisely predict the thermal behaviour of battery cells. However, due to the highly dynamic nature of EVs, high order computational models are difficult to be used in BTMS [19]. Hence, there is a pressing need for low-cost predictive models capable of quickly and accurately simulating the dynamic response of heat transfer within the battery module. Nonetheless, so far, most studies on battery cooling have been concerned with steady cases and little attention has been paid to development of dynamic models. The current work aims to fill this gap by using a computational model of unsteady heat transfer in a battery pack. It, then, determines the applicability of the linear predictive tools by evaluating the non-linearity of the dynamic response of heat convection in the battery pack.
The rest of this paper is structured as follows. The numerical and theoretical methods used for simulations and the subsequent postprocessing of the results are described in Section 2. The results are presented and discussed in Section 3 and Section 4 summarises the key findings of the study.

Problem configuration and assumptions
Sketches of a battery pack, the simulated model and a twodimensional schematic of a single cell section are shown in Fig. 1a, b and 1c, respectively. The model has been configured in a staggered cell configuration and has six primary cells. This number of cells was chosen based on the previous studies in this area [20,23,37]. Simulating a full battery module is a computationally expensive process, therefore, only a section of it was simulated, as shown in Fig. 1a, outlined in green. The battery cells have a radius R of 20 mm [19,20,28], a single cell section, shown in Fig. 1c, has a height of 50 mm, length of 100 mm and a depth of 100 mm [28]. The following assumptions are made throughout the proceeding analyses.
The flow is thermally and hydrodynamically fully developed. The flow is fully turbulent with a steady inlet velocity. For air, the fluid is an ideal gas, whereas, for water and nanofluids, it is simulated using a constantly density model. Gravitational effects and heat generation due to viscous resistance are ignored. Nanofluids are modelled as a single-phase liquid.
Further, the unsteady thermal heat flux was applied to the external surface of each battery cell once the model reached steady state conditions. In justifying this, it is noted that any heat generated in the battery can only be removed from the system once it has reached the surface of the battery cell. Therefore, the heat loss from each battery cell could be modelled by a surface heat flux [38e42].
In the following subsections the employed mathematical and computational models are explained. The numerical scheme is then validated rigorously against a range of existing experimental and numerical data.

Governing equations, boundary conditions and numerical model
Taking a Reynolds-Averaged Navier-Stokes approach to the modelling of turbulent flow, the continuity equation is expressed by Conservation of momentum for the fluid flow in the x, y and z direction are given by v vt where d ij is the Kronecker delta function defined by The conservation of energy is written as The boundary conditions applied to the simulated model include a no-slip boundary on the external surface of the battery cell walls. Further, the top, bottom, front, and back walls of the model have symmetry condition applied since the battery cells can be packed in the y and z-direction. Finally, the ambient and fluid inlet temperature were assumed to be 300 K with an outlet pressure of 0.1 MPa.
Reynolds number based on the battery cell diameter is defined as Nusselt number is expressed on the basis of the following equations.
Under steady conditions, the thermal boundary around each battery cell was modelled using a constant heat flux of 6500 W/m 2 . This value was set based on the current model geometry, 85 kW battery specification of the Tesla Model S [43] and an average battery cell efficiency of 85% [44]. It is worth noting that the expected air velocity of a vehicle being driven on outer city roads is 30 m/s (Re ¼ 75,000), which was also the chosen velocity for the air model. Due to the poor thermal conductivity and capacity of air, higher flow velocities are required to adequately cool down battery cells for optimal operation. Furthermore, the Reynolds number for water and nanofluids is 2300. These values of Re are kept constant throughout the study. Once the model reached steady-state conditions, the heat flux was modulated with time by superimposing a sinusoidal wave with varying amplitudes and frequencies The parameter a determines the amplitude of disturbance and is regarded as the amplitude of modulation/excitation in the proceeding sections. Further, the frequency is non-dimensionalised using the Strouhal number relation given in Eq. (8) and then normalised using Eq. (9).
The changes in dynamic viscosity and thermal conductivity of the Al 2 O 3 -Water nanofluid were calculated using relations given in Eqs. (10) and (11). These well-established models were developed by Brinkman and Maxwell-Garnetts, respectively [45,46].
where the concentration of nanoparticles is given by 4. Furthermore, during the nanofluid studies, 2.5% and 5.0% nanoparticle concentrations are used. The material properties for the coolants used in the current study can be seen in Table 1 with all the thermophysical properties evaluated at a temperature of 300 K. A parametric study was subsequently conducted where the coolant fluid, and the amplitude and frequency of the modulation were varied systematically. It should be noted that the heat flux on the surface of batteries was modulated at low frequencies (f 2HzÞ. This is because thermofluid systems are often most responsive to low forcing frequencies only [47,48]. The numerical simulations were conducted using the finitevolume based Computational Fluid Dynamics software, StarCCM þ v14.02-R8. An unsteady, three-dimensional, turbulent flow solver coupled with the energy equation was implemented within the model. A turbulent k-ε turbulence model was employed to accurately capture the flow field of the cooling fluid between the cells. Further, depending on the type of coolant, different fluid models were used. For air, the ideal gas model was employed, whereas, for water and nanofluids, the constant density model was utilised. An implicit unsteady solver was used to precisely model the fluid behaviour. The time-step was set to be three order of magnitude smaller than the physical time scale to accurately model the vortex shedding. Additionally, all the model utilised a secondorder discretisation scheme to further enhance their accuracy. The computational model was run in an HPC using a single compute node of an Intel Xeon 4830 (24-cores) due to the highmesh density and the complexity of the model. The simulated data were then exported to MATLAB 2019b for post processing.

Calculation of transfer functions
During this study, the inbuilt Fast Fourier Transform (FFT) functions in MATLAB were used to analyse the system response in the frequency domain. Prior to the application of FFT, the Nusselt number was normalised according to the following relation.
where NuðtÞ denotes Nusselt number at the current time and Nu is the time averaged Nusselt number. Furthermore, transfer functions were introduced to predict the dynamic response of batteries to perturbations in surface heat flux. Of course, this implies that the thermal system under investigation features linear dynamics. The amplitude of the transfer functions is calculated using aðf Þ ¼ jNormðNuðtÞÞ i j; (13) in which i specifies the battery cell. Also, cross correlation was employed to measure the time difference between the input and output signal, which is then non-dimensionalised: where ndd is the nondimensionalised delay, td is the crosscorrelation time delay, dl is distance from the inlet to the battery cell and i donates the battery cell. Moreover, the phase of the transfer function, measured in radians, was calculated using the time difference, td. That is However, during this study, the non-dimensionalised delay was used to show the phase of transfer functions.

Grid independency and validation
A polyhedral staggered mesh with prism layers around each battery cells was used to capture the heat transfer while minimising computational costs (see Fig. 1b). Multiple tests were carried out with varying mesh cell size to determine at which grid density the highest accuracy and lowest computational burden could be achieved. The cell size refers to the relative size control of all other values such as surface cell size, volume cell size and prism layer thickness. The grid independency was verified when the Nusselt number (Nu) of each battery cell fell within a variation band of one percent. These tests were performed with a Reynolds number of 5100. The outcomes are shown in Table 2. As the cell size is decreased, the mesh density increases, and the solution converges to within a margin of one percent. To balance the model accuracy and computational cost, a cell size of 0.0025 m was used throughout the current study.
The simulations were validated by comparing against the existing experimental and numerical data. First, the model parameters were changed so that the flow occurs over a single cylinder rather than a bundle of cylinders. Tables 3e5 depict the comparison between the simulated and experimental data. Table 3 shows that at low values of Re, the simulated results are in excellent agreement with the values obtained using the Churchill and Bernstein [49] empirical correlation. The error is higher at larger values of Reynolds number. Nonetheless, it is noted that it is still well within the error limits of these correlations which can reach over 20% [49]. Further, at the Reynolds number of 3900, the simulated Coefficient of Drag, C d , was compared with the experimental data, revealing an error of 1.4%.
From Table 4, the simulated Nu has a maximum error of 12.79% at a high value of Re and a minimum error of 1.91% at a low Re. However, comparison of the simulation results against the Churchill and Bernstein [49] correlation results in an error of less than 10% for all values of Re. Further, the simulations were repeated for different Re for comparison against experimental data. As seen in Table 5, for all values of Re, the error is less than five percent. Next, air flow across a bundle of tubes was compared against the empirical correlations by Grimson [51]. At a velocity of 30 m/s, the simulated model shown in Fig. 1b produced an average Nusselt number of 215.056, whereas, the correlation presented by Grimson produced an average Nusselt number of 233.05267, leading to an error of 8.4%.
The good agreement between the simulations, empirical correlations, and experimental data for both air and water confirms validity of the numerical analysis. Fig. 2 shows a flow chart illustrating the post processing steps. On completion of the numerical simulations, all data is exported to text files which can then be easily imported for post processing. During post processing, transfer function calculations and measures of non-linearity using Fourier Transforms and phase portraits are conducted. Further, if the latter is found to be less than 10%, transfer function analysis is conducted. These will be further discussed in the following section. To ensure about robustness of the analysis, all the developed post processing tools were validated against synthetic data.

Results and discussion
In this section the dynamic response of the battery cells to sinusoidal disturbances imposed on heat flux (see the details in Section 2) are investigated. It is noted that any arbitrary temporal disturbance can be decomposed into a series of sinusoidal components through Fourier transformation. Assuming linearity, the dynamic response of the system can be readily found by adding the system response to each of the sinusoids [54]. Hence, a knowledge of the system response to sinusoidal disturbances is essential for predicting the system response to all temporal fluctuations imposed on the battery heat release. For this to hold, the system should remain linear and therefore it is important to examine deviation of the dynamics of thermal system from a linear behaviour. The current section focuses on addressing these two issues. Fig. 3 shows the steady-state values of the heat convection coefficient, h, values for different coolants used in the current study.   The poor thermal conductivity of air leads to a low value of h as expected, whereas, water flow, with a significantly higher thermal conductivity, offers considerably higher convection coefficients. As foreseen from previous studies, the addition of nanoparticles leads to further improvements in the fluids heat transfer capabilities, leading to higher heat transfer coefficient values. This is further enhanced by increasing the concentration of nanoparticles in the fluid.
The spatiotemporal response of the fluid temperature field to a sinusoidal heat-flux disturbance imposed on the surface of each battery cell is shown in Fig. 4. This figure corresponds to a case with 30% amplitude of the fluctuations in the surface heat flux at forcing frequency of 1 Hz. Fig. 4 depicts a typical convective system in which the fluid temperature approaches that of the battery cell surface as it flows downstream towards the outlet. The temperature across rises when the sinusoid reaches its peak value. This is further illustrated by the cooling region near the inlet retreating and the downstream temperature increasing. As the sinusoid continues and reaches 180 , the cool region near the inlet advances, however, it does not equally represent the temperature field seen at 0 . This is due to the modulation frequency being 1 Hz, which means for a quarter segment of the sinusoid to be processed, only 0.25s are taken; an insufficient amount of time for the model to return to its previous state. However, at the sinusoids final segment, the cool regions near the inlet further advance as this is the trough of the sinusoidal disturbance. Further, this decrease in temperature across the model allows the temperature field to reset as the sinusoid returns to its starting position. Fig. 5 depicts the temporal response of the normalised Nusselt number and the corresponding response in the frequency domain, with air as the coolant and an amplitude of 10% and modulation frequency of 0.25 Hz. Fig. 5a shows that the Nusselt number response of the three examined battery cells is very closely represented by a sine wave, while there is an increase in the amplitude for those battery cells located farther from the inlet. This is due to an increase in the fluid temperature as it travels downstream requiring the need for a larger heat transfer coefficient in order to continue extracting the heat flux produced by the cells. It causes a larger variation in Nusselt number compared to cells located closer to the inlet, producing larger amplitudes. This behaviour is further Fig. 4. Spatiotemporal evolution of the temperature field exposed to a temporal sinusoidal disturbance imposed on the surface heat flux. Fluid Type : Air; a ¼ 30%; f ¼ 1:00Hz. confirmed by the spectrum of the time trace of Nusselt number shown in Fig. 5b. A classical sign of a linear system is the equality of the frequency of the response and that of excitation [54,55]. Although the spectrum in Fig. 5b contains two spikes, the second spike is comparatively insignificant; therefore, in this case, the heat transfer in the battery pack behaves as a linear system. Later, a more rigours method of assessing the system linearity will be put forward. Fig. 6 shows the temporal and spectral response of Nusselt number of different battery cells when water is the coolant and a disturbance with the amplitude of 10% and frequency of 0.25 Hz is imposed. The temporal responses shown in Figs. 5a and 6a differ significantly due to the large difference in the thermophysical properties of the two fluids. A general trend observed in Figs. 5 and 6 is the growth of Nusselt number amplitude for those cells that are located farther downstream of the coolant entry. Fig. 4 showed that the average coolant temperature increases as the outlet is approached, while the cells are set to lose a fixed value of heat flux. It follows that the downstream cells surrounded by hotter fluid should feature larger amplitudes of Nusselt number fluctuations. The extent of this behaviour is more significant in Fig. 5 compared to that in Fig. 6. This could be explained by noting the large  difference in the specific heat capacity of water and air, which results in much larger temperature increases in air in comparison with water. Unlike temporal responses in Figs. 5 and 6, the spectral responses are quite similar with a single large peak at the excitation frequency and a smaller, insignificant peak at the second harmonic of the excitation frequency. As with the previous case of air, the water case can be also approximated as a linear dynamic system. Further, the sign of another linear dynamic system can be seen in Fig. 7, where the coolant fluid was changed to a Al 2 O 3 -Water nanofluid with a nanoparticle concentration of 2.5%. Also, further increasing the nanoparticle concentration to 5.0% leads to the results similar to those shown in Fig. 7. Due to the similarity between the results shown for the cases of water and nanofluids, no further nanofluid results are presented.
A transfer function, providing information on the amplitude and phase of the dynamic response, can be used to predict the dynamics of linear systems [54]. The classical concept of a transfer function is regularly used for the linear systems with a single input and single output (SISO). In the current study, we consider the surface heat flux as the input, and the oscillating, normalised Nusselt number for each battery cell as the system output. Accordingly, six transfer functions are defined e one for each battery cell. The amplitude and the non-dimensional delay of the transfer functions were calculated for all battery cells and are presented in Figs. 8 and 9, respectively. Fig. 8a clearly shows that as the excitation frequency increases, the amplitude of the transfer function decreases, showing that the system is more responsive to lower excitation frequencies. Previous studies in which the inlet fluid was excited [48,56,57] rather than the thermal load on the system have also featured a stronger response at low frequencies. In general, excitation at lower frequencies tends to give the system larger time to respond to the  disturbances, leading to larger amplitudes. Further, as the thermal load disturbance frequency is increased, the non-dimensional delay decreases, which is in keeping with the findings of previous studies [56,57]. This is due to the available time between the disturbance and heat convection becoming significantly shorter with increasing frequency, causing the non-dimensional delay to diminish. As expected, the amplitude response of the cells located further from the inlet is larger than those located near the inlet. This can be explained by noting that the bulk temperature of the fluid flow increases as it travels downstream. Thus, as the fluid passes each battery cell, the temperature difference between the fluid and the subsequent battery cell decreases which, requires growth in the magnitude of convection coefficient and therefore larger amplitude of normalised Nusselt number. This is consistent with the behaviours shown in Figs. 5a, 6a and 7a in which the amplitudes of the normalised Nusselt number are significantly larger for battery cells near the outlet. The transfer function amplitude and the non-dimensional delay for the six battery cells with water as the coolant fluid are shown in Fig. 9a and b, respectively. From visual observation, Figs. 8 and 9 feature the same trend of a decrease in amplitude and nondimensional delay as the frequency increases. However, for the water case, the transfer function amplitude is far larger compared to the case of air. Further, at high frequencies, the non-dimensional delay is non-existent for water. This could be explained as follows. As shown in Fig. 3, the steady convection coefficients for water flows are significantly larger than those of air. Hence, under steady condition the surface temperature of the battery is smaller when air is used as the coolant. This implies that for a given disturbance in the heat flux, the water-cooled system needs to experience a larger fluctuation in the convection coefficient and Nusselt number, resulting in larger amplitude of the transfer function.
The foregoing transfer function findings were primarily based on the assumption of a linear system confirmed in Figs. 5e7, where an insignificant second peak is observed in the spectrum of the Nusselt number. However, by increasing the amplitude of modulation, some cases become either mildly non-linear e a small second peak is observed in the spectrume or strongly non-linear e where multiple large peaks appear. For example, a slightly nonlinear case can be seen in Fig. 10, where the amplitude of modulation was increased from 10% to 30% for air. It is clear that the normalised Nusselt number graph for battery 1 no longer represents a sine wave and the graph for battery 3 and 6 are starting to deviate from a regular sine wave, providing a clear indication of non-linearity. This inspection can be confirmed by the spectral response of the system seen in Fig. 10b, where a second peak and an insignificant third peak can be observed. Increasing the amplitude of modulation to 30% for the water and nanofluid cases led to the results like those shown in Figs. 6 and 7, and thus for brevity they are not shown here.
Increasing the amplitude of excitation to 60% with the coolant fluid being water leads to the appearance of signs of non-linearity. The time trace and spectrum of the Nusselt number for water with a modulation amplitude of 60% can be seen in Fig. 11a and b, respectively. From visual inspection, the response of normalised Nusselt number may seem to represent a regular sine wave, however, the FFT response shows two peaks e one large peak at the excitation frequency and another at double the excitation frequency. Further, a completely non-linear case can be seen in Fig. 12, where the coolant fluid is air and the amplitude of excitation is 60%. As seen in Fig. 12a, the graph for the normalised Nusselt number on cell 1 has completely deviated from a sinewave along with the graph for cell 3 and cell 6. This is further supported by the spectral analysis, where multiple large peaks are observed. For the case with air as the coolant with 30% amplitude modulation, three peaks occur. Yet, for the same coolant fluid, further increasing the amplitude to 60% leads to five peaks occurring in the spectral response of the system, with the largest peak being at the excitation frequency. This is a clear indication that the system has become strongly non-linear.
This finding is further supported by Figs. 13 and 14, which show the Lissajous patterns (phase portraits) of the normalised Nusselt number vs the normalised input signal calculated for the cases of air and water as the coolant fluids, respectively. Phase portraits have been used in the previous studies to further analyse the linearity of a dynamic system [48,58]. In Figs. 13a and 14a the oval shapes of the phase portrait are axisymmetric, corresponding to a linear case. In Fig. 14b, where the coolant fluid is water and the modulation amplitude is 30%, a slight divergence from axisymmetricity can be seen. However, this case can be still considered linear as the deviation is insignificant. This feature can no longer be seen in Fig. 13b, c and 14c indicating a clear departure from the linear response with Fig. 13c corresponding to a strongly non-linear system, for which the phase portrait becomes highly asymmetric. Unlike the linear cases shown in Figs. 5e7, prediction of the dynamics of non-linear systems are cumbersome and cannot be predicted using transfer functions. In many instances, the dynamics can be predicted only through high-order modelling, which is computationally expensive. Therefore, it is essential to identify conditions where the system dynamics are non-linear. This can be achieved by analysing the measure of non-linearity. The measure of non-linearity is used to evaluate the divergence of a system from linearity. This is given by Refs. [59,60].     n e n e À F d (16) where d nL is the measure of non-linearity, n e and F d are the Euclidean norm and the discrete Fourier transform spectrum of the normalised Nusselt number for each battery cell, respectively. For a truly linear system, Eq. (16) collapses to a value of 0 and for a completely non-linear system, a value of 1 is achieved. Therefore, any real system will have a value between 0 and 1. Fig. 15 shows the maximum measure of non-linearity (d nL ) at varying modulation amplitudes and frequencies for two different coolants, air and water. The figures for nanofluids at 2.5% and 5.0% concentration of nanoparticles are like those shown for water in Fig. 15b with minor differences and, therefore, are not shown. Fig. 15a shows that as the amplitude of modulation increases from 10% to 30%, there is a drastic increase in the value of measure of non-linearity. This further intensifies as the amplitude of modulation approaches 60% confirming the presence of a strongly nonlinear dynamic system at this amplitude. Doubling the excitation frequency from 0.25 to 0.5 Hz quickly decreased the measure of non-linearity regardless of the amplitude of modulation. This decrease continues as the frequency is further increased, however, the sharp decrease observed during the initial change in frequency is no longer present, rather a small gradual decrease can be seen. A similar trend but at a smaller scale is also present in Fig. 15b. Therefore, as expected, the system dynamics approach linearity at low amplitudes. In addition, linear behaviour is more likely to exist at higher forcing frequencies. Importantly, the tendency of air to feature nonlinear behaviour at high modulation amplitudes and low forcing frequencies is considerably higher than that of water and nanofluids.
The origin of non-linearity in the present problem is the interactions of the cooling fluid with the boundary layers and wake region around and behind each battery cell. It is well-established that such interactions are rather complex. The current results shown in Fig. 15 indicate that the long duration (low frequencies) modulation have enough time to sufficiently interact with the system generating a non-linear response. Conversely, short-term modulations (high frequencies) do not find the time required for the fluid dynamic interactions and hence present linear dynamics at a low modulation amplitude. This implies that in practice, long temporal disturbances with large magnitudes render nonlinear dynamics. However, short-term disturbances particularly with smaller magnitudes cause linear responses. Thus, a transfer function approach can be utilised to predict the dynamics of heat transfer as a computationally inexpensive predictive tool. However, at low frequencies and high amplitudes, the system dynamics cannot be interpreted using a simple transfer function, rather extensive high-order modelling for accurate prediction is required. Further, the significant difference between the linearity of air and water-cooled systems is likely to be related to their different thermophysical properties and particularly Prandtl number. Lower Prandtl number of air in comparison to water hinders convective heat transfer at fixed Reynolds number. This renders a stronger nonlinear response of the Nusselt number at any given amplitude of modulation. In closing, it is recalled that the current analysis neglects electrochemistry of the battery. Although it is speculated that the results are not considerably influenced by such assumption, this issue should be noted in any generalisation of the current results.

Conclusions
A major issue facing EVs is their battery performance, which is heavily impacted by the operational temperature. Batteries in EVs are exposed to a wide range of temporal scenarios in heat generation and transfer due to the existence of numerous vehicle manoeuvres and driving patterns. Therefore, thermal control of batteries requires prediction of their temperature. Such predictions are conventional for linear systems through utilisation of the classical transfer function approach. However, this becomes significantly more complicated for systems featuring nonlinear dynamics. Hence, it is important to determine under which conditions the battery cooling system features linear and nonlinear dynamic response to disturbances imposed on the battery heat generation. To achieve this, principle of Fourier transform was utilised noting that any temporal disturbance in battery heat release can be decomposed into a series of sinusoids with varying frequencies. The system response was, therefore, examined by numerical simulation of the response of the average Nusselt number on each battery cell to sinusoidal fluctuations in the heat flux. The results were then analysed in time and frequency domains. Further, a rigorous measure of nonlinearity was applied to the computational data to determine the regions at which the system response can be assumed linear. The key findings of this study are summarised as follows.
Dynamic response of heat convection from the battery cell remains nearly linear at low modulation amplitudes (10%). Under linear condition, a transfer function was developed for average Nusselt number which feature the characteristics of a conventional low-pass filter. At higher modulation amplitudes average Nusselt number features mildly to strongly nonlinear behaviour. Low frequency (long duration) temporal disturbances of the battery heat flux tend to generate strong nonlinear response. This can be attributed to the long-time available for completion of fluid dynamic and transport interactions. The classical transfer function approach can be utilised to predict the dynamics of battery cooling as long as the disturbance magnitude is relatively low (around 10% of the mean value) and its duration is short.
It remains as a future task to develop low-order models for predicting the dynamics of heat transfer with regards to non-linear systems. Implementation of the findings of this work into the battery thermal management system remains as a future task.

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