A new approach to the internal thermal management of cylindrical battery cells for automotive applications

Conventional cooling approaches that target either a singular tab or outer surface of common format cylindrical lithium-ion battery cells suffer from a high cell thermal resistance. Under an aggressive duty cycle, this resistance can result in the formation of large in-cell temperature gradients and high hot spot temperatures, which are known to accelerate ageing and further reduce performance. In this paper, a novel approach to internal thermal management of cylindrical battery cells to lower the thermal resistance for heat transport through the inside of the cell is investigated. The effectiveness of the proposed method is analysed for two common cylindrical formats when subject to highly aggressive electrical loading conditions representative of a high performance electric vehicle (EV) and hybrid electric vehicle (HEV). A mathematical model that captures the dominant thermal properties of the cylindrical cell is created and validated using experimental data. Results from the extensive simulation study indicate that the internal cooling strategy can reduce the cell thermal resistance by up to 67.8 ± 1.4% relative to single tab cooling, and can emulate the performance of a more complex pack-level double tab cooling approach whilst targeting cooling at a single tab. © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

A new approach to the internal thermal management of cylindrical battery cells for automotive applications 1

. Introduction
With their high energy and power density, lithium-ion batteries possess attractive characteristics for energy storage systems integrated within future hybrid (HEV), plug-in hybrid (PHEV) and full electric vehicles (EVs). Combined with the significant recent and forecast declines in the cost of this technology [1], hybridised and all electric automotive vehicles employing lithium-ion cells are witnessing increased market penetration rates. A recent 2015 report by KPMG [2] highlights the potential for electrified vehicles to be between 11 and 15% of new vehicle sales within the EU and China by 2025. Within the US, the market may comprise 16e20% of vehicles over the next 10 years. These predictions are comparable to those cited in Ref. [3]. The article collates a number of studies and concludes that by 2025, there will be in excess of 11 million EV sales worldwide, with approximately 6 million in North America (20% of new vehicle sales). Other reports [4,5] also forecast a substantial increases in market share for electric variants of vehicles within the commercial light duty vehicle market (LDV) over the coming decade.
Stringent requirements for automotive vehicles, as discussed within Ref. [6], include safety (i.e. thermal runaway avoidance for a battery based energy storage system), immediate performance and durability (i.e. drive range). These requirements present challenges for lithium-ion batteries as their performance and ageing rate is particularly sensitive to the thermal condition which must be controlled [7]. Low temperatures of circa < 15 C have been shown [8,9], to reduce the immediate available cell capacity and negatively impact power output from reduced mass transport rates. Higher temperatures, however, accelerate cell ageing due to the increased rates of unwanted parasitic side reactions that exist within the cell, mainly from the growing solid-electrolyte interface (SEI) layer which incurs a loss of cyclable lithium and increased capacity fade [10]. As such, to maximise performance whilst limiting ageing, a number of studies agree that the absolute operating temperature of lithium-ion cells should be kept within a range of circa 15e25 C regardless of the ambient temperature [8]. Temperature gradients that may form between cells and within the internal layers of individual cells is also of great concern. Yang et al. [11] demonstrated that temperature gradients between cells connected in a parallel string can exacerbate unbalanced discharging, with the ageing rate of the battery increasing linearly as the cell-to-cell temperature gradient increases. On the cell level, Fleckenstein et al. [12] concluded through simulation that temperature gradients throughout the internals of individual cylindrical cells cause inhomogeneities in cell current density, which in turn induces a local state of charge (SOC) imbalance within the cell. Troxler et al. [13] further investigated the effects of in-cell temperature gradients and noticed that the cell performance under a temperature gradient did not perform as if the cell was operating at the volume average temperature, but rather as if at a higher average temperature than the theoretical volume average. Owing to this, previous reports within the literature suggest that temperature gradients between lithium-ion cells in a pack and through the individual cell material should not exceed a maximum of 5 C [11,14e17]. Recently, however, Hunt et al. [18] demonstrated that not only do temperature gradients negatively affect cell performance and ageing, but the nature of the gradient can further exacerbate the effect and must be considered carefully by energy storage systems engineers tasked with aggregating individual components into complete battery modules or packs. Specifically, gradients perpendicular to the layers within the cell induced by surface cooling were found to accelerate the overall cell degradation rate compared to when the gradient is along each of the layers (i.e. in-plane) achieved from tab cooling.

Mathematical notation -Symbols and units
Careful consideration needs to be given to the design of the thermal management system to ensure that the method in which heat is added to or extracted from the cell does not induce an unfavourable temperature gradient that can compromise the benefit of absolute temperature control. Under electrical loads that involve relatively low current profiles, the cell volumetric heat generation rate is low and the formed temperature gradients between the external cell surface and core of the battery are of less concern [19]. However, the very poor effective cell thermal conductivity in the perpendicular direction (reported experimentally in the region of circa 0.25 W m À1 . K À1 [20e22]) can cause a large perpendicular temperature gradient to develop when the internal cell heat generation rate increases, for instance in hybrid electrical vehicle (HEV) applications [19]. Shah et al. reported that for a 26650 cylindrical cell, the gradient could exceed 20 C when subject to a steady state uniform heat generation of 6 W with aggressive active cooling on the outer surface [23]. The issue of internal temperature gradient is particularly problematic for cylindrical cells as the heat transfer pathway between the cell core and external surface is typically longer than for pouch cells, by circa: 44% for a common 18650 cell relative to a 10 mm thick pouch cell which can be cooled externally on both sides.
Due to the inherent lack of control over the increasing perpendicular temperature gradient associated with external cooling strategies under higher rates of volumetric heat generation, internal cooling strategies have been investigated to alleviate the thermal issues associated with the cell internals [7]. One such method, discussed by Sievers et al. [24], utilises a cooling channel to pass mineral oil through both the mandrel of a cylindrical cell and along its outer surface. Their simulations revealed that the internal temperature gradient could be reduced (for a set volumetric heat generation) upon increasing the size of the inner cooling channel at the cost of a decline in the cell volumetric energy density. More recently, Shah et al. [22] inserted a heat pipe into the mandrel of an experimental cylindrical cell where the condenser end protruded from the bottom tab. Under natural convection on the cell external radial surface with forced convection on the heat pipe condenser end, the temperature of the cell core was reduced by 20 C with the heat pipe when the cell heat generation rate was 1.62 W. Similar to the Sievers et al. solution, this method enabled the hot spot to shift from the core nearer to the outer surface of the cell, reducing both the magnitude of the temperature gradient and hot spot. However, the effectiveness of such a strategy relies on the presence of an external cooling mechanism to exist at the protruding heat pipe condenser section in addition to at the external radial surface of the cell, which would likely lead to additional complexities in the design of the thermal management system on the pack-level to accommodate cooling at multiple surfaces.
Contrary to surface cooling approaches, battery packs designed for tab cooling have the potential to be more compact given that cooling mechanisms in between cells are not required. Compared to surface cooled methods, tab cooling approaches for battery cells have the added benefit of providing an in-plane temperature gradient through the cells which can extend their service life [18]. However, one primary disadvantage of tab cooling is the presence of a longer heat transfer pathway given the length of conventional cells which can increase the magnitude of internal temperature gradients. Cooling both tabs as opposed to one can mitigate this issue, however, employing heat transfer mechanisms at both ends of the cell may potentially double the number of indirect cold plates and leakage sources.
In this paper, a new approach to internal thermal management of cylindrical battery cells is introduced to significantly enhance the thermal performance of a singular tab cooled approach, thus alleviating any pack-level complications associated with cooling at multiple surfaces of the cell. The thermal modelling development and description of the cooling strategy is discussed in Section 2, while validation of the baseline thermal model against experimental temperature measurements is discussed in Section 3. A case study analysis to rank the effectiveness of the internal cooling strategy against conventional tab cooling and outer surface cooling approaches, subject to the same aggressive electrical loading conditions, is conducted in Section 4. Further work and Conclusions are presented in Sections 5 and 6 respectively.

Model methodology
This section discusses the modelling methodology utilised for the thermal model, which includes the assumptions, battery heat generation characterisation and the numerical solution approach. Details of the analysed cooling strategies used in the case study analysis in Section 4 is also discussed.

Thermal model development
Homogenous bulk layer thermal models that use effective values for the cell heat capacity, density and anisotropic thermal conductivity have proven to be effective for thermal modelling of lithium-ion cells [25]. Examples of these approaches can be viewed in Refs. [21,23,26]. In these models, the complete cell (active material, electrolyte contact layer and metallic housing) is treated as one homogenous material displaying effective anisotropic thermal conductivity, heat capacity and density.
The newly proposed internal cooling method involves placing a heat pipe inside the mandrel of the battery cell, whereby both the condenser and evaporator ends of the heat pipe are connected to an aluminium heat spreader disc. This formed heat conduction network is termed the 'heat pipe system'. The top spreader disc in contact with the heat pipe evaporator acts to conduct heat through the top of the cell spiral roll (i.e. all the layers of the cell) into the heat pipe. The spreader disc in contact with the condenser end of the heat pipe acts to dissipate heat from the heat pipe across the base of the bottom tab which is actively cooled.
A schematic representation of the heat pipe cooled method is shown in Fig. 1(a). Conventional cooling approaches showing bottom tab cooling (b), radial/surface cooling (c) and double tab cooling (d) are also shown. Locations of the cell subject to a cooling source are represented by the presence of a heat transfer coefficient, these being h z0 at the bottom tab, h zL at the top tab and h r at the outer surface of the cell at radial location r ¼ R 0 .
The governing heat conduction equation for the cell composite, which includes the heat pipe and spreader discs, under the assumption of azimuthal symmetry is given by: Where r is the material density [kg m À3 ], T the local cell temperature [K], r the radial location [m], C p the material heat capacity [J kg À1 K À1 ], k r the perpendicular thermal conductivity in the 'r' direction and k z the axial thermal conductivity in the 'z' direction [W m À1 K À1 ]. This modelling approach assumes that the heat pipe behaves as a solid conductor which has an effective density, heat capacity and thermal conductivity [27,28]. The modelling implications employed for the heat pipe model have been used previously in related research [22,27]. The methodology behind the calculation of the effective heat pipe thermal conductivity is discussed in Section 2.2.
The solution of the governing heat conduction equation within the cell composite is achieved through the use of a finite difference approach using the widely employed, computationally efficient alternating direction implicit (ADI) method [29,30]. The general discretised form of the heat conduction equation in two dimensional (2-D) polar coordinates, where the thermal properties of the material may vary based on spatial location, is given below by the following ADI equations. These equations are valid for grid points between j ¼ 2:J and k ¼ 2:K as represented by the shaded area in the solution grid, see Fig. 1(e). The ADI equations in each direction with variable thermal conductivity properties are [29]: r direction: Where dz is the distance between nodes in the z direction [m] and dr the distance between nodes in the r direction [m]. The superscripts for the temperature and heat generation parameters refer to the time step n. For the thermal conductivity material properties the superscripts represent the direction in which the thermal conductivity is defined (i.e. z for axial and r for perpendicular). The subscripts for the physical properties refer to its location within the solution grid as defined by the spatial parameters j and k. Spatial variations of the thermal properties for the composite material is included through the use of the harmonic mean value for the thermal conductivity described by Patankar [31]. Use of the harmonic mean value as opposed to a linear interpolation between two grid points enables the interface thermal conductivity to physically capture steep changes in thermal conductivity gradients between two materials more accurately, and is recommended for use in modelling composite materials [31,32]. The harmonic mean values for the interface thermal conductivities are defined as: For polar coordinates the r terms in Eqs. (2) and (3) are discretised as [33]: The source term at the time interval nþ 1 2 is taken as a linear interpolation between the value for the source term at time interval nþ1 and time interval n: Newton's law of cooling is implemented at each external boundary condition: For the top of the cell: With the discretisation: For the bottom of the cell: With the discretisation: For the external radial surface: With the discretisation: For the centre of the cell at r ¼ 0 the insulation condition is applied due to symmetry: With the discretisation: The volumetric heat generation within the cell bulk material (not including heat pipe or spreader discs) is assumed to comprise of solely irreversible heat generation mechanisms, whereby entropic heating and other modes of heat generation [34] are ignored. This enables the heat generation term to be calculated from knowledge of the cell overpotential resistance as outlined in Ref. [35]. The irreversible heat generation term is given by: Where I is the cell current [A], R h the overpotential/internal resistance of the cell [U] and v c the volume of the bulk cell material [m 3 ]. The temperatures at each node and each time step is solved for using Matlab. For instances when there is no heat pipe or spreader discs, the spreader disc thickness is set to zero and the thermal conductivity of the heat pipe set to 2e-10 W m À1 K À1 to achieve an insulating condition at r ¼ R i .

Heat pipe physical properties
It is beyond the scope of the paper to discuss, in detail, the fundamental operation of a heat pipe. Heat pipes are an established technology that have been successfully deployed in a number of sectors [36], such as aerospace [37], automotive [16] and electronics [38]. A detailed description of their functionality is provided in Ref. [28] and a number of educational texts [39]. However, for completeness, a brief description is provided here. A heat pipe generally comprises 3 distinct layers, which include a wall shell layer, a wick structure and a vapour space. At the evaporator section, the vaporisation of the working fluid enables large quantities of heat to be absorbed isothermally. The formed vapour creates a pressure gradient within the heat pipe which channels the flow of vapour towards the condenser section where it is condensed. The condensed liquid is then pumped back towards the evaporator, via capillary forces created in the wick structure, to complete the cycle.

Estimation of physical properties
The thermal properties of the heat pipe used in the model are estimated through considering the design strategy of the heat pipe as outlined in Ref. [28]. The operating range of the heat pipe is chosen to identify a suitable working fluid. Given that lithium-ion cells prefer an operating range between 20 and 25 C, water is chosen as the working fluid.
For the wick structure a sintered copper wick is chosen with a pore size and permeability of 9Â 10 À6 m and 1.74Â 10 À12 m 2 respectively [28]. The thickness of the copper outer heat pipe wall is determined to withstand a 200 C bonding temperature to the spreader discs, whereby the thickness is calculated using the thin cylinder formula: Where t w is the thickness of the copper shell, [m], r hp the radius of the heat pipe [m], P the vapour pressure of water at 200 C [Pa] and U the permissible stress of annealed copper at 200 C which is taken as 18.5 MN/m 2 from Ref. [40]. For a 6 mm diameter heat pipe, t w is calculated as 0.3 mm. For simplicity, this thickness is taken as constant and used for smaller heat pipe diameters used in the simulation study. The size of the vapour space is taken as a design variable which sets the remaining wick layer thickness, thus dictating the power capability of the heat pipe. For correct operation of the heat pipe, the capillary pumping pressure must be equal to or greater than the sum of the pressure drops present in the heat pipe and is given by Ref. [28]: Where l e , l c and l a are the lengths of the evaporator, condenser and adiabatic section respectively. Assuming laminar flow of vapour and complete pressure recovery in the condenser then DP v is given by Ref. [41]: Where g is the gravitational acceleration [m s À2 ] and 4 the angle between the horizontal plane and the heat pipe. The maximum power of the heat pipe is determined through substituting the above expressions into Eq. (22) and solving for Q.
It is assumed that the thermal resistance for heat transfer through the vapour core is very small relative to other resistance terms, therefore the overall heat pipe thermal resistance network [41] is simplified giving an expression for the total heat pipe thermal resistance ðR hp ) as: Where R w is the thermal resistance through the copper wall and R wi the thermal resistance through the wick layer. The contact resistance between the heat pipe and cell material interfaces is neglected.
With reference to Fig. 1(a), the area of the condenser and evaporator sections are assumed to be those solely in contact with the spreader disc for determining l eff and R hp . The adiabatic section is taken as the remaining length of the heat pipe which passes through the cell mandrel. This relies on the assumption that the heat flux into the heat pipe from the cell material is very small relative to that from the heat flux at the aluminium spreader discs heat pipe interface. The validity of this assumption is tested further, through an extensive simulation study, discussed within Section 4.2.4.
The copper wall and wick thermal resistances are expressed as: Where r wi the radial location where the wick layer meets the copper wall [m], k c the thermal conductivity of the wall material (copper) [W m À1 K À1 ], k wi the effective thermal conductivity of the wick [W m À1 K À1 ] and l s the height of the aluminium spreader discs which is set as 2 mm. For heterogeneous copper sintered wicks, the Maxwell relation [42] is used to determine k wi which is given by: Where k l and k s is the thermal conductivity of the liquid water and solid particles (copper) respectively [W m À1 K À1 ] and ε the porosity of the wick, which is given as 0.52 based on the reported permeability and pore radius results of the sintered wick sample from Ref. [28]. The overall heat pipe thermal conductivity is calculated from the expression for R hp , assuming that it acts as a solid conductor given by: Where l hp is the length of the heat pipe [m] and k hp the overall effective heat pipe thermal conductivity [W m À1 K À1 ].

Experimental validation of the battery thermal model
The accuracy of the baseline thermal model (Section 2.1) without the heat pipe and spreader discs is compared to experimental surface data obtained from 4 new unaged 18650 type cells with a graphite anode (LiC 6 ) and Lithium Nickel Cobalt Aluminium Oxide (NCA) cathode. The nominal capacity of the test cells is 2.9 A h. The 4 cells are placed horizontally on the top tray within a climate chamber that circulates air to maintain its set point ambient temperature. A thermocouple is placed at the outer surface of each cell at mid height. Fig. 2 displays the test setup adopted prior to attachment of the thermocouples. Each of the 4 cells are subject to two separate electrical loading profiles, each performed at a different ambient temperature. The electrical loading profiles for both tests are shown in Fig. 3(b).
The physical properties for the 18650 cell used in the thermal model are displayed in Table 1. Given that a moderate degree of air circulation is present within the climate chamber, an h value of 25 W m À2 K À1 [43] is assumed at both the tabs and exterior surface of the cell in the thermal model. The properties of the mandrel in the thermal model are set with a thermal conductivity of 2e-10 W m À1 K À1 to provide the desired insulating condition. The spreader disc thickness is set as 0 mm.
The current pulse method [44] with a pulse length of 10s is applied at 4 SOC locations for the cell at both 25 C and 10 C to determine the overpotential resistance at these SOC levels. The average results over the 4 cells are shown in Fig. 3(a). Given that data points for <20% and >95% SOC were unavailable, the fitted polynomial is used to extrapolate the resistance values within these SOC regions for calculating the cell heat generation in the thermal model [19]. This is highlighted by the presence of vertical dashed lines in the simulation/experimental results comparison curves in Fig. 3(e) & (f), where between the lines the polynomial is interpolated between the data points.
The resistance polynomial is used to calculate the resistance value used in the thermal model at the given SOC level. Coulomb counting is used to estimate the cell SOC for a given nominal cell capacity and current profile via:  Fig. 3(c) and (d). The slight discrepancy in temperature between the cells are believed to be due to the increasing air convection occurring on the cells nearer the right of the chamber (in closer proximity to the air circulation inlet) leading to additional cooling. As a result of the variable heat transfer conditions experienced by the cells, validation of the thermal model is compared to the readings of Cell 2 which is more representative of the static h value of 25 W m À2 K À1 used in the thermal model. Fig. 3(e) and (f) highlight that there is good agreement between the baseline thermal model and T mid,surf data from Cell 2 for both test cases. The highest discrepancy occurs for the lower temperature test and corresponding drive cycle in the deep discharge region (20% < SOC), as seen to the right of the dashed vertical line. This is likely due to the heightened error associated with the polynomial extrapolation. The model gives a peak error of 11.3% which occurs for the 10 C test, with a mean absolute error of 1.94% across the whole 10 C test. Between the vertical lines bounded by resistance data points, the peak error and mean absolute error reduces to 3.15% and 1.11% respectively for the 10 C test.

Cooling case study analysis
In this section, the thermal model is utilised to investigate the thermal performance of two common cylindrical cell formats (18650 and 32113) as a function of cooling strategy when subject to two developed aggressive duty cycles. For the heat pipe cooled method, the calculated heat pipe physical properties and properties of the aluminium spreader disc material are input at the grid locations where these objects exist. In all instances, mesh resolution and time step parameters are chosen such that their values remain independent of the results. A steady state analysis is conducted to investigate the cell and overall thermal resistance for heat transfer between the cell hot spot and bulk heat transfer medium as a function of the convective heat transfer coefficient and location of cooling. An analysis into the effect of thermal conductivity of the heat pipe/solid conductor material on the cell thermal performance is also investigated. A full transient analysis using the full duty cycle profiles is undertaken to quantify the effect of cooling strategy choice on the magnitude of the temperature profiles. The thermal properties used for the cells considered in the simulation study are viewable in Table 1. In all instances, the initial cell temperature and temperature of the bulk heat transfer medium (T ∞ ) are set at 25 C.

Performance EV and HEV duty cycle
The current and calculated heat generation profiles for the performance EV cycle and HEV cycle are shown in Fig. 4. The methodology for construction of the EV duty cycle is beyond the scope of this paper and is discussed further in Ref. [49]. Both the EV and HEV cycle is obtained through a simulation with the IPG CarMaker software. For this case study a medium sized parallel hybrid vehicle is modelled for the HEV that represents a commercially available 5 door premium saloon. Table 2 presents the vehicle  parameters used. For the HEV, the battery control is specified such that the SOC window is constrained to between 70% and 20% SOC and the target SOC of the battery during operation is set to 50%, so as to avoid the limits of deep discharge and deep charge. The initial SOC of the battery at the beginning of the simulation is set to 70%. Electric driving is enabled above an SOC of 50% up to a speed limit of 100 km/h. Assistance through the electric motor is enabled at an SOC above 25% and becomes active if the vehicle torque demand exceeds either the maximum or optimum combustion engine torque. This vehicle model is used in conjunction with the speed-time trace of the US06 driving cycle to determine the HEV battery duty profile.
The full resistance curve at 25 C is used to calculate the heat generation from Eq. (20) given the SOC profile obtained from Eq. (33). For the HEV cell, as the SOC window is constrained to between 70 and 20% thereby avoiding large fluctuations in the resistance value, and given lack of a real 32113 cell to experimentally determine the resistance curve, a static value for the overpotential resistance of 4 mU for the cell has been assumed.

Steady state thermal analysis
The thermal resistance concept [50] is employed to rank the effectiveness of the cooling approaches as a function of cooling degree (i.e. 'h' value) applied at the associated cooling boundary surface(s). Thermal management strategies that are of an overall lower thermal resistance (R) can more effectively supply or remove heat into/out of the cell. For a transient analysis, a lower R value is indicative of a system which can more readily dampen out temperature oscillations [51] and would therefore be preferable for use within a system with large current fluctuations, e.g. in an HEV cell. The overall thermal resistance constant for heat transfer between the cell hot spot and bulk heat transfer medium is defined as: Where T max is the cell hot spot temperature [K], T ∞ the bulk heat

Calculated heat pipe parameters
The layer dimensions of the heat pipe which dictate the value for k hp are specified to meet the requirement for the heat pipe power Q. This requirement is set based on the time averaged heat generation values for both cell types. From Fig. 4, this is 3.11 W for the EV cycle and 5.80 W for the HEV cycle. Given that all the heat generated by the cell will not transmit into the heat pipe, the Q value for the 32113 cell is set as the maximum value obtainable with a 4 mm heat pipe which is 4.32 W. For the 18650 cell, the maximum Q value is set using a 3 mm heat pipe which is 3.60 W. Values are determined using water fluid properties from Ref. [28] at the lowest heat pipe operating temperature of 25 C to give the worst case. The 3 mm heat pipe is not deemed suitable for the 32113 cell given that the maximum obtainable Q is now only 2.03 W.
The determined values for R hp , k hp and the occupied volume fractions of each heat pipe layer relative to the overall heat pipe volume are shown in Table 3. Here, f w ,f wi and f v are the volume fractions of the copper wall, wick layer and vapour core layer respectively. The calculated values for k hp range from circa 8000e19,000 W m À1 K À1 and are comparable to similar analytical computed results reported in Ref. [41].

Thermal resistance as a function of cooling strategy
4.2.2.1. Overall thermal resistance. The results of the R value as a function of h and cooling strategy are displayed in Fig. 5(a) & (b). The number of h value data points used is shown in the curve for double tab cooling as an example. The effect of thermal conductivity of the solid conductor/heat pipe on the R value is shown in Fig. 5(c) & (d).
All cooling options involving a heat pipe or solid conductor include both the top and bottom 2 mm aluminium heat spreader disc.
For the 18650 cell, radial cooling enables the greatest heat removal for a given h value until h ¼ 520 W m À2 K À1 , where double tab cooling becomes more effective. This threshold increases to circa: h ¼ 1300 W m À2 K À1 for the case of the 3 mm internal heat pipe cell due to the much lower available heat transfer area and longer heat transfer pathway present associated with single tab cooling. Compared to tab cooling based methods, however, the R curve for radial cooling begins to level off much sooner and indicates that the thermal resistance for heat transfer through the cell material ðR cell Þ in the perpendicular direction begins to dominate the value of R at lower h values. The introduction of the heat pipe and spreader discs effectively lowers R cell , prolonging the region in which R is dominated by convective components of the resistance. Fully utilising the potential of the heat pipe system therefore requires a greater degree of heat removal at the base of the tab to address the higher heat fluxes created.
Upon transition towards the larger 32113 cell, the R value for double tab cooling equals that of radial cooling at a lower h value of 255 W m À2 K À1 . This is over half the value for the 18650 cell of 520 W m À2 K À1 and indicates that the effectiveness of tab cooling improves relative to radial cooling upon increasing the cell aspect ratio (i.e. cell diameter divided by height). Ideal cells for tab cooling would thus be short and stubby. Fig. 5(c) and (d) outline the impact of the effective thermal conductivity of the solid conductor on the R value. Increasing the thermal conductivity value prolongs the region in which R has a linear dependence with h on the log-log scale, signifying a lower R cell value. For a convective coefficient value of 3000 W m À2 K À1 (representing moderate liquid cooling with water [52,53]) for the 32113 cell, the overall thermal resistance decreases by 14.0% upon addition of a solid copper bar (k ¼ 398 W m À1 K À1 ) connected to the discs relative to single tab cooling without these additions. With a k hp value of 8225 W m À1 K À1 to simulate the heat pipe, R decreases by 49.1%. The much greater reduction in R clearly highlights the benefit of utilising a heat pipe over a simpler copper rod solution. For the 18650 cell at h ¼ 3000 W m À1 K À1 , the R value decreases by 20.3% with the copper rod and 51.5% with the 3 mm heat pipe with k hp of 9611 W m À1 K À1 .

4.2.2.2.
Cell material thermal resistance. For 1-D heat transfer cases (radial, singular tab & double tab cooling), the constant R cell can be expressed by: Where T s is the temperature of the bulk cell surface at the cooling source [K]. Knowledge of R cell is useful as T s is also the minimum cell temperature, therefore R cell can be used to calculate the magnitude of the cell temperature gradient ðDT max Þ for a given value of Q c . With the heat pipe system, the heat transfer is multidimensional and therefore T s (which is the temperature at the interface between the bottom spreader disc and bulk cell material) is not uniform. In this instance, the average temperature for T s across the interface is taken from the thermal model and Eq. (35) used to calculate R cell at logarithmic h value increments between 500 and 100,000 W m À2 K À1 . The average value of R cell over this range is taken and given to represent the constant value of R cell . Table 4 highlights the degree of reduction in R cell from the introduction of the heat pipe and spreader discs relative to single tab cooling without these additions. The error associated with R cell values used to estimate DT max for the heat pipe cases when compared to the direct DT max calculation from the thermal model is included. For the 18650 cell, a 3 mm heat pipe offers a 67.8 ± 1.4% reduction in R cell approaching the lower resistance offered by double tab cooling. Increasing the size of the heat pipe in the heat pipe system further past 3 mm offers diminishing reductions, indicating that larger heat pipes may not be worth the penalty of a reduced cell capacity as the size of the cell spiral roll would be reduced to accommodate the larger heat pipe. For the 32113 cell, the 6 mm heat pipe may be the most ideal given the diminishing reduction in R cell upon transition to the larger 8 mm size. Relative to a reference cell mandrel size of 3 mm, the 6 mm heat pipe results in a 2.67% reduction in cell capacity (which assumes that the loss in volume of cell bulk material directly relates to a loss in capacity).

Implications for heat transfer medium choice as a function of cooling strategy
Increasing the degree of convection is dependent on the fluid properties of the heat transfer medium, and may therefore dictate whether an air or liquid based battery thermal management system (BTMS) is required for a given cooling strategy and heat generation rate. An upper limit for the h value achieved with an air BTMS may be in the region of h ¼ 200 W m À2 K À1 [53], which would correspond to a situation where cross-flow of air exists over the cell to promote turbulence. An example of such a series type battery thermal management design is in Ref. [54]. However, as stated in Ref. [54], this system would consume a much greater amount of parasitic power to operate relative to a liquid cooled system due to the larger volumetric flowrates and hence higher pressure drop incurred by the air system [51]. Smaller mass flowrates to reduce parasitic power would be penalised through incurring a greater pack-level temperature gradient due to the inferior heat capacity of air. As such, such a high h value may not be practical to achieve with air cooling in a large sized EV battery with a long flow path.
Parallel type designs where air is blown axially through channels over the outer cell surface may benefit from improved pack level temperature uniformity [55] but suffer from a reduced h value relative to the cross-flow design. For a parallel type system, Kim and Pesaran [51] report an h value of 60 W m À2 K À1 , which from Fig. 5 would result in an R value of 9.01 K W À1 for the 18650 cell and 4.21 K W À1 for the 32113 cell. To achieve the same R value, Fig. 5(a) & (b) indicate that singular tab cooling would require an h value of 875 W m À2 K À1 and 680 W m À2 K À1 respectively for the 18650 and 32113 cell which would position it in the region of a water glycol system [51]. Tab cooling methods alone are thus ineffective with an air based BTMS and should only be considered when a liquid or evaporative [52] heat transfer pack-level strategy is required e.g. for packs that necessitate a large cooling power [56]. It is important to note that these heat transfer medium choices would require the use of an indirect thermal management approach, thus incurring additional thermal resistances between the base of the tab and the cooling/evaporator plate which is not modelled here. To ensure that the tab cooling method remains effective, care must be taken to minimise these resistances when implementing the external heat transfer mechanism.

Thermal performance for time averaged drive cycle heat generation
Steady state temperature contours throughout the 18650 cell subject to the time averaged heat generation of 3.11 W for the performance EV cycle are shown in Fig. 6. An h value of 60 W m À2 K À1 is chosen for the radial cooled case to reflect an air cooled parallel design and 875 W m À2 K À1 for the tab cooled cases (with and without heat pipe system) to reflect water-glycol cooling as discussed in Section 4.2.3. Fig. 6 visualises the effectiveness of the spreader discs and heat pipe in isothermalising the internal temperature gradients within the cell. Relative to radial cooling, the cell hot spot is shifted towards the outer surface slightly above the cell mid height. This enables the hot spot of the cell to be monitored at an accessible outer surface and is desirable from a control standpoint as estimations of the core temperature are not necessary. The total amount of heat entering the heat pipe, given the negative heat flux along its axial length at the heat pipe cell interface and heat pipe spreader disc interface is 1.52 W. Of this, the majority of the total heat (1.24 W) enters the heat pipe across the top spreader disc interface with the remaining 0.28 W entering across the heat pipecell material interface. This highlights the importance of the top spreader disc, as it acts as the dominant mechanism in transporting heat from the cell into the heat pipe. This may explain why the results in Ref. [22] did not show much benefit from transitioning from a solid copper rod to a heat pipe in cooling the core of the cell, as the poor perpendicular heat transfer pathway effectively hinders the supply of heat into the rod/heat pipe relative to its axial conduction. The much smaller heat of 0.28 W also supports the assumption used in the heat pipe effective length calculation (Eq. (25)) that the heat pipe evaporator section lies at the heat pipespreader disc interface, with the remaining contact area between the cell and heat pipe being relatively adiabatic. There still exists a more unfavourable perpendicular gradient between the layers of the cell with the heat pipe system. The effects of such a complex gradient on the overall cell degradation are not fully understood and would require further investigation (see Section 5). Nevertheless the magnitude of DT max with the heat pipe is reduced by a factor of 3 relative to both single tab and radial cooling. A summary of DT max , T max and volume averaged cell temperature (T vol ) values for each cooling strategy is shown in Table 5. The results for the 32113 cell are also presented when subject to the time averaged heat generation rate of 5.80 W. In this instance, the total heat input into the 6 mm heat pipe is 2.99 W. All heat pipe properties used are those displayed in Table 3.

Transient thermal analysis
A transient thermal analysis of the thermal performance of the cells for each corresponding aggressive duty cycle case is shown in Fig. 7 for the 32113 cell and Fig. 8 for the 18650 cell. The effective heat capacity and density of the heat pipe used in the transient analysis is calculated from the sum of the heat capacity/density properties of each material layer given its fraction of occupation within the heat pipe [27]: More aggressive values for h are chosen in this analysis, which for air cooling is 100 W m À2 K À1 e.g. to reflect a narrower channel hydraulic diameter for parallel cooling (at the cost of greater pressure drop and parasitic fan power requirement) and 3000 W m À2 K À1 for more aggressive liquid cooling. For double tab cooling, the h value is halved at the tabs. This is used to demonstrate that the heat pipe design has the potential to reach that of double tab cooling provided that the higher heat flux at the bottom tab can be addressed. From Fig. 7 for the radial cooled air case, T max reaches a peak of 50.8 C during the third loop of the US06 cycle with a corresponding perpendicular DT max of 17.7 C. For the single tab cooled case T max is slightly lower at 48.3 C but with a higher gradient of 18.9 C due to the increased degree of convection. It is expected that the tab cooled method will remain more thermally preferable over radial cooling given the reduced cell ageing rates associated with the direction of the in-plane gradient of similar magnitude [18].
Clearly, both methods result in high peaks for both T max and DT max and would be undesirable from an ageing standpoint.
For the double tab cooled case, T max peaks at 38.4 C with a DT max of 7.4 C. The corresponding values for T max and DT max achieved with the heat pipe system are 39.7 C and 9.7 C respectively, indicating similar thermal performance without the associated pack-level complications from multi surface cooling. The transient total power load on the heat pipe is also shown. As the peak power load reaches 6.8 W, the 6 mm heat pipe properties from Table 3 are not acceptable given that Q exceeds 4.3 W. Therefore, a thicker wick layer has been chosen to achieve a higher design Q value of 7.0 W at 25 C with a corresponding k hp value of 9861 W m À1 K À1 . For the EV performance cycle, similar observations are made. The lower thermal resistance values for double tab and heat pipe cooling dampen out the thermal oscillations relative to radial and single tab cooling, particularly during the deep discharge range where the heat generation rate spikes. Here, the transient heat pipe power load remains within the design limit of 3.6 W for the 3 mm heat pipe and therefore does not require modification of the wick structure.

Further work
Given the significant theoretical improvements in the thermal performance of the heat pipe system over singular tab cooling, further investigation should target the practical issues associated with incorporating the heat pipe and spreader discs into the cell internals. This should focus on identifying methods to avoid a potential short circuit given the electrical pathway created from the heat pipe and spreader discs, whilst avoiding additional resistances that can compromise the added benefit of utilising the axial conduction pathways within the cell. To further understand the effect of complex internal cell gradient on cell ageing, manufacture of a test cell specimen should be carried out together with an experimental analysis to quantify and compare the ageing rates of the heat pipe cooled method relative to that of conventional bottom tab cooling under identical cell electrical loading and external cooling conditions.

Conclusion
Relative to past applications of heat pipes for cooling cylindrical battery cells, the proposed heat pipe system discussed in this work utilises the more efficient axial heat conduction pathways present within the cell to further increase the rate of heat transfer from the cell into the heat pipe. This is achieved through connecting a 2 mm thick metallic (aluminium) heat spreader disc to both ends of the heat pipe (acting as the cell mandrel material) which directly contacts the top and bottom portions of the cell material. The thermal modelling highlights that the formed heat-conduction network enables a dramatic reduction in internal cell temperature gradient, owing to the reduced thermal resistance for heat transfer through the inside of the cell. This internal thermal management approach requires an external heat transfer mechanism to exist at only the base of the cell, which can simplify the pack-level thermal management design strategy. The heat pipe system also has the potential to be incorporated directly within the internals of individual battery cells for simpler integration within existing battery pack designs employing single tab cooling.
The immediate penalty of incorporating the heat pipe and spreader discs into the cell is a decrease in energy density and an increase in cell mass. Relative to the reference 18650 cell, the addition of the 3 mm heat pipe and 2 mm spreader discs reduces the cell energy density by 5.8% and increases the cell mass by 11.7%. For the 32113 cell, the cell energy density is reduced by 6.0% and the mass increased by 10.0% with the 6 mm heat pipe and 2 mm spreader discs.
As the heat pipe system increases the rate of heat transfer to the base of the cell, the degree of cooling applied at the base should be maximised to fully utilise the potential offered by the heat pipe system in minimising the cell temperature rise. This necessitates the use of either forced convection with a liquid heat transfer medium or evaporative heat transfer.