Probing Heterogeneity in Li-Ion Batteries with Coupled Multiscale Models of Electrochemistry and Thermal Transport using Tomographic Domains

This work presents a methodology for coupling two open-source modelling frameworks in a highly parallel fashion across multiple length scales to solve an electrical current and heat transport problem for commercial cylindrical lithium-ion batteries. The global current and heat transfer problems are formulated as resistor networks and solved using a ﬁ nite difference method on a network extracted from an X-ray tomogram of an MJ1 18650 battery. The electrochemistry governing the heat generation is solved at the local level using a physically parameterized model. Electrochemical models are solved for different regions of a spirally wound cylindrical cell in parallel, coupled via charge conservation at the current collectors in a “ battery of batteries ” fashion, similar to the concept of modelling a pack. Thermal connections between layers in the spiral winding are established and heat transport is solved globally in a two-dimensional fashion, allowing for the subsequent extension to three dimensions. Great heterogeneity in local current density is predicted by the model which is also found to have some temperature dependence with rami ﬁ cations for battery degradation.©2020The Author(s). Published on behalf of The Electrochemical Society by IOP Publishing Limited. This is an open access article distributed under the terms of the Creative Commons Attribution License

Lithium-ion batteries (Li-ion) are an integral component of a clean energy future. 1 However, with increasing demand for higher power and energy densities, it is becoming clear that traditional designs have their limitations from a thermal management perspective. [2][3][4][5] Internal temperature differences between the core and surface of a battery can lead to other non-uniformities in current flow, state-of-charge, particle stress and levels of degradation. [6][7][8][9][10][11][12] This is particularly the case for cylindrical cells with spirally wound electrode layers, termed the "jelly-roll," which are becoming more common for commercial applications with their adoption by Tesla. Additionally, under conditions of abuse, the Li-ion battery may experience temperature sufficient to trigger catastrophic uncontrolled reactions, termed thermal runaway. These events can trigger levels of degradation in the cell that are so severe that shortcircuiting and exothermic electrolyte decomposition occur, potentially leading to cell ruptures and explosions. Under these conditions, the thermally insulating nature of the materials used as electrodes and separators, and geometrical configurations that place large portions of the cell away from any heat sinks may exacerbate the thermal runaway and reduce the temperature range over which batteries may safely operate. [13][14][15][16][17][18][19][20] Many modeling approaches have aimed to capture and explain the key physical processes present within a Li-ion battery, having been inspired by the pioneering work of Newman and co-workers. [21][22][23] Atomistic models help understand the fundamental reaction pathways, crystallographic structures and intercalation processes at the nanoscale. Fully resolved direct numerical simulations, utilizing micron size resolution tomograms of detailed electrode structures, can be used to model sections of single electrodes and resolve the lithium concentrations and transport processes with high fidelity. 24 This micron scale is also well suited to investigate material properties and can inform constitutive relations at coarser length scales.
Overlapping with the electrode scale are the continuum models that treat the battery components as homogenized and volume-averaged domains, with transport governed by relations between ion diffusivity and bulk porosity, for example. 8,11,[25][26][27] Continuum models have been employed to investigate problems from multi-electrode cell scales up to pack-level problems. At the largest length scale, battery management system (BMS) models typically employ resistor-capacitor networks that can replicate the behavior of hundreds and even thousands of cells forming battery packs used in electric power trains and other systems. 28 Relatively few studies have attempted to combine models across different length scales in a direct way; instead, most rely on analytical relations to transfer information up the length scales. Often, this approach is suitable but assumes that information travels in one direction only. For transient problems, there is feedback between scales. For example, the generation of heat is best modelled at the particle level with careful consideration of the local state of lithiation in each electrode, as this determines the open-circuit voltage and also the changes in entropy. However, the global heat transfer problem is also determined by the rate at which the heat energy can be removed at the cell boundaries and transported between layers of the cell, so must be resolved at larger "entirebattery" level scales. The cells' outer boundary conditions are also dependent on pack designs and cell positioning within packs, adding an extra scale to consider. With definition of the appropriate boundary conditions at the particle level, it is possible to approximately model cell thermal transport over time with a lumped approach, but information is averaged and cell-level design features such as the placement of current collector tabs, that can transport both current and heat, are not resolved.
Harb and LaFollette combined a 1D electrochemical model employing porous electrode theory with a 2D conservation of charge model to determine the behavior of a spirally wound lead acid cell that included temperature dependence. 29 Their algorithm was an iterative two-step solution procedure, first calculating the heat sources from the electrochemical reactions under constant current and then the temperature distribution by integrating the energy equation over the global domain. Once temperature profiles were obtained, the cycle repeats and the concentrations and electrode parameters such as porosity were updated in a step-wise fashion. This iterative process was explicit and required a relatively small time-step to avoid oscillations in the solution. Song and Evans applied a similar methodology to study the temperature-dependent properties and profiles in a prismatic Li-ion cell stack. 30 They computed a 2D temperature profile and modelled each cell in the stack with a 1D electrochemical model assuming uniform and constant temperature between iterations of the heat problem with results comparing well with experiment.
Gerver and Meyers 31 developed a similar method to investigate the behavior of a Li-ion pouch cell in 3D. They modelled the current collector foils as a 2D network of resistors and the electrode and separator layers in-between resistor nodes as a system of individual 1D batteries connected in parallel: an implementation of a "battery of batteries" approach that is often used to model packs. The overall approach is similar to resistor-capacitor networks with tunable nonlinear resistance but instead of using fitting parameters, real physical parameters are used to determine the local electrochemical behavior. Gerver and Meyers's implementation method was to run the 1D model for a range of current densities as a pre-processing step to determine the I-V characteristic relations which were then used in the resistor networks. These characteristic equations were periodically recalibrated for changing states of charge and temperature. In later work, McCleary et al. used the same modelling framework to investigate current density and temperature distributions in spirally wound and prismatic cell configurations and investigated the influence of increasing the number of tabs and current collector thickness and also different cooling conditions. 32 Their simulations show that increasing the number of tabs results in more uniform current density and temperature distributions throughout the cells.
With advances in computing power and newly developed powerful and efficient asymptotically reduced methods for solving the battery problem, it is no longer necessary to rely on the calibrated curve method. Here, a methodology is presented that combines the precision and realism of the particle-level electrochemical modelling approach, with the versatility and efficiency of the continuum modelling approach. Open source software is used in both cases: PyBaMM for the local electrochemistry 33 and OpenPNM for the global current and thermal transport. 34 OpenPNM stands for open source pore network modelling software and has been used to model transport in various different porous materials, [35][36][37][38][39][40] in this case we utilize it's ability to model resistor networks for arbitrarily defined graphs. The method is demonstrated for a two-dimensional heattransfer problem resolving the "jelly-roll" electrode structures of spirally rolled cylindrical batteries using real structures obtained from tomogram data of a commercially available LG MJ1 18650 cell. A single discharge of the cell is performed and the effects of Crate, cooling conditions and the number of current collector tabs are explored.

Method
Electrochemistry-Model equations.-The electrochemistry is described by considering the limit in which the ratio of the distance between the current collectors (i.e. the local thickness) compared with the total spiral length is small, and the electrical conductivity of the current collectors is large compared to that of the electrode material. In this limit, charge transport in the current collectors is one-dimensional along the spiral length, whereas charge transport in the electrodes and separator is one-dimensional in the direction normal to the current collectors. This results in a so-called "N+1D" model which comprises a collection of 1D electrochemical models describing the local cell behavior coupled via an ND electrical model in the plane of the current collector. Calculations made in the current collector domain will be referred to as "global" and calculations made for the electrode and separator sections between current collectors will be referred to as "local." As the dimensions of the current collector are much longer along the length of the winding compared with the height of the battery a 1D approximation is valid, ignoring the fact that tabs do not span the entire length of the cylinder.
The behavior in the local electrochemical models is described by porous electrode theory, using the model of Doyle, Fuller and Newman (DFN). 22,41 The DFN model comprises equations for charge and mass transport in the solid and electrolyte, and also describes the electrochemical reactions that transfer charge between the active solid material and electrolyte. Heat transport occurs on the global (jelly roll) scale, and each of the 1D electrochemical models is held at a constant temperature, fixed by its location around the spiral, during each iteration. In certain operating regimes, for instance, low C-rate discharge, other simpler models, such as the Single Particle Model, may be sufficient to describe the local electrochemistry and provide improved computational performance. Such models can be derived from the DFN model via formal asymptotic analysis (e.g. 33 ).
Under the 1+1D assumption, the flow of current in the current collectors is primarily along the winding direction, z. The potential in the current collectors F  cc is assumed to depend only on the position around the spiral so that F = F   z t , .
where  i z t , cc ( ) is the current density in the current collectors,  L cc is the thickness of the current collector, s  cc is the electrical conductivity, and I z t , ( ) is the local through-cell current density. The potential is set to zero on the negative tab(s), and a total current I app is drawn from the positive tab(s). Elsewhere, the boundaries are assumed to be insulated as shown in Fig. 2b. This gives the boundary conditions: where + A tab is the area of the positive tabs. Current flow in the electrodes is assumed to be one-dimensional in the direction normal to the current collectors, x, and is also described by Ohm's law. Charge conservation then leads to ) is the potential in the solid electrode material, ) is the current in the solid electrode material, s  is the electrical conductivity,  a is the surface area per unit volume, and  j x t z , ; ( ) is the electrochemical current density that describes charge transfer between the solid and electrolyte. Boundary conditions for the electrode/current collector interface come from satisfying continuity of potential and current. At the electrode/separator interface there is no flow of electronic current: the charge must be transported through ionic current through the separator.
In the electrolyte, the ionic current is driven by potential and concentration gradients, and the effects of interacting species must be accounted for, resulting in a modified Ohm's law 33 ) is the potential in the electrolyte, i x t z , ; e ( ) is the current in the electrolyte, c e is the concentration of lithium-ions in the electrolyte, T z t , ( ) is the (local) temperature, s c T , e e ( ) is the electrolyte conductivity, e is the volume fraction, b is the Brüggeman coefficient, + t is the transference number, R is the universal gas constant, and F is Faraday's constant. At the electrode/separator interface the potential and current must be continuous. At the electrode/current collector interface, no charge is transferred via ionic current.
Mass transport in the electrolyte is governed by a reactiondiffusion equation 33 ) is the electrolyte diffusivity. At the electrode/ current collector interface the flux of ions is zero.
In the DFN, the active material particles are assumed to be spherical, with the mass transport described by Fick's law, and uniform lithium transfer across the outer surface. As a result, the concentration in the solid  c s is governed by ) is the solid diffusivity, r is the radial coordinate within the particle, and  R is the particle radius. The concentration within the particles is assumed to be initially uniform in space. In the DFN framework, there is a particle as each location, x, across the electrode and the solid concentration depends both on radial position and the position within the cell, i.e. =   c c r t x z , ; , . s s ( ) The electrochemical reactions are modelled using symmetric Butler-Volmer kinetics. The electrochemical current density is then given by 33 where the quantity  j 0 is the exchange current density, ) is the open-circuit potential. Note that the expressions for the electrochemical reactions are evaluated on the surface of the particle =  r R .
Global heat transfer-model equations.-Heat transport is calculated based on the assumption that the battery is a homogenous medium with volume averaged thermal properties. The thermophysical properties are taken from Chen et al. 42 and presented in Table I. The energy balance equation can be written as where r is the density, C p is the specific heat capacity at constant pressure,, k is the thermal conductivity which is considered to be anisotropic, as described later, and ¢ q are the heat sources: where DS is obtained from measurements of the change in open circuit potential of the cell w.r.t. temperature (dU/dT), including contributions from both anode and cathode. 43,44 A finite difference approach is used to solve the heat transport problem utilizing a heat resistor network employing Fourier's Law of conduction. 34 Key analytical relations and their temperature dependence.-The electrolyte conductivity is that of LiPF 6 in EC:EMC, which is a function of ion concentration with polynomial fit from Nyman et al. 45 : where c e is the concentration of Li ions per liter. The electrolyte diffusivity is taken from the same source and is also expressed as a function of ionic concentration: Both the electrolyte conductivity and diffusivity, as well as the electronic diffusivity and reaction rates, are assumed to follow Arrhenius temperature dependence. For example, electrolyte diffusivity is scaled as follows: where E a is the activation energy with assumed values presented in Table II, T ref is the reference temperature, 298.15 K, T is the local cell temperature and R is the universal gas constant. The equilibrium potentials and entropic changes for each electrode (graphite with silicon anode and NMC-811 cathode) are interpolated from the experimental data gathered by Sturm et al. for an MJ1 battery 44 and shown in Fig. 1.
X-ray CT.-Lab-based micro-CT was achieved through the use of an Xradia 520 VERSA X-ray instrument (Zeiss 520 VERSA, Carl Zeiss Inc., CA, USA). All scans were conducted with a stationary tungsten anode target on a copper substrate that produces a polychromatic cone-beam geometry with a characteristic intensity peak (W-Kα). around 58 keV, as a result of a tube voltage and power of 120 kV p and 10 W, respectively. The cell was immobilised using double-sided tape onto a flat-top 360°rotational stage. Radiographs were collected at a 0.4× magnification through 360°of rotation, 4500 in total, each with an exposure time of 10 s. This produced a total scan time of approximately 12.5 h per tomogram. The magnification produced an isotropic voxel size of 10.4 μm (after reconstruction). Sections of the battery were imaged sequentially by moving stage vertically in z between acquisitions. Representative slices from the upper, middle and lower part of the battery were selected from the central part of each scan to minimize beamhardening effects. After acquisition, the 2D radiographs were reconstructed into 3D tomograms using commercial software employing filtered back projection algorithms ("Reconstructor Scout-and-Scan," Carl Zeiss., CA, U.S.A.).
Image processing.-A custom image processing routine was devised using several functions from the Scikit-image 47 Python package to segment the image into the separate electrode layers of the jelly-roll and identify current collectors. This consists of the following steps: • Adjust and level the image intensity to account for radial beamhardening effects.
• Remove the outer features, including the can, by centering the image and applying a mask of zero outside a fixed radius.
• Apply a binary threshold to isolate features with high intensity.
As the cathode electrode material and anode current collector appear brightest, this clearly separates the image into two distinct brighter regions.
• Apply a binary dilation to close any erroneous small features inside the high-intensity regions.
• Skeletonize each region with a medial axis transform to find the central lines in regions corresponding to the current collectors.
• Remove the side branches of the skeleton which are image artefacts by applying a rank-sum to the skeleton to identify junctions and end-points and removing sections between them that are not part of the main branch.
• Apply a "dart-board" image mask with lines at fixed arc angles separated every 10°. Where these lines intersect the current collector lines a node or resistor junction is created for simulation.
• Create connections between nodes for electrical and thermal transport. Nodes connected in series define the current collector domain with equivalent resistors connecting the current collector radially in place of the batteries for the global current calculation.   Solution procedure.- Figure 3 outlines the algorithmic steps to solve the coupled electrochemical thermal problem. There are essentially three important submodels: the DFN models which are local and the Equivalent Circuit Model (ECM) and Heat Resistor Network (HRN) which are global. The DFN models the electrochemistry of the battery segments, the ECM is used for matching the applied global current with the combined local current densities of each segment and for specifying the boundary conditions of the DFN models which change for each successive step. The HRN is used for the temperature update using the heat sources from the DFNs. Both the ECM and HRN use the same network topology displayed in Fig. 2c and OpenPNM 34 is used to solve Ohm's Law and Fourier's Law of electrical and thermal conduction, respectively. Details of how this is done can be found in previous work. 37,48 As the boundary conditions and local temperature for the DFN are held fixed for each time step, the time step is kept small to ensure reasonable continuity. A single DFN model object is set up with local current, local temperature and current collector length as input parameters. The model is run many times in parallel with varying input parameters and boundary conditions representing the local behavior in the electrode sections that are segmented by the image analysis. Apart from the length of the current collector in each battery segment, the other physical dimensions for the model are the thickness of each component in the electrode layers and the cylinder height. The thicknesses are determined from image analysis and detailed in Table I and the cylinder length is set to be 65 mm as the battery is an MJ1 18650.
Significant potential or thermal gradients are not expected to occur in the 3rd (cylinder length) dimension because low currents are being simulated and it is assumed that the surface of the battery is evenly cooled, although uneven cooling is investigated. Even if tab cooling strategies were to be employed, then significant thermal gradients are not likely to occur because the cell casing is connected to the tabs for cylindrical cells. For an 18650, if the jellyroll electrodes are un-rolled, they can stretch to over 1 m. In comparison, a cylinder height of 65 mm is a much smaller distance over which to establish potential or thermal gradients. Even though the tabs are not typically extended for the full length of the cylinder, the model essentially assumes that they are and that current non-uniformity in this direction is assumed to be negligible when compared to the length around the jellyroll.
The DFN models are all initialized at the same state-of-charge and with a guess for the initial current flow through each model being equal to the terminal current scaled by the ratio of the current collector segment length to the total current collector length. This is only semi-realistic because nodes farther away from the tabs will experience a greater circuit resistance compared with those closer to the tabs and will, therefore, contribute less current initially. However, the initial step is only used to determine a guess for the equivalent resistance of the DFN models used in the ECM. The models are run for one very small time-step of 1 μs under these conditions and the local overpotentials are divided by the local currents to get the local equivalent resistance according to Ohm's law. These resistances are used to update the resistor network and this is solved by iteratively adjusting the global terminal voltage until the total current matches the desired applied terminal current. Terminals are specified with great flexibility by simply picking nodes in the resistor network to apply boundary conditions. Three different scenarios for terminal placement are investigated and detailed later in Table III.
With a better initial guess for the local currents, including the circuit resistance effects of the losses in the current collectors, the boundary conditions for each DFN are updated and they are run for one normal time step with fixed local current. This time-step is chosen to be relatively small (30 s for a current of 1 A) and is inversely scaled by the applied current so that the total charge transferred always remains small between steps. The total capacity of the battery is determined to be around 3.5 Ah and so for 1C discharge, approximately 420 time steps are required to simulate a full discharge from a fully charged state.
After stepping the electrochemical models with fixed local current the total heat generated is calculated and the global temperature is updated using an effective medium approximation. The same topological network is used with thermal conductivity calculated based on resistors-in-series for the electrode layers in the radial direction and resistors-in-parallel for the spiral length direction. Fick's 2nd law is integrated over the same time period with a single boundary condition applied to all the nodes on the outer casing layer of the model. This generates a temperature map from which to update the local temperatures for the next iteration of the electrochemical models.
The process repeats until one of the termination criteria is reached: either the total time requested has elapsed or the concentrations in the particles or terminal voltage has reached the limit of their allowed range. No limit on temperature is  specified for the present framework but would be a consideration for control engineers. The steps are summarized and pictured in Fig. 3. The DFN models are completely agnostic to their position in the resistor network as ion transport is assumed to be 1D. Therefore, they can be solved in parallel, which results in a highly efficient solution procedure. A full discharge using the DFN model for ∼1300 nodes completes 300 steps in approximately 30 min using a standard desktop computer.
Thermal material properties and boundary conditions.-For the heat transport problem, a homogenous medium with average density and specific heat capacity is assumed. Anisotropic thermal conductivity is assumed in the directions parallel and perpendicular to the current collector windings, termed the spiral direction and the radial directions. The current collectors have intrinsic conductivity several orders of magnitude greater than the surrounding components. However, this is counteracted somewhat by their small thickness and studies have found that the parallel effective conductivity is approximately one order of magnitude greater than the perpendicular direction. 49 The thermal diffusivities for the lumped model in the spiral and radial directions are calculated as follows: where the summation is over the battery components comprising the unit cell from the center of the anode current collector to the center of the cathode current collector. When evaluated with the parameters in Table I  Other battery parameters.-Whilst utilizing a specific jelly-roll geometry for the modelling domain, the focus of this paper is not to analyze a specific battery and evaluate its performance. Rather, the aim is to demonstrate the effectiveness of the coupling methodology in spatially resolving important parameters with high fidelity in an efficient manner to provide the foundation for future parametric analysis. Experimental work is required to characterize a specific battery as there are many important parameters to determine and for experimental validation, care must be taken to evaluate as many of these as possible. Typical values are found from the literature for batteries and chemistries as close to the MJ1 example cell as possible and so the results should be considered qualitative and not quantitative. Future studies combined with thorough experimental characterization of specific battery designs and chemistries are on-going.

Results and Discussion
Cases under investigation.-Cases are referred to alphabetically according to Table III: For all cases, the C-rate is varied between 0.5, 1.0 and 1.5 C, which equates to 1.75, 3.  50 but the details of which are out of the present scope of the study. It is our intention to couple this modelling framework with future CFD analysis providing more in-depth information on different cooling mechanisms.
The fraction of the outer surface area being cooled is also changed to 1/3 for the 1 tab cases (P-T) to simulate the effect of a cooling channel only running past one side of the battery, as often occurs with liquid cooling.
The effect of C-rate.-The MJ1 battery has a nominal capacity of 3.5 Ah, although this is dependent slightly on the C-rate, with higher rates resulting in lower battery voltages and so reaching the lower voltage cut-off at lesser discharge capacity. Figure 4 shows that for the three C-rates investigated in the present study, the capacity is relatively unchanged. The cooling conditions are the same for each C-rate and so with the same heat energy produced over shorter times at higher C-rates, but dissipated at the same rate, the battery temperature increases in proportion with the C-rate. There is also a minimum in the temperature profile at the lowest C-rate caused by the negative entropic coefficient in the anode around 0.5 SOC. For the higher C-rates the increased Joule heating compensates for this and the heating load is greater than the cooling load.
The particular battery chosen for the study is not a high power cell and so investigating higher C-rates was not advisable by the manufacturer, but will be performed in a future study using more suitable cells. The pattern for the variation in local current density is similar for each current but with higher rates producing a greater relative variation. The spatial position in the jellyroll for the minimal and maximal local current density changes over time. To demonstrate this, Fig. 5 provides a snapshot of the local current density at various states-ofcharge. In the base cases (A-E), the positive tab is set to the innermost node on the positive current collector and the negative tab is set to the outermost node. As the negative current collector is copper and the positive is aluminium, and both have the same thickness, the electrical resistance in the positive current collector is about twice that of the negative current collector. Therefore, the path that minimizes electron transport through the positive foil is the least resistive and so initially, as all parts of the battery are at the same temperature and SOC, the regions closer to the positive tab are favorably discharged and the local current density hotspots can be seen nearest the center of the jellyroll. In the equivalent circuit, the equivalent resistance of the battery is changed over time as the overpotential is dependent on both temperature and SOC, which both change locally during discharge. The innermost sections heat up more quickly than the outermost sections, both because there is more current produced here and therefore more heat, but also due to the innermost sections being furthest away from the cooled outer surface. The increase in temperature reduces the local resistance slightly due to increased ion diffusion and conductivity, and more favourable reaction kinetics based on the assumptions of Arrhenius-like dependence of these properties. However, the reduction in these overpotentials is counterbalanced by the faster reduction in local SOC and the associated reduction in OCV, which leads to an overall increase in local equivalent resistance and a shifting of current hot spots to the outer sections. The process happens gradually, at times in a wave-like ripple fashion, but sudden changes can also occur, shifting the maximum back to the inner sections when the equivalent battery resistance drops in outer sections and the resistance in the current collectors again dominates the profile of the current density (animations of this process are provided as Supplementary Information is available online at stacks.iop.org/JES/167/110538/mmedia). The difference between the minimum and maximum current density can be as great as 15% of the mean current density and this relative difference increases with C-rate. During a constant current discharge, the total amount of current being produced in each section of the battery should sum to the same amount, within certain tolerances. However, the local rates ebb and flow more severely for some sections compared with others and this could impact local degradation.
The effect of the number of tabs.-In the present study, increasing the number of tabs simply means that on average the current has to travel a shorter distance to enter or leave the cell. Numerically, this is achieved by setting the applied voltage and drawing current at intermediate nodes along the length of the jellyroll. The network is not physically changed to accommodate the additional material that tabs would require; this will be considered in future reports as early studies indicate that local geometric discontinuities can provide starting points for stressrelated degradation, 51,52 and indeed, a method based on real geometry will be the only way to effectively determine the electrochemical forces influencing such degradation. The capacity of the cell will also slightly reduce as there will generally be less active material per unit volume. These limitations aside, the next set of results provide interesting and useful information regarding the likely temperatures and current inhomogeneity that would occur in multi-tabbed cells.
Comparing cases A, F and K, the least cooled cases with 1, 2 and 5 tabs respectively, it can be seen in Fig. 6 that for the highest C-rate, a temperature difference of ∼8 K is found between the simple 1-tab case and the multi-tab cases, which behave similarly. The magnitude of current inhomogeneity is also clearly reduced for the multi-tab cases.
Comparing the snapshots of current density in Fig. 5, the base case, and Fig. 7, the case with an additional central tab on each current collector, it is clear that the position of the initial hot-spot for current density has shifted to a more central position. In both figures, the wave-like switching between relatively high and low current generation is clear but now the difference between the minimum and maximum local current density is never more than about 3%.
Although informative, it is somewhat difficult to appreciate the whole pattern of the variation in local current density over time from the snapshots. Imagining that the jellyroll has been unrolled and laid flat, with the innermost section on the left-hand side, the percentage deviation from the mean current density can be presented over the entire discharge period as a "current density map" (Fig. 8). It is useful to refer to greater than average current density as "overcurrent" and lesser as "under-current." Comparing cases A, F and K for all C-rates, it can be seen that bands appear along both axes. For case A, the only bands are along the time/capacity axis. They can be categorized approximately as one of three spatial regimes: "innerfavoured," where the over-current is nearer the centre of the jellyroll, approximately "evenly distributed" and "outer-favoured." Approximately half of the charge is transferred under the evenly distributed regime, which is a neutral result for potential battery degradation. However, if the battery were to be cycled with higher frequency around a fixed SOC with zero-mean current, one would expect the sections with the greatest over-current to experience the highest charge transfer and therefore degradation. This corresponds to the sections nearest the positive tab for the single tab case.
For cases F and K, with 2 and 5 tabs respectively, the spatial influence on the current density is clear, with bands appearing prominently along the x-axis (normalised position) as well as the y-axis (state-of-charge/time), as with case A. The increase in the number of tabs has the effect of smoothing out the current density variations in general, but there are still points in space and state-ofcharge where local variation between minimum and maximum could be as much as 20% for the 2-tab case and 6% for the 5-tab case. The base case can experience a relative difference in local current density of 40% near the end of discharge. This percentage difference appears to increase with total applied current as the scale bars in Fig. 8 clearly increase from left to right.
The effect of cooling conditions.-To assess the appropriateness of the thermal boundary conditions imposed in the model, an MJ1 battery was cycled at the investigated C-rates inside an environmental chamber at 25˚C. Three cycles were performed with an hour left between charging and discharging to allow the battery temperature to equilibrate with the environmental temperature. A thermocouple was placed on the outer can at half-height and the results are compared with the h = 28 [W m −2 K −1 ] simulation cases in Fig. 9. It can be seen that approximately the initial 1/3 of the temperature profile matches the simulations quite well. Afterwards, the experimental temperatures continue rising, whereas the simulation temperatures begin to fall slightly. Then towards the end of discharge where overpotentials are highest, the temperatures begin to converge again with simulation rising faster than the experiment. As mentioned previously, the OCV and entropic data was obtained for an MJ1 cell, but gathered from a previously published study by Sturm et al. 44 In fact, on analysis of the temperature profiles gathered by Sturm et al. there appears to be a similar discrepancy between their Figure 9. Comparison of experimental and simulated surface temperatures for MJ1 18650 battery using the base case number of tabs. Figure 10. Comparison of different heat transfer coefficients with increasing heat loss (5, 10, 28, 50 100 [W m −2 K −1 ] through Case A-E) at an applied current of 5.25 A with even cooling around all outer surfaces. As cooling increases, cell temperature decreases and terminal voltage decreases. The total amount of heating actually increases with increased overpotential but the cooling more than compensates for this. simulated results and the experimental values. An additional discrepancy could also be that there is some temperature dependence on the heat transfer coefficient which could be adjusted in the model but requires further analysis with a more sophisticated model for the environmental surroundings. The values for thermal conductivity and specific heat capacity in the model are also the literature values and will be experimentally determined in future studies. The maximum difference between the numerical and experimental values appears at roughly 50% SOC and can be attributed to the reversible entropic heating term associated with the anode. This behavior could also be dependent on the amount of silicon used in the anode and warrants further experimental evaluation-preferably with varying amounts of silicon. However, this work is out of the present scope.
The effect of changing the heat transfer coefficient to different constant values in order to simulate natural and forced convection cooling regimes can be seen in Fig. 10 for cases A-E. Increasing the heat transfer coefficient lowers the temperature of the battery but increases the thermal gradient inside the battery, which may be important to consider at higher C-rates. 53 There is little effect on the current density distributions with even cooling around the whole external surface. However, many battery cooling systems implemented packs do not feature even cooling around the whole circumference and this will be investigated in the following section.
The effect of uneven cooling.-All previous results have assumed uniform cooling around all outer surfaces of the battery. A more realistic scenario for thermal management in larger packs where a cooling channel runs in-between rows of batteries, heat rejection will be heterogeneous over the outer surface of each battery and likely affect individual cells in different ways depending on location in the pack.
The temperature profile for case T is shown in Fig. 11 and was produced with nearest neighbor interpolation without smoothing, so sharper interfaces are produced where the lateral thermal gradients nearest the cooled surface are higher. These are a feature of the coarseness of the thermal resistor network on the outer rolls and so for higher C-rates, the number of nodes may need to increase. However, as discussed previously, the thermal conductivity in the spiral direction is at least an order of magnitude higher than the radial direction and so gradients in this direction will always be smaller.
A comparison of the current density and temperature over the whole discharge period is shown in Fig. 12 for the different cooling cases. The battery temperature of cases C and T is roughly similar because the total amount of heat being taken from the battery is  approximately the same on a volumetric basis. However, the uneven cooling creates a larger difference between minimum and maximum temperature, as one would expect. The current density minima and maxima are approximately the same in all compared cases but the spatial position of these extrema is temperature dependent.
The effect of uneven cooling on the distribution of the current density is shown in Fig. 13 with the minimum and maximum depending most on the local SOC, as found previously with even cooling. A decrease in local temperature near the cooling surface has the effect of increasing the current inhomogeneity slightly, most clearly shown at 0.5 SOC, but at the investigated C-rates the greatest internal temperature difference is still only a few degrees and so the Arrhenius temperature dependence of the key transport variables such as the electrolyte diffusivity and also reaction rates will not vary significantly either.

Conclusions
For the first time, a coupled methodology for solving electrochemistry and thermal transport using tomographic domains at the cell level is presented. The software utilized for the procedure is all open source and written in Python which is accessible to the scientific community. The solution procedure is very efficient but further optimization is still possible and ongoing; however, it compares favorably to previous reports. 31,32 Image processing techniques are outlined to identify the current collectors in tomographic images of a single jellyroll with one inner and one outer tab. Nodes along the spiral length of the current collectors are defined for an equivalent circuit type approach to solving the global heat and current transport.
Interesting wave-like variations in local current density are observed and are attributed to the local equivalent resistance changes in the battery electrode sections. Whilst the model is parameterized for a low power, high energy cell, reasonably high temperatures are predicted to occur along with internal temperature gradients and the cooling strategy is demonstrated to have a pronounced effect on the average, extrema and location of temperature and current density hotspots. Further work to include degradation mechanisms and parameterize the model for higher power cells is intended for future work. Due to the relatively small internal temperature gradients limited by C-rate and cooling rate, the spatial-thermal dependence of battery states is not very prominent but is expected to increase for high C-rates and the case of more extreme cooling.
Because the model treats a battery as a collection of smaller batteries connected electrically in parallel (a battery of batteries), the same solution procedure can also be applied to packs of batteries with different considerations given to the thermal transfer between battery nodes and cooling systems. The boundary conditions for current collector tabs are varied and demonstrated to have a profound effect on the current density homogeneity with more tabs decreasing the internal resistance to current flow by reducing the path length. This, in turn, leads to lower cell temperatures with lesser Joule heating and it is hypothesized lower rates of degradation. The mechanical influence of the tabs is not considered in the present study and is left for future work but could be incorporated into the current image-based modelling framework. Beyond multiple tabs this framework creates a method for probing real behaviour based on real geometry which will be essential for understanding the role of geometric imperfections or stress related degradation.